# Escaping Spurious Local Minima of Low-Rank Matrix Factorization Through
Convex Lifting

Ching-pei Lee
Academia Sinica
Taipei, Taiwan
<EMAIL_ADDRESS>
&Ling Liang
Department of Mathematics
National University of Singapore
Singapore
<EMAIL_ADDRESS>Tianyun Tang
Department of Mathematics
National University of Singapore
Singapore
<EMAIL_ADDRESS>
&Kim-Chuan Toh
Department of Mathematics & Institute of Operations Research and Analytics
National University of Singapore
Singapore
<EMAIL_ADDRESS>
###### Abstract

This work proposes a rapid global solver for nonconvex low-rank matrix
factorization (MF) problems that we name MF-Global. Through convex lifting
steps, our method efficiently escapes saddle points and spurious local minima
ubiquitous in noisy real-world data, and is guaranteed to always converge to
the global optima. Moreover, the proposed approach adaptively adjusts the rank
for the factorization and provably identifies the optimal rank for MF
automatically in the course of optimization through tools of manifold
identification, and thus it also spends significantly less time on parameter
tuning than existing MF methods, which require an exhaustive search for this
optimal rank. On the other hand, when compared to methods for solving the
lifted convex form only, MF-Global leads to significantly faster convergence
and much shorter running time. Experiments on real-world large-scale
recommendation system problems confirm that MF-Global can indeed effectively
escapes spurious local solutions at which existing MF approaches stuck, and is
magnitudes faster than state-of-the-art algorithms for the lifted convex form.

##  1 Introduction

The matrix factorization (MF) problem that recovers a low-rank matrix $X$ from
noisy and partially observed data by minimizing a loss with respect to its
factorization $(W,H)$ such that $X=WH^{\top}$ finds thriving interest in
various applications, such as recommendation systems [42], phase retrieval
[11], and the Internet of things [39]. As one of the most popular and simple
nonconvex problems in machine learning and the simplest special case of neural
networks (with only one hidden layer and the activation function being the
identity map), lately it also receives growing theoretical interests in
studying both its optimization landscape [19, 48, 20, 32, 33, 17] and whether
simple methods for it can make $(W,H)$ efficiently converge to a local optimum
without being stuck at strict saddle points [34, 13, 35, 56]. In particular,
it has been shown in recent works like [19, 20, 48, 13, 56] that for a
quadratic or linear loss, either under some stringent and non-verifiable
assumptions on the global solutions as well as the observed data entries, or
provided that a high enough rank is given, methods for MF can find the global
optima with high probability. However, these theoretical results either need
to overestimate the rank to a certain degree, or require that the entries
observed are noiseless and decided by a uniform random sampling. The former
means that one needs to start from an unnecessarily high rank, in which case
the practical efficiency of matrix factorization approaches is by and large
eliminated, making such theoretically guaranteed approaches infeasible for
large-scale problems in real-world applications because the increased rank
will lead to prohibitively high cost per iteration. The latter requirement is
also impractical for two reasons. First, the uniform sampling assumption
clearly violates the real scenario in recommendation systems such that entries
recommended to users (and thus those of which we might get an observation) are
actively selected by mechanisms of an existing recommendation system. Second,
measurement errors and data corruption are omnipresent in real-world
applications, so the noiseless assumption can hardly be satisfied.

Recently, it has also been shown by [53] that there indeed exists a family of
MF problems whose number of spurious local minima grows exponentially with the
problem dimension, which validates that those theoretical guarantees above
could fail when the problem to solve fails some of the assumptions in those
theoretical analyses. In practice, it is unlikely that any method considers
the nonconvex MF formulation only could escape all such spurious local minima
to reach the global optima, even if such a method avoids all strict saddle
points. Moreover, those theoretical analyses all focus on very specific
problem cases to utilize properties that cannot be extended to general
problems, and thus a practically efficient and general approach for finding
global minima of low-rank matrix factorization problems is still in an urgent
need. We hence propose in this work an efficient and practical algorithmic
framework for MF to escape from such stationary points that are not global
optima by utilizing convex lifting steps. Our framework is a flexible one that
can be combined with any algorithm for MF to easily provide it with an
efficient escape from spurious stationary points and hence global convergence
to global optima, and we name this method MF-Global.

Alternative to the nonconvex MF problem that minimizes $(W,H)$, low-rank
matrix completion (MC) that directly minimizes their product $X$ can lead to a
convex formulation, so convergence to the global optima is guaranteed [10, 9,
41, 40], and such a convex formulation forms the cornerstone of the lifting
steps in our framework. However, such convergence guarantee for global
optimality comes with much higher computational and spatial costs. In
particular, when the matrix dimension is high, it is impractical to explicitly
compute and store $X$, and thus in implementation, one still needs to form
some factorization of $X$. To make such factorization of a reasonable size, we
must then enforce all iterates to be of low-rank, and thus the nuclear norm
regularization is needed. However, accurate computation of the nuclear norm
requires a partial singular value decomposition (SVD) of the positive singular
values, whose exact computation from scratch could be expensive. Although the
state of the art for convex MC [49, 23, 55] works around this difficulty by
using either the Lanczos or the power method to compute approximate SVDs in a
lower cost through matrix-vector products, these approximate SVDs are still
rather expensive, and the practical efficiency of such approaches for MC is
hardly comparable to that of methods for MF. In contrast, our framework
utilizes MF methods as much as possible, and only resorts to a convex lifting
step when we need to escape a spurious stationary point, so the encounters
with (approximate) SVDs is reduced to minimal. We also leverage on recent
advancements in numerical linear algebra to incorporate a limited-memory
Krylov subspace method [36] to make our computation of approximate SVDs much
more efficient than existing MC approaches. The resulting MF-Global is
therefore magnitudes faster than the state of the art for MC in our numerical
experiments with real-world large-scale recommendation system data.

The main goal of this work is to propose a highly efficient framework for MF
problems that absorbs the advantages of both the nonconvex but simpler MF form
and the convex MC form. Like methods for the convex MC, our framework
possesses guaranteed convergence to the global optima, while it fully utilizes
the simpler form of nonconvex MF without the expensive nuclear norm to obtain
better parallelism and a much lower per iteration cost. Partly inspired by
[54] that accelerates certain SDP problems by polynomial optimization, our MF-
Global alternates between solving the smooth MF problem and a convex lifting,
or escape, step that optimizes the convex MC problem with only one iteration
to escape from spurious stationary points in the nonconvex formulation and to
decide the rank of the matrices in the next MF round.

Another advantage of MF-Global over those for MF is the ability to find the
optimal rank without parameter search. Existing methods for MF need to pre-
specify the rank of $W$ and $H$ as a parameter in the optimization problem,
and the optimal value for $k$ is found through an exhaustive search such that
multiple MF problems with different ranks are solved, and thus efficiency of
MF methods over MC approaches is not as significant as it would appear due to
this additional parameter search time. On the other hand, through tools of
manifold identification [29, 22, 31], our method is guaranteed to quickly and
automatically fixate the rank at the MF steps to that of the global optimum to
which the iterates converge, and thus MF-Global is empirically much more
efficient than using a MF solver only through saving the expensive cost on
parameter search.

Contributions. The main contributions of this work is four-fold.

  1. 1.

_Global optimality guarantee with no assumption_ : We propose MF-Global, a
universally safeguarded approach for MF that can effectively escape all
spurious stationary points and hence possesses convergence guarantees to the
global optima, without any requirement on the data or the MF solver.

  2. 2.

_Minimizing the need for expensive SVDs_ : Through the steps that focus on the
easier-to-optimize MF form, MF-Global requires significantly fewer times of
computing partial SVDs than MC approaches, and thus achieves better efficiency
and parallelism.

  3. 3.

_Provable rank identification_ : We show that after a finite number of
iterations, all iterates of MF-Global automatically have the same rank as the
global optimum that they converge to. By knowing the optimal rank in advance,
MF-Global conducts approximate SVDs more efficiently than existing MC
approaches because computation of the singular vectors/values beyond the
optimal rank is avoided. MF-Global also totally removes the parameter search
time needed in MF approaches for tuning the rank.

  4. 4.

_Superior empirical performance_ : Numerical results show that on real-world
recommendation system data whose dimension is as high as a million by hundreds
of thousands, MF-Global is magnitudes faster than existing methods for MC, and
is also superior than the state of the art for MF by escaping spurious
stationary points at which the MF approach gets stuck at.

Organization. This paper is organized as follows. In Section 2, we describe
the problem formulation and details of MF-Global. Theoretical results are then
stated in Section 3 with their proofs in the appendix. Section 4 reviews
existing works and contrasts our approach with them. Experiments in Section 5
provides empirical evidence of the superior efficiency of the proposed MF-
Global over existing MF and MC solvers. Finally, Section 6 concludes this
work. Experiments and implementation details, additional experiments, and
proofs are in the appendices.

##  2 Algorithmic framework

Given the problem dimension $m,n>0$, we consider the following matrix
factorization problem:

$\min_{W\in\mathbb{R}^{m\times k},H\in\mathbb{R}^{n\times
k}}\left(F\left(W,H\right)\coloneqq
f\left(WH^{\top}\right)+\frac{\lambda}{2}\left({\left\|{W}\right\|}_{F}^{2}+{\left\|{H}\right\|}_{F}^{2}\right)\right),$
(MF)

where $\lambda>0$ and $k\in\mathbb{N}$ are pre-specified parameters and
${\left\|{\cdot}\right\|}_{F}$ is the Frobenius norm. We assume without loss
of generality $m\leq n$, which can easily be achieved by a matrix transpose.
On the other hand, the convex lifting formulation we consider is also known as
the matrix completion problem.

$\min\nolimits_{X\in\mathbb{R}^{m\times n}}\quad F\left(X\right)\coloneqq
f\left(X\right)+\lambda{\left\|{X}\right\|}_{*},$ (MC)

where ${\left\|{\cdot}\right\|}_{*}$ is the nuclear norm that sums up the
singular values of the given matrix. For Equation MF and Equation MC, we
assume the loss $f$ is convex and has $L$-Lipschitz-continuous gradient for
some $L>0$, and although $X$ is too large to be stored explicitly, we can
afford either forming $\nabla f(X)$ or computing $\nabla f(X)v$ for any given
$v$ efficiently. Therefore, both terms in Equation MF are smooth, so many
efficient and highly parallelizable algorithms can be applied to solve
Equation MF with a low per-iteration cost, but the nonconvex nature of
$F(W,H)$ makes it hard to find the global optima of Equation MF. On the other
hand, Equation MC is more costly to minimize because even evaluating the
objective value requires a round of expensive SVD, which is often the
bottleneck of methods for Equation MC.

The standard approach for solving Equation MC is the proximal gradient (PG)
method, whose exact form in the case of Equation MC with a given stepsize
$\alpha_{t}>0$ is as follows.

