---
title: 'Isserlis''s Estimator: Gaussian Moments & RMA'
url: https://www.emergentmind.com/topics/isserlis-s-estimator
type: topic
---

# Isserlis's Estimator: Gaussian Moments & RMA

“Isserlis’s estimator” denotes two distinct constructions in statistical literature. In the contemporary Gaussian moment-tensor setting, it is a plug-in estimator obtained by substituting the sample covariance into the Wick–Isserlis expansion for even-order Gaussian moments; in that role, it yields dimension-free, non-asymptotic improvements over the standard sample moment estimator for tensors of even order \(p>2\) [2507.06166]. In the classical errors-in-variables literature, the same label is often applied to the reduced major axis or standard major axis slope \(b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)\), although the invariant characterization of that functional form is attributed in the cited Bayesian treatment to Frisch (1934) and Samuelson (1942), not to Isserlis [1202.0957].

## 1. Terminological scope and dual usage

The term has a bifurcated usage. One usage is exact and theorem-driven: for Gaussian vectors, higher-order even moments are polynomial functionals of the covariance matrix \(\Sigma\), so estimation can proceed by first estimating \(\Sigma\) and then applying Isserlis’s theorem. This is the meaning developed in the study of Gaussian moment tensors [2507.06166].

The second usage belongs to invariant straight-line fitting with measurement error in both coordinates. There, the estimator often called “Isserlis’s estimator” is the reduced major axis, also called the standard major axis or geometric mean functional relationship, with slope
\[
b_{\mathrm{RMA}}=\operatorname{sign}(r)\frac{s_y}{s_x}.
\]
The same source explicitly notes that it does not cite Isserlis for this construction, and instead attributes the invariance characterization and the specific slope form to Frisch (1934) and Samuelson (1942) [1202.0957].

This indicates that the label is not univocal across subfields. In one context it names a covariance plug-in estimator justified by the Wick–Isserlis formula; in the other it refers to an interchange- and scale-invariant point estimate for a bivariate linear relationship under unspecified measurement-error variances.

## 2. Gaussian moment-tensor estimator from Wick–Isserlis expansions

Let \(X\in\mathbb{R}^d\) be Gaussian with law \(X\sim\mathcal{N}(\mu,\Sigma)\). The raw moment tensor of order \(p\) is
\[
M_p:=\mathbb{E}[X^{\otimes p}],
\]
with entries \((M_p)_{j_1,\dots,j_p}=\mathbb{E}[X_{j_1}\cdots X_{j_p}]\). The central moment tensor is
\[
\tilde M_p:=\mathbb{E}[(X-\mu)^{\otimes p}],
\]
with entries
\[
(\tilde M_p)_{j_1,\dots,j_p}=\mathbb{E}[(X_{j_1}-\mu_{j_1})\cdots(X_{j_p}-\mu_{j_p})].
\]
For zero-mean Gaussian vectors, all odd central moments vanish, \(\tilde M_p=0\) whenever \(p\) is odd, while for even \(p\) the tensor is completely determined by \(\Sigma\) through Isserlis’s theorem [2507.06166].

For centered Gaussian \(X\sim\mathcal{N}(0,\Sigma)\), Isserlis’s theorem gives, for even \(p\),
\[
\mathbb{E}[X_{j_1}\cdots X_{j_p}]
=
\sum_{\pi\in\mathcal{P}_2([p])}\prod_{(a,b)\in\pi}\Sigma_{j_a j_b},
\]
and the expectation is zero for odd \(p\). Here \(\mathcal{P}_2([p])\) is the set of all perfect matchings of \(\{1,\dots,p\}\), with cardinality \((p-1)!!\). For \(p=4\),
\[
\mathbb{E}[X_iX_jX_kX_\ell]
=
\Sigma_{ij}\Sigma_{k\ell}
+\Sigma_{ik}\Sigma_{j\ell}
+\Sigma_{i\ell}\Sigma_{jk}.
\]
Tensorially,
\[
\mathbb{E}[X^{\otimes p}]
=
\sum_{\pi\in\mathcal{P}_2([p])}\Sigma^{\sharp\pi},
\]
where \(\Sigma^{\sharp\pi}\) inserts covariance entries slotwise according to the pairing \(\pi\); the resulting tensor is symmetric in its \(p\) modes [2507.06166].

