ROOT   Reference Guide
Searching...
No Matches
TDecompChol.cxx
Go to the documentation of this file.
1// @(#)root/matrix:$Id$
2// Authors: Fons Rademakers, Eddy Offermann Dec 2003
3
4/*************************************************************************
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. * 9 * For the list of contributors see$ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12/** \class TDecompChol
13 \ingroup Matrix
14
15 Cholesky Decomposition class
16
17 Decompose a symmetric, positive definite matrix A = U^T * U
18
19 where U is a upper triangular matrix
20
21 The decomposition fails if a diagonal element of fU is <= 0, the
22 matrix is not positive negative . The matrix fU is made invalid .
23
24 fU has the same index range as A .
25*/
26
27#include "TDecompChol.h"
28#include "TMath.h"
29
31
32////////////////////////////////////////////////////////////////////////////////
33/// Constructor for (nrows x nrows) matrix
34
36{
37 fU.ResizeTo(nrows,nrows);
38}
39
40////////////////////////////////////////////////////////////////////////////////
41/// Constructor for ([row_lwb..row_upb] x [row_lwb..row_upb]) matrix
42
44{
45 const Int_t nrows = row_upb-row_lwb+1;
46 fRowLwb = row_lwb;
47 fColLwb = row_lwb;
48 fU.ResizeTo(row_lwb,row_lwb+nrows-1,row_lwb,row_lwb+nrows-1);
49}
50
51////////////////////////////////////////////////////////////////////////////////
52/// Constructor for symmetric matrix A . Matrix should be positive definite
53
55{
56 R__ASSERT(a.IsValid());
57
59 fCondition = a.Norm1();
60 fTol = a.GetTol();
61 if (tol > 0)
62 fTol = tol;
63
64 fRowLwb = a.GetRowLwb();
65 fColLwb = a.GetColLwb();
66 fU.ResizeTo(a);
67 fU = a;
68}
69
70////////////////////////////////////////////////////////////////////////////////
71/// Constructor for general matrix A . Matrix should be symmetric positive definite
72
74{
75 R__ASSERT(a.IsValid());
76
77 if (a.GetNrows() != a.GetNcols() || a.GetRowLwb() != a.GetColLwb()) {
78 Error("TDecompChol(const TMatrixD &","matrix should be square");
79 return;
80 }
81
83 fCondition = a.Norm1();
84 fTol = a.GetTol();
85 if (tol > 0)
86 fTol = tol;
87
88 fRowLwb = a.GetRowLwb();
89 fColLwb = a.GetColLwb();
90 fU.ResizeTo(a);
91 fU = a;
92}
93
94////////////////////////////////////////////////////////////////////////////////
95/// Copy constructor
96
98{
99 *this = another;
100}
101
102////////////////////////////////////////////////////////////////////////////////
103/// Matrix A is decomposed in component U so that A = U^T * U
104/// If the decomposition succeeds, bit kDecomposed is set , otherwise kSingular
105
107{
108 if (TestBit(kDecomposed)) return kTRUE;
109
110 if ( !TestBit(kMatrixSet) ) {
111 Error("Decompose()","Matrix has not been set");
112 return kFALSE;
113 }
114
115 Int_t i,j,icol,irow;
116 const Int_t n = fU.GetNrows();
117 Double_t *pU = fU.GetMatrixArray();
118 for (icol = 0; icol < n; icol++) {
119 const Int_t rowOff = icol*n;
120
121 //Compute fU(j,j) and test for non-positive-definiteness.
122 Double_t ujj = pU[rowOff+icol];
123 for (irow = 0; irow < icol; irow++) {
124 const Int_t pos_ij = irow*n+icol;
125 ujj -= pU[pos_ij]*pU[pos_ij];
126 }
127 if (ujj <= 0) {
128 Error("Decompose()","matrix not positive definite");
129 return kFALSE;
130 }
131 ujj = TMath::Sqrt(ujj);
132 pU[rowOff+icol] = ujj;
133 const Double_t inv_ujj = 1.0 / ujj;
134
135 if (icol < n-1) {
136 for (j = icol+1; j < n; j++) {
137 for (i = 0; i < icol; i++) {
138 const Int_t rowOff2 = i*n;
139 pU[rowOff+j] -= pU[rowOff2+j]*pU[rowOff2+icol];
140 }
141 }
142 for (j = icol+1; j < n; j++)
143 pU[rowOff + j] *= inv_ujj;
144 }
145 }
146
147 for (irow = 0; irow < n; irow++) {
148 const Int_t rowOff = irow*n;
149 for (icol = 0; icol < irow; icol++)
150 pU[rowOff+icol] = 0.;
151 }
152
154
155 return kTRUE;
156}
157
158////////////////////////////////////////////////////////////////////////////////
159/// Reconstruct the original matrix using the decomposition parts
160
162{
163 if (TestBit(kSingular)) {
164 Error("GetMatrix()","Matrix is singular");
165 return TMatrixDSym();
166 }
167 if ( !TestBit(kDecomposed) ) {
168 if (!Decompose()) {
169 Error("GetMatrix()","Decomposition failed");
170 return TMatrixDSym();
171 }
172 }
173
175}
176
177////////////////////////////////////////////////////////////////////////////////
178/// Set the matrix to be decomposed, decomposition status is reset.
179
181{
182 R__ASSERT(a.IsValid());
183
184 ResetStatus();
185 if (a.GetNrows() != a.GetNcols() || a.GetRowLwb() != a.GetColLwb()) {
186 Error("SetMatrix(const TMatrixDSym &","matrix should be square");
187 return;
188 }
189
191 fCondition = -1.0;
192
193 fRowLwb = a.GetRowLwb();
194 fColLwb = a.GetColLwb();
195 fU.ResizeTo(a);
196 fU = a;
197}
198
199////////////////////////////////////////////////////////////////////////////////
200/// Solve equations Ax=b assuming A has been factored by Cholesky. The factor U is
201/// assumed to be in upper triang of fU. fTol is used to determine if diagonal
202/// element is zero. The solution is returned in b.
203
205{
206 R__ASSERT(b.IsValid());
207 if (TestBit(kSingular)) {
208 Error("Solve()","Matrix is singular");
209 return kFALSE;
210 }
211 if ( !TestBit(kDecomposed) ) {
212 if (!Decompose()) {
213 Error("Solve()","Decomposition failed");
214 return kFALSE;
215 }
216 }
217
218 if (fU.GetNrows() != b.GetNrows() || fU.GetRowLwb() != b.GetLwb()) {
219 Error("Solve(TVectorD &","vector and matrix incompatible");
220 return kFALSE;
221 }
222
223 const Int_t n = fU.GetNrows();
224
225 const Double_t *pU = fU.GetMatrixArray();
226 Double_t *pb = b.GetMatrixArray();
227
228 Int_t i;
229 // step 1: Forward substitution on U^T
230 for (i = 0; i < n; i++) {
231 const Int_t off_i = i*n;
232 if (pU[off_i+i] < fTol)
233 {
234 Error("Solve(TVectorD &b)","u[%d,%d]=%.4e < %.4e",i,i,pU[off_i+i],fTol);
235 return kFALSE;
236 }
237 Double_t r = pb[i];
238 for (Int_t j = 0; j < i; j++) {
239 const Int_t off_j = j*n;
240 r -= pU[off_j+i]*pb[j];
241 }
242 pb[i] = r/pU[off_i+i];
243 }
244
245 // step 2: Backward substitution on U
246 for (i = n-1; i >= 0; i--) {
247 const Int_t off_i = i*n;
248 Double_t r = pb[i];
249 for (Int_t j = i+1; j < n; j++)
250 r -= pU[off_i+j]*pb[j];
251 pb[i] = r/pU[off_i+i];
252 }
253
254 return kTRUE;
255}
256
257////////////////////////////////////////////////////////////////////////////////
258/// Solve equations Ax=b assuming A has been factored by Cholesky. The factor U is
259/// assumed to be in upper triang of fU. fTol is used to determine if diagonal
260/// element is zero. The solution is returned in b.
261
263{
264 TMatrixDBase *b = const_cast<TMatrixDBase *>(cb.GetMatrix());
265 R__ASSERT(b->IsValid());
266 if (TestBit(kSingular)) {
267 Error("Solve()","Matrix is singular");
268 return kFALSE;
269 }
270 if ( !TestBit(kDecomposed) ) {
271 if (!Decompose()) {
272 Error("Solve()","Decomposition failed");
273 return kFALSE;
274 }
275 }
276
277 if (fU.GetNrows() != b->GetNrows() || fU.GetRowLwb() != b->GetRowLwb())
278 {
279 Error("Solve(TMatrixDColumn &cb","vector and matrix incompatible");
280 return kFALSE;
281 }
282
283 const Int_t n = fU.GetNrows();
284
285 const Double_t *pU = fU.GetMatrixArray();
286 Double_t *pcb = cb.GetPtr();
287 const Int_t inc = cb.GetInc();
288
289 Int_t i;
290 // step 1: Forward substitution U^T
291 for (i = 0; i < n; i++) {
292 const Int_t off_i = i*n;
293 const Int_t off_i2 = i*inc;
294 if (pU[off_i+i] < fTol)
295 {
296 Error("Solve(TMatrixDColumn &cb)","u[%d,%d]=%.4e < %.4e",i,i,pU[off_i+i],fTol);
297 return kFALSE;
298 }
299 Double_t r = pcb[off_i2];
300 for (Int_t j = 0; j < i; j++) {
301 const Int_t off_j = j*n;
302 r -= pU[off_j+i]*pcb[j*inc];
303 }
304 pcb[off_i2] = r/pU[off_i+i];
305 }
306
307 // step 2: Backward substitution U
308 for (i = n-1; i >= 0; i--) {
309 const Int_t off_i = i*n;
310 const Int_t off_i2 = i*inc;
311 Double_t r = pcb[off_i2];
312 for (Int_t j = i+1; j < n; j++)
313 r -= pU[off_i+j]*pcb[j*inc];
314 pcb[off_i2] = r/pU[off_i+i];
315 }
316
317 return kTRUE;
318}
319
320////////////////////////////////////////////////////////////////////////////////
321/// Matrix determinant det = d1*TMath::Power(2.,d2) is square of diagProd
322/// of cholesky factor
323
325{
326 if ( !TestBit(kDetermined) ) {
327 if ( !TestBit(kDecomposed) )
328 Decompose();
329 TDecompBase::Det(d1,d2);
330 // square det as calculated by above
331 fDet1 *= fDet1;
332 fDet2 += fDet2;
334 }
335 d1 = fDet1;
336 d2 = fDet2;
337}
338
339////////////////////////////////////////////////////////////////////////////////
340/// For a symmetric matrix A(m,m), its inverse A_inv(m,m) is returned .
341
343{
344 if (inv.GetNrows() != GetNrows() || inv.GetRowLwb() != GetRowLwb()) {
345 Error("Invert(TMatrixDSym &","Input matrix has wrong shape");
346 return kFALSE;
347 }
348
349 inv.UnitMatrix();
350
351 const Int_t colLwb = inv.GetColLwb();
352 const Int_t colUpb = inv.GetColUpb();
353 Bool_t status = kTRUE;
354 for (Int_t icol = colLwb; icol <= colUpb && status; icol++) {
355 TMatrixDColumn b(inv,icol);
356 status &= Solve(b);
357 }
358
359 return status;
360}
361
362////////////////////////////////////////////////////////////////////////////////
363/// For a symmetric matrix A(m,m), its inverse A_inv(m,m) is returned .
364
366{
367 const Int_t rowLwb = GetRowLwb();
368 const Int_t rowUpb = rowLwb+GetNrows()-1;
369
370 TMatrixDSym inv(rowLwb,rowUpb);
371 inv.UnitMatrix();
372 status = Invert(inv);
373
374 return inv;
375}
376
377////////////////////////////////////////////////////////////////////////////////
378/// Print class members .
379
381{
383 fU.Print("fU");
384}
385
386////////////////////////////////////////////////////////////////////////////////
387/// Assignment operator
388
390{
391 if (this != &source) {
393 fU.ResizeTo(source.fU);
394 fU = source.fU;
395 }
396 return *this;
397}
398
399////////////////////////////////////////////////////////////////////////////////
400/// Solve min {(A . x - b)^T (A . x - b)} for vector x where
401/// A : (m x n) matrix, m >= n
402/// b : (m) vector
403/// x : (n) vector
404
406{
408 Bool_t ok;
409 return ch.Solve(TMatrixD(TMatrixD::kTransposed,A)*b,ok);
410}
411
412////////////////////////////////////////////////////////////////////////////////
413/// Solve min {(A . x - b)^T W (A . x - b)} for vector x where
414/// A : (m x n) matrix, m >= n
415/// b : (m) vector
416/// x : (n) vector
417/// W : (m x m) weight matrix with W(i,j) = 1/std(i)^2 for i == j
418/// = 0 for i != j
419
420TVectorD NormalEqn(const TMatrixD &A,const TVectorD &b,const TVectorD &std)
421{
422 if (!AreCompatible(b,std)) {
423 ::Error("NormalEqn","vectors b and std are not compatible");
424 return TVectorD();
425 }
426
427 TMatrixD mAw = A;
428 TVectorD mBw = b;
429 for (Int_t irow = 0; irow < A.GetNrows(); irow++) {
430 TMatrixDRow(mAw,irow) *= 1/std(irow);
431 mBw(irow) /= std(irow);
432 }
434 Bool_t ok;
435 return ch.Solve(TMatrixD(TMatrixD::kTransposed,mAw)*mBw,ok);
436}
437
438////////////////////////////////////////////////////////////////////////////////
439/// Solve min {(A . X_j - B_j)^T (A . X_j - B_j)} for each column j in
440/// B and X
441/// A : (m x n ) matrix, m >= n
442/// B : (m x nb) matrix, nb >= 1
443/// mX : (n x nb) matrix
444
446{
449 ch.MultiSolve(mX);
450 return mX;
451}
452
453////////////////////////////////////////////////////////////////////////////////
454/// Solve min {(A . X_j - B_j)^T W (A . X_j - B_j)} for each column j in
455/// B and X
456/// ~~~
457/// A : (m x n ) matrix, m >= n
458/// B : (m x nb) matrix, nb >= 1
459/// mX : (n x nb) matrix
460/// W : (m x m) weight matrix with W(i,j) = 1/std(i)^2 for i == j
461/// = 0 for i != j
462/// ~~~
463
464TMatrixD NormalEqn(const TMatrixD &A,const TMatrixD &B,const TVectorD &std)
465{
466 TMatrixD mAw = A;
467 TMatrixD mBw = B;
468 for (Int_t irow = 0; irow < A.GetNrows(); irow++) {
469 TMatrixDRow(mAw,irow) *= 1/std(irow);
470 TMatrixDRow(mBw,irow) *= 1/std(irow);
471 }
472
475 ch.MultiSolve(mX);
476 return mX;
477}
#define b(i)
Definition RSha256.hxx:100
#define a(i)
Definition RSha256.hxx:99
bool Bool_t
Definition RtypesCore.h:63
int Int_t
Definition RtypesCore.h:45
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
constexpr Bool_t kTRUE
Definition RtypesCore.h:100
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:377
TVectorD NormalEqn(const TMatrixD &A, const TVectorD &b)
Solve min {(A .
#define R__ASSERT(e)
Definition TError.h:118
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:185
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
TMatrixTSym< Double_t > TMatrixDSym
TMatrixTRow< Double_t > TMatrixDRow
TMatrixT< Double_t > TMatrixD
Definition TMatrixDfwd.h:23
TVectorT< Double_t > TVectorD
Definition TVectorDfwd.h:23
Decomposition Base class.
Definition TDecompBase.h:34
Int_t GetRowLwb() const
Definition TDecompBase.h:73
Double_t fDet1
Definition TDecompBase.h:37
Double_t fDet2
Definition TDecompBase.h:38
void ResetStatus()
Definition TDecompBase.h:43
virtual Bool_t MultiSolve(TMatrixD &B)
Solve set of equations with RHS in columns of B.
Int_t fRowLwb
Definition TDecompBase.h:40
Double_t fTol
Definition TDecompBase.h:36
TDecompBase & operator=(const TDecompBase &source)
Assignment operator.
Double_t fCondition
Definition TDecompBase.h:39
Int_t fColLwb
Definition TDecompBase.h:41
void Print(Option_t *opt="") const override
Print class members.
virtual void Det(Double_t &d1, Double_t &d2)
Matrix determinant det = d1*TMath::Power(2.,d2)
Cholesky Decomposition class.
Definition TDecompChol.h:25
void Print(Option_t *opt="") const override
Print class members .
Bool_t Decompose() override
Matrix A is decomposed in component U so that A = U^T * U If the decomposition succeeds,...
TMatrixD fU
Definition TDecompChol.h:28
Int_t GetNrows() const override
Definition TDecompChol.h:43
virtual void SetMatrix(const TMatrixDSym &a)
Set the matrix to be decomposed, decomposition status is reset.
Bool_t Solve(TVectorD &b) override
Solve equations Ax=b assuming A has been factored by Cholesky.
TDecompChol & operator=(const TDecompChol &source)
Assignment operator.
const TMatrixDSym GetMatrix()
Reconstruct the original matrix using the decomposition parts.
void Det(Double_t &d1, Double_t &d2) override
Matrix determinant det = d1*TMath::Power(2.,d2) is square of diagProd of cholesky factor.
TMatrixDSym Invert()
Definition TDecompChol.h:60
Int_t GetNrows() const
void Print(Option_t *name="") const override
Print the matrix as a table of elements.
Int_t GetRowLwb() const
const TMatrixTBase< Element > * GetMatrix() const
Element * GetPtr() const
const Element * GetMatrixArray() const override
Definition TMatrixT.h:225
TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1) override
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition TObject.h:199
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:780
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:987
const Int_t n
Definition legend1.C:16
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:662
Bool_t AreCompatible(const TMatrixTBase< Element1 > &m1, const TMatrixTBase< Element2 > &m2, Int_t verbose=0)
Check that matrice sm1 and m2 areboth valid and have identical shapes .
void inv(rsa_NUMBER *, rsa_NUMBER *, rsa_NUMBER *)
Definition rsaaux.cxx:949