Improved Coresets and Sublinear Algorithms for Power Means in Euclidean Spaces
Year: 2022

Abstract:
In this paper, we consider the problem of finding high dimensional power means: given a set A of n points in R, find the point m that minimizes the sum of Euclidean distance, raised to the power z, over all input points. Special cases of problem include the well-known Fermat-Weber problem – or geometric median problem – where z = 1, the mean or centroid where z = 2, and the Minimum Enclosing Ball problem, where z =∞. We consider these problem in the big data regime. Here, we are interested in sampling as few points as possible such that we can accurately estimate m. More specifically, we consider sublinear algorithms as well as coresets for these problems. Sublinear algorithms have a random query access to the set A and the goal is to minimize the number of queries. Here, we show that Õ ( ε−z−3 ) samples are sufficient to achieve a (1+ε)-approximation, generalizing the results from Cohen, Lee, Miller, Pachocki, and Sidford [STOC ’16] and Inaba, Katoh, and Imai [SoCG ’94] to arbitrary z. Moreover, we show that this bound is nearly optimal, as any algorithm requires at least Ω ( ε−z+1 ) queries to achieve said approximation. The second contribution are coresets for these problems, where we aim to find find a small, weighted subset of the points which approximates cost of every candidate point c ∈ R up to a (1 ± ε) factor. Here, we show that Õ ( ε−2 ) points are sufficient, improving on the Õ ( dε−2 ) bound by Feldman and Langberg [STOC ’11] and the Õ ( ε−4 ) bound by Braverman, Jiang, Krauthgamer, and Wu [SODA 21].

text: We consider these problem in the big data regime. Here, we are interested in sampling as few points as possible such that we can accurately estimate m. More specifically, we consider sublinear algorithms as well as coresets for these problems. Sublinear algorithms have a random query access to the set A and the goal is to minimize the number of queries. Here, we show that Õ ( ε−z−3 ) samples are
sufficient to achieve a (1+ε)-approximation, generalizing the results from Cohen, Lee, Miller, Pachocki, and Sidford [STOC ’16] and Inaba, Katoh, and Imai [SoCG ’94] to arbitrary z. Moreover, we show that this bound is nearly optimal, as any algorithm requires at least Ω ( ε−z+1 ) queries to achieve said approximation.
The second contribution are coresets for these problems, where we aim to find find a small, weighted subset of the points which approximates cost of every candidate point c ∈ Rd up to a (1 ± ε) factor. Here, we show that Õ ( ε−2 )
points are sufficient, improving on the Õ ( dε−2 ) bound by Feldman and Langberg [STOC ’11]
and the Õ ( ε−4 ) bound by Braverman, Jiang, Krauthgamer, and Wu [SODA 21].

