---
title: Optimal Poisson Subsampling Algorithm
url: https://www.emergentmind.com/topics/optimal-poisson-subsampling-algorithm
type: topic
---

# Optimal Poisson Subsampling Algorithm

Searching arXiv for the cited paper and closely related optimal Poisson subsampling work to ground the article in current literature.
arxiv_search({"query":"2606.23363 optimal Poisson subsampling quantile regression large-scale longitudinal data", "max_results": 5})
An optimal Poisson subsampling algorithm is a subsample-construction and reweighted-estimation procedure for massive data in which each unit is independently retained with a designed Bernoulli probability, typically chosen to minimize an asymptotic mean squared error criterion under a fixed expected subsample size. In the large-scale longitudinal quantile-regression formulation, the algorithm combines subject-level Poisson inclusion, weighted smoothed quantile generalized estimating equations, induced smoothing, and a pilot-based approximation to A-optimal inclusion probabilities, thereby reducing computational burden while preserving consistency and asymptotic normality [2606.23363]. Related work has developed closely related Poisson-optimal designs for generalized estimating equations, quasi-likelihood estimation, M-estimation, and logistic regression, indicating that the method belongs to a broader family of unequal-probability subsampling schemes for statistically efficient inference under severe computational constraints [2508.20803], [2005.10435], [2209.04569], [1802.02698].

## 1. Statistical setting and target of inference

In the longitudinal quantile-regression setting, the model is specified at quantile level $\tau$ by
\[
Q_\tau(y_{ij} \mid \mathbf{x}_{ij}) = \mathbf{x}_{ij}^\top \boldsymbol{\beta}_0,
\]
where repeated measurements within subject induce within-subject correlation [2606.23363]. The estimating framework is quantile generalized estimating equations, with
\[
\mathbf{S}_n(\boldsymbol{\beta}) = \frac{1}{n} \sum_{i=1}^n \mathbf{X}_i^\top \boldsymbol{\Gamma}_i \mathbf{R}_i^{-1}\Big( \tau - I\{\mathbf{Y}_i \leq \mathbf{X}_i \boldsymbol{\beta}\} \Big),
\]
where $\mathbf{X}_i$ is the $(m \times p)$ covariate matrix for subject $i$, $\mathbf{Y}_i$ is the response vector, $\mathbf{R}_i$ is a working correlation matrix, and $\boldsymbol{\Gamma}_i$ is a diagonal matrix of conditional densities at zero [2606.23363].

The computational bottleneck is that direct evaluation of the full estimating equation becomes expensive when the number of subjects is large. The optimal Poisson subsampling algorithm addresses this by selecting a smaller, informative subset of subjects while preserving the full-sample target through weighting [2606.23363]. The same objective appears in other large-data settings: the proposed method in generalized estimating equations with diverging covariate dimension is explicitly designed to reduce storage and computation costs while retaining statistical accuracy, and the quasi-likelihood literature likewise treats Poisson subsampling as a device for replacing infeasible full-data optimization by weighted subsample estimation [2508.20803], [2005.10435].

A central feature of this literature is that “optimality” is criterion-dependent. In the quantile-regression formulation, the criterion is A-optimality, expressed through minimization of the trace of the asymptotic variance matrix of the subsample estimator [2606.23363]. In related frameworks, both A- and L-optimality appear, and more general invariant linear criteria have also been proposed for Poisson designs [2508.20803], [2304.03019].

## 2. Poisson inclusion mechanism and weighted estimating equations

The Poisson stage is defined at the subject level. For each subject $i$, let $\pi_i$ denote its inclusion probability and let $\delta_i \sim \mathrm{Bernoulli}(\pi_i)$ indicate whether the subject enters the subsample [2606.23363]. This independence structure is the defining feature of Poisson subsampling in the statistical literature: each unit is sampled independently, and the realized sample size is random with expectation controlled by $\sum_i \pi_i$ [2508.20803].

Because the quantile score contains an indicator and is therefore non-differentiable, the method uses induced smoothing for computational stability. Specifically, $I\{\mathbf{Y}_i \le \mathbf{X}_i \boldsymbol{\beta}\}$ is replaced by
\[
\Phi\Big( \frac{\mathbf{Y}_i - \mathbf{X}_i \boldsymbol{\beta}}{h} \Big),
\]
where $\Phi(\cdot)$ is the standard normal CDF [2606.23363]. The resulting weighted smoothed quantile GEE is
\[
\mathbf{S}_r^\Phi(\boldsymbol{\beta}) =
\frac{1}{n} \sum_{i=1}^n
\frac{\delta_i}{\pi_i}
\mathbf{X}_i^\top \boldsymbol{\Gamma}_i \mathbf{R}_i^{-1}
\left(
\Phi\Big( \frac{\mathbf{Y}_i - \mathbf{X}_i \boldsymbol{\beta}}{h} \Big)
-
(1-\tau)
\right).
\]
The estimator $\tilde{\boldsymbol{\beta}}$ solves $\mathbf{S}_r^\Phi(\boldsymbol{\beta}) = 0$ on the selected subset [2606.23363].

