Data subsampling for Poisson regression with pth-root-link

Authors: han cheng lie
alexander  munteanu
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 10:26:07

Abstract: We develop and analyze data subsampling techniques for Poisson regression, the standard model for count data $y\in\mathbb{N}$. In particular, we consider the Poisson generalized linear model with ID- and square root-link functions. We consider the method of \emph{coresets}, which are small weighted subsets that approximate the loss function of Poisson regression up to a factor of $1\pm\varepsilon$. We show $\Omega(n)$ lower bounds against coresets for Poisson regression that continue to hold against arbitrary data reduction techniques up to logarithmic factors. By introducing a novel complexity parameter and a domain shifting approach, we show that sublinear coresets with $1\pm\varepsilon$ approximation guarantee exist when the complexity parameter is small. In particular, the dependence on the number of input points can be reduced to polylogarithmic. We show that the dependence on other input parameters can also be bounded sublinearly, though not always logarithmically. In particular, we show that the square root-link admits an $O(\log(y_{\max}))$ dependence, where $y_{\max}$ denotes the largest count presented in the data, while the ID-link requires a $\Theta(\sqrt{y_{\max}/\log(y_{\max})})$ dependence. As an auxiliary result for proving the tightness of the bound with respect to $y_{\max}$ in the case of the ID-link, we show an improved bound on the principal branch of the Lambert $W_0$ function, which may be of independent interest. We further show the limitations of our analysis when $p$th degree root-link functions for $p\geq 3$ are considered, which indicate that other analytical or computational methods would be required if such a generalization is even possible.

Main paper: 
 introduction: 
 Random sampling is arguably one of the most popular approaches to reduce large amounts of data to save memory, runtime, and further downstream resources such as communication bandwidth and energy. In contrast, classic statistical learning theory often uses uniform sampling and provides only asymptotic approximation guarantees. These guarantees often require strict assumptions such as i.i.d. data and for model assumptions to be met exactly. However, data collected from real applications often violate these conditions: only finite samples are available, independence might not be satisfied, and the model may deviate from reality. When model fitting algorithms are applied to such data, we are not only interested in reducing the above-mentioned resource requirements, but also in providing rigorous worst-case guarantees on approximation.Arguably, the most popular approach is the Sensitivity Framework [26, 19] , which provides a generalpurpose importance sampling scheme that yields a weighted subsample -or coreset -that given a data matrix X approximates some loss function f (Xβ) within a factor (1 ± ε) for any query point β. This guarantee can be stated as follows: a significantly smaller subset K ⊆ X, k := |K| ≪ |X| together with corresponding weights w ∈ R k is a (1 ± ε)-coreset for X if it satisfiesEQUATIONWe point the interested reader to [37, 35, 32] for a gentle introduction and overview on coresets.Our aim is of course not only to obtain good approximation accuracy as stated in Equation ( 1), but also for the subsample achieving the bound to be sublinear in the input size. Unfortunately, most generalized linear models do not admit strongly sublinear data summaries with reasonable approximation guarantees [36, 31] . This also holds for the Poisson models considered in this paper.To go beyond the worst-case setting and enable meaningful data reduction, a natural approach is to parameterize the analysis with a quantity that captures the fit of data to the statistical model and quantifies the achievable size of succinct data summaries [36] . Another ingredient that is commonly used to tackle data reduction for generalized linear models is to relate their loss to ℓ p norms, for which ℓ p sensitivities or leverage scores yield viable importance sampling distributions [16, 12, 36, 34] .We provide the first rigorous analysis for (1 ± ε)-approximate data reduction for Poisson models:1. We show Ω(n) lower bounds against coresets for Poisson regression (Lemma 6.1), showing that changing the link function alone does not resolve the problem of bounding the complexity of coresets for the log-link, which is incompressible. Our lower bound extends to arbitrary data reduction techniques up to a log(n) factor (Lemma 6.2). 2. We introduce a novel complexity parameter ρ that captures the compressibility of data under a Poisson pth-root-link model (Equation ( 2)). This parameter corresponds naturally to the statistical model assumptions and establishes a relationship between these assumptions and an optimization perspective of the compressibility problem. 3. We conduct a parameterized analysis, showing that sublinear coresets exist under the statistically natural assumption of small ρ parameter (Theorem 3.8), and using a novel domain shift idea for their optimization (Theorem 4.2). 4. We prove a square root upper bound for the Lambert W 0 function over [-1/e, 0) (Lemma 6.3) that allows us to prove tight bounds for the slope of a linear lower envelope (Lemma 6.4). This justifies an Θ( √ y max )1 dependence for the ID-link, which is contrasted by an O(log (y max )) dependence for the square root-link. 5. We show the limitation of our domain shifting approach, showing that the error of this method cannot be bounded to give the required (1 + ε) approximation for p ≥ 3 (Lemma 6.5). This indicates a limitation to the common choices p ∈ {1, 2}, and suggests that different techniques than the ones we develop below may be needed to overcome this limitation.The general outline of our analysis follows the established method of sensitivity sampling [26, 19] . Several steps along this outline require novel ideas due to peculiarities of the Poisson loss function defined in Equations ( 3) and (4) below. A VC dimension bound of d 2 is easy to obtain by counting the number of arithmetic operations required to compare an individual loss to a given threshold as a measure of complexity [4] . A near-linear Õ(d/ε) bound is obtained by a more fine-grained analysis, by grouping and rounding the associated count data (respectively, sensitivity scores) to powers of (1 + ε) (resp. to powers of 2). This results in a surrogate loss function that admits group-wise linear VC dimension, by splitting the domain of each loss function at its unique global minimum into two regions such that the restriction of the function to each region is monotone, and connecting the resulting construction to hyperplane classifiers. Approximating the surrogate finally implies the desired (1 ± ε) approximation for the original loss as well.For bounding the sensitivities, the domain of the loss function g(z) is split into three intervals: 1) one interval consisting of 'moderate' values of z, such that the g(z) values are bounded above and below within constants and can be treated using simple uniform sampling; 2) one interval where z has large values, in which case g(z) is closely related to z p ; 3) one interval, where z is close to 0, in which case the negative logarithm dominates the loss.Tackling the interval in 2) requires relating the loss function g y (z) to z p . Specifically we would like to bound z p ≥ g y (z) ≥ z p /λ for sufficiently large z and for some value of λ. This requires special care, since the loss function g y (z) is translated polynomially towards larger values of z with growing y, but the minimum of g y (z) grows only logarithmically with y. Informally speaking, the loss function 'widens', and its minimum 'moves mainly to the right', so for large y, we would need a very flat lower bound, which requires large λ ≈ y (ignoring polylogarithmic terms). However, this is undesirable, since λ appears to be a crucial parameter for bounding the subsample size. Specifically, this would yield a i∈[n] y i = Ω(n) dependence in the coreset size. So instead, we bound the loss roughly as z p ≥ g y (z) ≥ (z -y 1/p ) p /λ, which amounts to translating the lower envelope with growing y as well. Additionally, we introduce a novel complexity parameter ρ which balances the translated (z -y 1/p ) p lower bound with the z p upper bound. We stress that these translation and balancing arguments are not artificial or just used to make calculations go through, but they are naturally consistent with the statistical model (see the discussion below Equation ( 2)). This proof strategy eventually captures the loss function within interval 2) more closely and yields sublinear bounds for λ as well.We remark that in contrast to µ-complexity in previous work [36, 34] , a bounded balancing complexity parameter ρ does not handle the asymmetry between intervals 2) and 3). Tackling the interval in 3) thus also requires completely new ideas, as the negative logarithm has an infinite asymptote at 0, which we exploit to prove Ω(n) lower bounds on subsample size in Lemma 6.1 and Lemma 6.2. Such asymptotes have not been mentioned or analyzed in previous related work on sensitivity sampling for GLMs. To circumvent the lower bound, we avoid this interval by introducing a novel domain shifting approach, requiring all feasible solutions to satisfy z > η for a suitable η > 0 for optimization. Choosing η in the order of ε, we can argue that the solution in the shifted domain is a (1 + ε) approximation. Avoiding the asymptote enables a coreset construction for the shifted domain.We believe that the domain shifting approach is necessary: if an instance consists of the extreme points on the convex hull, and all but a small (sublinear) number of points are separated by an ε distance to the boundary, then required structure is already in the data. But if non-extremal points are allowed to be arbitrarily close to the boundary, and we do not shift the domain, then we will not avoid high sensitivity points that are strictly inside the convex hull. Then the coreset size would necessarily depend on the distance of these non-extremal points to the boundary, and crucially on the number of points that are very close to the boundary of the convex hull, which again can be Ω(n).Exploiting the asymmetry between the intervals 2) and 3) where the loss exhibits z p and -log(z) growth respectively, we prove Ω(n) lower bounds on subsample size, by adapting known reduction techniques [31, 36, 41] . We also provide lower bounds on parameters used in our analysis, showing their tightness. In particular, the aforementioned λ slope parameter is of size Θ( y/ log(y)) in the case p = 1. The proof is conducted by an exact characterization of the tangent point between our linear lower envelope and the loss function. Since this requires balancing between z and log(z) and examining the log(y!) function, further bounds rely on Stirling's approximation and the principal branch of the Lambert W 0 function. Recent bounds in [38] for the Lambert W 0 function imply our square root upper bound but only a cubic root lower bound. We thus significantly tighten their bound in an appropriate interval. This result may be of independent interest, since the Lambert function cannot be expressed in terms of elementary functions and has many important applications in various fields. As a result, we obtain a matching square root lower bound on λ in our context.Classic work on data subsampling started with linear ℓ 2 regression [16] and was extended to linear ℓ p regression [12] . More recently, the study continued with non-linear transformations such as in generalized linear regression models. The first guaranteed finite subsample bounds for logistic regression appeared in [36] , while impossibility results for Poisson regression were given in [31] . Research on generalized linear models was continued for Probit regression [34] . Asymptotic properties of subsampling for generalized linear models, including for Poisson regression, were studied in [3, 27] .A finite-sample size result is given in [3, Theorem 5 ] that exhibits O( √ d) approximation error. [13] studied a sampling-based feature space reduction for a wide array of generalized linear models with additive errors. However, parts of their assumptions specifically do not apply to Poisson regression. 
