Fast Perfect $L_2$ Sampling

Introduction

The scale of modern data systems often makes it infeasible to store or repeatedly process an entire dataset. In many settings, such as network monitoring and large-scale database systems, data instead arrives as a continuous stream of updates. The streaming model of computation formalizes this setting by requiring algorithms to process the input using memory far smaller than the size of the data itself. In the standard turnstile model, the stream implicitly defines a high-dimensional frequency vector $x \in\mathbb{R}^n$ that is initialized to the all zeros vector and evolves through updates of the form $x_i \gets x_i + \Delta$, where $\Delta$ may be positive or negative. This model naturally captures dynamic settings in which data is both inserted and deleted, such as updates to database records or changes in network traffic.

A central goal in the streaming model is to design primitives that extract useful information from the frequency vector while respecting stringent space constraints. One particularly powerful primitive is $L_p$ sampling, which returns an index $i \in[n]$ with probability proportional to $(1 \pm\varepsilon) \cdot\frac{|x_i|^p}{\lVert x \rVert_p^p} + \frac{1}{\operatorname{poly}(n)}$. Informally, an $L_p$ sampler selects coordinates according to their contribution to the $p$-th frequency moment, biasing the sample toward larger entries as $p$ increases. The sampler is called perfect if $\varepsilon= 0$ and approximate otherwise. Because $L_p$ sampling captures structure in the frequency vector beyond global aggregates by extracting representative coordinates according to the geometry induced by the $L_p$ norm, $L_p$ sampling has emerged as a powerful primitive in a variety of sublinear and streaming algorithms. In particular, it has found applications in a broad spectrum of settings, including network monitoring and anomaly detection [11 HNG+07, 14, 44, 33], distributed computation frameworks [7, 46 JSTW19, 49, 4], database query processing tasks [20, 31 GM98, 5, 40, 18, 30, 21], and techniques for compact representation and summarization of large datasets [12, 9, 10, 1, 23, 34, 8, 22 MRWZ20, 29, 38, 16, 36, 32]. Despite this broad applicability, the complexity and behavior of $L_p$ sampling depend strongly on the choice of $p$, and different regimes require fundamentally different algorithmic ideas.

Among the family of $L_p$ samplers, the case $p=2$ occupies a particularly important position. From a geometric perspective, the $L_2$ norm underlies Euclidean structure and is therefore deeply connected to many algorithmic techniques based on linear algebra and sketching. Within database systems, the second frequency moment $F_2=\lVert x\rVert_2^2$ corresponds directly to self-join size, making $L_2$ sampling a natural primitive for approximate query processing, join estimation, and regression over joins without explicitly materializing them, where sampling-based techniques are commonly used [27]. Moreover, it is used in a number of subroutines such as streaming algorithms for $L_p$ norm estimation for $p>2$ [50, 48]. Thus, progress on $L_2$ sampling directly influences the complexity of a much broader class of problems. In particular, $p=2$ is the largest exponent for which polylogarithmic space $L_p$ sampling is possible, as $\widetilde{\Omega}(n^{1-2/p})$ space is necessary for $p>2$ [15, 48, 24].

Previous work. Despite its central role, constructing efficient $L_p$ samplers has been well-studied in the streaming literature. In the insertion-only setting, the classical reservoir sampling technique [47] gives an exact $L_1$ sampler using $O(\log n)$ bits of space with no distortion, both multiplicative and additive. However, beyond the $L_1$ case or outside insertion-only streams, the problem becomes significantly more complex and the question of whether such samplers exist was explicitly posed by [6], and subsequently resolved in the affirmative for $p\leq2$ in the approximate setting [37]. These initial constructions used $\operatorname{poly}(1/\varepsilon,\log n)$ space, and were later improved to near-optimal bounds by [2] and by [25], achieving roughly $O\left(\frac{1}{\varepsilon^{\max(1,p)}}\log^2 n\right)$ space for $p\in(0,2)$ and $O\left(\frac{1}{\varepsilon^2}\log^3 n\right)$ for $p=2$. While these results established near-optimal space complexity for approximate sampling, they left open finer questions regarding distortion and randomness requirements. This gap was later resolved by [28], who introduced the first perfect $L_p$ samplers for $p\in(0,2)$ with $\widetilde{O}(\log^2 n)$ space and for $p=2$ with $\widetilde{O}(\log^3 n)$ space.

