---
title: Average Gradient Outer Product (AGOP)
url: https://www.emergentmind.com/topics/average-gradient-outer-product
type: topic
---

# Average Gradient Outer Product (AGOP)

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) [2402.13728], [2006.03550]. 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 [2402.13728], [2605.15082], [2502.01594].

## 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\) [2006.03550]. 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)\) [2006.03550]. 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 [2402.13728].

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_{\mathrm{sym}}=\operatorname{sym}(G)=\frac{G+G^T}{2}.
\]
The formula
\[
G_s=G+G^T-G\circ I
\]
is shown to be incorrect as a gradient formula on \(\mathbb S^{n\times n}\); the correct statement is
\[
\nabla \phi_{\mathrm{sym}}(A)=\operatorname{sym}(\nabla \phi(A))
\]
with respect to the Frobenius geometry [1911.06491].

A closely related perspective appears in optimization. The EGOP of a scalar objective \(f:\mathbb R^d\to\mathbb R\) with respect to a sampling distribution \(\rho\) is
\[
\EGOP(f) \defeq \mathbb{E}_{\theta\sim \rho}\left[\nabla f(\theta)\nabla f(\theta)^{\T}\right],
\]
with empirical estimator
\[
\hat{P} = \frac{1}{M}\sum_{i=1}^M \nabla f(\theta_i)\nabla f(\theta_i)^\T.
\]
In that setting, the eigenvectors identify directions of average squared directional derivative, since
\[
\lambda_i = \mathbb{E}_{\theta\sim\rho}[\langle \nabla f(\theta), v_i\rangle^2 ].
\]
This suggests that AGOP is not only a representation-learning object but also a coordinate-selection object for basis-sensitive optimization methods [2502.01594].

## 2. Subspace recovery and sufficient dimension reduction

A foundational role of AGOP lies in multi-index regression and sufficient dimension reduction. If
\[
f(x)=g(V^\top x)
\quad\text{or}\quad
f^*(x)=h(Ux),
\]
then the gradient lies in the relevant low-dimensional subspace:
\[
\nabla f^*(x)=U^\top \nabla h(Ux).
\]
Consequently,
\[
\mathbb E[\nabla f^*(x)\nabla f^*(x)^\top]
=
U^\top \mathbb E[\nabla h(Ux)\nabla h(Ux)^\top]U,
\]
so its range lies in \(\operatorname{row}(U)\) [2605.15082]. Under mild assumptions, “the column space of the EGOP is exactly the relevant subspace” [2006.03550]. This is the basis of the classical dimension-reduction interpretation: if \(f\) does not vary along a direction \(v\), then \(v\) lies in the nullspace of the operator [2006.03550].