preliminaries and the poisson pth-root-link model: 
 Poisson regression models aim to predict a count variable Y ∈ N 0 using a generalized linear model with link function h : R → R, i.e., h(E(Y | x)) = xβ , where x = (1, x (1) , . . . , x (d-1) ) ∈ R d is a row vector, and β ∈ R d is a column vector carrying the model parameters that in particular include an intercept β 1 [30] . Common choices for h are the canonical log-link h(v) = ln(v), the ID-link h(v) = v, and the root-link h(v) = v 1/2 . The latter two can be cast into a general framework by introducing the pth-root-link h(v) = v 1/p , for any R ∋ p ≥ 1, where the ID-link and root-link correspond to p ∈ {1, 2} [10] . Subsampling for the log-link is not possible with the multiplicative (1 + ε) error guarantees that we aim for, since it entails preserving the exp(xβ) function [31] . We will also show impossibility results for the pth-root-link. However, we parameterize our analysis with a data-dependent parameter that reflects naturally how well the realized data distribution is captured by the Poisson model. This parameter is inspired from previous work [36] : we will refer to data X, y as being 'ρ-complex', if there exists a 0 < ρ < ∞ such that denoting by x j the j-th row of the data matrix X ∈ R n×d and by y j the j-th entry of the vector y ∈ N n 0 , it holds thatEQUATIONWe may interpret the parameter ρ as follows. The rate parameter of the predicted Poisson distribution is E(Y | x) = (xβ) p . Hence, the mean and variance of Y j given x j is (x j β) p , for each j = 1, . . . , n.To obtain the maximum likelihood estimator of β, we seek β such that for every j = 1, . . . , n, x j β is as close as possible to y 1/p j , since y 1/p j minimizes the jth summand g yj of the loss function specified in Equations ( 3) and (4) below. However, choosing β so that |x j β -y 1/p j | is small will imply that each summand in the numerator of Equation (2) will be close to y 1/p j . In this case, the variance of the (Y j ) n j=1 will not be captured effectively by the Poisson model, and ρ will be large. Thus, smaller values of ρ, i.e., values of ρ that are closer to 1, indicate that the true data distribution is better captured by the Poisson model. Thus the ρ parameter in Equation (2) plays a similar role of quantifying model fit as the µ w (X) parameter from [36, Section 2] ; see in particular the comments at the end of that section.Assuming that the value of ρ is small allows us to use the proximity of the negative log-likelihood to ℓ p norms, together with some novel optimization ideas involving a shifted domain. This yields the first provable finite and sublinear subsample size with rigorous (1 + ε) approximation guarantee. We focus on the special cases p ∈ {1, 2} since they are the most popular (in fact the only practical) alternatives to the intractable log-link [10, 31] . The ID-link has been used in epidemiology [40, 29] . The root-link function has been applied to forecasting for queueing systems [39] and to account for misspecification bias in maximum likelihood estimation [17] . When the estimated mean count of the data is zero, then the canonical log-link causes problems that can be avoided by using the root-link; see e.g. [28, Section 5.4 ]. We also discuss in our lower bounds section other choices for p, and show that for any natural number p ≥ 3 the bound implied by our novel shifting idea must fail. This bound is crucial to obtain our final approximation, indicating that other methods would be required to tackle a generalization for p ≥ 3, if this is even possible. . Given a set of i.i.d. observations expressed as the rows x i of a data matrix X ∈ R n×d with corresponding labels y ∈ N n 0 we can obtain a maximum likelihood estimate of the parameter β by minimizing the negative log-likelihood, which takes the formEQUATIONwhereg yi (x i β) := (x i β) p -py i log(x i β) + log(y i !).(4) For any pth-root-link, the loss function includes a log(xβ) term, which restricts the feasible set to all β such that for all x i , i ∈ [n] it holds that x i β > 0. We note that for summands corresponding to y i = 0, the function g 0 (z) simplifies to g 0 (z) = z p > 0 with well-known properties of the ℓ p norm. We thus focus on summands g yi for y i ∈ N below.For arbitrary y ∈ N, the function g y (z) = z p -py ln(z) + ln(y) on R >0 is strictly convex with first and second derivatives g ′ y (z) = pz p-1 -py z and g ′′ y (z) = p(p -1)z p-2 + py z 2 > 0 respectively. Thus g y (z) decreases on the interval z ∈ (0, y 1/p ), increases on z ∈ (y 1/p , ∞), and has a unique minimizer at z * = y 1/p with corresponding value y -y log(y) + log(y!) ≈ 1 2 log(y) + Θ(1) ≥ 1 by Stirling's approximation. We shall use the following lower bounds to capture the y-dependence. Lemma 2.1.It holds for all z ∈ R >0 , p ∈ [1, ∞), y ∈ N that g y (z) = z p -py log(z) + log(y!) ≥ max 1, 1 3 (1 + p log(z)) .The next two results bound the individual loss contributions from above and below by roughly a value of z p . For the lower bound, however, we note that as the value of y grows, the loss function is translated polynomially towards larger z values, since its minimum is attained at z * = y 1/p . However, by Lemma 2.1, and the properties above, the increase of g y (z * ) is only logarithmic in z * , and thus also logarithmic in y. Denote by λ the scaling parameter of the lower bound on g y given in Lemma 2.2. Then the logarithmic growth of g y (z * ) implies that we would need λ ≈ y/ log(y).Unfortunately, the value of λ will affect the coreset size, which is undesirable, since it can become linear simply due to large values of y. We thus require a sublinear dependence on y max , i.e., the largest value of y presented in the data. To this end, we shift the lower envelope by the minimizer z * = y 1/p . The value of λ can subsequently be bounded in a desirable way, but differs significantly depending on the value of p: in the case p = 1 we prove λ ∈ O( y/ log(y)) to be sufficient, while the case p = 2 even constant λ = 1 will suffice. In Section 6, we will show a separation by a superconstant and matching square root lower bound on the value of λ in the dominating case p = 1.Lemma 2.2. For any p ≥ 1 and y ∈ N it holds that z p ≥ g y (z) for z > y 1/p . If p = 1, then for some λ ∈ O( y/ log(y)), it holds that g y (z) ≥ (z-y 1/p ) p λ for z > y 1/p . Lemma 2.3. Let p ≥ 2 and y ∈ N, and λ = 1. Then g y (z) ≥ (z-y 1/p ) p λ for z > y 1/p .i∈[n] (Z β) i for Z from the construction of X 3. Compute b := max{1, 2 1/p • | min i (Z β) i |} 4. Define β := (b, β) ∈ R d+1 .• Construction of y:1. Compute λ := (Xβ) p ∈ R n 2. For each i = 1, . . . , n, draw Y i ∼ Poisson(λ i ), and store the resulting vector as y.All experiments were run on a commodity machine with Intel Core i7-7700K processor (4 cores, 4.2GHz, 32GB RAM) and took overall around 50 minutes to complete. The Python code of [34] was adapted to the Poisson regression setting. 4 We applied it with the appropriate p ∈ {1, 2} to the datasets with dimensions n = 100 000, d = 7 generated as detailed in the previous section.We compared our method to uniform sampling as a baseline, which is widely popular due to its simplicity and general applicability.We varied the reduced size between 50 and 600 in equal increments of size 50. For each reduced size and each method, we performed 201 independent repetitions.Our results are shown in Figure 1 . The red (Poisson with pth-root-link) and blue (uniform sampling) solid lines display the median approximation ratio across 201 independent repetitions for each reduced size. The shaded areas below and above the solid lines indicate 2× standard errors of the respective medians.For both p ∈ {1, 2} the results look widely similar, although the case p = 1 is slightly more distinctive. We thus focus our further description on the case p = 1.Our novel Poisson subsampling method generally outperforms uniform sampling, and their 2× standard error intervals are very narrow.The shaded blue area under the blue solid line indicates that the 2× standard error for uniform sampling is slightly more narrow than the 2× standard error for 1-Poisson regression, i.e., Poisson regression with ID-link. However, this error is optimistically calculated only on repetitions that succeed in providing a valid approximation when applied to the original full data.The shaded blue area above the blue solid line is unbounded, which indicates that some of the repetitions yield solutions that are infeasible for the original full data problem, and thus fail to give an appropriate approximation. Even in a few feasible cases, approximation ratios were 1.5 -2.5 × 10 Our method is presented in red and compared against uniform sampling, which is presented in blue. Solid lines indicate the median and shaded areas indicate ±2 standard errors around the median taken across 201 independent repetitions for each reduced size between 50 and 600 in equal increment steps of 50. For the blue shaded area below the blue solid line, only feasible repetitions were counted, while the blue shaded area above represents the unbounded standard error without this restriction. For some lower reduced sizes, even the median was infinite, which results in an interrupted blue solid line. This indicates that more than half of the repetitions gave infeasible results when using uniform sampling with low sample sizes, while our method never produced infeasible results.non-negligible and we note that the solid blue line was interrupted below a reduced size of 250 (respectively 150 for p = 2) meaning that even the median was infinite, indicating that more than half of the repetitions of uniform sampling were infeasible for the original problem.In contrast to that, our method produced feasible results in all repetitions, across all reduced sizes, confirming our discussion in Appendix F, and giving approximation ratios very close to 1.• The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. Justification: Computing resources and total runtime are detailed in Appendix G. We only used standard commodity hardware on which all experiments could be carried out in around 50 minutes.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
coreset construction: 
 We begin by summarizing some key aspects of the sensitivity framework. Formal definitions for the sensitivity framework (including sensitivities, the VC dimension, and the main subsampling theorem) are given in Appendix A. In the sensitivity framework, the goal is to obtain coresets using importance sampling techniques to approximate loss functions [26] . Given a loss function whose value depends on a collection of input points, the main idea of the framework is to measure the sensitivity of any input point in terms of its worst-case contribution to the loss function. More precisely, given a point x j , its sensitivity for the loss function of the form f (Xη) = j∈[n] g(x j η) is given byσ j = sup η g(x j η) f (Xη) .The main subsampling theorem, Proposition A.5, combines the sensitivities together with the theory of VC dimension. The idea is to sample points according to probabilities that are proportional to the sensitivities, in order to create an appropriately reweighted subsample of the initial collection of points. Suppose the total sensitivity S = j∈[n] σ j of the points and the VC dimension ∆ associated with a set system based on the summands g(x i η) in the loss function f (Xη) are bounded, and choose a failure probability δ. If the subsample is of size k = O( S ε 2 (∆ log(S) + log( 1 δ ))), then it is in fact a (1 + ε)-coreset [19] with probability at least 1 -δ. Unfortunately, it is often just as difficult to compute the exact sensitivities as it is to solve the original problem. The remedy is to exploit the fact that one does not need the exact sensitivities themselves: it suffices to use any upper bounds on the sensitivities, provided that the upper bounds are not too loose, since a larger upper bound will lead to a larger coreset. Thus, we can reduce the task of coreset construction to two tasks: control of the VC dimension and sensitivity estimation of the loss function. This is handled in the following sections.We prove two different bounds on the VC dimension. The first one is a simple quadratic bound of O(d2 ). Our proof in the appendix simply counts the number of operations required to compare the loss function to a given threshold. The VC dimension bound then follows from a standard result in the context of bounding the VC dimension of neural networks [4, Thm. 8.14] . Lemma 3.1. The VC dimension of the range space associated with the class of Poisson loss functions as in Equation (3) is bounded by∆(R F * ) ≤ O(d 2 ).It is noteworthy that the quadratic dependence is implied already only from one single application of the exponential function, which is sufficient but likely not necessary in our context. Hence, we show in the remainder a more refined near-linear bound of O( d ε log(n) log(y max )) = Õ(d/ε), while keeping the dependence on other input parameters-namely, on n and y max -logarithmic.To this end, we subdivide the set of input functions into groups of growing values of their response parameter y i , and of their sensitivity ς i in a geometric progression. By rounding these values in each group to their next power in the geometric progression, we obtain disjoint sets of functions that closely approximate the original weighted loss functions, and whose VC-dimension can be bounded in O(d). Since there is only a logarithmic number of groups in both progressions, we obtain the claimed VC dimension bound.Recall the responses (y i ) i are nonnegative integers. We define their largest value (for a given input) to be y max = max{y i | i ∈ [n]}. Therefore, they are naturally bounded between 0 ≤ y i ≤ y max for all i ∈ [n] and there are at most y max + 1 different values of y i . Also note that the sensitivity values are naturally bounded by 0 ≤ ς i ≤ 1, but since they are continuous, they must be bounded away from 0 in order for the geometric progression to end in a finite (logarithmic) number of steps. If we increase each sensitivity by 1/n then the total sensitivity grows only by a constant, since S ′ = i∈[n] (ς i + 1/n) = S + 1. For these reasons, we can thus assume that 1/n ≤ ς i ≤ 1. Now, we would like to increase the sensitivities even more to their next power of 2, which will clearly increase the total sensitivity by no more than 2S ′ = 2S + 2 ≤ 3S.(5) By our above upper and lower bounds on the sensitivities, this will result in O(log(n)) groups, where in each group all sensitivities are equal. Note that in Proposition A.5, the reweighting of points depends only on fixed terms except for the sensitivities. Thus, in each group, all weights are constant. We have the following bound that applies to each group and for any fixed y ∈ N 0 . Lemma 3.2. The VC dimension of the range space induced by the set of functionsF c = {g i (β) = c • g yi (x i β) | i ∈ [n]} with equal weight c ∈ R ≥0 , and equal y i = y ∈ N 0 for all i ∈ [n] satisfies ∆ (R Fc ) = O(d).For the values of y i , we proceed in a very similar way. However, unlike the sensitivities, a constant approximation as in Equation ( 5) provided by powers of 2 will not suffice. Instead, we group the values of y i into powers of (1 + ε) and round all y i that belong to the same group to the next larger power. We argue that this preserves a (1 ± O(ε)) approximation to the original loss function. Indeed, this claim even holds for each summand g yi if y i is large enough, as the following lemma shows. Lemma 3.3. Let y ≥ 8 and 1 ≥ ε > 0. Let y < y ′ ≤ (1 + ε)y i . Then for arbitrary z > 0 it holds that (1 -3ε)g y (z) ≤ g y ′ (z) ≤ (1 + 3ε)g y (z).A direct consequence of Lemma 3.3 is that any coreset for the rounded version is a coreset for the original loss function and vice versa, up to an additional (1 ± O(ε)) error. We can therefore work with the rounded version of the loss function, which yields better bounds for the VC dimension.A general theorem for bounding the VC dimension of the union or intersection of t range spaces, each of bounded VC dimension at most D, was given in [6] . Their result yields O(tD log(t)). Here, we give a bound of O(tD) for the special case that the range spaces are disjoint 2 . Lemma 3.4. Let F be any family of functions, and let F 1 , . . . , F t ⊆ F be nonempty sets that form a partition of F, i.e., their disjoint union satisfies ˙ i∈[t] = F. Let the VC dimension of the range space induced by F i be bounded by D for all i ∈ [t]. Then the VC dimension of the range space induced byF satisfies ∆ (R F ) ≤ tD.As a result of our previous partition into groups and the O(d) bound on each group, we obtain the desired result. Lemma 3.5. Let F be the set of functions in the Poisson model. We can round and group the values of y i and the associated sensitivities ς i to obtain F * such that each function in F * is weighted by0 < w i ∈ W := {u 1 , . . . , u t } for t ∈ O(ε -1 log(n) • log(y max )). The range space induced by F * satisfies ∆ (R F * ) ≤ O( d ε log(n) log(y max )).As a direct corollary of Lemmas 3.1 and 3.5, we obtain the following combined bound. Corollary 3.6. The VC dimension ∆(R F * ) of the range space associated with the class of Poisson loss functions as in Equation ( 3) is bounded by∆ (R F * ) ≤ O d • min d, log(n) log(y max ) ε .We split the loss function into two parts:EQUATIONWe will ignore the first sum, since we will see later in the main approximation of Section 4, that by shifting the hyperplanes defined by parameter vectors in the solution space, everything can be shifted to the second sum where x i β ≥ η. In this way, we preserve a (1 + ε)-approximation, if η = Θ(ε) is small enough. We note that this shifting technique still requires the extreme points on the convex hull Ext(X) to be maintained; we address this issue in Section 5. We will focus on bounding the sensitivities for the remaining points with x i β ≥ η. In the next lemma we require the concept of a well-conditioned-basis [12] . Let q ∈ {2, ∞} denote the dual norm of p ∈ {1, 2}, respectively. We say that U is a '(α, γ, p)-well-conditioned basis' for the column span ofEQUATIONLemma 3.7. Let X ∈ R n×d , y ∈ N n 0 be a ρ-complex dataset, i.e., Equation (2) holds. Let p ∈ {1, 2}. Let λ ≥ 1 be the slope parameter from either Lemma 2.2 or Lemma 2.3 depending on the value of p. Let γ be a conditioning parameter and η > 0 be arbitrary. Then the sensitivity for each x i with x i β > η is bounded by ς i ≤ λργ p ∥U i ∥ p p + 2/n. Their total sensitivity is bounded by S ≤ O ρd y max /log(y max ) + log log(1/η) for p = 1, and S ≤ O(ρd + log (y max ) + log log(1/η)) for p = 2.Putting all steps (VC dimension, total sensitivity) together into the sensitivity framework, Proposition A.5 yields the following computational result, where in particular ℓ p well-conditioning is established constructively using ℓ 2 subspace embeddings [9] , resp. using ℓ 1 spanning sets [44] . Theorem 3.8. Let X ∈ R n×d , y ∈ N n 0 be a ρ-complex dataset, i.e., Equation (2) holds. We can compute a weighted coreset(K, w) ∈ R k×d × R k ≥0 for the pth-root-link Poisson regression problem with p ∈ {1, 2} on D(η) := {β ∈ R d : ∀i, x i β > η}. The size of the coreset is bounded by k = Õ(ε -2 d • min{d, ε -1 log(n) log(y max )} • m), where m = ρd log log(d) y max /log(y max ) + log log(1/η) p = 1 ρd + log (y max ) + log log(1/η) p = 2. 
