Subsample-and-Aggregate Methods

• Subsample-and-Aggregate is a methodological template that divides large datasets into small subsamples, computes local statistics, and aggregates them to address memory, privacy, and scalability challenges.
• It leverages tuning parameters, where larger subsample sizes reduce bias and increased aggregation lowers variance, ensuring efficient computation even with limited resources.
• The approach underpins diverse applications—from subbagging and adaptive data analysis to time-series privacy—demonstrating practical benefits and theoretical rigor in modern inference.

Subsample-and-aggregate denotes a family of procedures in which a dataset is broken into smaller subsamples, a local statistic is computed on each subsample, and those local outputs are combined by averaging, noise addition, or another aggregation rule. In one canonical big-data formulation, for i.i.d. observations $Z_1,\dots,Z_N\sim F_Z$, $m_N$ subsamples of size $k_N\ll N$ are drawn and the aggregate estimator is

$$\hat\theta_{k_N,m_N}=\frac1{m_N}\sum_{j=1}^{m_N}\hat\theta_{k_N,s^{(j)}}.$$

This same pattern appears in subbagging for memory-constrained estimation (Zou et al., 2021), in block-based subagging for scalable inference (Politis, 2021), in adaptive data analysis through subsampling queries that output few bits (Blanc, 2023), in time-series differential privacy through subsampling and noisy re-stitching in the time domain (Koga et al., 2022), and in recent private aggregation schemes based on robust winsorization (Ramsay et al., 23 Jan 2025). Across these settings, subsample-and-aggregate is not a single algorithm but a methodological template for trading full-sample computation against memory, privacy, robustness, and generalization.

1. Canonical forms and notation

In the random-subsample formulation studied by Zou et al., the full dataset has size $N$, the parameter $\theta_0\in\mathbb R^d$ is defined by estimating equations $E[\psi_{\theta_0}(Z)]=0$, and the analyst draws $m_N$ subsets of size $k_N$ from the collection $S_N^{k_N}$ of all $m_N$0 subsets. Each subset is sampled uniformly without replacement within the subset and with replacement across subsets, so overlap between drawn subsamples is permitted. For each subsample $m_N$1, one computes the Z-estimator $m_N$2 solving $m_N$3, and then averages the resulting estimators. When all $m_N$4 subsets are used, the construction is a complete U-statistic; when only $m_N$5 subsets are used, it is an incomplete U-statistic of infinite order (Zou et al., 2021).

A distinct but related formulation replaces random subsets by non-random block subsamples. Politis considers $m_N$6, a subsample size $m_N$7, an overlap parameter $m_N$8, and the blocks

$m_N$9

The subagging estimator is then the simple average $k_N\ll N$0, where $k_N\ll N$1. The same block collection also supports subsampling-based distribution estimation through

$k_N\ll N$2

This shows that subsample-and-aggregate can be instantiated with either random or deterministic subsamples (Politis, 2021).

A third canonical description appears in the privacy literature: classical subsample-and-aggregate partitions the rows of a dataset into random subsets, computes a local statistic on each, then aggregates them with noise. The time-domain mechanism of Huang et al. takes an analogous view in the time dimension rather than the row dimension, replacing row subsampling by random selection of time indices (Koga et al., 2022). This suggests that the essential object is the decomposition-and-recombination pattern, not a fixed sampling geometry.

2. Statistical subbagging under memory constraints

The subbagging theory in Zou et al. is built around the asymptotic regime $k_N\ll N$3, $k_N\ll N$4, and $k_N\ll N$5, with regularity conditions C1–C4 on smoothness, moments, and nonsingularity. Defining

$k_N\ll N$6

their main asymptotic result states that if

$k_N\ll N$7

then

$k_N\ll N$8

Relative to the full-sample estimator, the asymptotic variance is inflated by the factor $k_N\ll N$9; when $$\hat\theta_{k_N,m_N}=\frac1{m_N}\sum_{j=1}^{m_N}\hat\theta_{k_N,s^{(j)}}.$$0, this inflation vanishes and the full-sample asymptotic variance is recovered (Zou et al., 2021).

