Title: LEARNING AND DATA SELECTION IN BIG DATASETS

Abstract: Finding a dataset of minimal cardinality to characterize the optimal parameters of a model is of paramount importance in machine learning and distributed optimization over a network. This paper investigates the compressibility of large datasets. More specifically, we propose a framework that jointly learns the input-output mapping as well as the most representative samples of the dataset (sufficient dataset). Our analytical results show that the cardinality of the sufficient dataset increases sublinearly with respect to the original dataset size. Numerical evaluations of real datasets reveal a large compressibility, up to 95%, without a noticeable drop in the learnability performance, measured by the generalization error.

1. INTRODUCTION:

During the last decade, new artificial intelligence methods have offered outstanding prediction performance on complex tasks including face and speech recognition (Zhao et al., 2003; Schalkwyk et al., 2010; Hinton et al., 2012) , autonomous driving (Michels et al., 2005) , and medicine (Kourou et al., 2015) . To achieve such amazing results, state-of-the-art machine learning methods often need to be trained on increasingly large datasets. For example, (MNIST) is a typical dataset for natural image processing of handwritten digits with more than 70,000 samples, and (MovieLens) 20M is a typical dataset for recommendation systems that includes more than 20,000,000 ratings. As we show throughout this paper, most of the samples in these datasets are redundant, carrying almost no additional information for the learning task. A fundamental open question in learning theory is how to characterize and algorithmically identify a small set of critical samples, hereafter called a small representative dataset, that best describes an unknown model. Studying its behavior around those critical samples helps us to better understand the unknown model as well as the inefficiencies of the sample acquisition process in the original dataset. For instance, in a multi-agent system, it may be enough to share these small representative datasets among the agents instead of the original big ones, leading to a significant reduction in power consumption and networking latency (Jiang et al., 2018) .

Experiment design (Sacks et al., 1989) or active learning (Settles, 2012) provides algorithmic approaches to obtain a minimal set of samples to be labeled by an "oracle" (e.g., a human annotator). Active learning is well-motivated in many modern machine learning applications where obtaining a new labeled training sample is expensive. The main components of active learning are a parameterized model, a measure of the model's uncertainty, and an acquisition function that decides based on the model's uncertainty the next sample to be labeled. This approach has several challenges including lack of scalability to high-dimensional data (Settles, 2012) (which has been partially addressed in some recent publications (Gal et al., 2017) ), lack of theoretical guarantees, and lack of formal uncertainty measure. More importantly, the acquisition function is usually greedy in the sense that it sequentially finds the new samples to be labeled one-by-one. Consequently, the resulting sub-sampled dataset may not necessarily be a small representative dataset due to the greedy nature of active learning.

In this paper, we investigate a different approach than experiment design or active learning. Instead of reducing the total labeling cost (like active learning), we focus on the following scenario: how to find the a small representative set of samples with cardinality K in a big dataset of labeled samples with cardinality N ( K)? Using optimization theory, we establish a scalable algorithm with provable theoretical guarantees to jointly find K most representative samples. We also show fundamental relations between K and N to guarantee any arbitrary learning performance. Our framework and algorithmic solution approaches can be very useful tools to better understand compressibility of the existing datasets and to improve distributed learning over a multi-agent systems and Internet-of-Things (see Appendix B). Moreover, further analysis of the representative dataset with respect to the original big dataset of an unknown model helps understand the sources of inefficiency in the current sampling process and how to improve it for other similar models.

The main research question of this work, and in particular our order analysis, is closely related to the sample complexity (Clarkson et al., 2012) , information-theoretic concepts of sampling (Jerri, 1977) like the Shannon-Nyquist sampling, compression (Cover & Thomas, 2012) , and compressive sensing (Donoho, 2006) when the function is sparse in some predetermined basis. All these methods address the following question: how many samples are required to reconstruct a function with a predefined error? We show that the size of such a compressed dataset grows sub-linearly with respect to the cardinality of the original dataset.

