# TMatrixD Inversion Issue

Hi all-
I have a covariance matrix that I need to invert to use for a Chi-2 calculation (X^2 = M^(T) * C^(-1) * M).
I first convert the covariance matrix from a TH2 to a TMatrixD. Then I invert it and assigned the inverted matrix back to another TH2. I then attempt to multiply the two TH2’s (doing standard matrix multiplication) to get the identity. What I get is nonsense. The TH2’s I’m multiplying are `hCovMatrixTotal_alt` and `hCovMatrixTotalInv`. I keep in mind that TH2’s are enumerated from the bottom-up for rows and from left to right for columns.

``````        //Need to alter covariance matrix to exclude the first two bins that have zero values
//Essentially deletes the whole first "offset" number of columns and rows of the cov. matrix above
int offset(2);
TH2D *hCovMatrixTotal_alt = new TH2D("Total Cov. Matrix for Ratio "+variable1+" to "+variable2, "Total Cov. Matrix for Ratio "+variable1+" to "+variable2, effBinNum-offset, 0, effBinNum-offset, effBinNum-offset, 0, effBinNum-offset);
for (int i(1); i < effBinNum+1-offset; i++){
for (int j(1); j < effBinNum+1-offset; j++){
hCovMatrixTotal_alt->SetBinContent(i, j, hCovMatrixTotal->GetBinContent(i+offset, j+offset));
}
}

//Find inverse of covariance matrix; kInverted uses LU decomposition by default(?)
TMatrixD matCovMatrixTotal(effBinNum-offset, effBinNum-offset);
//Matrix bin numbering different from histogram
for (int i(0); i < effBinNum-offset; i++){
for (int j(0); j < effBinNum-offset; j++){
matCovMatrixTotal(i,j) = hCovMatrixTotal_alt->GetBinContent(i+1,j+1);
}
}

TMatrixD matCovMatrixTotalInv(TMatrixD::kInverted, matCovMatrixTotal);

//Make histogram out of inverted TMatrixD
TH2D *hCovMatrixTotalInv = new TH2D("Total Cov. Matrix Inv. for Ratio "+variable1+" to "+variable2, "Total Cov. Matrix Inv. for Ratio "+variable1+" to "+variable2, effBinNum-offset, 0, effBinNum-offset, effBinNum-offset, 0, effBinNum-offset);
for (int i(0); i < effBinNum-offset; i++){
for (int j(0); j < effBinNum-offset; j++){
hCovMatrixTotalInv->SetBinContent(i+1, j+1, matCovMatrixTotalInv(i,j));
}
}

//Check to see if these two TH2's multiplied together makes the identity!
//Display first column of resultant matrix
std::vector<double> column(effBinNum-offset);
for (int i(1); i < effBinNum+1-offset; i++){
for (int j(1); j < effBinNum+1-offset; j++){
column[i-1] += hCovMatrixTotal_alt->GetBinContent(-(i-1)+effBinNum-offset,j)*hCovMatrixTotalInv->GetBinContent(-(j-1)+effBinNum-offset,1);
}
cout << column[i-1] << endl;
}

``````

This gives me:

``````-27.2628
9.35068
-0.697222
3.34982
1.84428
1.87325
2.51832
2.71394
3.14706
3.18361
2.26186
``````

if I replace the line
`column[i-1] += hCovMatrixTotal_alt->GetBinContent(-(i-1)+effBinNum-offset,j)*hCovMatrixTotalInv->GetBinContent(-(j-1)+effBinNum-offset,1);`
with
`column[i-1] += hCovMatrixTotal_alt->GetBinContent(i,j)*hCovMatrixTotalInv->GetBinContent(j,1);`
I essentially get the identity matrix up-side down (the first element printed is the zero element):

``````1
2.34188e-16
-1.09288e-16
4.97866e-16
2.37657e-16
-2.94903e-17
1.73472e-17
6.64399e-16
5.20417e-18
4.19803e-16
-1.62197e-16
8.04912e-16
``````

This second kind of iteration option would yield the identity, if both matrices were flipped upside-down.

Does using the kInverted variable flip the rows and columns assignment of the inverted matrix in any way? Or could someone offer advice on how to correctly do this calculation?

Best,
Andrew

Hi Andrew,

Can you please post a running script, so we can reproduce your issue. If you are using covariance matrix, which are symmetric and positive defined, I would recommend you using the TMatrixDSym class and the Cholesky decomposition, TDecompChol.
Otherwise, inversion can be difficult numerically, especially form matrix with high condition number.
What is the value returned from TDecompBase::GetCondition ?

Cheers

Lorenzo

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.