On the Convergence of Muon and Beyond

## Abstract

Abstract The Muon optimizer has demonstrated remarkable empirical success in handling matrix-structured parameters for training neural networks. However, a significant gap remains between its practical performance and theoretical understanding. Existing analyses show that the Muon variants achieve only a suboptimal iteration complexity of 𝒪 ​ ( T − 1 / 4 ) \mathcal{O}(T^{-1/4}) in stochastic non-convex settings, where T T denotes the number of iterations. To explore the theoretical limits of the Muon framework, we analyze two Momentum-based Variance-Reduced variants: a one-batch version (Muon-MVR1) and a two-batch version (Muon-MVR2). We provide the first rigorous proof that incorporating variance reduction enables Muon-MVR2 to attain the optimal iteration complexity of 𝒪 ~ ​ ( T − 1 / 3 ) \tilde{\mathcal{O}}(T^{-1/3}) , thereby matching the theoretical lower bound for this class of problems. Furthermore, our analysis establishes last-iterate convergence guarantees for Muon variants under the Polyak-Łojasiewicz (PŁ) condition. Extensive experiments on vision (CIFAR-10) and language (C4) benchmarks corroborate our theoretical findings on per-iteration convergence. Overall, this work offers the first proof of optimality for a Muon-style optimizer and clarifies the path toward developing more practically efficient, accelerated variants.

## 1 Introduction

The immense computational cost of pre-training Large Language Models (LLMs) has spurred a surge of research into novel optimization methods designed to enhance parameter efficiency and training stability hoffmann2022training; Liu2023SophiaAS; Chen2023SymbolicDO; vyas2025soap; pethick2025training; Yuan2024MARSUT. Among these, methods based on matrix orthogonalization have recently garnered significant attention from both academia and industry jordan2024muon; Liu2025MuonIS. In particular, the Muon optimizer has emerged as a notable milestone due to its impressive empirical performance Liu2025MuonIS; An2025ASGOAS; liu2025cosmos; shah2025practical. However, despite its practical success, the theoretical understanding of Muon’s underlying mechanisms has surprisingly lagged behind, with existing convergence analyses being fraught with limitations and even critical fallacies.

