Robust Regression Revisited: Acceleration and Improved Estimation Rates
Year: 2022

Abstract:
We study fast algorithms for statistical regression problems under the strong contamination model, where the goal is to approximately optimize a generalized linear model (GLM) given adversarially corrupted samples. Prior works in this line of research were based on the robust gradient descent framework of [PSBR20], a firstorder method using biased gradient queries, or the Sever framework of [DKK19], an iterative outlier-removal method calling a stationary point finder. We present nearly-linear time algorithms for robust regression problems with improved runtime or estimation guarantees compared to the state-of-the-art. For the general case of smooth GLMs (e.g. logistic regression), we show that the robust gradient descent framework of [PSBR20] can be accelerated, and show our algorithm extends to optimizing the Moreau envelopes of Lipschitz GLMs (e.g. support vector machines), answering several open questions in the literature. For the well-studied case of robust linear regression, we present an alternative approach obtaining improved estimation rates over prior nearly-linear time algorithms. Interestingly, our algorithm starts with an identifiability proof introduced in the context of the sum-of-squares algorithm of [BP21], which achieved optimal error rates while requiring large polynomial runtime and sample complexity. We reinterpret their proof within the Sever framework and obtain a dramatically faster and more sample-efficient algorithm under fewer distributional assumptions.

heading: 1 Introduction
text: Parameter estimation in generalized linear models (GLMs), such as linear and logistic regression problems, is among the most fundamental and well-studied statistical optimization problems. It serves as the primary workhorse in statistical studies arising from a variety of disciplines, ranging from economics [Smi12], biology [VGSM05], and the social sciences [Gor10]. Formally, given a link function γ : R2 → R and a dataset of covariates and labels {(Xi, yi)}i∈[n] ⊂ Rd × R drawn from an underlying distribution DXy , the problem of statistical (generalized linear) regression asks to
estimate θ? := argminθ∈Rd {
E (X,y)∼DXy
[γ(〈θ,X〉 , y)] } . (1)
For example, when γ(v, y) = 12 (v−y) 2, (1) corresponds to (statistical) linear regression. The problem (1) also has an interpretation as computing a maximum likelihood estimate for a parameterized distributional model for data generation, and indeed is only tractable under certain distributional
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
assumptions, since we only have access to samples from DXy rather than the underlying distribution itself (see e.g. [BP21] for tractability results in the linear regression setting).
However, in many modern settings these strong distributional assumptions fail to hold. In practically relevant settings, regression is often performed on massive datasets, where the data come from a poorly-understood distribution and have not been thoroughly vetted or cleaned of outliers. This has prompted the study of robust regression, i.e. regression under weak distributional or corruption assumptions. In this work, we study (1) in the strong contamination model. In this model, we assume the data points we receive are independently drawn from DXy, but that an arbitrary -fraction of the samples are then adversarially contaminated or replaced. The strong contamination model has garnered significant recent interest in the algorithmic statistics and learning communities for several reasons. Firstly, it is a flexible model of corruption and can be used to study both truly adversarial data poisoning attacks (where e.g. part of the dataset is sourced from malicious respondents), as well as model misspecification, where the generative DXy does not exactly satisfy our distributional assumptions, but is close in total variation to a distribution that does. Furthermore, a line of work building upon [DKK+16, LRV16] (discussed in our survey of prior work in Section 1.2) has achieved remarkable positive results for mean estimation and related problems under strong contamination, with statistical guarantees scaling independently of the dimension d. This dimension-free error promise is important in modern high-dimensional settings.

