Nonconvex Optimization Meets Low-Rank Matrix Factorization: An Overview

Introduction

Here, $f(\cdot)$ is a certain empirical risk function (e.g. Euclidean loss, negative log-likelihood) that evaluates how well a candidate solution fits the observations, and the set $\mathcal{C}$ encodes additional prior constraints, if any. This problem is often highly nonconvex and appears daunting to solve to global optimality at first sight. After all, conventional wisdom usually perceives nonconvex optimization as a computationally intractable task that is susceptible to local minima.

To bypass the challenge, one can resort to convex relaxation, an effective strategy that already enjoys theoretical success in addressing a large number of problems. The basic idea is to convexify the problem by, amongst others, dropping or replacing the low-rank constraint (1b) by a nuclear norm constraint , and solving the convexified problem in the full matrix space (i.e. the space of $\bm{M}$). While such convex relaxation schemes exhibit intriguing performance guarantees in several aspects (e.g. near-minimal sample complexity, stability against noise), its computational cost often scales at least cubically in the size of the matrix $\bm{M}$, which often far exceeds the time taken to read the data. In addition, the prohibitive storage complexity associated with the convex relaxation approach presents another hurdle that limits its applicability to large-scale problems.

This overview article focuses on provable low-rank matrix estimation based on nonconvex optimization. This approach operates over the parsimonious factorized representation (1b) and optimizes the nonconvex loss directly over the low-rank factors $\bm{L}$ and $\bm{R}$. The advantage is clear: adopting economical representation of the low-rank matrix results in low storage requirements, affordable per-iteration computational cost, amenability to parallelization, and scalability to large problem size, when performing iterative optimization methods like gradient descent. However, despite its wide use and remarkable performance in practice , the foundational understanding of generic nonconvex optimization is far from mature. It is often unclear whether an optimization algorithm can converge to the desired global solution and, if so, how fast this can be accomplished. For many nonconvex problems, theoretical underpinnings had been lacking until very recently.

2 Nonconvex optimization meets statistical models

Fortunately, despite general intractability, some important nonconvex problems may not be as hard as they seem. For instance, for several low-rank matrix factorization problems, it has been shown that: under proper statistical models, simple first-order methods are guaranteed to succeed in a small number of iterations, achieving low computational and sample complexities simultaneously (e.g. ). The key to enabling guaranteed and scalable computation is to concentrate on problems arising from specific statistical signal estimation tasks, which may exhibit benign structures amenable to computation and help rule out undesired “hard” instances by focusing on the average-case performance. Two messages deserve particular attention when we examine the geometry of associated nonconvex loss functions:

Basin of attraction. For several statistical problems of this kind, there often exists a reasonably large basin of attraction around the global solution, within which an iterative method like gradient descent is guaranteed to be successful and converge fast. Such a basin might exist even when the sample complexity is quite close to the information-theoretic limit .

Benign global landscape. Several problems provably enjoy benign optimization landscape when the sample size is sufficiently large, in the sense that there is no spurious local minima, i.e. all local minima are also global minima, and that the only undesired stationary points are strict saddle points .

These important messages inspire a recent flurry of activities in the design of two contrasting algorithmic approaches:

Two-stage approach. Motivated by the existence of a basin of attraction, a large number of works follow a two-stage paradigm: (1) initialization, which locates an initial guess within the basin; (2) iterative refinement, which successively refines the estimate without leaving the basin. This approach often leads to very efficient algorithms that run in time proportional to that taken to read the data.

Saddle-point escaping algorithms. In the absence of spurious local minima, a key challenge boils down to how to efficiently escape undesired saddle points and find a local minimum, which is the focus of this approach. This approach does not rely on carefully-designed initialization.

The research along these lines highlights the synergy between statistics and optimization in data science and machine learning. The algorithmic choice often needs to properly exploit the underlying statistical models in order to be truly efficient, in terms of both statistical accuracy and computational efficiency.

3 This paper

Understanding the effectiveness of nonconvex optimization is currently among the most active areas of research in machine learning, information and signal processing, optimization and statistics. Many exciting new developments in the last several years have significantly advanced our understanding of this approach for various statistical problems. This article aims to provide a thorough, but by no means exhaustive, technical overview of important recent results in this exciting area, targeting the broader machine learning, signal processing, statistics, and optimization communities.

The rest of this paper is organized as follows. Section 2 reviews some preliminary facts on optimization that are instrumental to understanding the materials in this paper. Section 3 uses a toy (but non-trivial) example (i.e. rank-1 matrix factorization) to illustrate why it is possible to solve a nonconvex problem to global optimality, through both local and global lenses. Section 4 introduces a few canonical statistical estimation problems that will be visited multiple times in the sequel. Section 5 and Section 6 review gradient descent and its many variants as a local refinement procedure, followed by a discussion of other methods in Section 7. Section 8 discusses the spectral method, which is commonly used to provide an initialization within the basin of attraction. Section 9 provides a global landscape analysis, in conjunction with algorithms that work without the need of careful initialization. We conclude the paper in Section 10 with some discussions and remarks. Furthermore, a short note is provided at the end of several sections to cover some historical remarks and provide further pointers.

4 Notations

Additionally, the standard notation $f(n)=O\left(g(n)\right)$ or $f(n)\lesssim g(n)$ means that there exists a constant $c>0$ such that $\left|f(n)\right|\leq c|g(n)|$, $f(n)\gtrsim g(n)$ means that there exists a constant $c>0$ such that $|f(n)|\geq c\left|g(n)\right|$, and $f(n)\asymp g(n)$ means that there exist constants $c_{1},c_{2}>0$ such that $c_{1}|g(n)|\leq|f(n)|\leq c_{2}|g(n)|$.

Preliminaries in optimization theory

We start by reviewing some basic concepts and preliminary facts in optimization theory. For simplicity of presentation, this section focuses on an unconstrained problem

The optimal solution, if it exists, is denoted by

When $f({\bm{x}})$ is strictly convex,Recall that $f({\bm{x}})$ is said to be strictly convex if and only if for any $\lambda\in(0,1)$ and ${\bm{x}},{\bm{y}}\in\mbox{dom}(f)$, one has $f(\lambda{\bm{x}}+(1-\lambda){\bm{y}})<\lambda f({\bm{x}})+(1-\lambda)f({\bm{y}})$ unless $\bm{x}=\bm{y}$. $\bm{x}_{\mathsf{opt}}$ is unique. But it may be non-unique when $f(\cdot)$ is nonconvex.

To solve (2), arguably the simplest method is (vanilla) gradient descent (GD), which follows the update rule

Here, $\eta_{t}$ is the step size or learning rate at the $t$th iteration, and $\bm{x}_{0}$ is the initial point. This method and its variants are widely used in practice, partly due to their simplicity and scalability to large-scale problems.

A central question is when GD converges fast to the global minimum $\bm{x}_{\mathsf{opt}}$. As is well-known in the optimization literature, GD is provably convergent at a linear rate when $f(\cdot)$ is (locally) strongly convex and smooth. Here, an algorithm is said to converge linearly if the error $\|\bm{x}_{t}-\bm{x}_{\mathsf{opt}}\|_{2}$ converges to 0 as a geometric series. To formally state this result, we define two concepts that commonly arise in the optimization literature.

With these definitions in place, we have the following standard result (e.g. ).

Suppose that $f$ is $\alpha$-strongly convex and $\beta$-smooth within a local ball $\mathcal{B}_{\zeta}(\bm{x}_{\mathsf{opt}}):=\left\{\bm{x}:\|\bm{x}-\bm{x}_{\mathsf{opt}}\|_{2}\leq\zeta\right\}$, and that $\bm{x}_{0}\in\mathcal{B}_{\zeta}(\bm{x}_{\mathsf{opt}})$. If $\eta_{t}\equiv{1}/{\beta}$, then GD obeys

The optimality of $\bm{x}_{\mathsf{opt}}$ indicates that $\nabla f(\bm{x}_{\mathsf{opt}})=\bm{0}$, which allows us to rewrite the GD update rule as

where $\bm{x}\left(\tau\right):=\bm{x}_{\mathsf{opt}}+\tau(\bm{x}_{t}-\bm{x}_{\mathsf{opt}})$. Here, the second line arises from the fundamental theorem of calculus [33, Chapter XIII, Theorem 4.2]. If $\bm{x}_{t}\in\mathcal{B}_{\zeta}(\bm{x}_{\mathsf{opt}})$, then it is self-evident that $\bm{x}(\tau)\in\mathcal{B}_{\zeta}(\bm{x}_{\mathsf{opt}})$, which combined with the assumption of Lemma 1 gives

Therefore, as long as $\eta_{t}\leq 1/\beta$ (and hence $\|\eta_{t}\nabla^{2}f(\bm{x}(\tau))\|\leq 1$), we have

This together with the sub-multiplicativity of $\|\cdot\|$ yields

A byproduct is: if $\bm{x}_{t}\in\mathcal{B}_{\zeta}(\bm{x}_{\mathsf{opt}})$, then the next iterate $\bm{x}_{t+1}$ also falls in $\mathcal{B}_{\zeta}(\bm{x}_{\mathsf{opt}})$. Consequently, applying the above argument recursively and recalling the assumption $\bm{x}_{0}\in\mathcal{B}_{\zeta}(\bm{x}_{\mathsf{opt}})$, we see that all GD iterates remain within $\mathcal{B}_{\zeta}(\bm{x}_{\mathsf{opt}})$. Hence, (7) holds true for all $t$. This immediately concludes the proof. ∎

This result essentially implies that: to yield $\varepsilon$-accuracy (in a relative sense), i.e. $\|\bm{x}_{t}-\bm{x}_{\mathsf{opt}}\|_{2}\leq\varepsilon\|\bm{x}_{\mathsf{opt}}\|_{2}$, the number of iterations required for GD — termed the iteration complexity — is at most