Specifically, the current theoretical exploration of Muon’s convergence faces three primary obstacles. First, existing analyses have not established convergence to a stationary point without reliance on problem dimension or batch size; the available bounds become valid only when the batch size is sufficiently large, so the basic requirement for a trustworthy optimizer remains unfulfilled Sato2025ConvergenceBA; sfyraki2025lions. Second, some analyses are predicated on flawed mathematical assumptions, such as the erroneous use of the inequality $\|\mathbf{S}^{-1}\|_{2}\leq 1/\|\mathbf{S}\|_{2}$ in a key part of their proof, which casts serious doubt on the validity of their convergence claims Li2025ANO. Finally, the most rigorous existing work analyzes both the standard Muon algorithm (like Algorithm [1](https://arxiv.org/html/2509.15816v3#alg1), Option MVR1 with $\gamma=0$) and its simplified variant , considering their convergence to a non-standard $\epsilon$-nuclear norm stationary point. This leaves the behavior of the Nesterov-Accelerated Muon algorithm (like Algorithm [1](https://arxiv.org/html/2509.15816v3#alg1), Option MVR1), as discussed in Liu2025MuonIS; Sato2025ConvergenceBA, and Variance-Reduction Muon (like Algorithm [1](https://arxiv.org/html/2509.15816v3#alg1), Option MVR2) under standard settings as an open question shen2025convergence.

To bridge this critical theoretical gap, this work establishes a rigorous and complete theoretical foundation for the Muon optimizer. Our main contributions are threefold:

- •
For general non-convex settings, we provide a convergence analysis for both the standard Muon (Algorithm [1](https://arxiv.org/html/2509.15816v3#alg1), Option MVR1 with $\gamma=0$) and the Muon-MVR1 algorithm (Algorithm [1](https://arxiv.org/html/2509.15816v3#alg1), Option MVR1).
- •
Regarding the optimal iteration complexity, we are the first to demonstrate that the Variance-Reduction Muon-MVR2 algorithm (Algorithm [1](https://arxiv.org/html/2509.15816v3#alg1), Option MVR2) can overcome the bottleneck of traditional momentum methods, achieving the optimal iteration complexity of $\tilde{\mathcal{O}}(T^{-1/3})$.
- •
Under the Polyak-Łojasiewicz (PL) condition, we prove that our proposed algorithms, Muon-MVR1, and Muon-MVR2, all achieve sublinear convergence rates. More precisely, we demonstrate that Muon-MVR1 converge at a rate of $\tilde{\mathcal{O}}(T^{-1/2})$, and Muon-MVR2 converges at an accelerated rate of $\tilde{\mathcal{O}}(T^{-2/3})$

Table [1](https://arxiv.org/html/2509.15816v3#S1.T1) summarizes the main contributions of our work and compares them with existing methods.

**Table 1: Comparison of Existing Muon-type Analyses with Ours.**
|  | Smooth^a | Stoc. Gradient<br>Estimator ^b | Batch<br>Size | Iteration<br>Complexity | Last-iterate<br>Conv. Rate |
| --- | --- | --- | --- | --- | --- |
| Li2025ANO | $L$ | MVR1($\gamma=0$) | $\mathcal{O}(1)$ | $\mathcal{O}(T^{-1/4})$ | ✗ |
| Sato2025ConvergenceBA | $L$ | MVR1($\gamma=0$) | $\mathcal{O}(\epsilon^{-1})$ | $\mathcal{O}(T^{-1})+\mathcal{O}(1)^{c}$ | ✗ |
| Sato2025ConvergenceBA | $L$ | MVR1 | $\mathcal{O}(\epsilon^{-1})$ | $\mathcal{O}(T^{-1})+\mathcal{O}(1)^{c}$ | ✗ |
| sfyraki2025lions | $L_{+}$ | MVR2 | $\mathcal{O}(\epsilon^{-1})^{d}$ | $\mathcal{O}(T^{-1/3})$ | ✗ |
| shen2025convergence | $L$ | MVR1($\gamma=0$) | $\mathcal{O}(1)$ | $\mathcal{O}(T^{-1/4})$ | ✗ |
| Ours | $L$ | MVR1($\gamma=0$) | $\mathcal{O}(1)$ | $\tilde{\mathcal{O}}(T^{-1/4})$ | $\tilde{\mathcal{O}}(T^{-1/2})$ |
| Ours | $L$ | MVR1 | $\mathcal{O}(1)$ | $\tilde{\mathcal{O}}(T^{-1/4})$ | $\tilde{\mathcal{O}}(T^{-1/2})$ |
| Ours | $L_{+}$ | MVR2 | $\mathcal{O}(1)$ | $\tilde{\mathcal{O}}(T^{-1/3})$ | $\tilde{\mathcal{O}}(T^{-2/3})$ |

$\bullet$ Organizations. The rest of the paper is organized as follows. Section [2](https://arxiv.org/html/2509.15816v3#S2) reviews existing variants of the Muon algorithm. Section [3](https://arxiv.org/html/2509.15816v3#S3) addresses iteration complexity and last-iterate convergence rate. Section [4](https://arxiv.org/html/2509.15816v3#S4) presents experimental results demonstrating the effectiveness of our method. Related work is reviewed in Section [5](https://arxiv.org/html/2509.15816v3#S5), and conclusions are drawn in Section [6](https://arxiv.org/html/2509.15816v3#S6).

$\bullet$ Notations. We denote scalars by non-bold letters (e.g., $a,A$), vectors in $\mathbb{R}^{d}$ by bold lowercase letters (e.g., $\mathbf{a}$), and matrices by bold uppercase letters (e.g., $\mathbf{A}$). The space $\mathbb{R}^{d}$ is endowed with the Euclidean inner product $\langle\mathbf{x},\mathbf{y}\rangle_{2}:=\mathbf{x}^{\top}\mathbf{y}$ and norm $\|\mathbf{x}\|_{2}$. For matrices, we employ the Frobenius inner product $\langle\mathbf{A},\mathbf{B}\rangle_{F}:=\text{tr}(\mathbf{A}^{\top}\mathbf{B})$ and the corresponding norm $\|\mathbf{A}\|_{F}$. The nuclear norm, denoted by $\|\mathbf{A}\|_{*}$, is defined as the sum of the singular values of the matrix, $\|\mathbf{A}\|_{*}=\sum_{i}\sigma_{i}(\mathbf{A})$. Throughout the paper, $[m]$ denotes the set of integers $\{1,2,\ldots,m\}$, and $\mathbb{N}$ denotes the set of non-negative integers. The model is parameterized by a matrix $\mathbf{X}\in\mathbb{R}^{m\times n}$. Without loss of generality, we assume $m\geq n$, so the rank of the matrix is at most $n$. The model is optimized by minimizing the empirical loss function $f(\mathbf{X}):=\frac{1}{N}\sum_{i\in[N]}f_{i}(\mathbf{X})$, where $N$ is the number of training data points, and $f_{i}(\mathbf{X})$ is the loss function for $\mathbf{X}\in\mathbb{R}^{m\times n}$ with respect to the $i$-th training data point $\mathbf{z}_{i}$ (for $i\in[N]$). Let $\xi$ be a random variable that is independent of $\mathbf{X}\in\mathbb{R}^{m\times n}$, and let $\mathbb{E}_{\xi}[\mathbf{X}]$ denote the expectation of a random variable $\mathbf{X}$ with respect to $\xi$.

## 2 Revisiting the Muon Algorithms

We consider the following optimization problem:

$$ $\min_{\mathbf{X}\in\mathbb{R}^{m\times n}}f(\mathbf{X}),\mathrm{~where~}f(\mathbf{X})=\mathbb{E}_{\xi\sim\mathcal{D}}[f(\mathbf{X};\xi)],$ (1) $$

where $f:\mathbb{R}^{m\times n}\rightarrow\mathbb{R}$ is the loss function, $\mathbf{X}$ denotes the decision variable, and $\xi$ represents a random variable(e.g., a training data sample) drawn from an unknown distribution $\mathcal{D}$. We assume that $f$ is differentiable and possibly nonconvex.

The Muon optimizer begins by computing a momentum-based variance-reduced gradient update, similar in spirit to ADAM Kingma2014AdamAM, STORM Cutkosky2019MomentumBasedVR, and SGD with Nesterov momentum sutskever2013importance. The momentum term $\mathbf{M}_{t}$ is then projected onto the set of orthogonal matrices. This orthogonalization step equalizes the singular values, ensuring that no principal component direction dominates the optimization. Finally, the resulting scaled orthogonal matrix is used to update the model parameters. The Muon algorithm is summarized in Algorithm [1](https://arxiv.org/html/2509.15816v3#alg1).

Figure: Algorithm 1 Muon

This orthogonalization step can be obtained from the Singular Value Decomposition (SVD) of $\mathbf{M}_{t}$:

$$ $\mathbf{M}_{t}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top},\qquad\mathbf{O}_{t}=\mathbf{U}\mathbf{V}^{\top}.$ $$

Equivalently, we have $\mathbf{O}_{t}=(\mathbf{M}_{t}\mathbf{M}_{t}^{\top})^{-1/2}\mathbf{M}_{t}$, which shows that orthogonalization reduces to computing an inverse square root rather than a full SVD. Since forming and decomposing $\mathbf{M}_{t}\mathbf{M}_{t}^{\top}$ remains costly, the popular Muon implementation jordan2024muon uses the quintic Newton–Schulz iteration to approximate the inverse square root. This recurrence converges in only a few steps (typically five), producing a numerically stable, rank-preserving orthogonalization of $\mathbf{M}_{t}$ that is nearly as accurate as SVD but far more efficient.

Algorithm [1](https://arxiv.org/html/2509.15816v3#alg1) incorporates two distinct strategies for momentum-based variance reduction, termed Muon-MVR1 and Muon-MVR2. These options present a fundamental trade-off between computational efficiency and theoretical rigor. MVR2 implements a principled variance reduction scheme at the cost of two gradient evaluations per step, while MVR1 serves as a computationally cheaper, single-gradient approximation. We detail both below.

$\blacktriangleright$ Option 1: Muon-MVR1 (One-batch Approximation). The first option, MVR1, augments the classical momentum update with a variance-reducing term that reuses the gradient from the previous step:

$$ $\displaystyle\mathbf{M}_{t}=\beta_{t}\mathbf{M}_{t-1}+(1-\beta_{t})\nabla f(\mathbf{X}_{t};\xi_{t})+\gamma\cdot\beta_{t}\cdot(\nabla f(\mathbf{X}_{t};\xi_{t})-\nabla f(\mathbf{X}_{t-1};\xi_{t-1})).$ (2) $$

The primary advantage of this formulation is its computational efficiency, as it requires only one stochastic gradient evaluation per iteration. This update rule is flexible:

(i) When $\gamma=0$, Rule ([2](https://arxiv.org/html/2509.15816v3#S2.E2)) degenerates to the standard exponential moving average (EMA) of gradients (Rule ([3](https://arxiv.org/html/2509.15816v3#S2.E3))), a stochastic gradient estimator widely used in optimizers like Adam Kingma2014AdamAM.

$$ $\mathbf{M}_{t}=\beta_{t}\mathbf{M}_{t-1}+(1-\beta_{t})\nabla f(\mathbf{X}_{t};\xi_{t}).$ (3) $$

(ii) With the specific parameterization ($\beta_{t}=\mu,\gamma=1-\beta_{t}$), it can be shown that Rule ([2](https://arxiv.org/html/2509.15816v3#S2.E2)) becomes mathematically equivalent to the update scheme of the Muon optimizer jordan2024muon; Liu2025MuonIS, given by:

$$ $\displaystyle\mathbf{C}_{t}$ $\displaystyle=\mu\mathbf{C}_{t-1}+\nabla f(\mathbf{X}_{t};\xi_{t}),$ (4) $\displaystyle\mathbf{M}_{t}$ $\displaystyle=\mu\mathbf{C}_{t}+\nabla f(\mathbf{X}_{t};\xi_{t}).$ $$

$\blacktriangleright$ Option 2: Muon-MVR2 (Two-batch Principled VR). The second option, MVR2, incorporates a more rigorous variance-reduction mechanism inspired by methods like SPIDER fang2018spider, STORM Cutkosky2019MomentumBasedVR, SVRG zhou2020stochastic, SUPER-Adam Huang2021SUPERADAMFA, and MARS Yuan2024MARSUT:

$$ $\displaystyle\mathbf{M}_{t}=\beta_{t}\mathbf{M}_{t-1}+(1-\beta_{t})\nabla f(\mathbf{X}_{t};\xi_{t})+\gamma\cdot\beta_{t}\cdot(\nabla f(\mathbf{X}_{t};\xi_{t})-\nabla f(\mathbf{X}_{t-1};\xi_{t})).$ (5) $$

The key distinction from MVR1 is the correction term. MVR2 subtracts the gradient computed on the previous parameters but with the current mini-batch, i.e., $\nabla f(\mathbf{X}_{t-1};\xi_{t})$. This modification is crucial as it is designed to directly cancel the variance introduced by the mini-batch $\xi_{t}$ Cutkosky2019MomentumBasedVR; Huang2021SUPERADAMFA; Yuan2024MARSUT. However, this theoretical benefit comes at the cost of requiring two gradient evaluations per step. MVR1 can be formally understood as a practical approximation of MVR2. The difference between their update rules is a single noise term, $\Delta^{\text{Noise}}_{t-1}=\nabla f(\mathbf{X}_{t-1};\xi_{t})-\nabla f(\mathbf{X}_{t-1};\xi_{t-1})$. Under the standard assumption of bounded variance ($\mathbb{E}_{\xi}[\|\nabla f(\mathbf{X};\xi)-\nabla f(\mathbf{X})\|_{F}^{2}]\leq\sigma^{2}$), the variance of this noise is well-controlled, satisfying $\mathbb{E}[\|\Delta^{\text{Noise}}_{t-1}\|_{F}^{2}]\leq 2\sigma^{2}$. While MVR1 is often sufficient in practice, this structural difference leads to fundamentally different theoretical guarantees. The principled variance cancellation in MVR2 allows our algorithm to achieve a superior iteration complexity of $\tilde{\mathcal{O}}(T^{-1/3})$, as we will formally establish in Theorem [3.2](https://arxiv.org/html/2509.15816v3#S3.Thmtheorem2).

## 3 Convergence Analysis

We begin our analysis by situating it within the established context of first-order stochastic optimization. The $\mathcal{O}(T^{-1/4})$ iteration complexity is a well-known bottleneck for methods like SGD, and recent work by shen2025convergence confirmed this limitation for the standard Muon algorithm. Accordingly, our first step in this section is to analyze the Muon-MVR1 variant and formally reaffirm this baseline iteration complexity in Theorem [3.1](https://arxiv.org/html/2509.15816v3#S3.Thmtheorem1). While the rate itself is standard, our primary contribution here lies in the novel and unified analytical framework we develop to prove it, as this framework will be instrumental for subsequent results.

Our analysis is conducted under a diminishing stepsize schedule, a standard setting that guarantees convergence to an optimal solution. Nevertheless, our theoretical framework is general enough to encompass the constant-stepsize setting, where it yields a tighter convergence bound devoid of the logarithmic factor $\ln(T)$.

To facilitate the analysis of convergence for Muon, we make the following assumptions:

###### Assumption 3.1 .

The function $f$ is bounded from below. There exists $f^{*}>-\infty$ such that $f(\mathbf{X})\geq f^{*}$, for all $\mathbf{X}\in\mathbb{R}^{m\times n}$.

###### Assumption 3.2 .

The function $f$ is $L$-smooth: $\|\nabla f(\mathbf{Y})-\nabla f(\mathbf{X})\|_{F}\leq L\|\mathbf{Y}-\mathbf{X}\|_{F}$.

###### Assumption 3.3 .

The function $f$ is $L$-smooth for any $\xi$: $\|\nabla f(\mathbf{Y};\xi)-\nabla f(\mathbf{X};\xi)\|_{F}\leq L\|\mathbf{Y}-\mathbf{X}\|_{F}$.

###### Assumption 3.4 .

The variance of unbiased stochastic gradient is finite. Specifically, there exists a constant $\sigma>0$ such that for all $\mathbf{X}\in\mathbb{R}^{m\times n}$, the following holds: $\mathbb{E}[\nabla f(\mathbf{X};\xi)]=\nabla f(\mathbf{X})$ and $\mathbb{E}\|\nabla f(\mathbf{X};\xi)-\nabla f(\mathbf{X})\|_{F}^{2}\leq\sigma^{2}$.

These assumptions are quite common Zhou2018OnTC; Chen2018OnTC; Huang2021SUPERADAMFA; Guo2021ANC; Li2023ConvergenceOA; Wang2023ClosingTG; 10586270; Yuan2024MARSUT.

###### Assumption 3.1 .

###### Assumption 3.2 .

###### Assumption 3.3 .

###### Assumption 3.4 .

### 3.1 Erogdic Convergence of Muon

In this subsection, we establish the ergodic convergence of Muon-MVR1 and Muon-MVR2.

#### 3.1.1 Option MVR1

We begin our analysis with Option MVR1, a straightforward implementation of one-batch momentum-based variance reduction. The following theorem establishes its ergodic convergence rate, demonstrating that the algorithm converges to a stationary point at a rate of $\tilde{\mathcal{O}}(T^{-1/4})$ for specific choices of learning rate and momentum schedules.

###### Theorem 3.1 .

Suppose Assumptions [3.1](https://arxiv.org/html/2509.15816v3#S3.Thmassumption1), [3.2](https://arxiv.org/html/2509.15816v3#S3.Thmassumption2), and [3.4](https://arxiv.org/html/2509.15816v3#S3.Thmassumption4) hold. Consider Algorithm [1](https://arxiv.org/html/2509.15816v3#alg1) with a learning rate of $\eta_{t}=t^{-3/4}$. The expected average squared norm of the gradient is bounded for the following options:

- 1.
For Option MVR1 ($\gamma=0$), with momentum parameter $\beta_{t}=1-t^{-1/2}$, the bound is given by:
$\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\|\nabla f(\mathbf{X}_{t})\|_{F}\leq\frac{f(\mathbf{X}_{1})-f^{*}}{T^{1/4}}+\frac{A_{1}\ln T+A_{2}}{T^{1/4}},$
where $A_{1}=2L^{-1}\sigma^{2}+4\sqrt{2}Ln+Ln+L/2$ and $A_{2}=4L^{-1}\sigma^{2}+4\sqrt{2}Ln+Ln+L/2$.
- 2.
For Option MVR1 ($\gamma_{t}=t^{-1/2}$), with momentum parameter $\beta_{t}=1-(t+1)^{-1/2}$, the bound is given by:
$\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\|\nabla f(\mathbf{X}_{t})\|_{F}\leq\frac{f(\mathbf{X}_{1})-f^{*}}{T^{1/4}}+\frac{A_{1}\ln T+A_{2}}{T^{1/4}},$
where $A_{1}=4L^{-1}\sigma^{2}+8\sqrt{2}Ln+Ln+L/2$ and $A_{2}=10L^{-1}\sigma^{2}+8\sqrt{2}Ln+Ln+L/2$.

See Appendix [B](https://arxiv.org/html/2509.15816v3#A2) for details.

###### Remark 3.1 .

Theorems [3.1](https://arxiv.org/html/2509.15816v3#S3.Thmtheorem1) establishes that as $T\to\infty$, the leading terms diminish to zero, and the algorithm converges to a neighborhood of a stationary point. The size of this neighborhood is determined by the stochastic gradient variance, the learning rate, and the momentum parameter. Notably, this convergence bound is free of any non-vanishing additive error term that depends on the dimension $n$ Sato2025ConvergenceBA. We absorb the dimensional dependence into the $\mathcal{O}$-notation to define $\tilde{\mathcal{O}}$, which yields the following complexity for the algorithm:

$$ $\min_{t=1,\dots,T}\mathbb{E}\|\nabla f(\mathbf{X}_{t})\|_{F}\leq\tilde{\mathcal{O}}(T^{-1/4}).$ $$

###### Theorem 3.1 .

###### Remark 3.1 .

#### 3.1.2 Option MVR2

While prior work has established a iteration complexity of $\tilde{\mathcal{O}}(T^{-1/3})$ for variance-reduction in non-convex settings for various algorithms, such as SGD fang2018spider; Cutkosky2019MomentumBasedVR; zhou2020stochastic and Adam Huang2021SUPERADAMFA; Yuan2024MARSUT, the theoretical convergence properties of this technique when applied to the Muon optimizer have remained an open question. We bridge this theoretical gap in Theorem [3.2](https://arxiv.org/html/2509.15816v3#S3.Thmtheorem2), which rigorously proves that Option MVR2 of Algorithm [1](https://arxiv.org/html/2509.15816v3#alg1) (Muon-MVR2) achieves the same $\tilde{\mathcal{O}}(T^{-1/3})$ iteration complexity in the general non-convex setting.

###### Theorem 3.2 .

Under Assumptions [3.1](https://arxiv.org/html/2509.15816v3#S3.Thmassumption1), [3.3](https://arxiv.org/html/2509.15816v3#S3.Thmassumption3), and [3.4](https://arxiv.org/html/2509.15816v3#S3.Thmassumption4), consider Algorithm [1](https://arxiv.org/html/2509.15816v3#alg1) with Option MVR2 and update rule ([5](https://arxiv.org/html/2509.15816v3#S2.E5)) parameters set as $\beta_{t}=1-\eta_{t}$, $\eta_{t}=t^{-2/3}$ for $t\geq 1$, and $\gamma=1$. Then, the following bounds hold:

$$ $\displaystyle\begin{aligned} \frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\|\nabla f(\mathbf{X}_{t})\|_{F}&\leq\frac{\sqrt{4\sigma^{2}+(16L^{2}n+8\sigma^{2})(1+\ln T)}}{T^{1/3}}+\frac{G}{T^{1/3}},\end{aligned}$ $$

where

$$ $\displaystyle G=f(\mathbf{X}_{1})-f^{*}+2\sigma^{2}L^{-1}+2Ln+\left(8Ln+4\sigma^{2}L^{-1}+\frac{Ln}{2}\right)(1+\ln T).$ $$

See Appendix [D](https://arxiv.org/html/2509.15816v3#A4) for details.

###### Remark 3.2 .

Theorem [3.2](https://arxiv.org/html/2509.15816v3#S3.Thmtheorem2) establishes a significant theoretical advantage. In contrast to prior work sfyraki2025lions that requires a computationally expensive growing batch size ($b=\mathcal{O}(T^{1/3})$) to control gradient variance, our result shows that Muon-MVR2 achieves the same $\tilde{\mathcal{O}}(T^{-1/3})$ iteration complexity in a more practical, fixed mini-batch setting. This is made possible by the specific hyperparameter schedule, where the learning rate $\eta_{t}=t^{-2/3}$ and the momentum parameter $\beta_{t}=1-\eta_{t}$ are critically coupled. This schedule ensures a careful balance between optimization progress and noise suppression.

###### Remark 3.3 .

As noted in Yuan2024MARSUT, a more sophisticated, adaptive setting for the variance-reduction parameter $\gamma$ can be employed. Specifically, by setting $\gamma=\gamma_{t}=1-\frac{A_{t}}{\beta_{t}}$, where $A_{t}$ is defined in Lemma [C.2](https://arxiv.org/html/2509.15816v3#A3.Thmlemma2), a key term in the analysis becomes $P_{t}=-\mathbb{E}\|\nabla f(\mathbf{X}_{t};\xi_{t})-\nabla f(\mathbf{X}_{t-1};\xi_{t})\|_{F}^{2}\cdot A_{t}^{2}$, which is strictly negative. This leads to a tighter convergence bound, as it effectively introduces an additional beneficial term into the recurrence. However, computing this adaptive $\gamma_{t}$ is often impractical as it depends on quantities that are difficult to estimate during training. Consequently, we adhere to the common and more practical approach of using a constant $\gamma\leq 1$.

###### Theorem 3.2 .

###### Remark 3.2 .

###### Remark 3.3 .

### 3.2 Non-Ergodic Convergence of Muon

In this subsection, we examine the performance of Muon-MVR1 and Muon-MVR2 under the setting where the non-convex objective functions satisfy the Polyak-Łojasiewicz (PL) condition. Our analysis is based on the following additional assumption:

###### Assumption 3.5 .

We assume the function $f$ is $\mu$-PL, i.e., $\|\nabla f(\mathbf{X}_{t})\|_{F}^{2}\geq 2\mu(f(\mathbf{X}_{t})-f^{*})$.

###### Remark 3.4 .

The PL condition has been widely employed in the convergence analysis of various first-order algorithms karimi2016linear; xie2020linear; li2022simple, though typically under restricted settings. Note that if $f(\cdot)$ is strongly convex, then it is necessarily convex and satisfies the PL condition. However, the converse does not hold in general; a counterexample is given by $f(x)=x^{2}+3\sin(2x)$.

###### Theorem 3.3 .

Suppose Assumptions [3.1](https://arxiv.org/html/2509.15816v3#S3.Thmassumption1), [3.2](https://arxiv.org/html/2509.15816v3#S3.Thmassumption2), [3.4](https://arxiv.org/html/2509.15816v3#S3.Thmassumption4), and [3.5](https://arxiv.org/html/2509.15816v3#S3.Thmassumption5) hold. Let $\{\mathbf{X}_{t}\}$ be the sequence of iterates generated by Algorithm [1](https://arxiv.org/html/2509.15816v3#alg1) with a step size of $\eta_{t}=t^{-3/4}$. We analyze the following two MVR1 schemes:

- 1.
Scheme 1 ($\gamma=0$): Using the update rule ([3](https://arxiv.org/html/2509.15816v3#S2.E3)) with momentum $\beta_{t}=1-t^{-1/2}$. For this case, we define the constant $\mathcal{A}_{1}$ as:
$\mathcal{A}_{1}=6L^{-1}\sigma^{2}+(8\sqrt{2}+2)Ln+L.$
- 2.
Scheme 2 ($\gamma_{t}=t^{-1/2}$): Using the update rule ([2](https://arxiv.org/html/2509.15816v3#S2.E2)) with parameters $\beta_{t}=1-(t+1)^{-1/2}$. For this case, we define the constant $\mathcal{A}_{2}$ as:
$\mathcal{A}_{2}=14L^{-1}\sigma^{2}+(16\sqrt{2}+2)Ln+L.$

Then, for either scheme $i\in\{1,2\}$, there exists a constant $T_{0}=2e^{4}$ such that for all iterations $T\geq T_{0}$, the expected suboptimality gap is bounded as follows:

$$ $\mathbb{E}[f(\mathbf{X}_{T+1})]-f^{*}\leq\frac{\mathcal{A}_{i}^{2}}{8\mu}\cdot\frac{(\ln T)^{2}}{T^{1/2}}.$ $$

###### Theorem 3.4 .

Suppose Assumptions [3.1](https://arxiv.org/html/2509.15816v3#S3.Thmassumption1), [3.2](https://arxiv.org/html/2509.15816v3#S3.Thmassumption2), [3.4](https://arxiv.org/html/2509.15816v3#S3.Thmassumption4), and [3.5](https://arxiv.org/html/2509.15816v3#S3.Thmassumption5) hold. Let $\{\mathbf{X}_{t}\}$ be the sequence of iterates generated by Algorithm [1](https://arxiv.org/html/2509.15816v3#alg1) with Option MVR2 (see update rule ([5](https://arxiv.org/html/2509.15816v3#S2.E5))). By setting the parameters to $\eta_{t}=t^{-2/3}$, $\beta_{t}=1-t^{-2/3}$, and $\gamma=1$, we define the constant $\mathcal{A}_{3}$ as :

$$ $\mathcal{A}_{3}=20L^{-1}\sigma^{2}+34Ln+L.$ $$

Then, there exists a constant $T_{0}=2e^{3}$ such that for all iterations $T\geq T_{0}$, the expected suboptimality gap is bounded as follows:

$$ $\mathbb{E}[f(\mathbf{X}_{T+1})]-f^{*}\leq\frac{2\mathcal{A}_{3}^{2}}{9\mu}\cdot\frac{(\ln T)^{2}}{T^{2/3}}.$ $$

Detailed proofs for Theorems [3.3](https://arxiv.org/html/2509.15816v3#S3.Thmtheorem3) and [3.4](https://arxiv.org/html/2509.15816v3#S3.Thmtheorem4) are provided in Appendix [F](https://arxiv.org/html/2509.15816v3#A6) and [G](https://arxiv.org/html/2509.15816v3#A7), respectively.

###### Remark 3.5 .

Our non-ergodic convergence proofs are unified by a recursive inequality on the potential function $\Delta_{t}=\mathbb{E}[f(\mathbf{X}_{t})]-f^{*}$:

$$ $\Delta_{t+1}\leq\Delta_{t}-\frac{\sqrt{2\mu}}{t^{p}}\sqrt{\Delta_{t}}+\Gamma_{t},$ $$

where $\sum_{i=1}^{t}\Gamma_{t}\leq\mathcal{O}(\ln t)$. The exponent $p\in(0,1)$ dictates the final convergence rate. For Theorem [3.3](https://arxiv.org/html/2509.15816v3#S3.Thmtheorem3), we establish this inequality with $p=3/4$, which leads to a rate of $\tilde{\mathcal{O}}(T^{-1/2})$. The superior variance control in Theorem [3.4](https://arxiv.org/html/2509.15816v3#S3.Thmtheorem4) enables a tighter analysis with $p=2/3$, which in turn yields the accelerated convergence rate of $\tilde{\mathcal{O}}(T^{-2/3})$.

###### Assumption 3.5 .

###### Remark 3.4 .

###### Theorem 3.3 .

###### Theorem 3.4 .

###### Remark 3.5 .

## 4 Experiments

In this section, we evaluate the performance of the Muon-variant optimizers on pretraining tasks. All experiments were conducted using four Ascend 910C (64GB) NPUs and four NVIDIA RTX 4090 (24GB) GPUs. Detailed experimental settings are provided in Appendix [H](https://arxiv.org/html/2509.15816v3#A8).

Figure: (a) Train and test accuracy curves
Refer to caption: x1.png

$\blacktriangleright$ ResNet18 on CIFAR10 Dataset. We evaluate our proposed methods on the CIFAR-10 image classification benchmark. We use a ResNet-18 he2016deep architecture and train for 100 epochs with a batch size of 128 and compare Muon-VR1 ($\gamma=0$), Muon-MVR1, and Muon-MVR2 against the widely-used baselines SGD and Adam. For each optimizer, the learning rate is selected via a grid search over the set $\{1\times 10^{-4},5\times 10^{-4},10^{-3},5\times 10^{-3},10^{-2},5\times 10^{-2},10^{-1}\}$. To ensure a robust comparison, all experiments are repeated over five different random seeds, and we report the mean results with one standard deviation shaded.
The training dynamics and performance are presented in Figure [1a](https://arxiv.org/html/2509.15816v3#S4.F1.sf1) and Figure [1b](https://arxiv.org/html/2509.15816v3#S4.F1.sf2). Muon-variant optimizers exhibit faster initial convergence in terms of training loss compared to SGD and Adam. More importantly, all Muon variants achieve lower final test error than the baselines. Among our proposed methods, Muon-MVR2 demonstrates the best performance, consistently securing the lowest test error and outperforming both the standard Muon-MVR1 ($\gamma=0$) and Muon-MVR1.

$\blacktriangleright$ LLaMA2 on C4 Dataset. To evaluate our optimizer in the context of large language model pre-training, we conduct experiments on the C4(Colossal Clean Crawled Corpus) dataset raffel2020exploring. We train a LLaMA2-130M model Touvron2023Llama2O from scratch and compare Muon-MVR1 ($\gamma=0$), Muon-MVR1, and Muon-MVR2 against the widely-used baseline AdamW loshchilov2017decoupled optimizer. For a fair and direct comparison, a fixed learning rate of $1\times 10^{-3}$ was used across all optimizers. All models were trained for 10,000 steps with a global batch size of 512. Performance is primarily assessed by tracking the training and validation loss throughout the process.
The training and validation loss curves are depicted in Figure [2](https://arxiv.org/html/2509.15816v3#S4.F2). The results indicate that Muon-MVR2 achieves a faster decrease and a lower final value in training loss compared to AdamW, Muon-MVR1 ($\gamma=0$), and Muon-MVR1. More critically, this improvement translates to better generalization performance; Muon-MVR2 consistently maintains a lower validation loss than the other Muon variants and demonstrates a significant advantage over the AdamW baseline. These findings suggest that the variance reduction component is particularly effective for stabilizing and accelerating convergence in large-scale language model training.

Figure: (a) Train loss curve
Refer to caption: x4.png

###### Remark 4.1 .

The precise, Variance-Reduction formulation ([5](https://arxiv.org/html/2509.15816v3#S2.E5)) is more rigorously aligned with stochastic optimization theory and, as demonstrated in our experiments, consistently achieves the highest validation accuracy. However, its computational overhead can be prohibitive in large-scale settings, as shown in Figure [1c](https://arxiv.org/html/2509.15816v3#S4.F1.sf3). We note that prior analysis in sfyraki2025lions; Sato2025ConvergenceBA suggests that increasing batch size can mitigate the impact of stochastic noise on convergence properties. Given that mini-batch training is the standard paradigm in practice, a trade-off between theoretical exactness and computational feasibility is often inherent. Therefore, echoing the sentiment in Yuan2024MARSUT, we posit that the approximate version serves as a pragmatic alternative when computational efficiency is a primary concern, as it incurs only a minimal degradation in performance. Conversely, in scenarios where maximizing validation accuracy is paramount, the precise formulation remains the recommended approach.

###### Remark 4.1 .

## 5 Related Work

### 5.1 The Evolution of Optimization Algorithms

The advancement of deep learning relies on first-order optimization methods. Nesterov1983AMF proposed the Momentum method, which leverages historical gradients to accelerate convergence. Subsequently, hinton2012lecture introduced RMSprop, enabling per-parameter adaptive learning rates. Kingma2014AdamAM then integrated these ideas in Adam, an optimizer that adapts using the first and second moments of the gradients, establishing it as a standard choice in the field.

Numerous variants have been proposed to address the limitations of Adam.  Reddi2018OnTC introduced AMSGrad to ensure a non-increasing learning rate, while  dozat2016incorporating created NAdam by incorporating Nesterov momentum. To improve regularization,  loshchilov2017decoupled developed AdamW, which decouples weight decay from the optimization step. Other methods focus on controlling the learning rate and variance:  Luo2019AdaptiveGM proposed AdaBound to dynamically clip learning rates;  Liu2019OnTV designed RAdam to rectify variance estimates in early training; and  Zhuang2020AdaBeliefOA created AdaBelief to adapt step sizes based on gradient belief. Recent developments continue this trend, such as  Xie2020AdaptiveID which decouples adaptation from momentum in the Adai framework, and  10586270 which introduced Adan with a novel Nesterov momentum estimation. The latest works from  Liang2024CautiousOI and  Yuan2024MARSUT further enhance efficiency with new masking and variance-reduction strategies, respectively. Additionally, the MGUP strategy proposed by chang2025mgup smooths zero masks into small non-zero values to alleviate potential non-convergence issues that the Cautious mask may cause in Adam.

Beyond Adam variants, research has explored other paradigms such as preconditioning methods that use parameter curvature. Gupta2018ShampooPS pioneered this direction with the Shampoo optimizer. Building on this work,  jordan2024muon proposed Muon, which adapts to curvature by orthogonalizing gradient momentum. Subsequent variants emerged, such as AdaMuon si2025adamuon which adds element-wise adaptivity, and COSMOS liu2025cosmos which integrates ideas from SOAP vyas2025soap for large model training. While these methods showed practical benefits, they often lacked convergence proofs. To bridge this theoretical gap for LMO-based methods, Gluon riabinin2025gluon introduces a novel layer-wise smoothness assumption, providing convergence guarantees that align with the practical implementations of optimizers like Muon and Scion pethick2025training. In the literature, other related preconditioning methods include ASGO An2025ASGOAS, PolarGrad lau2025polargrad, and AdaGO Zhang2025AdaGradMM, which introduces the Adagrad-Norm step size ward2020adagrad into a simplified version of Muon. Meanwhile, other high-performing optimizers not belonging to the Shampoo family also warrant attention, such as Sophia Liu2023SophiaAS, which improves second-moment estimation through efficient diagonal Hessian approximation and coordinate clipping, and Lion Chen2023SymbolicDO, a lightweight optimizer that only tracks momentum and uses the sign function to normalize updates.

### 5.2 Analysis of Muon

Recent theoretical analyses have clarified the mechanisms and convergence of the Muon optimizer. One central perspective frames modern optimizers as steepest descent methods under specific norm constraints; while this view does not directly account for momentum, it offers crucial insights into the update directions of optimizers like Muon and Adam bernstein2024old. Building on this, subsequent work has linked Muon with weight decay to the stochastic Frank-Wolfe method, identifying it as an instance operating under a spectral norm constraint sfyraki2025lions. From a constrained optimization viewpoint, Muon has also been characterized as a special case of Lion-$\mathcal{K}$, with its convergence to a KKT point proven in both deterministic and stochastic settings chen2025muon. These analyses, alongside broader research on norm-constrained stochastic conditional gradient methods, form Muon’s theoretical underpinnings pethick2025training. Beyond its formal framework, Muon’s implicit bias has become a key area of investigation. It has been demonstrated that Muon tends to converge toward max-margin solutions with respect to the spectral norm of the weight matrix, revealing a distinct implicit regularization and generalization preference compared to Adam fan2025implicit. A unifying preconditioning perspective, based on matrix polar decomposition, explains this behavior by distinguishing between curvature and gradient anisotropy. This decomposition clarifies the differing ways Muon and Adam handle various parameter types lau2025polargrad. The convergence analysis of Muon has also seen active development. Although rigorous proofs and iteration complexity in non-convex settings have been established under various smoothness assumptions shen2025convergence, the validity of some early results has been challenged. For instance, certain proofs were found to rely on incorrect mathematical inequalities, casting doubt on their conclusions Li2025ANO. Other analyses, while confirming the convergence of several Muon variants, report slow convergence rates or the need for stringent conditions to reach a stationary point Sato2025ConvergenceBA. On a practical level, especially for pre-training large-scale language models (LLMs), weight decay has been identified as an indispensable component for Muon. To enhance its scalability and utility, methods based on RMS analysis have also been proposed to effectively transfer learning rates from Adam to Muon Liu2025MuonIS.

## 6 Conclusion

In this work, we establish a rigorous theoretical foundation for the Muon optimizer, addressing the gap between its empirical success and formal analysis. We analyze two momentum-based variance-reduced variants of Muon: a one-batch version (Muon-MVR1) and a two-batch version (Muon-MVR2). We provide the first proof that Muon-MVR2 achieves the optimal iteration complexity of $\tilde{\mathcal{O}}(T^{-1/3})$ in stochastic nonconvex settings, matching the theoretical lower bound arjevani2023lower. Furthermore, we show that under the Polyak–Łojasiewicz condition, Muon-MVR1 and Muon-MVR2 attain sublinear nonergodic convergence rates of $\tilde{\mathcal{O}}(T^{-1/2})$ and $\tilde{\mathcal{O}}(T^{-2/3})$, respectively. Our theoretical results are validated by extensive experiments on the CIFAR-10 and C4 benchmarks, which confirm the practical acceleration and superior performance of Muon-MVR2 over its standard counterpart and other widely used optimizers. Overall, this research strengthens the Muon framework by providing robust theoretical guarantees and a practically effective new variant for deep learning training.

## Appendix

## Appendix A Lemmas for Theorem 3.1

### A.1 Lemma A.1

###### Lemma A.1 .

For Algorithm [1](https://arxiv.org/html/2509.15816v3#alg1), choosing an arbitrary parameter $\alpha>0$, we have the following inequality:

$$ $f(\mathbf{X}_{t+1})\leq f(\mathbf{X}_{t})-\eta_{t}\|\mathbf{M}_{t}\|_{F}+\frac{\eta_{t}\alpha}{2}\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}+\frac{\eta_{t}n}{2\alpha}+\frac{L\eta^{2}n}{2}.$ $$

###### Proof.

According to Assumption [3.2](https://arxiv.org/html/2509.15816v3#S3.Thmassumption2), we have the upper bound for the function value:

$$ $\displaystyle f(\mathbf{X}_{t+1})$ $\displaystyle\leq f(\mathbf{X}_{t})+\langle\nabla f(\mathbf{X}_{t}),\mathbf{X}_{t+1}-\mathbf{X}_{t}\rangle+\frac{L}{2}\|\mathbf{X}_{t+1}-\mathbf{X}_{t}\|_{F}^{2}$ $\displaystyle\leq f(\mathbf{X}_{t})-\eta_{t}\langle\nabla f(\mathbf{X}_{t}),\mathbf{O}_{t}\rangle+\frac{L\eta_{t}^{2}}{2}\|\mathbf{O}_{t}\|_{F}^{2}$ $\displaystyle\leq f(\mathbf{X}_{t})-\eta_{t}\langle\mathbf{M}_{t},\mathbf{O}_{t}\rangle-\eta_{t}\langle\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t},\mathbf{O}_{t}\rangle+\frac{L\eta_{t}^{2}}{2}\|\mathbf{O}_{t}\|_{F}^{2}.$ $$

Now we bound the three terms on the right-hand side respectively:

Descent Term: According to the definition of $\mathbf{O}_{t}$ and the norm property $\|\mathbf{M}_{t}\|_{*}\geq\|\mathbf{M}_{t}\|_{F}$, we have:

$$ $-\eta_{t}\langle\mathbf{M}_{t},\mathbf{O}_{t}\rangle=-\eta_{t}\langle\mathbf{M}_{t},\mathbf{U}_{r}\mathbf{V}_{r}^{\top}\rangle=-\eta_{t}\|\mathbf{M}_{t}\|_{*}\leq-\eta_{t}\|\mathbf{M}_{t}\|_{F}.$ $$

Cross Term: This is the key to eliminating the dimension-dependent error. We use Young’s inequality with a parameter ($ab\leq\frac{\alpha}{2}a^{2}+\frac{1}{2\alpha}b^{2}$) and the fact that $\|\mathbf{O}_{t}\|_{F}^{2}\leq\sum_{i=1}^{n}1=n$:

$$ $\displaystyle-\eta_{t}\langle\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t},\mathbf{O}_{t}\rangle$ $\displaystyle\leq\eta_{t}\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}\|\mathbf{O}_{t}\|_{F}$ $\displaystyle\leq\eta_{t}\left(\frac{\alpha}{2}\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}+\frac{1}{2\alpha}\|\mathbf{O}_{t}\|_{F}^{2}\right)$ $\displaystyle\leq\eta_{t}\left(\frac{\alpha}{2}\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}+\frac{n}{2\alpha}\right).$ $$

Quadratic Term:

$$ $\frac{L\eta_{t}^{2}}{2}\|\mathbf{O}_{t}\|_{F}^{2}\leq\frac{L\eta_{t}^{2}n}{2}.$ $$

Substituting these three bounds into the inequality for $f(\mathbf{X}_{t+1})$:

$$ $f(\mathbf{X}_{t+1})\leq f(\mathbf{X}_{t})-\eta_{t}\|\mathbf{M}_{t}\|_{F}+\frac{\eta_{t}\alpha}{2}\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}+\frac{\eta_{t}n}{2\alpha}+\frac{L\eta_{t}^{2}n}{2}.$ $$

This completes the proof.
∎

###### Lemma A.1 .

###### Proof.

### A.2 Lemma A.2

###### Lemma A.2 .

For Algorithm [1](https://arxiv.org/html/2509.15816v3#alg1) option MVR1 ($\gamma=0$), the accumulated error between the momentum term and the true gradient is bounded:

$$ $\displaystyle\mathbb{E}\left[\|\mathbf{M}_{t+1}-\nabla f(\mathbf{X}_{t+1})\|_{F}^{2}\right]\leq\beta_{t+1}\mathbb{E}\left[\|\mathbf{M}_{t}-\nabla f(\mathbf{X}_{t})\|_{F}^{2}\right]+\frac{\beta_{t+1}^{2}}{1-\beta_{t+1}}L^{2}\eta_{t}^{2}n+(1-\beta_{t+1})^{2}\sigma^{2}.$ $$

###### Proof.

First, we have

$$ $\displaystyle\|\mathbf{M}_{t+1}-\nabla f(\mathbf{X}_{t+1})\|_{F}^{2}$ $\displaystyle=\|\beta_{t+1}\mathbf{M}_{t}+(1-\beta_{t+1})\nabla f(\mathbf{X}_{t+1};\xi_{t+1})-\nabla f(\mathbf{X}_{t+1})\|_{F}^{2}$ $\displaystyle=\|\beta_{t+1}(\mathbf{M}_{t}-\nabla f(\mathbf{X}_{t}))+(1-\beta_{t+1})(\nabla f(\mathbf{X}_{t+1};\xi_{t+1})-\nabla f(\mathbf{X}_{t+1}))$ $\displaystyle\quad+\beta_{t+1}(\nabla f(\mathbf{X}_{t})-\nabla f(\mathbf{X}_{t+1}))\|_{F}^{2}$ $\displaystyle=\beta_{t+1}^{2}\|\mathbf{M}_{t}-\nabla f(\mathbf{X}_{t})\|_{F}^{2}+\beta_{t+1}^{2}\|\nabla f(\mathbf{X}_{t})-\nabla f(\mathbf{X}_{t+1})\|_{F}^{2}$ $\displaystyle\quad+(1-\beta_{t+1})^{2}\|\nabla f(\mathbf{X}_{t+1};\xi_{t+1})-\nabla f(\mathbf{X}_{t+1})\|_{F}^{2}$ $\displaystyle\quad+2\beta_{t+1}^{2}\langle\mathbf{M}_{t}-\nabla f(\mathbf{X}_{t}),\nabla f(\mathbf{X}_{t})-\nabla f(\mathbf{X}_{t+1})\rangle_{F}$ $\displaystyle\quad+2\beta_{t+1}(1-\beta_{t+1})\langle\mathbf{M}_{t}-\nabla f(\mathbf{X}_{t}),\nabla f(\mathbf{X}_{t+1};\xi_{t+1})-\nabla f(\mathbf{X}_{t+1})\rangle_{F}$ $\displaystyle\quad+2\beta_{t+1}(1-\beta_{t+1})\langle\nabla f(\mathbf{X}_{t})-\nabla f(\mathbf{X}_{t+1}),\nabla f(\mathbf{X}_{t+1};\xi_{t+1})-\nabla f(\mathbf{X}_{t+1})\rangle_{F}.$ $$

According to Assumption [3.4](https://arxiv.org/html/2509.15816v3#S3.Thmassumption4). Taking the expectation of its squared norm, and using the unbiasedness and independence of the stochastic gradient, we obtain:

$$ $\displaystyle\mathbb{E}[\|\mathbf{M}_{t+1}-\nabla f(\mathbf{X}_{t+1})\|_{F}^{2}]=$ $\displaystyle\beta_{t+1}^{2}\mathbb{E}[\|\mathbf{M}_{t}-\nabla f(\mathbf{X}_{t})\|_{F}^{2}]$ $\displaystyle\quad+(1-\beta_{t+1})^{2}\mathbb{E}[\|\nabla f(\mathbf{X}_{t+1};\xi_{t+1})-\nabla f(\mathbf{X}_{t+1})\|_{F}^{2}]$ $\displaystyle+\beta_{t+1}^{2}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})-\nabla f(\mathbf{X}_{t+1})\|_{F}^{2}]$ $\displaystyle+2\beta_{t+1}^{2}\mathbb{E}[\langle\mathbf{M}_{t}-\nabla f(\mathbf{X}_{t}),\nabla f(\mathbf{X}_{t})-\nabla f(\mathbf{X}_{t+1})\rangle].$ $$

Applying Young’s inequality with a parameter ($ab\leq\frac{\epsilon}{2}a^{2}+\frac{1}{2\epsilon}b^{2}$), we have

$$ $\langle\mathbf{M}_{t}-\nabla f(\mathbf{X}_{t}),\nabla f(\mathbf{X}_{t})-\nabla f(\mathbf{X}_{t+1})\rangle_{F}\leq\frac{\epsilon}{2}\|\mathbf{M}_{t}-\nabla f(\mathbf{X}_{t})\|_{F}^{2}+\frac{1}{2\epsilon}\|\nabla f(\mathbf{X}_{t})-\nabla f(\mathbf{X}_{t+1})\|_{F}^{2}.$ $$

Thus, we have:

$$ $\displaystyle\mathbb{E}\left[\|\mathbf{M}_{t+1}-\nabla f(\mathbf{X}_{t+1})\|_{F}^{2}\right]$ $\displaystyle\leq\beta_{t+1}^{2}(1+\epsilon)\mathbb{E}\left[\|\mathbf{M}_{t}-\nabla f(\mathbf{X}_{t})\|_{F}^{2}\right]$ $\displaystyle\quad+\beta_{t+1}^{2}\left(1+\frac{1}{\epsilon}\right)\mathbb{E}\left[\|\nabla f(\mathbf{X}_{t})-\nabla f(\mathbf{X}_{t+1})\|_{F}^{2}\right]$ $\displaystyle\quad+(1-\beta_{t+1})^{2}\mathbb{E}\left[\|\nabla f(\mathbf{X}_{t+1};\xi_{t+1})-\nabla f(\mathbf{X}_{t+1})\|_{F}^{2}\right].$ $$

According to Assumption [3.2](https://arxiv.org/html/2509.15816v3#S3.Thmassumption2),

$$ $\displaystyle\|\nabla f(\mathbf{X}_{t})-\nabla f(\mathbf{X}_{t+1})\|_{F}^{2}$ $\displaystyle\leq L^{2}\|\mathbf{X}_{t}-\mathbf{X}_{t+1}\|_{F}^{2}$ $\displaystyle=L^{2}\eta_{t}^{2}\|\mathbf{O}_{t}\|_{F}^{2}$ $\displaystyle\leq L^{2}\eta_{t}^{2}n.$ $$

Therefore:

$$ $\displaystyle\mathbb{E}\left[\left\|\mathbf{M}_{t+1}-\nabla f(\mathbf{X}_{t+1})\right\|_{F}^{2}\right]$ $\displaystyle\leq\beta_{t+1}^{2}(1+\epsilon)\mathbb{E}\left[\left\|\mathbf{M}_{t}-\nabla f(\mathbf{X}_{t})\right\|_{F}^{2}\right]$ $\displaystyle+\beta_{t+1}^{2}\left(1+\frac{1}{\epsilon}\right)L^{2}\eta_{t}^{2}n+(1-\beta_{t+1})^{2}\sigma^{2}.$ $$

Then, by letting $\epsilon:=\frac{1-\beta_{t}}{\beta_{t}}$, we have

$$ $\displaystyle\mathbb{E}\left[\|\mathbf{M}_{t+1}-\nabla f(\mathbf{X}_{t+1})\|_{F}^{2}\right]\leq\beta_{t+1}\mathbb{E}\left[\|\mathbf{M}_{t}-\nabla f(\mathbf{X}_{t})\|_{F}^{2}\right]+\frac{\beta_{t+1}^{2}}{1-\beta_{t+1}}L^{2}\eta_{t}^{2}n+(1-\beta_{t+1})^{2}\sigma^{2}.$ (6) $$

∎

###### Lemma A.2 .

###### Proof.

### A.3 Lemma A.3

###### Lemma A.3 .

Suppose that $\{E_{i},A_{i}\}$ are two nonnegative sequences. Assume $E_{t+1}\leq(1-\alpha_{t+1})E_{t}+A_{t+1}$ where $\alpha_{t}=t^{-p}$, $p\in(0,1]$. Then we have:

$$ $\alpha_{t}E_{t}\leq 2(E_{t}-E_{t+1}+A_{t+1}).$ $$

###### Proof.

We derive the following inequalities:

$$ $\displaystyle\alpha_{t}E_{t}-c\left(E_{t}-E_{t+1}+A_{t+1}\right)$ $\displaystyle\stackrel{{\scriptstyle(\bullet)}}{{\leq}}\alpha_{t}E_{t}-c\left(E_{t}+A_{t+1}\right)+c\cdot\left(E_{t}-\alpha_{t+1}E_{t}+A_{t+1}\right)$ $\displaystyle=E_{t}\left(\alpha_{t}-c\alpha_{t+1}\right)$ $\displaystyle=E_{t}\cdot(t+1)^{-p}\cdot\left(\left(\frac{t}{t+1}\right)^{-p}-c\right)$ $\displaystyle\stackrel{{\scriptstyle(\circ)}}{{\leq}}E_{t}\cdot(t+1)^{-p}\cdot(2-c)$ $\displaystyle\stackrel{{\scriptstyle(\star)}}{{\leq}}0,$ $$

where $(\bullet)$ follows from $E_{t+1}\leq(1-\alpha_{t+1})E_{t}+A_{t+1}$; $(\circ)$ is due to $(\frac{t}{t+1})^{-p}\leq 2^{p}\leq 2$; $(\star)$ is due to our choice $c=2$.
∎

###### Lemma A.3 .

###### Proof.

## Appendix B Proofs of Theorem 3.1

###### Proof.

According to Lemma [A.1](https://arxiv.org/html/2509.15816v3#A1.Thmlemma1), we have:

$$ $\displaystyle f(\mathbf{X}_{t+1})$ $\displaystyle\leq f(\mathbf{X}_{t})-\eta_{t}\|\mathbf{M}_{t}\|_{F}+\frac{\eta_{t}\alpha}{2}\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}+\frac{\eta_{t}n}{2\alpha}+\frac{L\eta_{t}^{2}n}{2}$ $\displaystyle\stackrel{{\scriptstyle(\circ)}}{{\leq}}f(\mathbf{X}_{t})-\eta_{t}\|\mathbf{M}_{t}\|_{F}+\frac{\eta_{t}^{2/3}}{2L}\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}+\frac{\eta_{t}^{4/3}Ln}{2}+\frac{L\eta_{t}^{2}n}{2}$ $\displaystyle\stackrel{{\scriptstyle(\star)}}{{\leq}}f(\mathbf{X}_{t})-\eta_{t}\|\mathbf{M}_{t}\|_{F}+\frac{\eta_{t}^{2/3}}{2L}\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}+Ln\eta_{t}^{4/3},$ $$

where $(\circ)$ by setting $\alpha=\frac{1}{\eta_{t}^{1/3}L}$; $(\star)$ follows from $\eta_{t}\leq 1$, we have $L\leq L/\eta_{t}^{1/3}$.

Thus, taking the expectation yields

$$ $\displaystyle\mathbb{E}[f(\mathbf{X}_{t+1})]$ $\displaystyle\leq\mathbb{E}[f(\mathbf{X}_{t})]-\eta_{t}\mathbb{E}[\|\mathbf{M}_{t}\|_{F}]+\frac{\eta_{t}^{2/3}}{2L}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}]+Ln\eta_{t}^{4/3}$ (7) $\displaystyle\stackrel{{\scriptstyle(\circ)}}{{\leq}}\mathbb{E}[f(\mathbf{X}_{t})]-\eta_{t}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})\|_{F}]+\eta_{t}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}]$ $\displaystyle+\frac{\eta_{t}^{2/3}}{2L}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}]+Ln\eta_{t}^{4/3}$ $\displaystyle\stackrel{{\scriptstyle(\star)}}{{\leq}}\mathbb{E}[f(\mathbf{X}_{t})]-\eta_{t}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})\|_{F}]+\frac{1}{2\epsilon}\eta_{t}^{2}+\frac{\epsilon}{2}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}]$ $\displaystyle+\frac{\eta_{t}^{2/3}}{2L}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}]+Ln\eta_{t}^{4/3}$ $\displaystyle\stackrel{{\scriptstyle(\bullet)}}{{=}}\mathbb{E}[f(\mathbf{X}_{t})]-\eta_{t}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})\|_{F}]$ $\displaystyle+\underbrace{\frac{L\eta_{t}^{4/3}}{2}+\left(\frac{\eta_{t}^{2/3}}{2L}+\frac{\eta_{t}^{2/3}}{2L}\right)\mathbb{E}[\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}]+Ln\eta_{t}^{4/3}}_{\Gamma_{t}},$ $$

where $(\circ)$ follows from the reverse triangle inequality $-\|\mathbf{M}_{t}\|_{F}\leq\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}-\|\nabla f(\mathbf{X}_{t})\|_{F}$; $(\star)$ applies Young’s inequality to the term $\eta_{t}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}]$; and $(\bullet)$ collects the residual terms into $\Gamma_{t}$ and sets $\epsilon=\eta_{t}^{2/3}/L$.

Next, we set $\eta_{t}=t^{-3/4}$ and $\beta_{t}=1-t^{-1/2},\alpha_{t}=t^{-1/2}$.

Case 1: $\gamma=0$.

By Lemma [A.2](https://arxiv.org/html/2509.15816v3#A1.Thmlemma2) inequality ([6](https://arxiv.org/html/2509.15816v3#A1.E6)), we have

$$ $\displaystyle\mathbb{E}\left[\|\mathbf{M}_{t+1}-\nabla f(\mathbf{X}_{t+1})\|_{F}^{2}\right]$ $\displaystyle\leq\beta_{t+1}\mathbb{E}\left[\|\mathbf{M}_{t}-\nabla f(\mathbf{X}_{t})\|_{F}^{2}\right]+\frac{\beta_{t+1}^{2}}{1-\beta_{t+1}}L^{2}\eta_{t}^{2}n+(1-\beta_{t+1})^{2}\sigma^{2}$ $\displaystyle\leq\beta_{t+1}\mathbb{E}\left[\|\mathbf{M}_{t}-\nabla f(\mathbf{X}_{t})\|_{F}^{2}\right]+\frac{L^{2}\eta_{t}^{2}n}{1-\beta_{t+1}}+(1-\beta_{t+1})^{2}\sigma^{2}.$ $$

Let $\mathbf{S}_{t+1}=\mathbf{M}_{t+1}-\nabla f(\mathbf{X}_{t+1})$. Thus, setting $\alpha_{t+1}=(t+1)^{-1/2}$, we observe the following relationship for $t\geq 1$:

$$ $\frac{\eta_{t}^{2}}{1-\beta_{t+1}}=\frac{t^{-3/2}}{(t+1)^{-1/2}}=\frac{\sqrt{t+1}}{t^{3/2}}\leq\frac{2\sqrt{2}}{t+1}.$ $$

This allows us to bound the expectation as follows:

$$ $\mathbb{E}\|\mathbf{S}_{t+1}\|_{F}^{2}\leq(1-\alpha_{t+1})\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}+\alpha_{t+1}^{2}(2\sqrt{2}L^{2}n+\sigma^{2}).$ $$

According to Lemma [A.3](https://arxiv.org/html/2509.15816v3#A1.Thmlemma3), by letting $A_{t+1}=\alpha_{t+1}^{2}(2\sqrt{2}L^{2}n+\sigma^{2})$, we have

$$ $\alpha_{t}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}\leq 2(\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}-\mathbb{E}\|\mathbf{S}_{t+1}\|_{F}^{2}+A_{t+1}).$ $$

Furthermore, since

$$ $\displaystyle\mathbb{E}\|\mathbf{S}_{1}\|$ $\displaystyle=\mathbb{E}\|\nabla f(\mathbf{X}_{1})-\mathbf{M}_{1}\|_{F}^{2}=\mathbb{E}\|\nabla f(\mathbf{X}_{1})-(1-\beta_{1})\nabla f(\mathbf{X}_{1};\xi_{1})\|_{F}^{2}$ $\displaystyle=\mathbb{E}\|\nabla f(\mathbf{X}_{1})-\nabla f(\mathbf{X}_{1};\xi_{1})\|_{F}^{2}\leq\sigma^{2}\quad\quad\quad(\text{since }\beta_{1}=0).$ $$

It follows that

$$ $\displaystyle\sum_{t=1}^{T}\alpha_{t}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}$ $\displaystyle\leq 2\sum_{t=1}^{T}(\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}-\mathbb{E}\|\mathbf{S}_{t+1}\|_{F}^{2}+A_{t+1})$ (8) $\displaystyle\leq 2\mathbb{E}\|\mathbf{S}_{1}\|_{F}^{2}+2(2\sqrt{2}L^{2}n+\sigma^{2})\sum_{t=1}^{T}\frac{1}{t+1}$ $\displaystyle\leq 2\sigma^{2}+2(2\sqrt{2}L^{2}n+\sigma^{2})(\ln T+1).$ $$

Thus,

$$ $\displaystyle\Gamma_{t}$ $\displaystyle=\frac{\eta_{t}^{2/3}}{L}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}+(L/2+Ln)\eta_{t}^{4/3}$ $\displaystyle=\frac{\alpha_{t}}{L}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}+(L/2+Ln)\alpha_{t}^{2}.$ $$

Next, we define $A_{1}=4\sqrt{2}Ln+Ln+2L^{-1}\sigma^{2}+L/2$ and $A_{2}=4\sqrt{2}Ln+Ln+4L^{-1}\sigma^{2}+L/2$.

$$ $\displaystyle\sum_{t=1}^{T}\Gamma_{t}$ $\displaystyle=\frac{1}{L}\sum_{t=1}^{T}\alpha_{t}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}+(L/2+Ln)\sum_{t=1}^{T}\alpha_{t}^{2}$ (9) $\displaystyle\stackrel{{\scriptstyle(\circ)}}{{\leq}}\frac{2\sigma^{2}+2(2\sqrt{2}L^{2}n+\sigma^{2})(\ln T+1)}{L}+(L/2+Ln)(\ln T+1)$ $\displaystyle\leq A_{1}\ln T+A_{2},$ $$

