Global Optimality in Low-rank Matrix Optimization

Zhihui Zhu, Qiuwei Li, Gongguo Tang, Michael B. Wakin

I Introduction

Consider the minimization of a general objective function f(X)f(\boldsymbol{X}) over all low-rank n×mn\times m matrices:

The bilinear nature of the parameterization renders the objective function of (2) nonconvex even when f(X)f(\boldsymbol{X}) is a convex function. Hence, the objective function in (2) can potentially have spurious local minima (i.e., local minimizers that are not global minimizers) or “bad” saddle points that prevent a number of iterative algorithms from converging to the global solution. By analyzing the landscape of nonconvex functions, several recent works have shown that the factored objective function h(U,V)h(\boldsymbol{U},\boldsymbol{V}) in certain matrix inverse problems has no spurious local minima .

We generalize this line of work by focusing on a general objective function f(X)f(\boldsymbol{X}) in the optimization (1), not necessarily a quadratic loss function coming from a matrix inverse problem. By focusing on a general objective function, we attempt to provide a unifying framework for low-rank matrix optimizations with the factorization approach. We provide a geometric analysis for the factored program (2) and show that, under certain conditions on f(X)f(\boldsymbol{X}), all critical points of the objective function h(U,V)h(\boldsymbol{U},\boldsymbol{V}) are well-behaved. Our characterization of the geometry of the objective function ensures that a number of iterative optimization algorithms converge to a global minimum.

The purpose of this paper is to analyze the geometry of the factored problem h(U,V)h(\boldsymbol{U},\boldsymbol{V}) in (2). In particular, we attempt to understand the behavior of all of the critical points of the objective function in the reformulated problem (2).

Before presenting our main results, we lay out the necessary assumptions on the objective function f(X)f(\boldsymbol{X}). As is known, without any assumptions on the problem, even minimizing traditional quadratic objective functions is challenging. For this purpose, we focus on the model where f(X)f(\boldsymbol{X}) is (2r,4r)(2r,4r)-restricted strongly convex and smooth, i.e., for any n×mn\times m matrices X,G\boldsymbol{X},\boldsymbol{G} with rank⁡(X)≤2r\operatorname{rank}(\boldsymbol{X})\leq 2r and rank⁡(G)≤4r\operatorname{rank}(\boldsymbol{G})\leq 4r, the Hessian of f(X)f(\boldsymbol{X}) satisfies

for some positive α\alpha and β\beta. A similar assumption is also utilized in [20, Conditions 5.3 and 5.4]. With this assumption on f(X)f(\boldsymbol{X}), we summarize our main results in the following informal theorem.

As guaranteed by Proposition 1 (in Section III), the (2r,4r)(2r,4r)-restricted strong convexity and smoothness property (3) ensures that X⋆\boldsymbol{X}^{\star} is the unique global minimum of (1). Theorem 1 then implies that we can recover the rank-r⋆r^{\star} global minimizer X⋆\boldsymbol{X}^{\star} of (1) by many iterative algorithms (such as the trust region method and stochastic gradient descent ) even from a random initialization. This is because 1) as guaranteed by Theorem 2, the strict saddle property ensures local search algorithms converge to a local minimum, and 2) there are no spurious local minima.

Since our main result only requires the (2r,4r)(2r,4r)-restricted strong convexity and smoothness property (3), aside from low-rank matrix recovery , it can also be applied to many other low-rank matrix optimization problems which do not necessarily involve quadratic loss functions. Typical examples include robust PCA , 1-bit matrix completion and Poisson principal component analysis (PCA) .

I-B Related Works

Compared with the original program (1), the factored form (2) typically involves many fewer variables (or variables with much smaller size) and can be efficiently solved by simple but powerful methods (such as gradient descent , the trust region method , and alternating methods ) for large-scale settings, though it is nonconvex. In recent years, tremendous effort has been devoted to analyzing nonconvex optimizations by exploiting the geometry of the corresponding objective functions. These works can be separated into two types based on whether the geometry is analysed locally or globally. One type of work analyzes the behavior of the objective function in a small neighborhood containing the global optimum and requires a good initialization that is close enough to a global minimum. Problems such as phase retrieval , matrix sensing , and semi-definite optimization have been studied.

Another type of work attempts to analyze the landscape of the objective function and show that it obeys the strict saddle property. If this particular property holds, then simple algorithms such as gradient descent and the trust region method are guaranteed to converge to a local minimum from a random initialization rather than requiring a good guess. We approach low-rank matrix optimization with general objective functions (1) via a similar geometric characterization. Similar geometric results are known for a number of problems including complete dictionary learning , phase retrieval , orthogonal tensor decomposition , and matrix inverse problems . Empirical evidence also supports using the factorization approach for estimating a low-rank PSD matrix from a set of rank-one measurements corrupted by arbitrary outliers and for recovering a dynamically evolving low-rank matrix from incomplete observations .

