Introduction

In several modern machine learning (ML) applications, such as ML security [BNJT10, BNL12, SKL17, DKK + 19a] and exploratory analysis of real datasets, e.g., in population genetics [RPW + 02, PLJD10, LAT + 08, DKK + 17], typical datasets contain a non-trivial fraction of arbitrary (or even adversarial) outliers. Robust statistics [HRRS86, HR09] is the subfield of statistics aiming to design estimators that are tolerant to a constant fraction of outliers, independent of the dimensionality of the data. Early work in this field, see, e.g., [Tuk60,Hub64,Tuk75] developed sample-efficient robust estimators for various basic tasks, alas with runtime exponential in the dimension.

During the past five years, a line of work in computer science, starting with [DKK + 16, LRV16], has developed the first computationally efficient robust high-dimensional estimators for a range of tasks. This progress has led to a revival of robust statistics from an algorithmic perspective (see, e.g., [DK19, DKK + 21] for surveys on the topic). In this work, we focus on high-dimensional estimation tasks in the presence of sparsity constraints. To rigorously study these problems, we need to formally define the model of data corruption. Throughout this work, we work with the following standard contamination model. Definition 1.1 (Strong Contamination Model, see [DKK + 16]). Given a parameter 0 < < 1/2 and a distribution family D on R d , the adversary operates as follows: The algorithm specifies a number of samples n, and n samples are drawn from some unknown D ∈ D. The adversary is allowed to inspect the samples, remove up to n of them and replace them with arbitrary points. This modified set of n points is then given as input to the algorithm. We say that a set of samples is -corrupted if it is generated by the above process.

High-dimensional robust statistics is algorithmically challenging because the natural optimization formulations of such tasks are typically non-convex. The recent line of work on algorithmic robust statistics has led to a range of sophisticated algorithms. In some cases, such algorithms require solving large convex relaxations, rendering them computationally prohibitive for large-scale problems. In other cases, they involve a number of hyper-parameters that may require careful tuning. Motivated by these shortcomings of known algorithms, recent work [CDGS20, ZJS20] established an intriguing connection between high-dimensional robust estimation and non-convex optimization. The high-level idea is quite simple: Even though typical robust statistics tasks lead to non-convex formulations, it may still be possible to leverage the underlying structure to show that standard first-order methods provably and efficiently reach near-optimal solutions. Indeed, [CDGS20, ZJS20] were able to prove such statements for robust mean estimation under natural distributional assumptions. Specifically, these works established that any (approximate) stationary point of a well-studied non-convex formulation for robust mean estimation yields a near-optimal solution for the underlying robust estimation task.

In this work, we continue this line of work with a focus on sparse estimation tasks. Leveraging sparsity in high-dimensional datasets is a fundamental problem of significant practical importance. Various formalizations of this problem have been investigated in statistics and machine learning for at least the past two decades (see, e.g., [HTW15] for a textbook on the topic). We focus on robust sparse mean estimation and robust sparse PCA. Sparse mean estimation is arguably one of the most fundamental sparse estimation tasks and is closely related to the Gaussian sequence model [Tsy08,Joh17]. The task of sparse PCA in the spiked covariance model, initiated in [Joh01], has been extensively investigated (see Chapter 8 of [HTW15] and references therein).

In the context of robust sparse mean estimation, we are given an -corrupted set of samples from a distribution with unknown mean µ ∈ R d where µ is k-sparse, and we want to compute a vector µ close to µ. In the context of robust sparse PCA (in the spiked covariance model), we are given an -corrupted set of samples from a distribution with covariance matrix I + ρvv T , where v ∈ R d is k-sparse and the goal is to approximate v. It is worth noting that for both problems, we have access to much fewer samples compared to the non-sparse case (roughly O(k 2 log d) instead of Ω(d)). Consequently, the design and analysis of optimization formulations for robust sparse estimation requires new ideas and techniques that significantly deviate from the standard (non-sparse) case.

Our Results and Contributions

We show that standard first-order methods lead to robust and efficient algorithms for sparse mean estimation and sparse PCA. Our main contribution is to propose novel (non-convex) formulations for these robust estimation tasks, and to show that approximate stationarity suffices for near-optimality. We establish landscape results showing that any approximate stationary point of our objective function yields a near-optimal solution for the underlying robust estimation task. Consequently, gradient descent (or any other methods converging to stationarity) can solve these problems.