if we initialize GD properly such that ${\bm{x}}_{0}$ lies in the local region $\mathcal{B}_{\zeta}(\bm{x}_{\mathsf{opt}})$. In words, the iteration complexity scales linearly with the condition number — the ratio $\beta/\alpha$ of smoothness to strong convexity parameters. As we shall see, for multiple problems considered herein, the radius $\zeta$ of this locally strongly convex and smooth ball $\mathcal{B}_{\zeta}(\bm{x}_{\mathsf{opt}})$ can be reasonably large (e.g. on the same order of $\|\bm{x}_{\mathsf{opt}}\|_{2})$.

2 Convergence under regularity conditions

Another condition that has been extensively employed in the literature is the Regularity Condition (RC) (see e.g. ), which accommodates algorithms beyond vanilla GD, as well as is applicable to possibly nonsmooth functions. Specifically, consider the iterative algorithm

$\bm{g}(\cdot)$ is said to obey the regularity condition $\mathsf{RC}(\mu,\lambda,\zeta)$ for some $\mu,\lambda,\zeta>0$ if

for all $\bm{x}\in\mathcal{B}_{\zeta}(\bm{x}_{\mathsf{opt}})$.

This condition basically implies that at any feasible point $\bm{x}$, the associated negative search direction $\bm{g}(\bm{x})$ is positively correlated with the error $\bm{x}-\bm{x}_{\mathsf{opt}}$, and hence the update rule (8) — in conjunction with a sufficiently small step size — drags the current point closer to the global solution. It follows from the Cauchy-Schwarz inequality that one must have $\mu\lambda\leq 1$.

It is worth noting that this condition does not require $\bm{g}(\cdot)$ to be differentiable. Also, when $\bm{g}(\cdot)=\nabla f(\cdot)$, it does not require $f(\cdot)$ to be convex within $\mathcal{B}_{\zeta}(\bm{x}_{\mathsf{opt}})$; see Fig. 1 for an example. Instead, the regularity condition can be viewed as a combination of smoothness and “one-point strong convexity” (as the condition is stated w.r.t. a single point ${\bm{x}}_{\mathsf{opt}}$) defined as follows

for all $\bm{x}\in\mathcal{B}_{\zeta}(\bm{x}_{\mathsf{opt}})$. To see this, observe that

where the second line follows from an equivalent definition of the smoothness condition [34, Theorem 5.8]. Combining (10) and (11) arrives at the regularity condition with $\mu=1/\beta$ and $\lambda=\alpha$.

Under the regularity condition, the iterative algorithm (8) converges linearly with suitable initialization.

Under $\mathsf{RC}(\mu,\lambda,\zeta)$, the iterative algorithm (8) with $\eta_{t}\equiv\mu$ and $\bm{x}_{0}\in\mathcal{B}_{\zeta}(\bm{x}_{\mathsf{opt}})$ obeys

Assuming that $\bm{x}_{t}\in\mathcal{B}_{\zeta}(\bm{x}_{\mathsf{opt}})$, we can obtain

where the first inequality comes from $\mathsf{RC}(\mu,\lambda,\zeta)$, and the last line arises if $0\leq\eta_{t}\leq\mu$. By setting $\eta_{t}=\mu$, we arrive at

which also shows that $\bm{x}_{t+1}\in\mathcal{B}_{\zeta}(\bm{x}_{\mathsf{opt}})$. The claim then follows by induction under the assumption that $\bm{x}_{0}\in\mathcal{B}_{\zeta}(\bm{x}_{\mathsf{opt}})$. ∎

In view of Lemma 2, the iteration complexity to reach $\varepsilon$-accuracy (i.e. $\|\bm{x}_{t}-\bm{x}_{\mathsf{opt}}\|_{2}\leq\varepsilon\|\bm{x}_{\mathsf{opt}}\|_{2}$) is at most

as long as a suitable initialization is provided.

3 Critical points

An iterative algorithm like gradient descent often converges to one of its fixed points . For gradient descent, the associated fixed points are (first-order) critical points or stationary points of the loss function, defined as follows.

A first-order critical point (stationary point) ${\bm{x}}$ of $f(\cdot)$ is any point that satisfies

Moreover, we call a point ${\bm{x}}$ an $\varepsilon$-first-order critical point, for some $\varepsilon>0$, if it satisfies $\|\nabla f({\bm{x}})\|_{2}\leq\varepsilon$.

A critical point can be a local minimum, a local maximum, or a saddle point of $f(\cdot)$, depending on the curvatures at / surrounding the point. Specifically, denote by $\nabla^{2}f({\bm{x}})$ the Hessian matrix at ${\bm{x}}$, and let $\lambda_{\min}(\nabla^{2}f({\bm{x}}))$ be its minimum eigenvalue. Then for any first-order critical point ${\bm{x}}$:

if $\nabla^{2}f({\bm{x}})\prec\bm{0}$, then ${\bm{x}}$ is a local maximum;

if $\nabla^{2}f({\bm{x}}))\succ\bm{0}$, then ${\bm{x}}$ is a local minimum;

$\lambda_{\min}(\nabla^{2}f({\bm{x}}))=0$, then ${\bm{x}}$ is either a local minimum or a degenerate saddle point;

$\lambda_{\min}(\nabla^{2}f({\bm{x}}))<0$, then ${\bm{x}}$ is a strict saddle point.

Another useful concept is second-order critical points, defined as follows.

A point $\bm{x}$ is said to be a second-order critical point (stationary point) if $\nabla f({\bm{x}})=\bm{0}$ and $\nabla^{2}f({\bm{x}})\succeq\bm{0}$.

Clearly, second-order critical points do not encompass local maxima and strict saddle points, and, as we shall see, are of particular interest for the nonconvex problems considered herein. Since we are interested in minimizing the loss function, we do not distinguish local maxima and strict saddle points.

A warm-up example: rank-1 matrix factorization

For pedagogical reasons, we begin with a self-contained study of a simple nonconvex matrix factorization problem, demonstrating local convergence in the basin of attraction in Section 3.1 and benign global landscape in Section 3.2. The analysis in this section requires only elementary calculations.

The aim is to find the best rank-$1$ approximation of $\bm{M}$. Clearly, this can be posed as the following problem:Here, the preconstant $1/4$ is introduced for the purpose of normalization and does not affect the solution.

We attempt to minimize the nonconvex function $f(\cdot)$ directly in spite of nonconvexity. This problem, though simple, plays a critical role in understanding the success of nonconvex optimization, since several important nonconvex low-rank estimation problems can be regarded as randomized versions or extensions of this problem.

To begin with, we demonstrate that gradient descent, when initialized at a point sufficiently close to the true optimizer (i.e. $\pm\sqrt{\lambda_{1}}\bm{u}_{1}$), is guaranteed to converge fast.

Consider the problem (14). Suppose $\eta_{t}\equiv\frac{1}{4.5\lambda_{1}}$ and $\|\bm{x}_{0}-\sqrt{\lambda_{1}}\bm{u}_{1}\|_{2}\leq\frac{\lambda_{1}-\lambda_{2}}{15\sqrt{\lambda_{1}}}$. Then the GD iterates (4) obey

By symmetry, Theorem 1 continues to hold if $\bm{u}_{1}$ is replaced by $-\bm{u}_{1}$.

In a nutshell, Theorem 1 establishes linear convergence of GD for rank-1 matrix factorization, where the convergence rate largely depends upon the eigen-gap (relative to the largest eigenvalue). This is a “local” result, assuming that a suitable initialization is present in the basin of attraction $\mathcal{B}_{\zeta}(\sqrt{\lambda_{1}}{\bm{u}}_{1})$. Its radius, which is not optimized in this theorem, is given by $\zeta=\frac{\lambda_{1}-\lambda_{2}}{15\lambda_{1}}\|\sqrt{\lambda_{1}}{\bm{u}}_{1}\|_{2}$. This also depends on the relative eigen-gap $(\lambda_{1}-\lambda_{2})/\lambda_{1}$.

The proof mainly consists of showing that $f(\cdot)$ is locally strongly convex and smooth, which allows us to invoke Lemma 1. The gradient and the Hessian of $f(\bm{x})$ are given respectively by

For notational simplicity, let $\bm{\Delta}:=\bm{x}-\sqrt{\lambda_{1}}\bm{u}_{1}$. A little algebra yields that if $\|\bm{\Delta}\|_{2}\leq\frac{\lambda_{1}-\lambda_{2}}{15\sqrt{\lambda_{1}}}$, then

We start with the smoothness condition. The triangle inequality gives

Next, it comes from the definition of $\bm{\Delta}$ and (17) that

Substitution into (16) yields a strong convexity lower bound:

where the last inequality is an immediate consequence of (18).

In summary, for all $\bm{x}$ obeying $\|\bm{\Delta}\|_{2}\leq\frac{\lambda_{1}-\lambda_{2}}{15\sqrt{\lambda_{1}}}$, one has

Applying Lemma 1 establishes the claim. ∎

The question then comes down to whether one can secure an initial guess of this quality. One popular approach is spectral initialization, obtained by computing the leading eigenvector of $\bm{M}$. For this simple problem, this already yields a solution of arbitrary accuracy. As it turns out, such a spectral initialization approach is particularly useful when dealing with noisy and incomplete measurements. We refer the readers to Section 8 for detailed discussions of spectral initialization methods.

2 Global optimization landscape

We then move on to examining the optimization landscape of this simple problem. In particular, what kinds of critical points does $f({\bm{x}})$ have? This is addressed as follows; see also [31, Section 3.3].

Consider the objective function (14). All local minima of $f(\cdot)$ are global optima. The rest of the critical points are either local maxima or strict saddle points.

