The streaming model is an abstraction of computing over massive data streams, which is a popular way of dealing with large-scale modern data analysis.
In this model, there is a stream of data points, one after the other.
A streaming algorithm is only allowed one pass over the data stream, and the goal is to perform some analysis during the stream while using as small space as possible.

Clustering problems (such as $k$-means and $k$-median) are fundamental unsupervised machine learning primitives, and
streaming clustering algorithms have been extensively studied in the past.
However, since data privacy becomes a central concern in many real-world applications, non-private clustering algorithms are not applicable in many scenarios.

In this work, we provide the first differentially private streaming algorithms for $k$-means and $k$-median clustering of $d$-dimensional Euclidean data points over a stream with length at most $T$ using $\poly(k,d,\log(T))$ space to achieve a constant multiplicative error and a $\poly(k,d,\log(T))$ additive error.
In particular, we present a differentially private streaming clustering framework which only requires an offline DP coreset or clustering algorithm as a blackbox.
By plugging in existing results from DP clustering [29, 42], we achieve (1) a $(1+\gamma)$-multiplicative approximation with $\tilde{O}_\gamma(\poly(k,d,\log(T)))$ space for any $\gamma>0$, and the additive error is $\poly(k,d,\log(T))$ or (2) an $O(1)$-multiplicative approximation with $\tilde{O}(k^{1.5} \cdot \poly(d,\log(T)))$ space and $\poly(k,d,\log(T))$ additive error.

In addition, our algorithmic framework is also differentially private under the continual release setting, i.e., the union of outputs of our algorithms at every timestamp is always differentially private.

§ INTRODUCTION

In real-world applications, a major challenge in dealing with large-scale data is that the entire datasets are too large to be stored in the computing system.
The need to address this challenge and the success of large-scale systems (such as Spark Streaming [44]) that process data in streams have driven the study of the streaming model, which is introduced by a seminal work [3].
In this model, there is a stream of data points.
During the stream, a data point arrives at each timestamp, and a streaming algorithm can only access these data points by a single pass.
The goal is to output an (approximate) solution to a problem with respect to the set of data points that are arrived while using as small space as possible .
If an algorithm is required to output at every timestamp when a new data point arrives, then it is called the continual release setting.
Otherwise, the algorithm only needs to output at the end of the stream, which is called non-continual release or the one-shot setting.
The streaming model has attracted a lot of attention from different areas in the past decades.
In particular, streaming clustering algorithms have been extensively studied by the clustering literature.

Clustering is an essential primitive in unsupervised machine learning, and its geometric formulations, such as $k$-means and $k$-median, have been studied extensively, e.g., [4, 11, 32, 12, 13, 5, 40, 36, 2].
Euclidean $k$-clustering problem is stated as follows.
Given a point set $P\subset \mathbb{R}^d$, the goal is to output a set of centers $C\subset\mathbb{R}^d$ with size $|C|=k$ such that the $k$-clustering cost $\textsf{cost}(C,P)=\sum_{p\in P}\min_{c\in C}\|p-c\|_2^z$ is minimized, where $z=1$ and $z=2$ correspond to $k$-median and $k$-means problem respectively.
There is a long list of work (e.g., a subset includes [32, 31, 14, 26, 8, 18, 20]) studying $k$-means and $k$-median problem in the streaming setting.
In the streaming $k$-clustering problem, there is a stream of points in $\mathbb{R}^d$, and the goal is to output a set of $k$ centers at each timestamp $t$, and minimize the $k$-clustering cost with respect to data points arrived before $t$.
The state-of-the-art result is achieved by [20] which uses $\tilde{O}\left(\frac{kd}{\gamma^2}\right)\cdot \min\left(\frac{1}{\gamma^z},k\right)\cdot \poly(\log\log T)$ space to obtain a $(1+\gamma)$-approximation with probability $0.9$ at the end of the stream.
However, none of these algorithms are private, which means they are not applicable when the dataset involves personal information and privacy is considered to be a major concern in real-world applications.

Differential privacy (DP) [21] has become the de facto standard for preserving data privacy due to its compelling privacy guarantees and mathematically rigorous definition.
DP $k$-clustering algorithms in the offline setting have been studied for years [38, 25, 28, 30, 6, 34, 39, 42, 29, 15], where the main focus is to improve the approximation ratio and achieve fast sequential running time.
DP $k$-clustering problem has also been studied in other computational models which are more relevant to large-scale computations such as sublinear-time [7] and massively parallel computing (MPC) [16, 17], and distributed computing setting [43, 10].
However, the landscape of the DP $k$-clustering problem in the streaming model is still mysterious.
In fact, to the best of our knowledge, there is no previously known DP streaming algorithm achieving $O(1)$-multiplicative error using sublinear space even in the one-shot setting.

In this work, we present the first DP streaming algorithms for Euclidean $k$-means and $k$-median clustering using $\poly(k,d,\log(T))$ space to achieve an $O(1)$-multiplicative error and a $\poly(k,d,\log(T))$-additive error.
In addition, our algorithms are DP under the continual release setting, which means that the union of all historical outputs of our algorithm is DP.
Note that any DP algorithm under the continual release setting is always DP under the one-shot setting.

§.§ Computation Model and Differential Privacy

In this section, we formally define differential privacy and the streaming model under continual release setting. The input is a stream of points $x_1,x_2,\cdots,x_T\in\mathbb{R}^d$, where each $x_i\in\mathbb{R}^d$ satisfies $\|x_i\|_2\leq \Lambda$, i.e., we assume all input points are within a ball of radius $\Lambda$.
Streams $\cS=(x_1,\ldots, x_T)$ and $\cS'=(x'_1,\ldots, x'_T)$ are neighboring if there exists at most one timestamp $t^* \in [T]$ for which $x_{t^*} \neq x'_{t^*}$ and $x_{t} = x'_{t}$ for all $t \neq t^*$.
In this paper, we consider every streaming algorithm $\cA$ under the continual release setting, i.e., the entire output of $\cA$ is $(s_1,s_2,\cdots,s_T)$ where $s_t$ is the output of $\cA$ at timestamp $t$ with respect to the data arrived no later than $t$.
If we do not specify the timestamp, then the output of $\cA$ indicates the entire output of $\cA$ over the stream.

A randomized algorithm $\cA$ is $(\eps,\delta)$-DP if for every pair of neighboring datasets $D \sim D'$, and for all sets $\cS$ of possible outputs, we have that $\Pr[\cA(D) \in \cS ] \leq e^\eps \Pr [\cA(D') \in \cS] +\delta  $. When $\delta=0$ we simply say that the algorithm is $\eps$-DP.

§.§ Our Results

We present a general framework for DP $k$-clustering in the streaming setting which utilizes an offline DP semicoreset (see semicoreset) algorithm as a black-box. Using existing results from the DP clustering literature, the cost of the resulting clustering output by our framework achieves (1) a $(1+\gamma)$-multiplicative error with space complexity having a $\poly(k)$ dependency — using the DP coreset algorithm from [29]  (see inf-final-ghazi), or (2) an $O(1)$-multiplicative error with space complexity having a $k^{1.5}$ dependency  — using the DP clustering algorithm from [42] (see inf-final-stemmer). We use the term semicoreset as a generalization for any coreset that has an additive error proportional to the optimal cost. Notably any DP coreset is also a DP semicoreset, and [17] show how to obtain a DP semicoreset from an offline DP $k$-clustering algorithm which we will use as a blackbox in the sequel to achieve better space bounds in the case of inf-final-stemmer.

We state our results for $k$-means in the sequel, but note that our results easily generalize to $k$-median. In the following results, we assume we are given a non-DP algorithm in the offline setting that can compute a $(1+\gamma)$-approximation to $k$-means (e.g. [41]). As is standard in DP clustering literature, we assume $\Lambda$ is an upper bound on the diameter of the space of input points.

Given dimension $d$, clustering parameter $k$, $(1+\gamma)$-approx non-DP coreset algorithm $(\gamma\in(0,0.5))$ [18], $\eps$-DP coreset algorithm from [29], and $(1+\gamma)$-approx non-DP $k$-means clustering algorithm [41]. Let $\cS$ be the stream of input points in Euclidean space  with length at most $T$. There exists a streaming algorithm $\cA'$ for $k$-means that outputs a set of $k$ centers $\cC_{\hat{\cY}}$ at every timestep $t \in [T]$ such that

* (Privacy) $\cA'$ is $3\eps$-DP under the continual release setting.

* (Accuracy) With high probability,

\begin{align}
\mathsf{cost}(\cC_{\hat{\cY}},\cS) &\leq (1+\gamma')\mathsf{cost}(\cC^{opt}_\cS, \cS) %\nonumber\\&
+ O_{\gamma}(1)\cdot V'(d,k,\eps,T,\Lambda) %+ O_{\gamma,\rho}(1) \cdot V''(d,k,\eps,T,\Lambda,\xi)
\end{align}

where $\cS$ denotes the set of all points given by the stream before timestamp $t$, $\cC^{opt}_\cS$ is the optimal solution, $\gamma' = 2\gamma + \gamma^2 + \frac{2(1+\gamma)^5}{(1-\gamma)^3} =O(\gamma) $, and  $V'(d,k,\eps,T,\Lambda) = \tilde{O}( \Lambda^2\cdot\poly( k, d,\frac{1}{\eps},\log(T)))$.

* (Space) $\cA'$ consumes $\tilde{O}_\gamma \left(2^{O_{\gamma}(d')}\cdot\poly(k,d,\frac{1}{\varepsilon},\log(T))\right)$

space where $d'=O_\gamma(\log k)$.

Given dimension $d$, clustering parameter $k$, $(1+\gamma)$-approx non-DP coreset algorithm $(\gamma\in(0,0.5))$ [18], $(\eps,\delta)$-DP k-means clustering algorithm from [42], and $(1+\gamma)$-approx non-DP $k$-means clustering algorithm [41]. Let $\cS$ be a stream of input points in Euclidean space with length at most $T$. There exists a streaming algorithm $\cA'$ for $k$-means that outputs a set of centers $\cC_{\hat{\cY}}$ at every timestep $t \in [T]$ such that

* (Privacy) $\cA'$ is $(5\eps,\delta)$-DP under the continual release setting.

* (Accuracy) With high probability,

\begin{align}
\mathsf{cost}(\cC_{\hat{\cY}},\cS)
&\leq O_{\gamma}(1) \cdot \mathsf{cost}(\cC^{opt}_\cS, \cS)+O_{\gamma}(1) \cdot V'(d,k,\eps,\delta, T,\Lambda)
\end{align}

where $\cS$ denotes the set of all points given by the stream before timestamp $t$, $\cC_{\cS}^{opt}$ is the optimal solution, and $V'(d,k,\eps,\delta,T,\Lambda) = \tilde{O}( \Lambda^2\cdot\poly(\log(T),\log(1/\delta),d,\frac{1}{\eps},k))$.

* (Space) $\cA'$ consumes $\tilde{O}(k^{1.5} \cdot \poly(\log(T),\log(\frac{1}{\delta}),d,\frac{1}{\eps}))$.

The notation $O_x(\cdot)$ ignores factors involving $x$.

The space complexity achieved by our framework depends on the additive error of the offline DP semicoreset (see Definition <ref> for the definition of semicoreset) algorithm used as a black-box. In the case of inf-final-stemmer, the DP $k$-means algorithm from [42] achieves an additive error with $k^{1.5}$ dependency whereas their DP coreset algorithm has $\poly(k)$ dependency. Thus we use their DP $k$-means algorithm to compute a DP semicoreset albeit with a slight loss in privacy to save on space (see clustering-to-semicoreset for the transformation). In the case of inf-final-ghazi, the additive error of the DP $k$-means clustering algorithm and DP coreset from [29] has the same dependency in $k$. Hence we use their DP coreset algorithm directly.

Finally, note that if there exists a DP semicoreset algorithm which achieves a better additive error in the offline setting, then our framework immediately gives a DP clustering algorithm with better space bounds in the continual release setting.

§.§ Comparison with Concurrent Work [35]

[35] concurrently released results on differentially private $k$-means clustering under continual observation.
However their work is mostly incomparable with ours: their main focus is to optimize the approximation ratio under the setting of continual observation, but their algorithm uses a large space.
In contrast, our main focus is to optimize the approximation ratio using as small space as possible, and the property that our algorithm is DP under continual observation (or equivalently, continual release) setting is a side outcome of our algorithm.

§.§ Technical Overview

Our techniques apply to both $k$-means and $k$-median clustering, but we assume we are dealing with $k$-means for simplicity. Before discussing our algorithm in more detail, we first outline the challenges to designing a DP $k$-means clustering algorithm in the continual release setting.

Vanilla Merge and Reduce approaches fail.
A natural idea is to follow a standard approach called Merge-and-Reduce [27] in the streaming clustering literature.
A coreset of a set of points is a weighted subset of the points such that the clustering cost of the coreset with respect to any $k$-centers $C$ is a good approximation of the clustering cost of the original point set with respect to $C$.
We refer readers to Definition <ref> for a formal definition of coreset.
According to [18], any point set has a $(1+\gamma)$-coreset with size $k\cdot \poly(1/\gamma)$.
Merge-and-Reduce approach suggests the following.
We maintain $C_1,C_2,\cdots,C_L$ for $L=\log(T)$ during the stream such that $C_1\cup C_2\cup\cdots\cup C_L$ is a $(1+\gamma)$-coreset of the points already arrived.
We choose some threshold $M=k\poly(\log(T))$.
When a new point $x$ comes, we make $C_1\gets C_1\cup \{x\}$.
Then we do the following check: if there is $i\in[L]$ such that $|C_i|\geq M$, we compute a $(1+1/\poly(\log(T)))$-coreset $C'$ of $C_i$, and we set $C_i\gets \emptyset$ and $C_{i+1}\gets C_{i+1}\cup C'$.
We repeat the above checking process until every $C_i$ has size $|C_i|\leq M$.
It is clear that at any time, $C_i$ is a coreset of a subset $S_i$ of input points and $S_1\cup S_2 \cup\cdots\cup S_L$ is the entire set of the input points.
By induction, we can show that $C_i$ must be a $(1+1/\poly(\log(T)))^i$-coreset of $S_i$, and thus $C_1\cup C_2\cup \cdots \cup C_L$ is a $(1+1/\poly(\log(T)))^L=(1+O(\log(T))/\poly(\log(T)))$-coreset of the entire input point set.

However, the above approach fails when using DP coreset for $C_1$.
The reason is that any DP coreset must introduce some additive error $\eta=\Omega(\Lambda^z)$ (see e.g., [29]).
In particular, let us follow the above Merge and Reduce process, but when $C_1\geq M$, we compute a DP coreset $C'$ of $C_1$ instead.
Then $C'$ is merged to $C_2$, and thus $C_2$ is a coreset of $S_2$ with some additive error $\eta$.
Now consider $\hat{C}$ is a corest of $\hat{S}$ with additive error $\hat{\eta}$ and  $\bar{C}$ is a corest of $\bar{S}$ with additive error $\bar{\eta}$.
Then even if we compute any non-DP coreset $\tilde{C}$ of $\hat{C}\cup\bar{C}$ with no relative error, $\tilde{C}$ is a coreset of $\hat{S}\cup\bar{S}$ with additive error $\bar{\eta}+\hat{\eta}$.
The above example implies that the additive error can actually accumulate.
Therefore, $C_i$ is a coreset of $S_i$ with additive error $2^{O(i)}\cdot \Omega(\Lambda^z)$ which can be as large as $T^{\Omega(1)}\cdot \Lambda^z$.

To overcome the above issue, we develop a novel method of DP Merge and Reduce: instead of running Merge and Reduce for all points together, we firstly partition the space $\mathbb{R}^d$ such that close points are likely to be in the same group, then we run Merge and Reduce in a DP way for points in each group.
In this way, we are able to show that the additive error introduced can always be bounded by the optimal clustering cost with some small multiplicative factors.
In the remainder of this section, we briefly discuss the high level ideas of our approach.

Our Approach.
For every timestamp $t \in [T]$, our algorithm
can be split into the following main steps — (1) Compute a set of centers $\cF$ in an online DP fashion that satisfies a bicriteria approximation[a bicriteria approximation algorithm in this case means that the algorithm may output a set of more than $k$ centers] to $k$-means. (2) Maintain DP semicoresets of the points assigned to centers in $\cF$ in parallel via multiple DP Merge-and-Reduce instances (3) Output the union of these semicoresets called $\hat{\cY}$. In a post-processing step — Compute a non-DP $k$-means $\rho$-approximation algorithm on $\hat{\cY}$.

Bicriteria Approximation.  The bicriteria approximation uses two main ingredients — quadtrees and heavy hitters. A quadtree creates a nested series of grids that partitions $\bbR^d$ and can be used to embed input points into a Hierarchically Separated Tree (HST) metric, which often makes computing $k$-means cost simpler. We use this embedding to map every input point to the center of a grid (or cell) at every quadtree level. For a fixed level, our goal is to approximately choose the $O(k)$ cells that have the most points, i.e., we want to find the “heaviest” cells in a DP fashion and store them as candidate centers in set $\cF$. We achieve this by hashing the cells into $O(k)$ substreams and running a continual release black-box DP heavy hitter algorithm from [24] on each hash substream. Since with large enough probability, the heavy cells will not collide, this achieves our goal. Note that since we need to do this over logarithmically many levels of the quadtree, we will end up with a bicriteria approximation.

