Matrix inversion

ROOT
1

Hello,

I have a problem with matrix inversion:

Matrix ist small (< 10x10) and quadratic.
I invert the matrix using TMatrix (which is, I think: typedef TMatrixT TMatrix) but it happens with a lot of other root-matrices
I multiply the inverted matrix to the original matrix and look at the result, that should be unity by definition.
It isn’t. The diagonal elements are any number and the off-diagonal elements are far off the Zero (let’s say around 1e-1 to 1e+20).

Now the usual “send me a simple running example” doesn’t work, since this error doesn’t occur in “small simple running programs”. It only occurs if there is a bit of an overhead. This leads to the assumption that it is a memory problem. But the memory consumed is not so much that it could be an out of memory problem (<100Mb of 4Gb) and the prog is not so large, that I lost track of allocated memory -I checked that more than once. And it should not happen with the first inversion but, let’s say with the hundredth.

Can you give me a reason why this happens and how this problem can be solved.

K.Ehrhardt

PS: I discussed that with some colleagues and they gave me a lot of advice:
-don’t use THilbertMatrix of any type, there inversion always fails,
-don’t use TMatrixT, there inversion fails,
-compile your program, don’t let it run in CINT
But these didn’t change anything (since I already used them).

2

Hi,

Let me first react to the advice you have heard:

Always perform your inversion in double precision ! It is has by definition
a higher precision than a float calculation. However, the accuracy of the
answer might still not be sufficient, see below my explanation about the
condition number of a matrix. What your colleagues might have meant is
that storing in a file a double is not meaningful because of the limited
accuracy of your numbers . For that case ROOT has the type

typedef double Double32_t;  //Double 8 bytes in memory, written as a 4 bytes float

I do not know what this statement means, it is something like “do not
invert a singular matrix, it always fails”.
The Hilbert matrix has a very large condition number as its size grows
and therefore its accuracy decreases, again see condtion number.

CINT runs the same compiled matrix routines, so what would the difference be ?

Coming to your problem …
Please have a look at the user’s guide chapter 14, pages 266-267:
ftp://root.cern.ch/root/doc/14LinearAlgebra.pdf

Here the concept of matrix condition number is explained and its
relation to the accuracy of an inverted matrix . I guess that this settles
the case whether to double or float in an inversion and the Hilbert
matrix case.

If I understand your problem correctly, your inversion problems occur
in your larger program . My advice would be to not use the Invert()
functionality but to go step-by-step with a TDecomp class. An extensive
example is given in the tutorial: tutorials/matrix/invertMatrix.C
Now you can calculate the condition number and stream your problematic
cases to a TFile and analyze these cases outside your program.

Eddy

3

Hello,

[quote]Always perform your inversion in double precision ! It is has by definition
a higher precision than a float calculation.[/quote]
Yes, I thought that too. But what does higher precision help, if it doesn’t work?

[quote]What your colleagues might have meant is
that storing in a file a double is not meaningful because of the limited
accuracy of your numbers [/quote]
Nope, we were just talking about inversion- the matrix is not written to file.

You haven’t used root to an extend that I have if you still ask that question. In principle your statement should be correct. In principle. But the libraries have not heard yet about that rule. And that doesn’t only count for matrices but a number of other routines and data structures.

[quote]If I understand your problem correctly, your inversion problems occur in your larger program . My advice would be to not use the Invert()
functionality but to go step-by-step with a TDecomp class. [/quote]
Ok, the complete inversion is done in the TDecomp-classes? Good, I’ll check that, but I won’t be surprised if the problem remains.

I have to stress, that the matrix was all right. I wrote the matrix to file and reread it with another program inverting it successfully with the same, the very same procedure that failed before.
BUT I’ve got approximately 10^10 matrix inversions to do, each counting on the result of the calculation before (don’t fix me on an order of magnitude more or less). Writing them all to file and solving them outside in a separate program is out of question.

K.Ehrhardt

4

I am afraid that part of my message did not come across .

Precision and accuracy are two different things . You can do an inaccurate
calculation with very high precision .

The matrix condition number can give you an estimate of the accuracy
of the matrix inversion given a certain precision.

I was suggesting that you dump those matrices into a file
for which the A*A^-1 has large off-diagonal components and do
diagnostics on those cases.

To analyze your problem further we will need a short running example !

The problem could be as simple as that you forgot that A.Invert()
does not only return the A^-1 but actually inverts A itself etc etc.

Eddy

5

[quote]I was suggesting that you dump those matrices into a file
for which the A*A^-1 has large off-diagonal components and do
diagnostics on those cases.[/quote]
I’m afraid, I misunderstood your suggestion. But it s still not applicable since I cannot stop the program, when some inversion failed, read the matrix in another program and start the original program at that point again. And that 10^10 times.
But, I did that already, and wrote that in my first post to that topic. It is the very same matrix, that is inverted in the small program and is not inverted in the large program.

hm. The problem is, that it works in a short running example… But it breaks, if there is some overhead.

[quote]The problem could be as simple as that you forgot that A.Invert()
does not only return the A^-1 but actually inverts A itself etc etc. [/quote]

I see, you want some code. Some time ago I had a look at the root-linear-algebra-classes. And I found out, the hard way, that some basic features were not so easy to use, like transposing the matrix (protected instead of public) and other things. That’s why I defined and implemented an own linear-algebra package, along with some other geometry stuff. The classes there are derived from TObject to have the streamer()-routine. Everything in these classes is calculated using these geometry/lin-algebra classes and the basic c-types.
The one exception is (or better was) the inversion, that went like that:

matrixNxM matrixNxM::operator-()const { Double_t tmp[numberOfLines*numberOfColumns]; for(int i=0;i<numberOfLines;i++) { for(int j=0;j<numberOfColumns;j++) tmp[i*numberOfColumns+j]=getValue(i,j); } TMatrixD mat(numberOfLines,numberOfColumns,tmp); mat.Invert(); //mat.InvertFast(); if(!mat.IsValid()) return matrixNxM(0,0); matrixNxM ret(numberOfColumns,numberOfLines); for(int i=0;i<numberOfLines;i++) { for(int j=0;j<numberOfColumns;j++) ret.setValue(i,j,mat[i][j]); } return ret; }
This class and the inversion routine defined by the overwritten unary operator-() is used in a various number of fitting routines, like line-fitting etc. What gave the problems was a kinfit, where there are 8 matrices defined, one having dimension of 20x20, the other ones much smaller. The problematic matrix was 6x6.
The test program used the same matrix class, the same precision, the same numbers, the same code, compiled at the same hour, running almost in parallel.
This was tried with several compilers (yes, every code linked was compiled with the same compiler) from gcc2.95 to gcc 4.1.2 and several root5 subversions both on 32 and 64 bit machines with different memory sizes. The result was surprisingly constant.

6

You are performing matrix operations with your own package and
ran into problems , you do not supply code/example . So how can the ROOT community help you :confused: .

We are not aware of any matrix functionality that is protected instead of
public .

Eddy