High-Dimensional Robust Mean Estimation in Nearly-Linear Time

Introduction

In this paper, we study the robust (or agnostic) setting when a constant $\epsilon<1/2$ fraction of our samples can be adversarially corrupted. We consider the following model of robust estimation (see, e.g., [DKK+16]) that generalizes other existing models, including Huber’s contamination model [Hub64]:

In the context of robust mean estimation studied in this paper, the goal is to output a hypothesis vector $\widehat{\mu}$ such that $\|\widehat{\mu}-\mu^{\star}\|_{2}$ is as small as possible. How do we estimate $\mu^{\star}$ in this regime? A moment’s thought reveals that the empirical mean inherently fails in the robust setting: even a single corrupted sample can arbitrarily compromise its performance. However, one can construct more sophisticated estimators that are provably robust. The information-theoretically optimal error for robustly estimating the mean of $\mathcal{N}(\mu^{\star},I)$ is $\Theta(\epsilon+\sqrt{d/N})$ [Tuk75, DG92, CGR15]. That is, when there are enough samples ($N=\Omega(d/\epsilon^{2})$) one can estimate the mean to accuracy $\Theta(\epsilon)$. Under different assumptions on the distribution of the good data, the optimal error guarantee may be different as well (see Section 1.2). However, the standard robust estimators (e.g., Tukey’s median [Tuk75]) require exponential time in the dimension $d$ to compute. On the other hand, a number of natural approaches (e.g., naive outlier removal, coordinate-wise median, geometric median, etc.) can only guarantee error $\Omega(\epsilon\sqrt{d})$ (see, e.g., [DKK+16, LRV16]), even in the infinite sample regime. That is, the performance of these estimators degrades polynomially with the dimension $d$, which is clearly unacceptable in high dimensions.

Recent work [DKK+16, LRV16] gave the first polynomial time robust estimators for a range of high-dimensional statistical tasks, including mean and covariance estimation. Specifically, [DKK+16] obtained the first robust estimators for the mean with dimension-independent error guarantees, i.e., whose error only depends on the fraction of corrupted samples $\epsilon$ but not on the dimensionality of the data. Since the dissemination of [DKK+16, LRV16], there has been a substantial number of subsequent works obtaining robust learning algorithms for a variety of unsupervised and supervised high-dimensional models. (See Section 1.3 for a summary of related work.)

Although the aforementioned works gave polynomial time robust learning algorithms for several fundamental learning tasks, these algorithms are at least a factor $d$ slower than their non-robust counterparts (e.g., the sample average for the case of mean estimation), hence are significantly slower in high dimensions. It is an important goal to design robust learning algorithms with near-optimal sample complexity that are also nearly as efficient as their non-robust counterparts. In particular, we propose the following broad question:

Can we design (nearly-)sample optimal robust learning algorithms — with dimension independent error guarantees — that run in nearly-linear time?

Here by nearly-linear time, we mean that the runtime is proportional to the size of the input, within poly-logarithmic in the input size and $\operatorname{poly}(1/\epsilon)$ factors. In addition to its potential practical implications, we believe that understanding the above question is of fundamental theoretical interest as it can elucidate the effect of the robustness requirement on the computational complexity of high-dimensional statistical learning/estimation.

For example, for the prototypical problem of robustly estimating the mean of a high-dimensional distribution, previous robust algorithms [DKK+16, LRV16, SCV18] have runtime at least $\Omega(Nd^{2})$ for constant $\epsilon$. Since the input size is $\Theta(Nd)$, we would like to obtain algorithms that run in time $\widetilde{O}(Nd)/\operatorname{poly}(\epsilon)$, where the $\widetilde{O}(\cdot)$ notation hides logarithmic factors in its argument. As the main contribution of this paper, we obtain such algorithms under different assumptions about the distribution of the good data. Our algorithms have optimal sample complexity, provide the information-theoretically optimal accuracy, and — importantly — run in time $\widetilde{O}(Nd)/\operatorname{poly}(\epsilon)$.

2 Our Results

It is well-known (see, e.g., [DKK+17]) that the optimal error guarantee under the assumptions of Theorem 1.2 is $\Omega(\epsilon\sqrt{\log 1/\epsilon})$, even in the infinite sample regime. Moreover, the sample complexity of the learning problem is known to be $\Omega(d/\delta^{2})$ even without corruptions. Thus, our algorithm has best possible error guarantee and sample complexity, up to constant factors. Prior work [DKK+16, DKK+17] gave algorithms with the same error and sample complexity guarantees, but with runtime $\Omega(Nd^{2})$, even for constant $\epsilon$. We note that for the very special case that $D=\mathcal{N}(\mu^{\star},I)$, an error of $O(\epsilon)$ is information-theoretically possible. However, as shown in [DKS17], any Statistical Query algorithm that runs in time $\operatorname{poly}(N)$ needs to have error $\Omega(\epsilon\sqrt{\log(1/\epsilon)})$. Our algorithm achieves this accuracy guarantee in nearly-linear time. See Section 1.3 for a detailed summary of previous work.

Theorem 1.2 handles the case that the covariance matrix of the good data distribution is known a priori. This is a somewhat limiting assumption. In our second main algorithmic result, we obtain a similarly robust algorithm under the much weaker assumption that the covariance matrix is unknown and bounded from above. Specifically, we show:

Similarly, the sample complexity of our algorithm is best possible within a logarithmic factor, even without corruptions; the $O(\sigma\sqrt{\epsilon})$ error guarantee is known to be information-theoretically optimal, up to constants, even in the infinite sample regime. Previous algorithms [DKK+17, SCV18] gave the same sample complexity and error guarantees, but again with significantly higher time complexities in high dimensions. Specifically, the iterative spectral algorithm of [DKK+17] has runtime $\Omega(N^{2}\cdot d)=\Omega(d^{3}/\epsilon^{2})$. See Section 1.3 for more detailed comparisons.

We note that an efficient algorithm for robust mean estimation under bounded covariance assumptions has been recently used as a subroutine [PSBR18, DKK+18b] to obtain robust learners for a wide range of supervised learning problems that can be phrased as stochastic convex programs. This includes linear and logistic regression, generalized linear models, SVMs (learning linear separators under hinge loss), and many others. The algorithm of Theorem 1.3 provides a faster implementation of such a subroutine, hence yields faster robust algorithms for all these problems.

