[Math] Numerical trace of inverse matrix from Cholesky

matricesmatrix inversenumerical linear algebratraces

This question was somewhat answered here: Fast trace of inverse of a square matrix. However, I feel like there was no complete answer wrt the Cholesky case.

I have the matrix $\Sigma=LL^T$. Is there a way of getting $Tr(\Sigma^{-1})$ without using the SVD? I'm guessing eigen decomposition is just as costly as SVD. I have already computer the lower Cholesky matrix $L$ for a previous computation.

The matrix is symmetric, positive definite and (unfortunately) dense.

Best Answer

Let $\alpha = \mathrm{trace}(LL^t)$. We are avoiding your notation $\Sigma$ since it often considered as a singular value matrix.

There are two options. The first method is to compute all singular values of $L$ without computing singular vectors. We note that the computation of singular vectors are expensive. Then the required trace is given by $$ \alpha = \sum_{i=1}^n \sigma_i^{-2} $$ where $\sigma_i$ are singular values of $L$. Of course, we assume that there are no tiny singular values (compared with $\sigma_1$) which are unsuitable for reciprocation. The positive definiteness of $LL^t$ precludes such worries.

The second option is to invert the matrix $L$ which can be done quite efficiently since it is triangular. We do not compute any singular values but they should not be tiny as previously. If there are, then they will have the effect of creating large diagonal elements in the inverse. In the worst case, there will be overflows in arithmetic operations. However, since $LL^t$ is assumed to be positive definite it is very unlikely. If $M = L^{-1}$ then the required trace is given by $$ \alpha = \mathrm{trace}(M^tM) = \mathrm{trace}(MM^t) = \|M\|_2^2 $$ where $\|.\|_{\mathrm F}$ is the Frobenius norm.

Additional notes (as requested):

When we compute the singular value decomposition of a matrix $A$ in the form $A = U\Sigma V^t$, the matrix $A$ is diagonalised using orthogonal transformations and these transformations are necessary to compute $U$ and $V$. There are many applications where the singular vectors are NOT required. In that case, we use orthogonal transformations to diagonalise the matrix but we do not accumulate those transformations to get $U$ and $V$. Standard software libraries have this option. As an example, please google "LAPACK". This is a public domain library in Fortran but now it has been translated to many leading languages. The LAPACK is the gold standard for the SVD.

The inversion of an upper (or lower) tridiagonal matrix $L$ is a standard problem taught at undergraduate level (or even sometimes at high school). We solve the equations $$ L x_i = e_i $$ for $i=1,\ldots,n$ where $e_i$ is the $i$th column of the $n \times n$ identity matrix $I$. Then if you form the matrix $X$ from the columns of $x_i$ then $X$ is the inverse of $L$. If $L$ is upper (lower) diagonal then it has the same shape as $L$.