Isserlis's Estimator: Gaussian Moments & RMA

• Isserlis's Estimator is a dual-use procedure that plugs sample covariance into Wick–Isserlis expansions for even-order Gaussian moments.
• It provides dimension-free non-asymptotic error bounds and improved sample complexity compared to standard sample moment estimators.
• In errors-in-variables models, it refers to the reduced major axis (RMA) slope, ensuring invariance under scaling and interchange of axes.

“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$ (Al-Ghattas et al., 8 Jul 2025). 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 (Leonard, 2012).

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 (Al-Ghattas et al., 8 Jul 2025).

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) (Leonard, 2012).

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

$b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$0

with entries

$b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$1

For zero-mean Gaussian vectors, all odd central moments vanish, $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$2 whenever $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$3 is odd, while for even $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$4 the tensor is completely determined by $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$5 through Isserlis’s theorem (Al-Ghattas et al., 8 Jul 2025).

For centered Gaussian $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$6, Isserlis’s theorem gives, for even $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$7,

$b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$8

and the expectation is zero for odd $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$9. Here $\Sigma$0 is the set of all perfect matchings of $\Sigma$1, with cardinality $\Sigma$2. For $\Sigma$3,

$\Sigma$4

Tensorially,

$\Sigma$5

where $\Sigma$6 inserts covariance entries slotwise according to the pairing $\Sigma$7; the resulting tensor is symmetric in its $\Sigma$8 modes (Al-Ghattas et al., 8 Jul 2025).

Given i.i.d. samples $\Sigma$9 from $\Sigma$0, define

$\Sigma$1

For even-order central moments, the Isserlis estimator is obtained by substituting $\Sigma$2 into Wick’s formula: $\Sigma$3 or, in tensor notation,

$\Sigma$4

For odd $\Sigma$5, it is set to $\Sigma$6 for a centered Gaussian target. For raw moments with $\Sigma$7, one reconstructs

$\Sigma$8

where $\Sigma$9 is computed by Wick’s formula for even $$b_{\mathrm{RMA}}=\operatorname{sign}(r)\frac{s_y}{s_x}.$$0 and set to zero for odd $$b_{\mathrm{RMA}}=\operatorname{sign}(r)\frac{s_y}{s_x}.$$1 (Al-Ghattas et al., 8 Jul 2025).

The contrast with the standard sample moment estimator is structural. The standard estimator

$$b_{\mathrm{RMA}}=\operatorname{sign}(r)\frac{s_y}{s_x}.$$2

targets the $$b_{\mathrm{RMA}}=\operatorname{sign}(r)\frac{s_y}{s_x}.$$3-th Gaussian chaos directly, whereas the Isserlis estimator targets covariance first and then maps it to the $$b_{\mathrm{RMA}}=\operatorname{sign}(r)\frac{s_y}{s_x}.$$4-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 (Al-Ghattas et al., 8 Jul 2025).

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

The analysis is carried out in the operator and entrywise maximum norms. For an order-$$b_{\mathrm{RMA}}=\operatorname{sign}(r)\frac{s_y}{s_x}.$$5 tensor $$b_{\mathrm{RMA}}=\operatorname{sign}(r)\frac{s_y}{s_x}.$$6,

$$b_{\mathrm{RMA}}=\operatorname{sign}(r)\frac{s_y}{s_x}.$$7

and

$$b_{\mathrm{RMA}}=\operatorname{sign}(r)\frac{s_y}{s_x}.$$8

The bounds depend on effective dimensions rather than ambient dimension $$b_{\mathrm{RMA}}=\operatorname{sign}(r)\frac{s_y}{s_x}.$$9: $X\in\mathbb{R}^d$0 The paper describes the resulting guarantees as dimension-free because they contain no explicit dependence on the ambient dimension (Al-Ghattas et al., 8 Jul 2025).

