Simple and Optimal Sublinear Algorithms for Mean Estimation

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.

Regarding proofs.

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.

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 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 good. The value $7/10$ is chosen to an arbitrary constant strictly greater than $1/2$.

Lemma 2.3. Event $\mathcal{E}$ holds with probability at least $1-\delta$.

The proof contains the following steps:

1. A sample mean derived from $1440 \varepsilon^{-1}$ samples is good with probability $0.9$.
2. If we have $50 \log \delta^{-1}$ subsamples of $1440 \varepsilon^{-1}$ points each, at least $7/10$ fraction of the mean will be good with probability $1-\delta$. We 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$.

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 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 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.

Next, 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 good.

The 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

$\mathbb{P}[|G| \leq 0.7 b \log \delta^{-1}] \leq \text{exp}(-\frac{2}{81} b \log \delta^{-1}).$

Thus, 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$.

We will add a more detailed version of the lemma and its proof to the final version of the paper.

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$.

We 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|$.

First, 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$.

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|$.

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|$.

From 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|$.

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:

1. For all $\hat \mu$ in $L_k$, we have that $|\hat \mu_{i,k} - \mu_k| \ge |\nu_{CWM,k} - \mu_k|$.
2. For all $\hat \mu$ in $R_k$, we have that $|\hat \mu_{i,k} - \mu_k| \ge |\nu_{CWM,k} - \mu_k|$.

We will update the final version of the paper accordingly.

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|$.

Regarding experiments.

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.

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.

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.

Our 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.

Regarding 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. The 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 for multivariate data, and it is not clear a priori what definition of multivariate median generalizes the median-of-means estimator for one-dimension.

Specifically, our contributions are:

• 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.

• 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.

• 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.

As 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.

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.

Comment. A sample size $m$ is introduced ... missing from the paper.

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.

Comment and Question. Three fixed datasets ... Gaussian/uniform sample; While this is ... much larger sample sizes.

Response. Please refer to the last response to reviewer otjC and tables therein.

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.

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).

Comment. Runtime barely grows ... sample size is important.

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.

Table 1: MNIST — Accuracies vs. Sample Size

Table 2: MNIST — Runtimes vs. Sample Size

Question. Last small question: line 179: what kind of non-optimality is spoken about here?

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 median-of-means tournaments. See the Lugosi and Mendelson survey for a detailed explanation.

Question. I believe the clarity of Lemma 2.3 ... infinitesimally close to it.

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.

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.

Question. The idea ... the algorithm itself.

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.

Comment. The algorithm looks like ... reject this nice paper.

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.

Comment. The experiments are very poor ... random numbers would do.

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$.

The 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).

Unfortunately, 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. We 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. We will add these new experimental results to the final version of the paper.

Table 3: Synthetic Dataset ($d = 200$) — Accuracies vs. Sample Size

Table 4: Synthetic Dataset ($d = 200$) — Runtimes vs. Sample Size