Update time. Despite resolving the space complexity, the algorithms of [28] and related work [5] suffer from a severe computational bottleneck, duplicating each coordinate $n^c$ times and assigning each copy an independent exponential weight. Consequently, a single update requires $\Omega(n^c)$ structural operations, so that the update time per element is polynomial in $n$, which can be prohibitively large in many settings. This sparked work on fast perfect $L_p$ sampling. For insertion-only streams, [41] achieved $O(1)$ update time for $p\leq1$ in the random oracle model. Subsequently, [43] introduced a continuous simulation framework for general turnstile streams, achieving polylogarithmic update time for $p<2$. Unfortunately, the framework of [43] does not extend to the critical case $p=2$, where the limiting behavior underlying the continuous simulation approach changes qualitatively and the techniques used for $p<2$ no longer apply (see § 3 for more details).

At first glance, this distinction may appear to be merely a change of regime. However, the situation is less clear from a complexity perspective. Indeed, while a combination of techniques from [24, 42] showed an $\Omega(\log^3 n)$ space lower bound for perfect $L_2$ sampling, perfect samplers for $p < 2$ are known using only $\widetilde{O}(\log^2 n)$ space [28, 43]. Thus, it is natural to ask whether the additional memory complexity of the $L_2$ regime fundamentally obstructs fast algorithms, or whether one can still achieve polylogarithmic update time while maintaining optimal space.

Our contributions. In this paper, we resolve the open question of whether perfect $L_2$ sampling admits both optimal space and fast update time in the turnstile model. We introduce a new approach for implementing perfect $L_2$ sampling without relying on the explicit polynomial duplication used in previous constructions to achieve the first perfect $L_2$ sampler with expected polylogarithmic update time and optimal space complexity.

We remark that the output $\perp$ corresponds to a failure event that occurs with probability at most $\delta$; this is standard in perfect sampling constructions and can be reduced arbitrarily by repetition without affecting the conditional sampling guarantee. Conceptually, our result shows that the additional complexity of the $L_2$ regime does not fundamentally obstruct efficient algorithms. Although the methods developed for $p < 2$ do not directly extend to $p = 2$, we show that this qualitative change in behavior need not incur an algorithmic penalty. Thus, despite occupying a distinct space complexity regime, perfect $L_2$ sampling remains achievable with both optimal space and efficient updates.

Motivation and Applications

Statistical indistinguishability. Sampling is a fundamental primitive for extracting information from massive datasets. In particular, $L_2$ sampling has been used as a subroutine for central streaming tasks such as heavy-hitters, norm and moment estimation, duplicate detection, and data summarization [37, 2, 25, 3, 28]. Since these summaries effectively represent the dataset itself, it is important that the samples accurately reflect the target distribution. Approximate samplers permit multiplicative distortion in the sampling probabilities. While this distortion may be small for an individual sample, the effect can accumulate when many samples are repeatedly combined or reused throughout a larger computational pipeline. Consequently, hidden statistical bias may propagate into downstream procedures and alter their behavior. Perfect $L_2$ samplers avoid this issue by reducing the additive distortion in the output distribution to an arbitrarily small inverse-polynomial quantity while preserving essentially the same asymptotic space complexity. For example, if an application extracts $n^c$ samples for some constant $c > 0$, then choosing an additive error of $n^{-(c+100)}$ ensures that the resulting joint distribution remains statistically indistinguishable from one obtained using truly exact samples. This makes perfect samplers particularly attractive as reusable black-box primitives.

Distributed databases and query processing. The distinction between approximate and perfect sampling becomes especially important in distributed database systems. Modern data processing architectures partition datasets across many machines, where local summaries are constructed independently and later combined into global estimates. Sampling provides a natural mechanism for constructing such summaries while avoiding communication of the raw data itself. However, local inaccuracies need not remain local. Even small distortions introduced by approximate samplers can accumulate when summaries are merged, reused, or incorporated into later stages of processing. This phenomenon is particularly relevant in relational settings involving repeated statistical estimation, joins, and self-join computations, where the output of one randomized procedure often becomes the input to another. Perfect samplers substantially reduce such accumulated bias and therefore preserve the fidelity of both local and global statistical summaries.

Information leakage and adaptive interaction. Approximate sampling probabilities may also reveal unintended information about the underlying dataset. Since the permitted multiplicative distortion can depend on global properties of the data, repeated observations of samples may allow an observer to infer information beyond the sampled index itself. For instance, suppose a sampler increases the total probability mass of a large set of coordinates by a factor of $(1+\varepsilon)$ whenever a property $\mathcal{P}$ holds and decreases it by $(1-\varepsilon)$ otherwise. Then an observer could estimate whether $\mathcal{P}$ holds simply by collecting sufficiently many samples and measuring aggregate frequencies. Perfect $L_2$ samplers mitigate this issue because their output distributions differ from the ideal distribution only by an arbitrarily small additive amount. As a result, repeated interactions with the sampler reveal essentially no additional information beyond that inherent to the target distribution itself.

Algorithmic and Technical Overview

