Title: CORESETS FOR ACCELERATING INCREMENTAL GRADIENT METHODS

Abstract: Many machine learning problems reduce to the problem of minimizing an expected risk. Incremental gradient (IG) methods, such as stochastic gradient descent and its variants, have been successfully used to train the largest of machine learning models. IG methods, however, are in general slow to converge and sensitive to stepsize choices. Therefore, much work has focused on speeding them up by reducing the variance of the estimated gradient or choosing better stepsizes. An alternative strategy would be to select a carefully chosen subset of training data, train only on that subset, and hence speed up optimization. However, it remains an open question how to achieve this, both theoretically as well as practically, while not compromising on the quality of the final model. Here we develop CRAIG, a method for selecting a weighted subset (or coreset) of training data in order to speed up IG methods. We prove that by greedily selecting a subset S of training data that minimizes the upper-bound on the estimation error of the full gradient, running IG on this subset will converge to the (near)optimal solution in the same number of epochs as running IG on the full data. But because at each epoch the gradients are computed only on the subset S, we obtain a speedup that is inversely proportional to the size of S. Our subset selection algorithm is fully general and can be applied to most IG methods. We further demonstrate practical effectiveness of our algorithm, CRAIG, through an extensive set of experiments on several applications, including logistic regression and deep neural networks. Experiments show that CRAIG, while achieving practically the same loss, speeds up IG methods by up to 10x for convex and 3x for non-convex (deep learning) problems.

1. INTRODUCTION:

Mathematical optimization lies at the core of training large-scale machine learning systems, and is now widely used over massive data sets with great practical success, assuming sufficient data resources are available. Achieving this success, however, also requires large amounts of (often GPU) computing, as well as concomitant financial expenditures and energy usage (Strubell et al., 2019) . Significantly decreasing these costs without decreasing the learnt system's resulting accuracy is one of the grand challenges of machine learning and artificial intelligence today (Asi & Duchi, 2019) .

Training machine learning models often reduces to the problem of optimizing a regularized empirical risk function. Given a convex loss l, and a µ-strongly convex regularizer r, one aims to find model parameter vector x * over the parameter space X that minimizes the loss f over the training data V :

EQUATION

where V = {1, . . . , n} is an index set of the training data, and functions f i : R d → R are associated with training examples (a i , y i ), where a i ∈ R d is the feature vector, and y i is the point i's label.

The standard Gradient Descent can find the minimizer of this problem, but requires repeated computations of the full gradient ∇f (x)-sum of the gradients over all training data points/functions i-and is therefore prohibitive for massive data sets. This issue is further exacerbated in case of deep neural networks where gradient computations (backpropagation) are expensive. Incremental Gradient (IG) methods, such as Stochastic Gradient Descent and its accelerated variants, including SGD with momentum (Qian, 1999), Adagrad (Duchi et al., 2011) , Adam (Kingma & Ba, 2014) , SAGA (Defazio et al., 2014) , and SVRG (Johnson & Zhang, 2013) iteratively estimate the gradient on random subsets/batches of training data. While this provides an unbiased estimate of the full gradient, the randomized batches introduce variance in the gradient estimate (Hofmann et al., 2015) , and therefore stochastic gradient methods are in general slow to converge (Johnson & Zhang, 2013; Defazio et al., 2014) . The majority of the work speeding up IG methods has thus primarily focused on reducing the variance of the gradient estimate (SAGA (Defazio et al., 2014) , SVRG (Johnson & Zhang, 2013) , Katysha (Allen-Zhu, 2017)) or more carefully selecting the gradient stepsize (Adagrad (Duchi et al., 2011 ), Adadelta (Zeiler, 2012) , Adam (Kingma & Ba, 2014) ).

However, the direction that remains largely unexplored is how to carefully select a small subset S ⊆ V of the full training data V , so that the model can only be trained on the subset S while still (approximately) converging to the globally optimal solution (i.e., the model parameters that would be obtained if training/optimizing on the full V ). If such a subset S can be quickly found, then this would directly lead to a speedup of |V |/|S| (which can be very large if |S| |V |) per epoch of IG.

There are four main challenges in finding such a subset S. First is that a guiding principle to select S is unclear. For example, selecting training points close to the decision boundary might allow the model to fine tune the decision boundary, while picking the most diverse set of data points would allow the model to get a better sense of the training data distribution. Second is that finding S must be fast as otherwise identifying the set S may take longer than the actual optimization, and so no overall speed-up would be achieved. Third is that finding a subset S is not enough. One also has to decide on an ordering over S and a gradient stepsize for each data point in S, as they affect the convergence. And last, while the method might work well empirically on some data sets, one also requires theoretical understanding and mathematical convergence guarantees.

Here we develop CoResets for Accelerating Incremental Gradient descent (CRAIG), for selecting a subset of training data points to speed up training of large machine learning models. Our key idea is to directly approximate the gradient. That is, we aim to find a weighted and ordered subset S of training data V that is representative of the full gradient of V . We prove that the subset S that minimizes an upper-bound on the error of estimating the full gradient of V maximizes a submodular facility location function. As a result, the subset S can be efficiently found using a fast greedy algorithm. A further benefit of our approach is that set S is created incrementally which induces a natural ordering over data in S. Thus, rather than processing data points in a random or arbitrary order, CRAIG processes them using in the order provided by the procedure, which we show further speeds up the convergence of the method.