$\displaystyle X_{t+1}$ $\displaystyle=\mbox{\rm
prox}_{{\alpha_{t}}\lambda{\left\|{\cdot}\right\|}_{*}}\left(X_{t}-{\alpha_{t}}\nabla
f\left(X_{t}\right)\right)$ (1)
$\displaystyle=\operatorname*{\arg\min}\nolimits_{Y}\,\big{(}Q^{{\alpha_{t}}}_{X_{t}}\left(Y\right)\coloneqq\langle\nabla
f(X_{t}),\,Y-X_{t}\rangle+{(2\alpha_{t})^{-1}}{\left\|{Y-X_{t}}\right\|}^{2}_{F}+\lambda{\left\|{Y}\right\|}_{*}\big{)},$

where $\langle A,\,B\rangle=\text{trace}(A^{\top}B)$ is the inner product of
matrices $A,B$, and given any function $g$, $\mbox{\rm prox}_{g}$ is its
proximal operator: $\mbox{\rm
prox}_{g}(x)\coloneqq\operatorname*{\arg\min}_{y}\\{{\left\|{x-y}\right\|}^{2}/2+g(y)\\}$.
When $g=\beta{\left\|{\cdot}\right\|}_{*}$ for $\beta>0$, it is known that for
any $\hat{Z}_{t}$ with $\mbox{\rm rank}(\hat{Z}_{t})=k_{t}$ and SVD
$\hat{Z}_{t}=U_{t}\mbox{\rm diag}(\hat{\Sigma}_{t})V_{t}^{\top}$, where
$\hat{\Sigma}_{t}\in\mathbb{R}^{k_{t}}_{+}$, $\mathbb{R}^{k_{t}}_{+}$ is the
nonnegative orthant of $\mathbb{R}^{k_{t}}$, and $U_{t}\in\mathbb{R}^{m\times
k_{t}},V_{t}\in\mathbb{R}^{n\times k_{t}}$ are orthonormal (matrices with
orthonormal columns), the proximal operation has a closed-form solution [30]

$\mbox{\rm
prox}_{\beta{\left\|{\cdot}\right\|}_{*}}\left(\hat{Z}_{t}\right)=U_{t}\mbox{\rm
diag}\left(\Sigma_{t}\right)V_{t}^{\top},\quad\Sigma_{t}\coloneqq[\hat{\Sigma}_{t}-\beta
e]_{+},$ (2)

where $e\coloneqq(1,\dotsc,1)$, and $[\cdot]_{+}$ is the projection onto
$\mathbb{R}^{k_{t}}_{+}$. Computing the exact SVD of $X_{t}-{\alpha_{t}}\nabla
f(X_{t})$ needed in Equation 2 is of cost $O(m^{2}n)$, which can be
prohibitively large when we are not even able to afford to explicitly form
$X$. Therefore, we will conduct inexact PG steps whose SVD part is approximate
but the truncation in Equation 2 using the approximate SVD is exact.

Our method MF-Global is an iterative process such that at the $t$th iteration
with the current iterate $X_{t}$, it conducts the following two operations in
order.

  1. 1.

MF phase: Solve Equation MF using a descent method to a certain stopping
tolerance, starting from the point $(\tilde{W}_{t},\tilde{H}_{t})$ obtained
from the SVD of $X_{t}$ computed in the inexact PG step, to get
$(W_{t},H_{t})$ satisfying a certain stopping criterion. We call such an
intermediate iterate $\tilde{X}_{t}\coloneqq W_{t}H_{t}^{\top}$.

  2. 2.

Convex lifting step: Conduct an inexact PG step at $\tilde{X}_{t}$ with error
measure $\epsilon_{t}>0$ to obtain the next iterate $X_{t+1}$ as follows.

$Z_{t}\coloneqq\tilde{X}_{t}-{\alpha_{t}}\nabla f(\tilde{X}_{t}),\quad
X_{t+1}\in\mbox{\rm
prox}^{\epsilon_{t}}_{{\alpha_{t}}\lambda{\left\|{\cdot}\right\|}_{*}}\left(Z_{t}\right)\coloneqq\big{\\{}Y\mid\min_{g\in\partial
Q^{{\alpha_{t}}}_{\tilde{X}_{t}}\left(Y\right)}{\left\|{g}\right\|}\leq\epsilon_{t}\big{\\}}.$
(3)

In the convex lifting step, standard PG approaches for Equation MC fixes
$\alpha_{t}\equiv\alpha\in(0,2/L)$. However, it is known that the convergence
speed of standard PG is slow, thus we follow [52] to use the Barzilai-Borwein
(BB) initialization [3] for the stepsize with linesearch to accelerate the
practical convergence of the convex lifting step. Given
$\alpha_{\max}\geq\alpha_{\min}>0$, we first compute

$\alpha_{t}^{\text{BB}}\coloneqq\max\big{\\{}\alpha_{\min},\min\big{\\{}\alpha_{\max},\frac{\langle
X_{t}-X_{t-1},\,\nabla f(X_{t})-\nabla
f(X_{t-1})\rangle}{{\left\|{X_{t}-X_{t-1}}\right\|}_{F}^{2}}\big{\\}}\big{\\}}.$
(4)

But different from [52], our backtracking linesearch does not search for an
$\alpha_{t}$ that gives sufficient descent. Instead, given
$\beta,\delta\in(0,1)$, we find the smallest nonnegative integer $i$ such that

$\alpha_{t}=\alpha_{t}^{\text{BB}}\beta^{i},\quad f\left(X_{t+1}\right)\leq
f\left(X_{t}\right)+\langle\nabla
f\left(X_{t}\right),\,X_{t+1}-X_{t}\rangle+\delta\alpha_{t}^{-1}{\left\|{X_{t+1}-X_{t}}\right\|}^{2}.$
(5)

From the $L$-Lipschitz continuity of $\nabla f,$ Equation 5 is satisfied for
any $\alpha_{t}\leq 2\delta/L$, so the backtracking procedure will terminate
in finite steps.

For the proximal part in Equation 3, we use the following procedure to
generate an iterate satisfying Equation 3. First, we compute a rank-$k_{t}$
approximate eigendecomposition of $Z_{t}Z_{t}^{\top}$ to get

$Z_{t}Z_{t}^{\top}\approx\hat{U}_{t}\mbox{\rm
diag}\left(\bar{\Sigma}^{t}\odot\bar{\Sigma}^{t}\right)\hat{U}_{t}^{\top}$ (6)

for $\bar{\Sigma}^{t}\in\mathbb{R}^{k_{t}}_{+}$ and $\hat{U}_{t}$ orthonormal,
where $\odot$ denotes entry-wise multiplication of two vectors. We then
conduct an exact SVD on $\hat{U}_{t}^{\top}Z_{t}\in\mathbb{R}^{k_{t}\times n}$
to obtain $\hat{U}_{t}^{\top}Z_{t}=\tilde{U}_{t}\mbox{\rm
diag}\big{(}\hat{\Sigma}^{t}\big{)}V_{t}^{\top}$ with cost
$O(k_{t}^{3}+nk_{t}^{2})$, which is much cheaper than SVD for $Z_{t}$ when
$k_{t}\ll m$. The approximate SVD is then obtained by

$Z_{t}\approx\hat{Z}_{t}\coloneqq U_{t}\mbox{\rm
diag}\big{(}\hat{\Sigma}^{t}\big{)}V_{t}^{\top},\quad
U_{t}\coloneqq\hat{U}_{t}\tilde{U}_{t}.$ (7)

Clearly, $U_{t}$ and $V_{t}$ are both orthonormal, so this is actually the
exact SVD of $\hat{Z}_{t}=\hat{U}_{t}\hat{U}_{t}^{\top}Z_{t}$. The inexact PG
step is then finished by applying Equation 2 to $\hat{Z}_{t}$ in Equation 7.
Note that we only store $(U_{t},\Sigma^{t},V_{t})$ without explicitly forming
$X_{t+1}$ throughout the process. Here we slightly abuse the notation to let
$\Sigma^{t}$ denote only the coordinates with nonzero values, and let $U_{t}$
and $V_{t}$ contain only the corresponding columns. This procedure is
summarized in Algorithm 1.

For the approximate eigendecomposition, as mentioned in Section 1, we cannot
directly form $\tilde{X}_{t}$ or $Z_{t}$ to conduct exact eigendecomposition
needed in Equation 2. Fortunately, given the factorized smaller matrices
$W_{t}$ and $H_{t}$, $\tilde{X}_{t}v$ and $Z_{t}v$ for any $v$ can be computed
efficiently. The state of the art for Equation MC hence utilizes this property
to use the power method to compute an approximate eigendecomposition of
$Z_{t}Z_{t}^{\top}$. Given an initial point $R_{t}$, the power method computes
$\hat{U}_{t}^{(0)}=\text{orth}(R_{t})$ and
$\hat{U}_{t}^{(i)}=\text{orth}(Z_{t}Z_{t}^{\top}\hat{U}_{t}^{(i-1)})$ for
$i>0$ where for any matrix $A$, $\text{orth}(A)$ is its orthogonalization, and
then takes $\hat{U}_{t}=\hat{U}_{t}^{(T_{t})}$ for some $T_{t}$. An
inefficient part of the power method is that $\hat{U}_{t}^{(j)}$ for $j<i-1$
is not used in obtaining $\hat{U}_{t}^{(i)}$. Thus, we adopt the limited-
memory Krylov subspace method of [36] instead to form the best approximation
$U_{t}^{(i)}$ using both $\text{orth}(Z_{t}Z_{t}^{\top}U_{t}^{(i-1)})$ and
$\\{\hat{U}_{t}^{(i-1)},\dotsc,U_{t}^{(i-M)}\\}$ for some $M\geq 0$ to improve
the slow convergence of the power method. When compared with the power method,
this approach takes slightly more storage and computation per iteration, but
lead to much improved convergence speed and thus greatly reduces the time cost
of the approximate eigendecomposition. For initializing $R_{t}$, we use all
columns of $\tilde{W}_{t}$ and $W_{t}$, and add in few columns truncated in
Equation 2 in previous iterations with random perturbation when no truncation
happens in Equation 2 and the rank of $\text{orth}(R_{t})$ stops increasing,
so we can ensure that eventually the approximate eigendecomposition approaches
the exact one and $\epsilon_{t}\rightarrow 0$ in Equation 3.

The MF stage requires an initial point $(\tilde{W}_{t},\tilde{H}_{t})$ from
$X_{t}$, and we know that for any $\tilde{W}_{t}\tilde{H}_{t}^{\top}=X_{t}$,
the value of $f$ remains the same, so the key is the part of the norms. The
following lemma justifies our usage of Equation MC as the convex lifting for
Equation MF, and provides the optimal $(\tilde{W}_{t},\tilde{H}_{t})$. The
first part of it is directly from [42].

######  Lemma 2.1.

Given any $X\in\mathbb{R}^{m\times n}$, we have

