Neural LoFi with Backward Coupling: A Spectral Theory of Cross-Layer Feature Learning

Neural LoFi with Backward Coupling: A Spectral Theory of Cross-Layer Feature Learning Yatin Dandi∗ Matteo Vilucchio∗ Hugo Tabanelli

Florent Krzakala

arXiv:2610.32849v1 [cs.LG] 23 Oct 2026

Information, Learning, and Physics Laboratory (IDEPHICS) EPFL, Lausanne, Switzerland

Abstract During joint training, backpropagation allows later layers to reshape earlier representations, a feedback mechanism absent from purely layer-wise learning. We extend Neural Low-Degree Filtering (Neural LoFi), a recently proposed layer-wise spectral surrogate for feature learning, to capture this interaction. Starting from the exact backpropagation gradients, we derive a local coupling between a neuron’s incoming weights and its corresponding next-layer column. In a controlled small-initialization regime with smooth activations and fixed surroundings, these dynamics amplify singular modes of a cross-layer correlation operator. Each growing mode pairs an earlier feature with an upper-layer response through shared neuron coefficients. This motivates a backward-coupling algorithm in which upper-layer sensitivities generate sample-dependent effective labels for earlier spectral steps, while corrected features refine upward transmission. A function-space variational principle characterizes the selected feature–response pairs by maximizing the label correlation of their product under kernel-norm constraints. For three-layer networks, we derive self-consistent continuous-time dynamics for the filtered hidden-layer features, the activation kernels, and the network output under the proposed backward coupling. In a hierarchical multi-index model, we establish a feature-recovery separation: backward coupling recovers lower-layer directions invisible to the original forward label-weighted statistic. The linear specialization recovers residual-dependent singular-mode dynamics, while paired modes provide a dynamical mechanism for linked low-rank covariance components consistent with the rainbow-network picture. Experiments with fully connected and convolutional networks on binary CIFAR-10 show improved accuracy over forward Neural LoFi at matched widths.

1

Introduction

Deep neural networks learn hierarchies of task-adapted representations through backpropagation (Rumelhart et al., 1986; LeCun et al., 2015). Yet a hierarchy in the forward computation does not imply a layer-wise order of learning. During joint training, features emerging in later layers reshape the learning signals received below, while changes in earlier representations alter what the upper layers can learn. Understanding this feedback requires going beyond how features are constructed successively. It raises a central question: How does learning in later layers change which features earlier layers can discover? This paper builds on Neural Low-Degree Filtering (Neural LoFi), introduced by Dandi et al. (2026b) as an explicit spectral construction for purely layer-wise learning. Given the current representation, each layer selects directions through low-degree correlations with the labels and ∗

Joint first authors.

1

passes them through a nonlinear random lift. Depth changes the geometry in which subsequent layers search for signal, making structure that is complex in the input accessible through lowdegree statistics of intermediate representations. However, this construction proceeds only forward: features learned above cannot refine the directions selected below. A useful upper representation may therefore become available only after the lower feature-selection step that could benefit from it has already been completed. Existing theories describe complementary aspects of feature learning. Neural-tangent limits characterize training with essentially fixed representations (Jacot et al., 2018; Lee et al., 2019); mean-field and dynamical-field frameworks describe evolving features in shallow and deep nonlinear networks (Mei et al., 2018; Chizat and Bach, 2018; Nguyen and Pham, 2020; Bordelon and Pehlevan, 2022; Hajjar et al., 2024). Feature-space mean-field models and tensor programs further describe feature learning at infinite width (Fang et al., 2021; Yang and Hu, 2021). Deep linear networks admit explicit analyses of coupled singular-mode dynamics (Saxe et al., 2014), while spectral analyses of quadratic and diagonal networks connect learned weight structure to generalization (Defilippis et al., 2026b). Backward feature correction studies how joint training with higher layers can refine lower representations (Allen-Zhu and Li, 2023). In this paper, we seek a tractable and interpretable spectral mechanism that identifies how an upper-layer response changes the "earlier" features. We extend Neural LoFi with a tractable and interpretable spectral surrogate that captures backward feature correction and, in addition, identifies a distinct effect of joint training: learning how lower-level cues selectively activate higher-level features. This couples the learning of lower representations to the learning of their transmission upward. Starting from the exact backpropagation gradients, we couple a neuron’s incoming row to the next-layer column that reads its activation. For a jointly trained block at small initialization, with smooth activations and fixed surroundings, we prove that the leading dynamics amplify singular modes of a cross-layer correlation operator. Each mode pairs an earlier feature with an upper response, with the same growing coefficients across intermediate neurons in both weight factors. The paired updates coordinate what a neuron detects with how its activation is transmitted upward: the incoming row acquires the lower feature, while the outgoing column couples it to the selected upper response. A similar mechanism has been seen to spontaneously emerge in the feed-forward blocks of modern LLMs (Geva et al., 2021; 2022) and we note additionally that with linear activations, the residual-driven row–column equations are exact and recover the singular-mode framework of Saxe et al. (2014). This mechanism motivates BackLoFi, a spectral correction of a forward LoFi network that visits layer boundaries from the top down. Upper-layer sensitivities define sample-dependent effective labels for earlier spectral steps; the corrected nonlinear activations then determine the corresponding upward transmission. The algorithm uses spectral selection and normalized corrections as a surrogate for the local dynamics. Its selection rule has a function-space interpretation: an earlier feature χ and b n [yψχ], subject to the kernel-norm constraints of their respective an upper response ψ maximize E spaces. The search for a lower feature thus depends on the response learned above as well as on the input representation and the labels. That backward feedback can make a previously invisible lower feature detectable, providing a drastic advantage to cross-layer training over layer-wise one. We illustrate this in a hierarchical multi-index model (Nichani et al., 2023; Wang et al., 2023; Fu et al., 2025; Dandi et al., 2026a; Tabanelli et al., 2026) with a visible anchor and additional hidden directions. In its Gaussian analysis, the original label-weighted second moment has no signal in the hidden block, while modulation by a sensitivity of the upper conditional predictor creates one. This illustrates a staircase learning mechanism through the backward pass. A second theoretical contribution is a closed description of coupled nonlinear feature learning in a three-layer network. For a continuous variant with fixed Gaussian lifts, continual forward filtering, 2

and readout gradient flow, we prove convergence to self-consistent feature and kernel dynamics. The evolving kernels determine subsequent feature selection, while an accumulated coupling retains earlier column corrections as lower features change. This provide an unifying description of how multi-layer kernel adaptivity Bach (2021) arise through cross-layer feature learning. It generalizes and extend the seminal work of Saxe et al. (2014) that was limited to linear networks. This also provides a dynamical perspective on the key–value organization identified in transformer feed-forward blocks by Geva et al. (2021; 2022): our local analysis shows how feature detectors and their associated output directions can be jointly selected and amplified during training. Further, this mechanism also allows to provide an explanation for the rainbow-network picture of (Guth et al., 2024): Shared neuron coefficients yield linked low-rank covariance components in adjacent weight matrices, providing a local dynamical mechanism predicting the rainbow structure. Finally, mechanistic experiments with fully connected and convolutional networks on binary CIFAR-10 show considerable improved accuracy over forward Neural LoFi at matched widths. They demonstrates that using upper-layer information to revise earlier representations is useful beyond the controlled theoretical setting, and that the predicted mechanism translate into real data scenarios. A public implementation, together with codes to reproduce the numerical illustrations, is available at https://github.com/IdePHICS/BackLoFi-Theory.

2

Neural Low-Degree Filtering with backward coupling

Backpropagation, and the guiding approximation — We first discuss the exact updates of a neuron’s incoming row and outgoing column, then use their local coupling to motivate our b n for spectral algorithm. Consider data Dn = {(xµ , yµ )}nµ=1 , with centered scalar labels, and write E p the empirical average. Let z0 : X → R 0 be a given input feature map; raw inputs are the special case z0 (x) = x. A depth-L network constructs hℓ = Wℓ zℓ−1 ,

1 zℓ = √ σ(hℓ ), pℓ

fb = ⟨aL , zL ⟩ + b.

(1)

Here Wℓ ∈ Rpℓ ×pℓ−1 and the activation acts coordinatewise. The sensitivity of the output to the preactivations, gℓ = ∂ fb/∂hℓ , obeys the backward recursion σ ′ (hL ) gL = √ ⊙ aL , pL

σ ′ (hℓ ) ⊤ gℓ = √ ⊙ (Wℓ+1 gℓ+1 ). pℓ

(2)

For the square loss, write ∆ = y − fb for the residual; gradient flow then satisfies b n [∆ gℓ z⊤ ]. Ẇℓ = E ℓ−1

(3)

Neuron j of layer ℓ owns two weight vectors: its incoming row, whose transpose we write wℓ,j (row j of Wℓ ), and its outgoing column ω ℓ+1,j (column j of Wℓ+1 ), through which the next layer reads its activation (Figure 1). By (2) and (3), their updates by gradient descent are #

"

′ (h ) ℓ,j bn ∆ σ √ ẇℓ,j = E ⟨ω ℓ+1,j , gℓ+1 ⟩zℓ−1 , pℓ

b n [∆ z gℓ+1 ]. ω̇ ℓ+1,j = E ℓ,j

(4)

Each update depends on the other vector: the outgoing column sets the sample-dependent gain driving the row, while the row’s activation determines the column update.

3

Layer-wise vs coupled training — In (Dandi et al., 2026b), the authors considered the layerwise “Neural LoFi regime” of the first equation in which the layers above are still untrained and vanishingly small. The factor ⟨ω ℓ+1,j , gℓ+1 (x)⟩, the effective readout of neuron j, then reduces to a sample-independent random constant āj , and the residual to the label. For an activation smooth near zero, σ ′ (hℓ,j ) ≃ σ ′ (0) + σ ′′ (0)hℓ,j gives  āj  ′ b n [y z ] + σ ′′ (0) E b n [y z z⊤ ] w ẇℓ,j ≃ √ σ (0) E ℓ−1 ℓ−1 ℓ−1 ℓ,j . pℓ

(5)

The rows evolve independently under the same label-weighted moments. Neural LoFi retains the first-moment direction and the leading second-moment eigenvectors by absolute eigenvalue, then applies a random lift. This leads to Alg.1 and weights of the form (9). We consider now a joint training procedure. With a learned upper representation, the effective readout becomes sample dependent. To isolate the effect of this coupling, we treat gℓ+1 and zℓ−1 as fixed. For a smooth activation with σ(0) = 0, we retain the constant derivative σ ′ (0) in the row update and the linear activation in the column update. Assuming σ ′ (0) > 0 gives σ ′ (0) e ⊤ ω ℓ+1,j , ẇℓ,j ≃ √ A pℓ ℓ

σ ′ (0) e ω̇ ℓ+1,j ≃ √ A ℓ wℓ,j , pℓ

b n [∆ gℓ+1 z⊤ ]. e ℓ := E A ℓ−1

(6)

Therefore, the same matrix acts in both directions and its singular pairs determine the growing modes. Close to initialization, under the initialization assumptions of App. B, the prediction is small, so ∆ ≃ y and the frozen coupling matrix becomes b n [y gℓ+1 z⊤ ]. Aℓ := E ℓ−1

(7)

Write si and (ui , vi ) for the singular values and pairs of Aℓ . We obtain the following approximation for the joint-evolution of ẇℓ,j and ⟨ω ℓ+1,j : Theorem 1 (Informal: amplification across adjacent layers). Under vanishing initialization and fixed surrounding layers, for small enough time, a singular mode (ui , vi , si ) of the frozen Aℓ , with non-zero initial component, has growing components: ′

√

′

√

′

√

′

√

σ (0)si t/ pℓ −σ (0)si t/ pℓ ⟨wℓ,j (t), vi ⟩ ≃ c+ + c− , ji e ji e σ (0)si t/ pℓ −σ (0)si t/ pℓ − c− , ⟨ω ℓ+1,j (t), ui ⟩ ≃ c+ ji e ji e

(8)

1 with amplitudes c± ji = 2 (⟨wℓ,j (0), vi ⟩ ± ⟨ω ℓ+1,j (0), ui ⟩) set by initialization. See App. B for details and error bounds.

Collecting the growing amplitudes of the linear approximation (6) over the neurons of layer ℓ, ⊤ ⊤ qi = (c+ ji )j , the growing parts have the rank-one directions qi vi in Wℓ and ui qi in Wℓ+1 , which we call a column spike. The same coefficients qi specify which intermediate neurons learn the feature and which columns transmit it. Neural LoFi with backward coupling — Our algorithm, BackLoFi (Algorithm 2), is motivated by this result and uses these singular modes to construct a spectral surrogate for cross-layer feature learning. On a fitted network, using labels in place of residuals is a design choice; the theorem justifies this replacement near initialization. The starting point is to use a first pass with forward Neural LoFi (Dandi et al., 2026b) (Algorithm 1), which builds the network one layer at a time from the input up. Given zℓ−1 , a filter 4

