# Robust Regression Revisited:
Acceleration and Improved Estimation Rates

Arun Jambulapati Stanford University<EMAIL_ADDRESS>Jerry Li Microsoft
Research<EMAIL_ADDRESS>Tselil Schramm Stanford University,
<EMAIL_ADDRESS>Kevin Tian Stanford University<EMAIL_ADDRESS>
###### Contents

  1. 1 Introduction
    1. 1.1 Our results
    2. 1.2 Prior work
    3. 1.3 Techniques
  2. 2 Preliminaries
    1. 2.1 Notation
    2. 2.2 Our statistical models
    3. 2.3 Linear regression
    4. 2.4 Regularity assumptions: Lipschitz and smooth stochastic optimization
    5. 2.5 Robustly decreasing the covariance operator norm
  3. 3 Linear regression
    1. 3.1 Filtering under $(\epsilon,\frac{\epsilon^{2}}{\alpha})$-goodness
    2. 3.2 Identifiability proof for linear regression
    3. 3.3 Halving the distance to $\theta^{\star}$
    4. 3.4 Last phase analysis
    5. 3.5 Full algorithm
  4. 4 Robust acceleration
    1. 4.1 Noisy gradient oracle
    2. 4.2 Halving the distance to $\theta^{\star}_{F}$
    3. 4.3 Full accelerated algorithm
  5. 5 Lipschitz generalized linear models
    1. 5.1 Noisy gradient oracle for the Moreau envelope
    2. 5.2 Accelerated optimization of the regularized Moreau envelope
  6. A Deferred proofs from Section 2
    1. A.1 Proof of Proposition 2
    2. A.2 Proof of Proposition 5

##  1 Introduction

Parameter estimation in generalized linear models, 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_ $\gamma:\mathbb{R}^{2}\to\mathbb{R}$ and a
dataset of covariates and labels
$\\{(X_{i},y_{i})\\}_{i\in[n]}\subset\mathbb{R}^{d}\times\mathbb{R}$ drawn
from an underlying distribution $\mathcal{D}_{Xy}$, the problem of statistical
(generalized linear) regression asks to

$\text{estimate
}\theta^{\star}:=\textup{argmin}_{\theta\in\mathbb{R}^{d}}\left\\{\operatorname*{\mathbb{E}}_{(X,y)\sim\mathcal{D}_{Xy}}\left[\gamma(\left\langle\theta,X\right\rangle,y)\right]\right\\}.$
(1)

For example, when $\gamma(v,y)=\frac{1}{2}(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 assumptions, since we only have access to samples from
$\mathcal{D}_{Xy}$ rather than the underlying distribution itself (see e.g.
[BP21] for tractability results in the linear regression setting).

However, in modern settings, these strong distributional assumptions may fail
to hold. In practically relevant settings, regression is often performed on
massive datasets, where the data comes from a poorly-understood distribution
and has not been thoroughly vetted or cleaned of outliers. This has prompted
the study of highly robust regression. In this work, we study the problem of
regression (1) in the _strong contamination model_. In this model, we assume
the data points we receive are independently drawn from $\mathcal{D}_{Xy}$,
but that an arbitrary $\epsilon$-fraction of the samples are then
adversarially contaminated or replaced. The strong contamination model has
recently drawn 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 $\mathcal{D}_{Xy}$ 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.

###  1.1 Our results

We give multiple _nearly-linear time algorithms_ 111Throughout, we reserve the
description “nearly-linear” for runtimes scaling linearly in the dataset size
$nd$, and polynomially in $\epsilon^{-1}$ and the condition number, up to a
polylogarithmic overhead in problem parameters. 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 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. In this work, we show that both approaches
can be sped up dramatically, and give two complementary types of algorithms
within these frameworks.

#### Robust acceleration.

Our first contribution is to demonstrate that within the noisy gradient
estimation framework for minimizing well-conditioned regression problems of
the form (1), an _accelerated_ rate of optimization can be achieved, answering
an open question asked by [PSBR20]. We demonstrate the following result for
smooth statistical regression problems, where we assume the uncorrupted data
is drawn from $\mathcal{D}_{Xy}$ with marginals $\mathcal{D}_{X}$ and
$\mathcal{D}_{y}$, $\mathcal{D}_{X}$ has support in $\mathbb{R}^{d}$, and
$\widetilde{O}$ hides polylogarithmic factors in problem parameters (cf.
Section 2.1 for technical definitions).

######  Theorem 1 (informal, see Theorem 7).

Suppose $\gamma:\mathbb{R}^{2}\times\mathbb{R}$ is such that
$\gamma_{y}(v):=\gamma(v,y)$ is convex and has (absolute) first and second
derivatives at most $1$ for all $y$ in the support of $\mathcal{D}_{y}$, and
$\mathcal{D}_{X}$ has second moment matrix $\boldsymbol{\Sigma}^{\star}\preceq
L\cdot\mathbf{I}$. For some $\mu\geq 0$, let $\kappa=\max(1,\frac{L}{\mu})$
and let

$\theta^{\star}_{\textup{reg}}:=\textup{argmin}_{\theta\in\mathbb{R}^{d}}\left\\{\operatorname*{\mathbb{E}}_{(X,y)\sim\mathcal{D}_{Xy}}\left\\{\gamma\left(\left\langle\theta,X\right\rangle,y\right)\right\\}+\frac{\mu}{2}\left\lVert\theta\right\rVert_{2}^{2}\right\\}$

be the solution to the true regularized statistical regression problem.222To
simplify our 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 regression 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. There is an algorithm that given $n:=\widetilde{O}(\frac{d}{\epsilon})$
$\epsilon$-corrupted samples from $\mathcal{D}_{Xy}$, for $\epsilon\kappa^{2}$
at most an absolute constant, runs in time $\widetilde{O}(nd\sqrt{\kappa})$
and obtains $\theta$ with
$\left\lVert\theta-\theta^{\star}_{\textup{reg}}\right\rVert_{2}=O\left(\sqrt{\frac{\kappa\epsilon}{\mu}}\right)$
with probability at least $1-\delta$.

A canonical example of a link function $\gamma$ satisfying the assumptions of
Theorem 1 is the logit function $\gamma(v,y)=\log(1+\exp(-vy))$, when the
labels $y$ are $\pm 1$. To contextualize Theorem 1, the earlier work [PSBR20]
obtains a similar statistical guarantee in its setting, using
$\widetilde{O}(\kappa)$ calls to a _noisy gradient oracle_ , which they
implement via a subroutine inspired by works on robust mean estimation. At the
time of its initial dissemination, nearly-linear time robust mean estimation
algorithms were not known; 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 mean
estimation techniques of [CDG19], and gave an algorithm running in time
$\widetilde{O}(nd\kappa\epsilon^{-6})$. Theorem 1 represents an improvement to
these results on two fronts: we apply tools from the mean-estimation algorithm
of [DHL19] to remove the $\text{poly}(\epsilon^{-1})$ runtime dependence for a
general class of regression problems, and we achieve an iteration count of
$\widetilde{O}(\sqrt{\kappa})$, matching the accelerated gradient descent
runtime of [Nes83] for smooth optimization in the non-robust setting.

We demonstrate the generality of our acceleration framework by demonstrating
that it applies to optimizing the _Moreau envelope_ for Lipschitz, but
possibly non-smooth, link functions $\gamma$; a canonical example of such a
function is the hinge loss $\gamma(v,y)=\max(0,1-vy)$ with $\pm 1$ labels,
used in training support vector machines. The Moreau envelope is a well-
studied smooth approximation to a non-smooth function which everywhere
additively approximates the original function if it is Lipschitz (see e.g.
[Sho97]), and in the non-robust setting many state-of-the-art rates for
Lipschitz optimization are known to be attained by accelerated optimization of
an appropriate Moreau envelope [TJNO20]. We show that even without explicit
access to the Moreau envelope, we can obtain approximate minimizers to it
through our robust acceleration framework.

######  Theorem 2 (informal, see Theorem 8).

Suppose $\gamma:\mathbb{R}^{2}\times\mathbb{R}$ is such that
$\gamma_{y}(v):=\gamma(v,y)$ is convex and has (absolute) first derivative at
most $1$ for all $y$ in the support of $\mathcal{D}_{y}$, and
$\mathcal{D}_{X}$ has bounded second moment matrix. For some $\mu,\lambda\geq
0$, let $\kappa=\max(1,\frac{1}{\lambda\mu})$ and let

$\displaystyle\theta^{\star}_{\textup{env}}:=\textup{argmin}_{\theta\in\mathbb{R}^{d}}\left\\{F^{\star}_{\lambda}(\theta)+\frac{\mu}{2}\left\lVert\theta\right\rVert_{2}^{2}\right\\}\textup{,
where
}F^{\star}_{\lambda}(\theta):=\inf_{\theta^{\prime}}\left\\{F^{\star}(\theta^{\prime})+\frac{1}{2\lambda}\left\lVert\theta-\theta^{\prime}\right\rVert_{2}^{2}\right\\}$
$\displaystyle\textup{ is the Moreau envelope of
}F^{\star}(\theta):=\operatorname*{\mathbb{E}}_{(X,y)\sim\mathcal{D}_{Xy}}\left\\{\gamma\left(\left\langle\theta,X\right\rangle,y\right)\right\\}.$

There is an algorithm that given $n:=\widetilde{O}(\frac{d}{\epsilon})$
$\epsilon$-corrupted samples from $\mathcal{D}_{Xy}$, for $\epsilon\kappa^{2}$
at most an absolute constant, runs in time
$\widetilde{O}(\frac{nd\sqrt{\kappa}}{\epsilon})$ and obtains $\theta$ with
$\left\lVert\theta-\theta^{\star}_{\textup{reg}}\right\rVert_{2}=O\left(\sqrt{\frac{\kappa\epsilon}{\mu}}\right)$
with probability at least $1-\delta$.

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.

The specific problem of robust linear regression is perhaps the most
ubiquitous example of statistical regression [KKM18, KKK19, DKS19, ZJS20,
CAT+20, BP21]. Amongst the algorithms developed for this problem, the only
nearly-linear time algorithm is the recent work of [CAT+20]. For an instance
of robust linear regression with noise variance bounded by $\sigma^{2}$ and
covariate second moment matrix
$\boldsymbol{\Sigma}^{\star}:=\operatorname*{\mathbb{E}}_{X\sim\mathcal{D}_{X}}[XX^{\top}]$,
the algorithms of [DKK+19, PSBR20, CAT+20] attain distance to the true
regression minimizer $\theta^{\star}$ scaling as $\sigma\kappa\sqrt{\epsilon}$
in the $\boldsymbol{\Sigma}^{\star}$ norm (the “Mahalanobis distance”) under a
bounded $4^{\textup{th}}$ moment assumption. We measure error in the
$\boldsymbol{\Sigma}^{\star}$ norm as it is scale invariant and the natural
norm in which to measure the underlying (quadratic) statistical regression
error.333Some prior works gave $\ell_{2}$ norm guarantees, which we have
translated for comparison. 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 directly applying our robust acceleration framework
leads to a similar estimation guarantee as [DKK+19, PSBR20, CAT+20] under the
same assumptions.

######  Theorem 3 (informal, see Theorem 6).

Suppose $\mathcal{D}_{X}$ is a $2$-to-$4$ hypercontractive distribution with
second moment matrix
$\boldsymbol{\Sigma}^{\star}=\operatorname*{\mathbb{E}}_{X\sim\mathcal{D}_{X}}[XX^{\top}]$
satisfying $\mu\cdot\mathbf{I}\preceq\boldsymbol{\Sigma}^{\star}\preceq
L\cdot\mathbf{I}$, and $y\sim\mathcal{D}_{y}$ is generated as
$\left\langle\theta^{\star},X\right\rangle+\delta$, for
$\delta\sim\mathcal{D}_{\delta}$ with variance at most $\sigma^{2}$. Let
$\kappa:=\frac{L}{\mu}$. There is an algorithm, $\mathsf{RobustAccel}$, that
given $n:=\widetilde{O}(\frac{d}{\epsilon})$ $\epsilon$-corrupted samples from
$\mathcal{D}_{Xy}$, for $\epsilon\kappa^{2}$ at most an absolute constant,
runs in time $\widetilde{O}(nd\sqrt{\kappa})$ and obtains $\theta$ with
$\left\lVert\theta-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}=O(\sigma\kappa\sqrt{\epsilon})$
with probability $1-\delta$.

We give a formal definition of $2$-to-$4$ hypercontractivity in Section 2.1;
as a lower bound of [BP21] shows, attaining estimation rates for robust linear
regression scaling polynomially in $\epsilon$ is impossible under only bounded
second moments, and such a $4^{\textup{th}}$ moment bound is the minimal
assumption under which robust estimation is known to be possible. Theorem 6
matches the distribution assumptions and error of [CAT+20], while obtaining an
accelerated runtime.

Interestingly, under the $4^{\textup{th}}$ moment bound used in Theorem 6,
[BP21] showed that the information-theoretically optimal rate of estimation in
the $\boldsymbol{\Sigma}^{\star}$ norm is _independent_ of $\kappa$, and
presented a matching upper bound under an analogous, but more stringent,
distributional assumption.444The algorithm of [BP21] requires
$\mathcal{D}_{X}$ to be _certifiably hypercontractive_ , an algebraic
condition frequently required by the sum-of-squares algorithmic paradigm to
apply to robust statistical estimation problems. We remark that thus far
robust linear regression algorithms have broadly fallen under two categories:
The first category (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(\sigma\kappa\sqrt{\epsilon})$ and the requirement that
$\epsilon\kappa^{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, and
attaining improved statistical performance compared to Theorem 3 while only
requiring $\epsilon\kappa=O(1)$.

######  Theorem 4 (informal, see Theorem 5).

Suppose $\mathcal{D}_{X}$ is a $2$-to-$4$ hypercontractive distribution with
second moment matrix
$\boldsymbol{\Sigma}^{\star}=\operatorname*{\mathbb{E}}_{X\sim\mathcal{D}_{X}}[XX^{\top}]$
satisfying $\mu\cdot\mathbf{I}\preceq\boldsymbol{\Sigma}^{\star}\preceq
L\cdot\mathbf{I}$, and $y\sim\mathcal{D}_{y}$ is generated as
$\left\langle\theta^{\star},X\right\rangle+\delta$, for
$\delta\sim\mathcal{D}_{\delta}$, a $2$-to-$4$ hypercontractive distribution
with variance at most $\sigma^{2}$. Let $\kappa:=\frac{L}{\mu}$. There is an
algorithm, $\mathsf{FastRegression}$, that given
$n:=\widetilde{O}((\frac{d^{2}}{\epsilon^{3}}+\frac{d}{\epsilon^{4}}))$
$\epsilon$-corrupted samples from $\mathcal{D}_{Xy}$, for $\epsilon\kappa$ at
most an absolute constant, uses $\widetilde{O}(\frac{1}{\epsilon})$ calls to
an empirical risk minimization routine555The empirical risk minimization
algorithm used is up to the practitioner; its runtime will never scale worse
than $\widetilde{O}(nd\sqrt{\kappa})$ by applying (non-robust) accelerated
gradient descent, but can be substantially better if recent advances in
stochastic gradient methods are used, e.g. [Zhu17]. and
$\widetilde{O}(\frac{nd}{\epsilon})$ additional runtime, and obtains $\theta$
with
$\left\lVert\theta-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}=O(\sigma\sqrt{\kappa\epsilon})$
with probability at least $\frac{9}{10}$.

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 $\kappa$) 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 framework. Our sample complexity for Theorem 4
dramatically improves that of [DKK+17]’s original linear regression algorithm
in the Sever framework for moderate $\epsilon$, which used
$\widetilde{O}(\frac{d^{5}}{\epsilon^{2}})$ samples (in addition to beating
their weaker error guarantee). We elaborate on these points further in Section
1.3.

###  1.2 Prior work

We give a general overview contextualizing our work in this section, and defer
the comparison of specific technical components we develop in this work to
relevant discussions.

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 full survey of this field is beyond our scope, and we defer a more
comprehensive overview to [DK19, Li18, Ste18].

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, 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 intense 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
$\epsilon^{-6}$ runtime overhead removed using more sophisticated mean
estimation techniques.

###  1.3 Techniques

We now describe the techniques we use to obtain the accelerated rates of
Theorems 1, 2, and 3 as well as the robust linear regression algorithm of
Theorem 4.

#### Robust acceleration.

Our robust acceleration framework is based on the following abstract
formulation of an optimization problem: there is an unknown function
$F^{\star}:\mathbb{R}^{d}\to\mathbb{R}$ with minimizer $\theta^{\star}$ which
is $L$-smooth and $\mu$-strongly convex, and we wish to estimate
$\theta^{\star}$, but our only mode of accessing $F^{\star}$ is through a
_noisy gradient oracle_ $\mathcal{O}_{\textup{ng}}$. Namely, for some
$\sigma$, $\epsilon$, we can query $\mathcal{O}_{\textup{ng}}$ at any point
$\theta\in\mathbb{R}^{d}$ with an upper bound
$R\geq\left\lVert\theta-\theta^{\star}\right\rVert_{2}$ and receive an
estimate $G(\theta)$ such that

$\left\lVert G(\theta)-\nabla
F^{\star}(\theta)\right\rVert_{2}=O\left(\sqrt{L\epsilon}\sigma+L\sqrt{\epsilon}R\right).$
(2)

In other words, we receive gradients perturbed by both fixed additive noise,
and multiplicative noise depending on the distance to $\theta^{\star}$. The
prior works [DKK+19, PSBR20, CAT+20] observed that by using tools from robust
mean estimation, appropriate noisy gradient oracles could be constructed for
the functions
$F^{\star}(\theta)=\operatorname*{\mathbb{E}}_{(X,y)\sim\mathcal{D}_{Xy}}[\sigma(\left\langle\theta,X\right\rangle-y)]$
arising from the distributional assumptions in Theorems 1, 2, and 3. Our first
contribution is speeding up the implementation of $\mathcal{O}_{\textup{ng}}$
to run in nearly-linear time $\widetilde{O}(nd)$, leveraging recent advances
by [DHL19] for robust mean estimation.

Our second, and more technically involved, contribution is demonstrating that
accelerated runtimes are achievable under the noisy gradient oracle access
model of (2). Designing accelerated algorithms under noisy gradient access is
an extremely well-studied problem, and there are both strong positive results
[d’A08, MS13, DG16, CDO18, MRJ19, BJL+19] as well as negative results [DGN14]
showing that under certain noise models, accelerated gradient descent may be
outperformed by unaccelerated methods. Indeed, it was asked (motivated by
these negative results) as an open question in [PSBR20] whether an accelerated
rate was possible under the noise model (2).

Our accelerated algorithm runs in logarithmically many phases, where we halve
the distance to the optimizer (while it is above a certain noise floor
depending on the additive error in (2)) in each phase. The subroutine we
design for implementing each phase is a robust accelerated “outer loop”
tolerant to noisy gradient access in the manner provided by our oracle
$\mathcal{O}_{\textup{ng}}$. By carefully balancing the accuracy of subproblem
solutions, the multiplicative error in our gradient estimates within the
accelerated outer loop, and the drift of the phase’s iterates (which may
venture further from $\theta^{\star}$ than our initial iterate upper bound),
we show that above the noise floor we can halve the distance to
$\theta^{\star}$ in $\widetilde{O}(\sqrt{\kappa})$ queries to
$\mathcal{O}_{\textup{ng}}$; recursing on this guarantee yields our complete
algorithm.