3 Related and Prior Work

Learning in the presence of outliers is an important goal in statistics and has been studied in the robust statistics community since the 1960s [Hub64]. After several decades of work, a number of sample-efficient and robust estimators have been discovered (see [HR09, HRRS86] for book-length introductions). For example, the Tukey median [Tuk75] is a sample-efficient robust mean estimator for various symmetric distributions [DG92, CGR15]. However, it is NP-hard to compute in general [JP78, AK95] and the many heuristics for computing it degrade in the quality of their approximation as the dimension scales [CEM+93, Cha04, MS10].

Until recently, all known computationally efficient high-dimensional estimators could only tolerate a negligible fraction of outliers, even for the simplest statistical task of mean estimation. Recent work in the theoretical computer science community [DKK+16, LRV16] gave the first efficient robust estimators for basic high-dimensional unsupervised tasks, including mean and covariance estimation. Since the dissemination of [DKK+16, LRV16], there has been a flurry of research activity on robust learning algorithms in both supervised and unsupervised settings [BDLS17, CSV17, DKK+17, DKS17, DKK+18a, SCV18, DKS18b, DKS18a, HL18, KSS18, PSBR18, DKK+18b, KKM18, DKS19, LSLC18, CDKS18].

For the specific task of robust mean estimation, [DKK+16] designs two related algorithmic techniques with similar sample complexities and error guarantees: a convex programming method and an iterative spectral outlier removal method (filtering). The former method inherently relies on the ellipsoid algorithm (leading to polynomial, yet impractical, runtimes), while the latter only requires repeated applications of power iteration to compute the highest eigenvalue-eigenvector of a covariance-like matrix. The total number of power iteration calls can be as large as $\Omega(d)$, for constant $\epsilon$, leading to runtimes of the form $\widetilde{\Omega}(Nd^{2})$. We note that the filter-based robust mean estimation algorithm, as presented in [DKK+16], applies to the sub-gaussian case (as in Theorem 1.2). A slight variant of the method [DKK+17] applies under second moment assumptions (as in Theorem 1.3).

The work [LRV16] gives a recursive dimension-halving technique with near-optimal accuracy, up to a logarithmic factor in the dimension. The aforementioned method requires computing the SVD of a second moment matrix $\Omega(\log d)$ times. Consequently, each iteration incurs runtime $\Omega(d^{3})$. Similarly, the robust mean estimation algorithm under bounded second moments in [SCV18] requires computing the SVD of a matrix multiple times, leading to $\Omega(d^{3})$ runtime.

4 Our Approach and Techniques

In this section, we provide a detailed outline of our algorithmic approach in tandem with a brief comparison to the most technically relevant prior work. To robustly estimate the unknown mean $\mu^{\star}$, we proceed as follows: Starting with an initial guess $\nu$, in a sequence of iterations we either certify that the current guess is close to the true mean $\mu^{\star}$ or refine our current guess with a new one that is provably closer to $\mu^{\star}$.

Our approach will try to minimize the weighted second order moment $\sum_{i=1}^{N}w_{i}(X_{i}-\nu)(X_{i}-\nu)^{\top}$ for all $w\in\Delta_{N,\epsilon}$, with the intended solution being assigning $1/|G|$ weight to all the good samples. This can be formalized as an SDP:

This SDP is similar to the convex program used in [DKK+16] but has some important conceptual differences that allow us to get a faster algorithm. The convex program in [DKK+16] is essentially this SDP with $\nu=\mu^{\star}$. However, of course one cannot solve it directly as we do not know $\mu^{\star}$. To overcome this difficulty, [DKK+16] designs a separation oracle, which roughly corresponds to finding a direction of large variance. The whole convex programming algorithm in [DKK+16] then relies on the ellipsoid algorithm and is therefore slow in high dimensions.

In contrast, we fix a guess $\nu$ for the true mean in the SDP. Even though this $\nu$ may not be correct, we will establish a win-win phenomenon: either $\nu$ is a good guess of $\mu^{\star}$ in which case we get a good set of weights, or $\nu$ is far from $\mu^{\star}$ and we can efficiently find a new guess $\nu^{\prime}$ that is closer to $\mu^{\star}$ by a constant factor.

More precisely, we will show that for any guess $\nu$ that is sufficiently close to the actual mean $\mu^{\star}$, the optimal value of the SDP is small. In this case, the weights $\{w_{i}\}_{i=1}^{N}$ computed by the SDP can be used to produce an accurate estimate of the mean: $\widehat{\mu}_{w}=\sum_{i=1}^{n}w_{i}X_{i}$ (see Lemma 3.2). Note that in this case the estimate $\widehat{\mu}_{w}$ will be more accurate than the current guess $\nu$. On the other hand, when the guess $\nu$ is far from $\mu^{\star}$, the optimal value of the SDP is large, and the optimal dual solution gives a certificate on why the second order moment $\sum_{i=1}^{N}w_{i}(X_{i}-\nu)(X_{i}-\nu)^{\top}$ cannot be small no matter how we re-weight the samples using $w\in\Delta_{N,\epsilon}$. Intuitively, the reason that the second moment matrix cannot have small spectral norm is because of the extra component $(\nu-\mu^{\star})(\nu-\mu^{\star})^{\top}$ in the expected second moment matrix. That is, the dual solution gives us information about $\nu-\mu^{\star}$ (Lemma 3.3).

So far, we have sketched our approach of reducing the algorithmic problem to solving a small number of SDPs. To get a fast algorithm, we need to solve the primal and dual SDPs in nearly-linear time. We achieve this by reducing them to covering/packing SDPs and using the solvers in [ALO16, PTZ16]. We note that these solvers rely on the matrix multiplicative weights update method (mirror descent), though we will not use this fact in our analysis. The main technical challenge here is that the approximate solutions to the reduced SDPs may violate some of the original constraints (specifically, the resulting $w$ may not be in $\Delta_{N,\epsilon}$). We show that our main arguments are robust enough to handle these mild violations.

