Efficiently calculating leave-one-out conditional multivariate normal distributions

Question

I have a multivariate normal distribution for vector $\mathbf{x}$ with mean vector $\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$. In my specific use-case, $\boldsymbol{\Sigma}$ is actually a correlation matrix. For the sake of calculating some leave-one-out log-likelihood values downstream, I need to efficiently calculate conditional distributions for each dimension (i.e. get the conditional distribution for $\mathbf{x_1}$ when setting the other values to a pre-specified vector $\mathbf{a}$). As shown on the relevant wikipedia page, these are:

$$\bar{\boldsymbol\mu} = \boldsymbol\mu_1 + \boldsymbol\Sigma_{12} \boldsymbol\Sigma_{22}^{-1} \left( \mathbf{a} - \boldsymbol\mu_2 \right) \\ \overline{\boldsymbol\Sigma} = \boldsymbol\Sigma_{11} - \boldsymbol\Sigma_{12} \boldsymbol\Sigma_{22}^{-1} \boldsymbol\Sigma_{21}$$

Where $\boldsymbol{\Sigma}_{11}$, $\boldsymbol{\Sigma}_{12}$, $\boldsymbol{\Sigma}_{21}$ and $\boldsymbol{\Sigma}_{22}$ represent sub-blocks of $\boldsymbol{\Sigma}$.

In my case I need every leave-one-out conditional distribution, so repeatedly computing $\boldsymbol\Sigma_{22}^{-1}$ becomes computationally burdensome. I think the complexity for one inverse (which is feasible) is $O((n-1)^3)$, so $n$ of those is $O(n(n-1)^3)$ which is too expensive. For my purposes, $n$ gets up to the neighborhood of ~5000.

This might be nothing more than a simple linear algebra problem. My instinct is to compute $\boldsymbol{\Sigma}^{-1}$ once at first, then repeatedly "subtract off" the influence from $\boldsymbol{\Sigma}_{11}$, $\boldsymbol{\Sigma}_{12}$, and $\boldsymbol{\Sigma}_{21}$ for each partitioning, but I can't see the path to doing that "subtraction".

A simple example of what I need in R for $n = 3$ and the first partitioning:

> Sigma = matrix(c(1, .15, .2, .15, 1, .15, .2, .15, 1), nrow = 3)
> Sigma
     [,1] [,2] [,3]
[1,] 1.00 0.15 0.20
[2,] 0.15 1.00 0.15
[3,] 0.20 0.15 1.00
> solve(Sigma) # What I have
           [,1]       [,2]       [,3]
[1,]  1.0579004 -0.1298701 -0.1920996
[2,] -0.1298701  1.0389610 -0.1298701
[3,] -0.1920996 -0.1298701  1.0579004
> solve(Sigma[-1,-1]) # What I need to calculate quickly without using solve() again
           [,1]       [,2]
[1,]  1.0230179 -0.1534527
[2,] -0.1534527  1.0230179
> # and so on for Sigma[-2,-2] and Sigma[-3,-3]

Any help is much appreciated.

Answer 1

Based on the Woodbury formula, this answer suggests the following code:

Sigma <- matrix(c(1, .15, .2, .15, 1, .15, .2, .15, 1), nrow = 3)

Sigma.inv <- solve(Sigma)               # what you have
(Sigma23.inv <- solve(Sigma[-1,-1]))    # this what is to be sped up

# ingredients to Woodbury
U <- cbind(c(1,rep(0,2)), c(0,Sigma[2:3,1])) 
V <- rbind(c(0,Sigma[1,2:3]), c(1,rep(0,2)))

solve(Sigma-U%*%V)[-1,-1] # effectively still inverting the submatrix, just a check that we get the same result

tmp <- solve(diag(2)-V%*%Sigma.inv%*%U) # a cheap 2x2 inversion
Sigma.inv + Sigma.inv%*%U%*%tmp%*%V%*%Sigma.inv # now entirely in terms of available big inverse plus a cheap 2x2 inversion