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 $n \times m$, 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!