Logo ROOT   6.07/09
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 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
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 
379 void TDecompChol::Print(Option_t *opt) const
380 {
381  TDecompBase::Print(opt);
382  fU.Print("fU");
383 }
384 
385 ////////////////////////////////////////////////////////////////////////////////
386 /// Assignment operator
387 
389 {
390  if (this != &source) {
391  TDecompBase::operator=(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 
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 
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 }
Int_t GetRowLwb() const
Definition: TMatrixTBase.h:132
virtual TMatrixTBase< Element > & UnitMatrix()
Make a unit matrix (matrix need not be a square one).
virtual void Det(Double_t &d1, Double_t &d2)
Matrix determinant det = d1*TMath::Power(2.,d2) is square of diagProd of cholesky factor...
static double B[]
virtual Bool_t MultiSolve(TMatrixD &B)
Solve set of equations with RHS in columns of B.
Element * GetPtr() const
Double_t fCondition
Definition: TDecompBase.h:42
Bool_t IsValid() const
Definition: TVectorT.h:89
Int_t fRowLwb
Definition: TDecompBase.h:43
const char Option_t
Definition: RtypesCore.h:62
Int_t GetNrows() const
Definition: TVectorT.h:81
void Print(Option_t *opt="") const
Print class members.
Decomposition Base class.
Definition: TDecompBase.h:36
#define R__ASSERT(e)
Definition: TError.h:98
double inv(double x)
For comparisons.
Definition: inv.h:58
TMatrixD fU
Definition: TDecompChol.h:32
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TArc * a
Definition: textangle.C:12
const Bool_t kFALSE
Definition: Rtypes.h:92
STL namespace.
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1)
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
Definition: TMatrixT.cxx:1210
static double A[]
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:739
const TMatrixDSym GetMatrix()
Reconstruct the original matrix using the decomposition parts.
Double_t fTol
Definition: TDecompBase.h:39
Element Norm1() const
Definition: TMatrixTBase.h:188
Bool_t IsValid() const
Definition: TMatrixTBase.h:157
Element GetTol() const
Definition: TMatrixTBase.h:139
TMatrixDSym Invert()
Definition: TDecompChol.h:64
TDecompBase & operator=(const TDecompBase &source)
Assignment operator.
virtual void SetMatrix(const TMatrixDSym &a)
Set the matrix to be decomposed, decomposition status is reset.
Int_t fColLwb
Definition: TDecompBase.h:44
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 .
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:925
TVectorT< Double_t > TVectorD
Definition: TVectorDfwd.h:24
Int_t GetColUpb() const
Definition: TMatrixTBase.h:136
TMatrixTRow< Double_t > TMatrixDRow
Int_t GetNrows() const
Definition: TMatrixTBase.h:134
Element * GetMatrixArray()
Definition: TVectorT.h:84
TMatrixT< Double_t > TMatrixD
Definition: TMatrixDfwd.h:24
TVectorD NormalEqn(const TMatrixD &A, const TVectorD &b)
Solve min {(A .
void Print(Option_t *name="") const
Print the matrix as a table of elements.
Cholesky Decomposition class.
Definition: TDecompChol.h:28
const double tol
TRandom2 r(17)
virtual Bool_t Solve(TVectorD &b)
Solve equations Ax=b assuming A has been factored by Cholesky.
Bool_t TestBit(UInt_t f) const
Definition: TObject.h:159
virtual const Element * GetMatrixArray() const
Definition: TMatrixT.h:223
void ResetStatus()
Definition: TDecompBase.h:46
const TMatrixTBase< Element > * GetMatrix() const
Int_t GetColLwb() const
Definition: TMatrixTBase.h:135
Linear Algebra Package.
Double_t fDet2
Definition: TDecompBase.h:41
Double_t fDet1
Definition: TDecompBase.h:40
Int_t GetRowLwb() const
Definition: TDecompBase.h:76
virtual Int_t GetNrows() const
Definition: TDecompChol.h:47
#define ClassImp(name)
Definition: Rtypes.h:279
Int_t GetLwb() const
Definition: TVectorT.h:79
double Double_t
Definition: RtypesCore.h:55
TMatrixTSym< Double_t > TMatrixDSym
Int_t GetNcols() const
Definition: TMatrixTBase.h:137
virtual Bool_t Decompose()
Matrix A is decomposed in component U so that A = U^T * U If the decomposition succeeds, bit kDecomposed is set , otherwise kSingular.
virtual void Det(Double_t &d1, Double_t &d2)
Matrix determinant det = d1*TMath::Power(2.,d2)
TDecompChol & operator=(const TDecompChol &source)
Assignment operator.
Int_t GetInc() const
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
const Bool_t kTRUE
Definition: Rtypes.h:91
void Print(Option_t *opt="") const
Print class members .
const Int_t n
Definition: legend1.C:16