[7] | 1 | /*! |
---|
| 2 | * \file |
---|
| 3 | * \brief Matrices in decomposed forms (LDL', LU, UDU', etc). |
---|
| 4 | * \author Vaclav Smidl. |
---|
| 5 | * |
---|
| 6 | * ----------------------------------- |
---|
| 7 | * BDM++ - C++ library for Bayesian Decision Making under Uncertainty |
---|
| 8 | * |
---|
| 9 | * Using IT++ for numerical operations |
---|
| 10 | * ----------------------------------- |
---|
| 11 | */ |
---|
| 12 | |
---|
| 13 | #ifndef DC_H |
---|
| 14 | #define DC_H |
---|
| 15 | |
---|
[262] | 16 | |
---|
[180] | 17 | #include "../itpp_ext.h" |
---|
[7] | 18 | |
---|
[219] | 19 | /*! |
---|
| 20 | \defgroup math Auxiliary math functions |
---|
| 21 | @{ |
---|
| 22 | */ |
---|
| 23 | |
---|
[7] | 24 | using namespace itpp; |
---|
| 25 | |
---|
[37] | 26 | //! Auxiliary function dydr; dyadic reduction |
---|
| 27 | void dydr( double * r, double *f, double *Dr, double *Df, double *R, int jl, int jh, double *kr, int m, int mx ); |
---|
| 28 | |
---|
| 29 | //! Auxiliary function ltuinv; inversion of a triangular matrix; |
---|
| 30 | //TODO can be done via: dtrtri.f from lapack |
---|
| 31 | mat ltuinv( const mat &L ); |
---|
| 32 | |
---|
[427] | 33 | /*! \brief Abstract class for representation of double symmetric matrices in square-root form. |
---|
[7] | 34 | |
---|
[37] | 35 | All operations defined on this class should be optimized for the chosen decomposition. |
---|
[7] | 36 | */ |
---|
[22] | 37 | class sqmat |
---|
| 38 | { |
---|
| 39 | public: |
---|
| 40 | /*! |
---|
[75] | 41 | * Perfroms a rank-1 update by outer product of vectors: \f$V = V + w v v'\f$. |
---|
[22] | 42 | * @param v Vector forming the outer product to be added |
---|
| 43 | * @param w weight of updating; can be negative |
---|
[7] | 44 | |
---|
[22] | 45 | BLAS-2b operation. |
---|
| 46 | */ |
---|
| 47 | virtual void opupdt ( const vec &v, double w ) =0; |
---|
[7] | 48 | |
---|
[22] | 49 | /*! \brief Conversion to full matrix. |
---|
| 50 | */ |
---|
[7] | 51 | |
---|
[168] | 52 | virtual mat to_mat() const =0; |
---|
[7] | 53 | |
---|
[75] | 54 | /*! \brief Inplace symmetric multiplication by a SQUARE matrix \f$C\f$, i.e. \f$V = C*V*C'\f$ |
---|
[22] | 55 | @param C multiplying matrix, |
---|
| 56 | */ |
---|
[26] | 57 | virtual void mult_sym ( const mat &C ) =0; |
---|
| 58 | |
---|
[75] | 59 | /*! \brief Inplace symmetric multiplication by a SQUARE transpose of matrix \f$C\f$, i.e. \f$V = C'*V*C\f$ |
---|
[26] | 60 | @param C multiplying matrix, |
---|
| 61 | */ |
---|
| 62 | virtual void mult_sym_t ( const mat &C ) =0; |
---|
[7] | 63 | |
---|
| 64 | |
---|
[22] | 65 | /*! |
---|
| 66 | \brief Logarithm of a determinant. |
---|
[7] | 67 | |
---|
[22] | 68 | */ |
---|
[26] | 69 | virtual double logdet() const =0; |
---|
[22] | 70 | |
---|
| 71 | /*! |
---|
[75] | 72 | \brief Multiplies square root of \f$V\f$ by vector \f$x\f$. |
---|
[22] | 73 | |
---|
| 74 | Used e.g. in generating normal samples. |
---|
| 75 | */ |
---|
[32] | 76 | virtual vec sqrt_mult (const vec &v ) const =0; |
---|
[22] | 77 | |
---|
| 78 | /*! |
---|
[75] | 79 | \brief Evaluates quadratic form \f$x= v'*V*v\f$; |
---|
[22] | 80 | |
---|
| 81 | */ |
---|
[32] | 82 | virtual double qform (const vec &v ) const =0; |
---|
[22] | 83 | |
---|
[75] | 84 | /*! |
---|
| 85 | \brief Evaluates quadratic form \f$x= v'*inv(V)*v\f$; |
---|
| 86 | |
---|
| 87 | */ |
---|
| 88 | virtual double invqform (const vec &v ) const =0; |
---|
| 89 | |
---|
[22] | 90 | // //! easy version of the |
---|
[7] | 91 | // sqmat inv(); |
---|
| 92 | |
---|
[22] | 93 | //! Clearing matrix so that it corresponds to zeros. |
---|
| 94 | virtual void clear() =0; |
---|
[7] | 95 | |
---|
[22] | 96 | //! Reimplementing common functions of mat: cols(). |
---|
| 97 | int cols() const {return dim;}; |
---|
[7] | 98 | |
---|
[427] | 99 | //! Reimplementing common functions of mat: rows(). |
---|
[22] | 100 | int rows() const {return dim;}; |
---|
| 101 | |
---|
[28] | 102 | //! Destructor for future use; |
---|
| 103 | virtual ~sqmat(){}; |
---|
[32] | 104 | //! Default constructor |
---|
| 105 | sqmat(const int dim0): dim(dim0){}; |
---|
[270] | 106 | //! Default constructor |
---|
| 107 | sqmat(): dim(0){}; |
---|
[22] | 108 | protected: |
---|
[33] | 109 | //! dimension of the square matrix |
---|
[22] | 110 | int dim; |
---|
[7] | 111 | }; |
---|
| 112 | |
---|
| 113 | |
---|
| 114 | /*! \brief Fake sqmat. This class maps sqmat operations to operations on full matrix. |
---|
| 115 | |
---|
| 116 | This class can be used to compare performance of algorithms using decomposed matrices with perormance of the same algorithms using full matrices; |
---|
| 117 | */ |
---|
[22] | 118 | class fsqmat: public sqmat |
---|
| 119 | { |
---|
| 120 | protected: |
---|
[33] | 121 | //! Full matrix on which the operations are performed |
---|
[22] | 122 | mat M; |
---|
| 123 | public: |
---|
| 124 | void opupdt ( const vec &v, double w ); |
---|
[168] | 125 | mat to_mat() const; |
---|
[26] | 126 | void mult_sym ( const mat &C); |
---|
| 127 | void mult_sym_t ( const mat &C); |
---|
[75] | 128 | //! store result of \c mult_sym in external matrix \f$U\f$ |
---|
[32] | 129 | void mult_sym ( const mat &C, fsqmat &U) const; |
---|
[75] | 130 | //! store result of \c mult_sym_t in external matrix \f$U\f$ |
---|
[32] | 131 | void mult_sym_t ( const mat &C, fsqmat &U) const; |
---|
[22] | 132 | void clear(); |
---|
[7] | 133 | |
---|
[26] | 134 | //! Default initialization |
---|
[270] | 135 | fsqmat(){}; // mat will be initialized OK |
---|
[26] | 136 | //! Default initialization with proper size |
---|
| 137 | fsqmat(const int dim0); // mat will be initialized OK |
---|
[22] | 138 | //! Constructor |
---|
| 139 | fsqmat ( const mat &M ); |
---|
[37] | 140 | //! Constructor |
---|
[183] | 141 | fsqmat ( const fsqmat &M, const ivec &perm ):sqmat(M.rows()){it_error("not implemneted");}; |
---|
| 142 | //! Constructor |
---|
[37] | 143 | fsqmat ( const vec &d ):sqmat(d.length()){M=diag(d);}; |
---|
[22] | 144 | |
---|
[28] | 145 | //! Destructor for future use; |
---|
| 146 | virtual ~fsqmat(){}; |
---|
| 147 | |
---|
| 148 | |
---|
[22] | 149 | /*! \brief Matrix inversion preserving the chosen form. |
---|
| 150 | |
---|
| 151 | @param Inv a space where the inverse is stored. |
---|
| 152 | |
---|
| 153 | */ |
---|
[433] | 154 | void inv ( fsqmat &Inv ) const; |
---|
[22] | 155 | |
---|
[26] | 156 | double logdet() const {return log ( det ( M ) );}; |
---|
[32] | 157 | double qform (const vec &v ) const {return ( v* ( M*v ) );}; |
---|
[75] | 158 | double invqform (const vec &v ) const {return ( v* ( itpp::inv(M)*v ) );}; |
---|
[37] | 159 | vec sqrt_mult (const vec &v ) const {mat Ch=chol(M); return Ch*v;}; |
---|
[22] | 160 | |
---|
[85] | 161 | //! Add another matrix in fsq form with weight w |
---|
| 162 | void add ( const fsqmat &fsq2, double w=1.0 ){M+=fsq2.M;}; |
---|
| 163 | |
---|
[37] | 164 | //! Access functions |
---|
| 165 | void setD (const vec &nD){M=diag(nD);} |
---|
| 166 | //! Access functions |
---|
[51] | 167 | vec getD (){return diag(M);} |
---|
| 168 | //! Access functions |
---|
[37] | 169 | void setD (const vec &nD, int i){for(int j=i;j<nD.length();j++){M(j,j)=nD(j-i);}} //Fixme can be more general |
---|
| 170 | |
---|
[85] | 171 | |
---|
[33] | 172 | //! add another fsqmat matrix |
---|
[22] | 173 | fsqmat& operator += ( const fsqmat &A ) {M+=A.M;return *this;}; |
---|
[33] | 174 | //! subtrack another fsqmat matrix |
---|
[22] | 175 | fsqmat& operator -= ( const fsqmat &A ) {M-=A.M;return *this;}; |
---|
[33] | 176 | //! multiply by a scalar |
---|
[22] | 177 | fsqmat& operator *= ( double x ) {M*=x;return *this;}; |
---|
| 178 | // fsqmat& operator = ( const fsqmat &A) {M=A.M; return *this;}; |
---|
[33] | 179 | //! print full matrix |
---|
[32] | 180 | friend std::ostream &operator<< ( std::ostream &os, const fsqmat &sq ); |
---|
[26] | 181 | |
---|
[7] | 182 | }; |
---|
| 183 | |
---|
[180] | 184 | /*! \brief Matrix stored in LD form, (commonly known as UD) |
---|
[33] | 185 | |
---|
[75] | 186 | Matrix is decomposed as follows: \f[M = L'DL\f] where only \f$L\f$ and \f$D\f$ matrices are stored. |
---|
[33] | 187 | All inplace operations modifies only these and the need to compose and decompose the matrix is avoided. |
---|
| 188 | */ |
---|
[427] | 189 | class ldmat: public sqmat |
---|
[22] | 190 | { |
---|
| 191 | public: |
---|
| 192 | //! Construct by copy of L and D. |
---|
| 193 | ldmat ( const mat &L, const vec &D ); |
---|
| 194 | //! Construct by decomposition of full matrix V. |
---|
[26] | 195 | ldmat (const mat &V ); |
---|
[180] | 196 | //! Construct by restructuring of V0 accordint to permutation vector perm. |
---|
| 197 | ldmat (const ldmat &V0, const ivec &perm):sqmat(V0.rows()){ ldform(V0.L.get_cols(perm), V0.D);}; |
---|
[22] | 198 | //! Construct diagonal matrix with diagonal D0 |
---|
| 199 | ldmat ( vec D0 ); |
---|
[26] | 200 | //!Default constructor |
---|
[22] | 201 | ldmat (); |
---|
[26] | 202 | //! Default initialization with proper size |
---|
| 203 | ldmat(const int dim0); |
---|
[180] | 204 | |
---|
[28] | 205 | //! Destructor for future use; |
---|
| 206 | virtual ~ldmat(){}; |
---|
| 207 | |
---|
[22] | 208 | // Reimplementation of compulsory operatios |
---|
[7] | 209 | |
---|
[22] | 210 | void opupdt ( const vec &v, double w ); |
---|
[168] | 211 | mat to_mat() const; |
---|
[26] | 212 | void mult_sym ( const mat &C); |
---|
| 213 | void mult_sym_t ( const mat &C); |
---|
[33] | 214 | //! Add another matrix in LD form with weight w |
---|
[22] | 215 | void add ( const ldmat &ld2, double w=1.0 ); |
---|
[26] | 216 | double logdet() const; |
---|
[32] | 217 | double qform (const vec &v ) const; |
---|
[75] | 218 | double invqform (const vec &v ) const; |
---|
[22] | 219 | void clear(); |
---|
[32] | 220 | vec sqrt_mult ( const vec &v ) const; |
---|
[7] | 221 | |
---|
[180] | 222 | |
---|
[22] | 223 | /*! \brief Matrix inversion preserving the chosen form. |
---|
| 224 | @param Inv a space where the inverse is stored. |
---|
| 225 | */ |
---|
[433] | 226 | void inv ( ldmat &Inv ) const; |
---|
[8] | 227 | |
---|
[75] | 228 | /*! \brief Symmetric multiplication of \f$U\f$ by a general matrix \f$C\f$, result of which is stored in the current class. |
---|
[33] | 229 | @param C matrix to multiply with |
---|
[26] | 230 | @param U a space where the inverse is stored. |
---|
[22] | 231 | */ |
---|
[32] | 232 | void mult_sym ( const mat &C, ldmat &U) const; |
---|
[12] | 233 | |
---|
[75] | 234 | /*! \brief Symmetric multiplication of \f$U\f$ by a transpose of a general matrix \f$C\f$, result of which is stored in the current class. |
---|
[33] | 235 | @param C matrix to multiply with |
---|
[26] | 236 | @param U a space where the inverse is stored. |
---|
| 237 | */ |
---|
[32] | 238 | void mult_sym_t ( const mat &C, ldmat &U) const; |
---|
[26] | 239 | |
---|
| 240 | |
---|
[75] | 241 | /*! \brief Transforms general \f$A'D0 A\f$ into pure \f$L'DL\f$ |
---|
[26] | 242 | |
---|
[75] | 243 | The new decomposition fullfills: \f$A'*diag(D)*A = self.L'*diag(self.D)*self.L\f$ |
---|
[22] | 244 | @param A general matrix |
---|
| 245 | @param D0 general vector |
---|
| 246 | */ |
---|
[26] | 247 | void ldform (const mat &A,const vec &D0 ); |
---|
[22] | 248 | |
---|
[32] | 249 | //! Access functions |
---|
| 250 | void setD (const vec &nD){D=nD;} |
---|
[33] | 251 | //! Access functions |
---|
| 252 | void setD (const vec &nD, int i){D.replace_mid(i,nD);} //Fixme can be more general |
---|
| 253 | //! Access functions |
---|
[32] | 254 | void setL (const vec &nL){L=nL;} |
---|
| 255 | |
---|
[98] | 256 | //! Access functions |
---|
| 257 | const vec& _D() const {return D;} |
---|
| 258 | //! Access functions |
---|
| 259 | const mat& _L() const {return L;} |
---|
| 260 | |
---|
[33] | 261 | //! add another ldmat matrix |
---|
[22] | 262 | ldmat& operator += ( const ldmat &ldA ); |
---|
[33] | 263 | //! subtract another ldmat matrix |
---|
[22] | 264 | ldmat& operator -= ( const ldmat &ldA ); |
---|
[33] | 265 | //! multiply by a scalar |
---|
[22] | 266 | ldmat& operator *= ( double x ); |
---|
| 267 | |
---|
[33] | 268 | //! print both \c L and \c D |
---|
[32] | 269 | friend std::ostream &operator<< ( std::ostream &os, const ldmat &sq ); |
---|
[22] | 270 | |
---|
| 271 | protected: |
---|
[75] | 272 | //! Positive vector \f$D\f$ |
---|
[22] | 273 | vec D; |
---|
[75] | 274 | //! Lower-triangular matrix \f$L\f$ |
---|
[22] | 275 | mat L; |
---|
| 276 | |
---|
[7] | 277 | }; |
---|
| 278 | |
---|
| 279 | //////// Operations: |
---|
[33] | 280 | //!mapping of add operation to operators |
---|
[22] | 281 | inline ldmat& ldmat::operator += ( const ldmat &ldA ) {this->add ( ldA );return *this;} |
---|
[33] | 282 | //!mapping of negative add operation to operators |
---|
[22] | 283 | inline ldmat& ldmat::operator -= ( const ldmat &ldA ) {this->add ( ldA,-1.0 );return *this;} |
---|
[7] | 284 | |
---|
[219] | 285 | /*! @} */ |
---|
| 286 | |
---|
[7] | 287 | #endif // DC_H |
---|