Efficient Multivariate Robust Mean Estimation Under Mean-Shift Contamination

We thank the reviewer for their effort and time in assessing our work. We respond to the points raised individually below:

(Experiments) We would like to emphasize that our primary contribution is to characterize the computational-statistical landscape for this fundamental learning task—in terms of error guarantee, sample complexity, and runtime. That being said, our algorithm is already fairly simple and potentially practical, as it uses a dimension-reduction procedure with each step essentially being an SVD of a reweighted covariance matrix. Also, while it is easy to generate synthetic Gaussian data, it is not clear what distribution on errors should be considered, making the meaningfulness of experiments questionable.

(Novelty) We would like to stress that our work develops a novel spectral dimension-reduction technique, rather than applying/adapting methods from prior work. To the best of our knowledge, there are two “spectral methods” proposed in prior work for robust mean estimation (in Huber’s model). Both of them are significantly different than our algorithm, as we discuss below:

• (Spectral filtering from Diakonikolas et. al (2016)) This method consists of iteratively computing the top eigenvector of the empirical covariance and removing points whose projected deviation from the mean along that direction is higher than a carefully chosen threshold. This has the effect of removing more inliers than outliers until the dataset becomes so clean that the empirical mean is accurate. Importantly, this iterative filtering approach cannot yield consistent estimation (for the mean-shift model that we study), essentially because it does not take into account the additional structure of the outliers. The algorithm we propose is fundamentally different. Rather than removing points based solely on the top eigenvector of the empirical covariance, our approach performs dimension reduction. It uses multiple eigenvectors of a reweighted covariance matrix as the subspace of the next iteration. In our method, each point is assigned an exponential weight, with carefully designed reweighting factors that are updated throughout the execution of the algorithm (see lines 126-144, second column and 175-186, first column for intuition on choosing the reweighting factor). Finally, once the problem’s dimension has been reduced sufficiently, we employ an inefficient estimator on the lower-dimensional data.

• (Dimension reduction from Lai et al. (2016)) Although the goal in that technique is also to reduce the dimension, as we explain in lines 662-671 in Appendix A of the paper and in our response to Reviewer BKkQ, the similarities are only superficial. The method from Lai et al. (2016) uses the standard covariance matrix to identify the subspace for the next iteration. This method also fails to leverage the special nature of the outliers. In contrast, our method reweights the space with an exponential function that encodes our prior knowledge on the outliers in order to truncate their effect in the second moment matrix. Although one could broadly view this reweighting as a kind of “soft filtering”, this is fundamentally different from the filtering of the previous paragraph: there, points were removed based solely on their deviation along a single direction, where here the weight assigned to each point is based on the norm of the vector. Moreover, as we comment in lines 662-671 of the paper and in our response to Reviewer BKkQ, our dimension reduction comes with a number of technical obstacles that are due to the nature of our reweighting and need to be carefully addressed.

(Dependence on $\alpha$) Theorem 1.2 holds for any $\alpha \leq 0.49$. As we reply to Reviewer oS19, the most interesting parameter regime in our opinion is when the fraction of outliers $\alpha$ is a positive constant. Essentially, this is the hardest setting for efficient estimation. In that regime, no sub-exponential time algorithm was known to achieve error lower than $\Omega(1)$, while our algorithm can obtain arbitrarily small error in sample-polynomial time; and its sample complexity is information-theoretically near-optimal (as follows from the lower bound of Kotekal & Gao, 2025). As we show in Appendix F, the parameter $0.49$ in the condition of Theorem 1.2 can be replaced by any other absolute constant strictly smaller than $½$ (and the error will degrade by an absolute constant factor). Furthermore, if one prefers to use the reparameterization $\alpha = ½ - c$, where $c$ is a parameter (rather than an absolute constant), the error degrades in a way that increases as $c$ approaches zero. We refer to the guarantee in Claim F.2 for the precise functional form of this dependence, which quantifies how the problem becomes more difficult as $\alpha$ increases.

We thank the reviewer for their time and for their positive assessment of our work. Below, we provide a clarification in response to their comment.

The introductory proof sketch mentions that we can assume $\lVert\mu\rVert = O(1)$ due to naive outlier removal. This means that even if $\lVert\mu\rVert$ is arbitrary, we can preprocess the data to reduce it to the case where $\lVert\mu\rVert = O(1)$; a step that is formally carried out in the subsequent sections containing the full proofs. The argument is quite simple: At the start of the algorithm (line 6 of Algorithm 1), we perform a preprocessing step to obtain a rough estimate $\hat{\mu}_0$ satisfying $\lVert\hat{\mu}_0 - \mu\rVert = O(1)$. This can be achieved using any black-box outlier removal estimator designed for the Huber model from prior work. We then draw a new dataset and replace each point $x$ in that dataset by $x - \hat{\mu}_0$. This transformation effectively shifts the data distribution so that its true mean becomes $\mu - \hat{\mu}_0$, which has norm $O(1)$. After the algorithm produces an $\epsilon$-approximation $\hat{\mu}_1$ for this shifted mean, we simply add back $\hat{\mu}_0$ to $\hat{\mu}_1$, obtaining an $\epsilon$-approximation of $\mu$ (line 26 of Algorithm 1).

We thank the reviewer for their effort and their positive assessment of our work. We respond to the comments below:

