Average Gradient Outer Product (AGOP)

• Average Gradient Outer Product (AGOP) is a symmetric positive semidefinite matrix that aggregates local input gradients to identify task-relevant directions and subspaces.
• It plays a key role in kernel methods and deep neural networks by enabling efficient feature extraction and sufficient dimension reduction.
• Estimation methods for AGOP range from direct empirical averaging and finite differences to smoothed variants, facilitating both representation recovery and adaptive optimization.

Searching arXiv for recent and foundational papers on Average Gradient Outer Product and closely related terms. Average Gradient Outer Product (AGOP) denotes a gradient-second-moment matrix of the form

$\AGOP(f, X) \triangleq \frac{1}{N} \sum_{c=1}^K \sum_{i=1}^N \frac{d f(x_{ci})}{d x} \frac{d f(x_{ci})}{d x}^\top,$

or, in scalar-output settings, its population analogue

$G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$

Across the literature, closely related terminology includes the expected gradient outer product (EGOP) and, for vector-valued outputs, the expected Jacobian outerproduct (EJOP) (Beaglehole et al., 2024, Trivedi et al., 2020). AGOP is a symmetric positive semidefinite matrix that aggregates local input sensitivity over a dataset or distribution, and its leading eigenspaces are used to identify task-relevant directions, central subspaces, learned metrics, and feature-learning mechanisms in kernel methods and deep networks (Beaglehole et al., 2024, Zhu et al., 14 May 2026, DePavia et al., 3 Feb 2025).

1. Definitions and core operator-theoretic structure

In the scalar case, the standard population object is the expected gradient outer product

$$\mathbb{E}_{X} G(X) \triangleq \mathbb{E}_X\!\left(\nabla f(X)\cdot \nabla f(X)^\top\right),$$

with $G(X)=\nabla f(X)\nabla f(X)^\top$ (Trivedi et al., 2020). In multiclass settings, the corresponding Jacobian-based generalization is

$$\mathbb{E}_{X}\!\left(\mathbf J_f(X)\mathbf J_f(X)^\top\right),$$

where

$$\mathbf J_f(x)\mathbf J_f(x)^\top = \sum_{k=1}^c \nabla f_k(x)\nabla f_k(x)^\top$$

for $f=(f_1,\dots,f_c)$ (Trivedi et al., 2020). The 2024 deep-learning literature uses the name AGOP for the empirical Jacobian version and describes it as an uncentered covariance matrix of gradients, since it averages $\nabla f\,\nabla f^\top$ directly without subtracting the mean gradient (Beaglehole et al., 2024).

The matrix is symmetric by construction. This is operationally important because if one differentiates a scalar objective $F(M)$ with respect to an AGOP-type matrix variable

$$M=\mathbb{E}[\nabla f(x)\nabla f(x)^\top],\qquad \hat M=\frac1n\sum_{i=1}^n \nabla f(x_i)\nabla f(x_i)^T,$$

then the correct matrix gradient on the space of symmetric matrices is the orthogonal projection of the ambient gradient onto the symmetric subspace: $G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$0 The formula

$G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$1

is shown to be incorrect as a gradient formula on $G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$2; the correct statement is

$G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$3

with respect to the Frobenius geometry (Srinivasan et al., 2019).

A closely related perspective appears in optimization. The EGOP of a scalar objective $G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$4 with respect to a sampling distribution $G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$5 is

$G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$6

with empirical estimator

$G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$7

In that setting, the eigenvectors identify directions of average squared directional derivative, since

$G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$8

This suggests that AGOP is not only a representation-learning object but also a coordinate-selection object for basis-sensitive optimization methods (DePavia et al., 3 Feb 2025).

2. Subspace recovery and sufficient dimension reduction

A foundational role of AGOP lies in multi-index regression and sufficient dimension reduction. If

$G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$9

then the gradient lies in the relevant low-dimensional subspace: $$\mathbb{E}_{X} G(X) \triangleq \mathbb{E}_X\!\left(\nabla f(X)\cdot \nabla f(X)^\top\right),$$0 Consequently,