heading: 1.1 Our results
text: We give multiple nearly-linear time algorithms1 for problem (1) under the strong contamination model, with improved statistical or runtime guarantees compared to the state-of-the-art. Prior algorithms for (1) under the strong contamination model in the literature typically followed one of two frameworks. The first, which we refer to as robust gradient descent, was pioneered by [PSBR20], and is based on reframing (1) as a problem where we have noisy gradient access to an unknown function we wish to optimize, coupled with the design of a noisy gradient oracle based on a robust mean estimation primitive. The second, which we refer to as Sever, originated in work of [DKK+19], and uses the guarantees of stationary point finders such as stochastic gradient descent to repeatedly perform outlier removal. Interestingly, we show that both approaches can be dramatically sped up, and give two complementary types of algorithms within these frameworks.
Robust acceleration. We demonstrate that under the noisy gradient estimation framework for solving well-conditioned regression problems (1), an accelerated rate of optimization can be achieved, answering an open question of [PSBR20]. We give the following result for smooth statistical regression, where we assume the uncorrupted data is drawn from DXy with marginals DX and Dy; throughout, Õ hides polylogarithmic factors in problem parameters and failure probabilities. Finally, we remove κ dependences in sample complexity statements, as they are subsumed by dependences due to assumed bounds on κ or κ2. Theorem 1 (informal, cf. Theorem 7, supplement). Suppose γ : R2 → R is such that γy(v) := γ(v, y) is convex and has (absolute) first and second derivatives at most 1 for all y in the support of Dy , and DX has second moment matrix Σ? LI. For some µ ≥ 0, let κ := max(1, Lµ ) and let
θ?reg := argminθ∈Rd {
E (X,y)∼DXy {γ (〈θ,X〉 , y)}+ µ 2 ‖θ‖22 } be the solution to the true regularized statistical regression problem.2 There is an algorithm that given n := Õ(d ) -corrupted samples from DXy, for κ 2 at most an absolute constant, obtains θ
with ∥∥θ − θ?reg∥∥2 = O (√κ µ ) with high probability in time Õ(nd√κ).
A canonical example of a link function γ satisfying assumptions of Theorem 1 is the logit function γ(v, y) = log(1 + exp(−vy)), when the labels y are ±1. To contextualize Theorem 1, [PSBR20]
1Throughout, we reserve the description “nearly-linear” for runtimes scaling linearly in the dataset size nd, and polynomially in −1 and the condition number, up to a polylogarithmic overhead in problem parameters.
2To simplify bounds and avoid estimation error for non-strongly convex statistical regression problems scaling with the initial search radius (which may be dimension-dependent), we focus on regularized problems. There is a substantial line of work on reductions between rates for strongly convex and convex smooth optimization in the non-robust setting, see e.g. [ZH16], and we defer an analogous exploration in the robust setting to future work.
obtains a similar statistical guarantee in its setting, using Õ(κ) calls to a noisy gradient oracle, implemented via a subroutine based on robust mean estimation. Since then, [CAT+20] showed that for the case of linear regression (see Theorem 3 for the formal setup, as the linear regression link function is not Lipschitz), the framework was amenable to nearly-linear time mean estimation techniques of [CDG19], and gave an algorithm with runtime Õ(ndκ −6). Theorem 1 represents an improvement to these results on two fronts: we apply tools from [DHL19] to remove the poly( −1) runtime dependence for a general class of regression problems, and we achieve an iteration count of Õ( √ κ), matching the accelerated runtime of [Nes83] for (non-robust) smooth optimization. We remark that the application of tools inspired by [DHL19] is fairly straightforward, and not a primary contribution of our work compared to the accelerated dependence on κ.
We demonstrate the generality of our acceleration framework by applying it to optimizing the Moreau envelope for Lipschitz, but possibly non-smooth, link functions γ; a canonical example is the hinge loss γ(v, y) = max(0, 1 − vy) with ±1 labels, used in training support vector machines. The Moreau envelope is an extremely well-studied smooth approximation, and is an everywhere additive approximation if the original function is Lipschitz (see e.g. [Sho97]). In the non-robust setting many state-of-the-art rates for Lipschitz optimization are attainable by accelerated optimization of an appropriate Moreau envelope [TJNO20]. We show that even without explicit access to the Moreau envelope, we can approximately minimize it with our robust acceleration framework.
Theorem 2 (informal, cf. Theorem 8, supplement). Suppose γ : R2 → R is such that γy(v) := γ(v, y) is convex and has (absolute) first derivative at most 1 for all y in the support of Dy , and DX has bounded second moment matrix. For some µ, λ ≥ 0, let κ = max(1, 1λµ ) and let
θ?env := argminθ∈Rd { F ?λ (θ) + µ
2 ‖θ‖22
} , where F ?λ (θ) := inf
θ′
{ F ?(θ′) + 1
2λ ‖θ − θ′‖22 } is the Moreau envelope of F ?(θ) := E
(X,y)∼DXy {γ (〈θ,X〉 , y)} .
There is an algorithm that given n := Õ(d ) -corrupted samples from DXy, for κ 2 at most an absolute constant, obtains θ with ∥∥θ − θ?reg∥∥2 = O (√κ µ ) with high probability in time Õ(nd√κ ).
To obtain this result, we give a nearly-linear time construction of a noisy gradient oracle for the Moreau envelope, which may be of independent interest; we note similar gradient oracle constructions in different settings have been developed in the optimization literature (see e.g. [CJJS21]).
Robust linear regression. Perhaps the most ubiquitous example of statistical regression, the specific problem of robust linear regression has received substantial attention in the literature (cf. Section 1.2). Amongst the algorithms developed for the variant of this problem that we study, the only nearlylinear time algorithm is the recent work of [CAT+20].3 For a robust linear regression problem with noise variance bounded by σ2 and covariate second moment matrix Σ? := EX∼DX [XX>], [DKK+19, PSBR20, CAT+20] attain distance to the true regression minimizer θ? scaling as σκ √ in the Σ? norm (the “Mahalanobis distance”)4 under a bounded 4th moment assumption. We give one result (Theorem 3) which improves the runtime of [DKK+19, PSBR20, CAT+20] under the noisy gradient descent framework, and one result (Theorem 4) which improves its estimation rate, under the Sever framework.
We first demonstrate that applying our robust acceleration framework leads to a similar estimation guarantee as [DKK+19, PSBR20, CAT+20] under the same assumptions, but with improved runtime.
Theorem 3 (informal, cf. Theorem 6, supplement). Suppose DX is 2-to-4 hypercontractive with second moment matrix Σ? = EX∼DX [XX>] satisfying µI Σ
? LI, and y ∼ Dy is generated as 〈θ?, X〉 + δ, for δ ∼ Dδ with variance at most σ2 independent of X . Let κ := Lµ . There is an algorithm that given n := Õ(d ) -corrupted samples fromDXy , for κ
2 at most an absolute constant, obtains θ with ‖θ − θ?‖Σ? = O(σκ √ ) with high probability in time Õ(nd √ κ).
3Another algorithm was recently given by [] under different assumptions, most notably a known second moment matrix.
4We measure error in the Σ? norm as it is scale invariant and the natural norm to measure the underlying (quadratic) statistical regression problem under; some prior works gave `2 norm guarantees, which we convert.
We give a formal definition of 2-to-4 hypercontractivity in Section 2. We remark that attaining estimation rates for robust linear regression scaling polynomially in is known to be impossible under only bounded second moments ([BP21]); the 4th moment bound we require is the minimal assumption known in the literature under which such robust estimation is possible.5 Theorem 3 matches the distribution assumptions and error of [CAT+20], while obtaining an accelerated runtime.
Interestingly, under the 4th moment bound used in Theorem 3, [BP21] showed that the informationtheoretically optimal rate of estimation in the Σ? norm is independent of κ, and presented a matching upper bound under an analogous, but more stringent, distributional assumption.6 However, thus far robust linear regression algorithms have broadly fallen under two categories. The first family (e.g. [KKM18, ZJS20, BP21]), based on the sum-of-squares paradigm for algorithm design, sacrifices practicality to obtain improved error rates by paying a large runtime and sample complexity overhead (as well as requiring stronger distributional assumptions). The second (e.g. [DKK+19, PSBR20, CAT+20]), which opts for more practical approaches to algorithm design, has been bottlenecked at Mahalanobis distance O(σκ √ ) and the requirement that κ2 = O(1).
We present a nearly-linear time method for robust linear regression overcoming this bottleneck for the first time amongst non-sum-of-squares algorithms, attaining improved statistical rates compared to Theorem 3 while only requiring κ = O(1).
Theorem 4 (informal, cf. Theorem 5, supplement). Suppose DX is 2-to-4 hypercontractive with second moment matrix Σ? = EX∼DX [XX>] satisfying µI Σ
? LI, and y ∼ Dy is generated as 〈θ?, X〉+ δ, for δ ∼ Dδ , a 2-to-4 hypercontractive distribution with variance at most σ2 independent of X . Let κ := Lµ . There is an algorithm that given n := Õ(( d2 + d 4 )) samples from DXy , for κ at most an absolute constant, uses Õ( 1 ) calls to an empirical risk minimization algorithm 7 and Õ(nd ) additional runtime, and obtains θ with ‖θ − θ?‖Σ? = O(σ √ κ ) with probability at least 910 .
This second algorithm does require more resources than that of Theorem 3: the sample complexity scales quadratically in d, and the runtime is never faster. Further, we make the slightly stronger assumption of hypercontractive noise for the uncorrupted samples. On the other hand, the improved dependence on the condition number in the error can be significant for distributions in practice, which may be far from isotropic. All told, Theorem 4 presents an intermediate tradeoff inheriting some statistical gains of the sum-of-squares approach (albeit still depending on κ) without sacrificing a nearly-linear runtime. Interestingly, we obtain Theorem 4 by reinterpreting an identifiability proof used in the algorithm of [BP21], and combining it with tools inspired by the Sever robust estimation framework. We note that our sample complexity dramatically improves that of [DKK+17]’s original linear regression algorithm in the Sever framework in terms of the dependence on d (but not ), which requires Õ(d 5
2 ) samples (in addition to beating their weaker error guarantee). We elaborate on these points further in Section 4; we believe it is an interesting open problem to understand if the worse sample complexity of Theorem 4 is truly necessary for the improved statistical guarantees.

heading: 1.2 Prior work
text: We give a general overview contextualizing our work in this section, and defer the comparison of technical components we develop in this work to the relevant sections. Specifically, in Section 3 we discuss prior work on acceleration with noisy gradients, and how it compares to our results.
The study of learning in the presence of adversarial noise is known as robust statistics, with a long history dating back over 60 years [Ans60, Tuk60, Hub64, Tuk75, Hub04]. Despite this, the first efficient algorithms with near-optimal error for many fundamental high dimensional robust statistics problems were only recently developed [DKK+16, LRV16, DKK+17]. Since these works, efficient robust estimators have been developed for a variety of more complex problems; a comprehensive overview may be found in [DK19, Li18, Ste18].
5This is true not only for efficient algorithms but even information-theoretically in the case of DX . 6The algorithm of [BP21] requires DX to be certifiably hypercontractive, an algebraic condition frequently required by the sum-of-squares algorithmic paradigm to apply to robust statistical estimation problems. 7The empirical risk minimization algorithm used is up to the practitioner; its runtime will never scale worse than Õ(nd √ κ) by applying (non-robust) accelerated gradient descent, but can be substantially better if recent advances in stochastic gradient methods are used, e.g. [Zhu17].
Our results sit within the line of work in this field on robust stochastic optimization. The first works which achieved dimension-independent error rates with efficient algorithms for the problems we consider in this paper are the aforementioned works of [PSBR20, DKK+19]. Similar problems were previously considered in [CSV17a, BDLS17]. In [CSV17a], the authors consider a setting where a majority of the data is corrupted, and the goal is to output a short list of hypotheses so that at least one is close to the true regressor. However, because most of their data is corrupted, they achieve weaker statistical rates; in particular, their techniques do not achieve vanishing error as the fraction of error goes to zero. In [BDLS17], the authors consider a somewhat different model with stronger assumptions on the structure of the functions. In particular, they assume that the uncorrupted covariates are Gaussian with identity covariance, and are primarily concerned with the case where the regressors are sparse. Their main goal is to achieve sublinear sample complexities by leveraging sparsity. We also remark that the algorithms in [CSV17a, BDLS17] are also much more cumbersome, requiring heavy-duty machinery such as black-box SDP solvers and cutting plane methods, and as a result are more computationally intensive than those considered in [PSBR20, DKK+17].
There has been a large body of subsequent work on the special case of robust linear regression [KKM18, KKK19, DKS19, ZJS20, CAT+20, BP21]; however, the majority of this line of work focuses on achieving improved error rates under additional distributional assumptions by using the sum-of-squares hierachy. As a result, their algorithms are likely impractical in high dimensions, and require large (albeit polynomial) sample complexity and runtime. Of particular interest to us is [CAT+20], who combine the framework of [PSBR20] with the robust mean estimation algorithm of [CDG19] to achieve nearly-linear runtimes in the problem dimension and the number of samples. Our Theorem 3 can be thought of as the natural accelerated version of [CAT+20], with an additional −6 runtime overhead removed using more sophisticated mean estimation techniques.
After the initial submission of this paper, we were made aware of another recent work [PJL20] preventing a different approach to robust regression based on covariate filtering. Their work focuses on the highly well-conditioned setting, where the true second moment matrix has constant condition number, whereas our work primarily aims to improve various rates in terms of dependence on the condition number κ in the ill-conditioned regime (i.e. κ 1).

heading: 2 Preliminaries
text: Notation. For d ∈ N we let [d] := {j | j ∈ N, 1 ≤ j ≤ d}. The `p norm of a vector is ‖·‖p. The all-ones vector (of appropriate dimension from context) is 1, and the identity matrix is I. The (solid) probability simplex is ∆n := {w ∈ Rn≥0, ‖w‖1 ≤ 1}. For v ∈ Rn and S ⊆ [n], we let vS ∈ Rn zero out coordinates [n] \ S. We call S1, S2 a bipartition of S if S1 ∩ S2 = ∅ and S1 ∪ S2 = S. For symmetric d× dA,B we write A B to mean B−A is positive semidefinite. For positive definite M, we define the induced norm ‖v‖M := √ v>Mv. We use ‖·‖op to mean the `2-`2 operator norm.
Functions. We say differentiable f : Rd → R is λ-Lipschitz in ‖·‖M if ‖∇f(θ)‖M−1 ≤ λ for all θ ∈ Rd, and that twice-differentiable f : Rd → R is L-smooth and µ-strongly convex in ‖·‖M if µM ∇2f(θ) LM, for all θ ∈ Rd. When M is not specified, we assume M = I.
Distributions. For w ∈ ∆n and a set of vectors X := {Xi}i∈[n], the empirical second moment matrix is denoted Covw (X) := ∑ i∈[n]
wi ‖w‖1 XiX > i . We say distributionD supported on Rd is 2-to-4
hypercontractive if for all v ∈ Rd, EX∼D[〈X, v〉4] ≤ O(1)EX∼D[〈X, v〉2]2. Corruption model. We provide provable guarantees for statistical problems captured by the following standard statistical model (the “strong contamination model”): given link function γ : R2 → R and a “true” underlying distributionDXy on Rd×R with marginalsDX ,Dy , a dataset {(X̃i, ỹi)}i∈[n] is independently sampled from DXy, and an arbitrary (unknown) subset B ⊂ [n] is replaced by arbitrary points in supp(DXy). For i ∈ B, we observe the corrupted (Xi, yi), and otherwise we observe (Xi, yi) ← (X̃i, ỹi). When |B| = n, we say the dataset is -corrupted. We define fi(θ) := γ(〈Xi, θ〉 , yi), gi(θ) := ∇fi(θ), and for w ∈ ∆n, Fw(θ) := ∑ i∈[n] wifi(θ).
Filtering. Our algorithms make frequent use of an algorithmic primitive we call “filtering” [ABL14, DKK+17, SCV18]. Here, we have an index set [n] with fixed unknown bipartition G∪B; intuitively, G is a “good” set of indices we wish to keep, and B is a “bad” set we would like to remove. Naively, one would expect that we would chooseG to be the set of remaining uncorrupted points. For technical reasons, we will often have to choose (large) subsets of the set of remaining uncorrupted points.
The filtering algorithm maintains weightsw ∈ ∆n it iteratively downweights using “scores” τ ∈ Rn≥0, with the goal of producing weights close to the uniform distribution onG. We callw ∈ ∆n c-saturated (or saturated for short if c = 0) if w ≤ 1n1 entrywise and ∥∥[ 1n1− w]G∥∥1 ≤ ∥∥[ 1n1− w]B∥∥1 + c. The following lemma (implicit in [DKK+17, CSV17b, Li18, Ste18]) is representative of this technique. Lemma 1. Suppose G ∪B is a bipartition of [n], and 〈wG, τ〉 ≤ 〈wB , τ〉 for c-saturated w ∈ ∆n, τ ∈ Rn≥0. Then if w′i ← (1− τi τmax )wi ∀i ∈ [n] and τmax := maxi∈[n]|wi 6=0 τi, w′ is also c-saturated.
In other words, we can keep weights saturated (retaining most of the weight on G) by repeatedly identifying scores whose empirical average is large due to B. We use a variety of different scores; one subroutine which follows as a consequence of this paradigm, which is implicit in [DHL19], is a way to rapidly decrease the operator norm of an empirical second moment matrix whose restriction to a “good” majority set is bounded. Many of our algorithms use this subroutine, stated here. Lemma 2 (cf. Proposition 5, supplement). There is an algorithm, FastCovFilter, taking inputs V := {vi}i∈[n] ∈ Rn×d, saturated w ∈ ∆n with respect to [n] = G ∪B with |B| = n, and R ≥ 0 with the promise that ‖ 1|G| ∑ i∈G viv > i ‖op ≤ R. Then, FastCovFilter returns saturated w′ ∈ ∆n
such that ‖ ∑ i∈[n] w ′ iviv > i ‖op = O(R) with high probability in time Õ(nd).
Organization. In Section 3, we describe a general framework for attaining accelerated optimization rates under a “noisy gradient oracle” model, which we use to obtain Theorems 1, 2, and 3 by efficiently constructing appropriate oracles. In Section 4, we overview techniques used in designing the nearly-linear time linear regression algorithm of Theorem 4. Due to space constraints, we defer proofs, a precise statement of our statistical models and assumptions, and an extended exposition of our techniques to the unabridged version of this paper in the supplementary material.

