Isserlis's theorem

{{about||the result in quantum field theory about products of creation and annihilation operators|Wick's theorem}}

In probability theory, Isserlis's theorem or Wick's probability theorem is a formula that allows one to compute higher-order moments of the multivariate normal distribution in terms of its covariance matrix. It is named after Leon Isserlis.

This theorem is also particularly important in particle physics, where it is known as Wick's theorem after the work of {{harvtxt|Wick|1950}}.{{cite journal

| last = Wick

| first = G.C.

| year = 1950

| title = The evaluation of the collision matrix

| journal = Physical Review

| volume = 80

| issue = 2

| pages = 268–272

| doi = 10.1103/PhysRev.80.268

| bibcode = 1950PhRv...80..268W

}} Other applications include the analysis of portfolio returns,{{cite journal

| last1 = Repetowicz

| first1 = Przemysław

| last2 = Richmond

| first2 = Peter

| year = 2005

| title = Statistical inference of multivariate distribution parameters for non-Gaussian distributed time series

| journal = Acta Physica Polonica B

| volume = 36

| issue = 9

| pages = 2785–2796

| url = http://th-www.if.uj.edu.pl/acta/vol36/pdf/v36p2785.pdf

| bibcode = 2005AcPPB..36.2785R

}} quantum field theory{{cite journal

| last1 = Perez-Martin

| first1 = S.

| last2 = Robledo

| first2 = L.M.

| year = 2007

| title = Generalized Wick's theorem for multiquasiparticle overlaps as a limit of Gaudin's theorem

| journal = Physical Review C

| volume = 76

| issue = 6

| pages = 064314

| doi = 10.1103/PhysRevC.76.064314

| arxiv = 0707.3365

| bibcode = 2007PhRvC..76f4314P

| s2cid = 119627477

}} and generation of colored noise.{{cite journal

| last = Bartosch

| first = L.

| year = 2001

| title = Generation of colored noise

| journal = International Journal of Modern Physics C

| volume = 12

| issue = 6

| pages = 851–855

| doi = 10.1142/S0129183101002012

| bibcode = 2001IJMPC..12..851B

| s2cid = 54500670

}}

Statement

If $(X_1,\dots, X_{n})$ is a zero-mean multivariate normal random vector, then $\operatorname{E} [\,X_1 X_2\cdots X_{n}\,] = \sum_{p\in P_n^2}\prod_{\{i,j\}\in p} \operatorname{E}[\,X_i X_j\,] = \sum_{p\in P_n^2}\prod_{\{i,j\}\in p} \operatorname{Cov}(\,X_i, X_j\,),$ where the sum is over all the pairings of $\{1,\ldots,n\}$ , i.e. all distinct ways of partitioning $\{1,\ldots,n\}$ into pairs $\{i,j\}$ , and the product is over the pairs contained in $p$ .{{Cite book|url=https://www.cambridge.org/core/books/gaussian-hilbert-spaces/658C87D5A0E7FB440FC34D82B08167FC|title=Gaussian Hilbert Spaces|last=Janson|first=Svante|date=June 1997|publisher=Cambridge Core|access-date=2019-11-30|doi=10.1017/CBO9780511526169|isbn=9780521561280}}{{cite journal

| last1 = Michalowicz

| first1 = J.V.

| last2 = Nichols

| first2 = J.M.

| last3 = Bucholtz

| first3 = F.

| last4 = Olson

| first4 = C.C.

| title = An Isserlis' theorem for mixed Gaussian variables: application to the auto-bispectral density

| journal = Journal of Statistical Physics

| year = 2009

| volume = 136

| issue = 1

| pages = 89–102

| doi = 10.1007/s10955-009-9768-3

| bibcode = 2009JSP...136...89M

| s2cid = 119702133

}}

More generally, if $(Z_1,\dots, Z_{n})$ is a zero-mean complex-valued multivariate normal random vector, then the formula still holds.

The expression on the right-hand side is also known as the hafnian of the covariance matrix of $(X_1,\dots, X_{n})$ .

= Odd case =

If $n=2m+1$ is odd, there does not exist any pairing of $\{1,\ldots,2m+1\}$ . Under this hypothesis, Isserlis's theorem implies that $\operatorname{E}[\,X_1 X_2\cdots X_{2m+1}\,] = 0.$