Recall that $\bm{x}$ is a critical point if and only if $\nabla f({\bm{x}})=({\bm{x}}{\bm{x}}^{\top}-{\bm{M}}){\bm{x}}=\bm{0}$, or equivalently,

As a result, ${\bm{x}}$ is a critical point if either it aligns with an eigenvector of ${\bm{M}}$ or ${\bm{x}}=\bm{0}$. Given that the eigenvectors of $\bm{M}$ obey ${\bm{M}}{\bm{u}}_{i}=\lambda_{i}{\bm{u}}_{i}$, by properly adjusting the scaling, we determine the set of critical points as follows

To further categorize the critical points, we need to examine the associated Hessian matrices as given by (16). Regarding the critical points $\pm\sqrt{\lambda_{k}}{\bm{u}}_{k}$, we have

With regards to the points $\{\pm\sqrt{\lambda_{1}}{\bm{u}}_{1}\}$, one has

and hence they are (equivalent) local minima of $f(\cdot)$;

For the points $\{\pm\sqrt{\lambda_{k}}{\bm{u}}_{k}\}_{k=2}^{n}$, one has

and therefore they are strict saddle points of $f(\cdot)$.

Finally, the critical point at the origin satisfies

and is hence either a local maxima (if $\lambda_{n}>0$) or a strict saddle point (if $\lambda_{n}=0$).

This result reveals the benign geometry of the problem (14) amenable to optimization. All undesired fixed points of gradient descent are strict saddles, which have negative directional curvature and may not be difficult to escape or avoid.

Formulations of a few canonical problems

For an article of this length, it is impossible to cover all nonconvex statistical problems of interest. Instead, we decide to focus on a few concrete and fundamental matrix factorization problems. This section presents formulations of several such examples that will be visited multiple times throughout this article. Unless otherwise noted, the assumptions made in this section (e.g. restricted isometry for matrix sensing, Gaussian design for phase retrieval) will be imposed throughout the rest of the paper.

Clearly, we cannot distinguish $\bm{X}_{\star}$ from $\bm{X}_{\star}\bm{H}$ for any orthonormal matrix $\bm{H}\in\mathcal{O}^{r\times r}$, as they correspond to the same low-rank matrix $\bm{X}_{\star}\bm{X}_{\star}^{\top}=\bm{X}_{\star}\bm{H}\bm{H}^{\top}\bm{X}_{\star}^{\top}$. This simple fact implies that there exist multiple global optima for (20), a phenomenon that holds for most problems discussed herein.

the true low-rank factors, where $\bm{M}_{\star}=\bm{U}_{\star}\bm{\Sigma}_{\star}\bm{V}_{\star}^{\top}$ stands for its singular value decomposition.

A useful property for the sensing operator that enables tractable algorithmic solutions is the Restricted Isometry Property (RIP), which says that the operator preserves approximately the Euclidean norm of the input matrix when restricted to the set of low-rank matrices. More formally:

holds simultaneously for all $\bm{X}$ of rank at most $r$.

As an immediate consequence, the inner product between two low-rank matrices is also nearly preserved if $\mathcal{A}$ satisfies the RIP. Therefore, $\mathcal{A}$ behaves approximately like an isometry when restricting its operations over low-rank matrices.

If an operator $\mathcal{A}$ satisfies the $2r$-RIP with RIP constant $\delta_{2r}<1$, then

holds simultaneously for all $\bm{X}$ and $\bm{Y}$ of rank at most $r$.

Notably, many random sensing designs are known to satisfy the RIP with high probability, with one remarkable example given below.

If the entries of $\bm{A}_{i}$ are i.i.d. Gaussian entries $\mathcal{N}(0,1)$, then $\mathcal{A}$ as defined in (23) satisfies the $r$-RIP with RIP constant $\delta_{r}$ with high probability as soon as $m\gtrsim(n_{1}+n_{2})r/\delta_{r}^{2}$.

2 Phase retrieval and quadratic sensing

This subsumes phase retrieval as a special case, and comes up in applications such as covariance sketching for streaming data , and phase space tomography under the name coherence retrieval . Here, we wish to solve

Clearly, this is equivalent to the matrix sensing problem by taking $\bm{A}_{i}=\bm{a}_{i}\bm{a}_{i}^{\top}$. Here and throughout, we assume i.i.d. Gaussian design as follows, a tractable model that has been extensively studied recently.

Suppose that the design vectors $\bm{a}_{i}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(\bm{0},\bm{I}_{n})$.

3 Matrix completion

The matrix completion problem then boils down to recovering $\bm{M}_{\star}$ from $\mathcal{P}_{\Omega}(\bm{M}_{\star})$ (or equivalently, from the partially observed entries of $\bm{M}_{\star}$) . This arises in numerous scenarios; for instance, in collaborative filtering, we may want to predict the preferences of all users about a collection of movies, based on partially revealed ratings from the users. Throughout this paper, we adopt the following random sampling model:

Each entry is observed independently with probability $0<p\leq 1$, i.e.

For the positive semidefinite case where $\bm{M}_{\star}=\bm{X}_{\star}\bm{X}_{\star}^{\top}$, the task can be cast as solving

When it comes to the more general case where $\bm{M}_{\star}=\bm{L}_{\star}\bm{R}_{\star}^{\top}$, the task boils down to solving

One parameter that plays a crucial role in determining the feasibility of matrix completion is a certain coherence measure, defined as follows .

As shown in , a low-rank matrix cannot be recovered from a highly incomplete set of entries, unless the matrix satisfies the incoherence condition with a small $\mu$.

Throughout this paper, we let $n:=\max\{n_{1},n_{2}\}$ when referring to the matrix completion problem, and set $\kappa=\sigma_{1}(\bm{M}_{\star})/\sigma_{r}(\bm{M}_{\star})$ to be the condition number of $\bm{M}_{\star}$.

4 Blind deconvolution (the subspace model)

To solve this problem, one seeks a solution to

The recovery performance typically depends on an incoherence measure crucial for blind deconvolution.

Let the incoherence parameter $\mu$ of ${\bm{h}}_{\star}$ be the smallest number such that

5 Low-rank and sparse matrix decomposition / robust principal component analysis

The goal is to separate $\bm{M}_{\star}$ and $\bm{S}_{\star}$ from the (possibly partial) entries of $\bm{\Gamma}_{\star}$. This problem is also known as robust principal component analysis , since we can think of it as recovering the low-rank factors of $\bm{M}_{\star}$ when the observed entries are corrupted by sparse outliers (modeled by $\bm{S}_{\star}$). The problem spans numerous applications in computer vision, medical imaging, and surveillance.

Similar to the matrix completion problem, we assume the random sampling model (31), where $\Omega$ is the set of observed entries. In order to make the problem well-posed, we need the incoherence parameter of ${\bm{M}}_{\star}$ as defined in Definition 7 as well, which precludes ${\bm{M}}_{\star}$ from being too spiky. In addition, it is sometimes convenient to introduce the following deterministic condition on the sparsity pattern and sparsity level of $\bm{S}_{\star}$, originally proposed by . Specifically, it is assumed that the non-zero entries of $\bm{S}_{\star}$ are “spread out”, where there are at most a fraction $\alpha$ of non-zeros per row / column. Mathematically, this means:

Local refinement via gradient descent

This section contains extensive discussions of local convergence analysis of gradient descent. GD is perhaps the most basic optimization algorithm, and its practical importance cannot be overstated. Developing fundamental understanding of this algorithm sheds light on the effectiveness of many other iterative algorithms for solving nonconvex problems.

In the sequel, we will first examine what standard GD theory (cf. Lemma 1) yields for matrix factorization problems; see Section 5.1. While the resulting computational guarantees are optimal for nearly isotropic sampling operators, they become highly pessimistic for most of other problems. We diagnose the cause in Section 5.2.1 and isolate an incoherence condition that is crucial to enable fast convergence of GD. Section 5.2.2 discusses how to enforce proper regularization to promote such an incoherence condition, while Section 5.3 illustrates an implicit regularization phenomenon that allows unregularized GD to converge fast as well. We emphasize that generic optimization theory alone yields overly pessimistic convergence bounds; one needs to blend computational and statistical analyses in order to understand the intriguing performance of GD.

To analyze local convergence of GD, a natural strategy is to resort to the standard GD theory in Lemma 1. This requires checking whether strong convexity and smoothness hold locally, as done in Section 3.1. If so, then Lemma 1 yields an upper bound on the iteration complexity. This simple strategy works well when, for example, the sampling operator is nearly isotropic. In the sequel, we use a few examples to illustrate the applicability and potential drawback of this analysis strategy.

For notational simplicity, this subsection focuses on the symmetric case where $\bm{A}_{i}=\bm{A}_{i}^{\top}$. The gradient update rule (4) for this problem reads

where $\mathcal{A}$ is defined in Section 4.1, and $\mathcal{A}^{*}$ is the conjugate operator of $\mathcal{A}$.

When the sensing matrices $[{\bm{A}}_{i}]_{1\leq i\leq m}$ are random and isotropic, (42) can be viewed as a randomized version of the rank-1 matrix factorization problem discussed in Section 3. To see this, consider for instance the case where the ${\bm{A}}_{i}$’s are drawn from the symmetric Gaussian design, i.e. the diagonal entries of $\bm{A}_{i}$ are i.i.d. $\mathcal{N}(0,1)$ and the off-diagonal entries are i.i.d. $\mathcal{N}(0,1/2)$. For any fixed $\bm{x}$, one has

which coincides with the warm-up example (15) by taking $\bm{M}=\bm{x}_{\star}\bm{x}_{\star}^{\top}$. This bodes well for fast local convergence of GD, at least at the population level (i.e. the case when the sample size $m\rightarrow\infty$).