where $(\circ)$ due to inequality ([8](https://arxiv.org/html/2509.15816v3#A2.E8)).

Therefore, we have

$$ $\displaystyle\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\|\nabla f(\mathbf{X}_{t})\|_{F}$ $\displaystyle=\frac{1}{T}\sum_{t=1}^{T}t^{3/4}\cdot t^{-3/4}\mathbb{E}\|\nabla f(\mathbf{X}_{t})\|_{F}$ $\displaystyle\leq\frac{1}{T}\sum_{t=1}^{T}t^{3/4}\cdot\eta_{t}\mathbb{E}\|\nabla f(\mathbf{X}_{t})\|_{F}$ $\displaystyle\leq\frac{T^{3/4}}{T}\sum_{t=1}^{T}f(\mathbf{X}_{t})-f(\mathbf{X}_{t+1})+\frac{T^{3/4}}{T}\sum_{t=1}^{T}\Gamma_{t}$ $\displaystyle\leq\frac{f(\mathbf{X}_{1})-f^{*}}{T^{1/4}}+\frac{A_{1}\ln T+A_{2}}{T^{1/4}}$ $\displaystyle=\tilde{\mathcal{O}}(T^{-1/4}).$ $$

Case 2: $\gamma\neq 0$.

We set $\eta_{t}=t^{-3/4}$ and $\beta_{t}=1-(t+1)^{-1/2},\alpha_{t}=t^{-1/2}$, and $\gamma=1-\beta_{t-1}=t^{-1/2}$, we first note that an equivalent form of Algorithm [1](https://arxiv.org/html/2509.15816v3#alg1) Option MVR1 is given by

$$ $\displaystyle\mathbf{C}_{t}$ $\displaystyle=\beta_{t-1}\mathbf{C}_{t-1}+(1-\beta_{t-1})\nabla f(\mathbf{X}_{t};\xi_{t})$ (10) $\displaystyle\mathbf{M}_{t}$ $\displaystyle=\beta_{t}\mathbf{C}_{t}+(1-\beta_{t})\nabla f(\mathbf{X}_{t};\xi_{t}).$ $$

In this case, then we have

$$ $\displaystyle\mathbb{E}\|\mathbf{C}_{t}-\nabla f(\mathbf{X}_{t})\|_{F}^{2}$ $\displaystyle=\mathbb{E}\|\beta_{t-1}\mathbf{C}_{t-1}+(1-\beta_{t-1})\nabla f(\mathbf{X}_{t};\xi_{t})-\nabla f(\mathbf{X}_{t})\|_{F}^{2}$ $\displaystyle=\mathbb{E}\|\beta_{t-1}(\mathbf{C}_{t-1}-\nabla f(\mathbf{X}_{t-1}))+(1-\beta_{t-1})(\nabla f(\mathbf{X}_{t};\xi_{t})-\nabla f(\mathbf{X}_{t}))$ $\displaystyle\quad+\beta_{t-1}(\nabla f(\mathbf{X}_{t-1})-\nabla f(\mathbf{X}_{t}))\|_{F}^{2}$ $\displaystyle\stackrel{{\scriptstyle(\circ)}}{{\leq}}\beta_{t-1}\mathbb{E}\|\mathbf{C}_{t-1}-\nabla f(\mathbf{X}_{t-1})\|_{F}^{2}+\frac{\beta_{t-1}^{2}L\eta_{t-1}^{2}n}{1-\beta_{t-1}}+(1-\beta_{t-1})^{2}\sigma^{2}$ $\displaystyle=(1-\alpha_{t})\mathbb{E}\|\mathbf{C}_{t-1}-\nabla f(\mathbf{X}_{t-1})\|_{F}^{2}+\alpha_{t}^{2}(2\sqrt{2}Ln+\sigma^{2}),$ $$

where ($\circ$) follows from Lemma [A.2](https://arxiv.org/html/2509.15816v3#A1.Thmlemma2).

Then, we define $\mathbf{S}_{t}^{\prime}=\|\mathbf{M}_{t}-\nabla f(\mathbf{X}_{t})\|_{F}^{2},\mathbf{S}_{t}=\|\mathbf{C}_{t}-\nabla f(\mathbf{X}_{t})\|_{F}^{2},A_{t+1}=\alpha_{t+1}^{2}(2\sqrt{2}Ln+\sigma^{2})$. Using the conclusion of Lemma [A.3](https://arxiv.org/html/2509.15816v3#A1.Thmlemma3), we have

$$ $\mathbb{E}\|\mathbf{S}_{t+1}\|_{F}^{2}\leq(1-\alpha_{t+1})\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}+A_{t+1}.$ $$

Furthermore, since

$$ $\displaystyle\mathbb{E}\|\mathbf{S}_{1}\|$ $\displaystyle=\mathbb{E}\|\nabla f(\mathbf{X}_{1})-\mathbf{C}_{1}\|_{F}^{2}=\mathbb{E}\|\nabla f(\mathbf{X}_{1})-(1-\beta_{0})\nabla f(\mathbf{X}_{1};\xi_{1})\|_{F}^{2}$ $\displaystyle=\mathbb{E}\|\nabla f(\mathbf{X}_{1})-\nabla f(\mathbf{X}_{1};\xi_{1})\|_{F}^{2}\leq\sigma^{2}\quad\quad\quad(\text{since }\beta_{0}=0).$ $$

Thus, as in inequality ([8](https://arxiv.org/html/2509.15816v3#A2.E8)), we have

$$ $\displaystyle\sum_{t=1}^{T}\alpha_{t}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}$ $\displaystyle\leq 2\sum_{t=1}^{T}(\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}-\mathbb{E}\|\mathbf{S}_{t+1}\|_{F}^{2}+A_{t+1})$ $\displaystyle\leq 2\mathbb{E}\|\mathbf{S}_{1}\|_{F}^{2}+2(2\sqrt{2}L^{2}n+\sigma^{2})\sum_{t=1}^{T}\frac{1}{t+1}$ $\displaystyle\leq 2\sigma^{2}+2(2\sqrt{2}L^{2}n+\sigma^{2})(\ln T+1).$ $$

Then,

$$ $\displaystyle\|\mathbf{M}_{t}-\nabla f(\mathbf{X}_{t})\|_{F}$ $\displaystyle=\|\beta_{t}\mathbf{C}_{t}+(1-\beta_{t})\nabla f(\mathbf{X}_{t};\xi_{t})-\nabla f(\mathbf{X}_{t})\|_{F}$ $\displaystyle=\|\beta_{t}(\mathbf{C}_{t}-\nabla f(\mathbf{X}_{t}))+(1-\beta_{t})(\nabla f(\mathbf{X}_{t};\xi_{t})-\nabla f(\mathbf{X}_{t}))\|_{F}.$ $$

From this, we can bound the squared norm using the inequality $\|\mathbf{A}+\mathbf{B}\|_{F}^{2}\leq 2\|\mathbf{A}\|_{F}^{2}+2\|\mathbf{B}\|_{F}^{2}$:

$$ $\displaystyle\|\mathbf{M}_{t}-\nabla f(\mathbf{X}_{t})\|_{F}^{2}\leq 2\beta_{t}^{2}\|\mathbf{C}_{t}-\nabla f(\mathbf{X}_{t})\|_{F}^{2}+2(1-\beta_{t})^{2}\|\nabla f(\mathbf{X}_{t};\xi_{t})-\nabla f(\mathbf{X}_{t})\|_{F}^{2}.$ $$

Thus, we have

$$ $\displaystyle\mathbb{E}\|\mathbf{M}_{t}-\nabla f(\mathbf{X}_{t})\|_{F}^{2}$ $\displaystyle\leq 2\beta_{t}^{2}\mathbb{E}\|\mathbf{C}_{t}-\nabla f(\mathbf{X}_{t})\|_{F}^{2}+2(1-\beta_{t})^{2}\mathbb{E}\|\nabla f(\mathbf{X}_{t};\xi_{t})-\nabla f(\mathbf{X}_{t})\|_{F}^{2}$ (11) $\displaystyle\leq 2\beta_{t}^{2}\mathbb{E}\|\mathbf{C}_{t}-\nabla f(\mathbf{X}_{t})\|_{F}^{2}+2(1-\beta_{t})^{2}\sigma^{2}$ $\displaystyle\leq 2\mathbb{E}\|\mathbf{C}_{t}-\nabla f(\mathbf{X}_{t})\|_{F}^{2}+2(1-\beta_{t})^{2}\sigma^{2}.$ $$

Thus,

$$ $\displaystyle\sum_{t=1}^{T}\alpha_{t}\mathbb{E}\|\mathbf{S}^{\prime}_{t}\|_{F}^{2}$ $\displaystyle\leq 2\sum_{t=1}^{T}\alpha_{t}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}+2\sigma^{2}\sum_{t=1}^{T}\frac{1}{t(t+1)}$ (12) $\displaystyle\stackrel{{\scriptstyle(\circ)}}{{\leq}}4\sigma^{2}+4(2\sqrt{2}L^{2}n+\sigma^{2})(\ln T+1)+2\sigma^{2}$ $\displaystyle=6\sigma^{2}+4(2\sqrt{2}L^{2}n+\sigma^{2})(\ln T+1),$ $$

where $(\circ)$ is due to $\sum_{t=1}^{T}\frac{1}{t(t+1)}=\sum_{t=1}^{T}(\frac{1}{t}-\frac{1}{t+1})=1-\frac{1}{T+1}\leq 1$.

Therefore,

$$ $\displaystyle\Gamma_{t}$ $\displaystyle=\frac{\eta_{t}^{2/3}}{L}\mathbb{E}\|\mathbf{S}^{\prime}_{t}\|_{F}^{2}+(L/2+Ln)\eta_{t}^{4/3}$ $\displaystyle=\frac{\alpha_{t}}{L}\mathbb{E}\|\mathbf{S}^{\prime}_{t}\|_{F}^{2}+(L/2+Ln)\alpha_{t}^{2}.$ $$

Then, we define $A_{1}=8\sqrt{2}Ln+Ln+4L^{-1}\sigma^{2}+L/2$ and $A_{2}=8\sqrt{2}Ln+Ln+10L^{-1}\sigma^{2}+L/2$. We have

$$ $\displaystyle\sum_{t=1}^{T}\Gamma_{t}$ $\displaystyle=\frac{1}{L}\sum_{t=1}^{T}\alpha_{t}\mathbb{E}\|\mathbf{S}^{\prime}_{t}\|_{F}^{2}+(L/2+Ln)\sum_{t=1}^{T}\alpha_{t}^{2}$ (13) $\displaystyle\stackrel{{\scriptstyle(\circ)}}{{\leq}}\frac{6\sigma^{2}+4(2\sqrt{2}L^{2}n+\sigma^{2})(\ln T+1)}{L}+(L/2+Ln)(\ln T+1)$ $\displaystyle\leq A_{1}\ln T+A_{2},$ $$

where ($\circ$) is due to inequality ([12](https://arxiv.org/html/2509.15816v3#A2.E12)).

Then, we have

$$ $\displaystyle\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\|\nabla f(\mathbf{X}_{t})\|_{F}$ $\displaystyle=\frac{1}{T}\sum_{t=1}^{T}t^{3/4}\cdot t^{-3/4}\mathbb{E}\|\nabla f(\mathbf{X}_{t})\|_{F}$ $\displaystyle\leq\frac{1}{T}\sum_{t=1}^{T}t^{3/4}\cdot\eta_{t}\mathbb{E}\|\nabla f(\mathbf{X}_{t})\|_{F}$ $\displaystyle\leq\frac{T^{3/4}}{T}\sum_{t=1}^{T}f(\mathbf{X}_{t})-f(\mathbf{X}_{t+1})+\frac{T^{3/4}}{T}\sum_{t=1}^{T}\Gamma_{t}$ $\displaystyle\leq\frac{f(\mathbf{X}_{1})-f^{*}}{T^{1/4}}+\frac{A_{1}\ln T+A_{2}}{T^{1/4}}$ $\displaystyle=\tilde{\mathcal{O}}(T^{-1/4}).$ $$

This completes the proof.
∎

###### Proof.

## Appendix C Lemmas for Theorem 3.2

### C.1 Lemma C.1

###### Lemma C.1 .

Let $\{A_{t}\}_{t\geq 1}$ and $\{B_{t}\}_{t\geq 1}$ be non-negative sequences satisfying the relation $A_{t+1}\leq(1-\varepsilon_{t+1})A_{t}+B_{t+1}$ for all $t\geq 1$. If we define the sequence $\varepsilon_{t}=t^{-p}$ for some constant $p\in(0,1]$, then for all $t\geq 1$, the following inequality holds:

$$ $\sqrt{\varepsilon_{t}}A_{t}\leq 4\left(\frac{A_{t}}{\sqrt{\varepsilon_{t}}}-\frac{A_{t+1}}{\sqrt{\varepsilon_{t+1}}}+\frac{B_{t+1}}{\sqrt{\varepsilon_{t+1}}}\right).$ $$

###### Proof.

First, we define the function $F(t,q)=\frac{1}{4}t^{-q}-t^{q}+(t+1)^{q}-(t+1)^{-q}$ for $t\geq 1$ and $q\in(0,1/2]$. To analyze its properties, let us define $g(t)=F(t,q)$ for a fixed $q\in(0,1/2]$ and $f(q)=F(t,q)$ for a fixed $t\geq 1$.

Case 1: For $t=1$.
We have $g(1)=\frac{1}{4}-1+2^{q}-\frac{1}{2^{q}}$. Since $q\in(0,1/2]$, this expression is bounded above by its value at $q=1/2$, yielding $g(1)\leq\frac{1}{4}-1+2^{1/2}-\frac{1}{2^{1/2}}<0$. Thus, for $t=1$, the inequality $g(t)\leq 0$ holds.

Case 2: For $t>1$.
For any given $t>1$ and $q\in(0,1/2]$, we examine the derivative of $f(q)$:

$$ $\displaystyle f^{\prime}(q)$ $\displaystyle=\left(\ln(t+1)(t+1)^{-q}-\frac{1}{4}\ln(t)t^{-q}\right)+(\ln(t+1)(t+1)^{q}-\ln(t)t^{q})$ (14) $\displaystyle\stackrel{{\scriptstyle(\circ)}}{{\geq}}\left(\ln(t+1)(t+1)^{-q}-\frac{1}{4}\ln(t)t^{-q}\right)$ $\displaystyle=\ln(t)t^{-q}\cdot\left\{\frac{\ln(t+1)}{\ln(t)}\cdot\left(\frac{t}{t+1}\right)^{q}-\frac{1}{4}\right\}$ $\displaystyle\stackrel{{\scriptstyle(\star)}}{{\geq}}\ln(t)t^{-q}\cdot\left\{1\cdot\sqrt{\frac{t}{t+1}}-\frac{1}{4}\right\}$ $\displaystyle\stackrel{{\scriptstyle(\bullet)}}{{\geq}}\ln(t)t^{-q}\cdot\left\{\sqrt{\frac{1}{2}}-\frac{1}{4}\right\}\geq 0,$ $$

where $(\circ)$ holds because $\ln(t+1)>\ln(t)$ and $(t+1)^{q}>t^{q}$ for all $t\geq 1$ and $q\in(0,1/2]$. The inequality $(\star)$ holds because $\frac{\ln(t+1)}{\ln(t)}>1$ for $t>1$, and the function $(\frac{t}{t+1})^{q}$ is decreasing in $q$, thus its minimum on $(0,1/2]$ is achieved at $q=1/2$. The final inequality $(\bullet)$ holds because $t\geq 1$ implies $\sqrt{\frac{t}{t+1}}\geq\sqrt{\frac{1}{2}}$, and $\sqrt{\frac{1}{2}}>1/4$.

Inequality ([14](https://arxiv.org/html/2509.15816v3#A3.E14)) implies that $f(q)$ is monotonically increasing with respect to $q$ on the interval $(0,1/2]$.

Next, we consider the boundary condition at $q=1/2$. Let $h(t):=f(\frac{1}{2})=\frac{1}{4}t^{-1/2}-t^{1/2}+(t+1)^{1/2}-(t+1)^{-1/2}$. It can be verified that $h(t)\leq 0$ for all $t>1$. Consequently, we have $F(t,\frac{1}{2})=f(\frac{1}{2})\leq 0$ for all $t>1$.

Finally, for all $t>1$ and $q\in(0,1/2]$, we have

$$ $F(t,q)\leq F(t,1/2)\leq 0,$ (15) $$

where the first inequality holds because $f(q)$ is monotonically increasing in $q$ on $(0,1/2]$ for any fixed $t>1$.

With this result, we can proceed as follows:

$$ $\displaystyle\sqrt{\varepsilon_{t}}A_{t}-4\left(\frac{A_{t}}{\sqrt{\varepsilon_{t}}}-\frac{A_{t+1}}{\sqrt{\varepsilon_{t+1}}}+\frac{B_{t+1}}{\sqrt{\varepsilon_{t+1}}}\right)$ $\displaystyle\stackrel{{\scriptstyle(\circ)}}{{\leq}}\sqrt{\varepsilon_{t}}A_{t}-4\left(\frac{A_{t}}{\sqrt{\varepsilon_{t}}}+\frac{B_{t+1}}{\sqrt{\varepsilon_{t+1}}}\right)+\frac{4}{\sqrt{\varepsilon_{t+1}}}\cdot\left(A_{t}-\varepsilon_{t+1}A_{t}+B_{t+1}\right)$ $\displaystyle=\sqrt{\varepsilon_{t}}A_{t}-4\frac{A_{t}}{\sqrt{\varepsilon_{t}}}+\frac{4}{\sqrt{\varepsilon_{t+1}}}(1-\varepsilon_{t+1})A_{t}$ $\displaystyle=A_{t}\cdot\left(\sqrt{\varepsilon_{t}}-\frac{4}{\sqrt{\varepsilon_{t}}}+\frac{4}{\sqrt{\varepsilon_{t+1}}}-4\sqrt{\varepsilon_{t+1}}\right)$ $\displaystyle\stackrel{{\scriptstyle(\star)}}{{=}}4A_{t}\cdot\left(\frac{1}{4}t^{-q}-t^{q}+(t+1)^{q}-(t+1)^{-q}\right)$ $\displaystyle\stackrel{{\scriptstyle(\bullet)}}{{\leq}}0.$ $$

Here, $(\circ)$ follows from the assumption $A_{t+1}\leq(1-\varepsilon_{t+1})A_{t}+B_{t+1}$. The equality $(\star)$ is obtained by substituting $\varepsilon_{t}=t^{-p}$ and setting $q=\frac{p}{2}$ (note that $p\in(0,1]$ implies $q\in(0,1/2]$). Finally, $(\bullet)$ is a direct consequence of our result in inequality ([15](https://arxiv.org/html/2509.15816v3#A3.E15)).
∎

###### Lemma C.1 .

###### Proof.

### C.2 Lemma C.2

###### Lemma C.2 .

For Algorithm [1](https://arxiv.org/html/2509.15816v3#alg1) option MVR2, let $\Delta_{t}=\nabla f(\mathbf{X}_{t+1};\xi_{t+1})-\nabla f(\mathbf{X}_{t};\xi_{t+1})$, $\delta_{t}=\nabla f(\mathbf{X}_{t})-\nabla f(\mathbf{X}_{t+1})$, $\mathbf{S}_{t}=\mathbf{M}_{t}-\nabla f(\mathbf{X}_{t})$, and $\mathbf{R}_{t+1}=\nabla f(\mathbf{X}_{t+1};\xi_{t+1})-\nabla f(\mathbf{X}_{t+1})$. Then we have the following inequality:

$$ $\mathbb{E}\|\mathbf{S}_{t+1}\|_{F}^{2}\leq\beta_{t+1}^{2}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}+2\beta_{t+1}^{2}L^{2}\mathbb{E}\|\mathbf{X}_{t+1}-\mathbf{X}_{t}\|_{F}^{2}+2(1-\beta_{t+1})^{2}\sigma^{2}+P_{t+1},$ $$

where

$$ $\displaystyle A_{t+1}$ $\displaystyle=\frac{B_{t+1}+\beta_{t+1}(\mathbb{E}\|\Delta_{t}\|_{F}^{2}-\|\mathbb{E}\Delta_{t}\|_{F}^{2})}{\mathbb{E}\|\Delta_{t}\|_{F}^{2}}$ $\displaystyle B_{t+1}$ $\displaystyle=(1-\beta_{t+1})\mathbb{E}\langle\Delta_{t},\mathbf{R}_{t+1}\rangle+\beta_{t+1}\mathbb{E}\langle\Delta_{t},\mathbf{S}_{t}\rangle$ $\displaystyle P_{t+1}$ $\displaystyle=\mathbb{E}\|\Delta_{t}\|_{F}^{2}(\beta_{t+1}(1-\gamma_{t+1})-A_{t+1})^{2}-\mathbb{E}\|\Delta_{t}\|_{F}^{2}A_{t+1}^{2}.$ $$

If we choose $\gamma_{t+1}=1-\frac{A_{t+1}}{\beta_{t+1}}$ or $\gamma_{t+1}=1$, then

$$ $\mathbb{E}\|\mathbf{S}_{t+1}\|_{F}^{2}\leq\beta_{t+1}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}+2\beta_{t+1}^{2}L^{2}\mathbb{E}\|\mathbf{X}_{t+1}-\mathbf{X}_{t}\|_{F}^{2}+2(1-\beta_{t+1})^{2}\sigma^{2}.$ $$

###### Proof.

First, we have

$$ $\displaystyle\mathbf{M}_{t+1}$ $\displaystyle=\beta_{t+1}\mathbf{M}_{t}+(1-\beta_{t+1})\nabla f(\mathbf{X}_{t+1};\xi_{t+1})+\gamma_{t+1}\cdot\beta_{t+1}(\nabla f(\mathbf{X}_{t+1};\xi_{t+1})-\nabla f(\mathbf{X}_{t};\xi_{t+1}))$ $\displaystyle=(1-\beta_{t+1})\nabla f(\mathbf{X}_{t+1};\xi_{t+1})+\beta_{t+1}(\mathbf{M}_{t}+\gamma_{t+1}(\nabla f(\mathbf{X}_{t+1};\xi_{t+1})-\nabla f(\mathbf{X}_{t};\xi_{t+1})))$ $\displaystyle=(1-\beta_{t+1})\nabla f(\mathbf{X}_{t+1};\xi_{t+1})+\beta_{t+1}(\mathbf{M}_{t}+\gamma_{t+1}\Delta_{t}).$ $$

Hence,

$$ $\displaystyle\mathbf{M}_{t+1}-\nabla f(\mathbf{X}_{t+1})$ $\displaystyle=(1-\beta_{t+1})\left(\nabla f(\mathbf{X}_{t+1};\xi_{t+1})-\nabla f(\mathbf{X}_{t+1})\right)$ $\displaystyle\quad+\beta_{t+1}\left(\mathbf{M}_{t}-\nabla f(\mathbf{X}_{t})+\delta_{t}+\gamma_{t+1}\Delta_{t}\right)$ $\displaystyle=(1-\beta_{t+1})\mathbf{R}_{t+1}+\beta_{t+1}(\mathbf{S}_{t}+\delta_{t}+\gamma_{t+1}\Delta_{t}).$ $$

Note that, according to Assumption [3.4](https://arxiv.org/html/2509.15816v3#S3.Thmassumption4):

$$ $\delta_{t}=\nabla f(\mathbf{X}_{t})-\nabla f(\mathbf{X}_{t+1})=-\mathbb{E}_{t+1}[\Delta_{t}].$ $$

Therefore,

$$ $\mathbf{S}_{t+1}=(1-\beta_{t+1})\mathbf{R}_{t+1}+\beta_{t+1}\mathbf{S}_{t}+\beta_{t+1}(\gamma_{t+1}\Delta_{t}-\mathbb{E}\Delta_{t}).$ $$

Thus,

$$ $\displaystyle\mathbb{E}\|\mathbf{S}_{t+1}\|_{F}^{2}$ $\displaystyle=\mathbb{E}\|(1-\beta_{t+1})\mathbf{R}_{t+1}+\beta_{t+1}\mathbf{S}_{t}+\beta_{t+1}(\gamma_{t+1}\Delta_{t}-\mathbb{E}\Delta_{t})\|_{F}^{2}$ $\displaystyle=\mathbb{E}\|(1-\beta_{t+1})\mathbf{R}_{t+1}+\beta_{t+1}\mathbf{S}_{t}+\beta_{t+1}((\gamma_{t+1}-1)\Delta_{t}+\Delta_{t}-\mathbb{E}\Delta_{t})\|_{F}^{2}$ $\displaystyle=\underbrace{\mathbb{E}\|(1-\beta_{t+1})\mathbf{R}_{t+1}+\beta_{t+1}\mathbf{S}_{t}+\beta_{t+1}(\Delta_{t}-\mathbb{E}\Delta_{t})\|_{F}^{2}}_{\text{Term A.1}}$ $\displaystyle\quad+\underbrace{\beta_{t+1}^{2}(\gamma_{t+1}-1)^{2}\mathbb{E}\|\Delta_{t}\|_{F}^{2}}_{\text{Term A.2}}$ $\displaystyle\quad+\underbrace{2\beta_{t+1}(\gamma_{t+1}-1)\mathbb{E}\langle\Delta_{t},(1-\beta_{t+1})\mathbf{R}_{t+1}+\beta_{t+1}\mathbf{S}_{t}+\beta_{t+1}(\Delta_{t}-\mathbb{E}\Delta_{t})\rangle}_{\text{Term A.3}}.$ $$

First, let’s consider Term A.1:

$$ A.1 $\displaystyle=\mathbb{E}\|(1-\beta_{t+1})\mathbf{R}_{t+1}+\beta_{t+1}\mathbf{S}_{t}+\beta_{t+1}(\Delta_{t}-\mathbb{E}\Delta_{t})\|_{F}^{2}$ $\displaystyle\stackrel{{\scriptstyle(\circ)}}{{=}}\mathbb{E}\|(1-\beta_{t+1})\mathbf{R}_{t+1}+\beta_{t+1}(\Delta_{t}-\mathbb{E}\Delta_{t})\|_{F}^{2}+\beta_{t+1}^{2}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}$ $\displaystyle\stackrel{{\scriptstyle(\star)}}{{\leq}}2(1-\beta_{t+1})^{2}\mathbb{E}\|\mathbf{R}_{t+1}\|_{F}^{2}+2\beta_{t+1}^{2}\mathbb{E}\|\Delta_{t}-\mathbb{E}\Delta_{t}\|_{F}^{2}+\beta_{t+1}^{2}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2},$ $$

where $(\circ)$ holds because Assumption [3.4](https://arxiv.org/html/2509.15816v3#S3.Thmassumption4) implies $\mathbb{E}[\mathbf{R}_{t+1}]=0$ and $\mathbb{E}[\Delta_{t}-\mathbb{E}\Delta_{t}]=0$, making the cross-terms with $\mathbf{S}_{t}$ zero; and $(\star)$ follows from the inequality $\|\mathbf{A}+\mathbf{B}\|_{F}^{2}\leq 2\|\mathbf{A}\|_{F}^{2}+2\|\mathbf{B}\|_{F}^{2}$ for any $\mathbf{A},\mathbf{B}\in\mathbb{R}^{m\times n}$.

Next, considering the sum of Terms A.2 and A.3, let us define

$$ $\displaystyle A_{t+1}$ $\displaystyle=\frac{B_{t+1}+\beta_{t+1}(\mathbb{E}\|\Delta_{t}\|_{F}^{2}-\|\mathbb{E}\Delta_{t}\|_{F}^{2})}{\mathbb{E}\|\Delta_{t}\|_{F}^{2}}$ $\displaystyle B_{t+1}$ $\displaystyle=(1-\beta_{t+1})\mathbb{E}\langle\Delta_{t},\mathbf{R}_{t+1}\rangle+\beta_{t+1}\mathbb{E}\langle\Delta_{t},\mathbf{S}_{t}\rangle.$ $$

Then we have

$$ $\displaystyle\text{A.2}+\text{A.3}$ $\displaystyle=\beta_{t+1}^{2}(\gamma_{t+1}-1)^{2}\mathbb{E}\|\Delta_{t}\|_{F}^{2}$ $\displaystyle+2(\gamma_{t+1}-1)\beta_{t+1}\mathbb{E}\langle\Delta_{t},(1-\beta_{t+1})\mathbf{R}_{t+1}+\beta_{t+1}\mathbf{S}_{t}+\beta_{t+1}(\Delta_{t}-\mathbb{E}\Delta_{t})\rangle$ $\displaystyle=\beta_{t+1}^{2}\mathbb{E}\|\Delta_{t}\|_{F}^{2}\left((\gamma_{t+1}-1)^{2}+2(\gamma_{t+1}-1)\frac{A_{t+1}}{\beta_{t+1}}\right)$ $\displaystyle=\beta_{t+1}^{2}\mathbb{E}\|\Delta_{t}\|_{F}^{2}\left(\gamma_{t+1}-1+\frac{A_{t+1}}{\beta_{t+1}}\right)^{2}-\beta_{t+1}^{2}\mathbb{E}\|\Delta_{t}\|_{F}^{2}\left(\frac{A_{t+1}}{\beta_{t+1}}\right)^{2}$ $\displaystyle=\mathbb{E}\|\Delta_{t}\|_{F}^{2}(\beta_{t+1}(1-\gamma_{t+1})-A_{t+1})^{2}-A_{t+1}^{2}\mathbb{E}\|\Delta_{t}\|_{F}^{2}$ $\displaystyle=P_{t+1}.$ $$

Therefore, combining the bounds:

$$ $\displaystyle\mathbb{E}\|\mathbf{S}_{t+1}\|_{F}^{2}$ $\displaystyle\leq 2(1-\beta_{t+1})^{2}\mathbb{E}\|\mathbf{R}_{t+1}\|_{F}^{2}+2\beta_{t+1}^{2}\mathbb{E}\|\Delta_{t}-\mathbb{E}\Delta_{t}\|_{F}^{2}+\beta_{t+1}^{2}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}+P_{t+1}$ $\displaystyle\stackrel{{\scriptstyle(\circ)}}{{\leq}}2(1-\beta_{t+1})^{2}\sigma^{2}+2\beta_{t+1}^{2}\mathbb{E}\|\Delta_{t}\|_{F}^{2}+\beta_{t+1}^{2}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}+P_{t+1}$ $\displaystyle\stackrel{{\scriptstyle(\star)}}{{\leq}}2(1-\beta_{t+1})^{2}\sigma^{2}+2\beta_{t+1}^{2}L^{2}\mathbb{E}\|\mathbf{X}_{t+1}-\mathbf{X}_{t}\|_{F}^{2}+\beta_{t+1}^{2}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}+P_{t+1},$ $$

where $(\circ)$ uses Assumption [3.4](https://arxiv.org/html/2509.15816v3#S3.Thmassumption4) (bounded variance, i.e., $\mathbb{E}\|\mathbf{R}_{t+1}\|_{F}^{2}\leq\sigma^{2}$) and the property that $\mathbb{E}\|\Delta_{t}-\mathbb{E}\Delta_{t}\|_{F}^{2}\leq\mathbb{E}\|\Delta_{t}\|_{F}^{2}$; and $(\star)$ follows from Assumption [3.3](https://arxiv.org/html/2509.15816v3#S3.Thmassumption3), which implies $\|\Delta_{t}\|_{F}^{2}\leq L^{2}\|\mathbf{X}_{t+1}-\mathbf{X}_{t}\|_{F}^{2}$.

Next, we set

$$ $\gamma_{t+1}=1-\frac{A_{t+1}}{\beta_{t+1}}\text{ or }\gamma_{t+1}=\gamma=1.$ (16) $$

Then,

$$ $P_{t+1}=-\mathbb{E}\|\Delta_{t}\|_{F}^{2}A_{t+1}^{2}\leq 0\text{ or }P_{t+1}=0.$ $$

This leads to the final result:

$$ $\displaystyle\mathbb{E}\|\mathbf{S}_{t+1}\|_{F}^{2}$ $\displaystyle\leq\beta_{t+1}^{2}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}+2\beta_{t+1}^{2}L^{2}\mathbb{E}\|\mathbf{X}_{t+1}-\mathbf{X}_{t}\|_{F}^{2}+2(1-\beta_{t+1})^{2}\sigma^{2}$ $\displaystyle\stackrel{{\scriptstyle(\circ)}}{{\leq}}\beta_{t+1}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}+2\beta_{t+1}^{2}L^{2}\mathbb{E}\|\mathbf{X}_{t+1}-\mathbf{X}_{t}\|_{F}^{2}+2(1-\beta_{t+1})^{2}\sigma^{2},$ $$

where $(\circ)$ follows from $\beta_{t+1}^{2}\leq\beta_{t+1}\leq 1$.
∎

###### Lemma C.2 .

###### Proof.

### C.3 Lemma C.3

###### Lemma C.3 .

Let $\{\mathbf{S}_{t}\}_{t\geq 1}$ be a sequence of matrices satisfying the recursive inequality $\mathbb{E}\|\mathbf{S}_{t+1}\|_{F}^{2}\leq(1-\eta_{t+1})\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}+2\eta_{t+1}^{2}(2L^{2}n+\sigma^{2})$ for some constants $L,n,\sigma^{2}>0$. If we set the step size $\eta_{t}=t^{-2/3}$ and $\gamma=1$, then we have the following upper bound on the time-averaged expectation:

$$ $\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}\leq\frac{4\sigma^{2}+8(2L^{2}n+\sigma^{2})(1+\ln T)}{T^{2/3}}.$ $$

###### Proof.

The proof begins with the recursive inequality derived from a preceding Lemma [C.2](https://arxiv.org/html/2509.15816v3#A3.Thmlemma2):

$$ $\mathbb{E}\|\mathbf{S}_{t+1}\|_{F}^{2}\leq(1-\eta_{t+1})\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}+2\eta_{t}^{2}L^{2}n+2\eta_{t+1}^{2}\sigma^{2},$ $$

where $\mathbf{S}_{t+1}=\nabla f(\mathbf{X}_{t+1})-\mathbf{M}_{t+1}$.
We have noticed the following facts

$$ $\frac{1}{t^{2/3}}\leq\frac{2}{(t+1)^{2/3}}.$ $$

Therefore , we have

$$ $\mathbb{E}\|\mathbf{S}_{t+1}\|_{F}^{2}\leq(1-\eta_{t+1})\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}+2\eta_{t+1}^{2}(2L^{2}n+\sigma^{2}).$ $$

Let $A_{t}=\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}$ and $B_{t+1}=2\eta_{t+1}^{2}(2L^{2}n+\sigma^{2})$. The inequality can be written as $A_{t+1}\leq(1-\eta_{t+1})A_{t}+B_{t+1}$. This structure allows us to apply a standard result Lemma [C.1](https://arxiv.org/html/2509.15816v3#A3.Thmlemma1) which yields:

$$ $\sqrt{\eta_{t}}A_{t}\leq 4\left(\frac{A_{t}}{\sqrt{\eta_{t}}}-\frac{A_{t+1}}{\sqrt{\eta_{t+1}}}+\frac{B_{t+1}}{\sqrt{\eta_{t+1}}}\right).$ $$

Then, we define $P_{t}=\frac{4A_{t}}{\sqrt{\eta_{t}}}=\frac{4\mathbb{E}\|\mathbf{S}_{t}\|^{2}_{F}}{\sqrt{\eta_{t}}}$. Substituting $P_{t}$ and the definition of $B_{t+1}$ into the inequality gives:

$$ $\sqrt{\eta_{t}}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}\leq P_{t}-P_{t+1}+\frac{4\cdot 2\eta_{t+1}^{2}(2L^{2}n+\sigma^{2})}{\sqrt{\eta_{t+1}}}=P_{t}-P_{t+1}+8\eta_{t+1}^{3/2}(2L^{2}n+\sigma^{2}).$ $$

