Nonconvex Optimization Meets Low-Rank Matrix Factorization: An Overview

2018·arXiv

Abstract

1 Introduction

Modern information processing and machine learning often have to deal with (structured) low-rank matrix factorization. Given a few observations about a matrix seeks a low-rank solution compatible with this set of observations as well as other prior constraints. Examples include low-rank matrix completion [1–3], 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 as with low-rank factors and , one attempts recovery by solving an optimization problem in the form of

Here, 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 [1–3, 11–13], 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. [16–26]). 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 [16–22].

• 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 [27–31].

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 denote its respectively. For any matrix 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 norm (i.e. the largest norm of the rows), and the entrywise norm (the largest magnitude of all entries), respectively. We denote by th largest singular value (resp. eigenvalue) of M, and let (resp. ) represent its jth row (resp. column). The condition number of M is denoted by . In addition, indicate 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 stands for the trace. The matrix denotes the identity matrix, and denotes the ith column of . For a linear operator A, denote by its adjoint operator. For example, if . We also let vec(Z) denote the vectorization of a matrix Z. The indicator function equals 1 when the event A holds true, and 0 otherwise. Further, the notation denotes the set of orthonormal matrices. Let A and B be two square matrices. We write ) if their difference is a positive definite (resp. positive semidefinite) matrix.

Additionally, the standard notation means that there exists a constant c > 0 such that means that there exists a constant c > 0 such that , and means that there exist constants

2 Preliminaries in optimization theory

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

The optimal solution, if it exists, is denoted by

When f(x) is strictly convex,is unique. But it may be non-unique when 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, is the step size or learning rate at the tth iteration, and is the initial point. This method and its variants are widely used in practice, partly due to their simplicity and scalability to large-scale problems.

A central question is when GD converges fast to the global minimum . As is well-known in the optimization literature, GD is provably convergent at a linear rate when is (locally) strongly convex and smooth. Here, an algorithm is said to converge linearly if the error converges 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 is said to be -strongly convex in a set B if

Definition 2 (Smoothness). A twice continuously differentiable function is said to be in a set B if

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

Lemma 1. Suppose that -strongly convex and -smooth within a local ball and that , then GD obeys

Proof of Lemma 1. The optimality of indicates that , which allows us to rewrite the GD update rule as

where . Here, the second line arises from the fundamental theorem of calculus [33, Chapter XIII, Theorem 4.2]. If , then it is self-evident that , which combined with the assumption of Lemma 1 gives

Therefore, as long as (and hence

This together with the sub-multiplicativity of

By setting , we arrive at the desired error contraction, namely,

A byproduct is: if , then the next iterate also falls in . Consequently, applying the above argument recursively and recalling the assumption , we see that all GD iterates remain within ) holds true for all t. This immediately concludes the proof.

Figure 1: An example of taken from [20]. Here, . It satisfies , but is clearly nonconvex.

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

if we initialize GD properly such that lies in the local region . 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 can be reasonably large (e.g. on the same order of

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 . In vanilla GD, , but can also incorporate several variants of GD; see Section 6. The regularity condition is defined as follows.

is said to obey the regularity condition

for all

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

It is worth noting that this condition does not require to be differentiable. Also, when it does not require to be convex within ; 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 ) defined as follows

for all . 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 Under the regularity condition, the iterative algorithm (8) converges linearly with suitable initialization.

, the iterative algorithm (8) with

Proof of Lemma 2. Assuming that , we can obtain

where the first inequality comes from , and the last line arises if . By setting arrive at

which also shows that . The claim then follows by induction under the assumption that

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) is any point that satisfies

Moreover, we call a point , if it satisfies

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

• if is a local maximum;

• if is a local minimum;

Figure 2: Illustration of the function

• is either a local minimum or a degenerate saddle point;

• saddle 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

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.

3 A warm-up example: rank-1 matrix factorization

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

Specifically, consider a positive semidefinite matrix (which is not necessarily low-rank) with eigenvalue decomposition . 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:

where is a degree-four polynomial and highly nonconvex. The solution to (14) can be expressed in closed form as the scaled leading eigenvector . See Fig. 2 for an illustration of the function f(x) when . 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 [36–42].