Recent work sharpens this picture for kernel regression. In the Boolean-hypercube multi-index model
\[
x^{(i)} \overset{\text{i.i.d.}{\sim} \tau_d,\qquad
y_i = f^*(x^{(i)})+\varepsilon_i,\qquad
f^*(x)=h(Ux),
\]
the empirical AGOP of a fitted kernel ridge regressor
\[
\widehat M := \frac1n\sum_{i=1}^n \nabla \hat f(x^{(i)})\nabla \hat f(x^{(i)})^\top
\]
is proved to recover the central subspace even when prediction remains inaccurate [2605.15082]. The central approximation theorem states that, under the stated assumptions and sample regime \(n=d^{p+\delta}\) with \(\delta\in(0,1)\),
\[
\|\widehat M-M_{\le p}\|_{\rm op}
=
O_{d,\mathbb P}\!\Big(
R_d(U)\,(\|f^*\|_{L_2}^2+\sigma_\varepsilon^2)
\Big),
\]
where \(M_{\le p}\) is the low-degree population AGOP and \(R_d(U)\) is the explicit rate term given in the paper [2605.15082]. If a latent truncated AGOP is nondegenerate, the top \(r\) eigenvectors of \(\widehat M\) consistently recover the central subspace via a Davis–Kahan argument [2605.15082].

The paper’s main conceptual claim is a separation between prediction and representation. If the target has degree \(p^*\), then accurate prediction by KRR requires
\[
n\asymp d^{p^*},
\]
but if the degree-\(\le p\) component already contains all subspace information, then subspace recovery occurs already at
\[
n=d^{p+\delta},\qquad \delta\in(0,1),
\]
which can be much smaller when \(p<p^*\) [2605.15082]. 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
\[
M :=\mathbb{E}_{\theta\sim N(0_d,\sigma_\theta^2 I_d)}
\big[\beta_h(\theta)\beta_h(\theta)^\top\big],
\]
with
\[
\beta_h(\theta)=\mathbb{E}_{Z\sim N(\theta,h^2I_d)}[\nabla g(Z)],
\]
and Monte Carlo approximation
\[
\widetilde M :=\frac1m\sum_{j=1}^m \beta_h(\theta_j)\beta_h(\theta_j)^\top.
\]
This construction is shown to recover the central mean subspace while attaining a parametric \(n^{-1/2}\)-type subspace rate under known \(P_X\) and suitable moment conditions [2312.15469]. 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:
\[
\hat M=\frac1n\sum_{i=1}^n \nabla \hat f(x^{(i)})\nabla \hat f(x^{(i)})^\top.
\]
In kernel ridge regression with kernel
\[
K(x,x')=g\!\left(\frac{\langle x,x'\rangle}{d}\right),
\]
the fitted predictor is
\[
\hat f(x)=K(x,X)(K(X,X)+\lambda I_n)^{-1}y,
\]
and its gradient is
\[
\nabla \hat f(x)
=
\frac1d\sum_{j=1}^n \alpha_j g'\!\left(\frac{\langle x,x^{(j)}\rangle}{d}\right)x^{(j)},
\qquad
\alpha=(K+\lambda I)^{-1}y.
\]
Evaluating these gradients at the training points yields the empirical AGOP used for eigenspace estimation [2605.15082].

For multiclass nonparametric classification, a rough estimator of the Expected Jacobian Outerproduct uses finite differences of a kernel estimate \(\bar f_{n,h}(x)\). The \((i,j)\)-entry of the Jacobian is estimated by
\[
\Delta_{t,i,j} f_n(x)=\frac{f_{n,i}(x+t e_j)-f_{n,i}(x-t e_j)}{2t},
\]
and the estimated EJOP is obtained by averaging \(\hat{\mathbf J}_f(x)\hat{\mathbf J}_f(x)^\top\) over the sample [2006.03550]. The paper proves consistency of the operator estimate, as well as consistency of its eigenvalues and eigenspaces [2006.03550].

The smoothed-gradient literature uses a different pipeline. By Stein’s identity,
\[
\mathbb{E}[g(Z)(Z-\theta)] = h^2 \mathbb{E}[\nabla g(Z)]
\]
for \(Z\sim N(\theta,h^2I_d)\), yielding
\[
\beta_h(\theta)=h^{-2}\mathbb{E}_{Z\sim N(\theta,h^2I_d)}[g(Z)(Z-\theta)].
\]
With density ratio
\[
\rho_h(x;\theta):=\frac{\varphi_h(x-\theta)}{p(x)},
\]
one has
\[
\mathbb{E}\big[\rho_h(X;\theta)\,Y\,(X-\theta)\big] = h^2 \beta_h(\theta).
\]
This allows unbiased estimation of a smoothed gradient functional from response-weighted first moments rather than from direct pointwise gradient estimation [2312.15469].