Now, we sum this inequality from $t=1$ to $T$:

$$ $\displaystyle\sum_{t=1}^{T}\sqrt{\eta_{t}}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}$ $\displaystyle\leq\sum_{t=1}^{T}\left(P_{t}-P_{t+1}+8\eta_{t+1}^{3/2}(2L^{2}n+\sigma^{2})\right)$ $\displaystyle=(P_{1}-P_{T+1})+8(2L^{2}n+\sigma^{2})\sum_{t=1}^{T}\eta_{t+1}^{3/2}.$ $$

Since $P_{T+1}\geq 0$, we can drop this term to simplify the bound. By setting the step size $\eta_{t}=t^{-2/3}$, we have $\eta_{t+1}^{3/2}=((t+1)^{-2/3})^{3/2}=(t+1)^{-1}$. The summation becomes:

$$ $\displaystyle\sum_{t=1}^{T}\sqrt{\eta_{t}}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}$ $\displaystyle\leq P_{1}+8(2L^{2}n+\sigma^{2})\sum_{t=1}^{T}\frac{1}{t+1}$ (17) $\displaystyle\leq P_{1}+8(2L^{2}n+\sigma^{2})\sum_{t=1}^{T}\frac{1}{t}$ $\displaystyle\stackrel{{\scriptstyle(\circ)}}{{\leq}}P_{1}+8(2L^{2}n+\sigma^{2})(1+\ln T),$ $$