There exists an $\eps$-DP algorithm that for every timestamp $t \in [T]$, computes a weighted set of $O({k}\log^2(k) \log(\Lambda) \log T)$ centers with $d^{O(1)}$-multiplicative error to the best $k$-means (or $k$-median) clustering and $\tilde{O}(\frac{k\rho\Lambda^2}{\eps}\cdot( d\log T)^{O(1)}))$-additive error in  $O(k \log^2\Lambda \log(k) \poly\left( \log \left({T}{\Lambda}k\right)\right))$ space.

Reducing Multiplicative Error. At this point, one idea is to assign input points to candidate centers in $\cF$ obtained from the bicriteria approximation in an online fashion and release a weighted DP semicoreset at every timestep. The postprocessing step of applying an offline non-DP clustering algorithm on this weighted DP coreset would ensure a DP set of $k$ centers, however the cost of clustering of the resulting set of centers incur a multiplicative error of $d^{O(1)}$ (see bicriteria), which can be quite large especially in high dimensional Euclidean space.

Inspired by [17] who recently used techniques introduced by [14] to obtain a near optimal DP Clustering algorithm in the MPC setting, we define rings $R_r$ ( where $1\leq r \leq \log \Lambda$) based on distance from the current centers in set $\cF$ and map input points to these rings. The main idea is to compute DP semicoresets for points in each ring and then by taking the union of the semicoresets over all rings, we can obtain a DP semicoreset for the stream of points seen so far. Since the DP semicoreset algorithms we use [29, 42] achieve constant multiplicative error, intuitively, this technique will result in constant multiplicative error of the resulting DP semicoreset. In particular, to compute DP semicoresets for points in each ring we run Merge and Reduce in a DP manner.

DP Merge and Reduce. We design an online differentially private variant of the celebrated Merge and Reduce framework [33, 1] to compute DP semicoresets for points in each ring  $R_r$ ( where $1\leq r \leq \log \Lambda$). An informal overview of the (non-DP) Merge and Reduce technique is — partition the input stream into blocks, compute a coreset for each block, take the union of the resulting coresets (merge step), and compute the coreset of the union (reduce step) in a tree-like fashion.
On a high-level, our framework computes DP semicoresets at the base level (of the tree) using a DP semicoreset algorithm (e.g. [42, 29]) and then computes coresets for subsequent levels using a non-DP coreset Algorithm (e.g. [19]). Thus using our DP Merge and Reduce framework, at any given timestep, the number of input points we actually need to store is roughly equal to the size of the block at the base level or the threshold (denoted as $M$). This is because as soon as the number of input points exceed $M$, we can compute a DP coreset for this block. We perform this check privately via the sparse vector technique [23]. Note that the size of the base level block controls the height of the tree, i.e., larger $M$ results in smaller height of the tree which in turn results in smaller additive error of the resulting DP coreset. However, we show that we can store a sublinear number of input points and achieve reasonable additive error for a fixed instance of DP Merge and Reduce (see dp-merge-reduce for a formal statement).

Charging Additive Error to Multiplicative Error. Although we show that a fixed instance of DP Merge and Reduce for a fixed ring releases a DP semicoreset that achieves reasonably small additive error $\approx \tilde{O}(\frac{N_r}{M}+\log(T))$ where $N_r$ is the number of points in ring $R_r$, the additive error of the union of these DP semicoresets over all rings can be large, e.g., linear in $T$. We use ideas similar to [17] to address this challenge — by using the properties of the bicriteria solution and carefully choosing threshold $M$ for the DP Merge and Reduce instances, we show that we can charge the large additive error obtained by summing up the cost of clustering over all rings to a small constant-factor blow up in the multiplicative error (see final-accuracy-main). In particular, we achieve an additive error that is $\approx \tilde{O}(\poly(\log T))$.

Finally, observe that the bicriteria solution $\cF$ can change over time, i.e., new centers can be added to $\cF$. Thus whenever a new center (or centers) are added to $\cF$, we redefine the rings (based on distance to the set $\cF$) and initiate new DP Merge and Reduce instances for each ring. See extend-cluster for more details.

§ PRELIMINARIES

Norms and heavy hitters. Let $p \geq 1$, the $\ell_p$-norm of a vector $\mathbf{x}=(x_1, \ldots, x_t)$ is defined as $\|\mathbf{x}\|_p = (\sum^t_{i=1} \vert x_i \vert^p)^{1/p}$. Given a multiset $\cS$, denote the frequency of an item $x$ appearing in $\cS$ as $f(x)$. We say that an item $x$ is an $\theta$-heavy hitter ($\theta$-HH for short) if $f(x) \geq \theta \|\cS \|_1$.

Differential Privacy. Streams $\cS=(x_1,\ldots, x_T)$ and $\cS'=(x'_1,\ldots, x'_T)$ are neighboring if there exists at most one timestamp $t^* \in [T]$ for which $x_{t^*} \neq x'_{t^*}$ and $x_{t} = x'_{t}$ for all $t \neq t^*$.

*  (Sequential) Let $M_i$ each provide $(\eps_i,\delta_i)$-DP. The sequence of $M_i(X)$ provides $(\sum_i\eps_i, \sum_i \delta_i)$-DP.

* (Parallel) Let $M_i$ each provide $(\eps,\delta)$-DP. Let $D_i$ be arbitrary disjoint subsets of input domain $D$. The sequence of $M_i(X \cap D_i)$ provides $(\eps,\delta)$-DP.

Let $\eps \geq 0, \gamma \in (0,0.5)$, there is an $\eps$-DP algorithm for the sum of the stream in the  continual release model. With probability $1-\xi$, the additive error of the output for every timestamp $t \in [T]$ is always at most $O(\frac{1}{\eps} \log^{2.5}(T) \log(\frac{1}{\xi}))$ and uses $O(\log T)$ space.

Let $\eps>0$, $\gamma_h \in (0,0.5)$, $0< \theta < 1$, $\xi \in (0,0.5)$. There is an $\eps$-DP algorithm in the streaming continual release model such that with probability at least $1-\xi$, it always outputs a set $H \subseteq \cU$ and a function $\hat{f} : H \to \bbR$ for every timestamp $t \in [T]$ such that

* $\forall a \in H$, $\hat{f}(a) \in (1\pm \gamma_h) \cdot f_a$ where $f_a$ is the frequency of $a$ in the stream $\cS = (a_1, a_2, \ldots, a_t)$

* $\forall a \in \cU$, if $f_a \geq \frac{1}{\eps \gamma_h} \poly\bigl(\log\bigl( \frac{T \cdot\vert \cU \vert }{\theta \xi \gamma_h }\bigr)\bigr)$ and $f^1_a \geq \theta  \|\cS \|_1$ then $a \in H$

* The size of $H$ is at most $O((\log(T/\xi)+\log \vert \cU \vert ) \cdot (\frac{1+\gamma_h}{1-\gamma_h}) \cdot \frac{1}{\theta})$

The algorithm uses $\frac{1 }{\gamma_h^2\theta^3} \poly \left( \log\left( \frac{T \cdot  \vert \cU \vert}{\xi \theta }\right)\right)$ space.

Clustering.  For points $x,y \in \bbR^d$, we let $d(x,y) = \|x-y\|_2$ be the Euclidean distance between $x$ and $y$. Given a set $\cC$, we define $d(x,\cC):= \min_{c \in \cC} d(x,c)$.

For a set of centers $C$, we define the cost of clustering for the set $\cS$ wrt $C$ as

$$cost(C, \cS) = \sum_{x \in \cS} d^z(x,C)$$

where $z=1$ for $k$-median, and $z=2$ for $k$-means.

Our goal in DP clustering is to produce a set of $k$ centers $C_\cS$ for input stream $\cS$ such that (1) $C_\cS$ is $(\eps,\delta)$-DP wrt $\cS$, and (2) $cost(C_\cS, \cS) \leq \alpha \cdot cost(C^{opt}_\cS, \cS) + \beta$.

Given a point set $P$ in $\bbR^d$, a weighted subset $Q \subseteq P$ is a $(\kappa,\eta)$-coreset of $P$ for $k$-clustering ($k$-means or $k$-median) if

$$ \frac{1}{\kappa}\cdot \textsf{cost}(C,P)-\eta \leq \textsf{cost}(C,Q) \leq \kappa \cdot \textsf{cost}(C,P)+\eta$$

for any set of $k$-centers $C \subseteq \bbR^d$.

In particular, given a point set $P$ in $\bbR^d$, a weighted subset $Q \subseteq P$ is a $(1+\gamma,\eta)$-coreset of $P$ for $k$-clustering ($k$-means or $k$-median) if

$$ (1-\gamma) \cdot \textsf{cost}(C,P)-\eta \leq \textsf{cost}(C,Q) \leq (1+\gamma) \cdot \textsf{cost}(C,P)+\eta$$

for any set of $k$-centers $C \subseteq \bbR^d$.

Given a point set $P$ in $\bbR^d$, a weighted subset $Q \subseteq P$ is a $(\kappa,\eta_1,\eta_2)$-semicoreset of $P$ for $k$-clustering ($k$-means or $k$-median) if

$$ \frac{1}{\kappa}\cdot \textsf{cost}(C,P)-\eta_1 \cdot \textsf{cost}(C^{opt}_P,P) -\eta_2 \leq \textsf{cost}(C,Q) \leq \kappa \cdot \textsf{cost}(C,P)+\eta_2 $$

for any set of $k$-centers $C \subseteq \bbR^d$.

Let $\cA$ be an $(\eps,\delta)$-DP (static) algorithm that takes as input a dataset $P$ of size at least $\Omega(k \cdot \eps^{-1} \log(n))$ contained in a ball $\cB(0,\Lambda)$ (with known center) and outputs a set of $k$ centers $\cC$ contained in $\cB$ such that $cost(\cC,\cX ) \leq O(1) cost(\cC^{opt},\cX)+V(n,d,k,\eps,\delta) \cdot \Lambda^p$. Then for any fixed $\eta>0$ there exists a $(3\eps,\delta)$-DP (sequential) algorithm $\cB$ on the dataset $\cY$ that outputs a set $\hat{\cP} \subset B$ such that for any set $\cC$ of size at most $k$.

\begin{align}
cost(\hat{P}, \cC) &\leq O(1) \cdot \left( cost(\cC^{opt}_P,P) + cost(\cC,P)+(V(n,d,k,\eps,\delta)+k \eps^{-1} \log(n) )\cdot \Lambda^p \right) \\
cost(P,\cC) &\leq O(1) \cdot \left( cost(\cC^{opt}_P,P) + cost(\cC,\hat{P}) +(V(n,d,k,\eps,\delta)+k \eps^{-1} \log(n))\cdot \Lambda^p \right)
\end{align}

In addition constructing the set requires only $O(k^{1+\eta} \cdot d \cdot \log(n))$ space and at most $O(k^\eta \cdot d \cdot \log(n))$ time per point in $P$ for some arbitrarily chosen constant $\eta>0$.

Any $(\kappa,\eta)$-coreset is also a $(\kappa,0,\eta)$-semicoreset. The output of any DP $k$-clustering algorithm can be made into a DP semicoreset via clustering-to-semicoreset.

There exists a non-DP $(1+\gamma)$-coreset of size $O(k\log k \cdot (\gamma^{-2-\max(2,z)}) \cdot 2^{O(z \log z)} \cdot \poly\log(\gamma^{-1}))$ for $(k,z)$-clustering in Euclidean spaces.

§ DIFFERENTIALLY PRIVATE CLUSTERING FRAMEWORK

Our Differentially Private Clustering Framework in the continual release setting is given by extend-cluster. There is a postprocessing step of applying a $\rho$-approximation non-DP clustering algorithm on the output of extend-cluster. Notably, our framework allows us to plug in any existing offline DP coreset or clustering algorithm to obtain a corresponding DP clustering algorithm in the continual release setting.

Let $r \in \bbR$. Ring $R_r$ for set $\cF$ contains the set of points $\{x_i\}_{i \in [T]}$ such that $2^{r-1} \leq d(x_i,\cF)< 2^{r}$.

Main Algorithm

[1]
Stream $\cS$ of points $x_1, \ldots, x_T \in \bbR^d$, Privacy parameters $\eps,\eps_1,\delta_1$

Initialize bicriteria solution $\cF=\emptyset$, and DP coreset $\hat{\cY}=\emptyset $
when new point $x_t$ arrives
$\cF_t \leftarrow $Update$(x_t)$ of $\mathsf{\textbf{BicriteriaDPCenters}}(\eps)$See dp-cluster-online
$\cF_t \neq \emptyset$ and $\vert \cF \cap \cF_t \vert < \vert \cF_t \vert $ New centers added to $\cF$ — need to redefine rings
$\cF \leftarrow \cF \cup \cF_t$
Run non-DP coreset algorithm on the semicoreset $\hat{\cY}$ of points so far and store resulting semicoreset
Delete existing input points (if any) in memory

$1 \leq r \leq \log \Lambda$, run the following in parallel
Initialize $\hat{\cY}_{r}=\emptyset$
Let $R_r$ represent the ring centered at $\cF$ (see ring-center)
Create new instance $\textbf{DP-Merge-Reduce}_{R_r}$ of $\textbf{DP-Merge-Reduce}(\eps,\eps_1,\delta_1)$
$ x_t \in R_r$
$\hat{\cY}_{r} \leftarrow$ Update$(x_t)$ of $\textbf{DP-Merge-Reduce}_{R_r}$ See dp-merge-reduce
$\hat{\cY}_{r} \leftarrow$ Update$(\perp)$ of $\textbf{DP-Merge-Reduce}_{R_r}$
Set $\hat{\cY} \leftarrow \hat{\cY} \cup (\cup_{r} \hat{\cY}_{r})$
Output the semicoreset computed by running a non-DP coreset algorithm on $\hat{\cY}$

$\vert \cF \cap \cF_t \vert = \vert \cF_t \vert$ No new centers — the rings are already defined wrt existing centers
$\cF \leftarrow \cF \cup \cF_t$
$1 \leq r \leq \log \Lambda$, run the following in parallel

Let $R_r$ represent the ring centered at $\cF$ (see ring-center)

$ x_t \in R_r$  $\textbf{DP-Merge-Reduce}_{R_r}$ has already been created
$\hat{\cY}_{r} \leftarrow$ Update$(x_t)$ of $\textbf{DP-Merge-Reduce}_{R_r}$
$\hat{\cY}_{r} \leftarrow$ Update$(\perp)$ of $\textbf{DP-Merge-Reduce}_{R_r}$
Set $\hat{\cY} \leftarrow \hat{\cY} \cup (\cup_{r} \hat{\cY}_{r})$
Output the semicoreset computed by running a non-DP coreset algorithm on $\hat{\cY}$

Main Algorithm. When a new point $x_t$ arrives, our algorithm does the following

*  Update the bicriteria solution $\cF$ (see dp-cluster-online)

* Create rings according to ring-center for the set $\cF$ and add $x_t$ to an existing ring[Notice that in the pseudo code extend-cluster the symbol $\perp$ represents an empty update that is effectively ignored. This is needed for technical reasons to ensure DP by avoiding the value of the input affecting the number of events in the sub-streams.]. For each ring $1 \leq r \leq \log(\Lambda)$ maintain an instance of $\textsf{DP-Merge-Reduce}_r$ (see dp-merge-reduce) which outputs a DP semicoreset per ring. If no new centers have been added to $\cF$ in this timestep, then instead of creating new rings, the algorithm adds $x_t$ to an existing ring (and corresponding  instance).

* Release the union of these DP semicoresets as $\hat{\cY}$.

If new centers have been added to $\cF$, then in order to keep our space usage small, before creating the new rings, we apply a $(1+\gamma)$-approximation non-DP coreset algorithm to the existing union of semicoresets of these rings as $\hat{\cY}$. We do the same thing, i.e., apply $(1+\gamma)$-approximation non-DP coreset algorithm before releasing $\hat{\cY}$ at every timestep in order to keep the size of $\hat{\cY}$ small.

Finally, in an offline postprocessing step, we apply a $\rho$-approximate non-DP clustering algorithm to $\hat{\cY}$.

For the sake of analysis, we split the entire stream of points into epochs dictated by the addition of new centers to the bicriteria solution $\cF$.
Let $T_1, \ldots, T_e$ be the epochs such that for a fixed $i$, the set of bicriteria centers $\cF$ over the timesteps $t \in T_i$ is fixed. Clearly $T_1 \cup \ldots \cup T_e = [T]$ and $T_1 \cap \ldots \cap T_e = \emptyset$.

We first state the theoretical guarantees of  and  as we need these results to state the guarantees of our Main Algorithm (see extend-cluster). The proofs for these statements can be found in final-acc-dp-hh and dp-merge-reduce.

Let $\cS:= \{x_1,\ldots, x_T\}$ be the stream of input points in Euclidean space. For $t\in [T]$, let $\cF_t$ be the set of centers until time step $t$. Let $cost(\cF, \cS):=\sum^T_{t=1} cost(\cF_t)$ where $cost(\cF_t) := \min_{f \in \cF_t} dist^2(x_t,f)$. There exists an algorithm  (see dp-cluster-online) that outputs a set of centers $\cF$ at every timestep $t\in [T]$ such that

* (Privacy)  is $\eps$-DP under the continual release setting.

* (Accuracy) With probability at least $1-1/k^{O(\poly(k,\log(\Lambda))}$,

\begin{align*}
cost(\cF, \cS) &\leq  O(d^3)cost(C^{opt}_{\cS}, \cS) + \tilde{O}\left(\frac{d^{2} \Lambda^2 k}{\eps} \poly\left(\log\left({T\cdot k \cdot\Lambda}\right)\right)\right)
\end{align*}

where $cost(C^{opt}_{\cS}, \cS)$ is the optimal $k$-means cost for $\cS$.

* (Space)  uses $O(k \log(\Lambda) \log^2(k) \poly\left( \log \left({T}{\Lambda}k\right)\right))$ space.

* (Size) $\cF$ has at most $O({k}\log^2(k) \log(\Lambda) \log T)$ centers.

Recall that our Main Algorithm runs multiple concurrent instances of  at any given timestep, so $N$ in the theorem statement below represents the number of points seen by a specific instance of . We use a non-DP coreset algorithm in  which we fix in the analysis to be the current state-of-the-art clustering result by nondp-coreset.

Let $0<\xi<1$, $T$ be the length of the entire stream, $\eps,\delta$ be privacy parameters, $M$ be an arbitrary parameter such that  $M>\frac{12}{\eps}\log(\frac{2T}{\xi})$, and $P$ be a sub-stream of non-empty points with length $N$.

Suppose we are given black-box access to an offline $(\eps_1,\delta_1)$-DP algorithm $\cA$ that computes a $(\kappa,\eta_1,\eta_2)$-semicoreset of $X \subseteq \bbR^d$ of size $SZ_\cA(N,k,d,\eps, \delta, \kappa,\eta_1,\eta_2, \xi_A)$ using space $S_\cA(N,k,d,\eps, \delta, \kappa,\eta_1,\eta_2, \xi_A)$ with failure probability $\xi_A$. And we are given black-box access to an offline non-DP algorithm $\cB$ (see nondp-coreset) that computes a $(1+\gamma)$-coreset in time $X \subseteq \bbR^d$ of size $SZ_{\cB}(N,k,d,\gamma,\xi_B)$ using space $S_\cB(N,k,d,\gamma, \xi_B)$ with failure probability $\xi_B$.
Then there exists an algorithm  (see dp-merge-reduce) in the streaming model such that

* (Privacy)  is $(\eps+\eps_1,\delta_1)$-DP under the continual release setting.

* (Accuracy) With probability $1-\xi_A-\xi_B-\xi$, the semicoreset released by  is  a $((1+\gamma)\kappa,(\frac{4N}{M}-1)(1+\gamma)\eta+\tilde{M})$-semicoreset of $P$. Where $\tilde{M}:= {M}+\frac{6}{\eps}\log(\frac{2T}{\xi})$.

* (Space)  requires $S_\cA(M, k,d,\eps,\delta,\kappa,\eta_1,\eta_2, \xi_A)+ S_\cB(SZ_\cA(M,k,d,\eps,\delta,\kappa,\eta_1,\eta_2, \xi_A),k,d, \gamma, \xi_B)+\lceil \log(2N/M) \rceil \cdot S_\cB(SZ_\cB(M,k,d,\gamma, \xi_B),k,d, \gamma, \xi_B) +3M/2$ space.

* (Size) The resulting coreset has size at most  $O(k\log k \cdot \gamma^{-4})$.

The parameter $M$ denotes the block size of the base level in  (see dp-merge-reduce). We treat $M$ as an arbitrary parameter in dp-merge-reduce, main-space and yhat-dp-merge-reduce-output-semicoreset. However we set $M$ to be a function of $\alpha$ (the multiplicative approximation error of the bicriteria solution), $\eta_2$ (the additive approximation error of the DP semicoreset algorithm) and a parameter $C_M$ in the proof of yhat-semicoreset (see M). We assign an appropriate value to $C_M$ and consequently $M$ in order to obtain an $O(1)$-multiplicative approximation in final-stemmer and a $(1+\gamma')$-approximation in final-ghazi.

Privacy. We first show that the output of our main algorithm is indeed differentially private.

Let the underlying DP semicoreset algorithm used as a blackbox be $(\eps_1,\delta_1)$-DP, then extend-cluster is $(\eps_1+2\eps,\delta_1)$-DP under the continual release setting.

We first observe that for a fixed input point $x_t$, we can view extend-cluster as two main steps (1) $x_t$ is first used to update the subroutine  which produces a set of candidate centers $\cF$ (2) $x_t$ is assigned to an appropriate ring $R_r$ (defined according to $\cF$) and added to the corresponding $_{R_r}$ instance of that ring which produces a DP semicoreset (3) Finally, the algorithm releases the union of these DP semicoresets over disjoint rings at timestep $t$.

Step (1) involves $x_t$ being processed by  which is $\eps$-DP under the continual release setting (by final-acc-dp-hh). In Step (2), since the rings partition the input space by definition, the corresponding $_{R_r}$ instances are disjoint. Since the semicoreset output by each $_{R_r}$ instance is $(\eps_1+\eps,\delta_1)$-DP under the continual release setting (by dp-merge-reduce), then by parallel composition (Item <ref> of composition), the total privacy over all semicoresets output by $_{R_r}$ instances (where $1 \leq r \leq \Lambda$) is $(\eps_1+\eps,\delta_1)$-DP. In Step (3), the union of DP semicoresets is DP by postprocessing.

Thus, each step is private under the continual release setting and note that each step is carried out in the same fixed timestep $t$.  By sequential composition (Item <ref> of composition) of the privacy for Steps (1) and (2), and postprocessing of Step (3), our claim follows.

Space and Size. We present the space (see main-space) and coreset size (see coreset-size-main) of extend-cluster in terms of the space and coreset size of the underlying DP semicoreset algorithm $\cA$ and non-DP coreset algorithm $\cB$ used in . We use the algorithm from [19] as our non-DP coreset algorithm $\cB$ whose guarantees are given by nondp-coreset.

extend-cluster consumes

\begin{align*}%\label{eq:space}
&\log(\Lambda) \cdot (S_\cA(M, k,d,\eps,\delta,\kappa,\eta_1,\eta_2,\xi_\cA)+ S_\cB(SZ_\cA(M,k,d,\eps,\delta,\kappa,\eta_1,\eta_2,\xi_\cA),k,d, \gamma, \xi_B)+3M/2) \nonumber\\&+\tilde{O}(k \log(T/M)\log(\Lambda)\log k \cdot \gamma^{-4}) +O({k}\log(k) \log^2(\Lambda) \log T)+O(k \log(\Lambda) \log^2(k) \poly\left( \log \left({T}{\Lambda}k\right)\right))
\end{align*}

The last term in the above claim is the total space used by  and the second last term is the space used to store the bicriteria solution $\cF$ (see dp-hh-cr). We focus on proving that the  instances consume the space specified by the first term of the above claim.

We sometimes abuse notation and omit the input for the non-DP coreset size $SZ_\cB$ as by nondp-coreset we know that $SZ_\cB(\cdot) = \tilde{O}(k \log(k) \cdot \gamma^{-4})$. For a fixed epoch $T_i$, recall that the set of centers $\cF_t$ is fixed for $t \in T_i$. According to extend-cluster, at timestep $t \in T_i$, we compute $ \log(\Lambda)$ instances of  in parallel. By dp-merge-reduce, since the space required to compute semicoreset $\hat{\cY}^{(t)}_{r}$ using  is
$S_\cA(M, k,d,\eps,\delta,\kappa,\eta_1,\eta_2,\xi_\cA)+ S_\cB(SZ_\cA(M,k,d,\eps,\delta,\kappa,\eta_1,\eta_2,\xi_\cA),k,d, \gamma,\xi_\cB)+\lceil \log(2N/M) \rceil \cdot S_\cB(SZ_\cB(\cdot),k,d, \gamma, \xi_\cB) +3M/2$, the total space at timestep $t \in T_i$ is

\begin{align}
& \sum^{\log(\Lambda)}_{r=1}  (S_\cA(M, k,d,\eps,\delta,\kappa,\eta_1,\eta_2,\xi_\cA)+ S_\cB(SZ_\cA(M,k,d,\eps,\delta,\kappa,\eta_1,\eta_2,\xi_\cA),k,d, \gamma,\xi_\cB)\nonumber\\&+\lceil \log(2 N^{(t)}_r/M) \rceil \cdot S_\cB(SZ_\cB(\cdot),k,d, \gamma,\xi_\cB) +3M/2)  \nonumber \\
&=\log(\Lambda) \cdot (S_\cA(M, k,d,\eps,\delta,\kappa,\eta_1,\eta_2,\xi_\cA)+ S_\cB(SZ_\cA(M,k,d,\eps,\delta,\kappa,\eta_1,\eta_2,\xi_\cA),k,d, \gamma,\xi_\cB)+3M/2) \nonumber \\&+ S_\cB(SZ_\cB(\cdot),k,d, \gamma,\xi_\cB) \sum^{\log(\Lambda)}_{r=1} \lceil \log(2 N^{(t)}_r/M) \rceil \label{eq:sum-r-other}
\end{align}

We focus on the last term in sum-r-other. In particular, since the space used by $\cB$ is linear in the coreset size, we have that $S_\cB(SZ_\cB(\cdot),k,d, \gamma)=\tilde{O}(k \log k \cdot \gamma^{-4})$. Next we simplify the adjoining sum

\begin{align}
&\sum^{\log(\Lambda)}_{r=1} \lceil \log(\frac{2N^{(t)}_r}{M}) \rceil\leq \sum^{\log(\Lambda)}_{r=1} (\log(\frac{2N^{(t)}_r}{M}) +1) = \log\left(\frac{2^{\log(\Lambda)}}{M^{\log(\Lambda)}}\cdot \prod^{\log(\Lambda)}_{r=1}N^{(t)}_r\right) + \log(\Lambda) \nonumber\\
&< \log\left(\frac{\Lambda \cdot 2^{\log(\Lambda)}}{M^{\log(\Lambda)}}\cdot T^{\log(\Lambda)}\right) = \log(\Lambda)(2+\log(T)-\log(M))
\end{align}

where we use the fact that the number of points in any ring $N_r$ cannot be larger than $T$ in the second last step.

extend-cluster releases a coreset of size at most $\tilde{O}(k \log(k)\cdot \gamma^{-4})$.

Consider a fixed epoch $T_i$. By dp-merge-reduce, the coreset $\hat{\cY}^{(t)}_{r}$ has size $\tilde{O}(k \log(k)\cdot \gamma^{-4})$ for $t\in T_i$. Since we run a non-DP coreset algorithm at the end of each epoch (whose theoretical guarantees are given by nondp-coreset), the size of the coreset at the end of the epoch is also $\tilde{O}(k \log(k)\cdot \gamma^{-4})$. Recall that an epoch is created every time a new center is added to $\cF$, therefore the total size of the coreset is given by $\vert \cF \vert \cdot\tilde{O}(k \log(k)\cdot \gamma^{-4})$. But since we apply another non-DP coreset algorithm to $\hat{\cY}$ before releasing it, the coreset size is $\tilde{O}(k \log(k)\cdot \gamma^{-4})$.

Accuracy. We analyze the accuracy of extend-cluster in the sequel. We first state the accuracy guarantee of the DP semicoreset $\hat{\cY}^{(t)}_{r}$ released by $\textsf{DP-Merge-Reduce}_r$ for each ring $R^{(t)}_r$ at timestep $t \in [T]$ in yhat-dp-merge-reduce-output-semicoreset.

Let $\hat{\cY}^{(t)}_{r}$ be the output of $\textsf{DP-Merge-and-Reduce}_r$ (see extend-cluster) and $N^{(t)}_{r}$ be the number of (non-empty) points in $R^{(t)}_r$ at timestep $t \in [T]$. Then $\hat{\cY}^{(t)}_{r}$ is a DP $((1+\gamma)\kappa,(1+\gamma)\eta_1, (\frac{4N^{(t)}_r}{M}-1)(1+\gamma)\eta_2+\tilde{M})$-semicoreset of $R^{(t)}_r$ at timestep $t$. Where  $\tilde{M}:= {M}+\frac{6}{\eps}\log(\frac{2T}{\xi})$. In other words, for any set of $k$ centers $\cC$, with probability $1-3\xi$,

\begin{align}\label{eq:yhat-dp-merge-reduce-output-semicoreset}
&\frac{(1-\gamma)}{\kappa}\cdot \mathsf{cost}(\cC,R^{(t)}_{r})-  (1-\gamma)\eta_1 \cdot \textsf{cost}(\cC^{opt}_{R^{(t)}_r},R^{(t)}_r)\cdot  (2^r)^2 - ((\frac{4N^{(t)}_{r}}{M}-1)(1+\gamma)\eta_2 {+\tilde{M}})\cdot (2^r)^2 \nonumber \\&\leq \mathsf{cost}(\cC,\hat{\cY}^{(t)}_{r})\nonumber\\ &\leq (1+\gamma)\kappa \cdot \mathsf{cost}(\cC,R^{(t)}_{r}) + ((\frac{4N^{(t)}_{r}}{M}-1)(1+\gamma)\eta_2 {+\tilde{M}}) \cdot (2^r)^2
\end{align}

The lemma follows from the accuracy guarantees of  (see dp-merge-reduce) and the fact that the additive cost incurred is proportional to the radius of the ring $R_r$ which is given by $2^r$.

Next, we show that the union of the DP semicoresets over all rings $\hat{\cY}$ is a semicoreset for the stream $\cS$.

Given dimension $d$, clustering parameter $k$, arbitrary parameter $C_M$, non-DP $(1+\gamma)$-coreset, $O(d^3)$-approximate bicriteria solution from final-acc-dp-hh, DP $(\kappa,\eta_1,\eta_2)$-semicoreset, and privacy parameter $\eps$. Let $\hat{\cY}$ be the output of extend-cluster for stream $\cS=\{x_1, \ldots, x_T\}$. Then for any set of $k$ centers $\cC$, we have that with probability $1-\frac{1}{T^2}-\frac{1}{k^{O(\poly(k,\log(\Lambda))}}$, the following holds

\begin{align}\label{eq:yhat-semicoreset}
&\frac{(1-\gamma)^2}{\kappa} \mathsf{cost}(\cC,\cS)- ((1-\gamma)^2\eta_1\Lambda^2 + 16 C_M(1-\gamma^2))  \cdot \mathsf{cost}(\cC^{opt}_\cS, \cS) -C_M(1-\gamma^2) \cdot V'(d,k,\eps,T,\Lambda) ) \nonumber \\&- (1-\gamma)V''(d,k,\eps,T,\Lambda) \nonumber \\&\leq \mathsf{cost}(\cC,\hat{\cY})\nonumber\\ &\leq (1+\gamma)^2\kappa \mathsf{cost}(\cC,\cS) +16 C_M(1+\gamma)^2\mathsf{cost}(\cC^{opt}_\cS, \cS) + C_M(1+\gamma)^2 \cdot V'(d,k,\eps,T,\Lambda) )+(1+\gamma)V''(d,k,\eps,T,\Lambda)
\end{align}

where $V'(d,k,\eps,T,\Lambda) = O\left(\frac{\Lambda^2 k^2}{d\eps} \log^2(\Lambda) \log^4(k) \log(T) \poly\log\left({T\cdot k \cdot\Lambda}\right)\right)$, $V''(d,k,\eps, T,\Lambda)= O({k}\Lambda^2\log^2(k) \log(\Lambda) \log T)\cdot (\frac{\alpha\eta_2}{C_M}+\tilde{O}(\frac{1}{\eps} \cdot \poly(\log(T\cdot k\cdot \Lambda))$, and $\alpha:=O(d^3)$.

Fix an epoch $T_i$. Let $\hat{\cY}^{(t)}$ be the output of extend-cluster at timestep $t \in T_i$. Since the centers in $\cF\vert_{T_i}$ and consequently the rings $R_r$ are fixed, we can compute the cost of $\hat{\cY}^{(t)}$ by summing yhat-dp-merge-reduce-output-semicoreset from yhat-dp-merge-reduce-output-semicoreset over all rings $R_r$ where $1 \leq r \leq \log(\Lambda)$ as follows:

\begin{align}\label{eq:bjr-yt-fixed-semicoreset}
&\frac{(1-\gamma)}{\kappa}\sum_r \mathsf{cost}(\cC,R^{(t)}_{r})-  (1-\gamma)\eta_1 \sum_r \textsf{cost}(\cC^{opt}_{R^{(t)}_r},R^{(t)}_r)\cdot  (2^r)^2 - (1+\gamma)\eta_2\sum_r(\frac{4N^{(t)}_{r}}{M}-1)\cdot (2^r)^2 -\sum_r\tilde{M}\cdot (2^r)^2 \nonumber \\&\leq \mathsf{cost}(\cC,\hat{\cY}^{(t)})\nonumber\\ &\leq (1+\gamma)\kappa \sum_r \mathsf{cost}(\cC,R^{(t)}_{r}) + (1+\gamma)\eta_2\sum_r(\frac{4N^{(t)}_{r}}{M}-1) +\sum_r \tilde{M} \cdot (2^r)^2
\end{align}

Next, we take the union of $\hat{\cY}^{(t)}$ over all timesteps in epoch $T_i$. Note that for a fixed epoch and fixed ring $R_r$, the additive error $\tilde{M}$ is incurred at most once since this error stems from privately testing if the number of points at the base level is larger than the block size $M$ (see dp-above-threshold) in $\textsf{DP-Merge-Reduce}_r$. Thus summing bjr-yt-fixed-semicoreset over all timesteps in $T_i$ gives us

\begin{align}\label{eq:bjr-yt-all-semicoreset}
&\frac{(1-\gamma)}{\kappa} \mathsf{cost}(\cC,\cS\vert_{T_i})-  (1-\gamma)\eta_1\Lambda^2  \textsf{cost}(\cC^{opt}_{\cS\vert_{T_i}},\cS\vert_{T_i})- (1+\gamma)\eta_2\sum_{t \in T_i}\sum_r(\frac{4N^{(t)}_{r}}{M}-1)\cdot (2^r)^2 -\tilde{M}\cdot\frac{4}{3}(\Lambda^2-1) \nonumber \\&\leq \mathsf{cost}(\cC,\hat{\cY}^{(t)}\vert_{t\in T_i})\nonumber\\ &\leq (1+\gamma)\kappa \mathsf{cost}(\cC,\cS\vert_{T_i}) + (1+\gamma)\eta_2\sum_{t \in T_i}\sum_r(\frac{4N^{(t)}_{r}}{M}-1) +\tilde{M} \cdot \frac{4}{3}(\Lambda^2-1)
\end{align}

where $\cS\vert_{T_i}$ denotes the input points in stream $\cS$ restricted to epoch $T_i$ and $\hat{\cY}^{(t)}\vert_{t\in T_i}$ is defined analogously.

If $x_i$ is assigned to the ring $R_{r}$, let $r(x_i):= 2^r$. Then

\begin{align}
\sum_{x_i \in \cS} r(x_i)^2 \leq \sum_{x_i \in \cS} 4d(x_i,\cF)^2 = 4\textsf{cost}(\cF,\cS) \leq 4\alpha \cdot \textsf{cost}(\cC^{opt}_\cS,\cS)+ V(d,k,\eps,T,\Lambda)
\end{align}

where $\alpha:=O(d^3)$ and $V(d,k,\eps,T,\Lambda)=  O\left(\frac{d^{2} \Lambda^2 k \log(\Lambda) \log^2(k)}{\eps} \poly\bigl(\log\bigl({T\cdot k \cdot  \Lambda}\bigr)\bigr)\right)$ are the multiplicative/additive factors from the bicriteria approximation (see final-acc-dp-hh).

The statement follows immediately from the definition of the ring $R_r$ (ring-center) and the accuracy guarantees of the bicriteria solution (final-acc-dp-hh).

Observe that $\sum_{t\in T_i}\sum_{r} N^{(t)}_{r}\cdot (2^r)^2 = \sum_{x \in \cS\vert_{T_i}} r(x)^2$, thus we can use radius-opt-cost to simplify the additive error in bjr-yt-all-semicoreset to obtain

\begin{align}\label{eq:bjr-yt-epoch-semicoreset}
&\frac{(1-\gamma)}{\kappa} \mathsf{cost}(\cC,\cS\vert_{T_i})-  (1-\gamma)\eta_1\Lambda^2  \textsf{cost}(\cC^{opt}_{\cS\vert_{T_i}},\cS\vert_{T_i})- \frac{4}{M}(1+\gamma)\eta_2 (4\alpha \cdot \mathsf{cost}(\cC^{opt}_{\cS \vert_{T_i}}, \cS\vert_{T_i}) +V(d,k,\eps,T,\Lambda) )-\tilde{M}\cdot\frac{4 \Lambda^2}{3} \nonumber \\&\leq \mathsf{cost}(\cC,\hat{\cY}^{(t)}\vert_{t\in T_i})\nonumber\\ &\leq (1+\gamma)\kappa \mathsf{cost}(\cC,\cS\vert_{T_i})  + \frac{4}{M}(1+\gamma)\eta_2 (4\alpha \cdot \mathsf{cost}(\cC^{opt}_{\cS \vert_{T_i}}, \cS\vert_{T_i}) +V(d,k,\eps,T,\Lambda) ) +\tilde{M} \cdot \frac{4\Lambda^2}{3}
\end{align}

Since a new epoch starts whenever a new center is added to the bicriteria solution $\cF$, and we run a $(1+\gamma)$-approx non-DP coreset algorithm on the coreset $\hat{\cY}$ at the beginning of a new epoch, denote this new coreset as $\hat{\cY}_{T_i}$. We state the guarantees of $\hat{\cY}_{T_i}$ below.

\begin{align}\label{eq:bjr-new-epoch-semicoreset}
&\frac{(1-\gamma)^2}{\kappa} \mathsf{cost}(\cC,\cS\vert_{T_i})-  (1-\gamma)^2\eta_1\Lambda^2  \textsf{cost}(\cC^{opt}_{\cS\vert_{T_i}},\cS\vert_{T_i})- \frac{4}{M}(1-\gamma^2)\eta_2 (4\alpha \cdot \mathsf{cost}(\cC^{opt}_{\cS \vert_{T_i}}, \cS\vert_{T_i}) +V(d,k,\eps,T,\Lambda) )\nonumber\\&-\frac{4 \Lambda^2 (1-\gamma)\tilde{M}}{3} \nonumber \\&\leq \mathsf{cost}(\cC,\hat{\cY}_{T_i})\nonumber\\ &\leq (1+\gamma)^2\kappa \mathsf{cost}(\cC,\cS\vert_{T_i}) + \frac{4}{M}(1+\gamma)^2\eta_2 (4\alpha \cdot \mathsf{cost}(\cC^{opt}_{\cS \vert_{T_i}}, \cS\vert_{T_i}) +V(d,k,\eps,T,\Lambda) )  +\frac{4\Lambda^2 (1+\gamma)\tilde{M}}{3}
\end{align}

Observe that the total number of epochs is bounded by the total number of centers in $\cF$ at the end of the stream. By final-acc-dp-hh we know that $\vert \cF \vert = O({k}\log^2(k) \log(\Lambda) \log T)$. Therefore, we can sum bjr-new-epoch-semicoreset over all epochs to have that for any set of $k$ centers $\cC$:

\begin{align}\label{eq:bjr-yt-epoch-all-semicoreset}
&\frac{(1-\gamma)^2}{\kappa} \mathsf{cost}(\cC,\cS)-  (1-\gamma)^2\eta_1\Lambda^2  \textsf{cost}(\cC^{opt}_{\cS},\cS)- \frac{4}{M}(1-\gamma^2)\eta_2 (4\alpha \cdot \mathsf{cost}(\cC^{opt}_\cS, \cS) +V^*(d,k,\eps,T,\Lambda) )\nonumber\\&-{\vert \cF \vert} \cdot \frac{4 \Lambda^2 (1-\gamma)\tilde{M}}{3} \nonumber \\&\leq \mathsf{cost}(\cC,\hat{\cY})\nonumber\\ &\leq (1+\gamma)^2\kappa \mathsf{cost}(\cC,\cS) + \frac{4}{M}(1+\gamma)^2\eta_2 (4\alpha \cdot \mathsf{cost}(\cC^{opt}_\cS, \cS) +V^*(d,k,\eps,T,\Lambda) )+{\vert \cF \vert} \cdot \frac{4\Lambda^2 (1+\gamma)\tilde{M}}{3}
\end{align}

where $V^*(d,k,\eps,T,\Lambda) = O\left(\frac{d^{2} \Lambda^2 k^2}{\eps} \log^2(\Lambda) \log^4(k) \log(T) \poly\left(\log\left({T\cdot k \cdot\Lambda}\right)\right)\right)$.

Grouping like terms together, we have

\begin{align}\label{eq:bjr-yt-epoch-all-semicoreset-final}
&\frac{(1-\gamma)^2}{\kappa} \mathsf{cost}(\cC,\cS)-  ((1-\gamma)^2\eta_1\Lambda^2 + \frac{16\alpha}{M}(1-\gamma^2)\eta_2)  \cdot \mathsf{cost}(\cC^{opt}_\cS, \cS) -\frac{4}{M}(1-\gamma^2)\eta_2 \cdot V^*(d,k,\eps,T,\Lambda) )\nonumber\\&-{\vert \cF \vert} \cdot \frac{4 \Lambda^2 (1-\gamma)\tilde{M}}{3} \nonumber \\&\leq \mathsf{cost}(\cC,\hat{\cY})\nonumber\\ &\leq (1+\gamma)^2\kappa \mathsf{cost}(\cC,\cS) + \frac{4}{M}(1+\gamma)^2\eta_2 (4\alpha \cdot \mathsf{cost}(\cC^{opt}_\cS, \cS) +V^*(d,k,\eps,T,\Lambda) )+{\vert \cF \vert} \cdot \frac{4\Lambda^2 (1+\gamma)\tilde{M}}{3}
\end{align}

We set

\begin{align}\label{eq:M}
M:= \frac{\alpha \eta_2 }{C_M}
\end{align}

where $C_M$ is a parameter chosen in the sequel (see param-M). Simplifying bjr-yt-epoch-all-semicoreset-final and taking a union bound over all rings and epochs, the desired claim in the theorem statement holds with probability $1- \log(\Lambda)\vert \cF \vert 3 \xi-\frac{1}{k^{O(\poly(k,\log(\Lambda))}}$. Thus we set $\xi:=  \frac{1}{3\vert \cF \vert \log(\Lambda)T^2}$.

Finally, final-accuracy-main gives the cost of clustering result for the output $\hat{\cY}$ after the offline postprocessing step is executed. Recall that the postprocessing step involves running a $\rho$-approximation non-DP clustering algorithm on $\hat{\cY}$. We present the clustering guarantee in terms of the $\rho$-approx non-DP clustering algorithm, the  non-DP $(1+\gamma)$-coreset algorithm, the DP $(\kappa, \eta_1,\eta_2)$-coreset algorithm and parameter $C_M$.

Given dimension $d$, clustering parameter $k$, arbitrary parameter $C_M$, non-DP $(1+\gamma)$-coreset, $O(d^3)$-approximate bicriteria solution from final-acc-dp-hh, DP $(\kappa,\eta_1,\eta_2)$-semicoreset, and privacy parameter $\eps$. Let $\cC_{\hat{\cY}}$ be the set of $k$ centers obtained from running the offline $\rho$-approx non-DP $k$-means algorithm on $\hat{\cY}$. Then,

\begin{align}
&\mathsf{cost}(\cC_{\hat{\cY}},\cS) \nonumber\\
&\leq \frac{\kappa}{(1-\gamma)^3}\cdot ((1+\gamma)^3\rho (\kappa + 16 C_M) +((1-\gamma)^3\eta_1\Lambda^2 + 16C_M(1-\gamma^2)(1-\gamma)))\mathsf{cost}(\cC^{opt}_{\cS},\cS) \nonumber\\&+ \frac{\kappa}{(1-\gamma)^3}\cdot(\rho C_M(1+\gamma)^3+ C_M(1-\gamma^2)(1-\gamma)) V'(d,k,\eps,T,\Lambda)\nonumber \\&+ \frac{\kappa}{(1-\gamma)^3}\cdot((1+\gamma)^2\rho+ (1-\gamma)^2) V''(d,k,\eps,T,\Lambda)
\end{align}

where $V'(d,k,\eps,T,\Lambda) = O\left(\frac{\Lambda^2 k^2}{d\eps} \log^2(\Lambda) \log^4(k) \log(T) \poly\log\left({T\cdot k \cdot\Lambda}\right)\right)$, $V''(d,k,\eps, T,\Lambda)= O({k}\Lambda^2\log^2(k) \log(\Lambda) \log T)\cdot \left(\frac{\alpha\eta_2}{C_M}+\tilde{O}(\frac{1}{\eps} \cdot \poly(\log(T\cdot k\cdot \Lambda)\right)$, and $\alpha:=O(d^3)$.

We first recall that our algorithm releases the semicoreset $\hat{\cY}$ after running a non-DP $(1+\gamma)$-coreset algorithm on it. Thus, the guarantee of the resulting semicoreset which we call $\hat{\cY}_{off}$ is

\begin{align}\label{eq:yhat-semicoreset-opt}
&\frac{(1-\gamma)^3}{\kappa} \mathsf{cost}(\cC,\cS)-  ((1-\gamma)^3\eta_1\Lambda^2 + 16C_M(1-\gamma^2)(1-\gamma))  \cdot \mathsf{cost}(\cC^{opt}_\cS, \cS) - C_M(1-\gamma^2)(1-\gamma) \cdot V'(d,k,\eps,T,\Lambda) ) \nonumber \\&- (1-\gamma)^2V''(d,k,\eps,T,\Lambda) \nonumber \\&\leq \mathsf{cost}(\cC,\hat{\cY}_{off})\nonumber\\ &\leq (1+\gamma)^3\kappa \mathsf{cost}(\cC,\cS) + 16 C_M(1+\gamma)^3\mathsf{cost}(\cC^{opt}_\cS, \cS) + C_M(1+\gamma)^3 \cdot V'(d,k,\eps,T,\Lambda) )+(1+\gamma)^2V''(d,k,\eps,T,\Lambda)
\end{align}

where $V'(d,k,\eps,T,\Lambda) = O\left(\frac{\Lambda^2 k^2}{d\eps} \log^2(\Lambda) \log^4(k) \log(T) \poly\log\left({T\cdot k \cdot\Lambda}\right)\right)$, $V''(d,k,\eps, T,\Lambda)= O({k}\Lambda^2\log^2(k) \log(\Lambda) \log T)\cdot (\frac{\alpha\eta_2}{C_M}+\tilde{O}(\frac{1}{\eps} \cdot \poly(\log(T\cdot k\cdot \Lambda))$, and $\alpha:=O(d^3)$.

Recall from the theorem statement that $\cC_{\hat{\cY}}$ is the set of $k$ centers obtained from running the offline $\rho$-approx non-DP $k$-means algorithm on $\hat{\cY}_{off}$ which gives us the following guarantee,

\begin{align}\label{eq:y-opt-sc}
\mathsf{cost}(\cC_{\hat{\cY}},\hat{\cY}_{off}) \leq \rho \cdot \mathsf{cost}(\cC^{opt}_{\hat{\cY}}, {\hat{\cY}_{off}})
\end{align}

Since yhat-semicoreset-opt is true for any set of $k$ centers, in particular it holds for $\cC_{\hat{\cY}}$, which gives us the following

\begin{align}\label{eq:opt-y-semicoreset}
&\frac{(1-\gamma)^3}{\kappa} \mathsf{cost}(\cC_{\hat{\cY}},\cS)-  ((1-\gamma)^3\eta_1\Lambda^2 + 16C_M(1-\gamma^2)(1-\gamma))  \cdot \mathsf{cost}(\cC^{opt}_\cS, \cS) -C_M(1-\gamma^2)(1-\gamma) \cdot V'(d,k,\eps,T,\Lambda) ) \nonumber \\&- (1-\gamma)^2V''(d,k,\eps,T,\Lambda) \nonumber \\&\leq \mathsf{cost}(\cC_{\hat{\cY}},\hat{\cY}_{off})\nonumber\\ &\leq (1+\gamma)^3\kappa \mathsf{cost}(\cC_{\hat{\cY}},\cS) + 16 C_M(1+\gamma)^3\mathsf{cost}(\cC^{opt}_\cS, \cS) + C_M(1+\gamma)^3 \cdot V'(d,k,\eps,T,\Lambda) )+(1+\gamma)^2V''(d,k,\eps,T,\Lambda)
\end{align}