In this study, we investigate compressibility of large datasets and develop a general framework for function approximation in which choosing the samples that best describe the function is done jointly with learning the function itself. We formulate a corresponding mixed integer non-linear program, and propose an iterative algorithm that alternates between a data selection step and a function approximation step. We show the convergence of the proposed algorithm to a stationary point of the problem, despite the combinatorial nature of the learning task. We then demonstrate that our algorithm outputs a small dataset of carefully chosen samples that solves a learning task, as accurately as if it were solved using original large dataset. Comprehensive numerical analyses on synthetic and real datasets reveal that our algorithm can significantly compress the datasets, by as much as 95%, with almost no noticeable penalty in the learning performance.

The rest of the paper is organized as follows. Section 2 presents the problem setting and our algorithmic solution approach. Section 3 provides main theoretical results. We apply our algorithms on synthetic and real datasets in Section 4, and then conclude the paper in Section 5. Due to lack of space, we have moved all the proofs and some applications to the appendix.

Notation: Normal font a or A, bold font a, and calligraphic font A denote scalar, vector, and set, respectively. |A| is the cardinality of set A. I is indicator function. a T denotes the transpose of a, and a 0 and a 2 are its l 0 and l 2 norms. 1 is a vector of all ones of proper size. For any integer N , [N ] denotes set {1, 2, . . . , N }.

2. SETTING AND SOLUTION APPROACH:

2.1 PROBLEM SETTING Consider input space X , output space Y, an unknown function f : X → Y from some function space F, an index set [N ] := {1, . . . , N } for N ∈ N + , and a dataset of training samples D = {(x i , f (x i ))} i∈ [N ] , where x i ∈ X . In a classical regression problem, we use the dataset D to learn f , namely find a function h : X → Y that has a minimal distance (for some distance measure, also called loss) to the true function f . Formally, for a given loss function : X × Y × Y → [0, ∞], the regression task solves the following empirical risk minimization problem:

(P 1) : h ∈ arg min h∈F 1 N i∈[N ] (x i , f (x i ), h(x i )) .

(1)

Here, we assume that h ∈ F. However, considering a different function class for h would not change the generality of our results.

In many applications (see Appendix B) one might not want to work with the entire dataset D, e.g., due to its large size, but rather with a small subset E ⊆ D, where possibly |E| |D|. We associate a binary variable z i with each training sample

(x i , f (x i )) such that z i = I{(x i , f (x i )) ∈ E}, representing sample (x i , f (x i )) being selected or dropped. Letting z = [z 1 , • • • , z N ] T ,

the novel problem of jointly learning h and selecting E can be formulated as

EQUATION

EQUATION

where constraint (2b) prevents overfitting when "generalizing" from E to D, and constraint (2c) prevents degenerate/trivial solutions to the problem (e.g., where E empty). We show later that K is a very important parameter that trades off the compression rate, defined as 1 -|E|/|D|, and the generalization error of learning with E.

Let i (h) := (x i , f (x i ), h(x i )) denote the loss corresponding to sample (x i , f (x i )). Followings are some assumptions used throughout the paper, which are prevalent in the learning literature.

Assumption 1. i (h) is continuous and convex in h.

Assumption 2. The original dataset D is duplicate-free, namely x m = x n for all m, n ∈ [N ] and x m , x n ∈ X . Moreover, for all h ∈ F, x m = x n implies m (h) = n (h).

Note that i (h) does not have to be smooth and differentiable in general. We stress the existence of a large family of loss functions, including L p spaces, for which both assumptions hold, as exemplified in Appendix C. Moreover, the convexity in Assumption 1 may also be relaxed at the expense of a weaker convergence property using the block successive upper-bound minimization framework (Razaviyayn et al., 2013) . Finally, if the dataset contains some duplicated samples, we can add insignificant perturbations to satisfy Assumption 2, without affecting the information content (Bertsekas, 1998) .

2.2. SOLUTION APPROACH:

(P 2) is a non-convex combinatorial optimization problem with coupling cost and constraint functions.

In the following, we provide a solution approach based on block-coordinate descent (BCD), which splits (P 2) into two subproblems: (P 2a) for data selection and (P 2b) for function approximation. Let h (k) and z (k) be the value of h and z at iteration k. BCD yields the following update rules:

EQUATION

(P 2b) : h (k+1) ∈ arg min h∈F g(h, z (k+1) ), s.t. g 1 (h) ≤ .