heading: 1 Introduction
text: Large data sets have shifted the focus of algorithm design. In the past, an algorithm might have been deemed feasible if its running time was polynomial in the input size and so a textbook fast algorithm can have time complexity for example quadratic. For truly gargantuan data sets, even linear time or nearly linear time algorithms could be considered too slow or requiring too much memory. This led to the emergence of the field of sublinear algorithms: How well can we solve a problem without reading the entire input?
Except for trivial problems, deterministic time sublinear algorithms do not exist. Our primary tool in designing sublinear algorithms is thus the following basic approach:
• Take a uniform sample of the input.
• Run an algorithm on the sample.
Hence, the performance of a sublinear algorithm is often measured in terms of its query complexity, i.e. the number of samples required such that we can extract a high quality solution in the second
∗Equal contribution.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
step above that generalizes to the entire input. Sublinear algorithms have close ties to questions in learning theory and estimation theory, where we are similarly interested in a quality to sample size tradeoff.
A perhaps very fundamental problem of primary importance in machine learning and data analysis is to efficiently estimate the parameters of a distribution. For example, given a distribution D, how many samples do we need to estimate the mean? Even such a simple and basic question has surprisingly involved answers and are still subject to ongoing research (Lugosi & Mendelson (2019); Lee & Valiant (2021)).
In this paper, we investigate the possibility of estimating power means in high dimensional Euclidean spaces. Specifically, given an arbitrary set of points A, we wish to determine the number of uniform queries S such that we can extract a power mean m with
cost(m) := ∑ p∈A ‖p−m‖z ≤ (1 + ε) · min µ ∑ p∈A ‖p− µ‖z,
where ‖p‖ denotes the Euclidean norm of a vector p. The power mean problem captures a number of important problems in computational geometry and multivariate statistics. For example, for z = 1, this corresponds to the Fermat-Weber problem also known as the geometric median. For z = 2, the problem is to determine the mean or centroid of the data set. Letting z →∞, we have the Minimum Enclosing Ball (MEB), where one needs to find the Euclidean sphere of smallest radius containing all input points.
For z > 2, the problem is not as well studied, but it still has many applications. First, higher powers allows us to interpolate between z = 2 and z →∞, which is interesting as the latter admits no sublinear algorithms2. Skewness (a measure of the asymmetry of the probability distribution of a real-valued random variable about its mean) and kurtosis (a measure of the ”tailedness” of the probability distribution) are the centralized moments with respect to the three and the four norms and are frequently used in statistics. The power mean is a way of estimating these values for multivariate distributions.
Another application is when dealing with non-Euclidean distances, such as the Hamming metric, coreset constructions for powers of z can be reduced to coreset constructions for powers 2z. So for example if we want the mean in Hamming space, we can reduce it to the z = 4 case /in squared Euclidean spaces Huang & Vishnoi (2020).
These problems are convex and thus can be approximated in the near-linear time efficiently via convex optimization techniques. However, aside from the mean (z = 2), doing so in a sublinear setting is challenging and to the best of our knowledge, only the mean and the geometric median (z = 1) are currently known to admit nearly linear time algorithms.
Our main result is:
Theorem 1. There exists an algorithm that, with query complexity O(ε−z−3 · polylog(ε−1) log2 1/δ), computes a (1 + ε) approximate solution to the high dimensional power mean problem with probability at least 1− δ. A key component in designing this algorithm is a novel analysis for coresets for these problems. Coresets are succinct summaries that approximate cost for any center solution c.
Theorem 2. For any set of points in d-dimensional Euclidean space, there exists an ε-coreset for the high dimensional power mean problem of size Õ(ε−2 · 2O(z)).
With the exception of the mean, previous coresets for these problems achieved boundsO(ε−2 ·2O(z) · min(d, ε−2)) (Cohen-Addad et al. (2021); Braverman et al. (2021); Feldman & Langberg (2011)), or had weaker guarantees such as merely approximating an optimal solution or requiring removal of outliers from the data set.
Comparing the bounds in Theorem 1 and Theorem 2, one may question whether the exponential dependency in z is necessary for computing an approximation. Indeed, previous sublinear algorithms
2To see this, we place n− 1 points at 0 and one point with probability 1/2 at 1 and with probability 1/2 at 0. Any 2 − ε approximation can distinguish between the two cases, but this clearly requires us to query Ω(n) points.
for both the geometric median and the mean had a query complexity of O(ε−2), and thus matched our coresets bounds. Unfortunately, we show that the exponential dependency in the power is indeed necessary even in a single dimension:
Theorem 3. For any ε > 0 and z, any algorithm that computes with probability more than 4/5 a (1 + ε)-approximation for a one-dimensional power mean has query complexity Ω(ε−z+1).
Hence, up to constants in the exponent, our sublinear algorithm is tight. Moreover, the algorithm is very simple to implement and performs well empircally.

heading: 1.1 Techniques
text: While stochastic gradient descent has been used for a variety of center-based problems (Clarkson et al. (2012); Cohen et al. (2016)), it is difficult to apply it for higher powers. Indeed, Cohen et al. (2016) remark in their paper that even for the mean3 (z = 2) their analysis does not work as the objective function is neither Lipschitz, nor strictly convex.
Hence, one needs to use new tools. A natural starting point is to use techniques from coresets, as they allow us to preserve most of the relevant information, using a substantially smaller number of points. Unfortunately, coresets have a drawback: the sampling distributions used to construct coresets is non-uniform and therefore difficult to use in a sublinear setting. The first step is therefore to design coresets from uniform sampling. To do this, we use and improve upon a technique originally introduced by Chen (2009). Chen showed that, given a sufficiently good initial solution q, one can partition the points into rings exponentially increasing radii such that the points cost the same, up to constant factors. Thereafter, taking a uniform sample of size Õ(d · ε−2) from each ring produces a coreset. Since there are at most O(log n) rings in the worst case, this yields a coreset of size Õ(d · ε−2 · log n). To realize these ideas in a sublinear setting, we are now faced with a number of challenges. First, rings that are particularly far from q may contain very few points. This makes is difficult for a sublinear algorithm to access them. Second, partitioning the points into rings depends on the cost of q. It is very simple to construct examples where estimating the cost of an optimal power center requires Ω(n) many queries. Finally, this analysis loses both factors log n and d, which we aim to avoid.
We improve and extend this framework in two ways. First, we show that it is sufficient to only consider O(log ε−1) many rings, which, in of itself, already removes the dependency on log n. Moreover, we show that it is possible to simply ignore any ring containing too few points, i.e. any ring with less than εz+O(1) · n points may be discarded. The intuition is that, while rings with few point may contribute significantly to the cost, these points do not influence the position of the optimal center by much. Thus, using a number of carefully chosen pruning steps, we show how to reduce both the problem of obtaining a sublinear algorithm as well as obtaining coresets to sampling from a select few rings containing many points.
The second improvement directly considers an improved analysis of sampling from rings. The standard way to prove a coreset guarantee is to show that using s log 1/δ samples we preserve the cost of a single solution with probability 1 − delta. If there exist T solutions then we set δ = 1/T and have obtained a coreset, which we call the ”naive” union bound. This works in certain cases such as finite metrics, but is insufficient if we have infinitely many solutions such as Euclidean spaces. The simplest way to improve over the naive union bound is to discretize the space and then apply the naive union bound on the discritization. In literature this is sometimes called an ε-net bound. This can be optimal or close to optimal for certain metrics, but so far these arguments have only lead to Õ(ε−2 min(d, ε−2)) bounds for coresets in Euclidean spaces.
Instead of applying a union bound over the discretization in ”one shot”, we apply a union bounded over a nested sequence of increasingly better discretizations. Essentially, instead of only using a set of centers Cε as a substitute for all centers is Rd, we use centers ch ∈ Ch for different values of h ∈ {1, . . . , log 1/ε}. In literature, this is known as chaining, see Nelson (2016) for a survey. We can then write, for any input point set P , ∑ p∈P cost(p, c) = ∑ p∈P ∑ h≥0 cost(p, c
h+1) − cost(p, ch), where cost(p, c0) is defined to be 0 and |cost(p, ch) − cost(p, c)| ≤ 2−hcost(p, c). We now only
3For the special case of the mean, Inaba et al. Inaba et al. (1994) observed that O(ε−2) samples are nevertheless sufficient.
apply the naive union bound for successive summands, i.e. we approximate ∑ p∈P cost(p, c
h+1)− cost(p, ch) up to an error of ε · ∑ p∈P cost(p, c).
Essentially, the idea is to use a sequence of solutions ch that approximate a candidate solution c with increasing accuracy as h → ∞. Approximating the cost of a candidate center c can then be written as a telescoping sum of solutions ch, i.e. cost(c) = ∑∞ h=0 cost(ch+1) − cost(ch), where ch is a solution that has distance to c of the order 2−h and cost(c0) := 0. The notion of distance between candidate solutions is perhaps most easily understood by imagining a solution to be a cost vector vc where the ith entry corresponds to the cost of the ith point of A when using c as a candidate center, i.e., vci = cost(pi, c). Informally, we consider two solutions c and c
′ to have distance at most 2−h if |‖pi − c‖z − ‖pi − c′‖z| ≤ 2−h · (‖pi − c‖z + ‖pi − c′‖z) for all points pi ∈ A.

heading: 1.2 Related Work
text: Sublinear Approximation for Clustering: A number of sublinear algorithm are known for clustering problems. For k-Median, under the constraint that the input space has a small diameter, a constant factor approximation is known Czumaj & Sohler (2007). Ben-David (2007) proposed a different set of conditions under which a sublinear algorithm for k-median and k-means exists. Other approximations, with different constraint are also known: for instance, Meyerson et al. (2004) give an algorithm achieving a O(1)-approximation in time poly(k/ε) for discrete metrics, when each cluster has size Ω(nε/k). For the 1-median problem, this assumption is always satisfied, and their algorithm gives a constant factor approximation in constant running time. The algorithm by Cohen et al. (2016) produces a (1 + ε)-approximation in time min(nd log3(n/ε), d/ε2) for Euclidean spaces of dimension d. Ding (2020, 2021) and Clarkson et al. (2012) showed how to obtain sublinear algorithms for the minimum enclosing ball problem assuming that either the algorithm is allowed to drop a fraction of the points, or with an additive error. For the unconstrained version of the problem, no sublinear time algorithm is possible. For the k-means problem, Bachem et al. (2016) showed how to approximate the k-means++ algorithm in sublinear time. To our knowledge, no algorithm is known for higher distance powers z.
Coreset Constructions: A coreset is a weighted subset of the input points, such that the cost of any clustering is the same on the coreset than on the input points, up to a (1 ± ε) factor. The arguably most widely studied problem is coresets for the k clustering problems with powers of distances. Following a long line of work Bachem et al. (2018); Becchetti et al. (2019); Braverman et al. (2021); Chen (2009); Cohen-Addad et al. (2021); Feldman & Langberg (2011); Feldman et al. (2020); Feng et al. (2021); Fichtenberger et al. (2013); Har-Peled & Kushal (2007); Har-Peled & Mazumdar (2004); Huang & Vishnoi (2020); Langberg & Schulman (2010); Sohler & Woodruff (2018), we now know that coreset of size Õ(kε−4 · 2O(z) ·min(k, ε−max z−2,0))) exist. In a low-dimensional regime, this may be improved to Õ(k2dε−22O(z)) and Õ(kdε−2z) Feldman & Langberg (2011). For the special case of clustering with a single center, this yields the state of the art Õ(ε−4) bound due to Braverman et al. (2021). For the Minimum Enclosing Ball problem and its generalizations, these algorithms yield no space saving, although sketches of weaker gaurantees of size O(ε−1) exist (Badoiu & Clarkson (2008)). Given a coreset, it is easy to compute a (1 + ε)-approximation in time independent of n: one just needs to find a (1 + ε)-approximation on the coreset points.
Coresets have also been studied for many other problems: we cite non-comprehensively fair clustering Cohen-Addad & Li (2019); Huang et al. (2019); Schmidt et al. (2019) determinant maximization Indyk et al. (2020), diversity maximization Ceccarello et al. (2018); Indyk et al. (2014) logistic regression Huggins et al. (2016); Munteanu et al. (2018), dependency networks Molina et al. (2018), or low-rank approximation Maalouf et al. (2019). The interested reader is referred to the recent surveys for more information Feldman (2020); Munteanu & Schwiegelshohn (2018).

heading: 1.3 Organization
text: We prove our key sampling result in Section 2. We follow up on that result by applying it to the sublinear setting (Section 3) and to coresets (Section 4). We conclude with a short experimental evaluation. Due to space constraints, all proofs are included in the supplementary material.

heading: 1.4 Preliminaries
text: We denote the Euclidean distance of a vector x by ‖x‖ := √∑d
i=1 x 2 i . Similarly, we define the Hamming norm ‖x‖1 = ∑d i=1 |xi|. For a set of points A, we say that ‖A‖0 is the distinct number
of points. Note that if A contains multiplicities |A| 6= ‖A‖0. We write Õ(x) to denote O(x · loga x), where a is any constant. For any set A and candidate solution c, we defined cost(A, c) =∑ p∈A ‖p − c‖z . If the set is clear from context, we simply write cost(c). We frequently use the following generalized triangle inequality, see Cohen-Addad & Schwiegelshohn (2017); Makarychev et al. (2019) for proofs and similar statements.
Lemma 1 (Triangle Inequality for Powers). Let a, b, c be an arbitrary set of points in a metric space with distance function d and let z be a positive integer. Then for any ε > 0
d(a, b)z ≤ (1 + ε)z−1d(a, c)z + ( 1 + ε
ε
)z−1 d(b, c)z
|d(a, b)z − d(a, c)z| ≤ ε · d(a, c)z + ( z + ε
ε
)z−1 d(b, c)z.
We also use the fact that uniform sampling is efficient to approximate the density: we review details on sampling in bounded VC dimension in the supplementary material.
Lemma 2 (Li et al. (2001)). Given a range space (X,R) with VC-dimension d, an constants ε, δ, η, and a uniform sample S ⊂ X of size at least O( 1η·ε−2 (d log 1/η+ log 1/δ)), we have for all ranges R ∈ R with |X ∩R| ≥ η · |X|∣∣∣∣ |X ∩R||X| − |S ∩R||S|
∣∣∣∣ ≤ ε · |X ∩R||X| and for all ranges R ∈ R with |X ∩R| ≤ η · |X|∣∣∣∣ |X ∩R||X| − |S ∩R||S|
∣∣∣∣ ≤ ε · η with probability at least 1− δ. The only range space we will consider is the one induced by Euclidean spheres centered around a single fixed point. The VC dimension induced of this range space is 2, which seems to be a well known fact, although we could not find a reference. For completeness, we added a short proof in the supplementary material.
First, we recall the commonly used coreset definition for clustering problems in Euclidean spaces.
Definition 1. Let A be a set of points in Rd. Then a set Ω is a strong (ε, z)-coreset if there exists a weight function w : Ω→ R+ and a constant κ such that for any point c∣∣∣∣∣∣ ∑ p∈A ‖p− c‖z − ( ∑ p∈S wp‖p− c‖z + κ) ∣∣∣∣∣∣ ≤ ε · ∑ p∈A′ ‖p− c‖z.
We say that Ω is a weak (ε, z)-coreset if for some α ∈ [0, 1] any point satisfying ∑ p∈Ω wp · ‖p − c′‖z ≤ (1 + α · ε)argmin c∈Rd ∑ p∈Ω wp · ‖p − c‖z also satisfies ∑ p∈A ‖p − c′‖z ≤ (1 + ε)argmin c∈Rd ∑ p∈A ‖p− c‖z.
The difference between the two notions is that strong coresets give a guarantee for all candidate centers, whereas the weak coreset guarantee only applies for the optimum. In an offline setting, there is no difference in the size for our construction for either guarantee. In the sublinear setting, we will be satisfied with a weak coreset as it can be obtained with a nearly optimal query complexity.

