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