(4) The data selection (D-)step (P 2a) is optimized for a given hypothesis h. Then, in the function approximation (F-)step (P 2b), the hypothesis is optimized for the updated compressed dataset z (k+1) . Next, we derive solutions to each step.

2.2.1. D-STEP:

Given h (k) , namely the value of h at iteration k, we let c

EQUATION

N ], for notation simplicity. Then, the first subproblem is written as

EQUATION

where we have used the fact that z T 1 ≥ K holds with equality; shown in Proposition 1. Thus, 1 T z = K can be removed from the denominator to equivalently1 write (5) as

EQUATION

)

Though combinatorial, the form in (P 2a) allows to easily derive its solution, as shown by Proposition 1, i.e., z k) with the smallest values. In other words, the optimal solution is obtained by selecting the smallest K elements of c (k) , and setting the corresponding indices of z to 1.

(k+1) i = I{i ∈ S (k) }, where S (k) is the set of K-indices in c (

2.2.2. F-STEP:

Given the updated data selection, z (k+1) , we use the fact that 1 T z (k+1) = K and rewrite (P 2b) as

EQUATION

)

Algorithm 1 Alternating Data Selection and Function Approximation (DF) k+1) by solving (P 2a) in (5), using Proposition 1 // F-step Update h (k+1) by solving (P 2b) in ( 7)

// Initialize z (1) = 1 for k = 1, 2, 3, . . . do // D-step Compute c (k) i , ∀i ∈ [N ] Update z (

Stop if |g(h (k+1) , z (k+1) ) -g(h (k) , z (k) )| < γ It follows from z (k+1) i

∈ {0, 1} that both the cost and constraints in (P 2b) consist of convex combinations of the loss, i , assumed convex (Assumption 1): Thus, (P 2b) is convex in h (Boyd & Vandenberghe, 2004, Chap. 3 ). In the following section, we show that there exist loss functions that lead to closed-form solutions.

Algorithm 1 summarizes our alternating data selection and function approximation steps. We establish the algorithm's convergence to a stationary point of (P 2) in Proposition 2.

2.3. SPECIAL CASE 1: LINEAR REGRESSION AND DATA SELECTION:

We specialize our approach to a linear regression problem with data selection where h(

x i ) = x T i w for w, x i ∈ R d , d being the dimension of each sample. Optimization problem (P 2) reduces to arg min w,z i∈[N ] z i 1 T z x T i w -f (x i ) 2 s.t. 1 N i∈[N ] x T i w -f (x i ) 2 ≤ , 1 T z ≥ K , z ∈ {0, 1} N .

The D-step is identical to (P 2a), which can be solved by Proposition 1. Given z (k+1) , the F-step reduces to the following quadratically-constrained quadratic programming:

w (k+1) = arg min w A (k+1) X T w -f (X) 2 2 , s.t. X T w -f (X) 2 2 ≤ ,

where

A (k+1) := diag z (k+1) 1 , • • • , z (k+1) N , X := [x 1 , • • • , x N ], and f (X) := [f (x 1 ), . . . , f (x N )] T .

The problem is convex and can be solved using standard Lagrangian techniques to yield

w (k+1) = X A (k+1) 2 + λI N X T -1 X A (k+1) + λI N f (X) ,

where λ ≥ 0 is a Lagrange multiplier that satisfies the complementary slackness condition and can be found using 1D search methods.

Note that computational complexity of w (k+1) is dominated by the matrix inversion, which is in the order of O(d 3 ) and does not scale with the size of the training set, N . However, if it is still considered significant in some applications, w (k+1) can be obtained using stochastic gradient-based methods whose computational complexity is O(d/α) for accuracy threshold α > 0.

2.4. SPECIAL CASE 2: ROBUST LEARNING:

Consider the following continuous relaxation of (P 2), where z ∈ {0,

1} N is relaxed to z ∈ [0 , 1] N , (P 3) : arg min h∈F ,z 1 1 T z i∈[N ] z i i (h) s.t. 1 N i∈[N ] i (h) ≤ , 1 T z ≥ K , z ∈ [0, 1] N .

Thus, z i can be seen as a non-negative weight assigned to sample (x i , f (x i )) in that training set, representing level of confidence in its quality (higher values of z i imply better qualities). From this perspective, the resulting problem becomes a robust learning problem, in presence of non-uniform sample quality: learning h, jointly with the best samples of the training set. Some applications include de-noising sensor measurements and outlier detection.

