# Numerical trace of inverse matrix from Cholesky

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.

• What about computing an upper triangular Cholesky factor instead? Then, from the shape of $\Sigma^{-1}=U^{-T}U^{-1}$, you get that the diagonal entries are the squares of the inverse of the diagonal of $U$. – Federico Poloni Apr 15 '14 at 6:20
• Oh, and for a SPD matrix, essentially the eigendecomposition is the SVD. – Federico Poloni Apr 15 '14 at 6:21
• Theoretically, yes. Unfortunately my matrix has some eigen values close to zero that sometimes it dips below zero. Basically numerical problems. For the first comment, I guess I was hoping not to use inverse of U (can do the same operation with L) – sachinruk Apr 15 '14 at 7:47
• So you wish to compute $\mbox{trace}(L^{-T}L^{-1})$ without actually computing $L^{-1}$? You can use quadrature + randomisation based algorithms to approximate this term, but I don't know if that's what you're after.... – Suvrit Apr 15 '14 at 10:05

## 1 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$.

• welcome! Though this is the beginning of an answer, it still lack a lot of details - what is the right and efficient way to invert a triangular matrix? Same goes for computing the eigenvalues without the eigen-vectors? As it is, the OP can't use your answer. – Amir Sagiv Mar 15 '17 at 8:27
• Amir Sagiv : I have revised my Answer so that almost anybody can understand. – Vini Mar 15 '17 at 10:51
• Additionally, it should be noted that $M$ is also lower triangular (because $L$ is), so that $||M||_{F}$ can be calculated on just the lower triangular portion of the matrix. – aepound May 23 '17 at 14:47