My stuff
Nonconvex Optimization Meets Low-Rank Matrix Factorization: An Overview

Substantial progress has been made recently on developing provably accurate and efficient algorithms for low-rank matrix factorization via nonconvex optimization. While conventional wisdom often takes a dim view of nonconvex optimization algorithms due to their susceptibility to spurious local minima, simple iterative methods such as gradient descent have been remarkably successful in practice. The theoretical footings, however, had been largely lacking until recently.

In this tutorial-style overview, we highlight the important role of statistical models in enabling efficient nonconvex optimization with performance guarantees. We review two contrasting approaches: (1) two-stage algorithms, which consist of a tailored initialization step followed by successive refinement; and (2) global landscape analysis and initialization-free algorithms. Several canonical matrix factorization problems are discussed, including but not limited to matrix sensing, phase retrieval, matrix completion, blind deconvolution, robust principal component analysis, phase synchronization, and joint alignment. Special care is taken to illustrate the key technical insights underlying their analyses. This article serves as a testament that the integrated consideration of optimization and statistics leads to fruitful research findings.

Modern information processing and machine learning often have to deal with (structured) low-rank matrix factorization. Given a few observations  y ∈ Rm about a matrix  M⋆ ∈ Rn1×n2 of rank r ≪ min{n1, n2}, oneseeks a low-rank solution compatible with this set of observations as well as other prior constraints. Examples include low-rank matrix completion [13], phase retrieval [4], blind deconvolution and self-calibration [5, 6], robust principal component analysis [7, 8], synchronization and alignment [9, 10], to name just a few. A common goal of these problems is to develop reliable, scalable, and robust algorithms to estimate a low-rank matrix of interest, from potentially noisy, nonlinear, and highly incomplete observations.

1.1 Optimization-based methods

Towards this goal, arguably one of the most popular approaches is optimization-based methods. By factorizing a candidate solution  M ∈ Rn1×n2as  LR⊤with low-rank factors  L ∈ Rn1×rand  R ∈ Rn2×r, one attempts recovery by solving an optimization problem in the form of


Here,  f(·)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 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 [13, 1113], and solving the convexified problem in the full matrix space (i.e. the space of 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 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 L and 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 [14, 15], 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.

1.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. [1626]). 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 [1622].

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 [2731].

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.

1.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.

1.4 Notations

It is convenient to introduce a few notations that will be used throughout. We use boldfaced symbols to represent vectors and matrices. For any vector  v, we let ∥v∥2, ∥v∥1 and ∥v∥0denote its  ℓ2, ℓ1, and ℓ0 norm,respectively. For any matrix  M, let ∥M∥, ∥M∥F, ∥M∥∗, ∥M∥2,∞, and ∥M∥∞stand for the spectral norm (i.e. the largest singular value), the Frobenius norm, the nuclear norm (i.e. the sum of singular values), the ℓ2/ℓ∞norm (i.e. the largest  ℓ2norm of the rows), and the entrywise  ℓ∞norm (the largest magnitude of all entries), respectively. We denote by  σj(M) (resp. λj(M)) the jth largest singular value (resp. eigenvalue) of M, and let  Mj,·(resp.  M·,j) represent its jth row (resp. column). The condition number of M is denoted by  κ(M). In addition,  M ⊤, M H and Mindicate the transpose, the conjugate transpose, and the entrywise conjugate of M, respectively. For two matrices A and B of the same size, we define their inner product as ⟨A, B⟩ = Tr(A⊤B), where Tr(·)stands for the trace. The matrix  Indenotes the  n × nidentity matrix, and eidenotes the ith column of  In. For a linear operator A, denote by  A∗ its adjoint operator. For example, if A maps X ∈ Rn×n to [⟨Ai, X⟩]1≤i≤m, then A∗(y) = �mi=1 yiAi. We also let vec(Z) denote the vectorization of a matrix Z. The indicator function  1Aequals 1 when the event A holds true, and 0 otherwise. Further, the notation  Or×r denotes the set of  r × rorthonormal matrices. Let A and B be two square matrices. We write  A ≻ B (resp. A ⪰ B) if their difference  A − Bis a positive definite (resp. positive semidefinite) matrix.

Additionally, the standard notation  f(n) = O (g(n)) or f(n) ≲ g(n)means that there exists a constant c > 0 such that  |f(n)| ≤ c|g(n)|, f(n) ≳ g(n)means that there exists a constant c > 0 such that  |f(n)| ≥ c |g(n)|, and  f(n) ≍ g(n)means that there exist constants  c1, c2 > 0 such that c1|g(n)| ≤ |f(n)| ≤ c2|g(n)|.

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(x) is strictly convex,1 xoptis unique. But it may be non-unique when  f(·)is nonconvex.

2.1 Gradient descent for locally strongly convex functions

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


Here,  ηtis the step size or learning rate at the tth iteration, and  x0is 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  xopt. As is well-known in the optimization literature, GD is provably convergent at a linear rate when  f(·)is (locally) strongly convex and smooth. Here, an algorithm is said to converge linearly if the error  ∥xt − xopt∥2converges to 0 as a geometric series. To formally state this result, we define two concepts that commonly arise in the optimization literature.

Definition 1 (Strong convexity). A twice continuously differentiable function  f : Rn �→ Ris said to be α-strongly convex in a set B if


Definition 2 (Smoothness). A twice continuously differentiable function  f : Rn �→ Ris said to be  β-smoothin a set B if


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

Lemma 1. Suppose that  f is α-strongly convex and  β-smooth within a local ball  Bζ(xopt) := {x : ∥x − xopt∥2 ≤ ζ},and that  x0 ∈ Bζ(xopt). If ηt ≡ 1/β, then GD obeys


Proof of Lemma 1. The optimality of  xoptindicates that  ∇f(xopt) = 0, which allows us to rewrite the GD update rule as


where  x (τ) := xopt + τ(xt − xopt). Here, the second line arises from the fundamental theorem of calculus [33, Chapter XIII, Theorem 4.2]. If  xt ∈ Bζ(xopt), then it is self-evident that  x(τ) ∈ Bζ(xopt), which combined with the assumption of Lemma 1 gives


Therefore, as long as  ηt ≤ 1/β(and hence  ∥ηt∇2f(x(τ))∥ ≤ 1), we have


This together with the sub-multiplicativity of  ∥ · ∥ yields


By setting  ηt = 1/β, we arrive at the desired  ℓ2error contraction, namely,


A byproduct is: if  xt ∈ Bζ(xopt), then the next iterate  xt+1also falls in  Bζ(xopt). Consequently, applying the above argument recursively and recalling the assumption  x0 ∈ Bζ(xopt), we see that all GD iterates remain within  Bζ(xopt). Hence, (7) holds true for all t. This immediately concludes the proof.


Figure 1: An example of  f(·)taken from [20]. Here,  f(x) = x2 if x ∈ [−6, 6] and f(x) = x2 + 1.5|x|(cos(|x| −6) − 1) if |x| > 6. It satisfies  RC(µ, λ, ζ) with some µ, λ > 0 and ζ = ∞, but is clearly nonconvex.

This result essentially implies that: to yield  ε-accuracy (in a relative sense), i.e.  ∥xt − xopt∥2 ≤ ε∥xopt∥2,the number of iterations required for GD — termed the iteration complexity — is at most


if we initialize GD properly such that  x0lies in the local region  Bζ(xopt). In words, the iteration complexity scales linearly with the condition number — the ratio  β/αof smoothness to strong convexity parameters. As we shall see, for multiple problems considered herein, the radius  ζof this locally strongly convex and smooth ball  Bζ(xopt)can be reasonably large (e.g. on the same order of  ∥xopt∥2).

2.2 Convergence under regularity conditions

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


for some general mapping  g(·) : Rn �→ Rn. In vanilla GD,  g(x) = ∇f(x), but  g(·)can also incorporate several variants of GD; see Section 6. The regularity condition is defined as follows.

Definition 3 (Regularity condition). g(·)is said to obey the regularity condition  RC(µ, λ, ζ) for some µ, λ, ζ >


for all  x ∈ Bζ(xopt).

This condition basically implies that at any feasible point x, the associated negative search direction g(x) is positively correlated with the error  x − xopt, 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  µλ ≤ 1.

It is worth noting that this condition does not require  g(·)to be differentiable. Also, when  g(·) = ∇f(·),it does not require  f(·)to be convex within  Bζ(xopt); 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  xopt) defined as follows


for all  x ∈ Bζ(xopt). 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  µ = 1/β and λ = α.Under the regularity condition, the iterative algorithm (8) converges linearly with suitable initialization.

Lemma 2. Under RC(µ, λ, ζ), the iterative algorithm (8) with  ηt ≡ µ and x0 ∈ Bζ(xopt) obeys


Proof of Lemma 2. Assuming that  xt ∈ Bζ(xopt), we can obtain


where the first inequality comes from  RC(µ, λ, ζ), and the last line arises if  0 ≤ ηt ≤ µ. By setting  ηt = µ, wearrive at


which also shows that  xt+1 ∈ Bζ(xopt). The claim then follows by induction under the assumption that x0 ∈ Bζ(xopt).


as long as a suitable initialization is provided.

2.3 Critical points

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

Definition 4 (First-order critical points). A first-order critical point (stationary point)  x of f(·)is any point that satisfies


Moreover, we call a point  x an ε-first-order critical point, for some ε > 0, if it satisfies  ∥∇f(x)∥2 ≤ ε.

A critical point can be a local minimum, a local maximum, or a saddle point of  f(·), depending on the curvatures at / surrounding the point. Specifically, denote by  ∇2f(x)the Hessian matrix at x, and let λmin(∇2f(x))be its minimum eigenvalue. Then for any first-order critical point x:

if  ∇2f(x) ≺ 0, then xis a local maximum;

if  ∇2f(x)) ≻ 0, then xis a local minimum;


Figure 2: Illustration of the function  f(x) = 14∥xx⊤ − M∥2F , where x = [x1, x2]⊤ and M =� 1 −0.5−0.5 1

 λmin(∇2f(x)) = 0, then xis either a local minimum or a degenerate saddle point;

 λmin(∇2f(x)) < 0, then x is a strictsaddle point.

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

Definition 5 (Second-order critical points). A point x is said to be a second-order critical point (stationary point) if  ∇f(x) = 0 and ∇2f(x) ⪰ 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.

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.

Specifically, consider a positive semidefinite matrix  M ∈ Rn×n(which is not necessarily low-rank) with eigenvalue decomposition  M = �ni=1 λiuiu⊤i . We assume throughout this section that there is a gap between the 1st and 2nd largest eigenvalues


The aim is to find the best rank-1 approximation of M. Clearly, this can be posed as the following problem:2


where  f(·)is a degree-four polynomial and highly nonconvex. The solution to (14) can be expressed in closed form as the scaled leading eigenvector  ±√λ1u1. See Fig. 2 for an illustration of the function f(x) when x ∈ R2. This problem stems from interpreting principal component analysis (PCA) from an optimization perspective, which has a long history in the literature of (linear) neural networks and unsupervised learning; see for example [3642].

We attempt to minimize the nonconvex function  f(·)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.

3.1 Local linear convergence of gradient descent

To begin with, we demonstrate that gradient descent, when initialized at a point sufficiently close to the true optimizer (i.e.  ±√λ1u1), is guaranteed to converge fast.

Theorem 1. Consider the problem (14). Suppose  ηt ≡ 14.5λ1and  ∥x0 − √λ1u1∥2 ≤ λ1−λ215√λ1. Then the GD iterates (4) obey


Remark 1. By symmetry, Theorem 1 continues to hold if  u1is replaced by  −u1.

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  Bζ(√λ1u1). Its radius, which is not optimized in this theorem, is given by  ζ = λ1−λ215λ1 ∥√λ1u1∥2. This also depends on the relative eigen-gap  (λ1 − λ2)/λ1.

Proof of Theorem 1. The proof mainly consists of showing that  f(·)is locally strongly convex and smooth, which allows us to invoke Lemma 1. The gradient and the Hessian of f(x) are given respectively by


For notational simplicity, let  ∆ := x − √λ1u1. A little algebra yields that if  ∥∆∥2 ≤ λ1−λ215√λ1 , then


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


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 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.

3.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(x) have? This is addressed as follows; see also [31, Section 3.3].

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

Proof of Theorem 2. Recall that x is a critical point if and only if  ∇f(x) = (xx⊤ −M)x = 0, or equivalently,


As a result, x is a critical point if either it aligns with an eigenvector of M or x = 0. Given that the eigenvectors of M obey  Mui = λiui, 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  ±√λkuk, we have


We can then categorize them as follows:

1. With regards to the points  {±√λ1u1}, one has


2. For the points  {±√λkuk}nk=2, one has


3. Finally, the critical point at the origin satisfies


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.

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.

4.1 Matrix sensing

Suppose that we are given a set of m measurements of  M⋆ ∈ Rn1×n2 of the form


where  {Ai ∈ Rn1×n2}is a collection of sensing matrices known a priori. We are asked to recover  M⋆— which is assumed to be of rank r — from these linear matrix equations [12, 43].

When  M⋆ = X⋆X⊤⋆ ∈ Rn×nis positive semidefinite with  X⋆ ∈ Rn×r, this can be cast as solving the least-squares problem


Clearly, we cannot distinguish  X⋆ from X⋆Hfor any orthonormal matrix  H ∈ Or×r, as they correspond to the same low-rank matrix  X⋆X⊤⋆ = X⋆HH⊤X⊤⋆ . This simple fact implies that there exist multiple global optima for (20), a phenomenon that holds for most problems discussed herein.


Similarly, we cannot distinguish  (L⋆, R⋆)from  (L⋆H, R⋆(H⊤)−1)for any invertible matrix  H ∈ Rr×r, as L⋆R⊤⋆ = L⋆HH−1R⊤⋆ . Throughout the paper, we denote by


the true low-rank factors, where  M⋆ = U⋆Σ⋆V ⊤⋆stands for its singular value decomposition.

In order to make the problem well-posed, we need to make proper assumptions on the sensing operator A : Rn1×n2 �→ Rm defined by:


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:

Definition 6 (Restricted isometry property [44]). An operator  A : Rn1×n2 �→ Rm is said to satisfy the r-RIP with RIP constant  δr < 1 if


holds simultaneously for all X of rank at most r.

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

Lemma 3 ([44]). If an operator A satisfies the 2r-RIP with RIP constant  δ2r < 1, then


holds simultaneously for all X and 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.

Fact 1 (RIP for Gaussian matrices [12, 43]). If the entries of  Aiare i.i.d. Gaussian entries N(0, 1), then A as defined in (23) satisfies the r-RIP with RIP constant  δrwith high probability as soon as  m ≳ (n1 + n2)r/δ2r.

4.2 Phase retrieval and quadratic sensing

Imagine that we have access to m quadratic measurements of a rank-1 matrix  M⋆ := x⋆x⊤⋆ ∈ Rn×n:


where  ai ∈ Rnis the design vector known a priori. How can we reconstruct  x⋆ ∈ Rn— or equivalently, M⋆ = x⋆x⊤⋆ — from this collection of quadratic equations about  x⋆? This problem, often dubbed as phase retrieval, arises for example in X-ray crystallography, where one needs to recover a specimen based on intensities of the diffracted waves scattered by the object [4, 4547]. Mathematically, the problem can be posed as finding a solution to the following program


More generally, consider the quadratic sensing problem, where we collect m quadratic measurements of a rank-r matrix M⋆ := X⋆X⊤⋆ with X⋆ ∈ Rn×r:


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


Clearly, this is equivalent to the matrix sensing problem by taking  Ai = aia⊤i. Here and throughout, we assume i.i.d. Gaussian design as follows, a tractable model that has been extensively studied recently.


4.3 Matrix completion

Suppose we observe partial entries of a low-rank matrix  M⋆ ∈ Rn1×n2of rank r, indexed by the sampling location set  Ω. It is convenient to introduce a projection operator  PΩ : Rn1×n2 �→ Rn1×n2such that for an input matrix  M ∈ Rn1×n2,


The matrix completion problem then boils down to recovering  M⋆ from PΩ(M⋆)(or equivalently, from the partially observed entries of  M⋆) [1, 13]. 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:

Assumption 2 (Random sampling in matrix completion). Each entry is observed independently with probability  0 < p ≤ 1, i.e.


For the positive semidefinite case where  M⋆ = X⋆X⊤⋆ , the task can be cast as solving


When it comes to the more general case where  M⋆ = L⋆R⊤⋆ , 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 [1].

Definition 7 (Incoherence for matrix completion). A rank-r matrix M⋆ ∈ Rn1×n2 with singular value decomposition (SVD)  M⋆ = U⋆Σ⋆V ⊤⋆is said to be  µ-incoherent if


As shown in [13], a low-rank matrix cannot be recovered from a highly incomplete set of entries, unless the matrix satisfies the incoherence condition with a small  µ.

Throughout this paper, we let  n := max{n1, n2}when referring to the matrix completion problem, and set  κ = σ1(M⋆)/σr(M⋆)to be the condition number of  M⋆.

4.4 Blind deconvolution (the subspace model)

Suppose that we want to recover two objects  h⋆ ∈ CKand  x⋆ ∈ CN— or equivalently, the outer product M⋆ = h⋆xH⋆ — from mbilinear measurements of the form


To explain why this is called blind deconvolution, imagine we would like to recover two signals  g ∈ Cmand  d ∈ Cmfrom their circulant convolution [5]. In the frequency domain, the outputs can be written as y = diag(ˆg) ˆd, where  ˆg(resp. ˆd) is the Fourier transform of g (resp. d). If we have additional knowledge that  ˆg = Ax⋆and ˆd = Bh⋆lie in some known subspace characterized by  A = [a1, · · · , am]Hand B = [b1, · · · , bm]H, then yreduces to the bilinear form (35). In this paper, we assume the following semi-random design, a common subspace model studied in the literature [5, 22].

Assumption 3 (Semi-random design in blind deconvolution). Suppose that  aji.i.d.∼ N�0, 12IN�+iN�0, 12IN�,and that  B ∈ Cm×K is formed by the first K columns of a unitary discrete Fourier transform (DFT) matrix F ∈ Cm×m.

To solve this problem, one seeks a solution to


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

Definition 8 (Incoherence for blind deconvolution). Let the incoherence parameter  µof  h⋆be the smallest number such that


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

Suppose we are given a matrix  Γ⋆ ∈ Rn1×n2that is a superposition of a rank-r matrix  M⋆ ∈ Rn1×n2and a sparse matrix  S⋆ ∈ Rn1×n2:


The goal is to separate  M⋆ and S⋆from the (possibly partial) entries of  Γ⋆. This problem is also known as robust principal component analysis [7, 8, 52], since we can think of it as recovering the low-rank factors of  M⋆when the observed entries are corrupted by sparse outliers (modeled by  S⋆). 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  Ωis the set of observed entries. In order to make the problem well-posed, we need the incoherence parameter of  M⋆as defined in Definition 7 as well, which precludes  M⋆from being too spiky. In addition, it is sometimes convenient to introduce the following deterministic condition on the sparsity pattern and sparsity level of  S⋆,originally proposed by [7]. Specifically, it is assumed that the non-zero entries of  S⋆are “spread out”, where there are at most a fraction  αof non-zeros per row / column. Mathematically, this means:

Assumption 4. It is assumed that  S⋆ ∈ Sα ⊆ Rn1×n2, where


For the positive semidefinite case where  M⋆ = X⋆X⊤⋆ with X⋆ ∈ Rn×r, the task can be cast as solving


The general case where  M⋆ = L⋆R⊤⋆ can be formulated similarly by replacing  XX⊤ with LR⊤ in (40) andoptimizing over  L ∈ Rn1×r, R ∈ Rn2×r and S ∈ Sα, namely,


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.

5.1 Computational analysis via strong convexity and smoothness

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.

5.1.1 Measurements that satisfy the RIP (the rank-1 case)

We begin with the matrix sensing problem (19) and consider the case where the truth has rank 1, i.e.  M⋆ =x⋆x⊤⋆ for some vector  x⋆ ∈ Rn. This requires us to solve


For notational simplicity, this subsection focuses on the symmetric case where  Ai = A⊤i . The gradient update rule (4) for this problem reads


where A is defined in Section 4.1, and  A∗ is the conjugate operator of A.

When the sensing matrices  [Ai]1≤i≤mare 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  Ai’s are drawn from the symmetric Gaussian design, i.e. the diagonal entries of  Aiare i.i.d. N(0, 1) and the off-diagonal entries are i.i.d. N(0, 1/2). For any fixed x, one has


which coincides with the warm-up example (15) by taking  M = x⋆x⊤⋆. This bodes well for fast local convergence of GD, at least at the population level (i.e. the case when the sample size  m → ∞).

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

Theorem 3 (GD for matrix sensing (rank-1)). Consider the problem (42), and suppose the operator (23) satisfies 4-RIP for RIP constant  δ4 ≤ 1/44. If ∥x0 − x⋆∥2 ≤ ∥x⋆∥2/12, then GD with  ηt ≡ 1/(3∥x⋆∥22) 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.

5.1.2 Measurements that satisfy the RIP (the rank-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 case3 (20)with  Ai = A⊤i, in which the gradient update rule (4) satisfies


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

Example 1. Suppose A∗A(·)is identity, then  f(·)reduces to


It can be shown that for any  Z ∈ Rn×r (see [26, 29]):


Think of the following example


for two unit vectors  u, v obeying u⊤v = 0 and any 0 < δ < 1. It is straightforward to verify that


which violates convexity. Moreover, this happens even when X is arbitrarily close to  X⋆(by taking  δ → 0).

Fortunately, the above issue can be easily addressed. The key is to recognize that: one can only hope to recover  X⋆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  HXis a classical problem called the orthogonal Procrustes problem [53]. 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  X⋆is a global minimizer of  f(·), and make the further homogeneity assumption  ∇f(X)H = ∇f(XH)for any orthonormal matrix  H ∈ Or×r— a common fact that arises in matrix factorization problems.

Lemma 4. Suppose that  f is β-smooth within a ball  Bζ(X⋆) := {X : ∥X − X⋆∥F ≤ ζ}, and that ∇f(X)H =∇f(XH)for any orthonormal matrix  H ∈ Or×r. Assume that for any  X ∈ Bζ(X⋆) and any Z,


If  ηt ≡ 1/β, then GD with  X0 ∈ Bζ(X⋆) obeys


Proof of Lemma 4. For notational simplicity, let


First, by definition of  dist(·, ·) we have


Next, we claim that a modified regularity condition  RC(1/β, α, ζ)(cf. Definition 3) holds for all  t ≥ 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  X(τ) := XtHt + τ(X⋆ − XtHt)for some  0 ≤ τ ≤ 1. If  XtHtlies within  Bζ(X⋆), then one has X(τ) ∈ Bζ(X⋆)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 tth iteration, provided that  XtHt ∈ Bζ(X⋆). Finally, since the initial point is assumed to fall within  Bζ(X⋆), weimmediately learn from (53) and induction that  XtHt ∈ Bζ(X⋆) for all t ≥ 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  ZHZ − X⋆, where one first adjusts the orientation of Z to best align with the global minimizer. To confirm that such restriction is sensible, we revisit Example 1. With proper rotation, one has  ZHZ = X⋆ and hence


which becomes strictly positive for  δ ≤ 1/3. In fact, if X is sufficiently close to  X⋆, 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 A satisfies the RIP for some sufficiently small RIP constant, then  ∇2f(·)is locally not far from  ∇2f∞(·)in Example 1, meaning that the condition (51) continues to hold for some  α > 0. This leads to the following result. It is assumed that the ground truth  M⋆has condition number  κ.

Theorem 4 (GD for matrix sensing (rank-r) [23, 54]). Consider the problem (20), and suppose the operator (23) satisfies the 6r-RIP with RIP constant  δ6r ≤ 1/10. Then there exist some universal constants  c0, c1 > 0such that if  dist2(X0, X⋆) ≤ σr(M⋆)/16, then GD with  ηt ≡ c0/σ1(M⋆) obeys


Remark 2. The algorithm (45) is referred to as Procrustes flow in [23].

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  κ: 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.

Remark 3 (Asymmetric case). Our discussion continues to hold for the more general case where  M⋆ = L⋆R⊤⋆ ,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


with  λa regularization parameter, e.g.  λ = 1/32as suggested in [23, 55].4Here, the regularization term ��L⊤L − R⊤R��2F is included in order to balance the size of the two low-rank factors L and R. If one applies GD to the regularized loss:


then the convergence rate for the symmetric case remains valid by replacing  X⋆ (resp. Xt) with X⋆ =�L⋆R⋆

(resp.  Xt =�LtRt

5.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 [18]. The name “Wirtinger flow” stems from the fact that Wirtinger calculus is used to calculate the gradient in the complex-valued case.

The associated sampling operator A (cf. (23)), unfortunately, does not satisfy the standard RIP (cf. Defini-tion 6) unless the sample size far exceeds the statistical limit; see [46, 48].5 We can, however, still evaluate the local strong convexity and smoothness parameters to see what computational bounds they produce. Recall that the Hessian of (27) is given by


Using standard concentration inequalities for random matrices [56, 57], one derives the following strong convexity and smoothness bounds [26, 58, 59].6

Lemma 5 (Local strong convexity and smoothness for phase retrieval). Consider the problem (27). There exist some constants  c0, c1, c2 > 0such that with probability at least  1 − O(n−10),


holds simultaneously for all  x obeying ∥x − x⋆∥2 ≤ c1∥x⋆∥2, provided that  m ≥ c0n log n.

This lemma says that  f(·)is locally 0.5-strongly convex and  c2n-smooth when the sample size  m ≳ n log n.The sample complexity only exceeds the information-theoretic limit by a logarithmic factor. Applying Lemma 1 then reveals that:

Theorem 5 (GD for phase retrieval (loose bound) [18]). Under the assumptions of Lemma 5, the GD iterates obey


with probability at least  1 − O(n−10), provided that  ηt ≡ 1/(c2n∥x⋆∥22) and ∥x0 − x⋆∥2 ≤ c1∥x⋆∥2.

This is precisely the computational guarantee given in [18], albeit derived via a different argument. The above iteration complexity bound, however, is not appealing in practice: it requires  O(n log 1ε)iterations to guarantee  ε-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.

5.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).

5.2.1 Restricted strong convexity and smoothness

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(·)is locally 0.5-strongly convex, but the smoothness parameter is exceedingly large (see Lemma 5). On closer inspection, those points x that are too aligned with any of the sampling vector  {ai}incur ill-conditioned Hessians. For instance, suppose  x⋆is a unit vector independent of {ai}. Then the point  x = x⋆ + δ aj∥aj∥2for some constant  δoften results in extremely large  x⊤∇2f(x)x.7This 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.

Lemma 6 (Restricted smoothness for phase retrieval [26]). Under the assumptions of Lemma 5, there exist some constants  c0, · · · , c3 > 0such that if  m ≥ c0n log n, then with probability at least  1 − O(mn−10),


holds simultaneously for all  x ∈ Rn obeying


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 x andthe 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 [19, 21, 22, 26, 6062]. In the sequel, we single out the matrix completion problem to illustrate this fact. The interested readers are referred to [19, 21, 54] for regularity conditions for matrix completion, to [60] for strong convexity and smoothness for quadratic sensing, to [26, 62] (resp. [22, 61]) for strong convexity and smoothness (resp. regularity condition) for blind deconvolution and demixing.

Lemma 7 (Restricted strong convexity and smoothness for matrix completion [26]). Consider the problem (32). Suppose that  n2p ≥ c0κ2µrn log nfor some large constant  c0 > 0. Then with probability  1 − O�n−10�, the Hessian obeys

vec (ZHZ − X⋆)⊤ ∇2f (X) vec (ZHZ − X⋆) ≥ 0.5σr(M⋆) ∥ZHZ − X⋆∥2F (64a)��∇2f (X)�� ≤ 2.5σ1(M⋆) (64b)

for all  Z (with HZ defined in (50)) and all Xsatisfying


where  ϵ ≤ c1/�κ3µr log2 nfor some constant  c1 > 0.8

This lemma confines attention to the set of points obeying


Given that each observed entry  Mi,jcan be viewed as  e⊤i Mej, the sampling basis relies heavily on the standard basis vectors. As a result, the above lemma is essentially imposing conditions on the incoherence between X and the sampling basis.

5.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 [16, 1922, 61, 63, 64], in the hope of improving computational guarantees. Specifically, one can regularize the loss function, by adding additional regularization term  G(·)to the objective function and designing the GD update rule w.r.t. the regularized problem


with  λ > 0the regularization parameter. For example:

Matrix completion (33): the following regularized loss has been proposed [16, 19, 65]


for some scalars  α1, · · · , α4 > 0. There are numerous choices of  G0; the one suggested by [19] is  G0(z) = max{z − 1, 0}2. With suitable learning rates and proper initialization, GD w.r.t. the regularized loss provably yields  ε-accuracy in  O(poly(n) log 1ε)iterations, provided that the sample size n2p ≳ µ2κ6nr7 log n [19].

Blind deconvolution (36): [22, 61, 66] suggest the following regularized loss


with  G0(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  ε-accuracy in  O(poly(m) log 1ε)iterations, provided that the sample size  m ≳ (K + N) log2 m [22].

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.

5.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 [26, 60, 62]. 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.

To illustrate this fact, we display in Fig. 3 a typical GD trajectory. The incoherence region — which enjoys local strong convexity and smoothness — is often a polytope (see the shaded region in the right panel of Fig. 3). The implicit regularization suggests that with high probability, the entire GD trajectory is constrained within this polytope, thus exhibiting linear convergence. It is worth noting that this cannot be derived from generic GD theory like Lemma 1. For instance, Lemma 1 implies that starting with a good initialization, the next iterate experiences  ℓ2error contraction, but it falls short of enforcing the incoherence condition and hence does not preclude the iterates from leaving the polytope.

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


where  x(⌧) := ⌧x + (1 � ⌧)xopt. The second line arises since  f(·)is  ↵-strongly convex, and the last line uses the smoothness condition.

where  x(⌧) := ⌧x + (1 � ⌧)xopt. The second line arises since  f(·)is  ↵-strongly convex, and the last line uses the smoothness condition.

where  x(⌧) := ⌧x + (1 � ⌧)xopt. The second line arises since  f(·)is  ↵-strongly convex, and the last line uses the smoothness condition.


where  x(⌧) := ⌧x + (1 � ⌧)xopt. The second line arises since  f(·)is  ↵-strongly convex, and the last line uses the smoothness condition.

where  x(⌧) := ⌧x + (1 � ⌧)xopt. The second line arises since  f(·)is  ↵-strongly convex, and the last line uses the smoothness condition.


Figure 3: The GD iterates and the locally strongly convex and smooth region (the shaded region). (Left) When this region is an  ℓ2ball, then standard GD theory implies  ℓ2convergence. (Right) When this region is a polytope, the implicit regularization phenomenon implies that the GD iterates still stay within this nice region.

Theorem 6 (GD for phase retrieval (improved bound) [26]). Under the assumptions of Lemma 5, the GD iterates with proper initialization (see, e.g., spectral initialization in Section 8.2) and  ηt ≡ 1/(c3∥x⋆∥22 log n)obey


for all  t ≥ 0with probability  1 − O(n−10). Here, c3, c4 > 0are some constants, and we assume  ∥x0 − x⋆∥2 ≤∥x0 + x⋆∥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 1ε)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.

Theorem 7 (GD for quadratic sensing [60]). Consider the problem (29). Suppose the sample size satisfies m ≥ c0nr4κ3 log nfor some large constant  c0, then with probability  1 − O(mn−10), the GD iterates with proper initialization (see, e.g., spectral initialization in Section 8.2) and  ηt = η ≡ 1/(c1(rκ + log n)2σ1(M⋆))obey


for all  t ≥ 0. Here, c1 > 0is some absolute constant.

This theorem demonstrates that vanilla GD converges within  O�max{r, log n}2 log 1ε�iterations for quadratic sensing of a rank-r matrix. This significantly improves upon the computational bounds in [59] which do not consider the incoherence region.

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


with  PΩdefined in (30). The theory for this update rule is:

Theorem 8 (GD for matrix completion [26, 67]). Consider the problem (32). Suppose that the sample size satisfies  n2p ≥ c0µ2r2n log nfor some large constant  c0 > 0, and that the condition number  κ of M⋆ = X⋆X⊤⋆is a fixed constant. With probability at least  1 − O�n−3�, the GD iterates (71) with proper initialization (see, e.g., spectral initialization in Section 8.2) satisfy


for all  t ≥ 0, with  HXtdefined in (50). Here,  c1 > 0is some constant, and  ηt ≡ c2/(κσ1(M⋆))for some constant  c2 > 0.

This theorem demonstrates that vanilla GD converges within  O�log 1ε�iterations. The key enabler of such a convergence rate is the property (72b), which basically implies that the GD iterates stay incoherent with the standard basis vectors. A byproduct is that: GD converges not only in Euclidean norm, it also converges in other more refined error metrics, e.g. the one measured by the  ℓ2/ℓ∞ norm.

The last example is blind deconvolution. To measure the discrepancy between any  z :=�hx�and z⋆ :=�h⋆x⋆�,we define


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


which enjoys the following theoretical support.

Theorem 9 (WF for blind deconvolution [26]). Consider the problem (36). Suppose the sample size  m ≥c0µ2 max{K, N}poly log mfor some large constant  c0 > 0. Then there is some constant  c1 > 0such 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  ηt ≡ c1 satisfy


For conciseness, we only state that the estimation error converges in  O�log 1ε�iterations. The incoherence conditions also provably hold across all iterations; see [26] for details. Similar results have been derived for the blind demixing case as well [62].

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 [6778]. The interested readers can find a general recipe using this argument in [26].

5.4 Notes

Provably valid two-stage nonconvex algorithms for matrix factorization were pioneered by Keshavan et al. [16, 65], 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 [16, 65] 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. [18], 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) [19], matrix sensing [23, 54], and blind deconvolution under subspace prior [22]. 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 [26] 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 [75].

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.

6.1 Projected gradient descent

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 C can be either convex or nonconvex. For many important sets 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.

6.1.1 Projection for computational benefits

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 [21, 63, 64], where in the positive semidefinite case (i.e.  M⋆ = X⋆X⊤⋆ ), one runs projected GD w.r.t. the loss function  f(·) in (32):


where  ηtis the step size and  PCdenote the Euclidean projection onto the set of incoherent matrices:


with  X0being 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, . . . , 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.

Theorem 10 (Projected GD for matrix completion [21, 64]). Suppose that the sample size satisfies  n2p ≥c0µ2r2n log nfor some large constant  c0 > 0, and that the condition number  κof  M⋆ = X⋆X⊤⋆is a fixed constant. With probability at least  1 − O�n−1�, the projected GD iterates (77) satisfy


for all  t ≥ 0, provided that  dist2(X0, X⋆) ≤ c3σr(M⋆)and  ηt ≡ η := c2/(µrσ1(M⋆))for some constant c1, c2, c3 > 0.


Remark 4. The results can be extended to the more general asymmetric case by applying similar modifications mentioned in Remark 3; see [63, 64].

6.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  x⋆ ∈ C. Therefore, it is natural to apply projection to enforce the desired structural constraints. One such example is sparse phase retrieval [79, 80], where it is known  a priori that x⋆in (26) is k-sparse, where  k ≪ n. If we have prior knowledge about  ∥x⋆∥1, then we can pick the constraint set C as follows to promote sparsity


as a sparse signal often (although not always) has low  ℓ1norm. With this convex constraint set in place, applying projected GD w.r.t. the loss function (27) can be efficiently implemented [81]. The theoretical guarantee of projected GD for sparse phase retrieval is given below.

Theorem 11 (Projected GD for sparse phase retrieval [79]). Consider the sparse phase retrieval problem where x⋆is k-sparse. Suppose that  m ≥ c1k log nfor some large constant  c1 > 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  ηt ≡ 1/(c2n∥x⋆∥22) and ∥x0 − x⋆∥2 ≤ ∥x⋆∥2/8.

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


This leads to a hard-thresholding operation, namely,  PC(x)becomes the best k-term approximation of x (obtained by keeping the k largest entries (in magnitude) of x and setting the rest to 0). The readers are referred to [82] for details. See also [80] 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  x0when  m ≲ 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  k2 log nunder i.i.d. Gaussian designs, unless k is sufficiently large or other structural information is available [48, 80, 8385]. All in all, the computational bottleneck for sparse phase retrieval lies in the initialization stage.

6.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 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 [20, 86, 87], 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 [24, 64, 88].

6.2.1 Truncation for computational and statistical benefits

We use phase retrieval to illustrate this benefit. All results discussed so far require a sample size that exceeds m ≳ n log n. When it comes to the sample-limited regime where  m ≍ 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 ≍ n, the gradient — which depends on 4th moments of  {ai}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  {∇fi(xt) :=�(a⊤i xt)2 − yi�aia⊤i xt}whose size deviate too much from the typical size. Specifically, the truncation rule proposed in [20] is:9


for some trimming criteria defined as


where  αlb, αub, αhare predetermined thresholds. This trimming rule — called Truncated Wirtinger flow — effectively removes the “heavy tails”, thus leading to much better concentration and hence enhanced performance.

Theorem 12 (Truncated GD for phase retrieval [20]). Consider the problem (27). With probability  1−O(n−10),the iterates (84) obey


for some constant  0 < ρ < 1, provided that  m ≥ c1n, ∥x0 − x⋆∥2 ≤ c2∥x⋆∥2and  ηt ≡ c3/∥x⋆∥22for some constants  c1, c2, c3 > 0.