We derive the desired expression in the theorem statement by starting with the LHS of opt-y-semicoreset and then applying y-opt-sc to obtain rho-cy-sc. By observing that the cost of $\hat{\cY}_{off}$ wrt its optimal set of centers $\cC^{opt}_{\hat{\cY}}$ is less than its cost wrt to $\cC^{opt}_{\cS}$, we obtain opt-sc. coreset-app is obtained by simply applying the RHS of the coreset relation from opt-y-semicoreset wrt $\cC^{opt}_{\cS}$. The final statement is derived by grouping like terms and simplifying simplify-sc.

\begin{align}
&\frac{(1-\gamma)^3}{\kappa} \mathsf{cost}(\cC_{\hat{\cY}},\cS)\\
&\leq \mathsf{cost}(\cC_{\hat{\cY}},\hat{\cY}_{off}) +  ((1-\gamma)^3\eta_1\Lambda^2 + 16C_M(1-\gamma^2)(1-\gamma))  \cdot \mathsf{cost}(\cC^{opt}_\cS, \cS) \nonumber \\&+ C_M(1-\gamma^2)(1-\gamma) \cdot V'(d,k,\eps,T,\Lambda) ) + (1-\gamma)^2V''(d,k,\eps,T,\Lambda) \\
&\leq \rho \mathsf{cost}(\cC^{opt}_{\hat{\cY}},\hat{\cY}_{off}) +  ((1-\gamma)^3\eta_1\Lambda^2 + 16C_M(1-\gamma^2)(1-\gamma))  \cdot \mathsf{cost}(\cC^{opt}_\cS, \cS) \nonumber \\&+ C_M(1-\gamma^2)(1-\gamma) \cdot V'(d,k,\eps,T,\Lambda) ) + (1-\gamma)^2V''(d,k,\eps,T,\Lambda) \label{eq:rho-cy-sc} \\
&\leq \rho \mathsf{cost}(\cC^{opt}_{\cS},\hat{\cY}_{off}) +  ((1-\gamma)^3\eta_1\Lambda^2 + 16C_M(1-\gamma^2)(1-\gamma))  \cdot \mathsf{cost}(\cC^{opt}_\cS, \cS) \nonumber \\&+ C_M(1-\gamma^2)(1-\gamma) \cdot V'(d,k,\eps,T,\Lambda) ) + (1-\gamma)^2V''(d,k,\eps,T,\Lambda) \label{eq:opt-sc}\\
&\leq \rho ((1+\gamma)^3\kappa \mathsf{cost}(\cC^{opt}_{\cS},\cS) + 16 C_M(1+\gamma)^3\mathsf{cost}(\cC^{opt}_\cS, \cS) + C_M(1+\gamma)^3 \cdot V'(d,k,\eps,T,\Lambda) \nonumber\\&+(1+\gamma)^2V''(d,k,\eps,T,\Lambda)) +  ((1-\gamma)^3\eta_1\Lambda^2 + 16C_M(1-\gamma^2)(1-\gamma))  \cdot \mathsf{cost}(\cC^{opt}_\cS, \cS) \nonumber \\&+ C_M(1-\gamma^2)(1-\gamma) \cdot V'(d,k,\eps,T,\Lambda) ) + (1-\gamma)^2V''(d,k,\eps,T,\Lambda) \label{eq:coreset-app}\\
&\leq ((1+\gamma)^3\rho (\kappa + 16 C_M) +((1-\gamma)^3\eta_1\Lambda^2 + 16C_M(1-\gamma^2)(1-\gamma)))\mathsf{cost}(\cC^{opt}_{\cS},\cS) \nonumber\\&+(\rho C_M(1+\gamma)^3+ C_M(1-\gamma^2)(1-\gamma)) V'(d,k,\eps,T,\Lambda)+((1+\gamma)^2\rho + (1-\gamma)^2) V''(d,k,\eps,T,\Lambda) \label{eq:simplify-sc}
\end{align}

