# Training Neural Networks at Any Scale

Thomas Pethick, Kimon Antonakopoulos, Antonio Silveti-Falls, Leena

Chennuru Vankadara, Volkan Cevher

<figure><img src="image_1.png" /></figure>

## To cite this version:

Thomas Pethick, Kimon Antonakopoulos, Antonio Silveti-Falls, Leena Chennuru Vankadara, Volkan Cevher. Training Neural Networks at Any Scale. 2025. hal-05420726

HAL Id: hal-05420726

https://inria.hal.science/hal-05420726v1

Preprint submitted on 17 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Training Neural Networks at Any Scale

Thomas Pethick*, Kimon Antonakopoulos*, Antonio Silveti-Falls†,
Leena Chennuru Vankadara‡, Volkan Cevher*

## INTRODUCTION

Deep learning thrives on the powerful interplay between data, architectures, and optimization. As the sizes of foundation models continuously grow in scale with 405B Llama models and 671B DeepSeek models, we see unprecedented performance and new transformative capabilities. However, this increase in scale also demands optimization methods that scale efficiently and reliably.

Traditionally, optimization theory for neural networks has treated training as a black-box optimization problem. By black-box optimization, we mean optimizing without taking into account possible structure present in the objective function or its gradients. While this generality has been convenient, it has limited our ability to exploit inherent characteristics of the problem that can significantly enhance performance. Recent advances [1]–[3], however, have demonstrated the power of jointly designing neural network architectures, initialization schemes, and optimization algorithms, while understanding their training regimes through the notion of feature learning.

Key insights have emerged from the study of infinite-width neural network limits, particularly through frameworks such as Neural Tangent Kernel (NTK) [4] and mean-field analyses [5]. A landmark development, called the Maximal Update Parameterization (μP) [6], showed that feature learning can be ensured when the network architecture, initialization, and layerwise stepsizes are carefully aligned. Remarkably, this alignment also results in a beneficial property called hyperparameter transfer, in which hyperparameter settings (like learning rate) remain optimal across model scales, significantly stabilizing the training.

Furthermore, optimization methods have recently evolved to be both more reliable and faster by explicitly leveraging the network's structure as a collection of layers. A notable example is the Shampoo optimizer [7], which topped the recent AlgoPerf benchmark hosted by MLCommons. Shampoo explicitly utilizes the network's structure as a collection of weight matrices by preconditioning the update to each layer's weight matrix separately, connecting elegantly to spectral descent methods dating back to

<footnote>*École Polytechnique Fédérale de Lausanne (EPFL), Laboratory for Information and Inference Systems (LIONS)</footnote>

<footnote>†CVN, CentraleSupélec, Université Paris-Saclay, Inria</footnote>

<footnote>‡Gatsby Computational Neuroscience Unit, University College London</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Different optimizers for training large language models (LLMs), illustrated with a GPT model having up to 3 billion parameters (plots from [1]). A modern optimizer like Scion (middle) exploits the geometry of the neural network to achieve faster convergence (right), while simultaneously being parameterized such that the optimizer becomes agnostic to scale (middle, all other parameters constant). This design leads to the transfer of hyperparameters, like learning rate, across model width, as demonstrated above. Performance is also boosted in large-scale training (3 billion parameters, right).</figcaption></figure>

[8]. As we explore in this article, the spectral geometry emerges naturally from analyzing how signals or features propagate through the network, and can be used for more effective training.

To this end, we present a contemporary perspective on deep learning optimization, highlighting key connections between optimization techniques and neural network architectures. Our goal is to equip both experts and practitioners with actionable insights on choosing and configuring optimizers in tandem with architectures to achieve scale-agnostic neural network training.

*Outline:* We begin by introducing a general framework for deriving optimization algorithms that leverage problem geometry through an appropriate choice of norm. Perhaps unexpectedly, we find that the classical stochastic Frank-Wolfe method for constrained optimization [9] emerges as a central component of this framework. Second, we examine scaling rules: *how to optimally adjust the training recipe as the problem size varies*. We analyze signal propagation in neural networks under scaling, focusing on how to balance training stability with meaningful feature evolution, particularly as network width increases.

Neural networks can learn features even in the infinite-width limit when a specific *parameterization*, i.e., a layerwise scaling of initialization variances and learning rates, is employed. Under common parameterizations, such as the Neural Tangent Parameterization (NTP) [4], infinitely wide networks operate in the so-called lazy regime and do not exhibit feature learning. In contrast, when the network satisfies the maximal update criterion [6], which promotes maximal feature learning across all layers, we empirically observe hyperparameter transfer (cf., Figure 1). In practice, this means that the optimal learning rate becomes independent of model width and a learning rate tuned on a smaller model can be directly applied to a larger one. We conclude with a discussion of current limitations and open challenges.
GEOMETRY-AWARE OPTIMIZATION METHODS

Adam, an adaptive optimizer based on prior methods like AdaGrad and RMSProp, has dominated deep learning by dynamically changing the geometry in which the updates to parameters are computed. Methods like this adapt to the local geometry on-the-fly using the gradients seen by the algorithm so far. This enables reliable training with little tuning. However, these algorithms also treat the neural network (NN) to be trained as a black-box with all parameters treated irrespective of their role or position in the architecture.

A different approach instead chooses an appropriate geometry to compute the updates to parameters prior to the optimization phase. This choice of geometry is encoded through a norm on the parameters and is deliberately made with the NN architecture to be trained in mind, in contrast to the previously described adaptive methods. In the following section, we lay out both of these algorithmic approaches under a unified algorithmic template involving steepest descent and conditional gradient methods.

Preliminaries

We begin by describing the mathematical formalism that we will use to discuss different algorithms and initializations. We frame training as solving the following non-convex optimization problem:

$$ \min_{x \in \mathcal{X}} f(x), \qquad (1) $$

where $f$ is $L$-smooth (i.e., its gradient $\nabla f$ is Lipschitz continuous with constant $L > 0$ in some general norm) and the domain $\mathcal{X}$ is unconstrained and given by $\mathbb{R}^p$, where $p$ is the number of NN parameters.

We assume access to some stochastic oracle $f(\cdot, \xi)$ instead of $f(\cdot)$ directly, where $\xi$ is a random variable, typically representing a minibatch of data. From an algorithmic perspective, we aim to solve (1) with iterative methods that use a stochastic first-order oracle (SFO). At each iteration, the algorithm can query a black-box mechanism at a given point $x$ that returns a stochastic estimate $g(x, \xi)$ of $\nabla f(x)$.

For instance, the workhorse of modern deep learning, AdamW [10], presented below without the bias correction proceeds with the following updates at each iteration $k$, given $\beta_{1,2} \in [0, 1)$:

$$ d^k = \beta_1 d^{k-1} + (1 - \beta_1)g^k, \quad v^k = \beta_2 v^{k-1} + (1 - \beta_2)(g^k \odot g^k), \quad x^{k+1} = (1 - \underbrace{s_k \lambda}_{=\lambda_k})x^k - \underbrace{s_k \gamma}_{=\gamma_k} \frac{d^k}{\sqrt{v^k + \varepsilon^2}}, $$

where multiplication $\odot$, division and square root are taken coordinatewise; $\lambda \in [0, 1)$ is the decoupled weight decay constant; $\gamma$ is the stepsize; $s_k$ is the learning rate scheduler; and $\varepsilon > 0$ is a small numerical constant. We take $d^0 = g^0$ and $v^0 = g^0 \odot g^0$ to avoid the bias correction. In the sequel, we refer to $\lambda_k$ and $\gamma_k$ as the scheduled weight decay and scheduled stepsize respectively.

Intuitively, AdamW builds a diagonal preconditioner on-the-fly that rescales coordinates by an exponential moving average (EMA) of their gradient energy via $\beta_2$. As opposed to directly using the stochastic gradients $g^k$, AdamW uses their (Polyak) momentum accumulation via $\beta_1$. AdamW relies on the weight
decay parameter $\lambda > 0$ to prevent the parameters from growing unbounded, which is observed to improve generalization.

To connect AdamW to a broader picture, we introduce the concept of “dual feedback,” denoted by $d^k$, which represents a key signal the optimizer uses at each iteration. In the deterministic setting, $d^k$ can simply be the gradient $\nabla f(x^k)$ while in the stochastic setting choosing $d^k$ can be more sophisticated, involving things like variance reduction. For now we will consider two main examples, in which $d^k$ is taken to be $g^k$ directly or a momentum accumulation of previous $g^k$; both can be defined recursively via $d^k = (1 - \alpha_k)d^{k-1} + \alpha_k g^k$ for some $\alpha_k \in (0, 1]$, which is equivalent to the estimator in Adam when $\alpha_k = 1 - \beta_1$. This momentum estimator reduces variance but is a biased estimate of $\nabla f(x^k)$.

Having established the dual feedback, we now introduce our general master algorithmic template:

$$x^{k+1} \in \arg \min_{x \in \mathcal{X}} \gamma_k \langle d^k, x \rangle + h_k(x), \quad (2)$$

where $\gamma_k > 0$ is the stepsize and $h_k$ is a convex function which facilitates a tractable solution of (2), i.e., a regularizer. The generality of this formulation allows us to unify several celebrated algorithms for training neural networks within the somewhat classical framework below with different choices of $d^k$ and $h_k$.

In what follows, we explore the flexibility of the scheme in (2) to adapt to the geometry of (1), both on-the-fly or a priori, through appropriate choice of the regularizer $h_k$. We will also discuss how the two other parts, the dual feedback $d^k$ and the stepsize $\gamma_k$, are essential for treating the stochastic case.

## On-the-fly adaptation to the geometry

We start with describing how AdamW adapts on-the-fly to the problem geometry. To make this explicit, we first consider the norm associated to the Mahalanobis inner product in a Euclidean space as follows: $\|x\|_{2,H} = \sqrt{\langle x, H^{1/2}x \rangle}$, where $H > 0$ is a positive define matrix. Then, given a sequence of such matrices $\{H_k\}$, we can define a sequence of regularizers that will play a unifying role for adaptation:

$$h_k(x) = \frac{1}{2} \|x - (1 - \lambda_k)x^k\|_{2, H_k}^2. \quad (3)$$

Using this family of regularizers, (2) reduces to the preconditioned stochastic gradient descent method with decoupled weight decay PreSGDW, given by the following generic template:

$$x^{k+1} = (1 - \lambda_k)x^k - \gamma_k H_k^{-1/2} d^k. \quad (\text{PreSGDW})$$

