Paper deep dive
Retrieval-Centric Deep Learning in Growing Nonparametric Neural Networks
Maximilian Schlegel, Rajai Nasser, Seijin Kobayashi, Yanick Schimpf, Oliver Sieberling, Robert Obryk, Kazuki Irie, João Sacramento, Johannes von Oswald
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%
Last extracted: 10/6/2026, 10:24:11 AM
Summary
The paper introduces Retrieval-Centric Deep Learning (RCDL), a paradigm where neural networks grow by storing key-value pairs for each data point instead of compressing data into fixed-size weight matrices. It derives functional gradient-based learning rules for kernelized attention layers using Radial Basis Function (RBF) and weighted softmax kernels, addressing the non-trivial optimization challenges of replacing linear attention with more powerful attention mechanisms. The authors demonstrate RCDL's efficiency on image classification and teacher-student tasks and establish a formal connection between advanced linear attention variants (MesaNet/DeltaNet) and conventional optimizers.
Entities (8)
Relation Signals (7)
Retrieval-Centric Deep Learning → emplokeskernel → Radial Basis Function Kernel
confidence 95% · develop functional gradient-based learning rules for kernelized attention layers, based on radial basis function (RBF) and softmax-like kernels
Retrieval-Centric Deep Learning → emplokeskernel → Weighted Softmax Kernel
confidence 95% · we introduce a 'weighted softmax kernel' and derive a learning algorithm for the corresponding growing NN.
Retrieval-Centric Deep Learning → demonstratesperformanceon → Image Classification
confidence 92% · Empirically, we demonstrate the promising performance and learning-efficiency of RCDL on image classification
Retrieval-Centric Deep Learning → demonstratesperformanceon → Teacher-Student Learning
confidence 92% · Empirically, we demonstrate the promising performance and learning-efficiency of RCDL on ... synthetic teacher-student learning tasks.
MesaNet → isvariantof → Linear Attention
confidence 90% · advanced LA variants, namely MesaNet/DeltaNet
DeltaNet → isvariantof → Linear Attention
confidence 90% · advanced LA variants, namely MesaNet/DeltaNet
Retrieval-Centric Deep Learning → uses → Linear Attention
confidence 90% · While Irie et al. ... have put forward this perspective from the classic duality expressing any linear layer ... as linear attention (LA) ... replacing LA by more powerful attention functions
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:We investigate a general-purpose layer for deep learning that, instead of compressing arbitrary-size training data into fixed-size weight matrices, stores a new pair of key-value representations for every data point during training, and retrieves and recombines these representations through an attention mechanism at inference time - resulting in a growing neural net (NN). While Irie et al. (arXiv:2202.05798) have put forward this perspective from the classic duality expressing any linear layer in a deep NN trained by gradient descent as linear attention (LA) over the training data points, replacing LA by more powerful attention functions, as they suggest, turns out to be non-trivial: we show that naively applying learning rules from the LA case to advanced kernels does not lead to principled optimization. Here we fill this gap and develop functional gradient-based learning rules for kernelized attention layers, based on radial basis function (RBF) and softmax-like kernels - establishing the principled "retrieval-centric deep learning" (RCDL) paradigm. Empirically, we demonstrate the promising performance and learning-efficiency of RCDL on image classification and synthetic teacher-student learning tasks. Moreover, we show that replacing LA in the dual form of NNs by advanced LA variants, namely MesaNet/DeltaNet, yields a formal connection to recently proposed optimizers for conventional fixed-size NNs, offering a novel perspective on deep learning optimization.
Tags
Links
- Source: https://arxiv.org/abs/2610.03858v1
- Canonical: https://arxiv.org/abs/2610.03858v1
Trouble viewing inline? Open PDF directly →
Full Text
196,174 characters extracted from source content.
Expand or collapse full text
† *All authors contributed equally. Correspondence: schlegelm,rajainasser,jvoswald@google.com, osieberl@mit.edu, kazuki.irie@yale.edu. Public code: https://github.com/paradigms-of-intelligence/rcdl Retrieval-Centric Deep Learning in Growing Nonparametric Neural Networks Maximilian Schlegel Rajai Nasser Seijin Kobayashi Yanick Schimpf Affiliation: Google, Paradigms of Intelligence Team, Zurich, Switzerland Oliver Sieberling Robert Obryk Affiliation: Google, Paradigms of Intelligence Team, Zurich, Switzerland Affiliation: MIT CSAIL, Cambridge, MA, USA Kazuki Irie João Sacramento Johannes von Oswald Affiliation: Google, Paradigms of Intelligence Team, Zurich, Switzerland Affiliation: Department of Computer Science and Wu Tsai InstituteYale University, New Haven, CT, USA Abstract We investigate a general-purpose layer for deep learning that, instead of compressing arbitrary-size training data into fixed-size weight matrices, stores a new pair of key-value representations for every data point during training, and retrieves and recombines these representations through an attention mechanism at inference time—resulting in a growing neural net (N). While Irie et al. [33] have put forward this perspective from the classic duality expressing any linear layer in a deep N trained by gradient descent as linear attention (LA) over the training data points, replacing LA by more powerful attention functions, as they suggest, turns out to be non-trivial: we show that naively applying learning rules from the LA case to advanced kernels does not lead to principled optimization. Here we fill this gap and develop functional gradient-based learning rules for kernelized attention layers, based on radial basis function (RBF) and softmax-like kernels—establishing the principled “retrieval-centric deep learning” (RCDL) paradigm. Empirically, we demonstrate the promising performance and learning-efficiency of RCDL on image classification and synthetic teacher-student learning tasks. Moreover, we show that replacing LA in the dual form of NNs by advanced LA variants, namely MesaNet/DeltaNet, yields a formal connection to recently proposed optimizers for conventional fixed-size NNs, offering a novel perspective on deep learning optimization. Table 1: Comparison of parametric vs. retrieval-centric deep learning paradigms. Model Parametric neural networks Growing nonparametric networks Building block Linear layer Attention layer Learnable state Fixed-size weight matrix Growing database of keys & values (for t train examples) ∈ℝdout×din W ^d_out× d_in =∈ℝdin×t,∈ℝdout×tD=\ K ^d_in× t, V ^d_out× t\ Inference = y= W x =s()=Attention(,,) y=s_D( x)=Attention( K, V, x) (∈ℝdin x ^d_in, ∈ℝdout y ^d_out) Illustration Forming deep models Compose linear layers and nonlinearities, e.g., σ(2σ(1))σ( W_2σ( W_1 x)) Compose attention layers and nonlinearities, e.g., σ(s2(σ(s1())))σ(s_D_2(σ(s_D_1( x)))) Learning objective Minimize a global loss ℒL (e.g., via gradient descent) Learning update Update the fixed weight matrix Append a new key-value pair (e.g., gradient-based) (e.g., Δ∝−∇ℒ=−∇tℒt⊤ W - _ WL=- _ y_tL\, x_t ) (e.g., t=t k_t= x_t, t∝−∇tℒ v_t - _ y_tL) Key question How to update the weight matrix? What key-value vectors to append? 1 Introduction Traditional deep learning uses gradient descent to compress arbitrary-size training datasets and/or experiences into fixed-size weight matrices of a deep neural network (N). A classic but lesser-known mathematical property thereof is that any linear layer in such a gradient descent-trained deep N admits a nonparametric dual form: Its exact input/output relation can also be achieved by an equivalent “key-value memory system” that stores all the training datapoints presented to the layer during training and the initial weights, and produces outputs through retrieval using “unnormalized dot-product attention” over the entire training experience stored in the form of “key/value pairs” [1, 33] (reviewed in Sec. 2). From a modern sequence modeling research perspective, however, such an unnormalized dot-product attention represents a form of “linear attention” [40, 62], which is well-known to be sub-optimal to process sequences, e.g., compared with the softmax attention [4] commonly used in the transformer [69]. This observation has raised a natural question: can we achieve more powerful deep learning systems by replacing the linear attention (LA) underlying deep NNs by a more powerful attention function [33, 37]? While this question has been left unanswered, the time is ripe to revisit it, especially given how recent research in sequence modeling has advanced LA and provided insights to interpret and improve LA through various angles, including its connection to fast weight programming and learning rules [60, 59, 35, 34, 36, 79, 77, 63], the underlying local optimization perspectives [70, 71, 72, 67, 6, 8, 7, 73], as well as traditional sequence processing tools such as gating [51, 78]. Here we pursue two approaches that aim to go beyond current deep learning by extending the LA mechanism underlying any modern deep NNs’ dual form. First, we advance the retrieval-centric view of the dual form by developing a theory for the “growing neural network” obtained by introducing a (non-linearizable) kernel to the dual form’s LA, and devise its learning algorithms (which prescribe the key and value representations to be stored for each data point) that optimize a given global objective function. We study the cases of positive definite kernels (e.g., the classic radial basis function (RBF) [61]) and the “weighted softmax” kernel we propose. In fact, while Irie et al. [33] have suggested to directly “insert” the softmax function in the corresponding linear attention to enhance retrieval from NNs’ training experience, devising a principled learning algorithm for such a system turns out to be highly non-trivial, as such an insertion invalidates the optimality of the key and value representations computed in the linear case. This is especially true for the softmax, for which we introduce a closely related “weighted softmax kernel” to formally derive a learning rule for the resulting retrieval-centric growing N. We conduct experiments using classic image classification datasets and teacher-student learning tasks, and demonstrate promising results in performance and sample efficiency compared to the standard deep NNs with fixed-size weight matrices. The second extension we pursue is to leverage recent advances in linear attention, namely the weight updates of DeltaNet [59, 77], and MesaNet [72]. In contrast to the first direction above, this approach gives rise to optimizers that improve conventional weight-based learning. We derive a formal connection to a recently proposed optimizer and second-order optimization, and empirically demonstrate the efficiency of the resulting optimizers, with competitive performance to Muon [39]. Finally, as a new framework that represents a departure from standard deep learning (see comparison in Table 1), retrieval-centric deep learning (RCDL) still faces many challenges. Perhaps the most central one is that model size grows at every training step when a nonlinear kernel is chosen. We finish our paper by presenting first results on a hybrid learning algorithm which manages capacity allocation dynamically, by combining standard (parametric) gradient updates with model growing (non-parametric) steps, triggered only when the inputs to a given layer are sufficiently distinct from past ones. We show empirically that this hybrid algorithm is able to match the performance of the non-parametric learner at a small fraction of its size. 2 Background: Dual Form of a Linear Layer Here we briefly review the dual form of a linear layer trained by gradient descent as a linear attention system [1, 33], which is our starting point. Let t, T, dind_in, doutd_out denote positive integers. Consider a linear layer inside an arbitrary complex/deep neural network, that transforms a layer input ∈ℝdin x ^d_in into a layer output ∈ℝdout y ^d_out as = y= W x where ∈ℝdout×din W ^d_out× d_in is a weight matrix. Given an arbitrary differentiable loss function ℒL and a sequence (1,…,T)( x_1,…, x_T) of T training inputs presented to the layer, with t∈ℝdin x_t ^d_in for 1≤t≤T1≤ t≤ T, we denote the corresponding (backpropagated) error signals as (1,…,T)( v_1,…, v_T) with t∈ℝdout v_t ^d_out defined by gradient descent, that is t=−ηt∇tℒ v_t=- _t _ y_tL where ηt∈ℝ>0 _t _>0 is the learning rate and ∇tℒ _ y_tL is the gradient of the loss w.r.t. the layer output t y_t corresponding to input t x_t. At training step t, gradient descent updates the weight matrix (WM), initialized with 0 W_0, as t=t−1−ηt(∇tℒ)⊗t=t−1+t⊗t W_t= W_t-1- _t( _ y_tL) x_t= W_t-1+ v_t x_t, where ⊗ denotes the outer product, that is t⊗t=tt⊤ v_t x_t= v_t x_t . We can unfold this equation to express the WM at step T as: T=0+∑t=1Ttt⊤ W_T= W_0+ _t=1^T v_t x_t . Thus, the layer output for an arbitrary input ∈ℝdin q ^d_in can be expressed as: =T=0+∑t=1Ttt⊤=0+∑t=1Ttϕ(t,) y= W_T q= W_0 q+ _t=1^T v_t x_t q= W_0 q+ _t=1^T v_tφ( x_t, q) (1) where the last term is a weighted sum of vectors t∈ℝdout v_t ^d_out with weights t⊤∈ℝ x_t q , which can be interpreted as an (unnormalized) “attention mechanism” by identifying t v_t as values, t x_t as keys, and q as the query; and ϕ(t,)=t⊤φ( x_t, q)= x_t q as attention weights [46, 69, 4]. To be more specific, this is (a simple form of) linear attention [40, 62]. This implies that the input/output relation of the linear layer can also be realized by a “growing nonparametric neural network”: Instead of storing the fixed-size weight matrix W as the system’s training memory, the system maintains a growing database by storing a new pair (tCLOSE( k_t, OPENt) v_t) of key-value vectors for each training data point—the layer input t=t∈ℝdin k_t= x_t ^d_in as the “key” and the corresponding error signal t=−ηt∇tℒ∈ℝdout v_t=- _t _ y_tL ^d_out as the “value” vector (in addition to the initial weight matrix 0 W_0). Inference is performed through the linear attention mechanism above (Eq. 1). Given the ubiquity of linear layers in deep neural networks, this dual view suggests that a fundamental mechanism linking training and inference in virtually any deep learning systems is based on (a cascade of) this simple linear attention mechanism. However, recent research on linear attention in sequence modeling indicates that such “simple linear attention” mechanisms are highly suboptimal. This points to the potential for fundamentally improving deep learning by enhancing the underlying linear attention mechanism. 3 Growing Nets with Non-Linear Attention With the goal of developing a general-purpose deep learning layer that can improve upon the conventional linear layer, we replace ϕφ in Eq. 1 by a more powerful kernel that is not reducible to a finite-dimensional linear map. Consequently, the resulting layer can no longer be converted into a linear layer with a fixed-size weight matrix. Instead, we obtain a “growing layer” in which a new pair of key/value representations is added to the system for each training datapoint. This raises a fundamental question: how should these key and value representations be chosen to implement a principled learning algorithm that minimizes a given loss function? Here we formally derive such algorithms in two classes of powerful attention functions. 3.1 Positive Definite Kernel Attention We first investigate the case where ϕφ in Eq. 1 is a positive definite kernel (PDK). By the Moore-Aronszajn theorem [2], any PDK admits a representation as an inner product of feature embeddings into a Hilbert space ℋϕH_φ. That is, there exists a map ψϕ:ℝdin→ℋϕ _φ:R^d_in _φ such that ϕ(,′)=⟨ψϕ(),ψϕ(′)⟩ℋϕφ( x, x )= _φ( x), _φ( x ) _H_φ. This allows us to view the nonparametric layer as a linear layer operating on the implicit feature space ℋϕH_φ with an infinite-dimensional weight matrix ∈ℝdout×∞ W ^d_out×∞. In Appendix A, we rigorously show that performing standard gradient descent on this implicit weight matrix is mathematically equivalent to appending a new key-value pair to the database at every step. Analogously to the analysis of the linear kernel case in Sec. 2, the gradient update of the implicit weight matrix is Δ∝(−ηt∇tℒ)ψϕ(t)⊤ W (- _t _ y_tL ) _φ( x_t) . This translates exactly to adding a new term (−ηt∇tℒ)ϕ(t,⋅) (- _t _ y_tL )φ( x_t,·) to the layer’s output function. Consequently, to build a growing layer based on a PDK, we retain the exact key-value definitions from the linear case: setting the key t=t k_t= x_t and the value t=−ηt∇tℒ v_t=- _t _ y_tL implements a principled loss-minimizing learning rule in the RKHS. In our experiments, we use the radial basis function (RBF) kernel—a classic PDK [61]. 3.2 Softmax and Weighted Softmax Attention While PDKs offer a robust theoretical foundation, modern deep learning relies heavily on softmax attention, which empirically outperforms unnormalized kernel attention in many sequence modeling tasks [69]. Here we study the case where we replace ϕφ in Eq. 1 by the softmax function (more formally, ϕ(t,)φ( x_t, q) becomes ϕ(t,,τ1≤τ≤T)φ( x_t, q;\ x_τ\_1≤τ≤ T) as it involves normalization). As it turns out, theoretically dealing with the standard softmax attention is highly non-trivial (as we discuss below; and in Appendix B.3). Instead, we introduce a “weighted softmax kernel” and derive a learning algorithm for the corresponding growing N. Let ⊂ℝdin×ℝdoutD ^d_in×R^d_out be a finite set of key-value pairs or “database”. This database induces the standard softmax attention function s:ℝdin→ℝdouts_D:R^d_in ^d_out defined as: s()=∑(,)∈e⟨,⟩/∑(,)∈e⟨,⟩.s_D( q)= _( k, v) e k, q v / _( k, v) e k, q \,. Consider the objective of making a small “corrective” adjustment to D to move s_D in a direction that minimizes a global loss. Ideally, we would achieve this by adding a new key-value pair to the database, analogous to the “growing layers” derived in Sec. 3.1. However, this is not possible in an infinitesimal manner for this softmax attention. Introducing any new key discontinuously alters the denominator, regardless of how small the associated value vector is. To enable the “smooth addition” of new key-value pairs, we introduce a scalar weight for each entry. This leads to the “weighted softmax attention” formulation, where the database D consists of a finite set of triples (w,,)∈ℝ+×ℝdin×ℝdout(w, k, v) ^+×R^d_in×R^d_out, and the retrieval function becomes: s()=∑(w,,)∈we⟨,⟩∑(w,,)∈we⟨,⟩.s_D( q)= _(w, k, v) we k, q v _(w, k, v) we k, q \,. (2) In the remainder of this section, we develop a theoretical framework that allows us to derive a “growing database” formulation which corresponds to performing a gradient descent in a space containing all weighted softmax functions. Now consider an infinite database containing a weight for every possible key-value pair. We formalize this by defining a weight function w:ℝdin×ℝdout→ℝ+w:R^d_in×R^d_out ^+, which induces the continuous softmax function sw:ℝdin→ℝdouts_w:R^d_in ^d_out as: sw()=∫ℝdin×ℝdoutw(,)e⟨,⟩∫ℝdin×ℝdoutw(,)e⟨,⟩.s_w( q)= _R^d_in×R^d_outw( k, v)e k, q v\,d k\,d v _R^d_in×R^d_outw( k, v)e k, q \,d k\,d v\,. (3) Appendix B.1 provides further remarks on the functional space of weighted softmax functions. The case of Gaussian mixtures. We analyze the case where w:ℝdin×ℝdout→ℝ+w:R^d_in×R^d_out ^+ represents the density of a Gaussian mixture. Let =(wi,i,i):1≤i≤n⊂ℝ+×ℝdin×ℝdoutD=\(w_i, k_i, v_i):1≤ i≤ n\ ^+×R^d_in×R^d_out be a finite set of n weighted key-value pairs, and let σk,σv>0 _k, _v>0. We define the Gaussian mixture density w_D as follows: w(,)=∑i=1nwi(2π)din+dout2σkdinσvdoute−‖−i‖22σk2−‖−i‖22σv2.w_D( k, v)= _i=1^n w_i(2π) d_in+d_out2 _k^d_in _v^d_oute^- \| k- k_i\|^22 _k^2- \| v- v_i\|^22 _v^2\,. This formulation corresponds to mixing n Gaussian densities centered at (i,i)1≤i≤n( k_i, v_i)_1≤ i≤ n with respective mixing weights (wi)1≤i≤n(w_i)_1≤ i≤ n. The covariance matrix Σ for each of the n Gaussian components is a block diagonal matrix given by: Σ=[Σk00Σv]=[σk2din00σv2dout]. = bmatrix _k&0\\ 0& _v bmatrix= bmatrix _k^2 I_d_in&0\\ 0& _v^2 I_d_out bmatrix\,. (4) Now, let us compute the softmax function sws_w_D induced by this mixture density. Substituting w_D into Eq. 3, we can show that (see the full derivation in Appendix B.2): sw() s_w_D( q) =∫ℝdin×ℝdoutw(,)e⟨,⟩∫ℝdin×ℝdoutw(,)e⟨,⟩ = _R^d_in×R^d_outw_D( k, v)e k, q v\,d k\,d v _R^d_in×R^d_outw_D( k, v)e k, q \,d k\,d v =∑i=1nwie⟨i,⟩+12σk2‖2i∑i=1nwie⟨i,⟩+12σk2‖2=∑i=1nwie⟨i,⟩i∑i=1nwie⟨i,⟩=s(). = _i=1^nw_ie k_i, q + 12 _k^2\| q\|^2 v_i _i=1^nw_ie k_i, q + 12 _k^2\| q\|^2= _i=1^nw_ie k_i, q v_i _i=1^nw_ie k_i, q =s_D( q)\,. This demonstrates that the continuous softmax function induced by the Gaussian mixture density w_D is exactly equivalent to the discrete weighted softmax function s_D defined by the database D. The shared variance-dependent term e12σk2‖2e 12 _k^2\| q\|^2 cancels out in the numerator and denominator, rendering the prediction independent of the specific variances σk _k and σv _v used in the mixture model. It is worth noting that the same cancellation holds for any general block-diagonal covariance matrix Σ (where k and v are independent), regardless of the specific correlations within the key or value dimensions. The functional space of Gaussian-modulated softmax functions. Now observing that Gaussian mixtures with a fixed covariance matrix Σ are absolutely continuous with respect to the centered Gaussian measure μ=(,Σ)μ=N(0, ), we can analyze the system using this background measure. We define the Gaussian-modulated softmax function for a relative weight function w as: sw()=∫ℝdin×ℝdoutw(,)e⟨,⟩μ(,)∫ℝdin×ℝdoutw(,)e⟨,⟩μ(,).s_w( q)= _R^d_in×R^d_outw( k, v)e k, q v\,dμ( k, v) _R^d_in×R^d_outw( k, v)e k, q \,dμ( k, v)\,. (5) This formulation relates to the original Lebesgue-density formulation via a Radon-Nikodym derivative. If the physical density of the database is p(,)p( k, v), the corresponding functional weight is w(,)=dpdμ(,)w( k, v)= dpdμ( k, v). This perspective reveals a crucial structural property: a Gaussian mixture in the physical space corresponds to a linear combination of exponential functions in the functional space. Specifically, the ratio of a shifted Gaussian (|,Σ)N(x|u, ) to the base measure (|,Σ)N(x|0, ) is proportional to e⟨,Σ−1⟩e , ^-1u . Consequently, adding a key-value pair to the database is mathematically equivalent to adding an exponential basis function to the weight w(,)w( k, v). This justifies searching for updates in the span of exponential functions during gradient descent. Let ℋσk,σvH_ _k, _v be the space of functions w:ℝdin×ℝdout→ℝw:R^d_in×R^d_out such that (,)∼(,Σ)[w(,)2 _( k, v) (0, )[w( k, v)^2 e⟨,⟩+⟨,′⟩]<∞,∀∈ℝdin,∀′∈ℝdout. e k, q + v, v ]<∞\,, ∀ q ^d_in\,,∀ v ^d_out\,. We endow ℋσk,σvH_ _k, _v with the L2(μ)L_2(μ) inner product: ⟨w,w′⟩ℋσk,σv=(,)∼(,Σ)[w(,)w′(,)]. w,w _H_ _k, _v=E_( k, v) (0, )[w( k, v)w ( k, v)]\,. Functional gradient descent update rule. We derive the update rule for the weight function w∈ℋσk,σvw _ _k, _v by performing functional gradient descent on a loss function ℒ(sw(),)L(s_w( q), y). For a single observation (t,t)( q_t, y_t), we compute the Fréchet derivative of the loss functional J(w)=ℒ(sw(t),t)J(w)=L(s_w( q_t), y_t) and identify the functional gradient ∇wJ _wJ via the Riesz representation theorem in ℋσk,σvH_ _k, _v. As detailed in Appendix C, the gradient ∇wJ _wJ takes a specific form involving the exponential kernel e⟨,t⟩e k, q_t and the error signal: ∇wJ(,)=e⟨,t⟩D(w)⟨∇^ℒ,−sw(t)⟩withD(w)=∫ℝdin×ℝdoutw(,)e⟨,t⟩μ(,). _wJ( k, v)= e k, q_t D(w) _ yL, v-s_w( q_t) \,\,with\,\,D(w)= _R^d_in×R^d_outw( k, v)e k, q_t \,dμ( k, v)\,. (6) D(w)D(w) is the denominator of the softmax functional induced by w. The functional update wt+1←wt−η∇wJw_t+1← w_t-η _wJ corresponds to adding a new term to the weight function, which is constructively equivalent to appending a new key-value tuple to the database (cf. Sec. C.4). Specifically, this establishes the update rule that adds a tuple (wnew,new,~new)(w_new, k_new, v_new) where the new key is a scaled version of the query, new=Σkt=σk2t k_new= _k q_t= _k^2 q_t, and the weight and value are determined by the prediction error and the local density. To ensure numerical stability, we store the weighted value vector ~=w v=w v, leading to the robust update rules: wnew∝η⟨t,t⟩,~new∝−ηΣvt∝−ησv2t,w_new η s_t, g_t , v_new -η _v g_t -η _v^2 g_t\,, (7) where t=∇^ℒ g_t= _ yL and t=swt(t) s_t=s_w_t( q_t). Refer to Appendix C for more details. In the case of mini-batch stochastic gradient descent with a batch of size B, the linearity of the gradient operator implies that the total functional gradient is the sum of individual gradients centered at distinct queries. Consequently, a single step corresponds to adding B distinct tuples to the database simultaneously, reflecting the non-parametric nature of the approach where model complexity grows with data. Remark 3.1. Given that new=σk2t k_new= _k^2 q_t, one can interpret σk2 _k^2 as the inverse temperature of the resulting softmax function. In practice, we can store new=t k_new= q_t and instead use a temperature hyperparameter T in the Softmax, as is typically done11 1 σk2 _k^2 relates to the temperature because increasing σk2 _k^2 scales the positions of new keys quadratically, causing them to spread apart faster than the linear expansion of the standard deviation of individual Gaussian components. Consequently, higher variance σk2 _k^2 acts as a separation parameter, reducing component overlap and sharpening the attention distribution., i.e., the final weighted-softmax function becomes: s()=∑(w,,~)∈e1Tdin⟨,⟩~∑(w,,~)∈e1Tdin⟨,⟩w.s_D( q)= _(w, k, v) e 1T d_in k, q v _(w, k, v) e 1T d_in k, q w\,. (8) In Appendix D, we provide further remarks regarding extensions of the learning algorithm above (corresponding to vanilla gradient descent) with weight decay and momentum. 3.2.1 Handling negative weights: clipping and gauge invariance One challenge with the update rule of equation 7 is that when ⟨t,t⟩<0 s_t, g_t <0 we get a negative weight wnew<0w_new<0. This may cause denominator cancellation and numerical instability (cf. Sec. C.5). A simple heuristic to enforce non-negativity is via clipping, wnew←max(0,wnew)w_new← (0,w_new). Another method for addressing the negative weights issue without clipping is via a translational symmetry of the softmax functional, which we call gauge-invariance (cf. Sec. C.7): If we shift all the value vectors in the database by ++ b, the output of the weighted softmax also shifts by ++ b. Hence, if we then subtract b from the output, we recover the original weighted softmax. Therefore, by equipping the softmax functional with an auxiliary global bias b, the model gains a translational symmetry between the values and the bias (cf. Proposition C.3). This means that the same function has a different representation for each offset (gauge). We can dynamically choose an offset frame—a “virtual gauge” t b_t along the gradient direction—that guarantees that the effective projection ⟨t,t+t⟩ g_t, s_t+ b_t is strictly positive, ensuring wnew>0w_new>0. Refer to Sec. C.7 for a detailed description of this method. 3.3 Sublinear Growth with Hybrid Non-Parametric and Parametric Learning Appending a new key-value tuple at every training step incurs a linear memory growth (T)O(T) with respect to the number of training steps T. To achieve scalable, sublinear growth, we introduce a hybrid learning rule that dynamically balances non-parametric expansion with parametric refinement: We treat the existing key-value pairs in the database as learnable parameters optimized via standard optimizers (e.g., Adam or AdamW), while conditionally growing the database using functional gradient descent only when a newly encountered key x is sufficiently distinct from the existing keys j k_j in the database, based on a threshold parameter τ, i.e., minj∥−j∥2>τ _j x- k_j _2>τ. By restricting memory growth to genuinely novel data while leveraging parametric optimization for representation tuning, this hybrid mechanism dramatically curbs memory footprint while preserving the performance advantages of the fully non-parametric model (see Figure 5). 4 Advanced Linear Attention Perspectives Sec. 3 focused on improving the LA underlying the dual form of NNs (Eq. 1) by introducing nonlinearity. Here we plug advanced LA models into it, and show how they lead to improved N optimizers. We focus on the mesa-layer [71, 72], which replaces the linear attention of Eq. 1 by the best possible linear key-value map, in the sense that it achieves the lowest (regularized) squared error over the keys and values seen up to time t: t W_t =argmin12∑τ=1t‖τ−τ‖2+λ2‖F2, = _ W 12 _τ=1^t\| v_τ- W x_τ\|^2+ λ2\| W\|^2_F\,, (9) where λ>0λ>0 is a static (scalar) hyperparameter which controls the strength of regularization, given by the squared Frobenius norm ‖F2\| W\|^2_F of the weights. Note that the mesa-layer considers the cumulative (in time) loss, and fully minimizes it. von Oswald et al. [72] empirically found that the mesa-layer outperforms both standard LA and the DeltaNet on various sequence modeling problems. Interestingly, as we show next, the dual view of a linear layer trained by gradient descent yields an approximate second-order optimizer once we replace linear attention (Eq. 1) by the mesa-layer (Eq. 9). The only remaining step is to determine what the values t v_t should be. In this section, we set them to be the gradient-corrected layer outputs: t=t−1t−ηt∇tℒ. v_t= W_t-1 x_t- _t _ y_tL. (10) An important property of the objective of the mesa-layer (Eq. 9) is that it can be solved analytically. We use this property to obtain our optimization algorithm, which we develop for the mini-batch setting—the current standard in deep learning. Before proceeding, we introduce additional notation to handle mini-batches. We denote by T the number of batches seen so far, each containing B data points. In this setting, T X_T is the din×Bd_in× B matrix that holds the current batch inputs in its columns, and T V_T is the dout×Bd_out× B matrix that holds the corresponding gradient-corrected outputs in its columns, following Eq. 10. Additionally, we define <T X_<T and <T V_<T as the history of all inputs and gradient-corrected outputs (resp.) observed thus far, excluding the current batch. Solving Eq. 9 analytically, we obtain a closed-form expression for the weights at iteration T i.e., T=(<T<T⊤+TT⊤)⋅(<T<T⊤+TT⊤+λ)−1. W_T=( V_<T X_<T + V_T X_T )·( X_<T X_<T + X_T X_T +λ I)^-1. It turns out that the resulting layer implements exactly a recently proposed optimizer [9, FOOF;], which has been shown to be highly efficient in deep learning. FOOF was derived as a simplification of KFAC [47], itself an approximate second-order optimization method. Intuitively, the optimization speed-up that FOOF (here, the mesa-layer) enjoys stems from the additional inverse matrix which is not present in linear attention. This additional factor applies a whitening transform to the inputs, which yields a full (second-order) Newton update for a shallow linear model, and an approximate second-order update for a deep neural network. We refer to Benzing [9] for an analysis of this optimizer, and to Appendix E for additional details. As an additional reference, we further consider the DeltaNet [59] case. As derived from the delta rule [75], this case corresponds to performing an online gradient descent step w.r.t. W on the layer-local, time-dependent loss penalizing the discrepancy between the prediction t W x_t and a target value t=−η∇tℒ v_t=-η _ y_tL: ℓt():=12‖t−t‖22=12‖t+η∇tℒ‖22 _t( W):= 12\| W x_t- v_t\|_2^2= 12\| W x_t+η _ y_tL\|_2^2 (11) with learning rate λ (corresponding to an online version of the mesa-layer case). This yields the update t=t−1+λ(−η∇tℒ−t−1t)t⊤ W_t= W_t-1+λ(-η _ y_tL- W_t-1 x_t) x_t . Appendix E provides additional details. 5 Experiments Here we conduct experiments to empirically validate theories arising from the two extended views of linear layers above: the newly derived growing neural networks and their learning rules (Sec. 3) and advanced optimizers (Sec. 4). Basic Settings. The baseline architecture is a residual MLP similar to the transformer’s MLP block [69]; the model has multiple layers with a block structure, where every block of the model updates the residual stream ∈ℝde e ^d_e following ←+2σ(1RMSNorm()), e← e+ W_2σ( W_1RMSNorm( e)), (12) with nonlinearity σ(⋅)σ(·) and RMSNorm prenormalization [81]. The first block is applied to the net input ∈ℝdx x ^d_x. The output layer is the standard softmax layer for classification or a linear layer for regression, using the cross-entropy or mean-squared error loss, respectively. To find a strong baseline, we perform intensive search over the choice of activation σ(⋅)σ(·), namely: ReLU [23, 50], GELU [31] and softmax [13]; as well as over the space of optimizers, including plain SGD, Adam [41], AdamW [45], RMSProp [32], and Muon [39]. These serve as the baselines for the MesaNet/DeltaNet-derived optimizers introduced in Sec. 4. We perform a fine-grained scan over hyperparameters, including model depth and nonlinearity type. All our results are averaged over 5 seeds, and shaded areas reflect stdev. unless noted otherwise (see Appendix F.1 for more details). Our nonlinear attention-based growing nets (Sec. 3) follow the same residual architecture, but perform: ←+g(RMSNorm()), e← e+g_D(RMSNorm( e)), (13) where g(⋅)g_D(·) can either be the weighted-softmax layer (Eq. 2) or the RBF layer (Eq. 1, with ϕφ an RBF). Figure 1: Teacher-student regression loss (mean ± s.e.m.) comparing: (blue, Weighted Softmax) our growing weighted softmax model; (red, Parametric Softmax) a very wide MLP baseline with width matching the final number of training samples using the softmax activation; and (gray, Reference MLP) an overparameterized student with matched architecture providing a lower bound. The final output layer is similar but without the residual connection. Importantly, note that the formal learning rules we derived for these growing nets correspond to training them using plain SGD (Sec. 3); derivations for more sophisticated optimizers are left for future work. Teacher-Student Experiments. We begin with experiments on a synthetic regression problem, where scalar target outputs y∗∈ℝy^* are generated by a held-out residual MLP teacher y∗()y^*( x) from 100-dimensional inputs ∈ℝ100 x ^100 [24]. The teacher N has 1 hidden layer with ReLU non-linearity, residual-stream de=100d_e=100, and random, fixed weights sampled from a LeCun normal distribution (fan-in). The teacher then projects the residual stream activation to an output dimension of 11. The student y()y( x) is trained through online SGD on the expected squared error loss ∼p()[(y∗()−y())2]E_ x p( x)[(y^*( x)-y( x))^2], by drawing new data samples at every training iteration. The inputs x are generated by first sampling from a standard normal distribution, and then projecting to the unit sphere. Appendix F provides further setup details. In this task, an overparameterized student—with the same architecture as the teacher, but larger width—is expected to perform very strongly [64]. In Fig. 1 we see that this is indeed the case; the MLP student (denoted by “Reference MLP” in the plot) outperforms all other models, and serves as the gold standard for this task. We then compare our Growing Weighted Softmax model (GWS; Sec. 3.2) with a fixed-size “parametric softmax” baseline that uses softmax activation (Eq. 12) and dimensions matching the full database size—closely mirroring the GWS model by treating 1 W_1 and 2 W_2 as learnable keys and values [26]. We compare a single-layer GWS with a single-block parametric softmax baseline (with two trainable weight matrices; the best model trained with Adam). We find that GWS significantly outperforms the parametric baseline. Despite the architectural mismatch with the teacher, and despite only having one trainable layer (comprising values and respective weights), our growing weighted softmax model can reach significantly lower loss values. Figure 2: (Left Plots) Test set accuracy on CIFAR-10/100: Growing Weighted Softmax & RBF, while based on plain SGD, significantly outperform the SGD-trained MLP and Adam/RMSProp and are competitive with MLPs trained with Muon. The Mesa/FOOF optimizer achieves high accuracy on CIFAR-10 and performs best on CIFAR-100. (Middle Right) Training loss for Growing Weighted Softmax & RBF decreases much faster than that of the other models. (Right) Test set accuracy on CIFAR-10 after a single pass on the training data. Growing Weighted Softmax & RBF are the most efficient learning algorithms. To place our results, we include the test accuracy of two classical (not single pass) methods: k-nearest-neighbor (kNN) and SVM (starred; details in Appendix F.10). CIFAR-10/100 Experiments. Next, we investigate the performance of the growing RCDL models on the classic CIFAR-10 and CIFAR-100 image classification tasks [42], using standard data augmentation and non-convolutional architectures. Fig. 2 reports the performance of several parametric baselines and the growing nonparametric models based on RBF and weighted softmax. Fig. 2/Left shows test accuracy on CIFAR-10/100. We observe that both the growing weighted softmax and RBF models (Sec. 3), whose learning rules are based on plain SGD, significantly outperform the baselines trained by plain SGD, as well as Adam, AdamW and RMSProp; and the GWS model slightly outperforms the RBF variant. Using more advanced optimizers, the parametric MLP eventually surpasses the SGD-based growing models: the MLP trained by the MesaNet-derived (FOOF) optimizer outperforms all methods on CIFAR-100 and performs on par with the growing variants on CIFAR-10 (see Appendix F.8 for further discussion on the DeltaNet-based one). Similarly, Muon proves to be a strong baseline on both datasets achieving top accuracy on CIFAR-10. These results highlight the importance of developing learning rules for growing models that correspond to these advanced optimizers. Fig. 2/Middle reports the evolution of training loss on CIFAR-10. Interestingly, the training loss decreases drastically faster for the two nonparametric models compared to all the other methods, while having good generalization performance as shown above; this may be a promising property for tasks that require both good generalization and strong memorization. Figure 3: (Left) The growing weighted softmax model trained on CIFAR-10 improves with depth, and displays a non-monotonic temperature behavior w.r.t. test set acc. For the optimal value of T=10T=10 and din=3072d_in=3072, 1Tdin⟨,⟩∈[−5.54,5.54] 1T d_in k, q ∈[-5.54,5.54] implying that the softmax layer is operating in the non-linear regime. (Right) This is further confirmed by linearizing post-training the “exp” nonlinearity of the layer, which results in a significant performance degradation. Figure 4: Comparison between parametric softmax (very wide MLP with softmax nonlin.) and our growing non-parametric weighted softmax model. Despite sharing the same architecture, the models behave very differently, both in terms of optimization and generalization: The parametric model displays a widening train-test gap, that is not present in its nonparametric weighted counterpart. We train the former for 2000, the latter for 50 epochs. Finally, Fig. 2/Right shows the learning efficiency of all methods when trained only on a single pass over non-augmented CIFAR-10 dataset. In this online regime, the growing methods outperform all parametric models including Muon and Mesa. We refer to Appendix G for further experimental results and Sec. G.2 for results using the gauge-invariance version of the GWS layer (described in Sec. 3.2.1). Figure 5: CIFAR-10 performance vs final capacity for various growth thresholds τ. The model retains high test acc. while exhibiting graceful performance degradation when compressing further. Analysis. Here we analyze how the GWS model’s test performance varies as a function of the number of layers and its temperature (Sec. 3.2), and verify that the model effectively leverages the nonlinear kernel. Fig. 4 shows that adding hidden layers improves the performance on CIFAR-10 (though performance saturates with 2 hidden layers), and that there is a non-monotonic dependence on temperature. We also confirm that the model is in the nonlinear regime by ablating the trained model by linearizing the exponential function (i.e., using the approximation exp(x)≈1+x (x)≈ 1+x) in the weighted softmax, post-training, which leads to a severe drop in performance. We further compare the training and generalization behaviors of the GWS model with that of the parametric softmax baseline. Figure 4 shows their training and test accuracy on CIFAR-10 as a function of the training progress. We observe that their behaviors are very different: despite being able to train a huge set of parameters at every iteration, the parametric baseline is slower to fit the training data, despite being highly prone to overfitting; we could not find a hyperparameter configuration that avoided a widening train-test performance gap. The GWS model has no such issues and reaches much higher test accuracy overall. Mixing retrieval-centric and parametric layers. It is also straightforward to mix RCDL layers and parametric layers within the same model. To illustrate this, we conduct a sanity check experiment on vision transformers (ViTs) [18] where we replace the feed-forward block in every ViT layer with an RCDL layer as a drop-in replacement, while keeping all other layers (e.g., self-attention) parametric and training them by gradient descent. Trained for a single epoch on CIFAR-10, the purely parametric ViT achieves 48.4%48.4\% while our mixed model achieves 54.4%54.4\% test set accuracy, a promising result opening up many avenues for future research. For more details, see Appendix F.11. Sub-linear growing. Finally, to mitigate the prohibitive memory and compute growth of vanilla RCDL layers, we provide preliminary results of the hybrid layer described in Sec. 3.3 which combines parametric learning with (conditional) non-parametric growth, again when training CIFAR-10 for one epoch. As shown in Figure 5, the memory growth can be made sublinear while preserving strong performance. Further experimental details can be found in Appendix F.4. 6 Discussion, Related Work, and Conclusion The RCDL concept has connections to many existing methods in machine learning, such as kernel machines [1, 12, 14, 61] (further related work can be found in Appendix H). In particular, the principled derivation of our update rules—particularly our novel rules for the weighted softmax kernel—closely parallels functional gradient boosting [48, 22]. While this perspective remains relatively underexplored in deep learning, prior work [29] derived an algorithm for constructing networks by iteratively adding differentiable modules trained to match the functional gradient of the loss; moreover, the same work presented a functional gradient-based rule for reproducing kernel Hilbert spaces, which essentially corresponds to our result of Sec. 3.1. We complement this early study with modern deep learning experiments, and with our novel (weighted) softmax attention. Furthermore, RCDL’s growing network is fundamentally a cascade of kernel machines [1, 12, 14, 61], with features in every layer trained in an end-to-end fashion through gradient descent—a core principle of deep learning. This has an interesting difference to Domingos [17]’s view on how an entire deep N trained by gradient descent can be approximated by a kernel machine. In contrast, the development of RCDL builds on the exact layer-wise equivalence/duality between a linear layer trained by gradient descent and a kernel machine [1, 33]. This might open new avenues to study optimization in NNs. RCDL may also open up a new path for continual learning. Unlike conventional weight-based learning which involves possibly irreversible compression during training, RCDL preserves past activations; we may potentially augment RCDL with some algorithms that revive old memories (e.g., within the forward pass of the net). Loosely, this property also aligns with the neurobiological findings on “silent memory” [58, 56, 25] where forgetting occurs due to retrieval failures (not because the actual memory content is lost) and some (optogenetic) procedure can revive the apparently lost memory. Limitations. RCDL has several open challenges. One of the biggest concerns is the memory and compute requirements that linearly grow with the number of training datapoints, which is generally prohibitive. While some of the key-value pruning methods from the transformer literature (e.g., sliding window or token merging [11]) may provide some insights, deriving a more principled solution for RCDL remains a future research topic. At the same time, the internet (and the corresponding data storage) itself continues to grow, and memory continues to be cheaper as hardware advances—ever growing RCDL may not be completely inconceivable in a distant future. Another notable limitation is the class of base optimization algorithms currently supported by our theory; while we covered some of the basic concepts such as weight decay and momentum to augment vanilla gradient descent, we do not know yet how to derive more advanced optimizers such as Adam [41] and Muon [39] for RCDL. Conclusion. Retrieval-centric deep learning (RCDL) aims to improve upon conventional deep learning by replacing the “linear attention” (LA) underlying the classic linear layer trained via gradient descent—a fundamental building block of deep learning—with more powerful attention mechanisms. While replacing LA with MesaNet/DeltaNet-inspired update rules yields efficient optimizers for conventional weight-based neural nets, replacing it with kernelized and softmax attention leads to unconventional nonparametric networks that learn by growing a key-value memory database instead of a fixed-size weight matrix. Using techniques from functional analysis and functional gradient methods, we derived gradient-based learning rules to train deep networks composed of the latter, and we showed empirically that nonparametric deep networks perform strongly on standard benchmarks. Our findings open a number of avenues for future work, such as investigating how to compress the growing nonparametric layers or integrate such layers into language models. References [1] M. A. Aizerman, E. M. Braverman, and L. I. Rozonoer (1964) Theoretical foundations of potential function method in pattern recognition. Automation and Remote Control 25 (6), p. 917–936. Cited by: §1, §2, §6, §6. [2] N. Aronszajn (1950) Theory of reproducing kernels. Transactions of the American mathematical society 68 (3), p. 337–404. Cited by: §3.1. [3] J. Ba, J. R. Kiros, and G. E. Hinton (2016) Layer normalization. Preprint arXiv:1607.06450. Cited by: §F.8.1. [4] D. Bahdanau, K. Cho, and Y. Bengio (2015) Neural machine translation by jointly learning to align and translate. In Int. Conf. on Learning Representations (ICLR), San Diego, CA, USA. Cited by: §1, §2. [5] A. Behrouz, Z. Li, Y. Deng, P. Zhong, M. Razaviyayn, and V. Mirrokni (2026) Memory caching: RNNs with growing memory. In Proc. Int. Conf. on Machine Learning (ICML), Seoul, South Korea. Cited by: Appendix H. [6] A. Behrouz, Z. Li, P. Kacham, M. Daliri, Y. Deng, P. Zhong, M. Razaviyayn, and V. Mirrokni (2025) Atlas: learning to optimally memorize the context at test time. Preprint arXiv:2505.23735. Cited by: §1. [7] A. Behrouz, M. Razaviyayn, P. Zhong, and V. Mirrokni (2025) Nested learning: the illusion of deep learning architectures. In Proc. Advances in Neural Information Processing Systems (NeurIPS), San Diego, CA, USA. Cited by: §1. [8] A. Behrouz, P. Zhong, and V. Mirrokni (2025) Titans: learning to memorize at test time. In Proc. Advances in Neural Information Processing Systems (NeurIPS), San Diego, CA, USA. Cited by: §1. [9] F. Benzing (2022) Gradient descent on neurons and its link to approximate second-order optimization. In Proc. Int. Conf. on Machine Learning (ICML), K. Chaudhuri, S. Jegelka, L. Song, C. Szepesvári, G. Niu, and S. Sabato (Eds.), Cited by: §4. [10] V. Berges, B. Oğuz, D. Haziza, W. Yih, L. Zettlemoyer, and G. Ghosh (2025) Memory layers at scale. In Proc. Int. Conf. on Machine Learning (ICML), Vancouver, Canada. Cited by: Appendix H. [11] D. Bolya, C. Fu, X. Dai, P. Zhang, C. Feichtenhofer, and J. Hoffman (2023) Token merging: your vit but faster. In Int. Conf. on Learning Representations (ICLR), Kigali, Rwanda. Cited by: §6. [12] B. E. Boser, I. Guyon, and V. Vapnik (1992) A training algorithm for optimal margin classifiers. In Proc. Annual ACM Conference on Computational Learning Theory (COLT), Pittsburgh, PA, USA, p. 144–152. Cited by: §6, §6. [13] J. S. Bridle (1990) Probabilistic interpretation of feedforward classification network outputs, with relationships to statistical pattern recognition. Neurocomputing: Algorithms, architectures and applications, p. 227–236. Cited by: §5. [14] C. J. Burges (1998) A tutorial on support vector machines for pattern recognition. Data mining and knowledge discovery 2 (2), p. 121–167. Cited by: §6, §6. [15] E. Butkus, K. G. Gupta, and N. Kriegeskorte (2026) Growing a neural network in breadth, depth, and time. Preprint arXiv:2605.25174. Cited by: Appendix H. [16] L. Chizat, E. Oyallon, and F. Bach (2019) On lazy training in differentiable programming. In Proc. Advances in Neural Information Processing Systems (NeurIPS), Vancouver, Canada. Cited by: Appendix H. [17] P. Domingos (2020) Every model learned by gradient descent is approximately a kernel machine. Preprint arXiv:2012.00152. Cited by: §6. [18] A. Dosovitskiy, L. Beyer, A. Kolesnikov, D. Weissenborn, X. Zhai, T. Unterthiner, M. Dehghani, M. Minderer, G. Heigold, S. Gelly, J. Uszkoreit, and N. Houlsby (2021) An image is worth 16x16 words: transformers for image recognition at scale. In Int. Conf. on Learning Representations (ICLR), Virtual only. Cited by: §5. [19] U. Evci, B. van Merrienboer, T. Unterthiner, F. Pedregosa, and M. Vladymyrov (2022) GradMax: growing neural networks using gradient information. In Int. Conf. on Learning Representations (ICLR), Virtual only. Cited by: Appendix H. [20] S. Eyuboglu, R. Ehrlich, S. Arora, N. Guha, D. Zinsley, E. Liu, W. Tennien, A. Rudra, J. Zou, A. Mirhoseini, et al. (2025) Cartridges: lightweight and general-purpose long context representations via self-study. Preprint arXiv:2506.06266. Cited by: Appendix H. [21] S. Fahlman and C. Lebiere (1989) The cascade-correlation learning architecture. In Proc. Advances in Neural Information Processing Systems (NIPS), Vol. 2. Cited by: Appendix H. [22] J. H. Friedman (2000) Greedy function approximation: a gradient boosting machine. Annals of Statistics 29, p. 1189–1232. Cited by: §6. [23] K. Fukushima (1969) Visual feature extraction by a multilayered network of analog threshold elements. IEEE Transactions on Systems Science and Cybernetics 5 (4), p. 322–333. Cited by: §5. [24] E. Gardner and B. Derrida (1989) Three unfinished works on the optimal storage capacity of networks. Journal of Physics A: Mathematical and General 22 (12). Cited by: §5. [25] S. J. Gershman, I. Fiete, and K. Irie (2025) Key-value memory in the brain. Neuron 113. Cited by: §6. [26] M. Geva, R. Schuster, J. Berant, and O. Levy (2021) Transformer feed-forward layers are key-value memories. In Proc. Conf. on Empirical Methods in Natural Language Processing (EMNLP), Punta Cana, Dominican Republic, p. 5484–5495. Cited by: §5. [27] A. Graves, G. Wayne, and I. Danihelka (2014) Neural turing machines. Preprint arXiv:1410.5401. Cited by: Appendix H. [28] A. Graves, G. Wayne, M. Reynolds, T. Harley, I. Danihelka, A. Grabska-Barwińska, S. G. Colmenarejo, E. Grefenstette, T. Ramalho, J. Agapiou, et al. (2016) Hybrid computing using a neural network with dynamic external memory. Nature 538 (7626), p. 471–476. Cited by: Appendix H. [29] A. Grubb and J. A. Bagnell (2010) Boosted backpropagation learning for training deep modular networks. In Proc. Int. Conf. on Machine Learning (ICML), J. Fürnkranz and T. Joachims (Eds.), Cited by: §6. [30] D. Hendrycks and K. Gimpel (2016) Gaussian Error Linear Units (GELUs). arXiv preprint arXiv:1606.08415. Cited by: §F.8.1. [31] D. Hendrycks (2016) Gaussian error linear units (GELUs). Preprint arXiv:1606.08415. Cited by: §5. [32] G. Hinton, N. Srivastava, and K. Swersky (2012) Neural networks for machine learning (lecture 6e). Note: Coursera, video lectures Cited by: §F.9, §5. [33] K. Irie, R. Csordás, and J. Schmidhuber (2022) The dual form of neural networks revisited: connecting test time predictions to training patterns via spotlights of attention. In Proc. Int. Conf. on Machine Learning (ICML), Baltimore, MD, USA. Cited by: §1, §1, §1, §2, §6, Abstract. [34] K. Irie, F. Faccio, and J. Schmidhuber (2022) Neural differential equations for learning to program neural nets through continuous learning rules. In Proc. Advances in Neural Information Processing Systems (NeurIPS), New Orleans, LA, USA. Cited by: §1. [35] K. Irie, I. Schlag, R. Csordás, and J. Schmidhuber (2021) Going beyond linear transformers with recurrent fast weight programmers. In Proc. Advances in Neural Information Processing Systems (NeurIPS), Virtual only. Cited by: §1. [36] K. Irie, I. Schlag, R. Csordás, and J. Schmidhuber (2022) A modern self-referential weight matrix that learns to modify itself. In Proc. Int. Conf. on Machine Learning (ICML), Baltimore, MA, USA, p. 9660–9677. Cited by: §1. [37] K. Irie (2026) The dual form of the perceptron (1964) meets deep learning and what that enables. Preprint. Cited by: §1. [38] A. Jacot, F. Gabriel, and C. Hongler (2018) Neural tangent kernel: convergence and generalization in neural networks. In Proc. Advances in Neural Information Processing Systems (NeurIPS), Montréal, Canada. Cited by: Appendix H. [39] K. Jordan, Y. Jin, V. Boza, J. You, F. Cesista, L. Newhouse, and J. Bernstein (2024) Muon: an optimizer for hidden layers in neural networks. External Links: Link Cited by: §F.9, §1, §5, §6. [40] A. Katharopoulos, A. Vyas, N. Pappas, and F. Fleuret (2020) Transformers are RNNs: fast autoregressive transformers with linear attention. In Proc. Int. Conf. on Machine Learning (ICML), Virtual only. Cited by: §1, §2. [41] D. P. Kingma and J. Ba (2015) Adam: a method for stochastic optimization. In Int. Conf. on Learning Representations (ICLR), San Diego, CA, USA. Cited by: §F.12, §F.9, §5, §6. [42] A. Krizhevsky (2009) Learning multiple layers of features from tiny images. Master’s Thesis, Computer Science Department, University of Toronto. Cited by: §F.1.1, §5. [43] J. Lee, L. Xiao, S. Schoenholz, Y. Bahri, R. Novak, J. Sohl-Dickstein, and J. Pennington (2019) Wide neural networks of any depth evolve as linear models under gradient descent. In Proc. Advances in Neural Information Processing Systems (NeurIPS), Vancouver, Canada. Cited by: Appendix H. [44] P. Lewis, E. Perez, A. Piktus, F. Petroni, V. Karpukhin, N. Goyal, H. Küttler, M. Lewis, W. Yih, T. Rocktäschel, S. Riedel, and D. Kiela (2020) Retrieval-augmented generation for knowledge-intensive NLP tasks. In Proc. Advances in Neural Information Processing Systems (NeurIPS), Virtual only. Cited by: Appendix H. [45] I. Loshchilov and F. Hutter (2019) Decoupled weight decay regularization. In Int. Conf. on Learning Representations (ICLR), New Orleans, LA, USA. Cited by: §D.3, §F.9, §5. [46] M. Luong, H. Pham, and C. D. Manning (2015) Effective approaches to attention-based neural machine translation. In Proc. Conf. on Empirical Methods in Natural Language Processing (EMNLP), Lisbon, Portugal. Cited by: §2. [47] J. Martens and R. B. Grosse (2015) Optimizing neural networks with kronecker-factored approximate curvature. In Proceedings of the 32nd International Conference on Machine Learning, ICML 2015, Lille, France, 6-11 July 2015, F. R. Bach and D. M. Blei (Eds.), JMLR Workshop and Conference Proceedings, Vol. 37, p. 2408–2417. External Links: Link Cited by: §4. [48] L. Mason, J. Baxter, P. L. Bartlett, and M. R. Frean (1999) Boosting algorithms as gradient descent. In Proc. Advances in Neural Information Processing Systems (NIPS), S. A. Solla, T. K. Leen, and K. Müller (Eds.), Cited by: §6. [49] S. Min, W. Shi, M. Lewis, X. Chen, W. Yih, H. Hajishirzi, and L. Zettlemoyer (2023) Nonparametric masked language modeling. In Findings of the Association for Computational Linguistics, Toronto, Canada, p. 2097–2118. Cited by: Appendix H. [50] V. Nair and G. E. Hinton (2010) Rectified linear units improve restricted boltzmann machines. In Proc. Int. Conf. on Machine Learning (ICML), Haifa, Israel, p. 807–814. Cited by: §5. [51] H. Peng, N. Pappas, D. Yogatama, R. Schwartz, N. A. Smith, and L. Kong (2021) Random feature attention. In Int. Conf. on Learning Representations (ICLR), Virtual only. Cited by: §1. [52] J. Platt (1991) A resource-allocating network for function interpolation. Neural Computation 3 (2), p. 213–225. Cited by: Appendix H. [53] A. Pritzel, B. Uria, S. Srinivasan, A. P. Badia, O. Vinyals, D. Hassabis, D. Wierstra, and C. Blundell (2017) Neural episodic control. In Proc. Int. Conf. on Machine Learning (ICML), Lille, France, p. 2827–2836. Cited by: Appendix H. [54] H. Ramsauer, B. Schäfl, J. Lehner, P. Seidl, M. Widrich, L. Gruber, M. Holzleitner, T. Adler, D. Kreil, M. K. Kopp, G. Klambauer, J. Brandstetter, and S. Hochreiter (2021) Hopfield networks is all you need. In Int. Conf. on Learning Representations (ICLR), Virtual only. Cited by: Appendix H. [55] H. Robbins and S. Monro (1951) A stochastic approximation method. The Annals of Mathematical Statistics 22 (3), p. 400–407. Cited by: §F.9. [56] D. S. Roy, S. Muralidhar, L. M. Smith, and S. Tonegawa (2017) Silent memory engrams as the basis for retrograde amnesia. Proc. of National Academy of Sciences (PNAS) 114 (46), p. E9972–E9979. Cited by: §6. [57] A. A. Rusu, N. C. Rabinowitz, G. Desjardins, H. Soyer, J. Kirkpatrick, K. Kavukcuoglu, R. Pascanu, and R. Hadsell (2016) Progressive neural networks. Preprint arXiv:1606.04671. Cited by: Appendix H. [58] T. J. Ryan, D. S. Roy, M. Pignatelli, A. Arons, and S. Tonegawa (2015) Engram cells retain memory under retrograde amnesia. Science 348 (6238), p. 1007–1013. Cited by: §6. [59] I. Schlag, K. Irie, and J. Schmidhuber (2021) Linear Transformers are secretly fast weight programmers. In Proc. Int. Conf. on Machine Learning (ICML), Virtual only. Cited by: §1, §1, §4. [60] J. Schmidhuber (1992) Learning to control fast-weight memories: an alternative to dynamic recurrent networks. Neural Computation 4 (1), p. 131–139. Cited by: §1. [61] B. Schölkopf and A. J. Smola (2002) Learning with kernels: support vector machines, regularization, optimization, and beyond. MIT press. Cited by: §1, §3.1, §6, §6. [62] Z. Shen, M. Zhang, H. Zhao, S. Yi, and H. Li (2018) Efficient attention: attention with linear complexities. Preprint arXiv:1812.01243. Cited by: §1, §2. [63] J. Siems, T. Carstensen, A. Zela, F. Hutter, M. Pontil, and R. Grazzi (2025) DeltaProduct: improving state-tracking in linear RNNs via householder products. In Proc. Advances in Neural Information Processing Systems (NeurIPS), San Diego, CA, USA. Cited by: §1. [64] B. Simsek, F. Ged, A. Jacot, F. Spadaro, C. Hongler, W. Gerstner, and J. Brea (2021) Geometry of the loss landscape in overparameterized neural networks: symmetries and invariances. In Proc. Int. Conf. on Machine Learning (ICML), Virtual only. Cited by: §5. [65] K. O. Stanley and R. Miikkulainen (2002) Evolving neural networks through augmenting topologies. Evolutionary computation 10 (2), p. 99–127. Cited by: Appendix H. [66] S. Sukhbaatar, A. Szlam, J. Weston, and R. Fergus (2015) End-to-end memory networks. In Proc. Advances in Neural Information Processing Systems (NIPS), Montréal, Canada, p. 2440–2448. Cited by: Appendix H. [67] Y. Sun, X. Li, K. Dalal, J. Xu, A. Vikram, G. Zhang, Y. Dubois, X. Chen, X. Wang, S. Koyejo, et al. (2025) Learning to (learn at test time): RNNs with expressive hidden states. In Proc. Int. Conf. on Machine Learning (ICML), Vancouver, Canada. Cited by: §E.1, §1. [68] A. M. Turing (1936) On computable numbers, with an application to the entscheidungsproblem. Proc. of London Mathematical Society. Cited by: Appendix H. [69] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin (2017) Attention is all you need. In Proc. Advances in Neural Information Processing Systems (NIPS), Long Beach, CA, USA, p. 5998–6008. Cited by: Appendix H, §1, §2, §3.2, §5. [70] J. von Oswald, E. Niklasson, E. Randazzo, J. Sacramento, A. Mordvintsev, A. Zhmoginov, and M. Vladymyrov (2023) Transformers learn in-context by gradient descent. In Proc. Int. Conf. on Machine Learning (ICML), Honolulu, HI, USA. Cited by: §1. [71] J. von Oswald, E. Niklasson, M. Schlegel, S. Kobayashi, N. Zucchet, N. Scherrer, N. Miller, M. Sandler, M. Vladymyrov, R. Pascanu, et al. (2023) Uncovering mesa-optimization algorithms in Transformers. Preprint arXiv:2309.05858. Cited by: §1, §4. [72] J. von Oswald, N. Scherrer, S. Kobayashi, L. Versari, S. Yang, M. Schlegel, K. Maile, Y. Schimpf, O. Sieberling, A. Meulemans, et al. (2026) MesaNet: sequence modeling by locally optimal test-time training. In Int. Conf. on Learning Representations (ICLR), Rio de Janeiro, Brazil. Cited by: §1, §1, §4, §4. [73] K. A. Wang, J. Shi, and E. B. Fox (2025) Test-time regression: a unifying framework for designing sequence models with associative memory. Preprint arXiv:2501.12352. Cited by: §1. [74] J. Weston, S. Chopra, and A. Bordes (2014) Memory networks. Preprint arXiv:1410.3916. Cited by: Appendix H. [75] B. Widrow and M. E. Hoff (1960) Adaptive switching circuits. In Proc. IRE WESCON Convention Record, Los Angeles, CA, USA, p. 96–104. Cited by: §4. [76] A. G. Wilson, Z. Hu, R. Salakhutdinov, and E. P. Xing (2016) Deep kernel learning. In Proc. Int. Conf. on Artificial Intelligence and Statistics (AISTATS), Cadiz, Spain, p. 370–378. Cited by: Appendix H. [77] S. Yang, J. Kautz, and A. Hatamizadeh (2025) Gated delta networks: improving Mamba2 with delta rule. In Int. Conf. on Learning Representations (ICLR), Vancouver, Canada. Cited by: §1, §1. [78] S. Yang, B. Wang, Y. Shen, R. Panda, and Y. Kim (2024) Gated linear attention transformers with hardware-efficient training. In Proc. Int. Conf. on Machine Learning (ICML), Vienna, Austria. Cited by: §1. [79] S. Yang, B. Wang, Y. Zhang, Y. Shen, and Y. Kim (2024) Parallelizing linear transformers with the delta rule over sequence length. In Proc. Advances in Neural Information Processing Systems (NeurIPS), Vancouver, Canada. Cited by: §E.1, §1. [80] J. Yoon, E. Yang, J. Lee, and S. J. Hwang (2018) Lifelong learning with dynamically expandable networks. In Int. Conf. on Learning Representations (ICLR), Vancouver, Canada. Cited by: Appendix H. [81] B. Zhang and R. Sennrich (2019) Root mean square layer normalization. In Proc. Advances in Neural Information Processing Systems (NeurIPS), Vancouver, Canada. Cited by: §5. Appendix A Derivation of the Growing Layer Algorithm for Positive Definite Kernels In this appendix, we provide the formal derivation of the learning algorithms for the growing kernel networks introduced in Section 3.1. We establish the correspondence between the “growing database” update rules and stochastic gradient descent (SGD) performed on an implicit, infinite-dimensional weight matrix defined in a Reproducing Kernel Hilbert Space (RKHS). A.1 The Growing Kernel Layer A growing kernel layer with input dimension dind_in and output dimension doutd_out is defined by a tuple (ϕ,)(φ,D), where: • ϕ:ℝdin×ℝdin→ℝφ:R^d_in×R^d_in is a positive definite kernel. • =(i,i)i=1MD=\( k_i, v_i)\_i=1^M is a finite database of key-value pairs, where i∈ℝdin k_i ^d_in and i∈ℝdout v_i ^d_out. The layer implements a function fϕ,:ℝdin→ℝdoutf_φ,D:R^d_in ^d_out defined as: fϕ,()=∑(,)∈ϕ(,).f_φ,D( q)= _( k, v) φ( k, q) v\,. (14) Unlike standard parametric layers with fixed dimensions, this layer “grows” by appending new key-value pairs to D. Deep growing kernel networks are constructed by composing such layers, interleaving them with element-wise non-linearities σ. For example, a two-layer network is represented as F()=σ(fϕ2,2(σ(fϕ1,1())))F( x)=σ(f_ _2,D_2(σ(f_ _1,D_1( x)))). A.2 The Implicit Weight Matrix View By the Moore-Aronszajn theorem, the positive definite kernel ϕφ is associated with a Reproducing Kernel Hilbert Space (RKHS) ℋϕH_φ and a feature map ψϕ:ℝdin→ℋϕ _φ:R^d_in _φ such that ϕ(,′)=⟨ψϕ(),ψϕ(′)⟩ℋϕφ( x, x )= _φ( x), _φ( x ) _H_φ. Since ℋϕH_φ is typically separable (e.g., for the RBF kernel), it is isomorphic to ℓ2(ℕ) _2(N). We can thus view ψϕ() _φ( x) as an infinite-dimensional vector. Substituting this inner product structure into Eq. 14, we obtain: fϕ,() f_φ,D( q) =∑(,)∈⟨ψϕ(),ψϕ()⟩ℋϕ = _( k, v) v _φ( k), _φ( q) _H_φ =(∑(,)∈ψϕ()⊤)ψϕ() = ( _( k, v) v _φ( k) ) _φ( q) =ϕ,ψϕ(), = W_φ,D _φ( q)\,, (15) where ϕ,=∑(,)∈ψϕ()⊤ W_φ,D= _( k, v) v _φ( k) is an implicit linear operator (or infinite matrix) mapping from ℋϕH_φ to ℝdoutR^d_out. This formulation reveals that a growing kernel layer is theoretically equivalent to a linear layer with a weight matrix ∈ℝdout×∞ W ^d_out×∞ acting on the feature space, where the weight matrix has a finite rank determined by the size of D. A.3 Derivation of the Learning Algorithm We now derive the update rule for the database D by minimizing a global loss function ℒL via Stochastic Gradient Descent (SGD) on the implicit weights W. Consider a deep kernel network F and a specific layer ℓ with implicit weight matrix (ℓ) W^( ), kernel ϕ(ℓ)φ^( ), and feature map ψ(ℓ)ψ^( ). Let hin(ℓ)()h_in^( )( x) denote the input to this layer and hout(ℓ)()=(ℓ)ψ(ℓ)(hin(ℓ)())h_out^( )( x)= W^( )ψ^( )(h_in^( )( x)) denote the output. Given a mini-batch of training data B=(b,b)b=1|B|B=\( x_b, y_b)\_b=1^|B|, the loss is =∑b=1|B|ℒ(b,b,(ℓ))J= _b=1^|B|L( x_b, y_b; W^( )). The gradient of the loss with respect to the implicit matrix (ℓ) W^( ) is computed via the chain rule. Let b(ℓ)e_b^( ) be the error signal (gradient of the loss with respect to the layer output) for sample b: b(ℓ):=∇hout(ℓ)ℒ(b,b)∈ℝdout(ℓ).e_b^( ):= _h_out^( )L( x_b, y_b) ^d_out^( )\,. (16) The gradient with respect to (ℓ) W^( ) is the outer product of the error signal and the layer input features: ∇(ℓ)=∑b=1|B|b(ℓ)ψ(ℓ)(hin(ℓ)(b))⊤. _ W^( )J= _b=1^|B|e_b^( )ψ^( )(h_in^( )( x_b)) \,. (17) Applying the SGD update rule t+1(ℓ)=t(ℓ)−η∇(ℓ) W_t+1^( )= W_t^( )-η _ W^( )J yields: t+1(ℓ)=t(ℓ)−η∑b=1|B|b(ℓ)ψ(ℓ)(hin(ℓ)(b))⊤. W_t+1^( )= W_t^( )-η _b=1^|B|e_b^( )ψ^( )(h_in^( )( x_b)) \,. (18) We now examine the function computed by the layer at step t+1t+1 for an arbitrary query q: ft+1(ℓ)() f_t+1^( )( q) =t+1(ℓ)ψ(ℓ)() = W_t+1^( )ψ^( )( q) =t(ℓ)ψ(ℓ)()+∑b=1|B|(−ηb(ℓ))ψ(ℓ)(hin(ℓ)(b))⊤ψ(ℓ)() = W_t^( )ψ^( )( q)+ _b=1^|B| (- _b^( ) )ψ^( )(h_in^( )( x_b)) ψ^( )( q) =ft(ℓ)()+∑b=1|B|new,bϕ(ℓ)(new,b,), =f_t^( )( q)+ _b=1^|B| v_new,bφ^( )( k_new,b, q)\,, (19) where new,b=hin(ℓ)(b) k_new,b=h_in^( )( x_b) and new,b=−ηb(ℓ) v_new,b=- _b^( ). This confirms that performing SGD on the implicit infinite-dimensional weight matrix is exactly equivalent to appending the current batch inputs (as keys) and their scaled negative gradients (as values) to the database D. A.4 Optimization Dynamics: Weight Decay and Momentum We extend the framework to include weight decay and momentum, deriving their dual forms in the growing database context. A.4.1 Weight Decay and Database Pruning Consider the L2L_2-regularized objective ℒreg()=ℒ()+γ2η‖HS2L_reg( W)=L( W)+ γ2η\| W\|_HS^2, where ∥⋅∥HS\|·\|_HS is the Hilbert-Schmidt (Frobenius) norm and γ∈(0,1)γ∈(0,1) is the decay factor. The SGD update becomes: t+1=(1−γ)t−η∇ℒ. W_t+1=(1-γ) W_t-η _ WL\,. (20) In the functional form f()=ψ()f( q)= Wψ( q), the scaling factor (1−γ)(1-γ) distributes linearly over the sum of key-value pairs. Thus, weight decay is implemented by scaling down the value vectors of all existing pairs in the database at each step: t+1=(,(1−γ))∣(,)∈t∪(new,new).D_t+1=\( k,(1-γ) v) ( k, v) _t\∪\( k_new, v_new)\\,. (21) Pruning: This update implies that the magnitude of a value vector stored at time t0t_0 decays as ‖(t)‖=(1−γ)t−t0‖(t0)‖\| v(t)\|=(1-γ)^t-t_0\| v(t_0)\|. This justifies a principled pruning strategy: once ‖(t)‖\| v(t)\| drops below a threshold ϵε, the pair’s contribution becomes negligible, and it can be removed from memory, bounding the database size. A.4.2 Momentum We adopt the standard momentum formulation where a velocity matrix t V_t accumulates gradients. Let t=−η∇ℒt G_t=-η _ WL_t be the update term at step t. The velocity update is t+1=μt+t V_t+1=μ V_t+ G_t, and the weight update is t+1=t+t+1 W_t+1= W_t+ V_t+1. Unrolling the recurrence, the weight matrix at time T is a linear combination of all past update terms τ G_τ (created at time τ<Tτ<T): T=∑τ=1Tc(T,τ)τ,with c(T,τ)=∑k=0T−τμk=1−μT−τ+11−μ. W_T= _τ=1^Tc(T,τ) G_τ\,, c(T,τ)= _k=0^T-τμ^k= 1-μ^T-τ+11-μ\,. (22) Since each τ G_τ corresponds to the batch of key-value pairs created at time τ, we do not need to store explicit velocity vectors. Instead, we store the creation timestamp tcreatet_create with each pair. During the forward pass at time T, we dynamically compute the coefficient c(T,tcreate)c(T,t_create): fT()=∑(,init,tcreate)∈(1−μT−tcreate+11−μ)ϕ(,)init.f_T( q)= _( k, v_init,t_create) ( 1-μ^T-t_create+11-μ )φ( k, q) v_init\,. (23) This allows exact momentum computation with no additional memory overhead per pair. Combining both techniques results in decaying the stored values for regularization while scaling the retrieval output for momentum. Appendix B Further Details and Remarks for the Softmax Case B.1 Remark: the functional space of weighted softmax functions In order to make the softmax function in Eq. 3 well-defined, we require that ∫ℝdin×ℝdoutw(,)e⟨,⟩+⟨,′⟩ _R^d_in×R^d_outw( k, v)e k, q + v, v dd<∞, d k\,d v<∞\,, ∀∈ℝdin,∀′∈ℝdout. ∀ q ^d_in\,,∀ v ^d_out\,. In other words, if we interpret w as the density of a measure μw _w on ℝdin×ℝdoutR^d_in×R^d_out, then we require that the moment-generating function of μw _w be well-defined and finite everywhere. Remark B.1. In the case of a finite database of weighted key-value pairs, if we have several tuples having the same key, then we can combine them into one tuple. For example, we can combine (w1,,1)(w_1, k, v_1) and (w2,,2)(w_2, k, v_2) and replace them with (w1+w2,,w11+w22w1+w2) (w_1+w_2, k, w_1 v_1+w_2 v_2w_1+w_2 ), and this will induce the same softmax function. Analogously, if we define the “marginalized” key-weight function w~:ℝdin→ℝ+ w:R^d_in ^+ as w~()=∫ℝdoutw(,), w( k)= _R^d_outw( k, v)d v\,, and the marginalized value function v~:ℝdin→ℝdout v:R^d_in ^d_out which associates a unique value for every key as v~():=1w~()∫ℝdoutw(,), v( k):= 1 w( k) _R^d_outw( k, v) vd v\,, then we can write sws_w in terms of w~ w and v~ v as follows: sw()=sw~,v~():=∫ℝdinw~()e⟨,⟩v~()∫ℝdinw~()e⟨,⟩.s_w( q)=s_ w, v( q):= _R^d_in w( k)e k, q v( k)\,d k _R^d_in w( k)e k, q \,d k\,. (24) B.2 Continuous softmax function induced by the Gaussian mixture density & Discrete weighted softmax function Here we provide a full derivation demonstrating that the continuous softmax function induced by the Gaussian mixture density w_D is exactly equivalent to the discrete weighted softmax function s_D defined by the database D. sw() s_w_D( q) =∫ℝdin×ℝdoutw(,)e⟨,⟩∫ℝdin×ℝdoutw(,)e⟨,⟩ = _R^d_in×R^d_outw_D( k, v)e k, q v\,d k\,d v _R^d_in×R^d_outw_D( k, v)e k, q \,d k\,d v =∑i=1nwi(,)∼((i,i),Σ)[e⟨,⟩]∑i=1nwi(,)∼((i,i),Σ)[e⟨,⟩] = _i=1^nw_iE_( k, v) (( k_i, v_i), )[e k, q v] _i=1^nw_iE_( k, v) (( k_i, v_i), )[e k, q ] =∑i=1nwi∼(i,σk2)[e⟨,⟩]∼(i,σv2)[]∑i=1nwi∼(i,σk2)[e⟨,⟩] = _i=1^nw_iE_ k ( k_i, _k^2 I)[e k, q ]E_ v ( v_i, _v^2 I)[ v] _i=1^nw_iE_ k ( k_i, _k^2 I)[e k, q ] =∑i=1nwie⟨i,⟩+12σk2‖2i∑i=1nwie⟨i,⟩+12σk2‖2 = _i=1^nw_ie k_i, q + 12 _k^2\| q\|^2 v_i _i=1^nw_ie k_i, q + 12 _k^2\| q\|^2 =e12σk2‖2∑i=1nwie⟨i,⟩ie12σk2‖2∑i=1nwie⟨i,⟩ = e 12 _k^2\| q\|^2 _i=1^nw_ie k_i, q v_ie 12 _k^2\| q\|^2 _i=1^nw_ie k_i, q =∑i=1nwie⟨i,⟩i∑i=1nwie⟨i,⟩=s(). = _i=1^nw_ie k_i, q v_i _i=1^nw_ie k_i, q =s_D( q)\,. B.3 Further discussion on the standard (unweighted) softmax case In Sec. 3.2 of the main text, we introduced the weighted softmax by noting that deriving a principled learning rule for the standard/unweighted softmax is non-trivial due to the discontinuous nature. Here we illustrate this difficulty by presenting an example of a naive derivation of a learning rule for the standard softmax case and showing where it breaks down. Model description. Consider a growing softmax layer that, at every training step t, maintains a key-value memory state consisting of two matrices t−1=[1,…,t−1]∈ℝdin×(t−1) K_t-1=[ k_1,…, k_t-1] ^d_in×(t-1) and t−1=[1,…,t−1]∈ℝdout×(t−1) V_t-1=[ v_1,…, v_t-1] ^d_out×(t-1) (both to be defined later), receives a layer input t∈ℝdin x_t ^d_in and produces a layer output t∈ℝdout y_t ^d_out through the softmax attention function: t q_t =t = x_t (25) t y_t =t−1softmax(t−1⊤t)=∑τ=1t−1αt−1,ττ = V_t-1softmax( K_t-1 q_t)= _τ=1^t-1 _t-1,τ v_τ (26) where αt−1,τ∈ℝ _t-1,τ denotes the attention score for the key τ k_τ using the query t q_t and the system’s key-value memory state from step t−1t-1 (i.e., before updating the system at step t). That is, t x_t is used as the query t q_t in the attention operation to retrieve from the key-value memory consisting of matrices t−1 K_t-1 and t−1 V_t-1. Then, the system’s memory state is updated as: t K_t =[t−1,t]∈ℝdin×t =[ K_t-1, k_t] ^d_in× t (27) t V_t =[t−1,t]∈ℝdout×t =[ V_t-1, v_t] ^d_out× t (28) where [][] denotes concatenation along the time dimension. By following the linear case (Sec. 2), let us set the key to be the layer input, i.e., t=t k_t= x_t. Now the question is how should we set t∈ℝdout v_t ^d_out so that this memory update corresponds to a principled learning algorithm for this system in a way that it optimizes a given objective function ℒL? Direct Derivation. For that, we first set the new value vector t∈ℝdout v_t ^d_out as an arbitrary free parameter vector, and consider the virtually updated system with the newly added key-value pair (t,)( x_t, v). After this update, the output of the system on the input t=t q_t= x_t becomes: tnew y_t^new =tsoftmax(t⊤t)=∑τ=1t−1αt,τnewτ+αt,tnew = V_tsoftmax( K_t q_t)= _τ=1^t-1 _t,τ^new v_τ+ _t,t^new v (29) that is, compared to the old system of Eq. 26, one extra term (the last term in Eq. 29) is introduced in the sum as we now have one more value vector, with newly computed attention scores αt,τnew _t,τ^new normalized over the new set of keys, including the attention score αt,tnew _t,t^new for the newly added key t x_t. Our goal is to improve the loss on input t q_t by adding the new value token t v_t to the system. For that, we can apply a step of gradient descent on loss ℒL w.r.t. t v_t with a learning rate ηt∈ℝ>0 _t _>0, which yields: t=old−ηt∂ℒt(tnew)∂ v_t= v^old- _t _t( y_t^new)∂ v (30) where two more variables, old v^old and the gradient term, are still to be determined; as we see below, this raises difficulties we cannot trivially overcome. First, what should old v^old be? Essentially, old v^old is the value for newly added t v_t that preserves the old system (representing the initial system before applying a step of gradient descent). For this, we can express the output of the new system tnew y_t^new (Eq. 29) as a function of that of the old one t y_t (Eq. 26) and search for the value of t v_t that makes them equal. The apparently tricky difference between them is the modified attention scores: αt,τnew _t,τ^new (Eq. 29) vs. αt−1,τ _t-1,τ (Eq. 26) for 1≤τ≤t−11≤τ≤ t-1. However, thankfully, we can notice a simple relation that connects the old attention scores before the update (αt−1,τ _t-1,τ) and the new ones (αt,τnew _t,τ^new) after the update. By denoting ZoldZ^old and ZnewZ^new the denominators (normalizers) of the attention scores when the input is t q_t, before and after the update, respectively, i.e., Zold=∑τ=1t−1eτ⊤t;Znew=∑τ=1teτ⊤t=Zold+et⊤t, Z^old= _τ=1^t-1e x_τ q_t\,;\,Z^new= _τ=1^te x_τ q_t=Z^old+e q_t q_t\,, (31) we can express the new attention scores as a function of the old one: αt,τnew=eτ⊤tZnew=ZoldZnew×eτ⊤tZold=ZoldZnewαt−1,τ, _t,τ^new= e x_τ q_tZ^new= Z^oldZ^new× e x_τ q_tZ^old= Z^oldZ^new _t-1,τ\,, (32) and additionally, the normalizer ratio can also be simplified, since Znew=Zold+et⊤tZ^new=Z^old+e q_t q_t, as follows: ZoldZnew Z^oldZ^new =Zold+et⊤t−et⊤tZnew=Znew−et⊤tZnew=1−et⊤tZnew=1−αt,tnew. = Z^old+e q_t q_t-e q_t q_tZ^new= Z^new-e q_t q_tZ^new=1- e q_t q_tZ^new=1- _t,t^new\,. (33) By inserting this ratio to Eq. 32, it yields: αt,τnew=(1−αt,tnew)αt−1,τ, _t,τ^new=(1- _t,t^new) _t-1,τ\,, (34) that is, the new attention scores for the old keys shrink by factor (1−αt,tnew)(1- _t,t^new). This allows us to express the new output tnew y_t^new as a function of the old one t y_t: tnew y_t^new =∑τ=1t−1αt,τnewτ+αt,tnew = _τ=1^t-1 _t,τ^new v_τ+ _t,t^new v (29) =∑τ=1t−1(1−αt,tnew)αt−1,ττ+αt,tnew=(1−αt,tnew)(∑τ=1t−1αt−1,ττ)+αt,tnew = _τ=1^t-1(1- _t,t^new) _t-1,τ v_τ+ _t,t^new v=(1- _t,t^new) ( _τ=1^t-1 _t-1,τ v_τ )+ _t,t^new v (35) =(1−αt,tnew)t+αt,tnew =(1- _t,t^new) y_t+ _t,t^new v (36) =t+αt,tnew(−t). = y_t+ _t,t^new( v- y_t)\,. (37) Here we can at least find t=old v_t= v^old such that the updated system reproduces the old system’s output on input t q_t, that is, tnew=t y_t^new= y_t, which, from Eq. 37 holds when t=old=t v_t= v^old= y_t, canceling the last term. However, this only preserves the mapping for the current query t q_t; as the new key t k_t alters the softmax denominator globally, the output for any other query will be modified, failing to achieve the global function preservation. Second, while the computation of the gradient term seems straightforward: ∂ℒt(tnew)∂=∂ℒt(tnew)∂tnew×∂tnew∂=(∇tnewℒt)αt,tnew, _t( y_t^new)∂ v= _t( y_t^new)∂ y_t^new× ∂ y_t^new∂ v=( _ y_t^newL_t) _t,t^new\,, (38) by differentiating Eq. 29 w.r.t. v, this requires access to the post-update gradient of the loss w.r.t. the new output (∇tnewℒt)( _ y_t^newL_t). While this creates a circular dependency in the general case, by setting t=t v_t= y_t as prescribed by the derivation above, ∇tnewℒt _ y_t^newL_t can be essentially replaced by ∇tℒt _ y_tL_t. Putting these together, we obtain the following update rule: t=old−ηt∂ℒt(tnew)∂=t−ηtαt,tnew∇tnewℒt, v_t= v^old- _t _t( y_t^new)∂ v= y_t- _t _t,t^new _ y_t^newL_t\,, (39) where αt,tnew=et⊤tZnew∈ℝ _t,t^new= e x_t x_tZ^new can be computed simply as we have Znew=Zold+et⊤tZ^new=Z^old+e x_t x_t for the normalizer, and ZoldZ^old is computed to compute t y_t. However, this is not a satisfactory derivation as the choice of old v^old is limited to the local function preservation on one data point. This derivation illustrates that a naive approach that consists in parameterizing the value vector as a free parameter and applying gradient descent on it does not work out of the box. Appendix C Functional Gradient Descent in the Weighted Softmax Space In this appendix, we rigorously derive the functional gradient of the loss with respect to the weight function w∈ℋσk,σvw _ _k, _v and show that the gradient descent update corresponds to adding a new key-value tuple to the database. C.1 Setup and Definitions Let ℋ=ℋσk,σvH=H_ _k, _v be the Hilbert space of weight functions w:ℝdin×ℝdout→ℝw:R^d_in×R^d_out equipped with the inner product ⟨f,g⟩ℋ=∫f(,)g(,)μ(,) f,g _H= f( k, v)g( k, v)dμ( k, v) and norm ‖w‖ℋ\|w\|_H, where μ=(,Σ)μ=N(0, ) is the centered Gaussian measure with block covariance Σ=diag(Σk,Σv)=diag(σk2,σv2) =diag( _k, _v)=diag( _k^2 I, _v^2 I). Recall the definition of the softmax functional sw()=∫ℝdin×ℝdoutw(,)e⟨,⟩μ(,)∫ℝdin×ℝdoutw(,)e⟨,⟩μ(,)=N(w)D(w),s_w( q)= _R^d_in×R^d_outw( k, v)e k, q v\,dμ( k, v) _R^d_in×R^d_outw( k, v)e k, q \,dμ( k, v)= N(w)D(w)\,, where we define the denominator functional D:ℋ→ℝD:H and the numerator operator N:ℋ→ℝdoutN:H ^d_out as follows: D(w) D(w) =∫ℝdin×ℝdoutw(,)e⟨,⟩μ(,), = _R^d_in×R^d_outw( k, v)e k, q \,dμ( k, v)\,, (40) N(w) N(w) =∫ℝdin×ℝdoutw(,)e⟨,⟩μ(,). = _R^d_in×R^d_outw( k, v)e k, q v\,dμ( k, v)\,. (41) Let Nj(w)N_j(w) denote the j-th component of the vector N(w)N(w). It is a scalar linear functional: Nj(w)=∫w(,)e⟨,⟩vjμ(,).N_j(w)= w( k, v)e k, q v_j\,dμ( k, v)\,. By defining the kernels KD(,)=e⟨,⟩,KN,j(,)=e⟨,⟩vj,K_D( k, v)=e k, q , K_N,j( k, v)=e k, q v_j\,, we can rewrite D and NjN_j as D(w) D(w) =∫ℝdin×ℝdoutw(,)KD(,)μ(,), = _R^d_in×R^d_outw( k, v)K_D( k, v)\,dμ( k, v)\,, (42) Nj(w) N_j(w) =∫ℝdin×ℝdoutw(,)KN,j(,)μ(,). = _R^d_in×R^d_outw( k, v)K_N,j( k, v)\,dμ( k, v)\,. (43) By the Riesz Representation Theorem, if the functions representing these kernels belong to ℋH, the functionals are bounded. We verify membership in ℋH by checking the condition μ[f(,)2e⟨,⟩+⟨,⟩]<∞E_μ[f( k, v)^2e k,u + v,z ]<∞ for all ,u,z. For KDK_D, we have ‖KD‖ℋ2=(,)∼μ[(e⟨,⟩)2e⟨,⟩+⟨,⟩]=[e⟨,2+⟩][e⟨,⟩].\|K_D\|_H^2=E_( k, v) μ [(e k, q )^2e k,u + v,z ]=E_ k[e k,2 q+u ]E_ v[e v,z ]\,. These are moment generating functions of Gaussian variables, which are finite everywhere. For KN,jK_N,j, we have ‖KN,j‖ℋ2=(,)∼μ[(e⟨,⟩vj)2e⟨,⟩+⟨,⟩]=[e⟨,2+⟩][vj2e⟨,⟩].\|K_N,j\|_H^2=E_( k, v) μ [(e k, q v_j)^2e k,u + v,z ]=E_ k[e k,2 q+u ]E_ v[v_j^2e v,z ]\,. The second term involves the second moment of a Gaussian under an exponential tilt, which is also finite. Thus, D(w)D(w) and N(w)N(w) are bounded linear operators. Specifically, there exist constants CN,CD<∞C_N,C_D<∞ such that for any δ∈ℋδ : ‖N(δ)‖2≤CN‖δ‖ℋ,|D(δ)|≤CD‖δ‖ℋ\|N(δ)\|_2≤ C_N\|δ\|_H, |D(δ)|≤ C_D\|δ\|_H (44) C.2 Fréchet Differentiability of sws_w We claim that the Fréchet derivative of sws_w at w, denoted dswds_w, is the linear operator acting on a perturbation δ∈ℋδ as: dsw(δ)=N(δ)D(w)−N(w)D(δ)D(w)2ds_w(δ)= N(δ)D(w)-N(w)D(δ)D(w)^2 (45) To prove this, we must show that the remainder term R(δ)=sw+δ()−sw()−dsw(δ)R(δ)=s_w+δ( q)-s_w( q)-ds_w(δ) satisfies ‖R(δ)‖=o(‖δ‖ℋ)\|R(δ)\|=o(\|δ\|_H) as ‖δ‖ℋ→0\|δ\|_H→ 0. Proof. Let us expand the difference sw+δ−sws_w+δ-s_w: sw+δ−sw s_w+δ-s_w =N(w)+N(δ)D(w)+D(δ)−N(w)D(w) = N(w)+N(δ)D(w)+D(δ)- N(w)D(w) (46) =D(w)(N(w)+N(δ))−N(w)(D(w)+D(δ))D(w)(D(w)+D(δ)) = D(w)(N(w)+N(δ))-N(w)(D(w)+D(δ))D(w)(D(w)+D(δ)) (47) =D(w)N(δ)−N(w)D(δ)D(w)2+D(w)D(δ) = D(w)N(δ)-N(w)D(δ)D(w)^2+D(w)D(δ) (48) Subtracting the candidate derivative dsw(δ)ds_w(δ): R(δ) R(δ) =D(w)N(δ)−N(w)D(δ)D(w)(D(w)+D(δ))−D(w)N(δ)−N(w)D(δ)D(w)2 = D(w)N(δ)-N(w)D(δ)D(w)(D(w)+D(δ))- D(w)N(δ)-N(w)D(δ)D(w)^2 (49) =[D(w)N(δ)−N(w)D(δ)](1D(w)(D(w)+D(δ))−1D(w)2) =[D(w)N(δ)-N(w)D(δ)] ( 1D(w)(D(w)+D(δ))- 1D(w)^2 ) (50) =[D(w)N(δ)−N(w)D(δ)](D(w)−(D(w)+D(δ))D(w)2(D(w)+D(δ))) =[D(w)N(δ)-N(w)D(δ)] ( D(w)-(D(w)+D(δ))D(w)^2(D(w)+D(δ)) ) (51) =−[D(w)N(δ)−N(w)D(δ)]D(δ)D(w)2(D(w)+D(δ)) =-[D(w)N(δ)-N(w)D(δ)] D(δ)D(w)^2(D(w)+D(δ)) (52) Now we bound the norm of R(δ)R(δ). Let w be fixed such that D(w)≠0D(w)≠ 0. For sufficiently small δ (specifically ‖δ‖ℋ<|D(w)|/(2CD)\|δ\|_H<|D(w)|/(2C_D)), we have |D(w)+D(δ)|≥|D(w)|−|D(δ)|≥|D(w)|−CD‖δ‖ℋ≥|D(w)|/2.|D(w)+D(δ)|≥|D(w)|-|D(δ)|≥|D(w)|-C_D\|δ\|_H≥|D(w)|/2\,. Furthermore, ‖D(w)N(δ)−N(w)D(δ)‖ \|D(w)N(δ)-N(w)D(δ)\| ≤|D(w)||N(δ)|+‖N(w)‖|D(δ)| ≤|D(w)|\|N(δ)\|+\|N(w)\||D(δ)| (53) ≤(|D(w)|CN+‖N(w)‖CD)‖δ‖ℋ ≤(|D(w)|C_N+\|N(w)\|C_D)\|δ\|_H (54) =C1‖δ‖ℋ. =C_1\|δ\|_H\,. (55) Substituting these back into the expression for R(δ)R(δ): ‖R(δ)‖ \|R(δ)\| ≤C1‖δ‖ℋ⋅|D(δ)||D(w)|2|D(w)+D(δ)| ≤ C_1\|δ\|_H·|D(δ)||D(w)|^2|D(w)+D(δ)| (56) ≤C1‖δ‖ℋ⋅CD‖δ‖ℋ|D(w)|2(|D(w)|/2) ≤ C_1\|δ\|_H· C_D\|δ\|_H|D(w)|^2(|D(w)|/2) (57) =2C1CD|D(w)|3‖δ‖ℋ2 = 2C_1C_D|D(w)|^3\|δ\|_H^2 (58) =O(‖δ‖ℋ2). =O(\|δ\|_H^2)\,. (59) Since ‖R(δ)‖≤C‖δ‖ℋ2\|R(δ)\|≤ C\|δ\|_H^2, it implies lim‖δ‖→0‖R(δ)‖δ‖=0 _\|δ\|→ 0 \|R(δ)\|\|δ\|=0. Thus, sws_w is Fréchet differentiable. ∎ C.3 Derivation of the Functional Gradient Let ℒ(^,)L( y, y) be a loss function that is differentiable with respect to the prediction y. By the chain rule for Fréchet derivatives, the composition J(w)=ℒ(sw(),)J(w)=L(s_w( q), y) is differentiable with respect to w, and its differential is: dJ(δ)=⟨∇^ℒ,dsw(δ)⟩dJ(δ)= _ yL,ds_w(δ) (60) Substituting dsw(δ)ds_w(δ): dJ(δ) dJ(δ) =⟨∇^ℒ,N(δ)D(w)−sw()D(δ)D(w)⟩ = _ yL, N(δ)D(w)-s_w( q) D(δ)D(w) (61) =1D(w)⟨∇^ℒ,N(δ)⟩−⟨∇^ℒ,sw()⟩D(w)D(δ). = 1D(w) _ yL,N(δ) - _ yL,s_w( q) D(w)D(δ)\,. (62) Using the integral definitions of N(δ)N(δ) and D(δ)D(δ) and the linearity of the inner product, we can rewrite this as an inner product in ℋH: dJ(δ) dJ(δ) =∫δ(,)[1D(w)⟨∇^ℒ,⟩e⟨,⟩−⟨∇^ℒ,sw()⟩D(w)e⟨,⟩]μ = δ( k, v) [ 1D(w) _ yL, v e k, q - _ yL,s_w( q) D(w)e k, q ]dμ (63) =∫δ(,)e⟨,⟩D(w)⟨∇^ℒ,−sw()⟩μ = δ( k, v) e k, q D(w) _ yL, v-s_w( q) dμ (64) =⟨δ,∇wJ⟩ℋ. = δ, _wJ _H\,. (65) By the Riesz representation theorem, the function ∇wJ _wJ inside the inner product is the unique functional gradient: ∇wJ(,)=e⟨,⟩D(w)⟨∇^ℒ,−sw()⟩. _wJ( k, v)= e k, q D(w) _ yL, v-s_w( q) \,. (66) C.4 Constructive Update Analysis We now analyze the impact of the functional update wt+1←wt−η∇wJw_t+1← w_t-η _wJ on the softmax functional swt+1()s_w_t+1( q) and demonstrate that it is equivalent to adding a specific key-value tuple to the database. The impact of this update is determined by the changes to the numerator and denominator functionals. For the denominator, the update adds a term: ΔD() D( q) =−η∫∇wJ(,)e⟨,⟩dμ(,) =-η _wJ( k, v)e k, q dμ( k, v) =−ηD(wt)∫e⟨,t+⟩⟨t,−t⟩dμ(,), =- ηD(w_t) e k, q_t+ q g_t, v- s_t dμ( k, v)\,, where t=∇^ℒ g_t= _ yL and t=swt(t) s_t=s_w_t( q_t). Since the measure μ is a centered Gaussian, μ[]=E_μ[ v]=0. Consequently, the term involving v integrates to zero, leaving only the scalar part: ΔD()=η⟨t,t⟩D(wt)μk[e⟨,t+⟩]. D( q)= η g_t, s_t D(w_t)E_ _k[e k, q_t+ q ]\,. The expectation is the moment generating function of (,Σk)N(0, _k) evaluated at t+ q_t+ q, which equals e12‖Σk1/2(t+)‖2e 12\| _k^1/2( q_t+ q)\|^2. Factoring this, we get a term proportional to e⟨Σkt,⟩e _k q_t, q . This corresponds precisely to the contribution of a new key new=Σkt k_new= _k q_t (or σk2t _k^2 q_t in the isotropic case). Similarly, for the numerator, the update adds a term involving v: ΔN() N( q) =−ηD(wt)∫e⟨,t+⟩⟨t,−t⟩dμ(,). =- ηD(w_t) e k, q_t+ q g_t, v- s_t vdμ( k, v)\,. The term ⟨t,−t⟩ g_t,- s_t v integrates to zero (zero mean). The remaining term involves the second moment μv[⊤]=ΣvE_ _v[ v v ]= _v. Thus: ΔN()=−ηD(wt)Σvtμk[e⟨,t+⟩]. N( q)=- ηD(w_t) _v g_tE_ _k[e k, q_t+ q ]\,. Comparing ΔN N and ΔD D, we see they share the same key kernel e⟨new,⟩e k_new, q . The equivalent added tuple (wnew,new,new)(w_new, k_new, v_new) is thus: new=Σkt,wnew∝η⟨t,t⟩,new∝−Σvt⟨t,t⟩. k_new= _k q_t, w_new η s_t, g_t , v_new - _v g_t s_t, g_t \,. This confirms that functional gradient descent in ℋσk,σvH_ _k, _v is constructively equivalent to growing the database with specific, tractable update rules. In order to avoid dividing by ⟨t,t⟩ s_t, g_t which could be small and might introduce numerical instabilities, we can adopt a representation where we store the weighted value vector ~=w v=w v instead of the value vector v, as in Eq. 8. In this case, the added tuple becomes new=Σkt,wnew∝η⟨t,t⟩,~new∝−ηΣvt. k_new= _k q_t, w_new η s_t, g_t , v_new -η _v g_t\,. If Σv=σv2 _v= _v^2 I and β=ησv2β=η _v^2, the weighted value vector becomes ~new∝−βt. v_new -β g_t\,. C.5 Negative Weights and Numerical Stability A critical inspection of the update rule for wneww_new reveals that its sign is determined by the inner product δt=⟨t,t⟩ _t= s_t, g_t . Since the learning rate η and normalization D(wt)D(w_t) are positive (assuming initialization with positive mass), wneww_new can be negative if the gradient of the loss forms an obtuse angle with the current prediction vector. This introduces the concept of “negative attention” or signed measures into the softmax mechanism. Negative weights act to reduce the denominator (partition function), effectively increasing the magnitude of the output vector s()s( q) in the direction of the value vectors. This occurs when the current prediction t s_t is negatively correlated with the gradient t g_t, indicating that gradient descent is attempting to scale up the prediction to minimize the loss. While considering negative weights might increase the expressive power of our model, they might introduce numerical instabilities as the denominator might become close to 0 if the negative and positive terms cancel. To ensure robustness, we propose enforcing non-negativity by clipping the new weight: wnew←max(0,wnew).w_new← (0,w_new)\,. (67) This clipping strategy can be theoretically justified as performing a descent step in a larger, decoupled function space ℋσk,σv2H^2_ _k, _v, where the numerator and denominator weight functions are optimized independently. Refer to Sec. C.6 for further details. In Sec. C.7, we propose another method for solving the negative weights issue. We leave the investigation of deriving an update rule that benefits from the potential expressive power of the negative weights while ensuring numerical stability for future research. C.6 Functional Gradient Descent in the Decoupled Weight Space In this appendix, we consider the functional gradient descent on the decoupled softmax functional defined in Eq. 68 below. We define the product Hilbert space ℋ2=ℋσk,σv⊕ℋσk,σvH^2=H_ _k, _v _ _k, _v containing pairs of weight functions (wN,wD)(w_N,w_D). The decoupled softmax functional is: swN,wD()=N(wN)D(wD)=∫wN(,)e⟨,⟩μ∫wD(,)e⟨,⟩μ.s_w_N,w_D( q)= N(w_N)D(w_D)= w_N( k, v)e k, q v\,dμ w_D( k, v)e k, q \,dμ\,. (68) We compute the Fréchet derivatives with respect to wNw_N and wDw_D independently. For wNw_N, the functional is linear in the numerator: ∇wNs=1D(wD)∇wNN. _w_Ns= 1D(w_D) _w_NN\,. The gradient of the loss J(wN,wD)J(w_N,w_D) with respect to wNw_N is: ∇wNJ _w_NJ =1D(wD)e⟨,t⟩⟨∇^ℒ,⟩. = 1D(w_D)e k, q_t _ yL, v \,. For wDw_D, the functional involves the inverse of the denominator: ∇wDs=−N(wN)D(wD)2∇wDD=−swN,wDD(wD)∇wDD. _w_Ds=- N(w_N)D(w_D)^2 _w_DD=- s_w_N,w_DD(w_D) _w_DD\,. The gradient of the loss with respect to wDw_D is: ∇wDJ _w_DJ =−1D(wD)e⟨,t⟩⟨∇^ℒ,swN,wD(t)⟩. =- 1D(w_D)e k, q_t _ yL,s_w_N,w_D( q_t) \,. We now map these updates to the growing database formulation. Updating wNw_N adds a term ΔwN∝−1D⟨t,⟩e⟨,t⟩ w_N - 1D g_t, v e k, q_t . In the “weighted value vector” representation (where we store ~ v instead of v), adding a tuple corresponds to adding a basis function to wNw_N. Specifically, if we add (~new,new)( v_new, k_new), we modify the numerator density by a term proportional to e⟨,Σk−1new⟩⟨,Σv−1~new⟩e k, _k^-1 k_new v, _v^-1 v_new . Matching terms, we find new=Σkt k_new= _k q_t and ~new∝−ηΣvt v_new -η _v g_t. This matches the numerator update derived in the coupled case. Updating wDw_D adds a term ΔwD∝1D⟨t,t⟩e⟨,t⟩ w_D 1D g_t, s_t e k, q_t . This corresponds to adding a basis function to the denominator weight density wDw_D. In the database, this means adding a tuple (wnew,new)(w_new, k_new). Matching terms, we find new=Σkt k_new= _k q_t and wnew∝η⟨t,t⟩w_new η s_t, g_t . This matches the denominator update derived in the coupled case. Therefore, the coupled update rule derived in Appendix A is a special case of the decoupled update where we enforce wN=wD=w_N=w_D=w initially but allow them to evolve. Crucially, since the optimization of wDw_D is independent of wNw_N in this formulation, applying a constraint wD≥0w_D≥ 0 (via clipping wnew=max(0,wnew)w_new= (0,w_new)) corresponds to a valid Descent step in ℋσk,σv⊕ℋσk,σvH_ _k, _v _ _k, _v which would be positively correlated with the gradient direction. This ensures numerical stability while maintaining the numerator’s ability to reduce error. C.7 Gauge-Invariance: Another Method for Solving the Negative Weights Issue As established in Sec. C.5, the constructive functional gradient update for the scalar weight of a newly added weighted key-value tuple is given by: wnew∝η⟨t,t⟩,w_new η s_t, g_t \,, (69) where t=swt(t) s_t=s_w_t( q_t) is the current prediction and t=∇^ℒ g_t= _ yL is the error signal. When the gradient forms an obtuse angle with the current prediction vector (i.e., ⟨t,t⟩<0 s_t, g_t <0), this rule mandates appending a negative weight (wnew<0w_new<0). While the non-negativity constraint can be enforced heuristically via clipping wnew←max(0,wnew)w_new← (0,w_new) (theoretically grounded in the decoupled weight space of Sec. C.6), clipping discards the denominator update entirely whenever ⟨t,t⟩<0 s_t, g_t <0. In this appendix, we develop an alternative, geometrically principled resolution to the negative weight obstruction: an augmented representation with gauge invariance. By introducing a global bias degree of freedom, we show that the underlying function space possesses an exact translation symmetry in value space. By dynamically selecting an appropriate coordinate frame (a “virtual gauge”) at each training step, we guarantee that the required weight update is strictly positive (wnew>0w_new>0) without clipping, while preserving the locality of the database update. C.7.1 Augmented Functional and Gauge Symmetry We introduce an auxiliary global bias degree of freedom ∈ℝdout b ^d_out into both the continuous and discrete weighted softmax functionals. Definition C.1 (Augmented Softmax Functionals). Let ∈ℝdout b ^d_out be a bias vector. 1. For a functional weight density w∈ℋσk,σvw _ _k, _v, the continuous augmented softmax output sw,:ℝdin→ℝdouts_w, b:R^d_in ^d_out is defined as: sw,():=sw()−=∫ℝdin×ℝdoutw(,)e⟨,⟩μ(,)∫ℝdin×ℝdoutw(,)e⟨,⟩μ(,)−.s_w, b( q):=s_w( q)- b= _R^d_in×R^d_outw( k, v)e k, q v\,dμ( k, v) _R^d_in×R^d_outw( k, v)e k, q \,dμ( k, v)- b\,. (70) 2. For a discrete database =(wi,i,~i)i=1ND=\(w_i, k_i, v_i)\_i=1^N of weighted key-value tuples with ~i=wii v_i=w_i v_i, the discrete augmented retrieval functional is defined as: s,():=s()−=∑i=1Ne⟨i,⟩~i∑i=1Nwie⟨i,⟩−.s_D, b( q):=s_D( q)- b= _i=1^Ne k_i, q v_i _i=1^Nw_ie k_i, q - b\,. (71) The standard model corresponds to a bias that is set to zero: s()=s,()s_D( q)=s_D,0( q). We now define a transformation of the representation under a spatial translation in value space. Definition C.2 (Shifted Density and Database). Let ∈ℝdout ξ ^d_out be an arbitrary shift vector. • We define the shifted weight density w_ ξ by translating the value coordinate: w(,):=w(,−)⋅exp(2⟨,⟩−‖22σv2).w_ ξ( k, v):=w( k, v- ξ)· ( 2 ξ, v -\| ξ\|^22 _v^2 )\,. (72) • In the discrete representation, we shift the value vectors as i↦i+ v_i v_i+ ξ which implies that the stored weighted values shift as ~i↦wi(i+)=~i+wi v_i w_i( v_i+ ξ)= v_i+w_i ξ. Hence, we define the shifted database D_ ξ as: :=(wi,i,~i+wi)i=1N.D_ ξ:= \(w_i, k_i, v_i+w_i ξ) \_i=1^N\,. (73) Proposition C.3 (Gauge Invariance). For any density w∈ℋσk,σvw _ _k, _v, discrete database D, bias vector ∈ℝdout b ^d_out, and shift vector ∈ℝdout ξ ^d_out, the augmented functional is invariant under simultaneous shifts of the values and the bias: sw,+()=sw,()ands,+()=s,().s_w_ ξ, b+ ξ( q)=s_w, b( q) s_D_ ξ, b+ ξ( q)=s_D, b( q)\,. (74) In particular, any augmented state (,)(D, b) is gauge-equivalent to a canonical zero-bias state: s,()=s−,().s_D, b( q)=s_D_- b,0( q)\,. (75) Proof. For the continuous case, note that sw() s_w_ ξ( q) =∫ℝdin×ℝdoutw(,)e⟨,⟩μ(,)∫ℝdin×ℝdoutw(,)e⟨,⟩μ(,) = _R^d_in×R^d_outw_ ξ( k, v)e k, q v\,dμ( k, v) _R^d_in×R^d_outw_ ξ( k, v)e k, q \,dμ( k, v) =∫ℝdin×ℝdoutw(,−)⋅exp(2⟨,⟩−‖22σv2)e⟨,⟩exp(−‖22σk2−‖22σv2)∫ℝdin×ℝdoutw(,−)⋅exp(2⟨,⟩−‖22σv2)e⟨,⟩exp(−‖22σk2−‖22σv2) = _R^d_in×R^d_outw( k, v- ξ)· ( 2 ξ, v -\| ξ\|^22 _v^2 )e k, q v\, (- \| k\|^22 _k^2- \| v\|^22 _v^2 )d kd v _R^d_in×R^d_outw( k, v- ξ)· ( 2 ξ, v -\| ξ\|^22 _v^2 )e k, q \, (- \| k\|^22 _k^2- \| v\|^22 _v^2 )d kd v =∫ℝdin×ℝdoutw(,−)e⟨,⟩exp(−‖22σk2−‖−‖22σv2)∫ℝdin×ℝdoutw(,−)e⟨,⟩exp(−‖22σk2−‖−‖22σv2) = _R^d_in×R^d_outw( k, v- ξ)e k, q v\, (- \| k\|^22 _k^2- \| v- ξ\|^22 _v^2 )d kd v _R^d_in×R^d_outw( k, v- ξ)e k, q \, (- \| k\|^22 _k^2- \| v- ξ\|^22 _v^2 )d kd v =∫ℝdin×ℝdoutw(,)e⟨,⟩(+)exp(−‖22σk2−‖22σv2)∫ℝdin×ℝdoutw(,)e⟨,⟩exp(−‖22σk2−‖22σv2) = _R^d_in×R^d_outw( k, v)e k, q ( v+ ξ)\, (- \| k\|^22 _k^2- \| v\|^22 _v^2 )d kd v _R^d_in×R^d_outw( k, v)e k, q \, (- \| k\|^22 _k^2- \| v\|^22 _v^2 )d kd v =∫ℝdin×ℝdoutw(,)e⟨,⟩(+)μ(,)∫ℝdin×ℝdoutw(,)e⟨,⟩μ(,) = _R^d_in×R^d_outw( k, v)e k, q ( v+ ξ)\,dμ( k, v) _R^d_in×R^d_outw( k, v)e k, q \,dμ( k, v) =sw()+. =s_w( q)+ ξ\,. Evaluating the augmented functional with shifted bias + b+ ξ, we get sw,+()=sw()−(+)=(sw()+)−(+)=sw()−=sw,().s_w_ ξ, b+ ξ( q)=s_w_ ξ( q)-( b+ ξ)=(s_w( q)+ ξ)-( b+ ξ)=s_w( q)- b=s_w, b( q)\,. (76) Similarly, for the discrete case, we have: s() s_D_ ξ( q) =∑i=1Ne⟨i,⟩(~i+wi)∑i=1Nwie⟨i,⟩ = _i=1^Ne k_i, q ( v_i+w_i ξ) _i=1^Nw_ie k_i, q =∑i=1Ne⟨i,⟩~i∑i=1Nwie⟨i,⟩+∑i=1Nwie⟨i,⟩∑i=1Nwie⟨i,⟩ = _i=1^Ne k_i, q v_i _i=1^Nw_ie k_i, q + ξ _i=1^Nw_ie k_i, q _i=1^Nw_ie k_i, q =s()+. =s_D( q)+ ξ\,. (77) Evaluating the augmented functional with shifted bias + b+ ξ, we get s,+()=s()−(+)=(s()+)−(+)=s()−=s,().s_D_ ξ, b+ ξ( q)=s_D_ ξ( q)-( b+ ξ)=(s_D( q)+ ξ)-( b+ ξ)=s_D( q)- b=s_D, b( q)\,. (78) Choosing =− ξ=- b yields s,()=s−,()s_D, b( q)=s_D_- b,0( q), completing the proof. ∎ C.7.2 Gauge Fixing for Strictly Positive Weight Updates Suppose at step t the model is stored in the canonical zero-bias frame (t)D^(t) with = b=0, producing prediction t=s(t)(t) s_t=s_D^(t)( q_t) on query t q_t and loss gradient t=∇^ℒ(t,t) g_t= _ yL( s_t, y_t). Instead of computing the functional gradient in the canonical frame (which risks wnew<0w_new<0), we perform a gauge transformation to view the database in a virtual gauge with bias t∈ℝdout b_t ^d_out. By Proposition C.3, mapping the canonical zero-bias state to a frame with bias t b_t requires shifting the values by =t ξ= b_t, resulting in the shifted virtual database: t(t)=(wi,i,~i+wit)i=1N.D_ b_t^(t)= \(w_i, k_i, v_i+w_i b_t) \_i=1^N\,. (79) Recall that st(t),t(t)=st(t)(t)−t=s(t)(t)+t−t=s(t)(t)s_D_ b_t^(t), b_t( q_t)=s_D_ b_t^(t)( q_t)- b_t=s_D^(t)( q_t)+ b_t- b_t=s_D^(t)( q_t). Now consider adding a tuple (wnew,new,~new′)(w_new, k_new, v _new) to t(t)D_ b_t^(t) in this frame. As derived in Sec. C.4, the functional gradient descent update on the weight density w yields a scalar mass update proportional to the projection of the gradient onto the effective prediction relative to the frame bias: wnew=η⟨t,st(t)(t)⟩=η⟨t,s(t)(t)+t⟩=η⟨t,t+t⟩.w_new=η g_t,s_D_ b_t^(t)( q_t) =η g_t,s_D^(t)( q_t)+ b_t =η g_t, s_t+ b_t \,. (80) Because the gauge choice t b_t is a free degree of freedom that leaves the function invariant, we can fix the gauge dynamically at each step to ensure that the inner product in equation 80 is strictly positive. Specifically, we choose: t=(−⟨t,t⟩/∥t∥2+αt)t, b_t= (- s_t, g_t /\| g_t\|^2+ _t ) g_t\,, (81) where αt>0 _t>0 is a scalar gauge parameter. Substituting equation 81 into equation 80 yields: wnew=η⟨t,t+(−⟨t,t⟩/∥t∥2+αt)t⟩=ηαt∥t∥2>0.w_new=η g_t, s_t+ (- s_t, g_t /\| g_t\|^2+ _t ) g_t =η _t\| g_t\|^2>0\,. (82) Thus, for any non-zero error gradient t g_t, the new weight is guaranteed to be strictly positive. The target weighted value vector ~new′ v _new in the virtual frame is directed along the negative error gradient, following standard functional descent on the numerator functional: ~new′=−βt, v _new=-β g_t\,, (83) where β>0β>0 is the value learning rate. The key is set as usual to new=Σkt=σk2t k_new= _k q_t= _k^2 q_t. Appending this tuple to t(t)D_ b_t^(t) yields the preliminary database prel,(t+1)D^prel,(t+1) in the virtual frame: prel,(t+1):=t(t)∪(wnew,new,~new′).D^prel,(t+1):=D_ b_t^(t)∪\(w_new, k_new, v _new)\\,. C.7.3 Canonical Return and Locality of the Database Update To maintain a standard, zero-bias database across training steps, we must transform prel,(t+1)D^prel,(t+1) back to the canonical gauge (= b=0). Shifting the bias from t b_t back to 0 corresponds to applying the shift =−t ξ=- b_t to all values. A naive implementation might suggest that shifting values by −t- b_t requires modifying every existing entry in the database. Remarkably, this is not the case: the gauge transformation yields a local update. Theorem C.4 (Locality of the Gauge-Invariant Update). The transition from the canonical database (t)D^(t) to (t+1)D^(t+1) consists solely of appending a single new tuple: (t+1)=(t)∪(wnew,new,~final),D^(t+1)=D^(t)∪ \(w_new, k_new, v_final) \\,, (84) where all previously stored tuples remain completely unaltered, and the canonical (weighted) value vector of the new tuple is given by: ~final=(−wnewαt−β+ηαt⟨t,t⟩)t. v_final= (-w_new _t-β+η _t s_t, g_t ) g_t\,. (85) Proof. Let (wi,i,~i)∈(t)(w_i, k_i, v_i) ^(t) be an arbitrary pre-existing tuple. 1. Forward transformation to virtual gauge: Under shift =+t ξ=+ b_t, the weighted value becomes ~i′=~i+wit v_i = v_i+w_i b_t. 2. Atom insertion: The new atom (wnew,new,~new′)(w_new, k_new, v _new) is added. Existing tuples are unaffected. 3. Return transformation to canonical gauge: Shifting the bias by −t- b_t transforms all weighted values via ~↦~−wt v v-w b_t. For any existing tuple i: ~i′=~i′−wit=(~i+wit)−wit=~i. v_i = v_i -w_i b_t=( v_i+w_i b_t)-w_i b_t= v_i\,. (86) All previously stored tuples return identically to their initial values. For the newly added atom, its weighted value in the virtual frame was ~new′=−βt v _new=-β g_t. Shifting it to the canonical frame yields: ~final v_final =~new′−wnewt = v _new-w_new b_t =−βt−wnew(−⟨t,t⟩/∥t∥2+αt)t =-β g_t-w_new (- s_t, g_t /\| g_t\|^2+ _t ) g_t =(−wnewαt−β+wnew⟨t,t⟩‖t‖2)t = (-w_new _t-β+w_new s_t, g_t \| g_t\|^2 ) g_t (87) =(−wnewαt−β+ηαt‖t‖2⟨t,t⟩‖t‖2)t = (-w_new _t-β+η _t\| g_t\|^2 s_t, g_t \| g_t\|^2 ) g_t (88) =(−wnewαt−β+ηαt⟨t,t⟩)t. = (-w_new _t-β+η _t s_t, g_t ) g_t\,. (89) Thus, no historical data points are modified, and the update is strictly local. ∎ C.7.4 Algorithm Summary The complete procedure is formalized in Algorithm 1. As mentioned in Remark 3.1, we store new←t k_new← q_t in the database instead of new←σk2t k_new← _k^2 q_t. Algorithm 1 Gauge-Invariant Weighted Softmax Growing Step 1: Input: Current database (t)D^(t), input query t q_t, loss gradient t=∇^ℒ(t,t) g_t= _ yL( s_t, y_t) 2: Hyperparameters: Weight learning rate η>0η>0, value learning rate β>0β>0, gauge parameter αt>0 _t>0, temperature T 3: 1. Compute prediction: 4: t←∑(w,,~)∈(t)e1Tdin⟨,t⟩~∑(w,,~)∈(t)e1Tdin⟨,t⟩w s_t← _(w, k, v) ^(t)e 1T d_in k, q_t v _(w, k, v) ^(t)e 1T d_in k, q_t w 5: 2. Compute strictly positive weight update: 6: wnew←ηαt‖t‖2w_new←η _t\| g_t\|^2 7: 3. Compute canonical weighted value: 8: ~final←(−wnewαt−β+ηαt⟨t,t⟩)t v_final← (-w_new _t-β+η _t s_t, g_t ) g_t 9: 4. Set new key: 10: new←t k_new← q_t 11: 5. Append new tuple to database: 12: (t+1)←(t)∪(wnew,new,~final)D^(t+1) ^(t)∪ \(w_new, k_new, v_final) \ C.7.5 Hyperparameter Choice that Yields Changes in the Order of the Gradient Length The weight of the added tuple is wnew=ηαt‖t‖2w_new=η _t\| g_t\|^2 which is quadratic in the length of the gradient. In order to ensure that we add updates that are in the order of the gradient length, we can set αt:=α‖t‖ _t:= α\| g_t\|. In this case, we get wnew=ηα‖t‖,w_new=ηα\| g_t\|\,, and ~final v_final =(−wnewαt−β+ηαt⟨t,t⟩)t = (-w_new _t-β+η _t s_t, g_t ) g_t =(−ηαt‖t‖2αt−β+ηαt⟨t,t⟩)t = (-η _t\| g_t\|^2 _t-β+η _t s_t, g_t ) g_t =(−ηα‖t‖‖t‖2α‖t‖−β+ηα‖t‖⟨t,t⟩)t = (-η α\| g_t\|\| g_t\|^2 α\| g_t\|-β+η α\| g_t\| s_t, g_t ) g_t =(−ηα2−β+ηα‖t‖⟨t,t⟩)t. = (-ηα^2-β+η α\| g_t\| s_t, g_t ) g_t\,. This yields the following algorithm. Algorithm 2 Gauge-Invariant Weighted Softmax Learning with Updates in the Order of the Gradient Length 1: Input: Current database (t)D^(t), input query t q_t, loss gradient t=∇^ℒ(t,t) g_t= _ yL( s_t, y_t) 2: Hyperparameters: Weight learning rate η>0η>0, value learning rate β>0β>0, gauge parameter α>0α>0, temperature T 3: 1. Compute prediction: 4: t←∑(w,,~)∈(t)e1Tdin⟨,t⟩~∑(w,,~)∈(t)e1Tdin⟨,t⟩w s_t← _(w, k, v) ^(t)e 1T d_in k, q_t v _(w, k, v) ^(t)e 1T d_in k, q_t w 5: 2. Compute strictly positive weight update: 6: wnew←ηα‖t‖w_new←ηα\| g_t\| 7: 3. Compute canonical weighted value: 8: ~final←(−ηα2−β+ηα‖t‖⟨t,t⟩)t v_final← (-ηα^2-β+η α\| g_t\| s_t, g_t ) g_t 9: 4. Set new key: 10: new←t k_new← q_t 11: 5. Append new tuple to database: 12: (t+1)←(t)∪(wnew,new,~final)D^(t+1) ^(t)∪ \(w_new, k_new, v_final) \ Appendix D Optimization Dynamics in the Growing Database Standard deep learning optimization techniques such as weight decay and momentum can be naturally incorporated into the growing database framework. Weight decay, which corresponds to L2L_2 regularization in the functional space, translates to a simple uniform decay of the scalar weights w and weighted values ~ v of all stored key-value pairs at each step. This mechanism not only regularizes the model but also introduces a “forgetting” factor where older memories gradually fade in influence. As we show in Sec. D.1, weight-decay corresponds to the update rule t= _t= ((1−γ)w,,(1−γ)~)∣(w,,~)∈t−1 \((1-γ)w, k,(1-γ) v) (w, k, v) _t-1\ ∪(wnew,new,~new). ∪\,\(w_new, k_new, v_new)\\,. (90) This means that the weights of old tuples decay exponentially with training steps, thereby providing a principled path towards “compression” of the database: Old tuples with negligible weights can be pruned with minimal impact on the function, potentially allowing the database to be maintained at a constant size. Momentum, typically implemented by maintaining a velocity vector, would naively require doubling the storage (storing a velocity database alongside the parameter database). However, we derive a more efficient formulation where momentum is realized analytically. By storing the creation timestamp tit_i of each key, we can compute the effective momentum contribution of a key during the forward pass via a time-dependent scaling coefficient: c(t,ti)=1−μt−ti+11−μ,c(t,t_i)= 1-μ^t-t_i+11-μ\,, (91) where μ is the momentum factor. As we show in Sec. D.2, this allows the growing network to benefit from accelerated convergence without the storage overhead of explicit velocity buffers. One can also combine both approaches. The resulting algorithm applies the geometric decay to the stored weights at every step to handle regularization, while applying the time-dependent momentum coefficient during the retrieval operation to simulate the accumulation of gradients. We provide detailed derivations of these dynamics in Sec. D.3. In this appendix, we derive the update rules for weight decay and momentum within the growing softmax framework. We assume the decoupled representation where the database stores (w,,~)(w, k, v), and the functions are wD()=∑wie⟨,i⟩w_D( k)=Σ w_ie k, k_i and wN()=∑~ie⟨,i⟩w_N( k)=Σ v_ie k, k_i . Note that here we use ~ v to denote the weighted value vector, consistent with the main text notation. D.1 Weight Decay Weight decay corresponds to minimizing a regularized objective: Jreg(wN,wD)=J(wN,wD)+γ2η(‖wN‖ℋ2+‖wD‖ℋ2).J_reg(w_N,w_D)=J(w_N,w_D)+ γ2η(\|w_N\|_H^2+\|w_D\|_H^2)\,. Here, we adopt the “decoupled weight” formulation of Sec. C.6. Nevertheless, following a “coupled weight” formulation as in equation 5 yields similar results. The functional gradient descent update for wDw_D becomes: wD,t←wD,t−1−η(∇wDJ+γηwD,t−1)=(1−γ)wD,t−1−η∇wDJ.w_D,t← w_D,t-1-η( _w_DJ+ γηw_D,t-1)=(1-γ)w_D,t-1-η _w_DJ\,. Similarly for the numerator function wNw_N: wN,t←(1−γ)wN,t−1−η∇wNJ.w_N,t←(1-γ)w_N,t-1-η _w_NJ\,. In the database representation, the function wD,t−1w_D,t-1 is represented by the sum of existing kernels with weights wiw_i. Multiplying the function by (1−γ)(1-γ) is equivalent to multiplying every coefficient wiw_i in the database by (1−γ)(1-γ). The same logic applies to the numerator function and its coefficients ~i v_i. Thus, the update rule for the database t−1D_t-1 is: t=((1−γ)w,,(1−γ)~)∣(w,,~)∈t−1∪(wnew,new,~new),D_t=\((1-γ)w, k,(1-γ) v) (w, k, v) _t-1\∪\(w_new, k_new, v_new)\\,, where the new tuple is derived from the gradients −η∇wDJ-η _w_DJ and −η∇wNJ-η _w_NJ as in Sec. C.4. This operation implies that the weight of a tuple added at time t0t_0 decays geometrically as (1−γ)t−t0(1-γ)^t-t_0 at step t. Consequently, older tuples will eventually have weights strictly bounded by some small ϵε. This suggests a natural pruning strategy: removing tuples where |w|<ϵ|w|<ε and ‖~‖<ϵ\| v\|<ε results in a negligible perturbation to the softmax function output. By setting a threshold, one can effectively bound the maximum size of the database, converting the growing nonparametric model into a fixed-capacity model with a dynamic memory. D.2 Momentum Momentum accelerates optimization by accumulating a velocity function ν∈ℋν . We apply momentum to the gradient updates for both wNw_N and wDw_D. Let’s focus on wDw_D (the derivation for wNw_N is identical). The update rules are: νt=μνt−1+∇wDJt,wD,t=wD,t−1−ηνt. _t=μ _t-1+ _w_DJ_t, w_D,t=w_D,t-1-η _t\,. Unrolling the recursion, the weight function at step t can be expressed as a linear combination of all past gradients: wD,t=wD,0−η∑j=1t(∑k=0t−jμk)∇wDJj=wD,0−η∑j=1t1−μt−j+11−μ∇wDJj.w_D,t=w_D,0-η _j=1^t ( _k=0^t-jμ^k ) _w_DJ_j=w_D,0-η _j=1^t 1-μ^t-j+11-μ _w_DJ_j\,. Recall that the gradient ∇wDJj _w_DJ_j at step j corresponds to adding a new kernel centered at j=σ2j k_j=σ^2 q_j with a specific weight Δwj w_j. The unrolled equation implies that the accumulated weight for the kernel at j k_j at current time t is simply the initial update weight Δwj w_j scaled by the coefficient: c(t,tj)=1−μt−tj+11−μ.c(t,t_j)= 1-μ^t-t_j+11-μ\,. Therefore, we can implement momentum without storing the velocity function νt _t explicitly. Instead, we store the creation time tjt_j for each tuple. During the retrieval (forward pass) at current time t, the contribution of the j-th tuple to the denominator is scaled by c(t,tj)c(t,t_j). The same scaling applies to the weighted value vector ~j v_j in the numerator. D.3 Combined Dynamics In this section, we provide a rigorous derivation of the update rule when both weight decay and momentum are active. We adopt the decoupled weight decay formulation [45], which is the standard for modern deep learning optimizers like AdamW and SGDW. D.3.1 Setup and Update Equations Let wt∈ℋw_t denote the function state at time t. Let ∇Jt∇ J_t denote the functional gradient of the loss at step t. The update dynamics are defined by the following system of recurrence relations: νt _t =μνt−1+∇Jt, =μ _t-1+∇ J_t\,, (92) wt w_t =(1−γ)wt−1−ηνt. =(1-γ)w_t-1-η _t\,. (93) where: • μ∈[0,1)μ∈[0,1) is the momentum coefficient. • γ∈[0,1)γ∈[0,1) is the weight decay factor per step. • η>0η>0 is the learning rate. • νt _t is the velocity function (initialized at ν0=0 _0=0). • w0w_0 is the initial function state. D.3.2 Derivation of the Closed-Form State We aim to express wtw_t as a direct summation of the history of gradients ∇Jjj=1t\∇ J_j\_j=1^t. For convenience, we define the effective momentum μ~ μ and effective learning rate η~ η as: μ~=μ1−γ,η~=η1−γ. μ= μ1-γ, η= η1-γ\,. First, we unroll the velocity term νt _t using equation 92: νk=∑j=1kμk−j∇Jj. _k= _j=1^kμ^k-j∇ J_j\,. (94) Next, we substitute this into the weight update equation 93. By unrolling the recurrence for wtw_t, we obtain: wt w_t =(1−γ)wt−1−ηνt =(1-γ)w_t-1-η _t =(1−γ)tw0−η∑k=1t(1−γ)t−kνk =(1-γ)^tw_0-η _k=1^t(1-γ)^t-k _k =(1−γ)tw0−η~∑k=1t(1−γ)t−k+1νk. =(1-γ)^tw_0- η _k=1^t(1-γ)^t-k+1 _k\,. (95) Substituting the expression for νk _k and rearranging the summation order (1≤j≤k≤t1≤ j≤ k≤ t): wt=(1−γ)tw0−η~∑j=1t∇Jj(∑k=jt(1−γ)t−k+1μk−j).w_t=(1-γ)^tw_0- η _j=1^t∇ J_j ( _k=j^t(1-γ)^t-k+1μ^k-j )\,. (96) We simplify the inner summation (t,j)C(t,j) as follows: (t,j) (t,j) =∑k=jt(1−γ)t−k+1μk−j = _k=j^t(1-γ)^t-k+1μ^k-j (97) =∑m=0t−j(1−γ)(t−j+1)−mμm = _m=0^t-j(1-γ)^(t-j+1)-mμ^m =(1−γ)t−j+1∑m=0t−j(μ1−γ)m =(1-γ)^t-j+1 _m=0^t-j ( μ1-γ )^m =(1−γ)t−j+1∑m=0t−jμ~m =(1-γ)^t-j+1 _m=0^t-j μ^m (98) =(1−γ)t−j+11−μ~t−j+11−μ~. =(1-γ)^t-j+1 1- μ^t-j+11- μ\,. (99) Substituting this back into the expression for wtw_t, we arrive directly at the compact closed-form solution: wt=(1−γ)tw0−η~∑j=1t[(1−γ)t−j+11−μ~t−j+11−μ~]∇Jj.w_t=(1-γ)^tw_0- η _j=1^t [(1-γ)^t-j+1 1- μ^t-j+11- μ ]∇ J_j\,. (100) Similarly to the implementation of momentum, we can store the creation time tjt_j for each tuple and then take it into account in the forward pass. Appendix E Additional technical details on advanced linear attention layers E.1 DeltaNet In (unnormalized) linear attention, each step writes a new key-value association into a fixed size state via an outer-product update. The delta rule instead updates the value associated with the current key by (partly) overwriting the previously stored value with write strength λ∈[0,1],tnew=λt+(1−λ)toldλ∈[0,1], v^new_t=λ v_t+(1-λ) v_t^old and t=t−1−toldt⊤+tnewt⊤ W_t= W_t-1- v_t^old k_t + v_t^new k_t . In our setting, we set the key to be the activation, i.e., t:=t k_t:= x_t and define the previously stored value for that key as the current prediction t−1t W_t-1 x_t. To compress a mapping from activations to (negative) gradients into a fixed size state as in vanilla gradient descent we set the target value to be the negative output-gradient scaled by a step size η, i.e., t:=−η∇tℒ v_t:=-η _ y_tL, which gives tnew=−λη∇tℒ+(1−λ)t−1t. v^new_t=-λη _ y_tL+(1-λ) W_t-1 x_t. This yields the recurrent update t=t−1(−λtt⊤)+(−λη∇tℒ)t⊤, W_t= W_t-1( I-λ x_t x_t )+(-λη _ y_tL) x_t , (101) which shows that a DeltaNet-like update can be interpreted as gradient descent preceded by a targeted forgetting mechanism along the current activation direction. Notably, the update itself corresponds to performing a gradient descent step w.r.t. W on the layer-local, time-dependent loss ℓt():=12‖t−t‖22=12‖t+η∇tℒ‖22 _t( W):= 12\| W k_t- v_t\|_2^2= 12\| W x_t+η _ y_tL\|_2^2 (102) with learning rate λ [79, 67], meaning the Delta-rule-based update implicitly performs gradient descent within the gradient descent weight update. In practice, the update can be implemented with little computational overhead by rewriting Eq. 101 to t=t−1+λ(−η∇tℒ−t−1t)t⊤, W_t= W_t-1+λ(-η _ y_tL- W_t-1 x_t) x_t , (103) caching the activation t−1t W_t-1 x_t from the forward pass and using it to form the corrected term −η∇tℒ−t−1t-η _ y_tL- W_t-1 x_t before computing the outer-product update. E.2 MesaNet: Problem formulation and connection to FOOF Next, we show the connection between FOOF and Mesa layer when applied to the database. We present a slightly modified formulation of the Mesa layer, with decay factor and an additional regularization center. We present the case of batch size 1, but the result can easily be extended to the general case. Consider the online estimation of a parameter matrix ∈ℝk×d W ^k× d from a stream of observations (t,t)( x_t, y_t) (key and value of the database) for t=1,…,Tt=1,…,T. The mesa layer seeks to minimize the following time-discounted objective function JT()J_T( W) at each step T: JT()=∑t=1TγT−t‖t−t‖2+ΛT‖−T‖F2J_T( W)= _t=1^Tγ^T-t\| y_t- W x_t\|^2+ _T\| W- C_T\|_F^2 (104) where: • γ∈[0,1]γ∈[0,1] is the forgetting factor. • ΛT>0 _T>0 is the regularization strength at time T. • T∈ℝk×d C_T ^k× d is the regularization center (prior mean) at time T, taken as the exponential moving average of past weights, defined as T=γΛT−1ΛTT−1+(1−γΛT−1ΛT)T−1 C_T=γ _T-1 _T C_T-1+(1-γ _T-1 _T) W_T-1, and 0= C_0=0. The loss in the main text is recovered by choosing γ=1γ=1 and Λt=λ∀t _t=λ∀ t. It can be shown that the optimal solution T W_T satisfies the normal equation: TΦT=T W_T _T= P_T (105) where ΦT _T and T P_T are the regularized correlation matrix and cross-correlation vector, defined respectively as: ΦT _T =∑t=1TγT−ttt⊤+ΛT = _t=1^Tγ^T-t x_t x_t + _T I (106) T P_T =∑t=1TγT−ttt⊤+ΛTT = _t=1^Tγ^T-t y_t x_t + _T C_T (107) E.3 Derivation of the Recursive Update We seek to express T W_T as a function of T−1 W_T-1 and the current observation (T,T)( x_T, y_T). First, we expand the recursive definitions of the sufficient statistics. E.3.1 Recursive Matrix Update Separating the term at t=Tt=T and substituting ΦT−1 _T-1: ΦT _T =γ(∑t=1T−1γT−1−ttt⊤)+TT⊤+ΛT =γ ( _t=1^T-1γ^T-1-t x_t x_t )+ x_T x_T + _T I (108) =γ(ΦT−1−ΛT−1)+TT⊤+ΛT =γ( _T-1- _T-1 I)+ x_T x_T + _T I =γΦT−1+TT⊤+(ΛT−γΛT−1) =γ _T-1+ x_T x_T +( _T-γ _T-1) I E.3.2 Recursive Vector Update Similarly for PTP_T: T P_T =γ(T−1−ΛT−1T−1)+TT⊤+ΛTT =γ( P_T-1- _T-1 C_T-1)+ y_T x_T + _T C_T (109) =γT−1+TT⊤+(ΛTT−γΛT−1T−1) =γ P_T-1+ y_T x_T +( _T C_T-γ _T-1 C_T-1) E.3.3 Sufficient Condition for the Delta Rule Assuming T−1=T−1ΦT−1−1 W_T-1= P_T-1 _T-1^-1, we substitute T−1=T−1ΦT−1 P_T-1= W_T-1 _T-1 into the expression for T P_T: T=T−1(γΦT−1)+TT⊤+(ΛTT−γΛT−1T−1) P_T= W_T-1(γ _T-1)+ y_T x_T +( _T C_T-γ _T-1 C_T-1) (110) From the matrix update step, we substitute γΦT−1=ΦT−TT⊤−(ΛT−γΛT−1)γ _T-1= _T- x_T x_T -( _T-γ _T-1) I: T=T−1(ΦT−TT⊤−(ΛT−γΛT−1))+TT⊤+(ΛTT−γΛT−1T−1) P_T= W_T-1 ( _T- x_T x_T -( _T-γ _T-1) I )+ y_T x_T +( _T C_T-γ _T-1 C_T-1) (111) Rearranging terms to isolate the prediction error: T=T−1ΦT+(T−T−1T)T⊤+ℛ P_T= W_T-1 _T+( y_T- W_T-1 x_T) x_T +R (112) where the residual term ℛR is defined as: ℛ=(ΛTT−γΛT−1T−1)−T−1(ΛT−γΛT−1)R=( _T C_T-γ _T-1 C_T-1)- W_T-1( _T-γ _T-1) (113) Multiplying by ΦT−1 _T^-1 from the right yields the update for T W_T: T=T−1+(T−T−1T)T⊤ΦT−1+ℛΦT−1 W_T= W_T-1+( y_T- W_T-1 x_T) x_T _T^-1+R _T^-1 (114) For the residual term ℛR to vanish, it is necessary and sufficient to impose the following consistency condition on the update of the regularization center T C_T: ΛTT=γΛT−1T−1+(ΛT−γΛT−1)T−1, _T C_T=γ _T-1 C_T-1+( _T-γ _T-1) W_T-1, (115) which is satisfied by definition. Under the above condition, and using as target at every timestep T=T−1T−η∇Tℒ y_T= W_T-1 x_T-η _ y_TL, applying the loss in Equation 104 to the database would correspond to having a state, ,Φ W, , which is updated following: ΦT _T =γΦT−1+TT⊤+(ΛT−γΛT−1) =γ _T-1+ x_T x_T +( _T-γ _T-1) I (116) T W_T =T−1−η∇T−1ℒ⋅ΦT−1 = W_T-1-η _ W_T-1L· _T^-1 (117) It now only remains to define a meaningful ΛT _T. E.4 Proposed Algorithms We identify two distinct regimes satisfying the consistency condition equation 115. E.4.1 Algorithm A: Static Regularization We assume a constant regularization strength ΛT=λ _T=λ for all T. Matrix Update: ΦT=γΦT−1+TT⊤+λ(1−γ) _T=γ _T-1+ x_T x_T +λ(1-γ) I Which corresponds exactly to the FOOF optimizer. E.4.2 Algorithm B: Cumulative Regularization The previous algorithm causes growing norm of Φ and a decrease in the importance of the prior, which may incur numerical and hyperparameter instabilities. We therefore propose the next algorithm, which we use in practice. We assume the regularization strength accumulates over time according to ΛT=λ+γΛT−1 _T=λ+γ _T-1. This is equivalent to optimizing the weighted averaged loss JT()=1∑t=1TγT−t∑t=1TγT−t‖t−t‖2+λ‖−T‖F2J_T( W)= 1 _t=1^Tγ^T-t _t=1^Tγ^T-t\| y_t- W x_t\|^2+λ\| W- C_T\|_F^2 (118) Matrix Update: ΦT=γΦT−1+TT⊤+λ _T=γ _T-1+ x_T x_T +λ I While keeping the above objective, we introduce additional tricks to stabilize numeric and hyperparameter across different layer widths, resulting in the following recursive update for a batch of B inputs T∈ℝd×B X_T ^d× B: ΓT _T =γΓT−1+1 =γ _T-1+1 (119) ΨT _T =γΓT−1ΓTΨT−1+1BdΓTTT⊤ = γ _T-1 _T _T-1+ 1Bd _T X_T X_T (120) T W_T =T−1−η∇T−1ℒ⋅(ΨT+λ)−1, = W_T-1-η _ W_T-1L·( _T+λ I)^-1, (121) while initializing Ψ0= _0=0, 0 W_0 to a standard weight initialization, and Γ0=0 _0=0. Appendix F Experimental Setup F.1 Tasks and Model Selection We run experiments on image classification and teacher-student regression tasks. For all experiments, we report the mean and standard deviation over 5 seeds. To produce the results in Figure 2, we run each model/optimizer with a comparable set of hyperparameters where they are shared and for each method add the needed, additional hyperparameters. For Figures 1, 4, and 4, we select a subset of the hyperparameters relevant to the produced result and again detail the selection in the respective model subsection. The considered ranges of hyperparameters are reported beneath in the respective model subsections. Generally, for each plot and for each method displayed in the respective plot, we select those hyperparameters that maximize/minimize the metric displayed in the plot. The only place where we share the selection across multiple plots is Fig. 2. There, the runs visualized in Fig. 2/Left and Fig. 2/Middle are tied in that the two metrics (Test Accuracy and Training Loss) are displayed for the same selection of runs that maximize the final test accuracy. In particular, looking for the hyperparameter that minimizes the training loss would yield a different selection of runs achieving even smaller training loss (but potentially worse generalization). In the following two subsections, we first briefly describe the data processing/generation pipeline used for classification and regression experiments. F.1.1 Image Classification CIFAR-10 and CIFAR-100 [42] are image classification datasets on 32×32×332× 32× 3 dimensional RGB images from 10 and 100 classes, respectively. Before handing these images to the network, we first apply standard normalization by subtracting the mean and dividing by the standard deviation per channel. Then, to increase the effective size of the training dataset, we apply standard image transformations, randomly sampled during training. These include reflect-padding (±4± 4 pixels) followed by a random 32×3232× 32 crop applied to all 3 channels and an epoch-wise horizontal flip. Test data is only normalized without applying augmentation. Before handing the images to our MLP network, we flatten them into a 30723072-dimensional vector. The notion of epoch in the respective experiments refers to the following: take the 50000 training samples from the respective dataset, apply a randomly chosen augmentation (following the above rules) to each sample, and ultimately shuffle the samples. Only for Fig. 2/Right, where only a single pass on the training data is performed, we use the non-augmented standard version of CIFAR-10. F.1.2 Teacher-Student Regression We generate synthetic regression tasks using a random teacher network. The inputs ∈ℝdin x ^d_in are sampled from a standard normal distribution and projected onto the unit sphere. The targets for the randomly sampled inputs are generated by a randomly initialized MLP teacher with L residual hidden blocks adding to a residual stream of width ded_e. Crucially, this teacher has nonlinearities thereby, upon initialization, implementing some nonlinear input-output(target) mapping. Particularly, the teacher first projects the input data to the hidden dimension by applying a linear layer and followed by a ReLU activation. Subsequently, L−1L-1 residual blocks, where each block consists of a linear layer and a ReLU activation, are applied to the initial residual stream activation. A final linear projection maps the last residual stream activation mapping to the output feature (the target for the regression task) of dimension doutd_out. All weights are randomly initialized (standard variance scaling init). In our experiments, we set din=100d_in=100, de=100d_e=100, L=2L=2, and dout=1d_out=1. F.2 Model Architecture Throughout the entire paper, and importantly shared across all methods, all experiments (unless explicitly stated otherwise) are run with a ResNet-MLP architecture. The inputs act as the initial residual stream activation, thereby governing the residual stream dimension de=dind_e=d_in. The network backbone consists of L residual blocks arranged in a pre-normalization configuration. Specifically, each block applies layer normalization, followed by a trainable module (a standard dense layer or one of the proposed alternatives: Delta, Mesa, RBF, or Weighted Softmax). In the case of a standard dense layer, Delta, or Mesa, a non-linear activation function is applied to the output of the block and the norm rescaling parameters are trained. Note that for RBF and Weighted Softmax, nonlinearities are by design included in the block. Also, for those we do not train the scale parameters of the norm but fix them to be unit. In all cases, the block’s output is added back to the residual stream. The final residual stream activation is processed by a final normalization (again scale parameters not trained for RBF and Softmax and trained for all others) before being mapped to the C class logits via a final output block (again, either a dense layer or one of the proposed methods). In the following subsections, we detail the various instantiations of the residual blocks with the methods described in this work. F.3 Weighted Softmax We evaluate Weighted Softmax kernel-based RCDL models on both synthetic teacher-student regression and image classification tasks, using the classic CIFAR-10 and CIFAR-100 benchmarks. For these experiments, we instantiate the blocks of the ResNet-MLP backbone architecture, described in Appendix F.2, with Weighted Softmax RCDL modules. Algorithm 3 Growing Weighted Softmax Attention 1: Input: Stream of query batches t,nn=1B\ q_t,n\_n=1^B, temperature T, value-learning rate β>0β>0, weight-learning rate η>0η>0, 0=1B W_0= 1B, 0∼(0,1din) K_0 (0, 1 d_in I), 0∼(0,1dout) V_0 (0, 1 d_out I) 2: Initialize: 0←(0,0,0)D_0←\( W_0, K_0, V_0)\ 3: for t=1,2,…t=1,2,… do 4: for n=1,…,Bn=1,…,B do 5: Compute prediction: t,n←∑(w,,~)∈t−1e1Tdin⟨,t,n⟩~∑(w,,~)∈t−1we1Tdin⟨,t,n⟩ s_t,n← _(w, k, v) _t-1e 1T d_in k, q_t,n \, v _(w, k, v) _t-1w\,e 1T d_in k, q_t,n 6: Compute gradient: t,n←∇ℒ(t,n) g_t,n← ( s_t,n) 7: Compute inner product: δt,n←⟨t,n,t,n⟩ _t,n← s_t,n, g_t,n 8: end for 9: Update database: t←t−1∪(max0,ηδt,n,t,n,−βt,n):n=1,…,BD_t _t-1∪\( \0,η\, _t,n\,\; q_t,n,\;-β\, g_t,n):n=1,…,B\ 10: end for A Weighted Softmax module maintains a growing database of key, value, weight tuples t=(t,t,t)=(i,i,wii=0t)D_t=( K_t, V_t, W_t)=(\ k_i, v_i,w_i\_i=0^t) of size (t+1)B(t+1)B at time step t, where B is the batch size. For numerical stability, we initialize this database with a randomly sampled batch at t=0t=0. Particularly, we set 0=(0,0,0)D_0=( K_0, V_0, W_0). Here, 0∈ℝB×de K_0 ^B× d_e and 0∈ℝB×dout V_0 ^B× d_out (dout=ded_out=d_e for the residual blocks) are sampled from a standard normal distribution and subsequently normalized along the feature dimension by root mean squared normalization. All entries of 0∈ℝB W_0 ^B are set to a constant weight w0=1Bw_0= 1B, where B is the batch size. At timestep t>0t>0, the database is grown following the steps outlined in Algorithm 3. Briefly, we apply the Weighted Softmax update rule but clip negative weights to 0 (cf. Appendix C.5). Moreover, for numerical stability, we drop the division by the normalization Dt−1D_t-1 for both the weight and the value update. Ultimately, we introduce learning rates η and β to scale weights and values independently. Table 2: Hyperparameters Nonparametric Growing Weighted Softmax for Teacher-Student experiments. Hyperparameter Value Model Single Weighted Softmax layer Data normalization RMSNorm (on input) Value learning rate β [0.1,0.3,1.0][0.1,0.3,1.0] Weight learning rate η [10.0,30.0,100.0][10.0,30.0,100.0] Softmax temperature T [0.1,1.0,3.0,10.0][0.1,1.0,3.0,10.0] Momentum [0.6,0.9][0.6,0.9] Batch size 6464 Epochs 200000200000 We report the hyperparameters considered to produce the Weighted Softmax numbers in Figure 1 in Table 2. Additionally, the dimensionality of the inputs and thus of the keys in the synthetic regression task is de=100d_e=100, Input=100Input=100. Table 3: Hyperparameters Growing Weighted Softmax for experiments on CIFAR-10 and CIFAR-100. Hyperparameter Value Dataset CIFAR-10, CIFAR-100 Input 32×32×332× 32× 3 images flattened to 30723072 Model Residual MLP Residual stream size (ded_e) 30723072 # residual blocks (L) [0,1,2][0,1,2] Normalization in blocks RMSNorm (pre-norm, per block) Final normalization RMSNorm before head Classifier head Weighted softmax-kernel RCDL layer Weight learning rate η [0.3,1.0,1.5,3.0][0.3,1.0,1.5,3.0] Value learning rate β [20.0,30.0,40.0,60.0][20.0,30.0,40.0,60.0] Softmax temperature T [1.0,3.0,10.0,15.0,30.0][1.0,3.0,10.0,15.0,30.0] Momentum [0.0,0.1,0.9,0.999][0.0,0.1,0.9,0.999] Batch size 6464 Epochs 5050 The hyperparameters considered to produce the Weighted Softmax numbers in Figure 2 and 4 are reported in Table 3. Both Figures display data from the same set of experiments. On top, for Figure 2/Right we also include the batch size in the sweeps and consider the values 4,16,324,16,32 and 6464. For the Weighted Softmax numbers reported in Figure 4 we again follow Table 3 but only consider models with L=0L=0 residual blocks. F.3.1 Gauge-Invariant Weighted Softmax The Gauge-Invariant Weighted Softmax module uses the same database, initialization and prediction as the Weighted Softmax module (Appendix F.3). Only the weight and value of each new tuple differ (Algorithm 2). As explained in Appendix C.7, we can choose the gauge offset t b_t so that the added weight wt=ηα∥t∥≥0w_t=ηα g_t ≥ 0 is non-negative by construction, so the clipping max0,⋅ \0,·\ of Algorithm 3 is no longer needed. As before, we drop the division by Dt−1D_t-1 in both updates. In practice we add ϵ=10−8ε=10^-8 to ∥t∥ g_t in the denominator for numerical stability. F.4 Hybrid Weighted Softmax Algorithm 4 Hybrid Weighted Softmax Attention 1: Input: Stream of query batches t,nn=1B\ q_t,n\_n=1^B, temperature T, growth threshold τ, maximal size N, value-learning rate β>0β>0, weight-learning rate η>0η>0, Adam learning rate α, 0=1B W_0= 1B, 0∼(0,1din) K_0 (0, 1 d_in I), 0∼(0,1dout) V_0 (0, 1 d_out I) 2: Initialize: 0←(0,0,0)D_0←\( W_0, K_0, V_0)\ 3: for t=1,2,…t=1,2,… do 4: for n=1,…,Bn=1,…,B do 5: Compute prediction: t,n←∑(w,,~)∈t−1e1Tdin⟨,t,n⟩~∑(w,,~)∈t−1we1Tdin⟨,t,n⟩ s_t,n← _(w, k, v) _t-1e 1T d_in k, q_t,n \, v _(w, k, v) _t-1w\,e 1T d_in k, q_t,n 6: Compute gradient: t,n←∇ℒ(t,n) g_t,n← ( s_t,n) 7: Compute inner product: δt,n←⟨t,n,t,n⟩ _t,n← s_t,n, g_t,n 8: end for 9: Select novel queries: t←n:min(⋅,,~)∈t−1‖t,n−‖2>τS_t←\n: _(·, k, v) _t-1\| q_t,n- k\|_2>τ\, keeping only the first N−|t−1|N-|D_t-1| elements if necessary 10: Parametric update: update all (w,,~)∈t−1(w, k, v) _t-1 with one Adam step (learning rate α) on 1B∑nℒ(t,n) 1B _nL( s_t,n) 11: Update database: t←t−1∪(max0,ηδt,n,t,n,−βt,n):n∈tD_t _t-1∪\( \0,η\, _t,n\,\; q_t,n,\;-β\, g_t,n):n _t\ 12: end for The Hybrid Weighted Softmax layer keeps the database D of the Weighted Softmax (Appendix F.3), but differs in two ways (Algorithm 4). First, all stored tuples are treated as parameters and trained by backpropagation with standard gradient descent optimization. Second, a query is only added if its Euclidean distance to the nearest stored key exceeds the growth threshold τ. The new tuple follows the Weighted Softmax rule of Algorithm 3, with δt,n _t,n clipped to [0,∞)[0,∞) for stability. 0D_0 consists of B tuples with zero keys and values and weight 1B 1B. Because queries are RMS-normalized, ‖2=din≈55.4\| q\|_2= d_in≈ 55.4, which sets the scale of τ. We define capacity as the total number of stored tuples after training, summed over all layers and divided by the maximum 4N4N. Table 4: Hyperparameters for the Hybrid Growing Weighted Softmax on CIFAR-10 (Figure 5). Hyperparameter Value Model L=3L=3 residual blocks + readout layer, no nonlinearity Embedding dimension ded_e 30723072 (flattened image, no input embedding) Normalization RMSNorm (no learned scale) before each layer and the readout Maximal width per layer N 5001650016 Growth threshold τ [35,40,45,50,55,60][35,40,45,50,55,60] Value learning rate β [10,30,60][10,30,60] Weight learning rate η 1.01.0 Softmax temperature T [1.0,5.0,10.0][1.0,5.0,10.0] Adam learning rate α [0.0003,0.001,0.003][0.0003,0.001,0.003] Batch size 1616 Epochs 11 Table 4 lists the hyperparameters for Figure 5. F.5 Parametric Softmax Here we describe the Parametric Softmax baseline that is referenced in Figures 1 and 4. The motivation for this baseline is to compare growing a (weighted) softmax kernel-based module with learning a fixed size key-value parameterization for a softmax kernel where all keys and values are treated as parameters i.e., are optimized jointly to minimize the loss. Unlike the growing module, this baseline uses a static architecture with fixed size parameters (keys, values and weights), pre-allocated with enough capacity to match the size of a fully grown database from the RCDL case from the beginning. We obtain predictions t−1 s_t-1 applying a softmax to this huge database, but instead of growing the database step-wise, update the entire key-value parameterization by gradient descent at every timestep. Essentially, this baseline can be seen as wide MLP with a softmax nonlinearity. Table 5: Hyperparameters Parametric Softmax for experiments on CIFAR-10. Hyperparameter Value Dataset CIFAR-10 Input 32×32×332× 32× 3 images flattened to 30723072 Model MLP Normalization in blocks RMSNorm (pre-norm, per block) WKeyW_Key shape (2560000, 3072) WValueW_Value shape (2560000, 10) Nonlinear Activation Softmax Learning rate α [0.003,0.01,0.03,0.1][0.003,0.01,0.03,0.1] Softmax temperature T [0.001,0.03,0.1,0.3,1.0,3.0,10.0][0.001,0.03,0.1,0.3,1.0,3.0,10.0] Batch size 6464 Epochs 20002000 We report the hyperparameters considered to produce Parametric Softmax numbers in CIFAR-10 experiments (cf. Figure 4) in Table 5. Table 6: Hyperparameters Parametric Softmax for Teacher-Student experiments. Hyperparameter Value Model Single Parametric Softmax layer Hidden size 1024000010240000 Data normalization RMSNorm (on input) Learning rate η [0.0003,0.001,0.003,0.01,0.03][0.0003,0.001,0.003,0.01,0.03] Temperature T [0.001,0.01,0.1,0.3,1.0,3.0,10.0,30.0,100.0][0.001,0.01,0.1,0.3,1.0,3.0,10.0,30.0,100.0] Optimizer Adam Momentum 0.850.85 Batch size 6464 Training steps 200000200000 For the synthetic Teacher-Student regression task presented in Figure 1, we adapt the "database size" of our fixed parametric representations to shapes (10240000,din)(10240000,d_in) for Key W_Key and (10240000,dout)(10240000,d_out) for Value W_Value. We report the remaining hyperparameters for this result in Table 6. We find that training the weights for a weighted version of softmax leads to severe training instabilities in this parametric setup and hence report the numbers for an unweighted softmax variant. F.6 RBF Algorithm 5 Growing RBF 1: Input: Stream of query batches t,nn=1B\ q_t,n\_n=1^B, temperature T, value-learning rate η>0η>0 2: Initialize: 0←∅D_0← 3: for t=1,2,…t=1,2,… do 4: for n=1,…,Bn=1,…,B do 5: Compute prediction: t,n←∑(,)∈t−1exp(−‖t,n−‖222Tdin) s_t,n← _( k, v) _t-1 \! (- \| q_t,n- k\|_2^22T d_in ) v 6: Compute gradient: t,n←∇ℒ(t,n) g_t,n← ( s_t,n) 7: end for 8: Update database: t←t−1∪(t,n,−ηt,n):n=1,…,BD_t _t-1∪\( q_t,n,\;-η\, g_t,n):n=1,…,B\ 9: end for We evaluate RBF kernel-based RCDL models on image classification tasks, using the classic CIFAR-10 and CIFAR-100 benchmarks. For these experiments, we instantiate the blocks of the ResNet-MLP backbone architecture, described in Appendix F.2, with RBF RCDL modules. In all RBF kernel-based RCDL modules, we predict and grow the database according to Algorithm 5. Table 7: Hyperparameters RBF for experiments on CIFAR-10 and CIFAR-100. Hyperparameter Value Dataset CIFAR-10, CIFAR-100 Input 32×32×332× 32× 3 images flattened to 30723072 Model Residual MLP Residual stream size (ded_e) 30723072 # residual blocks (L) [0,1,2][0,1,2] Normalization in blocks RMSNorm (pre-norm, per block) Final normalization RMSNorm before head Classifier head Weighted RBF-kernel RCDL layer Value learning rate η [0.3,1.0,3.0][0.3,1.0,3.0] Temperature T [0.5,1.0,1.5,5.0,15.0][0.5,1.0,1.5,5.0,15.0] Momentum [0.0,0.1,0.9,0.999][0.0,0.1,0.9,0.999] Batch size 6464 Epochs 5050 We report the hyperparameters considered to produce the RBF numbers in Figure 2 in Table 7. On top, for Figure 2/Right we also include the batch size in the sweeps and consider the values 4,16,324,16,32 and 6464. F.7 Mesa (FOOF) We evaluate the performance of the Mesa update rule introduced in Section 4 on CIFAR-10 and CIFAR-100 benchmarks. For these experiments, we instantiate the blocks of the ResNet-MLP backbone architecture, described in Appendix F.2, with vanilla dense linear layers. We train the parameters of the resulting network by applying the Mesa weight update detailed in Section 4. Table 8: Best Average Accuracy on CIFAR-10 and CIFAR-100 across different numbers of steps for the newton iteration algorithm for approximating the inverse. We find that generalization performance on CIFAR saturates with 1010 iterations. Mesa newton-approximation steps Dataset 5 10 Inf. CIFAR-10 0.6914 0.7130 0.7124 CIFAR-100 0.4149 0.4323 0.4319 Table 9: Hyperparameters Mesa for experiments on CIFAR-10 and CIFAR-100. Hyperparameter Value Dataset CIFAR-10, CIFAR-100 Input 32×32×332× 32× 3 images flattened to 30723072 Model Residual MLP Residual stream size (ded_e) 30723072 # residual blocks (L) [2,4][2,4] Normalization in blocks RMSNorm (pre-norm, per block) Final normalization RMSNorm before head MLP activation GELU Classifier head Linear (3072→10)(3072→ 10) Mesa-γ [0.9,0.99,1.0][0.9,0.99,1.0] Mesa-λ [0.001,0.01,0.1][0.001,0.01,0.1] Optimizer Adam Learning rate α [0.0001,0.0003,0.001,0.003,0.01][0.0001,0.0003,0.001,0.003,0.01] Momentum 0.850.85 Batch size 6464 Epochs 5050 We report the hyperparameters used to produce the numbers for the Mesa update rule in Figure 2 in Table 9. We implement the computation of the matrix inverse (cf. Appendix E) using the Newton-Schulz Iteration. We report the effect of iterations for this algorithm in Table 8. F.8 Delta Lastly, we evaluate the performance of the Delta rule variant introduced in Section 4 on the CIFAR-10 and CIFAR-100 benchmark. For these experiments, we instantiate the blocks of the ResNet-MLP backbone architecture, described in Appendix F.2, with vanilla dense linear layers. We train the parameters of the resulting network by applying the Delta weight update detailed in Section 4. Table 10: Hyperparameters Delta for experiments on CIFAR-10 and CIFAR-100. Hyperparameter Value Dataset CIFAR-10, CIFAR-100 Input 32×32×332× 32× 3 images flattened to 30723072 Model Residual MLP Residual stream size (ded_e) 30723072 # residual (L) [2,4][2,4] Normalization in blocks RMSNorm (pre-norm, per block) Final normalization RMSNorm before head MLP activation GELU Decay rate λ [0.0000003,0.000001,0.000003,0.00001,0.0001,0.0003][0.0000003,0.000001,0.000003,0.00001,0.0001,0.0003] Classifier head Linear (3072→10)(3072→ 10) Learning rate α [0.0001,0.0003,0.001,0.003,0.01][0.0001,0.0003,0.001,0.003,0.01] Batch size 6464 Epochs 5050 We report the hyperparameters used to produce the numbers for the Delta update rule in Figure 2 in Table 10. On top, for Figure 2/Right we also include the batch size in the sweeps and consider the values 4,16,324,16,32 and 6464. F.8.1 Additional Experimental Results: Extended training for deep networks with Delta regularization Figure 6: Training and test accuracy for a 10-layer residual MLP trained with SGD under varying DeltaNet decay rate λ. While λ=0λ=0 achieves near-perfect training accuracy quickly, test accuracy plateaus at around 67%. Modest DeltaNet-like decay rates significantly improve generalization, as indicated by the reduced train-test gap, which demonstrates the regularizing effect of the update rule. Curves are smoothed using a running average. Figure 7: Training and test accuracy for a 10-layer residual MLP trained with SGD under varying multiplicative decay rate β. Similarly to DeltaNet-like decay the training accuracy increases much slower, however, generalization improves only marginally. Curves are smoothed using a running average. Here, we present additional experimental results to evaluate the DeltaNet-like gradient descent update, which were obtained in a different experimental setting compared to the results in Section 5, to better assess the regularizing effect for deep neural networks. Specifically, we train a 10-layer ResNet MLP on Cifar-10 for 2000 epochs. Our architecture downprojects the preprocessed and flattened input to a hidden size of 512, and uses 2-layer MLPs with 4x expansion ratio, GELU activation [30] and a pre-block Layer norm [3]. We train with a batch size of 64 and implement the DeltaNet update in the batch setting by applying all decays within a batch simultaneously, that is, multiplying the linear layer by −λtt⊤ I-λ X_t X_t before applying the gradient update, where t∈ℝdin×B X_t ^d_in× B denotes the stacked input activations from batch t. In contrast to the derivation in Section 4, we reparametrize the update rule to align the decay rate with the learning rate scheduling (analogous to standard weight decay). Specifically, we update the weights as t=t−1(−ληttt⊤)−ηttt⊤, W_t= W_t-1( I-λ _t X_t X_t )- _t G_ Y_t X_t , (122) where t G_ Y_t contains the stacked output gradients, λ is the Delta rule decay rate and ηt _t is the learning rate, for which we use cosine annealing with linear warmup. A more detailed description of the experimental setup can be found in Table 11. We compare plain SGD, SGD with weight decay, and SGD with Delta rule decay when training for 20002000 epochs. Note that this is significantly longer than in other experiments (with a significantly deeper model), as it is computationally feasible and allows us to better understand the effect of regularization. As shown in Fig. 6, adding a DeltaNet-like decay mechanism significantly improves the final test accuracy (74% vs. 67%), far exceeding that of multiplicative weight decay (Fig. 7). While both mechanisms lead to a much slower increase in train accuracy, only the DeltaNet-like decay achieves superior generalization. We also tested combining both decay mechanisms (as in Gated DeltaNet), but this yielded no additional improvement. Additionally, we experimented with normalized DeltaNet updates as is standard for sequence modeling, i.e., using an input-dependent write strength of λη/‖22λη/\| x\|_2^2. However, this slightly degraded performance. A likely explanation is that faithfully following the DeltaNet derivation would require scaling the learning rate accordingly; doing so, however, leads to highly unstable optimization, since it makes the (effective) per-sample learning rate dependent on the activation norm. Table 11: Hyperparameters used for DeltaNet-GD experiments on extended training on Cifar-10 with deep nets. Hyperparameter Value Dataset CIFAR-10 Input 32×32×332× 32× 3 images flattened to 30723072 Model Residual MLP Input projection Linear (3072→512)(3072→ 512) Residual stream size (ded_e) 512512 # residual blocks 1010 Normalization in blocks LayerNorm (pre-norm, per block) Final normalization LayerNorm before head MLP expansion 4×4× MLP hidden dimension 20482048 MLP activation GELU Dropout 0.10.1 Classifier head Linear (512→10)(512→ 10) Optimizer SGD Learning rate α 1×10−21× 10^-2 Momentum 0.00.0 LR schedule Step-wise: 5 epoch linear warmup + cosine annealing Cosine schedule CosineAnnealingLR over remaining steps Batch size 6464 Epochs 20002000 F.9 MLP Baselines We complement the RCDL modules discussed in this work with strong MLP baselines on image classification tasks, using the classic CIFAR-10 and CIFAR-100 benchmarks. For these experiments, we instantiate the blocks of the ResNet-MLP backbone architecture, described in Appendix F.2, with vanilla dense linear layers. We train the parameters of the resulting network with either SGD [55], RMSProp [32], Adam [41], AdamW [45], and Muon [39]. Table 12: Hyperparameters MLP Baselines for experiments on CIFAR-10 and CIFAR-100. Hyperparameter Value Dataset CIFAR-10, CIFAR-100 Input 32×32×332× 32× 3 images flattened to 30723072 Model Residual MLP Residual stream size (ded_e) 30723072 # residual (L) [2,4][2,4] Normalization in blocks RMSNorm (pre-norm, per block) Final normalization RMSNorm before head MLP activation GELU Classifier head Linear (3072→10)(3072→ 10) Optimizer [SGD,RMSProp,Adam,AdamW,Muon][SGD,RMSProp,Adam,AdamW,Muon] Learning rate α [0.0001,0.0003,0.001,0.003,0.01][0.0001,0.0003,0.001,0.003,0.01] Weight decay (for AdamW) [0.0001,0.01][0.0001,0.01] Momentum 0.00.0 Batch size [16,64][16,64] Epochs 5050 We report the hyperparameters used to produce the numbers for the respective optimizer in Figure 2 in Table 12. On top, for Figure 2/Right we also include the batch size in the sweeps and consider the values 4,16,324,16,32 and 6464. Table 13: Hyperparameters MLP Baseline for Teacher-Student setup. Hyperparameter Value Model Residual MLP Residual stream size (ded_e) [256,512,1024,4096,8192][256,512,1024,4096,8192] # residual (L) 22 Normalization None MLP activation ReLU Optimizer SGD Learning rate [0.0003,0.001,0.003,0.01,0.03,0.1,0.3,1.0,3.0][0.0003,0.001,0.003,0.01,0.03,0.1,0.3,1.0,3.0] Momentum [0.0,0.9,0.999][0.0,0.9,0.999] Batch size [64,256,1024][64,256,1024] Training steps 200000200000 We report the hyperparameters used to produce the numbers for the reference MLP in the Teacher-Student setup in Figure 1 in Table 13. F.10 Nonparametric Baselines In Figure 2/Right we complement the RCDL methods with two nonparametric baselines: A kernel ridge regression (least squares (LS-)SVM) on the 10-dimensional one-hot targets and a k-Nearest Neighbour (kNN) implementation. F.10.1 Kernel ridge regression (least squares (LS-)SVM) For the LS-SVM we chose the RBF-kernel k(,′)=exp(−γ∥−′∥22)k( x, x )= (-γ x- x _2^2) inducing the Gram matrix ∈ℝn×n K ^n× n, with Ki,j=k(i,j)K_i,j=k( x_i, x_j), where i,j x_i, x_j are the channel-wise normalized RGB images, flattened into 30723072 dimensional vectors. We parametrize the bandwidth parameter following the median heuristic γ=γmultmediani,j∥i−j∥22γ= _multmedian_i,j x_i- x_j _2^2 and approximate the sample median on all pairwise distances of a subsample of 20482048 images. With the one-hot targets ∈ℝn×10 Y ^n× 10 we invoke a standard solver to compute =(+λn)−1∈ℝn×10 α=( K+λ I_n)^-1 Y ^n× 10 and evaluate f()=∑i=1nαik(,i)∈ℝ10f( x)= _i=1^n _ik( x, x_i) ^10 on the test samples x (predicting the argmax class). Table 14: Hyperparameters RBF kernel regression for experiments on CIFAR-10. Hyperparameter Value Dataset CIFAR-10 Input 32×32×332× 32× 3 images flattened to 30723072 Model RBF kernel Classification argmax on predicted output Optimizer Exact solver Bandwidth γmult _mult [0.25,0.5,1.0,2.0,3.0,5.0][0.25,0.5,1.0,2.0,3.0,5.0] Ridge regulariser λ [0.001,0.01,0.1,0.3,1.0][0.001,0.01,0.1,0.3,1.0] We report the considered hyperparameter ranges for producing the data-point displayed in Figure 2/Right in Table 14. F.10.2 k-Nearest Neighbors (kNN) We implement a kNN search using similarity score s(,i)s( x, x_i) on the channel-wise normalized and flattened test image and respective database entry ,i∈ℝ3072 x, x_i ^3072. We consider two similarity scores s(,i)=−∥−i∥22euclidean⊤iunnormalised dot product.s( x, x_i)= cases- x- x_i _2^2&euclidean\\ x x_i&unnormalised dot product. cases Denoting k()N_k( x) as the indices of the top-k neighbors (according to the similarity score) we then weight the labels of these entries either uniformly or proportional to the respective similarity score leading to one of the following weights wi=1uniform (regardless of similarity score)1∥−i∥2+10−5proportional, euclidean (inverse distance)exp(⊤i)∑j∈Nkexp(⊤j)proportional, dot (softmax over top-k),w_i= cases1&uniform (regardless of similarity score)\\ 1 x- x_i _2+10^-5&proportional, euclidean (inverse distance)\\ ( x x_i) _j∈ N_k ( x x_j)&proportional, dot (softmax over top-k), cases for i∈k()i _k( x). Using these weights, the prediction for an input x is given by taking the argmax over ∑i∈k()wi1[yi=c]∑i∈k()wi. _i _k( x)w_i1[y_i=c] _i _k( x)w_i. To obtain the number reported in Figure 2/Right we sweep the grid of the two similarity metrics, the two weighting mechanisms, and k∈5,10,15,20,25k∈\5,10,15,20,25\. F.11 Vision Transformer Experiments We investigate the performance of RCDL modules beyond ResNet architectures in classic vision Transformers (ViTs). We employ pre-LayerNorm transformer blocks consisting of standard softmax-attention and feed-forward networks (FFNs) and a linear readout on the CLS token. In the Weighted Softmax ViT, we replace each block’s FFN with a single Growing Weighted Softmax Module (Algorithm 3), all other components of the model are trained with standard parametric learning. Both are trained for a single epoch on CIFAR-10 (Appendix F.1.1), with hyperparameters in Tables 15 and 16. Table 15: Hyperparameters Standard ViT for single-epoch experiments on CIFAR-10. Hyperparameter Value Dataset CIFAR-10 Input 32×32×332× 32× 3 images, 4×44× 4 patches (6464 patch tokens + CLS token) Model ViT, 22 transformer blocks Embedding dimension ded_e 128128 Attention heads 88 Positional embedding learned Normalization in blocks LayerNorm (pre-norm, per block) + final LayerNorm FFN Dense (128→512128→ 512) → GELU → Dense (512→128512→ 128) FFN blowup factor [1,2,4,8,16,32][1,2,4,8,16,32] Readout Linear layer on CLS token Optimizer Adam Learning rate α [0.0003,0.001,0.003][0.0003,0.001,0.003] Batch size 6464 Epochs 11 Seeds 55 Table 16: Hyperparameters Weighted Softmax ViT for single-epoch experiments on CIFAR-10. Hyperparameter Value Dataset CIFAR-10 Input 32×32×332× 32× 3 images, 4×44× 4 patches (6464 patch tokens + CLS token) Model ViT, 22 transformer blocks, FFN replaced by growing GWS Embedding dimension ded_e 128128 Attention heads 88 Positional embedding learned Normalization in blocks LayerNorm (pre-norm, per block) + final LayerNorm Readout Linear layer on CLS token Optimizer (non-growing params.) Adam Learning rate (Adam) α [0.0003,0.001,0.003][0.0003,0.001,0.003] Weight learning rate η [0.1,0.3,1.0][0.1,0.3,1.0] Value learning rate β [10,30][10,30] Softmax temperature T [0.1,0.3,1.0,3.0][0.1,0.3,1.0,3.0] Value momentum [0.0,0.999][0.0,0.999] Batch size 6464 Epochs 11 Seeds 55 F.12 Computational considerations Growing Weighted Softmax, RBF-kernel based RCDL, and Parametric Weighted Softmax modules have a significant memory footprint. In the reported CIFAR-10,100 experiments, where each epoch consists of 5×1045× 10^4 training examples, we add 2.5×1062.5× 10^6 tuples (w,,)(w, v, k) (or (,)( v, k), respectively, for RBF-kernel based modules) to D in growing modules. For a hidden layer with din=dout=3072d_in=d_out=3072, implemented in float32float32, as used in CIFAR experiments, this results in a total memory footprint of 61.561.5GB (sum of memory footprint for both keys, values and weights). To realize the training we rely on parallelization and shard this database across 44 TPUs for 5050-epoch CIFAR experiments. For the teacher-student regression experiments, training for 200200 epochs adds 10710^7 tuples to D in growing modules. For a hidden layer with din=dout=100d_in=d_out=100, this results in a minimal memory footprint of approximately 88GB. As a result, this task does not require parallelization or database sharding and can be efficiently trained on a single device. However, note that the update is relatively cheap in growing as we do not have to compute gradients for these large databases. This is in contrast to Parametric Weighted Softmax modules, where we furthermore have to consider the additional footprint of the equally sized gradients, and if using advanced optimizers such as Adam [41], first and second moments, making training significantly more expensive and slower than by learning the Weighted Softmax module by growing. Appendix G Further Experimental Results G.1 Training loss on CIFAR-100 Complementary to Figure 2/Middle where we report the train loss on CIFAR-10, we here report the respective numbers for CIFAR-100 in Figure 8. Figure 8: Training loss for Weighted Softmax & RBF decreases much faster than that of the other models. G.2 Gauge-Invariant Weighted Softmax We compare the Gauge-Invariant Weighted Softmax (Algorithm 2) with the standard Weighted Softmax and RBF in the single-epoch setting of Fig. 2/Right: non-augmented CIFAR-10, batch size 1616, L=2L=2 residual blocks, de=3072d_e=3072, and RMSNorm without learned scale. For Weighted Softmax and RBF we use the best configurations from the respective sweeps: Weighted Softmax with T=6T=6, η=0.003η=0.003, β=3β=3; RBF with temperature 2020 and learning rate 1010. For the Gauge-Invariant variant we use T=5T=5, α=10α=10, η=0.03η=0.03, β=100β=100. Fig. 9 shows mean and standard deviation over seeds. All three methods reach similar final test accuracy (Weighted Softmax 59.4±0.4%59.4± 0.4\%, RBF 58.8±0.2%58.8± 0.2\%, Gauge-Invariant 58.4±0.4%58.4± 0.4\%). Guaranteeing non-negative weights by the choice of gauge, rather than by clipping, therefore performs comparably and removes the need for the clipping heuristic. Given this similar behavior, we decided for the results in the main text of the paper to only analyze standard weighted softmax RCDL layers. Figure 9: Test accuracy of non-parametric methods (Weighted Softmax, Gauge-Invariant WSM & RBF) as a function of the number of samples seen during training. Appendix H Further Related Work and Discussions Here we provide further related work and discussions. At a high-level, retrieval-centric deep learning (RCDL) also shares conceptual similarities with existing “differentiable memory models” such as Differentiable Neural Computer [27, 28], Memory Networks [74, 66] and Neural Episodic Control [53] (see also more recent work in the context of language models, such as Berges et al. [10], Eyuboglu et al. [20]); with a crucial difference that RCDL’s memory retrieval operates on the training datapoints. Its growing behavior, recruiting “resources” (new key/value allocations) as the system processes increasingly more training data (conceptually to infinity, with practical hardware limitations) is somewhat also reminiscent of Turing machines [68] as is the case for the transformer at runtime [69] (see also how Behrouz et al. [5] adopt a similar principle to recurrent nets). RCDL also differs from “deep kernel learning” [76] which uses a deep N as a learnable pre-processor to produce features for a single-layer Gaussian Process (GP); as well as from a special variant of modern Hopfield-lookup layers [54] that stores the static training data and retrieves their labels via softmax attention. This comparison also provides an interesting contrast to infinite-width deep learning theories, where NNs asymptotically reduce to a kernel machine, such as the neural tangent kernel (NTK; Jacot et al. [38], Lee et al. [43]); in that regime, a central question is whether deep NNs meaningfully learn representations or operate in a lazy training regime with fixed random features [16]. In contrast, RCDL learns representations in an N that is explicitly structured as a (deep) kernel machine—an intriguing future work may be to study RCDL through the lens of NTK/GP. The challenge of expanding model capacity has antecedents in early constructive algorithms [21, 52], which dynamically allocate new units to reduce error. In modern deep learning, this concept evolved into architectures like “progressive NNs” [57], which instantiate new columns for novel tasks to prevent forgetting. In fact, the idea of growing networks is rather common in continual learning (see, e.g., Yoon et al. [80]), as well as in architecture search/optimization (see, e.g., Stanley and Miikkulainen [65], Evci et al. [19], Butkus et al. [15]). RCDL offers a distinct path: unlike heuristics that add coarse-grained parametric modules, RCDL grows by accumulating fine-grained nonparametric key-value pairs, resulting in capacity that scales linearly with training data. Crucially, while this growth rate is fixed, the update rule itself is strictly principled: we do not add units arbitrarily, but rather derive the addition of each new key-value pair as a precise step of functional gradient descent on the global loss. “Retrieval” has also received recent popularity in the context of large language models (LLMs) [44, 49] as a method to complement the deficiency of LLMs’ knowledge stored in their fixed-size weight matrices. In a sense, RCDL integrates such a process within the network itself; store the entire internet in the N itself, and “compression” occurs at inference time through attention that recombines features—a slogan here is “store the entire internet in a growing neural net, compression can wait!”.