1
$\begingroup$

Assumptions

Let's start with the following assumptions. Assume the matrices $\Sigma$, $L$, and $S$ defined below are known in advance.

  • $\Sigma$ is a symmetric positive definite $n \times n$ matrix.
  • $\Sigma$ has Cholesky factorization $L L^T$ where $L$ is lower triangular.
  • Matrix $S$ is a subset of unique rows of the $n \times n$ identity matrix.
    • $S$ has $m$ rows and $n$ columns, where $m < n$.
    • Each row of $S$ has a single 1 and all the other elements equal to 0.
    • The rows of $S$ are orthonormal.
  • By construction, the product $S \Sigma$ is an $m \times n$ matrix whose rows are a unique subset of the rows of $\Sigma$.

Question

Is there a fast shortcut to get the Cholesky factor of $S \Sigma S^T$ using $L$ and $S$? Looking for something much quicker and easier to compute than the standard Cholesky decomposition algorithms which uses the existing Cholesky decomposition of $\Sigma$.

Context

I am asking because I want to compute the marginal density of a multivariate normal random variable using https://statproofbook.github.io/P/mvn-marg.html. I already have the Cholesky factor of the covariance of the full multivariate normal distribution, and my Stan code will run much faster if I do not need to compute an additional Cholesky factor for the marginal distribution.

$\endgroup$
8
  • $\begingroup$ I don't know how to do what you're asking, but if you're trying to compute the inverse of the population covariance for a subset of the random variables and you already have the inverse of the population covariance for all of the variables, these matrices are related via the Schur Complement $\endgroup$ Commented Mar 27, 2023 at 20:44
  • $\begingroup$ How large is that correlation matrix $\Sigma$ ? Unless this goes into the thousands it is hard to believe that the standard Cholesky algorithm is so slow that one does not want to perform it on a submatrix. $\endgroup$
    – Kurt G.
    Commented Mar 28, 2023 at 7:13
  • $\begingroup$ Thanks for the tip about the Schur complement. $\endgroup$
    – landau
    Commented Mar 28, 2023 at 15:09
  • $\begingroup$ As for the size of $\Sigma$, I expect it to only be around 4x4 on average, maybe 8x8 in some cases. The issue is not the size of the matrix, but rather the sheer number of Cholesky factorizations over the course of the MCMC algorithm. $\endgroup$
    – landau
    Commented Mar 28, 2023 at 15:09
  • $\begingroup$ The MCMC has at least 4000 iterations, and each iteration has a unique value of $\Sigma$. Within each iteration, each multivariate normal observation in the data has its own $S$ matrix, and there are upwards of 100 observations. Without a shortcut, that multiplies out to 400000 Cholesky factorizations of $S \Sigma S^T$ instead of just 4000, which I have observed seriously impacts the speed of my Stan code. $\endgroup$
    – landau
    Commented Mar 28, 2023 at 15:10

1 Answer 1

1
$\begingroup$

EDIT 2023-04-17

The answer below is wrong because $R \Sigma R$ is not lower-triangular in the general case. I leave it here for posterity. Down-votes may help ensure people know not to trust it.

Original post (which is incorrect)

Thanks for the nudges to explain more about my issue. The overarching problem is a statistical issue regarding multivariate normal distributions. Suppose random variable $y$ is a vector of length $n$ with a multivariate normal distribution with mean vector $\mu$ and positive-definite symmetric covariance matrix $\Sigma$. We denote this by:

$$ y \sim \text{MVN}(\mu, \Sigma) $$

My goal is to efficiently work with the marginal distribution of a subset of the components of $y$:

$$ S R y \sim \text{MVN} ( SR\mu, SR\Sigma R^TS^T) $$

In this problem, we assume we already know the Cholesky factorization $\Sigma = L L^T$. In order to quickly evaluate the density of the marginal distribution above, we want a shortcut to get the Cholesky factor of $SR\Sigma R^TS^T$ without having to run the standard factorization algorithm from scratch. This is because multi_normal_cholesky_lpdf() is much faster than multi_normal_lpdf() in Stan. Fortunately, a shortcut exists, as I prove below.