Alternatively, choosing $d^k = g^k$, $H_k = I$ (identity), and weight decay parameter $\lambda = 0$, we reach the stochastic gradient descent (SGD). By choosing $d^k$ as the momentum estimator with $\alpha_k = 1 - \beta_1$ and $H_k = \text{diag}(v^k + \varepsilon^2 I)$, we reach AdamW.
<table><caption>TABLE I
A CHRONOLOGICAL TAXONOMY OF SELECTED FIRST-ORDER OPTIMIZATION METHODS UNDER PRESGDW</caption><thead><tr><th>Algorithm</th><th>Reference</th><th>Dual feedback d<sup>k</sup></th><th>Preconditioner H<sub>k</sub><sup>1</sup></th></tr></thead><tbody><tr><td>SGD</td><td>[11]</td><td>g<sup>k</sup></td><td>[H<sub>k</sub>]<sub>ii</sub> = 1; ∀i ∈ [p]</td></tr><tr><td>Momentum (Polyak)</td><td>[11]</td><td>β<sub>1</sub>d<sup>k-1</sup> + (1 - β<sub>1</sub>)g<sup>k</sup></td><td>[H<sub>k</sub>]<sub>ii</sub> = 1; ∀i ∈ [p]</td></tr><tr><td>AdaGrad (diagonal)</td><td>[12], [13]</td><td>g<sup>k</sup></td><td>[H<sub>k</sub>]<sub>ii</sub> = [H<sub>k-1</sub>]<sub>ii</sub> + (g<sub>i</sub><sup>k</sup>)<sup>2</sup>; ∀i ∈ [p]</td></tr><tr><td>AdaGrad (full)</td><td>[12], [13]</td><td>g<sup>k</sup></td><td>H<sub>k</sub> = H<sub>k-1</sub> + g<sub>i</sub><sup>k</sup>(g<sub>i</sub><sup>k</sup>)<sup>T</sup></td></tr><tr><td>RMSProp</td><td>[14]</td><td>g<sup>k</sup></td><td>[H<sub>k</sub>]<sub>ii</sub> = β<sub>2</sub>[H<sub>k-1</sub>]<sub>ii</sub> + (1 - β<sub>2</sub>)(g<sub>i</sub><sup>k</sup>)<sup>2</sup>; ∀i ∈ [p]</td></tr><tr><td>Adam / AdamW</td><td>[10], [15]</td><td>β<sub>1</sub>d<sup>k-1</sup> + (1 - β<sub>1</sub>)g<sup>k</sup></td><td>[H<sub>k</sub>]<sub>ii</sub> = β<sub>2</sub>[H<sub>k-1</sub>]<sub>ii</sub> + (1 - β<sub>2</sub>)(g<sub>i</sub><sup>k</sup>)<sup>2</sup>; ∀i ∈ [p]</td></tr><tr><td>Stochastic ℓ<sub>∞</sub>-descent<sup>2</sup></td><td>[8], [16]</td><td>g<sup>k</sup></td><td>[H<sub>k</sub>]<sub>ii</sub> = (g<sub>i</sub><sup>k</sup>)<sup>2</sup>; ∀i ∈ [p]</td></tr><tr><td>Stochastic spectral descent<sup>3</sup></td><td>[8], [16]</td><td>G<sub>k</sub> = mat(g<sup>k</sup>)</td><td>R<sub>k</sub> ⊗ I; I ⊗ L<sub>k</sub>; or R<sub>k</sub><sup>1/2</sup> ⊗ L<sub>k</sub><sup>1/2</sup>: <span style="font-size: 2em; vertical-align: middle;">{</span> R<sub>k</sub> = G<sub>k</sub>G<sub>k</sub><sup>T</sup><br>L<sub>k</sub> = G<sub>k</sub><sup>T</sup>G<sub>k</sub><span style="font-size: 2em; vertical-align: middle;">}</span></td></tr><tr><td>SignSGD / Signum</td><td>[17]</td><td>g<sup>k</sup>/β<sub>1</sub>d<sup>k-1</sup> + (1 - β<sub>1</sub>)g<sup>k</sup></td><td>[H<sub>k</sub>]<sub>ii</sub> = (g<sub>i</sub><sup>k</sup>)<sup>2</sup>; ∀i ∈ [p]</td></tr><tr><td>Shampoo</td><td>[7]</td><td>G<sub>k</sub> = mat(g<sup>k</sup>)</td><td>R<sub>k</sub><sup>1/2</sup> ⊗ L<sub>k</sub><sup>1/2</sup>: <span style="font-size: 2em; vertical-align: middle;">{</span> R<sub>k</sub> = β<sub>2</sub>R<sub>k-1</sub> + (1 - β<sub>2</sub>)G<sub>k</sub>G<sub>k</sub><sup>T</sup><br>L<sub>k</sub> = β<sub>2</sub>L<sub>k-1</sub> + (1 - β<sub>2</sub>)G<sub>k</sub><sup>T</sup>G<sub>k</sub><span style="font-size: 2em; vertical-align: middle;">}</span></td></tr><tr><td>uScion / Muon (spectral)</td><td>[1], [2]</td><td>β<sub>1</sub>d<sup>k-1</sup> + (1 - β<sub>1</sub>)G<sub>k</sub></td><td>R<sub>k</sub> ⊗ I; I ⊗ L<sub>k</sub>; or R<sub>k</sub><sup>1/2</sup> ⊗ L<sub>k</sub><sup>1/2</sup>: <span style="font-size: 2em; vertical-align: middle;">{</span> R<sub>k</sub> = G<sub>k</sub>G<sub>k</sub><sup>T</sup><br>L<sub>k</sub> = G<sub>k</sub><sup>T</sup>G<sub>k</sub><span style="font-size: 2em; vertical-align: middle;">}</span></td></tr></tbody></table>

Table I summarizes several popular first-order optimization algorithms as specific instantiations of PreSGDW. Each method differs only in two essential modeling choices: (i) how the dual feedback $d^k$ is formed from the observed gradients, and (ii) what geometry, encoded via the preconditioner $H_k$, and hence the regularizer $h_k$ above, is used to compute the next iterate. These differences reveal how existing optimizers are related in the way they adapt to the geometry of the underlying problem.

Methods like SGD and Polyak momentum operate under a fixed Euclidean geometry with an identity preconditioner. In contrast, AdaGrad, RMSProp, and AdamW introduce a diagonal preconditioner that evolves with the observed gradients, endowing the parameter space with an anisotropic geometry that varies across iterations. From the table, it is clear that RMSProp is a generalization of diagonal AdaGrad with a tunable $\beta_2$ parameter in the computation of $H_k$, while Stochastic $\ell_\infty$-descent is the limit of RMSProp with $\beta_2 \to 0$. AdamW further expands this class of algorithms with momentum-based dual feedback, reducing dependence on (unbiased) minibatch stochastic gradients with $\beta_1$.

In general, going beyond diagonal preconditioners in PreSGDW can dramatically improve performance but this is especially resource heavy in deep learning, including memory and compute. However, we can

<footnote n="1">¹ Unless stated otherwise, we have $[H_k]_{ij} = 0; \forall i \neq j \in [p]$; ² The stepsize $\gamma_k$ needs to be scaled by $||d^k||_1$.</footnote>

<footnote n="3">³ The stepsize $\gamma_k$ needs to be scaled by $||d^k||_{S_1}$.</footnote>
obtain a practical compromise by maintaining {$H_k = R_k \otimes L_k$} with kronecker-factors {$R_k$} and {$L_k$} for right and left preconditioning, respectively. We can then treat the stochastic first-order dual feedback $d^k$ layerwise as a matrix $D_k = \mathrm{mat}(d^k)$ and rely on the identity $H_k^{-1/2}d^k = \mathrm{vec}(L_k^{-1/2}D_kR_k^{-1/2})$. Table I describes Stochastic spectral descent (SSD) and Shampoo as two algorithms that use this perspective. We see that SSD is the limit of Shampoo, with both forming the spectral analogs of stochastic $\ell_\infty$-descent and AdaGrad, respectively.

It is worth noting that the inclusion of the weight decay term $\lambda_k$ inside the regularizer (3) is ad-hoc; we write things this way just so that all methods fit the same mold in Table I. Indeed, it is introduced primarily to obtain the generic form PreSGDW that unifies preconditioning and decoupled weight decay within a single framework. In the sequel, we shift focus to a priori geometry design, where we will encode architectural knowledge directly into the optimizer through a norm that does not change across iterations of the algorithm. Also within this new framework, the weight decay $\lambda$ will play a more natural and prominent role.

## A priori adaptation to the geometry

When training NNs, the form of (1) has a particular, predictable structure induced by the network architecture and the input domain. This useful yet relatively generic information can be encoded geometrically via the construction of a norm on $\mathcal{X}$. The norm $\| \cdot \|$ together with its dual norm,

$$ \|u\|_* := \max_{\|x\| \le 1} \langle u, x \rangle, \qquad (\text{Dual Norm}) $$

can guide the choice of $h_k$ to get a priori adaptation of the optimizer to the new geometry. We will describe how to construct such a norm without an inner product structure in later sections; for now we start by introducing two prominent families of regularizers based on this non-Euclidean paradigm.

**Steepest descent:** One natural way to generalize SGD to a new geometry is to simply regularize with a non-Euclidean norm $h_k(x) = \frac{1}{2}\|x - x^k\|^2$. With this choice, the algorithmic template (2) reduces to the stochastic steepest descent (SD) method in a normed space [8], [16], which can be expressed concisely by defining the sharp-operator, $s^\sharp \in \arg\max_{x \in \mathcal{X}} \{\langle s, x \rangle - \frac{1}{2}\|x\|^2\}$, introduced in [18]:

$$ x^{k+1} = x^k - \gamma_k [d^k]^\sharp. \qquad (\text{Stochastic SD}) $$

Here, the sharp-operator maps the dual feedback $d^k$ to the primal space according to the chosen geometry.

Key for deep learning applications are geometries that induce so-called dense updates. An illustrative example is given by the $\ell_\infty$-vector norm, for which Stochastic SD reduces to a sign-based update, e.g., Stochastic $\ell_\infty$-descent (cf., Table I) [8]:

$$ x^{k+1} = x^k - \gamma_k \|d^k\|_1 \mathrm{sign}(d^k). \qquad (4) $$
The update is dense as each nonzero coordinate of $d^k$ produces the same magnitude change in $x^{k+1}$.

The matrix analog of the $\ell_\infty$-norm is the Schatten-$\infty$ matrix norm (also known as the spectral norm) which, if chosen to induce $h_k$, results in the stochastic spectral descent method [8]:

$$x^{k+1} = x^k - \gamma_k \|\sigma^k\|_1 U^k (V^k)^\top,$$

where the reduced singular value decomposition (SVD) is given as $d^k = U^k \text{diag}(\sigma^k)(V^k)^\top$.

For geometries arising from an inner product, which can represent both the updates of SGD and PreSGD, the sharp-operator is linear. In contrast, the sharp-operator induced by a non-Euclidean geometry is notably nonlinear as explicitly illustrated in Table I though the master template (2). Nonetheless, preconditioning ideas can be naturally extended to this non-linear setting, as discussed in [19].

We also note that it is possible to construct a non-Euclidean norm tailored to the task of training deep neural networks by bounding the weight matrices layerwise in the spectral norm [16]. We will later revisit and motivate this idea in the context of feature learning.

**Conditional Gradient Methods:** While deterministic steepest descent is guaranteed to decrease the function value $f(x^k)$ at every iteration, it comes at the cost of picking a specific stepsize which is a function of the Lipschitz modulus $L$ of the gradient (specifically $\gamma_k < 2/L$). In practice, this leads to slower training, since typically the global Lipschitz modulus is too pessimistic where the function is less volatile.

