On Coresets for Logistic Regression
Alexander Munteanu, Chris Schwiegelshohn, Christian Sohler, David P. Woodruff
Introduction
Scalability is one of the central challenges of modern data analysis and machine learning. Algorithms with polynomial running time might be regarded as efficient in a conventional sense, but nevertheless become intractable when facing massive data sets. As a result, performing data reduction techniques in a preprocessing step to speed up a subsequent optimization problem has received considerable attention. A natural approach is to sub-sample the data according to a certain probability distribution. This approach has been successfully applied to a variety of problems including clustering (Langberg & Schulman 2010; Feldman & Langberg 2011; Barger & Feldman 2016; Bachem et al. 2018), mixture models (Feldman et al. 2011; Lucic et al. 2016), low rank approximation (Cohen et al. 2017), spectral approximation (Alaoui & Mahoney 2015; Li et al. 2013), and Nyström methods (Alaoui & Mahoney 2015; Musco & Musco 2017).
The unifying feature of these works is that the probability distribution is based on the sensitivity score of each point. Informally, the sensitivity of a point corresponds to the importance of the point with respect to the objective function we wish to minimize. If the total sensitivity, i.e., the sum of all sensitivity scores , is bounded by a reasonably small value , there exists a collection of input points known as a coreset with very strong aggregation properties. Given any candidate solution (e.g., a set of centers for -means, or a hyperplane for linear regression), the objective function computed on the coreset evaluates to the objective function of the original data up to a small multiplicative error. See Sections 2 and 4 for formal definitions of sensitivity and coresets.
Our first contribution is an impossibility result: logistic regression has no sublinear streaming algorithm. Due to a standard reduction between coresets and streaming algorithms, this also implies that logistic regression admits no coresets or bounded sensitivity scores in general.
Our third contribution is an analysis of our sampling distribution for a parametrized class of instances we call -complex, placing our work in the framework of beyond worst-case analysis (Balcan et al. 2015; Roughgarden 2017). The parameter roughly corresponds to the ratio between the log of correctly estimated odds and the log of incorrectly estimated odds. The condition of small is justified by the fact that for instances with large , logistic regression exhibits methodological problems like imbalance and separability, cf. (Mehta & Patel 1995; Heinze & Schemper 2002). We show that the total sensitivity of logistic regression can be bounded in terms of , and that our sampling scheme produces the first coreset of provably sublinear size, provided that is small.
All proofs and additional plots from the experiments are in the appendices A and B, respectively.
Preliminaries and Problem Setting
In this paper we assume we have a very large number of observations in a moderate number of dimensions, that is, . In order to speed up the computation and to lower memory and storage requirements we would like to significantly reduce the number of observations without losing much information in the original data. A suitable data compression reduces the size to a sublinear number of data points while the dependence on and the approximation parameters may be polynomials of low degree. To achieve this, we design a so-called coreset construction for the objective function. A coreset is a possibly (re)weighted and significantly smaller subset of the data that approximates the objective value for any possible query points. More formally, we define coresets for the weighted logistic regression function.
-Complex Data Sets We will see in Section 3 that in general, there is no sublinear one-pass streaming algorithm approximating the objective function up to any finite constant factor. More specifically there exists no sublinear summary or coreset construction that works for all data sets. For the sake of developing coreset constructions that work reasonably well, as well as conducting a formal analysis beyond worst-case instances, we introduce a measure that quantifies the complexity of compressing a given data set.
weighted by is called -complex if .
We conjecture that computing the value of is hard. However, it can be approximated in polynomial time. It is not necessary to do so in practical applications, but we include this result for those who wish to evaluate whether their data has nice -complexity.
Lower Bounds
At first glance, one might think of taking a uniform sample as a coreset. We demonstrate and discuss on worst-case instances in Appendix C that this won’t work in theory or in practice. In the following we will show a much stronger result, namely that no efficient streaming algorithms or coresets for logistic regression can exist in general, even if we assume that the points lie in -dimensional Euclidean space. To this end we will reduce from the INDEX communication game. In its basic variant, there exist two players Alice and Bob. Alice is given a binary bit string and Bob is given an index . The goal is to determine the value of with constant probability while using as little communication as possible. Clearly, the difficulty of the problem is inherently one-way; otherwise Bob could simply send his index to Alice. If the entire communication consists of only a single message sent by Alice to Bob, the message must contain bits (Kremer et al. 1999).
A similar reduction also holds if Alice’s message consists of points forming a coreset. Hence, the following corollary holds.
We note that the proof can be slightly modified to rule out any finite additive error as well. This indicates that the notion of lightweight coresets with multiplicative and additive error (Bachem et al. 2018) is not a sufficient relaxation. Independently of our work Tolochinsky & Feldman 2018 gave a linear lower bound in a more general context based on a worst case instance to the sensitivity approach due to Huggins et al. 2016. Our lower bounds and theirs are incomparable; they show that if a coreset can only consist of input points it comprises the entire data set in the worst-case. We show that no coreset with can exist, irrespective of whether input points are used. While the distinction may seem minor, a number of coreset constructions in literature necessitate the use of non-input points, see (Agarwal et al. 2004) and (Feldman et al. 2013).
Sampling via Sensitivity Scores
The sensitivity of a point measures its worst-case importance for approximating the objective function on the entire input data set. Performing importance sampling proportional to the sensitivities of the input points thus yields a good approximation. Computing the sensitivities is often intractable and involves solving the original optimization problem to near-optimality, which is the problem we want to solve in the first place, as pointed out in (Braverman et al. 2016). To get around this, it was shown that any upper bound on the sensitivities also has provable guarantees. However, the number of samples needed depends on the total sensitivity, that is, the sum of their estimates , so we need to carefully control this quantity. Another complexity measure that plays a crucial role in the sampling complexity is the VC dimension of the range space induced by the set of functions under study.
Recently a framework combining the sensitivity scores with a theory on the VC dimension of range spaces was developed in (Braverman et al. 2016). For technical reasons we use a slightly modified version.
where each element of is sampled i.i.d. with probability from , denotes the weight of a function that corresponds to , and where is an upper bound on the VC dimension of the range space induced by that can be obtained by defining to be the set of functions where each function is scaled by .
Now we show that the VC dimension of the range space induced by the set of functions studied in logistic regression can be related to the VC dimension of the set of linear classifiers. We first start with a fixed common weight and generalize the result to a more general finite set of distinct weights.
We will see later how to bound the number of distinct weights by a logarithmic term in the range of the involved weights. It remains for us to derive tight and efficiently computable upper bounds on the sensitivities.
In the second case, the element under study is bounded by a constant. We consider two sub cases. If there are a lot of contributions, which are not too small, and thus cost at least a constant each, then we can lower bound the total cost by a constant times their total weight. If on the other hand there are many very small negative values, then this implies again that the cost is within a fraction of the total weight.
Combining both lemmas yields general upper bounds on the sensitivities that we can use as an importance sampling distribution. We also derive an upper bound on the total sensitivity that will be used to bound the sampling complexity.
We combine the above results into the following theorem.
holds, where .
Using this, we can show that the -complexity is not violated too much after one stage of sampling.
Let be a sampling and reweighting matrix according to Theorem 15 where parameter is replaced by . That is is the resulting reweighted sample when Theorem 15 succeeds on -complex input . Suppose that simultaneously Lemma 17 holds. Let
Then we have
Now we are ready to prove our theorem regarding the recursive subsampling algorithm.
Experiments
We ran a series of experiments to illustrate the performance of our coreset method. All experiments were run on a Linux machine using an Intel i7-6700, 4 core CPU at 3.4 GHz, and 32GB of RAM. We implemented our algorithms in Python. Now, we compare our basic algorithm to simple uniform sampling and to sampling proportional to the sensitivity upper bounds given by Huggins et al. 2016.
The exact QR-decomposition is rather slow on large data matrices. We thus optimized the running time of our approach in the following way. We used a fast approximation algorithm based on the sketching techniques of Clarkson & Woodruff 2013, cf. (Woodruff 2014). That leads to a provable constant approximation of the square root of the leverage scores with constant probability, cf. (Drineas et al. 2012), which means that the total sensitivity bounds given in our theory will grow by only a small constant factor. A detailed description of the algorithm is in the proof of Theorem 15.
The subsequent optimization was done for all approaches with the standard gradient based optimizer from the scipy.optimize http://www.scipy.org/ package.
Data Sets We briefly introduce the data sets that we used. The Webb Spam https://www.cc.gatech.edu/projects/doi/WebbSpamCorpus.html data consists of unigrams with features from web pages which have to be classified as spam or normal pages ( positive). The Covertype https://archive.ics.uci.edu/ml/datasets/covertype data consists of cartographic observations of different forests with features. The task is to predict the type of trees at each location ( positive). The KDD Cup ’99 http://kdd.ics.uci.edu/databases/kddcup99/kddcup99.html data comprises network connections with features and the task is to detect network intrusions ( positive).
For each data set, we ran all three subsampling algorithms for a number of thirty regular subsampling steps in the range . For each step, we present the mean relative error as well as the trade-off between mean relative error and running time, taken over twenty independent repetitions, in Figure 1. Relative running times, standard deviations and absolute values are presented in Figure 2 respectively in Table 1 in Appendix B.
Evaluation The accuracy of the QR-sampling distribution outperforms uniform sampling and the distribution derived from -means on all instances. This is especially true for small sampling sizes. Here, the relative error especially for uniform sampling tends to deteriorate. While -means sampling occasionally improved over uniform sampling for small sample sizes, the behavior of both distributions was similar for larger sampling sizes. The standard deviations had a similarly low magnitude as the mean values, where the QR method usually showed the lowest values.
The trade-off between the running time and relative errors shows a common picture for Webb Spam and Covertype. QR is nearly always more accurate than the other algorithms for a similar time budget, except for regions where the relative error is large, say above 5-10% while for larger time budgets, QR is better by a factor between - and drops more quickly towards . The conclusion so far could be that for a quick guess, say a -approximation, the competitors are faster, but to provably obtain a reasonably small relative error below 5%, QR outperforms its competitors. However, for KDD Cup ’99, QR always has a lower error than its competitors. Their relative errors remain above 15% or much worse, while QR never exceeds 22% and drops quickly below 4%. As a side note, our estimates for support our experimental findings, especially that KDD Cup ’99 seems more difficult to approximate than the others. The estimated values were for Webb Spam, for Covertype, and for KDD Cup ’99.
The relative running time for the QR-distribution was comparable to -means and only slightly higher than uniform sampling. However, it never exceeded a factor of two compared to its competitors and remained negligible compared to the full optimization task, see Figure 2 in Appendix B. The standard deviations were negligible except for the -means algorithm and the KDD Cup ’99 data set, where the uniform and -means based algorithms showed larger values. The QR method had much lower standard deviations. This indicates that the resulting coresets are more stable for the subsequent numerical optimization.
We note that the savings of all presented data reduction methods become even more significant when performing more time consuming data analysis tasks like MCMC sampling in a Bayesian setting, see e.g., (Huggins et al. 2016; Geppert et al. 2017).
Conclusions
Our experimental evaluation shows that our implementation of the basic algorithm outperforms uniform sampling as well as state of the art methods in the area of coresets for logistic regression while being competitive to both regarding its running time.
Acknowledgments
We thank the anonymous reviewers for their valuable comments. We also thank our student assistant Moritz Paweletz for implementing and conducting the experiments. This work was partly supported by the German Science Foundation (DFG) Collaborative Research Center SFB 876 "Providing Information by Resource-Constrained Analysis", projects A2 and C4 and by the ERC Advanced Grant 788893 AMDROMA.
References
Appendix A Proofs
Assume we had a streaming algorithm using space. We construct the following protocol for INDEX: Consider an instance of INDEX, i.e., Alice has a string and Bob has an index . We transform the instance into an instance for logistic regression. For each , Alice adds a point . Note that all of these points have unit Euclidean norm and hence any single point may be linearly separated from the others. All of Alice’s points have label . Alice summarizes the point set by running the streaming algorithm and sends a message containing the working memory of the streaming algorithm to Bob. Bob now adds the point for small enough with label . From the contents of Alice’s message and , Bob now obtains a solution to the logistic regression instance. Clearly, if Alice added and hence then the optimal solution will have cost at least , since there will be at least one misclassification. If, on the other hand, Alice did not add and hence , then the two point sets are linearly separable and the cost tends to . Distinguishing between these two cases, i.e. approximating the cost of logistic regression beyond a factor solves the INDEX problem.
To conclude the theorem, let us consider the space required to encode the points added by Alice. For the reduction to work, it is only important that any point added by Alice can be linearly separated from the others. This can be achieved by using bits per point, i.e., the space of Alice’s point set is at most . The space bound now follows from the lower bound of bits due to Kremer et al. 1999 for the INDEX problem. ∎
If we had a coreset construction with points, we have a protocol for INDEX: Alice computes a coreset for her point set defined in the proof of Theorem 4 and sends it to Bob. Bob computes an optimal solution on the union of the coreset and his point. This solves INDEX using communication, which contradicts the lower bound of Kremer et al. 1999. So Alice’s coreset cannot exist. ∎
(cf. Huggins et al. 2016) For all , we have
Now note that corresponds to the set of points that is shattered by the affine hyperplane classifier . We can conclude that
which means that the VC dimension of is since the VC dimension of the set of hyperplane classifiers is (Kearns & Vazirani 1994; Vapnik 1995). ∎
We partition the functions into disjoint classes having equal weights. Let , for . For the sake of contradiction, suppose . Then there exists a set of size that is shattered by the ranges of . Now consider the sets , for . Due to the disjointness property, each set must be shattered by the ranges induced by . But at least one of them must be as large as , which contradicts Lemma 10. Thus follows. ∎
Let , where is an orthonormal basis for the columnspace of . It follows from and monotonicity of that
Let and . Note that and . Also,
Thus if then
If on the other hand then . Thus
From Lemma 12 and Lemma 13 we have for each
From this, the second claim follows via the Cauchy-Schwarz inequality and using the fact that the Frobenius norm satisfies due to orthonormality of . We have
The algorithm computes the QR-decomposition of . Note that is an orthonormal basis for the columnspace of . We would like to use the upper bounds on the sensitivities from Lemma 14. Namely, to sample the input points proportional to the sampling probabilities However, to keep control of the VC dimension of the involved range space, we modify them to obtain upper bounds such that each value corresponds to but is rounded up to the closest power of two. It thus holds for all . The input points are sampled proportional to the sampling probabilities From Lemma 14 we know that .
In the proof of Theorem 9, the VC dimension bound is applied to a set of functions which are reweighted by . We denote this set of functions . Now note that the sensitivities satisfy
Also note that and are fixed values. Since the values are scaled to powers of two, by (2) there can be at most distinct values of . Putting this into Lemma 11, we have .
Putting all these pieces into Theorem 9 for error parameter and failure probability , we have that a reweighted random sample of size
is a coreset with probability as claimed.
It remains to prove the claims regarding streaming and running time. We can compute the QR-decomposition of in time , see (Golub & van Loan 2013). Once is available, we can inspect it row-by-row computing and give it as input together with to independent copies of a weighted reservoir sampler (Chao 1982), which takes time to collect all sampled non-zero entries. This gives a total running time of since the computations are dominated by the QR-decomposition.
This sums up to two passes over the data and a running time of . ∎
where . Note that since the weights are non-negative, sampling and reweighting does not change the sign of the entries. This implies for and that
The claim follows by folding the constant into . ∎
Recall, due to Lemma 18, the -complexity at the -th recursion level is upper bounded by . We thus apply Theorem 15 recursively times with parameter for . First we bound the approximation ratio, which is the product of the single stages. We have
Initially all weights are equal to one. So in the first application of Theorem 15 we have . This value might grow as the weights are reassigned. However, from Inequality (2) and the discussion below it follows, that the value of can grow only by a factor of in each recursive iteration. So it remains bounded by in all levels of our recursion. Its contribution to the lower order terms given in Theorem 15 is thus bounded by
The size of the data set at recursion level satisfies
for some constant . Solving the recursion until we reach we get the following bound on . We use that for our choice we have and .
We conclude that for some constant
To reduce this even further, note that in the final iteration we do not need to preserve the -complexity. We can thus apply Theorem 15 with the original approximation parameter to obtain a coreset as claimed of size
It remains to bound the failure probability. Note that we use a factor in the sampling sizes at all stages rather than . The failure probability at each stage is thus bounded by for by adjusting constants. We can thus take a union bound over the stages to get an error probability of at most
Now recall from Theorem 15 the two pass streaming algorithm whose running time was dominated by . We can thus bound the running time of the recursive algorithm for sufficiently large by
Regarding the number of passes, note that for any , after recursion steps, the leading term in the size of the coreset is as low as , after which we may arguably assume, that the coreset fits into memory. The algorithm thus takes streaming passes over the data before it turns to an internal memory algorithm. ∎
Appendix B Material for the experimental section
Appendix C Discussion of uniform sampling
As we have discussed in the lower bounds section 3, uniform sampling cannot help to build coresets of sublinear size for worst case instances. Actually this also holds for other techniques for solving logistic regression that rely on uniform subsampling, such as stochastic gradient descent (SGD).
Note that assuming and , we have , since by construction
This implies that the approximation ratio is , which turned out very large in the experiment above, cf. Figure 3.