1
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$ Oct 17, 2022 at 16:10
  • $\begingroup$ 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. $\endgroup$
    – zkxg
    Oct 17, 2022 at 16:22

1 Answer 1

0
$\begingroup$

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
$\endgroup$
2
  • 1
    $\begingroup$ 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) $\endgroup$
    – zkxg
    Oct 17, 2022 at 17:38
  • $\begingroup$ I am glad it is useful - and surprised about the need for the parentheses! $\endgroup$ Oct 18, 2022 at 4:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.