Data subsampling for Poisson regression with pth-root-link

In the last line of our abstract (lines 18-20) we state that we show the limitations of our analysis for $p$-th degree root link functions for $p\geq 3$, and state that these limitations show the need for other methods if one aims to generalize our approach to this range of values of $p$. We repeat this point in lines 62-65. We also indicated on the NeurIPS Paper Checklist that Section 6 contains lower bounds that indicate the limitations of constructing coresets for $p$-th power link Poisson models in general, and the limitations of our specific approach and analysis, for instance the limitations regarding input/complexity parameters such as $y_{\max}$. We also point to [23] for limitations regarding the canonical log link.

Our lower bounds on the parameters, together with linear VC dimension and linear sensitivity, and linear $\rho$ dependence leave no room for improvement if one uses the sensitivity framework. We also refer to the related work Section 1.3, where previous finite sample size results were either unbounded or $O(\sqrt{d})$ error, instead of our $1+\varepsilon$. We will add some comments on this.

We do not agree with the statement that Poisson regression with links other than the canonical log link is 'artificial'. The identity link was applied in epidemiology, see e.g. [a], or [b]. The root link function has been applied to forecasting for queueing systems [c], and to account for misspecification bias in maximum likelihood estimation [d]. 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. the paragraph "Use of the square root (sqrt) link function" in Section 5.4 of [e]. There are also informal discussions on stackexchange that shed more light on the choice of link function, which we unfortunately are not allowed to link in our response.

While the general outline follows the same structure as previous papers (in fact almost all papers based on the sensitivity framework share this outline), the details are quite different. First note that the change of link function changes the norms to relate to from $\ell_\infty$-norm to $\ell_1/\ell_2$. The reviewer seems to suggest that this choice of polynomial functions resolves the problem given previous literature, which is indeed not the case as our linear lower bounds prove. Indeed, the fact that even after chosing polynomial functions, the loss still depends on the convex hull (i.e. $\infty$-norm) is somewhat challenging in our analysis. Next, the 'balancedness' $\mu$ assumption in [24,25] is needed due to the asymmetry of logistic/probit losses. However, in our paper, the asymmetry of Poisson loss is taken care of by our novel domain shift idea to avoid large contributions. It is not taken care of by the $\rho$ parameter! Also, $\rho$ has a different interpretation that is natural for Poisson, but not applicable to logit nor probit. It balances the mean and the variance, and motivates 'moving' the polynomial lower envelopes. Without this, we stress that the slope $\lambda$ parameter would need to be $y_{\max}/\log(y_{\max})$ since with growing $y$ the loss function 'moves' polynomially along the $x\beta$ axis but only logarithmically along the $g_y(x\beta)$ axis. See our high level discussion in lines 83-89. Indeed, though not obvious from the current writeup, this induces a $\sum_{i \in[n]} y_i =\Omega(n)$ dependence. So it is a completely different difficulty that the supposedly 'not surprising' or 'straightforward' analysis addresses by the balancing condition. The most common point between $\rho$ and $\mu$ is the fact that both correspond to or arise naturally from their respective statistical modeling. Please also consider the further discussion and the separation between the $p\in\\{1,2\\}$ cases with respect to the $y_{\max}$ dependence. These are further points that were previously completely unexplored and led to novel and interesting non-trivial results regarding the important Lambert $W_0$ function as a byproduct. So while some basic parts were borrowed from previous literature that we cited, we kindly ask the reviewer to recognize as well our non-trivial extensions and modifications that go beyond the polynomial growth and $\ell_p$ space approximation of [24,25] and other previous literature.

