ON THE EXPECTED RUNNING TIME OF NONCONVEX OPTIMIZATION WITH EARLY STOPPING
Year: 2019

Abstract:
This work examines the convergence of stochastic gradient algorithms that use early stopping based on a validation function, wherein optimization ends when the magnitude of a validation function gradient drops below a threshold. We derive conditions that guarantee this stopping rule is well-defined and analyze the expected number of iterations and gradient evaluations needed to meet this criteria. The guarantee accounts for the distance between the training and validation sets, measured with the Wasserstein distance. We develop the approach for stochastic gradient descent (SGD), allowing for biased update directions subject to a Lyapunov condition. We apply the approach to obtain new bounds on the expected running time of several algorithms, including Decentralized SGD (DSGD), a variant of decentralized SGD, known as Stacked SGD, and the stochastic variance reduced gradient (SVRG) algorithm. Finally, we consider the generalization properties of the iterate returned by early stopping.

heading: 1 INTRODUCTION
text: This work considers the minimization of a differentiable and possible nonconvex objective function:
min x∈Rd f(x). (1)
A generally accepted success criteria for algorithms that use only first-order information is that an approximate stationary point is generated. These are points x ∈ Rd at which the function f has a small gradient. In a typical machine learning scenario, f is the average loss over a dataset of training examples, and the method of choice involves using some form of stochastic gradient method, for instance, stochastic gradient descent, or SGD (see Algorithm 1). The success of SGD in machine learning problems has led to many extensions of the algorithm, including variance-reduced and distributed variants (reviewed in Section 1.1).
A common approach to stopping optimization in practice is to use early stopping based on a validation function. In this scenario, a stopping criterion is periodically evaluated on the validation function, and the algorithm stops once this criterion is met. The validation and training sets often are disjoint (although this is not required in the present work). Although this approach is used frequently, there is little theoretical work on the runtime of nonconvex optimization using early stopping based on a validation function. In general, the runtime and performance will depend on several factors, including the relation between the validation and training functions, and the desired level of solution accuracy. In this work, the stopping criterion is that the algorithm has generated a point, which is an approximate stationary point for the validation function. Our analysis focuses on bounding the number of iterations and gradient evaluations used until the algorithm meets the stopping criteria. Formally, we consider the stopping time defined as the first time an iterate has the property of being an approximate stationary point for the validation function, and we derive upper bounds on the expected value of this stopping time. Using bounds on the Wasserstein distance between the training and validation sets, we also may derive a bound on the stationary gap of the training function at the resulting iterate. As an extension, we also describe how Wasserstein concentration bounds can be used to bound the stationarity gap with respect to the testing distribution to which both the training and validation datasets are drawn.
We apply our analysis to several settings, including stochastic gradient descent (SGD), stochastic variance reduced gradient (SVRG), decentralized SGD (DSGD) and stacked SGD (SSGD). The
result is new bounds on the expected number of Incremental First-order Oracle (IFO) calls needed to generate approximate stationary points for some known algorithms (SGD; DSGD; SVRG), as well as a new algorithm (SSGD).
Main contributions Our main contributions include:
– We present a non-asymptotic analysis of SGD with early stopping that leads to a bound on the expected number of gradient evaluations needed to find approximate stationary points of the training function (Corollary 3.5). The analysis allows for biases in the update direction, subject to a Lyapunov-type inequality on the error terms.
– We rigorously analyze the expected running time of two distributed SGD algorithms: Decentralized SGD, and Stacked SGD (Algorithm 2), a new decentralized form of SGD designed to exploit connectivity patterns consisting of a network of parameter-server type clusters (Corollary 4.4.)
– We apply the analysis to nonconvex SVRG to obtain a bound on the expected number of IFO calls for this algorithm (Corollary 6.2).
– We demonstrate how Wasserstein concentration bounds can be leveraged to bound the generalization performance of parameters returned by SGD with early stopping (Corollary 7.2).

heading: 1.1 RELATED WORK
text: The study of stochastic optimization goes back (at least) to the pioneering efforts of Robbins and Monro Robbins & Monro (1951). Subsequent developments include the ordinary differential equation (ODE) method Ljung (1977) and stochastic approximation Kushner & Clark (1978), which emphasizes the asymptotic behavior of the algorithms.
The subject of non-asymptotic performance guarantees has attracted interest as well, including lower bounds on algorithm performance Nemirovski & Yudin (1983). For a review of non-asymptotic guarantees for convex SGD, the reader may consult Nemirovski et al. (2009); Rakhlin et al. (2012); Bach & Moulines (2013). Extensions such as distributed Zinkevich et al. (2010) and asynchronous Zhang et al. (2015); Agarwal & Duchi (2011) variants of the algorithms also have been investigated. Many machine learning optimization problems involve an objective represented as a finite sum of functions, and, for this case, variance reduction techniques lead to improved rates of convergence over SGD Johnson & Zhang (2013); L. Roux et al. (2012); Defazio et al. (2014).
The randomized stochastic gradient (RSG) method Ghadimi & Lan (2013) uses randomization to obtain a non-asymptotic performance guarantee for SGD. The randomization technique has became a standard tool for analyzing optimization algorithms in the nonconvex setting Ghadimi & Lan (2016); Lian et al. (2017); Allen-Zhu (2018a;b); Reddi et al. (2016b;a); Zhang et al. (2016); Reddi et al. (2016b); Lei et al. (2017); Lian et al. (2015). Follow-up works have included analysis of nonconvex optimization in more sophisticated algorithmic settings, such as asynchronous Lian et al. (2015) and decentralized Lian et al. (2017) optimization. Analysis of variance reduced optimization has extended beyond convex functions beginning with an application to principal components analysis Shamir (2015) and later to general nonconvex functions Allen-Zhu & Hazan (2016); Reddi et al. (2016b;a); Lei et al. (2017). A highlight in this area is that the IFO complexity of SVRG for nonconvex functions is superior to that of RSG Allen-Zhu & Hazan (2016); Reddi et al. (2016a). Algorithms with better convergence rates than SVRG also have been developed Allen-Zhu (2018a;b).
We are particularly interested in results for biased SGD, e.g., Bertsekas & Tsitsiklis (2000) which considers the asymptotic convergence of biased SGD to stationary points. In this work, we are interested in a similar scenario but instead focus on the non-asymptotic behavior of the algorithm.
Several recent works also have explored the average amount of resources needed to reach a desired performance level in optimization. The expected running time of a stochastic trust region algorithm (STORM) is given in Blanchet et al. (2016). A similar methodology also has been used to analyze stochastic line search methods Paquette & Scheinberg (2018). Our convergence analysis is similar in spirit to these works as we also are interested in the expected amount of time or other resources required to meet the performance guarantee. However, the algorithms and assumptions in this work differ.
Other recent work has analyzed the theoretical aspects of early stopping. For instance, Duvenaud et al. (2016) developed an interpretation of early stopping in terms of variational Bayesian inference. Early stopping for a least squares problem in a reproducing kernel Hilbert space has been treated in Lin & Rosasco (2016). The implications of early stopping generalization were studied in Hardt et al. (2016). However, to our knowledge, this work is the first to analyze the algorithms’ runtimes when using a validation function for early stopping in nonconvex optimization.

