Paper deep dive
Convergence guarantees for Muon: New parameter regimes and generalizations
Arthur C. B. de Oliveira, Dhruv D. Jatkar, Guilherme S. Vicinansa, Eduardo D. Sontag
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 89%
Last extracted: 9/28/2026, 3:09:40 AM
Summary
This paper establishes the first asymptotic convergence guarantees for the Muon optimizer by analyzing it as a preconditioned Polyak heavy-ball method. The authors introduce a 'soft-sign' regularization variant (Soft Muon) that ensures a bounded preconditioner, enabling classical Lyapunov analysis to prove convergence to stationary points and linear convergence under the Polyak-Łojasiewicz condition. Additionally, they propose 'Muesterov', a Nesterov-based variant of Muon, and demonstrate similar convergence properties through theoretical proofs and numerical experiments on scalar and transformer-based tasks.
Entities (10)
Relation Signals (8)
Muesterov → isvariantof → Muon
confidence 95% · We formalize this idea by introducing Muesterov, a Nesterov-based variant of Muon...
Soft Muon → isvariantof → Muon
confidence 95% · The soft-sign formulation exposes Muon as a preconditioned heavy-ball method... Soft Muon (8)...
Soft Muon → providesconvergenceguaranteesfor → Muon
confidence 92% · In this paper, we establish the first asymptotic convergence guarantees for a soft-sign variation of the classical formulation of Muon.
Soft Muon → isinterpretedas → Preconditioned Polyak Heavy-Ball
confidence 90% · The soft-sign formulation exposes Muon as a preconditioned heavy-ball method...
Muon → uses → Newton-Schulz iteration
confidence 90% · Muon (Jordan et al., 2024) turns this principle into a practical algorithm by approximating the matrix sign via a small number of Newton-Schulz iterations...
Muesterov → uses → Nesterov Acceleration
confidence 90% · Muesterov, a Nesterov-based variant of Muon...
Soft Muon → satisfies → Polyak-Łojasiewicz Condition
confidence 85% · ...under a global Polyak-Łojasiewicz condition, that the sequence of function values converges linearly.
Muon → isappliedto → nanoGPT
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:In this paper, we establish the first asymptotic convergence guarantees for the Muon algorithm through a more accurate proxy for the Newton-Schultz iteration than the typical matrix sign function. We prove that, for appropriate choices of hyperparameters, the iterates satisfy $\lim_{k\to\infty}\|\nabla f(x_k)\|=0$, and, under a global Polyak-Łojasiewicz condition, that the sequence of function values converges linearly. The key insight is that the regularization, implicit in Muon's Newton-Schulz implementation, induces a bounded preconditioner, exposing Muon as a \emph{preconditioned Polyak heavy-ball} method and enabling a classical Lyapunov analysis. This observation naturally motivates applying the same preconditioning structure to the Nesterov gradient evaluation. We formalize this idea by introducing \emph{Muesterov}, a Nesterov-based variant of Muon, and prove that it enjoys the same convergence guarantees, extending the theoretical framework beyond the heavy-ball setting. Numerical experiments on a scalar cross-entropy problem corroborate the theory and illuminate the joint role of the learning rate and the Newton-Schulz regularizer in controlling convergence. Preliminary numerical simulations training the nanoGPT dataset provide intuition regarding the relevance of the observations in this paper to practical applications.
Tags
Links
- Source: https://arxiv.org/abs/2609.30546v1
- Canonical: https://arxiv.org/abs/2609.30546v1
Trouble viewing inline? Open PDF directly →
Full Text
82,680 characters extracted from source content.
Expand or collapse full text
Convergence guarantees for Muon: New parameter regimes and generalizations Arthur Castello B. de Oliveira Affiliation: Department of Electrical and Computer Engineering Affiliation: Northeastern University Affiliation: Boston, MA 02115 Email: a.castello@northeastern.edu Dhruv D. Jatkar Affiliation: Department of Electrical and Computer Engineering Affiliation: Northeastern University Affiliation: Boston, MA 02115 Email: jatkar.d@northeastern.edu Guilherme S. Vicinansa Affiliation: Department of Telecommunications and Control Affiliation: Escola Politecnica, University of São Paulo Affiliation: São Paulo, SP 05508-010 Email: guilherme.vicinansa@gmail.com Eduardo D. Sontag Affiliation: Department of Electrical and Computer Engineering and Department of Bioengineering Affiliation: Northeastern University Affiliation: Boston, MA 02115 Email: e.sontag@northeastern.edu Abstract In this paper, we establish the first assymptotic convergence guarantees for the Muon algorithm through a more accurate proxy for the Newton-Schultz iteration than the typical matrix sign function. We prove that, for appropriate choices of hyperparameters, the iterates satisfy limk→∞‖∇f(xk)‖=0 _k→∞\|∇ f(x_k)\|=0, and, under a global Polyak-Łojasiewicz condition, that the sequence of function values converges linearly. The key insight is that the regularization, implicit in Muon’s Newton-Schulz implementation, induces a bounded preconditioner, exposing Muon as a preconditioned Polyak heavy-ball method and enabling a classical Lyapunov analysis. This observation naturally motivates applying the same preconditioning structure to the Nesterov gradient evaluation. We formalize this idea by introducing Muesterov, a Nesterov-based variant of Muon, and prove that it enjoys the same convergence guarantees, extending the theoretical framework beyond the heavy-ball setting. Numerical experiments on a scalar cross-entropy problem corroborate the theory and illuminate the joint role of the learning rate and the Newton-Schulz regularizer in controlling convergence. Preliminary numerical simulations training the nanoGPT dataset provide intuition regarding the relevance of the observations in this paper to practical applications. 1 Introduction Training large transformer models has driven sustained interest in optimizers that go beyond the element-wise second-moment estimates of Adam (Kingma and Ba, 2015) and AdamW (Loshchilov and Hutter, 2019). A number of approaches exploit richer structure in the gradient or its second-order statistics. Methods such as K-FAC (Martens and Grosse, 2015) and Sophia (Liu et al., 2024) incorporate curvature information, while algorithms such as Shampoo (Gupta et al., 2018) and SOAP (Vyas et al., 2025) leverage matrix-valued preconditioning and second-moment structure to improve convergence. These methods can yield substantial gains, but typically require maintaining and manipulating additional matrix statistics. A computationally cheaper alternative is to normalise the gradient matrix directly. Bernstein and Newhouse (Bernstein and Newhouse, 2024) formalise this idea through modular duality, showing that for any layer-wise norm, the steepest descent direction is given by its dual map. In the case of the spectral (Schatten-∞) norm, this dual corresponds to the matrix sign function (Higham, 2008), sign(B)=UV⊤sign(B)=UV , where B=UΣV⊤B=U V is an SVD. This can be viewed as a matrix analogue of the element-wise sign used by methods such as signSGD (Bernstein et al., 2018) and Lion (Chen et al., 2023). Muon (Jordan et al., 2024) turns this principle into a practical algorithm by approximating the matrix sign via a small number of Newton-Schulz iterations, which can be implemented efficiently in low precision and parallelised on modern hardware. Subsequent work has studied improved matrix-sign iterations and their numerical properties (Amsel et al., 2025), as well as the empirical performance of Muon in large-scale language model training (Liu et al., 2025). Combined with momentum, this orthogonalisation step leads to updates that can be interpreted as a form of preconditioned momentum iteration, an observation made in several recent theoretical analyses (Li and Hong, 2025; Shen et al., 2025; Sato et al., 2025; Chang et al., 2025). Despite this progress, providing a classical convergence theory for Muon remains challenging. A central difficulty is that the exact matrix sign operator produces updates of essentially constant magnitude: for a rank-r matrix B, one has ‖sign(B)‖F=r\|sign(B)\|_F= r, so the induced preconditioner is not uniformly bounded. This prevents direct application of standard Lyapunov arguments for preconditioned methods. As a result, existing analyses typically establish finite-time or asymptotic stationarity guarantees – often in terms of averaged or expected gradient norms – under smoothness and stochastic assumptions (Li and Hong, 2025; Shen et al., 2025; Sato et al., 2025; Kim and Oh, 2026), with improved rates obtainable under variance reduction (Chang et al., 2025). These results, while valuable, do not yield classical full-sequence convergence guarantees for the iterates themselves. In this paper, we establish the first asymptotic convergence guarantees for a soft-sign variation of the classical formulation of Muon. Our key observation is that the regularization implicit in this variation of Muon’s Newton-Schulz implementation induces a uniformly bounded preconditioner, exposing it as a preconditioned Polyak heavy-ball method and enabling a classical Lyapunov analysis. Under suitable hyperparameter choices, we prove that the iterates satisfy limk→∞‖∇f(xk)‖=0 _k→∞\|∇ f(x_k)\|=0 and, under a global Polyak–Łojasiewicz condition, that the function values converge linearly. Motivated by this interpretation and the known advantages of Nesterov over Polyak heavy-ball (Lessard et al., 2016; Shi et al., 2021), we additionally introduce Muesterov, a Nesterov-based variant of Muon obtained by applying the same normalization structure to extrapolated gradients. We prove that Muesterov satisfies analogous convergence guarantees, and compare it to Muon on scalar synthetic simulations. Finally, we present synthetic and transformer-based training experiments illustrating how the learning rate and the Newton-Schulz regularizer jointly control the transition between oscillatory sign-like dynamics and stable convergence. In particular, the fixed-batch nanoGPT experiments illustrate that the same regularization mechanisms identified by the theory and synthetic simulations remain qualitatively visible in transformer training. All proofs are provided in the technical appendix, and all code used to generate the numerical results is available in Anonymous (2026). 2 Problem Setup and Preliminary Analysis The set of natural and real numbers are denoted ℕN and ℝR respectively. Let ∥⋅∥\|·\| and ⟨⋅,⋅⟩ ·,· denote the Frobenius norm and inner product if taken over matrices, and the ℓ2 _2 norm and inner product if taken over vectors. Given m,n∈ℕm,n and a differentiable function f:ℝm×n→ℝf:R^m× n , let f∗:=minx∈ℝm×nf(x)f^*:= _x ^m× nf(x) denote its minimum value, assumed to be well defined and finite throughout. A function f is said to be proper if the pre-images of compact sets are compact. Furthermore, a function f is said to be L-smooth for some L>0L>0 if for all x,y∈ℝm×nx,y ^m× n it holds that ‖∇f(x)−∇f(y)‖≤L‖x−y‖.\|∇\!f(x)-∇\!f(y)\|≤ L\|x-y\|. (1) Finally, a function f is said to satisfy a global Polyak-Łojasiewicz inequality (PŁI) with constant μ>0μ>0 if for every x∈ℝm×nx ^m× n it holds that ‖∇f(x)‖2≥μ(f(x)−f∗).\|∇\!f(x)\|^2≥μ(f(x)-f^*). (2) For any matrix x∈ℝm×nx ^m× n, let x=UΣV⊤x=U V be a valid singular value decomposition (SVD) of x. Then the matrix sign function (Higham, 2008) is given by sign(x):=UV⊤sign(x):=UV , and is well defined up to the choice of singular vectors corresponding to zero singular values. 2.1 Equivalent Muon formulations We now introduce three different formulations for Muon: the practical implementation as described by Jordan et al. (2024); an idealized sign-based model typically considered in the literature; and a regularized analytical proxy to be used for convergence analysis. We will discuss the equivalence between these three algorithms and why the distinctions matter. For n≥m∈ℕn≥ m , let f:ℝm×n→ℝf:R^m× n be an L-smooth function that satisfies a global Polyak-Łojasiewicz inequality with constant μ>0μ>0. Furthermore, let f be proper and bounded below, with minimum value f∗:=f(x∗)f^*:=f(x^*) attained at some finite point x∗∈ℝm×nx^* ^m× n. Consider, then, the following optimization problem. minimizex∈ℝm×nf(x). x ^m× nminimize f(x). (3) Many methods exist for solving (3), however Muon (Jordan et al., 2024) has rapidly gained widespread popularity for the training of the hidden-layers of transformer models. Muon (henceforth called canonical Muon, or simply Muon) proposes finding the solution of (3) by the following iteration Bk B_k =βBk−1+∇f(xk−1) =β B_k-1+∇\!f(x_k-1) (4a) Ok O_k =(Bk) =N(B_k) (4b) xk x_k =xk−1−αOk, =x_k-1-α O_k, (4c) where (x)N(x) denotes the application of a Newton–Schulz iteration, empirically tuned for Muon applications, to the normalized matrix x/(‖x‖F+δ)x/(\|x\|_F+δ). Prior work (Jordan et al., 2024; Bernstein and Newhouse, 2024; Amsel et al., 2025) shows that Muon’s Newton-Schulz algorithm aims at approximating the matrix sign function (as defined in (Higham, 2008)). In other words, if UkΣkVk⊤U_k _kV_k is a valid singular value decomposition of BkB_k, then Ok≈UkVk⊤=:sign(Bk)O_k≈ U_kV_k =: sign(B_k). From this observation, an idealized version of Muon (henceforth called ideal Muon) can be written as Hk H_k =βHk−1+(1−β)∇f(xk−1) =β H_k-1+(1-β)∇\!f(x_k-1) (5a) xk x_k =xk−1−αsign(Hk). =x_k-1-α sign(H_k). (5b) Furthermore, the ideal Muon is equivalent to the canonical Muon in the following sense. Proposition 1. Let (Bk)=sign(Bk)N(B_k)= sign(B_k). Also, for any x0∈ℝm×nx_0 ^m× n let x~k,B~k\ x_k, B_k\ and x¯k,H¯k\ x_k, H_k\ be the sequences generated by (4) and (5) respectively, each initialized at x0,B0\x_0,B_0\ and x0,H0\x_0,H_0\ with B0=(1−β)−1H0B_0=(1-β)^-1H_0. Then, for every k∈ℕk it holds that x~k=x¯k x_k= x_k and B~k=(1−β)−1H¯k B_k=(1-β)^-1 H_k. The proposition shows that the two different formulations induce identical trajectories in the parameter space, and up to a rescaling of the momentum. The (1−β)(1-β) factor is introduced to match the standard heavy-ball form, in which the gradient appears with weight (1−β)(1-β). This scaling does not affect Proposition 1, as it only rescales the momentum variable while leaving the induced updates unchanged. Crucially, it places (5) in the form of a preconditioned Polyak heavy-ball iteration, which will be exploited in the next section. Note that Ideal Muon does not converge asymptotically under constant α and β, since ‖sign(Hk)‖F=rank(Hk)\| sign(H_k)\|_F= rank(H_k) prevents vanishing step sizes near critical points. Existing analyses therefore establish ε -stationarity or averaged-gradient guarantees rather than pointwise convergence (Shen et al., 2025; Li and Hong, 2025). To address the issue of absence of vanishing step size near critical points, we introduce the soft-sign normalization. (x,ε):=x2+ε2S(x, ):= x x^2+ ^2 (6) for any ε>0 >0. The matrix-valued soft sign :ℝm×n→ℝm×nS:R^m× n ^m× n is, then, simply the scalar soft sign applied to the singular values of the input matrix. That is, let x=UΣV⊤x=U V be a valid SVD of a matrix x∈ℝm×nx ^m× n, then (x,ε):=U(Σ,ε)V⊤=(x⊤+ε2I)−1/2x,S(x, ):=US( , )V =(x + ^2I)^-1/2x, (7) which is well-defined for all x∈ℝm×nx ^m× n. This allows the formulation of the following iteration (henceforth called Soft Muon) Hk H_k =βHk−1+(1−β)∇f(xk−1) =β H_k-1+(1-β)∇\!f(x_k-1) (8a) xk x_k =xk−1−α(HkHk⊤+ε2I)−1/2Hk, =x_k-1-α(H_kH_k + ^2I)^-1/2H_k, (8b) which can also be shown to be equivalent to Muon in the following sense. Proposition 2. Let (Bk)=(Bk,ε)N(B_k)=S(B_k, ) for some ε>0 >0. Also, for any x0∈ℝm×nx_0 ^m× n let x~k,B~k\ x_k, B_k\ and x¯k,H¯k\ x_k, H_k\ be the sequences generated by (4) and (8) respectively, with each initialized at x0,B0\x_0,B_0\ and x0,H0\x_0,H_0\ with B0=(1−β)−1H0B_0=(1-β)^-1H_0. Further assume (8) generates its sequence with constant ε(1−β) (1-β) instead of ε . Then, for every k∈ℕk it holds that x~k=x¯k x_k= x_k and B~k=(1−β)−1H¯k B_k=(1-β)^-1 H_k. Note that neither Proposition 1 nor 2 assert that the Newton-Schulz polynomial is exactly equal to either the matrix sign or the soft sign. Rather, they show that once a particular surrogate for the Newton-Schulz map is chosen, the resulting canonical formulation can be rewritten exactly in the corresponding heavy-ball form. Unlike ideal Muon, the soft-sign formulation permits vanishing updates near critical points, enabling asymptotic convergence analysis, as we will show later in the paper. 2.2 Why so many Muons? (a) (b) (c) (d) Figure 1: Simulations for f(x)=x2f(x)=x^2 with α=0.1α=0.1 and β=0.5β=0.5. Columns (a)-(c) show the phase portrait (xk×Hkx_k× H_k) for Canonical Muon (4), Ideal Muon (5), and Soft Muon (8) respectively; (d) shows the loss across iterations. The first row uses ε=10−7 =10^-7 and δ=kδε10−7=484.87610−7δ=k_δ 10^-7=484.87610^-7, where δ denotes the regularizer inside the Newton-Schulz algorithm. To generate the second row, the values of δ and ε were increased until Canonical Muon and Soft Muon respectively started presenting a converging trajectory. Ideal Muon has no such regularizer and instead always settles into a periodic orbit around x∗=0x^*=0. Propositions 1 and 2 establish conditions in which Canonical Muon (4), Ideal Muon (5), and Soft Muon (8) generate identical parameter trajectories up to a rescaling of the momentum variable. Figure 1 illustrates how these formulations begin to differ once regularization becomes non-negligible. The first row corresponds to a weakly regularized regime, in which all three methods behave nearly identically and fail to converge asymptotically. In contrast, the second row shows that increasing δ and ε sufficiently induces vanishing step sizes near the optimum: Canonical and Soft Muon now converge toward x⋆=0x =0, while Ideal Muon continues to orbit the minimizer. Thus, the three formulations differ only through the effect of regularization – when the regularization is negligible they coincide, while sufficiently strong regularization restores asymptotic convergence. Although Canonical and Soft Muon exhibit the same qualitative transition, their trajectories may differ quantitatively for small values of δ and ε , reflecting the higher-order structure of the Newton-Schulz polynomial. Throughout the paper, we align the local behavior of the two formulations near the origin by matching their scalar first-order expansions, yielding the relation δ≈484.876ε=:kδεεδ≈ 484.876\, =:k_δ \, used in the simulations. Finally, both (5) and (8) can be interpreted as preconditioned variants of the Polyak heavy-ball method. The next section develops a convergence theory for preconditioned heavy-ball iterations and specializes it to Soft Muon. 3 Muon and the preconditioned Polyak heavy-ball The soft-sign formulation exposes Muon as a preconditioned heavy-ball method, allowing classical convergence analysis through the induced regularized preconditioner. The generic preconditioned heavy-ball iteration can be written as Hk H_k =βHk−1+(1−β)∇f(xk−1), =β H_k-1+(1-β)∇\!f(x_k-1), (9a) xk x_k =xk−1−αPk−1Hk, =x_k-1-α P_k^-1H_k, (9b) where Pk∈ℝm×mP_k ^m× m is a positive definite preconditioning matrix, and Pk−1HkP_k^-1H_k defines the descent direction. For this formulation, the following result can be stated. Theorem 1. Consider the preconditioned heavy-ball iteration presented in (9) for solving (3), and assume the preconditioning matrix PkP_k is uniformly bounded – i.e. there exist pmax>pmin>0p_ >p_ >0 such that pminI⪯Pk⪯pmaxIp_ I P_k p_ I for all k∈ℕk . Then there exist a constant β⋆∈(0,1)β ∈(0,1) and a positive function α⋆:[0,β⋆)→ℝ>0α :[0,β ) _>0 such that for every β∈[0,β⋆)β∈[0,β ) and every α∈(0,α⋆(β))α∈(0,α (β)), it holds that limk→∞‖∇f(xk)‖=0, _k→∞ \|∇\!f(x_k) \|=0, (10) that is, the method converges to a first-order stationary point in the classical sense. Furthermore, linear convergence follows naturally from this and the PŁ condition, as shown next. Corollary 1. Let f satisfy a global Polyak–Łojasiewicz inequality with constant μ>0μ>0, and let PkP_k be uniformly bounded. Then, for every β∈[0,β⋆)β∈[0,β ) and every α∈(0,α⋆(β))α∈(0,α (β)) (where β⋆β and α⋆(⋅)α (·) are as in Theorem 1), there exist constants C>0C>0 and ρ∈(0,1)ρ∈(0,1) such that f(xk)−f⋆≤Cρk,∀k∈ℕ.f(x_k)-f ≤ Cρ^k, ∀ k . (11) Theorem 1 and Corollary 1 provide a general convergence framework for preconditioned heavy-ball methods. In particular, they show that classical convergence guarantees hold whenever the preconditioning matrices remain uniformly bounded above and below. In the context of Muon, this condition is not immediate. While the regularization in the soft-sign formulation induces a natural lower bound on the preconditioner, a global upper bound does not follow directly from the assumptions on f. Establishing such a bound therefore becomes the key technical step required to apply the general theory. We will next show that the Soft Muon iteration (8) generates bounded trajectories under suitable conditions, which in turn implies uniform boundedness of the associated preconditioner. This allows us to invoke Theorem 1 and conclude convergence of the method. 3.1 Why soft Muon works The convergence results for preconditioned heavy-ball methods clarify why the soft-sign formulation is the right proxy for Muon’s Newton-Schulz iteration from an analysis perspective: the regularization parameter ε provides a uniform lower bound on the induced preconditioner. The lower bound alone is not sufficient to apply Theorem 1, since a uniform upper bound on PkP_k would also require a uniform bound on HkH_k. Such a bound is not available globally, because Pk=(HkHk⊤+ε2I)1/2P_k=(H_kH_k + ^2I)^1/2 grows with the momentum variable. The main purpose of this subsection is therefore to prove boundedness of the generated trajectory, from which the required upper bound follows. For that purpose, a key property of dε(Hk):=(HkHk⊤+ε2I)−1/2Hkd_ (H_k):=(H_kH_k + ^2I)^-1/2H_k is that ‖dε(H)‖<m\|d_ (H)\|< m for every H∈ℝm×nH ^m× n, so ‖xk−xk−1‖=α‖dε(Hk)‖<mα\|x_k-x_k-1\|=α\|d_ (H_k)\|< mα (12) holds unconditionally for any β∈[0,1)β∈[0,1). This bound is independent of both ‖Hk‖\|H_k\| and ‖∇f(xk)‖\|∇ f(x_k)\|, and is the mechanism that prevents arbitrarily large one-step excursions, allowing us to prove that for any initialization x0∈ℝm×nx_0 ^m× n, the sequence generated by (8) is bounded, circumventing the upper-bound issue. Theorem 2. Let xk,Hk\x_k,H_k\ be the sequence generated by (8) from the initial condition x0,∇f(x0)\x_0,∇\!f(x_0)\. Assume f is proper, L-smooth, has a unique global minimizer x⋆x , and satisfies a global Polyak–Łojasiewicz inequality with constant μ>0μ>0. Then there exists a positive function α⋆:[0,1)→ℝ>0α :[0,1) _>0 such that for every β∈[0,1)β∈[0,1) and for every α∈(0,α⋆(β))α∈(0,α (β)) the sequences xk\x_k\ and Hk\H_k\ are bounded. Once boundedness of the iterates is established, the preconditioner Pk=(HkHk⊤+ε2I)1/2P_k=(H_kH_k + ^2I)^1/2 is uniformly bounded along any trajectory. Hence the hypothesis of Theorem 1 hold, yielding following convergence result. Corollary 2. Consider the Soft Muon (8), then there exists α¯(β)>0 α(β)>0 such that for every β∈[0,β⋆)β∈[0,β ) and for every α∈(0,α¯(β))α∈(0, α(β)) the sequences xk\x_k\ and Hk\H_k\ are bounded and limk→∞‖∇f(xk)‖=0 _k→∞\|∇\!f(x_k)\|=0, where β⋆=1/2β =1/ 2 is the threshold from Theorem 1. Corollary 1 establishes linear convergence of Soft Muon under a global PŁ inequality, provided the learning rate and momentum parameter satisfy the bounds of Theorem 1. Together with Theorem 2, this shows that the regularization induced by the soft-sign formulation restores classical convergence guarantees despite the normalization-based update structure. Note that the constraint β<1/2β<1/ 2 does not reflect the empirical use of Muon, where β is typically picked at around 0.950.95. The constraint comes from the general result in Theorem 1, where convergence needs to be guaranteed for any bounded preconditioner, and does not mean that β>1/2β>1/ 2 will necessarily break convergence for Muon specifically. We believe that the specific structure of the Newton-Schulz iteration can be leveraged to tighten this constraint, and will investigate this gap in future works. Having established this connection between normalization-based updates and preconditioned heavy-ball dynamics, we next investigate whether the same perspective can be extended to Nesterov-type acceleration schemes. 4 A Nesterov-based Muon-type algorithm Nesterov acceleration naturally motivates investigating whether the normalization structure identified for Muon can be extended beyond heavy-ball dynamics, and whether similar convergence guarantees can be established. In light of this question, consider the preconditioned Nesterov algorithm below. yk y_k =xk+β(xk−xk−1), =x_k+β(x_k-x_k-1), (13a) xk+1 x_k+1 =yk−αPk−1∇f(yk), =y_k-α P_k^-1∇\!f(y_k), (13b) where α>0α>0, β∈[0,1)β∈[0,1) are the learning rate and momentum parameter, respectively, yky_k is the momentum-extrapolated point at which the gradient is evaluated, andPk−1∇f(yk)P_k^-1∇ f(y_k) defines the preconditioned descent direction. Analogously to the heavy-ball case, we establish convergence of the preconditioned Nesterov method under uniform boundedness of the preconditioner. Theorem 3. Consider the preconditioned Nesterov defined in (13) for solving (3), and assume the preconditioning matrix PkP_k is uniformly bounded – i.e. there exist pmax>pmin>0p_ >p_ >0 such that pminI⪯Pk⪯pmaxIp_ I P_k p_ I for all k∈ℕk . Then there exist a constant β∗∈(0,1)β^*∈(0,1) and a positive function α∗:[0,β∗)→ℝ>0α^*:[0,β^*) _>0 such that for every β∈[0,β∗)β∈[0,β^*) and every α∈(0,α∗(β))α∈(0,α^*(β)), it holds that limk→∞‖∇f(yk)‖=0. _k→∞\|∇ f(y_k)\|=0. (14) That is, the gradient evaluated at the extrapolated points converges to zero. Under a global Polyak–Łojasiewicz condition, this again implies linear convergence of the function values. Corollary 3. Let f satisfy a global Polyak–Łojasiewicz inequality with constant μ>0μ>0. Then, for every β∈[0,β⋆)β∈[0,β ) and every α∈(0,α⋆(β))α∈(0,α (β)) (where β∗β^* and α∗(⋅)α^*(·) are as in Theorem 3), there exist constants C>0C>0 and ρ∈(0,1)ρ∈(0,1) such that f(xk)−f⋆≤Cρk,∀k∈ℕ.f(x_k)-f ≤ Cρ^k, ∀ k . (15) Theorem 3 and Corollary 3 extend the convergence guarantees of Section 3 to the Nesterov setting. In particular, they show that preconditioned Nesterov methods exhibit the same qualitative behavior as their heavy-ball counterparts, provided the preconditioner remains uniformly bounded. We now specialize this framework to the preconditioning structure induced by Muon, leading to a Nesterov-based variant that we refer to as Muesterov. 4.1 The Muesterov algorithm Replacing the heavy-ball momentum in Soft Muon with a Nesterov extrapolation step leads to the following iteration, which we refer to as Muesterov. yk y_k =xk+β(xk−xk−1), =x_k+β(x_k-x_k-1), (16a) xk+1 x_k+1 =yk−α(GkGk⊤+ε2I)−1/2Gk, =y_k-α(G_kG_k + ^2I)^-1/2G_k, (16b) where Gk:=∇f(yk)G_k:=∇ f(y_k). Note that this preconditioned Nesterov formulation uses an analogous preconditioner to Muon, since (GkGk⊤+ε2I)−1/2Gk→sign(Gk)(G_kG_k + ^2I)^-1/2G_k→ sign(G_k) as ε→0 → 0. Moreover, the update can be implemented using the same Newton-Schulz iteration as in Muon, allowing the method to inherit the same computational advantages and parallelization structure. As in the heavy-ball case, the key step is to establish boundedness of the generated trajectory, which ensures uniform boundedness of the induced preconditioner. We provide this result in the following theorem. Theorem 4. Let xk,yk\x_k,y_k\ be the sequence generated by (16) from the initial condition x0,x0\x_0,x_0\. Assume f is proper, L-smooth, has a unique global minimizer x∗x^*, and satisfies a global PŁ inequality with constant μ>0μ>0. Then there exist β∗∈(0,1]β^*∈(0,1] and a positive function α∗:[0,β∗)→ℝ>0α^*:[0,β^*) _>0 such that for every β∈[0,β∗)β∈[0,β^*) and every α∈(0,α∗(β))α∈(0,α^*(β)) the sequences xk\x_k\ and yk\y_k\ remain bounded. Once boundedness is established, the preconditioner remains uniformly bounded along the trajectory, allowing us to invoke Theorem 3 to obtain convergence. Corollary 4. For the algorithm presented in (16), there exists β¯∈(0,1) β∈(0,1) and α¯:[0,β¯)→ℝ+ α:[0, β) _+, such that every β∈[0,β¯)β∈[0, β) and every α∈(0,α¯(β))α∈(0, α(β)) the sequences xk\x_k\ and yk\y_k\ are bounded and limk→∞‖∇f(yk)‖=0 _k→∞\|∇\!f(y_k)\|=0, i.e. the gradient evaluated at the extrapolated points converges to zero. We have thus extended the convergence framework developed for Muon to a Nesterov-based setting, showing that normalization-based preconditioning is compatible with accelerated momentum schemes. Crucially, Muesterov preserves the normalization-based structure of Muon, differing only in the point at which the gradient is evaluated. This provides both a theoretical justification for the proposed algorithm and a principled foundation for comparing Muon and Muesterov. 5 Simulations 5.1 Building intuition through a synthetic example (a) (b) (c) (d) Figure 2: Synthetic simulations for Canonical Muon (4) (first row, gray), Soft Muon (8) (second row, blue), and Muesterov (16) (third row, red). Column (a) depicts four training cost trajectories for α=10−3α=10^-3, β=0.5β=0.5, and increasing values of ε . Column (b) depicts the training cost for β=0.5β=0.5, ε=10−7 =10^-7, and increasing values of α. Columns (c) and (d) repeat the same simulations as column (b) but for much larger values of ε for the larger α. To build intuition for the role of the hyperparameters α and ε , we consider a simple scalar cross-entropy minimization problem for which closed-form expressions are available. This setting isolates the effect of normalization and regularization, and provides a transparent testbed for understanding the dynamics of Muon, Soft Muon, and Muesterov. Specifically, let the true label follow a Bernoulli distribution BpB_p with p=0.7p=0.7, and consider the parametric family qθ(X=1)=θq_θ(X=1)=θ, qθ(X=0)=1−θq_θ(X=0)=1-θ. The associated cross-entropy loss and its derivative are given by f(θ)=−plog(θ)−(1−p)log(1−θ)f(θ)=-p (θ)-(1-p) (1-θ), and f′(θ)=−pθ+1−p1−θf (θ)=- pθ+ 1-p1-θ, with unique minimizer θ∗=pθ^*=p. Furthermore, f satisfies a global PŁ inequality on (0,1)(0,1) with constant μ=2(p1/3+(1−p)1/3)3μ=2 (p^1/3+(1-p)^1/3 )^3. In this scalar setting, the soft-sign operator makes the effect of ε particularly transparent, (h,ε)=h/ε+O(h3)near h=0S(h, )=h/ +O(h^3)~near h=0, while (h,ε)≈sign(h)when |h|≫εS(h, )≈ sign(h)~when |h| . Figure 2 compares Canonical Muon, Soft Muon, and Muesterov under different choices of α and ε . For small regularization, all methods behave similarly and exhibit persistent oscillations characteristic of nearly sign-based updates. Increasing ε transitions Soft Muon and Muesterov toward a gradient-like regime with effective gain proportional to 1/ε1/ , stabilizing the dynamics and restoring convergence. Canonical Muon exhibits the same qualitative transition: for small δ, its trajectories are nearly indistinguishable from Soft Muon, while larger δ introduces quantitative differences due to the higher-order structure of the Newton-Schulz polynomial. The experiments further suggest that α and the normalization regularizer must be tuned jointly. Large learning rates combined with weak regularization produce sustained oscillations, whereas increasing ε (or δ) restores stable behavior and can substantially accelerate convergence relative to the nearly sign-based regime. Excessively large regularization, however, eventually slows convergence. Experiments also indicate potential advantage of using Muesterov over Muon, with the former reaching lower loss values faster than both Canonical and Soft Muon under identical hyperparameters in a simple 1D loss simulation. We next examine the problem of training GPT models, one of the primary practical applications of Muon. 5.2 Training nanoGPT Figure 3: Fixed-batch nanoGPT overfitting at 100×100× the baseline script-level learning rate (α=0.36α=0.36 instead of 0.00360.0036). The Muon optimizer trains the transformer blocks at 0.1α0.1α while AdamW remains at the baseline rate (0.00360.0036) on the output head. After the first minibatch is drawn it is reused at every step, so the plotted train loss is a fixed-batch diagnostic rather than a sampled training curve. Faint traces show the per-step values; bold traces show a 5050-step moving average. Top Panel. With Newton-Schulz regularizer δ=0.1δ=0.1, the trajectory exhibits sustained peak–trough oscillations superimposed on an overall descent. Bottom Panel. With δ=1δ=1, the same learning rate produces a single early overshoot followed by smooth, near-monotone descent toward a fixed point. The contrast is consistent with the soft-sign / Newton-Schulz analysis: larger δ reduces the effective gain of the normalized update and damps the oscillatory regime. The preceding experiments isolate, in a low-dimensional setting, how the learning rate α and the regularization in the soft-sign/Newton-Schulz normalization interact to produce either oscillatory behavior or convergence. We next investigate whether the same qualitative mechanism is visible in a practical transformer-training setting. Rather than performing a standard large-scale training benchmark, the goal of this experiment is to evaluate if intuition from theory is verifiable using a real model: we intentionally construct a controlled setting in which oscillatory behavior becomes visible, allowing us to directly examine the role of the Newton-Schulz regularizer in the practical Muon implementation. To achieve this goal, we edited the original script from (Keller and contributors, 2025) in the following ways: we increased the script-level learning rate from 0.00360.0036 to 0.360.36, with the AdamW learning rate for the output-head parameters kept at its baseline value, of 0.00360.0036; we also used a short 5050-iteration learning rate warmup, removed the warmdown, and fixed all seeds. Crucially, we also implemented the optimizer to be fixed to only one batch: after the initial batch is drawn, the same batch is reused throughout training. The plotted training loss should therefore be read as a fixed-batch loss, rather than as a standard minibatch-sampled training curve. Figure 3 shows that the smaller-regularization run develops visible oscillatory behavior. Although the loss decreases overall, the trajectory is repeatedly interrupted by sharp spikes followed by recovery. Since the same batch is reused throughout the run, these oscillations are not caused by changing minibatches. Rather, they are consistent with the regularization mechanism identified in the synthetic experiments and convergence analysis. Increasing the regularizer from δ=0.1δ=0.1 to δ=1δ=1 substantially smooths the trajectory. The run still has an initial transient, but the repeated large-amplitude peak-trough events are strongly reduced. This agrees with the soft-sign interpretation of the Newton-Schulz regularizer. A larger regularizer reduces the effective gain of the normalized update when the momentum matrix is small relative to the normalization scale. Thus, while this fixed-batch nanoGPT experiment is intentionally diagnostic rather than benchmark-oriented, it nevertheless exhibits the same qualitative α-regularization interaction predicted by the synthetic experiments and theoretical analysis. 6 Conclusion In this paper we established the first asymptotic convergence guarantees for Muon by interpreting the regularization implicit in its Newton-Schulz implementation as a bounded preconditioning mechanism. This perspective exposes Muon as a preconditioned Polyak heavy-ball method and enables a classical Lyapunov analysis under suitable hyperparameter choices. Motivated by the same interpretation, we introduced Muesterov, a Nesterov-based variant that preserves the normalization structure of Muon while enjoying analogous convergence guarantees. Beyond the theoretical results, the synthetic and transformer-training experiments suggest that the interaction between the learning rate and the Newton-Schulz regularizer plays a central role in controlling the transition between oscillatory sign-like dynamics and stable convergence. Understanding how these regularization mechanisms behave in practical large-scale training regimes remains an important direction for future work. References Amsel et al. (2025) N. Amsel, D. Persson, C. Musco, and R. M. Gower The polar express: optimal matrix sign methods and their application to the muon algorithm. arXiv preprint arXiv:2505.16932. Cited by: §1, §2.1. Anonymous (2026) Anonymous Code for “convergence guarantees for muon”. Note: https://anonymous.4open.science/r/2026_why_muon_works-217E/Anonymous repository Cited by: §1. Bernstein and Newhouse (2024) J. Bernstein and L. Newhouse Modular duality in deep learning. arXiv. External Links: Document, Link Cited by: §1, §2.1. Bernstein et al. (2018) J. Bernstein, Y. Wang, K. Azizzadenesheli, and A. Anandkumar SignSGD: compressed optimisation for non-convex problems. In International Conference on Machine Learning, p. 560–569. Cited by: §1. Chang et al. (2025) D. Chang, Y. Liu, and G. Yuan On the convergence of muon and beyond. arXiv preprint arXiv:2509.15816. Cited by: §1, §1. Chen et al. (2023) X. Chen, C. Liang, D. Huang, E. Real, K. Wang, H. Pham, X. Dong, T. Luong, C. Hsieh, Y. Lu, et al. Symbolic discovery of optimization algorithms. Advances in neural information processing systems 36, p. 49205–49233. Cited by: §1. Gupta et al. (2018) V. Gupta, T. Koren, and Y. Singer Shampoo: preconditioned stochastic tensor optimization. In International Conference on Machine Learning, p. 1842–1850. Cited by: §1. Higham (2008) N. J. Higham Functions of matrices: theory and computation. Society for Industrial and Applied Mathematics. External Links: ISBN 9780898717778, Link, Document Cited by: §1, §2.1, §2. Jordan et al. (2024) K. Jordan, Y. Jin, V. Boza, J. You, F. Cesista, L. Newhouse, and J. Bernstein Muon: an optimizer for hidden layers in neural networks. External Links: Link Cited by: §1, §2.1, §2.1, §2.1. Keller and contributors (2025) J. Keller and contributors modded-nanogpt. Note: https://github.com/kellerjordan/modded-nanogptGitHub repository. Commit: “Muon improvements” Cited by: §5.2. Kim and Oh (2026) G. Y. Kim and M. Oh Convergence of muon with newton-schulz. arXiv preprint arXiv:2601.19156. Cited by: §1. Kingma and Ba (2015) D. P. Kingma and J. Ba Adam: a method for stochastic optimization. In International Conference on Learning Representations, Cited by: §1. Lessard et al. (2016) L. Lessard, B. Recht, and A. Packard Analysis and design of optimization algorithms via integral quadratic constraints. SIAM Journal on Optimization 26 (1), p. 57–95. Cited by: §1. Li and Hong (2025) J. Li and M. Hong A note on the convergence of muon. arXiv preprint arXiv:2502.02900. Cited by: §1, §1, §2.1. Liu et al. (2024) H. Liu, Z. Li, D. Hall, P. Liang, and T. Ma Sophia: a scalable stochastic second-order optimizer for language model pre-training. In International Conference on Learning Representations, Cited by: §1. Liu et al. (2025) J. Liu, J. Su, X. Yao, Z. Jiang, G. Lai, Y. Du, Y. Qin, W. Xu, E. Lu, J. Yan, et al. Muon is scalable for LLM training. arXiv preprint arXiv:2502.16982. Cited by: §1. Loshchilov and Hutter (2019) I. Loshchilov and F. Hutter Decoupled weight decay regularization. In International Conference on Learning Representations, Cited by: §1. Martens and Grosse (2015) J. Martens and R. Grosse Optimizing neural networks with Kronecker-factored approximate curvature. In International Conference on Machine Learning, p. 2408–2417. Cited by: §1. Sato et al. (2025) N. Sato, H. Naganuma, and H. Iiduka Convergence bound and critical batch size of muon optimizer. arXiv preprint arXiv:2507.01598. Cited by: §1, §1. Shen et al. (2025) W. Shen, R. Huang, M. Huang, C. Shen, and J. Zhang On the convergence analysis of muon. arXiv preprint arXiv:2505.23737. Cited by: §1, §1, §2.1. Shi et al. (2021) B. Shi, S. S. Du, M. I. Jordan, and W. J. Su Understanding the acceleration phenomenon via high-resolution differential equations. Mathematical Programming 195, p. 79–148. Cited by: §1. Vyas et al. (2025) N. Vyas, D. Morwani, R. Zhao, I. Shapira, D. Brandfonbrener, L. Janson, and S. Kakade SOAP: improving and stabilizing Shampoo using Adam. In International Conference on Learning Representations, Cited by: §1. Appendix A Proofs of theoretical claims A.1 Proof of Proposition 1 Statement. Let (Bk)=sign(Bk)N(B_k)= sign(B_k). Also, for any x0∈ℝm×nx_0 ^m× n let x~k,B~k\ x_k, B_k\ and x¯k,H¯k\ x_k, H_k\ be the sequences generated by (4) and (5) respectively, each initialized at x0,B0\x_0,B_0\ and x0,H0\x_0,H_0\ with B0=(1−β)−1H0B_0=(1-β)^-1H_0. Then, for every k∈ℕk it holds that x~k=x¯k x_k= x_k and B~k=(1−β)−1H¯k B_k=(1-β)^-1 H_k. Proof. The proof is simply algebraic. First note that for any given sequence xk⊂ℝm×n\x_k\ ^m× n and for B0=(1−β)−1H0B_0=(1-β)^-1H_0 it holds that Hk=βkH0+(1−β)∑i=0k−1βk−1−i∇f(xi)=(1−β)(βkB0+∑i=0k−1βk−1−i∇f(xi))=(1−β)Bk,H_k=β^kH_0+(1-β) _i=0^k-1β^k-1-i∇\!f(x_i)=(1-β) (β^kB_0+ _i=0^k-1β^k-1-i∇\!f(x_i) )=(1-β)B_k, and thus sign(Hk)=sign((1−β)Bk)=sign(Bk) sign(H_k)= sign((1-β)B_k)= sign(B_k). Since both x~k x_k and x¯k x_k are initialized at the same point and have the same update law, it implies that x¯k,H¯k=x~k,(1−β)B~k\ x_k, H_k\=\ x_k,(1-β) B_k\, proving the claim. ∎ A.2 Proof of Proposition 2 Statement. Let (Bk)=(Bk,ε)N(B_k)=S(B_k, ) for some ε>0 >0. Also, for any x0∈ℝm×nx_0 ^m× n let x~k,B~k\ x_k, B_k\ and x¯k,H¯k\ x_k, H_k\ be the sequences generated by (4) and (8) respectively, with each initialized at x0,B0\x_0,B_0\ and x0,H0\x_0,H_0\ with B0=(1−β)−1H0B_0=(1-β)^-1H_0 and with constant ε(1−β) (1-β) in (8). Then, for every k∈ℕk it holds that x~k=x¯k x_k= x_k and B~k=(1−β)−1H¯k B_k=(1-β)^-1 H_k. Proof. The proof is, again, simply algebraic and follows almost identically to the previous proposition. First note that for any given sequence xk⊂ℝm×n\x_k\ ^m× n and for B0=(1−β)−1H0B_0=(1-β)^-1H_0 it holds that Hk=βkH0+(1−β)∑i=0k−1βk−1−i∇f(xi)=(1−β)(βkB0+∑i=0k−1βk−1−i∇f(xi))=(1−β)Bk,H_k=β^kH_0+(1-β) _i=0^k-1β^k-1-i∇\!f(x_i)=(1-β) (β^kB_0+ _i=0^k-1β^k-1-i∇\!f(x_i) )=(1-β)B_k, and thus (Hk,ε(1−β))=((1−β)Bk,ε(1−β))=(Bk,ε)S(H_k, (1-β))=S((1-β)B_k, (1-β))=S(B_k, ). Since both x~k x_k and x¯k x_k are initialized at the same point and have the same update law, it implies that x¯k,H¯k=x~k,(1−β)B~k\ x_k, H_k\=\ x_k,(1-β) B_k\, proving the claim. ∎ A.3 Proof of Theorem 1 Statement. Consider the preconditioned heavy-ball iteration presented in (9), and assume the preconditioning matrix PkP_k is uniformly bounded – i.e. there exist pmax>pmin>0p_ >p_ >0 such that pminI⪯Pk⪯pmaxIp_ I P_k p_ I for all k∈ℕk . Then there exist a constant β⋆∈(0,1)β ∈(0,1) and a positive function α⋆:[0,β⋆)→ℝ>0α :[0,β ) _>0 such that for every β∈[0,β⋆)β∈[0,β ) and every α∈(0,α⋆(β))α∈(0,α (β)), it holds that limk→∞‖∇f(xk)‖=0. _k→∞ \|∇\!f(x_k) \|=0. Proof. Let f⋆:=infxf(x)f := _xf(x) and define ek:=Hk−∇f(xk)e_k:=H_k-∇\!f(x_k), which satisfies e0=0e_0=0 by initialization. Henceforth adopt the shorthands Φk:=Φ(xk,ek) _k:= (x_k,e_k), Gk:=∇f(xk)G_k:=∇\!f(x_k), and dk:=Pk−1Hkd_k:=P_k^-1H_k, and define the Lyapunov function Φ(xk,ek):=f(xk)−f⋆+A‖ek‖2, (x_k,e_k):=f(x_k)-f +A \|e_k \|^2, (17) with A:=1/(4L)A:=1/(4L). We prove convergence by showing that for sufficiently small α and β, Φk _k is a Lyapunov function of the discrete-time system (9), i.e. that Φk+1−Φk<0 _k+1- _k<0. We first look at f(xk+1)−f(xk)f(x_k+1)-f(x_k). By L-smoothness of f, f(xk+1)≤f(xk)−α⟨Gk,dk⟩+Lα22‖dk‖2.f(x_k+1)≤ f(x_k)-α G_k,\,d_k + Lα^22 \|d_k \|^2. Writing Hk=Gk+ekH_k=G_k+e_k and using (1/pmax)I⪯Pk−1⪯(1/pmin)I(1/p_ )I P_k^-1 (1/p_ )I, ⟨Gk,dk⟩=⟨Gk,Pk−1Gk⟩+⟨Gk,Pk−1ek⟩≥1pmax‖Gk‖2−1pmin‖Gk‖‖ek‖, G_k,\,d_k = G_k,\,P_k^-1G_k + G_k,\,P_k^-1e_k ≥ 1p_ \|G_k \|^2- 1p_ \|G_k \| \|e_k \|, and ‖dk‖≤(1/pmin)‖Hk‖≤(1/pmin)(‖Gk‖+‖ek‖) \|d_k \|≤(1/p_ ) \|H_k \|≤(1/p_ )( \|G_k \|+ \|e_k \|), so f(xk+1)−f(xk)≤−(αpmax−Lα2pmin2)‖Gk‖2+αpmin‖Gk‖‖ek‖+Lα2pmin2‖ek‖2.f(x_k+1)-f(x_k)≤-\! ( αp_ - Lα^2p_ ^2 ) \|G_k \|^2+ αp_ \|G_k \| \|e_k \|+ Lα^2p_ ^2 \|e_k \|^2. (18) Next, consider ‖ek+1‖2−‖ek‖2 \|e_k+1 \|^2- \|e_k \|^2. From the preconditioned heavy-ball recursion (9) and the definition of eke_k, ek+1=βek+(Gk−Gk+1),e_k+1=β e_k+ (G_k-G_k+1 ), and by L-smoothness, ‖Gk+1−Gk‖≤Lα‖dk‖≤(Lα/pmin)(‖Gk‖+‖ek‖) \|G_k+1-G_k \|≤ Lα \|d_k \|≤(Lα/p_ )( \|G_k \|+ \|e_k \|). Applying (a+b)2≤2a2+2b2(a+b)^2≤ 2a^2+2b^2, A(‖ek+1‖2−‖ek‖2)≤4AL2α2pmin2‖Gk‖2−A(1−2β2−4L2α2pmin2)‖ek‖2.A ( \|e_k+1 \|^2- \|e_k \|^2 )≤ 4AL^2α^2p_ ^2 \|G_k \|^2-A\! (1-2β^2- 4L^2α^2p_ ^2 ) \|e_k \|^2. (19) Writing the iteration for Φk _k by adding (18) and (19) results in Φk+1−Φk≤−(αpmax−2Lα2pmin2)⏟=:CG‖Gk‖2+αpmin⏟=:Cc‖Gk‖ek‖−(1−2β24L−2Lα2pmin2)⏟=:Ce‖ek‖2, _k+1- _k≤- ( αp_ - 2Lα^2p_ ^2 )_=:\,C_G \|G_k \|^2+ αp_ _=:\,C_c \|G_k \| \|e_k \|- ( 1-2β^24L- 2Lα^2p_ ^2 )_=:\,C_e \|e_k \|^2, (20) where we used A=1/(4L)A=1/(4L) to collect the 2Lα2/pmin22Lα^2/p_ ^2 factor appearing identically in both the ‖Gk‖2 \|G_k \|^2 and ‖ek‖2 \|e_k \|^2 terms. Now Φk+1−Φk≤0 _k+1- _k≤ 0 whenever the quadratic form CG‖Gk‖2+Ce‖ek‖2−Cc‖Gk‖‖ek‖C_G \|G_k \|^2+C_e \|e_k \|^2-C_c \|G_k \| \|e_k \| is positive semi-definite, which by Sylvester’s criterion holds if and only if CG>0C_G>0, Ce>0C_e>0, and 4CGCe>Cc24C_GC_e>C_c^2. For fixed β∈[0,1)β∈[0,1), the condition CG>0C_G>0 reads α<pmin22Lpmax=:α1,α< p_ ^22Lp_ =: _1, and the condition Ce>0C_e>0 reads α<pmin2L1−2β22=:α2(β),α< p_ 2L 1-2β^22=: _2(β), with no constraint for β=0β=0 and α2(β)>0 _2(β)>0 for all β<1/2β<1/ 2. Since 4CGCe−Cc24C_GC_e-C_c^2 vanishes at α=0α=0, we factor out α>0α>0 and the discriminant condition 4CGCe>Cc24C_GC_e>C_c^2 is equivalent to Q(α)>0Q(α)>0, where Q(α):=16L2pmin4α3−8Lpmaxpmin2α2+4β2−3pmin2α+1−2β2Lpmax.Q(α):= 16L^2p_ ^4\,α^3- 8Lp_ p_ ^2\,α^2+ 4β^2-3p_ ^2\,α+ 1-2β^2Lp_ . (21) The constant term in (21) is strictly positive if and only if β2<1/2β^2<1/2, i.e. β<β⋆β<β where β⋆:=12∈(0,1).β := 1 2∈(0,1). (22) Since the leading coefficient of (21) is positive and the constant term is positive for β<β⋆β<β , we have Q(0)>0Q(0)>0 and Q(α)→−∞Q(α)→-∞ as α→−∞α→-∞, so Q has at least one negative root. Moreover, Q(α1)<0Q( _1)<0 and Q(α)→+∞Q(α)→+∞ as α→+∞α→+∞, hence Q has two positive roots α−(β)≤α+(β) _-(β)≤ _+(β) and is positive on (0,α−(β))(0, _-(β)). To verify that α−(β)<α1 _-(β)< _1, note that at α=α1α= _1 we have CG=0C_G=0, so 4CGCe−Cc2=−Cc2<04C_GC_e-C_c^2=-C_c^2<0, which implies Q(α1)<0Q( _1)<0. Hence α1∈[α−(β),α+(β)] _1∈[ _-(β), _+(β)] and α−(β)<α1 _-(β)< _1. By the same argument, Q(α2(β))<0Q( _2(β))<0, so α−(β)<α2(β) _-(β)< _2(β). Setting α⋆(β):=α−(β)α (β):= _-(β), all three conditions CG>0C_G>0, Ce>0C_e>0, and 4CGCe>Cc24C_GC_e>C_c^2 hold for every α∈(0,α⋆(β))α∈(0,α (β)). It remains to convert positive definiteness of the quadratic form into a bound with independent squared terms. Applying Young’s inequality with δ=Ce/CG>0δ= C_e/C_G>0 to the cross term gives Cc‖Gk‖‖ek‖≤Cc2CGCe‖Gk‖2+Cc2CeCG‖ek‖2,C_c \|G_k \| \|e_k \|≤ C_c2 C_GC_e \|G_k \|^2+ C_c2 C_eC_G \|e_k \|^2, and therefore Φk+1−Φk≤−c1‖Gk‖2−c2‖ek‖2, _k+1- _k≤-c_1 \|G_k \|^2-c_2 \|e_k \|^2, (23) where c1:=CG−(Cc/2)CG/Ce>0c_1:=C_G-(C_c/2) C_G/C_e>0 and c2:=Ce−(Cc/2)Ce/CG>0c_2:=C_e-(C_c/2) C_e/C_G>0, both of which are strictly positive since 4CGCe>Cc24C_GC_e>C_c^2 implies 2CGCe>Cc2 C_GC_e>C_c. Since Φk≥0 _k≥ 0 by definition, summing over k gives ∑k=0∞‖Gk‖2≤Φ0c1<∞, _k=0^∞ \|G_k \|^2≤ _0c_1<∞, and in particular ‖Gk‖2→0 \|G_k \|^2→ 0, which proves the statement. ∎ A.4 Proof of Corollary 1 Statement. Let f satisfy a global Polyak–Łojasiewicz inequality with constant μ>0μ>0, and let PkP_k be uniformly bounded. Then, for every β∈[0,β⋆)β∈[0,β ) and every α∈(0,α⋆(β))α∈(0,α (β)) (where β⋆β and α⋆(⋅)α (·) are as in Theorem 1), there exist constants C>0C>0 and ρ∈(0,1)ρ∈(0,1) such that f(xk)−f⋆≤Cρk,∀k∈ℕ.f(x_k)-f ≤ Cρ^k, ∀ k . (24) Proof. From the proof of Theorem 1, for every β∈[0,β⋆)β∈[0,β ) and every α∈(0,α⋆(β))α∈(0,α (β)), there exist constants c1,c2>0c_1,c_2>0 such that (23) holds. We now upper-bound Φk _k in terms of ‖Gk‖2 \|G_k \|^2 and ‖ek‖2 \|e_k \|^2. By the Polyak–Łojasiewicz inequality applied at xkx_k, f(xk)−f⋆≤1μ‖Gk‖2,f(x_k)-f ≤ 1μ \|G_k \|^2, and recalling the definition (17) of Φk _k, Φk≤1μ‖Gk‖2+A‖ek‖2≤C~(c1‖Gk‖2+c2‖ek‖2), _k≤ 1μ \|G_k \|^2+A \|e_k \|^2≤ C (c_1 \|G_k \|^2+c_2 \|e_k \|^2 ), (25) where C~:=max1/(μc1),A/c2∈(0,∞) C:= \1/(μ c_1),\,A/c_2\∈(0,∞). Combining (23) and (25), Φk+1−Φk≤−(c1‖Gk‖2+c2‖ek‖2)≤−1C~Φk. _k+1- _k≤- (c_1 \|G_k \|^2+c_2 \|e_k \|^2 )≤- 1 C\, _k. Therefore, defining δ:=1/C~>0δ:=1/ C>0, we obtain Φk+1≤(1−δ)Φk _k+1≤(1-δ) _k. Since Φk≥0 _k≥ 0 for all k, we may take ρ∈(0,1)ρ∈(0,1) such that Φk+1≤ρΦk _k+1≤ρ _k. Iterating gives Φk≤ρkΦ0. _k≤ρ^k _0. Finally, since f(xk)−f⋆≤Φkf(x_k)-f ≤ _k, it follows that f(xk)−f⋆≤Φ0ρk.f(x_k)-f ≤ _0ρ^k. Setting C:=Φ0C:= _0 concludes the proof. ∎ A.5 Proof of Theorem 2 Statement. Let xk,Hk\x_k,H_k\ be the sequence generated by (8) from the initial condition x0,∇f(x0)\x_0,∇\!f(x_0)\. Assume f is proper, L-smooth, has a unique global minimizer x⋆x , and satisfies a global Polyak–Łojasiewicz inequality with constant μ>0μ>0. Then for every β∈[0,1)β∈[0,1) there exists a positive function α⋆:[0,1)→ℝ>0α :[0,1) _>0 such that for every α∈(0,α⋆(β))α∈(0,α (β)) the sequences xk\x_k\ and Hk\H_k\ are bounded. Proof. For every C≥f(x0)−f⋆C≥ f(x_0)-f , fix an open neighbourhood D of x⋆x with D¯⊂int(DC) D (D_C), where DC:=f≤f⋆+CD_C:=\f≤ f +C\, and define rD:=dist(D¯,DCc)>0r_D:=dist( D,D_C^c)>0, γ(C):=minx∈DC∖D‖∇f(x)‖>0,Γ(C):=maxx∈DC‖∇f(x)‖<∞,γ(C):= _x∈ D_C D\|∇\!f(x)\|>0, (C):= _x∈ D_C\|∇\!f(x)\|<∞, (26) and α⋆(β,C):=min2γ(C)2(1−β)L(2βΓ(C)m+m(1−β)Γ(C)2+ε2),rD(C)/m>0.α (β,C):= \! \ 2γ(C)^2(1-β)L (2β (C) m+m(1-β) (C)^2+ ^2 ),\;r_D(C)/ m \>0. (27) Since α⋆(β,C)>0α (β,C)>0 for every C≥f(x0)−f⋆C≥ f(x_0)-f , define α⋆(β):=supC≥f(x0)−f⋆α⋆(β,C)>0α (β):= _C≥ f(x_0)-f α (β,C)>0. For any α∈(0,α⋆(β))α∈(0,α (β)), there exists C≥f(x0)−f⋆C≥ f(x_0)-f such that α<α⋆(β,C)α<α (β,C). Fix such C and write γ:=γ(C)γ:=γ(C), Γ:=Γ(C) := (C) for brevity. We prove by induction that xk∈DCx_k∈ D_C and ‖Hk‖≤Γ\|H_k\|≤ for all k≥0k≥ 0. The base case holds since x0∈DCx_0∈ D_C by definition of C and ‖H0‖=‖∇f(x0)‖≤Γ\|H_0\|=\|∇\!f(x_0)\|≤ . Assume x0,…,xk∈DCx_0,…,x_k∈ D_C and ‖Hj‖≤Γ\|H_j\|≤ for all j≤kj≤ k. Since e0=H0−∇f(x0)=0e_0=H_0-∇\!f(x_0)=0 and the lag error ej:=Hj−∇f(xj)e_j:=H_j-∇\!f(x_j) satisfies the recursion ej+1=βej+(∇f(xj)−∇f(xj+1))e_j+1=β e_j+(∇\!f(x_j)-∇\!f(x_j+1)), iterating with ‖∇f(xj+1)−∇f(xj)‖≤Lmα\|∇\!f(x_j+1)-∇\!f(x_j)\|≤ L mα (from L-smoothness and (12)) gives ‖ek‖≤Lmα∑j=0k−1βj<mLα1−β,\|e_k\|≤ L mα _j=0^k-1β^j< mLα1-β, (28) uniformly in k. The bound on ‖Hk+1‖\|H_k+1\| then follows directly: ‖Hk+1‖≤β‖Hk‖+(1−β)‖∇f(xk)‖≤Γ.\|H_k+1\|≤β\|H_k\|+(1-β)\|∇\!f(x_k)\|≤ . It remains to show xk+1∈DCx_k+1∈ D_C. Since xk+1=xk−αdε(Hk+1)x_k+1=x_k-α d_ (H_k+1), we consider two cases. If xk∈Dx_k∈ D, then ‖xk+1−xk‖=α‖dε(Hk+1)‖<αm≤rD\|x_k+1-x_k\|=α\|d_ (H_k+1)\|<α m≤ r_D by (12) and the choice α<α⋆(β,C)≤rDα<α (β,C)≤ r_D, so xk+1x_k+1 cannot exit DCD_C. If xk∈DC∖Dx_k∈ D_C D, then ‖∇f(xk)‖≥γ\|∇\!f(x_k)\|≥γ and by L-smoothness, f(xk+1)≤f(xk)−α⟨∇f(xk),dε(Hk+1)⟩+Lα2m2.f(x_k+1)≤ f(x_k)-α ∇\!f(x_k),\,d_ (H_k+1) + Lα^2m2. The step dε(Hk+1)d_ (H_k+1) involves Hk+1=∇f(xk)+hkH_k+1=∇\!f(x_k)+h_k, where the prospective lag hk:=Hk+1−∇f(xk)=βekh_k:=H_k+1-∇\!f(x_k)=β e_k satisfies ‖hk‖=β‖ek‖<βmLα1−β\|h_k\|=β\|e_k\|< β mLα1-β by (28). Using ‖Hk+1‖≤Γ\|H_k+1\|≤ and Hk+1=∇f(xk)+hkH_k+1=∇\!f(x_k)+h_k, we have ⟨∇f(xk),dε(Hk+1)⟩ ∇\!f(x_k),\,d_ (H_k+1) =⟨∇f(xk),(Hk+1Hk+1⊤+ε2I)−1/2Hk+1⟩ = ∇\!f(x_k),\,(H_k+1H_k+1 + ^2I)^-1/2H_k+1 =⟨∇f(xk),(Hk+1Hk+1⊤+ε2I)−1/2∇f(xk)⟩ = ∇\!f(x_k),\,(H_k+1H_k+1 + ^2I)^-1/2∇\!f(x_k) +⟨∇f(xk),(Hk+1Hk+1⊤+ε2I)−1/2hk⟩ + ∇\!f(x_k),\,(H_k+1H_k+1 + ^2I)^-1/2h_k ≥‖∇f(xk)‖2‖Hk+1‖2+ε2−1ε‖∇f(xk)‖‖hk‖ ≥ \|∇\!f(x_k)\|^2 \|H_k+1\|^2+ ^2- 1 \|∇\!f(x_k)\|\|h_k\| ≥γ2Γ2+ε2−Γε‖hk‖ ≥ γ^2 ^2+ ^2- \|h_k\| >γ2Γ2+ε2−βΓmLαε(1−β). > γ^2 ^2+ ^2- β mLα (1-β). Substituting, f(xk+1)≤f(xk)−α(γ2Γ2+ε2−βΓmLαε(1−β))+Lα2m2=:f(xk)−cα,f(x_k+1)≤ f(x_k)-α ( γ^2 ^2+ ^2- β mLα (1-β) )+ Lα^2m2=:f(x_k)-c_α, and cα>0c_α>0 holds whenever α<2γ2(1−β)L(2βΓm+m(1−β)Γ2+ε2),α< 2γ^2(1-β)L (2β m+m(1-β) ^2+ ^2 ), which is precisely the first term in (27). Therefore f(xk+1)≤f(xk)≤f⋆+Cf(x_k+1)≤ f(x_k)≤ f +C, i.e. xk+1∈DCx_k+1∈ D_C. By induction, xk∈DCx_k∈ D_C and ‖Hk‖≤Γ\|H_k\|≤ for all k≥0k≥ 0. Since f is proper, the sublevelset DCD_C is compact, implying that both sequences are bounded. ∎ A.6 Proof of Corollary 2 Statement. For the algorithm (8), for every β∈[0,β⋆)β∈[0,β ) there exists α¯(β)>0 α(β)>0 such that for every α∈(0,α¯(β))α∈(0, α(β)) the sequences xk\x_k\ and Hk\H_k\ are bounded and limk→∞‖∇f(xk)‖=0 _k→∞\|∇\!f(x_k)\|=0, where β⋆=1/2β =1/ 2 is the threshold from Theorem 1. Proof. Fix any β∈[0,β⋆)β∈[0,β ) and define α¯(β):=minα1⋆(β),α2⋆(β)>0 α(β):= \ _1 (β), _2 (β)\>0, where α1⋆(β) _1 (β) is as in Theorem 1 and α2⋆(β) _2 (β) is the threshold from Theorem 2. For any α∈(0,α¯(β))α∈(0, α(β)), Theorem 2 guarantees that xk\x_k\ and Hk\H_k\ remain in a compact set, so ‖Hk‖≤Γ\|H_k\|≤ uniformly and the eigenvalues of Pk=(HkHk⊤+ε2I)1/2P_k=(H_kH_k + ^2I)^1/2 lie in [ε,Γ2+ε2][ , ^2+ ^2] for all k≥0k≥ 0. We then conclude that PkP_k is bounded with pmin=εp_ = and pmax=Γ2+ε2p_ = ^2+ ^2, and Theorem 1 gives limk→∞‖∇f(xk)‖=0 _k→∞\|∇\!f(x_k)\|=0. ∎ A.7 Proof of Theorem 3 Statement. Consider the preconditioned Nesterov defined in (13), and assume the preconditioning matrix PkP_k is uniformly bounded – i.e. there exist pmax>pmin>0p_ >p_ >0 such that pminI⪯Pk⪯pmaxIp_ I P_k p_ I for all k∈ℕk . Then there exist a constant β∗∈(0,1)β^*∈(0,1) and a positive function α∗:[0,β∗)→ℝ>0α^*:[0,β^*) _>0 such that for every β∈[0,β∗)β∈[0,β^*) and every α∈(0,α∗(β))α∈(0,α^*(β)), it holds that limk→∞‖∇f(yk)‖=0. _k→∞\|∇ f(y_k)\|=0. Proof. Let f⋆:=infx∈ℝm×nf(x)f := _x ^m× nf(x) and define the function Φ(xk):=f(xk)−f⋆+A‖xk−xk−1‖2. (x_k):=f(x_k)-f +A\|x_k-x_k-1\|^2. (29) We will prove convergence by showing that Φ(⋅) (·) is a Lyapunov function of the discrete time system (13). Henceforth consider the shorthand Φk:=Φ(xk) _k:= (x_k), Gk:=∇f(yk)G_k:=∇\!f(y_k), and dk:=Pk−1Gkd_k:=P_k^-1G_k. Define further sk:=xk−xk−1s_k:=x_k-x_k-1 and then rewrite (13) as xk+1=xk+βsk−αdk, x_k+1=x_k+β s_k-α d_k, (30a) sk+1=βsk−αdk. s_k+1=β s_k-α d_k. (30b) To show that Φk _k is a Lyapunov function, we need to show that it decreases at every iteration, i.e. we need to show that Φk+1−Φk<0 _k+1- _k<0. Consider Φk+1−Φk=f(xk+1)−f(xk)+A(‖sk+1‖2−‖sk‖2). _k+1- _k=f(x_k+1)-f(x_k)+A(\|s_k+1\|^2-\|s_k\|^2). We will first look at the term ‖sk+1‖2−‖sk‖2\|s_k+1\|^2-\|s_k\|^2. Note that ‖sk+1‖2 \|s_k+1\|^2 =‖βsk−αdk‖2 =\|β s_k-α d_k\|^2 ≤β2‖sk‖2+α2‖dk‖2+2αβ‖sk‖‖dk‖. ≤β^2\|s_k\|^2+α^2\|d_k\|^2+2αβ\|s_k\|\|d_k\|. Next, consider f(xk+1)−f(xk)f(x_k+1)-f(x_k), and note that from L-smoothness of f and uniform boundedness of PkP_k we have that f(xk+1) f(x_k+1) ≤f(yk)−α⟨Gk,dk⟩+Lα22‖dk‖2 ≤ f(y_k)-α G_k,d_k + Lα^22\|d_k\|^2 ≤f(yk)−(αpmin−α2L2)⏟cα‖dk‖2, ≤ f(y_k)- (α p_ - α^2L2 )_c_α\|d_k\|^2, and note that cα>0c_α>0 for all 0<α<2pmin/L0<α<2p_ /L. Next, and again by L-smoothness of f, using the fact that yk=xk+βsky_k=x_k+β s_k, and Young’s inequality we get that f(yk) f(y_k) ≤f(xk)+β⟨∇f(xk⟩,sk)+L2β2∥sk∥2. ≤ f(x_k)+β ∇\!f(x_k ,s_k)+ L2β^2\|s_k\|^2. Then, note that ∇f(xk)=Gk+ek∇\!f(x_k)=G_k+e_k where ek=∇f(xk)−∇f(yk)e_k=∇\!f(x_k)-∇\!f(y_k). Then by L-smoothness of f and uniform boundedness of PkP_k ⟨∇f(xk⟩,sk) ∇\!f(x_k ,s_k) =⟨Gk+ek,sk⟩ = G_k+e_k,s_k =⟨Gk,sk⟩+⟨ek,sk⟩ = G_k,s_k + e_k,s_k ≤pmax‖dk‖‖sk‖+Lβ‖sk‖2. ≤ p_ \|d_k\|\|s_k\|+Lβ\|s_k\|^2. Hence f(yk) f(y_k) ≤f(xk)+βpmax‖dk‖‖sk‖+(Lβ2+L2β2)‖sk‖2 ≤ f(x_k)+β p_ \|d_k\|\|s_k\|+ (Lβ^2+ L2β^2 )\|s_k\|^2 =f(xk)+βpmax‖dk‖‖sk‖+3L2β2‖sk‖2. =f(x_k)+β p_ \|d_k\|\|s_k\|+ 3L2β^2\|s_k\|^2. Putting the current equations together results in f(xk+1)≤f(xk)+3L2β2‖sk‖2+βpmax‖dk‖‖sk‖−cα‖dk‖2.f(x_k+1)≤ f(x_k)+ 3L2β^2\|s_k\|^2+β p_ \|d_k\|\|s_k\|-c_α\|d_k\|^2. Writing the iteration for Φk _k then results in Φk+1−Φk _k+1- _k ≤(3L2β2+A(β2−1))⏟=:−Cs(α,β)∥sk∥2−(cα−Aα2)⏟=:Cd(α)∥dk∥2+β(pmax+2Aα)⏟=:Cc(α,β)∥sk∥dk∥. ≤ ( 3L2β^2+A(β^2-1) )_=:\,-C_s(α,β)\|s_k\|^2- (c_α-Aα^2 )_=:\,C_d(α)\|d_k\|^2+ β(p_ +2Aα)_=:\,C_c(α,β)\|s_k\|\|d_k\|. Now choose A=pmin4α.A= p_ 4α. Then 2Aα=pmin/22Aα=p_ /2, so Cc(α,β)=β(pmax+pmin/2)C_c(α,β)=β(p_ +p_ /2) is independent of α, and Cd(α)=cα−Aα2=α(3pmin4−Lα2)>0for all α∈(0,3pmin2L)=:(0,α1).C_d(α)=c_α-Aα^2=α ( 3p_ 4- Lα2 )>0 all α∈ (0, 3p_ 2L )=:(0, _1). Hence Φk+1−Φk≤0 _k+1- _k≤ 0 whenever the quadratic form Cs‖sk‖2+Cd‖dk‖2−Cc‖sk‖‖dk‖C_s\|s_k\|^2+C_d\|d_k\|^2-C_c\|s_k\|\|d_k\| is positive definite, which by Sylvester’s criterion holds if and only if Cs>0C_s>0 and 4CsCd−Cc2>04C_sC_d-C_c^2>0. For fixed β∈[0,1)β∈[0,1), the condition Cs>0C_s>0 reads explicitly α<pmin(1−β2)6Lβ2=:α2(β),α< p_ (1-β^2)6Lβ^2=: _2(β), (31) with no constraint for β=0β=0. Substituting Cc=β(pmax+pmin/2)C_c=β(p_ +p_ /2), Cs=pmin(1−β2)/(4α)−3Lβ2/2C_s=p_ (1-β^2)/(4α)-3Lβ^2/2, and Cd=α(3pmin/4−Lα/2)C_d=α(3p_ /4-Lα/2) into 4CsCd−Cc2>04C_sC_d-C_c^2>0 and expanding, the discriminant condition is equivalent to the quadratic inequality in α 3L2β2α2−Lpmin(1+8β2)2α+3pmin2(1−β2)4−β2(pmax+pmin2)2>0.3L^2β^2\,α^2- Lp_ (1+8β^2)2\,α+ 3p_ ^2(1-β^2)4-β^2\! (p_ + p_ 2 )^\!2>0. (32) The constant term in (32) is strictly positive if and only if β2<3pmin2/43pmin2/4+(pmax+pmin/2)2,β^2< 3p_ ^2/43p_ ^2/4+ (p_ +p_ /2 )^2, i.e. β<β∗β<β^* where β∗:=3pmin/23pmin2/4+(pmax+pmin/2)2∈(0,1).β^*:= 3p_ /2 3p_ ^2/4+ (p_ +p_ /2 )^2\;∈\;(0,1). (33) Since the leading coefficient of (32) is positive, the linear coefficient is negative, and the constant term is positive for β<β∗β<β^*, the quadratic opens upward and is positive at α=0α=0, so it has two positive roots α−(β)≤α+(β) _-(β)≤ _+(β) and is positive on (0,α−(β))(0, _-(β)). Define Δ(β):=L2pmin2(1+8β2)24−3L2β2[3pmin2(1−β2)4−β2(pmax+pmin2)2], (β):= L^2p_ ^2(1+8β^2)^24-3L^2β^2\! [ 3p_ ^2(1-β^2)4-β^2\! (p_ + p_ 2 )^\!2 ], (34) so that α−(β)=Lpmin(1+8β2)/2−Δ(β)6L2β2, _-(β)= Lp_ (1+8β^2)/2- (β)6L^2β^2, (35) and note that if β=0β=0 then the problem reduces to a preconditioned gradient descent which can easily be shown satisfy this Theorem. To see that α−(β)<α1 _-(β)< _1, evaluate (32) at α=α1=3pmin/(2L)α= _1=3p_ /(2L): 3L2β2⋅9pmin24L2−Lpmin(1+8β2)2⋅3pmin2L+3pmin2(1−β2)4−β2(pmax+pmin2)2=−β2(pmax+pmin2)2≤0,3L^2β^2· 9p_ ^24L^2- Lp_ (1+8β^2)2· 3p_ 2L+ 3p_ ^2(1-β^2)4-β^2\! (p_ + p_ 2 )^\!2=-β^2\! (p_ + p_ 2 )^\!2≤ 0, so the quadratic is non-positive at α1 _1, hence α−(β)<α1 _-(β)< _1 for all β∈(0,β∗)β∈(0,β^*), and Cd(α)>0C_d(α)>0 is automatically satisfied on (0,α−(β))⊂(0,α1)(0, _-(β))⊂(0, _1). Moreover, α−(β)<α2(β) _-(β)< _2(β) always holds for β<β∗β<β^*, since at α=α2(β)α= _2(β) one has Cs=0C_s=0, so 4CsCd−Cc2=−Cc2<04C_sC_d-C_c^2=-C_c^2<0, meaning α2(β) _2(β) lies above the root α−(β) _-(β), and hence Cs>0C_s>0 is automatically satisfied on (0,α−(β))(0, _-(β)). Setting α∗(β):=α−(β)α^*(β):= _-(β), all three conditions Cs>0C_s>0, Cd>0C_d>0, and 4CsCd>Cc24C_sC_d>C_c^2 hold for every α∈(0,α∗(β))α∈(0,α^*(β)). It remains to convert positive definiteness of the quadratic form into a bound with independent squared terms. Applying Young’s inequality with δ=Cd/Cs>0δ= C_d/C_s>0 to the cross term gives Cc‖sk‖‖dk‖≤CcCs2Cd‖sk‖2+CcCd2Cs‖dk‖2,C_c\|s_k\|\|d_k\|≤ C_c C_s2 C_d\|s_k\|^2+ C_c C_d2 C_s\|d_k\|^2, and therefore Cs∥sk∥2+Cd∥dk∥2−Cc∥sk∥dk∥≥Cs(Cs−Cc2Cd)⏟=:c1∥sk∥2+Cd(Cd−Cc2Cs)⏟=:c2∥dk∥2, C_s\|s_k\|^2+C_d\|d_k\|^2-C_c\|s_k\|\|d_k\|≥ C_s\! ( C_s- C_c2 C_d )_=:\,c_1\|s_k\|^2+ C_d\! ( C_d- C_c2 C_s )_=:\,c_2\|d_k\|^2, where c1,c2>0c_1,c_2>0 since 4CsCd>Cc24C_sC_d>C_c^2 implies 2CsCd>Cc2 C_sC_d>C_c. Therefore Φk+1−Φk≤−c1‖sk‖2−c2‖dk‖2. _k+1- _k≤-c_1\|s_k\|^2-c_2\|d_k\|^2. Since Φk≥0 _k≥ 0 by definition, summing over k gives ∑k=0∞‖sk‖2<∞,∑k=0∞‖dk‖2<∞. _k=0^∞\|s_k\|^2<∞, _k=0^∞\|d_k\|^2<∞. In particular, limk→∞‖dk‖=0 _k→∞\|d_k\|=0. Finally, by uniform boundedness of PkP_k, ‖Gk‖=‖Pkdk‖≤pmax‖dk‖→0,\|G_k\|=\|P_kd_k\|≤ p_ \|d_k\|→ 0, which proves (14). ∎ A.8 Proof of Corollary 3 Statement. Let f satisfy a global Polyak–Łojasiewicz inequality with constant μ>0μ>0. Then, for every β∈[0,β⋆)β∈[0,β ) and every α∈(0,α⋆(β))α∈(0,α (β)) (where β∗β^* and α∗(⋅)α^*(·) are as in Theorem 3), there exist constants C>0C>0 and ρ∈(0,1)ρ∈(0,1) such that f(xk)−f⋆≤Cρk,∀k∈ℕ.f(x_k)-f ≤ Cρ^k, ∀ k . (36) Proof. From the proof of Theorem 3, for every β∈[0,β∗)β∈[0,β^*) and every α∈(0,α∗(β))α∈(0,α^*(β)), there exist constants c1,c2>0c_1,c_2>0 such that Φk+1−Φk≤−c1‖sk‖2−c2‖dk‖2. _k+1- _k≤-c_1\|s_k\|^2-c_2\|d_k\|^2. (37) We now upper bound Φk _k in terms of ‖sk‖2\|s_k\|^2 and ‖dk‖2\|d_k\|^2. By L-smoothness of f and the identity yk=xk+βsky_k=x_k+β s_k, we have f(xk) f(x_k) ≤f(yk)+⟨Gk,xk−yk⟩+L2‖xk−yk‖2 ≤ f(y_k)+ G_k,x_k-y_k + L2\|x_k-y_k\|^2 =f(yk)−β⟨Gk,sk⟩+L2β2‖sk‖2 =f(y_k)-β G_k,s_k + L2β^2\|s_k\|^2 ≤f(yk)+β‖Gk‖‖sk‖+L2β2‖sk‖2. ≤ f(y_k)+β\|G_k\|\|s_k\|+ L2β^2\|s_k\|^2. Applying the PŁ inequality at x=ykx=y_k gives f(yk)−f⋆≤1μ‖Gk‖2.f(y_k)-f ≤ 1μ\|G_k\|^2. Therefore, f(xk)−f⋆≤1μ‖Gk‖2+β‖Gk‖‖sk‖+L2β2‖sk‖2.f(x_k)-f ≤ 1μ\|G_k\|^2+β\|G_k\|\|s_k\|+ L2β^2\|s_k\|^2. Next, by uniform boundedness of PkP_k we have ‖Gk‖=‖Pkdk‖≤pmax‖dk‖\|G_k\|=\|P_kd_k\|≤ p_ \|d_k\|, hence f(xk)−f⋆≤pmax2μ‖dk‖2+βpmax‖dk‖‖sk‖+L2β2‖sk‖2.f(x_k)-f ≤ p_ ^2μ\|d_k\|^2+β p_ \|d_k\|\|s_k\|+ L2β^2\|s_k\|^2. Applying Young’s inequality with parameter η>0η>0 to the cross term, βpmax‖dk‖‖sk‖≤η2‖dk‖2+β2pmax22η‖sk‖2,β p_ \|d_k\|\|s_k\|≤ η2\|d_k\|^2+ β^2p_ ^22η\|s_k\|^2, gives f(xk)−f⋆≤(pmax2μ+η2)‖dk‖2+(β2pmax22η+L2β2)‖sk‖2.f(x_k)-f ≤ ( p_ ^2μ+ η2 )\|d_k\|^2+ ( β^2p_ ^22η+ L2β^2 )\|s_k\|^2. Recalling the definition of Φk _k, we conclude that Φk≤(pmax2μ+η2)⏟=:Cd∥dk∥2+(A+β2pmax22η+L2β2)⏟=:Cs∥sk∥2. _k≤ ( p_ ^2μ+ η2 )_=:\,C_d\|d_k\|^2+ (A+ β^2p_ ^22η+ L2β^2 )_=:\,C_s\|s_k\|^2. (38) Combining (38) with (37) gives Φk+1−Φk _k+1- _k ≤−c1‖sk‖2−c2‖dk‖2 ≤-c_1\|s_k\|^2-c_2\|d_k\|^2 ≤−minc2Cd,c1Cs(Cd‖dk‖2+Cs‖sk‖2) ≤- \! \ c_2C_d, c_1C_s \\! (C_d\|d_k\|^2+C_s\|s_k\|^2 ) ≤−minc2Cd,c1CsΦk. ≤- \! \ c_2C_d, c_1C_s \ _k. Therefore, defining δ:=minc2/Cd,c1/Cs>0δ:= \c_2/C_d,\,c_1/C_s\>0, we obtain Φk+1≤(1−δ)Φk _k+1≤(1-δ) _k. Since Φk≥0 _k≥ 0 for all k, we may take ρ∈(0,1)ρ∈(0,1) such that Φk+1≤ρΦk _k+1≤ρ _k. Iterating gives Φk≤ρkΦ0. _k≤ρ^k _0. Finally, since f(xk)−f⋆≤Φkf(x_k)-f ≤ _k, it follows that f(xk)−f⋆≤Φ0ρk.f(x_k)-f ≤ _0ρ^k. Setting C:=Φ0C:= _0 concludes the proof. ∎ A.9 Proof of Theorem 4 Statement. Let xk,yk\x_k,y_k\ be the sequence generated by (16) from the initial condition x0,x0\x_0,x_0\. Assume f is proper, L-smooth, has a unique global minimiser x∗x^*, and satisfies a global PŁ inequality with constant μ>0μ>0. Then there exist β∗∈(0,1]β^*∈(0,1] and a positive function α∗:[0,β∗)→ℝ>0α^*:[0,β^*) _>0 such that for every β∈[0,β∗)β∈[0,β^*) and every α∈(0,α∗(β))α∈(0,α^*(β)) the sequences xk\x_k\ and yk\y_k\ are bounded. Proof. Throughout this proof, ∥⋅∥\|·\| and ⟨⋅,⋅⟩ ·,· denote the Frobenius norm and inner product on ℝn×mR^n× m, i.e. ‖A‖2=tr(A⊤A)\|A\|^2=tr(A A) and ⟨A,B⟩=tr(A⊤B) A,B =tr(A B). For every C≥f(x0)−f∗C≥ f(x_0)-f^*, fix an open neighborhood D of x∗x^* with ¯⊂int(C) D (D_C) where C:=f≤f∗+CD_C:=\f≤ f^*+C\, and define :=C∖K:=D_C D, r:=dist(¯,Cc)>0r_D:=dist( D,D_C^c)>0, γ(C):=minx∈‖∇f(x)‖>0,Γ(C):=maxx∈C‖∇f(x)‖<∞,γ(C):= _x \|∇ f(x)\|>0, (C):= _x _C\|∇ f(x)\|<∞, (39) and β∗(C):=γ(C)2/(Γ(C)mγ(C)2+ε2)1+γ(C)2/(Γ(C)mγ(C)2+ε2)∈(0,1).β^*(C):= γ(C)^2/\! ( (C) m γ(C)^2+ ^2 )1+γ(C)^2/\! ( (C) m γ(C)^2+ ^2 )\;∈\;(0,1). (40) Since β∗(C)∈(0,1)β^*(C)∈(0,1) for every C≥f(x0)−f∗C≥ f(x_0)-f^*, the supremum β∗:=supC≥f(x0)−f∗β∗(C)β^*:= _C\,≥\,f(x_0)-f^*β^*(C) (41) is well-defined and satisfies β∗∈(0,1]β^*∈(0,1]. Now fix any β∈[0,β∗)β∈[0,β^*) and define the set of admissible levels (β):=C≥f(x0)−f∗:β<β∗(C),C(β):= \C≥ f(x_0)-f^*\;:\;β<β^*(C) \, (42) which is non-empty by definition of β∗β^*: since β<β∗=supCβ∗(C)β<β^*= _Cβ^*(C), there exists at least one C∈(β)C (β). The remainder of the proof holds for every C∈(β)C (β); we write γ:=γ(C)γ:=γ(C) and Γ:=Γ(C) := (C) for brevity. Write sk:=xk−xk−1s_k:=x_k-x_k-1 and dε(y):=(∇f(y)∇f(y)⊤+ε2I)−1/2∇f(y)d_ (y):=(∇ f(y)∇ f(y) + ^2I)^-1/2∇ f(y). If ∇f(y)=UΣV⊤∇ f(y)=U V is the SVD with singular values σ1≥⋯≥σm≥0 _1≥·s≥ _m≥ 0, then dε(y)=U(Σ2+ε2I)−1/2ΣV⊤d_ (y)=U( ^2+ ^2I)^-1/2 V with singular values σi/σi2+ε2<1 _i/ _i^2+ ^2<1, so ‖dε(y)‖2=∑i=1mσi2σi2+ε2<m.\|d_ (y)\|^2= _i=1^m _i^2 _i^2+ ^2<m. (43) The recursion sk+1=βsk−αdε(yk)s_k+1=β s_k-α d_ (y_k) with s0=0s_0=0 and (43) give, by induction, ‖sk‖≤αm1−β=:s¯,‖yk−xk‖≤ρ:=αβm1−β,∀k≥0,\|s_k\|≤ α m1-β=: s, \|y_k-x_k\|≤ρ:= αβ m1-β, ∀\,k≥ 0, (44) unconditionally for any β∈[0,1)β∈[0,1). For y∈y , by direct computation using the SVD, ⟨∇f(y),dε(y)⟩=∑i=1mσi2σi2+ε2. ∇ f(y),d_ (y) = _i=1^m _i^2 _i^2+ ^2. (45) The function ϕ(t):=t/t+ε2φ(t):=t/ t+ ^2 is concave on [0,∞)[0,∞) with ϕ(0)=0φ(0)=0, hence subadditive, so ∑i=1mϕ(σi2)≥ϕ(∑i=1mσi2)=‖∇f(y)‖2‖∇f(y)‖2+ε2≥γ2γ2+ε2, _i=1^mφ( _i^2)≥φ\! ( _i=1^m _i^2 )= \|∇ f(y)\|^2 \|∇ f(y)\|^2+ ^2≥ γ^2 γ^2+ ^2, (46) where the last inequality uses that t↦t2/t2+ε2t t^2/ t^2+ ^2 is increasing and ‖∇f(y)‖≥γ\|∇ f(y)\|≥γ on K. By L-smoothness and (43), f(y−αdε(y)) f(y-α d_ (y)) ≤f(y)−α⟨∇f(y),dε(y)⟩+L2α2‖dε(y)‖2 ≤ f(y)-α ∇ f(y),d_ (y) + L2α^2\|d_ (y)\|^2 ≤f(y)−(αγ2γ2+ε2−Lm2α2) ≤ f(y)- ( αγ^2 γ^2+ ^2- Lm2α^2 ) =:f(y)−cα, =:f(y)-c_α, (47) where cα>0c_α>0 for all α∈(0,α1)α∈(0, _1) with α1(C):=2γ2/(Lmγ2+ε2)>0 _1(C):=2γ^2/(Lm γ^2+ ^2)>0. Next, we show yk∈Cy_k _C implies xk+1,yk+1∈Cx_k+1,y_k+1 _C. We consider two cases. Case 1: yk∈y_k . Since ‖xk+1−yk‖=‖αdε(yk)‖≤αm\|x_k+1-y_k\|=\|α d_ (y_k)\|≤α m and ‖yk+1−yk‖≤‖yk+1−xk+1‖+‖xk+1−yk‖≤ρ+αm=s¯\|y_k+1-y_k\|≤\|y_k+1-x_k+1\|+\|x_k+1-y_k\|≤ρ+α m= s, both xk+1x_k+1 and yk+1y_k+1 remain in CD_C provided s¯<r s<r_D, i.e. αm1−β<r(C). α m1-β<r_D(C). (48) Case 2: yk∈y_k . By (47), f(xk+1)≤f(yk)−cα≤f∗+Cf(x_k+1)≤ f(y_k)-c_α≤ f^*+C, so xk+1∈Cx_k+1 _C and ‖∇f(xk+1)‖≤Γ\|∇ f(x_k+1)\|≤ . By L-smoothness, f(yk+1)≤f(xk+1)+Γρ+L2ρ2f(y_k+1)≤ f(x_k+1)+ ρ+ L2ρ^2. For yk+1∈Cy_k+1 _C it suffices to have Γρ+L2ρ2<cα ρ+ L2ρ^2<c_α, since then f(yk+1)<f(xk+1)+cα≤f(yk)≤f∗+Cf(y_k+1)<f(x_k+1)+c_α≤ f(y_k)≤ f^*+C, where the last inequality in the chain comes from yk∈y_k . We therefore require Γρ+L2ρ2<cα. ρ+ L2ρ^2<c_α. (49) We next look into admissible ranges for α. Condition (48) holds for all α<(1−β)r(C)/mα<(1-β)r_D(C)/ m. Substituting ρ=αβm/(1−β)ρ=αβ m/(1-β) into (49) and rearranging, the condition holds whenever α<α2(β,C)α< _2(β,C) where α2(β,C):=2(1−β)2Lm(1−2β+2β2)(γ2γ2+ε2−Γβm1−β). _2(β,C):= 2(1-β)^2Lm(1-2β+2β^2) ( γ^2 γ^2+ ^2- β m1-β ). (50) Since C∈(β)C (β), we have β<β∗(C)β<β^*(C), which is precisely the condition that makes the bracket in (50) strictly positive, so α2(β,C)>0 _2(β,C)>0. Defining α∗(β,C):=minα1(C),α2(β,C),(1−β)r(C)m>0,α^*(β,C):= \! \ _1(C),\; _2(β,C),\; (1-β)r_D(C) m \>0, (51) all three terms of which are strictly positive for every C∈(β)C (β), both (48) and (49) hold for all α∈(0,α∗(β,C))α∈(0,α^*(β,C)). Since α∗(β,C)>0α^*(β,C)>0 for every C∈(β)C (β) and (β)C(β) is non-empty, define α∗(β):=supC∈(β)α∗(β,C)> 0.α^*(β):= _C\,∈\,C(β)α^*(β,C)\;>\;0. (52) For any α∈(0,α∗(β))α∈(0,α^*(β)), by definition of the supremum there exists C∈(β)C (β) such that α<α∗(β,C)α<α^*(β,C), and hence both (48) and (49) are satisfied for this C. Since C≥f(x0)−f∗C≥ f(x_0)-f^*, we have x0,y0∈Cx_0,y_0 _C. The above shows yk∈Cy_k _C implies xk+1,yk+1∈Cx_k+1,y_k+1 _C for all α∈(0,α∗(β))α∈(0,α^*(β)) and β∈[0,β∗)β∈[0,β^*). By induction, xk,yk∈Cx_k,y_k _C for all k≥0k≥ 0. Since CD_C is compact, the sequences are bounded. ∎ A.10 Proof of Corollary 4 Statement. For the algorithm presented in (16), there exists β¯∈(0,1) β∈(0,1) and α¯:[0,β¯)→ℝ+ α:[0, β) _+, such that every β∈[0,β¯)β∈[0, β) and every α∈(0,α¯(β))α∈(0, α(β)) the sequences xk\x_k\ and yk\y_k\ are bounded and limk→∞‖∇f(yk)‖=0 _k→∞\|∇\!f(y_k)\|=0. Proof. Define β¯:=minβ3∗,β4∗∈(0,1) β:= \β^*_3,β^*_4\∈(0,1) and α¯(β):=minα3∗(β),α4∗(β)>0 α(β):= \α^*_3(β),α^*_4(β)\>0, where βi∗β^*_i and αi∗(⋅)α^*_i(·) are as in Theorem i for i=3,4i=3,4. Fix any β∈[0,β¯)β∈[0, β) and α∈(0,α¯(β))α∈(0, α(β)). Since β<β4∗β<β^*_4 and α<α4∗(β)α<α^*_4(β), Theorem 4 guarantees that xk\x_k\ and yk\y_k\ remain in the compact set CD_C for some C∈(β)C (β). In particular, the eigenvalues of Pk=(∇f(yk)∇f(yk)⊤+ε2I)1/2P_k=(∇\!f(y_k)∇\!f(y_k) + ^2I)^1/2 lie in [ε,Γ(C)2+ε2][ , (C)^2+ ^2] for all k≥0k≥ 0, PkP_k is bounded with pmin=εp_ = and pmax=Γ(C)2+ε2p_ = (C)^2+ ^2. Since β<β3∗β<β^*_3 and α<α3∗(β)α<α^*_3(β), Theorem 3 then gives limk→∞‖∇f(yk)‖=0 _k→∞\|∇\!f(y_k)\|=0. ∎