Robust Sparse Estimation Tasks in High Dimensions

Introduction

In the last couple of decades, there has been a large amount of work in machine learning and statistics on how to exploit sparsity in high dimensional data analysis. Motivated by the ever-increasing quantity and dimensionality of data, the goal at a high level is to utilize the underlying sparsity of natural data to extract meaningful guarantees using a number of samples that is sublinear in the dimensionality of the data. In this paper, we will consider the unsupervised setting, where we have sample access to some distribution with some underlying sparsity, and our goal is to recover this distribution by exploiting this structure. Two natural and well-studied problems in this setting that attempt to exploit sparsity are sparse mean estimation and sparse PCA. In both problems, the shared theme is that we assume that one wishes to find a distinguished sparse direction of a Gaussian data set. However, the algorithms inspired by this line of work tend to be quite brittle—it can be shown that they fail when the model is slightly perturbed.

This connects to a major concern in high dimensional data analysis: that of model misspecification. At a high level, the worry is that our algorithms should be able to tolerate the case when our assumed model and the true model do not perfectly coincide. In the distributional setting, this (more or less) corresponds to the regime when a small fraction of our samples are adversarially corrupted. The study of these so-called robust estimators, i.e., estimators which work in the presence of such noise, is a classical subfield of statistics. Unfortunately, the classical algorithms for these problems fail to scale as the dimensionality of the problem grows—either the algorithms run in time which is exponential in the dimension, or the error guarantees for these algorithms degrade substantially as the dimension grows. In a flurry of recent work, we now know new algorithms which circumvent this “curse of dimensionality”: they run efficiently, and provide dimension independent error guarantees. However, these algorithms are unable to exploit any inherent sparsity in the problem.

Do the statistical gains (achievable by computationally efficient algorithms) for sparse estimation problems persist in the presence of noise?

More formally: Suppose we are asked to solve some estimation task given samples from some distribution $D$ with some underlying sparsity constraint (e.g. sparse PCA). Suppose now an $\varepsilon$-fraction of the samples are corrupted. Can we still solve the same sparse estimation problem? In this work, we initiate the study of such issues. Interestingly, new gaps between computational and statistical rates seem to emerge in the presence of noise. In particular, while the sparse mean estimation problem was previously quite simple to solve, the efficient algorithms which achieve the minimax rate for this problem break down in the presence of this adversarial noise. More concretely, it seems that the efficient algorithms which are robust to noise run into the same computational issues as those which plague sparse PCA. A very interesting question is whether this phenomenon is inherent to any computationally efficient algorithm.

We study the natural robust versions of two classical, well-studied statistical tasks involving sparsity, namely, sparse mean estimation, and sparse PCA.

Here, we get a set of $d$-dimensional samples from $\operatorname{\mathcal{N}}(\mu,I)$, where $\mu$ is $k$-sparse, and an $\varepsilon$-fraction of the points are corrupted adversarially. Our goal then is to recover $\mu$. Our main contribution is the following:

There is an efficient algorithm, which given a set of $\varepsilon$-corrupted samples of size $\widetilde{O}(\frac{k^{2}\log d}{\varepsilon^{2}})$ from $\operatorname{\mathcal{N}}(\mu,I)$ where $\mu$ is $k$-sparse, outputs a $\widehat{\mu}$ so that with high probability, $\|\widehat{\mu}-\mu\|_{2}\leq\varepsilon\sqrt{\log 1/\varepsilon}$.

Any efficient algorithm for robust sparse mean estimation needs $\widetilde{\Omega}(\frac{k^{2}\log d}{\varepsilon^{2}})$ samples.

In Appendix D we give some intuition for why it seems to be true. At a high level, it seems that any technique to detect outliers for the mean must look for sparse directions in which the variance is much larger than it should be; at which point the problem faces the same computational difficulties as sparse PCA. We leave closing this gap as an interesting open problem.

Robust sparse PCA

Here, we study the natural robust analogue of the spiked covariance model. Classically, two problems are studied in this setting. The detection problem is given as follows: given sample access to the distributions, we are asked to distinguish between $\operatorname{\mathcal{N}}(0,I)$, and $\operatorname{\mathcal{N}}(0,I+\rho vv^{T})$ where $v$ is a $k$-sparse unit vector. That is, we wish to understand if we can detect the presence of any sparse principal component. Our main result is the following:

Fix $\rho>0$, and let $\eta=O(\varepsilon\sqrt{\log 1/\varepsilon})$. If $\rho>\eta$, there is an efficient algorithm, which given a set of $\varepsilon$-corrupted samples of size $O(\frac{k^{2}\log d}{\rho^{2}})$ which distinguishes between $\operatorname{\mathcal{N}}(0,I)$, and $\operatorname{\mathcal{N}}(0,I+\rho vv^{T})$ with high probability.

The condition that $\varepsilon=\widetilde{O}(\rho)$ is necessary (up to log factors), as otherwise the problem is impossible information theoretically. Observe that this (up to log factors) matches the optimal rate for computationally efficient detection for sparse PCA without noise (under reasonable complexity theoretic assumptions, see [BR13, WBS16]), and so it seems that noise does not introduce an additional gap here. The recovery problem is similar, except now we want to recover the planted spike $v$, i.e. find a $u$ minimizing $L(u,v)=\frac{1}{\sqrt{2}}\|uu^{T}-vv^{T}\|$, which turns out to be the natural measure for this problem. For this, we show:

Fix $\varepsilon>0$ and $0<\rho=O(1)$, and let $\eta=O(\varepsilon\sqrt{\log 1/\varepsilon})$. There is an efficient algorithm, which given a set of $\varepsilon$-corrupted samples of size $O(\frac{k^{2}\log d}{\eta^{2}})$ from $\operatorname{\mathcal{N}}(0,I+\rho vv^{T})$, outputs a $u$ so that $L(u,v)=O\left(\frac{\eta}{\rho}\right)$ with high probability.

This rate is non-trivial—in particular, it provides guarantees for recovery of $v$ when the number of samples we take is at the detection threshold. Moreover, up to log factors, our rate is optimal for computationally efficient algorithms–[WBS16] gives an algorithm with rate roughly $O(\varepsilon/\rho)$, and show that this is necessary.

Techniques

We first introduce a simple way to describe the optimization problems used for solving sparse mean estimation and sparse PCA. This approach is very similar to the approach taken by [CRPW12] for solving under-determined linear systems. We observe that any set $\mathcal{S}$ in a Hilbert space naturally induces a dual norm $\|x\|^{*}_{\mathcal{S}}=\max_{y\in\mathcal{S}}|\langle x,y\rangle|,$ and that well-known efficient algorithms for sparse mean estimation and sparse PCA simply compute this norm, and the corresponding dual witness $y\in\mathcal{S}$ which maximizes this norm, for appropriate choices of $\mathcal{S}$. These norms give us a language to only consider deviations in directions we care about, which allows us to prove concentration bounds which are not true for more traditional norms.

We now describe our techniques for robust sparse mean estimation. Our starting point is the convex programming approach of [DKK+16]. We assign each sample point a weight, which morally corresponds to our belief about whether the point is corrupted, and we optimize these weights. In previous work of [DKK+16], the approach was to find weights so that the empirical covariance with these weights looked like the identity in spectral norm.

Unfortunately, such an approach fundamentally fails for us because the spectrum of the covariance will never concentrate for us with the number of samples we take. Instead, we utilize a novel connection to sparse PCA. We show that if instead we find weights so that the empirical covariance with these weights looks like the identity in the dual norm induced by a natural SDP for sparse PCA (in the noiseless setting), then this suffices to show that the trucnated empirical mean with these weights is close to the truth. We do so by convex programming. While we cannot explicitly write down the feasible set of weights, it is a convex set. Thus, by the classical theory of convex optimization, it suffices to give a separation oracle for this convex set to optimize over this set. We show that in fact the SDP for sparse PCA gives us such a separation oracle, if one is sufficiently careful to always work with sparsity preserving objects. This in turns suffices to allow us to (approximately) find a point in the desired feasible set of points, which we show suffices to recover the true mean.

We now turn to robust sparse PCA. We first consider the detection problem, which is somewhat easier technically. Here, we again use the dual norm induced by the SDP for sparse PCA. We show that if we can find weights on the samples (as before) so that the empirical covariance with these samples has minimal dual norm, then the value of the dual norm gives us a distinguisher between the spiked and non-spiked case. To find such a set of weights, we observe that norms are convex, and thus our objective is convex. Thus, as before, to optimize over this set it suffices to give a separation oracle, which again the SDP for sparse PCA allows us to do.

We now turn our attention to the recovery problem. Here, the setup is very similar, except now we simultaneously find a set of weights and an “explainer” matrix $A$ so that the empirical covariance with these weights is “maximally explained” by $A$, in a norm very similar to the one induced by the sparse PCA SDP. Utilizing that norms are convex, we show that this can be done via a convex program using the types of techniques described above, and that the top eigenvector of the optimal $A$ gives us the desired solution. While the convex program would be quite difficult to write down in one shot, it is quite easily expressible using the abstraction of dual norms.

2 Related Work

As mentioned previously, there has been a large amount of work on various ways to exploit sparsity for machine learning and statistics. In the supervised setting, perhaps the most well-known of these is compressive sensing and its variants (see [CW08, HTW15] for more details). We do not attempt to provide an exhaustive overview the field here. Other well-known problems in the same vein include general classes of linear inverse problems, see [CRPW12] and matrix completion ([CR12]).

The question of estimating a sparse mean is very related to a classical statistical model known as the Gaussian sequence model, and the reader is referred to [Tsy09, Joh11, Rig15] for in-depth surveys on the area. This problem has also garnered a lot of attention recently in various distributed and memory-limited settings, see [GMN14, SD15, BGM+16]. The study of sparse PCA was initiated in [Joh01] and since yielded a very rich algorithmic and statistical theory ([dEGJL07, dBG08, AW08, WTH09, JNRS10, ACCD11, LZ12, Ma13, BJNP13, CMW13, OMH14, GWL14, CRZ16, BMVX16, PWBM16]). In particular, we highlight a very interesting line of work [BR13, KNV15, MW15, WGL15, WBS16], which give evidence that any computationally efficient estimator for sparse PCA must suffer a sub-optimal statistical rate rate. We conjecture that a similar phenomenon occurs when we inject noise into the sparse mean estimation problem.

In this paper we consider the classical notion of corruption studied in robust statistics, introduced back in the 70’s in seminal works of [HR09, Tuk75, HRRS86]. Unfortunately, essentially all robust estimators require exponential time in the dimension to compute ([JP78, Ber06, HM13]). Subsequent work of [LT15, BD15] gave efficient SDP-based estimators for these problems which unfortunately had error guarantees which degraded polynomially with the dimension. However, a recent flurry of work ([DKK+16, LRV16, CSV16, DKK+17, DKS17, DKS16]) have given new, computationally efficient, robust estimators for these problems and other settings which avoid this loss, and are often almost optimal. Independent work of [DSS17] also considers the robust sparse setting. They give a similar result for robust mean estimation, and also consider robust sparse PCA, though in a somewhat different setting than we do, as well as robust sparse linear regression.

The questions we consider are similar to learning in the presence of malicious error studied in [Val85, KL93], which has received a lot of attention, particularly in the setting of learning halfspaces ([Ser03, KLS09, ABL14]). They also are connected to work on related models of robust PCA ([Bru09, CLMW11, LMTZ12, ZL14]). We refer the reader to [DKK+16] to a detailed discussion on the relationships between these questions and the ones we study.

Definitions

We will study the following contamination model:

As discussed in [DKK+16], this is a strong notion of sample corruption that is able to simulate previously defined notions of error. In particular, this can simulate (up to constant factors) the scenario when our samples do not come from $D$, but come from a distribution $D^{\prime}$ with total variation distance at most $O(\varepsilon)$ from $D$.

We may now formally define the algorithmic problems we consider.

there is a poly-time algorithm which outputs $\widehat{\mu}$ so that w.p. $1-\delta$, we have $\|\mu-\widehat{\mu}\|_{2}\leq O(\eta)$.

It is well-known that information theoretically, the best error one can achieve is $\Theta(\varepsilon)$, as achieved by Fact A.1. We show that it is possible to efficiently match this bound, up to a $\sqrt{\log 1/\varepsilon}$ factor. Interestingly, our rate differs from that in Fact A.1: our sample complexity is (roughly) $\widetilde{O}(k^{2}\log d/\varepsilon^{2})$ versus $O(k\log d/\varepsilon^{2})$. We conjecture this is necessary for any efficient algorithm.

Robust sparse PCA

We will consider both the detection and recovery problems for sparse PCA. We first focus detection problem for sparse PCA. Here, we are given $\rho>0$, and an $\varepsilon$-corrupted set of samples from a $d$-dimensional distribution $D$, where $D$ can is either $\operatorname{\mathcal{N}}(0,I)$ or $\operatorname{\mathcal{N}}(0,I+\rho vv^{T})$ for some $k$-sparse unit vector $v$. Our goal is to distinguish between the two cases, using as few samples as possible. It is not hard to show that information theoretically, $O(k\log d/\rho^{2})$ samples suffice for this problem, with an inefficient algorithm (see Appendix A). Our first result is that efficient robust sparse PCA detection is possible, at effectively the best computationally efficient rate:

Fix $\rho,\delta,\varepsilon>0$. Let $\eta=O(\varepsilon\sqrt{\log 1/\varepsilon})$. Then, if $\eta=O(\rho)$, and we are given a we are given a $\varepsilon$-corrupted set of samples from either $\operatorname{\mathcal{N}}(0,I)$ or $\operatorname{\mathcal{N}}(0,I+\rho vv^{T})$ for some $k$-sparse unit vector $v$ of size

then there is a polynomial time algorithm which succeeds with probability $1-\delta$ for detection.

It was shown in [BR13] that even without noise, at least $n=\Omega(k^{2}\log d/\varepsilon^{2})$ samples are required for any polynomial time algorithm for detection, under reasonable complexity theoretic assumptions. Up to log factors, we recover this rate, even in the presence of noise.

We next consider the recovery problem. Here, we are given an $\varepsilon$-corrupted set of samples from $\operatorname{\mathcal{N}}(0,I+\rho vv^{T})$, and our goal is to output a $u$ minimizing $L(u,v)$, where $L(u,v)=\frac{1}{\sqrt{2}}\|uu^{T}-vv^{T}\|$. For the recovery problem, we recover the following efficient rate:

Fix $\varepsilon,\rho>0$. Let $\eta$ be as in Theorem 2.2. There is an efficient algorithm, which given a set of $\varepsilon$-corrupted samples of size $n$ from $\operatorname{\mathcal{N}}(0,I+\rho vv^{T})$, where

