1 Introduction↩︎

The streaming model of computation has emerged as a central framework for studying algorithms that process massive amounts of sequential data. In this setting, an underlying dataset is implicitly defined by a stream of updates and an algorithm must maintain a compact summary to approximately answer queries or report relevant statistics about the data. Often, the algorithm is restricted to a single pass over the stream and must use space sublinear in both the stream length \(m\) and the universe size \(n\). This framework has proven remarkably effective for reasoning about large-scale data processing, and has led to a rich body of work spanning algorithm design, lower bounds, and applications.

In the classical streaming setting, a standard assumption in the literature is that the input stream is fixed in advance and independent of the algorithm’s internal randomness. However, this assumption is often violated in modern applications, where the input stream may be generated adaptively based on previous outputs of the algorithm. For instance, queries to a database may depend on previous answers, iterative scientific procedures may update their state using previously computed estimates, and feedback loops may arise in online recommendation platforms and financial systems. In such settings, the input stream may depend on the behavior of the algorithm itself (and on the internal randomness), thus raising the question of whether classical streaming guarantees continue to hold.

1.0.0.1 Adversarially robust streaming.

To capture this phenomenon, [1] introduced the adversarially robust streaming model, where future updates may be chosen adaptively based on previous outputs of the algorithm. Specifically, at each time \(t \in [m]\), the algorithm receives an update \((a_t, \Delta_t)\) to the underlying frequency vector \({\mathbf{x}}\in \mathbb{Z}^n\), where \(a_t \in [n]\) is an index and \(\Delta_t \in \mathbb{Z}\) denotes an increment or decrement to coordinate \(a_t\) of \({\mathbf{x}}\) by \(\Delta_t\). Concretely, the frequency vector \({\mathbf{x}}\) is defined implicitly by \(x_i^{(t)} = \sum_{t' \leq t: a_{t'} = i} \Delta_{t'}\). The algorithm must then output an estimate of a function \(f({\mathbf{x}}^{(t)})\) of the current prefix. Crucially, the algorithm must remain accurate at every step, even when each update \((a_t, \Delta_t)\) is chosen adaptively based on the entire transcript of previous updates and responses of the algorithm.

Informally, the model captures an interaction between a randomized algorithm \({\mathcal{A}}\) and an adaptive adversary that constructs a sequence of updates \(\{u_1,\ldots,u_m\}\), where each update \(u_t\) may depend on all previously observed updates and outputs. The interaction proceeds as follows:

1. For each \(t \in [m]\), the adversary selects an update \(u_t\), where the choice of \(u_t\) may be an arbitrary (possibly randomized) function of \(u_1,\ldots,u_{t-1}\) and \(Z_1,\ldots,Z_{t-1}\).

2. Upon receiving \(u_t\), the algorithm \({\mathcal{A}}\) updates its internal state to incorporate \(u_t\).

3. The algorithm then outputs an estimate \(Z_t\) of \(f({\mathbf{x}}^{(t)})\) prior to the next update.

Designing algorithms that remain accurate under such adaptivity, while still using limited space and a single pass, has emerged as a central challenge in recent work [1]–[33].

1.0.0.2 Norm and distinct element estimation.

A canonical class of problems in the streaming model is the estimation of the frequency moments of the underlying data. These moments capture a variety of fundamental statistics of the data stream and are formally defined as follows. For a frequency vector \({\mathbf{x}}\in \mathbb{Z}^n\), the \(p\)-th frequency moment is given by \(F_p({\mathbf{x}})=\sum_{i\in[n]}|x_i|^p\) for \(p>0\), while \(F_0({\mathbf{x}})=|\{i\in[n] : x_i \neq 0\}|\) denotes the number of distinct elements in the stream. These quantities capture a range of important properties of the data; for example, \(F_1\) corresponds to the total mass of the stream, \(F_2\) is closely related to the collision probability, and \(F_0\) measures the support size of the underlying vector. More broadly, the family of \(F_p\) moments for \(p>0\) provides a rich spectrum of statistics, and via the identity \(\|{\mathbf{x}}\|_p = (F_p({\mathbf{x}}))^{1/p}\), directly corresponds to the problem of estimating \(L_p\) norms.

These problems have been studied extensively in the classical (non-adaptive) streaming model since the seminal work of [34], with tight space bounds known for a wide range of parameters [34]–[39]. For the distinct elements problem corresponding to \(p =0\), it is known that \(\tilde{\Theta}\left(\frac{1}{\varepsilon^2}\log n\right)\) bits of space is both necessary [34], [37] and sufficient [38]. Similarly for \(p\in(0,2]\), it is known that for all \(\varepsilon\in(n^{-1/2+c},1)\) for any constant \(c>0\), \(\Theta\left(\frac{1}{\varepsilon^2}\log n\right)\) bits of space is both necessary [34], [37] and sufficient [34], [37], [39]. Finally, for \(p>2\), it is known that \(\tilde{\Theta}\left(\frac{1}{\varepsilon^2}\cdot n^{1-2/p}\right)\) bits of space is both necessary [40], [41] and sufficient [42], [43]. Similar bounds are known for \(p\in[0,2]\) and integer \(p>2\) in the adversarially robust model when the adversary is restricted to \(\Delta_t\ge 0\) for all \(t\ge 0\), i.e., the insertion-only setting [44].

In contrast, significantly less is understood in the adversarially robust setting in the presence of deletions. A long line of work has shown that linear sketching, which underlies essentially all known turnstile streaming algorithms in the non-adaptive setting [45]–[49], is fundamentally limited in this setting, and cannot achieve even a constant-factor approximation to these quantities unless the sketching dimension is \(r=\Omega(n)\), i.e., the sketch requires near-linear space [23], [26], [50]. This has led to a significant gap in our understanding of the power and limitations of turnstile streaming algorithms under adaptive updates.

Recent work has begun to close this gap by demonstrating that it is possible to design sublinear-space adversarially robust algorithms in the turnstile model. In particular, it is now known that a \((1+\varepsilon)\)-approximation for the \(L_2\) norm can be achieved using \(\mathop{\mathrm{poly}}\left(\frac{1}{\varepsilon},\log n\right)\) space [32]. For \(p\neq 2\) and for the distinct elements problem, however, only constant-factor approximations were previously known in sublinear space, leaving a significant gap in our understanding of the achievable approximation guarantees under adversarial updates.

1.1 Our Contributions↩︎

In this work, we significantly advance the understanding of adversarial robustness in the insertion-deletion streaming model. Our main results establish that, despite strong impossibility results for linear sketches [23], [26], [50], accurate estimation of a broad class of frequency moments is still achievable in sublinear space under adaptive updates. We further develop a conceptual connection between adversarial robustness and classical sketching, showing that these two models are, in a weak sense, equivalent up to polynomial losses. We also show that our robust \(F_p\) estimation primitives can be combined with standard metric and norm embeddings to obtain the first adversarially robust sublinear-space algorithms for a variety of other problems. These include geometric and clustering objectives such as Earth Mover Distance and \(k\)-median cost, as well as statistical measures like Shannon entropy.

1.1.0.1 \(F_p\) moment estimation for \(p \in [0,2]\).

Our first main result resolves the approximability of frequency moments and distinct elements in the adversarially robust turnstile model for the entire regime \(p \in [0,2]\). Prior to our work, \((1+\varepsilon)\)-approximation was only known for \(p=2\), while for all other values of \(p \in [0,2)\) (including the distinct elements problem \(F_0\)), only constant-factor approximations were known in sublinear space [32].

[thm:lp:small] provides a unified and nearly optimal guarantee for all \(p \in [0,2]\), matching the space complexity of classical (non-adaptive) streaming up to polynomial factors [34], [37]–[39]. In particular, [thm:lp:small] gives the first \((1+\varepsilon)\)-approximation for the distinct elements problem in the adversarially robust turnstile model using polylogarithmic space, thereby closing a central open question in the area.

1.1.0.2 Equivalence between robustness and sketching.

Our second main result provides a structural explanation for the power and limitations of adversarially robust streaming algorithms. As previously discussed, linear sketching is known to be inherently non-robust to adaptive inputs, yet it underlies essentially all optimal oblivious streaming algorithms. We show that this tension can be formalized via a general equivalence theorem.

A key conceptual ingredient underlying this result is again the central role of \(L_1\). Building on [thm:lp:small], which gives a robust and accurate estimator for the \(L_1\) norm, we leverage a foundational connection between sketching and metric embeddings. Specifically, a seminal result of [51] establishes that for finite-dimensional normed spaces, admitting an efficient classical linear sketch is mathematically equivalent to admitting a low-distortion embedding into \(L_{1-\varepsilon}\) (and into \(L_1\) for norms closed under sum-product). By this characterization, any normed space that admits an efficient oblivious sketch must embed into such spaces. This allows us to simulate arbitrary linear sketches in the adversarially robust setting by composing these guaranteed embeddings with our robust \(L_p\) estimation primitives. Conversely, any adversarially robust algorithm can be converted into a classical sketch by compressing its interaction transcript. Together, these reductions yield a general equivalence between the two models.

[thm:equivalence] establishes a weak equivalence between the adversarially robust and oblivious models: any progress in one model translates moderately to the other. In particular, [thm:equivalence] shows that the existence of efficient adversarially robust algorithms implies new linear sketches, while classical sketching lower bounds extend to the robust setting. Thus, this connection provides a unifying framework for understanding which problems admit efficient robust algorithms, and explains the inherent difficulty of extending linear sketches to the adversarial setting.

A central theme underlying both our algorithms and our structural results is the role of \(L_1\) as a unifying geometry. At a technical level, our \((1+\varepsilon)\)-approximation for \(F_p\) in [thm:lp:small] is obtained by reducing to robust \(F_1\) estimation via suitable transformations. More broadly, rather than \(L_1\) merely serving as a convenient target space for a wide range of embeddings arising in theoretical computer science, the equivalence of [51] demonstrates that \(L_1\) (and \(L_{1-\varepsilon}\)) embeddability is the fundamental and unavoidable mechanism for sketching norms. This perspective is powerfully reinforced by [thm:equivalence], whose proof proceeds by extracting such embeddings from classical sketches and simulating them in the adversarially robust setting. Taken together, these results establish a compelling triad—adversarial robustness, oblivious sketching, and \(L_1\) embeddability are intrinsically tied together—and suggest that robust \(L_1\) estimation is a fundamental primitive from which more general adversarially robust streaming algorithms can be derived.

1.1.0.3 Additional applications.

Building on this viewpoint, we obtain a number of additional applications by combining known embeddings into \(L_p\) with our robust estimation framework. We highlight a few representative examples below.

Our first application is to the Earth Mover Distance (EMD), a central metric for comparing distributions over discrete domains. It is well-known that EMD over the grid \([\Delta]^d\) admits an embedding into \(L_1\) with distortion \(\mathcal{O}\left(d \log \Delta\right)\). By combining this embedding with our robust \(F_1\) estimation primitives, we obtain an adversarially robust approximation algorithm for EMD in the turnstile model.

Our second application is to the \(k\)-median clustering objective, which can also be reduced to \(L_1\) via standard embeddings. Using our framework, we obtain adversarially robust algorithms for maintaining approximate clustering solutions on insertion-deletion streams.

Finally, we consider entropy estimation, which can be reduced to estimating frequency moments. Using our robust \(F_p\) estimators as a black-box primitive, we obtain accurate entropy estimates even in the presence of adversarial updates.

Beyond these examples, our framework applies more broadly to any problem that admits a low-distortion embedding into \(L_1\) or can be reduced to frequency moment estimation. This includes a variety of geometric, statistical, and combinatorial estimation tasks, highlighting the generality of our approach and the central role of \(L_p\) embeddings in the adversarially-robust streaming model.

1.2 Technical Overview↩︎

1.2.1 \(F_p\) Estimation Algorithm for \(p<2\)↩︎

A natural starting point for [thm:lp:small] might be to generalize the existing adversarially robust \((1+\varepsilon)\)-approximate \(L_2\) estimation algorithm that handles insertions and deletions.

1.2.1.1 Estimator-corrector-learner framework.

At a high level, the \(F_2\) algorithm of [32] departs from classical linear sketching by introducing fresh randomness during the execution of the algorithm and organizing the computation around a estimator-corrector-learner paradigm. The stream is conceptually partitioned into blocks, and within each block we decompose the underlying frequency vector as \({\mathbf{x}}={\mathbf{z}}-{\mathbf{q}}\), where \({\mathbf{z}}\) corresponds to the portion of the stream before the block. For convenience, we write \(-{\mathbf{q}}\) as the updates of the block, so that \(\|{\mathbf{x}}\|\) can be interpreted as a metric measuring the distance between \({\mathbf{z}}\) and \({\mathbf{q}}\).

The algorithm maintains multiple sketches that correspond to three ingredients, which play different roles in this decomposition:

1. The estimator produces a candidate approximation to \(\|{\mathbf{z}}-{\mathbf{q}}\|_2^2\) via \(\|{\mathbf{z}}-{\mathbf{z}}'\|_2^2+\|{\mathbf{z}}'-{\mathbf{q}}\|_2^2\), where \({\mathbf{z}}'\) is a current guess for the vector \({\mathbf{z}}\), which arrived before the current block and is thus unknown from the perspective of the sketches initialized at the beginning of the block.

2. The algorithm also maintains a corrector, which uses a separate sketch on \({\mathbf{z}}-{\mathbf{q}}\) initialized at the beginning of the data stream, to verify whether the estimator is accurate. If the estimate is correct, the algorithm proceeds using the output of the estimator; otherwise, the corrector overrides the estimator and outputs an estimate for \(\|{\mathbf{z}}- {\mathbf{q}}\|_2^2\), and activates the learner.

3. Upon detecting an incorrect estimate, the learner updates the internal iterate \({\mathbf{z}}'\), refining its approximation to \({\mathbf{z}}\) using information revealed by the current query \({\mathbf{q}}\), for which the estimator was not accurate.

The key idea is that each time that the estimator is not accurate, this exposes structural information about \({\mathbf{z}}\) to the algorithm, allowing us to update the guess \({\mathbf{z}}'\) to move closer to \({\mathbf{z}}\) in a controlled manner. After the update, subsequent estimators operate with the refined iterate, and the process continues.

A crucial property of the framework is that the number of learner updates is small. Each query \({\mathbf{q}}\) can only correspond to an incorrect output by the estimator if \({\mathbf{q}}\) is aligned with some direction of \({\mathbf{z}}\) not captured by the guess \({\mathbf{z}}'\). For example, when the estimator errs via \(\|{\mathbf{z}}-{\mathbf{z}}'\|_2^2+\|{\mathbf{z}}'-{\mathbf{q}}\|_2^2>(1+\varepsilon)\cdot\|{\mathbf{z}}-{\mathbf{q}}\|_2^2\), this implies \[\begin{align} \langle{\mathbf{z}}-{\mathbf{z}}',{\mathbf{z}}'-{\mathbf{q}}\rangle\ge\varepsilon\cdot\|{\mathbf{z}}-{\mathbf{z}}'\|_2\cdot\|{\mathbf{z}}'-{\mathbf{q}}\|_2. \end{align}\] Thus, it can be shown that updating \({\mathbf{z}}'\) toward \({\mathbf{q}}\) by setting \({\mathbf{z}}''={\mathbf{z}}'+\alpha({\mathbf{z}}'-{\mathbf{q}})\) for some carefully chosen \(\alpha\approx\frac{\varepsilon\cdot\|{\mathbf{z}}-{\mathbf{z}}'\|_2}{\|{\mathbf{z}}'-{\mathbf{q}}\|_2}\) makes quantifiable progress toward \({\mathbf{z}}\): \[\begin{align} \label{eq:progress:L2} \|{\mathbf{z}}-{\mathbf{z}}''\|_2^2\le(1-\varepsilon^2)\cdot\|{\mathbf{z}}-{\mathbf{z}}'\|_2^2. \end{align}\tag{1}\] Since the frequency vector has entries bounded by \(\mathop{\mathrm{poly}}(n)\), it follows that the iterate \({\mathbf{z}}'\) converges to \({\mathbf{z}}\) after only \(\mathcal{O}\left(\frac{1}{\varepsilon^2} \log n\right)\) updates to \({\mathbf{z}}'\). Therefore, the corrector sketch only needs to be accurate for a small number of adaptive queries, which can be ensured with only polylogarithmic overhead in the total space complexity. Then, the final insight is that robustly estimating the second term \(\|{\mathbf{z}}' - {\mathbf{q}}\|_2^2\) is a problem of precisely the same form as the original problem, and thus it can be recursively computed by implementing the estimator-corrector-learner framework with estimator \(\|({\mathbf{z}}' - {\mathbf{q}}_1) - {\mathbf{u}}\|_2^2 + \|{\mathbf{u}}- {\mathbf{q}}_2\|_2^2\), where \({\mathbf{q}}= {\mathbf{q}}_1 + {\mathbf{q}}_2\) is decomposed so that segments \({\mathbf{q}}_1\) and \({\mathbf{q}}_2\) can be sketched using independent sketches.

1.2.1.2 An alternative update rule.

We remark that while [32] analyzed the algorithm for the update rule \({\mathbf{z}}'' \gets {\mathbf{z}}' + \alpha({\mathbf{z}}' - {\mathbf{q}})\), one can show the same convergence as in 1 for a modified update rule, given as follows: let \({\mathcal{L}}= \{{\mathbf{q}}_i\}_{i = 1}^L\) be a list of queries on which the estimator was previously inaccurate. If the estimator is not accurate on the current suffix \({\mathbf{q}}\), we update \[{\mathbf{z}}'' := \mathop{\mathrm{Proj}}_{\mathop{\mathrm{Span}}({\mathcal{L}}\cup \{{\mathbf{q}}\})} ({\mathbf{z}}),\] noting that the projection can be computed by solving a regression problem in the sketch space. In fact, this update rule has the nice (but not strictly necessary) property that once the iterate is updated to \({\mathbf{z}}''\), the new estimator \(\|{\mathbf{z}}- {\mathbf{z}}''\|_2^2 + \|{\mathbf{z}}'' - {\mathbf{q}}\|_2^2\) is precisely equal to \(\|{\mathbf{z}}- {\mathbf{q}}\|_2^2\) by Pythagorean theorem.

1.2.1.3 Possible approaches for \(F_p\) estimation.

In order to generalize the framework of [32] for other \(L_p\) norms, we must determine an appropriate decomposition of \(\|{\mathbf{x}}\|_p^p = \|{\mathbf{z}}- {\mathbf{q}}\|_p^p \approx g({\mathbf{z}}-{\mathbf{z}}') + g({\mathbf{z}}' - {\mathbf{q}})\) for some sketchable function \(g: \mathbb{R}^n \rightarrow \mathbb{R}\) with the following properties:

• The prefix \({\mathbf{z}}- {\mathbf{z}}'\) and suffix \({\mathbf{z}}' - {\mathbf{q}}\) can be sketched independently. Moreover, for \({\mathbf{q}}= {\mathbf{q}}_1 + {\mathbf{q}}_2\), the function \(g\) should have the property that a sketch for \(g({\mathbf{z}})\) can subsequently be used to track \(g({\mathbf{z}}- {\mathbf{q}}_1)\). Note that this is needed in order to implement the recursion, as described above.

• This decomposition should allow us to argue that if \(g({\mathbf{z}}-{\mathbf{z}}') + g({\mathbf{z}}' -{\mathbf{q}}) \not \in (1\pm\varepsilon) \cdot \|{\mathbf{z}}- {\mathbf{q}}\|_p^p\), the algorithm can recover more information about \({\mathbf{z}}\) from the current query \({\mathbf{q}}\).

Recall that in the \(L_2\) algorithm of [32], the crucial insight was that due to the Pythagorean properties of Euclidean space and the squared \(L_2\) norm, an inaccurate estimator directly implies that the dot product between \({\mathbf{z}}-{\mathbf{z}}'\) and \({\mathbf{z}}'-{\mathbf{q}}\) must be significant, thus enabling the algorithm to make progress in learning \({\mathbf{z}}\). Unfortunately, it is not clear how to establish a similar claim for \(L_p\) norms. For instance, one might try to write down the binomial expansion of \(\|{\mathbf{z}}- {\mathbf{q}}\|_p^p\) for integer \(p\), but then observe that incorrect outputs by the estimator might be caused by numerous correlations between coordinate-wise powers of \({\mathbf{z}}-{\mathbf{z}}'\), \({\mathbf{z}}'-{\mathbf{q}}\), and \({\mathbf{z}}-{\mathbf{q}}\). Thus, even by considering a more complex potential function of \({\mathbf{z}}-{\mathbf{z}}'\), it is not at all clear how to show convergence.

1.2.1.4 A hard instance: the Hadamard code.

Our first observation is that the iterate \({\mathbf{z}}'\) maintained by the learner of [32] uses a very specific update rule \({\mathbf{z}}'\gets{\mathbf{z}}'+\alpha({\mathbf{z}}'-{\mathbf{q}})\) and so \({\mathbf{z}}'\) can be written as a linear combination of a list \({\mathbf{u}}_1,\ldots,{\mathbf{u}}_L\) of previous queries on which the estimator failed. Hence, \({\mathbf{z}}'\) is contained within the subspace \({\mathbf{U}}\) spanned by the vectors \({\mathbf{u}}_1,\ldots,{\mathbf{u}}_L\). Thus, the distance \(\|{\mathbf{z}}'-{\mathbf{z}}\|_2^2\) between the iterate \({\mathbf{z}}'\) and the target vector \({\mathbf{z}}\) is at least the distance between \({\mathbf{z}}\) and the closest vector \({\mathbf{v}}\) to \({\mathbf{z}}\) contained within the subspace \({\mathbf{U}}\), i.e., \({\mathbf{v}}=\mathop{\mathrm{argmin}}_{{\mathbf{y}}\in\mathbb{R}^L}\|{\mathbf{U}}{\mathbf{y}}-{\mathbf{z}}\|_2^2\). Therefore, we could similarly upper bound the number of iterations required for convergence if at each step, we simply take the closest vector \({\mathbf{v}}\) to \({\mathbf{z}}\) contained within the subspace \({\mathbf{U}}\).

We could try adapting this approach to \(L_p\) norms by maintaining the list of queries \({\mathbf{u}}_1,\ldots,{\mathbf{u}}_L\) on which the estimator failed and iteratively updating \({\mathbf{z}}'\) to be the closest vector to the subspace \({\mathbf{U}}\) in \(L_p\) distance. However, it turns out that this approach fails catastrophically. Suppose \({\mathbf{z}}=\mathbf{1}^n\) and each query \({\mathbf{q}}\in\{-1,+1\}^n\) is a codeword in the Hadamard code not equal to \({\mathbf{z}}\), so that \(\sum_{i=1}^n q_i=0\). Hence, once such a query is added, the maintained subspace \({\mathbf{U}}=\mathop{\mathrm{Span}}({\mathbf{u}}_1,\ldots,{\mathbf{u}}_L)\) is contained in \(\left\{{\mathbf{y}}\in\mathbb{R}^n:\sum_{i=1}^n y_i=0\right\}\). The update step sets \({\mathbf{z}}'\) to be the minimizer of \(\min_{{\mathbf{y}}\in {\mathbf{U}}} \|\mathbf{1}^n - {\mathbf{y}}\|_p^p\). For \(p>1\), the function \(t \mapsto |1-t|^p\) is strictly convex and so by Jensen’s inequality, for any \({\mathbf{y}}\in {\mathbf{U}}\), \[\frac{1}{n}\sum_{i=1}^n |1 - y_i|^p \ge \left|1 - \frac{1}{n}\sum_{i=1}^n y_i\right|^p = 1,\] where we used that \(\sum_i y_i = 0\). Moreover, equality holds if and only if all coordinates of \({\mathbf{y}}\) are equal, which together with the zero-sum constraint forces \({\mathbf{y}}= \mathbf{0}^n\). Thus the unique minimizer is \({\mathbf{z}}' = \mathbf{0}^n\). Note that while \(\mathbf{0}^n\) is also the minimizer for \(p=2\), the algorithm avoids this trap because the initial estimator evaluates to \(\|{\mathbf{z}}-{\mathbf{z}}'\|_2^2 + \|{\mathbf{z}}'-{\mathbf{q}}\|_2^2 = n+n=2n\), which matches the true distance \(\|{\mathbf{z}}-{\mathbf{q}}\|_2^2 = \frac{n}{2}|1-1|^2 + \frac{n}{2}|1-(-1)|^2 = 2n\), meaning the corrector never flags these queries. For \(p \in (1, 2)\), however, the true distance is \(\|{\mathbf{z}}-{\mathbf{q}}\|_p^p = n 2^{p-1}\), which is a constant-factor multiplicative separation from the estimate of \(2n\). Consequently, \({\mathbf{z}}'\) remains \(\mathbf{0}^n\) throughout the process, so that it never progresses toward \({\mathbf{z}}\), even while each query is continuously flagged by the corrector and added to the list.

1.2.1.5 Embedding into \(L_2\) and regularized regression.

We first recall the existence of isometric embeddings \(f\) onto Hilbert spaces [52]–[54], so that for any \(p\in[0,2)\), we have for all vectors \({\mathbf{u}},{\mathbf{v}}\in\mathbb{R}^n\), \[\|{\mathbf{u}}-{\mathbf{v}}\|_p^p=\|f({\mathbf{u}})-f({\mathbf{v}})\|_2^2,\qquad \|f(\mathbf{0})\|_2=0,\] so that for the \(p\)-th moment of the vector \({\mathbf{z}}-{\mathbf{q}}\), it suffices to estimate \(\|f({\mathbf{z}})-f({\mathbf{q}})\|_2^2\). In the discussion below, we let \(N\) denote the embedding dimension of \(f(\cdot)\). Ideally, if we could apply the embedding on the fly to each arriving element in the data stream, then this problem would directly reduce to \(F_2\) estimation. Unfortunately, these embeddings are non-linear mappings into a higher dimensional space and it is not clear how to maintain embedded vectors in the streaming setting.

In fact, it turns out that we only require the existence of such embeddings to design our robust \(F_p\) estimation algorithm by simulating the estimator-corrector-learner paradigm implicitly in the embedded space. As such, we now describe the algorithm as though we have access to embedded vectors, although the true algorithm will strictly operate on the original vectors. By analogy to the \(L_2\) algorithm, suppose the estimator for \(\|{\mathbf{z}}- {\mathbf{q}}\|_p^p = \|f({\mathbf{z}}) - f({\mathbf{q}})\|_2^2\) is \(\|f({\mathbf{z}}) - {\mathbf{v}}\|_2^2 + \|{\mathbf{v}}- f({\mathbf{q}})\|_2^2\), where \({\mathbf{v}}= \mathop{\mathrm{Proj}}_{{\mathcal{L}}} f({\mathbf{z}})\) is the iterate that approximates \(f({\mathbf{z}})\) and \(\mathcal{L} = \{f({\mathbf{q}}_i)\}\) is the list of previous queries on which the estimator was not accurate ¹. Let \(|{\mathcal{L}}| = L\). Since \({\mathbf{v}}\) is a projection onto the span of the previous hard queries \(f({\mathbf{q}}_i)\), we have that \({\mathbf{v}}= \sum_{i = 1}^{L} \alpha_i \cdot f({\mathbf{q}}_i)\) for some coefficients \(\alpha_i \in \mathbb{R}\). Importantly, note that \({\mathbf{v}}\) does not correspond to an explicit vector in the original space, and may not admit a pre-image under the embedding \(f\).

Now, we can estimate the quantity \(\|f({\mathbf{z}}) - {\mathbf{v}}\|_2^2\) in the streaming setting by first observing that \[\begin{align} \label{eq:expansion:first:term} \|f({\mathbf{z}}) - {\mathbf{v}}\|_2^2 = \|f({\mathbf{z}})\|_2^2 - 2\sum_{i = 1}^{L} \alpha_i \cdot \langle f({\mathbf{z}}), f({\mathbf{q}}_i) \rangle + \sum_{1 \leq i,j \leq L} \alpha_i \cdot \alpha_j \cdot \langle f({\mathbf{q}}_i), f({\mathbf{q}}_j) \rangle. \end{align}\tag{2}\]

Then, we can simply estimate all of these right-hand side terms as follows:

• \(\|f({\mathbf{z}})\|_2^2 = \|{\mathbf{z}}\|_p^p\) can be estimated using an \(F_p\) estimation sketch for \({\mathbf{z}}\).

• \(\langle f({\mathbf{z}}), f({\mathbf{q}}_i) \rangle\) can be estimated by noting that \[\langle f({\mathbf{z}}), f({\mathbf{q}}_i) \rangle = \frac{1}{2} \left(\|f({\mathbf{z}})\|_2^2 + \|f({\mathbf{q}}_i)\|_2^2 - \|f({\mathbf{z}}) - f({\mathbf{q}}_i) \|_2^2 \right).\]

  Therefore, we can estimate \(\|f({\mathbf{q}}_i)\|_2^2 = \|{\mathbf{q}}_i\|_p^p\) and \(\|f({\mathbf{z}}) - f({\mathbf{q}}_i)\|_2^2 = \|{\mathbf{z}}- {\mathbf{q}}_i\|_p^p\) using \(F_p\) estimation sketches on \({\mathbf{q}}_i\) and \({\mathbf{z}}- {\mathbf{q}}_i\) respectively, so we obtain an estimate of \(\langle f({\mathbf{z}}), f({\mathbf{q}}_i) \rangle\) up to \(\varepsilon\cdot \left(\|{\mathbf{z}}\|_p^p + \|{\mathbf{q}}_i\|_p^p + \|{\mathbf{z}}-{\mathbf{q}}_i\|_p^p \right)\) additive error. We remark that the estimate of \(\|{\mathbf{z}}- {\mathbf{q}}_i\|_p^p\) precisely corresponds to a previous output of our algorithm.

• By the same reasoning, terms \(\langle f({\mathbf{q}}_i), f({\mathbf{q}}_j) \rangle\) can be estimated using \(F_p\) estimation sketches for previous (fixed) queries \({\mathbf{q}}_i\), \({\mathbf{q}}_j\), and \({\mathbf{q}}_i - {\mathbf{q}}_j\).

As before, the corrector will use an \(F_p\) estimation sketch for \({\mathbf{z}}-{\mathbf{q}}\) to check whether the current estimator is accurate, and to trigger updates to the iterate \({\mathbf{v}}\).

1.2.1.6 Recursion.

At this point, it suffices to show how to recursively estimate \(\|{\mathbf{v}}-f({\mathbf{q}})\|_2^2\) in the next level of the algorithm. However, since the iterate \({\mathbf{v}}\) does not necessarily have a pre-image under \(f\), this significantly complicates the implementation of the recursion. To this end, we let \({\mathbf{q}}= {\mathbf{q}}_L + {\mathbf{q}}_{L+1}\) and expand this quantity as follows: \[\begin{align} \label{eq:recursion:expansion} \|{\mathbf{v}}- f({\mathbf{q}})\|_2^2 &= \left \| \sum_{i = 1}^{L}\alpha_i \cdot f({\mathbf{q}}_i) - f({\mathbf{q}}) \right\|_2^2 \\ &= \sum_{1 \leq i\leq j \leq L} \left( \alpha_i \cdot \alpha_j \cdot \langle f({\mathbf{q}}_i), f({\mathbf{q}}_j) \rangle \right) - 2 \sum_{i= 1}^{L} \alpha_i \cdot \langle f({\mathbf{q}}_i), f({\mathbf{q}}) \rangle + \|f({\mathbf{q}})\|_2^2\nonumber \end{align}\tag{3}\]

We describe how to compute each of the terms above. First, note that since past queries \({\mathbf{q}}_i \in {\mathcal{L}}\) are fixed, inner products of the form \(\langle f({\mathbf{q}}_i), f({\mathbf{q}}_j) \rangle\) can be computed robustly by simply invoking \(F_p\) estimation sketches on \({\mathbf{q}}_i\), \({\mathbf{q}}_j\), and \({\mathbf{q}}_i - {\mathbf{q}}_j\) (as described earlier). Thus, we focus on estimating terms which involve the evolving portion of the stream \({\mathbf{q}}= {\mathbf{q}}_L + {\mathbf{q}}_{L+1}\). In particular, by the properties of the embedding \(f\), the inner products \(\langle f({\mathbf{q}}_i), f({\mathbf{q}}) \rangle\) appearing in the second term can be expressed as

\[\begin{align} \langle f({\mathbf{q}}_i), f({\mathbf{q}}) \rangle &= \frac{1}{2}\left(\|f({\mathbf{q}}_i)\|_2^2 + \|f({\mathbf{q}})\|_2^2- \|f({\mathbf{q}}_i) - f({\mathbf{q}})\|_2^2 \right) \\ &= \frac{1}{2}\left(\|f({\mathbf{q}}_i)\|_2^2 + \|{\mathbf{q}}_L + {\mathbf{q}}_{L+1}\|_p^p - \| {\mathbf{q}}_i - ({\mathbf{q}}_L + {\mathbf{q}}_{L+1})\|_p^p \right) \\ & = \frac{1}{2} \left(\|f({\mathbf{q}}_i)\|_2^2 + \|f({\mathbf{q}}_L) - f(-{\mathbf{q}}_{L+1})\|_2^2 - \|f({\mathbf{q}}_i-{\mathbf{q}}_L) - f({\mathbf{q}}_{L+1})\|_2^2 \right). \end{align}\]

Likewise, the third term can be written as \(\|f({\mathbf{q}})\|_2^2 = \|{\mathbf{q}}_L + {\mathbf{q}}_{L+1}\|_p^p = \|f({\mathbf{q}}_L) - f(-{\mathbf{q}}_{L+1})\|_2^2\). At this point, note that estimating each of the above norms is a problem of precisely the same form as the original problem, and it may seem as though we must then recursively estimate each of these \(\mathcal{O}\left(L^2\right)\) norms arising from inner products in 3 . In fact, it turns out that with some careful algebraic manipulation, we can reduce the problem to only \(\textit{one}\) recursive call on the squared-norm \(\left \|\sum_{i = 1}^L \alpha_i \cdot f({\mathbf{q}}_i - {\mathbf{q}}_L) - f({\mathbf{q}}_{L+1})\right \|_2^2\), and all remaining terms can be robustly computed by directly using \(F_p\) estimation sketches on previous queries \({\mathbf{q}}_i - {\mathbf{q}}_{L}\) or \({\mathbf{q}}_i\) for \(i \in [L]\). ².

1.2.1.7 Bounding the additive error.

Let us examine the total error incurred by our estimates in the approach described above. Recall that for any \({\mathbf{w}}, {\mathbf{y}}\in \mathbb{R}^n\) and any \(\varepsilon' > 0\), each inner product \(\langle f({\mathbf{w}}), f({\mathbf{y}}) \rangle\) can be estimated up to error \(\varepsilon' \cdot \left(\|{\mathbf{w}}\|_p^p + \|{\mathbf{y}}\|_p^p + \|{\mathbf{w}}- {\mathbf{y}}\|_p^p \right)\). Now, suppose that the underlying frequency vector \({\mathbf{x}}= {\mathbf{z}}- {\mathbf{q}}\) satisfies \(\|{\mathbf{x}}\|_p^p \leq 10 A\) at all times for some fixed parameter \(A\) ³. Furthermore, suppose that the coefficients \(\alpha_i\) appearing in 2 and 3 are bounded by \(\mathcal{O}\left(1\right)\). Then, by setting \(\varepsilon' = \frac{\varepsilon}{L^2}\), this means that the error of estimating the first term \(\|f({\mathbf{z}}) - {\mathbf{v}}\|_2^2\) via dot products in 2 is roughly \(L^2 \cdot \mathcal{O}\left(\frac{\varepsilon}{L^2} \cdot A\right) \approx \mathcal{O}\left(\varepsilon\right) \cdot A\); by the same argument, we can also ensure that the error due to the second term \(\|{\mathbf{v}}- f({\mathbf{q}})\|_2^2\) is \(\mathcal{O}\left(\varepsilon\right) \cdot A\) as well.

However, this analysis relied on the assumption that the projection of \(f({\mathbf{z}})\) onto \(\mathop{\mathrm{Span}}({\mathcal{L}})\), given by \({\mathbf{v}}= \sum_{i = 1}^L \alpha_i \cdot f({\mathbf{q}}_i)\), has coefficients \(\alpha_i\) bounded by \(\mathcal{O}\left(1\right)\). In fact, it turns out that this is not necessarily true for arbitrary list vectors \({\mathbf{q}}_1,\ldots, {\mathbf{q}}_{L}\). To overcome this, we instead imagine that the list vectors are actually \({\mathbf{u}}_i=f({\mathbf{q}}_i)+\nu\cdot{\mathbf{e}}_{N+i} \in \mathbb{R}^{N + L}\), and let \({\mathcal{L}}' = \{{\mathbf{u}}_i\}_{i \in [L]}\). By scaling, suppose that \(\|{\mathbf{u}}_i\|_2 = 1\) for all \(i \in [L]\), and let \({\mathbf{U}}\) be the matrix whose columns are \({\mathbf{u}}_1,\ldots, {\mathbf{u}}_{L}\). First, we observe that \(\sigma_{\max}({\mathbf{U}}) \leq \|{\mathbf{U}}\|_F \leq \sqrt{L} \cdot \sqrt{\mathcal{O}\left(1\right) + \nu^2} = \mathcal{O}\left(\sqrt{L}\right)\). Moreover, since the \(L \times L\) identity matrix scaled by \(\nu\) is appended in the dimensions after \(N\), we have that \(\sigma_{\min}({\mathbf{U}}) \geq \nu\). Thus, it follows that the condition number of \({\mathbf{U}}\) is \(\kappa = \frac{\sqrt{L}}{\nu} = \frac{L \cdot H^2}{\varepsilon^2}\), and thus the projection of \(f({\mathbf{z}})\) (with \(0\)’s in additional dimensions) onto \(\mathop{\mathrm{Span}}({\mathbf{u}}_1,\ldots, {\mathbf{u}}_{L})\) will have coefficients satisfying \(|\alpha_i| \leq \kappa = \frac{\sqrt{L}}{\nu}\).

With the above observation in mind, we are now ready to define our modified update rule. Note that finding the projection of \(f({\mathbf{z}})\) onto the space spanned by \({\mathbf{u}}_1, \ldots, {\mathbf{u}}_L\) is equivalent to computing \(\mathop{\mathrm{argmin}}_{{\mathbf{u}}'\in\mathop{\mathrm{Span}}({\mathcal{L}}')}\|f({\mathbf{z}})-{\mathbf{u}}'\|_2^2\). However, since we appended the \(L \times L\) identity scaled by \(\nu\) in the added dimensions, this actually corresponds to solving the regularized regression \(\mathop{\mathrm{argmin}}_{{\mathbf{u}}' \in \mathop{\mathrm{Span}}({\mathcal{L}})} \left(\|f({\mathbf{z}}) - {\mathbf{u}}' \|_2^2 + \nu^2 \cdot \|\mathbf{\alpha}\|_2^2\right)\), where \(\alpha \in \mathbb{R}^L\) denotes the coefficient vector of \({\mathbf{u}}'\) in the basis \({\mathcal{L}}= \{f({\mathbf{q}}_i)\}_{i \in [L]}\). As before, all computations in this regularized regression take place in the sketch space by estimating dot products of the form \(\langle f({\mathbf{z}}), f({\mathbf{q}}_i)\rangle\) or \(\langle f({\mathbf{q}}_i), f({\mathbf{q}}_j)\rangle\) up to additive error \(\mathcal{O}\left(\varepsilon\right) \cdot A\). Finally, if \({\mathbf{u}}' = \sum_{ i =1}^L \alpha_i \cdot {\mathbf{u}}_i\) for \({\mathbf{u}}_i \in {\mathcal{L}}'\), we output the coefficients of \({\mathbf{u}}_1 + {\mathbf{u}}'\), which we use to represent the linear combination of our new iterate \({\mathbf{v}}\).

Overall, this ensures that the coefficients \(\alpha_i\) in the linear combinations of 2 and 3 are bounded, and thus we ensure that the total error accrued in estimating all \(\mathcal{O}\left(L^2\right)\) dot products is \(\mathcal{O}\left(\varepsilon\right) \cdot A\) by setting \(\varepsilon' = \frac{\varepsilon}{\kappa \cdot L^2}\) to be the error parameter for each sketch.

1.2.1.8 From additive error to multiplicative error.

We now describe how to convert a robust additive error algorithm \({\mathcal{A}}\) to a robust multiplicative error algorithm \({\mathcal{A}}'\) by stratifying the estimation process across multiple scales.

Recall that the algorithm above outputs estimates with error \(\varepsilon\cdot A\) when the underlying frequency vector \({\mathbf{x}}\) satisfies \(\|{\mathbf{x}}\|_p^p \leq 10 A\). To obtain multiplicative error, we initialize \(K = \mathcal{O}\left(\frac{1}{\varepsilon}\log n\right)\) independent copies \({\mathcal{A}}_1,\ldots, {\mathcal{A}}_K\) of the robust additive-error algorithm, where \({\mathcal{A}}_i\) achieves error \(\varepsilon\cdot (1+\varepsilon)^i\) when \(\|{\mathbf{x}}\|_p^p \leq 10 \cdot (1+\varepsilon)^i\). If we could maintain the property that input vectors to the robust algorithm \({\mathcal{A}}_i\) for scale \(S_i = [(1+\varepsilon)^i, (1+\varepsilon)^{i+1})\) satisfies \(\|{\mathbf{x}}\|_p^p = \Theta{((1+\varepsilon)^i)}\), then at scale \(S_i\), the additive error of the estimator is \(\varepsilon\cdot \Theta((1+\varepsilon)^i)\), which corresponds to a multiplicative \((1+\varepsilon)\)-approximation for all queries for which the \(F_p\) moment lies in \(\Theta((1+\varepsilon)^i)\). To determine the appropriate scale for a given query \({\mathbf{z}}-{\mathbf{q}}\), we begin at the top level \(K\). The estimator at this level provides an additive approximation of \(\varepsilon\cdot\Theta((1+\varepsilon)^K)\), which suffices to determine whether the true value is smaller than \(2\varepsilon\cdot\Theta((1+\varepsilon)^K)\). If this is the case, we can safely proceed to a lower level and repeat the process. A more fine-grained version of this iterative search allows the algorithm to iteratively zoom in on the correct level, ultimately yielding a multiplicative \((1+\varepsilon)\)-approximation across all scales of \(\|{\mathbf{z}}-{\mathbf{q}}\|_p^p\). So, if each of the algorithms \({\mathcal{A}}_i\) is accurate with high probability, \({\mathcal{A}}'\) returns a multiplicative \((1+\varepsilon)\) approximation.

Finally, it remains to show that this new algorithm \({\mathcal{A}}'\) is robust, i.e. each of the \({\mathcal{A}}_i\)’s continue to be accurate even when updates are adaptive. By contradiction, suppose that there exists a strategy \({\mathcal{S}}\) for the adversary to break a particular \({\mathcal{A}}_i\) with non-trivial probability. In that case, we can construct an adaptive attack on \({\mathcal{A}}_i\) by letting the adversary simulate the execution of the algorithm for the remaining scales \({\mathcal{A}}_1,\ldots, {\mathcal{A}}_{i-1}, {\mathcal{A}}_{i+1}, \ldots, {\mathcal{A}}_K\) by itself, and then execute strategy \({\mathcal{S}}\) to break \({\mathcal{A}}_i\). However, this contradicts the robustness of \({\mathcal{A}}_i\), so we conclude that \({\mathcal{A}}'\) is adversarially robust on streams of length \(m = \mathop{\mathrm{poly}}(n)\).

1.2.2 Sketching, EMD, and Clustering↩︎

1.2.2.1 Equivalence between robustness and sketching.

To formalize the connection between adversarial robustness and classical sketching, we use a well-known result by [51], who established that any finite-dimensional normed space \(\mathcal{X}=(\mathbb{R}^n, \|\cdot\|_{\mathcal{X}})\) admitting an efficient classical linear sketch must also admit a low-distortion linear embedding into \(L_p\) for \(p \in\left(\frac{2}{3},1\right)\) (and into \(L_1\) if the norm is closed under sum-product). Because this embedding is a linear map, we can immediately integrate it within our turnstile streaming framework. Specifically, given a static embedding matrix \({\mathbf{A}}\in \mathbb{R}^{M \times n}\), each stream update to the underlying vector \({\mathbf{x}}\) implicitly updates the coordinates of the corresponding vector \({\mathbf{A}}{\mathbf{x}}\) in the embedded space. We can then execute our adversarially robust \(L_p\) estimation algorithm directly on this virtual stream.

The major issue is that the continuous embedding from [51] is existential and not necessarily computable. To make this algorithmic, we operate in an advice model and prove that the continuous embedding can be discretized. Fortunately, by using standard techniques, it can be shown the target dimension of the virtual stream can be upper bounded by \(M = \mathcal{O}\left(n\log n\right)\), and that the entries of \({\mathbf{A}}\) require only \(\mathcal{O}\left(\log n\right)\) bits of precision. Since the space complexity of our robust \(L_p\) approximation algorithm depends only polylogarithmically on \(M\), we can run it on the transformed stream to establish a weak equivalence between classical sketching and adversarial robustness.

1.2.2.2 Derandomizing embeddings for EMD and clustering.

This approach directly gives robust algorithms for important tasks such as approximating the Earth Mover Distance (EMD) over the grid \([\Delta]^d\). In the classical oblivious setting, EMD is estimated by embedding the spatial distributions into \(L_1\) via randomly shifted quadtrees, achieving \(\mathcal{O}\left(d \log \Delta\right)\) distortion [55]–[57]. However, this randomized embedding \({\mathbf{A}}_s\) only provides a “for each” guarantee, so that an adaptive adversary could observe the algorithm’s outputs, deduce the random shift \({\mathbf{s}}\), and generate updates that force worst-case distortion.

To prevent adversarial attacks against the embedding, we instead completely derandomize the construction. Namely, we define a deterministic embedding \({\mathbf{A}}{\mathbf{x}}= \frac{1}{n} \bigoplus_{s \in [\Delta]^d} {\mathbf{A}}_s{\mathbf{x}}\) that concatenates the scaled linear embeddings induced by all possible \(\Delta^d\) grid shifts simultaneously. By linearity of expectation, the \(L_1\) norm of this concatenated vector exactly matches the expected value over the random shifts, providing a universal “for all” guarantee that holds against any adaptive adversary. While aggregating over all shifts \({\mathbf{s}}\in[\Delta]^d\) blows up the embedding dimension from \(\mathcal{O}\left(n\right)\) to \(M =\mathcal{O}\left(n^2\right)\), where \(n = \Delta^d\), our robust \(F_p\) estimation algorithm again only uses space polylogarithmic in \(M\).

We naturally extend this robust EMD primitive to the \(k\)-median clustering problem, which can be formulated as the Earth Mover Distance between the data distribution \(\mu\) and an assignment vector \(\nu_{{\mathcal{C}}}\) that routes mass to the nearest centers in a candidate set \({\mathcal{C}}\). Using our deterministic \(L_1\) embedding, we can dynamically maintain a robust sketch of the dataset. Because the embedding is deterministic and universally preserves distances, we can safely enumerate over a standard collection of candidate center sets and evaluate their robust costs on the fly, achieving the first robust algorithm for \(k\)-median clustering on turnstile streams.

2 Preliminaries↩︎

2.1 Notation↩︎

Throughout, we write \([n]\) to refer to the set of integers \(\{1,2,\ldots,n\}\) for any positive integer \(n\). The notation \(\mathop{\mathrm{poly}}(n)\) refers to some polynomial in \(n\) with fixed, but unspecified, degree. An event is said to happen with high probability if its likelihood is at least \(1 - \frac{1}{\mathop{\mathrm{poly}}(n)}\). We typically use normal font for scalars and bold font for vectors or matrices.

We typically use normal font for scalars and bold font for vectors or matrices. For vectors \({\mathbf{x}}, {\mathbf{y}}\in \mathbb{R}^n\), we write \(\langle {\mathbf{x}}, {\mathbf{y}}\rangle\) for the standard Euclidean inner product. For \(p > 0\), the \(L_p\) norm (or quasi-norm when \(p < 1\)) of \({\mathbf{x}}\in \mathbb{R}^n\) is defined as \[\|{\mathbf{x}}\|_p = \left(\sum_{i=1}^n |x_i|^p\right)^{1/p}.\] For a vector \({\mathbf{x}}\) and a subspace \(V\), \(\mathop{\mathrm{Proj}}_V({\mathbf{x}})\) denotes the orthogonal projection of \({\mathbf{x}}\) onto \(V\), i.e., the vector in \(V\) closest to \({\mathbf{x}}\) in Euclidean distance. For a set of vectors \(S\), \(\mathop{\mathrm{Span}}(S)\) denotes the linear span of \(S\), i.e., all linear combinations of vectors in \(S\).

We use \(\bigoplus\) to denote vector concatenation (direct sum). For vectors \({\mathbf{x}}^{(1)} \in \mathbb{R}^{n_1}, \ldots, {\mathbf{x}}^{(k)} \in \mathbb{R}^{n_k}\), we define \[\bigoplus_{i=1}^k {\mathbf{x}}^{(i)} = ({\mathbf{x}}^{(1)}, \ldots, {\mathbf{x}}^{(k)}) \in \mathbb{R}^{n_1 + \cdots + n_k}.\]

We use \(\circ\) to denote column-wise concatenation. For vectors \({\mathbf{a}}_1, \ldots, {\mathbf{a}}_k \in \mathbb{R}^n\), the matrix \({\mathbf{a}}_1 \circ \cdots \circ {\mathbf{a}}_k \in \mathbb{R}^{n \times k}\) is defined by taking \({\mathbf{a}}_j\) as its \(j\)-th column. For a matrix \({\mathbf{A}}\in \mathbb{R}^{n \times k}\), we write \({\mathbf{A}}^\dagger\) to denote its Moore-Penrose pseudoinverse, defined as the unique matrix satisfying \[{\mathbf{A}}{\mathbf{A}}^\dagger {\mathbf{A}}= {\mathbf{A}}, \quad {\mathbf{A}}^\dagger {\mathbf{A}}{\mathbf{A}}^\dagger = {\mathbf{A}}^\dagger, \quad ({\mathbf{A}}{\mathbf{A}}^\dagger)^\top = {\mathbf{A}}{\mathbf{A}}^\dagger, \quad ({\mathbf{A}}^\dagger {\mathbf{A}})^\top = {\mathbf{A}}^\dagger {\mathbf{A}}.\] When \({\mathbf{A}}\) has full column rank, this reduces to \({\mathbf{A}}^\dagger = ({\mathbf{A}}^\top {\mathbf{A}})^{-1} {\mathbf{A}}^\top\).

2.2 Preliminaries for Small Moment Estimation↩︎

We recall the existence of an embedding of \(F_p\) into \(F_2\).

Similarly, we describe an embedding for \(F_0\) into \(F_2\) as follows. Given a vector \({\mathbf{x}}\in\mathbb{R}^n\), suppose all coordinates of \(x\) are bounded within \([-m,-(m-1),\ldots,-1,0,1,\ldots,m-1,m]\). Now for each coordinate \(i \in [n]\) and each possible value \(a\) that coordinate \(i\) may take, introduce a distinct orthonormal basis vector \(e_{i,a}\) in a high-dimensional Euclidean space. Then, we define \[g({\mathbf{x}}) = \frac{1}{\sqrt{2}}\sum_{i=1}^n e_{i, x_i},\] so that each coordinate of \({\mathbf{x}}\) is represented by a one-hot vector indicating its value. Now for any two vectors \({\mathbf{x}},{\mathbf{y}}\) whose coordinates are within this range, we then have \[\begin{align} \|g({\mathbf{x}}) - g({\mathbf{y}}) \|_2^2 &= \frac{1}{2}\sum_{i=1}^n \| e_{i, x_i} - e_{i, y_i} \|_2^2\\ &= \frac{1}{2}\cdot 2 \cdot |\{ i : x_i \neq y_i \}|\\ &= \|{\mathbf{x}}- {\mathbf{y}}\|_0. \end{align}\]

However, note that for \(p=0\), we deterministically have \(\|g({\mathbf{x}})\|_2^2=\frac{n}{2}\) for all \({\mathbf{x}}\in\mathbb{R}^n\). To account for this, we further define the embedding by

\[f({\mathbf{x}}) \coloneq g({\mathbf{x}}) - g(\mathbf{0})\]

In particular, observe that this embedding preserves the \(L_0\) norm as before, and maps \(\mathbf{0}\) to \(\mathbf{0}\): \[\|f({\mathbf{x}}) - f({\mathbf{y}})\|_2^2 = \|{\mathbf{x}}- {\mathbf{y}}\|_0 \textrm{ and } \|f(\mathbf{0})\|_2^2 = 0\]

Thus, the one-hot encoding yields an isometric embedding of \((\mathbb{R}^n, \|\cdot\|_0)\) into \((\mathbb{R}^D, \|\cdot\|_2^2)\), where each differing coordinate contributes a fixed amount to the squared Euclidean distance. Hence, we have:

Next, we recall oblivious turnstile streaming algorithms for \(F_p\) moment estimation and distinct element estimation:

None

We now show the correctness of the algorithm \({\mathsf{DotEst}}\), which estimates the dot product of two vectors in the embedded space.

3 Robust \(F_p\) Algorithm for \(p\in[0,2)\)↩︎

In this section, we present our adversarially robust algorithm for \(F_p\) moment estimation on insertion-deletion streams, for \(p\in[0,2)\). Our goal is to generalize the estimator, corrector, and learner framework developed for \(F_2\) moment estimation in [32]. The key challenge is that \(L_p\) spaces for \(p \neq 2\) are no longer Hilbert spaces, i.e., there is no inner product that induces the \(L_p\) norm, so geometric reasoning based on orthogonality or the Pythagorean theorem does not apply. Instead, our analysis relies on the existence of an embedding of \(L_p^p\) into \(L_2^2\).

3.1 Subroutines for Small Moment Estimation↩︎

We next describe an algorithm for \(L_p\) norm estimation, with \(p\in[0,2)\). Utilizing [thm:embed], there exists an embedding \(f\) from \(F_p\) to \(F_2\), so that we have \[\|{\mathbf{x}}-{\mathbf{y}}\|_p^p=\|f({\mathbf{x}})-f({\mathbf{y}})\|_2^2=\|f({\mathbf{x}})\|_2^2-2\langle f({\mathbf{x}}),f({\mathbf{y}})\rangle+\|f({\mathbf{y}})\|_2^2.\]

Note that for all \(p \in [0,2)\) we have \(f(\mathbf{0})=\mathbf{0}\), so that \[\|f({\mathbf{x}})\|_2^2=\|f({\mathbf{x}})-f(\mathbf{0})\|_2^2=\|{\mathbf{x}}-\mathbf{0}\|_p^p=\|{\mathbf{x}}\|_p^p\] and thus \[\|{\mathbf{x}}-{\mathbf{y}}\|_p^p=\|{\mathbf{x}}\|_p^p-2\langle f({\mathbf{x}}),f({\mathbf{y}})\rangle+\|{\mathbf{y}}\|_p^p.\] Then it follows that we can use \(F_p\) estimation algorithms to estimate \(\|{\mathbf{x}}\|_p^p\) and \(\|{\mathbf{y}}\|_p^p\). Moreover, since \(\|{\mathbf{x}}-{\mathbf{y}}\|_p^p=\|f({\mathbf{x}})-f({\mathbf{y}})\|_2^2\), then we can further estimate \(2\langle f({\mathbf{x}}),f({\mathbf{y}})\rangle\) using \(F_p\) estimation algorithms, such as in \({\mathsf{DotEst}}\).

In this subsection, we design and analyze a robust algorithm which estimates the \(F_p\) moment up to \(\varepsilon\cdot A\) additive error when the underlying input vector satisfies \(\|{\mathbf{x}}\|_p^p = \mathcal{O}\left(A\right)\). Then, in 3.2 we show how to modify our \(F_p\) estimation algorithm to obtain the \((1+\varepsilon)\) multiplicative error guarantee.

3.1.0.1 Estimator, corrector, and learner.

As in [32], we have a learner that maintains an approximation \({\mathbf{v}}\) to \(f({\mathbf{z}})\). By analogy, we can update \({\mathbf{v}}\) by adding \(\varepsilon\cdot ({\mathbf{v}}+ f({\mathbf{q}}))\) each time that the estimator is incorrect on a query \({\mathbf{q}}\); however, in this setting, this update rule may result in large coefficients in the linear combination \({\mathbf{v}}= \sum_{i \in [L]} \alpha_i \cdot f({\mathbf{q}}_i)\). Thus, we modify the previous update rule as follows: we add the vector \({\mathbf{q}}\) to a list \({\mathcal{L}}_i\) and solve a kernel regression problem to find the vector \({\mathbf{v}}\) in the span of \({\mathcal{L}}_i\) that best approximates \({\mathbf{z}}\). Crucially, we embed all vectors into \(L_2\) space so that it suffices to solve an \(L_2\) regression problem in the kernel space, where we can compute the closed-form solution using a linear combination of dot products of either the form \(\langle f({\mathbf{z}}),f({\mathbf{q}}_i)\rangle\) or \(\langle f({\mathbf{q}}_i),f({\mathbf{q}}_j)\rangle\), which can be computed using \(\mathcal{O}\left(\frac{1}{\varepsilon^2} \log^3 n\right)\) bits of space by 1.

Formally, the procedure \({\mathsf{EmbedEst}}\) in 4 computes the squared distance \(\|f({\mathbf{z}})-{\mathbf{v}}\|_2^2\) in the kernel space using estimates of the norms and pairwise dot products. The procedure \({\mathsf{EstLevel}}\) in 2 then recursively estimates the \(F_p\) moment at each recursive level \(i\) by combining the embedded approximation from \({\mathsf{EmbedEst}}\) with the estimate from the previous level. Meanwhile, \({\mathsf{MaintainList}}\) in 3 maintains a list \({\mathcal{L}}_i\) at level \(i\), which includes each query vector \({\mathbf{q}}\) for which the corrector flags the output from the estimator as being inaccurate. When this occurs, the subroutine \({\mathsf{GetIterate}}\) in 5 computes \({\mathbf{v}}\) by finding a near-optimal approximation of \(f({\mathbf{z}})\) in the span of the list \({\mathcal{L}}_i\).

3.1.0.2 Recursion.

In the top-level of the algorithm, we use the \(\mathsf{DotEst}\) subroutine to estimate \(\|f({\mathbf{z}}) - {\mathbf{v}}\|_2^2\). Following the robust \(L_2\) algorithm of [32], we need to recursively compute the remaining quantity \(\|{\mathbf{v}}- f({\mathbf{q}})\|_2^2\) in the next level, where the evolving portion of the stream \({\mathbf{q}}= {\mathbf{q}}_L + {\mathbf{q}}_{L+1}\) is divided into blocks which must be sketched independently. Unfortunately, we do not have access to the sketch of \({\mathbf{v}}\) in the embedded space, and in fact, \({\mathbf{v}}\) may not have a pre-image in the original space. This significantly complicates our implementation of the recursion and related computations. In the following lemma, we show that all terms of \(\|{\mathbf{v}}- f({\mathbf{q}})\|_2^2\) can be computed robustly, and that it suffices to recursively compute only one squared-norm.

In order to implement the corrector at each level of the recursion, in 4 we show that the input to the recursion at each level has bounded norm.

3.1.0.3 Update rule and regularized regression.

Next, we discuss our modified update rule for the iterate \({\mathbf{v}}\) in each level of the algorithm. In particular, let \({\mathcal{L}}= \{f({\mathbf{q}}_i)\}_{i \in [L]}\) be the list of queries on which the algorithm was previously inaccurate. Then, we note that the \(L_2\) algorithm of [32] can be implemented by following the alternative update rule \({\mathbf{v}}= \arg \min_{{\mathbf{y}}\in \mathbb{R}^L}{\|{\mathbf{U}}{\mathbf{y}}- f({\mathbf{z}})\|_2^2}\), where \({\mathbf{U}}\) is the matrix with columns \(f({\mathbf{q}}_i)\). Then, by expressing \({\mathbf{v}}= \sum_{i = 1}^L \alpha_i \cdot f({\mathbf{q}}_i)\), we can approximate the estimator \(\|f({\mathbf{z}}) - {\mathbf{v}}\|_2^2 + \|{\mathbf{v}}- f({\mathbf{q}})\|_2^2\) by estimating the necessary inner products of the form \(\langle f({\mathbf{z}}), f({\mathbf{q}}_i) \rangle\) or \(\langle f({\mathbf{q}}_i), f({\mathbf{q}}_j) \rangle\) for \(i, j \in [L]\) via the \(\mathsf{DotEst}\) subroutine. For instance, for the first term of the estimator, we have

\[\|f({\mathbf{z}}) - {\mathbf{v}}\|_2^2 = \|f({\mathbf{z}})\|_2^2 - 2 \sum_{i = 1}^L \alpha_i \cdot \langle f({\mathbf{z}}), f({\mathbf{q}}_i) \rangle + \sum_{1 \leq i,j \leq L} \alpha_i \cdot \alpha_j \langle f({\mathbf{q}}_i), f({\mathbf{q}}_j) \rangle.\]

Similarly, by 3, the second term \(\|{\mathbf{v}}- f({\mathbf{q}})\|_2^2\) can be expanded into inner products involving the previous fixed vectors \(f({\mathbf{q}}_i)\) and \(f({\mathbf{z}})\), and a single squared norm on which we recurse. By 1, each inner product \(\langle f({\mathbf{z}}), f({\mathbf{q}}_i) \rangle\) can be estimated with additive error \(\varepsilon\cdot \left(\|{\mathbf{z}}\|_p^p + \|{\mathbf{q}}_i\|_p^p + \|{\mathbf{z}}- {\mathbf{q}}_i\|_p^p \right)\) by using \(F_p\) estimation sketches, and similarly, \(\langle f({\mathbf{q}}_i), f({\mathbf{q}}_j) \rangle\) can be estimated with additive error \(\varepsilon\cdot \left(\|{\mathbf{q}}_i\|_p^p + \|{\mathbf{q}}_j\|_p^p + \|{\mathbf{q}}_i - {\mathbf{q}}_j\|_p^p \right)\). By assumption, we have that \(\|{\mathbf{z}}\|_p^p = \Theta(A)\), \(\|{\mathbf{z}}- {\mathbf{q}}_i\|_p^p = \Theta(A)\) and all queries \({\mathbf{q}}_i\) also satisfy \(\|{\mathbf{q}}_i\|_p^p \leq \mathcal{O}\left(A\right)\), so each such inner product can be approximated up to \(\varepsilon\cdot A\) additive error. Therefore, in order to upper bound the total error incurred by our estimates via \(\mathsf{DotEst}\), we must ensure that the coefficients \(\alpha_i\) are bounded by \(\mathop{\mathrm{poly}}\left(\frac{1}{\varepsilon}, \log n \right)\). Since this is not true in general, we must modify our update rule, as we describe next.

To this end, we instead consider the vectors \(f({\mathbf{q}}_i)\) with added dimensions, i.e. let \({\mathbf{u}}_i = f({\mathbf{q}}_i) + \nu \cdot {\mathbf{e}}_{N + i} \in \mathbb{R}^{N+L}\) for \(i \in [L]\), and let \({\mathcal{L}}'\) be the list containing \({\mathbf{u}}_1, \ldots, {\mathbf{u}}_L\). Then, we observe that projecting \(f({\mathbf{z}})\) onto \({\mathcal{L}}'\) exactly corresponds to the regularized regression problem \(\mathop{\mathrm{argmin}}_{{\mathbf{u}}' \in \mathop{\mathrm{Span}}({\mathcal{L}})} \left(\|f({\mathbf{z}}) - {\mathbf{u}}' \|_2^2 + \nu^2 \cdot \|\mathbf{\alpha}\|_2^2\right)\), where \(\alpha \in \mathbb{R}^L\) denotes the vector of coefficients of the iterate. Let \({\mathbf{u}}' = \sum_{i=1}^L \alpha_i \cdot {\mathbf{u}}_i\) be the solution to this regularized regression problem; then, for a hard query \({\mathbf{q}}\), we update the current iterate \({\mathbf{v}}\) to \({\mathbf{v}}' = \sum_{ i = 1}^L \alpha_i \cdot f({\mathbf{q}}_i)\).

Finally, recall that in 3 we required that the sum of the coefficients satisfies \(\sum_{ i =1}^L \alpha_i = 1\). To ensure this, we simply impose an additional linear constraint and instead project \(f({\mathbf{z}}) - {\mathbf{u}}_1\) onto the affine subspace \({\mathbf{u}}_1 + \mathop{\mathrm{Span}}\{ {\mathbf{u}}_i - {\mathbf{u}}_1 : i = 2, \ldots,L \}\). The following lemmas allow us to bound the condition number of the vectors \({\mathbf{u}}_1,\ldots, {\mathbf{u}}_{L}\).

3.1.0.4 Algorithm analysis.

Next, we justify the correctness guarantees of our algorithm, assuming that the underlying frequency vector \({\mathbf{x}}\) satisfies \(\|{\mathbf{x}}\|_p^p \leq 10 A\) when the sketch is invoked. In the next lemma, we show that \(P_i\) in \(\mathsf{EstLevel}\) is an accurate approximation for \(\|f({\mathbf{z}}) - {\mathbf{v}}\|_2^2 + \nu^2 \|\alpha\|_2^2\).

Building on the previous lemma, we now argue that the output \(P_i + Q_i\) of \(\mathsf{EstLevel}\)\((i)\) accurately approximates our desired estimator \(S = \|f({\mathbf{z}}) - {\mathbf{v}}\|_2^2 + \|{\mathbf{v}}- f({\mathbf{q}})\|_2^2 + 2 \nu^2 \|\alpha\|_2^2\).

We now argue that \({\mathsf{EstLevel}}(i)\) returns an estimate within additive error \(\frac{\mathcal{O}\left(\varepsilon^2 A\right) \cdot i}{100H^2}\) of the \(F_p\) moment \(\|{\mathbf{z}}- {\mathbf{q}}\|_p^p\).

Next, we show that the estimator in block \(C_{i,j}\) of level \(i \in [H]\) is inaccurate at most \(\mathcal{O}\left(\frac{1}{\eta^2} \log n\right)\) times throughout the stream, where \(\eta = \frac{\varepsilon}{100 H}\) is the total allowable error in each of the \(H\) levels of the recursion. Recall that the estimator at level \(i\) is given by \(\|f({\mathbf{z}}) - {\mathbf{v}}\|_2^2 + \| {\mathbf{v}}- f({\mathbf{q}})\|_2^2 + 2 \nu^2 \|\alpha\|_2^2\), where \({\mathbf{z}}\) represents the previous frequency vector, \({\mathbf{v}}\) is the projection of \(f({\mathbf{z}})\) onto the list of previous queries on which the estimator had erred, and \({\mathbf{q}}\) contains the arriving updates. Since our estimator simulates the robust \(L_2\) algorithm by embedding into \(L_2^2\), progress of our algorithm will also be measured in the \(L_2\) space.

3.2 Achieving Multiplicative Error↩︎

In this section, we show how to convert an algorithm that estimates a function \(F(\cdot)\) up to additive \(\varepsilon\cdot (1+\varepsilon)^i\) error to an algorithm that achieves a multiplicative \((1+\varepsilon)\)-approximation guarantee. This allows us to complete our adversarially robust \(F_p\) moment estimation algorithm on insertion-deletion streams. Let \({\mathcal{A}}_1, \ldots, {\mathcal{A}}_{K}\) denote \(K = \mathcal{O}\left(\frac{1}{\varepsilon}\log n\right)\) independent copies of a robust additive-error algorithm, where \({\mathcal{A}}_i\) achieves additive \(\varepsilon\cdot (1+\varepsilon)^i\) error when the frequency vector \({\mathbf{x}}\) satisfies \(F({\mathbf{x}}) \leq 10 \cdot (1+\varepsilon)^i\).

None

In the next lemma, we show that this algorithm returns a \((1\pm \varepsilon)\) approximation to \(F({\mathbf{x}})\) and is robust to \(m = \mathop{\mathrm{poly}}(n)\) adaptive updates.

We now show correctness for our adversarially robust \(F_p\) estimation algorithm for \(p \in [0,2)\), which achieves a \((1 + \varepsilon)\)-approximation to \(F_p\) using polylogarithmic space.

4 Embedding Applications↩︎

4.1 Embedding into \(L_p\)↩︎

In this section, we give a number of additional applications of our adversarially robust \(L_p\) estimation algorithm for \(p\in[0,2)\). Unlike the applications in 4.3 which embed into \(L_2^2\), the applications in this section embed into \(L_p\) for \(p\in[0,2)\). To that end, recall that a result by [51] guarantees that any norm \({\mathcal{X}}= (\mathbb{R}^n, \|\cdot\|_{{\mathcal{X}}})\) satisfying triangle inequality and sub-multiplicativity of norms admits a classical sketch and embeds linearly into a measure space \(L_p(\mu)\) that solves the following problem:

However, the proof of 4 is existential and not necessarily algorithmically computable in general. Thus to efficiently evaluate this embedding, we utilize an advice model where the streaming algorithm is granted read-only access to a static advice string, i.e., the embedding matrix \({\mathbf{A}}\), that depends solely on the metric space \({\mathcal{X}}\) and dimension \(n\). In particular, we assume that stream updates arrive as discrete increments in the standard computational basis \({\mathbf{e}}_1, \ldots, {\mathbf{e}}_n\) of \(\mathbb{R}^n\) and that the metric space is well-conditioned with respect to this basis; otherwise, the norm of a simple unit stream update could exceed the representable word size. Formally:

Given that the distortion of the embedding of 4 is \(\mathcal{O}\left(sD\right)\), we cannot subsequently embed \(L_p\) into \(L_2^2\), as this could cause a blow-up of \((sD)^{\mathcal{O}\left(H\right)}\) across \(H\) different recursive levels. Instead, we observe that the embedding of 4 is crucially linear, so we can directly map the input stream to a virtual stream over \(\mathbb{R}^M\), where \(M\) is the dimension of the embedding. Our robust \(L_p\) estimation algorithm is then executed over the virtual stream. However, this requires some care about the magnitude of \(M\), as well as the magnitude of the resulting entries, which we now observe.

Thus, we have the following corollary to 4:

Putting together 1 with known equivalences between turnstile algorithms and linear sketches, we have the following:

Finally, we give a number of specific applications of 1 where the linear transformation \({\mathbf{A}}\) is explicit.

4.1.0.1 Robust graph total variation.

In graph signal processing, the graph total variation is a standard tool for quantifying the smoothness of a signal \({\mathbf{x}}\in \mathbb{R}^V\) defined on a graph \(G=(V,E)\). Fixing an arbitrary orientation of the edges, the corresponding incidence matrix \({\mathbf{B}}\in \{-1,0,1\}^{|E| \times |V|}\) maps vertex signals to edge-wise differences. The graph total variation norm can be written as \(\|{\mathbf{x}}\|_{\text{TV}} = \|{\mathbf{B}}{\mathbf{x}}\|_1\), i.e., the \(L_1\) norm of edge differences. The matrix \({\mathbf{B}}\) thus provides an isometric linear embedding from \((\mathbb{R}^V, \|\cdot\|_{\text{TV}})\) into \((\mathbb{R}^E, \|\cdot\|_1)\). Then applying 1, we have:

4.1.0.2 Robust trend filtering (\(k\)-th order total variation).

Consider a temporal signal \({\mathbf{x}}\in \mathbb{R}^n\), which can be viewed as a sequence of observations over time, e.g., measurements or prices. A natural way to quantify the variability of such a signal is through its discrete differences, where the first-order differences capture successive changes \(x_{i+1} - x_i\) and higher-order differences capture more refined notions of variation, such as curvature for \(k=2\) and higher-order fluctuations. Formally, the \(k\)-th order difference operator is represented by a matrix \({\mathbf{D}}^{(k)}\), obtained by iteratively applying a total of \(k\) times the first-order operator \({\mathbf{D}}^{(1)} \in \mathbb{R}^{(n-1)\times n}\), which has \(1\) on the main diagonal and \(-1\) on the superdiagonal, i.e., the diagonal that is directly above and to the right of the main diagonal. In this way, the \(k\)-th order total variation norm is defined as \(\|{\mathbf{x}}\|_{TV^{(k)}} = \|{\mathbf{D}}^{(k)}{\mathbf{x}}\|_1\), so that it is the sum of the absolute values of these higher-order differences. Intuitively, this quantity measures how irregular the signal is at scale \(k\). For example, \(k=1\) penalizes sharp jumps, while \(k=2\) penalizes changes in slope, encouraging piecewise linear structure. This notion forms the basis of trend filtering, where one seeks signals that balance fidelity to data with low higher-order variation. Note that since \({\mathbf{D}}^{(k)}\) can be explicitly computed, the \(k\)-th order total variation norm embeds explicitly and isometrically into \(L_1\). Thus by applying 1, we have:

4.1.0.3 Robust Earth Mover distance (EMD).

Earth Mover distance (EMD) is a common metric for comparing two distributions over a discrete grid \([\Delta]^d\), so that \(n = \Delta^d\). Intuitively, EMD quantifies the minimal “work” needed to transform one distribution into another, where work is the amount of mass moved times the distance the mass is moved. We interpret \(\|{\mathbf{x}}\|_{\mathop{\mathrm{EMD}}}\) as the Earth Mover Distance between the signed measure \({\mathbf{x}}\) and the zero vector \(\mathbf{0}\), i.e., \(\|{\mathbf{x}}\|_{\mathop{\mathrm{EMD}}} := \mathop{\mathrm{EMD}}({\mathbf{x}}, \mathbf{0})\). Equivalently, for two distributions \(\mu,\nu\), we have \(\mathop{\mathrm{EMD}}(\mu,\nu) = \|\mu - \nu\|_{\mathop{\mathrm{EMD}}}\). Perhaps the most well-used embeddings for EMD are based on randomly shifted quadtrees and give \(\mathcal{O}\left(d\log\Delta\right)\) distortion, but provide a “for each” guarantee. That is, for a fixed signal \({\mathbf{x}}\), the embedding preserves distances in expectation over its random choices. However, this randomized embedding is insufficient for our purposes in the adversarially robust streaming setting, as an adaptive adversary can exploit the particular random shift to create worst-case distortions.

Instead, consider the following construction of a fully deterministic embedding. Let \({\mathbf{A}}_{{\mathbf{s}}} : \mathbb{R}^n \to \mathbb{R}^K\) denote the linear embedding induced by a single quadtree shift \({\mathbf{s}}\), with \(K = \mathcal{O}\left(n\right)\). We define the deterministic, concatenated embedding \[{\mathbf{A}}{\mathbf{x}}= \frac{1}{n} \bigoplus_{s \in [\Delta]^d} {\mathbf{A}}_{{\mathbf{s}}}{\mathbf{x}},\] which aggregates the contributions from all possible shifts. By linearity of expectation, the \(L_1\) norm of \({\mathbf{A}}{\mathbf{x}}\) exactly matches the expected value over all shifts: \[\|{\mathbf{A}}{\mathbf{x}}\|_1 = \frac{1}{n} \sum_{{\mathbf{s}}\in [\Delta]^d} \|{\mathbf{A}}_{{\mathbf{s}}}{\mathbf{x}}\|_1 = \underset{{\mathbf{s}}}{\mathbb{E}}\left[\|{\mathbf{A}}_{{\mathbf{s}}}{\mathbf{x}}\|_1\right].\] This gives a deterministic, universal “for all” guarantee, so that for all \({\mathbf{x}}\in\mathbb{R}^n\), we have \[\|{\mathbf{x}}\|_{\mathop{\mathrm{EMD}}} \le \|{\mathbf{A}}{\mathbf{x}}\|_1 \le \mathcal{O}\left(d \log \Delta\right)\|{\mathbf{x}}\|_{\mathop{\mathrm{EMD}}}.\] Although the embedding dimension grows from \(K=\mathcal{O}\left(n\right)\) to \(M=nK=\mathcal{O}\left(n^2\right)\), the resulting \(L_1\) estimation algorithm only uses space polylogarithmic in the dimension. Hence by applying 1, we have:

4.1.0.4 Robust \(k\)-median clustering.

Next, we consider the \(k\)-median clustering problem over the grid \([\Delta]^d\), where a dataset is represented as a multiset of points, or equivalently, a mass vector \(\mu \in \mathbb{R}^{[\Delta]^d}\). The \(k\)-median objective seeks a set of \(k\) centers \({\mathcal{C}}\subseteq [\Delta]^d\) minimizing the total \(L_2\) distance from each point to its nearest center. Equivalently, for an assignment vector \(\nu_{{\mathcal{C}}}\) that routes each unit of mass in \(\mu\) to its closest center in \({\mathcal{C}}\), the \(k\)-median cost can be expressed as an Earth Mover Distance: \[\mathop{\mathrm{cost}}_k(\mu) = \min_{{\mathcal{C}}: |{\mathcal{C}}|=k} \mathop{\mathrm{EMD}}(\mu, \nu_{{\mathcal{C}}}).\] Thus, to approximate this quantity in the streaming setting, previous works often utilize the standard embedding of \(\mathop{\mathrm{EMD}}\) into \(L_1\) via randomly shifted quadtrees. Fix a shift \(s \in \mathbb{Z}^d\) and suppose without loss of generality that \(\Delta=2^\ell\) for some integer \(\ell>0\). Then for each scale \(t \in \{0,1,\ldots,\ell\}\), define a grid \(\mathcal{G}_{s,t}\) with side length \(2^t\). For a mass vector \(\mu\), let \(G_{s,t}\mu\) denote the vector recording the total mass in each cell of \(\mathcal{G}_{s,t}\). Then the embedding is defined by concatenation across scales: \[G_s \mu = (G_{s,0}\mu)\circ (2 \cdot G_{s,1}\mu)\circ \cdots \circ (2^t \cdot G_{s,t}\mu)\circ \cdots \circ (2^\ell \cdot G_{s,\ell}\mu).\] The key intuition is that the contribution of a unit of mass to the \(L_1\) norm \(\|G_s(\mu - \nu)\|_1\) is determined by the finest scale at which the mass assigned in \(\mu\) and \(\nu\) fall into different cells. Because coarser grids are weighted more heavily, this construction ensures that the \(L_1\) norm faithfully captures transportation cost. In particular, this embedding provides an upper bound on \(\mathop{\mathrm{EMD}}(\mu,\nu)\) up to a factor of \(\mathcal{O}\left(d \log \Delta\right)\).

As in the EMD setting, to obtain adversarial robustness we derandomize the construction by concatenating over all shifts, yielding a deterministic embedding \({\mathbf{A}}\) with a “For All” guarantee. Applying 1 to the embedded vector allows us to approximate \(\|{\mathbf{A}}(\mu - \nu_{{\mathcal{C}}})\|_1\) for all candidate center sets \({\mathcal{C}}\) simultaneously in sublinear space.

Finally, we recall that there exists a collection \({\mathcal{C}}\) of candidate sets of centers and assigned masses such that \(|{\mathcal{C}}|\le(\Delta m)^{\mathcal{O}\left(kd\right)}\), one of which can be used to find a constant-factor approximation to the optimal \(k\)-median clustering cost [56], [57]. For our purposes, it suffices to enumerate over this candidate set, though we remark that there are more efficient top-down search processes to find a good candidate set [56], [57]. Therefore, we have:

4.2 Applications to Entropy Estimation↩︎

In this section, we leverage our adversarially robust \(F_p\) estimation algorithms to obtain a black-box procedure for estimating entropy. Recall that given a frequency vector \({\mathbf{v}}\in\mathbb{R}^n\), the Shannon entropy is defined as \(H({\mathbf{v}})=-\sum_{i=1}^n v_i\log v_i\). We first recall the following relationship between the Shannon entropy and its exponentiation:

In light of this observation, it suffices to focus on computing a \((1+\varepsilon)\)-multiplicative approximation to \(h(v)=2^{H(v)}\). To that end, we recall to following structural property that facilitates entropy approximation.

The points \(\{y_0,\ldots,y_k\}\) from 12 can be constructed explicitly as described in [61]. Let \(\ell=\frac{1}{2(k+1)\log m}\) and define \[f(z)=\frac{(k^2\ell)z-\ell(k^2+1)}{2k^2+1}.\] Then each \(y_i\) is given by \(y_i=1+f\left(\cos\left(\frac{i\pi}{k}\right)\right)\), which allows the entire set to be computed in linear time. Finally, a \((1+\varepsilon)\)-multiplicative estimate of \(h(v)=2^{H(v)}\) is obtained from \(2^{P(0)}\), where \(P(x)\) is the degree-\(k\) polynomial interpolating the points \(\{(y_i, F_{y_i}(v))\}_{i=0}^k\). Crucially, since \(y_i\in(0,2)\) for each \(i\), we can apply our \(F_p\) estimation algorithm to get an additive \(\varepsilon\)-approximation to the Shannon entropy:

4.3 Embedding into \(L_2^2\)↩︎

In this section, we show that our algorithm for \(L_p\) estimation, \(p \in [0,2)\) can be extended to a broader class of functions. In particular, consider the task of estimating \(\left\lVert{\mathbf{x}}\right\rVert_g := \sum_i g(x_i)\) for some specific function \(g\). We show that if \(g\) has the form \(g(t) = f(t^2)\) for some Bernstein function \(f\) (4), then we can design an adversarially robust streaming algorithm which computes a \((1 \pm \varepsilon)\)-approximation to \(\left\lVert{\mathbf{x}}\right\rVert_g\) at each time in the adaptive stream, using \(\mathop{\mathrm{poly}}\left(\frac{1}{\varepsilon},\log n\right)\) bits of space. We note that this definition captures many interesting choices of the function \(g\), which are relevant for various applications.

• (Pseudo–Huber Loss): \(g_{\tau}(x) = \tau(\sqrt{1 + (x / \tau)^2} - 1)\)

• (Cauchy/Lorentzian Loss): \(g_{\tau}(x) = \log (1 + x^2 / \tau)\)

• (Generalized Charbonnier): \(g_{\tau} = (1 + x^2 / \tau)^\beta - 1\) for \(0 < \beta \le 1\).

• (Welsch/Leclerc Loss): \(g_{\tau}(x) = 1 - e^{-x^2 / \tau}\)

• (Geman–McClure loss): \(\displaystyle g_{\tau} =\frac{x^2}{x^2 + \tau}\)

We first recall the following classical result in [52].

By applying 8 and 9 coordinate-wise and then taking the sum over all coordinates, we have the following theorem. Without loss of generality, we can assume that \(\Phi(0) = 0\); otherwise, we can simply replace the embedding \(\Phi({\mathbf{x}})\) by \(\Phi({\mathbf{x}}) - \Phi(0)\) to satisfy this condition.

To utilize our algorithm for \(L_p\) estimation for \(p \in[0, 2)\), we also need a (non-robust) algorithm for estimating \(\left\lVert{\mathbf{x}}\right\rVert_g\). In the work of [32], the authors give the following lemma, which is based on the zero-one law of the form of the function \(g\) in [64].

Finally, by combining 10 and 13, we obtain adversarially robust turnstile streaming algorithms that compute a \((1+\varepsilon)\)-approximation for any \(g\)-sum \(\|{\mathbf{x}}\|_g\), where \(g(x) = h(x^2)\) for some Bernstein function \(h\) and \(h(0) = 0\).

Acknowledgements↩︎

Elena Gribelyuk and Huacheng Yu are supported in part by an NSF CAREER award CCF-2339942. Honghao Lin was supported in part by a Simons Investigator Award and a CMU Paul and James Wang Sercomm Presidential Graduate Fellowship. David P. Woodruff is supported in part by Office of Naval Research award number N000142112647 and a Simons Investigator Award. Samson Zhou is supported in part by NSF CCF-2335411 and gratefully acknowledges funding provided by the Oak Ridge Associated Universities (ORAU) Ralph E. Powe Junior Faculty Enhancement Award.

References↩︎

---

1. As we explain next, \({\mathbf{v}}\) can be computed via a kernel regression over these previously observed queries \(f({\mathbf{q}}_1),f({\mathbf{q}}_2),\ldots,f({\mathbf{q}}_L)\) on which the estimator failed. Crucially, this regression can be implemented using only estimates of norms and pairwise inner products of the form \(\langle f({\mathbf{q}}_i), f({\mathbf{q}}_j)\rangle\) and \(\langle f({\mathbf{z}}), f({\mathbf{q}}_i)\rangle\), which can themselves be obtained using \(F_p\) sketches.↩︎

2. We refer the reader to 3 for more details regarding how to estimate various terms in the expansion.↩︎

3. Note that for this condition, we require that any any time, the prefix of the stream has \(F_p\) moment bounded by \(A\). Specifically, \(\|{\mathbf{z}}\|_p^p \leq \mathcal{O}\left(A\right)\) and \(\|{\mathbf{z}}- {\mathbf{q}}_i\|_p^p \leq \mathcal{O}\left(A\right)\) for all previous queries \({\mathbf{q}}_i\). Also, this directly implies that \(\|{\mathbf{q}}_i\|_p^p \leq \mathcal{O}\left(A\right)\), so \(\|{\mathbf{q}}_i - {\mathbf{q}}_j\|_p^p \leq \mathcal{O}\left(A\right)\) as well.↩︎