One recurring theme in empirical deep learning is that we want to keep intermediate activations at a healthy size – we don’t want them to blow up or vanish.

This principle has inspired a bunch of design choices that are popular in empirical deep learning today – for example initialization schemes (e.g., Xavier and Kaiming He initialization), layer normalization, and gradient clipping.

Muon[1]Muon: An optimizer for hidden layers in neural networks [link]
Keller Jordan, Yuchen Jin, Vlado Boza, Jiacheng You, Franz Cesista, Laker Newhouse, and Jeremy Bernstein. 2024.
is inspired by this same idea, applied to optimization. The goal is to keep activations at every layer reasonably sized (we’ll define “reasonably sized” in a moment); these activations should be reasonably sized not just at initialization, but also throughout training.

Natural norms

We first need to define what we mean when we say that an activation vector vRd\mathbf{v} \in \mathbb{R}^d is “reasonably sized”.

One standard way to measure the size of a vector is the 2\ell_2 norm: v2=i=1dvi2\|\mathbf{v}\|_2 = \sqrt{\sum_{i=1}^d v_i^2}.

However, this norm is dimension-dependent: for example, a vector with all entries ±1\pm 1 has v2=d\|\mathbf{v}\|_2 = \sqrt{d}, which grows with dimension.

The natural norm for vectors: RMS norm

The root mean square (RMS) norm fixes this by averaging over entries instead of summing:

vRMS:=1di=1dvi2=1dv2.\|\mathbf{v}\|_{\text{RMS}} := \sqrt{\frac{1}{d} \sum_{i=1}^d v_i^2} = \frac{1}{\sqrt{d}} \|\mathbf{v}\|_2.

Now a ±1\pm 1 vector has vRMS=1\|\mathbf{v}\|_{\text{RMS}} = 1 regardless of dimension.

Sanity check with a Gaussian random vector
E[vRMS2]=E ⁣[(1dv2)2]=1dE[v22]=1dE ⁣[i=1dvi2]=1di=1dE[vi2]=1dd1=1,\begin{aligned} \mathbb{E}[\|\mathbf{v}\|_{\text{RMS}}^2] &= \mathbb{E}\!\left[\left(\frac{1}{\sqrt{d}} \|\mathbf{v}\|_2\right)^2\right] \\ &= \frac{1}{d}\mathbb{E}[\|\mathbf{v}\|_2^2] \\ &= \frac{1}{d}\mathbb{E}\!\left[\sum_{i=1}^d v_i^2\right] \\ &= \frac{1}{d}\sum_{i=1}^d \mathbb{E}[v_i^2] \\ &= \frac{1}{d} \cdot d \cdot 1 = 1, \end{aligned}

since each viN(0,1)v_i \sim \mathcal{N}(0, 1) has E[vi2]=Var(vi)=1\mathbb{E}[v_i^2] = \text{Var}(v_i) = 1.

Compare with the Euclidean norm:

E[v22]=E ⁣[i=1dvi2]=i=1dE[vi2]=d1=d,\mathbb{E}[\|\mathbf{v}\|_2^2] = \mathbb{E}\!\left[\sum_{i=1}^d v_i^2\right] = \sum_{i=1}^d \mathbb{E}[v_i^2] = d \cdot 1 = d,

which grows with dimension.

Intuitively, the RMS norm measures the scale of a typical entry of v\mathbf{v}, rather than accumulating across all entries like the 2\ell_2 norm does. This makes it the natural norm for dense vectors1: vRMS=Θ(1)\|\mathbf{v}\|_{\text{RMS}} = \Theta(1) exactly captures “the entries of v\mathbf{v} are Θ(1)\Theta(1)-sized”, independent of dimension.

The natural norm for matrices: RMS-to-RMS operator norm

What does it mean for a weight matrix WRdout×din\mathbf{W} \in \mathbb{R}^{d_\text{out} \times d_\text{in}} to be “well-behaved”? It should map reasonably-sized inputs to reasonably-sized outputs: if xRMS=Θ(1)\|\mathbf{x}\|_{\text{RMS}} = \Theta(1), then we want WxRMS=Θ(1)\|\mathbf{W}\mathbf{x}\|_{\text{RMS}} = \Theta(1).

The natural way to measure this is the RMS-to-RMS operator norm – the worst-case RMS stretch:

