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, the…