where $(\circ)$ follows from the harmonic series, $\sum_{t=1}^{T}\frac{1}{t}\leq 1+\ln T$.

Finally, we establish the bound on the time-averaged expectation. With our choice of $\eta_{t}=t^{-2/3}$, we have $\sqrt{\eta_{t}}=t^{-1/3}$. Therefore:

$$ $\displaystyle\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}$ $\displaystyle=\frac{1}{T}\sum_{t=1}^{T}t^{1/3}\cdot t^{-1/3}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}$ $\displaystyle=\frac{1}{T}\sum_{t=1}^{T}t^{1/3}\sqrt{\eta_{t}}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}$ $\displaystyle\leq\frac{T^{1/3}}{T}\sum_{t=1}^{T}\sqrt{\eta_{t}}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}$ $\displaystyle\text{(since }t^{1/3}\leq T^{1/3}\text{ for }t\leq T\text{)}$ $\displaystyle\stackrel{{\scriptstyle(\circ)}}{{\leq}}\frac{T^{1/3}}{T}\left(P_{1}+8(2L^{2}n+\sigma^{2})(1+\ln T)\right)$ $\displaystyle=\frac{P_{1}+8(2L^{2}n+\sigma^{2})(1+\ln T)}{T^{2/3}},$ $$

where $(\circ)$ follows from the inequality ([17](https://arxiv.org/html/2509.15816v3#A3.E17)).

Substituting $P_{1}=4\mathbb{E}\|\mathbf{S}_{1}\|_{F}^{2}/\sqrt{\eta_{1}}=4\mathbb{E}\|\mathbf{S}_{1}\|_{F}^{2}$ completes the main proof. To establish the final bound, we now analyze the initial term $P_{1}$. Given the definition $\mathbf{S}_{1}=\nabla f(\mathbf{X}_{1})-\mathbf{M}_{1}$, we have:

$$ $\displaystyle P_{1}$ $\displaystyle=4\mathbb{E}\|\nabla f(\mathbf{X}_{1})-\mathbf{M}_{1}\|_{F}^{2}$ $\displaystyle=4\mathbb{E}\|\nabla f(\mathbf{X}_{1})-(1-\beta_{1}+\gamma_{1}\beta_{1})\nabla f(\mathbf{X}_{1};\xi_{1})\|_{F}^{2}$ $\displaystyle=4\mathbb{E}\|\nabla f(\mathbf{X}_{1})-\nabla f(\mathbf{X}_{1};\xi_{1})\|_{F}^{2}$ $\displaystyle\stackrel{{\scriptstyle(\star)}}{{\leq}}4\sigma^{2}.$ $$

The final inequality $(\star)$ holds by setting the parameter $\gamma_{1}=\gamma=1$. This choice nullifies the first term, as $(\gamma_{1}-1)^{2}=0$, and simplifies the coefficient of the variance to $(1-\beta_{1}+\beta_{1})^{2}=1$.
∎

###### Lemma C.3 .

###### Proof.

## Appendix D Proofs of Theorem 3.2

###### Proof.

According to Assumption [3.3](https://arxiv.org/html/2509.15816v3#S3.Thmassumption3), and based on the descent lemma, we have

$$ $\displaystyle f(\mathbf{X}_{t+1})$ $\displaystyle\leq f(\mathbf{X}_{t})+\langle\nabla f(\mathbf{X}_{t}),\mathbf{X}_{t+1}-\mathbf{X}_{t}\rangle+\frac{L}{2}\|\mathbf{X}_{t+1}-\mathbf{X}_{t}\|_{F}^{2}$ $\displaystyle\leq f(\mathbf{X}_{t})+\langle\mathbf{M}_{t},\mathbf{X}_{t+1}-\mathbf{X}_{t}\rangle+\langle\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t},\mathbf{X}_{t+1}-\mathbf{X}_{t}\rangle+\frac{L}{2}\|\mathbf{X}_{t+1}-\mathbf{X}_{t}\|_{F}^{2}$ $\displaystyle\stackrel{{\scriptstyle(\circ)}}{{\leq}}f(\mathbf{X}_{t})-\eta_{t}\|\mathbf{M}_{t}\|_{*}+\langle\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t},\mathbf{X}_{t+1}-\mathbf{X}_{t}\rangle+\frac{L}{2}\|\mathbf{X}_{t+1}-\mathbf{X}_{t}\|_{F}^{2}$ $\displaystyle\leq f(\mathbf{X}_{t})-\eta_{t}\|\mathbf{M}_{t}\|_{*}+\frac{1}{2\alpha}\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}+\frac{\alpha+L}{2}\|\mathbf{X}_{t+1}-\mathbf{X}_{t}\|_{F}^{2}$ $\displaystyle\stackrel{{\scriptstyle(\star)}}{{\leq}}f(\mathbf{X}_{t})-\eta_{t}\|\mathbf{M}_{t}\|_{F}+\frac{\sqrt{\eta_{t}}}{2L}\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}+\frac{\frac{L}{\sqrt{\eta_{t}}}+L}{2}\|\mathbf{X}_{t+1}-\mathbf{X}_{t}\|_{F}^{2},$ $$