Here, $R$ is an $n \times n$ permutation matrix. In my case, $Ry$ is a column vector that puts the desired elements of the subset of $y$ at the front of the vector, while preserving the relative order of the desired components. To explain more about what I mean, suppose $y = [a \ \ b \ \ c \ \ d]^T$ and I want the marginal distribution of $[b \ \ c]^T$. Then, $Ry = [b \ \ c \ \ a \ \ d]^T$, with $b$ and $c$ at the front and $b$ before $c$.

Also, $S$ is an $m \times n$ matrix, where $m < n$ is the cardinality of the desired subvector $S R y$. Here, $S = [I_m \ \ 0_{m \times (n - m)}]$, where $I_m$ is the $m \times m$ identity matrix and $0_{m \times (n - m)}$ is the appropriately sized matrix of zeroes.

Now, let $L$ be the lower-triangular Cholesky factor of $\Sigma$ (with $\Sigma = L L^T$). We know $R L R^T$ is lower triangular because $R$ permutes the rows of $L$ in the exact same way as $R^T$ permutes the columns of $L$. (The previous sentence is definitely wrong, and it is straightforward to come up with small counterexamples.) In addition, $R^T = R^{-1}$ because $R$ is a permutation matrix and thus orthonormal. Thus:

$$ R L R^T(R L R^T)^T = R L R^T R L^T R^T = R L L^T R^T = R \Sigma R^T $$

That means $R L R^T$ is the Cholesky factor of $R \Sigma R^T$. Finally, the matrix $S$ is constructed such that $S R \Sigma R^T S^T$ is the upper-left $m \times m$ block of $R \Sigma R^T$ and $S R L R^T S^T$ is the upper-left $m \times m$ block of $R \Sigma R^T$. By the lemma below, this makes $S R L R^T S^T$ the lower-triangular Cholesky factor of $S R \Sigma R^T S^T$. Best of all, $S R L R^T S^T$ is super fast to compute because it is just a subset of the rows and columns of $L$.

Lemma

Define $\Sigma$ to be an $n \times n$ positive definite symmetric matrix with lower-triangular covariance matrix $L$. Write these matrices in block notation:

$$ \Sigma = \begin{bmatrix} A & B \\ C & D \end{bmatrix} \qquad L = \begin{bmatrix} X & O \\ Y & Z \end{bmatrix} $$

where $A$ and $X$ are $m \times m$ blocks $(m < n)$ and $O$ is the appropriately sized matrix of zeroes. Then, $X$ is the $m \times m$ lower-triangular Cholesky factor of $m \times m$ matrix $A$.

Proof of the Lemma

By construction, matrix $X$ is already lower-triangular. From matrix multiplication and the definition of the Cholesky factorization, we find:

$$ \begin{aligned} \Sigma &= L L^T \\ \begin{bmatrix} A & B \\ C & D \end{bmatrix} &= \begin{bmatrix} X & O \\ Y & Z \end{bmatrix} \begin{bmatrix} X^T & Y^T \\ O^T & Z^T \end{bmatrix} = \begin{bmatrix} XX^T + OO^T & XY^T + OZ^T \\ YX^T + ZO^T & YY^T + ZZ^T \end{bmatrix} = \begin{bmatrix} XX^T & XY^T + OZ^T \\ YX^T + ZO^T & YY^T + ZZ^T \end{bmatrix} \end{aligned} $$

Thus, $A = XX^T$.

$\endgroup$
4
  • $\begingroup$ How much does that speed up your code overall? $\endgroup$
    – Kurt G.
    Commented Mar 30, 2023 at 21:14
  • $\begingroup$ It’s running right now, and it looks to be at least 10x faster, but still pretty slow. $\endgroup$
    – landau
    Commented Mar 30, 2023 at 21:58
  • $\begingroup$ It could be some other problem in the MCMC, I am confident this exercise was worthwhile. $\endgroup$
    – landau
    Commented Mar 31, 2023 at 1:40
  • $\begingroup$ It surely was. Just being curious. 10 times faster sounds pretty good. Sorting algorithms are notorious bottlenecks. That's why I was asking. $\endgroup$
    – Kurt G.
    Commented Mar 31, 2023 at 6:38

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .