PaperlayerView exact revision on arXiv

Source-linked AI summary

Optimal Subsampling Algorithms for Big Data Regressions

Mingyao Ai, Jun Yu, Huiming Zhang, HaiYing Wang

arXiv:1806.06761v2stat.MEmath.STstat.CO

TL;DR

Massive-data GLM fitting remains computationally difficult because iterative Newton-Raphson optimization requires O(p2n) time per iteration. This paper develops optimal subsampling strategies based on A- and L-optimality, establishes their theoretical properties, and evaluates an adaptive two-step approach that outperforms uniform subsampling in reported simulations.

  • Problem

    Fitting generalized linear models on massive datasets is computationally difficult because Newton-Raphson optimization requires O(p2n) time per iteration.

  • Method

    The paper derives mV- and mVc-optimal subsampling probabilities under A- and L-optimality criteria, establishes asymptotic properties, and develops finite-sample analyses for the resulting estimators.

  • Results

    The two-step approach outperforms uniform subsampling across all four cases for both Poisson and negative binomial regression when r0/r ∈[0.1, 0.9].

  • Takeaways & Limitations

    Optimal subsampling provides a theoretically analyzed way to reduce data volume while approximating full-data GLM estimation, with reported confidence-interval performance uniformly better than uniform subsampling.

  • Takeaways & Limitations

    The method relies on moment and covariate conditions, including assumptions needed for consistency and the Lindeberg-Feller central limit theorem; extremely small sampling probabilities can also inflate estimator variance.

Abstract

from arXiv · show

To fast approximate maximum likelihood estimators with massive data, this paper studies the Optimal Subsampling Method under the A-optimality Criterion (OSMAC) for generalized linear models. The consistency and asymptotic normality of the estimator from a general subsampling algorithm are established, and optimal subsampling probabilities under the A- and L-optimality criteria are derived. Furthermore, using Frobenius norm matrix concentration inequalities, finite sample properties of the subsample estimator based on optimal subsampling probabilities are also derived. Since the optimal subsampling probabilities depend on the full data estimate, an adaptive two-step algorithm is developed. Asymptotic normality and optimality of the estimator from this adaptive algorithm are established. The proposed methods are illustrated and evaluated through numerical experiments on simulated and real datasets.

1. Introduction

Massive datasets make full-data GLM fitting computationally difficult, motivating subsampling strategies that reduce data volume while preserving estimation quality. The paper develops optimal and adaptive subsampling methods for GLMs.

  • Fitting GLMs on massive data requires O(p2n) time per Newton-Raphson iteration.
  • Subsampling reduces data volume and has been developed into leverage-based and optimal procedures for regression estimation.
  • Existing statistically guaranteed subsampling results remain limited for complex models, with most work focused on linear, logistic, or multinomial regression.
  • The paper derives general and optimal GLM subsampling methods, finite-sample properties, and an adaptive two-step algorithm evaluated on simulations and real data.

2. Preliminaries

The paper formulates GLMs through exponential-family conditional distributions and studies weighted estimators formed from with-replacement subsamples. Under regularity conditions, the subsample estimator consistently approximates the full-data MLE and has asymptotically normal approximation error.

  • 2.1 Models and Assumptions: GLMs model each conditional response distribution through an exponential-family form with natural parameter u(β^Tx_i).
  • 2.2 General Subsampling Algorithm and its Asymptotic Properties: The general algorithm samples r observations with replacement according to probabilities π_i and computes a weighted log-likelihood estimator.
  • 2.2 General Subsampling Algorithm and its Asymptotic Properties: Inverse-probability weighting makes the subsample estimator a weighted MLE and avoids the general bias of unweighted estimation when probabilities depend on responses.
  • 2.2 General Subsampling Algorithm and its Asymptotic Properties: Under assumptions including compactness, positive-definiteness, finite covariate moments, and sampling-probability conditions, the estimator is consistent for the full-data MLE.
  • 2.2 General Subsampling Algorithm and its Asymptotic Properties: The estimator converges to the full-data MLE at rate r^-1/2, and its approximation error is asymptotically normal in conditional distribution.

3. Optimal Subsampling Strategies

