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.