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

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.

1. INTRODUCTION:

This work considers the minimization of a differentiable and possible nonconvex objective function:

min x∈R d 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 ∈ R d 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).

1.1. RELATED WORK:

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.

2. PRELIMINARIES:

Let f : R q × R d → 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 : R q × R d → R is bounded from below by f * ∈ R, and the derivative of f with respect to x is L-Lipschitz continuous:

∀y ∈ R q , x 1 , x 2 ∈ R d , ∇ x f (y, x 1 ) -∇ x f (y, x 2 ) ≤ L x 1 -x 2 .

At times, we will make a distinction between the training function f T , which is used to calculate gradients, and a validation function f V used to decide when to stop training. Assumption 2.2. The training function f T is defined using a set Y T ⊆ R q of n T elements as

f T (x) = 1 n T y∈Y T f (y, x), and the validation function f V is defined using a set Y V ⊆ R q of n V elements as f V (x) = 1 n V y∈Y V 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 f T and f V will be similar when the datasets Y T and Y V are similar. Formally, the datasets Y T and Y V determine probability measures µ T and µ Y , defined as

µ T = 1 n T y∈Y T δ y and µ V = 1 n V y∈Y V δ y , respectively

, where δ y is the delta measure δ y (A) = 1 y∈A for all sets A. We can compare these measures using the Wasserstein distance as in:

For q ≥ 1, p ≥ 1, we denote by P p (R q ) the probability measures on R q with finite moments of order p. Recall that a coupling of probability measures µ 1 and µ 2 is a probability measure γ on R q × R q such that for all measurable sets A, γ(A × R q ) = µ 1 (A) and γ(R q × 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 P p (R q ), denoted by d p , is defined as:

EQUATION

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 ∈ R d , ∇f V (x) -∇f T (x) ≤ Gd 1 (µ 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 → ∇ x f (y, x) is a G-Lipschitz function, uniformly for all x. Consider the following example: Example 2.4. Suppose that g : R q × R d → R is a smooth function. Let h : R d → R be the function that applies the hyperbolic tangent function to each of its components: h(x) = (tanh(x 1 ), . . . , tanh(x d )), and define f (y, x) = g(y, h(x)). Further suppose that the training data are bounded: y ≤ J for all y ∈ Y T . Then Assumption 2.3 is satisfied because

y → ∇ x f (x, y) is a Lipschitz function with G = sup y ≤J, x ≤ √ d ∂ 2 g ∂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 ∈ R d 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), ∇ x f (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 x 1 ∈ R d 2: t = 1 3: while ∇f V (x t ) 2 > do 4: for n = t to t + m -1 do 5: x n+1 = x n -ηh n 6: end 7: t = t + m 8: end 9: return x t

This section details our analysis of SGD with early stopping shown on the right in Algorithm 1. Starting from an initial point x 1 , at each iteration, the parameter is updated with an approximate gradient h n 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 h t is a sum of two components, v t and ∆ t , that represent an unbiased gradient estimate and an error term, respectively:

h t = v t + ∆ t .

(3)

Let {F t } t≥0 be a filtration such that x 1 is F 0 -measurable, and for all t > 1, the variables v t , ∆ t are F t -measurable. Our assumptions on the v t are as follows.

Assumption 3.1. For any t ≥ 1, it holds that

EQUATION

EQUATION

Assumption 3.1 states that the update directions v t are valid approximations to the gradient ∇f T (x t ), 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 V 1 , V 2 , . . ., and U 1 , U 2 , . . . such that for all t ≥ 1 the variables V t and U t are F t -measurable, ∆ t 2 ≤ V t , and the V t satisfy the following Lyapunov-type inequality: For constants α ∈ [0, 1) and β ≥ 0,

EQUATION

)

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 U t terms, prevents the error from vanishing completely. Note that the assumption would be satisfied in the unbiased case by simply setting V t = 0.

We can now state our result on the expected number of iterations for SGD with early stopping: Proposition 3.3. Let {x t } 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 τ

EQUATION

)

Furthermore, it holds with probability 1 that ∇f T (x τ ( )

) 2 ≤ 2 + 2G 2 d 1 (µ 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 > G 2 d 1 (µ V , µ T ) 2 . Let c ∈ (0, 1) and let the step-size be

EQUATION

)

