Consistent Estimation for PCA and Sparse Regression with Oblivious Outliers
Year: 2021

Abstract:
We develop machinery to design efficiently computable and consistent estimators, achieving estimation error approaching zero as the number of observations grows, when facing an oblivious adversary that may corrupt responses in all but an α fraction of the samples. As concrete examples, we investigate two problems: sparse regression and principal component analysis (PCA). For sparse regression, we achieve consistency for optimal sample size n & (k log d)/α2 and optimal error rate O( √ (k log d)/(n · α2)) where n is the number of observations, d is the number of dimensions and k is the sparsity of the parameter vector, allowing the fraction of inliers to be inverse-polynomial in the number of samples. Prior to this work, no estimator was known to be consistent when the fraction of inliers α is o(1/log log n), even for (non-spherical) Gaussian design matrices. Results holding under weak design assumptions and in the presence of such general noise have only been shown in dense setting (i.e., general linear regression) very recently by d’Orsi et al. (dNS21). In the context of PCA, we attain optimal error guarantees under broad spikiness assumptions on the parameter matrix (usually used in matrix completion). Previous works could obtain non-trivial guarantees only under the assumptions that the measurement noise corresponding to the inliers is polynomially small in n (e.g., Gaussian with variance 1/n2). To devise our estimators, we equip the Huber loss with non-smooth regularizers such as the `1 norm or the nuclear norm, and extend d’Orsi et al.’s approach (dNS21) in a novel way to analyze the loss function. Our machinery appears to be easily applicable to a wide range of estimation problems. We complement these algorithmic results with statistical lower bounds showing that the fraction of inliers that our PCA estimator can deal with is optimal up to a constant factor.

text: We develop machinery to design efficiently computable and consistent estimators, achieving estimation error approaching zero as the number of observations grows, when facing an oblivious adversary that may corrupt responses in all but an α fraction of the samples. As concrete examples, we investigate two problems: sparse regression and principal component analysis (PCA). For sparse regression, we achieve consistency for optimal sample size n & (k log d)/α2 and optimal error rate O( √ (k log d)/(n · α2)) where n is the number of observations, d is the number of dimensions and k is the sparsity of the parameter vector, allowing the fraction of inliers to be inverse-polynomial in the number of samples. Prior to this work, no estimator was known to be consistent when the fraction of inliers α is o(1/log log n), even for (non-spherical) Gaussian design matrices. Results holding under weak design assumptions and in the presence of such general noise have only been shown in dense setting (i.e., general linear regression) very recently by d’Orsi et al. (dNS21). In the context of PCA, we attain optimal error guarantees under broad spikiness assumptions on the parameter matrix (usually used in matrix completion). Previous works could obtain non-trivial guarantees only under the assumptions that the measurement noise corresponding to the inliers is polynomially small in n (e.g., Gaussian with variance 1/n2). To devise our estimators, we equip the Huber loss with non-smooth regularizers such as the `1 norm or the nuclear norm, and extend d’Orsi et al.’s approach (dNS21) in a novel way to analyze the loss function. Our machinery appears to be easily applicable to a wide range of estimation problems. We complement these algorithmic results with statistical lower bounds showing that the fraction of inliers that our PCA estimator can deal with is optimal up to a constant factor.