heading: 3 Acceleration under noisy gradient oracle access
text: Our robust acceleration framework, which we use to prove Theorems 1, 2, and 3, addresses the following problem formulation: there is an unknown function F ? with minimizer θ? which is Lsmooth and µ-strongly convex, and we wish to estimate θ?, but our only mode of accessing F ? is through a noisy gradient oracle Ong. Namely, for some σ, , we can query Ong with some θ ∈ Rd and an upper bound R ≥ ‖θ − θ?‖2 and receive an estimate G(θ) such that
‖G(θ)−∇F ?(θ)‖2 = O (√ L σ + L √ R ) . (2)
In other words, we receive gradients perturbed by both fixed additive noise, and multiplicative noise depending on (a known upper bound on) distance to θ?. Prior works observed that by using tools from robust mean estimation, appropriate noisy gradient oracles could be constructed for the “mean” functions E(X,y)∼DXy [γ(〈X, θ〉 , y)] arising from the distributional assumptions in Theorems 1 and 3. EfficientOng construction. Our first contribution is speeding up the implementation ofOng to run in nearly-linear time Õ(nd), leveraging recent advances by [DHL19] for robust mean estimation. This
improves a similar construction by [CAT+20], running in time Õ(nd · −6) in the linear regression setting. To obtain this construction, we use the following lemma implicit in [DHL19].
Lemma 3 (Lemma 11, supplement). Letw ∈ ∆n be saturated with respect to bipartition [n] = G∪B, and let w?G := 1 |G|1G, w̃ := w ‖w‖1
. Then,∥∥∇Fw̃(θ)−∇Fw?G(θ)∥∥2 = O (√ ) (∥∥Covw?G{gi(θ)}i∈[n]∥∥ 12op + ∥∥Covw̃{gi(θ)}i∈[n]∥∥ 12op) . By demonstrating that our “uncorrupted” dataset has function gradients which are appropriately bounded under our distributional assumptions, we can apply FastCovFilter (Lemma 2) in conjunction with the above lemma to produce an weighting w̃ ∈ ∆n whose empirical gradient satisfies (2), in time Õ(nd). We give aOng construction for linear regression in the distributional model of Theorem 3 and a radiusless Ong construction for smooth GLMs (satisfying (2) with no dependence on R and σ = 1) in the distributional model of Theorem 1, in Corollaries 1 and 2 of the supplement. Finally, in the setting of Theorem 2 (Lipschitz GLMs), we give a reduction from constructing a noisy gradient oracle for Moreau envelopes to Õ( −1) queries of a noisy gradient oracle (cf. Corollary 3, supplement) using a characterization of the Moreau gradient as a solution to a regularized objective.
Noisy proximal oracle construction. For the remainder of this section, we assume access to Ong satisfying (2). Our second, and much more technically involved, contribution is demonstrating that acceleration is achievable under the noise model (2). Designing accelerated algorithms under noisy gradient access is extremely well-studied, and there are both strong positive results [d’A08, MS13, DG16, CDO18, MRJ19, BJL+19] as well as negative results [DGN14] showing under certain noise assumptions, accelerated gradient descent may be outperformed by unaccelerated methods. Motivated by these negative results, [PSBR20] posed the question of whether acceleration was possible under (2). We demonstrate the accelerated proximal point framework of [MS13] is amenable to such noise.
In the noiseless setting, accelerated proximal point algorithms (introduced by [Gül92]) are a reduction from optimization of a function F ? to iteratively solving proximal subproblems of the form
θ?θ̄ ← argminθ∈Rd { F ?(θ) + 1
2λ
∥∥θ − θ̄∥∥2 2 } , (3)
for some θ̄ ∈ Rd. By tuning λ we can trade off how efficiently we can solve the subproblem with how many times we need to solve them. In our accelerated algorithm, we will always set λ = 1L in proximal subproblems. We first use the (unaccelerated) robust gradient descent analysis of [PSBR20, CAT+20] to show that given a current upper bound R on ∥∥θ̄ − θ?∥∥ 2 , we can solve the proximal subproblems (3) to appropriate accuracy, yielding a “noisy proximal oracle.”
Lemma 4 (Definition 5, Proposition 7, supplement). There is an algorithm, NoisyProximalOracle, which given θ̄ ∈ Rd, R ≥ ∥∥θ̄ − θ?∥∥ 2 , and Ong giving noisy gradient estimates satisfying (2), returns θ̂ satisfying ‖θ̂ − θ? θ̄ ‖ = O(σ √ L + √ R), where θ? θ̄ is the minimizer to (3). The complexity of the algorithm is dominated by Õ(1) calls to a noisy gradient oracle Ong.
In other words, NoisyProximalOracle runs in Õ(nd) time (as expected, since (3) with λ = L−1 has condition number ≤ 2). We now discuss how to use NoisyProximalOracle in our main subroutine. Halving the radius with NoisyProximalOracle. Our accelerated algorithm runs in logarithmically many phases, each halving an upper bound on the distance to the optimizer (while above a certain noise floor, Ω(σ(κ µ ) 1 2 ), due to the additive error in (2)). We state one phase of our accelerated algorithm below as HalfRadiusAccel, assuming access to oracle Onp with the guarantees of Lemma 4.8
In Line 5, scalars {at, At} are parameters used by acceleration potential analyses, defined by recursions A0 = 0, At = La2t , At+1 = At + at+1. We remark HalfRadiusAccel is simply the proximal framework of [MS13], where we call NoisyProximalOracle in place of a (classical) proximal oracle in Line 6, and constrain Line 7. Our main technical innovation is to show the analysis of [MS13] is robust to guarantees of NoisyProximalOracle, by carefully balancing accumulated errors.
8Though our algorithm is randomized, we omit discussion of failure here for brevity. Our algorithm succeeds with high probability (with sample sizes and runtimes depending polylogarithmically on failure probabilities).
Algorithm 1 HalfRadiusAccel(θ̄, R,Onp) 1: Input: θ̄ ∈ Rd, R ≥ ∥∥θ̄ − θ?∥∥
2 with R = Ω(σ(κ µ ) 1 2 ), Onp satisfying guarantees of Lemma 4
2: Output: θ̂ ∈ Rd with ‖θ̂ − θ?‖2 ≤ 12R 3: T ← O( √ κ) for a sufficiently large constant, t← 0, θ0 ← θ̄, v0 ← θ̄ 4: while t < T do 5: yt ← AtAt+1 θt + at+1 At+1
vt 6: θt+1 ← Onp(yt, 3R) 7: vt+1 ← argminv∈B{at+1 〈yt − θt+1, v〉+ 12L ‖v − vt‖ 2 2}, where B := {v | ∥∥v − θ̄∥∥ 2 ≤ R} 8: t← t+ 1 9: end while
10: return θt
To make the analysis more modular, we first observe via characterizing proximal minimizers θ?yt (following (3)) that though iterates of Algorithm 1 may drift more than R away from θ? (due to error in Onp), they never drift too far. Formally, Lemma 13 of the supplement shows all yt satisfy ‖yt − θ?‖2 ≤ 3R, validating calls to Onp. We then give our main noisy potential bound. Lemma 5 (main potential bound, Lemma 14, supplement). For all 0 ≤ t ≤ T , define Et := F ?(θt)− F ?(θ?), Dt := 12 ‖vt − θ ?‖22, and Φt := AtEt +Dt. Then for all 0 ≤ t < T ,
Φt+1 − Φt ≤ O ( tσ √
L R+ t
√ R2 + t2σ2
L + t2 R2
) .
The potential Φt defined in Lemma 5 is typical in acceleration analyses, and indeed the (non-robust) analysis of [MS13] shows Φt+1 ≤ Φt under an exact proximal oracle. Since the coefficient of function error in the potential, At, grows as Ω( t 2
L ), it is straightforward (via strong convexity) to show HalfRadiusAccel halves the radius in O( √ κ) iterations classically. We demonstrate the analysis still goes through under the potential increase of Lemma 5 in the same number of iterations (up to constant factors), concluding analysis of HalfRadiusAccel. Our final algorithm iterates HalfRadiusAccel until reaching the noise floor, requiring Õ( √ κ) iterations dominated by calls to Ong.
Discussion. Our result builds upon recent studies of the noise tolerance of proximal point methods, e.g. [BJL+19], who used a model with fixed additive gradient error, and [CJJS21], who used noisy solutions to (3) with fixed additive error. Our model crucially tolerates both additive and multiplicative guarantees for both gradient and proximal solution estimates. An interesting open question is whether a similar result is obtainable without a proximal point abstraction, e.g. by extending [CDO18]. We view our result as a proof-of-concept that acceleration is possible under (2). We believe a unified study of acceleration under models encompassing (2) warrants further exploration, and defer it to interesting future work.