$$\mathbb{E}_{X} G(X) \triangleq \mathbb{E}_X\!\left(\nabla f(X)\cdot \nabla f(X)^\top\right),$$1

so its range lies in $$\mathbb{E}_{X} G(X) \triangleq \mathbb{E}_X\!\left(\nabla f(X)\cdot \nabla f(X)^\top\right),$$2 (Zhu et al., 14 May 2026). Under mild assumptions, “the column space of the EGOP is exactly the relevant subspace” (Trivedi et al., 2020). This is the basis of the classical dimension-reduction interpretation: if $$\mathbb{E}_{X} G(X) \triangleq \mathbb{E}_X\!\left(\nabla f(X)\cdot \nabla f(X)^\top\right),$$3 does not vary along a direction $$\mathbb{E}_{X} G(X) \triangleq \mathbb{E}_X\!\left(\nabla f(X)\cdot \nabla f(X)^\top\right),$$4, then $$\mathbb{E}_{X} G(X) \triangleq \mathbb{E}_X\!\left(\nabla f(X)\cdot \nabla f(X)^\top\right),$$5 lies in the nullspace of the operator (Trivedi et al., 2020).

Recent work sharpens this picture for kernel regression. In the Boolean-hypercube multi-index model

$$\mathbb{E}_{X} G(X) \triangleq \mathbb{E}_X\!\left(\nabla f(X)\cdot \nabla f(X)^\top\right),$$6

the empirical AGOP of a fitted kernel ridge regressor

$$\mathbb{E}_{X} G(X) \triangleq \mathbb{E}_X\!\left(\nabla f(X)\cdot \nabla f(X)^\top\right),$$7

is proved to recover the central subspace even when prediction remains inaccurate (Zhu et al., 14 May 2026). The central approximation theorem states that, under the stated assumptions and sample regime $$\mathbb{E}_{X} G(X) \triangleq \mathbb{E}_X\!\left(\nabla f(X)\cdot \nabla f(X)^\top\right),$$8 with $$\mathbb{E}_{X} G(X) \triangleq \mathbb{E}_X\!\left(\nabla f(X)\cdot \nabla f(X)^\top\right),$$9,

$G(X)=\nabla f(X)\nabla f(X)^\top$0

where $G(X)=\nabla f(X)\nabla f(X)^\top$1 is the low-degree population AGOP and $G(X)=\nabla f(X)\nabla f(X)^\top$2 is the explicit rate term given in the paper (Zhu et al., 14 May 2026). If a latent truncated AGOP is nondegenerate, the top $G(X)=\nabla f(X)\nabla f(X)^\top$3 eigenvectors of $G(X)=\nabla f(X)\nabla f(X)^\top$4 consistently recover the central subspace via a Davis–Kahan argument (Zhu et al., 14 May 2026).

The paper’s main conceptual claim is a separation between prediction and representation. If the target has degree $G(X)=\nabla f(X)\nabla f(X)^\top$5, then accurate prediction by KRR requires

$G(X)=\nabla f(X)\nabla f(X)^\top$6

but if the degree-$G(X)=\nabla f(X)\nabla f(X)^\top$7 component already contains all subspace information, then subspace recovery occurs already at

$G(X)=\nabla f(X)\nabla f(X)^\top$8

which can be much smaller when $G(X)=\nabla f(X)\nabla f(X)^\top$9 (Zhu et al., 14 May 2026). This suggests that AGOP can recover representational structure before full function approximation.

A related but distinct development replaces the classical EGOP by a smoothed variant. The Expected Smoothed Gradient Outer Product (ESGOP) is

$$\mathbb{E}_{X}\!\left(\mathbf J_f(X)\mathbf J_f(X)^\top\right),$$0

with

$$\mathbb{E}_{X}\!\left(\mathbf J_f(X)\mathbf J_f(X)^\top\right),$$1

and Monte Carlo approximation

$$\mathbb{E}_{X}\!\left(\mathbf J_f(X)\mathbf J_f(X)^\top\right),$$2