Our work is most closely related to certain recent works in low-rank matrix optimization. Bhojanapalli et al. showed that the low-rank, PSD matrix sensing problem has no spurious local minima and obeys the strict saddle property. Similar results were exploited for PSD matrix completion , PSD matrix factorization and low-rank, PSD matrix optimization problems with generic objective functions . Our work extends this line of analysis to general low-rank matrix (not necessary PSD or even square) optimization problems. Another closely related work considers the low-rank, non-square matrix sensing problem and matrix completion with the factorization approach . We note that our general objective function framework includes the low-rank matrix sensing problem as a special case (see Section III-C). Furthermore, our result covers both over-parameterization where r>r⋆r>r^{\star} and exact parameterization where r=r⋆r=r^{\star}. Wang et al. also considered the factored low-rank matrix minimization problem with a general objective function which satisfies the restricted strong convexity and smoothness condition. Their algorithms require good initializations for global convergence since they characterized only the local landscapes around the global optima. By categorizing the behavior of all the critical points, our work differs from in that we instead characterize the global landscape of the factored objective function.

This paper continues in Section II with formal definitions for strict saddles and the strict saddle property. We present the main results and their implications in matrix sensing, weighted low-rank approximation, and 1-bit matrix completion in Section III. The proof of our main results is given in Section IV. We conclude the paper in Section VI.

II Preliminaries

II-B Strict Saddle Property

We say x\boldsymbol{x} a critical point if the gradient at x\boldsymbol{x} vanishes, i.e., ∇h(x)=0\nabla h(\boldsymbol{x})={\bf 0}.

A critical point x\boldsymbol{x} is a strict saddle if the Hessian matrix evaluated at this point has a strictly negative eigenvalue, i.e., λmin⁡(∇2h(x))<0\lambda_{\min}(\nabla^{2}h(\boldsymbol{x}))<0.

A twice differentiable function satisfies the strict saddle property if each critical point either corresponds to a local minimum or is a strict saddle.

Intuitively, the strict saddle property requires a function to have a directional negative curvature at all critical points but local minima. This property allows a number of iterative algorithms such as noisy gradient descent and the trust region method to further decrease the function value at all the strict saddles and thus converge to a local minimum.

(informal) For a twice continuously differentiable objective function satisfying the strict saddle property, a number of iterative optimization algorithms (such as gradient descent and the the trust region method) can find a local minimum.

III Problem Formulation and Main Results

This paper considers the problem (1) of minimizing a general function f(X)f(\boldsymbol{X}) (over the set of low-rank matrices) which is assumed to have a low-rank critical point X⋆\boldsymbol{X}^{\star} with rank⁡(X⋆)=r⋆≤r\operatorname{rank}(\boldsymbol{X}^{\star})=r^{\star}\leq r such that ∇f(X⋆)=0\nabla f(\boldsymbol{X}^{\star})={\bf 0}. Because of the restricted strong convexity and smoothness condition (3), the following result establishes that if f(X)f(\boldsymbol{X}) has a critical point X⋆\boldsymbol{X}^{\star} with rank⁡(X⋆)≤r\operatorname{rank}(\boldsymbol{X}^{\star})\leq r, then it is the unique global minimum of (1).

Suppose f(X)f(\boldsymbol{X}) satisfies the (2r,4r)(2r,4r)-restricted strong convexity and smoothness condition (3) with positive α\alpha and β\beta. Assume X⋆\boldsymbol{X}^{\star} is a critical point of f(X)f(\boldsymbol{X}) with rank⁡(X⋆)=r⋆≤r\operatorname{rank}(\boldsymbol{X}^{\star})=r^{\star}\leq r. Then X⋆\boldsymbol{X}^{\star} is the global minimum of (1), i.e.,

and the equality holds only at X=X⋆\boldsymbol{X}=\boldsymbol{X}^{\star}.

First note that if X⋆\boldsymbol{X}^{\star} is a critical point of f(X)f(\boldsymbol{X}), then

where X~=tX⋆+(1−t)X\widetilde{\boldsymbol{X}}=t\boldsymbol{X}^{\star}+(1-t)\boldsymbol{X} for some t∈t\in. This Taylor expansion together with ∇f(X⋆)=0\nabla f(\boldsymbol{X}^{\star})={\bf 0} and (3) (both X~\widetilde{\boldsymbol{X}} and X′−X⋆\boldsymbol{X}^{\prime}-\boldsymbol{X}^{\star} have rank at most 2r2r) gives

Although the new variable W\boldsymbol{W} has much smaller size than X\boldsymbol{X} when r≪min⁡{n,m}r\ll\min\{n,m\}, the objective function in the factored problem (2) may have a much more complicated landscape due to the bilinear form about U\boldsymbol{U} and V\boldsymbol{V}. The reformulated objective function h(U,V)h(\boldsymbol{U},\boldsymbol{V}) could introduce spurious local minima or degenerate saddle points even when f(X)f(\boldsymbol{X}) is convex. Our goal is to guarantee that this does not happen.