We also provide theoretical analysis of CRAIG and prove the convergence of our approach. In particular, for a µ-strongly convex risk function and a subset S selected by CRAIG that estimates the full gradient by an error of at most , we prove that IG method with diminishing stepsize α k = α/k τ at epoch k (with 0 < τ < 1 and 0 < α), converges to an 2R /µ neighborhood of the optimal solution at rate O(1/ √ k). Here, R = min{d 0 , (rγ max C + )/µ} where d 0 is the initial distance to the optimum, C is an upper-bound on the norm of the gradients, r = |S|, and γ max is the largest weight for the elements in the subset obtained by CRAIG. Moreover, we prove that if in addition to the strong convexity, component functions have smooth gradients, IG with the same diminishing stepsize on subset S converges to a 2 /µ neighborhood of the optimum solution at rate O(1/k τ ).

The above convergence rates are the same as convergence rate of IG on V for a strongly convex risk (and smooth component functions), and therefore IG on S converges in the same number epochs as IG on the full V . But because every epoch only uses a subset S of the data, it requires fewer gradient computations and thus leads to a |V |/|S| speedup over traditional IG methods, while still (approximately) converging to the optimal solution. Furthermore, CRAIG only requires the knowledge of estimated gradient differences and does not involve any (exact) gradient calculations. Therefore, CRAIG can be used as a simple preprocessing step before IG starts and no additional storage or gradient calculations are required during IG, which makes CRAIG extremely practical. As such CRAIG can be used to speed up any existing IG methods, including IG, Adam, SAGA, SVRG as we show in the experiments section. and accuracy as the underlying incremental gradient descent methods, speeds up gradient methods by up to 10x for convex and 3x for non-convex loss functions. We also demonstrate that the deliberate ordering scheme of the CRAIG algorithm significantly improves convergence time.

2. RELATED WORK:

Convergence of IG methods has been long studied under various conditions (Zhi-Quan & Paul, 1994; Mangasariany & Solodovy, 1994; Bertsekas, 1996; Solodov, 1998; Tseng, 1998) , however IG's convergence rate has been characterized only more recently (see (Bertsekas, 2015) for a survey). In particular, Nedić & Bertsekas (2001) provides a O(1/ √ k) convergence rate for diminishing stepsizes α k per epoch k under a strong convexity assumption, and Gürbüzbalaban et al. (2015) proves a O(1/k τ ) convergence rate with diminishing stepsizes α k = Θ(1/k τ ) for τ ∈ (0, 1] under an additional smoothness assumption for the components. While these works provide convergence on the full dataset, our analysis provides the same convergence rates on subsets obtained by CRAIG.

It has been empirically observed that ordering of data significantly affects the convergence rate of IG. However, finding a favorable ordering for IG has been a long standing open question. Among the few results are that of (Recht & Re, 2012) showing that without-replacement random sampling improves convergence of IG for least means squares problem, and the very recent result of (Gurbuzbalaban et al., 2017) showing that a Random Reshuffling (RR) method with iterate averaging and a diminishing stepsize Θ(1/k τ ) for τ ∈ (1/2,1) converges at rate Θ(1/k 2τ ) with probability one in the suboptimality of the objective value, thus improving upon the Ω(1/k) rate of SGD. Contrary to the above randomized analysis, we propose the first deterministic ordering on the data points and empirically show that the ordering provided by CRAIG provides a significant speedup for the convergence of IG.

Techniques for speeding up SGD, are mostly focused on variance reduction techniques (Roux et al., 2012; Shalev-Shwartz & Zhang, 2013; Johnson & Zhang, 2013; Hofmann et al., 2015; Allen-Zhu et al., 2016) , and accelerated gradient methods when the regularization parameter is small (Frostig et al., 2015; Lin et al., 2015; Xiao & Zhang, 2014) . Very recently, Hofmann et al. (2015) ; Allen-Zhu et al. (2016) exploited neighborhood structure to further reduce the variance of stochastic gradient descent and improve its running time. Our CRAIG method and analysis are complementary to variance reduction and accelerated methods. CRAIG can be applied to these methods as well to speed them up (as we show in experiments).

3. CORESETS FOR INCREMENTAL GRADIENT DESCENT (CRAIG):

We proceed as follows: First, we define an objective function L for selecting an optimal set S of size r that best approximates the gradient of the full training dataset V of size n. Then, we show that L can be turned into a submodular function F and thus S can be efficiently found using a greedy algorithm. Crucially, we also show that the approximation error between the estimated and the true gradient can be efficiently minimized in a way that is independent of the actual optimization procedure and thus CRAIG can simply be used as a preprocessing step before the actual optimization starts.

Incremental gradient methods aim at estimating the full gradient i∈V ∇f i (x) over V by iteratively making a step based on the gradient of every function f i . Our key idea in CRAIG is that if we can find a small subset S such that the weighted sum of the gradients of its elements closely approximates the full gradient over V , we can apply IG only to the set S (with stepsizes equal to the weight of the elements in S), and we should still converge to the (approximately) optimal solution, but much faster. Specifically, our goal in CRAIG is to find the smallest subset S ⊆ V and corresponding per-element stepsizes γ j > 0 that approximate the full gradient with an error at most > 0 for all the possible values of the optimization parameters x ∈ X .1 

