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 S\mathfrak{S}, is bounded by a reasonably small value S{S}, 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 kk centers for kk-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.

∙\bullet 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.

∙\bullet Our third contribution is an analysis of our sampling distribution for a parametrized class of instances we call μ\mu-complex, placing our work in the framework of beyond worst-case analysis (Balcan et al. 2015; Roughgarden 2017). The parameter μ\mu roughly corresponds to the ratio between the log of correctly estimated odds and the log of incorrectly estimated odds. The condition of small μ\mu is justified by the fact that for instances with large μ\mu, 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 μ\mu, and that our sampling scheme produces the first coreset of provably sublinear size, provided that μ\mu 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, n≫dn\gg d. 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 o(n)o(n) data points while the dependence on dd 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.

μ\mu-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 μ\mu that quantifies the complexity of compressing a given data set.

XX weighted by ww is called μ\mu-complex if μw(X)≤μ\mu_{w}(X)\leq\mu.

We conjecture that computing the value of μ(X)\mu(X) 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 μ\mu-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 22-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 x∈{0,1}nx\in\{0,1\}^{n} and Bob is given an index i∈[n]i\in[n]. The goal is to determine the value of xix_{i} 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 Ω(n)\Omega(n) 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 o(n/log⁡n)o(n/\log n) 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 si≥ςis_{i}\geq\varsigma_{i} also has provable guarantees. However, the number of samples needed depends on the total sensitivity, that is, the sum of their estimates S=∑i=1nsi≥∑i=1nςi=SS=\sum\nolimits_{i=1}^{n}s_{i}\geq\sum\nolimits_{i=1}^{n}\varsigma_{i}=\mathfrak{S}, 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 RR is sampled i.i.d. with probability pj=sjSp_{j}=\frac{s_{j}}{S} from F\mathcal{F}, ui=Swjsj∣R∣u_{i}=\frac{Sw_{j}}{s_{j}|R|} denotes the weight of a function fi∈Rf_{i}\in R that corresponds to fj∈Ff_{j}\in\mathcal{F}, and where Δ\Delta is an upper bound on the VC dimension of the range space RF∗\mathfrak{R}_{\mathcal{F}^{*}} induced by F∗\mathcal{F}^{*} that can be obtained by defining F∗\mathcal{F}^{*} to be the set of functions fj∈Ff_{j}\in\mathcal{F} where each function is scaled by Swjsj∣R∣\frac{Sw_{j}}{s_{j}|R|}.

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 tt 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 μ\mu 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 ε′=ε/μ+1\varepsilon^{\prime}=\varepsilon/\sqrt{\mu+1}.

Using this, we can show that the μ\mu-complexity is not violated too much after one stage of sampling.

Let TT be a sampling and reweighting matrix according to Theorem 15 where parameter ε\varepsilon is replaced by ε/μ+1\varepsilon/\sqrt{\mu+1}. That is TDwXTD_{w}X is the resulting reweighted sample when Theorem 15 succeeds on μ\mu-complex input X,wX,w. Suppose that simultaneously Lemma 17 holds. Let

Then we have μ′≤(1+ε)μ.\mu^{\prime}\leq(1+\varepsilon)\mu.

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 350,000350,000 unigrams with 127127 features from web pages which have to be classified as spam or normal pages (61%61\% positive). The Covertype https://archive.ics.uci.edu/ml/datasets/covertype data consists of 581,012581,012 cartographic observations of different forests with 5454 features. The task is to predict the type of trees at each location (49%49\% positive). The KDD Cup ’99 http://kdd.ics.uci.edu/databases/kddcup99/kddcup99.html data comprises 494,021494,021 network connections with 4141 features and the task is to detect network intrusions (20%20\% positive).

For each data set, we ran all three subsampling algorithms for a number of thirty regular subsampling steps in the range k∈[⌊2n⌋,⌈n/16⌉]k\in[\lfloor 2\sqrt{n}\rfloor,\lceil n/16\rceil]. 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 kk-means on all instances. This is especially true for small sampling sizes. Here, the relative error especially for uniform sampling tends to deteriorate. While kk-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 1.51.5-33 and drops more quickly towards 00. The conclusion so far could be that for a quick guess, say a 1.11.1-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 μ\mu support our experimental findings, especially that KDD Cup ’99 seems more difficult to approximate than the others. The estimated values were 4.394.39 for Webb Spam, 1.861.86 for Covertype, and 35.1835.18 for KDD Cup ’99.