heading: 1 Introduction
text: Estimating information from structured data is a central theme in statistics that by now has found applications in a wide array of disciplines. On a high level, a typical assumption in an estimation problem is the existence of a –known a priori– family of probability distributions P : { β β ∈ Ω} over some spaceZ that are each indexed by some parameter β ∈ Ω. We then observe a collection
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
of n independent observations Z (Z1 , . . . , Zn) drawn from an unknown probability distribution β∗ ∈ P. The goal is to (approximately) recover the hidden parameter β∗.1
Oftentimes, real-world data may contain skewed, imprecise or corrupted measurements. Hence, a desirable property for an estimator is to be robust to significant, possibly malicious, noise perturbations on the given observations. Indeed, in the last two decades, a large body of work has been developed on designing robust algorithms (e.g see (BGN09; BBC11; MMYSB19)). However, proving strong guarantees often either demands strong assumptions on the noise model or requires that the fraction of perturbed observations is small. More concretely, when we allow the noise to be chosen adaptively, i.e., chosen dependently on the observations and hidden parameters, a common theme is that consistent estimators – estimators whose error tends to zero as the number of observations grows – can be attained only when the fraction of outliers is small.
In order to make vanishing error possible in the presence of large fractions of outliers, it is necessary to consider weaker adversary models that are oblivious to the underlying structured data. In recent years, a flurry of works have investigated oblivious noise models (CLMW11; ZLW+10; TJSO14; BJKK17; SZF18; SBRJ19; PJL21; dNS21). These results however are tailor-made to the specific models and problems. To overcome this limitation, in this paper we aim to provide a simple blueprint to design provably robust estimators under minimalistic noise assumptions for a large class of estimation problems. As a testbed for our blueprint, we investigate two well-studied problems:
Principal component analysis (PCA): Given a matrix Y : L∗ + N where L∗ ∈ n×n is an unknown parameter matrix and N is an n-by-n random noise matrix, the goal is to find an estimator L̂ for L∗ which is as close as possible to L∗ in Frobenius norm.
Sparse regression: Given observations (X1 , y1), . . . , (Xn , yn) following the linear model yi 〈Xi , β∗〉 + ηi where Xi ∈ d , β∗ ∈ d is the k-sparse parameter vector of interest (by k-sparse we mean that it has at most k nonzero entries) and η1 , . . . ηn is noise, the goal is to find an estimator β̂ for β∗ achieving small squared prediction error 1n ‖X(β̂ − β∗)‖2, where X is the matrix whose rows are X1 , . . . ,Xn .2
Principal component analysis. A natural way to describe principal component analysis under oblivious perturbations is that of assuming the noise matrix N to be an n-by-n matrix with a uniformly random set of α · n2 entries bounded by some small ζ > 0 in absolute value. In these settings, we may think of the α fraction of entries with small noise as the set of uncorrupted observations. Moreover, for ζ > 0, the fact that even for uncorrupted observations the noise is non-zero allows us to capture both gross sparse errors and small entry-wise noise in the measurements at the same time (for example if ζ 1 then the model captures settings with additional standard Gaussian noise). Remarkably, for ζ 0, Candès et al.’s seminal work (CLMW11) provided an algorithm that exactly recovers L∗ even for a vanishing fraction of inliers, under incoherence conditions on the signal L∗. The result was slightly extended in (ZLW+10) where the authors provided an algorithm recovering L∗ up to squared error O(ζ2 · n4) thus allowing polynomially small measurement noise ζ, but still failing to capture settings where standard Gaussian measurement noise is added to the sparse noise. Even for simple signal matrices L∗, prior to this work, it remained an open question whether consistent estimators could be designed in presence of both oblivious noise and more reasonable measurement noise (e.g., standard Gaussian).
Linear regression. Similarly to the context of principal component analysis, a convenient model for oblivious adversarial corruptions is that of assuming the noise vector η to have a random set of α · n entries bounded by ζ in absolute value (here all results of interest can be extended to any ζ > 0 just by scaling, so we consider only the case ζ 1). This model trivially captures the classical settings with Gaussian noise (a Gaussian vector with variance σ2 will have an α Θ( 1σ ) fraction of inliers) and, again, allows us to think of the α-fraction of entries with small noise as the set of uncorrupted observations. Early works on consistent regression focused on the regime with Gaussian design X1 , . . . ,Xn ∼ N(0,Σ) and deterministic noise η.3 (BJKK17) presented an estimator achieving error
1In this paper, we assume Ω ⊆ d andZ ⊆ D for some d ,D, denote random variables in bold face, and hide absolute constant factors with the notation O(·),Ω(·) , & , . and logarithmic factors with Õ(·), Ω̃(·).
2Our analysis also works for the parameter error ‖β∗ − β̂‖. 3For Gaussian design X one may also consider a deterministic noise model. This noise model is subsumed by the random noise model discussed here (if X is Gaussian). As shown in (dNS21), roughly speaking the
Õ(d/(α2 · n)) for any α larger than a fixed constant. (SBRJ19) extended this result by achieving comparable error rates even for a vanishing fraction of inliers α & 1/log log n. Assuming n & d2/α2 (TJSO14) proved that the Huber-loss estimator (Hub64) achieves optimal error rate O(d/α2 · n) even for polynomially small fraction of inliers α & √ d/n. This line of work culminated in (dNS21) which extended the result of (TJSO14) achieving the same guarantees with optimal sample complexity n & d/α2. For Gaussian design, similar result as (dNS21) can be extracted from independent work (PF20). Furthermore, the authors of (dNS21) extended these guarantees to deterministic design matrices satisfying only a spreadness condition (trivially satisfied by sub-Gaussian design matrices).
Huber loss was also analyzed in context of linear regression in (SZF18; DT19; SF20; PJL21). These works studied models that are different from our model, and these results do not work in the case α o(1). (SZF18) used assumptions on the moments of the noise, (DT19) and (SF20) studied the model with non-oblivious adversary, and the model from (PJL21) allows corruptions in covariates.
In the context of sparse regression, less is known. When the fraction of observations is constant α > Ω(1), (NTN11) provided a first consistent estimator. Later (DT19) and (SF20) improved that result (assuming α > Ω(1), but these results also work with non-oblivious adversary). In the case α o(1), the algorithm of (SBRJ19) (for Gaussian designs) achieves nearly optimal error convergence Õ ( (k log2 d)/α2n ) , but requires α & 1/log log n. More recently, (dNS21) presented an algorithm for standard Gaussian design X ∼ N(0, 1)n×d , achieving the nearly optimal error convergence Õ(k/(α2 · n)) for nearly optimal inliers fraction α & √ (k · log2 d)/n. Both (SBRJ19) and (dNS21) use an iterative process, and require a bound ‖β∗‖ 6 dO(1) (for larger ‖β∗‖ these algorithms also work, but the fraction of inliers or the error convergence is worse). The algorithm from (dNS21) , however, heavily relies on the assumption that X ∼ N(0, 1)n×d and appears unlikely to be generalizable to more general families of matrices (including non-spherical Gaussian designs).
Our contribution. We propose new machinery to design efficiently computable consistent estimators achieving optimal error rates and sample complexity against oblivious outliers. In particular, we extend the approach of (dNS21) to structured estimation problems, by finding a way to exploit adequately the structure therein. While consistent estimators have already been designed under more benign noise assumptions (e.g. the LASSO estimator for sparse linear regression under Gaussian noise), it was previously unclear how to exploit this structure in the setting of oblivious noise. One key consequence of our work is hence to demonstrate what minimal assumptions on the noise are sufficient to make effective recovery (in the sense above) possible. Concretely, we show
Oblivious PCA: Under mild assumptions on the noise matrix N and common assumption on the parameter matrix L∗ –traditionally applied in the context of matrix completion (NW12)– we provide an algorithm that achieves optimal error guarantees.
Sparse regression: Under mild assumptions on the design matrix and the noise vector –similar to the ones used in (dNS21) for dense parameter vectors β∗ – we provide an algorithm that achieves optimal error guarantees and sample complexity.
For both problems, our analysis improves over the state-of-the-art and recovers the classical optimal guarantees, not only for Gaussian noise, but also under much less restrictive noise assumptions. At a high-level, we achieve the above results by equipping the Huber loss estimator with appropriate regularizers. Our techniques closely follow standard analyses for M-estimators, but crucially depart from them when dealing with the observations with large perturbations. Furthermore, our analysis appears to be mechanical and thus easily applicable to many different estimation problems.