We can address this dependence on the Lipschitz constant by normalizing the gradient's magnitude and considers only its direction. An example is the Normalized SGD (NSGD) with decoupled weight decay:

$$x^{k+1} = (1 - \underbrace{\lambda_k}_{=s_k \lambda}) x^k - \underbrace{\gamma_k \frac{\nabla f(x^k)}{\|\nabla f(x^k)\|_2}}_{=s_k \gamma}, \quad (\text{NSGD})$$

which uses the $\ell_2$-norm normalization. A very related technique in practice is known as *gradient norm clipping*, which only normalizes the gradient when its $\ell_2$-norm exceeds a predefined threshold.

A natural generalization of NSGD to non-Euclidean norms follows from noticing that the normalization can be written in terms of a Linear Minimization Oracle (lmo)

$$\text{lmo}_\rho(d) \in \arg \min_{x \in \mathcal{D}} \langle d, x \rangle, \quad (\text{lmo})$$

where the constraint is a (possibly non-Euclidean) norm-ball $\mathcal{D} := \{x \mid \|x\| \le \rho\}$. When $\rho = 1$, we will omit it in the notation. Indeed, by choosing the $\ell_2$-norm ball with radius $\rho = \gamma/\lambda$ for $\mathcal{D}$, NSGD can be seen as an instance of the so-called Stochastic Conditional Gradient method (SCG) [1], [9], [20]–[22],

$$x^{k+1} = (1 - \lambda_k)x^k + \lambda_k \text{lmo}_\rho(d^k) = (1 - \lambda_k)x^k + \lambda_k \rho \text{lmo}(d^k), \quad (\text{SCG})$$

which requires $\lambda_k \in [0, 1]$ and $x^0 \in \mathcal{D}$ to ensure the iterates satisfy constraint $x^k \in \mathcal{D}$. We refer to these class of methods as conditional gradient methods since they are *conditioned* on the constraint. This
algorithm is also known as the stochastic Frank-Wolfe algorithm, where $\lambda_k$ is referred to as Frank-Wolfe stepsize. SCG can be used to provide numerical solutions to the constrained problem (1) with $\mathcal{X} \equiv \mathcal{D}$.

**A twist in perspective:** In addition to NSGD, the algorithms Lion, Signum, Muon, LARS and LAMB with decoupled weight decay can all be written as instances of the stochastic conditional gradient method (SCG) as shown in [1], [23]. Intriguingly, SCG does not have a “weight decay” parameter. However, in light of their Frank-Wolfe interpretation, we can interpret that the roles of the stepsize and weight decay in these algorithms become inverted. The update direction corresponds to the output of the linear minimization oracle lmo over a ball of radius $\rho = \gamma/\lambda$ which does not depend on the dual feedback’s magnitude. Consequently, the coefficient $\gamma$ in front of the normalized gradient no longer determines the actual progress along the optimization path as it merely sets the implicit radius of the feasible set. Instead, the contraction factor $(1 - \lambda_k)$, commonly seen in the context of weight decay, *and is interpreted as such*, plays the role of the *Frank-Wolfe stepsize*, dictating how far the iterate moves toward the boundary of the norm ball at each step. Viewed through this lens, the so-called, decoupled “weight decay” parameter $\lambda$ is, geometrically, the true stepsize of the algorithm, whereas the nominal stepsize $\gamma$ only defines the scale of the constraint by scaling the lmo.

**Unconstrained variant:** However, normalized gradient methods without weight decay also exist. Indeed, NSGD without weight decay can be seen as an instance of the unconstrained Stochastic Conditional Gradient method (uSCG) [1], which is also known as the continuous greedy method in submodular optimization:

$$x^{k+1} = x^k + \gamma_k \text{lmo}(d^k). \qquad (\text{uSCG})$$

The uSCG method can also be seen as a trust-region method with a linear surrogate objective or as the generalized Matching Pursuit algorithm with a set of atoms $\mathcal{A}$ whose convex hull is $\mathcal{D}$, since the main operation used in those algorithms corresponds to the lmo. Whereas steepest descent is an instance of (2) with a regularizer $h_k$ that enforces a soft constraint, uSCG on the other hand is obtained by enforcing a hard constraint through the regularizer choice $h_k(x) = \iota_{\mathcal{D}}(x - x^k)$ where $\iota_{\mathcal{D}}(x) := 0$ if $x \in \mathcal{D}$ and $+\infty$ otherwise. This hard constraint explains the stability of the method, since the update never leaves the trust-region around the previous iterate $x^k$.

<figure><img src="image_4.png" /><figcaption>Fig. 2. Visualization of the lmo for several different balls corresponding to $\ell_p$ norms for $p \in \{1, \frac{3}{2}, 2, 3, +\infty\}$.</figcaption></figure>
<caption>TABLE II</caption>

SPECIAL INSTANTIATIONS OF CONDITIONAL GRADIENT METHODS (USCG AND SCG) AND STOCHASTIC SD ACCORDING TO DIFFERENT CHOICES OF NORM. THE REDUCED SVD IS GIVEN AS $d = U \text{diag}(\sigma)V^\top$ (HERE $d \in \mathbb{R}^{p_2 \times p_1}$ IS A MATRIX).

<table><thead><tr><th></th><th>Method</th><th>α<sub>k</sub></th><th>Norm</th><th>-lmo(d) or [d]<sup>#</sup></th><th>Reference</th></tr></thead><tbody><tr><td rowspan="4">Conditional Gradient</td><td>Normalized SGD</td><td>1</td><td>Euclidean ||·||<sub>2</sub></td><td><sup>d</sup>/<sub>||d||<sub>2</sub></sub></td><td>[24]</td></tr><tr><td>Normalized SGD with momentum</td><td>[0, 1]</td><td>Euclidean ||·||<sub>2</sub></td><td><sup>d</sup>/<sub>||d||<sub>2</sub></sub></td><td>[25]</td></tr><tr><td>SignSGD</td><td>1</td><td>Max-norm ||·||<sub>∞</sub></td><td>sign(d)</td><td>[17]</td></tr><tr><td>Signum</td><td>[0, 1]</td><td>Max-norm ||·||<sub>∞</sub></td><td>sign(d)</td><td>[17]</td></tr><tr><td>Scion, Muon<sup>1</sup></td><td>[0, 1]</td><td>Layerwise Spectral ||·||<sub>S<sub>∞</sub></sub></td><td>UV<sup>⊤</sup></td><td>[1], [2]</td></tr><tr><td rowspan="10">Steepest descent</td><td>SGD</td><td>1</td><td>Euclidean ||·||<sub>2</sub></td><td>d</td><td>[11]</td></tr><tr><td>SGD with momentum</td><td>[0, 1]</td><td>Euclidean ||·||<sub>2</sub></td><td>d</td><td></td></tr><tr><td>PreSGDW (AdaGrad, RMSProp)</td><td>1</td><td>Mahalanobis ||·||<sub>2,H</sub></td><td>(H)<sup>-1/2</sup>d</td><td>[12]–[14]</td></tr><tr><td>PreSGDW with momentum (Adam/AdamW)</td><td>[0, 1]</td><td>Mahalanobis ||·||<sub>2,H</sub></td><td>(H)<sup>-1/2</sup>d</td><td>[10], [15]</td></tr><tr><td>Norm-scaled SignSGD</td><td>1</td><td>Max-norm ||·||<sub>∞</sub></td><td>||d||<sub>1</sub> sign(d)</td><td>[26]</td></tr><tr><td>Norm-scaled SignSGD with momentum</td><td>[0, 1]</td><td>Max-norm ||·||<sub>∞</sub></td><td>||d||<sub>1</sub> sign(d)</td><td></td></tr><tr><td>Stoch. Spectral descent</td><td>1</td><td>Layerwise Spectral ||·||<sub>S<sub>∞</sub></sub></td><td>||σ||<sub>1</sub> UV<sup>⊤</sup></td><td>[16]</td></tr><tr><td>Stoch. Spectral descent with momentum</td><td>[0, 1]</td><td>Layerwise Spectral ||·||<sub>S<sub>∞</sub></sub></td><td>||σ||<sub>1</sub> UV<sup>⊤</sup></td><td>[27]</td></tr></tbody></table>

To see how uSCG can in general be considered a normalized version of steepest descent, we point to the fact that the sharp operator and lmo can be defined in terms of each other. Specifically, we have that

$$[d]^# = -\|d\|_* lmo(d). \qquad (5)$$

From this relationship we see that the lmo can be expressed as the normalized sharp-operator.

*Examples:* We can also deduce from the previous relationship that uSCG with $l_\infty$-norm constraints leads to the SignSGD update [17]: $x^{k+1} = x^k - \gamma_k \text{sign}(d^k)$. This scheme was made popular in the distributed optimization literature due to the low communication cost of sending the bit-representable vector $\text{sign}(d^k)$ as oppose to the floating point vector $d^k$. Similarly, for the spectral norm, uSCG reduces to $x^{k+1} = x^k - \gamma_k U^k (V^k)^\top$. When $d^k$ is a momentum accumulation of past gradients, e.g., $d^k = (1 - \alpha_k) d^{k-1} + \alpha_k g^k$, we recover the primary operation used in the Muon implementation [2].