To obtain Theorem 2, we demonstrate that for Lipschitz functions $F^{\star}$
admitting a _radiusless_ noisy gradient oracle, i.e. one which satisfies (2)
with no dependence on $R$, we can further efficiently construct a noisy
gradient oracle for the Moreau envelope of $F^{\star}$ using projected
subgradient descent. This construction enables applying our robust accelerated
method to Lipschitz regression problems.

Our acceleration framework crucially tolerates both additive and
multiplicative guarantees for gradient estimation. While it is possible that
arguments of other noisy acceleration frameworks e.g. [CDO18] may be extended
to capture our gradient noise model, we give a self-contained derivation
specialized to our specific oracle access for convenience. We view our result
as a proof-of-concept that acceleration is possible under this noise model; we
believe a unified study of acceleration under noise models encompassing (2)
warrants further exploration, and defer it to interesting future work.

#### Robust linear regression.

For the special case of linear regression, as discussed earlier, the robust
optimization methods of [DKK+19, PSBR20, CAT+20] attain Mahalanobis distance
scaling as $\sigma\kappa\sqrt{\epsilon}$ from the true minimizer. Directly
plugging in deterministic conditions proven by [CAT+20] to hold under an
appropriate statistical model into our robust gradient descent framework, we
obtain a similar guarantee (Theorem 3) at an accelerated rate. In this
technical overview, we now focus on how we obtain the improvements of Theorem
4.

At a high level, prior works lose two factors of $\sqrt{\kappa}$ in their
error guarantees because of two _norm conversions_ from the
$\boldsymbol{\Sigma}^{\star}$ norm to the $\ell_{2}$ norm: one in gradient
space, and one in parameter space. Because we do not have access to the true
covariance $\boldsymbol{\Sigma}^{\star}$, it is natural to perform both
gradient estimation and the gradient descent procedure itself in the
$\ell_{2}$ norm. When the $\ell_{2}$ guarantees of both subroutines are
converted back to the $\boldsymbol{\Sigma}^{\star}$ norm, the error rate is
lossy by a factor of $\kappa$.

We give a different approach to robust linear regression which bypasses this
barrier in parameter space, saving a factor of $\sqrt{\kappa}$ in our error
rate. In particular, we measure progress of our parameter estimates entirely
in the $\boldsymbol{\Sigma}^{\star}$ norm in our analysis, which removes the
need for an additional norm conversion. Our starting point is the following
_identifiability proof_ guarantee of [BP21], which we slightly repurpose for
our needs. Let $w\in\Delta^{n}$ be entrywise less than
$\frac{1}{n}\mathbbm{1}$ such that $\left\lVert w_{G}\right\rVert_{1}\geq
1-4\epsilon$ where $G$ is our “uncorrupted” data, and let
$\theta\in\mathbb{R}^{d}$. We demonstrate in Proposition 6 that

$\left\lVert\theta-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}=O\left(\sigma\sqrt{\kappa\epsilon}+\sqrt{\left\lVert\textup{Cov}_{w}\left(\left\\{g_{i}(\theta)\right\\}_{i\in[n]}\right)\right\rVert_{\textup{op}}\frac{\epsilon}{\mu}}+\sqrt{\epsilon}\left\lVert\nabla
F_{w}(\theta)\right\rVert_{(\boldsymbol{\Sigma}^{\star})^{-1}}\right).$ (3)

In the above display, $g_{i}(\theta)$ is the empirical gradient of the squared
loss at our $i^{\text{th}}$ data point, $\textup{Cov}_{w}(\cdot)$ is the
empirical second moment matrix of its argument under the weighting $w$, and
$F_{w}$ is the empirical risk under $w$. The guarantee (3) suggests a natural
approach to estimation: if we can _simultaneously_ verify that $\theta$ is an
approximate minimizer to $F_{w}$, and that the empirical second moment of
gradients at $\theta$ (according to $w$) are small, then we have a proof that
$\theta$ and $\theta^{\star}$ are close.

Prior work by [BP21] used this approach to obtain a polynomial-time estimator
by solving a joint optimization problem in $(w,\theta)$, via an appropriate
semidefinite program relaxation. However, in designing near-linear time
algorithms, we cannot afford to use said relaxation. This raises a chicken-
and-egg issue: for fixed $w$, it is straightforward to make $\left\lVert\nabla
F_{w}(\theta)\right\rVert_{(\boldsymbol{\Sigma}^{\star})^{-1}}$ small, by
setting $\theta$ to the empirical risk minimizer (ERM) of $F_{w}$. Likewise,
for fixed $\theta$, known _filtering_ techniques rapidly decrease
$\|\textup{Cov}_{w}(\\{g_{i}(\theta)\\})_{i\in[n]}\|_{\textup{op}}$ while
preserving most of $w_{G}$, by using that the second moment restricted to
uncorrupted points has a small operator norm as a certificate for outlier
removal. However, performing either of these subroutines to guarantee one of
our sufficient conditions passes (small operator norm or gradient norm) may
adversely affect the quality of the other.

We circumvent this chicken-and-egg problem by introducing a _third_ potential,
namely the actual function value $F_{w}(\theta)$. In particular, notice that
the two subroutines we described earlier (downweighting $w$ or setting
$\theta$ to the ERM) both decrease this third potential. Our linear regression
algorithm is an alternating procedure which iteratively filters $w$ based on
the gradients at the current $\theta$ (to make the operator norm small), and
sets $\theta$ to the ERM of $F_{w}$ (to zero out the gradient norm). We
further show that if the ERM step does not make significant function progress
(our third potential), then it was not strictly necessary to make progress
according to (3), since the gradient norm was already small. This gives a
dimension-independent bound on the number of times we could have alternated,
via tracking function progress, yielding Theorem 4.

##  2 Preliminaries

We give the notation used throughout this paper in Section 2.1, and set up the
statistical model we consider in Section 2.2. In Section 2.3, we give the
deterministic regularity assumptions used by our regression algorithm in
Section 3. In Section 2.4, we give the deterministic regularity assumptions
used by our stochastic optimization algorithms in Sections 4 and 5. Finally,
in Section 2.5, we state a nearly-linear time procedure for robustly
decreasing the operator norm of the second moment matrix of a set of vectors.
Some proofs are deferred to the appendices.

###  2.1 Notation

#### General notation.

For $d\in\mathbb{N}$ we let $[d]:=\\{j\mid j\in\mathbb{N},1\leq j\leq d\\}$.
The $\ell_{p}$ norm of a vector argument is denoted
$\left\lVert\cdot\right\rVert_{p}$, where
$\left\lVert\cdot\right\rVert_{\infty}$ is the element with largest absolute
value; when the argument is a symmetric matrix, we overload this to mean the
Schatten-$p$ norm. The all-ones vector (of appropriate dimension from context)
is denoted $\mathbbm{1}$. The (solid) probability simplex is denoted
$\Delta^{n}:=\\{w\in\mathbb{R}^{n}_{\geq 0},\left\lVert w\right\rVert_{1}\leq
1\\}$. We use $\widetilde{O}$ to suppress logarithmic factors in dimensions,
distance ratios, the problem condition number $\kappa$, the inverse corruption
parameter $\epsilon^{-1}$, and the inverse failure probability. For
$v\in\mathbb{R}^{n}$ and $S\subseteq[n]$, we let $v_{S}\in\mathbb{R}^{n}$
denote $v$ with coordinates in $[n]\setminus S$ zeroed out. For a set $S$, we
call $S_{1}$, $S_{2}$ a _bipartition_ of $S$ if $S_{1}\cap S_{2}=\emptyset$
and $S_{1}\cup S_{2}=S$.

#### Matrices.

Matrices are denoted in boldface. We denote the zero and identity matrices (of
appropriate dimension) by $\mathbf{0}$ and $\mathbf{I}$. The $d\times d$
symmetric matrices are $\mathbb{S}^{d}$, and the $d\times d$ positive
semidefinite cone is $\mathbb{S}_{\geq 0}^{d}$. For
$\mathbf{A},\mathbf{B}\in\mathbb{S}^{d}$ we write
$\mathbf{A}\preceq\mathbf{B}$ to mean
$\mathbf{B}-\mathbf{A}\in\mathbb{S}_{\geq 0}^{d}$. The largest and smallest
eigenvalue and trace of a symmetric matrix are respectively denoted
$\lambda_{\textup{max}}(\cdot)$, $\lambda_{\textup{min}}(\cdot)$, and
$\textup{Tr}(\cdot)$. The inner product on $\mathbb{S}^{d}$ is
$\left\langle\mathbf{A},\mathbf{B}\right\rangle:=\textup{Tr}(\mathbf{A}\mathbf{B})$.
For positive definite $\mathbf{M}$, we define the induced norm $\left\lVert
v\right\rVert_{\mathbf{M}}:=\sqrt{v^{\top}\mathbf{M}v}$. We use
$\left\lVert\cdot\right\rVert_{\textup{op}}$ to mean the $\ell_{2}$-$\ell_{2}$
operator norm of a matrix; when the argument is symmetric, it is synonymous
with $\lambda_{\textup{max}}$, and otherwise is the largest singular value.

#### Functions.

The gradient and Hessian of a twice-differentiable function are denoted
$\nabla$ and $\nabla^{2}$. We say differentiable
$f:\mathbb{R}^{d}\to\mathbb{R}$ is $\lambda$-Lipschitz in a quadratic norm
$\left\lVert\cdot\right\rVert_{\mathbf{M}}$ if $\left\lVert\nabla
f(\theta)\right\rVert_{\mathbf{M}^{-1}}\leq\lambda$ for all
$\theta\in\mathbb{R}^{d}$. We say twice-differentiable
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is $L$-smooth and $\mu$-strongly
convex in $\left\lVert\cdot\right\rVert_{\mathbf{M}}$ if

$\mu\mathbf{M}\preceq\nabla^{2}f(\theta)\preceq L\mathbf{M},\text{ for all
}\theta\in\mathbb{R}^{d}.$

When $\mathbf{M}$ is not specified, we assume $\mathbf{M}=\mathbf{I}$ (i.e.
the norm in question is $\ell_{2}$). For any $\mathbf{M}$, smoothness and
strong convexity imply the following bounds for all
$\theta,\theta^{\prime}\in\mathbb{R}^{d}$,

$f(\theta)+\left\langle\nabla
f(\theta),\theta^{\prime}-\theta\right\rangle+\frac{\mu}{2}\left\lVert\theta^{\prime}-\theta\right\rVert_{\mathbf{M}}^{2}\leq
f(\theta^{\prime})\leq f(\theta)+\left\langle\nabla
f(\theta),\theta^{\prime}-\theta\right\rangle+\frac{L}{2}\left\lVert\theta^{\prime}-\theta\right\rVert_{\mathbf{M}}^{2}.$

It is well-known that $L$-smoothness of function $f$ implies $L$-Lipschitzness
of the function gradient $\nabla f$, i.e. $\left\lVert\nabla f(\theta)-\nabla
f(\theta^{\prime})\right\rVert_{\mathbf{M}^{-1}}\leq
L\left\lVert\theta-\theta^{\prime}\right\rVert_{\mathbf{M}}$ for all
$\theta,\theta^{\prime}\in\mathbb{R}^{d}$. For any $f$ which is $L$-smooth and
$\mu$-strongly convex in $\left\lVert\cdot\right\rVert_{\mathbf{M}}$, with
$\theta^{*}:=\textup{argmin}_{\theta\in\mathbb{R}^{d}}f(\theta)$, it is
straightforward to show

$\frac{1}{2L}\left\lVert\nabla f(\theta)\right\rVert_{\mathbf{M}^{-1}}^{2}\leq
f(\theta)-f(\theta^{\star})\leq\frac{1}{2\mu}\left\lVert\nabla
f(\theta)\right\rVert_{\mathbf{M}^{-1}}^{2}\text{ for all
}\theta\in\mathbb{R}^{d}.$

#### Distributions.

The multivariate Gaussian distribution with mean $\mu$ and covariance
$\boldsymbol{\Sigma}$ is denoted $\mathcal{N}(\mu,\boldsymbol{\Sigma})$. For
weights $w\in\mathbb{R}_{\geq 0}^{n}$ and a set of vectors
$\mathbf{X}:=\\{X_{i}\\}_{i\in[n]}$, we let

$\mu_{w}\left(\mathbf{X}\right):=\sum_{i\in[n]}\frac{w_{i}}{\left\lVert
w\right\rVert_{1}}X_{i},\;\textup{Cov}_{w,\bar{X}}\left(\mathbf{X}\right):=\sum_{i\in[n]}\frac{w_{i}}{\left\lVert
w\right\rVert_{1}}\left(X_{i}-\bar{X}\right)\left(X_{i}-\bar{X}\right)^{\top}.$

be the empirical mean and (centered) covariance matrix; when $\bar{X}$ is not
specified, it is the zeroes vector. Draws from the uniform distribution on
$\\{X_{i}\\}_{i\in[n]}$ are denoted $X\sim_{\textup{unif}}\mathbf{X}$. We say
distribution $\mathcal{D}$ supported on $\mathbb{R}^{d}$ is $2$-to-$4$
hypercontractive with parameter $C_{2\to 4}$ if for all $v\in\mathbb{R}^{d}$,

$\operatorname*{\mathbb{E}}_{X\sim\mathcal{D}}\left[\left\langle
X,v\right\rangle^{4}\right]\leq C_{2\to
4}\operatorname*{\mathbb{E}}_{X\sim\mathcal{D}}\left[\left\langle
X,v\right\rangle^{2}\right]^{2}.$

We will refer to this property as being $C_{2\to 4}$-hypercontractive for
short; by massaging the definition, we observe $C_{2\to 4}$-hypercontractivity
is preserved under linear transformations of the distribution.

#### Filtering.

We will make much use of the following algorithmic technique, which refer to
as _filtering_. In the filtering paradigm, we have an index set $[n]$, and a
fixed, unknown bipartition $[n]=G\cup B$, $G\cap B=\emptyset$. The set $G$ is
a “good” set of indices that we wish to keep, and the set $B$ is a set of
“bad” indices which we would like to remove. The algorithm maintains a set of
weights $w\in\Delta^{n}$ (with the goal of producing a weight vector which is
close to the uniform distribution on $G$). These weights are iteratively
updated according to “scores” $\tau\in\mathbb{R}^{n}_{\geq 0}$; the goal of
filtering is to assign large scores to coordinates in $B$ and small scores to
coordinates in $G$, so that the bad coordinates can be filtered out according
to their scores. Concretely, we use the following definition.

######  Definition 1 (saturated weights).

We say weights $w\in\Delta^{n}$ are $c$-_saturated_ with respect to the
bipartition $G\cup B=[n]$ if $w\leq\frac{1}{n}\mathbbm{1}$ entrywise, and

$\left\lVert\left[\frac{1}{n}\mathbbm{1}-w\right]_{G}\right\rVert_{1}\leq\left\lVert\left[\frac{1}{n}\mathbbm{1}-w\right]_{B}\right\rVert_{1}+c.$

If $c=0$, we refer to $w$ as simply _saturated_.

In words, $w$ is saturated if its difference from the uniform distribution has
more weight on $B$ than $G$ (in the context of our algorithm, if we have
started with uniform weights and produced a saturated $w$, then we have
removed more mass from $B$ than $G$). We allow for a “fudge factor” of an
additive $c$ to relax the above definition, which will come in handy in our
linear regression applications.

######  Definition 2 (safe scores).

Suppose $G\cup B$ is a bipartition of $[n]$, and suppose $w\in\Delta^{n}$ is a
set of weights. We call a set of scores
$\tau=\\{\tau_{i}\\}_{i\in[n]}\in\mathbb{R}^{n}_{\geq 0}$ _safe with respect
to $w$_ if it satisfies

$\left\langle w_{G},\tau\right\rangle\leq\left\langle
w_{B},\tau\right\rangle.$

The following simple lemma (implicit in prior works [DKK+17, CSV17b, Li18,
Ste18]) is the crux of the filtering paradigm, relating these two definitions.

######  Lemma 1.

Suppose $w\in\Delta^{n}$ is saturated, and $\tau\in\mathbb{R}^{n}_{\geq 0}$ is
safe with respect to $w$. Defining $w^{\prime}$ by

$w^{\prime}_{i}\leftarrow\left(1-\frac{\tau_{i}}{\tau_{\max}}\right)w_{i}\text{
for all }i\in[n],\text{ and }\tau_{\max}:=\max_{i\in[n]\mid w_{i}\neq
0}\tau_{i},$

then $w^{\prime}\in\Delta^{n}$ is also saturated.

###### Proof.

If $G\cup B=[n]$ is the bipartition with good coordinates $G$, then by
definition of safe scores,

$\left\lVert\left[w-w^{\prime}\right]_{G}\right\rVert_{1}=\frac{1}{\tau_{\max}}\left\langle
w_{G},\tau\right\rangle\leq\frac{1}{\tau_{\max}}\left\langle
w_{G},\tau\right\rangle\leq\left\lVert\left[w-w^{\prime}\right]_{B}\right\rVert_{1}.$

Now, since $w$ is saturated, and since $w^{\prime}\leq
w\leq\frac{1}{n}\mathbbm{1}$ by definition of saturation,

$\left\lVert\left[\tfrac{1}{n}\mathbbm{1}-w^{\prime}\right]_{G}\right\rVert_{1}=\left\lVert\left[\tfrac{1}{n}\mathbbm{1}-w\right]_{G}\right\rVert_{1}+\left\lVert\left[w-w^{\prime}\right]_{G}\right\rVert_{1}\leq\left\lVert\left[\tfrac{1}{n}\mathbbm{1}-w\right]_{B}\right\rVert_{1}+\left\lVert\left[w-w^{\prime}\right]_{B}\right\rVert_{1}=\left\lVert\left[\tfrac{1}{n}\mathbbm{1}-w^{\prime}\right]_{G}\right\rVert_{1}.$

∎

We will also frequently using the following simple fact.

######  Lemma 2.

Suppose $w\in\Delta^{n}$ is $c$-saturated with respect to bipartition
$[n]=G\cup B$, and suppose the bad set $|B|\leq\epsilon n$. Let
$\tilde{w}=\frac{w}{\left\lVert w\right\rVert_{1}}$ be the distribution with
probabilities proportional to $w$ and
$w^{\star}_{G}=\frac{1}{|G|}\mathbbm{1}_{G}$ be uniform over $G$. Then,
$\left\lVert\tilde{w}-w^{\star}_{G}\right\rVert_{1}\leq 6\epsilon+2c$.

###### Proof.

By the definition of saturation, since there is only $\epsilon$ mass to remove
from the coordinates of $B$ on $\frac{1}{n}\mathbbm{1}$, clearly $\left\lVert
w\right\rVert_{1}\geq 1-2\epsilon-c$. By the triangle inequality, we have

$\left\lVert\tilde{w}-w^{\star}_{G}\right\rVert_{1}\leq\left\lVert\tilde{w}-w\right\rVert_{1}+\left\lVert
w-\tfrac{1}{n}\mathbbm{1}\right\rVert_{1}+\left\lVert\tfrac{1}{n}\mathbbm{1}-w^{\star}_{G}\right\rVert_{1}.$