We attempt to minimize the nonconvex function 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. ), is guaranteed to converge fast.

Theorem 1. Consider the problem (14). Suppose and . Then the GD iterates (4) obey

Remark 1. By symmetry, Theorem 1 continues to hold if is replaced by

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 . Its radius, which is not optimized in this theorem, is given by . This also depends on the relative eigen-gap

Proof of Theorem 1. The proof mainly consists of showing that 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 . A little algebra yields that if

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

We can then categorize them as follows:

1. With regards to the points

2. For the points

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.

4 Formulations of a few canonical problems

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

4.1 Matrix sensing

Suppose that we are given a set of m measurements of of the form

where is a collection of sensing matrices known a priori. We are asked to recover — which is assumed to be of rank r — from these linear matrix equations [12, 43].

When is positive semidefinite with , this can be cast as solving the least-squares problem

Clearly, we cannot distinguish for any orthonormal matrix , as they correspond to the same low-rank matrix . 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 from for any invertible matrix , as . Throughout the paper, we denote by

the true low-rank factors, where 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 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 is said to satisfy the r-RIP with RIP constant

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

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 are i.i.d. Gaussian entries N(0, 1), then A as defined in (23) satisfies the r-RIP with RIP constant with high probability as soon as

4.2 Phase retrieval and quadratic sensing

Imagine that we have access to m quadratic measurements of a rank-1 matrix

where is the design vector known a priori. How can we reconstruct — or equivalently, — from this collection of quadratic equations about ? 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, 45–47]. 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-

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 . 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 of rank r, indexed by the sampling location set . It is convenient to introduce a projection operator such that for an input matrix

The matrix completion problem then boils down to recovering (or equivalently, from the partially observed entries of ) [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

For the positive semidefinite case where , the task can be cast as solving

When it comes to the more general case where , 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].

with singular value decomposition (SVD) 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 when referring to the matrix completion problem, and set to be the condition number of

4.4 Blind deconvolution (the subspace model)

Suppose that we want to recover two objects and — or equivalently, the outer product bilinear measurements of the form

To explain why this is called blind deconvolution, imagine we would like to recover two signals and from their circulant convolution [5]. In the frequency domain, the outputs can be written as , where (resp. ) is the Fourier transform of g (resp. d). If we have additional knowledge that and lie in some known subspace characterized by and B = reduces 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 and that is formed by the first K columns of a unitary discrete Fourier transform (DFT) matrix

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 be the smallest number such that

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

Suppose we are given a matrix that is a superposition of a rank-r matrix and a sparse matrix

The goal is to separate 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 when the observed entries are corrupted by sparse outliers (modeled by ). 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 as defined in Definition 7 as well, which precludes from being too spiky. In addition, it is sometimes convenient to introduce the following deterministic condition on the sparsity pattern and sparsity level of originally proposed by [7]. Specifically, it is assumed that the non-zero entries of 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

For the positive semidefinite case where , the task can be cast as solving

The general case where can be formulated similarly by replacing optimizing over

5 Local reﬁnement via gradient descent

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

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

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. for some vector . This requires us to solve

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

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