heading: 2 Results
text: Our estimators are based on regularized versions of the Huber loss. The regularizer we choose depends on the underlying structure of the estimation problem: We use `1 regularization to enforce sparsity in linear regression and nuclear norm regularization to enforce a low-rank structure in the context of PCA. More formally, the Huber penalty is defined as the function fh : → >0 such that underlying reason is that the Gaussianity of X allows one to obtain several other desirable properties "for free". For example, one could ensure that the noise vector is symmetric by randomly flipping the sign of each observation (yi ,Xi), as the design matrix will still be Gaussian. See Section 2 for a more in-depth discussion.
fh(t) : { 1 2 t 2 for |t | 6 h ,
h(|t | − h2 ) otherwise. (2.1)
where h > 0 is a penalty parameter. For X ∈ D , the Huber loss is defined as the function Fh(X) : ∑ i∈[D] fh(Xi). We will define the regularized versions in the following sections. For a matrix A, we use ‖A‖, ‖A‖nuc, ‖A‖F, ‖A‖max to denote its spectral, nuclear, Frobenius, maximum4 norms, respectively. For a vector v, we use ‖v‖ and ‖v‖1 to denote its `2 and `1 norms.

heading: 2.1 Oblivious principal component analysis
text: For oblivious PCA, we provide guarantees for the following estimator (ζ will be defined shortly):
L̂ B argmin L∈ n×n , ‖L‖max6ρ/n
( Fh(Y − L) + 100 √ n ( ζ + ρ/n ) ‖L‖nuc ) . (2.2)
Theorem 2.1. Let L∗ ∈ n×n be an unknown deterministic matrix and let N be an n-byn random matrix with independent, symmetrically distributed (about zero) entries and α : mini , j∈[n] { Ni j 6 ζ} for some ζ > 0. Suppose that rank(L∗) r and ‖L∗‖max 6 ρ/n. Then, with probability at least 1 − 2−n over N , given Y L∗ + N , ζ and ρ, the estimator Eq. (2.2) with Huber parameter h ζ + ρ/n satisfies L̂ − L∗ F 6 O (√rnα ) · (ζ + ρ/n) . We first compare the guarantees of Theorem 2.1 with the previous results on robust PCA (CLMW11; ZLW+10).5 The first difference is that they require L∗ to satisfy certain incoherent conditions.6 Concretely, they provide theoretical guarantees for r 6 O ( µ−1n(log n)−2 ) , where µ is the incoherence parameter. In certain regimes, such a constraint strongly binds with the eigenvectors of L∗, restricting the set of admissible signal matrices. Using the different assumption ‖L∗‖max 6 ρ/n (commonly used for matrix completion, see (NW12)), we can obtain nontrivial guarantees (i.e.
L̂ − L∗ F/‖L∗‖F → 0 as n →∞) even when the µ-incoherence conditions are not satisfied for any µ 6 n/log2 n, and hence the results (CLMW11; ZLW+10) cannot be applied. We remark that, without assuming incoherence, the dependence of the error on the maximal entry of L∗ is inherent (see Remark 3.2).
The second difference is that Theorem 2.1 provides a significantly better dependence on the magnitude ζ of the entry-wise measurement error. Specifically, in the settings of Theorem 2.1, (ZLW+10) showed that if L∗ satisfies the incoherence conditions, the error of their estimator is O ( n2ζ ) . If the entries of N are standard Gaussian with probability α (and hence ζ 6 O(1)), and the entries of L∗ are bounded by O(1), then the error of our estimator is O( √ rn/α), which is considerably better than O(n2) as in (ZLW+10). On the other hand, their error does not depend on the magnitude ρ/n of the signal entries, so in the extreme regimes when the singular vectors of L∗ satisfy the incoherence conditions but L∗ has very large singular values (so that the magnitude of the entries of L is significantly larger than n), their analysis provides better guarantees than Theorem 2.1.
As another observation to understand Theorem 2.1, notice that our robust PCA model also captures the classical matrix completion settings. In fact, any instance of matrix completion can be easily transformed into an instance of our PCA model: for the entries (i , j) that we do not observe, we can set Ni , j to some arbitrarily large value ±C(ρ, n) ρ/n, making the signal-to-noise ratio of the entry arbitrarily small. The observed (i.e. uncorrupted) Θ(α · n2) entries may additionally be perturbed by Gaussian noise with variance Θ(ζ2). The error guarantees of the estimator in Theorem 2.1 is
4For n × m matrix A, ‖A‖max maxi∈[n], j∈[m] |Ai j |. 5We remark that in (CLMW11) the authors showed that they can consider non-symmetric noise when the fraction of inliers is large α > 1/2. However for smaller fraction of inliers their analysis requires the entries of the noise to be symmetric and independent, so for α < 1/2, their assumptions are captured by Theorem 2.1.
6A rank-r n × n dimensional matrix M is µ-incoherent if its singular vector decomposition M : UΣV>
satisfies maxi∈[n]‖U>ei ‖2 6 µr n , maxi∈[n]‖V>ei ‖2 6 µr n and ‖UV>‖∞ 6 √ µr n .
O ( ( ρ/n + ζ )√ rn/α ) . Thus, the dependency on the parameters ρ, n , ζ, and r is the same as in matrix
completion and the error is within a factor of Θ( √
1/α) from the optimum for matrix completion. However, this worse dependency on α is intrinsic to the more general model considered and it turns out to be optimal (see Theorem 2.4). On a high level, the additional factor of Θ( √ 1/α) comes from the fact that in our PCA model we do not know which entries are corrupted. The main consequence of this phenomenon is that a condition of the form α & √ r/n appears inherent to achieve consistency. To get some intuition on why this condition is necessary, consider the Wigner model where we are given a matrix Y xx> + σW for a flat vector x ∈ {±1}n and a standard Gaussian matrix W . Note, that the entries of W fit our noise model for ζ 1 , ρ/n 1 , r 1 and α Θ(1/σ). The spectral norm of σW concentrates around 2σ √ n and thus it is information-theoretically impossible to approximately recover the vector x for σ 1/α ω( √ n) (see (PWBM16)).

heading: 2.2 Sparse regression
text: Our regression model considers a fixed design matrix X ∈ n×d and observations y : Xβ∗+η ∈ n where β∗ is an unknown k-sparse parameter vector and η is random noise with (|ηi | 6 1) > α for all i ∈ [n]. Earlier works (BJKK17), (SBRJ19) focused on the setting that the design matrix consists of i.i.d. rows with Gaussian distribution N(0,Σ) and the noise is η ζ + w where ζ is deterministic (α · n)-sparse vector and w is subgaussian. As in (dNS21), our results for a fixed design and random noise can, in fact, extend to yield the same guarantees for this early setting (see Theorem 2.3). Hence, a key advantage of our results is that the design X does not have to consist of Gaussian entries. Remarkably, we can handle arbitrary deterministic designs as long as they satisfy some mild conditions. Concretely, we make the following three assumptions, the first two of which are standard in the literature of sparse regression (e.g., see (Wai19), section 7.3):
1. For every column X i of X, X i 6 √n.
2. Restricted eigenvalue property (RE-property): For every vector u ∈ d such that7 usupp(β∗) 1 > 0.1 · ‖u‖1, we have 1n ‖Xu‖2 > λ · ‖u‖2 for some parameter λ > 0. 3. Well-spreadness property: For some (large enough) m ∈ [n] and for every vector u ∈ d
such that usupp(β∗) 1 > 0.1 · ‖u‖1 and for every subset S ⊆ [n] with |S | > n − m, it holds that ‖(Xu)S‖ > 12 ‖Xu‖.
Denote F2(β) : n∑
i 1
f2 ( yi − 〈Xi , β〉 ) , where Xi are the rows of X, and f2 is as in Eq. (2.1). We
devise our estimator for sparse regression and state its statistical guarantees below:
β̂ B arg min β∈ d
( F2(β) + 100 √ n log d · β 1) . (2.3) Theorem 2.2. Let β∗ ∈ d be an unknown k-sparse vector and let X ∈ n×d be a deterministic matrix such that for each column X i of X, ‖X i ‖ 6 √ n, satisfying the RE-property with λ > 0 and well-spreadness property with m & k log d λ·α2 (recall that n > m).
Further, let η be an n-dimensional random vector with independent, symmetrically distributed (about zero) entries and α mini∈[n] { ηi 6 1}. Then with probability at least 1 − d−10 over η, given X and y Xβ∗ + η, the estimator Eq. (2.3) satisfies
1 n X (β̂ − β∗) 2 6 O ( 1 λ · k log d α2 · n ) and
β̂ − β∗ 2 6 O ( 1 λ2 · k log d α2 · n ) .
There are important considerations when interpreting this theorem. The first is the special case η ∼ N(0, σ2)n , which satisfies our model for α Θ(1/σ). For this case, it is well known (e.g.,
7For a vector v ∈ d and a set S ⊆ [d], we denote by vS the restriction of v to the coordinates in S.
see (Wai19), section 7.3) that under the same RE assumption, the LASSO estimator achieves a prediction error rate of O( σ2λ · k log d n ) O( k log d λ·α2 ·n ), matching our result. Moreover, this error rate is essentially optimal. Under a standard assumption in complexity theory (NP 1 P/poly), the RE assumption is necessary when considering polynomial-time estimators (ZWJ14). Further, this also shows that the dependence on the RE constant seems unavoidable. Under mild conditions on the design matrix, trivially satisfied if the rows are i.i.d. Gaussian with covariance Σ whose condition number is constant, our guarantees are optimal up to constant factors for all estimators if k 6 d1−Ω(1) (e.g. k 6 d0.99), see (RWY11). This optimality also shows that our bound on the number of samples is best possible since otherwise we would not be able to achieve vanishing error. The (non-sparse version of) well-spreadness property was first used in the context of regression in (dNS21). In the same work the authors also showed that, under oblivious noise assumptions, some weak form of spreadness property is indeed necessary.
The second consideration is the optimal dependence on α: Theorem 2.2 achieves consistency as long as the fraction of inliers satisfies α ω( √ k log d/n). To get an intuition, observe that lower bounds for standard sparse regression show that, already for η ∼ N(0, σ · Idn), it is possible to achieve consistency only for n ω(σ2k log d) (if k 6 d1−Ω(1)). As for this η, the number of entries of magnitude at most 1 is O(n/σ) with high probability, it follows that for α Θ(1/σ) 6 O (√ (k log d)/n ) , no estimator is consistent.
To the best of our knowledge, Theorem 2.2 is the first result to achieve consistency under such minimalistic noise settings and deterministic designs. Previous results (BJKK17; SBRJ19; dNS21) focused on simpler settings of Gaussian design X and deterministic noise, and provide no guarantees for more general models. The techniques for Theorem 2.2 also extend to this case.
Theorem 2.3. Let β∗ ∈ d be an unknown k-sparse vector and let X be a n-by-d random matrix with i.i.d. rows X1 , . . .Xn ∼ N(0,Σ) for a positive definite matrix Σ. Further, let η ∈ n be a deterministic vector with α·n coordinates bounded by 1 in absolute value. Suppose that n & ν(Σ)·k log d
σmin(Σ)·α2 , where ν(Σ) is the maximum diagonal entry of Σ and σmin(Σ) is its smallest eigenvalue. Then, with probability at least 1 − d−10 over X , given X and y Xβ∗ + η, the estimator Eq. (2.3) satisfies
1 n X (β̂ − β∗) 2 6 O ( ν(Σ) · k log d σmin(Σ) · α2 · n ) and
β̂ − β∗ 2 6 O ( ν(Σ) · k log d σ2min(Σ) · α2 · n ) .
Even for standard Gaussian design X ∼ N(0, 1)n×d , the above theorem improves over previous results, which required sub-optimal sample complexity n & ( k/α2 ) · log d · log
β∗ . For non-spherical Gaussian designs, the improvement over state of the art (SBRJ19) is more serious: their algorithm requires α > Ω(1/log log n), while our Theorem 2.3 doesn’t have such restrictions and works for all α & √ ν(Σ) σmin(Σ) · k log d n ; in many interesting regimes α is allowed to be smaller than n −Ω(1). The
dependence on α is nearly optimal: the estimator is consistent as long as α > ω (√
ν(Σ) σmin(Σ) · k log d n
) ,
and from the discussion after Theorem 2.3, if α 6 O( √ (k log d)/n), no estimator is consistent.
Note that while we can deal with general covariance matrices, to compare Theorem 2.3 with Theorem 2.2 it is easier to consider Σ in a normalized form, when ν(Σ) 6 1. This can be easily achieved by scaling X . Also note that Theorem 2.2 can be generalized to the case ‖X i ‖ 6 √ νn for arbitrary ν > 0, and then the error bounds and the bound on m should be multiplied by ν.
The RE-property of Theorem 2.2 is a standard assumption in sparse regression and is satisfied by a large family of matrices. For example, with high probability a random matrix X with i.i.d. rows sampled from N(0,Σ), with positive definite Σ ∈ d×d whose diagonal entries are bounded by 1, satisfies the RE-property with parameter Ω(σmin(Σ)) for all subsets of [d] of size k (so for every possible support of β∗) as long as long as n & 1σmin(Σ) · k log d (see (Wai19), section 7.3.3). The wellspread assumption is satisfied for such X with for all sets S ⊂ [n] of size m 6 cn (for sufficiently small c) and for all subsets of [d] of size k as long as n & 1σmin(Σ) · k log d.

heading: 2.3 Optimal fraction of inliers for principal component analysis under oblivious noise
text: We show here that the dependence on α we obtain in Theorem 2.1 is information theroetically optimal up to constant factors. Concretely, let L∗ ,N , Y , α, ρ and ζ be as in Theorem 2.1, and let 0 < ε < 1 and 0 < δ < 1. A successful (ε, δ)-weak recovery algorithm for PCA is an algorithm that takes Y as input and returns a matrix L̂ such that
L̂ − L∗ F 6 ε · ρ with probability at least 1 − δ. It can be easily seen that the Huber-loss estimator of Theorem 2.1 fails to be a successful weakrecovery algorithm if α o( √ r/n) (for both cases ζ 6 ρ/n and ρ/n 6 ζ, we need α Ω( √ r/n).)
A natural question to ask is whether the condition α Ω( √
r/n) is necessary in general. The following theorem shows that if α o( √ r/n), then weak-recovery is information-theoretically impossible. This means that the (polynomially small) fraction of inliers that the Huber-loss estimator of Theorem 2.1 can deal with is optimal up to a constant factor. Theorem 2.4. There exists a universal constant C0 > 0 such that for every 0 < ε < 1 and 0 < δ < 1, if α : mini , j∈[n] [|Ni , j | 6 ζ] satisfies α < C0 · (1 − ε2)2 · (1 − δ) · √ r/n, and n is large enough, then it is information-theoretically impossible to have a successful (ε, δ)-weak recovery algorithm. The problem remains information-theoretically impossible (for the same regime of parameters) even if we assume that L∗ is incoherent; more precisely, even if we know that L∗ has incoherence parameters that are as good as those of a random flat matrix of rank r, the theorem still holds.

heading: 3 Techniques
text: To illustrate our techniques in proving statistical guarantees for the Huber-loss estimator, we first use sparse linear regression as a running example. Then, we discuss how the same ideas apply to principal component analysis. Finally, we also remark our techniques for lower bounds.

heading: 3.1 Sparse linear regression under oblivious noise
text: We consider the model of Theorem 2.2. Our starting point to attain the guarantees for our estimator Eq. (2.3), i.e., β̂ B arg minβ∈ d F2(β) + 100 √ n log d
β 1, is a classical approach for M-estimators (see e.g. (Wai19), chapter 9). For simplicity, we will refer to F2(β) as the loss function and to
β 1 as the regularizer. At a high level, it consists of the following two ingredients:
(I) an upper bound on some norm of the gradient of the loss function at the parameter β∗, (II) a lower bound on the curvature of the loss function (in form of a local strong convexity bound)
within a structured neighborhood of β∗. The structure of this neighborhood can roughly be controlled by choosing the appropriate regularizer.
The key aspect of this strategy is that the strength of the statistical guarantees of the estimator crucially depends on the directions and the radius in which we can establish lower bounds on the curvature of the function. Since these features are inherently dependent on the landscape of the loss function and the regularizer, they may differ significantly from problem to problem. This strategy has been applied successfully for many related problems such as compressed sensing or matrix completion albeit with standard noise assumptions.8 Under oblivious noise, (dNS21) used a particular instantiation of this framework to prove optimal convergence of the Huber-loss – without any regularizer – for standard linear regression. Such estimator, however, doesn’t impose any structure on the neighborhood of β∗ considered in (II) and thus, can only be used to obtain sub-optimal guarantees for sparse regression.
In the context of sparse regression, the above two conditions translate to: (I) an upper bound on the largest entry in absolute value of the gradient of the loss function at β∗, and (II) a lower bound on the curvature of F2 within the set of approximately k-sparse vectors9 close to β∗. We use this recipe to show that all approximate minimizers of F are close to β∗. While the idea of restricting to only approximately sparse directions has also been applied for the LASSO estimator in sparse
8The term "standard noise assumptions" is deliberately vague; as a concrete example, we will refer to (sub)-Gaussian noise distributions. See again chapter 9 of (Wai19) for a survey.
9We will clarify this notion in the subsequent paragraphs.
regression under standard (sub)-Gaussian noise, in the presence of oblivious noise, our analysis of the Huber-loss function requires a more careful approach.
More precisely, under the assumptions of Theorem 2.2, the error bound can be computed as O ( s‖G‖∗reg
κ
) , (3.1)
where G is a gradient of Huber loss at β∗, ‖·‖∗reg is a norm dual to the regularization norm (which is equal to ‖·‖max for `1 regularizer), s is a structure parameter, which is equal to √ k/λ, and κ is a restricted strong convexity parameter. Note that by error here we mean 1√ n ‖X(β̂ − β∗)‖.
Similarly, under under assumptions of Theorem 2.1, we get the error bound Eq. (3.1), where G is a gradient of Huber loss at L∗, ‖·‖∗reg is a norm dual to the nuclear norm (i.e. the spectral norm), structure parameter s is √ r, and κ is a restricted strong convexity parameter. For more details on the conditions of the error bound, see the supplementary material. Below, we explore the bounds on the norm of the gradient and on the restricted strong convexity parameter.
Bounding the gradient of the Huber loss. The gradient of the Huber-loss F2(·) at β∗ has the form ∇F2(β∗) ∑n i 1 f ′ 2 [ ηi ] · Xi , where Xi is the i-th row of X. The random variables f ′2 [ ηi ] , i ∈ [n], are independent, centered, symmetric and bounded by 2. Since we assume that each column of X has norm at most √ n, the entries of the row Xi are easily bounded by √ n. Thus, ∇F2(β∗) is a vector with independent, symmetric entries with bounded variance, so its behavior can be easily studied through standard concentration bounds. In particular, by a simple application of Hoeffding’s inequality, we obtain, with high probability, ∇F2(β∗) max maxj∈[d] ∑
i∈[n] f ′2
[ ηi ] · Xi j 6 O (√n log d) . (3.2) Local strong convexity of the Huber loss. Proving local strong convexity presents additional challenges. Without the sparsity constraint, (dNS21) showed that under a slightly stronger spreadness assumptions than Theorem 2.2, the Huber loss is locally strongly convex within a constant radius R ball centered at β∗ whenever n & d/α2. (This function is not globally strongly convex due to its linear parts.) Using their result as a black-box, one can obtain the error guarantees of Theorem 2.2, but with suboptimal sample complexity. The issue is that with the substantially smaller sample size of n > Õ(k/α2) that resembles the usual considerations in the context of sparse regression, the Huber loss is not locally strongly convex around β∗ uniformly across all directions, so we cannot hope to prove convergence with optimal sample complexity using this argument. To overcome this obstacle, we make use of the framework of M-estimators: Since we consider a regularized version of the Huber loss, it will be enough to show local strong convexity in a radius R uniformly across all directions which are approximately k-sparse. For this substantially weaker condition, Õ(k/α2) will be enough. More in details, for observations y Xβ∗ + η and an arbitrary u ∈ d of norm ‖u‖ 6 R, it is possible to lower bound the Hessian10 of the Huber loss at β∗ + u by:11
HF2(β∗ + u) n∑
i 1
f ′′2 [ (Xu)i − ηi ] · XiXiT n∑ i 1 1[|(Xu)i−ηi |62] · XiXi T
M(u) : n∑
i 1
1[|〈Xi ,u〉|61] · 1[|ηi |61] · XiXi T
As can be observed, we do not attempt to exploit cancellations between Xu and η. Let Q : { i ∈ [n] ηi 6 1} be the set of uncorrupted entries of η. Given that with high probability Q has size Ω(α ·n), the best outcome we can hope for is to provide a lower bound of the form 〈u ,M(u)u〉 > αn in the direction β̂ − β∗ . In (dNS21), it was shown that if the span of the measurement matrix X is well spread, then 〈u ,M(u)u〉 > Ω(α · n).
10The Hessian does not exist everywhere. Nevertheless, the second derivative of the penalty function f2 exists
as an L1 function in the sense that f ′2(b) − f ′ 2(a) 2 ∫ b a 1[|t |62]dt. This property is enough for our purposes.
11A more extensive explanation of the first part of this analysis can be found in (dNS21).
If the direction β̂− β∗ was fixed, it would suffice to show the curvature in that single direction through the above reasoning. However, β̂ depends on the unknown random noise vector η . Without the regularizer in Eq. (2.3), this dependence indicates that the vector β̂ may take any possible direction, so one needs to ensure local strong convexity to hold in a constant-radius ball centered at β∗. That is, min‖u‖6R λmin(M(u)) > Ω(α · n) . It can be shown through a covering argument (of the ball) that this bound holds true for n > d/α2. This is the approach of (dNS21).
The minimizer of the Huber loss follows a sparse direction. The main issue with the above approach is that no information concerning the direction β̂ − β∗ is used. In the settings of sparse regression, however, our estimator contains the regularizer ‖β‖1. The main consequence of the regularizer is that the direction β̂−β∗ is approximately flat in the sense ‖β̂−β∗‖1 6 O (√ k‖β̂ − β∗‖ ) .
The reason12 is that due to the structure of the objective function in Eq. (2.3) and concentration of the gradient Eq. (3.2), the penalty for dense vectors is larger than the absolute value of the inner product 〈∇F(β∗), β̂ − β∗〉 (which, as previously argued, concentrates around its (zero) expectation). This specific structure of the minimizer implies that it suffices to prove local strong convexity only in approximately sparse directions. For these set of directions, we carefully construct a sufficiently small covering set so that n > Õ(k/α2) samples suffice to ensure local strong convexity over it. Remark 3.1 (Comparison with LASSO). it is important to remark that while this approach of only considering approximately sparse directions has also been used in the context of sparse regression under Gaussian noise (e.g. the LASSO estimator), obtaining the desired lower bound is considerably easier in these settings as it directly follows from the restricted eigenvalue property of the design matrix. In our case, we require an additional careful probabilistic analysis which uses a covering argument for the set of approximately sparse vectors. As we see however, it turns out that we do not need any additional assumptions on the design matrix when compared with the LASSO estimator except for the well-spreadness property (recall that some weak version of well-spreadness is indeed necessary in robust settings, see (dNS21)).