${\left\|{X}\right\|}_{*}=\min_{W,H:WH^{\top}=X}\quad\frac{1}{2}\left({\left\|{W}\right\|}_{F}^{2}+{\left\|{H}\right\|}_{F}^{2}\right).$
(8)

Moreover, if the SVD of $X$ is
$X=USV^{\top}=\sum_{i=1}^{k}\sigma_{i}u_{i}v_{i}^{\top}$, where
$S=\mathrm{diag}(\sigma_{1},\dots,\sigma_{k})$ and $\sigma_{i}>0$ are the
singular values, $k=\mbox{\rm rank}(X)$,
$U=[u_{1},\dotsc,u_{k}]\in\mathbb{R}^{m\times k}$ and
$V=[v_{1},\dotsc,v_{k}]\in\mathbb{R}^{n\times k}$ are both orthonormal, the
minima of Equation 8 are exactly those

$\hat{W}\coloneqq\left[\sqrt{\sigma_{\tau(1)}}u_{\tau(1)},\dotsc,\sqrt{\sigma_{\tau(k)}}u_{\tau(k)}\right],\hat{H}\coloneqq\left[\sqrt{\sigma_{\tau(1)}}v_{\tau(1)},\dotsc,\sqrt{\sigma_{\tau(k)}}v_{\tau(k)}\right],$
(9)

where $\tau$ is any random permutation of $\\{1,\dots,k\\}$. Therefore, for
any global optimum $(W^{*},H^{*})$ of Equation MF, $X^{*}\coloneqq
W^{*}\left(H^{*}\right)^{\top}$ is also a global optimum to Equation MC,
provided that there is a global optimum $\hat{X}$ of Equation MC with
$\mbox{\rm rank}(\hat{X})\leq k$, and for any optimal solution $X^{*}$ of
Equation MC with SVD $X^{*}=U^{*}S^{*}(V^{*})^{\top}$,
$W^{*}=U^{*}(S^{*})^{1/2}$ and $H^{*}=V^{*}(S^{*})^{1/2}$, form a global
solution of Equation MF for any $k\geq\mbox{\rm rank}(X^{*})$.

The flow of our MF-Global is summarized in Algorithm 2.

1

2for  _$i=0,1,\dotsc$_ do

3 $\alpha_{t}\leftarrow\alpha_{t}^{\text{BB}}\beta^{i}$, and compute an
orthonormal matrix $\hat{U}_{t}\in\mathbb{R}^{m\times k_{t}}$ from a
top-$k_{t}$ approximate eigendecomposition of $Z_{t}Z_{t}^{\top}$ as in
Equation 6, for $Z_{t}$ defined in Equation 3 with the current $\alpha_{t}$,
using $R_{t}$ as the initial point for $\hat{U}_{t}$ and let
$k_{t}\coloneqq\mbox{\rm rank}(R_{t})$.

4 Compute the exact SVD of $\hat{U}_{t}^{\top}Z_{t}$ to get
$\tilde{U}_{t},\hat{\Sigma}_{t},V_{t}$, and then get $\Sigma_{t}$ from
$\hat{\Sigma}_{t}$ as in Equation 2.

5 if  _Equation 5 holds_ then Terminate

6

Algorithm 1 $(U_{t},\Sigma_{t},V_{t})=$ Inexact Proximal Gradient
$(\lambda,\alpha_{t}^{\text{BB}},\beta,\tilde{X}_{t},R_{t})$

input : $\lambda>0,{\alpha_{\max}\geq\alpha_{\min}>0}$, initial rank $k$,
nonnegative sequences $\\{\epsilon_{t}\\}$ with $\epsilon_{t}\rightarrow 0$

1

2Initialize $\tilde{W}_{0}\in\mathbb{R}^{m\times
k},\tilde{H}_{0}\in\mathbb{R}^{n\times k}$.

3for  _$t=0,\dotsc$_ do

4 Compute $(W_{t},H_{t})$ as an approximate solution to Equation MF with rank
$k$ (starting from $(\tilde{W}_{t},\tilde{H}_{t})$) satisfying
$F(W_{t},H_{t})\leq F(X_{t})$.

5 Compute $\alpha^{\text{BB}}_{t}$ by Equation 4 and obtain
$U_{t},\Sigma_{t},V_{t}$ for the SVD of $X_{t+1}$ by Algorithm 1 that
satisfies Equation 3, with $\tilde{X}_{t}=W_{t}H_{t}^{\top}$ and $R_{t}$
initialized using $\tilde{W}_{t}$, $W_{t}$ with small random perturbations.

6 $k\leftarrow\mbox{\rm rank}(\mbox{\rm diag}(\Sigma^{t}))$,
$\tilde{W}_{t+1}\leftarrow U_{t}\mbox{\rm diag}(\sqrt{\Sigma^{t}})$,
$\tilde{H}_{t+1}\leftarrow V_{t}\mbox{\rm diag}(\sqrt{\Sigma^{t}})$.

Algorithm 2 Algorithm Framework of MF-Global for Equation MF

##  3 Analysis

we first discuss the global convergence of MF-Global associated with
$\\{\epsilon_{t}\\}$ in Equation 3. Rank identification of Algorithm 2 is then
proven under a nondegeneracy condition, which shows that for any subsequence
$\\{X_{t_{i}}\\}_{i=0}^{\infty}$ of the iterates that converge to a solution
$X^{*}$, $\mbox{\rm rank}(X_{t_{i}})=\mbox{\rm rank}(X^{*})$ for all $i$ large
enough, so the optimal rank for Equation MF can be automatically decided on
the fly.

For analyzing the worst-case convergence, we define $\Omega^{*}$ as the
solution set to Equation MC and the optimal objective as $F^{*}$. The
following result establishes the global convergence of MF-Global.

######  Theorem 3.1.

Consider Equation MC with $f$ lower-bounded. Then $\Omega^{*}$ is compact and
$F$ is coercive. Moreover, let $\mbox{\rm
dist}(X,\Omega^{*})\coloneqq\min_{X^{*}\in\Omega^{*}}{\left\|{X-X^{*}}\right\|}$,
if $\nabla f$ is Lipschitz continuous and $\sum\epsilon_{t}^{2}<\infty$ in
Algorithm 2, then for any initialization of $\tilde{W}_{0},\tilde{H}_{0}$, we
always have $\mbox{\rm dist}(X_{t},\Omega^{*})\rightarrow 0$, there is at
least one limit point of $\\{X_{t}\\}$, any such limit point is a global
solution to Equation MC, and $F(X_{t})\rightarrow F^{*}$.

Although our global convergence is guaranteed by the inexact PG step, existing
analyses for inexact PG that utilizes the geometry of the iterates like those
in [15, 47, 24, 21] are not applicable, due to that the additional MF phase
could move the iterates arbitrarily in the level set. Therefore, another
contribution of this work is developing new proof techniques for obtaining
global convergence guarantee for general nonmonotone inexact PG combined with
other optimization steps.

Inspired by the proof techniques of [26], under the premise of convergence, we
show by the idea of manifold identification from partial smoothness [29, 22,
31] of the nuclear norm [16] that $\mbox{\rm rank}(X_{t})$ will eventually
become fixed even in the presence of the inexactness from the approximate SVD.

######  Theorem 3.2.

Consider a sequence of iterates $\\{X_{t}\\}$ generated by Algorithm 2 from
some starting point $X_{0}$ with $\epsilon_{t}\rightarrow 0$ in Equation 3.
Then the following hold.

  1. 1.

For a subsequence $\\{\tilde{X}_{t_{i}}\\}_{i}$ such that
$\tilde{X}_{t_{i}}\rightarrow X^{*}$ for some $X^{*}\in\Omega^{*}$,
$X_{t_{i}+1}\rightarrow X^{*}$ as well.

  2. 2.

For the same subsequence in the previous item, if $X^{*}$ satisfies the
nondegeneracy condition

$0\in\mbox{\rm relative interior}\left(\partial F\left(X^{*}\right)\right),$
(10)

then there is $i_{0}\geq 0$ such that $\mbox{\rm rank}(X_{t_{i}+1})=\mbox{\rm
rank}(X^{*})$ for all $i\geq i_{0}$.

Theorem 3.2 implies that when the full iterate sequence converge to a point
$X^{*}$, we get $k_{t}{\coloneqq\mbox{\rm rank}(X_{t})}=\mbox{\rm
rank}(X^{*})$ for all $t$ large enough, and thus through the help of Lemma
2.1, the optimal rank for Equation MF is automatically found. This also shows
that the convex lifting step can effectively make the iterates escape from
saddle points and spurious local minima resulted from an insufficient rank.

##  4 Related works

Methods for Equation MF. State-of-the-art practically efficient solvers for
Equation MF include (semi-asynchronous) parallel stochastic gradient (SG) [14]
and parallel block-coordinate descent (BCD) [57, 51] that fully exploit modern
multicore environments. The SG approach has good performance in terms of the
prediction criterion when suitable parameters are used, but as an optimization
approach it is not as efficient as the BCD approach. Moreover, it requires
more overall running time for tuning additional parameters. We thus use
polyMF-SS of [51] that conducts exact line search for the update step obtained
by BCD as our solver for Equation MF, and also as the baseline for comparison,
because it is shown to outperform [57]. Another popular approach for Equation
MF is Riemannian optimization, including [50, 37, 8, 59, 18], that explicitly
enforces a fixed-rank constraint on $X$ through retractions. The main drawback
of the Riemannian approaches is, just like other methods for Equation MF, that
they have no convergence guarantee to a local or global optimum. In addition,
the retractions could be more expensive, so it is also often observed in
practice that the Riemannian approaches are empirically slower than methods
directly working on $(W,H)$.

Methods for Equation MC. The convex problem Equation MC falls in the category
of regularized optimization, and many efficient algorithms in the optimization
literature are available. However, most methods for regularized optimization
focus on the scenario that the proximal operation can be conducted efficiently
but obtaining information of the smooth term is the major computational
bottleneck, which is apparently not the case for Equation MC. Practical
methods specifically designed for Equation MC all take into serious account
the expensive SVD involved in the nuclear norm [49, 23, 55], and they all
focused on the most popular setting that $f$ is a quadratic term. In this
case, high-order methods like proximal (quasi-)Newton is useless because the
subproblem has the same form as the original problem itself. Therefore, these
methods all consider first-order methods, including inexact PG and inexact
accelerated PG (APG) [38, 5, 6]. [49] used the Lanczos method to conduct
approximate SVD, and [23] proposed to apply the power method for approximate
SVD and to use the rank $k$ decided by their inexact PG to conduct another
convex optimization step with respect to a subproblem of dimension $k\times k$
after each PG step. The power method in [23] effectively uses the current
iterate as warmstart and is much more efficient than the Lanczos approach in
[49], but the additional convex optimization step turns out to be time-
consuming and hard to parallelize. To improve the efficiency of [49] and [23],
Yao et al. [55] combined the two approaches to conduct inexact APG using the
power method, and significantly improved the running time upon both, but the
part on the power method becomes the new bottleneck of [55]. Although the
major computation of matrix-matrix products at each iteration of the power
method is embarrassingly parallel, they are memory-heavy but computation-light
operations, so its parallelism is confined by the memory bandwidth and thus
quite low. Different from [23, 55], for approximate SVDs we use the limited-
memory Krylov subspace method [36] that also utilizes the warmstart from the
current iterate, but has faster convergence and better parallelism than the
power method because it balances computation-heavy and memory-heavy operations
better, so our PG step can also be conducted more efficiently.