When the sensing matrices are random and isotropic, (42) can be viewed as a randomized version of the rank-1 matrix factorization problem discussed in Section 3. To see this, consider for instance the case where the ’s are drawn from the symmetric Gaussian design, i.e. the diagonal entries of are 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 . This bodes well for fast local convergence of GD, at least at the population level (i.e. the case when the sample size

What happens in the finite-sample regime? It turns out that if the sensing operator satisfies the RIP (cf. Definition 6), then does not deviate too much from its population-level counterpart, and hence 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 , then GD with

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 casewith , in which the gradient update rule (4) satisfies

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

is identity, then reduces to

It can be shown that for any

Think of the following example

for two unit vectors . It is straightforward to verify that

which violates convexity. Moreover, this happens even when X is arbitrarily close to (by taking

Fortunately, the above issue can be easily addressed. The key is to recognize that: one can only hope to recover 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 is 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 is a global minimizer of , and make the further homogeneity assumption for any orthonormal matrix — a common fact that arises in matrix factorization problems.

Lemma 4. Suppose that -smooth within a ball for any orthonormal matrix . Assume that for any

If , then GD with

Proof of Lemma 4. For notational simplicity, let

First, by definition of

Next, we claim that a modified regularity condition (cf. Definition 3) holds for all

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 for some . If lies within , then one has 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 . Finally, since the initial point is assumed to fall within immediately learn from (53) and induction that . 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 , 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

which becomes strictly positive for . In fact, if X is sufficiently close to , 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 is locally not far from in Example 1, meaning that the condition (51) continues to hold for some . This leads to the following result. It is assumed that the ground truth 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 . Then there exist some universal constants such that if , then GD with

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 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. as suggested in [23, 55].Here, the regularization term 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

(resp.

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

Lemma 5 (Local strong convexity and smoothness for phase retrieval). Consider the problem (27). There exist some constants such that with probability at least

holds simultaneously for all , provided that

This lemma says that is locally 0.5-strongly convex and -smooth when the sample size 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 , provided that

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 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: 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 incur ill-conditioned Hessians. For instance, suppose is a unit vector independent of . Then the point for some constant often results in extremely large This simple instance suggests that: in order to ensure well-conditioned Hessians, one needs to preclude points that are too “coherent” with the sampling vectors, as formalized below.

Lemma 6 (Restricted smoothness for phase retrieval [26]). Under the assumptions of Lemma 5, there exist some constants such that if , then with probability at least

holds simultaneously for all

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 “the sampling vectors.

Going beyond phase retrieval, the notion of incoherence is not only crucial to control smoothness, but also plays a critical role in ensuring local strong convexity (or regularity conditions). A partial list of examples include matrix completion, quadratic sensing, blind deconvolution and demixing, etc [19, 21, 22, 26, 60–62]. 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 for some large constant . Then with probability , the Hessian obeys

for all satisfying

where for some constant

This lemma confines attention to the set of points obeying

Given that each observed entry can be viewed as , 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, 19–22, 61, 63, 64], in the hope of improving computational guarantees. Specifically, one can regularize the loss function, by adding additional regularization term to the objective function and designing the GD update rule w.r.t. the regularized problem

with the regularization parameter. For example:

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

for some scalars . There are numerous choices of ; the one suggested by [19] is . With suitable learning rates and proper initialization, GD w.r.t. the regularized loss provably yields -accuracy in iterations, provided that the sample size

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

with . It has been demonstrated that under proper initialization and step size, gradient methods w.r.t. the regularized loss reach -accuracy in iterations, provided that the sample size

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 error 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 . The second line arises since is -strongly convex, and the last line uses the smoothness condition.

where . The second line arises since is -strongly convex, and the last line uses the smoothness condition.

where . The second line arises since 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 ball, then standard GD theory implies convergence. (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 obey

for all with probability are some constants, and we assume

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 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 for some large constant , then with probability , the GD iterates with proper initialization (see, e.g., spectral initialization in Section 8.2) and obey

for all is some absolute constant.

This theorem demonstrates that vanilla GD converges within 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 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 for some large constant , and that the condition number is a fixed constant. With probability at least , the GD iterates (71) with proper initialization (see, e.g., spectral initialization in Section 8.2) satisfy

for all , with defined in (50). Here, is some constant, and for some constant

This theorem demonstrates that vanilla GD converges within 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

The last example is blind deconvolution. To measure the discrepancy between any 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 for some large constant . Then there is some constant such that with probability exceeding , the iterates (74) with proper initialization (see, e.g., spectral initialization in Section 8.2) and

For conciseness, we only state that the estimation error converges in 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 [67–78]. 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].

6 Variants of gradient descent

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

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. ), one runs projected GD w.r.t. the loss function

where is the step size and denote the Euclidean projection onto the set of incoherent matrices:

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

for i = 1, 2, . . . , 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 for some large constant , and that the condition number of is a fixed constant. With probability at least , the projected GD iterates (77) satisfy

for all , provided that and for some constant

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 . 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 in (26) is k-sparse, where . If we have prior knowledge about , then we can pick the constraint set C as follows to promote sparsity

as a sparse signal often (although not always) has low norm. 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 is k-sparse. Suppose that for some large constant . The projected GD iterates w.r.t. (27) and the constraint set (80) obey

with probability at least , provided that

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, 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 when (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 under i.i.d. Gaussian designs, unless k is sufficiently large or other structural information is available [48, 80, 83–85]. 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 . When it comes to the sample-limited regime where , 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 , the gradient — which depends on 4th moments of 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 whose size deviate too much from the typical size. Specifically, the truncation rule proposed in [20] is:

for some trimming criteria defined as

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

Theorem 12 (Truncated GD for phase retrieval [20]). Consider the problem (27). With probability the iterates (84) obey

for some constant , provided that and for some constants

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 . 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 is orderwise optimal;

• Optimal computational complexity: truncated WF yields accuracy in 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 by hard thresholding the residual, namely,

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

where ) 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 then becomes our estimate of the sparse outlier matrix 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 is a fixed constant. Suppose that the sample size and the sparsity of the outlier satisfy and for some constants . With probability at least , the iterates satisfy

for all , provided that . Here, are some constants, and for some constant

Remark 6. In the full data case, the convergence rate can be improved to for some constant

This theorem essentially says that: as long as the fraction of entries corrupted by outliers does not exceed , then the nonconvex algorithm described above provably recovers the true low-rank matrix in about 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 . Mathematically, the measurement model in the presence of outliers is given by

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

Here, only includes samples whose residual size does not deviate much from the

where 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 and such that if , , and , then with probability at least , the median-truncated GD iterates satisfy

6.3 Generalized gradient descent

In all the examples discussed so far, the loss function has been a smooth function. When 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, 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 often has better curvature around the truth, compared to the intensity-based loss function ] 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 such that if , then with high probability,

