ROOT   Reference Guide
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
134 if (icol < n-1) {
135 for (j = icol+1; j < n; j++) {
136 for (i = 0; i < icol; i++) {
137 const Int_t rowOff2 = i*n;
138 pU[rowOff+j] -= pU[rowOff2+j]*pU[rowOff2+icol];
139 }
140 }
141 for (j = icol+1; j < n; j++)
142 pU[rowOff+j] /= ujj;
143 }
144 }
145
146 for (irow = 0; irow < n; irow++) {
147 const Int_t rowOff = irow*n;
148 for (icol = 0; icol < irow; icol++)
149 pU[rowOff+icol] = 0.;
150 }
151
153
154 return kTRUE;
155}
156
157////////////////////////////////////////////////////////////////////////////////
158/// Reconstruct the original matrix using the decomposition parts
159
161{
162 if (TestBit(kSingular)) {
163 Error("GetMatrix()","Matrix is singular");
164 return TMatrixDSym();
165 }
166 if ( !TestBit(kDecomposed) ) {
167 if (!Decompose()) {
168 Error("GetMatrix()","Decomposition failed");
169 return TMatrixDSym();
170 }
171 }
172
174}
175
176////////////////////////////////////////////////////////////////////////////////
177/// Set the matrix to be decomposed, decomposition status is reset.
178
180{
181 R__ASSERT(a.IsValid());
182
183 ResetStatus();
184 if (a.GetNrows() != a.GetNcols() || a.GetRowLwb() != a.GetColLwb()) {
185 Error("SetMatrix(const TMatrixDSym &","matrix should be square");
186 return;
187 }
188
190 fCondition = -1.0;
191
192 fRowLwb = a.GetRowLwb();
193 fColLwb = a.GetColLwb();
194 fU.ResizeTo(a);
195 fU = a;
196}
197
198////////////////////////////////////////////////////////////////////////////////
199/// Solve equations Ax=b assuming A has been factored by Cholesky. The factor U is
200/// assumed to be in upper triang of fU. fTol is used to determine if diagonal
201/// element is zero. The solution is returned in b.
202
204{
205 R__ASSERT(b.IsValid());
206 if (TestBit(kSingular)) {
207 Error("Solve()","Matrix is singular");
208 return kFALSE;
209 }
210 if ( !TestBit(kDecomposed) ) {
211 if (!Decompose()) {
212 Error("Solve()","Decomposition failed");
213 return kFALSE;
214 }
215 }
216
217 if (fU.GetNrows() != b.GetNrows() || fU.GetRowLwb() != b.GetLwb()) {
218 Error("Solve(TVectorD &","vector and matrix incompatible");
219 return kFALSE;
220 }
221
222 const Int_t n = fU.GetNrows();
223
224 const Double_t *pU = fU.GetMatrixArray();
225 Double_t *pb = b.GetMatrixArray();
226
227 Int_t i;
228 // step 1: Forward substitution on U^T
229 for (i = 0; i < n; i++) {
230 const Int_t off_i = i*n;
231 if (pU[off_i+i] < fTol)
232 {
233 Error("Solve(TVectorD &b)","u[%d,%d]=%.4e < %.4e",i,i,pU[off_i+i],fTol);
234 return kFALSE;
235 }
236 Double_t r = pb[i];
237 for (Int_t j = 0; j < i; j++) {
238 const Int_t off_j = j*n;
239 r -= pU[off_j+i]*pb[j];
240 }
241 pb[i] = r/pU[off_i+i];
242 }
243
244 // step 2: Backward substitution on U
245 for (i = n-1; i >= 0; i--) {
246 const Int_t off_i = i*n;
247 Double_t r = pb[i];
248 for (Int_t j = i+1; j < n; j++)
249 r -= pU[off_i+j]*pb[j];
250 pb[i] = r/pU[off_i+i];
251 }
252
253 return kTRUE;
254}
255
256////////////////////////////////////////////////////////////////////////////////
257/// Solve equations Ax=b assuming A has been factored by Cholesky. The factor U is
258/// assumed to be in upper triang of fU. fTol is used to determine if diagonal
259/// element is zero. The solution is returned in b.
260
262{
263 TMatrixDBase *b = const_cast<TMatrixDBase *>(cb.GetMatrix());
264 R__ASSERT(b->IsValid());
265 if (TestBit(kSingular)) {
266 Error("Solve()","Matrix is singular");
267 return kFALSE;
268 }
269 if ( !TestBit(kDecomposed) ) {
270 if (!Decompose()) {
271 Error("Solve()","Decomposition failed");
272 return kFALSE;
273 }
274 }
275
276 if (fU.GetNrows() != b->GetNrows() || fU.GetRowLwb() != b->GetRowLwb())
277 {
278 Error("Solve(TMatrixDColumn &cb","vector and matrix incompatible");
279 return kFALSE;
280 }
281
282 const Int_t n = fU.GetNrows();
283
284 const Double_t *pU = fU.GetMatrixArray();
285 Double_t *pcb = cb.GetPtr();
286 const Int_t inc = cb.GetInc();
287
288 Int_t i;
289 // step 1: Forward substitution U^T
290 for (i = 0; i < n; i++) {
291 const Int_t off_i = i*n;
292 const Int_t off_i2 = i*inc;
293 if (pU[off_i+i] < fTol)
294 {
295 Error("Solve(TMatrixDColumn &cb)","u[%d,%d]=%.4e < %.4e",i,i,pU[off_i+i],fTol);
296 return kFALSE;
297 }
298 Double_t r = pcb[off_i2];
299 for (Int_t j = 0; j < i; j++) {
300 const Int_t off_j = j*n;
301 r -= pU[off_j+i]*pcb[j*inc];
302 }
303 pcb[off_i2] = r/pU[off_i+i];
304 }
305
306 // step 2: Backward substitution U
307 for (i = n-1; i >= 0; i--) {
308 const Int_t off_i = i*n;
309 const Int_t off_i2 = i*inc;
310 Double_t r = pcb[off_i2];
311 for (Int_t j = i+1; j < n; j++)
312 r -= pU[off_i+j]*pcb[j*inc];
313 pcb[off_i2] = r/pU[off_i+i];
314 }
315
316 return kTRUE;
317}
318
319////////////////////////////////////////////////////////////////////////////////
320/// Matrix determinant det = d1*TMath::Power(2.,d2) is square of diagProd
321/// of cholesky factor
322
324{
325 if ( !TestBit(kDetermined) ) {
326 if ( !TestBit(kDecomposed) )
327 Decompose();
328 TDecompBase::Det(d1,d2);
329 // square det as calculated by above
330 fDet1 *= fDet1;
331 fDet2 += fDet2;
333 }
334 d1 = fDet1;
335 d2 = fDet2;
336}
337
338////////////////////////////////////////////////////////////////////////////////
339/// For a symmetric matrix A(m,m), its inverse A_inv(m,m) is returned .
340
342{
343 if (inv.GetNrows() != GetNrows() || inv.GetRowLwb() != GetRowLwb()) {
344 Error("Invert(TMatrixDSym &","Input matrix has wrong shape");
345 return kFALSE;
346 }
347
348 inv.UnitMatrix();
349
350 const Int_t colLwb = inv.GetColLwb();
351 const Int_t colUpb = inv.GetColUpb();
352 Bool_t status = kTRUE;
353 for (Int_t icol = colLwb; icol <= colUpb && status; icol++) {
354 TMatrixDColumn b(inv,icol);
355 status &= Solve(b);
356 }
357
358 return status;
359}
360
361////////////////////////////////////////////////////////////////////////////////
362/// For a symmetric matrix A(m,m), its inverse A_inv(m,m) is returned .
363
365{
366 const Int_t rowLwb = GetRowLwb();
367 const Int_t rowUpb = rowLwb+GetNrows()-1;
368
369 TMatrixDSym inv(rowLwb,rowUpb);
370 inv.UnitMatrix();
371 status = Invert(inv);
372
373 return inv;
374}
375
376////////////////////////////////////////////////////////////////////////////////
377/// Print class members .
378
380{
382 fU.Print("fU");
383}
384
385////////////////////////////////////////////////////////////////////////////////
386/// Assignment operator
387
389{
390 if (this != &source) {
392 fU.ResizeTo(source.fU);
393 fU = source.fU;
394 }
395 return *this;
396}
397
398////////////////////////////////////////////////////////////////////////////////
399/// Solve min {(A . x - b)^T (A . x - b)} for vector x where
400/// A : (m x n) matrix, m >= n
401/// b : (m) vector
402/// x : (n) vector
403
405{
407 Bool_t ok;
408 return ch.Solve(TMatrixD(TMatrixD::kTransposed,A)*b,ok);
409}
410
411////////////////////////////////////////////////////////////////////////////////
412/// Solve min {(A . x - b)^T W (A . x - b)} for vector x where
413/// A : (m x n) matrix, m >= n
414/// b : (m) vector
415/// x : (n) vector
416/// W : (m x m) weight matrix with W(i,j) = 1/std(i)^2 for i == j
417/// = 0 for i != j
418
419TVectorD NormalEqn(const TMatrixD &A,const TVectorD &b,const TVectorD &std)
420{
421 if (!AreCompatible(b,std)) {
422 ::Error("NormalEqn","vectors b and std are not compatible");
423 return TVectorD();
424 }
425
426 TMatrixD mAw = A;
427 TVectorD mBw = b;
428 for (Int_t irow = 0; irow < A.GetNrows(); irow++) {
429 TMatrixDRow(mAw,irow) *= 1/std(irow);
430 mBw(irow) /= std(irow);
431 }
433 Bool_t ok;
434 return ch.Solve(TMatrixD(TMatrixD::kTransposed,mAw)*mBw,ok);
435}
436
437////////////////////////////////////////////////////////////////////////////////
438/// Solve min {(A . X_j - B_j)^T (A . X_j - B_j)} for each column j in
439/// B and X
440/// A : (m x n ) matrix, m >= n
441/// B : (m x nb) matrix, nb >= 1
442/// mX : (n x nb) matrix
443
445{
448 ch.MultiSolve(mX);
449 return mX;
450}
451
452////////////////////////////////////////////////////////////////////////////////
453/// Solve min {(A . X_j - B_j)^T W (A . X_j - B_j)} for each column j in
454/// B and X
455/// ~~~
456/// A : (m x n ) matrix, m >= n
457/// B : (m x nb) matrix, nb >= 1
458/// mX : (n x nb) matrix
459/// W : (m x m) weight matrix with W(i,j) = 1/std(i)^2 for i == j
460/// = 0 for i != j
461/// ~~~
462
463TMatrixD NormalEqn(const TMatrixD &A,const TMatrixD &B,const TVectorD &std)
464{
465 TMatrixD mAw = A;
466 TMatrixD mBw = B;
467 for (Int_t irow = 0; irow < A.GetNrows(); irow++) {
468 TMatrixDRow(mAw,irow) *= 1/std(irow);
469 TMatrixDRow(mBw,irow) *= 1/std(irow);
470 }
471
474 ch.MultiSolve(mX);
475 return mX;
476}
bool Bool_t
Definition: RtypesCore.h:63
int Int_t
Definition: RtypesCore.h:45
const Bool_t kFALSE
Definition: RtypesCore.h:101
const Bool_t kTRUE
Definition: RtypesCore.h:100
const char Option_t
Definition: RtypesCore.h:66
TVectorD NormalEqn(const TMatrixD &A, const TVectorD &b)
Solve min {(A .
ClassImp(TDecompChol)
#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:187
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 Float_t Float_t b
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
Definition: TMatrixTBase.h:123
void Print(Option_t *name="") const override
Print the matrix as a table of elements.
Int_t GetRowLwb() const
Definition: TMatrixTBase.h:121
const TMatrixTBase< Element > * GetMatrix() const
Int_t GetInc() const
Element * GetPtr() const
const Element * GetMatrixArray() const override
Definition: TMatrixT.h:222
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...
Definition: TMatrixT.cxx:1211
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition: TObject.h:187
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:696
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:893
const Int_t n
Definition: legend1.C:16
static double B[]
static double A[]
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition: TMath.h:659
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
auto * a
Definition: textangle.C:12