Simple and Optimal Sublinear Algorithms for Mean Estimation
Authors: Beatrice Bertolotti
Matteo Russo
Chris Schwiegelshohn
Sudarshan Shyam
Venue: NeurIPS 2025 poster ()
Keywords: Sublinear Algorithms
Mean Estimation
Learning
Primary area: theory
TL;DR: We present algorithms with optimal sample complexity and running time for estimating a high dimensional Euclidean mean.

Abstract: We study the sublinear multivariate mean estimation problem in $d$-dimensional Euclidean space. Specifically, we aim to find the mean $\mu$ of a ground point set $A$, which minimizes the sum of squared Euclidean distances of the points in $A$ to $\mu$. We first show that a multiplicative $(1+\varepsilon)$ approximation to $\mu$ can be found with probability $1-\delta$ using $O(\varepsilon^{-1}\log \delta^{-1})$ many independent uniform random samples, and provide a matching lower bound. Furthermore, we give two estimators with optimal sample complexity that can be computed in optimal running time for extracting a suitable approximate mean:
	1. The coordinate-wise median of $\log \delta^{-1}$ sample means of sample size $\varepsilon^{-1}$. As a corollary, we also show improved convergence rates for this estimator for estimating means of multivariate distributions.
	2.  The geometric median of  $\log \delta^{-1}$ sample means of sample size $\varepsilon^{-1}$. To compute a solution efficiently, we design a novel and simple gradient descent algorithm that is significantly faster for our specific setting than all other known algorithms for computing geometric medians.

In addition, we propose an order statistics approach that is empirically competitive with these algorithms, has an optimal sample complexity and matches the running time up to lower order terms.

We finally provide an extensive experimental evaluation among several estimators which concludes that the geometric-median-of-means-based approach is typically the most competitive in practice.

Simple and Optimal Sublinear Algorithms for Mean Estimation
Beatrice Bertolotti University of Pavia, Italy
Matteo Russo
Sapienza University of Rome, Italy mrusso@diag uniromal.... it
beatrice. beatrice. awr.:ooe.rrroott..eeeuuvvvaaaaar it
Chris Schwiegelshohn
Sudarshan Shyam
Aarhus University, Denmark
Aarhus University, Denmark shyan0cs... .aw.cc..auu... dk
schwiegelshhon@cs. .u.
Abstract
We study the sublinear multivariate mean estimation problem in d-dimensional Euclidean space. Specifically, we aim to find the mean js of ground point set A, which minimizes the sum squared Euclidean distances of the points in A We first show that a multipiccative (1 + E) approximation to can be found with probability 1 ô using O log d-1) many independent uniform random samples, and provide a matching lower bound. Furthermore, we give two estimators with optimal sample complexity that can be computed in optimal running time for extracting suitable approximate mean:
1. The coordinate- wise median log 5-1 sample means of sample size e- As corollary, we also show improved convergence rates for this estimator for estimating means of multivariate distributions.
The geometric median of log 5-1 sample means of sample size e-1_ To compute solution efficiently, algorithm that is significantly faster for our specific setting than all other known algorithms for computing geometric medians.
design novel and simple gradient descent
In addition, propose order statistics approach that is empirically competitive with these algorithms, has an optimal sample complexity and matches the running time up to lower order terms.
We finally provide extensive experimental evaluation among several estimators which concludes that the -aaeeerri-eeddia-----eeassssseee approach is typically the most competitive in practice.
Introduction
An extremely simple algorithmic paradigm is to sample subset of the input independently and uniformly at random and solve a problem on the sample. Such algorithms are called sublincar algorithms and form staple of large data analysis. The most desirable property of sublinear algorithm is that the query sample complexity is as small as possible. For most problems, rare for sublinear algorithm with only constant queries to produce multiplicative guarantees. Remarkably, finding the multivariate mean of a point set A in Euclidean space one of the cases where not only obtain multiplicative error with constant queries, but can do with arbitrary precision. Specifically, we are interested the objective
min
1p c
CERS
PEA
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