We can use Algorithm 1 to address (P 3), with a minor modification in the D-step. Note that relaxing the binary constraint implies that the D-step cannot be solved using the simple method of Proposition 1. However, we use the linear program relaxation of (P 2a),

EQUATION

which can be efficiently solved using standard linear program solvers. In Appendix A.6, we have shown that optimization problem ( 8) is equivalent to (P 2a), and therefore the linear program relaxation is optimal.

3. MAIN THEORETICAL RESULTS:

In this section, we present main results of this paper. Detailed proofs are available in Appendix.

Proposition 1 (Solution of (P 2a)). Under Assumption 2, we define an index sequence j such that for any

j m , j n ∈ [N ], c (k) jm < c (k) jn iff m < n, where c (k) i

is defined in Section 2.2.1. The solution of the data selection subproblem is

EQUATION

and z (k+1) 0 = K.

Proposition 1 implies that the optimal solution to the binary data selection subproblem is simple: evaluate the loss function for all training samples of D using h (k) , sort the values, and keep K data samples having the smallest losses. Next, we establish the convergence of our BCD-based algorithm to a stationary point of (P 2).

Proposition 2 (Convergence). Let {g(h (k) , z (k) )} k≥1 denote the sequence generated by the BCD updates in Algorithm 1. Then, this sequence monotonically decreases with each update, i.e.,

g(h (k) , z (k) ) ≥ g(h (k) , z (k+1) ) ≥ g(h (k+1) , z (k+1) ), k = 1, 2, 3, . . .

and converges to a stationary point of (P 2).

Proposition 3 (Computational Complexity). The D-step is run in O(N ), and the F-step has the same complexity as of (P 1), per iteration of Algorithm 1.

To analyze the asymptotic behavior of our approach, we make additional assumptions on the class of loss functions, and on F. In particular, we assume that belongs to the L p space and F is the space of L-Lipschitz functions defined on some compact support.

Proposition 4 (Sample Complexity). Assume that the samples are noiseless, and F is the set of L-Lipschitz functions defined on interval [0, T ] d . Consider optimization problem (P 2). Let

g (h , z ) := 1 T z -1 i∈[N ] z i (x i , f (x i ), h (x i )).

For any arbitrary constant δ > 0, the following holds:

When (x m , f (x m ), h (x m )) := |f (x m ) -h (x m )| 2 , g (h , z ) ≤ δ if K ≥ (1 + 2LT d/δ) d

, where • is the ceiling function.

Corollary 1 (Asymptotic Sample Complexity). Consider the assumptions of Proposition 4. Define compression ratio as CR := 1 -K/N . As N grows large: ∀δ > 0, ∃K ≤ N s.t. g (h , z ) < δ , and CR → 1 .

Proposition 4 and Corollary 1 imply that the sufficient dataset E (which can be used to learn any function f in class F with any arbitrary accuracy) is the output of a sub-linear sub-sampler (our Algorithm 1) of the original big dataset, namely K/N → 0 asymptotically. In other words, as we experimentally show in the next section, most of the existing big datasets are highly redundant, and the redundancy brings almost no additional gain for the accuracy of the learning task. Finding this sufficient dataset is of manageable complexity, as specified in Proposition 3.

4. EXPERIMENTAL RESULTS:

In this section, we first present two toy examples to illustrate (P 2) and algorithmic solution. We then focus on real databases and evaluate the effectiveness of our approach on finding the small representative dataset with a negligible loss in the generalization capability.

4.1. EXPERIMENTAL SETTING:

In every experiment, we select the input space X , output space Y, mapping f , training dataset D, test dataset T , and hypothesis class H. We then run Algorithm 1 to find the optimal compressed dataset E ⊆ D. We run this experiment for different values of CR.

To evaluate the true generalization capability of E , we run a conventional regression problem (P 1) using both D and E , find the corresponding optimal approximation function, and then evaluate their accuracy on T . We denote the normalized generalization error by e(D, T ),

EQUATION

)

where h is found by running (P 1) with D. When we find h by running (P 1) with E , (10) shows e(E , T ). When is the L 2 -norm, (10) reduces to the normalized square error, i∈T |f