EQUATION

Given such an S * and associated weights {γ} j , we are guaranteed that gradient updates on S * will be similar to the gradient updates on V regardless of the value of x.

Unfortunately, directly solving the above optimization problem is not feasible, due to two problems. Problem 1: Eq. ( 2) requires us to calculate the gradient of all the functions f i over the entire space X , which is too expensive and would not lead to overall speedup. In other words, it would appear that solving for S * is as difficult as solving problem (1), as it involves calculating i∈V ∇f i (x) for various x ∈ X . And Problem 2: even if calculating the normed difference between the gradients in Eq. ( 2) would be fast, as we discuss later finding the optimal subset S * in NP-hard. In the following, we address the above two challenges and discuss how we can quickly find a near-optimal subset S.

3.1. UPPER-BOUND ON THE ESTIMATION ERROR:

We first address Problem 1, i.e., how to quickly estimate the error/discrepancy of the weighted sum of gradients of functions f j associate with data points j ∈ S, vs the full gradient, for every x ∈ X .

Let S be a subset of r data points. Furthermore, assume that there is a mapping ς x : V → S that for every x ∈ X assigns every data point i ∈ V to one of the elements j in S, i.e., ς

x (i) = j ∈ S. Let C j = {i ∈ [n]|ς(i) = j} ⊆

V be the set of data points that are assigned to j ∈ S, and γ j = |C j | be the number of such data points. Hence, {C j } r j=1 form a partition of V . Then, for any arbitrary (single) x ∈ X we can write

EQUATION

Subtracting and then taking the norm of the both sides, we get an upper bound on the error of estimating the full gradient with the weighted sum of the gradients of the functionsf j forj ∈ S. I.e.,

EQUATION

where the inequality follows from the triangle inequality. The upper-bound in Eq. ( 5) is minimized when ς x assigns every i ∈ V to an element in S with most gradient similarity at x, or minimum Euclidean distance between the gradient vectors at x. That is:

EQUATION

The right hand side of Eq. ( 6) is minimized when S is the set of r medoids (exemplars) for all the components in the gradient space. So far, we considered upper-bounding the gradient estimation error at a particular x ∈ X . To bound the estimation error for all x ∈ X , we consider a worst-case approximation of the estimation error over the entire parameter space X . Formally, we define a distance metric d ij between gradients of f i and f j as the maximum normed difference between ∇f i (x) and ∇f j (x) over all x ∈ X :

EQUATION

Thus, by solving the following minimization problem, we obtain the smallest weighted subset S * that approximates the full gradient by an error of at most for all x ∈ X : 1

EQUATION

Note that Eq. ( 8) requires that the gradient error is bounded overX . However, we show (Appendix B) for several classes of convex problems, including linear regression, ridge regression, logistic regression, and regularized support vector machines (SVMs), the normed gradient difference between data points can be efficiently boundedly approximated by (Allen-Zhu et al., 2016; Hofmann et al., 2015) :

EQUATION

Note when x is bounded for all x ∈ X , i.e., max x∈X O( x ) < ∞, upper-bounds on the Euclidean distances between the gradients can be pre-computed. This is crucial, because it means that estimation error of the full gradient can be efficiently bounded independent of the actual optimization problem (i.e., point x). Thus, these upper-bounds can be computed only once as a pre-processing step before any training takes place, and then used to find the subset S by solving the optimization problem (8). We address upper-bounding the normed difference between gradients for deep models in Section 3.4.

Algorithm 1 CRAIG (CoResets for Accelerating Incremental Gradient descent)

Input: Set of component functions f i for i ∈ V = {1, • • • , n}.

Output: Subset S ⊆ V with corresponding per-element stepsizes {γ} j∈S , and an ordering σ.

1: S 0 ← ∅, s 0 = 0, i = 0. 2: while F (S) < L({s 0 }) -do 3: j ∈ argmax e∈V \Si-1 F (e|S i-1 )

4:

S i = S i-1 ∪ {j} 5: σ i = j 6: i = i + 1 7: end while 8: for j = 1 to |S| do 9: γ j = i∈V I j = argmin s∈S max x∈X ∇f i (x) -∇f s (x) 10: end for 3.2 THE CRAIG ALGORITHM

Optimization problem (8) produces a subset S of elements with their associated weights {γ} j∈S or per-element stepsizes that closely approximates the full gradient. Here, we show how to efficiently approximately solve the above optimization problem in order to find a near-optimal subset S.

The optimization problem ( 8) is NP-hard as it involves calculating the value of L(S) for all the 2 |V | subsets S ⊆ V . We show, however, to transform it into a submodular set cover problem, for which efficient approximation algorithms exist.

Formally, F is submodular if F (S ∪ {e}) -f (S) ≥ F (T ∪ {e}) -F (T )

, for any S ⊆ T ⊆ V and e ∈ V \ T . We denote the marginal utility of an element s w.r.t. a subset S as F (e|S) = F (S ∪ {e}) -F (S). Function F is called monotone if F (e|S) ≥ 0 for any e ∈ V \S and S ⊆ V . The submodular cover problem is defined as finding the smallest set S that achieves utility ρ. Precisely,

