# Robust Second-Order Nonconvex Optimization and
Its Application to Low Rank Matrix Sensing

Shuyao Li
University of Wisconsin-Madison
<EMAIL_ADDRESS>&Yu Cheng
Brown University
<EMAIL_ADDRESS>
Ilias Diakonikolas
University of Wisconsin-Madison
<EMAIL_ADDRESS>&Jelena Diakonikolas
University of Wisconsin-Madison
<EMAIL_ADDRESS>
Rong Ge
Duke University
<EMAIL_ADDRESS>&Stephen Wright
University of Wisconsin-Madison
<EMAIL_ADDRESS>

###### Abstract

Finding an approximate second-order stationary point (SOSP) is a well-studied
and fundamental problem in stochastic nonconvex optimization with many
applications in machine learning. However, this problem is poorly understood
in the presence of outliers, limiting the use of existing nonconvex algorithms
in adversarial settings.

In this paper, we study the problem of finding SOSPs in the strong
contamination model, where a constant fraction of datapoints are arbitrarily
corrupted. We introduce a general framework for efficiently finding an
approximate SOSP with _dimension-independent_ accuracy guarantees, using
$\widetilde{O}({D^{2}}/{\epsilon})$ samples where $D$ is the ambient dimension
and $\epsilon$ is the fraction of corrupted datapoints.

As a concrete application of our framework, we apply it to the problem of low
rank matrix sensing, developing efficient and provably robust algorithms that
can tolerate corruptions in both the sensing matrices and the measurements. In
addition, we establish a Statistical Query lower bound providing evidence that
the quadratic dependence on $D$ in the sample complexity is necessary for
computationally efficient algorithms.

##  1 Introduction

Learning in the presence of corrupted data is a significant challenge in
machine learning (ML) with many applications, including ML security [Bar+10,
BNL12, SKL17, Dia+19] and exploratory data analysis of real datasets, e.g., in
biological settings [Ros+02, Pas+10, Li+08, Dia+17]. The goal in such
scenarios is to design efficient learning algorithms that can tolerate a small
constant fraction of outliers and achieve error guarantees independent of the
dimensionality of the data. Early work in robust statistics [Ham+86, HR09]
gave sample-efficient robust estimators for various tasks (e.g., the Tukey
median [Tuk75] for robust mean estimation), alas with runtimes exponential in
the dimension. A recent line of work in computer science, starting with
[Dia+16, LRV16], developed the first computationally efficient robust
algorithms for several fundamental high-dimensional tasks. Since these early
works, there has been significant progress in algorithmic aspects of robust
high-dimensional statistics (see [DK19] and [DK23] for a comprehensive
overview).

In this paper, we study the general problem of smooth (with Lipschitz gradient
and Hessian) stochastic nonconvex optimization $\min_{x}\bar{f}(x)$ in the
outlier-robust setting, where
$\bar{f}(x):=\operatorname*{\mathbb{E}}_{A\sim\mathcal{G}}f(x,A)$ and
$\mathcal{G}$ is a possibly unknown distribution of the random parameter $A$.
We will focus on the following standard adversarial contamination model (see,
e.g., [Dia+16]).

######  Definition 1.1 (Strong Contamination Model).

Given a parameter $0<\epsilon<1/2$ and an inlier distribution $\mathcal{G}$,
an algorithm receives samples from $\mathcal{G}$ with $\epsilon$-contamination
as follows: The algorithm first specifies the number of samples $n$, and $n$
samples are drawn independently from $\mathcal{G}$. An adversary is then
allowed to inspect these samples, and replace an $\epsilon$-fraction of the
samples with arbitrary points. This modified set of $n$ points is said to be
$\epsilon$-corrupted, which is then given to the algorithm.

The stochastic optimization problem we consider is computationally intractable
in full generality — even without corruption — if the goal is to obtain
globally optimal solutions. At a high level, an achievable goal is to design
sample and computationally efficient robust algorithms for finding locally
optimal solutions. Prior work [Pra+20, Dia+19] studied outlier-robust
stochastic optimization and obtained efficient algorithms for finding
approximate first-order stationary points. While first-order guarantees
suffice for convex problems, it is known that in many tractable non-convex
problems, first-order stationary points may be bad solutions, but all second-
order stationary points (SOSPs) are globally optimal. This motivates us to
study the following questions:

> Can we develop a general framework for finding second-order stationary
> points in outlier-robust stochastic optimization?
>
> Can we obtain sample and computationally efficient algorithms for outlier-
> robust versions of tractable nonconvex problems using this framework?

In this work, we answer both questions affirmatively. We introduce a framework
for efficiently finding an approximate SOSP when $\epsilon$-fraction of the
functions are corrupted and then use our framework to solve the problem of
outlier-robust low rank matrix sensing.

In addition to the gradient being zero, a SOSP requires the Hessian matrix to
not have negative eigenvalues. The second-order optimality condition is
important because it rules out suboptimal solutions such as strict saddle
points. It is known that all SOSPs are globally optimal in nonconvex
formulations of many important machine learning problems, such as matrix
completion [GLM16], matrix sensing [BNS16], phase retrieval [SQW16], phase
synchronization [BBV16], dictionary learning [SQW17], and tensor decomposition
[Ge+15] (see also [WM22, Chapter 7]). However, the properties of SOSPs are
highly sensitive to perturbation in the input data. For example, it is
possible to create spurious SOSPs for nonconvex formulations of low rank
matrix recovery problems, even for a semi-random adversary that can add
additional sensing matrices but cannot corrupt the measurements in matrix
sensing [GC23] or an adversary who can only reveal more entries of the ground-
truth matrix in matrix completion [CG18]. Those spurious SOSPs correspond to
highly suboptimal solutions.

Finding SOSPs in stochastic nonconvex optimization problems in the presence of
arbitrary outliers was largely unaddressed prior to our work. Prior works
[Pra+20, Dia+19] obtained efficient and robust algorithms for finding first-
order stationary points with dimension-independent accuracy guarantees. These
works relied on the following simple idea: Under certain smoothness
assumptions, projected gradient descent with an approximate gradient oracle
efficiently converges to an approximate first-order stationary point.
Moreover, in the outlier-robust setting, approximating the gradient at a
specific point amounts to a robust mean estimation problem (for the underlying
distribution of the gradients), which can be solved by leveraging existing
algorithms for robust mean estimation.

Our work is the first to find approximate SOSPs with dimension-independent
errors in outlier-robust settings. Note that in standard non-robust settings,
approximate SOSPs can be computed using first-order methods such as perturbed
gradient descent [Jin+17, Jin+21]. This strategy might seem extendable to
outlier-robust settings through perturbed approximate gradient descent,
utilizing robust mean estimation algorithms to approximate gradients. The
approach in [Yin+19] follows this idea, but unfortunately their second-order
guarantees scale polynomially with dimension, even under very strong
distributional assumptions (e.g., subgaussianity). Our lower bound result
provides evidence that approximating SOSPs with dimension-independent error is
as hard as approximating _full_ Hessian, suggesting that solely approximating
the gradients is not sufficient. On a different note, [IPL23] recently
employed robust estimators for both gradient and Hessian in solving certain
convex stochastic optimization problems, which has a different focus than ours
and does not provide SOSPs with the guarantees that we achieve.

###  1.1 Our Results and Contributions

The notation we use in this section is defined in Section 2. To state our
results, we first formally define our generic nonconvex optimization problem.
Suppose there is a true distribution over functions
$f:\mathbb{R}^{D}\times\mathcal{A}\rightarrow\mathbb{R}$, where $f(x,A)$ takes
an argument $x\in\mathbb{R}^{D}$ and is parameterized by a random variable
$A\in\mathcal{A}$ drawn from a distribution $\mathcal{G}$. Our goal is to find
an $(\epsilon_{g},\epsilon_{H})$-approximate SOSP of the function
$\bar{f}(x):=\operatorname*{\mathbb{E}}_{A\sim\mathcal{G}}f(x,A)$.

######  Definition 1.2 ($\epsilon$-Corrupted Stochastic Optimization).

The algorithm has access to $n$ functions $(f_{i})_{i=1}^{n}$ generated as
follows. First $n$ random variables $(A_{i})_{i=1}^{n}$ are drawn
independently from $\mathcal{G}$. Then an adversary arbitrarily corrupts an
$\epsilon$ fraction of the $A_{i}$’s. Finally, the $\epsilon$-corrupted
version of $f_{i}(\cdot)=f(\cdot,A_{i})$ is sent to the algorithm as input.
The task is to find an approximate SOSP of the ground-truth average function
$\bar{f}(\cdot):=\operatorname*{\mathbb{E}}_{A\sim\mathcal{G}}f(\cdot,A)$.

######  Definition 1.3 (Approximate SOSPs).

A point $x$ is an $(\epsilon_{g},\epsilon_{H})$-approximate second-order
stationary point (SOSP) of $\bar{f}$ if
$\left\|\nabla\bar{f}(x)\right\|\leq\epsilon_{g}$ and
${\lambda_{\min}\left(\nabla^{2}\bar{f}(x)\right)\geq-\epsilon_{H}}$.

We make the following additional assumptions on $f$ and $\mathcal{G}$.

######  Assumption 1.4.

There exists a bounded region $\mathcal{B}$ such that the following conditions
hold:

  1. (i)

There exists a lower bound $f^{*}>-\infty$ such that for all
$x\in\mathcal{B}$, $f(x,A)\geq f^{*}$ with probability $1$.

  2. (ii)

There exist parameters $L_{D_{g}}$, $L_{D_{H}}$, $B_{D_{g}}$, and $B_{D_{H}}$
such that, with high probability over the randomness in $A\sim\mathcal{G}$,
letting $g(x)=f(x,A)$, we have that $g(x)$ is $L_{D_{g}}$-gradient Lipschitz
and $L_{D_{H}}$-Hessian Lipschitz over $\mathcal{B}$, and $\left\|\nabla
g(x)\right\|\leq B_{D_{g}}$ and $\left\|\nabla^{2}g(x)\right\|_{F}\leq
B_{D_{H}}$ for all $x\in\mathcal{B}$.

  3. (iii)

There exist parameters $\sigma_{g},\sigma_{H}>0$ such that for all
$x\in\mathcal{B}$,

$\left\|\operatorname{Cov}_{A\sim\mathcal{G}}(\nabla
f(x,A))\right\|_{\operatorname{op}}\leq\sigma_{g}^{2}$ and
$\left\|\operatorname{Cov}_{A\sim\mathcal{G}}(\operatorname{vec}(\nabla^{2}f(x,A)))\right\|_{\operatorname{op}}\leq\sigma_{H}^{2}$.

Note that the radius of $\mathcal{B}$ and the parameters $L_{D_{g}}$,
$L_{D_{H}}$, $B_{D_{g}}$, $B_{D_{H}}$ are all allowed to depend polynomially
on $D$ and $\epsilon$ (but not on $x$ and $A$).

Our main algorithmic result for $\epsilon$-corrupted stochastic optimization
is summarized in the following theorem. A formal version of this theorem is
stated as Theorem 3.1 in Section 3.

######  Theorem 1.5 (Finding an Outlier-Robust SOSP, informal).

Suppose $f$ satisfies Assumption 1.4 in a region $\mathcal{B}$ with parameters
$\sigma_{g}$ and $\sigma_{H}$. Given an arbitrary initial point
$x_{0}\in\mathcal{B}$ and an $\epsilon$-corrupted set of
$n=\widetilde{\Omega}\\!\left(D^{2}/\epsilon\right)$ functions where $D$ is
the ambient dimension, there exists a polynomial-time algorithm that with high
probability outputs an
$(O\\!\left(\sigma_{g}\sqrt{\epsilon}\right),O\\!\left(\sigma_{H}\sqrt{\epsilon}\right))$-approximate
SOSP of $\bar{f}$, provided that all iterates of the algorithm stay inside
$\mathcal{B}$.

Although the bounded iterate condition in Theorem 1.5 appears restrictive,
this assumption holds if the objective function satisfies a “dissipativity”
property, which is a fairly general phenomenon [Hal10]. Moreover, adding an
$\ell_{2}$-regularization term enables any Lipschitz function to satisfy the
dissipativity property [RRT17, Section 4]. As an illustrating example, a
simple problem-specific analysis shows that this bounded iterate condition
holds for outlier-robust matrix sensing by exploiting the fact that the matrix
sensing objective satisfies the dissipativity property.

In this paper, we consider the problem of outlier-robust symmetric low rank
matrix sensing, which we formally define below. We focus on the setting with
Gaussian design.

######  Definition 1.6 (Outlier-Robust Matrix Sensing).

There is an unknown rank-$r$ ground-truth matrix $M^{*}\in\mathbb{R}^{d\times
d}$ that can be factored into $U^{*}{U^{*}}^{\top}$ where
$U^{*}\in\mathbb{R}^{d\times r}$. The (clean) sensing matrices
$\\{A_{i}\\}_{i\in[n]}$ have i.i.d. standard Gaussian entries. The (clean)
measurements $y_{i}$ are obtained as $y_{i}=\langle
A_{i},M^{*}\rangle+\zeta_{i}$, where the noise
$\zeta_{i}\sim\mathcal{N}(0,\sigma^{2})$ is independent from all other
randomness. We denote the (clean) data generation process by
$(A_{i},y_{i})\sim\mathcal{G}_{\sigma}$. When $\sigma=0$, we have
$\zeta_{i}=0$ and we write $\mathcal{G}:=\mathcal{G}_{0}$ for this noiseless
(measurement) setting. In outlier robust matrix sensing, an adversary can
arbitrarily change any $\epsilon$-fraction of the sensing matrices and the
corresponding measurements. This corrupted set of $(A_{i},y_{i})$’s is then
given to the algorithm as input, where the goal is to recover $M^{*}$.

We highlight that in our setting, both the sensing matrices
$A_{i}\in\mathbb{R}^{d\times d}$ and the measurements $y_{i}\in\mathbb{R}$ can
be corrupted, presenting a substantially more challenging problem compared to
prior works (e.g., [Li+20a, Li+20]) that only allow corruption in $y_{i}$.

Let ${\sigma^{\star}_{1}}$ and ${\sigma^{\star}_{r}}$ denote the largest and
the smallest nonzero singular value of $M^{*}$ respectively. We assume
${\sigma^{\star}_{r}}$ and the rank $r$ are given to the algorithm, and we
assume that the algorithm knows a multiplicative upper bound $\Gamma$ of
${\sigma^{\star}_{1}}$ such that $\Gamma\geq 36{\sigma^{\star}_{1}}$ (a
standard assumption in matrix sensing even for non-robust settings [GJZ17,
Jin+17]). Let $\kappa=\Gamma/{\sigma^{\star}_{r}}$.

Our main algorithmic result for the low rank matrix sensing problem is
summarized in the following theorem. For a more detailed statement, see
Theorems 3.2 and 3.3 in Section 3.

######  Theorem 1.7 (Our Algorithm for Outlier-Robust Matrix Sensing).

Let $M^{*}\in\mathbb{R}^{d\times d}$ be the rank $r$ ground-truth matrix with
smallest nonzero singular value ${\sigma^{\star}_{r}}$. Let $\Gamma\geq
36\left\|M^{*}\right\|_{\operatorname{op}}$ and let
$\kappa=\Gamma/{\sigma^{\star}_{r}}$. There exists an algorithm for outlier-
robust matrix sensing, where an
$\epsilon=O\big{(}1/{(\kappa^{3}r^{3})}\big{)}$ fraction of samples from
$\mathcal{G}_{\sigma}$ as in Definition 1.6 gets arbitrarily corrupted, that
can output a rank-$r$ matrix $\widehat{M}$ such that
$\big{\|}\widehat{M}-M^{*}\big{\|}_{F}\leq\iota$ with probability at least
$1-\xi$, where $\iota>0$ is the error parameter:

  1. 1)

If $\sigma\geq r\Gamma$, then $\iota=O\\!\left(\sigma\sqrt{\epsilon}\right)$;

  2. 2)

If $\sigma\leq r\Gamma$, then
$\iota=O\\!\left(\kappa\sigma\sqrt{\epsilon}\right)$;

  3. 3)

If $\sigma=0$ (noiseless), then $\iota$ can be made arbitrarily small,
achieving exact recovery.

The algorithm uses
$n=\widetilde{O}\Big{(}\frac{d^{2}r^{2}+dr\log(\Gamma/\xi)}{\epsilon}\Big{)}$
samples and runs in time
$\operatorname{poly}(n,\kappa,\log({\sigma^{\star}_{r}}/\iota))$.

Finally, we complement our algorithmic results for outlier-robust matrix
sensing with a Statistical Query (SQ) lower bound, which provides strong
evidence that quadratic dependence on $d$ in the sample complexity is
unavoidable for efficient algorithms. A detailed statement of this result is
provided in Section 4.

###  1.2 Our Techniques

##### Outlier-robust nonconvex optimization.

To obtain our algorithmic result in the general nonconvex setting, we leverage
existing results on robust mean estimation [DKP20], which we use as a black
box to robustly estimate the gradient and the (vectorized) Hessian. We use
these robust estimates as a subroutine in a randomized nonconvex optimization
algorithm (described in Appendix A.2), which can tolerate inexactness in both
the gradient and the Hessian. With high probability, this algorithm outputs an
$(\epsilon_{g},\epsilon_{H})$-approximate SOSP, where $\epsilon_{g}$ and
$\epsilon_{H}$ depend on the inexactness of the gradient and Hessian oracles.

We remark that robust estimation of the Hessian is crucial to obtaining our
_dimension-independent_ approximation error result and is what causes $D^{2}$
dependence in the sample complexity (which is unavoidable for SQ algorithms as
discussed below). Notably, the only prior work on approximating SOSPs in the
outlier-robust setting [Yin+19] used robust mean estimation only on the
gradients and had sample complexity scaling linearly with $D$; however, they
can only output an order $(\sqrt{\epsilon},(\epsilon D)^{1/5})$-SOSP, which is
uninformative for many problems of interest, including the matrix sensing
problem considered in this paper, due to the dimensional dependence in the
approximation error.

##### Application to low rank matrix sensing.

Our main contribution on the algorithmic side is showing that our outlier-
robust nonconvex optimization framework can be applied to solve outlier-robust
matrix sensing with _dimension-independent_ approximation error, even
achieving exact recovery when the measurements are noiseless. We obtain this
result using the following geometric insights about the problem: We show that
the norm of the covariance of the gradient and the Hessian can both be upper
bounded by the sum of $\sigma^{2}$ and a function of the distance to the
closest optimal solution. We further prove that all iterates stay inside a
nice region using a “dissipativity” property [Hal10]), which says that the
iterate aligns with the direction of the gradient when the iterate’s norm is
large. This allows us to invoke Theorem 1.5 to obtain an approximate SOSP of
the ground-truth objective function.

We show that this approximate SOSP must be close to a global optimal solution.
Additionally, we establish a local regularity condition in a small region
around globally optimal solutions (which is similar to strong convexity but
holds only locally). This local regularity condition bounds below a measure of
stationarity, which allows us to prove that gradient descent-type updates
contract the distance to the closest global optimum under appropriate
stepsize. We leverage this local regularity condition to prove that the
iterates of the algorithm stay near a global optimum, so that the regularity
condition continues to hold, and moreover, the distance between the current
solution and the closest global optimum contracts, as long as it is larger
than a function of $\sigma$. Consequently, the distance-dependent component of
the gradient and Hessian covariance bound contracts as well, which allows us
to obtain more accurate gradient and Hessian estimates. While such a statement
may seem evident to readers familiar with linear convergence arguments, we
note that proving it is quite challenging, due to the circular dependence
between the distance from the current solution to global optima, the
inexactness in the gradient and Hessian estimates, and the progress made by
our algorithm.

The described distance-contracting argument allows us to control the
covariance of the gradient and Hessian, which we utilize to recover $M^{*}$
exactly when $\sigma=0$, and recover $M^{*}$ with error roughly
$O(\sigma\sqrt{\epsilon})$ when $0\neq\sigma\leq r\Gamma$. We note that the
$\sigma\sqrt{\epsilon}$ appears unavoidable in the $\sigma\neq 0$ case, due to
known limits of robust mean estimation algorithms [DKS17].

##### SQ lower bound.