This construction is shown to recover the central mean subspace while attaining a parametric $$\mathbb{E}_{X}\!\left(\mathbf J_f(X)\mathbf J_f(X)^\top\right),$$3-type subspace rate under known $$\mathbb{E}_{X}\!\left(\mathbf J_f(X)\mathbf J_f(X)^\top\right),$$4 and suitable moment conditions (Yuan et al., 2023). The paper explicitly positions ESGOP/ASGOP as a smoothed surrogate rather than a direct estimator of the unsmoothed AGOP.

3. Estimation methods and computational realizations

A basic empirical AGOP estimator averages gradient outer products at sample locations: $$\mathbb{E}_{X}\!\left(\mathbf J_f(X)\mathbf J_f(X)^\top\right),$$5 In kernel ridge regression with kernel

$$\mathbb{E}_{X}\!\left(\mathbf J_f(X)\mathbf J_f(X)^\top\right),$$6

the fitted predictor is

$$\mathbb{E}_{X}\!\left(\mathbf J_f(X)\mathbf J_f(X)^\top\right),$$7

and its gradient is

$$\mathbb{E}_{X}\!\left(\mathbf J_f(X)\mathbf J_f(X)^\top\right),$$8

Evaluating these gradients at the training points yields the empirical AGOP used for eigenspace estimation (Zhu et al., 14 May 2026).

For multiclass nonparametric classification, a rough estimator of the Expected Jacobian Outerproduct uses finite differences of a kernel estimate $$\mathbb{E}_{X}\!\left(\mathbf J_f(X)\mathbf J_f(X)^\top\right),$$9. The $$\mathbf J_f(x)\mathbf J_f(x)^\top = \sum_{k=1}^c \nabla f_k(x)\nabla f_k(x)^\top$$0-entry of the Jacobian is estimated by

$$\mathbf J_f(x)\mathbf J_f(x)^\top = \sum_{k=1}^c \nabla f_k(x)\nabla f_k(x)^\top$$1

and the estimated EJOP is obtained by averaging $$\mathbf J_f(x)\mathbf J_f(x)^\top = \sum_{k=1}^c \nabla f_k(x)\nabla f_k(x)^\top$$2 over the sample (Trivedi et al., 2020). The paper proves consistency of the operator estimate, as well as consistency of its eigenvalues and eigenspaces (Trivedi et al., 2020).

The smoothed-gradient literature uses a different pipeline. By Stein’s identity,

$$\mathbf J_f(x)\mathbf J_f(x)^\top = \sum_{k=1}^c \nabla f_k(x)\nabla f_k(x)^\top$$3

for $$\mathbf J_f(x)\mathbf J_f(x)^\top = \sum_{k=1}^c \nabla f_k(x)\nabla f_k(x)^\top$$4, yielding

$$\mathbf J_f(x)\mathbf J_f(x)^\top = \sum_{k=1}^c \nabla f_k(x)\nabla f_k(x)^\top$$5

With density ratio

$$\mathbf J_f(x)\mathbf J_f(x)^\top = \sum_{k=1}^c \nabla f_k(x)\nabla f_k(x)^\top$$6

one has

$$\mathbf J_f(x)\mathbf J_f(x)^\top = \sum_{k=1}^c \nabla f_k(x)\nabla f_k(x)^\top$$7

This allows unbiased estimation of a smoothed gradient functional from response-weighted first moments rather than from direct pointwise gradient estimation (Yuan et al., 2023).

In optimization, empirical EGOP estimation is straightforward: $$\mathbf J_f(x)\mathbf J_f(x)^\top = \sum_{k=1}^c \nabla f_k(x)\nabla f_k(x)^\top$$8 The eigendecomposition

$$\mathbf J_f(x)\mathbf J_f(x)^\top = \sum_{k=1}^c \nabla f_k(x)\nabla f_k(x)^\top$$9

defines an orthonormal reparameterization $f=(f_1,\dots,f_c)$0, which is then optimized by Adagrad or Adam in the EGOP eigenbasis (DePavia et al., 3 Feb 2025). The same paper also proposes block EGOP estimation,

$f=(f_1,\dots,f_c)$1