More recently, motivated by PG’s ability of manifold identification, [2]
proposed to alternate between an exact PG step and a Riemannian (truncated)
Newton step on the currently identified manifold for general regularized
optimization, and applied this method to a toy problem of Equation MC in their
experiments. Different from ours, their usage of manifold identification is
for showing that their method could obtain superlinear convergence, but their
algorithm is not feasible for large-scale problems considered in this work
because they considered exact PG only and require the explicit computation of
$X_{t}$.

Convergence of Inexact Methods The global convergence of our method is
achieved through the safeguard of inexact PG steps, but we do not require the
inexact PG step to always decrease the objective value. This feature combines
with the MF phase makes the analysis difficult. Existing analyses for inexact
PG [15, 47, 24, 21] utilize telescope sums of inequalities in the form of
${\left\|{X_{t+1}-X^{*}}\right\|}^{2}\leq\text{error
term}+{\left\|{X_{t}-X^{*}}\right\|}^{2}-\alpha\left(F(X_{t})-F(X^{*})\right)$
for any $X^{*}\in\Omega^{*}$ and some $\alpha>0$ to prove convergence and
rates. Therefore, those approaches cannot allow for alternating between
inexact PG and other steps because that will nullify the technique of
telescope sums. On the other hand, analyses compatible with other steps like
those in [46, 7, 27] require strict decreasing of the objective in the inexact
PG step, which imposes an additional burden.

The work of [54] that applies inexact PG to the SDP relaxation of polynomial
optimization problems is probably the closest to ours in that they allow an
alternative step and nonmonotone inexact PG. However, their alternative step
is only accepted when it decreases the objective by an absolute amount of some
$\epsilon>0$, which makes the alternative steps eventually always rejected
when the objective converges to the global optimum, but we do not have such
restrictions. One possible criticism for our analysis is the lack of a rate.
We emphasize that complexity analysis is easily attainable if additional
assumptions are imposed, but those worst-case rates anyway will be pessimistic
and faraway from the fast empirical convergence we observed, and we thus
prefer the current result with novel techniques.

Methods for escaping saddle points. There has recently been a thriving
interest in studying smooth optimization methods that can escape strict saddle
points with at least one negative eigenvalue in the Hessian [28, 25, 44, 12,
45, 1]. However, these methods are unable to deal with degenerate saddle
points where the smallest eigenvalue of the Hessian is exactly zero, and nor
could they deal with spurious local minima that might be arbitrarily far away
from the global ones. On the other hand, our method can handle all such
difficult cases appearing in Equation MF by resorting to convex lifting, and
Lemma 2.1 together with Theorem 3.1 show that our method indeed will converge
to the global optima, so it is more desirable for the matrix factorization
case. Moreover, these methods for escaping strict saddles are mainly of
theoretical interest but their empirical performance is not very impressive,
while our method is designed for practical large-scale usage and greatly
outperforms state-of-the-art methods for Equation MC and Equation MF by a
large margin on real-world data.

Absence of spurious stationary points for Equation MF. To cope with the
possible spurious local minima and degenerate saddle points of matrix
factorization, there are also growing interest in analyzing its optimization
landscape. Most such works consider the quadratic loss defined below, where
given an index set $\Omega$ of observed entries, $(P_{\Omega}(X))_{i,j}$
returns $X_{i,j}$ when $(i,j)\in\Omega$ and $0$ otherwise.

$f(X)={\left\|{P_{\Omega}\left(A-X\right)}\right\|}_{F}^{2}.$ (11)

[19, 48, 20] confined their analyses to Equation 11 and discussed only the
ideal case in which each $(i,j)$ belongs to $\Omega$ with a fixed probability
$p\in(0,1]$ and the observations are noiseless, and showed the absence of
spurious local minima, but these assumptions generally fail in practice. In
particular, for applications like recommendation systems, elements of $\Omega$
are already biased selections by an existing system and will never obey the
independent random assumption, and real-world data always contain noisy
observations and measurement corruption. [13, 56] focused on Equation 11 with
$\Omega=\\{1,\dotsc,m\\}\times\\{1,\dotsc,n\\}$ only, and their techniques are
not generalizable to other problems. Moreover, these two works discussed only
gradient descent, which is barely used (if any) in large-scale real
applications. Different from these works, we do not have any assumption on the
underlying data, as it has been shown recently in [53] that some problem
classes of matrix factorization indeed contain numerous spurious local minima,
yet our method can still find the global optima in the presence of non-strict
saddle points and spurious local optima with ultra-high practical efficiency.

##  5 Experiments

We conduct three sets of experiments to evaluate our algorithm. In particular,
we examine the rank of $X_{t}$ and compare MF-Global with solvers for Equation
MF and for Equation MC. We focus on the case of Equation 11 for its popularity
and use four public large-scale recommendation system data sets and a server
with eight cores. All experiments utilize all available cores for all
algorithms. Detailed data statistics and experiment settings are in Section
A.3. To compare different methods, we consider two criteria, one is the
relative objective value from the optimization perspective, and the other is
the relative RMSE from the task-oriented angle. Their technical definitions
are Section A.3. Here we only mention that the smaller is the better for both,
while their values at an optimum should both be $0$.

###  5.1 Stabilization of the rank

We first demonstrate numerically that MF-Global is able to identify the
optimal rank quickly. We check the rank of $X_{t}$ in Algorithm 2 over time
and the results are in Figure 1, where we use solid lines and dash lines to
respectively show the relative objective value and the rank of $X_{t}$.
Moreover, the gray line represents the rank at the optimum. We can see that
the rank of $X_{t}$ increases quickly at first, and then soon stabilizes at
the rank of the optimal solution in all cases.

|  |  |
---|---|---|---

(a) movielens100k

(b) movielens10m

(c) netflix

(d) yahoo-music

Figure 1: Rank and relative objective of the iterates of MF-Global over
running time.

###  5.2 Comparison with solving Equation MF only

We next compare our method with one that solves Equation MF only. The purpose
of this experiment is to show that solving Equation MF can indeed get stuck at
spurious stationary points, while MF-Global can efficiently escape from such
points. Note that using the relative objective alone is sufficient for this
comparison. To solve Equation MF, we directly use the solver in our MF phase,
namely, polyMF-SS of [51], with their original random initialization to avoid
starting from the origin, which is a known saddle point. We favor their method
to directly assign $k$ in Equation MF as the final rank of an optimal
solution, but we emphasize that in real-world applications, finding this $k$
will require additional efforts in parameter search. Comparisons of polyMF-SS
and MF-Global are shown in Figure 2.

We observe that in both iterations and running time, polyMF-SS is fast at
first, but then its convergence quickly slows down, suggesting that likely the
iterates are attracted to a saddle point or a spurious local minimum that is
strictly worse than the global optima in terms of objective function value. On
the other hand, MF-Global could be slightly slower at the early stage
sometimes, but it maintains a steady convergence towards the global optimum
throughout and always outperforms polyMF-SS eventually. Another observation is
that although we did not prove a convergence rate for Algorithm 2, in the
numerical experiments, the observed convergence speed is nearly linear.

Overall speaking, MF-Global is as efficient as running a solver for Equation
MF alone, while it provides additional advantages like rank identification and
guaranteed convergence to global optima.

movielens100k | movielens10m | netflix | yahoo-music
---|---|---|---
|  |  |
|  |  |

Figure 2: Comparison between polyMF-SS and MF-Global using the relative
objective. Top row: iterations. Bottom row: time (seconds).

###  5.3 Comparison with methods for Equation MC

Now that it is clear MF-Global is advantageous over solvers for Equation MF
alone, we proceed to compare it with the state-of-the-art solvers for Equation
MC. In particular, we compare it with the following:

  * •

Active-ALT [23]: This method alternates between conducting inexact PG and
solving a lower-dimensional convex subproblem. In the approximate SVD, they
use the power method with warmstart input from the output of the previous
iteration plus random columns as a safeguard.

  * •

AIS-Impute [55]: An inexact APG method that also uses the power method for
approximate SVD. They use the combination of the outputs of the previous
iteration and the iteration preceding it to warmstart the power method.

The inexact APG method in [49] is not included because the underlying APG part
is the same as that of AIS-Impute, but their approximate SVD using Lanczos is
less efficient, as shown by [55].

The results are shown in Figure 3. Clearly, MF-Global outperforms the other
two for Equation MC significantly on both criteria. Particularly, our results
exemplify that for reaching a satisfactory RMSE, MF-Global is magnitudes
faster than state-of-the-art solvers for Equation MC.

movielens100k | movielens10m | netflix | yahoo-music
---|---|---|---
|  |  |
|  |  |

Figure 3: Comparison between MF-Global and methods for Equation MC in running
time (seconds). Top row: relative objective. Bottom row: relative RMSE (time
in log scale).

##  6 Conclusions

In this work, we propose an efficient algorithm for finding global optima of
the nonconvex matrix factorization problem. Our method conducts aperiodic
convex lifting steps to effectively escape spurious stationary points of the
nonconvex problems, and exhibits extremely fast convergence in numerical
experiments. The proposed approach pushes forward the state of the art by
being magnitudes faster than existing methods, especially on large-scale
problems. We hope this work will inspire more research on safeguarding fast
nonconvex optimization approaches through convex lifting as well as
accelerating solvers for convex optimization through nonconvex steps. Based on
the results of this work, we will release an open-source package for the
matrix factorization problem soon.

