Calculating the log-determinant of large, sparse covariance matrices 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:


*

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

*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!
 A: 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$).
A: 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$.