main approximation result: 
 In the previous section, we developed a coreset for the sum of individual losses where x i β ≥ η, i.e., for the second sum of Equation ( 6). Since we cannot bound the remaining first sum where x i β < η, we choose to simply avoid it instead, by shifting each solution by η. Define for any η ≥ 0D(η) := {β ∈ R d : ∀i, x i β > η}. 3The original domain of optimization is D(0) and the shifted domain is D(η). Shifting the domain does not remove the need to store the extreme points on the convex hull of the input, since we need these points to determine the feasible domain D(η) during optimization. However, shifting removes the need to approximate the first sum in Equation ( 6) over points with unbounded sensitivity located in a small slab of width η within the convex hull.Since we have a (1 ± ε) coreset for D(η), we need to find a suitable choice for η and show that the optimizer(β ′ ) * ∈ D(η) is a (1 + O(ε)) approximation for the optimizer in the original domain β * ∈ D(0). We thus defineEQUATIONLemma 4.1. It holds for sufficiently small η > 0 thatf (Xβ * ) ≤ f (X β * ) ≤ (1 + O(η))f (Xβ * ).Now we combine the preceding results, namely the coreset and the domain shifting bound, for our main theorem. Theorem 4.2. Let ε ∈ (0, 1/14). Let (C, w) be a coreset according to Theorem 3.8. Letβ := argmin β∈D(ε) f w (Cβ), β * := argmin β∈D(0) f (Xβ). Then f (Xβ * ) ≤ f (X β) ≤ (1 + ε)f (Xβ * ). 
extreme points on the convex hull: 
 In the previous sections, we argued that for obtaining a (1 + ε) approximation, it suffices to calculate a coreset that is valid for all β ∈ D(ε), and then to minimize our loss function over β ∈ D(ε) using the coreset instead of the full data. Note that for the optimizer to stay in the feasible set D(ε), one must store the extreme points on the convex hull of the input points denoted by Ext(X). This is true even when ε = 0, i.e., even when the original function is considered and no shifting occurs. Note that there exist datasets such as our 'points on a circle' example considered in the lower bounds of Section 6, such that |Ext(X)| = n.There are several ways to either characterize |Ext(X)| for typical inputs in a sublinear way, or to approximate the convex hull by a smaller sublinear subset, called an ε-kernel, with an error of at most ε. Since these methods are usually relative to the diameter of the data, and since we need an additive error for our shifting approach, we first normalize the data to be within the unit ball. We note that this does not change the value of the loss function, since this involves scaling by a fixed value C ≥ 1, and since we can use the fact thatf (Xβ) = f X C (Cβ) ,as well as the one-to-one correspondence between any β and Cβ. Thus, we can run the algorithm on the rescaled data, obtain a good or optimal solution Cβ * , and rescale Cβ * to obtain the corresponding β * . Rescaling steps such as normalizing data to zero mean and unit variance are standard in statistical data analysis [22] .Smoothed complexity of the convex hull Instead of the worst case |Ext(X)| over all possible datasets X ∈ R n×d , in smoothed complexity we considersup X∈R n×d E Ξ∼N |Ext(X + Ξ)|,where N denotes the distribution over matrices Ξ with the same dimensions as X and with i.i.d. Gaussian entries with mean 0 and variance σ 2 , i.e., Ξ i,j ∼ N (0, σ 2 ) for all i ∈ [n], j ∈ [d] [11] . This is motivated by the fact that many datasets are recorded with measurement errors, which can often be assumed to be Gaussian. Specifically for the convex hull, [11, Chapter 4] showed that for normalized data in the unit cube, the supremum above is bounded byO log 3 2 d-1 (n) σ d + log d-1 (n) ,which is sublinear in n, though exponential in d.ε-kernels The purpose of ε-kernels is to approximate the extent of a point set up to an error of 1 -ε for any direction in R d based on a subset of the data. They were introduced by [1] and improved by [8] to optimal Θ(1/ε (d-1)/2 ) size, see the survey [2] . Since we assume our data to be normalized within the unit ball, this translates to an additive error that is bounded by ε. Thus, the boundary of the convex hull of the ε-kernel can be smaller than the boundary of the original convex hull by at most ε.Improvement for structured data It is known that ε-kernels can have size up to Ω(1/ε (d-1)/2 ) in the worst case. Beyond the worst case, the structure of the given data may allow for much smaller ε-kernels to exist, even in high dimensions. Motivated by this, [5] developed a 'greedy clustering' approach that produces an ε-kernel of size O(k min /ε 2 ), where k min = k min (X, ε) denotes the smallest possible size of a subset that gives the required ε-kernel guarantee for the original input X.Using ε-kernels for optimizing Poisson models The above options include the possibility to calculate the convex hull and rely on the sublinear smoothed complexity bound, if this is reasonable in the given context. Otherwise, if we have access to any of the above ε-kernel constructions, we shift the hyperplane away from the approximation of the convex hull, provided by the ε-kernel, by a distance of η = 2ε instead of just ε. As a result, the hyperplane will be shifted away from the original convex hull by at least ε and at most 2ε. Thus, we still compute a (1 + Θ(ε)) approximation in this way. 
lower bounds: 
 We complement our coreset constructions for the variants of the Poisson model by a series of lower bounds. Our bounds in Lemma 6.1 and Lemma 6.2 are specifically for the pth-root-link, and not for the log-link, which was studied in [31] . We use similar constructions of the bad dataset as those used frequently in previous literature, e.g., in [31, Theorem 6] , and in [36, 41, 23, 21] . However, each of these references require specific adaptations to the respective loss function that do not directly apply to our setting. Our arguments are thus adapted to the Poisson pth-root-link model to show that it does not admit sublinear coresets without imposing assumptions on the data or restricting the model.The hard instance consists of n equidistant points on the unit circle. Recall that in the Poisson regression formulation, every point has an additional intercept coordinate which is 1, and the corresponding parameter of β determines the affine translation; see Section 2. For every i ∈ [n], letx i = (1, cos( 2πi n ), sin( 2πi n )), and y i = 1. Recall that any feasible β must satisfyx i β > 0, ∀i ∈ [n].The hyperplane parameterized by β is thus always outside the point set and β points in the direction of the(x i ) i∈[n] .The idea for showing that any point has sensitivity 1 is that if β points to the center of the point set, and the hyperplane is translated to just 'touch' point x i , then x i β is arbitrarily close to 0, implying that the cost is arbitrarily large. All other points are sufficiently bounded away from the hyperplane, but also not too far away, so their cost is bounded. This implies that the sensitivity of point x i can be made arbitrarily close to 1. By symmetry of our construction, this holds for any point. Below, we prove an even more interesting statement, i.e., that no compression is possible below Ω(n) bits. While this statement appears to give a weaker Ω(n/ log(n)) bound against coresets, it in fact gives a stronger bound in some sense. This is because the bound holds against any possible data reduction algorithm and against any data structure that answers negative log-likelihood queries to within a small error, independent of what (possibly randomized) operations the data reduction algorithm performs. For example, the algorithm could subsample, it could select input points as coreset constructions do, or it could take linear combinations as in linear sketching. More generally, the bound holds against any sort of bit encoding that represents the reduced data. The reduction is based on the same data example of equidistant points on a circle embedded in d ≥ 3 dimensions. Lemma 6.2. Let Σ D be a data structure for D = [X, y] ∈ R n×d × R n , d ≥ 3, that approximates negative log-likelihood queries Σ D (β) for Poisson regression with the pth-root-link for p ∈ {1, 2}, such that for some φ ≥ 1 it holds that∀β ∈ R d : f (Xβ) ≤ Σ D (β) ≤ φ • f (Xβ). If φ < n 8 log(n) then Σ D requires Ω(n) bits of memory.Next, we prove that the parts of our analysis that are specific to the Poisson model are tight. In particular, the scale parameter λ is only a constant for p = 2, but in the case p = 1 we only have a y/log(y) upper bound. Since y max , i.e., the largest y, can potentially be very large, one may ask if we can do better. Our next result exactly characterizes the smallest possible parameter λ such that our linear lower envelope approximation is tangent to the actual loss function, in order to show a tight λ = Θ( y/log(y)) bound. This characterization of tangent points relies on properties of the principal branch W 0 of the Lambert function, which is defined by the equationW 0 (x)e W0(x) = x, for x ≥ -1/e.The only approximations that the final bounds obey follow from Stirling's approximation and from the following upper bound on the Lambert function W 0 by a square root function. This upper bound improves the recent cubic root bound of [38, Theorem 3.2] within a small region, and is crucial to obtaining a tight square root (rather than cubic root) lower bound on λ.Lemma 6.3. For all x ∈ [-1/e, 0), it holds that W 0 (x) ≤ 2(1 + ex) -1.The above bound on the W 0 function is novel and may be of independent interest. In our context, it allows us to prove the following tight bound on λ that resembles the same asymptotic upper bound as in Lemma 2.2, and establishes a matching lower bound. Lemma 6.4. Let y ∈ N be arbitrary, p = 1, and τ = y 1/p in the definition Equation (4) of g y . Let h λ (z) := (z-y 1/p ) p λ for z > 0. Then g y and h λ are tangent to each other if and only if λ = λ * (y) = (W 0 ( -y (y!) 1/y exp(2) ) + 1) -1 , in which case the unique tangent point is z * (y) = yλ * (y) λ * (y)-1 . In addition, λ * (y) = Θ( y max / log(y max )).The next lemma shows for p ≥ 3 that there exists no constant C such that the domain shifting approach that we developed in Section 4 yields a 1 + Cε error bound. As this error bound is a crucial sufficient condition for our main approximation results given in Lemma 4.1 and Theorem 4.2 to hold, the lemma suggests that different techniques may be needed. It indicates a limitation of our analysis to the most common values p ∈ {1, 2}, which are the main parameterizations considered in our work. Lemma 6.5. Let p ∈ N, p ≥ 3. Then there does not exist an absolute constant C ≥ 0 such that for all sufficiently small η > 0 and for all β ∈ D(0),EQUATION) 
