Matrix inversion problem

ROOT
1

Hello ROOT talkers,

I am doing a little soft to perform analytical calculations of the track parameters resolution doing vertex detector stand-alone tracking. I am having a problem of matrix inversion. For the moment I have to invert matrices no bigger than 12x12 (6 layers with 2 measurements). I was using the TDecompLU ROOT class by doing as follows,

TDecompLU lu_meas(A);
lu.SetTol(1.0e-20);
int nr = 0;
while(!lu.Invert(InvA) && nr < 10) {
lu.SetMatrix(A);
lu.SetTol(0.1*lu_meas.GetTol());
nr++;
}

where A is the matrix to invert and InvA is the inverted matrix. Both, A and InvA are TMatrixD objects.

I noticed the problem by checking the product A*InvA which is sometimes significantly different from the unit matrix: diagonal (non-diagonal) elements significantly different from 1 (0).

I was wondering if anyone know something about matrix inversion either in ROOT and/or C++.

Thanks in advance,

Alejandro

2

Hi,

If your matrix is symmetric and positive defined, like maybe is the case if it represents a covariance matrix, you can try to use the Choleski decomposition (TDecompLU), otherwise you might try the SVD decomposition (TDecompSVD).
But what is the condition number ? If it is high it is expected that the inversion will suffer from numerical error.

See also root.cern.ch/root/htmldoc/guide … mpositions

Lorenzo

3

Hi Alejandro,

I assume that your are trying to find the solution x in A x = b by doing
x = A^-1 b .

This is both from the point of view of accuracy and calculation time not the optimal
way to do it. There is no reason to calculate the inverted matrix, better use the
back-substitution method as given in TDecompXX::Solve(…).

-Eddy