heading: 2 Uniform Sampling Routine
text: In this section, we outline the proof of our basic sampling subroutine. We assume that we are given a point q, and a set of points R such that for all p, p′ ∈ R
‖p− q‖z ≤ 2 · ‖p′ − q‖z.
In the following sections, we refer to R as a ring. Under this assumption, the following claim holds.
Theorem 4. Let q and R be defined as above and let Ω be a uniform random sample consisting of s ∈ Õ(ε−2 · log 1/δ) points. Then with probability at least 1− δ, we have for all candidate centers c ∣∣∣∣∣∣cost(R, c)− ∑ p∈Ω |R| s · ‖p− c‖z
∣∣∣∣∣∣ ≤ ε · (cost(R, q) + cost(R, c)). We prove this theorem, first by proving the following slightly simpler result and then showing how to apply the result recursively to remove dependency in ‖R‖0:
Lemma 3. Let q and A be defined as above and let Ω be a uniform random sample consisting of s ∈ Õ(ε−2 ·log ‖R‖0 ·log 1/δ) points. Then with probability at least 1−δ, we have for all candidate centers c that ∣∣∣cost(R, c)−∑p∈Ω |R|s · ‖p− c‖z∣∣∣ ≤ ε · (cost(R, q) + cost(R, c)). To set up a chaining analysis, we require two ingredients: (1) a notion of nets and (2) a Gaussian process. We focus on the latter first. For any point c, let vc be the |R|-dimensional vector, henceforth called a cost vector, such that vpi = ‖pi − c‖z for some arbitrary but fixed ordering of the points in R. Note that ‖v‖1 = cost(R, c). Let pj ∈ Ω with j ∈ {1, . . . , s} be the jth point of the sample. Observe for any cost vector v induced by some center c, we have
EΩ [n s vpj ] = ∑ p∈R ‖p− c‖z = cost(R, c).
We now symmetrize the expectation. Let g1, . . . gs be standard Gaussian random variables, i.e. Gaussians with mean 0 and variance 1. Then we have for any collection of cost vectors N (see Appendix B.3 of Rudra & Wootters (2014))
EΩsup v∈N ∣∣∣∣∣ ∑s j=1 |R| s vpj − ‖v‖1 cost(R, q) + ‖v‖1 ∣∣∣∣∣ ≤ √2πEΩEgsupv∈N ∣∣∣∣∣ ∑s j=1 |R| s vpj · gj cost(R, q) + ‖v‖1 ∣∣∣∣∣ (1) Note that the first term in Equation 1 is the expected maximum deviation from the (normalized) expected cost of the sample, if the cost vectors are induced by centers. In other words, if
sup c∈Rd ∣∣∣∣∑sj=1 |R|s ‖pj−c‖z−cost(R,c)cost(R,q)+cost(R,c) ∣∣∣∣ ≤ ε, we have the desired coreset guarantee. Our cost vectors will not be induced by centers for technical reasons, but are in a well-defined sense close enough such that it will be enough to bound the deviation for these approximate cost vectors to obtain a bound for all center induced cost vectors. Introducing Gaussians is standard in this type of analysis as concentration bounds typically used for weighted Bernoulli random variables such as Hoeffding or
Bernstein are too weak. Our goal is therefore to show that EΩEg sup c∈Rd ∣∣∣∣∑sj=1 |R|s ‖pj−c‖z·gjcost(R,q)+cost(R,c) ∣∣∣∣ ≤ ε√2π . We now define the nets we will use.
Definition 2. Let R be a set of points and let q be a candidate solution. For β > 0, we say that a set of vectors Nβ is a β-net of R, if there exists a vector v′ ∈ Nβ such that for every point p ∈ R with ‖p− c‖ ≤ 8zε · ‖p− q‖, we have∣∣‖p− c‖z − v′p∣∣ ≤ β · (‖p− c‖z + ‖p− q‖z) . We need to show the following three properties:
1. By bounding EΩ sup v∈Nε/10 ∣∣∣∣∑sj=1 |R|s vpj−∑p∈Ri vpcost(R,q)+cost(R,c) ∣∣∣∣, we can also bound EΩ sup
c∈Rd ∣∣∣∣∑sj=1 |R|s ‖pj−c‖z−cost(R,c)cost(R,q)+cost(R,c) ∣∣∣∣ 2. There exist β-nets of size exp ( z2 log ‖R‖0 · β−2 · log 1ε·β ) .
3. Let v be a cost vector written as a telescoping sum v = ∑ i=1∞vi − vi−1 of cost vectors from 2i-nets. Let v′, v′′ be cost vectors from successive β and β/2 nets Nβ and Nβ/2, respectively. By bounding the term EΩEg sup
v′−v′′ ∣∣∣∑sj=1 |R|s |v′pj − v′′pj | · gj∣∣∣, we can also bound EΩEg sup
v∈Nε ∣∣∣∑sj=1 |R|s vpj · gj∣∣∣ We start with item 1 via the following lemma. Lemma 4. Let Nε/10 be an ε/10-net of R. Then if supv∈Nε/10 |∑sj=1 |R|s vpj−‖v‖1| cost(R,q)+cost(R,c) ≤ ε 10 we have
supc∈Rd |∑sj=1 |R|s ‖pj−c‖z−cost(R,c)| cost(R,q)+cost(R,c) ≤ ε for all c ∈ R d.
To prove item 2, we use terminal embeddings Elkin et al. (2017); Mahabadi et al. (2018); Narayanan & Nelson (2019), defined as follows.
Definition 3 (Terminal Embeddings). Let A be a set of points in Rd. A mapping f : Rd → Rk is a terminal embedding if for all p ∈ A and all c ∈ Rd
(1− ε) · ‖p− c‖2 ≤ ‖f(p)− f(c)‖2 ≤ (1 + ε) · ‖p− c‖2.
The guarantee of a terminal embedding is very similar to the guarantee of the famous Johnson Lindenstrauss lemma, but stronger in one crucial detail. A terminal embedding preserves not only the distances between the points of A but also the distance between an arbitrary point in Rd and any point of A. Despite this stronger guarantee, the target dimension of terminal embedding is in fact no worse than that of the Johnson Lindenstrauss lemma. Specifically:
Theorem 5 (Narayanan & Nelson (2019)). For any n point-set A ⊂ Rd, there exists a terminal embedding f : Rd → Rk with k ∈ γ · ε−2 log n for some absolute constant γ.
We now use the terminal embeddings to show that small nets exist.
Lemma 5. Let R and q be defined as above. Then for every β > 0, there exists a β-net of R of size at most exp(γ · z3β2 log ‖R‖0 · log ε−1), where γ is an absolute constant.
We now move onto item 3.
Lemma 6. Let R and q be defined as above and let Ω be a uniform sample consisting of s points. Then for any point c and s ≥ η·z328z ·ε−2 ·log ‖R‖0 ·log3 ε−1 for some absolute constant η, we have
EΩEg sup c∈Rd ∣∣∣∣∑sj=1 |R|s ‖pj−c‖z·gjcost(R,q)+cost(R,c) ∣∣∣∣ ≤ ε. Moreover, if s ≥ η · z328z · ε−2 · log ‖R‖0 · log4 ε−1 log 1/δ for some absolute constant η, then we have sup
c∈Rd ∣∣∣∣∑sj=1 |R|s ‖pj−c‖z·gjcost(R,q)+cost(R,c) ∣∣∣∣ ≤ ε, with probability at least 1− δ.
The proof of Lemma 3 is now a direct consquence of Lemma 6 and Equation 1. To finally prove Theorem 4, we apply Lemma 6 recursively. The result is essentially a special case of Theorem 3.1 from Braverman et al. (2021).
Lemma 7. Let R and q be defined as above. Suppose a uniform sample of size s ∈ Õ(Γ · log ‖R‖0 · log 1/δ) satisfies with probability at least 1− δ for all candidate centers c
|cost(R, c)− ∑ p∈Ω |R| s · ‖p− c‖z| ≤ ε · (cost(R, q) + cost(R, c)).
Then a uniform sample of size Õ(Γ · log 1/δ) achieves the same guarantee.

