Year: 2019

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.

heading: 1 INTRODUCTION
text: 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 :
x∗ ∈ argmin x∈X f(x), f(x) := ∑ i∈V fi(x) + r(x), fi(x) = l(x, (ai, yi)), (1)
where V = {1, . . . , n} is an index set of the training data, and functions fi : Rd → R are associated with training examples (ai, yi), where ai ∈ Rd is the feature vector, and yi 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{d0, (rγmaxC + )/µ} where d0 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.
We demonstrate the effectiveness of CRAIG via an extensive set of experiments using logistic regression (a convex optimization problem) as well as training neural networks (non-convex optimization problems). We show that CRAIG speeds up incremental gradient methods, including SGD, SAGA, SVRG, Adam, Adagrad, and NAG. In particular, CRAIG while achieving practically the same loss
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.

heading: 2 RELATED WORK
text: 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, Nedić & 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/k2τ ) 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).

heading: 3 CORESETS FOR INCREMENTAL GRADIENT DESCENT (CRAIG)
text: 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 ∇fi(x) over V by iteratively making a step based on the gradient of every function fi. 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
S∗ = argmin S⊆V,γj≥0 ∀j |S|, s.t. max x∈X ‖ ∑ i∈V ∇fi(x)− ∑ j∈S γj∇fj(x)‖ ≤ . (2)
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.
1Note 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.
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 fi 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 ∇fi(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.

heading: 3.1 UPPER-BOUND ON THE ESTIMATION ERROR
text: We first address Problem 1, i.e., how to quickly estimate the error/discrepancy of the weighted sum of gradients of functions fj 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 Cj = {i ∈ [n]|ς(i) = j} ⊆ V be the set of data points that are assigned to j ∈ S, and γj = |Cj | be the number of such data points. Hence, {Cj}rj=1 form a partition of V . Then, for any arbitrary (single) x ∈ X we can write∑
i∈V ∇fi(x) = ∑ i∈V ( ∇fi(x)−∇fςx(i)(x) +∇fς(i)(x) ) (3)
= ∑ i∈V ( ∇fi(x)−∇fςx(i)(x) ) + ∑ j∈S γj∇fj(x). (4)
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 functionsfj forj∈S. I.e.,
‖ ∑ i∈V ∇fi(x)− ∑ j∈S γj∇fj(x)‖ ≤ ∑ i∈V ‖∇fi(x)−∇fςx(i)(x)‖, (5)
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: ςx(i)∈argminj∈S‖∇fi(x)−∇fj(x)‖. Hence,
min S⊆V ‖ ∑ i∈V ∇fi(x)− ∑ j∈S γj∇fj(x)‖ ≤ ∑ i∈V min j∈S ‖∇fi(x)−∇fj(x)‖. (6)
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 dij between gradients of fi and fj as the maximum normed difference between ∇fi(x) and ∇fj(x) over all x ∈ X :
dij , max x∈X ‖∇fi(x)−∇fj(x)‖. (7)
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
S∗ = argmin S⊆V |S|, such that L(S) , ∑ i∈V min j∈S dij ≤ . (8)
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):
∀x, i, j ‖∇fi(x)−∇fj(x)‖ ≤ dij ≤ max x∈X O(‖x‖) · ‖ai − aj‖ = const. ‖ai − aj‖. (9)
Note when ‖x‖ is bounded for all x ∈ X , i.e., maxx∈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 fi for i ∈ V = {1, · · · , n}. Output: Subset S ⊆ V with corresponding per-element stepsizes {γ}j∈S , and an ordering σ.
1: S0 ← ∅, s0 = 0, i = 0. 2: while F (S) < L({s0})− do 3: j ∈ argmaxe∈V \Si−1 F (e|Si−1) 4: Si = Si−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 = argmins∈S maxx∈X ‖∇fi(x)−∇fs(x)‖
] 10: end for