We remark that W⋆\boldsymbol{W}^{\star} is still a global minimizer of the factored problem (5) since f(X)f(\boldsymbol{X}) achieves its global minimum over the low-rank set of matrices at X⋆\boldsymbol{X}^{\star} and g(W)g(\boldsymbol{W}) also achieves its global minimum at W⋆\boldsymbol{W}^{\star}. The regularizer g(W)g(\boldsymbol{W}) is applied to force the difference between the two Gram matrices of U\boldsymbol{U} and V\boldsymbol{V} to be as small as possible. The global minimum of g(W)g(\boldsymbol{W}) is , which is achieved when U\boldsymbol{U} and V\boldsymbol{V} have the same Gram matrices, i.e., when W\boldsymbol{W} belongs to

III-B Main Results

Our main argument is that, under certain conditions on f(X)f(\boldsymbol{X}), the objective function ρ(W)\rho(\boldsymbol{W}) has no spurious local minima and satisfies the strict saddle property. This is equivalent to categorizing all the critical points into two types: 1) the global minima which correspond to the global solution of the original convex problem (1) and 2) strict saddles such that the Hessian matrix ∇2ρ(W)\nabla^{2}\rho(\boldsymbol{W}) evaluated at these points has a strictly negative eigenvalue. We formally establish this in the following theorem, whose proof is given in the next section.

For any μ>0\mu>0, each critical point W=[UV]\boldsymbol{W}=\begin{bmatrix}\boldsymbol{U}\\ \boldsymbol{V}\end{bmatrix} of ρ(W)\rho(\boldsymbol{W}) defined in (5) satisfies

Equation (7) shows that any critical point W\boldsymbol{W} belongs to E\mathcal{E} for the objective function in the factored problem (5) with any positive μ\mu. This demonstrates the reason for adding the regularizer g(U,V)g(\boldsymbol{U},\boldsymbol{V}). Thus, any iterative optimization algorithm converging to some critical point of ρ(W)\rho(\boldsymbol{W}) results in a solution within E\mathcal{E}. Furthermore, the strict saddle property along with the lack of spurious local minima ensures that a number of iterative optimization algorithms find the global minimum.

The constants appearing in Theorem 3 are not optimized. We use μ≤116α\mu\leq\frac{1}{16}\alpha simply to include μ=116\mu=\frac{1}{16} which is utilized for the matrix sensing problem in . If the ratio between the restricted strong convexity and smoothness constants βα≤1.4\frac{\beta}{\alpha}\leq 1.4, then we can show that ρ(W)\rho(\boldsymbol{W}) has no spurious local minima and obeys the strict saddle property for any μ≤14α\mu\leq\frac{1}{4}\alpha (where μ=14\mu=\frac{1}{4} is utilized for the matrix sensing problem in ). In all cases, a smaller μ\mu yields a more negative constant in (8); see Section IV for more discussion on this. This implies that when the restricted strong convexity constant α\alpha is not provided a priori, one can always choose a small μ\mu to ensure the strict saddle property holds, and hence guarantee the global convergence of many iterative optimization algorithms.

The constant 1.51.5 for the dynamic range βα\frac{\beta}{\alpha} in Theorem 3 is also not optimized and it is possible to slightly relax this constraint with more sophisticated analysis. However, the following example involving weighted symmetric matrix factorization implies that the room for improving this constant is rather limited. Let

Now consider the following weighted low-rank matrix factorization:

whose gradient ∇h(U)\nabla h(\boldsymbol{U}) and Hessian ∇2h(U)\nabla^{2}h(\boldsymbol{U}) are given by:

and λ2=4a>0\lambda_{2}=4a>0. We conclude that this U\boldsymbol{U} is a strict saddle point when a<2a<2 and a spurious local minimum when a>2a>2. This weighted symmetric matrix factorization problem (9) satisfies the restricted strong convexity and smoothness condition (3) with constants α=∥Ω∥min⁡2=1\alpha=\|\boldsymbol{\Omega}\|_{\min}^{2}=1 and β=∥Ω∥max⁡2=1+a\beta=\|\boldsymbol{\Omega}\|_{\max}^{2}=1+a (where ∥Ω∥min⁡\|\boldsymbol{\Omega}\|_{\min} and ∥Ω∥max⁡\|\boldsymbol{\Omega}\|_{\max} represent the smallest and largest entries in Ω\boldsymbol{\Omega}; see Section III-C). Thus, we have a counter example which demonstrates the existence of spurious local minima when βα>3\frac{\beta}{\alpha}>3.

We finally remark that although Theorem 3 requires the additional regularizer (4), empirical evidence (see experiments in Section V) shows we can get rid of this regularizer for many iterative algorithms with random initialization.

We prove Theorem 3 in Section IV. Before proceeding, we present two stylized applications of Theorem 3 in matrix sensing and weighted low-rank approximation.

III-C Stylized Applications

We first consider the implication of Theorem 3 in the matrix sensing problem where

holds for any n×mn\times m matrix X\boldsymbol{X} with rank⁡(X)≤r\operatorname{rank}(\boldsymbol{X})\leq r.

Note that, in this case, the gradient of f(X)f(\boldsymbol{X}) at X⋆\boldsymbol{X}^{\star} is

