Recursive Feature Machine (RFM)

• Recursive Feature Machine (RFM) is a kernel method that recursively updates its metric using the Average Gradient Outer Product (AGOP) to uncover task-relevant directions.
• It bridges fixed-kernel approaches and neural networks by combining kernel regression with an iterative feature-learning mechanism that reshapes similarity relations during training.
• Variants of RFM improve robustness, scalability, and interpretability, though theoretical challenges and efficient kernel computations remain active research areas.

Recursive Feature Machine (RFM) denotes a class of kernel methods in which the feature geometry is not fixed in advance, but is updated recursively from the model’s own input sensitivities. In the RFM literature, the central object is the average gradient outer product (AGOP), a positive semidefinite matrix built from predictor gradients; its eigenspaces define task-relevant directions, and the learned matrix $M$ is then used as a Mahalanobis-type metric inside the kernel. RFMs are therefore positioned between fixed-kernel methods and neural networks: they retain kernel-regression structure, but acquire an explicit iterative feature-learning mechanism that reshapes similarity relations during training (Mallinar et al., 2024, Radhakrishnan et al., 2024, Beaglehole et al., 12 Aug 2025).

1. Core mechanism and algorithmic structure

An RFM alternates between predictor fitting and metric update. For a predictor $f:\mathbb{R}^d\to\mathbb{R}^c$, the AGOP is defined in this literature as

$$\mathrm{AGOP}(f;\{x^{(j)}\}_{j=1}^n) = \frac1n\sum_{j=1}^n \frac{\partial f(x^{(j)})}{\partial x} \left(\frac{\partial f(x^{(j)})}{\partial x}\right)^\top \in\mathbb{R}^{d\times d}.$$

This matrix summarizes directions in input space along which the current predictor varies most strongly. Large diagonal entries indicate coordinatewise sensitivity, while leading eigendirections identify low-dimensional structure used by the model (Mallinar et al., 2024).

In the standard kernel-RFM formulation, one starts from a positive definite or positive semidefinite matrix $M_0$, fits a kernel predictor under the current metric, computes AGOP, and replaces the metric by a matrix function of AGOP. A representative iteration used in the literature is

$$f^{(t)}(x)=k(x,X;M_t)\alpha,\qquad \alpha = k(X,X;M_t)^{-1}y,$$

followed by

$M_{t+1}=\bigl[\mathrm{AGOP}(f^{(t)})\bigr]^s,$

with $s>0$, often $s=\tfrac12$. The predictor remains a kernel expansion, but the kernel geometry changes because $M_t$ changes (Mallinar et al., 2024).

This learned matrix can be read as a supervised metric. In the Laplace form used in several papers,