where $(\circ)$ holds because, by the definition of $\mathbf{O}_{t}$ and the property of norms $\|\mathbf{M}_{t}\|_{*}\geq\|\mathbf{M}_{t}\|_{F}$, we have:

$$ $\langle\mathbf{M}_{t},\mathbf{X}_{t+1}-\mathbf{X}_{t}\rangle=-\eta_{t}\langle\mathbf{M}_{t},\mathbf{O}_{t}\rangle=-\eta_{t}\langle\mathbf{M}_{t},\mathbf{U}_{r}\mathbf{V}_{r}^{\top}\rangle=-\eta_{t}\|\mathbf{M}_{t}\|_{*}\leq-\eta_{t}\|\mathbf{M}_{t}\|_{F},$ $$

and $(\star)$ holds by setting $\alpha=\frac{L}{\sqrt{\eta_{t}}}$.
Thus, we have

$$ $\displaystyle\sum_{t=1}^{T}\eta_{t}\mathbb{E}\|\mathbf{M}_{t}\|_{F}$ $\displaystyle\leq\sum_{t=1}^{T}\left(\mathbb{E}[f(\mathbf{X}_{t})]-\mathbb{E}[f(\mathbf{X}_{t+1})]\right)+\sum_{t=1}^{T}\frac{\sqrt{\eta_{t}}}{2L}\mathbb{E}\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}$ $\displaystyle\quad+\sum_{t=1}^{T}\frac{\frac{L}{\sqrt{\eta_{t}}}+L}{2}\mathbb{E}\|\mathbf{X}_{t+1}-\mathbf{X}_{t}\|_{F}^{2}$ $\displaystyle\leq f(\mathbf{X}_{1})-f^{*}+\sum_{t=1}^{T}\frac{\sqrt{\eta_{t}}}{2L}\mathbb{E}\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}+\sum_{t=1}^{T}\frac{L(\eta_{t}^{3/2}+\eta_{t}^{2})}{2}n$ $\displaystyle\leq f(\mathbf{X}_{1})-f^{*}+\frac{1}{2L}\sum_{t=1}^{T}\sqrt{\eta_{t}}\mathbb{E}\|\mathbf{S}_{t}\|_{F}^{2}+\frac{Ln}{2}\sum_{t=1}^{T}t^{-1}+\frac{Ln}{2}\sum_{t=1}^{T}t^{-4/3}$ $\displaystyle\stackrel{{\scriptstyle(\circ)}}{{\leq}}f(\mathbf{X}_{1})-f^{*}+\frac{4\sigma^{2}+8(2L^{2}n+\sigma^{2})(1+\ln T)}{2L}$ $\displaystyle\quad+\frac{Ln}{2}(1+\ln T)+\frac{Ln}{2}\sum_{t=1}^{T}t^{-4/3}$ $\displaystyle\stackrel{{\scriptstyle(\star)}}{{\leq}}f(\mathbf{X}_{1})-f^{*}+\frac{2\sigma^{2}}{L}+2Ln+4(2Ln+\sigma^{2}L^{-1})(1+\ln T)+\frac{Ln}{2}(1+\ln T),$ $$

where $(\circ)$ uses Lemma [C.2](https://arxiv.org/html/2509.15816v3#A3.Thmlemma2); $(\star)$ follows from the fact that $\sum_{t=1}^{T}\frac{1}{t^{4/3}}\leq 4$.
Next, we let

$$ $G=f(\mathbf{X}_{1})-f^{*}+\frac{2\sigma^{2}}{L}+\left(8Ln+4\sigma^{2}L^{-1}+\frac{Ln}{2}\right)(1+\ln T)+2Ln.$ $$

Thus, we have

$$ $\displaystyle\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\|\mathbf{M}_{t}\|_{F}$ $\displaystyle\leq\frac{1}{T}\sum_{t=1}^{T}\frac{t^{2/3}}{t^{2/3}}\mathbb{E}\|\mathbf{M}_{t}\|_{F}$ $\displaystyle\leq\frac{T^{2/3}}{T}\sum_{t=1}^{T}\frac{1}{t^{2/3}}\mathbb{E}\|\mathbf{M}_{t}\|_{F}=\frac{1}{T^{1/3}}\sum_{t=1}^{T}\eta_{t}\mathbb{E}\|\mathbf{M}_{t}\|_{F}$ $\displaystyle\leq\frac{G}{T^{1/3}}.$ $$

Next, we have

$$ $\displaystyle\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}$ $\displaystyle\stackrel{{\scriptstyle(\circ)}}{{\leq}}\sqrt{\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}}$ $\displaystyle\stackrel{{\scriptstyle(\star)}}{{\leq}}\sqrt{\frac{4\sigma^{2}+8(2L^{2}n+\sigma^{2})(1+\ln T)}{T^{2/3}}}.$ $$

where $(\circ)$ uses Jensen’s inequality; $(\star)$ uses Lemma [C.3](https://arxiv.org/html/2509.15816v3#A3.Thmlemma3) by letting $\mathbf{S}_{t}=\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}$.
Thus, we have

$$ $\displaystyle\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\|\nabla f(\mathbf{X}_{t})\|_{F}$ $\displaystyle\leq\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}+\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\|\mathbf{M}_{t}\|_{F}$ $\displaystyle\leq\sqrt{\frac{4\sigma^{2}+8(2L^{2}n+\sigma^{2})(1+\ln T)}{T^{2/3}}}+\frac{G}{T^{1/3}}$ $\displaystyle=\mathcal{O}\left(\frac{\ln T}{T^{1/3}}\right).$ $$

This completes the proof.
∎

###### Proof.

## Appendix E Lemmas for Theorems 3.3 and 3.4

### E.1 Lemma E.1

###### Lemma E.1 .

Let $\{\Gamma_{t}\}_{t\geq 1}$ be a non-negative sequence whose partial sum $D_{t}=\sum_{t=1}^{T}\Gamma_{t}$ is bounded by $D_{t}\leq A_{1}\ln T+A_{2}$ for some positive constants $A_{1}$ and $A_{2}$. Define the function $f(t)=\frac{t^{1-p}}{\ln t}$ where $p\in(0,1)$. Fix $T_{0}=2e^{1/(1-p)}$, where $e$ is the base of the natural logarithm. Then, the following inequality holds for all $T>T_{0}$:

$$ $\sum_{t=T_{0}}^{T-1}\Gamma_{t}f(t)\leq\left(A_{1}+A_{2}\right)(T-1)^{1-p}.$ $$

###### Proof.

We choose an integer $T_{0}=2e^{1/(1-p)}>e^{1/(1-p)}$ large enough such that for all $t\geq T_{0}$, the function $f(t)$ is positive and monotonically increasing.
We use summation by parts to bound the sum $\sum_{t=T_{0}}^{T-1}\Gamma_{t}f(t)$:

$$ $\sum_{t=T_{0}}^{T-1}\Gamma_{t}f(t)=\sum_{t=T_{0}}^{T-1}(D_{t}-D_{t-1})f(t)=D_{T-1}f(T-1)-D_{T_{0}-1}f(T_{0})-\sum_{t=T_{0}}^{T-1}D_{t-1}(f(t)-f(t-1)).$ $$

Since $D_{t-1}\geq 0$ and $f(t)$ is monotonically increasing for $t\geq T_{0}$ (i.e., $f(t)-f(t-1)\geq 0$), the last term, $\sum_{t=T_{0}}^{T-1}D_{t-1}(f(t)-f(t-1))$, is non-negative. Therefore, we have:

$$ $\sum_{t=T_{0}}^{T-1}\Gamma_{t}f(t)\leq D_{T-1}f(T-1).$ $$

Substituting the given bound on the partial sums, $D_{T-1}\leq A_{1}\ln(T-1)+A_{2}$, we obtain:

$$ $\displaystyle\sum_{t=T_{0}}^{T-1}\Gamma_{t}f(t)$ $\displaystyle\leq(A_{1}\ln(T-1)+A_{2})\frac{(T-1)^{1-p}}{\ln(T-1)}$ $\displaystyle=\left(A_{1}+\frac{A_{2}}{\ln(T-1)}\right)(T-1)^{1-p}$ $\displaystyle\stackrel{{\scriptstyle(\circ)}}{{\leq}}(A_{1}+A_{2})(T-1)^{1-p},$ $$

where $(\circ)$ follows from $\ln(T-1)\geq 1$ for sufficiently large $T$.

Consequently, for any given $\epsilon>0$, the expression can be bounded by $(A_{1}+A_{2}+\epsilon)T^{1-p}$ for sufficiently large $T$.
∎

###### Lemma E.1 .

###### Proof.

### E.2 Lemma E.2

###### Lemma E.2 .

Suppose a positive sequence $\{\Delta_{t}\}_{t\geq 1}$ satisfies the following recursive inequality:

$$ $\Delta_{t+1}\leq\Delta_{t}-\frac{\sqrt{2\mu}}{t^{p}}\sqrt{\Delta_{t}}+\Gamma_{t},$ $$