Then

E[τ ( )] = O m 2 (1 + R/(1 -α)) (1 -c) c ( -G 2 d 1 (µ 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 y t ∈ Y T uniformly at random and setting v t = ∇ x f (y t , x t ) and that there

is a σ 2 v ≥ 0 such that ∀ y ∈ Y T , x ∈ R d , ∇ x f (y , x) -1 n T y∈Y T ∇ x f (y, x) 2 ≤ σ 2 v .

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 mn V + m 2 ( -G 2 d 1 (µ V , µ T ) 2 ) 2 + n V

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) ). 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).

4. STACKED SGD:

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 x i 1 (for workers), or

x i 1 (for communicators) 2: t = 1 3: while ∇f V 1 M M j=1 x j t 2 > do 4: for n = t to t + m -1 do 5:

if node i is a computation node then 6: 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 {F t } t≥0 as follows:

x i n+ 1 2 = x i n -ηv i n 7: x i n+1 = 1 K+1 x c(i) n+ 1 2 + K j=1 x j n+ 1 2 8: else (node i is a communication node) 9: x i n+ 1 2 = 1 M M j=1 x j n 10: x i n+1 = 1 K+1 x i n+ 1 2 + j∈c -1 (i) x j

F 0 = σ x i 1 1 ≤ i ≤ KM ∪ x i 1 1 ≤ i ≤ M , ∀t ≥ 1, F t = σ x i 1 , v i n 1 ≤ n ≤ t, 1 ≤ i ≤ KM ∪ x i 1 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,

EQUATION

EQUATION

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 i t be as defined in Line 6 of Algorithm 2. Suppose the step-size satisfies η < 1/(2LK). Define the variables V 1 , U 1 , V 2 , U 2 , . . . and the constants α, β as follows:

EQUATION

)

EQUATION

)

Then, for all t ≥ 1, it holds that

V t+1 ≤ αV t + U t and E[U t | F t-1 ] ≤ β.

In the preceding definitions, V t 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:

EQUATION

)

Let x t be the average of the communicator states at time t:

x t = 1 M M i=1 x i t . For > 0, define τ ( ) to be the stopping time τ ( ) = inf{n ≥ 1 | n ≡ 1 (mod m) and ∇f V ( x n ) 2 ≤ }. Then, E[τ ( )] = O m 2 (1 -c)c( -G 2 d 1 (µ 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 y j t ∈ Y T uniformly at random and setting v j t = ∇ x f (y j t , x j t ) and that there is a

σ 2 ≥ 0, such that ∀ y ∈ Y T , x ∈ R d , ∇ x f (y , x) -1 n T y∈Y T ∇ x f (y, x) 2 ≤ σ 2 v .

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(n V + mK) ( -G 2 d 1 (µ V , µ T ) 2 ) 2 + n V .

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.

5. DECENTRALIZED SGD:

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; a i,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 x t :

EQUATION

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 x 1 ∈ R d , initial parameters x i 1 . 2: t = 1 3: while ∇f V (x t ) 2 > do 4: for n = t to t + m -1 do 5: x i n+1 = M j=1 a i,j x j n -ηv i n 6: end 7: t = t + m 8: end 9: return x t

The intuitive justification for DSGD is that it may be more efficient compared to naive approaches to parallelizing SGD, since whenever a i,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 {F t } t≥0 as follows: 2 is assumed to satisfy ρ < 1.

F 0 = σ x i 1 1 ≤ i ≤ M , ∀t ≥ 1, F t = σ x i 1 , v i n 1 ≤ n ≤ t, 1 ≤ i ≤ M

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

EQUATION

EQUATION

)

For the analysis, we show that the sequence of averages x t 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 x i t 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 V 1 , U 1 , V 2 , U 2 , .

. . and the constants α, β as follows:

EQUATION

EQUATION

) α = (1 + ρ) 2ρ √ ρ + 1 2 2 (19c) β = η L √ 2 1 -ρ σ 2 v (19d)

Then for all t ≥ 1 it holds that