Given i.i.d. samples \(X_1,\dots,X_n\) from \(\mathcal{N}(\mu,\Sigma)\), define
\[
\hat\mu=\frac{1}{n}\sum_{t=1}^n X_t,
\qquad
\hat\Sigma=\frac{1}{n}\sum_{t=1}^n (X_t-\hat\mu)(X_t-\hat\mu)^\top.
\]
For even-order central moments, the Isserlis estimator is obtained by substituting \(\hat\Sigma\) into Wick’s formula:
\[
\big(\hat M_p^{\mathrm{Isserlis}}\big)_{j_1,\ldots,j_p}
=
\sum_{\pi\in\mathcal{P}_2([p])}\prod_{(a,b)\in\pi}\hat\Sigma_{j_a j_b},
\]
or, in tensor notation,
\[
\hat M_p^{\mathrm{Isserlis}}
=
\sum_{\pi\in\mathcal{P}_2([p])}\hat\Sigma^{\sharp\pi}.
\]
For odd \(p\), it is set to \(0\) for a centered Gaussian target. For raw moments with \(\mu\neq 0\), one reconstructs
\[
\hat M_p
=
\sum_{k=0}^p \binom{p}{k}\hat\mu^{\otimes(p-k)}\odot \hat{\tilde M}_k^{\mathrm{Isserlis}},
\]
where \(\hat{\tilde M}_k^{\mathrm{Isserlis}}\) is computed by Wick’s formula for even \(k\) and set to zero for odd \(k\) [2507.06166].

The contrast with the standard sample moment estimator is structural. The standard estimator
\[
\hat M_p^{\mathrm{SM}}=\frac{1}{n}\sum_{t=1}^n X_t^{\otimes p}
\]
targets the \(p\)-th Gaussian chaos directly, whereas the Isserlis estimator targets covariance first and then maps it to the \(p\)-th moment through Wick’s formula. The stated consequence is that the plug-in estimator removes higher-order chaos and reduces the problem to covariance estimation [2507.06166].

## 3. Norms, non-asymptotic bounds, and sample-complexity gains

The analysis is carried out in the operator and entrywise maximum norms. For an order-\(p\) tensor \(T\in\mathbb{R}^{d_1\times\cdots\times d_p}\),
\[
\|T\|_{\mathrm{op}}
:=
\sup_{v_k\in S^{d_k},\,1\le k\le p}
\langle T, v_1\otimes\cdots\otimes v_p\rangle,
\]
and
\[
\|T\|_\infty
:=
\max_{i_1,\dots,i_p}|T_{i_1,\dots,i_p}|.
\]
The bounds depend on effective dimensions rather than ambient dimension \(d\):
\[
r_2(\Sigma):=\frac{\mathrm{Tr}(\Sigma)}{\|\Sigma\|},
\qquad
r_{\max}(\Sigma):=\frac{(\mathbb{E}\|X\|_\infty)^2}{\|\Sigma\|_{\max}}.
\]
The paper describes the resulting guarantees as dimension-free because they contain no explicit dependence on the ambient dimension [2507.06166].

In the symmetric centered case \(X\sim\mathcal{N}(0,\Sigma)\) with even \(p\) and target \(T=\mathbb{E}[X^{\otimes p}]\), the expected operator-norm error satisfies
\[
\mathbb{E}\|\hat T_S-T\|_{\mathrm{op}}
\asymp_p
\|\Sigma\|^{p/2}
\left(
\sqrt{\frac{r_2(\Sigma)}{n}}
+
\frac{r_2(\Sigma)^{p/2}}{n}
\right),
\]
while
\[
\mathbb{E}\|\hat T_I-T\|_{\mathrm{op}}
\lesssim_p
\|\Sigma\|^{p/2}
\left(
\sqrt{\frac{r_2(\Sigma)}{n}}
+
\left(\frac{r_2(\Sigma)}{n}\right)^{p/2}
\right).
\]
Entrywise, the analogous statements are
\[
\mathbb{E}\|\hat T_S-T\|_{\infty}
\asymp_p
\|\Sigma\|_{\max}^{p/2}
\left(
\sqrt{\frac{r_{\max}(\Sigma)}{n}}
+
\frac{r_{\max}(\Sigma)^{p/2}}{n}
\right),
\]
and
\[
\mathbb{E}\|\hat T_I-T\|_{\infty}
\lesssim_p
\|\Sigma\|_{\max}^{p/2}
\left(
\sqrt{\frac{r_{\max}(\Sigma)}{n}}
+
\left(\frac{r_{\max}(\Sigma)}{n}\right)^{p/2}
\right).
\]