EQUATION

)

Although finding S * is NP-hard since it captures such well-known NP-hard problems as Minimum Vertex Cover, for many classes of submodular functions (Nemhauser et al., 1978; Wolsey, 1982) , a simple greedy algorithm is known to be very effective. The greedy algorithm starts with the empty set S 0 = ∅, and at each iteration i, it chooses an element e ∈ V that maximizes (e|S i-1 ), i.e., S i = S i-1 ∪ {argmax e∈V F (e|S i-1 ) -F (S i-1 )}. Greedy gives us a logarithmic approximation, i.e. |S| ≤ 1 + ln(max e F (e|∅)) |S * |. The computational complexity of the greedy algorithm is O(|V | • |S|). However, its running time can be reduced to O(|V |) using stochastic algorithms (Mirzasoleiman et al., 2015a) and further improved using lazy evaluation (Minoux, 1978) , and distributed implementations (Mirzasoleiman et al., 2015b; 2016) .

Given a subset S ⊆ V , the facility location function quantifies the coverage of the whole data set V by the subset S by summing the similarities between every i ∈ V and its closest element j ∈ S. Formally, facility location is defined as F f l (S) = i∈V max j∈S s i,j , where s i,j is the similarity between i, j ∈ V . The facility location function has been used in a number of applications, including scene and documents summarization (Simon et al., 2007; Lin & Bilmes, 2012) .

By introducing an auxiliary element s 0 we can turn L(S) in Eq. ( 8) into a monotone submodular facility location function,

EQUATION

where L({s 0 }) is a constant. In words, F measures the decrease in the estimation error associated with the set S versus the estimation error associated with just the auxiliary element. It is easy to see that for suitable choice of s 0 , maximizing F is equivalent to minimizing L. Therefore, we apply the greedy algorithm to approximately solve the following problem to get the subset S defined in 8:

EQUATION

At every step, the greedy algorithm selects an element that reduces the upper bound on the estimation error the most. In fact, the size of the smallest subset S that estimates the full gradient by an error of at most depends on the structural properties of the data. Intuitively, as long as the marginal gains of facility location are considerably large, we need more elements to improve our estimation of the full gradient. Having found S, the weight γ j of every element j ∈ S is the number of components that are closest to it in the gradient space, and are used as stepsize of element j ∈ S during IG. The pseudocode for CRAIG subset selection method is outlined in Algorithm 1.

3.3. ORDERING ON THE ELEMENTS OF THE SUBSET:

Notice that CRAIG creates subset S incrementally one element at a time, which produces a natural order σ to the elements in S. Adding the element with largest marginal gain j ∈ argmax e∈V F (e|S i-1 ) improves our estimation from the full gradient by an amount bounded by the marginal gain. Formally, at every step i, we have F (S i ) ≥ (1 -e -i/|S| )F (S * ), and hence

EQUATION

Intuitively, the first elements of the ordering contribute the most to provide a close approximation of the full gradient and the rest of the elements further refine the approximation. Hence, the first incremental gradient updates gets us close to x * , and the rest of the updates further refine the solution.

We show experimentally that processing data points in the order of S leads to faster convergence than when we consider elements in S in a random order. We defer the formal proof to future work.

3.4. APPLICATION OF CRAIG TO DEEP NETWORKS:

Incremental gradient methods, including SGD with momentum (Qian, 1999), Adam (Kingma & Ba, 2014) and Adagrad (Duchi et al., 2011) are widely used to train deep networks. As discussed, CRAIG selects a subset that closely approximates the full gradient, and hence can be also applied for speeding up training deep networks. The challenge here is that we cannot use inequality (9) to bound the normed difference between gradients for all x ∈ X and find the subset as a preprocessing step.

However, it has been shown that for neural networks, the variation of the gradient norms is mostly captured by the gradient of the loss w.r.t. the last layer (see Section 3.2 of (Katharopoulos & Fleuret, 2018) ) that is often not expensive or only slightly more expensive than calculating the loss. In many cases, where we have cross entropy loss with soft-max as the last layer, the gradient of the loss w.r. 

4. CONVERGENCE RATE ANALYSIS OF CRAIG:

The idea of CRAIG is to selects a subset that closely approximates the full gradient, and hence can be applied to speed up most IG variants as we show in our experiments. Here, we briefly introduce the original IG method, and then prove the convergence rate of IG applied to subsets found by CRAIG.

4.1. INCREMENTAL GRADIENT METHODS (IG):

Incremental gradient (IG) methods are core algorithms for solving Problem (1) and are widely used and studied. IG aims at approximating the standard gradient method by sequentially stepping along the gradient of the component functions f i in a cyclic order. Starting from an initial point x 1 0 ∈ R d , it makes k passes over all the n components. At every epoch k ≥ 1, it iteratively updates x k i-1 based on the gradient of f i for i = 1, • • • , n using stepsize α k > 0. Formally,

EQUATION

with the convention that x k+1 0 = x k n . Note that for a closed and convex subset X of R d , the results can be projected onto X , and the update rule becomes

EQUATION

where P X denotes projection on the set X ⊂ R d .