What happens in the finite-sample regime? It turns out that if the sensing operator satisfies the RIP (cf. Definition 6), then $\nabla^{2}f(\cdot)$ does not deviate too much from its population-level counterpart, and hence $f(\cdot)$ remains locally strongly convex and smooth. This in turn allows one to invoke the standard GD theory to establish local linear convergence.

Consider the problem (42), and suppose the operator (23) satisfies $4$-RIP for RIP constant $\delta_{4}\leq{1}/{44}$. If $\|\bm{x}_{0}-\bm{x}_{\star}\|_{2}\leq{\|\bm{x}_{\star}\|_{2}}/{12}$, then GD with $\eta_{t}\equiv 1/(3\|\bm{x}_{\star}\|_{2}^{2})$ obeys

This theorem, which is a deterministic result, is established in Appendix A. An appealing feature is that: it is possible for such RIP to hold as long as the sample size $m$ is on the order of the information-theoretic limits (i.e. $O(n)$), in view of Fact 1. The take-home message is: for highly random and nearly isotropic sampling schemes, local strong convexity and smoothness continue to hold even in the sample-limited regime.

1.2 Measurements that satisfy the RIP (the rank-r𝑟r case)

The rank-1 case is singled out above due to its simplicity. The result certainly goes well beyond the rank-1 case. Again, we focus on the symmetric caseIf the ${\bm{A}}_{i}$’s are asymmetric, the gradient is given by $\nabla f(\bm{X})=\frac{1}{2m}\sum_{i=1}^{m}(\langle\bm{A}_{i},\bm{X}\bm{X}^{\top}\rangle-y_{i})(\bm{A}_{i}+\bm{A}_{i}^{\top})\bm{X}$, although the update rule (45) remains applicable. (20) with $\bm{A}_{i}=\bm{A}_{i}^{\top}$, in which the gradient update rule (4) satisfies

At first, one might imagine that $f(\cdot)$ remains locally strongly convex. This is, unfortunately, not true, as demonstrated by the following example.

Suppose $\mathcal{A}^{*}\mathcal{A}(\cdot)$ is identity, then $f(\cdot)$ reduces to

for two unit vectors $\bm{u},\bm{v}$ obeying $\bm{u}^{\top}\bm{v}=0$ and any $0<\delta<1$. It is straightforward to verify that

which violates convexity. Moreover, this happens even when $\bm{X}$ is arbitrarily close to $\bm{X}_{\star}$ (by taking $\delta\rightarrow 0$).

Fortunately, the above issue can be easily addressed. The key is to recognize that: one can only hope to recover $\bm{X}_{\star}$ up to global orthonormal transformation, unless further constraints are imposed. Hence, a more suitable error metric is

a counterpart of the Euclidean error when accounting for global ambiguity. For notational convenience, we let

Finding $\bm{H}_{\bm{X}}$ is a classical problem called the orthogonal Procrustes problem .

With these metrics in mind, we are ready to generalize the standard GD theory in Lemma 1 and Lemma 2. In what follows, we assume that $\bm{X}_{\star}$ is a global minimizer of $f(\cdot)$, and make the further homogeneity assumption $\nabla f(\bm{X})\bm{H}=\nabla f(\bm{X}\bm{H})$ for any orthonormal matrix $\bm{H}\in\mathcal{O}^{r\times r}$ — a common fact that arises in matrix factorization problems.

If $\eta_{t}\equiv 1/\beta$, then GD with $\bm{X}_{0}\in\mathcal{B}_{\zeta}(\bm{X}_{\star})$ obeys

First, by definition of $\mathsf{dist}(\cdot,\cdot)$ we have

Next, we claim that a modified regularity condition $\mathsf{RC}(1/\beta,\alpha,\zeta)$ (cf. Definition 3) holds for all $t\geq 0$, that is,

If this claim is valid, then it is straightforward to adapt Lemma 2 to obtain

thus establishing the advertised linear convergence. We omit this part for brevity.

The rest of the proof is thus dedicated to justifying (52). First, Taylor’s theorem reveals that

where $\bm{X}(\tau):=\bm{X}_{t}\bm{H}_{t}+\tau(\bm{X}_{\star}-\bm{X}_{t}\bm{H}_{t})$ for some $0\leq\tau\leq 1$. If $\bm{X}_{t}\bm{H}_{t}$ lies within $\mathcal{B}_{\zeta}(\bm{X}_{\star})$, then one has $\bm{X}(\tau)\in\mathcal{B}_{\zeta}(\bm{X}_{\star})$ as well. We can then substitute the condition (51) into (54) to reach

which can be regarded as a modified version of the one point convexity (10). In addition, repeating the argument in (11) yields

which is a consequence of the smoothness assumption. Combining (55) and (56) establishes (52) for the $t$th iteration, provided that $\bm{X}_{t}\bm{H}_{t}\in\mathcal{B}_{\zeta}(\bm{X}_{\star})$. Finally, since the initial point is assumed to fall within $\mathcal{B}_{\zeta}(\bm{X}_{\star})$, we immediately learn from (53) and induction that $\bm{X}_{t}\bm{H}_{t}\in\mathcal{B}_{\zeta}(\bm{X}_{\star})$ for all $t\geq 0$. This concludes the proof. ∎

The condition (51) is a modification of strong convexity to account for global rotation. In particular, it restricts attention to directions of the form $\bm{Z}\bm{H}_{\bm{Z}}-\bm{X}_{\star}$, where one first adjusts the orientation of $\bm{Z}$ to best align with the global minimizer. To confirm that such restriction is sensible, we revisit Example 1. With proper rotation, one has $\bm{Z}\bm{H}_{\bm{Z}}=\bm{X}_{\star}$ and hence

which becomes strictly positive for $\delta\leq 1/3$. In fact, if $\bm{X}$ is sufficiently close to $\bm{X}_{\star}$, then the condition (51) is valid for (46). Details are deferred to Appendix C.

Further, similar to the analysis for the rank-1 case, we can demonstrate that if $\mathcal{A}$ satisfies the RIP for some sufficiently small RIP constant, then $\nabla^{2}f(\cdot)$ is locally not far from $\nabla^{2}f_{\infty}(\cdot)$ in Example 1, meaning that the condition (51) continues to hold for some $\alpha>0$. This leads to the following result. It is assumed that the ground truth $\bm{M}_{\star}$ has condition number $\kappa$.

Consider the problem (20), and suppose the operator (23) satisfies the $6r$-RIP with RIP constant $\delta_{6r}\leq{1}/{10}$. Then there exist some universal constants $c_{0},c_{1}>0$ such that if $\mathsf{dist}^{2}(\bm{X}_{0},\bm{X}_{\star})\leq\sigma_{r}(\bm{M}_{\star})/16$, then GD with $\eta_{t}\equiv c_{0}/\sigma_{1}(\bm{M}_{\star})$ obeys

The algorithm (45) is referred to as Procrustes flow in .

Three implications of Theorem 4 merit particular attention: (1) the quality of the initialization depends on the least singular value of the truth, so as to ensure that the estimation error does not overwhelm any of the important signal direction; (2) the convergence rate becomes a function of the condition number $\kappa$: the better conditioned the truth is, the faster GD converges; (3) when a good initialization is present, GD converges linearly as long as the sample size $m$ is on the order of the information-theoretic limits (i.e. $O(nr)$), in view of Fact 1.

Our discussion continues to hold for the more general case where ${\bm{M}}_{\star}={\bm{L}}_{\star}{\bm{R}}_{\star}^{\top}$, although an extra regularization term has been suggested to balance the size of the two factors. Specifically, we introduce a regularized version of the loss (21) as follows

then the convergence rate for the symmetric case remains valid by replacing ${\bm{X}}_{\star}$ (resp. ${\bm{X}}_{t}$) with ${\bm{X}}_{\star}={\scriptsize\begin{bmatrix}\bm{L}_{\star}\\ \bm{R}_{\star}\end{bmatrix}}$ (resp. ${\bm{X}}_{t}={\scriptsize\begin{bmatrix}\bm{L}_{t}\\ \bm{R}_{t}\end{bmatrix}}$) in the error metric.

1.3 Measurements that do not obey the RIP

There is no shortage of important examples where the sampling operators fail to satisfy the standard RIP at all. For these cases, the standard theory in Lemma 1 (or Lemma 2) either is not directly applicable or leads to pessimistic computational guarantees. This subsection presents a few such cases.

We start with phase retrieval (27), for which the gradient update rule is given by

This algorithm, also dubbed as Wirtinger flow, was first investigated in . The name “Wirtinger flow” stems from the fact that Wirtinger calculus is used to calculate the gradient in the complex-valued case.

Using standard concentration inequalities for random matrices , one derives the following strong convexity and smoothness bounds .While [58, Chapter 15.4.3] presents the bounds only for the complex-valued case, all arguments immediately extend to the real case.

Consider the problem (27). There exist some constants $c_{0},c_{1},c_{2}>0$ such that with probability at least $1-O(n^{-10})$,

holds simultaneously for all $\bm{x}$ obeying $\|\bm{x}-\bm{x}_{\star}\|_{2}\leq c_{1}\|\bm{x}_{\star}\|_{2}$, provided that $m\geq c_{0}n\log n$.

This lemma says that $f(\cdot)$ is locally 0.5-strongly convex and $c_{2}n$-smooth when the sample size $m\gtrsim n\log n$. The sample complexity only exceeds the information-theoretic limit by a logarithmic factor. Applying Lemma 1 then reveals that:

Under the assumptions of Lemma 5, the GD iterates obey

with probability at least $1-O(n^{-10})$, provided that $\eta_{t}\equiv 1/(c_{2}n\|\bm{x}_{\star}\|_{2}^{2})$ and $\|\bm{x}_{0}-\bm{x}_{\star}\|_{2}\leq c_{1}\|\bm{x}_{\star}\|_{2}$.