which implies that X⋆\boldsymbol{X}^{\star} is a critical point of f(X)f(\boldsymbol{X}). The Hessian quadrature form ∇2f(X)[Y,Y]\nabla^{2}f(\boldsymbol{X})[\boldsymbol{Y},\boldsymbol{Y}] for any n×mn\times m matrices X\boldsymbol{X} and Y\boldsymbol{Y} is given by

If A\mathcal{A} satisfies the 4r4r-restricted isometry property with constant δ4r\delta_{4r}, then f(X)f(\boldsymbol{X}) satisfies the (2r,4r)(2r,4r)-restricted strong convexity and smoothness condition (3) with constants α=1−δ4r\alpha=1-\delta_{4r} and β=1−δ4r\beta=1-\delta_{4r} since

for any rank-4r4r matrix Y\boldsymbol{Y}. Now, applying Theorem 3, we can characterize the geometry for the following matrix sensing problem with the factorization approach:

where g(U,V)g(\boldsymbol{U},\boldsymbol{V}) is the added regularizer defined in (4).

Suppose A\mathcal{A} satisfies the 4r4r-RIP with constant δ4r≤15\delta_{4r}\leq\frac{1}{5}, and set μ≤1−δ4r16\mu\leq\frac{1-\delta_{4r}}{16}. Then the objective function in (11) has no spurious local minima and satisfies the strict saddle property.

This result follows directly from Theorem 3 by noting that βα=1+δ4r1−δ4r≤1.5\frac{\beta}{\alpha}=\frac{1+\delta_{4r}}{1-\delta_{4r}}\leq 1.5 if δ4r≤15\delta_{4r}\leq\frac{1}{5}. We remark that Park et al. [19, Theorem 4.3] provided a similar geometric result for (11). Compared to their result which requires δ4r≤1100\delta_{4r}\leq\frac{1}{100}, our result has a much weaker requirement on the RIP of the measurement operator.

III-C2 Weighted Low-Rank Matrix Factorization

We now consider the implication of Theorem 3 in the weighted matrix factorization problem , where

Here Ω\boldsymbol{\Omega} is an n×mn\times m weight matrix consisting of positive elements and ∘\circ denotes the point-wise product between two matrices. In this case, the gradient of f(X)f(\boldsymbol{X}) at X⋆\boldsymbol{X}^{\star} is

which implies that X⋆\boldsymbol{X}^{\star} is a critical point of f(X)f(\boldsymbol{X}). The Hessian quadrature form ∇2f(X)[Y,Y]\nabla^{2}f(\boldsymbol{X})[\boldsymbol{Y},\boldsymbol{Y}] for any n×mn\times m matrices X\boldsymbol{X} and Y\boldsymbol{Y} is given by

Thus f(X)f(\boldsymbol{X}) satisfies the (2r,4r)(2r,4r)-restricted strong convexity and smoothness condition (3) with constants α=∥Ω∥min⁡2\alpha=\|\boldsymbol{\Omega}\|_{\min}^{2} and β=∥Ω∥max⁡2\beta=\|\boldsymbol{\Omega}\|_{\max}^{2} since

where ∥Ω∥min⁡\|\boldsymbol{\Omega}\|_{\min} and ∥Ω∥max⁡\|\boldsymbol{\Omega}\|_{\max} represent the smallest and largest entries in Ω\boldsymbol{\Omega}, respectively. Now we consider the following weighted matrix factorization problem:

where g(U,V)g(\boldsymbol{U},\boldsymbol{V}) is the added regularizer defined in (4). For an arbitrary weight matrix Ω\boldsymbol{\Omega}, it is proven that the weighted low-rank factorization can be NP-hard and has spurious local minima. When the elements in the weight matrix Ω\boldsymbol{\Omega} are concentrated, it is expected that (12) can be efficiently solved by a number of iterative optimization algorithms as it is close to an (unweighted) matrix factorization problem (where Ω\boldsymbol{\Omega} is a matrix of ones) which obeys the strict saddle property . The following result characterizes the geometric structure in the objection function of (12) by directly applying Theorem 3.

Suppose Ω\boldsymbol{\Omega} satisfies ∥Ω∥max⁡2∥Ω∥min⁡2≤1.5\frac{\|\boldsymbol{\Omega}\|_{\max}^{2}}{\|\boldsymbol{\Omega}\|_{\min}^{2}}\leq 1.5. Set μ≤∥Ω∥min⁡216\mu\leq\frac{\|\boldsymbol{\Omega}\|_{\min}^{2}}{16}. Then the objective function in (12) has no spurious local minima and satisfies the strict saddle property.

III-C3 1-bit Matrix Completion

for all (i,j)∈Ω(i,j)\in\Omega. Typical choices for qq include the logistic regression model where q(x)=ex1+exq(x)=\frac{e^{x}}{1+e^{x}} and the probit regression model where q(x)=1−Φ(−x/σ)=Φ(x/σ)q(x)=1-\Phi(-x/\sigma)=\Phi(x/\sigma). Here Φ\Phi is the cumulative distribution function (CDF) of a mean-zero Gaussian distribution with variance σ2\sigma^{2}. In , the authors attempt to recover X⋄\boldsymbol{X}^{\diamond} from the incomplete nonlinear measurements {Yij}(i,j)∈Ω\{Y_{ij}\}_{(i,j)\in\Omega} by minimizing the negative log-likelihood function

