Logo ROOT   6.21/01
Reference Guide
invertMatrix.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_matrix
3 /// \notebook -nodraw
4 /// This macro shows several ways to invert a matrix . Each method
5 /// is a trade-off between accuracy of the inversion and speed.
6 /// Which method to chose depends on "how well-behaved" the matrix is.
7 /// This is best checked through a call to Condition(), available in each
8 /// decomposition class. A second possibility (less preferred) would be to
9 /// check the determinant
10 ///
11 /// #### USAGE
12 ///
13 /// This macro can be executed with Cling or ACLIC
14 /// - via the interpretor, do
15 /// ~~~{.cpp}
16 /// root > .x invertMatrix.C
17 /// ~~~
18 /// - via ACLIC
19 /// ~~~{.cpp}
20 /// root > gSystem->Load("libMatrix");
21 /// root > .x invertMatrix.C+
22 /// ~~~
23 ///
24 /// \macro_output
25 /// \macro_code
26 ///
27 /// \author Eddy Offermann
28 
29 #include <iostream>
30 #include "TMath.h"
31 #include "TMatrixD.h"
32 #include "TMatrixDLazy.h"
33 #include "TVectorD.h"
34 #include "TDecompLU.h"
35 #include "TDecompSVD.h"
36 
37 void invertMatrix(Int_t msize=6)
38 {
39  if (msize < 2 || msize > 10) {
40  std::cout << "2 <= msize <= 10" <<std::endl;
41  return;
42  }
43  std::cout << "--------------------------------------------------------" <<std::endl;
44  std::cout << "Inversion results for a ("<<msize<<","<<msize<<") matrix" <<std::endl;
45  std::cout << "For each inversion procedure we check the maximum size " <<std::endl;
46  std::cout << "of the off-diagonal elements of Inv(A) * A " <<std::endl;
47  std::cout << "--------------------------------------------------------" <<std::endl;
48 
49  TMatrixD H_square = THilbertMatrixD(msize,msize);
50 
51  // ### 1. InvertFast(Double_t *det=0)
52  // It is identical to Invert() for sizes > 6 x 6 but for smaller sizes, the
53  // inversion is performed according to Cramer's rule by explicitly calculating
54  // all Jacobi's sub-determinants . For instance for a 6 x 6 matrix this means:
55  // \# of 5 x 5 determinant : 36
56  // \# of 4 x 4 determinant : 75
57  // \# of 3 x 3 determinant : 80
58  // \# of 2 x 2 determinant : 45 (see TMatrixD/FCramerInv.cxx)
59  //
60  // The only "quality" control in this process is to check whether the 6 x 6
61  // determinant is unequal 0 . But speed gains are significant compared to Invert() ,
62  // up to an order of magnitude for sizes <= 4 x 4
63  //
64  // The inversion is done "in place", so the original matrix will be overwritten
65  // If a pointer to a Double_t is supplied the determinant is calculated
66  //
67 
68  std::cout << "1. Use .InvertFast(&det)" <<std::endl;
69  if (msize > 6)
70  std::cout << " for ("<<msize<<","<<msize<<") this is identical to .Invert(&det)" <<std::endl;
71 
72  Double_t det1;
73  TMatrixD H1 = H_square;
74  H1.InvertFast(&det1);
75 
76  // Get the maximum off-diagonal matrix value . One way to do this is to set the
77  // diagonal to zero .
78 
79  TMatrixD U1(H1,TMatrixD::kMult,H_square);
80  TMatrixDDiag diag1(U1); diag1 = 0.0;
81  const Double_t U1_max_offdiag = (U1.Abs()).Max();
82  std::cout << " Maximum off-diagonal = " << U1_max_offdiag << std::endl;
83  std::cout << " Determinant = " << det1 << std::endl;
84 
85  // ### 2. Invert(Double_t *det=0)
86  // Again the inversion is performed in place .
87  // It consists out of a sequence of calls to the decomposition classes . For instance
88  // for the general dense matrix TMatrixD the LU decomposition is invoked:
89  // - The matrix is decomposed using a scheme according to Crout which involves
90  // "implicit partial pivoting", see for instance Num. Recip. (we have also available
91  // a decomposition scheme that does not the scaling and is therefore even slightly
92  // faster but less stable)
93  // With each decomposition, a tolerance has to be specified . If this tolerance
94  // requirement is not met, the matrix is regarded as being singular. The value
95  // passed to this decomposition, is the data member fTol of the matrix . Its
96  // default value is DBL_EPSILON, which is defined as the smallest number so that
97  // 1+DBL_EPSILON > 1
98  // - The last step is a standard forward/backward substitution .
99  //
100  // It is important to realize that both InvertFast() and Invert() are "one-shot" deals , speed
101  // comes at a price . If something goes wrong because the matrix is (near) singular, you have
102  // overwritten your original matrix and no factorization is available anymore to get more
103  // information like condition number or change the tolerance number .
104  //
105  // All other calls in the matrix classes involving inversion like the ones with the "smart"
106  // constructors (kInverted,kInvMult...) use this inversion method .
107  //
108 
109  std::cout << "2. Use .Invert(&det)" << std::endl;
110 
111  Double_t det2;
112  TMatrixD H2 = H_square;
113  H2.Invert(&det2);
114 
115  TMatrixD U2(H2,TMatrixD::kMult,H_square);
116  TMatrixDDiag diag2(U2); diag2 = 0.0;
117  const Double_t U2_max_offdiag = (U2.Abs()).Max();
118  std::cout << " Maximum off-diagonal = " << U2_max_offdiag << std::endl;
119  std::cout << " Determinant = " << det2 << std::endl;
120 
121  // ### 3. Inversion through LU decomposition
122  // The (default) algorithms used are similar to 2. (Not identical because in 2, the whole
123  // calculation is done "in-place". Here the original matrix is copied (so more memory
124  // management => slower) and several operations can be performed without having to repeat
125  // the decomposition step .
126  // Inverting a matrix is nothing else than solving a set of equations where the rhs is given
127  // by the unit matrix, so the steps to take are identical to those solving a linear equation :
128  //
129 
130  std::cout << "3. Use TDecompLU" << std::endl;
131 
132  TMatrixD H3 = H_square;
133  TDecompLU lu(H_square);
134 
135  // Any operation that requires a decomposition will trigger it . The class keeps
136  // an internal state so that following operations will not perform the decomposition again
137  // unless the matrix is changed through SetMatrix(..)
138  // One might want to proceed more cautiously by invoking first Decompose() and check its
139  // return value before proceeding....
140 
141  lu.Invert(H3);
142  Double_t d1_lu; Double_t d2_lu;
143  lu.Det(d1_lu,d2_lu);
144  Double_t det3 = d1_lu*TMath::Power(2.,d2_lu);
145 
146  TMatrixD U3(H3,TMatrixD::kMult,H_square);
147  TMatrixDDiag diag3(U3); diag3 = 0.0;
148  const Double_t U3_max_offdiag = (U3.Abs()).Max();
149  std::cout << " Maximum off-diagonal = " << U3_max_offdiag << std::endl;
150  std::cout << " Determinant = " << det3 << std::endl;
151 
152  // ### 4. Inversion through SVD decomposition
153  // For SVD and QRH, the (n x m) matrix does only have to fulfill n >=m . In case n > m
154  // a pseudo-inverse is calculated
155  std::cout << "4. Use TDecompSVD on non-square matrix" << std::endl;
156 
157  TMatrixD H_nsquare = THilbertMatrixD(msize,msize-1);
158 
159  TDecompSVD svd(H_nsquare);
160 
161  TMatrixD H4 = svd.Invert();
162  Double_t d1_svd; Double_t d2_svd;
163  svd.Det(d1_svd,d2_svd);
164  Double_t det4 = d1_svd*TMath::Power(2.,d2_svd);
165 
166  TMatrixD U4(H4,TMatrixD::kMult,H_nsquare);
167  TMatrixDDiag diag4(U4); diag4 = 0.0;
168  const Double_t U4_max_offdiag = (U4.Abs()).Max();
169  std::cout << " Maximum off-diagonal = " << U4_max_offdiag << std::endl;
170  std::cout << " Determinant = " << det4 << std::endl;
171 }
int Int_t
Definition: RtypesCore.h:41
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:725
TMatrixT.
Definition: TMatrixDfwd.h:22
LU Decomposition class.
Definition: TDecompLU.h:23
Single Value Decomposition class.
Definition: TDecompSVD.h:23
TMatrixT< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant.
Definition: TMatrixT.cxx:1396
THilbertMatrixT< Double_t > THilbertMatrixD
Definition: TMatrixDLazy.h:35
double Double_t
Definition: RtypesCore.h:55
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
TMatrixT< Element > & InvertFast(Double_t *det=0)
Invert the matrix and calculate its determinant, however upto (6x6) a fast Cramer inversion is used ...
Definition: TMatrixT.cxx:1410