V t+1 ≤ αV t + U t and E[U t | F t-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:

x i 1 = x j 1 for all 1 ≤ i, j ≤ M . Let c ≤ 1- √ ρ 4 √ 2 and define R = L √ 2σ 2 v /(1 -ρ)

Algorithm 4 SVRG with early stopping

1: input: Initial point x 1 m ∈ R d 2: for s = 1, 2, . . . do 3: x s+1 0 = x s m 4: g s+1 = 1 n T y∈Y T ∇f T (y, x s+1 0 ) 5: if g s+1 2 ≤ then return x s+1 0 6: for t = 0 to m -1 do 7:

Sample y s t uniformly at random from Y T 8:

v s t = ∇f (y s t , x s+1 t

) -∇f (y s t , x s+1 0 ) + g s+1 9:

x 

EQUATION

)

For > 0 define τ ( ) to be the first time the norm of the gradient of the validation function falls below ; That is, τ

EQUATION

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.

6. SVRG:

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 x s+1 t defined by Algorithm 4. Let ξ = 1/4 and suppose that the step-size is set to η = ξ/(Ln 2/3 T ) and the epoch length is m = n T /(3ξ) . For > 0, define τ ( ) to be the stopping time τ

( ) = inf s ≥ 1 ∇f T (x s+1 0 ) 2 ≤ . Then, E[τ ( )] ≤ 1 + (40Ln 2/3 T (f T (x 1 m ) -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 n T + 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((n 5/3 T / ) + n T ).

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 .

7. GENERALIZATION PROPERTIES:

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 µ : R d → 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 Y V and Y T . 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 R d , such that J = E µ y 3 1/3 < ∞. Then, there is a constant

κ d , such that E[d 2 (µ, µ V ) 2 ] ≤ κ d Jn -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 Y V is an empirical version of µ, and

y → ∇ x f (y, x) is uniformly G-Lipschitz. If x τ ( ) is the output of Algorithm 1, then E[ ∇f µ (x τ ( ) ) 2 ] ≤ 2 + 2G 2 κ d Jn -3/d V .

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 3 rd 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 n T . 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.

8. DISCUSSION:

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.

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

A PRELIMINARIES

Our analyses make use of a quadratic bound for the training function which follows from Assumption 2.1:

EQUATION

)

A.1 STOCHASTIC PROCESSES

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 x 1 and defines a sequence of random variables {x t (ω)} 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 {F t } t=0,1,... is an increasing sequence of σ-algebras, with the interpretation that F t represents the information available to an algorithm up to and including time t. A random variable x : Ω → R d is said to be F t 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 {F t } t=0,1,... . Suppose there is a number c < ∞ such that τ ≤ c with probability one. Let x 1 , x 2 , . . . be any sequence of random variables such that each x t is F t -measurable and

EQUATION

)

Proof. This is a consequence of the optional stopping theorem (Theorem 10.10 in Williams (1991) ).

Define S 0 = 0 and for t ≥ 1, let Proof. For convenience, define the random variables δ t for t = 1, 2, . . . as δ t = v t -∇f T (x t ). From equation 22, it holds that

S t = t i=1 (x i -E[x i | F i-1 ]). Then S 0 , S 1 , . . .

f T (x t+1 ) ≤ f T (x t ) -η∇f T (x t ) T (∇f T (x t ) + δ t + ∆ t ) + L 2 η 2 ∇f T (x t ) + δ t + ∆ t 2 .

Under review as a conference paper at ICLR 2020 Summing this over t = 1, . . . , m yields

f T (x m+1 ) ≤ f T (x 1 ) - m t=1 η∇f T (x t ) T (∇f T (x t ) + δ t + ∆ t ) + m t=1 L 2 η 2 ∇f T (x t ) + δ t + ∆ t 2 = f T (x 1 ) - m t=1 η 1 - L 2 η ∇f T (x t ) 2 - m t=1 η(1 -Lη)∇f T (x t ) T δ t + m t=1 L 2 η 2 δ t 2 - m t=1 η(1 -Lη)∇f T (x t ) T ∆ t + m t=1 L 2 η 2 ∆ t 2 + m t=1

Lη 2 δ T t ∆ t .

(24)

Note that, in general, for any numbers a, b it is the case that |ab| ≤ 1 2 a 2 + 1 2 b 2 . Then

EQUATION

)

Combining Equations equation 24, equation 25, and equation 26, we obtain