which results in a maximum likelihood (ML) estimate.

We note that FΩ,YF_{\Omega,\boldsymbol{Y}} is a convex function for both the logistic model and the probit model. The following result also establishes that FΩ,YF_{\Omega,\boldsymbol{Y}} satisfies the restricted strong convexity and smoothness condition if we observe full 1-bit measurements, i.e., Ω=[n]×[m]\Omega=[n]\times[m].

Then FΩ,YF_{\Omega,\boldsymbol{Y}} satisfies the restricted strong convexity and smoothness condition:

The proof of Lemma 1 is given in Appendix A. Now we consider the logistic regression model where q(x)=ex1+exq(x)=\frac{e^{x}}{1+e^{x}}.

Suppose Ω=[n]×[m]\Omega=[n]\times[m] and γ≤1.3\gamma\leq 1.3. Consider the logistic regression model where q(x)=ex1+exq(x)=\frac{e^{x}}{1+e^{x}}. Then FΩ,YF_{\Omega,\boldsymbol{Y}} satisfies the restricted strong convexity and smoothness condition with

Applying Lemma 1 with direct calculation gives

where q′(x)=ex(1+ex)2q^{\prime}(x)=\frac{e^{x}}{(1+e^{x})^{2}}. Now if we restrict ∥X∥∞≤1.3\|\boldsymbol{X}\|_{\infty}\leq 1.3, we have

Under the assumption that X⋄\boldsymbol{X}^{\diamond} is low-rank, a nuclear norm constraint is utilized in to force a low-rank solution. Corollary 3 implies that we can apply matrix factorization for 1-bit matrix recovery given that the elements of X\boldsymbol{X} are bounded. For the setting where Ω\Omega is only a subset of [n]×[m][n]\times[m], considered the 1-bit matrix completion problem with the rank constraint and established a stronger statistical recovery guarantee than that in . Empirical evidence (see and Section V-C) supports that matrix factorization also works for 1-bit matrix completion.

IV Proof of Theorem 3

In this section, we provide a formal proof of Theorem 3. The main argument involves showing that each critical point of ρ(W)\rho(\boldsymbol{W}) either corresponds to the global solution of (1) or is a strict saddle whose Hessian ∇2ρ(W)\nabla^{2}\rho(\boldsymbol{W}) has a strictly negative eigenvalue. Specifically, we show that W\boldsymbol{W} is a strict saddle by arguing that the Hessian ∇2ρ(W)\nabla^{2}\rho(\boldsymbol{W}) has a strictly negative curvature along Δ:=W−W⋆R\boldsymbol{\Delta}:=\boldsymbol{W}-\boldsymbol{W}^{\star}\boldsymbol{R}, i.e., [∇2ρ(W)](Δ,Δ)≤−τ∥Δ∥F2[\nabla^{2}\rho(\boldsymbol{W})](\boldsymbol{\Delta},\boldsymbol{\Delta})\leq-\tau\|\boldsymbol{\Delta}\|_{F}^{2} for some τ>0\tau>0. Here R\boldsymbol{R} is an r×rr\times r orthonormal matrix such that the distance between W\boldsymbol{W} and W⋆\boldsymbol{W}^{\star} rotated through R\boldsymbol{R} is as small as possible.

We first present some useful results. The (2r,4r)(2r,4r)-restricted strong convexity and smoothness assumption (3) implies the following isometry property, whose proof is given in Appendix B.

Suppose the function f(X)f(\boldsymbol{X}) satisfies the (2r,4r)(2r,4r)-restricted strong convexity and smoothness condition (3) with positive α\alpha and β\beta. Then for any n×mn\times m matrices Z,G,H\boldsymbol{Z},\boldsymbol{G},\boldsymbol{H} of rank at most 2r2r, we have

We remark that Lemma 2 is a variant of [19, Lemma 3.2]. While the result there requires the 4r4r-RIP condition of the objective function, our result depends on the (2r,4r)(2r,4r)-restricted strong convexity and smoothness condition. Our result is also slightly tighter than [19, Lemma 3.2].

If C=0\boldsymbol{C}={\bf 0}, then we have

We present one more useful result in the following Lemma.

Finally, we provide the gradient and Hessian expressions for ρ(W)\rho(\boldsymbol{W}). The gradient of ρ(W)\rho(\boldsymbol{W}) is given by

IV-B The Formal Proof

Any critical point W\boldsymbol{W} of ρ(W)\rho(\boldsymbol{W}) satisfies ∇ρ(W)=0\nabla\rho(\boldsymbol{W})={\bf 0}, i.e.,

Now we turn to prove the strict saddle property and that there are no spurious local minima.