The relative running time for the QR-distribution was comparable to kk-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 kk-means algorithm and the KDD Cup ’99 data set, where the uniform and kk-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 o(n/log⁡n)o(n/\log n) space. We construct the following protocol for INDEX: Consider an instance of INDEX, i.e., Alice has a string x∈{0,1}nx\in\{0,1\}^{n} and Bob has an index i∈[n]i\in[n]. We transform the instance into an instance for logistic regression. For each xj=1x_{j}=1, Alice adds a point pj=(cos⁡(jn),sin⁡(jn))p_{j}=(\cos(\frac{j}{n}),\sin(\frac{j}{n})). 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 11. 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 pi=(1−δ)⋅(cos⁡(in),sin⁡(in))p_{i}=(1-\delta)\cdot(\cos(\frac{i}{n}),\sin(\frac{i}{n})) for small enough δ>0\delta>0 with label −1-1. From the contents of Alice’s message and pip_{i}, Bob now obtains a solution to the logistic regression instance. Clearly, if Alice added pip_{i} and hence xi=1x_{i}=1 then the optimal solution will have cost at least ln⁡(2)\ln(2), since there will be at least one misclassification. If, on the other hand, Alice did not add pip_{i} and hence xi=0x_{i}=0, then the two point sets are linearly separable and the cost tends to 00. Distinguishing between these two cases, i.e. approximating the cost of logistic regression beyond a factor lim⁡x→0ln⁡(2)x\lim\limits_{x\rightarrow 0}\frac{\ln(2)}{x} 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 O(log⁡n)O(\log n) bits per point, i.e., the space of Alice’s point set is at most n′∈O(nlog⁡n)n^{\prime}\in O(n\log n). The space bound now follows from the lower bound of Ω(n)⊆Ω(n′/log⁡n)\Omega(n)\subseteq\Omega(n^{\prime}/\log n) bits due to Kremer et al. 1999 for the INDEX problem. ∎

If we had a coreset construction with o(n/log⁡n)o(n/\log n) 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 o(n)o(n) communication, which contradicts the lower bound of Kremer et al. 1999. So Alice’s coreset cannot exist. ∎

(cf. Huggins et al. 2016) For all G⊆FlogcG\subseteq\mathcal{F}^{c}_{log}, we have

Now note that {c⋅gi∈G∣xiβ≥g−1(r/c)}\{c\cdot g_{i}\in G\mid x_{i}\beta\geq g^{-1}(r/c)\} corresponds to the set of points that is shattered by the affine hyperplane classifier xi↦1{xiβ−g−1(r/c)≥0}x_{i}\mapsto\mathbf{1}_{\{x_{i}\beta-g^{-1}(r/c)\geq 0\}}. We can conclude that

which means that the VC dimension of RFlogc\mathfrak{R}_{\mathcal{F}^{c}_{log}} is d+1d+1 since the VC dimension of the set of hyperplane classifiers is d+1d+1 (Kearns & Vazirani 1994; Vapnik 1995). ∎

We partition the functions into tt disjoint classes having equal weights. Let Fi={wj⋅gj∈Flog∣wj=vi}F_{i}=\{w_{j}\cdot g_{j}\in\mathcal{F}_{log}\mid w_{j}=v_{i}\}, for i∈[t]i\in[t]. For the sake of contradiction, suppose Δ(RFlog)>t⋅(d+1)\Delta(\mathfrak{R}_{\mathcal{F}_{log}})>t\cdot(d+1). Then there exists a set GG of size ∣G∣>t⋅(d+1)|G|>t\cdot(d+1) that is shattered by the ranges of RFlog\mathfrak{R}_{\mathcal{F}_{log}}. Now consider the sets Fi∩GF_{i}\cap G, for i∈[t]i\in[t]. Due to the disjointness property, each set Fi∩GF_{i}\cap G must be shattered by the ranges induced by FiF_{i}. But at least one of them must be as large as ∣G∣t>t⋅(d+1)t=d+1\frac{|G|}{t}>\frac{t\cdot(d+1)}{t}=d+1, which contradicts Lemma 10. Thus Δ(RFlog)≤t⋅(d+1)∈O(dt)\Delta(\mathfrak{R}_{\mathcal{F}_{log}})\leq t\cdot(d+1)\in O(dt) follows. ∎

Let DwX=URD_{w}X=UR, where UU is an orthonormal basis for the columnspace of DwXD_{w}X. It follows from 0.5≤xiβ0.5\leq x_{i}\beta and monotonicity of gg that

Let K−={j∈[n]  ∣  xjβ≤−2}K^{-}=\{j\in[n]\;|\;x_{j}\beta\leq-2\} and K+={j∈[n]  ∣  xjβ>−2}K^{+}=\{j\in[n]\;|\;x_{j}\beta>-2\}. Note that g(−2)>1/10g(-2)>1/10 and g(xiβ)≤g(0.5)<1g(x_{i}\beta)\leq g(0.5)<1. Also, ∑j∈K−wj+∑j∈K+wj=W.\sum_{j\in K^{-}}w_{j}+\sum_{j\in K^{+}}w_{j}=\mathcal{W}.