In the symmetric centered case $X\in\mathbb{R}^d$1 with even $X\in\mathbb{R}^d$2 and target $X\in\mathbb{R}^d$3, the expected operator-norm error satisfies

$X\in\mathbb{R}^d$4

while

$X\in\mathbb{R}^d$5

Entrywise, the analogous statements are

$X\in\mathbb{R}^d$6

and

$X\in\mathbb{R}^d$7

Quantity	Standard sample moment	Isserlis estimator
Operator norm	$X\in\mathbb{R}^d$8	$X\in\mathbb{R}^d$9
Entrywise maximum norm	$X\sim\mathcal{N}(\mu,\Sigma)$0	$X\sim\mathcal{N}(\mu,\Sigma)$1
Consistency threshold	$X\sim\mathcal{N}(\mu,\Sigma)$2 or $X\sim\mathcal{N}(\mu,\Sigma)$3	$X\sim\mathcal{N}(\mu,\Sigma)$4 or $X\sim\mathcal{N}(\mu,\Sigma)$5

These bounds quantify a strict improvement. The standard estimator is consistent only when $X\sim\mathcal{N}(\mu,\Sigma)$6 in operator norm, or $X\sim\mathcal{N}(\mu,\Sigma)$7 entrywise, because the term $X\sim\mathcal{N}(\mu,\Sigma)$8 or $X\sim\mathcal{N}(\mu,\Sigma)$9 must be negligible. By contrast, the Isserlis estimator is consistent already when $p$0 or $p$1, inheriting the sample-covariance rates (Al-Ghattas et al., 8 Jul 2025).

The mechanism is summarized by the perturbation control

$p$2

This expresses the moment tensor as a smooth functional of $p$3, so its estimation error is governed by covariance perturbation rather than by direct $p$4-th-chaos fluctuations. For $p$5 and $p$6, the paper gives explicit operator- and entrywise-norm bounds of the same form, with the higher-order correction terms $p$7 and $p$8, respectively (Al-Ghattas et al., 8 Jul 2025).

4. Asymmetric tensors, assumptions, and computation

The same plug-in principle extends to asymmetric Gaussian tensors. Let

$p$9

be zero-mean Gaussian with block covariances

$M_p:=\mathbb{E}[X^{\otimes p}],$0

and target

$M_p:=\mathbb{E}[X^{\otimes p}],$1

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 (Al-Ghattas et al., 8 Jul 2025).

More precisely, provided $M_p:=\mathbb{E}[X^{\otimes p}],$2,

$M_p:=\mathbb{E}[X^{\otimes p}],$3

whereas the standard estimator satisfies a larger bound containing both a sum term and a product term. Entrywise, provided $M_p:=\mathbb{E}[X^{\otimes p}],$4,

$M_p:=\mathbb{E}[X^{\otimes p}],$5

The stated implication is that $M_p:=\mathbb{E}[X^{\otimes p}],$6 is consistent once $M_p:=\mathbb{E}[X^{\otimes p}],$7 dominates the largest effective dimension among blocks, while $M_p:=\mathbb{E}[X^{\otimes p}],$8 is consistent only when $M_p:=\mathbb{E}[X^{\otimes p}],$9 dominates the product term (Al-Ghattas et al., 8 Jul 2025).

The scope conditions are explicit. The results require i.i.d. Gaussian samples, exact Gaussianity, and even $$(M_p)_{j_1,\dots,j_p}=\mathbb{E}[X_{j_1}\cdots X_{j_p}]$$0 for nontrivial central moments. In the symmetric case, $$(M_p)_{j_1,\dots,j_p}=\mathbb{E}[X_{j_1}\cdots X_{j_p}]$$1 is assumed for the cleanest formulation; nonzero mean is handled by centering with $$(M_p)_{j_1,\dots,j_p}=\mathbb{E}[X_{j_1}\cdots X_{j_p}]$$2 and adding mean terms. For odd $$(M_p)_{j_1,\dots,j_p}=\mathbb{E}[X_{j_1}\cdots X_{j_p}]$$3, central moments vanish for centered Gaussians, and the Isserlis estimator is identically zero (Al-Ghattas et al., 8 Jul 2025).