We exhibit a hard instance of low rank matrix sensing problem to show that
quadratic dependence on the dimension in sample complexity is unavoidable for
computationally efficient algorithms. Our SQ lower bound proceeds by
constructing a family of distributions, corresponding to corruptions of low
rank matrix sensing, that are nearly uncorrelated in a well-defined technical
sense [Fel+17]. To achieve this, we follow the framework of [DKS17] which
considered a family of distributions that are rotations of a carefully
constructed one-dimensional distribution. The proof builds on [DKS19, Dia+21],
using a new univariate moment-matching construction which yields a family of
corrupted conditional distributions. These induce a family of joint
distributions that are SQ-hard to learn.

###  1.3 Roadmap

Section 2 defines the necessary notation and discusses relevant building
blocks of our algorithm and analysis. Section 3 introduces our framework for
finding SOSPs in the outlier-robust setting. Section 3.1 presents how to
extend and apply our framework to solve outlier-robust low rank matrix
sensing. Section 4 proves that our sample complexity has optimal dimensional
dependence for SQ algorithms. Most proofs are deferred to the supplementary
material due to space limitations.

##  2 Preliminaries

For an integer $n$, we use $[n]$ to denote the ordered set
$\\{1,2,\dots,n\\}$. We use $[a_{i}]_{i\in\mathcal{I}}$ to denote the matrix
whose columns are vectors $a_{i}$, where $\mathcal{I}$ is an ordered set. We
use $\mathbbm{1}_{E}(x)$ to denote the indicator function that is equal to $1$
if $x\in E$ and $0$ otherwise. For two functions $f$ and $g$, we say
$f=\widetilde{O}\\!\left(g\right)$ if $f=O\big{(}g\log^{k}(g)\big{)}$ for some
constant $k$, and we similarly define $\widetilde{\Omega}$.

For vectors $x$ and $y$, we let $\langle x,y\rangle$ denote the inner product
$x^{\top}y$ and $\left\|x\right\|$ denote the $\ell_{2}$ norm of $x$. For
$d\in\mathbb{Z}_{+}$, we use $I_{d}$ to denote the identity matrix of size
$d\times d$. For matrices $A$ and $B$, we use
$\left\|A\right\|_{\operatorname{op}}$ and $\left\|A\right\|_{F}$ to denote
the spectral norm and Frobenius norm of $A$ respectively. We use
$\lambda_{\max}(A)$ and $\lambda_{\min}(A)$ to denote the maximum and minimum
eigenvalue of $A$ respectively. We use $\operatorname{tr}(A)$ to denote the
trace of a matrix $A$. We use $\langle
A,B\rangle=\operatorname{tr}(A^{\top}B)$ to denote the entry-wise inner
product of two matrices of the same dimension. We use
$\operatorname{vec}(A)=[a_{1}^{\top},a_{2}^{\top},\dots,a_{d}^{\top}]^{\top}$
to denote the canonical flattening of $A$ into a vector, where
$a_{1},a_{2},\dots,a_{d}$ are columns of $A$.

######  Definition 2.1 (Lipschitz Continuity).

Let $\mathcal{X}$ and $\mathcal{Y}$ be normed vector spaces. A function
$h:\mathcal{X}\rightarrow\mathcal{Y}$ is $\ell$-Lipschitz if
$\left\|h(x_{1})-h(x_{2})\right\|_{\mathcal{Y}}\leq\ell\left\|x_{1}-x_{2}\right\|_{\mathcal{X}},$
$\forall x_{1},x_{2}.$

In this paper, when $\mathcal{Y}$ is a space of matrices, we take
$\left\|\cdot\right\|_{\mathcal{Y}}$ to be the spectral norm
$\left\|\cdot\right\|_{\operatorname{op}}$. When $\mathcal{X}$ is a space of
matrices, we take $\left\|\cdot\right\|_{\mathcal{X}}$ to be the Frobenius
norm $\left\|\cdot\right\|_{F}$; this essentially views the function $h$ as
operating on the vectorized matrices endowed with the usual $\ell_{2}$ norm.
When $\mathcal{X}$ or $\mathcal{Y}$ is the Euclidean space, we take the
corresponding norm to be the $\ell_{2}$ norm.

##### A Randomized Algorithm with Inexact Gradients and Hessians.

We now discuss how to solve the unconstrained nonconvex optimization problem
$\min_{x\in\mathbb{R}^{D}}f(x),$ where $f(\cdot)$ is a smooth function with
Lipschitz gradients and Lipschitz Hessians. The goal of this section is to
find an approximate SOSP as defined in Definition 1.3.

######  Proposition 2.2 ([LW23]).

Suppose a function $f$ is bounded below by $f^{*}>-\infty,$ has
$L_{g}$-Lipschitz gradient and $L_{H}$-Lipschitz Hessian, and its inexact
gradient and Hessian computations $\widetilde{g}_{t}$ and $\widetilde{H}_{t}$
satisfy
$\left\|\widetilde{g}_{t}-\nabla{f}(x_{t})\right\|\leq\frac{1}{3}\epsilon_{g}$
and
$\big{\|}\widetilde{H}_{t}-\nabla^{2}f(x_{t})\big{\|}_{\operatorname{op}}\leq\frac{2}{9}\epsilon_{H}$.
Then there exists an algorithm (Algorithm 2) with the following guarantees:

  1. 1.

(Correctness) If Algorithm 2 terminates and outputs $x_{n}$, then $x_{n}$ is a
$(\frac{4}{3}\epsilon_{g},\frac{4}{3}\epsilon_{H})$-approximate SOSP.

  2. 2.

(Runtime) Algorithm 2 terminates with probability $1$. Let
$C_{\epsilon}:=\min\left(\frac{\epsilon_{g}^{2}}{6L_{g}},\frac{2\epsilon_{H}^{3}}{9L_{H}^{2}}\right)$.
With probability at least $1-\delta$, Algorithm 2 terminates after $k$
iterations for

$k=O\Big{(}\frac{f(x_{0})-f^{*}}{C_{\epsilon}}+\frac{L_{H}^{2}L_{g}^{2}\epsilon_{g}^{2}}{\epsilon_{H}^{6}}\log\Bigl{(}\frac{1}{\delta}\Bigr{)}\Big{)}\;.$
(1)

The constants $1/3$ and $2/9$ are chosen for ease of presentation. For all
constructions of Hessian oracles in this paper, we take the straightforward
relaxation
$\big{\|}\widetilde{H}_{t}-\nabla^{2}f(x_{t})\big{\|}_{\operatorname{op}}\leq\big{\|}\widetilde{H}_{t}-\nabla^{2}f(x_{t})\big{\|}_{F}$
and upper bound Hessian inexactness using Frobenius norm. Proof of a
simplified version of Proposition 2.2 with a weaker high probability bound
that is sufficient for our purposes is provided in Appendix A.2 for
completeness.

##### Robust Mean Estimation.

Recent algorithms in robust mean estimation give dimension-independent error
in the presence of outliers under strong contamination model.

We use the following results, see, e.g., [DKP20], where the upper bound
$\sigma$ on the spectral norm of covariance matrix is unknown to the
algorithm.

######  Proposition 2.3 (Robust Mean Estimation).

Fix any $0<\xi<1$. Let $S$ be a multiset of $n=O((k\log
k+\log(1/\xi))/\epsilon)$ i.i.d. samples from a distribution on
$\mathbb{R}^{k}$ with mean $\mu_{S}$ and covariance $\Sigma$. Let
$T\subset\mathbb{R}^{k}$ be an $\epsilon$-corrupted version of $S$ as in
Definition 1.1. There exists an algorithm (Algorithm 3) such that, with
probability at least $1-\xi$, on input $\epsilon$ and $T$ (but not
$\left\|\Sigma\right\|_{\operatorname{op}}$) returns a vector $\widehat{\mu}$
in polynomial time so that
$\left\|\mu_{S}-\widehat{\mu}\right\|=O\big{(}\sqrt{\left\|\Sigma\right\|_{\operatorname{op}}\epsilon}\big{)}$.

Algorithm 3 is given in Appendix A.3. Proposition 2.3 states that for
Algorithm 3 to succeed with high probability,
$\widetilde{O}\\!\left({k}/{\epsilon}\right)$ i.i.d. samples need to be drawn
from a $k$-dimensional distribution of bounded covariance. State of the art
algorithms for robust mean estimation can be implemented in near-linear time,
requiring only a logarithmic number of passes on the data, see, e.g, [Dia+22,
CDG19, DHL19]. Any of these faster algorithms could be used for our purposes.

With the above results, the remaining technical component for applying the
robust estimation subroutine (Algorithm 3) in this paper is handling the
dependence across iterations.

Because we will run $\mathsf{RobustMeanEstimation}$ in each iteration of our
optimization algorithm, the gradients $\\{\nabla f_{i}(x_{t})\\}_{i=1}^{n}$
and Hessians $\\{\nabla^{2}f_{i}(x_{t})\\}_{i=1}^{n}$ can no longer be
considered as independently drawn from a distribution after the first
iteration. Although they are i.i.d. for fixed $x$, the dependence on previous
iterations through $x$ will break the independence assumption. Therefore, we
will need a union bound over all $x_{t}$ to handle dependence across different
iterations $t$. We deal with this technicality in Appendix B.

##  3 General Robust Nonconvex Optimization

In this section, we establish a general result that uses Algorithm 2 to obtain
approximate SOSPs in the presence of outliers under strong contamination. The
inexact gradient and inexact Hessian oracles are constructed with the robust
mean estimation subroutine (Algorithm 3). We consider stochastic optimization
tasks in Definition 1.2 satisfying Assumption 1.4. We construct the inexact
gradient and Hessian oracle required by Algorithm 2 as follows:

$\displaystyle\widetilde{g}_{t}$
$\displaystyle\leftarrow\mathbf{RobustMeanEstimation}(\\{\nabla
f_{i}(x_{t})\\}_{i=1}^{n},4\epsilon)$ $\displaystyle\widetilde{H}_{t}$
$\displaystyle\leftarrow\mathbf{RobustMeanEstimation}(\\{\nabla^{2}f_{i}(x_{t})\\}_{i=1}^{n},4\epsilon)$

Then we have the following guarantee:

######  Theorem 3.1.

Suppose we are given $\epsilon$-corrupted set of functions
$\\{f_{i}\\}_{i=1}^{n}$ for sample size $n$, generated according to Definition
1.2. Suppose Assumption 1.4 holds in a bounded region
$\mathcal{B}\subset\mathbb{R}^{D}$ of diameter $\gamma$ with gradient and
Hessian covariance bound $\sigma_{g}$ and $\sigma_{H}$ respectively, and we
have an arbitrary initialization $x_{0}\in\mathcal{B}$. Algorithm 2
initialized at $x_{0}$ outputs an $(\epsilon_{g},\epsilon_{H})$-approximate
SOSP for a sufficiently large sample with probability at least $1-\xi$ if the
following conditions hold:

  1. (I)

All iterates $x_{t}$ in Algorithm 2 stay inside the bounded region
$\mathcal{B}$.

  2. (II)

For an absolute constant $c>0$, it holds that $\sigma_{g}\sqrt{\epsilon}\leq
c\epsilon_{g}$ and $\sigma_{H}\sqrt{\epsilon}\leq c\epsilon_{H}$.

The algorithm uses $n=\widetilde{O}\big{(}D^{2}/\epsilon\big{)}$ samples,
where $\widetilde{O}\\!\left(\cdot\right)$ hides logarithmic dependence on
$D,\epsilon,L_{D_{g}},L_{D_{H}},B_{D_{g}},B_{D_{H}},\gamma/\sigma_{H},\gamma/\sigma_{g},$
and $1/\xi$. The algorithm runs in time polynomial in the above parameters.

Note that we are able to obtain dimension-independent errors $\epsilon_{g}$
and $\epsilon_{H}$, provided that $\sigma_{g}$ and $\sigma_{H}$ are dimension-
independent.

###  3.1 Low Rank Matrix Sensing Problems

In this section, we study the problem of outlier-robust low rank matrix
sensing as formally defined in Definition 1.6. We first apply the above
framework to obtain an approximate SOSP in Section 3.1.2. Then we make use of
the approximate SOSP to obtain a solution that is close to the ground-truth
matrix $M^{*}$ in Section 3.1.3; this demonstrates the usefulness of
approximate SOSPs.

####  3.1.1 Main results for Robust Low Rank Matrix Sensing

The following are the main results that we obtain in this section:

######  Theorem 3.2 (Main Theorem Under Noiseless Measurements).

Consider the noiseless setting as in Theorem 1.7 with $\sigma=0$. For some
sample size
$n=\widetilde{O}\big{(}{(d^{2}r^{2}+dr\log(\Gamma/\xi))}/{\epsilon}\big{)}$
and with probability at least $1-\xi$, there exists an algorithm that outputs
a solution that is $\iota$-close to $M^{*}$ in Frobenius norm in
$O\big{(}r^{2}\kappa^{3}\log(1/\xi)+\kappa\log({\sigma^{\star}_{r}}/\iota)\big{)}$
calls to the robust mean estimation subroutine (Algorithm 3).

This result achieves exact recovery of $M^{*}$, despite the strong
contamination of samples. Each iteration involves a subroutine call to robust
mean estimation. Algorithm 3 presented here is one simple example of robust
mean estimation; there are refinements [Dia+22, DHL19] that run in nearly
linear time, so the total computation utilizing those more efficient
algorithms indeed requires $\widetilde{O}\left(r^{2}\kappa^{3}\right)$ passes
of data (computed gradients and Hessians).

######  Theorem 3.3 (Main Theorem Under Noisy Measurements).

Consider the same setting as in Theorem 1.7 with $\sigma\neq 0$. There exists
a sample size
$n=\widetilde{O}\big{(}{(d^{2}r^{2}+dr\log(\Gamma/\xi))}/{\epsilon}\big{)}$
such that

  * •

if $\sigma\leq r\Gamma$, then with probability at least $1-\xi$, there exists
an algorithm that outputs a solution $\widehat{M}$ in
$\widetilde{O}(r^{2}\kappa^{3})$ calls to robust mean estimation routine 3,
with error
$\big{\|}\widehat{M}-M^{*}\big{\|}_{F}=O\\!\left(\kappa\sigma\sqrt{\epsilon}\right)$;

  * •

if $\sigma\geq r\Gamma$, then with probability at least $1-\xi$, there exists
a (different) algorithm that outputs a solution $\widehat{M}$ in one call to
robust mean estimation routine 3, with error
$\big{\|}\widehat{M}-M^{*}\big{\|}_{F}=O\\!\left(\sigma\sqrt{\epsilon}\right)$.

We prove Theorem 3.3 in Appendix C.4, and instead focus on the noiseless
measurements with $\sigma=0$ when we develop our algorithms in this section;
the two share many common techniques. In the remaining part of Section 3.1, we
use $\mathcal{G}_{0}$ in Definition 1.6 for the data generation process.

We now describe how we obtain the solution via nonconvex optimization.
Consider the following objective function for (uncorrupted) matrix sensing:

$\min_{\begin{subarray}{c}M\in\mathbb{R}^{d\times d}\\\
\operatorname{rank}(M)=r\end{subarray}}\frac{1}{2}\operatorname*{\mathbb{E}}_{(A_{i},y_{i})\sim\mathcal{G}_{0}}(\langle
M,A_{i}\rangle-y_{i})^{2}.$ (2)

We can write $M=UU^{\top}$ for some $U\in\mathbb{R}^{d\times r}$ to
reparameterize the objective function. Let

$f_{i}(U):=\frac{1}{2}\left(\langle UU^{\top},A_{i}\rangle-y_{i}\right)^{2}.$
(3)

We can compute

$\bar{f}(U):=\operatorname*{\mathbb{E}}_{(A_{i},y_{i})\sim\mathcal{G}_{0}}f_{i}(U)=\frac{1}{2}\operatorname{Var}\langle
UU^{\top}-M^{*},A_{i}\rangle=\frac{1}{2}\left\|UU^{\top}-M^{*}\right\|_{F}^{2}.$
(4)

We seek to solve the following optimization problem under the corruption model
in Definition 1.2:

$\min_{U\in\mathbb{R}^{d\times r}}\bar{f}(U).$ (5)

The gradient Lipschitz constant and Hessian Lipschitz constant of $\bar{f}$
are given by the following result.

######  Fact 3.4 ([Jin+17], Lemma 6).

For any $\Gamma>{\sigma^{\star}_{1}}$, $\bar{f}(U)$ has gradient Lipschitz
constant $L_{g}=16\Gamma$ and Hessian Lipschitz constant
$L_{H}=24\Gamma^{\frac{1}{2}}$ inside the region
$\\{U:\left\|U\right\|_{\operatorname{op}}^{2}<\Gamma\\}$.

####  3.1.2 Global Convergence to an Approximate SOSP

In this section, we apply our framework Theorem 3.1 to obtain global
convergence from an arbitrary initialization to an approximate SOSP, by
providing problem-specific analysis to guarantee that both Assumption 1.4 and
algorithmic assumptions (I) and (II) required by Theorem 3.1 are satisfied.

######  Theorem 3.5 (Global Convergence to a SOSP).

Consider the noiseless setting as in Theorem 1.7 with $\sigma=0$ and
$\epsilon=O\big{(}1/{(\kappa^{3}r^{3})}\big{)}$. Assume we have an arbitrary
initialization $U_{0}$ inside
$\\{U:\left\|U\right\|_{\operatorname{op}}^{2}\leq\Gamma\\}$. There exists a
sample size
$n=\widetilde{O}\left({(d^{2}r^{2}+dr\log(\Gamma/\xi))}/{\epsilon}\right)$
such that with probability at least $1-\xi$, Algorithm 2 initialized at
$U_{0}$ outputs a
$(\frac{1}{24}{\sigma^{\star}_{r}}^{3/2},\frac{1}{3}{\sigma^{\star}_{r}})$-approximate
SOSP using at most $O\\!\left(r^{2}\kappa^{3}\log(1/\xi)\right)$ calls to
robust mean estimation subroutine (Algorithm 3).

###### Proof of Theorem 3.5.

To apply Theorem 3.1, we verify Assumption 1.4 first. To verify (i), for all
$U$ and $A_{i}$, $f_{i}(U)=\frac{1}{2}\left(\langle UU^{\top},A_{i}\rangle-
y_{i}\right)^{2}\geq 0$ , so $f^{*}=0$ is a uniform lower bound. We verify
(ii) in Appendix C.2: conceptually, by Fact 3.4, $\bar{f}$ is gradient and
Hessian Lipschitz; both gradient and Hessian of $f_{i}$ are sub-exponential
and concentrate around those of $\bar{f}$. To check (iii), we calculate the
gradients and Hessians of $f_{i}$ in Appendix C.1.1 and bound their
covariances from above in Appendix C.1.2 and C.1.3. The result is summarized
in the following lemma. Note that the domain of the target function in
Algorithm 2 and Theorem 3.1 is the Euclidean space $\mathbb{R}^{D}$, so we
vectorize $U$ and let $D=dr$. The gradient becomes a vector in
$\mathbb{R}^{dr}$ and the Hessian becomes a matrix in $\mathbb{R}^{dr\times
dr}$.

######  Lemma 3.6 (Gradient and Hessian Covariance Bounds).

For all $U\in\mathbb{R}^{d\times r}$ with
$\left\|U\right\|_{\operatorname{op}}^{2}\leq\Gamma$ and $f_{i}$ defined in
Equation (3), it holds

$\displaystyle\left\|\operatorname{Cov}(\operatorname{vec}(\nabla
f_{i}(U)))\right\|_{\operatorname{op}}$ $\displaystyle\leq
8\left\|UU^{\top}-M^{*}\right\|_{F}^{2}\left\|U\right\|_{\operatorname{op}}^{2}\leq
32r^{2}\Gamma^{3}$ (6)
$\displaystyle\left\|\operatorname{Cov}(\operatorname{vec}(H_{i}))\right\|_{\operatorname{op}}$
$\displaystyle\leq
16r\left\|UU^{\top}-M^{*}\right\|_{F}^{2}+128\left\|U\right\|_{\operatorname{op}}^{4}\leq
192r^{3}\Gamma^{2}$ (7)

We proceed to verify the algorithmic assumptions in Theorem 3.1. For the
assumption (I), we prove the following Lemma in Appendix C.2 to show that all
iterates stay inside the bounded region in which we compute the covariance
bounds.

######  Lemma 3.7.

All iterates of Algorithm 2 stay inside the region
$\\{U:\left\|U\right\|_{\operatorname{op}}^{2}\leq\Gamma\\}$.