Our results provide new insights and techniques in designing and analyzing (non-convex) optimization formulations of robust estimation tasks. Our formulations and structural results immediately lead to simple and practical algorithms for robust sparse estimation. Importantly, the gradient of our objectives can be computed efficiently via a small number of basic matrix operations. In addition to their simplicity and practicality, our methods provably succeed under more general distributional assumptions compared to prior work.

For robust sparse mean estimation and robust sparse PCA, our landscape results require deterministic conditions on the original set of good samples. We refer to these conditions as stability conditions (Definitions 2.1 and 2.2, formally defined in Section 2). At a high level, they state that the first and second moments of a set of samples are stable when any -fraction of the samples are removed. These stability conditions hold with high probability for a set of clean samples drawn from natural families of distributions (e.g., subgaussian).

For robust sparse mean estimation, we establish the following result. Theorem 1.2 (Robust Sparse Mean Estimation). Let 0 < < 0 for some universal constant 0 and let δ > . Let G be a set of n samples that is (k, , δ)-stable (per Definition 2.1) w.r.

t. a distribution with unknown

There is an algorithm that on inputs S, k, , and δ, runs in polynomial time and returns a k-sparse vector

We emphasize that a key novelty of Theorem 1.2 is that the underlying algorithm is a first-order method applied to a novel non-convex formulation of the problem. The major advantage of our algorithm over prior work [BDLS17, DKK + 19b] is its simplicity, practicality, and the fact that it seamlessly applies to a wider class of distributions on the clean data.

As we will discuss in Section 3, when the ground-truth distribution D is subgaussian with unknown k-sparse mean µ ∈ R d and identity covariance, a set of n = Ω(kfoot_2 log d/ 2 ) samples drawn from D is (k, , δ)-stable (Definition 2.1) with high probability for δ = O( log(1/ )). It follows as an immediate corollary of Theorem 1.2 that, given an -corrupted set of samples, we can compute a vector µ that is O(δ) = O( log(1/ )) close to the true mean µ. This sample complexity matches the known computational-statistical lower bounds [DKS17, BB20]. More generally, one can relax the concentration assumption on the clean data and obtain qualitatively similar error guarantees.

Next we state our main result for robust sparse PCA. Theorem 1.3 (Robust Sparse PCA). Let 0 < ρ ≤ 1 and 0 < < 0 for some universal constant 0 . Let G be a set of n samples that is (k, , δ)-stable (as in Definition 2.2) w.r.t. a centered distribution with covariance Σ = I + ρvv , for an unknown k-sparse

There is an algorithm that on inputs S, k, and , runs in polynomial time and returns a unit vector

Interestingly, our algorithm for robust sparse PCA is a first-order method applied to a simple convex formulation of the problem. We view the existence of a convex formulation as an intriguing fact that, surprisingly, was not observed in prior work. As we will discuss in Section 4, when the ground-truth distribution D is centered subgaussian with covariance Σ = I + ρvv , for an unknown k-sparse unit vector v ∈ R d , a set of n = Ω(k 2 log d/ 2 ) samples drawn from D is (k, , δ)-stable (Definition 2.2) with high probability for δ = O( log(1/ )). Therefore, our algorithm outputs a vector that is O( log(1/ )/ρ) close to the true direction v. The sample complexity in this case nearly matches the computational-statistical lower bound of Ω(k 2 log d/ 2 ) [BR13] which holds even without corruptions. While the error guarantee of our algorithm is slightly worse compared to prior work [BDLS17, DKK + 19b] for Gaussian data (we get O( δ/ρ) rather than O(δ/ρ)), we note that our algorithm works for a broader family of distributions.

Prior Work on Robust Sparse Estimation. We provide a detailed summary of prior work for comparison. [BDLS17] obtained the first sample-efficient and polynomial-time algorithms for robust sparse mean estimation and robust sparse PCA. These algorithms succeed for Gaussian inliers and inherently use the ellipsoid method. The separation oracle required for the ellipsoid algorithm turns out to be another convex program -corresponding to an SDP to solve sparse PCA. As a consequence, the running time of these algorithms, while polynomially bounded, is impractically high. [LLC19] proposed an algorithm for robust sparse mean estimation via iterative trimmed hard thresholding, which can only tolerate a sub-constant fraction of corruptions. [DKK + 19b] gave iterative spectral robust algorithms for sparse mean estimation and sparse PCA. These algorithms are still quite complex and are only shown to succeed under Gaussian inliers.

Overview of Our Approach

In this section, we give an overview of our approach for robust sparse mean estimation. At a very high level, we assign a nonnegative weight to each data point and try to find a good set of (1 -)n samples. The constraint on the weight vector is that it represents at least a (fractional) set of (1 -)-portion of the input dataset. Formally, given n datapoints (X i ) n i=1 , the goal is to find a weight vector w ∈ R n such that µ w = i w i X i is close to the true mean µ. The constraint on w is that it belongs to

