Why is the following integration involving multivariate Gaussian distribution so slow and generating an error? Is there a better integration strategy? All that I'm doing is considering $(X_1,X_2,X_3,X_4) \sim N(\mu, \Sigma)$ and computing the (trivial expectation) $E[1]$, which obviously equals to $1$.

muvec = ConstantArray[0.1, 4];
sigmat = IdentityMatrix[4];
npdf[x_] := PDF[MultinormalDistribution[muvec, sigmat], x]; 
NIntegrate[ 1 * npdf[{x1, x2, x3, x4}], 
           {x1, -Infinity, Infinity}, 
           {x2, -Infinity, Infinity}, 
           {x3, -Infinity, Infinity}, 
           {x4, -Infinity, Infinity} ] // AbsoluteTiming

This results in the errors

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>

NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 1.000000567592556` and 0.000021505642234874004` for the integral and error estimates. >>

and the result

{12.521, 1.}

The result is correct, but it has taken far too long (in my opinion), and moreover, I'm a little bit puzzled by the error message (despite reading http://reference.wolfram.com/language/tutorial/NIntegrateIntegrationStrategies.html), and I've tried other integration methods and it's still relatively slow.

In the actual application I have in mind, I want to consider $h(y) := E[ g(y;X_1,X_2,X_3,X_4)] $, where $g$ is somewhat complex and I need to use the expression $h(y)$ over and over again, and hence if I can't even compute $E[1]$ quickly and accurately, I have little hopes for computing the more difficult expectation.

Edit I emphasize that a symbolic integration of the solution is not of interest. In particular, since the question at hand is a "toy example", of which the real application is computing a more complex expectation $E[g(X_1,\ldots, X_n)]$ of which the function $g$ is sufficiently complex that there's no hope of a symbolic integral solution.


3 Answers 3


Because the integrand is highly localized at the origin; e.g.,

Log[10, npdf[{6, 0, 0, 0}]] // N
(* -9.16177 *)

limiting the range of integration reduces run time by a factor of three and moderately improves accuracy.

muvec = ConstantArray[1/10, 4];
sigmat = IdentityMatrix[4];
npdf[x_] := PDF[MultinormalDistribution[muvec, sigmat], x];
NIntegrate[1*npdf[{x1, x2, x3, x4}], {x1, -6, 6}, {x2, -6, 6}, {x3, -6, 6}, 
    {x4, -6, 6}, PrecisionGoal -> 5] // AbsoluteTiming

Not nearly as good as the answer by ciao but perhaps helpful. (Note that only the error message NIntegrate::slwcon persists.)

By the way, I first tried several NIntegrate methods but with little improvement. The issue, it seems to me, is that the integrand becomes ever more localized as the dimension increases. 3D works fine, for instance. Also, if a Precision of, say, 3 is acceptable, big savings can be achieved: {0.275684, 1.00009}. Decreasing the range of integration in this case helps further, because it concentrates the integration more on the central peak: {0.241199, 1.00001}. Taking advantages of symmetries, if any, also could be helpful. Finally, a coordinate transformation to spread the central peak relative to the wings could be helpful.

  • $\begingroup$ Thanks! This does answer my original question about computing $E[1]$ using numerical integration techniques. The subsequent follow up on $E[g(\ldots)]$ is an afterthought. $\endgroup$
    – user32416
    Commented Aug 22, 2015 at 2:51

While the unknown g might or might not be amenable to this approach, it's really fast on g == 1:

  1, {x1, x2, x3, x4} \[Distributed] 
   MultinormalDistribution[muvec, sigmat]] // RepeatedTiming

(23000 times faster than NIntegrate.)

Then there's also NExpectation to try, too.

Ever since I came across Guess who it is's use of the Chebyshev polynomial proxy this morning in How to find numerically all roots of a function in a given range?, I've been itching for something to try it out on. It's really cool. It approximates a univariate function, but chebfun has a bivariate method; perhaps it is extensible. What I want to show here is how well this works on the following.

We approximate g over the interval from -10 to 10; outside that the Gaussian is so small that practically there is no contribution to the integral.

r = 10;
g = BesselJ[2, #] &;
n = 32;
cnodes = Rescale[N[Cos[Pi Range[0, n]/n], 30], {-1, 1}, {-r, r}];
cc = Sqrt[2/n] FourierDCT[g /@ cnodes, 1];
cc[[{1, -1}]] /= 2;

gcheb = Expand[cc.Table[ChebyshevT[n - 1, x1/r], {n, Length@cc}]] /. 
   c_ /; c == 0 :> 0;
   {x1, x2, x3, x4} \[Distributed] 
   MultinormalDistribution[muvec, sigmat]] // AbsoluteTiming