From (5) we also see a clear distinction between the lmo and the sharp operator. While the lmo is scale invariant (i.e., $lmo(a \cdot s) = lmo(s)$ for $a > 0$), the sharp operator is scale equivariant (since $[a \cdot s]^# = a[s]^#$ for $a \in \mathbb{R}$). The scale or affine invariance of conditional gradient methods makes the family of methods agnostic to the gradient size. In particular, lmo outputs can always be chosen to lie on the boundary, i.e., $\|lmo(s)\| = 1$, such that the magnitude of the update in uSCG can be predetermined.

<footnote n="1">¹ With non-Nesterov based momentum and ignoring the prescribed use of AdamW for some parameters.</footnote>
**Norm control of Conditional Gradient Methods:** Both uSCG and SCG provide a control of the norm of the iterates through $||x^k|| \le \rho \sum_{i=1}^k \gamma_i$ and $||x^k|| \le \rho$ for all $k \in \mathbb{N}$, respectively. Specifically, SCG solves the constrained problem

$$ \min_{x \in \mathcal{D}_{\rho}} f(x). \qquad (6) $$

Importantly, the control on the iterates is in the chosen norm. This norm control acts as a regularizer which turns out to be superior to placing an (Euclidean) $\ell_2$-regularization on the objective itself, which has led to the popularity of the AdamW optimizer, which replaces the $\ell_2$-regularization present in the Adam optimizer with decoupled weight decay.

In practice, controlling the norm is not only important for avoiding overfitting but also for providing numerical stability, especially for long training runs or under low-precision. Long training runs are particularly common for generative adversarial networks (GANs), diffusion models, and in the pre-training of large language models. We illustrate the importance of the norm control in Figure 3. A similar benefit of decoupled weight decay can be observed for spectral methods as shown in [1], [28], which was later used to train the 1T parameter Kimi-K2 model through Muon with decoupled weight decay.

**Frank-Wolfe vs. PyTorch weight decay:** In SCG, the update employs three distinct scalars: a learning-rate schedule multiplier $s_k$, a Frank-Wolfe stepsize $\lambda$ (with $\lambda_k = s_k \lambda$), and a constraint radius $\rho = \gamma/\lambda$. As we will see, under $\mu$-parameterization, we will choose scalings of the hyperparameters like the learning rate layerwise to zero-shot transfer them to larger models, which is known as $\mu$-transfer [29]. In the Frank-Wolfe view, the stepsize $\lambda$ is treated separately from the base radius $\rho$ (and schedule $s_k$), which will enable a clean hyperparameter transfer. PyTorch's default implementation of weight decay instead folds these into each other, and changing $lr$ for $\mu$-transfer unintentionally changes the effective decay:

$$ x^{k+1} = (1 - lr \cdot \text{weight\_decay})x^k + lr \text{lmo}(d^k). \qquad (\text{SCG-PyTorch}) $$

The fix is trivial: treat $1/\text{weight\_decay}$ as the constraint radius $\rho$ and enforce $lr_k = s_k \lambda/\text{weight\_decay}$. With this, PyTorch and the Frank-Wolfe formulation are equivalent and both $\mu$-transfer correctly.

For any positive $\text{weight\_decay}$, it is similarly possible to interpret many other normalized algorithmic updates with PyTorch weight decay added within the SCG framework (cf., Table II). In particular, adding the weight decay simply modifies the constraint by scaling the radius to get an effective constraint set $\mathcal{D}_{1/\text{weight\_decay}}$ whose radius is $\frac{1}{\text{weight\_decay}}$ and a modified stepsize $lr \cdot \text{weight\_decay}$.

**Momentum and Dual Averaging:** Momentum, as we have defined by

$$ d^k = (1 - \alpha_k) d^{k-1} + \alpha_k \nabla f(x^k, \xi_k), \qquad (7) $$
<figure><img src="image_5.png" /><figcaption>Fig. 3. (left) The validation accuracy of a Vision Transformer on ImageNet benefits substantially from the tighter norm control provided by an instance of SCG (Scion) compared with uSCG (Unconstrained Scion) [1]. (right) Momentum sweep on NanoGPT, where the stepsize $\gamma$ is retuned for each configuration. Performance reduces as momentum is disabled ($\alpha \to 0$).</figcaption></figure>

has a clear motivation in the non-Euclidean setting [1], [9], [25] since both the lmo and the sharp-operator may be biased. This can happen even if the SFO, $\nabla f(x^k, \xi_k)$ is unbiased, because of the nonlinearity of the lmo or the sharp-operator. With careful choice of vanishing $\alpha_k$, the momentum estimator acts as variance reduction and can be used to guarantee convergence of uSCG [1]. In this sense, momentum provides a practical alternative to increasing batch sizes, which may not otherwise be practical.

The empirical benefits of momentum are evidenced by its presence in countless popular optimizers for deep learning. This use of momentum is also well-motivated theoretically. For the generic template SCG and uSCG with a constant stepsize $\gamma = \Theta(\frac{1}{n^{3/4}})$ and vanishing momentum $\alpha_k = \frac{1}{\sqrt{k}}$ it is possible to show convergence to a first-order stationary point with the rate

$$O(\frac{1}{n^{1/4}} + \frac{L}{n^{3/4}}), \qquad (8)$$

for $\bar{x}^k$ sampled uniformly over the iterates $\{x^k\}_{k \in [n]}$ [1]. Without momentum ($\alpha_k = 1$) the algorithm would only be guaranteed to converge to a noise-dominated region [1].

What does the rate in (8) tell us? Since the $O(n^{-1/4})$ term asymptotically dominates, to see the benefit of a favorably Lipschitz constant $L$, we need to be in a regime where either the noise is small enough (i.e., by taking large batch sizes) or where the horizon $n$ is small. This aligns well with practice where methods exploiting non-Euclidean geometries, such as SignSGD, Adam, Muon and Scion, typically benefit from the regime of large batch size and short runs.

As mentioned, the momentum estimator plays a central role for establishing convergence by reducing the variance. This leads one to wonder if we can design more sophisticated estimators? The space of
possible dual feedback mechanism is indeed rich, which we illustrate with two prominent examples:

$$ \text{Lion :} \quad \begin{cases} \tilde{d}^k = (1 - \alpha_k) \tilde{d}^{k-1} + \alpha_k \nabla f(x^k, \xi_k), \\ d^k = (1 - \beta_k) \tilde{d}^{k-1} + \beta_k \nabla f(x^k, \xi_k). \end{cases} $$

$$ \text{STORM : } d^k = (1 - \alpha_k)d^{k-1} + \alpha_k \nabla f(x^k, \xi_k) + (1 - \alpha_k)(\nabla f(x^k, \xi_k) - \nabla f(x^{k-1}, \xi_k)). $$

The estimator in the Lion optimizer was discovered empirically through a symbolic search over algorithms. Notably, it is similar to the estimator in optimistic follow the regularized leader in online learning which uses $\nabla f(x^k, \xi_k) + \sum_{i=1}^k \nabla f(x^i, \xi_i)$. Lion with $\alpha_k = \beta_k$ is also commonly referred to as Nesterov momentum and is used, e.g., in Muon. Empirically, the Lion estimator has been observed to perform better than the momentum estimator in (7) (recovered with $\beta_k = 0$), while theoretically its mechanism remains to be explained.

On the other hand, the STORM estimator was originally discovered theoretically as a variance reduction mechanism for accelerating the convergence rate on nonconvex problems by improving the classical $O(n^{-1/4})$ rate e.g., appearing in (8) to $O(n^{-1/3})$. It was later used in the MARS optimizer, which used a stale version of STORM in order to avoid computing two gradients per iteration (i.e., $\nabla f(x^{k-1}, \xi_k)$ was replaced by the previously computed $\nabla f(x^{k-1}, \xi_{k-1})$). Further improvements could be expected in finite-sum problems (i.e., the multi-epochs setting) by for example integrating the SPIDER estimator as in SPIDER-FW [30].

Two themes in the design of dual feedback are averaging gradients and averaging iterates, which are equivalent in the Euclidean setting, since the lmo is then linear. In the non-Euclidean setting, the nonlinearity of the lmo means the two are no longer equivalent, and the role of each becomes more clear. As we shall see next, averaging the iterates is important in producing the final model that is used.

**Selecting a candidate solution:** An attentive reader might have noticed that the guarantee in (8) for the stochastic case holds only for the average over the trajectory. There are roughly two paradigms for providing guarantees for a final, candidate iterate: either (i) a new sequence of iterates is produced (e.g., by averaging), or (ii) the stepsize schedule is modified. Maybe surprisingly, the latter approach can also be viewed as a different way of averaging the iterates.

A classical approach under convexity is to output the Polyak-Ruppert average, $\hat{x}^n := \frac{1}{n} \sum_{k=1}^n x^k$, which is justified by Jensen's inequality. In the non-convex case, this is unfortunately no longer justified theoretically, although variants of this are still used. Similar to RMSProp, it is common in deep learning to instead output an exponentially moving average (EMA), $\hat{x}^{k+1} = \tau \hat{x}^{k+1} + (1-\tau)x^k$ where typically $\tau \in \{0.999, 0.9999\}$.
An alternative approach, when the horizon is known, is through the choice of the stepsize schedule. Popular choices includes the cosine decay and the linear decay, both of which importantly decay the stepsize to zero (i.e., $\gamma_n = 0$). For the linear decay $\gamma_k = \gamma(1 - k/n)$, the schedule is well-motivated in the convex setting, as it provably leads to last iterate guarantees for both non-smooth and stochastic problems [31], [32]. A practical advantage of this scheduling approach is that no additional memory is required, as would otherwise be the case when maintaining an averaged iterate like $\hat{x}^{k+1}$.

Remarkably, we can generalize and interpolate between these two seemingly distinct approaches [32]. The generalization is obtained through the online learning framework and leads to the so-called schedule-free method, which in addition to the iterates $x^k$ themselves maintains an average over the iterates $z^k$:

$$y^k = (1 - \beta_k)x^k + \beta_k z^k, \quad z^{k+1} = (1 - c_{k+1})z^{k+1} + c_{k+1}x^{k+1}$$

where $c_{k+1} = \frac{1}{k+1}$ and $x^{k+1}$ is produced by some base optimizer receiving the feedback $\nabla f(y^k, \xi_k)$, e.g., SCG. We recover the mentioned Polyak-Ruppert averaging with $\beta_k = 0$, Polyak momentum with $\beta_k = 1$, and the linear decay stepsize schedule with $\beta_k = k/n$ where $n$ is the horizon.

As suggested by the dependency on $n$, linear decay only guarantees good performance for the final iteration at the end of training, in contrast with a constant $\beta_k = \beta$, which provides an any-time guarantee. Empirically, it has been observed that $\beta \approx 0.9$ works well, but the role of $\beta$ is still not fully understood.

A reader familiar with optimization might recognize the schedule-free method is essentially a repa-rameterization of Nesterov momentum. Interestingly, this perspective provides a connection to recent work on Lookahead-based optimizers, which similarly wraps a base optimizer, but executes multiple inner iterations before applying the outer iteration. DiLoCo uses Nesterov momentum for the outer loop and has been observed to improve performance significantly even in the single-client case, despite being developed for the distributed (multi-client) case.

Apart from this * cooldown phase where the stepsize decays, there are typically one or two other phases in a scheduler. First, it is common to have a * stable phase before cooldown in which the stepsize is taken constant. Before this stable phase, an initial linear warmup of the stepsize for SGD and Adam is also common. Interestingly, the warmup phase does not seem to be necessary for optimizers with a priori adaptation to a favorable geometry, nor for long runs when using Adam. These common schedules are visualized in Figure 4.

<figure><img src="image_6.png" /><figcaption>Fig. 4. Typical learning rate schedules: linear and cosine decay.</figcaption></figure>
**Computational complexity:** The computational cost of steepest descent and conditional gradient methods depends on the specific norm that induces the lmo. The sharp-operator used in steepest descent is only as expensive as the lmo, since the additional dual norm scaling can be computed as $-\|d\|_* = \langle d, \text{lmo}(d) \rangle = \text{vec}(d)^{\top} \text{vec}(\text{lmo}(d))$ due to the definition of the dual norm and the optimality of lmo(d). When $\|\cdot\|$ is a vector norm, the computation of the lmo is typically free, as is the case for $\ell_\infty$-norm (SignSGD) and $\ell_2$-norm (Normalized SGD).

However, the cost of the lmo can be significant when considering matrix norms. In particular, let us focus on the lmo for the spectral norm, which is used in many modern optimizers (Table II):

$$U \operatorname{diag}(\sigma) V^{\top} \mapsto UV^{\top} \in \mathbb{R}^{p_2 \times p_1}$$

One way of computing the spectral lmo is through the (reduced) singular value decomposition (SVD), which in the worse case has a slow $O(p_2p_1 \cdot \min(p_2, p_1))$ time complexity. Efficient alternatives to computing $UV^\top$ have fortunately been studied for decades. We highlight three prominent families:

* *Matrix-multiplication based* routines such as Newton-Schultz and PolarExpress [33]. One major benefit is that they are GPU friendly and work in low-precision. Newton-Schultz is the solver used in Muon, where the coefficients of the solver were hand-tuned to find a coarse approximation fast. PolarExpress further optimized these coefficients analytically.

* *QR based* routines such as power iteration used in the Dion optimizer [34], QDWH and ZOLO-PD, which have favorable properties in distributed settings, since the QR decomposition can be parallelized across devices. One downside is that the QR decomposition requires full precision.

* *Sketching*, which is a randomized method that constructs a low-rank approximation to the SVD. This was the solver used originally in the stochastic spectral descent method (SSD) [16].

Sketching and QR based methods have the additional benefit of enabling variants with efficient communication between devices by being able to truncate $U$ and $V^\top$, since they are materialized. Another difference with the other methods is that, if the input is low-rank, then the output remains low-rank, which is in contrast with Newton-Schultz and PolarExpress for which the rank may increase to full.

## SCALING DEEP LEARNING MODELS

In the previous section we have covered how optimization methods can be designed to exploit the geometry of the optimization problem in general. In this section, we focus on neural networks and discuss how one can leverage theoretical tools to identify effective strategies to scale up deep learning models. We will see that such scaling rules often yield desirable practical benefits such as monotonic improvement in performance with scale and the transferrability of certain hyperparameters such as learning rate from
small to large scale models [29]. We will also see how considering the propagation of signals through
the neural networks guides our choice of norm for optimizing the network.

## Background

Deep learning models are increasingly deployed at large scales. This is currently driven by the emergence of empirical neural scaling laws, which show that model performance predictably improves with scale [35] and follows a simple power-law relationship with the amount of available compute. These scaling laws play a crucial role in the design of foundation models in practice.

Since training at such scales typically requires substantial financial and computational resources for even a single run, design choices are primarily informed by extrapolating from smaller models, allowing researchers to compare alternative design choices—such as Adam versus Shampoo—without exhaustive large-scale experiments. Beyond improving performance, increasing scale also appears to induce qualitatively novel behaviors, known as emergent capabilities, suggesting that scaling is not merely a quantitative phenomenon [37].

Crucially, however, scaling does not always yield better performing models, and in general arbitrary scaling can hurt model performance. For instance, [36] shows that if we increase the width of a Mamba model [38] according to the

default scaling rule used in practice, performance of the model deteriorates after a certain width threshold (see Figure 5 for an illustration). Examples such as these are plentiful in literature [29].

<figure><img src="image_7.png" /><figcaption>Fig. 5. Figure from [36]. Test loss against learning rate on Mamba with varying widths under standard scaling. Performance of the optimally tuned model starts to worsen after a certain width.</figcaption></figure>

To ensure improved model performance with scale, we need *principled scaling rules*.

As the available computational budget increases, the training configuration of a neural network can be adapted in countless ways. These include scaling dimensions such as network width, depth, dataset size, or batch size. A scaling rule prescribes how key aspects of training—such as the architecture, hyperparameters, and dataset—should be adjusted as the computational budget varies.

To scale neural networks effectively, it is essential to establish scaling rules that satisfy at least two key desiderata:

1) *Optimality*: The large-scale model should be highly performant relative to alternative scaling rules.
2) *Predictability*: Performance of large-scale models can be reliably inferred from small-scale models, facilitating efficient exploration of the design space.