In particular, observe that when $\eta=O(\rho)$, so when $\varepsilon=\widetilde{O}(\rho)$, this implies that we recover $v$ to some small constant error. Therefore, given the same number of samples as in Theorem 2.2, this algorithm begins to provide non-trivial recovery guarantees. Thus, this algorithm has the right “phase transition” for when it begins to work, as this number of samples is likely necessary for any computationally efficient algorithm. Moreover, our rate itself is likely optimal (up to log factors), when $\rho=O(1)$. In the non-robust setting, [WBS16] showed a rate of (roughly) $O(\varepsilon/\rho)$ with the same number of samples, and that any computationally efficient algorithm cannot beat this rate. We leave it as an interesting open problem to show if this rate is achievable or not in the presence of error when $\rho=\omega(1)$.

Preliminaries

In this section we provide technical preliminaries that we will require throughout the paper.

We will require the following (straightforward) preprocessing subroutine from [DKK+16] to remove all points which are more than $\widetilde{\Omega}(d)$ away from the true mean.

Let $X_{1},\ldots,X_{n}$ be an $\varepsilon$-corrupted set of samples from $\operatorname{\mathcal{N}}(\mu,I)$, and let $\delta>0$. There is an algorithm $\textsc{NaivePrune}(X_{1},\ldots,X_{n},\delta)$ which runs in $O(\varepsilon d^{2}n^{2})$ time so that with probability $1-\delta$, we have that (1) NaivePrune removes no uncorrupted points, and (2) if $X_{i}$ is not removed by NaivePrune, then $\|X_{i}-\mu\|_{2}\leq O(\sqrt{d\log(n/\delta)})$. If these two conditions happen, we say that NaivePrune has succeeded.

2 Concentration inequalities

In this section we give a couple of concentration inequalities that we will require in the remainder of the paper. These “per-vector” and “per-matrix” concentration guarantees are well-known and follow from (scalar) Chernoff bounds, see e.g. [DKK+16].

Then, with probability $1-\delta$, we have

Then, with probability $1-\delta$, we have:

We make the following observation. For any subset $I\subseteq[n]$, if we let $w^{I}$ be the vector whose $i$th coordinate is $1/|I|$ if $i\in I$ and otherwise, we have

The set $S_{n,\varepsilon}$ will play a key role in our algorithms. We will think of elements in $S_{n,\varepsilon}$ as weights we place upon our sample points, where higher weight indicates a higher confidence that the sample is uncorrupted, and a lower weight will indicate a higher confidence that the sample is corrupted.

Concentration for sparse estimation problems via dual norms

In this section we give a clean way of proving concentration bounds for various objects which arise in sparse PCA and sparse mean estimation problems. We do so by observing they are instances of a very general “meta-algorithm” we call dual norm maximization. This will prove crucial to proving the correctness of our algorithms for robust sparse recovery. While this may sound similar to the “dual certificate” techniques often used in the sparse estimation literature, these techniques are actually quite different.

Let $\mathcal{H}$ be a Hilbert space with inner product $\langle\cdot,\cdot\rangle$. Fix any set $S\subseteq\mathcal{H}$. Then the dual norm induced by $S$, denoted $\|\cdot\|_{S}^{*}$, is defined by $\|x\|_{S}^{*}=\sup_{y\in S}|\langle x,y\rangle|$. The dual norm maximizer of $x$, denoted $d_{S}(x)$, is the vector $d_{S}(x)=\operatorname*{arg\,max}_{v\in S}|\langle v,x\rangle|$.

We show in Appendix B.1 that existing well-known algorithms for sparse mean recovery and sparse PCA without noise can be naturally written in this fashion.

Another detail we will largely ignore in this paper is the fact that efficient algorithms for these problems can only approximately solve the dual norm maximization problem. However, we explain in Appendix B.2 why this does not affect us in any meaningful way. Thus, for the rest of the paper we will assume we have access to the exact maximizer, and the exact value of the norm.

We now show how the above concentration inequalities allow us to derive very strong concentration results for the dual norm maximization problem for $\mathcal{U}_{k}$ and $\mathcal{X}_{k}$. Conceptually, we view these concentration results as being the major distinction between sparse estimation and non-sparse estimation tasks. Indeed, these results are crucial for adapting the convex programming framework for robust estimation to sparse estimation tasks. Additionally, they allow us to give an easy proof that the $L_{1}$ relaxation works for sparse PCA.

Fix $\varepsilon,\delta>0$. Let $X_{1},\ldots,X_{n}\sim\operatorname{\mathcal{N}}(0,I)$, where

Then $\|\frac{1}{n}\sum_{i=1}^{n}X_{i}\|^{*}_{\mathcal{U}_{k}}\leq\varepsilon$.

Fix a set of $k$ coordinates, and let $S$ be the set of unit vectors supported on these $k$ coordinates. By Fact 3.2 and a net argument, one can show that for all $\delta$, given $n=\Omega\left(\frac{k+\log 1/\delta}{\varepsilon^{2}}\right)$, we have that

with probability $1-\delta$. The result then follows by setting $\delta^{\prime}=\binom{d}{k}^{-1}\delta$ and union bounding over all sets of $k$ coordinates. ∎

The second concentration bound, which bounds deviation in $\mathcal{X}_{k}$ norm, uses ideas which are similar at a high level, but requires a bit more technical work.

Fix $\varepsilon,\delta>0$. Let $X_{1},\ldots X_{n}\sim\operatorname{\mathcal{N}}(0,I)$, where

Let us first introduce the following definition.

Let $n=O\left(\frac{\min(d,k^{2})+\log\binom{d^{2}}{k^{2}}+\log 1/\delta}{\varepsilon^{2}}\right)$. Then, with probability $1-\delta$, the following holds:

We will also require the following structural lemma.

where each $Y_{i}$ is symmetric, have $\sum_{i=1}^{O(n^{2}/k^{2})}\|Y_{i}\|_{F}\leq 4$, and each $Y_{i}$ is $k^{2}$-sparse.

However, as written it is not clear that the $Y_{i}$’s must be symmetric, and indeed they do not have to be. The only real condition we needed was that the $Y_{i}$’s (1) had disjoint support, (2) summed to $X$, (3) are each $\Theta(k^{2})$ sparse (except potentially the last one), and (4) the largest entry of $Y_{i+1}$ is bounded by the smallest entry of $Y_{i}$. It should be clear that this can be done while respecting symmetry by doubling the number of $Y_{i}$, which also at most doubles the bound in the sum of the Frobenius norms. We omit the details for simplicity. ∎

where each $Y_{i}$ is symmetric, have $\sum_{i=1}^{O(d^{2}/k^{2})}\|Y_{i}\|_{F}\leq 4$, and each $Y_{i}$ is $k^{2}$-sparse. Thus,

where (a) follows since each $Y_{i}/\|Y_{i}\|_{F}$ satisfies the conditions in (4), and (b) follows from the bound on the sum of the Frobenius norms of the $Y_{i}$. ∎

We will require the following concentration inequalities for weighted sums of Gaussians, where the weights come from $S_{n,\varepsilon}$, as these objects will naturally arise in our algorithms. These bounds follow by applying the above bounds, then carefully union bounding over all choices of possible subsets of $\binom{n}{\varepsilon n}$ subsets. We need to be careful here since the number of things we are union bounding over increases as $n$ increases. We include the proofs in Appendix C.

