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