Near-Optimal Algorithms for Gaussians with Huber Contamination: Mean Estimation and Linear Regression
Authors: Ilias Diakonikolas
Daniel Kane
Ankit Pensia
Thanasis Pittas
Venue: NeurIPS 2023 poster (2023)
Keywords: robust statistics
high-dimensional inference
regression
nearly linear time algorithms
TL;DR: We provide algorithms for Gaussian mean estimation and linear regression under Huber contamination that achieve optimal error and sample complexity and have almost linear runtime

Abstract: We study the fundamental problems of Gaussian mean 
estimation and linear regression with Gaussian covariates 
in the presence of Huber contamination. Our main 
contribution is the design of the first sample near-optimal 
and almost linear-time algorithms with optimal error 
guarantees for both these problems. Specifically, for 
Gaussian robust mean estimation on $\mathbb R^d$ with 
contamination parameter $\epsilon \in (0, \epsilon_0)$ for a small 
absolute constant $\epsilon_0$, we give an 
algorithm with sample complexity $n = \tilde{O}(d/\epsilon^2)$ 
and almost linear runtime that approximates the target 
mean within $\ell_2$-error $O(\epsilon)$. 
This improves on 
prior work that achieved this error guarantee with 
polynomially suboptimal sample and time complexity. 
For robust linear 
regression, we give the first algorithm with sample 
complexity $n = \tilde{O}(d/\epsilon^2)$ and almost linear 
runtime that approximates the target regressor within 
$\ell_2$-error $O(\epsilon)$. This is the first polynomial 
sample and time algorithm achieving the optimal error 
guarantee, answering an open question in the literature. 
At the technical level, we develop a methodology that 
yields almost-linear time algorithms for multi-directional 
filtering that may be of broader interest.