• (Comparison with Lai, Rao, and Vempala 2016) As we discuss in detail in Appendix (lines 662–671), although Lai, Rao, and Vempala (2016) iteratively reduce the dimension using the centered second moment of the data, this approach can only achieve an error of order $\alpha$, at best (in fact, their algorithm additionally incurs error with a logarithmic dependence in the dimension). Our dimension-reduction technique is fundamentally different, with the key innovation being the use of exponential reweighting for each sample. This mitigates the impact of outliers (in our model) and enables consistent estimation. Implementing this approach introduces additional challenges that must be addressed: (i) our dimension reduction method requires the space to have at least $1/\epsilon^2$ dimensions (as opposed to just one dimension in Lai, Rao, and Vempala 2016); (ii) the parameter used in exponential reweighting influences both the centering applied in the computation of the centered second moment and the sample complexity. As explained in Section 1.2, the first attempt for setting this parameter correctly does not immediately allow for reduction all the way down to the desired level. Instead, our algorithm operates in two stages: the first stage performs the initial dimension reduction from the last sentence, and the second stage resets that parameter to a refined value that allows for further reduction.
• (Runtime) As the reviewer points out, the runtime of $2^{O(1/\epsilon^2)}$ is unavoidable, even in one dimension. This is because, even in the one-dimensional case, estimating the mean up to error $\epsilon$ information-theoretically requires $2^{\Omega(1/\epsilon^2)}$ samples, as shown in Kotekal & Gao (2025). The key conceptual point of our paper is that, prior to our work, it was unknown whether the high-dimensional problem incurs a significantly higher runtime (e.g., $2^{\Omega(d/\epsilon^2)}$). We rule out this possibility by showing that runtime $poly(d)2^{\Omega(1/\epsilon^2)}$ is enough (i.e., the dependence on the dimension is polynomial instead of exponential). Notably, this implies that setting $\epsilon = \Theta(1/\sqrt{\log d})$ allows for error tending to zero as $d \to \infty$) within an overall $\text{poly}(d)$ runtime. In contrast, every prior algorithm for the Huber model cannot achieve an error smaller than $\Omega(\alpha)$, which could be an absolute constant.

We thank the reviewer for their time, effort and their clarifying questions. We comment on the points mentioned and respond to the reviewers questions below:

(Identity covariance) As pointed out by the reviewer, our algorithm works under the assumption that the covariance of the inliers is known. Please note that our algorithm does not require the covariance to be exactly equal to the identity; it is sufficient for the covariance to be known. (As pointed out in lines 85-88, for known covariance, by applying an appropriate transformation to the samples, we can reduce the problem to the identity covariance case.) We study the known covariance case, as it is the most fundamental setting for which we demonstrate the main conceptual point of the paper, which is an examination of how more structured adversaries in robust statistics can alter the statistical and computational landscape of estimation. As discussed in Section 3, extending our algorithmic result to the case of unknown covariance is an interesting problem for future work. There are several nuances in defining the problem for the unknown covariance case. For the sake of space, we point the reviewer to lines 629-654 for a thorough discussion.

(Dependence on alpha). We discuss the two relevant points below:

• (Dependence in the main theorem). The result in Theorem 1.2 holds for any $\alpha < 0.49$ (where $0.49$ can be replaced by any constant near $1/2$, as shown in Appendix F). The main conceptual point of the theorem is that even when nearly half of the samples are outliers, the algorithm efficiently achieves any desired error $\epsilon$ with the stated sample complexity of roughly $d/\epsilon^{2 + o(1)} + 2^{O(1/\epsilon^2)}$. This sample complexity is tight in this setting. Any prior algorithm would either require exponential runtime or incur a significantly larger error of $\Omega(1)$ (see lines 72-80 of the paper where this is explained), whereas our result allows for arbitrarily small error. Regarding other regimes for $\alpha$, prior work on one dimension (Kotekal & Gao (2025)) has fully specified the sample complexity as both a function of $\alpha$ and $\epsilon$. It would be interesting to investigate the exact dependence in high-dimensions. However, as with most prior work on robust statistics, the most interesting regime is to obtain a result for when the fraction of outliers $\alpha$ is a positive constant. This is because this is the hardest scenario for estimation (with the smallest amount of information about the target mean).
• (Special cases for $\alpha$). Formally speaking, setting $\alpha = \epsilon$ is insufficient to obtain error $\epsilon$ using a Huber robust estimator, as any estimator for the Huber model incurs an error of at least $C\epsilon$ for some constant $C > 1$ (as established in the robust statistics literature). Ignoring this nuance regarding absolute constants, the regime $\alpha = \Theta(\epsilon)$ (and even more so $\alpha = 0$) corresponds to special cases where existing algorithms can be directly applied to achieve error $O(\epsilon)$—specifically, any algorithm designed for Huber contamination or the sample mean setting, respectively. Thus, these cases are orthogonal to our work as they are special cases that can be handled by existing algorithms. As mentioned earlier, the main question we focus on is obtaining a result for the case when $\epsilon \ll \alpha$, i.e., we aim to estimate with error significantly smaller than the fraction of outliers. This is an objective that was impossible to achieve with algorithms designed for the Huber model, and our algorithmic result for constant fraction of outliers has tight sample complexity.