Paper deep dive
HTMuon: Improving Muon via Heavy-Tailed Spectral Correction
Tianyu Pang, Yujie Fang, Zihang Liu, Shenyang Deng, Lei Hsiung, Shuhua Yu, Yaoqing Yang
Intelligence
Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 94%
Last extracted: 3/13/2026, 12:58:29 AM
Summary
HTMuon is a novel matrix-based optimizer for LLM training that improves upon Muon by incorporating Heavy-Tailed Self-Regularization (HT-SR) theory. By applying a power transform (p ∈ (0,1)) to the singular values of the momentum matrix, HTMuon induces heavier-tailed weight spectra, which correlates with better generalization and lower perplexity compared to the orthogonalized update rule of standard Muon.
Entities (5)
Relation Signals (3)
HTMuon → basedon → HT-SR theory
confidence 95% · Motivated by the Heavy-Tailed Self-Regularization (HT-SR) theory, we propose HTMuon.
HTMuon → improvesupon → Muon
confidence 95% · HTMuon consistently improves performance over state-of-the-art baselines and can also serve as a plug-in on top of existing Muon variants.
HTMuon → trains → LLaMA
confidence 90% · For example, on LLaMA pretraining on the C4 dataset, HTMuon reduces perplexity
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Muon has recently shown promising results in LLM training. In this work, we study how to further improve Muon. We argue that Muon's orthogonalized update rule suppresses the emergence of heavy-tailed weight spectra and over-emphasizes the training along noise-dominated directions. Motivated by the Heavy-Tailed Self-Regularization (HT-SR) theory, we propose HTMuon. HTMuon preserves Muon's ability to capture parameter interdependencies while producing heavier-tailed updates and inducing heavier-tailed weight spectra. Experiments on LLM pretraining and image classification show that HTMuon consistently improves performance over state-of-the-art baselines and can also serve as a plug-in on top of existing Muon variants. For example, on LLaMA pretraining on the C4 dataset, HTMuon reduces perplexity by up to $0.98$ compared to Muon. We further theoretically show that HTMuon corresponds to steepest descent under the Schatten-$q$ norm constraint and provide convergence analysis in smooth non-convex settings. The implementation of HTMuon is available at this https URL.
Tags
Links
- Source: https://arxiv.org/abs/2603.10067v1
- Canonical: https://arxiv.org/abs/2603.10067v1
Trouble viewing inline? Open PDF directly →
Full Text
94,320 characters extracted from source content.
Expand or collapse full text
HTMuon: Improving Muon via Heavy-Tailed Spectral Correction Tianyu Pang ∗1 , Yujie Fang ∗2 , Zihang Liu 3,4 , Shenyang Deng 1 , Lei Hsiung 1 , Shuhua Yu 5 , Yaoqing Yang 1 1 Dartmouth College 2 Microsoft 3 International Computer Science Institute 4 University of California, Berkeley 5 Meta Abstract Muon has recently shown promising results in LLM training. In this work, we study how to further improveMuon. We argue thatMuon’s orthogonalized update rule suppresses the emergence of heavy-tailed weight spectra and over-emphasizes the training along noise-dominated directions. Motivated by the Heavy-Tailed Self-Regularization (HT-SR) theory, we proposeHTMuon.HTMuonpreservesMuon’s ability to capture parameter interdependencies while producing heavier-tailed updates and inducing heavier-tailed weight spectra. Experiments on LLM pretraining and image classification show thatHTMuonconsistently improves performance over state-of-the-art baselines and can also serve as a plug-in on top of existing Muonvariants. For example, on LLaMA pretraining on the C4 dataset,HTMuonreduces perplexity by up to 0.98 compared toMuon. We further theoretically show thatHTMuoncorresponds to steepest descent under the Schatten-qnorm constraint and provide convergence analysis in smooth non-convex settings. The implementation of HTMuon is available at https://github.com/TDCSZ327/HTmuon. 1 Introduction Optimizers play a central role in training Large language models (LLMs). A well-designed optimizer can help LLMs learn more effectively from large-scale datasets (Semenov et al., 2025), leading not only to lower training loss but also to improved generalization (Foret et al., 2020; Kaddour et al., 2022). Over the decades, optimizers have often been developed from a vector-based view, where gradients are treated as flattened vectors during updates. Within this framework,Adam(Kingma, 2014) andAdamW(Loshchilov and Hutter, 2017) have become the default choices for training LLMs (Groeneveld et al., 2024; OLMo et al., 2024). A key reason is that they incorporate first- and second-order moment information in an element-wise manner, allowing them to adapt to the different statistical patterns (e.g., local geometry of the loss landscape) across individual coordinates. However, an overly element-wise update scheme often ignores interdependencies among coordinates (e.g., geometric correlations among parameters), which can limit the effectiveness of the optimizer. Although recent work (Liu et al., 2024b; Wang et al., 2025a; Zhang et al., 2024; Zhou et al., 2023) has gone beyondAdamand its variants by using layer-wise or module-wise learning rates to capture parameters’ relationships, whether these adjustments are sufficient to capture the complex coupling across these parameters remains to be seen. Recently,Muon(Jordan et al., 2024) has been proposed as a representative matrix-based optimizer * . Building on earlier matrix-based optimizers such asAdaGrad(Duchi et al., 2011) andShampoo(Gupta et al., 2018),Muonperforms preconditioning on the momentum matrix, and its update can be interpreted as an * Equal contribution. * Matrix-based optimizers use matrix-valued preconditioners, which can better capture geometry-induced dependencies among different parameters. 1 arXiv:2603.10067v1 [cs.LG] 10 Mar 2026 orthogonalization step, which effectively captures geometry interdependencies among parameters. Moreover, the orthogonalization update can be viewed as the steepest descent under the Schatten−∞norm constraint (Bernstein and Newhouse, 2024; Pethick et al., 2025), which leads to a competitive convergence rate and improved stability (Chen et al., 2025; Ma et al., 2026; Shen et al., 2025).Muonhas shown promising results (Liu et al., 2025a; Shah et al., 2025; Wang et al., 2025b; Wen et al., 2025) and has been adopted in large-scale LLM training, including Moonshot’s Kimi K2 (Team et al., 2025) and GLM-4.5 (Zeng et al., 2025). However, some studies have found that the improved performance ofMuonis inversely proportional to the model scale and the number of training steps (Semenov et al., 2025; Wen et al., 2025). The orthogonalization step inMuonupdate sets all singular values of the momentum matrix to one, which means it assigns the same weight to each singular-vector direction of the updates. Although this may helpMuondecrease the loss faster in the early phase of training (Shen et al., 2025), it is well known that directions associated with small singular values tend to be more noise-dominated (Chen et al., 2024; Defilippis et al., 2025; Liu et al., 2025c; Sharma et al., 2023). In Section 3.2, we conduct experiments to show that using uniform weights on all singular vector directions may be suboptimal. More importantly, the orthogonalization step makes the spectrum of the momentum update matrix light-tailed. This, in turn, makes the spectrum of learned weight matrices light-tailed. Here, the spectrum refers to the empirical spectral density (ESD) of the weight matrices. However, a line of work by Martin and Mahoney (2021); Martin et al. (2021) shows that well-trained neural networks—namely, networks that have extracted strong correlations in their weights through learning from data—tend to have heavy-tailed ESDs in their weight matrices (Hodgkinson et al., 2025; Kothapalli et al., 2025). Within a reasonable range, the degree of heavy-tailedness is strongly correlated with model quality (Martin et al., 2021). Motivated by these observations, they proposed the Heavy-Tailed Self-Regularization theory (HT-SR), which has been supported by extensive empirical studies and further theoretical analyses (Defilippis et al., 2025; Martin and Hinrichs, 2025; Yang et al., 2023). In Section 3.3, we empirically show that theMuonoptimizer leads to a less heavy-tailed weight ESD, and it performs worse than variants ofMuonthat induce more heavy-tailed spectra. Based on this, we claim that theMuonupdate rule, which sets all singular values of the momentum matrix to one, limits the final quality that the model can achieve. Motivated by this observation, we propose a matrix-based optimizerHTMuon.HTMuonaims to makeMuon’s momentum update more heavy-tailed, which can lead to more heavy-tailed weight ESDs, while preserving Muon’s advantage in capturing interdependencies among parameters. The design ofHTMuonis simple: we raise the singular values of the momentum matrix to the power ofp, wherep ∈(0,1). See Algorithm 3. Settingp= 1 reduces the method toSGDM(Sutskever et al., 2013), whose updates are independent across parameters. In contrast,p= 0 recoversMuon, which, as discussed earlier, leads to less heavy-tailed updates. Therefore, we choosep∈(0,1) inHTMuonto maintain the matrix-based ability to model parameter coupling, while producing updates that are more heavy-tailed than those ofMuon. Unless otherwise specified, we use p= 0.125 forHTMuonthroughout. We discuss the specific choice ofpin Section 5.7. To evaluate the empirical performance ofHTMuon, we study a range of models and tasks, including training LLaMA models (Touvron et al., 2023) on the C4 dataset, GPT-2 models on the OpenWebText dataset and training image classification models using the ResNet (He et al., 2016) and ViT (Dosovitskiy, 2020) architectures. All of these tasks are standard and have been widely used in prior literature (He et al., 2025; Yuan et al., 2024; Zhao et al., 2024). We compareHTMuonagainst commonly used baseline optimizers such asAdam,Muon,Cautious(Liang et al., 2024),GaLore(Zhao et al., 2024),Sophia(Liu et al., 2024a),Mars(Yuan et al., 2024),SOAP(Vyas et al., 2024),COSMOS(Liu et al., 2025b), and show thatHTMuonnot only achieves superior performance but can also serve as a plug-in module on top of existingMuonvariants likeNorMuon(Li et al., 2025),AdaMuon (Si et al., 2025) to further improve test performance. Moreover, our analysis reveals thatHTMuonproduces update and weight matrices that are more heavy-tailed than those ofMuonduring training, consistent with prior work on HT-SR theory. Finally, we present theoretical results showing thatHTMuonis equivalent to the steepest-descent method under the Schatten−qnorm constraint, generalizingMuon’s equivalence to steepest descent under a Schatten−∞norm constraint, this provides a matrix-version ofpbSGD(Zhou et al., 2020) and also placesHTMuonwithin the linear minimization oracle framework (Bernstein and Newhouse, 2024; Pethick et al., 2025). We also provide a convergence analysis, showing thatHTMuonmatches the sample 2 complexity upper bound of Muon and SGDM in smooth non-convex settings. To summarize, the main contributions of our work are the following: •We argue thatMuon’s orthogonalized update rule biases updates toward noise-dominated directions and suppresses the emergence of heavy-tailed ESDs in the model’s weight matrices, thereby limiting generalization performance according to HT-SR theory. • We proposeHTMuon, an optimizer that makesMuonupdates more heavy-tailed while preservingMuon’s strength in modeling parameter interdependence. We demonstrateHTMuon’s effectiveness in improving performance on both LLM pretraining and image classification tasks relative to state-of-the-art optimizers includingCautious(Liang et al., 2024),Mars(Yuan et al., 2024),SOAP(Vyas et al., 2024),NorMuon (Li et al., 2025),AdaMuon(Si et al., 2025),COSMOS(Liu et al., 2025b), etc. For example,HTMuonreduces perplexity by 0.92 compared toMuonwhen training LLaMA-60M on the C4 dataset, and by 0.98 when training LLaMA-135M on C4. Moreover,HTMuoncan be used as an add-on method in combination with existingMuonvariants to achieve further improvements. We also design two accelerated implementations ofHTMuon. Using these implementations, we show thatHTMuonoutperformsMuonon LLaMA-1B, highlighting its potential for training large-scale models. • We present theoretical results analyzingHTMuon. Specifically, we show that it is equivalent to steepest descent under a Schatten−qnorm constraint, and we provide a convergence analysis demonstrating thatHTMuonmatches the sample-complexity upper bounds ofMuonandSGDMin nonconvex smooth settings. These results supportHTMuon’s competitive convergence rate and improved training stability. 2 Related Work Pretraining Optimizers. A strong optimizer is critical for successful LLM training. For a long time,Adam (Kingma, 2014) andAdamW(Loshchilov and Hutter, 2017) have been the dominant choices. In recent years, a number of methods have emerged that build on theAdamfamily: some focus on reducing gradient variance (e.g.,MARS(Yuan et al., 2024) andCautious(Liang et al., 2024)), while others target memory efficiency (e.g.,Lion(Chen et al., 2023),Adam-mini(Zhang et al., 2024) andGaLore(Zhao et al., 2024)). However, theseAdamvariants remain vector-based optimizers and have difficulty capturing interdependencies among parameters. Some works (Wang et al., 2025a; Zhou et al., 2023) extendAdamby using layer- or module-wise learning rates to capture relationships among parameters at a coarse level. Recently,Muon(Jordan et al., 2024) has been proposed as a representative matrix-based optimizer. It uses the momentum matrix as a preconditioner, which better captures interdependencies among parameters and improves training stability. Muonhas shown promising results in LLM training (Liu et al., 2025a; Shah et al., 2025; Wen et al., 2025). Following this line, severalMuonvariants have emerged, includingNorMuon(Li et al., 2025) andAdaMuon(Si et al., 2025). Heavy-tailed Phenomenon. Heavy tailed weight spectrum in machine learning have been observed and studied in various forms. Among the most prominent empirical findings are a series of works by Martin and Mahoney (2021); Martin et al. (2021), which propose the HT-SR theory for deep neural networks. They observe that well-trained networks often exhibit heavy-tailed spectra in their weight matrices, and the degree of heavy-tailedness is strongly correlated with model quality. Furthermore, they quantify this behavior via power-law fits to the ESDs of weight matrices. From a theoretical perspective, a number of studies (Dandi et al., 2024; Defilippis et al., 2025; Hodgkinson et al., 2022; Kothapalli et al., 2025; Simon et al., 2023; Simsekli et al., 2020) have used tools such as stochastic differential equations, random matrix theory, eigenlearning framework, and approximate message passing to characterize power-law behavior in weight ESDs and its strong connections to generalization bounds and test loss. Empirically, treating HT-SR theory as a diagnostic tool, subsequent work has explored layer-wise hyperparameter scheduling (He et al., 2025; Liu et al., 2024b; Zhou et al., 2023) and layer-wise pruning (Hu et al., 2025; Lu et al., 2024) based on the power-law fit to 3 improve training performance and efficiency. In this work, motivated by HT-SR theory, we useHTMuonto make Muon updates more heavy-tailed to improve training. 3 Motivation for HTMuon In this section, we first briefly introduce HT-SR theory and then use two examples to motivate the importance of heavy-tailed spectra for the design of HTMuon. 3.1 Background on HT-SR Theory The HT-SR theory (Martin and Mahoney, 2021; Martin et al., 2021) offers a principled lens for analyzing the ESD of neural network weight matrices. Empirically, well-trained models tend to exhibit more heavy-tailed ESDs, and the degree of heavy-tailedness correlates with training quality. Based on prior work (Martin and Mahoney, 2021; Zhou et al., 2023), heavy-tailedness is quantified by fitting a power law (PL) to the ESD and using the resulting PL exponent α as the metric. Given a network withLlayers and weight matricesW l L l=1 of shapen× m, we compute the ESD by extracting the eigenvalues of the correlation matrixX l =W ⊤ l W l for each module. We fit a power law (PL) to the ESD in the formp(λ)∝ λ −α , λ min < λ < λ max ,wherep(λ) denotes the eigenvalue density within the specified range. The PL exponentαserves as a proxy for the degree of heavy-tailedness, smallerαindicates a more heavy-tailed weight ESD, and more heavy-tailed spectra are often associated with better model quality. 3.2 Unit Singular Values May Be Suboptimal In this subsection, we introduce an intriguing empirical finding: forMuonimplementation, the numerically accelerated version by Newton Schulz, which we callMuonNS(see Algorithm 1) can outperform the theoretically exact implementation MuonSVD (see Algorithm 2). LLaMA-60MLLaMA-135M 20.0 22.5 25.0 27.5 30.0 Perplexity ( ↓ ) Muon_NS Muon_SVD (a) MuonNS vs MuonSVD 10 2 10 1 10 0 Normalized Singular Values 10 0 10 1 10 2 Count Step 1 Step 9000 Step 19000 (b) MuonNS Update Spectrum 10 3 10 2 10 1 10 0 Normalized Singular Values 10 2 10 3 Step 1 Step 9000 Step 19000 (c) MuonSVD Update Spectrum Figure 1:MuonNSvs.MuonSVDon C4 dataset. (a) Validation perplexity for LLaMA-60M/135M trained: MuonNSconsistently achieves lower perplexity thanMuonSVD. Both Learning rates for 60M is 0.03 and for 135M is 0.02. (b)(c) Spectra of update matrices at steps 1/9000/19000 shown in different colors:MuonSVD enforces an exactly ”all-ones” spectrum in the update matrices;MuonNSstays close to one but retains noticeable deviations, implicitly down-weighting noise-dominated singular-vector directions and correlating with improved performance. The theoreticalMuonupdate in Algorithm 2 has a similar form to polar decomposition:O t =U t V ⊤ t = (M t M ⊤ t ) − 1 2 M t . This is equivalent to setting all singular values of the momentum matrixM t to one. In practice, however, computingU t andV t exactly viaSVDis expensive, soMuonuses aNewtonSchulziteration (Higham and Schreiber, 1990) to approximateO t . In Figure 1, we train LLaMA-60M and LLaMA-135M on C4 dataset usingMuonSVDandMuonNS, with identical hyperparameters for each model. Figure 1c 4 visualizes the singular-value distributions of theMuonupdate matrix at different training steps. Although theNewton-Schulziteration closely approximates the ”all-ones” singular spectrum, many singular values still deviate noticeably from one. Meanwhile, as shown in Figure 1a,MuonNSoutperformsMuonSVDon both LLaMA-60M and LLaMA-135M. This suggests that keeping the same weight for all singular-vector directions throughout training may not be the most helpful. In particular, directions associated with smaller singular values are often more noise-dominated (Chen et al., 2024; Defilippis et al., 2025; Liu et al., 2025c; Sharma et al., 2023; Wang et al., 2023); enforcing exactly equal weights across all directions yields a light-tailed update, which can make late-stage training more sensitive to noise and limit model capacity. In contrast, due to the polynomial used in theNewton-Schulzsteps,MuonNSimplicitly assigns smaller weights to noise-dominated singular-vector directions, which improves performance. This observation also supports the motivation for making Muon updates more heavy-tailed. Algorithm 1 Muon 1:Input: Initial weights W 0 ∈ R m×n , loss functionL, learning rateη, momentum parameterβ, weight decay λ. 2: Initialize M 0 ∈ R m×n ← 0 3: for t = 1, 2,... do 4:G t ←∇ W L(W t ) 5:M t ← βM t−1 + (1− β)G t 6: O t ← NewtonSchulz5(M t ) 7: s← q max 1, m n 8:W t+1 ← W t − ηλW t − ηsO t 9: end for Algorithm 2 MuonSVD 1:Input: Initial weights W 0 ∈ R m×n , loss functionL, learning rateη, momentum parameterβ, weight decay λ. 2: Initialize M 0 ∈ R m×n ← 0 3: for t = 1, 2,... do 4:G t ←∇ W L(W t ) 5:M t ← βM t−1 + (1− β)G t 6: U t , Σ t ,V ⊤ t ← SVD(M t ) 7: O t ← U t V ⊤ t 8: s← q max 1, m n 9:W t+1 ← W t − ηλW t − ηsO t 10: end for 3.3 From Updates to Weights: Muon Produces Lighter-Tailed Weight Spectra In this section, we show thatMuon’s orthogonalization update makes weight matrices ESDs less heavy-tailed, which, based on HT-SR theory, may limit the model performance. We train LLaMA-60M and LLaMA-135M on the C4 dataset usingMuonNSandCOSMOS, a variant ofMuon (Liu et al., 2025b). In Figure 2, we visualize the average fitted PL exponent ̄αfor the trained LLaMA models across layers under the two optimizers. We find that ̄αacross layers is higher forMuon, indicating thatMuon makes the weight ESD less heavy-tailed. We further observe thatMuonyields higher perplexity (PPL) than COSMOSon both models. Based on HT-SR theory, this indicates thatMuon’s orthogonalized update rule may limit the model’s final quality, which motivates us to makeMuonupdate more heavy-tailed, leading to more heavy-tailed weight matrices. 5 LLaMA-60MLLaMA-135M 4.0 4.2 4.4 4.6 4.8 5.0 ̄ α Muon COSMOS (a) Average PL ̄α LLaMA-60MLLaMA-135M 22 24 26 28 30 Perplexity ( ↓ ) Muon COSMOS (b) Perplexity Figure 2: (a)Average PL ̄αof weight ESDs for LLaMA-60M and LLaMA-135M trained on C4 withMuon andCOSMOS.Muonyields a higher mean ̄α, indicating less heavy-tailed spectra thanCOSMOS; (b)COSMOS outperforms Muon for LLaMA-60M and LLaMA-135M models on C4 datatset. 4 Methodology In this section, we introduce our methodHTMuon. Our design goal is to preserveMuon’s ability to capture parameter interdependencies while making its updates more heavy-tailed, which in turn leads to more heavy-tailed weight ESDs. We putHTMuonin Algorithm 3.HTMuonis similar toMuon, and the only difference is that we consider a different power transform of the momentum matrix’s singularpin line 7. It is easy to know that whenp= 1,HTMuonreduces toSGDM.SGDMcan be viewed as a vector-based optimizer and has limited ability to capture interdependencies among parameters. Whenp= 0,HTMuonreduces toMuon; as discussed in Section 3.3, this update makes weight-matrix ESDs less heavy-tailed. We therefore proposeHTMuonwithp ∈(0,1). We believeHTMuonhelps mitigate the design limitations ofMuonandSGDM: On the one hand, in this regime,HTMuonremains a matrix-based optimizer and thus retains the ability to capture parameter interdependencies. On the other hand,HTMuonupdates are more heavy-tailed thanMuonupdates, and we illustrate this with a simple example: suppose that during training, the singular values of the momentum matrixM t follows k =s 1 k −s , wheres 1 is the largest singular value of M t . By Lemma 4.1, the PL exponentαof theHTMuonupdate matrixO t is 1 + 1 2sp , it is easy to find that whenp→0, theαwill increase, this meansHTMuonupdate is more heavy-tailed thanMuon’s, which will in turn lead to more heavy-tailed weight ESDs. Lemma 4.1 (Proof in Appendix A.3). Suppose the singular values of matrixW ∈ R n×m follows k = s 1 k −s , 1≤ k ≤ m, we have PL exponent α of W satisfies α = 1 + 1 2s . 5 Experiments In this section, we evaluateHTMuonon various pretraining tasks, and compare it with state-of-the-art pretraining optimizers, includingAdam(Kingma, 2014),AdamW(Loshchilov and Hutter, 2017),SGDM(Sutskever et al., 2013),NorMuon(Li et al., 2025),AdaMuon(Si et al., 2025),Cautious(Liang et al., 2024),GaLore(Zhao et al., 2024),Sophia(Liu et al., 2024a),Mars(Yuan et al., 2024),SOAP(Vyas et al., 2024),COSMOS(Liu et al., 2025b). We provide detailed definitions of these optimizers in Appendix D. Unless otherwise specified, allMuonresults in our experiments use theMuonNSimplementation. We demonstrate thatHTMuonachieves superior performance than these optimizers across all datasets. In addition, we show thatHTMuoncan be used as an add-on method to combine with existingMuonvariant optimizers to achieve further improvements. Furthermore, to reduce the overhead ofSVDinHTMuon, we design two acceleration schemes that significantly cut the computational cost while still outperformingMuon. To ensure a fair comparison, we reproduce all results of the baseline methods with the codebases provided in previous papers. 6 Algorithm 3 HTMuon 1:Input: Initial weights W 0 ∈ R m×n , loss functionL, learning rateη, momentum parameterβ, weight decay λ, power p∈ (0, 1). 2: Initialize M 0 ∈ R m×n ← 0 3: for t = 1, 2,... do 4:G t ←∇ W L(W t ) 5:M t ← βM t−1 + (1− β)G t 6: U t , Σ t ,V ⊤ t ← SVD(M t ) 7: O t ← U t Σ p t V ⊤ t 8: s← q max 1, m n 9:W t+1 ← W t − ηλW t − ηsO t 10: end for Table 1: Comparison with dominant Pretraining Optimizers on LLaMA Models of Varying Sizes on the C4 dataset. Lower perplexity indicates better performance. Detailed hyperparameters are provides in Table 18, 19 and 20 in Appendix C. LLaMa-60M LLaMa-135M LLaMa-350M Adam32.2123.0117.11 AdamW31.8523.3316.96 Muon28.8022.2316.81 HTMuon27.8821.2516.79 5.1 Experimental Setup Tasks. We conduct experiments on various pretraining tasks, including: 1) LLM Pretraining, in which we pretrain LLaMA-family models on the C4 dataset (Raffel et al., 2020), GPT-2 family models on the OpenWebText dataset (Gokaslan and Cohen, 2019). 2) Image Classification, where we train ResNet models on the CIFAR-100 and CIFAR-10 datasets (Krizhevsky and Hinton, 2009) and train ViT-tiny on the ImageNet-1K dataset (Deng et al., 2009). Models. For LLM pretraining, we evaluateHTMuonon LLaMA-60M, LLaMA-135M, LLaMA-350M, LLaMA- 1B (Touvron et al., 2023) and GPT-2 small (125M) (Radford et al., 2019). For image classification, we use ResNet18, ResNet50 (He et al., 2016) and ViT-tiny (Dosovitskiy, 2020). Unless otherwise specified, we usep= 0.125 forHTMuonthroughout. We discuss the specific choice of pin Section 5.7. All baseline methods andHTMuonare carefully tuned for a fair comparison. For detailed hyperparameter settings, please refer to Appendix C. 5.2 LLM Pretraining In this section, we evaluate HTMuon against several baseline optimizers on C4 and OpenWebText datasets. HTMuon consistently outperformsMuon. In Table 1, We find thatHTMuonconsistently outperformsMuon across three model scales (60M, 135M, and 350M) on C4 dataset. In particular, on LLaMA-60M,HTMuon achieves a PPL that is 0.92 lower thanMuonand 4.33 lower thanAdam; on LLaMA-135M,HTMuonis 0.98 lower thanMuonand 2.08 lower thanAdamW. In Table 2, we show that on GPT-2 small,HTMuonachieves a PPL that is 0.26 lower thanMuonand 2.99 lower thanAdamW. In Figure 9 in Appendix B, we provide the training loss curves ofHTMuonandMuon. These results suggest thatHTMuon’s update rule effectively mitigates the limitations of Muon updates. 7 Table 2: Comparison with dominant Pretraining Optimizers on GPT-2 small on the OpenWebText dataset. Lower perplexity indicates better performance. Hyperparameter settings are provided in Appendix C. AdamW Muon HTMuon GPT-2 small25.1922.4622.20 HTMuonconsistently outperforms other variants ofMuon. In Figure 3, we find thatHTMuonconsistently outperformsMuonvariant optimizers likeAdaMuonandNorMuon. Specifically, on LLaMA-60M, our PPL is 0.29 lower than that ofNorMuon, the strongest baseline afterHTMuon; on LLaMA-135M, we outperformNorMuon by 0.74 PPL andAdaMuonby 1.18 PPL. Moreover, we observe that combiningHTMuonwithNorMuonyields further gains, reducing PPL by an additional 0.33 on LLaMA-60M and 0.14 on LLaMA-135M. MuonHTMuonAdaMuonNorMuon HTMuon+NorMuon 27.0 27.5 28.0 28.5 29.0 Perplexity ( ↓ ) (a) LLaMA-60M 21.0 21.5 22.0 22.5 (b) LLaMA-135M Figure 3: Comparison withMuonvariant optimizers on LLaMA-60M and 135M on C4 dataset. All optimizers are carefully tuned via grid search; detailed results and hyperparameter settings are provided in Table 7 and 19 in Appendix B and C. HTMuonconsistently outperforms state-of-the-art pretraining optimizers. In Figure 4, we compare HTMuonwith state-of-the-art optimizers on LLaMA-60M and LLaMA-135M.HTMuonconsistently outperforms all competing methods at both scales. For example, compared toCOSMOS(the second-best method),HTMuon effectively reduces PPL by 1.07 on LLaMA-60M and by 1.04 on LLaMA-135M. In Appendix B, we include more recent optimizers, including GaLore and Sophia. MuonHTMuonCautiousMARSSOAPCOSMOS 27.5 28.0 28.5 29.0 29.5 30.0 Perplexity ( ↓ ) (a) LLaMA-60M 21.0 21.5 22.0 22.5 23.0 23.5 (b) LLaMA-135M Figure 4: Comparison with state-of-the-art pretraining optimizers on LLaMA-60M and 135M on C4 dataset. All optimizers are carefully tuned via grid search; detailed results and hyperparameter settings are provided in Table 7 and 19 in Appendix B and C. 8 Table 3: Results on CIFAR-100 and CIFAR-10 dataset. Higher accuracy indicates better performance. CIFAR-100CIFAR-10 OptimizerResNet18ResNet50 ResNet18 ResNet50 SGDM77.8278.5795.1595.16 Muon77.9579.8595.3995.97 NorMuon77.8279.7895.5796.03 HTMuon78.2480.1695.6396.13 HTMuon + NorMuon78.5180.2295.4096.04 Table 4: Results on ImageNet-1K. Hyperparameter settings are provided in Appendix C. Adam Muon HTMuon ViT-Tiny67.6771.0271.16 5.3 Image Classification In this section, we compareHTMuonagainst several baseline optimizers on CIFAR and ImageNet-1K datasets. In Table 3, we find that, when training ResNet18 and ResNet50 on CIFAR-10 and CIFAR-100,HTMuon consistently achieves higher accuracy thanSGDM,Muon, andNorMuon. For example, relative toMuon,HTMuon improves accuracy by up to 0.31 % on CIFAR-100 and up to 0.24 % on CIFAR-10. Moreover, on CIFAR-100, combiningHTMuonwithNorMuonfurther improves the accuracy of both ResNet18 and ResNet50. In Table 4, we show that HTMuon also outperforms Muon and Adam when training ViT-tiny on ImageNet-1K. 5.4 More Efficient Implementations of HTMuon In this section, to reduce the overhead introduced bySVDinHTMuon, we consider two acceleration strategies: (i) ApplyHTMuonevery fixed number of training steps and use Muon otherwise. In Figure 5a and 5c, we evaluate applyingHTMuonon LLaMA-60M and LLaMA-135M every 5, 10, and 25 steps. Compared to usingHTMuonat every step, this strategy greatly reduces runtime overhead, making it competitive with Muon, while still outperformingMuon. We also observe that smaller intervals lead to better performance, further supporting the effectiveness of our optimizer. Based on these findings, we train LLaMA-1B on C4 usingHTMuonwith interval = 5. As shown in Table 5, HTMuon still outperforms Muon, highlighting its potential for large-scale training. (i) Replace theSVDwith a numerical iterative method for faster computation. In Algorithm 3, We note that theHTMuonupdateO t =U t Σ p t V ⊤ t admits the factorizationO t = (U t V ⊤ t ) (V t Σ p t V ⊤ t ).For theU t V ⊤ t part, this isMuon’s update; we can considerNewton-Schulzwith 5 steps. To computeV t Σ p t V ⊤ t efficiently, we apply theNSrootroutine toM ⊤ t M t in Algorithm 5, sinceM ⊤ t M t =V t Σ 2 t V ⊤ t , we have V t Σ p t V ⊤ t = (M ⊤ t M t ) p 2 . Accordingly,NSrootusesNewton-Schulziterations to approximate matrix square roots and applies successive square-root operations for [1− log 2 p] rounds. Combining these two, we obtain HTMuonNSin Algorithm 4. In Figure 5b and 5d, comparingHTMuonNSwithHTMuonwith interval= 1, we findHTMuonNSsubstantially reduces runtime overhead, while only slightly increases PPL. We also apply HTMuonNSon LLaMA-60M and LLaMA-135M every 5, 10, and 25 steps. We find that this further reduces runtime overhead while still maintaining better performance thanMuon, for example,HTMuonNS(interval = 5) achieves 0.59 s/step compared toMuon’s 0.51 s/step in Table 9(which only introduces a minor increase in overhead), Table 8 shows that HTMuonNS (interval = 5) still outperforms Muon by 0.38 PPL. We also evaluateHTMuon+NorMuonwith different intervals. As observed earlier, using larger intervals reduces runtime overhead while still outperforming NorMuon. Please refer to Figure 10 in Appendix B. 9 Algorithm 4 HTMuonNS 1:Input: Initial weights W 0 ∈ R m×n , loss functionL, learning rateη, momentum parameterβ, weight decay λ, power p. 2: Initialize M 0 ∈ R m×n ← 0 3: for t = 1, 2,... do 4:G t ←∇ W L(W t ) 5:M t ← βM t−1 + (1− β)G t 6: O ′ t ← NS5(M t ) 7: O ′ t ← NS 2proot(M ⊤ t M t ) 8: O t ← O ′ t O ′ t 9: s← q max 1, m n 10:W t+1 ← W t − ηλW t − ηsO t 11: end for Algorithm 5 NSroot 1: Input: Positive definite matrix X n×n , small constant ε, power p, nssteps T . 2: for t = 1, 2,..., [log 2 2 p ] do 3: α←∥X∥ F 4: X ← X ∥X∥ F +ε 5: Y ← X,Z ← I 6:for i = 1, 2,...,T do 7:Q← 3.0I −ZY 8:Y ← 0.5Y Q 9:Z ← 0.5QZ 10:end for 11: X ← √ αY , X ← X+X ⊤ 2 + εI 12: end for Table 5: We train LLaMA-1B with interval= 5 on C4. Hyperparameter settings are provided in Table 20 in Appendix C. Adam AdamW Muon HTMuon LLaMA-1B15.2215.1114.3314.17 Besides reporting per-step runtime overhead, we also report the wall-clock time forMuon,HTMuonand HTMuonNSwith intervals. To better report wall-clock time, we rerun training on NVIDIA RTX PRO 6000 GPUs: LLaMA-60M on C4 using 2 GPUs, and LLaMA-135M on C4 using 4 GPUs. We use the optimal hyperparameter settings reported in the Table 19. We summarize the PPL and wall-clock time ofMuon, HTMuon, HTMuonNS, HTMuon (Interval = 5), and HTMuonNS (Interval = 5) in the Table 10 in Appendix B. We observe that for LLaMA-60M,HTMuonNS(Interval = 5) achieves only a∼6% additional overhead yet outperformsMuonby 0.34 PPL. Moreover, the time forHTMuonNS(Interval = 5) to reach 28.84 PPL is 0.64 hour, compared to 0.67 hour forMuon. For LLaMA-135M,HTMuonNS(Interval = 5) incurs only∼11% additional overhead while outperformingMuonby 0.25 PPL; it reaches 22.27 PPL in 1.32 hour, compared to 1.40 hour for Muon. We also give the FLOPs analysis for Muon and HTMuonNS in Appendix B. We find that althoughHTMuoncould take longer thanMuon, the substantial gains achieved through heavy- tailed spectral correction suggest thatHTMuoncaptures the correct inductive bias for learning. Furthermore, we believe that efficient variants such asHTMuonNSwith interval-based updates provide a strong balance between performance and efficiency. 10 MuonInterval=1Interval=5Interval=10Interval=25Time(sec) 27.0 27.5 28.0 28.5 29.0 Perplexity ( ↓ ) 0.0 0.5 1.0 1.5 (a) HTMuon LLaMA-60M 27.0 27.5 28.0 28.5 29.0 0.0 0.5 1.0 1.5 (b) HTMuonNS LLaMA-60M 21.0 21.5 22.0 22.5 0.0 0.5 1.0 1.5 2.0 (c) HTMuon LLaMA-135M 21.0 21.5 22.0 22.5 0.0 0.5 1.0 1.5 2.0 Time (sec) (d) HTMuonNS LLaMA-135M Figure 5: we evaluate applyingHTMuonandHTMuonNSon LLaMA-60M and LLaMA-135M every 1, 5, 10, and 25 steps. We report the average per-step runtime overhead for all methods. Detailed results and hyperparameter settings are provided in Table 8 and 9 in Appendix B. 5.5Weight Spectrum Analysis HTMuonsuccessfully produces more heavy-tailed weight matrices. In Figure 6, we fit the layer-wise PLαfor models trained withMuonandHTMuon. Across both LLaMA and ResNet models, we find that HTMuonyields smallerαthanMuonfor most layers, leading to a lower average ̄α. Meanwhile,HTMuonachieves better performance, consistent with the HT-SR observation that better-trained models tend to have lowerα. For more PL α plots, please refer to Figure 11a and 11b in Appendix B. 11020304050 Layer index 5 10 15 α Muon ( ̄ α=4.786) HTMuon ( ̄ α=4.211) (a) LLaMA-60M on C4 11020304050607080 Layer index 5 10 15 α Muon ( ̄ α=4.750) HTMuon ( ̄ α=4.347) (b) LLaMA-135M on C4 15101520253035404550 Layer index 2 3 4 α Muon ( ̄ α=2.488) HTMuon ( ̄ α=2.295) (c) ResNet50 on CIFAR-100 15101520253035404550 Layer index 2 3 4 5 α Muon ( ̄ α=2.777) HTMuon ( ̄ α=2.364) (d) ResNet50 on CIFAR-10 Figure 6: Layer-wise PLαfor LLaMA and ResNet model weights trained withMuonandHTMuon. All models used for visualization are trained using each optimizer’s best-performing hyperparameter configuration. For hyperparameter configurations, please refer to Appendix C. Additional generalization metrics. In Figure 7, we visualize the layer-wise spectral norm and Frobenius norm for LLaMA models trained withMuonandHTMuon. We find that, for these models,HTMuonconsistently yields smaller norms thanMuon. This is consistent with prior findings (Miyato et al., 2018) that smaller spectral and Frobenius norms are often associated with better generalization. 5.6 Downstream Tasks To further evaluateHTMuon’s applicability to broader NLP scenarios, in this section, we evaluate LLaMA-1B model trained withHTMuonon 7 commonsense tasks using the lmevalharness framework. For all tasks we use the default prompt and perform zero-shot evaluation. As shown in the Table 6,HTMuonachieves the best test accuracy for 4 out of 7 commonsense tasks, and achieves the best average score across 7 tasks, significantly outperforming the second best optimizer (Muon) by 1.05. This demonstrates the effectiveness of HTMuon on wider NLP scenarios. 11 11020304050 Layer index 5 10 Spectral Norm Muon (Spectral Norm=7.091) HTMuon (Spectral Norm=4.037) (a) LLaMA-60M on C4 11020304050607080 Layer index 5 10 15 Spectral Norm Muon (Spectral Norm=9.798) HTMuon (Spectral Norm=4.903) (b) LLaMA-135M on C4 11020304050 Layer index 20 40 60 80 100 Frobenius Norm Muon (F_Norm=61.869) HTMuon (F_Norm=35.092) (c) LLaMA-60M on C4 11020304050607080 Layer index 50 100 150 Frobenius Norm Muon (F_Norm=100.665) HTMuon (F_Norm=50.382) (d) LLaMA-135M on C4 Figure 7: Layer-wise spectral norm and frobenius norm for LLaMA model weights trained withMuonand HTMuon. All models used for visualization are trained using each optimizer’s best-performing hyperparameter configuration. For hyperparameter configurations, please refer to Appendix C. Table 6: Zero-shot evaluation results (↑) on seven commonsense reasoning benchmarks for the LLaMA-1B model pretrained with different methods. Optimizer ARC-c ARC-e PIQA Hellaswag OBQA Winogrande BOOLQAvg. Adam21.9331.2765.1827.4917.4053.8362.1739.90 AdamW21.1632.3764.2527.2917.8052.4962.2339.66 Muon20.5632.6665.1331.5317.0051.2262.1440.03 HTMuon22.1833.1666.5933.7617.2052.7261.9341.08 5.7 Ablation Study Varying differentp. Sincepis an additional hyperparameter introduced byHTMuon, we perform a grid search overpon LLaMA models in Figure 8. We find thatp= 0.125 is a good choice, which is why we use p= 0.125 in most experiments. In Appendix B.4.2, we discuss more about the hyperparameter sensitivity for p, we show that although introducingpadds a tuning dimension and add hyperparameter search burden, p= 0.125 generalizes well across tasks and architectures and often can be optimal or near-optimal across tasks and architectures. We suggest that p = 0.125 can serve as a recommended default for other tasks. HTMuon HTMuon+NorMuon 00.10.1250.20.25 p 27 28 29 30 Perplexity ( ↓ ) (a) LLaMA-60M 00.10.1250.20.25 p 21.0 21.5 22.0 22.5 (b) LLaMA-135M Figure 8: We conduct grid search onp. We findp= 0.125 is a strong choice. Note thatp= 0 reduces to Muon. Detailed results provided in Table 14 in Appendix C. Varying learning rates. In Figure 12, we visualize performance across learning-rate grids for LLaMA models on C4 and ResNet models on CIFAR datasets. We find that, over most learning rates,HTMuon consistently outperformsMuonandNorMuon, demonstrating the stability ofHTMuon. For more details, please refer to Table 11, 12 and 13 in Appendix B. 12 Algorithm 6 HTMuonHT 1:Input: Initial weights W 0 ∈ R m×n , loss functionL, learning rateη, momentum parameterβ, weight decay λ, power α. 2: Initialize M 0 ∈ R m×n ← 0 3: for t = 1, 2,... do 4:G t ←∇ W L(W t ) 5:M t ← βM t−1 + (1− β)G t 6: U t , Σ t ,V t ⊤ ← SVD(M t ) 7:Σ t i ← i −α 8: O t ← U t Σ t V t ⊤ 9: s← q max 1, m n 10:W t+1 ← W t − ηλW t − ηsO t 11: end for Varying different ways to makeMuonupdate more heavy-tailed. To better prove makingMuon update more heavy-tailed can improve performance, we design another optimizerHTMuonHTin Algorithm 6. Specifically, at each step,HTMuonHTdirectly replaces the singular values of the momentum matrix with a heavy-tailed sequencei −α m i=1 . As shown in Table 17, althoughHTMuonHTis not as strong asHTMuon, it still outperformsMuonon LLaMA models. This shows that merely inducing a heavy-tailed spectrum already yields non-trival gains, and our algorithmHTMuonfurther show that heavy tail spectral correction on top of preserved spectral information is more effective, we put more discussion about the performance improvement attribution in Appendix B.4.3. In Figure 11c and 11d, we visualize the layer-wise PLαforHTMuonHT, and as intended, it yields lower α than Muon across layers. 6 Theoretical Analysis In this section, we provide theoretical understandings of HTMuon. Notation. For a vectorv, we denote itsl 2 norm as∥v∥ 2 . For a matrixA,B ∈ R m×n , we denote its spectral norm as∥A∥, Frobenius norm as∥A∥ F , nuclear norm as∥A∥ ∗ , where∥A∥ ∗ = P r i=1 σ i . Here,σ i r i=1 are singular values. Forp >1, we denote Schatten−pnorm as∥A∥ p , where∥A∥ p =( P r i=1 σ p i ) 1 p . We use⟨A,B⟩ to denote the inner product between A and B, i.e., ⟨A,B⟩ = Tr A ⊤ B . 6.1 Update Under the Schatten−q Norm Constraint We are here to provide a view of Schatten norm constrained optimization forHTMuon. To determine the descent direction at step t, we formulate the update as the following constrained minimization problem: min ∆W∈R m×n f (W t + ∆W )(1) subject to ∥∆W∥ q ≤ δ where∥·∥ q denotes the Schatten−qnorm. When we employ the first-order approximation for our optimization objective f (W t + ∆W )≈ f (W t ) + Tr(∇f (W t ) ⊤ ∆W )(2) It is straightforward to derive that optimizing this approximation is equivalent to the following objective: max ∆W Tr( f G t ⊤ ∆W )(3) subject to ∥∆W∥ q ≤ δ, 13 where f G t =−∇f (W t ) denotes the negative gradient. Theorem 6.1. Let f G t =UΣV ⊤ be theSVDof the negative gradient f G t =−∇f(W t ), where Σ = diag(σ 1 ,...,σ r ). Letpbe the conjugate exponent ofq, satisfying 1/p ′ + 1/q= 1. The explicit solution ∆W ∗ to the optimization problem 3 is given by: ∆W ∗ = δ 1 ∥ f G t ∥ p ′ /q p ′ U Σ p ′ −1 V ⊤ ,(4) where Σ p ′ −1 = diag(σ p ′ −1 1 ,...,σ p ′ −1 r ). We put proof in Appendix A.1. Based on Theorem 6.1, whenq ∈(2,∞),p ′ ∈(1,2), we can assume p=p ′ −1 and find the ∆W ∗ is exactly the update byHTMuon, which meansHTMuonis equivalent to steepest descent under a Schatten−qnorm constraint. It is well known thatMuoncan be viewed as the steepest descent under the Schatten−∞norm constraint (Bernstein and Newhouse, 2024; Pethick et al., 2025), so HTMuongeneralizesMuon’s equivalence to steepest descent under a Schatten−∞norm constraint and this will also benefit training (Davis and Drusvyatskiy, 2025). 6.2 Convergence Analysis Here we give the convergence analysis under a smooth non-convex setting for findingε-stationary points, i.e. ∥∇f (W )∥ ∗ ≤ ε. We consider the following stochastic optimization problem: min W∈R m×n f (W ) = E ξ [f (W ;ξ)],(5) where ξ is stochastic noise. We denote f ∗ = inf W f (W ), r = minm,n. We assume f ∗ >−∞. For ease of theoretical analysis, we rewrite HTMuon in a stochastic setting and present it in Algorithm 7. We assumefis anLFrobenius norm Lipchitz smooth function. We assume that∇f(W;ξ) is an unbiased stochastic estimator of the true gradient∇f(W) and has a bounded variance. We also assume thatM t ’s singular values are bounded by l. We then give assumptions below, some assumptions and definitions are adapted from Shen et al. (2025). Assumption 6.2 (F-norm Lipschitz smooth ). Definef:R m×n → RisLFrobenius norm Lipschitz smooth if for any W,W ′ ∈ R m×n , we have ∥∇f (W )−∇f (W ′ )∥ F ≤ L∥W −W ′ ∥ F . Assumption 6.3 (Bounded variance). We assume∇f(W;ξ) is an unbiased stochastic estimator of the true gradient ∇f (W ) and has a bounded variance, i.e. E[∇f (W ;ξ)] =∇f (W ), E∥∇f (W ;ξ)−∇f (W )∥ 2 F ]≤ σ 2 . Assumption 6.4 (Bounded singular value). We assume the largest singular value σ max (M t )≤ l. Definition 6.5. We say W is an ε-nuclear norm stationary point of f if ∥∇f (W )∥ ∗ ≤ ε. Theorem 6.6 (HTMuon). Under Assumptions 6.2, 6.3 and 6.4, if we applyHTMuonin Algorithm 7 with adaptive learning rateη t = ⟨M t ,ρ(M t )⟩ L∥ρ(M t )∥ 2 F ,ρ(M t ) =U t Σ p t V ⊤ t ,η=maxη t T t=1 , batch sizeB=Tandβ= 1 √ T , ∆ = f (W 0 )− f ∗ , we have 1 T T X t=1 E ∥∇f (W t )∥ 2 1+p ≤ 2Lr 1 s 1− ε ′ 1− ε ′ L + 1 Lε ′ 2σ T + 4σ 2 T 3 + 4rl 2 q L 2 η 2 T + 2Lr 1 s 1− ε ′ ∆ T , where s = 1+p 1−p and 0 < ε ′ < 1. We put proof in Appendix A.2. Based on Theorem 6.6 and Definition 6.5, whenB=T, the sample complexity upper bound isO(ε −4 ), which matchesMuon’s sample complexity in Appendix Theorem A.4 and SGDM’s sample complexity in Arjevani et al. (2023); Garrigos and Gower (2023). 14 7 Conclusion In this work, motivated by HT-SR theory, we introduceHTMuonoptimizer.HTMuonmaintainsMuon’s capacity to model parameter interdependencies, while yielding more heavy-tailed updates and promoting more heavy- tailed weight spectra. It delivers consistent gains in accuracy and training stability on LLM pretraining and image classification, and can be seamlessly integrated into existingMuonvariants. To reduce theSVD overhead inHTMuon, we introduce two acceleration implementations that further lower runtime cost while still outperformingMuon. Finally, we establish a theoretical connection to steepest descent under a Schatten-q norm constraint and provide convergence analysis for smooth non-convex settings. Limitations Despite achieving improvements in LLM pretraining and image classification,HTMuonhas some potential limitations. Due to limited computational resources, we have not evaluatedHTMuonon larger models than 1B or larger-scale datasets; we leave this to future work. In addition,HTMuonandHTMuonNSincur slightly higher per-step runtime overhead thanMuon. While usingHTMuonwith larger update intervals can reduce this cost, achieving stronger performance typically requires more frequentHTMuonupdates, which increases the overall runtime. Exploring better numerical methods to accelerate HTMuon is also part of future work. Acknowledgments We thank our colleagues and funding agencies. This work is supported by the DARPA AIQ and DARPA DIAL programs, the U.S. Department of Energy under Award Number DE-SC0025584, Dartmouth College, and Lambda AI. References Yossi Arjevani, Yair Carmon, John C Duchi, Dylan J Foster, Nathan Srebro, and Blake Woodworth. Lower bounds for non-convex stochastic optimization. Mathematical Programming, 199(1):165–214, 2023. Jeremy Bernstein and Laker Newhouse. Modular duality in deep learning. arXiv preprint arXiv:2410.21265, 2024. Da Chang, Yongxiang Liu, and Ganzhao Yuan. On the convergence of muon and beyond. arXiv preprint arXiv:2509.15816, 2025. Lei Chen, Joan Bruna, and Alberto Bietti. Distributional associations vs in-context reasoning: A study of feed-forward and attention layers. arXiv preprint arXiv:2406.03068, 2024. Lizhang Chen, Jonathan Li, and Qiang Liu. Muon optimizes under spectral norm constraints. arXiv preprint arXiv:2506.15054, 2025. Xiangning Chen, Chen Liang, Da Huang, Esteban Real, Kaiyuan Wang, Hieu Pham, Xuanyi Dong, Thang Luong, Cho-Jui Hsieh, Yifeng Lu, and Quoc V Le. Symbolic discovery of optimization algorithms. In Thirty-seventh Conference on Neural Information Processing Systems, 2023. URLhttps://openreview. net/forum?id=ne6zeqLFCZ. Yatin Dandi, Luca Pesce, Hugo Cui, Florent Krzakala, Yue M Lu, and Bruno Loureiro. A random matrix theory perspective on the spectrum of learned features and asymptotic generalization capabilities. arXiv preprint arXiv:2410.18938, 2024. Damek Davis and Dmitriy Drusvyatskiy. When do spectral gradient updates help in deep learning? arXiv preprint arXiv:2512.04299, 2025. 15 Leonardo Defilippis, Yizhou Xu, Julius Girardin, Emanuele Troiani, Vittorio Erba, Lenka Zdeborov ́a, Bruno Loureiro, and Florent Krzakala. Scaling laws and spectra of shallow neural networks in the feature learning regime. arXiv preprint arXiv:2509.24882, 2025. Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei. Imagenet: A large-scale hierarchical image database. In 2009 IEEE conference on computer vision and pattern recognition, pages 248–255. Ieee, 2009. Alexey Dosovitskiy. An image is worth 16x16 words: Transformers for image recognition at scale. arXiv preprint arXiv:2010.11929, 2020. John Duchi, Elad Hazan, and Yoram Singer. Adaptive subgradient methods for online learning and stochastic optimization. Journal of machine learning research, 12(7), 2011. Pierre Foret, Ariel Kleiner, Hossein Mobahi, and Behnam Neyshabur. Sharpness-aware minimization for efficiently improving generalization. arXiv preprint arXiv:2010.01412, 2020. Guillaume Garrigos and Robert M Gower. Handbook of convergence theorems for (stochastic) gradient methods. arXiv preprint arXiv:2301.11235, 2023. Aaron Gokaslan and Vanya Cohen.Openwebtext corpus.http://Skylion007.github.io/ OpenWebTextCorpus, 2019. Dirk Groeneveld, Iz Beltagy, Evan Walsh, Akshita Bhagia, Rodney Kinney, Oyvind Tafjord, Ananya Jha, Hamish Ivison, Ian Magnusson, Yizhong Wang, Shane Arora, David Atkinson, Russell Authur, Khyathi Chandu, Arman Cohan, Jennifer Dumas, Yanai Elazar, Yuling Gu, Jack Hessel, Tushar Khot, William Merrill, Jacob Morrison, Niklas Muennighoff, Aakanksha Naik, Crystal Nam, Matthew Peters, Valentina Pyatkin, Abhilasha Ravichander, Dustin Schwenk, Saurabh Shah, William Smith, Emma Strubell, Nishant Subramani, Mitchell Wortsman, Pradeep Dasigi, Nathan Lambert, Kyle Richardson, Luke Zettlemoyer, Jesse Dodge, Kyle Lo, Luca Soldaini, Noah Smith, and Hannaneh Hajishirzi. Olmo: Accelerating the science of language models. In Proceedings of the 62nd annual meeting of the association for computational linguistics (volume 1: Long papers), pages 15789–15809, 2024. Vineet Gupta, Tomer Koren, and Yoram Singer. Shampoo: Preconditioned stochastic tensor optimization. In International Conference on Machine Learning, pages 1842–1850. PMLR, 2018. Di He, Songjun Tu, Ajay Jaiswal, Li Shen, Ganzhao Yuan, Shiwei Liu, and Lu Yin. Alphadecay: Module-wise weight decay for heavy-tailed balancing in llms. arXiv preprint arXiv:2506.14562, 2025. Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016. Nicholas J Higham and Robert S Schreiber. Fast polar decomposition of an arbitrary matrix. SIAM Journal on Scientific and Statistical Computing, 11(4):648–655, 1990. Liam Hodgkinson, Umut Simsekli, Rajiv Khanna, and Michael Mahoney. Generalization bounds using lower tail exponents in stochastic optimizers. In International Conference on Machine Learning, pages 8774–8795. PMLR, 2022. Liam Hodgkinson, Zhichao Wang, and Michael W Mahoney. Models of heavy-tailed mechanistic universality. arXiv preprint arXiv:2506.03470, 2025. Yuanzhe Hu, Kinshuk Goel, Vlad Killiakov, and Yaoqing Yang. Eigenspectrum analysis of neural networks without aspect ratio bias. In Forty-second International Conference on Machine Learning, 2025. Keller Jordan, Yuchen Jin, Vlado Boza, Jiacheng You, Franz Cesista, Laker Newhouse, and Jeremy Bernstein. Muon: An optimizer for hidden layers in neural networks, 2024. URLhttps://kellerjordan.github. io/posts/muon/. 16 Jean Kaddour, Linqing Liu, Ricardo Silva, and Matt J Kusner. When do flat minima optimizers work? Advances in Neural Information Processing Systems, 35:16577–16595, 2022. Diederik P Kingma. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014. Vignesh Kothapalli, Tianyu Pang, Shenyang Deng, Zongmin Liu, and Yaoqing Yang. From spikes to heavy tails: Unveiling the spectral evolution of neural networks. Transactions on Machine Learning Research, 2025. Alex Krizhevsky and Geoffrey Hinton. Learning multiple layers of features from tiny images. 2009. Zichong Li, Liming Liu, Chen Liang, Weizhu Chen, and Tuo Zhao. Normuon: Making muon more efficient and scalable. arXiv preprint arXiv:2510.05491, 2025. Kaizhao Liang, Lizhang Chen, Bo Liu, and Qiang Liu. Cautious optimizers: Improving training with one line of code. arXiv preprint arXiv:2411.16085, 2024. Hong Liu, Zhiyuan Li, David Hall, Percy Liang, and Tengyu Ma. Sophia: A scalable stochastic second-order optimizer for language model pre-training, 2024a. URL https://arxiv.org/abs/2305.14342. Jingyuan Liu, Jianlin Su, Xingcheng Yao, Zhejun Jiang, Guokun Lai, Yulun Du, Yidao Qin, Weixin Xu, Enzhe Lu, Junjie Yan, Yanru Chen, Huabin Zheng, Yibo Liu, Shaowei Liu, Bohong Yin, Weiran He, Han Zhu, Yuzhi Wang, Jianzhou Wang, Mengnan Dong, Zheng Zhang, Yongsheng Kang, Hao Zhang, Xinran Xu, Yutao Zhang, Yuxin Wu, Xinyu Zhou, and Zhilin Yang. Muon is scalable for llm training. arXiv preprint arXiv:2502.16982, 2025a. Liming Liu, Zhenghao Xu, Zixuan Zhang, Hao Kang, Zichong Li, Chen Liang, Weizhu Chen, and Tuo Zhao. Cosmos: A hybrid adaptive optimizer for memory-efficient training of llms. arXiv preprint arXiv:2502.17410, 2025b. Zihang Liu, Yuanzhe Hu, Tianyu Pang, Yefan Zhou, Pu Ren, and Yaoqing Yang. Model balancing helps low-data training and fine-tuning. arXiv preprint arXiv:2410.12178, 2024b. Zihang Liu, Tianyu Pang, Oleg Balabanov, Chaoqun Yang, Tianjin Huang, Lu Yin, Yaoqing Yang, and Shiwei Liu. Lift the veil for the truth: Principal weights emerge after rank reduction for reasoning-focused supervised fine-tuning. arXiv preprint arXiv:2506.00772, 2025c. Ilya Loshchilov and Frank Hutter. Decoupled weight decay regularization. arXiv preprint arXiv:1711.05101, 2017. Haiquan Lu, Yefan Zhou, Shiwei Liu, Zhangyang Wang, Michael W Mahoney, and Yaoqing Yang. Alphaprun- ing: Using heavy-tailed self regularization theory for improved layer-wise pruning of large language models. Advances in neural information processing systems, 37:9117–9152, 2024. Jianhao Ma, Yu Huang, Yuejie Chi, and Yuxin Chen. Preconditioning benefits of spectral orthogonalization in muon, 2026. URL https://arxiv.org/abs/2601.13474. Charles H Martin and Christopher Hinrichs. Setol: A semi-empirical theory of (deep) learning. arXiv preprint arXiv:2507.17912, 2025. Charles H Martin and Michael W Mahoney. Implicit self-regularization in deep neural networks: Evidence from random matrix theory and implications for learning. Journal of Machine Learning Research, 22(165): 1–73, 2021. Charles H Martin, Tongsu Peng, and Michael W Mahoney. Predicting trends in the quality of state-of-the-art neural networks without access to training or testing data. Nature Communications, 12(1):4122, 2021. 17 Takeru Miyato, Toshiki Kataoka, Masanori Koyama, and Yuichi Yoshida. Spectral normalization for generative adversarial networks. arXiv preprint arXiv:1802.05957, 2018. Team OLMo, Pete Walsh, Luca Soldaini, Dirk Groeneveld, Kyle Lo, Shane Arora, Akshita Bhagia, Yuling Gu, Shengyi Huang, Matt Jordan, et al. 2 olmo 2 furious. arXiv preprint arXiv:2501.00656, 2024. Thomas Pethick, Wanyun Xie, Kimon Antonakopoulos, Zhenyu Zhu, Antonio Silveti-Falls, and Volkan Cevher. Training deep learning models with norm-constrained lmos. arXiv preprint arXiv:2502.07529, 2025. Alec Radford, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei, Ilya Sutskever, et al. Language models are unsupervised multitask learners. OpenAI blog, 1(8):9, 2019. Colin Raffel, Noam Shazeer, Adam Roberts, Katherine Lee, Sharan Narang, Michael Matena, Yanqi Zhou, Wei Li, and Peter J Liu. Exploring the limits of transfer learning with a unified text-to-text transformer. Journal of machine learning research, 21(140):1–67, 2020. Andrei Semenov, Matteo Pagliardini, and Martin Jaggi. Benchmarking optimizers for large language model pretraining. arXiv preprint arXiv:2509.01440, 2025. Ishaan Shah, Anthony M. Polloreno, Karl Stratos, Philip Monk, Adarsh Chaluvaraju, Andrew Hojel, Andrew Ma, Anil Thomas, Ashish Tanwer, Darsh J Shah, Khoi Nguyen, Kurt Smith, Michael Callahan, Michael Pust, Mohit Parmar, Peter Rushton, Platon Mazarakis, Ritvik Kapila, Saurabh Srivastava, Somanshu Singla, Tim Romanski, Yash Vanjani, and Ashish Vaswani. Practical efficiency of muon for pretraining. arXiv preprint arXiv:2505.02222, 2025. Pratyusha Sharma, Jordan T Ash, and Dipendra Misra. The truth is in there: Improving reasoning in language models with layer-selective rank reduction. arXiv preprint arXiv:2312.13558, 2023. Wei Shen, Ruichuan Huang, Minhui Huang, Cong Shen, and Jiawei Zhang. On the convergence analysis of muon. arXiv preprint arXiv:2505.23737, 2025. Chongjie Si, Debing Zhang, and Wei Shen. Adamuon: Adaptive muon optimizer.arXiv preprint arXiv:2507.11005, 2025. James B Simon, Madeline Dickens, Dhruva Karkada, and Michael R DeWeese. The eigenlearning framework: A conservation law perspective on kernel ridge regression and wide neural networks. Transactions on Machine Learning Research, 2023. Umut Simsekli, Ozan Sener, George Deligiannidis, and Murat A Erdogdu. Hausdorff dimension, heavy tails, and generalization in neural networks. Advances in Neural Information Processing Systems, 33:5138–5151, 2020. Ilya Sutskever, James Martens, George Dahl, and Geoffrey Hinton. On the importance of initialization and momentum in deep learning. In International conference on machine learning, pages 1139–1147. pmlr, 2013. Kimi Team, Yifan Bai, Yiping Bao, Guanduo Chen, Jiahao Chen, Ningxin Chen, Ruijue Chen, Yanru Chen, Yuankun Chen, Yutian Chen, et al. Kimi k2: Open agentic intelligence. arXiv preprint arXiv:2507.20534, 2025. Hugo Touvron, Thibaut Lavril, Gautier Izacard, Xavier Martinet, Marie-Anne Lachaux, Timoth ́e Lacroix, Baptiste Rozi`ere, Naman Goyal, Eric Hambro, Faisal Azhar, et al. Llama: Open and efficient foundation language models. arXiv preprint arXiv:2302.13971, 2023. Nikhil Vyas, Depen Morwani, Rosie Zhao, Mujin Kwun, Itai Shapira, David Brandfonbrener, Lucas Janson, and Sham Kakade. Soap: Improving and stabilizing shampoo using adam. arXiv preprint arXiv:2409.11321, 2024. 18 Jinbo Wang, Mingze Wang, Zhanpeng Zhou, Junchi Yan, and Lei Wu. The sharpness disparity principle in transformers for accelerating language model pre-training. arXiv preprint arXiv:2502.19002, 2025a. Shuche Wang, Fengzhuo Zhang, Jiaxiang Li, Cunxiao Du, Chao Du, Tianyu Pang, Zhuoran Yang, Mingyi Hong, and Vincent YF Tan. Muon outperforms adam in tail-end associative memory learning. arXiv preprint arXiv:2509.26030, 2025b. Zhichao Wang, Andrew Engel, Anand D Sarwate, Ioana Dumitriu, and Tony Chiang. Spectral evolution and invariance in linear-width neural networks. Advances in neural information processing systems, 36: 20695–20728, 2023. Kaiyue Wen, David Hall, Tengyu Ma, and Percy Liang. Fantastic pretraining optimizers and where to find them. arXiv preprint arXiv:2509.02046, 2025. Yaoqing Yang, Ryan Theisen, Liam Hodgkinson, Joseph E Gonzalez, Kannan Ramchandran, Charles H Martin, and Michael W Mahoney. Test accuracy vs. generalization gap: Model selection in nlp without accessing training or testing data. In Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, pages 3011–3021, 2023. Huizhuo Yuan, Yifeng Liu, Shuang Wu, Xun Zhou, and Quanquan Gu. Mars: Unleashing the power of variance reduction for training large models. arXiv preprint arXiv:2411.10438, 2024. Aohan Zeng, Xin Lv, Qinkai Zheng, Zhenyu Hou, Bin Chen, Chengxing Xie, Cunxiang Wang, Da Yin, Hao Zeng, Jiajie Zhang, et al. Glm-4.5: Agentic, reasoning, and coding (arc) foundation models. arXiv preprint arXiv:2508.06471, 2025. Yushun Zhang, Congliang Chen, Ziniu Li, Tian Ding, Chenwei Wu, Diederik P Kingma, Yinyu Ye, Zhi-Quan Luo, and Ruoyu Sun. Adam-mini: Use fewer learning rates to gain more. arXiv preprint arXiv:2406.16793, 2024. Jiawei Zhao, Zhenyu Zhang, Beidi Chen, Zhangyang Wang, Anima Anandkumar, and Yuandong Tian. Galore: Memory-efficient llm training by gradient low-rank projection. arXiv preprint arXiv:2403.03507, 2024. Beitong Zhou, Jun Liu, Weigao Sun, Ruijuan Chen, Claire J Tomlin, and Ye Yuan. pbsgd: Powered stochastic gradient descent methods for accelerated non-convex optimization. In IJCAI, pages 3258–3266, 2020. Yefan Zhou, Tianyu Pang, Keqin Liu, Michael W Mahoney, and Yaoqing Yang. Temperature balancing, layer-wise weight analysis, and neural network training. Advances in Neural Information Processing Systems, 36:63542–63572, 2023. 19 A Proofs Notation. For a vectorv, we denote itsl 2 norm as∥v∥ 2 . For a matrixA,B ∈ R m×n , we denote its spectral norm as∥A∥, Frobenius norm as∥A∥ F , nuclear norm as∥A∥ ∗ , where∥A∥ ∗ = P r i=1 σ i , here σ i r i=1 are singular values, forp >1, Schatten−pnorm as∥A∥ p , where∥A∥ p =( P r i=1 σ p i ) 1 p . We denote ⟨A,B⟩ = Tr A ⊤ B . A.1 Update Under the Schatten−q Norm Constraint Lemma A.1 (Von Neumann trace inequality). For matrixA ∈ R m×n , letσ i (A) denote thei-th largest singular value of A, and let r = min(m,n). Then for any A,B ∈ R m×n , we have ⟨A,B⟩ = Tr A ⊤ B ≤ r X i=1 σ i (A)σ i (B). Theorem A.2. Let f G t =UΣV ⊤ be theSVDof the negative gradient f G t =−∇f(W t ), where Σ = diag(σ 1 ,...,σ r ). Letpbe the conjugate exponent ofq, satisfying 1/p ′ + 1/q= 1. The explicit solution ∆W ∗ to the optimization problem 3 is given by: ∆W ∗ = δ 1 ∥ f G t ∥ p ′ /q p ′ U Σ p ′ −1 V ⊤ ,(6) where Σ p ′ −1 = diag(σ p ′ −1 1 ,...,σ p ′ −1 r ). Proof. First, by invoking Von Neumann’s trace inequality, we have the upper bound Tr( f G t ⊤ ∆W )≤ X i σ i ( f G t )σ i (∆W ) The equality holds if and only if ∆Wshares the same left and right singular vectors as f G t . Thus, the optimal ∆W must take the form: ∆W = U ΛV ⊤ ,where Λ = diag(λ 1 ,...,λ r ).(7) Here,λ i ≥0 represents the singular values of ∆W. The problem reduces to maximizing P i σ i λ i subject to ∥λ∥ q ≤ δ. By H ̈older’s inequality, we have the upper bound: X i σ i λ i ≤∥σ∥ p ′ ∥λ∥ q ≤∥σ∥ p ′ δ.(8) Equality is achieved if and only if λ q i ∝ σ p ′ i , which implies: λ i = cσ p ′ /q i ,(9) for some constant c≥ 0. Using the conjugate property 1/p ′ + 1/q = 1, the exponent simplifies to: p ′ q = p ′ 1− 1 p ′ = p ′ − 1.(10) Substituting this back yields the optimal singular values: λ i = c· σ p ′ −1 i ,(11) 20 wherecis a normalization constant. To determinec, we substitute these values into the active constraint ∥∆W∥ q = δ: X i (cσ p ′ −1 i ) q ! 1/q = c X i σ (p ′ −1)q i ! 1/q = δ.(12) Using the algebraic relation (p ′ −1)q=p ′ , the term involving singular values simplifies to the Schatten−p ′ norm of f G t . Solving for c, we obtain: c∥ f G t ∥ p ′ /q p ′ = δ =⇒ c = δ ∥ f G t ∥ p ′ /q p ′ .(13) Substituting c back into the expression ∆W = U ΛV ⊤ yields the stated closed-form solution.□ A.2 Convergence Analysis We consider the following stochastic optimization problem: min W∈R m×n f (W ) = E ξ [f (W ;ξ)], where ξ is stochastic noise. We denote f ∗ = inf W f (W ), r = minm,n. We assume f ∗ >−∞. We give the convergence analysis under a smooth non-convex setting for findingε-stationary points, i.e. ∥∇f (x)∥ ∗ ≤ ε. To better understand the proof, we give a stochastic version ofHTMuonin Algorithm 7 for convergence analysis. Algorithm 7 HTMuonStochastic 1:Input: Initial weights W 0 ∈ R m×n , loss functionL, learning rateη t , momentum parameterβ, batch size B, power p∈ (0, 1). 2: Initialize M 0 ∈ R m×n ← 0 3: for t = 1, 2,... do 4:Sample batch ξ t,i B i=1 uniformly 5:G t = 1 B P B i=1 ∇ W L (W t ;ξ t,i ) 6:M t ← βM t−1 + (1− β)G t 7: U t , Σ t ,V ⊤ t ← SVD(M t ) 8: O t ← U t Σ p t V ⊤ t 9:W t+1 ← W t − η t O t 10: end for Lemma A.3. Fort= 0,1,· ,T, M t andW t are generated by Algorithm 7. Consider thatM 0 =G 0 ,M t = βM t−1 + (1− β)G t ,under Assumptions 6.2 and 6.3, if we assumeη=maxη t T t=1 , under Assumption 6.4, we have E∥∇f (W t )−M t ∥ F ≤ s 1− β 1 + β σ √ B + β t σ √ B + √ rl p βLη 1− β , where G t = 1 B P B i=1 ∇f (W t ;ξ t,i ) and B is batch size. 21 Proof. We defineC 0 =∇f(W 0 ),C t =βC t−1 + (1−β)∇f(W t ) = (1−β) P t i=1 β t−i ∇f (W i ) +β t ∇f(W 0 ), under Assumptions 6.2 and 6.3 we note that E[∥∇f (W t −C t )∥ F ] =E[∥∇f (W t )− (βC t−1 + (1− β)∇f (W t ))∥ F ] =E[β∥∇f (W t −C t−1 )∥ F ] ≤E β∥∇f (W t−1 −C t−1 )∥ F +E β∥∇f (W t−1 ))−∇f (W t ))∥ F ≤E β∥∇f (W t−1 −C t−1 )∥ F +E βL∥W t−1 −W t ∥ F ≤E β∥∇f (W t−1 −C t−1 )∥ F +E βLη∥U t−1 Σ p t−1 V ⊤ t−1 ∥ F ≤E β∥∇f (W t−1 −C t−1 )∥ F + βLη √ rl p ≤β t ∥∇f (W 0 )−C 0 ∥ F + t X i=1 β i Lη √ rl p ≤ √ rl p βLη 1− β . And follow the same proof in Section B.1 in Shen et al. (2025), we have E[∥C t −M t ∥ F ]≤ s 1− β 1 + β σ √ B + β t σ √ B . So we have E∥∇f (W t )−M t ∥ F ≤E[∥∇f (W t −C t )∥ F ] + E[∥C t −M t ∥ F ] ≤ s 1− β 1 + β σ √ B + β t σ √ B + √ rl p βLη 1− β □ Theorem A.4 (Muon (Shen et al., 2025)). Under Assumptions 6.2 and 6.3, if we apply Muon withη t =η, then 1 T T−1 X t=0 E ∥∇f (W t )∥ ∗ ≤ E f (W 0 )− f (W T ) Tη + Lrη 2 + 2σ p r(1− β) p B(1 + β) + 2βσ √ r (1− β)T √ B + 2rηβL 1− β . Where B is batch size. Denote ∆ =f(W 0 )−f ∗ . If we setB= 1,η= q (1−β)∆ rTL , 1−β=min n √ L∆ σ √ T , 1 o , then 1 T T−1 X t=0 E ∥∇f (W t )∥ ∗ ≤ O 4 r r 2 L∆σ 2 T + r rL∆ T + √ rσ 2 √ L∆T ! . Thus, Muon can find an ε-nuclear norm stationary point of f with a complexity of O(r 2 Lσ 2 ∆ε −4 ). 22 Theorem A.5 (HTMuon: Theorem 6.6). Under Assumptions 6.2, 6.3 and 6.4, if we applyHTMuonin Algorithm 7 with adaptive learning rate η t = ⟨M t ,ρ(M t )⟩ L∥ρ(M t )∥ 2 F , ρ(M t ) = U t Σ p t V ⊤ t and batch size B = T , then 1 T T X t=1 E ∥∇f (W t )∥ 2 1+p ≤ 2Lr 1 s 1− ε ′ 1− ε ′ L + 1 Lε ′ 2(1− β)σ B(1 + β) + 2β 2 σ 2 (1− β 2 )BT + 2rl 2p β 2 L 2 η 2 (1− β) 2 + 4Lr 1 s 1− ε ′ f (W 0 )− f (W T ) T . Where s = 1+p 1−p and 0 < ε ′ < 1. Furthermore, if we define β = 1 √ T , ∆ = f (W 0 )− f ∗ , thus we have 1 T T X t=1 E ∥∇f (W t )∥ 2 1+p ≤ 2Lr 1 s 1− ε ′ 1− ε ′ L + 1 Lε ′ 2σ T + 4σ 2 T + 4rl 2p L 2 η 2 T + 4Lr 1 s 1− ε ′ ∆ T . Proof. Consider M t = U t Σ t V ⊤ t , we first define ρ(M t ) = U t Σ p t V ⊤ t . By Assumption 6.2, we have E f (W t )− f (W t+1 ) ≥ E η t ⟨∇f (W t ),ρ(M t )⟩− Lη 2 t 2 ∥ρ(M t )∥ 2 F ≥ E η t ⟨M t ,ρ(M t )⟩− Lη 2 t 2 ∥ρ(M t )∥ 2 F − η t ⟨M t −∇f (W t ),ρ(M t )⟩ . Here we consider adaptive learning rate η t = ⟨M t ,ρ(M t )⟩ L∥ρ(M t )∥ 2 F > 0, thus E f (W t )− f (W t+1 ) ≥ E ⟨M t ,ρ(M t )⟩ 2 L∥ρ(M t )∥ 2 F − ⟨M t ,ρ(M t )⟩ 2 2L∥ρ(M t )∥ 2 F − η t ∥M t −∇f (W t )∥ F ∥ρ(M t )∥ F ≥ E ⟨M t ,ρ(M t )⟩ 2 L∥ρ(M t )∥ 2 F − ⟨M t ,ρ(M t )⟩ 2 2L∥ρ(M t )∥ 2 F − ⟨M t ,ρ(M t )⟩ L∥ρ(M t )∥ F ∥M t −∇f (W t )∥ F By 2ab≤ ε ′ a 2 + b 2 ε ′ , where 0 < ε ′ < 1, we have E f (W t )− f (W t+1 ) ≥ E ⟨M t ,ρ(M t )⟩ 2 L∥ρ(M t )∥ 2 F − ⟨M t ,ρ(M t )⟩ 2 2L∥ρ(M t )∥ 2 F − ε ′ 2L ⟨M t ,ρ(M t )⟩ 2 ∥ρ(M t )∥ 2 F − 1 2Lε ′ ∥M t −∇f (W t )∥ 2 F ≥ E 1− ε ′ 2L ⟨M t ,ρ(M t )⟩ 2 ∥ρ(M t )∥ 2 F − ∥M t −∇f (W t )∥ 2 F 2Lε ′ By Holder’s Inequality, for p∈ (0, 1) and we set s = 1+p 1−p and m = 1+p 2p , we have ∥ρ(M t )∥ 2 F = r X i=1 σ 2p i ≤ r 1 s r X i=1 (σ 2p i ) m ! 1 m ≤ r 1 s r X i=1 σ 1+p i ! 2p p+1 . Thus we have ⟨M t ,ρ(M t )⟩ 2 ∥ρ(M t )∥ 2 F ≥ P r i=1 σ 1+p i 2 r 1 s P r i=1 σ 1+p i 2p p+1 = ∥M t ∥ 2 1+p r 1 s . Hence, we have E f (W t )− f (W t+1 ) ≥ E 1− ε ′ 2Lr 1 s ∥M t ∥ 2 1+p − ∥M t −∇f (W t )∥ 2 F 2Lε ′ 23 By Lemma A.3, we let η = maxη t T t=1 , we have E∥∇f (W t )−M t ∥ 2 F ≤ s 1− β 1 + β σ √ B + β t σ √ B + √ rl p βLη 1− β ! 2 ≤ 2(1− β)σ B(1 + β) + 2β 2t σ 2 B + 2rl 2p β 2 L 2 η 2 (1− β) 2 Thus we have E 1− ε ′ 2Lr 1 s ∥∇f (W t )∥ 2 1+p ≤ E 1− ε ′ Lr 1 s ∥∇f (W t )−M t ∥ 2 1+p + E 1− ε ′ Lr 1 s ∥M t ∥ 2 1+p ≤ E 1− ε ′ L ∥∇f (W t )−M t ∥ 2 F + E 1− ε ′ Lr 1 s ∥M t ∥ 2 1+p ≤ E 1− ε ′ L + 1 Lε ′ ∥∇f (W t )−M t ∥ 2 F + f (W t )− f (W t+1 ) Thus we have 1 T T X t=1 E ∥∇f (W t )∥ 2 1+p ≤ 1 T T X t=1 E 2Lr 1 s 1− ε ′ 1− ε ′ L + 1 Lε ′ ∥∇f (W t )−M t ∥ 2 F + 2Lr 1 s 1− ε ′ f (W 0 )− f (W T ) T ≤ 2Lr 1 s 1− ε ′ 1− ε ′ L + 1 Lε ′ 2(1− β)σ B(1 + β) + 2β 2 σ 2 (1− β 2 )BT + 2rl 2p β 2 L 2 η 2 (1− β) 2 + 2Lr 1 s 1− ε ′ f (W 0 )− f (W T ) T We define B = T,β = 1 √ T , ∆ = f (W 0 )− f ∗ , thus we have 1 T T X t=1 E ∥∇f (W t )∥ 2 1+p ≤ 2Lr 1 s 1− ε ′ 1− ε ′ L + 1 Lε ′ 2σ T + 2σ 2 (1− β 2 )T 3 + 2rl 2p L 2 η 2 T (1− β) 2 + 2Lr 1 s 1− ε ′ ∆ T . When β = 1 √ T , it is easy to get 1− β ≥ 1 2 , thus we have 1 T T X t=1 E ∥∇f (W t )∥ 2 1+p ≤ 2Lr 1 s 1− ε ′ 1− ε ′ L + 1 Lε ′ 2σ T + 4σ 2 T 3 + 4rl 2p L 2 η 2 T + 2Lr 1 s 1− ε ′ ∆ T . Based on Theorem 6.6 and Definition 6.5, whenB=T, the sample complexity upper bound isO(ε −4 ), which matchesMuon’s sample complexity in Theorem A.4 and (Chang et al., 2025) andSGDM’s sample complexity in Arjevani et al. (2023); Garrigos and Gower (2023). .□ A.3 Lemma for PL exponent α Lemma A.6. Suppose the singular values of matrix W ∈ R n×m follow s k = s 1 k −s , 1≤ k ≤ n, we have PL exponent α of W satisfies α = 1 + 1 2s . Proof. Since we haves k =s 1 k −s ,1≤ k ≤ n, which means eigenvaluesλ k =λ 1 k −2s ,1≤ k ≤ n. Here we suppose Λ is a random variable distributed according to the empirical distribution from these eigenvalues, we can see that the distribution function takes the following form: P(Λ > λ 1 k −2s ) = k n . 24 By changing variables λ 1 k −s = λ, we get the cumulative distribution function of Λ: P(Λ > λ)∼ λ − 1 2s . After that, we take the derivative with respect to λ, and we get the ESD: p(λ)∼ λ − ( 1 2s +1 ) . So we have α = 1 + 1 2s .□ B More Experimental Details In this section, we present more experimental details. B.1 Muon and HTMuon Algorithms In Algorithm 1 and 2, we give the implementations ofMuonNSandMuonSVD.MuonNScorresponds to the commonly used five-stepNewton Schulzimplementation ofMuon, whereasMuonSVDrepresents the theoretically exact SVD-based implementation. Empirically, we find that due to the numerical approximation inMuonNS, the singular values of the update matrix do not always equal 1 exactly. Such deviations may suppress weights along noise-dominated eigenvector directions, which could explain whyMuonNSperforms better than MuonSVD in practice. In Algorithm 4, we give the implementations ofHTMuonNS. As discussed in Section 5.4, we note that the HTMuon update O t = U t Σ p t V ⊤ t admits the factorization O t = (U t V ⊤ t ) (V t Σ p t V ⊤ t ). ForU t V ⊤ t part, this isMuon’s update, we can considerNewton-Schulzwith 5 steps. To compute the symmetric factor V t Σ p t V ⊤ t efficiently, we apply the NS root routine to M ⊤ t M t in Algorithm 5, since M ⊤ t M t = V t Σ 2 t V ⊤ t , we haveV t Σ p t V ⊤ t = (M ⊤ t M t ) p 2 . Accordingly,NSrootusesNewtonSchulziterations to approximate matrix square roots and applies successive square-root operations for [1− log 2 p] rounds. B.2 Training Loss Curve In Figure 9, we show the smoothed training loss curves for LLaMA-60M and LLaMA-135M on C4 dataset. We find in the later stages of training,HTMuonconsistently achieves a noticeably lower training loss than Muon, demonstrating the effectiveness and training stability. B.3 More Experiments Results In Table 7, we conduct evaluations on LLaMA-60M and LLaMA-135M on the C4 datatset, comparing HTMuonwith multipleMuon-variant optimizers as well as other state-of-the-art optimizers.HTMuonconsistently achieves the best results. Furthermore, combiningHTMuonwith otherMuon-variant optimizers will achieve better performance. Detailed hyperparameters are reported in Appendix C. In Table 8 and 9, we report the perplexity and average per-step runtime forHTMuon,HTMuonNS, and HTMuon+NorMuonunder different update intervals. We measure per-step runtime on 4 NVIDIA A6000 GPUs. We find that increasing the interval substantially reduces the average overhead ofHTMuon, while it still outperformsMuon. Under the same interval,HTMuonNSfurther reduces the average overhead and also remains better thanMuon. Similarly, forHTMuon+NorMuon, larger intervals make its average overhead approach that of NorMuon, while consistently outperformingNorMuon. For more visualizations, please refer to Figure 5 and 10. 25 40006000800010000 Steps 3.25 3.50 3.75 4.00 4.25 Training Loss Muon HTMuon (a) LLaMa-60M on C4 5000100001500020000 Steps 3.00 3.25 3.50 3.75 4.00 4.25 Training Loss Muon HTMuon (b) LLaMA-135M on C4 Figure 9: Training loss curves for LLaMA-60M and LLaMA-135M. Learning rate for both models and optimizers is 0.03. Both curves are smoothed via a simple moving average (uniform weights, window = 50). Table 7: We evaluateHTMuonon LLaMA-60M and LLaMA-135M, comparing it against severalMuonvariant optimizers and other state-of-the-art optimizers. Red indicates the best value, and blue denotes the second- best value. We find thatHTMuonconsistently outperforms these baselines. More hyperparameter details are provided in Appendix C. Muon HTMuon NorMuon HTMuon+NorMuon AdaMuon MARS SOAP Cautious COSMOS GaLore Sophia LLaMa-60M28.8027.8828.1727.5528.6729.1328.9629.7128.6234.2634.02 LLaMa-135M22.2321.2521.9921.1122.4322.4222.4323.1522.1525.1125.63 NorMuonInterval=1Interval=5Interval=10Interval=25Time(sec) 27.0 27.5 28.0 28.5 Perplexity ( ↓ ) 0.0 0.5 1.0 1.5 2.0 (a) HTMuon+NorMuon LLaMA-60M 21.0 21.5 22.0 22.5 0.0 0.5 1.0 1.5 2.0 Time (sec) (b) HTMuon+NorMuon LLaMA-135M Figure 10: we evaluate applyingHTMuon+NorMuonon LLaMA-60M and LLaMA-135M every 1, 5, 10, and 25 steps (while other steps applyingNorMuon) . We report the average per-step runtime overhead for all methods. Detailed results and hyperparameter settings are provided in Table 8 and 9 in Appendix C. In Table 10, we summarize the PPL and wall-clock time ofMuon,HTMuon,HTMuonNS,HTMuon(Interval = 5), andHTMuonNS(Interval = 5). To better report the wall-clock time, we rerun training on NVIDIA RTX PRO 6000 GPUs: LLaMA-60M on C4 using 2 GPUs, and LLaMA-135M on C4 using 4 GPUs. We find that althoughHTMuoncan take longer thanMuon, the substantial gains achieved through heavy-tailed spectral correction suggest thatHTMuoncaptures the correct inductive bias for learning. Furthermore, we believe that efficient variants such asHTMuonNSwith interval-based updates provide a strong balance between performance and efficiency. Here we also give the FLOPs analysis for Muon and HTMuon NS, we consider initial weights W 0 ∈ R m×n . 26 Table 8: We report the exact PPL obtained when training LLaMA-60M and LLaMA-135M usingHTMuon, HTMuonNS, andHTMuon+NorMuonunder different update intervals. ForMuonandNorMuon, interval updates are not applicable; therefore, we only report their PPL after training. LLaMA-60MLLaMA-135M Optimizer151025151025 Muon28.8022.23 HTMuon27.8828.3628.5828.6421.2522.0122.1322.20 HTMuonNS28.0328.4228.6728.7421.3522.0522.1422.20 NorMuon 28.1721.99 HTMuon+NorMuon27.5527.8528.0128.1221.1121.6121.8121.95 Table 9: We report the average per-step runtime overhead obtained when training LLaMA-60M and LLaMA- 135M usingHTMuon,HTMuonNS, andHTMuon+NorMuonunder different update intervals. We measure per-step runtime on 4 NVIDIA A6000 GPUs. ForMuonandNorMuon, interval updates are not applicable; therefore, we only report average per-step runtime overhead. LLaMA-60MLLaMA-135M Optimizer151025151025 Muon0.510.59 HTMuon1.260.670.610.541.700.820.700.63 HTMuonNS0.660.590.560.540.920.680.640.63 NorMuon0.540.60 HTMuon+NorMuon1.270.700.640.591.780.860.810.72 • FLOPs analysis for Muon(Algorithm 1): F Muon = 20mnr + O(mn). • FLOPs analysis for HTMuon NS(Algorithm 4): F HTMuonNS = 20mnr + 4mn 2 + 6LTn 3 + O(n 2 ), where r = min(m,n),L = h log 2 2 p i and T=ns steps in Algorithm 4. Compared toMuon, the additional FLOPs inHTMuonNSmainly come from theNewtonSchulziterations used for the fractional power approximation. In Figure 11a, 11b, we visualize the layer-wise PL exponentαfor ResNet18 on CIFAR-100 and CIFAR-10 dataset to further support that HTMuon makes the updated matrices more heavy-tailed than Muon. B.4 More Ablation Study B.4.1 Varying learning rates In Figure 12, Table 11, 12 and 13, we conduct learning grid search on LLaMA and ResNet models, our experiments show that for most learning rates,HTMuonoutperformsMuonandNorMuon, which exhibits our method’s effectiveness and robustness. 27 Table 10: We report the perplexity and wall-clock time of when training LLaMA-60M and LLaMA-135M usingMuon,HTMuon,HTMuonNS,HTMuon(Interval=5),HTMuonNS(Interval=5). We measure wall-clock time on 4 NVIDIA RTX PRO 6000 GPUs. LLaMA-60MLLaMA-135M OptimizerPPLTimePPLTime(hour) Muon28.840.6722.271.40 HTMuonNS(Interval=5)28.500.7122.021.55 HTMuon(Interval=5)28.361.0222.042.40 HTMuon NS28.060.8721.362.17 HTMuon27.872.4621.246.57 15101520 Layer index 2 3 4 5 6 α Muon ( ̄ α=2.879) HTMuon ( ̄ α=2.492) (a) ResNet18 on CIFAR-100 15101520 Layer index 1 2 3 4 5 α Muon ( ̄ α=2.888) HTMuon ( ̄ α=2.511) (b) ResNet18 on CIFAR-10 11020304050 Layer index 10 20 α Muon ( ̄ α=4.786) HTMuon_HT ( ̄ α=4.376) (c) HTMUONHT60M 11020304050607080 Layer index 5 10 α Muon ( ̄ α=4.750) HTMuon_HT ( ̄ α=3.809) (d) HTMUONHT135M Figure 11: (a)(b): Layer-wise PLαfor ResNet-18 model weight ESDs on CIFAR-100 and CIFAR-10 trained withMuonandHTMuon. (c)(d): Layer-wise PLαfor LLaMA model weight ESDs trained with MuonandHTMuonHT. All models used for visualization are trained using each optimizer’s best-performing hyperparameter configuration. For hyperparameter configurations, please refer to Appendix C. B.4.2 Varying different p In table 14, we give the exact results for ablation study onpin Figure 8, we test multiple values ofpon LLaMA-60M and LLaMA-135M on the C4 dataset, and find thatp= 0.125 is a robust choice. We therefore use p = 0.125 in most experiments. Here we discuss more about hyperparameter sensitivity forp. Although introducingpadds a tuning dimension and add hyperparameter search burden, also we only do grid search on LLaMA models to obtain the optimal p = 0.125, we will clarify that p = 0.125 can serve as a recommended default for other tasks. • p= 0.125 generalizes well across tasks and architectures. In our paper, we also usep= 0.125 for HTMuonwhen training GPT-2 small on OpenWebText dataset, ResNet-18 and ResNet-50 on CIFAR datasets, Tables 2, 12and 13 in our paper show that HTMuon consistently outperforms Muon. • p = 0.125 is often optimal or near-optimal across tasks and architectures. To explore whether p= 0.125 is also close to optimal on other tasks, we performed grid searches overpon the CIFAR and OpenWebText datasets in Table 15 and 16. We find that on CIFAR-10,p= 0.125 achieves the best performance for both ResNet-18 and ResNet-50, and on OpenWebText, p=0.125 also yields the best result for GPT-2 small. While on CIFAR-100 a larger value (e.g.,p= 0.2) can perform slightly better, this does not change our overall conclusion that p=0.125 is a strong and practical recommended default. B.4.3 Varying different ways to make Muon update more heavy-tailed In Algorithm 6, we explore another wayHTMuonHTto makeMuonupdates more heavy-tailed , which directly replaces the momentum matrix’s singular values with a pre-specified heavy-tailed distribution. We use α = 0.25 in Algorithm 6 for experiments. In Table 17, we find that, although HTMuonHT is not as strong as 28 MuonHTMuonNorMuon HTMuon+NorMuon 0.010.020.030.04 Learning rate 27 28 29 30 Perplexity ( ↓ ) (a) LLaMA-60M on C4 0.010.020.030.04 Learning rate 21 22 23 24 (b) LLaMA-135M on C4 0.0150.020.0250.03 Learning rate 77.0 77.5 78.0 78.5 Accuracy ( ↑ ) (c) ResNet18 on CIFAR-100 0.0150.020.0250.03 Learning rate 79.0 79.5 80.0 80.5 Accuracy ( ↑ ) (d) ResNet50 on CIFAR-100 0.010.0150.020.025 Learning rate 95.0 95.2 95.4 95.6 95.8 (e) ResNet18 on CIFAR-10 0.010.0150.020.025 Learning rate 95.6 95.8 96.0 96.2 (f) ResNet50 on CIFAR-10 Figure 12: We conduct learning-rate grid searches forMuon,NorMuon,HTMuon, andHTMuon+NorMuonacross various datasets and architectures, including LLaMA models on C4 and ResNet models on CIFAR-100 and CIFAR-10. We find thatHTMuonoutperformsMuonandNorMuonacross most learning rates, demonstrating the effectiveness and robustness of our method. HTMuon , it consistently outperformsMuon. In Figure 11c and 11d, we visualize the layer-wise PL exponentα for LLaMA-60M and LLaMA-135M, and observe thatHTMuonHTalso yields consistently smallerαthanMuon across layers. This suggests that the proposed heavy-tailed spectral correction can indeed improve Muon. Since some readers maybe confused that the performance improvement attribution of our methodHTMuon is from inducing heavy tails or preserving singular value information. Here we discuss more about the performance improvement attribution for our algorithmHTMuonbased onHTMuonHT. Our motivation for heavy-tail spectral correction is thatMuonenforces unit singular values for update matrices. This induces a light-tailed spectrum and can over-emphasize noise-dominated directions, which may limit performance. Guided by HT-SR theory, we therefore consider heavy-tailed spectral correction. We study two ways of introducing heavy-tailed spectra:HTMuonHT(which merely induces heavy tails) andHTMuon(which introduces heavy tail spectral correction on top of preserved spectral information). Below, we discuss the benefits of merely inducing heavy tails versus preserving singular value information. •Merely inducing a heavy-tailed spectrum already yields non-trivial gains. OurHTMuonHT ablation enforces a rigid, pre-specified heavy-tailed spectral distribution. In Table 12 of our paper, although this design discards the data-driven singular-value structure,HTMuonHTstill outperforms Muonin perplexity by 0.28 on LLaMA-60M and by 0.58 on LLaMA-135M. We would like to kindly emphasize that this is not ”poor performance.” First,Muonis a strong baseline: in Table 1,Muon outperforms the widely used Adam/AdamW by up to 3.21 PPL on LLaMA-60M and 135M. Under such a strong baseline, simply replacing the spectrum with a pre-defined heavy-tailed distribution that ignores data information yet still yields gains provides direct evidence that heavy-tail correction itself 29 Table 11: Detailed results for learning rate grid search on C4. Red indicates the best value, and blue denotes the second-best value. LLaMA-60MLLaMA-135M Optimizer0.010.020.030.040.010.020.020.04 Muon28.8028.9928.8829.9322.2322.5122.5923.73 NorMuon 28.6328.5928.1729.0721.9922.1722.0622.76 HTMuon28.4828.7527.8830.2222.9521.9621.2521.92 HTMuon+NorMuon30.1428.6327.5528.0322.7121.8221.1121.76 Table 12: Detailed results for learning rate grid search on CIFAR-100. Red indicates the best value, and blue denotes the second-best value. ResNet18ResNet50 Optimizer0.0150.020.0250.030.0150.020.0250.03 Muon77.1077.1677.9577.5379.1579.4179.6779.85 NorMuon76.7577.8277.6477.6379.1379.7879.7479.61 HTMuon77.1877.6578.2478.2179.5580.0380.1680.07 HTMuon+NorMuon77.5177.7878.5178.0279.6080.0379.9380.22 is effective. Second, the improvement ofHTMuonHTis not insignificant: as suggested by Table 2 in (Liu et al., 2025b) and (He et al., 2025) and Figure 2 in (Wen et al., 2025), after carefully tune the hyperparameters of the baselines on LLaMA/C4, an improvement of≥0.2 PPL overMuonis generally regarded as non-negligible. For example,COSMOSoutperformsMuonby 0.15 PPL for LLaMA-135M in (Liu et al., 2025b) andAlphaDecayoutperforms Adam by 0.11 PPL for LLaMA-1B in (He et al., 2025). Therefore, these results strongly indicate that merely inducing heavy-tail correction is already useful. Moreover, in Figure 1 of our motivating experiment, we observe thatMuonNSoutperformsMuonSVD, this further suggests that strictly enforcing a light-tailed, unit-singular-value spectrum may be overly restrictive, and that relaxing it via heavy-tail correction can be beneficial. •Heavy tail spectral correction on top of preserved spectral information is more effective. We acknowledge that the data-driven spectral information encoded in singular values is important, as it reflects meaningful structure learned from data. Our methodHTMuontherefore performs heavy-tail spectral correction based on the current spectral information: specifically, we raise singular values to the powerp(as motivated by Lemma 4.1) to give the spectrum a heavier tail thanMuon’s unit-spectrum, without discarding the learned geometry. This combination preserves directional information while attenuating noise-dominated directions more than signal-aligned ones compared toMuon, leading to stronger gains than using a rigid, pre-specified heavy-tailed spectrum alone. C Detailed Hyperparameters Since in practice ourMuonis a hybridMuonwithAdamWvariant(Jordan et al., 2024), we use the same hyperpa- rameter settings as AdamW for the small subset of parameters updated with AdamW. C.1 Hyperparameter settings for LLM pretraining For the C4 task, we adopt the LLaMA architecture specifications listed in Table 18 to ensure reproducibility and consistency with prior work. All model variants are trained with a uniform maximum sequence length of 30 Table 13: Detailed results for learning rate grid search on CIFAR-10. Red indicates the best value, and blue denotes the second-best value. ResNet18ResNet50 Optimizer0.010.0150.020.0250.010.0150.020.025 Muon95.1095.3995.2595.3695.895.9095.9795.90 NorMuon 95.3095.3595.5795.2495.7496.0395.9195.75 HTMuon95.2795.2395.6395.3996.1395.8496.0795.86 HTMuon+NorMuon95.3695.2595.4095.1796.0295.8596.0495.97 Table 14: Detailed results for ablation study onpon C4 dataset. Red indicates the best value, and blue denotes the second-best value. LLaMA-60MLLaMA-135M Optimizer00.10.1250.20.2500.10.1250.20.25 HTMuon28.8027.91 27.88 28.56 29.9722.23 21.33 21.25 21.66 22.43 HTMuon+NorMuon28.17 27.4827.5528.15 29.5421.99 21.21 21.11 21.3322.1 256, a batch size of 512, and an aggregate of 13K tokens per batch. For LLaMA-60M, we run all experiments on two NVIDIA L40 GPUs or two NVIDIA A6000 GPUs without gradient accumulation; for LLaMA-135M, we run all experiments on four NVIDIA L40 GPUs or four NVIDIA A6000 GPUs without gradient accumulation; for LLaMA-350M, we run all experiments on four NVIDIA RTX PRO 6000 Blackwell without gradient accumulation; for LLaMA-1B, we run all experiments on four NVIDIA RTX PRO 6000 Blackwell with gradient accumulation. For the OpenWebText task, we consider a uniform maximum sequence length of 1024 , a batch size of 480 and 10000 iterations for GPT-2 small. we run all experiments on four NVIDIA RTX PRO 6000 Blackwell with gradient accumulation=8. We do learning rate grid search on0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1forMuonandHTMuon, we set learning rate= 6e−4 forAdamW. We consider weight decay= 0.1 for all optimizers. We applyp= 0.125 for all LLM pretraining experiments. In our experiments, we apply Muon updates to the embedding and output layers during training. For more detailed hyperparameter settings, please refer to Table 19 and 20. C.2 Hyperparameter settings for Image classification For training on CIFAR datasets, we set batch size = 512 and we setp= 0.125. We run all the experiments on one NVIDIA L40 GPU. For more detailed hyperparameter settings, please refer to Table 21 and 22. For training on ImageNet-1K datasets, we set batch size = 1024 and we setp= 0.03125. We run all the experiments on one NVIDIA RTX PRO 6000 Blackwell. We set learning rate for0.003, 0.004, 0.005for Adam, Muon, HTMuon. D Baseline Optimizers In this section, we provide the algorithms for all optimizers evaluated in our study. We adopt the following notation:w t denotes the model parameters at stept,g t the corresponding gradient,ηthe learning rate,λ the weight decay,β 1 andβ 2 the moment decay rates,εa numerical stability constant,g norm the gradient norm, andmandvthe first and second moments, respectively. All operations are element-wise unless stated otherwise. 31 Table 15: Detailed results for ablation study onpforHTMuonon CIFAR datasets. Red indicates the best value. ResNet18ResNet50 Dataset00.10.1250.20.2500.10.1250.20.25 CIFAR-1095.39 95.14 95.6395.4695.4995.97 95.88 95.9595.9595.95 CIFAR-100 77.95 77.6578.24 78.58 78.3279.67 79.4080.16 80.30 80.27 Table 16: Detailed results for ablation study onpforHTMuonon OpenWebText dataset. Red indicates the best value. Optimizer00.10.1250.20.25 HTMuon22.46 22.40 22.20 22.49 22.83 Table 17: Detailed results forHTMuonHT. Red indicates the best value, and blue denotes the second-best value. We use α = 0.25 and learning rate = 0.03 in Algorithm 6 for experiments. Muon HTMuon NorMuon HTMuon+NorMuon HTMuonHT HTMuonHT+NorMuon LLaMa-60M28.827.8828.1727.5528.5228.02 LLaMa-135M22.2321.2521.9921.1121.6521.26 Table 18: Hyperparameters of LLaMA models. ParamsHiddenIntermediateHeadsBlocksStepsData amountBatch Size 60M51213768810K1B512 135M7682048121220K2B512 350M10242736162460K6B512 1B20485461322490K9B512 NorMuon(Li et al., 2025): AMuon-based optimizer that applies normalized and orthogonalized gradient updates with shape-aware scaling. We put the implementations in Algorithm 8. AdaMuon(Si et al., 2025): An adaptive variant ofMuonthat incorporates second-moment information into orthogonalized updates. We put the implementations in Algorithm 9. MARS (Yuan et al., 2024): A momentum-based adaptive optimizer included as a representative adaptive baseline. We put the implementations in Algorithm 10. SOAP(Vyas et al., 2024): An optimizer that applies stochastic orthogonalization or projection to gradient updates. We put the implementations in Algorithm 11. Cautious(Liang et al., 2024): An optimizer that modifies update application in a conservative manner based on gradient information. We put the implementations in Algorithm 12. COSMOS(Liu et al., 2025b): An optimizer that incorporates geometry-aware scaling into its parameter updates. We put the implementations in Algorithm 13. GaLore (Zhao et al., 2024): A memory-efficient optimizer that performs low-rank gradient projection to reduce optimizer-state and update costs, enabling large-model training under limited GPU memory. We put the implementations in Algorithm 14. Sophia (Liu et al., 2024a): : A second-order optimizer that uses a Hessian-based (diagonal) curvature estimate to precondition gradients and applies clipped updates for stability and efficiency in large-scale training. We put the implementations in Algorithm 15. 32 Table 19: Hyperparameters for LLaMA-60M and LLaMA-135M On C4 Dataset. The boldfaced values denote the optimal hyperparameters. LLaMA-60MLLaMA-135M OptimizerLRWDLRWD Adam1e-31e-51e-31e-5 AdamW1e-30.11e-30.1 Muon0.01,0.02,0.03,0.040.10.01,0.02,0.03,0.040.1 NorMuon0.01,0.02, 0.03,0.040.10.01,0.02,0.03,0.040.1 AdaMuon6e-4, 1e-3, 5e-3, 1e-20.1 1e-3, 5e-3, 1e-2, 2e-20.1 MARS1e-3, 2e-3, 3e-3, 4e-30.11e-3, 2e-3, 3e-3, 4e-30.1 SOAP2e-3, 3e-3, 4e-30.01 2e-3, 3e-3, 4e-30.01 Cautious1e-3, 2e-3, 3e-3, 4e-301e-3, 2e-3, 3e-3, 4e-30 COSMOS1e-3, 2e-3, 3e-3, 4e-30.11e-3, 2e-3, 3e-3, 4e-30.1 GaLore5e-3, 1e-2, 2e-2, 3e-2, 4e-205e-3, 1e-2, 2e-2, 3e-2, 4e-20 Sophia1e-4, 2e-4,3e-4, 4e-4, 5e-4, 6e-40.11e-4, 2e-4,3e-4, 4e-4, 5e-4, 6e-40.1 HTMuon0.01,0.02,0.03,0.040.10.01,0.02,0.03,0.040.1 HTMuon+NorMuon0.01,0.02,0.03,0.040.10.01,0.02,0.03,0.040.1 Table 20: Hyperparameters for LLaMA-350M and LLaMA-1B On C4 Dataset. The boldfaced values denote the optimal hyperparameters. LLaMA-350MLLaMA-1B OptimizerLRWDLRWD Adam1e-31e-56e-41e-6 AdamW 1e-30.16e-40.1 Muon0.0025, 0.005,0.01,0.0150.10.005,0.010.1 HTMuon0.0025, 0.005,0.01,0.0150.10.005,0.010.1 Table 21: Hyperparameters for ResNet18 and ResNet50 On CIFAR-100 Dataset. The boldfaced values denote the optimal hyperparameters. ResNet18ResNet50 OptimizerLRWDLRWD SGDM0.1,0.2,0.3,0.45e-40.1,0.2,0.3,0.45e-4 Muon0.015,0.02,0.025,0.030.10.015,0.02,0.025,0.030.1 NorMuon0.015,0.02,0.025,0.030.10.015,0.02,0.025,0.030.1 HTMuon0.015,0.02,0.025,0.030.10.015,0.02,0.025,0.030.1 HTMuon+NorMuon 0.015,0.02,0.025,0.030.10.015,0.02, 0.025, 0.030.1 Table 22: Hyperparameters for ResNet18 and ResNet50 On CIFAR-10 Dataset. The boldfaced values denote the optimal hyperparameters. ResNet18ResNet50 OptimizerLRWDLRWD SGDM0.1,0.2,0.3,0.45e-40.1,0.2,0.3,0.45e-4 Muon0.01,0.015,0.02,0.0250.10.01,0.015,0.02,0.0250.1 NorMuon 0.01,0.015,0.02,0.0250.10.01,0.015,0.02,0.0250.1 HTMuon 0.01,0.015,0.02,0.0250.10.01,0.015,0.02,0.0250.1 HTMuon+NorMuon0.01,0.015,0.02,0.0250.10.01,0.015,0.02,0.0250.1 33 Algorithm 8 NorMuon 1:Input: Initial weights W 0 ∈ R m×n , loss functionL, learning rateη, momentum parameters (β 1 ,β 2 ), perturbation parameter ε, weight decay λ. 2: Initialize momentum M 0 ∈ R m×n ← 0, second moment v 0 ∈ R m ← 0 3: for t = 1, 2,... do 4:G t ←∇ W L(W t ) 5:M t ← β 1 M t−1 + (1− β 1 )G t 6:O t ← NS5(M t ) 7:v t ← β 2 v t−1 + (1− β 2 ) mean cols (O t ⊙ O t ) 8:V t ← ExpandRows(v t )(V t ∈ R m×n ) 9: b O t ← O t ⊘ ( √ V t + ε) 10:ˆη ← 0.2η √ mn ∥ b O t ∥ F 11:W t+1 ← W t − ηλW t − ˆη b O t 12: end for Algorithm 9 AdaMuon 1:Input: Initial 2D-weights W 0 ∈ R n×m , loss functionL, learning rateη, weight decayλ, momentumβ, Newton–Schulz steps T , small constant ε. 2: Output: Updated weights W 3: Initialize first-momentum M 0 ← 0, second-momentum V 0 ← 0 4: for each iteration t = 1, 2,... do 5:Compute gradient: G t =∇ W t L(W t ) 6:Update first momentum: M t = β· M t−1 + G t 7:Compute sign-stabilized orthogonal direction: O t = Newton–Schulz(Sign(M t ), T ) 8:Update second momentum: V t = β· V t−1 + (1− β)· (O t ⊙ O t ) 9:Apply second momentum update: ˆ O t = O t ⊘ ( √ V t + ε· 1) 10:RMS-aligned: γ t = 0.2· √ mn ∥ ˆ O t ∥ F 11:Update weights: W t+1 = W t − η γ t ˆ O t + λW t 12: end for Algorithm 10 MARS 1: Hyperparameters: β 1 , β 2 , γ, ε, η, λ, g norm 2: State: m, v, g t−1 3: for t = 1, 2,... do 4: c t = g t + γ β 1 1−β 1 (g t − g t−1 ) 5:ˆc t = c t max1, g norm ∥c t ∥ 2 , 6: m t = β 1 m t−1 + (1− β 1 ) ˆc t 7: v t = β 2 v t−1 + (1− β 2 ) ˆc 2 t 8:ˆm t = m t 1−β t 1 ,ˆv t = v t 1−β t 2 9: w t+1 = w t − η ˆm t √ ˆv t +ε − η λw t . 10: end for 34 Algorithm 11 SOAP 1: Hyperparameters: β 1 , β 2 , μ, k, ε, blocksize,g norm 2: Partition all parameters into blocksize× blocksize 3: Update Rules for Each Block: 4: for t = 1, 2,... do 5:ˆg t = g t max1, g norm ∥g t ∥ 2 6:ˆg t = Q A ˆg t Q B 7: m t = β 1 m t−1 + (1− β 1 ) ˆg t , v t = β 2 v t−1 + (1− β 2 ) ˆg 2 t 8:ˆm t = m t 1−β t 1 , ˆv t = v t 1−β t 2 9: w t+1 = w t − η t Q ⊤ A ˆm t √ ˆv t +ε Q ⊤ B 10: G A = μG A + (1− μ) ˆg t ˆg ⊤ t + εI, G B = μG B + (1− μ) ˆg ⊤ t ˆg t + εI 11:if t mod k = 0 : Q A = QR(G A Q A ), Q B = QR(G B Q B ) 12: end for Algorithm 12 Cautious 1: Hyperparameters: β 1 , β 2 , ε, η, λ, g norm 2: State: m, v 3: for t = 1, 2,... do 4:ˆg t = g t max1, g norm ∥g t ∥ 2 5: m t = β 1 m t−1 + (1− β 1 ) ˆg t 6: v t = β 2 v t−1 + (1− β 2 ) ˆg 2 t 7:ˆm t = m t 1−β t 1 ,ˆv t = v t 1−β t 2 8: u t = ˆm t √ ˆv t +ε , s t = I u t · ˆg t > 0 9:ˆu t = u t ·s t mean(s t ) , w t+1 = w t − η ˆu t − η λw t 10: end for Algorithm 13 COSMOS 1:For anm× nlayerW, the algorithm maintain four matrices:U ∈ R n×r ,S ∈ R r×r ,V ∈ R m×r , and M ∈ R m×n per layer. 2: Input: Learning rateη, combination weightγ, projection rankr ≪ n, momentum parameters (β 1 ,β 2 ), perturbation parameter ε. For simplicity, we omit the initialization. 3: for t = 0,... do 4:Sample batch M t 5: G t ←∇ W φ M t (W t ) 6: M t ← β 1 M t−1 + (1− β 1 )G t 7: U t ← QR β 2 U t−1 S t−1 + (1− β 2 )G ⊤ t G t U t−1 8: S t ← U ⊤ t β 2 U t−1 S t−1 U ⊤ t−1 + (1− β 2 )G ⊤ t G t U t 9: V t ← β 2 V t−1 + (1− β 2 )(G t U t )⊙ (G t U t ) 10: A t = M t U t /(1− β t 1 ) p (V t + ε)/(1− β t 2 ) ! U ⊤ t 11: B t ← NORM NS5 M t − M t U t U ⊤ t ∥M t − M t U t U ⊤ t ∥ F 12: ̃ G t ← A t + γ· B t · √ m 13: W t+1 ← W t − η· NORM( ̃ G t )· √ m 14: end for 35 Algorithm 14 GaLore 1:For anm× nlayerW, the algorithm maintains a rank-rsubspace projector (left or right, depending on the shape) and Adam moments in the projected space. 2: Input: Learning rateη, scale factorα, projection rankr, subspace change frequencyT, Adam parameters (β 1 ,β 2 ), numerical constant ε. For simplicity, we omit the initialization. 3: for t = 0, 1,... do 4:Sample batch B t 5: G t ←∇ W φ B t (W t ) 6:if t mod T = 0 then 7:U, Σ,V ← SVD(G t ) 8:if m≤ n then 9:P t ← U [:,1:r] left projector 10:else 11:Q t ← V [:,1:r] right projector 12:end if 13:else 14:if m≤ n then 15:P t ← P t−1 reuse projector 16:else 17:Q t ← Q t−1 reuse projector 18:end if 19:end if 20:if m≤ n then 21:R t ← P ⊤ t G t project to compact space, R t ∈ R r×n 22:else 23:R t ← G t Q t project to compact space, R t ∈ R m×r 24:end if 25:UPDATE(R t ) by Adam in the projected space: 26: M t ← β 1 M t−1 + (1− β 1 )R t 27: V t ← β 2 V t−1 + (1− β 2 )(R t ⊙ R t ) 28: ˆ M t ← M t /(1− β t+1 1 ) 29: ˆ V t ← V t /(1− β t+1 2 ) 30: N t ← ˆ M t /( p ˆ V t + ε) 31:if m≤ n then 32: ̃ G t ← α· P t N t project back to original space 33:else 34: ̃ G t ← α· N t Q ⊤ t project back to original space 35:end if 36: W t+1 ← W t − η· ̃ G t 37: end for 36 Algorithm 15 Sophia 1:Input: Parametersθ 1 , learning ratesη t T t=1 , hyperparametersλ,γ,β 1 ,β 2 ,ε, curvature estimator choice Estimator∈Hutchinson, Gauss-Newton-Bartlett. 2: Set m 0 ← 0, v 0 ← 0, and h 1−k ← 0. 3: for t = 1,...,T do 4:Compute minibatch loss L t (θ t ). 5: g t ←∇L t (θ t ). 6: m t ← β 1 m t−1 + (1− β 1 )g t . 7:if t mod k = 1 then 8: ˆ h t ← Estimator(θ t ). 9:h t ← β 2 h t−k + (1− β 2 ) ˆ h t . 10:else 11:h t ← h t−1 . 12:end if 13: θ t ← θ t − η t λθ t weight decay 14: θ t+1 ← θ t − η t · clip m t maxγ·h t , ε , 1 . 15: end for 37