which is the convex hull of all uniform distributions over subsets S ⊆ [n] of size |S| = (1 -)n.

Let Σ w = i w i (X i -µ w )(X i -µ w ) denote the weighted empirical covariance matrix. It is well-known that if one can find w ∈ ∆ n, that minimizes the weighted empirical variance v Σ w v for all k-sparse unit vectors v, then µ w must be close to µ. Unfortunately, it is NP-Hard to find the sparse direction v with the largest variance. To get around this issue, [BDLS17] considered the following convex relaxation, minimizing the variance for convex combinations of sparse directions:

Given w, the optimal A can be found using semidefinite programming (SDP). [ZJS20] observed that any stationary point w of (1) gives a good solution for robust sparse mean estimation. However, solving (1) requires convex programming to compute the gradient in each iteration. As explained in the proceeding discussion, our approach circumvents this shortcoming, leading to a formulation for which each gradient can be computed using only basic matrix operations.

In this work, we propose and analyze the following optimization formulation:

where A F,k,k is the Frobenius norm of the k 2 entries of A with largest magnitude, with the additional constraint that these k 2 entries are chosen from k rows with k entries in each row.

We prove that any stationary point of f (w) yields a good solution for robust sparse mean estimation.

Here we provide a brief overview of our proof (see Section 3 for more details). Given a weight vector w, we show that if w is not a good solution, then moving toward w (the weight vector corresponding to the uniform distribution on the clean input samples) will decrease the objective value. Formally, we will show that, for any 0 < η < 1,

We can then take • F,k,k norm on both sides (after subtracting I) and show that the third term can be essentially ignored. If the third term were not there, we would have

Therefore, if w is a bad solution with f (w) much larger than f (w ), then w cannot be a stationary point because f decreases when we move from w to (1 -η)w + ηw .

Remark 1.4. The technical overview for robust sparse PCA follows a similar high-level approach, but is somewhat more technical. It is deferred to Section 4.

Roadmap.

In Section 2, we introduce basic notations and the deterministic stability conditions that we require on the good samples. We present our algorithms and analysis for robust sparse mean estimation in Section 3 and robust sparse PCA in Section 4. In Section 5, we evaluate our algorithm on synthetic datasets and show that it achieves good statistical accuracy under various noise models.

Preliminaries and Background

Notation. For a positive integer n, let [n] = {1, . . . , n}. For a vector v, we use v 0 , v 1 , v 2 , and v ∞ for the number of non-zeros, the 1 , 2 , and ∞ norm of v respectively. Let I be the identity matrix. For a matrix A, we use A 2 , A F , tr(A) for the spectral norm, Frobenius norm, and trace of A respectively. For two vectors x, y, let x y denote their inner product. For two matrices A, B, we use A • B = tr(A B) for their entrywise inner product. A matrix A is said to be positive semidefinite (PSD) if x Ax ≥ 0 for all x. We write A B iff (B -A) is PSD.

For a vector w ∈ R n , let diag(w) ∈ R n×n denote a diagonal matrix with w on the diagonal. For a matrix A ∈ R n×n , let diag(A) ∈ R n denote a column vector with the diagonal of A. For a vector v ∈ R d and a set S ⊆ [d], we write v S ∈ R d for a vector that is equal to v on S and zero everywhere else. Similarly, for a matrix A ∈ R d×d and a set S ⊆ ([d] × [d]), we write A S for a matrix that is equal to A on S and zero everywhere else.

For a vector v, we define v 2,k = max |S|=k v S 2 to be the maximum 2 -norm of any k entries of v. For a matrix A, we define A F,k 2 to be the maximum Frobenius norm of any k 2 entries of A. Moreover, we define A F,k,k to be the maximum Frobenius norm of any k 2 entries with the extra requirement that these entries must be chosen from k rows with k entries in each row. Formally,

Sample Reweighting Framework. We use n for the number of samples, d for the dimension, and for the fraction of corrupted samples. For sparse estimation, we use k for the sparsity of the ground-truth parameters. We use G for the original set of n good samples. We use S = G ∪ B for the input samples after the adversary replaced -fraction of G , where G ⊂ G is the set of remaining good samples and B is the set of bad samples (outliers) added by the adversary. Note that |G| = (1 -)n and |B| = n.