heading: 3.2 Principal component analysis under oblivious noise
text: A convenient feature of the approach in Section 3.1 for sparse regression, is that it can be easily applied to additional problems. We briefly explain here how to apply it for principal component analysis. We consider the model defined in Theorem 2.1. We use an estimator based on the Huber loss equipped with the nuclear norm as a regularizer to enforce the low-rank structure in our estimator
L̂ B argmin L∈ n×n , ‖L‖max6ρ/n
( Fζ+ρ/n(Y − L) + 100 √ n ( ζ + ρ/n ) ‖L‖nuc ) . (3.3)
In this setting, the gradient ∇Fζ+ρ/n(Y − L∗) is a matrix with independent, symmetric entries which are bounded (by ζ + ρ/n) and hence its spectral norm is O ( (ζ + ρ/n) √ n )
with high probability. Local strong convexity can be obtained in a similar fashion as shown in Section 3.1: due to the choice of the Huber transition point all entries with small noise are in the quadratic part of F. Moreover, the nuclear norm regularizer ensures that the minimizer is an approximately low-rank matrix in the sense that ‖M‖nuc 6 O (√ r‖M‖F ) . So again, it suffices to provide curvature of the loss function only on these subset of structured directions. Remark 3.2 (Incoherence vs. spikiness). Recall the discussion on incoherence in Section 2. If for every µ 6 n/log2 n matrix L∗ doesn’t satisfy µ-incoherence conditions, the results in (CLMW11; ZLW+10) cannot be applied. However, our estimator achieves error ‖L̂−L∗‖F/‖L∗‖F → 0 as n →∞. Indeed,let ω(1) 6 f (n) 6 o(log2 n) and assume ζ 0. Let u ∈ n be an f (n)-sparse unit vector whose nonzero entries are equal to 1/ √ f (n). Let v ∈ n be a vector with all entries equal to 1/ √ n. Then, uvT does not satisfy incoherence with any µ < n/ f (n). We have uvT F 1, and the error of
our estimator is O(1/(α √
f (n))), so it tends to zero for constant (or even some subconstant) α. Furthermore notice that the dependence of the error of Theorem 2.1 on the maximal entry of L∗ is inherent if we do not require incoherence. Indeed, consider L1 b · e1e1T for large enough b > 0
12This phenomenon is a consequence of the decomposability of the L1 norm, see the supplementary material.
and L2 e2e2T. For constant α, let |Ni j | be 1 with probability α/2, 0 with probability α/2 and b with probability 1 − α. Then given Y we cannot even distinguish between cases L∗ L1 or L∗ L2, and since ‖L1 − L2‖F > b, the error should also depend on b. Remark 3.3 (α vs α2: what if one knows which entries are corrupted?). As was observed in Section 2, the error bound of our estimator is worse than the error for matrix completion by a factor 1/ √ α. We observe similar effect in linear regression: if, as in matrix completion, we are given a randomly chosen α fraction of observations {( Xi , yi 〈Xi , β∗〉 + ηi )}n i 1 where η ∼ N(0, 1)
n , and since for the remaining samples we may not assume any bound on the signal-to-noise ratio, then this problem is essentialy the same as linear regression with αn observations. Thus the optimal prediction error rate is Θ( √ d/(αn)). Now, if we have y Xβ∗ + η, where η ∼ N(0, 1/α2), then with probability
Θ(α), |ηi | 6 1, but the optimal prediction error rate in this case is Θ( √
d/(α2n)). So in both linear regression and robust PCA, prior knowledge of the set of corrupted entries makes the problem easier.