The paper derives A- and L-optimal subsampling probabilities for GLMs, balancing estimation efficiency against computational cost. It also develops finite-sample accuracy results that relate approximation error to subsample size, dimension, covariate structure, and model fit.

  • A-optimality: A-optimal subsampling minimizes the asymptotic mean squared error of the subsample estimator.The resulting strategy is called mV-optimal.
  • Information matrices: Under mild conditions, the empirical information matrix can replace JX when computing mV-optimal probabilities because both converge to the same Fisher information limit.This replacement preserves the stated asymptotic result.
  • A-optimality: mV-optimal probabilities favor observations with larger directional derivatives because they provide greater A-optimal efficiency gains.The directional derivative measures the relative gain from adding an observation to the sample.
  • L-optimality: The mVc-optimal criterion reduces probability-computation cost from O(np^2) to O(np) by minimizing an upper bound through tr(Vc).The alternative criterion preserves Loewner-order improvements in the estimator covariance and is described as losing little estimation accuracy.
  • Finite-sample properties: Finite-sample excess-risk bounds quantify prediction accuracy for optimal subsampling at fixed n and r, with accuracy improving as r increases.The bounds expose effects from covariate structure, dimension, and goodness of fit.
  • Finite-sample properties: The error analysis suggests choosing r proportional to p, while also accounting for the condition number of the weighted design matrix.The condition number reflects covariate collinearity, and p captures the curse of dimensionality.

4. Practical Consideration and Implementation

The practical algorithm replaces unknown full-data quantities with a pilot estimate, thresholds residual-based probabilities for robustness, and then samples using approximated optimal probabilities. Under stated asymptotic conditions, the resulting estimator has established consistency-related bounds and asymptotic normality with the optimal variance form.

  • Practical implementation: The optimal probabilities depend on the unknown full-data MLE, so a pilot estimator is used to approximate them.The pilot estimator is based on a subsample of size r0 and substitutes for β_MLE in the probability formulas.
  • Practical implementation: Approximated optimal probabilities are proportional to absolute fitted residuals, assigning very small probabilities to observations with near-zero residuals.Observations whose responses exactly match fitted means may receive zero probability and never enter the subsample.
  • Practical implementation: A positive threshold δ replaces residual magnitudes in the probabilities, truncating weights to protect the weighted likelihood from domination by low-probability observations.The paper gives δ = 10^-6 as an example and notes that sufficiently small δ preserves efficiency approximately.
  • Two-step algorithm: The two-step algorithm first draws a uniform pilot subsample, then estimates probabilities, samples with replacement, and maximizes a weighted log-likelihood.The final subsample combines the pilot set with the observations drawn under the approximated optimal probabilities.
  • Theoretical guarantees: Under the stated assumptions and r0/r → 0, the two-step estimator has asymptotic properties including asymptotic normality with variance Vopt.The results are established as r0, r, and n diverge, conditional on the full data.

5. Numerical Studies

Numerical studies evaluate the proposed two-step subsampling methods for Poisson and negative binomial regression using eMSE, confidence intervals, allocation choices, tuning sensitivity, and computation time. The optimal-probability methods generally outperform uniform and leverage-based sampling, while performance depends on correlation, condition number, sample allocation, and model.

  • Simulation setup: The simulations use Poisson regression with n = 10,000, seven covariates, and four covariate-generation cases, including highly correlated covariates.Cases vary covariate independence, correlation, and marginal distributions.
  • Methods compared: The compared strategies include uniform sampling, leverage sampling, adjusted leverage sampling, and the proposed πmV and πmVc probabilities.The two-step implementation fixes r0 = 200 while varying the second-step size r.
  • Poisson regression: Across all four datasets, πmV and πmVc produce smaller eMSEs than uniform sampling and have similar performance when covariates are independent.Their performance can differ when covariates are highly correlated because πmVc reduces the impact of the data correlation structure.
  • Poisson regression: Condition number affects required sample size: Cases 1, 3, and 4 perform well with r = 100, whereas Case 2 has condition number ≈40 and needs a larger subsample.The other cases have condition number approximately 5.
  • Inference: The proposed methods provide similar and uniformly better confidence-interval performance than uniform sampling, with interval lengths decreasing as r increases.Intervals are longer in Case 2, supporting larger subsamples when the full-data condition number is higher.
  • NBR and tuning: For NBR with ν = 2, eMSEs are larger than in Poisson regression at the same r, while the two-step method outperforms uniform sampling for r0/r ∈ [0.1, 0.9].The allocation performs best when r0/r is around 0.2, and eMSE is not sensitive to δ when δ is not big.
  • Computational efficiency: All subsampling algorithms require significantly less computing time than the full-data approach, with πmVc using O(pn) time to compute probabilities.πmVc is faster than leverage-based methods and its computational advantage increases with dimension.

