How to deal with complicated Gaussian integrals in Mathematica?

As we know, for most Gaussian integrals, we can get the analytical result. Now I have many Gaussian integrals to treat, which have the following general form,

Integrate[
x1^n1 x2^n2 x3^n3 ... xd^
nd Exp[-{x1, x2, x3 ... xd}.matrix.{x1, x2, x3 ... xd} + c1 x1 +
c2 x2 ... cd xd], {x1, -Infinity, Infinity},
{x2, -Infinity, Infinity} ... {xd, -Infinity, Infinity}].

The matrix is a symmetrical positive definite matrix.

We can get the analytical integral result for this general form. However, it is a difficult task for Mathematica to deal with the similiar expressions, which takes so much time. I have written a package for treating this problem.

But, it is still not sufficient because I need to deal with lots of these kinds of integral expressions. So the computation time is important.

So have you ever seen the package for doing this kind of Gaussian integral? Or can you give me some suggestions about this problem? Thank you very much!

Thanks all of you. For the concept MultinormalDistribution, it is very helpful for me.

By the way, How can I transform the following integral expression to the standard form of MultinormalDistribution in a simple way?

vec1 := {r1x, r1y, r1z};
vec2 := {r2x, r2y, r2z};
vect := {r1x, r1y, r1z, r2x, r2y, r2z};
vcof := {c1, c2, c3, c4, c5, c6};
Integrate[Exp[-vect.mat.vect + vcof.vect] (vec1.vec2)^n (vec1.vec1)^
l (vec2.vec2)^m, {r1x, -Infinity, Infinity}, {r1y, -Infinity,
Infinity}, {r1z, -Infinity, Infinity}, {r2x, -Infinity,
Infinity}, {r2y, -Infinity, Infinity}, {r2z, -Infinity, Infinity}];

The mat is a symmetrical positive definite sparse matrix. And n, l, m are integers.
Thanks!

• Your SPD matrix is dense or sparse? Jun 15 '12 at 6:11
• Just like this one {{3, 0, 0, 1, 0, 0}, {0, 4, 0, 0, 1, 0}, {0, 0, 4, 0, 0, 1}, {2, 0, 0, 8, 0, 0}, {0, 2, 0, 0, 6, 0}, {0, 0, 2, 0, 0, 6}} // MatrixForm Maybe it can be called sparse matrix. Jun 16 '12 at 1:43
• Yes, it's sparse. You might be interested to find that Mathematica supports SparseArray[]. Jun 16 '12 at 1:49
• I only noticed just now that you added a follow-up question. See my updated answer, I hope that will be enough for you to put everything together.
– Jens
Jun 18 '12 at 1:45

What you have is a MultinormalDistribution. The quadratic and linear forms in the exponential can be rewritten in terms of $-\frac12(\vec{x}-\vec{\mu})^\top\Sigma^{-1}(\vec{x}-\vec{\mu})$ where $\vec{\mu}$ represents the mean and $\Sigma$ the covariance matrix, see the documentation.

With this, you can do integrals of the type given in the question by invoking Expectation, as in this example:

Expectation[
x^2 y^3, {x, y} \[Distributed]
MultinormalDistribution[{μ1, μ2},
{{σ1^2, ρ σ1 σ2},
{ρ σ1 σ2, σ2^2}}]]

The result is:

$\text{$\mu $1}^2 \text{$\mu $2}^3+\text{$\mu $2}^3 \text{$\sigma $1}^2+6 \text{$\mu $1} \text{$\mu $2}^2 \rho \text{$\sigma $1} \text{$\sigma $2}+3 \text{$\mu $1}^2 \text{$\mu $2} \text{$\sigma $2}^2+3 \text{$\mu $2} \text{$\sigma $1}^2 \text{$\sigma $2}^2+6 \text{$\mu $2} \rho ^2 \text{$\sigma $1}^2 \text{$\sigma $2}^2+6 \text{$\mu $1} \rho \text{$\sigma $1} \text{$\sigma $2}^3$

Edit

Regarding the normalization prefactor mentioned in Sjoerd's comment, you can use the fact that for any dimension $n$

$\iint\exp(-\frac{1}{2}\vec{z}^\top \Sigma^{-1}\vec{z})\mathrm dz^n = (2\pi)^{n/2}\sqrt{\det(\Sigma)}$

Hopefully these hints will be enough for you to fill in the missing linear-algebra steps to make the connection to your given matrix matrix.

Edit 2

In response to the comment by chris, for polynomials as prefactors one can also use the slightly simpler but equivalent form

Moment[
MultinormalDistribution[{μ1, μ2},
{{σ1^2, ρ σ1 σ2},
{ρ σ1 σ2, σ2^2}}],
{2, 3}]

This is the same example as above, with the powers of x and y appearing in the second argument. See the documentation for Moment.