heading: 3.3 Optimal fraction of inliers for principal component analysis under oblivious noise
text: In order to prove Theorem 2.4, we will adopt a generative model for the hidden matrix L∗: We will generate L∗ randomly but assume that the distribution is known to the algorithm. This makes the problem easier. Therefore, any impossibility result for this generative model would imply impossibility for the more restrictive model for which L∗ is deterministic but unknown.
We generate a random flat matrix L∗ using n · r independent and uniform random bits in such a way that L∗ is of rank r and incoherent with high probability. Then, for every constant 0 < ξ < 1, we find a distribution for the random noise N in such a way that the fraction of inliers satisfies α : [|Ni j | 6 ζ] Θ ( ξ √ r/n ) , and such that the mutual information between L∗ and Y L∗ + N can be upper bounded as I(L∗; Y ) 6 O(ξ · n · r). Roughly speaking the smaller ξ gets, the more independent L∗ and Y will be. Now using an inequality that is similar to the standard Fano-inequality but adapted to weak-recovery, we show that if there is a successful (ε, δ)-weak recovery algorithm for L∗ and N , then I(L∗; Y ) > Ω ( (1 − ε2)2 · (1 − δ) · n · r ) . By combining all these observations together, we can deduce that if ξ is small enough, it is impossible to have a successful (ε, δ)-weak recovery algorithm for L∗ and N .
Acknowledgments and Disclosure of Funding
The authors thank the anonymous reviewers for useful comments. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 815464).