By definition of $\tilde{w}$,
$\left\lVert\tilde{w}-w\right\rVert_{1}=1-\left\lVert w\right\rVert_{1}\leq
2\epsilon+c$. By saturation, $\left\lVert
w-\tfrac{1}{n}\mathbbm{1}\right\rVert_{1}\leq 2\epsilon+c$. Finally, since
$|G|\geq(1-\epsilon)n$,
$\left\lVert\tfrac{1}{n}\mathbbm{1}-w^{\star}_{G}\right\rVert_{1}\leq
2\epsilon$. Combining these pieces yields the claim. ∎

###  2.2 Our statistical models

In this paper, we provide provable guarantees for optimization problems
captured by the following statistical model.

######  Model 1 (stochastic optimization in the strong contamination model).

For $\mathcal{D}_{f}$ a distribution over functions
$f:\mathbb{R}^{d}\to\mathbb{R}$, our goal is to optimize
$\operatorname*{\mathbb{E}}_{f\sim\mathcal{D}_{f}}[f(\theta)]$. We are given
access to $n$ samples $\\{f_{i}\\}_{i\in[n]}$ produced as follows:

  1. 1.

Functions $\\{\tilde{f}_{i}\\}_{i\in[n]}$ are drawn independently from
$\mathcal{D}_{f}$.

  2. 2.

An arbitrary subset $B\subset[n]$ of the samples is replaced with arbitrary
functions from $\mathrm{supp}(\mathcal{D}_{f})$.

  3. 3.

For each $i\in[n]$, if $i\in B$ we observe $f_{i}$ as the corrupted sample,
otherwise, we observe $f_{i}=\tilde{f}_{i}$.

We call $B$ the _corrupted_ samples. When $|B|=\epsilon n$, we say
$\\{f_{i}\\}_{i\in[n]}$ is drawn _$\epsilon$ -corrupted from
$\mathcal{D}_{f}$_.

Throughout we use the convention that $G\cup B=[n]$ is the bipartition of
sample coordinates with $B$ the corrupted samples and $G$ the “good” samples.
For simplicity, we assume throughout that $\epsilon=\frac{|B|}{n}$ is smaller
than some globally fixed constant. We will also frequently use the notation
$g_{i}$ to mean $\nabla f_{i}$ for all $i\in[n]$. When $\mathcal{D}_{f}$ is
clear from context, we denote the “true average function” by

$F^{\star}(\theta):=\operatorname*{\mathbb{E}}_{f\sim\mathcal{D}_{f}}\left[f(\theta)\right].$

For $w\in\Delta^{n}$ and functions $\\{f_{i}\\}_{i\in[n]}$, the (unnormalized)
weighted empirical average function is

$F_{w}(\theta):=\sum_{i\in[n]}w_{i}f_{i}(\theta).$ (4)

We will use $w^{\star}_{G}$ to denote the uniform distribution over $G$,
$w^{\star}_{G}:=\frac{1}{|G|}\mathbbm{1}_{G}$, and we use $F_{G}$ as shorthand
for the function $F_{w^{\star}_{G}}$. The goal of robust parameter estimation
is to estimate the true optimizer, which we always denote by
$\theta^{\star}:=\textup{argmin}_{\theta\in\mathbb{R}^{d}}F^{\star}(\theta)$.
For example, the problem estimating the mean of $\mathcal{D}$ can be expressed
by choosing $f_{i}(\theta)=\frac{1}{2}\left\lVert\theta-
X_{i}\right\rVert_{2}^{2}$ for $X_{i}\sim\mathcal{D}$. In the uncorrupted
setting (i.e. $\epsilon=0$), a typical strategy (given reasonable regularity
assumptions on $\mathcal{D}_{f}$) is to choose the estimator
$\theta_{G}:=\textup{argmin}_{\theta\in\mathbb{R}^{d}}F_{G}(\theta)$. The
challenge is to obtain comparable estimation performance to $\theta_{G}$ which
is robust to an $\epsilon$-fraction of unknown corruptions.

Our focus in this paper is optimizing generalized linear models. In
particular, throughout we will work only with $\mathcal{D}_{f}$ of the
following form.

######  Model 2 (generalized linear model).

A generalized linear model is a distribution over functions
$f_{i}:\mathbb{R}^{d}\to\mathbb{R}$ which is defined by a joint distribution
$\mathcal{D}_{Xy}$ over pairs
$\\{X_{i},y_{i}\\}\in\mathbb{R}^{d}\times\mathbb{R}$ and a _link function_
$\gamma:\mathbb{R}^{2}\to\mathbb{R}$, so that samples
$f_{i}\sim\mathcal{D}_{f}$ are generated as

$f_{i}(\theta):=\gamma\left(\left\langle
X_{i},\theta\right\rangle,y_{i}\right),\quad\text{ for
}(X_{i},y_{i})\sim\mathcal{D}_{Xy}.$ (5)

Note that observing $\\{f_{i}\\}_{i\in[n]}$ is equivalent to observing the
dataset $\\{X_{i},y_{i}\\}_{i\in[n]}$ when $\gamma$ is known.

For instance, when $\gamma(v,y)=\frac{1}{2}(v-y)^{2}$, this is the problem of
(statistical) linear regression. Further, when
$\gamma(v,y)=\log(1+\exp(-vy))$, our problem is logistic regression, and when
$\gamma(v,y)=\max(0,1-vy)$, it is fitting a support vector machine. We refer
to the $X$ and $y$ marginals over $\mathcal{D}_{Xy}$ respectively by
$\mathcal{D}_{X}$ and $\mathcal{D}_{y}$, and we denote
$\boldsymbol{\Sigma}^{\star}:=\operatorname*{\mathbb{E}}_{X\sim\mathcal{D}_{X}}\left[XX^{\top}\right]$
when $\mathcal{D}_{X}$ is clear from context.

###  2.3 Linear regression

In Section 3 and (part of) Section 4, we develop algorithms for the well-
studied special case of the generalized linear model, Model 2 wherein
$\gamma(v,y)=\frac{1}{2}(v-y)^{2}$, i.e. a statistical variant of linear
regression. We obtain guarantees under the following model and regularity
assumptions for $\mathcal{D}_{Xy}$.

######  Model 3 (distributional regularity for linear regression).

Given distributions $\mathcal{D}_{X}$ and $\mathcal{D}_{\delta}$ over
$\mathbb{R}^{d},\mathbb{R}^{d}$ respectively, the distribution
$\mathcal{D}_{Xy}$ over $\mathbb{R}^{d}\times\mathbb{R}$ is sampled as
follows: for an underlying vector $\theta^{\star}\in\mathbb{R}^{d}$,
independently sample $X\sim\mathcal{D}_{X}$ and
$\delta\sim\mathcal{D}_{\delta}$, and set
$y\leftarrow\left\langle\theta^{\star},X\right\rangle+\delta$. Further,
$\mathcal{D}_{X}$ and $\mathcal{D}_{\delta}$ satisfy the following regularity
assumptions.

  1. 1.

For
$\boldsymbol{\Sigma}^{\star}=\operatorname*{\mathbb{E}}_{\mathcal{D}_{X}}[XX^{\top}]$
and $0<\mu<L$, we have $\mu\mathbf{I}\preceq\boldsymbol{\Sigma}^{\star}\preceq
L\mathbf{I}$.

  2. 2.

$\mathcal{D}_{X}$ is $C_{2\to 4}$-hypercontractive for a constant $C_{2\to
4}$.

  3. 3.

$\mathcal{D}_{\delta}$ is a $C_{2\to 4}$-hypercontractive distribution with
mean zero and variance $\leq\sigma^{2}$.

For $\\{(X_{i},y_{i})\\}_{i\in[n]}$ (in particular, overloading to include
$i\in B$), we use the notation $\delta_{i}:=y_{i}-\left\langle
X_{i},\theta^{\star}\right\rangle$. We also use the following notation when
discussing linear regression:

$f_{i}(\theta):=\frac{1}{2}\left(\left\langle X_{i},\theta\right\rangle-
y_{i}\right)^{2},\;g_{i}(\theta):=\nabla f_{i}(\theta)=X_{i}\left(\left\langle
X_{i},\theta\right\rangle-y_{i}\right).$ (6)

We will denote the condition number of $\boldsymbol{\Sigma}^{\star}$ by
$\kappa:=\frac{L}{\mu}$ throughout. Under Model 3, it is immediate from the
first-order optimality condition that for

$F^{\star}(\theta):=\operatorname*{\mathbb{E}}_{X,y\sim\mathcal{D}_{Xy}}\left[\frac{1}{2}\left(\left\langle
X,\theta\right\rangle-y\right)\right]^{2},$

the optimizer $\textup{argmin}_{\theta\in\mathbb{R}^{d}}F^{\star}(\theta)$ is
exactly $\theta^{\star}$.

In our setting, following the description in Section 2.2 we independently draw
$\\{(X_{i},y_{i})\\}_{i\in G}\sim\mathcal{D}_{Xy}$ for $|G|=(1-\epsilon)n$ and
observe $\\{(X_{i},y_{i})\\}_{i\in[n]}$ where $[n]=G\cup B$ and
$\\{(X_{i},y_{i})\\}_{i\in B}$ are arbitrarily chosen. We will frequently
refer to $\\{X_{i}\\}_{i\in[n]}$ as $\mathbf{X}\in\mathbb{R}^{n\times d}$.
Under Model 3, recent work [BP21] obtained the following results.

######  Proposition 1 ([BP21], Theorem 1.7, Theorem 1.9, Theorem 1.2).

For Models 1 and 3, the minimax optimal error rate for estimators
$\hat{\theta}$ is

$\left\lVert\hat{\theta}-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}=O\left(\sigma\epsilon^{\frac{3}{4}}\right).$

When the distribution $\mathcal{D}_{X}$ is further _certifiably
hypercontractive_ in the sum-of-squares proof system, there is a
$\textup{poly}(d)$-time estimator requiring $\textup{poly}(d)$ samples
achieving this rate with high probability. Moreover, without the
hypercontractivity condition in Model 3, even when $\mu,L=\Theta(1)$ it is
information-theoretically impossible to attain an error rate depending
polynomially on $\epsilon$.

The algorithmic result of [BP21] (and all known techniques with error rate
$\ll\sqrt{\epsilon}$) crucially requires that the distributions are sum-of-
squares certifiably hypercontractive, which is a stronger assumption than
(standard) hypercontractivity. There is evidence that the problem of
certifying $2\to 4$ hypercontractivity is computationally intractable in
general (under e.g. the small-set expansion hypothesis, see [BBH+12, BGG+19]).
Even for certifiably hypercontractive distributions, known algorithms require
use spectral estimators of higher-order moment matrices, and thus more samples
and increased runtime complexity. Hence, error $\sqrt{\epsilon}$ is a standing
barrier for fast algorithms.

In Section 3, whenever we discuss robust linear regression we operate in
Models 1 and 3. These assumptions imply that the data
$\\{X_{i},y_{i}\\}_{i\in[n]}$ will satisfy the following deterministic
conditions with probability $\frac{9}{10}$. For convenience we work with these
deterministic conditions directly in our proofs.

######  Assumption 1 (deterministic regularity for linear regression).

Let $\epsilon$ be sufficiently small, and let $r\in(0,\epsilon^{2})$. Assume
$\\{(X_{i},y_{i})\\}_{i\in[n]}\subset\mathbb{R}^{d}\times\mathbb{R}$ is
_$(\epsilon,r)$ -good for linear regression_ (or $(\epsilon,r)$-good if
context is clear), which means there is a partition $[n]=G\cup B$ with
$|G|\geq(1-\epsilon)n$ which satisfies:

  1. 1.

For any $w\in\Delta^{n}$ with $\left\lVert w_{G}\right\rVert_{1}\geq
1-4\epsilon$,
$\frac{1}{2}\boldsymbol{\Sigma}^{\star}\preceq\textup{Cov}_{w_{G}}(\mathbf{X})\preceq\frac{3}{2}\boldsymbol{\Sigma}^{\star}$.

  2. 2.

There is a constant $C_{\textup{est}}$ such that for all
$\theta\in\mathbb{R}^{d}$, there exists a $G^{\prime}\subseteq G$ satisfying
$|G^{\prime}|\geq(1-r)|G|$ such that for all $\epsilon$-saturated
$w\in\Delta^{n}$, if we let $\tilde{w}:=\frac{w_{G^{\prime}}}{\left\lVert
w_{G^{\prime}}\right\rVert_{1}}$,

$\displaystyle\left\lVert\nabla F_{\tilde{w}}(\theta)-\nabla
F^{\star}(\theta)\right\rVert_{2}\leq
C_{\textup{est}}\sqrt{L\epsilon}\left(\sigma+\left\lVert\theta-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\right),$
(7)
$\displaystyle\left\lVert\textup{Cov}_{\tilde{w}}\left(\left\\{g_{i}(\theta)\right\\}_{i\in
G^{\prime}}\right)\right\rVert_{\textup{op}}\leq
C_{\textup{est}}L\left(\left\lVert\theta-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}^{2}+\sigma^{2}\right).$
(8)

  3. 3.

There is a constant $C_{\textup{ub}}$ such that

$\sum_{i\in G}\frac{1}{|G|}f_{i}(\theta^{\star})=\frac{1}{2|G|}\sum_{i\in
G}\left(\left\langle X_{i},\theta^{\star}\right\rangle-y_{i}\right)^{2}\leq
C_{\textup{ub}}\sigma^{2}.$

We defer the proof of the following claim, which establishes the probabilistic
validity of Assumption 1 (up to adjusting constants in definitions) under the
statistical Models 1 and 3, to Appendix A.

######  Proposition 2.

Let $\alpha\geq 1$ and let $\epsilon>0$ be sufficiently small. Let
$\\{(X_{i},y_{i})\\}_{i\in[n]}\subset\mathbb{R}^{d}\times\mathbb{R}$ be an
$\epsilon$-corrupted set of samples from a distribution $\mathcal{D}_{Xy}$ as
in Model 3. Then, if

$n=O\left(\frac{d\alpha^{2}\log
d}{\epsilon^{4}}+\frac{d^{2}\alpha^{1.5}\log(d/\epsilon)}{\epsilon^{3}}\right)\;,$

the set $\\{(X_{i},y_{i})\\}_{i\in[n]}$ is
$(2\epsilon,\frac{\epsilon^{2}}{\alpha})$-good for linear regression with
probability at least $\frac{9}{10}$.

######  Remark 1.

We remark that the gaurantees of Proposition 2 may be recovered with
$n=O(d\log d/\epsilon^{4}+d^{2}\log d\log\frac{1}{\epsilon})$ samples when the
$X_{i}$ are further assumed to be subgaussian.

We observe that Assumption 1 implies the following useful bound.

######  Lemma 3.

Let $w\in\Delta^{n}$ be $\epsilon$-saturated with respect to bipartition
$[n]=G\cup B$, let $G^{\star}\subseteq G$ be the subset in Assumption 1.2
corresponding to $\theta^{\star}$, and let
$\tilde{w}=\frac{w_{G^{\star}}}{\left\lVert w_{G^{\star}}\right\rVert_{1}}$.
Let
$\theta_{\tilde{w}}=\textup{argmin}_{\theta\in\mathbb{R}^{d}}F_{\tilde{w}}(\theta)$
be the empirical minimizer of $F_{\tilde{w}}$. Then,

$\left\lVert\theta_{\tilde{w}}-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\leq
4C_{\textup{est}}\sigma\sqrt{\kappa\epsilon}.$

###### Proof.

Applying Assumption 1.2 with $\theta=\theta^{\star}$, we have that

$\left\lVert\nabla F_{\tilde{w}}(\theta^{\star})\right\rVert_{2}\leq
C_{\textup{est}}\sqrt{L\epsilon}\sigma.$

However, applying Assumption 1.1 on the weights $w_{G^{\star}}$ implies that
$F_{\tilde{w}}$ is $\frac{1}{2}$-strongly convex in the
$\boldsymbol{\Sigma}^{\star}$ norm (since $w_{G^{\star}}$ removes at most
$\epsilon^{2}$ mass from $w_{G}$) and minimized by $\theta_{\tilde{w}}$, and
hence by using consequences of strong convexity and
$(\boldsymbol{\Sigma}^{\star})^{-1}\preceq\frac{1}{\mu}\mathbf{I}$,

$\left\lVert\nabla
F_{\tilde{w}}(\theta^{\star})\right\rVert_{2}\geq\sqrt{\mu}\left\lVert\nabla
F_{\tilde{w}}(\theta^{\star})\right\rVert_{(\boldsymbol{\Sigma}^{\star})^{-1}}\geq\frac{\sqrt{\mu}}{4}\left\lVert\theta_{\tilde{w}}-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}.$

∎

Finally, in Section 4, when we develop an alternative approach to linear
regression based on the robust gradient descent framework, we require a
slightly weaker set of distributional assumptions and deterministic
implications, which we now state.

######  Model 4 (distributional regularity for linear regression, gradient
descent setting).

Given distributions $\mathcal{D}_{X}$ and $\mathcal{D}_{\delta}$ over
$\mathbb{R}^{d},\mathbb{R}^{d}$ respectively, the distribution
$\mathcal{D}_{Xy}$ over $\mathbb{R}^{d}\times\mathbb{R}$ is sampled as
follows: for an underlying vector $\theta^{\star}\in\mathbb{R}^{d}$,
independently sample $X\sim\mathcal{D}_{X}$ and
$\delta\sim\mathcal{D}_{\delta}$, and set
$y\leftarrow\left\langle\theta^{\star},X\right\rangle+\delta$. Further,
$\mathcal{D}_{X}$ and $\mathcal{D}_{\delta}$ satisfy the following regularity
assumptions.

  1. 1.

For
$\boldsymbol{\Sigma}^{\star}=\operatorname*{\mathbb{E}}_{\mathcal{D}_{X}}[XX^{\top}]$
and $0<\mu<L$, we have $\mu\mathbf{I}\preceq\boldsymbol{\Sigma}^{\star}\preceq
L\mathbf{I}$.

  2. 2.

$\mathcal{D}_{X}$ is $C_{2\to 4}$-hypercontractive for a constant $C_{2\to
4}$.

  3. 3.

$\mathcal{D}_{\delta}$ is a distribution with mean zero and variance
$\leq\sigma^{2}$.

The main difference between Model 3 and Model 4 is that the latter no longer
requires hypercontractive noise. This corresponds to the following
deterministic assumption.

######  Assumption 2 (deterministic regularity for linear regression, gradient
descent setting).

Let $\epsilon$ be sufficiently small. For every fixed
$\theta\in\mathbb{R}^{d}$, there is a partition $[n]=G_{\theta}\cup
B_{\theta}$ with $|G_{\theta}|\geq(1-\epsilon)n$ which satisfies: there is a
constant $C_{\textup{est}}$ such that for
$\tilde{w}:=\frac{w_{G_{\theta}}}{\left\lVert
w_{G_{\theta}}\right\rVert_{1}}$,

$\displaystyle\left\lVert\nabla F_{\tilde{w}}(\theta)-\nabla
F^{\star}(\theta)\right\rVert_{2}\leq
C_{\textup{est}}\sqrt{L\epsilon}\left(\sigma+\left\lVert\theta-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\right),$
(9)
$\displaystyle\left\lVert\textup{Cov}_{\tilde{w}}\left(\left\\{g_{i}(\theta)\right\\}_{i\in
G_{\theta}}\right)\right\rVert_{\textup{op}}\leq
C_{\textup{est}}L\left(\left\lVert\theta-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}^{2}+\sigma^{2}\right).$
(10)