IG with diminishing stepsizes converges at rate O(1/ √ k) for strongly convex sum function (Nedić & Bertsekas, 2001) . If in addition to the strong convexity of the sum function, every component function f i is smooth, IG with diminishing stepsizes α k = Θ(1/k s ), s ∈ (0, 1] converges at rate O(1/k s ) (Gürbüzbalaban et al., 2015) .

The convergence rate analysis of IG is valid regardless of order of processing the elements. However, in practice, the convergence rate of IG is known to be quite sensitive to the order of processing the functions (Bertsekas & Scientific, 2015; Gurbuzbalaban et al., 2017) . If problem-specific knowledge can be used to find a favorable order σ (defined as a permutation {σ 1 , • • • , σ n } of {1, 2, ..., n}), IG can be updated to process the functions according to this order, i.e.,

EQUATION

In general a favorable order is not known in advance, and a common approach is sampling the function indices with replacement from the set {1, 2, • • • , n} and is called the Stochastic Gradient Descent (SGD) method, a.k.a. the Robbins-Monro algorithm (Robbins & Monro, 1951 ) (also see (Bottou, 1998; Bertsekas, 2015; Nemirovski et al., 2009; Shalev-Shwartz & Srebro, 2008) ).

4.2. CONVERGENCE RATE OF IG ON SUBSETS FOUND BY CRAIG:

Next we analyze the convergence rate of IG applied to the weighted and ordered subset S found by CRAIG. In particular, we show that (1) applying IG to S converges to a close neighborhood of the optimal solution and that (2) this convergence happens at the same rate (same number of epochs) as IG on the full data. Formally, every step of IG on the subset becomes

EQUATION

Here, σ is a permutation of {1, 2, • • • , r}, and the per-element stepsize γ si for every function f si is the weight of the element s i ∈ S and is fixed for all epochs.

4.3. CONVERGENCE RATE FOR STRONGLY CONVEX FUNCTIONS:

We first provide the convergence analysis for the case where the sum function in Problem (1) is strongly convex, i.e. ∀x, y ∈ R d we have that f (y) ≥ f (x) + ∇f (x), y -x + µ 2 x -y 2 . Theorem 1. Assume that i∈V f i (x) is strongly convex, and S is a weighted subset of size r such that L(S) = i∈V min j∈S d ij ≤ . Then for the iterates {x k = x k 0 } generated by applying IG to S with per-epoch stepsize α k = α/k τ with α > 0 and τ ∈ [0, 1], we have

(i) if τ = 1, then x k -x * 2 ≤ 2 R/µ + r 2 γ 2 max C 2 /kµ, (ii) if 0 < τ < 1, then x k -x * 2 ≤ 2 R/µ, for k → ∞ (iii) if τ = 0, then x k -x * 2 ≤ (1 -αµ) k+1 x 0 -x * 2 + 2 R/µ + αr 2 γ 2 max C 2 /µ,

where C is an upper-bound on the norm of the component function gradients, i.e. max i∈V sup x∈X ∇f i (x) ≤ C, γ max = max j∈S γ j is the largest per-element step size, and R = min{d 0 , (rγ max C + )/µ}, where d 0 = x 0 -x * is the initial distance to the optimum x * .

All the proofs can be found in the Appendix. The above theorem shows that IG on S converges at the same rate O(1/ √ k) of IG on the entire data set V . However, compared to IG on V , the |V |/|S| speedup of IG on S comes at the price of getting an extra error term, 2 R/µ.

4.4. CONVERGENCE RATE FOR SMOOTH AND STRONGLY CONVEX FUNCTIONS:

If in addition to strong convexity of the expected risk, each component function has a Lipschitz gradient, i.e. ∀x ∈ X , i ∈ [n] we have ∇f i (x) -∇f i (y) ≤ β i x -y , then we get the following results about the iterates generated by applying IG to the weighted subset S returned by CRAIG.

Figure 1 : Loss residual and error rate of IG, SVRG, SAGA for Logistic Regression on Covtype data set with 581,012 data points. We compare performance of CRAIG (10% selected subset) vs. entire data set. We achieve the average speedup of 7x for achieving similar loss residual and error rate across the three optimization methods.

Theorem 2. Assume that i∈V f i (x) is strongly convex and let f i (x), i = 1, 2, • • • , n be convex and twice continuously differentiable component functions with Lipschitz gradients on X . Given a subset S such that L(S) = i∈V min j∈S d ij ≤ . Then for the iterates {x k = x k 0 } generated by applying IG to S with per-epoch stepsize α k = α/k τ with α > 0 and τ ∈ [0, 1], we have

(i) if τ = 1, then x k -x * ≤ 2 /µ + βCrγ 2 max /kµ, (ii) if 0 < τ < 1, then x k -x * ≤ 2 /µ, for k → ∞ (iii) if τ = 0, then x k -x * ≤ (1 -αµ) k x 0 -x * + 2 /µ + αβCrγ 2 max /µ,

where β = n i=1 β i is the sum of gradient Lipschitz constants of the component functions.

The above theorem shows that for τ > 0, IG applied to S converges to a 2 /µ neighborhood of the optimal solution, with a rate of O(1/k τ ) which is the same convergence rate for IG on the entire data set V . As shown in our experiments, in real data sets small weighted subsets constructed by CRAIG provide a close approximation to the full gradient. Hence, applying IG to the weighted subsets returned by CRAIG provides a solution of the same or higher quality compared to the solution obtained by applying IG to the whole data set, in a considerably shorter amount of time.