We present the main results from applying our DP clustering framework as final-stemmer and final-ghazi in the sequel. We first recall the existing DP clustering or coreset results that we use as a blackbox. Note that in final-stemmer, we need to convert the set of $k$ centers obtained from applying the DP clustering algorithm in dp-clustering-stemmer into a semicoreset via clustering-to-semicoreset.

There is an $(\eps,\delta)$-DP algorithm that takes a database $S$ containing $n$ points in the $d$-dimensional ball $\cB(0,\Lambda)$, and outputs with probability $1-\beta$, a $(\gamma,\eta)$-approximation for the $k$-means of $S$ where $\gamma=O(1)$ and $\eta={O}((k^{1.01} \cdot d^{0.51} + k^{1.5})\cdot \poly(\log(n),\log(\frac{1}{\beta}),\log(\frac{1}{\delta}),\frac{1}{\eps}) )$

Given dimension $d$, clustering parameter $k$, non-DP $(1+\gamma)$-coreset [18], $(\eps,\delta)$-DP clustering algorithm from dp-clustering-stemmer [42], and $(1+\gamma)$-approx non-DP Clustering algorithm.  Let $\cS$ be the stream of input points in Euclidean space  with length at most $T$. There exists a streaming algorithm $\cA'$ for $k$-means that outputs a set of $k$ centers $\cC_{\hat{\cY}}$ at every timestep $t \in [T]$ such that

* (Privacy) $\cA'$ is $(5\eps,\delta)$-DP under the continual release setting.

* (Accuracy) With probability $1-\frac{1}{T^2}-\frac{1}{k^{O(\poly(k,\log(\Lambda))}}$,

\begin{align}
\mathsf{cost}(\cC_{\hat{\cY}},\cS)
&\leq O_{\gamma}(1) \cdot \mathsf{cost}(\cC^{opt}_\cS, \cS)+O_{\gamma}(1) \cdot V'(d,k,T,\Lambda) + O_{\gamma}(1) \cdot V''(d,k,\eps,\delta, T,\Lambda)
\end{align}

where $\cS$ denotes the set of all points given by the stream before timestamp $t$, $\cC^{opt}_\cS$ is the optimal solution, $V'(d,k,\eps,T,\Lambda) = \tilde{O}(\frac{\Lambda^2 k^2}{d\eps} \poly\log(T))$, $V''(d,k,\eps,\delta,T,\Lambda)= \tilde{O}({k} \Lambda^2 \poly(\log(T),\log(\frac{1}{\delta}),d,\frac{1}{\eps},k))$.

* (Space) $\cA'$ consumes  $\tilde{O}((k^{1.01} \cdot d^{3.51} + k^{1.5}\cdot d^3)\poly(\log(T),\log(\frac{1}{\delta}),d,\frac{1}{\eps}))$

We first note that using clustering-to-semicoreset together with dp-clustering-stemmer gives us an algorithm $\cA$ that produces a $(\kappa,\eta_1,\eta_2)$-semicoreset where $\kappa=O(1), \eta_1=O(1)$ and $\eta_2={O}((k^{1.01} \cdot d^{0.51} + k^{1.5}+k)\cdot \poly(\log(T) \log(\frac{1}{\delta}),\frac{1}{\eps}) )$ and that we can use as a black-box in our DP clustering framework.
The privacy guarantee follows by observing that the resulting algorithm from clustering-to-semicoreset is $(3\eps,\delta)$-DP and using dp-main.
We choose $C_M=100$ and the accuracy guarantee follows from final-accuracy-main.

Recall that $M= \frac{\alpha \eta_2}{C_M}$ (see M) where $\alpha=O(d^3)$ by final-acc-dp-hh, $\eta_2={O}((k^{1.01} \cdot d^{0.51} + k^{1.5})\cdot \poly(\log(T) \log(\frac{1}{\delta}),\frac{1}{\eps}) )$ and $C_M=100$, thus $M= {O}((k^{1.01} \cdot d^{3.51} + d^3 \cdot k^{1.5})\cdot \poly(\log(T), \log(\frac{1}{\delta}),\frac{1}{\eps}) )$. From main-space, the space usage is given by

\begin{align}\label{eq:eq-space}
&\log(\Lambda) \cdot ({S_\cA(M, k,d,\eps,\delta,\kappa,\eta_1,\eta_2)}+ S_\cB(SZ_\cA(M,k,d,\eps,\delta,\kappa,\eta_1,\eta_2),k,d, \gamma)+3M/2) \nonumber\\&+\tilde{O}(k \log(T/M)\log(\Lambda)\log k \cdot \gamma^{-4}) +O({k}\log(k) \log^2(\Lambda) \log T)+O(k \log(\Lambda) \log^2(k) \poly\left( \log \left({T}{\Lambda}k\right)\right))
\end{align}

Note that $S_\cA(M,k,d)= O((M+k)d)+O(k^{1+\eta} \cdot d \cdot \log(M))$ for some arbitrarily chosen constant $\eta>0$ where the second term comes from the space usage of clustering-to-semicoreset, $SZ_\cA(\cdot) = k$ and $S_\cB (SZ_\cA(\cdot))= S_\cB (k) =\tilde{O}(k)$. The final space complexity is obtained by plugging in the value of $M$ and setting $\eta=0.01$.

There is a $2^{O_\gamma(d')}\poly(n)$-time $\eps$-DP algorithm that with probability 0.99 produces an $\left(1+\gamma, O_\gamma\left(\frac{k^2 2^{O_\gamma(d')}}{\eps} \poly\log n\right)\right)$-coreset for $k$-means (and $k$-median) where $d'=O_\gamma(\log k)$. The size of the coreset is $2^{O_\gamma(d')} \cdot \poly(k,\log n)$.

Given dimension $d$, clustering parameter $k$, non-DP $(1+\gamma)$-coreset [18], $\eps$-DP coreset from dp-coreset-ghazi [29], and $(1+\gamma)$-approx non-DP Clustering algorithm.  Let $\cS$ be the stream of input points in Euclidean space  with length at most $T$. There exists a streaming algorithm $\cA'$ for $k$-means that outputs a set of $k$ centers $\cC_{\hat{\cY}}$ at every timestep $t \in [T]$ such that

* (Privacy) $\cA'$ is $3\eps$-DP under the continual release setting.

* (Accuracy) With probability $1-\frac{1}{T^2}-\frac{1}{k^{O(\poly(k,\log(\Lambda))}}$,

\begin{align}
\mathsf{cost}(\cC_{\hat{\cY}},\cS) &\leq (1+\gamma')\mathsf{cost}(\cC^{opt}_\cS, \cS) \nonumber+((1+\gamma)^2 + \frac{(1+\gamma)^6}{(1-\gamma)^5})\cdot V'(d,k,T,\Lambda) \\&+ (\frac{(1+\gamma)(1-\gamma)^2}{(1-\gamma)^3}+\frac{(1+\gamma)^4}{(1-\gamma)^3}) \cdot V''(d,k,\eps,T,\Lambda)
\end{align}

where $\cS$ denotes the set of all points given by the stream before timestamp $t$, $\cC^{opt}_\cS$ is the optimal solution, $\gamma' = 2\gamma + \gamma^2 + \frac{2(1+\gamma)^5}{(1-\gamma)^3} $, $V'(d,k,\eps,T,\Lambda) = \tilde{O}(\frac{\Lambda^2 k^2}{d\eps} \poly\log(T))$, $V''(d,k,\eps,T,\Lambda)= \tilde{O}_\gamma (\frac{{k}^3 \Lambda^2 d^3}{\eps} \poly(\log(T)))$.

* (Space) $\cA'$ consumes $\tilde{O}\left(\poly\left(O_\gamma\left(\frac{d^3k^2 2^{O_\gamma(d')}}{\eps} \poly\log(T) \right) \right),k,d \right)$

First note that the $\eps$-DP coreset from dp-coreset-ghazi [29] is a $(\kappa,0,\eta)$-semicoreset (see semicoreset) where $\kappa=1+\gamma$ and $\eta = O_\gamma\left(\frac{k^2 2^{O_\gamma(d')}}{\eps} \poly\log n\right)$. The privacy guarantee follows from dp-main. The accuracy guarantee follows from final-accuracy-main where we choose $C_M=(1-\gamma)/16$.

The space usage follows from main-space, where we note that $S_\cA(M,d,k)= \poly(M,d,k)$ and total coreset size $SZ_\cA(\cdot) =2^{O_\gamma(d')} \cdot \poly(k,\log T)$, therefore $S_\cB (SZ_\cA(\cdot))= S_\cB (2^{O_\gamma(d')} \cdot \poly(k,\log T)) =\tilde{O}(2^{O_\gamma(d')} \cdot \poly(k,\log T))$. Recall that $M= \frac{\alpha \eta}{C_M}$ (see M) where $\alpha=O(d^3)$ by final-acc-dp-hh, $\eta= O_\gamma\left(\frac{k^2 2^{O_\gamma(d')}}{\eps} \poly\log (T)\right)$ by dp-coreset-ghazi where $d'=O_\gamma(\log k)$ . Thus $M=  O_\gamma\left(\frac{d^3k^2 2^{O_\gamma(d')}}{\eps} \poly\log T \right)$.

§ BICRITERIA APPROXIMATION IN CONTINUAL RELEASE SETTING

We describe our bicriteria approximation algorithm and analysis in more detail here.

[1]
Privacy parameter $\eps$, Stream $\cS$ of points $x_1, \ldots, x_T \in \bbR^d$
$\eps' := \frac{\eps}{4\log(\Lambda) \log^2(k)}$
Parallel quadtrees $Q_1,\ldots, Q_{\log(k)}$ such that:
each quadtree $Q_q$ has $\log(\Lambda)$ levels with the bottom level having grid size $\Theta(1)$