| Quantity | Standard sample moment | Isserlis estimator |
|---|---|---|
| Operator norm | \(\asymp_p \|\Sigma\|^{p/2}\left(\sqrt{r_2/n}+r_2^{p/2}/n\right)\) | \(\lesssim_p \|\Sigma\|^{p/2}\left(\sqrt{r_2/n}+(r_2/n)^{p/2}\right)\) |
| Entrywise maximum norm | \(\asymp_p \|\Sigma\|_{\max}^{p/2}\left(\sqrt{r_{\max}/n}+r_{\max}^{p/2}/n\right)\) | \(\lesssim_p \|\Sigma\|_{\max}^{p/2}\left(\sqrt{r_{\max}/n}+(r_{\max}/n)^{p/2}\right)\) |
| Consistency threshold | \(n\gg r_2(\Sigma)^{p/2}\) or \(n\gg r_{\max}(\Sigma)^{p/2}\) | \(n\gg r_2(\Sigma)\) or \(n\gg r_{\max}(\Sigma)\) |

These bounds quantify a strict improvement. The standard estimator is consistent only when \(n\gg r_2(\Sigma)^{p/2}\) in operator norm, or \(n\gg r_{\max}(\Sigma)^{p/2}\) entrywise, because the term \(r_2(\Sigma)^{p/2}/n\) or \(r_{\max}(\Sigma)^{p/2}/n\) must be negligible. By contrast, the Isserlis estimator is consistent already when \(n\gg r_2(\Sigma)\) or \(n\gg r_{\max}(\Sigma)\), inheriting the sample-covariance rates [2507.06166].

The mechanism is summarized by the perturbation control
\[
\frac{\|\hat M_p^{\mathrm{Isserlis}}-M_p\|}{\|M_p\|}
\le
\frac{p}{2}\cdot \frac{\|\hat\Sigma-\Sigma\|}{\|\Sigma\|}
\left(1+\frac{\|\hat\Sigma-\Sigma\|}{\|\Sigma\|}\right)^{p/2-1}.
\]
This expresses the moment tensor as a smooth functional of \(\Sigma\), so its estimation error is governed by covariance perturbation rather than by direct \(p\)-th-chaos fluctuations. For \(p=4\) and \(p=6\), the paper gives explicit operator- and entrywise-norm bounds of the same form, with the higher-order correction terms \((r_2(\Sigma)/n)^2\) and \((r_2(\Sigma)/n)^3\), respectively [2507.06166].

## 4. Asymmetric tensors, assumptions, and computation

The same plug-in principle extends to asymmetric Gaussian tensors. Let
\[
X=(X^{(1)},\dots,X^{(p)})
\]
be zero-mean Gaussian with block covariances
\[
\Sigma^{(j,k)}:=\mathbb{E}[X^{(j)}\otimes X^{(k)}],
\]
and target
\[
T:=\mathbb{E}[X^{(1)}\otimes\cdots\otimes X^{(p)}].
\]
In this setting the paper proves dimension-free operator- and entrywise-norm bounds in which the Isserlis estimator depends only on the largest block effective dimension, while the standard sample moment estimator carries an additional product term involving all blocks [2507.06166].