heading: 3 Sublinear Algorithm and Analysis
text: We first describe an algorithm that succeeds with constant probability. This probability can be amplified (non-trivially) using independent repetition. In the following we will use parameters β and η that depends on ε which we will specify later. We let A be the set of input points.
Algorithm 1 Sublinear Algorithm for Power Means 1. Sample a random point q. 2. Sample a set S of O(ε−z−3 · poly(ε−1)) points uniformly at random. 3. Compute the maximum distance d such that there exist 2/3 · ε · η · |S| points with distance at least d from q. Discard all points at distance greater than d. 4. Define rings Ri such that Ri ∩ S contains all the points at distance (d · 2−i, d · 2−i+1] from q, with i = {1, . . . , β}. 5. If |Ri ∩ S| < εη · |S|, remove all points in Ri ∩ S from S. 6. Define R̂i = n · |Ri∩S||S| . Weigh the points Ri ∩ S by R̂i |Ri∩S| .
7. Solve the problem on the (weighted) set S.
Our goal is to show that S satisfies a weak coreset guarantee. Specifically, we will show that S has the property that for all points A′ ⊆ A at distance at most d (defined in Algorithm 1) from q there exists weight function w : S → R+ such that S is a strong coreset for A′. We then show that a strong coreset for A′ is also weak coreset for A, i.e. we preserve the cost of all points in A′, and we preserve the optimum for A.
In order to call Theorem 4 on only few rings, we now show that the loss incurred by only considering A′ is indeed negligible.
Pruning Lemmas We first show that we can safely discard points that are sufficiently far away, as parameterized by γ. We recall the meaning of parameters: α is the approximation-factor of the initial solution, β is such that points at distance closer than 2βd can be merged to the center (see Lemma 9), η is such that rings with less than an εη-fraction of the points can be discarded (see Lemma 10)
Lemma 8. Suppose we are given an α-approximate center q. Let B(q, r) be the ball centered at q with radius r = 4 · ( 2α·OPT
n
)1/z . Then the following two statements hold.
1. Any α-approximate center is in B(q, r).
2. For any two points c, c′ ∈ B(q, r) and for any point p with ‖p − q‖z > γ · rz with γ > ε−z · (12z)z , we have ‖p− c‖z ≤ (1 + ε) · ‖p− c′‖z.
Unfortunately, we are not given a knowledge of OPT a priori. We therefore have to describe how to implement this lemma in a sublinear fashion. For this, we observe the following. First, the point q is an 2zα approximation with probability at least 1− 1/α. Second, the number of points that cost more than γ · 2zα · OPTn is at most 1/γ · n. Combining this observation with Lemma 2 ensures that we are not considering any points that are too far away. Since it is difficult to determine OPTn in a sublinear fashion, we will rely on additional pruning arguments for points that are close to q.
We proceed here in two steps. First, we consider the points that are very close to q. Second, we will show that the rings which contain too few points to be efficiently sampled have an overall negligible contribution to the cost.
Lemma 9. Suppose that q is an α-approximate solution. Let Anear ⊂ A be a set of points with cost at most (ε/(α5z))z · OPTn . Let Â = (1± ε)|Anear| Then for any candidate solution c we have∣∣∣∣∣∣Â · ‖q − c‖z − ∑ p∈Anear ‖p− c‖z ∣∣∣∣∣∣ ≤ ε/α ·  ∑ p∈Anear ‖p− c‖z + OPT
 . For rings with few points, we have the following.