First, note that as guaranteed by Proposition 1, X⋆\boldsymbol{X}^{\star} is the unique n×mn\times m matrix with rank at most rr. Also the gradient of f(X)f(\boldsymbol{X}) vanishes at X⋆\boldsymbol{X}^{\star} since (1) is an unconstraint optimization problem. Denote the set of critical points of ρ(W)\rho(\boldsymbol{W}) by

We separate C\mathcal{C} into two subsets:

satisfying C=C1∪C2\mathcal{C}=\mathcal{C}_{1}\cup\mathcal{C}_{2}. Since any critical point W\boldsymbol{W} satisfies (18), g(W)g(\boldsymbol{W}) achieves its global minimum at W\boldsymbol{W}. Also f(X)f(\boldsymbol{X}) achieves its global minimum at X⋆\boldsymbol{X}^{\star}. We conclude that W\boldsymbol{W} is the globally optimal solution of ρ\rho for any W∈C1\boldsymbol{W}\in\mathcal{C}_{1}. If we show that any W∈C2\boldsymbol{W}\in\mathcal{C}_{2} is a strict saddle, then we prove that there are no spurious local minima as well as the strict saddle property. Thus, the remaining part is to show that C2\mathcal{C}_{2} is the set of strict saddles.

To show that C2\mathcal{C}_{2} is the set of strict saddles, it is sufficient to find a direction Δ\boldsymbol{\Delta} along which the Hessian has a strictly negative curvature for each of these points. We construct Δ=W−W⋆R\boldsymbol{\Delta}=\boldsymbol{W}-\boldsymbol{W}^{\star}\boldsymbol{R}, the difference from W\boldsymbol{W} to its nearest global factor W⋆\boldsymbol{W}^{\star}, where

The following result (which is proved in Appendix E) states that Π1\Pi_{1} is strictly negative, while the remaining terms are relatively small, though they may be nonnegative:

where (i)(i) utilizes Lemmas 2 and 4, (ii)(ii) utilizes the following inequality (which is proved in Appendix F)

and (ii)(ii) holds because βα≤1.5\frac{\beta}{\alpha}\leq 1.5 and μ≤116α\mu\leq\frac{1}{16}\alpha. Thus, if X≠X⋆\boldsymbol{X}\neq\boldsymbol{X}^{\star}, [∇2ρ(X)](Δ,Δ)\left[\nabla^{2}\rho(\boldsymbol{X})\right](\boldsymbol{\Delta},\boldsymbol{\Delta}) is always negative. This implies that W\boldsymbol{W} is a strict saddle.

To complete the proof, we utilize Lemma 3 to further bound the last term in (23):

From (23), we observe that a smaller μ\mu yields a more negative bound on [∇2ρ(X)](Δ,Δ)\left[\nabla^{2}\rho(\boldsymbol{X})\right](\boldsymbol{\Delta},\boldsymbol{\Delta}). This can be explained intuitively as follows. First note that any critical point W\boldsymbol{W} satisfies (18) provided μ>0\mu>0, no matter how large or small μ\mu is. The Hessian information about g(W)g(\boldsymbol{W}) is represented by the terms Π3\Pi_{3} and Π4\Pi_{4}. We have

where the last line holds since for any r×rr\times r matrix A\boldsymbol{A},

Thus the Hessian of ρ\rho evaluated at any critical point W\boldsymbol{W} is a PSD matrixThis can also be observed since any critical point W\boldsymbol{W} is a global minimum point of ρ(W)\rho(\boldsymbol{W}), which directly indicates that ∇2ρ(W)⪰0\nabla^{2}\rho(\boldsymbol{W})\succeq{\bf 0}. instead of having a negative eigenvalue. In low-rank, PSD matrix optimization problems, the corresponding objective function (without any regularizer such as g(W)g(\boldsymbol{W})) is proved to have the strict saddle property . Therefore, h(W)h(\boldsymbol{W}) is also expected to have the strict saddle property, and so is ρ(W)\rho(\boldsymbol{W}) when μ\mu is small, i.e., the Hessian of g(W)g(\boldsymbol{W}) has little influence on the Hessian of ρ(W)\rho(\boldsymbol{W}) when μ\mu is small. Our results also indicate that when the restricted strict convexity constant α\alpha is not provided a priori, we can always choose a small μ\mu to ensure the strict saddle property of ρ(W)\rho(\boldsymbol{W}) is met, and hence we are guaranteed the global convergence of a number of local search algorithms applied to (5).

V Experiments

In this section, we present a set of experiments on matrix sensing, matrix completion, and 1-bit matrix completion to demonstrate the performance of iterative algorithms for low-rank matrix optimization. Unless noted otherwise, we denote the matrix factorization approach by NVX and use the minFunc packageSoftware available at https://www.cs.ubc.ca/∼\simschmidtm/Software/minFunc.html to perform the local search algorithms for the factored problem.

where the entries of each n×mn\times m matrix Yi\boldsymbol{Y}_{i} are independent and identically distributed (i.i.d.) normal random variables with zero mean and variance 1p\frac{1}{p} for i∈{1,2,…p}i\in\{1,2,\ldots p\}. For each pair of rr and the number of measurements, 10 Monte Carlo trials are carried out and for each trial, and we claim matrix recovery to be successful if the relative reconstruction error satisfies