heading: 2 PRELIMINARIES
text: Let f : Rq×Rd → R be a loss function whose value we denote by f(y, x). Intuitively, the variable y represents an input/output pair, and x represents the model parameters. Throughout, we shall assume the gradient of the objective function with respect to x is Lipschitz continuous (defined as follows). Assumption 2.1. The function f : Rq × Rd → R is bounded from below by f∗ ∈ R, and the derivative of f with respect to x is L-Lipschitz continuous:
∀y ∈ Rq, x1, x2 ∈ Rd, ‖∇xf(y, x1)−∇xf(y, x2)‖ ≤ L‖x1 − x2‖.
At times, we will make a distinction between the training function fT , which is used to calculate gradients, and a validation function fV used to decide when to stop training. Assumption 2.2. The training function fT is defined using a set YT ⊆ Rq of nT elements as fT (x) =
1 nT ∑ y∈YT f(y, x), and the validation function fV is defined using a set YV ⊆ Rq of nV
elements as fV (x) = 1nV ∑ y∈YV f(y, x).
To guarantee that the early stopping rule leads to a well-defined algorithm, we will assume a relation between the training and validation functions. Intuitively, the functions fT and fV will be similar when the datasets YT and YV are similar. Formally, the datasets YT and YV determine probability measures µT and µY , defined as µT = 1nT ∑ y∈YT δy and µV = 1 nV ∑ y∈YV δy, respectively, where δy is the delta measure δy(A) = 1y∈A for all sets A. We can compare these measures using the Wasserstein distance as in:
For q ≥ 1, p ≥ 1, we denote by Pp(Rq) the probability measures on Rq with finite moments of order p. Recall that a coupling of probability measures µ1 and µ2 is a probability measure γ on Rq × Rq such that for all measurable sets A, γ(A × Rq) = µ1(A) and γ(Rq × A) = µ2(A). Intuitively, a coupling transforms data distributed like µ1 into a dataset that is distributed according to µ2. The p-Wasserstein distance on Pp(Rq), denoted by dp, is defined as:
dp(µ1, µ2) = inf γ∈Γ(µ1,µ2)
( E
(x1, x2)∼γ
[ ‖x1 − x2‖p ])1/p , (2)
where Γ(µ1, µ2) is the set of all couplings of µ1 and µ2. For more details the reader is referred to Villani (2008). Assumption 2.3. There is a constant G ≥ 0 such that
∀x ∈ Rd, ‖∇fV (x)−∇fT (x)‖ ≤ Gd1(µV , µT ).
There are several cases in which this assumption will be satisfied. It is trivially satisfied if the training and validation sets are the same as µT = µV . As a consequence of the Kantorovich duality formula (Villani (2008), Remark 6.5), it is also satisfied if the function y 7→ ∇xf(y, x) is a G-Lipschitz function, uniformly for all x. Consider the following example: Example 2.4. Suppose that g : Rq × Rd → R is a smooth function. Let h : Rd → R be the function that applies the hyperbolic tangent function to each of its components: h(x) = (tanh(x1), . . . , tanh(xd)), and define f(y, x) = g(y, h(x)). Further suppose that the training data are bounded: ‖y‖ ≤ J for all y ∈ YT . Then Assumption 2.3 is satisfied because y 7→ ∇xf(x, y) is a Lipschitz function with G = sup‖y‖≤J,‖x‖≤√d ‖ ∂ 2g ∂x∂y (y, x)‖.
In our analyses the notion of success is that an algorithm generates an approximate stationary point: Definition 2.5. A point x ∈ Rd is an -approximate stationary point of f if ‖∇f(x)‖2 ≤ .
We measure the complexity of algorithms according to how many function value and gradient queries they make. Formally, an IFO) is defined as follows (Agarwal & Bottou, 2015): Definition 2.6. An IFO takes a parameter x and an input y and returns the pair (f(y, x),∇xf(y, x)).
In Appendix A.1, we briefly recall the notion of filtration, stopping times, and other concepts from stochastic processes that will be used in our analyses. We also refer readers to Williams (1991) for more details.
3 BIASED SGD
Algorithm 1 SGD with early stopping 1: input: Initial point x1 ∈ Rd 2: t = 1 3: while ‖∇fV (xt)‖2 > do 4: for n = t to t+m− 1 do 5: xn+1 = xn − ηhn 6: end 7: t = t+m 8: end 9: return xt This section details our analysis of SGD with early stopping shown on the right in Algorithm 1. Starting from an initial point x1, at each iteration, the parameter is updated with an approximate gradient hn using a step-size η. The gradient norm of the validation function is evaluated every m iterations, and the algorithm ends when the norm decreases below a threshold . We assume that the update direction ht is a sum of two components, vt and ∆t, that represent an unbiased gradient estimate and an error term, respectively:
ht = vt + ∆t. (3)
Let {Ft}t≥0 be a filtration such that x1 is F0-measurable, and for all t > 1, the variables vt,∆t are Ft-measurable. Our assumptions on the vt are as follows. Assumption 3.1. For any t ≥ 1, it holds that
E [vt −∇fT (xt) | Ft−1] = 0, (4) E [ ‖vt −∇fT (xt)‖2 | Ft−1 ] ≤ σ2v . (5)
Assumption 3.1 states that the update directions vt are valid approximations to the gradient∇fT (xt), and it also bounds the error in the approximation. For the variables ∆t we assume the following: Assumption 3.2. There is a sequence of random variables V1, V2, . . ., and U1, U2, . . . such that for all t ≥ 1 the variables Vt and Ut are Ft-measurable, ‖∆t‖2 ≤ Vt, and the Vt satisfy the following Lyapunov-type inequality: For constants α ∈ [0, 1) and β ≥ 0,
V1 ≤ β, (6) ∀t ≥ 2, Vt ≤ αVt−1 + Ut−1, (7) ∀t ≥ 1, E [Ut | Ft−1] ≤ β. (8)
Assumption 3.2 models a scenario where the gradient error dynamics is a combination of contracting and expanding behaviors. Contraction shrinks the error and is represented by a factor α. External noise, represented by the Ut terms, prevents the error from vanishing completely. Note that the assumption would be satisfied in the unbiased case by simply setting Vt = 0.
We can now state our result on the expected number of iterations for SGD with early stopping: Proposition 3.3. Let {xt}t≥1 be as in Algorithm 1. Let Assumptions 2.1, 2.2, 2.3, 3.1, and 3.2 hold. For > 0, let τ( ) be the stopping time τ( ) = inf{n ≥ 1 | n ≡ 1 (modm) and ‖∇fV (xn)‖2 ≤ }. Suppose that η ≤ 1L and − 2Lmησ2v − 2mβ/(1− α)−G2d1(µV , µT )2 > 0. Then
E[τ( )] ≤ ηG 2d1(µV , µT ) 2 + 2(fT (x1)− f∗) + η /m+ 2ηβ/(1− α) η /m− 2Lη2σ2v − 2ηβ/(1− α)− ηG2d1(µV , µT )2/m . (9)
Furthermore, it holds with probability 1 that ‖∇fT (xτ( ))‖2 ≤ 2 + 2G2d1(µV , µT )2.
This result can be strengthened by assuming a coupling between the step-size and the expansion bound β, then the result can be strengthened as demonstrated in the next corollary.
Corollary 3.4. Let Assumptions 2.1, 2.2, 2.3, 3.1, and 3.2 hold. In the context of Proposition 3.3, let the constant β be of the form β = ηR for some R ≥ 0, and suppose that > G2d1(µV , µT )2. Let c ∈ (0, 1) and let the step-size be
η = c ·min { 1
L , −G2d1(µV , µT )2 m(2Lσ2v + 2R/(1− α))
} . (10)
Then
E[τ( )] = O (
m2 (1 +R/(1− α)) (1− c) c ( −G2d1(µV , µT )2)2
) (11)
See the appendix for the complete formula, including lower order terms.
From this corollary we see that when β is proportional to the step-size η, the complexity bound is of the order O(1/ 2). As we shall see below, this coupling assumption is satisfied for Stacked SGD. For the case of using SGD to minimize a finite sum using unbiased, we obtain the following:
Corollary 3.5. Let Assumptions 2.1, 2.2, and 2.3 hold. Suppose each gradient estimate is obtained by selecting a data point yt ∈ YT uniformly at random and setting vt = ∇xf(yt, xt) and that there is a σ2v ≥ 0 such that ∀ y′ ∈ YT , x ∈ Rd, ∥∥∥∇xf(y′, x)− 1nT ∑y∈YT ∇xf(y, x)∥∥∥2 ≤ σ2v . If the step-sizes are defined according to Equation equation 20 using c = 1/2 then the expected number of IFO calls used by SGD before reaching an -approximate stationary point is
E [IFO ( )] = O (
mnV +m 2
( −G2d1(µV , µT )2)2 + nV ) See the appendix for the complete formula, including lower-order terms.
Note that when d(µV , µT ) is on the order of √ , this result states that the expected IFO complexity is O(1/( 2)). This can be compared with the RSG algorithm, where O(1/( 2)) iterations are sufficient for the expected (squared) norm of the gradient at a random iterate to be at most (Cor. 2.2 in Ghadimi & Lan (2013)).
4 STACKED SGD
Cluster 1
Cluster 2
Cluster 3
Cluster 4
In this section, we introduce SSGD for decentralized optimization, which involves a distributed system with two types of nodes: workers and communicators. The communication pattern is shown on the right for a hypothetical system consisting of 16 workers and 4 communicators. Workers (circles) are grouped into clusters, each containing a communicator (triangles).
Algorithm 2 depicts the steps of SSGD. For every m epochs, the gradient of the validation function is computed at the mean of the communication node parameters. The algorithm stops when the norm of this gradient drops below a threshold. This step is naturally carried out by one of the communication nodes because they can store the average computed during Line 9 of the algorithm. For workers, each iteration begins on Line 6 with one step of SGD. Then, on Line 7, an averaging step is performed to partially synchronize the model parameters within the local cluster. The steps for a communication node begin in Line 9 with an averaging step, among the other communication nodes from each cluster. On Line 10, there is a partial synchronization with the worker nodes in the cluster.
The key to the algorithm’s efficiency is that the averaging among communication nodes can happen in parallel with the gradient descent step as it occurs within the cluster. In the naive approach to parallelizing SGD, all nodes must block after each iteration for synchronization of their parameters.
We assume there are M ≥ 1 clusters, each containing K ≥ 1 computation nodes and 1 communication node. Thus, there are M(K + 1) total nodes. Given a worker node 1 ≤ i ≤ KM , we let
Algorithm 2 SSGD with early stopping 1: input: Node id i, initial parameters xi1 (for workers), or x̂ i 1 (for communicators)
2: t = 1 3: while ∥∥∥∇fV ( 1M ∑Mj=1 x̂jt)∥∥∥2 > do 4: for n = t to t+m− 1 do 5: if node i is a computation node then 6: xi
n+ 12 = xin − ηvin
7: xin+1 = 1
K+1
( x̂ c(i) n+ 12 + ∑K j=1 x j n+ 12 ) 8: else (node i is a communication node) 9: x̂i
n+ 12 = 1M
∑M j=1 x̂ j n
10: x̂in+1 = 1
K+1 ( x̂i n+ 12 + ∑ j∈c−1(i) x j n+ 12 ) 11: end if 12: end 13: t = t+m 14: end 15: return x̂t+ 12
4 8 16 32 64 128 256 Number of nodes
104
105
D at
a p
er se
co n
d
SGD
SSGD
4 8 16 32 64 128 256 Number of nodes
0
1
2
3
f T (x t)
SGD
SSGD
SSGD results for a climate modeling task. Top: throughput vs number of nodes.. Bottom: comparison of the loss of the algorithms.
c(i) ∈ {1, . . . ,M} denote the index of the communication node of the group containing worker i. For analysis, we define the filtration {Ft}t≥0 as follows:
F0 = σ ( { xi1 ∣∣ 1 ≤ i ≤ KM} ∪ {x̂i1 ∣∣ 1 ≤ i ≤M} ),
∀t ≥ 1, Ft = σ ( { xi1, v i n ∣∣ 1 ≤ n ≤ t, 1 ≤ i ≤ KM} ∪ {x̂i1 ∣∣ 1 ≤ i ≤M} ). We assume that the gradient estimates used in SSGD are unbiased and have bounded variance: Assumption 4.1. For any t ≥ 1 and 1 ≤ i ≤ KM,
E [ vit −∇fT (xit) | Ft−1 ] = 0, (12)
E [∥∥vit −∇fT (xit)∥∥2 | Ft−1] ≤ σ2v . (13)
The first step in our analysis is a bound on the dispersion of the parameters across the system. Proposition 4.2. Let Assumption 2.1, 2.2, and 4.1 hold, and let the variables x̂it be as defined in Line 6 of Algorithm 2. Suppose the step-size satisfies η < 1/(2LK). Define the variables V1, U1, V2, U2, . . . and the constants α, β as follows:
Vt = L2
M2 M∑ i=1 M∑ j=1 ∥∥∥x̂it − x̂jt∥∥∥2 , (14a) Ut = ηKL
(K + 5/8)
1
M2 M∑ i=1 M∑ j=1 ∥∥∥∥∥ ( 1 K ∑ k∈c−1(i) vkt −∇fT (x̂it) ) − ( 1 K ∑ k∈c−1(j) vkt −∇fT (x̂jt ) )∥∥∥∥∥ 2 (14b)
α = (K + 3/4)2
(K + 1)2 , (14c)
β = η · 4Lσ 2 v
(K + 5/8) . (14d)
Then, for all t ≥ 1, it holds that Vt+1 ≤ αVt + Ut and E[Ut | Ft−1] ≤ β.
In the preceding definitions, Vt represents the dispersion of the parameter values across different nodes. An averaging step tends to reduce the dispersion by a factor of α, while the independent gradient updates at each node may increase parameter dispersion by an amount β. This result allows us to model SSGD as a form of biased SGD, leading to the following: Proposition 4.3. Let Assumptions 2.1, 2.2, 2.3, and 4.1 hold. Assume that the initial parameters at every node are equal: xi1 = x j 1, for all 1 ≤ i, j ≤ KM and x̂i1 = x̂j1 for 1 ≤ i, j ≤ M . For some c < 12K , suppose that the step-size η is
η = c ·min { 1
L , −G2d1(µV , µT )2 m(2Lσ2v/K + 32Lσ 2 v(K + 1)/(K + 5/8))
} . (15)
Let x̂t be the average of the communicator states at time t: x̂t = 1M ∑M i=1 x̂ i t. For > 0, define τ( ) to be the stopping time τ( ) = inf{n ≥ 1 | n ≡ 1 (mod m) and ‖∇fV (x̂n)‖2 ≤ }. Then,
E[τ( )] = O (
m2 (1− c)c( −G2d1(µV , µT )2)2 ) .
Refer to the Appendix for the complete formula, including lower-order terms.
This leads to a bound on the expected IFO complexity to minimize a finite sum using SSGD: Corollary 4.4. Let Assumptions 2.1, 2.2, and 2.3 hold. Suppose each gradient estimate is obtained by selecting a data point yjt ∈ YT uniformly at random and setting vjt = ∇xf(yjt , xjt ) and that there is a σ2 ≥ 0, such that ∀ y′ ∈ YT , x ∈ Rd,
∥∥∥∇xf(y′, x)− 1nT ∑y∈YT ∇xf(y, x)∥∥∥2 ≤ σ2v . If the step-sizes are defined as in Equation equation 15 with c = 1/(4K), then the expected number of IFO calls used by SSGD before reaching an -approximate stationary point is
E [IFO( )] = O ( mK(nV +mK)
( −G2d1(µV , µT )2)2 + nV
) .
Experimental Result SSGD has been implemented to train a neural network model as part of research into spatio-temporal data analysis for climate research. The model is LSTNet, a neural network architecture that includes a convolutional neural network to extract short-term local dependency patterns from spatial variables and a recurrent neural network component to discover long-term patterns from time series trends Lai et al. (2018). LSTNet is trained for the prediction of solar radiation from past sensor measurements. Fig. 4 shows results from the experiment. The upper plot compares the throughput of SGD and SSGD as the number of nodes increases. We see that the performance of SGD degrades after 64 nodes, while SSGD maintains near-linear scalability. The lower plot shows that although SSGD involves biased gradients, it does not sacrifice accuracy, and yields similar prediction errors compared to SGD. See Appendix G.1 for more details about the experimental methodology.

