Reservoir-Sampling Distribution Estimation

• Reservoir-sampling-based distribution estimation is a framework that uses scalable, unbiased sampling to estimate distributions in high-dimensional or streaming-data scenarios.
• Techniques such as ReSWD and varoptₖ integrate importance weighting and variance-optimality to efficiently estimate sliced Wasserstein distances and subset-sums.
• The approach is applied in machine learning, computer vision, and network analysis, demonstrating empirical improvements in error reduction and computational speed.

Reservoir-sampling-based distribution estimation refers to a family of techniques that leverage the statistical properties of reservoir sampling—classically developed for scalable, unbiased sampling from data streams—in order to construct distribution estimators with provable variance guarantees and optimality properties. These methods are especially relevant in high-dimensional or streaming-data contexts, where memory and computation constraints preclude exact inference over the full data. Recent research has extended the classical reservoir sampling framework using importance weighting, variance-optimality criteria, and integration with measures such as sliced Wasserstein distances, enabling robust and scalable solutions to distribution matching and subset-sum estimation problems in statistics, machine learning, computer vision, and graphics.

1. Core Principles of Reservoir Sampling

Reservoir sampling is designed to maintain a sample of size $k$ from a potentially unbounded stream of data, such that every element thus far seen has a specified inclusion probability in the reservoir. For unweighted streams, classical one-pass algorithms (such as Vitter’s and Efraimidis–Spirakis methods) ensure uniform inclusion probability. For weighted streams, each item is assigned a nonnegative weight $w_i$, and inclusion probabilities must be proportional to $w_i$.

Weighted reservoir sampling generalizes the selection process such that the sampled reservoir constitutes a valid base for unbiased estimation. For each candidate element $\theta$ (e.g., a projection direction in distribution estimation), a random key $k = u^{1/w(\theta)}$ is generated, where $u$ is drawn uniformly from $(0,1)$. The reservoir then consists of the $k$ elements with the smallest keys, guaranteeing that each candidate survives in the reservoir with a probability proportional to $w(\theta)$ (Boss et al., 1 Oct 2025).

2. Sliced Wasserstein Distance Estimation via Reservoir Sampling

Sliced Wasserstein Distance (SWD), defined for probability measures $\mu,\nu$ on $\mathbb{R}^d$ as

$$\operatorname{SWD}_p(\mu,\nu) = \left( \int_{\theta\in S^{d-1}} W_p^p(\langle \theta,\mu\rangle,\langle \theta,\nu\rangle) \, d\theta \right)^{1/p}$$

where $W_p$ denotes the 1D Wasserstein distance of the projected distributions, is a scalable proxy for high-dimensional Wasserstein metrics. Monte Carlo (MC) estimators approach this integral by averaging over $L$ random projections $\theta_i \sim \mathrm{Uniform}(S^{d-1})$: $$\operatorname{SWD}_p(\mu,\nu) \approx \left( \frac{1}{L} \sum_{i=1}^{L} W_p^p(\pi_{\theta_i}\mu, \pi_{\theta_i}\nu) \right)^{1/p}$$ Despite unbiasedness, such MC estimators exhibit variance $\sim \operatorname{Var}[W_p^p] / L$, which may remain prohibitive in optimization or learning contexts.

The Reservoir SWD (ReSWD) estimator (Boss et al., 1 Oct 2025) integrates weighted reservoir sampling into the SWD estimation process. At each iteration, a reservoir $\mathcal{R}$ of $K$ projection directions is constructed using weighted sampling, with $w(\theta)$ set to the current 1D Wasserstein cost $D(\theta) = W_p(\pi_{\theta}\mu, \pi_{\theta}\nu)$. Self-normalized importance weights are computed to yield the estimator: $$\widehat{S}_p(\mu,\nu) = \sum_{i=1}^K w_i \cdot D(\theta_i)$$ with $w_i = [1/q(\theta_i)] / \sum_j [1/q(\theta_j)]$, where $q(\theta_i) \propto D(\theta_i)$ is the marginal inclusion probability under reservoir sampling.

3. Variance-Optimal Reservoir Sampling in Subset-Sum Estimation

Variance-optimal reservoir sampling, as formalized in the varopt$_k$ scheme (0803.0473), targets the problem of maintaining a reservoir of $k$ weighted items to enable unbiased estimation of the total weight $\sum_{i\in M} w_i$ for any subset $M$ of items. The algorithm maintains adjusted weights $\hat w_i$ such that the Horvitz–Thompson estimator over any $M$,

$\widehat S_M = \sum_{i \in S \cap M} \hat w_i,$

yields $E[\widehat S_M] = S_M$. The scheme is characterized by:

• Maintenance of a “threshold” $\tau_k$ determined by

$$\sum_{j \in S\cup\{\text{new}\}} \min\!\Bigl(1,\frac{\tilde w_j}{\tau}\Bigr)=k$$