*Scaling limits* serve as a powerful theoretical framework to derive such scaling rules and more generally to understand the behaviour of large scale models. Scaling limits are the limits of neural networks obtained by taking various scaling dimensions such as the width or depth of the network to infinity.

In principle, to achieve optimality, one can derive the scaling limits associated with different scaling rules and determine which ones yield the most favorable asymptotic properties, such as stability of the training dynamics and feature learning. Predictability, on the other hand, requires that models rapidly converge to their scaling limits, motivating an analysis of the rate of convergence under various scaling rules.

Given the complexity of deriving scaling limits for arbitrary scaling strategies, analyses often constrain the space of scaling rules under consideration. In this section, we primarily focus on scaling rules where only the width of the network is taken to infinity, through there is growing literature covering the increasing depth as well [39], [40].

Developing theoretical tools to analyze these scaling limits has been an ongoing challenge in deep learning theory. One influential approach in this direction is the Neural Tangent Kernel (NTK) framework [4], which allows us to study the behaviour of infinitely wide neural networks under a specific choice of width scaling, known as the Neural Tangent Kernel Parameterization. NTK theory has provided fundamental insights into deep learning dynamics. However, a key limitation of NTK analysis is its inability to capture *feature learning*. In this regime, the parameters remain close to their initialization, and training dynamics reduce to kernel regression with a fixed, data-independent kernel. Consequently, the underlying feature representation is determined solely by the architecture and initialization rather than learned from data, which restricts the relevance of NTK theory to practical deep learning settings.

Early attempts to circumvent this limitation focused on mean-field theory, particularly for two-layer networks [5]. More recently, the Maximal Update Parameterization (μP) has emerged as a principled generalization to deep architectures. Introduced by [6], μP is based on an expanded space of hyperparameters that allows for different learning rates and initialization variances in each layer of the neural network. This parameterization enables nontrivial feature evolution even in the infinite-width limit, thus bridging the gap between theory and practical deep learning.

The foundation of μP is the Tensor Programs framework [6], which provides a systematic way to analyze the asymptotic behavior of wide neural networks at both initialization and during training. Tensor Programs are defined through a structured sequence of operations, including matrix multiplications, non-linear outer products, and vector averages. This framework allows for a precise tracking of neural network
dynamics, leading to the discovery of μP as the unique scaling rule that guarantees feature learning for every layer, regardless of width. Furthermore, empirically, μP has been shown to enable the transferability of optimal hyperparameters like learning rate from small to large-scale models, a property that is highly desirable in large-scale deep learning applications.

A complementary perspective on feature learning is rooted in operator norms. The authors of [41] show that effective feature learning is achieved by scaling the spectral norm of weight matrices and their updates proportionally to the square root of the ratio between output dimension and input dimension. The spectral condition recovers μP and, importantly, motivates the non-Euclidean optimization methods covered in the previous section.

**Outline:** This section first introduces the abc-parameterization and identifies the Maximal Update Parameterization (μP). We then provide an intuitive understanding of μP through the spectral perspective by focusing on the infinite-width limit. From this perspective we will clearly see the implication for initialization, norm constraint and the importance of the low-rank structure of gradients. We close by discussing practical aspects of hyperparameter transfer, the connection with Lipschitz continuity, and the norm choice for the optimizer.

## Identifying optimal scaling limits

Consider the setting of an L-hidden-layer multilayer perceptron (MLP) defined iteratively via

$$h_1(z) := W_1 z, \quad o_\ell(z) := \Phi(h_\ell(z)), \quad h_{\ell+1}(z) := W_{\ell+1} o_\ell(z), \quad \Phi(z) := W_{L+1} o_L(z),$$

for inputs $z \in \mathbb{R}^{p_1}$ with weight matrices $W_1 \in \mathbb{R}^{p \times p_1}$, $W_\ell \in \mathbb{R}^{p \times p}$ for $\ell \in [2, L]$, and $W_{L+1} \in \mathbb{R}^{p_{L+1} \times p}$.
Given a loss function $\mathcal{L}(\Phi(z), y)$, we use the SGD update rule (SGD) to update the parameters of the MLP

$$W_{\ell}^{k+1} = W_{\ell}^{k} - \eta_{\ell} \nabla_{W_{\ell}^{k}} \mathcal{L}(\Phi(z), y). \qquad (\text{SGD})$$

To understand the scaling behaviour of MLPs trained under SGD, [6] introduced a natural class of width-scaling rules—abc-parameterizations—which allow layerwise scalings of the parameters, the variance of their initializations, and learning rates by introducing width-dependent scaling factors. Throughout, we will use $\Theta(\cdot)$ to represent the asymptotic order in the infinite width limit $p \to \infty$.

**Definition 1 (abc-parametrization [6]).** An abc-parametrization $\{a_\ell\}_{\ell \in [L+1]} \cup \{b_\ell\}_{\ell \in [L+1]} \cup \{c_\ell\}_{\ell \in [L+1]}$ parameterizes an MLP trained with SGD in the following way:

1) Parameterize the weight matrices $W_\ell$ of the MLP as $p^{-a_\ell} w_\ell$, for a trainable weight parameter $w_\ell$ of the same shape as $W_\ell$, for all $\ell \in [L+1]$.
2) Initialize the trainable weights as $w_{\ell} \sim \mathcal{N}(0, p^{-2b_{\ell}}I)$.

3) Train the weights using the SGD update rule with layerwise learning rates for some width independent $\eta$,

$$w_{\ell}^{k+1} = w_{\ell}^{k} - \eta p^{-c_{\ell}} \nabla_{w_{\ell}^{k}} \mathcal{L}(\Phi(z), y),$$

Most standard initialization schemes used in practice fall within this family of scaling rules. For instance, He/LeCun initializations are specified by $a_{\ell} = 0 \forall \ell \in [L+1]$, $b_1 = 0$ and $b_{\ell} = 1/2$ for $\ell > 1$, and $c_{\ell} = 0 \forall \ell \in [L+1]$, i.e.,

$$W_1 \sim \mathcal{N}(0, 1/p_1), W_{\ell} \sim \mathcal{N}(0, 1/p), \text{ for } \ell \in [2, L], W_{L+1} \sim \mathcal{N}(0, 1/p),$$

with a global learning rate $\eta_{\ell} = \eta$ for all $\ell \in [L+1]$. This parametrization, commonly referred to as the Standard Parametrization (SP), fails to yield width-independent feature evolution in the infinite-width limit when trained with Mean-Squared Error (MSE) loss.

To identify optimal scaling limits, one can analytically derive the infinite width limits of MLPs trained with SGD under this family of abc parameterizations and classify all possible limits according to their qualitative behaviour. To this end, the Tensor Program (TP) framework introduced by [6] provides a way to derive the infinite-width limits of neural networks under arbitrary abc-parameterizations.

As the width of hidden layers approaches infinity, Tensor Programs identify two distinct dynamical regimes among others:

(i) Kernel Regime: here, weight updates are small relative to initialization, and the network behaves like a kernel method where no features are learned. The limiting behavior is governed by a fixed kernel derived from the initialization. Formally, an abc-parametrization is in the kernel regime if there exists a positive semidefinite kernel $K$ such that, in the $p \to \infty$ limit,

$$\Phi^{k+1}(z) = \Phi^k(z) - \eta K(z, z^k) \mathcal{L}'(\Phi^k(z^k), y^k), \quad \forall k \ge 0. \quad (9)$$

(ii) Feature Learning Regime: weight updates are scaled to induce stable feature evolution, enabling the network to learn data-dependent representations. Formally, an abc-parametrization admits feature learning if, as $p \to \infty$, each entry of $\Delta o_L(z)$ is $\Theta(1)$ where $\Delta o_L(z)$ denotes the change in the last layer activation vector $o_L(z)$ after one step of SGD.

