# Faster computation of high-dimensional multivariate normal probabilities

My goal is to find a faster way to calculate something like

mvtnorm::pmvnorm(upper = rep(1,100))

that is, the tail probability of multivariate normal distribution with mean 0 and arbitrary covariance matrix. The upper bound is also arbitrary. pmvnorm implements three algorithms: GenzBretz (up to dimension 1000), Miwa (up to dimension 20) and TVPACK (2- and 3-dimensional problems for semi-infinite integration regions). For my case of dimension 100, I can only use GenzBretz and it is the golden standard. However, for each run of the above code, it would take 0.6 second and it is called 40 times in each of 600 iterations. Thus the whole algorithm takes a long time.

To speed up the computation, I have tried the following strategies:

1. Decrease the accuracy of GenzBretz. I think abseps (absolute error tolerance) is the key parameter, but the default seems already liberal and relaxing it does not improve the speed largely.

GenzBretz(maxpts = 25000, abseps = 0.001, releps = 0)

2. Find alternative functions. I found OpenMx::omxMnor is faster than pmvnorm in low dimensions, but it is much slower when the dimension becomes as high as 100. Another function is bayesm::ghkvec, which computes the so called GHK approximation using an importance sampling method. So the tricky part is that the result is not fixed, which causes problems in my algorithm.

3. Search the literature. I tested the approach described in this paper but found it is slower than pmvnorm.

4. Develop a customized approximation. I was suggested to use Laplace approximation and quadrature methods. But I did not find a readily available function for the multivariate case with semi-infinite bounds. It is also not promising in theory since there is no N data points.

So my question is, is there any function faster than pmvnorm when the dimension is 100, and can maintain a (fixed and) relatively accurate approximation for multivariate normal probabilities?

• You refer to 40 calls times 600 iterations. That suggests the possibility of capitalizing on some amount of precomputation to save time later or of reusing some computations. If you would like to solicit ideas along those lines, you would need to provide some information about how all those calculations might be interrelated (if at all). – whuber Sep 7 '16 at 22:53
• I second whuber's advice. In particular, what does "arbitrary covariance" really mean in your case? In addition to re-use of previous computations, you might also consider low-rank approximations of the covariance. That is, if you are already considering lowering error tolerances or using Laplace approximation, you may be able to get better results just truncating the SVD. This would depend on your covariance matrices though. – GeoMatt22 Sep 7 '16 at 23:10
• Thanks a lot @whuber for your suggestions! I have clustered data (40 clusters) and the multivariate normal probabilities are used in the likelihood function for each cluster. I use quasi-Newton method for maximum likelihood estimation. The upper bound and the covariance matrix are updated in each iteration though they might not change much after some iterations. I thought about reusing computations and will consider it more seriously. – Randel Sep 8 '16 at 14:12
• @GeoMatt22 Thanks a lot! My covariance matrix has values like 1.05 on the diagonal and 0.05 off diagonal. Your suggestion reminds that I may treat this matrix as diagonal and can use the product of univariate normal probabilities to approximate the multivariate normal probabilities. It seems to work in my case. – Randel Sep 8 '16 at 14:37