Given n samples X 1 , . . . , X n , we write X ∈ R d×n as the sample matrix where the i-th column is X i . For a weight vector w ∈ R n , we use µ w = Xw = i w i X i for the weighted empirical mean and

for the weighted empirical covariance. Let ∆ n, be the convex hull of all uniform distributions over subsets

In other words, every w ∈ ∆ n, corresponds to a fractional set of (1 -)n samples. We use w to denote the uniform distribution on G (the remaining good samples in S).

Deterministic Stability Conditions.

For robust sparse mean estimation and robust sparse PCA, we require the following conditions respectively. Definition 2.1 (Stability Conditions for Sparse Mean). A set of n samples G = (X i ) n i=1 is said to be (k, , δ)-stable (w.r.t. a distribution with mean µ) iff for any weight vector w ∈ ∆ n,2 , we have µ w -µ 2,k ≤ δ and Σ w -I F,k,k ≤ δ 2 / , where µ w and Σ w are the weighted empirical mean and covariance matrix respectively, and the • F,k,k norm is defined in Equation (2).

Definition 2.2 (Stability Conditions for Sparse PCA). A set of

First-Order Stationary Points. We give a formal definition of the notion of (approximate) firstorder stationary point that we use in this paper. Definition 2.3 (Approximate Stationary Points). Fix a convex set K and a differentiable function f . For γ ≥ 0, we say that x ∈ K is a γ-stationary point of f iff the following condition holds: For any unit vector u where x + αu ∈ K for some α > 0, we have u ∇f (x) ≥ -γ.

We note that the objective functions studied in this paper are not everywhere differentiable. This is because, taking the • F,k,k norm as an example, there can be ties in choosing the largest k 2 entries. When the function f is not differentiable, we use ∇f informally to denote an element of the sub-differential. We will show in Appendix C that, while f is not differentiable, it does have a nonempty subdifferential, as it can be written as the pointwise maximum of differentiable functions.

Robust Sparse Mean Estimation

In this section, we present our non-convex approach for robust sparse mean estimation. We will optimize the following objective, where • F,k,k is defined in Equation (2):

We will show that the objective function (3) has no bad stationary points (Theorem 3.1). In other words, every first-order stationary point of f yields a good solution for robust sparse mean estimation.

Our algorithm is stated in Algorithm 1. As a consequence of our landscape result (Theorem 3.1), we know that Algorithm 1 works no matter how we find a stationary point of f (because any stationary point works), so we intentionally did not specify how to find such a point. As a simple illustration, we show that (projected) gradient descent can be used to minimize f . The convergence analysis and iteration complexity are provided in Appendix C.

Algorithm 1: Robust sparse mean estimation.

Input: k > 0, 0 < < 0 , and an -corrupted set of samples (X i ) n i=1 drawn from a distribution with k-sparse mean µ.foot_3 Output: a vector µ that is close to µ.

1: Find a first-order stationary point w ∈ ∆ n, of the objective min w f (w) = Σ w -I F,k,k . 2: Return µ = (µ w ) Q where Q is a set of k entries of µ w with largest magnitude.

Formally, we first prove that Algorithm 1 can output a vector µ ∈ R d that is close to µ in • 2,k norm, as long as the good samples satisfies the stability condition in Definition 2.1. Theorem 3.1. Fix k > 0, 0 < < 0 , and δ > . Let G be a set of n samples that is (k, , δ)-stable (as in Definition 2.1) w.r.t. a distribution with unknown k-sparse

Once we have a vector µ w that is O(δ)-close to µ in • 2,k norm, we can guarantee that a truncated version of µ w (the output µ of Algorithm 1) is O(δ)-close to µ in the 2 -norm: Lemma 3.2. Fix two vectors x, y with x 0 ≤ k and x -y 2,k ≤ δ. Let z be a vector that keeps the k entries of y with largest absolute values and sets the rest to 0. We have x -z 2 ≤ √ 5δ.

Theorem 1.2 follows immediately from Theorem 3.1 and Lemma 3.2.

We can apply Theorem 1.2 to get an end-to-end result for subgaussian distributions. We show that the required stability conditions are satisfied with a small number of samples. Lemma 3.3. Fix k > 0 and 0 < < 0 . Let G be a set of n samples that are drawn i.i.d. from a subgaussian distribution with mean µ and covariance I. If n = Ω(k 2 log d/ 2 ), then with probability at least 1 -exp(-Ω(k 2 log d)), G is (k, , δ)-stable (as in Definition 2.1) for δ = O( log(1/ )).