to make the procedure tractable for large models (DePavia et al., 3 Feb 2025).

4. Feature learning in kernel methods and deep networks

AGOP has become a central explanatory object in feature-learning accounts of deep representation geometry. In deep networks, it is presented as the average of outer products of input-output gradients, and geometrically it measures the directions in input or feature space along which the predictor changes most (Beaglehole et al., 2024). Large eigenvalues correspond to directions of consistent sensitivity, while small eigenvalues correspond to directions the predictor mostly ignores (Beaglehole et al., 2024).

This interpretation is made algorithmic in Recursive Feature Machines (RFM) and Deep RFM. In the Deep RFM recursion, at layer $f=(f_1,\dots,f_c)$2,

$f=(f_1,\dots,f_c)$3

then the AGOP is computed as

$f=(f_1,\dots,f_c)$4

and the next-layer features are formed by

$f=(f_1,\dots,f_c)$5

The random feature map is usually

$f=(f_1,\dots,f_c)$6

with $f=(f_1,\dots,f_c)$7 Gaussian random and $f=(f_1,\dots,f_c)$8 ReLU (Beaglehole et al., 2024).

The paper argues that AGOP drives Deep Neural Collapse (DNC). Its Neural Feature Ansatz is

$f=(f_1,\dots,f_c)$9

Thus the Gram matrix $\nabla f\,\nabla f^\top$0 approximately matches the layerwise AGOP, implying that neural-network feature learning can be interpreted as approximately projection onto the AGOP eigenspaces (Beaglehole et al., 2024). The paper reports that, in standard DNNs, most within-class variability collapse comes from the $\nabla f\,\nabla f^\top$1 part of the SVD

$\nabla f\,\nabla f^\top$2

and that this singular structure is strongly correlated with AGOP (Beaglehole et al., 2024).

In Deep RFM, the same mechanism is explicit. The theorem

$\nabla f\,\nabla f^\top$3

shows contraction of the Gram matrix toward the collapsed class structure in a linearized-kernel high-dimensional setting (Beaglehole et al., 2024). Empirically, the paper states that the improvement in NC1 is “entirely due to $\nabla f\,\nabla f^\top$4,” while the random feature map can worsen NC1 (Beaglehole et al., 2024). This is presented as evidence that the AGOP step, not the random feature step, causes collapse.

A separate non-neural line of work uses AGOP to explain grokking-like emergence in modular arithmetic. Recursive Feature Machines iteratively fit a kernel regressor,

$\nabla f\,\nabla f^\top$5

then update the feature matrix by

$\nabla f\,\nabla f^\top$6

with $\nabla f\,\nabla f^\top$7 in the experiments (Mallinar et al., 2024). The paper reports that training loss is identically zero and train accuracy is $\nabla f\,\nabla f^\top$8 at every iteration, yet test accuracy stays near random before sharply transitioning to perfect test accuracy (Mallinar et al., 2024). This is used to argue that delayed generalization can result purely from feature learning through AGOP rather than from neural architecture or gradient descent.

The same work identifies a learned block-circulant feature structure. For modular addition and subtraction, the learned feature matrix takes the form

$\nabla f\,\nabla f^\top$9

where $F(M)$0 is an asymmetric circulant matrix and

$F(M)$1

for constants $F(M)$2 (Mallinar et al., 2024). The paper further reports high Pearson correlations between the square root of neural-network AGOP and the first-layer Neural Feature Matrix $F(M)$3: Add $F(M)$4, Sub $F(M)$5, Mul $F(M)$6, Div $F(M)$7 (Mallinar et al., 2024). This is presented as further evidence for AGOP as a mechanism of feature learning in neural networks.

5. Geometry, basis selection, and related outer-product constructions

Because AGOP is symmetric positive semidefinite, its eigendecomposition naturally defines anisotropic metrics. In EJOP-based multiclass classification, if

$F(M)$8

one may transform

$F(M)$9

equivalently inducing a Mahalanobis-type distance that emphasizes predictive directions (Trivedi et al., 2020). The paper reports improvements in nonparametric classification by using the estimated EJOP as a metric and as initialization for metric learning (Trivedi et al., 2020).