Maximal Update Parameterization (μP): Under the feature learning regime, there exists a unique abc-parameterization (up to an equivalence class) that enables maximally stable feature learning. Formally, for all $\ell \in [L]$, the update term $\Delta W_{\ell} o_{\ell-1}(z)$ has $\Theta(1)$ coordinates. In other words, weight updates in every layer induce a non-trivial effect on the network's output. As feature learning is considered an useful mechanism for generalization, one may argue that μP is an optimal scaling rule within this class.
For MLPs trained with SGD, μP can be specified by the parameter choices:

$$ 
\begin{array}{l}
a_{\ell} = 0, \quad \forall \ell \in [L+1], \quad b_1 = 0, \quad b_{\ell} = \frac{1}{2}, \quad \text{for } \ell \in [2, L], \quad b_{L+1} = 1, \\
c_1 = -1, \quad c_{\ell} = 0, \quad \text{for } \ell \in [2, L], \quad c_{L+1} = 1.
\end{array}
$$

This corresponds to the following initialization and learning rate scalings:

$$ \text{Initialization: } W_1 \sim \mathcal{N}(0, 1/p_1), \quad W_{\ell} \sim \mathcal{N}(0, 1/p) \quad \text{for } \ell \in [2, L], \quad W_{L+1} \sim \mathcal{N}(0, 1/p^2). \qquad (10) $$

$$ \text{Layerwise SGD learning rates: } \eta_1 = \eta p, \quad \eta_{\ell} = \eta, \quad \text{for } \ell \in [2, L], \quad \eta_{L+1} = \eta p^{-1}. $$

Alternative parameterizations that achieve maximally stable updates can be derived by leveraging the symmetries of abc-parameterizations (see [6]). Specifically, for any $\alpha \in \mathbb{R}$, shifting the parameters as

$$ a_{\ell} \rightarrow a_{\ell} + \alpha, \quad b_{\ell} \rightarrow b_{\ell} - \alpha, \quad c_{\ell} \rightarrow c_{\ell} - 2\alpha $$

does not affect the training dynamics. That is, the activations and the function computed by the network remain unchanged.

μP scaling offers several practical benefits. It often tends to improve performance and performance tends to monotonically improve with scale. Empirically, it has been shown that, under μP, many properties of wide neural networks including loss curves, pointwise predictions, pre-activation distributions, internal representations etc tend to be consistent across large widths compared to other parameterizations, such as the Neural Tangent Kernel (NTK) parameterization, reflecting better predictability.

Another particularly interesting practical benefit is that under μP, certain optimal hyperparameters such as learning rate can be transferred from small to large-scale models, which is far from standard practice. This also suggests that under μP, finite-sized models are potentially closer to their infinite-width limit.

*Deriving optimal scaling rules for other architectures and optimizers:* While our discussion has evolved around MLPs, the framework of Tensor Programs can be utilized to analytically derive infinite width limits and accordingly identify optimal<fnref n="1" /> scaling limits for a fairly broad class of architectures, such as ResNets and Transformer, and optimizers such as Adam. However, there are notable exceptions. Certain modern state space models, such as the Mamba architecture, cannot be represented as a Tensor Program. Rigorously establishing their infinite-width limits requires a significant generalization of the existing framework [36].

<footnote n="1">¹It is not clear a priori what the notion of optimality entails (see for instance [42]).</footnote>
Infinite-width limit & the spectral condition

Why are different layers treated differently in μP? To build intuition, we will see how feature learning
can follow from a spectral conditions on the weights and updates in the special case of MLPs and
normalized updates [41]. We will focus on the simplified case of linear MLPs, where the number of
assumptions needed is minimal. We define the initial hidden layer by $h_1(z) = W_1z$ and the remaining
layers by $h_\ell(z) = W_\ell h_{\ell-1}(z), \forall \ell \in 2, \dots, L+1$. We denote the global loss as $\mathcal{L}(h_L(z), y)$ where $\mathcal{L}$ is
the loss function and $y$ is a 1-hot encoded target vector. The RMS norm, $\|z\|_{\text{RMS}} := \frac{1}{\sqrt{p}}\|z\|_2$ for $z \in \mathbb{R}^p$,
will play a central role. We can redefine feature learning in terms of the RMS norm:

**Definition 2 (Feature Learning).** Let $\Delta h_{\ell}(z)$ denote the change in the layer-$\ell$ activations. Then we are said to be in the feature learning regime if the following holds for all $\ell \in [L+1]$:

(i) $\|h_{\ell}(z)\|_{\text{RMS}} = \Theta(1)$ (the forward pass is stable)

(ii) $\|\Delta h_{\ell}(z)\|_{\text{RMS}} = \Theta(1)$ (the update is non-vanishing and non-exploding)

The definition states that the typical element of the vectors should be of order $\Theta(1)$ as $p$ goes to infinity. If the entries were larger, then they would blow up with width, and thus become numerically unstable. If they were any smaller this could likewise lead to instability. Additionally, if $\|\Delta h_{\ell}(z)\|_{\text{RMS}}$ is not lower bounded, and thus potentially vanishing with width, the model may be in the less expressive Neural Tangent Kernel (NTK) regime where features do not evolve.

What we control in practice are not the preactivations $h_{\ell}(z)$ themselves, but rather the weight matrices $W_{\ell}$ that induce them. Fortunately, we can upper bound $\|h_{\ell}(z)\|_{\text{RMS}}$ through a norm control on $W_{\ell}$. Recall the operator norm defined as follows

$$
\|A\|_{\alpha \to \beta} := \max_{\omega \in \mathbb{R}^p, \omega \neq 0} \frac{\|A\omega\|_{\beta}}{\|\omega\|_{\alpha}} = \sup_{\|\omega\|_{\alpha}=1} \|A\omega\|_{\beta}. \quad (11)
$$

Directly from the definition, we see that if the input $\omega$ is bounded through $\|\omega\|_{\alpha} \le 1$, then the output $\|A\omega\|_{\beta}$ will be bounded when $\|A\|_{\alpha \to \beta}$ is bounded. Thus, $\|W_{\ell}h_{\ell-1}(z)\|_{\text{RMS}}$ can be upper bounded through an upper bound on $\|W_{\ell}\|_{\text{RMS} \to \text{RMS}}$, which is just a rescaled version of the Schatten-$\infty$ norm (a.k.a. the spectral norm), i.e., $\|W_{\ell}\|_{\text{RMS} \to \text{RMS}} = \sqrt{\frac{p_{\text{in}}}{p_{\text{out}}}} \|W_{\ell}\|_{S_{\infty}}$. The only additional requirement is that the original input to the network is bounded, i.e., $\|z\|_{\text{RMS}} < \infty$. A similar reasoning applies for upper bounding $\|\Delta h_{\ell}(z)\|_{\text{RMS}}$.

If we additionally knew that these bounds were tight, then feature learning could be restated solely in
terms of the weight matrices as follows:

**Definition 3 (Spectral condition).** The spectral condition requires that
(i) $||W_\ell||_{\text{RMS}\to\text{RMS}} = \Theta(1)$ (the weight matrices are properly initialized)

(ii) $||\Delta W_\ell||_{\text{RMS}\to\text{RMS}} = \Theta(1)$ (the updates are properly scaled at each iteration)

**Implication for initialization:** The spectral condition needs to hold for all the iterates throughout optimization. At initialization, condition (i) can be trivially satisfied by directly initializing $W_\ell$ to be a scaled semi-orthonormal matrix or through Gaussian initialization, $[W_\ell]_{ij} \sim N(0, \sigma_\ell^2)$, since random matrix theory tells us that the spectral norm scales as $||W_\ell||_{2\to2} \approx \sigma_\ell(\sqrt{p_{\ell-1}} + \sqrt{p_\ell})$. This motivates the following variance choice (cf., [41, Param. 1]):

$$ \sigma^2 = \Theta(\frac{1}{p_{\ell-1}} \min\{1, \frac{p_\ell}{p_{\ell-1}}\}). $$

For language modeling, where the vocabulary size (i.e., the output dimension for the last layer) is much larger than the network width (i.e., input dimension for the last layer), the variance requirement above implies that even the last layer can be initialized with $\sigma^2 = \Theta(\frac{1}{p_L})$, which is satisfied by e.g., the popular Kaiming initialization [43] otherwise typically associated with SP (see e.g., Figure 1, which successfully uses SP initialization for hyperparameter transfer).

**Low-rank structure of gradients:** How do we show that the upper bounds provided through the operator norm are tight, thereby proving that the spectral condition implies feature learning? The bounds can be particularly loose when the incoming vector $h_{\ell-1}(z)$ interacts with small singular values of the matrices $W_\ell$ and $\Delta W_\ell$. Fortunately, $h_{\ell-1}(z)$ turns out to be aligned with the top singular vectors of the matrices. The key insight is that the gradients happens to be rank-1 with right singular vector $h_{\ell-1}(z)$ when the batch size is one, so the alignment trivially holds.

Due to the rank-1 structure of the gradients, many update rules coincide for batch size one. For instance, in this case the Frobenius norm and spectral norm of the gradients are equivalent, so the two corresponding lmo-based updates are equivalent. Furthermore, both are equivalent to normalizing with the spectral norm. Thus, all three of these methods provide valid ways of ensuring that $||\Delta W_\ell||_{\text{RMS}\to\text{RMS}} = \Theta(1)$ in the single sample case.

For batch sizes greater than one, the equivalence between the methods no longer strictly holds and establishing alignment becomes more difficult. However, despite not being rank-1, it has been observed that the stable rank remains small even for large batch sizes in practice.

**The importance of layerwise treatment:** One central takeaway from $\mu$P and the preceding reasoning in the infinite width limit is that weights should be optimized differently depending on which layer they correspond to in the network (e.g., first vs last). Focusing on the spectral condition (Definition 3), which requires that $||\Delta W_\ell||_{S_\infty} = \Theta(\sqrt{\frac{p_\ell}{p_{\ell-1}}})$, there are three distinct cases:

(i) $||\Delta W_1||_{S_\infty} = \Theta(\sqrt{p_1})$ (the input scaling grows)
(ii) $||\Delta W_l||_{S_\infty} = \Theta(\sqrt{\frac{p_l}{p_{l-1}}})$ for all $l \in \{2, ..., L+1\}$ (the hidden layers are constant with width)

(iii) $||\Delta W_{L+1}||_{S_\infty} = \Theta(\frac{1}{\sqrt{p_L}})$ (the output scaling goes to zero)

The distinction is a simple consequence of the input and output dimension being kept fixed, but it has nonetheless stark implications. For example, while the hidden layers should have constant magnitude updates, we see that the scaling of the last layer updates should instead approach zero as $p_L \to \infty$.

This reasoning through the spectral condition highlights many of the key ingredients, such as the choice of initialization and the importance of the low-rank structure of gradients. However, the complete argument relies on a technical assumption stating that a certain cancellation of terms never occurs, requires additional assumptions to treat activation functions, and is restricted to MLPs and normalized updates. The Tensor Program framework provides a way to avoid these assumptions, and it is flexible enough to cover a large family of both networks and algorithms.

