I’m encountering some covariance matrices with small determinants that I need to invert. Root either gives me the “Singular matrix” error or an inverse matrix that is incorrect. (Mathematica can find a good approximate inverse.)
Here’s an example of my problem in Python. As you will see the determinant is of the order 10^-8, small but not zero.
Thanks for any ideas you have to fix or work around this problem.
I think that version worked because of the truncated values in the matrix. This version is more along the lines of what my C++ application is doing when it runs into problems. This version requires the attacted .root file.
[code]import os, sys, math, ROOT
def GetCovarianceMatrix( data ):
means = []
n = data.GetNrows()
d = data.GetNcols()
for i in range( n ):
for j in range ( d ):
try:
means[j] += data[i][j]/n
except IndexError:
means.append(data[i][j]/n)
covMat = ROOT.TMatrixD( d, d )
for i in range(n):
for j in range(d):
for k in range(d):
covMat[j][k] += (means[j] - data[i][j])*(means[k] - data[i][k])/(n - 1.0)
return covMat
def GetTMatrixFromTree( tree ):
vars = [“InvMass”, “mMet”, “m1metDphi”, “m1m2Dphi” , “MuCentrali”] #vars = [“InvMass”, “m1m2Dphi”, “mMet”]
mat = ROOT.TMatrixD(tree.GetEntries(), len(vars) )
print “Resized matrix . . .”
for i in range( 0, tree.GetEntries() ):
tree.GetEntry(i)
for j in range (len(vars)):
mat[i ][j] = tree.GetLeaf( vars[j] ).GetValue()
Using your tree.root file, I have created a script imatrix.C
reproducing Eddy’s style. I have run this script under linux and
windows with identical results (within machine precision).
I also ran imatrix.py (your script) on the two systems and I
obtain substantially different results.
I do not know if this is due to the different algorithms
in your respective code (no time to check).
See the results in the small tar file in attachement.
I’ve modified the C++ code to calculate the covariance matrix. It then attempts to invert the covariance matrix and gets the same results as the python script.
I ran your imatrix code (C++ version of your last submission)
with your root file and saw that the
covariance file coming from your GetCovarianceMatrix is :
5x5 matrix covMatrix is as follows
| 0 | 1 | 2 | 3 | 4 |
------------------------------------------------------------------
0 | 1.699 1.582 0.115 nan -0.1538
1 | 1.582 6.718 -0.4563 nan 1.501
2 | 0.115 -0.4563 0.5547 nan -0.2539
3 | nan nan nan nan nan
4 | -0.1538 1.501 -0.2539 nan 0.6157
Inverting this matrix will of course not work .
The reason is that you forget to initialize the array “means” to zero !
You also forget to do this with the covMat matrix .
OK - the C++ version I hastily translated from python this morning has problems. Can you initialize the values and try again or run this attached version? (I created the original python script in order to isolate the problem and make it easier for others to reproduce. Sorry this backfired.)
It’s not my real problem though. The real problem is that ROOT is not properly inverting the matrix covMatrix. (You already showed that it can invert the low precision version of covMatrix when the values in the matrix are converted strings rather than calculated double values. )
I’ve actually tried to isolate a problem which is occurring in data analysis. Some of the barely contributing backgrounds end up with small data samples and covariance matrices that won’t invert. This is one example.
Here are my results when I run both the C and the Python versions of the code.
So is this what one could expect ? Yes, it is .
The relevant quantity to see what the precision of an inversion will
be is the Condtion number (see User’s Guide) . (It is not the
determinant . Just take for example a unit matrix and scale the
diagonal with an arbitrary small number … )
The condition number of your matrix is : 1.56019e+17
The rule of thumb is that for double precision, a condition number
of 10^n will result in an accuracy of 15-n digits in the inverted matrix .
Mathematica is probably using some iterative scheme to calculate
the inversion . Our current package has not this option ,however, the
building blocks to set this up are there …