For practical computation, the paper proposes Newton–Raphson iteration,
\[
\boldsymbol{\beta}^{(t+1)} = \boldsymbol{\beta}^{(t)} -
\left[\bar{\mathbf{H}}_n\left(\boldsymbol{\beta}^{(t)}\right)\right]^{-1}
\bar{\mathbf{S}}_r^\Phi\left(\boldsymbol{\beta}^{(t)}\right),
\]
with termination when the parameter change falls below $10^{-4}$ [2606.23363]. This emphasis on iterative weighted estimation is shared by other Poisson-optimal procedures. In generalized estimating equations, the subsample-based weighted GEE is solved by iterative weighted least squares, and in quasi-likelihood estimation the weighted quasi-likelihood equation is solved on the Poisson subsample using inverse-probability correction [2508.20803], [2005.10435].

## 3. Optimal inclusion probabilities and the two-step construction

The key design problem is the choice of $\pi_i$. In the quantile-regression algorithm, the asymptotic variance has the form
\[
\mathbf{V} = \mathbf{H}_n^{-1} \mathbf{V}_c \mathbf{H}_n^{-1},
\]
with
\[
\mathbf{H}_n = -\frac{1}{n}\sum_{i=1}^n \mathbf{X}_i^\top \boldsymbol{\Gamma}_i \mathbf{R}_i^{-1} \boldsymbol{\Gamma}_i \mathbf{X}_i,
\]
and $\mathbf{V}_c$ containing the $1/\pi_i$ dependence that determines design efficiency [2606.23363]. Under the A-optimality criterion, the target is to minimize $\operatorname{tr}(\mathbf{V})$.

The paper defines
\[
z_i = \left\| \mathbf{H}_n^{-1} \mathbf{X}_i^\top \boldsymbol{\Gamma}_i \mathbf{R}_i^{-1} \right\|,
\]
and shows that the optimal Poisson subsampling probabilities are
\[
\pi_i = r \frac{z_i \wedge T}{\sum_{j=1}^n (z_j \wedge T)},
\]
where $r$ is the target expected subsample size and $T$ is a threshold ensuring $\pi_i \leq 1$ [2606.23363]. The structure is informative: larger values of $z_i$ correspond to subjects contributing more information about the regression coefficients and therefore receiving larger inclusion probability.

Because $\mathbf{H}_n$ and $\mathbf{R}_i$ are unknown, implementation proceeds through a two-step approximation. A small initial uniform subsample is used to estimate the regression parameter and the working correlation structure; these pilot estimates are then plugged into the formula for $z_i$ to obtain approximate optimal probabilities for the full dataset [2606.23363]. To prevent extremely small probabilities, the paper adopts a mixture with uniform sampling,
\[
\tilde{\pi}_i = (1-\rho)\pi_i^{\mathrm{opt}} + \rho \frac{1}{n},
\]
with tuning parameter $0 < \rho < 1$, typically $\rho = 0.2$ [2606.23363].

This clipped-probability architecture is not unique to quantile regression. In general M-estimation, the optimal Poisson probabilities also take a thresholded form involving $\|\dot m(Z_i,\hat{\boldsymbol{\theta}}_n)\| \wedge H$, and the same structural device appears in generalized estimating equations and quasi-likelihood estimation [2205.08588], [2508.20803], [2005.10435]. A plausible implication is that clipping is not merely a technical artifact of one model class, but a recurrent solution to the feasibility constraint $\pi_i \leq 1$ in optimal unequal-probability Poisson designs.

## 4. Large-sample theory and penalized estimation

The subsample estimator in the quantile-regression framework is $r^{-1/2}$-consistent for $\boldsymbol{\beta}_0$ [2606.23363]. Its asymptotic normality is stated as
\[
\mathbf{V}^{-1/2} \left( \tilde{\boldsymbol{\beta}} - \boldsymbol{\beta}_0 \right)
\xrightarrow{d} N(\mathbf{0}, \mathbf{I}),
\]
under standard, mild regularity conditions denoted $(\mathrm{C}1)$–$(\mathrm{C}7)$ [2606.23363]. The paper also states that the asymptotic variance is minimized by the proposed optimal sampling probabilities, so the design criterion is directly linked to the limiting law of the estimator.

