4
$\begingroup$

I am running code in R where I repeatedly evaluate the multivariate normal density of a high-dimensional object for thousands of sequential iterations. I am using dmvnorm() to do this. My understanding is that, under the hood, dmvnorm() computes the inverse of the covariance matrix to evaluate the MVN density. But I already compute the inverse of the covariance matrix elsewhere in my code for other purposes.

I don't see a way to pass the inverse of the covariance matrix directly to the dmvnorm() function, but are there any alternative functions that would allow me to do so? Or, is it possible for me to rewrite the dmvnorm() function myself using the Cholesky decomposition?

$\endgroup$
3
  • 2
    $\begingroup$ yes this would be a huge speedup; see search.r-project.org/CRAN/refmans/LaplacesDemon/html/… $\endgroup$ Jun 20, 2023 at 17:26
  • 3
    $\begingroup$ See also dmvn() in package mvn fast, which allows you to pass in the Cholesky decomposition of the covariance matrix if you already have it, otherwise it uses a Cholesky decomposition of the covariance matrix internally. $\endgroup$ Jun 20, 2023 at 17:49
  • 3
    $\begingroup$ Because dmvnorm is so simple -- its logarithm is a quadratic form -- why not just code it yourself? It's basically one line. Consider passing both the inverse and its log determinant as arguments so that you can precompute the normalizing factor, too. Work with the logarithms for as long as possible because in high dimensions overflow and underflow of the densities is almost unavoidable. $\endgroup$
    – whuber
    Jun 20, 2023 at 21:37

1 Answer 1

2
$\begingroup$

This version of dmvnorm requires you to supply the inverse covariance matrix (aka precision matrix) and the logarithm of the determinant of the precision matrix (which is needed to normalize the density). It returns the density at $x-\mu$ (optionally, the log density).

dmvnorm. <- function(x, Sigma.1, log.det.Sigma, mu = rep(0, length(x)), log.p = FALSE) {
  q <- -(x - mu) %*% Sigma.1 %*% (x - mu) / 2 + log.det.Sigma / 2 - log(2 * pi) * length(x)/2
  if(isTRUE(log.p)) q else exp(q)
}

You can see it is straightforward, because the log density for a $p$-vector $x$ is

$$\log f(x;\mu,\Sigma) = -\left((x-\mu)\Sigma^{-1}(x-\mu) - \log\det(\Sigma^{-1}) + p\log(2\pi) \right)/2$$

This differs from standard formulas only in replacing $\log\det \Sigma$ by $-\log\det(\Sigma^{-1}).$

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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