The difference between Moment and Expectation is that Moment is restricted to the expectation values of polynomials.

Edit 3

Before going on with the symbolic manipulations that I assumed are desired here, let me also point out that you can do your integrals pretty straightforwardly if your integrand contains no symbolic parameters. Then you just need to do a numerical integral by replacing Integrate with NIntegrate.

But now back to the symbolic part:

A follow-up question arose how to complete the square in the exponential to get to the standard form of the multinormal distribution, starting from a form like this:

$$\exp(\,\vec{x}^\top A\vec{x}+\vec{v}^\top\vec{x})$$

The matrix $A$ in the exponential is symmetric, $A^\top=A$, and positive definite. Therefore $A$ is invertible, and the inverse is symmetric,

$$\left(A^{-1}\right)^\top=A^{-1}$$

With this, you can verify

$$\left(\vec{x}+\frac{1}{2}A^{-1}\vec{v}\right)^\top A\left(\vec{x}+\frac{1}{2}A^{-1}\vec{v}\right)=\vec{x}^\top A\vec{x}+\vec{v}^\top\vec{x}+\frac{1}{4}\vec{v}^\top A^{-1}\vec{v}$$

by directly multiplying out the factors on the left. Therefore,

$$\exp(\vec{x}^\top A\vec{x}+\vec{v}^\top \vec{x})=\exp(\left(\vec{x}-\vec{\mu}\right)^\top A\left(\vec{x}-\vec{\mu}\right)-\frac{1}{4}\vec{v}^\top A^{-1}\vec{v})$$

where

$$\vec{\mu}\equiv-\frac{1}{2}A^{-1}\vec{v}$$

Compare this to the standard form of the Gaussian integral, and you see that in the notation of Mathematica's documentation

$$A \equiv -\frac{1}{2} \Sigma^{-1}$$

and our integral differs from the standard Gaussian one by a factor

$$\exp(-\frac{1}{4}\vec{v}^\top A^{-1}\vec{v})$$

Now we have all the pieces that are needed, except that you still have to calculate the inverse matrix $A^{-1} \equiv -2\Sigma$, using

Inverse[mat]

if I go back to your original notation where the matrix $A$ is called mat.

Edit 4:

In view of the other answers, I put together the above steps in a module so that my approach can be compared more easily to the alternatives. The result is quite compact and is not significantly slower than the fastest alternative (by @ybeltukov):

gaussMoment[fPre_, fExp_, vars_] :=
Module[{coeff, dist, ai, μ, norm},
coeff = CoefficientArrays[fExp, vars, "Symmetric" -> True];
ai = Inverse[2 coeff[]];
μ = -ai.coeff[];
dist = MultinormalDistribution[μ, -ai];
norm = 1/PDF[dist, vars] /. Thread[vars -> μ];
Simplify[
norm Exp[1/2 coeff[].μ + coeff[]] Distribute@
Expectation[fPre, vars \[Distributed] dist]]]

The normalization factor can be easily obtained from the PDF. I used the same approach as @ybeltukov to extract the matrix $A$ from the exponent, except that I added a factor of $2$ at that stage to prevent that factor from popping up twice at later points.

Here are some tests:

RepeatedTiming[
gaussMoment[(x^2 + x^4 + x^6), -(x - 1)^2, {x}]]

$$\left\{0.0023,\frac{223 \sqrt{\pi }}{8}\right\}$$

RepeatedTiming[
gaussMoment[(x^2 + x y) , -(x - a1)^2 - (y - a2)^2 - (x - y)^2, {x, y}]]

$$\left\{0.0014,\frac{\pi \left(4 \text{a1}^2+6 \text{a1} \text{a2}+2 \text{a2}^2+3\right) e^{-\frac{1}{3} (\text{a1}-\text{a2})^2}}{6 \sqrt{3}}\right\}$$

This is about four orders of magnitude faster than doing the integrals using plain Integrate.

One other big advantage (in addition to its simplicity and speed) is that it can deal with non-polynomial prefactors. Here is an example (it takes longer to run, but the other methods cannot do it at all):

gaussMoment[
Sin[2 Pi x] Cos[Pi x] , -(x - a1)^2 - (y - a2)^2 - (x -
y)^2, {x, y}]

$$\frac{i \pi e^{-\text{a1}^2-\frac{\text{a2}^2}{2}} \left(e^{\frac{1}{6} (2 \text{a1}+\text{a2}-i \pi )^2}-e^{\frac{1}{6} (2 \text{a1}+\text{a2}+i \pi )^2}+e^{\frac{1}{6} (2 \text{a1}+\text{a2}-3 i \pi )^2}-e^{\frac{1}{6} (2 \text{a1}+\text{a2}+3 i \pi )^2}\right)}{4 \sqrt{3}}$$