Lemmas 6.1 and 6.2 are not significantly different from the lower bound argument of Theorem 6 in [23]. Indeed, we had already stated in lines 96-97 that these results are adapted from previous literature, admittedly without specifying Lemmas 6.1 and 6.2 exactly. This was an oversight on our part and we will clarify the connection to Theorem 6 in [23]. While the construction of the bad dataset is similar to [23] and is widely standard across previous literature as in e.g. [25,27], [f], [g], each of these references as well as ours adapt the construction to the peculiarities of their loss function under study. In particular, [23] puts the hyperplane into the convex hull to detect the query point being 'far' on the expensive exponential side of the loss (other points are on the cheaper linear side). In our case, we put the hyperplane outside the convex hull, to detect the query point being 'close' to the hyperplane, so that its negative logarithmic loss dominates the polynomial loss of all other points. While this can in hindsight be considered a minor modification, with the same reasoning, we should reject a lot of previous NeurIPS, ICML, AAAI publications for building on the IJCAI paper [g] (and maybe even older references).

[a] I. C. Marschner (2010). "Stable Computation of Maximum Likelihood Estimates in Identity Link Poisson Regression", Journal of Computational and Graphical Statistics, 19(3), 666–683. DOI: 10.1198/jcgs.2010.09127

[b] Donna Spiegelman, Ellen Hertzmark (2005). "Easy SAS Calculations for Risk or Prevalence Ratios and Differences", American Journal of Epidemiology, Volume 162, Issue 3, Pages 199–200. DOI: 10.1093/aje/kwi188

[c] Haipeng Shen. Jianhua Z. Huang (2008). "Forecasting time series of inhomogeneous Poisson processes with application to call center workforce management", Ann. Appl. Stat. 2 (2) 601 - 623. DOI: 10.1214/08-AOAS164

[d] B. Efron (1992). "Poisson Overdispersion Estimates Based on the Method of Asymmetric Maximum Likelihood", Journal of the American Statistical Association, 87(417), 98–107. DOI: 10.1080/01621459.1992.10475180

[e] J. Maindonald, W. J. Braun (2010). "Data analysis and graphics using R: an example-based approach" (4th ed.). Cambridge University Press. DOI: 10.1017/CBO9781139194648

[f] J. Huggins, T. Campbell, T. Broderick (2016). "Coresets for scalable Bayesian logistic regression". Advances in neural information processing systems, 29. NeurIPS 2016.

[g] Sariel Har-Peled, Dan Roth, Dav Zimak (2007). "Maximum Margin Coresets for Active and Noise Tolerant Learning". 836-841, Proceedings of the 20th International Joint Conference on Artificial Intelligence. IJCAI 2007.

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.

1. We have difficulty understanding what exactly is lacking in the structure of the paper. We remind the reviewer of the structure of our paper: In Section 1, we provide an introduction that introduces the crucial concepts of a $(1+\varepsilon)$-coreset and sublinear bounds on the coreset size. In Section 1.1, we list our contributions. In Section 1.2, we describe the high level ideas of our techniques --- for almost one page --- and introduce the important ideas that underpin our results. In Section 1.3, we give a representative sample of related work. In Section 2, we introduce more important concepts, e.g. the complexity parameter in (2) and the loss functions in (3) and (4) that are essential for the rest of the paper. In Section 3, we again introduce important concepts relating to VC dimension and the loss function from Section 2. After more than six pages of structured introduction of concepts, our first novel result appears on page 7 in Theorem 3.8, which gives an upper bound on the coreset size. In Section 4, we state our main approximation result, and introduce a novel shifting idea to overcome a difficulty posed by our loss function. There still remain difficulties for the subsequent optimization. We address these difficulties in Section 5, using ideas based on previous literature. In Section 6, we give further context to the upper bound on the coreset size from Theorem 3.8. We do this by providing lower bounds on the coreset size.

The above shows that we introduced important concepts and presented our main results in a logical progression. Note also that we already summarized the contents of our paper in the abstract. We do not agree that a conclusion is absolutely necessary for a "well-constructed structure", but we would be willing to include a conclusion on a tenth page. For our initial submission, we had already moved as much content as possible to the supplementary and kept only the important results of our work in the main paper. Note that no other reviewer had issues with the structure of our paper.