To verify Theorem 3.1 (II), we let
$\epsilon_{g}=\frac{1}{32}{\sigma^{\star}_{r}}^{3/2},\epsilon_{H}=\frac{1}{4}{\sigma^{\star}_{r}}$
and $\sigma_{g}=8r\Gamma^{1.5},\sigma_{H}=16r^{1.5}\Gamma$. So if we assume
$\epsilon=O(1/{(\kappa^{3}r^{3}}))$, then for the absolute constant $c$ in
Theorem 3.1 it holds that

$\displaystyle\sigma_{g}\sqrt{\epsilon}\leq
c{\epsilon_{g}}\qquad\sigma_{H}\sqrt{\epsilon}\leq c{\epsilon_{H}}.$

Hence, Theorem 3.1 applies and Algorithm 2 outputs an
$(\epsilon_{g},\epsilon_{H})$-approximate SOSP with high probability in
polynomial time. To bound the runtime, since
$\bar{f}(U_{0})=1/2\left\|U_{0}U_{0}^{\top}-M^{*}\right\|_{F}^{2}=O(r^{2}\Gamma^{2})$
for an arbitrary initialization $U_{0}$ with
$\left\|U_{0}\right\|_{\operatorname{op}}^{2}<\Gamma$, the initial distance
can be bounded by $O(r^{2}\Gamma^{2})$. Setting
$L_{g}=16\Gamma,L_{H}=24\Gamma^{1/2},\bar{f}(U_{0})=O(r^{2}\Gamma^{2}),f^{*}=0$
and thus $C_{\epsilon}=O\big{(}{{\sigma^{\star}_{r}}^{3}}/{\Gamma}\big{)}$,
Proposition 2.2 implies that Algorithm 2 outputs a
$(\frac{1}{24}{\sigma^{\star}_{r}}^{3/2},\frac{1}{3}{\sigma^{\star}_{r}})$-approximate
second order stationary point $U_{SOSP}$ in
$O\big{(}r^{2}\kappa^{3}\log(1/\xi)\big{)}$ steps with high probability. ∎

####  3.1.3 Local Linear Convergence

In this section, we describe a local search algorithm that takes a
$(\frac{1}{24}{\sigma^{\star}_{r}}^{3/2},\frac{1}{3}{\sigma^{\star}_{r}})$-approximate
second-order stationary point as its initialization and achieves exact
recovery even in the presence of outliers.

Data: The initialization $U_{SOSP}$ is a
$(\frac{1}{24}{\sigma^{\star}_{r}}^{3/2},\frac{1}{3}{\sigma^{\star}_{r}})$-approximate
SOSP, corruption fraction is $\epsilon$, corrupted samples are
$\\{(A_{i},y_{i})\\}_{i=1}^{n}$, target distance to optima is $\iota$

Result: $U$ that is $\iota$-close in Frobenius norm to some global minimum

1 $\eta=1/\Gamma$, $U_{0}=U_{SOSP}$

2 for  _t = 0, 1, …_ do

3 $\widetilde{g}_{t}:=\mathbf{RobustMeanEstimation}(\\{\nabla
f_{i}(U_{t})\\}_{i=1}^{n},4\epsilon)$

4 $U_{t+1}\leftarrow U_{t}-\eta\widetilde{g}_{t}$

Algorithm 1 Local Inexact Gradient Descent

######  Theorem 3.8 (Local Linear Convergence).

Consider the same noiseless setting as in Theorem 1.7. Assume we already found
a
$(\frac{1}{24}{\sigma^{\star}_{r}}^{3/2},\frac{1}{3}{\sigma^{\star}_{r}})$-approximate
SOSP $U_{SOSP}$ of $\bar{f}$. Then there exists a sample size
$n=\widetilde{O}\left({dr\log(1/\xi)}/{\epsilon}\right)$ such that with
probability at least $1-\xi$, Algorithm 1 initialized at $U_{SOSP}$ outputs a
solution that is $\iota$-close to some global minimum in Frobenius norm after
$O\big{(}\kappa\log({\sigma^{\star}_{r}}/{\iota})\big{)}$ calls to robust mean
estimation subroutine (Algorithm 3). Moreover, all iterates $U_{t}$ are
$\frac{1}{3}{\sigma^{\star}_{r}}^{1/2}$-close to some global minimum in
Frobenius norm.

###### Proof Sketch.

First we use known properties of
$\bar{f}=\operatorname*{\mathbb{E}}_{(A_{i},y_{i})\sim\mathcal{G}_{0}}f_{i}$
from the literature [Jin+17] to show approximate SOSPs of $\bar{f}$ — in
particular our initialization $U_{SOSP}$ — are in a small neighborhood of the
global minima of $\bar{f}$. In that neighborhood, it was also known that
$\bar{f}$ satisfies some local regularity conditions that enable gradient
descent’s linear convergence.

However, the algorithm only has access to the inexact gradient from the robust
mean estimation subroutine (line 1 in Algorithm 1) and, therefore, we need to
establish linear convergence of inexact gradient descent. We achieve this with
the following iterative argument: As the iterate gets closer to global optima,
the covariance bound of sample gradient in Equation (6) gets closer to 0.
Because the accuracy of the robust mean estimation scales with the covariance
bound (see Proposition 2.3), a more accurate estimate of population gradient
$\nabla\bar{f}$ can be obtained via the robust estimation subroutine
(Algorithm 3). This, in turn, allows for an improved inexact gradient descent
step, driving the algorithm towards an iterate that is even closer to global
optima.

See Appendix C.3 for the complete proof. ∎

##  4 Statistical Query Lower Bound for Low Rank Matrix Sensing

Our general algorithm (Theorem 3.1) leads to an efficient algorithm for robust
low rank matrix sensing with sample complexity $O(d^{2}r^{2}/\epsilon)$
(Theorem 3.3). Interestingly, the sample complexity of the underlying robust
estimation problem — ignoring computational considerations — is
$\Theta(dr/\epsilon)$. The information-theoretic upper bound of
$O(dr/\epsilon)$ can be achieved by an exponential (in the dimension) time
algorithm (generalizing the Tukey median to our regression setting); see,
e.g., Theorem 3.5 in [Gao20].

Given this discrepancy, it is natural to ask whether the sample complexity
achieved by our algorithm can be improved via a different computationally
efficient method. In this section, we provide evidence that this may not be
possible. In more detail, we establish a near-optimal information-computation
tradeoff for the problem, within the class of Statistical Query (SQ)
algorithms. To formally state our lower bound, we require basic background on
SQ algorithms.

##### Basics on SQ Model.

SQ algorithms are a class of algorithms that, instead of access to samples
from some distribution $\mathcal{P}$, are allowed to query expectations of
bounded functions over $\mathcal{P}$.

######  Definition 4.1 (SQ Algorithms and $\mathsf{STAT}$ Oracle [Kea98]).

Let $\mathcal{P}$ be a distribution on $\mathbb{R}^{d^{2}+1}$. A Statistical
Query (SQ) is a bounded function $q:\mathbb{R}^{d^{2}+1}\rightarrow[-1,1]$.
For $\tau>0$, the $\mathsf{STAT}(\tau)$ oracle responds to the query $q$ with
a value $v$ such that
$\left|{v-\operatorname*{\mathbb{E}}_{X\sim\mathcal{P}}[q(X)]}\right|\leq\tau$.
An SQ algorithm is an algorithm whose objective is to learn some information
about an unknown distribution $\mathcal{P}$ by making adaptive calls to the
corresponding $\mathsf{STAT}(\tau)$ oracle.

In this section, we consider $\mathcal{P}$ as the unknown corrupted
distribution where $(A_{i},y_{i})$ are drawn. The SQ algorithm tries to learn
the ground truth matrix $M^{*}$ from this corrupted distribution; the goal of
the lower bound result is to show that this is hard.

The SQ model has the capability to implement a diverse set of algorithmic
techniques in machine learning such as spectral techniques, moment and tensor
methods, local search (e.g., Expectation Maximization), and several others
[Fel+17]. A lower bound on the SQ complexity of a problem provides evidence of
hardness for the problem. [Bre+21] established that (under certain
assumptions) an SQ lower bound also implies a qualitatively similar lower
bound in the low-degree polynomial testing model. This connection can be used
to show a similar lower bound for low-degree polynomials.

Our main result here is a near-optimal SQ lower bound for robust low rank
matrix sensing that applies even for rank $r=1$, i.e., when the ground truth
matrix is $M^{*}=uu^{\top}$ for some $u\in\mathbb{R}^{d}$. The choice of rank
$r=1$ yields the strongest possible lower bound in our setting because it is
the easiest parameter regime: Recall that the sample complexity of our
algorithm is $\widetilde{O}(d^{2}r^{2})$ as in Theorems 3.2 and 3.3, and the
main message of our SQ lower bound is to provide evidence that the $d^{2}$
factor is necessary for computationally efficient algorithms _even if_ $r=1$.

######  Theorem 4.2 (SQ Lower Bound for Robust Rank-One Matrix Sensing).

Let $\epsilon\in(0,1/2)$ be the fraction of corruptions and let $c\in(0,1/2)$.
Assume the dimension $d\in\mathbb{N}$ is sufficiently large. Consider the
$\epsilon$-corrupted rank-one matrix sensing problem with ground-truth matrix
$M^{*}=uu^{\top}$ and noise $\sigma^{2}=O(1)$. Any SQ algorithm that outputs
$\widehat{u}$ with $\left\|\widehat{u}-u\right\|=O(\epsilon^{1/4})$ either
requires $2^{\Omega(d^{c})}/d^{2-4c}$ queries or makes at least one query to
$\mathsf{STAT}\Big{(}e^{O\\!\left(1/\sqrt{\epsilon}\right)}/O\\!\left(d^{1-2c}\right)\Big{)}$.

In other words, we show that, when provided with SQ access to an
$\epsilon$-corrupted distribution, approximating $u$ is impossible unless
employing a statistical query of higher precision than what can be achieved
with a strictly sub-quadratic number (e.g., $d^{1.99}$) of samples. Note that
the SQ oracle
$\mathsf{STAT}\bigl{(}e^{O\big{(}1/\sqrt{\epsilon}\big{)}}/Obig(d^{1-2c}big)\bigr{)}$
can be simulated with $O(d^{2-4c})/e^{O(1/\epsilon)}$ samples, and this bound
is tight in general. Informally speaking, this theorem implies that improving
the sample complexity from $d^{2}$ to $d^{2-4c}$ requires exponentially many
queries. This result can be viewed as a near-optimal information-computation
tradeoff for the problem, within the class of SQ algorithms.

The proof follows a similar analysis as in [DKS19, Dia+21], using one-
dimensional moment matching to construct a family of corrupted conditional
distributions, which induce a family of corrupted joint distributions that are
SQ-hard to learn. We provide the details of the proof in Appendix D. Apart
from the formal proof, in Appendix E we also informally discuss the intuition
for why some simple algorithms that require $O(d)$ samples do not provide
dimension-independent error guarantees.

## Acknowledgements

Shuyao Li was supported in part by NSF Awards DMS-2023239, NSF Award
CCF-2007757 and the U. S. Office of Naval Research under award number
N00014-22-1-2348. Yu Cheng was supported in part by NSF Award CCF-2307106.
Ilias Diakonikolas was supported in part by NSF Medium Award CCF-2107079, NSF
Award CCF-1652862 (CAREER), a Sloan Research Fellowship, and a DARPA Learning
with Less Labels (LwLL) grant. Jelena Diakonikolas was supported in part by
NSF Award CCF-2007757 and by the U. S. Office of Naval Research under award
number N00014-22-1-2348. Rong Ge was supported in part by NSF Award
DMS-2031849, CCF-1845171 (CAREER) and a Sloan Research Fellowship. Stephen
Wright was supported in part by NSF Awards DMS-2023239 and CCF-2224213 and
AFOSR via subcontract UTA20-001224 from UT-Austin.

## References

  * [Bar+10] M. Barreno, B. Nelson, A.. Joseph and J.. Tygar  “The security of machine learning”  In _Machine Learning_ 81.2, 2010, pp. 121–148
  * [BBV16] Afonso S Bandeira, Nicolas Boumal and Vladislav Voroninski  “On the low-rank approach for semidefinite programs arising in synchronization and community detection”  In _Conference on learning theory_ , 2016, pp. 361–382  PMLR
  * [Ber+22] El Houcine Bergou, Youssef Diouane, Vladimir Kunc, Vyacheslav Kungurtsev and Clément W. Royer  “A subsampling line-search method with second-order results”  In _INFORMS Journal on Optimization_ 4.4, 2022, pp. 403–425
  * [BNL12] Battista Biggio, Blaine Nelson and Pavel Laskov  “Poisoning Attacks against Support Vector Machines”  In _Proceedings of the 29th International Coference on International Conference on Machine Learning_ Omnipress, 2012, pp. 1467–1474
  * [BNS16] Srinadh Bhojanapalli, Behnam Neyshabur and Nati Srebro  “Global optimality of local search for low rank matrix recovery”  In _Advances in Neural Information Processing Systems_ , 2016, pp. 3873–3881
  * [Bre+21] Matthew S Brennan, Guy Bresler, Sam Hopkins, Jerry Li and Tselil Schramm  “Statistical Query Algorithms and Low Degree Tests Are Almost Equivalent”  In _Proceedings of Thirty Fourth Conference on Learning Theory_ 134, Proceedings of Machine Learning Research  PMLR, 2021, pp. 774–774
  * [CDG19] Yu Cheng, Ilias Diakonikolas and Rong Ge  “High-Dimensional Robust Mean Estimation in Nearly-Linear Time”  In _Proceedings of the 30th ACM-SIAM Symposium on Discrete Algorithms (SODA)_ SIAM, 2019, pp. 2755–2771
  * [CG18] Yu Cheng and Rong Ge  “Non-convex matrix completion against a semi-random adversary”  In _Conference On Learning Theory_ , 2018, pp. 1362–1394  PMLR
  * [DHL19] Yihe Dong, Samuel Hopkins and Jerry Li  “Quantum entropy scoring for fast robust mean estimation and improved outlier detection”  In _Advances in Neural Information Processing Systems_ 32, 2019
  * [Dia+16] Ilias Diakonikolas, Gautam Kamath, Daniel M. Kane, Jerry Li, Ankur Moitra and Alistair Stewart  “Robust estimators in high dimensions without the computational intractability”  In _57th Annual IEEE Symposium on Foundations of Computer Science—FOCS 2016_ IEEE Computer Soc., Los Alamitos, CA, 2016, pp. 655–664
  * [Dia+17] Ilias Diakonikolas, Gautam Kamath, Daniel M. Kane, Jerry Li, Ankur Moitra and Alistair Stewart  “Being Robust (in High Dimensions) Can Be Practical”  In _Proceedings of the 34th International Conference on Machine Learning_ 70, Proceedings of Machine Learning Research  PMLR, 2017, pp. 999–1008
  * [Dia+19] Ilias Diakonikolas, Gautam Kamath, Daniel Kane, Jerry Li, Jacob Steinhardt and Alistair Stewart  “Sever: A Robust Meta-Algorithm for Stochastic Optimization”  In _Proceedings of the 36th International Conference on Machine Learning_ 97, Proceedings of Machine Learning Research  PMLR, 2019, pp. 1596–1606
  * [Dia+21] Ilias Diakonikolas, Daniel Kane, Ankit Pensia, Thanasis Pittas and Alistair Stewart  “Statistical Query Lower Bounds for List-Decodable Linear Regression”  In _Advances in Neural Information Processing Systems_ 34 Curran Associates, Inc., 2021, pp. 3191–3204
  * [Dia+22] Ilias Diakonikolas, Daniel M. Kane, Ankit Pensia and Thanasis Pittas  “Streaming Algorithms for High-Dimensional Robust Statistics”  In _Proceedings of the 39th International Conference on Machine Learning_ 162, Proceedings of Machine Learning Research  PMLR, 2022, pp. 5061–5117
  * [DK19] Ilias Diakonikolas and Daniel M Kane  “Recent advances in algorithmic high-dimensional robust statistics”  In _arXiv preprint arXiv:1911.05911_ , 2019
  * [DK23] Ilias Diakonikolas and Daniel M. Kane  “Algorithmic High-Dimensional Robust Statistics”  Cambridge University Press, 2023
  * [DKP20] Ilias Diakonikolas, Daniel M. Kane and Ankit Pensia  “Outlier Robust Mean Estimation with Subgaussian Rates via Stability”  In _Advances in Neural Information Processing Systems_ 33 Curran Associates, Inc., 2020, pp. 1830–1840
  * [DKS17] Ilias Diakonikolas, Daniel M. Kane and Alistair Stewart  “Statistical Query Lower Bounds for Robust Estimation of High-Dimensional Gaussians and Gaussian Mixtures”  In _2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS)_ , 2017, pp. 73–84
  * [DKS19] Ilias Diakonikolas, Weihao Kong and Alistair Stewart  “Efficient algorithms and lower bounds for robust linear regression”  In _Proceedings of the Thirtieth Annual ACM-SIAM Symposium on Discrete Algorithms_ SIAM, Philadelphia, PA, 2019, pp. 2745–2754
  * [Dur19] Rick Durrett  “Probability—theory and examples” 49, Cambridge Series in Statistical and Probabilistic Mathematics  Cambridge University Press, Cambridge, 2019
  * [Fel+17] V. Feldman, E. Grigorescu, L. Reyzin, S. Vempala and Y. Xiao  “Statistical Algorithms and a Lower Bound for Detecting Planted Cliques”  In _J. ACM_ 64.2, 2017, pp. 8:1–8:37
  * [Gao20] Chao Gao  “Robust regression via mutivariate regression depth”  In _Bernoulli_ 26.2, 2020, pp. 1139–1170
  * [GC23] Xing Gao and Yu Cheng  “Robust Matrix Sensing in the Semi-Random Model”  In _Proceedings of the 37th Conference on Neural Information Processing Systems (NeurIPS)_ , 2023
  * [Ge+15] Rong Ge, Furong Huang, Chi Jin and Yang Yuan  “Escaping from saddle points—online stochastic gradient for tensor decomposition”  In _Conference on Learning Theory_ , 2015, pp. 797–842
  * [GJZ17] Rong Ge, Chi Jin and Yi Zheng  “No Spurious Local Minima in Nonconvex Low Rank Problems: A Unified Geometric Analysis”  In _Proceedings of the 34th International Conference on Machine Learning_ 70, Proceedings of Machine Learning Research  PMLR, 2017, pp. 1233–1242
  * [GLM16] Rong Ge, Jason D Lee and Tengyu Ma  “Matrix completion has no spurious local minimum”  In _Advances in Neural Information Processing Systems_ , 2016, pp. 2973–2981
  * [Hal10] Jack K Hale  “Asymptotic behavior of dissipative systems”  American Mathematical Soc., 2010
  * [Ham+86] F.. Hampel, E.. Ronchetti, P.. Rousseeuw and W.. Stahel  “Robust statistics. The approach based on influence functions”  Wiley New York, 1986
  * [HR09] P.. Huber and E.. Ronchetti  “Robust statistics”  Wiley New York, 2009
  * [Hub64] Peter J. Huber  “Robust estimation of a location parameter”  In _Annals of Mathematical Statistics_ 35, 1964, pp. 73–101
  * [IPL23] Eirini Ioannou, Muni Sreenivas Pydi and Po-Ling Loh  “Robust empirical risk minimization via Newton’s method”  In _arXiv preprint arXiv:2301.13192_ , 2023
  * [Jin+17] Chi Jin, Rong Ge, Praneeth Netrapalli, Sham M. Kakade and Michael I. Jordan  “How to Escape Saddle Points Efficiently”  In _Proceedings of the 34th International Conference on Machine Learning_ 70, Proceedings of Machine Learning Research  PMLR, 2017, pp. 1724–1732
  * [Jin+21] Chi Jin, Praneeth Netrapalli, Rong Ge, Sham M. Kakade and Michael I. Jordan  “On nonconvex optimization for machine learning: gradients, stochasticity, and saddle points”  In _Journal of the ACM_ 68.2, 2021, pp. Art. 1129
  * [Kea98] M.. Kearns  “Efficient noise-tolerant Learning from Statistical Queries”  In _Journal of the ACM_ 45.6, 1998, pp. 983–1006
  * [Li+08] J.Z. Li, D.M. Absher, H. Tang, A.M. Southwick, A.M. Casto, S. Ramachandran, H.M. Cann, G.S. Barsh, M. Feldman, L.L. Cavalli-Sforza and R.M. Myers  “Worldwide human relationships inferred from genome-wide patterns of variation”  In _Science_ 319, 2008, pp. 1100–1104
  * [Li+20] Xiao Li, Zhihui Zhu, Anthony Man-Cho So and Rene Vidal  “Nonconvex robust low-rank matrix recovery”  In _SIAM Journal on Optimization_ 30.1, 2020, pp. 660–686
  * [Li+20a] Yuanxin Li, Yuejie Chi, Huishuai Zhang and Yingbin Liang  “Non-convex low-rank matrix recovery with arbitrary outliers via median-truncated gradient descent”  In _Information and Inference: A Journal of the IMA_ 9.2, 2020, pp. 289–325
  * [LRV16] K.. Lai, A.. Rao and S. Vempala  “Agnostic Estimation of Mean and Covariance”  In _focs2016_ , 2016, pp. 665–674
  * [LW23] Shuyao Li and Stephen J. Wright  “A randomized algorithm for nonconvex minimization with inexact evaluations and complexity guarantees”  In _arXiv preprint arXiv:2310.18841_ , 2023
  * [Pas+10] P. Paschou, J. Lewis, A. Javed and P. Drineas  “Ancestry Informative Markers for Fine-Scale Individual Assignment to Worldwide Populations”  In _Journal of Medical Genetics_ 47, 2010, pp. 835–847
  * [Pra+20] Adarsh Prasad, Arun Sai Suggala, Sivaraman Balakrishnan and Pradeep Ravikumar  “Robust estimation via robust gradient estimation”  In _Journal of the Royal Statistical Society. Series B. Statistical Methodology_ 82.3, 2020, pp. 601–627
  * [Ros+02] N. Rosenberg, J. Pritchard, J. Weber, H. Cann, K. Kidd, L.A. Zhivotovsky and M.W. Feldman  “Genetic structure of human populations”  In _Science_ 298, 2002, pp. 2381–2385
  * [RRT17] Maxim Raginsky, Alexander Rakhlin and Matus Telgarsky  “Non-convex learning via Stochastic Gradient Langevin Dynamics: a nonasymptotic analysis”  In _Proceedings of the 2017 Conference on Learning Theory_ 65, Proceedings of Machine Learning Research  PMLR, 2017, pp. 1674–1703
  * [SKL17] J. Steinhardt, P. Koh and P.. Liang  “Certified Defenses for Data Poisoning Attacks”  In _Advances in Neural Information Processing Systems 30_ , 2017, pp. 3520–3532
  * [SQW16] Ju Sun, Qing Qu and John Wright  “A geometric analysis of phase retrieval”  In _Information Theory (ISIT), 2016 IEEE International Symposium on_ , 2016, pp. 2379–2383  IEEE
  * [SQW17] Ju Sun, Qing Qu and John Wright  “Complete Dictionary Recovery Over the Sphere I: Overview and the Geometric Picture”  In _IEEE Trans. Inf. Theor._ 63.2, 2017, pp. 853–884
  * [Tuk75] J.W. Tukey  “Mathematics and picturing of data”  In _Proceedings of ICM_ 6, 1975, pp. 523–531
  * [WM22] John Wright and Yi Ma  “High-dimensional data analysis with low-dimensional models: Principles, computation, and applications”  Cambridge University Press, 2022
  * [Yin+19] Dong Yin, Yudong Chen, Ramchandran Kannan and Peter Bartlett  “Defending Against Saddle Point Attack in Byzantine-Robust Distributed Learning”  In _Proceedings of the 36th International Conference on Machine Learning_ 97, Proceedings of Machine Learning Research  PMLR, 2019, pp. 7074–7084

