Calculating the log-determinant of large, sparse covariance matrices

Question

As part of calculating the log-density of a very large multivariate normal distribution, I need to evaluate the following log-determinant: $$ f = \mathrm{ln}\left( \mathrm{det}(\mathbf{XAX}^\top + \mathbf{D}) \right), $$ where:

• $\mathbf{A}$ is $n \times n$, symmetric, sparse and positive-definite,
• $\mathbf{D}$ is $m \times m$, diagonal, sparse and positive-definite,
• $\mathbf{X}$ is $m \times n$, sparse and has positive entries.

The main issue here is that $m \approx 10^5$, and I've been struggling with memory/stability issues.

There are two parts to my question:

1. Does anyone have suggestions of how to approach a problem of this scale numerically? (particular implementations of algorithms ect)

2. Are there any tricks we can play to make the problem more tractable using the fact that $n \approx 10^3$, and therefore $ n \ll m$?

Regarding #2, I've been trying tricks like finding the determinant of $(\mathbf{XAX}^\top + \mathbf{D})^{-1}$ instead (since $\mathrm{det}(\mathbf{C}^{-1}) = \mathrm{det}(\mathbf{C})^{-1}$) and the applying the Woodbury identity, or decomposing $\mathbf{A} = \mathbf{LL}^{\top}$, but I haven't made much progress.

This has been causing me real problems, so any help is much appreciated, thanks!

Answer 1

Here's one approach that assumes all elements along the diagonal of $D$ are equal.

First, let's combine $X$ and $A$. Find some matrix $B$ such that $B^T B = A$. You can do this using the Cholesky decomposition, eigendecomposition, etc. (which is easy because $A$ is relatively small). Let $V = X B^T$, so $V V^T = X A X^T$. Now the problem is to compute $\log \det(V V^T + D)$.

One way to calculate the determinant is as the product of the eigenvalues, so we need to find the eigenvalues of $V V^T + D$.

Let $\lambda_1, \dots, \lambda_n$ denote the eigenvalues of $V V^T$. We've assumed that $D = cI$, where $c$ is a constant and $I$ is the identity matrix. Because of this special form for $D$, the eigenvalues of $V V^T + D$ are given by $\lambda_1 + c, \dots, \lambda_n + c$. Unfortunately, I don't think there's any convenient trick like this for more general forms of $D$.

Now we must find the eigenvalues of $V V^T$. $V$ has size $n \times m$ and, because $m \ll n$, $V V^T$ has at most $m$ nonzero eigenvalues. Instead of operating directly on $V V^T$ (size $n \times n$), we can find the eigenvalues using the singular value decomposition of $V$, which should be much more efficient. Let $s_1, \dots, s_m$ denote the singular values of $V$. Then:

$$\lambda_i = \left \{ \begin{array}{cl} s_i^2 & 1 \le i \le m \\ 0 & i > m \\ \end{array} \right .$$

Putting everything together, we have:

$$\begin{array}{lcl} \log \det(X A X^T + D) & = & \log \prod_{i=1}^n (\lambda_i + c) \\ & = & (n-m) \log c + \sum_{i=1}^m \log(s_i^2 + c) \end{array}$$

The main computational work is a Cholesky or eigendecomposition of $A$ (size $m \times m$) and a singular value decomposition of $V$ (size $n \times m$). Because $m \ll n$, this should be much more efficient than directly working with the matrix $X A X^T + D$ (size $n \times n$).

Answer 2

If we are able to decompose $\mathbf{A} = \mathbf{LL}^{\top}$, and let $\mathbf{V} = \mathbf{XL}$ then $$ \mathrm{det}(\mathbf{XAX}^\top + \mathbf{D}) = \mathrm{det}(\mathbf{VV}^\top + \mathbf{D}) $$ Now using Sylvester's determinant theorem we may write $$ \mathrm{det}(\mathbf{VV}^\top + \mathbf{D}) = \mathrm{det}(\mathbf{D})\mathrm{det}(\mathbf{V}^\top \mathbf{D}^{-1}\mathbf{V} + \mathbf{I}_n) $$ Calculating $\mathrm{det}(\mathbf{D})$ is trivial, so the problem has been reduced to finding the determinant of the $n \times n$ matrix $\mathbf{V}^\top \mathbf{D}^{-1}\mathbf{V} + \mathbf{I}_n$, which represents a significant reduction in complexity as $n \ll m$.