Combining Theorem 1.2 and Lemma 3.3, we know that given an -corrupted set of O(k 2 log d/ 2 ) samples drawn from a subgaussian distribution with k-sparse mean µ, the output of Algorithm 1 is O( log(1/ ))-close to µ in 2 -norm.

In the rest of this section, we will prove Theorem 3.1. Omitted proofs in this section are in Appendix A.

We start with some intuition on why we choose our objective function (3). We would like to design f (w) = g(Σ w -I) to satisfy the following properties:

1. g(Σ w -I) is an upper bound on v (Σ w -I)v for all k-sparse unit vectors v ∈ R d . This way, a small objective value implies that µ w -µ 2,k is small.

2. g(Σ w -I) is small for w (the uniform distribution on G). This guarantees that a good w exists.

3. Triangle inequality on g. This allows us to upper bound the objective value when we move w toward w by the sum of g(•) of each term on the right-hand side:

4. g(uu ) is close to g(vv ) where v keeps only the k largest entries of u. We want to approximate µ in • 2,k norm, so intuitively g(Σ w -I) should depend only on the largest k entries of (µ w -µ).

The last step uses ( 1 2 -4 ) Σ w -I 2 ≥ (4c 1 + 1) δ 2 which holds if ≤ 1/10 and c 2 2 ≥ 164 c 1 + 40.

It follows immediately that w cannot be a stationary point. Let u = w -w w -w 2 and h = η w -w 2 .

We have w

By Definition 2.3, we know w cannot be a γ-stationary point of f for some γ = O(n 1/2 δ 2 -3/2 ).

Robust Sparse PCA

We consider a spiked covariance model for sparse PCA. In this model, there is a direction v ∈ R d with at most k nonzero entries. The good samples are drawn from a ground-truth distribution with covariance Σ = I + ρvv , where ρ > 0 is a parameter that intuitively measures the strength of the signal. We consider the more interesting case when ρ ≤ 1 (if ρ is larger the problem becomes easier).

To solve the sparse PCA problem, we consider the following optimization problem, where

The objective function minimizes the Frobenius norm of the largest 2k 2 entries of a reweighted second-moment matrix M w . Note that f (w) is actually convex in w, because the matrix M w is linear in w and the • F,2k 2 norm is convex.

Let R be the support of vv . Intuitively, the k 2 entries in R could be large due to spiked covariance. By minimizing the norm of the largest 2k 2 entries, we hope to make the entries outside of R very small. Our algorithm is given in Algorithm 2.

Algorithm 2: Robust sparse PCA.

Input: k > 0, 0 < < 0 , and an -corrupted set of samples (X i ) n i=1 drawn from a distribution with covariance I + ρvv for a k-sparse unit vector v. Output: a vector u that is close to v.

1: Find a first-order stationary point w ∈ ∆ n, of the objective min w f (w) = M w -I F,2k 2 . 2: Let A = M w -I. Let Q be the k 2 entries of A with largest magnitude.

3: Return u = the top eigenvector of (A Q + A Q ).

Theorem 4.1. Let 0 < ρ ≤ 1, 0 < < 0 , and δ > . Let G be a set of n samples that is (k, , δ)stable (as in Definition 2.2) w.r.t. a centered distribution with covariance I + ρvv for an unknown k-sparse unit vector v ∈ R d . Let S = (X i ) n i=1 be an -corrupted version of G . Algorithm 2 outputs a vector u such that uu -vv F = O( δ/ρ).

Theorem 1.3 is an immediate corollary of Theorem 4.1.

We can apply Theorem 4.1 to get an end-to-end result for subgaussian distributions. Algorithm 2 requires the stability conditions (Definition 2.2) of the original good samples G . We show that these conditions are satisfied with a small number of samples. Lemma 4.2. Let 0 < ρ ≤ 1 and 0 < < 0 . Let D be a centered subgaussian distribution with covariance I + ρvv for a k-sparse unit vector v ∈ R d . Let G be a set of n = Ω(k 2 log d/δ 2 ) samples drawn from D. Then then with probability at least 1 -exp(-Ω(k 2 log d)), G is (k, , δ)stable (as in Definition 2.2) w.r.t. D for δ = O( log(1/ )).

Combining Theorem 4.1 and Lemma 4.2, given as input an -corrupted set of n = Ω(k 2 log d/ 2 ) samples drawn from a centered subgaussian distribution with covariance I + ρvv , Algorithm 2 returns a vector u with uu -vv F = O( log(1/ )/ρ).

We defer the proofs of Lemma 4.2 and Theorem 4.1 to Appendix B and give an overview of the proof of Theorem 4.1.