Dealing with the inverse of a positive definite symmetric (covariance) matrix?

Question

In statistics and its various applications, we often calculate the covariance matrix, which is positive definite (in the cases considered) and symmetric, for various uses. Sometimes, we need the inverse of this matrix for various computations (quadratic forms with this inverse as the (only) center matrix, for example). Given the qualities of this matrix, and the intended uses, I wonder:

What is the best, in terms of numerical stability, way to go about computing or using (let's say for quadratic forms or matrix-vector multiplication in general) this inverse? Some factorization that can come in handy?

Answer 1

The Cholesky factorization $C=R^TR$ leads to a Cholesky-like factorization of the inverse $C^{-1}=SS^T$ with the upper triangular matrix $S=R^{-1}$.

In practice, is best to keep the inverse factored. If $R$ is sparse then it is usually even better to keep $S$ implicit, as matrix-vector products $y=C^{-1}x$ can be computed by solving the two triangular systems $R^Tz=x$ and $Ry=z$.

Answer 2

A Cholesky factorization makes the most sense for the best stability and speed when you are working with a covariance matrix, since the covariance matrix will be positive semi-definite symmetric matrix. Cholesky is a natural here. BUT...

IF you intend to compute a Cholesky factorization, before you ever compute the covariance matrix, do yourself a favor. Make the problem maximally stable by computing a QR factorization of your matrix. (A QR is fast too.) That is, if you would compute the covariance matrix as

$$ C = A^{T} A $$

where $A$ has had the column means removed, then see that when you form $C$, it squares the condition number. So better is to form the QR factors of $A$ rather than explicitly computing a Cholesky factorization of $A^{T}A$.

$$ A = QR $$

Since Q is orthogonal,

$$ \begin{align} C &= (QR)^{T} QR \\ &= R^T Q^T QR \\ &= R^T I R \\ &= R^{T} R \end{align} $$

Thus we get the Cholesky factor directly from the QR factorization, in the form of $R^{T}$. If a $Q$-less QR factorization is available, this is even better since you don't need $Q$. A $Q$-less QR is a fast thing to compute, since $Q$ is never generated. It becomes merely a sequence of Householder transformations. (A column pivoted, $Q$-less QR would logically be even more stable, at the cost of some extra work to choose the pivots.)

The great virtue of using the QR here is it is highly numerically stable on nasty problems. Again, this is because we never had to form the covariance matrix directly to compute the Cholesky factor. As soon as you form the product $A^{T}A$, you square the condition number of the matrix. Effectively, you lose information down in the parts of that matrix where you originally had very little information to start with.

Finally, as another response points out, you don't even need to compute and store the inverse at all, but use it implicitly in the form of backsolves on triangular systems.

Answer 3

I did this for the first time recently, using suggestions from mathSE.

SVD was recommended by most I think, but I opted for the simplicity of Cholesky:

If matrix $M = A A^\top$, then I decompose $M$ to a triangular matrix $L$ using Cholesky, such that $M = L L^\top$. I then use backsubstitution or forwardsubstitution (depending on whether I choose L to be upper or lower triangular), to invert $L$, such that I have $L^{-1}$. From this, I can quickly calculate $M^{-1} = \left(L L^\top\right)^{-1} = L^{-\top}L^{-1}$.

---

Start with:

$M = A A^\top$, where $M$ is known and is implicitly symmetric and is also positive-definite.

Cholesky factorisation:

$M \rightarrow L L^\top$, where $L$ is square and non-singular

Back-substitution:

$L \rightarrow L^{-1}$, probably the fastest way to invert $L$ (don't quote me on that though)

Multiplication:

$M^{-1} = \left(L L^\top\right)^{-1} = L^{-\top} L^{-1}$

Notation used: Lower indices are rows, upper indices are columns and $L^{-\top}$ is the transpose of $L^{-1}$

---

My Cholesky algorithm (probably from Numerical Recipes or Wikipedia)

$L_i^j = \frac{M_i^j - M_i \cdot M_j}{M_i^i - M_i \cdot M_i}$

This can almost be done in-place (you only need temporary storage for the diagonal elements, an accumulator and some integer iterators).

---

My back-substitution algorithm (from Numerical Recipes, check their version as I may have made a mistake with the LaTeX markup)

$\left(L^{-1}\right)_i^j = \left\{\begin{array}{11} 1 / {L_i^i} & \mbox{if } i = j\\ \left(-L_i \cdot \left(L^{-T}\right)_j\right) / L_i^i & \mbox{otherwise} \end{array}\right.$

As $L^{-T}$ appears in the expression, the order that you iterate over the matrix is important (some parts of the result matrix depend on other parts of it that must be calculated beforehand). Check the Numerical Recipes code for a complete example in code. [Edit]: Actually, just check the Numerical Recipes example. I've over-simplified too much by using dot-products, to the point that the above equation has a cyclic dependency no matter what order you iterate...

Answer 4

If you know that the matrix has an inverse (i.e., if it is indeed positive definite) and if it isn't too large, then the Cholesky decomposition gives an appropriate means to characterize the inverse of a matrix.