(x i ) - h (x i )| 2 / i∈T |f (x i )| 2 .

This normalization is to have a fair comparison of the generalization error over various test datasets, which may have different norms. We say that our compressed dataset E performs as well as D when e(E , T ) is close to e(D, T ). Throughout the following experimental studies, we observe that the lower the compression ratio, the lower the gap |e(D, T ) -e(E , T )|.

However, for big datasets, we can substantially compress the dataset without any noticeable drop in the gap, indicating a high inefficiency in the data generation process.

4.2. ILLUSTRATIVE EXAMPLES:

In our first example, we pick a very smooth function f . Let X = [0, 8], F = H = Poly(10), where Poly(n) is polynomial functions of degree n. f (x) is given in figure 1 (a). Our original dataset D is 100 equidistant samples in X . We add i.i.d. Gaussian noise of standard deviation 0.05 to f (x). Figure 1 (a) shows an example of an optimal compressed dataset E computed for K = 12, along with the learned hypothesis h , which almost perfectly coincide with the true function f . Note that due to the random noise in D, the selected samples would be different in every run. To evaluate the true generalization capability of our algorithm, we run a Monte Carlo simulation for 100 random realizations of the noise on dataset, find E and h for each realization, and compute e(D, T ) and e(E , T ) from ( 10). Note that T is the set of 1000 equidistant examples in X . 

4.3. REAL DATASETS:

Motivated by the excellent performance of the proposed algorithm on simple syntectic data, in this section, we apply Algorithm 1 on real databases listed in Table 1 , available on Statlib (Sta) and UCI repositories (UCI). These databases have been extensively used in relevant machine learning and signal processing applications (Schapire, 1999; Aharon et al., 2006; Zhang & Li, 2010; Zhou et al., 2014; Tang et al., 2016; Chatterjee et al., 2017) .

For the learning task, and without loss of generality, we use the recently proposed extreme learning machine (ELM) architecture (Huang et al., 2006; 2012) , due to its implementation efficiency, good regression and classification performance, and convexity of the resulting optimization problem. An ELM typically uses a few hidden layers, each having many nodes, to project the input data vectors to high dimensional feature vectors. Then, a linear projection is used at the last layer to recover the output vector. An interesting property of ELM is the ability to use instances of random matrices in mapping to feature vectors (as opposed to usual deep neural networks in which these weight matrices should be optimized), and therefore we need to optimize only the weights of the last layer. In our implementation, we have used a single hidden layer, an instance of random matrix between input layer and hidden layer R, an element-wise rectifier linear unit (ReLU) function at each hidden node that returns σ(•) = max(•, 0), weights to the output layer W . Given dataset D with N samples and a binary selection vector z (k+1) , we use the following optimization problem in the F-step of Algorithm 1 at iterate k: 0% 0.0245 ± 0.0051 0.0301 ± 0.0056 0.1323 ± 0.0134 0.0082 ± 0.0007 0.0142 ± 0.0008 25% 0.0294 ± 0.0058 0.0345 ± 0.0071 0.1325 ± 0.0142 0.0083 ± 0.0007 0.0144 ± 0.0008 50% 0.0333 ± 0.0067 0.0323 ± 0.0080 0.1338 ± 0.0175 0.0085 ± 0.0008 0.0144 ± 0.0009 75% 0.0360 ± 0.0060 0.0374 ± 0.0076 0.1351 ± 0.0169 0.0086 ± 0.0010 0.0145 ± 0.0008 80% 0.0417 ± 0.0077 0.0382 ± 0.0076 0.1354 ± 0.0171 0.0086 ± 0.0009 0.0145 ± 0.0008 95% 0.0630 ± 0.0134 0.0538 ± 0.0122 0.1386 ± 0.0217 0.0088 ± 0.0012 0.0153 ± 0.0011

W (k+1) = arg min W 1 1 T z (k+1) i∈[N ] z (k+1) i f (x i ) -W σ(Rx i ) 2 2 + λ W 2 2 s.t. 1 N i∈[N ] f (x i ) -W σ(Rx i ) 2 2 ≤ ,