A perhaps surprising byproduct of our results is that a natural family of SDPs leads to asymptotically faster algorithms for robust mean estimation than the previous fastest spectral algorithm [DKK+16] for the most interesting parameter regime (corresponding to large dimension $d$ so that $d\gg\operatorname{poly}(1/\epsilon)$). We view this as an interesting conceptual implication of our results: in our setting, principled SDP formulations can lead to faster runtimes compared to spectral algorithms, by exploiting the additional structure of these SDPs. This phenomenon illustrates the value of obtaining a deeper understanding of such convex formulations.

5 Structure of This Paper

In Section 3, we describe our algorithmic approach for robust mean estimation and use it to obtain our algorithm for sub-gaussian distributions (thus establishing Theorem 1.2). In Section 4, we show that the corresponding SDPs can be solved in nearly-linear time. In Section 5, we adapt our approach from Section 3 to obtain our algorithm for robust mean estimation under bounded covariance assumptions (thus establishing Theorem 1.3). For the clarity of the presentation, some proofs have been deferred to an appendix.

Preliminaries

Throughout this paper, we use $D$ to denote the ground-truth distribution. We use $d$ for the dimension of $D$, $N$ for the number of samples, and $\epsilon$ for the fraction of corrupted samples. Let $G^{\star}$ be the original set of $N$ uncorrupted samples drawn from $D$. After the adversary corrupts an $\epsilon$-fraction of $G^{\star}$, we use $G\subseteq G^{\star}$ to denote the remaining set of good samples, and $B$ to denote the set of bad samples added by the adversary. Note that $G\cup B$ is the input given to the algorithm, and we have $G\subseteq G^{\star}$, $|G|\geq(1-\epsilon)N$, and $|B|\leq\epsilon N$.

Robust Mean Estimation for Known Covariance Sub-Gaussian Distributions

In this section, we will describe our algorithmic technique and give an algorithm establishing Theorem 1.2.

As we described in Section 1.4, our algorithm is going to make a guess $\nu$ for the actual mean $\mu^{\star}$, and try to certify its correctness by an SDP. In Section 3.1, we give the SDP formulation and describe the entire algorithm. In Section 3.2, we show that the optimal value of the primal/dual SDPs are closely related to the distance $\left\|\nu-\mu^{\star}\right\|_{2}$. When the current guess $\nu$ is close to $\mu^{\star}$, we show (Section 3.3) that the solution to the primal SDP is going to give an accurate estimate of $\mu^{\star}$. On the other hand, when the current guess $\nu$ is far, in Section 3.4 we analyze the dual solution and show how to find a new guess $\nu^{\prime}$ that is closer to $\mu^{\star}$. Finally, we combine these techniques and prove Theorem 1.2 in Section 3.5.

Intuitively, this SDP tries to re-weight the samples to minimize the second moment matrix $\sum_{i=1}^{N}w_{i}(X_{i}-\nu)(X_{i}-\nu)^{\top}$. The intended solution to this SDP is to assign weight $1/|G|$ on each of the good samples. This solution will have a small objective value whenever $\nu$ is close to $\mu^{\star}$.

When $\nu$ is far from $\mu^{\star}$, we need to consider the dual of (2). We will first derive the dual of (2). Note that the primal SDP is equivalent to the following:

Strong duality holds in our setting because the primal SDP admits a strictly feasible solution. The dual SDP can now be written as follows:

Observe that once we fix a dual solution $M$, it is easy to minimize the objective function over $w$: the minimum value is attained by assigning weight $w_{i}=\frac{1}{(1-\epsilon)N}$ to the smallest $(1-\epsilon)N$ inner products. Therefore, the dual SDP can be equivalently written as follows:

The dual SDP (3) certifies that there are no good weights that can make the spectral norm small. The intended solution for the dual is $M=yy^{\top}$, where $y=\frac{\nu-\mu^{\star}}{\left\|\nu-\mu^{\star}\right\|_{2}}$ is the direction between $\nu$ and $\mu^{\star}$. Note that when $M=yy^{\top}$, the value $(X_{i}-\nu)^{\top}M(X_{i}-\nu)$ is exactly the squared norm of the projection in the direction $y$. Intuitively, if we project the samples onto the direction of $y$, the mean of the good samples is going to be at distance $\left\|\nu-\mu^{\star}\right\|_{2}$, so even after removing the farthest $\epsilon$-fraction of the projected samples one cannot make the remaining values of $(X_{i}-\nu)^{\top}M(X_{i}-\nu)$ small. Of course, in general, the dual solution can be of rank higher than $1$, but we will show that any near-optimal dual solution must be close to rank $1$ later in Section 3.4.

To avoid dealing with the randomness of the good samples, we require the following deterministic conditions on the original set of $N$ good samples $G^{\star}$ (which hold with probability $1-\tau$) drawn from the sub-gaussian distribution. For all $w\in\Delta_{N,3\epsilon}$, we require the following conditions to hold for $\delta=c_{1}(\epsilon\sqrt{\log 1/\epsilon})$ and $\delta_{2}=c_{1}(\epsilon\log 1/\epsilon)$ for some universal constant $c_{1}$:

Intuitively, Equations (4) show that removing any $\epsilon$-fraction of good samples will not distort the mean and the covariance by too much. Equation (5) says that the good samples are not too far from the true mean.

We note that the above deterministic conditions are identical to the ones used in the convex programming technique of [DKK+16] to robustly learn the mean of $\mathcal{N}(\mu^{\star},I)$. We note that the proof of these concentration inequalities does not require the Gaussian assumption, and it directly applies to sub-Gaussian distributions with identity covariance. It follows from the analysis in [DKK+16] that after $N=\Omega(\delta^{-2}(d+\log(1/\tau)))$ samples, these conditions hold with probability at least $1-\tau$ on the set of good samples.

Throughout the rest of this section, we will assume that the above conditions are satisfied where we set the parameter $\tau$ to be a sufficiently small universal constant; selecting $\tau=1/30$ suffices for all our arguments.

We are now ready to present our algorithm (Algorithm 1) to robustly estimate the mean of known covariance sub-gaussian distributions.

