Optimal Poisson Subsampling Algorithm

• Optimal Poisson subsampling is a design-driven method that independently selects data points with tailored Bernoulli probabilities to minimize asymptotic variance.
• The algorithm employs weighted smoothed quantile generalized estimating equations and a pilot-based approximation to compute optimal inclusion probabilities efficiently.
• Empirical results demonstrate that this approach achieves lower mean squared error and improved computational efficiency compared to uniform subsampling in large-scale studies.

Searching arXiv for the cited paper and closely related optimal Poisson subsampling work to ground the article in current literature. arxiv_search({"query":"(Li et al., 22 Jun 2026) 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 (Li et al., 22 Jun 2026). 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 (Li et al., 28 Aug 2025, Yu et al., 2020, Fan et al., 2022, Wang, 2018).

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 (Li et al., 22 Jun 2026). 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 (Li et al., 22 Jun 2026).

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 (Li et al., 22 Jun 2026). 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 (Li et al., 28 Aug 2025, Yu et al., 2020).

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 (Li et al., 22 Jun 2026). In related frameworks, both A- and L-optimality appear, and more general invariant linear criteria have also been proposed for Poisson designs (Li et al., 28 Aug 2025, Imberg et al., 2023).

2. Poisson inclusion mechanism and weighted estimating equations

The Poisson stage is defined at the subject level. For each subject $i$, let $$Q_\tau(y_{ij} \mid \mathbf{x}_{ij}) = \mathbf{x}_{ij}^\top \boldsymbol{\beta}_0,$$0 denote its inclusion probability and let $$Q_\tau(y_{ij} \mid \mathbf{x}_{ij}) = \mathbf{x}_{ij}^\top \boldsymbol{\beta}_0,$$1 indicate whether the subject enters the subsample (Li et al., 22 Jun 2026). 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 $$Q_\tau(y_{ij} \mid \mathbf{x}_{ij}) = \mathbf{x}_{ij}^\top \boldsymbol{\beta}_0,$$2 (Li et al., 28 Aug 2025).

Because the quantile score contains an indicator and is therefore non-differentiable, the method uses induced smoothing for computational stability. Specifically, $$Q_\tau(y_{ij} \mid \mathbf{x}_{ij}) = \mathbf{x}_{ij}^\top \boldsymbol{\beta}_0,$$3 is replaced by

$$Q_\tau(y_{ij} \mid \mathbf{x}_{ij}) = \mathbf{x}_{ij}^\top \boldsymbol{\beta}_0,$$4

where $$Q_\tau(y_{ij} \mid \mathbf{x}_{ij}) = \mathbf{x}_{ij}^\top \boldsymbol{\beta}_0,$$5 is the standard normal CDF (Li et al., 22 Jun 2026). The resulting weighted smoothed quantile GEE is

$$Q_\tau(y_{ij} \mid \mathbf{x}_{ij}) = \mathbf{x}_{ij}^\top \boldsymbol{\beta}_0,$$6

The estimator $$Q_\tau(y_{ij} \mid \mathbf{x}_{ij}) = \mathbf{x}_{ij}^\top \boldsymbol{\beta}_0,$$7 solves $$Q_\tau(y_{ij} \mid \mathbf{x}_{ij}) = \mathbf{x}_{ij}^\top \boldsymbol{\beta}_0,$$8 on the selected subset (Li et al., 22 Jun 2026).

For practical computation, the paper proposes Newton–Raphson iteration,

$$Q_\tau(y_{ij} \mid \mathbf{x}_{ij}) = \mathbf{x}_{ij}^\top \boldsymbol{\beta}_0,$$9

with termination when the parameter change falls below $$\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),$$0 (Li et al., 22 Jun 2026). 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 (Li et al., 28 Aug 2025, Yu et al., 2020).

3. Optimal inclusion probabilities and the two-step construction

The key design problem is the choice of $$\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),$$1. In the quantile-regression algorithm, the asymptotic variance has the form

$$\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),$$2

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),$$3

and $$\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),$$4 containing the $$\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),$$5 dependence that determines design efficiency (Li et al., 22 Jun 2026). Under the A-optimality criterion, the target is to minimize $$\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),$$6.

The paper defines

$$\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),$$7

and shows that the optimal Poisson subsampling probabilities are

$$\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),$$8

where $$\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),$$9 is the target expected subsample size and $\mathbf{X}_i$0 is a threshold ensuring $\mathbf{X}_i$1 (Li et al., 22 Jun 2026). The structure is informative: larger values of $\mathbf{X}_i$2 correspond to subjects contributing more information about the regression coefficients and therefore receiving larger inclusion probability.

Because $\mathbf{X}_i$3 and $\mathbf{X}_i$4 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 $\mathbf{X}_i$5 to obtain approximate optimal probabilities for the full dataset (Li et al., 22 Jun 2026). To prevent extremely small probabilities, the paper adopts a mixture with uniform sampling,

$\mathbf{X}_i$6

with tuning parameter $\mathbf{X}_i$7, typically $\mathbf{X}_i$8 (Li et al., 22 Jun 2026).

This clipped-probability architecture is not unique to quantile regression. In general M-estimation, the optimal Poisson probabilities also take a thresholded form involving $\mathbf{X}_i$9, and the same structural device appears in generalized estimating equations and quasi-likelihood estimation (Wang et al., 2022, Li et al., 28 Aug 2025, Yu et al., 2020). A plausible implication is that clipping is not merely a technical artifact of one model class, but a recurrent solution to the feasibility constraint $(m \times p)$0 in optimal unequal-probability Poisson designs.

4. Large-sample theory and penalized estimation

The subsample estimator in the quantile-regression framework is $(m \times p)$1-consistent for $(m \times p)$2 (Li et al., 22 Jun 2026). Its asymptotic normality is stated as

$(m \times p)$3

under standard, mild regularity conditions denoted $(m \times p)$4–$(m \times p)$5 (Li et al., 22 Jun 2026). 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

$(m \times p)$6

where

$(m \times p)$7

with weights derived from the initial unpenalized subsampled estimator (Li et al., 22 Jun 2026). Minimization is carried out by a minorization-maximization Newton strategy, which modifies the gradient and Hessian to accommodate the non-smooth penalty (Li et al., 22 Jun 2026).

The asymptotic results for the penalized version include an oracle property: if $(m \times p)$8 and $(m \times p)$9, the procedure identifies the correct nonzero coefficients with probability tending to $i$0, and the nonzero coefficients retain asymptotic normality at the optimal rate (Li et al., 22 Jun 2026). 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 (Li et al., 22 Jun 2026). 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 (Li et al., 22 Jun 2026). 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 (Li et al., 22 Jun 2026).

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 (Wang et al., 2022). 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 (Fan et al., 2022). 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 (Li et al., 28 Aug 2025).

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 (Yu et al., 2020). 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 (Imberg et al., 2023).

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 $i$1, and practical large-scale implementations therefore use truncated Poisson batch samplers with capping and padding to satisfy fixed-shape hardware constraints (Chua et al., 2024). 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 (Dong et al., 8 May 2026).

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 (Comino-Trinidad et al., 2023, Flamich, 2023). 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.