Zhengyuan Zhou
email: zzhou@stern.nyu.edu

Panayotis Mertikopoulos
email: panayotis.mertikopoulos@imag.fr

Nicholas Bambos
email: bambos@stanford.edu

Yinyu Peter W Glynn
email: glynn@stanford.edu

PANAYOTIS MERTIKOPOULOS ♯ Nicholas Bambos ⋄

AND Peter W Glynn ⋄

Yinyu Ye
email: yinyu-ye@stanford.edu

Distributed stochastic optimization with large delays

Keywords: 2020 Mathematics Subject Classification. Primary 90C15, 90C26, secondary 90C25, 90C06 Distributed optimization, delays, stochastic gradient descent, stochastic approximation

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

With the advent of high-performance computing infrastructures that are capable of handling massive amounts of data, distributed stochastic optimization has become the predominant paradigm in a broad range of applications in operations research [START_REF] Dey | Analysis of sparse cutting planes for sparse milps with applications to stochastic milps[END_REF][START_REF] Feng | Online robust pca via stochastic optimization[END_REF][START_REF] Lee | Stochastic optimization models for rapid detection of viruses in cellphone networks[END_REF][START_REF] Lian | A comprehensive linear speedup analysis for asynchronous stochastic parallel optimization from zeroth-order to first-order[END_REF][START_REF] Ruszczyński | Stochastic programming models[END_REF][START_REF] Shapiro | A tutorial on stochastic programming[END_REF][START_REF] Uryasev | Stochastic Optimization: Algorithms and Applications[END_REF]. Starting with a series of seminal contributions by Tsitsiklis et al. [START_REF] Tsitsiklis | Distributed asynchronous deterministic and stochastic gradient optimization algorithms[END_REF], recent years have witnessed a commensurate surge of interest in the parallelization of first-order methods, ranging from ordinary (stochastic) gradient descent [START_REF] Agarwal | Distributed delayed stochastic optimization[END_REF][START_REF] Chaturapruek | Asynchronous stochastic convex optimization: the noise is in the noise and sgd don't care[END_REF][START_REF] Feyzmahdavian | An asynchronous mini-batch algorithm for regularized stochastic optimization[END_REF][START_REF] Lian | Asynchronous parallel stochastic gradient for nonconvex optimization[END_REF][START_REF] Mania | Perturbed iterate analysis for asynchronous stochastic optimization[END_REF][START_REF] Paine | Gpu asynchronous stochastic gradient descent to speed up neural network training[END_REF][START_REF] Recht | Hogwild: A lock-free approach to parallelizing stochastic gradient descent[END_REF], to coordinate/dual coordinate descent [START_REF] Avron | Revisiting asynchronous linear solvers: Provable convergence rate through randomization[END_REF][START_REF] Fercoq | Accelerated, parallel, and proximal coordinate descent[END_REF][START_REF] Liu | Asynchronous stochastic coordinate descent: Parallelism and convergence properties[END_REF][START_REF] Liu | An asynchronous parallel randomized kaczmarz algorithm[END_REF][START_REF] Tappenden | On the complexity of parallel coordinate descent[END_REF][START_REF] Tran | Scaling up stochastic dual coordinate ascent[END_REF], randomized Kaczmarz algorithms [START_REF] Liu | An asynchronous parallel randomized kaczmarz algorithm[END_REF], online methods [START_REF] Héliou | Gradient-free online learning in continuous games with delayed rewards[END_REF][START_REF] Hsieh | Multi-agent online optimization with delays: Asynchronicity, adaptivity, and optimism[END_REF][START_REF] Joulani | Delay-tolerant online convex optimization: Unified analysis and adaptive-gradient algorithms[END_REF][START_REF] Quanrud | Online learning with adversarial delays[END_REF], block coordinate descent [START_REF] Marecek | Distributed block coordinate descent for minimizing partially separable functions[END_REF][START_REF] Wang | Parallel and distributed block-coordinate frank-wolfe algorithms[END_REF][START_REF] Wright | Coordinate descent algorithms[END_REF], ADMM [START_REF] Hong | A distributed, asynchronous and incremental algorithm for nonconvex optimization: An admm approach[END_REF][START_REF] Zhang | Asynchronous distributed admm for consensus optimization[END_REF], and many others.

This popularity is a direct consequence of Moore's law of silicon integration and the commensurately increased distribution of computing power. For instance, in a typical supercomputer cluster, up to several thousands of "workers" (or sometimes tens of thousands) perform independent computations with little to no synchronization -as the cost of such coordination quickly becomes prohibitive in terms of overhead and energy spillage. Similarly, massively parallel computing grids and data centers (such as those of Google, Amazon, IBM or Microsoft) may house up to several million computing nodes and/or servers, all working asynchronously to execute a variety of different tasks. Finally, taking the concept of distributed computing to its logical extreme, volunteer computing grids (such as Berkeley's BOINC infrastructure or Stanford's folding@home project) essentially span the entire globe and harness the computing power of a vast, heterogeneous network of non-clustered nodes that receive and process computational requests in a non-concurrent fashion, rendering syncrhonization impossible. In this way, by eliminating the required coordination overhead, asynchronous operations become simultaneously more appealing (in physically clustered systems) and more scalable (in massively parallel and/or volunteer computing grids).

In this broad context, perhaps the most widely deployed method is distributed asynchronous stochastic gradient descent (DASGD) and its variants. In addition to its long history in mathematical optimization, DASGD has also emerged as one of the principal algorithmic schemes for training large-scale machine learning models. In "big data" applications in particular, obtaining first-order information on the underlying learning objective is a formidable challenge, to the extent that the only information that can be readily computed is an imperfect, stochastic gradient thereof [START_REF] Dean | Large scale distributed deep networks[END_REF][START_REF] Dean | Large scale distributed deep networks[END_REF][START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF][START_REF] Paine | Gpu asynchronous stochastic gradient descent to speed up neural network training[END_REF][START_REF] Zhang | Asynchronous stochastic gradient descent for dnn training[END_REF][START_REF] Zhang | Deep learning with elastic averaging sgd[END_REF]. This information is typically obtained from a group of computing nodes (or processors) working in parallel, and is then leveraged to provide the basis for a distributed descent step.

Depending on the specific computing architecture, the resulting DASGD scheme varies accordingly. More concretely, there are two types of distributed computing architectures that are common in practice: The first is a cluster-oriented, multi-core, shared memory architecture where different processors independently compute stochastic gradients and update a global model parameter [START_REF] Chaturapruek | Asynchronous stochastic convex optimization: the noise is in the noise and sgd don't care[END_REF][START_REF] Feyzmahdavian | An asynchronous mini-batch algorithm for regularized stochastic optimization[END_REF][START_REF] Lian | Asynchronous parallel stochastic gradient for nonconvex optimization[END_REF]. The second is a "master-worker" architecture used predominantly in computing grids (and, especially, volunteer computing grids): here, each worker node independently -and asynchronously -computes a stochastic gradient of the objective and sends it to the master; the master then updates the model's global parameter and sends out new computation requests [START_REF] Agarwal | Distributed delayed stochastic optimization[END_REF][START_REF] Lian | Asynchronous parallel stochastic gradient for nonconvex optimization[END_REF]. In both cases, DASGD is inherently susceptible to delays, a key impediment that is usually absent in centralized stochastic optimization settings. For instance, in a master-worker system, when a worker sends its gradient update to the master, the master may have already updated the model parameters several times (using updates from other workers), so the received gradient is already stale by the time it is received. In fact, even in the perfectly synchronized setting where all workers have the same speed and send input to the master in an exact round-robin fashion, there is still a constant delay that grows roughly proportionally to the number of workers in the system [START_REF] Agarwal | Distributed delayed stochastic optimization[END_REF].

This situation is exacerbated in volunteer computing grids: here, workers typically volunteer their time and resources following a highly erratic and inconstant update/work schedule, often being turned off and/or being used for different tasks for hours (or even days) on end. In such cases, there is no lower bound on the fraction of resources used by a worker to compute an update at any given time (this is especially true in heterogeneous computing grids such as BOINC and SimGrid), meaning in turn that there is no upper bound on the induced delays. This can also happen in parallel computing environments where many tasks with different priorities are executed at the same time across different machines and, likewise, even in multi-core infrastructures with a shared memory, memory-starved processors can become arbitrarily slow in performing gradient computations.

From a theoretical standpoint, the issue of delays and asynchronicities has been studied from the early days of distributed computing [START_REF] Bertsekas | Parallel and Distributed Computation: Numerical Methods[END_REF], and one of the principal results in the field is the subsequential convergence of DASGD with probability 1 when no constraints are present and when the observed delays grow moderately with time -i.e., sublinearly relative to a global clock [START_REF] Tsitsiklis | Problems in decentralized decision making and computation[END_REF][START_REF] Tsitsiklis | Distributed asynchronous deterministic and stochastic gradient optimization algorithms[END_REF]. In several current systems (for example, in volunteer computing grids), as slower workers become saturated and accumulate computation requests over time while new (and possibly faster) workers enter the system, delays can quickly add up and grow at a superlinear rate relative to the system's global timer. Further, in several applications, there are natural constraints imposed on a subset (or all) of the decision variables that represent the model parameters. In such contexts, the following questions remain open: How robust is the performance envelope of DASGD for constrained optimization under large delays and asynchronicities? Can this robustness be leveraged from a theoretical viewpoint in order to design new and more efficient algorithms?

1.1. Our Contributions and Related Work. Our aim in this paper is to establish the convergence of DASGD in the presence of large, superlinear delays, in as wide a class of objectives as possible and in the presence of constraints where efficient projection can be performed. To that end, we focus on the following classes of problems, where different convergence results can be obtained:

General non-convex objectives. We first consider the class of general smooth non-convex functions, with no structural assumption on the objectives. In this (difficult) case, Tsitsiklis et al. [START_REF] Tsitsiklis | Distributed asynchronous deterministic and stochastic gradient optimization algorithms[END_REF] showed that, under sublinear delays, DASGD converges to a level set of the objective which contains a critical point with probability 1; in particular, if every such point is a global minimizer (e.g., if the problem is pseudo-convex) and the method is run with an Ω(1/n) step-size schedule, DASGD converges to the problem's solution set. More recently, Lian et al. [START_REF] Lian | Asynchronous parallel stochastic gradient for nonconvex optimization[END_REF] derived an estimate for the rate of convergence of the surrogate length n -1 n k=1 E[∥∇ f (X k )∥ 2 2 ] as n → ∞ under the assumption that the delays affecting the algorithm are bounded. Our first contribution is to show that these assumptions on the delays are not needed: specifically, as we show in Theorem 1, by tuning the algorithm's step-size appropriately, it is possible to retain this convergence guarantee, even if delays grow as polynomials of arbitrary degree.

Variationally coherent objectives. Albeit directly applicable to general non-convex stochastic optimization, Theorem 1 only guarantees convergence to stationarity in the mean square sense; to ensure global optimality, stronger structural assumptions on the objectives must be imposed. The "gold standard" of such assumptions (and by far the most widely studied one) is convexity: in the context of distributed stochastic convex optimization, recent works by Agarwal and Duchi [START_REF] Agarwal | Distributed delayed stochastic optimization[END_REF] and Recht et al. [START_REF] Recht | Hogwild: A lock-free approach to parallelizing stochastic gradient descent[END_REF], have established convergence for DASGD under bounded delays for each of the two distributed computing architectures, while Chaturapruek et al. [START_REF] Chaturapruek | Asynchronous stochastic convex optimization: the noise is in the noise and sgd don't care[END_REF] and Feyzmahdavian et al. [START_REF] Feyzmahdavian | An asynchronous mini-batch algorithm for regularized stochastic optimization[END_REF] extended the bounded delays assumption to a setting with finite-mean i.i.d. delays. To go beyond this framework, we focus the class of mean variationally coherent optimization problems [START_REF] Zhou | Stochastic mirror descent in variationally coherent optimization problems[END_REF][START_REF] Zhou | On the convergence of mirror descent beyond stochastic convex programming[END_REF], which includes pseudo-, quasiand/ star-convex problems, as well as many other classes of functions with non-convex profiles. Our main result here is that, in such problems, the global state parameter X n of DASGD converges to a global minimum with probability 1, even when the delays between gradient updates and requests grow at a polynomial rate (and this, without any distributional assumption on how the underlying delays are generated).

To go beyond this framework, we focus on a class of unimodal problems, which we call variationally coherent, and which properly includes all pseudo-, quasi-and/ star-convex problems, as well as many other classes of functions with highly non-convex profiles. Our main result here may be stated as follows: in stochastic variationally coherent problems, provided a lazy descent scheme is used (akin to dual averaging) to mesh with the constraint set in the distributed procedure, the global state parameter X n of DASGD converges to a global minimum with probability 1, even when the delays between gradient updates and requests grow at a polynomial rate (and this, without any distributional assumption on how the underlying delays are generated).

This result extends the works mentioned above in several directions: specifically, it shows that 1. Convexity is not required to obtain global convergence results. 2. Constraints do not hinder almost sure convergence under a suitable lazy projection scheme. 3. The robustness of DASGD is guaranteed even under large, superlinear delays.

We find these outcomes particularly appealing because, coupled with the existing rich literature on the topic, they help explain and reaffirm the prolific empirical success of DASGD in large-scale machine learning problems, and offer concrete design insights for fortifying the algorithm's distributed implementation against delays and asynchronicities.

Techniques. Our analysis relies on techniques and ideas from stochastic approximation and (sub-)martingale convergence theory. A key feature of our approach is that, instead of focusing on the discrete-time algorithm, we first establish the convergence of an underlying, deterministic dynamical system by means of a particular energy (Lyapunov) function which is decreasing along continuous-time trajectories and "quasi-decreasing" along the iterates of DASGD. To control this gap, we connect the continuous-and discrete-time frameworks via the theory of asymptotic pseudotrajectories (APTs), as pioneered by Benaïm and Hirsch [START_REF] Benaïm | Asymptotic pseudotrajectories and chain recurrent flows, with applications[END_REF]. By itself, the APT method does not suffice to establish convergence under delays. However, if the step-size of the method is chosen appropriately (following a quasi-linear decay rate for polynomially growing delays), it is possible to leverage L p martingale tail convergence results to show that the problem's solution set is recurrent under DASGD. This, combined with the above, allows us to prove our core convergence results.

Even though the ordinary differential equation (ODE) approximation of discrete-time Robbins-Monro algorithms has been widely studied in control and optimization theory [START_REF] Ghadimi | Stochastic first-and zeroth-order methods for nonconvex stochastic programming[END_REF][START_REF] Kushner | Stochastic Approximation and Recursive Algorithms and Applications[END_REF], transferring the convergence guarantees of an ODE solution trajectory to a discretetime algorithm is a fairly subtle affair that must be done on a case-by-case basis. Further, even if this transfer is complete, the results typically have the nature of convergence-in-probability: almost-sure convergence is usually much harder to obtain [START_REF] Borkar | Stochastic Approximation: A Dynamical Systems Viewpoint[END_REF]. Specifically, exisiting stochastic approximation results cannot be applied to our setting because a) the non-invertibility of the projection map makes the underlying dynamical system on the problem's feasible region non-autonomous (so APT results do not apply); b) unbounded delays only serve to aggravate this issue, as they introduce a further disconnect between the DASGD algorithm and its continuous-time version. To control the discrepancy between discrete and continuous time requires a more fine-grained analysis, for which we resort to a sharper law of large numbers for L p -bounded martingales. Finally, we also mention that the recent work of Lian et al. [START_REF] Lian | A comprehensive linear speedup analysis for asynchronous stochastic parallel optimization from zeroth-order to first-order[END_REF] has also considered distributed zeroth-order methods (where only the function value, rather than the gradient is available) and used techniques that are different from gradient-based analyses.