More precisely, provided \(n\ge \max_k r_2(\Sigma^{(k)})\),
\[
\mathbb{E}\|\hat T_I-T\|_{\mathrm{op}}
\lesssim_p
\left(\prod_{k=1}^p \|\Sigma^{(k)}\|^{1/2}\right)
\left(\frac{\max_k r_2(\Sigma^{(k)})}{n}\right)^{1/2},
\]
whereas the standard estimator satisfies a larger bound containing both a sum term and a product term. Entrywise, provided \(n\ge \max_k r_{\max}(\Sigma^{(k)})\),
\[
\mathbb{E}\|\hat T_I-T\|_{\infty}
\lesssim_p
\left(\prod_{k=1}^p \|\Sigma^{(k)}\|_{\max}^{1/2}\right)
\left(\frac{\max_k r_{\max}(\Sigma^{(k)})}{n}\right)^{1/2}.
\]
The stated implication is that \(\hat T_I\) is consistent once \(n\) dominates the largest effective dimension among blocks, while \(\hat T_S\) is consistent only when \(n\) dominates the product term [2507.06166].

The scope conditions are explicit. The results require i.i.d. Gaussian samples, exact Gaussianity, and even \(p\) for nontrivial central moments. In the symmetric case, \(\mu=0\) is assumed for the cleanest formulation; nonzero mean is handled by centering with \(\hat\mu\) and adding mean terms. For odd \(p\), central moments vanish for centered Gaussians, and the Isserlis estimator is identically zero [2507.06166].

The computational bottleneck is combinatorial: the number of pairings is \((p-1)!!\). The paper therefore recommends exploiting symmetry, evaluating contractions such as
\[
\langle \hat M_p^{\mathrm{Isserlis}}, v_1\otimes\cdots\otimes v_p\rangle
\]
as sums of products of quadratic forms \(v_a^\top \hat\Sigma v_b\), and using the recurrence
\[
\mathbb{E}[X_{j_1}\cdots X_{j_p}]
=
\sum_{b=2}^p
\Sigma_{j_1j_b}\,
\mathbb{E}[X_{j_2}\cdots \widehat{X_{j_b}}\cdots X_{j_p}]
\]
for dynamic programming. Exhaustive pairing enumeration is described as typically feasible for moderate \(p\), such as \(p\le 8\) or \(10\), while contraction-based evaluation is preferable for larger \(p\). The same source notes that although the paper focuses on exact Gaussianity, numerical evidence suggests the advantage may persist more broadly under near-Gaussian misspecification [2507.06166].

## 5. The reduced major axis interpretation in errors-in-variables models

In a different statistical tradition, “Isserlis’s estimator” refers to the slope used in invariant straight-line fitting when both coordinates are measured with error. The underlying structural model has observations
\[
y_{1i}=\xi_{1i}+u_{1i},
\qquad
y_{2i}=\xi_{2i}+u_{2i},
\qquad i=1,\dots,n,
\]
with latent relation \(\xi_{2i}=\alpha+\beta \xi_{1i}\). In the structural normal model, \(\xi_{1i}\sim N(\mu_1,\tau^2)\), and the total errors \(u_1,u_2\) are independent, mean-zero normal with variances \(\sigma_1^2,\sigma_2^2\). The observed pairs are then bivariate normal with covariance
\[
\Sigma=
\begin{bmatrix}
\tau^2+\sigma_1^2 & \beta\tau^2 \\
\beta\tau^2 & \beta^2\tau^2+\sigma_2^2
\end{bmatrix}
\]
[1202.0957].

If \(S\) is the sample covariance matrix with divisor \(n-1\), define
\[
r=\frac{S_{12}}{(S_{11}S_{22})^{1/2}},
\qquad
s_x=S_{11}^{1/2},
\qquad
s_y=S_{22}^{1/2},
\qquad
l=\frac{s_y}{s_x}.
\]
The reduced major axis or standard major axis slope is
\[
b_{\mathrm{RMA}}=\operatorname{sign}(r)\,l.
\]
It depends only on \(r\) and \(s_y/s_x\), is scale invariant, and under interchange of axes maps to its reciprocal \(1/b_{\mathrm{RMA}}\). The same functional form can be written as the geometric mean of the two OLS slopes,
\[
b_{\mathrm{RMA}}
=
\operatorname{sign}(r)\sqrt{b_{Y|X}\,b_{X|Y}^{-1}},
\]
where
\[
b_{Y|X}=\frac{S_{12}}{S_{11}}=rl,
\qquad
b_{X|Y}^{-1}=\frac{S_{22}}{S_{12}}=\frac{l}{r}.
\]
Samuelson’s result, as summarized in the cited source, is that among point estimates depending only on \(r\) and \(s_y/s_x\) and respecting scale and interchange invariance, the only possibility is this geometric mean form [1202.0957].