This is precisely the computational guarantee given in , albeit derived via a different argument. The above iteration complexity bound, however, is not appealing in practice: it requires $O(n\log\frac{1}{\varepsilon})$ iterations to guarantee $\varepsilon$-accuracy (in a relative sense). For large-scale problems where $n$ is very large, such an iteration complexity could be prohibitive.

Phase retrieval is certainly not the only problem where classical results in Lemma 1 yield unsatisfactory answers. The situation is even worse for other important problems like matrix completion and blind deconvolution, where strong convexity (or the modified version accounting for global ambiguity) does not hold at all unless we restrict attention to a constrained class of decision variables. All this calls for new ideas in establishing computational guarantees that match practical performances.

2 Improved computational guarantees via restricted geometry and regularization

As emphasized in the preceding subsection, two issues stand out in the absence of RIP:

The smoothness condition may not be well-controlled;

Local strong convexity may fail, even if we account for global ambiguity.

We discuss how to address these issues, by first identifying a restricted region with amenable geometry for fast convergence (Section 5.2.1) and then applying regularized gradient descent to ensure the iterates stay in the restricted region (Section 5.2.2).

While desired strong convexity and smoothness (or regularity conditions) may fail to hold in the entire local ball, it is possible for them to arise when we restrict ourselves to a small subset of the local ball and / or a set of special directions.

Take phase retrieval for example: $f(\cdot)$ is locally $0.5$-strongly convex, but the smoothness parameter is exceedingly large (see Lemma 5). On closer inspection, those points $\bm{x}$ that are too aligned with any of the sampling vector $\{\bm{a}_{i}\}$ incur ill-conditioned Hessians. For instance, suppose $\bm{x}_{\star}$ is a unit vector independent of $\{\bm{a}_{i}\}$. Then the point $\bm{x}=\bm{x}_{\star}+\delta\frac{\bm{a}_{j}}{\|\bm{a}_{j}\|_{2}}$ for some constant $\delta$ often results in extremely large $\bm{x}^{\top}\nabla^{2}f(\bm{x})\bm{x}$. With high probability, $(\bm{a}_{j}^{\top}\bm{x})^{2}=\left(1-o(1)\right)\delta^{2}n,$ and hence the Hessian (cf. (60)) at this point $\bm{x}$ satisfies $\bm{x}^{\top}\nabla^{2}f(\bm{x})\bm{x}\geq\frac{3}{m}(\bm{a}_{j}^{\top}\bm{x})^{4}-O(\|\frac{1}{m}\sum\nolimits_{j}(\bm{a}_{j}^{\top}\bm{x}^{\star})^{2}\bm{a}_{j}\bm{a}_{j}^{\top}\|)=\left(3-o(1)\right)\delta^{4}{n^{2}}/{m}$, much larger than the strong convexity parameter when $m\ll n^{2}$. This simple instance suggests that: in order to ensure well-conditioned Hessians, one needs to preclude points that are too “coherent” with the sampling vectors, as formalized below.

Under the assumptions of Lemma 5, there exist some constants $c_{0},\cdots,c_{3}>0$ such that if $m\geq c_{0}n\log n$, then with probability at least $1-O(mn^{-10})$,

In words, desired smoothness is guaranteed when considering only points sufficiently near-orthogonal to all sampling vectors. Such a near-orthogonality property will be referred to as “incoherence” between $\bm{x}$ and the sampling vectors.

Going beyond phase retrieval, the notion of incoherence is not only crucial to control smoothness, but also plays a critical role in ensuring local strong convexity (or regularity conditions). A partial list of examples include matrix completion, quadratic sensing, blind deconvolution and demixing, etc . In the sequel, we single out the matrix completion problem to illustrate this fact. The interested readers are referred to for regularity conditions for matrix completion, to for strong convexity and smoothness for quadratic sensing, to (resp. ) for strong convexity and smoothness (resp. regularity condition) for blind deconvolution and demixing.

Consider the problem (32). Suppose that $n^{2}p\geq c_{0}\kappa^{2}\mu rn\log n$ for some large constant $c_{0}>0$. Then with probability $1-O\left(n^{-10}\right)$, the Hessian obeys

for all $\bm{Z}$ (with $\bm{H}_{\bm{Z}}$ defined in (50)) and all $\bm{X}$ satisfying

where $\epsilon\leq c_{1}/\sqrt{\kappa^{3}\mu r\log^{2}n}$ for some constant $c_{1}>0$.This is a simplified version of [26, Lemma 7] where we restrict the descent direction to point to ${\bm{X}}_{\star}$.

This lemma confines attention to the set of points obeying

Given that each observed entry $M_{i,j}$ can be viewed as $\bm{e}_{i}^{\top}\bm{M}\bm{e}_{j}$, the sampling basis relies heavily on the standard basis vectors. As a result, the above lemma is essentially imposing conditions on the incoherence between $\bm{X}$ and the sampling basis.

2.2 Regularized gradient descent

While we have demonstrated favorable geometry within the set of local points satisfying the desired incoherence condition, the challenge remains as to how to ensure the GD iterates fall within this set. A natural strategy is to enforce proper regularization. Several auxiliary regularization procedures have been proposed to explicitly promote the incoherence constraints , in the hope of improving computational guarantees. Specifically, one can regularize the loss function, by adding additional regularization term $G(\cdot)$ to the objective function and designing the GD update rule w.r.t. the regularized problem

with $\lambda>0$ the regularization parameter. For example:

Matrix completion (33): the following regularized loss has been proposed

for some scalars $\alpha_{1},\cdots,\alpha_{4}>0$. There are numerous choices of $G_{0}$; the one suggested by is $G_{0}(z)=\max\{z-1,0\}^{2}$. With suitable learning rates and proper initialization, GD w.r.t. the regularized loss provably yields $\varepsilon$-accuracy in $O(\mathsf{poly}(n)\log\frac{1}{\varepsilon})$ iterations, provided that the sample size $n^{2}p\gtrsim\mu^{2}\kappa^{6}nr^{7}\log n$ .

Blind deconvolution (36): suggest the following regularized loss

with $G_{0}(z):=\max\{z-1,0\}^{2}$. It has been demonstrated that under proper initialization and step size, gradient methods w.r.t. the regularized loss reach $\varepsilon$-accuracy in $O(\mathsf{poly}(m)\log\frac{1}{\varepsilon})$ iterations, provided that the sample size $m\gtrsim{(K+N)\log^{2}m}$ .

In both cases, the regularization terms penalize, among other things, the incoherence measure between the decision variable and the corresponding sampling basis. We note, however, that the regularization terms are often found unnecessary in both theory and practice, and the theoretical guarantees derived in this line of works are also subject to improvements, as unveiled in the next subsection (cf. Section 5.3).

Two other regularization approaches are also worth noting: (1) truncated gradient descent; (2) projected gradient descent. Given that they are extensions of vanilla GD and might enjoy additional benefits, we postpone the discussions of them to Section 6.

3 The phenomenon of implicit regularization

Despite the theoretical success of regularized GD, it is often observed that vanilla GD — in the absence of any regularization — converges geometrically fast in practice. One intriguing fact is this: for the problems mentioned above, GD automatically forces its iterates to stay incoherent with the sampling vectors / matrices, without any need of explicit regularization . This means that with high probability, the entire GD trajectory lies within a nice region that enjoys desired strong convexity and smoothness, thus enabling fast convergence.

In the sequel, we start with phase retrieval as the first example:

Under the assumptions of Lemma 5, the GD iterates with proper initialization (see, e.g., spectral initialization in Section 8.2) and $\eta_{t}\equiv 1/(c_{3}\|\bm{x}_{\star}\|_{2}^{2}\log n)$ obey

for all $t\geq 0$ with probability $1-O(n^{-10})$. Here, $c_{3},c_{4}>0$ are some constants, and we assume $\|\bm{x}_{0}-\bm{x}_{\star}\|_{2}\leq\|\bm{x}_{0}+\bm{x}_{\star}\|_{2}$.

In words, (68b) reveals that all iterates are incoherent w.r.t. the sampling vectors, and hence fall within the nice region characterized in Lemma 6. With this observation in mind, it is shown in (68a) that vanilla GD converges in $O(\log n\log\frac{1}{\varepsilon})$ iterations. This significantly improves upon the computational bound in Theorem 5 derived based on the smoothness property without restricting attention to the incoherence region.

Similarly, for quadratic sensing (29) where the GD update rule is given by

we have the following result, which generalizes Theorem 6 to the low-rank setting.

Consider the problem (29). Suppose the sample size satisfies $m\geq c_{0}nr^{4}\kappa^{3}\log{n}$ for some large constant $c_{0}$, then with probability $1-O(mn^{-10})$, the GD iterates with proper initialization (see, e.g., spectral initialization in Section 8.2) and $\eta_{t}=\eta\equiv 1/(c_{1}(r\kappa+\log{n})^{2}\sigma_{1}({\bm{M}}_{\star}))$ obey

for all $t\geq 0$. Here, $c_{1}>0$ is some absolute constant.

This theorem demonstrates that vanilla GD converges within O\big{(}\max\{r,\log n\}^{2}\log\frac{1}{\varepsilon}\big{)} iterations for quadratic sensing of a rank-$r$ matrix. This significantly improves upon the computational bounds in which do not consider the incoherence region.

The next example is matrix completion (32), for which the GD update rule reads

with $\mathcal{P}_{\Omega}$ defined in (30). The theory for this update rule is:

Consider the problem (32). Suppose that the sample size satisfies $n^{2}p\geq c_{0}\mu^{2}r^{2}n\log n$ for some large constant $c_{0}>0$, and that the condition number $\kappa$ of $\bm{M}_{\star}=\bm{X}_{\star}\bm{X}_{\star}^{\top}$ is a fixed constant. With probability at least $1-O\left(n^{-3}\right)$, the GD iterates (71) with proper initialization (see, e.g., spectral initialization in Section 8.2) satisfy