Notation. In this section, we will use $c_{1},\ldots,c_{7}$ to denote universal constants that are independent of $N$, $d$, and $\epsilon$. We will give a detailed description on how to set these constants in Appendix A.

2 Optimal Value of the SDPs

In this subsection, we will give upper and lower bounds on the optimal value of the SDPs (2) and (3). Recall that our high-level idea is to use the dual SDP to improve our guess $\nu$, until it is close enough to the true mean $\mu^{\star}$, and then solve the primal SDP to get a good set of weights. However, we cannot write an if statement based on $r=\left\|\nu-\mu^{\star}\right\|_{2}$ because we do not know $\mu^{\star}$.

In particular, when $r\geq c_{2}\beta$, we can simplify the above as

One feasible primal solution is to set $w_{i}=\frac{1}{|G|}$ for all $i\in G$ (and $w_{i}=0$ for all $i\in B$). Therefore,

Notice that $w$ can be viewed as a weight vector on $G^{\star}$ and we have $w\in\Delta_{N,\epsilon}$. This allows us to use Condition (4) in the second to last step.

One feasible dual solution is $M=yy^{\top}$ where $y=\frac{\mu^{\star}-\nu}{\left\|\mu^{\star}-\nu\right\|_{2}}$. The dual objective value is the mean of the smallest $(1-\epsilon)$-fraction of $\left((X_{i}-\nu)^{\top}M(X_{i}-\nu)\right)_{i=1}^{N}$, which is at least

This is because $|G|=(1-\epsilon)N$, the smallest $(1-\epsilon)N$ entries must include $S$, where $S$ is the smallest $(1-2\epsilon)N$ entries in $G$. Let $w^{\prime}_{i}=\frac{1}{|S|}$ for all $i\in S$ and $w^{\prime}_{i}=0$ otherwise. Note that $S\subset G$ and $|S|=(1-2\epsilon)N$, so $w^{\prime}$ can be viewed as a weight vector on $G^{\star}$ with $w^{\prime}\in\Delta_{N,2\epsilon}$. Therefore we have

3 When Primal SDP Has Good Solutions

In this section, we show that a good primal solution for any guess $\nu$ will give an accurate weighted empirical mean. Lemma 3.2 proves the contrapositive statement: if the weighted empirical mean $\widehat{\mu}_{w}$, with respect to weight-vector $w$, is far from the true mean, then no matter what our current guess $\nu$ is, $w$ cannot be a good solution to the primal SDP. More specifically, we show that the objective value of $w$ is at least $1+\Omega(\delta^{2}/\epsilon)$. Roughly speaking, we get a contribution of $1$ from the good samples and a contribution of $\Omega(\delta^{2}/\epsilon)$ from the bad samples.

We briefly explain why the bad samples contribute $\Omega(\delta^{2}/\epsilon)$. The empirical mean of the good samples is off by at most $\delta$ by Condition (4). Now if $\widehat{\mu}_{w}$ is far away from $\mu^{\star}$, the bad samples must shift the mean by more than $\Omega(\delta)$. Intuitively, if an $\epsilon$-fraction of the samples distort the mean by $\delta$, on average each of these sample contributes an error of $\delta/\epsilon$, which introduces a total error of $\epsilon(\delta/\epsilon)^{2}=\delta^{2}/\epsilon$ in the second moment matrix.

Fix any $w\in\Delta_{N,2\epsilon}$. If $\left\|\mu^{\star}-\nu\right\|_{2}\geq c_{5}\beta$, then because $w$ is feasible and by Lemma 3.1,

Therefore, for the rest of this proof, we can assume $\left\|\mu^{\star}-\nu\right\|_{2}<c_{5}\beta$.

We project the samples along the direction of $(\widehat{\mu}_{w}-\mu^{\star})$. Consider the unit vector $y=(\widehat{\mu}_{w}-\mu^{\star})/\left\|\widehat{\mu}_{w}-\mu^{\star}\right\|_{2}$. To bound from below the maximum eigenvalue, it is sufficient to show that

We first bound from below the contribution of the bad samples by $\Omega(\delta^{2}/\epsilon)$. By triangle inequality,

The last line follows from our choice of $y$, and the good samples satisfy Condition (4). By Cauchy-Schwarz, $\left(\sum_{i\in B}w_{i}\langle X_{i}-\nu,y\rangle^{2}\right)\left(\sum_{i\in B}w_{i}\right)\geq\left(\sum_{i\in B}w_{i}\langle X_{i}-\nu,y\rangle\right)^{2}\geq c_{6}^{2}\delta^{2}.$ Since $w_{B}\leq 2\epsilon$, we have $\sum_{i\in B}w_{i}\langle X_{i}-\nu,y\rangle^{2}\geq\frac{c_{6}^{2}}{2}(\delta^{2}/\epsilon)$.

We continue to lower bound the contribution of the good samples to the quadratic form by $1-O(\delta_{2})=1-O(\delta^{2}/\epsilon)$. This is because the true covariance matrix is $I$. By Condition (4),

Putting the good and bad samples together, we have $\sum_{i=1}^{N}w_{i}\langle X_{i}-\nu,y\rangle^{2}\geq 1-c_{7}\delta_{2}+\frac{c_{6}^{2}}{2}(\delta^{2}/\epsilon)=1+(\frac{c_{1}^{2}c_{6}^{2}}{2}-c_{1}c_{7})\beta^{2}\geq 1+c_{4}\beta^{2}$ as needed.

The constants in the proof are given in Appendix A. ∎

Lemma 3.2 guarantees that any good solution to the primal SDP gives a good set of weights. In other words, whenever we have a solution to the primal SDP whose objective value is at most $1+O(\beta^{2})$, we are done because the weighted empirical mean must be close to the true mean.

4 When Primal SDP Has No Good Solutions

We now deal with the other possibility: the primal SDP has no good solution. We will show that, in this case, we can move $\nu$ closer to $\mu^{\star}$ by solving the dual SDP (3), decreasing $\left\|\nu-\mu^{\star}\right\|_{2}$ by a constant factor.

