Hi Niladri,
Assuming that you want to invert the covariance as you seem to point out
in another forum entry
[url]Error in <TDecompChol::Decompose()>: matrix not positive definite
the idea is the following. Suppose you have positive definite matrix M,
define a new matrix M_sc, where
M_sc[i][j] = M[i][j]/sqrt(M[i][i]*M[j][j])
The matrix M_sc has now 1 on every diagonal position and off-diagonal <1and > -1,
so numerically a well-behaved matrix
Now invert M_sc => M_sc^-1 and repeat the scaling:
(M_sc^-1)_sc = M_sc^-1[i][j]/sqrt(M[i][i]*M[j][j])
Now the claim is M^-1 = (M_sc^-1)_sc
Proof:
M is positive definite => M = L^T L
define a diagonal matrix S with entry S[i][i] = 1/sqrt(M[i][i])
M_sc = S L^T L S
M_sc^-1 = ( S L^T L S)^-1 = S^-1 L^-1 (L^-1)^T (S^-1)
S M_sc^-1 S = S (S^-1 L^-1 (L^-1)^T (S^-1)) S = L^-1 (L^-1)^T = (L^T L)^-1 = M^-1
Implementation in ROOT:
const TVectorD diag = TMatrixDDiag(m);
const TMatrixDSym m_sc = m;
m.NormByDiag(diag,“D”);
TDecompChol chol(m_sc);
TMatrixDSym m_sc_inv = chol.Invert();
m_inv = m_sc_inv;
m_inv.NormByDiag(diag,“M”);