concluding remarks: 
 In Section 1.3, we recalled that previous finite sample size results had either unbounded or O( √ d) error instead of our (1 + ε) approximation. Our lower bounds on the parameters, together with linear VC dimension, linear sensitivity, and linear ρ dependence in our main quantitative bounds of Theorem 3.8, leave no room for improvement (up to polylogarithmic factors) if one uses a black-box application of the sensitivity framework. We remark that recent improvements on ℓ p sensitivity sampling [33] suggest that the dimension dependence can be improved to linear as well. Exploiting the fact that our coreset gives a guarantee for all β ∈ D(ε), it would be interesting to extend the statistical treatment to the Bayesian setting, inferring the distribution of parameters over this (sub-)domain, as was recently accomplished in the case of logistic and probit regression [14] .Definition A.1 (Coreset, cf. [19] ). Let X ∈ R n×d be a set of points {x 1 , . . . , x n }, weighted by w ∈ R n >0 . For any η ∈ R d , let the cost of η w.r.t. the point x i be described by a functionw i • f (x i η) mapping from R to (0, ∞). Thus, the cost of η w.r.t. the (weighted) set X is f w (Xη) = i w i • f (w i η). Then a set K ∈ R k×d , (re)weighted by u ∈ R k >0 is a (1 + ε)-coreset of X for the function f w if k ≪ n and ∀β ∈ R d : |f w (Xη) -f u (Kη)| ≤ ε • f w (Xη) .In our analysis we use sampling based on so-called sensitivity scores, the range space induced by the set of functions, and the VC-dimension. We define these notions next. Definition A.2 (Sensitivity, [26] ). Consider a family of functions F = {g 1 , . . . , g n } mapping from R d to [0, ∞) and weighted by w ∈ R n >0 . The sensitivity of g ℓ for the functionEQUATION)The total sensitivity is S = ℓ∈[n] ς ℓ .Definition A.3 (Range space; VC dimension). A range space is a pair R = (F, ranges), where F is a set and ranges is a family of subsets of F. The VC dimension ∆(R) of R is the size |G| of the largest subset G ⊆ F such that G is shattered by ranges, i.e., |{G ∩ R : R ∈ ranges}| = 2 |G| .Definition A.4 (Induced range space). Let F be a finite set of functions mapping fromR d to R ≥0 . For every x ∈ R d and r ∈ R ≥0 , let range F (x, r) = {f ∈ F : f (x) ≥ r}, and ranges(F) = {range F (x, r) : x ∈ R d , r ∈ R ≥0 }. Let R F = (F, ranges(F))be the range space induced by F.To construct coresets for Poisson models, we use a framework that combines sensitivity scores with the theory of VC dimension, originally proposed by [18, 7] . We use a more recent and slightly updated version as stated in the following theorem. w i f i (η) - f ∈R u i f i (η) ≤ ε f ∈F w i f i (η),where each element of R is sampled i.i.d. with probability p j = sj S from F, u i = Swj |R|sj denotes the weight of a function f i ∈ R that corresponds to f j ∈ F, and where ∆ is an upper bound on the VC dimension of the range space R F * induced by F * that can be obtained by defining F * to be the set of functions f j ∈ F where each function is scaled by ∈ R >0 , p ∈ [1, ∞), y ∈ N that g y (z) = z p -py log(z) + log(y!) ≥ max 1, 1 3 (1 + p log(z)) .Proof of Lemma 2.1. For y = 1 a direct calculation yields g y (z) ≥ g y (y Since g y has a unique minimizer at z * = y 1/p , it follows that g y (z) ≥ g y (y 1/p ) = y -y log(y) + log(y!) ≥ 9 10 + 1 2 log(y) > 1 .To prove the claimed bound generalizing to all z ∈ R >0 , we first argue thatEQUATIONwhich is equivalent to showing for arbitrary y ∈ N it holds that1 3 1 + log y + 1 y = 1 3 1 + log 1 + 1 y ≤ 1 3 + 1 3y ≤ 2 3 ≤ 9 10 ,where the first inequality follows from the inequality 1 + x ≤ e x for all values of x. We see that by monotonicity of the functions h(z) = 1 3 (1 + log(z p )) and the above properties of g y (z), we have for all z ≤ (y + 1) Equation ( 11)≥ 1 3 (1 + log(y + 1)) ≥ h(z),where the first inequality above follows from the inequality one line above Equation (11) . Finally for z ≥ (y + 1) 1/p the function g y grows at least as fast as the lower bound, sinceg ′ y (z) = pz p-1 - py z = p z p-1 - y z = p z p -y z ≥ p y + 1 -y z ≥ p 3z = h ′ (z).Lemma 2.2. For any p ≥ 1 and y ∈ N it holds that z p ≥ g y (z) for z > y 1/p . If p = 1, then for some λ ∈ O( y/ log(y)), it holds that g y (z) ≥ (z-y 1/p ) p λ for z > y 1/p . Proof of Lemma 2.2. We start with the upper bound. Using the fact that the assumption z > y 1/p implies log z > 0, and setting y i = 1 in Equation (4) yieldsg 1 (z) = z p -p log(z) + log 1 < z p .It remains to prove the claim for y ≥ 2. From Stirling's approximation we getlog(y!) ≤ y log(y) -y + 1 2 log(y) + 1 2 log(2π) + 1 12 .Now, note that since the derivative of the function y → 1 2 log y -y is strictly negative over the interval [2, ∞), it follows that the function itself is decreasing over the same interval. By rearranging terms and replacing y i in Equation ( 4) with an arbitrary y ≥ 2, it follows that g y (z) = z p -y log(z p ) + log(y!) ≤ z p -y log(y) + log(y!)EQUATION≤ z p + 1 2 log(2) -2 + 1 2 log(2π) + 1 12 = z p + 1 2 log(π) - 23 12 < z p ,where the last inequality holds since π < e 2 , and thus 1 2 log(π) < log(e) = 1. In the remainder, let p = 1. Let LB(y) := max 1, 1 3(1 + log(y))be the lower bound given in Lemma 2.1. Now, we want to prove that g y (z) ≥ z-y λ =: h(z) for λ ≥ 1 as small as possible.To this end, we consider the derivatives g ′ y (z) = 1 -y z and h ′ (z) = 1 λ , and find that1 - y z ≥ 1 λ ⇐⇒ z ≥ y 1 + 1 λ -1 .Next, we see that since h(y) = 0, we can guarantee h(z) ≤ LB(y) ≤ g y (z) for all z ∈ [y, y + ∆],where ∆ := λ • LB(y).To obtain a general lower bound on g y (z), we want both conditions to hold simultaneously, which is true whenevery λ -1 ≤ ∆ ⇐⇒ y ≤ λ(λ -1)LB(y)Solving for λ yields that g y (z) ≥ z-y λ holds for all λ ≥ Proof of Lemma 2.3. First we define for τ > 0 the function h λ , and compute its derivatives:EQUATIONSetting τ = y 1/p in Equation ( 12), we obtainh λ (y 1/p ) = (y 1/p -y 1/p ) p λ = 0 < 1 2 log(2πy) < g y (y 1/p ),where the second inequality follows from the fact that the minimizer of g y is y 1/p and from Stirling's approximation:1 2 log(2πy) + 1 12y + 1 < min z>0 g y (z) = y -y log y + log y!.Thus, g y (y 1/p ) ≥ h λ (y 1/p ), and a sufficient condition for g y (z) ≥ h λ (z) to hold for z > y 1/p is that g y (z) -h(z) is nondecreasing on z > y 1/p , i.e.0 ≤ g ′ y (z) -h ′ (z) = p z p-1 - y z - (z -y 1/p ) p-1 λ ⇐⇒ z p-1 ≥ y z + (z -y 1/p ) p-1 λ . Note that z > y 1/p =⇒ y z + (z -y 1/p ) p-1 λ ≤ y y 1/p + (z -y 1/p ) p-1λ , so a sufficient condition for g ′ y (z) -h ′ λ (z) to be nonnegative for all z > y 1/p isz p-1 ≥ (z -y 1/p ) p-1 λ + (y 1/p ) p-1 , ∀z > y 1/p .If λ ≥ 1, then a further sufficient condition for g ′ y (z) -h ′ λ (z) to be nonnegative for all z > y 1/p isEQUATION)For x ∈ R, let ⌊x⌋ := sup{p ∈ Z : p ≤ x} be the floor of x. We will prove that Equation ( 13) holds.We havez y 1/p p-1 = z y 1/p p-1-⌊p-1⌋ z y 1/p ⌊p-1⌋ = z y 1/p p-1-⌊p-1⌋ z y 1/p -1 + 1 ⌊p-1⌋ = z y 1/p p-1-⌊p-1⌋ ⌊p-1⌋ m=0 ⌊p -1⌋ m z y 1/p -1 m binomial thm., p ≥ 2 ≥ z y 1/p p-1-⌊p-1⌋ z y 1/p -1 ⌊p-1⌋ + 1 z y 1/p > 1 = z y 1/p p-1-⌊p-1⌋ z y 1/p -1 ⌊p-1⌋ + z y 1/p p-1-⌊p-1⌋ ≥ z y 1/p -1 p-1-⌊p-1⌋ z y 1/p -1 ⌊p-1⌋ + z y 1/p p-1-⌊p-1⌋ r ≥ 0 ⇒ d dx x r ≥ 0 > z y 1/p -1 p-1 + 1 z y 1/p > 1,as desired. Proof of Lemma 3.1. Because we use the coreset approach, we need to consider weighted subsets of the data, and thus weighted sums n i=1 w i g yi (x i β) that allow for different weights in Equation (3). In determining the VC dimension we range over all r ∈ R, β ∈ R d , and need to check whether w i g yi (x i β) > r. Note that the log factorial terms are independent of x i β. In particular, for each r ∈ R, there exists R ∋ s := r -w i log(y i !). We can thus instead count the number of operations required for evaluating each summandw i ((x i β) p -py i ln(x i β)) ≥ s.This can be rearranged to x i β ≤ exp (x i β) p -s/w i py i which can be accomplished using (at most) the following operations• 1 division to compute s/w i • d multiplications and d -1 additions to compute x i β • p -1 multiplications of x i β to compute (x i β) p • 1 subtraction to compute (x i β) p -s/w i • 1 multiplication to compute py i (0 for p = 1)• 1 division by py i • 1 exponential function evaluation to compute exp (xiβ) p -s/wi pyi • 1 operation to verify whether x i β ≤ exp (xiβ) p -s/wi pyi holds • 1 operation to output 0 or 1 depending on whether x i β ≤ exp (xiβ) p -s/wi pyi holds for a total of t = 2d + p + 5 operations, exactly q = 1 of which is an exponential function evaluation. Putting these quantities into the final conclusion of [4, Thm. 8.14] yields ∆(R F * ) ≤ (d(q + 1)) 2 + 11d(q + 1)(t + log 2 ((9d(q + 1)))= (2d) 2 + 22d(2d + p + 5 + log 2 (18d)) = O(d 2 ).Lemma 3.2. The VC dimension of the range space induced by the set of functionsF c = {g i (β) = c • g yi (x i β) | i ∈ [n]} with equal weight c ∈ R ≥0 , and equal y i = y ∈ N 0 for all i ∈ [n] satisfies ∆ (R Fc ) = O(d).Proof of Lemma 3.2. We split the functions g i into subfunctions g ≤y 1/p i (x i β) and g >y 1/p i (x i β) depending on whether x i β ≤ y 1/p i or x i β > y 1/p i . Since g i is minimized when x i β = y 1/p i , and due to strict convexity, the two subfunctions g ≤y 1/p i : (0, y1/p i ] → [g(y 1/p i ), ∞) and g >y 1/p i : (y 1/p i , ∞) → [g(y 1/p i ), ∞)are strictly monotonic and invertible on their respective ranges and domains. Now fix an arbitrary subset G ⊆ F c . Let Ω = R d × R ≥0 . We have the following bound:|{G ∩ R | R ∈ ranges(F c )}| = |{range G (β, r) | β ∈ R d , r ∈ R ≥0 }| = (β,r)∈Ω {{g i ∈ G | g i (β) ≥ r}} = (β,r)∈Ω {{g i ∈ G | c • g ≤y 1/p i (x i β) ≥ r ∨ c • g >y 1/p i (x i β) ≥ r}} ≤ (β,r)∈Ω {{g i ∈ G | x i β ≥ g -1 ≤y 1/p i (r/c)}} • (β,r)∈Ω {{g i ∈ G | x i β ≥ g -1 >y 1/p i (r/c)}} . (14)The inequality holds, since each non-empty set in the collection on the LHS satisfies either of the conditions of the sets in the collections on the RHS, or both, and is thus the union of two of the sets, one from each collection. It can thus comprise at most all unions obtained from combining any two of these sets. Now, note that both sets are of the form {g i ∈ G | x i β ≥ s 1 } where s 1 = g -1 ≤y 1/p i (r/c) maps any real r to a value of some subset of the reals s ∈ D ⊂ R as specified above. Extending the domain of s and x i β to the reals, we obtain exactly the points that are shattered by the affine hyperplane classifier x i → 1 {xiβ-s≥0} . The VC dimension of the set of hyperplane classifiers is d + 1 [24, 42] . The argument holds verbatim for s 2 = g -1>y 1/p i (r/c).We conclude the claimed bound on ∆(R Fc ) by showing that the above term Equation ( 14) is strictly less than 2 |G| for |G| = 10(d + 1). By a bound on the growth of the sets (see [6, 24] ), we have for this particular choice(β,r)∈Ω {{g i ∈ G | x i β ≥ g -1 ≤y 1/p i (r/c)}} • (β,r)∈Ω {{g i ∈ G | x i β ≥ g -1 >y 1/p i (r/c)}} ≤ {{g i ∈ G | x i β -s 1 ≥ 0} | β ∈ R d , s 1 ∈ R} • {{g i ∈ G | x i β -s 2 ≥ 0} | β ∈ R d , s 2 ∈ R} ≤ e|G| d + 1 (d+1) 2 < 30 (d+1) 2 = 2 2(d+1) log 2 (30) ≤ 2 10(d+1) = 2 |G|which implies that ∆(R F ℓ 1 ) < 10(d + 1).Lemma 3.3. Let y ≥ 8 and 1 ≥ ε > 0. Let y < y ′ ≤ (1 + ε)y i . Then for arbitrary z > 0 it holds that (1 -3ε)g y (z) ≤ g y ′ (z) ≤ (1 + 3ε)g y (z).Proof of Lemma 3.3. Recall thatg y ′ (z) = (z) p -py ′ log(z) (i) + log(y ′ !) (ii).We bound the two parts separately. The first part (i) is straightforward from the assumption on y ′ . We only need to distinguish two cases, depending on which either the upper or the lower bound of the assumption y < y ′ ≤ (1 + ε)y i applies: If log(z) ≥ 0, then there is nothing to prove since in that casez p -py ′ log(z) ≤ z p -py log(z),given the hypothesis that y < y ′ . Suppose that log(z) < 0. Thenz p -py ′ log(z) ≤ z p -(1 + ε)py log(z) ≤ (1 + 3ε)(z p -py log(z)),where the first inequality follows from the hypothesis that y ′ ≤ (1 + ε)y i and the second inequality follows since 1 + 3ε > 1 + ε > 1.For the second part (ii) we need some technical claims. Since y, y ′ ∈ N 0 , it must hold that y ′ ≥ y + 1. Then1 12y ′ ≤ 1 12(y + 1) ≤ 1 12y + 1 .Further, for any x ≥ e it holds thatlog((1 + ε)x) = log(x) + log(1 + ε) ≤ log(x) + ε • 1 ≤ (1 + ε) log(x),where the last inequality uses the condition that x ≥ e.With this in place, and using y ′ > y ≥ 8 > e 2 , we apply Stirling's approximation with a constant C := 1 2 log(2π) that is independent of y, y ′ , and that is common to both the upper and lower bounds from Stirling's approximation for the factorial function. We thus obtainlog(y ′ !) ≤ y ′ log(y ′ ) -y ′ + 1 2 log(y ′ ) + 1 12y ′ + C = y ′ (log(y ′ ) -1) + 1 2 log(y ′ ) + 1 12y ′ + C = y ′ (log(y ′ /e)) + 1 2 log(y ′ ) + 1 12y ′ + C ≤ (1 + ε) 2 y log(y/e) + (1 + ε) 1 2 log(y) + 1 12y + 1 + C ≤ (1 + ε) 2 y log(y/e) + 1 2 log(y) + 1 12y + 1 + C ≤ (1 + 3ε) log(y!)where the first inequality follows from the upper bound in Stirling's approximation for the factorial, the second inequality follows from the hypothesis that y ′ ≤ (1 + ε)y and the inequality 1 12y ′ ≤ 1 12y+1 proven above, and the last inequality follows from the lower bound in Stirling's approximation of the factorial and from the fact that (1 + ε) 2 = 1 + 2ε + ε 2 ≤ 1 + 3ε. The lower bound can be treated in a similar way. Overall, our claim follows.Lemma 3.4. Let F be any family of functions, and let F 1 , . . . , F t ⊆ F be nonempty sets that form a partition of F, i.e., their disjoint union satisfies ˙ i∈[t] = F. Let the VC dimension of the range space induced by F i be bounded by D for all i ∈ [t]. Then the VC dimension of the range space induced by F satisfies ∆ (R F ) ≤ tD.Proof of Lemma 3.4. We prove the claim by contradiction. To this end suppose the VC dimension for F is strictly larger than tD. Then there exists a set G of size |G| > tD that is shattered by the ranges of R G . Consider its intersections G i = G ∩ F i , i ∈ [t] with the sets F i . By their disjointness, each G i must be shattered by the ranges of R Fi . Note, that at least one G i must therefore satisfy |G i |/t > D, which contradicts the assumption that their VC dimension is bounded by D. Our claim thus follows. Lemma 3.5. Let F be the set of functions in the Poisson model. We can round and group the values of y i and the associated sensitivities ς i to obtain F * such that each function in F * is weighted by 0 < w i ∈ W := {u 1 , . . . , u t } for t ∈ O(ε -1 log(n) • log(y max )). The range space induced byF * satisfies ∆ (R F * ) ≤ O( d ε log(n) log(y max )).Proof of Lemma 3.5. We partition our input functions g yi , i ∈ [n] into disjoint sets with boundaries that increase in a geometric progression, depending on the sensitivities resp. weights, and on the response valuesG ij = {g y k | 8 • (1 + ε) i ≤ y k < ⌊8 • (1 + ε) i+1 ⌋, 2 j ς min ≤ ς k < 2 j+1 ς min }, i ∈ [0, O(log 1+ε (y max ))], j ∈ [0, O(log(n))].Additionally, we put the remaining values of y into a constant number of disjoint setsH yj = {g y k | y k = y, 2 j ς min ≤ ς k < 2 j+1 ς min }, y ∈ {0, 1, 2, . . . , 7}, j ∈ [0, O(log(n))].In particular, we note thatF = ˙ ij G ij ∪ ˙ yj H yjforms a partition of the whole function family, since the sets are disjoint and cover all functions by construction.Each member of a set is of the same form, i.e., after rounding the weights and y i , all members of a subset have equal y i = y, and they have equal weight. The assumptions are thus satisfied for invoking Lemma 3.2 to bound the VC dimension for each of the induced range spaces by O(d). By construction, the subsets are disjoint and their number is bounded byt = O(log 2 (n) • log 1+ε (y max )) = O log(n) • log(y max ) log(1 + ε) = O(ε -1 log(n) log(y max )).We can thus invoke Lemma 3.4 to obtain∆ (R F * ) ≤ O(dt) = O d ε log(n) log(y max ) = Õ d ε .Recall the condition Equation ( 7) for an (α, γ, p)-well-conditioned basis.Lemma 3.7. Let X ∈ R n×d , y ∈ N n 0 be a ρ-complex dataset, i.e., Equation (2) holds. Let p ∈ {1, 2}. Let λ ≥ 1 be the slope parameter from either Lemma 2.2 or Lemma 2.3 depending on the value of p. Let γ be a conditioning parameter and η > 0 be arbitrary. Then the sensitivity for each x i with x i β > η is bounded by ς i ≤ λργ p ∥U i ∥ p p + 2/n. Their total sensitivity is bounded by S ≤ O ρd y max /log(y max ) + log log(1/η) for p = 1, and S ≤ O(ρd + log (y max ) + log log(1/η)) for p = 2.Proof of Lemma 3.7. Case 1: y i = 0: We start with the special case where y i = 0. Recall that x i β > 0. Then,g yi (x i β) = (x i β) p ≤ ∥U i ∥ p p ∥Rβ∥ p q ≤ ∥U i ∥ p p (γ∥U Rβ∥ q ) p = γ p ∥U i ∥ p p ∥U Rβ∥ p p = γ p ∥U i ∥ p p ∥Xβ∥ p p = γ p ∥U i ∥ p p n j=1 g yj (x j β).Next, we consider 0 ̸ = y i ∈ N, and divide this into two sub-cases.Case 2(i):y i ∈ N, η < x i β ≤ y 1/p i :We start with the case 0 < η ≤ x i β ≤ y 1/p i . For g yi defined in Equation ( 4),g yi (x i β) := (x i β) p -py i log(x i β) + log(y i !)it holds using Lemma 2.1, the bounds on x i β and the monotonicity of g y (z) in the interval that 1 ≤ g yi (x i β) ≤ η p + y i log 1 η p + log(y i !) =: UB (y i ). Let G = {i ∈ [n] | 1 ≤ g yi (x i β) ≤ UB (y max )}. We subdivide G = ˙ l j=1 G j into disjoint sets G j = {i ∈ G | UB (y max ) • 2 -j < g yi (x i β) ≤ UB (y max ) • 2 -j+1 }, j ∈ {1, . . . ,ς i := sup β g yi (x i β) n i=1 g yi (x i β) ≤ sup β g yi (x i β) i∈Gj g yi (x i β) ≤ UB (y max ) • 2 -j+1 n j UB (y max ) • 2 -j = 2 n j Summing over i ∈ G yields i∈G ς i = l j=1 i∈Gj 2 n j = l j=1 2n j n j = 2l ≤ O log(y max ) + log log y max η ≤ O log(y max ) + log log 1 η .Case 2(ii): y i ∈ N, x i β > y 1/p i : Now we take care of the remaining region where x i β > y 1/p . In particular, by Lemmas 2.2 and 2.3 for some scaling λ = λ p that depends on p ∈ {1, 2}, we have that(z -y 1/p i ) p λ ≤ g yi (z) ≤ z p in that region.Let U R be a decomposition of X, so that x i β = U i Rβ, and U is again a p-well conditioned basis, in the sense of Equation ( 7). Now, using our assumption given in Equation ( 2), we have the following inequalitiesg yi (x i β) ≤ (x i β) p ≤ ∥U i ∥ p p ∥Rβ∥ p q ≤ γ p ∥U i ∥ p p ∥U Rβ∥ p p = γ p ∥U i ∥ p p n j=1 (x j β) p ≤ ργ p ∥U i ∥ p p n j=1 (x j β -y 1/p j ) p ≤ λργ p ∥U i ∥ p p n j=1 g yj (x j β) .Summing over all sensitivities, we get that the total sensitivity is bounded byS ≤ ρλ(αγ) p + O log(y max ) + log log y max η .For the first summand, there exists for each p ∈ [1, 2], a so-called Auerbach basis attaining α = d, γ = 1, see [43, Lemma 2.22] . For the special case p = 2, we even have that any orthonormal basis satisfiesα = √ d, γ = 1 since ∥U ∥ F = √ d, and for any z ∈ R d : ∥U z∥ 2 = √ z T U T U z = √ z T z = ∥z∥ 2 .Thus, in both cases we have suitable bases with (αγ) p = d.For p = 1 Lemma 2.2 yields λ ≤ O( y max /log y max ). With this, the overall bound simplifies to S ≤ O ρd y max /log y max + log log(1/η) .For p = 2 we have that λ = 1 suffices by Lemma 2.3. With this, the overall bound simplifies to S ≤ O(ρd + log (y max ) + log log(1/η)).Theorem 3.8. Let X ∈ R n×d , y ∈ N n 0 be a ρ-complex dataset, i.e., Equation (2) holds. We can compute a weighted coreset(K, w) ∈ R k×d × R k ≥0 for the pth-root-link Poisson regression problem with p ∈ {1, 2} on D(η) := {β ∈ R d : ∀i, x i β > η}. The size of the coreset is bounded by k = Õ(ε -2 d • min{d, ε -1 log(n) log(y max )} • m), where m = ρd log log(d) y max /log(y max ) + log log(1/η) p = 1 ρd + log (y max ) + log log(1/η) p = 2.Proof of Theorem 3.8. We put our bounds from Corollary 3.6 and Lemma 3.7 together into the main theorem of the sensitivity framework Proposition A.5. That is, we calculate the sensitivity upper bounds s i , take a sample according to the distribution p i = s i / n j=1 s j of the respective size, and reweight them accordingly. Then, the calculated bounds on the VC dimension and total sensitivity yield a bound on the required size, such that Proposition A.5 yields with constant probability that the reweighted subsample gives a (1 ± η)-approximation uniformly over β ∈ D(η).Since the Auerbach basis used in the sensitivity calculations of Lemma 3.7 can be expensive to compute depending on the value of p, we use more efficient approximation techniques here.In the case p = 2, we use a sparse oblivious ℓ 2 subspace embedding by [9] , which was explicitly proven in [34, Lemma 2.14] , to give a ( √ 2d, √ 2, 2)-well-conditioned basis. This is within absolute constant factors to the ( √ d, 1, 2)-conditioning of the Auerbach basis. Thus, the complexities given in Lemma 3.7 do not change in O-notation.In the case p = 1, we use a more recent technique introduced in [44] , called (α, γ, p)-well-conditioned spanning sets. This is a relaxation of the well-conditioned basis U given in Equation (7) , where U ∈ R n×s and z ∈ R s , are allowed to have slightly increased dimension s > d. We also note that we only need to bound norms of vectors of the form Xβ in the columnspan of the data matrix whose rank is bounded by d. We thus require the bounds to hold only for y ∈ R s , s > d that actually represent vectors Xβ in a different basis. Other aspects of Equation ( 7) remain unchanged.Our proof is nearly verbatim to [44, Theorem 1.11] . Their algorithm constructs a matrix R ∈ R d×s . We set U = XR ∈ R n×s . By [44, Lemma 4.1] , this can be done with s = O(d log log(d)), such that each column U (i) , for i ∈ [s] satisfies ∥U (i) ∥ p = 1, and for every ∥Xβ∥ p = 1 there exists a vector y ∈ R s , such that Xβ = U y and ∥y∥ 2 = O(1).For p = 1, this yields that U is an (α, γ, 1)-well-conditioned spanning set, where α = O(d log log(d)), and γ = O(1). For the bound on α, it holds that∥U ∥ 1 = s i=1 ∥U (i) ∥ 1 = s = O(d log log(d)).For the bound on γ, note that the Hölder dual for p = 1 is q = ∞. Now, it follows for every y ∈ R s that represents any Xβ as a linear combination of columns of U that,∥y∥ ∞ ≤ ∥y∥ 2 ≤ O(1) = O(1)∥Xβ∥ 1 = O(1)∥U y∥ 1 .As a consequence in the case p = 1, this computational result replaces the d factor from the Auerbach basis in the proof of Lemma 3.7 by a factor (αγ) p = O(d log log(d)), as we have claimed.Recall Equation ( 8):β * := argmin β ′ ∈D(η) f (Xβ ′ ), β * := argmin β∈D(0) f (Xβ). Lemma 4.1. It holds for sufficiently small η > 0 that f (Xβ * ) ≤ f (X β * ) ≤ (1 + O(η))f (Xβ * ).Proof of Lemma 4.1. Recall from Section 2 that we choose X to be the design matrix that includes an intercept, i.e., every row of X is of the form x = (1, x (1) , x (2) , . . . , x (d-1) ). Note that by definition we have that D(η) ⊂ D(0) is a proper subset. Thus, f (Xβ * ) ≤ f (X β * ) follows immediately.Next, define for every β ∈ D(0) its shifted version to be in one-to-one correspondence with a unique β ′ ∈ D(η) via the translationEQUATION) where e 1 is the first standard basis vector. Recall that we choose the first column of X to be (1, . . . , 1) ∈ R n . Now consider (β * ) ′ ∈ D(η) to be the shifted version of the global optimizer β * ∈ D(0). Since β * minimizes the loss function overD(η), it follows that f (X β * ) ≤ f (X(β * ) ′ ).Finally, we claim that there exists an absolute constant C ≥ 0 such that for all sufficiently small η > 0 we have for all β that Equation ( 9) holds, i.e.,f (Xβ ′ ) ≤ f (Xβ) + η p n + ηCf (Xβ).By summing the result of Lemma 2.1 over all n inputs, we have that f (Xβ * ) ≥ n. Applying our claim to β * thus yieldsf (X(β * ) ′ ) ≤ f (Xβ * ) + η p n + ηCf (Xβ * ) ≤ (1 + η + Cη)f (Xβ * ) = (1 + O(η))f (Xβ * ).This proves that Lemma 2.1 and Equation ( 9) imply the upper bound of the lemma.It remains to prove our claim of Equation ( 9) for p = 1, and p = 2 separately.Case p = 1: We havef (Xβ ′ ) = i∈[n] (x i β + η) -y i log(x i β + η) + log(y i !) ≤ i∈[n] x i β -y i log(x i β) + log(y i !) + η = f (Xβ) + ηn, which satisfies Equation (9) with C = 0. Case p = 2: We have f (Xβ ′ ) = i∈[n] (x i β) 2 + 2ηx i β + η 2 -2y i log(x i β + η) + 2y i log(x i β) -2y i log(x i β) + log y i ! = i∈[n] (x i β) 2 + 2ηx i β + η 2 -2y i log x i β + η x i β -2y i log(x i β) + log y i ! =f (Xβ) + i∈[n] 2ηx i β + η 2 -2y i log x i β + η x i β =f (Xβ) + i∈[n] 2ηx i β + η 2 -2y i log 1 + η x i β =f (Xβ) + nη 2 + 2η i∈[n] x i β - y i η log 1 + η x i β .Now using that log(1 + x) ≥ x 1+x , ∀x > -1, we bound the error for every y ∈ N by a function ϕ = ϕ y .EQUATIONThe first and second derivatives are given byϕ ′ (z) = 1 + y (z + δ) 2 ≥ 1, ϕ ′′ (z) = - 2y(z + δ) (z + δ) 4 < 0,from which we know that the function is monotonically increasing and concave. On the other hand, we know that g y (z) is convex, monotonically increasing on z ∈ [y 1/2 , ∞) and is bounded below by 1. We can thus show the claim for C = 3 by comparing the functions as well as their derivatives at z = y 1/2 + 1.First note that by monotonicity, we have for all z ≤ y 1/2 + 1 thatϕ(z) ≤ ϕ(y 1/2 + 1) = (y 1/2 + 1) 2 + η(y 1/2 + 1) -y y 1/2 + 1 + η = y + 2y 1/2 + 1 + ηy 1/2 + η -y y 1/2 + 1 + η = y 1/2 + 1 + η + (1 + η)y 1/2 y 1/2 + 1 + η = 1 + (1 + η)y 1/2 y 1/2 + 1 + η ≤ 1 + (1 + η)y 1/2 y 1/2 = 1 + 1 + η ≤ 3 • 1 ≤ 3g(z).In particular this holds for z = y 1/2 + 1 as well.It remains to compare the derivatives for the choice of z = y 1/2 + 1. We haveϕ ′ (y 1/2 + 1) = 1 + y (y 1/2 + 1 + δ) 2 ≤ 1 + y (y 1/2 ) 2 = 2,which implies that we also have3g ′ (y 1/2 + 1) = 3 • 2 (y 1/2 + 1) 2 -y y 1/2 + 1 = 3 • 2 y + 2y 1/2 + 1 -y y 1/2 + 1 = 3 • 2 2y 1/2 + 1 y 1/2 + 1 ≥ 3 • 2 > 2 ≥ ϕ ′ (y 1/2 + 1).This completes the proof of ϕ(z) ≤ 3g(z) for all z by convexity of g and concavity of ϕ.Overall, the claim of Equation ( 9) follows with C = 3(p -1), for both p ∈ {1, 2}, which also concludes the proof of the lemma. f (Xβ * ) ≤ f (X β) ≤ (1 + ε)f (Xβ * ).Proof of Theorem 4.2. We invoke Lemma 4.1 with η = ε to show that the shifted version of β * , i.e., β good := (β * ) ′ = β * + ηe 1 is a (1 + ε)-approximation and β good ∈ D(η). Thus, the optimizer (β ′ ) * ∈ D(η) cannot be worse than a (1 + ε)-approximation.The coreset construction of Theorem 3.8 works uniformly over D(η). It thus yields a coreset C ⊂ X of size k with weights w ∈ R k such that if we denote the weighted loss on the coreset byEQUATIONThen definingβ good := β * + ηe 1 , we have f (X β) ≥ f (Xβ * ) since D(η) ⊂ D(0). Moreover, assuming 0 < ε ≤ 1 2 we have f (X β) ≤ 1 1 -ε f w (C β) ≤ 1 1 -ε f w (Cβ good ) ≤ 1 + ε 1 -ε f (Xβ good ) ≤ (1 + ε) 2 1 -ε f (Xβ * ) ≤ (1 + 7ε)f (Xβ * )rescaling ε finishes our main result. Proof of Lemma 6.1. We first note that our construction can be embedded arbitrarily in d ≥ 3 dimensional spaces. For simplicity, we describe the construction for d = 3, where the first dimension corresponds to the affine translation and the other two describe the location in the 2-dimensional subspace.Recall that our point set is given byx i = (1, cos( 2πi n ), sin( 2πi n )), i ∈ [n], and y i = 1 for every i ∈ [n]. By symmetry of the construction, it suffices to analyze w.l.o.g. the sensitivity of point x n = (1, cos(2π), sin(2π)) = (1, 1, 0). Since the sensitivity is defined as the supremum over all β, it also suffices to find one β for which the sensitivity is arbitrarily close to 1. To this end, for a small η > 0 yet to be determined, we chooseβ = (1 + η, -cos(2π), -sin(2π)) T = (1 + η, -1, 0) T ,where 1 + η represents a translation term and (-1, 0) T represents a 'normal' term that lives within the 2-dimensional subspace mentioned above. This normal term defines a hyperplane H (which is in fact a line) within the 2-dimensional subspace. The normal points towards the center of the point set, and the hyperplane H is at distance exactlyx n β = 1 + η -1 = η from x n .A simple trigonometric calculation yields that the separation between x n and the neighboring points x 1 and x n-1 along the direction orthogonal to H is exactly 1 -cos(2π/n). Since n ≥ 8, it holds that (2π/n) 2 /3 ≤ 1 -cos(2π/n) ≤ (2π/n) 2 /2 by a second order Taylor series expansion of the cosine function. All other points are even farther away from x n than x 1 and x n-1 , and therefore also further from H. Also note that if we let x ′ n = (1, -1, 0) be the antipodal point of x n on the circle, we see that the distances of all points from H are less than x ′ n β = 1 + η + 1 = 2 + η < 3. Recall that for arbitrary p ≥ 1, the function g y defined in Equation ( 4) is minimized at y 1/p , which in this case equals y i = 1. We have that roughly half of the points are at distance at least 1 + η and distance at most 3 from H. By strict convexity, we have that g 1 is also strictly increasing on the interval [1, ∞). We can thus upper bound the contribution of each of these points by at mostEQUATIONFor the other half of the points (except x n ) the contribution is upper bounded by the loss that occurs closest to H. By strict convexity again, we have that g 1 is also strictly decreasing on the interval (0, 1]. We argued that the points are sufficiently separated, so we get that each of their contributions is bounded byEQUATIONNow, choosing η = exp(-n 2 ), we have that the cost of the point x n is lower bounded byg 1 (x n β) = g 1 (η) ≥ 1 exp(n 2 ) p + p log(exp(n 2 )) ≥ n 2 . =⇒ -log( 1 -τ 2 ) ≤ -log(-x) -1by Equation ( 20)⇐⇒ 1 + log(-x) ≤ log( 1 -τ 2 ).The last inequality is equivalent to1 ≤ log √ 1 -τ 2 -x ⇐⇒ e ≤ √ 1 -τ 2 -x ⇐⇒ -ex ≤ 1 -τ 2 . Now since ( √ 1 -τ 2 ) 2 = 1 -τ 2 ≤ 1 -τ 2 + ( τ 2 2 ) 2 = (1 -τ 2 2 ) 2 for every τ ∈ (-1, 0], it follows that -ex ≤ 1 - τ 2 2 ⇐⇒ τ 2 2 ≤ 1 + ex ⇐⇒ |W 0 (x) + 1| ≤ 2(1 + ex),where the rightmost equivalence above follows from the definition τ := -(W 0 (x)+1). The inequality on the right-hand side above implies the desired conclusion.Recall the function g yi defined in Equation ( 4).Lemma 6.4. Let y ∈ N be arbitrary, p = 1, and τ = y 1/p in the definition Equation (4) of g y . Let h λ (z) := (z-y 1/p ) p λ for z > 0. Then g y and h λ are tangent to each other if and only if λ = λ * (y) = (W 0 ( -y (y!) 1/y exp(2) ) + 1) -1 , in which case the unique tangent point is z * (y) = yλ * (y) λ * (y)-1 . In addition, λ * (y) = Θ( y max / log(y max )).Proof of Lemma 6.4. A point of tangency z of the curves g y and h λ is defined as a point where the functions agree and their derivatives agree. To identify the point where the derivatives of g y and h λ agree, we observeg ′ y (z) = 1 - y z = 1 λ = h ′ λ (z) ⇐⇒ 1 - 1 λ = y z ⇐⇒ z = λy λ -1 .Since g y (z) ≥ 1 2 log(2πy) and h λ (z) < 0 for z < y if τ = y, the tangent point cannot lie in the interval (0, y]. Hence, z > y must hold. Combining this observation with the equation z = λy λ-1 for the point where the derivatives agree, we conclude that λ > 1.Now suppose that z is a point where the functions g y and h λ agree:h λ (z) =g y (z) ⇐⇒ y log(z) -z + z λ = y λ + log(y!) ⇐⇒ y log(z) - z λ -1 λ = y λ + log(y!) ⇐⇒ log(z) -z λy λ -1 -1 = 1 λ + 1 y log(y!) ⇐⇒ z exp -z λy λ -1 -1 = exp 1 λ (y!) 1/y ⇐⇒ -z λy λ -1 -1 exp -z λy λ -1 -1 = - λy λ -1 -1 exp 1 λ (y!) 1/y .If z is a point of tangency, then by the equivalent condition above for g ′ y = h ′ λ , it follows that z = λy λ-1 . Substituting this into the last equation above yields-1 exp(-1) = - λy λ -1 -1 exp 1 λ (y!) 1/y ⇐⇒ - y (y!) 1/y exp(-1) = - λ -1 λ exp 1 λ ⇐⇒ - y (y!) 1/y exp(-2) = 1 λ -1 exp 1 λ -1 .The last equation above yields thatEQUATIONWe can further specify the branches of the Lambert W function as follows. By definition of the factorial, 1 ≤ y (y!) 1/y for all y ∈ N. On the other hand, by Stirling's approximation,EQUATIONHence, -exp(-1) < -y (y!) 1/y exp(-2) ≤ -exp(-2), y ∈ N. This implies that the argument of the Lambert W function in Equation ( 21) lies in the interval (-e -1 , -e -2 ]. By definition of the Lambert W function, this in turn implies that in Equation ( 21), we need to consider only k = 0 and k = -1, which means that1 λ * = 1 λ * (y) ∈ W k -y (y!) 1/y exp(2) + 1 : k ∈ {0, -1}For all x ∈ [-exp(-1), -exp(-2)], W 0 (x) ≥ W -1 (x), with equality holding only for x = -exp(-1). This follows from the definition of W 0 as the principal branch of the Lambert W function. Hence,λ * (y) = 1 W 0 -y (y!) 1/y exp(2) + 1 , y ∈ N,which proves the first statement of Lemma 6.4.Next, we show that λ * (y) = Θ( y max / log(y max )). We use a lower bound for the principal branch W 0 of the Lambert W function for negative arguments from [38, Theorem 3.2] :EQUATIONSince we showed above that the argument of the Lambert W function in Equation ( 21) lies in the interval (-e -1 , -e -2 ], we may apply Equation (23) .Combining Lemma 6.3 with Equation ( 23), implies that for all x ∈ [-e -1 , 0],√ ex + 1 -1 ≤ W 0 (x) ≤ 2(1 + ex) -1 ⇔ 1 2(1 + ex) ≤ 1 W 0 (z) + 1 ≤ 1 (1 + ex) .By Equation ( 22), we may thus set x = -y (y!) 1/y e -2 in the above inequalities, which yields (24) for all y ∈ N. Again by Equation ( 22),2 -1/2 -y (y!) 1/y e -1 + 1 -1/2 ≤ 1 W 0 -y (y!) 1/y e -2 + 1 = λ * (y) ≤ -y (y!) 1/y e -1 + 1 -1/21 -exp - 1 12y 2 + y - 1 2y log(2πy) < -y (y!) 1/y e -1 + 1 < 1 -exp - 1 12y 2 - 1 2ylog(2πy) . < 1 2 at y = 1, it suffices to consider the quantity 1 -exp(-x) for x ∈ [0, 1 2 ]. By Taylor's approximation, we haveEQUATIONand applying the lower bound in Equation (25) to the lower bound for -y (y!) 1/y e -1 + 1 below Equation ( 24) yields -y (y!) where the rightmost inequality follows since the function y → 6y log(2πy) is increasing on [1, ∞) and has a value strictly larger than 1 at y = 1. Applying the inequality above to the lower bound for λ * (y) in Equation ( 24) yieldsλ * (y) ≥ 1 2 y 2 y log(2πy) 1/2 = Ω y log(y), ∀y ∈ N.This completes the proof that λ * (y) = Θ( y max / log(y max )).Lemma 6.5. Let p ∈ N, p ≥ 3. Then there does not exist an absolute constant C ≥ 0 such that for all sufficiently small η > 0 and for all β ∈ D(0), β ′ := β + ηe 1 ∈ D(η) satisfiesEQUATIONProof of Lemma 6.5. First note thatf (Xβ ′ ) = i∈[n] (x i β + η) p -py i log(x i β + η) + log(y i !) = i∈[n] (x i β) p + η p + p-1 ℓ=1 p ℓ (x i β) ℓ η p-ℓ -py i log(x i β + η) + log(y i !) =f (Xβ) + η p n + i∈[n] p-1 ℓ=1 p ℓ (x i β) ℓ η p-ℓ -py i log x i β + η x i β ,where the last equation follows by the definition of f (Xβ), given the hypothesis on p. Define for every y ∈ N the auxiliary functionEQUATION)Then the inequality Equation ( 9) is equivalent toi∈[n] φ 1,yi,p (x i β) ≤ C i∈[n] g yi (x i β) ⇐⇒ 0 ≤ i∈[n]Cg yi (x i β) -φ 1,yi,p (x i β), which implies the following statement: there exists C, η * > 0 such that for every n ∈ N, (y i ) i∈[n] ∈N n , η ≤ η * , andEQUATIONproblem formulated over D(η). But as we prove in Appendix D, we can construct a coreset that holds for all β ∈ D(η), and D(η) contains at least one β that is a (1 + ε)-approximation for the optimal solution β * of the problem on the original domain D(0), and evaluated on the full dataset. These two parts are combined to prove that the final minimizer β ∈ D(η) optimized on the coreset is a (1 + O(ε))-approximation compared to the value of β * , when both are evaluated on the full dataset.In the other direction, no infeasible solution can become feasible because of the proper subset relation D(η) ⊂ D(0). Finally, note that for any data and any fixed η < ∞, both, D(η) and D(0) are non-empty, since they consist of all β that parameterize hyperplanes that put the convex hull of input points (respectively, the additive η-inflation of the convex hull of input points) in the positive open halfspace. Thus there always exist feasible solutions, which means that no instance can become completely infeasible by means of our methods.G Pseudocode, data and experimental resultsHere we give pseudocode for our coreset construction Algorithm 1 and for the subsequent optimization procedure Algorithm 2:Algorithm 1 Coreset algorithm for pth-root-link Poisson regression.Input: data X ∈ R n×d , Y ∈ N n 0 , number of rows k (see Theorem 3.8). Output: coreset C = (X ′ , Y ′ , w) ∈ R k ′ ×d × N k ′ × R k ′ with k ′ = k + |Ext(X)| rows 1: Let X CH be the extreme points Ext(X) on the convex hull of X or their ε-kernel approximation (cf. Section 5), let Y CH be their corresponding labels 2: Assign weight vector w CH = 1 corresponding to all points in Sketch the data to obtain X = ΠX (see [9] ) 9:X CH , Y CH 3: Let X = X \ X CH , Y = Y \ Y CH 4: if p = 1Calculate the QR decomposition of the sketch X = QR 10:Set Q = XR -1 (see [15, 34] ) 11: Approximate the ℓ p sensitivities by s i := ∥Q i ∥ p p + 1/n 12: Sample k rows of X and Y i.i.d. with probability p i = s i / j s j to obtain X core ∈ R k×d and their corresponding labels Y core ∈ N k 13: Set w core ∈ R k such that w core,j = 1/(kp i ) if sample j corresponds to row i 14: Concatenate X CH with X core to obtain X ′ 15: Concatenate Y CH with Y core to obtain Y ′ 16: Concatenate w CH with w core to obtain w 17: return C = (X ′ , Y ′ , w) Algorithm 2 Domain shift optimizer for pth-root-link Poisson regression.EQUATION). Output: (1 + ε)-approximate solution β. 1: Run Algorithm 1 with input X, Y and k as specified to obtain the coreset (X ′ , Y ′ , w). Let X CH be defined as in Algorithm 1. 2: Run any convex optimization algorithm to find the optimal solution β for the pth-root-linkPoisson regression objective on (X ′ , Y ′ , w) under the constraint that ∀x i ∈ X CH : x i β > ε (see Theorem 4.2) 3: return β.We generated for each p ∈ {1, 2} a dataset with dimensions n = 100 000, d = 7 with n labels corresponding to each point.• Construction of X:1. Start with 6 standard basis vectors (z i ) 6i=1 and add the all zero vector z 0 = 0 to be the extreme points on the convex hull. 2. Construct a matrix in R (n-7)×d with i.i.d. standard Gaussian entries. Translate and rescale the rows of this matrix so that the resulting rows lie in the interior of the convex hull of the points (z i ) 6 i=0 . Concatenate the matrix with the resulting to the matrix with rows given by (z i ) 6 i=0 generated by the first step. Call the resulting matrix Z ∈ R n×d . 3. Horizontally concatenate a column of length n consisting only of ones to the left of the matrix Z (add an intercept) to get X ∈ R n×(d+1) .• Construction of β:1. Draw one sample β ∼ 10 1/p • N (0, I d ), with d = 6 as above. 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Justification: Our research does not involve human subjects or human participants, and addresses all data-related concerns in the NeurIPS Code of Ethics. We do not anticipate any potentially harmful consequences of our research.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [NA] Justification: We consider the research in this paper to be foundational research that is not tied to any particular application or deployment. We do not anticipate any negative application unique to the research presented in this paper.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: The risks posed by the Poisson model analyzed in this paper are no greater than the risks posed by other statistical models used in regression, for example. No datasets were released as part of this paper.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: We only used existing code for which one of the authors has full rights to modify and reuse. The original source was appropriately cited.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 38th Conference on Neural Information Processing Systems (NeurIPS 2024).Õ(•) hides lower order log terms, analogous notation applies to all Landau symbols.The same bound and proof also appeared in[20] in a different context.Such domain restrictions do not lead to feasibility issues, as discussed in Appendix F due to page limitations.Our new code is available at https://github.com/Tim907/poisson-regression/. 
 
 Appendix:  