heading: 5 DECENTRALIZED SGD
text: In this section we show how our methodology can be applied to Decentralized SGD (DSGD) which is another variant of distributed SGD. Recently, DSGD was analyzed using randomization Lian et al. (2017). In this section we complement that analysis by studying the expected running time of the algorithm.
The steps of DSGD are shown in Algorithm 3. The procedure involves M > 0 worker nodes that participate in the optimization. A communication matrix a describes the connectivity among the workers; ai,j > 0 means that workers i and j will communicate after each gradient descent step. At each step of optimization, every node computes a weighted average of the parameters in its local neighborhood, as determined by the connectivity matrix. This is combined with a local gradient approximation to obtain the new parameter at the worker. The data that is returned by the algorithm (assuming the termination criteria is met) is the average of the parameters throughout the system, denoted xt:
xt = 1
M M∑ i=1 xi (16)
Every m epochs, the norm of gradient of the validation function is evaluated at the average parameter and the algorithm terminates when this norm falls below a threshold.
Algorithm 3 DSGD with early stopping 1: input: Initial point x1 ∈ Rd, initial parameters xi1. 2: t = 1 3: while ‖∇fV (xt)‖2 > do 4: for n = t to t+m− 1 do 5: xin+1 = M∑ j=1 ai,jx j n − ηvin 6: end 7: t = t+m 8: end 9: return xt The intuitive justification for DSGD is that it may be more efficient compared to naive approaches to parallelizing SGD, since whenever ai,j = 0 then the nodes i and j need not communicate. In Lian et al. (2017) those authors offer theoretical support for the superiority of DSGD. In the present work, our goal is to analyze the expected running time of DSGD as an example of our the abstract theory developed above may be applied in practice. We leave comparisons of the algorithms for future work. For the analysis, we define the filtration {Ft}t≥0 as follows: F0 = σ ( { xi1 ∣∣ 1 ≤ i ≤M} ),
∀t ≥ 1, Ft = σ ( { xi1, v i n ∣∣ 1 ≤ n ≤ t, 1 ≤ i ≤M} ) The connectivity matrix a is subject to the same conditions as in Lian et al. (2017): Assumption 5.1. The M ×M connectivity matrix a is symmetric and stochastic. The spectral gap, denoted by ρ and defined as ρ = (max{|λ2(a)|, |λM (a)|})2 is assumed to satisfy ρ < 1.
To make the proofs clear and concise, we make the assumption the parameters at each node are single real-numbers. That is, throughout this section we assume d = 1 in Assumption ??.
We also assume that the gradient estimates used at each worker at unbiased and have bounded variance. Assumption 5.2. For any t ≥ 1 and 1 ≤ i ≤M,
E [ vit −∇fT (xit) | Ft−1 ] = 0, (17)
E [ |vit −∇fT (xit)|2 | Ft−1 ] ≤ σ2v . (18)
For the analysis, we show that the sequence of averages xt for t = 1, 2, . . . can be modeled as being generated by a biased version of SGD, using the tools from Section 3. This involves showing that the distance between local parameter values xit and the system average can be controlled, as shown in the following. Proposition 5.3. Let Assumptions 2.1, 2.2, 5.1, and 5.2 hold. Suppose the step-size satisfies η ≤ (1−√ρ)/(4 √ 2L). Define the variables V1, U1, V2, U2, . . . and the constants α, β as follows:
Vt = L2
M M∑ i=1 |xit − xt|2, (19a)
Ut = 8 η 2 L
2(1 + ρ) M(1− ρ) M∑ i=1 |vin −∇f(xit)|2, (19b)
α = (1 + ρ)
2ρ
(√ ρ+ 1
2
)2 (19c)
β = η L √ 2
1− ρσ 2 v (19d)
Then for all t ≥ 1 it holds that Vt+1 ≤ αVt + Ut and E[Ut | Ft−1] ≤ β.
Using this result on the dispersion of the parameters, we can move to the main result on decentralized SGD. The result gives conditions that guarantee the expected time E[τ( )] is finite, and also bounds this time in terms of the problem data. Notably, it shows a dependence on ρ, which is the mixing rate of the connectivity matrix. Proposition 5.4. Let Assumptions 2.1, 2.2, 5.1, 2.3 and 5.2 hold. Assume that the initial parameters at every node are equal: xi1 = x j 1 for all 1 ≤ i, j ≤M . Let c ≤ 1−√ρ 4 √ 2 and define R = L √ 2σ2v/(1− ρ)
Algorithm 4 SVRG with early stopping 1: input: Initial point x1m ∈ Rd 2: for s = 1, 2, . . . do 3: xs+10 = x s m
4: gs+1 = 1nT ∑ y∈YT ∇fT (y, x s+1 0 ) 5: if ‖gs+1‖2 ≤ then return xs+10 6: for t = 0 to m− 1 do 7: Sample yst uniformly at random from YT 8: vst = ∇f(yst , xs+1t )−∇f(yst , xs+10 ) + gs+1 9: xs+1t+1 = x s+1 t − ηvst
10: end 11: end
10−6 10−5 10−4 10−3 10−2 10−1 100 ‖∇fT‖2 0.0
2.5
5.0
IF O
C al
ls
×105
SGD
SVRG
10−6 10−5 10−4 10−3 10−2 10−1 100 ‖∇fT‖2 0.0
2.5
5.0
IF O
C al
ls
×105
SGD
SVRG
IFO calls to approximate stationarity for MNIST (top) and CIFAR-10 (bot.).
and let α be as in Equation equation 19c. Let the step-size be η = c ·min { 1
L , −G2d1(µV , µT )2 m(2Lσ2v + 2R/(1− α))
} . (20)
For > 0 define τ( ) to be the first time the norm of the gradient of the validation function falls below ; That is, τ( ) = inf{n ≥ 1 | n ≡ 1 (mod m) and ‖∇fV (xn)‖2 ≤ }. Then
E[τ( )] = O (
m2 (1 +R/(1− α)) (1− c) c ( −G2d1(µV , µT )2)2
) (21)
Note that in the above result, the order of the convergence is the same as for regular SGD and stacked SGD. An interesting avenue for future work would be to explore whether it is possible to obtain bounds where the step-size condition does not depend on the epoch length m.