The tuning parameters $$\hat\theta_{k_N,m_N}=\frac1{m_N}\sum_{j=1}^{m_N}\hat\theta_{k_N,s^{(j)}}.$$1 and $$\hat\theta_{k_N,m_N}=\frac1{m_N}\sum_{j=1}^{m_N}\hat\theta_{k_N,s^{(j)}}.$$2 encode the central trade-off. Larger $$\hat\theta_{k_N,m_N}=\frac1{m_N}\sum_{j=1}^{m_N}\hat\theta_{k_N,s^{(j)}}.$$3 reduces the $$\hat\theta_{k_N,m_N}=\frac1{m_N}\sum_{j=1}^{m_N}\hat\theta_{k_N,s^{(j)}}.$$4 bias of each subsample estimator but increases memory use. Larger $$\hat\theta_{k_N,m_N}=\frac1{m_N}\sum_{j=1}^{m_N}\hat\theta_{k_N,s^{(j)}}.$$5 decreases the variance inflation factor $$\hat\theta_{k_N,m_N}=\frac1{m_N}\sum_{j=1}^{m_N}\hat\theta_{k_N,s^{(j)}}.$$6 but increases CPU time. Algorithm 1 in the paper requires $$\hat\theta_{k_N,m_N}=\frac1{m_N}\sum_{j=1}^{m_N}\hat\theta_{k_N,s^{(j)}}.$$7, for example $$\hat\theta_{k_N,m_N}=\frac1{m_N}\sum_{j=1}^{m_N}\hat\theta_{k_N,s^{(j)}}.$$8 with $$\hat\theta_{k_N,m_N}=\frac1{m_N}\sum_{j=1}^{m_N}\hat\theta_{k_N,s^{(j)}}.$$9, together with $N$0. Algorithm 2 adds first-order bias correction, reducing bias from $N$1 to $N$2, and thereby permits $N$3 with $N$4 while retaining $N$5-consistency.

The computational profile is explicitly favorable to out-of-core analysis. The procedure loads one subsample at a time, so memory is $N$6; runtime is $N$7, where $N$8 is the cost of solving the estimating equations on $N$9 points; subsampling itself costs $\theta_0\in\mathbb R^d$0; and the $\theta_0\in\mathbb R^d$1 solves can be performed independently. In the American airline example with $\theta_0\in\mathbb R^d$2 flight records, full-sample logistic regression required memory of approximately $\theta_0\in\mathbb R^d$3 GiB, whereas Algorithm 1 with $\theta_0\in\mathbb R^d$4, $\theta_0\in\mathbb R^d$5, and $\theta_0\in\mathbb R^d$6 used about $\theta_0\in\mathbb R^d$7 MiB and produced estimates almost identical to the full-sample fit, with slightly larger standard errors by approximately the $\theta_0\in\mathbb R^d$8 factor predicted by theory. The example makes explicit that subbagging is not only asymptotically justified but practically motivated by memory constraints.

3. Scalable subagging and subsampling-based inference

Subsampling was developed in the 1990s as a general method for estimating the sampling distribution of a statistic $\theta_0\in\mathbb R^d$9, especially for confidence intervals and hypothesis tests, and Politis revisits it in the big-data regime where even drawing one large random subsample can be costly (Politis, 2021). The key observation is that a carefully chosen deterministic set of block subsamples can replace random subset generation while preserving the inferential role of subsampling. Under the assumptions $E[\psi_{\theta_0}(Z)]=0$0, $E[\psi_{\theta_0}(Z)]=0$1, $E[\psi_{\theta_0}(Z)]=0$2 with $E[\psi_{\theta_0}(Z)]=0$3, and $E[\psi_{\theta_0}(Z)]=0$4, the distribution estimator $E[\psi_{\theta_0}(Z)]=0$5 converges in probability to the limiting distribution $E[\psi_{\theta_0}(Z)]=0$6 of $E[\psi_{\theta_0}(Z)]=0$7 at each continuity point $E[\psi_{\theta_0}(Z)]=0$8, and