Because $\operatorname{tr}(M)=1$, we can remove $1$ from both sides and get $\langle M,(\nu-\mu^{\star})(\nu-\mu^{\star})^{\top}\rangle\approx\left\|\nu-\mu^{\star}\right\|_{2}^{2}$. This condition implies that the top eigenvector of $M$ aligns approximately with $(\nu-\mu^{\star})$, which provides a good direction for us to move $\nu$.

The following lemma formalizes this intuition. Specifically, Lemma 3.3 shows that despite the error from solving the SDP approximately and the errors in the concentration inequalities, we can still use the top eigenvector of $M$ to move $\nu$ closer to $\mu^{\star}$.

We know $M\succeq 0$ and $\operatorname{tr}(M)=1$. Without loss of generality, we can assume $M$ is symmetric. Using Condition (4), we can prove that $\langle M,(\mu^{\star}-\nu)(\mu^{\star}-\nu)^{\top}\rangle\geq\frac{3}{4}\left\|\mu^{\star}-\nu\right\|_{2}^{2}$:

We will continue to show that the top eigenvector of $M$ aligns with $(\nu-\mu^{\star})$. Let $\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{d}\geq 0$ denote the eigenvalues of $M$, and let $v_{1},\ldots,v_{d}$ denote the corresponding eigenvectors. The conditions on $M$ implies that $\sum_{i=1}^{d}\lambda_{d}=1$. We decompose $(\mu^{\star}-\nu)$ and write it as $\mu^{\star}-\nu=\sum_{i=1}^{d}\alpha_{i}v_{i}$ where $\sum_{i=1}^{d}\alpha_{i}^{2}=\left\|\mu^{\star}-\nu\right\|_{2}^{2}$. Using these decompositions, we can rewrite $\langle M,(\mu^{\star}-\nu)(\mu^{\star}-\nu)^{\top}\rangle=\sum_{i=1}^{d}\lambda_{i}\alpha_{i}^{2}$.

The constants in the proof are given in Appendix A. ∎

5 Proof of Theorem 1.2

We are now ready to prove Theorem 1.2. This is mostly done by applying Lemmas 3.2 and 3.3 in appropriate scenarios. Because of the geometric improvement in Lemma 3.3, we will apply it at most a logarithmic number of times, and then the algorithm can terminate in the case of Lemma 3.2.

Let $\tau=1/30$. When $N=\Omega(d/\delta^{2})$, Condition (4) holds for the good samples with probability at least $1-\tau$, which is required in the proofs of Lemmas 3.1, 3.2, and 3.3.

We will use the empirical coordinate-wise median as our initial guess $\nu$. It is folklore that with high probability, the coordinate-wise median is within $O(\epsilon\sqrt{d})$ of the true mean $\mu^{\star}$. In Algorithm 1, whenever we update $\nu$ by Lemma 3.3, we move it closer to $\mu^{\star}$. Therefore, throughout the algorithm, the condition $\left\|\nu-\mu^{\star}\right\|_{2}\leq O(\epsilon\sqrt{d})$ always holds, which is required by Proposition 4.1.

Solving Primal/Dual SDPs in Nearly-Linear Time

By combining Lemmas 3.2 and 3.3 from Section 3, we know that we can make progress by either finding any solution to the primal SDP (2) with objective value at most $1+c_{4}\beta^{2}$, or by finding an approximately optimal solution to the dual SDP (3) whose objective value is at least $1+\frac{9}{10}c_{4}\beta^{2}$. This section is dedicated to proving Proposition 4.1, which shows that this can be done in time $\widetilde{O}(Nd)/\operatorname{poly}(\epsilon)$.

Previously, nearly-linear time SDP solvers were developed for packing/covering SDPs [ALO16, PTZ16]. At a high level, we first relate SDPs (2), (3) with a pair of packing/covering SDPs (6), (7), where we switch the objective function with some constraint and introduce an additional parameter $\rho>0$. Next, we show that to prove Proposition 4.1, it is sufficient to solve SDPs (6), (7) approximately for the correct value of $\rho$, and moreover, we can run binary search to find a suitable $\rho$. Finally, in Section 4.1, we show that our packing/covering SDPs (6), (7) can be solved in time $\widetilde{O}(Nd/\epsilon^{6})$. Note that these running times (specifically, the dependence on $\epsilon$) can be improved if better packing/covering SDP solvers are discovered. For example, [ALO16] mentioned the possibility of achieving a bound of $\widetilde{O}(Nd/\epsilon^{5})$ by combining their approach and the techniques from [WMMR15].

Consider the following packing SDP (6) and its dual covering SDP (7) with parameters $(\nu,\epsilon,\rho)$:

We will first show that the solutions of (6) and (7) are closely related to solutions of our original SDPs (2) and (3). Formally, the following lemma shows that if we (approximately) solve the packing/covering SDPs (6), (7) for some value of $\rho>0$ and the resulting objective values are close to $1$, then we can translate these solutions back to obtain solutions for SDPs (2) (3) with objective value roughly $1/\rho$.

We first construct a solution $w$ to SDP (2) with parameters $(\nu,2\epsilon)$ given $w^{\prime}$. Let $w=\frac{w^{\prime}}{\left\|w^{\prime}\right\|_{1}}$. Since $\left\|w^{\prime}\right\|_{1}\geq 1-\frac{\epsilon}{10}$, we know that $w\in\Delta_{N,2\epsilon}$ is feasible for SDP (2). SDP (6) guarantees that $\rho\sum_{i}w^{\prime}_{i}(X_{i}-\nu)(X_{i}-\nu)^{\top}\preceq I$, so the objective value of SDP (2) at $w$ satisfies $\lambda_{\max}\left(\sum_{i}w_{i}(X_{i}-\nu)(X_{i}-\nu)^{\top}\right)\leq\frac{1}{\rho(1-\epsilon/10)}$.

Next, we will construct a solution $M$ for the original dual SDP (3) with parameters $(\nu,\epsilon)$ given $(M^{\prime},y^{\prime})$. We will work with the following SDP that is equivalent to the dual SDP (3):