heading: 6 SVRG
text: We demonstrate how the method can applied to a variant of the SVRG Johnson & Zhang (2013) with early stopping, shown in Algorithm 4. Each epoch begins with a full gradient computation (Line 4), and then an inner loop runs for m steps. The first step of the inner loop is to choose a random data point (Line 7). Then, the update direction is computed (Line 8) and used obtain the next parameter (Line 9).
Combining some existing bounds for SVRG with our stopping time approach yields the following bound on the expected number of iterations until SVRG with early stopping terminates: Proposition 6.1. Let Assumptions 2.1 and 2.2 hold and consider the variables xs+1t defined by Algorithm 4. Let ξ = 1/4 and suppose that the step-size is set to η = ξ/(Ln2/3T ) and the epoch length is m = bnT /(3ξ)c. For > 0, define τ( ) to be the stopping time τ( ) = inf { s ≥ 1
∣∣ ‖∇fT (xs+10 )‖2 ≤ } . Then, E[τ( )] ≤ 1 + (40Ln2/3T (fT (x1m)− f∗))/ . Note that Proposition 6.1 counts the number of epochs until an approximate stationary point is generated. A bound on the number of IFO calls can be obtained by multiplying τ by the number of IFO calls per epoch, which is nT + 2m. This immediately leads to the following result: Corollary 6.2. Let Assumptions 2.1 and 2.2 hold and suppose the step-size η and epoch length m are defined as in Proposition 6.1. Then, the expected number of IFO calls until SVRG returns an approximate stationary point is E [IFO ( )] = O((n5/3T / ) + nT ).
This result may be compared with Cor. 4 of Reddi et al. (2016a), which concerns an upper bound on the IFO calls for the expected (squared) norm of the gradient at a randomly selected iterate to be less than . Our result concerns the expected number of IFO calls before the algorithm terminates with an iterate that is guaranteed to be an approximate stationary point with probability 1, a stronger property.
Fig. 6 illustrates the expected IFO complexity of SVRG and SGD. The top plot shows the results of an experiment using a simple logistic classifier on the MNIST dataset, and the bottom shows the result of training a one-layer neural net on the CIFAR-10 dataset. The error bands represent the standard deviation of the measurements over five independent runs. Appendix G.2 includes more details about the experimental methodology. In each case, SVRG is better at obtaining accurate solutions, while, for SGD, the expected IFO calls seem to become unbounded for sufficiently small .

heading: 7 GENERALIZATION PROPERTIES
text: Typically, the training and validation sets are made from independent samples of a test distribution µ, and it is of interest to estimate the model performance on samples from this test distribution. Define fµ : Rd → R as f(x) = Eµ[f(y, x)]. In this section, we derive a bound on E[‖∇fµ(xτ( ))‖2], where the expectation is not only over the variates generated by optimization, but also over the random choice of the datasets YV and YT . We will show how Wasserstein concentration bounds, which concern the average distance between µ and its empirical versions, can be used for this task. Consider the following from Dereich et al. (2013). Theorem 7.1 (Dereich et al. (2013), Special case of Theorem 1). Let d ≥ 3 and let µ be a measure on Rd, such that J = Eµ [ ‖y‖3 ]1/3 <∞. Then, there is a constant κd, such that E[d2(µ, µV )2] ≤ κdJn −3/d V .
This leads to a bound on the generalization performance of the iterates returned by SGD with early stopping, whose proof is deferred to the Appendix.
Corollary 7.2. Let the conditions of Proposition 3.3 hold. Further assume that J = Eµ [ ‖y‖3 ]1/3 < ∞, the validation set YV is an empirical version of µ, and y 7→ ∇xf(y, x) is uniformly G-Lipschitz. If xτ ( ) is the output of Algorithm 1, then E[‖∇fµ(xτ( ))‖2] ≤ 2 + 2G2κdJn−3/dV . Note that for strongly convex optimization, it is possible to obtain rates of convergence of the test error that are independent of the dimension Hsu & Sabato (2016). Corollary 7.2 is interesting as it accounts for data distribution properties (via the 3rd moment J) and does not depend on the number of iterations used in SGD. This result could be compared with Hardt et al. (2016), where the authors proved a bound on the generalization gap for function values in terms of the number of iterations T and samples in the training set nT . There, the bound is independent of d but increasing with T , while our bound is independent of the number of iterations in SGD. It will be interesting to determine if these two analyses can be combined.

heading: 8 DISCUSSION
text: This work presents an analysis of several stochastic gradient methods that use early stopping based on a validation function. We demonstrated that by blending existing analysis techniques with some basic tools related to stopping times, it is possible to bound the expected number of iterations and gradient evaluations to generate approximate stationary points. We also considered decentralized optimization and introduced a new algorithm, SSGD, that proved amenable to analysis in our framework. For SSGD, we obtained a convergence rate using our results for biased SGD, and experiments showed that the algorithm has favorable scaling properties compared to basic parallel SGD. Our application to SVRG demonstrated that the theoretical approach can be applied in various settings. We also considered the generalization properties of the output of early stopping. We hope these efforts inspire other works that investigate the theoretical and practical aspects of early stopping.

heading: APPENDIX: ON THE EXPECTED RUNNING TIME OF NONCONVEX OPTIMIZATION WITH EARLY STOPPING

heading: A PRELIMINARIES
text: Our analyses make use of a quadratic bound for the training function which follows from Assumption 2.1:
∀x, v ∈ Rn, fT (x+ v) ≤ fT (x) +∇fT (x)T v + L
2 ‖v‖2. (22)

heading: A.1 STOCHASTIC PROCESSES
text: The formal setting of a stochastic optimization algorithm involves a probability space (Ω,F ,P), consisting of a sample space Ω, a σ-algebra F of subsets of Ω and a probability measure P on the subsets of Ω that are in F . The algorithm takes an initial point x1 and defines a sequence of random variables {xt(ω)}t>1. Intuitively Ω represents the random data used by the algorithm, such as indices used to define mini-batches. For ease of notation we will omit the dependence of random variates in the algorithms on ω ∈ Ω. A filtration {Ft}t=0,1,... is an increasing sequence of σ-algebras, with the interpretation that Ft represents the information available to an algorithm up to and including time t. A random variable x : Ω→ Rd is said to be Ft measurable if it can be expressed in terms of the state of the algorithm up and including time t. A rule for stopping an algorithm is represented as a stopping time, which is a random variable τ : Ω→ {0, 1, . . . ,∞} with the property that the decision of whether to stop or continue at time n is made based on the information up to and including time n.
The following proposition will be used through out our analysis of the different algorithms.
Proposition A.1. Let τ be a stopping time with respect to a filtration {Ft}t=0,1,.... Suppose there is a number c < ∞ such that τ ≤ c with probability one. Let x1, x2, . . . be any sequence of random variables such that each xt is Ft-measurable and E[‖xt‖] <∞. Then
E [ τ∑ t=1 xt ] = E [ τ∑ t=1 E [xt | Ft−1] ] . (23)
Proof. This is a consequence of the optional stopping theorem (Theorem 10.10 in Williams (1991)). Define S0 = 0 and for t ≥ 1, let St = t∑ i=1
(xi − E[xi | Fi−1]). Then S0, S1, . . . is a martingale with respect to the filtration {Ft}t=0,1,..., and the optional stopping theorem implies E[Sτ ] = E[S0]. But E[S0] = 0, and therefore E[Sτ ] = 0, which is equivalent to Equation equation 23.