## Supplementary Material

Supplementary material is organized as follows. In Appendix A, we provide
useful auxiliary facts and relevant technical results from previous works.
Appendix B proves our result for general robust nonconvex optimization
(Theorem 3.1). Appendix C provides omitted computation and proofs for robust
low rank matrix sensing (Section 3.1). Appendix D proves our SQ lower bound
(Section 4) for the sample complexity of efficient algorithms for the outlier-
robust low rank matrix sensing problem, and Appendix E discusses the intuition
why some simple algorithms that violate our SQ lower bound fail.

##  Appendix A Technical Preliminaries

###  A.1 Notation and Auxiliary Facts

In the appendix, we use $\widetilde{O}\\!\left(\cdot\right)$ notation to
suppress, for conciseness, logarithmic dependences on all defined quantities,
even if they do not appear inside $\widetilde{O}\\!\left(\cdot\right)$. More
precise statements can be found in the corresponding main parts of the paper.
For matrix space $\mathbb{R}^{d\times r}$ and a set
$\mathcal{X}^{\star}\subset\mathbb{R}^{d\times r}$, we use
$\mathcal{P}_{\mathcal{X}^{\star}}(\cdot)$ to denote the Frobenius projection
onto $\mathcal{X}^{\star}$, i.e.,
$\mathcal{P}_{\mathcal{X}^{\star}}(U)=\operatorname*{\arg\min}_{Z\in\mathcal{X}^{\star}}\left\|U-Z\right\|_{F}^{2}$.
We use $\operatorname{dist}(U,\mathcal{X}^{\star})$ to denote
$\left\|U-\mathcal{P}_{\mathcal{X}^{\star}}(U)\right\|_{F}$.

######  Fact A.1.

For matrices $A,B$ with compatible dimensions,

$\displaystyle\left\|AB\right\|_{F}$
$\displaystyle\leq\left\|A\right\|_{\operatorname{op}}\left\|B\right\|_{F}$
$\displaystyle\left\|AB\right\|_{F}$
$\displaystyle\leq\left\|A\right\|_{F}\left\|B\right\|_{\operatorname{op}}$

For two matrices $A$ and $B$, let $A\otimes B$ denote the Kronecker product of
$A$ and $B$.

######  Fact A.2.

For matrices $A,B,C$ with compatible dimensions,
$\operatorname{vec}({ABC})=\left({C}^{\top}\otimes{A}\right)\operatorname{vec}({B})$.

######  Fact A.3.

$\left\|A\otimes
B\right\|_{\operatorname{op}}=\left\|A\right\|_{\operatorname{op}}\left\|B\right\|_{\operatorname{op}}$

We will frequently use the following fact about the mean and the variance of
the quadratic form for zero-mean multivariate normal distributions.

######  Fact A.4.

Let $X\sim\mathcal{N}(0,\Sigma)$ be a $k$-dimensional random variable and let
$G$ be a $k\times k$ real matrix; then

$\displaystyle\operatorname*{\mathbb{E}}[X^{\top}GX]$
$\displaystyle=\operatorname{tr}(G\Sigma)=\langle G,\Sigma\rangle$
$\displaystyle\operatorname{Var}[X^{\top}GX]$
$\displaystyle=\operatorname{tr}(G\Sigma(G+G^{\top})\Sigma)$

###  A.2 Simplified Proof for the Optimization Algorithm with Inexact
Gradients and Hessians

1 for  _$t=1,2,\dots$_ do

2 Obtain $\widetilde{g}_{t}$ from the inexact gradient oracle

3

4 if  _$\|\widetilde{g}_{t}\| >\epsilon_{g}$_ then

5 $x_{t+1}=x_{t}-\frac{1}{L}\widetilde{g}_{t}$

6 else

7 Obtain $\widetilde{H}_{t}$ from the inexact Hessian oracle

8 Compute smallest eigenvalue and its corresponding eigenvector
$(\widetilde{\lambda}_{t},\widetilde{p}_{t})$

9 if  _$\widetilde{\lambda}_{t} <-\epsilon_{H}$_ then

10 $\sigma_{t}=\pm 1$ with probability $\frac{1}{2}$

11 $x_{t+1}=x_{t}+\frac{2\epsilon_{H}}{L_{H}}\sigma_{t}\widetilde{p}_{t}$

12 else

13 return $x_{t}$

14

15

Algorithm 2 Nonconvex minimization with inexact gradients and Hessians [LW23]

We stated the following guarantee (Proposition 2.2) in Section 2:

######  Proposition A.5.

Suppose a function $f$ is bounded below by $f^{*}>-\infty,$ has
$L_{g}$-Lipschitz gradient and $L_{H}$-Lipschitz Hessian, and its inexact
gradient and Hessian computation $\widetilde{g}_{t}$ and $\widetilde{H}_{t}$
satisfy
$\left\|\widetilde{g}_{t}-\nabla{f}(x_{t})\right\|\leq\frac{1}{3}\epsilon_{g}$
and
$\big{\|}\widetilde{H}_{t}-\nabla^{2}f(x_{t})\big{\|}_{\operatorname{op}}\leq\frac{2}{9}\epsilon_{H}$.
Then there exists an algorithm (Algorithm 2) with the following guarantees:

  1. 1.

(Correctness) If Algorithm 2 terminates and outputs $x_{n}$, then $x_{n}$ is a
$(\frac{4}{3}\epsilon_{g},\frac{4}{3}\epsilon_{H})$-approximate second-order
stationary point.

  2. 2.

(Runtime) Algorithm 2 terminates with probability $1$. Let
$C_{\epsilon}:=\min\left(\frac{\epsilon_{g}^{2}}{6L_{g}},\frac{2\epsilon_{H}^{3}}{9L_{H}^{2}}\right)$.
With probability at least $1-\delta$, Algorithm 2 terminates after $k$
iterations for

$k=O\Big{(}\frac{f(x_{0})-f^{*}}{C_{\epsilon}}+\frac{L_{H}^{2}L_{g}^{2}\epsilon_{g}^{2}}{\epsilon_{H}^{6}}\log\Bigl{(}\frac{1}{\delta}\Bigr{)}\Big{)}\;.$

This section proves the correctness and a slightly weaker runtime guarantee in
terms of the dependence on $\delta$. We prove a $O(1/\delta)$ dependence
instead of $O(\log(1/\delta))$, i.e., with probability at least $1-\delta$,
Algorithm 2 terminates after $k$ iterations for

$k=O\Big{(}\frac{f(x_{0})-f^{*}}{\delta C_{\epsilon}}\Big{)}\;.$ (8)

###### Proof of correctness and a weaker runtime [LW23].

We prove the correctness first. From the stopping criteria, we have
$\left\|\widetilde{g}_{k}\right\|\leq\epsilon_{g}$ and
$\widetilde{\lambda}_{k}\geq-\epsilon_{H}$. Thus,

$\left\|\nabla
f(x_{k})\right\|\leq\left\|g_{k}\right\|+\frac{1}{3}\epsilon_{g}=\frac{4}{3}\epsilon_{g}.$

To bound $\lambda_{\min}(\nabla^{2}f(x_{k}))$, write
$U_{k}=\widetilde{H}_{k}-\nabla^{2}f(x_{k})$. By Hessian inexactness
condition,
$\lambda_{\min}(-U_{k})\geq-\left\|-U_{k}\right\|_{\operatorname{op}}\geq-\frac{2}{9}\epsilon_{H}$.

We use Weyl’s theorem to conclude that

$\lambda_{\min}(\nabla^{2}f(x_{k}))\geq\lambda_{\min}(\widetilde{H}_{k})+\lambda_{\min}(-U_{k})\geq-\epsilon_{H}-\frac{2}{9}\epsilon_{H}>-\frac{4}{3}\epsilon_{H},$

as required.

Now we analyze the runtime. We proceed by first establishing an expectation
bound and then use Markov’s inequality.

Write $v_{t}:=\widetilde{g}_{t}-\nabla f(x_{t})$. When the algorithm takes a
gradient step, we have $\|\widetilde{g}_{t}\|>\epsilon_{g}$. By gradient
inexactness condition, we have
$\|v_{t}\|\leq\frac{1}{3}\|\widetilde{g}_{t}\|$, so that $\|\nabla
f(x_{t})\|\geq\frac{2}{3}\|\widetilde{g}_{t}\|$. From Taylor’s Theorem, we
have

$\displaystyle f(x_{t+1})$
$\displaystyle=f\left(x_{t}-\frac{1}{L_{g}}\widetilde{g}_{t}\right)$
$\displaystyle\leq f(x_{t})-\frac{1}{L_{g}}\nabla
f(x_{t})^{\top}\widetilde{g}_{t}+\frac{L_{g}}{2}\cdot\frac{1}{L_{g}^{2}}\|\widetilde{g}_{t}\|^{2}$
$\displaystyle=f(x_{t})-\frac{1}{2L_{g}}\left(\|\nabla
f(x_{t})\|^{2}-\|v_{t}\|^{2}\right)$ $\displaystyle\leq
f(x_{t})-\frac{1}{2L_{g}}\left(\frac{4}{9}\|\widetilde{g}_{t}\|^{2}-\frac{1}{9}\|\widetilde{g}_{t}\|^{2}\right)$
$\displaystyle\leq f(x_{t})-\frac{1}{6L_{g}}\|\widetilde{g}_{t}\|^{2}$
$\displaystyle\leq f(x_{t})-\frac{1}{6L_{g}}\epsilon_{g}^{2},$ (9)

For the negative curvature step, recall that
$\big{\|}\widetilde{H}_{t}-\nabla^{2}f(x_{t})\big{\|}_{\operatorname{op}}\leq\frac{2}{9}\epsilon_{H}$.
It follows from Taylor’s theorem that

$\displaystyle f(x_{t+1})$
$\displaystyle=f(x_{t}+\frac{2\epsilon_{H}}{L_{H}}\sigma_{t}\widetilde{p}_{t})$
$\displaystyle\leq f(x_{t})+2\frac{\epsilon_{H}}{L_{H}}\nabla
f(x_{t})^{\top}\sigma_{t}\widetilde{p}_{t}+\frac{1}{2}\frac{4\epsilon_{H}^{2}}{L_{H}^{2}}\widetilde{p}_{t}^{\top}\nabla^{2}f(x_{t})\widetilde{p}_{t}+\frac{L_{H}}{6}\frac{8\epsilon_{H}^{3}}{L_{H}^{3}}$
$\displaystyle=f(x_{t})+\frac{2\epsilon_{H}^{2}}{L_{H}^{2}}\left(\widetilde{p}_{t}^{\top}\nabla^{2}f(x_{t})\widetilde{p}_{t}+\frac{2}{3}\epsilon_{H}\right)+2\frac{\epsilon_{H}}{L_{H}}\nabla
f(x_{t})^{\top}\sigma_{t}\widetilde{p}_{t}$

When the algorithm takes a negative curvature step, we have
$\widetilde{\lambda}_{t}<-\epsilon_{H}<0$, so by Hessian inexactness
condition, we have
$\|\widetilde{H}_{t}-\nabla^{2}f(x_{t})\|\leq\frac{2}{9}\epsilon_{H}\leq\frac{2}{9}|\widetilde{\lambda}_{t}|$.
It follows from the definition of operator norm and Cauchy-Schwarz that
$|\widetilde{p}_{t}^{\top}\widetilde{H}_{t}\widetilde{p}_{t}-\widetilde{p}_{t}^{\top}\nabla^{2}f(x_{t})\widetilde{p}_{t}|\leq\|\widetilde{H}_{t}\widetilde{p}_{t}-\nabla^{2}f(x_{t})\widetilde{p}_{t}\|\leq\frac{2|\widetilde{\lambda}_{t}|}{9}$,
so

$\widetilde{p}_{t}^{\top}\nabla^{2}f(x_{t})\widetilde{p}_{t}\leq\frac{2}{9}|\widetilde{\lambda}_{t}|+{p}_{t}^{\top}\widetilde{H}_{t}\widetilde{p}_{t}=\frac{2}{9}|\widetilde{\lambda}_{t}|+(-|\widetilde{\lambda}_{t}|)=-\frac{7}{9}|\widetilde{\lambda}_{t}|.$

We thus have

$\displaystyle f(x_{t+1})$ $\displaystyle\leq
f(x_{t})+\frac{2\epsilon_{H}^{2}}{L_{H}^{2}}\left(\widetilde{p}_{t}^{\top}\nabla^{2}f(x_{t})\widetilde{p}_{t}+\frac{2}{3}\epsilon_{H}\right)+2\frac{\epsilon_{H}}{L_{H}}\nabla
f(x_{t})^{\top}\sigma_{t}\widetilde{p}_{t}$ $\displaystyle\leq
f(x_{t})+\frac{2\epsilon_{H}^{2}}{L_{H}^{2}}\left(-\frac{7}{9}|\widetilde{\lambda}_{t}|+\frac{2}{3}\alpha_{t}\right)+2\frac{\epsilon_{H}}{L_{H}}\nabla
f(x_{t})^{\top}\sigma_{t}\widetilde{p}_{t}$ $\displaystyle\leq
f(x_{t})-\frac{2\epsilon_{H}^{3}}{9L_{H}^{2}}+2\frac{\epsilon_{H}}{L_{H}}\nabla
f(x_{t})^{\top}\sigma_{t}\widetilde{p}_{t}.$ (10)

We have the following result for expected stopping time of Algorithm 2. Here
the expectation is taken with respect to the random variables $\sigma_{t}$
used at the negative curvature iterations. For purposes of this and later
results, we define

$C_{\epsilon}:=\min\left(\frac{\epsilon_{g}^{2}}{6L_{g}},\frac{2\epsilon_{H}^{3}}{9L_{H}^{2}}\right).$
(11)

Assuming Lemma A.6, the runtime guarantee in Proposition A.5 follows from
Markov inequality. ∎

######  Lemma A.6.

Consider the same setting as in Proposition A.5. Let $T$ denote the iteration
at which Algorithm 2 terminates. Then $T<\infty$ almost surely and

$\operatorname*{\mathbb{E}}T\leq\frac{f(x_{0})-f^{*}}{C_{\epsilon}},$ (12)

where $C_{\epsilon}$ is defined in (11).

The proof of this result is given below. It constructs a supermartingale111A
supermartingale with respect to filtration
$\\{\mathcal{G}_{1},\mathcal{G}_{2},\dotsc\\}$ is a sequence of random
variables $\\{Y_{1},Y_{2},\dotsc\\}$ such that for all $k\in\mathbb{Z}_{+}$,
(i) $\operatorname*{\mathbb{E}}\left|{Y_{t}}\right|<\infty$, (ii) $Y_{t}$ is
$\mathcal{G}_{t}$-measurable, and (iii)
$\operatorname*{\mathbb{E}}(Y_{t+1}\,|\mathcal{G}_{t})\leq Y_{t}$. based on
the function value and uses a supermartingale convergence theorem and optional
stopping theorem to obtain the final result. A similar proof technique is used
in [Ber+22] but for a line-search algorithm. We collect several relevant facts
about supermartingales before proving the result.

First, we need to ensure the relevant supermartingale is well defined even
after the algorithm terminates, so that it is possible to let the index $t$ of
the supermartingale go to $\infty$.

######  Fact A.7 ([Dur19, Theorem 4.2.9]).

If $T$ is a stopping time and $X_{t}$ is a supermartingale, then
$X_{\min(T,t)}$ is a supermartingale.

The following supermatingale convergence theorem will be used to ensure the
function value converges, so that the algorithm terminates with probability 1.

######  Fact A.8 (Supermartingale Convergence Theorem, [Dur19, Theorem
4.2.12]).

If $X_{t}\geq 0$ is a supermartingale, then as $t\rightarrow\infty$, there
exists a random variable $X$ such that $X_{t}\rightarrow X\;a.s.$ and
$\operatorname*{\mathbb{E}}X\leq\operatorname*{\mathbb{E}}X_{0}$.

Finally, we will use the optional stopping theorem to derive the expected
iteration complexity. Note that we use a version of the optional stopping
theorem specific to nonnegative supermartingales that does not require uniform
integrability.

######  Fact A.9 (Optional Stopping Theorem, [Dur19, Theorem 4.8.4]).

If $X_{n}$ is a nonnegative supermartingale and $N\leq\infty$ is a stopping
time, then
$\operatorname*{\mathbb{E}}X_{0}\geq\operatorname*{\mathbb{E}}X_{N}$.

###### Proof of Lemma A.6.

We first construct a supermartingale based on function values. Since
$\operatorname*{\mathbb{E}}\left[\sigma_{t}\right]=0$, linearity of
expectation implies that
$\operatorname*{\mathbb{E}}\left[2\frac{\alpha_{t}}{M}\nabla
f(x_{t})^{\top}\sigma_{t}\widetilde{p}_{t}\middle|x_{t}\right]=0$. We
therefore have from (10) that