## References

  * Agarwal et al. [2017] Naman Agarwal, Zeyuan Allen-Zhu, Brian Bullins, Elad Hazan, and Tengyu Ma.  Finding approximate local minima faster than gradient descent.  In _Proceedings of the 49th Annual ACM SIGACT Symposium on Theory of Computing_ , pages 1195–1199, 2017.
  * Bareilles et al. [2020] Gilles Bareilles, Franck Iutzeler, and Jérôme Malick.  Newton acceleration on manifolds identified by proximal-gradient methods.  Technical report, 2020.  arXiv:2012.12936.
  * Barzilai and Borwein [1988] Jonathan Barzilai and Jonathan M Borwein.  Two-point step size gradient methods.  _IMA journal of numerical analysis_ , 8(1):141–148, 1988.
  * Bauschke and Combettes [2017] Heinz H. Bauschke and Patrick L. Combettes.  _Convex analysis and monotone operator theory in Hilbert spaces_.  Springer, second edition, 2017.
  * Beck and Teboulle [2009a] Amir Beck and Marc Teboulle.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  _SIAM Journal on Imaging Sciences_ , 2(1):183–202, 2009a.
  * Beck and Teboulle [2009b] Amir Beck and Marc Teboulle.  Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems.  _IEEE transactions on image processing_ , 18(11):2419–2434, 2009b.
  * Bonettini et al. [2016] Silvia Bonettini, Ignace Loris, Federica Porta, and Marco Prato.  Variable metric inexact line-search-based methods for nonsmooth optimization.  _SIAM Journal on Optimization_ , 26(2):891–921, 2016.
  * Cambier and Absil [2016] Léopold Cambier and P-A Absil.  Robust low-rank matrix completion by riemannian optimization.  _SIAM Journal on Scientific Computing_ , 38:S440–S460, 2016\.
  * Candes and Plan [2010] Emmanuel J Candes and Yaniv Plan.  Matrix completion with noise.  _Proceedings of the IEEE_ , 98(6):925–936, 2010\.
  * Candès and Recht [2009] Emmanuel J. Candès and Benjamin Recht.  Exact matrix completion via convex optimization.  _Foundations of Computational mathematics_ , 9(6):717–772, 2009.
  * Candès et al. [2015] Yonina C. Candès, Emmanuel J.and Eldar, Thomas Strohmer, and Vladislav Voroninski.  Phase retrieval via matrix completion.  _SIAM review_ , 57(2):225–251, 2015.
  * Carmon et al. [2018] Yair Carmon, John C. Duchi, Oliver Hinder, and Aaron Sidford.  Accelerated methods for nonconvex optimization.  _SIAM Journal on Optimization_ , 28(2):1751–1772, 2018.
  * Chen et al. [2019] Yuxin Chen, Yuejie Chi, Jianqing Fan, and Cong Ma.  Gradient descent with random initialization: Fast global convergence for nonconvex phase retrieval.  _Mathematical Programming_ , 176(1):5–37, 2019\.
  * Chin et al. [2015] Wei-Sheng Chin, Yong Zhuang, Yu-Chin Juan, and Chih-Jen Lin.  A fast parallel stochastic gradient method for matrix factorization in shared memory systems.  _ACM Transactions on Intelligent Systems and Technology_ , 6:2:1–2:24, 2015.
  * Combettes [2004] Patrick L Combettes.  Solving monotone inclusions via compositions of nonexpansive averaged operators.  _Optimization_ , 53(5-6):475–504, 2004.
  * Daniilidis et al. [2014] Aris Daniilidis, Dmitriy Drusvyatskiy, and Adrian S. Lewis.  Orthogonal invariance and identifiability.  _SIAM Journal on Matrix Analysis and Applications_ , 35(2):580–598, 2014.
  * Ding et al. [2022] Lijun Ding, Dmitriy Drusvyatskiy, and Maryam Fazel.  Flat minima generalize for low-rank matrix recovery.  Technical report, 2022.  arXiv:2203.03756.
  * Gao and Absil [2022] Bin Gao and P-A Absil.  A riemannian rank-adaptive method for low-rank matrix completion.  _Computational Optimization and Applications_ , 81:67–90, 2022.
  * Ge et al. [2016] Rong Ge, Jason D. Lee, and Tengyu Ma.  Matrix completion has no spurious local minimum.  In _Proceedings of the 30th International Conference on Neural Information Processing Systems_ , pages 2981–2989, 2016.
  * Ge et al. [2017] Rong Ge, Chi Jin, and Yi Zheng.  No spurious local minima in nonconvex low rank problems: A unified geometric analysis.  In _International Conference on Machine Learning_ , pages 1233–1242, 2017.
  * Hamadouche et al. [2022] Anis Hamadouche, Yun Wu, Andrew M. Wallace, and Joao F. C. Mota.  Sharper bounds for proximal gradient algorithms with errors.  Technical report, 2022.  arXiv:2203.02204.
  * Hare and Lewis [2004] Warren L. Hare and Adrian S. Lewis.  Identifying active constraints via partial smoothness and prox-regularity.  _Journal of Convex Analysis_ , 11(2):251–266, 2004.
  * Hsieh and Olsen [2014] Cho-Jui Hsieh and Peder Olsen.  Nuclear norm minimization via active subspace selection.  In _International Conference on Machine Learning_ , pages 575–583, 2014.
  * Jiang et al. [2012] Kaifeng Jiang, Defeng Sun, and Kim-Chuan Toh.  An inexact accelerated proximal gradient method for large scale linearly constrained convex SDP.  _SIAM Journal on Optimization_ , 22(3):1042–1064, 2012.
  * Jin et al. [2017] Chi Jin, Rong Ge, Praneeth Netrapalli, Sham M. Kakade, and Michael I. Jordan.  How to escape saddle points efficiently.  In _International Conference on Machine Learning_ , pages 1724–1732, 2017.
  * Lee [2020] Ching-pei Lee.  Accelerating inexact successive quadratic approximation for regularized optimization through manifold identification.  Technical report, 2020.
  * Lee and Wright [2019] Ching-pei Lee and Stephen J. Wright.  Inexact successive quadratic approximation for regularized optimization.  _Computational Optimization and Applications_ , 72:641–674, 2019.
  * Lee et al. [2019] Jason D. Lee, Ioannis Panageas, Georgios Piliouras, Max Simchowitz, Michael I. Jordan, and Benjamin Recht.  First-order methods almost always avoid strict saddle points.  _Mathematical programming_ , 176(1):311–337, 2019\.
  * Lewis [2002] Adrian S. Lewis.  Active sets, nonsmoothness, and sensitivity.  _SIAM Journal on Optimization_ , 13(3):702–725, 2002.
  * Lewis and Overton [1996] Adrian S. Lewis and Michael L. Overton.  Eigenvalue optimization.  _Acta numerica_ , 5:149–190, 1996.
  * Lewis and Zhang [2013] Adrian S. Lewis and Shanshan Zhang.  Partial smoothness, tilt stability, and generalized Hessians.  _SIAM Journal on Optimization_ , 23(1):74–94, 2013.
  * Li et al. [2019a] Qiuwei Li, Zhihui Zhu, and Gongguo Tang.  The non-convex geometry of low-rank matrix optimization.  _Information and Inference: A Journal of the IMA_ , 8(1):51–96, 2019a.
  * Li et al. [2019b] Xingguo Li, Junwei Lu, Raman Arora, Jarvis Haupt, Han Liu, Zhaoran Wang, and Tuo Zhao.  Symmetry, saddle points, and global optimization landscape of nonconvex matrix factorization.  _IEEE Transactions on Information Theory_ , 65(6):3489–3514, 2019b.
  * Li et al. [2018] Yuanzhi Li, Tengyu Ma, and Hongyang Zhang.  Algorithmic regularization in over-parameterized matrix sensing and neural networks with quadratic activations.  In _Conference On Learning Theory_ , pages 2–47, 2018.
  * Liu et al. [2021] Tianyi Liu, Yan Li, Song Wei, Enlu Zhou, and Tuo Zhao.  Noisy gradient descent converges to flat minima for nonconvex matrix factorization.  In _International Conference on Artificial Intelligence and Statistics_ , pages 1891–1899, 2021.
  * Liu et al. [2013] Xin Liu, Zaiwen Wen, and Yin Zhang.  Limited memory block krylov subspace optimization for computing dominant singular value decompositions.  _SIAM Journal on Scientific Computing_ , 35(3):A1641–A1668, 2013.
  * Mishra et al. [2013] Bamdev Mishra, Gilles Meyer, Francis Bach, and Rodolphe Sepulchre.  Low-rank optimization with trace norm penalty.  _SIAM Journal on Optimization_ , 23:2124–2149, 2013.
  * Nesterov [2013] Yurii Nesterov.  Gradient methods for minimizing composite functions.  _Mathematical Programming_ , 140(1):125–161, 2013\.
  * Pal [2010] Amitangshu Pal.  Localization algorithms in wireless sensor networks: Current approaches and future challenges.  _Netw. Protoc. Algorithms_ , 2(1):45–73, 2010\.
  * Recht [2011] Benjamin Recht.  A simpler approach to matrix completion.  _Journal of Machine Learning Research_ , 12(12), 2011.
  * Recht et al. [2010] Benjamin Recht, Maryam Fazel, and Pablo A. Parrilo.  Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization.  _SIAM review_ , 52(3):471–501, 2010.
  * Rennie and Srebro [2005] Jasson DM Rennie and Nathan Srebro.  Fast maximum margin matrix factorization for collaborative prediction.  In _Proceedings of the 22nd international conference on Machine learning_ , pages 713–719, 2005.
  * Rockafellar and Wets [2009] R. Tyrrell Rockafellar and Roger J.-B. Wets.  _Variational analysis_.  Springer Science & Business Media, 2009.
  * Royer and Wright [2018] Clément W. Royer and Stephen J. Wright.  Complexity analysis of second-order line-search algorithms for smooth nonconvex optimization.  _SIAM Journal on Optimization_ , 28(2):1448–1477, 2018.
  * Royer et al. [2020] Clément W Royer, Michael O’Neill, and Stephen J. Wright.  A Newton-CG algorithm with complexity guarantees for smooth unconstrained optimization.  _Mathematical Programming_ , 180(1):451–488, 2020\.
  * Scheinberg and Tang [2016] Katya Scheinberg and Xiaocheng Tang.  Practical inexact proximal quasi-Newton method with global complexity analysis.  _Mathematical Programming_ , 160(1-2):495–529, 2016.
  * Schmidt et al. [2011] Mark Schmidt, Nicolas Roux, and Francis Bach.  Convergence rates of inexact proximal-gradient methods for convex optimization.  In _Advances in Neural Information Processing Systems_ , pages 1458–1466, 2011.
  * Sun and Luo [2016] Ruoyu Sun and Zhi-Quan Luo.  Guaranteed matrix completion via non-convex factorization.  _IEEE Transactions on Information Theory_ , 62(11):6535–6579, 2016.
  * Toh and Yun [2010] Kim-Chuan Toh and Sangwoon Yun.  An accelerated proximal gradient algorithm for nuclear norm regularized linear least squares problems.  _Pacific Journal of optimization_ , 6(615-640):15, 2010.
  * Vandereycken [2013] Bart Vandereycken.  Low-rank matrix completion by riemannian optimization.  _SIAM Journal on Optimization_ , 23:1214–1236, 2013.
  * Wang et al. [2017] Po-Wei Wang, Chun-Liang Li, and J. Zico Kolter.  Polynomial optimization methods for matrix factorization.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 31, 2017.
  * Wright et al. [2009] Stephen J Wright, Robert D Nowak, and Mário AT Figueiredo.  Sparse reconstruction by separable approximation.  _IEEE Transactions on signal processing_ , 57(7):2479–2493, 2009.
  * Yalcin et al. [2022] Baturalp Yalcin, Haixiang Zhang, Javad Lavaei, and Somayeh Sojoudi.  Factorization approach for low-complexity matrix completion problems: Exponential number of spurious solutions and failure of gradient methods.  In _The 25th International Conference on Artificial Intelligence and Statistics_ , 2022.
  * Yang et al. [2021] Heng Yang, Ling Liang, Luca Carlone, and Kim-Chuan Toh.  An inexact projected gradient method with rounding and lifting by nonlinear programming for solving rank-one semidefinite relaxation of polynomial optimization.  Technical report, 2021.  arXiv:2105.14033.
  * Yao et al. [2018] Quanming Yao, James T. Kwok, Taifeng Wang, and Tie-Yan Liu.  Large-scale low-rank matrix learning with nonconvex regularizers.  _IEEE transactions on pattern analysis and machine intelligence_ , 41(11):2628–2643, 2018.
  * Ye and Du [2021] Tian Ye and Simon S. Du.  Global convergence of gradient descent for asymmetric low-rank matrix factorization.  In _Advances in Neural Information Processing Systems_ , volume 34, 2021.
  * Yu et al. [2014] Hsiang-Fu Yu, Cho-Jui Hsieh, Si Si, and Inderjit S. Dhillon.  Parallel matrix factorization for recommender systems.  _Knowledge and Information Systems_ , 41(3):793–819, 2014.
  * Yue et al. [2019] Man-Chung Yue, Zirui Zhou, and Anthony Man-Cho So.  A family of inexact SQA methods for non-smooth convex minimization with provable convergence guarantees based on the Luo–Tseng error bound property.  174(1-2):327–358, 2019.
  * Zhou et al. [2016] Guifang Zhou, Wen Huang, Kyle A Gallivan, Paul Van Dooren, and Pierre-Antoine Absil.  A riemannian rank-adaptive method for low-rank optimization.  _Neurocomputing_ , 192:72–80, 2016.

