Momentum in large-batch training: Polyak enlarges the critical batch size, Nesterov improves data efficiency

Momentum in large-batch training: Polyak enlarges the critical batch size, Nesterov improves data efficiency Jia-Nan Wang1,∗

Kairui Li1,∗

Lei Wu1,2,†

Peking University 2 AI for Science Institute, Beijing 3 University of Pennsylvania

jiananwang25@stu.pku.edu.cn, zixunh@wharton.upenn.edu kaisms@stu.pku.edu.cn, leiwu@math.pku.edu.cn

Abstract We study when and how momentum improves large-batch training in the one-pass regime, using power-law kernel regression as a tractable setting. We first characterize risk stability through the critical learning rate, defined as the largest learning rate for stable training, and crit crit crit ≂ min{1, B β (1 − ρ)}, where B is obtain ηsgd ≂ 1, ηpolyak ≂ min{1, B(1 − ρ)}, and ηnesterov the batch size, ρ is the momentum factor, and β > 1 is the capacity exponent. Within this admissible region, we derive scaling laws for the full risk dynamics, capturing the progression from an early transient, through power-law decay, to a noise floor. We then minimize the final-step risk over the admissible learning rates and momentum factors under a fixed data budget, yielding a three-regime batch-size phase diagram that reveals how the role of momentum changes with batch size. Notably, Polyak enlarges the critical batch size, the largest batch size preserving the best small-batch data-scaling exponent, thereby enabling greater parallelism without sacrificing data efficiency. In contrast, Nesterov achieves better data efficiency in the large batch regime because its look-ahead mechanism suppresses noise accumulation. Numerical experiments validate the predicted stability boundaries, risk dynamics, and batch-size phase diagram.

(a) Theory large batch

N=P=S

(b) Simulation ultra-large batch

N>P>S

small batch

N=P>S

Final excess risk

small batch

Data-scaling exponent

arXiv:2609.02728v1 [stat.ML] 2 Sep 2026

1

Zixun Huang3,∗

Nesterov Polyak SGD

large batch

ultra-large batch

10−1 6 × 10−2

Nesterov Polyak SGD

3 × 10−2 2 × 10−2

0

b1

b2

b3

Batch exponent b

1

100

101

Batch size B

102

103

Figure 1: Batch-size dependence of momentum benefits: Polyak enlarges the critical batch size, while Nesterov improves data efficiency. (a) Theoretical batch-size phase diagram from Theorem 6.3, showing the optimal data-scaling exponent as a function of the batch exponent b, defined by the scaling B = Db+o(1) . (b) Best-tuned final excess risk versus batch size B at fixed data budget D = 214 and (s, β) = (0.3, 4). The simulation recovers the same qualitative three-regime structure predicted by the theory. Experimental setup and additional results are provided in Section 7. ∗ †

Equal contribution. Corresponding author.

1

1

Introduction

Momentum has become a central mechanism in optimizer design for modern large-scale training, underlying popular optimizers such as Adam (Kingma & Ba, 2015) and Muon (Jordan et al., 2024; Liu et al., 2025). Classical theory views momentum primarily as an acceleration mechanism, quantified through reduced iteration complexity (Nesterov, 1983; Su et al., 2016; Wibisono et al., 2016; Wilson et al., 2021; Gitman et al., 2019; Velikanov et al., 2023; Wang et al., 2024; Liu & Belkin, 2018; Gupta et al., 2024). Yet iteration complexity alone does not fully capture the role of momentum in modern large-scale training, where the optimization regime itself changes as training resources, such as data and hardware memory, scale up. Modern large-scale training, particularly language-model pretraining, often operates in a one-pass training regime (Brown et al., 2020). In this regime, the batch size B controls not only hardware parallelism but also the optimization horizon: under a fixed data budget D, the number of sequential updates is D T = . B Larger batches improve hardware utilization and training throughput, but necessarily reduce the number of optimization steps available per data pass. This creates a fundamental tradeoff between hardware efficiency and data efficiency. Indeed, recent language-model studies show that the relative advantages of different optimizers can change substantially with batch size, with some optimizers exhibiting markedly stronger data efficiency in the large-batch regime (Marek et al., 2026; Shah et al., 2025; Semenov et al., 2025). These observations motivate a systematic characterization of the role of momentum across batch sizes. We therefore ask: Question: Can momentum enable larger batches without sacrificing data efficiency, or even improve data efficiency itself ? We address this question for two canonical momentum methods, Polyak (Polyak, 1964) and Nesterov (Nesterov, 1983), using power-law kernel (PLK) regression as a tractable learning problem. Its infinite-dimensional power-law spectrum spans a broad range of learning timescales, making the effects of momentum across spectral directions analytically accessible. The problem is characterized by two exponents: the capacity exponent β > 1, which controls the decay of the feature spectrum, and the source exponent s > 0, which controls how the target-function energy is distributed across the spectrum. Related power-law settings have been widely used to understand scaling laws in large-scale training across a variety of settings (Lin et al., 2024; Bahri et al., 2024; Bordelon et al., 2025; Li et al., 2026a; Wang et al., 2026; Li et al., 2026b). Building on the functional scaling law (FSL) framework of Li et al. (2026a), we extend the analysis to momentum methods. Our analysis of one-pass training follows the hierarchy AB −→ E(k, η, ρ, B) −→

min (η,ρ)∈AB

E(T, η, ρ, B) ,

T =

D . B

Here, AB denotes the admissible set of learning rates η and momentum factors ρ at batch size B, and E(k, η, ρ, B) denotes the expected excess risk. The three stages correspond, respectively, to 1) identifying the admissible hyperparameter region through stability analysis, 2) deriving scaling laws that characterize the risk dynamics, and 3) optimizing the final-step risk subject to both the admissibility constraint (η, ρ) ∈ AB and the one-pass constraint T = D/B. Our contributions.

Our main contributions are summarized as follows.

• Risk stability. We characterize risk stability through the critical learning rate η crit , the largest learning rate for which the stochastic dynamics remain stable (Section 4). In the 2

large-batch, high-momentum regime, we obtain the sharp scalings crit ηsgd ≂ 1,

crit ηpolyak ≂ min{1, B(1 − ρ)},

crit ηnesterov ≂ min{1, B β (1 − ρ)}.

(1)

Thus, stronger momentum requires a correspondingly larger batch to maintain stability. This compensation scales as B for Polyak but as B β for Nesterov. Since β > 1, Nesterov permits substantially higher momentum in large-batch training. • Scaling laws for risk dynamics. Within the risk-stable region, we derive sharp scaling laws for the risk dynamics of Polyak and Nesterov (Section 5). In the high-momentum regime, the risk transitions from an early exponential transient to power-law decay and eventually saturates at a noise floor. Writing ηρ := η/(1 − ρ), after the transient the dominant terms take the simple form Esgd (k) ≂ (kη)−s +

σ2η , B

Epolyak (k) ≂ (kηρ )−s +

σ 2 ηρ , B

while Nesterov retains the same accelerated signal term but reduces the noise floor to Enesterov (k) ≂ (kηρ )−s +

σ2 min{ηρ , ηρ1/β }. B