as long as

In comparison to Theorem 12, the generalized GD w.r.t. 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. 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 , the truncation step is not crucial and can be safely removed when dealing with the amplitude-based . A main reason is that for any fixed x, involves only the first and second moments of the (sub)-Gaussian random variables it exhibits much sharper measure concentration — and hence much better controlled gradient components — compared to the heavy-tailed , which involves fourth moments of . 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 given their pairwise relative phases. Alternatively, by setting , this problem reduces to estimating from — a matrix that encodes all pairwise phase differences . To account for the noisy nature of practical measurements, suppose that what we observe is , where W is a Hermitian matrix. Here, are 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 , where each variable can take m possible values, namely, . Suppose that estimation needs to be performed based on pairwise difference samples , where the ’s are i.i.d. noise and their distributions dictate the recovery limits. To facilitate computation, one strategy is to lift each discrete variable into a m-dimensional vector . 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

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 . 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 , we might assume , where denotes the Grassmann manifold which parametrizes all r-dimensional linear subspaces of the -dimensional space. In words, we are searching for a r-dimensional subspace L but ignores the global rotation. It is also assumed that to remove the global scaling ambiguity (otherwise is always equivalent to ). One might then try to minimize the loss function defined over the Grassmann manifold as follows

where

As it turns out, it is possible to apply GD to 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 w.r.t. the [100]:

where is the least-squares solution. The gradient on the Grassmann manifold, denoted by , is then given by

Let be its compact SVD, then the geodesic on the Grassmann manifold along the direction is given by

We can then update the iterates as

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

with . As can be verified, automatically stays on the unit sphere obeying 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 , 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 , provided that a proper initialization is given.

Other gradient descent approaches on manifolds include [100–106]. 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) [108–110], 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 is a subset of cardinality k selected uniformly at random. Here, is known as the mini-batch size. As one can expect, the mini-batch size k plays an important role in the trade-off between the computational cost per iteration and the convergence rate. A properly chosen mini-batch size will optimize the total computational cost given the practical constraints. Please see [90, 111–115] 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.

