[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" |
---|
[565] | 18 | #include "../bdmerror.h" |
---|
[7] | 19 | |
---|
[737] | 20 | namespace bdm { |
---|
[495] | 21 | |
---|
[7] | 22 | using namespace itpp; |
---|
| 23 | |
---|
[37] | 24 | //! Auxiliary function dydr; dyadic reduction |
---|
[477] | 25 | void dydr ( double * r, double *f, double *Dr, double *Df, double *R, int jl, int jh, double *kr, int m, int mx ); |
---|
[37] | 26 | |
---|
| 27 | //! Auxiliary function ltuinv; inversion of a triangular matrix; |
---|
| 28 | //TODO can be done via: dtrtri.f from lapack |
---|
[477] | 29 | mat ltuinv ( const mat &L ); |
---|
[37] | 30 | |
---|
[427] | 31 | /*! \brief Abstract class for representation of double symmetric matrices in square-root form. |
---|
[7] | 32 | |
---|
[37] | 33 | All operations defined on this class should be optimized for the chosen decomposition. |
---|
[7] | 34 | */ |
---|
[477] | 35 | class sqmat { |
---|
| 36 | public: |
---|
[1064] | 37 | /*! |
---|
| 38 | * Perfroms a rank-1 update by outer product of vectors: \f$V = V + w v v'\f$. |
---|
| 39 | * @param v Vector forming the outer product to be added |
---|
| 40 | * @param w weight of updating; can be negative |
---|
[7] | 41 | |
---|
[1064] | 42 | BLAS-2b operation. |
---|
| 43 | */ |
---|
| 44 | virtual void opupdt ( const vec &v, double w ) = 0; |
---|
[7] | 45 | |
---|
[1064] | 46 | /*! \brief Conversion to full matrix. |
---|
| 47 | */ |
---|
| 48 | virtual mat to_mat() const = 0; |
---|
[7] | 49 | |
---|
[1064] | 50 | /*! \brief Inplace symmetric multiplication by a SQUARE matrix \f$C\f$, i.e. \f$V = C*V*C'\f$ |
---|
| 51 | @param C multiplying matrix, |
---|
| 52 | */ |
---|
| 53 | virtual void mult_sym ( const mat &C ) = 0; |
---|
[7] | 54 | |
---|
[1064] | 55 | /*! \brief Inplace symmetric multiplication by a SQUARE transpose of matrix \f$C\f$, i.e. \f$V = C'*V*C\f$ |
---|
| 56 | @param C multiplying matrix, |
---|
| 57 | */ |
---|
| 58 | virtual void mult_sym_t ( const mat &C ) = 0; |
---|
[7] | 59 | |
---|
[1064] | 60 | /*! |
---|
| 61 | \brief Logarithm of a determinant. |
---|
[22] | 62 | |
---|
[1064] | 63 | */ |
---|
| 64 | virtual double logdet() const = 0; |
---|
[22] | 65 | |
---|
[1064] | 66 | /*! |
---|
| 67 | \brief Multiplies square root of \f$V\f$ by vector \f$x\f$. |
---|
[22] | 68 | |
---|
[1064] | 69 | Used e.g. in generating normal samples. |
---|
| 70 | */ |
---|
| 71 | virtual vec sqrt_mult ( const vec &v ) const = 0; |
---|
[22] | 72 | |
---|
[1064] | 73 | /*! |
---|
| 74 | \brief Evaluates quadratic form \f$x= v'*V*v\f$; |
---|
[22] | 75 | |
---|
[1064] | 76 | */ |
---|
| 77 | virtual double qform ( const vec &v ) const = 0; |
---|
[75] | 78 | |
---|
[1064] | 79 | /*! |
---|
| 80 | \brief Evaluates quadratic form \f$x= v'*inv(V)*v\f$; |
---|
[75] | 81 | |
---|
[1064] | 82 | */ |
---|
| 83 | virtual double invqform ( const vec &v ) const = 0; |
---|
[477] | 84 | |
---|
[22] | 85 | // //! easy version of the |
---|
[7] | 86 | // sqmat inv(); |
---|
| 87 | |
---|
[1064] | 88 | //! Clearing matrix so that it corresponds to zeros. |
---|
| 89 | virtual void clear() = 0; |
---|
[7] | 90 | |
---|
[1064] | 91 | //! Reimplementing common functions of mat: cols(). |
---|
| 92 | int cols() const { |
---|
| 93 | return dim; |
---|
| 94 | }; |
---|
[7] | 95 | |
---|
[1064] | 96 | //! Reimplementing common functions of mat: rows(). |
---|
| 97 | int rows() const { |
---|
| 98 | return dim; |
---|
| 99 | }; |
---|
[22] | 100 | |
---|
[1064] | 101 | //! Destructor for future use; |
---|
| 102 | virtual ~sqmat() {}; |
---|
| 103 | //! Default constructor |
---|
| 104 | sqmat ( const int dim0 ) : dim ( dim0 ) {}; |
---|
| 105 | //! Default constructor |
---|
| 106 | sqmat() : dim ( 0 ) {}; |
---|
[477] | 107 | protected: |
---|
[1064] | 108 | //! dimension of the square matrix |
---|
| 109 | int dim; |
---|
[7] | 110 | }; |
---|
| 111 | |
---|
| 112 | |
---|
| 113 | /*! \brief Fake sqmat. This class maps sqmat operations to operations on full matrix. |
---|
| 114 | |
---|
| 115 | This class can be used to compare performance of algorithms using decomposed matrices with perormance of the same algorithms using full matrices; |
---|
| 116 | */ |
---|
[477] | 117 | class fsqmat: public sqmat { |
---|
| 118 | protected: |
---|
[1064] | 119 | //! Full matrix on which the operations are performed |
---|
| 120 | mat M; |
---|
[477] | 121 | public: |
---|
[1064] | 122 | void opupdt ( const vec &v, double w ); |
---|
| 123 | mat to_mat() const; |
---|
| 124 | void mult_sym ( const mat &C ); |
---|
| 125 | void mult_sym_t ( const mat &C ); |
---|
| 126 | //! store result of \c mult_sym in external matrix \f$U\f$ |
---|
| 127 | void mult_sym ( const mat &C, fsqmat &U ) const; |
---|
| 128 | //! store result of \c mult_sym_t in external matrix \f$U\f$ |
---|
| 129 | void mult_sym_t ( const mat &C, fsqmat &U ) const; |
---|
| 130 | void clear(); |
---|
[7] | 131 | |
---|
[1064] | 132 | //! Default initialization |
---|
| 133 | fsqmat() {}; // mat will be initialized OK |
---|
| 134 | //! Default initialization with proper size |
---|
| 135 | fsqmat ( const int dim0 ); // mat will be initialized OK |
---|
| 136 | //! Constructor |
---|
| 137 | fsqmat ( const mat &M ); |
---|
[565] | 138 | |
---|
[1064] | 139 | /*! |
---|
| 140 | Some templates require this constructor to compile, but |
---|
| 141 | it shouldn't actually be called. |
---|
| 142 | */ |
---|
| 143 | fsqmat ( const fsqmat &M, const ivec &perm ) { |
---|
| 144 | bdm_error ( "not implemented" ); |
---|
| 145 | } |
---|
[565] | 146 | |
---|
[1064] | 147 | //! Constructor |
---|
| 148 | fsqmat ( const vec &d ) : sqmat ( d.length() ) { |
---|
| 149 | M = diag ( d ); |
---|
| 150 | }; |
---|
[22] | 151 | |
---|
[1064] | 152 | //! Destructor for future use; |
---|
| 153 | virtual ~fsqmat() {}; |
---|
[28] | 154 | |
---|
| 155 | |
---|
[1064] | 156 | /*! \brief Matrix inversion preserving the chosen form. |
---|
[22] | 157 | |
---|
[1064] | 158 | @param Inv a space where the inverse is stored. |
---|
[22] | 159 | |
---|
[1064] | 160 | */ |
---|
| 161 | void inv ( fsqmat &Inv ) const; |
---|
[22] | 162 | |
---|
[1064] | 163 | double logdet() const { |
---|
| 164 | return log ( det ( M ) ); |
---|
| 165 | }; |
---|
| 166 | double qform ( const vec &v ) const { |
---|
| 167 | return ( v* ( M*v ) ); |
---|
| 168 | }; |
---|
| 169 | double invqform ( const vec &v ) const { |
---|
| 170 | return ( v* ( itpp::inv ( M ) *v ) ); |
---|
| 171 | }; |
---|
| 172 | vec sqrt_mult ( const vec &v ) const { |
---|
| 173 | mat Ch = chol ( M ); |
---|
| 174 | return Ch*v; |
---|
| 175 | }; |
---|
[22] | 176 | |
---|
[1064] | 177 | //! Add another matrix in fsq form with weight w |
---|
| 178 | void add ( const fsqmat &fsq2, double w = 1.0 ) { |
---|
| 179 | M += fsq2.M; |
---|
| 180 | }; |
---|
[85] | 181 | |
---|
[1064] | 182 | //! Access functions |
---|
| 183 | void setD ( const vec &nD ) { |
---|
| 184 | M = diag ( nD ); |
---|
| 185 | } |
---|
| 186 | //! Access functions |
---|
| 187 | vec getD () { |
---|
| 188 | return diag ( M ); |
---|
| 189 | } |
---|
| 190 | //! Access functions |
---|
| 191 | void setD ( const vec &nD, int i ) { |
---|
| 192 | for ( int j = i; j < nD.length(); j++ ) { |
---|
| 193 | M ( j, j ) = nD ( j - i ); //Fixme can be more general |
---|
| 194 | } |
---|
| 195 | } |
---|
[37] | 196 | |
---|
[85] | 197 | |
---|
[1064] | 198 | //! add another fsqmat matrix |
---|
| 199 | fsqmat& operator += ( const fsqmat &A ) { |
---|
| 200 | M += A.M; |
---|
| 201 | return *this; |
---|
| 202 | }; |
---|
| 203 | //! subtrack another fsqmat matrix |
---|
| 204 | fsqmat& operator -= ( const fsqmat &A ) { |
---|
| 205 | M -= A.M; |
---|
| 206 | return *this; |
---|
| 207 | }; |
---|
| 208 | //! multiply by a scalar |
---|
| 209 | fsqmat& operator *= ( double x ) { |
---|
| 210 | M *= x; |
---|
| 211 | return *this; |
---|
| 212 | }; |
---|
[737] | 213 | |
---|
[1064] | 214 | //! cast to normal mat |
---|
| 215 | operator mat&() { |
---|
| 216 | return M; |
---|
| 217 | }; |
---|
[737] | 218 | |
---|
[22] | 219 | // fsqmat& operator = ( const fsqmat &A) {M=A.M; return *this;}; |
---|
[1064] | 220 | //! print full matrix |
---|
| 221 | friend std::ostream &operator<< ( std::ostream &os, const fsqmat &sq ); |
---|
| 222 | //!access function |
---|
| 223 | mat & _M ( ) { |
---|
| 224 | return M; |
---|
| 225 | }; |
---|
[26] | 226 | |
---|
[7] | 227 | }; |
---|
| 228 | |
---|
[583] | 229 | |
---|
[477] | 230 | /*! \brief Matrix stored in LD form, (commonly known as UD) |
---|
[33] | 231 | |
---|
[75] | 232 | Matrix is decomposed as follows: \f[M = L'DL\f] where only \f$L\f$ and \f$D\f$ matrices are stored. |
---|
[33] | 233 | All inplace operations modifies only these and the need to compose and decompose the matrix is avoided. |
---|
| 234 | */ |
---|
[477] | 235 | class ldmat: public sqmat { |
---|
| 236 | public: |
---|
[1064] | 237 | //! Construct by copy of L and D. |
---|
| 238 | ldmat ( const mat &L, const vec &D ); |
---|
| 239 | //! Construct by decomposition of full matrix V. |
---|
| 240 | ldmat ( const mat &V ); |
---|
| 241 | //! Construct by restructuring of V0 accordint to permutation vector perm. |
---|
| 242 | ldmat ( const ldmat &V0, const ivec &perm ) : sqmat ( V0.rows() ) { |
---|
| 243 | ldform ( V0.L.get_cols ( perm ), V0.D ); |
---|
| 244 | }; |
---|
| 245 | //! Construct diagonal matrix with diagonal D0 |
---|
| 246 | ldmat ( vec D0 ); |
---|
| 247 | //!Default constructor |
---|
| 248 | ldmat (); |
---|
| 249 | //! Default initialization with proper size |
---|
| 250 | ldmat ( const int dim0 ); |
---|
[28] | 251 | |
---|
[1064] | 252 | //! Destructor for future use; |
---|
| 253 | virtual ~ldmat() {}; |
---|
[7] | 254 | |
---|
[1064] | 255 | // Reimplementation of compulsory operatios |
---|
[7] | 256 | |
---|
[1064] | 257 | void opupdt ( const vec &v, double w ); |
---|
| 258 | mat to_mat() const; |
---|
| 259 | void mult_sym ( const mat &C ); |
---|
| 260 | void mult_sym_t ( const mat &C ); |
---|
| 261 | //! Add another matrix in LD form with weight w |
---|
| 262 | void add ( const ldmat &ld2, double w = 1.0 ); |
---|
| 263 | double logdet() const; |
---|
| 264 | double qform ( const vec &v ) const; |
---|
| 265 | double invqform ( const vec &v ) const; |
---|
| 266 | void clear(); |
---|
| 267 | vec sqrt_mult ( const vec &v ) const; |
---|
[8] | 268 | |
---|
[12] | 269 | |
---|
[1064] | 270 | /*! \brief Matrix inversion preserving the chosen form. |
---|
| 271 | @param Inv a space where the inverse is stored. |
---|
| 272 | */ |
---|
| 273 | void inv ( ldmat &Inv ) const; |
---|
[26] | 274 | |
---|
[1064] | 275 | /*! \brief Symmetric multiplication of \f$U\f$ by a general matrix \f$C\f$, result of which is stored in the current class. |
---|
| 276 | @param C matrix to multiply with |
---|
| 277 | @param U a space where the inverse is stored. |
---|
| 278 | */ |
---|
| 279 | void mult_sym ( const mat &C, ldmat &U ) const; |
---|
[26] | 280 | |
---|
[1064] | 281 | /*! \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. |
---|
| 282 | @param C matrix to multiply with |
---|
| 283 | @param U a space where the inverse is stored. |
---|
| 284 | */ |
---|
| 285 | void mult_sym_t ( const mat &C, ldmat &U ) const; |
---|
[26] | 286 | |
---|
[22] | 287 | |
---|
[1064] | 288 | /*! \brief Transforms general \f$A'D0 A\f$ into pure \f$L'DL\f$ |
---|
[32] | 289 | |
---|
[1064] | 290 | The new decomposition fullfills: \f$A'*diag(D)*A = self.L'*diag(self.D)*self.L\f$ |
---|
| 291 | @param A general matrix |
---|
| 292 | @param D0 general vector |
---|
| 293 | */ |
---|
| 294 | void ldform ( const mat &A, const vec &D0 ); |
---|
[98] | 295 | |
---|
[1064] | 296 | //! Access functions |
---|
| 297 | void setD ( const vec &nD ) { |
---|
| 298 | D = nD; |
---|
| 299 | } |
---|
| 300 | //! Access functions |
---|
| 301 | void setD ( const vec &nD, int i ) { |
---|
| 302 | D.replace_mid ( i, nD ); //Fixme can be more general |
---|
| 303 | } |
---|
| 304 | //! Access functions |
---|
| 305 | void setL ( const vec &nL ) { |
---|
| 306 | L = nL; |
---|
| 307 | } |
---|
[22] | 308 | |
---|
[1015] | 309 | //! Access functions |
---|
[1064] | 310 | const vec& _D() const { |
---|
| 311 | return D; |
---|
| 312 | } |
---|
[1015] | 313 | //! Access functions |
---|
[1064] | 314 | const mat& _L() const { |
---|
| 315 | return L; |
---|
| 316 | } |
---|
[1015] | 317 | //! Access functions |
---|
[1064] | 318 | vec& __D() { |
---|
| 319 | return D; |
---|
| 320 | } |
---|
[1015] | 321 | //! Access functions |
---|
[1064] | 322 | mat& __L() { |
---|
| 323 | return L; |
---|
| 324 | } |
---|
| 325 | void validate() { |
---|
[1079] | 326 | bdm_assert(L.rows()==D.length(),"Incompatible L and D in ldmat"); |
---|
[1064] | 327 | dim= L.rows(); |
---|
| 328 | } |
---|
[22] | 329 | |
---|
[1064] | 330 | //! add another ldmat matrix |
---|
| 331 | ldmat& operator += ( const ldmat &ldA ); |
---|
| 332 | //! subtract another ldmat matrix |
---|
| 333 | ldmat& operator -= ( const ldmat &ldA ); |
---|
| 334 | //! multiply by a scalar |
---|
| 335 | ldmat& operator *= ( double x ); |
---|
[22] | 336 | |
---|
[1064] | 337 | //! print both \c L and \c D |
---|
| 338 | friend std::ostream &operator<< ( std::ostream &os, const ldmat &sq ); |
---|
[477] | 339 | |
---|
| 340 | protected: |
---|
[1064] | 341 | //! Positive vector \f$D\f$ |
---|
| 342 | vec D; |
---|
| 343 | //! Lower-triangular matrix \f$L\f$ |
---|
| 344 | mat L; |
---|
[477] | 345 | |
---|
[7] | 346 | }; |
---|
| 347 | |
---|
| 348 | //////// Operations: |
---|
[33] | 349 | //!mapping of add operation to operators |
---|
[477] | 350 | inline ldmat& ldmat::operator += ( const ldmat & ldA ) { |
---|
[1064] | 351 | this->add ( ldA ); |
---|
| 352 | return *this; |
---|
[477] | 353 | } |
---|
[33] | 354 | //!mapping of negative add operation to operators |
---|
[477] | 355 | inline ldmat& ldmat::operator -= ( const ldmat & ldA ) { |
---|
[1064] | 356 | this->add ( ldA, -1.0 ); |
---|
| 357 | return *this; |
---|
[477] | 358 | } |
---|
[7] | 359 | |
---|
[495] | 360 | } |
---|
| 361 | |
---|
[7] | 362 | #endif // DC_H |
---|