where the second term in the objective function is the Tikonov regularization with parameter λ that alleviates the over-fitting problem. This convex quadratically constrained quadratic program can be efficiently solved by existing optimization toolboxes (Grant & Boyd, 2014) . Due to the randomness in R, we repeat experiments 100 times and report the mean value and standard deviation of the performance results.

Table 2 shows the regression performance of various databases, given in Table 1 . From this table, our approach can substantially compress the training datasets with a controlled loss on the generalization error. This compressibility increases by the dataset size and decreases by the input dimension.

For instance, 75% compression in Bodyfat results in a noticeable performance drop, while a big dataset like YearPredictionMSD can be compressed by 95% without a significant loss in the learning performance. Moreover, these results empirically show the sub-linear characteristic of the sufficient dataset for any fixed δ, namely CR → 1 as N → ∞. The results suggest that a small set of carefully selected samples are enough to run the learning task. As we have discussed in Appendix B, a similar preprocessing to identify and send only a small representative dataset could be inevitable to implement distributed learning over communication-constrained networks, e.g., intra-body networks or Internet-of-Things.

Further Discussions: In all the simulation experiments, we have observed a very fast convergence of the proposed Algorithm 1, usually after a few iterations in small datasets (e.g., Bodyfat) and a few tens of iterations in large datasets (e.g., Abalone) among D-step and F-step. Computational complexity of each step is characterized in Proposition 3. Both bigger datasets and larger N K correspond to larger search spaces and consequently slower convergence rate.

5. CONCLUSIONS:

We addressed the compressibility of large datasets, namely the problem of finding dataset of minimal cardinality (small representative dataset) for a learning task. We developed a framework that jointly learns the input-output mapping and the representative dataset. We showed that its cardinality increases sub-linearly with respect to that of the original dataset. While an asymptotic compressibility of almost 100% is available in theory, we have observed that real datasets may be compressed as much as 95% without any noticeable drop of the learning performance. These results challenge the efficiency and benefits of the existing approaches to create big datasets and serve as benchmark for distributed learning over communication-limited networks. Consider loss function

x = |f (x) -h(x)| 2 . When d = 1

, it is easy to see from Lemma 3 and figure A.1 that max si ≤ 4L 2 (x i -x i-1 ) 2 for every set S i , where x i -x i-1 is the measure of set S i . Now, since sets {S i }, i ∈ [K -1] have the same measure (defined based on equidistant grid points), we have

x i -x i-1 = T /(K -1), so E x x ≤ max i max si ≤ 4L 2 T 2 (K -1) 2 . By setting g (h , z ) ≤ E x x ≤ δ, we get K ≥ 1 + 2LT / √ δ.

For d > 1, we can define equidistant marks on every coordinate of X and define a grid of (K 1/d -1) d disjoint sets {S i } i , where we have assumed that K 1/d is an integer number to avoid unnecessary notation complications. The distance between two consecutive marks on every coordinate is T /(K 1/d -1), and therefore from Lemma 3

E x ≤ max i max si ≤ 2L T √ d K 1/d -1 2 . By setting g (h , z ) ≤ E x x ≤ δ, we get K ≥ 1 + 2LT d/δ d . This completes the proof. A.4 LEMMA 1 Assume i∈[N ] z (k+1) i = i∈[N ] z (k+1) ji = M ≥ K. For k = 1, if z (k+1) j1 = 1 the statement holds. If z (k+1) j1

= 0, then take any n for which z (k+1) jn = 1 and observe that the following inequality holds by definition of index set j:

i∈[N ] z (k+1) ji c (k) ji M = i∈[N ]\{n} z (k+1) ji c (k) ji + c (k) jn M ≥ i∈[N ]\{n} z (k+1) ji c (k) ji + c (k) j1 M , since c (k) j1 < c (k)

jn for any n. This completes the proof for k

= 1. For k = 2 ≤ K, if z (k+1) j2 = 1 the statement holds. If z (k+1) j2

= 0, then take any n ≥ 3 for which z

(k+1) jn = 1. Use z (k+1) j1 = 1 and observe that i∈[N ] z (k+1) ji c (k) ji M = c (k) j1 + c (k) jn + i∈[N ]\{1,n} z (k+1) ji c (k) ji M ≥ c (k) j1 + c (k) j2 + i∈[N ]\{1,n} z (k+1) ji c (k) ji M since c (k) j2 < c (k)