Decision: This paper studies the problem of constructing coresets for Poisson regression, obtaining sample complexity bounds that depend on a novel and natural complexity parameter of the problem. It builds on a growing line of work on constructing corsets for various “generalized linear models”, including, e.g., recent work on logistic regression. The paper has some clear draw backs, including a lack of experiments, limited new theoretical techniques, and it studies a non-canonical form of Poisson regression. Nevertheless, given the recent interest in coreset problems (e.g., due to their connection to agnostic active learning) we believe the paper will be of interest to many in the NeurIPS community, and all of the reviewers felt it should be accepted to the conference.

Reviews: review: paper_summary: paper_summary: For Poisson regression, where the outcome variable is a positive integer, the paper provides sublinear coresets under certain assumption on the data characterized by parameter $\rho$. The link function is the $p^{th}$ root link function with $p = 1 , p=2$.  Without any assumptions the authors show that no sublinear size coresets are possible, in fact, they show lower bound results for any other data reduction techniques also. For $p = 1 , p=2$ and the parameter $\rho$, the coresets are constructed using the standard sensitivity framework which involves calculating sensitivity upper bounds and bounding the VC-dimension.  For other values of $p$, the authors show that their technique and analysis is not sufficient to get an $\epsilon$-coreset.
strength_weakness: strength_weakness: The only main weakness I think is the lack of any experimental results in the paper. Some small set of proof-of-concept experiments would have further strengthened the paper. 