where we denote by X^\widehat{\boldsymbol{X}} the reconstructed matrix. Figure 1 displays the phase transition for factorized gradient descent starting from a random initialization, the singular value projection (SVP) method proposed in which requires a SVD in each iteration, and the convex approach which solves

We see that there are only negligible differences between the different approaches for matrix sensing; these approaches also have very similar performance guarantees when the Gaussian sensing operator A\mathcal{A} satisfies the RIP . We note that with or without the regularizer gg as defined in (4), local search algorithms have similar performance with random initialization. Hence, throughout all of the experiments, we simply discard the regularizer gg, but we stress that identical performance is observed if we have this regularizer gg.

V-B Matrix Completion

We compare the performance of the matrix factorization approach with SVP , the convex approach, and singular value thresholdingSoftware available at http://svt.stanford.edu/ (SVT) for matrix completion where we want to recover a low-rank matrix X⋆\boldsymbol{X}^{\star} from incomplete measurements {Xij⋆}(i,j)∈Ω\{X^{\star}_{ij}\}_{(i,j)\in\Omega}, where Ω⊂[n]×[m]\Omega\subset[n]\times[m]. Let PΩ\mathcal{P}_{\Omega} denote the projection onto the index set Ω\Omega. The convex approach (denoted by CVX) attempts to use the nuclear norm as a convex relaxation of the rankness and solves

Though PΩ\mathcal{P}_{\Omega} does not satisfy the rr-RIP (10) for all low-rank matrices X\boldsymbol{X}, it satisfies the RIP when restricted to low-rank incoherent matrices.

[48, Theorem 4.2] Without loss of generality, assume n≥mn\geq m. There exists a constant C≥0C\geq 0 such that for Ω∈[n]×[m]\Omega\in[n]\times[m] chosen according to the Bernouli model with density greater than Cu2r2log⁡n/δ2mCu^{2}r^{2}\log n/\delta^{2}m, with probability at least 1−e−nlog⁡n1-e^{-n\log n}, the RIP holds for all μ\mu-incoherent matrices X\boldsymbol{X} of rank at most rr.

Thus, if local search algorithms (such as gradient descent) start with a random initialization and the iterates remain incoherent, then Theorem 3 guarantees the global convergence of the matrix factorization approach with these algorithms. We note that this hypothesis is also required for SVP . Though we can add a regularizer for incoherence as in , empirical evidence supports this hypothesis that the iterates in gradient descent are incoherent.

In the first set of experiments, we set n=m=100n=m=100 and vary the rank rr from 11 to 3030. Similar to the setup for matrix sensing in Section V-A, we generate a rank-rr random matrix and randomly obtain pp entries, i.e., ∣Ω∣=p|\Omega|=p. Figure 3 displays the phase transition for gradient descent with a random initialization, SVP , singular value thresholding (SVT) , and the convex approach. As can been seen, the matrix factorization approach has similar phase transition to SVP, and is slightly better than SVT and the convex approach in terms of the number of measurements needed for successful recovery.

In the second set of experiments, we set r=5r=5 and p=3r(2n−r)p=3r(2n-r) (3 times the number of degrees of freedom within a rank-rr n×nn\times n matrix), and vary nn from 4040 to 51205120. We compare the time needed for the four approaches in Figure 4; our matrix factorization approach is much faster than the other methods. The time savings for the matrix factorization approach comes from avoiding performing the SVD, which is needed both for SVT and SVP in each iteration. We also observe that convex approach has the highest computational complexity and is not scalable (which is the reason that we only present its time for nn up to 640640).

V-C 1-bit Matrix Completion

In the last set of experiments, we compare the performance of the matrix factorization approach with the convex approachSoftware available at http://mdav.ece.gatech.edu/software/ in for 1-bit matrix completion. We first note that to make the recovery problem well-posed, a constraint on ∥X∥∞\|\boldsymbol{X}\|_{\infty} (the entry-wise maximum of the matrix X\boldsymbol{X}) is applied in to require that the matrix is not too “spiky”. Instead of using the constraint on ∥X∥∞\|\boldsymbol{X}\|_{\infty}, we add a smooth regularizer ∥X∥F2\|\boldsymbol{X}\|_{F}^{2} and turn to minimize the following objective function