*Verifying feature learning:* Suppose that we have parameterized our model and algorithm such that we should expect feature learning. How do we verify in practice that our updates are indeed in the feature learning regime, and not instead unstable like SP or in the lazy regime?

The fact that hyperparameters transfer across model sizes can be an indication, but it can be misleading, since transfer might be observed without all layers exhibiting feature learning. The appropriate test is instead the so-called coordinate check, which is also significantly cheaper to perform.

In essence, a coordinate check simply verifies the feature learning definition, which states that the coordinates of preactivations throughout the network should remain constant with width. This property can be verified by running a few steps of the algorithm across different model sizes, and either computing the $l_{\text{RMS}}$-norm or the normalized $l_1$-norm of the preactivations, which measures the

magnitude of the typical entry and the magnitude of the average entry, respectively. In practice, these two metrics are very similar. An example is provided in Figure 6 showing the clear distinction between Standard Parameterization (SP) and the Maximal Update Parameterization ($\mu$P).

<figure><img src="image_8.png" /><figcaption>Fig. 6. Coordinate check on a 3-layer GPT model using the AdamW optimizer with SP (top) and $\mu$P (bottom). Preactivations and logits remains constant with width for $\mu$P as opposed to SP for which they blow up.</figcaption></figure>
## Norm choice

We have seen how feature learning can be ensured through the spectral condition. In this section, we cover how the requirements of the spectral condition connect to a norm choice on the full network parameters, and how that choice leads to an instantiation of the generic algorithmic framework introduced in the previous section. Specifically, we will see that Lipschitz conditions can be shown for all layers of the neural network in a given norm, when the weights are bounded in the same norm. The associated Lipschitz constant can importantly be made independent of the width, which provides a non-asymptotic argument for hyperparameter transfer. Finally, we touch on some of the limitations with this approach, since these Lipschitz conditions (upper bounds) are themselves not sufficient to rule out the lazy regime.

**A norm on the entire network:** The spectral condition (Definition 3) requires that the RMS → RMS norm of each layer's weight matrix is controlled and the same for each layer's update. The upper bound of these requirements can be satisfied simultaneously by choosing the max-norm on the joint parameters $x = \{W_\ell\}_{\ell \in [L+1]}$:

$$ \|x\|_{\mathcal{X}} = \max_{\ell \in [L+1]} \frac{1}{\rho_{\ell}} \|W_{\ell}\|_{W_{\ell}} \quad (12) $$

where each norm $\| \cdot \|_{W_\ell}$ is a layer-dependent operator norm, e.g., RMS → RMS. With this choice, any algorithm that guarantees $\|x\|_{\mathcal{X}} \le 1$ will ensure that $\|W_\ell\|_{W_\ell} \le \rho_\ell$ with a layerwise radius $\rho_\ell$ for all $\ell \in [L+1]$. Note, that in contrast with the first part on optimization algorithms, which only covers vector and matrix cases, this choice places a product norm over a collection of matrices and vectors.

The lmo-based schemes such as SCG decompose layerwise under this norm since

$$ \mathrm{lmo}_{\|\cdot\|_{\mathcal{X}}}(d) = \{ \rho_1 \mathrm{lmo}_{\|\cdot\|_{W_1}}(d_1), \dots, \rho_{L+1} \mathrm{lmo}_{\|\cdot\|_{W_{L+1}}}(d_{L+1}) \} $$

for $d = \{d_1, \dots, d_{L+1}\}$ where $\mathrm{lmo}_{\|\cdot\|_{\mathcal{D}}}(d) := \arg \min_{\|x\|_{\mathcal{D}} \le 1} \langle d, x \rangle$. The layerwise update magnitudes are therefore dictated by the radii, i.e., $\| \mathrm{lmo}_{\|\cdot\|_{W_\ell}}(d_\ell) \|_{W_\ell} = \rho_\ell$ for all $\ell \in [L+1]$, which in turn guarantees that the upper bound in spectral condition (Definition 3) holds for all iterates of SCG as long as $\lambda_k$ and $\rho_\ell$ are chosen in a width-independent way (we are free to choose these parameters in practice). This layerwise decomposition is used in the Scion optimizer, which chooses $\max_\ell \|W_\ell\|_{\mathrm{op}_\ell}$ for some layer dependent operator norms (see Table III) that can all be made to satisfy the upper bounds in Definition 3 with appropriate scaling. This layerwise decomposition is also present in the LARS and LAMB optimizers, which instead choose the norm $\max_\ell \|W_\ell\|_F$.

This max-norm over layers has the benefit of decomposing the lmo computation. The steepest descent methods Stochastic SD on the other hand, require computing the dual norm $\|d\|_{X_*} = \sum_\ell \|d_\ell\|_{W_*}$, which entails multiplying each layer's update by the global sum of the layerwise norms.
<caption>TABLE III</caption>

EXAMPLE OPERATOR NORMS AND THE ASSOCIATED lmoS OF A MATRIX $A \in \mathbb{R}^{p_{\text{out}} \times p_{\text{in}}}$. THE REDUCED SVD IS GIVEN AS
$A = U \text{diag}(\sigma)V^{\top}$, sign ACTS ELEMENTWISE, $\text{col}_j(A) := A_{\cdot,j}$ AND $\text{row}_i(A) := A_{i,\cdot}$. THIS TABLE IS NOT EXHAUSTIVE.

<table><thead><tr><th></th><th>1 → RMS (ColNorm)</th><th>1 → ∞ (Sign)</th><th>RMS → RMS (Spectral)</th><th>RMS → ∞ (RowNorm)</th></tr></thead><tbody><tr><th>Norm</th><td>max<sub>j</sub> <sup>1</sup>&frasl;<sub>&radic;<span style="text-decoration:overline;">p<sub>out</sub></span></sub> || col<sub>j</sub>(A) ||<sub>2</sub></td><td>max<sub>i,j</sub> |A<sub>i,j</sub>|</td><td>&radic;<span style="text-decoration:overline;">p<sub>in</sub>/p<sub>out</sub></span> ||A||<sub>S∞</sub></td><td>max<sub>i</sub> &radic;<span style="text-decoration:overline;">p<sub>in</sub></span> || row<sub>i</sub>(A) ||<sub>2</sub></td></tr><tr><th>LMO</th><td>col<sub>j</sub>(A) ↦ &minus;&radic;<span style="text-decoration:overline;">p<sub>out</sub></span> <sup>col<sub>j</sub>(A)</sup>&frasl;<sub>|| col<sub>j</sub>(A) ||<sub>2</sub></sub></td><td>A ↦ &minus;sign(A)</td><td>A ↦ &minus;&radic;<span style="text-decoration:overline;">p<sub>out</sub>/p<sub>in</sub></span> UV<sup>T</sup></td><td>row<sub>i</sub>(A) ↦ &minus; <sup>1</sup>&frasl;<sub>&radic;<span style="text-decoration:overline;">p<sub>in</sub></span></sub> <sup>row<sub>i</sub>(A)</sup>&frasl;<sub>|| row<sub>i</sub>(A) ||<sub>2</sub></sub></td></tr></tbody></table>

More importantly, the choice of the maximum norm over layers leads to an update of all layers at a given iteration. This is in contrast with, e.g., the choice of the $l_1$-norm made in the otherwise prescient work by Flynn [44] on dualization of neural networks, which leads to a (block) coordinate-wise update rule which modifies only one layer per step. The choice of maximum norm is one of the key design choices in the modular norm [3], [40], which uses the max-norm to generate layerwise scalings of the stepsize that can lead to hyperparameter transfer. In [40], the modular norm is formed exactly by taking the max-norm over layerwise norms. However, its main practical use was for layerwise stepsize normalization, which corresponds to assigning a scaled Euclidean norm to each layer.

**Layerwise operator norms:** It turns out that constraining the weight matrices in RMS → RMS norm is also a way to ensure that the layers are Lipschitz in the same norm, under certain assumptions on the input and the loss function computed on the output [40]. This operator norm perspective is the key idea behind recent works that led to Muon and Scion [1], [3]. In [3] it was proposed to assign operator norms to each layer and compute the lmo with respect to them (although no consideration was given to stochasticity or boundedness of the parameters). In [1], it was shown that same principle of assigning operator norms can be extended and used with the stochastic conditional gradient algorithm, to guarantee bounded parameters. If we move beyond assigning the RMS norm to the input and output space, we can constrain our weight matrix in a different operator norm and ensure the resulting layer is Lipschitz with respect to the new norms. Specifically, [3] proposed the column-wise RMS norm (ColNorm) for the input layer, while [1] proposed the row-wise RMS norm (RowNorm) for the output layer (cf., Table III).

**Revisiting constrained optimization:** The Lipschitz condition naturally requires, as an assumption, that the parameters themselves are contained within a norm-ball described above. This provides additional motivation for the norm-ball constrained optimization formulation (6) from the optimization section. Due to the max-structure of the norm choice in (12), the stochastic conditional gradient method (SCG), a.k.a. stochastic Frank-Wolfe, decomposes layerwise,

$$W_{\ell}^{k+1} = (1 - s_k \lambda) W_{\ell}^k + s_k \lambda \text{ lmo}_{\|\cdot\|_{W_{\ell}} \le \rho_{\ell}}(d_{\ell}^k) = (1 - s_k \lambda) W_{\ell}^k + s_k \lambda \rho_{\ell} \text{ lmo}_{\|\cdot\|_{W_{\ell}} \le 1}(d_{\ell}^k) \quad \forall \ell \in [L+1],$$
<figure><img src="image_9.png" /><figcaption>Fig. 7. The optimal stepsize and constraint radius transfers, when the norm constrained is picked appropriately layerwise for a CNN (left, middle). On the right, we see improved performance when training a vision transformer to classify ImageNet using a geometry-aware optimizer like Scion (plot from [1])</figcaption></figure>

where we have decomposed the dual feedback (a.k.a. the gradient estimator) as $d^k = \{d_1^k, ..., d_{L+1}^k\}$. Note that $\rho_\ell$ is the effective radius of the norm-ball constraint and that it is defined layerwise.

In particular, we see that the radius (and thus the regularization) of the last ((L + 1)th) layer increases with width, i.e., $\rho_{L+1} = \frac{\rho}{\sqrt{\rho_L}}$ for some width-independent scaling factor $\rho$, with the spectral norm choice from Definition 3. This is in contrast with hidden layers which have $\rho_\ell = \rho \sqrt{\frac{\rho_\ell}{\rho_{\ell-1}}}$, which remains a constant in the limit.

We are now going to think about the Lipschitz conditions as a non-asymptotic version of $\mu$P. Empirically, this $\mu$P scaling of the spectral norm turns out to lead to hyperparameter transfer as illustrated in Figure 7 on an image classification task on the CIFAR10 dataset using a convolutional neural network (CNN). What can be observed is that both the stepsize $\lambda$ and the radius scaling $\rho$ transfer across width. If the radius (a.k.a. the inverse of the weight decay parameter) was kept the same across layers and the stepsize instead varied over layers, then transfer should not be expected for the weight decay.