heading: B ANALYSIS OF BIASED SGD

heading: PROOF OF PROPOSITION 3.3
text: Proof. For convenience, define the random variables δt for t = 1, 2, . . . as δt = vt−∇fT (xt). From equation 22, it holds that
fT (xt+1) ≤ fT (xt)− η∇fT (xt)T (∇fT (xt) + δt + ∆t) + L
2 η2‖∇fT (xt) + δt + ∆t‖2.
Summing this over t = 1, . . . ,m yields
fT (xm+1) ≤ fT (x1)− m∑ t=1 η∇fT (xt)T (∇fT (xt) + δt + ∆t)
+ m∑ t=1 L 2 η2‖∇fT (xt) + δt + ∆t‖2
= fT (x1)− m∑ t=1 η ( 1− L 2 η ) ‖∇fT (xt)‖2 − m∑ t=1 η(1− Lη)∇fT (xt)T δt
+ m∑ t=1 L 2 η2‖δt‖2 − m∑ t=1 η(1− Lη)∇fT (xt)T∆t + m∑ t=1 L 2 η2‖∆t‖2
+ m∑ t=1 Lη2δTt ∆t.
(24)
Note that, in general, for any numbers a, b it is the case that |ab| ≤ 12a2 + 12b2. Then
|δTt ∆t| ≤ ‖δt‖‖∆t‖ ≤ 1
2 ‖δt‖2 +
1 2 ‖∆t‖2 (25)
and
|∇fT (xt)T∆t| ≤ ‖∇fT (xt)‖‖∆t‖ ≤ 1
2 ‖∇fT (xt)‖2 +
1 2 ‖∆t‖2. (26)
Combining Equations equation 24, equation 25, and equation 26, we obtain
fT (xm+1) ≤ f(x1)− m∑ t=1 η ( 1− L 2 η − 1 2 (1− Lη) ) ‖∇fT (xt)‖2
− m∑ t=1 η(1− Lη)∇fT (xt)T δt + N∑ t=1 ( L 2 η2 + 1 2 Lη2 ) ‖δt‖2
+ N∑ t=1 ( L 2 η2 + 1 2 η(1− Lη) + Lη2 1 2 ) ‖∆t‖2
= f(x1)− m∑ t=1 η 1 2 ‖∇fT (xt)‖2 − m∑ t=1 η(1− Lη)∇fT (xt)T δt + m∑ t=1 Lη2‖δt‖2
+ m∑ t=1 1 2 η(1 + Lη)‖∆t‖2.
Rearranging terms and noting that fT (xm+1) ≥ f∗, this yields m∑ t=1 η 1 2 ‖∇fT (xt)‖2 ≤ fT (x1)− f∗ − m∑ t=1 η(1− Lη)∇fT (xt)T δt + m∑ t=1 Lη2‖δt‖2
+ m∑ t=1 1 2 η(1 + Lη)‖∆t‖2
≤ f(x1)− f∗ − m∑ t=1 η(1− Lη)∇fT (xt)T δt + m∑ t=1 Lη2‖δt‖2
+ m∑ t=1 ηVt.
(27)
where in the second inequality we used the assumptions that η ≤ 1/L and ‖∆t‖2 ≤ Vt. Based on our assumption that ‖∇fV (x)−∇fT (x)‖ ≤ Gd1(µV , µT ), it follows that
‖∇fV (x)‖2 ≤ 2G2d1(µV , µT )2 + 2‖∇fT (x)‖2. (28)
Also note that n∑ t=1 1t≡1 (mod m) = ⌈ n m ⌉ ≤ n m + 1. (29)
Combining equation 28 and equation 29 results in
τ( )∧n∑ t=1 1t≡1 (mod m)‖∇fV (xt)‖2 ≤ τ( )∧n∑ t=1 1t≡1 (mod m)G 2d1(µV , µT ) 2
+ τ( )∧n∑ t=1 1t≡1 (mod m)‖∇fT (xt)‖2
≤ G2d1(µV , µT )2 (
(τ( ) ∧ n) m + 1
) + τ( )∧n∑ t=1 ‖∇fT (xt)‖2.
(30)
For each n ≥ 1 define τ( )∧n to be the stopping time which is the minimum of τ( ) and the constant n. Applying Proposition A.1 and Assumption 4, it holds that
E τ( )∧n∑ t=1 ∇fT (xt)T δt  = 0 (31) and using Proposition A.1 with Assumption 5 gives
E τ( )∧n∑ t=1 ‖δt‖2  ≤ σ2vE[τ( ) ∧ n]. (32)
Next, according to conditions equation 6, and equation 7, it holds for any m ≥ 1 and with probability one that
m∑ t=1 Vt ≤ α m∑ t=1 Vt + m∑ t=1 Ut + β (33)
and by equation 8 together with Proposition A.1,
E τ( )∧n∑ t=1 Ut  ≤ E[τ( ) ∧ n]β. (34) Combining equation 33 and equation 34, then
E τ( )∧n∑ t=1 Vt  ≤ αE τ( )∧n∑ t=1 Vt + (E[τ( ∧ n)] + 1)β which, upon rearranging, results in
E τ( )∧n∑ t=1 Vt  ≤ (E[τ( ) ∧ n] + 1) β 1− α. (35)
Furthermore, by definition of τ ,
E τ( )∧n∑ t=1 1t≡1 (mod m)‖∇fV (xt)‖2  ≥ E (τ( )∧n)−1∑ t=1 1t≡1 (mod m)‖∇fV (xt)‖2 
≥ E[(τ( ) ∧ n)− 1] m
(36)
Here we used that n−1∑ t=1 1t≡1 (mod m) = ⌈ n− 1 m ⌉ ≥ n− 1 m
Combining equation 27, equation 30, equation 31, equation 32, equation 35 and equation 36 results in η
2m (E[τ( ) ∧ n]− 1) ≤η 2 G2d1(µV , µT ) 2
( E[τ( ) ∧ n]
m + 1
) + fT (x1)− f∗
+ Lη2σ2vE [τ ( ) ∧ n] + η β
1− α (E[τ( ) ∧ n] + 1) . This can be rearranged into( η
m − 2Lη2σ2v −
2ηβ 1− α − η m G2d1(µV , µT ) 2
) E[τ( ) ∧ n] ≤ ηG2d1(µV , µT )2
+ 2(fT (x1)− f∗) + 2η β 1− α + η m .
which in turn is equivalent to
E[τ( ) ∧ n] ≤ ηG 2d1(µV , µT ) 2 + 2(fT (x1)− f∗) + η /m+ 2ηβ/(1− α) η /m− 2Lη2σ2v − 2ηβ/(1− α)− ηG2d1(µV , µT )2/m . (37)
Note that the sequence of random variables {(τ( ) ∧ n)}n=1,2,... is monotone increasing, and converges pointwise to τ( ). Then the claimed relation equation 9 follows from equation 37 by the monotone convergence theorem.
Finally, by applying equation 28 with the roles of fV and fT reversed, and using the definition of τ( ), it follows that
‖∇fT (xτ( ))‖ ≤ 2‖∇fV (xτ( ))‖2 + 2G2d1(µV , µT )2
≤ 2 + 2G2d1(µV , µT )2.

heading: PROOF OF COROLLARY 3.4
text: Proof. According to the assumption on the step-size η (Equation equation 20), η [ ( −G2d1(µV , µT )2)/m− η(2Lσ2v + 2R/(1− α)) ] ≥ η(1− c)( −G2d1(µV , µT )2)/m (38) and
1 η ≤ L c + m(2Lσ2v + 2R/(1− α)) c ( −G2d1(µV , µT )2) . (39)
Combining these inequalities with the conclusion of Proposition 3.3 (relation equation 9) yields
E[τ( )] A ≤ ηG
2d1(µV , µT ) 2 + 2(fT (x1)− f∗) + η /m+ 2ηβ/(1− α) η( −G2d1(µV , µT )2)/m− 2Lη2σ2v − 2ηβ/(1− α) .
B = ηG2d1(µV , µT ) 2 + 2(fT (x1)− f∗) + η /m+ 2η2R/(1− α) η( −G2d1(µV , µT )2)/m− 2Lη2σ2v − 2η2R/(1− α) .
C ≤ ηG
2d1(µV , µT ) 2 + 2(fT (x1)− f∗) + η /m+ 2η2R/(1− α)
η(1− c) ( −G2d1(µV , µT )2)/m .
(40)
Step A was established by Proposition 3.3. Step B uses the assumption that β = ηR. Step C is an application of equation 38. Next, we will upper-bound the final inequality in three steps. First, using Inequality equation 39, we see that
2(fT (x1)− f∗) η(1− c) ( −G2d1(µV , µT )2)/m ≤ 2(fT (x1)− f ∗) (1− c) ( −G2d1(µV , µT )2)/m)
× ( L
c + m(2Lσ2v + 2R/(1− α)) c ( −G2d1(µV , µT )2) ) =
2m(fT (x1)− f∗) (1− c) c ( −G2d1(µV , µT )2) L
+ 2m2(fT (x1)− f∗) (1− c) c ( −G2d1(µV , µT )2)2 ( 2Lσ2v + 2R/(1− α) ) .
(41)
Next, 2η2R/(1− α)
η(1− c) ( −G2d1(µV , µT )2)/m = 2ηR/(1− α) (1− c) ( −G2d1(µV , µT )2)/m
≤ 2cR/(1− α) (1− c) ( −G2d1(µV , µT )2)/m × −G 2d1(µV , µT ) 2 m(2Lσ2v + 2R/(1− α))
= c (1− c) × 2R/(1− α)
(2Lσ2v + 2R/(1− α)) ≤ c
(1− c) . (42)
Finally, ηG2d1(µV , µT ) 2 + η /m
η(1− c) ( −G2d1(µV , µT )2)/m =
G2d1(µV , µT ) 2 + /m
(1− c) ( −G2d1(µV , µT )2)/m
= mcG2d1(µV , µT ) 2 + c
(1− c)c( −G2d1(µV , µT )2) .
(43)
Combining equation 40 with equation 41, equation 42 and equation 43, we find that
E[τ( )] ≤ 4m 2(fT (x1)− f∗)
( Lσ2v +R/(1− α) ) (1− c) c ( −G2d1(µV , µT )2)2
+ 2Lm(fT (x1)− f∗) +mcG2d1(µV , µT )2 + c
(1− c) c ( −G2d1(µV , µT )2) +
c 1− c . (44)