There are some small issues with writing. $x$ in $R^d$ is described in terms of its coordinates $1, x_1, x_2 \dots x_{d-1}$. However later $x_j$ is used to denote the $j^{th}$ vector in the dataset. 

The writing in the final part of the paper appears rushed, may be due to space constraints. The authors talk about Lambert functions however there is no introduction as to what they are. It is difficult for readers not familiar with the concept.
questions: questions: This may be basic question, maybe I am mission something, the feasible $\beta$ is such that $x_i\beta > 0$ is for the pth-root link or it will also hold in case of ID-link Poisson regression?
limitations: limitations: see weaknesses
paper_summary: paper_summary: This paper demonstrates the theoretical bound analysis on data subsampling for Poisson regression with ID-link and root-link functions based on coreset method and the leverage of $\ell_r$ norms to the loss function. The $\text { ‘ } \rho \text {-complex’ }$ is a novel meaningful parameter for data compressibility under the context of subsampling and is a good starting point of theoretical construction.
strength_weakness: strength_weakness: This paper lacks a well-constructed structure, as it does not clearly introduce important concepts, provide a detailed discussion of related work, or include a summary or conclusion. To improve clarity, the proposed data subsampling techniques for Poisson regression might be presented in a pseudo-code format. While numerical simulation is not necessary, it would be a beneficial addition.
questions: questions: The notation is a bit confusing for the Poisson regression models: is the discussion restricted to $p \in\{(1,2)\}$ cases so the notation in e.g., 119 line, are in scalar format? The notation should be more generalized as in Theorem 3.8.
limitations: limitations: The authors should consider restructuring the entire paper to clearly state the introduction, the proposed method along with its bound analysis, and the limitations and conclusion. Although the problem statement, proof, and results are meaningful, they lack clarity for the readers. Additionally, a major concern is the limited cases to which the analysis can be applied, specifically only for $p \in\{(1,2)\}$. This limitation significantly undermines the generalizability of the proposed method.
paper_summary: paper_summary: Coresets are a technique in efficient algorithms for data analysis in which a dataset is compressed into a weighted subset of its examples. Typically, coresets are constructed such that a loss function for a given optimization problem (e.g. linear regression, logistic regression, or Poisson regression for this work) is preserved up to a $(1\pm\epsilon)$ factor error. This work extends this line of work to handle the Poisson regression case.