Let $M=\frac{M^{\prime}}{\operatorname{tr}(M^{\prime})}$, $y=\frac{(1-\epsilon)N}{\rho\operatorname{tr}(M^{\prime})}y^{\prime}$, and $z=\frac{1}{\rho\operatorname{tr}(M^{\prime})}$. Note that $y$ is well-defined, because we always have $M^{\prime}\neq\bf{0}$, otherwise the objective value is at least $1/(1-\epsilon)>1$. Note that $(M,y,z)$ is a feasible solution to SDP (3): by the definition of $(M,y,z)$, the constraint $\rho X_{i}^{\top}M^{\prime}X_{i}+(1-\epsilon)Ny^{\prime}_{i}\geq 1$ translates to $\rho\operatorname{tr}(M^{\prime})X_{i}^{\top}MX_{i}+\rho\operatorname{tr}(M^{\prime})y_{i}\geq 1=\rho\operatorname{tr}(M^{\prime})z$, which is exactly $X_{i}^{\top}MX_{i}+y_{i}\geq z$. The objective value of $(M,y,z)$ is $z-\frac{\left\|y\right\|_{1}}{(1-\epsilon)N}=\frac{1-\left\|y^{\prime}\right\|_{1}}{\rho\operatorname{tr}(M^{\prime})}=1/\rho$. ∎

We will use binary search to find a suitable $\rho>0$, solve the packing/covering SDPs approximately, and then translate the solutions back using Lemma 4.2. The translated solutions will satisfy the conditions of Proposition 4.1. To make sure a suitable $\rho$ exists, we use the following lemma which shows the optimal value of SDPs (6), (7) is continuous and monotone in $\rho$.

Now we are ready to prove Proposition 4.1 by putting Lemmas 4.3 and 4.2 together.

We will define a target interval $[\rho_{1},\rho_{2}]$, such that solving SDPs (6), (7) for any parameters ($\rho\in[\rho_{1},\rho_{2}],\nu,\epsilon)$ will allow us to compute a pair of solutions $w^{\prime}$ and $(M^{\prime},y^{\prime})$ such that:

$w^{\prime}$ is a solution to packing SDP (6) with $\left\|w^{\prime}\right\|_{1}\geq 1-\frac{\epsilon}{10}$; and

$(M^{\prime},y^{\prime})$ is a solution to covering SDP (7) with $\operatorname{tr}(M^{\prime})+\left\|y^{\prime}\right\|_{1}\leq 1$.

It remains to show that we can find a suitable $\rho$ and solve the SDPs (6) (7) in time $\widetilde{O}(Nd/\epsilon^{6})$. We can find $\rho\in[\rho_{1},\rho_{2}]$ using binary search: if $\left\|w^{\prime}\right\|_{1}<1-\frac{\epsilon}{10}$ we will decrease $\rho$, and if $\operatorname{tr}(M^{\prime})+\left\|y^{\prime}\right\|_{1}>1$ we will increase $\rho$.

Finally, we bound from above the running time of the algorithm in this proposition. In each step of the binary search, we solve packing/covering SDPs (6), (7) for some $\rho$. We solve these SDPs to precision $(1\pm O(\epsilon))$ as required in this proof, which takes time $\widetilde{O}(Nd/\epsilon^{6})$, by Corollary 4.5 from Section 4.1. We repeat every use of Corollary 4.5 $O(\log\log(d/\epsilon))$ times, so that the failure probability is at most $1/10$ when we take a union bound over all $O(\log(d/\epsilon))$ iterations. Eventually, when we have a suitable $\rho$, we can convert the solution back to solutions for SDPs (2) (3) using Lemma 4.2. Therefore, the total running time is

In this subsection, we show how to solve packing/covering SDPs (6), (7) in time $\widetilde{O}(Nd/\epsilon^{6})$. It is known that positive (i.e., packing/covering) SDPs can be solved in nearly-linear time and poly-logarithmic number of iterations [JY11, ALO16, PTZ16]. Because SDPs (6), (7) are packing/covering SDPs, we can apply the positive SDP solvers in [PTZ16] directly (Corollary 4.5).

Let $A_{1},\ldots,A_{n}$ be $m\times m$ PSD matrices given in factorized form $A_{i}=C_{i}C_{i}^{\top}$. Consider the following pair of packing and covering SDPs:

An application of the above lemma yields the following corollary:

a $(1+O(\epsilon))$-approximate solution $w^{\prime}$ for the packing SDP (6) with parameters $(\nu,\epsilon,\rho)$; and

a $(1-O(\epsilon))$-approximate solution $(M^{\prime},y^{\prime})$ for the covering SDP (7) with parameters $(\nu,\epsilon,\rho)$.

The total number of non-zeros in all $C_{i}$’s is $q=N(d+1)=O(Nd)$, so by Lemma 4.4, we can solve SDPs (6), (7) in time $\widetilde{O}(Nd/\epsilon^{6})$ with probability $9/10$. Note that the dual solution should be maintained implicitly to avoid writing down an $(N+d)\times(N+d)$ matrix: the top-left block of the dual solution $Y$ is $M^{\prime}$, and the diagonals of the bottom-right block is $y^{\prime}$. ∎

Robust Mean Estimation under Second Moment Assumptions

In this section, we use the algorithmic ideas from Section 3 to establish Theorem 1.3. The algorithm in this case is similar to the one for the sub-gaussian case with some important differences, due to the different concentration properties in the two settings.

Note that it suffices to prove Theorem 1.3 under the assumption that $\sigma=1$, i.e., the covariance satisfies $\Sigma\preceq I$. This is without loss of generality: Given a distribution $D$ with $\Sigma\preceq\sigma^{2}I$, we can first divide every sample by $\sigma$, run the algorithm to learn the mean, and multiply the output by $\sigma$.

We will require a set of conditions to hold for the good samples (similar to Conditions (4) and (5) for sub-gaussian distributions in Section 1.3). We would like the set of good samples to have bounded variance in all directions. However, because we did not make any assumption on the higher moments, it may be possible for a few good samples to affect the empirical covariance too much. Fortunately, such samples have small probability and they do not contribute much to the mean, so we can remove them in the preprocessing step (see Remark 5.1).