heading: PROOF OF COROLLARY 3.5
text: Proof. If the algorithm runs until iteration τ( ), then the number of times that the full gradient of fV is calculated is dτ( )/me ≤ τ( )/m + 1, and the number of IFO calls for the training function is τ( )− 1. Therefore
IFO( ) ≤ ( τ( )
m + 1
) nV + (τ( )− 1) ≤ τ( ) (nV m + 1 ) + nV . (45)
Note that under our assumption on the gradient estimates vt, we are in the unbiased setting where R = 0. Combining equation 44 with R = 0 and equation 45, we obtain
E[IFO( )] ≤ (
4m2(fT (x1)− f∗)Lσ2v (1− c) c ( −G2d1(µV , µT )2)2
+ 2Lm(f(x1)− f∗) +mcG2d1(µV , µT )2 + c
(1− c) c ( −G2d1(µV , µT )2) +
c
1− c ) × (nV m + 1 ) + nV .
(46)
Using c = 1/2 and neglecting lower-order terms, then, E[IFO( )] = O ( mnV +m 2
( −G2d1(µV , µT )2)2 + nV
) .

heading: C ANALYSIS OF STACKED SGD

heading: PROOF OF PROPOSITION 4.2
text: Proof. Note that the variables inside the stacked algorithm satisfy the following identities: For all 1 ≤ i ≤MK,
x̂ c(i) t = x i t, (47)
and for all pairs communication nodes nodes 1 ≤ j ≤M , and 1 ≤ k ≤M , and t ≥ 1, x̂j t+ 12 = x̂kt+ 12 . (48) Let 1 ≤ i ≤M and 1 ≤ j ≤M be arbitrary indices of communication nodes. By the definitions in Algorithm 2, we can see that for any t ≥ 1,
‖x̂it+1 − x̂jt+1‖
A = ∥∥∥∥∥ 1K + 1  ∑ k∈c−1(i) xkt+ 12 − 1 K + 1  ∑ k∈c−1(j) xkt+ 12 ∥∥∥∥∥ B =
∥∥∥∥∥ 1K + 1  ∑ k∈c−1(i) ( xkt − ηvkt )− 1 K + 1  ∑ k∈c−1(j) ( xkt − ηvkt )∥∥∥∥∥ C =
∥∥∥∥∥ 1K + 1  ∑ k∈c−1(i) ( x̂it − ηvkt )− 1 K + 1  ∑ k∈c−1(j) ( x̂jt − ηvkt )∥∥∥∥∥ D =
∥∥∥∥∥ KK + 1 (x̂it − x̂jt)− 1K + 1 ∑ k∈c−1(i) ηvkt + 1 K + 1 ∑ k∈c−1(j) ηvkt ∥∥∥∥∥ E ≤ K K + 1 ∥∥∥x̂it − x̂jt∥∥∥+ η ∥∥∥∥∥∥ 1K + 1 ∑ k∈c−1(i) vkt − 1 K + 1 ∑ k∈c−1(j) vkt ∥∥∥∥∥∥ .
(49)
Step A follows from the definition of x̂it+1, x̂ j t+1 in Algorithm 2 and Equation equation 48. Step B follows from the definition of xk t+ 12
in Algo. 2. Step C follows from Equation equation 47. Step D follows from rearranging terms in the previous step, and noting that there are K workers in a cluster. Finally, Step E is simply the triangle inequality. For the second term on the right of the final inequality above,∥∥∥∥∥∥ 1K + 1 ∑
k∈c−1(i)
vkt − 1
K + 1 ∑ k∈c−1(j) vkt ∥∥∥∥∥∥ = K
K + 1 ∥∥∥∥∥∇f(x̂it)−∇f(x̂jt ) +  1 K ∑ k∈c−1(i) vkt −∇fT (x̂it) −  1 K ∑ k∈c−1(j) vkt −∇fT (x̂jt )
∥∥∥∥∥ ≤ KL K + 1 ‖x̂it − x̂jt‖
+ K
K + 1 ∥∥∥∥∥  1 K ∑ k∈c−1(i) vkt −∇fT (x̂it) −  1 K ∑ k∈c−1(j) vkt −∇fT (x̂jt ) ∥∥∥∥∥.
(50)
where the last inequality uses the Lipschitz gradient property (Assumption 2.1) and the triangle inequality. Combining equation 49 and equation 50 yields∥∥∥x̂it+1 − x̂jt+1∥∥∥ ≤ KK + 1(1 + Lη)∥∥∥x̂it − x̂jt∥∥∥
+ ηK
K + 1 ∥∥∥∥∥ ( 1 K ∑ k∈c−1(i) vkt −∇fT (x̂it) ) −  1 K ∑ k∈c−1(j) vkt −∇fT (x̂jt ) ∥∥∥∥∥. (51)
Note that for any k1 > 0 and all a, b we have |a+ b|2 ≤ (1 + k1)a2 + ( 1 + 1
k1
) b2. (52)
Combining equation 51 and equation 52 while using the assumption that η ≤ 1/(2LK) results in
‖x̂it+1 − x̂jt+1‖2 ≤ (1 + k1)K
2
(K + 1)2
( 1 + 1
2K )2 ∥∥∥x̂it − x̂jt∥∥∥2 + η (1 + k1)K
k1(K + 1)2 1 2L ∥∥∥∥∥  1 K ∑ k∈c−1(i) vkt −∇fT (x̂it) −  1 K ∑ k∈c−1(j) vkt −∇fT (x̂jt ) ∥∥∥∥∥ 2 . (53)
Let k1 = (K+3/4)2 (K+1/2)2 − 1. Then by the definition of α,
(1 + k1)K 2
(K + 1)2
( 1 + 1
2K
)2 = (1 + k1)(K + 1/2) 2
(K + 1)2 =
(K + 3/4)2 (K + 1/2)2 = α. (54)
Furthermore, 1 + k1 k1 ≤ (K + 1) 2 K/2 + 5/16 so
η (1 + k1)K k1(K + 1)2 1 2L ≤ η (K + 1)
2
(K/2 + 5/16)
K (K + 1)2 1 2L
= η K
(K/2 + 5/16)
1
2L
= η K
L(K + 5/8) .
(55)
Multiplying each side of equation 53 by L2/M2, summing the resulting inequality over i = 1, . . . ,M and j = 1, . . . ,M , and using the relations equation 54, equation 55 we see that for all t ≥ 1,
Vt+1 ≤ αVt + Ut. It remains to confirm that E[Ut | Ft−1] for all t ≥ 1.
Taking expectations in equation 14b, while noting that x̂c(k)t = x k t and applying the variance bound equation 13 along with the inequality ‖a+ b‖2 ≤ 2 ( ‖a‖2 + ‖b‖2 ) it holds for all t ≥ 1 that
E [Ut | Ft−1]
≤ η · KL (K + 5/8)
× 1 M2 M∑ i=1 M∑ j=1 E [∥∥∥∥∥  1 K ∑ k∈c−1(i) vkt −∇fT (x̂it) −  1 K ∑ k∈c−1(j) vkt −∇fT (x̂jt ) ∥∥∥∥∥ 2∣∣∣∣∣Ft−1 ]
= η KL
(K + 5/8)
1
M2 M∑ i=1 M∑ j=1 2 ( σ2v K + σ2v K ) = β.

