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 (!a.Decompose()) {
         cout << "Decomposition failed, matrix singular ?" << endl;
         cout << "condition number = " << = a.GetCondition() << endl;
      } else {
         a.Invert(b);
      }

For an extensive discussion, see the tutorial invertMatrix.C example.

Eddy Offermann
Last update: 25-Nov-2004 by Ilka Antcheva