A. Proofs

The proof establishes boundedness and conditional moment properties under the stated assumptions, then uses these facts to derive the required convergence statements.

  • A. Proofs: The subsample matrix estimator is conditionally unbiased, with E(ˇJ_X|F_n) = J_X.The proof applies conditional expectation and the tower property.
  • A. Proofs: Chebyshev’s inequality converts the conditional moment bounds into the convergence statement required by the lemma.The proof checks componentwise bounds before applying Chebyshev’s inequality.
  • A. Proofs: The remaining intermediate claims are obtained by analogous conditional-probability arguments under the same assumptions.The proof records these steps as direct consequences of the preceding bounds.
  • A. Proofs: Compactness and continuity imply that u(β^T x_i) and relevant derivatives are uniformly bounded.The argument invokes Assumption (H.1) and continuity on a compact parameter space.

A.1 Proof of Theorem 1

Theorem 1 shows that the subsample likelihood approximates the full-data likelihood closely enough for the subsample estimator to approach the full maximum likelihood estimator.

  • A.1 Proof of Theorem 1: As r →∞, n−1L*(β)−n−1L(β) →0 in conditional probability given F_n.The proof combines the preceding convergence result with compactness of the parameter space.
  • A.1 Proof of Theorem 1: The uniqueness of the full-data maximum likelihood estimator transfers the likelihood approximation into estimator consistency.The argument applies an argmax theorem to the continuous objective functions.
  • A.1 Proof of Theorem 1: Taylor expansion and bounded remainder terms provide the intermediate control needed for the estimator-rate result.The proof bounds derivatives and remainder terms using the regularity assumptions.
  • A.1 Proof of Theorem 1: The estimator difference satisfies ˜β − ˆβ_MLE = O_P|F_n(r^−1/2).Thus, conditionally on F_n, the subsample estimator approaches the full-data estimator at rate r^−1/2.

A.2 Proof of Theorem 2

Theorem 2 is proved by establishing a conditional central limit theorem for the subsampling contributions and then combining it with consistency and matrix convergence.

  • A.2 Proof of Theorem 2: Given F_n, the subsampling variables η_1,…,η_r are independent with mean zero and controlled variance.The variance is O_P(1) under the stated assumptions.
  • A.2 Proof of Theorem 2: The Lindeberg–Feller conditions hold in probability, yielding the required conditional convergence in distribution.The proof applies the conditional Lindeberg–Feller central limit theorem.
  • A.2 Proof of Theorem 2: The conditional limit is combined with earlier estimator and matrix results through Slutsky’s theorem.The proof explicitly invokes Lemma 1, Theorem 1, and Slutsky’s theorem.
  • A.2 Proof of Theorem 2: The equality condition in the Cauchy–Schwarz bound identifies when the associated optimality bound is attained.The proof uses the equality case of Cauchy–Schwarz to characterize the relevant condition.
  • A.2 Proof of Theorem 2: Defining 0/0 = 0 is equivalent to removing observations whose specified residual term is zero from V_c.This convention handles zero terms in the expression for V_c.

A.4 Proof of Theorems 5 and 6

The proofs of Theorems 5 and 6 use randomized matrix-product approximation and Frobenius-norm concentration to control subsample estimating equations and establish finite-sample results.

  • A.4 Proof of Theorems 5 and 6: The Monte Carlo algorithm samples matrix elements independently according to a discrete distribution and constructs a subsampled approximation to a matrix product.The sampling distribution determines which columns of A and corresponding rows of B are retained.
  • A.4 Proof of Theorems 5 and 6: Appropriate sampling probabilities make the reconstructed product CR a high-probability approximation to AB under a Frobenius-norm concentration inequality.Lemma 2 supplies the concentration result used by the later proofs.
  • A.4 Proof of Theorems 5 and 6: The maximum likelihood estimator satisfies an estimating equation expressed through the model’s diagonal derivative matrix.The proof defines the diagonal entries using ˙ψ(u(X^Tβ)).
  • A.4 Proof of Theorems 5 and 6: The subsampled estimator ˜β is defined through the sampled matrix S using probabilities π_mV, with the π_mVc case handled similarly.The proof then bounds the resulting estimating-equation error using matrix inequalities.
  • A.4 Proof of Theorems 5 and 6: Combining the concentration and perturbation bounds yields the stated theorem result with probability at least 1 − ε.The final step follows after applying the preceding inequalities.
  • A.4 Proof of Theorems 5 and 6: Weyl’s inequality and spectral-norm bounds control the difference between full-data and subsampled matrix terms.The argument bounds the perturbation through products of derivative and covariance-related matrix norms.