$E[\psi_{\theta_0}(Z)]=0$9

Thus $m_N$0 block subsamples suffice to estimate the distribution of the full-sample statistic at essentially the same computational order as computing the statistic once.

The same block family yields a subagging estimator whose behavior can be analyzed by bias-variance decomposition. When $m_N$1, the block statistics are independent and

$m_N$2

Under Assumption B, with $m_N$3 and $m_N$4 for $m_N$5, the mean-squared error satisfies

$m_N$6

This can be tuned either to match the full-sample rate $m_N$7 or, when the bias exponent $m_N$8 exceeds $m_N$9, to exceed the full-sample convergence rate. A plausible implication is that aggregation over smaller estimators can improve not only computational feasibility but also the statistical rate, provided the bias structure is favorable.

The paper also develops inference directly for the subagging estimator. Under a uniform $k_N$0-moment condition,

$k_N$1

and $k_N$2 can be estimated by the empirical variance of the block statistics. This yields confidence intervals centered at the subagged estimator rather than at the original $k_N$3. One common misconception is that random subsampling is essential to subagging; Politis’ construction shows that, at least in this framework, appropriately chosen non-random blocks are sufficient.

4. Adaptive data analysis and information-theoretic aggregation

In adaptive data analysis, the central difficulty is that an analyst may choose query $k_N$4 after observing the answers to queries $k_N$5, so standard nonadaptive concentration arguments no longer control generalization. Blanc studies a subsampling-based framework in which each query $k_N$6 operates on a random subsample of size $k_N$7 drawn without replacement from the dataset and returns a value in a small output space $k_N$8. The main information-theoretic statement is

$k_N$9

where $S_N^{k_N}$0, query $S_N^{k_N}$1 has arity $S_N^{k_N}$2, and $S_N^{k_N}$3. The paper then converts this mutual-information control into generalization guarantees for adaptively chosen test queries (Blanc, 2023).

Two mechanisms are developed in detail. For statistical queries $S_N^{k_N}$4, the mechanism repeatedly draws $S_N^{k_N}$5 points uniformly from the sample, flips Bernoulli coins with means $S_N^{k_N}$6, and returns their average. With

$S_N^{k_N}$7

the answers satisfy, with probability at least $S_N^{k_N}$8, the simultaneous error bound

$S_N^{k_N}$9

For approximate median queries, the sample is split into $m_N$00 disjoint blocks, each block is queried, and a binary-search procedure aggregates the resulting votes.

This framework differs from classical estimator averaging, but it remains recognizably subsample-and-aggregate: subsampling introduces limited dependence on any single data point, and aggregation converts many low-information local outputs into a global answer. The paper’s limitations are equally instructive. It assumes the queries are subsampling queries, requires $m_N$01 for the high-confidence $m_N$02 bounds, and studies generalization rather than runtime. A common misconception is that subsampling is merely a computational heuristic; in this setting it is the primary mechanism controlling information leakage and hence adaptivity-induced bias.

5. Time-domain subsample-and-aggregate under differential privacy

Huang et al. adapt subsample-and-aggregate to time-series data, where the sensitive object is not a row but an individual’s participation pattern across time. The data model observes $m_N$03 individuals over $m_N$04 discrete time steps, with participation vectors $m_N$05 and aggregate count signal

$m_N$06

Adjacent datasets differ in one individual’s entire participation vector, and the paper assumes that no individual can contribute at more than $m_N$07 time steps: $m_N$08 Without any assumption, the $m_N$09-sensitivity of the count vector is $m_N$10; under the participation bound it drops to $m_N$11. The paper then shows that random time subsampling reduces sensitivity further with high probability, because a Bernoulli($m_N$12) subsample of time indices is unlikely to retain all $m_N$13 active times of one individual. Specifically, for any $m_N$14,

$m_N$15

controls the event that the subsample contains more than $m_N$16 of an individual’s active times, and with probability at least $m_N$17 the subsampled count map has $m_N$18-sensitivity at most $m_N$19. By Hoeffding’s inequality, one may take