7 Beyond gradient methods

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

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 and are updated sequentially. Here, is 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 and are full-rank, the column space of matches perfectly with that of . Armed with this fact, the subsequent least squares problem (i.e. the update for ) is exact, in the sense that

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 . If we initialize leading left singular vectors of , 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 iterations, compared to 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 . Intuitively, if there is a singular value that is much larger than the remaining ones, then one can treat 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 and the signal

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 is an appropriate initial estimate, is the pseudo-inverse of , and 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 and such that if , then with probability at least , the estimates of AltMin (ER) satisfy

as long as

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 , AltMin proceeds as follows: for t = 1, 2, . . .

where 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, 122–125]:

where denotes 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 satisfiesafter iterations, provided that the sample complexity exceeds . 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 , with the proviso that the sample complexity exceeds . 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

where

Here, is the step size, and returns the best rank-r approximation of Z. Given that the iterates are always low-rank, one can store in 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 . If we initialize and adopt a step size then SVP achieves

Theorem 19 (SVP for matrix completion [76]). Consider the problem (32) and set . Suppose that for some large constant , and that the step size is set as . Then with probability exceeding , the SVP iterates achieve

as long as for some constant

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) [131–134], approximate message passing [135–137], block coordinate descent [19], coordinate descent [138], and conjugate gradient [139]. The readers are referred to these papers for detailed descriptions.

8 Initialization via spectral methods

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 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 be a symmetric matrix, whose eigenvalues are real-valued. In many cases, we only have access to a perturbed version ), 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 , with the proviso that the perturbation is sufficiently small in size. This was first established in the celebrated Davis-Kahan Theorem [140]. Specifically, let the eigenvalues of be partitioned into two groups

where . We assume that the eigen-gap between the two groups, , is strictly positive. For example, if and 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 (resp. ) an orthonormal matrix whose columns are the first r eigenvectors of ). In order to measure the distance between the two subspaces spanned by , 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

These results are referred to as the theorem because the distance metric is identical to are the so-called principal angles [141] between the two subspaces spanned by , respectively.

Furthermore, to deal with asymmetric matrices, similar perturbation bounds can be obtained. Suppose that ) consist respectively of the first r left and right singular vectors of (resp. Y ). Then we have the celebrated Wedin 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 denotes its adjoint. The generic version of the spectral method then proceeds by computing (i) two matrices and whose columns consist of the top-r left and right singular vectors of Y , respectively, and (ii) a diagonal matrix that contains the corresponding top-r singular values. In the hope that U, V and are reasonably reliable estimates of , respectively, we take

as estimates of the low-rank factors is known to be positive semidefinite, then we can also let be 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 can be well-controlled. A simple way to control arises when A satisfies the RIP in Definition 6.

Lemma 8. Suppose that 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

from which (123) follows.

In what follows, we first illustrate how to control the estimation error in the rank-1 case where . In this case, the leading eigenvector u of Y obeys

where (i) comes from Theorem 20, and (ii) holds if (which is guaranteed if according to Lemma 8). Similarly, we can invoke Weyl’s inequality and Lemma 8 to control the gap between the leading eigenvalue

where the last inequality makes use of (124), (125), and (??). This characterizes the difference between our estimateand the true low-rank factor

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 . Suppose A satisfies 2r-RIP with RIP constant for some sufficiently small constant . Then the spectral estimate (122) obeys

Remark 9. It is worth pointing out that in view of Fact 1, the vanilla spectral method needs to 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 . 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 with sensing matrices . To obtain an initial guess that is close to the truth , 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 as

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 is a reasonably good estimate of . 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. 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 is an approximated version of E[Y ], we expect ) to carry useful information about

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 theorem, one arrives at the following estimates:

Theorem 23 (Spectral method for phase retrieval [18]). Consider phase retrieval in Section 4.2, where is any given vector. Fix any , and suppose for some sufficiently large constant Then the spectral estimate (127) obeys

with probability at least

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 , where serves as an estimate of . In words, the diagonal entries of can be approximately viewed as the top-r eigenvalues of . The initial guess is then set as

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

Theorem 24 (Spectral method for quadratic sensing [18]). Consider quadratic sensing in Section 4.2, where . Fix any , and suppose for some sufficiently large constant . Then the spectral estimate (130) obeys