for all $t\geq 0$, with $\bm{H}_{\bm{X}_{t}}$ defined in (50). Here, $c_{1}>0$ is some constant, and $\eta_{t}\equiv{c_{2}}/(\kappa\sigma_{1}(\bm{M}_{\star}))$ for some constant $c_{2}>0$.

The last example is blind deconvolution. To measure the discrepancy between any $\bm{z}:=$$\begin{bmatrix}\bm{h}\\ \bm{x}\end{bmatrix}$ and $\bm{z}_{\star}:=$$\begin{bmatrix}\bm{h}_{\star}\\ \bm{x}_{\star}\end{bmatrix}$, we define

which accounts for unrecoverable global scaling and phase. The gradient method, also called Wirtinger flow (WF), is

which enjoys the following theoretical support.

Consider the problem (36). Suppose the sample size $m\geq c_{0}\mu^{2}\max\{K,N\}\mathsf{poly}\log m$ for some large constant $c_{0}>0$. Then there is some constant $c_{1}>0$ such that with probability exceeding $1-O(\min\{K,N\}^{-5})$, the iterates (74) with proper initialization (see, e.g., spectral initialization in Section 8.2) and $\eta_{t}\equiv c_{1}$ satisfy

For conciseness, we only state that the estimation error converges in O\big{(}\log\frac{1}{\varepsilon}\big{)} iterations. The incoherence conditions also provably hold across all iterations; see for details. Similar results have been derived for the blind demixing case as well .

Finally, we remark that the desired incoherence conditions cannot be established via generic optimization theory. Rather, these are proved by exploiting delicate statistical properties underlying the models of interest. The key technique is called a “leave-one-out” argument, which is rooted in probability and random matrix theory and finds applications to many problems . The interested readers can find a general recipe using this argument in .

4 Notes

Provably valid two-stage nonconvex algorithms for matrix factorization were pioneered by Keshavan et al. , where the authors studied the spectral method followed by (regularized) gradient descent on Grassmann manifolds. Partly due to the popularity of convex programming, the local refinement stage of received less attention than convex relaxation and spectral methods around that time. A recent work that popularized the gradient methods for matrix factorization is Candès et al. , which provided the first convergence guarantees for gradient descent (or Wirtinger flow) for phase retrieval. Local convergence of (regularized) GD was later established for matrix completion (without resorting to Grassmann manifolds) , matrix sensing , and blind deconvolution under subspace prior . These works were all based on regularity conditions within a local ball. The resulting iteration complexities for phase retrieval, matrix completion, and blind deconvolution were all sub-optimal, which scaled at least linearly with the problem size. Near-optimal computational guarantees were first derived by via a leave-one-out analysis. Notably, all of these works are local results and rely on proper initialization. Later on, GD was shown to converge within a logarithmic number of iterations for phase retrieval, even with random initialization .

Variants of gradient descent

This section introduces several important variants of gradient descent that serve different purposes, including improving computational performance, enforcing additional structures of the estimates, and removing the effects of outliers, amongst others. Due to space limitation, our description of these algorithms cannot be as detailed as that of vanilla GD. Fortunately, many of the insights and analysis techniques introduced in Section 5 are still applicable, which already shed light on how to understand and analyze these variants. In addition, we caution that all of the theory presented herein is developed for the idealistic models described in Section 4, which might sometimes not capture realistic measurement models. Practitioners should perform comprehensive comparisons of these algorithms on real data, before deciding on which one to employ in practice.

Projected gradient descent modifies vanilla GD (4) by adding a projection step in order to enforce additional structures of the iterates, that is

where the constraint set $\mathcal{C}$ can be either convex or nonconvex. For many important sets $\mathcal{C}$ encountered in practice, the projection step can be implemented efficiently, sometimes even with a closed-form solution. There are two common purposes for including a projection step: 1) to enforce the iterates to stay in a region with benign geometry, whose importance has been explained in Section 5.2.1; 2) to encourage additional low-dimensional structures of the iterates that may be available from prior knowledge.

Here, the projection is to ensure the running iterates stay incoherent with the sampling basis, a property that is crucial to guarantee the algorithm descends properly in every iteration (see Section 5.2.1). One notable example serving this purpose is projected GD for matrix completion , where in the positive semidefinite case (i.e. $\bm{M}_{\star}=\bm{X}_{\star}\bm{X}_{\star}^{\top}$), one runs projected GD w.r.t. the loss function $f(\cdot)$ in (32):

where $\eta_{t}$ is the step size and $\mathcal{P}_{\mathcal{C}}$ denote the Euclidean projection onto the set of incoherent matrices:

with ${\bm{X}}_{0}$ being the initialization and $c$ is a predetermined constant (e.g. $c=2$). This projection guarantees that the iterates stay in a nice incoherent region w.r.t. the sampling basis (similar to the one prescribed in Lemma 7), thus achieving fast convergence. Moreover, this projection can be implemented via a row-wise “clipping” operation, given as

for $i=1,2,\ldots,n$. The convergence guarantee for this update rule is given below, which offers slightly different prescriptions in terms of sample complexity and convergence rate from Theorem 8 using vanilla GD.

Suppose that the sample size satisfies $n^{2}p\geq c_{0}\mu^{2}r^{2}n\log n$ for some large constant $c_{0}>0$, and that the condition number $\kappa$ of $\bm{M}_{\star}=\bm{X}_{\star}\bm{X}_{\star}^{\top}$ is a fixed constant. With probability at least $1-O\left(n^{-1}\right)$, the projected GD iterates (77) satisfy

for all $t\geq 0$, provided that $\mathsf{dist}^{2}(\bm{X}_{0},\bm{X}_{\star})\leq c_{3}\sigma_{r}(\bm{M}_{\star})$ and $\eta_{t}\equiv\eta:={c_{2}}/(\mu r\sigma_{1}(\bm{M}_{\star}))$ for some constant $c_{1},c_{2},c_{3}>0$.

This theorem says that projected GD takes $O(\mu r\log\frac{1}{\varepsilon})$ iterations to yield $\varepsilon$-accuracy (in a relative sense).

The results can be extended to the more general asymmetric case by applying similar modifications mentioned in Remark 3; see .

1.2 Projection for incorporating structural priors

In many problems of practical interest, we might be given some prior knowledge about the signal of interest, encoded by a constraint set $\bm{x}_{\star}\in\mathcal{C}$. Therefore, it is natural to apply projection to enforce the desired structural constraints. One such example is sparse phase retrieval , where it is known a priori that $\bm{x}_{\star}$ in (26) is $k$-sparse, where $k\ll n$. If we have prior knowledge about $\|{\bm{x}}_{\star}\|_{1}$, then we can pick the constraint set $\mathcal{C}$ as follows to promote sparsity

Consider the sparse phase retrieval problem where ${\bm{x}}_{\star}$ is $k$-sparse. Suppose that $m\geq c_{1}k\log n$ for some large constant $c_{1}>0$. The projected GD iterates w.r.t. (27) and the constraint set (80) obey

with probability at least $1-O(n^{-1})$, provided that $\eta_{t}\equiv 1/(c_{2}n\|\bm{x}_{\star}\|_{2}^{2})$ and $\|\bm{x}_{0}-\bm{x}_{\star}\|_{2}\leq\|\bm{x}_{\star}\|_{2}/8$.

Another possible projection constraint set for sparse phase retrieval is the (nonconvex) set of $k$-sparse vectors ,

This leads to a hard-thresholding operation, namely, $\mathcal{P}_{\mathcal{C}}({\bm{x}})$ becomes the best $k$-term approximation of $\bm{x}$ (obtained by keeping the $k$ largest entries (in magnitude) of $\bm{x}$ and setting the rest to 0). The readers are referred to for details. See also for a thresholded GD algorithm — which enforces adaptive thresholding rather than projection to promote sparsity — for solving the sparse phase retrieval problem.

We caution, however, that Theorem 11 does not imply that the sample complexity for projected GD (or thresholded GD) is $O(k\log n)$. So far there is no tractable procedure that can provably guarantee a sufficiently good initial point $\bm{x}_{0}$ when $m\lesssim k\log n$ (see a discussion of the spectral initialization method in Section 8.3.3). Rather, all computationally feasible algorithms (both convex and nonconvex) analyzed so far require sample complexity at least on the order of $k^{2}\log n$ under i.i.d. Gaussian designs, unless $k$ is sufficiently large or other structural information is available . All in all, the computational bottleneck for sparse phase retrieval lies in the initialization stage.

2 Truncated gradient descent

Truncated gradient descent proceeds by trimming away a subset of the measurements when forming the descent direction, typically performed adaptively. We can express it as

where $\mathcal{T}$ is an operator that effectively drops samples that bear undesirable influences over the search directions.

There are two common purposes for enforcing a truncation step: (1) to remove samples whose associated design vectors are too coherent with the current iterate , in order to accelerate convergence and improve sample complexity; (2) to remove samples that may be adversarial outliers, in the hope of improving robustness of the algorithm .

We use phase retrieval to illustrate this benefit. All results discussed so far require a sample size that exceeds $m\gtrsim n\log n$. When it comes to the sample-limited regime where $m\asymp n$, there is no guarantee for strong convexity (or regularity condition) to hold. This presents significant challenges for nonconvex methods, in a regime of critical importance for practitioners.

To better understand the challenge, recall the GD rule (59). When $m$ is exceedingly large, the negative gradient concentrates around the population-level gradient, which forms a reliable search direction. However, when $m\asymp n$, the gradient — which depends on 4th moments of $\{\bm{a}_{i}\}$ and is heavy-tailed — may deviate significantly from the mean, thus resulting in unstable search directions.