Fix $\varepsilon\leq 1/2$ and $\delta\leq 1$, and fix $k\leq d$. There is a $\eta_{1}=O(\varepsilon\sqrt{\log 1/\varepsilon})$ so that for any $\eta>\eta_{1}$, if $X_{1},\ldots,X_{n}\sim\operatorname{\mathcal{N}}(0,I)$ and $n=\Omega\left(\frac{\min(d,k^{2})+\log\binom{d^{2}}{k^{2}}+\log 1/\delta}{\eta^{2}}\right)$, then

Fix $\varepsilon\leq 1/2$ and $\delta\leq 1$, and fix $k\leq d$. There is a $\eta=O(\varepsilon\sqrt{\log 1/\varepsilon})$ so that if $X_{1},\ldots,X_{n}\sim\operatorname{\mathcal{N}}(0,I)$ and $n=\Omega\left(\frac{\min(d,k^{2})+\log\binom{d^{2}}{k^{2}}+\log 1/\delta}{\eta^{2}}\right)$, then we have

A robust algorithm for robust sparse mean estimation

This section is dedicated to the description of an algorithm RecoverRobustSMean for robustly learning Gaussian sequence models, and the proof of the following theorem:

Fix $\varepsilon,\tau>0$. Let $\eta=O(\varepsilon\sqrt{\log 1/\varepsilon})$. Given an $\varepsilon$-corrupted set of samples of size $n$ from $\operatorname{\mathcal{N}}(\mu,I)$, where $\mu$ is $k$-sparse

then RecoverRobustSMean outputs a $\widehat{\mu}$ so that with probability $1-\tau,$ we have $\|\widehat{\mu}-\mu\|_{2}\leq O(\eta)$.

Our algorithm builds upon the convex programming framework developed in [DKK+16]. Roughly speaking, the algorithm proceeds as follows. First, it does a simple naive pruning step to remove all points which are more than roughly $\Omega(\sqrt{d})$ away from the mean. Then, for an appropriate choice of $\delta$, it will attempt to (approximately) find a point within the following convex set:

The main difficulty with finding a point in $C_{\tau}$ is that $\mu$ is unknown. A key insight of [DKK+16] is that it suffices to create an (approximate) separation oracle for the feasible set, as then we may use classical convex optimization algorithms (i.e. ellipsoid or cutting plane methods) to find a feasible point. In their setting (for a different $C_{\tau}$), it turns out that a simple spectral algorithm suffices to give such a separation oracle.

Our main contribution is the design of separation oracle for $C_{\tau}$, which requires more sophisticated techniques. In particular, we will ideas developed in analogy to hard thresholding and SDPs similar to those developed for sparse PCA to design such an oracle.

Throughout this section, we will condition on the following three deterministic events occurring:

2 The separation oracle

Our main result in this section is the description of a polynomial time algorithm RobustSMeanOracle and the proof of the following theorem of its correctness:

Fix $\varepsilon>0$ sufficiently small. Suppose that (7) and (8) hold. Let $w^{*}$ denote the set of weights which are uniform over the uncorrupted points. Then, there is a constant $1\leq c\leq 21$ so that RobustSMeanOracle satisfies:

(Completeness) If $w=w^{*}$, RobustSMeanOracle outputs “YES”.

Plugging these guarantees into an ellipsoid (or cutting-plane) method, we obtain the following:

Fix $\varepsilon>0$ sufficiently small. Suppose that (7) and (8) hold. There is an algorithm ApproxRecoverRobustSMean which queries RobustSMeanOracle at most $\operatorname{poly}(d,1/\varepsilon,\log 1/\delta)$ times, and so runs in time $\operatorname{poly}(d,1/\varepsilon,1/\delta)$ which outputs a $w^{\prime}$ so that $\|w-w^{\prime}\|_{\infty}\leq\varepsilon/(n\sqrt{d\log n/\delta})$, for some $w\in C_{c\tau}$.

Our separation oracle, formally described in Algorithm 1, proceeds as follows. Given $w\in S_{n,\varepsilon}$, it forms $\widehat{\mu}=\|\widehat{\mu}^{\prime}\|^{*}_{\mathcal{U}_{k}}\cdot d_{\mathcal{U}_{k}}(\widehat{\mu}^{\prime})$, where $\widehat{\mu}=\sum w_{i}X_{i}$. It then forms the matrix $\widehat{\Sigma}=\sum w_{i}(X_{i}-\widehat{\mu})(X_{i}-\widehat{\mu})^{T}$, and computes $A=d_{\mathcal{X}_{k}}(\widehat{\Sigma})$. The algorithm then checks if $\left|\langle A,\widehat{\Sigma}\rangle\right|>C$ for appropriately chosen threshold $C$. If it does not, the algorithm outputs “YES”. Otherwise, the algorithm outputs a separating hyperplane given by this matrix $A$.

We will require the following two lemmata:

Let $\omega_{1},\ldots,\omega_{m}$ be a set of non-negative weights that sum to 1. Let $a_{1},\ldots,a_{m}$ be any sequence of scalars. Then

Let $v=d_{\mathcal{U}_{k}}(u)$. Then since $vv^{T}\in\mathcal{X}_{k}$, we have that $(\|uu^{T}\|_{\mathcal{X}_{k}}^{*})\geq\langle vv^{T},uu^{T}\rangle=\langle u,v\rangle^{2}=(\|u\|_{\mathcal{U}_{k}}^{*})^{2}$. This proves the first inequality.

To prove the other inequality, we first prove the intermediate claim that $\sup_{M\in\mathcal{Y}_{k^{2}}}u^{T}Mu\leq(\|u\|_{\mathcal{U}_{k}}^{*})^{2}$, where $\mathcal{Y}_{k^{2}}$ is the set of symmetric matrices $M$ with at most $k^{2}$-non-zeroes satisfying $\|M\|_{F}=1$. Indeed, fix any $M\in\mathcal{Y}_{k}$. Let $S\subseteq[n]$ be the set of non-zeroes of $d_{\mathcal{U}_{k}}(u)$. This is exactly the set of the $k$ largest elements in $u$, sorted by absolute value. Let $P$ be the symmetric sparsity pattern respected by $M$. Fix an arbitrary bijection $\phi:P\setminus(S\times S)\to(S\times S)\setminus P$, and let $M^{\prime}$ be the following matrix:

Then we claim that $u^{T}Mu\leq u^{T}M^{\prime}u$. Indeed, we have

We can now prove the original lemma. By Lemma 4.4 we may write $A=\sum_{i=1}^{O(n^{2}/k^{2})}Y_{i}$ where each $Y_{i}$ is symmetric, $k^{2}$-sparse, and have $\sum_{i=1}^{O(n^{2}/k^{2})}\|Y_{i}\|_{F}\leq 4$. We therefore have

as claimed, where the second line follows from the arguments above. ∎

Let $w\in S_{n,\varepsilon}$, and let $\tau\geq\eta$. Assuming (7) and (8) hold, if $\left\|\sum_{i=1}^{n}w_{i}Y_{i}\right\|_{\mathcal{U}_{k}}^{*}\geq 3\tau$, then $\left\|\sum_{i=1}^{n}w_{i}Y_{i}Y_{i}^{T}-I\right\|_{\mathcal{X}_{k}}^{*}\geq\frac{\tau^{2}}{\varepsilon}$.