Near-Optimal Algorithms for Gaussians with Huber Contamination: Mean Estimation and Linear Regression
Ilaas Diakonikolas University of Wisconsin-Madiion ilias@cs. wisc visc
Daniel M. Kane
University of California, San Diego dakane@cs. ucsd. edu
Ankit Pensia IBM Research
Thanasis Pittas University of Wisconsin-Madisn pittasewisc. edu
ankit@ibm. com
Abstract
We study the fundamental problems of Gaussian mean estimation and linear regression with Gaussian covariates in the presence of Huber contamination. Our main contribution the design of the first sample near-optimal and almost linear-time algorithms with optimal error guarantees for both these problems. Specifically, for Gaussian robust mean estimation on Rd with contamination parameter E (0, eo) for small absolute constant 6o, we give an algorithm with sample complexity the target mean within eerrr O(e). This improves on prior work that achieved this error guarantee with polynomially suboptimal sample and time complexity. For robust linear regression, we give the first algorithm with sample complexity n = Ó(dee) and almost linear runtime that approximates the target regressor within -2-error O(e). This is the first polynomial sample and time algorithm achieving the optimal error guarantee, answering an open question in the literature. At the technical level, develop methodology that yields almost-linear time algorithms for multi-directional filtering that may be of broader interest.
Ódd/²) and almost linear runtime that approximates
Introduction
Modern machine learning systems operate with vast amounts of training data, which are difficult to carefully curate. Consequently, outliers have become a fixture modern training datasets. This contra- dicts the standard i.i.d. assumption of classical statistical theory and has spurred the development of ro- bust statistics, which seeks to develop algorithms that perform well in the presence of outliers (HR09; DK23J. The standard contamination model of outliers was first formalized by Huber [Hub64]:
Definition 1.1 (Huber Contamination Model). Given 0
1/2 and a distribution family D, the algorithu specifies E N and observes n i.i.d. samples from a distribution P = (1 - e)G + eB, where G is an unknown distribution D, and B is an arbitrary distribution. We say G is the distribution of inliers, B is the distribution of outliers, and P is an e-corrupted version of G.
The Huber contamination model has since served as bedrock for the development and evaluation of robust algorithms. Given that estimating properties Gaussian distributions prototypical task in statistics, Gaussian estimation under Huber contamination is, analogously, central problem in robust statistics. The univariate version of this problem addressed in Huber's 's work [Hub64]. In this paper, we focus estimating digh-simm Gaussian distributions under Huber contamination.
Authors are in alphabetical order.
37th Conference on Ncural Information Processing Systems (NcurIPS 2023).

---

Consider the fundamental problem of estimating the mean of d-dimensional isotropic Gaussian I). given samples from an e-corrupted distribution (Definition 1.1). [DKKLMS18] gave an = poly(d/e) sample and poly(n) (time algorithm for this problem, achieving the information- theoretically optimal error of ( Despite being polynomial in sample and time complexity, these guarantees are far from the linear sample complexity and linear runtime of the sample mean on uncontaminated data. While a number of works [CDG19; DHL19; DKPP22] have developed near-linear time Gaussian robust mean estimation algorithms, all these prior methods inherently Given the fundamental nature of Gaussian
suffer sub-optimal error guarantee of (lcc//loo//e)) robust mean estimation [DM22], these contrasting sets of results raise the question of whether is possible to achieve the best of both worlds. In other words:
Can we obtain (+)) error for Gaussian mean estimation with Huber contamination in near-linear time and sample complexity?
While mean estimation is the most basic unsupervised learning task, linear regression is arguably the most basic supervised learning task. Here we study the basic case of Gaussian covariates.
Definition 1.2 (Gaussian Linear Regression). Fix a > O and 3 E R.. Let G be the joint distribution of pairs (X, y), with X E R y E R, such that X ~ N (0, Ia) and 3 X + 17, where ~ N o²) independently of X. The goal of the algorithm is compute 3 such that ||3 3 is small.
The theeroi error for robust Gaussian linear regression with Huber contamination is @(re). However, all known polynomial time algorithms for this task incur higher error of ((oe//og(1/e) [DKS19; PJL20; CATJFB20; BP21], raising the following open question in [BP21]:
Can we obtain ((ae) error for Gaussian linear regression with Huber contamination in polynomial time and sample complexity?
Perhaps surprisingly, despite the extensive algorithmic progress on robust statistics over the past years [DKKLMS16; LRV16; DK19; DK23), these basic questions have remained open.
1.1 Our Results
Robust Mean Estimation We begin by stating our result for Gaussian robust mean estimation:
Theorem 1.3 (Almost Linear-Time Algorithm for Robust Mean Estimation). Let en be a sufficiently small positive constant. There is an algorithm that, given parameters e E (0. <a). e E (0, 1), E (0,1), and n > 3 (d + log(1 (1)) polylog(d/c) e-corrupted points from N(,,, La) (in the Huber contamination model per Definition 1.1), computes an estimate ÿt such that [| 1ll2 || O(e/c) with probability at least 8. Moreover, the algorithm runs in time (nd + 1/c² (ooyylooolyyy
Taking e to be small positive constant, our algorithm achieves the optimal asymptotic error (see, e.g., (CGR[C) up to constant factor with near-optimal sample complexity (CGR18) and almost-linear runtime. Moreover, the runtime of our algorithm near-optimal in the regime when e 2 d 2/c which is the main regime of interest." Moreover, we note that the algorithm continues to work for a wider family of distributions than Gaussians; See Remark 2.2.
Our techniques are also amenable to the streaming framework of [DKPP22), and Theorem 1.3 can be extended to a streaming algorithm (with polynomial time and sample complexity) that uses only Ö(d + poly(1/e) (1/e))-memory, -memory, which is optimal up to the additive poly(1/e) term (see Appendix E).
Robust Linear Regression We next consider robust Gaussian linear regression (Definition 1.2). As existing polynomial-time algorithms for robust Gaussian linear regression of O( log((1e)) (1/e)). we assume without loss of generality that the true regressor satisfies || = (oe (1/E)) (in fact, it can be ensured in nearly-linear time; see, c.g., [CATJFB20, Theorem 2.5]).
achieve error
Using a novel reduction of linear regression to robust mean estimation for Gaussians, we prove:
These prior algorithms continue work in the total variation/strong contamination model with error (l (lo//e)()/e)) Still, their underlying algorithmic approaches provably cannot obtain better eroor in the Huber model - the fundamental contamination model and the focus of our work; See Section 1.3 for details.
Using the expression n = d/e², the runtime is dominated by nid when > n -2/(4+o), Note that the sample complexity requirement directly implies that > 1/

---

Theorem 1.4 (Almost Linear-Time Algorithm for Robust Linear Regression). Let Eo be suffi- ciently small positive constant. Let G be the joint distribution over (X, y) following the Gaussian linear regression model (Definition 1.2) with unknown parameters a > 0 and 33 E I2 satis- fying |/1|| = Oge log(1/e)) There is an algorithm that, given e E co). C E (0.1) and n > (oly))ooy (d/)) i.idd labeled examples {x, y) from an e-corrupted version G, the al- gorithm returns a 3 such that ||3 - B12 = O(oe/c), with probability at least 0.9. Moreover, the algorithm runs in time (nd +
For a constant e > O, Theorem 1.4 has optimal asymptotic error [CGR16], near-optimal sample com- plexity [CGRI6] for constant failure probabilities, and near-optimal runtime (up to an additive factor of 1/e"). Thus, we provide the first polynomial time and sample algorithm (in fact, near-optimal time and sample complexity) for robust Gaussian linear regression with Huber contamination.
1.2 Our Techniques
Our first major result is near-optimal sample and almost-linear time algorithm for learning Gaussian mean to error (c) in the Huber contamination model. At a high level, our technique borrows ideas from the (c) error algorithm of [DKKLMS18 and the fast robust estimation techniques of [DKPP22]. We emphasize that a number of challenges need to be overcome in order to achieve this result, as elaborated below.
Roughly speaking, the fast algorithm of [DKPP22] works via careful implementation of the standard filtering algorithm [DKKLMS16; DKKLMS17: DK19]. If the empirical covariance matrix has no direction of large variance, this serves as certificate that the remaining outliers do not substantially affect the mean. If the empirical covariance matrix has any large eigenvalue, the algorithm projects all the points onto a (randomly chosen) direction of large variance and removes outliers in this direction. Moreover, if there is such a direction, careful analysis can be used to show that the filtering removes more outliers than inliers and that improves a carefully chosen potential function covariance matrix Σ', roughly tr((E' - -I) ). Repeating this procedure only polylog(d/e (d/e)) many times eventually yields sample set with no directions of large variance, whose sample mean is guaranteed to be close to the true mean. Crucially, to achieve this fast runtime, [DKPP22] must consider rundom directions of large variance obtained by applying suitable power of the (translated) empirical covariance matrix to a random vector. This is to ensure that the outliers cannot be arranged so that they will end up orthogonal to the directions considered.
based on the trace of an appropriate power of the (translated) empirical
Conceptually, the (cr--rrro algorithm in Huber's model [DKKLMS18] works via a more compli- cated filter whereby if the covariance matrix has k (for k = (2log((1cc) moderately large eigenvalues 1 1+ (()). the algorithm projects all points onto the subspace spanned by these k directions and re- moves those points whose projections are too far from the mean. The usage of multiple orthogonal directions permits better concentration bounds of Gaussians (namely, the Honnonn-rrigh inequality) and achieves stronger filtering that translates to smaller final error. In particular, the standard filter cannot reduce the leading eigenvalue to + o( log(1/e)), while this improved multi-directional filter reach the threshold + O(e) for the k-th largest eigenvalue. A brute-fore approach is then used to Iearn the projection of the mean onto the subspace of the & largest eigenvalues, while the sample mean is used an approximation to the mean in the orthogonal directions. 4. However, [DKKLMS] MS88] takes quadratic time, nd², to run because the subspace used to remove points is deterministic and outliers may be arranged in way that filtering does not remove the outliers that lie in the orthogonal subspace. A natural idea is to use insights from [DKPP22] to improve the algorithm in [DKKLMS18]
Unfortunately, combining these techniques is highly non-trivial. While [DKPP22] filters based on a single random direction, [DKKLMS 18] requires a multi-directional filter with (loog(1/))) many (nearly-) orthogonal directions. If the technique of [DKPP22] for producing random directions is called several times in order to produce the directions necessary, it may (if the sample covariance matrix is dominated by only few large cigenvalues) produce essentially the same vectors over and over. We deal with this by considering two complementary cases. In the first case, where the random
'We note that [DKPP22] uses JJI. matrix of polylog(d/e) size instead of single direction. However, that is superficial difference, and random direction also works; see, for example [DK23; DKPP23J. 'The brute- force- approach takes 20(k) time. Hence, & is chosen to be (log((1/)) to satisfy the filter requirement while ensuring the brute-force approach runs in polynomial time.
3

---

directions selected by the technique of [DKPP22] have small inner product, we show that we can simply take & = (log(d/e)) many random directions and perform the multi-directional filter of [DKKLMS] 18] in these directions. When this is not the case, the covariance matrix must have a small number of dominant eigenvectors. We take one of these eigenvectors and put aside (for now) and instead consider the projections of our original points on the orthogonal subspace. Via careful analysis, it can be shown that either of these steps will lead to a substantial drop in the potential function of [DKPP22]. This guarantees that after small number of filter steps, we are left with matrix with small eigenvalues, allowing us to use the sample mean to approximate the true mean, at least in the directions orthogonal to those put aside.
In the directions that have been put aside, we still need to compute the sample mean. Fortunately, this reduces to the case where the dimension is polylog(d/e). Although a brute-force approach will not run in polynomial time for such a large subspace, careful analysis of [DKKL 18] can be made to yield an appropriate runtime, because the "effective" dimension is now only polylog(d/e).
We now sketch our algorithm for robust linear regression with Gaussian covariates. Our algorithm for robust linear regression in the Huber model achieves near-optimal error of ((aa) by leveraging our above-described robust mean estimation algorithm. Prior our work, no polynomial time algorithm was known to achieve this optimal error guarantee. As an additional bonus, our algorithm has near near-optimal sample complexity and runs in almost- linear time. The basic idea here is to consider the distribution of X conditioned on the value of y (or more precisely, after rejection sampling based on y in a carefully chosen way). We show that this conditional distribution will be nearly-spherical Gaussian whose mean is proportional to B, the true regressor. By using our robust mean estimation algorithm on this conditional distribution, we can obtain our final estimate of f.
1.3 Related Work
Our work lies in the field of algorithmic robust statistics, an active area of research since [DKKLMS16; LRV16]; see [DK23] for a recent book on this topic and Appendix A for additional discussion.
Contamination Models The strong contamination model works as follows: a computationally unbounded adversary inspects all the samples and can replace e-fraction of samples with points of their choice [DKKLMS16; DK19]. In contrast, the Huber model is additive and preserves independence among data points. In the rest of this paragraph, we focus on the task of robust (isotropic) Gaussian mean estimation. Despite these differences between the contamination models, the ciemmation-thhorrtti error in both of these models is (e)). However, the computational landscape of the problem changes considerably. First, [DKKLMS18] was able to achieve this error with polynomial samples and time. However, under the strong contamination model, there is evidence in terms of statistical query lower bounds (and low-degree hardness using [BBHLT20] that all poly- nomial time algorithms must incur a larger error of Oi /log/1/////// /c)) [DKS17], matching the existing polynomial time algorithms [DKKLMS16], which necessitates that we consider Huber's model.5
Gaussianity and Identity Covariance assumptions First, even for Gaussians, knowing the co- variance matrix is crucial for computationally--ffiiiee algorithms." [DKKLMS18] that gets O(e) error runs in polynomial time for isotropic covariances, while the algorithm for unknown covariance matrix runs in quasi-polynomial time [DKKLMS18)." Thus, we focus our attention to (nearly) isotropic distributions in this work and improve the runtime of the algorithm from DKKLMS18) from a large polynomial to almost-linear. Looking beyond Gaussianity assumption, -itmationtthoorttccllyy to obtain Oee) error for arbitrary isotropic subgaussian distributions, even in single dimension and Huber contamination model [DK23]. Finally, our results can be extended to subset of symmetric isotropic subgaussian distributions; see Remark 2.2.
impossible
'The lower bound from [DKS17] continues to hold for total variation corruption model, which is stronger than the Huber model but weaker than the strong corruption model. 'For examplc, all polynomial-tin
(statistical query and low-degree testing) algorithms for Gaussian mean estimation must use (2(x²) samples get o(v2) error [DKS19; DKKPP22] as opposed to O(c) error achievable with (d) samples using a computationally inntioonnall---fffiiiinnt inefficient algorithm.
1IDKKMMS6 16] gets larger error D(e/log) e 1) for unknown covariance in polynomial time and samples.

---

Outlier-robust Estimators in Nearly-Linear Time Several recent works have developed nearly- linear time algorithmss for problems in robust statistics: mean estimation [CDG19; DL22; D9L19; LLV720; DKKLT21; DKPP22; CMY20), linear regression [CATJFB20), and PCA [JLT200 DKPP23].
2 Preliminaries
Notation We denote [n.:- ..). For w R4 1] and a distribution P. we use Pw to denote the weighted by version of P. i.e., the distribution with pdf (w.( w(x) Ex~p(w( (x)] We use μip, Ep for the mean and covariance of P. When the vector ja is clear from the context, we use Ep to denote the second moment matrix of P centered with respect to μa, i.e., pp == Ex~Pl(X (X ". We use I1 for l norm of vectors and tr(-) and || llop for trace and operator norm of matrices, respectively. For a subspace V. we use Hl to denote the orthogonal projection matrix that projects to V. use Projy(1) (1:) to denote the orthogonal projection of x onto V, and use V for the sub- space orthogonal to V. We use x S y to denote that x S Cy for some absolute constant C. We use the notation > b to mean that > Cb where C' is some sufficiently large constant. We use polylog() to denote quantity that is poly-logarithmic in its arguments and use O S2() and θ to hide such factors.
2.1 Goodness Condition
Our algorithm builds on the following (slightly modified) goodness condition from [DKKLMS18]
Definition 2.1 (e, a, Gooddesss,), We say that a distribution G on Rd is (e, a, k) good with respect to μ. E Rd,, the following conditions are satisfied:
(1) (Median) For all v E Sd-1, Prx~allo" (x- -p) <
< 1/2.
(2) For every weight function with Bxxallw(X)) 2 1 the following hold:
(2.a) (Mean) lING -Mll S a (lo((/a),
(2.b) (Covariance) |E - Iallop S a (1g(()/))
(2.c) (Concentration along nearly-orthogone and [|U Ullop < 2k/ (og) the degree-2 polynomial p(x) == |U(x - 11)1 satisfies Ex-o.. (p(X)1 (p(x) > 100k)] S e/log(1/e)
vectors) For any U E R with tr(U k
The first condition, (1), states that the median is good estimate in each direction. The next two conditions, (2.a) and (2.b), state that the mean and covariance of the inliers change by at most a lo/looo//// and a log(1/a). respectively, when a-fraction of the dataset is deleted; these two conditions with a = e have been extensively used in the strong contamination model for (sub)- Gaussian distributions but the resulting algorithms necessarily get stuck at error c/log error. Condition (2.c), is crucial in obtaining O(c) error because it provides stronger concentration for projections along k-dimensional (nearly-orthogonal) projections using Hanson-Wright inequality. In particular, it implies that at most O O (oog((1/ -fraction of inliers have p(x) >> 10k. In contrast, using only Conditions (2.a) and (2.b) would require the threshold be much larger, p(x) 2 k llg(k//)) which is insufficient for our application. We show in Appendix B that if G is set of n >> (d + (ogl (1)5) polylog(d/c)/? i.i.d. samples from (P, I), then satisfies the goodness condition mentioned above with probability 1 - 5; Prior work [DKKLMS18 8] had only shown poly (dde) bound for the sample complexity.
Remark 2.2. Our proof for the sample complexity directly extends to all distributions satisfying the following: (i) their linear projections are subgaussian and centrally symmetric with constant density around the median, and (ii) their quadratic projections satisfy Hanson Hanson-Wright inequality. In particular, this extension is crucial for deriving our results for robust linear regression.
Throughout the algorithm, we will assume that the inliers and the outliers satisfy the following:
Setting 2.3. Let E (0, <o) for a sufficiently small constant 6p, C' >> be a large enough constant, P be a uniform distribution on n points of Rd that has the mixture form (1 - - + eB for some G that is (c, a. k)-good (cf. Definition 2.1) with respect to p E Rd, where a c/ log(1/c) and k (log 2 (d/e)). Let w Rd > [0, 1] be a weight function with Ex~a [w(X)) 2 - a such that -Cd Ep. I Co 2cloo2(1/c)) (I/C)I.
The condition on Σp. in Setting 2.3 amounts to "warm start" for our algorithm, which can be obtained in cearly-lincar time by running the algorithm from [DKPP22] first (see Appendix B).
5

---

The following result gives a sufficient condition for the (weighted) sample mean using (2.a) and (2.b):
Lemma 2.4 (Certificate Lemma). In Setting 2.3, if Er 1lop S + X, then I|p,_ =l X vd.
In light of the certificate lemma, our goal will be to downweigh outliers until the top eigenvalue (in the subspace of interest) is most + O(e). That would make the bound above O(e).
We now explain the filtering procedure that we use. A filtering algorithm takes weights w(x) for each point and scores ((x) (that capture how much of an outlier the point is) the procedure updates the weights to wxx) such that: (i) it removes much more mass from outliers than inliers, and (ii) Ex~Bl [w (x)] (the contribution to the weighted average of scores by the outliers) is small after filtering. By considering appropriate scores T(x)) that use multi nurti-directional projections of the form of Condition (2.c) of Definition 2.1, we show the following in Appendix B.
Lemma 2.5 (Multi-Directional Filtering). Consider Setting 2.3. Given nearly-orthogonal matrix U E IRixd satisfying ||U Ullop x 2 tr U)/ log(1/e), there is an algorithm that reads e, the n points, their weights w(x), and returns weights w' in time ndk + polylog(d/e, |[ such that
(i) Ex~a|w (X) - w(() < (e/ 1og(1/) ()) Ex-alwi (xx (((X))((()))) wx))).
(ii) Exx=Blu (xx=B[w((X]]]U)
uP_) ) S Sttr(UTU)
2.2 Polynomial Time Algorithm
We now record (refined) guarantee of the algorithm from [DKKLMS18), which would be useful later on. This algorithm uses the multi-directional filter from Lemma B.27 with U set to be the top-k eigenvectors of the covariance matrix until the top-k eigenvalues are + O(c) for k ((log(1/)) Lemma 2.4 then guarantees that the empirical mean in the subspace of the (d - k) last eigenvalues is Ole)-e to p. Finally, the algorithm employs a brute force procedure that uses median and (Condition (1)) to estimate the mean in the remaining k-dimensional subspace.
Lemma 2.6 (Adapted From [DKKLMS18), Let T be a set of n i..dd samples from an e-corrupted version of (0, La) and let 8 E 1) be parameters. If (d + (og(1/8))ffe )olyyogg(d// /ool///) then there is algorithm that, when having as input T, 6, and r, after time Ond + (xF ), it outputs vector çi E Rd such that 1| pll2 Z c/r with probability 8.
Since the refined runtime above does not appear in [DKKLMS18] we provide proof in Appendix B. In particular, it is important that the runtime is nd2 + zir. Even though this is super-linear in the size of the input, nd, our main algorithm in the next section will apply Lemma 2.6 only after projecting the entire dataset to a subspace a much smaller dimension (polylog(d/e) in place of d before).
3 Almost-Linear Time Algorithm for Robust Mean Estimation
In this section, we describe the main ideas of the almost-linear time algorithm for Gaussian robust mean estimation that achieves Theorem 1.3.
As described earlier, our algorithm runs in two stages, both of which will run in almost-linear time. In the first stage, our goal is to find a low- -ow-diensional subspace V (of dimension at most polylog( (d/c)) and a vector that is (c) close to μ when projected in subspace y- In the second stage, we deploy the basic O(c) algorithm from 1DKKLMS188 MS188 on the input data after projecting it onto V "'and estimate μ in that subspace. By a refined analysis of their argument (Lemma 2.6), the second stage runs fast because dim(V) = (polyyoog( (d/e)) e)). By combining these two estimates we get an approximation of p in the whole Rd Since the second stage algorithm's 's analysis closely follows [DKKLMS18] in the reminder we focus on the first stage, below: (the probability of success can be boosted by standard arguments)
Theorem 3.1 (First Stage in Almost-Linear Time). Given a set of n samples, the uniform distribution on which has the form (1 - G + eB for some G satisfying Conditions (2.a) to (2.c) of Desinition 2.1 with respect to p E TRd and parameters = e/ log(1 /c) and k = (logg (n + d)). Algorithm takes as input the n points and e, and with probability 0.99, returns a vector ç and a subspace V such that The algorithm runs in time O(ddpolylog (d/)).
[|Projy= (īi 2) S and dim(V)=polylog(dd
Theorem 3.1 is realized by Algorithm 1. Without loss of generality, we assume Setting 2.3 holds at the beginning of the algorithm. Recall the high-level strategy that was outlined in Section 1.2: We want to iteratively downweigh points and add directions to V so that (i) the weight removed

---

from inliers is at most (oee/og(///)) (ii) the downweighted dataset along every direction in V has variance at most + O(e), (iii) and dim(V) S polylog(d/c). Having the first two, the certificate lemma (Lemma 2.4) implies that the empirical (weighted) mean is ( to ş along V The way to achieve (i) and (ii) in [DKKLMS 18] was by using the matrix U to be the top-k eigenvectors of the covariance matrix (also, the algorithm of [DKKLMS18] does not add directions to V until the very end; see Section 2.2)). As mentioned earlier, this approach runs in quadratic time. Our main technical insight to (a) randomize the choice of when safe to do so, and (b) when it is not safe, allow the algorithm to remove a direction by adding it to V/ at any time (as opposed to waiting until the end). We first describe the notation below for a more detailed overview.
Notation for Algorithm 1. For each round E [mmax] of our algorithm, we maintain weights Wt : IR + [0, 1] over the dataset, capturing the confidence the points being inliers, i.c. w(x) = represents outliers and wxx) = represents inliers. We also maintain a (low-dimensional) subspace Vi with the goal of making the covariance of the projections of the data on V/- be small the end. Let M E be the sample mean and covariance of the data after projected on V and weighted by wi, and define B, = Σ Hy (see Lines to 7 and 4 for precise definitions). We use the potential function
6 tr (r( (M. |M, ll to track the progress of our algorithm, where M, (B,') ) for p = log d. Observe that the potential function Φi ignores the contribution from the directions in V4.
Algorithm Robust Mean Estimation Under Huber Contamination (Stage 1)
Input: Parameter 10,1/2), uniform distribution over points that can be written as (1-e)G eB where G satisfies Definition 2.1 with appropriate parameters.
=
Output: An approximation the mean a subspace V and the orthogonal subspace V.
1: Let C be a sufficiently large constant, & C'log" (n + d), tmax = (log(d/c))
2: Initialize V1 < 0 and wh (ar) for all RRd.
3: for t = 1,.... tmax do
4: Let V4 be the subspace spanned by the vectors in V,, and Vi- be the perpendicular subspace.
5: Let P/ be the distribution /' re- ro-weighted by We. i.e., Pi wr(x) wx(x) (xxx][((X)]
6: Let 14: Σ, be the mean and covariance of Projy, (X) when P:
7:
Define B, = (Ex~ P [wxx)]]] (
(1 Cie) y where Iy is the orthogonal projec-
tion matrix for V, and M, == (B; ) for p log(d).
8:
Calculate N such that B |lop E [0.1/10] using power iteration. It x,, S Ce then return f4 and Vg,
b cf. Appendix B
9:
10:
Let g: == "rs.:"//(o.allMMMMM MM 2 > |M 2 2 1/12///3)
11:
Calculate an estimate Te such that 17: qel I tmox
12:
if q. S 1/kkkt then
D Case (cf. Section 3.1)
13:
for j E [K] do
14:
Vi M Fe. for Zt, (O,
15:
U, + [........ vr,,] i.e., the matrix with rows Uz for j E [k]. e, U.)
16:
W++1 < MULTI-DIRECTIONLLFI
D cf. Lemma B.27 D Case 2 (cf. Section 3.2) 2||2 for p' = C. log (dbnmxx)) (O.
17:
else
Ut < (B 2///2/10 e/I 3'2l2
& Power iteration
19:
V + 1/2 )
20:
Let μl, = Ex~p, [Projy, (rojj,(XX) (X) be the mean of P. after projection Vi.
21:
return μl, and V,.
We will show that in every iteration, the potential function decreases multiplicatively, i.e., de+1 S (1- p-lylogge//)) while removing at most O(e log(1/ (<))) fraction of inliers throughout the algorithm (so that we do not fall outside of Setting 2.3). Since o; at f 0 is at most poly ((d/<))) and the algorithm necessarily terminates when reaches below e" (because this implics that | llop S (0)) O(c) which would cause Line 9 to activate), then after fmax = polylog(d/e iterations the algorithm yields HP., such that [|Proyy (HP.. - #)]l2 O(c) by Lemma 2.4. Since each iteration is implementable in ((d)) time, the whole algorithm terminates in Ó(dd) time.
7

---

In the next two subsections, we explain how the algorithm decides whether to expand the subspace V. or to remove outliers. This decision is based on Line 12, which checks if two random M, 2, M, 2 for 2, z' -N(0,)) I) are nearly-orthogonal with reasonable probability.
3.1 Case 1: Many Large Eigenvalues
By construction, "Case 1" corresponds to the case where the rows of U, are nearly-orthogonal (with high probability); see Line 12. Thus, U, will be oethoynorthogonal, permitting the use of the multi-directional filtering algorithm from Lemma B.27.
The explanation above ensures that the multi-directional filtering procedure with random v; '' is correct. The reason that it is significantly faster than [DKKLMS 18) 88 is that the v,'s are now randomized along the top eigenvalues of B7 because they are of the form M, ||| In contrast, [DKKLMS1 sets v; 's deterministically equal to the top-k eigenvalues of B As outlined in the introduction, the random choice of e 's prevents the "adversary" to place the outliers in such way that the outliers in the orthogonal subspace are unaffected during filtering. In particular, filtering with the random 2y's 's reduces the contribution of outliers, not only along the exact top-k eigenvectors of B,", but also along all directions with variance comparable to the top eigenvector. We use the technical insights from [DKPP22] to show that the potential function decreases multiplicatively:
Claim 3.2. With high constant probability, for every round t that Line 12 succeeds, 01+1 S 0.996r. Moreover: E Ex~Gl (X)] 2 e/ (1/c throughout the algorithm's 's execution.
Proof Sketch. We first sketch how U, is valid for the check of Line 12 succeeds, with high probability, we have that for each t E [fmaax] rounds, the angles of every pair of rows of U, formed in line 15 have cosine at most 1/k2² Also, Hanson- Wright inequal- ity implies that with high probability, ||ve,s 112 X Ir(U, U.) log(k). Combining both of these with Gershgorin Discs Theorem and the choice of k. we have that IIU Uzllop I (U, U.)/ log(1/e). satisfying the requirements in Lemma B.27 and ensuring Ex~a(wen (X)] 2 - e/ log()/e).
filter, Lemma B.27. Since the
To show that Φi+1 reduces, we first use the following linear-algebraic result to relate Φ++1 with from [DKPP22]: Φ1+1 tr((M,)))) < d (tr(M, 1 M,)) the definition of B,", we have tr (M, B B - = Ep.,. [[ (z 11 - (r((M,)))) To upper bound this, we will use the guarantees of the nulti-dimeesional filter along with the goodness condi- sketch of MA with & ž log(n); see Appendix B) as follows: we show in the appendix that the contribution from inliers is small by Definition 2.1 and the contribution from outliers, EB, [, 1 is controlled by Lemma B.27. Combining the two aforementioned arguments with some algebraic manipulations we can formally show that:
tr(M,
क By
tions and the fact that | (= 10]l x [|U,(zz 10)113 (since U, is
Lemma 3.3 (Filtering Implication). After the filtering, tr(M, B MF) < 0.1 (lpptr(M,,)))))
Combining Lemma C.4 with Φt+1 S dd/2pp (tr(M, B žाा from earlier, and noting that /2/2p- d/2llogd < 3 and | [2 < (r((M+))) ti, we have 6.+1 S 0.99 desired.
3.2 Case 2: A Few Large Eigenvalues
Consider now the alternate case where 9: from Line 10 is large. Then the approach from the previous case is inapplicable because Mf will be highly correlated vectors. The following result formally proves that this happens only when the top eigenvalue of M, contributes significantly to its spectrum.
Claim 3.4. If Pr (0)
[M M, z')l M =|a M,
7] 2 a, then
J|M Il ((((F)
2 poly(a,7))
Proof Sketch. Let ge be the probability from the statement. Pretend for simplicity that | =13 is equal to its expectation tr( Applying Markov's 's inequality, we see that a < 9t S "
--2 El B[M+ 2/ that || | < tr(A)]|A| lop for a psd matrix A applied to M,- completes the proof.
Further using the standard fact
L
Since y are polylog(d/é (d/c) in our setting (Line 12), Claim 3.4 reveals that the top eigenvalue of (M,-) )2 is at least (dde). Thus, our algorithm can take a top eigenvector of M; and add it
8

---

to the subspace of directions to ignore, V4, i.e., project all data on the subspace perpendicular to u in the future iterations. Formally, we show that the potential decreases if u has large projection along M,:
Claim 3.5. Let V4+1 Vi U {u} for some unil vector Vi², then 6e41 S Φh. u (M² Fuu.
Proof Sketch. Going from step & to step (t + 1), the effect of adding u in the subspace V is that B = AB, A for the projection matrix A = - uu" Then by properties of the trace operator and Thirrigg inequality:
de+1 = tr( (B )PP) = (r A )2P) = tr((B, A )P))
tr( (B,) ) A 27) (B( )3 A) = (r( )2)) - u (B, ) = Φ u )
Claims 3.4 and 3.5 imply that if we use u equal to the top eigenvector of M then our proof will be completed. However, the algorithm cannot compute exact eigenvectors in almost-linear time. Fortunately, power iteration in Line 20 computes a fine enough approximation such that u' ) is a constant fraction of ||M lop. Thus, we obtain Φ1+1 S (1 - 1/polylog(d/e))) as desired.
4 Robust Linear Regression: Optimal Error In Almost-Linear Time
In this section, we focus our attention to linear regression in Huber contamination model and sketch the proof of Theorem 1.4. Recall that each inlier sample (X,)) y) is distributed according to the Gaussian linear regression model of Definition 1.2.
Algorithm 2 Robust Linear Regression Under Huber Contamination
Input: >>0 O, multiset S of Rd+1 containing pairs of the form (z, v) with E Rd, E R.
Output: A vector in Rd
1: Set y <- TRIMMEDMEAN( (y/ (1r, ES))
D o ( O(e log(1/e)))
2: Draw a uniformly at random from I-āy. Fyl-
3: Define the interval I := [a - l,a + f] for E: = 0,/ ///(l//(/)))
> Random choice of I
4: S' < {æ (z. E S. y E r}.
D Simulating G₁
5: B1 < ROBUSTMEAN("", 10c)
b Algorithm from Theorem 1.3
6: return 3 (03/0)
D Rescaling
Our main technical insight is a novel reduction to robust mean estimation. Existing reductions in the literature (see, c.g., [BDLS17] I7]) rely on the fact that Ely E[xX" 8+Xx]] = 3. However, it un- clear how estimate the mean of y.X up to error (c)) under Huber contamination because y.X is very far from Gaussian; e.g. y.X is subexponential as opposed to SuussGGuu and it is not even symmet- ric around 88, implying median might not even work along a direction. Thus, existing approaches using this methodology have error Mo log(1/ 1c)).. In contrast, our reduction reduces to robust (almost)- Gaussian mean estimation by using the conditional distribution of X given y = a under Definition 1.2.
Claim 4.1 (Conditional Distribution). Let (X,v) E Rd+1 follow Definition 1.2. Given a E TR, denote by Ga the distribution of X given y = a. Similarly, given an interval I C R. let G₁ represent the conditional distribution of X given y I. Define oz =² - -118112. Then, Go N(23,L4 11 23. (38))) and G₁ = PrivET S1 pa 0, 0, N2) (533, Ia 1 oy 333 where Φ(2; 0, v²)
3,
denotes the pdf of the N(0, at 2 R.
Since the conditional distribution G. above is just Gaussian whose mean is scaled version of B (and roughly isotropic covariance), one would idcally like to get (oorucco)) samples from Ga. then use the robust mean estimator from Theorem 1.3, and finally scale the result back appropriately.
Obviously,
impossible to simulate G2 using G since there is zero probability over the inliers that y is exactly equal to a in our dataset. Instead, we can simulate the conditional distribution given y E I for some interval I centered around a. The length of I needs to be carefully selected. On the one hand, I should be sufficiently narrow to ensure that the mean of G1 closely approximates the mean of G.. On the other hand, it should not be excessively narrow to avoid rejecting a significant portion of our samples. As we show later, an interval of length (ar/ log(1/e)) meets both of these criteria as long as ||:1|| S are log(1/e). Also, set of i.i.d. samples from G, for such a small interval satisfy the goodness condition with the correct parameters because G2 will be roughly isotropic.

---

Finally, it is critical to ensure that the fraction of outliers in the simulated samples from G. remains O(e). In fact, this may not be true depending on the choice of the interval I. For example, the adversary knows T, then outliers could all happen to have their labels inside 7, which would cause the contamination rate to blow up. To overcome this, we randomize the selection of I: Since the distributions of outliers is independent of 7, choosing the center of 7 to be random means that the bad event the previous sentence is now a small probability event.
4.1 Proof Sketch of Theorem 1.4
The algorithm is given in Algorithm 2. We now give the proof sketch of Theorem 1.4: We claim that the following three events hold simultaneously with probability at least 0.9 (8o, 31 and og are defined in Algorithm 2):
(i 1
s logg(((11)))
(ii) |a| 0.00016,
(iii) || pa, 12 S
(I)
Given the above, we first show that 3 -coe--loss to B. For simplicity, let us assume momentarily
of that 0g = oz Then by triangle inequality || - B|| I Ja] (113, Ha, + l|AG, -/ (112) x oye + oylloG, a
of A1|2- The second term arises from the fact that we use Gy instead of G. in the algorithm (since simulating samples from Ga is algorithmically impossible). Using the expression for Gy from Claim 4.1, we can upper bound this term by og Fyllec, 3 S Al|llla Since o,
|| A oc log(1/c) and f := / log(1/c), this expression is overall ((ac).
Now, we show why the events in (1) above hold. The first is simply the guarantee of the trimmed mean on subexponential distribution (cf. Appendix B) and the second holds because a ~ u([-5,,, G,y)).
The third is significantly more involved. We first claim that Line 4 approximately preserves (within log factors) the size of the dataset and maintains the ratio of inliers to outliers (within constants). Recall that we use G for the inlier distribution and B for the distribution of outliers.
Lemma 4.2. Consider the context of Theorem 1.4 and Algorithmm 2. First, for every possible choice of I that can be made in Line 3, pr y E 1] Z f/o.g Second, E₁ (x) E 1 1[4]] S Elayy where the outer expectation is taken with respect to the random choice of the center of I.
Proof Sketch. For inlicrs, y ~ (0,, o')) and thus the fraction of inliers in S that belong to S" is at least (l(/a,,) since I C 2- 20,,] with length f. For outliers, we may assume without loss of generality that y E [-a,- -t, +E] always (since otherwise y & I). Since / is independent of (X, v). we may treat (X, y) as fixed and only / as random. Thus, the probability of y E is at most the ratio of the length of I to the length of [-oy- E, oy f), which is at most O((/oy) as claimed.
O
Observe that the set S" = ((x, (x. y) E S,y E 1} in Line consists of i.i.d. points from (1 Gr + B,, where G, is the conditional distribution from Claim 4.1, By is the conditional
distribution of the outliers, and e = (I-e) Prcx).)--- Prxx)// Pr(x. [ S (1-e) P
<Prcx wiwllee))
<Pr(x.g)=mlvel
By Lemma 4.2, Ey = (e) and thus by Markov's inequality, with high constant probability we will have that e O(e). Finally, it remains to show that S' satisfies the goodness conditions as
dictated by Theorem 1.3. This requires technical effort and is deferred to Appendix D.2.
5 Discussion
In this paper, we provided fast algorithms for mean estimation and linear regression that achieve optimal eттor under Huber contamination model for isotropic Gaussian inlier distributions. Several open problems and avenues for improvement remain. First, the sample complexity of our linear regression algorithm is multiplicative in log(1/6), where 8 is the failure probability, as opposed sample complexity [CGR16]. More broadly, it is an
to additive in the iifomnation-theoree open problem to design algorithms that achieve similar optimal guarantees for other fundamental tasks: robust principal component analysis, sparse mean estimation, and covariance estimation. In particular, the best known algorithm for the covariance estimation (or mean estimation with unknown covariance) achieving the optimal error runs in quasi-polynomial time [DKKLMS18].
10

---

Acknowledgments
Ilias Diakonikolas is supported by NSF Medium Award CCF-2107079, NSF Award (CAREER), and DARPA Learning with Less Labels (LwLLL) grant. Daniel Kane is supported by NSF Medium Award CCF-2107547 and NSF Award CCF-1553288 (CAREER). Ankit Pensia was supported by NSF Awards CCF- 165286,,, and CCF-1841 190, and CCF-201 1255; The majority of this work was done while Ankit Pensia was at UW-Madison. Thanasis Pittas is supported by NSF Medium Award CCF-2107079 and NSF Award DMS-2023239 (TRIPODS
References
[ABHHRT72]
D. F. Andrews. P. J. Bickel, F. R. Hampel, P. J. Huber, W. H. Rogers, and J. W. Tukey. Robust Estimates of Location: Survey and Advances. Princeton, NJ, USA: Princeton University Press, 1972.
[BBHLT20]
M. Brennan, G. Bresler, S. B. Hopkins, J. Li, and T.Schramm. "Statistical Query Al- gorithms and Low-Degree Tests Are Almost Equivalent". In: CoRR abs/2009.06/1 (2020).
[BDLS17]
S. Balakrishnan, S. Du, J. Li, and Singh. "Computationally Efficient Robust Sparse Estimation in High Dimensions". In: Proc. 30th Annual Conference on Learning Theory (COLT). Vol. 65. 2017, pp. 1-44.
[Ber18]
D. Bernstein. Scalar, Vector, and Matrix Mathematics: Theory, Facts, and Formu-
las. Revised and expanded edition. Princeton: Princeton University Press, 2018.
[Bha13]
R. Bhatia. Matrix analysis. Vol. 169. Springer Science & Business Media, 2013.
[BLM13]
S. Boucheron, G. Lugosi, and P. Massart. Concentration Inequalities: Nonasymp-
totic Theory of Independence. Oxford University Press, 2013.
[BP21]
A. Bakshi and A. Prasad. "Robust Linear Regression: Optimal Rates in Polynomial Time". In: Proc. 53rd Anmual ACM Symposium on Theory of Computing (STOC).
ACM Press, 2021, pp. 102-115.
[CATJFB20]
Y. Cherapanamjeri, E. Aras, N. Tripuraneni, M. Jordan, N. Flammarion, and P. L. Bartlett. "Optimal Robust Linear Regression in Nearly Linear Time". In: abs/2007.08 (2020).
[CDG19]
Y. Cheng. Diakonikolas, and R. Ge. "High-Dimmenional Robust Mean Estimation in Nearly-Linean Time". In: Proc. 30th Annual Symposium on Discrete Algorithms (SODA). SIAM, 2019, pp. 2755-2771.
[CDGS20]
Y. Cheng, Diakonikolas, R. Ge, and M. Soltanolkotabi. "High-Dimensie Robust Mean Estimation via Gradient Descent". In: Proc. 37th International Conference on Machine Learning (ICML). 2020.
[CDGW19]
Y. Cheng, I. Diakonikolas, R. Ge, and D. P. Woodruff. "Faster Algorithms for High- Dimensional Robust Covariance Estimation". In: Proc. 32nd Annual Conference on Learning Theory (COLT). 2019,
[CDKGGS22]
Y. Cheng, I. Diakonikolas, D. M. Kane, R. Ge. S. Gupta, and M. Soltanolkotabi. "Outlier-Rouust Sparse Estimation via Non-Convex Optimization". In: Advances in Neural Information Processing Systems 35 (NeurIPS). 2022.
[CFB19]
Y. Cherapanamjeri, N. Flammarion, and P. L. Bartlett. "Fast Mean Estimation with Sub-Gaussian Rates". In: Proc. 32nd Anmual Conference on Learning Theory (COLT). 2019.
[CGRI6]
M. Chen, C. Gao, and Z. Ren. "A General Decision Theory for Huber's $VepsilonS- Contamination Model". In: Electronic Journal of Statistics 10.2 (2016), pp. 3752- 3774.
[CGR18]
M. Chen, C. Gao, and Z. Ren. "Robust Covariance and Scatter Matrix Estima- tion under Huber's Contamination Model". In: The Annals of Statistics 46 (2018), pp. 1932-1960.
[CMY20]
Y. Cherapanamjeri, S. Mohanty, and M. Yau. "List decodable mean estimation in nearly linear time". Proc. 61 IEFE Symposium on Foundations of Computer Science (FOCS). 2020.

---

[CTBJ22]
Y. Cherapanamjeri, N. Tripuraneni, P. L. Bartlett, and M. I. Jordan. "Optimal Mean Estimation without Variance". In: Proc. 35th Annual Conference on Learning
Theory (COLT). 2022.
[Dep20]
Depersin. "A Spectral Algorithm for Robust Regression with Subgaussian Rates". In: CoRR abs/2007. 0bss/2007.06072 (2020).
[DHL19]
Y. Dong. S. B. Hopkins, and J. Li. "Quantum Entropy Scoring for Fast Robust Mean Estimation and Improved Outlier Detection". In: Advances in Neural Information Processing Systems 32 (NeurIPS), 2019.
[DK19]
I. Diakonikolas and D. M. Kane. "Recent Advances in Algorithmic High- Dimensional Robust Statistics". In: CoRR abs/1911.05911 (2019).
[DK23]
Diakonikolas and D. M. Kane. Algorithmic High- High-DDiensional Robust Statistics. https h //sites ceaites.googg google. www/view/ ars book/. Cambridge University Press, 2023.
[DKKLMS16]
Diakonikolas, G. Kamath, D. M. Kane, J. Li, A. Moitra, and A. Stewart. "Robust Estimators in High Dimensions without the Computational Intractability". In: Proc. 57th IEEE Symposium on Foundations of Computer Science (FOCS). 2016, pp. 655- 664.
Diakonikolas, G. Kamath, D. M. Kane, J. Li, A. Moitra, and A. Stewart. "Be- ing Robust (in High Dimensions) Can Be Practical". In: Proceedings of the 34th International Conference on Machine Learning, ICMI. 2017. 2017, pp. 999-1008. [DKKLMS17]
[DKKLMS18]
I. Diakonikolas, G. Kamath, D. M. Kane, J. Li, A. Moitra, and A. Stewart. "Robustly Learning a Gaussian: Getting Optimal Error, Efficiently". In: Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Al Algorithmms SODA 2018. Full version available at https//arxiv..rr/aabe www..orr/abs/7704.3336603
2018, pp. 2683-2702.
[DKKLSS19]
Diakonikolas, G. Kamath, D. Kane, J. Li, J. Steinhardt, and A. Stewart. "Sever: A Robust Meta-Algorithm for Stochastic Optimization". In: Proceedings of the 36th International Conference on Machine Learning, ICML 2019. 2019, pp. 1596-1606.
[DKKLT21]
I. Diakonikolas, D. M. Kane, D. Kongsgaard, J. Li, and K. Tian. -List-decodable mean estimation in nearly-pca time". In: Advances in Neural Information Processing Systems 34 (NeurIPS). Vol. 34. 2021.
[DKKLT22]
Diakonikolas, D. M. Kane, D. Kongsgaard, J. Li, and K. Tian. "Clustering Mixture Models in Almost-Linest Time via List-Decodable Mean Estimation". In: Proc. 54th Annual ACM Symposium on Theory of Computing (STOC). 2022.
[DKKPP22]
Diakonikolas, D. M. Kane, S. Karmalkar, A. Pensia, and T. Pittas. "Robust Sparse Mean Estimation via Sum of Squares". In: Proc. 35th Annual Conference on Learning Theory (COLT). 2022.
[DKKPSI9]
Diakonikolas, D. M. Kane, S. Karmalkar, E. Price, and A. Stewart. "Outlier- Robust High-Dimensionn in Neural Information Processing Systems 32 (NeurIPS). 2019.
Sparse Estimation via Iterative Filtering". In: Advances
[DKLP22]
Diakonikolas, D. M. Kane, J. C. H. Lee, and A. Pensia. "Qutlier-obust Sparse Mean Estimation for Heavy-Tailed Distributions". In: Advances in Neural Informa- tion Processing Systems 35 (NeurIPS). 2022.
[DKP20]
I. Diakonikolas, D. M. Kane, and A. Pensia. "Outlier Robust Mean Estimation with Subgaussian Rates via Stability". In: Advances in Neural Information Processing Systems 33 (NeurIPS). 2020.
[DKPP22]
Diakonikolas, D. M. Kane, A. Pensia, and T. Pittas. "Streaming algorithms for high-dimensional robust statistics". In: International Conference on Machine Learning. PMLR. 2022, pp. 5061-5117.
[DKPP23]
Diakonikolas, D. M. Kane, A. Pensia, and T. Pittas. "Nearly-Linear Time and Streaming Algorithms for Outlier-Robust PCA". In: Proc. 40th International Con- ference on Machine Learning (ICML). 2023.
[DKS17]
Diakonikolas, D. M. Kane, and A. Stewart. "Statistical Query Lower Bounds for Robust Estimation of High-DDimnsiionl Gaussians and Gaussian Mixtures". In: 58th IEEE Annual Symposium on Foundations of Computer Science, FOCS 2017. Full version at http://arriv.or//ab///// 1.03473. 2017, 7, pp. 73 73-84.
12

---

[DKS19]
I. Diakonikolas, W. Kong, and A. Stewart. "Efficient Algorithms and Lower Bounds for Robust Linear Regression". In: Proc. 30th Anmual Symposium on Discrete Algorithms (SODA). 2019.
[DL22]
J. Depersin and G. Lecué. "Robust Sub-Gaussian Estimation of a Mean Vector in Nearly Linear Time". In: The Annals of Statisticy 50.1 (Feb. 2022), pp. 51 1-366.
[DM22]
A. S. Dalalyan and A. Minasyan. " "All-In-One Robust Estimator of the Gaussian Mcan". In: The Annals of Statistics 50 (2022).
[HLZ20]
S. B. Hopkins, J. Li, and F. Zhang. "Robust and Heavy-Tailed Mean Estimation Made Simple, via Regret Minimization". In: Advances in Neural Information Pro- cessing Systems 33 (NeurIPS). 2020.
[HR09]
P. J. Huber and E. M. Ronchetti. Robust Statistics. John Wiley & Sons, 2009.
[HR21]
L. Hu and O. Reingold. "Robust Mean Estimation on Highly Incomplete Data with Arbitrary Outliers". In: Proc. 24th International Conference on Artificial Intelligence and Statistics (AISTATS). 2021.
[Hub64]
Huber. "Robust Estimation of Location Parameter". In: The Annals of Mathe- matical Statistics 35,1 1 (Mar. 1964), pp. 73-101.
[JLT20]
A. Jambulapati, J. Li, and K. Tian. "Robust uub-gaussian principal compo- nent analysis and width-indeppnden schatten packing". In: Advances in Neu- ral Information Processing Systems 33 (NeurIPS) (2020). arxiv preprint at https://arxiv http://wwwi..orrgaaa////.0.. 06980.
[KKM18]
A. Klivans, P. K. Kothari, and R. Meka. "Efficient Algorithms for Outlier-Robust Regression" In: Proc. 31st Annual Conference on Learning Theory (COLT). 2018. W. Kong, R. Somani, S. Kakade, and S. Oh. "Robust Meta-learning for Mixed Linear Regression with Small Batches".. In: Advances in Neural Information Processing Systems 33 (NeurIPS). 2020.
[KSKO20]
[KSS18]
P. K. Kothari, J. Steinhardt, and D. Steurer. "Robust Moment Estimation and Improved Clustering via Sum of Squares". In: Proc. 50th Annual ACM Symposium on Theory of Computing (STOC). ACM Press, 2018, pp. 1035-1046.
[LLV720]
Z. Lei, K. Luh, P. Venkat, and F. Zhang. "A Fast Spectral Algorithm for Mean Esti- mation with Sub-Gaussian Rates". In: Proc. 33rd Annual Conference on Learning Theory (COLT). 2020.
[LM21]
G. Lugosi and S. Mendelson. "Robust Multivariate Mean Estimation: The Optimality of Trimmed Mean". In: The Annals of Statistics 49.1 (2021), pp. 393-410.
[LRV16]
K. A. Lai, A. B. Rao, and S. Vempala. "Agnostic Estimation of Mean and Covari- ance". In: Proceedings of FOCS' 16. 2016.
[LY20]
Li and Ye. "Robust Gaussian Covariance Estimation in Nearly-Matrix Multipli- cation Time". In: Advances in Neural Information Processing Systems 33 (NeurIPS). 2020.
[MM15]
C. Musco and C. Musco. "Randomized block krylov methods for stronger and faster approximate singular value decomposition". In: Advances in neural information processing systems 28 (2015).
[PBR19]
A. Prasad, S. Balakrishnan, and P. Ravikumar. "A Unified Approach to Robust Mean Estimation".. In: CoRR abs/1907. 00927 (July 2019).
[PJL20]
A. Pensia, V. Jog, and P. Loh. "Robust Regression with Covariate Filtering: Heavy Tails and Adversarial Contamination". In: CoRR abs/2009. 12976 (Sept. 2020). J. W. Tukey. A survey of sampling from contaminated distributions". In: Contribu- tions to probability and statistics (1960), pp. 448-485
[Tuk60]
[Ver18]
R. Vershynin. High-Dimensional Probability: An Introduction with Applications in Data Science. Cambridge University Press, 2018
[Wai19]
M. J. Wainwright. High- Dinennimeesional Statistics: A Non-Asymptotic Viewpoint. Cam- bridge University Press, 2019.
[XCMI3]
H. Xu, C. Caramanis, and S. Mannor. "Quutier-Rbuust PCA: The High-Dimensional Case", In: IEEE Transactions on Information Theory 59.1 (2013), pp. 546-572. B. Zhu, J. Jiao, and J. Steinhardt. "Robust Estimation via Generalized Quasi- Gradients". In: Information and Inference: Journal of the IMA 1 11.2 (2022),
[ZJS22]
pp. 581-636.
13

---

Appendix
Table of Contents
A Additional Related Work
14
Additional Preliminaries
15
B.1 Linear Algebraic Facts
15
B.2
Concentration of Measure Facts
15
B.3
Existing Robust Algorithms: Trimmed Mean and [DKPP22]
16
B.4
Goodness Condition and Its Sample Complexity
17
B.5
Filtering Procedure
22
B.6 Algorithm from [DKKLMS18]
23
Omitted Proofs from Section
25
C.1
Deterministic Conditions for Algorithm
26
C.2 Case
26
C,3 Case
28
C.4 Combining Everything Together
31
Omitted Proofs from Section
32
D.I Conditional Distribution of Covariates
32
D.2 Robustness Goodness Conditions to Spectral Noise ................................
33
D.3 Proof of Theorem 1.4
34
Adaptation of Algorithm 1 to the Streaming Setting
36
A Additional Related Work
A systematic study of robust statistics initiated in the 1960s ([Tuk60; Hub64; ABHHRT72),. We refer the reader to the book [HR09] for classical treatment of this topic. The theoretical investigation in these studies often focused the asymptotics of sample size going infinity. A non-asymptotic analysis of the sample complexity can be found in, for example, [CGR18; CGR16; LM21].
Algorithmic high-dimensional robust statistics has been an active area of research since the publication of [DKKLMS 16; LRV16]. Various high-dimensional inference tasks have been investigated the literature: (i) Unsupervised learning tasks include mean estimation [DKKLMM16; LRV16; KSS18;
CFB19; DKP20; HLZ20; CTBJ22: HR21J, sparse estimation [BDLS17: DKKPP22: DKLP22; DKKPS19], covariance estimation [CDGW19, 19; LY20J, and PCA [XCM13; KSKO20; JLT20], and (ii) Supervised learning tasks include linear regression [DKS19; KKM18; CATJFB20; PJL20; Dep20; BP21], and general stochastic optimization tasks [PBR19; DKKLSS19].
Within algorithmic robust statisties, our work is related a line of work that develops novel algorithmic approaches to robust inference. For example, approaches based on gradient-descent and matrix multiplicative weights update [CDGS20; ZJS22; HL.Z20; CDKGGS22). several recent works have developed fast algorithms, often nearly linear time opposed to simply polynomial time algorithms, for robust inference [CDG19; CDGW19, LY20; ZJS22; CDKGGS22; DHL19; CMY20; DKKLT21; DKKLT22).
particular,
14

---

B Additional Preliminaries
In this section, we state the results that will be used the proofs later. Appendices B.1 and B.2 contain facts pertaining to linear algebra and concentration of measure, respectively. Appendix B.3 records some results from both folklore and prior works in robust statistics that are needed this paper. Appendix B.4 focuses on the goodness condition and proves the samples complexity and key properties of good sets. Appendix B.5 formally states the guarantee of the multi-directional filtering procedure. Finally, Appendix B.6 presents the improved runtime guarantee of the algorithm in [DKKLMS18].
B.1 Linear Algebraic Facts
Fact B.1. Let A be PSD matrix. Then for any unit vector and 2 1, x A 2 (x Ax)
Fact B.2. If A E RRddd is symmetric and p 2 1, the Schatten norms of A satisfy the following:
||A p & |l S || [1lp
In the following we let C denote the set of complex numbers, though throughout the paper we will only work with real matrices.
Definition B.3 (Gershgorin Discs). For any complex n x n matrix A, for i E [n]] let (A = Eist (aegl and let GAA) U C aiil I [.)A Each disc {z E C: 12- and S
(A)) is called Gershgorin disc and their union G(A) is called the Gershgorin domain.
Fact B.4 (Gershgorin's Disc Theorem). For any complex n malrix A, all the eigenvalues of A belong to the Gershgorin domain G(A).
Fact B.5. AfA, B, C are symmetric d x d matrices with A 0 and B < C then tr(AB) S ((((AC).
The following two linear algebra results follow from the Lieb- Shib-nhiring inequality.
Fact B.6 (Lemma in [JLT20]). Let A, be positive semidefinite matrices with 3 and p E N. Then, tr (BP) S ir (Ap-1B).
Fact B.7 -Lie-Thhirinn Inequality; see, e.g., Problem IX.8.1 1 in [Bhal3]). Let A, be positive semidefinite matrices. For all & E N. tr((A))") S (B*)
Fact B.8. Let A be a PSD matrix. Then, ||||| S V A op Ali- In particular, Al I ||A]I1/1
then | < |///t
Fact B.9 (See, c.g., [Ber18, Proposition 11.10.34]). Let R and M two square matrices, then [RM]| < min [[] [M]
||||| [[]r and [RM] |
B.2 Concentration of Measure Facts
Fact B.10 (Power iteration). For any positive semidefinite matrix A E dddd and n. E (0, 1), if p > C (1)(d/(16)) for sufficiently large constant C, and u == A for = ~ (0, then
Pr [u Au/ 1ul 2 (1- App 2 8
Definition B.11 (Sub-Gaussian Gaussian and Sub-exponential Random Variables), A one-dimensional random variable Y is nab-Gaussian if || := =upp21 1---- 1/2 E [ is finite. We say that ||Y|], the sub- Gaussiau norm of Y. A random vector X in Rd is SabbGaussian if for every E 1d u0 x11,, is finite. The nab-Gaussian norm of the vector defined to be
|X
112 = sup
|
[ 1/
DESa-1
We call a centered one-dimensional random variable Y a (v, x)+ sub-exponential Elexp(XY)) < x² for all X < 1/a. We call [|Y||, == supp21 p=1 E||YP]] the sub-exponential norm of Y. We extend our definition to vectors as before.
Lemma B.12 (Properties of Sub-exponential Random Variables [Wai19; BLM13].. The class of sub-exponential random variables satisfy the following:
15

---

1. [Wail9, Proposition 2.9] Y is centered (v, a))+ sub-exponential random variable. then with probability 1 6, Y $ {vVlog(L/0) + log((1/0).
2. [BLM13, Theorem 2.3] If is centered random variable satisfying that for all 5 E (0, 1),
Y S <v/log(L/0) + log(1/6), then Y is (v',a')= sub-exponential with v/ S v + and S a.
3. [Wai19, Section 2.1.3] Let YI...., Yg be k centered independent (v. a) sub-exponential random variables. Then Ex.1 Y; is (uVli, a))+ sub-exponential random variable.
Fact B.13 (Hanson- Haanson-Wright Inequality). Let X (N, Then, for every 2 0:
(-0.1
12
K3 K2)1A119))
Pr| E AX| > 1 S exp
m-i min
K+
Fact B.14 (Johnson-Lindenstrrn Lindenstrauss Sketch). Fix set of a points 21..... In E Rd. Let g(z) == |a: b)ll be a polynomial for some A Rdddd and b E IRd, and let U be the (random) matrix having as rows the vectors u, Az/VK for 2; Ia) [[(( - 12. If C is sufficiently large constant and & C'log((n + d)/8). then with probability least 5, the the following holds:
1....,k. Define g(x)
1. 0.8g(z;) & ((x,)= & 1.2g(zz)) for every i E [n].
2. A S |U|P S 1.2
Fact B.15. For any d x d symmetric matrix A, we have that Var I AAz, A=') || If A is PSD matrix, then for any 3 > it holds Pr A Str (A)] 21 VeB.
Fact B.16 (VC incquality). Let FF be a class of Boolean functions with finite VC dimension VC(F) and let a probability distribution D over the domain of these functions. For a set S of n independent samples from D
VC(F)
log( 1/T)
sup Pr Pr(f(X) Pr [f(X) 2 SEF
71
with probability at least T.
B.3 Existing Robust Algorithms: Trimmed Mean and [DKPP22]
In this section, we state the guarantees of the trimmed mean and the fast robust mean estimation algorithm from [DKPP22].
We begin with the guarantee of the trimmed mean. The following holds even in the strong contami- nation model where adversary edit c-fraction of the samples arbitrarily. We will need it to estimate the (unknown) variance of the labels in our linear regression algorithm.
Fact B.17 (Univariate Trimmed Mean, see, e.g., [LM21]. Let co be sufficiently small absolute constant. There is an algorithm (trimmed mean) that, for every E (0, co) and univariate distribu- tion D that has a) -oorm at most o² (cf. Definition B.11), given a set of n > log (og(1/8)/ (c//)//(ee (clog" (1/c)) samples from D with corruption rate e, outpuls g such that T Ex=D[X] S d²e log (1e) with probability at least 5.
The following result is implicit in [DKPP22]. It will serve as preprocessing step to ensure small operator norm of the covariance matrix in our robust mean estimation algorithm.
Lemma B.18. There exists algorithm for which the following is true. S is a set of n points in Rdd containing a C S such that the uniform distribution on G satisfies conditions (2.a) and (2.b) of Definition 2.1 with
e/ log(1/c), then the algorithm w(x) E [0, for each r E such that the following hold with probability at least 0.999:
input S and outputs weights
1 Exec w(x) 2 1 logli/e) 1. [o|
2. Let Q be the discrete distribution on that assigns probability mass uxxx)/ Eves w(x) each x E S. Then, the top eigenvalue of the covariance matrix of Q is most 1 O(elog((1/e).
16

---

Moreover, the algorithm runs in time O(ndpolylog(d/e))
Proof. The algorithm that achieves the guarantee is a slightly modified version of Algorithm 1 in [DKPP22] where: (i) instead of outputting an estimate of the mean, the algorithm outputs the underlying weights that it has been updating and (ii) the "Downweightnggiller procedure is replaced by Algorithm 3 with 3 = log(11 (1/c,, i.e., the only difference is on the threshold that is used to decide when to stop filtering. Algorithm increases the threshold by the factor 3 (1e). which as we show in Lemma B.27 results in removing log(1/e) times more mass from the outliers than inliers. The increase the filtering threshold translates to having an extra log(1 /e) factor Item 2 (without this increased threshold, the variance along every direction would have been log(1 (o)(1))),, and the fact that this filter now removes log(1 (1/e) times more mass from the outliers than inliers translates to having divided by log(1/e) in Item 1.
B.4 Goodness Condition and Its Sample Complexity
We first restate the definition of the conditions. Recall that for a weight function w Re [0, 1] let G denote weighted by version of the distribution G. The matrix Eg,, below denotes Ex.G... [XX (x [] (as opposed to the covariance matrix Eg. Fxx~o..(X" px (Ga)
Definition 2.1 ( (e, at k) -ooodessss)) We say that a distribution G on Rdd is (e, g)oo with respect to μ E Rd, the following conditions are satisfied:
(1) (Median) For all E S Prx-allo (x- - -1) ž c] 1/2
(2) For every weight function with Ex~alw(XX (X)] 2
a, the following hold:
(2.a) (Mean) lI -Mll S alog((/aa)
(2.b) (Covariance) || - -allop I log(L/a)).
(2.c) (Concentration along nearly-orthogonal vectors) For any U E Rxx with U)( and | Ullop < 2k/ log(k/e), the degree-2 polynomial p(x) | - m))l satisfies Ex~G... (p()))
> 100k)] S e/log(1/c).
Remark B.19. We note that once the conditions of Detinition 2.1 are established, then they also hold for projections of the data on a subspace, something that we will need later on in our proofs. Concretely, let V be any subspace of Rd. Then Conditions (1), (2.a) and (2.b) continue to hold for the projected points Projy(z).. This is because all of these conditions concern properties of one-dimensional projections of the data. Regarding the third condition, if the rows of U belong in the subspace V (which will be the case in our analysis later on), the last condition also continues to hold for the points Projy (a:) since the (>))rojp>>]]3 ||U 11 = ||aar|] 12 where used that UTly U because the rows of U already live in the space V.
B.4.1 Sample Complexity
In our results, use the goodness condition from Definition 2.1 with parameters = /og/oo////// and = polylog(d/e). The following result implies that a set of ((++lo((/)))ppologggl i.i.d. samples from the inliers satisfies the desired stability condition.
Lemma B.20. Let eo be small enough positive constant, let e E (0, E0), and let 7 E (0, 1). Let S C Rdd be of n. independent samples from (p, for >> ((dlog(d/o)+ min(a? e²/ [//lo"(1/e))
If G denotes the uniform distribution over S, then with probability at least - G is (c, a, k)ogo
Proof. We prove the results for each condition separately.
Condition (1) As the VC dimension of Linear Threshold Functions (LTFs) d, the VC inequality (Theorem 16) implies that there exists sufficiently large universal constant C' such that for any unit vector and 7, |Pxx~s[p' (X (4) > 7 Prx= Tll s C(Vd/n+ + (log(1/T)/m)
with probability 1 1- T.
We next apply Taylor's Theorem to the complementary cumulative function of the Gaussian, (tt)= Prx=Ni1,,)(X t] (L/V2m) 1e e 222 The first two derivatives are d'())
17

---

-E -E1////V2 and (t) = te /2- and thus there exists & between 0 and such that 112 - t/V2 +
(22 I 1/2 1/ 1 (2 (2v2T)) Combining with the VC
inequality, with probability 1 - т, we have that
Pr T M) > - 7/V
C(vd/n+ + Vlog(1/7)/n). (2)
By our assumption that v d/n + v right hand side is less than 1/2.
(()) and thus taking 7 = O(e), we obtain that the
Conditions (2.a) and (2.b) Let G be the uniform distribution over a set of n samples drawn from (),, I). We recall the following result from the recent robust statisties literature:
Lemma B.21 ([DK19; DKKLMS16)). Let (0) for small enough positive constant 60. Let P be an O(1) D(1)-subgasssian distribution with mean μ and covariance Σ, where [E-Illop S a log(1/o). The uniform distribution over a set of n > (d + log(1 /8)) points drawn i..d. from P satisfies Conditions (2.a) and (2.b) from Definition 2.1 with probability 5.
Condition (2.c) First we argue that (. > 100tr(U U)]] < e/ cogg(1/)) for all the polynomials p(æ) described Condition (2.c). Fix such a p(x) and observe that Elp(a)] = tr(UTUp T U). Under the Gaussian distribution, by Hanson-Wright inequality (Fact B.13), we have that
Pr [[x- U - 11) tr(U U) > 2exp (-0.1min X-NluI)
t²
Let V/ = lo)(k/ Using our assumption
Ulpp S U)/KK Fact B.8 implies that JUTIIÉ S tr(U U)2/K' Plugging in 99tr(U
(UTU) U) above, we obtain that
Pr X-N(ma)
[[x -4) U - n1) 2 100tr(U (U U)]] - 2exp( 48k) = 26 2k:24
(4)
where we use the definition of K. Moreover, integrating the tail inequality Equation (3) implies that (X)) & (E(p(X)))² | S tr(U"UU) Altemaatiely, Equation (3) implies that with
probability 6,
(u()) 2 (X)) s 2(pu (x) (p))()) E (X)))
S / /2/K Vlog(1//0/ + (k/k') log(1/6)
The covering argument is based on the following discretization of the unit sphere.
Fact B.22 (Cover of the sphere). Let r 0. There exists a set of unit vectors of IRd, such that |C| < (1 + 2/r)" and for every E Sd-1 we have that minyec |yy - v/l2 Sr.
Let V be the set of all nearly orthogonal matrices:
V (U Ridd
TUTUl S 2tr(U 2 UU
Let V., C // be the cover of V for some 1) (0, such that for any U V, have matrix U such that |U-U'll S v). We now bound the cardinality of --oover Vny- Looking it as a vector in dk as flattened vector, there exists an noover of size 2(3/n)dk where we also use that we can ensure V,
using [Ver18, Exercise 4.2.9).
For any U E V. define Po(x) | and TV( (p(())(pp(()) > 100tr(U U)). Let Xi,..., X,, be n i.i.d. random variables from N/ I). We first claim that for N D. (x(X)- E[r)(X)) is (O(k2/k') O(b///")-sueexxpeeeiia random variable and Elru(()) S e. For the latter, the Sauchy-Schwand incquality along with Equation (3) applies that
E E (X) = Elpu (X) 100tr(U U\) < (U))
Exri / (x) 2 100tr(U UU)]
18

---

which is less than e². For the former, we have that with probability 1 - 8.
(u(x) E[ro(X)] ( 2 Elpo(X))) & Vk2/w 1og(1/8) +
Thus, Lemma B. implies that TU (z) E is (O(k2/k) O(k/k)) subexponential random variable. Applying the concentration properties of the mean of independent subexponenial random variables (Lemma B. 12),,, we obtain that for any fixed U E Vi with probability 1 -
Σ L To(x S +
122 (og(1/r)
k llg((//T)
(5)
k'
k'n
Taking union bound over every U Vn, obtain the following: with probability 1
T,
k2 log(1/T)
A klog(1/T)
dk log(1/n)
dk²
log(1/n)
YU' e Vn= (z) e
+
+
k'n
kn
k'n
k'n
(6)
Let U E and U' E VV, be arbitrary such that | Ue I n. The goal now is show that E 1 (z,) and Σi=1 I (x;) are close to each other and thus the former bounded similarly n to (6). Define the difference of the polynomials Ap(x) == (2) (z). We will also use the notation [n| to denote the averaging over [n] samples. First, we have that
E [To( = E (Pu(()) > 100k)]
E (c;) (x;) > 100k - Ap(x:)) + E (pr(()) > 100k)]
i~[n]
1~[1
= E
> 100k Ap(x;)) + App())]]]
We start with the second term of (7).
Ap(x)]]
= E
IK U U' الااد|
i~ [n]
i~(n]
< -UU-UU
U
E | L~|n]
ddn
where we used two things: First that by Gaussian norm concentration (e.g., Wrightt inequality with A = combined with Lemma 1.1
log (og(1/6) va
log(1/0) = O(d),
Σ 1 < d
(8)
n
n
n log(1/0) and, second, we used that
JU- -TU-UUT U UT (U (U- U') UU (U- U` T
< 2 max (||U|lapp U'llap) ||U U'li S 2dn-
We now bound the first term in (7),
E in]
IPU
(pv)
> 100k Ap(x:))
E
(x:
( > 100k Ap(x;) $ s)]
E pu (p:) 100k Ap(x:)) > k)]
i~
E
Pu (1:)
994)] + E
Ip (Ap(a:) > A)]
(9)
dow ml
The first term is bounded by (
There is small difference in the constant in front of k. but the same proof gives similar quantitative bounds.
19

---

Choose 17 = 0.001c /d². Then second term in (9), (x.in]lpoov:::) > k)] is zero unless there exists an 2; in the sample set with ||= > 10vd. This relies on the fact that Ap(x) v ||=:|||| U|Ir I 2d1||>=11 which becomes less than k if ||0]| v (00dandd = 0.001 e/P.
But the probability of the event ||a|||| > 10vd is exponentially small:
Pr X1 N1
B : |X > 10v S ne ne-10000 4/100
(10)
by Gaussian norm concentration and a union bound.
Therefore taking 1) = 0.001 /de, we have that with high probability the scores from V and V, are close up to e Combining this with Equation (6) and the union bound over when n > d log(111)))) we obtain that that with probability least 1 - -6/2 / ne -d/100 for all U E V:
n n
k2 log(1/T) k'n +
A llg((1/))
13 log(//( dke log(d/c)
Tu(x)
Wn
Wn
k'n
for a sufficiently constant C, and c E (o) a sufficiently small constant cp, all of the terms above become less than c/ log(1//)) We also need ne-d/100 to be most 5/2. Using the same
Thus, as soon as n = 1 (ddddggdd)))
kee (1/e
before, the quantity -ddd/100 is at most 8/2 if >> poly 1og(e-18-1) We can assume without loss of generality that the last condition is always true can artificially augment the dimension by padding the data with Gaussian coordinates until the resulting dimension becomes d' > poly log(s (og(<-18-1) and it is easy to see that once the goodness conditions get established for the augmented data, they continue to hold for the original data.
B.4.2 Helper Lemmata Related to Goodness Condition
We now record implications of the goodness conditions that will use in analysis of the main algorithm. The first lemma below provides certification that when the operator norm of the empirical covariance is suitably bounded, the empirical mean closely approximates the true mean. This allows us to halt the outlier filtering process and return the resulting vector.
Lemma B.23 (Certificate Lemma (restated)). Let 0 < <<<<< 1/4 and 8 (0, P (1 c)G eB be mixture of distributions, where G satisfies Conditions (2.a) and (2.b) of Definition 2.1 with respect to μ E Rd. Let wxx) be such that Ex~c[w(X)] > 1 a. Assume that the the top eigenvalue of Ep. less than + Then,
llur, 1112 S (ov(l//)//) + Vx + v ar log(1/o)
Proof. Let P
c (Ex=B[w(X)]. (xxxXxxXXXxX/X
Denote
Pu(
w(x) P(x) Ex-p (owpp((XX (x)
w(x)B(x)/ Ex=n[w(X)],, (z))
w(x G(x) (x) Ex-a [w( x)] the weighted by w(x) versions of the distributions P, B, G and
denote by Σp.. E EB.. Eg. their covariance matrices. We can write:
Σp = PEB. + + - -D)EG. + + () (BG) - (BB.).
(1))
Let be the vector in the direction of μ.G.. PBB.. Since the top eigenvalue of Er. is less than 1 + X, we obtain the following:
1+ v
p)v EG., p(1 () (MB. HG.)))
21 (1-p)) (1 a log(L/a)) p(1 (v (HB.. - MG,))2
where the second step uses Equation (11) and the last step that G satisfies Conditions (2.a) and (2.b) of Definition 2.1. The expression above implies the following:
p a log(//))
(u (pB.. - pG_)))
p(1 p(1-p)
We can now bound the error [|rr.. - ul12 as follows:
|. nll2 lo [NG ()) p(uB.. PG.)
20

---

S 10 (MG. rl) + p/o (UB. - -HGW!
|| pll2 + p/c (UB. MGG..
A p a log(1/a) S (og//o/(///+++ + vPV
1 1-p
where the last inequality uses that G satisfies Conditions (2.a) and (2.b). We now use bounds on to simplify the terms. Recall that p = x (xx))//xxxxlpleeee/ (xxpl[((X) S e/(1 a). As < 1/2, we get that < 2e. In addition, note that p < 1/2
l1r. Hll2 S av log(/a/ + Vec+et Vaclog(1/a))
Recall the the goodness conditions are phrased in terms of the deviation from the true mean. However, the true mean is unknown the algorithm and the algorithm will use the empirical mean as an approximation, which introduces. additional errors in our bounds. To analyze these emrors, require the following lemmata.
Lemma B.24 (see, e.g., Lemma 2.13 in (DKPP22)). Let Rd - [0, 1] such that Ex~G(w Z 1 a and let G be distribution satisfying Conditions (2.a) and (2.b) of Definition 2.1 with respect to f IRd with parameters (c, a, k) 9 For any matrix U Rmxd and any vector E Rd, we have that
[||((X - b)] = | (1 log(//a) + |0( b)11 20 2ogeooloooo//
| bll2
X-G.
Lemma B.25. Let Rdd - [0, 1] such that Ex~G [w( (X) 1 a and let G be an (c, a, k)- good distribution with respect μ. E Rdd in Definition 2.1. For any matrix U Ikxxd, if we
define the polynomial T(z) := ||U( b) [I([U((r b)ll > 100tr(U U)) and assume that
(U(p-b]13 b)13 < u(UUU), U), then it holds
2c
E [x(X) <
Ir(U U) 21 b)ll
X XmG.
log(1/e
Proof. We have that
UU - b)112 < 2/ /011 21| b)1
2 --113 uU)
where the first step uses the triangle inequality combined with the inequality (a + b)2 < 2a + 218. Using the above, we can write
E S2 E [|((r - 10 ([ (1-0x=100) b]13 XNG
49Ir(U + 2 b)11
XNG.
S (e/log(1/e) 2 b)ll3
where the last line uses the goodness condition for G.
Lastly, the following lemma shows that the matrix U used in the last goodness condition consists of nearly orthogonal vectors, then the matrix satisfies the near-orthogonality condition mentioned there.
Lemma B.26. Let us for E [k] be vectors in Id and define U = [ui..... ux]. If every i E [k]] with i # j, it holds ((u, my)l S | Ullop 2tr(U U)/ (k/ log
and max; ||u, < log k (og"tr(UTU) then
Proof. We note that by the Gershgorin circle theorem (Fact B.4),
||U T Ullop max
||
12 + Σ
(u,
u,)l
ve
k
This lemma does the last goodness condition; thus the parameter does not appear the conclusion.
21

---

11 ka 1 @ulll|luu l
max kk
||u
k2
||u.112 - 2 U)
- max sE[R]
k/ log(k)
B.5 Filtering Procedure
In this section, we formally describe the filtering procedure (Algorithm 3) and prove its formal guarantee (Lemma 2.5).
Algorithm 3 Down-weighting Filter
Input: Distribution on points, weight w(x) for each point x, score (xx) for every point with the guarantee ((r) < r, threshold T > 0. parameters > 0, Liax E Z Output: New weights w/ (x)).
2:
3: Initialize (x(xx) w(x).
4: max f- r/(e()).
5: for I,..., lmax do
Can be implemented in Ollog max) with Binary Search
6: if Exxxplw' (xx T3 then
7: w w'(x) 6 (1) (z)/r))
8: else return w'.
9:
return w'.
Lemma B.27 (Filtering Guarantee). Ler P c)G + eB be mixture of distributions supported
on points and 3 > 1. If e) Ex Ex-clu(X)F(XX < T, || S r and Emax > r v/T then the new weights w (x) satisfy:
1. (1 - Ex-Exoolla w'(X)) < /ー Ex~B (X) - w(())
2. Ex-plu(X)F(X)) S
The algorithum can be implemented in O log(r///TT) time.
Proof. The weight removed from inliers is:
(1 - E) E [x( w w(X)] 1 - r
E [x())((X) S//r
X~G
The weight removed from outliers is:
E X~B
[w(X) w' = E [x()F((X)
X~B
[wo [x(X)F(X)) - (1 - c) X
[w (x))(()))
X-P
T(B 1)
Regarding runtime, for any L r/(eT) have that (xxr(wx(XX)(
Ex-p (XX(XXX)(X)XXXXXXX)) < r < T, where we used the inequality de or S 1/ a). ot
Lemma 2.5 (Multi-Direetional Filtering), Consider Setting 2.3. Given a nearly-orthogonal matrix
U Rkxd satisfying Ullop 2 tr t)) log(1 [/)), there is algorithm that reads e, the such that
points, their weights w(x), and returns weights w' in time ndk: + polylog(d/c, | U
(i) Ex~o|u(X) w/ (X)] (e/ log(1/c)) ExaBle Ex~e(w(X) ( wx)).
(ii) e <Ex~Blu' (X) _p. ) S etr(U U).
22

---

Proof. Let g(z) = uPm>1 and r(x) : g(z)11p(x)) > 100tr(U'U)) U)). Denote Bp.. =
Dp Using the goodness conditions and the certificate lemma:
E
1 (X))((())
- log(1/c) |0 + 10( - -rr
(12)
X~G
log(1/c) + 11 - p. 12 ||U|
(by Lemma B.25)
2e
+ + +0 (lope) |
(13)
log(L/c)
log(1/c) + +c llop |
(14)
log(1/e) + c)pooyyoo((1/)) |
(By Setting 2.3)
|
(15)
log(1/c)
where the first inequality uses Lemma B.25 (this uses that |0 pr...113 N tr(UUU which
follows by the certificate lemma and that BB || < 1), the second inequality uses |(x =)ll tr(U U(z px 44) and then the fact tr(AB) I tr(A) [B for A U U.B (x pz T)T, the third line uses the certificate lemma (Lemma 2.4).
We can thus apply the filtering lemma (Lemmaa B.27) with T equal the right hand side of Equa- tion (15), and 8 = log(1/c) to obtain that (xxpp/ww(X)rr) S TB & U| In more detail,
E (x)()))( X~B
sc [w' (X)) 100tr(U (U))
X-B
w (X) (X)) + utr(UUU)
$ T log(L//) + etr U UU)
(by Lemma B.27)
x U)
This completes the proof of part (ii). Part (i) follows directly from Lemma B.27.
B.6 Algorithm from [DKKLMS18]
We now describe (a simplified version of) the algorithm from [DKKLMM18] and state its improved runtime.
23

---

Algorithm 4 Optimal Error Robust Mean Estimation Under Huber Contamination
1: Input: Parameters e, r, and sets S1, S2 C R where for each i 1, the uniform distribution
S2 is of the form (1 e)G, cB; for some G2 satisfying Definition 2.1 with appropriate parameters.
2: Output: A vector in IRd
3: Let k := r log(1/e), and P be the uniform distribution on S1.
4: Calculate naive estimate ja s.t. || /l22 R for I O ( log(//))) e.g., geometric mean
5: Initialize t <- 0 and (x) 1[I مll S 2R) for all
Rd
6: Denote by P. the weighted by w version of P. by Σ, the covariance of Pt.
7: Let vr.1,...;; Ut,d be the eigenvector decomposition of Σ, and de,1 ... Xt, the correspond-
2
ing eigenvalues.
8: while Ei=1 Ae 21+ A1+ for x = ed do
9:
_ Define the polynomial gx(x) Ue(x pt where U, = [....1/// /o.i/Vi,... .v....,, Ut //k/ Define the score function (x) = (1( (g) 100/r)
10:
Wt+1 < Filter P. c, w,, T²)
b Algorithm
12:
t t 1, and update vectors v4,1.... U2
Let V sp.n( } and y u vaj.
14: & Proju- (Ht).
P2 f BBrutFooceeSs V)
16: return μ1 + P2.
Lemma 2.6 (Adapted From [DKKLMS18), Let be set of i..dd samples from an e-corrupted
version of N (0, Ia) and let r, (0, 1) be parameters. If (og())))) (oo//// /polsoo//o then there is an algorithm that, when having as input T, E, and r, after time Ö(nd² + कه7) ), it outputs vector µ E Rd such that 1| ull2 S e/r with probability - 8.
Proof Sketch. The algorithm consists of two parts. The first part (while loop of Line 8) finds weights for the points Wy (x) and subspace V span (aan(u11-. (i....,,2 Uk such that the weighted empirical mean approximates the true mean in that subspace within éi-distance O(e/r), and the second part runs brute force procedure (Lemma B.28) to approximate the true mean in the remaining k-dimensional subspace. The brute force procedure runs in time exponential to the dimension of that subspace, but since k: = r log(1 /e), it becomes polynomial in 1/c. Concretely, the runtime of the brute force step is given in Lemma B.28. The dataset S2 that cardinality: since the dimension of the subspace k, n' = (kk/²²) samples in S2 are enough satisfy the goodness conditions. Thus, the runtime from Lemma B.28 with replaced by n' is 2 "
for this brute force step is of smaller
In the remainder, sketch the runtime bound for the first stage of the algorithm, by arguing that the number of iterations is (d). If we let our algorithm compute exact eigenvalues and eigenvectors, each iteration needs nd time, resulting in total of Ó(dd') for the runtime of the while loop of Line 8. This can improved Ó(nd²) by calculating approximate top k-eigenvectors (see, e.g., [MM15]).
To show that the number of iterations Ö(d), we will show that every time that the filtering procedure of Line is used, removes least (2(e/)) mass from the points' weights. Then, given that the outliers themselves have mass at most (e)) the algorithm would stop after Ö(d) iterations. The mass removed per iteration is (by definition of the filtering procedure)
(X) Wc+ (X)) Ex~p(w (X) X~P
E
max (x))00 Ta (x:)
The denominator is upper bounded by max (x) (x(x)= < Ō(d) by the prepossessing of Line 3, 10 For the numerator,
[w [X) x
E
E X~p
(x (X) < 100/r))
X X~P
10 Since the norm of the samples from the Gaussian distribution is tightly concentrated, the filtering procedure removes ole/ log(1/e)) fraction of inliers.

---

The first term is at least (1 A) because of the condition in Line 8. We also claim that the second term
most (1+)/2). These two would imply that Ex~p(w,(X] -X-) 1 = ?2()) = (2()), which implies that each iteration removes D D(e/d) fraction of points, completing the proof. The remaining
claim for bounding Ox>>l[w(XXeee (X (g) (X) < 100/r)] is similar to [DKPP23, Lemma 2.14
(iv)]: We first note that E Ex-a (X) ) (ge (X) & 100/r)] = 1- O(e/ log(1/e)), because at most (cc/lo(1//c) (og(t//))-fraction -frcciinn of inliers have gr (z) > 100/r (cf. goodness condition (2.c) in Definition 2.1).
Then,
E (X) S 1000r)]-- Xap
- (1
E X
S 100/r)]
+ E
wu (X) 1 100/r)] -
The first term is upper bounded by X/4 by the goodness condition (2.b) along with the certificate lemma (Lemma B.23). For the second term we can use the bound
E (X) < 100/r)] - 1 S (100/r ((er))<(MMM)))) 1a~B
-
We will use the following brute force algorithm to estimate the mean in a low-dimensional subspace. Lemma B.28. There exists algorithm, BRUTEFORCE, that given set of n samples in I& from e-cormupted version of a distribution G satisfying the first condition in Definition 2.1 and a subspace V of dimension k. runs in time O (dnk (n n + (k)((k)))0(k) and finds vector ja E TR such that [[Proiv(w m O
Proof. The algorithm first projects the points to the k-dimensional subspace. then computes an n-cover C of the k-dimensional unit sphere, which by Fact B.22 consists of at most (4/r)* vectors, and for every vector u E C, it finds an estimate fue for the u-projection of the mean satisfying Ju" j4-Ma pul I for r OEe) being the same as in the right hand side of the first goodness condition in Definition 2.1 (because of that goodness condition, letting uu be just the median the u-direction achieves y-error). It then solves the linear program
find x s.t. |w x pul < 27 Vu EC
Let E cd be in the direction of x μ and arg min ||I u1l2- Then,
|| nll2 To (x n)] |u (x n)) 1(u u) (x 14]
< w (x - -mu) + |u (wo - m) + ||I - ull2 ||x - 1ll2
47 + 7 || Hll2
(16)
which establishes that ||x- wl| < 47/(1
Regarding runtime, the projection of the data takes Oddnk) time, solving the median in each direction takes ((44n))" n) time. For the LP, the separation oracle that checks all the constraints exhaustively needs time (4/1)** and the LP can be solved by (poly(k) log(R/r)) calls to that separation oracle, where R, r are bounds on the volume of the feasible set 1 vol( {z e R* |u. x Hal S 71 Yu. E C}) < R. For the lower bound we have that the volume of the feasible set is least the volume of k-diiennsonna ball of radius 7 because μ feasible solution and + yv remains feasible for any unit vector v. For the upper bound, we have that the feasible set is included in a k-dimensional ball radius 4ry((1 - n)) by Equation (16). Thus the runtime for the LP is (poly(k) log( (Rg())r))) = (ppol((k)). Finally, choosing 1 1/4, simplifies the total runtime to O (dnk + (n log + poly (k)))
C Omitted Proofs from Section 3
We begin by stating additional conditions that must hold throughout the algorithm's execution. Subsequently, in Appendices C.2 and C.3, we present comprehensive analysis of the two cases in our main algorithm. Finally, Appendix C.4, we combine all the components conclude the proof of Theorem 1.3.
25

---

C.1 Deterministic Conditions for Algorithm
Condition C.1 (Deterministic Conditions for Algorithm 1). Consider the notation in Algorithm 1. For all t E [fmas], the following hold:
(i) Spectral norm of B A E [0.1] B lap, 10 B0 [lap]-
(ii) Frobenius norm ull E [0.8k] M 1i 1.2k | M; [1]]
(iii Scores: For all support(P), | (z Me /
M, (x - pe)1 [0.8k, 1.2k].
log tr U.). (iv) Maximum norm: maxje(4) |/. 112 -
(v) Probability estimate: 1 qel - 0.1/ (k2 tmax).
(vi) U, almost orthogonal: For every pair us uj of distinct rows of U, it holds Ku,, uyبل
Uy Iallll 1/ll///
(vii) Approximate cigenvector: u, from Line 19 satisfies uf Bi u, 2 (1 - 1/p)||
Lemma C.2. Condition C.1 is satisfied with probability at least 0.9.
Proof. The spectral norm condition holds by using Fact B.10 with 1 = 0.9, 5 0.000 tmax and union bound over the fmax iterations of the algorithm.
The next two conditions hold with probability 0.999 by using Fact B.14 with A = M Vh and a similar union bound.
Regarding condition (iv), with probability 0.999 we have the following regarding the vectors Ur from line 15
||0, 12 SosSVosss /lo((k/mxx) (M, ) IIF + log(kkl
)?lpop
log(k)aax))
log(A) (k)
S log(kk [|M ||
S
| 11 -
1
k
k
where the first step uses Soanson-Wright inequality with probability of failure 1/ m11000ktmax) the second step uses the fact || < tr (((A)) ||) I tr(A)² applied with = (M, )2, the third step
is due to condition (ii), and the last one is because fmax & poly(k) for the choice of values that have made for max and k.
For Condition (v), if Tt is the empirical average of ( 2, M > over (k. m instantiations of 2, ~ (O0 I), then the condition holds by a standard Chernoff bound.
Condition (vi) holds by the fact that the algorithm uses U, only in the case that q. S 1/(22tmaxi (max). By condition (v), in that case we have that (r,,,'~N(o,ri) (o, (M, 2, M, 2²)l > |M, 2||2|| </112/12) 2
0.9/(k2t maxx) Thus, by union bound, with high constant probability, for every iteration all pairs of rows of U, will have angle with cosine in 1-1/k2,
Regarding Condition (vii): By Fact 1.10 with 77 1/pp 1/ log(d) and 8 = 1/tmax (so that can do union bound over all iterations of the algorithm), have that if p' > (dt then n, B 2 (1 1/p) B llop with high constant probability throughout the algorithm.
C.2
Case
This section contains the details missing from Section 3.1. We restate and prove the following:
Claim C.3 (Claim 3.2 restated). Assuming that Condition C.I is satisfied, for every iteration t that the condition of Line is satisfied, Φ141 S 0.99%). Moreover, Bxxg[w,(X) 2 --e //og(1/)) and B = 0 throughou the algorithm's execution.
The reminder of the subsection is dedicated to the proof of the claim. In particular, we focus on the claim + S 0.99é,, since the invariant Ex~o[w:(X) 1- e/ log(1 /e) follows by the guarantee of the filter (Lemma 2.5) rather directly: First, note that by Lemma B.2 B.26, the deterministic
26

---

conditions (vi) and (iv) from Condition C.1 imply that the goodness condition (2.c) is applicable
to the polynomial ||,,(x: -1) where U, the matrix used in Line 15; thus, Lemma 2.5 is
applicable. Second, by that lemma, we have that every time that (x(xx) changes to W+++ (æ) it holds
that Ex-G(w: (X) - (+++(()) (X) < (e/ (1)) Ex-B (x)(() 204+1 (X)]. If t. denotes the final
iteration, this means that E G (W X) w2- (X)] < (c/log(1 (()) Ex~B (wo (X) we (X)) -
e/ (1/e), thus Ex~a[w> (X)) 2 - e/ log(e/e) throughout the execution.
We also note, that Ex~c//u [X]] 2 1 c/ log(1/c) implies that B, / 0, and B, = 0. This is
because, by condition (2.b) of the goodness condition, it holds Σ, - (1 - O(e)) Also, we have already shown that Ex-a (wy [))
e/ (1/))/)). Thus B, == (Ex- (wo (X)) "
C (1 O -1- (1- C 0, if C1 is large enough.
Regarding the potential function dr, we will establish the following series of inequalities, for c being some constant of the form c = O(1/C), where C' is the constant used in Algorithm I (we can assume < 0.0001 by selecting C' appropriately large in the algorithm):
1141 tr((M )3) |
3p
2p
2p
d (
M |M++ 112p+1
(by Fact B.2)
d=st [IM
2p+1
10+1
- 11
1011
$ d (tr(M, B M,
(by Fact B.6)
$ d [c]
|lop B
22y
(by Lemma C.4 below)
$ + c क
2u
([ 12
B
12p
1p
2n
(by Fact B.2)
- 0.99 B, 12gp
0.99t...
(since p log(d) and 0.0001)
Where the application of Fact B.6 uses that B2+1 = B₁: this is by definition of the covariance matrix Y) Y which shows that downweighting points can only make B, the matrix smaller in PSD order (note that B, has the covariance matrix multiplied by (Ex~ (x)x))))())))))
1
as 2(Ex-P(wx(X)]] (x)))² Ex,y Px.y~plwe(h
remains to establish the main inequality that we used above:
Lemma C.4 (Filtering Implication (Lemma C.4 restated)). In the context of Algorithm 1, and assuming that Condition C.I holds, for every round E [fmax] that the algorithm enters Line 13, it holds (M, M²) S e' |Mo f for some d' O(1/C) (where C is the constant appearing in the algorithm's 's pseudocode).
Proof. First, note that
E+t
E X~P++1
(Proiva, (X) - (tt1) (X) - =
E X~P++1
(Projvy (X) - - (X) -x)-14 |
Ex-rlwitt X~p [wee ( (X)ov) (X) - 1 (Proive (X) nt)
E
(17)
To avoid the repetitive notation of Projys, (X) in the remainder of the proof, we introduce the shorthand notation X denote the projection of the random variable X 1/41 To avoid confusion for later on, we also note that since we are in the case that the algorithm enters Line 13, V4+1 Vit. i.c., no change is done to the subspace and in this subsection it does not matter which of the two subspaces we use in our notation.
We decompose tr (M, B B,L, into the contribution from inliers and outliers: Using the definition of B = Ex-Plux41 (X)] (1 U and Equation (17), we have that
27

---

tr(M Bi+ )+))
(18)
$ E [X)] E (0 (X- mt (L
(Ie) 14 MF)
(19)
X-p
x
< (1 -)) E (X) - 11)
+
E
[x)+
14 )1
(20)
x
X
-(1 (1-CGe))((M))))
1.25c
$ E (X 14 )12) +
E w
[ - -1. )131
(21)
- -1 -11-Cic)IM;
(22)
where the first inequality uses the definition of B₁, Fact B.5 and Equation (17), for the last term in Equation (19) we used that Ny; M, M which can be shown as follows: recall that using
then since Hy is idempotent, multiplying that expression from any side does not change it. For Equation (21) we used 1 - < 1 and the Johnson-Lindenstraut sketch guarantee (condition (iii) of Condition C.1).
our definitions MA-- (B²)²-
[I B, P
We start by upper bounding the combined contribution of the first and last term. We will need our assumption that G (the uniform distribution on inlier samples) satisfies conditions (2.b) and (2.a) from Definition 2.1 with a = / (og(1/e)), order to apply Lemma B.24 with U M7 and b = Nt:
E Mee 1 (X) M, (X! pe I1 (1 Cie) M lR X~G
(23)
$ + 0(c)) M, 11 + |IM, (p' - -M4 O(())M 1Ell 1e 112 - -(1 MM 11
(24)
< (1+0 |M, |M (lope) ( (VIB+ - ([
11
(25)
S [O
(26)
where Equation (24) uses Lemma B.24, Equation (25) uses [M p MM
(z )T as well as the certificate lemma (Lemma 2.4) which gives bound where it can be checked
) ) and then the fact tr(AB) S tr(A)][B op for A M, M, B = (a+ )
that the dominant term is O / B, lap), and Equation (26) uses that e < {OI llop due to Line 9 of Algorithm 1 and condition (i) of Condition C.I.
We now focus on upper bounding the middle term of Equation (21), where we will use the guarantee of filtering, in particular the second part of Lemma 2.5. By Lemma B.26, the deterministic conditions (vi) (pairwise angles close to 90 degrees) and (iv) from Condition C.1 imply that the goodness condition (2.c) is applicable to the polynomial |,,(( - 11)l where U, is the matrix used Line 15. In other words, the requirement of Lemma 2.5 is satisfied and the lemma is applicable to case. We thus obtain that
1
E X E (X) n² ) S
& | S <|MI S S lB |M,
where the second inequality uses the Johnson-I -indons--indssss guarantee (condition (iii) of Condi- tion 1) and the last inequality uses that e < ( |lop which is due to Line 9 of Algorithm 1 and Condition (i) of Condition C.I.
C,3
Case
This section contains the details missing from Section 3.2. We begin by proving the formal version of Claim 3.4 below:
Since we are working with the projection of the data to V7- this lemma requires that the goodness conditions hold for the projected data, which we have noted before holds (see Remark B.1)).
28

---

2M =1) ||M M/ 2'12 > 7) 2 a, then
Claim C.5 (Formal version of Claim 3.4). If Pr. Pro,s
[IM 1. tr((M P)
2 7
Proof. For some parameter > 0, define the following events for random variables 2, 2. & {2:
||M (r(()))))), and &g = {=: > )r( We have the following series
of inequalities:
(M, 2, M, 2
a E
Pr
>
(0.
M M 2/ 11
(M, M x/)
Pr
M A + Pr
(10.)
/|2 z' ||2
2, 2²)
Pr
AA4 &
+2v65
(27)
(.0
tr((M)) )2)6
(M 2, M, 2²) (r((M
Pr
?>
2V86
(0.1)
1 Eg (O. 2
22
(28)
?
(2(22)
+ 2Vee,
(29)
7 (r( )2)262
where Equation (27) is obtained by union bound along with the second part of Fact B.15, Equa- tion (28) obtained by squaring and using Markov's inequality, and Equation (29) follows by using the first part of Fact B.15.
Picking the value d a² /1000 and reorganizing Equation (29) implies the lower bound [M,,) Pelle/t((MM))) )²) 2 2 Combining that with the upper bound A I Ar( for any PSD matrix A (Fact B.8), we obtain:
II ) tr((M+))) tr((M,))2
((M+)
S
(30)
We now turn the second claim, using the same notation in the pseudocode. The requirement that B, is positive semi-definite (PSD) has already been established. As shown in Claim C.3, whenever the weights are changed, we have proven that B, remains PSD.
Claim C.6 (Claim 3.5 restated). Assume that B, is PSD. Let V V U (we) for some unit vector We Vit, then de+1 6:- u, )M
Proof. We first examine the effect that projecting everything to the space perpendicular to u, has on the matrix B
Claim C.7. Define A: := I u.u. We have the following:
(i) II 141 I A,
(ii) E A Σ A.
(iii) B = A, A.
Proof. We prove each of these separately.
29

---

Proof of (i) We first recall the definition of the orthogonal projection matrix: Given any basis et.1..... Dgr F of the subspace Vi (r is the dimension of the subspace), if A, the d matrix that has the ut, ss as its columns, then the matrix that performs the orthogonal projection to V is II == (.)AAA))
Consider the basis that is orthonormal and starts with 1/ = ue (the vector mentioned in the claim's is statement; such a basis exists since E V, Then, II, = A,A, = u4 + Ei= Vt.-ug,. Now, let | be the matrix that projects to the subspace - we). Since uy itself belongs to Vi=, an orthonormal basis of V4+1 is 01,2 vs,,, and thus if we define Act1 == [V,,2..... ,t,r then the orthogonal projection matrix for Vi+a is I = A = Ei= V, We claim that this is the same as A, A,: Indeed, A,I, A, = (A,A))
V4+ == {- E V/==
and A, (I 2 A is the matrix with zeroes in the first column and -....-... Mer in the rest of them, thus (A,A) ( A,) = - = "0,,%
Proof of (ii) Let A, be exactly in the previous paragraph. For the specific orthonormal basis 1t,1.-..; Uir E IRd of V (where r is the dimension of that subspace) where the first vector has been chosen to be 11 U4, define the d x matrix A, [01,1, U1,2,-... pt.r] [Mee U2,2. -.. Mt.l
Also, let be the mean of Pt. Then,
A A.
E X~Pr
A.(Projy, (X - -()) (X- pe)) A
E
A (X - (x) - 144) II, A.]
X-P.
E
A (X Mi)(X - mO) A.A A.]
(31)
X-P,
E
AAIIATI (X - -M)(X pe) A141
X-P.
E X-P.
(Proivi, (X - (p))(Po (X -
(32)
for the last line used the expressions for A,, Aet1 that we gave in the previous paragraph to conclude that A = 1,
Proof of (iii) The proof follows by the definition of B = (Ex~plu,(X])) 2 - (1
(e) " the previous two claims, and linearity.
We now show how to complete the proof given Claim C.7.
ΦL+1 == tr (BB/+1)22) )20)
tr ((A,B, 32)
(by Claim C.7)
(BB/ A.)
(by properties of A, and trace; see below)
S tr (BB A
(by Lieb Tiir-Thiring incquality (Fact B.7)
tr ((B+))PA))
(A, : uuu, is idempotent) (by definition of A,)
tr ((B))²² (I wen,.)
(BB+) )27) tr ((B,) 2 (P
tr (BB )r)) uu (B,(BB 2 (2.))
(by cyclic property of trace)
A A, and observe that: (i) A A, since A, I "nin and (ii) using the cyclic property of trace, after applying trace to both sides move the first A, to the end.
where for the third step one can expand out (A,B) 2,) A,B! A,B, A
30

---

Finally, we show how the above two claims imply that Φ1+1 - (1 polylog(e/d))o, for every
iteration that Line 19 of the algorithm is executed:
0++1 Φ 2 u. (B: u
(by Claim C.6)
(u B 3p
(by Fact B.1)
2p
B, llop
(by condition (vii) of Condition C.I)
| op
= || 1
(2)
2N
(by Claim 3.4 with a, 7 = polylog(d/c)
polylog(d/)) (d/c)
polylog(d/e)
C.4 Combining Everything Together
We now describe in more detail the algorithm that achieves Theorem 1.3 and put together the previous parts to derive the guarantee and runtime of the main theorem.
The final algorithm is outlined Algorithm 5.
Algorithm 5 Full Algorithm
1: Input: e E (0, 1/2), sets S1, S2 C IRdd datasets, where for each = 1,2 the uniform distribution S, is of the form (1 G)G; B, for some G; satisfying Definition 2.1 with appropriate parameters.
2: Output: A vector in Rdd
3: Run the algorithm from Lemma B.18 on S and let wzz) for z E Si be the weights that it outputs. Let P be the distribution assigning mass wxx)/ (xx)/>rEs) w(x) to every x E Si.
4: Run Algorithm 1 on input P. Let V and μn be the subspace and the vector that it outputs.
5: Project all points from S2 to V1, i.e. S2 f (Projy, (x) : x E S2).
6: Run Algorithm& 4 on S2. Let #2 be the vector that it outputs.
7: return μ1 + P2.
Proof of Theorem 1.3. The set S, is comprised of i...d. points from the e-corrupted version of N(s, I) and its size is chosen to be a large enough multiple of e log(l//)) Lemma B.20 the uniform distribution on S₁ is (e, a, -)ooo with a e/ log(1/e) and k polylog(d/es (d/e). The preprocessing step of Line 3 ensures that the eigenvalues of the covariance matrix are most + cpolylog(1/c) which part of Setting 2.3 and utilized in the analysis of the next steps. The preprocessing step of Line 3 runs in nearly-lineat time (ndpolylog(d/e))
so that by
The step of Line 4 has been analyzed in Section 3 and Appendix C using the potential argument. We have demonstrated that in each iteration of the algorithm, the potential function s; == tr( (r(())))) undergoes multiplicative reduction: Φ++1 (1 polylog(e/d)jó. Initially, the potential is naively bounded by do poly(d/e)
With just polylog(d/4 (d/e)) iterations, the potential decreases to o, I d' resulting in [B, llop O(e) and leading to the termination of Algorithm with [|Projy: (Mi - 22 = (e)).
Each step of the algorithm can be implemented in ((ndpolylog(d/c) time, for example, computing (B-)logd d by iteratively multiplying by B-. Utilizing the special form of covariance matrices enables each multiplication to be calculated efficiently in O(nd) time.
Since Algorithm 1 only adds one direction per iteration the subspace V1, the dimension of Vi will be d' = Oppolyoogdd// (d/c)). We now use a smaller sized second dataset S2 that consists of 2dd + log(1 (o))p))llllgggggg (d/c) samples and run Algorithm 4 on S2. The runtime of that algorithm is given in Lemma 2.6, and since we are using d' = O(polylog(d/c) for the dimension, the runtime Finally, since the output #2 approximates μ in the subspace Vi within error O(e), combining μI and μ22 yields an estimate with [| + u22 - =ll2 S [|Projy (w) - =2]l2 + |[Projy, (N2 4) 2 O(e).
of that algorithm is O "olyoooollg((////
31

---

D Omitted Proofs from Section
In this section, we provide the details omitted from Section Appendix D.1 proves the properties of the conditional distribution of the covariates conditioned on the responses being in interval. Appendix D.2 then uses these properties to show that the samples from this distribution satisfy the goodness condition. Finally, we provide the details missing from the proof sketch in the main body in Appendix D.3.
D.1
Conditional Distribution of Covariates
We first begin by proving Claim 4.1.
Claim 4.1 (Conditional Distribution). Let (X, E gd+t1 follow Definition 1.2. Given E R.
denote by Ga the distribution of X given y = a. Similarly, given an interval I C R, let Gr represent the conditional distribution of X given y I. Define oz o² - || Then, Go 3 33 and G PFiyET ,, é(a'; 0, oN la ,3 )d', where o(z; 0² denotes the pdf of the NN0, 1²) at E R.
Proof. Using Fact D.1 with y1 = X, 1/2 U.11 12 0, E11 Ia, E12 B, E21 =3T, Ez og =0= 0² + ||312.
Fact D.1. If
111 N
P1
11
then y1 v1lu2 N(a, E), with p pi Σza
12
1/2
Σ
E12 Eng (y22
U2) and Σ11
Σ12E
Ez1-
For the distribution G₁, the claim follows by definition of conditional probability and law of total probability
Px)vrr(>) Sa'er Pxly (a))yyaa)daa
Prly I
and then the expression for the pdf of Ga that showed earlier.
D.1.1 Properties of Conditional Distribution
In Section we mainly focused on the case where the underlying distribution is (1, La). For the purposes of our robust regression algorithm, we show in this section that the goodness conditions are satisfied by the distribution G, and get corollary a robust mean estimator for pG,,
Lemma D.2 (Properties of G1). Let I ja,b) b) be an interval of length f. Further suppose that max [6]) < 20,, Let X ~ G1 and let and E be the mean and covariance of X. Then:
1. p - 00/3//0 < e118112/o, and Illll S 411,812//r,
2. [ Illpp $ /913112/0?
3. |X /llvz - + 13112/00
Proof. For any z E I, we know that the mean of Gg is equal to 28/0² by Claim 4.1. Since G is a convex combination of G2 for E 1, follows that μ is convex combination of aB/o and 6,3 Therefore, the mean μ satisfies || u3//2]11 < 1b al11///////
Similarly, for any 2 E I, Ew~G. [WW = Ia 33 + - 38 by Claim 4.1. Thus the second
moment X ~ Gy is equal to Ex-G, [XX Ex~G|[XX]]
= I. 33 +1-+(11o00 33T for a E [0, 1].
Thus, have that
Σ - Tallop
E XNGi
[XX]] 1411"
I
op
32

---

1
1 + of 2111112 + max(a² baxx(a,)>)1111133 1181//////
3 of S 91
S max (2 bb
The proof of the last claim is based on the fact that each component of the mixture Gy Gaussian (and thus sub-Gaussian) For single one-dimensional Gaussian Y ~ (N o20, the sub-Gaussian norm is constant, particular || 12/T..
1/p
X~G1 [|v (x ma,)p1
E
E E XnGa
[ (x -c, [0]
max E X-G,
[|v (x pG,))
[| (x mG.)) 1/p + Il - MG,ll2
max E
aet X-G,
S 22/T e | lo
Hence, (x - pG:))l S 22/T + || o for every direction E Sd-1
L
D.2 Robustness of Goodness Conditions to Spectral Noise
We show that large set of i.i.d. samples from the conditional distribution G2 satisfies the goodness condition provided that the interval is close to the origin and the norm of [ is small.
Lemma D.3. Let = ba be an interval of length Let eo be a small enough constant. Let E (0, Further suppose that max([al, (6]) < 20y || S de log(1/e), and a 2 e.
k (d (oog(//)))) then with min(a?.e ee/ log=(1/e))
Let S be set of n i..d points from the distribution G If >> probability 1 T, the uniform distribution S satisfies (c, a, -gogoddssss condition (Definition 2. with respect to PG,,
Proof. We show this for each one of the conditions separately.
Condition (1) We first calculate the true probability of u TT(X (X- mG,) being larger than The distribution of u - (HG) is NN ā.) where o. := (v B)² Using Taylor approximation for the edf of the projected points.
t
Pr IV (X - (aa)) 5, S X~G.
22T
2v2n
Thus,
max Pr To (X (GG)) 1ō, + max #ES4-1
||GG, pG [1]
XmG1
ae1
max Pr (oT(X (X (G tō max IING. mGxl12]
ael XNG.
aef
max Pr [v [uT(X- (X MC.) > tō
ael XwGa
1
2
2
2
Applying VC inequality as in Lemma B.20, we obtain that the empirical probability over a set of size (d log(1//))/ samples will satisfy that with probability 1 T,
1 1
1
max Pr T (X - (0G,),,) #ES4 X
max o. + max ||MG, - HG, 1ll s ae!
+ c.
vEsd
/2m
2/2x
The expression on the right is less than 1/2 for 2 c. Thus, we obtain that Condition (1) is satisfied with 71 S e max ||GG. MG,l12 5 fB/o
33

---

Conditions (2.a) and (2.b) The proof of these two conditions for the Gaussian case only uses its concentration. From Lemma D.2, we have that || HGG pl S 2 + é / which is less than constant when l S o and || 2 a.
Hence, identically to Lemma B,21, we have that if n > (d + log(1 (1)(///))//a, for any set S' C S with |S|| 2 (1 a)]ss, the sample mean As and the sample second moment Eg satisfy that Ilss pa, |ll S a (lglo(1//)) and |S - Es, lop X a log(1/a).
Using the triangle inequality, we have that ||S - Iloo < S oop + |Es - Σa, llop X log(1/a) + e² log (1/e) S a log(1/a), where we used Lemma D.2 and the assumptions on || and a. Thus, Conditions (2.a) and (2.b) are satisfied.
Condition (2.c) We will follow the same notation as in Lemma B.20 and claim 4.1. Observe that the proof of Lemma B.20 relied on Equation (4). Our goal will be to show that similar right tail concentration inequality is satisfied for the distribution G2. Since G2 is a convex combination of G2 for 2 E [a, b), suffices to show that each G2 satisfies Equation (4) after small changes.
Let X ~ G.. Then x = E// y μa for ~ IO, and μæ = 23/0² and Σ, = I - 3,3 /og.
We will show that the polynomial (x) Ux p1)ll satisfies similar right tail concentration inequality as ||Uy|] 122 in Lemma B.20. Applying triangle inequality,
(u( | 11)11 - 2 !/2 v12 + 2|| - (1)))
< 2 2 v + 2 11 - 13
By Lemma D.2, |M= - - 1 is at most e 113/1//0 S under the assumptions. Thus, we obtain the following deterministic expression:
pu(z) s U /2 v1 2tr(U U).
(33)
The expectation of ||E!!" y is ( U U = tr(U" UE.). which less than
since Σ 1 I. Combining this with Equation (33), obtain
P(po(x)= 2 100tr(U T U)) S U 1/2 wl
EI /y Z 48tr(U U)).
(34)
Finally, the right tail of ||UΣ! /v13 - EI| 1/2 2 depends on / Uu! I|
and
|x! U UUU!/" llop by the Hanson Wright inequality. Both of these expressions are again bounded U llop. respectively follows:
by ||U UIlF and I
| UUT!! <
|
11op
1lop
|ITUUFF
and
| UI 121 llop
E
op Σ 1op
S
||
Ullop
where we use | olop I 1 and Fact B.9. Hence, applying Hanson Wright inequality to || 1/2 y/1 along with the above relations, Equation (34) implies that pu (x) satisfies Equation (4) with (c/k)² in the right hand side. The rest of the proof goes through similar to Lemma B.20.
Thus, Theorem 3.1 is applicable and we obtain the following algorithm.
Corollary D.4. Let Eg be a small enough constant. Let e E (0, c0). c E (0,1), (0, 1). Let S Rd be an e-corrupted set of n >> d+log(1/7) polylog(d/e points from the distribution G3 of Claim 4.1 with 8 S a/ log(1/e) and || S de (le). Then, the Algorithm 5 produces an estimate t such that || - μa, 2 with probability at least 1 in time (nd + 2 d, 1/e, 1/T).
D.3 Proof of Theorem 1.4
We now provide the details missing from Section 4. The proof sketch there relied heavily on the events defined in Equation (1), restated below:
(i Jo -og1 < ofe log(///)
(ii |a| > 0.00016,
(ii) 11, pG, 12 S E/cc (35)
34

---

The details missing from the proof sketch in Section 4 are: (i) the events in Equation (35) hold, the events Equation (35) suffice for the algorithm's output be correct, and the proof of Lemma 4.2. We already explained why (i) and (ii) hold. The third holds by Corollary D.4.
We begin by formally showing that if the events in Equation (1) hold, then the algorithm's output is correct.
Claim D.5. Consider the setting and notation of Theorem 1.4 and Section 4. Let E (0, 1). Assuming that the following three events hold
(i) 0 ولل S ofe log(1/e).
(ii) |a] > 0.00016,,
(iii) |3, - pG, ll S e/c.
Then, ||3 - B|1 x ae/c.
Proof. We want show that given Items (i) to (iii), we have that ||3 - S oe/c, where = (G)//a)3 is the algorithm's output vector. Using the triangle inequality,
6
G2
oy
3 8
B1
31
PG,
+
PG,
3
(36)
We upper bound each term by O(ae/c) separately. For the first term we have that
02 B1
62
3, - na,
PG
la|
2
0g
B1 No.
(by Item (ii) above)
Tye/C
(by Item (iii)
OyE/C
(by Item (i) and 1/2)
+ 11811))ecc
(o) = 0² + ||31|3
(a oe wowooll////olcol logg S oe/c.
(by assumption
|| x lo l/c)
(e < 1/2)
For the second term of Equation (36), the triangle inequality implies that
02 y
o
+ 1 o
PG,
B
<
NG,
8
||MG,
1l2.
aa
2
We bound the first term above as follows:
o² |
PG,
PG.
- 3
a
$ Syy
a 8 PG,
σ²
(by Item (i))
2
2||3112
0gy
(by Lemma D.2)
|||81|2
< oe.
(by definition of f, Item (i), and assumption || oe oelog(1///))
The second term can be controlled as follows:
o22
| |l2 H oye log
(Items (i) and (ii)
lal
(a + 1)118112
X log(1
(by Lemma D.2)
a3
35

---

o log (0g()))),
(a ta+250 Jg by definition of a,
2A 23 e
(ogg(/)))l33))
(by Item (i))
ac.
(by assumption || 2 ac aclog(1/c))
Finally, we provide the proof of Lemma 4.2 below.
Lemma 4.2. Consider the context of Theorem 1.4 and Algorithm 2. First, for every possible choice of
that can be made in Line Pr(x g y E 1] Z f/0.g Second, E, (r( E [[]]] S E/o where the outer expectation is taken with respect to the random choice of the center of T.
Proof. We start with the first claim.
Pr y E [ 1]= (X.)~G
Pr W E
(0.00)
( E Do((: ofjdy
(olv; 0, of)
denotes the pdf of (O og))
VER
( E T)((; df
yel-1.10,,1100
[ Iel--loo
(0.2/ax) 0.28/oy,
(v E fjdy
where we the second line uses that I g-1.lo y,,oy] because its center belongs in [-ō ōgy and its length is = /,/ log(1/e), and ō, is very close to σ, (Item (i)), and < 1. In the third line we used that the pdf of (N(0, of) pointwise bigger than a small multiple 1/o,, in the range [-1.1o7 I.
We now prove the second claim. First,
Pr yy 1 Pr (Xy))B
I, E 1- l,ō +]]+ Pr ly I, $ (-ō, L, E]]
(X.y)BB
(Xy)a
(37)
s Pr y 0- - l.o, +[] (Xx)=B
(38)
< Pr ly e (Xa)>B
E 6 1, {
(39)
Now let B' be the conditional distribution B given y E 1-ō, - L, ō. + {] for saving space. Recall that the interval / itself random interval. If we take expectation with respect to I of both sides of Equation (39) we have that
Pr [y E S 1.X.
E
1g
E
E (a.g)~B' [
Bln(y mi]
(x.g)~
l
E
(x.y)~
2.20 y
2.20
where we changed the order of expectations, and used that given fixed point E 1-ō,, - E, + E. the probability that this y being hit by our random interval I is most the length e of divided by the the length of [-ō,- - e, G, + f] (which as we have shown before is subset of [-1.1o,, 1. In the last inequality we relate og and oy once more using Item (i).
E Adaptation of Algorithm 1 to the Streaming Setting
In this section, we give a brief sketch of how the robust mean estimation algorithm of Theorem 1.3 be adapted to the streaming setting using the techniques of [DKPP22). The adaptation follows closely [DKPP22] and we thus omit the details.
We recall the data access model for the streaming setting:
36

---

Definition E.1 (Single-Pass Streaming Model). Let S be a fixed set. In the one-pass streaming model, the elements of S are revealed one at a time to the algorithm, and the algorithm allowed a single pass over these points.
The samples that form the set S are still samples from an r-contaminated version of ( I) as before. The adaptation of our main result to that model is the following:
Theorem E.2 (Almost Linear-Time and Streaming Algorithm for Robust Mean Estimation). Let co be sufficiently small positive constant. There is algorithm that, given parameters с E (0,co). c C (0, 1), 5 E (0, 1), reads a stream of poly(d/e) points from an c-corrupted version Nf Ia) (cf. Definition 1.1), computes an estimate ji such that [| - 14/22 = O(e/c) with probability at least 1 - 8. Moreover, the algorithm runs in time (nd + 1/e2+t (polylo)) (d/c)), and uses additional memory O(d+ + poly((/c).
We also restrict our attention to Algorithm 1, which is the first stage of the overall algorithm. This is because the second stage, Algorithm 4, is only utilized when the dimension has been reduced to just (od/loo////)),..)) thus one can simulate the offline algorithm in the streaming setting by storing the whole dataset, which has size ([pooylog(d/e)({) (by "size" we mean the number of scalar numbers that the dataset consists of).
Main Differences from [DKPP22] For the streaming setting, we change our viewpoint and instead of denoting by P (1 e)G + eB the uniform distribution over fixed dataset, we now denote by P the underlying data distribution itself (in our case G = N(,, I) or more generally any distribution satisfying the goodness conditions of Definition 2. I). Consequently, the quantities E, B, M, Σ, : B, M, mentioned in the pseudocode and in our proofs, are all population-Ievel quantities. Instead of having direct access to them (as in Section 3), the algorithm is free to the stream of fresh i.i.d. samples from P and build estimators to approximate them. The estimators needed for this are the same as in [DKPP22], for which we will give overview later on. Apart from this difference, the algorithm from Section is easy to implement in the streaming setting with low memory: The operation of multiplying a vector 2 by an empirical covariance matrix (x (z 14) T (which was used frequently in Section 3) can be implemented in the streaming setting without storing the entire covariance matrix, if one calculates first the inner product (z-))) in the sum and then multiply that scalar result with (x 1--- 4u). Second, the weight function wr (x(2) that the algorithm maintains for re-weighting points is now defined over the entire Rd (and not just over fixed dataset, since there is no such dataset anymore). The algorithm can store a representation of we by storing the matrices U, and how many times Line of the filter has been executed (having these in memory then given any point x E IRd the algorithm can calculate wx(x)) Since the number of iterations (max
and the size of each U, is Od (olylog(d/e)). the memory needed to store the weights is (ld polylolog( overall.
hopooylog(d//
The proof of correctness of the resulting streaming algorithm follows the same arguments before, but needs to account extra errors for the approximations taking place. Instead of giving a complete formal proof, state the new deterministic conditions that are needed in Condition E.3. Then, discuss why each of them is needed and how it can be obtained.
Condition E.3 (Deterministic Conditions for Algorithm 6). Let Scover be the cover set of [DKPP22, Lemma 4.9] For all t [fmua]], the following hold:
(i) Scores: For every x E Scover. it holds |M 1m) |.(x - m. (
(ii) T = I For every w Rd + [0, 1], the algorithm has access to an estimator where xx) = ||U,(z - (4 (I
f(w) for the quantity E. (x-r(w(X))((C
He]1 10, U.)) such that f(w)= > T2/2 whenever (x-) > T..
(iii) Spectral norm of B, X E [0.1] Ilop- B/ llop].
(iv) Frobenius norm U [0.8k] M || 1.2k M, [[].
(v) Maximum norm: ma || |13 <
log k tr (U, U,).
(vi) Probability estimate: 1g: -el S 0.11 tmax))
37

---

(vii) U, almost orthogonal: For every pair uj új of distinct rows of U, it holds K2,, uzحو -
[ 1/ll/yyl////
(viii) Approximate eigenvector: The vector in from Line 19 satisfies u, B, uu 2 (1
1/p)||Bt
Condition This condition used instead Condition (iii) in Condition C.I to go from Equa- tion (20) to Equation (21). The issue with the old condition condition is that in the streaming setting there no fixed dataset for which we can require the approximation to hold. Instead, can the technique from [DKPP22] and define appropriate cover Seover for which Ex-B (X)
1 (see Section 4.2.1 of [DKPP22] for more details the cover argument). Finally, we remark that the additional term [M, E is due to the fact that the matrix U, is calculated Line 17 based on approximations of M, and not M, itself (see [DKPP22, Lemma 4.15)). That error term be tolerated in the steps of the correctness proof that we provided in Appendix C.
اا) as Ex~s(wt+4 [
Condition (ii) Since the stopping condition of the filtering procedure (Algorithm 3) depends population-Ieved quantities, it needs to be estimated within sufficient accuracy. The estimator for that and the correctness of the resulting filter has been analyzed in [DKPP22].
Condition (iii) An approximate top eigenvector for Condition (iii) is obtained by multiplying a random Gaussian vector by log d fresh comple-baseee estimators of B; (i.e., multiply (0, by Mr, where M,- as Line 14). This estimator has been analyzed in [DKPP22, page 24).
Conditions (v) and (vi) These follow similarly to the corresponding conditions of Appendix C. For Condition (vi), the estimator now will use fresh estimates M2 of the form of Line 14, instead of the population-level M,
Condition (viii) This condition can be obtained by [DKPP23, Lemma 4.3] and [DKPP22, Lemma 4.15].
38

---

Algorithm 6 Robust Mean Estimation Under Huber Contamination In Streaming Model (Stage 1)
Input: Parameter e (0, 1/2), stream of i.i.d. points from the distribution P (1 c)G + cB,
where G' satisfies Definition 2.1 with appropriate parameters.
Output: An approximation the mean a subspace V and the orthogonal subspace V.
Let C' be a sufficiently large constant, & C'log" (n + d), mmax (log(d/c))
2: Initialize V1 < 9 and wi (x) for all RRd.
3: for t = 1,..., tmax do
4: Let V4 be the subspace spanned by the vectors in V,, and Vi- be the perpendicular subspace.
5: Let Pt be the distribution /' re-we re-weighted by wy. i.e., Pi (x))) = w; (,)x) Ex- [wx())]
6: Let μF; Σj be the mean and covariance of Projy, (X) when X ~ Pt.
7:
Define B, (Ex~rlwy 2 - (1 C where Iy is the orthogonal pro- jection matrix for V and M, = (B, ) for p log(d). Calculate N such that \ | llop [0.1/10] using power iteration.
8:
& cf. Appendix B
If X, < Ce then return fe, and Ve,
9:
10:
Let q: == Pr..:"~N(o.j) MM MM > 2 x 122/1²)
11:
Calculate an estimate Te such that 17: 4el S tmox
12:
if q - 1/(k2t tmax) then
Case (cf. Section 3.1)
13:
for j E [K] do
Let B B,yee for l E [log d] sample-based versions of B
Let M = - log B. b...e*
Ut Mr Zt, for 21, NO,,
17:
U, < [vr.I...., vt,&] i.e. the matrix with rows ut for j E [k].
18:
W/+11 MULTT
w, U.)
cf. (DKPP22, Algorithm 4] D Case 2 (cf. Section 3.2) Power iteration
19:
else
Find a vector ue such that u, B 2 (1 1/P) llop
21:
V4+1 + V, U }
22:
Let μty = Ex~p, [Projy, (X)) be the mean of Pt after projection to V4.
23: return /,y, and V.,
39


Reviews:
[{"review_id": "douyaN5NKz", "reviewer": "Reviewer_rxyU", "summary": "This paper solves an important problem in robust mean estimation and robust linear regression, that is, can we achieve optimal error rate with nearly linear sample size and runtime when the inliers are standard Gaussian. All the previous works that have linear runtime suffer from an additional $\\sqrt{\\log(1/\\epsilon)}$ factor. To derive such results, they develop a new analysis for the classic filtering technique. Moreover, they develop a novel reduction for the Gaussian linear regression to mean estimation problem, which is also of independent interest.", "questions": "N/A", "limitations": "N/A", "rating": 5, "confidence": 3, "soundness": 4, "presentation": 2, "contribution": 3, "strengths": "- This study solves a long-standing problem in the robust mean estimation and robust linear regression, that is, can we close the gap between the information-theoretic lower bound and the guarantees for algorithms that enjoy linear runtime. This problem, though is nearly the simplest setting, is very fundamental in this field, and I feel like the analysis developed in this paper can be adapted to other tasks.\n\n- The proof technique is very strong. The reduction for Gaussian linear regression is quite novel.\n\n- I like the intuition illustrated in Section 1.2.", "weaknesses": "- In most of the literature listed in this paper, people consider the stronger outlier model, i.e., the strong contamination model. My question is, why the authors can only handle Huber's model? What is the main difficulty in extending to the strong contamination model? \n\n- Some notations are not consistent: the pdf function is denoted by both $P(x)$ (line 154) and $p(x)$ (line 173).\n\n- It would be great if there are some preliminary simulations to support the theory, like plotting the runtime v.s. dimension. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their effort and feedback. \n\nRegarding the summary of the paper, we would like to highlight that the algorithm for linear regression is in fact the *first polynomial algorithm* obtaining error $O(\\sigma \\epsilon)$, i.e., not only attains nearly-linear runtime and nearly-optimal sample complexity, but prior to this, no polynomial-time and sample complexity algorithm was known to achieve the $O(\\sigma \\epsilon)$ error (cf. lines 82-83).\n\nWe now proceed by addressing the reviewer’s questions below:\n\n**(Huber’s vs strong contamination)** There is an inherent difference in the guarantees that can be obtained between the two contamination models. While the information theoretic optimal error is $\\Theta(\\epsilon)$ for both contamination models (and achievable by exponential time algorithms with linear sample complexity), there is a statistical query lower bound providing evidence that it is computationally hard to reach error $o(\\epsilon\\sqrt{\\log(1/\\epsilon)})$ in the strong contamination model (or even total variation model); see lines 146-148 in the introduction. This computational lower bound necessitates that we consider the simpler Huber contamination.\n\n**(Notation)** Throughout the paper, we have consistently used capital letters (e.g., $P(x)$) to denote probability density functions (pdfs). The instance $p(x)$ mentioned in line 173 does not represent a pdf; instead, it denotes the polynomial defined in the previous line.\n\n**(Simulations)** While we acknowledge the reviewer's interest in experiments and the importance of developing practical algorithms, the primary contribution of our work is to complete the computational and statistical landscape for the problem in terms of error guarantee, sample complexity and runtime. As NeurIPS invites theoretical ML papers (such as ours),  we believe that our work should not be negatively evaluated for the lack of experiments, especially when no other nearly-linear time estimator with optimal guarantees was previously known. \n\n\nIn the light of these responses, we kindly request the reviewer to reconsider their score. \n"}, {"review_id": "ngz4kzOROj", "reviewer": "Reviewer_bgYs", "summary": "This paper presents an algorithm to estimate the mean of a Huber contaminated d-dimensional Gaussian with sample complexity of n=O(d/\\epsilon^2)  and a run time that is near-linear. They also provide an algorithm for robust linear regression with sample complexity of n=O(d/\\epsilon^2) and a run time that is again near-linear.  The results improve upon the existing poly(d/\\epsilon) algorithms that achieve \\Theta(\\epsilon) error. ", "questions": "In Appendix B (which forms the core of the proof for theorem 1.3), what is the core insight that helps obtain the final sample complexity? The paper claims that combining the techniques from [DKPP22] and [Dia+18] is highly non-trivial. However, The core proofs in sections B.4.1 do not contain any novel proof techniques apart from the ones already present in [DK19] and [Dia+16]. If you could offer more insight into this in the main paper, that would be good. ", "limitations": "This is a theoretical paper with no negative societal impact. ", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is well-written and the flow is good. Although the paper is highly theoretical, the writing is friendly. \n\nThe results of the paper seem to answer an important open question and the paper claims to achieve a sample complexity that is close the the mean estimation problem in the uncontaminated case, which is good. ", "weaknesses": "The paper seems to combine the techniques present in two existing papers to obtain their result. Although the authors clarify that such a combination is highly non-trivial, the math used in the proof section is very dense for a thorough verification of this claim.\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the positive feedback on the contributions and the writing. \n\nHowever, it appears that the reviewer touched upon and subtly intertwined two distinct aspects of the algorithm: \n(i) combining the algorithmic techniques used in [DKPP22, Dia+18] that is developed in Section 3, and\n(ii) the sample complexity for ensuring the deterministic conditions (discussed in Appendix B.4.1 and its relation to [DK19, Dia+16]).\n Specifically, we would like to highlight that (i) is the main contribution of our work and forms the core of Theorem 1.3. Thus, to minimize any potential confusion, we would like to provide clarifications on both of these points below.\n\n**(Algorithmic techniques from [DKPP22, Dia+18])** Our main contribution is achieving the optimal error in nearly-linear runtime, which requires a non-trivial combination of the algorithm components from [DKPP22, Dia+18]. Combining the techniques from [DKPP22, Dia+18] is the main focus of our paper which we argue that is non-trivial, both from an algorithmic and an analytical viewpoint. We highlight the roadblocks in combining these techniques in Lines 116-121. Our solution is then outlined in Lines 122-134, with the details presented in Section 3. We stress again that this solution is the main contribution of our paper (and orthogonal to [Dia+16,DK19]).\n\n**(Sample complexity for deterministic conditions)** The proposed algorithm above succeeds under a set of (novel) deterministic conditions (given in Definition 2.1). The last of these conditions (Condition 2.c), which is the key to obtaining the optimal error, is novel, and thus requires a new proof, provided in Appendix B.4.1 (only the conditions 2.a and 2.b come from [DK19, Dia+16]). In addition, we show that Condition 2.c holds with linear sample complexity, which should be contrasted with the (large) polynomial sample complexity in [Dia+18] for the corresponding condition there. That said, we stress again that this is not the main result/focus of our paper (and is thus deferred to Appendix). The reviewer is also asking about the core insight for linear sample complexity of Definition 2.1, which we summarize next:\n\n**(Core insight for linear sample complexity of Definition 2.1)** We focus our attention to Condition 2.c, which is the novel deterministic condition, and provide a very high level intuition for why its sample complexity should be (nearly) linear. Observe that this condition asks for a uniform concentration along all matrices $U$ of dimension $d \\times k$, which are objects of $dk$ parameters. Since $k$ is logarithmic in $d$ and $1/\\epsilon$, the total number of parameters is nearly linear in $d$. Informally, as the number of parameters of the objects that we ask uniform concentration for is roughly related to the sample complexity, this eventually results in a nearly linear sample complexity. (More formally, we construct a finite cover for these matrices whose size is exponential with the exponent being nearly linear. Since Gaussians have exponential concentration, the resulting sample complexity is logarithmic in the cover size, and thus nearly linear.)\n"}, {"review_id": "Sk6XnejplW", "reviewer": "Reviewer_BEmV", "summary": "Gaussian estimation and linear regression are fundamental statistical tasks, and the Huber contamination model is one of the most well-studied models for robust statistics. This paper works to design algorithms for these tasks with near-optimal sample complexity, error, and almost linear running time. The main result is for Gaussian robust mean estimation with near-optimal sample complexity of $\\tilde{O}(d/\\varepsilon^2)$, almost linear running time, and near-optimal $\\ell_2$ error of $O(\\varepsilon)$. (Previously there have only been algorithms with (i) near-optimal error but polynomially-suboptimal sample complexity and running time, or (ii) near-linear runtime but suboptimal error.) The paper obtains similar guarantees for robust Gaussian linear regression with Huber contamination: resolving the open problem of whether there exists an algorithm with $O(\\sigma \\varepsilon)$ error in polynomial time and sample complexity. Applicability to streaming is also discussed.", "questions": "Could you provide discussion on implications for more general covariances?\n\n\"reminder\" -> \"remainder\" line 233", "limitations": "The limitations are appropriately discussed.", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "The strengths will read similarly to the summary, in that the obtained results are for fundamental problems that mostly speak for themselves. The tasks of Gaussian mean estimation and linear regression with Huber contamination are core to algorithmic robust statistics, and it is highly desirable to have algorithms with near-optimal error, sample complexity, and near-linear running time simultaneously. \n\nThe main result of this paper is for Gaussian mean estimation with Huber contamination. Grossly oversimplifying, the result works by leveraging ideas related to [DKPP22] and [Dia+18]. [DKPP22] obtains near-linear time but suboptimal error guarantees. It obtains a fast running time by techniques with which it filters in random directions of large variance. [Dia+18] obtains near-optimal error with slower running time by employing a stronger filtering technique that deterministically filters with respect to the subspace for the k-largest eigenvalues. These techniques are not obvious how to use together, as trying to sample for random, strong filtering directions may just consistently produce the same directions, and deterministic filtering will naturally incur a slowdown. In summary, this work’s primary technical contribution is showing how to obtain the guarantees of strong filtering techniques without incurring the slower runtime that is naively associated with them.", "weaknesses": "One could argue that the very high-level approach for the main result of mean estimation is not too surprising in how it blends the intuitions of [Dia+18] and [DKPP22]. Still, there are significant technical obstacles to this approach and it is perhaps even nicer to see how techniques from this general body of work are flexible enough to coalesce into a new, important result.\n\nProviding some discussion on the relationship to settings with non-identity covariance might be beneficial.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the positive evaluation. We respond to their question about general covariance matrices below.\nFirst, we can restrict to the case that the covariance matrix is unknown to the algorithm; otherwise, using an appropriate whitening of the data can reduce to the identity covariance case.  For this case, we note the following:\n\n**(Information-theoretic Error)** For unknown covariance Gaussians, the information theoretic optimal for the Euclidean error is $\\Theta(\\epsilon \\sqrt{\\lVert \\Sigma \\rVert_{op}} )$. Moreover, this can be achieved with linear in $d$ sample complexity by (computationally-inefficient) algorithms. See, e.g., [DK23, Lemma 1.9 and Proposition 1.20]. \n\n\n**(Computationally-Efficient Error)** However, any computationally-efficient (statistical query) algorithm that achieves error $o( \\sqrt{\\epsilon} \\sqrt{\\lVert \\Sigma \\rVert_{op}} )$ requires $\\Omega(d^2)$ samples [Dia+22a, Theorem 6.11 with $k = \\sqrt{d}$]. With $d^2$ samples, one can robustly learn the covariance (without the extraneous $\\log(1/\\epsilon)$ factor) using [Dia+18] and reduce it to the known covariance setting. Crucially, the algorithm in [Dia+18] runs in quasipolynomial time and speeding up this algorithm is an important problem but beyond the scope of our work."}, {"review_id": "Oy3aemhHDF", "reviewer": "Reviewer_qzWP", "summary": "The authors study here a classical problem in statistics: the estimation of Gaussian means and linear regression with Gaussian covariates in the presence of Huber contamination.\n\nThis paper proposes a novel polylog time algorithm to compute the mean of a Gaussian in an $\\epsilon$-corrupted, high-dimensional setting. This algorithm is obtained as a concatenation of two algorithms presented in the literature. They apply this algorithm in Gaussian linear regression and analyze this setting.", "questions": "I believe actual experimentation is lacking: Is it possible to test the algorithm and see its performance for some distribution choices? While it is clear that the distribution is Gaussian, are there some constraints on the other distribution?\n", "limitations": "The paper is theoretical is nature, so there are no limit", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The paper seems technically sound and addresses some relevant questions in the field. The approach is based on an algorithm previously known in the literature [Dia+18], enhancing it with ideas from [DKPP22]", "weaknesses": "While the paper provides a solid theoretical analysis of the time and sample complexity for these algorithms to work, it could be more concerned with seeing it applied. A section could have been dedicated to applying the algorithm in some cases where the parameter values are relevant to some applications.\n\nAdditionally, the work is focused on proving that the algorithm performs as claimed, but it needs to address the limitation of the analysis of the algorithm. The conditions and assumptions are not clearly stated in the introduction of the problem.\n\nThe first part of the paper is straightforward and gives a clear introduction to understanding the questions addressed by the paper. It is pedagogical to explain the position in the current literature. The technical part of explaining the algorithm (specifically Section 3.2) takes more work to follow.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their time and feedback. We respond to the individual points below:\n\n* **(Simulations)** We refer to the response given to reviewer rxyU regarding the same matter.\n* **(``While it is clear that the distribution is Gaussian, are there some constraints on the other distribution?’’)** We are not entirely sure which is the other distribution that the reviewer is referring to. If the reviewer asks if there are any constraints on the distribution of the outliers, then the answer is no: The outliers could be sampled from any (arbitrary) distribution. If the the reviewer is asking whether the inlier distribution should necessarily be Gaussian, then we provide the response below:\n    * **(Distributional choices of inliers)** Regarding testing the algorithm for different distribution choices of inliers, we would like to point out that the Gaussianity assumption is somewhat critical for achieving error $O(\\epsilon)$. Even for univariate sub-Gaussian distributions, the adversary can corrupt the distribution (in the Huber contamination model) such that every algorithm incurs error $\\Omega(\\epsilon \\sqrt{\\log(1/\\epsilon)})$; see [DK19, page 4]. \n    * **(Possible extensions)**  As outlined above, it is impossible to relax the Gaussianity assumption to general sub-gaussian distributions for inliers. Nonetheless, we can extend our theoretical results to inliers following subgaussian distributions which are centrally symmetric (to ensure correctness of median) and satisfy Hanson-Wright inequality (to ensure that multi-directional filter works).\n\n* **(``The conditions and assumptions are not clearly stated in the introduction of the problem’’)** The introduction provides an informal overview of the problem, and Section 1.1. provides complete formal statements for the two results with all conditions and assumptions stated. We are not sure which are the missing conditions that the reviewer is referring to.\n\n* **(``Position in the current literature’’)** This is mentioned in the related work section; cf. Lines 144-151. Given the extra space in the camera ready version, we could mention it earlier in the paper as well. In summary, the paper is positioned as follows: \n    * For subgaussian distributions, there are existing nearly-linear runtime robust mean estimation algorithms getting the information-theoretic optimal error of $O(\\epsilon \\sqrt{\\log(1/\\epsilon)})$ in the strong contamination model.\n    * For Gaussian distributions, these algorithms obviously continue to work, but the information-theoretic optimal error is now $O(\\epsilon)$ and existing Statistical Query lower bounds show that achieving that error is computationally hard under the strong (or even TV-distance) contamination model. Thus, one needs to restrict to the Huber contamination model.\n    *  For Gaussians under Huber contamination, [Dia+18] gave the first polynomial time-and-sample algorithm achieving the optimal $O(\\epsilon)$ error. In the present paper, we improve the runtime and sample complexity to be nearly linear. Additionally, we also obtain similar guarantees for robust linear-regression, where prior to our results, not even polynomial algorithms were known for $O(\\sigma \\epsilon)$ error.\n"}, {"review_id": "eX5m4s2kMI", "reviewer": "Reviewer_7D1F", "summary": "In this paper, the authors give an improved algorithm for estimating the mean of a multivarite Gaussian random variable under the Huber contamination model. Diakonikolas et al. Have given the first algorithm in this model that achieves the information theoretical optimal $\\ell_2$ loss of $\\Theta(\\varepsilon)$ using a polynomial number of samples when the noise rate is $\\varepsilon$. Pensia et al. have given a fast streaming algorithm for the same problem with an $\\ell_2$ loss of $\\Theta(\\varepsilon\\log \\varepsilon^{-1})$ with a sample complexity of $d^2\\varepsilon^{-2}$ where $d$ is the dimension. The authors carefully combine the ideas in the above two papers to achieve the sample complexity of $d\\varepsilon^{-1}$ which is optimal. \n\nAt a high level, a mean vector $\\mu’$ and a subspace $V$ is first found such that $\\mu’$ and $\\mu$ are close in $\\ell_2$ norm in $V^{\\bot}$. Moreover, $V$ has dimension at most $\\mathrm{\\poly}\\log(d\\varepsilon^{-1})$ such that a high sample-complexity inspired from DIakonikolas et al. is affordable. The finding of this subspace $V$ uses ideas from the Pensia et al. Paper by looking at a random direction where the eigenvalue is large as opposed to looking at the top k eigenvalues as in the Diakonikolas paper which incurse the worse loss. \n\nThe paper also gives a linear-sample algorithm for robust Gaussian linear regression by reducing it to the robust mean estimation problem of a certain conditional distribution of the multivariate Gaussian. This algorithm is slightly surprising to me as any continuous distribution has zero probability of exactly seeing a particular item in the sample space and therefore conditional sampling in infeasible. However, the authors manage to achieve their reduction by carefully considering a large enough ball that allows their reduction to go through. \n\nI have not carefully checked the mathematical details of the paper.\n \n\n\n ", "questions": "A question to the authors: all the papers in this line of work achieves a recovery rate of say $\\varepsilon/1000$ which looks quite lossy to me. Univarite robust estimation achives close to $\\varepsilon$ rate, or smaller constants as do several other papers in statistics such as: [arxiv:2002.01432]. So, my question to the authors is: is there any hope that the constant can be improved, or this filtering technique has some inherent barrier?  \n", "limitations": "None.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The robust statistics in general and the multivariate robust mean estimation, in particular, have received wide attention recently. Therefore, the improved sample-optimal algorithm derived in this paper should be very interesting to the community.", "weaknesses": "None.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their positive comments and feedback. \n\n**(Constants in the breakdown point)** We would like to point out that constants could indeed be optimized, but that work falls outside of the primary aim of our paper. The filtering techniques that we use and their variants have been used to get better constants and breakdown points. In particular, [Dia+18] already takes special care to improve the constant in the error guarantee. Moreover, [arxiv:2002.01432] in its core uses similar filtering techniques to obtain better breakdown points. Finally, optimal break-down point can be achieved using filtering-based algorithms, see, for example, [ZJS22, , Exercise 2.10 of DK23, Appendix E of HLZ20]. As the primary focus on this paper is demonstrating that the nearly-linear runtime is possible, we did not optimize constants.\n"}]

Author rebuttal:
We thank the reviewers for their time and effort in providing feedback. We are encouraged by the positive comments, and that all the reviewers appreciated the paper for the following: (i) **novelty** (rxyU, BEmV, qzWP), **importance** (rxyU,bgYs,BEmV, 7D1F) (iii) **good** and **friendly** writing (bgYs), and (iv) **surprising results** (7D1F). We address the individual questions and comments by the reviewers separately.


Meta-review: The paper tackles the fundamental problem of Gaussian mean estimation under Huber contamination. The authors propose the first method that linear time method that is sample optimal.

The reviewers were positive about this paper. While there is some variance between the reviewer scores, some of the reasoning for the lower scores (such as lack of experiments for this theory paper) was not defended nor elaborated in the discussion period.

Given the importance of the contributions and the overall scores, I recommend acceptance.
Decision: Accept (poster)