$\displaystyle\operatorname*{\mathbb{E}}\left[f(x_{t+1})\middle|x_{t}\right]$
$\displaystyle\leq f(x_{t})-\frac{2\epsilon_{H}^{3}}{9L_{H}^{2}}.$

By combining with the (deterministic) first-order decrease estimate (9), we
have

$\operatorname*{\mathbb{E}}\left[f(x_{t+1})\middle|x_{t}\right]\leq
f(x_{t})-\min\left(\frac{\epsilon_{g}^{2}}{6L_{g}},\frac{2\epsilon_{H}^{3}}{9L_{H}^{2}}\right)=f(x_{t})-C_{\epsilon}.$

Consider the stochastic process $M_{t}:=f(x_{t})+tC_{\epsilon}$. We have

$\displaystyle\operatorname*{\mathbb{E}}\left[M_{t+1}\middle|x_{t}\right]$
$\displaystyle=\operatorname*{\mathbb{E}}\left[f(x_{t+1})+(t+1)C_{\epsilon}\middle|x_{t}\right]$
$\displaystyle\leq\operatorname*{\mathbb{E}}\left[f(x_{t})-C_{\epsilon}+(t+1)C_{\epsilon}\middle|x_{t}\right]$
$\displaystyle=\operatorname*{\mathbb{E}}\left[f(x_{t})+tC_{\epsilon}\middle|x_{t}\right]=M_{t}.$

We need to select a filtration to define the supermartingale $M_{t}$. We view
$x_{t}$ as random variables defined with respect to $\sigma_{i},i\leq k$.
Since $M_{t}$ is expressed as a function of $x_{t}$ only, we define the
filtration $\\{\mathcal{G}_{t}\\}$ to be the filtration generated by $x_{t}$,
and it naturally holds that $M_{t}$ is $\mathcal{G}_{t}$-measurable for all
$k$ and
$\operatorname*{\mathbb{E}}\left[M_{t+1}\middle|\mathcal{G}_{t}\right]=\operatorname*{\mathbb{E}}\left[M_{t+1}\middle|x_{t}\right]$.
Hence, $\\{M_{t}\\}$ is a supermartingale with respect to filtration
$\\{\mathcal{G}_{t}\\}$.

Let $T$ denote the iteration at which our algorithm stops. Since the decision
to stop at iteration $t$ depends only on $x_{t}$, we have
$\\{T=t\\}\in\mathcal{G}_{t}$, which implies $T$ is a stopping time.

We will use the supermartingale convergence theorem (Fact A.8) to show that
$T<+\infty$ almost surely, since the function value cannot decrease
indefinitely as it is bounded by $f^{*}$ from below. To apply Fact A.8, we
need to let $t\rightarrow\infty$, so we need to transform $\\{M_{t}\\}$ to
obtain a supermartingale $\\{Y_{t}\\}$that is well defined even after the
algorithm terminates.

It follows from Fact A.7 that $Y_{t}:=M_{\min(t,T)}$ is also a
supermartingale. Since $Y_{t}\geq f^{*}$, it follows from the supermartingale
convergence theorem (Fact A.8) applied to $Y_{t}-f^{*}$ that $Y_{t}\rightarrow
Y_{\infty}$ almost surely for some random variable $Y_{\infty}$ with
$\operatorname*{\mathbb{E}}Y_{\infty}\leq\operatorname*{\mathbb{E}}Y_{0}=\operatorname*{\mathbb{E}}M_{0}=f(x_{0})<\infty$.
Hence $\mathbb{P}[Y_{\infty}=+\infty]=0$. On the other hand, as
$t\rightarrow\infty$, we have $tC_{\epsilon}\rightarrow\infty$, so
$T=+\infty\implies Y_{t}=M_{t}\geq
f^{*}+tC_{\epsilon}\rightarrow\infty\implies Y_{\infty}=+\infty$. Therefore we
have $\mathbb{P}[T<+\infty]=1$.

We can then apply the optional stopping theorem (Fact A.9) to $Y_{t}-f^{*}$.
It follows that

$f^{*}+\operatorname*{\mathbb{E}}T\cdot
C_{\epsilon}\leq\operatorname*{\mathbb{E}}f(x_{T})+\operatorname*{\mathbb{E}}T\cdot
C_{\epsilon}=\operatorname*{\mathbb{E}}\left[M_{T}\right]=\operatorname*{\mathbb{E}}\left[Y_{T}\right]\leq\operatorname*{\mathbb{E}}[Y_{0}]=\operatorname*{\mathbb{E}}\left[M_{0}\right]=f(x_{0}),$

where the first equality uses $T<+\infty$ almost surely and the last
inequality is Fact A.9. By reorganizing this bound, we obtain
$\operatorname*{\mathbb{E}}T\leq\frac{f(x_{0})-f^{*}}{C_{\epsilon}}$, as
desired. ∎

###  A.3 Robust Mean Estimation—Omitted background

A class of algorithms that robustly estimate quantities in the presence of
outliers under strong contamination model (Definition 1.1) are based on the
stability of samples:

######  Definition A.10 (Definition of Stability [Dia+17]).

Fix $0<\epsilon<1/2$ and $\delta>\epsilon$. A finite set
$S\subset\mathbb{R}^{k}$ is $(\epsilon,\delta)$-stable with respect to mean
$\mu\in\mathbb{R}^{k}$ and $\sigma^{2}$ if for every $S^{\prime}\subset S$
with $|S^{\prime}|\geq(1-\epsilon)|S|$, the following conditions hold: (i)
$\left\|\mu_{S^{\prime}}-\mu\right\|\leq\sigma\delta$, and (ii)
$\left\|\bar{\Sigma}_{S^{\prime}}-\sigma^{2}I\right\|_{\operatorname{op}}\leq\sigma^{2}\delta^{2}\epsilon$,
where $\mu_{S^{\prime}}$ and $\bar{\Sigma}_{S^{\prime}}$ denote the empirical
mean and empirical covariance over the set $S^{\prime}$ respectively.

Data: $0<\epsilon<1/2$ and $T$ is an $\epsilon$-corrupted set

Result: $\widehat{\mu}$ with
$\left\|\widehat{\mu}-\mu_{S}\right\|=O\big{(}\sqrt{\left\|\Sigma\right\|_{\operatorname{op}}\epsilon}\big{)}$

1 Initialize a weight function $w:T\rightarrow\mathbb{R}_{\geq 0}$ with
$w(x)=1/|T|$ for all $x\in T$;

2 while  _$\left\|w\right\|_{1}\geq 1-2\epsilon$_ do

3 $\mu(w):=\frac{1}{\left\|w\right\|_{1}}\sum_{x\in T}w(x)x$;

4 $\Sigma(w):=\frac{1}{\left\|w\right\|_{1}}\sum_{x\in
T}w(x)(x-\mu(w))(x-\mu(w))^{\top}$;

5 Compute the largest eigenvector $v$ of $\Sigma(w)$ ;

6 $g(x):=|v^{\top}(x-\mu(w))|^{2}$;

7 Find the largest threshold $t$ such that $\Sigma_{x\in T:g(x)\geq
t}w(x)\geq\epsilon$;

8 $f(x):=g(x)\mathbbm{1}{\\{g(x)\geq t\\}}$ ;

9 $w(x)\leftarrow w(x)\left(1-\frac{f(x)}{\max_{y\in T:w(y)\neq
0}f(y)}\right)$ ;

10

return $\mu(w)$

Algorithm 3 $\mathsf{RobustMeanEstimation}$ with unknown covariance bound

On input the corrupted version of a stable set, Algorithm 3 has the following
guarantee.

######  Proposition A.11 (Robust Mean Estimation with Stability [DKP20,
Theorem A.3]).

Let $T\subset\mathbb{R}^{k}$ be an $\epsilon$-corrupted version of a set $S$,
where $S$ is $(C\epsilon,\delta)$-stable with respect to $\mu_{S}$ and
$\sigma^{2}$, and where $C>0$ is a sufficiently large constant. Algorithm 3 on
input $\epsilon$ and $T$ (but not $\sigma$ or $\delta$) returns
(deterministically) a vector $\widehat{\mu}$ in polynomial time so that
$\left\|\mu_{S}-\widehat{\mu}\right\|=O(\sigma\delta)$.

Proposition A.11 requires the uncorrupted samples to be stable. With a large
enough sample size, most independent and identically distributed samples from
a bounded covariance distribution are stable. The remaining samples can be
treated as corruptions.

######  Proposition A.12 (Sample Complexity for Stability [DKP20, Theorem
1.4]).

Fix any $0<\xi<1$. Let $S$ be a multiset of $n$ independent and identically
distributed samples from a distribution on $\mathbb{R}^{k}$ with mean $\mu$
and covariance $\Sigma$. Then, with probability at least $1-\xi$, there exists
a sample size $n=O\\!\left(\frac{k\log k+\log(1/\xi)}{\epsilon}\right)$ and a
subset $S^{\prime}\subseteq S$ such that
$\left|S^{\prime}\right|\geq\left(1-\epsilon\right)n$ and $S^{\prime}$ is
$\left(2\epsilon,\delta\right)$-stable with respect to $\mu$ and
$\left\|\Sigma\right\|_{\operatorname{op}}$, where
$\delta=O\\!\left(\sqrt{\epsilon}\right)$.

###### Proof of Proposition 2.3.

Proposition A.12 implies that for i.i.d. samples from a $k$-dimensional
bounded covariance distribution to contain a stable set of more than
$(1-\epsilon)$-fraction of samples with high probability, we need
$\widetilde{O}\\!\left(\frac{k}{\epsilon}\right)$ samples. We refer to the
remaining $\epsilon$-fraction of samples as unstable samples. Since the
adversary corrupts an $\epsilon$-fraction of clean samples $S$, the input set
$T$ can be considered as a $2\epsilon$-corrupted version of a stable set, if
we view the unstable samples as corruptions. Therefore, Proposition A.11
applies and gives the desired error guarantee. ∎

##  Appendix B Omitted Proofs in General Robust Nonconvex Optimization

##### Sample Size for Successful Robust Mean Estimation

For each iteration $t$ in Algorithm 2 and Algorithm 1, we would like to
robustly estimate the mean of a set of corrupted points
$\\{\operatorname{vec}{(\nabla f(x_{t}))}\\}_{i=1}^{n}$ and/or
$\\{\operatorname{vec}{(\nabla^{2}f(x_{t}))}\\}_{i=1}^{n}$ where the inliers
are drawn from a distribution with bounded covariance
$\Sigma\preceq\sigma^{2}I$. For fixed $x_{t}$, inliers are drawn independently
and identically distributed (i.i.d.), but the dependence across all iterations
$t$ is allowed.

######  Theorem B.1.

Let $\mathcal{X}\subset\mathbb{R}^{l}$ be a closed and bounded set with radius
at most $\gamma$ (i.e., $\left\|x\right\|\leq\gamma,\;\forall
x\in\mathcal{X}$). Let $p^{*}$ be a distribution over functions
$h:\mathcal{X}\rightarrow\mathbb{R}^{k}$. Let $\xi\in(0,1)$ be the failure
probability. Suppose $h$ is $L$-Lipschitz and uniformly bounded, i.e., there
exists $B>0$ such that $\left\|h(x)\right\|\leq B\;\forall x\in\mathcal{X}$
almost surely. Assume further that for each $x\in\mathcal{X}$ we have
$\left\|\operatorname{Cov}_{h\sim
p^{*}}(h(x))\right\|_{\operatorname{op}}\leq\sigma^{2}$. Let $S$ be a multiset
of $n$ i.i.d. samples $\\{h_{i}\\}_{i=1}^{n}$ from $p^{*}$. Then there exists
a sample size

$n=O\left(\frac{k\log k+l\log\left(\frac{\gamma
LB}{\sigma\epsilon\xi}\right)}{\epsilon}\right)$

such that with probability $1-\xi$, it holds that for each $x\in\mathcal{X}$,
there exists a subset $S^{\prime}\subset S$ with
$|S^{\prime}|\geq(1-2\epsilon)n$ (potentially different subsets $S^{\prime}$
for different $x$) such that $\\{h_{i}(x)\\}_{i\in S^{\prime}}$ is
$(2\epsilon,\delta)$-stable with respect to $\mu$ and $\sigma^{2}$, where
$\delta=O\\!\left(\sqrt{\epsilon}\right)$.

To prove this theorem, we work with an easier version of stability condition
specialized to samples from a bounded covariance distribution.

######  Claim B.2 (Claim 2.1 in [DKP20]).

(Stability for bounded covariance) Let $R\subset\mathbb{R}^{k}$ be a finite
multiset such that $\left\|\mu_{R}-\mu\right\|\leq\sigma\delta$, and
$\left\|\bar{\Sigma}_{R}-\sigma^{2}I\right\|_{\operatorname{op}}\leq\sigma^{2}\delta^{2}/\epsilon$
for some $0\leq\epsilon\leq\delta$. Then $R$ is
$\left(\Theta(\epsilon),\delta^{\prime}\right)$ stable with respect to $\mu$
and $\sigma^{2}$, where $\delta^{\prime}=O(\delta+\sqrt{\epsilon})$.

###### Proof of Theorem B.1.

Given Claim B.2, it suffices to show that with probability $1-\xi$, for each
$x\in\mathcal{X}$ there exists a subset $S^{\prime}\subset S$ with
$|S^{\prime}|\geq(1-2\epsilon)n$ such that

$\displaystyle\left\|\operatorname*{\mathbb{E}}_{i\in
S^{\prime}}h_{i}(x)-\operatorname*{\mathbb{E}}_{h\sim p^{*}}h(x)\right\|$
$\displaystyle\leq O\\!\left(\sigma\delta\right)$ (13)
$\displaystyle\left\|\operatorname*{\mathbb{E}}_{i\in
S^{\prime}}(h_{i}(x)-\operatorname*{\mathbb{E}}_{h\sim
p^{*}}h(x))(h_{i}(x)-\operatorname*{\mathbb{E}}_{h\sim
p^{*}}h(x))^{\top}-\sigma^{2}I\right\|_{\operatorname{op}}$ $\displaystyle\leq
O\\!\left(\sigma^{2}\delta^{2}/\epsilon\right)$ (14)

By Proposition A.12, for each $x\in\mathcal{X}$, with probability
$1-\left(\frac{\xi\gamma BL}{\sigma\sqrt{\epsilon}}\right)^{l}\geq
1-\xi/\left(\frac{\gamma BL}{\sigma\sqrt{\epsilon}}\right)^{l}$, there exists
a subset $S^{\prime}\subset S$ with $|S^{\prime}|\geq(1-2\epsilon)n$ such that
(13) and (14) hold.

We proceed with a net argument. Up to a multiplicative error that can be
suppressed by $O\\!\left(\cdot\right)$, if Equation (13) holds for some
$x\in\mathcal{X}$, it also holds for all other $x^{\prime}$ in a ball of
radius $\sigma\delta/L$ because $h(\cdot)$ is $L$-Lipschitz. Similarly, (14)
is equivalent to

$|\operatorname*{\mathbb{E}}_{i\in
S^{\prime}}((h_{i}(x)-\operatorname*{\mathbb{E}}_{h\sim
p^{*}}h(x))^{\top}v)^{2}-\sigma^{2}v^{\top}v|\leq
O\\!\left(\sigma^{2}\delta^{2}/\epsilon\right)\qquad\forall
v\in\mathbb{R}^{k}.$

Since $h(\cdot)$ is $L$-Lipschitz and uniformly bounded by $B$, if (14) holds
for some $x\in\mathcal{X}$, it also holds for all $x^{\prime}$ in a ball of
radius $\frac{\sigma^{2}\delta^{2}}{\epsilon BL}$. Therefore, it suffices for
Equation (13) and (14) to hold for a $\tau$-net of $\mathcal{X}$, where for
$\delta=O\\!\left(\sqrt{\epsilon}\right)$ we have
$\tau=\min\left(\frac{\sigma\delta}{L},\frac{\sigma^{2}\delta^{2}}{\epsilon
BL}\right)=\Omega\left(\frac{\sigma\sqrt{\epsilon}}{BL}\right)$. An
$\Omega\left(\frac{\sigma\sqrt{\epsilon}}{BL}\right)$-net of $\gamma$-radius
ball in $\mathbb{R}^{l}$ has size $O\Big{(}\left(\frac{\gamma
BL}{\sigma\sqrt{\epsilon}}\right)^{l}\Big{)}$. Taking a union bound over this
net completes the proof. ∎

##### Proof of Main Theorem

We establish a slightly enhanced version of Theorem 3.1 by including the last
sentence as an additional component to the original theorem.

######  Theorem B.3.

Suppose we are given $\epsilon$-corrupted set of functions
$\\{f_{i}\\}_{i=1}^{n}$ for sample size $n$, generated according to Definition
1.2. Suppose Assumption 1.4 holds in a bounded region
$\mathcal{B}\subset\mathbb{R}^{D}$ of radius $\Gamma$ with gradient and
Hessian covariance bound $\sigma_{g}$ and $\sigma_{H}$ respectively, and we
have an arbitrary initialization $x_{0}\in\mathcal{B}$. For some
$n=\widetilde{O}\big{(}D^{2}/\epsilon\big{)}$, Algorithm 2 initialized at
$x_{0}$ outputs an $(\epsilon_{g},\epsilon_{H})$-approximate SOSP in
polynomial time with high probability if the following conditions hold:

  1. (I)

All iterates $x_{t}$ in Algorithm 2 stay inside the bounded region
$\mathcal{B}$.

  2. (II)

For an absolute constant $c>0$, it holds that $\sigma_{g}\sqrt{\epsilon}\leq
c\epsilon_{g}$ and $\sigma_{H}\sqrt{\epsilon}\leq c\epsilon_{H}$.

Moreover, there exists an absolute constant $C_{est}$ such that for each
iteration $t$, the gradient oracle $\widetilde{g}_{t}$ and Hessian oracle
$\widetilde{H}_{t}$ satisfy
$\left\|\widetilde{g}_{t}-\nabla{f}(x_{t})\right\|\leq
C_{est}\sigma_{g}\sqrt{\epsilon}$ and
$\big{\|}\widetilde{H}_{t}-\nabla^{2}f(x_{t})\big{\|}_{F}\leq
C_{est}\sigma_{H}\sqrt{\epsilon}$.

And we recall that the we construct the gradient and Hessian oracles in
Algorithm 2 in the following way:

$\displaystyle\widetilde{g}_{t}$
$\displaystyle\leftarrow\mathbf{RobustMeanEstimation}(\\{\nabla
f_{i}(x_{t})\\}_{i=1}^{n},4\epsilon)$ $\displaystyle\widetilde{H}_{t}$
$\displaystyle\leftarrow\mathbf{RobustMeanEstimation}(\\{\nabla^{2}f_{i}(x_{t})\\}_{i=1}^{n},4\epsilon)$

If we ignore that the dependence between $\\{\nabla f_{i}(x_{t})\\}$
introduced via $x_{t}$, this theorem follows directly from Proposition 2.3.
Here we fix this technicality via a union bound over all $x_{t}$ with Theorem
B.1.

###### Proof.

By Proposition 2.2, it suffices to check that gradient and Hessian inexactness
condition is satisfied. That is, for all iterations $t$, it holds that

$\displaystyle\left\|\widetilde{g}_{t}-\nabla{f}(x_{t})\right\|$
$\displaystyle\leq\frac{1}{3}\epsilon_{g},$ (15)
$\displaystyle\big{\|}\operatorname{vec}(\widetilde{H}_{t})-\operatorname{vec}(\nabla^{2}f(x_{t}))\big{\|}$
$\displaystyle\leq\frac{2}{9}\epsilon_{H}.$ (16)

Here we use the fact that
$\big{\|}\operatorname{vec}(\widetilde{H}_{t})-\operatorname{vec}(\nabla^{2}f(x_{t}))\big{\|}=\big{\|}\widetilde{H}_{t}-\nabla^{2}f(x_{t})\big{\|}_{F}\geq\big{\|}\widetilde{H}_{t}-\nabla^{2}f(x_{t})\big{\|}_{\operatorname{op}}$.