where $\delta=c_{1}\sqrt{\epsilon}$ and $\delta_{2}=c_{1}$ for some universal constants $c_{1}$. It follows from Lemma A.18 of [DKK+17] that these conditions will be satisfied with high constant probability after $N=\Omega((d\log d)/\epsilon)$ samples.

The algorithm will be almost identical to Algorithm 1. The only difference is in the “if” statement, where we need a different threshold to decide if the current primal SDP solution is good (or equivalently, whether our guess of $\nu$ is close enough to $\mu^{\star}$).

We are now ready to present our algorithm (Algorithm 2) to robustly estimate the mean of bounded covariance distributions. Algorithm 2 is almost identical to Algorithm 1. The differences are highlighted in bold font. In the if statement, we need a different threshold to decide whether the current primal SDP solution is good (or equivalently, whether our guess of $\nu$ is close enough to $\mu^{\star}$). We use Proposition 5.5 to solve the SDPs. In addition, we need to replace Lemmas 3.2 and 3.3 with Lemmas 5.3 and 5.4, but this merely changes our analysis and has no impact on the algorithm.

We wish to throw away samples that are too far from $\mu^{\star}$. However, we do not know $\mu^{\star}$, so we run the following preprocessing step. We start with an $\epsilon$-corrupted set of $2N$ samples drawn from $D$ (the true distribution) and partition them into two sets of $N$ samples, $S_{1}$ and $S_{2}$. We first compute the coordinate-wise median $\widetilde{\mu}$ of $S_{1}$. Notice that $\left\|\widetilde{\mu}-\mu^{\star}\right\|_{2}\leq O(\epsilon\sqrt{d})$ with high probability. We will use $\widetilde{\mu}$ to throw away samples that are too far from $\mu^{\star}$.

Let $\mathcal{B}$ be a ball of radius $O(\sqrt{d/\epsilon})$ around $\widetilde{\mu}$. Let $G^{\star}_{2}$ denote the original set of uncorrupted samples corresponds to $S_{2}$. By the bounded-covariance assumption, we know that with high probability, $(1-O(\epsilon))$-fraction of the samples in $G^{\star}_{2}$ are in $\mathcal{B}$. Thus, if we change all samples in $S_{2}\setminus\mathcal{B}$ to $\widetilde{\mu}$, the resulting set $S^{\prime}$ is an $O(\epsilon)$-corrupted set of samples drawn from $D$, and all samples in $S^{\prime}$ are not too far from $\mu^{\star}$. Therefore we can use $S^{\prime}$ as the input to Algorithm 2. In the rest of this section, we abuse notation and use $\epsilon$ to denote the fraction of corrupted samples in $S^{\prime}$.

Notation. In this section, we use $c_{1},\ldots,c_{6}$ to denote universal constants. They can be chosen in a way that is similar to how we set constants for Section 3 in Appendix A.

In particular, when $\left\|\mu^{\star}-\nu\right\|_{2}\geq c_{2}\beta$, we have

We take the same feasible primal/dual solutions as in the proof of Lemma 3.1. We get different upper/lower bounds because we use Conditions (8) in this section.

Consider a feasible primal solution $w$ with $w_{i}=\frac{1}{|G|}$ for all $i\in G$ and $w_{i}=0$ otherwise.

One feasible dual solution is $M=yy^{\top}$ where $y=\frac{\mu^{\star}-\nu}{\left\|\mu^{\star}-\nu\right\|_{2}}$. Let $S$ denote the $(1-2\epsilon)N$ good samples with smallest $(X_{i}-\nu)^{\top}M(X_{i}-\nu)$. Let $w^{\prime}_{i}=\frac{1}{(1-2\epsilon)N}$ for all $i\in S$ and $w^{\prime}_{i}=0$ otherwise.

2 When Primal SDP Has Good Solutions

We prove that if the weighted empirical mean is far away from the true mean, then the value of SDP (2) must be large.

The next lemma is similar to Lemma 3.2. The same intuition still holds: if $\epsilon$-fraction of the samples distort the mean by $\Omega(\delta)$, then they must introduce $\Omega(\delta^{2}/\epsilon)$ error to the second moment matrix. Note that because the true covariance matrix is no longer $I$, we cannot say anything about the contribution of the good samples.

Let $w\in\Delta_{N,2\epsilon}$ denote the optimal primal solution. For $y=(\widehat{\mu}_{w}-\mu^{\star})/\left\|\widehat{\mu}_{w}-\mu^{\star}\right\|_{2}$, we have

By the Cauchy-Schwarz inequality and the fact that $w_{B}\leq 2\epsilon$, we get that $\sum_{i\in B}w_{i}\langle X_{i}-\nu,y\rangle^{2}\geq\frac{c_{6}^{2}}{2}(\delta^{2}/\epsilon)$. We conclude the proof by observing that

3 When Primal SDP Has No Good Solutions

We show that when the primal SDP has no good solutions, we can solve the dual (approximately) and the dual will allow us to move $\nu$ closer to $\mu^{\star}$ by a constant factor. The next lemma is similar to Lemma 3.3. The first part of the argument changes slightly because we are using Condition (8), while the second part (the geometric argument) is identical to that of Lemma 3.3.

We know that $M\succeq 0$, $\operatorname{tr}(M)=1$. Without loss of generality, we can assume $M$ is symmetric. By Condition (8), we can prove $\langle M,(\mu^{\star}-\nu)(\mu^{\star}-\nu)^{\top}\rangle\geq\frac{3}{4}\left\|\mu^{\star}-\nu\right\|_{2}^{2}$ as follows.

Therefore, we have a matrix whose inner product with $(\mu^{\star}-\nu)(\mu^{\star}-\nu)^{\top}$ is approximately maximized, this implies that the top eigenvector of $M$ aligns with $(\nu-\mu^{\star})$. We omit the rest of the proof because the geometric analysis is identical to that of Lemma 3.3. ∎

4 Proof of Theorem 1.3

