arXiv:2605.08871v1 [math.OC] 9 May 2026
Rennala MVR: Improved Time Complexity for Parallel Stochastic Optimization via Momentum-Based Variance Reduction Zhirayr Tovmasyan KAUST [email protected]
Artavazd Maranjyan KAUST [email protected]
Peter Richtárik KAUST [email protected]
Abstract Large-scale machine learning models are trained on clusters of machines that exhibit heterogeneous performance due to hardware variability, network delays, and system-level instabilities. In such environments, time complexity rather than iteration complexity becomes the relevant performance metric for optimization algorithms. Recent work by Tyurin and Richtárik [2023] established the first time complexity analysis for parallel first-order stochastic optimization, proposing Rennala SGD as a time-optimal method for smooth nonconvex optimization. However, Rennala SGD is fundamentally a modification of SGD, and variance reduction techniques are known to improve the iteration complexity of SGD. In this work, we investigate whether variance reduction can also improve time complexity in heterogeneous systems. We show that, under a mean-squared smoothness assumption, variance reduction can improve time complexity in relevant parameter regimes. To this end, we propose Rennala MVR, a variance-reduced extension of Rennala SGD based on momentum-based variance reduction, and analyze its oracle and time complexity. We establish lower bounds for time complexity under these assumptions. On a stochastic quadratic benchmark, experiments with the exact method support the theory, while neural-network experiments with a practical inexact variant show similar empirical gains over Rennala SGD.
1
Introduction
Modern machine learning models have grown to a scale that makes training on a single machine impractical. As a result, training is performed in distributed environments spanning multiple machines with heterogeneous computation speeds, communication delays, and system variability. In such settings, the efficiency of optimization algorithms is determined not only by the number of iterations required for convergence, but also by their time complexity, which is critically affected by system-level heterogeneity. Traditionally, optimization methods for training machine learning models have been analyzed through the lens of iteration complexity, which measures how many algorithmic steps are needed to reach a target accuracy. While this metric has been instrumental for understanding algorithmic efficiency in idealized settings, it becomes insufficient in distributed and heterogeneous environments. In practice, two methods with similar iteration complexity can exhibit vastly different runtimes due to idle time, synchronization overhead, and other system-level bottlenecks. This mismatch has motivated a Preprint.
Table 1: Comparison of time and oracle complexities of our method with the state-of-the-art parallel first-order method Rennala SGD, together with our derived lower bounds. We consider a system of n workers, where worker i requires τi seconds to compute a stochastic gradient; without loss of generality, we assume τ1 ≤ τ2 ≤ · · · ≤ τn . Both time and oracle complexities are reported for achieving an ε–stationary point in the nonconvex setting. The oracle complexity counts the total number of stochastic gradient evaluations used in the optimizer updates over the entire training process. The stochastic gradients computed by the workers are assumed to be unbiased with bounded variance 2 σ 2 (Assumption 2.3), i.e., Eξ∼D [∥∇f (x; ξ) − ∇f (x)∥ ] ≤ σ 2 for all x ∈ Rd . We denote by 0 ∗ ∗ ∆ := f (x )−f the initial suboptimality gap, where f = inf x f (x) (Assumption 2.2). Here, L is the smoothness constant of f , i.e., ∥∇f (x) − ∇f (y)∥2 ≤ L2 ∥x − y∥2 , and L̄ denotes the mean-squared smoothness constant (Assumption 2.4), defined by Eξ∼D ∥∇f (x; ξ) − ∇f (y; ξ)∥2 ≤ L̄2 ∥x − y∥2 . All stated complexities hide universal constant factors. Our method outperforms Rennala SGD in both oracle and time complexity for sufficiently small ε and when L̄ = O(L); see Section 5.2. Although our method does not always match the lower bound in terms of time complexity, it matches the lower bound in oracle complexity and, in certain regimes, can also achieve optimal time complexity. A detailed discussion is provided in Section 6.
Rennala SGD
L∆ min ε m∈[n]
[Tyurin and Richtárik, 2023] Rennala MVR (new)
(Theorem 5.2) Lower Bound (new) (Theorem 6.7)
Oracle Complexity (†)
Time Complexity
Algorithm
L̄∆ σ √ ε + ε
min
P m i=1
m −1 P 1 σ
m∈[n]
i=1
√
τi
L̄∆ √ + 1 min σ ε m∈[n]
ε
1 τi
−1 2 σ ε
+m
L∆σ 2 ε2
m −1 P 1 σ2 + m + min τi ε +m m∈[n]
P m i=1
1 τi
−1 2 σ ε
2 L̄∆σ + σε ε3/2
i=1
+m
(†) Oracle complexity counts the total number of stochastic gradient evaluations used in the optimizer updates. Define T (B) = min
2 L̄∆σ + σε (‡) ε3/2
m∈[n]
P
m 1 i=1 τi
−1
(B +
m) By Lemma D.1 2T (B) upper bounds the time required to obtain B stochastic gradients asynchronously. Therefore, if a method has time complexity KT (B), then its oracle complexity is obtained by multiplying K by the number of stochastic gradient evaluations used per iteration. In particular, for Rennala SGD this gives KB, while for Rennala MVR it gives B0 + 2KB, since each iteration uses B gradient pairs, i.e. 2B stochastic gradient evaluations, and the initialization costs B0 . (‡) This oracle complexity matches the known lower bound on oracle complexity established by Arjevani et al. [2022].
shift toward time complexity as a more faithful measure of algorithmic performance in real-world distributed systems. Recent work by Tyurin and Richtárik [2023] derived the first lower bounds on time complexity for parallel first-order methods with smooth nonconvex losses. They matched this lower bound with a method called Rennala SGD, which is a modification of classical stochastic gradient descent (SGD) [Robbins and Monro, 1951]. This line of work highlights a fundamental principle: achieving optimal performance in distributed systems often requires redesigning classical algorithms to account for system heterogeneity. Beyond SGD, variance reduction techniques have been shown to substantially improve optimization efficiency in the nonconvex regime under iteration complexity analysis. Methods such as STORM [Cutkosky and Orabona, 2019] (which uses momentum variance reduction (MVR)) or SNVRG [Zhou et al., 2020] reduce gradient noise and achieve stronger convergence guarantees than vanilla SGD under additional assumptions. However, despite their strong theoretical and empirical advantages, the time complexity of variance-reduced methods in heterogeneous distributed systems remains unexplored. This raises a natural question: if variance-reduced methods outperform SGD in terms of iteration complexity, can they also achieve better time complexity, and how should they be modified to do so? In this work, we initiate a systematic study of time complexity for distributed variance-reduced optimization. We propose Rennala MVR, a variance-reduced extension of Rennala SGD, and show that it achieves better time complexity under the stronger mean-squared smoothness assumption (Assumption 2.4). 1.1
Contributions
Our contributions are summarized as follows: 2
• We introduce Rennala MVR (Algorithm 1), a variance-reduced distributed optimization method for heterogeneous systems under the mean-squared smoothness assumption. • We derive upper bounds on its iteration and time complexity (Section 5), clarifying when variance reduction can yield time complexity gains over Rennala SGD. • We establish a new lower bound on achievable time complexity in this regime (Section 6). • We extend the time complexity analysis to the arbitrarily varying worker speeds setting (Section A). • We provide empirical results (Section 7) validating our theoretical findings and demonstrating practical performance improvements. 1.2
Related Work
Time complexity analysis for parallel stochastic optimization has recently emerged as a principled alternative to iteration complexity, especially in heterogeneous environments where the time per iteration depends on the workers’ speeds. Early modern treatments of time complexity for asynchronous methods—which avoid global synchronization at the cost of using stale gradients—include works by Mishchenko et al. [2022], Koloskova et al. [2022], Alahyane et al. [2025]. Asynchronous optimization itself dates back to classical work such as Tsitsiklis et al. [1986], and was later popularized in machine learning by Recht et al. [2011], Agarwal and Duchi [2011]; see Assran et al. [2020] for a survey. Building on these developments, Tyurin and Richtárik [2023] initiated a general time complexity framework for parallel stochastic optimization and showed that several existing asynchronous methods are not time-optimal under their model. This led to the Rennala SGD method, which attains an optimal time complexity. Subsequent works refined and extended these ideas and produced time-optimal variants of asynchronous SGD, including [Maranjyan et al., 2025c, Maranjyan and Richtárik, 2026]. Further extensions and related developments include [Maranjyan et al., 2025b, Wu and Luo, 2026], and a comprehensive overview appears in the work of Maranjyan [2025]. Among these, Maranjyan et al. [2025c] is particularly relevant to our work: their Ringmaster ASGD method is an asynchronous SGD with the same optimal time complexity as Rennala SGD. Beyond the fixed-computation-time setting of Assumption 2.1, several papers extend the framework of Tyurin and Richtárik [2023] by considering alternative time models and settings; see, e.g., [Tyurin et al., 2024a,b, Tyurin and Richtárik, 2024, Maranjyan et al., 2025a]. The second line of work most relevant to this paper is variance reduction. Classical variance-reduced methods—including SVRG [Johnson and Zhang, 2013, Allen-Zhu and Yuan, 2016], S2GD [Konečný and Richtárik, 2015], SAGA [Defazio et al., 2014], SPIDER [Fang et al., 2018], SARAH [Nguyen et al., 2017], JacSketch [Gower et al., 2018], L-SVRG [Kovalev et al., 2019], SNVRG [Zhou et al., 2020], PAGE [Li et al., 2021], and STORM [Cutkosky and Orabona, 2019]—provide improved oracle complexity guarantees in stochastic nonconvex optimization. Despite their strong theory, they have been less common in large-scale deep learning practice; one explanation is discussed by Defazio and Bottou [2019]. More recently, Yuan et al. [2025] revisited variance reduction for training large language models and proposed a practical modification inspired by STORM [Cutkosky and Orabona, 2019], reporting speedups in LLM training.
2
Problem Setup
We consider the nonconvex optimization problem minimize {f (x) := Eξ∼D [f (x; ξ)]} , x∈Rd
(1)
where f (x; ξ) is the loss function evaluated on a data sample ξ drawn from distribution D, and the model is parameterized by x ∈ Rd with d denoting the parameter dimensionality. We consider a distributed learning setting with n workers, where each worker i has access to the same data distribution D. This setting is common in data centers with either unified memory or uniformly partitioned data. Following the fixed computation model [Mishchenko et al., 2022], we formalize the heterogeneous computation times as follows: 3
Assumption 2.1 (Fixed Computation Model). Each worker i requires τi seconds to compute one stochastic gradient ∇f (x; ξ). Without loss of generality, we assume 0 < τ1 ≤ τ2 ≤ · · · ≤ τn . We assume instantaneous communication (zero latency) between workers and the server in both directions. This is the standard modeling assumption in prior work on time complexity for distributed stochastic optimization [Mishchenko et al., 2022, Koloskova et al., 2022, Tyurin and Richtárik, 2023, Maranjyan et al., 2025c], and we adopt it here in order to make a direct comparison with these results. We stress that this is a modeling simplification rather than a claim that communication is negligible in practice. Explicitly modeling communication costs leads to a substantially richer theoretical problem, since one must then also specify which communication-reduction mechanisms are allowed, such as compression, quantization, sparsification, or local updates. These questions are important, but largely orthogonal to the present paper, whose goal is to understand how variance reduction affects time complexity under heterogeneous worker speeds within the standard theoretical model. 2.1
Assumptions
We make the following standard assumptions: Assumption 2.2 (Lower boundedness). There exists f ∗ > −∞ such that f (x) ≥ f ∗ for all x ∈ Rd . We define ∆ := f (x0 ) − f ∗ , where x0 is the starting point of the optimization methods. Assumption 2.3. For every ξ, the function f (x; ξ) is differentiable with respect to its first argument x. Moreover, the stochastic gradients are unbiased and have bounded variance σ 2 ≥ 0, that is, Eξ∼D [∇f (x; ξ)] = ∇f (x), ∀x ∈ Rd , Eξ∼D ∥∇f (x; ξ) − ∇f (x)∥2 ≤ σ 2 , ∀x ∈ Rd . Assumption 2.4 (Mean-squared smoothness). There exists L̄ > 0 such that for all x, y ∈ Rd , Eξ∼D ∥∇f (x; ξ) − ∇f (y; ξ)∥2 ≤ L̄2 ∥x − y∥2 . This assumption is stronger than classical L–smoothness of f , as shown in the following lemma. Lemma 2.5 (Proof in Section B.1). Mean-squared smoothness (Assumption 2.4) implies that f is L̄–smooth, i.e. ∥∇f (x) − ∇f (y)∥ ≤ L̄ ∥x − y∥ , ∀x, y ∈ Rd . Mean-squared smoothness also implies expected similarity, a centered gradient-difference bound closely related to the notion of Hessian variance introduced by Szlendak et al. [2021]. Lemma 2.6 (Proof in Section B.2). Mean-squared smoothness (Assumption 2.4) implies expected similarity with the same constant L̄: h i 2 2 E ∥∇f (x; ξ) − ∇f (y; ξ) − (∇f (x) − ∇f (y))∥ ≤ L̄2 ∥x − y∥ , ∀ x, y ∈ Rd . Under these assumptions, our objective is to find an ε–stationary point: a (possibly random) point x satisfying E[∥∇f (x)∥2 ] ≤ ε.
3
Background and Motivation
Under Assumptions 2.2 and 2.3, combined with L–smoothness of f (instead of Assumption 2.4), the SGD method with iterations xk+1 = xk − γ∇f (xk ; ξ k ) and i.i.d. ξ k ∼ D achieves optimal oracle 2 complexity [Ghadimi and Lan, 2013, Arjevani et al., 2022] of O(L∆/ε + σ L∆/ε2 ). Oracle complexity measures the total number of stochastic gradient evaluations; on a single machine, this corresponds to the iteration complexity. However, achieving time-optimal performance with n parallel machines requires effectively parallelizing SGD. Tyurin and Richtárik [2023] addressed this question by proposing the Rennala SGD method, which is a minibatch extension of SGD. Instead of computing a single gradient, Rennala SGD collects B gradients and takes a step using their average. The key to achieving optimal time complexity is 4
distributing the batch of size B across all available machines and computing it asynchronously (as 2 quickly as possible). By setting B = max{1, σ /ε}, Rennala SGD attains the lower bound on time complexity. Importantly, the oracle complexity—the total number of gradient evaluations—remains 2 O(σ L∆/ε2 ), matching the lower bound up to a constant factor (see the last column in Table 1). When we additionally invoke Assumption 2.4, which is stronger than L–smoothness of f (Lemma 2.5), we can achieve better oracle complexity on a single machine using variance reduction methods such as STORM [Cutkosky and Orabona, 2019] or SNVRG [Zhou et al., 2020]. These methods improve 2 the oracle complexity to O(L̄∆σ/ε3/2 + σ /ε) and are optimal [Arjevani et al., 2022]. This motivates our central question: does time complexity also improve under the stronger Assumption 2.4, and which method can we use to surpass Rennala SGD’s time complexity? We focus on the STORM algorithm [Cutkosky and Orabona, 2019] and try to make use of its MVR technique because it has favorable properties and is easier to parallelize across multiple machines.
4
Rennala MVR
We first recall the MVR technique of Cutkosky and Orabona [2019], the main technique behind their STORM algorithm. Given a starting point x0 ∈ Rd , an initial gradient estimator g 0 , a momentum parameter p ∈ (0, 1], and a stepsize γ > 0, The MVR update at each iteration k takes the following form: xk+1 = xk − γg k , g
k+1
= ∇f (x
k+1
(2) k
k
k
k
; ξ ) + (1 − p)(g − ∇f (x ; ξ )),
(3)
k
where ξ is sampled i.i.d. from D. Note that setting p = 1 recovers standard SGD. A distinctive feature of MVR is that it uses two stochastic gradients per iteration, evaluated at xk and at xk+1 using the same sample ξ k . To extend this idea to the distributed parallel setting, we follow the same paradigm as in Rennala SGD: instead of forming a single stochastic gradient, the server continuously aggregates stochastic gradients computed by workers and constructs minibatch estimators. In our case, each worker computes two gradients—one at xk and one at xk+1 —for each sampled data point, and the server forms two minibatch averages from the first B arrivals. The resulting method is summarized in Algorithm 1. Algorithm 1 Rennala MVR 1: Input: initial point x0 ∈ Rd (stored on both the server and the workers), stepsize γ > 0,
minibatch size B ∈ {1, 2, . . .}, momentum parameter p ∈ (0, 1], initial batch size B0
2: Compute initial gradient estimator g 0 asynchronously using all workers with batch size B0 3: for k = 0, . . . , K − 1 do 4: Update the model: xk+1 := xk − γg k Broadcast xk+1 to workers (workers keep xk cached), then workers compute gradients at both 5:
xk and xk+1 6: Initialize g − := 0, g + := 0 and b := 0 7: while b < B do 8: Gradients ∇f xk ; ξik,b and ∇f xk+1 ; ξik,b arrive from worker ik,b k,b k,b k,b 9: g − ← g − + ∇f xk ; ξik,b g + ← g + + ∇f xk+1 ; ξik,b k,b 10: Worker ik,b immediately begins computing new gradients at xk and xk+1 11: b←b+1 12: end while+ − 13: g k+1 := gB + (1 − p)(g k − gB ) 14: end for
4.1
Algorithm Description
The algorithm begins by obtaining an initial gradient estimator g 0 (for example, a minibatch average computed at x0 ). At iteration k, the server first updates the model via xk+1 = xk − γg k , and then 5
ensures that the workers have access to the pair (xk , xk+1 ). In practice, it suffices to broadcast only the new point xk+1 since xk was broadcast in the previous iteration and can be cached by the workers. Workers operate asynchronously and continuously. Each time a worker samples ξ, it computes the gradient pair ∇f (xk ; ξ), ∇f (xk+1 ; ξ) and sends this pair to the server. Thus, one arrival to the server consists of two stochastic gradients corresponding to the same sample ξ. The server collects the first B such arrivals, that is, B gradient pairs, and forms the sums g − and g + from the gradients evaluated at xk and xk+1 , respectively. Equivalently, each iteration uses B gradient pairs, or 2B stochastic gradient evaluations. Once these B arrivals have been collected, the server constructs the next estimator as g− g+ k k+1 + (1 − p) g − , g = B B and then proceeds to the next iteration. 4.2
Connection to Rennala SGD +
When p = 1, the update simplifies to g k+1 = g /B , so the method reduces to minibatch SGD as a special case, using only gradients at the current point xk+1 . In this regime, the additional computation of gradients at xk (i.e., the g − term) is unnecessary. If we therefore modify the worker routine to compute only ∇f (xk+1 ; ξ), the resulting asynchronous minibatch SGD implementation coincides with Rennala SGD [Tyurin and Richtárik, 2023].
5
Theoretical Results
This section presents the main theoretical guarantees for Rennala MVR (Algorithm 1). We first establish its iteration complexity and the resulting oracle complexity—the total number of stochastic gradient evaluations computed by the algorithm. We then derive a time complexity bound, which is the more informative metric in parallel settings with heterogeneous workers. 5.1
Iteration Complexity
We begin with an iteration complexity result under our standard assumptions. Theorem 5.1 (Iteration Complexity; Proof in Section D.1). Under Assumptions 2.2–2.4, let the stepsize in Rennala MVR (Algorithm 1) beγ = 1/4L̄. Fix ε > 0 and assume ε < σ 2 and ε < 2L̄∆. √ 2 √ Choose B = ⌈6σ/ ε⌉, p = ε/σ and B0 = 6σ /ε . Then, K−1
i 1 X h 2 E ∇f (xk ) ≤ε, K k=0
for K≥
24∆L̄ σ +√ . ε ε
Our method can be viewed as a minibatch variant of the momentum variance-reduction (MVR) mechanism used in STORM [Cutkosky and Orabona, 2019]. The analysis is correspondingly simpler here, since we focus on the MVR component and do not require the additional algorithmic features of STORM. A direct corollary of Theorem 5.1 is an oracle complexity bound. The algorithm runs for K iterations; each iteration computes 2B stochastic gradients, and the initialization costs B0 stochastic gradients. Hence, the total oracle complexity (i.e., the number of stochastic gradient evaluations) is L̄∆σ σ 2 + . (4) B0 + 2B · K = O ε ε3/2 This matches the lower bound for this problem class established by Arjevani et al. [2022], and therefore Rennala MVR is oracle-optimal (up to constant factors) under Assumptions 2.2–2.4. 6
For comparison, in the classical L–smooth setting (replacing mean-squared smoothness Assumption 2.4), Rennala SGD achieves oracle complexity L∆σ 2 , O ε2 which is also optimal in that setting [Arjevani et al., 2022]; see Table 1. Nevertheless, this dependence on ε is worse than the oracle complexity of Rennala MVR in (4). 5.2
Time Complexity
We now turn to time complexity. In parallel optimization with heterogeneous workers, iteration complexity alone can be misleading, since the time per iteration depends on how quickly workers return gradients. Theorem 5.2 (Time Complexity; Proof in Section D.2). Under the assumptions and parameter choices of Theorem 5.1, the time complexity of Rennala MVR (Algorithm 1) is given by !−1 !−1 m m 2 X X L̄∆ σ σ σ 1 1 √ + m + min T = O +√ min +m . ε τ ε ε m∈[n] i=1 τi ε m∈[n] i=1 i Define T (B) = min m∈[n]
m X 1 i=1
τi
!−1 (B + m) ,
which upper bounds the time required to obtain B stochastic gradients asynchronously (Lemma D.1). Using this notation, the bound in Theorem 5.2 can be written more compactly as 2 L̄∆ σ σ σ TMVR = O +√ ·T √ +T . ε ε ε ε Comparison with Rennala SGD.
For Rennala SGD, the corresponding time complexity is 2 L∆ σ TSGD = O ·T . ε ε
The key distinction is the per-iteration waiting time. Each iteration of Rennala MVR waits for roughly σ/√ε gradient pairs and therefore costs on the order of T (σ/√ε) time, whereas each iteration of 2 Rennala SGD costs T (σ /ε), which is always at least as large. This smaller per-iteration time is the main source of potential speedups for Rennala MVR. On the other hand, the dominant iteration count for Rennala MVR (for sufficiently small ε) scales as L̄∆/ε, which can be larger than the L∆/ε scaling of Rennala SGD since L ≤ L̄. Consequently, when L̄ = O(L), the reduction in per-iteration time dominates and Rennala MVR enjoys an improved overall time complexity compared to Rennala SGD. Comparison with the lower bound. We next relate Theorem 5.2 to the time lower bound proved in the following section. The lower bound takes the form (Theorem 6.7) 2 L̄∆ σ √ +1 ·T T =Ω . ε σ ε First, note that the implied oracle complexity (recall: the number of stochastic gradient evaluations) is 2 L̄∆ σ L̄∆σ σ 2 √ +1 · = 3/2 + , ε ε σ ε ε which matches the oracle complexity of Rennala MVR from Theorem 5.1 (and Table 1). In other words, Rennala MVR computes the minimal number of stochastic gradients up to constant factors. The remaining gap is a batching/synchronization issue. The lower bound corresponds to producing 2 √ gradients in batches of size σ /ε, while Rennala MVR uses smaller batches of size σ/ ε. Smaller 7
batches may require more frequent synchronization (i.e., more frequent termination of ongoing gradient computations), which can increase the elapsed time even when the total number of computed gradients is optimal. Despite this gap in general, the upper bound can match the lower bound in regimes where T
σ2 ε
σ ≈ √ ·T ε
σ √ ε
. √
For example, this occurs in the homogeneous worker case with n = σ/ ε and τi = τ for all i ∈ [n]. Overall, the fact that the current time upper bound for Rennala MVR does not match the lower bound highlights a key distinction between variance-reduced methods and plain SGD under heterogeneous compute-time models: optimal oracle complexity does not automatically imply optimal time complexity. Closing this gap appears to require a method that can operate with much larger effective 2 batch sizes (on the order of σ /ε) while preserving the desired convergence rate. Within the current proof approach and stepsize restrictions, it is unclear how to achieve this via a simple modification of Rennala MVR, which suggests that fundamentally different algorithmic ideas may be needed. This issue is discussed in more detail in Section C.
6
Lower Bound on Time Complexity
In this section we formalize a framework for proving lower bounds on the time required to find an ε–stationary point in stochastic nonconvex optimization with heterogeneous parallel workers. Following Tyurin and Richtárik [2023], we specify function class F, an oracle class O, and an algorithm class A together with a time-based interaction protocol. We then define a minimax notion of time complexity and state our lower bound. Compared to the protocol in Tyurin and Richtárik [2023], our setting requires one additional feature: the algorithm may request variable batch sizes up to a maximum of B points per interaction (in particular, our Rennala MVR (Algorithm 1) asks for two stochastic gradients). 6.1
Protocol
We consider n oracles (workers) running in parallel. At each interaction, the algorithm returns (i) which oracle to query, (ii) the time at which it requests the reply, and (iii) a batch of up to B query points. The protocol is given below in Algorithm 2. In Protocol 2, ski denotes the internal Protocol 2 Time Multiple Oracles Protocol (variable batch size up to B) function f ∈ F, algorithm A ∈ ((O1 , . . . , On ), (D1 , . . . , Dn )) ∈ O(f ) 2: s0i = (0, 0, 0, 0) for all i ∈ [n] 3: t0 = 0 4: for k = 0, 1, 2, . . . do 5: (tk+1 , ik+1 , bk , X k ) = Ak (G1 , . . . , Gk ) , Gk+1 ) = Oik+1 (tk+1 , bk , X k , skik+1 , ξ k+1 ), 6: (sk+1 ik+1 7: sk+1 = skj ∀j ̸= ik+1 j 8: end for
1: Input:
A,
oracles and distributions
ξ k+1 ∼ Dik+1
state of oracle i at interaction k (defined precisely below). At interaction k, the algorithm outputs an oracle index ik+1 ∈ [n], a time tk+1 ≥ tk , a batch size bk ∈ [B], and a tuple of query points X k = (xk,1 , . . . , xk,B ) ∈ (Rd )B , with the convention xk,j = 0 for j > bk . The oracle returns a tuple Gk+1 = (g k+1,1 , . . . , g k+1,B ) ∈ (Rd )B in which only the first bk entries may be nonzero (the remaining entries are again 0 by convention). Setting B = 1 recovers the protocol of Tyurin and Richtárik [2023]. 8
6.2
Oracle Model
Each oracle Oi has an internal state encoding whether it is idle or busy and, if busy, which batch it is currently processing. We represent the state as si = (st , sX , sb , sq ) ∈ R≥0 × (Rd )B × {0, . . . , B} × {0, 1}, where sq = 0 indicates that the oracle is idle, and sq = 1 indicates that it is busy computing a batch of size sb started at time st for the stored query tuple sX (with the convention that (sX )j = 0 for j > sb ). Formally, define Oτfi (·),B : R≥0 × {1, . . . , B} × (Rd )B × (R≥0 × (Rd )B × {0, 1, . . . , B} × {0, 1}) × Di |{z} |{z} | {z } | {z } | {z } time
batch size
points
input state
randomness
→ (R≥0 × (Rd )B × {0, 1, . . . , B} × {0, 1}) × (Rd )B . {z } | {z } | output state
gradients
It is defined by if sq = 0, ((t, X, b, 1), 0), Oτfi (·),B (t, b, X, (st , sX , sb , sq ), ξ) = ((st , sX , sb , 1), 0), if sq = 1 and t < st + τi (sb ), (5) ((0, 0, 0, 0), G), if sq = 1 and t ≥ st + τi (sb ), where 0 denotes the all-zero element of (Rd )B , and G ∈ (Rd )B is given by ∇f ((sX )j ; ξ), j ≤ sb , Gj = 0, j > sb . Thus, if the oracle is idle, it starts computing the gradients for the submitted batch (of size b) and becomes busy. While busy, it returns the all-zero tuple until the computation finishes. Once t reaches the completion time st + τi (sb ), the oracle returns the stored batch gradients and resets to idle. In particular, each completed reply uses a single random sample ξ ∼ Di shared across the sb gradients in that reply. Relative to Tyurin and Richtárik [2023], our oracle model differs in two ways. First, we allow B > 1 (variable batch sizes). Second, the runtime is governed by a batch-time function τi (·) rather than a constant: τi (sb ) is the time required for worker i to compute sb stochastic gradients. We impose the following natural monotonicity property. Assumption 6.1 (Batch time). For each oracle i ∈ [n], the batch-time function τi (·) is nondecreasing: for any 1 ≤ k ≤ ℓ, τi (k) ≤ τi (ℓ). 6.3
Algorithm Class
At each interaction, the algorithm can use all previously received replies {G1 , . . . , Gk } to select the next oracle, the request time, the batch size, and the query points. We also impose the natural constraint tk+1 ≥ tk (the algorithm cannot query in the past). Definition 6.2 (Algorithm class). An algorithm A = {Ak }∞ k=0 is a sequence of mappings such that for each k ≥ 0, k Ak : (Rd )B → R≥0 × [n] × {1, . . . , B} × (Rd )B , and if (tk+1 , ·) = Ak (G1 , . . . , Gk ) and (tk , ·) = Ak−1 (G1 , . . . , Gk−1 ), then tk+1 ≥ tk . Moreover, if Ak (·) = (tk+1 , ik+1 , bk , X k ), then X k = (xk,1 , . . . , xk,B ) satisfies xk,j = 0 for all j > bk . We denote the class of all such algorithms by A. As in Tyurin and Richtárik [2023], we restrict attention to zero-respecting algorithms. Definition 6.3 (Zero-respecting). Consider an execution of Protocol 2. For each interaction r ≥ 0, let (tr+1 , ir+1 , br , X r ) 9
be the algorithm output, where X r = (x(r,1) , . . . , x(r,B) ) ∈ (Rd )B . Let Grcomp = (g (r,1) , . . . , g (r,B) ) ∈ (Rd )B denote the gradient tuple eventually returned when the batch submitted at interaction r is completed, with the convention g (r,k) = 0 for all k > br . The algorithm is zero-respecting if for all r ≥ 0 and all k ∈ [B], [ [ ′ support x(r,k) ⊆ support g (s,k ) . s<r k′ ∈[B]
We denote the class of all zero-respecting algorithms by Azr . 6.4
Function and Oracle Classes
We define the following. Definition 6.4 (Function class F∆,L̄ ). We define F∆,L̄ as the set of all L̄–smooth functions f : Rd → R such that f (0) − inf x f (x) ≤ ∆. 2
,L̄,B Definition 6.5 (Oracle class). For any f ∈ F∆,L̄ , the class Oτσ1 (·),...,τ (f ) returns a collection n (·)
of oracles Oi = Oτfi (·),B and distributions Di for all i ∈ [n], where stochastic gradients ∇f (·; ξ) are unbiased and σ 2 -variance-bounded (Assumption 2.3) and satisfy mean-squared smoothness (Assumption 2.4) with constant L̄. The oracle Oτfi (·),B is given by (5). 6.5
Time Complexity Measure
Let P[F∆,L̄ ] denote the set of all probability distributions over F∆,L̄ . We define the minimax time complexity as follows. Definition 6.6. For a given batch-budget B and problem parameters (∆, L̄, σ 2 ), we define the (minimax) time complexity as the smallest time t by which a zero-respecting algorithm can ensure an ε–stationary point in expectation, uniformly over all admissible oracle models and hard instance distributions: mtime (B, ∆, L̄, σ 2 ) :=
sup
sup
2
,L̄,B O∈Oτσ (·),...,τ 1
n (·)
inf inf
PF ∈P[F∆,L̄ ] A∈Azr
n
t≥0 E
inf
(r,k)
∇f (x
2
)
o ≤ε ,
(r,k)∈St
where, for each interaction r ≥ 0, we denote by cr the completion time of the batch submitted at interaction r, and define St := { (r, k) ∈ N0 × [B] | 1 ≤ k ≤ br , cr ≤ t }. We can now state our lower bound on the time complexity under the above protocol and classes. Theorem 6.7 (Proof in Section E.1). Fix ∆ > 0, L̄ > 0, σ 2 > 0, 0 < ε < c′ L̄∆, an integer B ≥ 1, and n workers with batch-time functions {τi (·)}ni=1 that satisfy Assumption 6.1. Write τi := τi (1) and assume 0 < τ1 ≤ · · · ≤ τn . Let the chain constants ∆0 , ℓ1 , γ∞ be as in Lemma E.3 and the estimator constants ς, ℓ̄1 as in Lemma E.4. Define j L∆ k 2ες 2 ℓ1 √ ℓ1 √ := := := p := min , 1 , L L̄ p , λ 2ε , T . σ2 L 2 ∆0 ℓ1 ε ℓ̄1 2
,L̄,B Then there exist f ∈ F∆,L̄ and an oracle class O ∈ Oτσ1 (·),...,τ (f ) such that, under Protocol 2, n (·)
2
mtime (B, ∆, L̄, σ ) ≥ c ·
√
L̄∆ min { ε/σ, 1} +1 ε 10
min m∈[n]
m X 1 i=1
τi
!−1
σ2 +m . ε
Rennala SGD avg(10 seeds) (\# 1): ° = 0:0625, B = 100
Rennala SGD avg(10 seeds) (\# 1): ° = 0:25, B = 80
Rennala SGD avg(10 seeds) (\# 2): ° = 0:03125, B = 100
101
10
Rennala MVR avg(10 seeds) (\# 3): ° = 0:5, S = 200, p = 0:001, S0 = 40000, ® = 1:0
0
0.2
0.4
0.6
0.8
Runtime (seconds)
(a) τi =
1.0 1e6
Rennala MVR avg(10 seeds) (\# 3): ° = 0:125, S = 20, p = 0:001, S0 = 400, ® = 1:0
krf(x t )k 2
krf(x t )k 2
krf(x t )k 2
0.0
Rennala MVR avg(10 seeds) (\# 2): ° = 0:5, S = 60, p = 0:005, S0 = 3600, ® = 1:0
10−3 10−4
10−5
10−6
Rennala SGD avg(10 seeds) (\# 3): ° = 0:25, B = 80
10−2
10−4
10−5
0
10−1
10−3
−4
10
Rennala MVR avg(10 seeds) (\# 3): ° = 0:5, S = 60, p = 0:001, S0 = 3600, ® = 1:0
10−2
10−3
Rennala SGD avg(10 seeds) (\# 2): ° = 0:5, B = 200 Rennala MVR avg(10 seeds) (\# 1): ° = 0:5, S = 60, p = 0:001, S0 = 3600, ® = 1:0
Rennala MVR avg(10 seeds) (\# 2): ° = 0:25, S = 40, p = 0:001, S0 = 1600, ® = 1:0
10−1
10−2
10
101
Rennala SGD avg(10 seeds) (\# 3): ° = 0:5, B = 200 Rennala MVR avg(10 seeds) (\# 1): ° = 0:5, S = 40, p = 0:001, S0 = 1600, ® = 1:0
Rennala MVR avg(10 seeds) (\# 2): ° = 0:125, S = 60, p = 0:001, S0 = 3600, ® = 1:0
10−1
Rennala SGD avg(10 seeds) (\# 1): ° = 0:125, B = 40
Rennala SGD avg(10 seeds) (\# 2): ° = 0:25, B = 100
101
Rennala SGD avg(10 seeds) (\# 3): ° = 0:03125, B = 80 Rennala MVR avg(10 seeds) (\# 1): ° = 1, S = 200, p = 0:001, S0 = 40000, ® = 1:0
100
10−5
0.0
0.2
0.4
0.6
0.8
Runtime (seconds)
√ i
1.0 1e6
0.0
0.2
0.4
0.6
0.8
Runtime (seconds)
(b) Mixture delays
1.0 1e6
(c) Uniform delays
Figure 1: Comparison of Rennala MVR and Rennala SGD on the stochastic quadratic benchmark with 10 workers under three delay models. Rennala SGD (\# 1): ° = 0:0625, B = 1
Rennala SGD (\# 2): ° = 0:25, B = 10
Rennala SGD (\# 2): ° = 0:25, B = 5
Rennala SGD (\# 3): ° = 0:125, B = 5
Rennala SGD (\# 3): ° = 0:125, B = 5
Rennala MVR (\# 1): ° = 1:0, S = 20, p = 0:01, S0 = 400, ® = 0:025
Rennala MVR (\# 1): ° = 1:0, S = 20, p = 0:1, S0 = 400, ® = 0:025
Rennala MVR (\# 3): ° = 1:0, S = 40, p = 0:005, S0 = 40, ® = 0:025
10−1 10−2 10
−3
100
Rennala SGD (\# 2): ° = 0:125, B = 5 Rennala SGD (\# 3): ° = 0:0625, B = 1
Rennala MVR (\# 2): ° = 1:0, S = 20, p = 0:05, S0 = 20, ® = 0:025 Rennala MVR (\# 3): ° = 0:5, S = 10, p = 0:2, S0 = 10, ® = 0:001
10−1
Rennala MVR (\# 1): ° = 1:0, S = 20, p = 0:05, S0 = 20, ® = 0:025
Training loss f(x t )
Rennala MVR (\# 2): ° = 1:0, S = 20, p = 0:01, S0 = 400, ® = 0:01
100
Rennala SGD (\# 1): ° = 0:25, B = 10
Rennala SGD (\# 1): ° = 0:25, B = 5
Training loss f(x t )
Training loss f(x t )
101
Rennala MVR (\# 2): ° = 1:0, S = 20, p = 0:1, S0 = 400, ® = 0:025 Rennala MVR (\# 3): ° = 1:0, S = 20, p = 0:05, S0 = 400, ® = 0:025
10
0
10−2
10−4 0
2000
4000
6000
8000
10000
0
2000
Runtime (seconds)
(a) τi =
√ i
4000
6000
8000
10000
0
2000
(b) Mixture delays
4000
6000
8000
10000
Runtime (seconds)
Runtime (seconds)
(c) Uniform delays
Figure 2: Comparison of the inexact variant of Rennala MVR and Rennala SGD on asynchronous neural-network training over a MNIST subset with 10 workers under three delay models. This figure reports experiments with the practical α-parameterized inexact variant from Section 7.2, not Algorithm 1 exactly.
7
Experiments
We evaluate our method against Rennala SGD [Tyurin and Richtárik, 2023] on two benchmarks: (i) a stochastic tridiagonal quadratic problem, where we test the exact method analyzed in the paper, and (ii) asynchronous neural-network training, where we use an inexact practical variant inspired by MARS [Yuan et al., 2025]. The first benchmark is theory-aligned, while the second is intended as an empirical test outside the scope of our current analysis. 7.1
Quadratic Benchmark
We consider the quadratic objective f (x) =
1 ⊤ x Ax − b⊤ x, 2
with dimension d = 100, where A ∈ Rd×d and b ∈ Rd are given by 2 −1 0 −1 . . .. .. 1 1 0 −1 A= b= , . . .. .. 4 4 .. . . −1 0 0 −1 2 √ We initialize at x0 = ( d, 0, . . . , 0). Each worker queries the same unbiased stochastic gradient oracle ∇f (x; ζ) = ∇f (x) + ζ, 2 where ζ ∼ N (0, σadd Id ) with σadd = 0.1. That is, each stochastic gradient is formed by adding isotropic Gaussian noise to the exact gradient. Distributed setup and tuning. For the quadratic experiments shown in Figures 1a to 1c, we model heterogeneity through fixed per-worker delays. We use n = 10 workers and a simulated time budget 11
of 106 . For each hyperparameter setting, we average over 10 independent runs with different √ per-node random seeds. We consider three delay models: the canonical square-root profile τi = i, randomly permuted across workers; delays sampled independently from the uniform distribution on [1, 10n]; and a three-peak Gaussian mixture, obtained by sampling three peak locations in [1, 10n], assigning each worker to one of the peaks, drawing its delay from the corresponding Gaussian component, and clipping the result to [1, 10n]. For both Rennala SGD and Rennala MVR, we tune the stepsize over γ ∈ {2j : j = −15, . . . , 2}, and the minibatch size over B ∈ {1, 5, 10, 20, 40, 60, 80, 100, 200}. For Rennala MVR, we additionally tune the momentum parameter over p ∈ {0.001, 0.005, 0.01, 0.05, 0.1, 0.2, 0.5, 0.9}, and the initialization batch size over B0 ∈ {B, B 2 }. We report the stationarity measure ∥∇f (x)∥2 as a function of time. In all plots, we show the three best-performing hyperparameter configurations for each method, where performance is evaluated by the median value of the reported metric over the final 1% of the time horizon. Results. The quadratic experiments in Figures 1a to 1c evaluate the exact method analyzed in this paper. Across all three delay models, the best-tuned configurations of Rennala MVR reach lower values of ∥∇f (x)∥2 than Rennala SGD in the 10-worker regime. This provides empirical support for the time complexity advantage predicted by our theory. 7.2
Neural Network Experiments
We also evaluate the practical behavior of Rennala MVR on asynchronous neural-network training. We train a two-layer ReLU network with hidden dimension 200 on a 3000-example MNIST subset using 10 heterogeneous workers, local batch size 4, and a simulated time budget of 104 . We again consider the square-root, uniform, and three-peak mixture delay models. Unlike the quadratic benchmark in Section 7.1, this experiment does not evaluate Algorithm 1 exactly. Instead, in the neural-network setting we study an inexact variant of Rennala MVR inspired by MARS [Yuan et al., 2025]. The motivation is practical: while the present network is small, for larger models the exact variant introduces additional runtime overhead, whereas the inexact variant reuses previous stochastic-gradient information and is therefore computationally more attractive. This choice is also consistent with the empirical findings of MARS, where the inexact variant was observed to perform comparably to the exact one in practice. Accordingly, the neural-network experiments should be interpreted as empirical evidence about the practical behavior of variance reduction in this setting, rather than as a direct validation of the theory developed for the exact method. We additionally introduce a scaling parameter α ∈ (0, 1] and consider the update xk+1 = xk − γg k , g k+1 = (1 − p)g k + p∇f (xk+1 ; ξ k ) + α(1 − p) ∇f (xk+1 ; ξ k ) − ∇f (xk ; ξ k−1 ) . The role of α is to scale the correction term. Recall, the update of the gradient estimator in our algorithm is given by (2) g k+1 = (1 − p)g k + p∇f (xk+1 ; ξ k ) + (1 − p) ∇f (xk+1 ; ξ k ) − ∇f (xk ; ξ k ) , where α = 1 and the last stochastic gradient is computed with the same current data point ξ k . This modification is motivated by the fact that pure momentum-based variance reduction often behaves poorly in neural-network training, whereas introducing a small correction coefficient can substantially improve stability and empirical performance, as also observed in MARS and related follow-up work. We tune the stepsize γ, the batch size B, the momentum parameter p, and the interpolation parameter α over the grids γ ∈ {2j : j = −10, . . . , 0}, B ∈ {1, 5, 10, 20, 40, 60}, p ∈ {0.001, 0.005, 0.01, 0.05, 0.1, 0.2}, α ∈ {0.001, 0.005, 0.01, 0.025}, 12
and set the initialization batch size to B0 ∈ {B, B 2 }. Here, B denotes the batch size used to construct the stochastic gradient estimator, and B0 denotes the batch size used for its initialization. For these neural-network runs, we report the training loss as a function of time. For visualization, the main figures show the best-performing hyperparameter configurations from the search grid, while Appendix F provides additional sensitivity plots. Across the tested parameter ranges, performance is relatively insensitive to the auxiliary parameters p and α, which suggests that, in practice, the dominant tuning burden remains the stepsize and batch size, as in Rennala SGD. Results. The neural-network experiments in Figures 2a to 2c evaluate the α-parameterized inexact variant of Rennala MVR introduced in this subsection rather than Algorithm 1 exactly. Although this setting lies outside the scope of our present analysis, it exhibits the same qualitative trend across all three delay models: the inexact variant achieves lower training loss than Rennala SGD under the square-root, mixture, and uniform delay models. These results suggest that the practical benefits of variance reduction may extend beyond the exact theory-covered setting.
8
Conclusion
This work initiates time complexity analysis for variance-reduced methods in heterogeneous parallel stochastic optimization. Focusing on the MVR mechanism in STORM [Cutkosky and Orabona, 2019], and working under the stronger mean-squared smoothness assumption, we proposed Rennala MVR and established its iteration, oracle, and time complexity guarantees. In particular, we showed that variance reduction can improve not only oracle complexity, but also time complexity in heterogeneous parallel environments, yielding improved guarantees over Rennala SGD [Tyurin and Richtárik, 2023], the previously time-optimal method in the standard smoothness regime, in relevant parameter regimes. More broadly, our results highlight a key conceptual message: in parallel stochastic optimization, oracle complexity is not sufficient to understand performance. Indeed, methods based on momentumbased variance reduction are known to be oracle-optimal in the single-worker setting under stronger smoothness assumptions [Cutkosky and Orabona, 2019, Arjevani et al., 2022]; however, as our results and Table 1 show, this does not automatically translate into time optimality in the parallel heterogeneous-worker setting. Time complexity behavior also depends on how optimization progress interacts with the worker-time model and the batch-size/time tradeoff. In this sense, a central novelty of the paper lies in the time complexity theory needed to understand this mechanism in the heterogeneous parallel regime. To support this perspective, we also established a new lower bound on the achievable time complexity under the same model. Together with our upper and lower bounds, this provides, to the best of our knowledge, the first theoretical understanding of the time complexity behavior of momentum-based variance reduction in heterogeneous distributed systems. At the same time, a gap between the current upper and lower bounds remains. Closing this gap appears to be nontrivial, and likely requires fundamentally different ideas rather than a routine extension of existing variance-reduction analyses.
Acknowledgments and Disclosure of Funding The research reported in this publication was supported by funding from King Abdullah University of Science and Technology (KAUST): i) KAUST Baseline Research Scheme, ii) CRG Grant ORFS-CRG12-2024-6460, and iii) Center of Excellence for Generative AI, under award number 5940.
References Alekh Agarwal and John C Duchi. Distributed delayed stochastic optimization. Advances in Neural Information Processing Systems, 24, 2011. (Cited on page 3.) Abdelkrim Alahyane, Céline Comte, Matthieu Jonckheere, and Éric Moulines. Optimizing asynchronous federated learning: A delicate trade-off between model-parameter staleness and update frequency. arXiv preprint arXiv:2502.08206, 2025. (Cited on page 3.) Zeyuan Allen-Zhu and Yang Yuan. Improved SVRG for non-strongly-convex or sum-of-non-convex objectives. In Maria Florina Balcan and Kilian Q. Weinberger, editors, Proceedings of The 33rd International Conference
13
on Machine Learning, volume 48 of Proceedings of Machine Learning Research, pages 1080–1089, New York, New York, USA, 20–22 Jun 2016. PMLR. URL https://proceedings.mlr.press/v48/allen-zhub16. html. (Cited on page 3.) Yossi Arjevani, Yair Carmon, John C Duchi, Dylan J Foster, Nathan Srebro, and Blake Woodworth. Lower bounds for non-convex stochastic optimization. Mathematical Programming, pages 1–50, 2022. (Cited on pages 2, 4, 5, 6, 7, 13, 27, 28, 30, and 31.) Mahmoud Assran, Arda Aytekin, Hamid Reza Feyzmahdavian, Mikael Johansson, and Michael G Rabbat. Advances in asynchronous parallel and distributed optimization. Proceedings of the IEEE, 108(11):2013– 2031, 2020. (Cited on page 3.) Ashok Cutkosky and Francesco Orabona. Momentum-based variance reduction in non-convex SGD. Advances in Neural Information Processing Systems, 32, 2019. (Cited on pages 2, 3, 5, 6, and 13.) Aaron Defazio and Leon Bottou. On the ineffectiveness of variance reduced optimization for deep learning. In Advances in Neural Information Processing Systems, volume 32, 2019. (Cited on page 3.) Aaron Defazio, Francis Bach, and Simon Lacoste-Julien. SAGA: A fast incremental gradient method with support for non-strongly convex composite objectives. Advances in Neural Information Processing Systems, 27, 2014. (Cited on page 3.) Cong Fang, Chris Junchi Li, Zhouchen Lin, and Tong Zhang. SPIDER: Near-optimal non-convex optimization via stochastic path integrated differential estimator. In NeurIPS Information Processing Systems, 2018. (Cited on page 3.) Saeed Ghadimi and Guanghui Lan. Stochastic first-and zeroth-order methods for nonconvex stochastic programming. SIAM Journal on Optimization, 23(4):2341–2368, 2013. (Cited on page 4.) Robert M. Gower, Peter Richtárik, and Francis Bach. Stochastic quasi-gradient methods: Variance reduction via Jacobian sketching, 2018. URL https://arxiv.org/abs/1805.02632. (Cited on page 3.) Rie Johnson and Tong Zhang. Accelerating stochastic gradient descent using predictive variance reduction. Advances in Neural Information Processing Systems, 26, 2013. (Cited on page 3.) Anastasiia Koloskova, Sebastian U Stich, and Martin Jaggi. Sharper convergence guarantees for asynchronous SGD for distributed and federated learning. Advances in Neural Information Processing Systems, 35: 17202–17215, 2022. (Cited on pages 3 and 4.) Jakub Konečný and Peter Richtárik. Semi-stochastic gradient descent methods, 2015. URL https://arxiv. org/abs/1312.1666. (Cited on page 3.) Dmitry Kovalev, Samuel Horvath, and Peter Richtarik. Don’t jump through hoops and remove those loops: SVRG and Katyusha are better without the outer loop, 2019. URL https://arxiv.org/abs/1901.08689. (Cited on page 3.) Zhize Li, Hongyan Bao, Xiangliang Zhang, and Peter Richtárik. PAGE: A simple and optimal probabilistic gradient estimator for nonconvex optimization. In International Conference on Machine Learning, pages 6286–6295. PMLR, 2021. (Cited on page 3.) Artavazd Maranjyan. First Provably Optimal Asynchronous SGD for Homogeneous and Heterogeneous Data. PhD thesis, 2025. (Cited on page 3.) Artavazd Maranjyan and Peter Richtárik. Ringleader ASGD: The first asynchronous SGD with optimal time complexity under data heterogeneity. In The Fourteenth International Conference on Learning Representations, 2026. URL https://openreview.net/forum?id=5wqTal0EuC. (Cited on page 3.) Artavazd Maranjyan, Omar Shaikh Omar, and Peter Richtárik. Mindflayer SGD: Efficient parallel SGD in the presence of heterogeneous and random worker compute times. In The 41st Conference on Uncertainty in Artificial Intelligence, 2025a. (Cited on page 3.) Artavazd Maranjyan, El Mehdi Saad, Peter Richtárik, and Francesco Orabona. ATA: Adaptive task allocation for efficient resource management in distributed machine learning. In International Conference on Machine Learning, 2025b. (Cited on page 3.) Artavazd Maranjyan, Alexander Tyurin, and Peter Richtárik. Ringmaster ASGD: The first asynchronous SGD with optimal time complexity. In International Conference on Machine Learning, 2025c. (Cited on pages 3 and 4.)
14
Konstantin Mishchenko, Francis Bach, Mathieu Even, and Blake E Woodworth. Asynchronous SGD beats minibatch SGD under arbitrary delays. Advances in Neural Information Processing Systems, 35:420–433, 2022. (Cited on pages 3 and 4.) Yurii Nesterov. Lectures on Convex Optimization, volume 137. Springer, 2018. (Cited on page 22.) Lam Nguyen, Jie Liu, Katya Scheinberg, and Martin Takáč. SARAH: A novel method for machine learning problems using stochastic recursive gradient. In The 34th International Conference on Machine Learning, 2017. (Cited on page 3.) Benjamin Recht, Christopher Re, Stephen Wright, and Feng Niu. HOGWILD!: A lock-free approach to parallelizing stochastic gradient descent. Advances in Neural Information Processing Systems, 24, 2011. (Cited on page 3.) Herbert Robbins and Sutton Monro. A stochastic approximation method. Annals of Mathematical Statistics, 22 (3):400–407, 1951. (Cited on page 2.) Rafał Szlendak, Alexander Tyurin, and Peter Richtárik. Permutation compressors for provably faster distributed nonconvex optimization. arXiv preprint arXiv:2110.03300, 2021. (Cited on page 4.) John Tsitsiklis, Dimitri Bertsekas, and Michael Athans. Distributed asynchronous deterministic and stochastic gradient optimization algorithms. IEEE Transactions on Automatic Control, 31(9):803–812, 1986. (Cited on page 3.) Alexander Tyurin. Tight time complexities in parallel stochastic optimization with arbitrary computation dynamics. In The Thirteenth International Conference on Learning Representations, 2025. URL https: //openreview.net/forum?id=cUN8lJB4rD. (Cited on pages 17 and 18.) Alexander Tyurin and Peter Richtárik. Optimal time complexities of parallel stochastic optimization methods under a fixed computation model. In Advances in Neural Information Processing Systems, volume 36, pages 16515–16577, 2023. (Cited on pages 1, 2, 3, 4, 6, 8, 9, 11, 13, 25, 27, 28, 29, 30, and 33.) Alexander Tyurin and Peter Richtárik. On the optimal time complexities in decentralized stochastic asynchronous optimization. Advances in Neural Information Processing Systems, 37, 2024. (Cited on page 3.) Alexander Tyurin, Kaja Gruntkowska, and Peter Richtárik. Freya PAGE: First optimal time complexity for large-scale nonconvex finite-sum optimization with heterogeneous asynchronous computations. Advances in Neural Information Processing Systems, 37, 2024a. (Cited on page 3.) Alexander Tyurin, Marta Pozzi, Ivan Ilin, and Peter Richtárik. Shadowheart SGD: Distributed asynchronous SGD with optimal time complexity under arbitrary computation and communication heterogeneity. Advances in Neural Information Processing Systems, 37, 2024b. (Cited on page 3.) Yidong Wu and Luo Luo. Optimal asynchronous stochastic nonconvex optimization under heavy-tailed noise, 2026. URL https://arxiv.org/abs/2601.19379. (Cited on page 3.) Huizhuo Yuan, Yifeng Liu, Shuang Wu, Zhou Xun, and Quanquan Gu. MARS: Unleashing the power of variance reduction for training large models. In Forty-second International Conference on Machine Learning, 2025. URL https://openreview.net/forum?id=NrcKQ3ASLZ. (Cited on pages 3, 11, and 12.) Dongruo Zhou, Pan Xu, and Quanquan Gu. Stochastic nested variance reduction for nonconvex optimization. Journal of Machine Learning Research, 21(103):1–63, 2020. URL http://jmlr.org/papers/v21/ 18-447.html. (Cited on pages 2, 3, and 5.)
15
Contents 1
2
Introduction
1
1.1
Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
1.2
Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
Problem Setup
3
2.1
4
Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
Background and Motivation
4
4
Rennala MVR
5
4.1
Algorithm Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
4.2
Connection to Rennala SGD . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
5 Theoretical Results
6
7
8
6
5.1
Iteration Complexity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
5.2
Time Complexity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
Lower Bound on Time Complexity
8
6.1
Protocol . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
6.2
Oracle Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
6.3
Algorithm Class . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
6.4
Function and Oracle Classes . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
6.5
Time Complexity Measure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
Experiments
11
7.1
Quadratic Benchmark . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
7.2
Neural Network Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
Conclusion
13
A Time Complexity Under Arbitrarily Varying Worker Speeds
17
A.1 Universal Computation Model . . . . . . . . . . . . . . . . . . . . . . . . . . . .
17
A.2 Time Complexity in the Universal Model . . . . . . . . . . . . . . . . . . . . . .
17
B Proof of Lemmas
18
B.1 Proof of Lemma 2.5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
B.2 Proof of Lemma 2.6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
B.3 Expected Similarity and Bounded Variance for the Minibatch Case . . . . . . . . .
19
C Why the Current Analysis Cannot Close the Time-Complexity Gap
20
D Proofs of Upper Bounds
21 16
D.1 Proof of Theorem 5.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
D.2 Proof of Theorem 5.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
25
E Proofs of Lower Bound
27
E.1 Proof of Theorem 6.7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . F Additional experiments
A
28 34
Time Complexity Under Arbitrarily Varying Worker Speeds
The fixed-computation model in Assumption 2.1 assumes that every worker processes gradients at a constant rate. This assumption is often too restrictive in practice: worker speeds may change over time due to interruptions, network interference, or temporary hardware slowdowns. To capture such effects, we now replace the constants {τi }ni=1 with time-varying computation rates. A.1
Universal Computation Model
We adopt the universal computation model of Tyurin [2025]. Assumption A.1 (Universal Computation Model). For each worker i ∈ [n], there exists a nonnegative function pi : R+ → R+ that is continuous almost everywhere. For any 0 ≤ T1 ≤ T2 , the number of stochastic gradients completed by worker i during the interval [T1 , T2 ] is $Z % T2 Ni (T1 , T2 ) := pi (s) ds . T1
Moreover, Assumption A.1 contains Assumption 2.1 as a special case: if pi (s) = 1/τi for all s ≥ 0, then T2 − T1 Ni (T1 , T2 ) = . τi A.2
Time Complexity in the Universal Model
Unlike the fixed-time model, the universal model does not generally lead to a closed-form expression for the total time. The natural guarantee is therefore implicit and expressed through the times at which enough gradients have been completed to finish the initialization and each subsequent iteration. For the universal-model upper bound, we assume that at each iteration boundary all workers are effectively reset to the idle state, so unfinished computations from the previous iteration do not contribute to the next one. Theorem A.2 (Time Complexity Under the Universal Computation Model). Assume the setup and parameter choices of Theorem 5.1, and define 24∆L̄ σ Kε := +√ . ε ε Let 2 6σ 6σ , B= √ , B0 = ε ε k and define the completion times {TMVR }k≥0 by ) ( n Z t X 0 pi (s) ds ≥ B0 , TMVR := min t ≥ 0 : i=1
( k := min TMVR
k−1 t ≥ TMVR :
0
$ Z % ) 1 t pi (s) ds ≥ B , k−1 2 TMVR i=1
n X
17
k ≥ 1.
Kε Then Rennala MVR (Algorithm 1) returns an ε–stationary point within at most TMVR seconds.
Proof. By Theorem 5.1, under the stated parameter choice, Rennala MVR reaches an ε–stationary point after at most Kε iterations. It therefore remains to upper bound the time needed to complete the initialization and these Kε iterations. 0 0 By definition of TMVR , by time TMVR the workers have collectively produced at least B0 stochastic 0 0 gradients at x . Hence the initialization step used to construct g 0 is completed no later than TMVR . k−1 k−1 Now fix k ≥ 1 and suppose that iteration k starts at time TMVR . Over the interval [TMVR , t], worker i can complete $Z % t pi (s) ds k−1 TMVR
stochastic gradients. One arrival in Rennala MVR at iteration k is a complete pair ∇f (xk−1 ; ξ), ∇f (xk ; ξ) , and therefore requires two stochastic-gradient computations on the same worker. Consequently, by time t, worker i can contribute at least % $ Z 1 t pi (s) ds k−1 2 TMVR complete gradient pairs. Summing over workers, the server has received at least $ Z % n X 1 t pi (s) ds k−1 2 TMVR i=1 k k pairs by time t. By the definition of TMVR , this quantity is at least B when t = TMVR . Hence k iteration k finishes no later than TMVR .
Applying the above argument inductively for k = 1, . . . , Kε , we conclude that all Kε iterations finish Kε by time TMVR . Combining this with Theorem 5.1 proves the claim. Comparison with Rennala SGD. The corresponding universal-model bound for Rennala SGD follows from Tyurin [2025, Theorem 5.3], specialized to our notation. In the regime ε < σ 2 KSGD considered in Theorem 5.1, Rennala SGD reaches an ε–stationary point within TSGD seconds, where 24L∆ KSGD := ε and 0 := 0, TSGD ( $Z % ) n t X σ2 k−1 k := min t ≥ TSGD TSGD : , k ≥ 1. pi (s) ds ≥ k−1 ε TSGD i=1 Thus, the universal model preserves the same qualitative trade-off as the fixed-time model. Each √ iteration of Rennala MVR waits for B = Θ(σ/ ε) complete pairs, whereas each iteration of Rennala SGD waits for Θ(σ 2 /ε) stochastic gradients. For sufficiently small ε, the number of stochastic gradients that must be completed before each Rennala MVR update √ is therefore asymptotically smaller. The trade-off is that Rennala MVR requires Kε = Θ(L̄∆/ε + σ/ ε) update rounds, whereas Rennala SGD requires KSGD = Θ(L∆/ε) rounds. Since the universal-model bounds are implicit, no unconditional ordering between the two total times follows without additional structure on the rate functions {pi }ni=1 . Nevertheless, whenever L̄ = O(L) and the reduction in per-round waiting time dominates, the same mechanism that improves the fixed-time bound can also favor Rennala MVR under time-varying worker speeds.
B
Proof of Lemmas
This appendix provides complete proofs of all lemmas referenced in the main text, including both those stated explicitly and auxiliary results needed for the theorem proofs. 18
B.1
Proof of Lemma 2.5
Let us first restate the lemma. Lemma 2.5. Mean-squared smoothness (Assumption 2.4) implies smoothness of f with constant L̄, i.e. ∥∇f (x) − ∇f (y)∥ ≤ L̄ ∥x − y∥ ,
∀ x, y ∈ Rd .
Proof. Using Jensen’s inequality, we have 2
2
∥∇f (x) − ∇f (y)∥ = ∥E [∇f (x; ξ) − ∇f (y; ξ)]∥ h i 2 2 ≤ E ∥∇f (x; ξ) − ∇f (y; ξ)∥ ≤ L̄2 ∥x − y∥ . Hence, taking square roots, ∥∇f (x) − ∇f (y)∥ ≤ L̄∥x − y∥ .
B.2
Proof of Lemma 2.6
Again restating the lemma first. Lemma 2.6. Mean-squared smoothness (Assumption 2.4) implies expected similarity with the same constant L̄: h i 2 2 E ∥∇f (x; ξ) − ∇f (y; ξ) − (∇f (x) − ∇f (y))∥ ≤ L̄2 ∥x − y∥ , ∀ x, y ∈ Rd .
Proof. This follows by doing some simple algebra E ∥∇f (x; ξ) − ∇f (y; ξ) − (∇f (x) − ∇f (y))∥2 = E ∥∇f (x; ξ) − ∇f (y; ξ)∥2 − ∥∇f (x) − ∇f (y)∥2 ≤ E ∥∇f (x; ξ) − ∇f (y; ξ)∥2 ≤ L̄2 ∥x − y∥2 .
B.3
Expected Similarity and Bounded Variance for the Minibatch Case
We now prove another lemma that will be used in the iteration complexity analysis. Lemma B.1. The following bounds hold: 2 B B X X 1 1 E ∇f (xk+1 ; ξ k,j ) − ∇f (xk+1 ) + ∇f (xk ) − ∇f (xk ; ξ k,j ) B j=1 B j=1 ≤
L̄2 k+1 2 x − xk , B
and
2
B X
1 E ∇f (xk ; ξ k,j ) − ∇f (xk ) B j=1
19
σ2 . ≤ B
Proof. Using Lemma 2.6 and independence of ξ k,j : 2 B B X X 1 1 E ∇f (xk+1 ; ξ k,j ) − ∇f (xk+1 ) + ∇f (xk ) − ∇f (xk ; ξ k,j ) B j=1 B j=1 B
=
i 1 X h 2 E ∇f (xk+1 ; ξ k,j ) − ∇f (xk+1 ) + ∇f (xk ) − ∇f (xk ; ξ k,j ) 2 B j=1 B
1 X 2 k+1 L̄2 k+1 2 2 ≤ 2 − xk = x − xk . L̄ x B j=1 B Similarly, using Assumption 2.3 and independence: 2 B B i 1 X h 2 1 X E ∇f (xk ; ξ k,j ) − ∇f (xk ) = 2 E ∇f (xk ; ξ k,j ) − ∇f (xk ) B j=1 B j=1 B
≤
C
1 X 2 σ2 σ = . 2 B j=1 B
Why the Current Analysis Cannot Close the Time-Complexity Gap
Our time lower bound (cf. Theorem 6.7) implies, up to universal constants, a requirement of the form 2 L̄∆ σ √ +1 T mtime (B, ∆, L̄, σ 2 ) ≥ Ω , (6) ε σ ε where T (B) = min m∈[n]
m X
!−1 τi−1
(B + m)
i=1
is an upper bound on the time needed to obtain B stochastic gradients asynchronously. On the other hand, the time upper bound we obtain in this paper for Rennala MVR has the schematic structure Tup ≲ T (B0 ) + K T (B),
(7)
where B0 is the initialization batch size, B is the per-iteration batch size, and K is the number of server updates (iterations). If one aims to align (7) with the lower-bound scaling (6) uniformly over heterogeneous systems (within our timing model), then the structure of T (·) can force a specific scaling of (B, K); a simple way to see this is to consider the homogeneous case τi ≡ τ , for which T (B) = min m∈[n]
τ (B + n) τ τ (B + m) = = B+τ . m n n
In particular, even in this “linear” regime, T (B) has a nonzero per-round overhead τ , and hence τ K T (B) = KB + τ K, n while 2 2 L̄∆ σ L̄∆ τ σ2 τ L̄∆ σ L̄∆ √ +1 T √ +1 √ +1 √ +1 . = +τ = +τ ε n ε n σ ε ε σ ε σ ε σ ε Thus, matching the lower-bound scaling in this example forces the additive contributions to be of the same order, which yields L̄∆ K ≍ √ +1, (8) σ ε 20
2
and then matching the leading linear-in-B terms yields B ≍ σε (with B0 of the same order as well). However, our iteration-complexity bound (11) shows that (8) cannot be achieved by Rennala MVR under the stepsize restriction (10) used in our analysis. Let us enforce the batching suggested by the 2 lower bound and set B0 = B = c σε for a sufficiently large constant c. Then (11) gives K−1
i 2∆ 2pσ 2 1 2σ 2 1 X h 2 E ∇f (xk ) ≤ + + . K γK B pK B0 k=0
With B = B0
2 = cσ /ε, the noise terms become
2pσ 2 2p = ε B c and
2 1 2σ 2 = ε. pK B0 cpK By choosing c large enough (and, say, any constant p ≤ 1), these two terms can be made ≤ ε/3 2∆ provided K ≳ 1/p. The remaining requirement comes from the descent term γK . Crucially, 2 the stepsize √ constraint (10) enforces γ = O(1/L̄) even when B is as large as σ /ε (the factor (1 − p)/ pB only makes γ smaller). Therefore, to make ε 2∆ ≤ , γK 3 we must have
L̄∆ . ε √ In the high-noise regime σ > ε, the lower-bound-implied target K≥Ω
K≍
(9)
L̄∆ √ σ ε
√
is smaller than L̄∆/ε by a factor σ/ ε. Hence, (9) rules out the simultaneous choice (8) within our current upper-bound analysis. In summary, while the lower bound (6) suggests that optimal time would be attained by producing gradients in batches of size σ 2 /ε and performing only √ Θ L̄∆/(σ ε) update rounds, the smoothness-limited descent mechanism behind (11) forces Ω(L̄∆/ε) update rounds (under (10)) regardless of how large we take B. This creates a gap relative to (6) for the particular Rennala MVR analysis developed in this paper, which explains why our time upper bound does not match the lower bound in full generality. Although one can argue that the remaining gap may be due to our proof technique, we believe the bounds are reasonably tight: in particular, when we ignore the timing model and count only oracle queries, our analysis recovers the optimal oracle complexity of MVR.
D
Proofs of Upper Bounds
In this section, we prove the upper bounds on the iteration and time complexity. D.1
Proof of Theorem 5.1
We start with the iteration complexity. Theorem 5.1 (Iteration Complexity of Algorithm 1). Let Assumptions 2.2, 2.3, and 2.4 hold. Fix parameters p ∈ (0, 1], integers B0 ≥ 1 and B ≥ 1. Consider Algorithm 1 with the modification that the initialization g 0 uses an unbiased estimate of ∇f (x0 ) using a minibatch of size B0 , while for all k ≥ 0 the iteration-k minibatch has size B. Choose the stepsize γ≤
1
1−p 2L̄ 1 + √ pB
21
.
(10)
Then for any K ≥ 1, K−1
i 1 X h 2∆ 2pσ 2 1 2σ 2 2 E ∇f (xk ) ≤ + + , K γK B pK B0
(11)
k=0
where ∆ := f (x0 ) − f ∗ . Fix ε > 0 and assume ε < σ 2 and ε < 2L̄∆. Choose 1 , γ= 4L̄
√
6σ 2 B0 := , ε
ε p := , σ
B :=
6pσ 2 ε
6σ = √ , ε
(12)
then the iterates of Algorithm 1 satisfy K−1
i 1 X h 2 E ∇f (xk ) ≤ ε, K k=0
for
K≥
24∆L̄ σ +√ . ε ε
Proof. Since Assumption 2.4 implies that f is L̄–smooth (see Lemma 2.5), the standard smoothness inequality (e.g., [Nesterov, 2018]) gives
f (xk+1 ) ≤ f (xk ) + ⟨∇f (xk ), xk+1 − xk ⟩ +
L̄ k+1 2 x − xk 2
L̄ k+1 2 x − xk 2 γ 1 L̄ k+1 2 2 2 2 g k − ∇f (xk ) − ∇f (xk ) − xk+1 − xk + x − xk 2 2γ 2 L̄ 1 γ 2 2 2 ∇f (xk ) − − xk+1 − xk + g k − ∇f (xk ) . 2γ 2 2
= f (xk ) + ⟨∇f (xk ), −γg k ⟩ + γ 2 γ = f (xk ) − 2 = f (xk ) +
Subtract f ∗ and take E [ · | Fk ], where Fk := σ(x0 , . . . , xk , g 0 , . . . , g k ): γ 2 E f (xk+1 ) − f ∗ | Fk ≤ f (xk ) − f ∗ − ∇f (xk ) 2 γ k 1 L̄ 2 2 − − xk+1 − xk + g − ∇f (xk ) . 2γ 2 2 22
h i 2 We need to control the variance term above. So, let us study E g k+1 − ∇f (xk+1 ) Fk . Using bias–variance decomposition: h i 2 E g k+1 − ∇f (xk+1 ) Fk 2 B B X X 1 1 Fk = E ∇f (xk+1 ; ξ k,j ) − ∇f (xk+1 ) + (1 − p) g k − ∇f (xk ; ξ k,j ) B j=1 B j=1 2
= (1 − p)2 g k − ∇f (xk ) 2 B B X X 1 1 ∇f (xk+1 ; ξ k,j ) − ∇f (xk+1 ) + (1 − p) ∇f (xk ) − ∇f (xk ; ξ k,j ) + E B j=1 B j=1
Fk
2
= (1 − p)2 g k − ∇f (xk ) " ! B 1 X k+1 k,j k+1 +E p ∇f (x ; ξ ) − ∇f (x ) B j=1 # ! 2 B B 1 X 1 X k+1 k,j k+1 k k k,j Fk ∇f (x ; ξ ) − ∇f (x ) + ∇f (x ) − ∇f (x ; ξ ) + (1 − p) B j=1 B j=1 2 B 2 1 X ≤ (1 − p)2 g k − ∇f (xk ) + 2p2 E ∇f (xk+1 ; ξ k,j ) − ∇f (xk+1 ) Fk B j=1
B B 1 X 1 X + 2(1 − p)2 E ∇f (xk+1 ; ξ k,j ) − ∇f (xk+1 ) + ∇f (xk ) − ∇f (xk ; ξ k,j ) B j=1 B j=1
2
Fk .
Using Lemma B.1, we get h 2 E g k+1 − ∇f (xk+1 )
i 2 2 2(1 − p)2 L̄2 k+1 2 2p σ 2 Fk ≤ (1−p)2 g k − ∇f (xk ) + x − xk . + B B Define the Lyapunov function: γ 2 Ψk = f (xk ) − f ∗ + g k − ∇f (xk ) . 2p We have: 1 γ L̄ γ k 2 2 2 ∇f (xk ) − xk+1 − xk + g − ∇f (xk ) E [ Ψk+1 | Fk ] ≤ f (xk ) − f ∗ − − 2 2γ 2 2 γ 2(1 − p)2 L̄2 k+1 2p2 σ 2 2 2 + (1 − p)2 g k − ∇f (xk ) + + x − xk 2p B B L̄ γ(1 − p)2 L̄2 γ 1 2 2 ∇f (xk ) − − − xk+1 − xk = f (xk ) − f ∗ − 2 2γ 2 pB γ γpσ 2 2 p + (1 − p)2 g k − ∇f (xk ) + . + 2p B Using p + (1 − p)2 ≤ 1, we get γ 2 E [ Ψk+1 | Fk ] ≤ Ψk − ∇f (xk ) − 2
1 L̄ γ(1 − p)2 L̄2 − − 2γ 2 pB
Next to ensure that the coefficient before xk+1 − xk 1 L̄ ≥ 4γ 2
and
2
xk+1 − xk
2
+
γpσ 2 . B
is non-negative, it is sufficient to have
1 γ(1 − p)2 L̄2 ≥ , 4γ pB 23
which is equivalent to γ≤ Hence, taking γ ≤
1 2L̄
γ≤
and
1 1−p 2L̄ √ pB 2
ensures 1 − L̄ − γ(1−p) L̄ 1 1−p 2γ 2 pB 2L̄ 1+ √ pB
2
.
≥ 0. Thus, we obtained
γ γpσ 2 2 ∇f (xk ) + . 2 B Taking full expectation and telescoping from k = 0 to K − 1 yields K−1 i γpσ 2 i X γ h γ h 2 k 2 − E ∇f (x ) ≤ ∆ + E g 0 − ∇f (x0 ) . 2 B 2p E [ Ψk+1 | Fk ] ≤ Ψk −
k=0
Hence, K−1
i 1 X h 2∆ 2pσ 2 1 2σ 2 2 E ∇f (xk ) ≤ + + . K γK B pK B0 k=0
Fix ε > 0. We assume ε < σ 2 . Otherwise, we are in the low-noise regime, hence we can set p = 1 and use SGD with minibatch size B, which reaches the target accuracy in O L̄∆ stochastic ε gradients. Also assume ε < 2L̄∆. Otherwise, since ∥∇f (x0 )∥2 ≤ 2L̄∆ ≤ ε, the initial point x0 is already ε–stationary. Recall that we have established K−1
i 1 X h 2∆ 2pσ 2 1 2σ 2 2 ≤ E ∇f (xk ) + + . K γK B pK B0 k=0
To make sure the right hand side is bounded by ε, we bound each term on the right-hand side by ε/3. We choose the initialization batch size 6σ 2 B0 := , ε
(13)
which gives 1 ε 1 2σ 2 · . ≤ pK B0 pK 3 To ensure this term is at most ε/3, we require 1 . p
K ≥ We set
6pσ 2 B := , ε
(14)
(15)
which directly ensures 2pσ 2 ε ≤ . B 3 We choose
√
ε , σ
(16)
6σ B= √ . ε
(17)
p := Substituting (16) into (15) yields
24
Moreover, this choice of p gives √ 1−p 1 1−p 1−p ε √ ≤p =√ · · p σ pB 6 p · 6pσ 2 /ε 1 1 = √ (1 − p) ≤ √ . 6 6 The above bound ensures that the step size γ :=
1 4L̄
(18)
satisfies the required condition, since 1 1−p ≤ 2L̄ 1 + √ < 4L̄ . 2L̄ 1 + √ pB 6 With γ = 1/4L̄, bounding the first term by ε/3 requires K ≥ Combining with (14) and recalling p =
√
24∆L̄ . ε
(19)
ε/σ , we take
K ≥
24∆L̄ σ +√ . ε ε
(20)
The parameter choices (13), (17), (16), (18), and (20) together guarantee K−1
i 1 X h 2 E ∇f (xk ) ≤ ε. K k=0
D.2
Proof of Theorem 5.2
Before proving the theorem, we first prove a simple lemma. A result of this type was used in [Tyurin and Richtárik, 2023, Theorem 7.5]; here we provide a simpler proof. Lemma D.1. Suppose we have n workers with computation times as in Assumption 2.1. Let B ∈ {1, 2, . . .}. Consider any iteration boundary at which the server starts collecting new stochastic gradients. Then, the time needed to collect B stochastic gradients is at most !−1 m X 1 2T (B) = 2 min (B + m) . τ m∈[n] i=1 i Proof. Fix t > 0. For worker i, in the worst case the worker finishes and sends a gradient just before the iteration starts. Therefore, the first interval of length τi after the iteration starts may not produce a new gradient usable in this iteration. Hence, within t seconds from the iteration start, worker i can contribute at least t − τi t ni (t) := max , 0 = max − 1, 0 τi τi new gradients. Set ∗
j = arg min j∈[n]
j X 1 i=1
25
τi
!−1 (B + j) ,
then, since max{u, 0} ≥ u, we have ni (t) ≥ n X
∗
ni (t) ≥
i=1
j X
j k t τi
− 1, and therefore ∗
∗
ni (t) ≥
i=1
j X t
τi
i=1
−1
≥
j X t i=1
τi
∗
−2
= t
j X 1
τ i=1 i
− 2j ∗ ,
where we used ⌊a⌋ ≥ a − 1 for all a ≥ 0. Thus, if t satisfies ∗
t
j X 1
τ i=1 i
− 2j ∗ ≥ B,
then by time t the server can collect at least B gradients from workers 1, . . . , j ∗ . Choosing −1 ∗ j X 1 t = 2 (B + j ∗ ) τ i=1 i makes the left-hand side equal to 2B ≥ B, concluding the proof. Theorem 5.2 (Time Complexity of Algorithm 1). Assume the setup of Theorem 5.1 and our distributed time model: there are n workers with per-sample computation times 0 < τ1 ≤ · · · ≤ τn . Assume that on worker i, computing the pair ∇f (xk ; ξ), ∇f (xk+1 ; ξ) takes 2τi seconds. Then, the initialization (collecting B0 single gradients) takes at most !−1 m X 1 2 · min (B0 + m) τ m∈[n] i=1 i
(21)
seconds, and each iteration (collecting B gradient-pairs) takes at most !−1 m X 1 4 · min (B + m) τ m∈[n] i=1 i
(22)
seconds. In particular, with the parameter choice of Theorem 5.1, Algorithm 1 returns an ε–stationary point within !−1 !−1 m m X X 1 1 (B0 + m) + 4K · min (B + m) (23) TMVR ≤ 2 · min τ τ m∈[n] m∈[n] i=1 i i=1 i seconds, where 6σ 2 B0 = , ε
6σ B= √ , ε
K≥
24∆L̄ σ +√ . ε ε
Thus, TMVR = O
σ L̄∆ √ + ε ε
min m∈[n]
m X 1 i=1
τi
!−1
!−1 ! m X σ 1 σ2 √ + m + min +m . τ ε ε m∈[n] i=1 i
Proof. The initialization bound follows from Lemma D.1 applied with B = B0 and computation times τi . For each iteration, the server collects B arrivals, and each arrival produced by worker i is a gradientpair and takes 2τi seconds. Applying Lemma D.1 with computation times 2τi gives that one iteration takes at most !−1 !−1 m m X X 1 1 2 min (B + m) = 4 min (B + m). 2τi τ m∈[n] m∈[n] i=1 i=1 i Summing the initialization time and K iteration times yields the desired bound. 26
E
Proofs of Lower Bound
We start by introducing several definitions and notations that will be used throughout this section. For a vector x ∈ Rd , let support(x) := {i ∈ [d] | xi ̸= 0} and x≥i := (xi , . . . , xd ) ∈ Rd−i+1 . For α ∈ [0, 1) define the progress progα (x) := max{ i ∈ {0, 1, . . . , d} | |xi | > α },
with x0 ≡ 1 .
(24)
Our goal is to construct functions and stochastic oracles such that, under zero-respecting algorithms, each completed oracle reply can activate at most one new coordinate. Definition E.1 (First-order zero-chain). A differentiable function F : RT → R is a (first-order) zero-chain if, for all x ∈ RT , prog0 ∇F (x) ≤ prog0 (x) + 1. (25) Intuition. In the noiseless case g(x, ξ) ≡ ∇F (x), (25) implies that a zero-respecting algorithm can reveal at most one new coordinate per completed oracle reply, so its progress is at most linear in the number of replies. Definition E.2 (Probability-p zero-chain). A stochastic mapping g : RT × D → RT is a probability-p zero-chain if, for all x ∈ RT , Pξ prog0 g(x, ξ) = prog 1 (x) + 1 ≤ p, (26) 4 Pξ prog0 g(x, ξ) > prog 1 (x) + 1 = 0. (27) 4
Setup We use the same deterministic chain function as in the works by Tyurin and Richtárik [2023] and Arjevani et al. [2022]. Define FT : RT → R by FT (x) := −Ψ(1) Φ(x1 ) +
T X
Ψ(−xi−1 ) Φ(−xi ) − Ψ(xi−1 ) Φ(xi ) ,
i=2
where 0,
t ≤ 12 , 1 1 exp 1 − (2t−1)2 , t > 2 , Z t 1 2 √ Φ(t) = e e− 2 τ dτ.
Ψ(t) =
−∞
Lemma E.3 (Properties of FT , cf. Lemma 2). There are absolute constants ∆0 = 12, ℓ1 = 152, γ∞ = 23 such that: 1. FT (0) − inf x FT (x) ≤ ∆0 · T . 2. ∇FT is ℓ1 –Lipschitz in ℓ1 . 3. For all x, ∥∇FT (x)∥∞ ≤ γ∞ . 4. prog0 ∇FT (x) ≤ prog1/2 (x) + 1. 5. If prog1 (x) < T , then ∥∇FT (x)∥ ≥ ∇prog0 (x)+1 FT (x) > 1. Next we define the estimator. Let Γ : R → R be smooth, non-decreasing, and Lipschitz, with Γ(t) = 1 for t ≥ 12 .
Γ(t) = 0 for t ≤ 41 , 27
(28)
Define, for each i,
Θi (x) := Γ 1 − Γ |x≥i |
(29) 2
so that 1{i > prog1/4 (x)} ≤ Θi (x) ≤ 1{i > prog1/2 (x)}.
(30)
A concrete choice is obtained by the “integrated bump”: Rt t ≤ 41 or t ≥ 12 , 0, Λ(τ ) dτ 1/4 Γ(t) = R 1/2 , Λ(t) = 1 ′ ) dτ ′ , 14 < t < 21 . exp − Λ(τ 1 1 1/4 100 (t − 4 ) ( 2 − t)
(31)
This Γ satisfies: Γ ∈ C ∞ , 0 ≤ Γ′ (t) ≤ 6, and |Γ′′ (t)| ≤ 128. Define the smoothed estimator ḡT (x, ξ) i := ∇i FT (x) νi (x, ξ),
νi (x, ξ) := 1 + Θi (x)
ξ p
−1 ,
ξ ∼ Bernoulli(p). (32)
Lemma E.4 (Mean-squared smooth estimator (Lemma 4 in Arjevani et al. [2022])). ḡT is a probability-p zero-chain, is unbiased for ∇FT , and there exist constants ς = 23 and ℓ̄1 = 328 such that for all x, y ∈ RT , h i 1−p 2 E ḡT (x, z) − ∇FT (x) ≤ ς2 , p h i ℓ̄ 2 2 E ḡT (x, z) − ḡT (y, z) ≤ 1 ∥x − y∥2 . p E.1
Proof of Theorem 6.7
Let us first restate the theorem before proving it. Theorem 6.7. Fix ∆ > 0, L̄ > 0, σ 2 > 0, 0 < ε < c′ L̄∆, an integer B ≥ 1, and n workers with batch-time functions {τi (·)}ni=1 that satisfy Assumption 6.1. Write τi := τi (1) and assume 0 < τ1 ≤ · · · ≤ τn . Let the chain constants ∆0 , ℓ1 , γ∞ be as in Lemma E.3 and the estimator constants ς, ℓ̄1 as in Lemma E.4. Define 2ες 2 ℓ1 √ p := min , 1 , L := L̄ p (≤ L̄), 2 σ ℓ̄1 $ % ℓ1 √ L∆ λ := 2ε, T := . L 2 ∆0 ℓ1 ε 2
,L̄,B Then there exist f ∈ F∆,L̄ and an oracle class O ∈ Oτσ1 (·),...,τ (f ) such that, under Protocol 2, n (·)
2
mtime (B, ∆, L̄, σ ) ≥ c ·
√
L̄∆ min { ε/σ, 1} +1 ε
min m∈[n]
m X 1 i=1
τi
!−1
σ2 +m . ε
Proof. Before presenting the proof, we briefly outline the argument. We first follow the construction by Arjevani et al. [2022], and then derive the time complexity lower bound using the time protocol analysis from Tyurin and Richtárik [2023]. We work under Protocol 2 with batch sizes bounded by B. By the definition of zero-respecting algorithms, for any interaction r ≥ 1 and any k ∈ [B], [ [ ′ support x(r,k) ⊆ support g (s,k ) . s<r k′ ∈[B]
28
At r = 1 the union on the right is empty, hence x(1,k) = 0 for all k ∈ [B]. Fix parameters ∆ > 0, accuracy ε > 0, and let L ≤ L̄ be chosen below. Define x ℓ1 √ ∆ Lλ2 L∆ , λ= 2ε, T = f (x) = FT = . ℓ1 λ L ∆0 (Lλ2 /ℓ1 ) 2∆0 ℓ1 ε
(33)
By Lemma E.3 (i),(ii), f ∈ F∆,L̄ . Moreover, ∇f (x) = Lλ ℓ1 ∇FT (x/λ). Let the Bernoulli parameter p ∈ (0, 1] be chosen below. Define x Lλ Lλ ∇f (x; ξ) j = ḡT (x/λ, ξ) j = 1 + Θj (x/λ) pξ − 1 ∇ j FT ℓ1 ℓ1 λ = ∇j f (x) 1 + Θj (x/λ) pξ − 1 , ξ ∼ Bernoulli(p). (34) By Lemma E.4, we have E[∇f (x; ξ)] = ∇f (x), and E ∇f (x; ξ) − ∇f (x)
2
≤
Lλ ℓ1
2
ς2
1−p . p
Choosing 2ες 2 p = min , 1 σ2
(35)
makes the variance ≤ σ 2 . Again by Lemma E.4, E ∇f (x; ξ) − ∇f (y; ξ)
2
≤
Thus taking ℓ1 √ ℓ1 L= L̄ p = L̄ · min ℓ̄1 ℓ̄1
ℓ̄1√ L ℓ1 p
2
∥x − y∥2 .
( √ ) ς 2ε , 1 ≤ L̄ σ
(36)
ensures that ∇f belongs to the oracle class in Definition 6.5. We now derive the time lower bound following Tyurin and Richtárik [2023]. From Lemma E.3 (v) we have that if prog1 (u) < T , then ∥∇FT (u)∥ > 1. Using the monotonicity of progα in α and prog0 (x/λ) = prog0 (x), we obtain, for any x ∈ RT , prog0 (x) < T =⇒ ∥∇f (x)∥2 =
Lλ ∇FT (x/λ) ℓ1
2 2
= 2ε ∥∇FT (x/λ)∥ > 2ε.
Equivalently, ∥∇f (x)∥2 > 2ε 1{prog0 (x) < T }
for all x ∈ RT .
From Lemma E.3(iv), we also have Lλ ∇FT (x/λ) = prog0 ∇FT (x/λ) ℓ1 ≤ prog1/2 (x/λ) + 1 ≤ prog1/4 (x/λ) + 1.
prog0 (∇f (x)) = prog0
Thus, for indices i > prog1/4 (x/λ) + 1 we have ∇i f (x) = 0, and using (34), [∇f (x; ξ)]i = 0,
if i > prog1/4 (x/λ) + 1.
Moreover, by (30), Θi (x/λ) = 1 for i = prog1/4 (x/λ) + 1, and hence ξ [∇f (x; ξ)]i = ∇i f (x) , p
if i = prog1/4 (x/λ) + 1.
29
(37)
Since prog1/4 (x/λ) ≤ prog0 (x), the only way to activate a new coordinate i > prog0 (x) is to have ξ = 1, and this can happen only at i = prog1/4 (x/λ) + 1 ≤ prog0 (x) + 1. Therefore, ∇f is a probability-p zero-chain in the sense of Definition E.2. Now let {x(r,k) }r≥1,k∈[B] be the query points of a zero-respecting algorithm A ∈ Azr interacting with the oracles under Protocol 2, and let g (r,k) = ∇f (x(r,k) ; ξ (r) ) denote the completed stochastic gradients associated with interaction r. By the zero-respecting property and the probability-p zero-chain structure, the process maxs≤r, k∈[B] prog0 (x(s,k) ) can increase by at most one per completed oracle reply, and this increase occurs only when ξ (r) = 1, which happens with probability at most p (the same ξ (r) is used for all k ∈ [B]). Moreover, since τi (1) ≤ τi (k) for any batch size k ≤ B, replacing the batch-time function by the smaller delay τi (1) can only make the oracles faster. Therefore, any lower bound proved with delays τi = τi (1) also applies to the original model. Therefore, all assumptions of Lemma D.2 in Tyurin and Richtárik [2023] are satisfied (with prog(x) ≡ prog0 (x) and delays τi ), and we can invoke it directly. In particular, for any δ ∈ (0, 1) and any time t satisfying !−1 m X 1 T 1 1 1 t ≤ + log +m , min (38) 24 2 δ m∈[n] i=1 τi p we have, with probability at least 1 − δ, that prog0 (x(r,k) ) < T for all queried points (r, k) whose replies have been received by time t. Combining this with (37) yields inf
∥∇f (x(r,k) )∥2 > 2ε
with probability at least 1 − δ,
(r,k)∈St
and hence
E
(r,k)
∥∇f (x
inf
)∥
2
≥ 2ε(1 − δ).
(r,k)∈St
Choosing δ = 12 gives
E
(r,k)
∥∇f (x
inf
2
)∥
>ε
(r,k)∈St
whenever (38) holds with δ = 21 . Finally, substituting T from (33) and recalling the choice of L in (36), we obtain that it is necessary to have !−1 √ m X 1 L̄∆ p 1 1 1 −1 , (39) min +m t ≥ 24 m∈[n] τ p 2 2∆0 ℓ̄1 ε i=1 i to ensure
E
(r,k)
∥∇f (x
inf
)∥
2
≤ ε.
(r,k)∈St √ L̄∆ p
√ L̄∆ p
As in Arjevani et al. [2022] we consider the cases 2∆0 ℓ̄1 ε ≥ 5 and 2∆0 ℓ̄1 ε < 5 separately. In the former case, using ⌊x⌋/2 − 1 ≥ x/4, we get !−1 √ m X 1 L̄∆ p 1 1 t ≥ min +m 24 8∆0 ℓ̄1 ε m∈[n] i=1 τi p √ !−1 m X 1 L̄∆ min{ σ2ες , 1} 1 σ2 ≥ min + m . 24 τ 2ες 2 m∈[n] 8∆0 ℓ̄1 ε i=1 i Choosing c′ = 1/(40ℓ̄1 ∆0 ) implies ε ≤ L̄∆ 8 , hence the conditions of Lemma E.5 hold, and combining both bounds yields the desired result. j √ k L̄∆ p In the latter case, 12 2∆0 ℓ̄1 ε − 1 ≤ 1, and having ε<
L̄∆ 40ℓ̄1 ∆0 30
precludes the option p = 1. Hence, the right-hand side in (39) is smaller than the lower bound in Lemma E.5 up to a universal constant. This completes the proof. Here we prove the lemma that was used in the proof above. Lemma E.5. Assume ε ≤ L̄∆/8. Consider Protocol 2 with n workers and batch-time functions {τi (·)}ni=1 satisfying Assumption 6.1. Write τi := τi (1) and assume 0 < τ1 ≤ · · · ≤ τn . Then there exist functions {fs }s∈{−1,+1} ⊂ F∆,L̄ and, for each s ∈ {−1, +1}, a collection of oracles 2
,L̄,B and distributions ((O1 , . . . , On ), (D1 , . . . , Dn )) ∈ Oτσ1 (·),...,τ (fs ) such that for every algorithm n (·) A ∈ A run under the protocol, for any time t ≥ 0, σ2 L̄∆ (r,k) 2 max Es inf ∥∇fs (x )∥ ≥ min , , (40) 64 N (t) 8 s∈{−1,+1} (r,k)∈St
where Es denotes expectation under instance s, St is the set of query-point indices whose corresponding gradients have been returned by time t, and n X t N (t) := , (41) τi i=1 2
σ with the convention σ 2 /(64 N (t)) := +∞ when N (t) = 0. Consequently, if t satisfies N (t) ≤ 64ε , then (r,k) 2 max Es inf ∥∇fs (x )∥ ≥ ε. (42) s∈{−1,+1}
(r,k)∈St
In particular, any algorithm satisfying max Es s
(r,k)
inf (r,k)∈St
∥∇fs (x
2
)∥
≤ε
must have t ≥ c0 min
m∈[n]
m X 1 i=1
τi
!−1
σ2 +m , ε
(43)
for some absolute constant c0 . Proof. The argument is the same as Arjevani et al. [2022, Lemma 11], with the sample budget T replaced by the maximal number of independent oracle draws available by time t, namely N (t). p Fix r ∈ 0, 2∆/L̄ and define, for s ∈ {−1, +1}, σ2 L̄ s Pξ := N rs, 2 , f (x, ξ) := ∥x∥2 − 2ξx1 + r2 , fs (x) := Eξ∼Pξs [f (x, ξ)]. 2 L̄ Let θs := (rs, 0, . . . , 0). Then fs (x) =
L̄ ∥x − θs ∥2 , 2
so fs is L̄–smooth and fs (0) − inf x fs (x) = L̄2 r2 ≤ ∆, hence fs ∈ F∆,L̄ . For each worker i ∈ [n], set Di := Pξs and let the stochastic gradient mapping be ∇f (x; ξ) := ∇x f (x, ξ) = L̄(x − ξe1 ), so that Eξ [∇f (x; ξ)] = ∇fs (x), h i 2 Eξ ∥∇f (x; ξ) − ∇fs (x)∥ = σ 2 , h i 2 2 Eξ ∥∇f (x; ξ) − ∇f (y; ξ)∥ = L̄2 ∥x − y∥ . 2
,L̄,B Thus, the corresponding oracles Oi = Oτfi (·),B (defined in (5)) belong to Oτσ1 (·),...,τ (fs ). n (·)
31
In (5), each completed reply uses a single sample ξ ∼ Di shared across the entire returned batch. Hence, a completed reply contributes at most one independent sample, regardless of the chosen batch size b ≤ B. Since τi (·) is nondecreasing (Assumption 6.1), using b > 1 cannot increase the number of independent samples received by time t and can only increase completion times. Therefore, for the purpose of lower bounding the obtainable information by time t, we may restrict attention to unit batches and work with τi = τi (1), which yields the sample budget N (t) in (41). Let S be uniform on {−1, +1}, indicating which instance is selected. Conditioned on S = s, each completed reply reveals one independent draw from Pξs (indeed, from any returned gradient g = L̄(x − ξe1 ) we recover ξ = x1 − g1 /L̄ exactly). By time t, worker i can complete at most ⌊t/τi ⌋ such replies, hence the algorithm can receive at most N (t) independent samples. Define, for s ∈ {−1, +1}, As :=
∇fs (x(r,k) ) .
inf (r,k)∈St
Define Ŝ :=
arg min As′ ,
with ties broken arbitrarily.
(44)
s′ ∈{−1,+1}
If Ŝ ̸= S, then AŜ ≤ AS , and hence 2AS ≥ A1 + A−1 ≥ inf (∥∇f1 (x)∥ + ∥∇f−1 (x)∥) x∈Rd
= L̄ inf (∥x − θ1 ∥ + ∥x − θ−1 ∥) ≥ L̄∥θ1 − θ−1 ∥ = 2rL̄ . x∈Rd
Hence, AS ≥ rL̄ whenever Ŝ ̸= S, and therefore E [AS ] ≥ rL̄ · P(Ŝ ̸= S),
(45)
where P is over the randomness of S, the oracle, and the algorithm. Write Ps for the law of the information available by time t under instance s. Any estimator of S based on this information has error at least 1 P(Ŝ ̸= S) ≥ 1 − ∥P1 − P−1 ∥TV . 2 By Pinsker’s inequality, r
1 DKL (P1 ∥P−1 ) . 2 Since the transcript by time t is a measurable function of at most N (t) i.i.d. samples ξ1 , . . . , ξN (t) with ξj ∼ Pξs under instance s, the data-processing inequality yields DKL (P1 ∥P−1 ) ≤ DKL (Pξ1 )⊗N (t) (Pξ−1 )⊗N (t) = N (t) DKL (Pξ1 ∥Pξ−1 ) . ∥P1 − P−1 ∥TV ≤
Moreover, DKL (Pξ1 ∥Pξ−1 ) = DKL
σ2 σ2 2r2 L̄2 N r, 2 N −r, 2 = . σ2 L̄ L̄
Thus, 1 P(Ŝ ̸= S) ≥ 2
! p rL̄ N (t) 1− , σ
rL̄ E[AS ] ≥ 2
! p rL̄ N (t) 1− . σ
and by (45),
Now set
( r := min
σ p , 2L̄ N (t) 32
r
2∆ L̄
) ,
p p := +∞ when N (t) = 0. Then r ≤ 2∆/L̄ ensures fs (0)−inf fs ≤ ∆, interpreting σ/(2 pL̄ N (t)) and r ≤ σ/(2L̄ N (t)) makes the parenthesis at least 1/2, so ( ) r rL̄ σ L̄∆ E[AS ] ≥ , = min p . 4 8 8 N (t) P Since E[AS ] = 21 s∈{±1} Es [As ], we get ( ) r L̄∆ σ (r,k) max Es inf ∥∇fs (x )∥ ≥ min p , . 8 s∈{±1} (r,k)∈St 8 N (t) Finally, for each s, Jensen’s inequality yields 2 Es inf ∥∇fs (x(r,k) )∥2 = Es A2s ≥ (Es [As ]) , (r,k)∈St
and squaring the previous bound gives (40). The implication (42) follows immediately from ε ≤ L̄∆/8 and N (t) ≤ σ 2 /(64ε). It remains to obtain (43). Set Sε := σ 2 /(64ε). If Sε < 14 , then ε > σ 2 /16 and any algorithm with (r,k) 2 max Es inf ∥∇fs (x )∥ ≤ ε s
(r,k)∈St
must have St ̸= ∅, hence t ≥ τ1 . Moreover, !−1 2 m X 1 1 1 σ σ2 min +m ≤ τ1 + 1 < τ1 , 384 m∈[n] i=1 τi ε 384 ε so (43) holds. Assume now that Sε ≥ 41 . Define τn+1 := ∞ and !−1 m X 1 jε⋆ := inf m ∈ [n] Sε < τm+1 . τ i=1 i Define
⋆ −1 jε X 1 , t1 := Sε τ i i=1
t2 := min
m∈[n]
m X 1 i=1
τi
!−1 (Sε + m) .
By [Tyurin and Richtárik, 2023, Lemma D.7] (applied with S = Sε ), we have t1 ≤ t2 ≤ 6t1 . Moreover, t1 < τjε⋆ +1 , and hence ⌊t1 /τi ⌋ = 0 for all i ≥ jε⋆ + 1, so ⋆
⋆
N (t1 ) =
jε X t1 i=1
τi
≤
jε X t1 i=1
τi
⋆
= t1
jε X 1 i=1
τi
= Sε .
Since N (·) is nondecreasing, for any t ≤ t2 /6 we have t ≤ t1 and thus N (t) ≤ N (t1 ) ≤ Sε . Therefore, any algorithm with maxs Es [inf k∈St ∥∇fs (x(k) )∥2 ] ≤ ε must satisfy t ≥ t2 /6, i.e. !−1 m X 1 1 σ2 t ≥ min +m . 6 m∈[n] i=1 τi 64ε 2 σ2 1 σ Finally, since 64ε + m ≥ 64 + m for all m ∈ [n], we obtain (43). ε
33
Rennala SGD avg(10 seeds) (\# 1): ° = 0:0625, B = 400
Rennala SGD avg(10 seeds) (\# 1): ° = 0:5, B = 200
101
Rennala SGD avg(10 seeds) (\# 2): ° = 0:125, B = 400
100
Rennala MVR avg(10 seeds) (\# 1): ° = 0:5, S = 400, p = 0:001, S0 = 160000, ® = 1:0
10
Rennala MVR avg(10 seeds) (\# 3): ° = 0:125, S = 80, p = 0:001, S0 = 6400, ® = 1:0
0
10
Rennala MVR avg(10 seeds) (\# 2): ° = 0:5, S = 40, p = 0:01, S0 = 1600, ® = 1:0 Rennala MVR avg(10 seeds) (\# 3): ° = 0:5, S = 60, p = 0:005, S0 = 3600, ® = 1:0
10
10−5
10−4
10−4
10−6
10−5
10−5
10−7 0.0
0.2
0.4
0.6
0.8
1.0 1e6
Runtime (seconds)
(a) τi =
Rennala MVR avg(10 seeds) (\# 3): ° = 0:5, S = 80, p = 0:005, S0 = 6400, ® = 1:0
10−3
10−3
10−4
Rennala MVR avg(10 seeds) (\# 2): ° = 0:5, S = 60, p = 0:01, S0 = 3600, ® = 1:0
−1
10−2
10−2
10−3
Rennala SGD avg(10 seeds) (\# 3): ° = 0:25, B = 100 Rennala MVR avg(10 seeds) (\# 1): ° = 0:5, S = 60, p = 0:005, S0 = 3600, ® = 1:0
0
krf(x t )k 2
krf(x t )k 2
krf(x t )k 2
10−1
10−2
Rennala SGD avg(10 seeds) (\# 2): ° = 0:125, B = 40
101
Rennala SGD avg(10 seeds) (\# 3): ° = 0:5, B = 100 Rennala MVR avg(10 seeds) (\# 1): ° = 0:5, S = 40, p = 0:005, S0 = 1600, ® = 1:0
Rennala MVR avg(10 seeds) (\# 2): ° = 0:25, S = 400, p = 0:001, S0 = 160000, ® = 1:0
10−1
Rennala SGD avg(10 seeds) (\# 1): ° = 0:5, B = 200
Rennala SGD avg(10 seeds) (\# 2): ° = 0:25, B = 100
101
Rennala SGD avg(10 seeds) (\# 3): ° = 0:03125, B = 100
0.0
0.2
0.4
0.6
0.8
1.0 1e6
Runtime (seconds)
√ i
0.0
0.2
0.4
0.6
0.8
1.0 1e6
Runtime (seconds)
(b) Mixture delays
(c) Uniform delays
Figure 3: Comparison of Rennala MVR and Rennala SGD on the stochastic quadratic benchmark with 100 workers under three delay models.
F
Additional experiments
We additionally performed experiments on the same stochastic quadratic benchmark with a larger number of workers, namely n = 100. The corresponding plots in Figure 3 show the same qualitative trend as in the n = 10 setting: Rennala MVR consistently outperforms Rennala SGD across all considered delay models. We also examine the sensitivity of Rennala MVR to its additional hyperparameters. For the quadratic benchmark, the only extra parameter is p. The heatmaps in Figure 4 show that, over the tested range, the method is relatively insensitive to the precise choice of p. This is encouraging from a practical perspective, since it suggests that the additional flexibility of Rennala MVR does not translate into a substantial tuning burden. We perform a similar analysis for the neural-network experiments, where the inexact variant introduces both p and the scaling parameter α. We first plot the full heatmaps for fixed values of α; see Figure 5. These plots exhibit a pattern similar to that observed in the quadratic case, namely a relatively weak dependence on p. We therefore additionally minimize over p and plot the resulting heatmaps as functions of the stepsize and α; see Figure 6. The resulting plots indicate that, within the tested range, the method is also relatively insensitive to the choice of α, which further supports its practical viability. Rennala-Mvr hyperparameter sensitivity on quadratic benchmark S=1 best score=4.807e-06
S=40 best score=1.937e-06
S=80 best score=1.504e-06
0.9
S=200 best score=8.693e-07
7
4.0
8 7
0.5
7
6
3.5
6
0.2
6
3.0
5
5
4
3
3
1.5
0.01
2
2
1.0
log 10 (score=best)
3
4
log 10 (score=best)
2.0
0.05
4
log 10 (score=best)
2.5
0.1
log 10 (score=best)
S0 = S p
5
2
0.005 1
0.5
1
1
0.001 0.0 -15
-14
-13
-12
-11
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
1
0
2
-15
-14
-13
-12
-11
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
1
0
2
-15
-14
-13
-12
-11
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
1
0
2
-15
-14
-13
-12
-11
-10
-9
-8
-7
-6
-5
-4
log 2 °
log 2 °
log 2 °
log 2 °
S=1 best score=4.807e-06
S=40 best score=1.157e-06
S=80 best score=9.459e-07
S=200 best score=5.385e-07
-3
-2
-1
0
1
2
8
0.9 4.0
8
7 7
0.5
3.5
7
6 6
0.2
3.0
6 5
4
5
4
3
1.5
log 10 (score=best)
3
log 10 (score=best)
2.0
0.05
4
log 10 (score=best)
0.1
log 10 (score=best)
S0 = S^2 p
5 2.5
3
0.01 2
2
1.0
2
0.005 1
0.5
1
1
0.001 0.0 -15
-14
-13
-12
-11
-10
-9
-8
-7
-6
log 2 °
-5
-4
-3
-2
-1
0
1
2
0 -15
-14
-13
-12
-11
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
1
2
0 -15
log 2 °
-14
-13
-12
-11
-10
-9
-8
-7
-6
log 2 °
-5
-4
-3
-2
-1
0
1
2
0 -15
-14
-13
-12
-11
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
1
2
log 2 °
Figure 4: Sensitivity of exact Rennala MVR on the stochastic quadratic benchmark under square-root √ delays τi = i. Each heatmap shows the performance criterion over the (γ, p) grid for a fixed choice of B and B0 .
34
Rennala-Mvr NN hyperparameter sensitivity | criterion=loss | alpha=0.001 S=1 best score=8.736e-02
S=5 best score=2.853e-02
S=10 best score=1.646e-02
Rennala-Mvr NN hyperparameter sensitivity | criterion=loss | alpha=0.005
S=20 best score=7.069e-03
S=40 best score=1.332e-02
S=60 best score=2.811e-02
S=1 best score=8.983e-02
3.5
3.5
3.5
S=10 best score=1.567e-02
3.5
0.1
S=40 best score=1.306e-02
3.5
S=60 best score=2.672e-02
3.5
3.5
3.5
3.0
3.0
2.5
2.5
1.5
1.5
2.0
1.5
2.0
1.5
2.0
1.5
1.0
1.0
1.0
1.0
0.5
0.5
0.5
0.5
0.5
0.5
0.0
0.0
0.0
0.0
0.0
0.0
0.05
S0 = S p
2.0
1.0
log 10 (score=best)
2.0
1.0
2.0
1.5
0.01
0.005
2.0
1.5
2.0
1.5
2.0
1.5
2.0
1.5
1.0
1.0
1.0
1.0
1.0
0.5
0.5
0.5
0.5
0.5
0.0
0.0
0.0
0.0
0.0
log 10 (score=best)
3.0
2.5
log 10 (score=best)
3.0
2.5
log 10 (score=best)
3.0
2.5
log 10 (score=best)
3.0
2.5
log 10 (score=best)
3.0
2.5
log 10 (score=best)
3.0
2.5
log 10 (score=best)
3.0
2.5
log 10 (score=best)
3.0
2.5
log 10 (score=best)
3.0
2.5
0.001
2.0
1.5
1.0
0.5
0.0
0.001
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
log 2 °
log 2 °
log 2 °
log 2 °
log 2 °
log 2 °
log 2 °
log 2 °
log 2 °
log 2 °
log 2 °
log 2 °
S=1 best score=8.736e-02
S=5 best score=2.471e-02
S=10 best score=1.612e-02
S=20 best score=8.624e-03
S=40 best score=2.045e-02
S=60 best score=5.647e-02
S=1 best score=8.983e-02
S=5 best score=2.426e-02
S=10 best score=1.247e-02
S=20 best score=6.359e-03
S=40 best score=1.804e-02
S=60 best score=5.890e-02
-2
-1
0
1.0
0.005
1.0
1.0
1.0
1.0
0.5
0.5
0.5
0.5
0.5
0.5
0.0
0.0
0.0
0.0
0.0
0.0
3.5
3.5
3.5
3.5
3.0
3.0
3.0
3.0
2.5
2.5
2.5
2.5
2.5
2.0
2.0
1.5
1.0
0.005
0.001
2.0
1.5
1.0
2.0
1.5
1.0
2.0
1.5
1.0
log 10 (score=best)
3.5
3.0
2.5
log 10 (score=best)
3.5
3.0
1.5
0.01
1.0
log 10 (score=best)
0.1
0.05
log 10 (score=best)
0.01
2.0
1.5
S0 = S^2 p
2.5
2.0
1.5
log 10 (score=best)
3.5
3.0
2.5
2.0
1.5
log 10 (score=best)
3.5
3.0
2.5
2.0
1.5
log 10 (score=best)
3.5
3.0
2.5
2.0
1.5
log 10 (score=best)
3.5
3.0
2.5
2.0
1.5
log 10 (score=best)
3.5
3.0
2.5
log 10 (score=best)
3.5
3.0
2.0
1.5
1.0
0.5
0.5
0.5
0.5
0.5
0.0
0.0
0.0
0.0
0.0
log 10 (score=best)
0.2
0.1
0.05
log 10 (score=best)
0.2
1.0
0.5
0.0
0.001
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
log 2 °
log 2 °
S=1 best score=8.403e-02
S=5 best score=2.960e-02
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
log 2 °
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
log 2 °
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
log 2 °
log 2 °
log 2 °
log 2 °
S=40 best score=1.239e-02
S=60 best score=2.608e-02
S=1 best score=8.525e-02
S=5 best score=2.304e-02
-2
-1
0
-10
(a) α = 0.001 S=10 best score=1.487e-02
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
log 2 °
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
log 2 °
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
log 2 °
-5
-4
-3
-2
-1
0
log 2 °
(b) α = 0.005
Rennala-Mvr NN hyperparameter sensitivity | criterion=loss | alpha=0.01
Rennala-Mvr NN hyperparameter sensitivity | criterion=loss | alpha=0.025
S=20 best score=5.548e-03
S=10 best score=1.358e-02
S=20 best score=3.477e-03
S=40 best score=1.262e-02
S=60 best score=2.215e-02
0.2
2.0
1.5
2.0
1.5
1.0
1.0
1.0
0.5
0.5
0.5
0.5
0.5
0.5
0.0
0.0
0.0
0.0
0.0
0.0
0.05
0.001
3.5
3.0
3.0
3.0
3.0
3.0
3.0
2.5
2.5
2.5
2.5
2.5
1.5
0.005
3.5
2.5
2.0
0.01
3.5
2.0
1.5
2.0
1.5
2.0
1.5
2.0
1.5
1.0
1.0
1.0
1.0
1.0
0.5
0.5
0.5
0.5
0.5
0.0
0.0
0.0
0.0
0.0
2.0
1.5
log 10 (score=best)
1.5
1.0
3.5
log 10 (score=best)
2.0
1.0
3.5
log 10 (score=best)
3.0
2.5
1.5
log 10 (score=best)
3.0
2.5
2.0
log 10 (score=best)
3.0
2.5
1.5
log 10 (score=best)
3.0
2.5
2.0
log 10 (score=best)
3.0
2.5
1.0
3.5
0.1
log 10 (score=best)
3.5
log 10 (score=best)
3.5
S0 = S p
3.5
3.0
1.5
0.005
3.5
2.5
2.0
0.01
log 10 (score=best)
S0 = S p
0.05
3.5
log 10 (score=best)
3.5
0.1
log 10 (score=best)
0.2
1.0
0.5
0.0
0.001
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
log 2 °
log 2 °
log 2 °
log 2 °
log 2 °
log 2 °
log 2 °
log 2 °
log 2 °
log 2 °
log 2 °
log 2 °
S=1 best score=8.403e-02
S=5 best score=2.488e-02
S=10 best score=1.516e-02
S=20 best score=4.777e-03
S=40 best score=2.015e-02
S=60 best score=5.444e-02
S=1 best score=8.525e-02
S=5 best score=2.430e-02
S=10 best score=9.902e-03
S=20 best score=4.074e-03
S=40 best score=1.677e-02
S=60 best score=5.023e-02
-2
-1
0
0.005
0.5
0.5
0.5
0.5
0.5
0.5
0.0
0.0
0.0
0.0
0.0
0.0
3.5
3.5
3.5
3.5
3.5
3.0
3.0
3.0
3.0
3.0
2.5
2.5
2.5
2.5
2.5
2.5
2.0
2.0
1.5
1.0
0.005
0.001
2.0
1.5
1.0
2.0
1.5
1.0
2.0
1.5
1.0
log 10 (score=best)
3.5
3.0
1.5
0.01
log 10 (score=best)
0.1
0.05
1.0
log 10 (score=best)
2.0
1.5
1.0
log 10 (score=best)
2.5
2.0
1.5
1.0
S0 = S^2 p
3.5
3.0
2.5
2.0
1.5
1.0
log 10 (score=best)
3.5
3.0
2.5
2.0
1.5
1.0
log 10 (score=best)
3.5
3.0
2.5
2.0
1.5
1.0
log 10 (score=best)
3.5
3.0
2.5
2.0
log 10 (score=best)
3.5
3.0
2.5
log 10 (score=best)
3.5
3.0
1.5
0.01
log 10 (score=best)
0.1
0.05
2.0
1.5
1.0
0.5
0.5
0.5
0.5
0.5
0.0
0.0
0.0
0.0
0.0
log 10 (score=best)
0.2
log 10 (score=best)
0.2
S0 = S^2 p
S=20 best score=1.100e-02
3.5
3.0
1.5
0.005
3.5
2.5
2.0
0.01
3.5
log 10 (score=best)
S0 = S p
0.05
S0 = S^2 p
S=5 best score=2.757e-02
0.2
3.5
0.1
log 10 (score=best)
0.2
1.0
0.5
0.0
0.001
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
log 2 °
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
log 2 °
-5
-4
-3
-2
-1
0
-10
-9
-8
log 2 °
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
log 2 °
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
log 2 °
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
log 2 °
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
log 2 °
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
log 2 °
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
log 2 °
(c) α = 0.01
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
log 2 °
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
log 2 °
-6
-5
-4
-3
-2
-1
0
log 2 °
(d) α = 0.025
Figure 5: Sensitivity √ of inexact Rennala MVR on asynchronous neural-network training under squareroot delays τi = i. Each panel corresponds to a fixed value of α and shows the performance over the (γ, p) grid for different choices of B and B0 .
Rennala-Mvr NN hyperparameter sensitivity | criterion=loss | min over p S=1 best score=3.997e-03, p=0.2
S=5 best score=9.723e-04, p=0.005
S=10 best score=7.193e-04, p=0.01
S=20 best score=6.996e-04, p=0.05
S=40 best score=1.145e-03, p=0.01
S=60 best score=4.835e-03, p=0.05
4.0
4.0
3.5
3.5
3.5
3.0
3.0
3.0
3.0
S0 = S ®
2.0 1.5
2.5 2.0 1.5
2.5 2.0 1.5
2.5 2.0 1.5
2.5 2.0 1.5
log 10 (score=best)
4.0
3.5
3.0
2.5
log 10 (score=best)
4.0
3.5
3.0
log 10 (score=best)
4.0
3.5
log 10 (score=best)
4.0
log 10 (score=best)
0.01
2.5 2.0 1.5
log 10 (score=best)
0.025
0.005 1.0
1.0
1.0
1.0
1.0
1.0
0.5
0.5
0.5
0.5
0.5
0.5
0.0
0.0
0.0
0.0
0.0
0.0
0.001
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
log 2 °
log 2 °
log 2 °
log 2 °
log 2 °
log 2 °
S=1 best score=3.997e-03, p=0.2
S=5 best score=6.603e-04, p=0.005
S=10 best score=4.730e-04, p=0.01
S=20 best score=1.589e-04, p=0.01
S=40 best score=3.559e-04, p=0.01
S=60 best score=1.901e-03, p=0.01
-1
0
4.0
4.0
3.5
3.5
3.5
3.0
3.0
3.0
3.0
S0 = S^2 ®
2.5 2.0 1.5
2.5 2.0 1.5
2.5 2.0 1.5
2.5 2.0 1.5
2.5 2.0 1.5
log 10 (score=best)
4.0
3.5
3.0
log 10 (score=best)
4.0
3.5
3.0
log 10 (score=best)
4.0
3.5
log 10 (score=best)
4.0
log 10 (score=best)
0.01
2.5 2.0 1.5
log 10 (score=best)
0.025
0.005 1.0
1.0
1.0
1.0
1.0
1.0
0.5
0.5
0.5
0.5
0.5
0.5
0.0
0.0
0.0
0.0
0.0
0.0
0.001
-10
-9
-8
-7
-6
-5
log 2 °
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
log 2 °
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
-10
log 2 °
-9
-8
-7
-6
-5
log 2 °
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
log 2 °
-4
-3
-2
-1
0
-10
-9
-8
-7
-6
-5
-4
-3
-2
-1
0
log 2 °
Figure 6: Sensitivity √ of inexact Rennala MVR on asynchronous neural-network training under squareroot delays τi = i after minimizing over p. Each heatmap shows the performance as a function of the stepsize and α for a fixed choice of B and B0 .
35