# How to speed up a high dimensional numerical integration of Gaussian-form integrand [duplicate]

My integral has the following form:

$$\int \sin(\theta) d\theta d\phi\left| \int d^3p_1 d^3p_2d^3p_3 \frac{f(p,\theta,\phi)}{g(p,\theta,\phi)}e^{-\frac{1}{2}(p^\mathrm{T} .A.p)}\right |^2,$$ where $p=(p_1,p_2,p_3)^\mathrm{T}$ and $f$ and $g$ are polynomials of $p$,and $A$ is a diagonal matrix.

Because of the complexity of the integrand, there is no analytic expression of such integration. In order to use NIntegrate, I have to change the integration to this form:

$$\int \sin(\theta)d\theta d\phi \int d^3p_1 d^3p_2d^3p_3 d^3p'_1 d^3p'_2d^3p'_3\frac{f(p,\theta,\phi)f^*(p',\theta,\phi)}{g(p,\theta,\phi)g^*(p',\theta,\phi)}e^{-\frac{1}{2}(p^\mathrm{T} .A.p)}e^{-\frac{1}{2}(p'^\mathrm{T} .A.p')},$$

The problem is the fact that the numerical integration is very slow even if the integrand has a relatively simple form:

Timing@
NIntegrate[
E^-(Qx^2 + Qy^2 + Qz^2 + Qx2^2 + Qy2^2 + Qz2^2) *
(Sin[θ]^2 - Cos[ϕ] Sin[θ] (Qx + Qx2) -
Sin[θ] Sin[ϕ] (Qy + Qy2) + (Qx Qx2 + Qy Qy2)) Sin[θ],
{Qx, -∞, ∞}, {Qy, -∞, ∞}, {Qz, -∞, ∞}, {Qx2, -∞, ∞}, {Qy2, -∞, ∞}, {Qz2, -∞, ∞},
{θ, 0, π}, {ϕ, 0, 2 π}
]

and Mathematica returns

NIntegrate::slwcon : Numerical integration converging too slowly...
NIntegrate::eincr: The global error of the strategy GlobalAdaptive...

(*{72.7813,259.935}*)

How can I speed up this integration? Are there any options for NIntegrate that might help in my case?

• This problem bears some resemblance to 92102, which may be of some help. See also the Wolfram tutorial on Advanced Numerical Integration. – bbgodfrey Aug 26 '15 at 3:20
• I posted an answer using the method described in my answer here: How to deal with complicated gaussian integrals in Mathematica?. You could also use NExpectation if numerical integration is really necessary. – Jens Sep 8 '15 at 1:40
• @bbgodfrey Thanks. I've seen that post before asking this one. I've even walking through the Advanced Numerical Integration long before. After some physical considerations, I've simplified the question to [this] one(mathematica.stackexchange.com/q/92684/6468). I think mathematical transformations have to be done before direct numerical integration. I think it is partially solved in that link. I've asked too contrieved example in this post. And I forgot to close this question after that post. Sorry about this. – luyuwuli Sep 9 '15 at 11:54

I don't consider myself an expert on numerical integration, but I can give you a few hints. One way to speed up integrations is to ask yourself which accuracy you require. You can set this with the option

PrecisionGoal -> # of digits

If you use your results to produce a plot, an accuracy of 2 digits is usually more than sufficient. Another possibility might be to rewrite the integrand into a more suitable form. For numerical evaluation of integrals it is generally preferable to avoid subtractions of large numbers (even though I believe that Mathematica takes already care of this through some preprocessing of the integrand). In the current case a possibility would be the use of spherical coordinates for the d^3Q d^3Q2 integrals.

integrand =E^-(Qx^2 + Qy^2 + Qz^2 + Qx2^2 + Qy2^2 + Qz2^2)*(Sin[θ]^2 -
Cos[ϕ] Sin[θ] (Qx + Qx2) -
Sin[θ] Sin[ϕ] (Qy + Qy2) + (Qx Qx2 +
Qy Qy2)) Sin[θ]/
R^2 /. {-Qx2^2 - Qy2^2 - Qz2^2 -> - R2^2} /. {-Qx^2 - Qy^2 -
Qz^2 -> -R1^2 } /. {Qx -> R1 Cos[theta] Sin[phi],
Qy -> R1 Sin[theta] Sin[phi], Qz -> R1 Cos[phi],
Qx2 -> R2 Cos[theta2] Sin[phi2], Qy2 -> R2 Sin[theta2] Sin[phi2],
Qz2 -> R2 Cos[phi2] } /. {R1 -> x* R, R2 -> (1 - x) *R}

with an additional jacobian factor of

jac =  R^7 x^2 (1 - x)^2 * Sin[phi]*Sin[phi2]

I used these coordinates in my following attempts, but I am actually not sure whether they give much of an improvement in this case. For NIntegrate without options it seems to make no difference.

NIntegrate[integrand * R^7 x^2 (1 - x)^2 * Sin[phi]*Sin[phi2], {R, 0,
Infinity}, {x, 0, 1}, {theta, 0, 2 Pi}, {theta2, 0, 2 Pi}, {phi,
0, Pi}, {phi2, 0, Pi}, {θ, 0, Pi}, {ϕ, 0,
2 Pi}] // AbsoluteTiming

it takes 64.389507 sec. and gives you the following warning

"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 259.727 and 0.97042 for the integral and error \ estimates"

With PrecisionGoal -> 2, the integral requires only 8.4 sec. and gives you as a result 255.435 (without any warning), which is not 259.727, but agrees within the required amount of digits. The integration can be further accelerated, if one uses MonteCarlo methods, which are generally suitable for higher dimensional integration problems. In the current case you find

NIntegrate[integrand * R^7 x^2 (1 - x)^2 * Sin[phi]*Sin[phi2], {R, 0,
Infinity}, {x, 0, 1}, {theta, 0, 2 Pi}, {theta2, 0, 2 Pi}, {phi,
0, Pi}, {phi2, 0, Pi}, {θ, 0, Pi}, {ϕ, 0, 2 Pi},
Method -> "MonteCarlo", MaxPoints -> 50000000,
PrecisionGoal -> 2] // AbsoluteTiming

yields as a result 260.314 in 2.607920 sec. Mathematica provides a few refined MonteCarlo methods, see http://reference.wolfram.com/language/tutorial/NIntegrateOverview.html, which are for sure worth to take a look at. Requesting on the other hand more than 2 digits, leads in this case to an enormous increase in computation time and/or failure to achieve the required precision. With

NIntegrate[integrand * R^7 x^2 (1 - x)^2 * Sin[phi]*Sin[phi2], {R, 0,
Infinity}, {x, 0, 1}, {theta, 0, 2 Pi}, {theta2, 0, 2 Pi}, {phi,
0, Pi}, {phi2, 0, Pi}, {θ, 0, Pi}, {ϕ, 0, 2 Pi},
MaxRecursion -> 10,
Method -> {GlobalAdaptive, MaxErrorIncreases -> 10000},
PrecisionGoal -> 4] // AbsoluteTiming

I still remain after 308.220384 sec. of integration time with an error estimate of .3305 which does not yet reach the required 4 digits. Playing around with MaxRecursion and MaxErrorIncreases might eventually lead to the required accuracy. In case of the MonteCarlo method, an increase of MaxPoints helps

NIntegrate[integrand * jac, {R, 0, Infinity}, {x, 0, 1}, {theta, 0,
2 Pi}, {theta2, 0, 2 Pi}, {phi, 0, Pi}, {phi2, 0, Pi}, {θ,
0, Pi}, {ϕ, 0, 2 Pi}, Method -> "MonteCarlo",
MaxPoints -> 100000000, PrecisionGoal -> 3] // AbsoluteTiming

provides in 266.21365 sec. the result 259.392 with the required accuracy.

• Thanks. However, in the real problems, after physical simplifications, the integrand is like this. There are hypersingular points in the integrand. And even in that post it's a simplied version, there are actually three quadratic polynomials in the denominator. the integration variables is 9. So Monte-carlo rules can't be used. Even though I haven't tried the sperical basis, I would doubt that this transformation would work in that case. – luyuwuli Sep 9 '15 at 12:08

Observe first that the Gaussian can also be written as

FullSimplify[Pi^3 PDF[
MultinormalDistribution[{0, 0, 0, 0, 0, 0},
1/2 IdentityMatrix[6]], {Qx, Qy, Qz, Qx2, Qy2, Qz2}]]

(* ==> E^(-Qx^2 - Qx2^2 - Qy^2 - Qy2^2 - Qz^2 - Qz2^2) *)

Therefore, we can calculate the integral as

Pi^3 Integrate[
Expectation[(Sin[θ]^2 -
Cos[ϕ] Sin[θ] (Qx + Qx2) -
Sin[θ] Sin[ϕ] (Qy + Qy2) + (Qx Qx2 +
Qy Qy2)) Sin[θ], {Qx, Qy, Qz, Qx2, Qy2,
Qz2} \[Distributed]
MultinormalDistribution[{0, 0, 0, 0, 0, 0},
1/2 IdentityMatrix[6]]], {θ, 0, Pi}, {ϕ, 0,
2 Pi}]

(* ==> (8 Pi^4)/3 *)

N[%]

(* ==> 259.758 *)

This calculation is almost instantaneous.

• Thanks and sorry for the late reply. However, in the actual problems, the integrand has some complicated denominators which maybe equal to zero on some hypersurface. I tried the contrived simple example and found that NIntegrate is too slow for this, let alone the complicated one. Based on some physical simplification, I've asked a related question and partially solve the problems. I saw your answer before, but NExpectation deeply is the same as NIntegrate which can be confirmed by trying the example in Edit in that post. – luyuwuli Sep 9 '15 at 11:45