Algorithm 1 Forward Neural Low-Degree Filtering (Dandi et al., 2026b) Require: Data {(xµ , yµ )}nµ=1 , depth L, ranks kℓ , widths pℓ ; z0 = x, p0 = d. 1: for ℓ = 1, . . . , L do b n [y z ], C b n [y z z⊤ ] bℓ ← E cℓ ← E 2: m ▷ label-weighted moments ℓ−1 ℓ−1 ℓ−1 b ℓ by |λ| ] cℓ /∥m cℓ ∥ | leading kℓ − 1 eigenvectors of C Fℓ ← [ m √ p ×k ⊤ ℓ ℓ 4: zℓ ← σ(Rℓ Fℓ zℓ−1 )/ pℓ , Rℓ ∈ R Gaussian 5: end for 6: Fit the readout (aL , b) on zL ; return fb(x) = ⟨aL , zL (x)⟩ + b 3:

▷ filter ▷ random lift

Algorithm 2 Neural Low-Degree Filtering with Backward Coupling (BackLoFi) Require: Forward Neural LoFi network (Algorithm 1) Wℓ = Rℓ F⊤ ℓ , ℓ = 1, . . . , L, with readout (aL , b); spike ranks rℓ , scale β > 0, passes P . 1: for t = 1, . . . , P do ′ 2: Compute hℓ , zℓ for all ℓ by (1); gL ← σ√(hpLL ) ⊙ aL ▷ exact top sensitivity 3: for ℓ = L − 1, . . . , 1 do ▷ boundaries, top-down b n [y gℓ+1 z⊤ ] 4: (Uℓ , Vℓ ) ← leading rℓ singular pairs of Aℓ ← E ℓ−1 ℓ+1 5: sℓ ← U⊤ g ▷ scores ℓ 6: Q0ℓ ← Wℓ Vℓ ▷ spike columns at order zero ← 7: Fℓ ← [Fℓ | Vℓ Γrow ▷ rows of layer ℓ ℓ ], Rℓ ← [Rℓ | Rℓ ] 8: Recompute hℓ , zℓ ▷ features of the corrected rows b n [z (ys )⊤ ] 9: Qℓ ← E ▷ correlations with the modified labels ℓ ℓ col 10: Fℓ+1 ← [Fℓ+1 | Qℓ Γℓ ], Rℓ+1 ← [Rℓ+1 | Uℓ ] ▷ columns of layer ℓ + 1
σ ′ (hℓ ) ⊤ ℓ+1 0 ℓ b ▷ first-order signal (12) 11: g ← √pℓ ⊙ En [Wℓ+1 g ] + Qℓ sℓ 12: end for 13: Recompute zL ; refit (aL , b) 14: end for 15: return fb(x) = ⟨aL , zL (x)⟩ + b Fℓ ∈ Rpℓ−1 ×kℓ collects the kℓ directions amplified in (5): the leading label weighted moments. A Gaussian lift Rℓ ∈ Rpℓ ×kℓ then distributes them over the pℓ neurons, Wℓ = Rℓ F⊤ ℓ .

(9)

After the forward pass, a backward pass visits boundaries from the top down, then refits the readout. Boundary ℓ couples the rows of Wℓ to the columns of Wℓ+1 . It receives the sensitivity updated above and reads the representation below, which has not yet changed during the pass. First, retain the rℓ leading singular pairs of Aℓ : Vℓ = [vi ] selects lower features and Uℓ = [ui ] selects upper responses (line 4). These modes and responses are frozen during the boundary update and recomputed at the next pass. Second, install the lower directions with fresh Gaussian lift columns Rℓ← . Theorem 1 motivates coupling their upward transmission to the same intermediate neurons. The algorithm estimates this transmission from the corrected activations. After the row update, recompute zℓ and set ℓ+1 sℓ = U⊤ , ℓ g

e ℓ = y sℓ , y

⊤ b n [z y Qℓ = E ℓ e ℓ ].

(10)

Each column of Qℓ is a first moment under a modified label. Equivalently, Uℓ Q⊤ ℓ is the labelweighted column gradient projected onto the retained upper directions. Computing Qℓ after the row 5

backward feedback

j

incoming row

paired spikes

outgoing column

Figure 1: Backward coupling and paired spikes. The incoming row (blue) and outgoing column (red) meet at neuron j. Upper-layer sensitivity changes the feature learned below, whose activation changes upward transmission. Growing modes form paired rank-one spikes with shared neuron coefficients. The small matrices illustrate a mode concentrated on neuron j: an incoming row and its matching outgoing column.

step includes features absent from the forward filter and uses the corrected nonlinear activations without linearizing them. The two updates are ⊤ Wℓ ← Wℓ + Rℓ← Γrow ℓ Vℓ ,

⊤ Wℓ+1 ← Wℓ+1 + Uℓ Γcol ℓ Qℓ ,

(11)

implemented sequentially at lines 7–10. The new directions must be added perturbatively so that the modes of Aℓ , selected on the current network, remain informative after the update. Since Wℓ = Rℓ F⊤ ℓ , this can be done either by adding the small low-rank terms of (11) or, equivalently, by appending columns to the filter and the lift. The positive diagonal matrices Γrow and Γcol ℓ ℓ normalize correction amplitudes using a common scale β (Appendix F). Fresh row lifts and separate normalization modify the theorem’s paired gradient flow while preserving coupling through the activations. Finally, propagate the sensitivity to the boundary below. At the top, gL is exact. Below, the exact sensitivity gℓ applies the chain rule (2) through the corrected Wℓ+1 . A local estimate is ⊤ also appealing. Since the column update (11) adds Uℓ Γcol ℓ Qℓ to Wℓ+1 , the new term to product col ⊤ gℓ+1 is Q Γ s : of the upper sensitivity, the new columns transmit only the scores. The Wℓ+1 ℓ ℓ ℓ first-order estimate keeps this term, with its loading replaced by a linearized estimate Q0ℓ from the pre-correction rows, and replaces the pre-existing contribution by its training mean:  σ ′ (hℓ )  b ⊤ gℓ = √ ⊙ En [Wℓ+1 gℓ+1 ] + Q0ℓ sℓ , pℓ

Q0ℓ ∝ Wℓ Vℓ .

(12)

The matrices in this estimate are saved at the start of the pass; the lower gate uses the corrected activations. Each column of Q0ℓ is rescaled to match the installed column’s coordinate RMS. Both the mean substitution and the pre-correction loading are additional approximations. Appendix F specifies the normalization, the derivation, and the alternative signals. Subsequent passes repeat the sweep on the corrected network. Over the rainbow — A first consequence is that this provide a concrete mechanism for the rainbow network of Guth et al. (2024). First, a forward Neural Lofi network is a rainbow network albeit a simple one. Indeed each layer is a linear dimension reduction Fℓ followed by white random features Rℓ , its weight covariance Fℓ F⊤ ℓ has low rank, and the dependence between layers passes only through the filters, the random lifts being independent. In their scattering-network experiments, Guth et al. (2024, Fig. 10) observe amplification of selected covariance directions that largely preserves the initial neuron configuration. Theorem 1 gives a local mechanism for this behavior: its growing modes amplify fixed, initialization-dependent neuron patterns qi . Backward coupling then goes beyond this picture: Shared coefficients link the weight covariances of consecutive layers: if their covariance is Cq , the selected mode matrices Vℓ 6

and Uℓ produce Vℓ Cq Vℓ⊤ in incoming rows and Uℓ Cq U⊤ ℓ in outgoing columns. They also couple the weights themselves, which the rainbow model takes to be independent across layers up to a rotation: the column of Wℓ+1 that reads neuron j carries the coefficient qji tied to its row, and the lift of layer ℓ + 1 acquires the deterministic columns ui next to its random ones. In this picture it is joint training that paints the rainbow: the layers above choose which covariance a layer grows, the growth starts from the initialization and keeps it, and the covariances of consecutive layers come out tied to each other through the neurons they share. Coupled modes in linear networks — The coupling has an exact counterpart in the onehidden-layer linear network f (x) = W2 W1 x, studied by Saxe et al. (2014, Sections 1.3 and 2). For unnormalized linear activations, whitened inputs and squared-loss gradient flow, aligned positive b n [yx⊤ ] = P si ui v ⊤ obey modes of the input–target correlation E i i αi := ∥W1 vi ∥,

βi := ∥W2⊤ ui ∥,

α̇i = βi (si − αi βi ),

β̇i = αi (si − αi βi ).

(13)

The growing components scale as esi t while αi βi ≪ si ; their product then saturates at si . For deeper networks with aligned singular directions, each layer’s update contains the product of the mode amplitudes in all other layers. With D weight matrices and equal positive amplitudes across layers, the end-to-end mode strength mi satisfies 2−2/D

ṁi = D mi

(si − mi ).

(14)

For D = 2, this gives exponential growth at small mi , followed by saturation. For D > 2, the additional powers of mi change the early amplification law. Appendix D derives both cases. Theorem 1 identifies the local counterpart for nonlinear networks: a small adjacent row and column amplify modes of a cross-layer correlation matrix determined by the surrounding features. With scalar targets the linear correlation has rank at most one; input-dependent nonlinear sensitivities can generate multiple cross-layer modes. Backward coupling thus lets the response learned above alter feature selection below, while the corrected activations determine transmission upward. Section 3 characterizes this selection in function space and follows the resulting kernel evolution.

3

Lessons from learning with backward coupling

Backward coupling has a function-space description parallel to, but richer than, Neural LoFi: a selection principle at fixed representation, followed by an evolution of the representation kernels. A relevance–complexity trade-off for pairs — Neural LoFi selects a feature whose square corb n [yφ2 ]| (Dandi et al., 2026b). Backward coupling relates maximally with the label, max∥φ∥H =1 |E ℓ−1

instead selects a pair: let Kℓ−1 (x, x′ ) = ⟨zℓ−1 (x), zℓ−1 (x′ )⟩ and Gℓ+1 (x, x′ ) = ⟨gℓ+1 (x), gℓ+1 (x′ )⟩, g with RKHSs Hℓ−1 and Hℓ+1 . The leading singular pair of Aℓ gives ψ = ⟨u, gℓ+1 ⟩ and χ = ⟨v, zℓ−1 ⟩ solving b n [y ψ χ]. s1 (Aℓ ) = max E (15) ∥ψ∥Hg

ℓ+1

≤1, ∥χ∥Hℓ−1 ≤1

The RKHS constraints favor features that are simple to express in each layer’s geometry. Within these constraints, (15) selects a lower feature χ and an upper response ψ whose product is maximally correlated with y; neither feature needs to be predictive on its own. The paired row and column updates learn the lower cue in the same neurons that transmit the upper response, allowing the cue to modulate where that response is expressed. For example, a fur-texture cue could be coupled to 7

an upper response associated with an eye or ear detector, strengthening that response where fur is present when their conjunction is informative about the label. This selective activation refines where higher-level features are expressed, and thereby the geometry of the representation. Exact theory for Kernel dynamics with three layers — Joint training changes both the representations and the sensitivities that guide feature selection. Theorem 2 describes their coupled evolution in a three-layer network and its sequential infinite-width limit for a continuous version of backward coupling. In this version, the same selected lower perturbation is added to every first-layer filtered coordinate, forward Neural LoFi updates the second-layer features as the lower kernel changes, and gradient flow trains the final readout. The row correction, column correction and readout gradient update all use the same step size dt before taking the continuous-time limit. Appendix C gives the algorithm, initialization and assumptions. We assume that a forward Neural LoFi pass has already filtered the features at layers 1 and 2 to fixed dimensions k1 and k2 , respectively. During the subsequent evolution, these filtered features ⋆ (x) = F⊤ z (x) ∈ Rk1 and h⋆ (x) = F⊤ z (x) ∈ Rk2 . Let r ∼ N (0, I ) and η ∼ N (0, I ) are h1,t k1 k2 1,t 0,t 2,t 2,t 1,t be independent Gaussian vectors that specify random projections of the filtered coordinates in the first and second hidden layers, respectively. At finite widths p1 , p2 , the first-layer lift consists of p1 i.i.d. draws r1 , . . . , rp1 , and the second-layer lift of p2 i.i.d. draws η1 , . . . , ηp2 , sampled independently of the first collection. Each draw indexes one neuron and remains fixed while the filtered features evolve. After taking the continuous-time limit, we let p1 → ∞ and then p2 → ∞. Averages over the sampled neurons become Gaussian expectations over r and η. We use these vectors to index lower and upper neurons, respectively. To retain the column updates added during training, we define the coupling function Ct (η, r) as their total contribution from lower neuron r to upper neuron η, with C0 = 0. Its action on the current lower activations gives the correction to the upper preactivation. The preactivations h1,t , h2,t and activations z1,t , z2,t are therefore ⋆ h1,t (r, x) = r⊤ h1,t (x),

zj,t = σ(hj,t ),

⋆ h2,t (η, x) = η ⊤ h2,t (x) + Er [Ct (η, r)z1,t (r, x)].

(16)

Let w3,t (η) be the readout coefficient of upper neuron η. We further define a response function St (η, s) that records this neuron’s contribution to the output sensitivity at a preactivation s ∈ R. This response and the scalar output h3,t are given by St (η, s) = w3,t (η)σ ′ (s),

h3,t (x) = Eη [w3,t (η)z2,t (η, x)].

Evaluating the response at the current preactivation gives gt2 (η, x) = St (η, h2,t (η, x)). The feature kernels K1,t , K2,t and sensitivity kernel Gt average products at x, x′ over the same Gaussian lift row: K1,t (x, x′ ) = Er [z1,t (r, x)z1,t (r, x′ )], K2,t (x, x′ ) = Eη [z2,t (η, x)z2,t (η, x′ )],

Gt (x, x′ ) = Eη [gt2 (η, x)gt2 (η, x′ )].

The limiting evolution closes in the filtered features and scalar output together with the functions Ct (η, r) and St (η, s) defined above.

8

Theorem 2 (Informal: coupled dynamics in three layers). Under the assumptions of Appendix C, the discrete updates converge under the sequential limits dt → 0, p1 → ∞ and p2 → ∞ to the following set of self-consistent equations, defined over the training inputs. At each time t, let (ψt , χt ) denote the maximizer of (15) for the current sensitivity kernel Gt and the input kernel K0 , with value ρt > 0. The filtered features, coupling and response satisfy ⋆ ḣ⋆1,t,a = χt (a = 1, . . . , k1 ), h2,t = FwdLoFi2 (K1,t , y), 1b 2 b ∂t Ct (η, r) = E n [yχt (X)gt (η, X)]En [yψt (X)z1,t (r, X)], ρt

(17)

b n [(y − h3,t (X))z2,t (η, X)], ∂t St (η, s) = σ ′ (s)E

(18)

where FwdLoFi2 (K1,t , y) gives the forward Neural LoFi features at layer 2 from the current kernel and labels. The output satisfies b n [K2,t (x, X)(y − h3,t (X))] + Eη [g 2 (η, x)η]⊤ ḣ⋆ (x) ḣ3,t (x) = E t 2,t b n [K1,t (x, X)yψt (X)] + ψt (x)E

+ Eη,r [gt2 (η, x)Ct (η, r)ż1,t (r, x)].

(19)

Theorem 2 extends the layer-wise kernel constructions of forward Neural LoFi and AGOP-based Deep Recursive Feature Machines (Radhakrishnan et al., 2024; Beaglehole et al., 2024; Dandi et al., 2026b) to a continuous evolution in which upper and lower representations are revised together. The learned upper response changes which lower feature is selected, learning that feature then changes the upper kernel and the response itself. Thus the space in which features are selected evolves with training, in contrast with selection of hierarchical features from a prescribed family (Bach, 2008). Following this evolution requires only the current features, output, accumulated coupling Ct and response St . The equations are autonomous in these quantities, with earlier column updates retained in Ct , without the history integrals and two-time response functions of the dynamical mean-field formulation of Bordelon and Pehlevan (2022). Because the corrections are selected by (15), this evolution can be studied through the signal-tonoise ratio of correlations between layers. An upper response can make a lower feature detectable through its correlation with the modified label yψt , even when the feature had little correlation with y alone. For the current kernels, an informative leading population singular pair is reliably selected when its gap dominates the empirical operator error (Yu et al., 2015). Controlling this error along the trajectory would extend the emergence predictions of Neural LoFi (Dandi et al., 2026b) from successive layers to joint learning, using kernel concentration and spectral perturbation theory to predict when an upper response enables a new lower feature to emerge.

4

Mechanistic Numerical Experiments

We now illustrate the considerable advantages of backward feature correction with BackLofi over Lofi and layer-wise training in different situations. Synthetic setting: The Deep Staircase — We first consider a class of synthetic hierarchical models (Nichani et al., 2024; Wang et al., 2023; Fu et al., 2025; Tabanelli et al., 2026; Dandi et al., 2026b; Wortsman-Zurich et al., 2026), adapted here to isolate the role of backward information. The construction is also related to staircase models (Abbe et al., 2021; 2022; 2023; Dandi et al., 2024a;b; Troiani et al., 2025; Tabanelli et al., 2025), with one key difference: after the visible features have 9

been learned, the next useful directions become accessible only through a backward signal. We refer to this mechanism as the deep staircase. Specifically, we consider n samples (xµ , yµ ) with Gaussian inputs xµ ∼ N (0, Id ) and scalar outputs yµ ∈ R, and define the sample exponent α = log n/ log d. Let Φ2 (x) ∈ RD2 (d) , with D2 (d) = d(d + 1)/2, denote the complete normalized quadratic Hermite representation of x. The first level contains d1 = pU + pV + pW dense quadratic features, hU = AU Φ2 (x),

hV = AV Φ2 (x),

hW = AW Φ2 (x),

(20)

where the three blocks have distinct roles. The block hU is the visible anchor, while hV and hW are relevant to the target but are constructed to be hidden to the ordinary forward population statistic. The second level combines two quadratic coordinates of the visible block, Q=

⊤ B h − tr(B ) hU Q U Q √ , 2

Z=

⊤ B h − tr(B ) hU Z U Z √ , 2

S = sign(hV⊤ ChW ).

(21)

We choose BQ and BZ symmetric, normalized and orthogonal, with a sign-flip symmetry that leaves Q invariant while sending Z to −Z. The raw output is Q + ZS √ + b0 , yraw = tanh 2 



(22)

We finally standardize yraw of mean µy and variance σy2 to y = (yraw − µy )/σy with zero mean and unit variance.. Schematically, x −→ h(1) =

h

U |{z}

, hV , hW |{z}

−→ h(2) = Q(hU ) + Z(hU )S(hV , hW ) −→ y = g(h(2) ). (23)

|{z}

visible hidden hidden

The precise teacher distributions and normalizations are given in Appendix √ E. The mechanism is particularly transparent after decomposing the target. Set t± = (Q ± Z)/ 2 + b0 , the target reads a(Q, Z) =

1 tanh(t+ ) + tanh(t− ) − 2µy , 2σy

b(Q, Z) =

1 tanh(t+ ) − tanh(t− ) , 2σy

y = a(Q, Z) + b(Q, Z)S.

(24) (25)

By construction, a(Q, −Z) = a(Q, Z) and b(Q, −Z) = −b(Q, Z) Thus a carries the visible part of the hierarchy, while bS contains the hidden interaction. These parities also summarize what the ordinary forward hierarchy can learn. At the first level, the label-weighted covariance has a nonzero population component on hU , while its entire hV ⊕ hW block cancels. Forward LoFi therefore recovers the visible first-level subspace. Once hU is available, the upper Forward step acts on its quadratic representation: the first moment identifies Q, while the second-order statistic, after removing the Q direction, identifies Z. Hence the forward staircase reaches the visible hierarchy b (1) = (h b U , 0, 0) −→ h b (2) = (Q, b Z), b x −→ h

(26)

but it still cannot access the hidden hV , hW branch. The backward pass provides the missing stepping stone. Specializing the cross-layer estimator of Section 2 to the first boundary, h i A = E y g(2) (x)z0 (x)⊤ , (27) 10

1

LoFi BackLoFi GD layer-wise GD end-to-end

0.6 0.6

1

Hidden V; W overlap

LoFi BackLoFi

Visible U overlap

Test MSE

1

0.75 0.5 0.25 0

1.5

2

2.5

3

3.5

1.5

® = log n=log d

2

2.5

3

3.5

1 0.75 0.5 0.25 0 1.5

® = log n=log d

2

2.5

3

3.5

® = log n=log d

Figure 2: Backward coupling unlocks the deep staircase. Performance as a function of the sample exponent α = log n/ log d. Left: Test MSE. Center: Overlap with the visible first-level block U . Right: Overlap with the hidden blocks V, W . Both spectral methods recover the anchor, the backward correction also recovers the hidden blocks. The prediction panel additionally compares layer-wise and end-to-end gradient descent. where g(2) is the sensitivity of the fitted upper predictor. A retained left singular mode u defines the sample-dependent score s = u⊤ g(2) . For the present teacher, the informative score has a a non-zero linear component along Z. The resulting Forward–Backward–Forward mechanism can therefore be summarized as Forward

Forward

y7→ys

backward

(1)

Forward

(2)

b U −−−−−→ (Q, b b b b . b Z) b − x −−−−−→ h −−−−−→ s(x) ∋ η Zb −−−−→ h −−−−→ h B = hU , hV W − B

(28)

Here the subscript B denotes the representation obtained after the backward correction. The backward score turns the previously invisible hidden interaction into a nonzero lower-level population signal: h i h i ⊤ ⊤ E y hV hW = 0, E y s h V hW ̸= 0. (29) Moreover, this cancellation persists under any scalar transformation of the label: g : R → R, analogoushto the Statistical Query staircase settings in Troiani et al. (2025); Damian et al. (2025b), i ⊤ namely, E g(y) hV hW = 0. The corrected lower pass can therefore recover the missing joint hidden subspace, after which an ordinary Forward-LoFi pass is rerun on the corrected representation. Figure 2 exhibits the same qualitative separation for gradient descent and its spectral surrogate: joint training and BackLoFi recover hidden directions missed by their layer-wise counterparts. The corresponding sample scales make the staircase quantitative. In the regime pU , pV , pW ≍ dϵ , 0 < ϵ < 1, the visible lower subspace appears at nU ≍ d2+ϵ ,

nQ,Z ≍ d2ϵ ,

nForwardLoFi ≍ d2+ϵ ,

(30)

while, conditional on an informative backward score, recovery of the hidden subspace again requires d2+ϵ samples. The additional bottleneck is learning the score itself. Appendix E.4 derives a d−1/2 contribution to the cross-layer matrix from the quadratic feature contraction, leading to the sample scale nscore ≍ d3 , nBackLoFi ≍ max{d2+ϵ , d3 } = d3 . (31) For the numerical choice ϵ = 1/2, this predicts the two characteristic scales αU = 52 and αBack = 3. The d3 prediction is conditional on the d−1/2 informative cross-layer gap. Appendix E gives the complete derivation and the assumptions underlying these scalings.

11

Convolutional network

Fully connected network Ridge RF FwdLoFi

BackLoFi (1 ord) BackLoFi (exact) GD

20

Test Error (%)

Test Error (%)

25

15

st

15 12

RF FwdLoFi BackLoFi (1st ord)

12

BackLoFi (exact) GD

10 8 7 6

10 103

n

104

103

n

104

Figure 3: Backward coupling improves on forward Neural LoFi and approaches gradient descent. Test error vs number of training samples n on binary CIFAR-10 (animal vs. vehicle), at matched width, for a fully connected network (left) and a convnet (right). Ridge is a linear regression on the pixels, RF the random-features network of the same architecture, FwdLoFi forward Neural LoFi, BackLoFi Algorithm 2 with the first-order or the exact backward signal, and GD full-batch gradient descent on the same architecture. Curves are means over seeds and bands one std error; Spectral configurations and GD hyperparameters are selected on mean test error (see App. H). Image classification — We now test whether the coupling isolated in Section 2 accounts for what joint training learns on real data. We run Algorithm 2 on a fully connected and on a convolutional multi-layer network on binary CIFAR-10, and compare it at matched width with forward Neural LoFi, with the random-features network of the same architecture, and with full-batch gradient descent (Figure 3; settings in Appendix H). One can read Fig. 3 in three steps, each holding for both architectures. First, BackLoFi improves on forward LoFi at every training-set size and at equal width. The information that the layers above return to a boundary selects directions that the forward filter, built from the same representation, did not find. The gain is largest when data are scarce: on the convolutional network, BackLoFi trained on 5,000 images comes within 0.15 points of forward Neural LoFi trained on 50,000. The replacement of features builds up over a few backward passes, until about half of the forward columns have been replaced (see Appendix G). Second, the exact and the first-order backward signals give nearly the same curves. Of the sensitivity of the layers above, the first-order signal keeps only the scores sℓ , one number per selected mode. What the layers above pass down is therefore captured by these rℓ numbers per sample. Third, BackLoFi follows gradient descent on the same architecture. It matches it up to n = 5,000 on the fully connected network and is ahead of it on the convolutional one, by two points at n = 500. Gradient descent pulls ahead only on the largest training sets, where BackLoFi still recovers about half of the gap between forward Neural LoFi and joint training. This is consistent with joint training using more than a few coupled modes per boundary once data are abundant. The same mechanism appeared in the deep staircase above, where directions invisible to the forward statistic become visible through the response of the layers above. Over most of the sample range, then, what joint training adds to layer-wise learning is a low-rank coupling between the rows of a layer and the columns that read them. Installing a few such modes

12

per boundary, with no gradient step and using only the scores of the layers above, reproduces it.

5

Conclusion

Backward Neural LoFi extends layer-wise spectral learning by coupling lower features to upper responses: information learned above changes which directions are selected below. The resulting feature–response pairs admit a variational characterization, while their local dynamics provide a mechanism for joint amplification and structured weight covariances. For a continuous three-layer variant, we prove convergence to a self-consistent evolution of the nonlinear features, kernels, and output. The hierarchical experiment illustrates hidden-feature discovery, and the image experiments show improvements over forward learning at matched width. Controlling repeated corrections beyond the local regime and establishing finite-sample guarantees for the implemented algorithm with learned upper scores remain open. We believe these results help to open new directions into the understanding of the unreasonable efficiency of learning with deep networks (LeCun et al., 2015).

Acknowledgments We would like to thank Surya Ganguli, Bruno Loureiro, Andrew Saxe, Matthieu Wyart and Lenka Zdeborová for discussions. We acknowledge funding from the Swiss National Science Foundation grants OperaGOST (grant number 200021_200390) and DSGIANGO (grant number 225837). This work was supported by the Simons Collaboration on the Physics of Learning and Neural Computation via the Simons Foundation grant (#1257412) and by the EPFL AI Center PhD Fellowship Program 2026.

References Emmanuel Abbe, Enric Boix-Adsera, Matthew S Brennan, Guy Bresler, and Dheeraj Nagaraj. The staircase property: How hierarchical structure can guide deep learning. In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan, editors, Advances in Neural Information Processing Systems, volume 34, pages 26989–27002. Curran Associates, Inc., 2021. URL https://proceedings.neurips.cc/paper_files/paper/2021/file/e2db7186375992e 729165726762cb4c1-Paper.pdf. Emmanuel Abbe, Enric Boix Adsera, and Theodor Misiakiewicz. The merged-staircase property: a necessary and nearly sufficient condition for SGD learning of sparse functions on two-layer neural networks. In Po-Ling Loh and Maxim Raginsky, editors, Proceedings of Thirty Fifth Conference on Learning Theory, volume 178 of Proceedings of Machine Learning Research, pages 4782–4887. PMLR, 02–05 Jul 2022. Emmanuel Abbe, Enric Boix-Adserà, and Theodor Misiakiewicz. SGD learning on neural networks: leap complexity and saddle-to-saddle dynamics. In Gergely Neu and Lorenzo Rosasco, editors, Proceedings of Thirty Sixth Conference on Learning Theory, volume 195 of Proceedings of Machine Learning Research, pages 2552–2623. PMLR, 12–15 Jul 2023. Zeyuan Allen-Zhu and Yuanzhi Li. Backward feature correction: How deep learning performs deep (hierarchical) learning. In Proceedings of Thirty Sixth Conference on Learning Theory, volume 195 of Proceedings of Machine Learning Research, pages 4598–4598, 2023. URL https: //proceedings.mlr.press/v195/allen-zhu23a.html. 13

Francis Bach. The quest for adaptivity. Machine Learning Research Blog, 155, 2021. Francis R. Bach. Exploring large feature spaces with hierarchical multiple kernel learning. In Advances in Neural Information Processing Systems, volume 21, pages 105–112, 2008. URL https://papers.nips.cc/paper_files/paper/2008/hash/a4a042cf4fd6bfb47701cbc8a16 53ada-Abstract.html. Jinho Baik, Gérard Ben Arous, and Sandrine Péché. Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices. Annals of probability, 33(5):1643–1697, 2005. Lorenzo Bardone and Sebastian Goldt. Sliding down the stairs: how correlated latent variables accelerate learning with neural networks. In International Conference on Machine Learning, volume 235, pages 3024–3045, 2024. Daniel Beaglehole, Adityanarayanan Radhakrishnan, Parthe Pandit, and Mikhail Belkin. Mechanism of feature learning in convolutional neural networks. arXiv preprint arXiv:2309.00570, 2023. Daniel Beaglehole, Peter Súkeník, Marco Mondelli, and Mikhail Belkin. Average gradient outer product as a mechanism for deep neural collapse. In Advances in Neural Information Processing Systems, volume 37, 2024. URL https://proceedings.neurips.cc/paper_files/paper/202 4/hash/ec045a5ca2d8cfc528591b4c34296370-Abstract-Conference.html. Eugene Belilovsky, Michael Eickenberg, and Edouard Oyallon. Decoupled greedy learning of cnns. In International Conference on Machine Learning, pages 736–745, 2020. Gérard Ben Arous, Reza Gheissari, and Aukosh Jagannath. Online stochastic gradient descent on non-convex losses from high-dimensional inference. Journal of Machine Learning Research, 22 (106):1–51, 2021. Gérard Ben Arous, Cédric Gerbelot, and Vanessa Piccolo. Stochastic gradient descent in high dimensions for multi-spiked tensor pca. Communications on Pure and Applied Mathematics, 79 (10):2291–2369, 2026. doi: https://doi.org/10.1002/cpa.70056. URL https://onlinelibrary. wiley.com/doi/abs/10.1002/cpa.70056. Enric Boix-Adserà, Neil Rohit Mallinar, James B. Simon, and Mikhail Belkin. FACT: A firstprinciples alternative to the neural feature ansatz for how networks learn representations. In The Fourteenth International Conference on Learning Representations, 2026. URL https: //openreview.net/forum?id=j4964wtJMz. Fabrizio Boncoraglio, Vittorio Erba, Emanuele Troiani, Yizhou Xu, Florent Krzakala, and Lenka Zdeborová. Single-head attention in high dimensions: A theory of generalization, weights spectra, and scaling laws. arXiv preprint arXiv:2509.24914, 2025. Tony Bonnaire, Giulio Biroli, and Chiara Cammarota. The role of the time-dependent hessian in high-dimensional optimization. Journal of Statistical Mechanics: Theory and Experiment, 2025 (8):083401, 2025. Blake Bordelon and Cengiz Pehlevan. Self-consistent dynamical field theory of kernel evolution in wide neural networks. In Advances in Neural Information Processing Systems, volume 35, pages 32240–32256, 2022. URL https://arxiv.org/abs/2205.09653. Joan Bruna. The multiscale single-index model: A stylized model for hierarchical feature learning. arXiv preprint arXiv:2607.03347, 2026. URL https://arxiv.org/abs/2607.03347. 14

Francesco Cagnetta, Leonardo Petrini, Umberto M. Tomasini, Alessandro Favero, and Matthieu Wyart. How deep neural networks learn compositional data: The random hierarchy model. Physical Review X, 14(3):031001, 2024. doi: 10.1103/PhysRevX.14.031001. URL https: //arxiv.org/abs/2307.02129. Taehun Cha, Daniel Beaglehole, Adityanarayanan Radhakrishnan, and Donghun Lee. The weight gram matrix captures sequential feature linearization in deep networks. arXiv preprint arXiv:2605.06258, 2026. Lénaïc Chizat and Francis Bach. On the global convergence of gradient descent for over-parameterized models using optimal transport. In Advances in Neural Information Processing Systems, volume 31, 2018. Alex Damian, Loucas Pillaud-Vivien, Jason Lee, and Joan Bruna. Computational-statistical gaps in gaussian single-index models. In The Thirty Seventh Annual Conference on Learning Theory, pages 1262–1262. PMLR, 2024. Alex Damian, Jason Lee, and Joan Bruna. The generative leap: Tight sample complexity for efficiently learning gaussian multi-index models. In D. Belgrave, C. Zhang, H. Lin, R. Pascanu, P. Koniusz, M. Ghassemi, and N. Chen, editors, Advances in Neural Information Processing Systems, volume 38, pages 28276–28311. Curran Associates, Inc., 2025a. doi: 10.52202/085713-0 953. Alex Damian, Jason D. Lee, and Joan Bruna. The generative leap: Tight sample complexity for efficiently learning gaussian multi-index models. In The Thirty-ninth Annual Conference on Neural Information Processing Systems, 2025b. URL https://openreview.net/forum?id=5X6PL4906S. Yatin Dandi, Florent Krzakala, Bruno Loureiro, Luca Pesce, and Ludovic Stephan. How two-layer neural networks learn, one (giant) step at a time. Journal of Machine Learning Research, 25(349): 1–65, 2024a. URL http://jmlr.org/papers/v25/23-1543.html. Yatin Dandi, Emanuele Troiani, Luca Arnaboldi, Luca Pesce, Lenka Zdeborová, and Florent Krzakala. The benefits of reusing batches for gradient descent in two-layer networks: Breaking the curse of information and leap exponents. In Ruslan Salakhutdinov, Zico Kolter, Katherine Heller, Adrian Weller, Nuria Oliver, Jonathan Scarlett, and Felix Berkenkamp, editors, Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pages 9991–10016. PMLR, 21–27 Jul 2024b. Yatin Dandi, Luca Pesce, Lenka Zdeborová, and Florent Krzakala. The computational advantage of depth in learning high-dimensional hierarchical targets. Advances in Neural Information Processing Systems, 38:160705–160764, 2026a. Yatin Dandi, Matteo Vilucchio, Luca Arnaboldi, Hugo Tabanelli, and Florent Krzakala. Deep learning as Neural Low-Degree Filtering: A spectral theory of hierarchical feature learning. arXiv preprint arXiv:2605.13612, 2026b. URL https://arxiv.org/abs/2605.13612. Leonardo Defilippis, Florent Krzakala, Bruno Loureiro, and Antoine Maillard. Optimal scaling laws in learning hierarchical multi-index models. arXiv preprint arXiv:2602.05846, 2026a. Leonardo Defilippis, Yizhou Xu, Julius Girardin, Vittorio Erba, Emanuele Troiani, Lenka Zdeborová, Bruno Loureiro, and Florent Krzakala. Scaling laws and spectra of shallow neural networks in the feature learning regime. In The Fourteenth International Conference on Learning Representations, 2026b. URL https://openreview.net/forum?id=Q3yLIIkt7z. 15

Samuel F Edwards and Raymund C Jones. The eigenvalue spectrum of a large symmetric random matrix. Journal of Physics A: Mathematical and General, 9(10):1595–1603, 1976. Vittorio Erba, Emanuele Troiani, Lenka Zdeborová, and Florent Krzakala. The nuclear route: Sharp asymptotics of erm in overparameterized quadratic networks. In The Thirty-ninth Annual Conference on Neural Information Processing Systems, 2025. Cong Fang, Jason Lee, Pengkun Yang, and Tong Zhang. Modeling from features: a mean-field framework for over-parameterized deep neural networks. In Proceedings of Thirty Fourth Conference on Learning Theory, volume 134 of Proceedings of Machine Learning Research, pages 1887–1936, 2021. URL https://proceedings.mlr.press/v134/fang21a.html. Hengyu Fu, Zihao Wang, Eshaan Nichani, and Jason D. Lee. Learning hierarchical polynomials of multiple nonlinear features with three-layer networks. In Proceedings of the International Conference on Learning Representations (ICLR), 2025. Wulfram Gerstner, Marco Lehmann, Vasiliki Liakoni, Dane Corneil, and Johanni Brea. Eligibility traces and plasticity on behavioral time scales: experimental support of neohebbian three-factor learning rules. Frontiers in neural circuits, 12:53, 2018. Mor Geva, Roei Schuster, Jonathan Berant, and Omer Levy. Transformer feed-forward layers are key-value memories. In Proceedings of the 2021 conference on empirical methods in natural language processing, pages 5484–5495, 2021. Mor Geva, Avi Caciularu, Kevin Wang, and Yoav Goldberg. Transformer feed-forward layers build predictions by promoting concepts in the vocabulary space. In Proceedings of the 2022 conference on empirical methods in natural language processing, pages 30–45, 2022. Florentin Guth, Brice Ménard, Gaspar Rochette, and Stéphane Mallat. A rainbow in deep network black boxes. Journal of Machine Learning Research, 25(350):1–59, 2024. URL https://jmlr.o rg/papers/v25/23-1573.html. Karl Hajjar, Lénaïc Chizat, and Christophe Giraud. Training integrable parameterizations of deep neural networks in the infinite-width limit. Journal of Machine Learning Research, 25(196):1–130, 2024. URL https://www.jmlr.org/papers/v25/21-1260.html. Bernd Illing, Wulfram Gerstner, and Johanni Brea. Biologically plausible deep learning—but how far can we go with shallow networks? Neural Networks, 118:90–101, 2019. Bernd Illing, Jean Ventura, Guillaume Bellec, and Wulfram Gerstner. Local plasticity rules can learn deep representations using self-supervised contrastive predictions. In Advances in Neural Information Processing Systems, volume 34, pages 30365–30379, 2021. Arthur Jacot, Franck Gabriel, and Clément Hongler. Neural tangent kernel: Convergence and generalization in neural networks. In Advances in Neural Information Processing Systems, volume 31, 2018. Daniel Kunin, Giovanni Luca Marchetti, Feng Chen, Dhruva Karkada, James B. Simon, Michael R. DeWeese, Surya Ganguli, and Nina Miolane. Alternating gradient flows: A theory of feature learning in two-layer neural networks. In Advances in Neural Information Processing Systems, 2025.

16

Julien Launay, Iacopo Poli, François Boniface, and Florent Krzakala. Direct feedback alignment scales to modern deep learning tasks and architectures. In Advances in Neural Information Processing Systems, volume 33, pages 9346–9360, 2020. Yann LeCun, Yoshua Bengio, and Geoffrey Hinton. Deep learning. Nature, 521(7553):436–444, 2015. doi: 10.1038/nature14539. Dong-Hyun Lee, Saizheng Zhang, Asja Fischer, and Yoshua Bengio. Difference target propagation. In Machine Learning and Knowledge Discovery in Databases, pages 498–515. Springer, 2015. Jaehoon Lee, Lechao Xiao, Samuel S. Schoenholz, Yasaman Bahri, Roman Novak, Jascha SohlDickstein, and Jeffrey Pennington. Wide neural networks of any depth evolve as linear models under gradient descent. In Advances in Neural Information Processing Systems, volume 32, 2019. Timothy P. Lillicrap, Daniel Cownden, Douglas B. Tweed, and Colin J. Akerman. Random synaptic feedback weights support error backpropagation for deep learning. Nature Communications, 7: 13276, 2016. Yue M Lu and Gen Li. Phase transitions of spectral initialization for high-dimensional non-convex estimation. Information and Inference: A Journal of the IMA, 9(3):507–541, 2020. Antoine Maillard, Bruno Loureiro, Florent Krzakala, and Lenka Zdeborová. Phase retrieval in high dimensions: Statistical and computational phase transitions. Advances in Neural Information Processing Systems, 33:11071–11082, 2020. Giovanni Luca Marchetti, Daniel Kunin, Adele Myers, Francisco Acosta, and Nina Miolane. Sequential group composition: A window into the mechanics of deep learning. In International Conference on Machine Learning, 2026. Charles H Martin and Michael W Mahoney. Implicit self-regularization in deep neural networks: Evidence from random matrix theory and implications for learning. Journal of Machine Learning Research, 22(165):1–73, 2021. Song Mei, Andrea Montanari, and Phan-Minh Nguyen. A mean field view of the landscape of two-layer neural networks. Proceedings of the National Academy of Sciences, 115(33):E7665–E7671, 2018. doi: 10.1073/pnas.1806579115. Marco Mondelli and Andrea Montanari. Fundamental limits of weak recovery with applications to phase retrieval. In Conference On Learning Theory, pages 1445–1450. PMLR, 2018. Andrea Montanari and Zihao Wang. Phase transitions for feature learning in neural networks. arXiv preprint arXiv:2602.01434, 2026. Phan-Minh Nguyen and Huy Tuan Pham. A rigorous framework for the mean field limit of multilayer neural networks, 2020. URL https://arxiv.org/abs/2001.11443. Eshaan Nichani, Alex Damian, and Jason D. Lee. Provable guarantees for nonlinear feature learning in three-layer neural networks. In Advances in Neural Information Processing Systems, volume 36, pages 10828–10875, 2023. Eshaan Nichani, Alex Damian, and Jason D Lee. Provable guarantees for nonlinear feature learning in three-layer neural networks. Advances in Neural Information Processing Systems, 36, 2024.

17

Arild Nøkland. Direct feedback alignment provides learning in deep neural networks. In Advances in Neural Information Processing Systems, volume 29, 2016. Arild Nøkland and Lars Hiller Eidnes. Training neural networks with local error signals. In International Conference on Machine Learning, pages 4839–4850, 2019. Erkki Oja. Simplified neuron model as a principal component analyzer. Journal of mathematical biology, 15(3):267–273, 1982. Adityanarayanan Radhakrishnan, Daniel Beaglehole, Parthe Pandit, and Mikhail Belkin. Mechanism for feature learning in neural networks and backpropagation-free machine learning models. Science, 383(6690):1461–1467, 2024. doi: 10.1126/science.adi5639. Adityanarayanan Radhakrishnan, Mikhail Belkin, and Dmitriy Drusvyatskiy. Linear recursive feature machines provably recover low-rank matrices. Proceedings of the National Academy of Sciences, 122(13):e2411325122, 2025. doi: 10.1073/pnas.2411325122. Yunwei Ren, Yatin Dandi, Florent Krzakala, and Jason D. Lee. Provable learning of random hierarchy models and hierarchical shallow-to-deep chaining. In Proceedings of Thirty Ninth Conference on Learning Theory, volume 336 of Proceedings of Machine Learning Research, pages 5556–5597. PMLR, 2026a. URL https://proceedings.mlr.press/v336/ren26a.html. Yunwei Ren, Zihao Wang, and Jason D. Lee. Algorithmic separation between constant-depth and logarithmic-depth neural networks. arXiv preprint arXiv:2607.25200, 2026b. URL https: //arxiv.org/abs/2607.25200. Valentina Ros, Gérard Ben Arous, Giulio Biroli, and Chiara Cammarota. Complex energy landscapes in spiked-tensor and simple glassy models: Ruggedness, arrangements of local minima, and phase transitions. Physical Review X, 9(1):011003, 2019. David E. Rumelhart, Geoffrey E. Hinton, and Ronald J. Williams. Learning representations by back-propagating errors. Nature, 323(6088):533–536, 1986. doi: 10.1038/323533a0. Terence D Sanger. Optimal unsupervised learning in a single-layer linear feedforward neural network. Neural networks, 2(6):459–473, 1989. Stefano Sarao Mannelli, Giulio Biroli, Chiara Cammarota, Florent Krzakala, and Lenka Zdeborová. Who is afraid of big bad minima? analysis of gradient-flow in spiked matrix-tensor models. Advances in neural information processing systems, 32, 2019. Stefano Sarao Mannelli, Giulio Biroli, Chiara Cammarota, Florent Krzakala, Pierfrancesco Urbani, and Lenka Zdeborová. Marvels and pitfalls of the langevin algorithm in noisy high-dimensional inference. Physical Review X, 10(1):011057, 2020. Andrew M. Saxe, James L. McClelland, and Surya Ganguli. Exact solutions to the nonlinear dynamics of learning in deep linear neural networks. arXiv:1312.6120, 2014. URL https: //arxiv.org/abs/1312.6120. Hugo Tabanelli, Pierre Mergny, Lenka Zdeborová, and Florent Krzakala. Computational thresholds in multi-modal learning via the spiked matrix-tensor model. arXiv preprint arXiv:2506.02664, 2025.

18

Hugo Tabanelli, Yatin Dandi, Luca Pesce, and Florent Krzakala. Deep learning of compositional targets with hierarchical spectral methods. In Proceedings of the 43rd International Conference on Machine Learning, 2026. URL https://icml.cc/virtual/2026/poster/61650. Emanuele Troiani, Yatin Dandi, Leonardo Defilippis, Lenka Zdeborová, Bruno Loureiro, and Florent Krzakala. Fundamental computational limits of weak learnability in high-dimensional multi-index models. In Yingzhen Li, Stephan Mandt, Shipra Agrawal, and Emtiyaz Khan, editors, Proceedings of The 28th International Conference on Artificial Intelligence and Statistics, volume 258 of Proceedings of Machine Learning Research, pages 2467–2475. PMLR, 03–05 May 2025. Roman Vershynin. High-Dimensional Probability: An Introduction with Applications in Data Science. Cambridge University Press, 2018. URL https://webapps.math.uci.edu/~rvershyn/papers /HDP-book/HDP-1.pdf. Cédric Villani. Topics in Optimal Transportation, volume 58 of Graduate Studies in Mathematics. American Mathematical Society, 2003. Matteo Vilucchio, Yatin Dandi, Matéo Pirio Rossignol, Cédric Gerbelot, and Florent Krzakala. Asymptotics of non-convex generalized linear models in high-dimensions: A proof of the replica formula. arXiv preprint arXiv:2502.20003, 2025. Zihao Wang, Eshaan Nichani, and Jason D Lee. Learning hierarchical polynomials with three-layer neural networks. arXiv preprint arXiv:2311.13774, 2023. Per-Åke Wedin. Perturbation bounds in connection with singular value decomposition. BIT Numerical Mathematics, 12:99–111, 1972. Garrett G. Wen, Hong Hu, Yue M. Lu, Zhou Fan, and Theodor Misiakiewicz. When does Gaussian equivalence fail and how to fix it: Non-universal behavior of random features with quadratic scaling. arXiv preprint arXiv:2512.03325, 2025. URL https://arxiv.org/abs/2512.03325. Arie Wortsman-Zurich, Hugo Tabanelli, Yatin Dandi, Florent Krzakala, and Bruno Loureiro. Scaling laws from sequential feature recovery: A solvable hierarchical model, 2026. URL https: //arxiv.org/abs/2605.14567. Greg Yang and Edward J. Hu. Tensor Programs IV: Feature learning in infinite-width neural networks. In Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pages 11727–11737, 2021. URL https: //proceedings.mlr.press/v139/yang21c.html. Yi Yu, Tengyao Wang, and Richard J. Samworth. A useful variant of the Davis–Kahan theorem for statisticians. Biometrika, 102(2):315–323, 2015. doi: 10.1093/biomet/asv008. URL https://arxiv.org/abs/1405.0680. Bohan Zhang, Zihao Wang, Hengyu Fu, and Jason D. Lee. Neural networks learn generic multi-index models near information-theoretic limit. arXiv preprint arXiv:2511.15120, 2025.

19

A

Further related works

We collect here the lines of work that the main text touches only briefly, organized by their relation to backward coupling. Feature correction by the layers above — The closest antecedent is the backward feature correction of Allen-Zhu and Li (2023), where jointly trained upper layers refine the representation learned below. Backward coupling isolates a narrower and more explicit effect: an upper response, entering as a modified label, changes which directions the lower spectral step can see at all. The same reversal of the forward order appears in a different guise in the alternating-gradient-flow analyses of Kunin et al. (2025), where features of a two-layer network activate one after the other through plateaus and sharp drops, and in the sequential learning of compositional tasks of Marchetti et al. (2026). Multi-pass learning and staircases — A single low-degree statistic of the input misses targets with information or leap exponent above two (Ben Arous et al., 2021; 2026; Damian et al., 2024; Troiani et al., 2025). Repeated passes over the data change this: reusing batches transforms the effective label and lets gradient descent escape these constraints (Dandi et al., 2024b), and staircasetype targets are learned one correlation at a time, each step being unlocked by the previous one (Abbe et al., 2021; 2023; Bardone and Goldt, 2024; Tabanelli et al., 2025). The generative-exponent picture (Damian et al., 2025a) quantifies what label transformations can achieve. BackLoFi sits in this family with a specific choice of transformation, the sample-wise response of the fitted layers above, and a specific place where it acts, the spectral step of an earlier layer. Hierarchical learning and depth — Dandi et al. (2026a) study Gaussian hierarchical singleand multi-index targets, showing how successive feature learning reduces the effective dimension and yields a sample-complexity advantage over shallow networks. In a related compositional setting, Tabanelli et al. (2026) recover latent features with layer-wise spectral estimators and establish sample-complexity gains over shallow kernel methods. Wortsman-Zurich et al. (2026) derive recovery thresholds for individual features with power-law strengths: as the sample size increases, successively weaker features become detectable, producing power-law decay of the prediction error. For a multiscale single-index model, Bruna (2026) proves recovery by online SGD on a correlation objective with all layers evolving on the same timescale. For discrete compositional data, Cagnetta et al. (2024) introduce the Random Hierarchy Model and relate the emergence of representations invariant to equivalent feature combinations to the detectability of correlations between input patches and labels. Ren et al. (2026a) prove efficient learning of these hierarchies under nondegeneracy assumptions by training a deep convolutional network layer by layer. For Boolean targets with hierarchical Fourier structure, Ren et al. (2026b) give a layer-wise coordinate-descent algorithm for logarithmic-depth networks and prove an approximation separation from constant-depth networks with polynomial width, regular activations, and controlled spectral norms. Hessian and spectral learning — Label-weighted second-order operators drive spectral initializations and recovery in single- and multi-index estimation (Lu and Li, 2020; Mondelli and Montanari, 2018; Maillard et al., 2020; Vilucchio et al., 2025), and describe the early dynamics of shallow networks through Hessian-like objects (Bonnaire et al., 2025; Zhang et al., 2025; Montanari and Wang, 2026; Defilippis et al., 2026a); the selection of spectral outliers against a bulk follows the

20

BBP/EJ mechanism (Baik et al., 2005; Edwards and Jones, 1976), as in spiked matrix–tensor models (Ros et al., 2019; Sarao Mannelli et al., 2019; 2020). Neural LoFi (Dandi et al., 2026b) applies one such operator per layer. The cross-layer matrix Aℓ of this paper correlates two representations through the labels; its leading singular pairs select the coupled features. Spectral structure of trained weights — Trained networks show low-rank, non-random structure in their weights: spikes over a bulk in quadratic and diagonal networks (Defilippis et al., 2026b), nuclear-norm-type capacity control in overparameterized quadratic networks (Erba et al., 2025), similar spectra in attention layers (Boncoraglio et al., 2025), and structured weight spectra across architectures (Martin and Mahoney, 2021). The rainbow model (Guth et al., 2024) organizes these observations into layers of random features with learned low-rank covariances, assumed independent across layers up to a rotation. Section 3 shows what backward coupling adds: linked covariances and a neuron-by-neuron dependence between the rows of one layer and the columns of the next. Feature learning from gradients of a fitted predictor — Recursive feature machines update the feature metric with the average gradient outer product (AGOP) of a fitted predictor (Radhakrishnan et al., 2024; 2025; Beaglehole et al., 2023; Boix-Adserà et al., 2026; Cha et al., 2026). These methods also feed information from the output back into the features, through the gradient of a predictor that must first be fitted. In backward coupling, the upper sensitivity gℓ+1 plays a similar role: together with y, it weights a first moment of the fixed lower representation, producing the cross-layer matrix Aℓ . Local and backpropagation-free training — Layer-wise and local training rules avoid an exact backward pass: feedback alignment and its direct variants replace the transposed weights by random matrices (Lillicrap et al., 2016; Nøkland, 2016; Launay et al., 2020), local losses train layers separately (Nøkland and Eidnes, 2019; Belilovsky et al., 2020), target propagation passes layer-wise targets (Lee et al., 2015), and label-modulated Hebbian rules approximate spectral operators (Oja, 1982; Sanger, 1989; Gerstner et al., 2018; Illing et al., 2019; 2021). The backward signal of BackLoFi is itself low-rank and approximate, carried by a few scores per sample, in the spirit of these rules.

B

Coupled mode amplification: proof of Theorem 1

We now give the full formal statement and assumptions for Theorem 1. Theorem B.1 specifies the initialization scales and training interval, and bounds the error between gradient flow and the coupled linear dynamics used in the main text. It also shows that the incoming row and matching outgoing column share the same exponentially growing term when projected onto paired singular vectors.

B.1

Assumptions and formal statement

We use the network notation of Section 2. Fix a neuron j in layer ℓ < L, with incoming row wℓ,j and matching outgoing column ω ℓ+1,j . We train this pair jointly and keep all other weights frozen. With x = zℓ−1 , the next preactivation is 1 hℓ+1 (x) = hℓ+1,0 (x) + √ ω ℓ+1,j σ(⟨wℓ,j , x⟩). pℓ Here hℓ+1,0 is the fixed contribution of all neurons in layer ℓ except j. 21

(32)

Assumption 1 (Joint training of one pair in a fixed network). The dataset, inputs x, preactivation contributions hℓ+1,0 (x), depth and widths are fixed. The selected pair (wℓ,j , ω ℓ+1,j ) follows equalb n [(y − f )2 ]; all other weights remain frozen during the entire rate Euclidean gradient flow for 12 E training interval. The subsequent layers have no biases and the output intercept is zero. Assumption 2 (Activation). The activation is C 3 , with bounded derivatives on the relevant bounded preactivation sets, σ(0) = 0, and σ ′ (0) > 0. Assumption 3 (Layer-wise vanishing initialization). Initialize wℓ,j (0) = αℓ wℓ,j and ω ℓ+1,j (0) = αℓ+1 ω ℓ+1,j . The frozen subsequent weights are Wk = αk Wk for k = ℓ + 2, . . . , L, and aL = αL+1 a0 . The normalized initial vectors wℓ,j , ω ℓ+1,j , matrices Wk and readout a0 are fixed. The scales are positive and satisfy αℓ −→ 0, αk+1 = o(αk ) (k = ℓ, . . . , L). (33) This is the scale hierarchy of Dandi et al. (2026b, Assumption 1), applied to the selected pair and the frozen weights above it. To obtain the linear approximation, evaluate the upper sensitivity at hℓ+1,0 (x), with the selected neuron’s contribution removed. This gives the fixed sensitivity and correlation matrix g0 (x) =

∂f , ∂hℓ+1 hℓ+1 =hℓ+1,0 (x)

b n [yg0 x⊤ ]. Aℓ = E

(34)

√ The activation derivative in the row and column updates contributes the factor σ ′ (0)/ pℓ . The largest growth rate of their linear approximation is therefore σ ′ (0) λ = √ ∥Aℓ ∥2 . pℓ Both Aℓ and λ depend on the initialization scales; they remain fixed during training. The sensitivity contains one weight factor from each subsequent layer, including the readout, so Q its size is bounded by a constant times L+1 k=ℓ+2 αk . We assume that averaging against the labels and inputs leaves a correlation of this same order: λ≥c

L+1 Y

αk

for some fixed c > 0 and all sufficiently small scales.

k=ℓ+2

This condition excludes cancellation of the leading correlation that drives the coupled growth.

22

(35)

Theorem B.1 (Formal version of Theorem 1). Under Assumptions 1–3 and (35), let the linear approximation solve σ ′ (0) σ ′ (0) lin lin lin ẇℓ,j = √ A⊤ ω̇ lin Aℓ wℓ,j ℓ ω ℓ+1,j , ℓ+1,j = √ pℓ pℓ from the same initialization. For any fixed M > 0 and sufficiently small initialization scales, the following bounds hold uniformly for 0 ≤ t ≤ T with λT ≤ M : ∥wℓ,j (t)∥ ≤ CM αℓ ,

lin ∥wℓ,j (t) − wℓ,j (t)∥ ≤ CM αℓ2 ,

∥ω ℓ+1,j (t)∥ ≤ CM αℓ ,

2 ∥ω ℓ+1,j (t) − ω lin ℓ+1,j (t)∥ ≤ CM αℓ .

(36)

Here CM is independent of the initialization scales and time. Hence, the weights move within their initialization scale αℓ , along the linearized trajectories with the error bounded as the smaller ± 1 scale O(αℓ2 ). For each singular triple, Aℓ vi = si ui and A⊤ ℓ ui = si vi , write cji = 2 (⟨wℓ,j (0), vi ⟩ ± ⟨ω ℓ+1,j (0), ui ⟩). Then ′

√

′

√

′

√

′

√

σ (0)si t/ pℓ −σ (0)si t/ pℓ ⟨wℓ,j (t), vi ⟩ = c+ + c− + ϵw ji (t), ji e ji e σ (0)si t/ pℓ −σ (0)si t/ pℓ ⟨ω ℓ+1,j (t), ui ⟩ = c+ − c− + ϵωji (t), ji e ji e

(37)

with both errors satisfying the approximation-error bounds in (36). Remark (the decaying component). The decaying component comes from the difference lin (t), v ⟩ between the initial row and column projections. For the linearized trajectories, let a(t) = ⟨wℓ,j i and b(t) = ⟨ω lin (t), u ⟩. Adding and subtracting their evolution equations gives i ℓ+1,j d σ ′ (0)si (a + b) = √ (a + b), dt pℓ

d σ ′ (0)si (a − b) = − √ (a − b). dt pℓ

− Thus, for si > 0, their sum grows while their difference decays. The coefficients c+ ji and cji are half the initial sum and difference, respectively. The decaying term is needed to recover both initial projections; it vanishes exactly when they are equal. The two weights share the growing amplitude c+ ji , as stated in Theorem 1. Since the outgoing ± column is initialized at scale αℓ+1 = o(αℓ ), cji = 12 αℓ ⟨wℓ,j , vi ⟩ + o(αℓ ). At t = 0, their leading contributions add in the row projection and cancel in the column projection, leaving its smaller initial value. The column starts at scale αℓ+1 = o(αℓ ) but can grow to scale αℓ during this interval. On the shorter interval λT ≤ M αℓ+1 /αℓ , it remains at its own initialization scale: ∥ω ℓ+1,j (t)∥ ≤ CM αℓ+1 . Both approximation errors are then bounded by CM αℓ αℓ+1 , as shown in the proof below.

Remark (gradient descent). We use gradient flow for convenience. For full-batch gradient √ descent with step size ∆t, the linear mode factors after m steps are (1 ± ∆t σ ′ (0)si / pℓ )m . The same remainder estimates and discrete Grönwall give the corresponding bounds relative to these linear iterates. If ∆t λ < 1, replacing the factors by their exponentials introduces relative error O(∆t λ2 T ). This error vanishes for fixed ∆t, since λ2 T ≤ M λ → 0.

B.2

Proof

Proof sketch. Lemma 1 bounds the prediction and the change in upper sensitivity as the trained pair moves. Using these bounds, Proposition 1 controls the error in replacing the gradient updates by the fixed linear system of Theorem B.1. Lemma 2 solves this system and bounds how much it 23

amplifies an error introduced at an earlier time. In the final proof, we integrate the remainder with this amplification factor and show that the weights remain small throughout [0, T ]. All constants C below are uniform in the initialization scales and time; they may depend on the fixed data, architecture, normalized seeds, derivative bounds and the constant in (35). Constants denoted CM may additionally depend on M . Write ρ = (∥wℓ,j ∥2 + ∥ω ℓ+1,j ∥2 )1/2 . We establish the local estimates for ρ ≤ r∗ , with a fixed small radius r∗ ≤ 1, before showing that the trajectory stays in this ball. Lemma 1 (Stability of the upper sensitivity). Under the assumptions of Theorem B.1, for sufficiently small scales and ρ ≤ r∗ , uniformly over the training inputs, ∥gℓ+1 − g0 ∥ ≤ Cλρ2 ,

∥gℓ+1 ∥ ≤ Cλ,

|f | ≤ Cλ.

(38)

The upper preactivation changes by at most Cρ2 , and λ = o(αℓ ). Proof. By σ(0) = 0 and the bounded inputs, ∥hℓ+1 (x) − hℓ+1,0 (x)∥ ≤ C∥ω ℓ+1,j ∥ ∥wℓ,j ∥ ≤ Cρ2 .

(39)

Thus the upper preactivations stay in a fixed bounded neighborhood when ρ ≤ r∗ . Since σ(0) = 0 and σ ′ is bounded on the relevant intervals, ∥zk ∥ ≤ Cαk ∥zk−1 ∥ for each frozen layer k = ℓ + 2, . . . , L. Q Together with |f | ≤ CαL+1 ∥zL ∥ and the boundedness of zℓ+1 , this gives |f | ≤ C L+1 k=ℓ+2 αk . Every chain-rule term for the first and second derivatives with respect to hℓ+1 contains each subsequent weight matrix and the readout at least once. The factors σ ′ and σ ′′ are bounded, and any additional weight factors are bounded as the scales vanish. Therefore |f | +

L+1 Y ∂f ∂2f + ≤ C αk ≤ Cλ, 2 ∂hℓ+1 ∂hℓ+1 2 k=ℓ+2

where the last inequality uses (35). The mean-value theorem and (39) give the sensitivity difference in (38). The same derivative bound at hℓ+1,0 (x), together with the definition of Aℓ , gives λ≤C

L+1 Y

αk = o(αℓ ),

k=ℓ+2

by the initialization hierarchy. This estimate controls the errors from freezing the upper sensitivity and replacing the residual by the label in the two gradient updates. Proposition 1 (Linear coupling of a row and its matching column). Under the assumptions of Lemma 1, for sufficiently small scales and ρ ≤ r∗ , σ ′ (0) σ ′ (0) ẇℓ,j − √ A⊤ ω , ω̇ − √ Aℓ wℓ,j ℓ+1,j ℓ+1,j pℓ ℓ pℓ

!

≤ Cλ(ρ2 + αℓ ρ).

(40)

Proof. The exact gradient updates in (4) are 1 b ′ ⊤ ℓ+1 ẇℓ,j = √ E )x], n [(y − f )σ (hℓ,j )(ω ℓ+1,j g pℓ

24

1 b ℓ+1 ω̇ ℓ+1,j = √ E σ(hℓ,j )]. n [(y − f )g pℓ

(41)

Since |hℓ,j | ≤ Cρ, Taylor’s theorem gives σ ′ (hℓ,j ) = σ ′ (0) + O(ρ),

σ(hℓ,j ) = σ ′ (0)hℓ,j + O(ρ2 ).

Replacing the activation and its derivative by these leading terms, the sensitivity by g0 , and the residual by y gives the linear system in Theorem B.1. By Lemma 1, the Taylor remainders contribute at most Cλρ2 , freezing the sensitivity contributes at most Cλρ3 , and replacing the residual contributes at most Cλ2 ρ. Using ρ ≤ 1 and λ = o(αℓ ) gives (40). The next lemma gives the modes of this linear approximation and the growth bound used to integrate the remainder. Lemma 2 (Amplification of the linear singular modes). Every solution of the linear system in Theorem B.1 satisfies, for 0 ≤ s ≤ t ≤ T with λT ≤ M , lin lin lin ∥(wℓ,j (t), ω lin ℓ+1,j (t))∥ ≤ CM ∥(wℓ,j (s), ω ℓ+1,j (s))∥.

The projection identities (37) hold for these linear trajectories with both errors equal to zero. Proof. The linear system is generated by the symmetric matrix !

σ ′ (0) 0 A⊤ ℓ , √ pℓ Aℓ 0 whose eigenvalues lie in [−λ, λ]. Since λ(t − s) ≤ M , the scalar evolution along each eigendirection is bounded by a constant depending only on M , proving the norm bound. Projecting the row and column equations onto a singular pair and then taking their sum and difference gives ′

√

± ±σ (0)si t/ pℓ lin ⟨wℓ,j (t), vi ⟩ ± ⟨ω lin . ℓ+1,j (t), ui ⟩ = 2cji e

(42)

Adding and subtracting these identities gives (37). Both weights have the same growing amplitude, and their decaying components have opposite signs. For si = 0, both projections are constant. Proof of Theorem B.1. We combine the gradient approximation from Proposition 1 with the growth bound in Lemma 2 to control the error over [0, T ]. The initial norm is at most Kαℓ for a fixed K > 0, so Lemma 2 bounds the norm of the linear solution by KM αℓ throughout [0, T ]. Let E(t) be the norm of the difference between the trained pair and its linear approximation. Integrating the remainder bound and using the bounded linear evolution gives, while ρ ≤ r∗ , E(t) ≤ CM λ

Z t

[ρ(s)2 + αℓ ρ(s)] ds.

(43)

0

Suppose ρ(s) ≤ 2KM αℓ for 0 ≤ s ≤ t. Since λt ≤ M , this gives E(t) ≤ CM λt αℓ2 ≤ CM αℓ2 . For sufficiently small scales, ρ(t) ≤ KM αℓ + E(t) < 32 KM αℓ and 2KM αℓ < r∗ . The bound assumed on ρ holds initially and cannot reach equality before T , by this strict improvement. Thus the flow stays in the bounded neighborhood and continues to T , proving (36). To obtain the bounds on the shorter interval, integrate the linear column equation and use the bound on the linear row: ∥ω lin ℓ+1,j (t)∥ ≤ Cαℓ+1 + CM λt αℓ . If λT ≤ M αℓ+1 /αℓ , this is at most CM αℓ+1 , while the preceding error estimate gives E(t) ≤ CM αℓ αℓ+1 . Adding the error to the linear column bound proves the stated bound for the trained column. Finally, Lemma 2 gives the projections of the linear trajectory. Projection onto a unit singular vector cannot increase the error, so the same bound proves (37). 25

C

Proof of Theorem 2

We start by defining the continuous-time variant of Algorithm 2 studied in Theorem 2. At each iteration, the row correction, column correction and readout gradient update are each multiplied by the same step size dt. We then take dt → 0. A common positive learning-rate multiplier can be absorbed into time, whereas unequal multipliers change the relative dynamics. Algorithm 3 specifies these updates for a three-layer network. Algorithm 3 Three-layer BackLoFi: updates for the continuous-time limit Require: Forward Neural LoFi features h1⋆ , h2⋆ and readout w3 ; step size dt; duration T . 1: Set C = 0 and keep both Gaussian lifts fixed. 2: for each time step of size dt, up to time T do 3:

Select the coupled features. b n [y g 2 ⊗ z0 ]. Compute g 2 = w3 σ ′ (h2 ) and the leading singular pair (u, v) of E

Set ψ = ⟨u, g 2 ⟩, χ = ⟨v, z0 ⟩; hold these scores fixed for this step. 4:

Correct the first layer. h⋆1,a ← h⋆1,a + dt χ for every a.

5:

Refresh forward LoFi. Recompute z1 , K1 and set h2⋆ ← FwdLoFi2 (K1 , y).

6:

Add the column spike. Using the updated z1 , set b n [yψz1 ], q←E

7:

C ← C + dt u ⊗ q.

Train the readout. Recompute h2 , z2 , h3 , then set b n [(y − h3 )z2 ]. w3 ← w3 + dt E

8: end for 9: return h3 = ⟨w3 , z2 ⟩.

C.1

Full description of the Theorem and the underlying order parameters

We now give the full description of Theorem 2 and the order parameters that allow the equations to close. We study the continuous-time limit of Algorithm 3, followed by the width limits p1 → ∞ and p2 → ∞. Filtered features and accumulated coupling. The input map z0 is fixed, with kernel ⋆ (x) ∈ Rk1 and h⋆ (x) ∈ Rk2 describe the K0 (x, x′ ) = ⟨z0 (x), z0 (x′ )⟩. The filtered features h1,t 2,t ⋆ = FwdLoFi (K , y) and two hidden layers before their Gaussian lifts. The first starts from h1,0 1 0 ⋆ = FwdLoFi (K , y). follows the lower correction; the second is refreshed as h2,t 2 1,t Let r ∼ ν = N (0, Ik1 ) and η ∼ µ = N (0, Ik2 ) be the fixed lift rows. The scalar function Ct (η, r) records the total column coupling between these rows, starting from C0 = 0. Its action is (Ct f )(η) = Er [Ct (η, r)f (r)].

(44)

The hidden preactivations and activations are therefore ⋆ h1,t (r, x) = r⊤ h1,t (x),

z1,t = σ(h1,t ),

⋆ h2,t (η, x) = η ⊤ h2,t (x) + (Ct z1,t (·, x))(η),

z2,t = σ(h2,t ).

(45)

Thus earlier columns act on the current lower features through Ct z1,t . The arguments of Ct have dimensions k2 and k1 , independent of network width. Neuron inner products are Gaussian 26

expectations, or empirical averages at finite width; input coefficient inner products are Euclidean. √ This normalization absorbs the factors 1/ pj used for the feature vectors in the main text. Output and sensitivity. The scalar output is h3,t = ⟨w3,t , z2,t ⟩. To track how readout learning changes the sensitivity, retain the response St (η, s) = w3,t (η)σ ′ (s). Evaluating it at the current preactivation gives gt2 (η, x) = St (η, h2,t (η, x)). (46) The output h3,t alone does not determine this sensitivity when the hidden features change; St retains ⋆ , h⋆ , C , S , h the readout’s response to that change. Together, h1,t t t 3,t determine the feature and 2,t sensitivity kernels: K1,t (x, x′ ) = Er [z1,t (r, x)z1,t (r, x′ )], K2,t (x, x′ ) = Eη [z2,t (η, x)z2,t (η, x′ )],

Gt (x, x′ ) = Eη [gt2 (η, x)gt2 (η, x′ )].

(47)

Selection and continuous updates. At each time, (ψt , χt ) maximizes (15) for the current kernels Gt , K0 . The singular-vector equations become b n [Gt (x, X)yχt (X)], ρt ψt (x) = E

b n [K0 (x, X)yψt (X)], ρt χt (x) = E

(48)

where both features have unit RKHS norm and ρt > 0 is the achieved correlation. Algorithm 3 gives ḣ⋆1,t,a (x) = χt (x), a = 1, . . . , k1 , 1b 2 b ∂t Ct (η, r) = E n [yχt (X)gt (η, X)]En [yψt (X)z1,t (r, X)], ρt b n [(y − h3,t (X))z2,t (η, X)]. ∂t St (η, s) = σ ′ (s)E

(49) (50) (51)

The first equation corrects the lower filter. The second accumulates its paired column spikes, using 2 b b ut = ρ−1 t En [yχt gt ] and qt = En [yψt z1,t ]. The third is the final-layer gradient-flow update, with s held fixed. The output evolves by (19). Forward filter. The forward map FwdLoFiℓ (K, y) is Algorithm 1 expressed in feature space. Let HK be the RKHS of the current kernel K. Following Dandi et al. (2026b), its features are selected by b n [yf ], f0 ∈ arg max E ∥f ∥HK =1

fa ∈

arg max ∥f ∥HK =1 ⟨f,fb ⟩HK =0, 1≤b<a

b n [yf 2 ] , E

a = 1, . . . , kℓ − 1.

(52)

The first coordinate maximizes first-order correlation; the remaining coordinates successively maximize second-order correlation. Orthogonality is measured in HK among the second-order features. The map returns FwdLoFiℓ (K, y) = (f0 , f1 , . . . , fkℓ −1 ). At initialization, choose the sign of each second-order forward feature and of χ to have positive empirical inner product with a fixed reference function. At later steps, use its previous value as the reference, and change the sign of ψ together with that of χ.

27

Assumption 4 (Data, lifts and initialization). The training set is finite, and the numbers of filtered coordinates k1 , k2 are fixed. The lift laws are ν = N (0, Ik1 ) and µ = N (0, Ik2 ), or their finite empirical versions. All lift rows are held fixed in time. Initially C0 = 0. Choose h3,0 in the span of the initial features z2,0 on the training inputs, and set S0 by (55). Assumption 5 (Activation). The activation is C 2 , and σ, σ ′ , σ ′′ are bounded. Assumption 6 (Spectral gaps). On [0, T ], the label-weighted first moments used by the forward filters have norms bounded away from zero. The retained eigenvalues are uniformly separated in absolute value from zero, from one another, and from the discarded spectrum. The leading coupled singular value stays uniformly separated from zero and from the remaining singular values. Assumption 7 (Sequential width limit). The two lifts are sampled from independent nested Gaussian sequences. The data, ranks and initial prediction vector are the same across widths. The same reference functions initialize the signs at every width, with nonzero inner product with the limiting features. The initial upper Gaussian-population kernel matrix is positive definite. The population trajectory satisfies Assumption 6 on [0, T ]. Theorem C.1 (Precise version of Theorem 2). Under Assumptions 4–6, the current-state equations have a unique local solution. On [0, T ], the piecewise-linear interpolants of Algorithm 3 converge to it uniformly as dt → 0, with error OT (dt). Under Assumption 7, the continuous flows also converge almost surely, first as p1 → ∞ at fixed admissible p2 , then as p2 → ∞, to the solution of (17)–(19), with expectations taken over both Gaussian lifts. Convergence is uniform on [0, T ] for Ct , St in the fixed norms specified in (59), and on the data and [0, T ] for the filtered features, output, selected scores and current kernels K1,t , K2,t , Gt . Here an admissible fixed upper sample means that the mixed Gaussian-lower/empirical-upper flow realizes the prescribed initial predictions, has a stable initial kernel inverse, and satisfies Assumption 6 on [0, T ]. The width statement concerns sufficiently large upper widths, for which these conditions hold.

C.2

Proof of Theorem C.1

Proof sketch. To obtain a closed evolution and compare networks of different widths, we represent the accumulated columns by a kernel Ct (η, r) on the fixed Gaussian label spaces. Each selected feature pair contributes a rank-one increment, and Ct z1,t records the accumulated action on the current lower features. Differentiating Ct z1,t therefore retains both newly added columns and the action of earlier columns on moving features. Together with the readout update, this gives the output equation and the evolving sensitivity. The spectral gaps ensure that the selected features vary continuously with the state, yielding a locally Lipschitz flow. Each iteration of Algorithm 3 approximates (49)–(51) over an interval of length dt with error O(dt2 ), giving an OT (dt) error on [0, T ]. Finally, we compare empirical and Gaussian lifts on the same probability spaces. Stability of the flow lets us average the first lift and then the second, while retaining the dependence of learned features on the sampled neurons. We start by establishing the continuity of the maximizing feature pair with respect to the evolving order parameters. Lemma 3 (Continuity of the maximizing feature pair). Fix the training data, input feature map z0 and lift laws µ, ν. Under Assumptions 5 and 6, the selected singular pair (u, v) and scores (ψ, χ) depend locally Lipschitz-continuously on the current order parameters (h1⋆ , C, S), with 28

b n [∥f (X)∥2 ] h2⋆ = FwdLoFi2 (K1 , y) and signs chosen as above. Use the empirical norm ∥f ∥2n = E for the filtered features and scores, the Hilbert–Schmidt norm for C, and the L2 (µ; D) norm for S, where D = span{σ ′ } and ∥bσ ′ ∥D = |b|. The statement holds on bounded sets with ∥S(η, ·)∥D ≤ M almost everywhere, for fixed M < ∞.

Proof. The first-layer features z1 (r, x) = σ(r⊤ h1⋆ (x)) depend Lipschitz-continuously on h1⋆ : the derivative of σ is bounded and r has finite second moment. Their label-weighted first and second moments inherit this local regularity. Assumption 6 keeps the first moment nonzero and separates the retained eigenvalues of the second-moment operator. Normalizing the first moment and applying eigenvector perturbation bounds therefore give local Lipschitz continuity of the forward features h2⋆ = FwdLoFi2 (K1 , y). For two sets of order parameters (h1⋆ , C, S) and (h̄1⋆ , C̄, S̄), let h2 , h̄2 be the corresponding upper preactivations. Equation (45) and the bound ∥Cf ∥2 ≤ ∥C∥HS ∥f ∥2 give ∥h2 − h̄2 ∥2 ≤ L ∥h1⋆ − h̄1⋆ ∥n + ∥C − C̄∥HS .

The L2 norms on preactivations and sensitivities average over the upper label and training inputs. Since g 2 (η, x) = S(η, h2 (η, x)), the corresponding sensitivities satisfy ∥g 2 − ḡ 2 ∥2 ≤ ∥σ ′ ∥∞ ∥S − S̄∥L2 (µ;D) + M ∥σ ′′ ∥∞ ∥h2 − h̄2 ∥2 . b n [y g 2 ⊗ z0 ] is bounded These estimates give local Lipschitz continuity of g 2 in (h1⋆ , C, S). Since A = E and linear in g 2 , A∗ A also depends locally Lipschitz-continuously on these order parameters. Its leading eigenvalue ρ2 is positive and simple. Fixing the sign of its unit eigenvector v by positive overlap with a reference, the equations

A∗ Av = ρ2 v,

∥v∥2 = 1

have an invertible linearization. On v⊥ , the inverse of ρ2 I − A∗ A has norm at most 1/(ρ2 − ρ22 ), where ρ2 is the second singular value of A (zero when A has rank one). Assumption 6 keeps this eigenvalue difference bounded away from zero. The implicit function theorem gives smooth dependence of ρ and v on A, and hence local Lipschitz dependence on the order parameters. The identities u = Av/ρ, ψ = ⟨u, g 2 ⟩ and χ = ⟨v, z0 ⟩ give the same regularity for the left singular vector and both maximizing functions. We next verify that the column updates accumulate in Ct as defined in (44), and derive the evolution of their action Ct z1,t on the changing lower features. Proposition 2 (Accumulated column coupling). Let µ, ν be the Gaussian lift laws or their empirical versions. Suppose ut ∈ L2 (µ) and qt ∈ L2 (ν) are continuous in t. The column update Ċt = ut ⊗ qt , C0 = 0, has a unique continuously differentiable Hilbert–Schmidt solution, whose scalar kernel is given in L2 (µ ⊗ ν) by Z t

Ct (η, r) =

uτ (η)qτ (r) dτ.

(53)

0

It maps functions of the lower label r to functions of the upper label η through (44). For any continuously differentiable lower-feature curve z1,t , ∂t (Ct z1,t ) = Ċt z1,t + Ct ż1,t .

29

Proof. A rank-one update has kernel ut (η)qt (r) and Hilbert–Schmidt norm ∥ut ∥L2 (µ) ∥qt ∥L2 (ν) . Integrating in this Hilbert space gives the unique solution and the bound ∥Ct ∥HS ≤

Z t 0

∥uτ ∥L2 (µ) ∥qτ ∥L2 (ν) dτ.

The kernel therefore belongs to L2 (µ ⊗ ν). Since applying a Hilbert–Schmidt operator to an L2 (ν) function is continuous and bilinear, the product rule gives the displayed derivative. Combining this representation with the readout update gives the closed response and output equations. Proposition 3 (Closed evolution of the response and output). Under Assumptions 4–6, let ⋆ , C , w ) solve the lower-feature and column equations (49)–(50) and final-layer gradient flow (h1,t t 3,t under squared loss, with the prescribed initialization. Then St (η, s) = w3,t (η)σ ′ (s) satisfies (51), and h3,t = ⟨w3,t , z2,t ⟩ satisfies (19). These derivatives, and those of the feature maps and their kernels, ⋆ , h⋆ , C , S , h ). are determined by the current values of (h1,t t t 3,t 2,t Proof. The update rates are continuous in the order parameters by Lemma 3. Thus any solution of their integral equations is continuously differentiable in time. We derive the response and output identities for such a solution; its existence and uniqueness are established in the proof of Theorem C.1 below. b n [y g 2 ⊗ z0 ]. For a Write the readout as w3,t ∈ L2 (µ). The cross-layer operator is At = E t b n [yχt g 2 ]; applying the upper feature map gives nonzero singular pair, At vt = ρt ut gives ut = ρ−1 E t t b n [Gt (x, X)yχt (X)]. The adjoint equation similarly gives ρt χt (x) = E b n [K0 (x, X)yψt (X)]. ρt ψt (x) = E Each singular vector lies in the span of its feature map and therefore has the minimum coefficient norm among representations of the selected feature, which equals the RKHS norm of that feature. In particular, the unit vectors ut , vt give ∥ψt ∥HGt = ∥χt ∥HK0 = 1. Every lower filtered coordinate follows (49). The remaining continuous updates are Ċt = ut ⊗ qt ,

b n [(y − h3,t (X))z2,t (·, X)], ẇ3,t = E 2 b ut = ρ−1 t En [yχt (X)gt (·, X)].

b n [yψt (X)z1,t (·, X)], qt = E

(54)

⋆ (x) + (C z (x))(η). Substituting the two representers Here (u ⊗ q)b = u⟨q, b⟩, and h2,t (η, x) = η ⊤ h2,t t 1,t into Ċt gives (50). The loading qt fits the modified label on current lower activations, and ut selects the upper response. To initialize the response from h3,0 , choose coefficients b0 on the training inputs satisfying

b n [K2,0 (x, X)b0 (X)] = h3,0 (x) E

for every training input x.

b n [b0 (X)z2,0 (η, X)] realizes h3,0 and lies A solution exists by Assumption 4. The readout w3,0 (η) = E in the span of the initial activation profiles. Any other readout with the same predictions differs by a function orthogonal to that span. Hence w3,0 has minimum norm and is independent of the chosen coefficient solution. The corresponding initial response is b n [b0 (X)z2,0 (η, X)]. S0 (η, s) = σ ′ (s)E

(55)

Holding s fixed and applying the readout gradient update gives (51); evaluating at h2,t gives the current sensitivity in (46). Its total derivative therefore contains both readout learning and the 30

changing activation derivative: b n [(y − h3,t (X))z2,t (η, X)] ∂t gt2 (η, x) = σ ′ (h2,t (η, x))E

+ ∂s St (η, h2,t (η, x)) ḣ2,t (η, x).

(56)

The first term changes the response even with the hidden features held fixed. The second changes it as the upper preactivation moves, with ∂s St (η, s) = w3,t (η)σ ′′ (s). Both update Gt , hence the next selected scores (ψt , χt ), which in turn update the lower filter and the column coupling. Since ∂t St is σ ′ (s) times a function of η, the evolution preserves the product form St (η, s) = w3,t (η)σ ′ (s) of the initialization. Proposition 2 represents all installed columns by the current kernel Ct . The accumulated columns act on the current lower features through Ct z1,t . Differentiating this product gives ⋆ ḣ2,t (η, x) = η ⊤ ḣ2,t (x) + (Ċt z1,t (x))(η) + (Ct ż1,t (x))(η), ⋆ ż1,t (r, x) = σ ′ (r⊤ h1,t (x))χt (x)

k1 X

ra .

(57)

a=1

The three terms in the first line refresh the forward filter, install the current column, and propagate lower-feature motion through previously installed columns. The second line follows because every ⋆ = χ P r . The change in the first-layer kernel filtered coordinate has the same velocity: r⊤ ḣ1,t t a a then refreshes the second-layer forward filter: K̇1,t (x, x′ ) = Er ż1,t (r, x)z1,t (r, x′ ) + z1,t (r, x)ż1,t (r, x′ ) , 



⋆ ḣ2,t = DK FwdLoFi2 (K1,t , y)[K̇1,t ],

where DK FwdLoFi2 is the derivative of the forward filter with respect to its kernel. Thus all three derivatives are determined by the current filtered features and their prescribed velocity. Differentiate h3,t = ⟨w3,t , z2,t ⟩ to obtain b n [K2,t (x, X)(y − h3,t (X))] + Eη [g 2 (η, x)ḣ2,t (η, x)], ḣ3,t (x) = E t b n [K1,t (x, X)yψt (X)]. Eη [gt2 (η, x)(Ċt z1,t (x))(η)] = ψt (x)E

(58)

b n [K1,t (x, X)yψt (X)]. Substituting (57) The second identity uses ⟨gt2 , ut ⟩ = ψt and ⟨qt , z1,t (x)⟩ = E 2 gives the four terms of (19), including Eη [gt (Ct ż1,t )] for earlier columns. These derivatives also hold for Gaussian lifts. The readout equation keeps w3,t uniformly bounded over neurons on finite intervals; the formulas for ut , qt and Ct then give |ḣ1,t (r, x)| ≤ CT ∥r∥ and |ḣ2,t (η, x)| ≤ CT (1 + ∥η∥). Bounded activation derivatives and finite Gaussian second moments therefore justify differentiating the neuron expectations by dominated convergence. Conversely, the output equation and the response update give identical derivatives for h3,t and ⟨w3,t , z2,t ⟩. Their equal initialization therefore ⋆ = FwdLoFi (K , y), the variational pair and the preserves equality. Since gt2 = St (η, h2,t ) and h2,t 2 1,t representers in (54) are determined by the current state, as are the feature and kernel derivatives.

Proof of Theorem C.1. We now combine Lemma 3 and Propositions 2 and 3 to prove Theorem C.1. Proposition 2 gives the evolution of the accumulated columns and their action on the changing lower features; Proposition 3 gives the closed response and output equations. Lemma 3 supplies the local Lipschitz bounds used to establish existence and uniqueness, control the time-discretization error, and pass to the two successive width limits. The filtered features give the Gaussian projections r⊤ h1⋆ (x) and η ⊤ h2⋆ (x). The coupling C maps the lower activation z1 (·, x) ∈ L2 (ν) to the upper-preactivation correction Cz1 (·, x) ∈ L2 (µ), while 31

S(η, ·) gives the sensitivity evaluated at the corrected upper preactivation. Their convergence is therefore naturally formulated on the Hilbert spaces of functions of these Gaussian variables: we measure C in Hilbert–Schmidt norm and S in L2 (µ; D), with the response norm defined in Lemma 3. Together with the filtered features on the training inputs, this gives the state space Rnk1 × HS(L2 (ν), L2 (µ)) × L2 (µ; D)

(59)

for (h1⋆ , C, S). Here HS denotes Hilbert–Schmidt operators: the squared norm of C is the integral of its squared kernel over both neuron labels. Since S(η, s) = w3 (η)σ ′ (s), the chosen response norm satisfies ∥S∥2L2 (µ;D) = Eη [|w3 (η)|2 ]. The same identity holds for differences of two responses. We can therefore prove existence, uniqueness and convergence for the coefficient function w3 , then recover S bypmultiplying by the fixed function σ ′ . Since ∥χt ∥HK0 = 1, the reproducing property gives |χt (x)| ≤ K0 (x, x). Hence ⋆ ∥ḣ1,t ∥n =

p

k1 ∥χt ∥n ≤

q

b n [K0 (X, X)]. k1 E

The fixed input kernel and filter rank make this bound uniform in width. The canonical initial readout is a finite linear combination of bounded activations, and its update satisfies ∥ẇ3,t ∥∞ ≤ ∥σ∥∞ (∥y∥n + ∥σ∥∞ ∥w3,t ∥2 ). Lemma 3 gives local Lipschitz continuity of the right-hand sides of (49)–(51) in (h1⋆ , C, S), while the response coefficients are bounded. Integrating these right-hand sides from the prescribed initial state gives a contraction on bounded trajectories near the initial state over a sufficiently short interval: the displayed readout bound keeps the coefficients bounded, and the Lipschitz constant is multiplied by the interval length. This proves local existence and uniqueness. Grönwall bounds ∥w3,t ∥2 and hence ∥w3,t ∥∞ on finite intervals, so the solution extends while Assumption 6 holds. Proposition 3 then identifies its output with (19). We now show that Algorithm 3 converges to these equations as dt → 0, with the lifts fixed. Denote the changes over one iteration by ∆h1⋆ , ∆C, ∆S. The lower-feature update is exactly dt times (49). The column and readout updates use the features obtained after the earlier updates in that iteration. Those features differ from their values at the start of the iteration by O(dt). The Lipschitz bounds therefore give ⋆ ∆h1⋆ = dt ḣ1,t ,

∆C = dt ∂t Ct + OT (dt2 ), ∆S = dt ∂t St + OT (dt2 ), where the derivatives are the right-hand sides of (49)–(51), evaluated at the start of the iteration. The errors are measured in (59). Integrating the continuous equations over the same interval gives these first-order increments with OT (dt2 ) remainders. Let Ek be the error between the algorithm and the continuous solution at time k dt, measured in (59). While the iterates remain in a neighborhood of the continuous trajectory on which Assumption 6 holds uniformly, the error satisfies Ek+1 ≤ (1 + LT dt)Ek + CT dt2 ,

E0 = 0.

Summing this recurrence gives maxk dt≤T Ek = OT (dt). For sufficiently small dt, this bound keeps every iterate in that neighborhood up to time T . Linear interpolation gives the same uniform error 32

between steps. The readout bound above holds for the discrete updates as well. Continuity of the forward filter and Lemma 3 transfer this convergence to the features, output, kernels and selected scores. We next take the width limits. The learned features depend on the sampled lifts, so we compare the finite-width and Gaussian equations along their evolving trajectories. For the lower lift, write the sampled rows as r1 , . . . , rp1 . Partition the Gaussian label space into measurable cells V1 , . . . , Vp1 , each of probability 1/p1 , chosen to minimize p1 Z X

∥r − ri ∥2 dν(r).

i=1 Vi

Set T1,p1 (r) = ri for r ∈ Vi . This is the quadratic optimal transport map from the Gaussian law to the empirical distribution of the sampled rows. The Gaussian density and finite second moments ensure that such a minimizing partition exists (Villani, 2003, Theorem 2.12). Construct T2,p2 in the same way for the upper lift. To see why the transport cost tends to zero, divide a bounded box B into finitely many cells of diameter at most ε. In each cell, the fraction of sampled rows converges almost surely to its Gaussian probability by the strong law. Assigning as much Gaussian mass as possible to sampled rows within each cell leaves a vanishing unassigned mass inside B, with squared distance at most ε2 on the assigned mass. For the remaining assignments, use ∥r − ri ∥2 ≤ 2∥r∥2 + 2∥ri ∥2 . Inside B, squared norms are Rbounded and the unassigned mass vanishes. Outside B, the sampled second moment converges to B c ∥r∥2 dν by the strong law. Since the chosen partition minimizes the cost, lim sup ∥T1,p1 (r) − r∥2L2 (ν) ≤ ε2 + 4 p1 →∞

Z

∥r∥2 dν(r)

almost surely.

Bc

Letting the cell diameter tend to zero and the box grow proves convergence. The same argument applies to the upper lift, giving δ1,p1 := ∥T1,p1 (r) − r∥L2 (ν) → 0,

δ2,p2 := ∥T2,p2 (η) − η∥L2 (µ) → 0

almost surely. Extend each sampled neuron profile to be constant on its cell. Because every cell has probability 1/pj , Gaussian integrals of these extensions equal the original empirical averages. Extend Ct on products of lower and upper cells, and St on upper cells. These extensions preserve the norms and all update equations. The finite-width and limiting solutions can therefore be compared in the same Hilbert space (59). At fixed order parameters, replacing r by T1,p1 (r) changes z1 by at most CT δ1,p1 in L2 . The first-layer kernel and refreshed forward filter satisfy the same bound. Replacing the upper lift changes the upper preactivation by at most CT (δ1,p1 + δ2,p2 ). Lemma 3 gives this bound for the resulting change in the selected pair and hence in the update equations. The constants are uniform in width when the response coefficients and lift moments are bounded and Assumption 6 holds uniformly. The initial responses also converge: their coefficients solve the kernel linear system above (55), whose limiting matrix is invertible by assumption. Let E(t) denote the distance between the two solutions after this identification, in the norm of (59). Integrating the difference of their update equations and applying Grönwall gives

sup E(t) ≤ CT E(0) + δ1,p1 + δ2,p2 . t≤T

33

(60)

Choose a neighborhood of the limiting trajectory on which Assumption 6 holds uniformly. The estimate applies until the finite-width solution leaves this neighborhood. Its right-hand side tends to zero, so for sufficiently large widths the solution stays in the neighborhood throughout [0, T ]. To obtain the sequential limit, fix p2 and an admissible set of upper lift rows. The two flows compared in (60) use these same rows, so the upper-lift error is zero. As p1 → ∞, the remaining error tends to zero almost surely. The limit therefore uses Gaussian expectations in the first layer and empirical averages over the fixed rows in the second. We then let p2 → ∞ and apply the same estimate to this intermediate evolution, obtaining the equations with Gaussian expectations in both layers. The positive-definite initial limiting kernel and stability under the spectral-gap assumption ensure that all sufficiently large upper samples are admissible almost surely. Continuity of the forward filter and the maximizing pair gives convergence of the filtered features and selected scores. Cauchy–Schwarz then gives convergence of the output and kernels from their expressions as inner products. Together with (60), these bounds establish almost-sure convergence uniformly on [0, T ] in the stated norms. This completes the proof of Theorem C.1.

D

The exact linear reduction

We first derive (13) for the one-hidden-layer linear network f (x) = W2 W1 x. Take unnormalized linear activations, unit-rate squared-loss gradient flow, whitened inputs E[xx⊤ ] = I, and possibly vector-valued targets. Write S = E[yx⊤ ]. Direct differentiation gives Ẇ1 = W2⊤ A∆ ,

A∆ = S − W2 W1 ,

Ẇ2 = A∆ W1⊤ .

(61)

Thus the incoming row and its matching outgoing column obey the two actions in (6) exactly, with the residual input–output correlation A∆ as coupling matrix. The prediction term produces saturation as the network fits the target. P Let S = i si ui vi⊤ and choose orthonormal hidden directions qi . The aligned initialization W1 =

X

αi qi vi⊤ ,

W2 =

X

i

βi ui qi⊤

(62)

i

is preserved: substituting it into (61) and using orthogonality gives α̇i = βi (si − αi βi ),

d 2 (α − βi2 ) = 0. dt i

β̇i = αi (si − αi βi ),

(63)

For a balanced positive seed, mi = αi βi satisfies ṁi = 2mi (si − mi ), hence mi (t) =

si . 1 + (si /mi (0) − 1)e−2si t

When mi (0) ≪ si , each factor initially grows as esi t , and their product approaches si . Alignment is needed for independent scalar modes; an arbitrary initialization also has cross-mode interactions. Without alignment, linearizing (61) at zero still gives, for each row–column pair, d w1,j dt ω 2,j

!

=

0 S⊤ S 0

!

!

w1,j . ω 2,j

√ The discarded terms are cubic. The block matrix has eigenvectors (vi , ±ui )/ 2 with eigenvalues ±si , which proves the shared exponential amplification near zero for arbitrary initial directions. This connects the linear dynamics to Theorem 1; the nonlinear surrogate additionally freezes gates and replaces the residual by a label-based approximation. 34

More than two weight matrices. For f (x) = WD · · · W1 x, aligned singular directions reQ 1 2 duce each mode to layer amplitudes a1i , . . . , aDi . Writing mi = D k=1 aki , its loss is 2 (si − mi ) . Differentiating gives the mode equations of Saxe et al. (2014, Section 2): ȧki = (si − mi )

Y

ari ,

ṁi = (si − mi )

r̸=k

D X

2

 Y 

k=1

ari  .

(64)

r̸=k

Equal positive initialization across layers is preserved by these equations. Then every aki = 1/D mi ,which gives (14). For D = 2, this is the balanced two-factor equation ṁi = 2mi (si − mi ). For 2−2/D D > 2, the early growth law is ṁi ≃ Dsi mi . The additional powers arise from multiplying the contributions of all layers. The constant-rate exponential approximation for two weight matrices therefore does not describe the simultaneous growth of all layers in a deeper network.

E

The Deep Staircase: A Theoretical Setting

We analyze here the deep-staircase model of Section 4. We first define the complete teacher, then specialize the forward and backward spectral operations of Algorithms 1 and 2 to this setting. We derive the population mechanism separating the forward and label-corrected statistics, and the resulting sample-complexity scalings. For simplicity, we do not provide full finite-sample proofs. These derivations can be turned into rigorous recovery guarantees by combining the population calculations below with standard concentration and spectral-perturbation arguments, as in Tabanelli et al. (2026); Wortsman-Zurich et al. (2026)

E.1

Teacher and high-dimensional regime

Target. Let x ∼ N (0, Id ) and let Φ2 (x) ∈ RD2 (d) , D2 (d) = d(d + 1)/2, denote the complete normalized quadratic Hermite representation of x. We draw 



AU   A⋆ =  AV  ∈ Rd1 ×D2 (d) , AW

d1 = pU + pV + pW ,

iid





(A⋆ )ij ∼ N 0, D2 (d)−1 ,

and define the first-level features hU = AU Φ2 (x),

hV = AV Φ2 (x),

hW = AW Φ2 (x).

(65)

This prior makes the planted directions dense and delocalized in the quadratic feature space. In the high-dimensional regime below, A⋆ A⊤ ⋆ → Id1 in operator norm. At the second level, let BQ , BZ ∈ RpU ×pU be symmetric and define ⊤ B h − tr(B )
hU h⊤ BZ hU − tr(BZ ) Q U Q √ √ , Z(hU ) = U , S = sign hV⊤ ChW . (66) 2 2 We assume that BQ and BZ are symmetric, have unit Frobenius norm and are Frobenius orthogonal, and that there exists an orthogonal transformation T such that T ⊤ BQ T = BQ and T ⊤ BZ T = −BZ . Thus, in the Gaussian latent reference, Q and Z are centered, unit-variance and uncorrelated, with d (Q, Z) = (Q, −Z). We take pV = pW = ph , let C ∈ Rph ×ph be an independent hidden coupling matrix. The raw target is   √ yraw = tanh (Q + ZS)/ 2 + b0 . (67)

Q(hU ) =

We finally standardize yraw of mean µy and variance σy2 to y = (yraw − µy )/σy with zero mean and unit variance. 35

High-dimensional regime.

For the high-dimensional regime we take

pU , pV , pW ≍ dϵ ,

d1 ≍ dϵ ,

0 < ϵ < 1,

n ≍ dα .

(68)

For the hidden-subspace sample complexity result we set pV = pW = ph and assume that C is full rank with balanced singular values −1/2

c− ph

−1/2

≤ σj (C) ≤ c+ ph

0 < c− ≤ c+ < ∞.

,

(69)

This assumption fixes the scale ∥C∥F = Θ(1) while excluding hidden directions whose interaction strength vanishes asymptotically relative to the others.

E.2

BackLoFi in the deep-staircase setting

Student architecture and Hermite specialization. We first relate the estimator used below to the general Neural LoFi architecture. In the numerical model, the raw input is first mapped through a frozen nonlinear random feature map z0 (x) = σ(W0 x). Each subsequent Neural LoFi layer has the reduce–lift form hℓ⋆ = Fℓ⊤ zℓ−1 ,

hℓ = Rℓ hℓ⋆ ,

zℓ = σ(hℓ ),

where Fℓ is estimated from the data and Rℓ is a fixed random lift. Thus h1⋆ plays the role of the first teacher level (hU , hV , hW ), while the next filtered coordinates h2⋆ represent second-level features such as Q and Z. For the theoretical analysis, we replace all the random feature map by its complete degree-two Hermite counterpart, for instance, for the first level we consider z0 (x) = Φ2 (x). The numerical estimator keeps the random feature maps. Forward pass.

For a scalar sample weight ξ and a representation z, define h

i

h

i

b n (ξ − ξ)zz ¯ ⊤ . Cbz [ξ] = E

b n (ξ − ξ)z ¯ , b z [ξ] = E m

(70)

We denote by mz [ξ] and Cz [ξ] the corresponding population quantities obtained by replacing the empirical average and empirical mean by their expectations. The first Forward step acts on z0 = Φ2 (x) and, in the present setting, uses Cbz0 [y] to recover the visible hU subspace. Once this representation has been learned, the upper Forward step acts on z1 = Φ2 (hU ) the quadratic Hermite representation of the recovered visible coordinates and uses mz1 [y] and Cz1 [y] respectively its first and second Hermite channels to estimate the two relevant second-level directions Q and Z respectively. Their precise population components are computed in the next subsection. Finally, therecovered  second-level coordinates are passed through the last (2) frozen random-feature map, z2 = σ R2 h , and a scalar ridge readout is fitted, bn b out = arg min E w w



⊤

y − w z2

2 

36

+ β∥w∥2 ,

⊤ b out fb = w z2 .

Backward score. is

In the general BackLoFi notation, the cross-layer operator at the first boundary h

i

b n y g2 (x)z0 (x)⊤ = UA ΣV ⊤ A=E A

where g2 is the sensitivity obtained by differentiating the fitted upper predictor with respect to its preactivation coordinates. In the Hermite specialization, after fitting the upper Hermite predictor fbH . The cross-layer SVD defined above retains left singular vectors that define the sample-dependent score si (x) = u⊤ (71) i g(x). In the present setting we retain only the row, or effective-label, channel of BackLoFi. The right singular vectors are not directly installed and no paired column correction is performed. The information returned to the lower spectral estimator is therefore entirely contained in the scores si . Second Forward pass. For each retained score, consider the effective label y si , we apply the same lower covariance estimator as in the first Forward Cbr0 [y si ] pass. After projecting out the already recovered visible subspace, the leading directions of this operator provide candidates for the missing hV , hW subspace. If several score channels are retained, their projected covariance operators are combined as specified in Algorithm 4. The selected directions are then appended to b (1) . At the second level, we the visible hU subspace to form the corrected first-level representation h b (1) set zB [y] and CzB [y]. The resulting 1 = Φ2 (h ) and apply again the Forward LoFi moments mzB 1 1 second-level representation is finally used to refit the predictor against the original labels.

E.3

The staircase mechanism in population

√ + b0 and define Target decomposition and parity. For the population proof only, set t± = Q±Z 2

a(Q, Z) =

1 tanh(t+ ) + tanh(t− ) − 2µy , 2σy

b(Q, Z) =

1 tanh(t+ ) − tanh(t− ) . 2σy

(72)

Then the standardized target satisfies exactly y = a(Q, Z) + b(Q, Z)S,

(73)

where a is even in Z while b is odd in Z. Their expansions around (Q, Z) = (0, 0) make the first relevant components explicit:    tanh′ (b0 ) tanh′′ (b0 )  2 √ Q+ Q + Z 2 + O ∥(Q, Z)∥3 , 4σy 2 σy ′ ′′   tanh (b0 ) tanh (b0 ) b(Q, Z) = √ Z+ QZ + O ∥(Q, Z)∥3 , 2σy 2 σy

a(Q, Z) = a0 +

(74)

with a0 = (tanh(b0 ) − µy )/σy . The corresponding expansion in hU follows immediately, since both Q and Z are quadratic features of hU . Lower-level Forward operator. We first analyze the population counterpart of the lower Forward estimator, since h(1) = (hU , hV , hW ) = A⋆ Φ2 (x), its restriction to the planted first-level coordinates is h

i

⊤ A⋆ Cz0 [y]A⊤ A⊤ ⋆ = A⋆ E y Φ2 (x)Φ2 (x) ⋆ =E

37

h

i

a(Q, Z) + b(Q, Z)S h(1) h(1)⊤ .

Algorithm 4 Covariance-row BackLoFi for the deep staircase Require: Data {(xµ , yµ )}nµ=1 , first-level dimensions (pU , pV , pW ), retained upper rank k2 . bz [y]. Retain its leading 1: Form the prescribed lower representation z0 = σ(W0 x) and compute C 0

pU directions by absolute eigenvalue and denote their projector by PbU . b U be the corresponding visible coordinates and form 2: Let h b U ). z1 = Φ2 (h b z1 [y] to estimate Q, then project this direction Apply upper Forward LoFi: use the direction of m out of the second-order statistic and retain its leading eigenvector to estimate Z. 3: Fit the upper predictor fbH and compute its fitted sensitivity g. 4: Compute the cross-layer SVD h

i

b n y g z⊤ = UA ΣV ⊤ . Ab = E 0 A 5: for I ∈ {{1}, {2}, {1, 2}} do 6: 7: 8: 9:

for i ∈ I do b Set si = u⊤ i g and form Cz0 [ysi ]. end for Set PbU⊥ = I − PbU and form TI =

X

PbU⊥ Cbz0 [ysi ]PbU⊥

2

.

i∈I

Retain the leading pV + pW eigenvectors of TI and append them to the protected visible b (1) . subspace, giving the corrected first-level representation h B 11: Form the corrected upper quadratic representation 10:

(1)

b zB 1 = Φ2 (hB ). 12:

Rerun Forward LoFi with the original labels y: retain the normalized direction selected b zB [y], then retain k2 − 1 additional directions from the corresponding deflated secondby m 1

(2)

b . Hermite/Wick operator. Denote the resulting upper coordinates by h B 13: Refit the scalar predictor against the original labels y by ridge regression on h

(2)

(2)

i

b ; Φ2 (h b ) . h B B

14: end for 15: Select I by validation error, and return the corresponding predictor and corrected first-level

representation.

Hence the U, U , V, V , W, W and V, W blocks of this matrix are exactly the corresponding population blocks of the operator diagonalized by Forward LoFi. In the Gaussian latent limit, the three first-level blocks are independent standard Gaussian vectors. In the decomposition above, a(Q, Z) is even in Z, while b(Q, Z) is odd. Since the output is centered and the hidden sign S is balanced, E[a(Q, Z)] = 0, while the symmetry Z 7→ −Z gives E[b(Q, Z)] = 0. Hence the complete hidden

38

block of the first Forward operator vanishes: "

⊤ hV hV⊤ hV hW E y ⊤ hW hV⊤ hW hW

!#

"

⊤ hV hV⊤ hV hW = E[a(Q, Z)] IpV +pW + E[b(Q, Z)] E S ⊤ hW hV⊤ hW hW {z } | {z } | =0

!#

= 0.

(75)

=0 by parity

The mixed visible-hidden blocks vanish by the same sign symmetries. By contrast, the visible block is nonzero: h i ⊤ AU Cz0 [y]A⊤ = E a(Q, Z)h h (76) U U ̸= 0, U accordingly to the expansion of a(Q, Z) in 74. Thus the first Forward operator contains signal on hU while its entire hV ⊕ hW block vanishes. Upper-level Forward operator. Once the visible subspace has been recovered, the next Forward step acts on z1 = Φ2 (hU ). Averaging over the hidden blocks removes the b(Q, Z)S contribution and gives EhV ,hW [y | hU ] = a(Q, Z). Since Q and Z are quadratic coordinates of hU , they are linear directions of z1 . In population, the expansion (74) enters directly into the two Forward statistics, at the leading order it yields: mz1 [y] = E[a(Q, Z)z1 ] =

tanh′ (b0 ) √ E[Qz1 ], 2 σy | {z }

(77)

∝ BQ

h

i

Cz1 [y] = E a(Q, Z)z1 z⊤ 1 =

′′

i tanh (b0 ) h 2 E (Q + Z 2 ))z1 z⊤ 1 . 4σy | {z }

(78)

spikes along BQ and BZ

The first moment identifies the weight BQ associated with Q. We then restrict the pure secondHermite operator to the orthogonal complement of BQ , so that its leading informative direction is BZ . Projecting these two recovered weights back onto z1 gives the scalar coordinates Q and Z. We therefore retain the two-dimensional upper representation h(2) = (Q, Z) ∈ R2 . These are intermediate features of the same scalar hierarchy, which are subsequently combined by the upper readout and provide the two channels used by the backward pass. Upper readout. After recovering h(2) = (Q, Z), the Forward student applies the final randomfeature map z2 = σ(R2 h(2) ) and fits a linear ridge readout. To see the population components selected by this fit, let r⊤ i = (ri,Q , ri,Z ) denote the i-th row of R2 . The corresponding random-feature coordinate is z2,i = σ(ri,Q Q + ri,Z Z) , whose low-order Hermite expansion is h

2 z2,i = α0 + α1 (ri,Q Q + ri,Z Z) + α2 ri,Q H2 (Q) +

√

Since the visible regression function a(Q, Z) is even in Z, 2 E[a(Q, Z)z2,i ] = α1 ri,Q E[aQ] + α2 ri,Q E[aH2 (Q)] 2 E[aH2 (Z)] + · · · , + α2 ri,Z

39

i

2 2 ri,Q ri,Z QZ + ri,Z H2 (Z) + · · · .

whereas the linear Z and QZ components vanish by parity. In the Gaussian RF limit, the population feature covariance is block-diagonal across Hermite sectors, so the ridge fit only rescales these components. Consequently, the first relevant terms of the fitted predictor have the form fb(Q, Z) = γQ Q + γQQ H2 (Q) + γZ H2 (Z) + · · · . Backward score. Consequently, the fitted upper sensitivity with respect to the second-level coordinates satisfies ! √ γQ + 2 γ Q + · · · QQ b √ g(Q, Z) = ∇(Q,Z) f (Q, Z) = . (79) 2 γZ Z + · · · The relevant distinction is their parity in Z: the first response is even in Z, whereas the second is odd and starts linearly in Z. Since the hidden coefficient b(Q, Z) is itself odd in Z, only the latter can produce a non-zero hidden coupling after backward modulation. Indeed, let AZ denote the cross-layer row associated with the second component of (79). We identify quadratic Hermite coefficient vectors with their symmetric-matrix representation and keep the same notation for both. We approximate the latent moments by replacing (hU , hV , hW ) with independent standard Gaussian vectors (u, v, w) of the same dimensions. The actual latent blocks are quadratic functions of the same input x; their joint law is not exactly Gaussian. We sketch the reason for this approximation below. In Eu , Q, Z are evaluated at hU = u; in Ev,w , S = sign(v⊤ Cw), with C fixed. The second-order V, W contraction gives #

"

xx⊤ − Id AZ = E a(Q, Z) + b(Q, Z)S gZ (Q, Z) √ 2 p p W V XX √ = 2 Eu [b(Q, Z)gZ (Q, Z)] Ev,w [S va wb ] (AV,a AW,b + AW,b AV,a ) + RZ ,

(80)

a=1 b=1

where RZ collects the other contractions and the error from replacing the dense quadratic latents by their Gaussian reference. The product of two normalized quadratic features satisfies the exact identity D

hV,a hW,b = tr(AV,a AW,b ) + AV,a AW,b + AW,b AV,a , xx⊤ − Id |

{z

second Hermite order

E F

+ ProjH4 (hV,a hW,b ) .

}

fourth Hermite order

|

{z

(81)

}

The constant term has zero projection onto the quadratic Hermite features. For the leading response √ gZ = 2γZ Z with γZ ̸= 0, the displayed √ contraction is nonzero: strict monotonicity of tanh gives b(Q, Z)Z > 0 for Z ̸= 0, so Eu [bgZ ] = 2γZ Eu [bZ] ̸= 0. Its size is calculated in Section E.4. If A = UA ΣVA⊤ and a retained left singular vector ui = (ui,Q , ui,Z ) has ui,Z ̸= 0, its score has the form si (Q, Z) = u⊤ ηi ̸= 0 (82) i g(Q, Z) = si,even (Q, Z) + ηi Z + · · · , Sensitivity through the random lift. The sensitivity used below is written with respect to the filtered second-level coordinates h(2) = (Q, Z), whereas the general BackLoFi construction uses the sensitivity with respect to the neuron preactivations. Here these preactivations are R2 h(2) . If gneu = ∇R2 h(2) fb,

40

then the chain rule gives g := ∇h(2) fb = R2⊤ gneu ,

h

i

h

i

⊤ ⊤ A = E y g z⊤ 0 = R2 E y gneu z0 .

(83)

Thus the two-dimensional cross-layer operator considered below is the neuron-level operator restricted to the subspace generated by the random lift. At finite width, orthonormalizing this subspace amounts to replacing g by (R2⊤ R2 )−1/2 g. For an iid Gaussian lift of fixed rank two, 1 ⊤ R R2 −→ I2 , p 2 −1/2 asymptotically, which does not so this only multiplies the score u⊤ i g by the common factor p change the recovered lower eigenspace. We therefore use the simpler reduced-coordinate sensitivity g below.

Second Forward pass. For each retained score si in (82), the backward correction simply reruns the lower Forward covariance with the effective label ysi . In planted coordinates, h

i

(1) (1)⊤ A⋆ Cz0 [ysi ]A⊤ . ⋆ = E ysi h h

(84)

Its V, W block is, at leading order in the odd component of the score, h

i

h

i

⊤ ⊤ E ysi hV hW = E[b(Q, Z)si (Q, Z)] E S hV hW ̸= 0

|

{z

ηi E[b(Q,Z)Z]+···

(85)

}

The even part of si does not contribute because b(Q, Z) is odd in Z. Thus every retained mode with ηi ̸= 0 produces a covariance spike on the hidden V, W subspace, which was absent from the ordinary Forward operator. When several scores are retained, their projected covariance operators are combined using the same squared-covariance aggregation as in the estimator, before extracting the hidden eigenspace. Since C has full rank, the hidden restriction of Cz0 [ysi ] has range hV ⊕ hW . Its leading pV + pW eigendirections therefore recover the joint hidden first-level subspace. We do not need to identify the V and W bases separately: the subspace is recovered up to an orthogonal change of coordinates. Projecting the data onto the visible and recovered hidden directions yields the corrected first-level representation
b (1) = h bU , h bV W . h B Corrected upper Forward pass. We then resume the ordinary Forward hierarchy using the original labels y and form   b (1) . zB 1 = Φ2 hB At population level, the relevant quadratic coordinates in this dictionary are Q, Z, and T = hV⊤ ChW . The coefficients of mzB [y] along the relevant quadratic directions satisfy 1

E[yQ] ̸= 0,

E[yZ] = 0,

E[yT ] = 0.

We retain this direction the Q direction and then apply the covariance estimator CzB [y] operator 1 on its orthogonal complement. The new feature made available by the backward correction is the nonzero coupling 



E[yZT ] = E a(Q, Z) + b(Q, Z)S ZT = E[b(Q, Z)Z] E|T | > 0, 41

(86)

where S = sign(T ). Hence the Q-deflated upper operator has a nontrivial component on span{Z, T }. Under a spectral-gap assumption, the corrected upper Forward step retains a rank-k2 representation (2) hB ∈ Rk2 whose informative population span contains Q together with span{Z, T }. The scalar  (2) (2)  predictor is then refitted against the original labels y using a ridge readout on hB ; Φ2 (hB ) Since the retained upper span contains both Z and T , this quadratic readout contains the interaction ZT , providing a non-zero low-degree prediction channel for the hidden branch.

E.4

Finite-sample scaling and the learned-score bottleneck

We now translate the population gaps of the previous subsection into finite-sample scales. The purpose is not to prove a sharp concentration theorem for the complete learned-score pipeline, but to identify the effective dimension and spectral gap controlling each stage. Throughout, D0 = D2 (d) ≍ d2 ,

pU , pV , pW ≍ dϵ ,

0 < ϵ < 1.

We suppress logarithmic factors from truncating the Hermite features and constants independent of d and n. Matrix concentration and subspace recovery. We consider the weighted covariances defined in (70): h i M = Cz [ξ] = E (ξ − Eξ)zz⊤ , D = dim(z). The lower Forward pass uses Cz0 [y], and the backward correction uses Cz0 [ysi ]. Both act on z0 = Φ2 (x), so D = D2 (d) = d(d + 1)/2. The upper Forward covariance is Cz1 [y], with z1 = Φ2 (hU ) c=C bz [ξ] be the empirical matrix from (70), and D = D2 (pU ) = pU (pU + 1)/2. In each case, let M b c. and let P and P project onto a target eigenspace of M and the corresponding eigenspace of M For independent samples, when the truncated centered summands have operator norm and pere sample matrix variance O(D) and the truncation bias is negligible, the matrix Bernstein inequality (Vershynin, 2018, Theorem 5.4.1) gives s c − M ∥op = O eP  ∥M



D D + . n n

(87)

For the truncation and concentration argument for Hermite features, see Wen et al. (2025, Appendix F.2.1, Lemmas F.4–F.5), in particular the proof of Lemma F.5. In the regime n ≫ D considered below, the square-root term dominates. If the target eigenspace is separated by a population gap ∆, the Davis–Kahan theorem (Yu et al., 2015) then gives the sufficient condition, up to logarithmic factors, D n≫ 2 (88) ∆ for consistent recovery of its projector. We use the analogous statement for the cross-layer matrix. If A is a rectangular population b corresponding left singular vectors separated by a operator, Ab its empirical estimate, and u and u singular-value gap ∆A , Wedin (1972) gives the same sufficient scaling, n≫

D , ∆2A

(89)

where D is the effective dimension controlling ∥Ab − A∥op . These estimates apply directly to fixed representations and scores. Reusing the samples that fitted them requires a separate dependence argument. 42

Visible Forward recovery. The lower Forward covariance acts on Φ2 (x) ∈ RD0 . From the −1/2 population calculation above, the smallest visible spike satisfies ∆U ≍ pU . Applying (88) therefore gives nU ≫ D0 pU ≍ d2+ϵ . (90) This is the scale at which the Forward estimator recovers the visible U -subspace. Its hidden population block remains zero, so this argument does not reveal the V, W directions. Upper Forward recovery.

Once U has been recovered, the upper problem lives in DU = D2 (pU ) ≍ p2U

dimensions. Assuming that the first-moment coefficient selecting Q and the Q-deflated Wick gap selecting Z remain bounded away from zero, the same effective-dimension argument gives nQ,Z ≫ p2U ≍ d2ϵ .

(91)

Since 2ϵ < 2 + ϵ, this stage is not a bottleneck. Hidden recovery given an informative score. Let ph = pV = pW . From the population V, W block derived above, the nonzero hidden eigenvalues have the same order as the singular values of C. Under the balanced-singular-value condition (69), the smallest hidden population gap scales as −1/2 ∆H ≍ ph . The operator again acts in the D0 -dimensional quadratic input space, so (88) yields nV,W ≫ D0 ph ≍ d2+ϵ .

(92)

Learned-score bottleneck. The practical method must learn the informative score from the cross-layer matrix h i b y g(x)Φ2 (x)⊤ . Ab = E We now compute the size of the V, W contraction in (81). For this calculation, we replace the actual latent blocks (hU , hV , hW ) by independent standard Gaussian vectors (u, v, w) of the same dimensions. We first explain why this is a natural approximation and what it establishes. Fix a teacher row and write A for its symmetric-matrix representation. Diagonalizing A gives d x⊤ Ax − tr A 1 X √ √ = λr (A)(ξr2 − 1), hA (x) = 2 2 r=1

iid

ξr ∼ N (0, 1).

For the dense Gaussian teacher draws, ∥A∥2F → 1 and ∥A∥op = OP (d−1/2 ). Thus the total variance tends to one, while the largest summand has variance OP (d−1 ). The Lindeberg central limit theorem then gives a standard Gaussian limit. The same argument applies to every fixed linear combination of teacher rows. Moreover, Ex [hAi hAj ] = tr(Ai Aj ) → δij . Consequently, conditional on typical teacher draws, any fixed collection of latent coordinates converges jointly, over the randomness of x, to independent standard Gaussians. This is also a degree-two instance of the Gaussian approximation criterion in Wen et al. (2025, Theorem 4.1). This argument explains the independent Gaussian variables used in the coefficient calculation. It does not yet bound the replacement error when all pU + pV + pW ≍ dϵ coordinates enter Q, Z, S. In particular, convergence in distribution alone does not show that the resulting error in the cross-layer row is smaller than d−1/2 . That error remains part of RZ in (80). 43

We therefore compute the following coefficients in the Gaussian approximation: evaluate Q, Z at u, set S = sign(v⊤ Cw), and keep C fixed. Then Hab = Ev,w [S va wb ],

βZ = Eu [b(Q, Z)gZ (Q, Z)],

Kd =

X

Hab (AV,a AW,b + AW,b AV,a ).

a,b

Under the norm-preserving identification of quadratic Hermite coefficient √ vectors with symmetric matrices, the displayed contribution to the Z row in (80) is AVZ W = 2 βZ Kd . Conditional on C, H is fixed when averaging over the independent V, W teacher rows. Conditioning on v gives r

Ew [Sw | v] =

2 C ⊤v , π ∥C ⊤ v∥

r

⟨H, C⟩F =

2 Ev ∥C ⊤ v∥. π

The singular-value bounds in (69) and Cauchy–Schwarz therefore imply r

2 c− Ev ∥v∥ ≤ ∥H∥F ≤ 1. √ π c+ ph

For the upper bound, the products va wb are orthonormal in Gaussian L2 , so Bessel’s inequality gives P √ 2 2 a,b (E[Sva wb ]) ≤ E[S ] = 1. Since E∥v∥/ ph → 1, these bounds show ∥H∥F = Θ(1) directly for the prescribed matrix C. Each teacher row, viewed as a symmetric matrix A, has Gaussian coefficients of variance D0−1 in an orthonormal Frobenius basis. Writing EA for expectation over these draws, EA [Aij Akl ] =

δik δjl + δil δjk , 2D0

EA [A2 ] =

Id . d

For two independent rows A, B, Gaussian fourth moments give EA ∥AB + BA∥2F = 2EA tr(A2 B 2 ) + 2EA tr(ABAB) 2 1 d 2 = + + 2 =: νd = + O(d−2 ). d D0 D0 d

(93)

In EA ∥Kd ∥2F , all cross terms between distinct row pairs vanish by independence and centering. Hence 96νd EA ∥Kd ∥2F = νd ∥H∥2F , VarA (∥Kd ∥2F ) ≤ ∥H∥4F . (94) D0 For the variance bound, rotate the independent Gaussian row families using an SVD of H, so that P Kd = j hj (Aj Bj + Bj Aj ), where hj are the singular values of H. With F = ∥Kd ∥2F , its gradient is ∇Aj F = 2hj (Kd Bj + Bj Kd ), and similarly for Bj . The Gaussian Poincaré inequality, with coefficient variance D0−1 , gives 



X 16 96 VarA F ≤ EA F h2j (∥Aj ∥2F + ∥Bj ∥2F ) ≤ EA [F ]∥H∥2F . D0 D 0 j

The last inequality follows by conditioning on either row family and using the Gaussian fourthmoment bound E[XY ] ≤ 3E[X]E[Y ] for nonnegative quadratic forms X, Y of the same centered Gaussian vector. Since the variance in (94) is O(d−1 ) times the squared mean, Chebyshev’s inequality yields √ ∥H∥F |βZ |∥H∥F √ ∥Kd ∥F = ( 2 + oP (1)) √ , ∥AVZ W ∥2 = (2 + oP (1)) . (95) d d 44

Stage Visible U recovery Upper Q, Z recovery Learned backward score Hidden V, W recovery given the score Corrected upper Forward pass

Predicted sample scale d2+ϵ d2ϵ d3 d2+ϵ d2ϵ

Table 1: Stagewise effective-dimension scales for the quadratic staircase. The d3 score-selection scale uses the contraction size derived in (95) and requires this contribution to remain spectrally isolated in the complete cross-layer matrix. Thus the displayed contraction has size ΘP (d√−1/2 ) when |βZ | = Θ(1). For the retained linear √ response gZ = 2γZ Z, this coefficient is βZ = 2γZ Eu [bZ]. In the numerical specialization, (Q, Z) √ ′ converges to two independent standard Gaussians and Eu [bZ] → E[tanh (G + b0 )]/( 2σy ) > 0, where G ∼ N (0, 1) and σy denotes its limiting value. If |γZ | = Θ(1), this gives the stated scale. The rank-one matrix containing only this Z row has a singular value and a gap to zero of order d−1/2 . For the complete cross-layer matrix, using this gap also requires control of RZ and separation from the other singular values. The calculation above does not supply those bounds. If the displayed contribution remains spectrally isolated at this scale and the fitted-source concentration bound holds, then the right feature dimension D0 ≍ d2 , (89) and ∆A ≍ d−1/2 give ncross ≫

D0 ≍ d3 . ∆2A

(96)

Hence the transition near α = 3 is naturally associated with learning the backward score rather than with the subsequent effective-label covariance. Corrected upper Forward pass. After the lower correction, a new upper Forward estimator acts on a representation of effective quadratic dimension D2 (d1 ) ≍ d21 ≍ d2ϵ . For a fixed number of retained upper directions and nonvanishing population gaps, its natural sample scale is therefore again d2ϵ . We leave the exact corrected upper rank unspecified until the corresponding numerical architecture is fixed. Scaling summary. Table 1 collects the resulting scales. For 0 < ϵ < 1, the learned-score route is therefore predicted to be controlled by n

o

nBackLoFi ≍ max d2+ϵ , d3 = d3 .

(97)

For the numerical choice ϵ = 1/2, the corresponding transitions are 5 αU = , 2

αBack = 3.

The d2+ϵ scales follow from the population gaps established above combined with the effectivedimension concentration rule (87). By contrast, the d3 prediction is based on the contraction size derived in (95). Its application to the complete fitted operator still requires control of the remaining contractions and the resulting singular-value gap. 45

Algorithm 5 One pass of Algorithm 2, as implemented. Require: Network (Fℓ , Rℓ , cℓ )ℓ≤L and readout (aL , b); rank fraction f , scale β; signal ∈ {first-order, exact}. 1: Compute zℓ , Dℓ by (98); gL ← DL aL ; keep a copy (F0ℓ , Rℓ0 ) of the network 2: for ℓ = L − 1, . . . , 1 do b n [y gℓ+1 z⊤ ]; s ← U⊤ gℓ+1 3: (U, V) ← top rℓ singular pairs of Aℓ = E ▷ frozen scores ℓ−1 4: I ← the rℓ columns of Fℓ of lowest importance (100); τℓ ← mean coordinate RMS of the kept fwd and read columns 5: Fℓ [:, I] ← vi βτℓ / RMS⟨vi , zℓ−1 ⟩; Rℓ [:, I] ← fresh N (0, 1) ▷ row step ⊤ b eℓ , Dℓ ; Q ← En [z eℓ (ys) ] 6: Recompute z 7: J ← the rℓ columns of Fℓ+1 of lowest importance; τℓ+1 as above √ eℓ ⟩; Rℓ+1 [:, J] ← U 8: Fℓ+1 [:, J] ← qi β pℓ+1 τℓ+1 / RMS⟨qi , z ▷ column step 9: if ℓ > 1 then
  b n 1 F0 (R0 )⊤ gℓ+1 + 1 Q0 s , q0 ∝ 1 R0 (F0 )⊤ vi 10: first-order: gℓ ← Dℓ E i ℓ+1 ℓ cℓ+1 ℓ+1 cℓ+1 cℓ ℓ ⊤ gℓ+1 through the corrected W exact: gℓ ← Dℓ Wℓ+1 ℓ+1 , per sample 12: end if 13: end for 14: Recompute zL ; refit (aL , b) by ridge regression

11:

F

The implemented algorithm

This appendix gives Algorithm 2 exactly as it is run for the experiments, with the choices that the main text leaves open. The values of every hyperparameter are listed in Appendix H. Parameterization. A layer is a filter Fℓ ∈ Rpℓ−1 ×Kℓ followed by a Gaussian lift Rℓ ∈ Rpℓ ×Kℓ with i.i.d. N (0, 1) entries, hℓ =

1 Rℓ F⊤ ℓ zℓ−1 , cℓ

σ(hℓ ) zℓ = √ , pℓ

Dℓ =

diag(σ ′ (hℓ )) , √ pℓ

b n ∥F⊤ z ∥2 , c2ℓ = E ℓ ℓ−1 2

(98)

so that Wℓ = Rℓ F⊤ ℓ /cℓ . The scale cℓ is fitted with the forward network and never refitted. A bℓ cutoff kℓ stores Kℓ = kℓ + 1 columns, the normalized first moment followed by kℓ eigenvectors of C deflated against it. Each column carries a tag, fwd, read or pair, according to whether it was installed by forward Neural LoFi, by a row step or by a column step. The rank of every boundary is rℓ = max{1, round(f k)} for a rank fraction f . The widths never change: every append replaces rℓ existing columns. Scales. The only free number is β, the size of an appended coordinate in units of a forward one. It sets the amplitude of every column appended by the row step and by the column step. b n ξ 2 )1/2 for the root mean square over the training set, and let τ be the mean of Write RMS ξ = (E ℓ RMS⟨fa , zℓ−1 ⟩ over the surviving columns fa of Fℓ that carry a random lift (fwd and read). The row step appends γirow vi to Fℓ and the column step appends γicol qi to Fℓ+1 , with √ β pℓ+1 τℓ+1 β τℓ row col γi = , γi = . (99) eℓ ⟩ RMS⟨vi , zℓ−1 ⟩ RMS⟨qi , z In the notation of (11), where the factor 1/c of (98) is absorbed, this reads Γrow = diag(γirow )/cℓ ℓ col and Γℓ = diag(γicol )/cℓ+1 . An appended read coordinate thus has root mean square βτℓ . A paired 46

√ column carries the extra factor pℓ+1 because its lift column ui has unit norm, while a forward column is lifted by N (0, 1) entries; with it, the paired column contributes to a neuron of layer ℓ + 1 as much as β forward coordinates. Each γi is set per column, which equalizes the amplitudes of the retained modes. The columns qi0 of the first-order signal are rescaled to the same target as γicol qi . We use β = 0.25 for the fully connected network and β = 4 for the convolutional one, each chosen by a local scan at n = 104 . Which columns are replaced. The widths are fixed, so each append replaces rℓ existing columns: the row step replaces the set I of columns of Fℓ , and the column step the set J of columns of Fℓ+1 . In both cases we remove the rℓ columns of lowest importance Ω, and the rule treats forward, read and paired columns alike at every pass. Importance is defined as in forward Neural LoFi, where a direction matters when gradient flow amplifies its coefficients in the rows, but it is computed on the current network, whose layers above are no longer random. For a column a of layer m it is built in three steps. b a be the unit direction of the column and ha = ⟨v b a , zm−1 ⟩ its coordinate. 1. Coefficients. Let v Neuron j reads it with the coefficient cj,a = (Rm )ja ∥va ∥/cm , and sees the effective label y ⟨ω j , gm+1 ⟩, where ω j is the column of Wm+1 that reads it. b a , with the label in place of the residual and the gate 2. Rate. Project the row update (4) on v b n hσ ′ (h)). This gives b n σ ′ (h), E expanded to first degree, σ ′ (h) ≃ dm + em h with (dm , em ) = (E ċj,a ∝ em λj (a) cj,a + dm ℓj (a), with the drivers λj , ℓj given in (100).

3. Importance. Take the root mean square of this rate over the neurons of the layer: b n [y g b h2a ] , λj (a) = ω j , E

b n [y g b ha ] , ℓj (a) = ω j , E

Ω(a) = RMSj em λj (a) cj,a + dm ℓj (a) .

(100)

A column whose coefficients no neuron would move can thus be removed at the least cost. In the rate, the term in λj is the covariance channel, multiplicative in the coefficient, and the term in ℓj b a by the the spike channel, additive; only the diagonal of the driver is kept, since the drive of v other directions a neuron reads averages out over neurons. The backward signal enters through b = ḡ + UU⊤ (gm+1 − ḡ), its mean plus its sample dependence along the retained modes of the g boundary above, the only part that reaches the neurons with a coherent sign. When the backward b n [yh2 ], and Ω orders the covariance b = ḡ, λj (a) is ⟨ω j , ḡ⟩ E signal carries no sample dependence, g a directions by the signed Rayleigh quotient that forward Neural LoFi uses, weighted by how strongly b ≡ 1. the neurons read them. For the top layer ω j = (aL )j and g Convolutional network. The network is conv1 (3 × 3, no filter), conv2 (3 × 3, filter kc over channels) with pooling, and a fully connected layer (filter k3 ). At the boundary conv2 | FC, A pairs the FC sensitivity with the spatial mean of conv1’s output; the row step replaces channel directions of conv2, each with a fresh random 3 × 3 kernel, and the column step replaces FC filter columns by Q computed on the√pooled map, with all H ′ W ′ lift rows of a paired channel equal to ui and the target divided by H ′ W ′ to match a forward channel. Conv1 has no filter, so the boundary conv1 | conv2 performs the column step only, with locations treated as samples and ui placed on the centre tap. The backward signal is carried through the pooling by its adjoint, and the importance is computed per filter and tap. Every statistic is streamed over the training set.

47

Numerics. The top rℓ singular pairs are computed exactly when rℓ + 32 exceeds half the smaller dimension of Aℓ and by a randomized SVD otherwise (oversampling 32, four power iterations, fixed seed). Statistics are chunked averages over the training set, accumulated in double precision.

F.1

Variants tested and not retained

We think it is useful to report approaches that did not work, so that others need not try them again. The following variants were implemented and run on the fully connected network, and none was retained. Gate-only signal (gℓ = Dℓ 1): no gain over forward Neural LoFi, as expected, since its readdirections stay in the span of the forward filter (Section 2). Mean-field and exact signals: the same test error as the first-order signal within seed noise; the first-order signal is kept because it needs no propagation through the installed blocks. √ Common column scale (the mean or the largest paired column set to β pℓ+1 τℓ+1 , relative amplitudes kept): no improvement over forward Neural LoFi; only the per-column scale did improve. √ Scale τ over all columns: paired columns carry p-larger coordinates, so including them inflates τ by a factor of about 6 per pass, and the later appends with it. b n [yh2 ]|): backward columns score high, so passes Weakest swap (the lowest Rayleigh quotients |E a accumulate them and consume the forward columns.

Replace-own swap (each pass rewrites the slots of the previous one): no gain beyond the first pass. eℓ before the swap): the same Update of the whole block ℓ + 1 (forward columns refitted on z test error, within 0.1 points at every n.

Dense interpolation (W ← (1 − α)W + α estimate, α = 0.1, no swap): the same test error, within 0.1 points at every n.

G

Additional numerical simulations

G.1

Number of passes at fixed width

Figure 4 shows the test error of the fully connected network against the number of passes of Algorithm 2 at fixed width. With rℓ = 5%, 7%, 10% and 12% of kℓ the error decreases with the passes to the same minimum, reached when about half of the filter columns have been replaced (pass 4 at 12%, pass 5 at 10%, pass 7–8 at 7%, pass 9–10 at 5%); further passes exhaust the forward columns and the error rises quickly. With rℓ = 3% the error is still decreasing slowly after ten passes, and with rℓ = 2% it stays at the forward level.

48

Test Error (%)

n = 1k r = 2% r = 3% r = 5% r = 7% r = 10% r = 12%

19 18

n = 10k

n = 5k

n = 50k

16 14

16 14

12

14

17 0

2

4

6

8 10

0

2

Pass

4

6

8 10

0

Pass

2

4

6

8 10

0

2

Pass

4

6

8 10

Pass

Figure 4: Test error of the fully connected network against the number of passes of Algorithm 2 at fixed width, for six rank fractions rℓ /kℓ and four training-set sizes n. Pass 0 is forward Neural LoFi (dotted line). Each point is the best filter width kℓ for its n, rank fraction and pass, selected by the mean over 10 seeds; bands are one standard error. Settings in Appendix H.

G.2

What successive passes select

Figure 5 shows, for each pass, the fraction of its new read and write directions that lies in the span of those installed by the earlier passes. For the read directions it stays above 0.97, against at most 0.02 for a random subspace: each pass selects again the same read subspace and installs it with fresh lifts. Successive passes therefore resample random features within the same subspace. The error in Figure 4 consequently depends on the fraction of columns replaced and is insensitive to the per-pass step rℓ . Boundary 2 | 3

Boundary 1 | 2

Fraction in span

1.0 0.8 0.6 read, r = 2% write, r = 2% read, r = 5%

0.4 2

write, r = 5% read, r = 10% write, r = 10%

3

4

5

6

Pass

2

3

4

5

6

Pass

Figure 5: Fraction of the read (full) and write (dashed) directions selected at a pass inside the span of those installed by the earlier passes, ∥Vt⊤ Q<t ∥2F /rℓ , at the two boundaries of the fully connected network (kℓ = 512, n = 10,000, mean and standard error over 3 seeds). A random subspace gives at most 0.02.

G.3

Damped updates

The fixed-width version of Algorithm 2 used in Figure 3 replaces at each pass the rℓ forward columns of lowest importance by the corrected ones, so that the passes consume the forward columns (Appendix G.1). Figure 6 tests a rule that consumes nothing: the layers are kept dense and the two

49

updates of (11) are applied as damped averages, ⊤ Wℓ ← (1 − α) Wℓ + α Rℓ← Γrow ℓ Vℓ ,

⊤ Wℓ+1 ← (1 − α) Wℓ+1 + α Uℓ Γcol ℓ Qℓ ,

(101)

with α = 0.1, the same estimates, no column replaced and the readout refitted after each pass. With rℓ = 25% of kℓ the error decreases monotonically over about 25 passes to 18.7%, 16.9%, 13.6%, 12.6% and 11.0% at n = 500, 1,000, 5,000, 10,000 and 50,000, within 0.2 points of the minimum of the swap (18.7%, 16.8%, 13.4%, 12.4% and 11.1%) and never significantly below it, then rises slowly, by 0.15–0.25 points at pass 50; with rℓ = 10% the descent is slower (minimum at passes 30–48 for n ≥ 5,000) and stops 0.2 points higher. The damped rule reaches essentially the same minimum error as the swap over more passes, followed by a gradual rise. This is consistent with
α the identity (1 − α)W + αE = (1 − α) W + 1−α E : the network is positively homogeneous and the ridge readout absorbs the factor 1 − α, so the damped rule is the update (11) with the scale β replaced by βα/(1 − α), applied at every pass. Both methods use the same direction-selection rule. damped, r = 10%

FwdLoFi

n = 1k

n = 500

n = 10k

n = 5k 14.5

19.5

Test Error (%)

swap (best r)

damped, r = 25%

n = 50k

13.5

17.5

12.0 14.0

19.0

13.0 11.5

17.0 12.5

13.5

18.5 0

25

Pass

50

0

25

50

11.0 0

25

Pass

Pass

50

0

25

Pass

50

0

25

50

Pass

Figure 6: Test error of the fully connected network against the number of passes of the damped update (101) with α = 0.1, at two rank fractions rℓ /kℓ , and of the fixed-width swap of Figure 3 at its best rank fraction (dashed), for five training-set sizes n. Pass 0 is forward Neural LoFi (dotted line). Each point is the best filter width kℓ for its n, rank fraction and pass, selected by the mean over 10 seeds; bands are one standard error. Settings in Appendix H.

G.4

Gradient-descent training curves

Figures 7 and 8 show the test error of full-batch gradient descent along training for the learning rates that do not diverge, on the fully connected and on the convolutional network. On the fully connected network the rate 0.1 is the best at every n. On the convolutional network the rate 0.003 is the best up to n = 5,000 and the largest stable rate, 0.005, from n = 10,000 upwards; the rate 0.001 is behind at every n. At large n the error is still decreasing at the last step on both networks, so the gradient-descent curves of Figure 3 are limited by the step budget there (3,000 steps on the fully connected network, 5,000 on the convolutional one).

50

BackLoFi

Test Error (%)

n = 1k

GD, lr = 0.01

GD, lr = 0.03

GD, lr = 0.1

n = 10k

n = 5k

n = 50k

30 20 25

20

15 15

20

15 101

102

103

10 101

Step

102

103

101

102

Step

103

101

Step

102

103

Step

Figure 7: Test error of full-batch gradient descent on the fully connected network against the step, for three learning rates; mean and standard error over 5 seeds. Dotted: BackLoFi of Figure 3. GD, lr = 0.001

BackLoFi

Test Error (%)

n = 1k

GD, lr = 0.003

GD, lr = 0.005

n = 10k

n = 5k

n = 50k

20

20

15

20

10 101

102

Step

103

10

10

10 101

102

103

101

Step

102

Step

103

101

102

103

Step

Figure 8: Test error of full-batch gradient descent on the convolutional network against the step, for three learning rates at every n; mean and standard error over 3 seeds for the rates 0.003 and 0.005, one seed for 0.001 at n ≥ 5,000. Dotted: BackLoFi of Figure 3.

H

Hyperparameters of the figures

The image experiments share one protocol: binary CIFAR-10 (vehicles +1, animals −1), a training sub-sample of size n, the 104 test images, a ridge readout with leave-one-out penalty among 50 log-spaced values in [10−6 , 106 ], test error as the fraction of wrong signs, means over seeds with every selection made on the seed mean. Algorithmic terms are defined in Appendix F. Figure 2. We use d = 100 and the finite-dimensional realization (pU , pV , pW ) = (6, 5, 5), corresponding to the ϵ = 1/2 scaling. The sixteen first-level quadratic directions are drawn from √a D (100) 2 Gaussian ensemble and jointly orthonormalized in R . For the second level,√BQ = I6 / 6, while BZ is a randomly rotated balanced contrast with three eigenvalues +1/ 6 and three √ eigenvalues −1/ 6. The hidden coupling C is drawn as a Haar orthogonal matrix scaled by √ 1/ 5. We use b0 = 1/2 and no additive label noise. Labels are centered and standardized using the training statistics, which are then kept fixed for validation and test. The spectral estimators follow Algorithm 4. The frozen lower random-feature map has width 65,536, unit-norm Gaussian weights, and activation σ≥2 (t) = ReLU(t) − 2t − √12π The lower filter retains 16 directions, of which the leading 6 define the protected visible anchor. The upper 51

Forward step uses the complete quadratic Hermite lift of these 16 coordinates, of dimension 136, and retains 16 directions. Within BackLoFi algorithm, after correction, the upper Forward step is rerun with the original labels. The final affine ridge readout is fitted on the linear and quadratic upper features, with relative regularization selected from {10−8 , 10−6 , 10−4 , 10−2 , 1} by validation MSE. For the gradient-descent baselines, we use the same frozen lower random-feature map followed by two trainable hidden layers of widths (512, 512) with ReLU activations and an affine scalar readout. We use AdamW with learning rate in {3 × 10−4 , 10−3 }, batch size min(n, 1024), weight decay 10−4 , gradient-norm clipping at 5, and 4000 updates. End-to-end GD trains both hidden layers jointly. For layer-wise GD, the budget is split into 2000 updates for the first layer and a temporary readout, followed by 2000 updates for the second layer and final readout with the first layer frozen. The best validation checkpoint and learning rate are selected separately for each sample size and seed. We use α ∈ {1.5, 1.6, . . . , 3.5} with n = round(100α ), two seeds {0, 1}, 16,384 validation samples and 32,768 test samples. Curves report the mean and one standard error over the two seeds. Figure 3. n ∈ {500, 1k, 2k, 5k, 10k, 20k, 50k}. Left: L = 3 fully connected ReLU layers, pℓ = 5,000, 10 seeds. FwdLoFi: kℓ ∈ {32, . . . , 2,048} (powers of two), best per n. BackLoFi (1st ord): fixed width, rℓ /kℓ ∈ {5, 10, 15, 20, 25}%, P ≤ 4 passes, β = 0.25, best (kℓ , rℓ , P ) per n; exact: rℓ /kℓ ∈ {10, 20, 25}%, kℓ ≥ 64. RF: identity filters. Ridge: on the 3,072 pixels. GD: same widths, full batch, MSE on ±1 labels, no momentum or weight decay, 3,000 steps, rates {0.01, 0.03, 0.1, 0.2}; curve at rate 0.1 (best on average over n) and its best recorded step, 5 seeds. Right: conv 3 × 3 with 4,096 channels (no filter), conv 3 × 3 with 4,096 channels (filter kc ), 4 × 4 max pooling, fully connected 5,000 (filter k3 ), ReLU, 5 seeds. FwdLoFi: kc ∈ {64, . . . , 1,024}, k3 ∈ {32, . . . , 1,024}. BackLoFi (1st ord): k3 ≤ 512, ranks round(f min(kc , k3 )) and round(f kc ) at the two boundaries, f ∈ {5, 10, 15, 20, 25}%, P ≤ 4, β = 4; exact: kc , k3 ∈ {128, . . . , 1,024}, f ∈ {10, 20, 25}%, P ≤ 3. RF: identity filters. GD: same widths, full batch, MSE on ±1 labels, no momentum or weight decay, 5,000 steps, bfloat16 autocast (a replica without it at n = 500 gives the same curve within the seed noise), rates {0.001, 0.003, 0.005} (0.01 diverges at the third step for every n); 3 seeds for 0.003 and 0.005 at every n and for 0.001 up to n = 2k, one seed above; the curve uses, at each n, the best of the rates run with three seeds, selected on the seed mean, at its best recorded step (0.003 up to n = 5k, 0.005 above). Shaded bands are one standard error over seeds. Backward signal: compared on the grid of the exact arm, the first-order signal of (12) and the exact signal backpropagated through the corrected network give the same error on the fully connected network within the standard error at every n; on the convolutional network the first-order signal is lower by about 0.3 points for 2,000 ≤ n ≤ 10,000 (0.4 at n = 500) and equal at n = 50,000. Reference values at n = 50,000: fully connected, RF 13.9%, FwdLoFi 12.3%, BackLoFi 11.1%, GD 9.5%; convolutional, RF 8.4%, FwdLoFi 7.8%, BackLoFi 6.7%, GD 5.7%. Figure 4. Fully connected network and BackLoFi (1st ord) of Figure 3, with rℓ /kℓ ∈ {2, 3, 5, 7, 10, 12}%, P = 1, . . . , 10, β = 0.25, kℓ ∈ {32, . . . , 2,048} (best per point), n ∈ {1k, 5k, 10k, 50k}, 10 seeds. Figure 5. Fully connected network, first-order signal, kℓ = 512, n = 10k, rℓ /kℓ ∈ {2, 5, 10}%, 6 passes, β = 0.25; 3 seeds. Figure 6. Fully connected network of Figure 3 with dense layers and the damped updates (101), first-order signal, α = 0.1, rℓ /kℓ ∈ {5, 10, 25}% (5% not shown: it stays within 0.1 points of 52

FwdLoFi for n ≤ 1k and descends more slowly than 10% above), P = 1, . . . , 50, β = 0.25, kℓ ∈ {32, . . . , 2,048} (best per point), n ∈ {500, 1k, 5k, 10k, 50k}, 10 seeds. The forward scale τ entering Γ is measured with the forward-phase filters on the current features; both updates are applied at each boundary in turn, so a hidden layer between two boundaries is damped twice per pass; a layer whose pre-activation scale drifts by more than a factor 8 would be rescaled by a power of two (never triggered). Swap reference: the 4-pass arm of Figure 3 at its best rℓ for each n. Figure 7. GD runs of Figure 3, rates {0.01, 0.03, 0.1}, n ∈ {1k, 5k, 10k, 50k}, 5 seeds. Figure 8. Convolutional network of Figure 3 with trainable layers of the same widths, random initialization, full batch (gradient accumulated over the whole training set), MSE on ±1 labels, no momentum or weight decay, 5,000 steps, bfloat16 autocast; rates {0.001, 0.003, 0.005} (0.01 diverges) at every n; seeds as in the caption.

53

Related documents

• Neural LoFi with Backward Coupling: A Spectral Theory of Cross-Layer Feature Learning#1448365
• How Finite-Rank Bottlenecks Shape the Low-Rank Adaptation Landscape#1387621
• Support vector machines for spectral unmixing#223453
• NObSP: Functional Decomposition of Neural Networks via Oblique Subspace Projections#1011540
• Don't Drop Dropout: Optimizing Layer Sparsity for Efficient LLM Training and Inference#660858