WRMSRMS:=maxx0WxRMSxRMS.\|\mathbf{W}\|_{\text{RMS} \to \text{RMS}} := \max_{\mathbf{x} \neq \mathbf{0}} \frac{\|\mathbf{W}\mathbf{x}\|_{\text{RMS}}}{\|\mathbf{x}\|_{\text{RMS}}}.

This is the natural spectral norm. It relates to the standard spectral norm W\|\mathbf{W}\|_* (largest singular value) by a dimensional factor:

WRMSRMS=dindoutW.\|\mathbf{W}\|_{\text{RMS} \to \text{RMS}} = \sqrt{\frac{d_\text{in}}{d_\text{out}}} \cdot \|\mathbf{W}\|_*.
Relating RMS-to-RMS and spectral norm
WxRMSxRMS=1doutWx21dinx2=dindoutWx2x2.\frac{\|\mathbf{W}\mathbf{x}\|_{\text{RMS}}}{\|\mathbf{x}\|_{\text{RMS}}} = \frac{\frac{1}{\sqrt{d_\text{out}}} \|\mathbf{W}\mathbf{x}\|_2}{\frac{1}{\sqrt{d_\text{in}}} \|\mathbf{x}\|_2} = \sqrt{\frac{d_\text{in}}{d_\text{out}}} \cdot \frac{\|\mathbf{W}\mathbf{x}\|_2}{\|\mathbf{x}\|_2}.

Taking the max over x\mathbf{x}:

WRMSRMS=dindoutW.\|\mathbf{W}\|_{\text{RMS} \to \text{RMS}} = \sqrt{\frac{d_\text{in}}{d_\text{out}}} \cdot \|\mathbf{W}\|_*.

The spectral scaling condition

Let’s recall what we’re trying to achieve. We want every layer’s weight matrix W\mathbf{W} to map reasonably-sized inputs to reasonably-sized outputs – that is, WRMSRMS=Θ(1)\|\mathbf{W}\|_{\text{RMS} \to \text{RMS}} = \Theta(1). Using the relationship we derived above, this is equivalent to:

W=Θ ⁣(doutdin).\|\mathbf{W}\|_* = \Theta\!\left(\sqrt{\frac{d_\text{out}}{d_\text{in}}}\right).

We want this to hold not just at initialization, but throughout training. After a weight update WW+ΔW\mathbf{W} \to \mathbf{W} + \Delta \mathbf{W}, the change in output for input x\mathbf{x} is Δy=ΔWx\Delta \mathbf{y} = \Delta \mathbf{W} \, \mathbf{x}. By definition of the operator norm:

ΔyRMS=ΔWxRMSΔWRMSRMSxRMS.\|\Delta \mathbf{y}\|_{\text{RMS}} = \|\Delta \mathbf{W} \, \mathbf{x}\|_{\text{RMS}} \leq \|\Delta \mathbf{W}\|_{\text{RMS} \to \text{RMS}} \cdot \|\mathbf{x}\|_{\text{RMS}}.

So if the input is reasonably sized (xRMS=Θ(1)\|\mathbf{x}\|_{\text{RMS}} = \Theta(1)) and we want the output change to also be reasonably sized (ΔyRMS=Θ(1)\|\Delta \mathbf{y}\|_{\text{RMS}} = \Theta(1)), we need ΔWRMSRMS=Θ(1)\|\Delta \mathbf{W}\|_{\text{RMS} \to \text{RMS}} = \Theta(1), i.e.:

ΔW=Θ ⁣(doutdin).\|\Delta \mathbf{W}\|_* = \Theta\!\left(\sqrt{\frac{d_\text{out}}{d_\text{in}}}\right).

Together, these two requirements form the spectral scaling condition from Yang et al. [2]A spectral condition for feature learning [link]
Greg Yang, James B. Simon, and Jeremy Bernstein. arXiv preprint. 2024.
: both the weights and their updates should have Θ(1)\Theta(1) RMS-to-RMS operator norm, or equivalently, spectral norm Θ(dout/din)\Theta(\sqrt{d_\text{out}/d_\text{in}}).

Deriving Muon