In this section, we provide an overview of our algorithmic framework and the technical novelties required to achieve polylogarithmic update time for perfect $L_2$ sampling, c.f., Figure 1. To contextualize our approach, it is instructive to first review the mechanics and limitations of previous sampling frameworks before explaining how we overcome them.

Exponential random variables and duplication. The standard approach [28, 43, 48] for perfect $L_p$ sampling involves scaling each stream coordinate $x_i$ by an inverse exponential random variable $e_i^{-1/p}$. Due to the max-stability property of exponential random variables [2], the probability that coordinate $i \in[n]$ yields the largest scaled value $\frac{|x_i|}{e_i^{1/p}}$ is exactly $\frac{|x_i|^p}{\lVert x\rVert_p^p}$. To recover this maximum element in a data stream, the coordinate is identified using a heavy-hitters data structure such as CountSketch, coupled with a statistical test to ensure the maximum is reliably isolated from the background noise. However, the probability of passing this test is heavily correlated with the identity of the maximum element itself, distorting the sampling probabilities by an $\Omega(1)$ multiplicative bias that destroys the “perfect” sampling guarantee. To handle this issue, [28] duplicated each coordinate $i \in[n]$ a total of $k = n^{C}$ times, scaling each duplicate by its own independent exponential random variable. This flattens the background noise so that the statistical test becomes nearly independent of the identity of the maximum element, recovering the desired perfect sampling distribution. Unfortunately, explicitly routing $n^{C}$ duplicate coordinates into the sketch requires a prohibitively high update time of $\Omega(n^{C})$.

Algorithm Pipeline: Fast Perfect $L_2$ Sampling

1. Explicit Head: Conceptually duplicate each update $n^{C}$ times, but explicitly generate only the top $\tau= \Theta(\log n)$ exponential weights.

2. Tail Simulation: Compress the remaining copies into a variance $V_{\mathrm{tail}}$, sampled via fast quadrature of its 1-stable continuous CDF.

3. Sketch Injection: Update a dense Gaussian CountSketch, injecting the head directly and the aggregate tail as a single $\mathcal{N}(0,V_{\mathrm{tail}})$.

4. Recovery: Run a statistical test and extract the maximum coordinate if it succeeds.

Fig. 1: Overview of the framework, reducing update time to an expected $\operatorname{poly}(\log n)$.

Algorithmic framework: explicit heads and simulated tails. Our goal is to preserve the exact statistical behavior of this duplication process while avoiding its explicit realization. Our approach is based on a simple structural observation: although the duplicated coordinates are numerous, only a small number contribute significantly to the sketch on any given update. We therefore separate the duplicated copies into two components. For a coordinate $i \in[n]$, the largest $\tau= \Theta(\log n)$ duplicated coordinates are generated and maintained explicitly (the “head”), while the remaining $m = k - \tau$ coordinates are aggregated into a single random contribution (the “tail”).

Crucially, this decomposition integrates seamlessly into the sketching framework. Following the approach of [43], we utilize a dense Gaussian CountSketch. Because Gaussian random variables are exactly 2-stable, the aggregate contribution of the uninstantiated tail duplicates in each sketch row behaves identically to a single Gaussian variable $\mathcal{N}(0,V_{\mathrm{tail}})$. Rather than explicitly materializing the $m$ tail coordinates, we can compress their contribution into a single sampled variance parameter $V_{\mathrm{tail}}$. Thus, as detailed in Algorithm 1, the algorithm maintains only the explicit heavy coordinates $\mathcal{H}_i$, while the remaining copies are implicitly injected into the sketch through the sampled tail variance.

The explicit heavy coordinates $\mathcal{H}_i$ and the necessary truncation parameters can be generated efficiently using standard properties of exponential spacings, as shown in Algorithm 2.

Sampling the continuous tail distribution and the $L_2$ barrier. The core algorithmic challenge lies in generating the variance parameter $V_{\mathrm{tail}}$. We simulate the sum of the tail variables $V = \sum_{j=1}^{m} e_{i,j}$ as a single random draw from its limiting continuous distribution, $\mathcal{Z}$. We perform this via inverse transform sampling (Algorithm 3): we draw $U \sim\operatorname{Unif}(0,1)$ and compute $z^*$ such that $F_{\mathcal{Z}}(z^*) \approx U$, where $F_{\mathcal{Z}}$ is the Cumulative Distribution Function (CDF) of the limiting tail.

To evaluate $F_{\mathcal{Z}}(x)$ efficiently inside the binary search, we must rely on numerical inversion of the Gil-Pelaez integral. This requires obtaining the characteristic function of the tail distribution.

Algorithm 1 Perfect $L_2$ stream update processing

Input: Stream update $(i, \Delta)$

Output: Update the dense CountSketch associated with coordinate $i$