*Limitations:* The Lipschitz condition complements the infinite limit results, by, on the one hand, providing non-asymptotic results but, on the other hand, being unable to provide guarantees for feature learning. The second issue arises since Lipschitz conditions inherently only require an argument regarding upper bounds, which we already saw manifest itself when making the choice of the max norm over layers. Without a matching lower bound, we cannot guarantee that the lazy training regime is avoided.

It should also be noted that the Lipschitz conditions might be loose even in terms of their constants (ignoring their width dependency). In practice, this means that the computed constants are not used to determine Lipschitz-dependent hyperparameters, which are instead chosen through hyperparameter search.
## OUTLOOK

We have covered modern optimizers for deep learning, their geometries, and how to make them agnostic to the size of the model. Many questions remain open. While we have primarily focused on the invariance to width, there is the much broader question of how to scale optimally and predictably: How do we generally choose, e.g. the number of iterations, batch size, stepsize, model width, and model depth as a function of our compute budget? Embracing the non-Euclidean perspective, how can we build even faster algorithms? How can we refine the constrained choice for a broad class of neural network architectures? At massive scale, system-level concerns become the main bottleneck. For example, how to distribute the spectral lmo becomes a challenge, although research in this direction is active [27], [33], [34]. These are exciting questions that will require cross-disciplinary thinking—we hope that with this review article we have managed to pique the interest of the signal processing community.

## ACKNOWLEDGMENTS

We thank Jeremy Bernstein, Saleh Soltan, and Abhishek Kumar for helpful discussion. This work was supported as part of the Swiss AI Initiative by a grant from the Swiss National Supercomputing Centre (CSCS) under project ID a06 on Alps. This work was supported by the Swiss National Science Foundation (SNSF) under grant number 200021_205011. This work was supported by Hasler Foundation Program: Hasler Responsible AI (project number 21043). Research was sponsored by the Army Research Office and was accomplished under Grant Number W911NF-24-1-0048. ASF was supported by a public grant from the Fondation Mathématique Jacques Hadamard. This work was supported by the Gatsby Charitable Foundation, GAT 3850 and GAT 4058.

All code is available at https://github.com/LIONS-EPFL/scion

## REFERENCES

[1] T. Pethick, W. Xie, K. Antonakopoulos, Z. Zhu, A. Silveti-Falls, and V. Cevher, "Training deep learning models with norm-constrained LMOs," in *International Conference on Machine Learning*, 2025. 1, 2, 5, 7, 8, 9, 10, 11, 24, 25

[2] K. Jordan, Y. Jin, V. Boza, Y. Jiacheng, F. Cecista, L. Newhouse, and J. Bernstein, "Muon: An optimizer for hidden layers in neural networks," 2024. 1, 5, 9

[3] J. Bernstein and L. Newhouse, "Modular duality in deep learning," in *International Conference on Machine Learning*, 2025. 1, 24

[4] A. Jacot, F. Gabriel, and C. Hongler, "Neural tangent kernel: Convergence and generalization in neural networks," in *Advances in Neural Information Processing Systems*, 2018. 1, 2, 16
[5] L. Chizat, E. Oyallon, and F. Bach, "On lazy training in differentiable programming," in *Advances in Neural Information Processing Systems*, 2019. 1, 16

[6] G. Yang and E. J. Hu, "Tensor programs iv: Feature learning in infinite-width neural networks," in *International Conference on Machine Learning*, 2021. 1, 2, 16, 17, 18, 19

[7] V. Gupta, T. Koren, and Y. Singer, "Shampoo: Preconditioned stochastic tensor optimization," in *International Conference on Machine Learning*, 2018. 1, 5

[8] D. Carlson, V. Cevher, and L. Carin, "Stochastic spectral descent for restricted boltzmann machines," in *Artificial Intelligence and Statistics*, 2015. 2, 5, 6, 7

[9] A. Mokhtari, H. Hassani, and A. Karbasi, "Stochastic conditional gradient methods: From convex minimization to submodular maximization," *Journal of Machine Learning Research*, 2020. 2, 7, 11

[10] I. Loshchilov and F. Hutter, "Decoupled weight decay regularization," in *International Conference on Learning Representations*, 2019. 3, 5, 9

[11] B. T. Polyak, *Introduction to Optimization*. New York: Optimization Softw., Inc., 1987. 5, 9

[12] J. Duchi, E. Hazan, and Y. Singer, "Adaptive subgradient methods for online learning and stochastic optimization," *Journal of Machine Learning Research*, 2011. 5, 9

[13] H. B. McMahan and M. Streeter, "Adaptive bound optimization for online convex optimization," *Conference on Learning Theory*, 2010. 5, 9

[14] G. Hinton, N. Srivastava, and K. Swersky, "Neural networks for machine learning lecture 6a overview of mini-batch gradient descent," 2012. 5, 9

[15] D. Kingma and J. Ba, "Adam: a method for stochastic optimization," in *International Conference on Learning Representations*, 2014. 5, 9

[16] D. Carlson, Y.-P. Hsieh, E. Collins, L. Carin, and V. Cevher, "Stochastic spectral descent for discrete graphical models," *IEEE Journal of Selected Topics in Signal Processing*, 2016. 5, 6, 7, 9, 14

[17] J. Bernstein, Y.-X. Wang, K. Azizzadenesheli, and A. Anandkumar, "signsgd: Compressed optimisation for non-convex problems," in *International Conference on Machine Learning*. PMLR, 2018, pp. 560–569. 5, 9

[18] Y. Nesterov, "Efficiency of coordinate descent methods on huge-scale optimization problems," *SIAM Journal on Optimization*, vol. 22, no. 2, pp. 341–362, 2012. 6

[19] D. E. Carlson, E. Collins, Y.-P. Hsieh, L. Carin, and V. Cevher, "Preconditioned spectral descent for deep learning," in *Proceedings of the 28th International Conference on Neural Information Processing Systems*, 2015, pp. 2971–2979. 7

[20] M. Frank, P. Wolfe et al., "An algorithm for quadratic programming," *Naval research logistics quarterly*, 1956. 7

[21] K. L. Clarkson, "Coresets, sparse greedy approximation, and the Frank-Wolfe algorithm," *ACM Trans. Algorithms*, 2010. 7

[22] M. Jaggi, "Revisiting frank-wolfe: Projection-free sparse convex optimization," in *International Conference on Machine Learning*, 2013. 7

[23] L. Chen, B. Liu, K. Liang et al., "Lion secretly solves a constrained optimization: As lyapunov predicts," in *International Conference on Learning Representations*, 2024. 8

[24] E. Hazan, K. Levy, and S. Shalev-Shwartz, "Beyond convexity: Stochastic quasi-convex optimization," in *Advances in Neural Information Processing Systems*, 2015. 9

[25] A. Cutkosky and H. Mehta, "Momentum improves normalized SGD," in *International Conference on Machine Learning*, 2020. 9, 11

[26] J. A. Kelner, Y. T. Lee, L. Orecchia, and A. Sidford, "An almost-linear-time algorithm for approximate max flow in
undirected graphs, and its multicombmodity generalizations,” in *ACM-SIAM symposium on Discrete algorithms*. SIAM, 2014. 9

[27] T. T.-K. Lau, Q. Long, and W. Su, “Polargrad: A class of matrix-gradient optimizers from a unifying preconditioning perspective,” *arXiv preprint arXiv:2505.21799*, 2025. 9, 26

[28] J. Liu, J. Su, X. Yao, Z. Jiang, G. Lai, Y. Du, Y. Qin, W. Xu, E. Lu, J. Yan et al., “Muon is scalable for LLM training,” *arXiv preprint arXiv:2502.16982*, 2025. 10

[29] G. Yang, E. Hu, I. Babuschkin, S. Sidor, X. Liu, D. Farhi, N. Ryder, J. Pachocki, W. Chen, and J. Gao, “Tuning large neural networks via zero-shot hyperparameter transfer,” in *Advances in Neural Information Processing Systems*, 2021. 10, 15

[30] A. Yurtsever, S. Sra, and V. Cevher, “Conditional gradient methods via stochastic path-integrated differential estimator,” in *International Conference on Machine Learning*. PMLR, 2019, pp. 7282–7291. 12

[31] M. Zamani and F. Glineur, “Exact convergence rate of the last iterate in subgradient methods,” *SIAM Journal on Optimization*, 2025. 13

[32] A. Defazio, X. Yang, H. Mehta, K. Mishchenko, A. Khaled, and A. Cutkosky, “The road less scheduled,” in *Advances in Neural Information Processing Systems*, 2024. 13

[33] N. Amsel, D. Persson, C. Musco, and R. M. Gower, “The polar express: Optimal matrix sign methods and their application to the muon algorithm,” *arXiv preprint arXiv:2505.16932*, 2025. 14, 26

[34] K. Ahn, B. Xu, N. Abreu, Y. Fan, G. Magakyan, P. Sharma, Z. Zhan, and J. Langford, “Dion: Distributed orthonormalized updates,” *arXiv preprint arXiv:2504.05295*, 2025. 14, 26

[35] J. Hoffmann, S. Borgeaud, A. Mensch, E. Buchatskaya, T. Cai, E. Rutherford, D. de Las Casas, L. A. Hendricks, J. Welbl, A. Clark et al., “Training compute-optimal large language models,” in *Advances in Neural Information Processing Systems*, 2022. 15

[36] L. C. Vankadara, J. Xu, M. Haas, and V. Cevher, “On feature learning in structured state space models,” in *Advances in Neural Information Processing Systems*, 2024. 15, 19

[37] J. Wei, Y. Tay, R. Bommasani, C. Raffel, B. Zoph, S. Borgeaud, D. Yogatama, M. Bosma, D. Zhou, D. Metzler et al., “Emergent abilities of large language models,” *Transactions on Machine Learning Research*, 2022. 15

[38] A. Gu and T. Dao, “Mamba: Linear-time sequence modeling with selective state spaces,” in *First conference on language modeling*, 2024. 15

[39] N. Dey, B. C. Zhang, L. Noci, M. Li, B. Bordelon, S. Bergsma, C. Pehlevan, B. Hanin, and J. Hestness, “Don’t be lazy: CompleteP enables compute-efficient deep transformers,” *arXiv preprint arXiv:2505.01618*, 2025. 16

[40] T. Large, Y. Liu, M. Huh, H. Bahng, P. Isola, and J. Bernstein, “Scalable optimization in the modular norm,” in *Advances in Neural Information Processing Systems*, 2024. 16, 24

[41] G. Yang, J. B. Simon, and J. Bernstein, “A spectral condition for feature learning,” *arXiv preprint arXiv:2310.17813*, 2023. 17, 20, 21

[42] M. Haas, J. Xu, V. Cevher, and L. C. Vankadara, “$\mu P^2$: Effective sharpness aware minimization requires layerwise perturbation scaling,” in *Advances in Neural Information Processing Systems*, 2024. 19

[43] K. He, X. Zhang, S. Ren, and J. Sun, “Delving deep into rectifiers: Surpassing human-level performance on imagenet classification,” in *Proceedings of the IEEE International Conference on Computer Vision (ICCV)*, 2015. 21

[44] T. Flynn, “The duality structure gradient descent algorithm: Analysis and applications to neural networks,” *arXiv:1708.00523*, 2017. 24