In this section we prove Theorem 1.3 (Correctness and Runtime of Algorithm 2). By combining Lemmas 5.3 and 5.4, we can make progress by either finding a solution to the primal SDP (2) with objective value at most $c_{4}\beta^{2}$, or finding an approximately optimal solution to the dual SDP (3) whose objective value is at least $0.9c_{4}\beta^{2}$. This next proposition shows that this can be done in time $\widetilde{O}(Nd)/\operatorname{poly}(\epsilon)$.

We omit the proof of Proposition 5.5 because its proof is almost identical to the proof of Proposition 4.1. The only difference is that the ratio between the objective values of the desired primal/dual solutions is now $\frac{0.9c_{4}\beta^{2}}{c_{4}\beta^{2}}=0.9$, instead of $\frac{1+0.9c_{4}\beta^{2}}{1+c_{4}\beta^{2}}$ as in Proposition 4.1. The problem of computing a desired pair of solutions becomes easier since the gap is larger.

Theorem 1.3 follows directly from Lemmas 5.3, 5.4, and Proposition 5.5. The running time analysis is identical to that of Theorem 1.2, we can move our guess $\nu$ at most $O(\log(d/\epsilon))$ times, and for each guess we invoke Proposition 5.5 to obtain a good primal or dual solution. The overall running time is $\widetilde{O}(Nd\log(1/\tau)/\epsilon^{6})$.

We note that in Algorithm 2, we are interested in whether the optimal value of SDPs (2) and (3) is at least $c_{4}$ or at most $0.9c_{4}$. Moreover, when we improve our guess $\nu$ using Lemma 5.4, the dual solution $M$ needs to be $0.95$-approximately optimal. Therefore, we only need to solve SDPs (2) and (3) to precision $(1-\epsilon^{\prime})$ for some constant $\epsilon^{\prime}$. However, we do not know how to solve these SDPs directly in $O(Nd)$ time, and when we reduce them to packing/covering SDPs, the constraint on $w$ becomes the objective function in SDP (6), and we must solve SDP (6) more precisely to precision $1-\frac{\epsilon}{10}$ (see, e.g., Lemma 4.2). This is the only reason that we need to pay $\operatorname{poly}(1/\epsilon)$ in the running time of Algorithm 2.

Conclusions and Future Directions

In this paper, we studied the problem of robust high-dimensional mean estimation for structured distribution families in the presence of a constant fraction of corruptions. As our main technical contribution, we gave the first algorithms with dimension-independent error guarantees for this problem that run in nearly-linear time. We hope that this work will serve as the starting point for the design of faster algorithms for high-dimensional robust estimation.

A number of natural directions suggest themselves: Do our techniques generalize to robust covariance estimation? We believe so, but we have not explored this direction in the current work. Can we obtain nearly-linear time robust algorithms for other inference tasks under sparsity assumptions [BDLS17] (e.g., for robust sparse mean estimation or robust sparse PCA)? Can we speed-up the convex programs obtained via the SoS hierarchy in this setting [HL18, KSS18]?

The running time of our algorithms is $\widetilde{O}(Nd)/\operatorname{poly}(\epsilon)$, i.e., it is nearly-linear when the fraction of corruptions $\epsilon$ is constant. Can we avoid the extraneous $\operatorname{poly}(1/\epsilon)$ dependence in the runtime? We believe progress in this direction is attainable. Note that solving a single covering SDP to multiplicative accuracy $(1+\epsilon)$ incurs a $\operatorname{poly}(1/\epsilon)$ slowdown. Is it possible to reframe the underlying optimization problem so that a constant factor multiplicative accuracy suffices? Alternatively, is it possible to speed-up the iterative filtering technique of [DKK+16]? Exploring alternate certificates of robustness may be a promising avenue towards these goals.

Acknowledgments

We thank Alistair Stewart for useful discussions.

References

Appendix A Setting Constants in Section 3

In this section, we describe how to set universal constants $c_{1},\ldots,c_{7}$ in Section 3. The constants are set in the following order: $c_{1}$, $c_{2}$, $c_{4}$, $c_{5}$, $c_{7}$, $c_{6}$, and $c_{3}$. In this order, every $c_{i}$ only depends on the constants set before it, and there are only lower bounds on the value of $c_{i}$, so we can set $c_{i}$ to a sufficiently large constant. Note that $c_{3}$ is the last constant we choose, and our guarantee at the end of the day is to output some hypothesis vector $\widehat{\mu}$ that is close to the true mean $\mu^{\star}$: $\left\|\widehat{\mu}-\mu^{\star}\right\|_{2}\leq c_{3}\delta$.

Recall that in Section 3, $0<\epsilon<1/3$, $\delta=c_{1}\epsilon\sqrt{\ln(1/\epsilon)}$, $\delta_{2}=c_{1}\epsilon\ln(1/\epsilon)$, and $\beta=\sqrt{\epsilon\ln(1/\epsilon)}$.

The constant $c_{1}$ appears in the concentration bounds for the good samples (Condition (4)), and it is related to the constants in Chernoff bounds and Hanson-Wright inequality. We can set $c_{1}$ to be any constant that Condition (4) holds with the right sample complexity.

If we use the dual solution, we know that $r\geq c_{2}\beta$. If we use the primal solution, we have $r\leq c_{5}\beta$. We choose $c_{5}$ where $c_{5}\geq c_{2}$ and $0.9c_{5}^{2}\geq c_{4}$ as needed in the proof of Lemma 3.2.

In the proof of Lemma 3.2, the constants $c_{6}$ and $c_{7}$ appear when we argue that the bad samples contribute at least $\Omega(\delta^{2}/\epsilon)$ to the second-moment, and the good samples contribute at least $1-O(\delta^{2}/\epsilon)$. We choose $c_{7}$ such that $c_{7}\geq 1+\frac{2c_{5}\beta}{\sqrt{\ln(1/\epsilon)}}$, and $c_{6}$ such that $\frac{c_{1}^{2}c_{6}^{2}}{2}\geq c_{4}+c_{1}c_{7}$. Finally, because the good samples shift the mean by at most $\delta$, if the empirical mean is off by more than $c_{3}\delta$ then most of the error are from the bad samples. We choose $c_{3}$ so that $c_{3}\geq c_{6}+1+\frac{2c_{5}\sqrt{\epsilon}}{c_{1}}$.