##  Appendix A Implementation and Experiment Details

This section supplies details in our implementation and the experiments not
described in the main paper.

###  A.1 Implementation details

We now provide implementation details of Algorithm 2.

Initialization. For initializing, it is known that the origin is a saddle
point for Equation MF with Equation 11. Therefore, we start from a point near
the origin with a pre-specified rank $k_{0}$ for $\tilde{W}_{0}$ and
$\tilde{H}_{0}$.

Rank of $R_{t}$ for the Krylov subspace method. As mentioned in the main
paper, we use the columns of $\tilde{W}_{t}$ and $W_{t}$ plus some small
random perturbation of columns truncated in previous iterations. The
truncation step in Equation 2 attempts to decrease the rank of $X_{t}$, while
the combination of $\tilde{W}_{t}$ and $W_{t}$ gives $R_{t}$ a rank that tends
to be higher than that of $X_{t}$. In particular, if columns of
$\tilde{W}_{t}$ and $W_{t}$ are linearly independent, then the rank of $R_{t}$
is twice of that of $X_{t}$. Thus in the fastest case, the rank of $X_{t}$
doubles at each iteration. Adding back the previously truncated columns to
$R_{t}$ then ensures that the rank of $R_{t}$ is indeed strictly larger than
that of $X_{t}$ even if the linear independence condition fails. When the rank
of $X_{t}$ stabilizes, we eventually maintain the rank of $R_{t}$ to be
$\mbox{\rm rank}(X_{t})+1$ to prevent unnecessary computation.

Parallelization. The parallelization of the Krylov subspace method is purely
that of the matrix operations, and it is directly achieved by the inherent
parallelization of matrix operations in MATLAB, while we do no additional
explicit steps to improve the parallelism. For the MF step, the algorithm is
implemented in C/C++ and parallelization is achieved through OpenMP by working
simultaneously on multiple elements of one row of $W_{t}$ and one row of
$H_{t}$.

Computation of Equation 5 The key difficulty in computing the right-hand side
of Equation 5 is the norm term when $X_{t}$ are not formed explicitly. This is
achieved through the observation that

${\left\|{X_{t+1}-X_{t}}\right\|}^{2}=\langle
X_{t+1}-X_{t},\,X_{t+1}-X_{t}\rangle=\langle X_{t+1},\,X_{t+1}\rangle-2\langle
X_{t+1},\,X_{t}\rangle+\langle X_{t},\,X_{t}\rangle,$

and the inner products can be computed efficiently through
$U_{t},\Sigma_{t},V_{t}$.

###  A.2 Selection of $R_{t}$

For the approximate eigendecomposition in Equation 6, the key to the
efficiency and asymptotic convergence to the exact eigendecomposition of both
the power method and the limited-memory Krylov space method is to have a
proper input $R_{t}$ for warmstarting. It is well-known that the optimality
condition for Equation MC is

$X=\mbox{\rm
prox}_{\alpha\lambda{\left\|{\cdot}\right\|}_{*}}\left(X-\alpha\nabla
f(X)\right)$ (12)

for any $\alpha>0$. Therefore, let $Z\coloneqq X-\alpha\nabla f(X)$, for
Equation 2, Equation 12 indicates that the singular vectors of $X$ are exactly
those of $Z$ corresponding to the singular values not truncated to $0$. Hence,
for a sequence $\\{X_{t}\\}$ that approaches the solution set of Equation MC,
we can hope that the singular vectors of $X_{t}$ will also eventually
stabilize. Through this intuition, we take the singular vectors of $X_{t}$ to
form a part of the warmstart input to our approximate eigendecomposition
procedure in computing Equation 6. However, if we only use such singular
vectors, the rank of $X_{t}$ will be monotonically decreasing, in which case
exact eigendecomposition will never be recovered. Therefore, we will need a
mechanism to increase the rank of $X_{t}$. As noted in Lemma 2.1, we know that
the output of approximately solving Equation MF should be close to the
singular vectors of $\tilde{X}_{t}$ (up to column-wise scaling). Therefore, we
will use $U_{t}$ from Equation 7 and $W_{t}$ from the output of the MF phase
to form

$\hat{R}_{t}\coloneqq\text{orth}\left(\left[U_{t},W_{t}\right]\right)$ (13)

to construct the warmstart input to the approximate eigendecomposition in
obtaining $X_{t+1}$ from $Z_{t}$. To further guarantee that
$\epsilon_{t}\rightarrow 0$ in Equation 3, we need to ensure that the rank of
$R_{t}$ is sufficiently large, and that $R_{t}$ will approach the singular
vectors corresponding to the singular values not truncated out. Ideally, we
hope that the output of our approximate eigendecomposition will be exactly all
the eigenvalues or singular values that are retained nonzero, plus the largest
one that is truncated out in Equation 2 so that we know all nonzero singular
values are found. Therefore, we will add in one column with randomness to
$\hat{R}_{t}$ to form $R_{t}$ whenever the rank of $X_{t}$ and $X_{t-1}$ are
the same, and there is at most one singular truncated out in the inexact PG
step at the $(t-1)$th iteration. Utilizing this idea, we retain the singular
vector $u$ that corresponds to the largest truncated singular value in the
latest iteration that such a truncation took place, and compute its projection
$u_{t}$ to $\text{ker}\left(R_{t}^{\top}\right)$, where $\text{ker}(A)$
denotes the null space of $A$ for any matrix $A$. (It is possible and
acceptable that $u_{t}=0$.) The warmstart input is then formed by

$R_{t}\coloneqq\text{orth}\left(\left[R_{t},u_{t}+\xi_{t}\right]\right),\quad\xi_{t}\in\text{ker}\left([R_{t},u_{t}]^{\top}\right),\quad{\left\|{\xi_{t}}\right\|}\leq\psi_{t},$
(14)

where $\xi_{t}$ is a random vector and $\\{\psi_{t}\\}$ is a sequence such
that $\psi_{t}\downarrow 0$.

###  A.3 Experiment details

Environment. Experiments on the first three datasets are conducted on an
Amazon AWS EC2 c6i.4xlarge instance with an 8-core Intel Xeon Ice Lake
processor and 32GB memory, while for the largest yahoo-music that consumes
more memory, an m6i.4xlarge instance with the same processor and 64GB memory
is used. All algorithms are implemented in MATLAB and C/C++, and our
experiments are run in MATLAB 2022b.

Datasets. The datasets used are listed in Table 1. They are obtained through
the following links.

  * •

movielens100k: https://www.kaggle.com/prajitdatta/movielens-100k-dataset. For
the training/test split, we used the split from ua.

  * •

movielens10m: https://www.kaggle.com/smritisingh1997/movielens-10m-dataset.
For the training/test split, we used the split from ra.

  * •

Netflix: https://www.kaggle.com/netflix-inc/netflix-prize-data.

  * •

Yahoo-musc: the R2 one at
https://webscope.sandbox.yahoo.com/catalog.php?datatype=r.

The only preprocessing we did was to transpose the data matrices when
necessary to conform to our assumption of $m\leq n$. For all data sets, we use
their original training/test split. For $\lambda$, we follow the values
provided in [23] obtained through cross-validation, while the final $k$ is the
rank of a global optimum, obtained by running our algorithm with the given
$\lambda$ till the objective cannot be further improved.

Table 1: Data statistics for matrix completion.  Data set | $m$ | $n$ | $|\Omega|$ | $|\Omega_{\text{test}}|$ | $\lambda$ | final $k$
---|---|---|---|---|---|---
movielens100k | 943 | 1682 | 90570 | 9430 | 15 | 68
movielens10m | 65133 | 71567 | 9301274 | 698780 | 100 | 50
netflix | 17770 | 2649429 | 99072112 | 1408395 | 300 | 68
yahoo-music | 624961 | 1000990 | 252800275 | 4003960 | 10000 | 52

Evaluation criteria. To evaluate the algorithms, we first run MF-Global till
the objective cannot be further improved, and take the obtained output $X^{*}$
as the reference point for the optimal solution. The first criterion we
consider is the relative objective

$\frac{F(X)-F^{*}}{F^{*}}$

for any point $X$, where $F^{*}\coloneqq F(X^{*})$ is the optimal objective
value. The second measure is the relative root mean squared error (RMSE),
which is computed as

$\frac{\text{RMSE}(X)-\text{RMSE}(X^{*})}{\text{RMSE}(0)-\text{RMSE}(X^{*})},$

where given the test set $\Omega_{\text{test}}$, the RMSE at $X$ is computed
by

$\text{RMSE}(X)\coloneqq\sqrt{\frac{{\left\|{P_{\Omega_{\text{test}}}\left(X-A\right)}\right\|}^{2}_{F}}{|\Omega_{\text{test}}|}}.$

Settings of our method. In all experiments, we use the following simple
setting for our method. We start from $k_{0}=8$ in initializing Algorithm 2 to
fully utilize all cores available for better parallelism from the beginning.
For simplicity in computing the speedup, this is also the setting in the
speedup experiment when we actually have fewer cores available so that the
speedup is indeed measured under the same amount of computation. For switching
between the MF phase and the convex lifting phase, we run three epochs of
polyMF-SS at each MF phase, where one epoch is counted as one sweep over the
training data. In our setting for Equation 5, we use $\beta=0.5$ and
$\delta=0.99$ throughout.