In adaptive optimization, the same eigenspace logic appears in parameter space. If

$$M=\mathbb{E}[\nabla f(x)\nabla f(x)^\top],\qquad \hat M=\frac1n\sum_{i=1}^n \nabla f(x_i)\nabla f(x_i)^T,$$0

then optimization is performed on

$$M=\mathbb{E}[\nabla f(x)\nabla f(x)^\top],\qquad \hat M=\frac1n\sum_{i=1}^n \nabla f(x_i)\nabla f(x_i)^T,$$1

The paper shows that, for isotropic $$M=\mathbb{E}[\nabla f(x)\nabla f(x)^\top],\qquad \hat M=\frac1n\sum_{i=1}^n \nabla f(x_i)\nabla f(x_i)^T,$$2,

$$M=\mathbb{E}[\nabla f(x)\nabla f(x)^\top],\qquad \hat M=\frac1n\sum_{i=1}^n \nabla f(x_i)\nabla f(x_i)^T,$$3

so in the EGOP eigenbasis the expected outer product becomes diagonal (DePavia et al., 3 Feb 2025). The same paper defines a stable-rank quantity

$$M=\mathbb{E}[\nabla f(x)\nabla f(x)^\top],\qquad \hat M=\frac1n\sum_{i=1}^n \nabla f(x_i)\nabla f(x_i)^T,$$4

and argues that strong EGOP spectral decay can substantially improve the convergence behavior of basis-sensitive adaptive optimizers (DePavia et al., 3 Feb 2025).

AGOP also has a broader relation to other outer-product-based second-order constructions. In diffusion models, the diffusion Fisher

$$M=\mathbb{E}[\nabla f(x)\nabla f(x)^\top],\qquad \hat M=\frac1n\sum_{i=1}^n \nabla f(x_i)\nabla f(x_i)^T,$$5

is shown to admit an exact outer-product span-space characterization: $$M=\mathbb{E}[\nabla f(x)\nabla f(x)^\top],\qquad \hat M=\frac1n\sum_{i=1}^n \nabla f(x_i)\nabla f(x_i)^T,$$6 in the Dirac setting, with an analogous integral formula in the general case (Wang et al., 29 May 2025). That work explicitly states that it is not a direct AGOP identity, but an exact pointwise outer-product decomposition of a Hessian/Fisher object that is “analogous to AGOP in spirit” (Wang et al., 29 May 2025).

A more indirect connection appears in neural-network derivative structure. For feedforward networks, the per-sample input gradient is

$$M=\mathbb{E}[\nabla f(x)\nabla f(x)^\top],\qquad \hat M=\frac1n\sum_{i=1}^n \nabla f(x_i)\nabla f(x_i)^T,$$7

so the per-sample AGOP contribution is immediately $$M=\mathbb{E}[\nabla f(x)\nabla f(x)^\top],\qquad \hat M=\frac1n\sum_{i=1}^n \nabla f(x_i)\nabla f(x_i)^T,$$8 (Bakker et al., 2018). That paper does not define AGOP explicitly, but it shows that per-sample derivatives in fully connected and recurrent networks have an outer-product structure, which makes AGOP-like constructions natural and computationally accessible in those architectures (Bakker et al., 2018).

6. Variants, misconceptions, and scope of the concept

AGOP is not a single method but a family of closely related operators. The literature distinguishes at least four recurring variants.

Variant	Representative form	Typical setting
EGOP / population AGOP	$$M=\mathbb{E}[\nabla f(x)\nabla f(x)^\top],\qquad \hat M=\frac1n\sum_{i=1}^n \nabla f(x_i)\nabla f(x_i)^T,$$9	Scalar regression, active subspaces, optimization
Empirical AGOP	$G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$00	KRR, RFM, feature learning
EJOP	$G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$01	Multiclass/vector-valued outputs
Smoothed AGOP variants	ESGOP / ASGOP	SDR with fast rates

One recurring misconception concerns symmetry-aware differentiation. Since AGOP-type matrices are symmetric by construction, matrix optimization over AGOP variables often proceeds on the symmetric subspace. The correct matrix gradient is