We’ve established that ΔW\Delta \mathbf{W} should have Θ(1)\Theta(1) RMS-to-RMS operator norm. That is, ΔWRMSRMSα\|\Delta \mathbf{W}\|_{\text{RMS} \to \text{RMS}} \leq \alpha for some constant α\alpha. Since WRMSRMS=din/doutW\|\mathbf{W}\|_{\text{RMS} \to \text{RMS}} = \sqrt{d_\text{in}/d_\text{out}} \cdot \|\mathbf{W}\|_*, this is equivalent to ΔWη\|\Delta \mathbf{W}\|_* \leq \eta where η:=αdout/din\eta := \alpha \sqrt{d_\text{out}/d_\text{in}}.

Given this budget, we want to decrease the loss as much as possible. To first order, the change in loss is WL,ΔW\langle \nabla_{\mathbf{W}} \mathcal{L}, \Delta \mathbf{W} \rangle, where A,B:=ijAijBij\langle \mathbf{A}, \mathbf{B} \rangle := \sum_{ij} A_{ij} B_{ij} is the entrywise inner product between matrices. So we want to solve:

minΔW  WL,ΔWsubject toΔWη.\min_{\Delta \mathbf{W}} \; \langle \nabla_{\mathbf{W}} \mathcal{L}, \Delta \mathbf{W} \rangle \quad \text{subject to} \quad \|\Delta \mathbf{W}\|_* \leq \eta.

We can write the gradient in its SVD:

WL=i=1rσiuivi,\nabla_{\mathbf{W}} \mathcal{L} = \sum_{i=1}^r \sigma_i \mathbf{u}_i \mathbf{v}_i^\top,

with singular values σ1σ2σr>0\sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_r > 0.

Plugging in the SVD, and using linearity of inner product, we get:

WL,ΔW  =  i=1rσiuivi,  ΔW  =  i=1rσiuivi,ΔW.\langle \nabla_{\mathbf{W}} \mathcal{L}, \Delta \mathbf{W} \rangle \;=\; \left\langle \sum_{i=1}^r \sigma_i \mathbf{u}_i \mathbf{v}_i^\top, \; \Delta \mathbf{W} \right\rangle \;=\; \sum_{i=1}^r \sigma_i \, \langle \mathbf{u}_i \mathbf{v}_i^\top, \Delta \mathbf{W} \rangle.

Note that only the projection of ΔW\Delta \mathbf{W} onto the singular directions {uivi}\{\mathbf{u}_i \mathbf{v}_i^\top\} appears in the objective – any orthogonal component doesn’t help decrease the loss, but can only increase ΔW\|\Delta \mathbf{W}\|_*. Thus, without loss of generality, we can write ΔW=i=1rciuivi\Delta \mathbf{W} = \sum_{i=1}^r c_i \mathbf{u}_i \mathbf{v}_i^\top.

This is an SVD of ΔW\Delta \mathbf{W} with singular values ci\lvert c_i \rvert, so ΔW=maxici\|\Delta \mathbf{W}\|_* = \max_i \lvert c_i \rvert, and the optimization reduces to:

minc1,,cr  i=1rσicisubject tomaxiciη.\min_{c_1, \ldots, c_r} \; \sum_{i=1}^r \sigma_i \, c_i \quad \text{subject to} \quad \max_i \lvert c_i \rvert \leq \eta.

Since all σi>0\sigma_i > 0, the minimum is achieved by setting ci=ηc_i = -\eta for all ii, giving:

ΔW=ηi=1ruivi=ηUV.\Delta \mathbf{W} = -\eta \sum_{i=1}^r \mathbf{u}_i \mathbf{v}_i^\top = -\eta \, \mathbf{U}\mathbf{V}^\top.

This is the Muon update: take the gradient’s SVD, replace all singular values with 1, and then scale by η-\eta.

Comparison with gradient descent

What if we had constrained the Frobenius norm instead of the spectral norm? That is, what if we solved:

minΔW  WL,ΔWsubject toΔWFη.\min_{\Delta \mathbf{W}} \; \langle \nabla_{\mathbf{W}} \mathcal{L}, \Delta \mathbf{W} \rangle \quad \text{subject to} \quad \|\Delta \mathbf{W}\|_F \leq \eta.

If we treat ΔW\Delta \mathbf{W} as a flat vector of entries, the Frobenius norm is just the 2\ell_2 norm, and this is asking: which direction decreases the loss most per unit 2\ell_2 step? That’s the definition of the gradient. So the solution is ΔWWL\Delta \mathbf{W} \propto -\nabla_{\mathbf{W}} \mathcal{L}: gradient descent.