f T (x m+1 ) ≤ f (x 1 ) - m t=1 η 1 - L 2 η - 1 2 (1 -Lη) ∇f T (x t ) 2 - m t=1 η(1 -Lη)∇f T (x t ) 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 (x 1 ) - m t=1 η 1 2 ∇f T (x t ) 2 - m t=1 η(1 -Lη)∇f T (x t ) T δ t + m t=1 Lη 2 δ t 2 + m t=1 1 2 η(1 + Lη) ∆ t 2 .

Rearranging terms and noting that f T (x m+1 ) ≥ f * , this yields

EQUATION

)

where in the second inequality we used the assumptions that η ≤ 1/L and ∆ t 2 ≤ V t . Based on our assumption that ∇f

EQUATION

EQUATION

Combining equation 28 and equation 29 results in

EQUATION

)

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

EQUATION

and using Proposition A.1 with Assumption 5 gives

EQUATION

Next, according to conditions equation 6, and equation 7, it holds for any m ≥ 1 and with probability one that

EQUATION

and by equation 8 together with Proposition A.1,

EQUATION

Combining equation 33 and equation 34, then

EQUATION

Furthermore, by definition of τ ,

EQUATION

Here we used that 

n-1 t=1 1 t≡1 (mod m) = n -1 m ≥ n -

2m (E[τ ( ) ∧ n] -1) ≤ η 2 G 2 d 1 (µ V , µ T ) 2 E[τ ( ) ∧ n] m + 1 + f T (x 1 ) -f * + Lη 2 σ 2 v E [τ ( ) ∧ n] + η β 1 -α (E[τ ( ) ∧ n] + 1) .

This can be rearranged into

η m -2Lη 2 σ 2 v - 2ηβ 1 -α - η m G 2 d 1 (µ V , µ T ) 2 E[τ ( ) ∧ n] ≤ ηG 2 d 1 (µ V , µ T ) 2 + 2(f T (x 1 ) -f * ) + 2η β 1 -α + η m .

which in turn is equivalent to

EQUATION

)

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 f V and f T reversed, and using the definition of τ ( ), it follows that

∇f T (x τ ( ) ) ≤ 2 ∇f V (x τ ( ) ) 2 + 2G 2 d 1 (µ V , µ T ) 2 ≤ 2 + 2G 2 d 1 (µ V , µ T ) 2 .

PROOF OF COROLLARY 3.4

Proof. According to the assumption on the step-size η (Equation equation 20),

EQUATION

Combining these inequalities with the conclusion of Proposition 3.3 (relation equation 9) yields

EQUATION

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

EQUATION

Next,

EQUATION

Combining equation 40 with equation 41, equation 42 and equation 43, we find that

EQUATION

PROOF OF COROLLARY 3.5

Proof. If the algorithm runs until iteration τ ( ), then the number of times that the full gradient of f V is calculated is τ ( )/m ≤ τ ( )/m + 1, and the number of IFO calls for the training function is τ ( ) -1. Therefore

EQUATION

Note that under our assumption on the gradient estimates v t , we are in the unbiased setting where R = 0. Combining equation 44 with R = 0 and equation 45, we obtain

EQUATION

Using c = 1/2 and neglecting lower-order terms, then,

E[IFO( )] = O mn V + m 2 ( -G 2 d 1 (µ V , µ T ) 2 ) 2 + n V .

C ANALYSIS OF STACKED SGD PROOF OF PROPOSITION 4.2

Proof. Note that the variables inside the stacked algorithm satisfy the following identities: For all

EQUATION

and for all pairs communication nodes nodes 1 ≤ j ≤ M , and 1 ≤ k ≤ M , and t ≥ 1,

EQUATION

)

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,

EQUATION

Step A follows from the definition of x i t+1 , x j t+1 in Algorithm 2 and Equation equation 48.

Step B follows from the definition of x k t+ 1 2 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,

EQUATION

where the last inequality uses the Lipschitz gradient property (Assumption 2.1) and the triangle inequality. Combining equation 49 and equation 50 yields

EQUATION

Note that for any k 1 > 0 and all a, b we have

EQUATION

)

Combining equation 51 and equation 52 while using the assumption that η ≤ 1/(2LK) results in

EQUATION

)

