1.1 MVN density computation and random number generation

The two computationally intensive steps in evaluating \(\log f(x)\) are computing \(\log|\Sigma|\), and \(z^\top z\), without explicitly inverting \(\Sigma\) or repeating mathematical operations. But once we have \(z\), computing the dot product \(z^\top z\) is cheap. How to perform these steps efficiently in practice depends on whether the covariance matrix \(\Sigma\), or the precision matrix \(\Sigma^{-1}\) is available. For both cases, we start by finding a lower triangular matrix root: \(\Sigma=LL^\top\) or \(\Sigma^{-1}=\Lambda\Lambda^\top\). Since \(\Sigma\) and \(\Sigma^{-1}\) are positive definite, we will use the Cholesky decomposition, which is the unique matrix root with all positive elements on the diagonal.

With the Cholesky decomposition in hand, we compute the log determinant of \(\Sigma\) by adding the logs of the diagonal elements of the factors. \[\begin{aligned} \label{eq:logDet} \log|\Sigma|= \begin{cases} \phantom{-}2\sum_{m=1}^M\log L_{mm}&\text{ when $\Sigma$ is given}\\ -2\sum_{m=1}^M\log \Lambda_{mm}&\text{ when $\Sigma^{-1}$ is given} \end{cases}\end{aligned}\]

Having already computed the triangular matrix roots also speeds up the computation of \(z^\top z\). If \(\Sigma^{-1}\) is given, \(z=\Lambda^\top(x-\mu)\) can be computed efficiently as the product of an upper triangular matrix and a vector. When \(\Sigma\) is given, we find \(z\) by solving the lower triangular system \(Lz=x-\mu\). The subsequent \(z^\top z\) computation is trivially fast.

The algorithm for simulating \(X\sim\mathbf{MVN}\!\left(\mu,\Sigma\right)\) also depends on whether \(\Sigma\) or \(\Sigma^{-1}\) is given. As above, we start by computing the Cholesky decomposition of the given covariance or precision matrix. Define a random variable \(Z\sim\mathbf{MVN}\!\left(0,I_M\right)\), and generate a realization \(z\) as a vector of \(M\) samples from a standard normal distribution. If \(\Sigma\) is given, then evaluate \(x=Lz+\mu\). If \(\Sigma^{-1}\) is given, then solve for \(x\) in the triangular linear system \(\Lambda^\top\left(x-\mu\right)=z\). The resulting \(x\) is a sample from \(\mathbf{MVN}\!\left(\mu,\Sigma\right)\). The following shows that this approach is correct:

\[\begin{aligned} \mathbf{E}\!\left(X\right)&=\mathbf{E}\!\left(LZ+\mu\right)=\mathbf{E}\!\left(\Lambda^\top Z+\mu\right)=\mu\\ \mathbf{cov}\!\left(X\right)&= \mathbf{cov}\!\left(LZ+\mu\right)=\mathbf{E}\!\left(LZZ^\top L^\top\right)=LL^\top=\Sigma\\ \mathbf{cov}\!\left(X\right)&=\mathbf{cov}\!\left((\Lambda^\top)^{-1}Z+\mu\right)=\mathbf{E}\!\left((\Lambda^\top)^{-1}ZZ^\top\Lambda^{-1}\right) =(\Lambda^\top)^{-1}\Lambda^{-1}=(\Lambda\Lambda^\top)^{-1}=\Sigma \end{aligned}\]

These algorithms apply when the covariance/precision matrix is either sparse or dense. When the matrix is dense, the computational complexity is \(\mathcal{O}\!\left(M^3\right)\) for a Cholesky decomposition, and \(\mathcal{O}\!\left(M^2\right)\) for either solving the triangular linear system or multiplying a triangular matrix by another matrix (Golub and Van Loan 1996). Thus, the computational cost grows cubically with \(M\) before the decomposition step, and quadratically if the decomposition has already been completed. Additionally, the storage requirement for \(\Sigma\) (or \(\Sigma^{-1}\)) grows quadratically with \(M\).

1.2 Sparse matrices in R

The Matrix package (Bates and Maechler 2017) defines various classes for storing sparse matrices in compressed formats. The most important class for our purposes is dsCMatrix, which defines a symmetric matrix, with numeric (double precision) elements, in a column-compressed format. Three vectors define the underlying matrix: the unique nonzero values (just one triangle is needed), the indices in the value vector for the first value in each column, and the indices of the rows in which each value is located. The storage requirements for a sparse \(M\times M\) symmetric matrix with \(V\) unique nonzero elements in one triangle are for \(V\) double precision numbers, \(V+M+1\) integers, and some metadata. In contrast, a dense representation of the same matrix stores \(M^2\) double precision values, regardless of symmetry and the number of zeros. If \(V\) grows more slowly than \(M^2\), the matrix becomes increasingly sparse (a smaller percentage of elements are nonzero), with greater efficiency gains from storing the matrix in a compressed sparse format.

2.1 An example

Suppose we want to generate samples from an MVN approximation to the posterior distribution of our example model from Section 1.2. sparseMVN includes functions to simulate data for the example (binary.sim()), and to compute the log posterior density (binary.f()), gradient (binary.grad()), and Hessian (binary.hess()). trustOptim::trust.optim() in trustOptim is a nonlinear optimizer that estimates the curvature of the objective function using a sparse Hessian.

The calls to trustOptim::trust.optim() return the posterior mode, and the Hessian evaluated at the mode. These results serve as the mean and the negative precision of the MVN approximation to the posterior distribution of the model.

We can then sample from the posterior using an MVN approximation, and compute the MVN log density for each sample.

samples <- rmvn.sparse(R, pm, CH, prec=TRUE)
logf <- dmvn.sparse(samples, pm, CH, prec=TRUE)

The ability to accept a precision matrix, rather than having to invert it to a covariance matrix, is a valuable feature of sparseMVN. This is because the inverse of a sparse matrix is not necessarily sparse. In the following chunk, we invert the Hessian, and drop zero values to maintain any remaining sparseness. Note that there are 10,404 total elements in the Hessian.

Matrix::nnzero(H)
Hinv <- Matrix::drop0(solve(H))
Matrix::nnzero(Hinv)

Nevertheless, we should check that the log densities from [dmvn.sparse] correspond to those that we would get from dmvnorm.

logf_dense <- mvtnorm::dmvnorm(samples, pm, as.matrix(Hinv), log=TRUE)
all.equal(logf, logf_dense)