The computational bottleneck is combinatorial: the number of pairings is $$(M_p)_{j_1,\dots,j_p}=\mathbb{E}[X_{j_1}\cdots X_{j_p}]$$4. The paper therefore recommends exploiting symmetry, evaluating contractions such as

$$(M_p)_{j_1,\dots,j_p}=\mathbb{E}[X_{j_1}\cdots X_{j_p}]$$5

as sums of products of quadratic forms $$(M_p)_{j_1,\dots,j_p}=\mathbb{E}[X_{j_1}\cdots X_{j_p}]$$6, and using the recurrence

$$(M_p)_{j_1,\dots,j_p}=\mathbb{E}[X_{j_1}\cdots X_{j_p}]$$7

for dynamic programming. Exhaustive pairing enumeration is described as typically feasible for moderate $$(M_p)_{j_1,\dots,j_p}=\mathbb{E}[X_{j_1}\cdots X_{j_p}]$$8, such as $$(M_p)_{j_1,\dots,j_p}=\mathbb{E}[X_{j_1}\cdots X_{j_p}]$$9 or $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$00, while contraction-based evaluation is preferable for larger $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$01. 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 (Al-Ghattas et al., 8 Jul 2025).

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

$b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$02

with latent relation $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$03. In the structural normal model, $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$04, and the total errors $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$05 are independent, mean-zero normal with variances $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$06. The observed pairs are then bivariate normal with covariance

$b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$07

(Leonard, 2012).

If $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$08 is the sample covariance matrix with divisor $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$09, define

$b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$10

The reduced major axis or standard major axis slope is

$b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$11

It depends only on $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$12 and $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$13, is scale invariant, and under interchange of axes maps to its reciprocal $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$14. The same functional form can be written as the geometric mean of the two OLS slopes,

$b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$15

where

$b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$16

Samuelson’s result, as summarized in the cited source, is that among point estimates depending only on $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$17 and $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$18 and respecting scale and interchange invariance, the only possibility is this geometric mean form (Leonard, 2012).

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) (Leonard, 2012).

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

The structural normal model does not identify $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$19 without additional information. Reiersøl (1950) showed that the sampling distribution identifies $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$20 but not $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$21, 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 (Leonard, 2012).

The Bayesian construction in the cited work is designed to respect the same invariances while acknowledging non-identifiability. Writing the scale-invariant slope as $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$22, the prior

$b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$23

is rotationally invariant, equivalently uniform on the angle $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$24. After integrating out $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$25, the marginal posterior is

$b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$26

and depends on the data only through $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$27 and $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$28. With the Cauchy prior on $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$29, the posterior is invariant to interchange and scaling of coordinates. The posterior mass lies between the two OLS slopes $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$30 and $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$31, and diminishes rapidly outside that interval (Leonard, 2012).

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 $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$32 showed broader coverage accuracy across a range of $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$33 than intervals centered on interchange-invariant point estimators. The reported examples include $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$34 coverage for the posterior versus approximately $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$35 for geometric mean, bisector, and orthogonal regression when $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$36 and $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$37, and $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$38 versus approximately $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$39 when $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$40, $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$41, and $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$42 (Leonard, 2012).

The same source introduces the R package leiv for computing the posterior density $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$43, 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 $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$44 with shortest 95% interval $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$45; in the milk-fat comparison, the posterior median is approximately $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$46 with shortest 95% interval $b_{\mathrm{RMA}}=\operatorname{sign}(r)(s_y/s_x)$47 (Leonard, 2012).

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)$48, which is natural under interchange and scaling symmetries but does not resolve the non-identifiability of the structural slope under unspecified error variances.