The main difference between Assumption 1 and Assumption 2 is that the latter
provides gradient bounds using a different set $G_{\theta}$ for each $\theta$
(as opposed to the former, which uses the same set $G$ for all $\theta$). The
upshot is that the corresponding required sample complexity is lower.

######  Proposition 3 ([CAT+20], Lemma 5.5).

Let $\epsilon>0$ be sufficiently small. Let
$\\{(X_{i},y_{i})\\}_{i\in[n]}\subset\mathbb{R}^{d}\times\mathbb{R}$ be an
$O(\epsilon)$-corrupted set of samples from a distribution $\mathcal{D}_{Xy}$
as in Model 2. Then if

$n=O\left(\frac{d\log(d/\epsilon)}{\epsilon}\right),$

Assumption 2 holds with probability at least $\frac{9}{10}$.

###  2.4 Regularity assumptions: Lipschitz and smooth stochastic optimization

In Section 4, we develop algorithms for the special case of Model 2 when all
$\gamma_{y_{i}}(v):=\gamma(v,y_{i})$, as viewed as a function of its first
variable, satisfy

$\left|\gamma_{y_{i}}^{\prime}(v)\right|\leq
1,\;0\leq\gamma_{y_{i}}^{\prime\prime}(v)\leq 1\text{ for all
}v\in\mathbb{R},\;i\in[n].$

In other words, all $\gamma_{y_{i}}$ are $1$-smooth and $1$-Lipschitz. A
canonical example is when all $y_{i}=\pm 1$ are positive or negative labels,
and $\gamma$ is the logistic loss function, $\gamma(v,y)=\log(1+\exp(-vy))$.
In this setting, we will focus on approximating the (population) regularized
optimizer,

$\theta^{\star}_{\textup{reg}}:=\textup{argmin}_{\theta\in\mathbb{R}^{d}}\left\\{F^{\star}(\theta)+\frac{\mu}{2}\left\lVert\theta\right\rVert_{2}^{2}\right\\},\text{
where
}F^{\star}(\theta):=\operatorname*{\mathbb{E}}_{f\sim\mathcal{D}_{f}}\left[f(\theta)\right].$
(11)

Here, $\mu\in\mathbb{R}_{\geq 0}$ controls the amount of regularization, and
is used to introduce some amount of strong convexity in the problem. Following
Section 2.2, the distribution $\mathcal{D}_{f}$ over sampled functions is
directly dependent on a dataset distribution, $\mathcal{D}_{Xy}$, through the
relationship in Model 2. Concretely, we make the following regularity
assumptions about the distribution $\mathcal{D}_{Xy}$ and its induced
$\mathcal{D}_{f}$.

######  Model 5 (distributional regularity for smooth GLMs).

$\mathcal{D}_{Xy}$, supported on $\mathbb{R}^{d}\times\mathbb{R}$, its
marginals $\mathcal{D}_{X}$, $\mathcal{D}_{y}$, and its induced function
distribution $\mathcal{D}_{f}$, have the following properties.

  1. 1.

Letting the second moment matrix of $\mathcal{D}_{X}$ be
$\boldsymbol{\Sigma}^{\star}$ and $0<L$, we have
$\boldsymbol{\Sigma}^{\star}\preceq L\mathbf{I}$.

  2. 2.

There is a link function $\gamma:\mathbb{R}^{2}\to\mathbb{R}$, such that for
all $y$ in the support of $\mathcal{D}_{y}$, $\gamma_{y}(v):=\gamma(v,y)$
satisfies $\left|\gamma_{y}^{\prime}(v)\right|\leq 1$ and
$0\leq\gamma^{\prime\prime}_{y}(v)\leq 1$ for all $v\in\mathbb{R}$.

  3. 3.

The distribution of $f\sim\mathcal{D}_{f}$ is generated as follows: for
$(X,y)\sim\mathcal{D}_{Xy}$, $f(\theta)=\gamma(\left\langle
X,\theta\right\rangle,y)$.

In Section 5, we further develop algorithms for the special case of Model 2
when all $\gamma_{y_{i}}(v):=\gamma(v,y_{i})$ satisfy only a first-derivative
bound,

$\left|\gamma^{\prime}_{y_{i}}(v)\right|\leq 1\text{ for all
}v\in\mathbb{R},i\in[n].$

In other words, all $\gamma_{y_{i}}$ are $1$-Lipschitz (but possibly non-
smooth). A canonical example is when all $y_{i}=\pm 1$ are positive or
negative labels, and $\gamma$ is the support vector machine loss function
(hinge loss), $\gamma(v,y)=\max(0,1-vy)$. In this setting, we will focus on
approximating the (population) regularized optimizer of the _Moreau envelope_
,

$\displaystyle\theta^{\star}_{\textup{env}}:=\textup{argmin}_{\theta\in\mathbb{R}^{d}}\left\\{F^{\star}_{\gamma}\left(\theta\right)+\frac{\mu}{2}\left\lVert\theta\right\rVert_{2}^{2}\right\\},\text{
where
}F^{\star}(\theta)=\operatorname*{\mathbb{E}}_{f\sim\mathcal{D}_{f}}[f(\theta)],$
$\displaystyle\text{and
}F^{\star}_{\lambda}(\theta):=\inf_{\theta^{\prime}}\left\\{F^{\star}(\theta^{\prime})+\frac{1}{2\lambda}\left\lVert\theta-\theta^{\prime}\right\rVert_{2}^{2}\right\\}.$

The Moreau envelope is extremely well-studied [PB14], and can be viewed as a
smooth approximation to a non-smooth function. We choose to focus on
optimizing the Moreau envelope because it is amenable to acceleration
techniques, trading off approximation for smoothness.

######  Model 6 (distributional regularity for Lipschitz GLMs).

$\mathcal{D}_{Xy}$, supported on $\mathbb{R}^{d}\times\mathbb{R}$, its
marginals $\mathcal{D}_{X}$, $\mathcal{D}_{y}$, and its induced function
distribution $\mathcal{D}_{f}$, have the following properties.

  1. 1.

Letting the second moment matrix of $\mathcal{D}_{X}$ be
$\boldsymbol{\Sigma}^{\star}$ and $0<L$, we have
$\boldsymbol{\Sigma}^{\star}\preceq L\mathbf{I}$.

  2. 2.

There is a link function $\gamma:\mathbb{R}^{2}\to\mathbb{R}$, such that for
all $y$ in the support of $\mathcal{D}_{y}$, $\gamma_{y}(v):=\gamma(v,y)$
satisfies $\left|\gamma_{y}^{\prime}(v)\right|\leq 1$ for all
$v\in\mathbb{R}$.

  3. 3.

The distribution of $f\sim\mathcal{D}_{f}$ is generated as follows: for
$(X,y)\sim\mathcal{D}_{Xy}$, $f(\theta)=\gamma(\left\langle
X,\theta\right\rangle,y)$.

In other words, Model 6 is Model 5 without the smoothness assumption. Under
the weaker Model 6, [DKK+19] showed that we can make the following simplifying
deterministic assumptions about our observed dataset (which also extends to
Model 5, as it a superset of conditions).

######  Assumption 3 (deterministic regularity for Lipschitz regression).

The set $\\{(X_{i},y_{i})\\}_{i\in[n]}\subset\mathbb{R}^{d}\times\mathbb{R}$,
and the link function $\gamma:\mathbb{R}^{2}\to\mathbb{R}$, have the following
properties, for $[n]=G\cup B$.

  1. 1.

Letting $w^{\star}_{G}:=\frac{1}{|G|}\mathbbm{1}_{G}$,
$\textup{Cov}_{w^{\star}_{G}}(\mathbf{X})\preceq\frac{3}{2}L\mathbf{I}$.

  2. 2.

There is a constant $C_{\textup{est}}$ such that for all
$\theta\in\mathbb{R}^{d}$ and all saturated $w\in\Delta^{n}$, letting
$\tilde{w}:=\frac{w_{G}}{\left\lVert w_{G}\right\rVert_{1}}$,

$\displaystyle\left\lVert\nabla F_{\tilde{w}}(\theta)-\nabla
F^{\star}(\theta)\right\rVert_{2}\leq C_{\textup{est}}\sqrt{L\epsilon},$
$\displaystyle\left\lVert\textup{Cov}_{\tilde{w}}\left(\left\\{g_{i}(\theta)\right\\}_{i\in
G}\right)\right\rVert_{\textup{op}}\leq C_{\textup{est}}L.$

######  Proposition 4 ([DKK+19], Proposition C.3).

Let $\epsilon>0$ be sufficiently small. Let
$\\{(X_{i},y_{i})\\}_{i\in[n]}\subset\mathbb{R}^{d}\times\mathbb{R}$ be an
$\epsilon$-corrupted set of samples from a distribution $\mathcal{D}_{Xy}$ as
in Model 6. Then, if

$n=O\left(\frac{d\log\left(d/\epsilon\right)}{\epsilon}\right),$

Assumption 3 holds with probability at least $\frac{9}{10}$.

Finally, we make the useful observation that $F^{\star}$ under Model 6 is also
Lipschitz.

######  Lemma 4.

Under Model 6,
$F^{\star}=\operatorname*{\mathbb{E}}_{f\sim\mathcal{D}_{f}}[f]$ is
$\sqrt{L}$-Lipschitz.

###### Proof.

We wish to prove $\left\lVert\nabla
F^{\star}(\theta)\right\rVert_{2}\leq\sqrt{L}$ for all
$\theta\in\mathbb{R}^{d}$. By nonnegativity of covariance,

$\displaystyle\operatorname*{\mathbb{E}}_{f\sim\mathcal{D}_{f}}\left[\left(\nabla
f(\theta)\right)\left(\nabla
f(\theta)\right)^{\top}\right]\succeq\left(\operatorname*{\mathbb{E}}_{f\sim\mathcal{D}_{f}}\left[\nabla
f(\theta)\right]\right)\left(\operatorname*{\mathbb{E}}_{f\sim\mathcal{D}_{f}}\left[\nabla
f(\theta)\right]\right)^{\top}=\left(\nabla
F^{\star}(\theta)\right)\left(\nabla F^{\star}(\theta)\right)^{\top}.$

Since the right-hand side of the above display is rank-one, $\left\lVert\nabla
F^{\star}(\theta)\right\rVert_{2}^{2}$ is at most the largest eigenvalue of
the gradient second moment matrix, so it suffices to show the left-hand side
is $\preceq L\mathbf{I}$: assuming $f\sim\mathcal{D}_{f}$ is associated with
$(X,y)\sim\mathcal{D}_{Xy}$,

$\displaystyle\operatorname*{\mathbb{E}}_{f\sim\mathcal{D}_{f}}\left[\left(\nabla
f(\theta)\right)\left(\nabla
f(\theta)\right)^{\top}\right]=\operatorname*{\mathbb{E}}_{X,y\sim\mathcal{D}_{Xy}}\left[X\left(\gamma^{\prime}\left(\left\langle
X,\theta\right\rangle,y\right)\right)^{2}X^{\top}\right]\preceq\boldsymbol{\Sigma}^{\star}\preceq
L\mathbf{I}.$

∎

###  2.5 Robustly decreasing the covariance operator norm

In this section, we describe a procedure, $\mathsf{FastCovFilter}$, which
takes as input a set of vectors
$\mathbf{V}:=\\{v_{i}\\}_{i\in[n]}\in\mathbb{R}^{n\times d}$ such that for an
(unknown) bipartition $[n]=G\cup B$, it is promised that
$\operatorname*{\mathbb{E}}_{i\sim_{\textup{unif}}G}\left[v_{i}v_{i}^{\top}\right]$
is bounded in operator norm by $R$. Given this promise,
$\mathsf{FastCovFilter}$ takes as input a set of saturated weights
$w\in\Delta^{n}$ and performs a sequence of safe weight removals to obtain a
new saturated $w^{\prime}\in\Delta^{n}$, such that
$\sum_{i\in[n]}w^{\prime}_{i}v_{i}v_{i}^{\top}$ is bounded in operator norm by
$O(R)$. Moreover, the procedure runs in nearly-linear time in the description
size of $\mathbf{V}$. We formally describe the guarantees of
$\mathsf{FastCovFilter}$ here as Proposition 5, and defer the proof to
Appendix A.

######  Proposition 5.

There is an algorithm, $\mathsf{FastCovFilter}$ (Algorithm 7), taking inputs
$\mathbf{V}:=\\{v_{i}\\}_{i\in[n]}\in\mathbb{R}^{n\times d}$, saturated
weights $w\in\Delta^{n}$ with respect to bipartition $[n]=G\cup B$ with
$|B|=\epsilon n$, $\delta\in(0,1]$, and $R\geq 0$ with the promise that

$\left\lVert\sum_{i\in
G}\frac{1}{|G|}v_{i}v_{i}^{\top}\right\rVert_{\textup{op}}\leq R.$

Then, with probability at least $1-\delta$, $\mathsf{FastCovFilter}$ returns
saturated $w^{\prime}\in\Delta^{n}$ such that

$\left\lVert\sum_{i\in[n]}w^{\prime}_{i}v_{i}v_{i}^{\top}\right\rVert_{\textup{op}}\leq
5R.$

The runtime of $\mathsf{FastCovFilter}$ is

$O\left(nd\log^{3}(n)\log\left(\frac{n}{\delta}\right)\right).$

An algorithm with similar guarantees to $\mathsf{FastCovFilter}$ is implicit
in the work [DHL19], but we give a self-contained exposition in this work for
completeness. Our approach in designing $\mathsf{FastCovFilter}$ is to use a
matrix multiplicative weights-based potential function, along with the
filtering approach suggested by Lemma 1, to rapidly decrease the quadratic
form of the empirical second moment matrix in a number of carefully-chosen
directions, and argue this quickly decreases the potential. We remark that
this potential-based approach was also used in the recent work [DKK+21].

##  3 Linear regression

Throughout this section, we operate under Models 1 and 3, and Assumption 1.
Namely, there is a dataset $\\{X_{i},y_{i}\\}_{i\in[n]}$ and an unknown
bipartition $[n]=G\cup B$, such that $\\{X_{i},y_{i}\\}_{i\in G}$ were draws
from a distribution $\mathcal{D}_{Xy}$, and we wish to estimate

$\theta^{\star}:=\textup{argmin}_{\theta\in\mathbb{R}^{d}}F^{\star}(\theta),\text{
where
}F^{\star}(\theta):=\operatorname*{\mathbb{E}}_{X,y\sim\mathcal{D}_{Xy}}\left[\frac{1}{2}\left(\left\langle
X,\theta\right\rangle-y\right)^{2}\right].$

We follow notation (4), (6) in this section, and define $G^{\star}\subseteq G$
as the subset given by Assumption 1.2 for the true minimizer $\theta^{\star}$.
We also define $B^{\star}:=[n]\setminus G^{\star}$, so $B^{\star}\supseteq B$.

We begin in Section 3.1 with a preliminary on filtering under a weaker
assumption than the safety condition in Definition 2; in particular,
Assumption 1 is not quite compatible with Definition 2 because $G^{\prime}$
can change based on $\theta$, but will not affect saturation by more than
constants. In Section 3.2, we then state a general “identifiability proof”
showing that controlling certain quantities such as the operator norm of the
gradient covariances and near-optimality of a current estimate $\theta$, with
respect to some weights $w\in\Delta^{n}$, suffices to bound closeness of
$\theta$ and $\theta^{\star}$. This identifiability proof is motivated by an
analogous argument in [BP21], and will guide our algorithm design. Next, we
give a self-contained oracle which rapidly halves the distance to
$\theta^{\star}$ outside of a sufficiently large radius in Section 3.3, and
analyze the final phase (once this radius is reached) in Section 3.4. We put
the pieces together to give our main result and full algorithm in Section 3.5.

###  3.1 Filtering under $(\epsilon,\frac{\epsilon^{2}}{\alpha})$-goodness

Throughout this section, we will globally fix a value (for a sufficiently
large constant)

$\alpha=O\left(\max\left(1,\epsilon\kappa\log\frac{R_{0}}{\sigma}\right)\right).$
(12)

In this definition, $R_{0}$ is an initial distance bound on
$\left\lVert\theta_{0}-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}$
we will provide to our final algorithm (see the statement of Theorem 5). We
operate under Assumption 1 such that our dataset is
$(\epsilon,\frac{\epsilon^{2}}{\alpha})$-good, which inflates the sample
complexity of Proposition 2 by a factor depending on $\alpha$.

At a high level, this technical complication is to ensure that throughout the
algorithm we never remove more than $3\epsilon$ mass from $G$, and that our
weights are always $\epsilon$-saturated with respect to $G\cup B$, allowing
for inductive application of Assumption 1. Formally, we demonstrate the
following (simple) generalization of Lemma 1, using safety definitions on
different sets.

######  Lemma 5.

Suppose Assumption 1 holds with $r=\frac{\epsilon^{2}}{\alpha}$. Consider any
algorithm which performs the following weight removal.

  1. 1.

$w_{0}\leftarrow\frac{1}{n}\mathbbm{1}$

  2. 2.

For $0\leq t<T$:

    1. (a)

$[w_{t+1}]_{i}\leftarrow\left(1-\frac{\tau_{i}}{\tau_{\max}}\right)[w_{t}]_{i}$
for all $i\in[n]$ and $\tau_{\max}:=\max_{i\in[n]\mid w_{i}\neq 0}\tau_{i}$,
for $\\{\tau_{i}\\}_{i\in[n]}$ safe with respect to $w_{t}$ and a bipartition
$G^{\prime}_{t}$, $[n]\setminus G^{\prime}_{t}$ where $G^{\prime}_{t}\subseteq
G$, $|G^{\prime}_{t}|\geq(1-\frac{\epsilon^{2}}{\alpha})|G|$.

Suppose the number of distinct sets $G^{\prime}_{t}$ throughout the algorithm
is bounded by $\frac{\alpha}{2\epsilon}$. Then throughout the algorithm,
$w_{t}$ is $\epsilon$-saturated with respect to the bipartition $G\cup B$.

###### Proof.

Consider some distinct set $G^{\prime}$. The proof of Lemma 1 demonstrates
that under these assumptions, in every iteration $t$ using $G^{\prime}$ for
weight removal, the amount of mass removed from $G^{\prime}$ is less than the
amount removed from $[n]\setminus G^{\prime}$. Moreover, the amount of mass
removed from $G$ in these iterations can be at most the amount removed from
$G^{\prime}$, plus the weight assigned to the _entire difference_ $G\setminus
G^{\prime}$, which is at most a $\frac{\epsilon^{2}}{\alpha}$ fraction.
Similarly, the amount of mass removed from $B$ is at least the amount of mass
removed from $[n]\setminus G^{\prime}$, minus their set difference, which is
again at most $\frac{\epsilon^{2}}{\alpha}$. Combining over all distinct sets,
this deviation is at most $\epsilon$. ∎