The framework is extended to penalized estimation through adaptive LASSO. The penalized estimating equation is
\[
\mathbf{S}_r^P(\boldsymbol{\beta}) =
\mathbf{S}_r^\Phi(\boldsymbol{\beta}) -
\lambda \tilde{\mathbf{w}} \circ \mathrm{sgn}(\boldsymbol{\beta}),
\]
where
\[
\tilde{\mathbf{w}} = (|\tilde{\beta}_1|^{-\gamma}, \ldots, |\tilde{\beta}_p|^{-\gamma}),
\]
with weights derived from the initial unpenalized subsampled estimator [2606.23363]. Minimization is carried out by a minorization-maximization Newton strategy, which modifies the gradient and Hessian to accommodate the non-smooth penalty [2606.23363].

The asymptotic results for the penalized version include an oracle property: if $\sqrt{r}\lambda \rightarrow 0$ and $r\lambda \rightarrow \infty$, the procedure identifies the correct nonzero coefficients with probability tending to $1$, and the nonzero coefficients retain asymptotic normality at the optimal rate [2606.23363]. This places the optimal Poisson subsampling algorithm within the class of high-dimensional inference procedures that combine design-based efficiency with support recovery.

## 5. Efficiency comparisons and empirical behavior

The empirical evidence reported for longitudinal quantile regression is consistent across simulated and real data. Relative to uniform Poisson subsampling, the optimal design yields uniformly lower mean squared error, with biases that are similar or smaller and standard deviations that are consistently lower [2606.23363]. The efficiency gains persist under working-correlation misspecification, while computation remains overwhelmingly faster than fitting the model on the full sample, even though optimal subsampling is slightly slower than uniform subsampling because probabilities must be calculated [2606.23363]. In the CHARLS application for depression scores, the optimal subsample and the full sample yield similar variable-selection results, with substantially improved computational efficiency for the subsampling approach [2606.23363].

The broader comparative literature supplies an important asymptotic interpretation. In a general target-function framework, Poisson subsampling and subsampling with replacement have the same variance structure when the subsampling ratio tends to zero, but for larger subsampling fractions the Poisson estimator has strictly smaller asymptotic variance [2205.08588]. This explains why Poisson designs are repeatedly recommended when the subsample is not negligible relative to the full dataset.

Efficiency gains also arise from the weighting stage, not only from the sampling stage. In capture–recapture M-estimation, an initial uniform sample followed by a second Poisson sample can be combined with empirical likelihood weighting rather than inverse probability weighting. The resulting ELW estimator is always no greater than the IPW estimator in asymptotic variance, can incorporate auxiliary information, and leads to more efficient sampling plans and more economical sample sizes for prespecified precision [2209.04569]. In generalized estimating equations with growing dimension, both A-optimal and L-optimal Poisson probabilities substantially outperform uniform Poisson sampling and remain effective even when the working correlation matrix is misspecified [2508.20803].

Taken together, these results support a precise interpretation of “optimal Poisson subsampling”: the method is not merely a faster approximation to full-data fitting, but a design-and-estimation scheme in which the sampling law and the estimating equation are jointly tuned to reduce the asymptotic covariance of the final estimator.

## 6. Generalizations, distinctions, and scope of the term

The optimal Poisson subsampling algorithm is best understood as a family of criterion-driven Bernoulli inclusion schemes rather than a single invariant formula. In distributed quasi-likelihood estimation, optimal Poisson probabilities are derived under A- and L-optimality and can be computed block-by-block, which removes the need to calculate all probabilities at once and makes the method suitable for data stored across multiple locations [2005.10435]. In the general theory of optimal subsampling design, Poisson and multinomial schemes are treated under a unified covariance calculus, and new invariant linear criteria are proposed to achieve near D-optimal efficiency at much lower computational demand than D-optimality itself [2304.03019].

The term also requires separation from other “Poisson subsampling” usages. In differentially private stochastic gradient descent, Poisson subsampling means that each example is independently included in a minibatch with probability $b/n$, and practical large-scale implementations therefore use truncated Poisson batch samplers with capping and padding to satisfy fixed-shape hardware constraints [2411.04205]. That setting optimizes privacy amplification rather than asymptotic estimation efficiency. Indeed, Balanced Iteration Subsampling was shown to achieve stronger privacy amplification than Poisson subsampling and to be optimal at both extremes of the noise spectrum, so “optimal” in privacy-sensitive training does not identify Poisson subsampling as the best scheme [2605.07072].

The label should also be distinguished from Poisson-disk subsampling for point-cloud decimation, where the objective is geometric uniformity rather than statistical inference, and from greedy Poisson rejection sampling, where a Poisson process is used for one-shot channel simulation and coding rather than for selecting a weighted statistical subsample [2311.17604], [2305.15313]. Within statistical learning and inference, however, the dominant meaning remains clear: an optimal Poisson subsampling algorithm is an independently randomized, unequal-probability subsampling design coupled with a weighted estimator, constructed so that a prescribed covariance-based criterion is minimized subject to a fixed computational budget.

Source: https://www.emergentmind.com/topics/optimal-poisson-subsampling-algorithm