TMatrixD not correctly inverted due to rounding errors

ROOT
1

As the title suggests, I have some problems with a TMatrixD that does not behave completely as I like it. It needs to be inverted and as a test I can simply multiply the inverse with the original. Please find the matrix, its inverse and the product below.

How can I calculate the inverse with full efficiency?

Any help you can provide is most welcome!

Cheers,
Alex

TMatrixD mInvertCovariance = (*mCovariance);
mInvertCovariance.Invert();

mCovariance->Print();
mInvertCovariance.Print();

mInvertCovariance *= (*mCovariance);
mInvertCovariance.Print();
2

Hi,
Have you tried other decompositions ?
See root.cern.ch/how/how-invert-matrix

Also what is the condition number of the matrix ?

Lorenzo

3

Check also the User’s Guide for explanation of condition number:

root.cern.ch/root/htmldoc/guide … ion-number

and also have a look at

tutorials/matrix/invertMatrix.C

4

Since your matrix carries “Covariance” in its name, I assume that it is symmetric and positive-definite. In that case you should be using a Cholesky decomposition (TDecompChol is ROOT’s implementation) to calculate the inverse.

In an ideal world

  1. one formulates the problem such that one never inverts matrices. Instead, one uses the decompositions to solve linear 1D problems, and
  2. one doesn’t use the covariance matrices for calculations but their square roots (these are the matrices corresponding to the triangular matrix H in the definition of the Cholesky decomposition, C=HH^t), as the square roots are better conditioned in general. Most formulas that typically appear can be decomposed into formulas on square roots

E.g., in the matrix generalization of the Pythagorean theorem which goes like
C3 = C1+C2 with H1, H2, H3 the corresponding square roots, i.e. this equation is equivalent to
H3 H3^t=H1 H1^t + H2 H2^t, which we can rewrite as (sorry for the awkward notation with the double transpose, but I don’t see a LaTeX input option which would allow me to write columns; the objects in parentheses are block matrices; th dot denotes matrix multiplication)
H1 H1^t + H2 H2^t = (H1^t,H2^t) . (H1^t,H2^t)^t
Now (H1^t,H2^t)^t can be brought into upper triangle form (a triangle of the same dimensions as H1, followed by an equal number of zero rows) by an orthogonal transformation O, O^t = O^-1. Call the triangle H3. That is (I’m wrting the transposed because I cannot write columns)
(H1^t,H2^t) = (H3^t, 0).O^t,
and thus
H1 H1^t + H2 H2^t = (H3^t, 0).O^t.O.(H3^t,0)^t = H3^t.H3,
and we have obtained the square root of C3 using square roots only (notice how the dimensionality magically reduces in the second equality, first the orthogonal transformation pops away, then the zero entries do the rest). In practice the orthogonal transformation needs not be calculated explicitly, one can calculate H3 directly using a sequence of householder transformations. Finally, if one needs inverses somewhere, inverting triangular matrices is dirt cheap.

Another example would be a similarity transformation (say, error propagation with a Jacobian): instead of J.C.J^t one calculates J.H (or, if need be, the triangular J.H.O with O an orthogonal transformation).

Sorry for indulging, but I found this a very satisfying extension of my knowledge of linear algebra 8)

5

Hi Tobi,

Thanks for that nice explanation.

Each TDecomp{Chol,LU,QR,SVD,…) has a “Solve” method that calculates the solution
of a linear equation without inverting the matrix !

Actually, the Invert() method of the decomposition class uses Solve with a unit matrix.

-Eddy

6

[quote=“Eddy Offermann”]Each TDecomp{Chol,LU,QR,SVD,…) has a “Solve” method that calculates the solution
of a linear equation without inverting the matrix !

Actually, the Invert() method of the decomposition class uses Solve with a unit matrix.[/quote]Hi Eddy,

absolutely, it’s good that you point this out. Just to clarify (to our audience :smiley: ), what I was trying to say was that actually calculating the inverse is wasteful in most all circumstances, because one usually is interested in a vector, and not the inverse matrix itself. E.g. in a typical linear problem
A x = b
the solution is (existence of the inverse provided)
x = A^-1 b

But of course, one never needs to calculate A^-1. Once the matrix is decomposed one simply solves one time, instead of actually evaluating the second equation. Doing that would involve solving once for each column of the (inverse) matrix (i.e. “[one] uses Solve with a unit matrix”), followed by the matrix multiplication. Of course this is more expensive and numerically less stable because of the superfluous intermediate calculation. Similar considerations extend to least-squares problems.