# Speeding up computation of Gaussian integral

Suppose $$x$$ comes from standard normal in $$d$$ dimensions. Are there NIntegrate tricks I can use to speed up the following Gaussian integral?

$$f(a)=E_x\left[\frac{x_1^2}{\|x\|^2}\log (1-a \|x\|^2)^2\right]$$

Summarizing tricks from answers

• convert to polar coordinates
• use symmetry to get rid of $$x_1$$ term so it's just integral over $$r$$
• use integration bounds like {var,start,singuarlity,end} to manually specify singularity at $$r=1/a$$

The following code runs into slowness/numeric issues beyond $$d=3$$, I suspect Method->Automatic is not picking the best integration algorithm.

(* Returns Gaussian PDF in variables x[1],x[2],...,x[d] *)
normalDensity[d_] := Module[{xvec, dist},
Clear[x];
xvec = Array[x, d];
dist =
If[d == 1, NormalDistribution[],
MultinormalDistribution[IdentityMatrix[d]]];
If[d == 1, PDF[dist, First@xvec], PDF[dist, xvec]]
];

(* Computes Gaussian expectation of f in d dimensions *)
gaussianExpectation[f_, d_] := (
pdf = normalDensity[d];
xvec = Array[x, d];
bounds = {#, -\[Infinity], \[Infinity]} & /@ xvec;
NIntegrate @@ Join[{f*pdf}, bounds]
);

setupProblem := (
xvec = Array[x, d];
pdf = normalDensity[d];
normSquared = Total[xvec*xvec];
f[a_] =  (x[1]^2)/normSquared Log[(1 - a normSquared)^2] +
Log[d]/d
);

d = 3
setupProblem[d];
gaussianExpectation[f[0.552], d]  (* NIntegrate: Numerical integration converging too slowly *)


This is part of root finding algorithm to solve for $$f(a)=0$$, so the integral needs to be fast enough for this to work

obj[a_?NumericQ] := gaussianExpectation[f[a], d]
obj[0.552]
FindRoot[obj[a], {a, 0.5512}]  (* takes a long time *)

• is normSquared supposed to be less than 1? Or does it range from 0 to Infinity? Aug 24, 2022 at 0:56
• I think if you change to spherical coordinates, you can reduce the integration to Block[{d = 3, a = 0.552}, NIntegrate[(2^(-d/2) E^(-(r^2/2)) r^2 (Log[d] + Log[(-1 + a r^2)^2])) / Gamma[1 + d/2], {r, 0, 1/Sqrt[a], Infinity}] ] Aug 24, 2022 at 1:32
• @MichaelE2 the 1d integral is super-fast, but a bit mysterious to me how you turned $x_1/\|x\|$ term into polar coordinates, should this approach also work for $d=4$? Aug 24, 2022 at 16:24
• The integral over ${\Bbb R}^d$ of $g(x)=f(\|x\|)\,x_i^2/\|x\|^2$ is the same for all $i$ and $\sum x_i^2/\|x\|^2=1$; therefore the integral of $g(x)$ is equal to the integral of $f(\|x\|)/d$, for any $d$. Aug 24, 2022 at 20:49
• See "User-Specified Singularities" for {r, 0, 1/Sqrt[a], Infinity}. Aug 24, 2022 at 20:54

## 5 Answers

Note: All of the results below use the initial definition of $$f(a)$$ with

$$f(a)=E_x\left[\frac{x_1^2}{\|x\|^2}\log (1-a \|x\|^2)^2\right]$$

and the additional term Log[d]/d used later in the OP's question is not included.

If the $$X_i$$ ($$i=1,2,\ldots,d$$) have independent standard normal distributions, then $$X_i^2$$ has a chisquare distribution with 1 degree of freedom. Further, the sum of independent chisquare random variables is also a chisquare with the sum of the associated degrees of freedom. Therefore, the two random variables in the model $$X_1^2$$ and $$X_1^2+X_2^2+\cdots+X_d^2$$ (normSquared) can be written just involving 2 independent chisquare random variables: $$Z_1=X_1^2$$ and $$Z_2=X_2^2+\cdots+X_d^2$$ where $$Z_1$$ has a chisquare distribution with 1 degree of freedom and and $$Z_2$$ has a chisquare distribution with d-1 degrees of freedom.

In short any for any $$d>1$$, only a 2-dimensional integral will be necessary to estimate the desired mean.

The joint distribution of $$W_1=Z_1$$ and $$W_2=Z_1+Z_2$$ (=normSquared) is given by

dist = TransformedDistribution[{z1, z1 + z2},
{z1 \[Distributed] ChiSquareDistribution[1],
z2 \[Distributed] ChiSquareDistribution[d - 1]}]


with a pdf

PDF[dist, {w1, w2}]


So the desired mean can be calculated as follows:

f[a_?NumericQ, d_?IntegerQ] :=
NIntegrate[(2^(-d/2) E^(-w2/2) (-w1 + w2)^(1/2 (-3 + d)))/
(Sqrt[π] Sqrt[w1] Gamma[1/2 (-1 + d)])*
(w1/w2) 2 Log[Abs[1 - a w2]], {w2, 0, ∞}, {w1, 0, w2}]

AbsoluteTiming[f[2, 3]]
(* {0.0780974, 0.68992} *)

AbsoluteTiming[f[2, 10]]
(* {0.128119, 0.565409} *)


Additional simplifications

The double integral above can be simplified to a single integral. The single integrand is found to be

Integrate[(2^(-d/2) E^(-w2/2) (-w1 + w2)^(1/2 (-3 + d)))/
(Sqrt[π] Sqrt[w1] Gamma[1/2 (-1 + d)]) (w1/w2) 2 Log[Abs[1 - a w2]],
{w1, 0, w2}, Assumptions -> w2 > 0 && a ∈ Reals && d > 1 && d ∈ Integers]
(* (2^(1 - d/2) E^(-w2/2) w2^(1/2 (-2 + d)) Log[Abs[1 - a w2]])/(d Gamma[d/2]) *)


So the mean can be found with

f[a_?NumericQ, d_?IntegerQ] := NIntegrate[(2^(1 - d/2) E^(-w2/2) w2^(1/2 (-2 + d)) Log[Abs[1 - a w2]])/
(d Gamma[d/2]), {w2, 0, ∞}]

AbsoluteTiming[f[2, 3]]
(* {0.0163033, 0.68992} *)

AbsoluteTiming[f[2, 10]]
(* {0.0147712, 0.565409} *)


Even more simplifications:

For even numbered values of $$d$$, there is further simplification (and possibly a single general formula for all even values of $$d$$ but I haven't worked on that). Using the above single integrand but now using Integrate instead of NIntegrate results in the following for $$d=2,4,6,8,10$$:

t = Table[{d, Integrate[(2^(1 - d/2) E^(-w2/2) w2^(1/2 (-2 + d)) Log[Abs[1 - a w2]])/(d Gamma[d/2]),
{w2, 0, ∞}, Assumptions -> a ∈ Reals] // FullSimplify},
{d, 2, 10, 2}];
TableForm[t, TableHeadings -> {None, {"d", "Expectation"}}]


Closed-form solutions are also found for $$d=3$$ but for higher odd numbered values of $$d$$ closed-form solutions seem to be only for $$a<0$$.

• Thanks, that's indeed fast enough to compute for high d. BTW, Michael's trick (doesn't seem to be documented) of using {a,0,1/a,Infinity} as Integration bounds helps with the accuracy here as well Aug 24, 2022 at 18:48
• @MichaelE2 's comment with the simplification ($f(||x||)/d$) allows for a much simpler approach than above. Essentially $||x||$ has a chisquare distribution with d degrees of freedom. (This avoids dealing with the joint distribution approach I used above.) I'll write-up a cleaner response by sometime this weekend.
– JimB
Aug 25, 2022 at 14:11
• Thanks! useful to figuring out for which $a$ the iteration $w=w-ax\langle w,x\rangle$ converges with Gaussian $x$, surprisingly non-trivial Aug 25, 2022 at 20:00
• btw, background for this problem is given here Apr 21, 2023 at 5:53

Depending on the desired precision and what "a long time" means, why not estimate the mean by taking lots of samples?

d = 3;
nsim = 1000000;
SeedRandom[12345];
xvec = RandomVariate[NormalDistribution[0, 1], {nsim, d}];
normSquared = (# . #) & /@ xvec;
r = xvec[[All, 1]]^2/normSquared;
f[a_?NumericQ] := Mean[2 r Log[Abs[1 - a normSquared]]]

AbsoluteTiming[f /. a -> 2]
(* {0.0207494, 0.687993} *)


An approximate 95% confidence interval for the expectation is given by

f[2] + {-1, 1} 1.96*StandardDeviation[2 r Log[Abs[1 - a normSquared]] /. a -> 2]/Sqrt[nsim]
(* {0.685372, 0.690614} *)


For d = 10

AbsoluteTiming[f[2]]
(* {0.023773, 0.565282} *)

• aka a Monte Carlo! My +1 Aug 24, 2022 at 13:46

This is a cleaner approach to getting closed-form solutions to $$f(a)$$ and speedier numerical solutions when closed-form solutions don't appear to exist.

Based on the proof outlined by @MichaelE2 in their comments

$$E_x\left[\frac{x_1^2}{\|x\|^2}\log (1-a \|x\|^2)^2\right]=E_x\left[\frac{x_2^2}{\|x\|^2}\log (1-a \|x\|^2)^2\right]=\cdots=E_x\left[\frac{x_d^2}{\|x\|^2}\log (1-a \|x\|^2)^2\right]$$

and

$$\sum_{i=1}^d \frac{x_i^2}{\|x\|^2} =\frac{\sum_{i=1}^d {x_i^2}}{\sum_{i=1}^d {x_i^2}}= 1$$

Therefore it follows that

$$f(a)=\frac{1}{d}\sum_{i=1}^d E_x\left[\frac{x_i^2}{\|x\|^2}\log (1-a \|x\|^2)^2\right]=\frac{1}{d}E_x\left[\sum_{i=1}^d \frac{x_i^2}{\|x\|^2}\log (1-a \|x\|^2)^2\right]=E_x\left[\log (1-a \|x\|^2)^2\right]/d$$

So we can deal with that "simpler" expectation.

Because all of the independent $$x_i$$ random variables have a standard normal distribution (mean of 0 and variance of 1), $$z=\|x\|^2=\sum_{i=1}^d x_i^2$$ has a chisquare distribution with $$d$$ degrees of freedom.

We can define the expectation function as

f[a_, d_] := Expectation[2 Log[Abs[1 - a z]]/d,
z \[Distributed] ChiSquareDistribution[d],
Assumptions -> a ∈ Reals] // FullSimplify


(Note: I stole the change of Log[(1 - a z)^2] to 2 Log[Abs[1 - a z]] from one of the answers to one of your previous questions.)

Examples:

f[a,2]
(* -E^(-(1/2)/a) ExpIntegralEi[1/(2 a)] *)

f[a,3]


[

f[a, 10]
(* (1/(1920 a^4))(2 a (1 + 2 a (5 + 4 a (9 + 50 a))) -
(1 + 8 a (1 + 6 a (1 + 4 a + 8 a^2))) E^(-(1/2)/a) ExpIntegralEi[1/(2 a)]) *)


Closed-form solutions appear to exist for even $$d$$ but for odd $$d \geq 5$$ a closed-form solution seems only to exist for $$a<0$$.

For those instances where no closed-form solution exists, the following could be used:

fn[a_, d_] := NExpectation[2 Log[Abs[1 - a z]]/d,
z \[Distributed] ChiSquareDistribution[d]] // Chop

fn[1, 5]
(* 0.39344 *)

• cool, even faster! BTW, constraints of the problem also imply that a>0 Aug 25, 2022 at 23:48

I have closed form for the integral (See answer of user: JimB "Even more simplifications:"): $$\int_0^{\infty } \frac{2^{1-\frac{d}{2}} e^{-\frac{\text{w2}}{2}} \text{w2}^{\frac{1}{2} (-2+d)} \log (| 1-a \text{w2}| )}{d \Gamma \left(\frac{d}{2}\right)} \, d\text{w2}=\\\frac{2 \pi \cot \left(\frac{d \pi }{2}\right)}{d}-\frac{\pi \cot \left(\frac{d \pi }{2}\right) \Gamma \left(\frac{d}{2},\frac{1}{2 a}\right)}{\Gamma \left(1+\frac{d}{2}\right)}-\frac{2 \, _2F_2\left(1,1;2,2-\frac{d}{2};-\frac{1}{2 a}\right)}{a (-2+d) d}+\frac{\log (4)}{d}+\frac{2 \log (a)}{d}+\frac{2 \psi ^{(0)}\left(\frac{d}{2}\right)}{d}$$ Only works for d=1,3,5,7,9,...,but for even d Mathematica give me: "Indeterminate" !!!

ANSWER = Simplify[InverseMellinTransform[
First@Integrate[
MellinTransform[(
2^(1 - d/2) E^(-w2/2) w2^(1/2 (-2 + d)) Log[Abs[1 - a w2]])/(
d Gamma[d/2]), a, s], {w2, 0, \[Infinity]},
Assumptions -> {a \[Element] Reals,
d \[Element] PositiveIntegers}], s, a,
Assumptions -> -1 < Re[s] < 0],
Assumptions -> {a \[Element] Reals,
d \[Element] PositiveIntegers}] // FunctionExpand

(*(2 \[Pi] Cot[(d \[Pi])/2])/d - (\[Pi] Cot[(d \[Pi])/2] Gamma[d/2, 1/(
2 a)])/Gamma[1 + d/2] - (
2 HypergeometricPFQ[{1, 1}, {2, 2 - d/2}, -(1/(2 a))])/(
a (-2 + d) d) + Log[4]/d + (2 Log[a])/d + (2 PolyGamma[0, d/2])/d*)

AbsoluteTiming[N[ANSWER /. a -> 2 /. d -> 3]]
(*{0.0027664, 0.68992}*)

f[a_?NumericQ, d_?IntegerQ] := NIntegrate[(2^(1 - d/2) E^(-w2/2) w2^(1/2 (-2 + d)) Log[
Abs[1 - a w2]])/(d Gamma[d/2]), {w2, 0, \[Infinity]},
Method -> "LocalAdaptive"]
AbsoluteTiming[f[2, 3]]
(*{0.0324361, 0.68992}*)


Workaround for even d:

 \$MaxExtraPrecision = 100;

{With[{e = 10^-30}, N[(ANSWER /. a -> 2 /. d -> 10 + e), 50]],
With[{e = 10^-30}, N[(ANSWER /. a -> 2 /. d -> 10 - e), 50]]}(*For a=2,d=10*)
{0.56540897717124285593865218296006923023978008859065, 0.56540897717124285593865218296013455689940774731420}


In the plot we see any discontinuity :

Plot[(ANSWER /. a -> 2), {d, 0, 20}]


• +1 Very nice! (And I think there is a closed-form for all evens. )
– JimB
Aug 26, 2022 at 15:05

Note that numeric integration gets very expensive computationally in high dimensions. As an example, if you want 100 sample points per dimension, for a 10-dimensional function you would already need $$10^{20}$$ sample points. In other words, simply sampling points to do numeric integration is in general not a good strategy for higher dimensions.

As written in Michael E2's comment, it seems you can switch to polar coordinates and then your example only depends on the radius. This should massively help to speed up computations. In general trying to reduce the dimensions of your problem is the best strategy to speed up numeric integration in high dimensions.

• Your answer would be more appropriate as a comment. Aug 24, 2022 at 12:08
• @bbgodfrey OP asked how to speed up their computation. This provides a possible answer to that question and some background why it works. To me Michael E2's comment is not appropriate as a comment because it should be an answer instead. Aug 24, 2022 at 12:12