Review summaries:
This paper studies sparse regression and principal component analysis under  oblivious perturbations.   Main contributions: 1. The authors propose estimators that achieve the optimal error guarantees for PCA and sparse regression under oblivious perturbations separately, by minimizing the Huber loss function with some regularization. For sparse regression the estimator also achieves the optimal sample complexity.  2. This machinery has the potential to be applied to other estimation problems.  Post rebuttal: I acknowledge that I have read the authors responses and other reviewers' comments.

The paper studies algorithms for sparse regression and principle component analysis using Huber loss in the presence of overwhelming oblivious noise. In particular, the fraction of inliers, $\alpha$, is going to $0$ but it is still possible to achieve consistency as the noise model is benign (oblivious, symmetric noise). The paper obtains optimal error rates and dependence on $\alpha$ for a wide range of distribution assumptions. These rates are achieved by using regularized Huber loss (nuclear norm for the matrix case and $\ell_1$-norm for regression case).  

The authors study PCA and sparse regression in the presence of an oblivious adversary. The observation model is as follows: the observations may have symmetric, independent, additive noise. Each noise term is bounded by $\zeta$ with probability at least $\alpha$ which can be much smaller than 1. The authors show that the estimators which minimize Huber loss with appropriate $\ell_1$ penalty achieve optimal error rates.