heading: 4 Improved statistical rates for robust linear regression
text: We now overview the algorithm that achieves Theorem 4, and describe its correctness proof. We define the following, strong deterministic condition on the uncorrupted points. Assumption 1 (deterministic regularity for linear regression). Let be sufficiently small, and let r ∈ (0, 2). Assume {(Xi, yi)}i∈[n] ⊂ Rd × R is ( , r)-good for linear regression (or ( , r)-good if context is clear), which means there is a partition [n] = G ∪B with |G| ≥ (1− )n which satisfies:
1. For any w ∈ ∆n with ‖wG‖1 ≥ 1− 4 , 1 2Σ ? CovwG(X) 32Σ ?.
2. There is a constant Cest such that for all θ ∈ Rd and all -saturated w ∈ ∆n, there exists a G′ ⊆ G satisfying |G′| ≥ (1− r)|G| so that if we let w̃ := wG′‖wG′‖1 ,
‖∇Fw̃(θ)−∇F ?(θ)‖2 ≤ Cest √ L (σ + ‖θ − θ?‖Σ?) , (4)∥∥Covw̃ ({gi(θ)}i∈G)∥∥op ≤ CestL(‖θ − θ?‖2Σ? + σ2) . (5)
3. There is a constantCub such that ∑ i∈G 1 |G|fi(θ ?) = 12|G| ∑ i∈G (〈Xi, θ?〉 − yi) 2 ≤ Cubσ2.
The following claim establishes that Assumption 1 (up to constants in definitions) holds with good probability under the statistical model in Theorem 4. We give the proof in the supplement. Proposition 1 (Proposition 2, supplement). Let α ≥ 1 and > 0 be sufficiently small. Let {(Xi, yi)}i∈[n] ⊂ Rd × R be -corrupted samples from DXy as in Theorem 4. Then, if
n = O
( dα2 log d
4 + d2α log(d/ ) 2
) ,
the set {(Xi, yi)}i∈[n] is (2 , 2 α )-good for linear regression with probability at least 9 10 .
Through this section, we will assume that our set of points is (2 , 2 α )-good, where we choose α = O ( max ( 1, κ log R0σ )) . Here, R0 is an initial distance bound on ‖θ0 − θ?‖Σ? . This assumption inflates the sample complexity of Proposition 1 by a factor depending on α.
Filtering under goodness. Enforcing Assumption 1, and in particular, Assumption 1.2 holds inductively introduces some additional technical difficulties for our algorithm. However, our choice of α ensures that we never remove more than 3 mass from G, and that our weights are always -saturated with respect to G ∪B, allowing for inductive application of Assumption 1. This holds because we show that our algorithm uses at most α2 distinct subsets of G for filtering. Each distinct set satisfies Lemma 1 up to a 2 2
α additive discrepancy, whose accumulation loses an additive .
Algorithms through identifiability. We now describe our algorithm, assuming Assumption 1 holds throughout. We begin with the following (slight) reinterpretation of an identifiability proof in [BP21]. Proposition 2 (Proposition 6, supplement). Let w ∈ ∆n be -saturated with respect to bipartition [n] = G ∪B, and let θ ∈ Rd. Assuming κ is sufficiently small,
‖θ − θ?‖Σ? = O (√ ( σ √ κ+ √∥∥∥Covw ({gi(θ)}i∈[n])∥∥∥op · µ−1 ) + ‖∇Fw(θ)‖(Σ?)−1 ) .
This shows if ‖∇Fw(θ)‖(Σ?)−1 and ‖Covw({gi(θ)}i∈[n])‖op are simultaneously bounded, then we obtain a distance bound to θ?. We can bound the former at fixed w by setting θ to minimize Fw, and the latter at fixed θ by filtering w via FastCovFilter. This introduces a chicken-and-egg problem of obtaining both bounds at once; we accomplish this via an alternating approach, and argue fast termination using a third potential monotone under both subroutines, the function value Fw(θ).
Halving the radius. We design a procedure, HalfRadiusLinReg, with the following guarantee. Suppose w̄ ∈ ∆n is -saturated and θ̄ ∈ Rd, and we knowR = Ω(σ) with the promise ∥∥θ̄ − θ?∥∥ Σ? ≤ R. HalfRadiusLinReg then returns a new θ with ‖θ − θ?‖Σ? ≤ 1 2R. We do so by using saturated weights to guide a potential analysis via Proposition 2. Before stating HalfRadiusLinReg, we give subroutines achieving our earlier stated goals. The first is an approximate optimizer. Definition 1. We call OERM a γ-approximate ERM oracle if on input F : Rd → R it returns a point θ̂ such that F (θ̂)− F (θ?F ) ≤ γ, for θ?F := argminθ∈RdF (θ).
The second controls the initial error, which we use to yield distance bounds via strong convexity. Lemma 6 (Lemma 6, supplement). FunctionFilter takes as input -saturated w ∈ ∆n, θ ∈ Rd, and R ≥ ‖θ − θ?‖Σ? , and produces -saturated w′ ∈ ∆n with Fw′(θ) ≤ 2Cub(σ2 +R2), in time Õ(nd).
FunctionFilter is a straightforward implementation of the filtering paradigm using Assumption 1.3, and we formally define it in the supplementary material. We can now state (a simplified version of) HalfRadiusLinReg, which iteratively filters weights initially set to uniform. It calls FunctionFilter to preprocess, and then alternates calls to OERM and FastCovFilter, which respectively decrease the gradient norm and operator norm in Proposition 1. It finally terminates when the function value decrease between consecutive iterations is small, whence we conclude an OERM call was not actually necessary for Proposition 2 to imply progress. Lemma 7 (Lemma 7, supplement). HalfRadiusLinReg correctly returns (w, θ) such that w ∈ ∆n is -saturated and ‖θ − θ?‖Σ? ≤ 1 2R in κ calls to OERM and Õ(ndκ) additional time.
Algorithm 2 HalfRadiusLinReg(X, y, w̄, θ̄, R,OERM)
1: Input: Dataset X = {Xi}i∈[n] ∈ Rn×d, y = {yi}i∈[n] ∈ Rn, satisfying Assumption 1, - saturated w̄ ∈ ∆n with respect to bipartition [n] = G ∪B with X>diag (w̄) X 8LI, θ̄ ∈ Rd
with ∥∥θ̄ − θ?∥∥
Σ? ≤ R for R = Ω(σ), O(σ
2
κ )-approximate ERM oracle OERM. 2: Output: Saturated w with respect to G ∪B, and θ with ‖θ − θ?‖Σ? ≤ 1 2R. 3: t← 0, w(0) ← FunctionFilter(w, θ̄, R), ∆0 ←∞ 4: while ∆t = Ω(R 2
κ ) do 5: θ(t) ← OERM(Fw(t)) 6: w(t+1) ← FastCovFilter({gi ( θ(t) ) }i∈[n], w(t), O(LR2))
7: ∆t+1 ← Fw(t+1) ( θ(t) ) − Fw(t+1) ( θ(t+1) ) 8: t← t+ 1 9: end while
10: return ( w(t), θ(t−1) ) We prove Lemma 7 by using the proof strategy outlined earlier: first, strong convexity, our initial function error bound, and monotonicity of function error ensure that ∥∥θ(t) − θ?∥∥ Σ?
= O(R) throughout. Next, if ∆t is ever too small, the gradient norm in Proposition 2 was small even before we ran the loop, so we could have terminated with the prior θ(t). Finally, the whole algorithm runs O(κ) times, since each loop decreases function error (initially at O(R2)) substantially.
Putting it all together. We iteratively apply Lemma 7 O(log R0σ ) times, at which point we are within distanceO(σ) in the Σ? norm. To achieve distanceO(σ √ κ ), the natural bottleneck of Proposition 2, we define LastPhase, a variant of Algorithm 2 with a more stringent termination condition in Line 4, and use a similar argument to show it iterates O( 1 ) times before concluding. Finally, we verify that throughout we only filtered against at most O(κ log R0σ + 1 ) distinct sets, where the summands are due to calls to HalfRadiusLinReg and LastPhase respectively. This fulfills the premise of using α2 distinct sets for the α defined at the beginning of the section, and concludes our analysis.

heading: Acknowledgments
text: KT is supported by NSF Grant CCF-1955039 and the Alfred P. Sloan Foundation.

Review summaries:
The paper develops fast algorithms for robust stochastic optimization in the presence of strong contamination model, where the objective function belongs to one of the following: linear regression, (regularized) generalized linear models, and (regularized) Moureau Envelope. The main contributions of the paper are as follows: 1. Building on the gradient descent framework of [PSBR20], the paper develops accelerated gradient descent method using a noisy gradient oracle. This allows them to get better dependence on the condition number in the all of the mentioned models (with similar statistical rate as existing works). 2. Inspired by the identifiability proof of [BP21], the paper gives a practical algorithm for linear regression that has better statistical rates than existing works (at the cost of additional sample complexity). This is the first work that achieves this rate without certifiable hypercontractivity.