VE
क denotes the Euclidean norm. It is well known that the mean of the point set A minimizes this expression. We say that i (1 + appproxmmmae mean, the cost of c within (1 + s) factor of the cost of the optimal mean.
where || =
The exact sample complexity for obtaining approximate means is still open problem. Inaba et al. [1994] showed that E uniform and independent samples are necessary and sufficient for the empirical mean to be (1 E) -approximate mean on expectation. To extend this to high probability guarantee, that is, a success probability of 5, O samples are sufficient. Unfortunately, this bound is tight when using the empirical mean (see Appendix E). Cohen-Addad al. [2021] improved the dependence on the probability of failure by achieving a sample complexity of Ó(e~" log 8-1), at the cost of a significantly worse dependency on é. In the same paper, the authors also gave a lower bound of (? 2) for any sublinear algorithm succeeding with constant probability, showing that the result by Inaba et al. [1994] is optimal in the constant success probability regime. Subsequently, Musco et al. [2022] presented an algorithm with O ds="d d. log(ds- + log 8-1) log 8-1) queries, which combines the tight dependency on from Inaba al. [1994] with the improved dependency on 5 from Cohen-Addad et al. [2021], at the cost of a dependency on the dimension. Recently, Woodruff and Yasuda [2024] achieved almost best of all worlds trade-offs with sample complexity of O( (log= + log 5- log 1-).
1.1 Our Contributions
Our main contributions are three algorithms that estimate the mean of a point set and that have optimal sample complexity. In addition, the running times linear in the sample complexity for one and almost linear for the other two and thus optimal. In the first step of the algorithms,, see also Algorithm 1, we draw O(s log 8-1) points uniformly at random. We then split the sample into Olog 8-1) subsamples of size O(s-1), compute the empirical mean for each subsample, and run an aggregation procedure, i.e., a subroutine AGGREGATE(? ()) that combines the sample set produced to output an estimate of the mean.
Algorithm 1 MEANESTIMATE(C,,
for = .,....... log 8-1 do Sample S; points independently and uniformly at random with |S:| = as-1
Compute the sample mean ju; [3:] Epes,
Output AGGREGATE((2,.. (A10og6-1)
The sample complexity of this algorithm, m = b blog 8-1 as turns out to be optimal, as shown in Section F. Phrased as generalization question for bounding the excess risk given independent samples S from an underlying arbitrary but fixed distribution, we show an optimal convergence rate of ④
log
The remaining question is to find good solution (i.e., run AGGREGATE) quickly. Clearly, the sample complexity is also a lower bound for the running time of any algorithm. Achieving this rumning time and particular, breaking through the Ollog 8-1) barrier inherent in all previous algorithms, is substantial challenge, even for computing any constant approximation. It is fairly easy to show that most of the sample means are good estimators for the true mean. Selecting one of the good estimators is nevertheless difficult, as we have no way of estimating the cost and thus no easy way of distinguishing successful estimators from unsuccessful ones. Instead, we propose three algorithm to perform such a task.
Coordinate- Sose-Mediin-of-Meene and Geometric-MMedianof---ee wise-median -f-n--of-nn of-means and the geometric median-o of means are natural generalizations of the median- of-means estimator in single dimension. It turns out both have optimal sample complexity, up to constant factors. The former can be straightforwardly computed in linear time. Achieving a good running time for the latter is more challenging. Black box applications of gradient descent methods do not offer fast running times to compute suitable medians in our setting, but we show that a relatively simple gradient descent (see Algorithm 2) finds suitable solution in linear time.
Both the coordinate-

---

Using theoretical benchmarks, it is not possible to determine the more favorable estimator and it is not too difficult to construct instances where one outperforms the other. We therefore conduct an extensive experimental evaluation on real world data sets to measure the performance of these two estimators, along with previously proposed algorithms and our third contribution, the MinSumSelect algorithm.
MinSumSelect. Our other algorithm, given in Algorithm 3 (in Appendix B), uses a carefully chosen order statistic combined with clustering ideas. Specifically, the observation is that the sum of distances of the t closest points allows us to distinguish between a successful estimate and an unsuccessful estimate, for an appropriate choice of t. Here, the running time improves to O((==L + log? 8 8-1) dolog 0 1) for any constant y > 0. The algorithm also be parameterized to achieve a similar running time to the gradient descent algorithm, at the cost of increasing the sample size by poly (log log 8-1) factors.
1.2 Related Work
We covered the related work specific sublinear algorithms for estimating the mean earlier, but there are several other closely related topics. Given a distribution, a classical question is to design good estimators for some underlying statistic. For high-dimensional means, there has been substantial attention recently, see Lee and Valiant [2022], Lugosi and Mendelson [2019, 2024] and references therein. There are common notions with our work, as the distance to the true mean and tail bounds for the estimators are an important quantity in setting as well. Nevertheless, in this line of work, the running time, which is a very important focus for sublinear algorithms, is not a primary concern.
Another interesting line of work focuses on data dependent mean estimation, where the bounds derived are dependent on the quality of the data [Dang al., 2023]. There are limited regimes in which estimators with better than Lub-Gaussian estimators are possible. This is especially interesting given the recent efforts in the direction of beyond worst-case analysis of algorithms.
In a more general setting, one could ask for approximate power means minimizing the objective min Epe ||p cll for some positive integer z. For > 2, this problem was first posed by Cohen-Addad et al. [2021] who showed that Ō(e-=-3 log" 8-1) samples were sufficient for high success probability while 12(==*+1) samples were necessary for constant success probabil- ity. Woodruff and Yasuda [2024] improved this even further, obtaining sample complexity of O(E==11 -log-= + log 5 (oogg~2))_0 > 1-) Perhaps surprisingly. the case z = 1, also known as the geo- metric median requires ( - log d )ー2) many samples, see Chen and Derezinski [2021], Parulekar et al. [2021] for lower bounds and Cohen al. [2016] for an optimal algorithm. Finally, the special case + do corresponds to the minimum enclosing ball problem. The solution is highly susceptible to outliers, thus any sublinear algorithm must either be allowed to discard a fraction of the input, or rely on additive guarantees (see Clarkson al. [2012] and Ding [2020] for examples).
The sample complexity for the k-means problem has been extensively studided in the context of generalization bounds. Here, we are given an arbitrary but fixed distribution 'D typically supported on the unit Euclidean sphere, we aim to find a set of centers C' minimizing J, min Ip al TT dp] Following a long line of work, the best currently known learning rates with a sample size n for
this problem are the order O (v |blo + log 6-1 [S]
see Bartlett et al. [1998], Clemençcon [2011], Cohen-Addad et al. [2025], Fefferman al. [2016]. Klochkov et al. [2021], Linder [2000], Narayanan and Mitter [2010], Pollard and references therein. Interestingly, the best known lower bounds on the (V [ss]]
[S]
learning rates for arbitrary values of are at least [2023].
see Bartlett et al. [1998], Bucarelli et al.
Another related line of work is on coresets. For center-based objectives, a coreset approximates the cost of any given center with respect to the sum of distances raised to some power, where : = 2 is the power related to the 1-means objective. While there are some exceptions that obtain very small coreset sizes via other techniques, see Afshani and Schwiegelshohn [2024], Braverman et al. [2022], Huang et al. [2023], Maalouf et al. [2021, 2022], the sample complexity sensitivity sampling has by far been the most widely studied complexity measure in this line of work, see Bansal et al. [2024], Braverman et al. [2021], Cohen- Addad et al. [2021], Feldman and Langberg [2011], Langherg and Schulman [2010].
3

---

Further Comments and Limitations. We assume that the samples are drawn independently from data set, and that there are no corruptions. Efficiency and scalability are limitations, the present work only has optimal results in that regard. Regarding notions of privacy: Private mean estimation an important and well studied problem, but if we wish purely multiplicative approximation, as we aim to do here, no privacy guarantees are possible. Regarding fairness notions, sublinear algorithms cannot be individually fair in a meaningful way, as we never access the entire data set. Other clustering fairness notions, such as representational fairness, do not apply for single cluster.
2 Preliminaries
For a set of points C IRd whose cardinality is |A| = n, we use (( := 1o2 Epe A A to denote the mean, using f if A is clear from the context. Denote by OPT = min ΣpEA lpp - cll In addition, for set of points S sampled uniformly at random from A. we call (iS) the empirical estimator of p(A), or simply ğ if S is clear from context. A useful relationship between the goodness of fi and the size of |S] is established via the following identity.
Lemma 2.1 (Lemma of Inaba et al. [1994]). E [IIx(A) A(S))l?
OPT
[ST
At the heart of most algorithms for Euclidean means lies the following (arguably folklore) identity. Lemma 2.2 (High Dimensional Mean-Variance Variance Decomposition). For any set of points C R and any Rd, have Epen [p clI2 Epes [p nl|² + n -| cll².
<OPT, which we will often use throughout this paper. For space reasons, some proof's are omitted in the main body and deferred to the appendix.
Lemma 2.2 implies that fi is (1 + =) appoximation if and only if /| - اlر S V
For the proofs in the next sections, require basic probabilistic lemma. We say that an empirical =OPT We denote the set of good sample means by
mean is good, if [/ all $ where i v G, and define event
log o- }
=
G
10
The following lemma is straightforward application of the Chernoff bound. Lemma 2.3. With probability event & holds, for sufficiently large absolute constants and b.
Proof. We set the values of the constants to be 1440, b 50 and r r= 11 V eOFT (for reasons
which will be clear in the proof). First, we show that sample mean of as- samples is good with probability at least 0.9. This follows from Lemma 2.1 and Markov's inequality.
[lA 1ul2 < ,² (- (-cOPT)
eOPT
121n ]
2 0.9.
Applying the standard Chernoff bound (with b 50), get
P
Gl
S
blogo-1
S exp
(-
2
b bllo 8-1
S
81
where the is chosen so that
26 8i
<
-1
3 Coordinate-Wiss Median-of-Means
In this section, we show that the coordinate-wir median-of-means an optimal estimator and thus leads to a linear time algorithm for mean estimation. That is, at the end of Algorithm 1 we simply

---

use AGGREGATE(#1.... Iog procedure. We first do so the context of sublinear algorithms and then extend the result the case of distributional mean estimation in high dimensions.
8-1) COORD WOSRDDISEEE 1ng as an aggregation
3.1 CCordinate-Wiss Median-of -fenias in the Sublinear Model
We prove the following theorem:
Theorem 3.1. Algorithm run with NOWM == IOORDWISMEEIANNII... foottt as an aggregation routine oulputs a (1 -appproxmmaeee Euclidean mean with probability 6. The running time is O (== .olo )-)).
Proof. Consider VCWM, the coordinate wise median of all the sample means. For each coordinate let La and R3 be the sets of sample means such that f S VCWM,A and fu, 2 VCWM, respectively. We know that Lg, Rg have cardinality at least blog 5 Depending on whether VCWM,K > P4 or not, 2 show that at least one of the following statements is true:
1. For all jú in Tak, have that jit, soal |CCWM,K 1
2. For all ji in Rg we have that [ mmll |UCWM, 1A
We have Lg U R&l = blog which the total number of sample means. First, from the definitions of La and Rg. we observe that there points (such that ت사 PCWMLK) which belong to both. Along with the definition of the coordinate-we wise median, we have L. R 2 blog 8-1/2. Note that if there are no points such that Aa = MCWMAk then | = R& = b log -/2..
We provide case analysis as to why at least 1/5 fraction of the samples are good and satisfy |ic, - mal 2 |CWM, Mk]-
First, we note that when the good event & holds, at least 1/5 fraction of the means are good both the sets Lg and Rg. This follows from the fact that at least 7/10 fraction of the means good and the number of good means in L,/ Rg is at least 7/10- 1/2 1/5.
Case I: #l S VCWM,. In this case, we have for all means in RR. ( ukl 2 |UWMM, 1kkl-
Case II: μn VCWMA- In this case, we have for all means in Lis lpi. uol 2 |CWM, pkl-
From the above statements, we have for at least one La and RR. we have that at least 1/5 fraction of the means are good and satisfy JĀA, M Z |PCWM, 14..
From Lemma 1, we know that at least TO probability at least - 6. Hence, we infer that at least 7 satisfy Íūc, - pul 2 |PCWM.k - pu]. On average.
7 fraction of the means are good (i.e., llA #l S r) with of the sample means are good and
Σ
liL. 2
MCWM, M4
G
SEG
Summing over all coordinates & gives
PCWM,K & mk/2
VCWM All
Σ
kE[d]
TEG
kE[4
Interchanging the sums the left-hand side, get
5
5
||CWWM
1111 S
like 1ul1 < 5r?.
G
LEG KE[d]
VEG
The theorem follows from the definition
+ COPT and the running time follows by rec-
ognizing that we take O(log 6-1) sample means, each composed of O(e-1) many points in dimensions.
D
5

---

3.2 Learning the Mean of a High-Dimensional Distribution
In this section, we relate the results in the literature with respect mean estimation of high dimen- sional multivariate probability distributions and place our results in the context of these works.
We first describe the setting. The goal is to estimate the mean μ of probability distribution D in R A standard assumption is the existence of a covariance matrix Σ. We want the estimate to be close to the the mean μ and the Euclidean distance to the mean is used as the objective. Given N i.i.d. samples X1, x2,-.., n ~ D, the goal is to construct an estimator y = v(X1- X2...., ......,''' Xx) such that, for any confidence interval å E [0, 1J. with probability least 1 6, we have || -u--1h vil s(N, 6, 2.). The function eNN, 6, 2) is referred to as the rate of convergence of error henceforth.
The coordinate- wise median-of-meane while not optimal, is simple and natural extension of the algorithm for the 1-dimension. A simple analysis, found in several works, see Lugosi and Mendelson [2019] and Minsker [2015], applies union bound over all dimensions and concludes that if the mean is concentrated along every dimension, it will be concentrated in general. This analysis proves a
convergence rate of the order O TR(E) log(d8 1) N
where TR(2)) the trace of the covariance
matrix. Somewhat misleadingly, the dependency log d is often implied to be necessary.
The analysis of the coordinate- wise median-of-meant for sublinear algorithms be adapted to yield convergence rate for the mean estimation problem. Specifically, the coordinate-wise median-of-means estimator has the following. dimension- -fee convergence rate.
Theorem 3.2. Let X2.X..... NN be independent samples from distribution DD with variance μ. and covariance matrix E, then the coordinate-wit wise median of means estimator VCWM, with probability at least 1 - 8, satisfies
(R(E) logō-1
|PWM
rbll S 40
N
The proof is an adaptation of the proof for our setting and the details are deferred to the appendix.
4 Geometric Median-of-Means and Gradient Descent
In the previous section, showed that the coordinate wise median of the sample means is a (1 + e)- approximation of the true mean and can be computed linear time in the number of samples. However, one could generate instances where the coordinate-wise median is consistently a worse approximation to the true mean than the geometric median is, and vice-versa. This is shown in Appendix F. Therefore, we prove that the geometric median-of-means estimator also has an optimal sample complexity and give a gradient descent algorithm (Algorithm 2) for computing a sufficiently good estimate of it efficiently. It important to note that Algorithm 2 does not always converge the geometric median of the points. Take, for example. the triangle with coordinates (0, 1). (0, -1), (1. 0). The geometric median is the point (333 0). If we initialize the algorithm at the origin, the algorithm
never makes any progress and stays at the origin.
The gradient of the geometric median objective Eper |lp - ell at any point q with respect to some reference point set (q)) = Epep V,(4), where Vp(g) : - gradient from point p (excluding cocaoe points g to make well defined).
2-2 is the contribution to the
Algorithm 2 FASTGD(P)
Compute VCWM = COORDW
for 1,..,, do
Compute the gradient ((G-11) - pep Te--1-P
G1-1PP
Project all p3 E onto C--1 V(cy-1), denoting the projection by Pi,c Compute the I-dimensional median ej of the Pt.c,. 's along the line cj-1 ((c-1)
Output eT
The proof that the geometric median-of-ment estimator is optimal is algorithmic, i.e. it will follow as result of the analysis of Algorithm 2. We initialize Algorithm 2 with P {pi,.... ... Pb and

---

the coordinate-wise median VOWM of P as the starting estimate. The goal is to prove the following theorem:
Theorem 4.1. Algorithm I run with Algorithm 2 aggregation routine outputs (1 + E)- approximate Euclidean mean with probability 1 - if T 2 Y for some absolute constant The running time is O (E dloo 6 6)
The proof of Theorem 4.1 is articulated in two main steps. First, we show that, given a candidate estimate e, is very far from jr, the gradient of the geometric median of means evaluated at с, points in direction of Second, we must determine good updated solution along this gradient. Line-sweeps second-order methods all possible candidates, but computationally expensive. Instead, perform and update by selecting the median-of-meant along the gradient. Third, we show that once the algorithm determines a solution close μ. the gradient updates will remain close thereafter.
Fast convergence to close estimates from far ones. In order to estimate the mean quickly with gradient descent, we will make of the following geometric median gradient properties whenever our current estimate Cj is far from j..
Step 1:
Lemma 4.2. Conditioned on event E, we have:
(a) For any point q Id with l|g اll > 10r, then [V
3 bllog
10
(b) If, at iteration j. Algorithm 2 chooses a point ej such that ||c, (ell 2 10r, then I|c--1 S
TO ||0; - nl|--
?
Proof. We first show that every vector corresponding good mean is approximately in the same direction from Let á be good sample mean. Consider the triangle {q. p, f and let a, and 3 be the angle in q and ğ respectively. We show an upper-bound for cos (a;). Using the law of sines.
sin(as)
(3:)
1
114 Mlll I
| As
11 Mlll
We get sin(a.) <
14-p4 To which implies cos(a;) 2 V1- (a.)
0.99. To
bound (g)). we consider the contribution of good and bad means separately. We write
(4)) Lig Vi(g) Vo(2) + Vz(g), where Va, VB respectively denote good and bad points contributions the gradient evaluated at q (see Figure 1).
Vo(()) =
Σ
v.(a)
N
Σ
cos(a:)
N
7 10
0.99
bog 8og
all
jec good
d: tii good
where the first inequality holds because projecting the gradient onto any direction only decreases the norm. Conversely, the contribution given by bad sample means to the gradient norm is at most
Σ
v.(g)
S
blog 8-1
10
A. good
where the inequality follows from the upper-bound on the number of bad means when & holds and the triangle inequality. Part (a) follows by combining both bounds: We have, due to the triangle inequality, an overall lower bound on the norm of the gradient of
Σ v.(a)
Σ
Vi(g)
Σ
vi(g)
3 >
blo 8
10
8: As good
8: As nOt good
For part (b), similarly to the above, we focus on the triangle {g.9 V..g VG V }. Let 3 be the angle in vertex g (see Figure with q instead of e,) and let 0 be the angle vertex q Va VB. Using the bounds of Vaal and Rll and the law of sines, we get
sin(
(0)
V
V
Vell
7

---

which implies that sin(3) -
V -
Next, upper-bound ||c ju]| to prove part (b). Consider the projection of μ onto e Va(c)) (the green line in Figure I), call it y. Now, let us consider the projection of y onto сj (cy). call proj, proj,(y). (y). Since the good samples (including both good and bad) lies within the convex hull of the good projected means. This holds because, under any projection, the median must be between two good points (analogously to the one-dimensional median-of-mean estimator). This implies that, for Cj+1, we have
in the majority, the median cj+1 of the projected sample means
||c--1 - μ1| - ||os+1 - yll + 1y mlll S ||oz+1 - proij prrro (v) - yll +
Sr
-in(())- xlll
sin ||: ا|l
2r
- 10
|1e pl|
where the third inequality follows from the fact that |le - yll < |e, will (given that Gj p
and - y are respectively the hypotenuse and the leg of the right-angle triangle (cy, 14, y) and Ile; - All 10r.
We wish to briefly remark on the sample complexity of the geometric median-of-mean estimator and generalization bounds. Optimality follows from part (a) of Lemma 4.2, as if the gradient of the geometric median PGM is O, ||GM اlll s 10r must hold.
Corollary 4.3. For sufficiently large absolute constants a and b, the geometric median VGM == cll of b 1 many sample means consisting of as- many points is arg min E/ 1/As CER'
(1+ + E) apprrxxiiee mean with probability at least 1 8.
Remark 4.4. The geometric median cannot be computed exactly, can merely approximate it. The quality of the choice of approximation is important here. Using event & by itself. it is not clear that even a 2-approximate geometric median will recover a (1+ -aproxxiio mean.
Indeed, if assume that a single it, has distance V (1+6 + 8-8--00 (0-00-0)--aap appproximme median. In Appendix E. we will show that this assumption is warranted for worst case instances. The fastest algorithm, for finding a (1 + -approxximete geometric median by Cohen et al. [2016] runs in time O(d log )) In our setting, = b blog 8-1. The resulting runming time of Ollog* 8-1 d) is thus substantially slower that the other three algorithms presented in this paper.
the optimal mean, we must compute 1:OPT
Step 2: Once the estimate is close, it remains close.
Lemma 4.5. Conditioned on event E, if, at iteration j. Algorithm chooses a point oj such that I|c, xl| 10r, then [|cz+1-
- w|l < 11r.
Proof. Consider a point ej chosen by Algorithm 2 such that ||c, - ال| I 10r. Its gradient vector is V(c,) and Algorithm projects all other sample means onto C2 - Vc;,) and takes the median to be C++1- For every mean f., we denote its projection by proj(ū). We make two observations. First, the set of proj(ft of all good points lie bounded line segment of length at most 2r (see Figure 1). Second, due to the event ε, we have that the median chosen lies within the bounded line segment.
From the second observation, we have ||++1 proj(n)] < r, which is the radius of the good ball (see Figure 1). Applying triangle inequality, we get
||oz+1
اll ||c/11 pr) + [|pro((μ) All S + ||c, - =lll S 11r,
where the second inequality follows from (r) being a projection on line passing through ej.
Proof of Theorem 4.1. We condition on event E. In Algorithm 2, the initial value FCWM is taken to be the coordinate-wise median of all the 6 blog 8-1 sample means. From Theorem 3.1, we know that
||CCM - μ1ll S /5r. Using part (b) of Lemma 4.2, we have that as long as ||e; - za|| 2 10r, it holds that [|cz+1 - ا|| S 3. ||e, - pll-- As the quantity ||c, اull decreasing exponentially. it follows that for some absolute constant 7y there exists an iteration y, such that ||e, ااه < 10r. In the next
8

---

-Vis
Va(e)
Figure 1: Good empirical means are represented in green and bad ones in red. The ball of good means centered at f has radius r. Projection of all the good means lic on the bounded line segment of length most 2r.
iterations of the algorithm, from Lemma 4.5, we know that for the terminating value of the algorithm, or satisfies ||er - اl| S 11r, which yields the desired approximation ratio.
Concerning the the running time, we first need to compute b blog 8-1 sample means from as many uniform samples, and then, in each of the y iterations, there are b log 8 sample means to project onto the gradient vector, an operation that takes Odd log 8-1) time, with d being the dimension of the ambient space. To compute the median of all projections the initial coordinate-wise median VCWM can be done using linear-time rank select procedure (see [Cormen et al., 2009, Chapter 9.3]). For each iteration, as well as the initialization, this therefore takes (dd log 8-1) time. Thus, the overall running time is O (e-1 dloog 8 8-1).
every iteration, as well as finding
5 Experimental Evaluation
We conclude this paper with a short experimental evaluation. Code base and results can be found https: /ttpa://gitte gigg.uu....... com/ matteorusso sublinear carried out on the Google Colab default CPU. Our goal is to understand the practical viability of various estimators, such as the three optimal estimators included here, as well the previous results from Cohen- Addad et al. [2021]. Woodruff and Yasuda [2024] and the empirical mean of the samples. In particular, we wish to understand if the gradient updates from Algorithm improve over an initialization given by the coordinate- wise median. It is not difficult to generate instances where one of these algorithms outperforms the other, but ultimately the empirically best estimator can only be determined via experiments on real-world data.
sublier __ean_ estimation. The experiments were
Algorithms. The reference algorithms for fast sublinear mean estimation is the 1-means coreset- based CSS algorithm Coden-Addad et al. [2021], the oetive-regression ptreset-based WY algorithm Woodruff and Yasuda [2024], and our Algorithms and 3. In addition, we also used the standard empirical mean as a baseline. Each of these algorithms are given set of m points sampled uniformly random from the underlying point set and compared with respect to running time and quality of the solution.
We implement both MINSUMSELECTAdE FASTGD, As described in Section 4, we take initial guess the coordinate-wist median of the computed sample means. We also simply compare against the coordinate- wise median initialization step of FASTGD, referred as COORDWISEMEDIAN,
Data Sets, Setup and Results. We test our algorithms against the benchmarks mentioned above on the following datasets: MNIST and Fashion-MNIST, both of which are composed of 60,000 points, each with 784 features. We also consider Cover' Type, composed of 581,012 points, each with 54 features. We describe the gist of the results here, more details such as numerical figures are given in Appendix G. For each sample size m E {10, 15, 20, 25, 30, 100, 200, 500, 1000, 2000, 5000, 10000}, times and report averages and variances across runs.
we repeat the execution of every algorithm We plot the accuracy (i.e., ALG/OrT), and runing time function of the sample size for the various algorithms. Due to space constraints, we only report results on MNIST in this section (see Figure 2). The behaviour of accuracy and runtime as a function of sample size is similar across all datasets.

---

MOST (ub 1500C
198ST: 2t 108300 750
NIO adlect
(a) Accuracy vs. sample size (in log-scale)
(b) Runtime vs. sample size
Figure 2: MNIST Dataset: Accuracy and runtime sample size.
Our observations from the experimental evaluation real-world datasets provide few insights on the proposed methods. Among all algorithms, MINSUMSELECT and COORDWISEMEDI emerge as the fastest, consistently achieving superior runtime performance across all datasets. Their efficiency highlights their suitability for large-scale applications where computational speed is critical. In terms of accuracy, both UMSLEET, COORDWISEMEDI
and FASTGD perform comparably, with FASTGD showing a slight edge in most cases. Notably, FASTGD almost never underperforms relative to COORD WISEMEDIAN, despite both offering the same theoretical guarantees. This suggests that, in practice, applying a few gradient-base -based refinement steps toward the geometric median, starting from the coordinate- wise median, can yield measurable improvements.
Interestingly, despite having the weakest theoretical guarantees, the empirical mean consistently ranks among the best algorithms in terms of accuracy. Since it can also be computed extremely quickly via the closed form, it is also by far the fastest among the candidate algorithms. While there exist examples where the empirical mean does not provide good and robust estimation, these examples do not seem to appear on any of the data sets we tried and may not be common in practice, see the appendix for a more detailed discussion. Nevertheless, the empirical mean is not consistently the best estimator, which goes to the CSS-estimator, which also requires extensive filtering and preprocessing and is among the slowest available methods. The WY algorithm performs worse in terms of accuracy. While it is not known whether the stated sample complexity upper bounds for CSS and WY are tight, these observation may indicate that these algorithms are not tightly analyzed. Additionally, they are considerably slower than our algorithms, regardless of sample size, due to inherent overhead in pre- and postprocessing steps. However, optimization very fast for these algorithms because they simply output the mean.
In summary, the algorithms proposed in this paper are always competitive in terms of accuracy with the other methods, at least when given moderately large sample size, while being among the fastest methods available. While the refinement via FASTGD can offer improvements practice, the best candidate algorithms seem to be either COORDWISEMEDIAY or taking the empirical mean.
33
Crse
Wide iniidee
1.
200
200
2008
2006
2000
25
200
1000
2000
500
1960
(a) Accuracy vs. sample size (in log-scalc)
(b) Runtime vs. sample size
Figure 3: Fashion-MNIST Dataset: Accuracy and runtime against sample size.
10

---

Acknowledgements
We thank the anonymous reviewer from a previous version of this paper for pointing an improvement of analysis of the coordinate- wise median-of-men supported by the FAIR (Future Artificial Intelligence Research) project PE0000013, the NextGen- erationEU program within the PNRR-PE-AT scheme (M4C2, investment 1.3, line Artificial Intelligence), the PNRR MUR project IR0000013 3-So and by the MUR PRIN grant 2022FKNE5K (Learning in Markets and Society). Chris Schwiegelshohn
estimator. Matteo Russo was partially
supported by Google Research Award and by the Independent Research Fund Denmark (DFF) under a Sapere Aude Re- search Leader grant No 1051 10000106.. Sudarshan Shyam was partially supported by the Independent Research Fund Denmark (DFF) under grant 2032-00185B.
References
Afshani and C. Schwiegelshonn. Optimal coresets for low-dimensional geometric median. In Forty-first International Conference on Machine Learning, ICML 2024, Vienna, Austria, July 21-27, 2024. OpenReview.net, 2024. URL https: htttt::t/ppenrreee http:://///peeeeieeeeem http:://pppp..ee.oee..eeeeerraaa mttf:foruuummtdddaiiiikku
N. Bansal, V. Cohen-Addad, M. Prabhu, D. Saulpic, and C. Schwiegelshohn, Sensitivity sampling for k-means: Worst case and stability optimal coreset bounds. In 65th IEEE Annual Symposium
Foundations of Computer Science, FOCS 2024, Chicago, IL, USA, October 27-30, 2024, pages 1707-1723. IEEE, 2024. doi: N0.1109/FOCS612666 2024.00006 URL https::/ddo org/10. 1109/F0CS61266 2024.00106.
Bartlett, T. Linder, and G. Lugosi. The minimax distortion redundancy in empirical quantizer design. IEEE Transactions Information Theory, 44(5):1802-1811 1998. doi: 10. 10.1109/18.705660
V. Braverman, S. H. Jiang, R. Krauthgamer, and X. Wu. Coresets for clustering in excluded graphs and beyond. In D. Marx, editor, Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms, SODA 2021, Virtual Conference, January 10 13, 2021, pages 2679-2696. SIAM, 2021.
V. Braverman, V. Cohen-Addad, S. H. Jiang, R. Krauthgamer, C. Schwiegelshohn, M. B. Toftrup, and X. Wu. The power of uniform sampling for coresets. In 63rd IEEE Annual Symposium on Foundations of Computer Science, FOCS 2022, Denver, CO, USA, October 31 November 2022, pages 462-473. IEEE, 2022. doi: 10099FF005545722222000000 org/oogg//// 10. 101111109FF555 /ttp///OC555552225855
URL https:
M. S. Bucarelli, M. F. Larsen, C. Schwiegelshohn, and M. Toftrup. On generalization bounds for projective clustering. In A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, editors, Advances in Neural Information Processing Systems 36: Annual Conference on Neural Information Processing Systems 2023, NeurIPS 2023, New Orleans, LA, USA, December 10 16, 2023, 2023.
O. Catoni. Challenging the empirical mean and empirical variance: a deviation study. Ann. Inst. Henri Poincaré Probab. Stat., 48(4): 2012. ISSN 0246-0203.1778-0111 doi: 10.1214/ 11-AIHP454. URL https:://// / /tt////oi.ioo /ttp::////ooo.o./////////////ooa 1214/11-ATHP454, AIHP454.
X. Chen and M. Derezinski. Query complexity of least absolute deviation regression via robust. uniform convergence. In M. Belkin and S. Kpotufe, editors, Conference on Learning Theory, COLT 2021, 15-19 August 2021, Boulder, Colorado, USA, volume 134 of Proceedings of Machine Learning Research, pages 11 1449. PMLR, 2021.
K. L. Clarkson, E. Hazan, and D. Woodruff. Sublinear optimization for machine learning. ACM, 59(5):23:1 2012.
S. Clemengcon. On u-processes and clustering performance. J. Shawe-Taylor R. Zemel, P. Bartlett, F. Pereira, and K. Weinberger, editors, Advances Neural Information Processing Systems, volume 24, Curran Associates, Inc., 2011.
11

---

M. B. Cohen, Y. T. Lee. G. L. Miller, J. Pachocki, and A. Sidford. Geometric median in nearly linear time. In Wichs and Y. Mansour, editors, Proceedings of the 48th Annual ACM SIGACT Symposium Theory of Computing, STOC 2016, Cambridge, MA, USA, June 18-21, 2016, pages 9-21. ACM, 2016.
V. Cohen-Adder D. Saulpic, and C. Schwiegelshohn, Improved coresets and sublinear algorithms for power means in euclidean spaces. In M. Ranzato, A. Beygelzimer, Y. N. Dauphin, P. Liang, and J. W. Vaughan, editors, Advances in Neural Information Processing Systems 34: Annual Conference on Neural Information Processing Systems 2021, NeurlPS 2021, December 6-14, 2021, virtual, pages 21085-21098, 2021.
V. Cohen-Addad, Lattanzi, and Schwiegelshohn. Almost optimal PAC learning for k-means. In M. Koucký and N. Bansal, editors, Proceedings of the 57th Annual ACM Symposium on Theory of Computing, STOC 2025, Prague, Czechia, June 23-27, 2025, pages 2019-2030, ACM, 2025. doi: URL hatpps://dee ott.dii.org//ie /ttp///.111//37775580 37788180.
10.1145/37178233371188
H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein. Introduction to Algorithms, 3rd Edition. MIT Press, 2009.
Dang. J. C. H. Lee, M. Song, and P. Valiant. Optimality in mean estimation: Beyond worst-case, beyond sub-gaussian, and beyond ++a moments. In NeurlPS, 2023.
H. Ding. A sub-linear time framework for geometric optimization with outliers in high dimensions. In 28th Annual European Symposium on Algorithms, ESA 2020, September 7-9, 2020, Pisa, Italy (Virtual Conference), pages 38:1-38:21, 2020.
C. Fefferman, S. Mitter, and H. Narayanan. Testing the manifold hypothesis. Journal of the American Mathematical Society, 2944)99331, 2016.
D. Feldman and M. Langberg. A unified framework for approximating and clustering data. In Proceedings of the 43rd ACM Symposium on Theory of Computing, STOC 2011, San Jose, CA, USA, 6-8 June 2011, pages 569-578, 2011.
L. Huang, R. Huang, Z. Huang, and X. Wu. On coresets for clustering in small dimensional euclidean spaces. In A. Krause, E. Brunskill, K. Cho, B. Engelhardt, S. Sabato, and J. Scarlett, editors, International Conference USA, volume 202 of Proceedings of Machine Learning Research, pages 13891-13915. PMLR, 2023. URL https:tt//// htth::://ppoccedingg mlr. wwwps/eee000 / hhtp///haaaa......m/ html. wttp:/gs..........e/eeeeeee
Machine Learning. ICML 2023, 23-29 July 2023, Honolulu, Hawaii,
M. Inaba, N. Katoh, and H. Imai. Applications of weighted voronoi diagrams and randomization to muriance-based k- k-clustering (extended abstract). In Proceedings of the Tenth Anmual Symposium Computational Geometry, Stony Brook, New York, USA, June 6-8, 1994, pages 332-339, 1994.
Y. Klochkov, A. Kroshnin, and N. Zhivotovskiy. Robust k- k-means clustering for distributions with two moments. The Annals of Statistics, 49(4) 2230, 2021.
M. Langberg and L. Schulman. Universal e-approximators for integrals. In Proceedings of the Twenty-First Anmual ACM- Symposium on Discrete Algorithmms, SODA 2010, Austin, Texas, USA, January 17-19, 2010, pages 598-607, 2010. doi: 10.1137/1978166111337733.8 URL http://dx. doi. oo.. 10. 1137 01.3771..78111111 97816119737775.. 50.
J. C. H. Lee and P. Valiant. Optimal sub-gaussian mean estimation in very high dimensions. In M. Braverman, editor, 13th Innovations in Theoretical Computer Science Conference, ITCS 2022, January 31 February 3, 2022, Berkeley, CA. USA, volume 215 of LIPlcs, pages 98: 18:1- 1-9:211. Schloss Dagstuhl Leibniz Zentrum für Informatik, 2022. doi: IS 2022.99.. URL https:: / /ooi. org/10 4230/LIPIcs ITCS. CTCC..0222.9 2022.98.
Linder. On the training distortion vector quantizers. IEEE Transactions on Information Theory, 46(4):1617-1623, 2000. doi: 011100988/880000
G. Lugosi and S. Mendelson. Mean estimation and regression under heavy-tailed distributions: survey. Found. Comput. Math., 19(5):115--1190 2019.
12

---

G. Lugosi and S. Mendelson. Multivariate mean estimation with direction-dependent accuracy. J. Eur. Math. Soc. (JEMS). 26(6):2211 1-22477, 2024, ISSN 1435-9855,1335-98665 doi: 10.4171/jems/322. URL https:://dd...
A. Maalouf, I. Jubran, and D. Feldman. Introduction to coresets: Approximated mean. CoRR, 10/2111.03040 2021.
A. Maalouf, I. Jubran, and D. Feldman. Fast and accurate -eastemean-squares solvers for high dimensional data. IEEE Trans. Pattern Anal. Mach. Intell., 41112):9777-9994, 2022.
Minsker. Geometric median and robust estimation in banach spaces. Bermoulli, 21(4):2308-233 URL https: attps://www... jstor org/staable43559552.
2015. doi:
C. Musco, C. Musco, D. P. Woodruff, and T. Yasuda. Active linear regression for lp norms and beyond. In 2022 IEEE 63rd Annual Symposium on Foundations of Computer Science (FOCS), pages 744-753, 2022.
H. Narayanan and S. K. Mitter. Sample complexity of testing the manifold hypothesis. In J. D. Lafferty, C. K. Williams, Shawe-Taylor, R. S. Zemel, and A. Culotta, editors, Advances in Neural Information Processing Systems 23: 24th Annual Conference Neural Information Processing Systems 2010. Proceedings of meeting held 6-9 December 2010, Vancouver, British Columbia, Canada, pages 1786-179% Curran Associates, Inc., 2010.
A. Parulekar, A. Parulekar, and E. Price. L1 Regression with Lewis Weights Subsampling. In M. Wootters and I.. Sanità, editors, Approximation, Randomization, and Combinatorial Opti- mization. Algorithms and Techniques (APPROX/RANDOM 2021), volume 207 of Leibniz Interna- tional Proceedings in Informatics (LIPIcs), pages 49:1-49:22 Dagstuhl, Germany, 2021. Schloss Dagstuhl Leibniz-Zentrum für Informatik. ISBN 978-3-95977-207-55
D. Pollard. Central Limit Theorem for k-Means Clustering. The Annals of Probability, 10(4):919- 926.
D. P. Woodruff and T. Yasuda. Coresets for multiple lp regression. In Forty-first International Conference on Machine Learning, ICML 2024, Vienna, Austria, July 21-27, 2024, 2024.
13

---

A Preliminaries
Lemma A.1 (High Dimensional Mean- Variance Decomposition). For any set of points C I and any E R. have ΣPEA |p -01²
EpEA Ip ul + n || cll².
Proof. We have
Σ
c12
Σ
|lp + c||
PEA
PEA
Σ (Ilp 1112 + || -cl - 2(p 11) T (m c))
|pp 11112 + +n |/00 cll?
PEA
where the last equality follows from Epea (p- n1)
B Mean Estimation via Order Statistics
In this section, we leverage order statistics of the candidate means so as to allow for a quicker aggregation of suitable candidate mean.
Algorithm 3 MINSUMSELECT((P.)
Input: Set of points Pı....PIP) recursion depth if i or IP| 1 then W +- P
else
Split P arbitrarily into v Pliszze clusters {PI.... P Virl
W <
for each P, do W + W U (MINSUMSELLECT(P. 1)}
Output COMPUTEW
Algorithm 4 COMPUTEWINNER(?
for each Pj E P do
Let p/, be the 7 Pl-coosest point to Pj TO
Compute Dj
1 llp py ||
Output arg minp,eP Dj
The main result is the following:
Theorem B.1. Algorithum 1 run with Algorithm 3 aggregation routine outputs (1 + E) approximation with probability 1 5, using sumple size O (100'=-1 log 8-1), and running in time
o (110''s + +log² '8-1)
"d')dlog' dloog
for any non-negative integer
The key observation is that good estimate f4 minimizing the sum of distances of the th means closest to βy. In a nutshell, any mean minimizing such a must be close to sufficiently many successful estimates. Unfortunately, naive implementation of this idea takes time (log" 8-1 d), proven in the following lemma.
always be identified via COMPUTEWINNER
Lemma B.2. COMPUTEWINNERR takes time P| d).
14

---

Proof. We compute all pairwise distances between the points in P, which takes time O(1P d). To compute D,, we first have to find the To -llossett point to Pj. which takes time O(IP)) with sufficiently good rank select procedure, see [Cormen et al., 2009, Chapter 9.3]. Thereafter, summing up all the distances takes time P() per Pj E P.
Nevertheless, these scores yield an improved running time by arbitrarily partitioning the estimates into vPP groups, finding good estimate in every group via the truncated sum statistic each in time [P|.d, dd, for total rumning time of | d, and then selecting the best estimate via another truncated sum statistic on all estimates returned by the groups in time Pl. The final algorithm consists of applying this idea recursively.
To prove that MINSUMSELECT outputs a good solution, we require parameterized notion of a successful empirical mean. We say that a mean M4 7-good, if [|u; - ulll S 7V COPT
straightforward application of the Chernoff bound guarantees us that all but a very small fraction of the points are 7-good. Assuming this, the following lemma determines the quality of the computed solution.
Lemma B.3. Let P be set of means such that at least (1) of the means are y-good.
Then MINSUMSELECT((1) i) returns mean that is 5 r-good.
We prove this lemma by induction. We will use the following lemma in both the base case and the inductive step.
Lemma B.4. Given a set of means P, suppose that at least io Pl are y-good. Then, the estimate returned by COMPUTEWINNER is 5ryggood
Proof. First, let fi E P be any 7-good estimator. We know that the 1o Pl-closest estimators to ji
have distance at most y =OPT, which likewise implies that me Dj S To 1PIP1227 27y EOPT Now
suppose that juy is the estimator returned by the algorithm. Let G(āy) be the set of 7-good estimators
among the 7 Pl closest to Pa.. By assumption, we have [G(p)) I0 Pl. which implies I0
Σ (LEG(A)
1AA ixll -
7 Pl 2%
сOPT 7
y 1llly 5
EOP
10
m
which gives
2OPT
min (.E((A))
|1a - imll S 2
Therefore,
7
EOPT
EOPT
eOPT
[Aa - All S min Axll + |1
fu (G(())
[|A
m|l
S
271
2
n
7V
n
57
Proof of Lemma B.3. We proceed with the induction starting from i = 0.
Base Case, For the base case = 0, MINSUMSELECC((P.0) only calls COMPUTEWINNER, Thus, this case holds due Lemma B.4.
Inductive Step. Let PI.... P Viri be the clusters of P. computed when first calling
3 (+1 m PP| means that are not 7-good.
MINSUMSELECT((P.) By assumption, we have most This implies that the number of clusters with more than most
3 (10) PI means that are not y-good is at
(
P
VPP.
3
P
10
10
Denote this set by B(P) and let G(P) the remaining clusters. For each P, E B(P) G(P), let Øy returned by MINSUMSELECTTPP, 1). If P. G(P), may use the inductive hypothesis
15

---

which states that the mean ßi returned by MINSUMSELECT((J.- is 5'y-good. Since at least 10 /Pl of the thus computed means 5'oo may Lemma B.4 which shows that the final mean returned by COMPUTEWINNER (P, EB(P)UG())) (ū;}) is 7-good which concludes
the proof.
To conclude the proof, we require two arguments, First, we must show that, for appropriate choice of and b, at least
fraction of the means are 1-good with probability at least 1 1- 6, allowing us to use Lemma B.3. Second, we will argue the running time. The first a simple application of the Chernoff bound.
(10)
i+1
Lemma B.5. For a 2 2 25
(
i+1 and b = 3, at least (1) of the means are
5- 1-++1) -(+1)-oood with probability 1 1
{ if is not 5- good We have 16 - - 1. :OPT due to Lemma
Proof. Let X;
0
if is 5 (6+1) good
2.1. This implies due to Markov's inequality that P[X; = 1] S 25
which, by our choice of a, is
(10) 1t))t11 Thus, by the Chernoff bound
less than
blg 5-1
3
i+1
blog 5-1
25
X;
2
b blog 8-1
I TP
Σ
X 2 2.
log 8-1
10
a
25
b log 8-1)
I exp
3a
which is at most 8 by our choice of
The approximation guarantee now follows from Lemma B.3 and Lemma B.5. What remains to be shown is the running time. Since we split collection of means into vł clusters, the running time consists of the time required to recursively run the algorithm on the vł clusters each of size vł and which takes O(t d) time. Thus, starting with an instance of size Pol E 6logg511 ). we can solve for the recursion
consolidating via COMPUTEWINNER
VIPI-ITV/I + IPI |Pp
if 2 |
T(PP))
if S P
which yields the desired running time Oii lo 8-1 + (log 8 1 dd) = 0(11g ; 'd).
We next formalize this and complete the proof of Theorem B.I.
Proof of Theorem B.I. Throughout this proof, assume a 2 2 +1 (20)
i+1 and 4 3.
We first argue correctness, then rumning time. Let P denote the entire set of means passed to the MINSUM algorithm. By Lemma B.3 and Lemma B.5 and with our choices of a and
MINSUMSHLECT(P. returns a 1-good mean, that is (1 +) approxximetee mean with probability least 1 6.
What remains to be shown is the running time. Denote by |P/ = blog 5-1 the initial set of sample means. MINSUMSELECT(PP, computes v\Pl children, each of which recursively calls
MINSUMSELECT. Consolidating via COMPUTEWINNER takes time O O(V(P)- d) = O(PP) d) due
Lemma B.2. Thus the overall recursion takes time
(IPI-T(VDD |
[P) + +PP d if P0
T(P))
| < PP
which solves running time of O(i LP/ d 1+2- d) O((b log / d).
The time to compute the initial set of means is O (as- of and and some overestimation of the constants, leads to an overall running time of
-1 bd lollg811 8-1). which with our choice
O ((100'e-1 + (1og&-1)² .d log -dlog--1).
16

---

Remark B.6. For any constant choice of recursion depth, the sample complexity only increases by constants. We did not attempt optimize the constants, but the exponential dependency on is not avoidable. If we were to prioritize the running time over the sample complexity, can set a recursion depth i = log log log 8-1, which achieves running time and sample complexity of O( log 5 (log log 8- d).
Furthermore, if the recursion depth is shallow (i.e., for small values of i), the recursion can be improved via better choice of cluster size. Specifically, in the case i 1, that is with just one set of children, Theorem B.1 yields running time of ((e-1 + vlogo-1).ddae .logg If we instead choose (log 8 8-1)2/s nuny clusters, each consisting of Vlog6-1 many estimators, we obtain a running time of O (c-1 (/log-1) ed .loog (l Similar improvements for other values i are also possible, but these improvements become increasingly irrelevant compared the bounds given in Theorem B.1 as gets larger.
C Learning the Mean of a High-Dimensional Distribution
Theorem 3.2. Let , Xx be independent samples from distribution 'D with variance μ and covariance matrix E, then the coordinate-wist wise median least 1 8, satisfies
caddian-offan of means estimator VCwM, with probability at
TR Σ log 8-1
|PCWMM - 111l < 40V
N
Proof. The N samples are partitioned into 8log 8-1 subsamples of size 8log 8-1 each. The algorithm computes the empirical mean of each subsample and returns the coordinate-win median of the set of empirical means.
We prove the analog of Lemma 2.3 for this case, showing that a large fraction of sample means lie close to the true mean. For empirical mean of log samples, we have B [ STa(E)
N
We call an empirical mean good if || - all < where 13 TR and let G denotes the set of good means. Using TR(2.) = E||| - μl| and the Markov's 's inequality, we have Plβ G]
8 0.95. Applying the Chernoff bound, get that with probability least 1 8, we have G log 5-1
Consider VCWM. the coordinate-we wise median of all the sample means. For each coordinate k. let Lx and R2 be the sets of sample means such that üi, I VCWM, and fú, 2 MCWM, respectively. We know that Lk. R2 have cardinality at least blog 6-1 Depending on whether VOWM,K > f4 or not, at pkl-
least for one of Lk, RR. we have that | - uil N PCWM,
We know that at least 7 fraction of the means good (i.e., I| μll S r) with probability at of the sample means good and satisfy
least 1 5. Hence, infer that least [04. /4al 2 |.wMM, - /4k)- On average.
PCWM, - mkl
lin. - -m41 1x
G
IEG
Summing over all coordinates l: gives
Σ
PCWM,K 14k]
|CWM - /4]
5
5
kE[d]
iEG
kEId
Interchanging the sums on the left-hand side, we get
||CCMM
<
Σ
PPik - pul
Σ
114 - ull2 s 55.².
LEG RE(d)
IEG
The theorem follows from the definition of 13
TR(X)) 8-1881
17

---

D Generalization Bounds
We place our results the context of generalization bounds for clustering problems. In this setting. we are given an arbitrary but fixed distribution D> supported on the unit Euclidean ball Byg The cost of point with respect D is defined as costp(c) Speny lpp cll F(pldp. The risk is defined as R == argmin co)) Given independent samples S drawn from D, we wish to compute an
estimate è with cost R costp(c) such that the excess risk R R is minimized. The cost of distribution D is in a limiting sense the average cost of point set with all points living in By Since the average cost is at most the squared radius (i.e. 1), obtaining a (1 + E) approximate solution also yields a solution that has excess risk of at most :OPT S E, which we then rewrite in terms of sample size [S] as O log [ST
". This discussion is summarized the following corollary.
Corollary D.1. Given set of independent samples S drawn from some underlying arbitrary but fixed distribution supported on Bg, the geometric median-of menns estimator has an excess risk for the least squared distances objective of
log 8
R - R
<
[8]
with probability 8 for some absolute constant 7 >>
Note that the learning rate given by Corollary D.1 stands in contrast to learning rates that are achievable for other center-based problems such indeed the geometric median with objective
costp (c)
||p - cl P(pldp,
DEBd
or for the k-means problem with objective
costp(C)) = JE
min |pp eec
c|2 PPpjddp
for k-center set C. As mentioned the related work section, all of these objectives require learning rates of at least V ST ignoring problem specific parameters.
E Lower Bounds
We complement our algorithmic results with matching lower bounds.
E.1 A High Probability Lower Bound for
Theorem E.1. Any sublinear algorithm that outputs a (1 + -approxxima Euclidean mean with probability least ô must sample least ((s
1 log 8 8-1) many points.
The idea is to generate two instances that sublinear algorithm for approximating means has to distinguish between and then bound the number of samples required to distinguish between such distributions. The first instance is virtually identical to the one used by Cohen- Cohen-Addad Addad et al. [2021]. places n points at and en points at 1. Thus. the optimal mean is I+E A routine calculation via en Lemma 2.2 shows that the optimal cost I+2 The second instance places all points
Since is
not a sufficiently good approximate mean for the first instance, a sublinear algorithm can distinguish between the two instances, which yields lower bound on the sample complexity.
The two instances lie within the unit Euclidean ball. By normalizing the number of points, it also yields a distribution supported on unit Fuclidean ball, which implies that the generalization bounds given in Corollary D.1 are sharp.
Proof. We give two instances. The first instance places n points at 0 and en points 1. Thus, the optimal mean is placed - A routine calculation via Lemma 2.2 shows that the optimal cost is OPT The second instance places all points at 0.
18

---

First we argue that a sublinear algorithm can distinguish between these two instances. Any approxi- mate mean for the second instance must output O. If we output 0 for the first instance, we incur a cost of en = =m(1+=) = OPT + E OPT, hence any algorithm improving over a (1 + E) approximation for the former cannot output O. Thus, a necessary condition to distinguish between the two instances is for the algorithm to sample least one point at 1.
1+2
Suppose we sample m points. Our goal is to show that for m Q2 all of the sampled points are drawn from 0 is less than ó for the first instance. Let X; denote the indicator variable of the ith sampled point. We have
1 log 8 1). the probability that
m
8 2 PI i,
X;
0] = P 0
1 1+
1
exp
In
1 1+2 E exp (-m (n(1 + t))
> exp (-mE),
where the final inequality follows from the Mercator series In + =) = Σ (-1) Thus, for any m < g-1 log 5-g we do not achieve the desired failure probability. Conversely, m E S1s==loggg 8-1) for an algorithm to succeed.
E.2 The Empirical Mean is not Good Estimator
We briefly show that the arguably most natural algorithm that outputs the empirical mean of a subsampled point set has a substantial increase in sample complexity in the high success probability regime and therefore also a substantial increase in running time compared to the algorithms presented in this paper. The result is folklore, but included for completeness.
Theorem E.2. For all E, 8 E 1), there exists an instance such that M(c-18-1) independently sampled points are required for the empirical mean to be a (1 + -aappoxximatee Euclidean mean with probability at least 5 5.
The proof is based on reducing the problem of finding good approximate mean with high probability to giving an example distribution for which the Chebychev inequality is tight.
Proof. Consider a sample [S]. with |S] being specified later. We generate an instance as follows. We place _ل = fraction of the points each at - /S and V and the remaining 1 - [ fraction
at 0. Then the optimal solution places the mean at 0 and has average cost (I5)vE)2 T = 1. By
Lemma 2.2, this implies that the empirical mean म 1Sr Epes iss better than + - apappooxiee
mean if and only if [ < Ve. We set the failure probability, that is the probability P[lA| > vel 8.
Then P[A] Z VE] is at least the probability that we sample exactly one point from either - [S] VE or V and the remaining |S| points from O. Using Bernoulli's inequality and the density function of the binomial distribution, we therefore have
1
(1-
1
1
6 = P| N
SS|
S
1
Sl
S
T (1- [ for |S] then yields |S] E 2(e 18-1) 1, as desired.
Solving
Although the empirical mean is the fastest among all mean estimators (it can be computed in closed form), the above construction shows that it is not a robust estimator of the true mean-even in one dimension. Indeed, for heavy-tailed distributions, the best available tail bound is Chebyshev's inequality, which is exponentially weaker. In fact, Catoni [2012] (similarly to the construction above) shows that there exist distributions for which Chebyshev's inequality tight. In contrast, the empirical mean performs well when estimating the mean of subgaussian distributions, where the associated tail bounds are sharply exponential. In practice, data distributions are often not as adversarial as those constructed to demonstrate the failure of the empirical mean (e.g., Catoni [2012]). This explains why the empirical mean frequently observed to perform well in real-world scenarios, both in terms of computational efficiency and estimation accuracy.
19

---

F Comparison of Geometric Median and the Coordinate-wise median
In this section, we give two instances which demonstrate that neither of the coordinate-win median or the geometric median is strictly better than the other for the problem of mean-estimatio Proposition F.1. Let > 0 be fixed. There exist instances such that, with probability at O()),
1og6-1 2(t-e-1)² sample means of E- sampled points
(1) The coordinate wise median of K
coincides with the true mean, but the geometric median of the same K sample means does not;
(2) The geometric median of 8 sample means of sampled points is a better approxi- mation to the true mean than the coordinate-wise median is.
Proof. For part (I), we construct instance for which the coordinate- wise median is better estimator than the geometric median. Consider the uniform distribution on the d-dimensional sim- plex where d = O 28 8-22 log" 8-1). The value of d ensures that with probability at least - 6, each of the - log 8 vertices sampled are distinct. We show that in this case the coordinate- wise median is a better aggregation procedure than the geometric median. The coordinate me- dian VCWM the origin, while the true mean is l (1/d, 1/d, 1/d). 'he geometric dian can be shown to lie at the empirical mean of all the samples due to symmetry, which is log(8-))..... 0, ...... 0) where there are exactly log (o-1) non-zero coor- that || MCWM || S || puM|l. To conclude, note that
JGM (e/ log (og((--),,=
(6-))= // /
dinates. Comparing the distance, with probability at least 1 5, the coordinate-wir wise wise median is a better estimator than the geometric median.
For part (2), we wish to prove that there exists an instance in which the geometric median better than the coordinate- wise median. Consider the uniform distribution the d-dimensional simplex. where = cs_ We choose c to be a constant large enough so that in every subsample, the probability of a vertex sampled is $ (we only require it to be bounded away from 5) As in the last example, we have that the actual mean (1/d, 1/d,... .,. 1/d). We observe that the coordinate-wise median Vewm is the origin with high probability. While - O, the geometric median tends towards true mean μ.. D
G Further Experimental Evaluation
We present here further experimental evaluation and specifically we plot the accuracy vs. sample size and runtime vs. sample size, as well as the accuracy vs. dimension and runtime vs. dimension.
We first describe some summary statistics for the datasets we consider:
Dataset MNIST
Shape
Mean Std-Dev
Min
Max
(60,000; 784)
0.1306
0.3081
0.0000
1.0000
Dashion-MNIST
(60,000; 784)
0.2860
0.3530
0.0000
1.0000
CoverType
(581,012; 54)
0.4567
0.4981
0.0000
1.0000
Table 1: Summary statistics for datasets.
Accuracy and Runtime Sample Size. We plot how accuracy and runtime vary as function of the sample size. First, we plot this relationship in log-scale (Figures 4-5).
As they visually compress differences and obscure linear growth, especially when the slope is small, we replot accuracy vs. sample size in linear scale (Figure 6). The trends are indeed mildly increasing and consistent with near-linear growth, albeit with small slope. We believe the soft growth observed on the original dataset is largely due to implementation-lee NumPy's vectorized operations can lead to non-obvious runtime behavior, as such operations benefit from low-level optimizations (e.g.. memory locality, multi-threaded backends) and often incur fixed overheads.
factors. particular, our use of
We also add further experimentation for larger sample sizes (Figure 7). Specifically, we generated synthetic datasct from a 900-dimmnsio multivariate Gaussian distribution and extended the sample
20

---

MOST oub 1000C.7441
196S7: 14 180300 750
(a) Accuracy vs. sample size (in log-scale)
(b) Runtime vs. sample size
Figure 4: MNIST Dataset: Accuracy and runtime against sample size.
CONITTOON (54) 063112.51
Coartope- Outh (33101244)
Cb
13
(a) Accuracy vs. sample size (in log-scale)
(b) Runtime vs. sample size
Figure 5: CoverType Dataset: Accuracy and runtime against sample size,
size up to 5 10 The results are line with theory: They show that the error decreases sharply as the sample size grows and that the runtime grows linearly with sample size. Unfortunately, we were unable to go beyond this sample size due to memory limitations in our environment, where process exhausts the available RAM.
Accuracy and Runtime vs. Dimension. The experiments here illustrate how accuracy and runtime vary with dimensionality (Figure 8). Specifically, we generated synthetic dataset of 10, 000 samples drawn from multivariate Gaussian distributions with dimensions ranging from 10 to 1000. As predicted by theory, the errors of both the coordinate-win median and the fast gradient descent estimator remain low and stable across dimensions. Interestingly, the CSS algorithm shows rapid improvement the dimension increases, stabilizing around dimension 500. As expected, the empirical mean achieves the best accuracy since the data are drawn i.i.d. from a Gaussian distribution. In terms of runtime, observe roughly linear growth with dimension, in line with theory.
21

---

PABE (nd) 155000 7941
1.40
135
CSS
130
200,90
1.0
FRETEN-HUSTT (80000.25) 2011
Ocend bote Froden
:520122541
Dces Podee
CSS
Figure 6: MNIST, Fashion-MNIST and CoverType Datasets: Accuracy against sample size (in linear scale).
22

---

Random
144
Ocerd Packi
2.24
140
180
(a) Accuracy sample size (in linear scale)
Runtine vs Sample Size (Grouped for Samples 10.0001
11
10:0000
roson
Sample
(b) Runtime vs. sample size
Figure 7: Synthetic Dataset: Accuracy and runtime against sample size (in linear scale).
Randon Demoee Accenny F1 (EOROL Blmanstie
Mieura 14:200
6090
0.20
1210
00.00
0.21
0090
031
(a) Accuracy dimension (in lincar scale)
(b) Runtime vs. dimension (in linear scale)
Figure 8: Synthetic Dataset: Accuracy and runtime against dimension.
23

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: We provide both verbal and visual intuition as well as formal proofs of all claims stated in the paper. We also give an in-depth experimental evaluation of the proposed algorithms.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: Please refer the Limitations subsection before the preliminaries.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (c.g.. independence assumptions, noiseless settings, model sell-specfiiaatinn asymptotic approximations only holding locally). ). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-tex, system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
24

---

Answer: [Yes]
Justification:
Please refer to Sections 3, 4, and Appendices B, E.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proof's can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Question:
Answer: [Yes]
Justification: We provide full description of how to re re-produce the experiments and include the code that is set up to automatically run all of them.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, descríbing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e-g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing verifying the results.
5. Open access to data and code
25

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification:
Code included in supplementary and anonymous link available in main body.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.. cc/
httpg//ddde//ddeeubmiasssse sisiinncny))) for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https: / /nips. httpppulllluddeeeeeeeeppeppppeee
onPolicy) for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data,
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Question:
Answer: [Yes]
Justification: We provide several paragraphs describing how the experiments were set up and what the results mean. Please refer Section 5.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification; We report, for each sample size run the experiments, averages and variances of runtimes and accuracies arising from 50 executions of every algorithm. These are the most relevant statistical measures in our setting.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
26

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not show in tables figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification:
We include the hardware specifications for the experiments.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics https: / hww.rps.u...///e
Answer: [Yes]
Justification:
We study a basic theoretical question that has been studied before.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- cration due to laws or regulations their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: This work concerns efficient mean estimation in high dimensions and therefore does not have negative societal impacts.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
27

---

If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
No such risk posed.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g.. code, data, models), used in the paper, properly credited and the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [NA]
Justification: No existing assets were used.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
28

---

If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com/dataseta has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. If this information is available online, the authors encouraged to reach out the asset's s creators.
13. New assets
Question: Are assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: Our code is well-documented and should re-produce experiments out of the box (once the datasets are downloaded).
Guidelines:
The answer NA means that the paper does not release assets.
Researchers should communicate the details of the dataset/code/mooel part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: Not applicable.
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
29

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology. scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: Not applicable.
Guidelines:
The answer NA means that the method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy (https:/ /ttpar////.u...//eeeee ceaurips..ccCoofeeeeneeeee cc//omfferences//225/1/.. for what should or should be described.
30


Reviews:
[{"review_id": "7BBntJdreD", "reviewer": "Reviewer_otjC", "summary": "The paper suggests techniques to compute an approximation to the mean of a given n points in an Euclidean space. The approximation is (1+epsilon) multiplicative factor, using only O(log(1/delta)/epsilon) i.i.d samples. The state of the art is  (1/(epsilon*delta)) is tight, but assuming that the algorithm is just \"compute the mean of the sample\". Instead, the authors suggest aggregation techniques and median approaches. Experiments on some datasets is also provided, as well as code.", "questions": "The idea of using the median of small samples for computing the global mean seems very related to the famous paper:\nhttps://www.tau.ac.il/~nogaa/PDFS/amsz4.pdf\n\nPlease provide some comparison, at least for the algorithm itself.", "limitations": "No", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strong:\n-A fundamental problem in learning and statistics.\n- Open code.\n- Provable results.\n\nWeak:\n- The algorithm looks like a random version of a paper with a very similar goal that was appeared recently in NeuRIPS:\nhttps://arxiv.org/abs/1906.04705\nAlthough there are not so many papers regarding this problem, the paper is not cited and it seems that the authors are unaware of the result which essentially computes the mean by multiple aggregations. I am afraid that this fact alone suffices to reject this nice paper.\n\n- The experiments are very poor. Why MNIST which is considered small and a \"toy dataset\"? Why not datasets of billions of points if the running time is sub-linear? Even random numbers would do..", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for the useful comments, questions and suggestions. Moreover, we would like to highlight that, in accordance with the NeurIPS 2025 policy, we are not allowed to include any links in our rebuttal response, nor are we permitted to upload images—either in the rebuttal or to the already submitted repository. We are thus providing tables containing the same information as plots. We regret that you are unable to view the updated plots, and we kindly ask for your understanding and trust given these constraints.\n\n**Question.** The idea ... the algorithm itself.\n\n**Response.** The idea of median-of-means estimator is indeed old. Well-known references include (Nemirovski and Yudin, 1983; Jerrum et al., 1986 including the paper the reviewer has cited - Alon et al., 1996). Our work focuses on adapting the idea for higher dimensions. There are indeed many options for this including coordinate-wise median, geometric median among others. The proofs of convergence and running time we provide are novel to the best of our knowledge. Specifically, the best known algorithm for computing (approximate) geometric median (Cohen et al. 2016) leads to a much worse running time of the algorithm.\n\n**Comment.** The algorithm looks like ... reject this nice paper.\n\n**Response.** The paper in reference, though dealing with problems of a  similar flavor, has no direct connection to the problem we aim to solve. The paper aims at computing loss-less summaries for the mean. Sublinear algorithms necessarily require some form of approximation to the mean. The difference is also highlighted in the running time: their algorithm has complexity $O(nd + d^4 \\log n)$. It is unlikely it could be used at any stage of the mean aggregation. But even if it could be used, its running time would be significantly worse that what we are aiming for.  However, we will add a reference and explain the differences in the final version.\n\n**Comment.** The experiments are very poor ... random numbers would do.\n\n**Response.** Following the Reviewers otjC and 1w8V’s suggestions, we ran additional experiments with larger sample sizes. Specifically, we generated a synthetic dataset from a 200-dimensional multivariate Gaussian distribution and extended the sample size up to $5 \\times 10^5$. \n\nThe results are in line with theory: They show that the error decreases sharply as the sample size grows and that the runtime grows linearly with sample size (Tables 3 and 4).\n\nUnfortunately, we were unable to go beyond this sample size due to memory limitations in our environment, where the process exhausts the available RAM. We will add a note in the final version explaining this limitation. We hope that the extended results are sufficient to address the reviewers’ concerns.\nWe have also generated two plots illustrating how accuracy and runtime vary with dimensionality: We elaborated more on this aspect in the response to Reviewer sRHN. \nWe will add these new experimental results to the final version of the paper.\n\n\n**Table 3: Synthetic Dataset ($d = 200$) — Accuracies vs. Sample Size**\n\n|    $m$ |     EM |            EM $\\sigma$ |    FGD |          FGD $\\sigma$ |    CWM |          CWM $\\sigma$ |    MSS |          MSS $\\sigma$ |    CSS |           CSS $\\sigma$ |     WT |           WT $\\sigma$ |\n| -----: | -----: | ---------------------: | -----: | --------------------: | -----: | --------------------: | -----: | --------------------: | -----: | ---------------------: | -----: | --------------------: |\n|     10 | 1.0000 |  $1.17 \\times 10^{-1}$ | 1.1274 | $3.42 \\times 10^{-2}$ | 1.1276 | $3.49 \\times 10^{-2}$ | 1.2528 | $2.42 \\times 10^{-2}$ | 1.0000 | $1.11 \\times 10^{-16}$ | 1.2587 | $2.19 \\times 10^{-2}$ |\n|     15 | 1.0000 |  $8.88 \\times 10^{-2}$ | 1.0620 | $1.22 \\times 10^{-2}$ | 1.0671 | $1.31 \\times 10^{-2}$ | 1.1427 | $1.29 \\times 10^{-2}$ | 1.0000 |                 $0.00$ | 1.1442 | $1.34 \\times 10^{-2}$ |\n|     20 | 1.0000 |  $9.42 \\times 10^{-2}$ | 1.0549 | $1.05 \\times 10^{-2}$ | 1.0589 | $9.72 \\times 10^{-3}$ | 1.1195 | $1.29 \\times 10^{-2}$ | 1.0000 |                 $0.00$ | 1.1250 | $1.17 \\times 10^{-2}$ |\n|     25 | 1.0000 |  $6.08 \\times 10^{-2}$ | 1.0399 | $5.84 \\times 10^{-3}$ | 1.0393 | $7.14 \\times 10^{-3}$ | 1.0873 | $8.39 \\times 10^{-3}$ | 1.0019 |  $8.89 \\times 10^{-4}$ | 1.0915 | $8.92 \\times 10^{-3}$ |\n|     30 | 1.0000 |  $2.72 \\times 10^{-2}$ | 1.0316 | $5.53 \\times 10^{-3}$ | 1.0299 | $4.48 \\times 10^{-3}$ | 1.0672 | $6.72 \\times 10^{-3}$ | 1.0013 |  $9.67 \\times 10^{-5}$ | 1.0708 | $5.76 \\times 10^{-3}$ |\n|    100 | 1.0000 |  $6.47 \\times 10^{-2}$ | 1.0088 | $1.32 \\times 10^{-3}$ | 1.0091 | $1.18 \\times 10^{-3}$ | 1.0302 | $2.62 \\times 10^{-3}$ | 1.0005 |  $3.69 \\times 10^{-5}$ | 1.0313 | $2.94 \\times 10^{-3}$ |\n|    200 | 1.0000 |  $3.85 \\times 10^{-2}$ | 1.0058 | $7.01 \\times 10^{-4}$ | 1.0057 | $6.43 \\times 10^{-4}$ | 1.0194 | $2.18 \\times 10^{-3}$ | 1.0003 |  $2.54 \\times 10^{-5}$ | 1.0202 | $1.97 \\times 10^{-3}$ |\n|    500 | 1.0000 |  $1.66 \\times 10^{-2}$ | 1.0022 | $2.49 \\times 10^{-4}$ | 1.0022 | $2.41 \\times 10^{-4}$ | 1.0095 | $7.99 \\times 10^{-4}$ | 1.0001 |  $1.03 \\times 10^{-5}$ | 1.0102 | $8.94 \\times 10^{-4}$ |\n|   1000 | 1.0000 |                 $0.00$ | 1.0011 | $9.79 \\times 10^{-5}$ | 1.0011 | $9.97 \\times 10^{-5}$ | 1.0048 | $4.04 \\times 10^{-4}$ | 1.0001 |  $6.40 \\times 10^{-6}$ | 1.0051 | $4.65 \\times 10^{-4}$ |\n|   2000 | 1.0000 |  $5.44 \\times 10^{-2}$ | 1.0006 | $6.73 \\times 10^{-5}$ | 1.0006 | $7.67 \\times 10^{-5}$ | 1.0028 | $2.76 \\times 10^{-4}$ | 1.0001 |  $3.81 \\times 10^{-6}$ | 1.0031 | $2.75 \\times 10^{-4}$ |\n|   5000 | 1.0000 |  $9.93 \\times 10^{-3}$ | 1.0002 | $2.73 \\times 10^{-5}$ | 1.0002 | $2.63 \\times 10^{-5}$ | 1.0013 | $1.07 \\times 10^{-4}$ | 1.0001 |  $2.04 \\times 10^{-6}$ | 1.0014 | $1.47 \\times 10^{-4}$ |\n|  10000 | 1.0000 | $1.17 \\times 10^{-10}$ | 1.0001 | $1.21 \\times 10^{-5}$ | 1.0001 | $1.04 \\times 10^{-5}$ | 1.0007 | $6.10 \\times 10^{-5}$ | 1.0001 |  $1.59 \\times 10^{-6}$ | 1.0008 | $7.20 \\times 10^{-5}$ |\n| 100000 | 1.0000 | $3.14 \\times 10^{-18}$ | 1.0000 | $1.43 \\times 10^{-6}$ | 1.0000 | $1.27 \\times 10^{-6}$ | 1.0001 | $6.64 \\times 10^{-6}$ | 1.0000 |  $6.10 \\times 10^{-7}$ | 1.0001 | $1.00 \\times 10^{-5}$ |\n| 500000 | 1.0000 | $3.14 \\times 10^{-18}$ | 1.0000 | $2.92 \\times 10^{-7}$ | 1.0000 | $2.52 \\times 10^{-7}$ | 1.0000 | $2.00 \\times 10^{-6}$ | 1.0000 |  $3.98 \\times 10^{-7}$ | 1.0000 | $2.00 \\times 10^{-6}$ |\n\n**Table 4: Synthetic Dataset ($d = 200$) — Runtimes vs. Sample Size**\n\n|    $m$ |      FGD |          FGD $\\sigma$ |      CWM |          CWM $\\sigma$ |      MSS |          MSS $\\sigma$ |      CSS |          CSS $\\sigma$ |       WT |           WT $\\sigma$ |\n| -----: | -------: | --------------------: | -------: | --------------------: | -------: | --------------------: | -------: | --------------------: | -------: | --------------------: |\n|     10 | 0.023158 | $2.01 \\times 10^{-3}$ | 0.000424 | $1.98 \\times 10^{-4}$ | 0.000495 | $6.10 \\times 10^{-5}$ | 0.000229 | $1.10 \\times 10^{-5}$ | 0.000466 | $2.60 \\times 10^{-5}$ |\n|     15 | 0.023287 | $2.97 \\times 10^{-3}$ | 0.000434 | $3.57 \\times 10^{-4}$ | 0.000663 | $3.94 \\times 10^{-4}$ | 0.000339 | $2.44 \\times 10^{-4}$ | 0.000701 | $5.38 \\times 10^{-4}$ |\n|     20 | 0.016261 | $4.85 \\times 10^{-3}$ | 0.000245 | $5.00 \\times 10^{-5}$ | 0.000349 | $2.80 \\times 10^{-5}$ | 0.000164 | $1.80 \\times 10^{-5}$ | 0.000338 | $3.30 \\times 10^{-5}$ |\n|     25 | 0.014721 | $2.97 \\times 10^{-3}$ | 0.000273 | $6.80 \\times 10^{-5}$ | 0.000385 | $7.80 \\times 10^{-5}$ | 0.000167 | $8.00 \\times 10^{-6}$ | 0.000392 | $9.60 \\times 10^{-5}$ |\n|     30 | 0.015118 | $3.30 \\times 10^{-3}$ | 0.000248 | $2.50 \\times 10^{-5}$ | 0.000353 | $2.10 \\times 10^{-5}$ | 0.000173 | $1.20 \\times 10^{-5}$ | 0.000364 | $3.60 \\times 10^{-5}$ |\n|    100 | 0.017072 | $3.41 \\times 10^{-3}$ | 0.000379 | $6.40 \\times 10^{-5}$ | 0.000439 | $1.90 \\times 10^{-5}$ | 0.000282 | $4.30 \\times 10^{-5}$ | 0.000625 | $2.20 \\times 10^{-5}$ |\n|    200 | 0.019220 | $3.32 \\times 10^{-3}$ | 0.000522 | $6.40 \\times 10^{-5}$ | 0.000655 | $4.50 \\times 10^{-5}$ | 0.000479 | $5.80 \\times 10^{-5}$ | 0.001034 | $9.30 \\times 10^{-5}$ |\n|    500 | 0.021631 | $3.98 \\times 10^{-3}$ | 0.000835 | $5.20 \\times 10^{-5}$ | 0.001078 | $1.04 \\times 10^{-4}$ | 0.000850 | $4.90 \\times 10^{-5}$ | 0.002033 | $1.54 \\times 10^{-4}$ |\n|   1000 | 0.022246 | $2.99 \\times 10^{-3}$ | 0.001404 | $7.10 \\times 10^{-5}$ | 0.001841 | $4.40 \\times 10^{-4}$ | 0.001629 | $2.26 \\times 10^{-4}$ | 0.003663 | $3.48 \\times 10^{-4}$ |\n|   2000 | 0.025202 | $3.42 \\times 10^{-3}$ | 0.002930 | $5.24 \\times 10^{-4}$ | 0.003299 | $3.50 \\times 10^{-4}$ | 0.002832 | $2.04 \\times 10^{-4}$ | 0.007801 | $1.19 \\times 10^{-3}$ |\n|   5000 | 0.047089 | $6.56 \\times 10^{-3}$ | 0.010890 | $9.59 \\times 10^{-4}$ | 0.008105 | $5.00 \\times 10^{-4}$ | 0.010210 | $1.04 \\times 10^{-3}$ | 0.018438 | $1.52 \\times 10^{-3}$ |\n|  10000 | 0.042100 | $3.96 \\times 10^{-3}$ | 0.016434 | $9.72 \\times 10^{-4}$ | 0.016744 | $1.03 \\times 10^{-3}$ | 0.022541 | $7.60 \\times 10^{-4}$ | 0.037705 | $2.83 \\times 10^{-3}$ |\n| 100000 | 0.256733 | $2.73 \\times 10^{-2}$ | 0.227310 | $2.16 \\times 10^{-2}$ | 0.223843 | $1.91 \\times 10^{-2}$ | 0.339804 | $2.20 \\times 10^{-2}$ | 0.483566 | $7.13 \\times 10^{-2}$ |\n| 500000 | 1.321969 | $1.16 \\times 10^{-1}$ | 1.278243 | $1.11 \\times 10^{-1}$ | 1.277424 | $1.12 \\times 10^{-1}$ | 1.695665 | $1.05 \\times 10^{-1}$ | 2.782799 | $3.10 \\times 10^{-1}$ |"}, {"review_id": "96Yzr3qO68", "reviewer": "Reviewer_1w8V", "summary": "The authors present several novel methods to estimate a multivariate mean of a ground data, including theoretical proofs for optimality of sampling complexity. The methods are built as variants of of an algorithm that computes several sample means and is parametrized by an aggregation function for said means. The authors prove theoretical guarantees of approximation optimality and asymptotic relative runtime optimality of those methods, and evaluate them against MNIST/F-MNIST datasets, comparing to other baseline methods.", "questions": "I believe the clarity of Lemma 2.3 should be adjusted, as per discussion before.\n\nWhile this is a more theoretical paper (or even more importantly because of that), I believe an experimentation on synthetic data could be crucial here. In the simplest case, the same experiments from MNIST could be run on a sample from e.g. a Gaussian distribution, perhaps also testing limits with much larger sample sizes.\n\nOne question I have about Algorithm 2: is there any statement about algorithm's asymptotic convergence to true geometric median (in relation to T)? Unless I missed it, perhaps it is good to clarify whether it's an open question, or perhaps it is easy to construct a set of points that has a unique geometric median, but the method does not go infinitesimally close to it.\n\nI believe addressing these issues would bring the paper to a higher quality.", "limitations": "Yes.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "The paper is written in a very clear English and good formulations. The authors strive for theoretical soundness, and for the major part provide thorough theoretical justification for their claims.\n\nWhile the methods, as the paper's title suggest, are simple, this adds more clarity to their analysis. Moreover, the fast gradient method described for one of the algorithms, that essentially replaces a gradient step with a calculation of one-dimensional median along the gradient direction, is a rather neat idea. \n\nPerhaps because of limited space, some proofs in the main paper feel too condensed and hop over several steps, reducing the clarity. One particular place is the beginning of Lemma 2.3's proof that starts with the phrase \"Recall that the values of the constants were set as a = 1440, b = 50\". In all fairness, the values a and b were barely introduced; they first mentioned inline in Algorithm 1, and it's trivial how they appear from the big-O of sample sizes. However neither these particular values 1440 and 50 are mentioned before, nor is it clear how exactly they are used. \n\nAnother example of a jump in the proof, to me, is in line 164 and the next inequality. Perhaps a better way to write the inequality, while keeping it short, would be to include an part in-between, like [average across G] $\\ge$ [average across all means] $\\ge$ RHS. But this is a minor detail, and I believe the proof is correct (small correction, line 162: should the referenced lemma be 2.3?).\n\nA bit larger concerns are related to experemental evaluation and method significance. One strength is the amount of different methods/baselines. However:\n* A sample size $m$ is introduced. But how does it relate to $a$, $b$, $\\delta$ or $\\varepsilon$? I believe I can see that in the experiment code, but the description is missing from the paper.\n* Three fixed datasets are tested, where it is also shown that the empirical mean outperforms most methods. I strongly believe that testing on synthetic data sampled from more clean distributions would have been necessary for any conclusions, e.g. a large ground Gaussian/uniform sample. Also, not sure if accuracy should be plotted in log scale, or whether the accuracy valus in large sample sizes are not meaningfully distinct.\n* Runtime barely grows with the sample size. Unless I'm misunderstanding what sample size is, should not it be essentially linear for introduced methods? Could the issue lie in vectorized numpy implementations with bottlenecks around the calls? Wouldn't it make sense to try an aforementioned synthetic dataset with a large sample size? And if it is not supposed to be linear, then again a clarification about sample size is important.\n\nLast small question: line 179: what kind of non-optimality is spoken about here?", "was_revised": true, "final_justification": "I believe the additional experiments, clarification to existing figures and the addition of some missing definitions improve the quality of the paper. \n\nBased on authors' reply and how I imagine the new version, I readjust my grading to \"Borderline Accept\". For a full \"Accept\" rating, one would have to verify the final version of the paper (which is currently not available due to review process restrictions) with updated plots to also reevaluate the significance/quality of the experimental results. The rest of the paper is solid.", "rebuttal": "We would like to thank the reviewer for the valuable comments and the raised questions. Due to NeurIPS 2025 restrictions, we’re unable to share links or images in the rebuttal or repository. We’ve included equivalent information in tabular form and kindly appreciate your understanding of these constraints.\n\n**Comment.** A sample size $m$ is introduced ... missing from the paper.\n\n**Response.** The sample size $m$ is set to $ab \\varepsilon^{-1} \\log \\delta^{-1}$. To clarify, this is the result of having $b\\log \\delta^{-1}$ subsamples of size $a \\varepsilon^{-1}$ each. We thank the reviewer for the comment and will add it to the paper explicitly. \n\n**Comment and Question.** Three fixed datasets ... Gaussian/uniform sample; While this is ... much larger sample sizes.\n\n**Response.** Please refer to the last response to reviewer otjC and tables therein.\n\n**Comment.** Also, not sure if accuracy should be plotted in log scale, or whether the accuracy values in large sample sizes are not meaningfully distinct.\n\n**Response.** Following the reviewer's suggestion, we plotted the accuracy in linear scale for both MNIST and the randomly generated dataset. Indeed, the accuracy values for large sample sizes essentially coincide (Tables 1 and 7).\n    \n**Comment.** Runtime barely grows ... sample size is important.\n\n**Response.** We thank the reviewer for pointing this out. The runtime plots in the paper are shown in log-scale, which can visually compress differences and obscure linear growth, especially when the slope is small. When replotted in linear scale, the trends are indeed mildly increasing and consistent with near-linear growth, albeit with a small slope. We also repeated the experiments on the randomly generated synthetic dataset referenced above, and in that case the runtime growth is visibly linear. While we cannot upload plots, we have provided tables summarizing the above information (Tables 1 and 2 below, and 3 and 4 in reviewer otjC's last response). We believe the soft growth observed on the original dataset is largely due to implementation-level factors. In particular, our use of numpy's vectorized operations can lead to non-obvious runtime behavior, as such operations benefit from low-level optimizations (e.g., memory locality, multi-threaded backends) and often incur fixed overheads. As a result, runtime may appear sublinear over moderate input sizes, especially when memory access. Lastly, we agree that applying the methods to significantly larger datasets (as we did in the synthetic experiments) helps confirm the expected linear scaling more clearly, and we are happy to include these linear-scale plots and additional clarifications in the final version of the paper.\n\n**Table 1: MNIST — Accuracies vs. Sample Size**\n\n| $m$   | EM       | EM $\\sigma$     | FGD      | FGD $\\sigma$    | CWM      | CWM $\\sigma$    | MSS      | MSS $\\sigma$    | CSS      | CSS $\\sigma$    | WT       | WT $\\sigma$     |\n| ----- | -------- | -------- | -------- | -------- | -------- | -------- | -------- | -------- | -------- | -------- | -------- | -------- |\n| 10    | 1.098528 | 0.028635 | 1.155126 | 0.034458 | 1.175392 | 0.037500 | 1.392344 | 0.063439 | 1.001193 | 0.000251 | 1.050220 | 0.109802 |\n| 15    | 1.089786 | 0.022853 | 1.091481 | 0.021414 | 1.089386 | 0.019651 | 1.175579 | 0.041341 | 1.001193 | 0.000290 | 1.183068 | 0.040979 |\n| 20    | 1.050944 | 0.017004 | 1.071097 | 0.016169 | 1.079405 | 0.015386 | 1.150386 | 0.033795 | 1.001149 | 0.000299 | 1.166793 | 0.038011 |\n| 25    | 1.042582 | 0.012022 | 1.058092 | 0.013737 | 1.056672 | 0.011956 | 1.108955 | 0.023705 | 1.001170 | 0.000305 | 1.120348 | 0.024574 |\n| 30    | 1.035054 | 0.009981 | 1.042930 | 0.009194 | 1.045296 | 0.009135 | 1.091521 | 0.016239 | 1.001194 | 0.000290 | 1.097396 | 0.022739 |\n| 100   | 1.009560 | 0.002123 | 1.011862 | 0.002859 | 1.012232 | 0.002303 | 1.041629 | 0.011196 | 1.001229 | 0.000242 | 1.042088 | 0.011341 |\n| 200   | 1.004750 | 0.000965 | 1.007189 | 0.001480 | 1.007697 | 0.001541 | 1.022038 | 0.005945 | 1.001251 | 0.000228 | 1.025499 | 0.005615 |\n| 500   | 1.001954 | 0.000498 | 1.002961 | 0.000648 | 1.002704 | 0.000572 | 1.010661 | 0.002279 | 1.001248 | 0.000221 | 1.012654 | 0.002951 |\n| 1000  | 1.000921 | 0.000239 | 1.001301 | 0.000301 | 1.001286 | 0.000261 | 1.005201 | 0.001392 | 1.001310 | 0.000216 | 1.006786 | 0.001536 |\n| 2000  | 1.000493 | 0.000131 | 1.000700 | 0.000114 | 1.000717 | 0.000105 | 1.002898 | 0.000617 | 1.001290 | 0.000185 | 1.003791 | 0.000942 |\n| 5000  | 1.000180 | 0.000054 | 1.000294 | 0.000072 | 1.000270 | 0.000071 | 1.001383 | 0.000327 | 1.001360 | 0.000165 | 1.001820 | 0.000490 |\n| 10000 | 1.000089 | 0.000027 | 1.000142 | 0.000027 | 1.000151 | 0.000031 | 1.000692 | 0.000117 | 1.001340 | 0.000169 | 1.000996 | 0.000299 |\n\n**Table 2: MNIST — Runtimes vs. Sample Size**\n\n| $m$   | FGD      | FGD $\\sigma$    | CWM      | CWM $\\sigma$    | MSS      | MSS $\\sigma$    | CSS      | CSS $\\sigma$    | WT       | WT $\\sigma$     |\n| ----- | -------- | -------- | -------- | -------- | -------- | -------- | -------- | -------- | -------- | -------- |\n| 10    | 0.331120 | 0.025323 | 0.303907 | 0.021514 | 0.305031 | 0.022829 | 0.752144 | 0.047994 | 0.498704 | 0.101033 |\n| 15    | 0.328824 | 0.029853 | 0.305248 | 0.023860 | 0.304551 | 0.019328 | 0.750669 | 0.052915 | 0.486801 | 0.095939 |\n| 20    | 0.321185 | 0.022370 | 0.313300 | 0.026494 | 0.304787 | 0.021680 | 0.753925 | 0.049201 | 0.469229 | 0.060214 |\n| 25    | 0.320447 | 0.024325 | 0.305469 | 0.018924 | 0.312867 | 0.026952 | 0.758950 | 0.052323 | 0.454810 | 0.045884 |\n| 30    | 0.324699 | 0.023393 | 0.306246 | 0.021164 | 0.307380 | 0.021195 | 0.754369 | 0.050649 | 0.484448 | 0.061920 |\n| 100   | 0.332200 | 0.029719 | 0.305576 | 0.022594 | 0.304266 | 0.021536 | 0.753291 | 0.047923 | 0.471533 | 0.064969 |\n| 200   | 0.339868 | 0.032710 | 0.307593 | 0.021598 | 0.306124 | 0.020822 | 0.753950 | 0.052971 | 0.472828 | 0.057820 |\n| 500   | 0.341225 | 0.032963 | 0.311128 | 0.020642 | 0.309915 | 0.021704 | 0.757828 | 0.039416 | 0.476949 | 0.061728 |\n| 1000  | 0.340164 | 0.029454 | 0.313098 | 0.025072 | 0.308668 | 0.021177 | 0.757372 | 0.050639 | 0.478237 | 0.062958 |\n| 2000  | 0.347602 | 0.029652 | 0.320347 | 0.026740 | 0.317822 | 0.022311 | 0.760508 | 0.050913 | 0.498500 | 0.058959 |\n| 5000  | 0.360183 | 0.035123 | 0.320959 | 0.022485 | 0.322097 | 0.022088 | 0.754888 | 0.051201 | 0.496606 | 0.064727 |\n| 10000 | 0.369807 | 0.032692 | 0.336506 | 0.023159 | 0.345213 | 0.029987 | 0.759944 | 0.047348 | 0.526814 | 0.062030 |\n\n**Question.** Last small question: line 179: what kind of non-optimality is spoken about here?\n\n**Response.** The coordinate-wise median is not optimal with respect to the convergence rate. The rate is of the form $\\sqrt{\\frac{Tr(\\Sigma) \\log \\delta^{-1}}{n}}$ instead of the optimal $\\sqrt{\\frac{Tr(\\Sigma)}{n}} + \\sqrt{\\frac{\\lambda_{\\max} \\log \\delta^{-1}}{n}}$. There are other algorithms which achieve optimality in the convergence rate like the \\emph{median-of-means} tournaments. See the Lugosi and Mendelson survey for a detailed explanation.\n     \n**Question.** I believe the clarity of Lemma 2.3 ... infinitesimally close to it.\n\n**Response.** Algorithm 2 does not always converge to the geometric median of the points. Take, for example, the triangle with coordinates $(0,1),  (0,-1), (1,0)$. The geometric median is the point $(1/\\sqrt 3,0)$. If we initialize the algorithm at the origin, the algorithm never makes any progress and stays at the origin. We wish to reiterate that our goal was only to be sufficiently close to the geometric median. As mentioned in Remark 4.4, the state of the art geometric median algorithms have a significantly higher running time than what we are aiming for. Since even the best possible algorithm converging to the geometric median will run in time $\\Omega(nd \\log \\varepsilon^{-1})$, a running time not achievable with current algorithms and a running time that is already superlinear in our setting, we can only hope to be sufficiently close."}, {"review_id": "UgF3FMJbFt", "reviewer": "Reviewer_Vydg", "summary": "The paper considers the following question: given $n$ points in the $d$-dimensional Euclidean space, we want to find a $(1+\\varepsilon)$-approximation algorithm for finding the mean (where the objective is the $1$-means cost). The goal is to randomly samply certain points and use a fast algorithm to find the approximate mean (with probability at least $1-\\delta$). There were known lower bounds implying lower bound of $\\Omega(1/\\varepsilon \\cdot \\log(1/\\delta))$. There are almost matching upper bounds (with an extra log factor). The paper gives a simple algorithm which samples $\\Omega(1/\\varepsilon \\cdot \\log(1/\\delta))$ points and uses the well-known median of means trick to show that this gives an optimal sample complexity bound. They also give a gradient descent based algorithm which seems to work well on real data sets.", "questions": "1. What is novel in techniques here? The median of means seems to be a standard idea. \n2. What about $\\ell_p$ norms in general? Do the ideas here extend to these more general versions?", "limitations": "yes.", "rating": 3, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths: \n\n1. Simple algorithm achieving optimal sample complexity bound for the 1-means problem. \n2. Validation of algorithms by experiments. \n\nWeakness:\n\n1. Technique is well known and has been applied to many related problems in streaming. \n2. Prior works gave almost optimal sample complexity (upto log error) for all $\\ell_p$ norms, whereas the idea here works only for $\\ell_2$ norm. \n3. Improvement over prior work is only $\\log(1/\\varepsilon)$.", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer’s comments and the opportunity to clarify the contributions and novelty of our work. We respond to the raised questions and comments in a unified manner.\n\nOur primary goal is to achieve optimal runtime in addition to optimal sample complexity for robust mean estimation in streaming settings. While prior works have indeed addressed sample complexity (often for all $\\ell_p$ norms), their algorithms are generally not optimized for runtime. In contrast, our results focus specifically on the $\\ell_2$ norm because our runtime improvements rely on structural properties and lemmas (e.g., Lemmas 2.1 and 2.2) that are unique to this norm. Extending such runtime guarantees to general $\\ell_p$ norms remains an interesting open direction.\n\nRegarding the techniques: although median-of-means is a classical idea in robust estimation, our contributions lie in how we combine and extend known techniques to obtain new algorithmic and theoretical insights. \nThe median-of-means estimator is designed for the one-dimensional mean estimation problem. Extending the idea to higher dimensions is non-trivial. There is no standard notion of a median\nfor multivariate data, and it is not clear a priori what definition of multivariate\nmedian generalizes the median-of-means estimator for one-dimension. \n\nSpecifically, our contributions are:\n\n- We introduce a gradient-descent-based algorithm that starts from the coordinate-wise median-of-means and efficiently converges to the true mean. This process provably runs in linear time, and to the best of our knowledge, this is the first such demonstration in this setting.\n\n- We show that the coordinate-wise median-of-means achieves dimension-free sample complexity in a distributional mean estimation setting—an aspect that, to our knowledge, was overlooked in prior literature.\n\n- In addition, we propose a novel estimator based purely on order statistics, which offers a distinct alternative for mean estimation and runs in near-optimal time.\n\nAs for the claimed improvement: there is also an improvement in by the order of $\\log \\delta^{-1}$ in both sample complexity and running time. While the improvement in terms of the precision $\\varepsilon$ may be more modest (even if it is optimal now and was not before), it can be argued that we obtain a quadratic improvement in both running time and sample complexity in terms of the failure probability. Moreover, breaking the $\\log^2(1/\\delta)$ barrier has been a notable challenge in the literature, and our method addresses this without sacrificing either sample complexity or computational efficiency in other parameters."}, {"review_id": "5mcddbfFII", "reviewer": "Reviewer_sRHN", "summary": "The paper discusses sublinear mean estimation and presents three different algorithms. For each of these algorithms, running time analysis and optimality discussions are provided. Several experiments comparing proposed algorithms with existing algorithms are provided.", "questions": "Majority of my questions are about the proofs and experiments.\n\nQ1. In line 143, it is stated that earlier values of constants were set as $a = 1440, b = 50$, but neither of these appear to be discussed in detail prior to this line. Is some text missing or can some more context be shared as to why these constants were chosen?\n\nQ2. Overall I found the proof of Lemma 2.3 to be quite vague and would prefer it to be explained in more detail, particularly because it is one of the cornerstone results of the paper. For instance, why do we get $\\|\\hat{\\mu}_i-\\mu\\|^2 > r^2$ with probability at most 0.1, and where does this 0.1 come from? How does $-2/81$ appear in line 145-146?\n\nQ3. In line 160 it is stated that \"$L_k, R_k$ have cardinality at least $b\\log \\delta^{-1}/2$.\" What does this mean when $L_k + R_k = b\\log \\delta^{-1}$? Following this logic, can you explain in a bit more detail how we can infer at least $1/5$ of the samples satisfy $|\\hat{\\mu}_{i,k}-\\mu_k|\\geq |\\nu_k - \\mu_k|$ from this?\n\nQ4. Following Q3, can you clarify the logic in lines 160-161? In particularly, since $L_k, R_k$ are sets, the sentence \"at least for one of $L_k,R_k$, we have $|\\hat{\\mu}_{i,k}-\\mu_k|\\geq |\\nu_k - \\mu_k|$\" doesn't make a lot of sense.\n\nQ5. Following Q4., because of Q4 it's quite difficult to judge if the sentence \"Depending on whether $\\nu_k>\\mu_k$ or not, at least for one of $L_k,R_k$, we have that $|\\hat{\\mu}_{i,k}-\\mu_k|\\geq |\\nu_k - \\mu_k|$.\" is logical. Wonder if you can help clarify it with some examples, for example say if $\\mu_k = 0, \\hat{\\mu}_k = 2, \\nu_k = 3$. In this case, the left hand side is 2 while the right hand side is 3 and the inequality doesn't hold.\n\nQ6. As briefly mentioned in an earlier section of the review, it appears from the experiments that the empirical mean performs consistently better than proposed algorithms, and is likely fast as well; can you share more insights on how this may be happening?\n\nQ7. Wonder if there are more experiments that look at how the algorithms perform under different data dimension, say, varying between 1- 1000 instead of a fixed value like 784 or 54.\n\nAdditional comments regarding formatting and typos:\n\n1. Line 65: \"three algorithm\"\n2. Line 109: \"studided\"\n3. Line 135: \"appoximation\"\n4. Line 191: \"proof for the our setting\"\n5. Line 313: \"runing time\"\n6. Line 330: \"these observation\"", "limitations": "The authors adequately addressed limitations regarding uniform sampling and data corruption.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 4, "originality": 3, "strengths_and_weaknesses": "1. Perhaps most importantly, the new proposed algorithms are optimal in many aspects such as scalability and efficiency, which is arguably the highlight of the paper. There are some follow up questions regarding proof details which are deferred to the \"Questions\" section of the review, but if shown to be correct, the results are fairly interesting.\n2. Overall, the paper is well-organized and discusses background, key concepts, and contributions well. It is clear that the authors spent a lot of effort polishing the paper. Several suggestions regarding typos or formatting issues are deferred to a later section.\n3. The experiments are relatively weak. For instance, it is clear that the empirical mean consistently outperforms proposed algorithms in terms of accuracy. While the empirical mean is not included in *runtime* comparisons, it should be significantly faster than the proposed algorithms, if I understand correctly. Perhaps some additional experiments or discussion would be helpful to address these observations.\n4. Following point 3., I think the result of Theorem 3.2 regarding learning the mean of high-dimensional distributions is also quite interesting but there doesn't appear to be any experiments that support the dimension-free result. Might be helpful if some extra experiments or discussions can be made there.", "was_revised": true, "final_justification": "My original concerns about details of the proofs were addressed sufficiently. Overall I think the paper has good quality and it will improve with further revisions. After engaging with other reviewers, I agree about some concerns regarding novelty. Based on this I'm suggesting a score of 4.", "rebuttal": "We thank the reviewer for the helpful feedback. As per NeurIPS 2025 policy, we cannot include links or images in the rebuttal or repository. We’ve therefore provided tables (in earlier responses) with the same information as the plots and appreciate your understanding given these limitations.\n\n**Regarding proofs.**\n\n**Q1 Response.** We apologize for the confusion. The sentence has a typo and should instead be \"We set the values of the constants to be \\dots \". In the answer to the next question, please find a detailed response including motivations behind the values of the constants.\n\n**Q2 Response.** For the algorithms, the first step is to sample $b \\log \\delta^{-1}$ subsamples of size $a \\varepsilon^{-1}$ each and compute their means. We want our final estimate to be within $\\|\\hat \\mu - \\mu \\| \\leq \\sqrt{\\frac{\\varepsilon OPT}{n}}$ distance of the true mean, as it implies a $(1+\\varepsilon)$-approximation. Define $r = \\frac{1}{11} \\sqrt { \\frac{\\varepsilon OPT}{n}}$ and  we call a mean \\emph{good} if $\\|\\hat \\mu - \\mu\\| \\leq r$. Further, we define the good event $\\mathcal{E}$ to be when at least $7/10$ fraction of the means are \\emph{good}. The value $7/10$ is chosen to an arbitrary constant strictly greater than $1/2$.\n\n**Lemma 2.3.** Event $\\mathcal{E}$ holds with probability at least $1-\\delta$.\n\nThe proof contains the following steps:\n1. A sample mean derived from $1440 \\varepsilon^{-1}$ samples is \\emph{good} with probability $0.9$.\n2. If we have $50 \\log \\delta^{-1}$ subsamples of $1440 \\varepsilon^{-1}$ points each, at least $7/10$ fraction of the mean will be \\emph{good} with probability $1-\\delta$.\nWe now give a detailed version of the proof with an explanation of the constants as well. For now, the values of constants are set as follows - $a = 1440, b = 50$.\n\n**Proof.** From Lemma 2.1 in the paper, we know that $ \\mathbb{E} [ \\| \\mu - \\hat \\mu \\|^2 ]= \\frac{1}{|S|} \\frac{OPT}{n} = \\frac{1}{a \\log \\varepsilon^{-1}} \\frac{OPT}{n}$. The Markov's inequality can be used to bound the probability that a single mean is \\emph{not good} - $\\mathbb{P}[\\| \\mu - \\hat \\mu \\|^2 > r^2] \\leq \\frac{r^2}{\\mathbb{E} [ \\| \\mu - \\hat \\mu \\|^2 ]} \\leq 0.1$. Conversely, the probability that a single mean is \\emph{good} is at least 0.9. We require this probability to be strictly greater that $0.7$ in order to be able to say that the good event $\\mathcal{E}$ holds with high probability. Specifically, the value of $a$ is set to $1440$ to achieve this.\n\nNext, we use the (multiplicative) Chernoff bound to show that if we have $b \\log \\delta^{-1}$ means, then with high probability, at least 0.7 fraction of the means are \\emph{good}.  \n\nThe Chernoff bound says $\\mathbb{P}[X < (1 - \\theta)\\mu] \\leq \\text{exp}(-\\theta^2\\mu/2)$, where $X$ is the sum of independent random variables and the expected value $\\mu = \\mathbb{E}[X]$. Setting $\\mu \\geq 0.9 b \\log \\delta^{-1}$, $\\theta = 2/9$, we get\n\n$\\mathbb{P}[|G| \\leq 0.7 b \\log \\delta^{-1}] \\leq \\text{exp}(-\\frac{2}{81} b \\log \\delta^{-1}).$\n\nThus, the constant $-2/81$ is the result of applying the Chernoff bound by setting the parameters appropriately. We choose the value of $b$ so that $-(2b)/81 < -1$ and the probability bound achieved is at most $\\delta$. Therefore, we set $b=50$.\n\nWe will add a more detailed version of the lemma and its proof to the final version of the paper.\n\n**Q3 Response.** We have $|L_k \\cup R_k| = b \\log \\delta^{-1}$ which is the total number of sample means. Recall that $L_k$ consists of the sample means such that $\\hat \\mu_k \\leq \\nu_{CWM,k}$. And $R_k$ is the set of sample means such that $\\hat \\mu_k \\geq \\nu_{CWM,k}$. The key detail is that there are points (such that $\\hat \\mu_k = \\nu_{CWM,k}$) which belong to both. Along with the definition of the coordinate-wise median, we have $|L_k|, |R_k| \\geq b \\log \\delta^{-1}/2$. Note that if there are no points such that $\\hat \\mu_k = \\nu_{CWM,k}$, then $|L_k|= |R_k| = b \\log \\delta^{-1}/2$.\n\nWe provide a detailed case analysis as to why at least $1/5$ fraction of the samples are good and satisfy $|\\hat \\mu_{i,k} - \\mu_k| \\ge |\\nu_{CWM,k} - \\mu_k|$.\n\nFirst, we note that when the good event $\\mathcal{E}$ holds,  at least $1/5$ fraction of the means are *good* in both the sets $L_k$ and $R_k$. This follows from the fact that at least $7/10$ fraction of the means are good and the number of good means in $L_k$/$R_k$ is at least $7/10 - 1/2 = 1/5$.\n\n**Case I:** $\\mu_k \\leq \\nu_{CWM,k}$. In this case, we have for all means in $R_k$, $|\\hat \\mu_{i,k} - \\mu_k| \\ge |\\nu_{CWM,k} - \\mu_k|$.\n\n**Case II.** $\\mu_k > \\nu_{CWM,k}$. In this case, we have for all means in $L_k$, $|\\hat \\mu_{i,k} - \\mu_k| \\ge |\\nu_{CWM,k} - \\mu_k|$.\n\nFrom the above statements, we have for at least one of $L_k$ and $R_k$, we have that at least $1/5$ fraction of the means are good and satisfy $|\\hat \\mu_{i,k} - \\mu_k| \\ge |\\nu_{CWM,k} - \\mu_k|$.\n \n**Q4 Response.** We apologize for the typo. We want to say that the claim is true for all means in $L_k$ or $R_k$. The formal statement is as follows. At least one of the following statements is true:\n1.  For all $\\hat \\mu$ in $L_k$, we have that $|\\hat \\mu_{i,k} - \\mu_k| \\ge |\\nu_{CWM,k} - \\mu_k|$.\n2. For all $\\hat \\mu$ in $R_k$, we have that $|\\hat \\mu_{i,k} - \\mu_k| \\ge |\\nu_{CWM,k} - \\mu_k|$.\n\nWe will update the final version of the paper accordingly.\n\n**Q5 Response.** Please refer to our answer to Q3 for the detailed proof. In the example given by the reviewer, $\\nu_{CWM,k} = 3$, implies that there are equal number of points with $\\hat \\mu_k < 3$ and $\\hat \\mu_k > 3$. So, if there exists a sample mean such that $\\hat \\mu_k = 2 < 3$, then there will also exist at least one mean such that $\\hat \\mu_k > 3$. And for that mean, which is in $R_k$, we have $|\\hat \\mu_{i,k} - \\mu_k| \\ge |\\nu_{CWM,k} - \\mu_k|$.\n\n**Regarding experiments.**\n\n**Q6 Response.** As the reviewer correctly notes, the empirical mean is the fastest among all mean estimators, as it can be computed in closed form. Despite this computational advantage, it is not a robust estimator of the true mean—even in one dimension. That said, the empirical mean performs well when estimating the mean of subgaussian distributions, where the associated tail bounds are sharply exponential. In contrast, for heavy-tailed distributions, the best available tail bound is Chebyshev’s inequality, which is exponentially weaker. In fact, Catoni (2012) shows that there exist distributions for which Chebyshev’s inequality is tight. Nevertheless, in practice, data distributions are often not as adversarial as those constructed to demonstrate the failure of the empirical mean (e.g., Catoni (2012)). This explains why the empirical mean is frequently observed to perform well in real-world scenarios, both in terms of computational efficiency and estimation accuracy.\n\n**Q7 Response.** We thank the reviewer for the suggestion. We have run two experiments: They illustrate how accuracy and runtime vary with dimensionality. Specifically, we generated a synthetic dataset of 10,000 samples drawn from multivariate Gaussian distributions with dimensions ranging from 10 to 1000. As predicted by theory, the errors of both the coordinate-wise median and the fast gradient descent estimator remain low and stable across dimensions. Interestingly, the CSS algorithm shows rapid improvement as the dimension increases, stabilizing around dimension 500. As expected, the empirical mean achieves the best accuracy since the data are drawn i.i.d. from a Gaussian distribution. In terms of runtime, we observe a roughly linear growth with dimension, in line with theory. Unfortunately, we could not insert these tables due to characters constraints, but we will add these new experimental results (and the corresponding plots) to the final version of the paper."}]

Author rebuttal:

Meta-review: This paper studies the multivariate mean estimation problem. It shows an aggregate estimator that averages mean estimates on bootstrap samples (with independent and uniform sampling) can output an $(1+\varepsilon)$-approximate mean with probability at least $1 - \delta$ by using  $O(\varepsilon^{-1} \log \delta^{-1})$ sampled points. This is an improvement of the classic  $O(\varepsilon^{-1} \delta^{-1})$ sample complexity and the SOTA $O(\varepsilon^{-1} \log \delta^{-1} (\log \varepsilon^{-1} + \log \delta^{-1}))$ sample complexity. The paper further presents three efficient algorithms that run the aggregate estimator in $O(\varepsilon^{-1} \log \delta^{-1})$ time. Experimental results on a real-world data set shows one of the proposed estimators achieve the best accuracy-runtime tradeoff, while the other two have competitive performance compared to the SOTA method. 

--- 

The paper makes a plausible improvement on SOTA sample complexity for mean estimation, and provides efficient estimators. 

What holds reviewers Vydg and sRHN back is the limited technical novelty. Reviewer Vydg also points out the improvement only applies to the $\ell_2$ norm while prior results apply to $\ell_p$. Another notable weakness is the experiment section. Concerns include the limited scale of the data set (reviewers otjC & 1w8V & sRHN), insufficient comparison with empirical mean as a baseline (reviewer sRHN) and missing verification of the dimension-free property (reviewer sRHN). 
 
Despite those concerns, most reviewers are on the positive side. I think the work is decent and the paper is well written. Its contribution may have a relatively narrowed scope, but its merits outweigh this limitation. Authors should incorporate all reviewers' feedback in revision.
Decision: Accept (poster)