5. EXPERIMENTS:

In our experimental evaluation we wish to address the following questions: (1) How do loss and accuracy of IG applied to the subsets returned by CRAIG compare to loss and accuracy of IG applied to the entire data; (2) How small is the size of the subsets that we can select with CRAIG and still get a comparable performance to IG applied to the entire data;

(3) how does the ordering affect the performance of IG on the subset; and (4) how well does CRAIG scale to large data sets, and extend to non-convex problems. To this end, we apply CRAIG to several convex and non-convex problems. In our experiments, we report the run-time as the wall-clock time for subset selection with CRAIG, plus minimizing the loss using IG or other optimizers with the specified learning rates. For the classification problems, we separately select subsets from each class while maintaining the class ratios in the whole data, and apply IG to the union of the subsets.

5.1. CONVEX EXPERIMENTS:

In our convex experiments, we apply CRAIG to IG, as well as variance reduction methods SVRG (Johnson & Zhang, 2013) , and SAGA (Defazio et al., 2014) , that try to reduce the variance of SGD either based on computations of full gradients at pivot points, or on keeping per data point corrections Figure 1 compares training loss residual and test error rate of IG, SVRG, and SAGA on the subsets of size 10% of covtype selected by CRAIG (with corresponding per-element stepsizes) to that of IG, SVRG, and SAGA on the entire data set. We used a constant learning rate of α for SAGA and SVRG, and α/ √ k for IG, where α = 10 -3 for Covtype and Ijcnn, and 3 × 10 -5 for SensIT. It can be seen that subsets obtained by CRAIG achieve a similar loss and error rate as the entire data sets, but much faster. In particular, we obtained a speedup of 8.3x, 8.3x, 4.5x from applying IG, SVRG and SAGA on the subsets of size 10% obtained by CRAIG.

Figure 2 top row, compares the speedup achieved by CRAIG to reach a similar loss residual as that of IG after 50 epochs for subsets of size 5% to 25%. The bottom row compares the L2-norm of the distance to optimal solution (estimated by running IG for a long time) for IG applied to the subsets of various size obtained by CRAIG. We observe that while smaller subsets provide a larger speedup, IG on smaller subsets may converge to a point farther away from the optimal solution.

Finally, Figure 3 shows the loss residual vs time for IG when it processes the elements of the subsets according to the ordering obtained by CRAIG compared to random permutations of the same subsets. We observe that the greedy ordering significantly improves the rate of convergence of IG. We apply CRAIG to several popular methods for training neural networks, including SGD, Nesterov Accelerated Gradient (NAG) method (Nesterov, 1983) , Adagrad (Duchi et al., 2011) , and Adaptive Moment Estimation (Adam) (Kingma & Ba, 2014) . Momentum is a method that helps accelerate SGD in the relevant direction and dampens oscillations. NAG adapts momentum updates to the slope of our error function and speed up SGD. Adagrad adapts the learning rate to the parameters, performing smaller updates (i.e. low learning rates) for parameters associated with frequently occurring features, and larger updates (i.e. high learning rates) for parameters associated with infrequent features. Adam computes individual adaptive learning rates for different parameters from estimates of first and second moments of the gradients. Fig. 4a shows loss, and error rate for training a 1-layer neural net on MNIST. Fig. 4b show similar quantities for training ResNet-56 on CIFAR10. For both problems, we used a constant learning rate of 10 -2 . Here, we apply CRAIG to select a subset of 30%-40% of the data at the beginning of every epoch and train only on the selected subset with the corresponding per-element stepsizes. Interestingly, in addition to achieving a speedup of 2x to 3x for training neural networks, the subsets selected by CRAIG provide a better generalization performance compared to models trained on the entire data set.

6. CONCLUSION:

We developed a method, CRAIG, for selecting an ordered subset (coreset) of data points with their corresponding per-element stepsizes to speed up iterative gradient (IG) methods. In particular, we showed that weighted subsets that minimize the upper-bound on the estimation error of the full gradient, maximize a submodular facility location function. Hence, we can obtain an ordered subset of data points with their corresponding learning rates using a greedy algorithm. We showed that IG on subsets S returned by CRAIG converges at the same rate as IG on the entire data set V , while providing a |V |/|S| speedup. In our set of experiments, we showed that various IG methods, including SAGA, SVRG, NAG, Adagrad and Adam runs up to 10x faster on convex and up to 3x on non-convex problems on subsets found by CRAIG while achieving practically the same trining loss and test error. Finally, we empirically demonstrated the effect of the ordering found by the greedy algorithm on the convergence rate of IG methods.

A CONVERGENCE RATE ANALYSIS:

We firs proof the following Lemma which is an extension of the [Chung et al. (1954) , Lemma 4]. Lemma 3. Let u k ≥ 0 be a sequence of real numbers. Assume there exist k 0 such that

u k+1 ≤ (1 - c k )u k + e k p + d k p+1 , ∀k ≥ k 0 ,

where e > 0, d > 0, c > 0 are given real numbers. Then