jn for any n > 2. We can use the same arguments recursively to prove that z

(k+1) jm = 1 for any m ≤ K. A.5 LEMMA 2 Assume i∈[N ] z (k+1) ji = M ≥ K. By Lemma 1, z (k+1) ji = 1 for all i ≤ K. We should show that i∈[K] c (k) ji K ≤ i∈[K] c (k) ji + N i=K+1 z (k+1) ji c (k) ji M for any z (k+1) ji

. This is clearly true as the left-hand-side is the average of the K smallest values of the loss function on dataset of size N . In particular,

i∈[K] c (k) ji + N i=K+1 z (k+1) ji c (k) ji M ≥ i∈[K] c (k) ji + M -K c (k) j K + c (k) j K + . . . + c (k) j K M = i∈[K] c (k) ji K + ≥0 (M -K) Kc (k) j K -(M -K) i∈[K] c (k) ji M K (a) ≥ i∈[K] c (k) ji K ,

where (a) holds as Kc (k)

j K ≥ i∈[K] c (k)

ji . This completes the proof.

A.6 OPTIMALITY OF (8):

To prove the optimality of (8), recall that 1 T z = K in (P 2a) is of the form Az = B, where A is a totally unimodular matrix, and B is an integer. Thus, optimization problem ( 8) is equivalent to (P 2a), and the linear program relaxation is optimal.

B APPLICATIONS TO MULTI-AGENT SYSTEMS AND INTERNET-OF-THINGS:

Consider a network of agents with some abstract form of very limited networking capacity (the so-called communication-limited networks in optimization and control literature (Smith et al., 2018; Magnússon et al., 2018; Nedic & Bertsekas, 2001; Tsitsiklis & Luo, 1987) ). The limitation can be, among others, due to low-power operation of agents (sensor nodes) in Internet-of-Things or low channel capacity in harsh wireless environments, like intra-body networks. Consequently, we have a very limited transmission rate among various agents. Each agent has a local dataset, e.g., some measurements, and they should all share the datasets to jointly run an optimization/learning task within a short deadline. Due to the tight deadline and limited communication capability, we cannot send all entries of the datasets to all other agents by the deadline. An important question here is to decide, locally at every agent, which data should be shared (data selection). To address this problem, one may consider the existence of an oracle that has access to all the datasets, finds the minimal representative dataset, and then informs all the agents what to share. Clearly, this oracle is not practical, rather it gives a theoretical benchmark on the performance of various solution approaches and on the cardinality of the minimal representative dataset. Our (P 2) models the oracle, and our algorithmic solution approaches characterize the solution of the oracle. We can also get useful insights on the "optimal" sub-sampling per agent, which can be exploited to develop practical solution approaches for the original problem. Therefore, our results are of paramount importance for the interesting problem of low-latency learning and inference over a communication-limited network.

C ADDITIONAL EXAMPLES:

The following example shows the generality of Assumptions 1 and 2.

Example 4. Let P denote the space of polynomial functions on R, f (x) = e x (∈ P), h(x) = N -1 n=0 x n /n! (∈ P) be the first N (< ∞) terms of the Taylor expansion of f (x), and n (h) := |f (x n ) -h(x n )| 2 . n (h) is compatible with Assumption 1. Moreover, for almost any x n , x m ∈ R

:

Two problems are equivalent, if the optimal solution to one can be obtained from the other, and viceversa(Boyd & Vandenberghe

, 2004).

Figures and Tables:

Figure 1(b) reports the average generalization error against the CR. As expected, the higher the compression ratio the higher the generalization error. The tail drop of this function is at least as a fast as a double exponential function of CR, implying that for a wide range of CR values, the extra generalization error |e(E , T ) -e(D, T )| is negligible. In particular, in the example of Figure 1(b) with N = 100, 70% compression leads to only 6 × 10 -5 extra generalization error.