heading: 3.2 THE CRAIG ALGORITHM
text: 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,
S∗ = argmin S⊆V |S|, such that F (S) ≥ ρ. (10)
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 S0 = ∅, and at each iteration i, it chooses an element e ∈ V that maximizes 4(e|Si−1), i.e., Si = Si−1 ∪ {argmaxe∈V F (e|Si−1) − F (Si−1)}. Greedy gives us a logarithmic approximation, i.e. |S| ≤ ( 1 + ln(maxe 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 Ffl(S) = ∑ i∈V maxj∈S si,j , where si,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 s0 we can turn L(S) in Eq. (8) into a monotone submodular facility location function, F (S) = L({s0})− L(S ∪ {s0}), (11) where L({s0}) 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 s0, 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:
S∗ = argmin S⊆V |S|, such that F (S) ≥ L({s0})− . (12)
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.

heading: 3.3 ORDERING ON THE ELEMENTS OF THE SUBSET
text: 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 ∈ argmaxe∈V F (e|Si−1) improves our estimation from the full gradient by an amount bounded by the marginal gain. Formally, at every step i, we have F (Si) ≥ (1− e−i/|S|)F (S∗), and hence
‖ ∑ i∈V ∇fi(x)− ∑ j∈S γj∇fj(x)‖ ≤ cnt− (1− e−i/|S|)L(S∗). (13)
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.

heading: 3.4 APPLICATION OF CRAIG TO DEEP NETWORKS
text: 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.t. the i-th input to the soft-max is simply pi − yi, where pi are logits (dimension p− 1 for p classes) and y is the one-hot encoded label. In this case, CRAIG does not need any backward pass or extra storage. That is, CRAIG can be applied at the beginning of every epoch to find a subset for that epoch. Note that, although CRAIG needs an additional O(|V | · |S|) complexity (or O(|V |) using stochastic greedy) to find the subset S at the beginning of every epoch, this complexity does not involve any (exact) gradient calculations and is negligible compared to the cost of backpropagations performed during the epoch. Hence, as we show in the experiments CRAIG is practical and salable.

heading: 4 CONVERGENCE RATE ANALYSIS OF CRAIG
text: 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.

heading: 4.1 INCREMENTAL GRADIENT METHODS (IG)
text: 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 fi in a cyclic order. Starting from an initial point x10 ∈ Rd, it makes k passes over all the n components. At every epoch k ≥ 1, it iteratively updates xki−1 based on the gradient of fi for i = 1, · · · , n using stepsize αk > 0. Formally,
xki = x k i−1 − αk∇fi(xki−1), i = 1, 2, · · · , n, (14)
with the convention that xk+10 = x k n. Note that for a closed and convex subset X of Rd, the results can be projected onto X , and the update rule becomes xki = PX (x k i−1 − αk∇fi(xki−1)), i = 1, 2, · · · , n, (15)
where PX denotes projection on the set X ⊂ Rd. IG with diminishing stepsizes converges at rate O(1/ √ k) for strongly convex sum function (Nedić & Bertsekas, 2001). If in addition to the strong convexity of the sum function, every component function fi is smooth, IG with diminishing stepsizes αk = Θ(1/ks), s ∈ (0, 1] converges at rate O(1/ks) (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.,
xki = x k i−1 − αk∇fσi(xki−1), i = 1, 2, · · · , n. (16)
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)).

heading: 4.2 CONVERGENCE RATE OF IG ON SUBSETS FOUND BY CRAIG
text: 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
xki = x k i−1 − αkγsσi∇fsσi (x k i−1), i = 1, 2, · · · , r, si ∈ S, |S| = r. (17)
Here, σ is a permutation of {1, 2, · · · , r}, and the per-element stepsize γsi for every function fsi is the weight of the element si ∈ S and is fixed for all epochs.

heading: 4.3 CONVERGENCE RATE FOR STRONGLY CONVEX FUNCTIONS
text: We first provide the convergence analysis for the case where the sum function in Problem (1) is strongly convex, i.e. ∀x, y ∈ Rd we have that f(y) ≥ f(x) + 〈∇f(x), y − x〉+ µ2 ‖x− y‖ 2.
Theorem 1. Assume that ∑ i∈V fi(x) is strongly convex, and S is a weighted subset of size r such
that L(S) = ∑ i∈V minj∈S dij ≤ . Then for the iterates {xk=xk0} generated by applying IG to S with per-epoch stepsize αk = α/kτ with α > 0 and τ ∈ [0, 1], we have
(i) if τ = 1, then ‖xk − x∗‖2 ≤ 2 R/µ+ r2γ2maxC2/kµ,
(ii) if 0 < τ < 1, then ‖xk − x∗‖2 ≤ 2 R/µ, for k →∞
(iii) if τ = 0, then ‖xk − x∗‖2 ≤ (1− αµ)k+1‖x0 − x∗‖2 + 2 R/µ+ αr2γ2maxC2/µ,
where C is an upper-bound on the norm of the component function gradients, i.e. maxi∈V supx∈X ‖∇fi(x)‖ ≤ C, γmax = maxj∈S γj is the largest per-element step size, and R = min{d0, (rγmaxC + )/µ}, where d0 = ‖x0 − 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/µ.

heading: 4.4 CONVERGENCE RATE FOR SMOOTH AND STRONGLY CONVEX FUNCTIONS
text: 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 ‖∇fi(x)−∇fi(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.
Theorem 2. Assume that ∑ i∈V fi(x) is strongly convex and let fi(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 minj∈S dij ≤ . Then for the iterates {xk = xk0} generated by applying IG to S with per-epoch stepsize αk = α/kτ with α > 0 and τ ∈ [0, 1], we have
(i) if τ = 1, then ‖xk − x∗‖ ≤ 2 /µ+ βCrγ2max/kµ,
(ii) if 0 < τ < 1, then ‖xk − x∗‖ ≤ 2 /µ, for k →∞
(iii) if τ = 0, then ‖xk − x∗‖ ≤ (1− αµ)k‖x0 − x∗‖+ 2 /µ+ αβCrγ2max/µ,
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.

heading: 5 EXPERIMENTS
text: 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.

heading: 5.1 CONVEX EXPERIMENTS
text: 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
in memory. We apply L2-regularized logistic regression: fi(x) = ln(1 + exp(−xTaiyi)) + 0.5λxTx to classify the following three datasets from LIBSVM: (1) covtype.binary including 581,012 data points of 54 dimensions, (2) SensIT including 78,823 training and 19,705 test data points of dimension 100, and (3) Ijcnn1 including 49,990 training and 91,701 test data points of 22 dimensions. As covtype do not come with labeled test data, we randomly split the training data into halves to make the training/test split (training and set sets are consistent for different methods). For the convex experiments, λ is set to 10−5.
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.

heading: 5.2 NON-CONVEX EXPERIMENTS
text: Our non-convex experiments involves applying CRAIG to train the following two neural networks: (1) Our smaller network is a fully-connected hidden layer of 100 nodes and ten softmax output nodes; sigmoid activation and L2 regularization with λ = 0.0001 and mini-batches of size 10 on MNIST dataset of handwritten digits containing 60,000 training and 10,000 test images. (2) Our large neural network is ResNet-56 for CIFAR10 with convolution, average pooling and dense layers with softmax outputs and L2 regularization with λ = 2× 10−4 CIFAR 10 includes 50,000 training and 10,000 test images from 10 classes, and we used mini-batches of size 128. Both MNIST and CIFAR10 data sets are normalized into [0, 1] by division with 255.
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.

heading: 6 CONCLUSION
text: 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.

heading: A CONVERGENCE RATE ANALYSIS
text: We firs proof the following Lemma which is an extension of the [Chung et al. (1954), Lemma 4]. Lemma 3. Let uk ≥ 0 be a sequence of real numbers. Assume there exist k0 such that
uk+1 ≤ (1− c
k )uk +
e
kp +
d
kp+1 , ∀k ≥ k0,
where e > 0, d > 0, c > 0 are given real numbers. Then
uk ≤ (dk−1 + e)(c− p+ 1)−1k−p+1 + o(k−p+1) for c > p− 1, p ≥ 1 (18) uk = O(k
−c log k) for c = p− 1, p > 1 (19) uk = O(k
−c) for c < p− 1, p > 1 (20) (21)
Proof. Let c > p − 1 and vk = kp−1uk − dk(c−p+1) − e
c−p+1 . Then, using Taylor approximation (1 + 1k ) p = (1 + pk ) + o( 1 k ) we can write
vk+1 = (k + 1) p−1uk+1 −
d (k + 1)(c− p+ 1) − e c− p+ 1
(22)
≤ kp−1(1 + 1 k
)p−1 (
(1− c k )uk + e kp + d kp+1
) − d
(k + 1)(c− p+ 1) − e c− p+ 1
(23)
= kp−1uk ( 1− c− p+ 1
k + o(
1 k ) ) + e k ( 1 + p− 1 k + o( 1 k ) )
(24)
+ d
k2
( 1 +
p− 1 k + o( 1 k ) ) − d (k + 1)(c− p+ 1) − e c− p+ 1
(25)
= ( vk +
d
k(c− p+ 1) +
e
c− p+ 1
)( 1− c− p+ 1
k + o(
1 k ) )
(26)
+ e
k
( 1 +
p− 1 k + o( 1 k ) ) + d k2 ( 1 + p− 1 k + o( 1 k ) )
(27)
− d (k + 1)(c− p+ 1) − e c− p+ 1
(28)
= vk
( 1− c− p+ 1
k + o(
1 k ) ) + d/(c− p+ 1) k(k + 1) + e(p− 1) k2 + d(p− 1) k3 + o( 1 k2 ) (29)
Note that for vk, we have ∞∑ k=0 ( 1− c− p+ 1 k + o( 1 k ) ) =∞ and (d/(c− p+ 1) k(k + 1) + e(p− 1) k2 + d(p− 1) k3 + o( 1 k2 ) )( 1− c− p+ 1 k + o( 1 k ) )−1 → 0.
Therefore, limk→∞ vk ≤ 0, and we get Eq. 18. For p = 1, we have uk ≤ ec . Hence, uk converges into the region u ≤ ec , with ratio 1− c k .
Moreover, for p− 1 ≥ c we have
vk+1 = uk+1(k + 1) c ≤ [ (1− c
k )uk +
e
kp +
d
kp+1
] kc ( 1 + c
k +
c2
2k2 + o(
1 k2 ) )
(30)
= ( 1− c 2
2k2 + o(
1 k2 ) ) vk +
d
kp−c+1
( 1 +O( 1 k ) ) + e kp−c ( 1 + c k +O( 1 k2 ) )
(31)
≤ vk + e′
kp−c (32)
for sufficiently large k. Summing over k, we obtain that vk is bounded for p − 1 > c (since the series ∑∞ k=1(1/k α) converges for α > 1) and vk = O(log k) for p = c + 1 (since ∑k i=1(1/i) = O(log k)).
In addition, based on [Chung et al. (1954), Lemma 5] for uk ≥ 0, we can write
uk+1 ≤ (1− c
ks )uk +
e
kp +
d kt , 0 < s < 1, s ≤ p < t. (33)
Then, we have
uk ≤ e
c
1
kp−s + o(
1
kp−s ). (34)
A.1 CONVERGENCE RATE FOR STRONGLY CONVEX FUNCTIONS
PROOF OF THEOREM 1
We now provide the convergence rate for strongly convex functions building on the analysis of Nedić & Bertsekas (2001). For non-smooth functions, gradients can be replaced by sub-gradients.
Let xk = xk0 . For every IG update on subset S we have
‖xkj − x∗‖2 = ‖xkj−1 − αkγj∇fj(xkj−1)− x∗‖2 (35) = ‖xkj−1 − x∗‖2 − 2αkγj∇fj(xkj−1)(xkj−1 − x∗) + α2k‖γj∇fj(xkj−1)‖2 (36) ≤ ‖xkj−1 − x∗‖2 − 2αk(fj(xkj−1)− fi(x∗)) + α2k‖γi∇fj(xkj−1)‖2. (37)
Adding the above inequalities over elements of S we get ‖xk+1 − x∗‖2 ≤ ‖xk − x∗‖2 − 2αk ∑ j∈S (fi(x k j−1)− fj(x∗)) + α2k ∑ j∈S ‖γj∇fj(xkj−1)‖2 (38)
= ‖xk − x∗‖2 − 2αk ∑ j∈S (fj(xk)− fi(x∗))
+ 2αk ∑ j∈S (fj(x k j−1)− fj(xk)) + α2k ∑ j∈S ‖γj∇fj(xkj−1)‖2 (39)
Using strong convexity we can write ‖xk+1 − x∗‖2 ≤ ‖xk − x∗‖2 − 2αk (∑ j∈S γj∇fj(x∗) · (xk − x∗) + µ 2 ‖xk − x∗‖2 ) + 2αk
∑ j∈S (fj(x k j−1)− fj(xk)) + α2k ∑ j∈S ‖γj∇fj(xkj−1)‖2 (40)
Using Cauchy–Schwarz inequality, we know | ∑ j∈S γj∇fj(x∗) · (xk − x∗)| ≤ ‖ ∑ j∈S γj∇fj(x∗)‖ · ‖xk − x∗‖. (41)
Hence,
− ∑ j∈S γj∇fj(x∗) · (xk − x∗) ≤ ‖ ∑ j∈S γj∇fj(x∗)‖ · ‖xk − x∗‖. (42)
From reverse triangle inequality, and the facts that S is chosen in a way that ‖ ∑ i∈V ∇fi(x∗) −∑
j∈S γj∇fj(x∗)‖ ≤ , and that ∑ i∈V ∇fi(x∗) = 0 we have ‖ ∑ j∈S γj∇fj(x∗)‖ ≤
‖ ∑ i∈V ∇fi(x∗)‖+ = . Therefore
‖ ∑ j∈S γj∇fj(x∗)‖ · ‖xk − x∗‖ ≤ · ‖xk − x∗‖ (43)
For a continuously differentiable function, the following condition is implied by strong convexity condition
‖xk − x∗‖ ≤ 1 µ ‖ ∑ j∈S γj∇fj(xk)‖. (44)
Assuming gradients have a bounded norm maxx∈X , j∈V ‖∇fj(x)‖ ≤ C, and the fact that
∑ j∈S γj = n
we can write ‖ ∑ j∈S γj∇fj(xk)‖ ≤ n · C. (45)
Thus for initial distance ‖x0 − x∗‖ = d0, we have ‖xk − x∗‖ ≤ min(n · C, d0) = R (46)
Putting Eq. 42 to Eq. 46 together we get ‖xk+1 − x∗‖2 ≤ (1− αkµ)‖xk − x∗‖2 + 2αk R/µ
+ 2αk ∑ j∈S (fj(xj−1,k)− fj(xk)) + α2krγ2maxC2. (47)
Now, from convexity of every fj for j ∈ S we have that fj(xk)− fj(xkj−1) ≤ ‖γj∇fj(xk)‖ · ‖xkj−1 − xk‖. (48)
In addition, incremental updates gives us
‖xkj−1 − xk‖ ≤ αk j−1∑ i=1 ‖γi∇fi(xki−1)‖ ≤ αk(j − 1)γmaxC. (49)
Therefore, we get 2αk ∑ j∈S (fj(xk)− fj(xkj−1)) + α2krγ2maxC2
≤ 2αk r∑ i=1 γmaxC · αk(j − 1)γmaxC + α2krγ2maxC2 (50)
= α2kr 2γ2maxC 2 (51) Hence, ‖xk+1 − x∗‖2 ≤ (1− αkµ)‖xk − x∗‖2 + 2αk R/µ+ α2kr2γ2maxC2. (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 = r2γ2maxC 2.
For s = 0, where we have a constant step size αk = α ≤ 1µ , we get
‖xk+1 − x∗‖2 ≤ (1−αµ)k+1‖x0 − x∗‖2
+ 2 α R k∑ j=0 (1− αµ)j/µ+ α2r2γ2maxC2 k∑ j=0 (1− αµ)j (53)
Since ∑k j=0(1− αµ)j ≤ 1 αµ , we get
‖xk+1 − x∗‖2 ≤ (1− αµ)k+1‖x0 − x∗‖2 + 2 R/µ+ αr2γ2maxC2/µ, (54) and therefore,
‖xk+1 − x∗‖2 ≤ (1− αµ)k+1‖xk − x∗‖2 + 2 R/µ+ αr2γ2maxC2/µ. (55)
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 xk+1 = xk − αk( ∑ j∈S γj∇fj(xk)− ek) (56)
ek = ∑ j∈S γi(∇fj(xk)−∇fj(xkj−1)) (57)
Building on the analysis of Gürbüzbalaban et al. (2015), for convex and twice continuously differentiable function, we can write∑
j∈S γj∇fj(xk)− ∑ j∈S γj∇fj(x∗) = Ark(xk − x∗) (58)
where Ark = ∫ 1 0 ∇2f(x∗ + τ(xk − x∗))dτ is average of the Hessian matrices corresponding to the r (weighted) elements of S on the interval [xk, x∗].
From Eq. 58 we have∑ i∈V (∇fi(xk)−∇fi(x∗))− ∑ j∈S γj(∇fj(xk)−∇fj(x∗)) = Ak(xk − x∗)−Ark(xk − x∗), (59)
where Ak is average of the Hessian matrices corresponding to all the n component functions on the interval [xk, x∗]. Taking norm of both sides and noting that ∑ i∈V fi(x∗) = 0 and hence
‖ ∑ j∈S γjfj(x∗)‖ ≤ , we get
‖(Ak −Ark)(xk − x∗)‖ = ‖ (∑ i∈V ∇fi(xk)− ∑ j∈S γj∇fj(xk) ) + ∑ j∈S γjfj(x∗)‖ ≤ 2 , (60)
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 ∇fi(xk)− ∑ j∈S γj∇fj(xk)‖ ≤ .
Substituting Eq. 58 into Eq. 56 we obtain
xk+1 − x∗ = (I − αkArk)(xk − x∗) + αkek (61)
Taking norms of both sides, we get
‖xk+1 − x∗‖ ≤ ‖(I − αkArk)(xk − x∗‖) + αk‖ek‖ (62)
Now, we have
‖(I − αkArk)(xk − x∗)‖ = ‖I(xk − x∗)− αkArk(xk − x∗)‖ (63) = ‖I(xk − x∗)− αk(Ark −Ak)(xk − x∗)− αkAk(xk − x∗)‖ (64) ≤ ‖(I − αkAk)(xk − x∗)‖+ αk‖(Ak −Ark)(xk − x∗)‖ (65) ≤ ‖(I − αkAk)(xk − x∗)‖+ 2αk (66)
Substituting into Eq. 62, we obtain
‖xk+1 − x∗‖ ≤ ‖I − αkAk‖ · ‖xk − x∗‖+ 2αk + αk‖ek‖ (67)
From strong convexity of ∑ i∈V fi(x), and gradient smoothness of each component fi(x) we have
µIn ∑ i∈V ∇2fi(x), Ak βIn, x ∈ X , (68)
where β = ∑ i∈V βi In addition, from the gradient smoothness of the components we can write
‖ek‖ ≤ ∑ j∈S γjβj‖xk − xkj ‖ (69)
≤ ∑ j∈S γjβj j−1∑ i=1 ‖xki−1 − xki ‖ (70) ≤ ∑ j∈S γjβjαk j−1∑ i=1 ‖γi∇fi(xki )‖ (71)
≤ αkβCrγ2max, (72)
where in the last line we used |S| = r. Therefore,
‖xk+1 − x∗‖ ≤ max(‖1− αkµ‖, ‖1− αkβ‖)‖xk − x∗‖+ 2αk + α2kβCrγ2max (73) ≤ (1− αkµ)‖xk − x∗‖+ 2αk + α2kβCrγ2max if αkβ ≤ 1. (74)
For 0 < s ≤ 1, the theorem follows by applying Lemma 3 to Eq. 73 with c = µ, e = 2 , d = βCrγ2max. For s = 0, where we have a constant step size αk = α ≤ 1β , we get
‖xk+1 − x∗‖ ≤ (1− αµ)k+1‖xk − x∗‖+ 2 k∑ i=0 (1− αµ)i + α2 k∑ i=0 (1− αµ)iβCrγ2max (75)
≤ (1− αµ)k+1‖xk − x∗‖+ 2 /µ+ αβCrγ2max/µ, (76) ≤ (1− αµ)k+1‖xk − x∗‖+ 2 /µ+ Crγ2max/µ, (77)
where the inequality in Eq. 76 follows since ∑k i=0(1− αµ)i ≤ 1 αµ .

heading: B NORM OF THE DIFFERENCE BETWEEN GRADIENTS
text: For ridge regression fi(x) = 12 (〈ai, x〉 − yi) 2 + λ2 ‖x‖ 2, we have ∇fi(x) = ai(〈ai, x〉 − yi) + λx. Therefore,
‖∇fi(x)−∇fj(x)‖ = (‖ai − aj‖.‖x‖+ ‖yi − yj‖)‖aj‖ (78)
For ‖ai‖ ≤ 1, and yi = yj we get
‖∇fi(x)−∇fj(x)‖ ≤ ‖ai − aj‖O(‖x‖) (79)
For reguralized logistic regression with y ∈ {−1, 1}, we have ∇fi(x) = yi/(1 + eyi〈ai,x〉). For yi = yj we get
‖∇fi(x)−∇fj(x)‖ = e‖ai−aj‖.‖x‖ − 1
1 + e−〈ai,x〉 ‖aj‖. (80)
For ‖ai‖ ≤ 1, using Taylor approximation ex ≤ 1 + x, and noting that 11+e−〈ai,x〉 ≤ 1 we get
‖∇fi(x)−∇fj(x)‖ ≤ ‖ai − aj‖.‖x‖ 1 + e−〈ai,x〉 ‖aj‖ ≤ ‖ai − aj‖O(‖x‖). (81)
For classification, we require yi = yj , hence we can select subsets from each class and then merge the results. On the other hand, in ridge regression we also need |yi − yj | 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

Review summaries:
This paper proposes a novel extension to SGD/incremental gradient methods called CRAIG. The algorithm selects a subset of datapoints to approximate the training loss at the beginning of each epoch in order to reduce the total amount of time necessary to solve the empirical risk minimization problem. In particular, the algorithm formulates a submodular optimization problem based on the intuition that the gradient of the problem on the selected subset approximates the gradient of the true training loss up to some tolerance. Each datapoint in the subset is a medoid and assigned a weight corresponding to the number of datapoints in the full set that are assigned to that particular datapoint. A greedy algorithm is employed to approximately solve the subproblem. Theory is proven based on based on an incremental subgradient method with errors. Experiments demonstrate significant savings in time for training both logistic regression and small neural networks.

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.