# Efficiently calculating leave-one-out conditional multivariate normal distributions

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.

• I vote to keep open - if anything, the question should imho go to math exchange and not SO, as the underlying problem (on which my answer draws directly) is a mathematical one, viz the Woodbury formula. Commented Oct 17, 2022 at 16:10
• Thanks. I asked here because I wasn't sure if the linear algebra update method was the correct approach to this (maybe there's some obscure statistical one-liner for this that someone on CrossValidated would know?), but if a mod wants to move it that's fine.
– zkxg
Commented Oct 17, 2022 at 16:22

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

• Brilliant, this seems to do the trick. I get a ~14x speedup with n = 1000. The one note I would add to this is that it seems necessary to add some strategic parentheses in the last line so that R does the multiplication efficiently i.e. Sigma.inv + Sigma.inv %*% U %*% tmp %*% (V %*% Sigma.inv)
– zkxg
Commented Oct 17, 2022 at 17:38
• I am glad it is useful - and surprised about the need for the parentheses! Commented Oct 18, 2022 at 4:30