$$K(x,x')=\exp\!\left(-\gamma\sqrt{(x-x')^TM(x-x')}\right),$$

so recursive updates change which directions dominate pairwise similarity. The literature also uses Mahalanobis Gaussian and quadratic kernels. A fixed kernel baseline is often defined as the same kernel with $f:\mathbb{R}^d\to\mathbb{R}^c$0 recursive updates; in one empirical study, the “baseline” is precisely the Laplacian kernel at iteration $f:\mathbb{R}^d\to\mathbb{R}^c$1 (Shen et al., 2024, Gupta et al., 2023).

A recurring interpretation is that RFM supplies kernel machines with an explicit feature-learning step analogous, in a limited sense, to hidden-layer adaptation in neural networks. Several papers connect this to the Neural Feature Ansatz, under which AGOP is correlated with matrices of the form $f:\mathbb{R}^d\to\mathbb{R}^c$2 in trained feedforward networks (Pasternak, 1 Jul 2026).

2. Linear RFM, spectral regularization, and the IRLS connection

The most developed theory in the supplied literature concerns the linear specialization, lin-RFM, for sparse linear regression and low-rank matrix recovery. In matrix sensing or completion, with observations $f:\mathbb{R}^d\to\mathbb{R}^c$3, lin-RFM trains a minimum-Frobenius-norm linear predictor in transformed coordinates and updates the transform via AGOP. In this regime, AGOP becomes exactly a Gram matrix: if $f:\mathbb{R}^d\to\mathbb{R}^c$4, then $f:\mathbb{R}^d\to\mathbb{R}^c$5 (Radhakrishnan et al., 2024).

The lin-RFM update can be written as

$f:\mathbb{R}^d\to\mathbb{R}^c$6

and

$f:\mathbb{R}^d\to\mathbb{R}^c$7

After reparameterization by $f:\mathbb{R}^d\to\mathbb{R}^c$8, fixed points of the iteration are first-order critical points of a constrained spectral optimization problem

$f:\mathbb{R}^d\to\mathbb{R}^c$9

where $$\mathrm{AGOP}(f;\{x^{(j)}\}_{j=1}^n) = \frac1n\sum_{j=1}^n \frac{\partial f(x^{(j)})}{\partial x} \left(\frac{\partial f(x^{(j)})}{\partial x}\right)^\top \in\mathbb{R}^{d\times d}.$$0 is induced by the matrix function $$\mathrm{AGOP}(f;\{x^{(j)}\}_{j=1}^n) = \frac1n\sum_{j=1}^n \frac{\partial f(x^{(j)})}{\partial x} \left(\frac{\partial f(x^{(j)})}{\partial x}\right)^\top \in\mathbb{R}^{d\times d}.$$1 (Radhakrishnan et al., 2024).

This yields an explicit bridge to classical iterative reweighting. For $$\mathrm{AGOP}(f;\{x^{(j)}\}_{j=1}^n) = \frac1n\sum_{j=1}^n \frac{\partial f(x^{(j)})}{\partial x} \left(\frac{\partial f(x^{(j)})}{\partial x}\right)^\top \in\mathbb{R}^{d\times d}.$$2, lin-RFM is a reparameterization of IRLS with

$$\mathrm{AGOP}(f;\{x^{(j)}\}_{j=1}^n) = \frac1n\sum_{j=1}^n \frac{\partial f(x^{(j)})}{\partial x} \left(\frac{\partial f(x^{(j)})}{\partial x}\right)^\top \in\mathbb{R}^{d\times d}.$$3

At $$\mathrm{AGOP}(f;\{x^{(j)}\}_{j=1}^n) = \frac1n\sum_{j=1}^n \frac{\partial f(x^{(j)})}{\partial x} \left(\frac{\partial f(x^{(j)})}{\partial x}\right)^\top \in\mathbb{R}^{d\times d}.$$4, the induced penalty corresponds to nuclear-norm or $$\mathrm{AGOP}(f;\{x^{(j)}\}_{j=1}^n) = \frac1n\sum_{j=1}^n \frac{\partial f(x^{(j)})}{\partial x} \left(\frac{\partial f(x^{(j)})}{\partial x}\right)^\top \in\mathbb{R}^{d\times d}.$$5-type behavior; at $$\mathrm{AGOP}(f;\{x^{(j)}\}_{j=1}^n) = \frac1n\sum_{j=1}^n \frac{\partial f(x^{(j)})}{\partial x} \left(\frac{\partial f(x^{(j)})}{\partial x}\right)^\top \in\mathbb{R}^{d\times d}.$$6, the fixed-point objective becomes log-determinant-like; for $$\mathrm{AGOP}(f;\{x^{(j)}\}_{j=1}^n) = \frac1n\sum_{j=1}^n \frac{\partial f(x^{(j)})}{\partial x} \left(\frac{\partial f(x^{(j)})}{\partial x}\right)^\top \in\mathbb{R}^{d\times d}.$$7, the method becomes even more aggressively rank promoting. The paper’s central claim is that feature learning in this linear RFM sense generalizes well-studied IRLS dynamics rather than constituting a wholly separate mechanism (Radhakrishnan et al., 2024).

The same work also proves an exact Neural Feature Ansatz statement for deep linear networks under balanced initialization and continuous-time gradient flow, and uses it to motivate “deep lin-RFM,” in which multiple AGOP-inspired matrices emulate the implicit bias of deeper linear architectures. On the computational side, it gives an SVD-free implementation for the $$\mathrm{AGOP}(f;\{x^{(j)}\}_{j=1}^n) = \frac1n\sum_{j=1}^n \frac{\partial f(x^{(j)})}{\partial x} \left(\frac{\partial f(x^{(j)})}{\partial x}\right)^\top \in\mathbb{R}^{d\times d}.$$8 case and reports scaling to matrices with millions of missing entries; on a $$\mathrm{AGOP}(f;\{x^{(j)}\}_{j=1}^n) = \frac1n\sum_{j=1}^n \frac{\partial f(x^{(j)})}{\partial x} \left(\frac{\partial f(x^{(j)})}{\partial x}\right)^\top \in\mathbb{R}^{d\times d}.$$9 rank-$M_0$0 completion task with $M_0$1 observations, the implementation reached $M_0$2 test error in $M_0$3 seconds, compared with $M_0$4 seconds and $M_0$5 seconds for the cited deep linear baselines (Radhakrishnan et al., 2024).

3. Generalization phenomena: grokking, symmetry, and feature scaling

A major strand of RFM research uses the method as a controlled model of feature learning itself. On modular arithmetic, RFM exhibits grokking: every iteration performs exact kernel interpolation, so training loss is identically zero and training accuracy is $M_0$6, yet test accuracy remains near zero for many iterations before rapidly jumping to $M_0$7. In the reported $M_0$8 experiments trained for $M_0$9 iterations, test metrics stay almost unchanged for roughly $$f^{(t)}(x)=k(x,X;M_t)\alpha,\qquad \alpha = k(X,X;M_t)^{-1}y,$$0 iterations, test loss begins improving around iteration $$f^{(t)}(x)=k(x,X;M_t)\alpha,\qquad \alpha = k(X,X;M_t)^{-1}y,$$1, and perfect test accuracy follows a few iterations later. The proposed explanation is not optimization progress, but gradual feature learning via AGOP (Mallinar et al., 2024).

The learned matrices in these algebraic tasks acquire highly structured forms. For modular addition and subtraction, the final $$f^{(t)}(x)=k(x,X;M_t)\alpha,\qquad \alpha = k(X,X;M_t)^{-1}y,$$2 is reported to be block-circulant,

$$f^{(t)}(x)=k(x,X;M_t)\alpha,\qquad \alpha = k(X,X;M_t)^{-1}y,$$3

with circulant off-diagonal block $$f^{(t)}(x)=k(x,X;M_t)\alpha,\qquad \alpha = k(X,X;M_t)^{-1}y,$$4; multiplication and division reveal the same structure after discrete-log reindexing. The paper further presents theoretical evidence that such block-circulant features implement the Fourier Multiplication Algorithm, previously proposed as the generalizing solution learned by neural networks on the same tasks (Mallinar et al., 2024).

A related group-theoretic study sharpens this picture. On algebraic tasks with one-hot encoded inputs, RFM reportedly generalizes only when the training set breaks a nontrivial symmetry of the target. For modular addition with $$f^{(t)}(x)=k(x,X;M_t)\alpha,\qquad \alpha = k(X,X;M_t)^{-1}y,$$5, holding out all fixed points under a reflection subgroup yields essentially $$f^{(t)}(x)=k(x,X;M_t)\alpha,\qquad \alpha = k(X,X;M_t)^{-1}y,$$6 test accuracy even though the training fraction is around $$f^{(t)}(x)=k(x,X;M_t)\alpha,\qquad \alpha = k(X,X;M_t)^{-1}y,$$7. By contrast, starting from that non-generalizing split, moving one random point from train to test can restore $$f^{(t)}(x)=k(x,X;M_t)\alpha,\qquad \alpha = k(X,X;M_t)^{-1}y,$$8 held-out accuracy for the Gaussian kernel on addition, whereas removing symmetric reflection pairs does not help. The authors interpret the learned AGOP as encoding permutation representations of subgroup elements, so that successful RFM generalization becomes a form of symmetry recovery rather than mere interpolation (Bernal et al., 31 Mar 2026).

A third empirical phenomenon concerns ambient dimension. In a synthetic regression study, RFMs were evaluated while appending irrelevant Gaussian coordinates to data whose target depends only on the first few inputs. The reported RFM test MSE follows a decrease–increase–decrease pattern: a sharp decrease up to about $$f^{(t)}(x)=k(x,X;M_t)\alpha,\qquad \alpha = k(X,X;M_t)^{-1}y,$$9, an increase until around $M_{t+1}=\bigl[\mathrm{AGOP}(f^{(t)})\bigr]^s,$0, and then a shallower decrease from there to $M_{t+1}=\bigl[\mathrm{AGOP}(f^{(t)})\bigr]^s,$1. The same qualitative pattern appears across dataset sizes, label-noise levels, and a second target function, while the fixed Laplacian-kernel baseline exhibits earlier turning points around $M_{t+1}=\bigl[\mathrm{AGOP}(f^{(t)})\bigr]^s,$2 and $M_{t+1}=\bigl[\mathrm{AGOP}(f^{(t)})\bigr]^s,$3. The authors explicitly state that the similarity to double descent is phenomenological rather than established; no interpolation-threshold theorem is given (Gupta et al., 2023).

4. Robustness and scalability extensions

Several later works modify the basic recursion to address regimes in which standard RFMs underperform. The most explicit such variant is K-Inverse-RFM, introduced for modular arithmetic, GCD, and related mathematical tasks under label corruption, imbalance, and compact Chinese Remainder Theorem (CRT) encodings. Its recursion alters three components at once: label projection,

$M_{t+1}=\bigl[\mathrm{AGOP}(f^{(t)})\bigr]^s,$4

ridge-regularized kernel regression,

$M_{t+1}=\bigl[\mathrm{AGOP}(f^{(t)})\bigr]^s,$5

and an inverse-loss-weighted moving average of recent AGOP matrices,

$M_{t+1}=\bigl[\mathrm{AGOP}(f^{(t)})\bigr]^s,$6

The emphasized case $M_{t+1}=\bigl[\mathrm{AGOP}(f^{(t)})\bigr]^s,$7 motivates the name “K-Inverse.” The thesis argues that standard RFMs are limited by weak multiclass feature sharing, and that projecting labels through powers of the kernel matrix makes classes interact in a shared kernel-feature coordinate system (Pasternak, 1 Jul 2026).

Empirically, K-Inverse-RFM is reported to bridge an average of $M_{t+1}=\bigl[\mathrm{AGOP}(f^{(t)})\bigr]^s,$8 of the gap between standard RFM and neural networks on noisy modular addition and multiplication, to improve over standard RFM by as much as $M_{t+1}=\bigl[\mathrm{AGOP}(f^{(t)})\bigr]^s,$9 percentage points in test accuracy, and to remain comparable to feedforward neural networks up to about $s>0$0 label noise. On imbalanced data, it is reported to yield a $s>0$1 absolute jump in test accuracy over standard RFM at divisible-by-3 input weight $s>0$2, and on CRT inputs it outperforms both standard RFM and neural networks in $s>0$3 of $s>0$4 encoding-task pairs (Pasternak, 1 Jul 2026).

For large-scale tabular learning, xRFM localizes RFM through a supervised tree. Internal nodes fit lightweight RFM-style split models, compute AGOP, extract a top eigenvector $s>0$5, and split samples at the median of $s>0$6; each leaf then trains a local RFM on at most $s>0$7 samples. The paper states training complexity $s>0$8 and inference complexity $s>0$9 for fixed leaf size. On the TALENT benchmark, xRFM is reported to achieve the best performance across $s=\tfrac12$0 regression datasets and to be competitive to the best methods across $s=\tfrac12$1 classification datasets while outperforming GBDTs overall; it is also reported to outperform standard RFM and kernel ridge regression on TALENT regression with Wilcoxon $s=\tfrac12$2 (Beaglehole et al., 12 Aug 2025).

5. Interpretability and application domains

Because the learned matrix $s=\tfrac12$3 is explicit, RFM research repeatedly treats it as both a training mechanism and an interpretation object. In quantitative structure–property relationship modeling, RFM was applied to solubility prediction with MACCS keys, Morgan fingerprints, and a custom hybrid fingerprint (HF). The paper defines a local score

$s=\tfrac12$4

and a global score by dataset averaging. It reports that Matérn and Laplace kernels performed best among the tested kernel families; on AqSolDB with MACCS, the reported test RMSEs were $s=\tfrac12$5 for Matérn, $s=\tfrac12$6 for Laplace, $s=\tfrac12$7 for Rational Quadratic, and $s=\tfrac12$8 for Gaussian. The same study claims robust local and global explanations, substantial top-$s=\tfrac12$9 overlap with permutation importance and SHAP on ESOL and FreeSolv, and stability after removing redundant features from the custom fingerprint representation (Shen et al., 2024).

RFM has also been used as a probe for hidden states of frozen generative models. In MusicRFM, layerwise RFMs are trained on mean-pooled MusicGen-large hidden states for symbolic music concepts, AGOP eigendecomposition produces orthogonal concept directions, and the top eigenvector at selected layers is injected back into the residual stream during autoregressive decoding. The steering rule is additive,

$M_t$0

with optional layer weighting, time-varying schedules, and Bernoulli gating. In the main note-control result, target-note accuracy is reported to increase from $M_t$1 to $M_t$2, while text-prompt adherence remains within approximately $M_t$3 of the unsteered baseline over much of the tested range (Zhao et al., 21 Oct 2025).

A closely related activation-space application identifies refusal subspaces in aligned LLMs. An RFM-AGOP procedure with probe-informed initialization,

$M_t$4

and exponential moving-average AGOP updates is used to recover a multi-dimensional refusal subspace from harmful-versus-harmless residual activations. The paper reports that the subspace can be identified in seconds, that larger models require multiple directions rather than a single linear refusal vector, and that for Qwen 3 8B at least three directions are needed to exceed $M_t$5 attack success rate under ablation (Winninger, 2 Jul 2026).

6. Limits, unresolved questions, and acronym ambiguity

Despite the breadth of applications, the supplied literature repeatedly emphasizes incomplete theory. The feature-scaling study does not derive a mechanism for its decrease–increase–decrease curve and does not provide diagnostics such as AGOP spectra or learned-metric condition numbers; the resemblance to double descent is explicitly described as suggestive rather than established (Gupta et al., 2023). The K-Inverse-RFM thesis argues that standard AGOP is only a linear input reweighting and leaves open whether stronger AGOP variants are needed for broader corrupted-data robustness (Pasternak, 1 Jul 2026). Steering papers likewise report strong empirical utility while stating that the relationship between RFM-derived subspaces and alternatives such as refusal cones or sparse autoencoder features remains unclear (Winninger, 2 Jul 2026).

Kernel scaling remains a general practical constraint. Standard RFM repeatedly rebuilds kernel matrices and recomputes AGOP, so large-$M_t$6 workloads motivate either structural modifications such as xRFM or task-specific adaptations. This suggests that future RFM research will likely continue to divide between theory of AGOP-driven geometry learning and systems work on local, sparse, or approximate implementations (Beaglehole et al., 12 Aug 2025).

The acronym itself is also ambiguous outside this literature. In scientific computing, “RFM” can mean Random Feature Method, a collocation least-squares solver for PDEs that is unrelated to Recursive Feature Machines (2207.13380). In robotics, “RFM-SLAM” denotes Relative Feature Measurements-based SLAM, where “RFM” refers to relative feature measurements rather than feature-learning kernels (Agarwal et al., 2016). Within machine learning, however, Recursive Feature Machine has acquired a specific technical meaning: a kernel machine whose metric is recursively updated from AGOP so that representation learning is performed explicitly, analyzably, and, in many applications, interpretably.