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!

  • $\begingroup$ There are some nice tricks if all diagonal entries of D are equal. Any chance that's the case? $\endgroup$ – user20160 Oct 13 '17 at 9:48
  • $\begingroup$ @user20160 Ideally they'd be allowed to be different, but I'd be quite happy to make the approximation that $\mathbf{D} = c \mathbf{I}$ if it helps get to a solution! $\endgroup$ – CBowman Oct 13 '17 at 10:26

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$).

| cite | improve this answer | |

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$.

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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