Comparison with polyMF-SS. In the comparison with polyMF-SS, we also examined
the behavior in terms of iterations. For polyMF-SS, one iteration is counted
as one sweep over the training data (or one epoch), while for our method, we
count either one convex lifting step or one iteration of polyMF-SS as one
iteration in Figure 2.

##  Appendix B Proofs

###  B.1 Proof of Lemma 2.1

The first result of Lemma 2.1 is well-known, and here we provide a proof for
completeness.

###### Proof.

We first prove Equation 8. Denoted by $v^{*}$ the optimal value of the right-
hand side of Equation 8. By construction, we see that $(\hat{W},\hat{H})$ is a
feasible point of the minimization problem. It follows that

$v^{*}\leq\frac{1}{2}\left(\|{\hat{W}}\|_{F}^{2}+\|{\hat{H}}\|_{F}^{2}\right)={\left\|{X}\right\|}_{*}.$

On the other hand, for any $(W,H)$ satisfying $X=WH^{\top}$, it holds that

${\left\|{X}\right\|}_{*}=\mathrm{tr}(S)=\mathrm{tr}\left(U^{\top}WH^{\top}V\right)\leq{\left\|{U^{\top}W}\right\|}_{F}{\left\|{V^{\top}H}\right\|}_{F}\leq\frac{1}{2}\left({\left\|{W}\right\|}_{F}^{2}+{\left\|{H}\right\|}_{F}^{2}\right).$
(15)

Therefore, we obtain ${\left\|{X}\right\|}_{*}\leq v^{*}$. Then, Equation 8
holds true.

We next show that all optimal solutions for the minimization problem in
Equation 8 take the form as Equation 9. Given any $(W,H)$ satisfying
$X=WH^{\top}$ and
${\left\|{X}\right\|}_{*}=\frac{1}{2}\left({\left\|{W}\right\|}_{F}^{2}+{\left\|{H}\right\|}_{F}^{2}\right)$,
by inequality Equation 15, we conclude that $W^{\top}U$ and $H^{\top}V$ are
linearly dependent and it holds that
${\left\|{W^{\top}U}\right\|}_{F}={\left\|{H^{\top}V}\right\|}_{F}$. As a
consequence, we see that $W^{\top}U=H^{\top}V$. Moreover, it follows that

$S=U^{\top}XV=U^{\top}WH^{\top}V=U^{\top}WW^{\top}U,\quad
S=V^{\top}X^{\top}U=V^{\top}HW^{\top}U=V^{\top}HH^{\top}V.$

Since $S$ is a diagonal matrix, we see that $(W,H)$ has the form in Equation
9, as desired.

Last, we prove that Equation MF and Equation MC are equivalent under the rank
condition. Let $(W^{*},H^{*})$ be an optimal solution of Equation MF, and
denote $X^{*}:=W^{*}(H^{*})^{\top}$ with $\mathrm{rank}(X^{*})=k$. Suppose
that problem Equation MC has an optimal solution $\hat{X}$ satisfying
$\mbox{\rm rank}{(\hat{X})}\leq k$. Then, on the one hand, it holds that
$F(\hat{X})\leq F(X^{*})$. On the other hand,

$\displaystyle F(\hat{X})$
$\displaystyle=f(\hat{X})+\lambda{\left\|{X}\right\|}_{*}$
$\displaystyle=f(\hat{X})+\min_{WH^{\top}=\hat{X}}\frac{1}{2}\left({\left\|{W}\right\|}_{F}^{2}+{\left\|{H}\right\|}_{F}^{2}\right)$
$\displaystyle=\min_{WH^{\top}=\hat{X}}f(WH^{\top})+\frac{1}{2}\left({\left\|{W}\right\|}_{F}^{2}+{\left\|{H}\right\|}_{F}^{2}\right),$

which implies that $F(\hat{X})\geq F(X^{*})$. Therefore,
$F(\hat{X})=F(X^{*})$. Thus, the proof is completed. ∎

###  B.2 Proof of Theorem 3.1

###### Proof.

The coerciveness of $F$ follows directly from the fact that the nuclear norm
is coercive and that $f$ is lower-bounded. This implies that the level sets of
$F$ are bounded, and thus so is $\Omega^{*}$.

From the assumption on $\epsilon_{t}^{2}$, we have

$\infty>C^{2}\coloneqq\sum_{t=0}^{\infty}\epsilon_{t}^{2}.$ (16)

From Equation 3, we have that there exists $g_{t}\in\mathbb{R}^{m\times n}$
such that $\|g_{t}\|\leq\epsilon_{t}$ and

$g_{t}\in\nabla
f\left(\tilde{X}_{t}\right)+\frac{1}{\alpha_{t}}\left(X_{t+1}-\tilde{X}_{t}\right)+\partial\left(\lambda\|X_{t+1}\|_{*}\right).$
(17)

From the termination criterion of backtracking Equation 5 and the convexity of
$\|.\|_{*}$, we have the following two inequalities:

$\displaystyle f(X_{t+1})$ $\displaystyle\leq
f\left(\tilde{X}_{t}\right)+\langle\nabla
f(\tilde{X}_{t}),\,X_{t+1}-\tilde{X}_{t}\rangle+\frac{\delta}{\alpha_{t}}\|X_{t+1}-\tilde{X}_{t}\|^{2},$
(18) $\displaystyle\lambda\|X_{t+1}\|_{*}$
$\displaystyle\leq-\langle\xi,\,\tilde{X}_{t}-X_{t+1}\rangle+\lambda\|\tilde{X}_{t}\|_{*},\quad\forall\xi\in\partial\left(\lambda\|X_{t+1}\|_{*}\right).$
(19)

By combining Equations 17, 18 and 19, we get the following inequality.

$\displaystyle F(X_{t+1})$ $\displaystyle\leq F(\tilde{X}_{t})+\langle
g_{t},\,X_{t+1}-\tilde{X}_{t}\rangle-\frac{1-\delta}{\alpha_{t}}\|X_{t+1}-\tilde{X}_{t}\|^{2}$
$\displaystyle\leq
F(\tilde{X}_{t})+\|g_{t}\|\|X_{t+1}-\tilde{X}_{t}\|-\frac{1-\delta}{\alpha_{\max}}\|X_{t+1}-\tilde{X}_{t}\|^{2}$
$\displaystyle\leq
F(\tilde{X}_{t})+\epsilon_{t}\|X_{t+1}-\tilde{X}_{t}\|-\Gamma\|X_{t+1}-\tilde{X}_{t}\|^{2}$
$\displaystyle\leq
F(X_{t})+\epsilon_{t}\|X_{t+1}-\tilde{X}_{t}\|-\Gamma\|X_{t+1}-\tilde{X}_{t}\|^{2},$
(20)

where $\Gamma\coloneqq(1-\delta)/\alpha_{\max}$, and the last inequality is
from that the MF phase does not increase the objective value. Equation 20
implies that

$\Gamma\|X_{t+1}-\tilde{X}_{t}\|^{2}\leq
F(X_{t})-F(X_{t+1})+\epsilon_{t}\|X_{t+1}-\tilde{X}_{t}\|.$ (21)

By summing Equation 21 from $t=0$ to $t=k,$ we have that

$\displaystyle\Gamma\sum_{t=0}^{k}\|X_{t+1}-\tilde{X}_{t}\|^{2}$
$\displaystyle\leq
F(X_{0})-F(X_{k+1})+\sum_{t=0}^{k}\epsilon_{t}\|X_{t+1}-\tilde{X}_{t}\|$ (22)
$\displaystyle\leq
F(X_{0})-F^{*}+\sum_{t=0}^{k}\epsilon_{t}\|X_{t+1}-\tilde{X}_{t}\|$
$\displaystyle\leq
F(X_{0})-F^{*}+\sqrt{\sum_{t=0}^{k}\epsilon_{t}^{2}}\sqrt{\sum_{t=0}^{k}\|X_{t+1}-\tilde{X}_{t}\|^{2}}$
(23)
$\displaystyle\stackrel{{\scriptstyle\lx@cref{creftype~refnum}{eq:C}}}{{\leq}}F(X_{0})-F^{*}+C\sqrt{\sum_{t=0}^{k}\|X_{t+1}-\tilde{X}_{t}\|^{2}},$
(24)

where Equation 23 is from the Cauchy-Schwarz inequality. By applying the
quadratic formula to Equation 24, we obtain that

$\sqrt{\sum_{t=0}^{k}\|X_{t+1}-\tilde{X}_{t}\|^{2}}\leq\frac{C+\sqrt{C^{2}+4\Gamma\left(F(X_{0})-F^{*}\right)}}{2\Gamma}.$
(25)

This implies that

$\sum_{t=0}^{k}\epsilon_{t}\|X_{t+1}-\tilde{X}_{t}\|\leq\sqrt{\sum_{t=0}^{k}\epsilon_{t}^{2}}\sqrt{\sum_{t=0}^{k}\|X_{t+1}-\tilde{X}_{t}\|^{2}}\leq\frac{C^{2}+C\sqrt{C^{2}+4\Gamma\left(F(X_{0})-F^{*}\right)}}{2\Gamma}.$
(26)

Combining Equation 26 and Equation 22, we have that $\\{F(X_{t})\\}$ is an
upper-bounded sequence. Then from the coerciveness of $F$, the sequence
$\\{X_{t}\\}$ is bounded, and thus it has at least one limit point.

Now, let us assume that $\mbox{\rm dist}(X_{t},\Omega^{*})$ does not approach
zero for contradiction. Then there exists $\sigma>0$ and a subsequence
$\\{X_{t_{k}}\\}_{k}$ such that $\mbox{\rm dist}(X_{t_{k}},\Omega)>\sigma>0$
for any $k$. Since $\\{X_{t_{k}}\\}$ is also a bounded sequence, we have that
there exists a subsubsequence $\\{X_{\ell_{k}}\\}\subseteq\\{X_{t_{k}}\\}$
such that

$X_{\ell_{k}}\rightarrow\widetilde{X^{*}}\notin\Omega^{*}.$ (27)

From Equation 17, we have that

$g_{t}+\nabla f(X_{t+1})-\nabla
f(\tilde{X}_{t})-\frac{1}{\alpha_{t}}(X_{t+1}-\tilde{X}_{t})\in\partial
F(X_{t+1}).$ (28)

From that $\sum_{t=0}^{\infty}\epsilon^{2}_{t}<\infty,$ we have
$\epsilon_{t}\rightarrow 0$ and hence $g_{t}\rightarrow 0.$ From Equation 25,
we have that $X_{t+1}-\tilde{X}_{t}\rightarrow 0.$ This together with the
Lipschitz continuity of $\nabla f$ and the boundedness of $\alpha_{t}$ implies
that $\nabla f(X_{t+1})-\nabla f(\tilde{X}_{t})\rightarrow 0$ and
$\alpha_{t}^{-1}(X_{t+1}-\tilde{X}_{t})\rightarrow 0.$ These results together
imply that