u k ≤ (dk -1 + e)(c -p + 1) -1 k -p+1 + o(k -p+1 ) for c > p -1, p ≥ 1 (18) u k = O(k -c log k) for c = p -1, p > 1 (19) u k = O(k -c ) for c < p -1, p > 1 (20) (21) Proof. Let c > p -1 and v k = k p-1 u k - d k(c-p+1

) -e c-p+1 . Then, using Taylor approximation

EQUATION

EQUATION

EQUATION

Note that for v k , we have

∞ k=0 1 - c -p + 1 k + o( 1 k ) = ∞ and d/(c -p + 1) k(k + 1) + e(p -1) k 2 + d(p -1) k 3 + o( 1 k 2 ) 1 - c -p + 1 k + o( 1 k ) -1 → 0.

Therefore, lim k→∞ v k ≤ 0, and we get Eq. 18. For p = 1, we have u k ≤ e c . Hence, u k converges into the region u ≤ e c , with ratio 1 -c k . Moreover, for p -1 ≥ c we have

EQUATION

for sufficiently large k. Summing over k, we obtain that v k is bounded for p -1 > c (since the series

∞ k=1 (1/k α ) converges for α > 1) and v k = O(log k) for p = c + 1 (since k i=1 (1/i) = O(log k)).

Assuming gradients have a bounded norm max x∈X , j∈V ∇f j (x) ≤ C, and the fact that j∈S γ j = n we can write

EQUATION

Thus for initial distance x 0 -x * = d 0 , we have

EQUATION

46) Putting Eq. 42 to Eq. 46 together we get

EQUATION

)

Now, from convexity of every f j for j ∈ S we have that

EQUATION

In addition, incremental updates gives us

EQUATION

Therefore, we get

2α k j∈S (f j (x k ) -f j (x k j-1 )) + α 2 k rγ 2 max C 2 ≤ 2α k r i=1 γ max C • α k (j -1)γ max C + α 2 k rγ 2 max C 2 (50) = α 2 k r 2 γ 2 max C 2 (51) Hence, x k+1 -x * 2 ≤ (1 -α k µ) x k -x * 2 + 2α k R/µ + α 2 k r 2 γ 2 max C 2 .

(52) where γ max is the size of the largest cluster, and C is the upperbound on the gradients. For 0 < s ≤ 1, the theorem follows by applying Lemma 3 to Eq. 52, with c = µ, e = 2 R/µ, and d = r 2 γ 2 max C 2 . For s = 0, where we have a constant step size α

EQUATION

) and therefore,

EQUATION

A.2 CONVERGENCE STRONGLY CONVEX AND SMOOTH COMPONENT FUNCTIONS PROOF OF THEOREM 2

IG updates for cycle k on subset S can be written as

EQUATION

EQUATION

Building on the analysis of Gürbüzbalaban et al. (2015) , for convex and twice continuously differentiable function, we can write

EQUATION

where

A r k = 1 0 ∇ 2 f (x * + τ (x k -x * )

)dτ is average of the Hessian matrices corresponding to the r (weighted) elements of S on the interval [x k , x * ].

From Eq. 58 we have i∈V

EQUATION

where A k is average of the Hessian matrices corresponding to all the n component functions on the interval [x k , x * ]. Taking norm of both sides and noting that i∈V f i (x * ) = 0 and hence j∈S γ j f j (x * ) ≤ , we get

EQUATION

where is the estimation error of the full gradient by the weighted gradients of the elements of the subset S, and we used i∈V ∇f i (x k ) -j∈S γ j ∇f j (x k ) ≤ . Substituting Eq. 58 into Eq. 56 we obtain

EQUATION

Taking norms of both sides, we get

EQUATION

Now, we have

EQUATION

Substituting into Eq. 62, we obtain

EQUATION

From strong convexity of i∈V f i (x), and gradient smoothness of each component f i (x) we have 

EQUATION

EQUATION

where in the last line we used |S| = r. Therefore,

EQUATION

For 0 < s ≤ 1, the theorem follows by applying Lemma 3 to Eq. 73 with c = µ, e = 2 , d = βCrγ 2 max . For s = 0, where we have a constant step size α k = α ≤ 1 β , we get

x k+1 -x * ≤ (1 -αµ) k+1 x k -x * + 2 k i=0 (1 -αµ) i + α 2 k i=0

(1 -αµ) i βCrγ 2 max (75)

EQUATION

EQUATION

where the inequality in Eq. 76 follows since k i=0 (1 -αµ) i ≤ 1 αµ .

B NORM OF THE DIFFERENCE BETWEEN GRADIENTS:

For ridge regression f i (x) = 1 2 ( a i , x -y i ) 2 + λ 2 x 2 , we have ∇f i (x) = a i ( a i , x -y i ) + λx. Therefore, ∇f i (x) -∇f j (x) = ( a i -a j . x + y i -y j ) a j (78)

For a i ≤ 1, and y i = y j we get

EQUATION

For reguralized logistic regression with y ∈ {-1, 1}, we have ∇f i (x) = y i /(1 + e yi ai,x ). For y i = y j we get ∇f i (x) -∇f j (x) = e ai-aj . x -1 1 + e -ai,x a j .

For a i ≤ 1, using Taylor approximation e x ≤ 1 + x, and noting that 1 1+e -a i ,x ≤ 1 we get ∇f i (x) -∇f j (x) ≤ a i -a j . x 1 + e -ai,x a j ≤ a i -a j O( x ).