It will be straightforward to verify that throughout this section, all weight
removals we perform will be of the form in Lemma 5, and that we never perform
weight removals with respect to more than $\frac{\alpha}{2\epsilon}$ distinct
sets. Hence, we will always assume that any weights $w$ we discuss are
$\epsilon$-saturated with respect to the bipartition $G\cup B$, and thus has
$\left\lVert w_{G}\right\rVert_{1}\geq 1-3\epsilon$, allowing for application
of Assumption 1. Finally, we remark we sometimes will apply Assumption 1.1
with weight vectors $w_{G^{\star}}$ instead of $w_{G}$; since their difference
is $O(\epsilon^{2})<\epsilon$, this is a correct application for any
$\left\lVert w_{G}\right\rVert_{1}\geq 1-3\epsilon$.

###  3.2 Identifiability proof for linear regression

In this section, we give an identifiability proof similar to that appearing in
[BP21] which demonstrates, for a given weight-parameter estimate pair
$(w,\theta)\in\Delta^{n}\times\mathbb{R}^{d}$, verifiable technical conditions
on this pair which certify a bound on
$\left\lVert\theta-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}$.
In short, Proposition 6 will show that if both of the quantities

$\left\lVert\nabla
F_{w}(\theta)\right\rVert_{(\boldsymbol{\Sigma}^{\star})^{-1}},\;\left\lVert\textup{Cov}_{w}\left(\left\\{g_{i}(\theta)\right\\}_{i\in[n]}\right)\right\rVert_{\textup{op}}$
(13)

are simultaneously controlled, then we obtain a distance bound to
$\theta_{G^{\star}}$.

######  Proposition 6.

Let $w\in\Delta^{n}$ be $\epsilon$-saturated with respect to the bipartition
$[n]=G\cup B$, and let $\theta\in\mathbb{R}^{d}$. Assuming $\epsilon\kappa$ is
sufficiently small, there is a universal constant $C_{\textup{id}}$ such that

$\left\lVert\theta-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\leq
C_{\textup{id}}\left(\sqrt{\epsilon}\left(\sigma\sqrt{\kappa}+\sqrt{\frac{\left\lVert\textup{Cov}_{w}\left(\left\\{g_{i}(\theta)\right\\}_{i\in[n]}\right)\right\rVert_{\textup{op}}}{\mu}}\right)+\left\lVert\nabla
F_{w}(\theta)\right\rVert_{(\boldsymbol{\Sigma}^{\star})^{-1}}\right).$

###### Proof.

Let $G^{\prime}$ be the set promised by Assumption 1.2 for the point $\theta$.
Throughout this proof, we define $G_{\theta}:=G^{\prime}\cap G^{\star}$, and
we let $w^{\star}_{\theta}:=\frac{1}{|G_{\theta}|}\mathbbm{1}_{G_{\theta}}$ be
the uniform weights on $G_{\theta}$. Finally, let $\hat{\theta}$ minimize
$F_{w^{\star}_{\theta}}$. By applying Lemma 3 on the weights
$w^{\star}_{\theta}$, we have that
$\|\hat{\theta}-\theta^{\star}\|_{\boldsymbol{\Sigma}^{\star}}=O(\sigma\sqrt{\kappa\epsilon})$.
Thus, in the remainder of the proof we focus on bounding
$\|\theta-\hat{\theta}\|_{\boldsymbol{\Sigma}^{\star}}$ by the required
quantity.

Let $\mathcal{C}(w^{\star}_{\theta},\tilde{w})$ supported on $[n]\times[n]$ be
an optimal _coupling_ between $i\sim w^{\star}_{\theta}$ and $j\sim\tilde{w}$,
where $\tilde{w}:=\frac{w}{\left\lVert w\right\rVert_{1}}$ is the distribution
proportional to $w$; we denote the coupling as $\mathcal{C}$ for short. For a
pair $(i,j)\in[n]\times[n]$ sampled from $\mathcal{C}$, we let
$\mathbbm{1}_{i=j}$ be the indicator of the event $i=j$, and similarly define
$\mathbbm{1}_{i\neq j}$. From the total variation characterization of
coupling, we have by Lemma 2 that

$\operatorname*{\mathbb{E}}_{i,j\sim\mathcal{C}}\left[\mathbbm{1}_{i\neq
j}\right]=\left\lVert\tilde{w}-w^{\star}_{\theta}\right\rVert_{1}\leq
9\epsilon.$

Here we used that $\left\lVert
w^{\star}_{\theta}-w^{\star}_{G}\right\rVert_{1}=O(\epsilon^{2})$, where
$w^{\star}_{G}$ is the uniform distribution on $G$, as in Lemma 2. Now, let
$v=\theta-\hat{\theta}$. We have by Assumption 1.1 that

$\operatorname*{\mathbb{E}}_{i\sim w^{\star}_{\theta}}\left[\left\langle
v,X_{i}\right\rangle\left\langle
X_{i},\theta-\hat{\theta}\right\rangle\right]=\left\langle\theta-\hat{\theta},\operatorname*{\mathbb{E}}_{i\sim
w^{\star}_{\theta}}\left[X_{i}X_{i}^{\top}\right](\theta-\hat{\theta})\right\rangle\geq\frac{1}{2}\left\lVert\theta-\hat{\theta}\right\rVert_{\boldsymbol{\Sigma}^{\star}}^{2}.$
(14)

On the other hand,

$\displaystyle\operatorname*{\mathbb{E}}_{i\sim
w^{\star}_{\theta}}\left[\left\langle v,X_{i}\right\rangle\left\langle
X_{i},\theta-\hat{\theta}\right\rangle\right]$
$\displaystyle=\operatorname*{\mathbb{E}}_{i\sim
w^{\star}_{\theta}}\left[\left\langle v,X_{i}\right\rangle\left(\left\langle
X_{i},\theta\right\rangle-
y_{i}\right)\right]+\operatorname*{\mathbb{E}}_{i\sim
w^{\star}_{\theta}}\left[\left\langle
v,X_{i}\right\rangle\left(y_{i}-\left\langle
X_{i},\hat{\theta}\right\rangle\right)\right]$
$\displaystyle=\operatorname*{\mathbb{E}}_{i\sim
w^{\star}_{\theta}}\left[\left\langle
v,g_{i}(\theta)\right\rangle\right]-\operatorname*{\mathbb{E}}_{i\sim
w^{\star}_{\theta}}\left[\left\langle
v,g_{i}(\hat{\theta})\right\rangle\right]=\operatorname*{\mathbb{E}}_{i\sim
w^{\star}_{\theta}}\left[\left\langle v,g_{i}(\theta)\right\rangle\right].$

Here, we used that $\operatorname*{\mathbb{E}}_{i\sim
w^{\star}_{\theta}}\left[g_{i}(\hat{\theta})\right]=\nabla
F_{w^{\star}_{\theta}}(\hat{\theta})=0$ by definition of $\hat{\theta}$.
Continuing,

$\displaystyle\operatorname*{\mathbb{E}}_{i\sim
w^{\star}_{\theta}}\left[\left\langle v,X_{i}\right\rangle\left\langle
X_{i},\theta-\hat{\theta}\right\rangle\right]$
$\displaystyle=\operatorname*{\mathbb{E}}_{i,j\sim\mathcal{C}}\left[\left\langle
v,g_{i}(\theta)\right\rangle\mathbbm{1}_{i=j}\right]+\operatorname*{\mathbb{E}}_{i,j\sim\mathcal{C}}\left[\left\langle
v,g_{i}(\theta)\right\rangle\mathbbm{1}_{i\neq j}\right]$ (15)
$\displaystyle=\operatorname*{\mathbb{E}}_{j\sim\tilde{w}}\left[\left\langle
v,g_{j}(\theta)\right\rangle\right]-\operatorname*{\mathbb{E}}_{i,j\sim\mathcal{C}}\left[\left\langle
v,g_{j}(\theta)\right\rangle\mathbbm{1}_{i\neq
j}\right]+\operatorname*{\mathbb{E}}_{i,j\sim\mathcal{C}}\left[\left\langle
v,g_{i}(\theta)\right\rangle\mathbbm{1}_{i\neq j}\right]$
$\displaystyle\leq\left\langle v,\nabla
F_{\tilde{w}}\left(\theta\right)\right\rangle+3\sqrt{\epsilon}\left(\operatorname*{\mathbb{E}}_{j\sim\tilde{w}}\left[\left\langle
v,g_{j}(\theta)\right\rangle^{2}\right]^{\frac{1}{2}}+\operatorname*{\mathbb{E}}_{i\sim
w^{\star}_{\theta}}\left[\left\langle
v,g_{i}(\theta)\right\rangle^{2}\right]^{\frac{1}{2}}\right).$

In the last line, we used Cauchy-Schwarz and
$\operatorname*{\mathbb{E}}_{i,j\sim\mathcal{C}}[\mathbbm{1}_{i\neq
j}^{2}]\leq 9\epsilon$ to deal with the second and third terms. To bound the
term corresponding to $i\sim w^{\star}_{\theta}$,

$\displaystyle\operatorname*{\mathbb{E}}_{i\sim
w^{\star}_{\theta}}\left[\left\langle
v,g_{i}(\theta)\right\rangle^{2}\right]\leq\left\lVert
v\right\rVert_{2}^{2}\left\lVert\textup{Cov}_{w^{\star}_{\theta}}\left(\left\\{g_{i}(\theta)\right\\}_{i\in[n]}\right)\right\rVert_{\textup{op}}\leq
1.1C_{\textup{est}}\kappa\left\lVert
v\right\rVert_{\boldsymbol{\Sigma}^{\star}}^{2}\left(\left\lVert\theta-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}^{2}+\sigma^{2}\right).$

The first inequality is by definition of
$\left\lVert\cdot\right\rVert_{\textup{op}}$, and the second used $\left\lVert
v\right\rVert_{2}^{2}\geq\frac{1}{\mu}\left\lVert
v\right\rVert_{\boldsymbol{\Sigma}^{\star}}^{2}$ and Assumption 1.2, since
$G_{\theta}$ is a subset of the set $G^{\prime}$ promised by Assumption 1.2
(where we adjusted by a constant for normalization). We similarly arrive at
the bound

$\operatorname*{\mathbb{E}}_{j\sim\tilde{w}}\left[\left\langle
v,g_{j}(\theta)\right\rangle^{2}\right]\leq\frac{\left\lVert
v\right\rVert_{\boldsymbol{\Sigma}^{\star}}^{2}}{\mu}\left\lVert\textup{Cov}_{w}\left(\left\\{g_{i}(\theta)\right\\}_{i\in[n]}\right)\right\rVert_{\textup{op}}.$

Now plugging the above displays into (15) and combining with (14), as well as
using
$\left\lVert\theta-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}^{2}+\sigma^{2}=O(\|\theta-\hat{\theta}\|_{\boldsymbol{\Sigma}^{\star}}^{2}+\sigma^{2})$
by the triangle inequality and our earlier bound on
$\|\theta^{\star}-\hat{\theta}\|_{\boldsymbol{\Sigma}^{\star}}$,

$\displaystyle\left\lVert\theta-\hat{\theta}\right\rVert_{\boldsymbol{\Sigma}^{\star}}^{2}$
$\displaystyle\leq\left\lVert\theta-\hat{\theta}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\left\lVert\nabla
F_{\tilde{w}}(\theta)\right\rVert_{(\boldsymbol{\Sigma}^{\star})^{-1}}+O\left(\sqrt{\epsilon\kappa}\right)\left(\sigma\left\lVert\theta-\hat{\theta}\right\rVert_{\boldsymbol{\Sigma}^{\star}}+\left\lVert\theta-\hat{\theta}\right\rVert_{\boldsymbol{\Sigma}^{\star}}^{2}\right)$
$\displaystyle+O\left(\sqrt{\epsilon}\right)\left(\left\lVert\theta-\hat{\theta}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\sqrt{\frac{\left\lVert\textup{Cov}_{w}\left(\left\\{g_{i}(\theta)\right\\}_{i\in[n]}\right)\right\rVert_{\textup{op}}}{\mu}}\right).$

In the first line, we used Cauchy-Schwarz to bound the term $\left\langle
v,\nabla F_{\tilde{w}}(\theta)\right\rangle$. Dividing through by
$\|\theta-\hat{\theta}\|_{\boldsymbol{\Sigma}^{\star}}$, and using that
$\epsilon\kappa$ is sufficiently small, yields the conclusion. ∎

Proposition 6 suggests a natural approach. On the one hand, if $\theta$ is an
(approximate) minimizer to $F_{w}$, the first quantity in (13) will be small.
On the other hand, for a fixed $\theta\in\mathbb{R}^{d}$, we can filter on $w$
using our subroutine $\mathsf{FastCovFilter}$ until the second quantity in
(13) is small. The main challenge is accomplishing both bounds simultaneously.
To this end, we show that the number of times we have to repeat this process
of filtering and then computing an approximate minimizer is bounded, using
$F_{w}$ as a potential function; by preprocessing $F_{w}$ so that is smooth at
all points, if $\theta$ is an approximate minimizer for the $F_{w}$ attained
_after_ filtering on gradients at $\theta$, then we can exit the subroutine.
Otherwise, we argue we make substantial function progress by calling an
empirical risk minimizer to terminate quickly. We make this strategy formal in
the following sections.

###  3.3 Halving the distance to $\theta^{\star}$

In this section, we design a procedure, $\mathsf{HalfRadiusLinReg}$, with the
following guarantee. Suppose that we have $\epsilon$-saturated
$\bar{w}\in\Delta^{n}$ and a parameter $\bar{\theta}\in\mathbb{R}^{d}$, as
well as scalar $R$ with the promise

$\left\lVert\bar{\theta}-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\leq
R.$

The goal of $\mathsf{HalfRadiusLinReg}$ is to return a new $\theta$ with
$\left\lVert\theta-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\leq\frac{1}{2}R$;
we require that $R=\Omega(\sigma)$ for a sufficiently large constant in this
section. We do so by using saturated weights to guide a potential analysis,
crucially using Proposition 6. Before stating $\mathsf{HalfRadiusLinReg}$, we
require two additional helper tools. The first is an approximate optimization
procedure.

######  Definition 3.

We call $\mathcal{O}_{\textup{ERM}}$ a $\gamma$-approximate ERM oracle if on
input $F:\mathbb{R}^{d}\rightarrow\mathbb{R}$ it returns a point
$\hat{\theta}$ such that $F(\hat{\theta})-F(\theta^{\star}_{F})\leq\gamma$,
for $\theta^{\star}_{F}:=\textup{argmin}_{\theta\in\mathbb{R}^{d}}F(\theta)$.

The second controls the initial error, which we use to yield distance bounds
via strong convexity.

######  Lemma 6.

There is an algorithm, $\mathsf{FunctionFilter}$, which takes as input
$\epsilon$-saturated $w\in\Delta^{n}$, $\theta\in\mathbb{R}^{d}$, and
$R\geq\left\lVert\theta-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}$,
and produces $\epsilon$-saturated $w^{\prime}\in\Delta^{n}$ such that
$F_{w^{\prime}}(\theta)\leq 2C_{\textup{ub}}(\sigma^{2}+R^{2})$, in time
$O(nd+n\log\frac{D}{R^{2}})$, where $D$ is a bound on the largest
$\frac{1}{2}(\left\langle X_{i},\theta\right\rangle-y_{i})^{2}$ for any
nonzero $w_{i}$.

###### Proof.

Define $\tau_{i}:=f_{i}(\theta)=\frac{1}{2}(\left\langle
X_{i},\theta\right\rangle-y_{i})^{2}$. Assumption 1.1 shows that
$F_{w_{G^{\star}}}$ is $\frac{3}{2}$-smooth in the
$\boldsymbol{\Sigma}^{\star}$ norm, since its Hessian is exactly $\left\lVert
w_{G^{\star}}\right\rVert_{1}\textup{Cov}_{w_{G^{\star}}}(\mathbf{X})\preceq\frac{3}{2}\boldsymbol{\Sigma}^{\star}$.
Moreover, letting $\hat{\theta}$ minimize $F_{w_{G^{\star}}}$, by Lemma 3, the
triangle inequality and the assumed bound $R$,

$\left\lVert\theta-\hat{\theta}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\leq
R+\left\lVert\theta^{\star}-\hat{\theta}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\leq
R+\sigma.$

Hence, assuming $C_{\textup{ub}}$ is large enough and $R\geq\sigma$,
smoothness demonstrates

$\sum_{i\in G^{\star}}w_{i}\tau_{i}=F_{w_{G^{\star}}}(\theta)\leq
F_{w_{G^{\star}}}(\hat{\theta})+\frac{3}{4}\left\lVert\theta-\hat{\theta}\right\rVert_{\boldsymbol{\Sigma}^{\star}}^{2}\leq
C_{\textup{ub}}(\sigma^{2}+R^{2}),$

Next, letting $\tau_{\max}:=\max_{i\in[n]\mid w_{i}\neq 0}\tau_{i}$ and
$K\in\mathbb{N}\cup\\{0\\}$ be smallest such that

$\sum_{i\in[n]}\left(1-\frac{\tau_{i}}{\tau_{\max}}\right)^{K}w_{i}\tau_{i}\leq
2C_{\textup{ub}}(\sigma^{2}+R^{2}),$

each of the first $K$ weight removals according to the scores $\tau$ are safe
with respect to $G^{\star}\cup B^{\star}$, and Lemma 5 implies we can output
$\epsilon$-saturated
$w^{\prime}_{i}=\left(1-\frac{\tau_{i}}{\tau_{\max}}\right)^{K}w_{i}$
entrywise. It remains to binary search for $K$; given access to the scores
$\tau$, checking if a given $K$ passes the above display takes $O(n)$ time. We
can upper bound $K$ by the following inequality:

$\sum_{i\in[n]}\left(1-\frac{\tau_{i}}{\tau_{\max}}\right)^{K}w_{i}\tau_{i}\leq\sum_{i\in[n]}\exp\left(-\frac{K\tau_{i}}{\tau_{\max}}\right)w_{i}\tau_{i}\leq\frac{1}{eK}\sum_{i\in[n]}w_{i}\tau_{\max}\leq\frac{\tau_{\max}}{K}.$

Here the second inequality used $x\exp(-Cx)\leq\frac{1}{eC}$ for all
nonnegative $x$, $C$, where we chose $C=\frac{K}{\tau_{\max}}$ and
$x=\tau_{i}$. Hence, $K=O(\frac{D}{R^{2}})$. The runtime follows as computing
scores takes time $O(nd)$. ∎

We remark that every time we use $\mathsf{FunctionFilter}$ is a weight removal
of the form in Lemma 5, which is safe with respect to the bipartition
$G^{\star}\cup B^{\star}$. This accounts for one distinct set throughout.

Algorithm 1
$\mathsf{HalfRadiusLinReg}(\mathbf{X},y,\bar{w},\bar{\theta},R,D,\mathcal{O}_{\textup{ERM}},\delta)$

1: Input: Dataset
$\mathbf{X}=\left\\{X_{i}\right\\}_{i\in[n]}\in\mathbb{R}^{n\times d}$,
$y=\\{y_{i}\\}_{i\in[n]}\in\mathbb{R}^{n}$, satisfying Assumption 1,
$\epsilon$-saturated $\bar{w}\in\Delta^{n}$ with respect to bipartition
$[n]=G\cup B$ such that
$\mathbf{X}^{\top}\textbf{{diag}}\left(\bar{w}\right)\mathbf{X}\preceq
8L\mathbf{I}$, $\bar{\theta}\in\mathbb{R}^{d}$ with
$\left\lVert\bar{\theta}-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\leq
R$, $\delta\in(0,1)$, $O(\frac{\sigma^{2}}{\kappa})$-approximate ERM oracle
$\mathcal{O}_{\textup{ERM}}$,  $\displaystyle D$
$\displaystyle\geq\max_{i\in[n]}\frac{1}{2}(\left\langle
X_{i},\bar{\theta}\right\rangle-y_{i})^{2}.$ (16)

