# Efficient calculation of the multivariate normal density function

The formula for the multivariate normal density function in the standard form contains $\Sigma^{-1}$ and the determinant of $\Sigma$, which are not very computationally friendly. Is it possible to calculate the density in some other way, that is more suitable for computer implementation?

EDIT: To be precise I am after the log of the density function, if it helps.

-
why do you say $\Sigma$ is not computer friendly? – mathemagician Jan 17 '13 at 13:13
@mathemagician, calculating determinants and inverse matrices is not computer friendly in general, because they are slow and numerically unstable (that's what I have been told at least!) – Grzenio Jan 17 '13 at 13:15
It is not that bad if you use good software packages. It could all go wrong tho if the matrix has both very large and very small numbers, as a process called floating point numbers is used. I don't know any further about this, but I use R and it works fine, even with very large matrices. – mathemagician Jan 17 '13 at 13:25

http://gallery.rcpp.org/articles/dmvnorm_arma/

They use an alternative implementation of the function dmvnorm() from package mvtnorm() in R. It's a quite good time reduction, take a look at this example I just ran:

> x=matrix(0,1,40)
> m=rep(0,40)
> s=diag(40)
> dmvnorm(x,m,s)
[1] 1.087433e-16
> dMvn(x,m,s)
[1] 1.087433e-16
> benchmark(dmvnorm(x,m,s), dMvn(x,m,s), replications=10000)[,1:4]
test replications elapsed relative
2    dMvn(x, m, s)        10000    1.35    1.000
1 dmvnorm(x, m, s)        10000    5.07    3.756


I have a C implementation of dmvnorm that is taking more than 95% of the total running time. Now, I will implement this dMvn() and I hope that the bottleneck can tone down. The function is:

dMvn <- function(X,mu,Sigma) {
k <- ncol(X)
rooti <- backsolve(chol(Sigma),diag(k))