Thus, Polyak preserves the SGD signal–noise tradeoff while operating at the enlarged effective learning rate ηρ , whereas Nesterov further suppresses noise accumulation at large ηρ through the curvature-adaptive damping induced by its look-ahead mechanism. • Data scaling and batch-size phase diagram. Let B = Db+o(1) . Optimizing the remaining hyperparameters yields three characteristic batch exponents b1 < b2 < b3 and a three-regime phase diagram (Section 6). For b ⩽ b1 , SGD, Polyak, and Nesterov achieve the same data-scaling exponent. For b1 < b < b3 , Nesterov outperforms Polyak, which in turn outperforms SGD, with Nesterov attaining its best data-scaling rate at b = b2 . For b ⩾ b3 , Polyak and Nesterov again share the same exponent and both outperform SGD. Here b1 and b3 mark the critical batch exponents of SGD and Polyak, respectively, showing that Polyak preserves the best small-batch scaling over a substantially wider batch-size range, while Nesterov further improves the achievable data-scaling rate within the large batch regime. Figure 1 provides a complete view of how the relative advantage of the three methods evolves with batch size. Finally, Section 7 empirically validates these theoretical predictions. Why momentum helps large-batch training. The batch-size dependence of momentum can be understood through the per-sample effective learning rate (ELR), ( η, SGD, ηeff ηsample := , ηeff = B ηρ := η/(1 − ρ), momentum. Ignoring the early transient and using T = D/B, SGD and Polyak obey the same data-level signal–noise tradeoff, E(D) ≂ (Dηsample )−s + σ 2 ηsample . The difference is the admissible range of the per-sample ELR. Stability condition (1) gives ηsample ≲ 1/B for SGD, whereas Polyak allows ηsample ≲ 1. Thus, as the batch size grows, the admissible per-sample ELR of SGD shrinks as 1/B, while Polyak can keep it at order one by increasing ηρ proportionally with B. Polyak therefore extends the SGD signal–noise tradeoff to substantially larger batch sizes without changing its optimum, explaining both the enlarged critical batch size and the unchanged best data-scaling exponent. 3

Nesterov changes the picture more fundamentally. Its look-ahead mechanism suppresses noise accumulation at large effective learning rates, improving the signal–noise tradeoff itself and thereby achieving better data efficiency in the large-batch regime. At ultra-large batch sizes, however, the limited number of sequential updates T = D/B becomes the dominant bottleneck, and the two momentum methods recover the same signal-limited scaling.

2

Related work

Theory of momentum. Classical analyses explain momentum through accelerated convergence, Lyapunov functions, and continuous-time analysis (Nesterov, 1983; Su et al., 2016; Wibisono et al., 2016; Wilson et al., 2021; Shi et al., 2022; Kovachki & Stuart, 2021). Under stochastic training, subsequent work has studied convergence, stability, and asymptotic behavior (Vidyasagar, 2025; Liu et al., 2020; Yuan et al., 2016; Gitman et al., 2019; Velikanov et al., 2023; Zhang et al., 2024). However, existing theory does not characterize how the role of momentum changes with batch size under a fixed data budget. We instead treat batch size as a scaling variable and compare Polyak and Nesterov across the resulting training regimes. Scaling-law theory. Scaling laws characterize how learning performance scales with resources such as data, model, and compute (Hestness et al., 2017; Kaplan et al., 2020; Hoffmann et al., 2022; Grattafiori et al., 2024; DeepSeek-AI et al., 2026). A growing theoretical literature has derived such laws in linear, kernel, and related tractable models (Sharma & Kaplan, 2022; Hutter, 2021; Maloney et al., 2022; Wei et al., 2022; Jain et al., 2024; Michaud et al., 2023; Nam et al., 2024; Atanasov et al., 2026; Bahri et al., 2024; Dohmatob et al., 2024; Bordelon et al., 2024; Lin et al., 2024; Paquette et al., 2024; Bordelon et al., 2025; Yan et al., 2026; Lin et al., 2025; Sous & Winer, 2026; Li et al., 2026a; Zhang et al., 2025; Wang et al., 2026; Li et al., 2026b; Ding et al., 2025). Functional scaling laws further characterize risk dynamics through signal learning and noise accumulation (Li et al., 2026a). We extend this framework to momentum methods. Closely related, Ferbach et al. (2026a) study momentum scaling laws in noiseless power-law random features with fixed batch size, showing that Polyak preserves the SGD scaling exponents while dimension- and time-adapted Nesterov-type methods can improve them. We instead consider constant-order label noise and scale the batch size under a fixed data budget, obtaining a batch-size phase diagram for canonical Polyak and Nesterov momentum. Stability of gradient-based methods. Stability characterizes the admissible hyperparameter region and has been studied extensively for gradient-based methods, with connections to implicit bias, generalization, and optimization geometry (Zhu et al., 2019; Wu et al., 2018; Nar & Sastry, 2018; Ma & Ying, 2021; Wu et al., 2022; Wu & Su, 2023; Chemnitz & Engel, 2025; Mulayoff et al., 2021; Nacson et al., 2023; Qiao et al., 2024; Kaur et al., 2023; Wang & Wu, 2023). For momentum methods, its joint dependence on batch size and momentum remains less understood. We derive sharp risk stability boundaries for Polyak and Nesterov, revealing a much stronger batch-size compensation for high momentum in Nesterov. Concurrent work. Morwani et al. (2026), which appeared while this manuscript was in preparation, studies batch-size tradeoffs for stochastic momentum in finite-dimensional linear regression and derives lower bounds on contraction rates. These bounds suggest that heavy-ball may retain data efficiency over a larger batch-size range, but matching upper bounds are left open. In contrast, we derive sharp scaling laws for the full risk dynamics and, after optimizing the learning rate and momentum factor, obtain the data-scaling exponent as an explicit function of batch size. This yields a complete batch-size phase diagram, including a large-batch regime in which Nesterov outperforms Polyak.

4

3

Preliminaries

Notation. In this paper, we use ≂ to denote equivalence up to a constant factor, and ≲ (resp. ≳) indicates inequality up to a constant factor. For two non-negative functions f, g, we write f (x) ≲ g(x) (resp. f (x) ≳ g(x)) if there exists a constant C > 0 such that f (x) ⩽ Cg(x) (resp. g(x) ⩽ Cf (x)) for all x under consideration. We write f (x) ≂ g(x) if both f (x) ≲ g(x) and f (x) ≳ g(x) hold. Throughout the paper, ⟨·, ·⟩ and ∥·∥2 denote the standard inner product and the associated 2-norm, respectively.

3.1

Power-law kernel regression

We consider a power-law kernel (PLK) regression on a Hilbert space H with the input x ∼ D, and the observed label is y = ⟨ϕ(x), θ ∗ ⟩+ϵ, where ϵ ∼ N (0, σ 2 ) is independent of x. Here we assume the feature function ϕ(x) to be Gaussian: ϕ(x) ∼ N (0, H), where H = Ex∼D [ϕ(x) ⊗ ϕ(x)] denotes the feature covariance operator. Let {λj }j⩾1 denote the eigenvalues of H, arranged in decreasing order. We learn the target function using the linear student  model f (x; θ) = ⟨ϕ(x), θ⟩ 1 2 by minimizing the population risk: R(θ) = 2 Ex,y (f (x; θ) − y) . A direct calculation gives 2 R(θ) = E(θ) + σ2 , where E(θ) = 12 ⟨θ − θ ∗ , H(θ − θ ∗ )⟩ is the excess risk. Since SGD, Polyak momentum, and Nesterov momentum are orthogonally equivariant, we work in the eigenbasis of the covariance operator without loss of generality; see Appendix A.2. Assumption 3.1 (Diagonal covariance). The covariance operator is diagonal: λ1 ⩾ λ2 ⩾ · · · > 0.

H = diag(λ1 , λ2 , . . .),

(2)

Throughout, we focus on the label-noise-dominated setting with a constant-order noise level. Assumption 3.2 (Constant-order label noise). Assume σ ≂ 1. Assumption 3.3 (Source and capacity conditions). There exist β > 1 and s > 0 such that λj ≂ j −β ,

λj |θj∗ |2 ≂ j −(1+sβ) .

Here β is the capacity exponent, which controls the decay rate of the feature spectrum, and s is the source exponent, which describes how the target energy decays across spectral directions. We also refer to s as the relative difficulty of the task: smaller s means that more target energy lies in low-curvature, small-eigenvalue directions, making the problem harder to learn. We remark that similar assumptions have been widely used in the analysis of kernel methods (Caponnetto & Vito, 2005; Caponnetto & De Vito, 2007; Spigler et al., 2020; Bordelon et al., 2020; Maloney et al., 2022). Our work builds upon and extends this line of research.

3.2

Optimization methods

At each step, we draw a fresh mini-batch Bk = {(xk,i , yk,i )}B i=1 and denote the mini-batch 1 P 2 gradient by g(θ; Bk ) = ∇θ ( 2B (x,y)∈Bk (f (x; θ) − y) ). SGD updates as follows SGD:

θk+1 = θk − ηg(θk ; Bk ),

(3)

where η is the learning rate. Polyak momentum, also called the heavy-ball method, is given by Polyak:

θk+1 = θk − ηg(θk ; Bk ) + ρ(θk − θk−1 ),

(4)

where ρ ∈ [0, 1) is the momentum factor. Nesterov momentum differs from Polyak momentum by evaluating the gradient at a look-ahead point: Nesterov:

θkla = θk + ρ(θk − θk−1 ) 5

(look ahead)

θk+1 = θk − ηg(θkla ; Bk ) + ρ(θk − θk−1 ),

(5)

where θkla is the look-ahead point obtained by extrapolating the current iterate along the previous update direction. This look-ahead mechanism distinguishes Nesterov from Polyak. For clarity, we decompose g(θ; Bk ) = ∇R(θ) + ξk (θ), where ∇R(θ) denotes the population gradient and ξk (θ) represents the stochastic gradient noise satisfying E[ξk (θ) | θ] = 0,

E[ξk (θ) ⊗ ξk (θ) | θ] =

1 Σ(θ), B

where Σ(θ) denotes the gradient-noise covariance of batch size 1. We next recall a noise–curvature alignment property that will be used throughout our analysis. Proposition 3.4 (Noise–curvature alignment). For any unit vector ν ∈ H, h i
1 E |⟨ξk (θk ), ν⟩|2 θk ≂ E(θk ) + σ 2 ⟨ν, Hν⟩ . B

(6)

Here, ⟨ν, Hν⟩ is the curvature along ν. The two terms in Proposition 3.4 have distinct origins and roles. The term proportional to E(θk ) arises from mini-batch sampling, depends on the current iterate, and vanishes at the optimum, acting as multiplicative noise. The term proportional to σ 2 arises from label noise and persists at the optimum, acting as additive noise. Both are curvature-aligned but affect the dynamics differently: multiplicative noise feeds back into the dynamics and determines stability, whereas additive noise determines the asymptotic noise floor. This noise–curvature alignment has been observed in neural networks (Zhu et al., 2019; Wu et al., 2022) and established rigorously in toy models (Wang & Wu, 2023). We use this constitutive relation throughout our analysis; see Appendix A.4 for a proof.

Regime-dependent damping in momentum methods

Consider the full-batch Polyak on a one-dimensional quadratic function f (x) = λx2 /2, where λ > 0 denotes the curvature. Since ∇f (x) = λx, Polyak updates: xk+1 = xk − ηλxk + ρ(xk − xk−1 ).

(7)

1

Momentum

3.3

= (1 fixed

)2

Underdamped

The corresponding characteristic polynomial is r2 − (1 + ρ − ηλ)r + ρ = 0. The relevant transition from real to complex Overdamped √ roots occurs at ηλ = (1 − ρ)2 . 00 1 For fixed η and ρ, different curvature directions can exhibit Effective curvature two qualitatively different damping behaviors. √ • Overdamped regime: λη < (1 − ρ)2 . In this regime, the characteristic roots r± are ηλ real with the dominant root satisfying r+ ≈ 1 − 1−ρ when 1 − ρ = o(1). Hence, for sufficiently large k, it holds that   ηλ k k k xk = c+ r+ + c− r− ≈ c+ r+ ≈ c+ exp − k , 1−ρ where c± depend on the initialization. Hence, its update behaves like gradient descent with the momentum-rescaled learning rate ηρ = η/(1 − ρ). √ • Underdamped regime: λη > (1 − ρ)2 . In this regime, the characteristic roots √ form a complex-conjugate pair with modulus ρ. The iterate therefore has the form   xk = ρk/2 c1 cos(kω) + c2 sin(kω) , where c1 , c2 , and ω depend on the initialization and curvature. Thus, the iterate exhibits oscillatory convergence to zero, with the amplitude decaying as ρk/2 . 6

The key difference between the two regimes lies in how curvature affects the dynamics. In the overdamped regime, the decay rate ηλ/(1 − ρ) is proportional to the curvature, so flatter directions converge more slowly. In the underdamped regime, the oscillation frequency depends on the curvature, but the amplitude envelope ρk/2 does not. Nesterov look-ahead induces curvature-adaptive damping. On the same quadratic objective, Nesterov updates using the gradient at the look-ahead point xla k = xk + ρ(xk − xk−1 ): xk+1 = xk − ηλ[xk + ρ(xk − xk−1 )] + ρ(xk − xk−1 ) = xk − ηλxk + ρ(1 − ηλ)(xk − xk−1 ),

(8)

which is equivalent to Polyak with the momentum factor ρ(1 − ηλ). For an underdamped direction, it accelerates convergence: ρk/2 → [ρ(1 − ηλ)]k/2 , which damps sharper directions more rapidly. By contrast, in sufficiently flat directions where ηλ ≪ 1, we have ρ(1 − ηλ) ≈ ρ, and Nesterov retains approximately the same dynamics as Polyak. Remark 3.5 (Spectral√decomposition into damping regimes). For a fixed pair (η, ρ), the damping (1− ρ)2

threshold λcut := partitions the spectrum into two parts. Directions with λj > λcut are η underdamped, whereas directions with λj < λcut are overdamped. Hence, in PLK regression, the risk dynamics can be analyzed by treating the underdamped spectral head and overdamped spectral tail separately and then combining their contributions.

4

Risk stability and critical learning rates

In this section, we characterize the maximal learning rates that ensure stable training and study their dependence on the momentum factor and batch size. We measure stability directly through the population risk. Definition 4.1 (Risk stability). Consider a stochastic optimization algorithm for minimizing a population risk R(·), producing model parameters {θk }k⩾0 . We say that the algorithm is risk stable if, for every admissible initialization with finite population risk, sup E[R(θk )] < ∞. k⩾0

Definition 4.2 (Critical learning rate). Consider a stochastic optimization algorithm A parameterized by a learning rate η and additional hyperparameters h. For fixed h, let crit ηA (h) := sup {η > 0 : A is risk stable} .

Hence, the critical learning rate is the largest learning rate compatible with risk stability. crit . In our setting, h = B When the dependence on h is clear from context, we simply write ηA for SGD and h = (B, ρ) for Polyak and Nesterov. The following theorem characterizes how the critical learning rate scales with batch size and momentum; the proof is deferred to Appendix B. Theorem 4.3 (Critical learning-rate scaling). For sufficiently large B and sufficiently small 1 − ρ, the critical learning rates satisfy crit ηsgd ≂ 1,

(9)

crit ηpolyak ≂ min{1, B(1 − ρ)}, crit ηnesterov ≂ min{1, B β (1 − ρ)},

where the implicit constants are independent of B and ρ.

7

(10) (11)

Joint dependence on momentum and batch size. Theorem 4.3 reveals a direct coupling among the learning rate, momentum factor, and batch size in determining training stability. Increasing ρ toward one reduces the critical learning rate, whereas increasing the batch size crit crit enlarges it. In the high-momentum regime, ηpolyak ≂ B(1 − ρ), ηnesterov ≂ B β (1 − ρ), before reaching their order-one ceilings. Thus, the batch-size effect on stability is substantially stronger for Nesterov, with a B β dependence rather than the B dependence of Polyak. This difference originates from the curvature adaptivity induced by Nesterov’s look-ahead mechanism. Multiplicative-noise feedback governs stability through the renewal kernel. For all three methods, Appendix BPshows that the risk
dynamics admit, up to constant factors, a 2 renewal equation ℓk ≂ hk + k−1 j=0 Kk−1−j ℓj + σ , where ℓk is the risk moment, hk is the free response, and the renewal kernel K describes how stochastic perturbations propagate through the optimizer dynamics. Additive and multiplicative noise play different roles: additive label noise acts only as an external forcing, whereas multiplicative noise feeds the current risk back through the renewal kernel, producing the feedback term K ∗ ℓ. Risk stability is therefore determined by the strength of this feedback, characterized by ∥K∥1 < 1; see Appendix B.1. Since K is independent of the label-noise scale σ, the stability condition in Theorem 4.3 is also independent of σ.

5

Scaling laws for risk dynamics across damping regimes

In this section, we derive sharp scaling laws for the full risk dynamics of Polyak and Nesterov, explicitly characterizing their joint dependence on the hyperparameters (η, ρ, B). The analysis follows how increasing ρ drives transitions in damping and stability: from a fully overdamped regime with SGD-like dynamics, to a mixed-damping regime, and eventually to instability. As discussed in Section 3.3, for a single eigendirection with curvature λ, the overdamped–underdamped √ transition occurs exactly at (1 − ρ)2 = ηλ. For the scaling analysis below, where ρ → 1, we use the equivalent high-momentum form ηλ ≂ (1 − ρ)2 . This motivates the curvature cutoff λcut ≂

(1 − ρ)2 , η

Iover = {j : λj ≲ λcut } ,

Iunder = {j : λj ≳ λcut } .

Hence, the sets Iover and Iunder consist of flat overdamped directions and sharp underdamped directions, respectively. Since the PLK spectrum contains arbitrarily small eigenvalues, overdamped directions always remain present. As ρ increases, the first underdamped directions appear when λcut ≲ λ1 , corresponding to √ 1 − ρ ≂ η. Combining this damping threshold with the risk stability conditions from Section 4 yields three regimes: √ ( fully overdamped: 1 − ρ ≳ η, 1, Polyak, √ q mixed damping: η/B ≲ 1 − ρ ≲ η, q= β, Nesterov. unstable: 1 − ρ ≲ η/B q , √ Thus, the damping threshold 1 − ρ ≂ η and the risk stability threshold 1 − ρ ≂ η/B q delimit the three regimes. For reference, the SGD scaling law for PLK regression is E[E(θk )] ≂ (ηk)−s +

σ2 η, B

kη ≳ 1.

(12)

This follows from Li et al. (2026a, Theorem 5.2); see Appendix C for details. The first term, (ηk)−s , describes signal learning, with the decay rate determined by the source exponent s: 8

larger s corresponds to an easier target and faster learning. The second term, σ 2 η/B, captures stochastic noise accumulation and determines the asymptotic noise floor. The overdamped modes evolve with the rescaled learning rate ηρ :=

η , 1−ρ

which plays the role of the effective learning rate for momentum. The following theorem establishes the scaling law for Polyak momentum, whose proof is deferred to Appendix D. Theorem 5.1 (Scaling law for Polyak). Assume η ≲ 1 and k ≳ max(ηρ−1 , (1 − ρ)−1 ). Then  σ2  √  (ηρ k)−s + ηρ , 1 − ρ ≳ η, B E[E(θk )] ≂ (13)  σ2 η √  −s P (k) + (ηρ k) + ηρ , ≲ 1 − ρ ≲ η, B B where P (k) denotes the underdamped transient and satisfies 0 ⩽ P (k) ≲ ρk . If 1 − ρ ≲ η/B, the dynamics are not risk stable. Theorem 5.1 identifies three regimes as ρ increases: fully overdamped SGD-like dynamics with effective learning rate ηρ , mixed damping with an additional transient P (k), and eventually risk instability. In the mixed-damping regime, the two risk contributions have a direct spectral interpretation. The power-law term (ηρ k)−s arises from the overdamped flat modes in Iover , whereas P (k) comes from the underdamped sharp modes in Iunder and decays with envelope ρk . Increasing ρ therefore accelerates flat-mode learning through ηρ , while slowing the relaxation of sharp modes. Since ρk ≈ exp(−(1 − ρ)k), the latter occurs on the timescale   1 η τρ := log . 1−ρ (1 − ρ)2 Once the underdamped transient becomes negligible, the Polyak risk law reduces to the SGD law under the replacement η 7→ ηρ , leaving the signal–noise tradeoff unchanged. The distinction lies in the admissible effective learning rate. By Theorem 4.3, SGD requires η ≲ 1, whereas Polyak allows ηρ ≲ B, so its admissible range grows linearly with the batch size. As we show in Section 6, this enlarged range increases the critical batch size but does not improve the best data-scaling exponent. For Nesterov, the same damping structure persists, but its look-ahead mechanism modifies both the underdamped transient and noise accumulation in the mixed-damping regime. Theorem 5.2 (Scaling law for Nesterov). Assume η < 1 and k ≳ max(ηρ−1 , (1 − ρ)−1 ). Then  2  √ (η k)−s + σ η ,  1 − ρ ≳ η, ρ ρ B E[E(θk )] ≂ (14)  σ2 η √ 1/β  N (k) + (ηρ k)−s + min{ηρ , ηρ }, ≲ 1 − ρ ≲ η, B Bβ where N (k) denotes the underdamped transient and satisfies 0 ⩽ N (k) ≲ min{1, (ηk)−s }ρk . If 1 − ρ ≲ η/B β , the dynamics are not risk stable. The proof is deferred to Appendix E. In the fully overdamped regime, Nesterov and Polyak obey the same scaling law, both behaving like SGD with effective learning rate ηρ . Their difference in the mixed-damping regime stems from the curvature-adaptive damping induced by Nesterov’s look-ahead mechanism, which damps sharper directions more strongly; see Eq. (8). 9

This suppresses the underdamped transient from P (k) ≲ ρk

N (k) ≲ min{1, (ηk)−s }ρk ,

to

while leaving the asymptotic power-law signal term unchanged at (ηρ k)−s . The same mechanism also suppresses stochastic noise accumulation. When ηρ ≫ 1, the noise floor scales as σ2 ηρ B

σ 2 1/β η B ρ

for Polyak,

for Nesterov.

For asymptotic data scaling, this reduction in noise accumulation is the key effect.

6

Optimal data scaling across batch-size regimes

We now optimize the final-step risk using the scaling laws derived above, subject to the one-pass constraint T = D/B and the stability constraints. We focus on the large-batch setting, allowing B to grow with the data budget D rather than remain fixed. Definition 6.1 (Batch-size exponent). For a sequence of training problems with D → ∞, the batch size has exponent b ∈ [0, 1) if B = Db+o(1) ,

log B . D→∞ log D

equivalently

b = lim

Definition 6.2 (Batch-size-dependent data-scaling exponent). Fix b ∈ [0, 1), so that B = Db+o(1) and T = D/B. For an optimization method with admissible hyperparameter region AB , define its optimal data-scaling exponent r(b) by inf (η,ρ)∈AB

E[E(θT )] = D−r(b)+o(1) .

Here, ρ is omitted for SGD. For the three methods considered, we write rsgd (b), rpolyak (b), and rnesterov (b). A larger r(b) indicates better data efficiency at batch-size exponent b. Define the three transition exponents b1 =

s , s+1

b2 =

2s + 1/β , 2s + 1 + 1/β

b3 =

2s + 1 . 2s + 2

(15)

Since β > 1 and s > 0, 0 < b1 < b2 < b3 < 1. The following theorem gives the optimal data-scaling exponents across batch-size regimes; the proof is deferred to Appendix F. Theorem 6.3 (Batch-size-dependent data scaling). For B = Db+o(1) with b ∈ [0, 1),    s ,  s , 0 ⩽ b < b , 0 ⩽ b < b3 , 1 s+1 s+1 rsgd (b) = rpolyak (b) = s(1 − b), b ⩽ b < 1, 2s(1 − b), b ⩽ b < 1, 1 3 and rnesterov (b) =

 s   s+1 ,  

s(1+(β−1)b)

sβ+1     2s(1 − b),

0 ⩽ b < b1 , ,

b1 ⩽ b < b 2 , b2 ⩽ b < 1.

This theorem shows that the achievable data efficiency depends sharply on batch-size scaling. The resulting phase structure is most clearly seen in Figure 1 (left), which visualizes the three 10

data-scaling exponents as functions of the batch exponent b and makes their relative ordering across regimes immediate. The three transition exponents play distinct roles: b1 is the critical batch exponent of SGD and marks the onset of Nesterov’s improvement, b2 maximizes the Nesterov data-scaling exponent, and b3 is the critical batch exponent of Polyak. • Small batch (b ⩽ b1 ). All three methods achieve the same optimal data-scaling exponent: rnesterov (b) = rpolyak (b) = rsgd (b). • Large batch (b1 < b < b3 ). The three methods separate: rnesterov (b) > rpolyak (b) > rsgd (b), with Nesterov attaining its best data-scaling exponent at b = b2 . • Ultra-large batch (b ⩾ b3 ). Polyak and Nesterov share the same exponent and outperform SGD: rnesterov (b) = rpolyak (b) > rsgd (b). Thus, Polyak enlarges the batch-size range over which the best small-batch data scaling is preserved, whereas Nesterov further improves the achievable data-scaling exponent in the large-batch regime.

7

Numerical validations

In this section, we numerically validate our theoretical predictions. Experimental details and additional results are provided in Appendix G. Risk stability boundaries. We validate the stability predictions of Theorem 4.3 for (s, β) = crit crit ≂ B and ηnesterov ≂ B β before (3, 2). Figure 2 (left) fixes ρ and varies B, confirming ηpolyak crit remains order one. Figure 2 (middle) fixes B = 32 and varies 1 − ρ, saturation, while ηsgd confirming the predicted linear dependence of both momentum stability boundaries on 1 − ρ. Additional results are provided in Appendix G.1. Risk-dynamics scaling laws. We next test the full risk dynamics predicted by Theorems 5.1 and 5.2. For fixed (η, ρ, B) and (s, β) = (0.3, 4), Figure 2 (right) verifies three key predictions: an underdamped transient whose transition timescale has dominant (1 − ρ)−1 dependence, a common power-law signal decay (ηρ k)−s , and distinct asymptotic noise floors. In particular, once the transient becomes negligible, Polyak and Nesterov follow the same power law, while Nesterov saturates at a substantially lower noise floor. These are consistent with our theoretical predictions. Additional results across (s, β) and hyperparameter settings are provided in Appendix G.2. Batch-size phase diagram at a fixed data budget. We first test the three-regime structure predicted by Theorem 6.3 at a finite data budget. We fix D = 214 and, for each batch size B, tune (η, ρ) to minimize the final excess risk. Figure 1 (right) exhibits a progression across the three regimes. At small batch sizes, SGD, Polyak, and Nesterov achieve identical risk. As B increases, SGD deteriorates first, while Polyak maintains the small-batch performance over a wider range, reflecting its enlarged critical batch size. Nesterov goes further: in the large-batch regime, its risk continues to decrease beyond the best level attained by SGD and Polyak, before eventually increasing at ultra-large batch sizes. Thus, Polyak enlarges the data-efficient batch-size range, while Nesterov further improves data efficiency. 11

Critical learning rates

10−2

10−3

103

SGD Prediction Polyak Prediction Nesterov Prediction

22

23

24

25

26

27

Batch Size B

28

29

10 1

10 2 10 2

210

SGD Prediction Polyak Prediction Nesterov Prediction

29

28

27

26

25

Momentum Damping 1

102

10−1

10−2

Polyak Nesterov

10−2

24

Transition-time scaling

Transition step

10−1

Loss

10

Scaling law predictions

100

100

0

Stability Boundary

Stability Boundary η

Critical learning rates

1−ρ

Transition phase

Scaling law phase

100

104

102

Step

106

Figure 2: Numerical validation of stability and risk-dynamics scaling laws. Left: Critical learning rate versus batch size at fixed momentum. Markers show numerical estimates and dashed curves the theoretical predictions. Middle: Critical learning rate versus 1 − ρ at fixed batch size, confirming the predicted linear scaling before the order-one ceiling. Right: Polyak and Nesterov risk dynamics, with empirical results shown by thin translucent lines and theoretical predictions by thick solid lines. The theory captures the underdamped transient, power-law decay, and noise floor. The inset confirms the dominant
(1 − ρ)−1 dependence of the predicted transition timescale τρ ≂ (1 − ρ)−1 log η/(1 − ρ)2 , with fitted log-log slopes −1.00 for Polyak and −1.01 for Nesterov.

Data-scaling exponents across batch-size regimes. Having verified the global phase structure, we next test the predicted data-scaling exponents at representative points in each regime. We continue with (s, β) = (0.3, 4) and choose B ≂ 1, B ≂ D0.4 , and B ≂ D0.8 , representing small batch, large batch, and ultra-large batch, respectively. For each data budget D, we set (η, ρ) according to the theoretically prescribed scalings. Figure 3 shows that the empirical slopes closely match Theorem 6.3: all three methods share the same exponent for small batch; for large batch, Nesterov outperforms Polyak, which outperforms SGD; and for ultra-large batch, Polyak and Nesterov again share the same exponent and both outperform SGD. These results quantitatively validate the batch-dependent data-scaling exponents underlying the phase diagram. Experimental details and additional settings are provided in Appendix G.4. Small batch

Large batch

SGD Polyak Nesterov Prediction

6 × 10−2

Nesterov Prediction

Loss

6 × 10−2

Loss

Loss

10−1

Ultra-large batch SGD Polyak

4 × 10−2

10−1 9 × 10−2 8 × 10−2

3 × 10−2

7 × 10−2

104

Data D

105

104

Data D

105

SGD Polyak Nesterov Prediction

104

105

Data D

Figure 3: Data-scaling exponents across the three batch-size regimes. We choose representative batch-size scalings B ≂ 1 (left), B ≂ D0.4 (middle), and B ≂ D0.8 (right), and vary the data budget D. Markers show empirical final excess risks, while dashed lines indicate the scaling exponents predicted by Theorem 6.3. For these representative choices, the predicted ordering is recovered: Nesterov = Polyak = SGD at small batch, Nesterov > Polyak > SGD at large batch, and Nesterov = Polyak > SGD at ultra-large batch.

8

Conclusion and discussion

We developed a unified analysis of Polyak and Nesterov momentum in large-batch, one-pass training, characterizing their risk stability, risk dynamics, and optimal data scaling. A central conclusion is that the benefit of momentum depends critically on the batch size. At small batch sizes, neither method improves the data-scaling exponent of vanilla SGD. As the batch size grows, however, the two methods provide distinct benefits: Polyak preserves the small-batch

12

optimal data-scaling exponent over a substantially wider range of batch sizes than SGD, whereas Nesterov further improves the data-scaling exponent in the large-batch regime. At ultra-large batch sizes, training becomes limited by the number of sequential updates, and the two momentum methods again exhibit the same signal-limited scaling. Implications for momentum scheduling. Our scaling analysis allows the batch size B, learning rate η, and momentum factor ρ to scale with the data budget D, while keeping them fixed within each training run. A complete theory of time-varying schedules is therefore beyond the present analysis. Nevertheless, our stability results already have direct implications for momentum scheduling. Recent work has increasingly explored time-varying momentum schedules, while also observing that increasing momentum during training can become unstable under stochastic gradients and may require warmup or additional damping (Pagliardini et al., 2025; Yarotsky & Velikanov, 2025; Yarotsky, 2026; Ferbach et al., 2026b). Our results provide a complementary quantitative perspective: for the classical Polyak and Nesterov updates, the stability boundary predicts how long an increasing-momentum schedule can remain stable and how this horizon scales with the batch size. A particularly instructive example is the classical scaling 1 − ρk ≂ k −1 , which underlies Nesterov acceleration in deterministic optimization (Nesterov, 1983). Although Theorem 4.3 is derived for constant momentum, applying its stability boundary locally along this slowly varying schedule predicts     B Bβ η ≲ min 1, for Polyak, η ≲ min 1, for Nesterov. (16) k k Thus, for fixed η and B, the schedule may remain stable for a substantial fraction of training but eventually cross the stochastic stability boundary as ρk → 1. The corresponding instability horizons are predicted to scale as ∗ ≂ kpolyak

B , η

∗ knesterov ≂

Bβ . η

Prior work has observed that increasing the batch size can delay the instability of aggressive momentum schedules (Yarotsky & Velikanov, 2025); to our knowledge, the explicit batch-size scaling of this delayed-instability horizon, and its distinction between Polyak and Nesterov, has not been characterized previously. Fixed η, ρk = 1 − 1/(k + 1) The right panel illustrates this delayed-instability SGD 10 phenomenon for PLK regression with (s, β) = (0.3, 4) Polyak Nesterov under the schedule ρk = 1 − 1/(k + 1); experimental details are provided in Appendix G.5. Both methods 10− remain stable and benefit from the increasing momentum before eventually becoming unstable. Consistent with 10− the predicted horizons, Nesterov remains stable for substantially more iterations than Polyak because its 10− stability horizon scales as B β rather than B. This allows 10 10 10 10 Nesterov to exploit the increasing-momentum schedule over Step a longer training interval and reach a lower risk before instability occurs. This example highlights a constraint on momentum scheduling that is absent from deterministic acceleration theory: the benefit of increasing momentum depends not only on how strongly it accelerates optimization, but also on how long the resulting trajectory remains inside the stochastic stability region. It also suggests a simple tuning principle. Momentum, learning rate, and batch size should be scheduled jointly so that 1 − ρk remains above the stochastic stability scale, namely 1 − ρk ≳ ηk /Bk for Polyak and 1 − ρk ≳ ηk /Bkβ for Nesterov. Increasing momentum can therefore be supported either by increasing the batch size or by 0

Loss

1

2

3

0

13

2

4

6

decreasing the learning rate, suggesting a direct route from the present stability theory to practical momentum-schedule design. Outlook. Our results suggest that learning rate, batch size, and momentum should be designed jointly rather than tuned in isolation. The stability conditions determine how aggressively momentum can be increased at a given learning rate and batch size, while the scaling laws quantify the resulting tradeoff between faster signal learning and stronger noise accumulation. Extending this framework to time-varying ηk , Bk , and ρk could therefore provide principled guidance for jointly scaling and scheduling these hyperparameters in large-scale training (Li et al., 2026a,b; Wang et al., 2025, 2026). We expect such a theory to provide a practical route toward designing more stable and data-efficient training schedules, rather than relying primarily on empirical tuning.

Use of AI We used AI tools for drafting and language editing throughout the manuscript. In Sections B–F, they were also used to accelerate algebraic calculations and asymptotic estimates. The proof strategies and mathematical arguments were developed by the authors, who verified all AI-assisted calculations and estimates.

Acknowledgement Lei Wu is supported in part by the National Natural Science Foundation of China (NSFC12522120, NSFC92470122, and NSFC12288101). Kairui Li is partially supported by the elite undergraduate training program of School of Mathematical Sciences in Peking University. We thank Sanxing Cao, Zilin Wang, and Jinbo Wang for helpful discussions.

References Alexander Atanasov, Jacob A. Zavatone-Veth, and Cengiz Pehlevan. Scaling and renormalization in high-dimensional regression. Journal of Statistical Mechanics: Theory and Experiment, 2026. (cited on page 4) Yasaman Bahri, Ethan Dyer, Jared Kaplan, Jaehoon Lee, and Utkarsh Sharma. Explaining neural scaling laws. Proceedings of the National Academy of Sciences, 121(27), 2024. doi: 10.1073/pnas.2311878121. (cited on pages 2 and 4) Blake Bordelon, Abdulkadir Canatar, and Cengiz Pehlevan. Spectrum dependent learning curves in kernel regression and wide neural networks. In International Conference on Machine Learning, pp. 1024–1034. PMLR, 2020. (cited on page 5) Blake Bordelon, Alexander Atanasov, and Cengiz Pehlevan. A dynamical model of neural scaling laws. In International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pp. 4345–4382, 2024. (cited on page 4) Blake Bordelon, Alexander Atanasov, and Cengiz Pehlevan. How feature learning can improve neural scaling laws. In International Conference on Learning Representations, 2025. (cited on pages 2 and 4)

Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al. Language models are few-shot learners. Advances in Neural Information Processing Systems, 33:1877–1901, 2020. (cited on page 2) 14

Andrea Caponnetto and Ernesto De Vito. Optimal rates for the regularized least-squares algorithm. Foundations of Computational Mathematics, 7:331–368, 2007. (cited on page 5) Andrea Caponnetto and Ernesto De Vito. Fast rates for regularized least-squares algorithm. 2005. (cited on page 5) Dennis Chemnitz and Maximilian Engel. Characterizing dynamical stability of stochastic gradient descent in overparameterized learning. Journal of Machine Learning Research, 26 (134):1–46, 2025. (cited on page 4) DeepSeek-AI, Anyi Xu, Bangcai Lin, Bing Xue, et al. DeepSeek-V4: Towards highly efficient million-token context intelligence. arXiv preprint arXiv:2606.19348, 2026. (cited on page 4) Shihong Ding, Haihan Zhang, Hanzhen Zhao, and Cong Fang. Scaling law for stochastic gradient descent in quadratically parameterized linear regression. arXiv preprint arXiv:2502.09106, 2025. (cited on page 4) Elvis Dohmatob, Yunzhen Feng, Pu Yang, Francois Charton, and Julia Kempe. A tale of tails: Model collapse as a change of scaling laws. In International Conference on Machine Learning, 2024. (cited on page 4) William Feller. An introduction to probability theory and its applications, volume 2. John Wiley & Sons, 1991. (cited on page 23) Damien Ferbach, Katie Everett, Gauthier Gidel, Elliot Paquette, and Courtney Paquette. Dimension-adapted momentum outscales sgd. Advances in Neural Information Processing Systems, 38:112780–112977, 2026a. (cited on page 4) Damien Ferbach, Courtney Paquette, Gauthier Gidel, Katie Everett, and Elliot Paquette. Logarithmic-time schedules for scaling language models with momentum. arXiv preprint arXiv:2602.05298, 2026b. (cited on page 13) Igor Gitman, Hunter Lang, Pengchuan Zhang, and Lin Xiao. Understanding the role of momentum in stochastic gradient methods. In Advances in Neural Information Processing Systems, volume 32, pp. 9633–9643, 2019. (cited on pages 2 and 4) Aaron Grattafiori, Abhimanyu Dubey, Abhinav Jauhri, et al. The LLaMA 3 herd of models. arXiv preprint arXiv:2407.21783, 2024. (cited on page 4) Kanan Gupta, Jonathan W. Siegel, and Stephan Wojtowytsch. Nesterov acceleration despite very noisy gradients. In Advances in Neural Information Processing Systems, volume 37, pp. 20694–20744, 2024. (cited on page 2) Joel Hestness, Sharan Narang, Newsha Ardalani, Gregory Diamos, Heewoo Jun, Hassan Kianinejad, Mostofa Ali Patwary, Yang Yang, and Yanqi Zhou. Deep learning scaling is predictable, empirically. arXiv preprint arXiv:1712.00409, 2017. (cited on page 4) Jordan Hoffmann, Sebastian Borgeaud, Arthur Mensch, Elena Buchatskaya, Trevor Cai, Eliza Rutherford, Diego de Las Casas, Lisa Anne Hendricks, Johannes Welbl, Aidan Clark, Tom Hennigan, Eric Noland, Katie Millican, George van den Driessche, Bogdan Damoc, Aurelia Guy, Simon Osindero, Karen Simonyan, Erich Elsen, Oriol Vinyals, Jack W. Rae, and Laurent Sifre. Training compute-optimal large language models. In Advances in Neural Information Processing Systems, volume 35, pp. 30016–30030, 2022. (cited on page 4) Marcus Hutter. Learning curve theory. arXiv preprint arXiv:2102.04074, 2021. (cited on page 4)

15

Ayush Jain, Andrea Montanari, and Eren Sasoglu. Scaling laws for learning with real and surrogate data. In Advances in Neural Information Processing Systems, volume 37, 2024. (cited on page 4) Keller Jordan, Yuchen Jin, Vlado Boza, Jiacheng You, Franz Cesista, Laker Newhouse, and Jeremy Bernstein. Muon: An optimizer for hidden layers in neural networks, 2024. URL https://kellerjordan. github. io/posts/muon, 6(3):4, 2024. (cited on page 2) Jared Kaplan, Sam McCandlish, Tom Henighan, Tom B. Brown, Benjamin Chess, Rewon Child, Scott Gray, Alec Radford, Jeffrey Wu, and Dario Amodei. Scaling laws for neural language models. arXiv preprint arXiv:2001.08361, 2020. (cited on page 4) Simran Kaur, Jeremy Cohen, and Zachary Chase Lipton. On the maximum Hessian eigenvalue and generalization. In Proceedings on “I Can’t Believe It’s Not Better! – Understanding Deep Learning Through Empirical Falsification” at NeurIPS 2022 Workshops, volume 187 of Proceedings of Machine Learning Research, pp. 51–65, 2023. (cited on page 4) Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. International Conference on Learning Representations, 2015. (cited on page 2)

In

Nikola B. Kovachki and Andrew M. Stuart. Continuous time analysis of momentum methods. Journal of Machine Learning Research, 22(17):1–40, 2021. (cited on page 4) Binghui Li, Fengling Chen, Zixun Huang, Lean Wang, and Lei Wu. Functional scaling laws in kernel regression: Loss dynamics and learning rate schedules. Advances in Neural Information Processing Systems, 38:101211–101269, 2026a. (cited on pages 2, 4, 8, 14, 33, and 46) Binghui Li, Zilin Wang, Fengling Chen, Shiyang Zhao, Ruiheng Zheng, and Lei Wu. Optimal learning-rate schedules under functional scaling laws: Power decay and warmup-stable-decay. arXiv preprint arXiv:2602.06797, 2026b. (cited on pages 2, 4, and 14) Licong Lin, Jingfeng Wu, Sham M. Kakade, Peter L. Bartlett, and Jason D. Lee. Scaling laws in linear regression: Compute, parameters, and data. In Advances in Neural Information Processing Systems, volume 37, 2024. (cited on pages 2 and 4) Licong Lin, Jingfeng Wu, and Peter L. Bartlett. Improved scaling laws in linear regression via data reuse. In Advances in Neural Information Processing Systems, volume 38, 2025. (cited on page 4)

Chaoyue Liu and Mikhail Belkin. Accelerating SGD with momentum for over-parameterized learning. arXiv preprint arXiv:1810.13395, 2018. (cited on page 2) Jingyuan Liu, Jianlin Su, Xingcheng Yao, Zhejun Jiang, Guokun Lai, Yulun Du, Yidao Qin, Weixin Xu, Enzhe Lu, Junjie Yan, Yanru Chen, Huabin Zheng, Yibo Liu, Shaowei Liu, Bohong Yin, Weiran He, Han Zhu, Yuzhi Wang, Jianzhou Wang, Mengnan Dong, Zheng Zhang, Yongsheng Kang, Hao Zhang, Xinran Xu, Yutao Zhang, Yuxin Wu, Xinyu Zhou, and Zhilin Yang. Muon is scalable for LLM training. arXiv preprint arXiv:2502.16982, 2025. (cited on page 2) Yanli Liu, Yuan Gao, and Wotao Yin. An improved analysis of stochastic gradient descent with momentum. In Advances in Neural Information Processing Systems, volume 33, 2020. (cited on page 4)

Chao Ma and Lexing Ying. On linear stability of SGD and input-smoothness of neural networks. In Advances in Neural Information Processing Systems, volume 34, 2021. (cited on page 4)

16

Alexander Maloney, Daniel A. Roberts, and James Sully. A solvable model of neural scaling laws. arXiv preprint arXiv:2210.16859, 2022. (cited on pages 4 and 5) Martin Marek, Sanae Lotfi, Aditya Somasundaram, Andrew Wilson, and Micah Goldblum. Small batch size training for language models: When vanilla SGD works, and why gradient accumulation is wasteful. Advances in Neural Information Processing Systems, 38: 148837–148862, 2026. (cited on page 2) Eric J. Michaud, Ziming Liu, Uzay Girit, and Max Tegmark. The quantization model of neural scaling. In Advances in Neural Information Processing Systems, volume 36, 2023. (cited on page 4)

Depen Morwani, Alexandru Meterez, Pranav Nair, and Sham Kakade. Compute efficiency and serial runtime tradeoffs for stochastic momentum methods. arXiv preprint arXiv:2606.19179, 2026. (cited on page 4) Rotem Mulayoff, Tomer Michaeli, and Daniel Soudry. The implicit bias of minima stability: A view from function space. In Advances in Neural Information Processing Systems, volume 34, 2021. (cited on page 4) Mor Shpigel Nacson, Rotem Mulayoff, Greg Ongie, Tomer Michaeli, and Daniel Soudry. The implicit bias of minima stability in multivariate shallow ReLU networks. In International Conference on Learning Representations, 2023. (cited on page 4) Yoonsoo Nam, Nayara Fonseca, Seok Hyeong Lee, and Ard Louis. An exactly solvable model for emergence and scaling laws in the multitask sparse parity problem. In Advances in Neural Information Processing Systems, volume 37, 2024. (cited on page 4) Kamil Nar and S. Shankar Sastry. Step size matters in deep learning. In Advances in Neural Information Processing Systems, volume 31, pp. 3436–3444, 2018. (cited on page 4) Yurii Nesterov. A method for solving the convex programming problem with convergence rate o (1/k 2 ). In Dokl Akad Nauk SSSR, volume 269, pp. 543, 1983. (cited on pages 2, 4, and 13) Matteo Pagliardini, Pierre Ablin, and David Grangier. The AdEMAMix optimizer: Better, faster, older. In International Conference on Learning Representations, volume 2025, pp. 64715–64757, 2025. (cited on page 13) Elliot Paquette, Courtney Paquette, Lechao Xiao, and Jeffrey Pennington. 4+3 phases of compute-optimal neural scaling laws. arXiv preprint arXiv:2405.15074, 2024. (cited on page 4) Boris T Polyak. Some methods of speeding up the convergence of iteration methods. USSR computational mathematics and mathematical physics, 4(5):1–17, 1964. (cited on page 2) Dan Qiao, Kaiqi Zhang, Esha Singh, Daniel Soudry, and Yu-Xiang Wang. Stable minima cannot overfit in univariate ReLU networks: Generalization by large step sizes. In Advances in Neural Information Processing Systems, volume 37, 2024. (cited on page 4) Andrei Semenov, Matteo Pagliardini, and Martin Jaggi. Benchmarking optimizers for large language model pretraining. arXiv preprint arXiv:2509.01440, 2025. (cited on page 2) Ishaan Shah, Anthony M Polloreno, Karl Stratos, Philip Monk, Adarsh Chaluvaraju, Andrew Hojel, Andrew Ma, Anil Thomas, Ashish Tanwer, Darsh J Shah, et al. Practical efficiency of Muon for pretraining. arXiv preprint arXiv:2505.02222, 2025. (cited on page 2) Utkarsh Sharma and Jared Kaplan. Scaling laws from the data manifold dimension. Journal of Machine Learning Research, 23(9):1–34, 2022. (cited on page 4) 17

Bin Shi, Simon S. Du, Michael I. Jordan, and Weijie J. Su. Understanding the acceleration phenomenon via high-resolution differential equations. Mathematical Programming, 195(1–2): 79–148, 2022. doi: 10.1007/s10107-021-01681-8. (cited on page 4) John Sous and Michael Winer. Asymmetric scaling laws from sparse features. arXiv preprint arXiv:2605.23591, 2026. (cited on page 4) Stefano Spigler, Mario Geiger, and Matthieu Wyart. Asymptotic learning curves of kernel methods: empirical data versus teacher–student paradigm. Journal of Statistical Mechanics: Theory and Experiment, 2020(12):124001, 2020. (cited on page 5) Weijie Su, Stephen Boyd, and Emmanuel J. Candès. A differential equation for modeling Nesterov’s accelerated gradient method: Theory and insights. Journal of Machine Learning Research, 17(153):1–43, 2016. (cited on pages 2 and 4) Maksim Velikanov, Denis Kuznedelev, and Dmitry Yarotsky. A view of mini-batch SGD via generating functions: Conditions of convergence, phase transitions, benefit from negative momenta. In International Conference on Learning Representations, 2023. (cited on pages 2 and 4)

Mathukumalli Vidyasagar. Convergence of momentum-based optimization algorithms with time-varying parameters. arXiv preprint arXiv:2506.11904, 2025. (cited on page 4) Jinbo Wang, Mingze Wang, Zhanpeng Zhou, Junchi Yan, Lei Wu, et al. The sharpness disparity principle in transformers for accelerating language model pre-training. arXiv preprint arXiv:2502.19002, 2025. (cited on page 14) Jinbo Wang, Binghui Li, Zhanpeng Zhou, Mingze Wang, Jiaqi Zhang, Xunliang Cai, and Lei Wu. Fast catch-up, late switching: Optimal batch size scheduling via functional scaling laws. In International Conference on Learning Representations, volume 2026, pp. 7310–7340, 2026. (cited on pages 2, 4, and 14) Mingze Wang and Lei Wu. A theoretical analysis of noise geometry in stochastic gradient descent. arXiv preprint arXiv:2310.00692, 2023. (cited on pages 4 and 6) Runzhe Wang, Sadhika Malladi, Tianhao Wang, Kaifeng Lyu, and Zhiyuan Li. The marginal value of momentum for small learning rate SGD. In International Conference on Learning Representations, 2024. (cited on page 2) Alexander Wei, Wei Hu, and Jacob Steinhardt. More than a toy: Random matrix models predict how real-world neural representations generalize. In International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pp. 23549–23588, 2022. (cited on page 4) Andre Wibisono, Ashia C. Wilson, and Michael I. Jordan. A variational perspective on accelerated methods in optimization. Proceedings of the National Academy of Sciences, 113 (47):E7351–E7358, 2016. doi: 10.1073/pnas.1614734113. (cited on pages 2 and 4) Ashia C. Wilson, Ben Recht, and Michael I. Jordan. A Lyapunov analysis of accelerated methods in optimization. Journal of Machine Learning Research, 22(113):1–34, 2021. (cited on pages 2 and 4)

Lei Wu and Weijie J. Su. The implicit regularization of dynamical stability in stochastic gradient descent. In International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pp. 37656–37684, 2023. (cited on page 4)

18

Lei Wu, Chao Ma, and Weinan E. How SGD selects the global minima in over-parameterized learning: A dynamical stability perspective. In Advances in Neural Information Processing Systems, volume 31, pp. 8279–8288, 2018. (cited on page 4) Lei Wu, Mingze Wang, and Weijie J. Su. The alignment property of SGD noise and how it helps select flat minima: A stability analysis. In Advances in Neural Information Processing Systems, volume 35, pp. 4680–4693, 2022. (cited on pages 4 and 6) Tingkai Yan, Haodong Wen, Binghui Li, Kairong Luo, Wenguang Chen, and Kaifeng Lyu. Larger datasets can be repeated more: A theoretical analysis of multi-epoch scaling in linear regression. In International Conference on Learning Representations, 2026. (cited on page 4) Dmitry Yarotsky. Corner gradient descent. In International Conference on Learning Representations, volume 2026, pp. 156292–156327, 2026. (cited on page 13) Dmitry Yarotsky and Maksim Velikanov. SGD with memory: fundamental properties and stochastic acceleration. In International Conference on Learning Representations, volume 2025, pp. 79673–79715, 2025. (cited on page 13) Kun Yuan, Bicheng Ying, and Ali H. Sayed. On the influence of momentum acceleration on online learning. Journal of Machine Learning Research, 17(192):1–66, 2016. (cited on page 4) Haihan Zhang, Yuanshi Liu, Qianwen Chen, and Cong Fang. The optimality of (accelerated) SGD for high-dimensional quadratic optimization. arXiv preprint arXiv:2409.09745, 2024. (cited on page 4) Hanlin Zhang, Depen Morwani, Nikhil Vyas, Jingfeng Wu, Difan Zou, Udaya Ghai, Dean Foster, and Sham M. Kakade. How does critical batch size scale in pre-training? In International Conference on Learning Representations, 2025. (cited on page 4) Zhanxing Zhu, Jingfeng Wu, Bing Yu, Lei Wu, and Jinwen Ma. The anisotropic noise in stochastic gradient descent: Its behavior of escaping from sharp minima and regularization effects. In International Conference on Machine Learning, pp. 7654–7663, 2019. (cited on pages 4 and 6)

19

Appendix Table of Contents A Technical preliminaries

21

A.1 Equivalence to kernel regression . . . . . . . . . . . . . . . . . . . . . . . . . .

21

A.2 Rotational equivariance of SGD and momentum methods . . . . . . . . . . .

21

A.3 Gaussian fourth-moment identities . . . . . . . . . . . . . . . . . . . . . . . .

22

A.4 Noise–curvature alignment . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

B Risk stability analysis

23

B.1 Renewal formulation of risk stability . . . . . . . . . . . . . . . . . . . . . . .

23

B.2 Risk stability of SGD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

B.3 Risk stability of Polyak . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

B.4 Risk stability of Nesterov . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

C SGD scaling law

33

D Polyak scaling law D.1 Signal–noise decomposition and risk recursion . . . . . . . . . . . . . . . . . .

33 33

D.2 Overdamped and underdamped modes decomposition . . . . . . . . . . . . .

36

D.3 Label noise calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

D.4 Uniform boundedness of expected risk . . . . . . . . . . . . . . . . . . . . . .

39

D.5 Risk dynamics in the fully overdamped regime . . . . . . . . . . . . . . . . .

40

D.6 Risk dynamics in the mixed-damping regime

43

. . . . . . . . . . . . . . . . . .

E Nesterov scaling law

46

E.1 Signal–noise decomposition and risk recursion . . . . . . . . . . . . . . . . . .

F

46

E.2 Overdamped and underdamped modes decomposition . . . . . . . . . . . . .

50

E.3 Label noise calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

E.4 Uniform boundedness of expected risk . . . . . . . . . . . . . . . . . . . . . .

53

E.5 Risk dynamics in the fully overdamped regime . . . . . . . . . . . . . . . . .

54

E.6 Risk dynamics in the mixed-damping regime

57

. . . . . . . . . . . . . . . . . .

Optimal data scaling across batch-size regimes

62

F.1 SGD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

62

F.2 Polyak . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

63

F.3 Nesterov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

64

G Experimental details and additional results

66

G.1 Risk stability boundaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

66

G.2 Risk-dynamics scaling laws . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67

G.3 Batch-size phase diagram at a fixed data budget . . . . . . . . . . . . . . . .

68

G.4 Data-scaling exponents across batch-size regimes . . . . . . . . . . . . . . . .

68

G.5 Delayed instability under increasing momentum . . . . . . . . . . . . . . . . .

69

20

A

Technical preliminaries

A.1

Equivalence to kernel regression

We briefly relate the PLK regression in Section 3.1 to standard kernel regression. Given a positive semidefinite kernel K : X × X → R, kernel regression learns over its associated reproducing kernel Hilbert space (RKHS) HK . Given an input distribution D, under mild conditions, Mercer’s theorem ensures that the kernel admits the following eigendecomposition ′

K(x, x ) =

∞ X

λj ej (x)ej (x′ ),

j=1

where {ej }j⩾1 are orthonormal in L2 (D). The standard capacity condition assumes λj ≂ j −β . Target regularity is commonly described through the interpolation spaces   ∞ ∞ 2 X  X aj r , HK = aj ej : < ∞   λrj j=1

j=1

r . with a source condition of the form f ∗ ∈ HK Our formulation is the corresponding feature-map representation. Under Assumption 3.1, 1/2 define ej := ϕj /λj . Then {ej }j⩾1 is orthonormal in L2 (D), and ∞ X

′

Kϕ (x, x ) =

′

ϕj (x)ϕj (x ) =

j=1

∞ X

λj ej (x)ej (x′ ).

j=1

Thus, Assumption 3.3, λj ≂ j −β , is precisely the standard capacity condition. For the target function, f∗ =

∞ X

θj∗ ϕj =

j=1

∞ X

1/2

a∗j ej ,

a∗j := λj θj∗ .

j=1

Our task-difficulty assumption gives |a∗j |2 ≂ j −1 λsj . Hence, for any δ ∈ (0, s), ∞ X |a∗j |2

λs−δ j=1 j

≂

∞ X

j −1 λδj ≂

j=1

∞ X

j −1−βδ < ∞.

j=1

s−δ Therefore, f ∗ ∈ HK for every δ ∈ (0, s). Thus, β corresponds to the standard capacity ϕ exponent, while s characterizes the source regularity of the target.

A.2

Rotational equivariance of SGD and momentum methods

We show that Assumption 3.1 is without loss of generality. Since H = E[ϕ(x)⊗ϕ(x)] is positive, self-adjoint, and trace class, it admits a spectral decomposition H = U ΛU ∗ ,

Λ = diag(λ1 , λ2 , . . .),

for some unitary operator U . Define the rotated variables θ ′ = U ∗ θ,

θ ∗′ = U ∗ θ ∗ ,

ϕ′ (x) = U ∗ ϕ(x).

Unitary invariance of the inner product gives f ′ (x; θ ′ ) = ⟨ϕ′ (x), θ ′ ⟩ = ⟨ϕ(x), θ⟩ = f (x; θ), 21

and hence, for every mini-batch Bk , g′ (θ ′ ; Bk ) = U ∗ g(θ; Bk ). It follows immediately that the SGD and Polyak recursions are preserved under θk′ = U ∗ θk . For Nesterov, the look-ahead point transforms as ′ θkla′ = U ∗ θkla = θk′ + ρ(θk′ − θk−1 ),

so its recursion is preserved as well. Finally, H′ = E[ϕ′ (x) ⊗ ϕ′ (x)] = U ∗ HU = Λ. Thus, SGD, Polyak, and Nesterov are equivariant under unitary changes of basis, while their predictions and risks are unchanged. We may therefore work in the eigenbasis of H, where the covariance operator is diagonal, without loss of generality.

A.3

Gaussian fourth-moment identities

We record several Gaussian fourth-moment identities used repeatedly below. Lemma A.1 (Gaussian fourth-moment identities). Let x ∼ N (0, H) in a separable Hilbert space H, where H is a positive trace-class covariance operator. For any u, v ∈ H,   (17) E ⟨x, u⟩2 ⟨x, v⟩2 = ⟨u, Hu⟩⟨v, Hv⟩ + 2⟨u, Hv⟩2 . Equivalently, for any self-adjoint trace-class operator A, E[(x ⊗ x)A(x ⊗ x)] = 2HAH + HTr(HA).

(18)

Proof. The directional identity follows directly from Isserlis’ formula. For any u, v ∈ H, the random variables ⟨x, u⟩ and ⟨x, v⟩ are jointly Gaussian, and hence   E ⟨x, u⟩2 ⟨x, v⟩2 = ⟨u, Hu⟩⟨v, Hv⟩ + 2⟨u, Hv⟩2 . P We next prove the operator identity. Let A = j aj ej ⊗ ej be the spectral decomposition of the self-adjoint trace-class operator A. For arbitrary u, v ∈ H, Isserlis’ formula gives ⟨u, E[(x ⊗ x)A(x ⊗ x)] v⟩ X   = aj E ⟨u, x⟩⟨ej , x⟩2 ⟨x, v⟩ j

=2

X

aj ⟨u, Hej ⟩⟨ej , Hv⟩ + ⟨u, Hv⟩

j

X

aj ⟨ej , Hej ⟩

j

= 2⟨u, HAHv⟩ + ⟨u, Hv⟩Tr(HA). Since this holds for all u, v ∈ H, E[(x ⊗ x)A(x ⊗ x)] = 2HAH + HTr(HA).

Corollary A.2 (Empirical covariance fluctuation). Let b := 1 H B

B X

xb ⊗ x b ,

b − H, Σ := H

i.i.d.

xb ∼ N (0, H).

b=1

Then, for any self-adjoint trace-class operator A, E[ΣAΣ] =

1 (HAH + HTr(HA)) . B 22

(19)

Proof. Write Σ = B −1

PB

b=1 Zb , where Zb := xb ⊗ xb − H are independent and centered. Hence

E[ΣAΣ] =

1 1 E[ZAZ] = (E[(x ⊗ x)A(x ⊗ x)] − HAH) , B B

and the result follows from Lemma A.1.

A.4

Noise–curvature alignment

P We prove Proposition 3.4. For a mini-batch of size B, write ξk (θ) = B1 B i=1 ζk,i (θ), where ζk,i (θ) are conditionally independent, centered single-sample gradient noises. Hence, for any unit vector ν ∈ H,    1  E |⟨ξk (θ), ν⟩|2 θ = E |⟨ζ(θ), ν⟩|2 θ . B ∗ Let δ := θ − θ
and ϕ := ϕ(x). Since y = ⟨ϕ, θ ∗ ⟩ + ϵ, the single-sample gradient is gx,y (θ) = ⟨ϕ, δ⟩ − ϵ ϕ, while ∇R(θ) = Hδ. Therefore, ζ(θ) = ⟨ϕ, δ⟩ϕ − Hδ − ϵϕ. Using the independence of ϵ and ϕ, E[ϵ] = 0, and E[ϕϕ] = H, we obtain     E |⟨ζ(θ), ν⟩|2 θ = E ⟨ϕ, δ⟩2 ⟨ϕ, ν⟩2 − ⟨δ, Hν⟩2 + σ 2 ⟨ν, Hν⟩. Applying Lemma A.1 gives   E ⟨ϕ, δ⟩2 ⟨ϕ, ν⟩2 = ⟨δ, Hδ⟩⟨ν, Hν⟩ + 2⟨δ, Hν⟩2 . Thus,  
E |⟨ζ(θ), ν⟩|2 θ = σ 2 + ⟨δ, Hδ⟩ ⟨ν, Hν⟩ + ⟨δ, Hν⟩2 .

(20)

The last term is nonnegative and, by Cauchy–Schwarz in the H-inner product, ⟨δ, Hν⟩2 ⩽ ⟨δ, Hδ⟩⟨ν, Hν⟩. Since E(θ) = 21 ⟨δ, Hδ⟩, Eq. (20) implies 
 E |⟨ζ(θ), ν⟩|2 θ ≂ E(θ) + σ 2 ⟨ν, Hν⟩. Combining this with the 1/B variance reduction from mini-batching and setting θ = θk yields  
1 E |⟨ξk (θk ), ν⟩|2 θk ≂ E(θk ) + σ 2 ⟨ν, Hν⟩, B which proves Proposition 3.4.

B

Risk stability analysis

B.1

Renewal formulation of risk stability

We first recall a general boundedness criterion for renewal equations (Feller, 1991). Consider a nonnegative sequence (xk )k⩾0 satisfying xk = ak +

k−1 X

Kk−1−j xj ,

j=0

where ak ⩾ 0 is the input and K = (Kk )k⩾0 is the renewal kernel. Its total mass ∥K∥1 := P K k⩾0 k measures the overall feedback strength of the recursion. 23

Lemma B.1 (Boundedness of a renewal equation). Suppose 0 < a := inf k⩾0 ak ⩽ supk⩾0 ak =: ā < ∞. Then sup xk < ∞ ⇐⇒ ∥K∥1 < 1. k⩾0

Proof. Suppose first that ∥K∥1 < 1, and let Mk := max0⩽j⩽k xj . Then xk ⩽ ā + ∥K∥1 Mk−1 , which implies   ā sup xk ⩽ max x0 , < ∞. 1 − ∥K∥1 k⩾0 Conversely, suppose that (xk ) is bounded and let L := lim inf k→∞ xk . For every fixed N , nonnegativity gives N X xk ⩾ a + Km xk−1−m . m=0

Taking lim inf k→∞ yields L ⩾ a + L P K m⩾0 m ⩾ 1.

PN

m=0 Km . Letting N → ∞, boundedness is impossible if

For each method A ∈ {sgd, polyak, nesterov}, we will derive a nonnegative risk moment (A) ℓk , equal to 2E[E(θk )] for SGD and Polyak and comparable to it for Nesterov, satisfying (A)

ℓk

(A)

= hk

+

k−1 X

(A)

(A)

Kk−1−j ℓj

+ σ2 .

(21)

j=0 (A)

(A)

Setting xk = ℓk + σ 2 and ak = hk + σ 2 reduces (21) to the renewal equation above. Hence, (A) once the coordinate-wise dynamics are stable and hk is bounded, Lemma B.1 gives the global feedback condition X (A) κA (η, ρ, B) := ∥K (A) ∥1 = Kk < 1. k⩾0

The kernel K (A) describes how perturbations propagate through the optimizer dynamics and depends on the method and hyperparameters (η, ρ, B), but not on the label-noise variance σ 2 . Additive label noise enters as the constant external forcing σ 2 , whereas multiplicative noise feeds the current risk back through K (A) . The following subsections derive K (A) , compute κA , and combine κA < 1 with the corresponding coordinate-wise stability condition to obtain the critical learning rates.

B.2

Risk stability of SGD

The SGD error recursion, including label noise, is b k )ek + η ζbk . ek+1 = (I − η H

(22)

Theorem B.2 (Risk stability of SGD). Let the Hessian be a positive trace-class operator with spectrum H = diag(λ1 , λ2 , . . .), λ1 = λmax ⩾ λ2 ⩾ · · · > 0, and suppose σ > 0. Then SGD is risk stable if and only if η<

2B (B + 1)λmax

(23)

X

ηλi < 1. 2B − (B + 1)ηλi

(24)

and κsgd (η) :=

i⩾1

24

Proof. Define the covariance operator by Ck := E[ek ⊗ ek ]. We work with its spectral coordinates. For every finite-risk initialization, the relevant risk moment is X   ℓk := Tr(HCk ) = λi E |⟨ei , ek ⟩|2 = 2E[E(θk )]. i⩾1

In particular, ℓ0 < ∞. The covariance identities below may be justified first on finite spectral truncations and then passed to the limit by monotone convergence. Expanding (22), using the vanishing cross terms, and applying Lemma A.1 gives Ck+1 = (I − ηH)Ck (I − ηH) +

η2 HCk H + HTr(HCk ) + σ 2 H . B

(25)

Define the coordinate variance by ui,k := ⟨ei , Ck ei ⟩. Taking the diagonal entries of (25) yields ui,k+1 = ai ui,k +

η2 λ i ℓk + σ 2 , B

ai := (1 − ηλi )2 +

η2 2 λ . B i

(26)

We first determine the local condition. Since     1 1 − ai = ηλi 2 − 1 + ηλi , B the coordinatewise contraction condition ai < 1

for every i

is equivalent to (23). If this condition fails, then (26) contains an eigendirection satisfying ui,k+1 ⩾ ui,k +

η2σ2 λi , B

and therefore the risk is unbounded. Hence the local condition is necessary. Assume henceforth that (23) holds. Iterating (26), multiplying the i-th equation by λi , and summing over all eigendirections gives ℓk = hk +

k−1 X

Kk−1−s ℓs + σ 2 ,

s=0

where hk :=

X

λi aki ui,0 ,

Km :=

i⩾1

η2 X 2 m λi ai . B i⩾1

The local condition gives 0 ⩽ ai < 1,

0 ⩽ hk ⩽ ℓ0 .

Moreover, X

Km =

m⩾0

X η 2 X λ2i ηλi = = κsgd (η). B 1 − ai 2B − (B + 1)ηλi i⩾1

i⩾1

Lemma B.1 therefore gives sup ℓk < ∞

⇐⇒

k⩾0

which completes the proof.

25

κsgd (η) < 1,

(27)

Corollary B.3 (Critical learning rate of SGD). There exists a unique critical learning rate satisfying 2B crit 0 < ηsgd < (B + 1)λmax and

crit λ ηsgd i

X

crit λ 2B − (B + 1)ηsgd i i⩾1

= 1.

(28)

The SGD iteration is risk stable if and only if crit 0 < η < ηsgd .

Moreover, 2B (B + 1)λmax + Tr(H)

crit ⩽ ηsgd ⩽ min



2B 2B , (B + 1)λmax Tr(H)

 .

(29)

Consequently, under the fixed-exponent power-law assumptions λi ≂ i−β , β > 1, we have crit ηsgd ≂ 1.

Proof. The renewal mass κsgd (η) is continuous and strictly increasing on the local stability interval. It vanishes at η = 0 and diverges as η approaches the local boundary, because the denominator of the leading eigendirection tends to zero. Hence there is a unique solution of κsgd (η) = 1, which gives the stated critical learning rate and stability interval. At the critical point, every denominator in (28) is at most 2B. Therefore 1 ⩾ which yields 2B crit ⩽ ηsgd . Tr(H)

crit ηsgd 2B Tr(H),

The local condition gives the other upper bound in (29). Conversely, every denominator is at least crit λmax . 2B − (B + 1)ηsgd Using the critical equation then gives 1⩽

crit Tr(H) ηsgd crit λ 2B − (B + 1)ηsgd max

,

and rearranging yields the lower bound in (29). Finally, under λi ≂ i−β with β > 1, we have λmax ≂ 1 and Tr(H) ≂ 1. The two-sided crit ≂ 1. bounds therefore imply ηsgd

B.3

Risk stability of Polyak

The Polyak error recursion, including label noise, is b k )ek + ρ(ek − ek−1 ) + η ζbk . ek+1 = (I − η H

(30)

Writing D = (1 + ρ)I − ηH, Equation (30) becomes ek+1 = Dek − ρek−1 − η(Hk − H)ek + η ζbk . Theorem B.4 (Risk stability of Polyak). Let the Hessian be a positive trace-class operator and assume H = diag(λ1 , λ2 , . . .),

λ1 = λmax ⩾ λ2 ⩾ · · · > 0, 26

0 ⩽ ρ < 1,

σ > 0.

Then Polyak momentum is risk stable if and only if η< and κpolyak (η, ρ) :=

2B(1 − ρ2 ) λmax [B(1 − ρ) + (1 + ρ)]

(31)

ηλi (1 + ρ) < 1. 2 2B(1 − ρ ) − ηλi [B(1 − ρ) + (1 + ρ)] i⩾1

X

(32)

Proof. Define the coordinate second moments Uk := E[ek ⊗ ek ],

Pk := E[ek ⊗ ek−1 ],

Qk := E[ek−1 ⊗ ek ], Rk := E[ek−1 ⊗ ek−1 ]. The risk moment relevant to risk stability is ℓk := Tr(HUk ) = 2E[E(θk )]. Expanding (30), using the vanishing cross terms, and applying Lemma A.1 gives  Uk+1 = DUk D − ρDPk − ρQk D + ρ2 Rk     
 η2   HUk H + HTr(HUk ) + σ 2 H , +  B Pk+1 = DUk − ρQk ,      Qk+1 = Uk D − ρPk ,     Rk+1 = Uk .

(33)

In the eigenbasis of the Hessian, define ui,k := ⟨ei , Uk ei ⟩,

vi,k := ⟨ei , Pk ei ⟩,

ri,k := ⟨ei , Rk ei ⟩.

The diagonal entries of the two cross-covariance operators coincide because Qk = P∗k . Set di := 1 + ρ − ηλi , and

Xi,k := (ui,k , vi,k , ri,k )⊤ , 

η 2 di +

Mi :=  



2

B

b := (1, 0, 0)⊤ ,

λ2i

di 1

−2ρdi −ρ 0

ρ2

  0 . 0

Taking the diagonal entries of (33) yields Xi,k+1 = Mi Xi,k +

η2 λi b ℓk + σ 2 . B

(34)

After removing the global risk-coupling term, the characteristic polynomial of the isolated block is     η2 2 η2 2 3 2 2 2 pi (µ) = µ + ρ − di − λi µ + ρ di − λi − ρ µ − ρ3 . B B The Jury stability criterion shows that all roots lie strictly inside the unit disk if and only if 2B(1 − ρ2 ) − ηλi [B(1 − ρ) + (1 + ρ)] > 0. Requiring this inequality in every eigendirection gives (31). If the local condition fails, the isolated second-moment system has a mode on or outside the unit circle. Because the label 27

noise forcing enters its current-error coordinate directly and has strictly positive variance, the corresponding coordinate risk is unbounded. Thus, the local condition is necessary. Assume henceforth that the local condition holds, and define gi,m := b⊤ Mm i b. This quantity is nonnegative. Iterating (34), taking the first component, multiplying by the corresponding eigenvalue, and summing over all eigendirections gives ℓk = hk +

k−1 X

Kk−1−s ℓs + σ 2 ,

(35)

s=0

where hk :=

X

λi b⊤ Mki Xi,0 ,

Km :=

η2 X 2 λi gi,m . B i⩾1

i⩾1

The strict local condition gives the uniform power bound supi,k ∥Mki ∥ < ∞. Positivity of the u +r initial coordinate covariance also gives |vi,0 | ⩽ i,0 2 i,0 . Together, these estimates imply X 0 ⩽ hk ⩽ C λi (ui,0 + ri,0 ) < ∞. i⩾1

A direct calculation gives
2 (1 + ρ) B 1 + ρ, −2ρd , ρ i . b⊤ (I − Mi )−1 = ηλi {2B(1 − ρ2 ) − ηλi [B(1 − ρ) + (1 + ρ)]}

(36)

Consequently, X

Km =

m⩾0

η2 X 2 ⊤ λi b (I − Mi )−1 b = κpolyak (η, ρ). B i⩾1

Lemma B.1 therefore gives sup ℓk < ∞

⇐⇒

κpolyak (η, ρ) < 1,

k⩾0

which completes the proof. Corollary B.5 (Critical learning rate of Polyak). There exists a unique critical learning rate satisfying 2B(1 − ρ2 ) crit < 0 < ηpolyak λmax [B(1 − ρ) + (1 + ρ)] and the critical equation X

crit λ (1 + ρ) ηpolyak i

i⩾1

crit λ [B(1 − ρ) + (1 + ρ)] 2B(1 − ρ2 ) − ηpolyak i

= 1.

(37)

The Polyak iteration is risk stable if and only if crit 0 < η < ηpolyak .

Moreover, 2B(1 − ρ2 ) crit ⩽ ηpolyak , λmax [B(1 − ρ) + (1 + ρ)] + (1 + ρ)Tr(H)   2B(1 − ρ2 ) 2B(1 − ρ) crit ηpolyak ⩽ min , . λmax [B(1 − ρ) + (1 + ρ)] Tr(H) Under the power-law assumption λi ≂ i−β , β > 1, when B ≫ 1, 1 − ρ ≪ 1, we have crit ηpolyak ≂ min{1, B(1 − ρ)}.

28

(38)

Proof. The renewal mass κpolyak (η, ρ) is continuous and strictly increasing throughout the local stability interval. It vanishes at the origin and diverges at the local boundary. Theorem B.4 therefore gives a unique critical learning rate and the stated stability interval. At the critical point, every denominator in (37) is at most 2B(1 − ρ2 ) and hence crit ηpolyak ⩽

2B(1 − ρ) , Tr(H)

while the local condition gives the other upper bound in (38). Conversely, every denominator crit λ is at least 2B(1 − ρ2 ) − ηpolyak max [B(1 − ρ) + (1 + ρ)] . Substitution into the critical equation gives crit (1 + ρ)Tr(H) ηpolyak . 1⩽ crit λ 2B(1 − ρ2 ) − ηpolyak max [B(1 − ρ) + (1 + ρ)] Rearranging this inequality gives the lower bound in (38). Finally, the power-law spectrum and the high-momentum regime give λmax ≂ 1,

Tr(H) ≂ 1,

1 − ρ2 ≂ 1 − ρ.

Consequently, both sides of (38) scale as B(1 − ρ) ≂ min{1, B(1 − ρ)}. B(1 − ρ) + 1

B.4

Risk stability of Nesterov

Define the look-ahead error ela k := (1 + ρ)ek − ρek−1 . The Nesterov error recursion, including label noise, is b k )ela + η ζbk . ek+1 = (I − η H k

(39)

With the deterministic transition D = I − ηH, the recursion becomes la b ek+1 = Dela k − η(Hk − H)ek + η ζk .

Theorem B.6 (Risk stability of Nesterov). Let the Hessian be a positive trace-class operator and assume H = diag(λ1 , λ2 , . . .),

λ1 = λmax ⩾ λ2 ⩾ · · · > 0,

0 ⩽ ρ < 1,

σ > 0.

Define   ∆nesterov (x) := 2B(1 − ρ2 ) + 4Bρ2 + (B − 1)ρ − (B + 1) x − (B + 1)ρ(1 + 2ρ)x2 . Then Nesterov momentum is risk stable if and only if ∆nesterov (ηλmax ) > 0 and κnesterov (η, ρ) :=

X ηλi [1 + ρ + ρ(1 + 2ρ)ηλi ] i⩾1

∆nesterov (ηλi )

29

(40)

< 1.

(41)

Proof. Define la Vk := E[ela k ⊗ ek ],

Uk := E[ek ⊗ ek ],

Qk := E[ek−1 ⊗ ek ],

Pk := E[ek ⊗ ek−1 ],

Rk := E[ek−1 ⊗ ek−1 ].

Expanding (39), using the vanishing cross terms, and applying Lemma A.1 gives 
η2    HVk H + HTr(HVk ) + σ 2 H , U = DV D + k+1 k   B     Vk = (1 + ρ)2 Uk − ρ(1 + ρ)(Pk + Qk ) + ρ2 Rk ,
P = D (1 + ρ)U − ρQ ,  k+1 k k 
   Qk+1 = (1 + ρ)Uk − ρPk D,     Rk+1 = Uk .

(42)

The renewal equation below naturally controls the look-ahead loss moment la 2 ℓla k := Tr(HVk ) = E∥ek ∥H .

This is equivalent to controlling the risk at the actual iterate. Define its loss moment by ℓk := Tr(HUk ) = 2E[E(θk )]. Then 2 2 ℓla k ⩽ 2(1 + ρ) ℓk + 2ρ ℓk−1 .

(43)

Conversely, since D commutes with H and H is trace class, the first equation in (42) gives
2 ℓk+1 ⩽ Cη,ρ,H,B ℓla (44) k +σ . Therefore sup ℓla k <∞

⇐⇒

sup E[E(θk )] < ∞. k

k

Let ui,k := ⟨ei , Uk ei ⟩,

pi,k := ⟨ei , Pk ei ⟩,

ri,k := ⟨ei , Rk ei ⟩,

and vi,k := ⟨ei , Vk ei ⟩. The cross-covariance operators satisfy Qk = P∗k . Therefore, vi,k = (1 + ρ)2 ui,k − 2ρ(1 + ρ)pi,k + ρ2 ri,k . Set di := 1 − ηλi ,

ai := d2i +

b := (1, 0, 0)⊤ , and

η2 2 λ , B i

Xi,k := (ui,k , pi,k , ri,k )⊤ ,

c := (1 + ρ)2 , −2ρ(1 + ρ), ρ2

  ai (1 + ρ)2 −2ai ρ(1 + ρ) ai ρ2 −di ρ 0 . Mi :=  di (1 + ρ) 1 0 0

Then vi,k = c⊤ Xi,k ,

ℓla k =

X i⩾1

30

λi c⊤ Xi,k ,

⊤

,

and the diagonal part of (42) becomes
η2 2 λi b ℓla k +σ . B

Xi,k+1 = Mi Xi,k +

(45)

The characteristic polynomial of the isolated block is     pi (µ) = µ3 − ai (1 + ρ)2 − di ρ µ2 + ai ρ di (1 + ρ)2 − ρ µ − ai di ρ3 . The Jury stability criterion shows that all roots lie strictly inside the unit disk if and only if ∆nesterov (ηλi ) > 0. The polynomial is a concave quadratic satisfying ∆nesterov (0) = 2B(1 − ρ2 ) > 0, and it has exactly one positive root. Imposing the coordinate condition in every eigendirection is therefore equivalent to (40). If it fails, the nondegenerate label-noise forcing excites a nonstable isolated mode, and the corresponding coordinate risk is unbounded. Assume henceforth that the local condition holds. Iterating (45), multiplying the i-th coordinate contribution by λi , and summing over all eigendirections gives ℓla k = hk +

k−1 X

2 Kk−1−s ℓla s +σ ,

(46)

s=0

where hk :=

X

λi c⊤ Mki Xi,0 ,

Km :=

η2 X 2 ⊤ m λi c Mi b. B i⩾1

i⩾1

Both sequences are nonnegative: their summands are look-ahead variances generated by isolated coordinate systems. The local blocks are uniformly power bounded, and positivity of the initial coordinate covariance gives X 0 ⩽ hk ⩽ C λi (ui,0 + ri,0 ) < ∞. i⩾1

A direct calculation gives ⊤

B 1 + ρ + ρ(1 + 2ρ)ηλi , −2ρ(1 + ρ), ρ2 [1 + ρ(1 − ηλi )] = . ηλi ∆nesterov (ηλi )

−1

c (I − Mi ) Its first entry yields

X

Km =

m⩾0

η2 X 2 ⊤ λi c (I − Mi )−1 b = κnesterov (η, ρ). B i⩾1

Lemma B.1, together with (43)– (44), gives sup E[E(θk )] < ∞

⇐⇒

k⩾0

which completes the proof.

31

κnesterov (η, ρ) < 1,

(47)

Corollary B.7 (Critical learning rate of Nesterov). There exists a unique critical learning rate satisfying
crit κnesterov ηnesterov ,ρ = 1 before the local stability boundary. The Nesterov iteration is risk stable if and only if crit 0 < η < ηnesterov .

Under the power-law and high-momentum assumptions λi ≂ i−β , β > 1, B ≫ 1, 1 − ρ ≪ 1, crit ηnesterov ≂ min{1, B β (1 − ρ)}.

(48)

Proof. First consider the local condition. The polynomial is concave and satisfies ∆nesterov (0) > 0. It has a unique positive root, denoted by x+ (B, ρ), and hence the local condition is equivalent to ηλmax < x+ (B, ρ). In the high-momentum regime, ρ → 1− ,

B ≫ 1, the root remains of constant order. Indeed, x+ (B, ρ) ≂ 1,

x+ (B, 1) =

4B − 2 4 −→ . 3(B + 1) 3

Together with the power-law normalization λmax ≂ 1, this gives the order-one local ceiling η ≲ 1. We next consider the global renewal mass κnesterov (η, ρ). Each summand is continuous and strictly increasing up to the local boundary, and the leading eigendirection diverges at that boundary. Hence Theorem B.6 gives a unique solution of the critical equation and the stated stability interval. On every fixed subinterval strictly inside the local stability region, κnesterov (η, ρ) ≂

1 X ηλi . B 1 − ρ + ηλi i⩾1

For the power-law spectrum, the spectral sum satisfies  η   1 − ρ, X ηλi 1/β ≂  η 1 − ρ + ηλi   i⩾1  , 1−ρ

η ≲ 1 − ρ, η ≳ 1 − ρ.

η η At a global critical point, the first regime would require 1−ρ ≂ B, contradicting 1−ρ ≲ 1 when the batch size is large. Therefore the global critical point lies in the second regime, where  1/β η ≂ B. 1−ρ

Thus the global condition permits learning rates up to η ≂ B β (1 − ρ). Combining this scale with the order-one local ceiling yields (48). 32

C

SGD scaling law

Theorem C.1 (SGD scaling law (Li et al., 2026a)). Under Assumptions 3.1 and 3.3, if the stability condition η ≲ 1 holds, then for k ≳ 1/η, the expected excess risk satisfies: E[E(θk )] ≂ (ηk)−s +

ησ 2 . B

(49)

Proof. Li et al. (2026a, Theorem 5.2) establishes that for a finite-width model of capacity M , the expected excess risk under constant learning rate η and batch size B scales as: E[E(θk )] ≂ M −sβ + (ηk)−s +

ησ 2 . B

Our formulation (Section 3.1) considers the infinite-dimensional Hilbert space H. Taking the limit M → ∞, the approximation error M −sβ → 0. Under the matching noise covariance structure (Proposition 3.4), the result immediately follows.

D

Polyak scaling law

In this section, we establish Theorem 5.1 using the Functional Scaling Law (FSL) framework of Li et al. (2026a). Our starting point is the following expected excess risk: k−1 X E[E(θk )] ≂ Spolyak (k) + Kpolyak (τ ) | {z } τ =0 | {z } signal learning

(50)

forgetting kernel

{z

| k−1

≂ Spolyak (k) +

η2 σ 2 + E[E(θk−1−τ )] B }

noise accumulation k−1 2 X

η2σ2 X η Kpolyak (τ ) + Kpolyak (τ )E[E(θk−1−τ )] . B B τ =0 τ =0 | {z } | {z } additive label noise

(51)

multiplicative noise

Our analysis proceeds in three steps: • First, Section D.1 derives the general discrete risk recursion (51) from an exact coordinate-wise impulse-response representation. Section D.2 then analyzes the characteristic roots of the second-order dynamics and separates the spectrum into overdamped and underdamped modes. • Second, Sections D.3 and D.4 deal with the stochastic terms in the recursion. Section D.3 computes the additive label noise, while Section D.4 shows that, under the risk stability condition, the multiplicative noise contribution can be absorbed into the label noise term up to constant factors. • Finally, Sections D.5 and D.6 characterize the deterministic signal learning term Spolyak (k) in the fully overdamped and mixed-damping regimes, and therefore derive the scaling law expression.

D.1

Signal–noise decomposition and risk recursion

Recall the parameter error ek := θk − θ ∗ . Substituting the gradient decomposition into the Polyak update (4), we obtain the error recursion:
ek+1 = (1 + ρ)I − ηH ek − ρek−1 − ηξk (θk ). 33

In the diagonal eigenbasis of the population Hessian H = diag(λ1 , λ2 , . . .), the j-th coordinate of the parameter error, denoted as ek (j), obeys the following second-order scalar difference equation: ek+1 (j) = (1 + ρ − ηλj )ek (j) − ρek−1 (j) − ηξk (j), where ξk (j) is the j-th coordinate of the stochastic gradient noise vector ξk (θk ). The homogeneous characteristic polynomial associated with the deterministic part of this coordinate-wise recursion is y 2 − (1 + ρ − ηλj )y + ρ = 0. For each coordinate j, let yj,+ and yj,− denote the two roots of this polynomial. We define the discrete Green function (impulse response) for the j-th coordinate as Gk (j) :=

k+1 k+1 yj,+ − yj,−

yj,+ − yj,−

,

k ⩾ 0,

with the natural continuous limiting definition when yj,+ = yj,− (i.e., critical damping). By viewing the stochastic noise sequence −ηξk (j) as an external forcing term, the stochastic parameter error trajectory can be decomposed exactly into a deterministic signal trajectory and a discrete convolution of the gradient noise. This exact impulse-response representation is formalized in the following proposition. Proposition D.1 (Impulse-response representation). For any step k ⩾ 0, the coordinate-wise error trajectory of Polyak satisfies: ek (j) = edet k (j) − η

k−1 X

Gk−1−m (j)ξm (j),

(52)

m=0

where edet k (j) is the unperturbed deterministic trajectory solving: det det edet k+1 (j) = (1 + ρ − ηλj )ek (j) − ρek−1 (j), ∗ det ∗ with initializations edet 0 (j) = e0 (j) = −θj and e−1 (j) = e−1 (j) = −θj (assuming θ0 = θ−1 = 0).

Proof. Define the stochastic residual error ẽk (j) := ek (j)−edet k (j). Subtracting the deterministic trajectory from the stochastic parameter error yields the inhomogeneous system: ẽk+1 (j) − (1 + ρ − ηλj )ẽk (j) + ρẽk−1 (j) = −ηξk (j),

(53)

subject to the initial conditions: ẽ0 (j) = ẽ−1 (j) = 0. We introduce the discrete linear difference operator Lj : (Lj u)k := uk+1 − (1 + ρ − ηλj )uk + ρuk−1 . By the definition of the discrete Green function Gk (j): (Lj G(j))k = 0

(∀k ⩾ 0),

with

G−1 (j) = 0,

G0 (j) = 1.

Let the impulse response to a point source at step m ⩾ 0 be defined as Hk (j; m) := Gk−1−m (j), which satisfies: ( 1, k = m, (Lj H(j; m))k = δk,m , where δk,m = 0, k ̸= m. 34

Applying the superposition principle to the linear operator Lj under the forcing term −ηξk (j) in (53) yields: ẽk (j) =

k−1 X

− ηξm (j) Hk (j; m)

m=0 k−1 X

= −η

Gk−1−m (j)ξm (j).

m=0

Substituting the definition of ẽk (j) back into the equation completes the proof: k−1 X

ek (j) = edet k (j) − η

Gk−1−m (j)ξm (j).

m=0

To quantify the excess risk, we define the deterministic bias (the signal component) and the convolution kernel as follows: X 1X 2 det Spolyak (k) := Kpolyak (k) := λ2j Gk (j)2 . λj edet k (j) = E(θk ), 2 j⩾1

j⩾1

The following proposition establishes that the expected excess risk decomposes, up to constant factors, into the deterministic signal bias, the accumulated additive label noise, and the accumulated multiplicative gradient noise. Proposition D.2 (Risk recursion for Polyak). Under the gradient-noise covariance assumption, the expected excess risk obeys the two-sided recursion: k−1

E[E(θk )] ≂ Spolyak (k) +

k−1

η2σ2 X η2 X Kpolyak (τ ) + Kpolyak (τ )E[E(θk−1−τ )]. B B τ =0

(54)

τ =0

Equivalently, the right-hand side provides both rigorous upper and lower bounds for the expected excess risk up to universal constant factors. Proof. By substituting the Green function representation (52) into the definition of the excess risk, we analyze the second moment for each error coordinate: 2

ek (j) =

k−1 X

edet k (j) − η

!2 Gk−1−m (j)ξm (j)

.

m=0

Because the stochastic gradient noise is conditionally mean-zero (E[ξm (j) | θm ] = 0) and independent across distinct time steps, all cross-terms between the deterministic trajectory and the noise, as well as cross-time noise products, vanish upon taking the expectation. Hence, 2 2 E[ek (j)2 ] = edet k (j) + η

k−1 X

Gk−1−m (j)2 E[ξm (j)2 ].

m=0

Multiplying by 12 λj and summing over all coordinates j ⩾ 1 gives the expected excess risk: k−1

E[E(θk )] ≂

1X η2 X X 2 λj edet (j) + λj Gk−1−m (j)2 E[ξm (j)2 ]. k 2 2 j⩾1

m=0 j⩾1

35

According to the property of curvature-aligned noise defined in Eq. (6), for the coordinate basis uj , the variance satisfies the following: E[ξm (j)2 ] ≂

1 λj E[E(θm )] + σ 2 , B

where the factor 1/B accounts for the mini-batching variance reduction. Substituting this into the previous equation yields: E[E(θk )] ≂

k−1
X 2 1X η2 X 2 2 σ + E[E(θm )] λj Gk−1−m (j)2 λj edet (j) + k 2 B j⩾1

= Spolyak (k) +

m=0

j⩾1

k−1

k−1

m=0

m=0

η2σ2 X η2 X Kpolyak (k − 1 − m) + Kpolyak (k − 1 − m)E[E(θm )]. B B

Finally, performing a change of variables by setting τ = k − 1 − m yields: k−1

k−1

τ =0

τ =0

η2σ2 X η2 X E[E(θk )] ≂ Spolyak (k) + Kpolyak (τ ) + Kpolyak (τ )E[E(θk−1−τ )]. B B The two-sided equivalence follows directly from the two-sided comparability in the noise covariance assumption and the non-negativity of Spolyak (k) and Kpolyak (k).

D.2

Overdamped and underdamped modes decomposition

Since both the deterministic error signal edet k (j) and the discrete Green function Gk (j) obey the same homogeneous difference equation, their dynamic behaviors are simultaneously dictated by the same characteristic roots. Depending on the sign of the corresponding discriminant, these roots transition from real to complex conjugate pairs, driving the system from monotonic decay to oscillatory transients. The following proposition rigorously establishes the damping threshold, defines the crossover index Jmax , and partitions the spectrum into two distinct global regimes. Proposition D.3 (Exact spectral decomposition and crossover boundary). For each coordinate j, the characteristic roots of the Polyak recursion are given by p (1 + ρ − ηλj ) ± ∆j yj,± = , 2 where the discriminant is ∆j := (1 + ρ − ηλj )2 − 4ρ. √ The critical transition ∆j = 0 corresponds exactly to the spectral threshold ηλj = (1 − ρ)2 , which is asymptotically equivalent to ηλj ≂ (1 − ρ)2 . This threshold classifies the individual eigen-directions into: √ 1. Overdamped Modes (ηλj < (1 − ρ)2 , i.e., ∆j > 0): The roots are strictly real. √ 2. Underdamped Modes (ηλj > (1− ρ)2 , i.e., ∆j < 0): The roots are complex conjugates. In the mixed-damping regime, the boundary index Jmax separating the underdamped modes from the overdamped modes is defined as the crossover index satisfying ηλJmax ≂ (1 − ρ)2 :  Jmax ≂

η (1 − ρ)2

1/β .

(55)

Consequently, based on the maximum eigenvalue λmax of the population Hessian H, the algorithm globally operates in one of two regimes: 36

• Fully Overdamped Regime (ηλmax ≲ (1 − ρ)2 ): All spectral modes are overdamped. • Mixed-Damping Regime (ηλmax ≳ (1 − ρ)2 ): The spectrum bifurcates into high-curvature underdamped modes (j ⩽ Jmax ) and flat overdamped modes (j > Jmax ). Proof. The phase transition between real and complex roots occurs exactly at the boundary where the discriminant vanishes (∆j = 0). Solving the quadratic equation (1 + ρ − ηλj )2 = 4ρ for the relevant stable branch yields the exact threshold: √ √ ηλj = 1 + ρ − 2 ρ = (1 − ρ)2 . √ By multiplying the numerator and denominator by (1+ ρ)2 , we can rewrite this exact threshold √ (1−ρ)2 √ 2 . Since 1 + ρ ≂ 2 for all ρ ∈ (0, 1), we obtain the fundamental scaling as: ηλj = (1+ ρ) equivalence for the damping boundary: ηλj ≂ (1 − ρ)2 . To locate the exact crossover boundary in the mixed-damping regime, we apply the power-law spectrum λj ≂ j −β . Setting the critical damping boundary condition ηλJmax ≂ (1 − ρ)2 , we 2 −β −β obtain: ηJmax ≂ (1 − ρ)2 ⇐⇒ Jmax ≂ (1−ρ) . Taking the −1/β-th power of both sides yields: η  Jmax ≂

η (1 − ρ)2

1/β ,

which establishes that the underdamped modes are exactly indexed by j ⩽ Jmax and the overdamped modes by j > Jmax . √ For an overdamped mode where ηλj < (1 − ρ)2 , the discriminant is strictly positive (∆j > 0), and the roots yj,+ and yj,− are real and distinct. √ For an underdamped mode where ηλj > (1 − ρ)2 , the discriminant is strictly negative √ (∆j < 0). The roots become a complex conjugate pair yj,± = ρe±iωj , where the exact phase 1+ρ−ηλ angle ωj ∈ (0, π) satisfies: cos ωj = 2√ρ j . Finally, to determine the global behavior of the system, we examine the maximum eigenvalue λmax . If ηλmax ≲ (1−ρ)2 , then every coordinate satisfies ∆j ⩾ 0 in the scaling sense, placing the system entirely in the Fully Overdamped Regime. If ηλmax ≳ (1 − ρ)2 , the spectral modes with λj ≳ (1−ρ)2 /η form a non-empty underdamped subset, while the remaining tail λj ≲ (1−ρ)2 /η remains overdamped, placing the system globally in the mixed-damping regime.

D.3

Label noise calculation

Before bounding the full multiplicative risk, we must isolate the purely additive label-noise contribution, which establishes the noise floor for both regimes. Lemma D.4 (Label-noise contribution). There exist constants C, c− , c+ > 0 such that if   1 1−ρ k ⩾ C max , , (56) 1−ρ η then the finite-time label-noise contribution satisfies k−1

c−

σ2η η2σ2 X σ2η ⩽ Kpolyak (τ ) ⩽ c+ . B(1 − ρ) B B(1 − ρ) τ =0

In particular, k−1

η2σ2 X σ2 η Kpolyak (τ ) ≂ . B B 1−ρ τ =0

37

(57)

Proof. We first compute the saturated coordinate Green energy exactly. Fix a coordinate j and write xj := ηλj , Aj := 1 + ρ − xj . The Green function satisfies Gk+1 (j) = Aj Gk (j) − ρGk−1 (j),

G0 (j) = 1,

G−1 (j) = 0.

Assume the coordinate is stable. Define ∞ X

Sj :=

2

Gτ (j) ,

Cj :=

τ =0

∞ X

Gτ (j)Gτ −1 (j),

τ =0

where G−1 (j) = 0. Multiplying the recursion by Gk (j) and summing over k ⩾ 0 gives Cj = Aj Sj . Aj Sj − ρCj , hence Cj = 1+ρ Next square the recursion and sum over k ⩾ 0:
Sj − 1 = A2j + ρ2 Sj − 2Aj ρCj . Substituting Cj = Aj Sj /(1 + ρ) gives "

# 2ρ 2A j 1 = Sj 1 − A2j − ρ2 + 1+ρ  1−ρ = Sj (1 + ρ)2 − (1 + ρ − xj )2 1+ρ 1−ρ = Sj xj (2(1 + ρ) − xj ) . 1+ρ

Therefore

∞ X

Gτ (j)2 =

τ =0

1+ρ . (1 − ρ)ηλj (2(1 + ρ) − ηλj )

(58)

We now pass from the saturated sum to the finite-time sum. The root formula for Gk (j) implies the following tail bound: there exist universal constants CG , cG > 0 such that    X ∞ ∞ X ηλj 2 Gτ (j)2 . (59) Gτ (j) ⩽ CG exp −cG k min 1 − ρ, 1−ρ τ =0

τ =k

This is the standard Polyak relaxation estimate. P Because λj ≂ j −β and β > 1, the spectral mass is finite: j⩾1 λj < ∞. P ∗ Choose a fixed integer J∗ , depending only on the spectrum constants and β, such that Jj=1 λj ⩾ P 1 λ . j⩾1 j 2 Since J∗ is fixed, λJ∗ is a positive constant independent of η, ρ, k, B. Thus, after increasing the universal constant C in (56) if necessary, condition (56) implies that for every 1 ⩽ j ⩽ J∗ ,   ηλj log(2CG ) ⩾ . k min 1 − ρ, 1−ρ cG By (59), for all j ⩽ J∗ ,

∞ X

∞

Gτ (j)2 ⩽

τ =0

τ =k

Therefore

k−1 X τ =0

1X Gτ (j)2 . 2

∞

1X Gτ (j)2 , Gτ (j) ⩾ 2 2

τ =0

38

1 ⩽ j ⩽ J∗ .

(60)

For all j, the trivial upper bound k−1 X

Gτ (j)2 ⩽

τ =0

∞ X

Gτ (j)2

(61)

τ =0

holds. We now insert these bounds into the finite-time label-noise term. Since all terms are nonnegative, we can exchange the τ - and j-sums: k−1

k−1

η2σ2 X 2 X η2σ2 X Kpolyak (τ ) = λj Gτ (j)2 . B B τ =0

j⩾1

τ =0

For the lower bound, use (60) and (58) on the first J∗ modes: k−1

J

τ =0

j=1

∗ 1 η2σ2 X 1+ρ η2σ2 X Kpolyak (τ ) ⩾ λ2j B 2 B (1 − ρ)ηλj (2(1 + ρ) − ηλj )

J

∗ λj 1 ησ 2 (1 + ρ) X . = 2 B(1 − ρ) 2(1 + ρ) − ηλj

j=1

Since 2(1 + ρ) − ηλj ⩽ 2(1 + ρ),

1+ρ 1 ⩾ . 2(1 + ρ) − ηλj 2

Thus k−1

J

τ =0

j=1

∗ η2σ2 X 1 ησ 2 X Kpolyak (τ ) ⩾ λj . B 4 B(1 − ρ)

P ∗ λj is a positive constant, so this gives the lower bound in (57) for a suitable The finite sum Jj=1 constant c− > 0. For the upper bound, use (61) and (58) over all modes: k−1

λj η2σ2 X ησ 2 (1 + ρ) X Kpolyak (τ ) ⩽ B B(1 − ρ) 2(1 + ρ) − ηλj τ =0

j⩾1

⩽

ησ 2 (1 + ρ)

X 1 λj . B(1 − ρ) 2(1 + ρ) − ηλ1 j⩾1

The last factor is finite, and the denominator is bounded below by the stability margin. Hence the upper bound in (57) holds for a suitable constant c+ > 0.

D.4

Uniform boundedness of expected risk

As long as the additive label noise is strictly positive (σ > 0), the scaling law derivation can be simplified. We can prove that the expected risk is uniformly bounded, which allows the multiplicative gradient noise to be completely absorbed by the additive label noise floor. Lemma D.5 (Uniform boundedness of expected risk). Suppose σ > 0. Under the stability condition η ≲ 1 and η ⩽ cB(1 − ρ) for a sufficiently small universal constant c > 0, the expected excess risk of Polyak is uniformly bounded for all k ⩾ 0: E[E(θk )] ⩽ M, where M > 0 is a constant independent of k, η, B, and ρ. 39

Proof. We proceed by induction on k. By Proposition D.2, there exist universal constants C1 , C2 , C3 > 0 such that for any k ⩾ 1: E[E(θk )] ⩽ C1 Spolyak (k) + C2

k−1

k−1

τ =0

τ =0

η2σ2 X η2 X Kpolyak (τ ) + C3 Kpolyak (τ )E[E(θk−1−τ )]. B B

(62)

First, Propositions D.7 and D.9 in the following subsections establish that the deterministic bias monotonically decays or is uniformly bounded by its initialization across all regimes. Thus: Spolyak (k) ⩽ Ch2 E(θ0 ) =: M1 . Second, applying

P∞

CK τ =0 Kpolyak (τ ) ⩽ η(1−ρ) from Lemma D.4 bounds the label noise: k−1

C2

CK η η2σ2 X Kpolyak (τ ) ⩽ C2 σ 2 . B B(1 − ρ) τ =0

Under η ⩽ cB(1 − ρ), this is uniformly bounded by C2 σ 2 CK c =: M2 . Assuming inductively E[E(θm )] ⩽ M for m < k, the multiplicative noise term yields: k−1

η2 X η2 C3 Kpolyak (τ )E[E(θk−1−τ )] ⩽ C3 B B τ =0



CK η(1 − ρ)

 M = C3 CK

η M ⩽ C3 CK cM. B(1 − ρ)

Choosing c ⩽ 2C31CK limits this term to M/2. Substituting into (62) yields: E[E(θk )] ⩽ M1 + M2 +

M . 2

Setting M = 2(M1 + M2 ) guarantees E[E(θk )] ⩽ M . The base case k = 0 holds trivially since E[E(θ0 )] = Spolyak (0) ⩽ M1 ⩽ M . This completes the induction. Corollary D.6 (Absorption of multiplicative noise). For σ > 0, under the stability condition η ≲ min{1, B(1 − ρ)}, the multiplicative gradient noise is strictly bounded by the additive label noise up to a constant factor: ! k−1 k−1 η2 X M η2σ2 X Kpolyak (τ )E[E(θk−1−τ )] ⩽ 2 Kpolyak (τ ) . B σ B τ =0

τ =0

Consequently, the multiplicative noise can be absorbed, and the expected excess risk is asymptotically equivalent to the sum of the deterministic bias and the label noise: k−1

E[E(θk )] ≂ Spolyak (k) +

η2σ2 X Kpolyak (τ ). B τ =0

D.5

Risk dynamics in the fully overdamped regime

In the fully overdamped regime, the deterministic error trajectory decays monotonically, yielding a polynomial bias. Proposition D.7 (Deterministic risk, fully overdamped regime). Under Assumptions 3.1 and 3.3, in the fully overdamped regime η ≲ (1 − ρ)2 , we have:  1, k ≲ 1−ρ η ,   s Spolyak (k) ≂ 1−ρ 1−ρ  ηk , k≳ η . 40

Proof. In the fully overdamped regime, we can assume: ηλ1 ⩽ c0 (1 − ρ)2 , with a sufficiently small universal constant c0 > 0. det det ∗ Let the deterministic coordinate multiplier hk (j) := edet k (j)/e0 (j) = −ek (j)/θj satisfy hk+1 (j) = (1 + ρ − ηλj )hk (j) − ρhk−1 (j),

h1 (j) = 1 − ηλj .

h0 (j) = 1,

The characteristic roots are rj,+ , rj,− =

1 + ρ − ηλj ±

p (1 + ρ − ηλj )2 − 4ρ . 2

We first give explicit root bounds. Since η ⩽ λc01 (1−ρ)2 and λj ⩽ λ1 , we have 0 ⩽ ηλj ⩽ c0 (1−ρ)2 for all j. Put r = 1 − u. The characteristic polynomial evaluated at 1 − u is (1 − u)2 − (1 + ρ − ηλj )(1 − u) + ρ = ηλj − (1 − ρ + ηλj )u + u2 . Choosing c0 > 0 sufficiently small, the two evaluations 2  ηλj ηλj + > 0, ηλj − (1 − ρ + ηλj ) 2(1 − ρ) 2(1 − ρ) and 2ηλj ηλj − (1 − ρ + ηλj ) + 1−ρ



2ηλj 1−ρ

2 <0

hold uniformly for 0 < ηλj ⩽ c0 (1 − ρ)2 . Therefore the slow root satisfies 1−

2ηλj ηλj ⩽ rj,+ ⩽ 1 − . 1−ρ 2(1 − ρ)

(63)

The fast root satisfies rj,+ rj,− = ρ = 1 − (1 − ρ). Since rj,+ ⩾ 1 − 2c0 (1 − ρ), after decreasing c0 if needed, 1−ρ 0 ⩽ rj,− ⩽ 1 − . (64) 2 Now write the solution as k k hk (j) = Aj rj,+ + Bj rj,− . From h0 = 1 and h1 = 1 − ηλj , we obtain the coefficients Aj =

1 − ηλj − rj,− , rj,+ − rj,−

Bj = 1 − Aj =

rj,+ − (1 − ηλj ) . rj,+ − rj,−

Using (63)–(64), and taking c0 sufficiently small, these coefficients obey uniform bounds: 1 ⩽ Aj ⩽ 2, 2

|Bj | ⩽ 2.

(65)

Equivalently, applying the initial conditions, the exact solution can be expressed in the fractional form:     k−1 k−1 k − rk rj,+ rj,− −rj,+ + rj,− + (1 − ηλj ) rj,+ j,− hk (j) = . (66) rj,+ − rj,− P To evaluate the deterministic bias Spolyak (k) = j⩾1 λj (θj∗ )2 hk (j)2 , we divide the time k into three successive stages. 1 . Stage 1: k ≲ 1−ρ In this extremely short-time regime, taking the limits rj,+ ≈ rj,− ≈ r ≈ 1 and using r2 ≈ ρ in (66) yields the heuristic trajectory: hk (j) ≃ −ρ(k − 1)rk−2 + (1 − ηλj )krk−1 = k ((1 − ηλj )r − ρ) rk−2 + ρrk−2 . 41

Since (1−ηλj )r −ρ = O(1−ρ), the first term is bounded by O(k(1−ρ)) ≲ 1, yielding hk (j) ≃ 1. Rigorously, for the leading modes (j = O(1)), since ηλj ⩽ c0 (1 − ρ)2 , the slow root satisfies 1 rj,+ ⩾ 1 − 2c0 (1 − ρ). For k ≲ 1−ρ , we have C

k rj,+ ⩾ (1 − 2c0 (1 − ρ)) 1−ρ ≂ 1.

Since the dynamics are unconditionally stable (hk (j) ⩽ 1), it follows that hk (j) ≂ 1 for the leading modes. Consequently, the bias is dominated by these non-decayed modes: X Spolyak (k) ≂ λj (θj∗ )2 ≂ 1. j⩾1 1 Stage 2: 1−ρ ≲ k ≲ 1−ρ η .

k  1 k ⩽ 1 − 1−ρ For k ≳ 1−ρ , the fast root decays exponentially since rj,− ≪ 1. The multiplier is 2 dominated by the slow root: hk (j) ≃

1 − ηλj − rj,− k k r ≃ rj,+ . rj,+ − rj,− j,+

For the leading modes with λj = Θ(1), within the time window k ≲ 1−ρ η , the exponent is bounded by:     ηλj k ηλj k rj,+ ≂ 1 − c ≂ exp −c k ≂ 1, 1−ρ 1−ρ ηλ

where we used the condition 1−ρj k ≲ λj ⩽ λ1 = O(1). Thus, the leading modes have not yet experienced significant exponential decay, preserving hk (j) ≂ 1. As a result, the deterministic bias remains unchanged: Spolyak (k) ≂ 1. Stage 3: k ≳ 1−ρ η . In this long-time regime, the exponential decay of the slow root becomes significant across the spectrum. Combining (63), (64), and (65), there are universal constants c1 , c2 , c3 , C1 > 0 such that, for all k ⩾ C1 (1 − ρ)−1 ,     ηλj ηλj c1 exp −c2 (67) k ⩽ |hk (j)| ⩽ C1 exp −c3 k . 1−ρ 1−ρ Substituting (67) into the deterministic bias gives     X X ηλj ηλj ∗ 2 ∗ 2 c k ⩽ Spolyak (k) ⩽ C k . λj (θj ) exp −C λj (θj ) exp −c 1−ρ 1−ρ j⩾1

(68)

j⩾1

By Assumption 3.3, λj (θj∗ )2 ≂ j −1−βs ,

λj ≂ j −β .

Thus it remains to bound sums of the form   X ηk −β −1−βs j exp −c j . 1−ρ j⩾1

ηk The condition on k implies 1−ρ ⩾ C.

For the lower bound, set J :=



ηk 1−ρ

1/β  .

ηk −β For j ⩾ J, we have 1−ρ j ⩽ 1, so

ηk −β exp(−C 1−ρ j ) ⩾ e−C . Therefore,     X X ηk −s ηk −β j −1−βs exp −C j ⩾ e−C j −1−βs ⩾ c . 1−ρ 1−ρ j⩾1

j⩾J

42

ηk −β For the upper bound, the integral comparison with u = c 1−ρ x gives     Z ∞ X ηk −β ηk −β −1−βs −1−βs x exp −c ⩽C dx j exp −c j x 1−ρ 1−ρ 1 j⩾1

 =C

ηk 1−ρ

−s Z c ηk

1−ρ

s−1 −u

u

e

 du ⩽ C

0

ηk 1−ρ

−s .

Combining the upper and lower bounds in (68) yields   1−ρ s Spolyak (k) ≂ . ηk Synthesizing the three stages completes the proof. Theorem D.8 (Scaling law, fully overdamped regime). Under Assumptions 3.1 and 3.3, and the stability condition η ≲ min{1, B(1 − ρ)}, for the fully overdamped regime η ≲ (1 − ρ)2 , when k ≳ 1−ρ η , the expected excess risk satisfies:   1−ρ s σ2η E[E(θk )] ≂ + . ηk B(1 − ρ) Proof. By Corollary D.6, the expected risk is completely determined by the deterministic bias and the label noise: k−1 η2σ2 X E[E(θk )] ≂ Spolyak (k) + Kpolyak (τ ). B τ =0

For k ≳ 1−ρ η , we substitute the deterministic bias from Proposition D.7:   1−ρ s Spolyak (k) ≂ . ηk We substitute the label noise floor from Lemma D.4: k−1

η2σ2 X σ2η Kpolyak (τ ) ≂ . B B(1 − ρ) τ =0

Summing these two terms directly yields the scaling law.

D.6

Risk dynamics in the mixed-damping regime

In the mixed-damping regime, the spectrum is partitioned into underdamped and overdamped  1/β η modes by the crossover index Jmax ≂ (1−ρ) . Consequently, the deterministic bias 2 decomposes into an overdamped polynomial tail and an underdamped oscillatory transient, denoted by Sover (k) and Sunder (k). Due to destructive interference among the underdamped modes, Sunder (k) lacks a strictly positive pointwise lower bound. Proposition D.9 (Deterministic bias, mixed-damping regime). Under Assumptions 3.1 and 3.3, for the mixed-damping regime η ≳ (1 − ρ)2 , the deterministic bias decomposes into Spolyak (k) = Sover (k) + Sunder (k), where the overdamped tail satisfies the two-sided bounds:   2 s 1  (1−ρ) , k ≲ 1−ρ , η   Sover (k) ≂ s 1  1−ρ , k ≳ 1−ρ , ηk and Sunder (k) represents the underdamped transient whose properties are formally established in Lemma D.10. 43

∗ Proof. Let the coordinate multiplier hk (j) := −edet k (j)/θj . Then it satisfies

hk+1 (j) = (1 + ρ − ηλj )hk (j) − ρhk−1 (j),

h0 (j) = 1,

h1 (j) = 1 − ηλj .

The characteristic discriminant is negative precisely when √ √ (1 − ρ)2 < ηλj < (1 + ρ)2 . Thus the comparison between η and (1 − ρ)2 determines whether a coordinate is overdamped or underdamped. In the mixed-damping regime η ≳ (1 − ρ)2 , the crossover index is defined by ηλJmax ≂ (1 − ρ)2 . We decompose the deterministic bias into the overdamped tail (j > Jmax ) and the underdamped leading block (j ⩽ Jmax ). For overdamped coordinates (j > Jmax ) in the stable small-step range, the root comparison used in Proposition D.7 gives constants c1 , C1 > 0 such that   ηλj |hk (j)| ⩽ C1 exp −c1 k . (69) 1−ρ Moreover, if ηλj ⩽ c1 (1 − ρ)2 and k ⩾ C1 (1 − ρ)−1 , the slow root dominates the fast root and   ηλj |hk (j)| ⩾ c1 exp −C1 k . 1−ρ P To evaluate the overdamped contribution Sover (k) = 12 j>Jmax λj (θj∗ )2 hk (j)2 , we discuss two cases for k: 1 Case 1: k ≲ 1−ρ . In this regime, for all overdamped modes j > Jmax , we have ηλj ≲ (1 − ρ)2 . Consequently, the 2 ηλ 1 exponent satisfies 1−ρj k ≲ (1−ρ) 1−ρ 1−ρ = 1. The exponential decay is bounded away from zero, so |hk (j)| ≂ 1. Therefore, the sum is dominated by its lower limit:  s X (1 − ρ)2 −1−βs −βs Sover (k) ≂ j ≂ (Jmax ) ≂ . η j>Jmax

1 . Case 2: k ≳ 1−ρ

Since λj ≂ j −β , the time condition implies that every j ⩾



Cηk 1−ρ

1/β

satisfies ηλj ⩽ c1 (1 − ρ)2 ,

after increasing C if necessary. Using λj (θj∗ )2 ≂ j −1−βs , we get   X 1−ρ s −1−βs . Sover (k) ⩾ c j ⩾c ηk 1/β

(70)

j⩾(Cηk/(1−ρ))

Conversely, by (69) and the integral comparison,     X ηk −β 1−ρ s −1−βs . Sover (k) ⩽ C j exp −c1 j ⩽C 1−ρ ηk j⩾1

Therefore the overdamped tail contributes exactly the polynomial term, up to constants. The remaining deterministic bias comes from the underdamped modes, defined as 1P Sunder (k) := 2 j⩽Jmax λj (θj∗ )2 hk (j)2 , which will be analyzed in Lemma D.10. Lemma D.10 (Properties of the underdamped transient). Let Sunder (k) denote the deterministic bias contributed by the underdamped modes (j ⩽ Jmax ). Under the conditions of Proposition D.9, for every step k ⩾ 0, we have: 0 ⩽ Sunder (k) ⩽ Cρk . 44

Proof. For j ⩽ Jmax , write the two roots as √

ρeiωj ,

√

ρe−iωj ,

where

cos ωj =

1 + ρ − ηλj . √ 2 ρ

Solving the two initial conditions gives the exact identity   1 − ρ − ηλj k/2 cos(kωj ) + √ hk (j) = ρ sin(kωj ) . 2 ρ sin ωj Set Aj :=

(71)

1 − ρ − ηλj . √ 2 ρ sin ωj

On the separated underdamped band Cmom (1 − ρ)2 ⩽ ηλj ⩽ cst /2, one has 4ρ sin2 ωj = ηλj − (1 −

√

√ ρ)2 (1 + ρ)2 − ηλj ⩾ cηλj ,

and (1 − ρ − ηλj )2 ⩽ (1 − ρ)2 , hence |Aj | ⩽ C. Thus (71) immediately gives the upper bound Sunder (k) =

1 X λj (θj∗ )2 hk (j)2 ⩽ Cρk , 2 j⩽Jmax

P because j⩾1 λj (θj∗ )2 < ∞. The lower bound is trivially zero due to the potential for destructive interference of the oscillatory terms at specific iterations. Theorem D.11 (Scaling law, mixed-damping regime). Under Assumptions 3.1 and 3.3, and 2 the stability n condition o η ≲ min{1, B(1 − ρ)}, for the mixed-damping regime η ≳ (1 − ρ) with 1 k ≳ max 1−ρ , 1−ρ , the expected excess risk satisfies: η  E[E(θk )] ≂ Sunder (k) +

1−ρ ηk

s +

σ2η , B(1 − ρ)

(72)

where the properties of the underdamped transient Sunder (k) are established in Lemma D.10. Proof. By Corollary D.6, the multiplicative noise is absorbed by the additive label noise floor, yielding: k−1 η2σ2 X E[E(θk )] ≂ Sunder (k) + Sover (k) + Kpolyak (τ ). (73) B τ =0 n o 1 By Lemma D.4, for k ≳ max 1−ρ , 1−ρ , the label noise summation captures the saturated η floor: k−1 η2σ2 X σ2η Kpolyak (τ ) ≂ . B B(1 − ρ) τ =0

1 By Proposition D.9, since k ≳ 1−ρ , the overdamped tail reaches its terminal polynomial decay phase:   1−ρ s Sover (k) ≂ . ηk

Substituting these components directly into (73) establishes the equivalence (72). 45

E

Nesterov scaling law

As for Polyak momentum, we establish Theorem 5.2 using the Functional Scaling Law (FSL) framework of Li et al. (2026a). Our starting point is the following expected excess risk: E[E(θk )] ≂ Snesterov (k) + | {z } signal learning

k−1 X τ =0

η2 Knesterov (τ ) σ 2 + E[E(θk−τ )] | {z } B

forgetting kernel

{z

| k−1

≂ Snesterov (k) +

(74)

}

noise accumulation k−1 2 X

η2σ2 X η Knesterov (τ ) + Knesterov (τ )E[E(θk−τ )] . B B τ =0 τ =0 {z } | {z } | additive label noise

(75)

multiplicative noise

Our analysis proceeds in three steps. First, Section E.1 derives the general discrete risk recursion (75) from an exact coordinate-wise impulse-response representation. Section E.2 then analyzes the characteristic roots of the second-order dynamics and separates the spectrum into overdamped and underdamped modes. Second, Sections E.3 and E.4 deal with the stochastic terms in the recursion. Section E.3 computes the additive label-noise contribution, which yields a lower noise floor than Polyak; Section E.4 shows that, under the risk stability condition, the multiplicative noise contribution can be absorbed into the label noise term up to constant factors. Finally, Sections E.5 and E.6 characterize the deterministic signal learning term Snesterov (k) in the fully overdamped and mixed-damping regimes, and therefore derive the scaling law expression.

E.1

Signal–noise decomposition and risk recursion

Recall the parameter error ek := θk − θ ∗ . For Nesterov, the gradient is evaluated at the look-ahead point θkla = θk + ρ(θk − θk−1 ). Let la ∗ ela k := θk − θ

denote the look-ahead error. Substituting the gradient decomposition into the Nesterov update (5), the error recursion in vector form is: la ek+1 = (I − ηH)ela k − ηξk (θk ).

In the diagonal eigenbasis of H, since ela k = (1+ρ)ek −ρek−1 , the j-th coordinate of the Nesterov error recursion satisfies: ek+1 (j) = (1 + ρ)(1 − ηλj )ek (j) − ρ(1 − ηλj )ek−1 (j) − ηξkla (j), where ξkla (j) is the j-th coordinate of the stochastic gradient noise ξk (θkla ) evaluated at the look-ahead point. Its homogeneous characteristic polynomial is y 2 − (1 + ρ)(1 − ηλj )y + ρ(1 − ηλj ) = 0. Let yj,+ and yj,− be the two roots of y 2 − (1 + ρ)(1 − ηλj )y + ρ(1 − ηλj ) = 0. 46

Define Gnesterov (j) := k

k+1 k+1 yj,+ − yj,−

yj,+ − yj,−

,

k ⩾ 0,

with the usual limiting definition at a repeated root. Then ek (j) = edet k (j) − η

k−1 X

la Gnesterov k−1−m (j)ξm (j),

(76)

m=0

where edet k (j) is the deterministic Nesterov error trajectory. Proposition E.1 (Impulse-response representation). The Green function is the unique solution of Gnesterov (j) = (1+ρ)(1−ηλj )Gnesterov (j)−ρ(1−ηλj )Gnesterov (j), k+1 k k−1

Gnesterov (j) = 1, 0

Gnesterov (j) = 0. −1

∗ det ∗ Moreover, if edet 0 (j) = −θj and e1 (j) = −(1 − ηλj )θj , then

edet k (j) = −

k − (1 − ηλ − y k (1 − ηλj − yj,− )yj,+ j j,+ )yj,− ∗ θj . yj,+ − yj,−

Proof. Fix a coordinate j. The inhomogeneous scalar recursion is ek+1 (j) = (1 + ρ)(1 − ηλj )ek (j) − ρ(1 − ηλj )ek−1 (j) − ηξkla (j). (j) satisfying The impulse response for this second-order recursion is the sequence Gnesterov k (j), (j)−ρ(1−ηλj )Gnesterov (j) = (1+ρ)(1−ηλj )Gnesterov Gnesterov k−1 k k+1

Gnesterov (j) = 1, 0

Gnesterov (j) = 0. −1

la (j) is inserted at time m, then its contribution to e la Indeed, if an impulse −ηξm m+1 (j) is −ηξm (j), la its contribution to em+2 (j) is −η(1 + ρ)(1 − ηλj )ξm (j), and subsequent contributions obey the la same homogeneous recursion. Hence the contribution at time k is −ηGnesterov k−1−m (j)ξm (j) for 0 ⩽ m ⩽ k − 1. Summing over all previous impulses gives

ek (j) = edet k (j) − η

k−1 X

la Gnesterov k−1−m (j)ξm (j),

m=0

where edet k (j) solves the homogeneous recursion with the same initialization. It remains to compute the closed form of Gnesterov (j). If yj,+ ̸= yj,− , then the general k k k . The initial solution of the homogeneous recurrence is a linear combination of yj,+ and yj,− conditions Gnesterov (j) = 0 and Gnesterov (j) = 1 give −1 0 Gnesterov (j) = k

k+1 k+1 − yj,− yj,+

yj,+ − yj,−

,

k ⩾ 0.

k . Thus If yj,+ = yj,− , this formula is interpreted by continuity, giving Gnesterov (j) = (k + 1)yj,+ k the displayed formula is exactly the unique Green function of the Nesterov coordinate recursion. Now consider the deterministic trajectory. It satisfies det det edet k+1 (j) = (1+ρ)(1−ηλj )ek (j)−ρ(1−ηλj )ek−1 (j),

∗ edet 0 (j) = −θj ,

∗ edet 1 (j) = −(1−ηλj )θj .

k and y k . Solving For yj,+ ̸= yj,− , write the deterministic solution as a linear combination of yj,+ j,− the two-by-two linear system determined by the two displayed initial conditions gives

edet k (j) = −

k − (1 − ηλ − y k (1 − ηλj − yj,− )yj,+ j j,+ )yj,− ∗ θj . yj,+ − yj,−

The repeated-root case follows by taking the same continuous limit. This proves the claim. 47

To quantify the excess risk, we define the deterministic bias (the signal component) and the convolution kernel as follows: X
2
2 1X Snesterov (k) := λj edet (j) , K (k) := λ2j Gnesterov (j) . nesterov k k 2 j⩾1

j⩾1

Because the stochastic gradient is evaluated at the look-ahead point, define 1 ∥2 . Lk := ∥ela 2 k H

ela k := ek + ρ(ek − ek−1 ),

Lemma E.2 (Nesterov look-ahead risk control). Assume the stable small-step condition ηλ1 ⩽ cst , where cst > 0 is sufficiently small. Then the look-ahead risk admits the coordinate expansion Lm =

1X λj ((1 + ρ)em (j) − ρem−1 (j))2 . 2 j⩾1

Moreover, the deterministic look-ahead risk is comparable to the next deterministic Nesterov risk: 2 1X  det λj (1 + ρ)edet (j) − ρe (j) ≂ Snesterov (m + 1). m m−1 2 j⩾1

For the stochastic iterates, the corresponding expectation-level control is c EEm+1 − C

η2σ2 ⩽ ELm ⩽ C EEm+1 , B

where the constants are independent of m, η, ρ, B. Proof. The coordinate expansion follows directly from ela m = (1+ρ)em −ρem−1 and the definition of the H-norm. For the deterministic comparison, the deterministic Nesterov update gives, for each coordinate j,   det det edet (j) = (1 − ηλ ) (1 + ρ)e (j) − ρe (j) . j m+1 m m−1 Since 0 < 1 − ηλj ⩽ 1 and 1 − ηλj ⩾ 1 − cst , we have 

det (1 + ρ)edet m (j) − ρem−1 (j)

2

2 ≂ edet m+1 (j) .

Multiplying by λj and summing over j ⩾ 1 yields the deterministic claim. For the stochastic comparison, the actual Nesterov update in coordinate j is la em+1 (j) = (1 − ηλj ) ((1 + ρ)em (j) − ρem−1 (j)) − ηξm (j).

Conditioning on the past, the noise is mean zero at the look-ahead point. Therefore the mixed term vanishes and X EEm+1 ≂ λj (1 − ηλj )2 E ((1 + ρ)em (j) − ρem−1 (j))2 j⩾1

+ η2

X

2 la λj E ξm (j) .

j⩾1

The first term is comparable to ELm , because 1 − ηλj is bounded above and below by positive constants. For the second term, Proposition 3.4 gives X j⩾1

la λ j E ξm (j)

2

≲

X 2 1 2 σ + ELm λj . B j⩾1

48

P Since j⩾1 λ2j < ∞ and cst is chosen small, the term proportional to η 2 ELm /B is absorbed into the comparable first term. Thus  2 2 η σ EEm+1 ≂ ELm + O , B which gives c EEm+1 − C

η2σ2 ⩽ ELm ⩽ C EEm+1 . B

This proves the lemma. Proposition E.3 (Nesterov risk recursion). Under the covariance relation in Proposition 3.4, the expected excess risk satisfies the two-sided recursion k−1

k−1

τ =0

τ =0

η2σ2 X η2 X EEk ≂ Snesterov (k) + Knesterov (τ ) + Knesterov (τ )ELk−1−τ . B B

(77)

Moreover, Lemma E.2 gives the required control of the look-ahead risk by the next-step risk. In particular, under the stable small-step condition, c EEm+1 − C

η2σ2 ⩽ ELm ⩽ C EEm+1 . B

Proof. Substitute the Nesterov impulse-response representation (76) into the coordinate-wise excess risk. For every coordinate j, ek (j) = edet k (j) − η

k−1 X

la Gnesterov k−1−m (j)ξm (j).

m=0 la (j) is conditionally mean zero given the past and the look-ahead point ela , and the The noise ξm m fresh mini-batches are independent across time. Therefore all mixed deterministic–stochastic and cross-time stochastic terms vanish after taking expectation. Hence

E[ek (j)2 ] = edet k (j)

2

+ η2

k−1 X

2
2 la Gnesterov k−1−m (j) E ξm (j) ,

m=0

up to the same universal constants in the covariance comparison of Proposition 3.4. Multiplying by λj and summing over j ⩾ 1 gives EEk ≂

X

λj edet k (j)

2

+ η2

k−1 X X

2
2 la λj Gnesterov k−1−m (j) E ξm (j) .

m=0 j⩾1

j⩾1

Since Nesterov evaluates the stochastic gradient at the look-ahead point ela m , Proposition 3.4 gives the coordinate-wise variance relation la E ξm (j)

2

≂

1 λj σ 2 + ELm . B

Substituting this into the previous display yields k−1

EEk ≂

X j⩾1

λj edet k (j)

2

+

X 2 nesterov
2 η2 X 2 σ + ELm λj Gk−1−m (j) . B m=0

j⩾1

49

Using the definitions of Snesterov and Knesterov , and then changing variables τ = k − 1 − m, we obtain k−1 k−1 η2σ2 X η2 X EEk ≂ Snesterov (k) + Knesterov (τ ) + Knesterov (τ )ELk−1−τ . B B τ =0

τ =0

This proves the claimed two-sided recursion. It remains to control the look-ahead risk. This is exactly Lemma E.2, which expands Lm in coordinates and then uses the Nesterov update from step m to step m + 1 to compare the look-ahead risk with the next-step risk.

E.2

Overdamped and underdamped modes decomposition

nesterov (j) obey Since both the deterministic error signal edet k (j) and the discrete Green function Gk the same homogeneous difference equation, their behavior is dictated by the same characteristic roots. Depending on the sign of the discriminant, these roots transition from real roots to a complex-conjugate pair. The following proposition establishes the exact Nesterov damping threshold, defines the crossover index Jmax , and partitions the spectrum into the two global regimes.

Proposition E.4 (Exact spectral decomposition and crossover boundary). For each coordinate j, the characteristic roots of the Nesterov recursion are q (1 + ρ)(1 − ηλj ) ± ∆nesterov j yj,± = , 2 where ∆nesterov := (1 + ρ)2 (1 − ηλj )2 − 4ρ(1 − ηλj ). j In the stable small-step range 0 < ηλj < 1, the critical transition ∆nesterov = 0 corresponds j exactly to (1 − ρ)2 ηλj = , (1 + ρ)2 which is asymptotically equivalent to ηλj ≂ (1 − ρ)2 . Thus the individual eigen-directions are classified as follows:
1. Overdamped Modes ηλj < (1 − ρ)2 /(1 + ρ)2 : the roots are real.
2. Underdamped Modes ηλj > (1 − ρ)2 /(1 + ρ)2 : the roots form a complex-conjugate pair. In the mixed-damping regime, the boundary index Jmax separating the underdamped modes from the overdamped modes is defined by the crossover relation ηλJmax ≂ (1 − ρ)2 :  Jmax ≂

η (1 − ρ)2

1/β .

(78)

Consequently, based on the maximum eigenvalue λmax , Nesterov operates in one of two global regimes:
• Fully Overdamped Regime ηλmax ≲ (1 − ρ)2 : all spectral modes are overdamped.
• Mixed-Damping Regime ηλmax ≳ (1 − ρ)2 : the high-curvature modes j ⩽ Jmax are underdamped, whereas the flat modes j > Jmax are overdamped.

50

Proof. In the stable small-step range, 1 − ηλj ∈ (0, 1). The condition ∆nesterov < 0 is j (1 + ρ)2 (1 − ηλj )2 − 4ρ(1 − ηλj ) < 0. 2

(1−ρ) 4ρ 4ρ Since 1 − ηλj > 0, this is equivalent to 1 − ηλj < (1+ρ) 2 , and hence ηλj > 1 − (1+ρ)2 = (1+ρ)2 . Therefore the equality case gives the exact damping boundary. Since (1 + ρ)2 ≂ 1 for 0 ⩽ ρ < 1, this boundary has the scaling form ηλj ≂ (1 − ρ)2 . To locate the crossover index in the mixed-damping regime, use λj ≂ j −β and impose −β ηλJmax ≂ (1 − ρ)2 . Then ηJmax ≂ (1 − ρ)2 , which gives

 Jmax ≂

η (1 − ρ)2

1/β .

Because the eigenvalues are non-increasing, the modes j ⩽ Jmax lie above the damping threshold and are underdamped, while the modes j > Jmax lie below the threshold and are overdamped. Examining whether this crossover occurs below the leading eigenvalue yields the two global regimes stated above.

E.3

Label noise calculation

The convolution estimates in Appendix D.3 will be used throughout this section. It remains to identify the Nesterov label-noise floor and to replace the look-ahead risk in the multiplicative-noise recursion. Lemma E.5 (Nesterov label-noise contribution). Assume the stable small-step condition ηλ1 ⩽ cst , where cst > 0 is a sufficiently small universal constant. There exist constants C, c− , c+ > 0 such that if   1 1−ρ k ⩾ C max , , 1−ρ η then σ2 c− min B

(

η 1−ρ

1

β

η , 1−ρ

)

k−1

η2σ2 X σ2 ⩽ Knesterov (τ ) ⩽ c+ min B B

(

τ =0

η 1−ρ

1

β

η , 1−ρ

) .

Hence the Nesterov label-noise floor is k−1

η2σ2 X σ2 Knesterov (τ ) ≂ min B B τ =0

(

η 1−ρ

1

β

η , 1−ρ

) .

Proof. We first compute the saturated Green energy in one eigendirection. The Nesterov Green function satisfies Gnesterov (j) = (1+ρ)(1−ηλj )Gnesterov (j)−ρ(1−ηλj )Gnesterov (j), k+1 k k−1

Gnesterov (j) = 1, 0

Gnesterov (j) = 0. −1

Multiplying this recursion by Gnesterov (j), summing over k ⩾ 0, and then squaring the recursion k and summing over k ⩾ 0 gives the standard second-order energy identity ∞ X k=0

Gnesterov (j) k

2

= 1 − ρ(1 − ηλj )



1 + ρ(1 − ηλj ) .
2 1 + ρ(1 − ηλj ) − (1 + ρ)2 (1 − ηλj )2

The second factor in the denominator simplifies as
2 1 + ρ(1 − ηλj ) − (1 + ρ)2 (1 − ηλj )2 = ηλj (1 + (1 + 2ρ)(1 − ηλj )) . 51

Under ηλ1 ⩽ cst , the remaining factors are bounded above and below by constants except for 1 − ρ + ηλj and ηλj . Therefore ∞ X

Gnesterov (j) k

2

≂

k=0

1 . ηλj (1 − ρ + ηλj )

Consequently the saturated label-noise contribution is ∞

λj η 2 σ 2 X 2 X nesterov
2 ησ 2 X Gk (j) ≂ λj B B 1 − ρ + ηλj j⩾1

k=0

=

j⩾1 2 X σ

B

λj . λ + 1−ρ η j⩾1 j

Here Nesterov differs from Polyak: the effective root product is ρ(1 − ηλj ), so the saturated spectral weight contains λj /(λj + (1 − ρ)/η), rather than simply λj . We next estimate the spectral sum. If η/(1 − ρ) ⩽ 1, then (1 − ρ)/η ≳ 1, and hence λj η X η ≂ λj ≂ . 1−ρ 1−ρ 1−ρ λ + η j⩾1 j j⩾1

X

−β , the number of indices If η/(1 − ρ) ⩾ 1, then the spectral cutoff is λj ≂ 1−ρ η . Since λj ≂ j  1 β η above this cutoff is 1−ρ . The modes above the cutoff each contribute a constant amount to the summand, while the tail below the cutoff contributes the same order:

λj



λj ≂ min λ + 1−ρ η j⩾1 j

(

X

1−ρ ≂ 1−ρ η λj <

η 1−ρ

1

β

.

η

Thus, in all cases, X

η 1−ρ

1

β

η , 1−ρ

) .

It remains to pass from the saturated energy to the finite-time sum. The root formula for (j) gives a uniform tail estimate: after increasing C if necessary, the time condition Gnesterov k   1 1−ρ k ⩾ C max , 1−ρ η ensures that the finite partial sum captures a fixed positive fraction of the saturated Green energy on the spectral region responsible for the lower bounds above. When η/(1 − ρ) ⩽ 1, this region may be taken to be a fixed finite block of leading eigendirections, whose relaxation time is of order (1 − ρ)/η. When η/(1 − ρ) ⩾ 1, it may be taken to be the block λj ≳ (1 − ρ)/η, whose relaxation time is at most of order (1 − ρ)−1 . Therefore k−1 X τ =0

Knesterov (τ ) ≂

X

λ2j

j⩾1

∞ X

Gnesterov (j) τ

2

τ =0

up to constants for the purpose of the two-sided label-noise estimate. Combining this finite-time comparison with the saturated computation above proves the lemma.

52

E.4

Uniform boundedness of expected risk

As long as the additive label noise is strictly positive (σ > 0), we can prove that the expected risk is uniformly bounded, which allows the multiplicative gradient noise to be completely absorbed by the additive label noise floor. Lemma E.6 (Uniform boundedness of expected risk). Suppose σ > 0. Under the stability condition η ≲ 1 and η ⩽ cB β (1 − ρ) for a sufficiently small universal constant c > 0, the expected excess risk of Nesterov is uniformly bounded for all k ⩾ 0: E[E(θk )] ⩽ M, where M > 0 is a constant independent of k, η, B, and ρ. Proof. We proceed by induction on k. By Proposition E.3, there exist universal constants C1 , C2 , C3 > 0 such that for any k ⩾ 1: k−1

E[E(θk )] ⩽ C1 Snesterov (k) + C2

k−1

η2σ2 X η2 X Knesterov (τ ) + C3 Knesterov (τ )E[E(θk−τ )]. B B τ =0

(79)

τ =0

First, Propositions E.8 and E.10 establish that the deterministic bias monotonically decays or is uniformly bounded by its initialization across all regimes. Thus: Snesterov (k) ⩽ Ch2 E(θ0 ) =: M1 . Second, applying Lemma E.5, the infinite sum of the convolution kernel satisfies: ( ) 1 ∞ β η2 X CK η η Knesterov (τ ) ⩽ min , . B B 1−ρ 1−ρ τ =0

Therefore, the label noise contribution is universally bounded by: k−1

CK η2σ2 X Knesterov (τ ) ⩽ C2 σ 2 C2 B B



τ =0

η 1−ρ

1

β

.

1/β  η Under the stability condition η ⩽ cB β (1 − ρ), we have 1−ρ ⩽ c1/β B. Thus, the label noise is bounded by C2 σ 2 CK c1/β =: M2 . To handle the multiplicative noise term, we must first isolate the τ = 0 component, which contains the current step E[E(θk )]. Note that X
2 X 2 Knesterov (0) = λ2j Gnesterov (j) = λj ⩽ κ0 < ∞. 0 j⩾1

j⩾1 2

The τ = 0 term is therefore bounded by C3 ηB κ0 E[E(θk )]. Since the local stability condition 2 implies η ≲ 1, we can choose the learning rate prefactor sufficiently small such that C3 ηB κ0 ⩽ 12 . Subtracting this term from both sides of (79) gives: k−1

k−1

τ =0

τ =1

η2σ2 X 1 η2 X E[E(θk )] ⩽ C1 Snesterov (k) + C2 Knesterov (τ ) + C3 Knesterov (τ )E[E(θk−τ )]. (80) 2 B B

53

Now, we assume inductively that E[E(θm )] ⩽ M for all m < k. For the remaining multiplicative noise summation (τ ⩾ 1 =⇒ k − τ < k), we can apply the induction hypothesis: ! k−1 ∞ η2 X η2 X Knesterov (τ )E[E(θk−τ )] ⩽ C3 Knesterov (τ ) M C3 B B τ =1 τ =0  1 β 1 η ⩽ C3 CK M B 1−ρ ⩽ C3 CK c1/β M. By choosing c sufficiently small such that c1/β ⩽ 4C31CK , this term is limited to M/4. Substituting the bounds M1 , M2 , and M/4 into (80) yields: 1 M E[E(θk )] ⩽ M1 + M2 + . 2 4 Setting M = 4(M1 + M2 ) guarantees that 1 M M M E[E(θk )] ⩽ + = =⇒ E[E(θk )] ⩽ M. 2 4 4 2 The base case k = 0 holds trivially since E[E(θ0 )] = Snesterov (0) ⩽ M1 < M . This completes the induction. Corollary E.7 (Absorption of multiplicative noise). For σ > 0, under the stability condition η ≲ min{1, B β (1 − ρ)}, the multiplicative gradient noise is strictly bounded by the additive label noise up to a constant factor: ! k−1 k−1 η2 X M η2σ2 X Knesterov (τ )E[E(θk−τ )] ⩽ 2 Knesterov (τ ) . B σ B τ =0

τ =0

Consequently, the multiplicative noise can be absorbed, and the expected excess risk is asymptotically equivalent to the sum of the deterministic bias and the label noise: k−1

E[E(θk )] ≂ Snesterov (k) +

η2σ2 X Knesterov (τ ). B τ =0

E.5

Risk dynamics in the fully overdamped regime

In the fully overdamped regime, every characteristic root is real and the deterministic bias has no oscillatory component. Proposition E.8 (Deterministic bias, fully overdamped regime). Under Assumptions 3.1 and 3.3, in the fully overdamped regime ηλ1 ⩽ c0 (1 − ρ)2 with a sufficiently small universal constant c0 > 0, we have:  1, k ≲ 1−ρ η , Snesterov (k) ≂  1−ρ s 1−ρ  ηk , k≳ η . Proof. For each coordinate, write the deterministic multiplier as nesterov edet (j)θj∗ . k (j) = −hk

The deterministic Nesterov coordinate recursion gives hnesterov (j) = (1+ρ)(1−ηλj )hnesterov (j)−ρ(1−ηλj )hnesterov (j), k+1 k k−1 54

hnesterov (j) = 1, 0

hnesterov (j) = 1−ηλj . 1

Its characteristic roots are yj,+ , yj,− =

(1 + ρ)(1 − ηλj ) ±

p (1 + ρ)2 (1 − ηλj )2 − 4ρ(1 − ηλj ) . 2

In the fully overdamped regime, ηλj ⩽ c0 (1 − ρ)2 for every j. To locate the two roots, evaluate the characteristic polynomial at 1 − u. A direct expansion gives (1 − u)2 − (1 + ρ)(1 − ηλj )(1 − u) + ρ(1 − ηλj )
= ηλj − 1 − ρ + (1 + ρ)ηλj u + u2 . Taking c0 > 0 sufficiently small, the two evaluations ηλj + 2(1 − ρ)



2ηλj ηλj − 1 − ρ + (1 + ρ)ηλj + 1−ρ



ηλj − 1 − ρ + (1 + ρ)ηλj

ηλj 2(1 − ρ)

and 2ηλj 1−ρ

2 >0

2 <0

hold uniformly for all 0 < ηλj ⩽ c0 (1 − ρ)2 . Hence the slow root satisfies 1−

2ηλj ηλj ⩽ yj,+ ⩽ 1 − . 1−ρ 2(1 − ρ)

Moreover, yj,+ yj,− = ρ(1 − ηλj ). Using the lower bound on yj,+ , and decreasing c0 if necessary, the fast root obeys 0 ⩽ yj,− ⩽ 1 −

1−ρ . 2

The solution of the two-point initial value problem is hnesterov (j) = k

1 − ηλj − yj,− k yj,+ − (1 − ηλj ) k yj,+ + yj,− . yj,+ − yj,− yj,+ − yj,−

The root bounds above imply 1 − ηλj − yj,− 1 ⩽ ⩽ 2, 2 yj,+ − yj,−

yj,+ − (1 − ηλj ) ⩽ Cc0 . yj,+ − yj,−

Indeed, yj,+ − yj,− ≂ 1 − ρ, while 1 − ηλj − yj,− ≂ 1 − ρ, and |yj,+ − (1 − ηλj )| ≲

ηλj ⩽ c0 (1 − ρ). 1−ρ

Thus the slow component has a uniformly positive coefficient and the fast component has a coefficient of order c0 .
2 P To evaluate the deterministic bias Snesterov (k) = 12 j⩾1 λj (θj∗ )2 hnesterov (j) , we divide k the time k into three successive stages. 1 . Stage 1: k ≲ 1−ρ In this extremely short-time regime, for the leading modes (j = O(1)), since ηλj ⩽ c0 (1 − ρ)2 , 1 the slow root satisfies yj,+ ⩾ 1 − 2c0 (1 − ρ). For k ≲ 1−ρ , we have C

k ⩾ (1 − 2c0 (1 − ρ)) 1−ρ ≂ 1. yj,+

55

Since the dynamics are unconditionally stable (hnesterov (j) ⩽ 1), it follows that hnesterov (j) ≂ 1 k k for the leading modes. Consequently, the bias is dominated by these non-decayed modes: X Snesterov (k) ≂ λj (θj∗ )2 ≂ 1. j⩾1 1 Stage 2: 1−ρ ≲ k ≲ 1−ρ η .

k  1 k ⩽ 1 − 1−ρ , the fast root decays exponentially since yj,− ≪ 1. The multiplier is For k ≳ 1−ρ 2 dominated by the slow root: hnesterov (j) ≃ k

1 − ηλj − yj,− k k y ≃ yj,+ . yj,+ − yj,− j,+

For the leading modes with λj = Θ(1), within the time window k ≲ 1−ρ η , the exponent is bounded by:     ηλj ηλj k k k ≂ 1, ≂ exp −c yj,+ ≂ 1 − c 1−ρ 1−ρ ηλ

where we used the condition 1−ρj k ≲ λj ⩽ λ1 = O(1). Thus, the leading modes have not (j) ≂ 1. As a result, the yet experienced significant exponential decay, preserving hnesterov k deterministic bias remains unchanged: Snesterov (k) ≂ 1. Stage 3: k ≳ 1−ρ η . In this long-time regime, the exponential decay of the slow root becomes significant across the spectrum. For the slow root, the preceding bounds give     ηλj ηλj k exp −C k ⩽ yj,+ ⩽ exp −c k . 1−ρ 1−ρ For the fast root, since ηλj /(1 − ρ) ⩽ c0 (1 − ρ), the ratio of the fast component to the slow component is bounded by   ηλj k ⩽ Cc0 exp (−c(1 − ρ)k) , Cc0 exp −c(1 − ρ)k + C 1−ρ after decreasing c0 . Therefore, for all k ⩾ C(1 − ρ)−1 , the fast component is at most a fixed fraction of the slow component. Consequently there exist constants c, C > 0 such that     ηλj ηλj nesterov c exp −C k ⩽ hk (j) ⩽ C exp −c k . 1−ρ 1−ρ Substituting this coordinate estimate into the deterministic bias gives     X X ηλj ηλj ∗ 2 ∗ 2 k ⩽ Snesterov (k) ⩽ C λj (θj ) exp −c k . c λj (θj ) exp −C 1−ρ 1−ρ j⩾1

j⩾1

By Assumption 3.3, λj (θj∗ )2 ≂ j −1−βs , It remains to estimate X j⩾1

j

−1−βs



λj ≂ j −β .

 ηk −β exp −c j . 1−ρ

56

ηk The condition k ≳ 1−ρ η gives 1−ρ ⩾ C.

For the lower bound, take the tail j ⩾



ηk 1−ρ

1/β 

. On this tail the exponential factor is

bounded below by a positive constant, hence     X ηk −s ηk −β −1−βs ⩾c j . j exp −C 1−ρ 1−ρ j⩾1

For the upper bound, an integral comparison gives     Z ∞ X ηk −β ηk −β −1−βs −1−βs x exp −c j exp −c j ⩽C x dx 1−ρ 1−ρ 1 j⩾1   ηk −s . ⩽C 1−ρ Combining the spectral upper and lower bounds yields   1−ρ s Snesterov (k) ≂ . ηk Synthesizing the three stages completes the proof. Theorem E.9 (Scaling law, fully overdamped regime). Under Assumptions 3.1 and 3.3, assume the stability condition η ≲ min{1, B β (1−ρ)} holds. For the fully overdamped regime η ≲ (1−ρ)2 with k ≳ 1−ρ η , the expected excess risk satisfies:  EEk ≂

1−ρ ηk

s +

σ2η . B(1 − ρ)

Proof. By Corollary E.7, the expected risk is completely determined by the deterministic bias and the label noise: k−1 η2σ2 X E[E(θk )] ≂ Snesterov (k) + Knesterov (τ ). B τ =0

For k ≳ 1−ρ η , we substitute the deterministic bias from Proposition E.8:  Snesterov (k) ≂

1−ρ ηk

s .

We substitute the label noise floor from Lemma E.5. In the fully overdamped regime η ≲ (1−ρ)2 , η we have 1−ρ ≲ 1 − ρ ≪ 1. Since β > 1, the minimum is attained by the linear term: k−1

η2σ2 X σ2 Knesterov (τ ) ≂ min B B

(

τ =0

η 1−ρ

1

β

η , 1−ρ

) =

σ2η . B(1 − ρ)

Summing these two terms directly yields the exact scaling law.

E.6

Risk dynamics in the mixed-damping regime

In the mixed-damping regime, the spectrum is partitioned into underdamped and overdamped  1/β η modes by the crossover index Jmax ≂ (1−ρ) . Consequently, the deterministic bias 2 decomposes into an overdamped polynomial tail and an underdamped oscillatory transient, denoted by Sover (k) and Sunder (k). Due to destructive interference among the underdamped modes, Sunder (k) lacks a strictly positive pointwise lower bound. 57

Proposition E.10 (Deterministic bias, mixed-damping regime). Under Assumptions 3.1 and 3.3, for the mixed-damping regime η ≳ (1 − ρ)2 , the deterministic bias decomposes into Snesterov (k) = Sover (k) + Sunder (k), where the overdamped tail satisfies the two-sided bounds:   2 s 1  (1−ρ) , k ≲ 1−ρ , η Sover (k) ≂  1−ρ s 1  , k ≳ 1−ρ , ηk and Sunder (k) represents the underdamped transient whose properties are formally established in Lemma E.11. nesterov (j)θ ∗ . The multiplier satisfies Proof. For each coordinate, write edet j k (j) = −hk

(j) = (1+ρ)(1−ηλj )hnesterov (j)−ρ(1−ηλj )hnesterov (j), hnesterov k+1 k k−1

hnesterov (j) = 1, hnesterov (j) = 1−ηλj . 0 1

The roots are yj,+ , yj,− =

(1 + ρ)(1 − ηλj ) ±

p (1 + ρ)2 (1 − ηλj )2 − 4ρ(1 − ηλj ) . 2

The discriminant is negative precisely when ηλj >

(1 − ρ)2 . (1 + ρ)2

Thus the comparison between η and (1 − ρ)2 determines whether a coordinate is overdamped or underdamped. In the mixed-damping regime η ≳ (1 − ρ)2 , the crossover index is defined by ηλJmax ≂ (1 − ρ)2 . We decompose the deterministic bias into the overdamped tail (j > Jmax ) and the underdamped leading block (j ⩽ Jmax ). For overdamped coordinates in the stable small-step range, the root comparison used in Proposition E.8 gives constants c1 , C1 > 0 such that   ηλj nesterov (j)| ⩽ C1 exp −c1 |hk k . 1−ρ Moreover, if ηλj ⩽ c1 (1 − ρ)2 and k ⩾ C1 (1 − ρ)−1 , the slow root dominates the fast root and   ηλj nesterov |hk (j)| ⩾ c1 exp −C1 k . 1−ρ 2

ηλ

1 1 , the exponent satisfies 1−ρj k ≲ (1−ρ) For k ≲ 1−ρ 1−ρ 1−ρ = 1. The exponential decay is bounded away from zero, so |hnesterov (j)| ≂ 1. The sum is dominated by its lower limit at the k crossover:  s X X (1 − ρ)2 2 −1−βs −βs Sover (k) = λj (θj∗ )2 hnesterov (j) ≂ j ≂ (J ) ≂ . max k η j>Jmax

j>Jmax

1 For k ≳ 1−ρ , applying λj ≂ j −β , the integral comparison yields the strict upper bound:

Sover (k) ⩽ C

X j>Jmax

j

−1−βs



ηk −β exp −2c1 j 1−ρ



 ⩽C

1−ρ ηk

s .

 s ηk 1/β By restricting the sum to j ⩾ (C 1−ρ ) , we obtain the matching lower bound c 1−ρ . ηk The remaining deterministic bias comes from the underdamped modes, defined as P Sunder (k) := 12 j⩽Jmax λj (θj∗ )2 hnesterov (j)2 , which is analyzed in Lemma E.11. k 58

Lemma E.11 (Properties of the underdamped transient). Let Sunder (k) denote the deterministic bias contributed by the underdamped modes (j ⩽ Jmax ). Under the conditions of Proposition E.10, define the piecewise envelope Eenv (k) by:  ρk , k ≲ η1 ,        1 1 ρk   , ≲k≲ ,  s  η (1 − ρ)2  (ηk)  s−1 Eenv (k) := 1 (1 − ρ)2 1 −2− β2 k   , pnesterov , ≲ k ≲ η 1/β (1 − ρ)  2  ηk η (1 − ρ)    s+ β1    (1 − ρ)2 −2− β2  k  , pnesterov , k ≳ η 1/β (1 − ρ) η   (1−ρ)2 4ρ2 where pnesterov := ρ 1 − (1+ρ) = (1+ρ) 2 2 < ρ. Then for every step k ⩾ 0, Sunder (k) satisfies: 0 ⩽ Sunder (k) ⩽ CEenv (k). Proof. For j ⩽ Jmax , write the two characteristic roots as: p q q (1 + ρ) 1 − ηλj iωj −iωj . ρ(1 − ηλj )e , ρ(1 − ηλj )e , where cos ωj = √ 2 ρ Solving the initial conditions gives the exact identity for the coordinate multiplier: " # p
(1 − ρ) 1 − ηλ j k/2 (j) = ρ(1 − ηλj ) cos(kωj ) + hnesterov sin(kωj ) . √ k 2 ρ sin ωj √ (1−ρ) 1−ηλj Set Aj := 2√ρ sin ωj and define tan ψj = Aj . The squared multiplier is bounded by the envelope:
2
k (j) = ρ(1 − ηλj ) (1 + A2j ) cos2 (kωj − ψj ). hnesterov k 2

On the separated underdamped band ηλj ⩾ (1−ρ) , we have 4ρ sin2 ωj ≂ ηλj , which guarantees (1+ρ)2 |Aj | ⩽ C. Defining the positive spectral weights:
k 1 wj (k) := λj (θj∗ )2 ρ(1 − ηλj ) (1 + A2j ), 2 the total underdamped bias decomposes into: X Sunder (k) = wj (k) cos2 (kωj − ψj ). j⩽Jmax

Evaluation of the static mass envelope Mw (k). We define the static mass: J

Mw (k) :=

max
k 1 X λj (θj∗ )2 ρ(1 − ηλj ) (1 + A2j ). 2

j=1

Since |Aj | ⩽ C, we have 1 + A2j ≂ 1. Substituting λj (θj∗ )2 ≂ j −1−βs and λj ≂ j −β , the static mass is asymptotically equivalent to: Mw (k) ≂ ρk

JX max

j −1−βs (1 − ηλj )k .

j=1

59

We evaluate Mw (k) across four distinct temporal stages: Stage 1 (k ≲ η1 ): For all modes j ∈ [1, Jmax ], we have ηλj ⩽ ηλ1 ≲ η. Since k ≲ η1 , the product satisfies kηλj ≲ 1, which implies: (1 − ηλj )k ≂ 1,

∀ j ∈ [1, Jmax ].

Consequently, the summation is dominated by the leading indices (j ∼ 1): Mw (k) ≂ ρk

JX max

j −1−βs ≂ ρk

j=1

∞ X

j −1−βs ≂ ρk .

j=1

1 k −β suppresses Stage 2 ( η1 ≲ k ≲ (1−ρ) 2 ): In this regime, kη ≳ 1, and (1 − ηλj ) ≂ exp −ckηj the leading modes. The dominant spectral mass is concentrated at the continuous peak jk ≂ 1 (ηk)1/β . Since k ≲ (1−ρ) 2 , this peak lies strictly within the summation range: jk ≂ (ηk)1/β ≲



η (1 − ρ)2

1/β ≂ Jmax .

Applying Euler–Maclaurin integral comparison with the change of variables u = ckηx−β : Mw (k) ≂ ρk

Z Jmax 1

ρk = β(ckη)

x−1−βs exp −ckηx−β dx Z ckη us−1 e−u du. s −β ckηJmax

−β Since the integration limits satisfy ckηJmax ≂ k(1−ρ)2 ≲ 1 and ckη ≳ 1, the integral is bounded above and below by universal positive constants:

Z ckη −β ckηJmax

s−1 −u

u

e

Z ∞ du ≂

us−1 e−u du = Γ(s) ≂ 1.

0

Therefore, we obtain: Mw (k) ≂

ρk . (ηk)s

1 1 1/β moves past J Stage 3 & 4 (k ≳ (1−ρ) max . 2 ): When k ≳ (1−ρ)2 , the continuous peak (ηk)

Inside the summation interval j ⩽ Jmax , the decay factor (1 − ηλj )k is slowest at the boundary 2 j = Jmax , where ηλJmax = (1−ρ) =: uc . We define the boundary decay factor: (1+ρ)2 (1 − ρ)2 pnesterov := ρ(1 − uc ) = ρ 1 − (1 + ρ)2 

 =

4ρ2 . (1 + ρ)2

Let m := Jmax − j ⩾ 0. Linearizing the eigenvalues near the boundary gives:
−β−2 2 ηλj − uc = η(λJmax −m − λJmax ) = ∆x m + O ηJmax m , where the discrete frequency spacing is: −β−1 ∆x := ηβJmax ≂η



(1 − ρ)2 η

60

1+ β1 = (1 − ρ)

2+ β2 − β1

η

.



−1−βs Factoring out the boundary coefficient Jmax

≂

(1−ρ)2 η

s+ 1

β

, the sum is governed by the

geometric series over m:  Mw (k) ≂  =

(1 − ρ)2 η (1 − ρ) η

s+ β1

pknesterov

Jmax X−1

exp −ck∆x m

m=0

 1 2 s+ β

pknesterov

1 − exp −ck∆x Jmax
1 − exp −ck∆x

! .

Since k∆x Jmax ≂ kuc ≂ k(1 − ρ)2 ≳ 1, the numerator satisfies 1 − exp −ck∆x Jmax ≂ 1, yielding: !  s+ β1 (1 − ρ)2 1
. Mw (k) ≂ pknesterov (81) η 1 − exp −ck∆x 1 1/β (1 − ρ) Stage 3 ( (1−ρ) 2 ≲ k ≲ η

−2− β2

2+ β2 − β1

(81) is small: k∆x = k(1−ρ)

η

 Mw (k) ≂

): In this regime, the exponent in the denominator of

≲ 1. Using the Taylor expansion 1−e−u ≂ u for u ∈ (0, 1]: s+ β1

(1 − ρ)2 η

pknesterov



1 ck∆x



! 1 s+ β1 β η (1 − ρ)2 pknesterov ≂ 2+ 2 η k(1 − ρ) β  s−1 1 (1 − ρ)2 k = pnesterov . ηk η 

2+ 2 − 1

−2− 2

β ): Here, the exponent is large: k∆ β η β ≳ 1. Stage 4 (k ≳ η 1/β (1 − ρ) x = k(1
− ρ) Thus, the denominator in (81) is bounded away from zero: 1 − exp −ck∆x ≂ 1. The discrete mode j = Jmax (m = 0) dominates the sum:

 Mw (k) ≂

(1 − ρ)2 η

s+ β1

k pnesterov .

Combining the four stages establishes that Mw (k) ≂ Eenv (k) uniformly across all k ⩾ 0. Since cos2 (·) ⩽ 1, we immediately obtain the pointwise upper bound: X Sunder (k) ⩽ wj (k) = Mw (k) ⩽ CEenv (k). j⩽Jmax

The lower bound is trivially zero due to the potential for destructive interference of the oscillatory terms at specific iterations. Theorem E.12 (Scaling law, mixed-damping regime). Under Assumptions 3.1 and 3.3, assuming the stability condition η o≲ min{1, B β (1 − ρ)}, for the mixed-damping regime n 1 η ≳ (1 − ρ)2 with k ≳ max 1−ρ , 1−ρ , the expected excess risk satisfies: η  E[E(θk )] ≂ Sunder (k) +

1−ρ ηk

s

σ2 + min B

(

η 1−ρ

1

β

η , 1−ρ

) ,

(82)

where the properties of the underdamped transient Sunder (k) are established in Lemma E.11. 61

Proof. By Corollary E.7, the multiplicative noise is absorbed by the additive label noise floor, yielding: k−1 η2σ2 X E[E(θk )] ≂ Sunder (k) + Sover (k) + Knesterov (τ ). (83) B τ =0 n o 1−ρ 1 By Lemma E.5, for k ≳ max 1−ρ , η , the label noise summation captures the saturated floor: ) ( 1 k−1 β σ2 η η η2σ2 X Knesterov (τ ) ≂ . min , B B 1−ρ 1−ρ τ =0

1 By Proposition E.10, since k ≳ 1−ρ , the overdamped tail reaches its terminal polynomial decay phase:   1−ρ s Sover (k) ≂ . ηk Substituting these components explicitly into (83) establishes the equivalence (82).

F

Optimal data scaling across batch-size regimes

This appendix proves the batch-size-dependent scaling laws stated in Theorem 6.3. Throughout this appendix, we work in the noisy setting σ > 0 considered in Section 6 and impose the data-budget constraint D B ≂ Db , T = ≂ D1−b , 0 ⩽ b < 1. B Thus, the batch exponent b is prescribed, and only the remaining hyperparameters are optimized. We use the transition exponents defined in Section 6: s 2s + 1/β 2sβ + 1 2s + 1 b1 = , b2 = = , b3 = . s+1 2s + 1 + 1/β 2sβ + β + 1 2s + 2 For the momentum methods, we parameterize the momentum gap and learning rate by 1 − ρ ≂ D−m and η ≂ D−c , where m, c ⩾ 0.

F.1

SGD

For SGD, write η ≂ D−c , where c ⩾ 0. By the SGD risk scaling law in Section 5, σ2η ≂ D−s(1−b−c) + D−(b+c) . B Therefore, for a fixed batch exponent b, the achievable decay exponent is min{s(1 − b − c), b + c}. If 0 ⩽ b < b1 , the two terms can be balanced by taking c = b1 − b ⩾ 0. Since b1 = s/(s + 1), this gives s s(1 − b − c) = b + c = b1 = . s+1 Hence rsgd (b) = s/(s + 1) for 0 ⩽ b < b1 . If b ⩾ b1 , the unconstrained balancing choice would require c = b1 − b ⩽ 0. Because the learning rate is not allowed to grow polynomially with D, we must have c ⩾ 0. Moreover, b + c ⩾ b ⩾ b1 , so the signal term is the slower-decaying term. It is optimized by taking c = 0, which yields rsgd (b) = s(1 − b) for b1 ⩽ b < 1. Consequently,   s , 0 ⩽ b < b1 , rsgd (b) = s + 1 s(1 − b), b ⩽ b < 1. E[E(θT )] ≂ (ηT )−s +

1

Thus, b1 is precisely the largest batch exponent for which SGD preserves the data-optimal rate s/(s + 1). 62

F.2

Polyak

We next prove the Polyak part of Theorem 6.3. We retain T ≂ D1−b , 1−ρ
≂ D−m , and η ≂ D−c
, where m, c ⩾ 0. We require b + m < 1, so that ρT ⩽ exp −(1 − ρ)T = exp −Θ(D1−b−m ) is smaller than any polynomial power of D. The stability condition η ≲ B(1 − ρ) becomes b + c − m ⩾ 0. Moreover, ηρ T =

ηρ η = ≂ D−(b+c−m) . B B(1 − ρ)

ηT ≂ D1−b−c+m , 1−ρ

The damping threshold η ≂ (1 − ρ)2 corresponds to c = 2m. We consider the fully overdamped and mixed-damping regimes separately. Fully overdamped regime.

If c ⩾ 2m, Theorem 5.1 gives

E[E(θT )] ≂ (ηρ T )−s +

σ 2 ηρ ≂ D−s(1−b−c+m) + D−(b+c−m) . B

Therefore, the polynomial decay exponent in the fully overdamped regime is min{s(1 − b − c + m), b + c − m}. Mixed-damping regime.

If c ⩽ 2m, Theorem 5.1 gives

σ 2 ηρ . B

Since b + m < 1, the transient P (T ) ≲ ρT ⩽ exp −(1 − ρ)T = exp −Θ(D1−b−m ) is smaller than any polynomial power of D. Hence E[E(θT )] ≂ P (T ) + (ηρ T )−s +

E[E(θT )] ≂ D−s(1−b−c+m) + D−(b+c−m) , and the polynomial decay exponent is again min{s(1 − b − c + m), b + c − m}. Thus, although the two regimes have different transient dynamics, they lead to the same polynomial optimization problem. The signal and additive-noise powers are balanced when s(1 − b − c + m) = b + c − m

⇐⇒

b+c−m=

s = b1 . s+1

If 0 ⩽ b < b1 , this balance is attained by taking m = 0 and c = b1 − b. Since c ⩾ 2m, this choice lies in the fully overdamped regime and yields rpolyak (b) = s/(s + 1). If b1 ⩽ b < b3 , the same balance is attained by taking c = 0 and m = b − b1 . For b > b1 , this choice satisfies c ⩽ 2m and therefore lies in the mixed-damping regime. Moreover, b+c−m = b1 , while the transient condition becomes b+m = 2b−b1 < 1. Since (1+b1 )/2 = (2s+1)/(2s+2) = b3 , this condition holds precisely when b < b3 . Hence rpolyak (b) = s/(s + 1) for 0 ⩽ b < b3 . Now suppose that b ⩾ b3 . Since c ⩾ 0 and m < 1 − b, every admissible choice satisfies b + c − m > 2b − 1 ⩾ b1 , so the signal term is the slower-decaying term and s(1 − b − c + m) < 2s(1 − b). This upper bound is approached by taking c = 0 and m ↑ 1 − b, which lies in the mixed-damping regime and gives s(1 − b − c + m) ↑ 2s(1 − b),

b + c − m ↓ 2b − 1.

Moreover, 2b − 1 ⩾ 2s(1 − b) holds exactly when b ⩾ b3 , so the signal term is indeed dominant. Consequently, rpolyak (b) = 2s(1 − b) for b3 ⩽ b < 1, where the endpoint value is understood in the supremal sense represented by the o(1) term in Definition 6.2. Combining the two regimes gives   s , 0 ⩽ b < b3 , rpolyak (b) = s + 1 2s(1 − b), b ⩽ b < 1. 3

63

Thus, Polyak preserves the optimal data-scaling exponent s/(s + 1) up to the critical batch exponent b3 . Compared with SGD, it extends the critical batch size from Db1 to Db3 , while leaving the optimal small-batch exponent unchanged.

F.3

Nesterov

We finally prove the Nesterov part of Theorem 6.3. We retain B ≂ Db , T ≂ D1−b , 1−ρ ≂ D−m , and η ≂ D−c , where m, c ⩾ 0. We require b + m < 1, so that the Nesterov transient N (T ) is smaller than any polynomial power of D. The stability condition η ≲ B β (1 − ρ) becomes βb + c − m ⩾ 0, while the damping threshold η ≂ (1 − ρ)2 corresponds to c = 2m. We optimize the fully overdamped and mixed-damping regimes separately and then compare their optimal decay exponents. Fully overdamped regime.

If c ⩾ 2m, Theorem 5.2 gives

E[E(θT )] ≂ (ηρ T )−s +

σ 2 ηρ ≂ D−s(1−b−c+m) + D−(b+c−m) . B

Thus, the decay exponent is min{s(1 − b − c + m), b + c − m}. Since c ⩾ 2m implies c − m ⩾ 0, we have b + c − m ⩾ b. If 0 ⩽ b < b1 , the signal and additive-noise powers can be balanced by taking m = 0 and c = b1 − b, which gives b + c − m = b1 and decay exponent s/(s + 1). If b ⩾ b1 , then b + c − m ⩾ b ⩾ b1 , so the signal term is limiting. Its exponent is maximized by taking m = c = 0, which gives s(1 − b). Therefore, the optimal exponent in the fully overdamped regime is   s , 0 ⩽ b < b1 , full rnesterov (b) = s + 1 s(1 − b), b ⩽ b < 1. 1

In particular, the fully overdamped Nesterov dynamics have the same batch-size-dependent scaling exponent as SGD. Mixed-damping regime.

If c ⩽ 2m, Theorem 5.2 gives

E[E(θT )] ≂ N (T ) + (ηρ T )−s +

n o σ2 min ηρ , ηρ1/β . B

Since b + m <
1, the first term N (T ) ≲ min{1, (ηT )−s }ρT ⩽ exp −(1 − ρ)T exp −Θ(D1−b−m ) is smaller than any polynomial power of D. Hence (  1/β ) 2 σ η η E[E(θT )] ≂ D−s(1−b−c+m) + min , . B 1−ρ 1−ρ

=

The signal term has decay power s(1 − b − c + m). Since η/(1 − ρ) ≂ Dm−c , the additive-noise decay power is  b + c − m, c ⩾ m, b + c − m , c ⩽ m. β We optimize these two branches separately.

64

The branch c ⩾ m.

In this branch, the decay exponent is min{s(1 − b − c + m), b + c − m}.

Because c−m ⩾ 0, the same argument as in the fully overdamped regime applies. For 0 ⩽ b < b1 , the balance b + c − m = b1 is feasible in the mixed-damping regime. For example, one may take m = b1 − b and c = 2(b1 − b), for which b + m = b1 < 1. For b ⩾ b1 , the optimal choice satisfies c − m = 0 and gives exponent s(1 − b). Therefore,   s , 0 ⩽ b < b1 , mixed, c⩾m rnesterov (b) = s + 1 s(1 − b), b ⩽ b < 1. 1

The branch c ⩽ m.

In this branch, the decay exponent is   c−m . min s(1 − b − c + m), b + β

If 0 ⩽ b < b1 , then at m − c = 0 we have b < s(1 − b). Increasing m − c increases the signal power but decreases the additive-noise power, so the optimal choice is m − c = 0, which gives decay exponent b. If b ⩾ b1 , balancing the signal and additive-noise powers gives
β (s + 1)b − s c−m s(1 − b − c + m) = b + , m−c= . β sβ + 1 For a fixed value of m − c, the least restrictive choice for the transient condition is c = 0, so we take
β (s + 1)b − s c = 0, m= . sβ + 1 This choice satisfies the stability condition because
βs 1 + (β − 1)b βb + c − m = > 0. sβ + 1 The signal and additive-noise powers are then equal to
s 1 + (β − 1)b m s(1 − b + m) = b − = . β sβ + 1 The transient condition holds precisely when b+m<1

⇐⇒

b<

2sβ + 1 2s + 1/β = = b2 . 2sβ + β + 1 2s + 1 + 1/β

Therefore, for b1 ⩽ b < b2 , the optimal exponent on the branch c ⩽ m is
s 1 + (β − 1)b . sβ + 1 If b ⩾ b2 , the preceding balance is no longer compatible with b + m < 1. Since c ⩾ 0 and m < 1 − b, every admissible choice satisfies s(1 − b − c + m) < 2s(1 − b). This upper bound is approached by taking c = 0 and m ↑ 1 − b. Under this choice, s(1 − b + m) ↑ 2s(1 − b),

b− 65

m (β + 1)b − 1 → . β β

At b = b2 , the two limiting powers agree, while for b > b2 , (β + 1)b − 1 > 2s(1 − b). β Thus, the signal term is limiting and the supremal decay exponent is 2s(1 − b). Consequently,   b, 0 ⩽ b < b1 ,   
 s 1 + (β − 1)b mixed, c⩽m rnesterov (b) = , b1 ⩽ b < b2 ,  sβ + 1     2s(1 − b), b2 ⩽ b < 1. Taking the larger exponent between the two additive-noise branches gives the optimal exponent in the mixed-damping regime:  s  , 0 ⩽ b < b1 ,  s + 1 
 mixed rnesterov (b) = s 1 + (β − 1)b , b1 ⩽ b < b2 ,  sβ + 1     2s(1 − b), b2 ⩽ b < 1. Indeed, the branch c ⩾ m is optimal for b < b1 , whereas the branch c ⩽ m is optimal for b > b1 . Finally, optimizing over the fully overdamped and mixed-damping regimes gives n o full mixed rnesterov (b) = max rnesterov (b), rnesterov (b) . For 0 ⩽ b < b1 , the two regimes attain the same exponent s/(s + 1). For b1 < b < b2 , the mixed-damping regime exponent s(1 + (β − 1)b)/(sβ + 1) is strictly larger than the fully overdamped exponent s(1 − b). For b ⩾ b2 , the mixed-damping regime exponent 2s(1 − b) is also larger than the fully overdamped exponent s(1 − b). We therefore obtain  s  , 0 ⩽ b < b1 ,   s + 1
 rnesterov (b) = s 1 + (β − 1)b , b1 ⩽ b < b2 ,  sβ + 1     2s(1 − b), b2 ⩽ b < 1. Thus, Nesterov agrees with SGD and Polyak in the small-batch regime, improves its data-scaling exponent once b exceeds b1 , and reaches its maximal exponent at b = b2 . For b ⩾ b2 , its exponent becomes 2s(1 − b), which coincides with the Polyak exponent in the ultra-large-batch regime.

G

Experimental details and additional results

Our theoretical analysis is formulated in the infinite-dimensional PLK setting. For numerical experiments, we truncate the system to the first d = 4096 eigendirections, so that, for example, H = diag(λ1 , . . . , λd ) ∈ Rd×d . Unless otherwise stated, all experiments use d = 4096.

G.1

Risk stability boundaries

Simulation setup. We run with θ−1 = θ0 = 0 for T = 2000 iterations. A run is classified as stable if its average risk over the final 20 iterations is below 1; it is classified as unstable if the risk blows up or exceeds 1012 . For each setting and random seed, the empirical stability boundary is the largest stable learning rate on the search grid. We repeat each scan over five independent seeds and report the mean boundary. 66

• Stability versus batch size. For Figure 2 (left), we use (s, β) = (3, 2), fix 1 − ρ = 2−24 for the momentum methods, and scan B ∈ {20 , 21 , . . . , 212 } and 50 logarithmically spaced learning rates between 2−10 and 8. The dashed curves are post-hoc fits to the predicted scalings: a constant for SGD, B(1 − ρ) for Polyak, and B β (1 − ρ) followed by a constant plateau for Nesterov. The dotted line marks the split used for the Nesterov fit. • Stability versus momentum damping. For Figure 2 (middle), we use (s, β) = (3, 2), σ = 0, and B = 32, vary 1 − ρ ∈ {2−24 , 2−22 , . . . , 2−2 }, and scan 50 logarithmically spaced learning rates between 2−22 and 8. Since SGD is independent of ρ, its boundary is constant across the damping grid. The dashed curves are post-hoc fits to the predicted form η crit = min{L, c(1 − ρ)} for Polyak and Nesterov, and to a constant for SGD.

100

10−1

10−2

10−3

SGD Prediction Polyak Prediction Nesterov Prediction

22

23

24

25

26

27

Batch Size B

28

29

210

(a) s = 0.5, β = 1.5

100

10−1

10−2

10−3

SGD Prediction Polyak Prediction Nesterov Prediction

22

23

24

25

26

27

Batch Size B

28

29

210

Stability Boundary η

Figure 4 shows additional results for different choices of (s, β).

Stability Boundary η

Stability Boundary η

Additional results.

100

10−1

10−2

10−3

(b) s = 0.2, β = 2

SGD Prediction Polyak Prediction Nesterov Prediction

22

23

24

25

26

27

Batch Size B

28

29

210

(c) s = 0.5, β = 2.5

Figure 4: Additional stability boundaries across (s, β).

G.2

Risk-dynamics scaling laws

Simulation setup. We simulate PLK regression initialized with θ−1 = θ0 = 0. For Figure 2 (right), we set (s, β) = (0.3, 4), σ = 1, B = 32, η = 0.03, and ρ = 0.98. Each method is run for 2 × 107 iterations. Let E(k) denote the population excess risk averaged over the 100 runs. Scaling-law fits. Theorems 5.1 and 5.2 determine the functional forms of the risk dynamics up to multiplicative constants. We therefore introduce Θ = (C1 , C2 , C3 ) ∈ R3+ and fit the following ansatz separately for each method. For Polyak, Fpolyak (k; Θ) = C1 ρk + C2 (ηρ k)−s + C3

σ2 ηρ , B

ηρ :=

η , 1−ρ

while for Nesterov, Fnesterov (k; Θ) = C1 min{1, (ηk)−s }ρk + C2 (ηρ k)−s + C3

σ2 min{ηρ , ηρ1/β }. B

The first terms use the transient envelopes from Theorems 5.1 and 5.2, so the fits test the predicted scaling forms rather than their multiplicative constants. For each method, the coefficients are fitted by minimizing the mean squared relative error: 1 X Θ∈R3+ |Ifit | min

k∈Ifit



F (k; Θ) − E(k) E(k)

2 ,

Ifit := {k ⩾ 1 : ηρ k ⩾ 1},

where F = Fpolyak or Fnesterov for the corresponding method.

67

Additional results. We repeat the same experiment for (s, β) ∈ {(0.3, 6), (0.5, 4), (0.5, 6)}, keeping σ = 1, B = 32, η = 0.03, and ρ = 0.98. Each run uses 5 × 105 iterations, and the reported risk is averaged over 100 independent runs. Figure 5 shows that the predicted scaling forms consistently track the empirical risk dynamics across these choices of (s, β). Scaling law predictions

Scaling law predictions 100

Polyak Nesterov

10−1

100

101

102

104

103

105

Polyak Nesterov

10−1

10−2

10−2

10−2

Scaling law predictions

100

Polyak Nesterov

Loss

10−1

Loss

Loss

100

101

100

102

Step

104

103

105

100

101

102

Step

(a) (s, β) = (0.3, 6)

104

103

105

Step

(b) (s, β) = (0.5, 4)

(c) (s, β) = (0.5, 6)

Figure 5: Additional validation of the risk-dynamics scaling laws for Polyak and Nesterov. Thin translucent lines show empirical excess risks, while thick solid curves show the fitted scaling-law predictions. All experiments use σ = 1, B = 32, η = 0.03, and ρ = 0.98, run for 5 × 105 iterations, and report averages over 100 independent runs.

G.3

Batch-size phase diagram at a fixed data budget

We simulate PLK regression with (s, β) = (0.3, 4) and σ = 0.1, with all methods initialized at θ−1 = θ0 = 0. We fix the total data budget at D = 214 . For each batch size, we search locally around the hyperparameter scalings predicted by Theorem 6.3. Each candidate is evaluated over 50 independent runs. Let W = max{1, ⌈0.01T ⌉}. We score each candidate by the mean excess risk over the final W iterations and across runs: 50

T X

1 X 1 E= 50 W

(r)

E(θk ).

k=T −W +1

r=1

1.2 1.1 1.0 0.9 0.8 0.7 0.6 0.5 0.4

0.7

SGD Polyak Nesterov

101

102

Batch size

(a) s = 0.1, β = 2

103

0.6

Loss

Loss

Loss

For each method and batch size, Figure 1 (right) reports the minimum E over the corresponding search grid. Figure 6 shows additional results on the batch-size phase diagram across (s, β).

0.5 SGD Polyak Nesterov

0.4 101

102

Batch size

0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.02

103

(b) s = 0.1, β = 4

SGD Polyak Nesterov

101

102

Batch size

103

(c) s = 0.3, β = 6

Figure 6: Additional experiments on the batch-size phase diagram under different (s, β).

G.4

Data-scaling exponents across batch-size regimes

In Figure 3, we simulate PLK regression with (s, β) = (0.3, 4), σ = 0.1, and θ−1 = θ0 = 0, using data budgets D ∈ {2000, 4000, 8000, 16000, 32000, 64000, 128000}. For a prescribed batch-size exponent b, we set B = round(CB Db ) and T = ⌊D/B⌋, where CB is an order-one constant. 68

For each method and D, we tune the multiplicative constants around the learning-rate and momentum scalings prescribed by Theorem 6.3, and report the best stable configuration. Each run is evaluated by averaging the risk over the final 1% of optimization steps, and each point is averaged over 10 independent runs. We evaluate one representative batch-size exponent from each of the three regimes. • Small batch. We take bsmall = 0. For all three methods, we search around η ≂ D−s/(s+1) ; for Polyak and Nesterov, we additionally use 1 − ρ ≂ 1. 2sβ+1 • Large batch. We take blarge = b2 − δ, where b2 = 2sβ+β+1 and δ > 0 is small. SGD uses −(b η ≂ 1. For Polyak, we search around η ≂ 1 and 1 − ρ ≂ D large −s/(s+1)) . For Nesterov, we search around η ≂ 1 and

1 − ρ ≂ D−cnesterov ,

cnesterov = blarge −

sβ − (β − 1)blarge . sβ + 1

• Ultra-large batch. We choose b3 < bultra < 1. SGD uses η ≂ 1, while for both momentum methods we search around η ≂ 1 and 1 − ρ ≂ D−(1−bultra −δ) , with δ > 0 chosen so that T (1 − ρ) → ∞. Additional results. Small-batch

5 × 10−1

4 × 10−1

3 × 10−1

104

Data D

8 × 10−1

Ultra large-batch

6 × 10−1

Loss

6 × 10−1

Loss

Large-batch

SGD Polyak Nesterov Prediction

Loss

7 × 10−1

Figures 7 and 8 show additional data-scaling results across (s, β).

7 × 10−1

SGD Polyak Nesterov Prediction

105

SGD Polyak Nesterov Prediction

104

105

104

Data D

105

Data D

Figure 7: Data scaling with (s, β) = (0.1, 4).

Small-batch

10−1

4 × 10−2

Data D

105

Ultra large-batch

2.4 × 10−1

Loss

2.2 × 10−1

10−1 6 × 10−2

104

2.6 × 10−1

2 × 10−1

Loss

Loss

2 × 10−1

Large-batch

3 × 10−1

SGD Polyak Nesterov Prediction

2 × 10−1 1.8 × 10−1 1.6 × 10−1

SGD Polyak Nesterov Prediction

1.4 × 10−1 10

4

10

5

Data D

SGD Polyak Nesterov Prediction

104

105

Data D

Figure 8: Data scaling with (s, β) = (0.3, 2).

G.5

Delayed instability under increasing momentum

We simulate PLK regression with (s = 0.3, β = 4). We set the learning rate to η = 0.05 and the batch size to B = 4. To isolate the effect of multiplicative noise, we set the additive noise to σ = 0. For Polyak and Nesterov, we use the time-varying momentum schedule ρk = 1 −

1 . k+1

Each method is run for at most 5 × 105 iterations over 50 independent runs, and the figure reports the median risk.

69

Related documents

• Accelerating Compound LLM Training Workloads with Maestro#175196
• Pallas: A Proactive KV Cache Migration Framework for LLM Inference in AI-RAN#622793
• Pallas: A Proactive KV Cache Migration Framework for LLM Inference in AI-RAN#459680
• Hardware-Algorithm Co-Optimization of Early-Exit Neural Networks for Multi-Core Edge Accelerators#1156239
• Deep-submicron embedded processor architectures for high-performance, low-cost and low-power#456709