where $p\in(0,1)$ and $\mu>0$ are constants. The non-negative noise sequence $\{\Gamma_{t}\}_{t\geq 1}$ satisfies the condition from Lemma [E.1](https://arxiv.org/html/2509.15816v3#A5.Thmlemma1) that $\sum_{t=1}^{T}\Gamma_{t}=A_{1}\ln T+A_{2}$. Then, the sequence $\Delta_{t}$ converges to 0 with the following rate:

$$ $\Delta_{t}\leq\frac{\mathcal{A}^{2}(1-p)^{2}}{C^{2}}\cdot\frac{(\ln t)^{2}}{t^{2(1-p)}},$ $$

where $\mathcal{A}:=A_{1}+A_{2}$ and $C=\frac{\sqrt{2\mu}}{2}$.

###### Proof.

We define the variable $G_{t}=\sqrt{\Delta_{t}}$.

Since $\Delta_{t}>0$, it follows that $G_{t}>0$. The function $f(x)=\sqrt{x}$ is concave, which implies $\sqrt{y}-\sqrt{x}\leq\frac{1}{2\sqrt{x}}(y-x)$ for any $x,y>0$. By setting $y=\Delta_{t+1}$ and $x=\Delta_{t}$, we obtain:

$$ $G_{t+1}-G_{t}\leq\frac{\Delta_{t+1}-\Delta_{t}}{2G_{t}}.$ $$

Substituting the original recurrence relation $\Delta_{t+1}-\Delta_{t}\leq-\frac{\sqrt{2\mu}}{t^{p}}\sqrt{\Delta_{t}}+\Gamma_{t}$ and noting that $\sqrt{\Delta_{t}}=G_{t}$, we have:

$$ $G_{t+1}-G_{t}\leq\frac{1}{2G_{t}}\left(-\frac{\sqrt{2\mu}}{t^{p}}G_{t}+\Gamma_{t}\right).$ $$

Letting $C=\frac{\sqrt{2\mu}}{2},T_{0}=2e^{1/(1-p)}$, we arrive at the core recurrence for $G_{t}$:

$$ $G_{t+1}\leq G_{t}-\frac{C}{t^{p}}+\frac{\Gamma_{t}}{2G_{t}}$ (18) $$

For $p\in(0,1)$, we will prove that $G_{t}\leq\frac{\mathcal{A}(1-p)}{C}\cdot\frac{\ln t}{t^{1-p}}$ by contradiction.

Assume the proposition is false. This means that for any constant $A>0$, there exist infinitely many time steps $t$ such that $G_{t}>B_{t}$, where $B_{t}=A\cdot\frac{\ln t}{t^{1-p}}$.

We choose a constant $A$ large enough (the specific condition will be derived later). By our assumption, there must exist an arbitrarily large time $T$ such that $G_{T}>B_{T}$. Consider an interval $[T_{0},T-1]$ and assume that $G_{t}>B_{t}$ for all $t$ in this interval.

Summing the inequality ([18](https://arxiv.org/html/2509.15816v3#A5.E18)) from $t=T_{0}$ to $T-1$ gives:

$$ $G_{T}-G_{T_{0}}\leq\sum_{t=T_{0}}^{T-1}(G_{t+1}-G_{t})\leq\sum_{t=T_{0}}^{T-1}\left(-\frac{C}{t^{p}}+\frac{\Gamma_{t}}{2G_{t}}\right).$ $$

Over this interval, since $G_{t}>B_{t}=A\frac{\ln t}{t^{1-p}}$, we have $\frac{1}{G_{t}}<\frac{1}{B_{t}}=\frac{t^{1-p}}{A\ln t}$. Substituting this bound yields:

$$ $G_{T}<G_{T_{0}}-\underbrace{C\sum_{t=T_{0}}^{T-1}\frac{1}{t^{p}}}_{V_{1}}+\underbrace{\frac{1}{2A}\sum_{t=T_{0}}^{T-1}\frac{\Gamma_{t}t^{1-p}}{\ln t}}_{V_{2}}.$ $$

Next, we analyze the asymptotic behavior of the two sums, $V_{1}$ and $V_{2}$.

For the negative drift term $V_{1}$, we use an integral approximation:

$$ $\sum_{t=T_{0}}^{T-1}\frac{1}{t^{p}}\geq\int_{T_{0}}^{T}x^{-p}dx=\frac{1}{1-p}(T^{1-p}-T_{0}^{1-p}).$ $$

The leading term of $V_{1}$ is thus on the order of $\frac{C}{1-p}T^{1-p}$.

For the positive noise term $V_{2}$, we use the condition on the partial sums of $\Gamma_{t}$, $S_{t}=\sum_{k=1}^{t}\Gamma_{k}\leq A_{1}\ln t+A_{2}$. Applying summation by parts (in a manner analogous to Lemma [E.1](https://arxiv.org/html/2509.15816v3#A5.Thmlemma1)), the dominant behavior of the sum in $V_{2}$ is given by $S_{T-1}\cdot\frac{(T-1)^{1-p}}{\ln(T-1)}$. This can be bounded by $((A_{1}+A_{2})\ln(T-1))\frac{(T-1)^{1-p}}{\ln(T-1)}=(A_{1}+A_{2})(T-1)^{1-p}$. For simplicity in asymptotics, let us denote the effective constant as $\mathcal{A}:=A_{1}+A_{2}$. The leading behavior of $V_{2}$ is thus on the order of $\frac{\mathcal{A}}{2A}T^{1-p}$.

Substituting the leading terms of these bounds back into the inequality for $G_{T}$:

$$ $G_{T}<G_{T_{0}}-\frac{C}{1-p}T^{1-p}+\frac{\mathcal{A}}{2A}T^{1-p}=G_{T_{0}}+T^{1-p}\left(\frac{\mathcal{A}}{2A}-\frac{C}{1-p}\right).$ $$

We now choose the constant $A$ such that the coefficient of the leading term is negative:

$$ $\frac{\mathcal{A}}{2A}-\frac{C}{1-p}<0\implies A>\frac{\mathcal{A}(1-p)}{2C}.$ $$

With an $A$ satisfying this condition, for a sufficiently large $T$, the negative term proportional to $T^{1-p}$ will dominate, forcing the right-hand side to become negative. This contradicts the fact that $G_{T}=\sqrt{\Delta_{T}}$ must be positive.

This contradiction proves that our assumption that $G_{t}$ can remain above $B_{t}$ for an arbitrarily long interval is false. Therefore, there must exist a time $T_{A}$ such that for all $t>T_{A}>T_{0}=2e^{1/(1-p)}$, $G_{t}\leq B_{t}$. This establishes that $G_{t}=\mathcal{O}\left(\frac{\ln t}{t^{1-p}}\right)$.

Finally, by choosing the constant $A=\frac{\mathcal{A}(1-p)}{C}$ , squaring both sides of this result yields the convergence rate for $\Delta_{t}$:

$$ $\displaystyle\Delta_{t}=G_{t}^{2}$ $\displaystyle\leq\frac{\mathcal{A}^{2}(1-p)^{2}}{C^{2}}\cdot\frac{(\ln t)^{2}}{t^{2(1-p)}}$ $\displaystyle=\mathcal{O}\left(\frac{(\ln t)^{2}}{t^{2(1-p)}}\right).$ $$

∎

###### Lemma E.2 .

###### Proof.

## Appendix F Proofs of Theorem 3.3

###### Proof.

By Theorem [3.1](https://arxiv.org/html/2509.15816v3#S3.Thmtheorem1) inequality ([7](https://arxiv.org/html/2509.15816v3#A2.E7)), we have

$$ $\mathbb{E}[f(\mathbf{X}_{t+1})]\leq\mathbb{E}[f(\mathbf{X}_{t})]-\eta_{t}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})\|_{F}]+\Gamma_{t}$ (19) $$

where

$$ $\Gamma_{t}=\frac{\eta_{t}^{2/3}}{L}\mathbb{E}\|\mathbf{M}_{t}-\nabla f(\mathbf{X}_{t})\|_{F}^{2}+(L/2+Ln)\eta_{t}^{4/3}.$ $$

Case 1: $\gamma=0$.

Let $\eta_{t}=t^{-3/4},\beta_{t}=1-t^{-1/2},\alpha_{t}=t^{-1/2}$. According to Theorem [3.1](https://arxiv.org/html/2509.15816v3#S3.Thmtheorem1) inequality ([9](https://arxiv.org/html/2509.15816v3#A2.E9)), we have

$$ $\displaystyle\sum_{t=1}^{T}\Gamma_{t}\leq A_{1}\ln T+A_{2},$ $$

where $A_{1}=2L^{-1}\sigma^{2}+4\sqrt{2}Ln+Ln+L/2$ and $A_{2}=4L^{-1}\sigma^{2}+4\sqrt{2}Ln+Ln+L/2$.

Case 2: $\gamma\neq 0$.

Let $\eta_{t}=t^{-3/4},\beta_{t}=1-(1+t)^{-1/2},\alpha_{t}=t^{-1/2}$. According to Theorem [3.1](https://arxiv.org/html/2509.15816v3#S3.Thmtheorem1) inequality ([13](https://arxiv.org/html/2509.15816v3#A2.E13)), we have

$$ $\displaystyle\sum_{t=1}^{T}\Gamma_{t}\leq A_{1}\ln T+A_{2},$ $$

where $A_{1}=4L^{-1}\sigma^{2}+8\sqrt{2}Ln+Ln+L/2$ and $A_{2}=10L^{-1}\sigma^{2}+8\sqrt{2}Ln+Ln+L/2$.

Then, we define $\Delta_{t}=\mathbb{E}[f(\mathbf{X}_{t})]-f^{*}$ . Applying Assumption [3.5](https://arxiv.org/html/2509.15816v3#S3.Thmassumption5) to inequality ([19](https://arxiv.org/html/2509.15816v3#A6.E19)), we have

$$ $\displaystyle\Delta_{t+1}$ $\displaystyle\leq\Delta_{t}-\eta_{t}\sqrt{2\mu\Delta_{t}}+\Gamma_{t}$ $\displaystyle=\Delta_{t}-\frac{\sqrt{2\mu}}{t^{3/4}}\sqrt{\Delta_{t}}+\Gamma_{t}.$ $$

By Lemma [E.2](https://arxiv.org/html/2509.15816v3#A5.Thmlemma2), let $G_{t}=\sqrt{\Delta_{t}}$, $C=\sqrt{2\mu}/2$ and $p=3/4$.

$$ $G_{t+1}\leq G_{t}-\frac{C}{t^{3/4}}+\frac{\Gamma_{t}}{2G_{t}}.$ $$

Setting $p=3/4$, for any constant $A>\frac{(A_{1}+A_{2})(1-p)}{2C}=\frac{A_{1}+A_{2}}{8C}$, there exists a time step $T_{A}$ such that for all $t>T_{A}$, we have $G_{t}\leq A\cdot\frac{\ln t}{t^{1/4}}$. By choosing $A=\frac{\mathcal{A}}{4C}$, we satisfy the condition $A>\frac{\mathcal{A}}{8C}$, which ensures this bound holds. Thus, we obtain the convergence rate for $\Delta_{t}$:

$$ $\Delta_{t}\leq\frac{\mathcal{A}^{2}(1-3/4)^{2}}{C^{2}}\cdot\frac{(\ln t)^{2}}{t^{2(1-3/4)}}=\frac{\mathcal{A}^{2}}{8\mu}\cdot\frac{(\ln t)^{2}}{t^{1/2}}$ $$

We can now summarize the results for both scenarios. In either case ($\gamma=0$ or $\gamma\neq 0$), the analysis yields the same asymptotic convergence rate. For any number of iterations $T>T_{0}=2e^{4}$, the analysis leads to the same asymptotic upper bound:

$$ $\mathbb{E}[f(\mathbf{X}_{T+1})]-f^{*}\leq\frac{\mathcal{A}^{2}}{8\mu}\cdot\frac{(\ln T)^{2}}{T^{1/2}}=\mathcal{O}\left(\frac{(\ln T)^{2}}{T^{1/2}}\right).$ $$

The specific definition of $\mathcal{A}$ is as follows:
(i) Case 1 ($\gamma=0$): $\mathcal{A}=6L^{-1}\sigma^{2}+8\sqrt{2}Ln+2Ln+L$. (ii) Case 2 ($\gamma\neq 0$): $\mathcal{A}=14L^{-1}\sigma^{2}+16\sqrt{2}Ln+2Ln+L$.

This completes the proof.
∎

###### Proof.

## Appendix G Proofs of Theorem 3.4

###### Proof.

According to Theorem [3.2](https://arxiv.org/html/2509.15816v3#S3.Thmtheorem2), we have

$$ $\displaystyle f(\mathbf{X}_{t+1})\leq f(\mathbf{X}_{t})-\eta_{t}\|\mathbf{M}_{t}\|_{F}+\frac{\sqrt{\eta_{t}}}{2L}\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}+\frac{\frac{L}{\sqrt{\eta_{t}}}+L}{2}\|\mathbf{X}_{t+1}-\mathbf{X}_{t}\|_{F}^{2}.$ $$

Since $\eta_{t}\leq 1$, we have $\frac{L}{\sqrt{\eta_{t}}}\geq L$. Thus, taking the expectation yields

$$ $\displaystyle\mathbb{E}[f(\mathbf{X}_{t+1})]$ $\displaystyle\leq\mathbb{E}[f(\mathbf{X}_{t})]-\eta_{t}\mathbb{E}[\|\mathbf{M}_{t}\|_{F}]+\frac{\sqrt{\eta_{t}}}{2L}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}]+\frac{L}{\sqrt{\eta_{t}}}\mathbb{E}[\|\mathbf{X}_{T+1}-\mathbf{X}_{t}\|_{F}^{2}]$ $\displaystyle\leq\mathbb{E}[f(\mathbf{X}_{t})]-\eta_{t}\mathbb{E}[\|\mathbf{M}_{t}\|_{F}]+\frac{\sqrt{\eta_{t}}}{2L}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}]+Ln\eta_{t}^{3/2}$ $\displaystyle\stackrel{{\scriptstyle(\circ)}}{{\leq}}\mathbb{E}[f(\mathbf{X}_{t})]-\eta_{t}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})\|_{F}]+\eta_{t}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}]$ $\displaystyle+\frac{\sqrt{\eta_{t}}}{2L}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}]+Ln\eta_{t}^{3/2}$ $\displaystyle\stackrel{{\scriptstyle(\star)}}{{\leq}}\mathbb{E}[f(\mathbf{X}_{t})]-\eta_{t}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})\|_{F}]+\frac{1}{2\epsilon}\eta_{t}^{2}+\frac{\epsilon}{2}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}]$ $\displaystyle+\frac{\sqrt{\eta_{t}}}{2L}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}]+Ln\eta_{t}^{3/2}$ $\displaystyle\stackrel{{\scriptstyle(\bullet)}}{{=}}\mathbb{E}[f(\mathbf{X}_{t})]-\eta_{t}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})\|_{F}]$ $\displaystyle+\underbrace{\frac{L\eta_{t}^{3/2}}{2}+\left(\frac{\epsilon}{2}+\frac{\sqrt{\eta_{t}}}{2L}\right)\mathbb{E}[\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}^{2}]+Ln\eta_{t}^{3/2}}_{\Gamma_{t}},$ $$

where $(\circ)$ follows from the reverse triangle inequality $-\|\mathbf{M}_{t}\|_{F}\leq\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}-\|\nabla f(\mathbf{X}_{t})\|_{F}$; $(\star)$ applies Young’s inequality to the term $\eta_{t}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}\|_{F}]$; and $(\bullet)$ collects the residual terms into $\Gamma_{t}$. By setting $\epsilon=\frac{\sqrt{\eta_{t}}}{L}$, we can simplify the expression for $\Gamma_{t}$:

$$ $\Gamma_{t}=\frac{\sqrt{\eta_{t}}}{L}\mathbb{E}[\|\mathbf{S}_{t}\|_{F}^{2}]+(L/2+Ln)\eta_{t}^{3/2},$ $$

where $\mathbf{S}_{t}=\nabla f(\mathbf{X}_{t})-\mathbf{M}_{t}$. The sum of these terms can be bounded. Let $A_{1}=17Ln+8L^{-1}\sigma^{2}+L/2$ and $A_{2}=17Ln+12L^{-1}\sigma^{2}+L/2$. We have

$$ $\displaystyle\sum_{t=1}^{T}\Gamma_{t}$ $\displaystyle=\frac{1}{L}\sum_{t=1}^{T}\sqrt{\eta_{t}}\mathbb{E}[\|\mathbf{S}_{t}\|_{F}^{2}]+\sum_{t=1}^{T}(L/2+Ln)\eta_{t}^{3/2}$ $\displaystyle\leq\frac{1}{L}\sum_{t=1}^{T}\sqrt{\eta_{t}}\mathbb{E}[\|\mathbf{S}_{t}\|_{F}^{2}]+(L/2+Ln)\sum_{t=1}^{T}\frac{1}{t}$ $\displaystyle\stackrel{{\scriptstyle(\circ)}}{{\leq}}\frac{4\sigma^{2}+8(2L^{2}n+\sigma^{2})(1+\ln T)}{L}+(L/2+Ln)(1+\ln T)$ $\displaystyle=A_{1}\ln T+A_{2},$ $$