Problem Setup

Let X be a subset of R d and let (Ω, F, P) be some underlying (complete) probability space. Throughout this paper, we focus on the following stochastic optimization problem:

minimize f (x) subject to x ∈ X , (Opt)
where the objective function f : X → R is of the form

f (x) = E[F (x; ω)] = Ω F (x; ω) d P(ω) (1)
for some random function 1 F : X × Ω → R. Using standard optimization terminologies, (Opt) is called an unconstrained stochastic optimization problem if X = R d , and is called a constrained stochastic optimization problem otherwise. In this paper, we study both unconstrained and constrained stochastic problems under smooth objectives. Specifically, we make the following regularity assumptions for the rest of the paper, which are standard in the literature:

Assumption 1. We assume the following:

(1) F (x; ω) is differentiable in x for P-almost all ω ∈ Ω.

(2) ∇ F (x; ω) has 2 finite second moment, that is,

sup x∈X E[∥∇F (x; ω)∥ 2 2 ] < ∞. (3) ∇F (x; ω) is Lipschitz continuous in the mean: E[∇F (x; ω)] is Lipschitz on X . Remark 1.
Assumptions 1 and 2 together imply that f is differentiable, because finite second moment (by Statement 2) implies finite first moment: E[∥∇F (x; ω)∥ 2 ] < ∞ for all x ∈ X ; and hence the expectation E[∇F (x; ω)] exists. By a further application of the dominated convergence theorem, we have ∇f

(x) = ∇ E[F (x; ω)] = E[∇F (x; ω)].
Additonally, Assumption 3 implies that ∇f is Lipschitz continuous. In the deterministic optimization literature, f is sometimes called L-smooth, where L is the Lipschitz constant.

One important class of motivating applications that can be cast in the current stochastic optimization problem (Opt) is empirical risk minimization (ERM) in machine learning. As is well-known in the distributed optimization/learning literature [START_REF] Agarwal | Distributed delayed stochastic optimization[END_REF][START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF][START_REF] Lian | Asynchronous parallel stochastic gradient for nonconvex optimization[END_REF][START_REF] Paine | Gpu asynchronous stochastic gradient descent to speed up neural network training[END_REF][START_REF] Zhang | Asynchronous stochastic gradient descent for dnn training[END_REF], the expectation in Eq. ( 1) contains as a special case the common machine learning objectives of the form 1 N N i=1 f i (x), where each f i (x) is the loss associated with the i-th training sample. This setup corresponds to ERM without regularization. With regularization, ERM takes the form 1 N N i=1 f i (x) + r(x), where r(•) is a regularizer (typically convex and known), which is again a special case of (Opt). Other related examples that are also special cases of (Opt) include the objective N i=1 v i f i (x), which are standard in curriculum learning: v i are weights (between 0 and 1) generated from a learned curriculum that indicates how much emphasis each loss f i should be given.

In the large-scale data setting (N is very large), such problems are typically solved in practice using stochastic gradient descent (SGD) on a distributed computing architecture. SGD 3 is widely used primarily because in many applications, stochastic gradient, computed by first drawing a sub-batch of data and then computing the average of the gradients on that sub-batch, is the only type of information that is computationally feasible and practically convenient to obtain 4 . Further, such problems are generally solved on a distributed computing architecture because: 1) Computing gradients is typically the computational bottleneck. Consequently, having multiple processors compute gradients in parallel can harness the available computing power in modern distributed systems. 2) The data (which determine 1 It is understood that a random function here means that F (x; •) : Ω → R is a measureable function for each x ∈ X .

2 It is understood here that the gradient ∇ F (x; ω) is only taken with respect to x: no differential structure is assumed on Ω. 3 In vanilla SGD (i.e. centralized/single-processor setting), an iid sample of the gradient of F at the current iterate is used to make a descent step (with an appropriate projection made in the constrained optimization case).

4

For instance, Google's Tensorflow system does automatic differentiation on samples for neural networks (i.e. f i (•)'s are parametrized neural networks). the individual cost functions f i 's) may simply be too large to fit on a single machine; and hence a distributed system is necessary even from a storage perspective.

With the above background, our goal in this paper is to study and establish theoretical convergence guarantees for applying stochastic gradient descent (SGD) to solve (Opt) on a distributed computing architecture. Two common distributed computing architectures that are widely delpoyed in practice (see also Fig. 1):

(1) Master-worker system. This architecture is mostly used in data-centers and parallel computing grids (each computing node is a single machine, virtual or physical). (2) Multi-processor system with shared memory. This architecture is mostly used in multi-core machines or GPU computing: in the former, each processor is a CPU, while in the latter, each processor is a GPU. In the next two subsections (Sections 2.1 and 2.2), we describe the standard procedure of parallelizing SGD on each of the two distributed computing architectures. Although running SGD on these two architectures have some differences, in Section 2.3 we give a meta algorithmic description, called distributed asynchronous stochastic gradient descent (DASGD) that unifies these two parallelizations on the same footing.

2.1. SGD on Master-Worker Systems. Here we consider the first distributed computing architecture: the master-worker system. The standard way of deploying stochastic gradient descent in such systems -and that which we adopt here -is for the workers to asychronously compute stochastic gradients and then send them to the master, 5 while the master updates the global state of the system and pushes the updated state back to the workers ( [START_REF] Agarwal | Distributed delayed stochastic optimization[END_REF][START_REF] Lian | Asynchronous parallel stochastic gradient for nonconvex optimization[END_REF]). This process is presented in Algorithm 1.

Due to the distributed nature of the master-worker system, a gradient received by the master on any given iteration can be stale. As a simple example, consider a fully coordinated update scheme where each worker sends the computed gradient to and receives the updated iterate from the master following a round-robin schedule. In this case, each worker's gradient is received with a delay exactly equal to K -1 (K is the number of workers in the system), because by the time the master receives worker K's computed gradient, the master has already applied K -1 gradient updates from workers 1 to K -1 (and since the schedule is round-robin, this delay of K -1 is true for any one of the K workers).

However, delays can be much worse since we allow full asynchrony: workers can compute and send (stochastic) gradients to the master without any coordinated schedule. In the 5 As alluded to before, in machine learning applications, this is done by sampling a subset of the training data, computing the gradient for each datapoint and averaging over all datapoints in the sample. Workers:

(a) Receive iterate (b) Compute an i.i.d. stochastic gradient (at the received iterate) (c) Send the computed gradient to master 5: until end asynchronous setting, fast workers (workers that are fast in computing gradients) will cause disproprotionately large delays to gradients produced by slow workers (workers that are slow in computing gradients): when a slower worker has finished computing a gradient, a fast worker may have already computed and communicated many gradients to the master. Since the master updates the global state of the system (the current iterate), one can gain a clearer representation of this scheme by looking at the master's update. This is given in Section 2.3.

SGD on Multi-Processor Systems with Shared Memory.

Here we consider the second distributed computing architecture: multi-processor system with shared memory. In this architecture, all processors can access a global, shared memory, which holds all the data needed for computing a (stochastic) gradient, as well as the current iterate (the global state of the system). The standard way of deploying stochastic gradient descent in such systems ( [START_REF] Chaturapruek | Asynchronous stochastic convex optimization: the noise is in the noise and sgd don't care[END_REF][START_REF] Lian | Asynchronous parallel stochastic gradient for nonconvex optimization[END_REF]) is for each processor to independently and asychronously read the current global iterate, compute a stochastic gradient 6 , and then update the global iterate in the shared memory. This process is given Algorithm 2:

Algorithm 2. Running SGD on a Multi-Processor System with Shared Memory Require: K processors and global (shared) memory.

1: The initial iterate in the global memory. The key difference from Algorithm 1 is that there is no central entity that updates the global state; instead, each processor can both read the global state and update it. Since each processor is performing the operations asynchronously, different processors may be reading the same global iterate at the same time. Further, the delays in this case is again caused by the heterogeneity across different processors: if a processor is slow in computing gradients, then by the time it finishes computing its gradient, the global iterate has been updated by other, faster processors many times over, thereby causing its own gradient stale. Here, we also adopt a common assumption that updating the global iterate is an atomic operation (and hence no two processors will be updating the global iterate at the same time). This is justified 7 because performing gradient update is a simple arithemtic operation, and hence takes negligible time compared to reading data and computing a stochastic gradient, which is the main computational bottleneck in practice. However, despite a cheap computation, performing the whole gradient update (typically achieved via locking) does have overhead. In particular, a less stringent model would be only updating one coordinate at a time: this is known as an inconsistent write/read model since different processors are updating different components of the global parameters and hence can read in elements of different ages. For simplicity, we do not consider this case, as our focus in this paper is on delays. See [START_REF] Lian | Asynchronous parallel stochastic gradient for nonconvex optimization[END_REF] and [START_REF] Liu | An asynchronous parallel randomized kaczmarz algorithm[END_REF] for lucid discussions and analyses on this model. Finally, one can gain a clearer picture of this update scheme by tracking the update to the global iterate in the shared memory. This is given in Section 2.3.

DASGD:

A Unifying Algorithmic Representation. In this subsection, we present a unified algorithmic description, aptly called distributed asynchronous stochastic gradient descent (DASGD), that formally captures both Algorithm 1 and Algorithm 2, where their differences are reflected in the assumptions of the meta algorithm's parameters. We start with the unconstrained case, see Algorithm 3.

X n+1 = X n -α n+1 ∇F (X s(n) , ω n+1 );
X n = pr X (Y n ); 4: Y n+1 = Y n -α n+1 ∇F (X s(n) , ω n+1 ); 5: n ← n + 1; 6: until end 7: return solution candidate X n
In more detail, n is a global counter and is incremented every time an update occurs to the current solution candidate X n (the global iterate): in the master-worker systems, the master updates it; in the multi-processor systems, each processor updates it. Since there are delays in both systems, the gradient applied to the current iterate X n can be a gradient associated with a previous time step. This fact is abstractly captured by Line 3 in Algorithm 4. In full 7 On a related note, we also note that our analysis can be further extended to cases where only one variable or a small block of variables are being updated at a time. We omit this discussion because the resulting notation is quite onerous, and will obscure the main ideas behind the already complex theoretical framework developed here.

generality, we will write s(n) for the iteration from which the gradient received at time n originated. In other words, the delay associated with iteration s(n) is n -s(n), since it took n -s(n) iterations for the gradient computed on iteration s(n) to be received at stage n. Note that s(n) is always no larger than n; and if n = s(n), then there is no delay in iteration n.

Now, the difference between the two distributed computing archecitures is reflected in the assumption of s(n). Specifically, in the master-worker systems, each s(•) is a one-toone function 8 , because no two workers will ever receive same iterates from the master per Algorithm 1. On the other hand, in multi-processor systems, s(n) can be the same for different n's (since different processors may read the current iterate at the same time); however, it is easy to observe that the same s will appear at most K times for different n's, since there are K processors in total. As an important note, our analysis is agnostic to whether s(n) is one-to-one or not. Consequently, in establishing theoretical guarantees for the meta algorithm DASGD, we obtain the same guarantees for both architectures simultaneously.

Notation-wise, we will write d n for the delay required to compute a gradient requested at iteration n. This gradient is received at iteration n + d n . Following this notation, the delay for a gradient received at n is d s(n) = n -s(n). Note also we have chosen the subscript associated with ω to be n + 1: we can do so because ω n 's are iid (and hence the indexing is irrelevant). Finally, in constrained optimization case (where X is a strict subset of R d ), projection must be performed. This results in DASGD with projection 9 , which is formally given in Algorithm 4.

General Nonconvex Objectives

In this section, we take X = R d (i.e. unconstrained optimization) and consider general non-convex objectives. Note that for a general non-convex objective where no further structural assumption (e.g. convexity) is made, convergence to an optimal solution (even a local optimal solution) cannot be expected, and will not hold in general, even in the absence of both noise and delays (i.e. single machine deterministic optimization). In such cases, the best one can hope for, which is also the standard metric to determine the stability of the algorithm, is that the gradient vanishes in the limit 10 .

3.1. Delay Assumption. Our goal here is to establish convergence guarantees (in mean square) of DASGD for a general non-convex objective in the presence of delays. In fact, a large family of unbounded delay processes can be tolerated. We state our main assumption regarding the delays and step-sizes:

Assumption 2.
The gradient delay process d n and the step-size sequence α n of DASGD (Algorithms 3 and 4) satisfy one of the following conditions:

(1) Bounded delays: sup n d n ≤ D for some positive number D and

∞ n=1 α 2 n < ∞, ∞ n=1 α n = ∞.
(2) Sublinearly growing delays: d n = O(n p ) for some 0 < p < 1 and α n ∝ 1/n. 8 Except initially if all the workers have the same initial point. 9 This tpe of projection is technically known as lazy projection. 10 An alternative phrase that is commonly used is that the criticality gap vanishes. This is also colloquially referred to as convergence to a stationary point/critical point in the machine learning community. Note further that convergence to second-order-stationary points can be achieved by stochastic gradient descent (and its variants) under weaker assumptions (than convexity): Lipschitz Hessian and strict saddle point property. See [START_REF] Ge | Escaping from saddle points?online stochastic gradient for tensor decomposition[END_REF][START_REF] Jin | How to escape saddle points efficiently[END_REF][START_REF] Lee | Gradient descent converges to minimizers[END_REF] for this line of work.

(3) Linearly growing delays:

d n = O(n) and α n ∝ 1/(n log n).
(4) Polynomially growing delays: d n = O(n q ) for some q ≥ 1 and α n ∝ 1/(n log n log log n).

Note that as delays get larger, we need to use less aggressive step-sizes. This is to be expected, because the larger the delays, the more "averaging" one needs perform in order to remove the staleness that is caused by the delays; and smaller step-sizes correspond to averaging over a longer horizon. This is a one of the important insights that occur throughout the paper. Another thing to note that is the Assumption 2 also highlights the quantitative relationship between the class of delays and the class of step-sizes. For instance, when the delays increase from a linear rate to a polynomial rate, only a factor of 1 log log n needs to be added (which is effectively a constant). From a practical standpoint, this means that a step-size on the order of 1/(n log n) will be a good model-agnostic choice and more-or-less sufficient for almost all delay processes.

3.2.

Controlling the Tail Behavior of Second Moments. We now turn to establish the theoretical convergence guarantees of DASGD for general non-convex objectives. Our first step lies in controlling the tail behavior of the second moments of the gradients that are generated from DASGD. By leveraging the Lipschitz continuity of the gradient, its telescoping sum, appropriate conditionings and a careful analysis of the interplay between delays and step-sizes, we show that (next lemma) a particularly weighted tail sum of the second moments are vanishingly small in the limit (see appendix for a detailed proof).

Lemma 1.

Under Assumptions 1 and 2, if

inf x∈X f (x) > -∞, then ∞ n=0 α n+1 E[∥ ∇ f (X n )∥ 2 2 ] < ∞. (2)
Remark 2. Since [START_REF] Bertsekas | Neuro-dynamic programming[END_REF]). Note that the converse is not true: when a subsequence of

