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

Question

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?

Answer 1

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.

Answer 2

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.