Thus if ∑j∈K+wj≥12W\sum_{j\in K^{+}}w_{j}\geq\frac{1}{2}\mathcal{W} then

If on the other hand ∑j∈K+wj<12W\sum_{j\in K^{+}}w_{j}<\frac{1}{2}\mathcal{W} then ∑j∈K−wj≥12W\sum_{j\in K^{-}}w_{j}\geq\frac{1}{2}\mathcal{W}. Thus

From Lemma 12 and Lemma 13 we have for each ii

From this, the second claim follows via the Cauchy-Schwarz inequality and using the fact that the Frobenius norm satisfies ∥U∥F=∑i∈[n],j∈[d]∣Uij∣2=d\|U\|_{F}=\sqrt{\sum\nolimits_{i\in[n],j\in[d]}|U_{ij}|^{2}}=\sqrt{d} due to orthonormality of UU. We have

The algorithm computes the QR-decomposition DwX=QRD_{w}X=QR of DwXD_{w}X. Note that QQ is an orthonormal basis for the columnspace of DwXD_{w}X. 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 si∑j=1nsj=∥Qi∥2+wi/W∑j=1n(∥Qj∥2+wj/W).\frac{s_{i}}{\sum\nolimits_{j=1}^{n}s_{j}}=\frac{\|Q_{i}\|_{2}+w_{i}/\mathcal{W}}{\sum\nolimits_{j=1}^{n}(\|Q_{j}\|_{2}+w_{j}/\mathcal{W})}. However, to keep control of the VC dimension of the involved range space, we modify them to obtain upper bounds si′s_{i}^{\prime} such that each value si′/wi{s_{i}^{\prime}}/{w_{i}} corresponds to si/wi{s_{i}}/{w_{i}} but is rounded up to the closest power of two. It thus holds si≤si′≤2sis_{i}\leq s^{\prime}_{i}\leq 2s_{i} for all i∈[n]i\in[n]. The input points are sampled proportional to the sampling probabilities pi=si′/∑j=1nsj′.p_{i}={s^{\prime}_{i}}/{\sum\nolimits_{j=1}^{n}s^{\prime}_{j}}. From Lemma 14 we know that S′=∑j=1nsj′≤2S∈O(μnd)S^{\prime}=\sum\nolimits_{j=1}^{n}s^{\prime}_{j}\leq 2S\in O(\mu\sqrt{nd}).

In the proof of Theorem 9, the VC dimension bound is applied to a set of functions which are reweighted by S′wisi′k\frac{S^{\prime}w_{i}}{s^{\prime}_{i}k}. We denote this set of functions Flog\mathcal{F}_{log}. Now note that the sensitivities satisfy

Also note that kk and S′S^{\prime} are fixed values. Since the values si′/wi{s_{i}^{\prime}}/{w_{i}} are scaled to powers of two, by (2) there can be at most O(log⁡nwmax⁡wmin⁡)⊆O(log⁡(ωn))O(\log\frac{nw_{\max}}{w_{\min}})\subseteq O(\log(\omega n)) distinct values of S′wisi′k\frac{S^{\prime}w_{i}}{s^{\prime}_{i}k}. Putting this into Lemma 11, we have Δ(RFlog)∈O(dlog⁡(ωn))\Delta(\mathfrak{R}_{\mathcal{F}_{log}})\in O(d\log(\omega n)).

Putting all these pieces into Theorem 9 for error parameter ε∈(0,1/2)\varepsilon\in(0,1/2) and failure probability η=n−c\eta=n^{-c}, we have that a reweighted random sample of size

is a (1±ε)(1\pm\varepsilon) coreset with probability 1−1/nc1-1/n^{c} as claimed.

It remains to prove the claims regarding streaming and running time. We can compute the QR-decomposition of DwXD_{w}X in time O(nd2)O(nd^{2}), see (Golub & van Loan 2013). Once QQ is available, we can inspect it row-by-row computing ∥Qi∥2+wi/W\|Q_{i}\|_{2}+w_{i}/\mathcal{W} and give it as input together with xix_{i} to kk independent copies of a weighted reservoir sampler (Chao 1982), which takes O(nnz⁡(X))O(\operatorname{\mathtt{nnz}}(X)) time to collect all sampled non-zero entries. This gives a total running time of O(nd2)O(nd^{2}) since the computations are dominated by the QR-decomposition.

This sums up to two passes over the data and a running time of O(nnz⁡(X)log⁡n+poly⁡(d)log⁡n)O(\operatorname{\mathtt{nnz}}(X)\log n+\operatorname{poly}(d)\log n). ∎