A.5 Proof of Theorem 7

The proof of Theorem 7 analyzes the adaptive two-step estimator under the initial subsampling strategy and establishes its asymptotic behavior using conditional convergence arguments and regularity assumptions.

  • Proof setup: The proof begins with the πmV case, defining the first- and second-step log-likelihoods used in the two-step algorithm.The πmVc case is treated analogously.
  • Supporting lemma: Under Assumptions (H.1)–(H.4), Lemma 3 provides the conditional asymptotic result needed for Theorem 7.The proof repeatedly invokes Lemma 3 together with the first-stage subsampling conditions.
  • Rate control: The proof controls remainder terms and inverse-information differences using bounds derived from the assumptions and Lemmas 1 and 3.These bounds yield the stated OP(r^-1) and OP|Fn(r^-1/2) rates for relevant terms.
  • Assumptions: The initial subsampling probabilities satisfy the required conditions, allowing Theorem 2 to be applied to the first-stage estimator.The parameter space is compact and the full-data maximum likelihood estimator is the unique global maximum of the continuous convex objective.
  • Conclusion: As r0/r → 0, the desired asymptotic result follows after combining the lemmas and noting the required conditional convergence.The proof concludes by applying the established rate relations.

A.6 Proof of Theorem 8

The proof of Theorem 8 establishes conditional asymptotic normality for the adaptive subsample estimator by verifying moment and Lindeberg–Feller conditions, then combining the resulting limit with Slutsky’s theorem.

  • Proof setup: The proof starts with the π̃mVc case and states that the π̃mV case follows by analogous modifications.The adaptive probabilities replace the full-data estimate with the pilot estimator β̃0.
  • Conditional CLT: The variables used in the conditional central limit argument are i.i.d. with zero mean and a stated variance.The proof derives the needed bounds from assumptions (H.1), (H.3), and (H.5).
  • Conditional CLT: The Lindeberg–Feller conditions hold in probability, yielding the conditional central limit theorem given Fn and β̃0.The proof invokes the Lindeberg–Feller theorem after verifying the relevant conditions.
  • Approximation: The difference between inverse information matrices is OP|Fn(r^-1/2), based on Lemma 3 and Theorem 7.This controls the replacement of the pilot-based quantity by the corresponding target quantity.
  • Conclusion: Slutsky’s theorem completes the proof of the desired asymptotic result.The preceding conditional limit and approximation results are combined at the final step.

B. Additional Simulation Results

Additional simulations examine how the allocation between pilot and main subsamples and the total subsample size affect estimation accuracy across Poisson and negative binomial regression settings.

  • Subsample allocation: The two-step approach performs best when r0/r is around 0.2, motivating comparisons across total sample sizes at that fixed ratio.The simulations vary r0 + r while holding r0/r = 0.2.
  • Methods and models: The additional experiments compare the proposed method with uniform, leverage, and adjusted leverage subsampling for Poisson and negative binomial regression.The corresponding results are reported in Figures 6 and 7.
  • Allocation comparison: The proposed methods are slightly better than using a fixed pilot size r0 = 200, but the improvement is not significant.This comparison is based on the results shown in Figures 6 and 7.
  • Idealized setting: An idealized simulation with the full-data MLE available varies r over 10p, 15p, 20p, 25p, and 30p and compares against uniform subsampling.The setup is intended to isolate pilot-subsampling effects.
  • Covariate conditioning: With fixed covariate condition number one, r = 10p and r = 20p have similar performance, while eMSE increases as the condition number grows.These results are reported for the additional Poisson regression simulation cases.
Loading 1806.06761v2…