We proceed to apply Theorem B.1 and consider the function
$h_{g}(x,A)=\operatorname{vec}(\nabla_{x}f(x,A))$ and
$h_{H}(x)=\operatorname{vec}(\nabla_{xx}^{2}f(x,A))$. We consider
$h_{H}:\mathbb{R}^{D}\rightarrow\mathbb{R}^{D^{2}}$ first. By assumption (I),
all iterates never leave the bounded region $\mathcal{B}$ with high
probability; we condition on this event in the remaining analysis and it
follows that
$\left\|\operatorname{Cov}_{A\sim\mathcal{G}}(\operatorname{vec}(\nabla^{2}f(x_{t},A)))\right\|_{\operatorname{op}}\leq\sigma_{H}^{2}$
for all iterations $t$.

Assumption 1.4 (ii) posits that with high probability $h_{H}(\cdot,A)$ is
$L_{D_{H}}$-Lipschitz and its $\ell_{2}$-norm is bounded above by $B_{D_{H}}$.
We further condition on this event. Let $\xi\in(0,1)$. Since
$\left\|\operatorname{Cov}_{A\in\mathcal{G}}(h_{H}(x,A))\right\|_{\operatorname{op}}\leq\sigma_{H}$
by Assumption 1.4 (iii), there exists a sample size

$n=O\left(\frac{D^{2}\log D+D\log\left(\frac{\gamma
L_{D_{H}}B_{D_{H}}}{\sigma_{H}\epsilon\xi}\right)}{\epsilon}\right),$

such that with probability $1-\xi$, it holds that for each $x\in\mathcal{X}$,
there exists a $(2\epsilon,O\\!\left(\sqrt{\epsilon}\right))$-stable subset of
size at least $(1-2\epsilon)n$ in the sense of Definition A.10 (potentially
different subsets $S^{\prime}$ for different $x$). Therefore, conditioning on
the event that for all $x$ there exists a stable subset, Proposition A.11
implies that there exists an absolute constant $C_{est}$ such that
$\big{\|}\operatorname{vec}(\widetilde{H}_{t})-\operatorname{vec}(\nabla^{2}f(x_{t}))\big{\|}\leq
C_{est}{\sigma_{H}\sqrt{\epsilon}}$ for all $t$. By assumption (II) with a
sufficiently small constant $c$, we have
$C_{est}{\sigma_{H}\sqrt{\epsilon}}\leq 2\epsilon_{H}/9$, and therefore
$\big{\|}\operatorname{vec}(\widetilde{H}_{t})-\operatorname{vec}(\nabla^{2}f(x_{t}))\big{\|}\leq\frac{2}{9}\epsilon_{H}.$

Now we established that for Equation (16) to hold for all iterations $t$ with
high probability, the sample complexity is required to be at least
$n=\widetilde{O}\\!\left(D^{2}/\epsilon\right)$. A similar argument implies
that $\widetilde{O}\\!\left(D/\epsilon\right)$ samples are needed for Equation
(15) to hold for all iterations $t$ with high probability, which is dominated
by $\widetilde{O}\\!\left(D^{2}/\epsilon\right)$. ∎

##  Appendix C Omitted Proofs in Low Rank Matrix Sensing

This section provides omitted computation and proofs for robust low rank
matrix sensing (Section 3.1), establishing Theorem 1.7. Appendix C.1 computes
the gradient and Hessian of $f_{i}$ defined in Equation (3) and their
covariance upper bounds under Definition 1.6 with noiseless measurements
($\sigma=0$). Appendix C.2 proves the global convergence of our general robust
nonconvex optimization framework (Theorem 3.1) applied to the noiseless robust
low rank matrix sensing problem, which leads to an approximate SOSP, and
Appendix C.3 proves the local linear convergence to a global minimum from this
approximate SOSP. Finally, Appendix C.4 discusses the case of noisy
measurements ($\sigma\neq 0$) under Definition 1.6.

###  C.1 Omitted Computation in Low Rank Matrix Sensing

####  C.1.1 Sample Gradient and Sample Hessian

##### Summary of Results

Consider the noiseless setting as in Theorem 1.7 with $\sigma=0$. We now
compute the gradient and the Hessian of $f_{i}:\mathbb{R}^{d\times
r}\rightarrow\mathbb{R}$. We firstly summarize the results. The gradient of
$f_{i}$ at point $U$, denoted by $\nabla f_{i}(U)$, is a matrix in
$\mathbb{R}^{d\times r}$ given by

$\nabla f_{i}(U)=\langle UU^{\top}-M^{*},A_{i}\rangle(A_{i}+A_{i}^{\top})U.$

The Hessian of $f_{i}$ at point $U$, denoted by $(H_{i})_{U}$, is a linear
operator $(H_{i})_{U}:\mathbb{R}^{d\times r}\rightarrow\mathbb{R}^{d\times
r}$. The operator $(H_{i})_{U}$ acting on a matrix $Y\in\mathbb{R}^{d\times
r}$ gives:

$(H_{i})_{U}(Y)=\langle
UU^{\top}-M^{*},A_{i}\rangle(A_{i}+A_{i}^{\top})Y+\langle
Y,(A_{i}+A_{i}^{\top})U\rangle(A_{i}+A_{i}^{\top})U.$

We want to identify this linear map with a matrix $H_{i}(U)$ of dimension
$\mathbb{R}^{dr\times dr}$. We vectorize both the domain and the codomain of
$(H_{i})_{U}$ so that it can be represented by a matrix.

Let $I_{r}$ and $I_{d}$ denote the identity matrices of dimension $r\times r$
and $d\times d$ respectively.

$H_{i}(U)=\langle UU^{\top}-M^{*},A_{i}\rangle
I_{r}\otimes(A_{i}+A_{i}^{\top})+\operatorname{vec}\left((A_{i}+A_{i}^{\top})U\right)\operatorname{vec}\left((A_{i}+A_{i}^{\top})U\right)^{\top}.$

In this paper, we sometimes abuse the notation and use $\nabla^{2}f_{i}(U)$ to
refer to either $H_{i}(U)$ or $(H_{i})_{U}$, but its precise meaning should be
clear from its domain.

##### Computation of Gradients and Hessians

Let $(D^{k}f_{i})_{U}$ denote the $k$-th order derivative of $f_{i}$ at point
$U\in\mathbb{R}^{d\times r}$ and $\mathcal{L}(X,Y)$ denote the space of all
linear mappings from $X$ to $Y$. We consider higher-order derivatives as
linear maps:

$\displaystyle(Df_{i})_{U}$ $\displaystyle:\mathbb{R}^{d\times
r}\rightarrow\mathbb{R}$ $\displaystyle(D^{2}f_{i})_{U}$
$\displaystyle:\mathbb{R}^{d\times
r}\rightarrow\mathcal{L}(\mathbb{R}^{d\times r},\mathbb{R})$

We identify $(Df_{i})_{U}$ with the matrix $\nabla f_{i}(U)$ such that

$(Df_{i})_{U}(Z)=\langle\nabla f_{i}(U),Z\rangle$

and identify $(D^{2}f_{i})_{U}$ with a linear operator
$(H_{i})_{U}:\mathbb{R}^{d\times r}\rightarrow\mathbb{R}^{d\times r}$ such
that

$(D^{2}f_{i})_{U}(Y)(Z)=\langle(H_{i})_{U}(Y),Z\rangle.$

Since $f_{i}$ is differentiable, applying its derivative at $U$ to a matrix
$Z$ gives the corresponding directional derivative at $U$ for the direction
$Z$.

$\displaystyle(Df_{i})_{U}(Z)$ $\displaystyle=\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!missing}{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!t}f_{i}(U+tZ)\Bigr{|}_{t=0}$
$\displaystyle=\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!missing}{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!t}\frac{1}{2}(\langle(U+tZ)(U+tZ)^{\top},A_{i}\rangle-
y_{i})^{2}\Bigr{|}_{t=0}$ $\displaystyle=\left(\langle UU^{\top},A_{i}\rangle-
y_{i}\right)\left(\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!missing}{\mathop{\mathrm{\mathstrut d}}\nolimits\\!t}\langle
UU^{\top}+t^{2}ZZ^{\top}+t(ZU^{\top}+UZ^{\top}),A_{i}\rangle\right)\Bigr{|}_{t=0}$
$\displaystyle=\left(\langle UU^{\top},A_{i}\rangle-
y_{i}\right)\langle\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!missing}{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!t}\left(t^{2}ZZ^{\top}+t(ZU^{\top}+UZ^{\top})\right),A_{i}\rangle\Bigr{|}_{t=0}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\left(\langle
UU^{\top},A_{i}\rangle-y_{i}\right)\langle ZU^{\top}+UZ^{\top},A_{i}\rangle$
$\displaystyle=\left(\langle UU^{\top},A_{i}\rangle-y_{i}\right)\langle
A_{i},ZU^{\top}\rangle+\langle A_{i},UZ^{\top}\rangle$
$\displaystyle=\left(\langle UU^{\top},A_{i}\rangle-y_{i}\right)\langle
A_{i}U,Z\rangle+\langle U^{\top}A_{i},Z^{\top}\rangle$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}\left(\langle
UU^{\top},A_{i}\rangle-y_{i}\right)\langle(A_{i}+A_{i}^{\top})U,Z\rangle.$

From (a) we conclude $(Df_{i})_{U}$ is the functional $Z\mapsto\left(\langle
UU^{\top},A_{i}\rangle-y_{i}\right)\langle ZU^{\top}+UZ^{\top},A_{i}\rangle$;
it would be easier to calculate the second derivative from this form. From
(b), using $y_{i}=\langle A_{i},M^{*}\rangle,$ we obtain the following closed-
form expression for the gradient $\nabla f_{i}(U)$:

$\nabla f_{i}(U)=\langle UU^{\top}-M^{*},A_{i}\rangle(A_{i}+A_{i}^{\top})U.$

To calculate the second derivative of $f_{i}$ at point $U$, we study the
variation of its first derivative $U\mapsto(Df_{i})_{U}$ at direction $Y$.
Recall that the second derivative lives in the space of linear functionals
$\mathcal{L}(\mathbb{R}^{d\times r},\mathbb{R})$ and we use
$Z\in\mathbb{R}^{d\times r}$ to denote the input of this linear functional.

$\displaystyle(D^{2}f_{i})_{U}(Y)$
$\displaystyle=\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!missing}{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!t}(Df_{i})_{U+tY}\Bigr{|}_{t=0}$
$\displaystyle=\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!missing}{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!t}\left\\{Z\mapsto\left(\langle(U+tY)(U+tY)^{\top},A_{i}\rangle-
y_{i}\right)\langle
Z(U+tY)^{\top}+(U+tY)Z^{\top},A_{i}\rangle\right\\}\Bigr{|}_{t=0}$
$\displaystyle=\left(\langle UU^{\top},A_{i}\rangle-
y_{i}\right)\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!missing}{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!t}\left\\{Z\mapsto\langle
Z(U+tY)^{\top}+(U+tY)Z^{\top},A_{i}\rangle\right\\}\Bigr{|}_{t=0}+$
$\displaystyle\qquad\qquad\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!missing}{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!t}\left\\{\langle(U+tY)(U+tY)^{\top},A_{i}\rangle-
y_{i}\right\\}\Bigr{|}_{t=0}\left\\{Z\mapsto\langle
ZU^{\top}+UZ^{\top},A_{i}\rangle\right\\}$ $\displaystyle=\left(\langle
UU^{\top},A_{i}\rangle-y_{i}\right)\\{Z\mapsto\langle
ZY^{\top}+YZ^{\top},A_{i}\rangle\\}+$ $\displaystyle\qquad\qquad\langle
UY^{\top}+YU^{\top},A_{i}\rangle\left\\{Z\mapsto\langle
ZU^{\top}+UZ^{\top},A_{i}\rangle\right\\}$
$\displaystyle=Z\mapsto\left(\langle UU^{\top},A_{i}\rangle-
y_{i}\right)\langle ZY^{\top}+YZ^{\top},A_{i}\rangle+\langle
UY^{\top}+YU^{\top},A_{i}\rangle\langle ZU^{\top}+UZ^{\top},A_{i}\rangle$
$\displaystyle=Z\mapsto(\langle UU^{\top},A_{i}\rangle-
y_{i})\langle(A_{i}+A_{i}^{\top})Y,Z\rangle+\langle
Y,(A_{i}+A_{i}^{\top})U\rangle\langle(A_{i}+A_{i}^{\top})U,Z\rangle.$

Recall that we identify $(D^{2}f_{i})_{U}$ with a linear operator
$(H_{i})_{U}:\mathbb{R}^{d\times r}\rightarrow\mathbb{R}^{d\times r}$ such
that

$(D^{2}f_{i})_{U}(Y)(Z)=\langle(H_{i})_{U}(Y),Z\rangle.$

Note that $(H_{i})_{U}:\mathbb{R}^{d\times r}\rightarrow\mathbb{R}^{d\times
r}$ viewed as a fourth-order tensor is difficult to write down in a closed
form, so we vectorize both the domain and the codomain of $(H_{i})_{U}$ so
that it can be represented by a matrix.

The operator $(H_{i})_{U}$ acting on matrix $Y\in\mathbb{R}^{d\times r}$ gives

$(H_{i})_{U}(Y)=\langle
UU^{\top}-M^{*},A_{i}\rangle(A_{i}+A_{i}^{\top})Y+\langle
Y,(A_{i}+A_{i}^{\top})U\rangle(A_{i}+A_{i}^{\top})U.$

We identify this linear map with a matrix $H_{i}(U)$ of dimension
$\mathbb{R}^{dr\times dr}$, acting on the vectorized version of $Y$. As a
shorthand, write $y=\operatorname{vec}(Y)$ and $B_{i}=A_{i}+A_{i}^{\top}$.

Let $I_{r}$ and $I_{d}$ denote the identity matrix of dimension $r\times r$
and $d\times d$ respectively. Then

$\displaystyle H_{i}(U)\,y$ $\displaystyle=\operatorname{vec}\left(\langle
UU^{\top}-M^{*},A_{i}\rangle(A_{i}+A_{i}^{\top})Y+\langle
Y,(A_{i}+A_{i}^{\top})U\rangle(A_{i}+A_{i}^{\top})U\right)$
$\displaystyle=\langle
UU^{\top}-M^{*},A_{i}\rangle\operatorname{vec}\left(B_{i}YI_{r}\right)+\operatorname{vec}\left(B_{i}U\right)\operatorname{vec}(B_{i}U)^{\top}\operatorname{vec}(Y)$
$\displaystyle=\langle UU^{\top}-M^{*},A_{i}\rangle I_{r}\otimes
B_{i}\operatorname{vec}\left(Y\right)+\operatorname{vec}\left(B_{i}U\right)\operatorname{vec}(B_{i}U)^{\top}\operatorname{vec}(Y)$
$\displaystyle=\langle UU^{\top}-M^{*},A_{i}\rangle(I_{r}\otimes
B_{i})y+\operatorname{vec}\left(B_{i}U\right)\operatorname{vec}(B_{i}U)^{\top}y.$

Hence

$H_{i}(U)=\langle UU^{\top}-M^{*},A_{i}\rangle I_{r}\otimes
B_{i}+\operatorname{vec}\left(B_{i}U\right)\operatorname{vec}(B_{i}U)^{\top}.$

####  C.1.2 Gradient Covariance Bound

######  Lemma C.1 (Equation (6)).

Consider the noiseless setting as in Theorem 1.7 with $\sigma=0$. For all
$U\in\mathbb{R}^{d\times r}$ with
$\left\|U\right\|_{\operatorname{op}}^{2}\leq\Gamma$, it holds that

$\left\|\operatorname{Cov}(\operatorname{vec}(\nabla
f_{i}(U)))\right\|_{\operatorname{op}}\leq
8\left\|UU^{\top}-M^{*}\right\|_{F}^{2}\left\|U\right\|_{\operatorname{op}}^{2}\leq
32r^{2}\Gamma^{3}.$ (17)

Recall that $\Gamma\geq 36{\sigma^{\star}_{1}}$, so we have the following
trivial bound that is used frequently

$\left\|M^{*}\right\|_{\operatorname{op}}={\sigma^{\star}_{1}}\leq\Gamma/36<\Gamma.$

###### Proof of Equation (6).

Recall that

$\nabla f_{i}(U)=\langle UU^{\top}-M^{*},A_{i}\rangle(A_{i}+A_{i}^{\top})U.$

We frequently encounter the following calculations.

######  Lemma C.2.

Let $P,Q\in\mathbb{R}^{d\times d}$ be given and let $A_{i}$ have i.i.d.
standard Gaussian entries. Then

$\displaystyle\operatorname*{\mathbb{E}}[\langle P,A_{i}\rangle\langle
Q,A_{i}\rangle]$ $\displaystyle=\langle P,Q\rangle,$ (18)
$\displaystyle\operatorname{Var}[\langle P,A_{i}\rangle\langle
Q,A_{i}\rangle]$ $\displaystyle\leq
2\left\|P\right\|_{F}^{2}\left\|Q\right\|_{F}^{2}.$ (19)

###### Proof.

By Definition 1.6, $X:=\operatorname{vec}(A_{i})$ is a standard Gaussian
vector with identity covariance $\Sigma=I_{d^{2}}$. Let
$a=\operatorname{vec}(P),b=\operatorname{vec}(Q)$, and $G=ab^{\top}$. Fact A.4
implies that

$\displaystyle\operatorname*{\mathbb{E}}\langle P,A_{i}\rangle\langle
Q,A_{i}\rangle$ $\displaystyle=\langle ab^{\top},I_{dr}\rangle=\langle
a,b\rangle=\langle P,Q\rangle$ $\displaystyle\operatorname{Var}[\langle
P,A_{i}\rangle\langle Q,A_{i}\rangle]$
$\displaystyle=\operatorname{tr}(G\Sigma(G+G^{\top})\Sigma)=\operatorname{tr}(G(G+G^{\top}))$
$\displaystyle=\operatorname{tr}(ab^{\top}(ab^{\top}+ba^{\top}))$
$\displaystyle=b^{\top}a\operatorname{tr}(ab^{\top})+b^{\top}b\operatorname{tr}(aa^{\top})$
$\displaystyle=\langle
a,b\rangle^{2}+\left\|a\right\|^{2}\left\|b\right\|^{2}$ $\displaystyle\leq
2\left\|a\right\|^{2}\left\|b\right\|^{2}\text{\qquad(Cauchy-Schwarz)}$
$\displaystyle=2\left\|P\right\|_{F}^{2}\left\|Q\right\|_{F}^{2}.\qed$

To compute the mean and the variance of the gradient $\nabla f_{i}(U)$, it is
more convenient to work with $\langle\nabla f_{i}(U),Z\rangle$ for some
$Z\in\mathbb{R}^{d\times r}$ with $\left\|Z\right\|_{F}=1$.

$\displaystyle\langle\nabla f_{i}(U),Z\rangle$ $\displaystyle=\langle
UU^{\top}-M^{*},A_{i}\rangle\langle(A_{i}+A_{i}^{\top})U,Z\rangle$
$\displaystyle=\langle UU^{\top}-M^{*},A_{i}\rangle\langle
A_{i}+A_{i}^{\top},ZU^{\top}\rangle$ $\displaystyle=\langle
UU^{\top}-M^{*},A_{i}\rangle\langle ZU^{\top}+UZ^{\top},A_{i}\rangle.$

Lemma C.2 implies that the expectation

$\displaystyle\operatorname*{\mathbb{E}}\langle\nabla f_{i}(U),Z\rangle$
$\displaystyle=\langle UU^{\top}-M^{*},ZU^{\top}+UZ^{\top}\rangle$
$\displaystyle=\langle\left((UU^{\top}-M^{*})+(UU^{\top}-M^{*})^{\top}\right)U,Z\rangle$
$\displaystyle=\langle 2\left(UU^{\top}-M^{*}\right)U,Z\rangle.$

By linearity of expectation, we conclude

$\operatorname*{\mathbb{E}}\nabla f_{i}(U)=2(UU^{\top}-M^{*})U.$ (20)

The variance bound can also be obtained via Lemma C.2

$\displaystyle\operatorname{Var}\langle\nabla f_{i}(U),Z\rangle$
$\displaystyle\leq
2\left\|UU^{\top}-M^{*}\right\|_{F}^{2}\left\|ZU^{\top}+UZ^{\top}\right\|_{F}^{2}\leq
8\|UU^{\top}-M^{*}\|_{F}^{2}\left\|U\right\|_{\operatorname{op}}^{2}\|Z\|_{F}^{2},$

where the last inequality comes from Fact A.1.