To stabilize the search directions, one strategy is to trim away those gradient components \{\nabla f_{i}({\bm{x}}_{t}):=\big{(}(\bm{a}_{i}^{\top}\bm{x}_{t})^{2}-y_{i}\big{)}\bm{a}_{i}\bm{a}_{i}^{\top}\bm{x}_{t}\} whose size deviate too much from the typical size. Specifically, the truncation rule proposed in is:Note that the original algorithm proposed in is designed w.r.t. the Poisson loss, although all theory goes through for the current loss.

where $\alpha_{\text{lb}}$, $\alpha_{\text{ub}}$, $\alpha_{h}$ are predetermined thresholds. This trimming rule — called Truncated Wirtinger flow — effectively removes the “heavy tails”, thus leading to much better concentration and hence enhanced performance.

Consider the problem (27). With probability $1-O(n^{-10})$, the iterates (84) obey

for some constant $0<\rho<1$, provided that $m\geq c_{1}n$, $\|\bm{x}_{0}-\bm{x}_{\star}\|_{2}\leq c_{2}\|\bm{x}_{\star}\|_{2}$ and $\eta_{t}\equiv c_{3}/\|\bm{x}_{\star}\|_{2}^{2}$ for some constants $c_{1},c_{2},c_{3}>0$.

for all ${\bm{x}}$ within a local ball around ${\bm{x}}_{\star}$, where $\lambda,\mu\asymp 1$. See for details.

In comparison to vanilla GD, the truncated version provably achieves two benefits:

Optimal sample complexity: given that one needs at least $n$ samples to recover $n$ unknowns, the sample complexity $m\asymp n$ is orderwise optimal;

Optimal computational complexity: truncated WF yields $\varepsilon$ accuracy in O\big{(}\log\frac{1}{\varepsilon}\big{)} iterations. Since each iteration takes time proportional to that taken to read the data, the computational complexity is nearly optimal.

At the same time, this approach is particularly stable in the presence of noise, which enjoys a statistical guarantee that is minimax optimal. The readers are referred to for precise statements.

2.2 Truncation for removing sparse outliers

In many problems, the collected measurements may suffer from corruptions of sparse outliers, and the gradient descent iterates need to be carefully monitored to remove the undesired effects of outliers (which may take arbitrary values). Take robust PCA (40) as an example, in which a fraction $\alpha$ of revealed entries are corrupted by outliers. At the $t$th iterate, one can first try to identify the support of the sparse matrix $\bm{S}_{\star}$ by hard thresholding the residual, namely,

Here, $c>1$ is some predetermined constant (e.g. $c=3$), and the operator $\mathcal{H}_{l}(\cdot)$ is defined as

where $A_{j,\cdot}^{(l)}$ (resp. $A_{\cdot,k}^{(l)}$) denotes the $l$th largest entry (in magnitude) in the $j$th row (resp. column) of $\bm{A}$. The idea is simple: an entry is likely to be an outlier if it is simultaneously among the largest entries in the corresponding row and column. The thresholded residual $\bm{S}_{t+1}$ then becomes our estimate of the sparse outlier matrix $\bm{S}_{\star}$ in the $(t+1)$-th iteration. With this in place, we update the estimate for the low-rank factor by applying projected GD

where $\mathcal{C}$ is the same as (78) to enforce the incoherence condition. This method has the following theoretical guarantee:

Assume that the condition number $\kappa$ of $\bm{M}_{\star}=\bm{X}_{\star}\bm{X}_{\star}^{\top}$ is a fixed constant. Suppose that the sample size and the sparsity of the outlier satisfy $n^{2}p\geq c_{0}\mu^{2}r^{2}n\log n$ and $\alpha\leq c_{1}/(\mu r)$ for some constants $c_{0},c_{1}>0$. With probability at least $1-O\left(n^{-1}\right)$, the iterates satisfy

for all $t\geq 0$, provided that $\mathsf{dist}^{2}(\bm{X}_{0},\bm{X}_{\star})\leq c_{3}\sigma_{r}(\bm{M}_{\star})$. Here, $0<c_{2},c_{3}<1$ are some constants, and $\eta_{t}\equiv{c_{4}}/(\mu r\sigma_{1}(\bm{M}_{\star}))$ for some constant $c_{4}>0$.

In the full data case, the convergence rate can be improved to $1-c_{2}$ for some constant $0<c_{2}<1$.

This theorem essentially says that: as long as the fraction of entries corrupted by outliers does not exceed $O(1/\mu r)$, then the nonconvex algorithm described above provably recovers the true low-rank matrix in about $O(\mu r)$ iterations (up to some logarithmic factor). When $r=O(1)$, it means that the nonconvex algorithm succeeds even when a constant fraction of entries are corrupted.

Another truncation strategy to remove outliers is based on the sample median, as the median is known to be robust against arbitrary outliers . We illustrate this median-truncation approach through an example of robust phase retrieval , where we assume a subset of samples in (26) is corrupted arbitrarily, with their index set denoted by $\mathcal{S}$ with $|\mathcal{S}|=\alpha m$. Mathematically, the measurement model in the presence of outliers is given by

The goal is to still recover $\bm{x}_{\star}$ in the presence of many outliers (e.g. a constant fraction of measurements are outliers).

It is obvious that the original GD iterates (59) are not robust, since the residual

can be perturbed arbitrarily if $i\in\mathcal{S}$. Hence, we instead include only a subset of the samples when forming the search direction, yielding a truncated GD update rule

Here, $\mathcal{T}_{t}$ only includes samples whose residual size $|r_{t,i}|$ does not deviate much from the median of $\{|r_{t,j}|\}_{1\leq j\leq m}$:

where $\mathsf{median}(\cdot)$ denotes the sample median. As the iterates get close to the ground truth, we expect that the residuals of the clean samples will decrease and cluster, while the residuals remain large for outliers. In this situation, the median provides a robust means to tell them apart. One has the following theory, which reveals the success of the median-truncated GD even when a constant fraction of measurements are arbitrarily corrupted.

Consider the problem (89) with a fraction $\alpha$ of arbitrary outliers. There exist some constants $c_{0},c_{1},c_{2}>0$ and $0<\rho<1$ such that if $m\geq c_{0}n\log n$, $\alpha\leq c_{1}$, and $\|\bm{x}_{0}-\bm{x}_{\star}\|_{2}\leq c_{2}\|\bm{x}_{\star}\|_{2}$, then with probability at least $1-O(n^{-1})$, the median-truncated GD iterates satisfy

3 Generalized gradient descent

In all the examples discussed so far, the loss function $f(\cdot)$ has been a smooth function. When $f(\cdot)$ is nonsmooth and non-differentiable, it is possible to continue to apply GD using the generalized gradient (e.g. subgradient) . As an example, consider again the phase retrieval problem but with an alternative loss function, where we minimize the quadratic loss of the amplitude-based measurements, given as

We can simply execute GD w.r.t. the generalized gradient:

Consider the problem (27). There exist some constants $c_{0},c_{1},c_{2}>0$ and $0<\rho<1$ such that if $m\geq c_{0}n$ and $\eta_{t}\equiv c_{2}$, then with high probability,

as long as $\|\bm{x}_{0}-\bm{x}_{\star}\|_{2}\leq\|\bm{x}_{\star}\|_{2}/10$.

4 Projected power method for constrained PCA

Many applications require solving a constrained quadratic maximization (or constrained PCA) problem:

where $\mathcal{C}$ encodes the set of feasible points. This problem becomes nonconvex if either $\bm{L}$ is not negative semidefinite or if $\mathcal{C}$ is nonconvex. To demonstrate the value of studying this problem, we introduce two important examples.

Phase synchronization . Suppose we wish to recover $n$ unknown phases $\phi_{1},\cdots,\phi_{n}\in[0,2\pi]$ given their pairwise relative phases. Alternatively, by setting $(x_{\star})_{i}=e^{\phi_{i}}$, this problem reduces to estimating $\bm{x}_{\star}=[(x_{\star})_{i}]_{1\leq i\leq n}$ from $\bm{x}_{\star}\bm{x}_{\star}^{\mathsf{H}}$ — a matrix that encodes all pairwise phase differences $(x_{\star})_{i}(x_{\star})_{j}^{\mathsf{H}}=e^{\phi_{i}-\phi_{j}}$. To account for the noisy nature of practical measurements, suppose that what we observe is $\bm{L}=\bm{x}_{\star}\bm{x}_{\star}^{\mathsf{H}}+\sigma\bm{W}$, where $\bm{\bm{W}}$ is a Hermitian matrix. Here, $\{W_{i,j}\}_{i\leq j}$ are i.i.d. standard complex Gaussians. The quantity $\sigma$ indicates the noise level, which determines the hardness of the problem. A natural way to attempt recovery is to solve the following problem

Joint alignment . Imagine we want to estimate $n$ discrete variables $\{(x_{\star})_{i}\}_{1\leq i\leq n}$, where each variable can take $m$ possible values, namely, $(x_{\star})_{i}\in\{1,\cdots,m\}$. Suppose that estimation needs to be performed based on pairwise difference samples $y_{i,j}=x_{i}-x_{j}+z_{i,j}~{}\mathsf{mod}~{}m$, where the $z_{i,j}$’s are i.i.d. noise and their distributions dictate the recovery limits. To facilitate computation, one strategy is to lift each discrete variable $x_{i}$ into a $m$-dimensional vector $\bm{x}_{i}\in\{\bm{e}_{1},\cdots,\bm{e}_{m}\}$. We then introduce a matrix $\bm{L}$ that properly encodes all log-likelihood information. After simple manipulation (see for details), maximum likelihood estimation can be cast as follows