where $(\circ)$ holds by Lemma [C.3](https://arxiv.org/html/2509.15816v3#A3.Thmlemma3) inequality ([17](https://arxiv.org/html/2509.15816v3#A3.E17)).

Next, let $\Delta_{t}=\mathbb{E}[f(\mathbf{X}_{t})]-f^{*}$. By applying the Assumption [3.5](https://arxiv.org/html/2509.15816v3#S3.Thmassumption5), $\|\nabla f(\mathbf{X}_{t})\|_{F}^{2}\geq 2\mu(f(\mathbf{X}_{t})-f^{*})$, and Jensen’s inequality to the main recurrence, we obtain

$$ $\displaystyle\Delta_{t+1}$ $\displaystyle\leq\Delta_{t}-\eta_{t}\mathbb{E}[\|\nabla f(\mathbf{X}_{t})\|_{F}]+\Gamma_{t}$ $\displaystyle\leq\Delta_{t}-\eta_{t}\sqrt{2\mu\mathbb{E}[f(\mathbf{X}_{t})-f^{*}]}+\Gamma_{t}$ $\displaystyle\leq\Delta_{t}-\eta_{t}\sqrt{2\mu\Delta_{t}}+\Gamma_{t}.$ $$

Setting the step size $\eta_{t}=1/t^{2/3}$, the recurrence becomes

$$ $\Delta_{t+1}\leq\Delta_{t}-\frac{\sqrt{2\mu}}{t^{2/3}}\sqrt{\Delta_{t}}+\Gamma_{t}.$ $$

This recurrence is in the form required by Lemma [E.2](https://arxiv.org/html/2509.15816v3#A5.Thmlemma2). Let $G_{t}=\sqrt{\Delta_{t}}$, $C=\sqrt{2\mu}/2$ and $p=2/3$. According to Lemma [E.2](https://arxiv.org/html/2509.15816v3#A5.Thmlemma2), for any constant $A>\frac{(A_{1}+A_{2})(1-p)}{2C}=\frac{A_{1}+A_{2}}{6C}$, there exists a time $T_{A}$ such that for all $t>T_{A}$, $G_{t}\leq A\cdot\frac{\ln t}{t^{1/3}}$. We choose a specific constant $A=\frac{\mathcal{A}}{3C}$, where $\mathcal{A}:=A_{1}+A_{2}=20L^{-1}\sigma^{2}+34Ln+L$. Squaring both sides of the inequality for $G_{t}$ establishes the convergence rate for $\Delta_{t}=G_{t}^{2}$:

$$ $\Delta_{t}\leq A^{2}\cdot\frac{(\ln t)^{2}}{t^{2/3}}=\frac{\mathcal{A}^{2}}{9C^{2}}\frac{(\ln t)^{2}}{t^{2/3}}=\frac{2\mathcal{A}^{2}}{9\mu}\frac{(\ln t)^{2}}{t^{2/3}}$ $$

Therefore, for any number of iterations $T>T_{0}=2e^{3}$, the expected suboptimality has the following asymptotic bound:

$$ $\displaystyle\mathbb{E}[f(\mathbf{X}_{T+1})]-f^{*}\leq\frac{2\mathcal{A}^{2}}{9\mu}\cdot\frac{(\ln T)^{2}}{T^{2/3}}=\mathcal{O}\left(\frac{(\ln T)^{2}}{T^{2/3}}\right).$ $$

This completes the proof.
∎

###### Proof.

## Appendix H Experimental Details

### H.1 Training on CIFAR10

The ResNet18 model he2016deep undergoes pretraining on the CIFAR-10 dataset with comprehensive hyperparameter specifications provided in Table [2](https://arxiv.org/html/2509.15816v3#A8.T2). For each optimizer, the learning rate is selected via a grid search over the set $\{1\times 10^{-4},5\times 10^{-4},10^{-3},5\times 10^{-3},10^{-2},5\times 10^{-2},10^{-1}\}$. To ensure a robust comparison, all experiments are repeated over five different random seeds, and we report the mean results with one standard deviation shaded.

**Table 2: Hyperparameters used for training ResNet18 on CIFAR10**
|  | SGD | Adam | Muon-Variant |
| --- | --- | --- | --- |
| Model Size | 42.7M |  |  |
| Training Epochs | 100 |  |  |
| Batch Size | 128 |  |  |
| Learning Rate | 0.1 | 0.01 | 0.05 |
| Learning Rate Scheduling | cosine to 10% |  |  |
| Numerical precision | float32 |  |  |
| Weight Decay | 0.01 |  |  |
| $(\beta_{1},\beta_{2})$ | ✗ | (0.9,0.999) | (0.9,-) |
| Gamma | ✗ | ✗ | 0.1 |

### H.2 Pretraining on C4

The LLaMA2-130M modelTouvron2023Llama2O was pre-trained on the C4 (Colossal Clean Crawled Corpus) dataset raffel2020exploring. Identical learning rates and scheduling protocols were systematically implemented across all optimizers during the training process. Comprehensive experimental specifications are tabulated in Table [3](https://arxiv.org/html/2509.15816v3#A8.T3). We use the Muon implementation from Moonlight(^1^11[https://github.com/MoonshotAI/Moonlight](https://github.com/MoonshotAI/Moonlight)).

**Table 3: Hyperparameters used for training LLaMA2-130M on C4**
|  | AdamW | Muon-Variant |
| --- | --- | --- |
| Model Size | 130M |  |
| Hidden Size | 768 |  |
| Head | 12 |  |
| Depth | 12 |  |
| Training Steps | 10000 |  |
| Warmup Steps | 1000 |  |
| Maximum Length | 1024 |  |
| Batch Size | 512 |  |
| Learning Rate | 1e-3 |  |
| Warmup Scheduling | linear from 8e-6 |  |
| Learning Rate Scheduling | cosine to 10% |  |
| Numerical precision | bfloat16 |  |
| Weight Decay | 0.1 |  |
| $(\beta_{1},\beta_{2})$ | (0.9,0.999) | (0.95,-) |
| Gamma | ✗ | 0.05 |

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 4
presentation: 2
questions: - When is Muon-MVR2 the better choice? Fig1 A/B seems to show it's better per-epoch, but it appears slower than Muon-MVR1 in wall-clock time on Fig 1C. Is there a point after N number of epochs where its faster convergence overcomes its 2x per-step cost? A larger experiment is needed to examine this possibility. 

- How well does Muon-MVR2 compare against variance-reduced optimizers? Such as SPIDER[3] and SVRG? 

- The theoretical contributions of the paper largely depend on perfect orthogonalization. But in practice this is an approximation (Newton-Schultz) for efficiency. Does using the approximation break these theoretical guarantees? 

[3] Fang, Cong, et al. "Spider: Near-optimal non-convex optimization via stochastic path-integrated differential estimator." Advances in neural information processing systems 31 (2018).

[4] Gower, Robert M., et al. "Variance-reduced methods for machine learning." Proceedings of the IEEE 108.11 (2020): 1968-1983.
soundness: 3
strengths: - The main contribution of this paper is it's rigorous theoretical proof of optimality for the Muon-class optimizer. There has been some recent interest in understanding the performance of Muon [1], and aiming to close the gap between Muon's practical usability versus its theoretical understanding is valuable. 

- The proposed Muon-MVR2 achieves its optimal convergence using a constant batch size $O(1)$. This is a significant improvement over prior work [2] that needed large initial batch sizes. The last-iterate non-ergodic convergence analysis is also rigorous, showing $O(T^{-2/3})$. 

- There are experiments in both vision and language domains, although the main contribution of the paper remains its theoretical clarification of the Muon framework. There is a clearly expressed tradeoff between the "cheaper" Muon-MVR1 variant and the theoretically backed MVR2 variant. 



[1] Shen, Wei, et al. "On the convergence analysis of muon." arXiv preprint arXiv:2505.23737 (2025).

[2] Sfyraki, Maria-Eleni, and Jun-Kun Wang. "Lions and muons: Optimization via stochastic frank-wolfe." arXiv preprint arXiv:2506.04192 (2025).
summary: This paper presents analysis of two momentum-based variance-reduced Muon-MVR1 and Muon-MVR2 algorithms, which achieve sublinear convergence rates. In addition, the paper also provides convergence analysis details for the standard Muon and Muon-MVR1 methods. Muon-MVR2 variant offers some additional analysis on achieving optimal iteration complexity of $O(T^{-1/3})$ for non-convex optimization. Primary comparisons are between Muon-MVR2 versus the theoretical lower bound. The paper presents some small-scale experiments (such as with CIFAR-10 and LLaMA2), and practical results imply that those looking for the best theoretical convergence rate should use the Muon-MVR2 variant.
weaknesses: - The most concerning point is that the proposed Muon-MVR2 variant requires 2 gradient evaluations per step. This practically introduces almost double the amount of computational overhead. 

- This practical-vs-theory tradeoff is shown in the paper's own experiments. The authors transparently note that the overhead can be prohibitive, especially visible in Fig 1C where the theoretically slower MVR1 is actually faster in practice per wall-clock time. Since the more practically usable MVR1 variant is stuck at the suboptimal $O(T^{-1/4})$ convergence rate, this means most people will likely use this lower-cost variant. Making the implications of the paper's main headline results less impactful. 

- Experiments are small-scale. CIFAR-10 is an outdated and small dataset, but more critically LLaMA2-130M is only trained for 10,000 steps. At this very small scale, it is unclear how the benefits outweigh the costs of Muon-MVR2. The variance reduction technique also introduces new hyperparameters ($\gamma$), which becomes another expensive step to tune in any practical setting.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: 1.	Why wasn't the wall-clock time shown in Figure 2? While O-notation is important, what is practically significant is essentially the wall-clock time.

2.	Muon is primarily known as an optimizer frequently used for LLMs. Therefore, I question the choice of using an older model from 2016 (ResNet-18) specifically for image classification experiments. If the goal was to demonstrate its effectiveness in the vision domain as well, why weren't large-scale datasets like ImageNet/MS COCO and more recent vision models such as CLIP or DINO employed? The chosen combination of ResNet18 & CIFAR-10 might inadvertently lead readers to suspect that the method failed on more modern and larger-scale vision tasks, and thus only results from this simpler setup were presented.
soundness: 4
strengths: The primary strength is the first proof that Muon-MVR2 achieves optimal iteration complexity in the stochastic non-convex function. The paper doesn’t stop at presenting theoretical results. It validates the algorithm through experiments on CIFAR-10 and C4 datasets.
summary: This paper aims to establish a rigorous theoretical foundation for the Muon optimizer, which has recently gained attention for its remarkable empirical performance. The authors point out the limitations of existing Muon analyses and propose two variants to address these issues.
weaknesses: While it's commendable that this study used a large dataset like C4, the testing was conducted on only one model, LLaMA2-130M. Although the paper's main contribution is theoretical, the lack of evaluation on more than one model and more than one dataset remains its biggest drawback.

---

confidence: 4
contribution: 2
initial_rating: 4
presentation: 3
questions: - Technical novelty: how do the proofs differ materially from normalized‑SGD‑with‑momentum (Cutkosky & Mehta, 2020) and prior Muon analyses; why are these not cited? Please also discuss connections to the OptML MVR notes and parameter‑agnostic results (Yang et al., 2023).
- Positioning vs. prior work: add a careful comparison to Sfyraki & Wang (2025) and include Kovalev (2025); specify which assumptions/estimators.
-  Baselines: add MARS‑AdamW on CIFAR‑10 and C4.
- Empirical details: clarify if Muon is used only for matrix‑shaped parameters while AdamW is used elsewhere; list parameter groups and per‑group learning rates; report weight decay values.
- LLaMA2/C4 setup: justify using a single LR for all optimizers or include at least a small per‑optimizer LR sweep.
- Constants: discuss the dependence on $L,\sigma,n$ in Theorems 3.1–3.2 and whether any factors are tight or can be improved.
soundness: 3
strengths: - The core idea is reasonable and timely: MARS-style variance reduction has worked well in practice, and combining MVR with Muon yields provable benefits, making the approach promising for practical impact.
- The analysis of a more practical, one-batch MVR variant appears novel at the algorithmic level; however, I am uncertain about the theoretical novelty/difficulty of this result.
- The main results are stated cleanly and appear correct.
- The algorithmic setup is presented clearly: MVR1 vs. MVR2, schedules, and update rules.
- The empirical section covers reasonable benchmarks (CIFAR-10/ResNet-18, LLaMA2-130M on C4) and broadly supports the claimed advantages, especially for MVR2.
summary: This work integrates Momentum Variance Reduction (MVR)—including a MARS‑style $\gamma$ factor—into Muon. It analyzes three variants (standard Muon/MVR1 with $\gamma=0$, MVR1 with $\gamma_t=t^{-1/2}$, and MVR2 with two gradients on the same mini‑batch), proving ergodic rates $\tilde{\mathcal O}(T^{-1/4})$ for MVR1 and $\tilde{\mathcal O}(T^{-1/3})$ for MVR2, plus PŁ last‑iterate rates $\tilde{\mathcal O}(T^{-1/2})$ (MVR1) and $\tilde{\mathcal O}(T^{-2/3})$ (MVR2). Experiments on CIFAR‑10/ResNet‑18 and LLaMA2‑130M/C4 are provided.
weaknesses: - Despite its practical motivation, the paper’s theoretical analysis appears to offer limited novelty. The proof technique closely mirrors that of normalized SGD with momentum (Cutkosky & Mehta, 2020)—which is, surprisingly, not cited—and is similar to prior analyses of Muon’s convergence. Moreover, MVR in this context has already been discussed (e.g., the OptML notes: https://optmlclass.github.io/notes/optforml_notes.pdf), and parameter-agnostic results under comparable assumptions (e.g., Yang et al., 2023) are not addressed. As a result, beyond analyzing a more practical $\gamma$-factor variant, the contribution seems incremental: with a decreasing schedule, $\gamma_t$, the analysis does not improve the bound over $\gamma=0$ (it essentially recovers the SGD rate). The paper should add a clear, dedicated discussion of technical novelty, explicitly situating its results relative to these works and clarifying what is genuinely new in the Muon-specific setting.

- The $\gamma$‑factor with a decreasing schedule ($\gamma_t=t^{-1/2}$) recovers the same $\tilde{\mathcal O}(T^{-1/4})$ rate as $\gamma=0$, so the analysis does not show a theoretical gain from the MARS‑like correction.
- Related‑work positioning: compares against relatively weak baselines and omits important works such as Kovalev (2025); needs a more detailed comparison to Sfyraki & Wang (2025).
- Optimality beyond the $T$‑exponent is unclear: the dependence on $L,\sigma,n$ (especially for Theorem 3.2) is not discussed relative to known parameter‑agnostic results.
- Missing important baseline: MARS‑AdamW should be included.
- Missing/unclear experimental details: whether Muon is combined with AdamW for some parameter groups, and the separate learning rates; absence of ablations on $\gamma$, weight decay, and batch size; on C4, a single fixed LR is used for all optimizers, which is hard to justify.
- “Nesterov‑accelerated” wording for Option MVR1 is incorrect for the presented update.

---
- Cutkosky, A., & Mehta, H. “Momentum improves normalized SGD.” ICML, 2020.
- Yang, Junchi, et al. “Two sides of one coin: the limits of untuned SGD and the power of adaptive methods.” NeurIPS 36, 2023.
- Kovalev, Dmitry. “Understanding gradient orthogonalization for deep learning via non-Euclidean trust‑region optimization.” arXiv:2503.12645, 2025.
- OptML notes on MVR in this context: https://optmlclass.github.io/notes/optforml_notes.pdf￼

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 2
questions: The main improvement of Muon lies in its orthogonalization step. Could the authors elaborate on why orthogonalization improves optimization dynamics over standard SGD? Is there a formal argument or empirical evidence that the equalization of singular values accelerates convergence?
soundness: 3
strengths: 1. **Clear theoretical contribution.** The paper provides a rigorous convergence analysis for Muon-MVR, bridging the gap between Muon’s empirical success and previously incomplete theory. The Lyapunov-style recursion and variance-control arguments are clean and easy to follow.

2. **Well-structured and consistent analysis.** The logical progression from assumptions to lemmas and theorems is smooth, and the proofs are well-motivated and technically precise. The framework appears general enough to extend to other momentum–variance-reduction methods.

3. **Readable exposition.** Despite the mathematical density, the writing is accessible and the notation—aside from minor issues—remains consistent. The paper is pleasant to read, and the proofs are presented with notable clarity.
summary: This paper analyzes the theoretical convergence of the Muon optimizer and introduces two momentum-based variance-reduced variants, Muon-MVR1 and Muon-MVR2. The authors show that Muon-MVR2 attains the optimal stochastic non-convex iteration complexity of $\tilde{\mathcal{O}}(T^{-1/3})$ and establish last-iterate convergence under the Polyak–Łojasiewicz condition. The analysis is elegant and unified: it extends the stochastic momentum framework to the orthogonalized Muon setting, with proofs that are logically consistent and mathematically transparent. Empirical results on CIFAR-10 and C4 further support the theoretical claims.
weaknesses: 1. **Overstated claims in presentation.** The statement that the method “overcomes the bottleneck of traditional momentum methods, achieving the optimal iteration complexity” is exaggerated. Sfyraki & Wang (2025) already analyze a Muon+VR variant (under a stochastic Frank–Wolfe formulation) achieving the same $\tilde{\mathcal{O}}(T^{-1/3})$ rate. The novelty here lies primarily in the unconstrained SGD-style analysis, not in attaining a new theoretical complexity.

2. **Ambiguity in Table 1.** The superscripts $c$ and $d$ are undefined. Their meaning should be clearly explained in the caption or a footnote.

3. **Formula inconsistency (Eq. (4)).** As written,
   $$
   M_t = \mu M_{t-1} + \nabla f(X_t, \xi_t) + \mu \big(\nabla f(X_t, \xi_t) - \nabla f(X_{t-1}, \xi_{t-1})\big)
   $$
   is not equivalent to Eq. (2). It should instead be
   $$
   C_t = \mu C_{t-1} + (1-\mu)\nabla f(X_t,\xi_t), \quad
   M_t = \mu C_t + (1-\mu)\nabla f(X_t,\xi_t).
   $$

4. **Possible overclaim in Remark 3.2.** Although the analysis avoids a growing-batch assumption, the comparison with Sfyraki & Wang (2025) is not fully aligned, as their problem setting, objective, and convergence metric differ. Characterizing this as a “significant theoretical advantage” is therefore overstated—especially without empirical evidence that fixed-batch training achieves the claimed $\tilde{\mathcal{O}}(T^{-1/3})$ scaling.

5. **Experimental limitations.** The empirical section primarily plots standard training curves and does not directly test the predicted scaling laws. Moreover, there is no comparison against other Muon-type optimizers. As a result, the experiments feel empirical rather than diagnostic.

6. **Typographical and constant issues.**
   - (1) Page 15 (L799): should read $\epsilon = \frac{1-\beta_{t+1}}{\beta_{t+1}}$.
   - (2) Page 15 (L801), Eq. (6): the right-hand side contains an extra “$t$”.
   - (3) Page 22 (L1187): the second term’s constant should be $2\eta_{t+1}^2(4L^2 n + \sigma^2)$; assuming $\eta_t \le 2\eta_{t+1}$, we have $2\eta_t^2 L^2 n \le 8\eta_{t+1}^2 L^2 n$, which affects multiple subsequent constants.

---

confidence: 4
contribution: 2
initial_rating: 2
presentation: 3
questions: - I understand that $\gamma$ is used to control the variance reduction term, as in [1]. However, this scheme is different from the originally proposed STORM scheme in [4]. I also see that in practice $\gamma$ is set to a very small value, which effectively reduces the impact of the variance reduction. Can the authors explain the motivation for including $\gamma$?
- Can the authors elaborate on how equation (4) is obtained from Algorithm 1?



[1] MARS: Unleashing the Power of Variance Reduction for Training Large Models - Yuan et al.

[4] Momentum-Based Variance Reduction in Non-Convex SGD - Cutkosky and Orabona
soundness: 2
strengths: - The paper is for the most part clearly written and its structure is easy to follow.
- The last-iterate convergence results under the PL condition are interesting.
- The convergence guarantees of the approximate variance reduction scheme (MVR1) without the additive error-term improve previous results.
- The experimental results showcase the improved performance of the proposed algorithms alongside other popular optimizers.
summary: This work studies the convergence properties of the Muon optimizer and two variants under suitable assumptions. In particular, for the stochastic non-convex setting, the authors provide convergence guarantees for the standard Muon algorithm and demonstrate improved guarantees for two variance-reduced versions. Additionally, the convergence of Muon and its two variants is examined under the Polyak-Lojasiewicz condition. Finally, the authors present experimental results on vision and language tasks, evaluating Muon and its two variants alongside other popular optimizers.
weaknesses: - The proposed algorithmic variants of Muon are not novel. Prior work [1] has proposed Shampoo (a generalization of Muon) with both one-batch and two-batch momentum variance reduction and also uses the additional $\gamma$ parameter to control the variance-reduction term. In addition, the algorithm in [1] uses weight decay. However, it is important to note that [1] does not provide convergence guarantees for their algorithm.
- The convergence of standard Muon in the non-convex setting has been studied in previous works, as also mentioned in the paper. For example, [2] shows an iteration complexity of $O(T^{-1/4})$ with a batch size of $O(1)$. In contrast, this works shows an iteration complexity of $O(\log(T)/ T^{1/4})$, that is with an additional $\log(T)$ factor. Additionally, [3] also shows that the stochastic first-order oracle calls in Muon are $O( T^{-1/4})$. However, [3] studies Muon with weight decay.
- The convergence of the variant MVR2 (with $\gamma =1$ and additional weight decay) in the non-convex setting was also studied in a previous work [3]. This work claims that their advantage lies in the use of a batch size of $O(1)$, as opposed to [3]. However, [3] demonstrates an interaction complexity of $O(T^{-1/3})$ using a large batch size only in the first iteration and for all subsequent iterations the batch size is set to $O(1)$. Therefore, the amortized average batch size is also $O(1)$. Furthermore, the complexity in [3] does not involve a logarithmic factor, as opposed to this work, which demonstrates a complexity of $O(\log(T)/ T^{1/3})$. Additionally, Assumption 3.3 of this work is stronger than Assumption 2 in [3].
- This work claims that the MVR1 variant is equivalent to Muon with Nesterov momentum for a specific choice of parameters (equation (4)). I do not see how this is possible, though I may be wrong.
- Some parts of the non-ergodic section are a bit unclear. For example, $\Gamma_t$ is not specified in Remark 3.5.
- The introduction of an additional hyperparameter $\gamma$ makes tuning even more tedious in practice.
- The optimizers in the experiments use weight decay, which is not supported by the theoretical results.
- It is not specified how the variance parameter $\gamma$, momentum parameter $\beta$ and weight decay parameter were tuned in the experiments.
- One minor presentation issue is that standard Muon is referred to as Muon-MVR1 with $\gamma=0$, which is a bit misleading, since in that case no variance reduction is applied.


[1] MARS: Unleashing the Power of Variance Reduction for Training Large Models - Yuan et al.

[2] Convergence Bound and Critical Batch Size of Muon Optimizer - Sato et al.

[3] Lions and Muons: Optimization via Stochastic Frank-Wolfe - Sfyraki and Wang


Decision: Reject