How to Invert a Matrix ?

The most straightforward way to invert a matrix is to invoke the Invert(Double_t &det=0) function or the appropriate constructor: 

   TMatrixD a(...);
   a.Invert();
or: 
   TMatrixD b(kInvert,a);
This functionality is available for both general and symmetric matrices. There is a slight variant of the Invert function available, InvertFast(Double_t &det=0). For matrix sizes <= (6x6), the Cramer algorithm will be applied which is faster but less accurate. The inversion can also be achieved in a more controlled way by using directly the decomposition classes:

 

Name Matrix Type Comment
 TDecompLU  general
 
 TDecompQRH  general  
 TDecompSVD  general  Can manipulate singular matrix.
 TDecompBK  symmetric  
 TDecompChol  symmetric  Matrix should also be positive definite.
 TDecompSparse  sparse  

If the required matrix type is general, it can of course also handle symmetric matrices.

The following example shows how to check whether the matrix is singular before attempting to invert it. 

   TDecompLU lu(a);
   TMatrixD b;
   if (!lu.Decompose()) {
      cout << "Decomposition failed, matrix singular ?" << endl;
      cout << "condition number = " << = a.GetCondition() << endl;
   } else {
      lu.Invert(b);
   }
For an extensive discussion, see the tutorial invertMatrix.C example.

  • Printer-friendly version
  • PDF version