In optimization, empirical EGOP estimation is straightforward:
\[
\hat{P} = \frac{1}{M}\sum_{i=1}^M \nabla f(\theta_i)\nabla f(\theta_i)^\T,
\qquad \theta_i \overset{\text{i.i.d.}{\sim} \rho.
\]
The eigendecomposition
\[
\hat P = V\Lambda V^\top
\]
defines an orthonormal reparameterization \(\tilde f(\tilde\theta)=f(V\tilde\theta)\), which is then optimized by Adagrad or Adam in the EGOP eigenbasis [2502.01594]. The same paper also proposes block EGOP estimation,
\[
\EGOP^{(i)} \defeq \frac{1}{M}\sum_{k=1}^M \nabla_{S_i} f(\theta_k)\nabla_{S_i} f(\theta_k)^\T,
\]
to make the procedure tractable for large models [2502.01594].

## 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 [2402.13728]. Large eigenvalues correspond to directions of consistent sensitivity, while small eigenvalues correspond to directions the predictor mostly ignores [2402.13728].

This interpretation is made algorithmic in Recursive Feature Machines (RFM) and Deep RFM. In the Deep RFM recursion, at layer \(l\),
\[
\alpha_l =Y(k_l(X_l, X_l)+\mu I)^{-1},
\qquad
f_l(\cdot) = \alpha_lk_l(\cdot, X_l),
\]
then the AGOP is computed as
\[
M_l = \sum_{c,i=1}^{K,n} \nabla f_l(x^l_{ci}) \nabla f_l(x^l_{ci})^\top,
\]
and the next-layer features are formed by
\[
X_{l+1} \xleftarrow[]{} \Phi_l(M_l^{1/2}X_l).
\]
The random feature map is usually
\[
\Phi_l(\cdot)=\sigma(W_l \cdot),
\]
with \(W_l\) Gaussian random and \(\sigma\) ReLU [2402.13728].

The paper argues that AGOP drives Deep Neural Collapse (DNC). Its Neural Feature Ansatz is
\[
W_l^\top W_l\, \propto\, \frac{1}{N} \sum_{c=1}^K \sum_{i=1}^n \frac{d f(x_{ci})}{d x^l} \frac{d f(x_{ci})}{d x^l}^\top.
\]
Thus the Gram matrix \(W_l^\top W_l\) approximately matches the layerwise AGOP, implying that neural-network feature learning can be interpreted as approximately projection onto the AGOP eigenspaces [2402.13728]. The paper reports that, in standard DNNs, most within-class variability collapse comes from the \(S_lV_l^\top\) part of the SVD
\[
W_l=U_lS_lV_l^\top,
\]
and that this singular structure is strongly correlated with AGOP [2402.13728].

In Deep RFM, the same mechanism is explicit. The theorem
\[
\|\tilde{X}_{l+1}^\top \tilde{X}_{l+1} - yy^\top\| \leq (1 - \epsilon) \|\tilde{X}_{l}^\top \tilde{X}_{l} - yy^\top\| + O(\lambda_k^2 \lambda_{\Phi}^{-2})
\]
shows contraction of the Gram matrix toward the collapsed class structure in a linearized-kernel high-dimensional setting [2402.13728]. Empirically, the paper states that the improvement in NC1 is “entirely due to \(M_l^{1/2}\),” while the random feature map can worsen NC1 [2402.13728]. 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,
\[
f^{(t)}(x) = k(x, X ; M_t)\alpha,
\qquad
\alpha = k(X, X; M_t)^{-1} y,
\]
then update the feature matrix by
\[
M_{t+1} = [G(f^{(t)})]^{s},
\]
with \(s=\frac12\) in the experiments [2407.20199]. The paper reports that training loss is identically zero and train accuracy is \(100\%\) at every iteration, yet test accuracy stays near random before sharply transitioning to perfect test accuracy [2407.20199]. 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
\[
M^* = \begin{pmatrix} A  & C\tran\\ C & A \end{pmatrix},
\]
where \(C\) is an asymmetric circulant matrix and
\[
A = c_1 I + c_2 \mathbf 1 \mathbf 1^\top
\]
for constants \(c_1,c_2\) [2407.20199]. The paper further reports high Pearson correlations between the square root of neural-network AGOP and the first-layer Neural Feature Matrix \(W_1^\top W_1\): Add \(0.955\), Sub \(0.942\), Mul \(0.924\), Div \(0.929\) [2407.20199]. 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
\[
M=VDV^\top,
\]
one may transform
\[
x \mapsto D^{1/2}V^\top x,
\]
equivalently inducing a Mahalanobis-type distance that emphasizes predictive directions [2006.03550]. The paper reports improvements in nonparametric classification by using the estimated EJOP as a metric and as initialization for metric learning [2006.03550].

In adaptive optimization, the same eigenspace logic appears in parameter space. If
\[
\EGOP(f)=V\Lambda V^\top,
\]
then optimization is performed on
\[
\tilde f(\tilde\theta)\defeq f(V\tilde\theta).
\]
The paper shows that, for isotropic \(\rho\),
\[
\mathbb{E}_{\theta\sim \rho}[\nabla \tilde f(\theta)\nabla \tilde f(\theta)^\T]
=
V^\T \mathbb{E}_{\theta\sim \rho}[\nabla f(\theta)\nabla f(\theta)^\T] V,
\]
so in the EGOP eigenbasis the expected outer product becomes diagonal [2502.01594]. The same paper defines a stable-rank quantity
\[
{\sr}_f \defeq \frac{\sum_{i=1}^d \sqrt{\lambda_i(\EGOP(f))}}{\sqrt{\lambda_{\max}(\EGOP(f))}},
\]
and argues that strong EGOP spectral decay can substantially improve the convergence behavior of basis-sensitive adaptive optimizers [2502.01594].

AGOP also has a broader relation to other outer-product-based second-order constructions. In diffusion models, the diffusion Fisher
\[
F_t(x_t,t):= -\nabla_{x_t}^2 \log q_t(x_t,t)
\]
is shown to admit an exact outer-product span-space characterization:
\[
F_t(x_t,t) = \frac{1}{\sigma_t^2}I - \frac{\alpha_t^2}{\sigma_t^4} \left[ \sum_i w_i y_i y_i^\top - \left(\sum_i w_i y_i\right)\left(\sum_i w_i y_i\right)^\top \right]
\]
in the Dirac setting, with an analogous integral formula in the general case [2505.23264]. 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” [2505.23264].

A more indirect connection appears in neural-network derivative structure. For feedforward networks, the per-sample input gradient is
\[
\frac{\partial f}{\partial v^{\left(0\right),i}} = \frac{\partial f}{\partial p^j} u^j_l \alpha^{\left(n,0\right),l}_i \equiv \zeta_i,
\]
so the per-sample AGOP contribution is immediately \(\zeta\zeta^\top\) [1810.03798]. 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 [1810.03798].

## 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** | \(\mathbb E[\nabla f(x)\nabla f(x)^\top]\) | Scalar regression, active subspaces, optimization |
| **Empirical AGOP** | \(\frac1n\sum_i \nabla \hat f(x_i)\nabla \hat f(x_i)^\top\) | KRR, RFM, feature learning |
| **EJOP** | \(\mathbb E[\mathbf J_f(X)\mathbf J_f(X)^\top]\) | 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
\[
\nabla_{\mathbb S}J(M)=\operatorname{sym}(G),
\]
not the doubled off-diagonal expression
\[
G+G^T-G\circ I
\]
sometimes quoted in statistics or control texts [1911.06491]. For AGOP objectives involving \(\log\det M\), \(\operatorname{tr}(CM)\), spectral penalties, or covariance-like regularizers, this distinction determines the correct steepest-descent direction in Frobenius geometry [1911.06491].

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
\[
\mathbf{W}_i^{*}=\sum_{m=1}^M \mathbf{x}_{i,m}\mathbf{g}_{i+1,m}^T
\]
and approximates it by a subset of rank-one terms [2110.09164]. 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 [2110.09164].

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
\[
E\left( \frac{\partial m(\bm X)}{\partial \bm X}\frac{\partial m(\bm X)}{\partial \bm X^\top} \right)
\]
for the conditional mode function \(m(\bm X)=\mathrm M(y\mid \bm X)\), and averages local modal gradient outer products
\[
\frac{1}{n}\sum_{j=1}^n \hat{\bm b}_{(j)}\hat{\bm b}_{(j)}^\top
\]
to recover a central subspace under skew-error mechanisms [2407.20738]. 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 [2006.03550], [2312.15469], [2502.01594].

## 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 [2605.15082]. 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 [2402.13728]. 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 [2407.20199]. In optimization, EGOP provides a principled orthonormal basis for reparameterizing basis-sensitive adaptive methods such as Adagrad and Adam [2502.01594].

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 \(P_X\), motivating smoothed surrogates such as ESGOP/ASGOP [2312.15469]. Some theoretical results are asymptotic and model-specific, as in Deep RFM’s high-dimensional kernel linearization analysis [2402.13728]. In optimization, the benefits of EGOP reparameterization depend on spectral decay and on the up-front cost of estimating and diagonalizing the empirical matrix [2502.01594]. 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 [1810.03798].

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 [2006.03550], [2402.13728], [2605.15082], [2502.01594].

Source: https://www.emergentmind.com/topics/average-gradient-outer-product