Observe that the $w_{i}/w^{b}$ are a set of non-negative weights summing to $1$. Hence, by Lemma 5.4, we have

Let $A=uu^{T}$. Observe that $A\in\mathcal{X}_{k}$. Then the above inequality is equivalent to the statement that

and together these two inequalities imply that

as claimed. The final inequality follows from the definition of $\eta$, and since $4>2$. ∎

Completeness follows from (8). We will now show soundness. Suppose $w\not\in C_{21\eta}$. We wish to show that we will output a separating hyperplane. From the description of the algorithm, this is equivalent to showing that $\|\widehat{\Sigma}-I\|_{\mathcal{X}_{k}}\geq 20\eta$. Let $\widehat{\mu}=\sum_{i=1}^{n}w_{i}X_{i}$, and let $\Delta=\mu-\widehat{\mu}$. By elementary manipulations, we may write

where (a) follows since $\sum_{i=1}^{n}w_{i}Y_{i}=\Delta$ by definition, (b) follows from a triangle inequality, and (c) follows from Lemma 5.5. If $\|\Delta\|_{\mathcal{U}_{k}}\leq\sqrt{\eta/2}$, then the RHS is at least $21\eta$ since the second term is at most $\eta$, and the first term is at least $21\eta$ since we assume that $w\not\in C_{21\eta}$. Conversely, if $\|\Delta\|_{\mathcal{U}_{k}}\geq\sqrt{\eta/2}$, then by Proposition 5.6, we have $\|\sum_{i=1}^{n}w_{i}Y_{i}Y_{i}-I\|_{\mathcal{X}_{k}}\geq\|\Delta\|_{\mathcal{X}_{k}}^{2}/(6\varepsilon)>48\|\Delta\|_{\mathcal{X}_{k}}^{2}$ as long as $\varepsilon\leq 1/288$. This implies that the RHS is at least $40\|\Delta\|_{\mathcal{X}_{k}^{2}}\geq 20\eta$, as claimed.

Hence by the triangle inequality and Lemma 5.5, we have

If $\|\Delta\|_{\mathcal{U}_{k}}^{*}\leq\sqrt{\eta/2}$, then this follows since the quantity on the RHS is at least $20\eta$ by assumption, and the quantity on the LHS is at most $17\eta$ by (10). If $\|\Delta\|_{\mathcal{U}_{k}}^{*}\geq\sqrt{\eta/2}$, then by Proposition 5.6, the RHS of (11) is at least $\left(\|\Delta\|_{\mathcal{U}_{k}}^{*}\right)^{2}/(3\varepsilon)$, which dominates the LHS as long as $\|\Delta\|_{\mathcal{U}_{k}}^{*}\geq\eta$ and $\varepsilon\leq 1/288$, which completes the proof. ∎

3 Putting it all together

We now have the ingredients to prove our main theorem. Given what we have, our full algorithm RecoverRobustSMean is straightforward: first run NaivePrune, then run ApproxRecoverRobustSMean on the pruned points to output some set of weights $w$. We then output $\|\widehat{\mu}\|_{\mathcal{U}_{k}}d_{\mathcal{U}_{k}}(\widehat{\mu})$. The algorithm is formally defined in Algorithm 2.

Let us condition on the event that (6), (7), and (8) all hold simultaneously. As previously mentioned, when $n=\Omega\left(\frac{\min(k^{2},d)+\log\binom{k^{2}}{d^{2}}+\log 1/\delta}{\eta^{2}}\right)$ these events simultaneously happen with probability at least $1-O(\delta)$. For simplicity of exposition, let us assume that NaivePrune does not remove any points. This is okay since if it succeeds, it never removes any good points, so if it removes any points, it can only help us. Moreover, since it succeeds, we know that $\|X_{i}-\mu\|_{2}\leq O(\sqrt{d\log(n/\delta)})$ for all $i\in[n]$. By Corollary 5.3, we know that there is some $w\in C_{21\eta}$ so that $\|w-w^{\prime}\|_{\infty}\leq\varepsilon/(n\sqrt{d\log n/\delta})$. We have

by Proposition 5.6. We now show that this implies that if we let $\mu^{\prime}=\|\widehat{\mu}\|^{*}_{\mathcal{U}_{k}}d_{\mathcal{U}_{k}}(\widehat{\mu})$, then $\|\mu^{\prime}-\mu\|_{2}\leq O(\eta)$. Let $S$ be the support of $\mu^{\prime}$, and let $T$ be the support of $\mu$. Then we have

Observe that $\sum_{i\in S\cap T}(\mu^{\prime}_{i}-\mu_{i})^{2}+\sum_{i\in S\setminus T}(\mu_{i}^{\prime})^{2}\leq\left(\|\widehat{\mu}-\mu\|_{\mathcal{U}_{k}}^{*}\right)^{2}$, since $\mu$ was originally nonzero on the entries in $S\setminus T$. Moreover, for all $i\in T\setminus S$ and $j\in S\setminus T$, we have $(\mu^{\prime}_{i})^{2}\leq(\mu^{\prime}_{j})^{2}$. Thus we have

Therefore we have $\|\mu^{\prime}-\mu\|_{2}^{2}\leq 3\left(\|\widehat{\mu}-\mu\|_{\mathcal{U}_{k}}^{*}\right)^{2}$, which implies that $\|\mu^{\prime}-\mu\|_{2}\leq O(\eta)$, as claimed. ∎

An algorithm for robust sparse PCA detection

In this section, we give an efficient algorithm for detecting a spiked covariance matrix in the presence of adversarial noise. Our algorithm is fairly straightforward: we ask for the set of weights $w\in S_{n,\varepsilon}$ so that the empirical second moment with these weights has minimal deviation from the identity in the dual $\mathcal{X}_{k}$ norm. We may write this as a convex program. Then, we check the value of the optimal solution of this convex program. If this value is small, then we say it is $\operatorname{\mathcal{N}}(0,I)$. if this value is large, then we say it is $\operatorname{\mathcal{N}}(0,I+\rho vv^{T})$. We refer to the former as Case 1 and the latter as Case 2. The formal description of this algorithm is given in Algorithm.

We first show that the algorithm presented above can be efficiently implemented. Indeed, one can show that by taking the dual of the SDP defining the $\|\cdot\|^{*}_{\mathcal{X}_{k}}$ norm, this problem can be re-written as an SDP with (up to constant factor blowups) the same number of constraints and variables, and therefore we may solve it using traditional SDP solver techniques.

Alternatively, one may observe that to optimize Algorithm 4 via ellipsoid or cutting plane methods, it suffices to, given $w\in S_{n,\varepsilon}$, produce a separating hyperplane for the constraint (12). This is precisely what dual norm maximization allows us to do efficiently. It is straightforward to show that the volume of $S_{n,\varepsilon}\times\mathcal{X}_{k}$ is at most exponential in the relevant parameters. Therefore, by the classical theory of convex optimization, (see e.g. [CITE]), for any $\xi$, we may find a solution $w^{\prime}$ and $\gamma^{\prime}$ so that $\|w^{\prime}-w^{*}\|_{\infty}\leq\xi$ and $\gamma^{\prime}$ so that $|\gamma-\gamma^{\prime}|<\xi$ for some exact minimizer $w^{*}$, where $\gamma$ is the true value of the solution, in time $\operatorname{poly}(d,n,1/\varepsilon,\log 1/\xi)$,