∞ n=0 α n+1 = ∞, Lemma 1 implies that lim inf n→∞ E[∥ ∇ f (X n )∥ 2 2 ] = 0 (see Appendix A of
E[∥ ∇ f (X n )∥ 2
2 ] converges to 0, the sum in Equation ( 5) need not be finite. As a simple example, consider α n+1 = 1 n , and

E[∥ ∇ f (X n )∥ 2 2 ] = 1 n , if n = 2 k 1, otherwise. (3)
Then the subsequence on indicies 2 n converges to 0, but the sum still diverges. Consequently, Lemma 1 is stronger than subsequence convergence.

3.3.

Bounding the Successive Differences. However, Lemma 1 is still not strong enough to guarantee that

lim n→∞ E[∥ ∇ f (X n )∥ 2 2 ] = 0.
This is because the convergent sum given in Equation ( 5) only limits the tail growth somewhat, but not completely. To demonstrate this point, let c t be the following boolean variable:

c n = 1, if n contains the digit 9 in its decimal expansion, 0, otherwise. (4)
For instance, c 9 = 1, c 11 = 0. Now define α n+1 = 1 n , and

E[∥ ∇ f (X n )∥ 2 2 ] = 1 n , if c n = 1 1, if c n = 0.
As first-year Berkeley Math PhDs delightfully found out during-or painfully found out after -their qualifying exam,

∞ n=1 α n+1 E[∥ ∇ f (X n )∥ 2 2 ] < ∞ (see Problem 1.3.24 in [12]), even though the limit E[∥ ∇ f (X n )∥ 2 2 ] does not exist. This indicates that to obtain convergence of E[∥ ∇ f (X n )∥ 2
2 ], we need to impose more stringent conditions to ensure its sufficient tail decay. Note that one issue that is revealed by the above counter-example is that the distance between two successive terms is always bounded away from 0, no matter how large n is. This obviously makes it possible for convergence to occur: a necessary condition for convergence is that the difference converges to 0. Note that intuitively, this pathological case cannot occur for gradient descent, because the step-size is shrinking to 0, hence making the successive difference converge to 0 (at least in expectation). Consequently, we can bound the difference between every two successive terms in terms of a decreasing sequence that is converging to 0. This ensures that E[∥ ∇ f (X n )∥ 2 2 ] cannot change two much from iteration to iteration. Further, the change between two successive terms will be vanishing. This result is formalized in the following lemma (the proof is given in the appendix):

Lemma 2.
Under Assumptions 1 and 2, there exists a constant C > 0 such that for every n,

E[∥ ∇ f (X n+1 )∥ 2 2 ] -E[∥ ∇ f (X n )∥ 2 2 ] ≤ Cα n+1 .
3.4. Main Convergence Result. Putting the above two characterizations together, we obtain:

Theorem 1.
Let X n be the DASGD iterates generated from Algorithm 3. Under Assumptions 1 and 2, if

inf x∈X f (x) > -∞, then lim n→∞ E[∥ ∇ f (X n )∥ 2 2 ] = 0. (5)
Remark 3. Three remarks are in order here. First, note that the condition inf x∈X f (x) > -∞ means that the optimization problem (Opt) has a solution. This is necessary, because otherwise, a stationary point may not exist in the first place, and DASGD (or even simple gradient descent) may continue decrease the objective value ad infinitum. Note that since f is smooth, a minimum point is necessarily a stationary point. Second, the above convergence is a fairly strong characterization of the fact that the gradient vanishes. In particular, it means that the gradient of the DASGD iterates converges to 0 in mean square. Consequently, this implies that the norm of the gradient vanishes in expectation, and that the gradient converges to 0 with high probability. Note further that if we strengthen Assumption 1 to require that all stochastic gradients are bounded almost surely (as opoosed to just bounded in second moments), then a similar analysis ensures almost sure convergence of ∥ ∇ f (X n )∥ 2 . We omit the details due to space limitation. Finally, that Theorem 1 is a direct consequence of Lemmas 1 and 2 is a simple excercise in elementary series theory (in particular, Lemma A.4).

Variationally Coherent Problems

In this section, our goal is to establish global optimality convergence guarantees of DASGD in as wide a class of optimization problems as possible. Since global convergence cannot be expected to hold for all non-convex optimization problems (even without delays). a standard structural assumption to make in the existing literature on the objectives (even in the no-delay case) is convexity. Here we instead consider a much broader class of stochastic optimization problems than convex problems. Further, we allow for constrained optimization; in particular, X is assumed to be a convex and compact subset of R d throughout the section. We first discuss the class of objectives and then present global convergence results and their analyses in the subsequent two subsections. We focus on the class of mean variationally coherent optimization problems [START_REF] Mertikopoulos | Optimistic mirror descent in saddle-point problems: Going the extra (gradient) mile[END_REF][START_REF] Mertikopoulos | Learning in games with continuous action sets and unknown payoff functions[END_REF][START_REF] Zhou | On the convergence of mirror descent beyond stochastic convex programming[END_REF][START_REF] Zhou | Distributed asynchronous optimization with unbounded delays: How slow can you go[END_REF], defined here as follows:

Assumption 3.
The optimization problem (Opt) is called variationally coherent in the mean if E[⟨∇ F (x; ω), x -x * ⟩] ≥ 0, (VC) for all x * ∈ X * and all x ∈ X with equality only if x ∈ X * . Note that we do not impose the "if" condition for equality: if x ∈ X * , then the equality may or may not hold. By Assumption 1, we can interchange expectation and differentiation in (VC) to obtain ⟨∇f (x), x -x * ⟩ ≥ 0, (6) for all x ∈ X , x * ∈ X * . As a result, mean variational coherence can be interpreted as an averaged coherence condition for the deterministic optimization problem with objective f (x). Mean variationally coherent optimization problems include convex programs, pseudoconvex programs, non-degenerate quasi-convex programs and star-convex programs as special cases. For instance, if g is star-convex, then ⟨∇f (x), x -x * ⟩ ≥ f (x) -f (x * ) for all x ∈ X , x * ∈ X * . This is easily seen to be a special case of variational coherence because ⟨∇f (x), x -x * ⟩ ≥ f (x) -f (x * ) ≥ 0, with the last inequality strict unless x ∈ X * . Note that star-convex functions contain convex functions as a subclass (but not necessarily pseudo/quasi-convex functions). See [START_REF] Zhou | Stochastic mirror descent in variationally coherent optimization problems[END_REF][START_REF] Zhou | On the convergence of mirror descent beyond stochastic convex programming[END_REF] for a more detailed discussion on why these various classes of optimization problems are special cases. Fig. 2 also provides two more elaborate examples of a variationally coherent optimization problem that are not quasi-convex. Put together, these examples indicate that variationally coherent objectives can have highly non-convex profiles. Nevertheless, for the class of variationally coherent functions, it is possible to establish almost sure convergence guarantees (which we do so in the subsequent sections). 4.1. Deterministic Analysis: Convergence to Global Optima. To streamline our presentation and build intuition along the way, we will begin with the deterministic case in this subsection, where there is no randomness in the calculation of a gradient update. In this case, DASGD boils down to distributed asynchronous gradient descent (DAGD), as illustrated in Algorithm 5: x n = pr X (y n );

4:

y n+1 = y n -α n+1 ∇f (x s(n) ); 5:
n ← n + 1; 6: until end 7: return solution candidate x n 4.1.1. Energy Function. We start with an energy function that measures on how "optimal" the dual variable y is: the smaller the energy (i.e. the closer it is to 0), the better the dual variable.

Definition 1.

Let x * ∈ X * . Define the energy function E : R d → R as follows:

E(y) = inf x * ∈X * E x * (y), where E x * (y) = ∥x * ∥ 2 2 -∥pr X (y)∥ 2 2 + 2⟨y, pr X (y) -x * ⟩. ( 7
)
Note that one can think of E x * (y) as the energy of y with respect to a fixed optimal solution x * , while E(y) is the best (smallest) energy for a given y. We next characterize a few of its useful properties.

Lemma 3.

For all y ∈ R d , we have:

(1) E(y) ≥ 0 with equality if and only if pr X (y) ∈ X * .

(

) Let {y n } ∞ n=1 be a given sequence. If lim n→∞ E(y n ) = 0, then pr X (y n ) → X * as 11 n → ∞. 2
Remark 4. The proof is given in the appendix, but it is helpful to make a few quick remarks. The first statement justifies the terminology of "energy", as E(y) is always non-negative. This energy function will also be the tool we use to establish an important component of the global convergence result. Further, given that E(y) ≥ 0, it should also be clear that when pr X (y) ∈ X * ⊂ X , we can choose a particular x * = pr X (y) such that E(y) = 0 (but that E(y) = 0 implies pr X (y) ∈ X * is less obvious). Statement 2 of the lemma provides us with a way to establish convergence to optimal solutions. If we can show that E(y n ) → 0, then x n = pr X (y n ) → X * . Nevertheless, as we shall see later, unlike many other Lyapunov 11 Following the convention in point-set topology, a sequence sn converges to a set S if dist(sn, S) → 0, with dist(•, •) being the standard point-to-set distance function: dist(sn, S) ≜ inf s∈S dist(sn, s), where dist(sn, s) is the Euclidean distance between sn and s.

functions in optimization, E(y n ) does not decrease monotonically; consequently, it is difficult to directly establish E(y n ) → 0. In fact, a finer-grained analysis is required to characterize the convergence behavior of E(y n ) (see Section 4.1.2).

We also assume that the converse of Statement 2 of the above lemma holds:

Assumption 4. If lim n→∞ E(y n ) = 0, then pr X (y n ) → X * as n → ∞.
Assumption 4 can be seen as a primal-dual analogue of the reciprocity conditions for the Bregman divergence. This assumption usually holds (e.g. when the feasible set X is a polytope), unless X is pathological.

Lemma 4.

Fix any x * ∈ X * .

(

) ∥ pr X (y) -ŷ∥ 2 2 -∥ pr X (ŷ) -ŷ∥ 2 2 ≤ ∥y -ŷ∥ 2 2 , for any y, ŷ ∈ R d . (2) E x * (y + ∆y) -E x * (y) ≤ 2⟨∆y, pr X (y) -x * ⟩ + ∥∆y∥ 2 1
2 , for any y, ∆y ∈ R d . Remark 5. The first statement of the above lemma serves as an important intermediate step in proving the second statement, and is established by leveraging the envelop theorem and several important properties of Euclidean projection. To see that this is not trivial, consider the quantity ∥ pr X (y) -ŷ∥ 2 -∥ pr X (ŷ) -ŷ∥ 2 , which we know by triangle's inequality satisfies:

∥ pr X (y) -ŷ∥ 2 -∥ pr X (ŷ) -ŷ∥ 2 ≤ ∥ pr X (y) -pr X (ŷ)∥ 2 ≤ ∥y -ŷ∥ 2 , (8)
where the last inequality follows from the fact that projection is a non-expansive map. However, this inequality is not sufficient for our purposes because in quantifying the perturbation E(y + ∆y) -E(y), we also need the squared term ∥∆y∥ 2 2 , which is not easily obtainable from Equation [START_REF] Bertsekas | Nonlinear Programming[END_REF]. In fact, a tighter analysis is needed to establish that ∥y -ŷ∥ 2 2 is an upper bound on ∥ pr X (y) -ŷ∥ 2 2 -∥ pr X (ŷ) -ŷ∥ 2 2 . 4.1.2. Main Convergence Result. An intermediate result, interesting on its own and useful also as a heavy-lifting tool for our convergence analysis is then provided by the following technical result:

Proposition 1.
Under Assumptions 1 to 4, DAGD admits a subsequence x n k that converges to X * as k → ∞.

We highlight the main steps below and refer the reader to the appendix for the details:

(1

) Letting b n = ∇ f (x s(n) ) -∇ f (x n ),
we can rewrite the gradient update in DAGD as:

y n+1 = y n -α n+1 ∇ f (x s(n) ) = y n -α n+1 ∇ f (x n ) -α n+1 {∇ f (x s(n) ) -∇ f (x n )} = y n -α n+1 (∇ f (x n ) + b n ). (9)
Recall here that s(n) denotes the previous iteration count whose gradient becomes available only at the current iteration n. By bounding the magnitude of b n using the delay sequence through a careful analysis, we establish that under any one of the conditions in Assumption 2, lim n→∞ ∥b n ∥ 2 = 0. The analysis here, particularly the one for the last three conditions, reveals the following pattern: as the magnitude of the delays gets larger and larger in the order of growth, one needs to use a more conservative step-size sequence in order to mitigate the damage done by the stale gradient information. Intuitively, smaller step-sizes are more helpful in larger delays because they carry a better "amortization" effect that makes DAGD more tolerant to delays.

(2) With the defintion of b n , DAGD can be written as:

x n = pr X (y n ),

y n+1 = y n -α n+1 (∇ f (x n ) + b n ). ( 10
)
We then use the energy function to study the behavior of y n and x n . More specifically, we look at the quantity E(y n+1 ) -E(y n ) and, using Lemma 4, we bound this onestep change using the step size α n , the b n sequence and the defining quantity ⟨∇ f (x n ), x n -x * ⟩ of a variationally coherent function (as well as another term that will prove inconsequential). We then telescope on E(y n+1 ) -E(y n ) to obtain an upper bound for E(y n+1 ) -E(y 0 ). Since the energy function is always non-negative, E(y n+1 ) -E(y 0 ) is at least -E(y 0 ) for every n. Then, utilizing the fact that b n converges to 0 and that ⟨∇ f (x n ), x n -x * ⟩ is always positive (unless the iterate is exactly an optimal solution, in which case it is 0), we show that the upper bound will approach -∞ if X n only enters N (X * , ϵ), an open ϵ-neighborhood of X * , a finite number of times (for an arbitrary ϵ > 0). This generates an immediate contradiction, and thereby establishes that X n will get arbitrarily close to X * for an infinite number of times. This then implies that there exists a subsequence of DAGD iterates that converges to the solution set of (Opt), i.e, x n k → X * as k → ∞.

Theorem 2. Under Assumptions 1 to 4, the global state variable x n of DAGD (Algorithm 5) converges to the solution set X * of (Opt).

We give an outline of the proof below, referring to the appendix for the details. Fix a δ > 0. Since x n k → x * , as k → ∞ (per Proposition 1), we have E(y n k ) → 0 as k → ∞ per Lemma 3. So we can pick an n that is sufficiently large and E(y n ) < δ. Our goal here is to show that for n large enough, once E(y n ) < δ, it will stay so forever: E(y m ) < δ, ∀m ≥ n. Note in particular that the above statement is not true for all n, but only for n large enough.

However, the behavior of E(y n ) is not very regular: it can certainly increase from iteration to iteration for any n. Nevertheless, we can precisely quantify how large this increment (if any) can be. This leads us to break it down to two cases:

(1) Case 1: E(y n ) < δ/2.

(2) Case 2: δ/2 ≤ E(y n ) < δ. For Case 1, we show in the appendix that

E(y n+1 ) -E(y n ) ≤ 2BC 4 α n+1 + 2α 2 n+1 (C 2 + B 2 ), (11)
for suitable constants B and C's. Now for n sufficiently large, we can make the right-hand arbitrarily small, and in particular, smaller than δ/2. This means E(y n+1 ) ≤ E(y n ) + δ 2 < δ. Consequently, in this case, the energy stays within the δ bound in the next iteration.

For Case 2, we show in the appendix that

E(y n+1 ) -E(y n ) ≤ -2α n+1 a 2 -α n+1 (C 2 + B 2 ) , ( 12
)
where a is a positive constant that depends only on δ. Again, since n is sufficiently large, we can make a 2 -α n+1 (C 2 + B 2 ) positive, thereby making the right-hand side negative. Consequently, E(y n+1 ) < E(y n ) < δ. Hence, again, the energy stays within the δ bound in the next iteration.

The key conclusion from the above is that, for large enough n, once E(y n ) is less than δ, E(y n+1 ) is less than δ as well and so are all the iterates afterwards. Since δ is arbitrary, it follows E(y n ) → 0, and therefore x n → X * by Lemma 3.

Stochastic Analysis:

Almost Sure Convergence to Global Optima. Having established deterministic global convergence of DAGD, we now proceed to study stochastic global convergence of DASGD. Compared to the deterministic analysis, the stochastic case is much more involved because randomness can lead to very volatile behavior in the presence of delays; in particular, the simple approach employed in Theorem 2 (to establish that once E(y n ) is less than δ, it will always remain so) no longer works. To deal with both delays and noise, a much more sophisticated analysis framework needs to be developed, which requires several news ideas. To streamline our presentation, we break the theoretical development into four subsections, each comprising an important component and step of the overall analysis. 4.2.1. Recurrence of DASGD. Our first step lies in generalizing Proposition 1 to the stochastic case. Specifically, in the presence of noise, we show that the iterates of DASGD visit any neighorhood of X * infinitely often almost surely.

Proposition 2.

Under Assumptions 1 to 4, DASGD admits a subsequence X n k that converges to X * almost surely: X n k → X * with probability 1 as k → ∞.

We outline the two main steps of the proof below, referring the reader to the appendix for the details.

(1) We begin by rewriting the gradient update step in DASGD as:

Y n+1 -Y n α n+1 = -∇F (X s(n) , ω n+1 ) = -∇ f (X n ) -[∇ f (X s(n) ) -∇ f (X n )] -[∇F (X s(n) , ω n+1 ) -∇ f (X s(n) )]. (13)
Letting

B n = ∇ f (X s(n) ) -∇ f (X n ) and U n+1 = ∇F (X s(n) , ω n+1 ) -∇ f (X s(n) )
, we can rewrite the DASGD update as

Y n+1 = Y n -α n+1 {∇ f (X n ) + B n + U n+1 }. ( 14
)
We then establish the following two facts in this step. First, we verify that

n r=0 U n+1 is a martingale adapted to Y 1 , Y 2 . . . , Y n+1 , where {U n+1 } ∞ n=0 is a L 2 - bounded martingale difference sequence. Second, we show that lim n→∞ ∥B n ∥ 2 = 0, a.s..
The second claim is done by first giving an upper bound on ∥B n ∥ 2 by writing

∇ f (X s(n) ) -∇ f (X n ) as a sum of one-step changes (∇ f (X s(n) ) -∇ f (X s(n)+1 ) + ∇ f (X s(n)+1 ) -• • • + ∇ f (X n-1 ) -∇ f (X n ))
and analyzing each such successive change. We then break that upper bound into two parts, one deterministic and one stochastic. For the deterministic part, the same analysis in the proof of Proposition 1 yields convergence to 0.

The stochastic part turns out to be the tail of a martingale. By leveraging the property of the step-size and a crucial property of martingale differences (two martingale differences at different time steps are uncorrelated), we establish that said martingale is L 2 -bounded. Then, by applying a version of Doob's martingale convergence theorem, it follows that said martingale converges almost surely to a limit random variable with finite second moment (and hence almost surely finite). Consequently, writing the tail as a difference between two terms (each of which converges to the same limit variablewith probability 1), we conclude that the tail converges to 0 (a.s.).

(2) The full DASGD update may then be written as

X n = pr X (Y n ) Y n+1 = Y n -α n+1 [∇ f (X n ) + B n + U n+1 ]. (15)
As in Step 2 of the proof of Proposition 1, we again bound the one-step change of the energy function E(Y n+1 ) -E(Y n ) and then telescope the differences. The two distinctions from the determinstic case are: 1) Everything is now a random variable. 2) We have three terms: in addition to the random gradient ∇ f (X n ) and the random drift B n , we also have a martingale term U n+1 . Since B n converges to 0 almost surely (as shown in the previous step), its effect can be shown to be negligible. Futher, the analysis utilizes law of large numbers for martingale as well as Doob's martingale convergence theorem to bound the effect of the various martingale terms and to establish that the final dominating term converges to -∞ almost surely (which generates a contradiction since the energy function is always positive) unless a subsequence X n k converges almost surely to X * . ■ Even though recurrence, which can be seen as the counterpart of Proposition 1, holds as per the above proposition, the random iterates X n are much more irregular than their deterministic counterpart x n in DAGD. To deal with this complexity, we work with and characterize the sample trajectories "generated" (to be made precise later) by X n (rather than individual iterates X n ). To work towards this general direction, we first push the DASGD update into a determinstic ordineary differential equation (ODE), as explained in the next subsection. 4.2.2. Mean-Field Approximation of DASGD. We can rewrite the DASGD update as:

X n = pr X (Y n ) Y n+1 = Y n -α n+1 {∇ f (X n ) + B n + U n+1 }. ( 16
)
Written in this way, DASGD can be viewed as a discretization of the "mean-field" ODE

x = pr X (y), ẏ = -∇ f (x). ( 17
)
The intuition is that this ODE provides a "mean" approximation of the DASGD update, because in [START_REF] Feng | Online robust pca via stochastic optimization[END_REF], the noise term U n+1 has 0 mean, and the term B n converges to 0 (and therefore has negilible effect in the long run). Thes leaves only the term Y n+1 = Y nα n+1 ∇ f (X n ), which can be seen as a Euler discretization of the ODE. (Of course, that Equation ( 17) is a good-enough approximation of Equation ( 16) for global almost sure convergence purposes here will be rigorously justified later.) Next, writing Equation ( 17) solely in terms of y yields ẏ = -∇ f (pr X (y)). Since ∇ f and pr X are both Lipschitz continuous and X is a compact set, the composition ∇ f • pr X is itself Lipschitz continuous and bounded. Standard results from the theory of dynamical systems then show that (17) admits a unique global solution y(t) for any initial condition y(0). On the other hand, since pr X is not a one-to-one map, it is not invertible; consequently, there need not exist a unique solution trajectory for x(t). By this token, the rest of our analysis will focus on the trajetory of y(t).

With the guarantee of the existence and uniqueness of the y trajectory, let P : R + × R d → R d be the semiflow 12 of (17), i.e., P (t, y 0 ) denotes the state of (17) at time t when the initial condition is y 0 . In other words, when viewed as a function of time, P (•, y 0 ) is the solution trajectory to ẏ = -∇ f (pr X (y)). It is worth pointing out that writing it in this double-argument form also allows us to interpret P as a function of the initial condition: for a fixed t, P (t, •) gives different states at t when the ODE starts from different initial conditions (in particular, P (0, y) = y). Both views will be useful later.

It turns out that with the energy function introduced here, E(P (t, y)) is always nonincreasing. Furthermore, it is also decreasing at a meaingful rate. We end this subsection with a "sufficient decrease" property of the mean dynamics (17) (the proof given in the appendix due to space limitation):

Lemma 5.
With notation as above, we have:

(1) If pr X (P (t, y)) / ∈ X * , then E(P (t, y)) is strictly decreasing in t for all y ∈ R d .

(2) For all δ > 0, there exists some T ≡ T (δ) > 0 such that, for all t ≥ T , we have

sup y {E(P (t, y)) -E(y) : E(P (t, y)) ≥ δ/2} ≤ -δ/2. ( 18
)
Lemma 5 essentially says E(P (t, y)) is strictly and uniformly (across all y) decreasing at a non-vanishing rate. More specifically, to give some intuition of the second part of Lemma 5, note that E(P (t, y)) -E(y) is the energy change at time t when starting at y. The first part says this quantity is always non-negative. While the second part says provided E(P (t, y)) ≥ δ/2, the decrease in energy will be at least δ 2 , no matter what the initial point y is. If E(P (t, y)) ≥ δ/2 does not hold, that means the energy at time t is already really small (i.e. E(P (t, y)) < δ/2). Consequently, the mantra of the above lemma can be stated succinctly as follows: either the energy is already close to 0, or the energy will decrease towards 0.

In fact, by some additional analysis, one can further show 13 that Lemma 5 implies P (t, y) → X * , ∀y as t → ∞. Now, despite being an interesting result on its own (which establishes that the continuous dynamics of DASGD converges to X * ), it is still some distance away from our final desideratum: our goal is to establish (almost sure) convergence of the discrete-time process in DASGD. So unless we can somehow relate the discrete-time iterates to the continuous-time trajectories, the convergence of DASGD is still uncertain. We fulfill this taks in the next subsection.

4.2.3.

Relating DASGD Iterates to ODE Trajectories. Our goal here is to establish a quantitative connection between the DASGD iterates and the ODE trajectory studied in the previous subsection. Our general idea is that if we show the trajectory generated by the discrete-time iterates of DASGD is path-by-path "close" to the continuous-time trajectory P (t, y), then likely almost-sure convergence of the DASGD iterates can be guaranteed as well.

To be more specific, there are two things that need to be more precisely defined from the preceding high-level discussion. First, what does it mean to be a trajectory generated by the discrete-time iterates of DASGD? Second, what does it mean for this trajectory to be "close" to the ODE trajectory? In general, the answers to these questions can vary depend on the specific goal at hand. In the current context, our goal is to establish global almost sure convergence (a very strong result). Consequently, we need to choose the answers rather judiciously: on the one hand, the answers must be stringent enough to ensure global almost sure convergence in the end (for instance, for the second question, a fairly strong notion of "closeness" is needed); on the other hand, they must also be flexible enough to fit in the current context. 13 Although this is an interesting conclusion, we do not prove it here because we are mainly concerned with establishing convergence of the DASGD iterates, rather than the ODE solution trajectory.

As it turns out, the answer to the first question is rather intuitive: (perhaps) the simplest way to generate a continuous trajectory from a sequence of discrete points is the affine interpolation: connect the iterates Y 0 , Y 1 , . . . , Y n at times 0, α 1 , . . . , n=1 r=1 α r . We call this curve the affine interpolation curve of DASGD and denote it by A(t). Note that A(t) is a random curve because the DASGD iterates Y 0 , Y 1 , . . . , Y n are random. To avoid confusion, we summarize the three different objects discussed so far:

(1) The DASGD iterates Y 0 , Y 1 , . . . , Y n .

(2) The affine interpolation curve A(t) of Y n .

(3) The flow P (t, y) of the ODE ( 17). The answer to the second question lies in the notion of an asymptotic pseudotrajectory (APT) , a concept introduced by Benaïm and Hirsch [START_REF] Benaïm | Asymptotic pseudotrajectories and chain recurrent flows, with applications[END_REF] and Benaïm and Schreiber [START_REF] Benaïm | Ergodic properties of weak asymptotic pseudotrajectories for semiflows[END_REF]. Specifically, in our current context, a continuous curve s(t) is considered close to ODE solution P (t, y) if the following holds:

Definition 2.
A continuous function s : R + → R d is an APT for P if for every T > 0,

lim t→∞ sup 0≤h≤T d(s(t + h), P (h, s(t))) = 0, (19)
where

d(•, •) is the Euclidean metric 14 in R d .
Intuitively, the definition matches exactly the naming: s is an APT for P if, for sufficiently large t, the flow lines of P remain arbitrarily close to s(t) over a time window of any (fixed) length. More precisely, for each fixed T > 0, one can find a large enough t 0 , such that for all t > t 0 , the curve s(t + h) approximates the trajectory P (h, s(t)) on the interval h ∈ [0, T ] with any predetermined degree of accuracy.

With this definition in place, to push through the agenda, we need to establish that A(t), the affine interpolation curve of the DASGD iterates, is an APT for the flow P (t, y) induced by the ODE [START_REF] Fercoq | Accelerated, parallel, and proximal coordinate descent[END_REF]. More precisely, we establish that A(t) is an APT for the flow P (t, y) almost surely, because as mentioned before, A(t) is a random curve.

Lemma 6.

Let A(t) be the random affine interpolation curve generated from the DASGD iterates. Then A(t) is an APT of P (t, y) almost surely.

Note that this result means any resulting affine interpolation curve of DASGD is close to the ODE trajectory. This also forms the basis for reasoning convergence on a path-by-path scale. However, more work still remains to be done because, unfortunately, the condition that A(t) is an APT for P almost surely does not guarantee that the discrete-time iterates of DASGD converge to X * (for many counterexamples in general dynamical systems, see Benaïm [START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF]). In other words, the notion of APT is not sharp enough to ensure direct convergence result. We fulfill this final gap in the next subsection.

Main Convergence Result.

Even though A(t) being an APT for P almost surely does not itself guarantee that the discrete-time iterates of DASGD converge to X * , we can use the energy function to further sharpen this result. Specifically, we use E(A(t)) to further control the behavior of the affine interpolation curve. In fact, the advantage of working with the affine interpolation curve A(t) is that once we show E(A(t)) is bounded by some δ almost surely from some point onwards, then we know E(Y n ) is bounded by δ almost surely (also from some point onwards): this is because the discrete-time iterates and the affine curve 14 As should be obvious from the definition, APTs can be defined more generally in metric spaces in exactly the same way. coincide at discrete time points. Consequently, we focus on boudning E(A(t)), which will enable us to obtain our main convergence result:

Theorem 3.
Under Assumptions 1 to 4, the global state variable X n of DASGD (Algorithm 4) converges (a.s.) to the solution set X * of (Opt).

Again, we only give an outline of the proof below. By Proposition 2, X n gets arbitrarily close to X * infinitely often. Thus, it suffices to show that, if X n ever gets ϵ-close to X * , all the ensuing iterates are ϵ-close to X * (a.s.). The way we show this "trapping" property is to use the energy function. Specifically, we consider E(A(t)) and show that no matter how small ϵ is, for all sufficiently large t 0 , if E(A(t 0 )) is less than ϵ for some t 0 , then E(A(t)) < ϵ, ∀t > t 0 . This would then complete the proof because A(t) actually contains all the DASGD iterates, and hence if E(A(t)) < ϵ, ∀t > t 0 , then E(Y n ) < ϵ for all sufficiently large n. Furthermore, since A(t) contains all the iterates, the hypothesis that " if E(A(t 0 )) is less than ϵ for some t 0 " will be satisfied due to Proposition 2.

We expand on one more layer of detail and defer the rest into appendix. To obtain control E(A(t)), we control two things: the energy on the ODE path E(P (t, y)) and the discrepancy between E(P (t, y)) and E(A(t)). The former can be made arbitrarily small as a result of Lemma 5 (we have a direct handle on how the ODE path would behave). The latter can also be made arbitrarily small as a result of Lemma 6: since A(t) is an APT for P , the two paths are close. Therefore, the discrepancy between E(P ) and E(A) should also be vanishingly small. Consequently, since E(A(t)) = E(P (t, y)) + {E(A(t)) -E(P (t, y))}, and both terms on the right can be made arbitrarily small, so can E(A(t)) be made arbitrarily small.

Discussion

We end the paper with a short simulation discussion that reveals an interesting practical observation. Specifically, we test the convergence of Algorithm 4 against a Rosenbrock test function with d = 1001 degrees of freedom, a standard non-convex global optimization benchmark [START_REF] Rosenbrock | An automatic method for finding the greatest or least value of a function[END_REF]. Specifically, we consider the objective

f Ros (x) = 1000 i=1 [1000(x i+1 -x 2 i ) 2 + (1 -x i ) 2 ], (20)
with x i ∈ [0, 2], i = 1, . . . , 1001. The global minimum of f Ros is located at (1, . . . , 1), at the end of a very thin and very flat parabolic valley which is notoriously difficult for first-order methods to traverse [START_REF] Rosenbrock | An automatic method for finding the greatest or least value of a function[END_REF]. Since the minimum of the Rosenbrock function is known, (VC) is easily checked over the problem's feasible region.