Let k 1 = (K+3/4) 2 (K+1/2) 2 -1. Then by the definition of α,

EQUATION

Furthermore,

EQUATION

Multiplying each side of equation 53 by L 2 /M 2 , 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,

V t+1 ≤ αV t + U t .

It remains to confirm that E[U t | F t-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 [U t | F t-1 ] ≤ η • KL (K + 5/8) × 1 M 2 M i=1 M j=1 E   1 K k∈c -1 (i) v k t -∇f T ( x i t )   -   1 K k∈c -1 (j) v k t -∇f T ( x j t )   2 F t-1 = η KL (K + 5/8) 1 M 2 M i=1 M j=1 2 σ 2 v K + σ 2 v K = β.

PROOF OF PROPOSITION 4.3

Proof. Note that according to Line 9 of SSGD,

x t = x 1 t+ 1 2 = . . . = x M t+ 1 2 (56)

and Lines 7 and 10 mean that for all j with c(j) = i,

EQUATION

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 i t+ 1 2 + j∈c -1 (i) x j t+ 1 2   B =   1 K + 1 x t + 1 M M i=1 1 K + 1 j∈c -1 (i) x j t -ηv j t   C =   1 K + 1 x t + 1 M M i=1 1 K + 1 j∈c -1 (i) ( x i t -ηv j t )   D = 1 K + 1 x t + K K + 1 x t - 1 M M i=1 1 K + 1 j∈c -1 (i) ηv j t E = x t -η   1 M M i=1 1 K + 1 j∈c -1 (i) v j t   .

For Step A, note that x t is the average of the x i t , and then use definition of the x i t from Line 9 of the SSGD algorithm. Step B follows from equation 56 and the definition of the x i t+ 1 2 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 i t .

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) v j t   = x t -η K K + 1   1 M M i=1 1 K j∈c -1 (i) v j t   = x t - ηK K + 1 (v t + ∆ t ) ,

where v t and ∆ t are

EQUATION

and

EQUATION

Based on the definition of v t in Equation equation 58 and on Assumption 4.1, for all t ≥ 1 it holds that

E [v t -∇f T ( x t ) | F t-1 ] = 0, E v t -∇f T ( x t ) 2 | F t-1 ≤ σ 2 v K .

Thus Assumption 3.1 is confirmed. Next we consider the Lyapunov condition of Assumption 3.2. Let the variables U t , V t for t ≥ 1 and the constants α, β be defined as in equation 14a-equation 14d.

By our assumption that each node has the same initial state, then 

EQUATION

)

According to Corollary 3.4, then, a step-size of

EQUATION

Dropping the lower order terms, we see that

EQUATION

)

PROOF OF COROLLARY 4.4

Proof. If SSGD runs until iteration τ ( ), then number of times that the full gradient of f V is calculated is τ ( )/m ≤ τ ( )/m + 1, and the number of IFO calls for the training function is (τ ( ) -1)K. Therefore

EQUATION

Combining equation 67 and equation 70, then,

y n+1 ≤ a -a ∞ + ηL √ 8 y n + η 8 M i=1 |v i n -∇f (x i n )| 2 ≤ √ ρ + ηL √ 8 y n + η 8 M i=1 |v i n -∇f (x i n )| 2

In the second step we have applied Assumption equation 5.1 and Lemma D.1. Squaring this equation, for any k 1 > 0 it holds that

EQUATION

) Let k 1 = 1-ρ 2ρ (in which case 1 + 1 k1 = 1+ρ 1-ρ ).

Multiplying each side of 71 by L 2 /M and noting that

EQUATION

It follows from the variance bound in Assumption 5.2 that

E [U t | F t-1 ] ≤ 8 η 2 L 2 1 + ρ 1 -ρ σ 2 v (73) Using the assumption that η ≤ 1- √ ρ 4L √

2 , then equation 72 and equation 73 become, respectively,

V t+1 ≤ (1 + ρ) 2ρ √ ρ + 1 2 2 V t + U t and E [U t | F t-1 ] ≤ √ 2 η(1 - √ ρ) L 1 + ρ 1 -ρ σ 2 v ≤ η L √ 2 1 -ρ σ 2 v PROOF OF PROPOSITION 5.4

Proof. To begin, note that the system average x t satisfies the recursion