As mentioned in Section B.2, neither approach will in general give exact solutions, however, both can achieve inverse polynomial accuracy in the parameters in polynomial time. We will ignore these issues of numerical precision throughout the remainder of this section, and assume we work with exact $\gamma$.

Observe that in general it may be problematic that we don’t have exact access to the minimizer $w^{*}$, since some of the $X_{i}$ may be unboundedly large (in particular, if it’s corrupted) in norm. However, we only use information about $\gamma$. Since $\gamma$ lives within a bounded range, and our analysis is robust to small changes to $\gamma$, these numerical issues do not change anything in the analysis.

2 Proof of Theorem 2.2

We now show that Algorithm 4 provides the guarantees required for Theorem 2.2. We first show that if we are in Case 1, then $\gamma$ is small:

Let $\rho,\delta>0$. Let $\varepsilon,\eta$ be as in Theorem 2.2. Let $X_{1},\ldots,X_{n}$ be an $\varepsilon$-corrupted set of samples from $\operatorname{\mathcal{N}}(0,I)$ of size $n$, where $n$ is as in Theorem 2.2. Then, with probability $1-\delta$, we have $\gamma\leq\rho/2$.

Let $w$ be the uniform weights over the uncorrupted points. Then it from Theorem 4.2 that $\|\sum_{w}\sum_{i=1}^{n}w_{i}(X_{i}X_{i}^{T}-I)\|^{*}_{\mathcal{X}_{k}}\leq O(\eta)$ with probability $1-\delta$. Since $w\in S_{n,\varepsilon}$, this immediately implies that $\gamma\leq O(\rho)$. By setting constants appropriately, we obtain the desired guarantee. ∎

We now show that if we are in Case 2, then $\gamma$ must be large:

Let $\rho,\delta>0$. Let $\varepsilon,\eta,n$ be as in Theorem 2.2. Let $X_{1},\ldots,X_{n}$ be an $\varepsilon$-corrupted set of samples from $\operatorname{\mathcal{N}}(0,I)$ of size $n$. Then, with probability $1-\delta$, we have $\gamma\geq(1-\varepsilon)\rho-(2+\rho)\eta$. In particular, for $\varepsilon$ sufficiently small, and $\eta=O(\rho)$, we have that $\gamma>\rho/2$.

It thus suffices to show that $|v^{T}\Sigma^{1/2}N\Sigma^{1/2}v|<(1+\rho)\eta$. Since $v$ is an eigenvector for $\Sigma$ with eigenvalue $1+\rho$, we have that $\Sigma^{1/2}v=\sqrt{\rho+1}\cdot v$ and thus

Lemmas 6.1 and 6.2 together imply the correctness of DetectRobustSPCA and Theorem 2.2.

An algorithm for robust sparse PCA recovery

In this section, we prove Theorem 2.3. We give some intuition here. Perhaps the first naive try would be to simply run the same SDP in (12), and hope that the dual norm maximizer gives you enough information to recover the hidden spike. This would more or less correspond to the simplest modification SDP of the sparse PCA in the non-robust setting that one could hope gives non-trivial information in this setting. However, this cannot work, for the following straightforward reason: the value of the SDP is always at least $O(\rho)$, as we argued in Section 6. Therefore, the noise can pretend to be some other sparse vector $u$ orthogonal to $v$, so that the covariance with noise looks like $w^{g}(I+\rho vv^{T})+w^{g}\rho uu^{T}$, so that the value of the SDP can be minimized with the uniform set of weights. Then it is easily verified that both $vv^{T}$ and $uu^{T}$ are dual norm maximizers, and so the dual norm maximizer does not uniquely determine $v$.

To circumvent this, we simply add an additional slack variable to the SDP, which is an additional matrix in $\mathcal{X}_{k}$, which we use to try to maximally explain away the rank-one part of $I+\rho vv^{T}$. This forces the value of the SDP to be very small, which allows us to show that the slack variable actually captures $v$.

Our algorithms and analyses will make crucial use of the following convex set, which is a further relaxation of $\mathcal{X}_{k}$:

Our algorithm, given formally in Algorithm 4, will be the following. We solve a convex program which simultaneously chooses a weights in $S_{n,\varepsilon}$ and a matrix $A\in\mathcal{W}_{k}$ to minimize the $\mathcal{W}_{k}$ distance between the sample covariance with these weights, and $A$. Our output is then just the top eigenvector of $A$.

This algorithm can be run efficiently for the same reasons as explained for DetectRobustSPCA. For the rest of the section we will assume that we have an exact solution for this problem. As before, we only use information about $A$, and since $A$ comes from a bounded space, and our analysis is robust to small perturbations in $A$, this does not change anything.

2 More concentration bounds

Before we can prove correctness of our algorithm, we require a couple of concentration inequalities for the set $\mathcal{W}_{k}$.

Fix $\varepsilon,\delta>0$. Let $X_{1},\ldots,X_{n}\sim\operatorname{\mathcal{N}}(0,I)$, where $n$ is as in Theorem 4.2. Then with probability $1-\delta$

Let $\widehat{\Sigma}$ denote the empirical covariance. Observe that $\mathcal{W}_{k}\subseteq\bigcup_{i=0}^{\infty}2^{-i}\mathcal{X}_{2^{i+1}k}$. Moreover, for any $i$, by Theorem 4.2, if we take

then $|\langle M,\widehat{\Sigma}\rangle|\leq\varepsilon$ for all $M\in 2^{-i}\mathcal{X}_{2^{i+1}k}$ with probability $1-\delta/2$. In particular, if we take

samples, then for any $i$, we have $|\langle M,\widehat{\Sigma}\rangle|\leq\varepsilon$ for all $M\in 2^{-1}\mathcal{X}_{2^{i+1}k}$ with probability at least $1-\delta^{2^{2i}}/2$. By a union bound over all these events, since $\sum_{i=0}^{\infty}\delta^{2^{2i}}\leq 2\delta$, we conclude that if we take $n$ to be as above, then $|\langle M,\widehat{\Sigma}\rangle|\leq\varepsilon$ for all $M\in\bigcup_{i=0}^{\infty}2^{-i}\mathcal{X}_{2^{i+1}k}$ with probability $1-\delta$. Since $\mathcal{W}_{k}$ is contained in this set, this implies that $\|\widehat{\Sigma}-\Sigma\|^{*}_{\mathcal{W}_{k}}\leq O(\varepsilon)$ with probability at least $1-\delta$, as claimed. ∎

By the same techniques as in the proofs of Theorems 4.5 and 4.6, we can show the following bound. Because of this, we omit the proof for conciseness.

Fix $\varepsilon,\delta>0$. Let $X_{1},\ldots,X_{n}\sim\operatorname{\mathcal{N}}(0,I)$ where $n$ is as in Theorem 4.6. Then there is an $\eta=O(\varepsilon\sqrt{\log 1/\varepsilon})$ so that

3 Proof of Theorem 2.3

In the rest of this section we will condition on the following deterministic event happening:

where $\eta=O(\varepsilon\log 1/\varepsilon)$. By Corollary 7.2, this holds if we take

The rest of this section is dedicated to the proof of the following theorem, which immediately implies Theorem 2.3.