Lemma 10. Suppose that q that is an α-approximate solution. Let Rcheap the union of rings with |Ri ∩ A| < ε · η · n and with radius at most 4(γ · α·OPTn )
1/z , where γ is given by Lemma 8. Then, for any candidate solution c∑
p∈Rcheap∩A ‖p− c‖z ≤ ε · 4z+1 · β · η · γ · α · ∑ p∈A ‖p− c‖z.
While we will defer the exact parameterization to later, observe that for β ∈ O(log 1/ε), α ∈ 2O(z) and η ∈ O(εz) this entire sum can be bounded by O(ε) ·OPT. The next lemma states that we can use these pruning results in a sublinear fashion.
Lemma 11. Let q be a point that is an α-approximation and let S be a uniform sample consisting of O ( α · η−1 · ε−3polylog(ε−1 · δ−1) ) points. Then with probability at least 1− δ for all rings Ri,
|Ri ∩A| − ε ·max (n · η, |Ri ∩A|) ≤ |Ri ∩ S| · n |S| ≤ |Ri ∩A|+ ε ·max (n · η, |Ri ∩A|) .
Furthermore, let d as in the algorithm, i.e., such that 23 · ε · η · |S| ≤ ∣∣S \ (B(q, d) ∩ S) ∣∣. Then
d < ( 3(εη)−1 · α·OPTn )1/z Since we do not know r, we also do not know which of the rings with i < 0, if any, satisfy 2i+1 ·d >
(ε/α5z) · (η/3α)1/z . However, the maximum number of these rings is at most log( 5αzε · ( 3α η )1/z ), which will turn out to be of the order O(log ε−1). For all of the rings that are not light, i.e. we cannot discard or snap to q, we now use Theorem 4. To ensure that we can call Theorem 4, we invoke Lemma 2 for every ring.
Probability Amplification While the aforementioned algorithm is guaranteed to produce a (1+ε)approximation with constant probability, amplifying this is non-trivial. Indeed, when running the algorithm multiple times, it is not clear how to distinguish a successful run from an unsuccessful one. The main issue in amplifying the probability lies in the initial solution q, as any invocation of Lemma 2 or Theorem 4 allows us to control the failure probability. The simplest way to achieve a success probability 1 − δ is to condition on ‖q − c‖z ≤ δ · OPTn . Unfortunately, this makes ε dependent on δ, which significantly increases the sampling complexity.
Instead, we use the following algorithm. We sample m ∈ O(log 1/δ) points q1, . . . qm uniformly at random. For each point, we additionally sample O(ε−2(log 1/ε + log 1/δ)) points Sqi . For qi, let di denote the minimum radius such that the points in n|Sqi |
· |B(qi, di) ∩ Sqi | > n2 . We output the point with minimal di.
The following lemma shows that this point is, with probability at least 1− δ, a 8z approximation. Lemma 12. Given query access to A, we can identify with probability 1 − δ a 8z-approximate solution using O(ε−2(log 1/ε+ log 1/δ) log 1/δ) samples.
To achieve an overall success probability of 1 − δ, we only need to sample from non-cheap rings. Thereafter, a high probability bound can be obtained by applying Theorem 4 applied to all Ri. The range space induced by rings centered around a single point q has constant VC dimension. Hence, Lemma 2 guarantees that a constant size sample will allow us to distinguish cheap rings from noncheap ones.