EQUATION

)

Define the variables v t and ∆ t , for t ≥ 1, as

v t = ∇f (x t ) + 1 M M i=1 v i t -∇f (x i t ) ∆ n = 1 M M i=1 ∇f (x i t ) -∇f (x t )

Then we can express the recursion equation 74 as

x t+1 = η (v t + ∆ t )

It holds trivially that for all s ≥ 1,

m-1 t=0 ∇f T (x s+1 t ) 2 = f T (x s m ) -f T (x s+1 m ) γ + δ s (79)

and by Proposition E.1, for all s ≥ 1, 

E[δ s | F s-1 ] = m-1 t=0 E ∇f T (x s+1 t ) 2 | F s-1 - f T (x s m ) -E[f T (x s+1 m ) | F s-1 ] γ ≤ 0.

EQUATION

Rearranging terms and noting that f T (x q+1 m ) ≥ f * results in

EQUATION

It follows that

EQUATION

For r ≥ 1, let τ ( ) ∧ r be the stopping time which is the minimum of τ ( ) and the constant value r. 

E[τ ( ) ∧ n] ≤ f T (x 1 m ) -f * γ + 1.

Applying the monotone convergence theorem, then,

E[τ ( )] ≤ f T (x 1 m ) -f * γ + 1.

Next, specialize η and m to η = ξ/(Ln 2/3 ) and m = n/(3ξ) with ξ = 1/4. Then by Proposition E.1, γ ≥ 1/(40Ln 2/3 ). Therefore, E[τ ( )] ≤ 40Ln 2/3 (f T (x 1 m )f * ) + 1.

F GENERALIZATION ANALYSIS:

PROOF OF COROLLARY 7.2

Proof. Under our assumption on the Lipschitz property of y → ∇ x f (y, x), it holds that ∇f µ (x τ ( ) ) ≤ ∇f V (x τ ( ) ) + Gd 1 (µ, µ V ).

Squaring and taking expectations, while noting that d 1 ≤ d 2 (see Villani (2008) , Remark 6.6),

E[ ∇f µ (x τ ( ) ) 2 ] ≤ 2E[ ∇f V (x τ ( ) ) 2 ] + 2G 2 E[d 2 (µ, µ V ) 2 ].

Using the Wasserstein concentration bound from Theorem 7.1 and the definition of τ ( ), we obtain 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.

G.2 SVRG EXPERIMENTS:

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.

:



Figures and Tables:

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)|})

approximate stationarity for MNIST (top) and CIFAR-10 (bot.). and let α be as in Equation equation 19c. Let the step-size be

is a martingale with respect to the filtration {F t } t=0,1,... , and the optional stopping theorem implies E[S τ ] = E[S 0 ]. But E[S 0 ] = 0, and therefore E[S τ ] = 0, which is equivalent to Equation equation 23. B ANALYSIS OF BIASED SGD PROOF OF PROPOSITION 3.3

V 1 = 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

, equation 84, and equation 85 yieldsγ E[(τ ( ) ∧ n) -1] ≤ f T (x 1 m )f * Rearranging terms in the above yields

[ ∇f µ (x τ ( ) ) 2 ] ≤ 2 + 2G 2 κ d JnThe neural network model used for these experiments is LSTNetLai 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.

Combining equation 27, equation 30, equation 31, equation 32, equation 35 and equation 36 results in η

Table content: <table/>

Peer review (iclr2020):
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. While the use of a validation function is useful for early stopping, it introduces additional cost. While bounds on the expected number of iterations are derived for several variants of SGD, it seems that most arguments are adapted from the existing analysis to take into account the validation function. In Corollary 3.4 and Corollary 3.5, the bound is an increasing function of m. This suggests that the best choice would be m=1. However, in this case, one needs to calculate the gradient of the validation function at each iteration, which may wastes a lot of computation. The authors consider constant step sizes. In practice, step sizes are often decreasing along the optimization. Can the analysis be extended to cover the case with decreasing step sizes? In eq (30), there is a missing factor of 2. There is a required in the results. Therefore, to achieve a high accuracy we need to be small. How many numbers of sample size to make small? This has an influence on the computational cost. ---------------------- After rebuttal: The authors do not respond. I would like to keep my original score.