Inverting TMatrixD matricies with small elements

ROOT
1

Hello,

I am having trouble inverting TMatrixD matricies with small elements.  For an example:

cplager@PointyIII> root
  *******************************************
  *                                         *
  *        W E L C O M E  to  R O O T       *
  *                                         *
  *   Version   4.02/00  17 December 2004   *
  *                                         *
  *  You are welcome to visit our Web site  *
  *          http://root.cern.ch            *
  *                                         *
  *******************************************

FreeType Engine v2.1.3 used to render TrueType fonts.
Compiled for win32gcc with thread support.

CINT/ROOT C/C++ Interpreter version 5.15.159, Nov 14 2004
Type ? for help. Commands must be C++ statements.
Enclose multiple statements between { }.
root [0] TMatrixD m(1,1)
root [1] m(0,0) = 5.58332e-22;
root [2] m.Invert()
(class TMatrixD)284111368
root [3] m(0,0)
Fatal in <TMatrixD>: IsValid() violated at line 204 of `include/TMatrixD.h'
aborting

On a linux machine with an older version of root, I get slightly more information, but the same result:

cplager@b0pcucla10> root
  *******************************************
  *                                         *
  *        W E L C O M E  to  R O O T       *
  *                                         *
  *   Version   4.00/08   1 December 2004   *
  *                                         *
  *  You are welcome to visit our Web site  *
  *          http://root.cern.ch            *
  *                                         *
  *******************************************

FreeType Engine v2.1.3 used to render TrueType fonts.
Compiled for linux with thread support.

CINT/ROOT C/C++ Interpreter version 5.15.138, May 23 2004
Type ? for help. Commands must be C++ statements.
Enclose multiple statements between { }.
root [0] TMatrixD m(1,1)
root [1] m(0,0) = 5.58332e-22;
Inroot [2] m.Invert()
mError in <InvertLU>: LU[0,0]=5.5833e-22 < 2.2204e-16
((class TMatrixD)149216040
root [3] m(0,0)
Fatal in <TMatrixD>: IsValid() violated at line 179 of `include/TMatrixD.h'
aborting

Note that we can’t upgrade our version of root as it is tied to the rest of our production.

Is there a way around this? (Is it fixed in newer versions of root?)

Cheers,
Charles

2

Hi Charles,

Your problem is caused by the fact that the algoritm believes that your
matrix is singular . As a consequence, the matrix is made invalid and
no further operation is allowed on it, like operator [] , otherwise resulting
in an assert .

Generally, it is a good idead to check that the inversion suceeded through
a .IsValid() check .

The decision whether a matrix is singular is made by checking some
matrix values (in case of a standard matrix inversion a LU decomposition
is used and the important matrix values are the diagonals) against the
matrix tolerance which has a default value of epsilon=2.2e-16 .

In your case the matrix is just scaled by some very small value, so the
matrix is not singular . To prevent the algorithm from giving false
alarm , you can do 2 things :

  1. Change the matrix tolerance through a SetTol call , so in your
    case m.SetTol(1.e-23)
  2. scale the matrix before calling Invert() and scale back.

Indeed, the error message disappeared in version 4.02/00 but is
fortunately back again with a clearer message :

::Error(“TDecompLU::InvertLU”,“matrix is singular, %d diag elements < tolerance of %.4e”,nrZeros,tol);

Eddy