The same paper is explicit about the status of the name. Although many sources refer to the RMA/SMA or geometric mean functional relationship as “Isserlis’s estimator,” the paper does not cite Isserlis, and instead attributes the invariance characterization and specific slope form to Frisch (1934) and Samuelson (1942) [1202.0957].

## 6. Non-identifiability, Bayesian alternatives, and applied interpretation

The structural normal model does not identify \(\beta\) without additional information. Reiersøl (1950) showed that the sampling distribution identifies \(\Sigma\) but not \(\beta\), and the paper emphasizes that the RMA/SMA slope is therefore an invariant point estimate rather than a generally identified estimator of the structural slope unless further conditions hold, such as specified error-variance information or special symmetry [1202.0957].

The Bayesian construction in the cited work is designed to respect the same invariances while acknowledging non-identifiability. Writing the scale-invariant slope as \(\tilde\beta=\beta/l\), the prior
\[
p(\tilde\beta)=\frac{1}{\pi}\frac{1}{1+\tilde\beta^2}
\]
is rotationally invariant, equivalently uniform on the angle \(\theta=\arctan(\tilde\beta)\). After integrating out \(\mu_1,\alpha,\tau^2,\sigma_1^2,\sigma_2^2\), the marginal posterior is
\[
p(\beta\mid y)
=
\frac{p(\beta)\,J(\beta,\nu,r,l)}
{\int_{-\infty}^{\infty} p(\beta)\,J(\beta,\nu,r,l)\,d\beta},
\qquad \nu=n-1,
\]
and depends on the data only through \(r\) and \(l\). With the Cauchy prior on \(\tilde\beta\), the posterior is invariant to interchange and scaling of coordinates. The posterior mass lies between the two OLS slopes \(b_{Y|X}=rl\) and \(b_{X|Y}^{-1}=l/r\), and diminishes rapidly outside that interval [1202.0957].

This Bayesian treatment also changes the interpretation of the point estimate often called Isserlis’s estimator. The paper states that the “Isserlis/RMA” point estimator itself is unbiased only under special symmetry conditions on errors, and cautions against relying on point estimators without specifying error variances. In its simulation study of 90% intervals from 1000 datasets, shortest posterior probability intervals from \(p(\beta\mid y)\) showed broader coverage accuracy across a range of \((\sigma_1,\sigma_2)\) than intervals centered on interchange-invariant point estimators. The reported examples include \(96.6\%\) coverage for the posterior versus approximately \(88\%\) for geometric mean, bisector, and orthogonal regression when \(n=100\) and \(\sigma_1=\sigma_2=0.2\), and \(42.1\%\) versus approximately \(0\%\) when \(n=100\), \(\sigma_1=1.0\), and \(\sigma_2=0.05\) [1202.0957].

The same source introduces the R package `leiv` for computing the posterior density \(p(\beta\mid y)\), posterior summaries, shortest credible intervals, and density plots. Its applications include Zellner’s artificial data, the Faber–Jackson relation in astronomy, and a method-comparison study of milk fat. In the astronomy example, the posterior median is reported as approximately \(3.6\) with shortest 95% interval \((1.8,6.1)\); in the milk-fat comparison, the posterior median is approximately \(0.972\) with shortest 95% interval \((0.953,0.991)\) [1202.0957].

Taken together, the two literatures attach the same label to conceptually different estimators. In Gaussian moment-tensor estimation, Isserlis’s estimator is a covariance plug-in construction that eliminates higher-order chaos and inherits covariance-level concentration. In errors-in-variables line fitting, the name is used for the invariant reduced major axis slope \(b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)\), which is natural under interchange and scaling symmetries but does not resolve the non-identifiability of the structural slope under unspecified error variances.

Source: https://www.emergentmind.com/topics/isserlis-s-estimator