To evaluate the performance of this factorization approach on 1-bit matrix completion, we generate n×rn\times r matrices U⋄\boldsymbol{U}^{\diamond} and V⋄\boldsymbol{V}^{\diamond} with entries drawn i.i.d. from a uniform distribution on [−12,12][-\frac{1}{2},\frac{1}{2}] and construct a random n×nn\times n matrix X⋄\boldsymbol{X}^{\diamond} with rank rr. Similar to the setup in , the matrix is then scaled so that ∥X⋄∥=1\|\boldsymbol{X}^{\diamond}\|=1. We obtain 1-bit observations {Yi,j}(i,j)∈Ω\{Y_{i,j}\}_{(i,j)\in\Omega} by adding Gaussian noise of variance σ2\sigma^{2} and recording the sign of the resulting value (15), where the subset of indices Ω\Omega is chosen at random with E⁡∣Ω∣=p\operatorname{E}|\Omega|=p. We compare the performance of the factorization approach and the convex approach over a range of different values of nn, pp, rr or σ\sigma. Figures 5(a)-(d) show the normalized squared Frobenius norm of the error ∥X^−X⋄∥∥X⋄∥F2\frac{\|\widehat{\boldsymbol{X}}-\boldsymbol{X}^{\diamond}\|}{\|\boldsymbol{X}^{\diamond}\|_{F}^{2}} (where X^\widehat{\boldsymbol{X}} denotes the reconstructed matrix) and average the results over 10 draws of Monte Carlo trials. We observe that matrix factorization approach has slightly better performance than the convex approach for 1-bit matrix completion . Note that this phenomenon (the factorization approach having better performance) is also observed in . We repeat these experiments but obtaining 1-bit observations with the logistic regression model where g(x)=ex1+exg(x)=\frac{e^{x}}{1+e^{x}} for (15) and display the results in Figure 6.

VI Conclusion

This paper considers low-rank matrix optimization on general (nonsymmetric and rectangular) matrices with general objective functions. By focusing on general objective functions, we provide a unifying framework for low-rank matrix optimizations with the factorization approach. Although the resulting optimization problem is not convex, we show that the reformulated objection function has a simple landscape: there are no spurious local minima and any critical point not being a local minimum is a strict saddle such that the Hessian evaluated at this point has a strictly negative eigenvalue. These properties guarantee that a number of iterative optimization algorithms (such as gradient descent and the trust region method) will converge to the global optimum from a random initialization.

Appendix A Proof of Lemma 1

We compute the partial derivative of FΩ,YF_{\Omega,\boldsymbol{Y}} in terms of Xi,jX_{i,j} as

Appendix B Proof of Proposition 2

This proof follows similar steps to the proof of [51, Lemma 2.1]. First note that the bilinear form [∇2f(Z)](G,H)=∑i,j,k,l∂2f(Z)∂Zij∂ZklGijHkl[\nabla^{2}f(\boldsymbol{Z})](\boldsymbol{G},\boldsymbol{H})=\sum_{i,j,k,l}\frac{\partial^{2}f(\boldsymbol{Z})}{\partial\boldsymbol{Z}_{ij}\partial\boldsymbol{Z}_{kl}}\boldsymbol{G}_{ij}\boldsymbol{H}_{kl} implies [∇2f(Z)](G,H)[\nabla^{2}f(\boldsymbol{Z})](\boldsymbol{G},\boldsymbol{H}) is invariant under all scalings for both G\boldsymbol{G} and H\boldsymbol{H}, i.e.,

Now suppose both G\boldsymbol{G} or H\boldsymbol{H} are nonzero. By the scaling invariance property of both sides in (3), we assume ∥G∥F=∥H∥F=1\|\boldsymbol{G}\|_{F}=\|\boldsymbol{H}\|_{F}=1 without loss of generality. Note that the (2r,4r)(2r,4r)-restricted strong convexity and smoothness condition (3) implies

Appendix C Proof of Lemma 2

It follows from (19) and (20) that any critical point W\boldsymbol{W} satisfies

On the other hand, we give an upper bound on the right hand side of (29):

Appendix D Proof of Lemma 3

When C≠0\boldsymbol{C}\neq{\bf 0}, the proof follows directly from the following results.

If C=0\boldsymbol{C}={\bf 0}, then we have

Appendix E Proof of (22)

We prove the upper bounds for the four terms as follows.

Bounding term Π1\Pi_{1}: Utilizing the fact that ΔU=U−U⋆R\boldsymbol{\Delta}_{\boldsymbol{U}}=\boldsymbol{U}-\boldsymbol{U}^{\star}\boldsymbol{R} and ΔV=V−V⋆R\boldsymbol{\Delta}_{\boldsymbol{V}}=\boldsymbol{V}-\boldsymbol{V}^{\star}\boldsymbol{R}, we have

where (i)(i) follows from (19) and (20), (ii)(ii) utilizes ∇f(X⋆)=0\nabla f(\boldsymbol{X}^{\star})={\bf 0}, and (iii)(iii) follows by using the (2r,4r)(2r,4r)-restricted strict convexity property (3):

where the first line follows from the integral form of the mean value theorem for vector-valued functions, and the second line uses the fact that both tX+(1−t)X⋆t\boldsymbol{X}+(1-t)\boldsymbol{X}^{\star} and X−X⋆\boldsymbol{X}-\boldsymbol{X}^{\star} have rank at most 2r2r, and the (2r,4r)(2r,4r)-restricted strong convexity of the Hessian ∇2f(⋅)\nabla^{2}f(\cdot).

Bounding term Π2\Pi_{2}: By the smoothness condition (3), we have

Appendix F Proof of (24)

To show (24), expanding the left hand side of (24), it is equivalent to show

Thus, we obtain (24) by noting that the above equation is equivalent to

References