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