$g_{t}+\nabla f(X_{t+1})-\nabla
f(\tilde{X}_{t})-\frac{1}{\alpha_{t}}(X_{t+1}-\tilde{X}_{t})\rightarrow 0.$
(29)

From Equations 27, 28 and 29 and the outer semi-continuity of $\partial F$
(see [43, Proposition 8.7] and [4, Proposition 20.37]) in Equation 28, we have
that

$0\in\partial F\big{(}\widetilde{X^{*}}\big{)}.$

From the convexity of $F$, we get $\widetilde{X^{*}}\in\Omega^{*}$,
contradicting Equation 27. Therefore, we conclude $\mbox{\rm
dist}(X_{t},\Omega^{*})\rightarrow 0.$ This also implies that any limit point
of $\\{X_{t}\\}$ must lie in $\Omega^{*}$.

From the definition of subdifferential, the convexity of $F$, and Equation 28,
we see that for any $X^{*}\in\Omega^{*}$,

$\displaystyle 0$ $\displaystyle\leq
F\left(X_{t+1}\right)-F\left(X^{*}\right)$ $\displaystyle\leq\langle
g_{t}+\nabla f\left(X_{t+1}\right)-\nabla
f\big{(}\tilde{X}_{t}\big{)}-\frac{1}{\alpha_{t}}\big{(}X_{t+1}-\tilde{X}_{t}\big{)},\,X_{t+1}-X^{*}\rangle$
$\displaystyle\leq{\left\|{g_{t}+\nabla f\big{(}X_{t+1}\big{)}-\nabla
f\big{(}\tilde{X}_{t}\big{)}-\frac{1}{\alpha_{t}}\big{(}X_{t+1}-\tilde{X}_{t}\big{)}}\right\|}{\left\|{X_{t+1}-X^{*}}\right\|}.$
(30)

By the boundedness of $\\{X_{t}\\}$, ${\left\|{X_{t+1}-X^{*}}\right\|}$ is
upper bounded, while the first norm term in Equation 30 approaches to $0$ as
shown in Equation 29. We thus conclude from the sandwich lemma that
$F(X_{t})\rightarrow F^{*}$. ∎

###  B.3 Proof of Theorem 3.2

Before proving Theorem 3.2, we will need to define partly smooth functions
first. As we consider the nuclear norm only, which is convex, we can use the
following simplified definition.

######  Definition B.1 (Partly smooth [29, 22]).

A convex function $\Psi$ is partly smooth at a point $x^{*}$ relative to a set
$\mathcal{M}$ containing $x^{*}$ if $\partial\Psi(x^{*})\neq\emptyset$ and:

  1. 1.

Around $x^{*}$, $\mathcal{M}$ is a $\mathcal{C}^{2}$-manifold and
$\Psi|_{\mathcal{M}}$ is $\mathcal{C}^{2}$.

  2. 2.

The affine span of $\partial\Psi(x)$ is a translate of the normal space to
$\mathcal{M}$ at $x^{*}$.

  3. 3.

$\partial\Psi$ is continuous at $x^{*}$ relative to $\mathcal{M}$.

Loosely speaking, this means that $\Psi|_{\mathcal{M}}$ is smooth at $x^{*}$,
but the value of $\Psi$ changes drastically along directions leaving
$\mathcal{M}$ around $x^{*}$. We will use the partial smoothness of the
nuclear norm to prove Theorem 3.2.

###### Proof of Theorem 3.2.

  1. 1.

Let us denote the exact solution of Equation 1 at the $t_{i}$-th iteration
given $\tilde{X}_{t_{i}}$ as $\hat{X}_{t_{i}+1}$, and the real update we use
from Equation 3 as $X_{t_{i}+1}$. Following the proof of [58, Proposition 1],
we have from the optimality of $X^{*}$, which implies $\mbox{\rm
prox}_{\alpha\lambda{\left\|{\cdot}\right\|}_{*}}(X^{*}-\alpha\nabla
f(X^{*}))=X^{*}$, that

$\displaystyle 0$
$\displaystyle\leq{\left\|{\hat{X}_{t_{i}+1}-\tilde{X}_{t_{i}}}\right\|}$
$\displaystyle={\left\|{\mbox{\rm
prox}_{\alpha_{t_{i}}\lambda{\left\|{\cdot}\right\|}_{*}}\left(\tilde{X}_{t_{i}}-\alpha_{t_{i}}\nabla
f(\tilde{X}_{t_{i}}\right)-\tilde{X}_{t_{i}}+\left(X^{*}-\mbox{\rm
prox}_{\alpha_{t_{i}}\lambda{\left\|{\cdot}\right\|}_{*}}(X^{*}-\alpha_{t_{i}}\nabla
f(X^{*}))\right)}\right\|}$
$\displaystyle\leq{\left\|{\tilde{X}_{t_{i}}-X^{*}}\right\|}+{\left\|{\mbox{\rm
prox}_{\alpha_{t_{i}}\lambda{\left\|{\cdot}\right\|}_{*}}\left(\tilde{X}_{t_{i}}-\alpha_{t_{i}}\nabla
f(\tilde{X}_{t_{i}})\right)-\mbox{\rm
prox}_{\alpha_{t_{i}}\lambda{\left\|{\cdot}\right\|}^{*}}\left(X^{*}-\alpha_{t_{i}}\nabla
f(X^{*})\right)}\right\|}$
$\displaystyle\leq{\left\|{\tilde{X}_{t_{i}}-X^{*}}\right\|}+{\left\|{\left(\tilde{X}_{t_{i}}-\alpha_{t_{i}}\nabla
f(\tilde{X}_{t_{i}})\right)-\left(X^{*}-\alpha_{t_{i}}\nabla
f(X^{*})\right)}\right\|}$ (31) $\displaystyle\leq
2{\left\|{\tilde{X}_{t_{i}}-X^{*}}\right\|}+\alpha_{t_{i}}{\left\|{\nabla
f(\tilde{X}_{t_{i}})-\nabla f(X^{*})}\right\|}$
$\displaystyle\leq(2+L\alpha_{\max}){\left\|{\tilde{X}_{t_{i}}-X^{*}}\right\|}\rightarrow
0,$ (32)

where Equation 31 is from the nonexpansiveness of the proximal operation of
any convex function and Equation 32 is from the upper bound in Equation 4 and
the assignment of $\alpha_{t}$ in Equation 5. Therefore, Equation 32 leads to

${\left\|{\hat{X}_{t_{i}+1}-\tilde{X}_{t_{i}}}\right\|}\rightarrow 0.$ (33)

On the other hand, $Q^{\alpha_{t_{i}}}_{\tilde{X}_{t_{i}}}(\cdot)$ is
$(\alpha_{t_{i}})^{-1}$-strongly convex for all $i$, so the minimizer is
unique and by definition it is $\hat{X}_{t_{i}+1}$. Moreover, through
straightforward computation, we know that such strong convexity of
$Q^{\alpha_{t_{i}}}_{\tilde{X}_{t_{i}}}$ and the condition Equation 3 imply
that

$\displaystyle 0$
$\displaystyle\leq\frac{1}{2\alpha_{t_{i}}}{\left\|{X_{t_{i}+1}-\hat{X}_{t_{i}+1}}\right\|}^{2}$
$\displaystyle\leq
Q^{\alpha_{t_{i}}}_{\tilde{X}_{t_{i}}}\left(X_{t_{i}+1}\right)-Q^{\alpha_{t_{i}}}_{\tilde{X}_{t_{i}}}\left(\hat{X}_{t_{i}+1}\right)$
$\displaystyle\leq\frac{\alpha_{t_{i}}}{2}\min_{g_{t_{i}}\in\partial
Q^{\alpha_{t_{i}}}_{\tilde{X}_{t_{i}}}\left(X_{t_{i}+1}\right)}{\left\|{g_{t_{i}}}\right\|}^{2}$
$\displaystyle\leq\frac{\alpha_{t_{i}}}{2}\epsilon_{t_{i}}^{2}\rightarrow 0.$
(34)

The limit is obtained from that $\epsilon_{t}\rightarrow 0$ and that
$\alpha_{t}$ is upper-bounded by $\alpha_{\max}$ from Equation 4. By combining
Equations 33 and 34, it is clear that

${\left\|{X_{t_{i}+1}-X^{*}}\right\|}\rightarrow
0,\quad{\left\|{X_{t_{i}+1}-\tilde{X}_{t_{i}}}\right\|}\rightarrow 0,$ (35)

and the former proves the desired result.

  2. 2.

It is well-known that $\lambda{\left\|{X}\right\|}_{*}$ is partly smooth at
every $X$ relative to the manifold

$\mathcal{M}_{X}\coloneqq\left\\{Y\in\mathbb{R}^{m\times n}\mid\mbox{\rm
rank}(Y)=\mbox{\rm rank}(X)\right\\},$ (36)

see, for example, [16] or [2, Example 2.4]. Therefore, the result of the
second item of Theorem 3.2 is equivalent to
$X_{t_{i}+1}\in\mathcal{M}_{X^{*}}$ for all $i$ large enough. As
$\epsilon_{t}\rightarrow 0$, we see that all conditions of [26, Theorem 1] are
satisfied, and therefore the conclusion follows.

∎

##  Appendix C Additional Experiments

We present some additional experimental results in this section.

###  C.1 Performance in terms of iterations

We present comparison between different methods in terms of iteration count.
For AIS-Impute and Active-ALT, we use their original defnitions for counting
iterations. For MF-Global, one sweep over the training data in the MF phase is
counted as one iteration, while one convex lifting step is also counted as one
iteration. This result is shown in Figure 4, and although each method might
have different per-iteration costs, we observe a trend similar to that in the
comparison in running time that MF-Global is much faster than approaches for
Equation MC.

movielens100k | movielens10m | netflix | yahoo-music
---|---|---|---
|  |  |
|  |  |

Figure 4: Comparison between MF-Global and methods for Equation MC in
iterations. Top row: relative objective. Bottom row: relative RMSE (iteration
count in log scale).

###  C.2 Speedup

We examine in Figure 5 the speedup of different methods with respect to the
number of cores. The speedup of each method is independently measured by

$\text{Speedup}(x)=\frac{\text{Running time of $x$ cores}}{\text{Running time
of $1$ core}}.$

We skip the small dataset movielens100k in this experiment because no method
achieves observable speedup on this data, as it is too small for the benefits
of multicore parallelization to compensate the additional overhead of
initializing multithread processes.

We can see that methods for Equation MC have barely any speedup at all, so
they are less suitable for large-scale problems. On the other hand, although
MF-Global sacrificed some speedup when compared with polyMF-SS alone, it can
still utilize the multiple cores available in modern computing environments to
have speedup increasing with number of cores used.

|  |
---|---|---

(a) movielens10m

(b) netflix

(c) yahoo-music

Figure 5: Speedup comparison of algorithms.