heading: 4 Improved Coreset Constructions
text: Our algorithm is as follows. First, we compute a point q that is a reasonably good approximation to the optimum4. In the following, let α = cost(q)OPT . Let ∆ = cost(q) n be the average cost of the input points when clustering them to q. We now partition the points into rings, defined as follow:
Ri = {p ∈ A | ( ε 2z )z · 1α · 2i ·∆ ≤ ‖p− q‖z ≤ ( ε2z )z · 1α · 2i+1 ·∆}. Let RM = ⋃log( 2zε )3zi=1 Ri. For each ring Ri with 1 ≤ i ≤ log ( 2z ε )3z , we sample a subset Si of s points uniformly at random, where each point is weighted by |Ri|s . We weigh q by the number of points in A \ RM . Finally, we set κ = ∑ i>log( 2zε ) 3z ∑ p∈Ri ‖p− q‖ z . We claim, for s ∈ Õ(ε−2) that for the weights defined
above, ⋃log( 2zε )3z i=1 Si ∪ {q} together with the constant κ is a coreset.
The analysis for every ring Ri ∈ RM is merely an application of Theorem 4. For the remaining points, we use the following lemma.
4We described an option in detail for the sublinear algorithm. In the interest of keeping the presentation succinct, we defer to that part of the paper and omit further discussion.
Lemma 13. Let q, RM and κ be defined as above. Then for any point c ∈ Rd, we have
|cost(A, c)− (cost(RM , c) + κ+ |A \RM | · ‖q − c‖2)| ≤ ε · cost(A, c).

heading: 5 Experimental Evaluation
text: While we can prove that the algorithm can compute a good solution in constant time for every constant ε and z, even for moderately small ε (e.g. ε = 1/2) the sampling complexity becomes quite large for even small values of z, as indeed our lower bound shows is necessarily the case. Our experiments therefore aim at evaluating the performance of the sublinear algorithm on realistic, not necessarily worst case data sets.
As baseline algorithm, we implemented a simple version of a batched gradient descent. Since all considered objectives are convex, we can expect such an algorithm to find a good solution in a reasonable time. The sublinear algorithm ran Algorithm 1 before calling the batched gradient descent. The code can be found at https://github.com/DaSau/power-mean We selected two data sets from the UCI repository Dua & Graff (2017), both of which are under the Creative Commons license. The first data set is the 3D Road Network data set from Kaul et al. (2013). It consists of elevation information with the attributes longitude, latitude and altitude. The total number of points is 434,874. For this data set, we considered all powers from z = 3 to 7. The second data set is the USCensus data set, consisting of the records from a 1990 census. The total size of the data set was 2,458,285 samples, each with 68 attributes. For this data set, we considered the powers z = 3, 4, 5.
The results essentially confirmed that the sublinear algorithm succeeded in finding a good candidate solution in a fraction of the time as batch gradient descent for essentially all considered problems.
A more extensive discussion can be found in the supplementary material.

heading: 6 Conclusion and Future Work
text: We gave sublinear algorithms and coresets for any power of means. Our bounds are nearly tight for the sublinear algorithms and we conjecture the coreset bound to be optimal, up to polylog factors.
The most immediate open question is whether our results generalize to coresets for k-clustering objectives. It seems likely that coresets of size O(k2/ε2) are achievable using our techniques. Improving on either this bound or theO(kε−2−z) from Cohen-Addad et al. (2021) is arguable the most important open problem in coresets.
In terms of sublinear algorithms, there is still a sub-optimal dependency on the ε by a factor ε−O(1). Obtaining tight bounds would be interesting. Finally, it is also an interesting open question whether sublinear algorithms for `p with p > 2 exist. It is known that for these spaces, no coreset that is independent of d can exist, even for the mean or the median. Is it nevertheless possible to obtain a sublinear algorithm that is independent of d?
Acknowledgments and Disclosure of Funding The work of David Saulpic is [partially] funded by the grant ANR-19-CE48-0016 from the French National Research Agency (ANR).

Review summaries:
Given n points in R^d, the goal is to compute a center x in R^d that minimizes the sum (dist(p,x))^z of distance to the power of z over every input point p. The main (sub-linear time) algorithm uniformly sample ~1/eps^z points and use smart processing to obtain a (1+epsilon) approximation to the desired center x. Nice application is differential privacy. Lower bounds are also provided. More generally, the authors suggest such coreset that approximates every query center.   

The paper studies sublinear algorithms for Power means i.e., given a set of $n$ points $a_1,\ldots, a_n$, compute a center $c$ that minimizes (upto $1+\varepsilon$) factor the objective $\sum_{i = 1}^n d(a_i, c)^z$ for any $z \ge 1$. The authors show that $\varepsilon^{-z-3}$ uniformly random samples is enough to achieve a $1+\varepsilon$ approximation. They also construct coresets of size $\varepsilon^{-2}$.

The paper gives sub-linear time algorithms for approximating ($1+\epsilon$ approximation) the centroid of a set of points in $d$-dimensional real space, under a distance function which is a power $z$ of the Euclidean distance. So, in particular, for $z=1$ the problem is estimating the Fermat-Weber point, and for $z=2$ the problem is estimating the mean. The algorithms use a random sample of the point set. Moreover, the paper also gives a probabilistic construction of a coreset of size independent of the dimension $d$ that estimates well the cost of any candidate center. The algorithms employ an interesting modification of an algorithm by Chen, as well as dimension reduction and a probabilistic analysis using chaining.