Remark 5. In this case, the truncated gradient is clearly not smooth, and hence we need to resort to the regularity condition (see Definition 3) with  g(x) = ∇ftr(x). Specifically, the proof of Theorem 12 consists of showing


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 ≍ nis orderwise optimal;

Optimal computational complexity: truncated WF yields  εaccuracy in  O�log 1ε�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 [20] for precise statements.

6.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  αof revealed entries are corrupted by outliers. At the tth iterate, one can first try to identify the support of the sparse matrix  S⋆by hard thresholding the residual, namely,


Here, c > 1 is some predetermined constant (e.g. c = 3), and the operator  Hl(·)is defined as


where  A(l)j,· (resp. A(l)·,k) denotes the lth largest entry (in magnitude) in the jth row (resp. column) of 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  St+1then becomes our estimate of the sparse outlier matrix  S⋆in the (t + 1)-th iteration. With this in place, we update the estimate for the low-rank factor by applying projected GD


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

Theorem 13 (Nonconvex robust PCA [64]). Assume that the condition number  κof  M⋆ = X⋆X⊤⋆is a fixed constant. Suppose that the sample size and the sparsity of the outlier satisfy  n2p ≥ c0µ2r2n log nand α ≤ c1/(µr)for some constants  c0, c1 > 0. With probability at least  1 − O�n−1�, the iterates satisfy


for all  t ≥ 0, provided that  dist2(X0, X⋆) ≤ c3σr(M⋆). Here,  0 < c2, c3 < 1are some constants, and ηt ≡ c4/(µrσ1(M⋆))for some constant  c4 > 0.

Remark 6. In the full data case, the convergence rate can be improved to  1−c2for some constant  0 < c2 < 1.

This theorem essentially says that: as long as the fraction of entries corrupted by outliers does not exceed O(1/µr), then the nonconvex algorithm described above provably recovers the true low-rank matrix in about O(µ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 [24, 88]. We illustrate this median-truncation approach through an example of robust phase retrieval [24], where we assume a subset of samples in (26) is corrupted arbitrarily, with their index set denoted by  S with |S| = αm. Mathematically, the measurement model in the presence of outliers is given by


The goal is to still recover  x⋆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 ∈ S. Hence, we instead include only a subset of the samples when forming the search direction, yielding a truncated GD update rule


Here,  Ttonly includes samples whose residual size  |rt,i|does not deviate much from the  median of {|rt,j|}1≤j≤m:


where  median(·)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.

Theorem 14 (Median-truncated GD for robust phase retrieval [24]). Consider the problem (89) with a fraction αof arbitrary outliers. There exist some constants  c0, c1, c2 > 0and  0 < ρ < 1such that if  m ≥ c0n log n, α ≤ c1, and  ∥x0 − x⋆∥2 ≤ c2∥x⋆∥2, then with probability at least  1 − O(n−1), the median-truncated GD iterates satisfy


6.3 Generalized gradient descent

In all the examples discussed so far, the loss function  f(·)has been a smooth function. When  f(·)is nonsmooth and non-differentiable, it is possible to continue to apply GD using the generalized gradient (e.g. subgradient) [89]. 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


Clearly,  famp(x)is nonsmooth, and its generalized gradient is given by, with a slight abuse of notation,


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


This amplitude-based loss function  famp(·)often has better curvature around the truth, compared to the intensity-based loss function  f(·) defined in (29); see [87, 90, 91] for detailed discussions. The theory is given as follows.

Theorem 15 (GD for amplitude-based phase retrieval [90]). Consider the problem (27). There exist some constants  c0, c1, c2 > 0 and 0 < ρ < 1such that if  m ≥ c0n and ηt ≡ c2, then with high probability,


as long as  ∥x0 − x⋆∥2 ≤ ∥x⋆∥2/10.

In comparison to Theorem 12, the generalized GD w.r.t.  famp(·)achieves both order-optimal sample and computational complexities. Notably, a very similar theory was obtained in [87] for a truncated version of the generalized GD (called Truncated Amplitude Flow therein), where the algorithm also employs the gradient update w.r.t.  famp(·)but discards high-leverage data in a way similar to truncated GD discussed in Section 6.2. However, in contrast to the intensity-based loss  f(·) defined in (29), the truncation step is not crucial and can be safely removed when dealing with the amplitude-based  famp(·). A main reason is that for any fixed x, famp(·)involves only the first and second moments of the (sub)-Gaussian random variables  {a⊤i x}. As such,it exhibits much sharper measure concentration — and hence much better controlled gradient components — compared to the heavy-tailed  f(·), which involves fourth moments of  {a⊤i x}. This observation in turn implies the importance of designing loss functions for nonconvex statistical estimation.

6.4 Projected power method for constrained PCA

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


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

Phase synchronization [92, 93]. Suppose we wish to recover n unknown phases  φ1, · · · , φn ∈ [0, 2π]given their pairwise relative phases. Alternatively, by setting  (x⋆)i = eφi, this problem reduces to estimating  x⋆ = [(x⋆)i]1≤i≤nfrom  x⋆xH⋆— a matrix that encodes all pairwise phase differences (x⋆)i(x⋆)Hj = eφi−φj. To account for the noisy nature of practical measurements, suppose that what we observe is  L = x⋆xH⋆ + σW, where W is a Hermitian matrix. Here,  {Wi,j}i≤jare i.i.d. standard complex Gaussians. The quantity  σ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 [10, 94]. Imagine we want to estimate n discrete variables  {(x⋆)i}1≤i≤n, where each variable can take m possible values, namely,  (x⋆)i ∈ {1, · · · , m}. Suppose that estimation needs to be performed based on pairwise difference samples  yi,j = xi − xj + zi,j mod m, where the  zi,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  xiinto a m-dimensional vector  xi ∈ {e1, · · · , em}. We then introduce a matrix L that properly encodes all log-likelihood information. After simple manipulation (see [10] for details), maximum likelihood estimation can be cast as follows


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

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) [10, 96], 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 C (e.g. the ones in phase synchronization and joint alignment), (97) is equivalent to projected GD when the step size  ηt → ∞.

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  σ ≲�n/ log n. This is information theoretically optimal up to some log factor [71, 92].

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 [10].

6.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 [98, 99].

To fix ideas, take matrix completion as an example. When factorizing  M⋆ = L⋆R⊤⋆, we might assume L⋆ ∈ G(n1, r), where  G(n1, r)denotes the Grassmann manifold which parametrizes all r-dimensional linear subspaces of the  n1-dimensional space10. In words, we are searching for a r-dimensional subspace L but ignores the global rotation. It is also assumed that  L⊤⋆ L⋆ = Irto remove the global scaling ambiguity (otherwise  (cL, c−1R)is always equivalent to  (L, R) for any c ̸= 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(·)over the Grassmann manifold by moving along the geodesics; here, a geodesic is the shortest path between two points on a manifold. See [98] 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(·)w.r.t. the  tth iterate Lt in the Euclidean space[100]:


where ˆRt = argminR��PΩ(M⋆ − LtR⊤)��2Fis the least-squares solution. The gradient on the Grassmann manifold, denoted by  ∇GF(·), is then given by


Let  −∇GF(Lt) = ˜Ut ˜Σt ˜V ⊤tbe its compact SVD, then the geodesic on the Grassmann manifold along the direction  −∇GF(Lt)is given by


We can then update the iterates as


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


with  σ := ∥∇GF(Lt)∥2. As can be verified,  Lt+1automatically stays on the unit sphere obeying  L⊤t+1Lt+1 = 1.One of the earliest provable nonconvex methods for matrix completion — the OptSpace algorithm by Keshavan et al. [16, 65] — performs gradient descent on the Grassmann manifold, tailored to the loss function:


where  L ∈ G(n1, r) and R ∈ G(n2, r), with some additional regularization terms to promote incoherence (see Section 5.2.2). It is shown by [16, 65] that GD on the Grassman manifold converges to the truth with high probability if  n2p ≳ µ2κ6r2n log n, provided that a proper initialization is given.

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

6.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) [108110], 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  Ωt ∈ {1, . . . , m}is a subset of cardinality k selected uniformly at random. Here,  k ≥ 1is 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 [90, 111115] for the application of SGD in phase retrieval (which has an interesting connection with the Kaczmarz method), and [116] for its application in matrix factorization.

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)


7.1 Alternating minimization

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


where  Rtand  Ltare updated sequentially. Here,  L0is an appropriate initialization. For many problems discussed here, both (107a) and (107b) are convex problems and can be solved efficiently.

7.1.1 Matrix sensing

Consider the loss function (21). In each iteration, AltMin proceeds as follows [17]: for t = 1, 2 . . . ,


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