(81)

For classification, we require y i = y j , hence we can select subsets from each class and then merge the results. On the other hand, in ridge regression we also need |y i -y j | to be small. Similar results can be deduced for other loss functions including square loss, smoothed hinge loss, etc.

Assuming x is bounded for all x ∈ X , upper-bounds on the euclidean distances between the gradients can be pre-computed. 

C ADDITIONAL EXPERIMENTS:



:

Note that in the worst case we may need |S * | ≈ |V | to approximate the gradient. However, as we show in experiments, in practice we find that a small subset is sufficient to accurately approximate the gradient.

Figures and Tables:

Figure 2: (top) Speedup of CRAIG applied to to get similar loss residual as IG after 50 epoch, and (bottom) distance to the optimal solution vs various subset sizes for (left) Covtype, (middle) SensIT, and (right) Ijcnn1. Smaller subsets provides larger speedups, but may converge farther away from the optimal solution.

Figure 3: Loss residual of CRAIG for Logistic Regression on (left) Covtype, (middle) SensIT, and (right) Ijcnn, where we process the points of the subset according to the ordering provided by CRAIG vs three random permutations of the same subset. Notice convergence is significantly faster when we process the points in CRAIG order.

Figure 4: Training loss and test accuracy of CRAIG applied to (a) SGD on MNIST with a 1-layer neural network, and (b) SGD, Adam, Adagrad, NAG, on CIFAR-10 with ResNet-56. CRAIG provides 2x to 3x speedup and a better generalization performance.

where β = i∈V β i In addition, from the gradient smoothness of the components we can write e k ≤ j∈S γ j β j x k -

Figure 5: Loss residual vs. time for IG on 10 subsets of size 10%, 20%, 30%, ..., 100% selected by CRAIG, random subsets, and random subsets weighted by |V |/|S|. Stepsizes are tuned for every subset separately, by preferring smaller training loss from a large number of parameter combinations for two types of learning scheduling: exponential decay η(t) = η 0 a t/n with parameters η 0 and a to adjust and t-inverse η(t) = η 0 (1 + b t/n ) -1 with η 0 and b to adjust.

Figure 6: Training loss and test accuracy for SGD applied to full MNIST vs. subsets of size 60% selected by CRAIG and random subsets of size 60%. Both the random subsets and the subsets found by CRAIG change at the beginning of every epoch.

Peer review (iclr2020):
The paper proposes a theoretically founded method to generate subsets of a dataset, together with corresponding sample weights in a way that the average gradient of the subset is at most epsilon far from the average gradient of the full dataset. Given such a subset, the authors provide theoretical guarantees for convergence to an epsilon neighborhood of the optimum for strongly convex functions. The proposed algorithm to create such a subset is a greedy algorithm that relies on parameter independent similarities between samples, namely similarity scores that are true regardless of the current value of the function parameters. Although I find the approach interesting, I have three main concerns with the proposed method. 1. The experimental setup is lacking significant information, baselines and baseline tuning (see below for more in depth comments). 2. The proposed upper bound which has been used for a similar purpose by [1] becomes nonsensical in high dimensions and although for [1] this would mean sampling with a non optimal sampling distribution for CRAIG it means converging very far from the optimum. What are the values of epsilon that you observe in practice? 3. I do not see how CRAIG would be applied to deep learning. The argument in section 3.4 is that the variance of the gradient norm is captured by the gradient of the last layer or last few layers, however this is true given the parameters of the neural network. The gradients can change arbitrarily after a very small number of parameter updates as shown by [2]. Experimental setup ---------------------------- For the case of the convex problems, the learning rate is not tuned independently for each method. Even more importantly the stepsizes of CRAIG are all numbers larger than 1 so the expected learning rate is multiplied by the average step size. This makes it difficult to understand whether the speedup is due to a larger learning rate or due to CRAIG. Similarly for figure 3 the result could be due to a non decreasing step size because of gamma_j while for CRAIG gamma_j are ordered in decreasing order. In addition, there is no experimental analysis of the epsilon bound and the actual difference of the gradients for the subset and the full dataset. There are also no baselines that use a subset to train. A comparison with a baseline that uses 1. random subset or 2. a subset selected via importance sampling from [1] would contribute towards understanding the particular benefits of CRAIG. Regarding the neural network experiments: 1. There is no explicit definition of the similarity function used for the case of neural networks. If we assume based on 3.4 that the algorithm requires an extra forward pass in the beginning of every epoch there should be visible steps in Figure 3 where time passes but the loss doesn*t move. 2. 2000 seconds and 80% accuracy on MNIST points towards a mistake on the implementation of the training. On a laptop CPU it takes ~15s per epoch and achieves ~95% test accuracy from the first epoch for the neural network described. 3. Similarly 80% accuracy on CIFAR10 is sufficiently low for Resnet-56 to be alarming. [1] Zhao, Peilin, and Tong Zhang. *Stochastic optimization with importance sampling for regularized loss minimization.* international conference on machine learning. 2015. [2] Defazio, Aaron, and Léon Bottou. *On the ineffectiveness of variance reduced optimization for deep learning.* arXiv preprint arXiv:1812.04529 (2018).