2: Output: With probability $\geq 1-\delta$, saturated $w$ with respect to
$G\cup B$, and $\theta$ with
$\left\lVert\theta-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\leq\frac{1}{2}R.$

3: $t\leftarrow 0$,
$w^{(0)}\leftarrow\mathsf{FunctionFilter}(w,\bar{\theta},R,D)$,
$\Delta_{0}\leftarrow\infty$

4: while $\Delta_{t}>\frac{R^{2}}{512\kappa C_{\textup{id}}^{2}}$ do

5: $\theta^{(t)}\leftarrow\mathcal{O}_{\textup{ERM}}(F_{w^{(t)}})$

6:
$w^{(t+1)}\leftarrow\mathsf{FastCovFilter}\left(\left\\{g_{i}\left(\theta^{(t)}\right)\right\\}_{i\in[n]},w^{(t)},\frac{\delta}{O(\kappa)},40C_{\textup{est}}C_{\textup{ub}}LR^{2}\right)$

7: $\Delta_{t+1}\leftarrow
F_{w^{(t+1)}}\left(\theta^{(t)}\right)-F_{w^{(t+1)}}\left(\theta^{(t+1)}\right)$

8: $t\leftarrow t+1$

9: end while

10: return $\left(w^{(t)},\theta^{(t-1)}\right)$

######  Lemma 7.

$\mathsf{HalfRadiusLinReg}$ is correct, i.e. if its preconditions are met, it
successfully returns $(w,\theta)$ such that $w\in\Delta^{n}$ is
$\epsilon$-saturated and
$\left\lVert\theta-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\leq\frac{1}{2}R$.
It runs in $O(\kappa)$ calls to $\mathcal{O}_{\textup{ERM}}$, plus

$O\left(n\log\left(\frac{D}{R^{2}}\right)+nd\kappa\log^{3}(n)\log\left(\frac{n\kappa}{\delta}\right)\right)\textup{
additional time.}$

###### Proof.

We discuss correctness and runtime separately.

Correctness. The first step of our correctness proof is to show that
throughout the algorithm,

$\left\lVert\theta^{(t)}-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\leq
6\sqrt{C_{\textup{ub}}}R.$ (17)

To see this, consider iteration $t$ and suppose $w^{(t)}$ is
$\epsilon$-saturated. Let $G^{\prime}\subseteq G$ be the set promised by
Assumption 1.2 for the pair $(w^{(t)},\theta^{(t)})$, and let
$G_{t}=G^{\prime}\cap G^{\star}$. At the beginning of the algorithm we applied
$\mathsf{FunctionFilter}$ (see Lemma 6), and the minimum value of $F_{w}$ is
monotone nonincreasing as $w$ is decreasing, and $\mathcal{O}_{\textup{ERM}}$
only decreases function value (else Line 4 would fail), so since
$R\geq\sigma$,

$F_{w_{G_{t}}^{(t)}}\left(\theta^{(t)}\right)\leq
F_{w^{(t)}}\left(\theta^{(t)}\right)\leq 4C_{\textup{ub}}R^{2}.$

On the other hand, $F_{w_{G_{t}}^{(t)}}$ is $\frac{1}{3}$-strongly convex in
the $\boldsymbol{\Sigma}^{\star}$ norm by Assumption 1.1 (adjusting for the
normalization factor), and hence letting $\theta^{\star}_{t}$ be the minimizer
of $F_{w_{G_{t}}^{(t)}}$, we have

$4C_{\textup{ub}}R^{2}\geq F_{w_{G_{t}}^{(t)}}\left(\theta^{(t)}\right)\geq
F_{w_{G_{t}}^{(t)}}\left(\theta^{(t)}\right)-F_{w_{G_{t}}^{(t)}}\left(\theta^{\star}_{t}\right)\geq\frac{1}{6}\left\lVert\theta^{(t)}-\theta^{\star}_{t}\right\rVert_{\boldsymbol{\Sigma}^{\star}}^{2}\implies\left\lVert\theta^{(t)}-\theta^{\star}_{t}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\leq
5\sqrt{C_{\textup{ub}}}R.$

Finally, by using Lemma 3 on $\theta^{\star}_{t}$, we have
$\left\lVert\theta^{\star}_{t}-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\leq
4C_{\textup{est}}\sigma\sqrt{\kappa\epsilon}\leq R$. From this and the above
display, the triangle inequality yields (17). Now, (17) with Assumption 1.2
shows that for all $t$,

$\left\lVert\textup{Cov}_{\frac{1}{|G_{t}|}\mathbbm{1}_{G_{t}}}\left(\left\\{g_{i}(\theta^{(t)})\right\\}_{i\in
G_{t}}\right)\right\rVert_{\textup{op}}\leq
40C_{\textup{est}}C_{\textup{ub}}LR^{2}.$

We argue later in this proof that there are at most $O(\kappa)$ loops
throughout the algorithm. This shows all calls to $\mathsf{FastCovFilter}$ are
safe with respect to the bipartition $G_{t}$, $[n]\setminus G_{t}$ (adjusting
the definition of $\epsilon$ by a constant in Proposition 5), and thus taking
a union bound the algorithm succeeds with probability at least $1-\delta$.
Condition on this for the remainder of the proof.

The success of all calls to $\mathsf{FastCovFilter}$ implies that in every
iteration $t$ (following Proposition 5),

$\left\lVert\textup{Cov}_{w^{(t+1)}}\left(\left\\{g_{i}(\theta^{(t)})\right\\}_{i\in[n]}\right)\right\rVert_{\textup{op}}=O\left(LR^{2}\right).$
(18)

Next, in any iteration where $\Delta_{t+1}\leq\frac{R^{2}}{512\kappa
C_{\textup{id}}^{2}}$, we claim that

$\left\lVert\nabla
F_{w^{(t+1)}}\left(\theta^{(t)}\right)\right\rVert_{\left(\boldsymbol{\Sigma}^{\star}\right)^{-1}}\leq\frac{R}{4C_{\textup{id}}}.$
(19)

This is because $F_{w^{(t+1)}}$ is $8\kappa$-smooth in the
$\boldsymbol{\Sigma}^{\star}$ norm, since
$\nabla^{2}F_{w^{(t+1)}}=\mathbf{X}^{\top}\textbf{{diag}}\left(w^{(t+1)}\right)\mathbf{X}\preceq
8L\mathbf{I}$ by assumption, and $8\kappa\boldsymbol{\Sigma}^{\star}\succeq
8L\mathbf{I}$ by assumption, so by the guarantees of
$\mathcal{O}_{\textup{ERM}}$,

$\displaystyle\frac{1}{16\kappa}\left\lVert\nabla
F_{w^{(t+1)}}\left(\theta^{(t)}\right)\right\rVert_{\left(\boldsymbol{\Sigma}^{\star}\right)^{-1}}^{2}$
$\displaystyle\leq
F_{w^{(t+1)}}\left(\theta^{(t)}\right)-\min_{\theta\in\mathbb{R}^{d}}F_{w^{(t+1)}}\left(\theta\right)$
$\displaystyle\leq\Delta_{t+1}+O\left(\frac{\sigma^{2}}{\kappa}\right)\leq\frac{R^{2}}{256\kappa
C_{\textup{id}}^{2}}.$

Rearranging indeed yields (19). Now by combining (18) and (19) in Proposition
6, we see that if $\Delta_{t+1}\leq\frac{R^{2}}{512\kappa
C_{\textup{id}}^{2}}$ in an iteration, we obtain the desired
$\left\lVert\theta^{(t)}-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\leq\frac{1}{2}R$.

Runtime. We first observe that the loop in Lines 4-9 of Algorithm 1 can only
run $O(\kappa)$ times. This is because $\mathsf{FunctionFilter}$ decreases the
initial function value until it is $O(R^{2})$, and every loop decreases the
function value by $\Omega(\frac{R^{2}}{\kappa})$. Hence, the cost of the whole
algorithm is one call to $\mathsf{FunctionFilter}$, and $O(\kappa)$ calls to
$\mathcal{O}_{\textup{ERM}}$, $\mathsf{FastCovFilter}$, and two function value
computations, which fit into the allotted runtime budget by Lemma 6 and
Proposition 5.

∎

Finally, we remark that throughout Algorithm 1, we only filtered weight based
on $\mathsf{FunctionFilter}$ and $\mathsf{FastCovFilter}$, with respect to at
most $O(\kappa)$ distinct sets (in the manner described by Lemma 5): the sets
$\\{G_{t}\\}$ used in correctness calls to $\mathsf{FastCovFilter}$, and the
set $G^{\star}$ for the one call to $\mathsf{FunctionFilter}$.

###  3.4 Last phase analysis

In this section, we give a slight variant of Algorithm 1 which applies when
$R\leq C_{\textup{lp}}\sigma$ for a universal constant $C_{\textup{lp}}$. It
will have a somewhat more stringent termination condition, because we require
the gradient term in Proposition 6 to be $O(\sigma\sqrt{\epsilon\kappa})$, but
otherwise is identical.

Algorithm 2
$\mathsf{LastPhase}(\mathbf{X},y,\bar{w},\bar{\theta},D,\mathcal{O}_{\textup{ERM}},\delta)$

1: Input: Dataset
$\mathbf{X}=\left\\{X_{i}\right\\}_{i\in[n]}\in\mathbb{R}^{n\times d}$,
$y=\\{y_{i}\\}_{i\in[n]}\in\mathbb{R}^{n}$, satisfying Assumption 1,
$\epsilon$-saturated $\bar{w}\in\Delta^{n}$ with respect to bipartition
$[n]=G\cup B$ such that
$\mathbf{X}^{\top}\textbf{{diag}}\left(\bar{w}\right)\mathbf{X}\preceq
8L\mathbf{I}$, $\bar{\theta}\in\mathbb{R}^{d}$ with
$\left\lVert\bar{\theta}-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\leq
C_{\textup{lp}}\sigma$, $\delta\in(0,1)$, $O(\sigma^{2}\epsilon)$-approximate
ERM oracle $\mathcal{O}_{\textup{ERM}}$,  $\displaystyle D$
$\displaystyle\geq\max_{i\in[n]}\frac{1}{2}(\left\langle
X_{i},\bar{\theta}\right\rangle-y_{i})^{2}.$ (20)

2: Output: With probability $\geq 1-\delta$, saturated $w$ with respect to
$G\cup B$, and $\theta$ with
$\left\lVert\theta-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}=O\left(\sigma\sqrt{\kappa\epsilon}\right).$

3: $t\leftarrow 0$,
$w^{(0)}\leftarrow\mathsf{FunctionFilter}(w,\bar{\theta},C_{\textup{lp}}\sigma,D)$,
$\Delta_{0}\leftarrow\infty$

4: while $\Delta_{t}>\sigma^{2}\epsilon$ do

5: $\theta^{(t)}\leftarrow\mathcal{O}_{\textup{ERM}}(F_{w^{(t)}})$

6:
$w^{(t+1)}\leftarrow\mathsf{FastCovFilter}\left(\left\\{g_{i}\left(\theta^{(t)}\right)\right\\}_{i\in[n]},w^{(t)},O\left(\delta\epsilon\right),40C_{\textup{est}}C_{\textup{ub}}C_{\textup{lp}}^{2}L\sigma^{2}\right)$

7: $\Delta_{t+1}\leftarrow
F_{w^{(t+1)}}\left(\theta^{(t)}\right)-F_{w^{(t+1)}}\left(\theta^{(t+1)}\right)$

8: $t\leftarrow t+1$

9: end while

10: return $\left(w^{(t)},\theta^{(t-1)}\right)$

######  Lemma 8.

$\mathsf{LastPhase}$ is correct, i.e. if its preconditions are met, it
successfully returns $(w,\theta)$ such that $w\in\Delta^{n}$ is
$\epsilon$-saturated and
$\left\lVert\theta-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}=O\left(\sigma\sqrt{\kappa\epsilon}\right)$.
It runs in $O(\frac{1}{\epsilon})$ calls to $\mathcal{O}_{\textup{ERM}}$, plus

$O\left(n\log\left(\frac{D}{R^{2}}\right)+\frac{nd}{\epsilon}\log^{3}(n)\log\left(\frac{n}{\delta\epsilon}\right)\right)\textup{
additional time.}$

###### Proof.

On the correctness side, the analysis is nearly identical to Lemma 7; the same
logic applies to yield an analogous bound to (17), which shows that all
iterates are within $O(\sigma)$ from the minimizer, so all calls to
$\mathsf{FastCovFilter}$ are correct. This implies that (analogous to (18))
the gradient operator norm is always bounded by $O(L\sigma^{2})$. Similarly,
since the threshold for termination is when the function decrease is
$\sigma^{2}\epsilon$, we have that on the terminating iteration,

$\frac{1}{16\kappa}\left\lVert\nabla
F_{w^{(t+1)}}\left(\theta^{(t)}\right)\right\rVert_{\left(\boldsymbol{\Sigma}^{\star}\right)^{-1}}^{2}=O(\sigma^{2}\epsilon),$

and combining this with the operator norm bound in Proposition 6 yields the
conclusion. On the runtime side, the analysis is the same as Lemma 7, except
that there are now $O(\frac{1}{\epsilon})$ iterations. ∎

Again, we remark here that throughout Algorithm 2, we filtered weight with
respect to at most $O(\epsilon^{-1})$ distinct sets (in the manner described
by Lemma 5).

###  3.5 Full algorithm

Before we give our full algorithm, we require some preliminary pruning
procedures on the dataset.

######  Lemma 9.

Under Assumption 1, for all $i\in G$, $\left\lVert
X_{i}\right\rVert_{2}\leq\sqrt{2Ln}$.

###### Proof.

Suppose otherwise for some $i\in G$. Then, since
$\textup{Cov}_{w_{G}}\left(\mathbf{X}\right)\succeq\frac{1}{|G|}X_{i}X_{i}^{\top}\succeq\frac{1}{n}X_{i}X_{i}^{\top}$,
$\textup{Cov}_{w_{G}}\left(\mathbf{X}\right)$ has an eigenvalue larger than
$\frac{3}{2}L$ (certified by $\frac{X_{i}}{\left\lVert
X_{i}\right\rVert_{2}}$), contradicting Assumption 1.1. ∎

We next give a bound on $D$ required by Algorithms 1 and 2, assuming a bound
on
$\left\lVert\bar{\theta}-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}$.

######  Lemma 10.

Suppose all $\\{X_{i}\\}_{i\in[n]}$ satisfy $\left\lVert
X_{i}\right\rVert_{2}\leq\sqrt{2Ln}$, and we have a bound
$\left\lVert\bar{\theta}-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\leq
R$. Under Assumption 1, it suffices to set $D$ in $\mathsf{HalfRadiusLinReg}$
or $\mathsf{LastPhase}$ to

$\displaystyle D$ $\displaystyle=2nC_{\textup{ub}}\sigma^{2}+2\kappa nR^{2}.$
(21)

###### Proof.

First, under Assumption 1 we have that for all $i\in G$,

$\left|\left\langle X_{i},\theta^{\star}\right\rangle-
y_{i}\right|\leq\sqrt{2nC_{\textup{ub}}\sigma^{2}}.$

Applying Cauchy-Schwarz, we have that for all $i\in G$, since
$\left\lVert\bar{\theta}-\theta^{\star}\right\rVert_{2}\leq\frac{1}{\sqrt{\mu}}\left\lVert\bar{\theta}-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\leq\frac{R}{\sqrt{\mu}}$,

$\left|\left\langle X_{i},\bar{\theta}\right\rangle-
y_{i}\right|\leq\sqrt{2nC_{\textup{ub}}\sigma^{2}}+\left\lVert
X_{i}\right\rVert_{2}\left\lVert\bar{\theta}-\theta^{\star}\right\rVert_{2}\leq\sqrt{2nC_{\textup{ub}}\sigma^{2}}+\sqrt{2\kappa
n}R.$

∎

Finally, we give our full algorithm for regression, $\mathsf{FastRegression}$,
below.

Algorithm 3
$\mathsf{FastRegression}(\mathbf{X},y,\theta_{0},R_{0},\mathcal{O}_{\textup{ERM}},\delta)$

1: Input: Dataset
$\mathbf{X}=\left\\{X_{i}\right\\}_{i\in[n]}\in\mathbb{R}^{n\times d}$,
$y=\\{y_{i}\\}_{i\in[n]}\in\mathbb{R}^{n}$, satisfying Models 1 and 3, and
satisfying Assumptions 1, $\theta_{0}\in\mathbb{R}^{d}$ with
$\left\lVert\theta_{0}-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\leq
R_{0}$, $\delta\in(0,1)$, $O(\sigma^{2}\epsilon)$-approximate ERM oracle
$\mathcal{O}_{\textup{ERM}}$.

2: Output: With probability $\geq 1-\delta$, $\theta$ with
$\left\lVert\theta-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}=O\left(\sigma\sqrt{\kappa\epsilon}\right).$

3: Remove all $(X_{i},y_{i})$ with $\left\lVert
X_{i}\right\rVert_{2}>\sqrt{2nL}$, $n\leftarrow$ new dataset size

4:
$w\leftarrow\mathsf{FastCovFilter}(\mathbf{X},\frac{1}{n}\mathbbm{1},\frac{\delta}{3},\frac{3}{2}L,2nL)$,
$R\leftarrow R_{0}+4C_{\textup{est}}\sigma\sqrt{\kappa\epsilon}$,
$\theta\leftarrow\theta_{0}$

5: $T\leftarrow O(\log\frac{R_{0}}{\sigma})$ for a sufficiently large constant

6: while $R>C_{\textup{lp}}\sigma$ do

7: $D\leftarrow$ value in (21) with current setting of $R$

8: Remove all $(X_{i},y_{i})$ not satisfying bound (16), $n\leftarrow$ new
dataset size

9:
$(w,\theta)\leftarrow\mathsf{HalfRadiusLinReg}(\mathbf{X},y,w,\theta,R,D,\mathcal{O}_{\textup{ERM}},\frac{\delta}{3T})$

10: $R\leftarrow\frac{1}{2}R$

11: end while

12: $D\leftarrow$ value in (21) with current setting of $R$

13: Remove all $(X_{i},y_{i})$ not satisfying bound (16), $n\leftarrow$ new
dataset size

14:
$(w,\theta)\leftarrow\mathsf{HalfRadiusLinReg}(\mathbf{X},y,w,\theta,D,\mathcal{O}_{\textup{ERM}},\frac{\delta}{3})$

15: return $\theta$

######  Theorem 5.

In Models 1 and 3, under Assumption 1 with $r=\frac{\epsilon^{2}}{\alpha}$ for
$\alpha$ as in (12), supposing $\epsilon\kappa$ is sufficiently small, given
$\theta_{0}\in\mathbb{R}^{d}$ and
$\left\lVert\theta_{0}-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}\leq
R_{0}$, $\mathsf{FastRegression}$ returns $\theta$ with
$\left\lVert\theta-\theta^{\star}\right\rVert_{\boldsymbol{\Sigma}^{\star}}=O(\sigma\sqrt{\kappa\epsilon})$
with probability at least $1-\delta$. The algorithm runs in
$O\left(\frac{1}{\epsilon}+\kappa\log\frac{R_{0}}{\sigma}\right)$ calls to a
$O(\sigma^{2}\epsilon)$-approximate ERM oracle, and