This also follows from the fact that $-X=(-X_1,\dots,-X_n)$ has the same distribution as $X$ , which implies that $\operatorname{E}[\,X_1 \cdots X_{2m+1}\,]=\operatorname{E}[\,(-X_1) \cdots (-X_{2m+1})\,]=-\operatorname{E}[\,X_1 \cdots X_{2m+1}\,] = 0$ .

= Even case =

In his original paper,{{cite journal |last=Isserlis |first=L. |year=1918 |title=On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables |url=https://zenodo.org/record/1431593 |journal=Biometrika |volume=12 |issue=1–2 |pages=134–139 |doi=10.1093/biomet/12.1-2.134 |jstor=2331932 |authorlink=Leon Isserlis}} Leon Isserlis proves this theorem by mathematical induction, generalizing the formula for the $4^{\text{th}}$ order moments,{{cite journal |last=Isserlis |first=L. |year=1916 |title=On Certain Probable Errors and Correlation Coefficients of Multiple Frequency Distributions with Skew Regression |url=https://zenodo.org/record/1431585 |journal=Biometrika |volume=11 |issue=3 |pages=185–190 |doi=10.1093/biomet/11.3.185 |jstor=2331846 |authorlink=Leon Isserlis}} which takes the appearance

: $\operatorname{E}[\,X_1 X_2 X_3 X_4\,] =
\operatorname{E}[X_1X_2]\,\operatorname{E}[X_3X_4] +
\operatorname{E}[X_1X_3]\,\operatorname{E}[X_2X_4] +
\operatorname{E}[X_1X_4]\,\operatorname{E}[X_2X_3].$

If $n=2m$ is even, there exist $(2m)!/(2^{m}m!) = (2m-1)!!$ (see double factorial) pair partitions of $\{1,\ldots,2m\}$ : this yields $(2m)!/(2^{m}m!) = (2m-1)!!$ terms in the sum. For example, for $4^{\text{th}}$ order moments (i.e. $4$ random variables) there are three terms. For $6^{\text{th}}$ -order moments there are $3\times 5=15$ terms, and for $8^{\text{th}}$ -order moments there are $3\times5\times7 = 105$ terms.

= Example =

We can evaluate the characteristic function of gaussians by the Isserlis theorem: $E[e^{-iX}] = \sum_k \frac{(-i)^k}{k!} E[X^k] = \sum_k \frac{(-i)^{2k}}{(2k)!} E[X^{2k}] = \sum_k \frac{(-i)^{2k}}{(2k)!} \frac{(2k)!}{k!2^k} E[X^{2}]^k = e^{-\frac 12 E[X^2]}$

Proof

Since both sides of the formula are multilinear in $X_1, ..., X_n$ , if we can prove the real case, we get the complex case for free.

Let $\Sigma_{ij} = \operatorname{Cov}(X_i, X_j)$ be the covariance matrix, so that we have the zero-mean multivariate normal random vector $(X_1, ..., X_n) \sim N(0, \Sigma)$ . Since both sides of the formula are continuous with respect to $\Sigma$ , it suffices to prove the case when $\Sigma$ is invertible.

Using quadratic factorization $-x^T\Sigma^{-1}x/2 + v^Tx - v^T\Sigma v/2 = -(x-\Sigma v)^T\Sigma^{-1}(x-\Sigma v)/2$ , we get

$\frac{1}{\sqrt{(2\pi)^n\det\Sigma}}\int e^{-x^T\Sigma^{-1}x/2 + v^Tx} dx = e^{v^T\Sigma v/2}$

Differentiate under the integral sign with $\partial_{v_1, ..., v_n}|_{v_1, ..., v_n=0}$ to obtain

$E[X_1\cdots X_n] = \partial_{v_1, ..., v_n}|_{v_1, ..., v_n=0}e^{v^T\Sigma v/2}$ .

That is, we need only find the coefficient of term $v_1\cdots v_n$ in the Taylor expansion of $e^{v^T\Sigma v/2}$ .

If $n$ is odd, this is zero. So let $n = 2m$ , then we need only find the coefficient of term $v_1\cdots v_n$ in the polynomial $\frac{1}{m!}(v^T\Sigma v/2)^m$ .

Expand the polynomial and count, we obtain the formula. $\square$

Generalizations

= Gaussian integration by parts =

An equivalent formulation of the Wick's probability formula is the Gaussian integration by parts. If $(X_1,\dots X_{n})$ is a zero-mean multivariate normal random vector, then

$\operatorname{E}(X_1 f(X_1,\ldots,X_n))=\sum_{i=1}^{n} \operatorname{Cov}(X_1, X_i)\operatorname{E}(\partial_{X_i}f(X_1,\ldots,X_n)).$ This is a generalization of Stein's lemma.

The Wick's probability formula can be recovered by induction, considering the function $f:\mathbb{R}^n\to\mathbb{R}$ defined by $f(x_1,\ldots,x_n)=x_2\ldots x_n$ . Among other things, this formulation is important in Liouville conformal field theory to obtain conformal Ward identities, BPZ equations{{Cite journal|last1=Kupiainen|first1=Antti|last2=Rhodes|first2=Rémi|last3=Vargas|first3=Vincent|date=2019-11-01|title=Local Conformal Structure of Liouville Quantum Gravity|journal=Communications in Mathematical Physics|volume=371|issue=3|pages=1005–1069|doi=10.1007/s00220-018-3260-3|issn=1432-0916|bibcode=2019CMaPh.371.1005K|arxiv=1512.01802|s2cid=55282482}} and to prove the Fyodorov-Bouchaud formula.{{cite journal|last=Remy|first=Guillaume|title=The Fyodorov–Bouchaud formula and Liouville conformal field theory|journal=Duke Mathematical Journal|year=2020|volume=169|doi=10.1215/00127094-2019-0045|arxiv=1710.06897|s2cid=54777103}}

= Non-Gaussian random variables =

For non-Gaussian random variables, the moment-cumulants formula{{cite journal |last1=Leonov |first1=V. P. |last2=Shiryaev |first2=A. N. |title=On a Method of Calculation of Semi-Invariants |journal=Theory of Probability & Its Applications |date=January 1959 |volume=4 |issue=3 |pages=319–329 |doi=10.1137/1104031 |url=https://epubs.siam.org/doi/abs/10.1137/1104031|url-access=subscription }} replaces the Wick's probability formula. If $(X_1,\dots X_{n})$ is a vector of random variables, then $\operatorname{E}(X_1 \ldots X_n)=\sum_{p\in P_n} \prod_{b\in p} \kappa\big((X_i)_{i\in b}\big),$ where the sum is over all the partitions of $\{1,\ldots,n\}$ , the product is over the blocks of $p$ and $\kappa\big((X_i)_{i\in b}\big)$ is the joint cumulant of $(X_i)_{i\in b}$ .

= Uniform distribution on the unit sphere =

Consider $X = (X_1,\dots,X_d)$ uniformly distributed on the unit sphere $S^{d-1}$ , so that $\|X\|=1$ almost surely. In this setting, the following holds.

If $n$ is odd,

$\operatorname{E}\bigl[X_{i_1}\,X_{i_2}\,\cdots\,X_{i_n}\bigr] \;=\; 0.\!$

If $n = 2k$ is even,

$\operatorname{E}\bigl[X_{i_1}\,\cdots\,X_{i_{2k}}\bigr]
\;=\;
\frac{1}{d\,\bigl(d+2\bigr)\bigl(d+4\bigr)\cdots\bigl(d + 2k - 2\bigr)}
\sum_{p \in P_{2k}^2}
\prod_{\{r,s\}\in p} \delta_{\,i_r,i_s},$

where $P_{2k}^2$ is the set of all pairings of $\{1,\ldots,2k\}$ , $\delta_{i,j}$ is the Kronecker delta.

Since there are $|P_{2k}^2|=(2k - 1)!!$ delta-terms, we get on the diagonal:

$\operatorname{E}[\,X_1^{2k}\,]
\;=\;
\frac{(2k - 1)!!}{d\,\bigl(d+2\bigr)\bigl(d+4\bigr)\cdots\bigl(d + 2k - 2\bigr)}.$

Here, $(2k - 1)!!$ denotes the double factorial.

These results are discussed in the context of random vectors and irreducible representations in the work by Kushkuley (2021).{{cite arXiv | last=Kushkuley | first=Alexander | title=A Remark on Random Vectors and Irreducible Representations | eprint=2110.15504 | year=2021| class=math.PR }}