(*  {0.081965, 0.0988631125184331759104}  *)

There are coefficients of the form 0``23.443340842020458, and it saves time with the integration to remove them; hence the replacement c_ /; c == 0 :> 0. The use of n = 32 was inherited from @Guess; the precision of cnodes was adjusted from 20 to 30 after examining gcheb, but it made little difference in the result.

The first NIntegrate misses the mark widely. Increasing MinRecursion helps, but it takes much, much longer, even though it recognizes a convergence problem.

  BesselJ[2, x1]*npdf[{x1, x2, x3, x4}], {x1, -Infinity, 
   Infinity}, {x2, -Infinity, Infinity}, {x3, -Infinity, 
   Infinity}, {x4, -Infinity, Infinity}] // AbsoluteTiming

(*  {105.429, 0.00890915}  *)

  BesselJ[2, x1]*npdf[{x1, x2, x3, x4}], {x1, -Infinity, 
   Infinity}, {x2, -Infinity, Infinity}, {x3, -Infinity, 
   Infinity}, {x4, -Infinity, Infinity}, 
  MinRecursion -> 3] // AbsoluteTiming


(*  {1563.69, 0.0988631}  *)
  • $\begingroup$ Thanks. For $g \equiv 1$, this is indeed a great approach (I suspect the same as symbolic integration). But for general $g$, it seems like NExpectation collapses back to NIntegrate anyway. $\endgroup$
    – user32416
    Commented Aug 22, 2015 at 2:53
  • $\begingroup$ Nice to see someone excited about this like I am. :) I've seen the paper on the bivariate Chebyshev approximant; there's still a lot to digest before I can try implementing it myself (unless you beat me to it ;) ). BTW: NIntegrate[] is capable of (tensor product) Clenshaw-Curtis quadrature, so at least that one is able to avoid an ad hoc choice of n like mine. $\endgroup$ Commented Aug 22, 2015 at 5:59
  • $\begingroup$ +1 for link to most interesting answer... $\endgroup$
    – ciao
    Commented Aug 22, 2015 at 6:33
  • $\begingroup$ @Guesswhoitis. Thanks for introducing it to the site. (I had missed your earliest answer, somehow.) Students arrived on campus yesterday and I was too busy to do any real work on it. Boyd's simplified overview makes choosing $n$ seem simple, but I didn't have time to look into it. $\endgroup$
    – Michael E2
    Commented Aug 22, 2015 at 12:01
  • $\begingroup$ Thanks, @ciao. It's a fairly simple, if advanced, technique and seems quite powerful. I hope people will appreciate Guess's other answers. My first thought was, "I wonder if this is what NSolve does?" But WRI claims the breakthrough in NSolve was due to their own research.. $\endgroup$
    – Michael E2
    Commented Aug 22, 2015 at 12:07

About 30X faster integration, no errors...

Clear[x, x1, x2, x3, x4]
muvec = ConstantArray[1/10, 4];
sigmat = IdentityMatrix[4];
pdf = PDF[MultinormalDistribution[muvec, sigmat], {x1, x2, x3, x4}] //
  1*pdf, {x1, -Infinity, Infinity}, {x2, -Infinity, Infinity}, 
         {x3, -Infinity, Infinity}, {x4, -Infinity, Infinity}] // AbsoluteTiming
  • $\begingroup$ Thanks for the answer but I don't think this is getting at the spirit of the problem. In particular, you are using symbolic integration and where I'm explicitly considering numerical integration (i.e. the real problem I have in mind is not to simply compute $E[1]$ but rather $E[g(X_1,...X_n)]$ of which there's no chance for a symbolic evaluation. Hence, I'm most interested in understanding what is the best strategy to even numerically integrate $E[1]$. $\endgroup$
    – user32416
    Commented Aug 22, 2015 at 0:02
  • 1
    $\begingroup$ @user32416 Then update your question with appropriate details for readers, else it's tilting and windmills... $\endgroup$
    – ciao
    Commented Aug 22, 2015 at 0:11
  • $\begingroup$ @user32416 and what is g? $\endgroup$
    – ciao
    Commented Aug 22, 2015 at 0:19
  • $\begingroup$ It's a highly problem specific function and I don't think that specifying it here is of general interest to the community. But we can simply think of $g$ as a continuous function in all the arguments. Thanks for asking! $\endgroup$
    – user32416
    Commented Aug 22, 2015 at 0:22
  • 1
    $\begingroup$ @user, in general, without at least knowing the structure of what's multiplying the Gaussian, it's hard to give good advice. Cubatures are hard, dude. Coming up with efficient ones more so. That's why any opportunity to exploit structure is relentlessly pursued. $\endgroup$ Commented Aug 22, 2015 at 0:32

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.