Example 2. Consider the case where A is identity (i.e. A(M) = 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  L⊤⋆ L0and  L⊤0 L0are full-rank, the column space of  R⋆matches perfectly with that of  R1. Armed with this fact, the subsequent least squares problem (i.e. the update for  L1) is exact, in the sense that L1R⊤1 = M⋆ = L⋆R⊤⋆ .

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

Theorem 16 (AltMin for matrix sensing [17]). Consider the problem (21) and suppose the operator (23) satisfies 2r-RIP with RIP constant  δ2r ≤ 1/(100κ2r). If we initialize  L0 by the rleading left singular vectors of  A∗(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  κ; that is, it obtains  ε-accuracy within  O(log(1/ε))iterations, compared to O(κ log(1/ε))iterations for GD. In addition, the requirement on the RIP constant depends quadratically on  κ, leading to a sub-optimal sample complexity. To address this issue, Jain et al. [17] further developed a stage-wise AltMin algorithm, which only requires  δ2r = O(1/r2). Intuitively, if there is a singular value that is much larger than the remaining ones, then one can treat  M⋆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 [17] for details.

7.1.2 Phase retrieval

Consider the phase retrieval problem. It is helpful to think of the amplitude measurements as bilinear measurements of the signs  b = {bi ∈ {±1}}1≤i≤mand the signal  x⋆, namely,


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


where we abuse the notation by letting


Therefore, by applying AltMin to the loss function f(b, x), we obtain the following update rule [118, 119]: for each t = 1, 2, . . .


where  x0is an appropriate initial estimate,  A†is the pseudo-inverse of  A := [a1, · · · , am]⊤, and √y :=[√yi]1≤i≤m. The step (110b)can again be efficiently solved using the conjugate gradient method [117]. This is exactly the Error Reduction (ER) algorithm proposed by Gerchberg and Saxton [45, 120] in the 1970s. Given a reasonably good initialization, this algorithm converges linearly under the Gaussian design.

Theorem 17 (AltMin (ER) for phase retrieval [119]). Consider the problem (27). There exist some constants c0, · · · , c4 > 0and  0 < ρ < 1such that if  m ≥ c0n, then with probability at least  1 − c2 exp(−c3m), the estimates of AltMin (ER) satisfy


as long as  ∥x0 − x⋆∥2 ≤ c4∥x⋆∥2.

Remark 7. AltMin for phase retrieval was first analyzed by [118] but for a sample-splitting variant; that is, each iteration employs fresh samples, which facilitates analysis but is not the version used in practice. The theoretical guarantee for the original sample-reuse version was derived by [119].

In view of Theorem 17, alternating minimization, if carefully initialized, achieves optimal sample and computational complexities (up to some logarithmic factor) all at once. This in turn explains its appealing performance in practice.

7.1.3 Matrix completion

Consider the matrix completion problem in (33). Starting with a proper initialization  (L0, R0), AltMin proceeds as follows: for t = 1, 2, . . .


where  PΩis defined in (30). Despite its popularity in practice [121], a clean analysis of the above update rule is still missing to date. Several modifications have been proposed and analyzed in the literature, primarily to bypass mathematical difficulty:

Sample splitting. Instead of reusing the same set of samples across all iterations, this approach draws a fresh set of samples at every iteration and performs AltMin on the new samples [17, 122125]:


where  Ωtdenotes the sampling set used in the tth iteration, which is assumed to be statistically independent across iterations. It is proven in [17] that under an appropriate initialization, the output satisfies��M⋆ − LtR⊤t��F ≤ εafter  t ≳ log(∥M⋆∥F/ε)iterations, provided that the sample complexity exceeds  n2p ≳ µ4κ6nr7 log n log(r∥M⋆∥F/ε). Such a sample-splitting operation ensures statistical independence across iterations, which helps to control the incoherence of the iterates. However, this necessarily results in undesirable dependency between the sample complexity and the target accuracy; for example, an infinite number of samples is needed if the goal is to achieve exact recovery.

Regularization. Another strategy is to apply AltMin to the regularized loss function in (67) [19]:


In [19], it is shown that the AltMin without resampling converges to  M⋆, with the proviso that the sample complexity exceeds  n2p ≳ µ2κ6nr7 log n. Note that the subproblems (113a) and (113b) do not have closed-form solutions. For properly chosen regularization functions, they might be solved using convex optimization algorithms.

It remains an open problem to establish theoretical guarantees for the original form of AltMin (112). Meanwhile, the existing sample complexity guarantees are quite sub-optimal in terms of the dependency on r and  κ, and should not be taken as an indicator of the actual performance of AltMin.

7.2 Singular value projection

Another popular approach to solve (106) is singular value projection (SVP) [76, 126, 127]. In contrast to the algorithms discussed so far, SVP performs gradient descent in the full matrix space and then applies a partial singular value decomposition (SVD) to retain the low-rank structure. Specifically, it adopts the update rule




Here,  ηtis the step size, and  Pr(Z)returns the best rank-r approximation of Z. Given that the iterates  Mtare always low-rank, one can store  Mtin a memory-efficient manner by storing its compact SVD.

The SVP algorithm is a popular approach for matrix sensing and matrix completion, where the partial SVD can be calculated using Krylov subspace methods (e.g. Lanczos algorithm) [117] or the randomized linear algebra algorithms [128]. The following theorem establishes performance guarantees for SVP.

Theorem 18 (SVP for matrix sensing [126]). Consider the problem (106) and suppose the operator (23) satisfies 2r-RIP for RIP constant  δ2r ≤ 1/3. If we initialize  M0 = 0and adopt a step size  ηt ≡ 1/(1 + δ2r),then SVP achieves


Theorem 19 (SVP for matrix completion [76]). Consider the problem (32) and set  M0 = 0. Suppose that n2p ≥ c0κ6µ4r6n log nfor some large constant  c0 > 0, and that the step size is set as  ηt ≡ 1. Then with probability exceeding  1 − O�n−10�, the SVP iterates achieve


as long as  t ≥ c1 log(∥M⋆∥/ε)for some constant  c1 > 0.

Theorem 18 and Theorem 19 indicate that SVP converges linearly as soon as the sample size is sufficiently large.

Further, the SVP operation is particularly helpful in enabling optimal uncertainty quantification and inference for noisy matrix completion (e.g. constructing a valid and short confidence interval for an entry of the unknown matrix). The interested readers are referred to [78] for details.

7.3 Further pointers to other algorithms

A few other nonconvex matrix factorization algorithms have been left out due to space, including but not limited to normalized iterative hard thresholding (NIHT) [129], atomic decomposition for minimum rank approximation (Admira) [130], composite optimization (e.g. prox-linear algorithm) [131134], approximate message passing [135137], block coordinate descent [19], coordinate descent [138], and conjugate gradient [139]. The readers are referred to these papers for detailed descriptions.

The theoretical performance guarantees presented in the last three sections rely heavily on proper initialization. One popular scheme that often generates a reasonably good initial estimate is called the spectral method. Informally, this strategy starts by arranging the data samples into a matrix Y of the form


where  Y⋆represents certain large-sample limit whose eigenspace / singular subspaces reveal the truth, and  ∆captures the fluctuation due to the finite-sample effect. One then attempts to estimate the truth by computing the eigenspace / singular subspace of Y , provided that the finite-sample fluctuation is well-controlled. This simple strategy has proven to be quite powerful and versatile in providing a “warm start” for many nonconvex matrix factorization algorithms.

8.1 Preliminaries: matrix perturbation theory

Understanding the performance of the spectral method requires some elementary toolkits regarding eigenspace / singular subspace perturbations, which we review in this subsection.

To begin with, let  Y⋆ ∈ Rn×n be a symmetric matrix, whose eigenvalues are real-valued. In many cases, we only have access to a perturbed version  Y of Y⋆ (cf. (115)), where the perturbation  ∆is a “small” symmetric matrix. How do the eigenvectors of Y change as a result of such a perturbation?

As it turns out, the eigenspace of Y is a stable estimate of the eigenspace of  Y⋆, with the proviso that the perturbation is sufficiently small in size. This was first established in the celebrated Davis-Kahan  sin ΘTheorem [140]. Specifically, let the eigenvalues of  Y⋆be partitioned into two groups


where  1 ≤ r < n. We assume that the eigen-gap between the two groups,  λr(Y⋆) − λr+1(Y⋆), is strictly positive. For example, if  Y⋆ ⪰ 0and has rank r, then all the eigenvalues in the second group are identically zero.

Suppose we wish to estimate the eigenspace associated with the first group. Denote by  U⋆ ∈ Rn×r(resp.  U ∈ Rn×r) an orthonormal matrix whose columns are the first r eigenvectors of  Y⋆ (resp. Y). In order to measure the distance between the two subspaces spanned by  U⋆ and U, we introduce the following metric that accounts for global orthonormal transformation


Remark 8. The bound we present in (117) is in fact a simplified, but slightly more user-friendly version of the original Davis-Kahan inequality. A more general result states that, if  λr(Y⋆) − λr+1(Y ) > 0, then


These results are referred to as the  sin Θtheorem because the distance metric  distp(U, U⋆)is identical to max1≤i≤r | sin θi|, where {θi}1≤i≤rare the so-called principal angles [141] between the two subspaces spanned by  U and U⋆, respectively.

Furthermore, to deal with asymmetric matrices, similar perturbation bounds can be obtained. Suppose that  Y⋆, Y , ∆ ∈ Rn1×n2 in (115). Let U⋆ (resp. U) and V⋆ (resp. V) consist respectively of the first r left and right singular vectors of  Y⋆(resp. Y ). Then we have the celebrated Wedin  sin ΘTheorem [142] concerning perturbed singular subspaces.


In addition, one might naturally wonder how the eigenvalues / singular values are affected by the perturbation. To this end, Weyl’s inequality provides a simple answer:


In summary, both eigenspace (resp. singular subspace) perturbation and eigenvalue (resp. singular value) perturbation rely heavily on the spectral norm of the perturbation  ∆.

8.2 Spectral methods

With the matrix perturbation theory in place, we are positioned to present spectral methods for various low-rank matrix factorization problems. As we shall see, these methods are all variations on a common recipe.

8.2.1 Matrix sensing

We start with the prototypical problem of matrix sensing (19) as described in Section 4.1. Let us construct a surrogate matrix as follows


where A is the linear operator defined in (23) and  A∗ denotes its adjoint. The generic version of the spectral method then proceeds by computing (i) two matrices  U ∈ Rn1×rand  V ∈ Rn2×rwhose columns consist of the top-r left and right singular vectors of Y , respectively, and (ii) a diagonal matrix  Σ ∈ Rr×r that contains the corresponding top-r singular values. In the hope that U, V and  Σare reasonably reliable estimates of U⋆, V⋆ and Σ⋆, respectively, we take


as estimates of the low-rank factors  L⋆ and R⋆ in (22). If M⋆ ∈ Rn×n is known to be positive semidefinite, then we can also let  U ∈ Rn×rbe a matrix consisting of the top-r leading eigenvectors, with  Σbeing a diagonal matrix containing all top-r eigenvalues.

Why would this be a good strategy? In view of Section 8.1, the three matrices U, V and  Σbecome reliable estimates if  ∥Y − M⋆∥can be well-controlled. A simple way to control  ∥Y − M⋆∥arises when A satisfies the RIP in Definition 6.

Lemma 8. Suppose that  M⋆is a rank-r matrix, and assume that A satisfies 2r-RIP with RIP constant


Proof of Lemma 8. Let x, y be two arbitrary vectors, then


where the last inequality is due to Lemma 3. Using a variational characterization of  ∥ · ∥, we have


from which (123) follows.

In what follows, we first illustrate how to control the estimation error in the rank-1 case where  M⋆ =λ⋆u⋆u⊤⋆ ⪰ 0. In this case, the leading eigenvector u of Y obeys


where (i) comes from Theorem 20, and (ii) holds if  ∥Y − M⋆∥ ≤ σ1(M⋆)/2(which is guaranteed if  δ2 ≤ 1/2according to Lemma 8). Similarly, we can invoke Weyl’s inequality and Lemma 8 to control the gap between the leading eigenvalue  λ of Y and λ⋆:


where the last inequality makes use of (124), (125), and (??). This characterizes the difference between our estimate√λuand the true low-rank factor √λ⋆u⋆.

Moving beyond this simple rank-1 case, a more general (and often tighter) bound can be obtained by using a refined argument from [23, Lemma 5.14]. We present the theory below. The proof can be found in Appendix D.

Theorem 22 (Spectral method for matrix sensing [23]). Fix  ζ > 0. Suppose A satisfies 2r-RIP with RIP constant  δ2r < c0√ζ/(√rκ)for some sufficiently small constant  c0 > 0. Then the spectral estimate (122) obeys


Remark 9. It is worth pointing out that in view of Fact 1, the vanilla spectral method needs  O(nr2) samplesto land in the local basin of attraction (in which linear convergence of GD is guaranteed according to Theorem 4).

As discussed in Section 5, the RIP does not hold for the sensing matrices used in many problems. Nevertheless, one may still be able to show that the leading singular subspace of the surrogate matrix Y contains useful information about the truth  M⋆. In the sequel, we will go over several examples to demonstrate this point.

8.2.2 Phase retrieval

Recall that the phase retrieval problem in Section 4.2 can be viewed as a matrix sensing problem, where we seek to recover a rank-1 matrix  M⋆ = x⋆x⊤⋆ with sensing matrices  Ai = aia⊤i . To obtain an initial guess  x0that is close to the truth  x⋆, we follow the recipe described in (121) by estimating the leading eigenvector u and leading eigenvalue  λof a surrogate matrix


The initial guess is then formed as11


Unfortunately, the RIP does not hold for the sensing operator in phase retrieval, which precludes us from invoking Theorem 22. There is, however, a simple and intuitive explanation regarding why  x0is a reasonably good estimate of  x⋆. Under the Gaussian design, the surrogate matrix Y in (126) can be viewed as the sample average of m i.i.d. random rank-one matrices�yiaia⊤i�1≤i≤m. When the number of samples m is large, this sample average should be “close” to its expectation, which is,


The best rank-1 approximation of E[Y ] is precisely  3x⋆x⊤⋆ . Now that Yis an approximated version of E[Y ], we expect  x0 in (127) to carry useful information about  x⋆.

The above intuitive arguments can be made precise. Applying standard matrix concentration inequalities [57] to the surrogate matrix in (126) and invoking the Davis-Kahan  sin Θtheorem, one arrives at the following estimates:

Theorem 23 (Spectral method for phase retrieval [18]). Consider phase retrieval in Section 4.2, where  x⋆ ∈Rn is any given vector. Fix any  ζ > 0, and suppose  m ≥ c0n log nfor some sufficiently large constant  c0 > 0.Then the spectral estimate (127) obeys


with probability at least  1 − O(n−2).

8.2.3 Quadratic sensing

An argument similar to phase retrieval can be applied to quadratic sensing in (29), recognizing that the expectation of the surrogate matrix Y in (126) now becomes


The spectral method then proceeds by computing U (which consists of the top-r eigenvectors of Y ), and a diagonal matrix  Σwhose ith diagonal value is given as  (λi(Y ) − σ)/2, where  σ = 1m�mi=1 yiserves as an estimate of  ∥X⋆∥2F. In words, the diagonal entries of  Σcan be approximately viewed as the top-r eigenvalues of 12�Y − ∥X⋆∥2F In�. The initial guess is then set as


for estimating the low-rank factor  X⋆. The theory is as follows.

Theorem 24 (Spectral method for quadratic sensing [18]). Consider quadratic sensing in Section 4.2, where X⋆ ∈ Rn×r. Fix any  ζ > 0, and suppose  m ≥ c0nr4 log nfor some sufficiently large constant  c0 > 0. Then the spectral estimate (130) obeys


with probability at least  1 − O(n−2).

8.2.4 Blind deconvolution

The blind deconvolution problem introduced in Section 4.4 has a similar mathematical structure to that of phase retrieval. Recall the sensing model in (35). Instead of reconstructing a symmetric rank-1 matrix, we now aim to recover an asymmetric rank-1 matrix  h⋆xH⋆with sensing matrices  Ai = biaHi (1 ≤ i ≤ m). Following (121), we form a surrogate matrix


Let  u, v, and σdenote the leading left singular vector, right singular vector, and singular value, respectively. The initial guess is then formed as


This estimate provably reveals sufficient information about the truth, provided that the sample size m is sufficiently large.

Theorem 25 (Spectral method for blind deconvolution [22, 26]). Consider blind deconvolution in Section 4.4. Suppose  h⋆satisfies the incoherence condition in Definition 8 with parameter  µ, and assume  ∥h⋆∥2 = ∥x⋆∥2.For any  ζ > 0, if m ≥ c0ζ−2µ2K log2 mfor some sufficiently large constant  c0 > 0, then the spectral estimate (131) obeys


with probability at least  1 − O(m−10).

8.2.5 Matrix completion

Turning to matrix completion as introduced in Section 4.3, which is another instance of matrix sensing with sensing matrices taking the form of


Then the measurements obey  ⟨Ai,j, M⋆⟩ = 1√p(M⋆)i,j. Following the aforementioned procedure, we can form a surrogate matrix as


Notably, the scaling factor in (132) is chosen to ensure that  E[Y ] = M⋆. We then construct the initial guess for the low-rank factors  L0 and R0in the same manner as (122), using Y in (133).

As�Ai,j1{(i,j)∈Ω}�is a collection of independent random matrices, we can use the matrix Bernstein inequality [57] to get a high-probability upper bound on the deviation  ∥Y − E[Y ]∥. This in turn allows us to apply the matrix perturbation bounds to control the accuracy of the spectral method.

Theorem 26 (Spectral method for matrix completion [19, 21, 26]). Consider matrix completion in Section 4.3. Fix  ζ > 0, and suppose the condition number  κ of M ⋆ is a fixed constant. There exist a constant  c0 > 0 suchthat if  np > c0µ2r2 log n, then with probability at least  1 − O(n−10), the spectral estimate (122) obeys


8.3 Variants of spectral methods

We illustrate modifications to the spectral method, which are often found necessary to further enhance sample efficiency, increase robustness to outliers, and incorporate signal priors.

8.3.1 Truncated spectral method for sample efficiency

The generic recipe for spectral methods described above works well when one has sufficient samples compared to the underlying signal dimension. It might not be effective though if the sample complexity is on the order of the information-theoretic limit.

In what follows, we use phase retrieval to demonstrate the underlying issues and how to address them12.Recall that Theorem 23 requires a sample complexity  m ≳ n log n, which is a logarithmic factor larger than the signal dimension n. What happens if we only have access to  m ≍ nsamples, which is the information-theoretic limit (order-wise) for phase retrieval? In this more challenging regime, it turns out that we have to modify the standard recipe by applying appropriate preprocessing before forming the surrogate matrix Y in (121).

We start by explaining why the surrogate matrix (126) for phase retrieval must be suboptimal in terms of sample complexity. For all  1 ≤ j ≤ m,


Under the Gaussian design,�yi/∥x⋆∥22�is a collection of i.i.d.  χ2 random variables with 1 degree of freedom. It follows from well-known estimates in extreme value theory [143] that


Recall that  E[Y ] = 2x⋆x⊤⋆ + ∥x⋆∥22In has a bounded spectral norm, then (135) implies that, to keep the deviation between Y and E[Y ] well-controlled, we must at least have


This condition, however, cannot be satisfied when we have linear sample complexity  m ≍ n. This explains why we need a sample complexity  m ≳ n log nin Theorem 23.

The above analysis also suggests an easy fix: since the main culprit lies in the fact that  maxi yiis unbounded (as  m → ∞), we can apply a preprocessing function  T (·)to  yito keep the quantity bounded. Indeed, this is the key idea behind the truncated spectral method proposed by Chen and Candès [20], in which the surrogate matrix is modified as




for some predetermined truncation threshold  γ. The initial point  x0is then formed by scaling the leading eigenvector of  YTto have roughly the same norm of  x⋆, which can be estimated by  σ = 1m�mi=1 yi. This isessentially performing a trimming operation, removing any entry of  yithat bears too much influence on the leading eigenvector. The trimming step turns out to be very effective, allowing one to achieve order-wise optimal sample complexity.

Theorem 27 (Truncated spectral method for phase retrieval [20]). Consider phase retrieval in Section 4.2. Fix any  ζ > 0, and suppose  m ≥ c0nfor some sufficiently large constant  c0 > 0. Then the truncated spectral estimate obeys


with probability at least  1 − O(n−2).

Subsequently, several different designs of the preprocessing function have been proposed in the literature. One example was given in [87], where


where  γ is the (cm)-largest value (e.g.  c = 1/6) in {yj}1≤j≤m. In words, this method only employs a subset of design vectors that are better aligned with the truth  x⋆. By properly tuning the parameter c, this truncation scheme performs competitively as the scheme in (137).

8.3.2 Truncated spectral method for removing sparse outliers

When the samples are susceptible to adversarial entries, e.g. in the robust phase retrieval problem (89), the spectral method might not work properly even with the presence of a single outlier whose magnitude can be arbitrarily large to perturb the leading eigenvector of Y . To mitigate this issue, a median-truncation scheme was proposed in [24, 88], where


for some predetermined constant  γ > 0. By including only a subset of samples whose values are not excessively large compared with the sample median of the samples, the preprocessing function in (139) makes the spectral method more robust against sparse and large outliers.

Theorem 28 (Median-truncated spectral method for robust phase retrieval [24]). Consider the robust phase retrieval problem in (89), and fix any  ζ > 0. There exist some constants  c0, c1 > 0such that if  m ≥ c0n andα ≤ c1, then the median-truncated spectral estimate obeys


with probability at least  1 − O(n−2).

The idea of applying truncation to form a spectral estimate is also used in the robust PCA problem (see Section 4.5). Since the observations are also potentially corrupted by large but sparse outliers, it is useful to first clean up the observations  (Γ⋆)i,j = (M⋆)i,j + (S⋆)i,jbefore constructing the surrogate matrix as in (133). Indeed, this is the strategy proposed in [64]. We start by forming an estimate of the sparse outliers via the hard-thresholding operation  Hl(·)defined in (86), as


where c > 0 is some predetermined constant (e.g.  c = 3) and PΩis defined in (30). Armed with this estimate, we form the surrogate matrix as


One can then apply the spectral method to Y in (141) (similar to the matrix completion case). This approach enjoys the following performance guarantee.

Theorem 29 (Spectral method for robust PCA [64]). Suppose that the condition number  κof  M⋆ = L⋆R⊤⋆is a fixed constant. Fix  ζ > 0. If the sample size and the sparsity fraction satisfy  n2p ≥ c0µr2n log nand α ≤ c1/(µr3/2)for some large constant  c0, c1 > 0, then with probability at least  1 − O�n−1�,


8.3.3 Spectral method for sparse phase retrieval

Last but not least, we briefly discuss how the spectral method can be modified to incorporate structural priors. As before, we use the example of sparse phase retrieval to illustrate the strategy, where we assume  x⋆is k-sparse (see Section 6.1.2). A simple idea is to first identify the support of  x⋆, and then try to estimate the nonzero values by applying the spectral method over the submatrix of Y on the estimated support. Towards this end, we recall that  E[Y ] = 2x⋆x⊤⋆ + ∥x⋆∥22In, and therefore the larger ones of its diagonal entries are more likely to be included in the support. In light of this, a simple thresholding strategy adopted in [80] is to compare  Yi,iagainst some preset threshold  γ:


The nonzero part of  x⋆is then found by applying the spectral method outlined in Section 8.2.2 to  Y ˆS =1

�mi=1 yiai, ˆSa⊤i, ˆS, where  ai, ˆSis the subvector of  aicoming from the support  ˆS. In short, this strategy provably leads to a reasonably good initial estimate, as long as  m ≳ k2 log n; the complete theory can be found in [80]. See also [82] for a more involved approach, which provides better empirical performance.

8.4 Precise asymptotic characterization and phase transitions for phase retrieval

The Davis-Kahan and Wedin  sin ΘTheorems are broadly applicable and convenient to use, but they usually fall short of providing the tightest estimates. For many problems, if one examines the underlying statistical models carefully, it is often possible to obtain much more precise performance guarantees for spectral methods.

In [144], Lu and Li provided an asymptotically exact characterization of spectral initialization in the context of generalized linear regression, which subsumes phase retrieval as a special case. One way to quantify the quality of this eigenvector is via the squared cosine similarity


which measures the (squared) correlation between the truth  x⋆and the initialization  x0. The result is this:

Theorem 30 (Precise asymptotic characterization [144]). Consider phase retrieval in Section 4.2, and let m/n = αfor some constant  α > 0. Under mild technical conditions on the preprocessing function T , the leading eigenvector  x0/∥x0∥2 of YT in (136) obeys


as  n → ∞. Here,  αc > 0is a fixed constant and  ρ∗(·)is a fixed function that is positive when  α > αc. Furthermore,  λ1(YT ) − λ2(YT )converges to a positive constant iff  α > αc.

Remark 10. The above characterization was first obtained in [144], under the assumption that  T (·)is nonnegative. Later, this technical restriction was removed in [145]. Analytical formulas of  αc and ρ∗(α) areavailable for any given  T (·), which can be found in [144, 145].

The asymptotic prediction given in Theorem 30 reveals a phase transition phenomenon: there is a critical sampling ratio  αcthat marks the transition between two very contrasting regimes.

An uncorrelated phase takes place when the sampling ratio  α < αc. Within this phase,  ρ(x⋆, x0) → 0,meaning that the spectral estimate is uncorrelated with the target.

A correlated phase takes place when  α > αc. Within this phase, the spectral estimate is strictly better than a random guess. Moreover, there is a nonzero gap between the 1st and 2nd largest eigenvalues of YT, which in turn implies that  x0can be efficiently computed by the power method.


Figure 4: Comparisons of 3 designs of the preprocessing function for estimating a  64 × 64 cameramanimage from phaseless measurements under Poisson noise. (Red curve) the subset scheme in (138), where the parameter  γhas been optimized for each fixed  α; (Green curve) the function  T ∗α (·)defined in (144); (Blue curve) the uniformly optimal design  T ∗(·) in (145).

The phase transition boundary  αcis determined by the preprocessing function  T (·). A natural question arises as to which preprocessing function optimizes the phase transition point. At first glance, this seems to be a challenging task, as it is an infinite-dimensional functional optimization problem. Encouragingly, this can be analytically determined using the asymptotic characterizations stated in Theorem 30.

Theorem 31 (Optimal preprocessing [145]). Consider phase retrieval in the real-valued setting. The phase transition point satisfies  αc > α∗ = 1/2for any  T (·). Further,  α∗can be approached by the following preprocessing function


The value  α∗ is called the  weak recovery threshold. When α < α∗, no algorithm can generate an estimate that is asymptotically positively correlated with  x⋆. The function  T ∗α (·)is optimal in the sense that it approaches this weak recovery threshold.

Another way to formulate optimality is via the squared cosine similarity in (142). For any fixed sampling ratio  α > 0, we seek a preprocessing function that maximizes the squared cosine similarity, namely,


The following theorem [146] shows that the fixed function  T ∗(·)is in fact uniformly optimal for all sampling ratio  α. Therefore, instead of using (144) which takes different forms depending on  α, one should directly use  T ∗(·).

Theorem 32 (Uniformly optimal preprocessing). Under the same settings of Theorem 31, we have  T optα (·) =T ∗(·), where T ∗(·)is defined in (145).

Finally, to demonstrate the improvements brought by the optimal preprocessing function, we show in Fig. 4 the results of applying the spectral methods to estimate a  64 × 64 cameramanimage from phaseless measurements under Poisson noise. It is evident that the optimal design significantly improves the performance of the method.

8.5 Notes

The idea of spectral methods can be traced back to the early work of Li [147], under the name of Principal Hessian Directions for general multi-index models. Similar spectral techniques were also proposed in [16, 17, 148], for initializing algorithms for low-rank matrix completion. Regarding phase retrieval, Netrapalli et al. [118] used this method to address the problem of phase retrieval, the theoretical guarantee of which was tightened in [18]. Similar guarantees were also provided for the randomly coded diffraction pattern model in [18]. The first order-wise optimal spectral method was proposed by Chen and Candès [20], based on the truncation idea. This method has multiple variants [87, 88], and has been shown to be robust against noise. The precise asymptotic characterization of the spectral method was first obtained in [144]. Based on this characterization, [145] determined the optimal weak reconstruction threshold for spectral methods.

Finally, the spectral method has been applied to many other problems beyond the ones discussed here, including but not limited to community detection [149, 150], phase synchronization [9], joint alignment [10], ranking from pairwise comparisons [73, 151, 152], tensor estimation [153156]. We have to omit these due to the space limit.

A separate line of work aims to study the global geometry of a loss function  f(·)over the entire parameter space, often under appropriate statistical models of the data. As alluded by the warm-up example in Section 3.2, such studies characterize the critical points and geometric curvatures of the loss surface, and highlight the (non-)existence of spurious local minima. The results of the geometric landscape analysis can then be used to understand the effectiveness of a particular optimization algorithm of choice.

9.1 Global landscape analysis

In general, global minimization requires one to avoid two types of undesired critical points: (1) local minima that are not global minima; (2) saddle points. If all critical points of a function  f(·)are either global minima or strict saddle points, we say that  f(·)has benign landscape. Here, we single out strict saddles from all possible saddle points, since they are easier to escape due to the existence of descent directions.13

Loosely speaking, nonconvexity arises in these problems partly due to “symmetry”, where the global solutions are identifiable only up to certain global transformations. This necessarily leads to multiple indistinguishable local minima that are globally optimal. Further, saddle points arise naturally when interpolating the loss surface between two separated local minima. Nonetheless, in spite of nonconvexity, a large family of problems exhibit benign landscape. This subsection gives a few such examples.

9.1.1 Two-layer linear neural network

A straightforward instance that has already been discussed is the warm-up example in Section 3. It can be slightly generalized as follows.

Example 3 (Two-layer linear neural network [36]). Given arbitrary data {xi, yi}mi=1, xi, yi ∈ Rn, we wish to fit a two-layer linear network (see Fig. 5) using the quadratic loss:


In particular, when X = Y , Example 3 reduces to rank-r matrix factorization [or principal component analysis (PCA)], an immediate extension of the rank-1 warm-up example. When  X ̸= Y, Example 3 is precisely the canonical correlation analysis (CCA) problem. This explains why both PCA and CCA, though highly nonconvex, admit efficient solutions.


Figure 5: Illustration of a 2-layer linear neural network.

9.1.2 Matrix sensing and rank-constrained optimization

Moving beyond PCA and CCA, a more nontrivial problem is matrix sensing in the presence of RIP. The analysis for this problem, while much simpler than other problems like phase retrieval and matrix completion — is representative of a typical strategy for analyzing problems of this kind.

Theorem 33 (Global landscape for matrix sensing [30, 31]). Consider the matrix sensing problem (20). If A satisfies  2r-RIP with δ2r < 1/10, then:

(All local minima are global): for any local minimum  X of f(·), it satisfies  XX⊤ = M⋆;

(Strict saddles): for any critical point X that is not a local minimum, it satisfies  λmin�∇2f(X)�≤−4σr(M⋆)/5.

Proof of Theorem 33. For conciseness, we focus on the rank-1 case with  M⋆ = x⋆x⊤⋆ and show that all local minima are global. The complete proof can be found in [30, 31].

Consider any local minimum  x of f(·). This is characterized by the first-order and second-order optimality conditions


Without loss of generality, assume that  ∥x − x⋆∥2 ≤ ∥x + x⋆∥2.

A typical proof idea is to demonstrate that: if x is not globally optimal, then one can identify a descent direction, thus contradicting the local optimality of x. A natural guess of such a descent direction would be the direction towards the truth, i.e.  x − x⋆. As a result, the proof consists in showing that: when the RIP constant is sufficiently small, one has


unless  x = x⋆. Additionally, the value of (147) is helpful in upper bounding  λmin(∇2f(x))if x is a saddle point. See Appendix E for details.

We take a moment to expand on this result. Recall that we have introduced a version of strong convexity and smoothness for  f(·)when accounting for global orthogonal transformation (Section 5.1.2). Another way to express this is through a different parameterization


which clearly satisfies  g(XX⊤) = f(X). It is not hard to show that: in the presence of the RIP, the Hessian ∇2g(·)is well-conditioned when restricted to low-rank decision variables and directions. This motivates the following more general result, stated in terms of certain restricted well-conditionedness of  g(·). One of the advantages is its applicability to more general loss functions beyond the squared loss.

Theorem 34 (Global landscape for rank-constrained problems [31]). Let g(·) be a convex function. Suppose


admits a solution  M⋆with  rank(M⋆) = r < n. Assume that for all M and D with  rank(M) ≤ 2rand rank(D) ≤ 4r,


holds with  β/α ≤ 3/2. Then the function  f(X) = g(XX⊤), where X ∈ Rn×r, has no spurious local minima, and all saddle points of  f(·)are strict saddles.

Remark 11. This result continue to hold if  rank(M) ≤ 2˜r and rank(D) ≤ 4˜r for some ˜r > r(although the restricted well-conditionedness needs to be valid w.r.t. this new rank) [31]. In addition, it has also been extended to accommodate nonsymmetric matrix factorization in [31].

9.1.3 Phase retrieval and matrix completion

Next, we move on to problems that fall short of restricted well-conditionedness. As it turns out, it is still possible to have benign landscape, although the Lischiptz constants w.r.t. both gradients and Hessians might be much larger. A typical example is phase retrieval, for which the smoothness condition is not well-controlled (as discussed in Lemma 5).

Theorem 35 (Global landscape for phase retrieval [27]). Consider the phase retrieval problem (27). If the sample size  m ≳ n log3 n, then with high probability, there is no spurious local minimum, and all saddle points of  f(·)are strict saddles.

Further, we turn to the kind of loss functions that only satisfy highly restricted strong convexity and smoothness. In some cases, one might be able to properly regularize the loss function to enable benign landscape. Here, regularization can be enforced in a way similar to regularized gradient descent as discussed in Section 5.2.2. In the following, we use matrix completion as a representative example.

Theorem 36 (Global landscape for matrix completion [29, 159, 160]). Consider the problem (32) but replace f(·)with a regularized loss


with  λ > 0a regularization parameter, and  G0(z) = max{z − α, 0}4 for some α > 0. For properly selected  αand  λ, if the sample size  n2p ≳ n max{µκr log n, µ2κ2r2}, then with high probability, all local minima X of freg(·) satisfies XX⊤ = X⋆X⊤⋆ , and all saddle points of  freg(·)are strict saddles.

Remark 12. The study of global landscape in matrix completion was initiated in [29]. The current result in Theorem 36 is taken from [160].

9.1.4 Over-parameterization

Another scenario that merits special attention is over-parametrization [161, 162]. Take matrix sensing and phase retrieval for instance: if we lift the decision variable to the full matrix space (as opposed to using the low-rank decision variable), the resulting optimization problem is


where  M⋆is the true low-rank matrix. Here,  Ai = aia⊤i and M⋆ = x⋆x⊤⋆ for phase retrieval. As it turns out, under minimal sample complexity,  fop(·)does not have spurious local minima even though we over-parametrize the model significantly; in other words, enforcing the low-rank constraint is not crucial for optimization.

Theorem 37 (Over-parametrized matrix sensing and phase retrieval). Any local minimum  X⋆of the above function  fop(·) obeys X⋆X⊤⋆ = M⋆, provided that the set


is a singleton  {M⋆}.

The singleton property assumed in Theorem 37 has been established, for example, for the following two problems:

matrix sensing (the positive semidefinite case where  M⋆ ⪰ 0), as long as A satisfies 5r-RIP with a RIP constant  δ5r ≤ 1/10 [163]; a necessary and sufficient condition was established in [164];

phase retrieval, as long as  m ≥ c0nfor some sufficiently large constant  c0 > 0; see [165, 166].

As an important implication of Theorem 37, even solving the over-parametrized optimization problem allows us to achieve perfect recovery for the above two problems. This was first observed by [161] and then made precise by [162]; they showed that running GD w.r.t. the over-parameterized loss  fop(·)also (approximately) recovers the truth under roughly the same sample complexity, provided that a near-zero initialization is adopted.

Proof of Theorem 37. Define the function


which satisfies  g(XX⊤) = fop(X). Consider any local minimum  X⋆of  fop(·). By definition, there is an ϵ-ball Bϵ(X⋆) := {X | ∥X − X⋆∥F ≤ ϵ}with small enough  ϵ > 0 such that


Now suppose that  X⋆ (resp. X⋆X⊤⋆ ) is not a global solution of  fop(·) (resp. g(·)), then there exists a  ˜X ˜X⊤(resp. ˜X) arbitrarily close to  X⋆X⊤⋆ (resp. X⋆) such that


For instance, one can take ˜X ˜X⊤ = (1 − ζ)X⋆X⊤⋆ + ζX⋆X⊤⋆ with ζ ↘ 0. This implies the existence of an  ˜Xthat contradicts (152), thus establishing the claim.

9.1.5 Beyond low-rank matrix factorization

Finally, we note that benign landscape has been observed in numerous other contexts beyond matrix factorization. While they are beyond the scope of this paper, we briefly mention two important cases based on chronological order:

Dictionary learning [167169]. Observe a data matrix Y that can be factorized as Y = AX, where A is a square invertible dictionary matrix and X encodes the sparse coefficients. The goal is to learn A. It is shown that a certain smooth approximation of the  ℓ1loss exhibits benign nonconvex geometry over a sphere.

M-estimator in statistical estimation [170]. Given a set of independent data points  {x1, · · · , xn}, this work studied when the empirical loss function inherits the benign landscape of the population loss function. The results provide a general framework for establishing uniform convergence of the gradient and the Hessian of the empirical risk, and cover several examples such as binary classification, robust regression, and Gaussian mixture models.

9.1.6 Notes

The study of benign global landscapes dates back to the works on shallow neural networks in the 1980s [36] with deterministic data, if not earlier. Complete dictionary learning analyzed by Sun et al. [167, 168] was perhaps the first “modern” example where benign landscape is analyzed, which exploits measure concentration of random data. This work inspired the line of global geometric analysis for many aforementioned problems, including phase retrieval [27, 91], matrix sensing [30, 31, 159, 171, 172], matrix completion [29, 31, 159, 160], and robust PCA [159]. For phase retrieval, [27] focused on the smooth squared loss in (27). The landscape of a more “robust” nonsmooth formulation  f(x) = 1m�mi=1��(a⊤i x)2 − (a⊤i x⋆)2��has been studied by [91]. The optimization landscape for matrix completion was pioneered by [29], and later improved by [159, 160]. In particular, [160] derived a model-free theory where no assumptions are imposed on  M⋆, which accommodates, for example, the noisy case and the case where the truth is only approximately low-rank. The global landscape of asymmetric matrix sensing / completion holds similarly by considering a loss function regularized by the term  g(L, R) :=��L⊤L − R⊤R��2F[159, 173] or by the term  g(L, R) := ∥L∥2F + ∥R∥2F[31]. Theorem 34 has been generalized to the asymmetric case in [40]. Last but not least, we caution that, many more nonconvex problems are not benign and indeed have bad local minima; for example, spurious local minima are common even in simple neural networks with nonlinear activations [174, 175].

9.2 Gradient descent with random initialization

For many problems described above with benign landscape, there is no spurious local minima, and the only task is to escape strict saddle points and to find second-order critical points, which are now guaranteed to be global optima. In particular, our main algorithmic goal is to find a second-order critical point of a function exhibiting benign geometry. To make the discussion more precise, we define the functions of interest as follows

Definition 9 (strict saddle property [176, 177]). A function  f(·)is said to satisfy the  (ε, γ, ζ)-strict saddle property for some  ε, γ, ζ > 0, if for each x at least one of the following holds:


In words, this property says that: every point either has a large gradient, or has a negative directional curvature, or lies sufficiently close to a local minimum. In addition, while we have not discussed the strong gradient condition in the preceding subsection, it arises for most of the aforementioned problems when x is not close to the global minimum.

A natural question arises as to whether an algorithm as simple as gradient descent can converge to a second-order critical point of a function satisfying this property. Apparently, GD cannot start from anywhere; for example, if it starts from any undesired critical point (which obeys  ∇f(x) = 0), then it gets trapped. But what happens if we initialize GD randomly?

A recent work [178] provides the first answer to this question. Borrowing tools from dynamical systems theory, it proves that:

Theorem 38 (Convergence of GD with random initialization [179]). Consider any twice continuously differ-entiable function f that satisfies the strict saddle property (Definition 9). If  ηt < 1/β with βthe smoothness parameter, then GD with a random initialization converges to a local minimizer or  −∞almost surely.

This theorem says that for a broad class of benign functions of interest, GD — when randomly initialized — never gets stuck in the saddle points. The following example helps develop a better understanding of this theorem.

Example 4. Consider a quadratic minimization problem


The GD rule is  xt+1 = xt − ηtAxt. If ηt ≡ η < 1/∥A∥, then


Now suppose that  λ1(A) ≥ · · · ≥ λn−1(A) > 0 > λn(A), and let E+be the subspace spanned by the first  n − 1eigenvectors. It is easy to see that 0 is a saddle point, and that  xt → 0only if  x0 ∈ E+. In fact, as long as x0contains a component outside  E+, then this component will keep blowing up at a rate  1 + η|λn(A)|. Giventhat  E+ is (n − 1)-dimensional, we have  P(x0 ∈ E+) = 0. As a result,  P(xt → 0) = 0.

The above theory has been extended to accommodate other optimization methods like coordinate descent, mirror descent, the gradient primal-dual algorithm, and alternating minimization [179181]. In addition, the above theory is generic and accommodates all benign functions satisfying the strict saddle property.

We caution, however, that almost-sure convergence does not imply fast convergence. In fact, there exist non-pathological functions such that randomly initialized GD takes time exponential in the ambient dimension to escape saddle points [182]. That being said, it is possible to develop problem-specific theory that reveals much better convergence guarantees. Once again, we take phase retrieval as an example.

Theorem 39 (Randomly initialized GD for phase retrieval [75]). Consider the problem (27), and suppose that m ≳ npoly log m. The GD iterates with random initialization  x0 ∼ N(0, ∥x⋆∥22n In) and ηt ≡ 1/(c3∥x⋆∥22) obey


with probability  1 − O(n−10). Here, c3, c4 > 0are some constants,  T0 ≲ log n, and we assume  ∥x0 − x⋆∥2 ≤∥x0 + x⋆∥2.

To be more precise, the algorithm consists of two stages:

 When 0 ≤ t ≤ T0 ≲ log n: this stage allows GD to find and enter the local basin surrounding the truth, which takes time no more than O(log n) steps. To explain why this is fast, we remark that the signal strength  |⟨xt, x⋆⟩|in this stage grows exponentially fast, while the residual strength��xt − |⟨xt,x⋆⟩|∥x⋆∥22 x⋆��2does not increase by much.

 When t > T0: once the iterates enters the local basin, the  ℓ2estimation error decays exponentially fast, similar to Theorem 6. This stage takes about  O(log 1ε)iterations to reach  ε-accuracy.

Taken collectively, GD with random initialization achieves  ε-accuracy in  O�log n + log 1ε�iterations, making it appealing for solving large-scale problems. It is worth noting that the GD iterates never approach or hit the saddles; in fact, there is often a positive force dragging the GD iterates away from the saddle points.

Furthermore, there are other examples for which randomly initialized GD converges fast; see [183, 184] for further examples. We note, however, that the theoretical support for random initialization is currently lacking for many important problems (including matrix completion and blind deconvolution).

9.3 Generic saddle-escaping algorithms

Given that gradient descent with random initialization has only been shown to be efficient for specific examples, it is natural to ask how to design generic optimization algorithms to efficiently escape saddle points for all functions with benign geometry (i.e. those satisfying the strict saddle property in Definition 9). To see why this is hopeful, consider any strict saddle point x (which obeys  ∇f(x) = 0). The Taylor expansion implies


for any  ∆sufficiently small. Since x is a strict saddle, one can identify a direction of negative curvature and further decrease the objective value (i.e.  f(x + ∆) < f(x)). In other words, the existence of negative curvatures enables efficient escaping from undesired saddle points.

Many algorithms have been proposed towards the above goal. Roughly speaking, the available algorithms can be categorized into three classes, depending on which basic operations are needed: (1) Hessian-based algorithms; (2) Hessian-vector-product-based algorithms; (3) gradient-based algorithms.

Remark 13. Caution needs to be exercised as this categorization is very rough at best. One can certainly argue that many Hessian-based operations can be carried out via Hessian-vector products, and that Hessian-vector products can be computed using only gradient information.

Before proceeding, we introduce two notations that are helpful in presenting the theory. First, following the convention [185], a point x is said to be an  ε-second-order critical point if


Next, the algorithms often require the loss functions to be sufficiently smooth in the following sense:

Definition 10 (Hessian Lipschitz continuity). A function  f(·) is L2-Hessian-Lipschitz if, for any  x1, x2, onehas


9.3.1 Hessian-based algorithms

This class of algorithms requires an oracle that returns  ∇2f(x)for any given x. Examples include the trust-region method and the cubic-regularized Newton method.

Trust-region methods (TRM). At the tth iterate, this method minimizes a quadratic approximation of the loss function in a local neighborhood, where the quadratic approximation is typically based on a second-order Taylor expansion at the current iterate [186, 187]. Mathematically,


for some properly chosen radius  ∆. By choosing the size of the local neighborhood — called the “trust region” — to be sufficiently small, we can expect the second-order approximation within this region to be reasonably reliable. Thus, if  xthappens to be a strict saddle point, then TRM is still able to find a descent direction due to the negative curvature condition. This means that such local search will not get stuck in undesired saddle points.

Cubic regularization. This method attempts to minimize a cubic-regularized upper bound on the loss function [185]:


Here,  L2is taken to be the Lipschitz constant of the Hessians (see Definition 10), so as to ensure that the objective function in (157) majorizes the true objective  f(·). While the subproblem (157) is nonconvex and may have local minima, it can often be efficiently solved by minimizing an explicitly written univariate convex function [185, Section 5], or even by gradient descent [188].

Both methods are capable of finding an  ε-second-order stationary point (cf. (154)) in  O(ε−1.5)iterations [185, 189], where each iteration consists of one gradient and one Hessian computations.

9.3.2 Hessian-vector-product-based algorithms

Here, the oracle returns  ∇2f(x) ufor any given x and u. In what follows, we review the method by Carmon et al. [190] which falls under this category. There are two basic subroutines that are carried out using Hessian-vector products.

Negative-Curvature-Descent, which allows one to find, if possible, a direction that decreases the objective value. This is achieved by computing the eigenvector corresponding to the smallest eigenvalue of the Hessian.

Almost-Convex-AGD, which invokes Nesterov’s accelerated gradient method [191] to optimize an almost convex function superimposed with a squared proximity term.

In each iteration, based on the estimate of the smallest eigenvalue of the current Hessian, the algorithm decides whether to move along the direction computed by Negative-Curvature-Descent (so as not to get trapped in saddle points), or to apply Almost-Convex-AGD to optimize an almost convex function. See [190] for details. This method converges to an  ε-second-order stationary point in  O(ε−7/4 log 1ε)steps, where each step involves computing a Hessian-vector product.

Another method that achieves about the same computational complexity is Agarwal et al. [192]. This is a fast variant of cubic regularization. The key idea is to invoke, among other things, fast multiplicative approximations to accelerate the subproblem of the cubic-regularized Newton step. The interested readers shall consult [192] for details.

9.3.3 Gradient-based algorithms

Here, an oracle outputs  ∇f(x)for any given x. These methods are computationally efficient since only first-order information is required.

Ge et al. [176] initiated this line of work by designing a first-order algorithm that provably escapes strict saddle points in polynomial time with high probability. The algorithm proposed therein is a noise-injected version of (stochastic) gradient descent, namely,


where  ζtis some noise sampled uniformly from a sphere, and  ∇f(xt)can also be replaced by a stochastic gradient. The iteration complexity, however, depends on the dimension polynomially, which grows prohibitively as the ambient dimension of  xtincreases. Similar high iteration complexity holds for [193] which is based on injecting noise to normalized gradient descent.

The computational guarantee was later improved by perturbed gradient descent, proposed in [194]. In contrast to (158), perturbed GD adds noise to the iterate before computing the gradient, namely,


with  ζtuniformly drawn from a sphere. The crux of perturbed GD is to realize that strict saddles are unstable; it is possible to escape them and make progress with a slight perturbation. It has been shown that perturbed GD finds an  ε-second-order stationary point in  O(ε−2)iterations (up to some logarithmic factor) with high probability, and is hence almost dimension-free. Note that each iteration only needs to call the gradient oracle once.

Moreover, this can be further accelerated via Nesterov’s momentum-based methods [191]. Specifically, [195] proposed a perturbed version of Nesterov’s accelerated gradient descent (AGD), which adds perturbation to AGD and combines it with an operation similar to Negative-Curvature-Descent described above. This accelerated method converges to an  ε-second-order stationary point in  O(ε−7/4)iterations (up to some logarithmic factor) with high probability, which matches the computational complexity of [190, 192].

There are also several other algorithms that provably work well in the presence of stochastic / incremental first-order and / or second-order oracles [196198]. Interested readers are referred to the ICML tutorial [199] and the references therein. Additionally, we note that it is possible to adapt many of the above mentioned saddle-point escaping algorithms onto manifolds; we refer the readers to [168, 187, 200].

9.3.4 Caution

Finally, it is worth noting that: in addition to the dependency on  ε, all of the above iteration complexities also rely on (1) the smoothness parameter, (2) the Lipschitz constant of the Hessians, and (3) the local strong convexity parameter. These parameters, however, do not necessarily fall within the desired level. For instance, in most of the problems mentioned herein (e.g. phase retrieval, matrix completion), the nonconvex objective function is not globally smooth, namely, the local gradient Lipschitz constant might grow to infinity as the parameters become unbounded. Motivated by this, one might need to further assume that optimization is performed within a bounded set of parameters, as is the case for [170]. Leaving out this boundedness issue, a more severe matter is that these parameters often depend on the problem size. For example, recall that in

Table 1: A selective summary of problems with provable nonconvex solutions


phase retrieval, the smoothness parameter scales with n even within a reasonably small bounded range (see Lemma 5). As a consequence, the resulting iteration complexities of the above saddle-escaping algorithm might be significantly larger than, say,  O(ε−7/4)for various large-scale problems.

Given that nonconvex matrix factorization is a rapidly growing field, we expect that many of the results reviewed herein may be improved or superseded. Nonetheless, the core techniques and insights reviewed in this article will continue to play a key role in understanding nonconvex statistical estimation and learning. While a paper of this length is impossible to cover all recent developments, our aim has been to convey the key ideas via representative examples and to strike a balance between useful theory and varied algorithms. Table 1 summarizes representative references on the canonical problems reviewed herein, with a pointer to their corresponding convex approaches for completeness. It should be noted that most of these results are obtained under random data models, which may not be satisfied in practice (e.g., in matrix completion the sampling patterns can be non-uniform). This means one should not take directly the performance guarantees as is before checking the reasonability of assumptions, and indeed, the performance may drop significantly without carefully modifying the nonconvex approaches on a case-by-case basis to specific problems.

Caution needs to be exercised that the theory derived for nonconvex paradigms may still be sub-optimal in several aspects. For instance, the number of samples required for recovering / completing a rank-r matrix often scales at least quadratically with the rank r in the current nonconvex optimization theory, which is outperformed by convex relaxation (whose sample complexity typically scales linearly with r). The theoretical dependency on the condition number of the matrix also falls short of optimality.

Due to the space limits, we have omitted developments in several other important aspects, most notably the statistical guarantees vis-à-vis noise. Most of our discussions (and analyses) can be extended without much difficulty to the regime with small to moderate noise; see, for example, [2022, 65, 77] for the stability of gradient methods in phase retrieval, matrix completion, and blind deconvolution. Encouragingly, the nonconvex methods not only allow to control the Euclidean estimation errors in a minimax optimal manner, they also provably achieve optimal statistical accuracy in terms of finer error metrics (e.g. optimal entrywise error control in matrix completion [26]). In addition, there is an intimate link between convex relaxation and nonconvex optimization, which allows us to improve the stability guarantees of convex relaxation via the theory of nonconvex optimization; see [77, 78] for details.

Furthermore, we have also left out several important nonconvex problems and methods in order to stay focused, including but not limited to (1) blind calibration and finding sparse vectors in a subspace [97, 220, 221]; (2) tensor completion, decomposition and mixture models [156, 222225]; (3) parameter recovery and inverse problems in shallow and deep neural networks [184, 226230]; (4) analysis of Burer-Monteiro factorization to semidefinite programs [55, 231233]. The interested readers can consult [234] for an extensive list of further references. We would also like to refer the readers to an excellent recent monograph by Jain and Kar [235] that complements our treatment. Roughly speaking, [235] concentrates on the introduction of generic optimization algorithms, with low-rank factorization (and other learning problems) used as special examples to illustrate the generic results. In contrast, the current article focuses on unveiling deep statistical and algorithmic insights specific to nonconvex low-rank factorization.

Finally, we conclude this paper by singling out several exciting avenues for future research:

investigate the convergence rates of randomly initialized gradient methods for problems without Gaussian-type measurements (e.g. matrix completion and blind deconvolution);

characterize generic landscape properties that enable fast convergence of gradient methods from random initialization;

relax the stringent assumptions on the statistical models underlying the data; for example, a few recent works studied nonconvex phase retrieval under more physically-meaningful measurement models [201, 202];

develop robust and scalable nonconvex methods that can handle distributed data samples with strong statistical guarantees;

study the capability of other commonly used optimization algorithms (e.g. alternating minimization) in escaping saddle points;

characterize the optimization landscape of constrained nonconvex statistical estimation problems like non-negative low-rank matrix factorization (the equality-constrained case has recently been explored by [236]);

develop a unified and powerful theory that can automatically accommodate the problems considered herein and beyond.

The authors thank the editors of IEEE Transactions on Signal Processing, Prof. Pier Luigi Dragotti and Prof. Jarvis Haupt, for handling this invited overview article, and thank the editorial board for useful feedbacks of the white paper. The authors thank Cong Ma, Botao Hao, Tian Tong and Laixi Shi for proofreading an early version of the manuscript.

To begin with, simple calculation reveals that for any  z, x ∈ Rn, the Hessian obeys (see also [30, Lemma 4.3])


With the notation (23) at hand, we can rewrite


where the last line uses the symmetry of  Ai.To establish local strong convexity and smoothness, we need to control  z⊤∇2f(x)zfor all z. The key ingredient is to show that:  z⊤∇2f(x)zis not too far away from


a function that can be easily shown to be locally strongly convex and smooth. To this end, we resort to the RIP (see Definition 6). When A satisfies 4-RIP, Lemma 3 indicates that


where the last line holds if  ∥x − x⋆∥2 ≤ ∥x⋆∥2. Similarly,


As a result, if  δ4is small enough, then putting the above results together implies that:  z⊤∇2f(x)zis sufficiently close to g(x, z). A little algebra then gives (see Appendix B)


for all z, which provides bounds on local strong convexity and smoothness parameters. Applying Lemma 1 thus establishes the theorem.

Without loss of generality, we assume that  ∥z∥2 = ∥x⋆∥2 = 1. When  ∥x − x⋆∥2 ≤ ∥x⋆∥2, some elementary algebra gives


In addition,


Putting these together with (160), (162) and (163), we have


If  δ4 ≤ 1/44 and ∥x − x⋆∥2 ≤ 1/12, then we arrive at


Here, we demonstrate that when X is sufficiently close to  X⋆, then the objective function of  f∞(·)of (46) exhibits the modified form of strong convexity in the form of (51). Set  V = ZHZ − X⋆. In view of (47), it suffices to show


for all X sufficiently close to  X⋆and all Z. To this end, we first observe a simple perturbation bound (the proof is omitted)


for some universal constant  c1 > 0, provided that  ∥X − X⋆∥Fis small enough. We then turn attention to g(X⋆, V ):


The last line holds by recognizing that  X⊤⋆ ZHZ ⪰ 0 [53, Theorem 2], which implies that  Tr(X⊤⋆ ZHZX⊤⋆ ZHZ) ≥0 and Tr�X⊤⋆ ZHZX⊤⋆ X⋆�≥ 0. Thus,


as long as  ∥X − X⋆∥F ≤ σr(X⊤⋆ X⋆)∥V ⊤V ∥2F2c1∥X⋆∥F·∥V ∥F. In summary,


We start with the rank-r PSD case, where we denote by  X0the initial estimate (i.e.  X0 = L0 = R0 for thiscase) and  X⋆the ground truth. Observe that


where the second line follows since  X0X⊤0 is the best rank-r approximation of Y (and hence  ∥X0X⊤0 − Y ∥ ≤∥X⋆X⊤⋆ − Y ∥), and the last inequality follows from Lemma 8. A useful lemma from [23, Lemma 5.4] allows us to directly upper bound  dist2(X0, X⋆)by the Euclidean distance between their low-rank counterparts.

Lemma 9 ([23]). For any U, X ∈ Rn×r, we have


As an immediate consequence of the above lemma, we get


where the last line follows from (165). Therefore, as long as  δ2r ≪ √ζ/(rκ), we have


In fact, [127] shows a stronger consequence of the RIP: one can improve the left-hand-side of (165) to ∥X0X⊤0 − M⋆∥F, which allows relaxing the requirement on the RIP constant to  δ2r ≪ √ζ/(√rκ) [23]. Theasymmetric case can be proved in a similar fashion.

Before proceeding, we first single out an immediate consequence of the first-order optimality condition (146a) that will prove useful. Specifically, for any critical point x,


This fact basically says that any critical point of  f(·)stays very close to the truth in the subpace spanned by this point. To verify (147), we observe the identity


where the last line follows from (146a) with a little algebra. This combined with (25) gives


We then need to show that (168) is negative. To this end, we quote an elementary algebraic inequality [30, Lemma 4.4]


This taken collectively with (168) and (167) yields


which is negative unless  x = x⋆. This concludes the proof.

Proof of Claim (167). To see why (167) holds, observe that


Therefore, for any Z, using (25) we have


[1] E. J. Candès and B. Recht, “Exact matrix completion via convex optimization,” Foundations of Computational Mathematics, vol. 9, no. 6, pp. 717–772, 2009.

[2] M. A. Davenport and J. Romberg, “An overview of low-rank matrix recovery from incomplete observations,” IEEE Journal of Selected Topics in Signal Processing, vol. 10, no. 4, pp. 608–622, 2016.

[3] Y. Chen and Y. Chi, “Harnessing structures in big data via guaranteed low-rank matrix estimation: Recent theory and fast algorithms via convex and nonconvex optimization,” IEEE Signal Processing Magazine, vol. 35, no. 4, pp. 14 – 31, 2018.

[4] Y. Shechtman, Y. C. Eldar, O. Cohen, H. N. Chapman, J. Miao, and M. Segev, “Phase retrieval with application to optical imaging: a contemporary overview,” IEEE signal processing magazine, vol. 32, no. 3, pp. 87–109, 2015.

[5] A. Ahmed, B. Recht, and J. Romberg, “Blind deconvolution using convex programming,” IEEE Transactions on Information Theory, vol. 60, no. 3, pp. 1711–1732, 2014.

[6] S. Ling and T. Strohmer, “Self-calibration and biconvex compressive sensing,” Inverse Problems, vol. 31, no. 11, p. 115002, 2015.

[7] V. Chandrasekaran, S. Sanghavi, P. Parrilo, and A. Willsky, “Rank-sparsity incoherence for matrix decomposition,” SIAM Journal on Optimization, vol. 21, no. 2, pp. 572–596, 2011.

[8] E. J. Candès, X. Li, Y. Ma, and J. Wright, “Robust principal component analysis?” Journal of the ACM, vol. 58, no. 3, pp. 11:1–11:37, 2011.

[9] A. Singer, “Angular synchronization by eigenvectors and semidefinite programming,” Applied and computational harmonic analysis, vol. 30, no. 1, pp. 20–36, 2011.

[10] Y. Chen and E. Candès, “The projected power method: An efficient algorithm for joint alignment from pairwise differences,” Communications on Pure and Applied Mathematics, vol. 71, no. 8, pp. 1648–1714, 2018.

[11] M. Fazel, “Matrix rank minimization with applications,” Ph.D. dissertation, Stanford University, 2002.

[12] B. Recht, M. Fazel, and P. A. Parrilo, “Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization,” SIAM review, vol. 52, no. 3, pp. 471–501, 2010.

[13] E. Candès and T. Tao, “The power of convex relaxation: Near-optimal matrix completion,” IEEE Transactions on Information Theory, vol. 56, no. 5, pp. 2053 –2080, May 2010.

[14] Y. Koren, R. Bell, and C. Volinsky, “Matrix factorization techniques for recommender systems,” Computer, no. 8, pp. 30–37, 2009.

[15] P. Chen and D. Suter, “Recovering the missing components in a large noisy low-rank matrix: Application to SFM,” IEEE transactions on pattern analysis and machine intelligence, vol. 26, no. 8, pp. 1051–1063, 2004.

[16] R. H. Keshavan, A. Montanari, and S. Oh, “Matrix completion from a few entries,” IEEE Transactions on Information Theory, vol. 56, no. 6, pp. 2980 –2998, June 2010.

[17] P. Jain, P. Netrapalli, and S. Sanghavi, “Low-rank matrix completion using alternating minimization,” in Proceedings of the forty-fifth annual ACM symposium on Theory of computing. ACM, 2013, pp. 665–674.

[18] E. Candès, X. Li, and M. Soltanolkotabi, “Phase retrieval via Wirtinger flow: Theory and algorithms,” Information Theory, IEEE Transactions on, vol. 61, no. 4, pp. 1985–2007, 2015.

[19] R. Sun and Z.-Q. Luo, “Guaranteed matrix completion via non-convex factorization,” IEEE Transactions on Information Theory, vol. 62, no. 11, pp. 6535–6579, 2016.

[20] Y. Chen and E. Candès, “Solving random quadratic systems of equations is nearly as easy as solving linear systems,” Communications on Pure and Applied Mathematics, vol. 70, no. 5, pp. 822–883, 2017.

[21] Y. Chen and M. J. Wainwright, “Fast low-rank estimation by projected gradient descent: General statistical and algorithmic guarantees,” arXiv preprint arXiv:1509.03025, 2015.

[22] X. Li, S. Ling, T. Strohmer, and K. Wei, “Rapid, robust, and reliable blind deconvolution via nonconvex optimization,” Applied and Computational Harmonic Analysis, 2018.

[23] S. Tu, R. Boczar, M. Simchowitz, M. Soltanolkotabi, and B. Recht, “Low-rank solutions of linear matrix equations via Procrustes flow,” in International Conference Machine Learning, 2016, pp. 964–973.

[24] H. Zhang, Y. Chi, and Y. Liang, “Provable non-convex phase retrieval with outliers: Median truncated Wirtinger flow,” in International conference on machine learning, 2016, pp. 1022–1031.

[25] P. Netrapalli, U. Niranjan, S. Sanghavi, A. Anandkumar, and P. Jain, “Non-convex robust PCA,” in Advances in Neural Information Processing Systems, 2014, pp. 1107–1115.

[26] C. Ma, K. Wang, Y. Chi, and Y. Chen, “Implicit regularization in nonconvex statistical estimation: Gradient descent converges linearly for phase retrieval, matrix completion and blind deconvolution,” arXiv preprint arXiv:1711.10467, accepted to Foundations of Computational Mathematics, 2017.

[27] J. Sun, Q. Qu, and J. Wright, “A geometric analysis of phase retrieval,” Foundations of Computational Mathematics, pp. 1–68, 2018.

[28] ——, “Complete dictionary recovery using nonconvex optimization,” in Proceedings of the 32nd International Conference on Machine Learning, 2015, pp. 2351–2360.

[29] R. Ge, J. D. Lee, and T. Ma, “Matrix completion has no spurious local minimum,” in Advances in Neural Information Processing Systems, 2016, pp. 2973–2981.

[30] S. Bhojanapalli, B. Neyshabur, and N. Srebro, “Global optimality of local search for low rank matrix recovery,” in Advances in Neural Information Processing Systems, 2016, pp. 3873–3881.

[31] Q. Li, Z. Zhu, and G. Tang, “The non-convex geometry of low-rank matrix optimization,” Information and Inference: A Journal of the IMA, 2018, in press.

[32] S. Bubeck, “Convex optimization: Algorithms and complexity,”  Foundations and Trends R⃝ in MachineLearning, vol. 8, no. 3-4, pp. 231–357, 2015.

[33] S. Lang, “Real and functional analysis,” Springer-Verlag, New York,, vol. 10, pp. 11–13, 1993.

[34] A. Beck, First-Order Methods in Optimization. SIAM, 2017, vol. 25.

[35] Y. Nesterov, Introductory lectures on convex optimization: A basic course. Springer Science & Business Media, 2013, vol. 87.

[36] P. Baldi and K. Hornik, “Neural networks and principal component analysis: Learning from examples without local minima,” Neural networks, vol. 2, no. 1, pp. 53–58, 1989.

[37] P. F. Baldi and K. Hornik, “Learning in linear neural networks: A survey,” IEEE Transactions on neural networks, vol. 6, no. 4, pp. 837–858, 1995.

[38] B. Yang, “Projection approximation subspace tracking,” IEEE Transactions on Signal processing, vol. 43, no. 1, pp. 95–107, 1995.

[39] X. Li, Z. Wang, J. Lu, R. Arora, J. Haupt, H. Liu, and T. Zhao, “Symmetry, saddle points, and global geometry of nonconvex matrix factorization,” arXiv preprint arXiv:1612.09296, 2016.

[40] Z. Zhu, Q. Li, G. Tang, and M. B. Wakin, “Global optimality in low-rank matrix optimization,” IEEE Transactions on Signal Processing, vol. 66, no. 13, pp. 3614–3628, 2018.

[41] Z. Zhu, D. Soudry, Y. C. Eldar, and M. B. Wakin, “The global optimization geometry of shallow linear neural networks,” arXiv preprint arXiv:1805.04938, 2018.

[42] R. A. Hauser, A. Eftekhari, and H. F. Matzinger, “PCA by determinant optimisation has no spurious local optima,” in Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining. ACM, 2018, pp. 1504–1511.

[43] E. J. Candès and Y. Plan, “Tight oracle inequalities for low-rank matrix recovery from a minimal number of noisy random measurements,” IEEE Transactions on Information Theory, vol. 57, no. 4, pp. 2342–2359, 2011.

[44] E. J. Candès, “The restricted isometry property and its implications for compressed sensing,” Comptes Rendus Mathematique, vol. 346, no. 9, pp. 589–592, 2008.

[45] J. R. Fienup, “Phase retrieval algorithms: a comparison.” Applied optics, vol. 21, no. 15, pp. 2758–2769, 1982.

[46] E. J. Candès, T. Strohmer, and V. Voroninski, “Phaselift: Exact and stable signal recovery from magnitude measurements via convex programming,” Communications on Pure and Applied Mathematics, vol. 66, no. 8, pp. 1241–1274, 2013.

[47] K. Jaganathan, Y. C. Eldar, and B. Hassibi, “Phase retrieval: An overview of recent developments,” arXiv preprint arXiv:1510.07713, 2015.

[48] Y. Chen, Y. Chi, and A. Goldsmith, “Exact and stable covariance estimation from quadratic sampling via convex programming,” IEEE Transactions on Information Theory, vol. 61, no. 7, pp. 4034–4059, July 2015.

[49] T. Cai and A. Zhang, “ROP: Matrix recovery via rank-one projections,” The Annals of Statistics, vol. 43, no. 1, pp. 102–138, 2015.

[50] L. Tian, J. Lee, S. B. Oh, and G. Barbastathis, “Experimental compressive phase space tomography,” Optics Express, vol. 20, no. 8, p. 8296, 2012.

[51] C. Bao, G. Barbastathis, H. Ji, Z. Shen, and Z. Zhang, “Coherence retrieval using trace regularization,” SIAM Journal on Imaging Sciences, vol. 11, no. 1, pp. 679–706, 2018.

[52] N. Vaswani, T. Bouwmans, S. Javed, and P. Narayanamurthy, “Robust subspace learning: Robust PCA, robust subspace tracking, and robust subspace recovery,” IEEE Signal Processing Magazine, vol. 35, no. 4, pp. 32–55, 2018.

[53] J. M. Ten Berge, “Orthogonal Procrustes rotation for two or more matrices,” Psychometrika, vol. 42, no. 2, pp. 267–276, 1977.

[54] Q. Zheng and J. Lafferty, “A convergent gradient descent algorithm for rank minimization and semidefi-nite programming from random linear measurements,” arXiv preprint arXiv:1506.06081, 2015.

[55] S. Bhojanapalli, A. Kyrillidis, and S. Sanghavi, “Dropping convexity for faster semi-definite optimization,” in Conference on Learning Theory, 2016, pp. 530–582.

[56] J. A. Tropp, “An introduction to matrix concentration inequalities,”  Foundations and Trends R⃝ inMachine Learning, vol. 8, no. 1-2, pp. 1–230, 2015.

[57] R. Vershynin, “Introduction to the non-asymptotic analysis of random matrices,” in Compressed Sensing, Y. C. Eldar and G. Kutyniok, Eds. Cambridge University Press, 2012, pp. 210–268.

[58] M. Soltanolkotabi, “Algorithms and theory for clustering and nonconvex quadratic programming,” Ph.D. dissertation, Stanford University, 2014.

[59] S. Sanghavi, R. Ward, and C. D. White, “The local convexity of solving systems of quadratic equations,” Results in Mathematics, vol. 71, no. 3-4, pp. 569–608, 2017.

[60] Y. Li, C. Ma, Y. Chen, and Y. Chi, “Nonconvex matrix factorization from rank-one measurements,” in Proceedings of International Conference on Artificial Intelligence and Statistics (AISTATS), vol. 89. PMLR, 16–18 Apr 2019, pp. 1496–1505.

[61] S. Ling and T. Strohmer, “Regularized gradient descent: a non-convex recipe for fast joint blind deconvolution and demixing,” Information and Inference: A Journal of the IMA, vol. 8, no. 1, pp. 1–49, 2018.

[62] J. Dong and Y. Shi, “Nonconvex demixing from bilinear measurements,” IEEE Transactions on Signal Processing, vol. 66, no. 19, pp. 5152–5166, 2018.

[63] Q. Zheng and J. Lafferty, “Convergence analysis for rectangular matrix completion using Burer-Monteiro factorization and gradient descent,” arXiv preprint arXiv:1605.07051, 2016.

[64] X. Yi, D. Park, Y. Chen, and C. Caramanis, “Fast algorithms for robust PCA via gradient descent,” in Advances in neural information processing systems, 2016, pp. 4152–4160.

[65] R. H. Keshavan, A. Montanari, and S. Oh, “Matrix completion from noisy entries,” The Journal of Machine Learning Research, vol. 11, pp. 2057–2078, July 2010.

[66] W. Huang and P. Hand, “Blind deconvolution by a steepest descent algorithm on a quotient manifold,” SIAM Journal on Imaging Sciences, vol. 11, no. 4, pp. 2757–2785, 2018.

[67] J. Chen, D. Liu, and X. Li, “Nonconvex rectangular matrix completion via gradient descent without ℓ2,∞regularization,” arXiv:1901.06116, 2019.

[68] C. Stein, “A bound for the error in the normal approximation to the distribution of a sum of dependent random variables,” in Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability, 1972.

[69] L. H. Chen, L. Goldstein, and Q.-M. Shao,  Normal approximation by Stein’s method.Springer Science & Business Media, 2010.

[70] N. El Karoui, “On the impact of predictor geometry on the performance on high-dimensional ridgeregularized generalized robust regression estimators,” Probability Theory and Related Fields, pp. 1–81, 2015.

[71] Y. Zhong and N. Boumal, “Near-optimal bounds for phase synchronization,” SIAM Journal on Optimization, vol. 28, no. 2, pp. 989–1016, 2018.

[72] P. Sur, Y. Chen, and E. J. Candès, “The likelihood ratio test in high-dimensional logistic regression is asymptotically a rescaled chi-square,” Probability Theory and Related Fields, pp. 1–72, 2017.

[73] Y. Chen, J. Fan, C. Ma, and K. Wang, “Spectral method and regularized MLE are both optimal for top-K ranking,” Annals of Statistics, vol. 47, no. 4, pp. 2204–2235, August 2019.

[74] E. Abbe, J. Fan, K. Wang, and Y. Zhong, “Entrywise eigenvector analysis of random matrices with low expected rank,” arXiv preprint arXiv:1709.09565, 2017.

[75] Y. Chen, Y. Chi, J. Fan, and C. Ma, “Gradient descent with random initialization: Fast global convergence for nonconvex phase retrieval,” Mathematical Programming, vol. 176, no. 1-2, pp. 5–37, July 2019.

[76] L. Ding and Y. Chen, “The leave-one-out approach for matrix completion: Primal and dual analysis,” arXiv preprint arXiv:1803.07554, 2018.

[77] Y. Chen, Y. Chi, J. Fan, C. Ma, and Y. Yan, “Noisy matrix completion: Understanding statistical guarantees for convex relaxation via nonconvex optimization,” arXiv preprint arXiv:1902.07698, 2019.

[78] Y. Chen, J. Fan, C. Ma, and Y. Yan, “Inference and uncertainty quantification for noisy matrix completion,” arXiv preprint arXiv:1906.04159, 2019.

[79] M. Soltanolkotabi, “Structured signal recovery from quadratic measurements: Breaking sample complexity barriers via nonconvex optimization,” IEEE Transactions on Information Theory, vol. 65, no. 4, pp. 2374–2400, 2019.

[80] T. T. Cai, X. Li, and Z. Ma, “Optimal rates of convergence for noisy sparse phase retrieval via thresholded Wirtinger flow,” The Annals of Statistics, vol. 44, no. 5, pp. 2221–2251, 2016.

[81] J. Duchi, S. Shalev-Shwartz, Y. Singer, and T. Chandra, “Efficient projections onto the  ℓ1-ball for learning in high dimensions,” in International conference on Machine learning, 2008, pp. 272–279.

[82] G. Wang, L. Zhang, G. B. Giannakis, M. Akçakaya, and J. Chen, “Sparse phase retrieval via truncated amplitude flow,” IEEE Transactions on Signal Processing, vol. 66, no. 2, pp. 479–491, 2018.

[83] X. Li and V. Voroninski, “Sparse signal recovery from quadratic measurements via convex programming,” SIAM Journal on Mathematical Analysis, vol. 45, no. 5, pp. 3019–3033, 2013.

[84] S. Oymak, A. Jalali, M. Fazel, Y. C. Eldar, and B. Hassibi, “Simultaneously structured models with application to sparse and low-rank matrices,” IEEE Transactions on Information Theory, vol. 61, no. 5, pp. 2886–2908, 2015.

[85] G. Jagatap and C. Hegde, “Fast, sample-efficient algorithms for structured phase retrieval,” in Advances in Neural Information Processing Systems, 2017, pp. 4917–4927.

[86] R. Kolte and A. Özgür, “Phase retrieval via incremental truncated Wirtinger flow,” arXiv preprint arXiv:1606.03196, 2016.

[87] G. Wang, G. B. Giannakis, and Y. C. Eldar, “Solving systems of random quadratic equations via truncated amplitude flow,” IEEE Transactions on Information Theory, vol. 64, no. 2, pp. 773–794, 2018.

[88] Y. Li, Y. Chi, H. Zhang, and Y. Liang, “Nonconvex low-rank matrix recovery with arbitrary outliers via median-truncated gradient descent,” arXiv preprint arXiv:1709.08114, 2017.

[89] F. H. Clarke, “Generalized gradients and applications,” Transactions of the American Mathematical Society, vol. 205, pp. 247–262, 1975.

[90] H. Zhang, Y. Zhou, Y. Liang, and Y. Chi, “A nonconvex approach for phase retrieval: Reshaped Wirtinger flow and incremental algorithms,” Journal of Machine Learning Research, vol. 18, no. 141, pp. 1–35, 2017.

[91] D. Davis, D. Drusvyatskiy, and C. Paquette, “The nonsmooth landscape of phase retrieval,” arXiv preprint arXiv:1711.03247, 2017.

[92] N. Boumal, “Nonconvex phase synchronization,” SIAM Journal on Optimization, vol. 26, no. 4, pp. 2355–2377, 2016.

[93] H. Liu, M.-C. Yue, and A. Man-Cho So, “On the estimation performance and convergence rate of the generalized power method for phase synchronization,” SIAM Journal on Optimization, vol. 27, no. 4, pp. 2426–2446, 2017.

[94] Y. Chen, C. Suh, and A. J. Goldsmith, “Information recovery from pairwise measurements,” IEEE Transactions on Information Theory, vol. 62, no. 10, October 2016.

[95] I. Waldspurger, A. d’Aspremont, and S. Mallat, “Phase recovery, maxcut and complex semidefinite programming,” Mathematical Programming, vol. 149, no. 1-2, pp. 47–81, 2015.

[96] X.-T. Yuan and T. Zhang, “Truncated power method for sparse eigenvalue problems,” Journal of Machine Learning Research, vol. 14, no. Apr, pp. 899–925, 2013.

[97] Y. Li, K. Lee, and Y. Bresler, “Blind gain and phase calibration for low-dimensional or sparse signal sensing via power iteration,” in Sampling Theory and Applications (SampTA), 2017 International Conference on. IEEE, 2017, pp. 119–123.

[98] A. Edelman, T. A. Arias, and S. T. Smith, “The geometry of algorithms with orthogonality constraints,” SIAM journal on Matrix Analysis and Applications, vol. 20, no. 2, pp. 303–353, 1998.

[99] P.-A. Absil, R. Mahony, and R. Sepulchre, Optimization algorithms on matrix manifolds. Princeton University Press, 2009.

[100] L. Balzano and S. J. Wright, “Local convergence of an algorithm for subspace identification from partial data,” Foundations of Computational Mathematics, pp. 1–36, 2014.

[101] W. Dai, E. Kerman, and O. Milenkovic, “A geometric approach to low-rank matrix completion,” IEEE Transactions on Information Theory, vol. 58, no. 1, pp. 237–247, 2012.

[102] K. Wei, J.-F. Cai, T. F. Chan, and S. Leung, “Guarantees of Riemannian optimization for low rank matrix recovery,” SIAM Journal on Matrix Analysis and Applications, vol. 37, no. 3, pp. 1198–1222, 2016.

[103] ——, “Guarantees of Riemannian optimization for low rank matrix completion,” arXiv preprint arXiv:1603.06610, 2016.

[104] B. Vandereycken, “Low-rank matrix completion by Riemannian optimization,” SIAM Journal on Optimization, vol. 23, no. 2, pp. 1214–1236, 2013.

[105] A. Uschmajew and B. Vandereycken, “On critical points of quadratic low-rank matrix optimization problems,” 2018.

[106] J.-F. Cai and K. Wei, “Solving systems of phaseless equations via Riemannian optimization with optimal sampling complexity,” arXiv preprint arXiv:1809.02773, 2018.

[107] ——, “Exploiting the structure effectively and efficiently in low rank matrix recovery,” arXiv preprint arXiv:1809.03652, 2018.

[108] H. Robbins and S. Monro, “A stochastic approximation method,” in Herbert Robbins Selected Papers. Springer, 1985, pp. 102–109.

[109] A. Nemirovski, A. Juditsky, G. Lan, and A. Shapiro, “Robust stochastic approximation approach to stochastic programming,” SIAM Journal on optimization, vol. 19, no. 4, pp. 1574–1609, 2009.

[110] L. Bottou, F. E. Curtis, and J. Nocedal, “Optimization methods for large-scale machine learning,” SIAM Review, vol. 60, no. 2, pp. 223–311, 2018.

[111] K. Wei, “Solving systems of phaseless equations via Kaczmarz methods: A proof of concept study,” Inverse Problems, vol. 31, no. 12, p. 125008, 2015.

[112] G. Li, Y. Gu, and Y. M. Lu, “Phase retrieval using iterative projections: Dynamics in the large systems limit,” in 53rd Annual Allerton Conference on Communication, Control, and Computing, 2015, pp. 1114–1118.

[113] Y. Chi and Y. M. Lu, “Kaczmarz method for solving quadratic equations,” IEEE Signal Processing Letters, vol. 23, no. 9, pp. 1183–1187, 2016.

[114] Y. S. Tan and R. Vershynin, “Phase retrieval via randomized Kaczmarz: Theoretical guarantees,” Information and Inference: A Journal of the IMA, vol. 8, no. 1, pp. 97–123, 2018.

[115] V. Monardo, Y. Li, and Y. Chi, “Solving quadratic equations via amplitude-based nonconvex optimization,” in 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2019, pp. 5526–5530.

[116] C. Jin, S. M. Kakade, and P. Netrapalli, “Provable efficient online matrix completion via non-convex stochastic gradient descent,” in Advances in Neural Information Processing Systems, 2016, pp. 4520– 4528.

[117] L. N. Trefethen and D. Bau III, Numerical linear algebra. SIAM, 1997, vol. 50.

[118] P. Netrapalli, P. Jain, and S. Sanghavi, “Phase retrieval using alternating minimization,” IEEE Transactions on Signal Processing, vol. 18, no. 63, pp. 4814–4826, 2015.

[119] I. Waldspurger, “Phase retrieval with random gaussian sensing vectors by alternating projections,” IEEE Transactions on Information Theory, vol. 64, no. 5, pp. 3301–3312, May 2018.

[120] R. W. Gerchberg, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik, vol. 35, p. 237, 1972.

[121] T. Hastie, R. Mazumder, J. D. Lee, and R. Zadeh, “Matrix completion and low-rank SVD via fast alternating least squares.” Journal of Machine Learning Research, vol. 16, pp. 3367–3402, 2015.

[122] R. H. Keshavan, “Efficient algorithms for collaborative filtering,” Ph.D. dissertation, Stanford University, 2012.

[123] M. Hardt, “Understanding alternating minimization for matrix completion,” in Foundations of Computer Science (FOCS), 2014 IEEE 55th Annual Symposium on. IEEE, 2014, pp. 651–660.

[124] M. Hardt and M. Wootters, “Fast matrix completion without the condition number,” in Proceedings of The 27th Conference on Learning Theory, 2014, pp. 638–678.

[125] T. Zhao, Z. Wang, and H. Liu, “Nonconvex low rank matrix factorization via inexact first order oracle.” Advances in Neural Information Processing Systems, 2015.

[126] P. Jain, R. Meka, and I. S. Dhillon, “Guaranteed rank minimization via singular value projection,” in Advances in Neural Information Processing Systems, 2010, pp. 937–945.

[127] S. Oymak, B. Recht, and M. Soltanolkotabi, “Sharp time–data tradeoffs for linear inverse problems,” IEEE Transactions on Information Theory, vol. 64, no. 6, pp. 4129–4158, 2018.

[128] N. Halko, P.-G. Martinsson, and J. A. Tropp, “Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions,” SIAM Revieweview, vol. 53, no. 2, pp. 217–288, 2011.

[129] J. Tanner and K. Wei, “Normalized iterative hard thresholding for matrix completion,” SIAM Journal on Scientific Computing, vol. 35, no. 5, pp. S104–S125, 2013.

[130] K. Lee and Y. Bresler, “Admira: Atomic decomposition for minimum rank approximation,” IEEE Transactions on Information Theory, vol. 56, no. 9, pp. 4402–4416, 2010.

[131] J. C. Duchi and F. Ruan, “Solving (most) of a set of quadratic equalities: Composite optimization for robust phase retrieval,” arXiv preprint arXiv:1705.02356, accepted to Information and Inference: A Journal of the IMA, 2017.

[132] J. Duchi and F. Ruan, “Stochastic methods for composite optimization problems,” arXiv preprint arXiv:1703.08570, 2017.

[133] V. Charisopoulos, Y. Chen, D. Davis, M. DÃaz, L. Ding, and D. Drusvyatskiy, “Low-rank matrix recovery with composite optimization: good conditioning and rapid convergence,” arXiv preprint arXiv:1904.10020, 2019.

[134] V. Charisopoulos, D. Davis, M. Díaz, and D. Drusvyatskiy, “Composite optimization for robust blind deconvolution,” arXiv preprint arXiv:1901.01624, 2019.

[135] D. L. Donoho, M. Gavish, and A. Montanari, “The phase transition of matrix recovery from gaussian measurements matches the minimax mse of matrix denoising,” Proceedings of the National Academy of Sciences, vol. 110, no. 21, pp. 8405–8410, 2013.

[136] J. Ma, J. Xu, and A. Maleki, “Optimization-based AMP for phase retrieval: The impact of initialization and  ℓ_2-regularization,” arXiv preprint arXiv:1801.01170, 2018.

[137] ——, “Approximate message passing for amplitude based optimization,” arXiv preprint arXiv:1806.03276, 2018.

[138] W.-J. Zeng and H.-C. So, “Coordinate descent algorithms for phase retrieval,” arXiv preprint arXiv:1706.03474, 2017.

[139] S. Pinilla, J. Bacca, and H. Arguello, “Phase retrieval algorithm via nonconvex minimization using a smoothing function,” IEEE Transactions on Signal Processing, vol. 66, no. 17, pp. 4574–4584, 2018.

[140] C. Davis and W. M. Kahan, “The rotation of eigenvectors by a perturbation. iii,” SIAM Journal on Numerical Analysis, vol. 7, no. 1, pp. 1–46, 1970.

[141] ˙A. Björck and G. H. Golub, “Numerical methods for computing angles between linear subspaces,” Mathematics of Computation, vol. 27, no. 123, pp. 579–594, Jul. 1973.

[142] P.-Å. Wedin, “Perturbation bounds in connection with singular value decomposition,” BIT Numerical Mathematics, vol. 12, no. 1, pp. 99–111, 1972.

[143] T. S. Ferguson, A Course in Large Sample Theory. Chapman & Hall/CRC, 1996, vol. 38.

[144] Y. M. Lu and G. Li, “Phase transitions of spectral initialization for high-dimensional nonconvex estimation,” arXiv preprint arXiv:1702.06435, 2017.

[145] M. Mondelli and A. Montanari, “Fundamental limits of weak recovery with applications to phase retrieval,” arXiv preprint arXiv:1708.05932, 2017.

[146] W. Luo, W. Alghamdi, and Y. M. Lu, “Optimal spectral initialization for signal recovery with applications to phase retrieval,” IEEE Transactions on Signal Processing, vol. 67, no. 9, pp. 2347–2356, May 2019.

[147] K.-C. Li, “On principal Hessian directions for data visualization and dimension reduction: Another application of Stein’s lemma,” J. Am. Stat. Assoc, vol. 87, no. 420, pp. 1025–1039, 1992.

[148] D. Achlioptas and F. McSherry, “Fast computation of low-rank matrix approximations,” Journal of the ACM, vol. 54, no. 2, p. 9, 2007.

[149] K. Rohe, S. Chatterjee, and B. Yu, “Spectral clustering and the high-dimensional stochastic blockmodel,” The Annals of Statistics, vol. 39, no. 4, pp. 1878–1915, 2011.

[150] E. Abbe, “Community detection and stochastic block models: Recent developments,” Journal of Machine Learning Research, vol. 18, no. 177, pp. 1–86, 2018.

[151] S. Negahban, S. Oh, and D. Shah, “Rank centrality: Ranking from pairwise comparisons,” Operations Research, vol. 65, no. 1, pp. 266–287, 2016.

[152] Y. Chen and C. Suh, “Spectral MLE: Top-k rank aggregation from pairwise comparisons,” in International Conference on Machine Learning, 2015, pp. 371–380.

[153] A. Montanari and N. Sun, “Spectral algorithms for tensor completion,” Communications on Pure and Applied Mathematics, 2016.

[154] B. Hao, A. Zhang, and G. Cheng, “Sparse and low-rank tensor estimation via cubic sketchings,” arXiv preprint arXiv:1801.09326, 2018.

[155] A. Zhang and D. Xia, “Tensor SVD: Statistical and computational limits,” IEEE Transactions on Information Theory, 2018.

[156] C. Cai, G. Li, H. V. Poor, and Y. Chen, “Nonconvex low-rank symmetric tensor completion from noisy data,” 2019.

[157] Y. N. Dauphin, R. Pascanu, C. Gulcehre, K. Cho, S. Ganguli, and Y. Bengio, “Identifying and attacking the saddle point problem in high-dimensional non-convex optimization,” in Advances in neural information processing systems, 2014, pp. 2933–2941.

[158] A. Anandkumar and R. Ge, “Efficient approaches for escaping higher order saddle points in non-convex optimization,” in Conference on Learning Theory, 2016, pp. 81–102.

[159] R. Ge, C. Jin, and Y. Zheng, “No spurious local minima in nonconvex low rank problems: A unified geometric analysis,” in International Conference on Machine Learning, 2017, pp. 1233–1242.

[160] J. Chen and X. Li, “Memory-efficient kernel PCA via partial matrix sampling and nonconvex optimization: a model-free analysis of local minima,” arXiv preprint arXiv:1711.01742, 2017.

[161] S. Gunasekar, B. E. Woodworth, S. Bhojanapalli, B. Neyshabur, and N. Srebro, “Implicit regularization in matrix factorization,” in Advances in Neural Information Processing Systems, 2017, pp. 6151–6159.

[162] Y. Li, T. Ma, and H. Zhang, “Algorithmic regularization in over-parameterized matrix sensing and neural networks with quadratic activations,” in Conference On Learning Theory, 2018, pp. 2–47.

[163] A. Kyrillidis and A. Kalev, “Implicit regularization and solution uniqueness in over-parameterized matrix sensing,” arXiv preprint arXiv:1806.02046, 2018.

[164] M. Wang, W. Xu, and A. Tang, “A unique nonnegative solution to an underdetermined system: From vectors to matrices,” IEEE Transactions on Signal Processing, vol. 59, no. 3, pp. 1007–1016, 2011.

[165] L. Demanet and P. Hand, “Stable optimizationless recovery from phaseless linear measurements,” Journal of Fourier Analysis and Applications, vol. 20, no. 1, pp. 199–221, 2014.

[166] E. J. Candès and X. Li, “Solving quadratic equations via PhaseLift when there are about as many equations as unknowns,” Foundations of Computational Mathematics, vol. 14, no. 5, pp. 1017–1026, 2014.

[167] J. Sun, Q. Qu, and J. Wright, “Complete dictionary recovery over the sphere I: Overview and the geometric picture,” IEEE Transactions on Information Theory, vol. 63, no. 2, pp. 853–884, 2017.

[168] ——, “Complete dictionary recovery over the sphere II: Recovery by Riemannian trust-region method,” IEEE Transactions on Information Theory, vol. 63, no. 2, pp. 885–914, 2017.

[169] Y. Zhai, Z. Yang, Z. Liao, J. Wright, and Y. Ma, “Complete dictionary learning via  ℓ4-norm maximization over the orthogonal group,” arXiv preprint arXiv:1906.02435.

[170] S. Mei, Y. Bai, and A. Montanari, “The landscape of empirical risk for nonconvex losses,” The Annals of Statistics, vol. 46, no. 6A, pp. 2747–2774, 2018.

[171] R. Y. Zhang, S. Sojoudi, and J. Lavaei, “Sharp restricted isometry bounds for the inexistence of spurious local minima in nonconvex matrix recovery,” arXiv preprint arXiv:1901.01631, 2019.

[172] R. Zhang, C. Josz, S. Sojoudi, and J. Lavaei, “How much restricted isometry is needed in nonconvex matrix recovery?” in Advances in neural information processing systems, 2018, pp. 5586–5597.

[173] D. Park, A. Kyrillidis, C. Carmanis, and S. Sanghavi, “Non-square matrix sensing without spurious local minima via the burer-monteiro approach,” in Artificial Intelligence and Statistics, 2017, pp. 65–74.

[174] P. Auer, M. Herbster, and M. K. Warmuth, “Exponentially many local minima for single neurons,” in Advances in neural information processing systems, 1996, pp. 316–322.

[175] I. Safran and O. Shamir, “Spurious local minima are common in two-layer ReLU neural networks,” in International Conference on Machine Learning, 2018, pp. 4430–4438.

[176] R. Ge, F. Huang, C. Jin, and Y. Yuan, “Escaping from saddle points-online stochastic gradient for tensor decomposition.” in Conference on Learning Theory (COLT), 2015, pp. 797–842.

[177] J. Sun, “When are nonconvex optimization problems not scary?” Ph.D. dissertation, 2016.

[178] J. D. Lee, M. Simchowitz, M. I. Jordan, and B. Recht, “Gradient descent only converges to minimizers,” in Conference on Learning Theory, 2016, pp. 1246–1257.

[179] J. D. Lee, I. Panageas, G. Piliouras, M. Simchowitz, M. I. Jordan, and B. Recht, “First-order methods almost always avoid strict saddle points,” Mathematical Programming, pp. 1–27.

[180] M. Hong, M. Razaviyayn, and J. Lee, “Gradient primal-dual algorithm converges to second-order stationary solution for nonconvex distributed optimization over networks,” in International Conference on Machine Learning, 2018, pp. 2014–2023.

[181] Q. Li, Z. Zhu, and G. Tang, “Alternating minimizations converge to second-order optimal solutions,” in International Conference on Machine Learning, 2019, pp. 3935–3943.

[182] S. S. Du, C. Jin, J. D. Lee, M. I. Jordan, A. Singh, and B. Poczos, “Gradient descent can take exponential time to escape saddle points,” in Advances in Neural Information Processing Systems, 2017, pp. 1067–1077.

[183] D. Gilboa, S. Buchanan, and J. Wright, “Efficient dictionary learning with gradient descent,” ICML Workshop on Modern Trends in Nonconvex Optimization for Machine Learning, 2018.

[184] A. Brutzkus and A. Globerson, “Globally optimal gradient descent for a Convnet with gaussian inputs,” in International Conference on Machine Learning, 2017, pp. 605–614.

[185] Y. Nesterov and B. T. Polyak, “Cubic regularization of Newton method and its global performance,” Mathematical Programming, vol. 108, no. 1, pp. 177–205, 2006.

[186] A. R. Conn, N. I. Gould, and P. L. Toint, Trust region methods. SIAM, 2000, vol. 1.

[187] P.-A. Absil, C. G. Baker, and K. A. Gallivan, “Trust-region methods on Riemannian manifolds,” Foundations of Computational Mathematics, vol. 7, no. 3, pp. 303–330, 2007.

[188] Y. Carmon and J. C. Duchi, “Gradient descent efficiently finds the cubic-regularized non-convex newton step,” arXiv preprint arXiv:1612.00547, 2016.

[189] F. E. Curtis, D. P. Robinson, and M. Samadi, “A trust region algorithm with a worst-case iteration complexity of  O(ϵ−3/2)for nonconvex optimization,” Mathematical Programming, vol. 162, no. 1-2, pp. 1–32, 2017.

[190] Y. Carmon, J. C. Duchi, O. Hinder, and A. Sidford, “Accelerated methods for nonconvex optimization,” SIAM Journal on Optimization, vol. 28, no. 2, pp. 1751–1772, 2018.

[191] Y. E. Nesterov, “A method for solving the convex programming problem with convergence rate  o(1/k2),”in Dokl. Akad. Nauk SSSR, vol. 269, 1983, pp. 543–547.

[192] N. Agarwal, Z. Allen-Zhu, B. Bullins, E. Hazan, and T. Ma, “Finding approximate local minima faster than gradient descent,” in Proceedings of the 49th Annual ACM SIGACT Symposium on Theory of Computing. ACM, 2017, pp. 1195–1199.

[193] K. Y. Levy, “The power of normalization: Faster evasion of saddle points,” arXiv preprint arXiv:1611.04831, 2016.

[194] C. Jin, R. Ge, P. Netrapalli, S. M. Kakade, and M. I. Jordan, “How to escape saddle points efficiently,” in International Conference on Machine Learning, 2017, pp. 1724–1732.

[195] C. Jin, P. Netrapalli, and M. I. Jordan, “Accelerated gradient descent escapes saddle points faster than gradient descent,” in Conference On Learning Theory, 2018, pp. 1042–1085.

[196] Z. Allen-Zhu, “Natasha 2: Faster non-convex optimization than SGD,” in Advances in Neural Information Processing Systems, 2018, pp. 2675–2686.

[197] Y. Xu, J. Rong, and T. Yang, “First-order stochastic algorithms for escaping from saddle points in almost linear time,” in Advances in Neural Information Processing Systems, 2018, pp. 5530–5540.

[198] Z. Allen-Zhu and Y. Li, “Neon2: Finding local minima via first-order oracles,” in Advances in Neural Information Processing Systems, 2018, pp. 3716–3726.

[199] Z. Allen-Zhu, “Recent advances in stochastic convex and non-convex optimization,” ICML tutorial, 2017. [Online]. Available: http://people.csail.mit.edu/zeyuan/topics/icml-2017.htm

[200] N. Boumal, P.-A. Absil, and C. Cartis, “Global rates of convergence for nonconvex optimization on manifolds,” IMA Journal of Numerical Analysis, vol. 39, no. 1, pp. 1–33, 2018.

[201] Q. Qu, Y. Zhang, Y. Eldar, and J. Wright, “Convolutional phase retrieval,” in Advances in Neural Information Processing Systems, 2017, pp. 6086–6096.

[202] T. Bendory, Y. C. Eldar, and N. Boumal, “Non-convex phase retrieval from STFT measurements,” IEEE Transactions on Information Theory, vol. 64, no. 1, pp. 467–484, 2018.

[203] R. Kueng, H. Rauhut, and U. Terstiege, “Low rank matrix recovery from rank one measurements,” Applied and Computational Harmonic Analysis, vol. 42, no. 1, pp. 88–116, 2017.

[204] J. A. Tropp, “Convex recovery of a structured signal from independent random linear measurements,” in Sampling Theory, a Renaissance. Springer, 2015, pp. 67–101.

[205] Z. Yang, L. F. Yang, E. X. Fang, T. Zhao, Z. Wang, and M. Neykov, “Misspecified nonconvex statistical optimization for sparse phase retrieval,” Mathematical Programming, pp. 1–27, 2019.

[206] T. Zhang, “Phase retrieval using alternating minimization in a batch setting,” Applied and Computational Harmonic Analysis, 2019.

[207] X. Zhang, L. Wang, Y. Yu, and Q. Gu, “A primal-dual analysis of global optimality in nonconvex low-rank matrix recovery,” in International conference on machine learning, 2018, pp. 5857–5866.

[208] D. Gross, “Recovering low-rank matrices from few coefficients in any basis,” IEEE Transactions on Information Theory, vol. 57, no. 3, pp. 1548–1566, March 2011.

[209] S. Negahban and M. Wainwright, “Restricted strong convexity and weighted matrix completion: Optimal bounds with noise,” The Journal of Machine Learning Research, vol. 98888, pp. 1665–1697, May 2012.

[210] V. Koltchinskii, K. Lounici, and A. B. Tsybakov, “Nuclear-norm penalization and optimal rates for noisy low-rank matrix completion,” The Annals of Statistics, vol. 39, no. 5, pp. 2302–2329, 2011.

[211] Y. Chen, “Incoherence-optimal matrix completion,” IEEE Transactions on Information Theory, vol. 61, no. 5, pp. 2909–2923, 2015.

[212] S. Ling and T. Strohmer, “Fast blind deconvolution and blind demixing via nonconvex optimization,” in Sampling Theory and Applications (SampTA), 2017 International Conference on. IEEE, 2017, pp. 114–118.

[213] A. Ganesh, J. Wright, X. Li, E. Candès, and Y. Ma, “Dense error correction for low-rank matrices via principal component pursuit,” in International Symposium on Information Theory, 2010, pp. 1513–1517.

[214] Y. Chen, A. Jalali, S. Sanghavi, and C. Caramanis, “Low-rank matrix recovery from errors and erasures,” IEEE Transactions on Information Theory, vol. 59, no. 7, pp. 4324–4337, 2013.

[215] Y. Cherapanamjeri, K. Gupta, and P. Jain, “Nearly optimal robust matrix completion,” in International Conference on Machine Learning, 2017, pp. 797–805.

[216] Q. Gu, Z. Wang, and H. Liu, “Low-rank and sparse structure pursuit via alternating minimization,” in Artificial Intelligence and Statistics, 2016, pp. 600–609.

[217] J.-F. Cai, T. Wang, and K. Wei, “Spectral compressed sensing via projected gradient descent,” SIAM Journal on Optimization, vol. 28, no. 3, pp. 2625–2653, 2018.

[218] Y. Chen and Y. Chi, “Robust spectral compressed sensing via structured matrix completion,” IEEE Transactions on Information Theory, vol. 60, no. 10, pp. 6576–6601, 2014.

[219] G. Tang, B. N. Bhaskar, P. Shah, and B. Recht, “Compressed sensing off the grid,” IEEE Transactions on Information Theory, vol. 59, no. 11, pp. 7465–7490, 2013.

[220] V. Cambareri and L. Jacques, “A non-convex blind calibration method for randomised sensing strategies,” in 2016 4th International Workshop on Compressed Sensing Theory and its Applications to Radar, Sonar and Remote Sensing (CoSeRa). IEEE, 2016, pp. 16–20.

[221] Q. Qu, J. Sun, and J. Wright, “Finding a sparse vector in a subspace: Linear sparsity using alternating directions,” in Advances in Neural Information Processing Systems, 2014, pp. 3401–3409.

[222] X. Yi, C. Caramanis, and S. Sanghavi, “Alternating minimization for mixed linear regression,” in International Conference on Machine Learning, 2014, pp. 613–621.

[223] R. Ge and T. Ma, “On the optimization landscape of tensor decompositions,” in Advances in Neural Information Processing Systems, 2017, pp. 3653–3663.

[224] G. B. Arous, S. Mei, A. Montanari, and M. Nica, “The landscape of the spiked tensor model,” accepted to Communications on Pure and Applied Mathematics, 2017.

[225] Q. Li and G. Tang, “Convex and nonconvex geometries of symmetric tensor factorization,” in Asilomar Conference on Signals, Systems, and Computers, 2017, pp. 305–309.

[226] K. Zhong, Z. Song, P. Jain, P. L. Bartlett, and I. S. Dhillon, “Recovery guarantees for one-hidden-layer neural networks,” in International Conference on Machine Learning, 2017, pp. 4140–4149.

[227] Y. Li and Y. Yuan, “Convergence analysis of two-layer neural networks with ReLU activation,” in Advances in Neural Information Processing Systems, 2017, pp. 597–607.

[228] P. Hand and V. Voroninski, “Global guarantees for enforcing deep generative priors by empirical risk,” in Conference On Learning Theory, 2018, pp. 970–978.

[229] H. Fu, Y. Chi, and Y. Liang, “Local geometry of one-hidden-layer neural networks for logistic regression,” arXiv preprint arXiv:1802.06463, 2018.

[230] J. Fan, C. Ma, and Y. Zhong, “A selective overview of deep learning,” arXiv preprint arXiv:1904.05526, 2019.

[231] S. Burer and R. D. Monteiro, “A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization,” Mathematical Programming, vol. 95, no. 2, pp. 329–357, 2003.

[232] S. Bhojanapalli, N. Boumal, P. Jain, and P. Netrapalli, “Smoothed analysis for low-rank solutions to semidefinite programs in quadratic penalty form,” in Conference On Learning Theory, 2018, pp. 3243–3270.

[233] A. S. Bandeira, N. Boumal, and V. Voroninski, “On the low-rank approach for semidefinite programs arising in synchronization and community detection,” in Conference on Learning Theory, 2016, pp. 361–382.

[234] [Online]. Available: http://sunju.org/research/nonconvex/

[235] P. Jain and P. Kar, “Non-convex optimization for machine learning,”  Foundations and Trends R⃝ inMachine Learning, vol. 10, no. 3-4, pp. 142–336, 2017.

[236] Q. Li, Z. Zhu, G. Tang, and M. B. Wakin, “The geometry of equality-constrained global consensus problems,” in ICASSP, May 2019, pp. 7928–7932.

Designed for Accessibility and to further Open Science