Fix $\varepsilon,\delta,$ and let $\eta$ be as in (14). Assume that (14) holds. Let $\widehat{v}$ be the output of $\textsc{RecoveryRobustSPCA}(X_{1},\ldots,X_{n},\varepsilon,\delta,\rho)$. Then $L(\widehat{v},v)\leq O(\sqrt{(1+\rho)\eta/\rho})$.

Our proof proceeds in a couple of steps. Let $\Sigma=I+\rho vv^{T}$ denote the true covariance. We first need the following, technical lemma:

Let $M\in\mathcal{W}_{k}$. Then $\Sigma^{1/2}M\Sigma^{1/2}\in(1+\rho)\mathcal{W}_{k}$.

Clearly, $\Sigma^{1/2}M\Sigma^{1/2}\succeq 0$. Moreover, since $\Sigma^{1/2}=I+(\sqrt{1+\rho}-1)vv^{T}$, we have that the maximum value of any element of $\Sigma^{1/2}$ is upper bounded by $\sqrt{1+\rho}$. Thus, we have $\|\Sigma^{1/2}M\Sigma^{1/2}\|_{1}\leq(1+\rho)\|M\|_{1}$. We also have

since $\|M\|\leq 1$. Thus $\Sigma^{1/2}M\Sigma^{1/2}\in(1+\rho)\mathcal{W}_{k}$, as claimed. ∎

Let $w^{*},A^{*}$ be the output of our algorithm. We first claim that the value of the optimal solution is quite small:

Indeed, if we let $w$ be the uniform set of weights over the good points, and we let $A=vv^{T}$, then by (14), we have

where $\|N\|^{*}_{\mathcal{X}_{k}}\leq\eta$, and $\Sigma=I+\rho vv^{T}$. Thus we have that

We now show that this implies the following:

Now, since $vv^{T}\in\mathcal{W}_{k}$, the above implies that

which by a further triangle inequality implies that

Since $0\leq v^{T}A^{*}v\leq 1$ (since $A\in\mathcal{X}_{k}$) and $B$ is PSD, this implies that in fact, we have

Hence $v^{T}A^{*}v\geq 1-(2+3\rho)\eta/\rho$, as claimed. ∎

which by a further triangle inequality implies that

We now show this implies the following intermediate result:

By Lemma 7.6, we have that $v^{T}A^{*}v=\lambda_{1}(v^{T}u)^{2}+v^{T}A_{1}v\geq 1-\gamma$. In particular, this implies that $(v^{T}u)^{2}\geq(1-2\gamma)/\lambda_{1}\geq 1-3\gamma$, since $1-\gamma\leq\lambda\leq 1$. ∎

We now wish to control the spectrum of $B$. For any subsets $S,T\subseteq[d]$, and for any vector $x$ and any matrix $M$, let $x_{S}$ denote $x$ restricted to $S$ and $M_{S,T}$ denote the matrix restricted to the rows in $S$ and the columns in $T$. Let $I$ be the support of $u$, and let $J$ be the support of the largest $k$ elements of $v$.

Observe that the condition (15) immediately implies that

Lemma 7.8 and (16) together imply that $\|v_{I}v_{I}^{T}-u_{I}u_{I}^{T}\|\leq O(\gamma)$. The desired bound then follows from a reverse triangle inequality. ∎

We now show this implies a bound on $B_{J\setminus I,J\setminus I}$:

$\|B_{J\setminus I,J\setminus I}\|\leq O(\rho\gamma)$.

Suppose $\|B_{J\setminus I,J\setminus I}\|\geq C\gamma$ for some sufficiently large $C$. Since $u$ is zero on $J\setminus I$, (15) implies that

for some universal $c$. By a triangle inequality, this implies that $\|v_{J\setminus I}\|_{2}^{2}=\|v_{J\setminus I}v_{J\setminus I}^{T}\|\geq(C-c)\gamma$. Since $v$ is a unit vector, this implies that $\|v_{I}\|_{2}^{2}\leq 1-(C-c)\gamma$, which for a sufficiently large $C$, contradicts Corollary 7.9. ∎

We now invoke the following general fact about PSD matrices:

Suppose $M$ is a PSD matrix, written in block form as

Suppose furthermore that $\|C\|\leq\xi$ and $\|E\|\leq\xi$. Then $\|M\|\leq O(\xi)$.

It is easy to see that $\|M\|\leq O(\max(\|C\|,\|D\|,\|E\|))$. Thus it suffices to bound the largest singular value of $D$. For any vectors $\phi,\psi$ with appropriate dimension, we have that

which immediately implies that the largest singular value of $D$ is at most $(\|A\|+\|B\|)/2$, which implies the claim. ∎

Therefore, Lemmas 7.8 and 7.10 together imply:

$\|v_{I\cup J}v_{I\cup J}^{T}-u_{I\cup J}u_{I\cup J}^{T}\|\leq O(\gamma)\;.$

Observe (15) immediately implies that $\|\rho(v_{I\cup J}v_{I\cup J}^{T}-u_{I\cup J}u_{I\cup J}^{T})+B_{I\cup J,I\cup J}\|\leq O(\rho\gamma)$, since $|I\cup J|\leq 2k$. Moreover, Lemmas 7.8 and 7.10 with Lemma 7.11 imply that $\|B_{I\cup J,I\cup J}\|\leq O(\rho\gamma)$, which immediately implies the statement by a triangle inequality. ∎

Finally, we show this implies $\|vv^{T}-u_{J}u_{J}^{T}\|\leq O(\gamma)$, which is equivalent to the theorem.

We will in fact show the slightly stronger statement, that $\|uu^{T}-v_{J}v_{J}^{T}\|_{F}\leq O(\gamma)$. Observe that since $uu^{T}-vv^{T}$ is rank 2, Corollary 7.12 implies that $\|v_{I\cup J}v_{I\cup J}^{T}-u_{I\cup J}u_{I\cup J}^{T}\|_{F}\leq O(\gamma)$, since for rank two matrices, the spectral and Frobenius norm are off by a constant factor. We have

by Corollary 7.12. Moreover, we have that

since $J\times J$ contains the $k^{2}$ largest entries of $uu^{T}$. This completes the proof. ∎

Acknowledgements

The author would like to thank Ankur Moitra for helpful advice throughout the project, and Michael Cohen for some surprisinglyIs it really surprising though? useful conversations.

References

Appendix A Information theoretic estimators for robust sparse estimation

This section is dedicated to the proofs of the following two facts:

there is an (inefficient) algorithm which outputs $\widehat{\mu}$ so that with probability $1-\delta$, we have $\|\mu-\widehat{\mu}\|_{2}\leq O(\varepsilon)$. Moreover, up to logarithmic factors, this rate is optimal.

Fix $\rho,\delta>0$. Suppose that $\rho=O(1)$. Then, there exist universal constants $c,C$ so that: (a) if $\varepsilon\leq c\rho$, and we are given a $\varepsilon$-corrupted set of samples from either $\operatorname{\mathcal{N}}(0,I)$ or $\operatorname{\mathcal{N}}(0,I+\rho vv^{T})$ for some $k$-sparse unit vector $v$ of size

then there is an (inefficient) algorithm which succeeds with probability $1-\delta$ for the detection problem. Moreover, if $\varepsilon\geq C\rho$, then no algorithm succeeds with probability greater than $1/2$, and this statistical rate is optimal.

The rates in Facts A.1 and A.2 are already known to be optimal (up to log factors) without noise. Thus in this section we focus on proving the upper bounds, and the lower bounds on error.