$O\left(\left(nd\log^{3}(n)\log\left(\frac{nR_{0}}{\sigma\delta\epsilon}\right)\right)\left(\frac{1}{\epsilon}+\kappa\log\frac{R_{0}}{\sigma}\right)\right)\textup{
additional time.}$

###### Proof.

First, correctness of Lines 3, 8, and 13 of the algorithm follow from Lemmas 9
and 10. Also, Line 4 ensures that throughout the algorithm we have
$\mathbf{X}^{\top}\textbf{{diag}}\left(w\right)\mathbf{X}\preceq
8L\mathbf{I}$, so the preconditions of $\mathsf{HalfRadiusLinReg}$ and
$\mathsf{LastPhase}$ are met. Finally, the initial setting of $R$ is correct
by Lemma 3 and the assumed bound $R_{0}$. The correctness and runtime then
follow from applying Lemma 7 $T$ times and Lemma 8 once. The failure
probability comes from union bounding over the one call to
$\mathsf{FastCovFilter}$, the $T$ calls to $\mathsf{HalfRadiusLinReg}$, and
the one call to $\mathsf{LastPhase}$.

Finally, for completeness we check that the promise of Section 3.1 is kept by
Algorithm 3. There are at most $\log(\frac{R_{0}}{\sigma})$ calls to
$\mathsf{HalfRadiusLinReg}$, and one call to $\mathsf{LastPhase}$. Combined,
this accounts for at most
$O(\frac{1}{\epsilon}+\kappa\log\frac{R_{0}}{\sigma})$ distinct sets we
filtered with respect to, in the manner described by Lemma 5. For a
sufficiently large $\alpha$ in (12), this is indeed at most
$\frac{\alpha}{2\epsilon}$ distinct sets. ∎

We conclude this section by noting that the guarantees of Proposition 2 imply
the sample complexity required for Algorithm 3 to succeed with probability at
least $\frac{9}{10}-\delta$ is
$n=\widetilde{O}(\frac{d}{\epsilon^{4}}+\frac{d^{2}}{\epsilon^{3}})$.

##  4 Robust acceleration

In this section, we give a general-purpose algorithm for solving statistical
optimization problems with a finite condition number $\kappa$ under the strong
contamination model. We study the following abstract problem: we wish to
minimize a function $F:\mathbb{R}^{d}\rightarrow\mathbb{R}$ which is
$L$-smooth and $\mu$-strongly convex with minimizer $\theta^{\star}_{F}$, but
we only have black-box access to $F$ through a _noisy gradient oracle_
$\mathcal{O}_{\textup{ng}}$. In particular, we can query
$\mathcal{O}_{\textup{ng}}$ at any point $\theta\in\mathbb{R}^{d}$ with a
parameter $R\geq\left\lVert\theta-\theta^{\star}_{F}\right\rVert_{2}$ and
receive $G(\theta)$ such that for a universal constant $C_{\textup{ng}}$,

$\left\lVert G(\theta)-\nabla F(\theta)\right\rVert_{2}\leq
C_{\textup{ng}}\left(\sqrt{L\epsilon}\sigma+L\sqrt{\epsilon}R\right).$

Notably, our algorithm is _accelerated_ , running in a number of iterations
depending on $\sqrt{\kappa}$ rather than $\kappa$. It applies to both the
regression setting of Section 2.3 and the smooth stochastic optimization
setting of Section 2.4. We demonstrate in Section 4.1 how to build a noisy
gradient oracle for regression and smooth stochastic optimization settings. We
then build in Section LABEL:ssec:proxproblem a simple subroutine based on the
robust gradient descent framework of [PSBR20] to approximately solve
_regularized subproblems_ encountered by our final algorithm. We put the
pieces together and give our complete algorithm in Sections 4.2 and 4.3.
Throughout we assume $\epsilon\kappa^{2}<1$ is sufficiently small.

###  4.1 Noisy gradient oracle

In this section, we build noisy gradient oracles for the problems in Sections
2.3 and 2.4. We now give a formal definition below; the remaining sections
will access $F$ through this abstraction.

######  Definition 4 (Noisy gradient oracle).

We call $\mathcal{O}_{\textup{ng}}(\theta,R)$ a $(L,\sigma,\delta)$-_noisy
gradient oracle_ for $F:\mathbb{R}^{d}\rightarrow\mathbb{R}$ with minimizer
$\theta^{\star}_{F}$ if on query $\theta\in\mathbb{R}^{d}$ and given
$R\geq\left\lVert\theta-\theta^{\star}_{F}\right\rVert_{2}$, with probability
$\geq 1-\delta$ it returns $G(\theta)$ satisfying for a universal constant
$C_{\textup{ng}}$,

$\left\lVert G(\theta)-\nabla F(\theta)\right\rVert_{2}\leq
C_{\textup{ng}}\left(\sqrt{L\epsilon}\sigma+L\sqrt{\epsilon}R\right).$

If the returned $G(\theta)$ always satisfies the stronger bound $\left\lVert
G(\theta)-\nabla F(\theta)\right\rVert_{2}\leq
C_{\textup{ng}}\sqrt{L\epsilon}\sigma$, we call $\mathcal{O}_{\textup{ng}}$ a
$(L,\sigma,\delta)$-_radiusless noisy gradient oracle_.

Before developing our implementations, we state a useful identifiability
result relating gradient estimation to controlling operator norms of gradient
second moments for finite sum functions.

######  Lemma 11.

Suppose $F_{G}(\theta)=\frac{1}{|G|}\sum_{i\in G}f_{i}(\theta)$ for some
functions $\\{f_{i}\\}_{i\in G}$, and let $w\in\Delta^{n}$ be saturated with
respect to bipartition $[n]=G\cup B$. For $\tilde{w}:=\frac{w}{\left\lVert
w\right\rVert_{1}}$ and $w^{\star}_{G}=\frac{1}{|G|}\mathbbm{1}_{G}$, we have

$\left\lVert\operatorname*{\mathbb{E}}_{i\sim w^{\star}_{G}}\left[\nabla
f_{i}(\theta)\right]-\operatorname*{\mathbb{E}}_{j\sim\tilde{w}}\left[\nabla
f_{j}(\theta)\right]\right\rVert_{2}\leq\sqrt{24\epsilon}\left(\left\lVert\textup{Cov}_{w^{\star}_{G}}\left(\left\\{\nabla
f_{i}(\theta)\right\\}_{i\in[n]}\right)\right\rVert_{\textup{op}}^{\frac{1}{2}}+\left\lVert\textup{Cov}_{\tilde{w}}\left(\left\\{\nabla
f_{i}(\theta)\right\\}_{i\in[n]}\right)\right\rVert_{\textup{op}}^{\frac{1}{2}}\right).$

###### Proof.

Throughout this proof, we let $\mathcal{C}$ supported on $[n]\times[n]$ be an
optimal coupling between $i\sim w^{\star}_{G}$ and $j\sim\tilde{w}$, and
follow notation from Proposition 6. For some unit vector $v\in\mathbb{R}^{d}$,
we have

$\displaystyle\left\langle v,\operatorname*{\mathbb{E}}_{i\sim
w^{\star}_{G}}\left[\nabla
f_{i}(\theta)\right]-\operatorname*{\mathbb{E}}_{j\sim\tilde{w}}\left[\nabla
f_{j}(\theta)\right]\right\rangle$
$\displaystyle=\operatorname*{\mathbb{E}}_{i,j\sim\mathcal{C}}\left[\left\langle
v,\nabla f_{i}(\theta)-\nabla f_{j}(\theta)\right\rangle\right]$
$\displaystyle=\operatorname*{\mathbb{E}}_{i,j\sim\mathcal{C}}\left[\left\langle
v,\nabla f_{i}(\theta)-\nabla f_{j}(\theta)\right\rangle\mathbbm{1}_{i\neq
j}\right]$
$\displaystyle\leq\sqrt{6\epsilon}\operatorname*{\mathbb{E}}_{i,j\sim\mathcal{C}}\left[\left\langle
v,\nabla f_{i}(\theta)-\nabla
f_{j}(\theta)\right\rangle^{2}\right]^{\frac{1}{2}}$
$\displaystyle\leq\sqrt{24\epsilon}\left(\operatorname*{\mathbb{E}}_{i\sim
w^{\star}_{G}}\left[\left\langle v,\nabla
f_{i}(\theta)\right\rangle^{2}\right]^{\frac{1}{2}}+\operatorname*{\mathbb{E}}_{j\sim\tilde{w}}\left[\left\langle
v,\nabla f_{j}(\theta)\right\rangle^{2}\right]^{\frac{1}{2}}\right).$

The conclusion follows from choosing $v$ to be in the direction of
$\operatorname*{\mathbb{E}}_{i\sim w^{\star}_{G}}\left[\nabla
f_{i}(\theta)\right]-\operatorname*{\mathbb{E}}_{j\sim\tilde{w}}\left[\nabla
f_{j}(\theta)\right]$, and using the definition of the operator norm. ∎

Lemma 11 implies that for approximating gradients of functions $F$ which are
“closely approximated” by an (unknown) finite sum function $F_{G}$, it
suffices to find a weighting $\tilde{w}$ such that the operator norm of
$\textup{Cov}_{\tilde{w}}$ applied to gradients is bounded. We now demonstrate
applications of this strategy to linear regression and smooth stochastic
optimization.

######  Corollary 1.

Consider a robust linear regression instance where we have sample access to
datasets $\mathbf{X}=\\{X_{i}\\}_{i\in[n]}\in\mathbb{R}^{n\times d}$ and
$y=\\{y_{i}\\}_{i\in[n]}\in\mathbb{R}^{n}$ under Models 1, 4, with sample size
$n$ corresponding to Proposition 3. For

$F(\theta)=F^{\star}(\theta)=\operatorname*{\mathbb{E}}_{X,y\sim\mathcal{D}_{Xy}}\left[\frac{1}{2}(\left\langle
X_{i},\theta\right\rangle-y_{i})^{2}\right],$

we can construct a $(L,\sigma,\delta)$-noisy gradient oracle for $F$ in
$O\left(nd\log^{3}(n)\log^{2}\left(\frac{n}{\delta}\right)\right)$ time, using
$O(\log\frac{1}{\delta})$ queries of samples from Proposition 3.

###### Proof.

We first demonstrate how to construct a noisy gradient oracle with success
probability $\geq\frac{8}{10}$. The algorithm is as follows: first, sample a
dataset under Models 1, 4, according to Proposition 3. Then, at point
$\theta\in\mathbb{R}^{d}$, with probability $\frac{9}{10}$ Assumption 2 gives
us a set $G:=G_{\theta}$ (where we drop the subscript for simplicity, as this
proof only discusses a single $\theta$) with $|G|=(1-\epsilon)$ such that (10)
holds. If we have the promise
$\left\lVert\theta-\theta^{*}\right\rVert_{2}\leq R$, let $w\in\Delta^{n}$ be
the output of

$\mathsf{FastCovFilter}\left(\\{g_{i}(\theta)\\}_{i\in[n]},\frac{1}{n}\mathbbm{1},\frac{9}{10},C_{\textup{est}}\left(L\sigma^{2}+LR^{2}\right)\right),\text{
where }g_{i}(\theta):=X_{i}\left(\left\langle X_{i},\theta\right\rangle-
y_{i}\right).$

where $\mathsf{FastCovFilter}$ (Algorithm 7) is the algorithm of Proposition
5. We then output
$\operatorname*{\mathbb{E}}_{j\sim\tilde{w}}\left[g_{j}(\theta)\right]$, where
$\tilde{w}=\frac{w}{\left\lVert w\right\rVert_{1}}$. The runtime is
$O(nd\log^{4}n)$ from the bottleneck operation of running
$\mathsf{FastCovFilter}$. The assumptions of $\mathsf{FastCovFilter}$, namely
a bound on
$\textup{Cov}_{w^{\star}_{G}}\left(\\{g_{i}(\theta)\\}_{i\in[n]}\right)$, are
satisfied by Assumption 2.2. Guarantees of $\mathsf{FastCovFilter}$ and Lemma
11 then imply

$\left\lVert\operatorname*{\mathbb{E}}_{i\sim
w^{\star}_{G}}\left[g_{i}(\theta)\right]-\operatorname*{\mathbb{E}}_{j\sim\tilde{w}}\left[g_{j}(\theta)\right]\right\rVert_{2}=O\left(\sqrt{L\epsilon}\sigma+L\sqrt{\epsilon}R\right).$

The conclusion follows from Assumption 2.2 which bounds
$\left\lVert\operatorname*{\mathbb{E}}_{i\sim
w^{\star}_{G}}\left[g_{i}(\theta)\right]-\nabla F(\theta)\right\rVert_{2}$.

We now describe how to boost the success probability, by calling our sample
access $T=O(\log\frac{1}{\delta})$ times. Let $G^{\star}:=\nabla F(\theta)$ be
the true gradient, and run the procedure described above $T$ times, producing
$\\{G_{t}\\}_{t\in[T]}$, such that each $G_{t}$ satisfies $\left\lVert
G_{t}-G^{\star}\right\rVert_{2}\leq
C(\sqrt{L\epsilon}\sigma+L\sqrt{\epsilon}R)$ with probability at least
$\frac{4}{5}$, for some constant $C$. By standard binomial concentration, with
probability at least $1-\delta$, at least $\frac{3}{5}$ of the
$\\{G_{t}\\}_{t\in[T]}$ will satisfy this bound; call such a satisfying
$G_{t}$ “good.” We return any $G_{t}$ which is within distance
$2C(\sqrt{L\epsilon}\sigma+L\sqrt{\epsilon}R)$ from at least $\frac{3}{5}$ of
the gradient estimates. Note this will never return any $G_{t}$ with
$\left\lVert
G_{t}-G^{\star}\right\rVert_{2}>4C(\sqrt{L\epsilon}\sigma+L\sqrt{\epsilon}R)$,
since the triangle inequality implies this $G_{t}$ will miss all the good
estimates, a contradiction since there is at most a $\frac{2}{5}$ fraction
which is not good. Thus, this procedure satisfies the requirements with
$C_{\textup{ng}}=4C$; the additional runtime overhead is
$O(\log^{2}(\frac{1}{\delta}))$ distance comparisons between our gradient
estimates. ∎

######  Corollary 2.

Consider a robust Lipschitz (not necessarily smooth) stochastic optimization
instance where we have sample access to datasets
$\mathbf{X}=\\{X_{i}\\}_{i\in[n]}\in\mathbb{R}^{n\times d}$ and
$y=\\{y_{i}\\}_{i\in[n]}\in\mathbb{R}^{n}$ under Models 1, 2, 6 with sample
size $n$ corresponding to Proposition 4. For

$F(\theta)=F^{\star}(\theta)+\frac{\mu}{2}\left\lVert\theta\right\rVert_{2}^{2}=\operatorname*{\mathbb{E}}_{f\sim\mathcal{D}_{f}}\left[f(\theta)\right]+\frac{\mu}{2}\left\lVert\theta\right\rVert_{2}^{2},$

we can construct a $(L,1,\delta)$-radiusless noisy gradient oracle for $F$ in
$O\left(nd\log^{3}(n)\log^{2}\left(\frac{n}{\delta}\right)\right)$ time, using
$O(\log\frac{1}{\delta})$ queries of samples from Proposition 4.

###### Proof.

The proof follows identically to that of Corollary 1, where we use Assumption
3.2 in place of Assumption 2.2. ∎

In most of Sections 4.2 and 4.3, we will no longer discuss any specifics of
the unknown function $F:\mathbb{R}^{d}\to\mathbb{R}$ we wish to optimize,
except that it is $L$-smooth, $\mu$-strongly convex, has minimizer
$\theta^{\star}_{F}$, and supports a noisy gradient oracle
$\mathcal{O}_{\textup{ng}}$. We will apply Corollaries 1 and 2 to derive
concrete rates and sample complexities for specific applications at the
conclusion of this section.

###  4.2 Halving the distance to $\theta^{\star}_{F}$

In this section, we give a subroutine used in our full algorithm, which halves
the distance to $\theta^{\star}_{F}$, the minimizer of $F$, outside of a
sufficiently large radius. Suppose that we have an initial point
$\bar{\theta}\in\mathbb{R}^{d}$, as well as a sufficiently large scalar $R\geq
0$ with the promise that (for a universal constant $C_{\textup{lp}}$)

$\left\lVert\bar{\theta}-\theta^{\star}_{F}\right\rVert_{2}\leq R,\text{ where
}R\geq C_{\textup{lp}}\sigma\sqrt{\frac{\kappa\epsilon}{\mu}}.$ (22)

We begin by stating a standard lemma from convex analysis, following from
first-order optimality.

######  Lemma 12 (Proximal three-point inequality).

Let $f$ be a convex function, and let $\mathcal{S}$ be a convex set. For any
point $y\in\mathcal{S}$, define
$\textup{Prox}(y):=\textup{argmin}_{x\in\mathcal{S}}\left\\{f(x)+\left\lVert
x-y\right\rVert_{2}^{2}\right\\}$. Then if $x=\textup{Prox}(y)$,

$\left\langle\nabla f(x),x-u\right\rangle\leq\frac{1}{2}\left(\left\lVert
y-u\right\rVert_{2}^{2}-\left\lVert x-u\right\rVert_{2}^{2}-\left\lVert
x-y\right\rVert_{2}^{2}\right),\text{ for all }u\in\mathcal{S}.$

We now state a procedure, $\mathsf{HalfRadiusAccel}$, which returns a new
$\theta$ with
$\left\lVert\theta-\theta^{\star}_{F}\right\rVert_{2}\leq\frac{1}{2}R$. In its
statement, we define a sequence of scalars $\\{a_{t},A_{t}\\}_{0\leq t<T}$
given by the recursions

$A_{0}=0,\;A_{t}=3a_{t}^{2},\;A_{t+1}=A_{t}+a_{t+1}.$ (23)

The following fact is well-known (see for instance Chapter 2.2 of [Nes03]).

######  Fact 1.

For all $0\leq t<T$, $A_{t}=\Theta(t^{2})$ and $a_{t}=\Theta(t)$.

Algorithm 4
$\mathsf{HalfRadiusAccel}(\bar{\theta},R,\mathcal{O}_{\textup{ng}},\delta)$

1: Input: $\mathcal{O}_{\textup{ng}}$, a $(L,\sigma,\frac{\delta}{T})$-noisy
gradient oracle for $L$-smooth, $\mu$-strongly convex
$F:\mathbb{R}^{d}\to\mathbb{R}$ with minimizer $\theta^{\star}_{F}$,
$\bar{\theta}\in\mathbb{R}^{d}$, $R\in\mathbb{R}_{\geq 0}$ satisfying (22),
and $T=O(\sqrt{\kappa})$, $\delta\in(0,1)$

2: Output: With probability $\geq 1-\delta$, $\theta\in\mathbb{R}^{d}$ with
$\left\lVert\theta-\theta^{\star}_{F}\right\rVert_{2}\leq\frac{1}{2}R$