While the most popular choice of the log-link does not admit small coresets, the authors study the case of the ID and square root link functions and show that for these log-link functions, Poisson regression admits small coresets, assuming the smallness of a complexity parameter $\rho$ which they introduce building on prior work on logistic regression.
strength_weakness: strength_weakness: The first main weakness of this work is that there are many modifications and assumptions that must be made to the Poisson regression in order to obtain small coresets. As the authors mention, the canonical log-link for Poisson regression does not admit small coresets, and this is a standard fact in the literature (e.g. formalized in [23], as the authors cite). Thus, the authors instead consider a polynomial link function, but the authors do not provide sufficient arguments as to why this is still an interesting problem to study, either from a theory or practice perspective. Indeed, once the growth rate of the function is polynomial and one assumes a balancedness condition, then small coresets can be constructed using known techniques in theory [24]. While this type of result may still be interesting if the regression problem is widely studied, Poisson regression with noncanonical links seems too artificial of a setting to be interesting for a wide audience.

The second main weakness is that the techniques largely are a straightforward generalization of techniques from [23, 24, 25] and the conclusions do not seem surprising given the prior work.
questions: questions: How are Lemmas 6.1 and 6.2 different from the lower bound argument of Theorem 6 in [23]?
limitations: limitations: See weaknesses.
paper_summary: paper_summary: the paper considers sampling for Poisson regression, making explicit the dependence of the size of coreset on various parameters, in particular the effect of different link functions. The paper has a theoretical flavor
strength_weakness: strength_weakness: Frankly I am not familiar with this particular approach. One possible weakness may be that it is not clear how sharp this result is. It is also not clear what characteristics of Poisson regression make the study entirely different from many other regression models. I submit this particular review early so the chair may decide whether to look for additional reviewers.
questions: questions: none
limitations: limitations: seems there is little discussion on limitation