So the key difference between gradient descent and Muon stems from the norm constraint. Gradient descent effectively constrains the Frobenius norm, which treats the weight matrix as an unstructured vector of numbers; Muon effectively constrains the spectral norm, which measures how the matrix acts on inputs.

The spectral norm only constrains the largest singular value of ΔW\Delta \mathbf{W} – once you’ve spent your budget η\eta on the top singular direction, the remaining directions are free. Muon takes advantage of this by stepping equally in every singular direction. Gradient descent can’t do this: under the Frobenius norm, every direction draws from a shared budget (ci2η2\sum c_i^2 \leq \eta^2), so stepping more in one direction means stepping less in another. The optimal allocation sets ciσic_i \propto -\sigma_i, concentrating on the largest singular directions at the expense of the smaller ones.

Making it practical: Newton-Schulz

The (ideal) Muon update requires computing UV\mathbf{U}\mathbf{V}^\top from the gradient WL=UΣV\nabla_{\mathbf{W}} \mathcal{L} = \mathbf{U}\boldsymbol{\Sigma} \mathbf{V}^\top. Computing the full SVD at every step is too expensive. But we don’t need the full SVD – we just need the map UΣVUV\mathbf{U}\boldsymbol{\Sigma} \mathbf{V}^\top \mapsto \mathbf{U}\mathbf{V}^\top.

Turns out there are some algorithms that do exactly this – Newton-Schulz iterations. The algorithm works iteratively: at each step, it applies an operation to the matrix that acts independently on each singular value while preserving the singular vectors, and the operation is chosen so that all singular values converge to 1. In practice, 5-10 iterations suffice, and each iteration is just a few matrix multiplications – fast on GPUs. See these blog posts[3, 1]Deriving Muon [link]
Jeremy Bernstein. 2025.
Muon: An optimizer for hidden layers in neural networks [link]
Keller Jordan, Yuchen Jin, Vlado Boza, Jiacheng You, Franz Cesista, Laker Newhouse, and Jeremy Bernstein. 2024.
for more details.

The Muon update rule

Putting it all together, the Muon update for a weight matrix WRdout×din\mathbf{W} \in \mathbb{R}^{d_\text{out} \times d_\text{in}} is:

WWαdoutdinNewtonSchulz(WL).\mathbf{W} \leftarrow \mathbf{W} - \alpha \cdot \sqrt{\frac{d_\text{out}}{d_\text{in}}} \cdot \text{NewtonSchulz}(\nabla_{\mathbf{W}} \mathcal{L}).

Let’s unpack this:

  • NewtonSchulz(WL)\text{NewtonSchulz}(\nabla_{\mathbf{W}} \mathcal{L}) approximately orthogonalizes the gradient: UΣVUV\mathbf{U}\boldsymbol{\Sigma} \mathbf{V}^\top \mapsto \mathbf{U}\mathbf{V}^\top. This has spectral norm 11.
  • The dout/din\sqrt{d_\text{out}/d_\text{in}} factor scales the RMS-to-RMS learning rate α\alpha to the spectral norm budget: recall η=αdout/din\eta = \alpha \sqrt{d_\text{out}/d_\text{in}}.

Note that the dout/din\sqrt{d_\text{out}/d_\text{in}} factor absorbs the layer dimensions into the update, which helps α\alpha transfer across weight matrices of different shapes and sizes.

References

References are listed in order of citation.

  1. Muon: An optimizer for hidden layers in neural networks [link]
    Keller Jordan, Yuchen Jin, Vlado Boza, Jiacheng You, Franz Cesista, Laker Newhouse, and Jeremy Bernstein. 2024.
  2. A spectral condition for feature learning [link]
    Greg Yang, James B. Simon, and Jeremy Bernstein. arXiv preprint. 2024.
  3. Deriving Muon [link]
    Jeremy Bernstein. 2025.

Footnotes

  1. A dense vector is one where every entry contributes a comparable amount to the squared norm. By contrast, a sparse vector has only a constant number of non-negligible entries, regardless of dimension – e.g., a one-hot encoding vector. For sparse vectors, the ordinary 2\ell_2 norm is already dimension-independent (e.g., a one-hot vector has ei2=1\|\mathbf{e}_i\|_2 = 1, regardless of dimension), so no correction is needed. The RMS norm is the natural norm specifically for dense vectors. Throughout this note we focus on dense vectors, since hidden activations in a transformer can be thought of as dense; the embedding layer (which maps sparse inputs to dense activations) requires separate treatment.