• Eviction and adjustment procedures guaranteeing that $|S|=k$ at all times.
• Zero total sum variance: $\operatorname{Var}[\widehat S_{[n]}]=0$.
• Optimal minimization of

$V_m = E_{|M|=m}[\operatorname{Var}[\widehat S_M]]$

across all subset sizes $m$ among all possible schemes with $k$ samples.

The algorithm is efficient, offering $O(\log k)$ per-element update complexity and supporting merge operations for distributed or parallel data streams (0803.0473).

4. Unbiasedness, Variance Reduction, and Theoretical Guarantees

Reservoir-sampling-based estimators are provably unbiased. For SWD estimation, self-normalized importance sampling ensures

$$E[\widehat S_p] = E_{\theta\sim\mathrm{Uniform}}[W_p(\pi_{\theta} \mu, \pi_{\theta} \nu)]$$

while the empirical variance is reduced relative to plain MC estimation—empirically by up to 20–30% for a fixed number of projections (Boss et al., 1 Oct 2025).

In subset-sum estimation, varopt$_k$ yields strictly minimal average variance $V_m$ for all $m$ and supports tight worst-case bounds, such as

$$\operatorname{Var}[\hat w_i] \leq w_i W / k, \quad \operatorname{Var}[\widehat S_M] \leq (W/k) S_M$$

where $W = \sum_i w_i$. Notably, covariance terms between adjusted weights are identically zero.

A key feature is composability in distributed contexts: varopt$_k$ reservoirs can be merged (using adjusted weights and threshold recomputation) to yield the same statistical guarantees as if the complete data stream had been processed sequentially (0803.0473).

5. Algorithms: Pseudocode Structure and Computational Complexity

The ReSWD update algorithm at optimization step $t$ involves:

• (Optional) Time-decay of old keys: $k_i \gets k_i \cdot \exp(-\text{age}/\tau)$.
• Drawing $M$ new directions $\theta$.
• Computing costs $D(\theta)$ and keys $k(\theta)$ for the union of old and new directions.
• Retaining $K$ directions with smallest keys.
• Computing inclusion probabilities $q_i$ and self-normalized weights $w_i$.
• Performing an effective sample size (ESS) check to trigger reservoir resets if necessary.

Complexity per update is $O((K+M) n \log n)$ for data of size $n$ and $O((K+M)\log K)$ for key sorting. In practice, setting $K+M\approx L$ (the MC sample size) ensures similar asymptotic costs as plain SWD, with modest overhead for structural maintenance (Boss et al., 1 Oct 2025).

In varopt$_k$, reservoir updates and merges are executed in $O(\log k)$, with amortized constant-time operations possible under specific implementations (0803.0473).

6. Empirical Performance and Applications

Empirical evaluation in (Boss et al., 1 Oct 2025) demonstrates that ReSWD achieves the lowest mean $W_1$ error in synthetic 3D-to-3D distribution matching at marginal computational cost (e.g., $0.622 \times 10^3$ $W_1$ error at $1.9\,$ms/step vs $0.670$–$0.733$ for baselines at $1.03\,$ms/step). In vision and graphics tasks such as color correction and diffusion guidance, ReSWD delivers measurable improvements in error metrics (e.g., RMSE reduced from $0.34\rightarrow0.31$, PSNR increased from $24.30\rightarrow24.64$) and efficiency (guidance for SD3.5 Large and Turbo yields $2$–$4\times$ speedup, $30$–$45\%$ lower $W_2$ color-distance) (Boss et al., 1 Oct 2025).

In streaming data applications, varopt$_k$ is widely used for network traffic analysis, streaming subset-sum estimation, and distributed statistics, offering strict statistical optimality and efficient support for parallel and mergeable computation (0803.0473).

7. Practical Considerations and Limitations

Reservoir size ($K$ in ReSWD, $k$ in varopt$_k$) and the number of new candidates per iteration ($M$) govern the balance between memory usage and adaptivity. Ablation studies in (Boss et al., 1 Oct 2025) suggest $K=64$, $M=8$ delivers robust performance. Decay parameters allow adaptation to nonstationary data; effective sample size reset heuristics prevent estimator degeneration by enforcing periodic full redraws.

ReSWD’s variance advantages decrease when $\mu\approx\nu$ (costs $D(\theta)$ nearly uniform); in this regime, ReSWD reverts to standard SWD. For extension beyond linear projections (e.g., learned kernels), no improvement was observed over random projections, attributed to the overwhelming search space dimensionality (Boss et al., 1 Oct 2025).

A plausible implication is that while reservoir-sampling-based estimation offers theoretically optimal, unbiased, and adaptive estimators for both streaming subset-sums and distribution matching objectives, its practical impact depends on task-specific parameterization and the structure of the data or distributions involved.