Comments: comment: I have read the authors' response and have no further comments.
Thank you for the additional clarifications, these are very helpful, and the use of the intercept indeed makes everything make sense. I will raise my score to accept.
Regarding feasibility: for the ID and square root link (in fact for any $p$th-root-link), recall that the loss function includes a $\log(x\beta)$ term, which restricts the feasible region to $\beta$ such that for all $x_i, i\in[n]$ it holds that $x_i \beta > 0$. So $\beta\in D(0)$ is the natural domain, induced by the model. This restriction is not our choice.

The domain shift idea restricts the domain even further to $\beta\in D(\eta)\subset D(0)$, for $\eta > 0$. Clearly, some solutions that are feasible in the problem formulated over $D(0)$ are no longer feasible in the problem formulated over $D(\eta)$. But as we prove, we can construct a coreset that holds for all $\beta \in D(\eta)$, and $D(\eta)$ contains at least one $\beta$ that is a $(1+\varepsilon)$-approximation for the optimal solution $\beta^*$ of the problem on the original domain $D(0)$, and evaluated on the full dataset. These two parts are combined to prove that the final minimizer $\tilde\beta\in D(\eta)$ evaluated on the coreset is a $(1+O(\varepsilon))$-approximation for $\beta^*$. 

In the other direction, no infeasible solution can become feasible because of the proper subset relation $D(\eta)\subset D(0)$.

Note that for any data and any fixed $\eta<\infty$, both $D(\eta)$ and $D(0)$ are non-empty, since they consist of all $\beta$ that parameterize hyperplanes that put the convex hull of input points (respectively, the additive $\eta$-inflation of the convex hull of input points) in the positive open halfspace. So there always exist feasible solutions, which means that no instance can become completely infeasible by our methods.

Regarding the question of whether "small coresets [can] still exist ... just not possible for such instances?": We think that small coresets can `still exist' only under assumptions that make the problem of finding a coreset trivial. This is our reasoning: if an instance consists of the extreme points on the convex hull, and all but a small (sublinear) number of points are separated by an $\varepsilon$ distance to the boundary, then indeed the domain shift would not be necessary since the resulting domain structure is already in the data. But if the non-extremal points are allowed to get arbitrarily close to the boundary, and if we do **not** shift the domain, then we will not avoid high sensitivity points that are strictly inside the convex hull. Then the coreset size would necessarily depend on the distance of the non-extremal points to the boundary, and crucially on the number of points that are very close to the boundary of the convex hull, which again can be constructed to be $\Omega(n)$.
Thank you again for your thoughtful comments and questions, and for appreciating our work more than before. Of course we will add the additional references to motivate the model. 

Regarding $\rho$ and $\lambda$: The short answer is yes, the circle lower bound gives $\Omega(n)$ even in the case where $\rho$ is small.

We try to give a more detailed explanation of the role of $\rho$ (and $\lambda$):
Recall in logistic regression [25] the (individual) loss function satisfies $z \leq g(z) \leq 2z$ for sufficiently large $z$. I.e. the upper and lower bounds are balanced up to a factor $2$ independent of $\mu$ (but $\mu$ is used to introduce balance between large (positive) and small (negative) contributions to avoid the circle lower bound).

Now let us take a look at Poisson with ID link (arguments are similar for $p=2$). We would like to bound $z/\lambda \leq g_y(z) \leq z$ for sufficiently large $z$ for some value of $\lambda$. Now $\lambda$ would need to be roughly $y$ (ignoring logarithmic factors) because with growing $y$, the loss function widens, and its minimum moves mainly to the 'right', so we would need a very flat lower bound, so large $\lambda\approx y$. This would unfortunately give $\sum_{i \in[n]} y_i =\Omega(n)$ dependence in coreset size (not reflected in the circle lower bound).

So instead, we additionally shift the lower bound: $(z-y)/\lambda \leq g_y(z) \leq z$ for sufficiently large $z$. This allows sublinear $\lambda$ and thus avoids $\Omega(n)$ (restricted to the large $z$). At the same time we need to relate the lower $(z-y)$ and upper bound $z$, which is the role of $\rho$. This shift and balancing assumption is not artificial or just used to make the calculations go through, but it is naturally consistent with the statistical model (see discussion below Equation (2)). 


Please note that the steps we described in the previous two paragraphs **only** help us to achieve the analogue in the Poisson regression setting of the bounds $z \leq g(z) \leq 2z$ that hold in the logistic regression setting, which is **just** for relating to the $\ell_1$-norm. This is consistent with the fact that the Poisson regression setting is more difficult than the logistic regression setting and indicates that the additional steps we have taken in our paper are nontrivial.

All the above does **not** handle the asymptote near zero, which is the main source of hardness leveraged in the circle $\Omega(n)$ lower bound. The latter difficulty is later tackled by the domain shift idea. 

[To be continued in an additional comment.]
We agree to changing the title to: "Data subsampling for Poisson regression with $p$th-root-link".
Please let us know if you have other suggestions in mind.

We have already chosen $X$ to be the design matrix that includes an intercept. This means that the first column of $X$ consists of only ones. We stated this explicitly in line 119: $x=(1,x^{(1)},x^{(2)},\ldots,x^{(d-1)})$ (notation slightly changed in response to other reviewers). We will add a reminder of that fact in the proof of Lemma 4.1, where the reviewer pointed us.
I have read the rebuttal.  I will keep my score
I have read the response and have no comment
Another comment: perhaps the authors could consider reflecting the fact that only $p$-th root links are considered in the title, since this paper does not give a solution to the canonical Poisson regression problem.

Possibly dumb question: In the proof of Lemma 4.1, why is $\beta'\in D(\eta)$? it seems like if $X\beta \geq 0$ (pointwise), then it does not follow that $X\beta' = X\beta + \eta Xe_1\geq \eta$ depending on the value of $X$ on the first column.
Thank you for the rebuttal. The additional references on the fact that the identity and root link functions are actually used in practice are extremely valuable to me, please consider including it in the draft. I have raised my score based on that.

I still don't quite see your point about the balancedness parameter $\rho$ not being the driving force behind the coreset bound. In the logistic regression setting, having a small $\mu$ parameter allows the coreset to be small, while the circle lower bound instance that is always used does not have a small $\mu$ parameter and thus has a $\Omega(n)$ lower bound on the coreset size [25]. Are you saying that the circle lower bound gives an $\Omega(n)$ lower bound despite a small $\rho$ parameter?

The domain shift idea seems to a technique of avoiding a pathological $\epsilon$-neighborhood of 0 in the loss function to facilitate the analysis, which still contains a $(1+\epsilon)$-approximate minimizer. Indeed, this has not been needed in prior works on sensitivity sampling for GLMs since the loss functions have no asymptotes in prior work. I did not realize this on my first read. Some follow-up questions on feasibility issues: When using the canonical log link, we never have to worry about anything like this since the log link ensures that the expectation of the GLM will be positive. How is this handled for the identity and square root link? And can't the domain restriction turn a feasible instance into an infeasible one? If so, can small coresets still exist and this work's solution avoids them, or are small coresets just not possible for such instances? This should be considered and discussed more carefully. **These are very interesting questions, and I will consider increasing my score to a weak accept depending on how the authors answer these questions. Note that I don't expect all of these questions to be fully answered, as long as a thoughtful discussion is provided.**

Tangential point: I found the discussion of the $\lambda$ parameter to be confusing/not well-explained. I now understand it to be a parameter which parameterizes a lower bound on the loss function.

Re: Lemmas 6.1 and 6.2, I now realize I missed that these lower bounds were for the p-th root versions rather than for the log link version as in [23]. I agree that these lower bounds need to be slightly modified for each coreset problem being studied. The question was not to suggest that including this lower bound is a weakness (although it wouldn't be a stand-alone contribution either), I just wanted to clarify that it couldn't have just been cited from [23].