Figure2illustrates the performance of our optimization problem on a less smooth function than that of figure 1(a). In particular, we consider the function of figure 2(a) which is from F = Poly(15), X = [0, 1], N = 100 equidistant samples from X in D, and 1000 equidistant samples from X in T . We observe a large compression of the dataset in figure2(b), where only |e(E , T ) -e(D, T )| = 2.3 × 10 -5 extra true generalization error after 60% compression. Moreover, we can see the double exponential tail of e(E , T ), which implies that only a small dataset of a few carefully chosen samples are enough to have a learning with a sufficiently good generalization capability. However, for a fixed error gap |e(E , T ) -e(D, T )|, functions having more variation are less compressible, since more samples are needed to maintain the same error gap.

Figure1: Learning compressed data set E and optimal hypothesis h with a dataset of size N = 100. Function f ∈ Poly(10) + noise. Selected samples in E (⊂ D) are denoted by circles. In the example of (a), f and h are visually indistinguishable. The higher the compression ratios the higher the generalization error. e(E , T ) behaves like a double exponential function of CR. For a wide range of CR values there is almost no loss compared to e(D, T ).

Figure A.1: Illustration of the functional class F. Input space X is divided into disjoint sets {Si}. F is the set of all L-Lipschitz functions passing through samples/marks {xi} i∈[K] . All functions f, h ∈ F lie in the dashed red parallelograms. The slopes of these parallelograms are ±L. Three possible functions are shown in the figure.

Databases for regression task.

Table content: <table><tr><td>Database</td><td colspan="3"># Training samples # Test samples Input dimension</td></tr><tr><td>Bodyfat</td><td>168</td><td>84</td><td>14</td></tr><tr><td>Housing</td><td>337</td><td>169</td><td>13</td></tr><tr><td>Space-ga</td><td>2,071</td><td>1,036</td><td>6</td></tr><tr><td>YearPredictionMSD</td><td>463,715</td><td>51,630</td><td>90</td></tr><tr><td>Power Consumption</td><td>1,556,445</td><td>518,814</td><td>9</td></tr></table>

Regression performance on real databases. The reported values are "average ± standard deviation" of the normalized true generalization error. "CR" stands for compression ratio. The values on row CR = 0% corresponds to e(D, T ).

Table content: <table><tr><td>CR</td><td>Bodyfat</td><td>Housing</td><td>Space-ga</td><td>YearPredictionMSD Power Consumption</td></tr></table>

Peer review (iclr2019):
 This paper proposes an algorithm for data compression which learns the most representative samples of a dataset by solving a mixed integer program. The paper then finds a stationary point of the program by a greedy method. Empirical results show that compressing a large fraction of samples using the method does not cause a noticeable drop in learning performance. The algorithm and theory in this paper are clearly presented. The idea however is not quite novel. The impact of reducing samples on learning problems have been studied in literature, e.g. [1][2][3]. Algorithms like sketching [2] and importance sampling [3] have been shown both in theory and practice that can reduce dataset sample size without a big drop in performance. Besides, both [2][3] requires less computation time than solving the original problem. However, the proposed algorithm takes more computation than minimizing the loss. Some applications of the proposed algorithm are presented in appendix B. However, it would be better to also show the learning or optimization problems in these applications since they can help us understand how to apply the proposed algorithm in practice. Some experiments on these applications are also preferred. For experiments, both synthetic data and large-scale real data are tested. However, it would be better comparing with other algorithms like [1][2][3] to show the novelty of the proposed algorithms. In section 4.2, the synthetic data experiments are conducted using only one constructed dataset and multiple random noise. Generating a random dataset each time will show that the performance of the algorithm is not a special case for some specific datasets and will also give insights on the robustness of the algorithm. The proposed method can be applied to both regression and classification problems. Some experiments on classification would help better exploring this algorithm. The experiments on reduced sample size and learning performance are comprehensive. However, I would like to see experiments on running time, since running time represents the key performance metric in real world applications. Some questions: In Section 2.4: Please explain what robust learning in your case is. In your experiments, how you choose your parameter epsilon? What is the average loss in your experiments? Average loss is more commonly used in practice. In your experiments, what if you pick these selected samples and train them? Data compression is also well studied in learning theory [4] (Part IV, Section 30). It is also interesting to know if the proposed algorithm can reach the theoretical bound. [1] Understanding Black-box Predictions via Influence Functions [2] sketching as a tool for numerical linear algebra [3] Randomized algorithms for matrices and data [4] understanding machine learning from theory to algorithms