$0 \leq \ell \leq \log(\Lambda)$ and $1 \leq q \leq \log(k)$
Initialize $\mathbf{DPFindCenters}_{\ell,q}$ of $\mathbf{DPFindCenters}(\eps')$
Set of candidate centers $\cF := \emptyset$

$0 \leq \ell \leq \log(\Lambda)$ and $1 \leq q \leq \log(k)$

$\hat{\cF}_t \leftarrow \mathbf{Update}(x_t) \text{ of } \mathbf{DPFindCenters}_{\ell,q}$See dp-find-center-online
$\cF \leftarrow \cF \cup \hat{\cF}_t$
Output $\cF$

[1]
Privacy parameter $\eps'$,
Stream of points $x_1, \ldots, x_T$

Set of candidate centers $\hat{\cF}_t$ at every timestep $t \in [T]$
$\eps' \leftarrow \eps'$ where $\eps' := \frac{\eps}{4\log(\Lambda) \log^2(k)}$
Hash function $h:[2^\ell] \to [w]$ s.t. $\forall \text{ cells } \bfc, \forall j \in [w], \Pr[h(\bfc)=j] = \frac{1}{w}$
$\hat{T}_{1}=0, \ldots, \hat{T}_{w}=0$DP Count for the size of hash bucket
Initialize $\textbf{BM}_1,\ldots,\textbf{BM}_w$ of $\textbf{BinaryMechanism}(T,\eps')$See binarymechanism [22]
Initialize $\textbf{DP-HH}_1,\ldots,\textbf{DP-HH}_w$ of $\textbf{DP-HH}(T,\eps')$See dp-hh-cr [24]
Initialize $\hat{\cF}_t=\emptyset$
Let $x_t$ be mapped to cell $\bfc^*$ Recall DPFindCenters is initialized per level $\ell$ of quadtree instance $q$
$p=1,\ldots,L$, where $L:=\log (k^2)$ run in parallel
$j \in [w]$ Update the DP count of each hash bucket

$j = h(\bfc^*)$
$\hat{T}_{j} \leftarrow \textbf{Update}(1)\text{ of }\textbf{BM}_{j}$
$\hat{T}_{j} \leftarrow \textbf{Update}(0)\text{ of }\textbf{BM}_{j}$
$j\in[w]$ Update the DP HHs of each hash bucket
$\hat{f}, H \leftarrow$ Update($\bfc^*$) of $\textbf{DP-HH}_{j}$
$\hat{f}, H \leftarrow$ Update($\perp$) of $\textbf{DP-HH}_j$
$\bfc \in H$
$\hat{f}(\bfc) \geq \frac{\theta}{1000} \cdot \hat{T}_{h(\bfc)}$
Add centerpoint of $\bfc$ to $\hat{\cF}_t$ as a center
Return $\hat{\cF}_t$

Algorithm. Our bicriteria approximation algorithm is given by dp-cluster-online which initializes $\log(k)$ parallel instances of randomly shifted quadtrees. Each input point $x_t$ is assigned to a cell in every level of every quadtree. For a fixed quadtree $1 \leq q \leq \log(k)$ and fixed level $0 \leq \ell \leq \Lambda$, the subroutine  (see dp-find-center-online) returns a candidate set of centers $\hat{\cF}_t$ which is added to the current set of candidate set of centers $\cF$.

The  subroutine (see dp-find-center-online) finds the approximate heaviest $O(k)$ cells in a fixed level of a fixed quadtree. It achieves this by first hashing the cell containing the current point to a bucket, note that there are $w:= O(k)$ many buckets. For each hash bucket $j \in [w]$, the algorithm maintains a continual release $\theta$-heavy hitter instance $\textsf{DP-HH}_j$.  We use the $\ell_1$-heavy hitter algorithm from [24] as $\textsf{DP-HH}$ — it returns a set $H$ of $\theta$-heavy hitters and their approximate counts $\hat{f}(\bfc)$ for all $\bfc \in H$. Since we are storing the centerpoints of all the cells marked as heavy hitters as candidate centers, we need to ensure that we do not store too many false positives, i.e., cells whose counts are much smaller than $\theta$ times the length of the hashed substream. To address this challenge, we have an additional pruning step that eliminates any cell $\bfc$ whose approximate count is less than  $\Theta(\theta) \hat{T}_{h(\bfc)}$ where $\hat{T}_{h(\bfc)}$ denotes the DP count of each hash bucket $j \in [w]$ at timestep $t \in [T]$. We keep track of $\hat{T}_{j}$ via an instance of the Binary Mechanism [22] denoted as $\textbf{BM}_j$ for each $j\in [w]$. Finally, only the centerpoints of cells that pass this pruning step are added as candidate centers to the set $\hat{\cF}_t$.

§.§ Proof of final-acc-dp-hh

is $\eps$-DP under the continual release setting.

For a fixed timestep $t$, an input point $x_t$ is passed to an instance of  (dp-find-center-online) where it is assigned to a specific cell for a specific level of the quadtree, and cells at the same level are disjoint. The cell containing $x_t$ is then hashed to a bucket which maintains a  instance.

There are $\log (\Lambda)$ levels per quadtree, thus point $x_t$ is a member of $\log(\Lambda)$ cells in total. As there are $2 \log^2(k)$ parallel processes — considering $\log(k)$ quadtrees and $\log(k^2)$ parallel processes (see Line <ref> of dp-find-center-online) per quadtree — a single point participates in $2\log \Lambda \log^2(k)$ total calls to .
Since  is DP under the continual release setting (by dp-hh-cr), assigning each  instance a privacy budget of $\frac{\eps}{4\log(\Lambda) \log^2(k)}$ preserves $(\eps/2)$-DP by sequential composition (Item <ref> of composition).

Next, the algorithm uses the Binary Mechanism (which we know is DP by binarymechanism) with a privacy budget of $\frac{\eps}{4\log(\Lambda) \log^2(k)}$ to keep track of the size of each hash substream $\hat{T}_j$ $\forall j \in [w]$. Since the input cells (and corresponding points within cells) are disjoint in each substream due to hashing, this preserves $\frac{\eps}{4\log(\Lambda) \log^2(k)}$-DP by parallel composition (Item <ref> of composition) which over $2\log(\Lambda) \log^2(k)$ parallel processes preserves $(\eps/2)$-DP by sequential composition (Item <ref> of composition).

Finally, since only the centerpoints of cells identified as heavy hitters by  are released, these centers maintain $\eps$-DP by sequential composition (Item <ref> of composition).

With probability at least $1-1/k^{O(\poly(k,\log(\Lambda))}$,

\begin{align*}
cost(\cF, \cS) &\leq  O(d^3)cost(C^{opt}_{\cS}, \cS) + \tilde{O}\left(\frac{d^{2} \Lambda^2 k}{\eps} \poly\left(\log\left({T\cdot k \cdot\Lambda}\right)\right)\right)
\end{align*}

where $cost(C^{opt}_{\cS}, \cS)$ is the optimal $k$-means cost for $\cS$.

We first state some geometric properties regarding the cells within the quadtree construction.
[17]
Let $\cB$ be an $\ell_\infty$ ball of radius $r$ contained in $[-\Lambda, \Lambda]^d$ (it forms a $d$-dimensional cube with each side length $2r$). Then for a randomly shifted quadtree and any level $\ell$ with grid size at least $r' \geq 2r$, $\cB$ is split by the grid in each dimension $j \in [d]$ independently with probability $\frac{2r}{r'}$.

Let $C^{opt}_\cS =\{c_1, \ldots, c_k\}$ be the optimal set of $k$ centers for the input set $\cS=\{x_1, \ldots, x_T\}$. For any radius, define $n_r$ as the number of points $x \in \cS$ such that $d(x,C^{opt}_\cS)\geq r$.
Note that the opt cost of $k$-means (and $k$-median) is given by $\sum_{p \in \bbZ} 2^{2p} \cdot n_{2^p}$ and $ \sum_{p \in \bbZ} 2^{p} \cdot n_{2^p}$ (up to an $O(1)$-approximation).

Fix some radius $r=2^p$ where $p \in \bbZ$ and consider a randomly shifted grid of size $20rd$. The following lemma characterizes cells containing $\cup^k_{i=1} \cB(c_i,r)$ with respect to the grid size.

[17]
$\cup^k_{i=1} \cB(c_i,r)$ is contained in at most $4k$ cells of grid length $20rd$ by the corresponding level of the quadtree with probability at least 1/2.

Let $G_\ell$ where $0 \leq \ell \leq \log(\Lambda)$ be the set of $4k$ good cells of length $20rd$ (equivalently $\ell_2$-radius of $10rd^{3/2}$) at level $\ell$. Let the number of points in $\cS$ uncovered by $G_\ell$ be $n_{G_\ell}$. Observe that by balls-cells, since $G_\ell$ contains $\cup^k_{i=1} \cB(c_i,r)$ with probability at least 1/2, we have that $n_{G_\ell}\leq n_r$. It follows that

\begin{align}\label{eq:len-radius-opt}
&\sum^{\log(\Lambda)}_{\ell=0}(\text{grid length at level }\ell)^2 \cdot n_{G_\ell} \nonumber\\
&\leq  O(d^{3})\sum_{p \in \bbZ\ :\ r=2^p \leq \Lambda} r^2 \cdot n_r  \leq O(d^3) \cdot  cost(C^{opt}_{\cS}, \cS)
\end{align}

Observe that we can define a one-one mapping between the level $\ell$ and the radius $r$, i.e., the radius $r$ (ranging from $1$ to $\Lambda$) maps to the grid length of a cell which is at most $\Lambda/2^\ell$ (level $\ell$ ranges from $\log (\Lambda)$ to 0). Since the grid length of a cell in $G_\ell$ at level $\ell$ is $20rd$ which maps to $20d \frac{\Lambda}{2^\ell}$, we can replace the leftmost term in the expression above as follows

\begin{align}\label{eq:len-radius-opt-1}
O(d^2)\sum^{\log(\Lambda)}_{\ell=0}(\Lambda/2^\ell)^2  n_{G_\ell} &\leq   O(d^{3})\sum_{p \in \bbZ\ :\ r=2^p \leq \Lambda} r^2 \cdot n_r
\leq O(d^3) \cdot cost(C^{opt}_{\cS}, \cS)
\end{align}

Recall that we define $\cF_t$ as the set of centers until time step $t$. For a fixed level $\ell$, let the set of cells the algorithm  marks as heavy at timestep $t$ at level $\ell$ as $H_{\ell,t}$. Note that although there is an extra pruning step in  after the cells are marked heavy by , we do not account for this here — as if a cell is an $\alpha$-HH and marked heavy by , and it survives the pruning step, it will still be an $\alpha$-HH. Then,

\begin{align*}
cost(\cF_t) \leq O(d^{2})\sum^{\log(\Lambda)}_{\ell =0}(\Lambda/2^\ell)^2 \cdot \mathbbm{1}[x_t \text{ uncovered by }H_{\ell,t}]
\end{align*}

Observe that

\begin{align}\label{eq:cost-f}
&cost(\cF) \nonumber\\
&= \sum^T_{t=1} cost(\cF_t) \nonumber\\
&\leq O(d^{2})\sum^{\log(\Lambda)}_{\ell =0}(\Lambda/2^\ell)^2 \cdot \sum^T_{t=1}\mathbbm{1}[x_t \text{ uncovered by }H_{\ell,t}]
\end{align}

For a fixed level $\ell$, with probability at least $1-\frac{12}{k}$,

\begin{align}\label{eq:hli-upperbd}
\sum^T_{t=1} \mathbbm{1}[x_t \text{ uncovered by }H_{\ell,t}] &\leq  (1+\theta)n_{G_\ell}
+ \frac{4k \log^2 \Lambda \log k}{\eps \eta} \poly\left(\log\left(\frac{T\cdot 2^\ell}{\theta \xi \gamma_h }\right)\right)&
\end{align}

Observe that

\begin{align*}
&\sum^T_{t=1} \mathbbm{1}[x_t \text{ uncovered by }H_{\ell,t}] \nonumber\\
&=\sum^T_{t=1} (\mathbbm{1}[(x_t \text{ uncovered by }H_{\ell,t})\land(x_t \text{ uncovered by }G_{\ell})]\nonumber+\mathbbm{1}[(x_t \text{ uncovered by }H_{\ell,t})\land(x_t \text{ covered by }G_{\ell})]) \nonumber\\
&=\sum^T_{t=1} \mathbbm{1}[(x_t \text{ uncovered by }H_{\ell,t})\land(x_t \text{ uncovered by }G_{\ell})] \nonumber+\sum^T_{t=1} \mathbbm{1}[(x_t \text{ uncovered by }H_{\ell,t})\land(x_t \text{ covered by }G_{\ell})] %\label{eq:sum-uncoverh-coverg}
\end{align*}

The first sum in the above expression can be upper bounded by $n_{G_\ell}$, thus it remains to bound the second sum. In order to bound the second sum, we will need some properties of good cells that are hashed to buckets in , the proofs of these claims can be found in bicriteria. In the sequel, we denote $N_{\ell,\bfc}$ as the number of points in the cell $\bfc$ at level $\ell$ and the hash buckets as $\cB_j$ where $j \in [w]$. For simplicity, we consider the number of hash buckets $w:= 40k$. We first show that for any good cell $\bfc$, it is unlikely that the bucket it is hashed to contains another good cell $\bfc'\neq \bfc$.

Let $\bfc\in G_\ell$, then with probability at least $1/2$, for any $\bfc' \in G_\ell$ such that $\bfc'\neq \bfc$, we have that $ h(\bfc') \neq h(\bfc)$.

In the next claim we give a bound on the size of the hash bucket denoted as $\cB_j$ where $j \in [w]$ in terms of the size of a good cell that is hashed to it and $n_r$.

For each $\bfc \in G_\ell$, suppose the hash bucket $\cB_j$ where $j \in [w]$, contains only one good cell which is $\bfc$. Let $N_{\ell,\bfc}:= y$. Then with probability at least 1/2, $\vert \cB_j \vert \leq 2(y+\frac{n_{G_\ell}}{40k})$.

Note that since the hashing procedure is run $\log(k^2)$ times in parallel, we can boost the success probabilities in the above claims to be $1-1/k^2$.

Observe that for a fixed hash bucket $\cB_j$, any cell $\bfc$ such that $N_{\ell,\bfc} \geq \theta \cdot 2(y+\frac{n_{G_\ell}}{40k})$ qualifies as a $\theta$-heavy hitter since $N_{\ell,\bfc} \geq \theta \cdot 2(y+\frac{n_{G_\ell}}{40k}) \geq \theta \vert \cB_j\vert$ (by hash). In particular, for good cell $\bfc_y$ such that $N_{\ell,\bfc_y} = y$, if $\bfc_y$ is an $\theta$-HH then $y \geq \theta \cdot 2(y+\frac{n_{G_\ell}}{40k})$ or $ y \geq \frac{\theta n_{G_\ell}}{20k}$. We formalize this intuition in the claim below where we use the accuracy guarantees of  given by dp-hh-cr to characterize the good cells that are reported as $\theta$-HHs.

Let $\bfc \in G_\ell$. If $N_{\ell,\bfc} \geq \frac{\theta n_{G_\ell}}{20k}$, and $N_{\ell,\bfc} \geq  \frac{2\log(\Lambda) \log^2(k)}{\eps \gamma_h} \poly( \log(\frac{T\cdot k \cdot 2^\ell}{\theta \gamma_h}))$, then with probability at least $1-\frac{12}{k}$, $\bfc$ is reported as an $\theta$-heavy hitter by .

Finally, we give an upper bound for the number of points that are covered by good cells but for which  fails to report as heavy.

With probability $1-\frac{12}{k}$,

\begin{align*}
&\sum^T_{t=1} \mathbbm{1}[(x_t \text{ uncovered by }H_{\ell,t})\land(x_t \text{ covered by }G_{\ell})] \leq\theta n_{G_\ell} + \frac{8k \log(\Lambda) \log^2(k)}{\eps \gamma_h} \poly\left(\log\left(\frac{T\cdot k \cdot  2^\ell}{\theta \gamma_h }\right)\right)
\end{align*}

Thus by combining uncovered-hh-covered-good with our observation about the first sum being upper bounded by $n_{G_\ell}$ in the decomposition of $\sum^T_{t=1} \mathbbm{1}[x_t \text{ uncovered by }H_{\ell,t}]$, we obtain our desired statement in hh-uncover.

Note that we have shown hh-uncover is true with probability at least $1-\frac{12}{k}$, for a fixed level. Since we have $\log(\Lambda)$ many levels in a specific quadtree, and $\log(k)$ many quadtree instances in parallel — we can boost our probability of success to be sufficiently high. It remains to bound the total $k$-means cost for the set of centers $\cF$ output by our algorithm. Combining len-radius-opt, len-radius-opt-1 and cost-f along with hh-uncover, we obtain the following.

\begin{align*}\small
&cost(\cF, \cS)
\\&= \sum^T_{t=1} cost(\cF_t, \cS) \\
&\leq O(d^{2})\sum^{\log(\Lambda)}_{\ell =0}(\Lambda/2^\ell)^2 \cdot \sum^T_{t=1}\mathbbm{1}[x_t \text{ uncovered by }H_{\ell,t}] \\
&\leq O(d^{2})\sum^{\log(\Lambda)}_{\ell =0}(\Lambda/2^\ell)^2 \cdot ((1+\theta)n_{G_\ell}+  \frac{8k \log(\Lambda) \log^2(k)}{\eps \gamma_h} \poly\bigl(\log\bigl(\frac{T\cdot k \cdot 2^\ell}{\theta  \gamma_h }\bigr)\bigr)) \\
&= O(d^{2})\sum^{\log(\Lambda)}_{\ell =0}(\Lambda/2^\ell)^2 \cdot (1+\theta)n_{G_\ell}+O(d^{2})\sum^{\log(\Lambda)}_{\ell =0}(\Lambda/2^\ell)^2 \cdot  \frac{8k \log(\Lambda) \log^2(k)}{\eps \gamma_h} \poly\bigl(\log\bigl(\frac{T\cdot k \cdot 2^\ell}{\theta \gamma_h }\bigr)\bigr) \\
&\leq O(d^{3})(1+\theta) \sum_{p \in \bbZ\ :\ r=2^p \leq \Lambda} r^2 \cdot n_r + O(d^{2})\sum^{\log(\Lambda)}_{\ell =0}(\Lambda/2^\ell)^2 \cdot \frac{8k \log(\Lambda) \log^2(k)}{\eps \gamma_h} \poly\bigl(\log\bigl(\frac{T\cdot k \cdot 2^\ell}{\theta \gamma_h }\bigr)\bigr) \\
&\leq O(d^3)(1+\theta) cost(C^{opt}_\cS,\cS) + O(d^{2})\sum^{\log(\Lambda)}_{\ell =0}(\Lambda/2^\ell)^2 \cdot  \frac{8k \log(\Lambda) \log^2(k)}{\eps \gamma_h}\poly\bigl( \log\bigl(\frac{T\cdot k \cdot 2^\ell}{\theta  \gamma_h }\bigr)\bigr)\\
&=O(d^3)(1+\theta) cost(C^{opt}_\cS,\cS) + O\left(d^{2} \Lambda^2 \frac{k \log(\Lambda) \log^2(k)}{\eps \gamma_h} \poly\bigl(\log\bigl(\frac{T\cdot k \cdot  \Lambda}{\theta \gamma_h}\bigr)\bigr)\right)
\end{align*}

Finally, we can set $\theta$ (threshold for HHs) and $\gamma_h$ (approximation factor for frequency of a cell marked as heavy from dp-hh-cr) to appropriate constants. The accuracy claim follows.

uses $O(k \log(\Lambda) \log^2(k) \poly\left( \log \left({T}{\Lambda}k\right)\right))$ space.

We analyze the total space usage for  in dp-find-center-online as this dominates space usage for the entire algorithm. From dp-hh-cr, one instance of  uses $\poly\left( \log \left({T}{\Lambda}k\right)\right)$. Since we run  on $O(k)$ many hash substreams, and $2\log(\Lambda) \log^2(k)$ parallel processes, the total space is $O(k \log(\Lambda) \log^2(k) \poly\left( \log \left({T}{\Lambda}k\right)\right))$.

[Size of $\cF$ in final-acc-dp-hh]
For all $j\in [w]$, suppose the size of the hash bucket $\cB_j$ is $\vert \cB_j \vert = \Omega( \frac{\log(\Lambda) \log^3(k) }{\eps} \log^{2.5}T )$ then with high probability, the total number of historical heavy hitters at the end of the stream is $O({k}\log^2(k) \log(\Lambda) \log T)$.

In other words, $\cF$ has atmost $O({k}\log^2(k) \log(\Lambda) \log T)$ centers.

The algorithm runs independent instances of the  algorithm for each bucket of each level in each instantiation of the quadtree, thus it is sufficient to first show that for a fixed quadtree $Q$, a fixed level $\ell$, and a fixed bucket $\cB_j$ where $j\in [w]$, the total number of historical heavy hitters is at most $O(\frac{(1+\gamma_h)}{{\theta}}\log T)$.

Let the timestamps of points that end up in $\cB_j$ be $t_i=2^{i}$, where $0 \leq i \leq \log(T)$. Let the state of the hash bucket at time step $t$ be $\cB^{(t)}_j$. We set the failure probability in binarymechanism as $\xi:=\frac{1}{k^2}$. From binarymechanism we know that with probability $1-\frac{1}{k^2}$, the DP count of the hash bucket $\hat{T}_{j}$ at timestep $t$ has additive error $O(\frac{\log(\Lambda) \log^3(k)}{\eps} \log^{2.5}(T))$. Thus for a fixed timestamp $t$, if $\vert \cB^{(t)}_j \vert = \Omega(\frac{\log(\Lambda)\log^3(k)}{\eps} \log^{2.5}(T) )$, then we can see that a cell $\bfc$ is added to $\cF_t$ only if $\hat{f}(\bfc) > \frac{2\theta}{1000} \vert \cB^{(t)}_j \vert$. Recall from Condition 1 of dp-hh-cr that $\hat{N}_{\ell,\bfc} \in (1\pm \gamma_h) N_{\ell, \bfc}$. Thus if $\bfc \in \cF_t$ and $\vert \cB^{(t)}_j \vert = \Omega(\frac{\log(\Lambda)\log^3(k)}{\eps} \log^{2.5}(T) )$ then it must be the case that with probability $1-\frac{1}{k^2}$, $N_{\ell, \bfc} \geq \frac{{2\theta}}{1000(1+\gamma_h)}\| \cB^{(i)}_j \|_1 \geq \frac{{2\theta}}{1000(1+\gamma_h)} {t_{i-1}}$.

Now, suppose for a contradiction, that the number of heavy hitters between $t_{i-1}$ and $t_i$ is at least $\frac{1000(1+\gamma_h)}{{\theta}}$. Then for each such cell $\bfc$, we have that $N_{\ell,\bfc}\geq \frac{{2\theta}}{1000(1+\gamma_h)} {t_{i-2}}$. Since there are at least $\frac{2000(1+\gamma_h)}{{2\theta}}$ such cells, this implies that the total number of points between $t_{i-1}$ and $t_i$ is $\geq \frac{2\theta}{1000(1+\gamma_h)} {t_{i-2}}\frac{2000(1+\gamma_h)}{{2\theta}} = 2^i = t_i$, which is a contradiction. Thus there must be at most $\frac{1000(1+\gamma_h)}{{\theta}}$ cells marked as heavy hitters between consecutive intervals, and since there are $\log(T)$ such intervals, we have that the total number of historical $\ell_1$ heavy hitters for a fixed bucket is $ O(\frac{(1+\gamma_h)}{{\theta}}\log T)$.

Boosting the success probability over $O(k)$ buckets, $\log(k^2)$ parallel processes, $\log(\Lambda)$ quadtree levels, and $\log(k)$ parallel processes of the quadtree instantiation, accounting for the additional number of historical HHs, and taking $\theta$ and $\gamma_h$ as appropriate constants, we obtain the claim as stated.

§ DP MERGE AND REDUCE ALGORITHM

We give a differentially-private variant of the widely-known Merge and Reduce framework [33, 1, 27] that is used to efficiently release a coreset for a stream of points. The main idea behind the Merge and Reduce technique is to partition the input stream into blocks, compute a coreset for each block, take the union of the resulting coresets (merge step), and compute the coreset of the union (reduce step). The merging and reducing of the coresets is done in a tree-like fashion. In order to reduce the error introduced by merging, the number of levels in the tree must be small. On a high-level, our framework computes coresets at the base level (of the tree) using a DP semicoreset Algorithm $\cA$ (e.g. [42, 29]) and then computes coresets for subsequent levels using a non-DP coreset Algorithm $\cB$ (e.g. [19]).

First, we show that the semicoreset definition (i.e., semicoreset) satisfies the Merge and Reduce properties, i.e., the union of semicoresets is a semicoreset and the coreset of a union of semicoresets is a valid semicoreset for the underlying points.

Algorithm $\mathsf{DP\text{-}Merge\text{-}Reduce}$

[1]

Point $x'_t$ to be added to coreset, DP parameters $\eps,\eps_1,\delta_1$

Block size $M$Value of $M$ depends on our choice of DP semicoreset Algorithm. See final-stemmer and final-ghazi for exact values
Size of $P_0$ as ${p}_0=0$
We only store $P_0$ in the actual algorithm, the rest of the $P_i$'s are only used in the accuracy analysis, so we do not run the following line in the actual algorithm. $N_r$ denotes the number of points in ring $R_r$.
$P_i=\emptyset$ for all $i=1, \ldots, u$ where $u= \lfloor \log(2{N}_{r}/M)\rfloor+1$
$\gamma_i:= \gamma/Ci^2$ for $i=1, \ldots, u-1$

$x'_t \neq \perp$

Insert $x'_t$ into $P_0$
$p_0 = p_0+1$update the size of $P_0$

$\textbf{LevelZero-AboveThreshold}(p_0, M, \eps) \to \top$ See dp-above-threshold

$Q_1 \leftarrow \cA(P_0,k,d, \eps,\delta,\kappa , \eta_1,\eta_2)$ $\cA$ computes a $(\eps_1,\delta_1)$-DP $(\kappa, \eta_1,\eta_2)$-semicoreset, e.g., [42],[29]
$P_1 = P_1 \cup P_0$
Delete elements in $P_0$, and set ${p}_0=0$
$\exists\ Q'_i$ such that $Q_i$ and $Q'_i$ are both at level $i$
Merge-Reduce computes $(1+\gamma_{i})$-coreset $Q_{i+1}$ of $Q_i\cup Q'_i$ and deletes $Q_i,Q'_i$
$Q_{i+1}=\textbf{Merge\text{-}Reduce}(\cB, Q_i,Q'_i,\gamma_i)$ $\cB$ computes a non-DP $(1+\gamma_i)$-coreset, e.g., [19]
$P_{i+1}=P_{i+1}\cup P_i$ Lines <ref> and <ref> are for analysis purpose only
Delete elements in $P_i$
return $\textbf{Merge\text{-}Reduce}(\cB, \{Q_i\}_{ i \leq u}, \gamma/3)$

return $\emptyset$
return $\emptyset$

Algorithm $\textsf{LevelZero-AboveThreshold}$

[1]
Size of set $P_0$ denoted as $p_0$, Block size $M$, Privacy parameter $\eps$
$\top$ if the size of set $P_0$ exceeds noisy threshold and $\perp$ otherwise
$\hat{M} = M+\Lap(2/\eps)$
$\nu = \Lap(4/\eps)$
$p_0+{\nu} \geq \hat{M}$

return $\top$
return $\bot$

* (Merge) If $Q$ is a $(1+\gamma,\eta_1,\eta_2)$-semicoreset of $P$, $Q'$ is a $(1+\gamma,\eta_1,\eta_2)$-semicoreset of $P'$ and $P,P'$ are disjoint, then $Q \cup Q'$ is a $(1+\gamma, \eta_1,2\eta_2)$-semicoreset of $P \cup P'$.

* (Reduce) If $R$ is a $(1+\gamma)$-coreset of $Q\cup Q'$, then $R$ is a $((1+\gamma)^2,(1-\gamma)\eta_1, (1+\gamma)2\eta_2)$-semicoreset of $P$.

We first prove the merge property. For any set of $k$-centers $C \subseteq \bbR^d$,

\begin{align*}
\textsf{cost}(C,Q\cup Q') &\leq \textsf{cost}(C,Q) + \textsf{cost}(C,Q') \\
&\leq (1+\gamma)\textsf{cost}(C,P) + \eta_2 + (1+\gamma)\textsf{cost}(C,P') + \eta_2 \\
&\leq (1+\gamma)\textsf{cost}(C,P \cup P') + 2\eta_2
\end{align*}

\begin{align*}
\textsf{cost}(C,Q\cup Q') &\geq \textsf{cost}(C,Q) + \textsf{cost}(C,Q') \\
&\geq (1-\gamma)\textsf{cost}(C,P) - \eta_1 \textsf{cost}(C^{opt}_P,P)-\eta_2 + (1-\gamma)\textsf{cost}(C,P') - \eta_1 \textsf{cost}(C^{opt}_{P'},P')-\eta_2 \\
&\geq (1-\gamma)\textsf{cost}(C,P \cup P') - \eta_1(\textsf{cost}(C^{opt}_P,P)+\textsf{cost}(C^{opt}_{P'},P'))- 2\eta_2\\
&\geq (1-\gamma)\textsf{cost}(C,P \cup P') - \eta_1\textsf{cost}(C^{opt}_{P\cup P'},P\cup P')- 2\eta_2
\end{align*}

Next we prove the Reduce property below. For any set of $k$-centers $C \subseteq \bbR^d$,

\begin{align*}
\textsf{cost}(C,R) &\leq (1+\gamma)\textsf{cost}(C,Q \cup Q') \\
&\leq (1+\gamma)((1+\gamma)\textsf{cost}(C,P \cup P') + 2\eta_2 )\\
&\leq (1+\gamma)^2\textsf{cost}(C,P \cup P') + (1+\gamma)2\eta_2  \\
%   &\leq (1+\gamma)^2\textsf{cost}(P \cup P') + (1+\gamma)2\eta\\
\end{align*}

\begin{align*}
\textsf{cost}(C,R) &\geq (1-\gamma)\textsf{cost}(C,Q \cup Q') \\
&\geq (1-\gamma)( (1-\gamma)\textsf{cost}(C,P \cup P') - \eta_1\textsf{cost}(C^{opt}_P,P\cup P')- 2\eta_2 )\\
&\geq (1-\gamma)^2\textsf{cost}(C,P \cup P') - (1-\gamma)\eta_1\textsf{cost}(C^{opt}_P,P\cup P')-(1-\gamma)2\eta_2  \\
%   &\leq (1+\gamma)^2\textsf{cost}(P \cup P') + (1+\gamma)2\eta\\
\end{align*}

Observe that at any timestep the total set of input points seen so far (denoted as $P$) is partitioned into subsets $P_0,\ldots, P_u$ where some $P_i$'s can be empty and $u=\lfloor \log(2{N}/M)\rfloor+1$. Note that we simulate this step of partitioning $P$ into $P_i$'s in dp-merge-reduce solely for the analysis. It is not necessary to store $P_1, \ldots, P_u$ explicitly in the actual algorithm.

We first prove some claims about .

[Privacy]  is $\eps$-DP under the continual release setting.

Recall that  checks whether the size of the set $P_0$ is above a certain threshold. In other words, it checks whether the total count of the group of elements $x_i \in P_0$ is above the given threshold. Once a positive response is returned by , elements in $P_0$ are deleted and the process is repeated with a new group of elements. This algorithm is equivalent to grouping a stream of counts in [24]. In particular, the proof of privacy is identical except that in our algorithm we do not release the actual noisy counts but instead we just release a positive/negative response to the same query.

For all $t\in [T]$, with probability $1-\xi$, we have that

* $\vert \nu \vert < \frac{4}{\eps}\log(\frac{2T}{\xi})$

* $\vert \hat{M}-M \vert < \frac{2}{\eps} \log(\frac{2T}{\xi})$

This follows from standard application of tail bounds for Laplace distribution and union bound over all $t \in [T]$.

If $p_0 + \nu \geq \hat{M}$ then with probability $1-\xi$, we have that $p_0 \geq M/2$.

We can simplify $p_0 + \nu \geq \hat{M}$ by applying the noise bounds from noise-above-threshold. Thus with probability at least $1-\xi$, we have that $p_0 \geq {M} - \frac{6}{\eps} \log(\frac{2T}{\xi})$. Finally using the assumption about $M>\frac{12}{\eps}\log(\frac{2T}{\xi})$, the statement follows.

The number of levels is given by $u=\lceil \log(2N/M)\rceil+1$.

By p0-range, $p_0 \geq M/2$, thus the total number of blocks at level 0 is $\leq \frac{2{N}}{M}$. The statement follows.

§.§ Proof of dp-merge-reduce

framework is $(\eps_1+\eps,\delta_1)$-DP under the continual release setting.

First, by dp-above-threshold,  is $\eps$-DP under the continual release setting. Since the semicoresets computed by $\cA$ at level 1 are $(\eps_1,\delta_1)$-DP, we can release these DP semicoresets as they are computed. Subsequent computations on these DP semicoresets via non-DP algorithm $\cB$ preserve DP by postprocessing. Thus by sequential composition (Item <ref> of composition) the claim follows.

With probability at least $1-\xi-\xi_A-\xi_B$,  framework releases a $((1+\gamma)\kappa,(1-\gamma)\eta_1, (\frac{4N}{M}-1)(1+\gamma)\eta_2+\tilde{M})$-semicoreset of $P$. Where $\tilde{M}:= {M}+\frac{6}{\eps}\log(\frac{2T}{\xi})$.

Recall that $P$ is partitioned by $P_1,\ldots, P_u$. We first prove the following claim about the coreset for a non-empty subset $P_r \subseteq P$.

Suppose $P_r$ is non-empty. Then $Q_r$ is a $((1+\gamma/3)\kappa,(1-\gamma/3)\eta_1,(1+\gamma/3)2^{r-1}\eta_2 )$-semicoreset of $P_r$.

We will first prove the claim that $Q_r$ is a $(\prod^{r-1}_{j=1}(1+\gamma_j)\kappa,\prod^{r-1}_{j=1}(1-\gamma_j)\eta_1, \prod^{r-1}_{j=1}(1+\gamma_j)2^{r-1}\eta_2 )$-semicoreset for $P_r$ where $r\geq 2$ by induction. Note that $Q_0=P_0$ for $p_0 + \nu < \hat{M}$. For $p_0 + \nu \geq \hat{M}$, $P_1 = P_1 \cup P_0$. Since we apply DP semicoreset algorithm $\cA$ to $P_1$, the resulting coreset $Q_1$ is a $(\kappa,\eta_1,\eta_2)$-semicoreset for ${P}_1$.

Base Case.
By merge-reduce-property, $Q_1 \cup Q'_1$ is a $(\kappa, \eta_1, 2\eta_2)$-semicoreset for $P_1 \cup P'_1$ and $Q_2$ is a $(\kappa(1+\gamma_1), (1+\gamma_1)\eta_1, (1+\gamma_1)2\eta_2)$-semicoreset of $Q_1\cup Q'_1$.  Here we use the notation $Q'_i$ or $P'_i$ to differentiate between sets or coresets at the same level $i$.

Inductive Hypothesis. Suppose the claim is true for $r=i$, i.e., $Q_i$ is a $(\kappa \prod^{i-1}_{j=1}(1+\gamma_j),\prod^{i-1}_{j=1}(1-\gamma_j)\eta_1,\prod^{i-1}_{j=1}(1+\gamma_j)2^{i-1}\eta_2 )$-semicoreset for $P_i$ and $Q'_i$ is a $(\kappa \prod^{i-1}_{j=1}(1+\gamma_j),\prod^{i-1}_{j=1}(1-\gamma_j)\eta_1, \prod^{i-1}_{j=1}(1+\gamma_j)2^{i-1}\eta_2 )$-semicoreset for $P'_i$.

By merge-reduce-property, the Merge step implies

\begin{align}
\mathsf{cost}(C,Q_i \cup Q'_i) &\leq \kappa\prod^{i-1}_{j=1}(1+\gamma_j)\mathsf{cost}(C,P_i \cup P'_i)+\prod^{i-1}_{j=1}(1+\gamma_j)2^{i} \eta_2 \\
\mathsf{cost}(C,Q_i \cup Q'_i) &\geq \kappa\prod^{i-1}_{j=1}(1-\gamma_j)\mathsf{cost}(C,P_i \cup P'_i) - \prod^{i-1}_{j=1}(1-\gamma_j)\eta_1\mathsf{cost}(C^{opt}_{P_i \cup P'_i},P_i \cup P'_i) - \prod^{i-1}_{j=1}(1-\gamma_j)2^{i} \eta_2
\end{align}

The Reduce step implies that the resulting $(1+\gamma_{i+1})$-coreset of $Q_i \cup Q'_i$ denoted as $Q_{i+1}$ is such that

\begin{align}
\mathsf{cost}(Q_{i+1}) &\leq \kappa\prod^{i}_{j=1}(1+\gamma_j)\mathsf{cost}(P_{i+1})+\prod^{i}_{j=1}(1+\gamma_j)2^{i}\eta_2\\
\mathsf{cost}(Q_{i+1}) &\geq \kappa\prod^{i}_{j=1}(1-\gamma_j)\mathsf{cost}(P_{i+1})- \prod^{i-1}_{j=1}(1-\gamma_j)\eta_1\mathsf{cost}(C^{opt}_{P_{i+1}},P_{i+1}) -\prod^{i}_{j=1}(1-\gamma_j)2^{i}\eta_2
\end{align}

Finally, provided $c$ is large enough, we have:

\begin{align}\label{eq:gamma-gamma'}
&\prod^i_{j=1} (1+\gamma_j) \leq \prod^i_{j=1} \exp(\frac{\gamma}{c j^2}) = \exp(\frac{\gamma}{c}\sum^i_{j=1} \frac{1}{j^2}) \leq  \exp(\frac{\gamma}{c}\cdot \frac{\pi^2}{6}) \leq 1+\gamma/3  \\
&\prod^i_{j=1} (1-\gamma_j) \leq \prod^i_{j=1} \exp(-\frac{\gamma}{c j^2}) = \exp(-\frac{\gamma}{c}\sum^i_{j=1} \frac{1}{j^2}) \leq  \exp(-\frac{\gamma}{c}\cdot \frac{\pi^2}{6}) \leq 1-\gamma/3
\end{align}

The statement follows.

Finally, we release a $(1+\gamma/3)$-coreset of $\cup_{i \leq u}Q_i$ of $\cup_{i\leq u}P_i$ for all non-empty $P_i$ which by similar arguments as above is a $((1+\gamma)\kappa, (1-\gamma)\eta_1, (1+\gamma)\eta_2 (2^{u}-1))$-semicoreset of $P$. The statement follows by plugging in the value for the total number of levels $u$.

Note that if ${p}_0+\nu < \hat{M}$ then we do not release anything. So we have to account for an additional additive error of ${M}+\frac{6}{\eps}\log(\frac{2T}{\xi})$ w.p. $1-\xi$ in this case. If ${p}_0+\nu \geq \hat{M}$, then we proceed by computing a DP coreset using Algorithm $\cA$ with failure probability $\xi_A$. We also compute all the coresets past the first level using Algorithm $\cB$ and a failure probability of $\xi_B/2u$, where $u$ is the number of levels. Thus for a fixed run of , by a union bound, the total failure probability for this part is at most $\xi_B$.

framework uses $S_\cA(M, k,d,\eps,\delta,\kappa,\eta_1,\eta_2, \xi_A)+ S_\cB(SZ_\cA(M,k,d,\eps,\delta,\kappa,\eta_1,\eta_2, \xi_A),k,d, \gamma, \xi_B)+\lceil \log(2N/M) \rceil \cdot S_\cB(SZ_\cB(M,k,d,\gamma, \xi_B),k,d, \gamma, \xi_B) +3M/2$ space.

First we need an upper bound for the size of the blocks at level 0, i.e., $p_0$ which is given by the contrapositive statement of the claim below.

If $p_0 \geq \frac{3M}{2}$ then with probability $1-\xi$, ${p}_0+\nu \geq \hat{M}$.

By applying the noise bounds from noise-above-threshold and using the assumption that $p_0 \geq \frac{3M}{2}$, we have that with probability $1-\xi$,

\begin{align*}
p_0+\nu \geq \frac{3M}{2} - \frac{4}{\eps}\log(\frac{2T}{\xi}) > \frac{3M}{2} - \frac{M}{3} = M+\frac{M}{6} > M+\frac{2}{\eps}\log(\frac{2T}{\xi}) > \hat{M}
\end{align*}

Thus, we only need to store at most
$3M/2$ points for the block of input points at level 0, plus the additional space required to execute the coreset construction. Note that the semicoreset computation of level 0 blocks consumes space $S_\cA(M, k,d,\eps,\delta,\kappa,\eta_1,\eta_2)$ and since the largest coreset construction wrt non-DP algorithm $\cB$ is the union of at most $u$ coresets that is reduced to a single coreset, and the largest input to non-DP coreset algorithm is the resulting DP semicoreset size — the additional storage size is at most $S_\cB(SZ_\cA(M,k,d,\eps,\delta,\kappa,\eta_1,\eta_2),k,d, \gamma)+(u-1) \cdot S_\cB(SZ_\cB(M,k,d,\gamma),k,d, \gamma)$. Note that the resulting coreset size for non-DP coreset algorithm $\cB$ is independent of the input set size.

The resulting semicoreset has size at most $\tilde{O}(k\log (k) \cdot \gamma^{-4})$.

Since $Q_1 \leftarrow \cA(P_0,k,d, \eps,\delta,\kappa , \eta)$, size of $Q_1$ is $SZ_\cA(M,k,d,\eps,\delta,\kappa,\eta_1,\eta_2)$. Now $Q_2, \ldots, Q_u$ are obtained by running non-DP algorithm $\cB$. In order to simplify our notation, we invoke the state of the art non-DP coreset algorithm as $\cB$ (see nondp-coreset). Thus the semicoreset size of $Q_i$ where $2 \leq i \leq u$ is $\tilde{O}(k\log k \cdot \gamma^{-4})$. Finally we take the union of the semicoresets $\cup_{1 \leq i \leq u}Q_i$ which has size at most  $SZ_\cA(M,k,d,\eps,\delta,\kappa,\eta_1,\eta_2)+ (u-1)O(k\log k \cdot \gamma^{-4})$ and then apply non-DP algorithm $\cB$ to obtain a semicoreset of size at most $\tilde{O}(k\log k \cdot \gamma^{-4})$.

[1]
Pankaj K. Agarwal, Sariel Har-Peled, and Kasturi R. Varadarajan.
Approximating extent measures of points.
J. ACM, 51(4):606–635, 2004.

[2]
Sara Ahmadian, Ashkan Norouzi-Fard, Ola Svensson, and Justin Ward.
Better guarantees for k-means and euclidean k-median by primal-dual
SIAM Journal on Computing, 49, 2020.

[3]
Noga Alon, Yossi Matias, and Mario Szegedy.
The space complexity of approximating the frequency moments.
Journal of Computer and system sciences, 58(1):137–147, 1999.

[4]
Vijay Arya, Naveen Garg, Rohit Khandekar, Adam Meyerson, Kamesh Munagala, and
Vinayaka Pandit.
Local search heuristic for k-median and facility location problems.
In STOC, 2001.

[5]
Pranjal Awasthi, Avrim Blum, and Or Sheffet.
Stability yields a ptas for k-median and k-means clustering.
In FOCS, 2010.

[6]
Maria-Florina Balcan, Travis Dick, Yingyu Liang, Wenlong Mou, and Hongyang
Differentially private clustering in high-dimensional euclidean
In Proceedings of the 34th International Conference on Machine
Learning, ICML 2017, Sydney, NSW, Australia, 6-11 August 2017, volume 70
of Proceedings of Machine Learning Research, pages 322–331. PMLR,

[7]
Jeremiah Blocki, Elena Grigorescu, and Tamalika Mukherjee.
Differentially-private sublinear-time clustering.
In 2021 IEEE International Symposium on Information Theory
(ISIT), pages 332–337. IEEE, 2021.

[8]
Vladimir Braverman, Dan Feldman, Harry Lang, and Daniela Rus.
Streaming coreset constructions for m-estimators.
In Approximation, Randomization, and Combinatorial Optimization.
Algorithms and Techniques (APPROX/RANDOM 2019). Schloss
Dagstuhl-Leibniz-Zentrum fuer Informatik, 2019.

[9]
T.-H. Hubert Chan, Elaine Shi, and Dawn Song.
Private and continual release of statistics.
ACM Trans. Inf. Syst. Secur., 14(3):26:1–26:24, 2011.

[10]
Alisa Chang, Badih Ghazi, Ravi Kumar, and Pasin Manurangsi.
Locally private k-means in one round.
In International Conference on Machine Learning, pages
1441–1451. PMLR, 2021.

[11]
Moses Charikar, Sudipto Guha, Éva Tardos, and David B Shmoys.
A constant-factor approximation algorithm for the k-median problem.
Journal of Computer and System Sciences, 65, 2002.

[12]
Ke Chen.
On k-median clustering in high dimensions.
In SODA, 2006.

[13]
Ke Chen.
A constant factor approximation algorithm for k-median clustering
with outliers.
In SODA, 2008.

[14]
Ke Chen.
On coresets for k-median and k-means clustering in metric and
euclidean spaces and their applications.
SIAM J. Comput., 39(3):923–947, 2009.

[15]
Edith Cohen, Haim Kaplan, Yishay Mansour, Uri Stemmer, and Eliad Tsfadia.
Differentially-private clustering of easy instances.
In Marina Meila and Tong Zhang, editors, Proceedings of the 38th
International Conference on Machine Learning, volume 139 of Proceedings
of Machine Learning Research, pages 2049–2059. PMLR, 18–24 Jul 2021.

[16]
Vincent Cohen-Addad, Alessandro Epasto, Silvio Lattanzi, Vahab Mirrokni, Andres
Munoz Medina, David Saulpic, Chris Schwiegelshohn, and Sergei Vassilvitskii.
Scalable differentially private clustering via hierarchically
separated trees.
In KDD, page 221–230, New York, NY, USA, 2022. Association
for Computing Machinery.

[17]
Vincent Cohen-Addad, Alessandro Epasto, Vahab Mirrokni, Shyam Narayanan, and
Peilin Zhong.
Near-optimal private and scalable $ k $-clustering.
In Advances in Neural Information Processing Systems, 2022.

[18]
Vincent Cohen-Addad, Kasper Green Larsen, David Saulpic, and Chris
Towards optimal lower bounds for k-median and k-means coresets.
In Proceedings of the 54th Annual ACM SIGACT Symposium on Theory
of Computing, pages 1038–1051, 2022.

[19]
Vincent Cohen-Addad, David Saulpic, and Chris Schwiegelshohn.
A new coreset framework for clustering.
In STOC '21: 53rd Annual ACM SIGACT Symposium on Theory of
Computing, Virtual Event, Italy, June 21-25, 2021, pages 169–182. ACM,

[20]
Vincent Cohen-Addad, David P. Woodruff, and Samson Zhou.
Streaming euclidean k-median and k-means to a (1 + eps)-approximation
with $o_{k,eps}$(log n) memory words.
In FOCS, 2023.

[21]
Cynthia Dwork, Frank McSherry, Kobbi Nissim, and Adam D. Smith.
Calibrating noise to sensitivity in private data analysis.
J. Priv. Confidentiality, 7, 2016.

[22]
Cynthia Dwork, Moni Naor, Toniann Pitassi, and Guy N. Rothblum.
Differential privacy under continual observation.
In Leonard J. Schulman, editor, Proceedings of the 42nd ACM
Symposium on Theory of Computing, STOC 2010, Cambridge, Massachusetts, USA,
5-8 June 2010, pages 715–724. ACM, 2010.

[23]
Cynthia Dwork and Aaron Roth.
The algorithmic foundations of differential privacy.
Found. Trends Theor. Comput. Sci., 9(3-4):211–407, 2014.

[24]
Alessandro Epasto, Jieming Mao, Andres Muñoz Medina, Vahab Mirrokni,
Sergei Vassilvitskii, and Peilin Zhong.
Differentially private continual releases of streaming frequency
moment estimations.
CoRR, abs/2301.05605, 2023.

[25]
Dan Feldman, Amos Fiat, Haim Kaplan, and Kobbi Nissim.
Private coresets.
In Proceedings of the 41st Annual ACM Symposium on Theory of
Computing, STOC 2009, Bethesda, MD, USA, May 31 - June 2, 2009, pages
361–370. ACM, 2009.

[26]
Dan Feldman and Michael Langberg.
A unified framework for approximating and clustering data.
In Proceedings of the forty-third annual ACM symposium on Theory
of computing, pages 569–578, 2011.

[27]
Dan Feldman, Melanie Schmidt, and Christian Sohler.
Turning big data into tiny data: Constant-size coresets for k-means,
pca, and projective clustering.
SIAM Journal on Computing, 49(3):601–657, 2020.

[28]
Dan Feldman, Chongyuan Xiang, Ruihao Zhu, and Daniela Rus.
Coresets for differentially private k-means clustering and
applications to privacy in mobile sensor networks.
In Proceedings of the 16th ACM/IEEE International Conference
on Information Processing in Sensor Networks, IPSN 2017, Pittsburgh, PA,
USA, April 18-21, 2017, 2017.

[29]
Badih Ghazi, Ravi Kumar, and Pasin Manurangsi.
Differentially private clustering: Tight approximation ratios.
In Advances in Neural Information Processing Systems 33: Annual
Conference on Neural Information Processing Systems 2020, NeurIPS 2020,
December 6-12, 2020, virtual, 2020.

[30]
Anupam Gupta, Katrina Ligett, Frank McSherry, Aaron Roth, and Kunal Talwar.
Differentially private combinatorial optimization.
In Proceedings of the Twenty-First Annual ACM-SIAM Symposium
on Discrete Algorithms, SODA 2010, Austin, Texas, USA, January 17-19,
2010, pages 1106–1125. SIAM, 2010.

[31]
Sariel Har-Peled and Akash Kushal.
Smaller coresets for k-median and k-means clustering.
In Proceedings of the twenty-first annual symposium on
Computational geometry, pages 126–134, 2005.

[32]
Sariel Har-Peled and Soham Mazumdar.
On coresets for k-means and k-median clustering.
In Proceedings of the thirty-sixth annual ACM symposium on
Theory of computing, pages 291–300, 2004.

[33]
Sariel Har-Peled and Soham Mazumdar.
On coresets for k-means and k-median clustering.
In Proceedings of the 36th Annual ACM Symposium on Theory of
Computing, Chicago, IL, USA, June 13-16, 2004, pages 291–300. ACM, 2004.

[34]
Zhiyi Huang and Jinyan Liu.
Optimal differentially private algorithms for k-means clustering.
In Jan Van den Bussche and Marcelo Arenas, editors, Proceedings
of the 37th ACM SIGMOD-SIGACT-SIGAI Symposium on Principles of Database
Systems, Houston, TX, USA, June 10-15, 2018, pages 395–408. ACM, 2018.

[35]
Max Dupré la Tour, Monika Henzinger, and David Saulpic.
Differential privacy for clustering under continual observation,

[36]
Shi Li and Ola Svensson.
Approximating k-median via pseudo-approximation.
SIAM Journal on Computing, 45, 2016.

[37]
Frank D McSherry.
Privacy integrated queries: an extensible platform for
privacy-preserving data analysis.
In Proceedings of the 2009 ACM SIGMOD International Conference
on Management of data, pages 19–30, 2009.

[38]
Kobbi Nissim, Sofya Raskhodnikova, and Adam D. Smith.
Smooth sensitivity and sampling in private data analysis.
In David S. Johnson and Uriel Feige, editors, Proceedings of the
39th Annual ACM Symposium on Theory of Computing, San Diego, California,
USA, June 11-13, 2007, pages 75–84. ACM, 2007.

[39]
Kobbi Nissim and Uri Stemmer.
Clustering algorithms for the centralized and local models.
In Algorithmic Learning Theory, ALT 2018, 7-9 April 2018,
Lanzarote, Canary Islands, Spain, volume 83 of Proceedings of Machine
Learning Research, pages 619–653. PMLR, 2018.

[40]
Rafail Ostrovsky, Yuval Rabani, Leonard J Schulman, and Chaitanya Swamy.
The effectiveness of lloyd-type methods for the k-means problem.
Journal of the ACM (JACM), 59, 2012.

[41]
David Saulpic, Vincent Cohen-Addad, and Andreas Emil Feldmann.
Near-linear time approximations schemes for clustering in doubling
In David Zuckerman, editor, 60th IEEE Annual Symposium on
Foundations of Computer Science, FOCS 2019, Baltimore, Maryland, USA,
November 9-12, 2019, pages 540–559. IEEE Computer Society, 2019.

[42]
Uri Stemmer and Haim Kaplan.
Differentially private k-means with constant multiplicative error.
In Advances in Neural Information Processing Systems 31: Annual
Conference on Neural Information Processing Systems 2018, NeurIPS 2018,
December 3-8, 2018, Montréal, Canada, 2018.

[43]
Chang Xia, Jingyu Hua, Wei Tong, and Sheng Zhong.
Distributed k-means clustering guaranteeing local differential
Computers & Security, 90:101699, 2020.

[44]
Matei Zaharia, Tathagata Das, Haoyuan Li, Scott Shenker, and Ion Stoica.
Discretized streams: An efficient and $\{$Fault-Tolerant$\}$ model
for stream processing on large clusters.
In 4th USENIX Workshop on Hot Topics in Cloud Computing
(HotCloud 12), 2012.

§ PROOFS FOR BICRITERIA APPROXIMATION

Proof of good-cell-collide.

For any good cell $\bfc' \neq \bfc$, define
$$X_{\bfc'}= \begin{cases} 1 & \text{if }\bfc'\text{ collides with }\bfc\\ 0 &\text{otherwise}\end{cases}$$

Now for a fixed $\bfc' \neq \bfc$, the expected number of collisions with $\bfc$ is given by $\E_h[X_{\bfc'}] = \Pr[X_{\bfc'}=1] =\Pr[h(\bfc)=h(\bfc')] = \frac{1}{w} = \frac{1}{40k}$. Thus the total number of collisions for $\bfc'$ with $\bfc$, in expectation, is given by $\E_h[\sum_{\bfc'} X_{\bfc'}] \leq 4k/w = 1/10 \leq 1/2$.

Therefore, by Markov, $\Pr[\sum_{\bfc'} X_{\bfc'}>1] \leq 1/2$. The claim follows.

Proof of hash.

\begin{align}
\E_h \vert \cB_j \vert &= \E_h [ y + \sum_{\bfc \not\in G_\ell} N_{\ell, \bfc}]\\
&=y + \E_h[\sum_{\bfc \not\in G_\ell} N_{\ell, \bfc}]\\
&=y+ \frac{n_{G_\ell}}{40k}
\end{align}

The claim follows by Markov inequality.

Proof of good-hh.

We set the failure probability of  (see dp-hh-cr) to be $\xi:=\frac{1}{k^2}$.
For $\bfc \in G_\ell$, define $\cE_\bfc := \cE_1 \wedge \cE_2 \wedge \cE_3$ where $\cE_1$ is the event that  algorithm is correct on all instances, $\cE_2$ is the event that there are no collisions between $\bfc$ and other good cells, and $\cE_3$ is the event that there exists a hash bucket that contains only $\bfc$, and if $N_{\ell,\bfc}:= y$, then the size of that hash bucket is $\leq 2(y+\frac{n_{G_\ell}}{40k})$. We know that $\Pr[\cE_1] \geq 1-\xi= 1-\frac{1}{k^2}$,
by the accuracy guarantee of  algorithm, $\Pr[\cE_2] \geq 1-\frac{1}{k^2}$ by good-cell-collide (and boosting the success probability), and $\Pr[\cE_3] \geq 1-\frac{1}{k^2}$, by hash (and boosting the success probability). Thus by a union bound, we have that for a fixed $\bfc \in G_\ell$, $\Pr[\cE_\bfc] \geq 1-\frac{3}{k^2}$.
By taking a union bound over all $4k$ good cells, with probability at least $1-\frac{12}{k}$, the claim holds.

Proof of uncovered-hh-covered-good.

Recall that $N^{(t)}_{\ell,\bfc}$ is the number of points in cell $\bfc$ at time step $t$. Because $x_t$ is covered in $G_\ell$, this means that we only need to care about the points covered by a good cell $\bfc \in G_\ell$ but $\bfc \not\in H_{\ell,t}$. Using good-hh, we know that if $N^{(t)}_{\ell,\bfc} \geq \frac{\theta n_{G_\ell}}{20k} $, and $N^{(t)}_{\ell,\bfc} \geq  \frac{2\log(\Lambda) \log^2(k)}{\eps \eta} \poly(\log(\frac{T\cdot k \cdot 2^\ell}{\theta \eta }))$, then with probability at least $1-\frac{12}{k}$, we have that $\bfc \in H_{\ell,t}$.

If $\bfc \in G_\ell$, and $\bfc \not \in H_{\ell,t}$, this means that either (1) $N^{(t)}_{\ell,\bfc} <  \frac{2\log(\Lambda) \log^2(k)}{\eps \eta} \poly(\log(\frac{T\cdot k \cdot 2^\ell}{\theta \eta }))$, or, (2) $\bfc$ is not an $\theta$-HH. For case (1), since there are $4k$ such good cells, the total number of uncovered points in such cells are $\frac{4k \cdot 2\log(\Lambda) \log^2(k)}{\eps \eta} \poly(\log(\frac{T\cdot k \cdot 2^\ell}{\theta  \eta }))$.
For case (2), this means that $N^{(t)}_{\ell,\bfc} < \frac{\theta n_{G_\ell}}{20k} $.
Again, since there are $4k$ such good cells, the total number of uncovered points in such cells are $< \frac{\theta n_{G_\ell}}{20k} \cdot 4k < \theta n_{G_\ell} $.