2. The goal of a literature review in an 8-page paper is not to be detailed or exhaustive. Instead, the goal is to point the reader to representative examples of relevant literature, and to explain the differences between these examples and our present work, and this is exactly what we have done on lines 109-117. If the reviewer can point to specific examples of literature that they feel should have been included in the review, then we invite the reviewer to state these examples and to indicate why these suggestions are not adequately represented by the references we cite in lines 109-117.

3. We are willing to provide pseudocode in the supplementary material.

4. We agree that numerical results that provide proof of concept are not necessary but could be beneficial. We now have results of this kind and are willing to report them in the supplementary material. Please refer to our response to reviewer SQxM for more details regarding numerical results.

We suspect that the confusion is due to an issue with our notation: on line 119, "$x=(1,x_1,\ldots,x_{d-1})\in\mathbb{R}^{d}$" indicates a row vector, whereas in Theorem 3.8, $x_i$ on line 277 indicates a row vector. We will change the notation of coordinates of $x\in\mathbb{R}^d$ to $x=(1,x^{(1)},x^{(2)},\ldots,x^{(d-1)})$ on line 119 and consistently write $x_i$ to denote the $i$-th row vector in the data matrix $X\in\mathbb{R}^{n\times d}$ as in Theorem 3.8.

The reviewer describes the problem statement, proof, and results as lacking clarity, and also as being meaningful. It is difficult for us to understand how both these statements can be true, since a necessary condition for a statement to be meaningful is that it is clear. We invite the reviewer to give one or two specific examples of results from our paper and to explain exactly what they do not understand about those results. Without this information, the reviewer's feedback appears vague and generic.

We point out that the other reviewers who are familiar with the area of coresets had no issues with the clarity of our problem statements, proofs, or results.

Our paper is about constructing coresets of suitable size for Poisson regression, and we achieve this goal for some link functions that satisfy two criteria: 1) they can be treated using a combination of known coreset construction methods and our novel techniques and results; and 2) they have been and continue to be used in statistics (see response to reviewer HZmN). At no point do we claim generalizability of the proposed method to all link functions or all $p$. For this reason, the reviewer's comment appears to be finding issue with aspects of the paper that are not relevant to our stated goals. On the contrary, our lower bound results (together with [23]) indicate the restriction to $p\in\\{1,2\\}$, and exclude $p\geq 3$, in particular the limiting case $p=\infty$ that corresponds to the canonical log-link. According to the reviewing guidelines for NeurIPS: "authors should be rewarded rather than punished for being up front about the limitations of their work" (see lines 18-20, and 62-65 for early mentions of the limitations). Moreover $p\in\\{1,2\\}$ are important standard choices in statistics (see response to reviewer HZmN).

We will describe some proof-of-concept results that have been obtained since we submitted the paper for review. We have generated 6-dimensional data that consists of the vertices of a regular simplex (as extreme points on the convex hull) and $n$ further points from a normal distribution rescaled to be in the convex hull. Our resulting approximation factors are close to $1$ even for only $20-50$ subsamples except for a few repetitions where it fails, and gives factors $2-7$. In contrast, even for $500$ samples, uniform sampling fails even to produce feasible $\beta$ in one third of the cases, and even when results are feasible and the extreme points are included, some approximation ratios are in the order of $1.5-2.5 \times 10^9$ which are clearly explained by missing points that are very close to the boundary of the convex hull, thus causing huge errors. We will elaborate more on this and add corresponding numerical proof-of-concept in the supplementary material.

We agree that the notation can be improved and will change the notation of coordinates of $x\in\mathbb{R}^d$ to $x=(1,x^{(1)},x^{(2)},\ldots,x^{(d-1)})$. We will continue to write $x_j$ to denote the $j$-th vector in the dataset.

If by 'final part of the paper' the reviewer was referring to Section 6, then it is true that we had to remove some text from this section due to space constraints. In particular, we will give the definition of the Lambert function in the text and give some intuition as to why it is necessary to consider this function. Some intuition is already included in the last paragraph of Section 1.2.

For ID-link Poisson regression, the set of feasible $\beta$ also consists of $\beta$ such that $x_i\beta >0$. This is because the loss function involves $\log (x_i\beta)$ for the ID link as well.