$G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$02

not the doubled off-diagonal expression

$G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$03

sometimes quoted in statistics or control texts (Srinivasan et al., 2019). For AGOP objectives involving $G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$04, $G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$05, spectral penalties, or covariance-like regularizers, this distinction determines the correct steepest-descent direction in Frobenius geometry (Srinivasan et al., 2019).

A second misconception is to conflate AGOP with any outer-product-based approximation. For example, approximate outer product gradient descent with memory represents the minibatch weight gradient as

$G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$06

and approximates it by a subset of rank-one terms (Hernandez et al., 2021). This is directly related to outer-product decompositions of backpropagation, but it is not an AGOP method in the usual sense because the object being approximated is the weight-gradient matrix itself, not a covariance-like matrix of input gradients with themselves (Hernandez et al., 2021).

A third distinction concerns the regression target. Classical OPG methods for sufficient dimension reduction target the conditional mean via averaged outer products of mean gradients. Local Modal OPG (LMOPG) instead estimates

$G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$07

for the conditional mode function $G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$08, and averages local modal gradient outer products

$G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$09

to recover a central subspace under skew-error mechanisms (Li et al., 2024). This suggests that the outer-product-of-gradients principle extends beyond mean regression, but the target subspace then changes accordingly.

A plausible implication is that AGOP is best viewed as a geometric template rather than a single estimator: choose a target function or predictor, define local derivatives in the relevant space, form a second-moment operator by averaging their outer products, and use the resulting spectrum to identify directions of relevance. The exact meaning of “relevance” depends on the model class, output type, smoothing scheme, and sampling distribution (Trivedi et al., 2020, Yuan et al., 2023, DePavia et al., 3 Feb 2025).

7. Significance and current directions

Several contemporary strands of research converge on AGOP as a unifying operator between statistical dimension reduction, kernel feature learning, deep representation geometry, and basis-adaptive optimization.

In kernel regression, AGOP now has a theorem-level role as a representation-recovery mechanism that can succeed in lower-sample regimes than full prediction (Zhu et al., 14 May 2026). In Deep RFM and related deep-learning analyses, AGOP is treated as the data-dependent feature-learning object that shapes learned representations and can drive Deep Neural Collapse (Beaglehole et al., 2024). In non-neural grokking studies, AGOP-driven feature updates are used to isolate emergence-like delayed generalization from gradient-descent dynamics and from training-loss effects (Mallinar et al., 2024). In optimization, EGOP provides a principled orthonormal basis for reparameterizing basis-sensitive adaptive methods such as Adagrad and Adam (DePavia et al., 3 Feb 2025).

At the same time, the literature is careful about limitations. Classical AGOP estimators based on pointwise gradient estimation can have nonparametric convergence rates under general $G(f):=\mathbb E[\nabla f(x)\nabla f(x)^\top].$10, motivating smoothed surrogates such as ESGOP/ASGOP (Yuan et al., 2023). Some theoretical results are asymptotic and model-specific, as in Deep RFM’s high-dimensional kernel linearization analysis (Beaglehole et al., 2024). In optimization, the benefits of EGOP reparameterization depend on spectral decay and on the up-front cost of estimating and diagonalizing the empirical matrix (DePavia et al., 3 Feb 2025). And in neural architectures, the clean per-sample outer-product derivative structure emphasized for fully connected and recurrent models does not hold in the same way for convolutional layers (Bakker et al., 2018).

The cumulative picture is that AGOP is a symmetric PSD operator encoding average local sensitivity. In statistics, it recovers low-dimensional predictive structure; in kernels, it updates geometry; in deep learning, it provides a concrete mechanism for feature learning; and in optimization, it supplies a data-adapted eigenbasis for coordinate selection. This suggests that AGOP occupies a central place among derivative-based second-moment methods, with current work extending it across output types, smoothing schemes, architectures, and computational regimes (Trivedi et al., 2020, Beaglehole et al., 2024, Zhu et al., 14 May 2026, DePavia et al., 3 Feb 2025).