with probability at least

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 with sensing matrices . Following (121), we form a surrogate matrix

Let 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 satisfies the incoherence condition in Definition 8 with parameter , and assume For any for some sufficiently large constant , then the spectral estimate (131) obeys

with probability at least

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 . Following the aforementioned procedure, we can form a surrogate matrix as

Notably, the scaling factor in (132) is chosen to ensure that . We then construct the initial guess for the low-rank factors in the same manner as (122), using Y in (133).

Asis a collection of independent random matrices, we can use the matrix Bernstein inequality [57] to get a high-probability upper bound on the deviation . 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 , and suppose the condition number is a fixed constant. There exist a constant that if , then with probability at least , 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 themRecall that Theorem 23 requires a sample complexity , which is a logarithmic factor larger than the signal dimension n. What happens if we only have access to samples, 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

Under the Gaussian design,is a collection of i.i.d. random variables with 1 degree of freedom. It follows from well-known estimates in extreme value theory [143] that

Recall that 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 . This explains why we need a sample complexity in Theorem 23.

The above analysis also suggests an easy fix: since the main culprit lies in the fact that is unbounded (as ), we can apply a preprocessing function to to 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

where

for some predetermined truncation threshold . The initial point is then formed by scaling the leading eigenvector of to have roughly the same norm of , which can be estimated by essentially performing a trimming operation, removing any entry of that 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 , and suppose for some sufficiently large constant . Then the truncated spectral estimate obeys

with probability at least

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

where -largest value (e.g. . In words, this method only employs a subset of design vectors that are better aligned with the truth . 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 . 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 . There exist some constants such that if , then the median-truncated spectral estimate obeys

with probability at least

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 before 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 defined in (86), as

where c > 0 is some predetermined constant (e.g. 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 is a fixed constant. Fix . If the sample size and the sparsity fraction satisfy and for some large constant , then with probability at least

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 is k-sparse (see Section 6.1.2). A simple idea is to first identify the support of , 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 , 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 against some preset threshold

The nonzero part of is then found by applying the spectral method outlined in Section 8.2.2 to

, where is the subvector of coming from the support . In short, this strategy provably leads to a reasonably good initial estimate, as long as ; 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 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 and the initialization . The result is this:

Theorem 30 (Precise asymptotic characterization [144]). Consider phase retrieval in Section 4.2, and let for some constant . Under mild technical conditions on the preprocessing function T , the leading eigenvector

as . Here, is a fixed constant and is a fixed function that is positive when . Furthermore, converges to a positive constant iff

Remark 10. The above characterization was first obtained in [144], under the assumption that is nonnegative. Later, this technical restriction was removed in [145]. Analytical formulas of available for any given , 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 that marks the transition between two very contrasting regimes.

• An uncorrelated phase takes place when the sampling ratio . Within this phase, meaning that the spectral estimate is uncorrelated with the target.

• A correlated phase takes place when . 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 , which in turn implies that can be efficiently computed by the power method.

Figure 4: Comparisons of 3 designs of the preprocessing function for estimating a image 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 defined in (144); (Blue curve) the uniformly optimal design

The phase transition boundary is determined by the preprocessing function . 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 for any . Further, can be approached by the following preprocessing function

The value is called the , no algorithm can generate an estimate that is asymptotically positively correlated with . The function 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 , we seek a preprocessing function that maximizes the squared cosine similarity, namely,

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

Theorem 32 (Uniformly optimal preprocessing). Under the same settings of Theorem 31, we have 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 image 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 [153–156]. We have to omit these due to the space limit.

9 Global landscape and initialization-free algorithms

A separate line of work aims to study the global geometry of a loss function 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 are either global minima or strict saddle points, we say that 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.

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.

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

• (All local minima are global): for any local minimum , it satisfies

• (Strict saddles): for any critical point X that is not a local minimum, it satisfies

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

Consider any local minimum . This is characterized by the first-order and second-order optimality conditions

Without loss of generality, assume that

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. . As a result, the proof consists in showing that: when the RIP constant is sufficiently small, one has

unless . Additionally, the value of (147) is helpful in upper bounding 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 when accounting for global orthogonal transformation (Section 5.1.2). Another way to express this is through a different parameterization

which clearly satisfies . It is not hard to show that: in the presence of the RIP, the Hessian 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 . 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])

