Can I get a Cholesky decomposition from the inverse of a matrix?

Question

I have the inverse of a giant covariance matrix from which I'd like to draw random instances.

The way I know how to do this is to do a Cholesky decomposition on the covariance matrix and use it to transform a vector of independent Gaussians. So the straightforward process is to invert the giant matrix and then do the Cholesky decomposition.

Since Cholesky decomposition can be used to get inverses in the first place, it seems like there might be a way to get what I need without the inversion step (which takes a long time).

Question: Is there a way to use Cholesky decomposition to do this without inverting the matrix?

Answer 1

You can avoid inverting the matrix by generating draws by means of the eigendecomposition method. According to this method, the draws are generated by doing this product: $$ (V D)^\top X^\top \,, $$

where $V$ is the eigenvectors of the matrix, $D$ is a diagonal matrix containing the square roots of the eigenvalues and $X$ is a matrix containing draws from the standard univariate $N(0,1)$ distribution.

It is straightforward to adapt this method and avoid recovering the original covariance matrix by using these results: 1) the eigenvalues of a matrix $A$ are the reciprocal of the eigenvalues of its inverse $A^{-1}$; 2) the eigenvectors of a matrix A are also eigenvectors of its inverse $A^{-1}$.

---

Example:

Let's say that the original covariance matrix is the following:

$$ A = \left[ \begin{array}{rrr} 1 & 0.8 & -0.4 \\ 0.8 & 2 & 0.3 \\ -0.4 & 0.3 & 3 \end{array} \right] \,. $$ But you have the inverse of this matrix, $B=A^{-1}$: $$ A^{-1}= B = \left[ \begin{array}{rrr} 1.699 & -0.725 & 0.299 \\ -0.725 & 0.817 & -0.178 \\ 0.299 & -0.178 & 0.391 \end{array} \right] \,. $$

The eigendecomposition method based on the original matrix $A$ yields the following covariance matrix for a sample of draws:

A <- rbind(c(1,0.8,-0.4), c(0.8,2,0.3), c(-0.4,0.3,3))
e1 <- eigen(A, symmetric=TRUE)
set.seed(1)
X <- matrix(rnorm(5000*ncol(A)), ncol=ncol(A))
draws1 <- t(e1$vectors %*% sqrt(diag(e1$values)) %*% t(X))
draws1.cov <- cov(draws1)
draws1.cov
#           [,1]      [,2]       [,3]
#[1,]  0.9765023 0.8030752 -0.3970233
#[2,]  0.8030752 1.9941052  0.3229827
#[3,] -0.3970233 0.3229827  3.1689348

Using the matrix that you have (the inverse of A), you just need to invert the eigenvalues:

B <- solve(A)
e2 <- eigen(B, symmetric=TRUE)
e2$values <- 1/e2$values
draws2 <- t(e2$vectors %*% sqrt(diag(e2$values)) %*% t(X))
draws2.cov <- cov(draws2)
draws2.cov
#           [,1]      [,2]       [,3]
#[1,]  0.9765023 0.8030752 -0.3970233
#[2,]  0.8030752 1.9941052  0.3229827
#[3,] -0.3970233 0.3229827  3.1689348

A covariance matrix closely matching the original one is obtained and we didn't need to invert B in order to recover the original covariance matrix A.

---

A small simulation to check the validity of this approach:

set.seed(3)
niter <- 1000
m <- matrix(0, nrow=ncol(A), ncol=ncol(A))
for (i in seq_len(niter))
{
    X <- matrix(rnorm(5000*ncol(A)), ncol=ncol(A))
    draws2 <- t(e2$vectors %*% sqrt(diag(e2$values)) %*% t(X))
    m <- m + cov(draws2)
}
m/niter
# average covariance matrix
#           [,1]      [,2]       [,3]
#[1,]  1.0005129 0.7995872 -0.4005644
#[2,]  0.7995872 1.9993231  0.2990850
#[3,] -0.4005644 0.2990850  2.9957277
# original covariance matrix 'A'
#    [,1] [,2] [,3]
#[1,]  1.0  0.8 -0.4
#[2,]  0.8  2.0  0.3
#[3,] -0.4  0.3  3.0

We can see that the covariance matrix of the draws are on average very close to the original covariance matrix A.