Since $\operatorname{Var}\langle\nabla
f_{i}(U),Z\rangle=\operatorname{Var}[\operatorname{vec}(\nabla
f_{i}(U))^{\top}\operatorname{vec}(Z)]=\operatorname{vec}(Z)^{\top}\operatorname{Cov}(\operatorname{vec}(\nabla
f_{i}(U)))\operatorname{vec}(Z)$, we have

$\displaystyle\left\|\operatorname{Cov}(\operatorname{vec}(\nabla
f_{i}(U)))\right\|_{\operatorname{op}}$ $\displaystyle\leq
8\left\|UU^{\top}-M^{*}\right\|_{F}^{2}\left\|U\right\|_{\operatorname{op}}^{2}$
$\displaystyle\leq
8\Gamma\left(\left\|UU^{\top}\right\|_{F}+\left\|M^{*}\right\|_{F}\right)^{2}$
$\displaystyle\leq
8\Gamma(r\Gamma+r{\sigma^{\star}_{1}})^{2}=32r^{2}\Gamma^{3}.\qed$

####  C.1.3 Hessian Covariance Bound

######  Lemma C.3 (Equation (7)).

Consider the noiseless setting as in Theorem 1.7 with $\sigma=0$. For all
$U\in\mathbb{R}^{d\times r}$ with
$\left\|U\right\|_{\operatorname{op}}^{2}\leq\Gamma$, it holds that

$\left\|\operatorname{Cov}(\operatorname{vec}(H_{i}))\right\|_{\operatorname{op}}\leq
16r\left\|UU^{\top}-M^{*}\right\|_{F}^{2}+128\left\|U\right\|_{\operatorname{op}}^{4}\leq
192r^{3}\Gamma^{2}.$ (21)

###### Proof.

As a shorthand, we write $B_{i}=A_{i}+A_{i}^{\top}$. Recall that

$H_{i}(U)=\langle UU^{\top}-M^{*},A_{i}\rangle I_{r}\otimes
B_{i}+\operatorname{vec}\left(B_{i}U\right)\operatorname{vec}\left(B_{i}U\right)^{\top}.$

Let $W$ be a $dr\times dr$ matrix that has Frobenius norm 1. We work with
$\langle H_{i}(U),W\rangle$ as we did in the gradient covariance bound
calculation.

By Young’s inequality for $L^{2}$ random variables (i.e.,
$\left\|\mathcal{A}+\mathcal{B}\right\|^{2}_{L^{2}}\leq
2\left\|\mathcal{A}\right\|^{2}_{L^{2}}+2\left\|\mathcal{B}\right\|^{2}_{L^{2}}$),

$\displaystyle\operatorname{Var}\langle H_{i}(U),W\rangle$ $\displaystyle\leq
2\operatorname{Var}\underbrace{[\langle UU^{\top}-M^{*},A_{i}\rangle\langle
I_{r}\otimes
B_{i},W\rangle]}_{\mathcal{A}}+2\operatorname{Var}\underbrace{\left[\langle\operatorname{vec}\left(B_{i}U\right)\operatorname{vec}(B_{i}U)^{\top},W\rangle\right]}_{\mathcal{B}}.$

We first bound the term $\mathcal{A}$. Note $I_{r}\otimes B_{i}$ consists of
$r$ copies of $B_{i}$ in the diagnal

$I_{r}\otimes B_{i}=\left[\begin{array}[]{c|c|c}B_{i}&\cdots&O\\\
\hline\cr\vdots&\ddots&\vdots\\\ \hline\cr O&\cdots&B_{i}\end{array}\right].$

We partition $W$ into $r^{2}$ submatrices of dimension $d\times d$.

$W=\left[\begin{array}[]{c|c|c}W_{11}&\cdots&W_{1r}\\\
\hline\cr\vdots&\ddots&\vdots\\\ \hline\cr
W_{r1}&\cdots&W_{rr}\end{array}\right].$

Then

$\displaystyle\operatorname{Var}[\langle UU^{\top}-M^{*},A_{i}\rangle\langle
I_{r}\otimes B_{i},W\rangle]$ $\displaystyle=\operatorname{Var}[\langle
UU^{\top}-M^{*},A_{i}\rangle\langle B_{i},\sum_{i=1}^{r}W_{ii}\rangle]$
$\displaystyle=\operatorname{Var}[\langle UU^{\top}-M^{*},A_{i}\rangle\langle
A_{i},\sum_{i=1}^{r}(W_{ii}+W_{ii}^{\top})\rangle]$ $\displaystyle\leq
2\left\|UU^{\top}-M^{*}\right\|_{F}^{2}\left\|\sum_{i=1}^{r}(W_{ii}+W_{ii}^{\top})\right\|_{F}^{2},$

where the last inequality uses Lemma C.2.

Write
$a_{i}=\operatorname{vec}(W_{ii})\in\mathbb{R}^{d^{2}},b_{i}=\operatorname{vec}(W_{ii}^{\top}),c_{i}=a_{i}+b_{i}$.
Then $\left\|a_{i}\right\|=\left\|b_{i}\right\|=\left\|W_{ii}\right\|_{F}$ and

$\displaystyle\left\|\sum_{i=1}^{r}(W_{ii}+W_{ii}^{\top})\right\|_{F}^{2}$
$\displaystyle=\left\|\sum_{i=1}^{r}c_{i}\right\|^{2}$
$\displaystyle=\sum_{j=1}^{d^{2}}\left(\sum_{i=1}^{r}c_{ij}\right)^{2}\leq\sum_{j=1}^{d^{2}}r\sum_{i=1}^{r}c_{ij}^{2}\qquad\text{(Cauchy-
Schwarz)}$
$\displaystyle=r\sum_{i=1}^{r}\left\|c_{i}\right\|^{2}=r\sum_{i=1}^{r}\left\|a_{i}+b_{i}\right\|^{2}$
$\displaystyle\leq
2r\left(\sum_{i=1}^{r}(\left\|a_{i}\right\|^{2}+\left\|b_{i}\right\|^{2})\right)$
$\displaystyle=4r\sum_{i=1}^{r}\left\|W_{ii}\right\|_{F}^{2}\leq
4r\left\|W\right\|_{F}^{2}=4r.$

Therefore,

$\displaystyle\operatorname{Var}[\langle UU^{\top}-M^{*},A_{i}\rangle\langle
I_{r}\otimes B_{i},W\rangle]$ $\displaystyle\leq
2\left\|UU^{\top}-M^{*}\right\|_{F}^{2}\left\|\sum_{i=1}^{r}(W_{ii}+W_{ii}^{\top})\right\|_{F}^{2}$
$\displaystyle\leq
8r\left\|UU^{\top}-M^{*}\right\|_{F}^{2}\left\|W\right\|_{F}^{2}=8r\left\|UU^{\top}-M^{*}\right\|_{F}^{2}.$

To bound the term $\mathcal{B}$, we apply Fact A.4 with
$\Sigma=\operatorname{Cov}(\operatorname{vec}(B_{i}U))\in\mathbb{R}^{dr\times
dr}$:

$\displaystyle\operatorname{Var}\left[\langle\operatorname{vec}\left(B_{i}U\right)\operatorname{vec}(B_{i}U)^{\top},W\rangle\right]$
$\displaystyle=\operatorname{Var}\left[\operatorname{vec}(B_{i}U)^{\top}W\operatorname{vec}(B_{i}U)\right]$
$\displaystyle=\operatorname{tr}(W\Sigma(W+W^{\top})\Sigma)$
$\displaystyle=\langle W\Sigma,\Sigma W\rangle+\langle W\Sigma,\Sigma
W^{\top}\rangle$ $\displaystyle\leq\left\|W\Sigma\right\|_{F}\left\|\Sigma
W\right\|_{F}+\left\|W\Sigma\right\|_{F}\left\|\Sigma W^{\top}\right\|_{F}$
$\displaystyle\leq\left\|W\right\|_{F}\left\|\Sigma\right\|_{\operatorname{op}}\left\|\Sigma\right\|_{\operatorname{op}}\left\|W\right\|_{F}+\left\|W\right\|_{F}\left\|\Sigma\right\|_{\operatorname{op}}\left\|\Sigma\right\|_{\operatorname{op}}\left\|W\right\|_{F}$
$\displaystyle=2\left\|W\right\|_{F}^{2}\left\|\Sigma\right\|_{\operatorname{op}}^{2}=2\left\|\Sigma\right\|_{\operatorname{op}}^{2},$

where the last inequality comes from Fact A.1.

To bound $\left\|\Sigma\right\|_{\operatorname{op}}^{2}$, we use Fact A.2
again. Since
$\operatorname{vec}(B_{i}U)=\operatorname{vec}(I_{d}B_{i}U)=(U^{\top}\otimes
I_{d})\operatorname{vec}(B_{i})$, we have

$\operatorname{Cov}(\operatorname{vec}{(B_{i}U)})=(U^{\top}\otimes
I_{d})\operatorname{Cov}(\operatorname{vec}(B_{i}))(U\otimes I_{d}).$ (22)

Let $P\in\mathbb{R}^{d^{2}\times d^{2}}$ be the permutation that maps
$\operatorname{vec}(A_{i})$ to $\operatorname{vec}(A_{i}^{\top})$. Then

$\displaystyle\operatorname{Cov}(\operatorname{vec}(B_{i}))$
$\displaystyle=\operatorname{Cov}(\operatorname{vec}(A_{i})+\operatorname{vec}(A_{i}^{\top}))$
$\displaystyle\leq
2I_{d^{2}}+2\operatorname{Cov}(\operatorname{vec}(A_{i}),\operatorname{vec}(A_{i}^{\top}))$
(Young’s inequality)
$\displaystyle=2I_{d^{2}}+2\operatorname{Cov}(\operatorname{vec}(A_{i}),P\operatorname{vec}(A_{i}))$
$\displaystyle=2I_{d^{2}}+2P\operatorname{Cov}(\operatorname{vec}(A_{i}),\operatorname{vec}(A_{i}))$
$\displaystyle=2I_{d^{2}}+2P.$

Since $P$ is a permutation, we have
$\left\|\operatorname{Cov}(\operatorname{vec}(B_{i}))\right\|_{\operatorname{op}}\leq
4$.

It follows from Equation (22) and Fact A.3 that
$\left\|\operatorname{Cov}(\operatorname{vec}(B_{i}U))\right\|_{\operatorname{op}}\leq
4\left\|U\right\|_{\operatorname{op}}^{2}$ and therefore

$\operatorname{Var}\left[\langle\operatorname{vec}\left(B_{i}U\right)\operatorname{vec}(B_{i}U)^{\top},W\rangle\right]\leq
64\left\|U\right\|_{\operatorname{op}}^{4}\left\|W\right\|_{F}^{2}=64\left\|U\right\|_{\operatorname{op}}^{4}.$

We conclude that for all $W\in\mathbb{R}^{dr\times dr}$,

$\operatorname{Var}\langle H_{i},W\rangle\leq
16r\left\|UU^{\top}-M^{*}\right\|_{F}^{2}+128\left\|U\right\|_{\operatorname{op}}^{4},$

hence

$\left\|\operatorname{Cov}(\operatorname{vec}(H_{i}))\right\|_{\operatorname{op}}\leq
16r\left\|UU^{\top}-M^{*}\right\|_{F}^{2}+128\left\|U\right\|_{\operatorname{op}}^{4}\leq
192r^{3}\Gamma^{2},$

where the last inequality uses $r\geq 1$ and
$\left\|U\right\|_{\operatorname{op}}^{2}\leq\Gamma$. ∎

###  C.2 Global Convergence in Low Rank Matrix Sensing—Omitted Proofs

We firstly verify Assumption 1.4 (ii). We recall the entire assumption below:

######  Assumption C.4 (Assumption 1.4).

There exists a bounded region $\mathcal{B}$ such that the function $f$
satisfies:

  1. (i)

There exists a lower bound $f^{*}>-\infty$ such that for all
$x\in\mathcal{B}$, $f(x,A)\geq f^{*}$ with probability $1$.

  2. (ii)

There exist parameters $L_{D_{g}},L_{D_{H}},B_{D_{g}},B_{D_{H}}$ such that
with high probability over the randomness in $A$, $f(\cdot,A)$ is
$L_{D_{g}}$-gradient Lipschitz and $L_{D_{H}}$-Hessian Lipschitz, and the
$\ell_{2}$-norm of gradient and the Frobenius norm of Hessian of $f(\cdot,A)$
are upper bounded by $B_{D_{g}}$ and $B_{D_{H}}$ respectively.

  3. (iii)

There exist parameters $\sigma_{g},\sigma_{H}>0$ such that

$\left\|\operatorname{Cov}_{A\sim\mathcal{G}}(\nabla
f(x,A))\right\|_{\operatorname{op}}\leq\sigma_{g}$ and
$\left\|\operatorname{Cov}_{A\sim\mathcal{G}}(\operatorname{vec}(\nabla^{2}f(x,A)))\right\|_{\operatorname{op}}\leq\sigma_{H}$.

###### Verifying Assumption 1.4(ii).

We take the bounded region $\mathcal{B}$ to be the region
$\\{U:\left\|U\right\|_{\operatorname{op}}^{2}\leq\Gamma\\}$. Let
$f_{i}(U)=\frac{1}{2}\left(\langle UU^{\top},A_{i}\rangle-y_{i}\right)^{2}$ be
the cost function corresponding to clean samples. We need to verify that
$\nabla f_{i}(U)$ and $\nabla^{2}f_{i}(U)$ are
$\operatorname{poly}(dr\Gamma/\epsilon)$-Lipschitz and bounded within
$\mathcal{B}$ with high probability so that the violated samples constitute at
most an $\epsilon$-fraction.

We discuss gradient-Lipschitzness as an example, focusing on the constant
$L_{D_{g}}$. The rest of conditions can be checked similarly. We drop the
subscript $L:=L_{D_{g}}$ in the following analysis for conciseness of the
notation.

Since $\nabla f_{i}(X)-\nabla f_{i}(Y)=\langle
XX^{\top}-M^{*},A_{i}\rangle(A_{i}+A_{i}^{\top})(X-Y)+\langle
XX^{\top}-YY^{\top},A_{i}\rangle(A_{i}+A_{i}^{\top})Y$, it suffices that
$\left\|A_{i}\right\|_{F}^{2}\leq 2d^{2}+3/\epsilon$ for $\nabla f_{i}(\cdot)$
to be $\operatorname{poly}(dr\Gamma/\epsilon)$-Lipschitz. Since
$\left\|A_{i}\right\|_{F}^{2}$ follows chi-square distribution with degree of
freedom $d^{2}$, by Laurent-Massart bound, we have

$\mathbb{P}\\{\left\|A_{i}\right\|_{F}^{2}-d^{2}\geq
d^{2}+3/\epsilon\\}\leq\exp(-1/\epsilon).$

This completes the verification of Assumption 1.4(ii) for $L_{D_{g}}$. We
proceed to discuss why this high probability result implies that at most an
$\epsilon$-fraction of clean samples violates gradient Lipschitzness. By
Chernoff’s inequality, the probability that more than $\epsilon$-fraction of
uncorrupted $A_{i}$’s fail to satisfy $\left\|A_{i}\right\|_{F}^{2}\leq
2d^{2}+3/\epsilon$ is less than

$\exp\left(-n\left((1-\epsilon)\log\frac{1-\epsilon}{1-\exp(-1/\epsilon)}+\epsilon\log\frac{\epsilon}{\exp(-1/\epsilon)}\right)\right)=O\\!\left(\frac{1}{\epsilon}\exp(-n)\right),$

which is a small if $n=\widetilde{O}\left(\frac{1}{\epsilon}\right)$ is
sufficiently large. ∎

######  Corollary C.5.

Consider the noiseless setting as in Theorem 1.7 with $\sigma=0$. There exists
a sample size $n=\widetilde{O}\left(\frac{d^{2}r^{2}}{\epsilon}\right)$ such
that with $n$ samples, all subroutine calls to $\mathsf{RobustMeanEstimation}$
(Algorithm 3) to estimate population gradients $\nabla\bar{f}(U_{t})$ and
population Hessians $\nabla^{2}\bar{f}(U_{t})$ succeed with high probability.
In light of Equations (6) and (7), this means that there exists a large enough
constant $C_{est}$ such that with high probability, for all iterations $t$, we
have

$\displaystyle\left\|\widetilde{g}_{t}-\nabla\bar{f}(U_{t})\right\|_{F}\leq
2C_{est}\left\|U_{t}U_{t}^{\top}-M^{*}\right\|_{F}\left\|U_{t}\right\|_{\operatorname{op}}$
(23)
$\displaystyle\left\|\widetilde{H}_{t}-\nabla^{2}\bar{f}(U_{t})\right\|_{F}\leq
4r^{1/2}\left\|U_{t}U_{t}^{\top}-M^{*}\right\|_{F}+16\left\|U_{t}\right\|_{\operatorname{op}}^{2}.$
(24)

###### Proof.

This follows directly from the last sentence of Theorem B.3, with
$D=d^{2}r^{2}$ and $\sigma_{g},\sigma_{H}$ given by Equations (6) and (7). ∎

We next prove Lemma 3.7, establishing that all iterates stay inside a bounded
region in which the covariance bounds are valid. The analysis utilizes a
“dissipativity” property [Hal10]), which says that the iterate aligns with the
direction of the gradient when the iterate’s norm is large. For the gradient
step, the gradient will reduce the norm of the iterate. For the negative
curvature step, dissipativity property provides a lower bound on the gradient
when the iterate’s norm is large, but we show that Algorithm 2 only takes a
negative curvature step when this lower bound is violated, therefore the
iterate’s norm must be small and negative curvature steps with a fixed and
small stepsize cannot increase the iterate’s norm by too much.

###### Proof of Lemma 3.7.

Recall Algorithm 2 consists of inexact gradient descent steps of size
$1/L_{g}=\frac{1}{16\Gamma}$ and randomized inexact negative curvature steps
of size

$\frac{2\epsilon_{H}}{L_{H}}=\frac{{\sigma^{\star}_{r}}}{36\Gamma}\leq\frac{\Gamma^{1/2}}{36}.$
(25)

We proceed to use induction to prove the following for Algorithm 2:

  1. 1.

Suppose at step $\tau$ we run the negative curvature step to update the
iterate from $U_{\tau}$ to $U_{\tau+1}$, then
$\left\|U_{\tau}\right\|_{\operatorname{op}}\leq\frac{1}{2}\Gamma^{1/2}$

  2. 2.

$\left\|U_{t}\right\|_{\operatorname{op}}\leq\Gamma^{1/2}$ for all $t\geq 0$.

First we consider the inexact gradient steps. We denote the gradient
inexactness by $e_{t}=\nabla\bar{f}(U_{t})-\widetilde{g}_{t}$. Recall that

$\left\|e_{t}\right\|_{F}\leq
4rC_{est}{\sigma^{\star}_{1}}\left\|U_{t}\right\|_{\operatorname{op}}\sqrt{\epsilon}$

according to Equation (23) in Lemma C.5 whenever
$\left\|U_{t}\right\|_{\operatorname{op}}\leq\Gamma^{1/2}$, which is true
according to the induction hypothesis. Therefore, we have

$\displaystyle\left\|U_{t+1}\right\|_{\operatorname{op}}$
$\displaystyle=\left\|U_{t}-\frac{1}{16\Gamma}\widetilde{g}_{t}\right\|_{\operatorname{op}}=\left\|U_{t}-\frac{1}{16\Gamma}(2(U_{t}U_{t}-M^{*})U_{t}-e_{t})\right\|_{\operatorname{op}}$
$\displaystyle\leq\left\|U_{t}-\frac{1}{8\Gamma}U_{t}U_{t}^{\top}U_{t}\right\|_{\operatorname{op}}+\frac{1}{8\Gamma}\left\|M^{*}U_{t}+\frac{1}{2}e_{t}\right\|_{\operatorname{op}}$
$\displaystyle\leq\max_{i}\\{\sigma_{i}(U_{t})-\frac{1}{8\Gamma}\sigma_{i}^{3}(U_{t})\\}+\frac{1}{8\Gamma}({\sigma^{\star}_{1}}+2rC_{est}{\sigma^{\star}_{1}}\sqrt{\epsilon})\left\|U_{t}\right\|_{\operatorname{op}}$