More examples of constrained PCA include an alternative formulation of phase retrieval , sparse PCA , and multi-channel blind deconvolution with sparsity priors .

To solve (96), two algorithms naturally come into mind. The first one is projected GD, which follows the update rule

Another possibility is called the projected power method (PPM) , which drops the current iterate and performs projection only over the gradient component:

While this is perhaps best motivated by its connection to the canonical eigenvector problem (which is often solved by the power method), we remark on its close resemblance to projected GD. In fact, for many constrained sets $\mathcal{C}$ (e.g. the ones in phase synchronization and joint alignment), (97) is equivalent to projected GD when the step size $\eta_{t}\rightarrow\infty$.

As it turns out, the PPM provably achieves near-optimal sample and computational complexities for the preceding two examples. Due to the space limitation, the theory is described only in passing.

Phase synchronization. With high probability, the PPM with proper initialization converges linearly to the global optimum, as long as the noise level $\sigma\lesssim\sqrt{n/\log n}$. This is information theoretically optimal up to some log factor .

Joint alignment. With high probability, the PPM with proper initialization converges linearly to the ground truth, as long as certain Kullback-Leibler divergence w.r.t. the noise distribution exceeds the information-theoretic threshold. See details in .

5 Gradient descent on manifolds

In many problems of interest, it is desirable to impose additional constraints on the object of interest, which leads to a constrained optimization problem over manifolds. In the context of low-rank matrix factorization, to eliminate global scaling ambiguity, one might constrain the low-rank factors to live on a Grassmann manifold or a Riemannian quotient manifold .

To fix ideas, take matrix completion as an example. When factorizing ${\bm{M}}_{\star}={\bm{L}}_{\star}{\bm{R}}_{\star}^{\top}$, we might assume ${\bm{L}}_{\star}\in\mathcal{G}(n_{1},r)$, where $\mathcal{G}(n_{1},r)$ denotes the Grassmann manifold which parametrizes all $r$-dimensional linear subspaces of the $n_{1}$-dimensional spaceMore specifically, any point in $\mathcal{G}(n,r)$ is an equivalent class of a $n\times r$ orthonormal matrix. See for details. . In words, we are searching for a $r$-dimensional subspace ${\bm{L}}$ but ignores the global rotation. It is also assumed that ${\bm{L}}_{\star}^{\top}{\bm{L}}_{\star}={\bm{I}}_{r}$ to remove the global scaling ambiguity (otherwise $(c\bm{L},c^{-1}\bm{R})$ is always equivalent to $(\bm{L},\bm{R})$ for any $c\neq 0$). One might then try to minimize the loss function defined over the Grassmann manifold as follows

As it turns out, it is possible to apply GD to $F(\cdot)$ over the Grassmann manifold by moving along the geodesics; here, a geodesic is the shortest path between two points on a manifold. See for an excellent overview. In what follows, we provide a very brief exposure to highlight its difference from a nominal gradient descent in the Euclidean space.

We start by writing out the conventional gradient of $F(\cdot)$ w.r.t. the $t$th iterate ${\bm{L}}_{t}$ in the Euclidean space :

for some properly chosen step size $\eta_{t}$. For the rank-1 case where $r=1$, the update rule (101) can be simplified to

with $\sigma:=\|\nabla_{\mathcal{G}}F({\bm{L}}_{t})\|_{2}$. As can be verified, ${\bm{L}}_{t+1}$ automatically stays on the unit sphere obeying ${\bm{L}}_{t+1}^{\top}{\bm{L}}_{t+1}=1$.

One of the earliest provable nonconvex methods for matrix completion — the $\mathsf{OptSpace}$ algorithm by Keshavan et al. — performs gradient descent on the Grassmann manifold, tailored to the loss function:

where ${\bm{L}}\in\mathcal{G}(n_{1},r)$ and ${\bm{R}}\in\mathcal{G}(n_{2},r)$, with some additional regularization terms to promote incoherence (see Section 5.2.2). It is shown by that GD on the Grassman manifold converges to the truth with high probability if $n^{2}p\gtrsim\mu^{2}\kappa^{6}r^{2}n\log n$, provided that a proper initialization is given.

Other gradient descent approaches on manifolds include . See for an extensive overview of recent developments along this line.

6 Stochastic gradient descent

Many problems have to deal with an empirical loss function that is an average of the sample losses, namely,

When the sample size is large, it is computationally expensive to apply the gradient update rule — which goes through all data samples — in every iteration. Instead, one might apply stochastic gradient descent (SGD) , where in each iteration, only a single sample or a small subset of samples are used to form the search direction. Specifically, the SGD follows the update rule

where $\Omega_{t}\in\{1,\ldots,m\}$ is a subset of cardinality $k$ selected uniformly at random. Here, $k\geq 1$ is known as the mini-batch size. As one can expect, the mini-batch size $k$ plays an important role in the trade-off between the computational cost per iteration and the convergence rate. A properly chosen mini-batch size will optimize the total computational cost given the practical constraints. Please see for the application of SGD in phase retrieval (which has an interesting connection with the Kaczmarz method), and for its application in matrix factorization.

Beyond gradient methods

Gradient descent is certainly not the only method that can be employed to solve the problem (1). Indeed, many other algorithms have been proposed, which come with different levels of theoretical guarantees. Due to the space limitation, this section only reviews two popular alternatives to gradient methods discussed so far. For simplicity, we consider the following unconstrained problem (with slight abuse of notation)

To optimize the core problem (106), alternating minimization (AltMin) alternates between solving the following two subproblems: for $t=1,2\ldots,$

where ${\bm{R}}_{t}$ and ${\bm{L}}_{t}$ are updated sequentially. Here, ${\bm{L}}_{0}$ is an appropriate initialization. For many problems discussed here, both (107a) and (107b) are convex problems and can be solved efficiently.

Consider the loss function (21). In each iteration, AltMin proceeds as follows : for $t=1,2\ldots,$

Each substep consists of a linear least-squares problem, which can often be solved efficiently via the conjugate gradient algorithm . To illustrate why this forms a promising scheme, we look at the following simple example.

Consider the case where $\mathcal{A}$ is identity (i.e. $\mathcal{A}(\bm{M})=\bm{M}$). We claim that given almost any initialization, AltMin converges to the truth after two updates. To see this, we first note that the output of the first iteration can be written as

As long as both ${\bm{L}}_{\star}^{\top}{\bm{L}}_{0}$ and ${\bm{L}}_{0}^{\top}{\bm{L}}_{0}$ are full-rank, the column space of ${\bm{R}}_{\star}$ matches perfectly with that of ${\bm{R}}_{1}$. Armed with this fact, the subsequent least squares problem (i.e. the update for $\bm{L}_{1}$) is exact, in the sense that $\bm{L}_{1}\bm{R}_{1}^{\top}=\bm{M}_{\star}=\bm{L}_{\star}\bm{R}_{\star}^{\top}$.

With the above identity example in mind, we are hopeful that AltMin converges fast if $\mathcal{A}$ is nearly isometric. Towards this, one has the following theory.

Consider the problem (21) and suppose the operator (23) satisfies $2r$-RIP with RIP constant $\delta_{2r}\leq 1/(100\kappa^{2}r)$. If we initialize ${\bm{L}}_{0}$ by the $r$ leading left singular vectors of $\mathcal{A}^{*}({\bm{y}})$, then AltMin achieves

In comparison to the performance of GD in Theorem 4, AltMin enjoys a better iteration complexity w.r.t. the condition number $\kappa$; that is, it obtains $\varepsilon$-accuracy within $O(\log(1/\varepsilon))$ iterations, compared to $O(\kappa\log(1/\varepsilon))$ iterations for GD. In addition, the requirement on the RIP constant depends quadratically on $\kappa$, leading to a sub-optimal sample complexity. To address this issue, Jain et al. further developed a stage-wise AltMin algorithm, which only requires $\delta_{2r}=O(1/r^{2})$. Intuitively, if there is a singular value that is much larger than the remaining ones, then one can treat $\bm{M}_{\star}$ as a (noisy) rank-1 matrix and compute this rank-1 component via AltMin. Following this strategy, one successively applies AltMin to recover the dominant rank-1 component in the residual matrix, unless it is already well-conditioned. See for details.

1.2 Phase retrieval

Consider the phase retrieval problem. It is helpful to think of the amplitude measurements as bilinear measurements of the signs $\bm{b}=\{b_{i}\in\{\pm 1\}\}_{1\leq i\leq m}$ and the signal ${\bm{x}}_{\star}$, namely,

This leads to a simple yet useful alternative formulation for the amplitude loss minimization problem

Therefore, by applying AltMin to the loss function $f(\bm{b},{\bm{x}})$, we obtain the following update rule : for each $t=1,2,\ldots$

where ${\bm{x}}_{0}$ is an appropriate initial estimate, $\bm{A}^{\dagger}$ is the pseudo-inverse of ${\bm{A}}:=[\bm{a}_{1},\cdots,\bm{a}_{m}]^{\top}$, and $\sqrt{{\bm{y}}}:=[\sqrt{y_{i}}]_{1\leq i\leq m}$. The step (110b) can again be efficiently solved using the conjugate gradient method . This is exactly the Error Reduction (ER) algorithm proposed by Gerchberg and Saxton in the 1970s. Given a reasonably good initialization, this algorithm converges linearly under the Gaussian design.

Consider the problem (27). There exist some constants $c_{0},\cdots,c_{4}>0$ and $0<\rho<1$ such that if $m\geq c_{0}n$, then with probability at least $1-c_{2}\exp(-c_{3}m)$, the estimates of AltMin (ER) satisfy

as long as $\|\bm{x}_{0}-\bm{x}_{\star}\|_{2}\leq c_{4}\|\bm{x}_{\star}\|_{2}$.