Our techniques for proving the upper bounds go through the technique of agnostic hypothesis selection via tournaments. Specifically, we use the following lemma:

Let $\mathcal{M}_{A}$ be the set of distributions $\{\operatorname{\mathcal{N}}(\mu^{\prime},I)\}$, where $\mu^{\prime}$ ranges over the set of $k$-sparse vectors so that each coordinate of $\mu^{\prime}$ is an integer multiple of $\varepsilon/(10\sqrt{d})$, and so that $\|\mu^{\prime}-\mu\|_{2}\leq A$. We then have:

There exists a $\operatorname{\mathcal{N}}(\mu^{\prime},I)=D\in\mathcal{M}_{A}$ so that $\|\mu-\mu^{\prime}\|_{2}\leq O(\varepsilon)$. Moreover, $|\mathcal{M}_{A}|\leq\binom{d}{k}\cdot(10A\sqrt{d}/\varepsilon)^{k}$.

The first claim is straightforward. We now prove the second claim. For each possible set of $k$ coordinates, there are at most $(10A\sqrt{d}/\varepsilon)^{k}$ vectors supported on those $k$ coordinates with each coordinate being an integer multiple of $\varepsilon/(10\sqrt{d})$ with distance at most $A$ from any fixed vector. Enumerating over all $\binom{d}{k}$ possible choices of $k$ coordinates yields the desired answer. ∎

The estimator is given as follows: first, run $\textsc{NaivePrune}(X_{1},\ldots,X_{n},\delta)$ to output some $\mu_{0}$ so that with probability $1-\delta$, we have $\|\mu_{0}-\mu\|_{2}\leq O(\sqrt{d\log n/\delta})$. Round each coordinate of $\mu_{0}$ so that it is an integer multiple of $\varepsilon/(10\sqrt{d})$. Then, output the set of distributions $\mathcal{M^{\prime}}=\{\operatorname{\mathcal{N}}(\mu^{\prime\prime},I)\}$, where $\mu^{\prime\prime}$ is any $k$-sparse vector with each coordinate being an integer multiple of $\varepsilon/(10\sqrt{d})$, with $\|\alpha\|_{2}\leq O(\sqrt{d\log n/\delta})$. With probability $1-\delta$, we have $\mathcal{M^{\prime}}\subseteq\mathcal{M}_{O(\sqrt{d\log n/\delta})}$. By Claim A.6, applying Lemma A.5 to this set of distributions yields that we will select, with probability $1-\delta$, a $\mu^{\prime}$ so that $\|\mu-\mu^{\prime}\|_{2}\leq O(\varepsilon)$. By Claim A.6, this requires

samples, which simplifies to the desired bound, as claimed.

A.2 Proof of Upper Bound in Fact A.2

Our detection algorithm is given as follows. We let $\mathcal{N}$ be an $O(1)$-net over all $k$-sparse unit vectors, and we apply Lemma A.5 to the set $\{\operatorname{\mathcal{N}}(0,I+\rho uu^{T})\}_{u\in\mathcal{N}}\cup\{\operatorname{\mathcal{N}}(0,I)\}$. Clearly, we have:

By Fact A.4 and the guarantees of Lemma A.5, by an appropriate setting of parameters, if we have

samples, then with probability $1-\delta$ we will output $\operatorname{\mathcal{N}}(0,I)$ if and only if the true model is $\operatorname{\mathcal{N}}(0,I)$. This proves the upper bound.

Appendix B Omitted Details from Section 4

In this section we will briefly review well-known non-robust algorithms for sparse mean recovery and for sparse PCA, and write them using our language.

Recall that in the (non-robust) sparse mean estimation problem, one is given samples $X_{1},\ldots,X_{n}\sim\operatorname{\mathcal{N}}(\mu,I)$ where $\mu$ is $k$-sparse. The goal is then to recover $\mu$. It turns out the simple thresholding algorithm ThresholdMean given in Algorithm 5 suffices for recovery:

The correctness of this algorithm follows from the following folklore result, whose proof we shall omit for conciseness:

Fix $\varepsilon,\delta>0$, and let $X_{1},\ldots,X_{n}$ be samples from $\operatorname{\mathcal{N}}(\mu,I)$, where $\mu$ is $k$-sparse and

Then, with probability $1-\delta$, if $\widehat{\mu}^{\prime}$ is the output of ThresholdMean, we have $\|\widehat{\mu}^{\prime}-\widehat{\mu}\|_{2}\leq\varepsilon$.

To write this in our language, observe that

where $\widehat{\mu}=\frac{1}{n}\sum_{i=1}^{n}X_{i}$.

In various scenarios, including recovery of a spiked covariance, one may envision the need to take $k$-sparse eigenvalues a matrix $A$, that is, vectors which solve the following non-convex optimization problem:

However, this problem is non-convex and cannot by solved efficiently. This motivates the following SDP relaxation of (18): First, one rewrites the problem as

The work of [dEGJL07] shows that this indeed detects the presence of a spike (but at an information theoretically suboptimal rate).

Finally, by definition, for any PSD matrix $A$, if $X$ is the solution to (20) with input $A$, we have $X=d_{\mathcal{X}_{k}}(A)$.

B.2 Numerical precision

In general, we cannot find closed form solutions for $d_{\mathcal{X}_{k}}(A)$ in finite time. However, it is well-known that we can find these to very high numerical precision in polynomial time. For instance, using the ellipsoid method, we can find an $M^{\prime}$ so that $\|M^{\prime}-d_{\mathcal{X}_{k}}(A)\|_{\infty}\leq\varepsilon$ in time $\operatorname{poly}(d,\log 1/\varepsilon)$. It is readily verified that if we set $\varepsilon^{\prime}=\operatorname{poly}(\varepsilon,1/d)$ then the numerical precision of the answer will not effect any of the calculations we make further on. Thus for simplicity of exposition we will assume throughout the paper that given any $A$, we can find $d_{\mathcal{X}_{k}}(A)$ exactly in polynomial time.

Appendix C Omitted Proofs from Section 4

Fix $n$ as in Theorem 4.5, and let $\delta_{1}=\binom{n}{\varepsilon n}^{-1}\delta$. By convexity of $S_{n,\varepsilon}$ and the objective function, it suffices to show that with probability $1-\delta$, the following holds:

By Corollary 4.1, this occurs with probability $1-O(\delta)$.

Fix any $I\subseteq[n]$ so that $|I|=(1-\varepsilon)n$. By Corollary 4.1 applied to $I^{c}$, we have that there is some universal constant $C$ so that as long as

then with probability $1-\delta^{\prime}$,

Since $\log\binom{n}{\varepsilon n}=\Theta(n\varepsilon\log 1/\varepsilon)$, (22) is equivalent to the condition that

Let $\alpha=O(\sqrt{\log 1/\varepsilon})$. By our choice of $\eta$, we have that $0\leq\varepsilon-\frac{\varepsilon\log 1/\varepsilon}{\eta^{2}}\leq\varepsilon/(2C)$, and by an appropriate setting of constants, since by our choice of $n$ we have

we have that (23) holds with probability $1-\delta^{\prime}$. Thus by a union bound over all $\binom{n}{\varepsilon n}$ choices of $I$ so that $|I|=(1-\varepsilon)n$, we have that except with probability $1-\delta$, we have that (23) holds simultaneously for all $I$ with $|I|=(1-\varepsilon)n$. The desire result then follows from this and (21), and a union bound. ∎