Since the function $t\mapsto t-\frac{1}{8\Gamma}t^{3}$ is increasing in
$[0,\sqrt{\frac{8\Gamma}{3}}]$ and
$\left\|U_{t}\right\|_{\operatorname{op}}\leq\Gamma^{1/2}\leq\sqrt{\frac{8\Gamma}{3}}$
by induction hypothesis, the maximum is taken when $i=1$, hence

$\displaystyle\left\|U_{t+1}\right\|_{\operatorname{op}}$
$\displaystyle\leq\left\|U_{t}\right\|_{\operatorname{op}}-\frac{1}{8\Gamma}\left\|U_{t}\right\|_{\operatorname{op}}^{3}+\frac{1}{8\Gamma}({\sigma^{\star}_{1}}+2rC_{est}{\sigma^{\star}_{1}}\sqrt{\epsilon})\left\|U_{t}\right\|_{\operatorname{op}}$
$\displaystyle\leq\left\|U_{t}\right\|_{\operatorname{op}}-\frac{1}{8\Gamma}(\left\|U_{t}\right\|_{\operatorname{op}}^{2}-{\sigma^{\star}_{1}}-2rC_{est}{\sigma^{\star}_{1}}\sqrt{\epsilon})\left\|U_{t}\right\|_{\operatorname{op}}$
$\displaystyle\leq\left\|U_{t}\right\|_{\operatorname{op}}-\frac{1}{8\Gamma}(\left\|U_{t}\right\|_{\operatorname{op}}^{2}-2{\sigma^{\star}_{1}})\left\|U_{t}\right\|_{\operatorname{op}},$
(26)

where the last inequality uses
$\epsilon=O\\!\left(\frac{1}{\kappa^{3}r^{3}}\right)=O\\!\left(\frac{1}{r^{2}}\right)$.
We split into two cases now:

Case $\left\|U_{t}\right\|_{\operatorname{op}}>\frac{1}{2}\Gamma^{1/2}$.
Recall that $\Gamma\geq 36{\sigma^{\star}_{1}}$, hence

$\displaystyle\left\|U_{t+1}\right\|_{\operatorname{op}}$
$\displaystyle\leq\left\|U_{t}\right\|_{\operatorname{op}}-\frac{1}{8\Gamma}(\left\|U_{t}\right\|_{\operatorname{op}}^{2}-2{\sigma^{\star}_{1}})\left\|U_{t}\right\|_{\operatorname{op}}$
$\displaystyle\leq\left\|U_{t}\right\|_{\operatorname{op}}-\frac{1}{8\Gamma}(\frac{\Gamma}{4}-2{\sigma^{\star}_{1}})\frac{\Gamma^{1/2}}{2}$
$\displaystyle\leq\left\|U_{t}\right\|_{\operatorname{op}}-\frac{1}{8\Gamma}(\frac{\Gamma}{4}-\frac{\Gamma}{18})\frac{\Gamma^{1/2}}{2}$
$\displaystyle\leq\left\|U_{t}\right\|_{\operatorname{op}}-\frac{1}{96}\Gamma^{1/2},$

where the second inequality is because
$t\mapsto(t^{2}-2{\sigma^{\star}_{1}})t$ is increasing on
$\left[\sqrt{\frac{2{\sigma^{\star}_{1}}}{3}},+\infty\right]$ and
$\left\|U\right\|_{\operatorname{op}}\geq\frac{1}{2}\Gamma^{1/2}\geq
3{\sigma^{\star}_{1}}$. This case captures the dissipativity condition
satisfied by the matrix sensing problem, which says that the iterate aligns
with the direction of the gradient when the iterate’s norm is large, so
descending along the gradient decreases the norm of the iterate.

Case $\left\|U_{t}\right\|_{\operatorname{op}}\leq\frac{1}{2}\Gamma^{1/2}$.
From Equation (C.2) we know that if
$\left\|U_{t}\right\|_{\operatorname{op}}^{2}\geq
2\left\|M^{*}\right\|_{\operatorname{op}}$, it always holds that
$\left\|U_{t+1}\right\|_{\operatorname{op}}\leq\left\|U_{t}\right\|_{\operatorname{op}}\leq\frac{1}{2}\Gamma^{1/2}$.
For $\left\|U_{t}\right\|_{\operatorname{op}}^{2}\leq
2\left\|M^{*}\right\|_{\operatorname{op}}$, we have
$\left\|U_{t}\right\|_{\operatorname{op}}\leq\sqrt{2}{\sigma^{\star}_{1}}^{1/2}$,
and therefore

$\displaystyle\left\|U_{t+1}\right\|_{\operatorname{op}}$
$\displaystyle\leq\left\|U_{t}\right\|_{\operatorname{op}}-\frac{1}{8\Gamma}(\left\|U_{t}\right\|_{\operatorname{op}}^{2}-2{\sigma^{\star}_{1}})\left\|U_{t}\right\|_{\operatorname{op}}$
$\displaystyle\leq\left\|U_{t}\right\|_{\operatorname{op}}+\frac{1}{8\Gamma}(\sqrt{2}{\sigma^{\star}_{1}}^{1/2}+\left\|U_{t}\right\|_{\operatorname{op}})((\sqrt{2}{\sigma^{\star}_{1}}^{1/2}-\left\|U_{t}\right\|_{\operatorname{op}}))\left\|U_{t}\right\|_{\operatorname{op}}$
$\displaystyle\leq\left\|U_{t}\right\|_{\operatorname{op}}+\frac{2\sqrt{2}{\sigma^{\star}_{1}}^{1/2}}{8\Gamma}(\sqrt{2}{\sigma^{\star}_{1}}^{1/2}-\left\|U_{t}\right\|_{\operatorname{op}})\sqrt{2}{\sigma^{\star}_{1}}^{1/2}$
$\displaystyle=\left\|U_{t}\right\|_{\operatorname{op}}+\frac{{\sigma^{\star}_{1}}}{2\Gamma}(\sqrt{2}{\sigma^{\star}_{1}}-\left\|U_{t}\right\|_{\operatorname{op}})$
$\displaystyle\leq\sqrt{2}{\sigma^{\star}_{1}}^{1/2}\text{\qquad(because the
above expression is increasing in
$\left\|U_{t}\right\|_{\operatorname{op}}$)}$
$\displaystyle\leq\frac{\sqrt{2}}{6}\Gamma^{1/2}<\Gamma^{1/2},$

where in the second last inequality we use
${\sigma^{\star}_{1}}\leq\Gamma/36$.

In summary, for all $t\in\mathbb{N}$,

$\displaystyle\left\|U_{t}\right\|_{\operatorname{op}}>\frac{1}{2}\Gamma^{1/2}$
$\displaystyle\implies\left\|U_{t}-\frac{1}{16\Gamma}\widetilde{g}_{t}\right\|_{\operatorname{op}}\leq\left\|U_{t}\right\|_{\operatorname{op}}-\frac{1}{96}\Gamma^{1/2}$
(27)
$\displaystyle\left\|U_{t}\right\|_{\operatorname{op}}\leq\frac{1}{2}\Gamma^{1/2}$
$\displaystyle\implies\left\|U_{t}-\frac{1}{16\Gamma}\widetilde{g}_{t}\right\|_{\operatorname{op}}\leq\frac{1}{2}\Gamma^{1/2}$
(28)

In particular, for gradient step regime
$U_{t+1}:=U_{t}-\frac{1}{16\Gamma}\widetilde{g}_{t}$, we know the iterate
$U_{t+1}$ stays inside the region
$\\{U:\left\|U\right\|_{\operatorname{op}}^{2}\leq\Gamma\\}$ provided that
$U_{t}$ is inside the region.

We proceed to analyze negative curvature steps, which only happen if the
inexact gradient is small
$\left\|\widetilde{g}_{t}\right\|_{F}\leq\epsilon_{g}=\frac{1}{32}{\sigma^{\star}_{r}}^{3/2}$.
Note that it follows from Equation (27) that

$\displaystyle\left\|\widetilde{g}_{t}\right\|_{F}\leq\frac{1}{32}{\sigma^{\star}_{r}}^{3/2}$
$\displaystyle\implies\left\|\widetilde{g}_{t}\right\|_{\operatorname{op}}\leq\frac{1}{32}{\sigma^{\star}_{r}}^{3/2}\leq\frac{1}{32}\Gamma^{3/2}$
$\displaystyle\implies\left\|U_{t}-\frac{1}{16\Gamma}\widetilde{g}_{t}\right\|_{\operatorname{op}}\geq\left\|U_{t}\right\|_{\operatorname{op}}-\frac{1}{512}\Gamma^{1/2}>\left\|U_{t}\right\|_{\operatorname{op}}-\frac{1}{96}\Gamma^{1/2}$
$\displaystyle\implies\left\|U_{t}\right\|_{\operatorname{op}}\leq\frac{1}{2}\Gamma^{1/2}.$

Therefore, if at step $\tau$ we run the inexact negative curvature step to
update the iterate from $U_{\tau}$ to $U_{\tau+1}$, then
$\left\|U_{\tau}\right\|_{\operatorname{op}}\leq\frac{1}{2}\Gamma^{1/2}$.
Recall in Equation (25) the negative curvature stepsize
$\frac{2\epsilon_{H}}{L_{H}}\leq\frac{\Gamma^{1/2}}{36}<\frac{\Gamma^{1/2}}{2}$,
so for the negative curvature update

$U_{\tau+1}=U_{\tau}+\frac{2\epsilon_{H}}{L_{H}}\sigma_{\tau}\widetilde{p}_{\tau}$

where $\sigma_{\tau}=\pm 1$ with probability 1/2 and $p_{\tau}$ is a vector
with $\left\|p_{\tau}\right\|_{F}=1$, we have

$\displaystyle\left\|U_{\tau+1}\right\|_{\operatorname{op}}$
$\displaystyle\leq\left\|U_{\tau}\right\|_{\operatorname{op}}+\left\|\frac{2\epsilon_{H}}{L_{H}}\sigma_{\tau}\widetilde{p}_{\tau}\right\|_{\operatorname{op}}$
$\displaystyle\leq\left\|U_{\tau}\right\|_{\operatorname{op}}+\left\|\frac{2\epsilon_{H}}{L_{H}}\sigma_{\tau}\widetilde{p}_{\tau}\right\|_{F}<\frac{1}{2}\Gamma^{1/2}+\frac{1}{2}\Gamma^{1/2}=\Gamma^{1/2}.$

Finally, the initialization satisfies
$\left\|U_{0}\right\|_{\operatorname{op}}^{2}\leq\Gamma$ by assumption, so the
induction is complete. ∎

###  C.3 Local Linear Convergence in Low Rank Matrix Sensing—Omitted Proofs

Recall that Theorem 3.5 gives us a
$(\frac{1}{24}{\sigma^{\star}_{r}}^{3/2},\frac{1}{3}{\sigma^{\star}_{r}})$-approximate
SOSP as the initialization for Algorithm 1. We now explain why approximate
SOSP is useful for the local search in the low rank matrix sensing problems.

######  Definition C.6 (Strict Saddle Property).

Function $f(\cdot)$ is a $(\epsilon_{g},\epsilon_{H},\zeta)$-strict saddle
function if for any $U\in\mathbb{R}^{r\times d}$, at least one of following
properties holds:

  1. (a)

$\left\|\nabla f(U)\right\|>\epsilon_{g}$.

  2. (b)

$\lambda_{\min}(\nabla^{2}f(U))<-\epsilon_{H}$.

  3. (c)

$U$ is $\zeta$-close to $\mathcal{X}^{\star}$ — the set of local minima;
namely, $\operatorname{dist}(U,\mathcal{X}^{\star})\leq\zeta$.

Approximate SOSPs are close to the set of local minima for strict saddle
functions:

######  Proposition C.7.

Let $f$ be a $(\epsilon_{g},\epsilon_{H},\zeta)$-strict saddle function and
$\mathcal{X}^{\star}$ be the set of its local minima. If $U_{SOSP}$ is a
$(\epsilon_{g},\epsilon_{H})$-approximate SOSP of $f$, then
$\operatorname{dist}(U_{SOSP},\mathcal{X}^{\star})\leq\zeta$.

###### Proof.

Since $U_{SOSP}$ is a $(\epsilon_{g},\epsilon_{H})$-approximate SOSP of $f$,
we have $\left\|\nabla f(U)\right\|\leq\epsilon_{g}$ and
$\lambda_{\min}(\nabla^{2}f(U))\geq-\epsilon_{H}$, i.e., (a) and (b) in
Definition C.6 fail to hold. Therefore, (c) holds, i.e.,
$\operatorname{dist}(U,\mathcal{X}^{\star})\leq\zeta$. ∎

Some regularity conditions only hold in a local neighborhood of
$\mathcal{X}^{\star}$. We focus on the following definition.

######  Definition C.8 (Local Regularity Condition).

In a $\zeta$-neighborhood of the set of local minima $\mathcal{X}^{\star}$
(that is, $\operatorname{dist}(U,\mathcal{X}^{\star})\leq\zeta$), the function
$f(\cdot)$ satisfies an $(\alpha,\beta)$-regularity condition if for any $U$
in this neighborhood:

$\langle\nabla
f(U),U-\mathcal{P}_{\mathcal{X}^{\star}}(U)\rangle\geq\frac{\alpha}{2}\left\|U-\mathcal{P}_{\mathcal{X}^{\star}}(U)\right\|^{2}+\frac{1}{2\beta}\left\|\nabla
f(U)\right\|^{2}.$ (29)

Local regularity condition is a weaker condition than strong convexity, but
both conditions would allow a first-order algorithm to obtain local linear
convergence.

Recall that $\bar{f}(U)=\frac{1}{2}\left\|UU^{\top}-M^{*}\right\|_{F}^{2}$.
The global minima of $\bar{f}$ solve the so-called symmetric low rank matrix
factorization problem and the corresponding function value is zero. Let
$TDT^{\top}$ be the singular value decomposition of the real symmetric matrix
$M^{*}$. We observe $U^{*}=TD^{1/2}$ is a global optimum of $\bar{f}$.

The function $\bar{f}$ satisfies the local regularity condition and the strict
saddle property, and all of its local minima are global minima [Jin+17].

######  Fact C.9 ([Jin+17], Lemma 7).

For $\bar{f}$ defined in (4), all local minima are global minima. The set of
global minima is characterized by
$\mathcal{X}^{\star}=\\{U^{*}R|RR^{\top}=R^{\top}R=I\\}$. Furthermore,
$\bar{f}(U)$ satisfies:

  1. 1.

$(\frac{1}{24}({\sigma^{\star}_{r}})^{3/2},\frac{1}{3}{\sigma^{\star}_{r}},\frac{1}{3}({\sigma^{\star}_{r}})^{1/2})$-strict
saddle property, and

  2. 2.

$(\frac{2}{3}{\sigma^{\star}_{r}},10{\sigma^{\star}_{1}})$-regularity
condition in the $\frac{1}{3}({\sigma^{\star}_{r}})^{1/2}$-neighborhood of
$\mathcal{X}^{\star}$ in Frobenius norm.

Similar to Corollary C.5, we need a theorem to say that all calls to
$\mathsf{RobustMeanEstimation}$ (Algorithm 3) are successful. For the local
search (Algorithm 1), only inexact gradient oracles are required. Since
gradients have dimension $dr$, the sample complexity needed for all robust
estimates of gradient to be accurate with high probability is
$\widetilde{O}\left(\frac{dr}{\epsilon}\right)$.

######  Corollary C.10.

Consider the noiseless setting as in Theorem 1.7 with $\sigma=0$. There exists
a sample size $n=\widetilde{O}\left(\frac{dr}{\epsilon}\right)$ such that with
$n$ samples, all subroutine calls to $\mathsf{RobustMeanEstimation}$
(Algorithm 3) to estimate gradients $\nabla\bar{f}(U_{t})$ succeed with high
probability. In light of Equation (6), this means that there exists a large
enough constant $C_{est}$ such that with high probability which is suppresed
by $\widetilde{O}(\cdot)$, for all iterations $t$, we have

$\displaystyle\left\|\widetilde{g}_{t}-\nabla\bar{f}(U_{t})\right\|_{F}\leq
2C_{est}\left\|U_{t}U_{t}^{\top}-M^{*}\right\|_{F}\left\|U_{t}\right\|_{\operatorname{op}}.$
(30)

###### Proof.

The proof follows exactly the same line of argument as Corollary C.5 and
Theorem B.3, which discusses the sample complexity required for robust mean
estimation of both gradients and Hessians to be successful. Since the Hessian
has dimension $d^{2}r^{2}$, the sample complexity in Corollary C.5 is
dominated by $\widetilde{\Omega}\\!\left(d^{2}r^{2}/\epsilon\right)$. Here the
gradient has dimension $dr$, so we have the required sample complexity
$\widetilde{\Omega}\\!\left(dr/\epsilon\right)$. ∎

###### Proof of Theorem 3.8.

Let $e_{t}=\widetilde{g}_{t}-\nabla\bar{f}(U_{t})$ and we rewrite the inexact
gradient descent update as

$U_{t+1}=U_{t}-\eta(\nabla\bar{f}(U_{t})+e_{t}).$

It follows from Corollary C.10 that, with high probability and for some large
constant $C_{est}$, all calls to $\mathbf{RobustMeanEstimation}$ are
successful, i.e., for all iterations $t$,

$\displaystyle\left\|e_{t}\right\|_{F}\leq
2C_{est}\Gamma^{1/2}\left\|U_{t}U_{t}^{\top}-M^{*}\right\|_{F}\sqrt{\epsilon}.$
(31)

We condition on this event for the rest of the analysis.

It is straightforward to check that
$\frac{1}{3}{\sigma^{\star}_{r}}^{1/2}$-neighborhood of $\mathcal{X}^{\star}$
in Frobenius norm is inside the region
$\\{U:\left\|U\right\|_{\operatorname{op}}^{2}\leq\Gamma\\}$. By strict saddle
property in Fact C.9, we know from Proposition C.7 that $U_{0}=U_{SOSP}$ is
$\frac{1}{3}{\sigma^{\star}_{r}}^{1/2}$-close to $\mathcal{X}^{\star}$ in
Frobenius norm. We assume for the sake of induction that $U_{t}$ is
$\frac{1}{3}{\sigma^{\star}_{r}}^{1/2}$-close to $\mathcal{X}^{\star}$.

Let $\mathcal{P}_{\mathcal{X}^{\star}}(U)$ be the Frobenius projection of
$U\in\mathbb{R}^{d\times r}$ onto $\mathcal{X}^{\star}$. By Fact C.9, for a
given $U\in\mathbb{R}^{d\times r}$, there exists a rotation $R_{U}$ so that
$\mathcal{P}_{\mathcal{X}^{\star}}(U)=U^{*}R_{U}$. Therefore,
$M^{*}=U^{*}{U^{*}}^{\top}=U^{*}R_{U}R_{U}^{\top}{U^{*}}^{\top}=\mathcal{P}_{\mathcal{X}^{\star}}(U)\mathcal{P}_{\mathcal{X}^{\star}}(U)^{\top}$.
We have

$\displaystyle UU^{\top}-M^{*}$
$\displaystyle=UU^{\top}-U\mathcal{P}_{\mathcal{X}^{\star}}(U)^{\top}+U\mathcal{P}_{\mathcal{X}^{\star}}(U)^{\top}-\mathcal{P}_{\mathcal{X}^{\star}}(U)\mathcal{P}_{\mathcal{X}^{\star}}(U)^{\top}$
$\displaystyle=U(U-\mathcal{P}_{\mathcal{X}^{\star}}(U))^{\top}+(U-\mathcal{P}_{\mathcal{X}^{\star}}(U))\mathcal{P}_{\mathcal{X}^{\star}}(U)^{\top}.$

Since $\left\|U_{t}\right\|_{\operatorname{op}}\leq\Gamma^{1/2}$ and
$\left\|\mathcal{P}_{\mathcal{X}^{\star}}(U)\right\|_{\operatorname{op}}\leq\Gamma^{1/2}$,
it follows from (31) that

$\left\|e_{t}\right\|_{F}\leq
2C_{est}\Gamma^{1/2}\sqrt{\epsilon}\left\|U_{t}U_{t}^{\top}-M^{*}\right\|_{F}\leq
4C_{est}\Gamma\sqrt{\epsilon}\left\|U_{t}-\mathcal{P}_{\mathcal{X}^{\star}}(U_{t})\right\|_{F}.$

We proceed to derive the exponential decrease of the distance to global
minima.