# 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? –  Guess who it is. 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. –  simpler Jun 16 '12 at 1:43
Yes, it's sparse. You might be interested to find that Mathematica supports SparseArray[]. –  Guess who it is. Jun 16 '12 at 1:49
Thanks, I will learn it. –  simpler Jun 16 '12 at 1:54
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 \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 \Sigma$, using

Inverse[mat]


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

-
Looks OK, but isn't there a need for normalization of the original integral? –  Sjoerd C. de Vries 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. –  chris 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! –  simpler Jun 18 '12 at 2:16

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[3]}

Integrate[Exp[-(x - 1)^2 - y^2 - x y], {x, -∞, ∞}, {y, -∞, ∞}] // AbsoluteTiming

{0.685393, (2 E^(1/3) π)/Sqrt[3]}


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[0]];
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, σ[[#[[1]], #[[1]]]],
σ[[#[[1]], #[[2]]]]] &@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[3])}

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[3])}

-