where ε′≤εμ+1\varepsilon^{\prime}\leq\frac{\varepsilon}{\mu+1}. Note that since the weights are non-negative, sampling and reweighting does not change the sign of the entries. This implies for η+=∣∥(TX′β)+∥1−∥(X′β)+∥1∣\eta^{+}=|\|(TX^{\prime}\beta)^{+}\|_{1}-\|(X^{\prime}\beta)^{+}\|_{1}| and η−=∣∥(TX′β)−∥1−∥(X′β)−∥1∣\eta^{-}=|\|(TX^{\prime}\beta)^{-}\|_{1}-\|(X^{\prime}\beta)^{-}\|_{1}| that max⁡{η+, η−}≤η++η−=∣∥TX′β∥1−∥X′β∥1∣≤ε′∥X′β∥1.\max\{\eta^{+},\,\eta^{-}\}\leq\eta^{+}+\eta^{-}=|\|TX^{\prime}\beta\|_{1}-\|X^{\prime}\beta\|_{1}|\leq\varepsilon^{\prime}\|X^{\prime}\beta\|_{1}.

The claim follows by folding the constant 14\frac{1}{4} into ε\varepsilon. ∎

Recall, due to Lemma 18, the μ′\mu^{\prime}-complexity at the ii-th recursion level is upper bounded by μ(1+ε)i\mu(1+\varepsilon)^{i}. We thus apply Theorem 15 recursively l=log⁡log⁡nl=\log\log n times with parameter εi=ε2lμ+1(1+ε)i\varepsilon_{i}=\frac{\varepsilon}{2l\sqrt{\mu+1}(1+\varepsilon)^{i}} for i∈{0…l−1}i\in\{0\ldots l-1\}. 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 ω=1\omega=1. This value might grow as the weights are reassigned. However, from Inequality (2) and the discussion below it follows, that the value of ω\omega can grow only by a factor of 2n2n in each recursive iteration. So it remains bounded by ω≤(2n)log⁡log⁡n\omega\leq(2n)^{\log\log n} in all levels of our recursion. Its contribution to the lower order terms given in Theorem 15 is thus bounded by O(log⁡((2n)1+log⁡log⁡n)))⊆O(log⁡nlog⁡log⁡n).O\left(\log((2n)^{1+\log\log n}))\right)\subseteq O\left(\log n\log\log n\right).

The size of the data set at recursion level i+1i+1 satisfies

for some constant C>1C>1. Solving the recursion until we reach n0=nn_{0}=n we get the following bound on nln_{l}. We use that for our choice l=log⁡log⁡nl=\log\log n we have 2l=log⁡n2^{l}=\log n and n2−l=2log⁡n2l=2n^{2^{-l}}=2^{\frac{\log n}{2^{l}}}=2.

We conclude that for some constant C′>CC^{\prime}>C

To reduce this even further, note that in the final iteration we do not need to preserve the μ\mu-complexity. We can thus apply Theorem 15 with the original approximation parameter ε\varepsilon to obtain a coreset as claimed of size

It remains to bound the failure probability. Note that we use a log⁡n\log n factor in the sampling sizes at all stages rather than log⁡ni\log n_{i}. The failure probability at each stage is thus bounded by 1nc′\frac{1}{n^{c^{\prime}}} for c′=c+1>2c^{\prime}=c+1>2 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 O(nnz⁡(X)log⁡ni+poly⁡(d)log⁡ni)O(\operatorname{\mathtt{nnz}}(X)\log n_{i}+\operatorname{poly}(d)\log n_{i}). We can thus bound the running time of the recursive algorithm for sufficiently large C>1C>1 by

Regarding the number of passes, note that for any η>0\eta>0, after log⁡(1η)\log(\frac{1}{\eta}) recursion steps, the leading term in the size of the coreset is as low as n2−log⁡1η=nηn^{2^{-\log\frac{1}{\eta}}}=n^{\eta}, after which we may arguably assume, that the coreset fits into memory. The algorithm thus takes 2log⁡(1η)2\log(\frac{1}{\eta}) 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 β0=0\beta_{0}=0 and β1≫(1+ε)\beta_{1}\gg(1+\varepsilon), we have f(Xβ)>2nβ1≫(1+ε)mf(X\beta)>2n\beta_{1}\gg(1+\varepsilon)m, since by construction

This implies that the approximation ratio is f(Xβ)f(Xβ^)≥2nβ1m=2nβ12n+2⟶n→∞β1\frac{f(X\beta)}{f(X\hat{\beta})}\geq\frac{2n\beta_{1}}{m}=\frac{2n\beta_{1}}{2n+2}\overset{n\rightarrow\infty}{\longrightarrow}\beta_{1}, which turned out very large in the experiment above, cf. Figure 3.