3: $T\leftarrow O\left(\sqrt{\kappa}\right)$ for a sufficiently large
constant, $t\leftarrow 0$, $\theta_{0}\leftarrow\bar{\theta}$,
$v_{0}\leftarrow\bar{\theta}$

4: while $t<T$ do

5:
$y_{t}\leftarrow\frac{A_{t}}{A_{t+1}}\theta_{t}+\frac{a_{t+1}}{A_{t+1}}v_{t}$

6: $g_{t}\leftarrow\mathcal{O}_{\textup{ng}}(y_{t},2R)$

7:
$\theta_{t+1}\leftarrow\textup{argmin}_{\theta\in\mathbb{B}}\left\\{\frac{1}{3L}\left\langle
g_{t},\theta\right\rangle+\frac{1}{2}\left\lVert\theta-
y_{t}\right\rVert_{2}^{2}\right\\}$, where
$\mathbb{B}:=\\{\theta\mid\left\lVert\theta-\bar{\theta}\right\rVert_{2}\leq
R\\}$

8:
$v_{t+1}\leftarrow\textup{argmin}_{v\in\mathbb{B}}\left\\{\frac{a_{t+1}}{L}\left\langle
g_{t},v\right\rangle+\frac{1}{2}\left\lVert
v-v_{t}\right\rVert_{2}^{2}\right\\}$

9: $t\leftarrow t+1$

10: end while

11: return $\theta_{t}$

We remark throughout we assume that Lines 7 and 8 are implemented exactly for
simplicity; it is straightforward to verify from the proof of Lemma 13 that it
suffices to implement these steps to inverse-polynomial precision in problem
parameters. This can be done by a standard binary search (see e.g. Proposition
8 of [CJJ+20]), and is not the bottleneck operation compared to calling
$\mathcal{O}_{\textup{ng}}$.

We give the main technical lemma of this section, which shows a potential
bound on iterates of $\mathsf{HalfRadiusAccel}$. Our proof is based on a
standard analysis of accelerated methods by [AO17].

######  Lemma 13.

In every iteration $0\leq t\leq T$, define

$E_{t}:=F(\theta_{t})-F(\theta^{\star}_{F}),\;D_{t}:=\frac{1}{2}\left\lVert
v_{t}-\theta^{\star}_{F}\right\rVert_{2}^{2},\;\Phi_{t}:=\frac{A_{t}}{L}E_{t}+D_{t}.$

Then, for all $0\leq t<T$, for a universal constant $C_{\textup{pot}}$,

$\Phi_{t+1}-\Phi_{t}\leq
C_{\textup{pot}}\left(t\sigma\sqrt{\frac{\epsilon}{L}}R+t\sqrt{\epsilon}R^{2}+t^{2}\sigma^{2}\frac{\epsilon}{L}+t^{2}\epsilon
R^{2}\right).$

###### Proof.

Throughout this proof Lemma 12 will be applied to the set
$\mathcal{S}=\mathbb{B}$. Fix an iteration $t$. We observe that $\theta_{t}$
and $v_{t}$ lie in $\mathbb{B}$ by the constraints on Lines 7 and 8: therefore
$y_{t}\in\mathbb{B}$. As
$\left\lVert\bar{\theta}-\theta^{\star}_{F}\right\rVert_{2}\leq R$, the
triangle inequality yields $\left\lVert
y_{t}-\theta^{\star}_{F}\right\rVert_{2}\leq\left\lVert\bar{\theta}-\theta^{\star}_{F}\right\rVert_{2}+\left\lVert\bar{\theta}-y_{t}\right\rVert_{2}\leq
2R$: thus by the guarantee of $\mathcal{O}_{\textup{ng}}$ we have, for some
$C_{\textup{ng}}=O(1)$ (adjusting the definition of $C_{\textup{ng}}$ by a
constant)

$\left\lVert g_{t}-\nabla
F\left(\theta^{\star}_{t+1}\right)\right\rVert_{2}\leq
C_{\textup{ng}}\left(\sqrt{L\epsilon}\sigma+L\sqrt{\epsilon}R\right).$

For convenience, we define
$\rho=C_{\textup{ng}}\left(\sqrt{L\epsilon}\sigma+L\sqrt{\epsilon}R\right)$.
We define the helper function

$\mathsf{Prog}(y;g)=\min_{\theta\in\mathbb{B}}\left\\{\left\langle
g,\theta-y\right\rangle+\frac{3L}{2}\left\lVert\theta-y\right\rVert_{2}^{2}\right\\},$

and observe

$\displaystyle\mathsf{Prog}(y_{t};g_{t})$
$\displaystyle=\min_{\theta\in\mathbb{B}}\left\\{\left\langle g_{t},\theta-
y_{t}\right\rangle+\frac{3L}{2}\left\lVert\theta-
y_{t}\right\rVert_{2}^{2}\right\\}$ (24)
$\displaystyle\begin{subarray}{c}(i)\\\
=\end{subarray}\frac{3L}{2}\left\lVert\theta_{t+1}-y_{t}\right\rVert^{2}+\left\langle
g_{t},\theta_{t+1}-y_{t}\right\rangle$
$\displaystyle=\left(\frac{L}{2}\left\lVert\theta_{t+1}-y_{t}\right\rVert^{2}+\left\langle\nabla
F(y_{t}),\theta_{t+1}-y_{t}\right\rangle\right)+L\left\lVert\theta_{t+1}-y_{t}\right\rVert_{2}^{2}-\left\langle\nabla
F(y_{t})-g_{t},\theta_{t+1}-y_{t}\right\rangle$
$\displaystyle\begin{subarray}{c}(ii)\\\
\geq\end{subarray}F(\theta_{t+1})-F(y_{t})+L\left\lVert\theta_{t+1}-y_{t}\right\rVert_{2}^{2}-\left\langle\nabla
F(y_{t})-g_{t},\theta_{t+1}-y_{t}\right\rangle$
$\displaystyle\begin{subarray}{c}(iii)\\\
\geq\end{subarray}F(\theta_{t+1})-F(y_{t})-\frac{1}{4L}\left\lVert\nabla
F(y_{t})-g_{t}\right\rVert_{2}^{2}\geq
F(\theta_{t+1})-F(y_{t})-\frac{\rho^{2}}{4L}.$

Here $(i)$ is by the optimality of $\theta_{t+1}$, $(ii)$ holds via the
$L$-smoothness of $F$, and $(iii)$ follows from Young’s inequality $\langle
a,b\rangle-\frac{1}{2}\left\lVert
a\right\rVert_{2}^{2}\leq\frac{1}{2}\left\lVert b\right\rVert_{2}^{2}$ and our
bound on $\left\lVert\nabla F(y_{t})-g_{t}\right\rVert_{2}$. Next, we note

$\frac{a_{t+1}}{L}\left\langle
g_{t},v_{t+1}-\theta^{\star}_{F}\right\rangle\leq\frac{1}{2}\left(\left\lVert
v_{t}-\theta^{\star}_{F}\right\rVert_{2}^{2}-\left\lVert
v_{t+1}-\theta^{\star}_{F}\right\rVert_{2}^{2}-\left\lVert
v_{t}-v_{t+1}\right\rVert_{2}^{2}\right)$ (25)

by the proximal three-point inequality Lemma 12 on Line 8. Define

$\widetilde{y}_{t}=\frac{A_{t}}{A_{t+1}}\theta_{t}+\frac{a_{t+1}}{A_{t+1}}v_{t+1},$

and note that
$y_{t}-\widetilde{y}_{t}=\frac{a_{t+1}}{A_{t+1}}\left(v_{t}-v_{t+1}\right)$
and $\widetilde{y}_{t}\in\mathbb{B}$ by convexity. Consequently,

$\displaystyle\frac{a_{t+1}}{L}\left\langle
g_{t},v_{t}-\theta^{\star}_{F}\right\rangle$
$\displaystyle=\frac{a_{t+1}}{L}\left\langle
g_{t},v_{t}-v_{t+1}\right\rangle+\frac{a_{t+1}}{L}\left\langle
g_{t},v_{t+1}-\theta^{\star}_{F}\right\rangle$
$\displaystyle\begin{subarray}{c}(i)\\\
\leq\end{subarray}\frac{a_{t+1}}{L}\left\langle
g_{t},v_{t}-v_{t+1}\right\rangle+\frac{1}{2}\left(\left\lVert
v_{t}-\theta^{\star}_{F}\right\rVert_{2}^{2}-\left\lVert
v_{t+1}-\theta^{\star}_{F}\right\rVert_{2}^{2}-\left\lVert
v_{t}-v_{t+1}\right\rVert_{2}^{2}\right)$
$\displaystyle\begin{subarray}{c}(ii)\\\
=\end{subarray}\frac{A_{t+1}}{L}\left\langle
g_{t},y_{t}-\widetilde{y}_{t}\right\rangle-\frac{A_{t+1}^{2}}{2a_{t+1}^{2}}\left\lVert
y_{t}-\widetilde{y}_{t}\right\rVert_{2}^{2}+D_{t}-D_{t+1}$
$\displaystyle\begin{subarray}{c}(iii)\\\
=\end{subarray}\frac{A_{t+1}}{L}\left(\left\langle
g_{t},y_{t}-\widetilde{y}_{t}\right\rangle-\frac{3L}{2}\left\lVert
y_{t}-\widetilde{y}_{t}\right\rVert_{2}^{2}\right)+D_{t}-D_{t+1}$
$\displaystyle\begin{subarray}{c}(iv)\\\
\leq\end{subarray}-\frac{A_{t+1}}{L}\mathsf{Prog}(y_{t};g_{t})+D_{t}-D_{t+1}$
$\displaystyle\begin{subarray}{c}(v)\\\
\leq\end{subarray}-\frac{A_{t+1}}{L}\left(F(\theta_{t+1})-F(y_{t})-\frac{\rho^{2}}{4L}\right)+D_{t}-D_{t+1}.$

Here $(i)$ uses (25), $(ii)$ uses the definition of $\widetilde{y}_{t}$,
$(iii)$ uses that $A_{t+1}=3a_{t+1}^{2}$, $(iv)$ uses the definition of
$\mathsf{Prog}$, and $(v)$ uses our lower bound on
$\mathsf{Prog}(y_{t};g_{t})$ (24). Finally, by convexity of $F$ we have

$\displaystyle\frac{a_{t+1}}{L}\left(F(y_{t})-F(\theta^{\star}_{F})\right)$
$\displaystyle\leq\frac{a_{t+1}}{L}\left\langle\nabla
F(y_{t}),y_{t}-\theta^{\star}_{F}\right\rangle$
$\displaystyle=\frac{a_{t+1}}{L}\left\langle\nabla
F(y_{t}),y_{t}-v_{t}\right\rangle+\frac{a_{t+1}}{L}\left\langle\nabla
F(y_{t}),v_{t}-\theta^{\star}_{F}\right\rangle$
$\displaystyle\begin{subarray}{c}(i)\\\
=\end{subarray}\frac{A_{t}}{L}\left\langle\nabla
F(y_{t}),\theta_{t}-y_{t}\right\rangle+\frac{a_{t+1}}{L}\left\langle\nabla
F(y_{t}),v_{t}-\theta^{\star}_{F}\right\rangle$
$\displaystyle\leq\frac{A_{t}}{L}\left(F(\theta_{t})-F(y_{t})\right)+\frac{a_{t+1}}{L}\left\langle
g_{t},v_{t}-\theta^{\star}_{F}\right\rangle+\frac{a_{t+1}}{L}\left\langle\nabla
F(y_{t})-g_{t},v_{t}-\theta^{\star}_{F}\right\rangle$
$\displaystyle\begin{subarray}{c}(ii)\\\
\leq\end{subarray}\frac{A_{t}}{L}\left(F(\theta_{t})-F(y_{t})\right)+\frac{a_{t+1}}{L}\left\langle
g_{t},v_{t}-\theta^{\star}_{F}\right\rangle+\frac{a_{t+1}}{L}\left\lVert\nabla
F(y_{t})-g_{t}\right\rVert_{2}\left\lVert
v_{t}-\theta^{\star}_{F}\right\rVert_{2}$
$\displaystyle\begin{subarray}{c}(iii)\\\
\leq\end{subarray}\frac{A_{t}}{L}\left(F(\theta_{t})-F(y_{t})\right)+\frac{a_{t+1}}{L}\left\langle
g_{t},v_{t}-\theta^{\star}_{F}\right\rangle+\frac{2a_{t+1}\rho R}{L}$

where $(i)$ used
$y_{t}-v_{t}=\frac{A_{t}}{a_{t+1}}\left(\theta_{t}-y_{t}\right)$, $(ii)$ used
the Cauchy-Schwarz inequality, and $(iii)$ used that $\left\lVert
v_{t}-\theta^{\star}_{F}\right\rVert_{2}\leq\left\lVert
v_{t}-\bar{\theta}\right\rVert_{2}+\left\lVert\bar{\theta}-\theta^{\star}_{F}\right\rVert_{2}\leq
2R$. Combining the above two equations and rearranging,

$\displaystyle\frac{a_{t+1}}{L}\left(F(y_{t})-F(\theta^{\star}_{F})\right)+\frac{A_{t+1}}{L}\left(F(\theta_{t+1})-F(y_{t})\right)$
$\displaystyle\leq\frac{A_{t}}{L}\left(F(\theta_{t})-F(y_{t})\right)$
$\displaystyle+D_{t}-D_{t+1}+\frac{2a_{t+1}\rho
R}{L}+\frac{A_{t+1}\rho^{2}}{4L^{2}}.$

Adding $\frac{A_{t}}{L}\left(F(y_{t})-F(\theta^{\star}_{F})\right)$ to both
sides, we obtain

$\Phi_{t+1}-\Phi_{t}\leq\frac{2a_{t+1}\rho
R}{L}+\frac{A_{t+1}\rho^{2}}{4L^{2}}.$

Finally, applying Fact 1, we see that this potential increase is indeed
bounded as

$\displaystyle\frac{2a_{t+1}\rho R}{L}+\frac{A_{t+1}\rho^{2}}{4L^{2}}$
$\displaystyle=O\left(\frac{t\rho
R}{L}+\frac{t^{2}\rho^{2}}{L^{2}}\right)=O\left(t\sigma\sqrt{\frac{\epsilon}{L}}R+t\sqrt{\epsilon}R^{2}+t^{2}\sigma^{2}\frac{\epsilon}{L}+t^{2}\epsilon
R^{2}\right).$

∎

Finally, we are ready to analyze the output of $\mathsf{HalfRadiusAccel}$.

######  Lemma 14.

With probability at least $1-\delta$, the output $\theta_{T}$ of
$\mathsf{HalfRadiusAccel}$ satisfies

$\left\lVert\theta_{T}-\theta^{\star}_{F}\right\rVert_{2}\leq\frac{1}{2}R.$

###### Proof.

The failure probability comes from union bounding over the failures of calls
to $\mathcal{O}_{\textup{ng}}$, so we discuss correctness assuming all calls
succeed. By telescoping Lemma 13 over $T$ iterations, we have

$\Phi_{T}\leq\Phi_{0}+C_{\textup{pot}}\sum_{t\in[T]}\left(t\sigma\sqrt{\frac{\epsilon}{L}}R+t\sqrt{\epsilon}R^{2}+t^{2}\sigma^{2}\frac{\epsilon}{L}+t^{2}\epsilon
R^{2}\right).$

It is clear from definition that $\Phi_{0}\leq\frac{1}{2}R^{2}$, so it remains
to bound all other terms. By examination,

$\displaystyle\sum_{t\in[T]}t\sigma\sqrt{\frac{\epsilon}{L}}R=O\left(\sigma\kappa\sqrt{\frac{\epsilon}{L}}R\right)=O\left(R^{2}\right),$
$\displaystyle\sum_{t\in[T]}t\sqrt{\epsilon}R^{2}=O\left(\kappa\sqrt{\epsilon}R^{2}\right)=O\left(R^{2}\right),$
$\displaystyle\sum_{t\in[T]}t^{2}\sigma^{2}\frac{\epsilon}{L}=O\left(\kappa^{1.5}\sigma^{2}\frac{\epsilon}{L}\right)=O\left(R^{2}\right),$
$\displaystyle\sum_{t\in[T]}t^{2}\epsilon R^{2}=O\left(\kappa^{1.5}\epsilon
R^{2}\right)=O\left(R^{2}\right).$

Each of the above lines follows from $\epsilon\kappa^{2}$ being sufficiently
small, and the lower bound in (22). Thus for a large enough value of
$C_{\textup{lp}}$ in (22), we have that $\Phi_{T}\leq R^{2}$. Since
$\Phi_{T}=A_{T}E_{T}+D_{T}\geq A_{T}E_{T}$, choosing a sufficiently large
value of $T=\sqrt{\kappa}$ combined with Fact 1 yields

$E_{T}\leq\frac{R^{2}}{A_{T}}\leq\frac{LR^{2}}{8\kappa}=\frac{\mu R^{2}}{8}.$

The conclusion follows from strong convexity of $F$, which implies
$E_{T}\geq\frac{\mu}{2}\left\lVert\theta_{T}-\theta^{\star}_{F}\right\rVert_{2}^{2}$.

A note on constants. To check there are no conflict of interests hidden in the
constants of this proof, note first that conditional on the bound at time $T$
being at most $R^{2}$, the number of iterations $T$ can be chosen solely as a
function of the constants in Fact 1. From this point, the constant
$C_{\textup{lp}}$ in (22) can be chosen to ensure that the potential bound is
indeed $R^{2}$. ∎

###  4.3 Full accelerated algorithm

We conclude this section with a statement of a complete accelerated algorithm,
and its applications.

Algorithm 5
$\mathsf{RobustAccel}(\theta_{0},R_{0},\mathcal{O}_{\textup{ng}},\delta)$

1: Input: $\mathcal{O}_{\textup{ng}}$, a $(L,\sigma,\frac{\delta}{N})$-noisy
gradient oracle for $L$-smooth, $\mu$-strongly convex
$F:\mathbb{R}^{d}\to\mathbb{R}$ with minimizer $\theta^{\star}_{F}$ for
$N=O\left(\sqrt{\kappa}\log\left(\frac{R_{0}\sqrt{L}}{\sigma\sqrt{\epsilon}}\right)\right)$,
$\theta_{0}\in\mathbb{R}^{d}$ with
$\left\lVert\theta_{0}-\theta^{\star}_{F}\right\rVert_{2}\leq R_{0}$,
$\delta\in(0,1)$

2: Output: With probability $\geq 1-\delta$, $\theta\in\mathbb{R}^{d}$ with
$\left\lVert\theta-\theta^{\star}_{F}\right\rVert_{2}=O\left(\sigma\sqrt{\frac{\kappa\epsilon}{\mu}}\right).$
(26)

3: $T\leftarrow
O\left(\log\left(\frac{R_{0}\sqrt{\mu}}{\sigma\sqrt{\kappa\epsilon}}\right)\right)$
for a sufficiently large constant, $t\leftarrow 0$

4: while $t<T$ do

5:
$\theta_{t+1}\leftarrow\mathsf{HalfRadiusAccel}(\theta_{t},R_{t},\mathcal{O}_{\textup{ng}},\frac{\delta}{T})$

6: $R_{t+1}\leftarrow\frac{1}{2}R_{t}$

7: $t\leftarrow t+1$

8: end while