For our numerical experiments, we considered a) a synchronous update schedule as a baseline; and b) an asynchronous master-worker framework with random delays that scale as d n = Θ(n). In both cases, Algorithm 4 was run with a decreasing step-size of the form α n ∝ 1/(n log n) and stochastic gradients drawn from a standard multivariate Gaussian distribution (i.e., zero mean and identity covariance matrix).

Our results are shown in Fig. 3. Starting from a random (but otherwise fixed) initial condition, we ran S = 10 5 realizations of DASGD (with and without delays). We then plotted a randomly chosen trajectory ("test sample" in Fig. 3), the sample average, and the min/max over all samples at every update epoch. For comparison purposes, we also plotted the value of the so-called "ergodic average"

Xn = n k=1 α k X k n k=1 α k , ( 21
) • • • • • • • • • • • • • • • • • •••••••• •• •• • • •••••• • ••• ••• •••• •• ••••• ••••••••••••• •••••••• ••••• •••• ••• ••• •• •• •• • • • • • • • • • • • • • • •• △ △ △ △ △ △ △ △ △ △ △ △ △△△ △△△△△ △△△△△△△ △ △△△△△△△△△△ △△△△△△△△△△△△△△△△△△△ △△△△△△△△△△△△△△△△△△△△△△△ △△△△△△△△△△△△ △△△△△△△△ △△△△△ △△△△ △△△ △ ���� ������ • ������ ������� △ ������� ������� ������ ���/��� � � �� �� ��� ��� ���� �� -� ����� ����� ����� � ��������� �������� �����

����� ����������� (�������� ���� �=��� �������)

(a) Convergence with no delays between gradient updates which is often used in the analysis of DASGD in the convex case (see e.g., 1). Even though this averaging leads to very robust convergence rate estimates in the convex case, we see here that it performs worse than the worst realization of DASGD. The reason for this is the lack of convexity: due to the ridges and talwegs of the Rosenbrock function, Jensen's inequality fails dramatically to produce an improvement over X n (and, in fact, causes delays as it causes X n to deviate from its gradient path). Consequently, this simple simulation indicates that establishing convergence of the iterate X n itself is not only theoretically stronger (and hence more difficult) than convergence of the ergodic average, but also more practically relevant.

• • • • • • • • • • • • • • •• •• •• •• • • •• • • • • • • • • • • • • •• • • • • • • • • • •• • • •• • •• • •• •• •• ••• ••• •••• ••••••• ••••••••• ••••••• •••• ••• ••• •• •• •• •• • • • • • △ △ △ △ △ △ △ △ △ △ △△ △△△ △△△△ △ △△△ △△ △△ △ △ △ △△ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △

A. Auxiliary Results

We collect here in one place all the auxiliary results in the existing literature that will be used in our proofs subsequently. The first one is a well-known characterization of convex sets and the projection operator given in [START_REF] Nesterov | Introductory Lectures on Convex Optimization: A Basic Course[END_REF]:

Lemma A.1.
Let X be a compact and convex subset of R d . Then for any x ∈ X , y ∈ R d :

⟨pr X (y) -x, pr X (y) -y⟩ ≤ 0. (A.1)
The second one is an L p -bounded martingale convergence theorem given in [START_REF] Hall | Martingale Limit Theory and Its Application[END_REF]:

Lemma A.2.
Let S n be a martingale adapted to the filtration S n . If for some p ≥ 1,

sup n≥0 E[|S n | p ] < ∞, then S n converges almost surely to a limiting random variable S ∞ with E[|S ∞ | p ] < ∞. Remark 6. Note that E[|S ∞ | p ] < ∞ for p ≥ 1 obviously implies S ∞ is finite almost surely.
The third one is the classical envelope theorem (see [START_REF] Carter | Foundations of Mathematical Economics[END_REF]).

∇ V (θ) = ∇ f (x * (θ), θ). (A.3)
The fourth one is an elementary sequence result (see [START_REF] Bertsekas | Nonlinear Programming[END_REF]).

Lemma A.4. Let a n , b n be two non-negative sequences such that

∞ n=1 a n = ∞, ∞ n=1 a n b n < ∞. If there exists a real number K > 0 such that |b n+1 -b n | ≤ Ka n . Then, lim n→∞ b n = 0.
The fifth one is law of large numbers for martingales given in [START_REF] Hall | Martingale Limit Theory and Its Application[END_REF]:

Lemma A.5. Let M n = n k=0 d k be a martingale adapted to (F n ) ∞ n=0 and let (u n ) ∞ n=0 be a nondecreasing sequence of positive numbers with lim n→∞ u n = ∞. If ∞ n=0 u -p n E[|d k | p |F n ] < ∞ for some p ∈ [1, 2] (a.s.), then: lim n→∞ M n u n = 0 (a.s.) (A.4)
Finally, we recall the standard notion of semiflow. Definition 3. A semiflow P on a metric space (M, d) is a continuous map P : R + ×M → M :

(t, x) → P t (x),
such that the semi-group properties hold: P 0 = identity, P t+s = P t •P s for all (t, s) ∈ R + ×R + .

Remark 7. A standard way to induce a semiflow is via an ODE. Specifically, if F : R m → R m is a continuous function and if the following ODE has a unique solution trajectory for each initial point x ∈ R m :

dx dt = F (x),
x(0) = x, then P t (x) defined by the solution trajectory x(t) ∈ R m as follows is a semiflow: P t (x) ≜ x(t) with x(0) = x. We say P defined in this way is the semiflow induced by the corresponding ODE.

B. Technical Proofs

B.1. General Nonconvex Objectives.

B.1.1. Proof of Lemma 1. Proof: We start by rewriting the delayed gradient update X n+1 = X n -α n+1 ∇F (X s(n) , ω n+1 ) in Algorithm 3 in two forms as follows, both of which will be used subsequently:

X n+1 = X n -α n+1 ∇ f (X s(n) ) + ∇ F (X s(n) , ω n+1 ) -∇ f (X s(n) ) = X n -α n+1 ∇ f (X n ) + ∇ f (X s(n) ) -∇ f (X n ) + ∇ F (X s(n) , ω n+1 ) -∇ f (X s(n) ) . (B.1) Denoted B ≜ sup x∈X E[∥∇F (x, ω)∥ 2 2 ] per Assumption 1. Since ∇ f (x) = E[∇ F (x, ω)
] is Lipschitz per Assumption 1, letting L be the Lipschitz constant, we have:

f (X n+1 ) -f (X n ) ≤ ⟨∇ f (X n ), X n+1 -X n ⟩ + L 2 ∥X n+1 -X n ∥ 2 2 = -α n+1 ⟨∇ f (X n ), ∇F (X s(n) , ω n+1 )⟩ + L 2 ∥α n+1 ∇F (X s(n) , ω n+1 )∥ 2 2 = -α n+1 ⟨∇ f (X n ), ∇ f (X n ) + ∇ f (X s(n) ) -∇ f (X n ) + ∇ F (X s(n) , ω n+1 ) -∇ f (X s(n) )⟩ + L 2 α 2 n+1 ∥ ∇ f (X s(n) ) + ∇ F (X s(n) , ω n+1 ) -∇ f (X s(n) )∥ 2 2 = -α n+1 ∥ ∇ f (X n )∥ 2 2 -α n+1 ⟨∇ f (X n ), ∇ f (X s(n) ) -∇ f (X n )⟩ -α n+1 ⟨∇ f (X n ), ∇ F (X s(n) , ω n+1 ) -∇ f (X s(n) )⟩ + L 2 α 2 n+1 ∥ ∇ f (X s(n) ) + ∇ F (X s(n) , ω n+1 ) -∇ f (X s(n) )∥ 2 2 ≤ -α n+1 ∥ ∇ f (X n )∥ 2 2 + α n+1 ∥ ∇ f (X n )∥ 2 ∥ ∇ f (X s(n) ) -∇ f (X n )∥ 2 -α n+1 ⟨∇ f (X n ), ∇ F (X s(n) , ω n+1 ) -∇ f (X s(n) )⟩ + L 2 α 2 n+1 2∥ ∇ f (X s(n) )∥ 2 2 + 2∥ ∇ F (X s(n) , ω n+1 ) -∇ f (X s(n) )∥ 2 2 ≤ -α n+1 ∥ ∇ f (X n )∥ 2 2 + √ Bα n+1 ∥ ∇ f (X s(n) ) -∇ f (X n )∥ 2 -α n+1 ⟨∇ f (X n ), ∇ F (X s(n) , ω n+1 ) -∇ f (X s(n) )⟩ ≤ Lα 2 n+1 √ B + ∥ ∇ F (X s(n) , ω n+1 ) -∇ f (X s(n) )∥ 2 2 ,
where in the last inequality follows because by Jensen's inequality, we have:

∥ ∇ f (X n )∥ 2 2 ≤ sup x∈X ∥ E[∇F (x, ω)]∥ 2 2 ≤ sup x∈X E[∥∇F (x, ω)∥ 2 2 ] ≤ B.
Denote the filtration generated by X 0 , X 1 , . . . , X n to be F n . We take the expectation of both sides of Equation (B.2) and obtain:

E[f (X n+1 ) -f (X n )] ≤ -α n+1 E[∥ ∇ f (X n )∥ 2 2 ] + √ Bα n+1 E[∥ ∇ f (X s(n) ) -∇ f (X n )∥ 2 ] -α n+1 E[⟨∇ f (X n ), ∇ F (X s(n) , ω n+1 ) -∇ f (X s(n) )⟩] + L √ Bα 2 n+1 + Lα 2 n+1 E[∥ ∇ F (X s(n) , ω n+1 ) -∇ f (X s(n) )∥ 2 2 ] = -α n+1 E[∥ ∇ f (X n )∥ 2 2 ] + √ Bα n+1 E[∥ ∇ f (X s(n) ) -∇ f (X n )∥ 2 ] -α n+1 E E ⟨∇ f (X n ), ∇ F (X s(n) , ω n+1 ) -∇ f (X s(n) )⟩ F n + L √ Bα 2 n+1 + Lα 2 n+1 E[∥ ∇ F (X s(n) , ω n+1 ) -∇ f (X s(n) )∥ 2 2 ] = -α n+1 E[∥ ∇ f (X n )∥ 2 2 ] + √ Bα n+1 E[∥ ∇ f (X s(n) ) -∇ f (X n )∥ 2 ] -α n+1 E ⟨∇ f (X n ), E ∇ F (X s(n) , ω n+1 ) -∇ f (X s(n) ) F n ⟩ + L √ Bα 2 n+1 + Lα 2 n+1 E[∥ ∇ F (X s(n) , ω n+1 ) -∇ f (X s(n) )∥ 2 2 ] = -α n+1 E[∥ ∇ f (X n )∥ 2 2 ] + √ Bα n+1 E[∥ ∇ f (X s(n) ) -∇ f (X n )∥ 2 ] + L √ Bα 2 n+1 + Lα 2 n+1 E[∥ ∇ F (X s(n) , ω n+1 ) -∇ f (X s(n) )∥ 2 2 ] ≤ -α n+1 E[∥ ∇ f (X n )∥ 2 2 ] + √ Bα n+1 E[∥ ∇ f (X s(n) ) -∇ f (X n )∥ 2 ] + L √ Bα 2 n+1 + 4LBα 2 n+1 ,
where the third equality follows from

E ∇ F (X s(n) , ω n+1 ) -∇ f (X s(n) ) F n = 0, since ω n+1 is independent of F n , and the last inequality follows from E[∥ ∇ F (X s(n) , ω n+1 ) - ∇ f (X s(n) )∥ 2 2 ] ≤ sup x∈X E[∥ ∇ F (x, ω n+1 ) -∇ f (x)∥ 2 2 ] ≤ 4B. Since ∇ f is L-Liptchiz contin- uous, we have: ∥ ∇ f (X s(n) ) -∇ f (X n )∥ 2 ≤ L∥X s(n) -X n ∥ 2 = L X s(n) -X s(n)+1 + X s(n)+1 -X s(n)+2 + • • • + X n-1 -X n 2 = L n-1 r=s(n) X r -X r+1 2 = L n-1 r=s(n) α r+1 ∇F (X s(r) , ω r+1 ) 2 ≤ L n-1 r=s(n) α r+1 ∇F (X s(r) , ω r+1 ) 2 . (B.2)
Taking the expectation of both sides of the above equation then yields:

E[∥ ∇ f (X s(n) ) -∇ f (X n )∥ 2 ] ≤ L n-1 r=s(n) α r+1 E[∥∇F (X s(r) , ω r+1 )∥ 2 ] ≤ L n-1 r=s(n) α r+1 sup x∈X E[∥∇F (x, ω r+1 )∥ 2 ] (B.3) ≤ L B n-1 r=s(n) α r+1 , (B.4)
where

B ≜ sup x∈X E[∥∇F (x, ω r+1 )∥ 2 ] < ∞ per Remark 1. Combining Equation (B.3) with Equation (B.
2) then yields:

E[f (X n+1 )] -E[f (X n )] ≤ -α n+1 E[∥ ∇ f (X n )∥ 2 2 ] + L B√ Bα n+1 n-1 r=s(n) α r+1 + L √ Bα 2 n+1 + 4LBα 2 n+1 . (B.5)
Telescoping Equation (B.5) then yields:

-∞ < inf x∈X f (x) -f (X 0 ) ≤ E[f (X T +1 )] -f (X 0 ) = T n=0 E[f (X n+1 )] -E[f (X n )] ≤ - T n=0 α n+1 E[∥ ∇ f (X n )∥ 2 2 ] + L B√ B T n=0 α n+1 n-1 r=s(n) α r+1 + (L √ B + 4LB) T n=0 α 2 n+1 . (B.6)
Taking T → ∞, the above equation yields:

-∞ < - ∞ n=0 α n+1 E[∥ ∇ f (X n )∥ 2 2 ] + L B√ B ∞ n=0 α n+1 n-1 r=s(n) α r+1 + (L √ B + 4LB) ∞ n=0 α 2 n+1 . (B.7)
Note that in all cases in Assumption 2, we have

∞ n=0 α 2 n+1 < ∞. We next proceed to bound ∞ n=0 α n+1 n-1
r=s(n) α r+1 and show that it is finite in each of the cases in Assumption 2.

(1) In the bounded delays case, since sup n d n ≤ D, it follows that s(n) + D ≥ n and hence:

∞ n=0 α n+1 n-1 r=s(n) α r+1 ≤ ∞ n=0 α n+1 n r=n-1-D α r+1 ≤ (D + 1) ∞ n=0 α n+1 max r∈{n-1-D,...,n} α r+1 ≤ (D + 1) ∞ n=0 max r∈{n-1-D,...,n} α r+1 • max r∈{n-1-D,...,n} α r+1 = (D + 1) ∞ n=0 max r∈{n-1-D,...,n} α 2 r+1 ≤ (D + 1) ∞ n=0 n r=n-1-D α 2 r+1 = (D + 1) 2 ∞ n=0 α 2 n+1 < ∞,
where all the terms α r are defined to be 0 when r drops below 0. (2) In the sublinearly growing delays case, since d n = O(n p ), it follows that s(n) + Ks p (n) ≥ n for some positive number K, which further implies that s(n

) + Ks(n) ≥ s(n) + Ks p (n) ≥ n, thereby leading to s(n) ≥ n K+1
. Consequently, we have:

∞ n=0 α n+1 n-1 r=s(n) α r+1 ≤ ∞ n=0 1 n n r=s(n) 1 r ≤ ∞ n=0 1 n Ks p (n) s(n) = K ∞ n=0 1 n s p-1 (n) ≤ K ∞ n=0 1 n ( n K + 1 ) p-1 ≤ (K + 1) -p ∞ n=0 n p-2 < ∞. (B .8)
(3) In the linearly growing delays case, since d n = O(n), it follows that s(n) + Ks(n) ≥ n for some positive number K and hence again s(n) ≥ n K+1 . Consequently, we have:

∞ n=0 α n+1 n-1 r=s(n) α r+1 ≤ ∞ n=0 1 n log n n r=s(n) 1 r log r ≤ ∞ n=0 1 n log n s(n)+Ks(n) r=s(n) 1 r log r ≤ ∞ n=0 1 n log n s(n)+Ks(n) s(n) 1 r log r dr = ∞ n=0 1 n log n log log(s(n) + Ks(n)) log s(n) = ∞ n=0 1 n log n log log(K + 1) + log s(n) log s(n) = ∞ n=0 1 n log n log(1 + log(K + 1) log s(n) ) ≤ ∞ n=0 1 n log n log(K + 1) log s(n) ≤ ∞ n=0 1 n log n log(K + 1) log n -log(K + 1) ∼ K ∞ n=0 1 n(log n) 2 < ∞.
(B.9) (4) In the polynomially growing delays case, since d n = O(n q ), it follows that s(n) + Ks q (n) ≥ n for some positive number K. Note that in this case, s(n) = Ω(n

1 2q ), because otherwise, s(n) + Ks q (n) = O(n 1 2q + Kn 1 2 ) = o(n), which is a contradiction. Consequently, we have: ∞ n=0 α n+1 n-1 r=s(n) α r+1 ≤ ∞ n=0 1 n log n log log n n r=s(n) 1 r log r log log r ≤ ∞ n=0 1 n log n log log n s(n)+Ks q (n) r=s(n) 1 r log r log log r ≤ ∞ n=0 1 n log n log log n s(n)+Ks q (n) s(n) 1 r log r log log r dr ≤ ∞ n=0 1 n log n log log n log log((K + 1) log s(n) + a log s(n)) log log s(n) ≤ ∞ n=0 1 n log n log log n log(1 + log(K + 1 + a) log log s(n) ) ≤ ∞ n=0 1 n log n log log n log(K + 1 + a) log log s(n) ≤ ∞ n=0 1 n log n log log n log(K + 1 + a) log log n 1 2q K ′ ∼ K ∞ n=0 1 n log n(log log n) 2 < ∞. (B.10)
Consequently, in each of the above 4 cases, we have

∞ n=0 α n+1 n-1 r=s(n) α r+1 < ∞. Therefore, Equation (B.7) yields -∞ < - ∞ n=0 α n+1 E[∥ ∇ f (X n )∥ 2 2 ] + L B√ B ∞ n=0 α n+1 n-1 r=s(n) α r+1 + (L √ B + 4LB) ∞ n=0 α 2 n+1 ≤ - ∞ n=0 α n+1 E[∥ ∇ f (X n )∥ 2 2 ] + C, (B.11)
for some finite constant C. Reversing the above inequality immediately yields the result:

∞ n=0 α n+1 E[∥ ∇ f (X n )∥ 2 2 ] < ∞. ■ B.1.2. Proof of Lemma 2. Proof:
We first recall a useful fact: for any two vectors a, b and any finite-dimensional vector norm ∥ • ∥,

(∥a∥ + ∥b∥)(∥a∥ -∥b∥) ≤ ∥a + b∥∥a -b∥. (B.12)
Using this fact, we can expect bound the term in question as follows:

E[∥∇f (X n+1 )∥ 2 2 ] -E[∥∇f (X n )∥ 2 2 ] = E ∥∇f (X n+1 )∥ 2 + ∥∇f (X n )∥ 2 ∥∇f (X n+1 )∥ 2 -∥∇f (X n )∥ 2 ≤ E ∥∇f (X n+1 )∥ 2 + ∥∇f (X n )∥ 2 • ∥∇f (X n+1 )∥ 2 -∥∇f (X n )∥ 2 ≤ E ∇f (X n+1 ) + ∇f (X n ) 2 • ∇f (X n+1 ) -∇f (X n ) 2 ≤ E 2 sup x∈X ∥∇f (x)∥ 2 • ∇f (X n+1 ) -∇f (X n ) ≤ 2 √ BE ∇f (X n+1 ) -∇f (X n ) 2 ≤ 2 √ BE L X n+1 -X n 2 = 2L √ BE X n+1 -X n 2 = 2L √ BE α n+1 ∇F (X s(n) , ω n+1 ) 2 ≤ 2L √ Bα n+1 sup x∈X E ∇F (x, ω n+1 ) 2 = 2L √ Bα n+1 sup x∈X E ∇F (x, ω n+1 ) 2 2 ≤ 2L √ Bα n+1 sup x∈X E ∇F (x, ω n+1 ) 2 2 ≤ 2L √ Bα n+1 √ B = 2LBα n+1 ,
where the first inequality is an application of Jensen's inequality, the second inequality follows from Equation (B.12) the fifth inequality follows from Lipschitz continuity and the second-to-last inequality follows from another application of Jensen's inequality and that the squre root function is concave. ■ B.2. Variationally Coherent Objectives. We define the following constants that will be handy later:

(1)

C 1 = sup x∈X E[∥∇F (x; ω)∥ 2 ]. (2)
C 2 = sup x∈X E[∥∇F (x; ω)∥ 2 2 ]. (3) C 3 is the Lipschitz constant for ∇ f (x)(= E[∇F (x; ω)]) : ∥ ∇ f (x) -∇ f (x ′ )∥ 2 ≤ C 3 ∥x -x ′ ∥ 2 . (B.13) (4) C 4 = sup x,x ′ ∈X ∥x -x ′ ∥ 2 .
B.2.1. Proof of Lemma 3. Per the definition of the energy function, we have:

E x * (y) -∥ pr X (y) -x * ∥ 2 2 = ∥x * ∥ 2 2 -∥ pr X (y)∥ 2 2 + 2⟨y, pr X (y) -x * ⟩ -∥ pr X (y)∥ 2 2 -2⟨pr X (y), x * ⟩ + ∥x * ∥ 2 2 = -2∥ pr X (y)∥ 2 2 + 2⟨y -pr X (y), pr X (y) -x * ⟩ + 2⟨pr X (y), pr X (y) -x * ⟩ + 2⟨pr X (y), x * ⟩ = 2⟨y -pr X (y), pr X (y) -x * ⟩ ≥ 0, (B.14)
where the last inequality follows from Lemma A.1. Consequently, E x * (y) ≥ ∥ pr X (y)-x * ∥ 2 2 ≥ 0. Taking the infimum over X * yields E(y) ≥ 0.

Next, we establish that E(y) = 0 if and only if pr X (y) ∈ X * . The if part is already established in Remark 4. It suffices to show that E(y) = 0 implies pr X (y) ∈ X * . To see this, observe that if E x * (y) = 0, we must have ∥ pr X (y) -x * ∥ 2 2 = 0, therefore implying pr X (y) = x * . Since E x * (y) is a continuous function of x * for each fixed y, and since X * is a compact set, inf x * ∈X * E x * (y) must be achieved by some z * ∈ X * . Namely E(y) = E z * (y). Consequently, by the preceding observation, E(y) = 0 implies pr X (y) = z * ∈ X * .

For the second statement, suppose on the contrary E(y n ) → 0 but pr X (y n ) does not converge to X * . Then there must exist a subsequence n k such that pr X (y n k ) is bounded away from X * . Denoting by N (X * , ϵ) the ϵ-open ball around X * (i.e. N (X * , ϵ) ≜ {x ∈ R d | dist(x, X * ) < ϵ}). Then for the subsequence n k , there must exist some positive ϵ such that pr X (y n ) remains in X ∩ N c (X * , ϵ). Since X ∩ N c (X * , ϵ) is an intersection of two closed sets, it is itself closed; it is also bounded since X is bounded: hence it is compact. Further, since for each fixed x * , E x * (y) is a continuous function of pr X (y), E x * (y) must achieve the minimum value on the compact set X ∩ N c (X * , ϵ), where the minimum value a x * is positive per the first statement:

E x * (y) ≥ a x * > 0, ∀ pr X (y) ∈ X ∩ N c (X * , ϵ), ∀x * ∈ X *
Finally, since E x * (y) is continuous in x * , it must achieve the minimum value (over x * ) on the compact set X * , where the minimum value a (corresponding to some E x * (y)) must again be positive:

E(y) = inf x * ∈X * E x * (y) ≥ inf x * ∈X * a x * = a > 0. Consequently, E(y n k ) ≥ a > 0 since pr X (y n k ) ∈ X ∩ N c (X * ,
∥ pr X (y) -ŷ∥ 2 2 -∥ pr X (ŷ) -ŷ∥ 2 2 = ∥ pr X (y) -y + y -ŷ∥ 2 2 -∥ pr X (ŷ) -ŷ∥ 2 2 = ∥y -ŷ∥ 2 2 + ∥ pr X (y) -y∥ 2 2 + 2⟨pr X (y) -y, y -ŷ⟩ -∥ pr X (ŷ) -ŷ∥ 2 2 = ∥y -ŷ∥ 2 2 -∥ pr X (ŷ) -ŷ∥ 2 2 -∥ pr X (y) -y∥ 2 2 -2⟨y -pr X (y), ŷ -y⟩ . (B.15)
Now define the function f (x, y) = ∥x -y∥ 2 2 . It follows easily that the solution to the problem max x∈X ∥x -y∥ 2 2 is x * (y) = pr X (y). Consequently, by Lemma A.3, V (y) = f (x * (y), y) is a differential function in y and its derivative can be computed explicitly as follows:

∇ V (y) = ∇ f (x * (y), y) = 2(y -x * (y)) = 2(y -pr X (y)). (B.16)
Futher, since for each x ∈ X , f (x, y) is a convex function in y, and taking the maximum preserves convexity, we have V (y) is also a convex function in y. This means that V (ŷ) -V (y) -⟨∇ V (y), ŷ -y⟩ ≥ 0.

By Equation (B. [START_REF] Feng | Online robust pca via stochastic optimization[END_REF]) and that V (y) = f (x * (y), y) = ∥ pr X (y) -x∥ 2 2 , the above equation becomes:

∥ pr X (ŷ) -ŷ∥ 2 2 -∥ pr X (y) -y∥ 2 2 -2⟨y -pr X (y), ŷ -y⟩ ≥ 0. Consequently, Equation (A.2) then immediately yields: ∥ pr X (y) -ŷ∥ 2 2 -∥ pr X (ŷ) -ŷ∥ 2 2 ≤ ∥y -ŷ∥ 2 2 . We now prove the second part. Expanding using the definition of the Lyapunov function (and skipping some tedious algebra in between), we have:

E x * (y + ∆) -E x * (y) = ∥x * ∥ 2 2 -∥ pr X (y + ∆)∥ 2 2 + 2⟨y + ∆, pr X (y + ∆) -x * ⟩ -∥x * ∥ 2 2 -∥ pr X (y)∥ 2 2 + 2⟨y, pr X (y) -x * ⟩ = ∥ pr X (y)∥ 2 2 -∥ pr X (y + ∆)∥ 2 2 -2⟨y, pr X (y) -x * ⟩ + 2⟨y + ∆, pr X (y + ∆) -x * ⟩ = ∥ pr X (y)∥ 2 2 -∥ pr X (y + ∆)∥ 2 2 + 2⟨y, pr X (y + ∆) -pr X (y)⟩ + 2⟨∆, pr X (y) -x * + pr X (y + ∆) -pr X (y)⟩ = 2⟨∆, pr X (y) -x * ⟩ + 2⟨y + ∆, pr X (y + ∆) -pr X (y)⟩ + ∥ pr X (y)∥ 2 2 -∥ pr X (y + ∆)∥ 2 2 = 2⟨∆, pr X (y) -x * ⟩ + ∥ pr X (y) -(y + ∆)∥ 2 2 -∥ pr X (y + ∆) -(y + ∆)∥ 2 2 ≤ 2⟨∆, pr X (y) -x * ⟩ + ∥∆∥ 2 2 (B.17)
where the last equality follows from completing the squares and the last inequality follows from the first part of the lemma. ■ B.2.3. Proof of Proposition 1. We provide the details for all the steps.

(1) Defining b n = ∇ f (x s(n) ) -∇ f (x n ), we can rewrite the gradient update in DAGD as:

y n+1 = y n -α n+1 ∇ f (x s(n) ) = y n -α n+1 ∇ f (x n ) -α n+1 {∇ f (x s(n) ) -∇ f (x n )} = y n -α n+1 (∇ f (x n ) + b n ). (B.18)
Recall here once again that s(n) denotes the previous iteration count whose gradient becomes available only at the current iteration n. To establish the claim, we start by expanding b n as follows:

∥b n ∥ 2 = ∥ ∇ f (x s(n) ) -∇ f (x n )∥ 2 ≤ C 3 ∥x s(n) -x n ∥ 2 = C 3 ∥ pr X (y s(n) ) -pr X (y n )∥ 2 ≤ C 3 ∥y s(n) -y n ∥ 2 ≤ C 3 ∥y s(n) -y s(n)+1 ∥ 2 + ∥y s(n)+1 -y s(n)+2 ∥ 2 + • • • + ∥y n-1 -y n ∥ 2 = C 3 n-1 r=s(n) ∥α r+1 ∇ f (x s(r) )∥ 2 ≤ C 3 sup x∈X ∥ ∇ f (x)∥ 2 n-1 r=s(n) α r+1 = C 3 V max n-1 r=s(n) α r+1 . (B.19)
We now consider two cases, depending on whether the delays are bounded or not. where the limit approaching 0 follows from lim n→∞ α n = 0, which itself is a consequence of Assumption 2. This implies lim n→∞ C 3 V max n-1 r=s(n) α r+1 = 0 and consequently, lim n→∞ ∥b n ∥ 2 = 0. (b) We consider each of the three conditions in turn.

When α n-1 = 1 n and d n = o(n), we have n-s(n) ≤ Ko(s(n)) for some universal constant K > 0, which means n ≤ s(n) + Ko(s(n)). Consequently, we have:

0 < n-1 r=s(n) α r+1 = s(n)+Ko(s(n)) r=s(n) α r ≤ s(n)+Ko(s(n)) s(n) 1 r dr = log(s(n) + Ko(s(n))) -log s(n) = log Ko(s(n)) + s(n) s(n) → log 1 = 0 as n → ∞, (B.21)
where the last limit follows from s

(n) → ∞ as n → ∞ because n ≤ s(n) + Ko(s(n)).
When α n-1 = 1 n log n and d n = O(n), it is easy to verify (by integration) that this particular choice of sequence satisfies

∞ n=1 α 2 n < ∞, ∞ n=1 α n = ∞. Since d n = O(n),
we have n -s(n) ≤ Ks(n) for some universal constant K > 0, which means n ≤ s(n) + Ks(n). Consequently, we have:

0 < n-1 r=s(n) α r+1 = s(n)+Ks(n) r=s(n) α r ≤ s(n)+Ks(n) s(n) 1 r log r dr = log log(s(n) + Ks(n)) log s(n) = log log(K + 1) + log s(n) log s(n) → 0 as n → ∞, (B.22)
where the last limit follows from s(n) → ∞ as n → ∞ because n ≤ s(n)+Ks(n).

When α n-1 = 1 n log n log log n and d n = O(n a ), a > 1, it is again easy to verify (by integration) that this particular choice of sequence satisfies Assumption 2. Since d n = O(n a ), we have n -s(n) ≤ Ks(n) a for some universal constant K > 0, which means n ≤ s(n) + Ks(n) a . Consequently, we have:

0 < n-1 r=s(n) α r+1 = s(n)+Ks(n) a r=s(n) α r ≤ s(n)+Ks(n) a s(n) 1 r log r log log r dr = log log log(s(n) + Ks(n) a ) log log s(n) ≤ log log log(s(n) a + Ks(n) a ) log log s(n) < log log((K + 1) log s(n) + a log s(n)) log log s(n) = log log(K + 1 + a) + log log s(n) log log s(n) → 0 as n → ∞, (B.23)
where the last limit follows from the fact that s(n) → ∞ as n → ∞ (again, because n ≤ s(n) + Ks(n) a ). ( 2) With the defintion of b n , DAGD can be written as:

x n = pr X (y n ),

y n+1 = y n -α n+1 (∇ f (x n ) + b n ). (B.24)
To prove the claim, we start by fixing an arbitrary x * ∈ X * and applying Lemma 4 to bound the energy change in a single gradient update as follows:

E x * (y n+1 ) -E x * (y n ) ≤ 2⟨y n+1 -y n , pr X (y n ) -x * ⟩ + ∥y n -y n+1 ∥ 2 2 = -2⟨α n+1 (∇ f (x n ) + b n ), x n -x * ⟩ + ∥y n -y n+1 ∥ 2 2 .
(B.25) Now telescoping the above inequality yields:

E x * (y n+1 ) -E x * (y 0 ) = n r=0 {E x * (y r+1 ) -E x * (y r )} ≤ n r=0 {-2α r+1 ⟨∇ f (x r ) + b r , x r -x * ⟩ + α 2 r+1 ∥ ∇ f (x r ) + b r ∥ 2 2 } ≤ -2 n r=0 α r+1 {⟨∇ f (x r ), x r -x * ⟩ -∥b r ∥ 2 ∥x r -x * ∥ 2 } + 2 n r=0 α 2 r+1 {∥ ∇ f (x r )∥ 2 2 + ∥b r ∥ 2 2 } ≤ -2 n r=0 α r+1 {⟨∇ f (x r ), x r -x * ⟩ -C 4 ∥b r ∥ 2 } + 2 n r=0 α 2 r+1 {C 2 + B}, (B.26)
where the last inequality follows from the fact that b n 's must be bounded (since lim n→∞ ∥b n ∥ 2 = 0) and hence let B ≜ sup n ∥b n ∥ 2 . By Assumption 2, we have 2 n r=0 α 2 r+1 {C 2 + B} = B < ∞. Now fix any positive number ϵ. Assume for contradiction purposes x n only enters N (x * , ϵ) a finite number of times and let t 1 be the last time this occurs. This means that for all r > t 1 , x r is outside the open set N (X * , ϵ). Therefore, since a continuous function always achieves its minimum on a compact set, we have: ⟨∇ f (x r ), x r -x * ⟩ ≥ min x∈X -N (X * ,ϵ) ⟨∇ f (x), x -x * ⟩ ≜ a > 0, ∀r > t 1 (note that a depends on ϵ). Further, since b r → 0 as r → ∞, pick t 2 such that ∥b r ∥ 2 < a 2C4 , ∀r ≥ t 2 . Denoting t = max(t 1 , t 2 ), we continue the chain of inequalities in Equation (B.26) below:

-E x * (y 0 ) ≤ E x * (y n+1 ) -E x * (y 0 ) ≤ -2 t r=0 α r+1 {⟨∇ f (x r ), x r -x * ⟩ -C 4 ∥b r ∥ 2 } -2 n r=t+1 α r+1 {⟨∇ f (x r ), x r -x * ⟩ -C 4 ∥b r ∥ 2 } + 2 n r=0 α 2 r+1 {C 2 + B} ≤ -2 t r=0 α r+1 {⟨∇ f (x r ), x r -x * ⟩ -C 4 ∥b r ∥ 2 } -2 n r=t+1 α r+1 {a -C 4 ∥b r ∥ 2 } + B ≤ 2C 4 t r=0 α r+1 |b r ∥ 2 + B -2 n r=t+1 α r+1 {a - a 2 } = B -a n r=t+1 α r+1 → -∞, as n → ∞
where the first inequality follows from the energy function always being positive (Lemma 3), the second-to-last inequality follows from variational coherence and the limit on the last line follows from Assumption 2 and that B ≜ 2C 4 t r=0 α r+1 |b r ∥ 2 +B is just some finite constant. This yields an immediate contradiction and the claim is therefore established. ■ B.2.4. Proof of Theorem 2. Fix a given δ > 0. Since α n → 0, b n → 0 as n → ∞, for any a > 0, we can pick an N large enough (depending on δ and a) such that ∀n ≥ N , the following three statements all hold:

2BC 4 α n+1 + 2α 2 n+1 (C 2 + B 2 ) ≤ δ 2 , C 4 ∥b n ∥ 2 ≤ a 2 , α n+1 (C 2 + B 2 ) < a 2 . (B.27)
We show that under either of the following two (exhaustive) possibilities, if E(y n ) is less than δ, E(y n+1 ) is less than δ as well, where n ≥ N .

(1) Case 1:

E(y n ) < δ 2 . (2) Case 2: δ 2 ≤ E(y n ) < δ. Under Case 1, it follows from Equation (B.25): E x * (y n+1 ) -E x * (y n ) ≤ -2α n+1 ⟨∇ f (x n ) + b n , x n -x * ⟩ + α 2 n+1 ∥ ∇ f (x n ) + b n ∥ 2 2 ≤ -2α n+1 ⟨b n , x n -x * ⟩ + α 2 n+1 ∥ ∇ f (x n ) + b n ∥ 2 2 ≤ 2α n+1 ∥b n ∥ 2 ∥x n -x * ∥ 2 + 2α 2 n+1 (C 2 + B 2 ) ≤ 2BC 4 α n+1 + 2α 2 n+1 (C 2 + B 2 ) ≤ δ 2 , (B.28)
where the second inequality follows from variational coherence. Taking the infimum of x * over X * then yields: E(y n+1 )-E(y n ) ≤ δ 2 . This then implies that E(x * , y n+1 ) ≤ E(x * , y n )+ δ 2 < δ. Under Case 2, Eq. (B.25) readily yields:

E x * (y n+1 ) -E x * (y n ) ≤ -2α n+1 ⟨∇ f (x n ) + b n , x n -x * ⟩ + α 2 n+1 ∥ ∇ f (x n ) + b n ∥ 2 2 = -2α n+1 ⟨∇ f (x n ), x n -x * ⟩ -2α n+1 ⟨b n , x n -x * ⟩ + α 2 n+1 ∥ ∇ f (x n ) + b n ∥ 2 2 ≤ -2α n+1 a + 2α n+1 ∥b n ∥ 2 ∥x n -x * ∥ 2 + 2α 2 n+1 (C 2 + B 2 ) ≤ -2α n+1 a -C 4 ∥b n ∥ 2 -α n+1 (C 2 + B 2 ) ≤ -2α n+1 a - a 2 -α n+1 (C 2 + B 2 ) = -2α n+1 a 2 -α n+1 (C 2 + B 2 ) < 0, (B.29)
where the second inequality follows from ⟨∇ f (x n ), x n -x * ⟩ ≥ a under Case 2 15 . Taking the infimum of x * over X * then yields: E(y n+1 ) ≤ E(y n ). This then implies that E(y n+1 ) ≤ E(y n ) < δ 2 . Consequently, putting the above two cases together therefore yields E(y n+1 ) < δ. ■ B.2.5. Proof of Proposition 2. Using the definitions introduced in the main text, we rewrite the gradient update in DASGD as:

Y n+1 = Y n -α n+1 {∇ f (X n ) + B n + U n+1 }. (B.30) (1) To see that n r=0 U n+1 is a martingale adapted to Y 0 , Y 1 . . . , Y n+1 , first note that, by defintion, B n is adapted to Y 0 , Y 1 . . . , Y n (since X n is a deterministic function of Y n ) and Y n+1 , Y n , B n together determine U n+1 .
We then check that their first moments are bounded:

E[∥ n r=0 ∥U r+1 ∥ 2 ] ≤ n r=0 E[∥U r+1 ∥ 2 ] = n r=0 E[∥∇F (X s(r) , ω r+1 ) -∇ f (X s(r) )∥ 2 ] ≤ n r=0 E[∥∇F (X s(r) , ω r+1 )∥ 2 ] + E[∥ ∇ f (X s(r) )∥ 2 ] ≤ n r=0 sup x∈X E[∥∇F (x, ω)∥ 2 ] + sup x∈X ∥ ∇ f (x)∥ 2 = n r=0 sup x∈X E[∥∇F (x, ω)∥ 2 ] + sup x∈X ∥ E[∇F (x, ω)]∥ 2 ≤ n r=0 2 sup x∈X E[∥∇F (x, ω)∥ 2 ] = n r=0 2C 1 = 2(n + 1)C 1 < ∞, (B.31)
where the last inequality follows from Jensen's inequality (since ∥ • ∥ 2 is a convex function). Finally, the martingale property holds because:

E[ n r=0 U r+1 | Y 1 , . . . , Y n+1 ] = E[∇F (X s(n) , ω n+1 ) -∇ f (X s(n) ) | Y 1 , . . . , Y n ] + n-1 r=0 U r+1 = n-1 r=0 U r+1 ., since ω n+1 is iid. Therefore, n r=0
U r+1 is a martingale, and U n+1 is a martingale difference sequence adapted to Y 0 , Y 1 . . . , Y n+1 . 15 Here a is a constant that only depends on δ: if E(x * , y) ≥ δ 2 , then by part 2 of Lemma 3, pr X (y) must be outside an ϵ-neighborhood of x * , for some ϵ > 0. On this neighborhood, the strictly positive continuous function ⟨∇ f (x), x -x * ⟩ must achieve a minimum value a > 0.

Next, we show that lim n→∞ ∥B n ∥ 2 = 0, a.s.. By definition, we can expand B n as follows:

∥B n ∥ 2 = ∥ ∇ f (X s(n) ) -∇ f (X n )∥ 2 ≤ C 3 ∥X s(n) -X n ∥ 2 = C 3 ∥ pr X (Y s(n) ) -pr X (Y n )∥ 2 ≤ C 3 ∥Y s(n) -Y n ∥ 2 = C 3 Y s(n) -Y s(n)+1 + Y s(n)+1 -Y s(n)+2 + • • • + Y n-1 -Y n 2 = C 3 n-1 r=s(n) Y r -Y r+1 2 = C 3 n-1 r=s(n) α r+1 ∇F (X s(r) , ω r+1 ) 2 = C 3 n-1 r=s(n) α r+1 ∇ f (X s(r) ) + ∇F (X s(r) , ω r+1 ) -∇ f (X s(r) ) 2 = C 3 n-1 r=s(n) α r+1 ∇ f (X s(r) ) + n-1 r=s(n) α r+1 U r+1 2 ≤ C 3 n-1 r=s(n) α r+1 ∥ ∇ f (X s(r) )∥ 2 + C 3 n-1 r=s(n) α r+1 U r+1 2 ≤ C 3 C 1 n-1 r=s(n) α r+1 + C 3 ∥ n-1 r=s(n) α r+1 U r+1 ∥ 2 = C 3 C 1 n-1 r=s(n) α r+1 + C 3 ∥ n-1 r=0 α r+1 U r+1 - s(n)-1 r=0 α r+1 U r+1 ∥ 2 , (B.32)
where the first inequality follows from ∇ f being Liptichz-continuous (Assumption 3) and the second inequality follows from pr X is a non-expansive map.

By the same analysis as in the deterministic case, the first part of the last line of Equation (B.32) converges to 0 (under each one of the conditions on step-size and delays in Assumption 2):

lim n→∞ C 3 C 1 n-1 r=s(n) α r+1 = 0.
(B.33)

We then analyze the limit of ∥

n-1 r=0 α r+1 U r+1 - s(n)-1 r=0 α r+1 U r+1 ∥ 2 . Define: M n = n-1 r=0 α r+1 U r+1 .
Since U r+1 's are martingale differences, M n is a martingale. Further, in each of the three conditions,

∞ n=1 α 2 n < ∞. This implies that M n is an L 2 -bounded martingale because: sup n E[∥M n ∥ 2 2 ] = sup n E[ n-1 r=0 α r+1 U r+1 2 2 ] = sup n E[⟨ n-1 r=0 α r+1 U r+1 , n-1 r=0 α r+1 U r+1 ⟩] = sup n E[ i,j ⟨α i+1 U i+1 , α j+1 U j+1 ⟩] = sup n n-1 r=0 E[⟨α r+1 U r+1 , α r+1 U r+1 ⟩] = sup n n-1 r=0 α 2 r+1 E[∥U r+1 ∥ 2 2 ] ≤ sup n 4C 2 n-1 r=0 α 2 r+1 ≤ 4C 2 ∞ r=0 α 2 r+1 < ∞, (B.34)
where the last inequality in the second line follows from the martingale property as follows:

E[⟨α i+1 U i+1 , α j+1 U j+1 ⟩] = α i+1 α j+1 E[⟨U i+1 , U j+1 ⟩] = α i+1 α j+1 E[E[⟨U i+1 , U j+1 ⟩ | Y 0 , Y 1 , . . . , Y i+1 ]] = α i+1 α j+1 E[⟨U i+1 , E[U j+1 | Y 0 , Y 1 , . . . , Y i+1 ]⟩] = α i+1 α j+1 E[⟨U i+1 , 0⟩] = 0, (B.35)
where the second equality follows from the tower property (and without loss of generality, we have assumed i < j, the third equality follows from U i+1 is adapted to Y 0 , Y 1 , . . . , Y i+1 and the second-to-last equality follows from U n+1 is a martingale difference. Consequently, all the cross terms in the second line of Equation (B.34) are 0. Therefore, by Lemma A.2, by taking p = 2 lim n→∞ M n = M ∞ , a.s., where M ∞ has finite second-moment. Further, since in all three cases s(n) → ∞ as n → ∞ (because there is at most a polynomial lag between s(n) and n), we have

lim n→∞ M s(n) = M ∞ , a.s.. Therefore lim n→∞ n-1 r=0 α r+1 U r+1 - s(n)-1 r=0 α r+1 U r+1 = lim n→∞ M n -M s(n) = 0, a.s.,
thereby implying:

lim n→∞ C 3 n-1 r=0 α r+1 U r+1 - s(n)-1 r=0 α r+1 U r+1 2 = 0. (B.36)
Combining Equation (B.33) and Equation (B.36) yields lim n→∞ ∥B n ∥ 2 = 0, a.s.. (2) The full DASGD update is then:

X n = pr X (Y n ) (B.37) Y n+1 = Y n -α n+1 {∇ f (X n ) + B n + U n+1 }. (B.38)
We now bound the one-step change of the energy function E(X * , Y n+1 ) -E(X * , Y n ) (which is now a random quantity) and then telescope the differences. Pick an arbitrary x * ∈ X * and apply Lemma 4, we have:

E x * (Y n+1 ) -E x * (Y n ) ≤ -2α n+1 ⟨∇ f (X n ) + B n + U n+1 , X n -x * ⟩ + ∥Y n -Y n+1 ∥ 2 2 = -2α n+1 ⟨∇ f (X n ) + B n + U n+1 , X n -x * ⟩ + α 2 n+1 ∥ ∇ f (X n ) + B n + U n+1 ∥ 2 2 ≤ -2α n+1 ⟨∇ f (X n ) + B n + U n+1 , X n -x * ⟩ + 3α 2 n+1 ∥ ∇ f (X n )∥ 2 2 + ∥B n ∥ 2 2 + ∥U n+1 ∥ 2 2 ≤ -2α n+1 ⟨∇ f (X n ) + B n + U n+1 , X n -x * ⟩ + 3α 2 n+1 C 2 + ∥B n ∥ 2 2 + ∥U n+1 ∥ 2 2 . (B.39)
For contradiction purposes assume X n enters N (X * , ϵ) only a finite number of times with positive probability. By starting the sequence at a later index if necessary, we can without loss of generality X n never enters N (x * , ϵ) with positive probability. Then on this event (of X n never entering N (X * , ϵ)), we have ⟨∇ f (X n ), X n -x * ⟩ ≥ a > 0 as before. Telescoping Equation (B.39) then yields:

-∞ < -E x * (Y 0 ) ≤ E x * (Y n+1 ) -E x * (Y 0 ) = n r=0 {E x * (Y r+1 ) -E x * (Y r )} ≤ -2 n r=0 α n+1 ⟨∇ f (X n ) + B n + U n+1 , X n -x * ⟩ + 3 n r=0 α 2 n+1 C 2 + ∥B n ∥ 2 2 + ∥U n+1 ∥
sup n E[ n r=0 α 2 n+1 ∥U n+1 ∥ 2 2 ] ≤ sup n n r=0 α 2 n+1 E[∥U n+1 ∥ 2 2 ] ≤ sup n n r=0 α 2 n+1 E[∥U n+1 ∥ 2 2 ] = sup n n r=0 α 2 n+1 E[∥∇F (X s(n) , ω n+1 ) -∇ f (X s(n) )∥ 2 2 ]
≤ 2 sup In addition, α n+1 ⟨U n+1 , X n -x * ⟩ is a martingale difference that is L 2 bounded because α n+1 is square summable and

E[∥⟨U n+1 , X n -x * ⟩∥ 2 2 ] ≤ E[∥U n+1 ∥ 2 2 ∥X n -x * ∥ 2 2 ] ≤ C 4 E[∥U n+1 ∥ 2 2 ] ≤ 4C 4 C 2 < ∞.
Consequently, by applying Lemma A.5 with p = 2 and u n = n r=1 α r+1 (which is not summable), law of large number for martingales therefore implies: where the last inequality is strict unless pr X (y(t)) = x(t) = x * . Take the infimum over x * then yields the result: in particular, if pr X (y(t)) = x(t) / ∈ X * , then d dt E x * (y(t)) ≤ -ϵ < 0, ∀x * ∈ X * , hence yielding the strict part of the inequality. Note also that even though pr X (y(t)) -y(t) is not differentiable, ∥ pr X (y(t)) -y(t)∥ 2 2 is; and in computing its derivative, we applied the envelope theorem as given in Lemma A.3.

For the second claim, consider all y that satisfy E(P (t, y)) > δ 2 . Fix any x * ∈ X * . By the monotonicity property in the first part of the lemma, it follows that E x * (P (s, y)) > δ 2 , ∀0 ≤ s ≤ t. Consequently, P (s, y) must be outside some ϵ neighborhood of X * for 0 ≤ s ≤ t, for otherwise, E x * (P (t, y)) would be 0 for at least some x * ∈ X * , which is a contradiction to E(P (t, y)) > δ 2 . This means that there exists some positive constant a(δ) such that ∀0 ≤ s ≤ t, ∀x * ∈ X * : First, from previous analysis, we already know that U n+1 is a martingale difference sequence with sup n E[∥U n+1 ∥ 2 2 ] < ∞ and that α n is a square-summable sequence. Recall also lim n→∞ B n = 0 (a.s.). Now, by combining Proposition 4.1 and Proposition 4.2 in Benaïm [START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF], we obtain the following sufficient condition for APT:

The affine interpolation curve of the iterates generated by the difference equation Y n+1 = Y n -α n+1 {G(X n ) + U n+1 } is an APT for the solution to the ODE ẏ = -G(y) if the following three conditions all hold:

(1) G is Lipschitz continuous and bounded. 16 .

(2) U n+1 is a martingale difference sequence with sup n E[∥U n+1 ∥ p 2 ] < ∞ and ∞ n=0 α

1+ p 2 n+1 < ∞ for some p > 2.
By using a similar analysis as in Benaïm [START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF] (which we omit here due to space limitation), we can show that if the following conditions all hold:

(1) G is Lipschitz continuous and bounded.

(2) lim n→∞ B n = 0 (a.s.).

(3) U n+1 is a martingale difference sequence with

sup n E[∥U n+1 ∥ p 2 ] < ∞ and ∞ n=0 α 1+ p 2
n+1 < ∞ for some p > 2, then, the affine interpolation curve of the iterates generated by the difference equation Y n+1 = Y n -α n+1 {G(X n ) + B n + U n+1 } is an APT for the solution to the ODE ẏ = -G(y).

Take p = 2 and recall ∇ f (pr X (•)) is Lipschitz continuous and bounded. The above list of three conditions are thus all verified, thereby yielding the result. B.3.1. Proof of Theorem 3. By Proposition 2, Y n will get arbitrarily close to X * infinitely often. It then suffices to show that, after long enough iterations, if Y n ever gets ϵ-close to X * , all the ensuing iterates will be ϵ-close to X * almost surely. The way we show this "trapping" property is to use the energy function. Specifically, we consider E(x * , A(t)) and show that no matter how small ϵ is, for all sufficiently large t, if E(x * , A(t 0 )) is less than ϵ for some t 0 , then E(x * , A(t)) < ϵ, ∀t > t 0 . This would then complete the proof because A(t) actually contains all the DASGD iterates, and hence if E(x * , A(t)) < ϵ, ∀t > t 0 , then E(x * , Y n ) < ϵ for all sufficiently large n. Furthermore, since A(t) contains all the iterates, the hypothesis that " if E(x * , A(t 0 )) is less than ϵ for some t 0 " will be satisfied due to Prop 2.

We now flesh out more details of the proof. Fix any ϵ > 0. Since A(t) is an asymptotic pseudotrajectory for P , we have: Consequently, for any δ > 0, there exists some τ (δ, T ) such that ∥Y (t + h) -P (h, Y (t))∥ 2 < δ for all t ≥ τ and all h ∈ [0, T ]. We therefore have the following chain of inequalities: where in the last step we have choosen δ small enough such that C 4 δ + 1 2 δ 2 = ε 2 . Now by Proposition 2, there exists some τ 0 such that E x * (A(τ 0 )) < ε 2 . Our goal is to establish that E x * (A(τ 0 + h)) < ε for all h ∈ [0, ∞). To that end, partition the [0, ∞) into disjoint time intervals of the form [(n -1)T ε , nT ε ) for some appropriate T ε . By Lemma 5, we have: E x * (P (h, A(τ 0 ))) ≤ E x * (P (0, A(τ 0 ))) = E x * (A(τ 0 )) < ε 2 for all h ≥ 0. (B.48)

E x * (A(t + h)) = E x
Consequently:

E x * (A(τ 0 + h)) < E x * (P (h, A(τ 0 ))) + ε 2 < ε 2 + ε 2 = ε, (B.49)
where the last inequality is a consequence of (B.48). Now, assume inductively that Eq. (B.49) holds for all h ∈ [(n -1)T ε , nT ε ) for some n ≥ 1. Then, for all h ∈ [(n -1)T ε , nT ε ), we have:

E x * (A(τ 0 + T ε + h)) < E x * (P (T ε , A(τ 0 + h))) + ε 2 ≤ max ε 2 , E x * (A(τ 0 + h)) - ε 2 + ε 2 ≤ ε 2 + ε 2 = ε. (B.50)
Consequently, Eq. (B.49) holds for all h ∈ [nT ε , (n + 1)T ε ). This completes the induction. Taking the infimum over x * then completes our proof.

Figure 1 .

 1 Figure 1. Two commonly used distributed computing architecutures: (a) masterworker (left) and multi-processorwith shared memory (right).

2 : repeat 3 :

 23 (a) Each proceesor reads the current global iterate. (b) Each processor reads data from memory and computes a stochastic gradient. (c) Each processor updates the global iterate. 4: until end

Algorithm 3 .

 3 Distributed asynchronous stochastic gradient descentRequire: Initial state X 0 ∈ R d , step-size sequence α n 1: n ← 0; 2: repeat 3:

Figure 2 .

 2 Figure 2. Examples of variationally coherent objectives: on the top row, the Beale function f (x1, x2) = (1.5 -x1 + x1x2) 2 + (2.25 -x1 + x1x 2 2 ) 2 + (2.625 -x1 + x1x 3 2 ) 2 over the benchmark domain [-4, 4] × [-4, 4]; on the bottom row, the polar example f (r, θ) = (3 + sin(5θ) + cos(3θ))r 2 (5/3 -r) over the unit ball 0 ≤ r ≤ 1. In both figures, the black curves indicate a sample trajectory of DASGD with linearly growing delays.

Algorithm 5 .

 5 Distributed Asynchronous Gradient Descent (DAGD) Require: Initial state y 0 ∈ R d , step-size sequence α n 1: n ← 0 2: repeat 3:

Figure 3 .

 3 Figure 3. Value convergence in a non-convex stochastic optimization problem with d = 1001 degrees of freedom.

Lemma A. 3 .. 2 )

 32 Let f : R n × R m → R be a continuously differentiable function. Let U be a compact set and consider the problem max x∈U f (x, θ). (ALet x * : O → R m be a continuous function defined on an open set O ⊂ R m such that for each θ ∈ O, x * (θ) solves the problem in Equation A.2. Define V : R m → R where V (θ) = f (x * (θ), θ). Then V (θ) is differentiable on O and:

  (a) If {α n } ∞ n=1 and d n ≤ D, ∀n satisfy Assumption 2, then d s(n) = n -s(n) ≤ D. α r ≤ D max r∈{n-D,...,n} α r → 0 as n → ∞, (B.20)

2 ]2C 2 3 n r=0 α 2 α

 2232 ∥ ∇ F (X s(n) , ω n+1 )∥ 2 2 ] + E[∥ ∇ f (X s(n) )∥ 2 α 2 n+1 < ∞.Consequently, by martingale convergence theorem (Lemma A.2 by taking p = 1), n+1 ∥U n+1 ∥ 2 2 → R, a.s., for some random variable R that is almost surely finite (in fact E[|R|] < ∞). (c) Since ∥B n ∥ 2 converges to 0 almost surely, its average also converges to 0 almost surely:n+1 ⟨B n , X n -x * ⟩ n r=1 α r+1 = 0, a.s., since |⟨B n , X n -x * ⟩| ≤ ∥B n ∥ 2 ∥X n -x * ∥ 2 ≤ C 4 ∥B n ∥ 2 .

∞α

  n=0 n+1 ⟨U n+1 , X n -x * ⟩ n r=1 α r+1 = 0, a.s.Combining the above two limits, we havelim n→∞ n r=0 α n+1 ⟨B n + U n+1 , X n -x * ⟩ n r=0 α r+1 = 0, a.s. Consequently, -n r=0 α n+1 a + ⟨B n + U n+1 , X n -x * ⟩ = -{ n r=0 α n+1 } a + n r=0 αn+1⟨Bn+Un+1,Xn-x * ⟩ n r=0 αr+1→ -∞, as n → -∞.

2 . (B. 43 )B. 3 .

 2433 d ds E x * (P (s, y)) = -⟨∇ f (x(s)), x(s) -x * ⟩ ≤ -a(δ). (B.42)Consequently, pick T (δ) = δ 2a(δ) , Equation (B.42) implies that for any t > T (δ):E x * (P (t, y)) ≤ E x * (P (T (δ), y)) ≤ E x * (y) -T (δ)a(δ) ≤ E x * (y) -δTaking the supremum over x * ∈ X * then yields:E(P (t, y)) ≤ E(y) -δ 2 . (B.44)Since Equation (B.44) is true for any y, taking sup over y establishes the claim. ■ Proof of Lemma 6.

lim t→∞ sup 0≤h≤T ∥Y

 0≤h≤T (t + h) -P (h, Y (t))∥ 2 = 0. (B.[START_REF] Shapiro | A tutorial on stochastic programming[END_REF]

  ϵ), ∀k. However, on this subsequence, the energy function still converges to 0 by assumption: E(y n k ) → 0, which immediately yields a contradiction. The claim is hence established. ■ B.2.2. Proof of Lemma 4. We first prove the first claim. By expanding it, we have:

  + ⟨B n + U n+1 , X n -x * ⟩ + 3We justify the last-line limit of Equation (B.40) by looking at each of its components in turn: (a) Since n r=0 α 2 n+1 < ∞, and per the previous step, lim n→∞ ∥B n ∥ 2 2 = 0, a.s., we have 3∞ r=0 α 2 n+1 C 2 + ∥B n ∥ 2 2 = C, a.s., for some constant C < ∞. (b) n r=0 α 2 n+1 ∥U n+1 ∥ 22 is submartingale that is L 1 bounded since:

				2
				2
	≤ -2	n	α n+1 a n	α 2 n+1 C 2 + ∥B n ∥ 2 2 + ∥U n+1 ∥ 2 2
		r=0	r=0
	→ -∞ a.s. as n → ∞.
				(B.40)

  B.2.6. Proof of Lemma 5. The first claim follows by computing the derivative of the energy function with respect to time (for notational simplicity, here we just use y(t) to denote P (t, y)):

	d dt	E x * (y(t)) =	d dt	∥x

* ∥ 2 2 -∥ pr X (y(t))∥ 2 2 + 2⟨y(t), pr X (y(t)) -x * ⟩ = d dt -∥ pr X (y(t)) -y(t)∥ 2 2 + ∥y(t)∥ 2 2 + 2⟨y(t), -x * ⟩ = 2⟨ ẏ(t), pr X (y(t)) -y(t)⟩ + 2⟨y(t), ẏ(t)⟩ + 2⟨ ẏ(t), -x * ⟩ = -2⟨∇ f (x(t)), x(t) -y(t)⟩ -2⟨∇ f (x(t)), y(t)⟩ -2⟨∇ f (x(t)), -x * ⟩ = -⟨∇ f (x(t)), x

(

t) -x * ⟩ ≤ 0, (B.41)

This is again done by sampling a subset of the training data in the global memory and computing the gradient at the iterate for each datapoint and averaging over all the comptued gradients in the sample.

See Appendix A for a more rigorous definition. Furthermore, R + is the set of non-negative reals.

This condition is also sufficient for the existence and uniqueness of the ODE solution

Acknowledgments

The authors are grateful to an associate editor and two anonymous referees for their suggestions and remarks, and to John Tsitsiklis for his insightful input on previous work. The first author in particular enjoyed the conversations with John on this topic during the winter of 2019 when visiting MIT.

Zhengyuan Zhou was partially supported by IBM Goldstine fellowship. P. Mertikopoulos was partially supported by the COST Action CA16228 "European Network for Game Theory" (GAMENET). P. Mertikopoulos also received financial support from the French National Research Agency (ANR) in the framework of the "Investissements d'avenir" program (ANR-15-IDEX-02), the LabEx PERSYVAL (ANR-11-LABX-0025-01), MIAI@Grenoble Alpes (ANR-19-P3IA-0003), and the grants ORACLESS (ANR-16-CE33-0004) and ALIAS (ANR-19-CE48-0018-01).