$m_N$20

so the effective sensitivity behaves like $m_N$21 up to logarithmic terms (Koga et al., 2022).

The paper also studies filter-plus-subsample mechanisms. A circulant linear time-invariant filter $m_N$22 with row $m_N$23-norm $m_N$24 is applied to the count signal, and the filtered series is then subsampled in time. A matrix-concentration argument based on Tropp’s theorem yields a high-probability sensitivity bound of the form

$m_N$25

Given such a bound, the Gaussian mechanism adds i.i.d. $m_N$26 noise with

$m_N$27

Algorithmically, the mechanism filters, subsamples time indices, adds Gaussian noise to the retained coordinates, and interpolates deterministically back to length $m_N$28.

The conceptual link to classical subsample-and-aggregate is explicit in the paper: the subsample step randomly selects time indices, the local statistic is the filtered count at each sampled time, and the aggregate step adds Gaussian noise coordinatewise and re-stitches the result. Empirically, on all three real datasets in the study, the subsample-only mechanism achieved the lowest mean absolute error among the reported methods: on PeMS traffic flow, Gaussian $m_N$29, DFT $m_N$30, and subsample-only $m_N$31; on Gowalla, $m_N$32, $m_N$33, and $m_N$34; on Foursquare, $m_N$35, $m_N$36, and $m_N$37. On synthetic data, performance improved as the sampling frequency increased, while the baselines remained flat, and with additional observation noise the filter-plus-subsample variant outperformed subsample-only. These results show that subsample-and-aggregate can be instantiated along the time axis rather than the sample axis.

6. Private aggregation, robustness, and the modern view of aggregation

Recent work emphasizes that the aggregation step is itself a major design variable. Chao et al. consider a corrupted sample $m_N$38, partition it into $m_N$39 disjoint subsets of size $m_N$40, compute a non-private statistic $m_N$41 on each subsample, and release

$m_N$42

where $m_N$43 is a differentially private mean aggregator (Ramsay et al., 23 Jan 2025). Their proposed aggregator is the private modified winsorized mean (PMWM). In the univariate setting, PMWM first estimates the $m_N$44- and $m_N$45-quantiles of the $m_N$46 via a differentially private quantile algorithm, then winsorizes the $m_N$47 into the resulting interval, averages, and finally adds either Laplace noise for pure DP or Gaussian noise for zCDP. Under adversarial contamination and the paper’s support and grid conditions, the estimator is minimax optimal, up to constants and $m_N$48 factors, over broad heavy-tailed and contaminated distribution classes.

For subsample-and-aggregate itself, the paper derives a finite-sample deviation bound that isolates the role of the subsample estimator’s bias and robustness. If $m_N$49 estimates $m_N$50, then with probability at least $m_N$51,

$m_N$52

The paper extracts two explicit insights from this bound: the optimal choice of subsamples depends on the bias of the estimator computed on the subsamples, and the rate of convergence of the subsample-and-aggregate estimator depends on the robustness of the estimator computed on the subsamples. If $m_N$53 is unbiased and its variance behaves like $m_N$54, then the stochastic term is $m_N$55; since $m_N$56, achieving the overall $m_N$57 rate requires choosing $m_N$58 large enough that the subsample bias is also $m_N$59.

Taken together, the recent literature corrects several persistent misconceptions. Subsample-and-aggregate does not require a plain arithmetic mean as the combiner: the aggregation may be a simple average, a noisy Gaussian release, or a private modified winsorized mean (Ramsay et al., 23 Jan 2025). It does not require random subsets in every implementation: deterministic blocks can be sufficient for both aggregation and inference (Politis, 2021). And it is not restricted to computational acceleration: in adaptive data analysis it controls mutual information (Blanc, 2023), while in time-series privacy it directly reduces sensitivity (Koga et al., 2022). A plausible implication is that the modern interpretation of subsample-and-aggregate is best understood as a modular architecture whose statistical, computational, and privacy properties are determined jointly by three choices: how the subsamples are formed, what local statistic is computed on each subsample, and how those local outputs are aggregated.