• Looks OK, but isn't there a need for normalization of the original integral? Jun 15 '12 at 8:28
• I would think this is in fact the n1,nd derivative w.r.t. I c1,I c2.. cn of the characteristic function of a MultinormalDistribution; so mathematica can provide you with the Characteristic function and you can differentiate it. Jun 15 '12 at 20:42
• @chris Absolutely - and you can go even further and use the built-in command Moment for this, too. But Expectation is a bit more general in that it allows arbitrary prefactors, not just polynomials.
– Jens
Jun 15 '12 at 20:49
• For the part Edit 3, I got it. Thanks, Jens! Jun 18 '12 at 2:16
• I think it's worth mentioning that the slowness of Moment that would potentially affect this answer appears to have been fixed in Mathematica version 10.3. See also How to efficiently find moments of a multinormal distribution
– Jens
Jan 30 '16 at 1:34

If there is no pre-exponential factor then integral has a simple answer $$\int \exp({\bf x}^TA{\bf x}+{\bf b}^T{\bf x}+c)d{\bf x} = \frac{\pi^{n/2}}{\sqrt{\det(-A)}}\exp\left(c - \frac{1}{4}{\bf b}^T A^{-1}{\bf b}\right)$$

it can be implemented in Mathematica very compactly

gaussianIntegralExp[expr_, vars_] := π^(Length[vars]/2)
Exp[#1 - LinearSolve[#3, #2].#2/4]/Sqrt[Det[-#3]] & @@
CoefficientArrays[expr, vars, "Symmetric" -> True] // Simplify

gaussianIntegralExp[-(x - 1)^2 - y^2 - x y, {x, y}] // AbsoluteTiming
{0.000685, (2 E^(1/3) π)/Sqrt}
Integrate[Exp[-(x - 1)^2 - y^2 - x y], {x, -∞, ∞}, {y, -∞, ∞}] // AbsoluteTiming
{0.685393, (2 E^(1/3) π)/Sqrt}

Great speedup!

Pre-exponential factor can be treated as moment of multinormal distribution. One can use Moment and MultinormalDistribution as in the Jens' answer but it is slow (especially for large moments). I describe this problem as a separate question and post my own answer in which I describe the moment function

gaussianIntegral[expr_, vars_List] :=
Module[{coeff, exp, μ, a, b, c, rul, factor, σ, moment},
{coeff, exp} = {expr/Exp[#], #} &[Plus @@ Cases[FactorList[expr], {E^x_, p_} :> p x]];
{c, b, a} = CoefficientArrays[exp, vars, "Symmetric" -> True];
μ = LinearSolve[a, b]/2;
rul = Cases[#, x_ /; EvenQ@Total@First[x]] &@
CoefficientRules[coeff /. Thread[vars -> vars - μ], vars];
If[rul == {}, Return];
factor = π^(Length[vars]/2) Exp[c - μ.b/2]/Sqrt@Det[-a];
If[MatchQ[rul, {{0 ..} -> _}], Return[rul[[1, 2]] factor]];
σ = -Inverse[a]/2;
moment[{0 ...}] = 1;
moment[x_List] /; Total[x] == 2 := If[Length[#] == 1, σ[[#[], #[]]],
σ[[#[], #[]]]] &@Position[x, y_?Positive][[All, 1]];
moment[x_List] := moment[x] =
Module[{i = LengthWhile[x, # == 0 &] + 1, x2 = x}, x2[[i]]--;
Sum[If[x2[[j]] > 0, x2[[j]] σ[[i, j]] moment@MapAt[# - 1 &, x2, j], 0],
{j, i, Length[x]}]];
Total[#2 moment[#1] & @@@ rul] factor
]

Description:

1. Extract the pre-factor and the exponent.
2. Shift pre-factor by μ (the mean of the distribution)
3. If there is only odd moments then return 0.
4. If the pre-factor is a constant then return known result.
5. Calculate moments (see here for the description of moment).

For the nice syntax highlighting add

SyntaxInformation[gaussianIntegral] = {"ArgumentsPattern" -> {_, _},
"LocalVariables" -> {"Solve", 2}};

Examples:

gaussianIntegral[(x^2 + x^4 + x^6) Exp[-(x - 1)^2], {x}] // AbsoluteTiming
{0.003386, (223 Sqrt[π])/8}
Integrate[(x^2 + x^4 + x^6) Exp[-(x - 1)^2], {x, -∞, ∞}] // AbsoluteTiming
{0.389498, (223 Sqrt[π])/8}
gaussianIntegral[(x^2 +x y) Exp[-(x - a1)^2 - (y - a2)^2 - (x - y)^2], {x, y}]
// AbsoluteTiming // FullSimplify
{0.007893, ((3 + 2 (a1 + a2) (2 a1 + a2)) E^(-(1/3) (a1 - a2)^2) π)/(6 Sqrt)}
Integrate[(x^2 + x y) Exp[-(x - a1)^2 - (y - a2)^2 - (x - y)^2], {x, -∞, ∞}, {y, -∞, ∞}]
// AbsoluteTiming
{8.053473, ((3 + 2 (a1 + a2) (2 a1 + a2)) E^(-(1/3) (a1 - a2)^2) π)/(6 Sqrt)}

I'm a little bit late to this party, but I had written this function for another question that turned out not to need it, so I'll put this here.

My strategy is to straightforwardly calculate the integral via this formula

\begin{align} \int f(\vec x) &\, \exp\left( - \frac 1 2 \sum_{i,j=1}^{n}A_{ij} x_i x_j \right) d^nx = \\ & \sqrt{(2\pi)^n\over \det A} \, \left. \exp\left({1\over 2}\sum_{i,j=1}^{n}(A^{-1})_{ij}{\partial \over \partial x_i}{\partial \over \partial x_j}\right)f(\vec{x})\right|_{\vec{x}=0} \end{align}

I had some help from @Jens and @DanielLichtblau in writing a function for the exponential differential operator. Here is the function for the multidimensional Gaussian integral $f e^g$

gaussianIntegral[f_, g_, vars_] :=
Module[{expD, fdegree, nvars, amatrix, ainverse, jvector, rest, f2,
fnew},
fdegree =
Exponent[# /. ((# -> \[FormalX] RandomReal[]) & /@
Variables[#]), \[FormalX]] &@f;
nvars = Length@vars;
{rest, jvector,
amatrix} = ({1, 1, -2} CoefficientArrays[g, vars,
"Symmetric" -> True]);
Which[
Exponent[# /. ((# -> \[FormalX] RandomReal[]) & /@
Variables[#]), \[FormalX]] &@g != 2,
"Non Gaussian integral",
Det[amatrix] <= 0,
"A matrix is not positive definite, divergent integral",
True,
ainverse = Inverse[amatrix];
fnew = f /. Inner[Rule, vars, vars + jvector.ainverse, List];
expD = Module[{dop},
dop =
Function[{f2},
1/2 Sum[ainverse[[i, j]] D[
f2, {vars[[i]], 1}, {vars[[j]], 1}], {i, nvars}, {j,
nvars}]];
(1/Range[0, fdegree]!).NestList[dop, fnew, fdegree] /.
];
Sqrt[(2 π)^nvars/Det[amatrix]]
expD Exp[1/2 jvector.ainverse.jvector + rest]
]
];

This method is a fair bit slower than @ybeltukov's, (below I try a degree 6 polynomial and the timing is similar, but for a degree 30 polynomial my function takes ~10 times longer.

f1 = With[{n = 2},
Expand[x^Range[0, n - 1].RandomInteger[10, n]*
y^Range[0, n - 1].RandomInteger[10, n]]*
z^Range[0, n - 1].RandomInteger[10, n]]
(* (6 + 42 x + 6 y + 42 x y) (4 + 4 z) *)

With[{f = f1, g = x - x^2 - 7 y - x y - y^2 - 7 z^2 + 2 z y},
{gaussianIntegral[f, g, {x, y, z}] // AbsoluteTiming,
gaussianIntegralYbeltukov[f Exp[g], {x, y, z}] // AbsoluteTiming}]
(* {{0.001136, -((2603616 E^(398/17) π^(3/2))/(
4913 Sqrt))}, {0.004434, -((2603616 E^(398/17) π^(3/2))/(
4913 Sqrt))}} *)

You see they get the same result, however @ybeltukov's function has a couple points that mine addresses. Firstly, if you try to integrate a function with a covariance matrix that isn't positive definite, an error is thrown,

With[{f = x + y, g = 3*x - x^2 - 7*y - x*y - y^2 - 3*x*y}, {gaussianIntegral[f, g, {x, y}],
gaussianIntegral2[f*Exp[g], {x, y}], Integrate[f*Exp[g], {x, -Infinity, Infinity}, {y, -Infinity, Infinity}]}]
(* {"Covariance matrix is not positive definite, divergent integral",
(2*I*Pi)/(3*Sqrt*E^(71/6)),
Integrate[(-(1/2))*E^(49/4 + x*(17 + 3*x))*Sqrt[Pi]*(7 + 2*x)
, {x, -Infinity, Infinity}]} *)

Secondly, I get an error for ybeltukov's function when approximate numbers are used,

With[{f = 2.0 x - y^2, g = x^2 + y^2},
{gaussianIntegral[f, g, {x, y}] // AbsoluteTiming,
gaussianIntegral2[f Exp[g], {x, y}] // AbsoluteTiming}]

During evaluation of Power::infy: Infinite expression 1/0 encountered. >>

(* {{0.000606, 1.5708}, {0.031715, ComplexInfinity}} *)

Maybe someone will find this function useful