heading: PROOF OF PROPOSITION 4.3
text: Proof. Note that according to Line 9 of SSGD,
x̂t = x̂ 1 t+ 12 = . . . = x̂Mt+ 12 (56)
and Lines 7 and 10 mean that for all j with c(j) = i,
x̂it = x j t (57)
Using these equations, together with the definitions in the SSGD algorithm, we obtain that
x̂t+1 A =
1
M M∑ i=1 1 K + 1 x̂it+ 12 + ∑ j∈c−1(i) xj t+ 12  B =  1 K + 1 x̂t + 1 M M∑ i=1 1 K + 1 ∑ j∈c−1(i) ( xjt − ηvjt
) C =  1 K + 1 x̂t + 1 M M∑ i=1 1 K + 1 ∑ j∈c−1(i) (x̂it − ηvjt )
 D = ( 1
K + 1 x̂t +
K
K + 1 x̂t −
1
M M∑ i=1 1 K + 1 ∑ j∈c−1(i) ηvjt
)
E = x̂t − η  1 M M∑ i=1 1 K + 1 ∑ j∈c−1(i) vjt  . For Step A, note that x̂t is the average of the x̂it, and then use definition of the x̂it from Line 9 of the SSGD algorithm. Step B follows from equation 56 and the definition of the xi
t+ 12 from Line 5 of
SSGD. Step C follows from equation 57. Step D follows by rearranging terms in the previous step, and again noting the definition of x̂t is the average of the x̂it. Step E follows by grouping terms.
Continuing, we can express this recursion as an approximate form of gradient descent:
x̂t+1 = x̂t − η  1 M M∑ i=1 1 K + 1 ∑ j∈c−1(i) vjt  = x̂t − η K
K + 1  1 M M∑ i=1 1 K ∑ j∈c−1(i) vjt  = x̂t − ηK
K + 1 (vt + ∆t) ,
where vt and ∆t are
vt = ∇fT (x̂t) + 1
M M∑ i=1  1 K  ∑ j∈c−1(i) vjt −∇fT (x̂it)  (58)
and
∆t = 1
M ( M∑ i=1 ∇fT (x̂it) ) −∇fT (x̂t). (59)
Based on the definition of vt in Equation equation 58 and on Assumption 4.1, for all t ≥ 1 it holds that
E [vt −∇fT (x̂t) | Ft−1] = 0,
E [ ‖vt −∇fT (x̂t)‖2 | Ft−1 ] ≤ σ 2 v
K .
Thus Assumption 3.1 is confirmed. Next we consider the Lyapunov condition of Assumption 3.2. Let the variables Ut, Vt for t ≥ 1 and the constants α, β be defined as in equation 14a-equation 14d.
Then by Assumption 2.1, ‖∆t‖2 = ∥∥∥∥∥ 1M M∑ i=1 ∇fT (x̂it)−∇f(x̂t) ∥∥∥∥∥ 2
≤ L2 ∥∥∥∥∥ 1M M∑ i=1 ( x̂it − x̂t )∥∥∥∥∥ 2
= L2 ∥∥∥∥∥∥ 1M M∑ i=1 x̂it − 1M M∑ j=1 x̂jt ∥∥∥∥∥∥ 2
= L2 ∥∥∥∥∥∥ 1M2 M∑ i=1 M∑ j=1 ( x̂it − x̂jt )∥∥∥∥∥∥ 2
≤ L2 1 M2 M∑ i=1 M∑ j=1 ∥∥∥x̂it − x̂jt∥∥∥2 = Vt.
The second line uses the Lipschitz gradient property, and the second to last line follows from Jensen’s inequality (Section 6.6, Williams (1991)).
By our assumption that each node has the same initial state, then V1 = 0, hence Inequality equation 6 holds. The Inequalities equation 7 and equation 8 are established by Prop 4.1. According to the definition of β (Equation equation 14d, we may write β = ηR where R = 4Lσ2v/(K + 5/8).
Note that 1/(1− α) = 2(K + 1)2/(K + 5/8). Therefore, R
1− α = 8Lσ 2 v(K + 1) 2/(K + 5/8)2 ≤ 32Lσ2v (60)
According to Corollary 3.4, then, a step-size of η = c ·min { 1
L , −G2d1(µV , µT )2 m(2Lσ2v/K + 16Lσ 2 v(K + 1) 2/(K + 5/8)2) } leads to
E[τ( )] ≤ 4m 2(fT (x1)− f∗)
( Lσ2v/K + 32Lσ 2 v ) (1− c) c ( −G2d1(µV , µT )2)2
+ 2Lm(fT (x1)− f∗) +mcG2d1(µV , µT )2 + c
(1− c) c ( −G2d1(µV , µT )2) +
c 1− c . (61)
Dropping the lower order terms, we see that E[τ( )] = O (
m2 (1− c)c( −G2d1(µV , µT )2)2 ) . (62)

heading: PROOF OF COROLLARY 4.4
text: Proof. If SSGD runs until iteration τ( ), then number of times that the full gradient of fV is calculated is dτ( )/me ≤ τ( )/m + 1, and the number of IFO calls for the training function is (τ( )− 1)K. Therefore
IFO( ) ≤ ( τ( )
m + 1
) nV + (τ( )− 1)K ≤ τ( ) (nV m +K ) + nV . (63)
Next, note that (1− c)c = (1− 14K ) 14K = 4K−1(4K)2 , which implies 1
(1− c)c = 16K2 4K − 1 ≤ 16K 3 . (64)
Combining equation 62, equation 63, and equation 64 we see that E [IFO( )] = O ( mK(nV +mK)
( −G2d1(µV , µT )2)2 + nV
) .

heading: D ANALYSIS OF DECENTRALIZED SGD
text: The following result is a restatement of Lemma 5 of Lian et al. (2017). Lemma D.1. Under Assumption 5.1, the matrix limk→∞ ak is well-defined and has entries a∞i,j = 1 M for 1 ≤ i, j ≤M . Furthermore, for all k ≥ 1, bound on the spectral gap implies ‖a∞ − ak‖2 ≤ ρk.

heading: PROOF OF PROPOSITION 5.3
text: Proof. For ease of notation, for each t ≥ 1 let yt and zt be the M -dimensional vectors with components defined as
yit = x i t − xt, (65)
zit = v i t −
1
M M∑ j=1 vjt . (66)
Let a∞ be the M ×M matrix with entries a∞i,j = 1M (see Lemma D.1). Then according to Line 5 of Algorithm 3, the yt satisfy the recursion yt+1 = (a− a∞)yt + ηzt. (67) Note that zt can be expressed as
zit = ∇f(xit)−∇f(xt)
+ vit −∇f(xit)− 1
M M∑ j=1 (vjt −∇f(xjt )) + 1 M M∑ j=1 (∇f(xjt )−∇f(xt)) (68)
Using the Lipschitz property (Assumption 2.1 ) then,
|zin| ≤ L|xin − xn|+ |vin −∇f(xin)|+ 1
M M∑ j=1 |vjn −∇f(xjn)|+ L M M∑ j=1 |xjn − xn| (69)
Squaring and summing Equation equation 69, M∑ i=1 |zin|2 = M∑ i=1 L|xin − xn|+ |vin −∇f(xin)|+ 1M M∑ j=1 |vjn −∇f(xjn)|+ L M M∑ j=1 |xjn − xn| 2
= L24 M∑ i=1 |xin − xn|2 + 4 M∑ i=1 |vin −∇f(xin)|2
+ 4
M M∑ i=1 M∑ j=1 |vjn −∇f(xjn)|2 + 4L2 M M∑ i=1 M∑ j=1 |xjn − xn|2
= L28‖yn‖2 + 8 M∑ i=1 |vin −∇f(xin)|2
Taking square roots on each sides of this equation yields
‖zn‖ ≤ √√√√L28‖yn‖2 + 8 M∑ i=1 |vin −∇f(xin)|2
≤ √ L28‖yn‖+ √√√√8 M∑ i=1 |vin −∇f(xin)|2 (70)
Combining equation 67 and equation 70, then,
‖yn+1‖ ≤ ( ‖a− a∞‖+ ηL √ 8 ) ‖yn‖+ η √√√√8 M∑ i=1 |vin −∇f(xin)|2
≤ (√ ρ+ ηL √ 8 ) ‖yn‖+ η √√√√8 M∑ i=1 |vin −∇f(xin)|2
In the second step we have applied Assumption equation 5.1 and Lemma D.1. Squaring this equation, for any k1 > 0 it holds that
‖yn+1‖2 ≤ (1 + k1) (√ ρ+ ηL √ 8 )2 ‖yn‖2 + η2 ( 1 + 1
k1
) 8 M∑ i=1 |vin −∇f(xin)|2. (71)
Let k1 = 1−ρ2ρ (in which case 1 + 1 k1 = 1+ρ1−ρ ). Multiplying each side of 71 by L 2/M and noting that Vt = L2 M ‖yt‖2, it follows that
Vt+1 ≤ (1 + ρ)
2ρ
(√ ρ+ ηL √ 8 )2 Vt + Ut (72)
It follows from the variance bound in Assumption 5.2 that
E [Ut | Ft−1] ≤ 8 η2 L2 1 + ρ
1− ρσ 2 v (73)
Using the assumption that η ≤ 1− √ ρ
4L √ 2 , then equation 72 and equation 73 become, respectively,
Vt+1 ≤ (1 + ρ)
2ρ
(√ ρ+ 1
2
)2 Vt + Ut
and
E [Ut | Ft−1] ≤ √ 2 η(1−√ρ)L1 + ρ 1− ρσ 2 v
≤ η L √ 2
1− ρσ 2 v

heading: PROOF OF PROPOSITION 5.4
text: Proof. To begin, note that the system average xt satisfies the recursion
xt+1 = xt + η
M M∑ i=1 vit (74)
Define the variables vt and ∆t, for t ≥ 1, as
vt = ∇f(xt) + 1
M M∑ i=1 ( vit −∇f(xit) ) ∆n = 1
M M∑ i=1 ( ∇f(xit)−∇f(xt) ) Then we can express the recursion equation 74 as
xt+1 = η (vt + ∆t)
We will show that this can be interpreted as a form of biased SGD and therefore we may apply Corollary 3.4. For the unbiased component vt, observe that
E [vt −∇fT (xt) | Ft−1] = E [ 1
M M∑ i=1 (vit −∇f(xit)) | Ft−1 ] = 0 (75)
and
E [ ‖vt −∇fT (xt)‖2 | Ft−1 ] ≤ E
[ 1
M M∑ i=1 |vit −∇f(xit)|2 ] = σ2v (76)
For the bias term, note that
‖∆t‖2 ≤ L2
M M∑ i=1 |xit − xt|2 = Vt
Assumption equation 3.1 follows from equation 75 and equation 76, while Assumption equation 3.2 follows from Proposition 5.3. The result then follows from Corollary 3.4.

heading: E ANALYSIS OF SVRG
text: For the analysis of SVRG, define the filtration {Ft}t=0,1,... as follows. F0 = σ(x1m) and for all s ≥ 1, Fs = σ ({ x1m } ∪ { ijt ∣∣∣ 0 ≤ t ≤ m− 1, 1 ≤ j ≤ s}) . We will leverage some prior results concerning the behavior of SVRG. The following is adapted from Reddi et al. (2016a).
Proposition E.1. Let Assumptions 2.1 and 2.2 hold. Let β > 0 and define the constants cm, cm−1, . . . , c0 as follows: cm = 0, and for 0 ≤ t ≤ m−1, let ct = ct+1(1+ηβ+2η2L2)+η2L3. Define Γt for 0 ≤ t ≤ m− 1 as Γt = η − ct+1ηβ − η2L− 2ct+1η2. Suppose that the step-size η and the analysis constant β are chosen so that Γt > 0 for 0 ≤ t ≤ m − 1, and set γ = inf0≤t<m Γt. Then for all s ≥ 1,
m−1∑ t=0 E[‖∇fT (xs+1t )‖2 | Fs−1] ≤ fT (x s m)− E[fT (xs+1m ) | Fs−1] γ . (77)
Furthermore, if η is of the form η = ξ/(Ln2/3) for some ξ ∈ (0, 1) and if the epoch length is set to m = bn/(3ξ)c, then there is a value for β such that γ ≥ ν(ξ)
Ln2/3 where ν(ξ) is a constant dependent
only on ξ. In particular, if ξ = 1/4 then
γ ≥ 1 40Ln2/3 . (78)
Proof. The proof of equation 77 follows from nearly the same reasoning used to establish Equation 10 in (Section B, Reddi et al. (2016a)), the only difference being that conditional expectations replace expectations in all of the relevant formulas. The details are left to the reader.
Formula equation 78 is a numerical inequality whose proof can be derived from the proof of Theorem 3 given in in (Appendix B, Reddi et al. (2016a)).

heading: PROOF OF PROPOSITION 6.1
text: Proof. First, note that τ( ) is a well-defined stopping time with respect to the filtration {Fs}s=0,1,.... For s = 1, 2, . . . define the random variables δs as
δs = m−1∑ t=0 ‖∇fT (xs+1t )‖2 − fT (x s m)− fT (xs+1m ) γ
It holds trivially that for all s ≥ 1, m−1∑ t=0 ‖∇fT (xs+1t )‖2 = fT (x s m)− fT (xs+1m ) γ + δs (79)
and by Proposition E.1, for all s ≥ 1,
E[δs | Fs−1] = m−1∑ t=0 E [∥∥∇fT (xs+1t )∣∣2 | Fs−1]− fT (xsm)− E[fT (xs+1m ) | Fs−1]γ
≤ 0. (80)
Summing Equation equation 79 over s = 1, . . . , q yields
q∑ s=1 m−1∑ i=0 ‖∇fT (xs+1i )‖2 = fT (x 1 m)− fT (xq+1m ) γ + q∑ s=1 δs, (81)
Rearranging terms and noting that fT (xq+1m ) ≥ f∗ results in
γ q∑ s=1 m−1∑ i=0 ‖∇fT (xs+1i )‖2 ≤ fT (x1m)− f∗ + γ q∑ s=1 δs. (82)
It follows that
γ q∑ s=1 ‖∇fT (xs+10 )‖2 ≤ fT (x1m)− f∗ + γ q∑ s=1 δs. (83)
For r ≥ 1, let τ( ) ∧ r be the stopping time which is the minimum of τ( ) and the constant value r. Applying Proposition A.1 together with Equation 80, it holds that
E τ( )∧r∑ s=1 δs  ≤ 0 (84) Furthermore, by definition of τ ,
E τ( )∧r∑ s=1 ‖∇fT (xs+10 )‖2  ≥ E (τ( )∧r)−1∑ s=1 ‖∇fT (xs+10 )‖2  ≥E (τ( )∧r)−1∑ s=1  = E[(τ( ) ∧ r)− 1]. (85)
Combining equation 83, equation 84, and equation 85 yields
γ E[(τ( ) ∧ n)− 1] ≤ fT (x1m)− f∗
Rearranging terms in the above yields
E[τ( ) ∧ n] ≤ fT (x 1 m)− f∗ γ + 1.
Applying the monotone convergence theorem, then,
E[τ( )] ≤ fT (x 1 m)− f∗ γ + 1.
Next, specialize η and m to η = ξ/(Ln2/3) and m = bn/(3ξ)c with ξ = 1/4. Then by Proposition E.1, γ ≥ 1/(40Ln2/3). Therefore,
E[τ( )] ≤ 40Ln 2/3(fT (x 1 m)− f∗) + 1.

heading: F GENERALIZATION ANALYSIS

heading: PROOF OF COROLLARY 7.2
text: Proof. Under our assumption on the Lipschitz property of y 7→ ∇xf(y, x), it holds that ‖∇fµ(xτ( ))‖ ≤ ‖∇fV (xτ( ))‖+Gd1(µ, µV ).
Squaring and taking expectations, while noting that d1 ≤ d2 (see Villani (2008), Remark 6.6), E[‖∇fµ(xτ( ))‖2] ≤ 2E[‖∇fV (xτ( ))‖2] + 2G2E[d2(µ, µV )2].
Using the Wasserstein concentration bound from Theorem 7.1 and the definition of τ( ), we obtain
E[‖∇fµ(xτ( ))‖2] ≤ 2 + 2G2κdJn−3/dV .

heading: G EXPERIMENTAL METHODOLOGY

heading: G.1 SSGD EXPERIMENTS
text: The neural network model used for these experiments is LSTNet Lai et al. (2018)with CUDA-aware MPI and extended to use the Stacked SGD training method. The objective function is the squared error between the prediction and the true sensor measurement, averaged over the dataset of training instances.
Our experiments compared SSGD (Algorithm 2) and SGD (1). For SSGD, each cluster has 4 worker nodes and 1 communication node. Hence, K = 4 and M is varied from 1 to 64. Since one physical node has four GPU devices, each physical node in the HPC environment is modeled a single local cluster in SSGD. The scalability of SSGD is compared with the basic parallel implementation of SGD, where all node synchronize with all-reduce collective communication call after each parameter update.
The experiment is conducted in a high performance computing environment that is equipped with 108 nodes. Each node has a dual Intel Xeon E5-2695v4 (Broadwell) CPU, four NVidia K40 GPUs, SAS-based local storage, and 256 GB of memory. The nodes are inter-connected with a non-blocking Infiniband EDR fabric.

heading: G.2 SVRG EXPERIMENTS
text: For the first experiment in this section, we trained a neural network with no hidden layer (i.e., a logistic classifier) for the MNIST classification task. For the second experiment, we trained a neural network with one hidden layer of n = 100 nodes on the CIFAR-10 classification task. The activation function in the hidden layer was the logistic function σ(x) = (1 + e−x)−1. In both cases, the objective function is the average cross entropy loss.
Our experiments compared SVRG (Algorithm 4) and SGD (Algorithm 1). The only modification to the algorithms was that mini-batches were used to make the gradient estimates. For each algorithm we determined the values of the learning rate and mini-batch size using a validation set. For the learning rate, we searched over the values η ∈ {0.001, 0.01, 0.1, 1.0}. For the mini-batch size, we considered values in {32, 64, 128}. The values used for full training were determined by running the training procedure for several epochs and evaluating the resulting model on a held-out portion of the training dataset. The parameters that gave best performance on this held-out dataset were used for full training. Using the found settings, we ran five independent runs of training, and report the mean and confidence bands (representing one standard deviation) for the expected IFO calls.

Review summaries:
In this paper, the authors consider stochastic optimization in the setting where a validation function is used to guide the termination of the algorithm. In more details, the algorithm terminates if the gradient of the validation function at an iterate is smaller than a threshold. In this framework, the authors consider several variants of SGD, including distributed variant and SVRG, for each of which the authors study the expected number of iterations for a prescribed accuracy under an assumption between the training and validation set.

This paper proposes an optimization approach in which the optimizer computes the gradient on a given function yet uses another to decide a stopping time. Conceptually those functions are empirical errors on train and validation folds in the most common setting, although the authors seem to use other settings later in the paper to consider decentralized optimization schemes. The authors introduce a bound on the Wasserstein distance between the train and validation distributions in their analysis which plays a crucial role in their results. The authors use these results to motivate variants of existing optimization algorithms.