admits a solution with . Assume that for all M and D with and

holds with . Then the function , has no spurious local minima, and all saddle points of are strict saddles.

Remark 11. This result continue to hold if (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 , then with high probability, there is no spurious local minimum, and all saddle points of 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 with a regularized loss

with a regularization parameter, and . For properly selected and , if the sample size , then with high probability, all local minima X of , and all saddle points of 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 is the true low-rank matrix. Here, for phase retrieval. As it turns out, under minimal sample complexity, 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 of the above function , provided that the set

is a singleton

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

• matrix sensing (the positive semidefinite case where ), as long as A satisfies 5r-RIP with a RIP constant ]; a necessary and sufficient condition was established in [164];

• phase retrieval, as long as for some sufficiently large constant

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 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 . Consider any local minimum of . By definition, there is an with small enough

Now suppose that ) is not a global solution of ), then there exists a (resp. ) arbitrarily close to ) such that

For instance, one can take . This implies the existence of an that 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 [167–169]. 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 loss exhibits benign nonconvex geometry over a sphere.

• M-estimator in statistical estimation [170]. Given a set of independent data points , 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 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 , 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 [159, 173] or by the term [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 is said to satisfy the -strict saddle property for some , 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 ), 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 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

Now suppose that be the subspace spanned by the first eigenvectors. It is easy to see that 0 is a saddle point, and that only if . In fact, as long as contains a component outside , then this component will keep blowing up at a rate that -dimensional, we have . As a result,

The above theory has been extended to accommodate other optimization methods like coordinate descent, mirror descent, the gradient primal-dual algorithm, and alternating minimization [179–181]. 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 . The GD iterates with random initialization

with probability are some constants, , and we assume

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

• : 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 in this stage grows exponentially fast, while the residual strengthdoes not increase by much.

• : once the iterates enters the local basin, the estimation error decays exponentially fast, similar to Theorem 6. This stage takes about iterations to reach

Taken collectively, GD with random initialization achieves -accuracy in 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 ). 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. ). 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 -Hessian-Lipschitz if, for any has

9.3.1 Hessian-based algorithms

This class of algorithms requires an oracle that returns 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 happens 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, is 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 . 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 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 for 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 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 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 is some noise sampled uniformly from a sphere, and 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 increases. 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 uniformly 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 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 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 [196–198]. 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, for various large-scale problems.

10 Concluding remarks

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, [20–22, 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, 222–225]; (3) parameter recovery and inverse problems in shallow and deep neural networks [184, 226–230]; (4) analysis of Burer-Monteiro factorization to semidefinite programs [55, 231–233]. 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.

Acknowledgment

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.

A Proof of Theorem 3

To begin with, simple calculation reveals that for any , 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 To establish local strong convexity and smoothness, we need to control for all z. The key ingredient is to show that: is 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 . Similarly,

As a result, if is small enough, then putting the above results together implies that: is 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.

B Proof of Claim (164)

C Modiﬁed strong convexity for (46)

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

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

for some universal constant , provided that is small enough. We then turn attention to

The last line holds by recognizing that , Theorem 2], which implies that

as long as . In summary,

D Proof of Theorem 22

We start with the rank-r PSD case, where we denote by the initial estimate (i.e. case) and the ground truth. Observe that

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

As an immediate consequence of the above lemma, we get

where the last line follows from (165). Therefore, as long as

In fact, [127] shows a stronger consequence of the RIP: one can improve the left-hand-side of (165) to , which allows relaxing the requirement on the RIP constant to asymmetric case can be proved in a similar fashion.

E Proof of Theorem 33 (the rank-1 case)

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 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 . This concludes the proof.

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

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

References

[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,” Learning, 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,” Machine 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 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, 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-, 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 -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 -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 -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 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 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,” Machine 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