Finding Low-Rank Solutions via Non-Convex Matrix Factorization, Efficiently and Provably

Introduction

Specific instances of (1), where the solution is assumed low-rank, appear in several applications in diverse research fields; a non-exhaustive list includes factorization-based recommender systems , multi-label classification tasks , dimensionality reduction techniques , density matrix estimation of quantum systems , phase retrieval applications , sensor localization and protein clustering tasks, image processing problems , as well as applications in system theory , just to name a few. Thus, it is critical to devise user-friendly, efficient and provable algorithms that solve (1), taking into consideration such near low-rank structure of $X^{\star}$.

While, in general, imposing a low-rankness may result in an NP-hard problem, (1) with a rank-constraint can be solved in polynomial-time for numerous applications, where $f$ has specific structure. A prime—and by now well-known—example of this is the matrix sensing/matrix completion problem (we discuss this further in Section 1.1). There, $f$ is a least-squares objective function and the measurements satisfy the appropriate restricted isometry/incoherence assumptions. In such a scenario, the optimal low-rank $X^{\star}$ can be recovered in polynomial time, by solving (1) with a rank-constraint , or by solving its convex nuclear-norm relaxation, as in .

In view of the above and although such algorithms have attractive convergence rates, they directly manipulate the $n\times n$ variable matrix $X$, which in itself is computationally expensive in the high-dimensional regime. Specifically, each iteration in these schemes typically requires computing the top-$r$ singular value/vectors of the matrix. As $n$ scales, these computational demands at each iteration can be prohibitive.

Note that characterizations (2) and (1) are equivalent in the case $\text{rank}(X^{\star})=r$.Here, by equivalent, we mean that the set of global minima in (2) contains that of (1). It remains an open question though whether the reformulation in (2) introduces spurious local minima in the factored space for the majority of $f$ cases. Observe that such parameterization leads to a very specific kind of non-convexity in $f$. Even more importantly, proving convergence for these settings becomes a harder task, due to the bi-linearity of the variable space.

Our contributions. While the computational gains are apparent, such bi-linear reformulations $X=UV^{\top}$often lack theoretical guarantees. Only recently, there have been attempts in providing answers to when and why such non-convex approaches perform well in practice, in the hope that they might provide a new algorithmic paradigm for designing faster and better algorithms; see .

As we describe below and in greater detail in Section 1.2, our work is more general, addressing important settings that could not (as far as we know) be treated by the previous literature. Our contributions can be summarized as follows:

We study a gradient descent algorithm on the non-convex formulation given in (2) for non-square matrices. We call this Bi-Factored Gradient Descent (BFGD). Recent developments (cited above, and see Section 1.2 for further details) rely on properties of $f$ for special cases , and their convergence results seem to rely on this special structure. In this work, we take a more generic view of such factorization techniques, closer to results in convex optimization. We provide local convergence guarantees for general smooth (and strongly convex) $f$ objectives.

In particular, when $f$ is only (restricted) smooth, we show that a simple lifting technique leads to a local sublinear rate convergence guarantee, using results from that of the square and PSD case . Moreover, we provide a simpler and improved proof than , which requires a weaker initial condition.

When $f$ is both (restricted) strongly convex and smooth, results from the PSD case do not readily extend. In such cases, of significant importance is the use of a regularizer in the objective, that restricts the geometry of the problem at hand. Here, we improve upon —where such a regularizer was used only for the cases of matrix sensing/completion and robust PCA—and solve a different formulation that lead to local linear rate convergence guarantees. Our proof technique proves a significant generalization: using any smooth and strongly convex regularizer on the term $(U^{\top}U-V^{\top}V)$, with optimum at zero, one can guarantee linear convergence.

Our theory is backed up with extensive experiments, including affine rank minimization (Section 6.3), compressed noisy image reconstruction from a subset of image pixels (Section 6.4), and 1-bit matrix completion tasks (Section 6.5). Overall, our proposed scheme shows superior recovery performance, as compared to state-of-the-art approaches, while being $(i)$ simple to implement, $(ii)$ scalable in practice and, $(iii)$ versatile to various applications.

In this section, we describe some applications that can be modeled as in (2). The list includes criteria with $(i)$ smooth and strongly convex objective $f$ (e.g., quantum state tomography from a limited set of observations and compressed image de-noising), and $(ii)$ just smooth objective $f$ (e.g., 1-bit matrix completion and logistic PCA). For all cases, we succinctly describe the problem and provide useful references on state-of-the-art approaches; we restrict our discussion on first-order, gradient schemes. Some discussion regarding recent developments on factorized approaches is deferred to Section 1.2. Section 6 provides specific configuration of our algorithm, for representative tasks, and make a comparison with state of the art.

Matrix sensing (MS) problems have gained a lot of attention the past two decades, mostly as an extension of Compressed Sensing to matrices; see . The task involves the reconstruction of an unknown and low-rank ground truth matrix $X^{\star}$, from a limited set of measurements. The assumption on low-rankness depends on the application at hand and often is natural: e.g., in background subtraction applications, $X^{\star}$ is a collection of video frames, stacked as columns, where the “action” from frame to frame is assumed negligible ; in (robust) principal component analysis , we intentionally search for a low-rank representation of the data; in linear system identification, the low rank $X^{\star}$ corresponds to a low-order linear, time-invariant system ; in sensor localization, $X^{\star}$ denotes the matrix of pairwise distances with rank-dependence on the, usually, low-dimensional space of the data ; in quantum state tomography, $X^{\star}$ denotes the density state matrix of the quantum system and $X^{\star}$ is designed to be rank-1 (pure state) or rank-$r$ (almost pure state), for $r$ relatively small .

In a non-factored form, MS is expressed via the following criterion:

Critical assumption for $\mathcal{A}$ that renders (3) a polynomially solvable problem, is the restricted isometry property (RIP) for low-rank matrices :

A linear map $\mathcal{A}$ satisfies the $r$-RIP with constant $\delta_{r}$, if

It turns out linear maps that satisfy Definition 1.1 also satisfy the (restricted) smoothness and strong convexity assumptions ; see Theorem 2 in and Section 2 for their definition.

State-of-the-art approaches. The most popularized approach to solve this problem is through convexification: show that the nuclear norm $\|\cdot\|_{*}$ is the tightest convex relaxation of the non-convex $\text{rank}(\cdot)$ constraint and algorithms involving nuclear norm have been shown to be effective in recovering low-rank matrices. This leads to:

Efficient implementations include Augmented Lagrange Multiplier (ALM) methods , convex conic solvers like the TFOCS software package and, convex proximal and projected first-order methods . However, due to the nuclear norm, in most cases these methods are binded with full SVD computations per iteration, which constitutes them impractical in large-scale settings.

From a non-convex perspective, algorithms that solve (3) in a non-factored form include SVP and Randomized SVP algorithms , Riemannian Trust Region Matrix Completion algorithm (RTRMC) , ADMiRA and the Matrix ALPS framework .

In all cases, algorithms admit fast linear convergence rates. Moreover, the majority of approaches assumes a first-order oracle: information of $f$ is provided through its gradient $\nabla f(X)$. For MS, $\nabla f(X)=-2\mathcal{A}^{*}\left(y-\mathcal{A}(X)\right)$, which requires $O({\rm T}_{\text{map}})$ complexity, where ${\rm T}_{\text{map}}$ denotes the time required to apply linear map (or its adjoint $\mathcal{A}^{*}$) $\mathcal{A}$. Formulations (3)-(5) require at least one top-$r$ SVD calculation per iteration; this translates into additional $O(mnr)$ complexity.

Motivation for factorizing (3). Problem (3) can be factorized as follows:

For this case and assuming a first-order oracle, the gradient of $f$ with respect to $U$ and $V$ can be computed respectively as $\nabla_{U}f(UV^{\top}):=\nabla f(X)V$ and $\nabla_{V}f(UV^{\top}):=\nabla f(X)^{\top}U$, respectively. This translates into $2\cdot O({\rm T}_{\text{map}}+mnr)$ time complexity. However, one avoids performing any SVD calculations per iteration, which in practice is considered a great computational bottleneck, even for moderate $r$ values. Thus, if there exist linearly convergent algorithms for (6), intuition indicates that one could obtain computational gains.

1.2 Logistic PCA and Low-Rank Estimation on Binary Data

Finding a low-rank approximation of binary matrices has gain a lot of interest recently, due to the wide appearance of categorical responses in real world applications . While regular linear principal component analysis (PCA) is still applicable for binary or categorical data, $(i)$ the way data are pre-processed (e.g., centering data before applying PCA), and/or $(ii)$ the least-squares objective in PCA, constitute it a natural choice for real-valued data, where observations are assumed to follow a Gaussian distribution. propose generalized versions of PCA for other type of data sets: In the case of binary data, this leads to Logistic Principal Component Analysis (Logistic PCA), where each binary data vector is assumed to follow the multivariate Bernoulli distribution, parametrized by the principal components that live in a $r$-dimensional subspace. Moreover, collaborative filtering on binary data and network sign prediction tasks have shown that standard least-squares loss functions perform poorly, while generic logistic loss optimization shows more interpretable and promising results.

where we assume independence among entries of $Y$. The negative log-likelihood for log-odds parameter $X$ is given by:

Assuming a compact, i.e., low-rank, representation for the latent variable $X$, we end up with the following optimization problem:

observe that the objective criterion is just a smooth convex loss function.

State-of-the-art approaches.Here, we note that proposes a slightly different way to generalize PCA than , based on a different interpretation of Pearson’s PCA formulation . The resulting formulation looks for a weighted projection matrix $UU^{\top}$ (instead of $UV^{\top})$, where the number of parameters does not increase with the number of samples and the application of principal components to new data requires only one matrix multiplication. For this case, the authors in propose, among others, an alternating minimization technique where convergence to a local minimum is guaranteed. Even for this case though, our framework applies. In , the authors consider the problem of sign prediction of edges in a signed network and cast it as a low-rank matrix completion problem: In order to model sign inconsistencies between the entries of binary matrices, the authors consider more appropriate loss functions to minimize, among which is the logistic loss. The proposed algorithmic solution follows (stochastic) gradient descent motions; however, no guarantees are provided. utilizes logistic PCA for collaborative filtering on implicit feedback data (page clicks and views, purchases, etc.): to find a local minimum, an alternating gradient descent procedure is used—further, the authors use AdaGrad to adaptively update the gradient descent step size, in order to reduce the number of iterations for convergence. A similar alternating gradient descent approach is followed in , with no known theoretical guarantees.

Motivation for factorizing (7). Following same arguments as before, in logistic PCA and logistic matrix factorization problems, we often assume that the observation binary matrix is generated as the sign operation on a linear factored model: $\texttt{sign}(UV^{T})$. Parameterized by the latent factors $U,V$, we obtain the following optimization criterion:

where $U_{i}$, $V_{j}$ represent the $i$-th and $j$-th row of $U$ and $V$, respectively.

2 Related Work

As it is apparent from the discussion above, this is not the first time such transformations have been considered in practice. Early works on principal component analysis and non-linear estimation procedures , use this technique as a heuristic; empirical evaluations show that such heuristics work well in practice . further popularized these ideas for solving SDPs: their approach embeds the PSD and linear constraints into the objective and applies low-rank variable re-parameterization. While the constraint considered here is of different nature—i.e, rank constraint vs. PSD constraint—the motivation is similar: in SDPs, by representing the solution as a product of two factor matrices, one can remove the positive semi-definite constraint and thus, avoid computationally expensive projections onto the PSD cone.

We provide an overview of algorithms that solve instances of (2); for discussions on methods that operate on $X$ directly, we defer the reader to for more details; see also Table 1 for an overview of the discussion below. We divide our discussion into two problem settings: $(i)$ $X^{\star}$ is square and PSD and, $(ii)$ $X^{\star}$ is non-square.

Several recent works have studied (9). For the special case where $f$ is a least-squares objective for an underlying linear system, and propose gradient descent schemes that function on the factor $U$. Both studies employ careful initialization (performing few iterations of SVP for the former and, using a spectral initialization procedure—as in —for the latter) and step size selection, in order to prove convergence.Recently, and proved that $UU^{\top}$ factorization introduces no spurious local minima for the cases of matrix completion and sensing, respectively: random initialization eventually leads to convergence to the optimal $X^{\star}$ (or close to $X^{\star}$ in the approximately low rank case). The extension of these results to the non-square matrix sensing setting can be found in . However, their analysis is designed only for least-squares instances of $f$. Some results and discussion on their step size selection/initialization and how it compares with this work are provided in Section 6.

The work of proposes a first-order algorithm for (9), where $f$ is more generic. The algorithmic solution proposed can handle additional constraints on the factors $U$; the nature of these constraints depends on the problem at hand.Any additional constraints should satisfy the faithfulness property: a constraint set $\mathcal{C}$ is faithful if for each $U\in\mathcal{C}$, within some bounded radius from optimal point, we are guaranteed that the closest (in the Euclidean sense) rotation of optimal $U^{\star}$ lies within $\mathcal{U}$. The authors present a broad set of exemplars for $f$—matrix completion and sensing, as well as sparse PCA, among others. For each problem, a set of assumptions need to be satisfied; i.e., faithfulness, local descent, local Lipschitz and local smoothness conditions; see for more details. Under such assumptions and with proper initialization, one can prove convergence with $O(\nicefrac{{1}}{{\varepsilon}})$ or $O(\log(\nicefrac{{1}}{{\varepsilon}}))$ rate, depending on the nature of $f$, and for problems that even fail to be locally convex.

proposes Factored Gradient Descent (FGD) algorithm for (9). FGD is also a first-order scheme; key ingredient for convergence is a novel step size selection that can be used for any $f$, as long as it is (restricted) gradient Lipschitz continuous; when $f$ is further (restricted) strongly convex, their analysis lead to faster convergence rates. Using proper initialization, this is the first paper that provably solves (9) for general convex functions $f$ and under common convex assumptions. An extension of these ideas to some constrained cases can be found in .

To summarize, most of these results guarantee convergence—up to linear rate—on the factored space, starting from a “good” initialization point and employing a carefully selected step size.

Non-square $X^{\star}$. propose AltMinSense, an alternating minimization algorithm for matrix sensing and matrix completion problems. This is one of the first works to prove linear convergence in solving (2) for the MS model. improves upon for the case of reasonably well-conditioned matrices. Their algorithm handles problem cases with bad condition number and gaps in their spectrum . Recently, extended the Procrustes Flow algorithm to the non-square case, where gradient descent, instead of exact alternating minimization, is utilized. extended the first-order method of for matrix completion to the rectangular case. All the studies above focus on the case of least-squares objective $f$.

generalize the results in : the authors show that, under common incoherence conditions and sampling assumptions, most first-order variants (e.g., gradient descent, alternating minimization schemes and stochastic gradient descent, among others) indeed converge to the low-rank ground truth $X^{\star}$. Both the theory and the algorithm proposed are restricted to the matrix completion objective.

Recently, —based on the inexact first-order oracle, previously used in —proved that linear convergence is guaranteed if $f(UV^{\top})$ is strongly convex over either $U$ and $V$, when the other is fixed. While the technique applies for generic $f$ and for non-square $X$, the authors provide algorithmic solutions only for matrix completion / matrix sensing settings.E.g., in the gradient descent case, the step size proposed depends on RIP constants and it is not clear what a good step size would be in other problem settings. Furthermore, their algorithm requires QR-decompositions after each update of $U$ and $V$; this is required in order to control the notion of inexact first order oracle.

Preliminaries

Given a matrix $X$, we denote its best rank-$r$ approximation with $X_{r}$; $X_{r}$ can be computed in polynomial time via the SVD. For our discussions from now on and in an attempt to simplify our notation, we denote the optimum point we search for as $X^{\star}_{r}$, both $(i)$ in the case where we intentionally restrict our search to obtain a rank-$r$ approximation of $X^{\star}$—while $\text{rank}(X^{\star})>r$—and $(i)$ in the case where $X^{\star}\equiv X^{\star}_{r}$, i.e., by default, the optimum point is of rank $r$.

An important issue in optimizing $f$ over the factored space is the existence of non-unique possible factorizations for a given $X$. Since we are interested in obtaining a low-rank solution in the original space, we need a notion of distance to the low-rank solution $X^{\star}_{r}$ over the factors. Among infinitely many possible decompositions of $X^{\star}_{r}$, we focus on the set of “equally-footed” factorizations :

Given a pair $(U,V)$, we define the distance to $X^{\star}_{r}$ as:

Assumptions. We consider applications that can be described $(i)$ either by restricted strongly convex functions $f$ with gradient Lipschitz continuity, or $(ii)$ by convex functions $f$ that have only Lipschitz continuous gradients.Our ideas can be extended in a similar fashion to the case of restricted strong convexity . We state these standard definitions below.

The Factored Gradient Descent (FGD) algorithm. Part of our contributions is inspired by the work of , where the FGD algorithm is proposed. For completeness, we describe here the problem they consider and the proposed algorithm. considers the problem:

and proposes the following first-order recursion for its solution:

Key property in their analysis is the positive semi-definiteness of the feasible space. For a proper initialization and step size, show sublinear and linear convergence rates towards optimum, depending on the nature of $f$.

The Bi-Factored Gradient Descent (BFGD) Algorithm

In this section, we provide an overview of the Bi-Factored Gradient Descent (BFGD) algorithm for two problem settings in (1): $(i)$ $f$ being a $L$-smooth convex function and, $(ii)$ $f$ being $L$-smooth and $\mu$-strongly convex. For both cases, we assume a good initialization point $X_{0}=U_{0}V_{0}^{\top}$ is provided; for a discussion regarding initialization, see Section 5. Given $X_{0}$ and under proper assumptions, we further describe the theoretical guarantees that accompany BFGD.

As introduced in Section 1, BFGD is built upon non-convex gradient descent over each factor $U$ and $V$, written as

When $f$ is convex and smooth, BFGD follows exactly the motions in (13); in the case where $f$ is also strongly convex, BFGD is based on a different set of recursions, which we discuss in more detail in the rest of the section.

In , the authors describe a simple technique to transform problems similar to (1) into problems where we look for a square and PSD solution. The key idea is to lift the problem and introduce a stacked matrix of the two factors, as follows:

Following this idea, one can utilize algorithmic solutions designed only to work on square and PSD-based instances of (1), where $f$ is just $L$-smooth. Here, we use the Factored Gradient Descent (FGD) algorithm of on the $W$-space, as follows:

Then, it is easy to verify the following remark:

A natural question is whether this reduction gives a desirable convergence behavior. Since FGD solves for a different function $\hat{f}$ from the original $f$, the convergence analysis depends also on $\hat{f}$. When $f$ is convex and smooth, we can rely on the result from .

If $f$ is convex and $L$-smooth, then $\hat{f}$ is convex and $\tfrac{L}{2}$-smooth.

where the first inequality follows from the $L$-smoothness of $f$. ∎

Based on the above proposition, we use FGD to solve (2) with $\hat{f}$: its procedure is exactly (13), with a different step size:

While one can rely on the sublinear convergence analysis from , we provide a new guarantee with a weaker initial condition; see Section 4.

2 Using BFGD when f𝑓f is L𝐿L-Smooth and Strongly Convex

Assume $f$ function satisfies both properties in Definitions 2.1 and 2.2. In this case, we cannot simply rely on the lifting technique as above since $\hat{f}$ is clearly not strongly convex. Instead, we consider a slight variation, where we appropriately regularize the objective and force the solution pair $(\widehat{U},\widehat{V})$ to be “balanced”. This regularization is based on the set of optimal pairs $(U^{\star},V^{\star})$ in $\mathcal{X}^{\star}_{r}$, as defined in (10). In particular, given $\mathcal{X}^{\star}_{r}$, the equivalent optimization problem that “forces” convergence to balanced $(U^{\star},V^{\star})$ is as follows:

$g$ is convex and minimized at zero point; i.e., $\nabla g(0)=0$.

$g$ is $\mu_{g}$-strongly convex and $L_{g}$-smooth.

The necessity of the regularizer. As we show next, the theoretical guarantees of BFGD heavily depend on the condition number of the pair $(U^{\star},V^{\star})$ the algorithm converges to. In particular, one of the requirements of BFGD is that every estimate $U_{t}$ (resp. $V_{t}$) be “relatively close” to the convergent point $U^{\star}$ (resp. $V^{\star}$), such that their distance $\|U_{t}-U^{\star}\|_{F}$ is bounded by a function of $\sigma_{r}(U^{\star})$, for all $t$. Then, it is easy to observe that, for arbitrarily ill-conditioned $(U^{\star},V^{\star})\in\mathcal{X}_{r}^{*}$, such a condition might not be easily satisfied by BFGD per iterationEven if $UV^{\top}$ is close to $U^{\star}{V^{\star}}^{\top}$, the condition numbers of $U$ and $V$ can be much larger than the condition number of $UV^{\top}$., unless we “force” the sequence of estimates $(U_{t},V_{t}),~{}\forall t,$ to converge to a better conditioned pair $(U^{\star},V^{\star})$. This is the key role of regularizer $g$: it guarantees putative estimates $U_{t}$ and $V_{t}$ are not too ill-conditioned, per iteration.

An example of $g$ is the Frobenius norm (weighted by $\mu/2$), as proposed in . Other examples are sums of element-wise (at least) $\mu_{g}$-strongly convex and (at most) $L_{g}$-gradient Lipschitz functions (of the form $g(X)=\sum_{i,j}g_{ij}(X_{ij})$) with the optimum at zero. However, any other user-friendly $g$ can be selected, as long as it satisfies the above conditions. We show in this paper that any such regularizer results provably in convergence, with attractive convergence rate; see Section 6.2 for a toy example where the addition of $g$ leads to faster convergence rate in practice.

The BFGD algorithm. BFGD is a first-order, gradient descent algorithm for (16), that operates on the factored space $(U,V)$ in an alternating fashion. Principal components of BFGD is a proper step size selection and a “decent” initialization point. BFGD can be considered as the non-squared extension of FGD algorithm in , which is specifically designed to solve problems as in (2), for $U=V$ and $m=n$. The key differences with FGD though, other than the necessity of a regularizer $g$, are:

Our analysis leads to provable convergence results in the non-square case. Such a result cannot be trivially obtained from .

The main recursion followed is different in the two schemes: in the non-squared case, we update the left and right factors $(U,~{}V)$ with a different rule, according to which:

The parameter $\lambda>0$ is arbitrarily chosen.

Due to this new update rule, a slightly different and proper step size selection should be devised for BFGD. Our step size is selected as follows:

Compared to the step size proposed in (which is of the same form with (15)), our analysis drops the dependence to $\|\nabla f(\cdot)\|_{2}$ at the denominator. This is leads to a faster computed step size and highlights the non-necessity of this term for proof of convergence, i.e., the presence of the $\|\nabla f(\cdot)\|_{2}$ term is sufficient but not necessary.

The scheme is described in Algorithm 1. As shown in the next section, constant step size (17) is sufficient to lead to attractive convergence rates for BFGD, for $f$ $L$-smooth and $\mu$-strongly convex.

Local Convergence for BFGD

This section includes the main theoretical guarantees of BFGD, both for the cases of just smooth $f$, and $f$ being smooth and strongly convex. To provide such local convergence results, we assume that there is a known “good” initialization which ensures the following.

Assumption A1. Define $\kappa=\tfrac{\max\{L,L_{g}\}}{\min\{\mu,\mu_{g}\}}$ where $\mu_{g}$ and $L_{g}$ are the strong convexity and smoothness parameters of $g$, respectively. Then, we assume we are provided with a “good” initialization point $X_{0}=U_{0}V_{0}^{\top}$ such that:

For our analysis, we will use the following step size assumptions:

While these step sizes are different than the ones we use in practice, there is a constant-fraction connection between $\widehat{\eta}$ and $\eta$.

Let $(U_{0},V_{0})$ be such that Assumption A1 is satisfied. Then, (18) holds if (17) is satisfied, and (19) holds if (15) is satisfied.

The proof is provided in the Appendix A. By this lemma, our analysis below is equivalent—up to constants—to that if we were using the original step size $\eta$ of the algorithm. However, for clarity reasons and ease of exposition, we use $\widehat{\eta}$ below.

For the case of strongly convex $f$, both Assumption A1 and the step size depends on the strong convexity and smoothness parameters of $g$. When $\mu$ and $L$ are known a priori, this dependency can be removed since one can choose $g$ such that at least $\mu$-restricted strongly convex and at most $L$-smooth. Then, $\kappa$ becomes the condition number of $f$, and the step size depends only on $L$.

The following theorem proves that, under proper initialization, BFGD admits linear convergence rate, when $f$ is both $L$-smooth and $\mu$-restricted strongly convex.

for every $t\geq 0$, where the contraction parameter $\gamma_{t}$ satisfies:

The proof is provided in Section B. The theorem states that if $X^{\star}$ is (approximately) low-rank, the iterates converge to a close neighborhood of $X^{\star}_{r}$.

The above result can also be expressed w.r.t. the function value $f(UV^{\top})$, as follows:

Under the same initial condition with Theorem 4.2, Algorithm (1) satisfies the following recursion w.r.t. the distance of function values from the optimal point:

2 Local Sublinear Convergence

In Section 3.1, we showed that a lifting technique can reduce our problem (2) to a rank-constrained semidefinite program, and applying FGD from is exactly BFGD (13). While the sublinear convergence guarantee of FGD can also be applied to our problem, we provide an improved result.

Theorem 4.4 guarantees a local sublinear convergence with a looser initial condition. While requires $\min_{R\in O(r)}\left\|{W-W^{\star}R}\right\|_{F}\leq\frac{\sigma_{r}^{2}(W^{\star})}{100\sigma_{1}^{2}(W^{\star})}\cdot\sigma_{r}(W^{\star})$, our result requires that the initial distance to the $W^{\star}$ is merely a constant factor of $\sigma_{r}(W^{\star})$.

Initialization

In this section, we present initialization procedures for the case where $f$ is strongly convex and smooth. Our main theorem guarantees linear convergence in the factored space given that the initial point $(U_{0},V_{0})$ is within a ball around the closest target factors, with radius $O(\kappa^{-1/2}\sigma_{r}(X^{\star}_{r})^{1/2})$. To find such a solution, we propose an extension of the initialization in .

Consider an initial solution $U_{0}V_{0}^{\top}$ which is the best rank-$r$ approximation of

Combined with Lemma 5.14 in , which transforms a good initial solution from the original space to the factored space, the following corollary gives one sufficient condition for global convergence of BFGD with the SVD of (21) as initialization.

where $A_{0}\Sigma_{0}B_{0}$ is the SVD of $-\frac{1}{L}\nabla f(0)$ satisfies the initial condition of Theorem 4.2.

Corollary 5.2 requires weaker conditions than in order Theorem 4.2 be transformed to global guarantees. While our theoretical results can only guarantee global convergence for a well-conditioned problem ($\kappa$ close to one), we show in the experiments that the algorithm performs well in practice where the sufficient conditions are yet to be satisfied.

Experiments

In this section, we first provide comparison results regarding the actual computational complexity of SVD and matrix-matrix multiplication procedures; while such comparison is not thoroughly complete, it provides some evidence about the gains of optimizing over $U,V$ factors, in lieu of SVD-based rank-$r$ approximations. We also describe a toy example that highlights the effect of the regularizer $g$ in convergence rates, for strongly convex and smooth $f$. Next, we provide extensive results on the performance of BFGD, as compared with state of the art, for the following problem settings: $(i)$ affine rank minimization, where the objective is smooth and (restricted) strongly convex, $(ii)$ image denoising/recovery from a limited set of observed pixels, where the problem can be cast as a matrix completion problem, and $(iii)$ 1-bit matrix completion, where the objective is just smooth convex. In all cases, the task is to recover a low rank matrix from a set of observations, where our machinery naturally applies.

Figure 1 (left panel) shows execution time results for the algorithms under comparison, as a function of the dimension $m$. Rank $r$ is fixed to $r=100$. While both SVD and matrix multiplication procedures are known to have $O(m^{2}r)$ complexity, it is obvious that the latter on dense matrices is at least two-orders of magnitude faster than the former. In Table 2, we also report the approximation guarantees of some faster SVD subroutines, as compared to svds: while irblablk seems to be faster, it returns a very rough approximation of the singular values, when $r$ is relatively large. Similar findings are depicted in Figures 1 (middle and right panel).

2 The Role of Regularizer g𝑔g

In this section, we provide a simple example that illustrates the role of the regularizer $g$. As discussed in Section 3.2, $g$ forces our algorithm to converge to a well-conditioned factorization of $X^{\star}$. This regularizer not only enables us to control and guarantee convergence of BFGD, but also provides a better convergence rate, as we know next.

Figure 2 (left panel) shows the convergence behavior of BFGD, when $f$ and $f+g$ is used, with an ill-conditioned initial point $(U_{0},V_{0})$. It is obvious from the convergence plot that adding the regularizer results into faster convergence to an optimum. This difference in convergence rate is due to dependency on the condition numbers of $U^{\star}$ and $V^{\star}$ that the algorithm converges to. As shown in Figure 2 (right panel), the algorithm converges to a well-conditioned factorization of $X^{\star}$, while the condition number is not forced to decrease when there is no regularizer.

3 Affine Rank Minimization Using Noiselet Linear Maps

List of algorithms. We compare the following state-of-the-art algorithms: $(i)$ the Singular Value Projection (SVP) algorithm , a non-convex, projected gradient descent algorithm for (3), with constant step size selection (we study the case where $\mu=\nicefrac{{1}}{{3}}$, as it is the one that showed the best performance in our experiments), $(ii)$ the SparseApproxSDP extension to non-square cases for (5) in , based on , where a putative solution is refined via rank-1 updates from the gradientSparseApproxSDP in avoids computationally expensive operations per iteration, such as full SVDs. In theory, at the $r$-th iteration, these schemes guarantee to compute a $\tfrac{1}{r}$-approximate solutio, with rank at most $r$, i.e., achieves a sublinear rate., $(iii)$ the matrix completion algorithm in , which we call GuaranteedMCWe note that the original algorithm in is designed for the matrix completion problem, not the matrix sensing problem here., where the objective is (6), $(iv)$ the Procrustes Flow algorithm in for (6), and $(v)$ the BFGD algorithm.The algorithm in assumes step size that depends on RIP constants, which are NP-hard to compute; since no heuristic is proposed, we do not include this algorithm in the comparison list.

Implementation details. To properly compare the algorithms in the above list, we preset a set of parameters that are common. In all experiments, we fix the number of observations in $y$ to $p=C\cdot n\cdot r$, where $n\geq m$ in our cases, and for varying values of $C$. All algorithms in comparison are implemented in a Matlab environment, where no mex-ified parts present, apart from those used in SVD calculations; see below.

In all algorithms, we fix the maximum number of iterations to $T=4000$, unless otherwise stated. We use the same stopping criteria for the majority of algorithms as:

where $X_{t},~{}X_{t-1}$ denote the current and the previous estimates in the $X$ space and $\rm{tol}:=5\cdot 10^{-6}$. For SVD calculations, we use the lansvd implementation in PROPACK package . For fairness, we modified all the algorithms so that they exploit the true rank $r$; however, we observed that small deviations from the true rank result in relatively small degradation in terms of the reconstruction performance.In case the rank of $X^{\star}$ is unknown, one has to predict the dimension of the principal singular space. The authors in , based on ideas in , propose to compute singular values incrementally until a significant gap between singular values is found. For a more recent discussion on how to efficiently estimate the numerical rank of a matrix, refer to

In the implementation of BFGD, we set $g$ to be $\tfrac{1}{16}\cdot\|U^{\top}U-V^{\top}V\|_{F}^{2}$, as suggested in , for ease of comparison. Moreover, for our implementation of Procrustes Flow, we set the constant step size as $\mu:=\tfrac{2}{187}\cdot\{\tfrac{1}{\|U_{0}\|_{F}^{2}},\tfrac{1}{\|V_{0}\|_{F}^{2}\|}\}$, as suggested in . We use the implementation of , with random initialization (unless otherwise stated) and regularization type soft, as suggested by their implementation. In , we require an upper bound on the nuclear norm of $X^{\star}$; in our experiments we assume we know $\|X^{\star}\|_{*}$, which requires a full SVD calculation. Moreover, for our experiments, we set the curvature constant for the SparseApproxSDP implementation to its true value $C_{f}=1$.

For initialization, we consider the following settings: $(i)$ random initialization, where $X_{0}=U_{0}V_{0}^{\top}$ for some randomly selected $U_{0}$ and $V_{0}$ such that $\|X_{0}\|_{F}=1$, and $(ii)$ specific initialization, as suggested in each of the papers above. Our specific initialization is based on the discussion in Section 5, where $X_{0}=\mathcal{P}_{r}(-\frac{1}{L}\nabla f(0))$. Algorithms SVP, SparseApproxSDP and the solver in work with random initialization. For the initialization phase of , we consider two cases: $(i)$ the condition number $\kappa$ is known, where according to Theorem 3.3 in , we require $T_{\text{init}}:=\lceil 3\log(\sqrt{r}\cdot\kappa)+5\rceil$ SVP iterationsObserve that setting $T_{\text{init}}=1$ leads to spectral method initialization and the algorithm in for non-square cases, given sufficient number of samples., and $(ii)$ the condition number $\kappa$ is unknown, where we use Lemma 3.4 in .

Results using random initialization. Figure 3 depicts the convergence performance of the above algorithms w.r.t. total execution time. Top row corresponds to the case $m=n=1024$, bottom row to the case $m=2048,~{}n=4096$. For all cases, we fix $r=50$; from left to right, we decrease the number of available measurements, by decreasing the constant $C$. BFGD shows the best performance, compared to the rest of the algorithms. It is notable that BFGD performs better than SVP, by avoiding SVD calculations and employing a better step size selection.If our step size is used in SVP, we get slightly better performance, but not in a universal manner. For this setting, GuaranteedMC converges to a local minimum while SparseApproxSDP and Procrustes Flow show close to sublinear convergence rate.

To further show how the performance of each algorithms scales as dimension increases, we provide aggregated results in Tables 3-4. Observe that BFGD is one order of magnitude faster than the rest non-convex factorization algorithms. Table 5 shows the median time per iteration, spent by each algorithm, for both problem instances and $C=3$. Observe that SVP requires one order of magnitude more time to complete one iteration, mostly due to the SVD step. In stark contrast, all factorization-based approaches spend less time per iteration, as was expected by the discussion in Section 6.1.

Results using specific initialization. recently proved that random initialization is sufficient to lead to the optimum $X^{\star}$ for matrix sensing problems, while operating on the factors for the case where $X^{\star}\succeq 0$, i.e., $X^{\star}$ is square. We conjecture similar results can be proved for the non-square case, but we still consider specific initializations for completeness. In this case, we study the effect of initialization in the convergence performance of each algorithm. To do so, we focus only on the factorization-based algorithms: Procrustes Flow, GuaranteedMC, and BFGD. We consider two problem cases: $(i)$ all these schemes use our initialization procedure, and $(ii)$ each algorithm uses its own suggested initialization procedure. The results are depicted in Tables 6-7, respectively.

Using our initialization procedure for all algorithms, we observe that both Procrustes Flow and GuaranteedMC schemes can compute an approximation $\widehat{X}$ such that $\tfrac{\|\widehat{X}-X^{\star}\|_{F}}{\|X^{\star}\|_{F}}>10^{-1}$. In contrast, our approach achieves a solution $\widehat{X}$ that is close to the stopping criterion, i.e., $\tfrac{\|\widehat{X}-X^{\star}\|_{F}}{\|X^{\star}\|_{F}}\approx 10^{-6}$.

Using different initialization schemes per algorithm, the results are depicted in Table 7. We remind that GuaranteedMC is designed for matrix completion tasks, where the linear operator is a selection mask of the entries. Observe that Procrustes Flow’s performance improves significantly by using their proposed initialization: the idea is to perform SVP iterations to get to a good initial point; then switch to non-convex factored gradient descent for low per-iteration complexity. However, this initialization is computationally expensive: Procrustes Flow might end up performing several SVP iterations. This can be observed e.g., in the case $m=n=1024,~{}r=5$ and comparing the results in Tables 6-7: for this case, Procrustes Flow performs $T=4000$ iterations when our initialization is used and spends $\sim 200$ seconds, while in Table 7 it performs $T\ll 4000$ iterations, at least $20\%$ of them using SVP, and consumes $\sim 2000$ seconds.

4 Image Denoising as Matrix Completion Problem

In this example, we consider the matrix completion setting for an image denoising task: In particular, we observe a limited number of pixels from the original image and perform a low rank approximation based only on the set of measurements—similar experiments can be found in . We use real data images: While the true underlying image might not be low-rank, we apply our solvers to obtain low-rank approximations.

Figures 4-6 depict the reconstruction results for three image cases. In all cases, we compute the best $100$-rank approximation of each image (see e.g., the top middle image in Figure 4, where the full set of pixels is observed) and we observe only the $35\%$ of the total number of pixels, randomly selected—a realization is depicted in the top right plot in Figure 4. Given a fixed common tolerance level and the same stopping criterion as before, the top rows of Figures 4-6 show the recovery performance achieved by a range of algorithms under consideration—the peak signal to noise ration (PSNR), depicted in dB, corresponds to median values after 10 Monte-Carlo realizations. Our algorithm shows competitive performance compared to simple gradient descent schemes as SVP and Procrustes Flow, while being a fast and scalable solver. Table 8 contains timing results from 10 Monte Carlo random realizations for all image cases.

5 1-bit Matrix Completion

Similar to classic matrix completion results, we assume $\Omega$ is chosen uniformly at random, e.g., we assume $\Omega$ follows a binomial model, as in . Two natural choices for $\sigma$ function are: $(i)$ the logistic regression model, where $\sigma(x)=\tfrac{e^{x}}{1+e^{x}}$, and $(ii)$ the probit regression model, where $\sigma(x)=1-\Phi(-x/\sigma)$ for $\Phi$ being the cumulative Gaussian distribution function. Both models correspond to different noise assumptions: in the first case, noise is modeled according to the standard logistic distribution, while in the second case, noise follows standard Gaussian assumptions. Under this model, propose two convex relaxation algorithmic solutions to recover $X^{\star}$: $(i)$ the convex maximum log-likelihood estimator under nuclear norm and infinity norm constraints:

and $(ii)$ the the convex maximum log-likelihood estimator under only nuclear norm constraints. In both cases, $f(X)$ satisfies the expression in (7). proposes a spectral projected-gradient descent method for both these criteria; in the case where only nuclear norm constraints are present, SVD routines compute the convex projection onto norm balls, while in the case where both nuclear and infinity norm constraints are present, propose a alternating-direction method of multipliers (ADMM) solution, in order to compute the joint projection onto these sets.

Figure 7 depicts the recovery performance of BFGD, as compared to variants of (27) in . We consider their performance over different noise levels w.r.t. the normalized Frobenius norm distance $\tfrac{\|\widehat{X}-X^{\star}\|_{F}}{\|X^{\star}\|_{F}}$. As noted in , the performance of all algorithms is poor when $\sigma$ is too small or too large, while in between, for moderate noise levels, we observe better performance for all approaches.

By default, in all problem settings, we observe that the estimate of (27) is not of low rank: to compute the closest rank-$r$ approximation to that, we further perform a debias step via truncated SVD. The effect of the debias step is better illustrated in Figure 7, focusing on the differences between left and right plot: without such step, BFGD has a better performance in terms of $\tfrac{\|\widehat{X}-X^{\star}\|_{F}}{\|X^{\star}\|_{F}}$, within the “sweet” range of noise levels, compared to the convex analog in (27). Applying the debias step, both approaches have comparable performance, with that of (27) being slightly better.

Perhaps somewhat surprisingly, the performance of BFGD, in terms of estimating the correct sign pattern of the entries, is better than that of , even with the debias step. Figure 8 (left panel) illustrates the observed performances for various noise levels.

Finally, we study the performance of the algorithms under consideration as a function of the number of measurements, for fixed settings of dimensions $m=n=200$ and noise level $\sigma=0.244$. By the discussion above, such noise level leads to good performance from all schemes. We considered matrices $X^{\star}$ with rank $r\in\left\{3,~{}5,~{}10\right\}$ and generate $p=C\cdot n^{2}$, over a wide range of $0<C<1$. Figure 8 (right panel) shows the performance of BFGD and the approach for (27) in , in terms of the relative Frobenius norm of the error. All approaches do poorly when there are only $p<0.35\cdot n^{2}$ measurements, since this is near the noiseless information-theoretic limit. For higher numbers of measurements, the non-convex approach in BFGD returns more reasonable solutions and outperforms convex approaches, taking advantage of the prior knowledge on low-rankness of the solution.

Recommendation system using the MovieLens data set. We compare 1-bit matrix completion solvers on the 100k MovieLens data set. To do so, we repeat the experiment in Section 4.3 of : we use the MovieLens 100k, which consists of 100k movie ratings, from 1000 users on 1700 movies. Each user entry denotes the movie rating, ranging from 1 to 5. To convert this data set into 1-bit measurements, we convert these ratings to binary observations by comparing each rating to the average rating for the entire data set (which is approximately 3.5), according to . To evaluate the performance of the algorithms, we assume part of the observed ratings as unobserved (5k of them) and check if the estimate of $X^{\star}$, $\widehat{X}$, predicts the sign of these ratings. We perform ML estimation using logistic function $\sigma(x)=\tfrac{e^{x}}{1+e^{x}}$ in $f$.

We compare the following algorithms: $(i)$ the spectral projected gradient descent (SPG) implementation of (27) in for 1-bit matrix completion, $(ii)$ the standard matrix completion implementation TFOCS , where we observe the unquantized data set (actual values)Using TFOCS, we set the regularizer $\mu=10^{-3}$ as the parameter value that returned the best recovery results over a wide range of $\mu$ values., $(iii)$ BFGD for various values of rank parameter $r$. The results are shown Table 9 over 10 Monte Carlo realizations (i.e., we randomly selected 5k ratings as test sets and solved the problem for different runs of the algorithms). The values in Table 9 denote the accuracy in predicting whether the unobserved ratings are above or below the average rating of 3.5. BFGD shows competitive performance, compared to convex approaches. Moreover, setting the parameter $r$ is an “easier” and more intuitive task: our algorithm administers precise control on the rankness of the solution, which might lead to further interpretation of the results. Convex approaches lack of this property: the mapping between the regularization parameters and the number of rank-1 components in the extracted solution is highly non-linear. At the same time, BFGD shows much faster convergence to a good solution, which constitutes it a preferable algorithmic solution for large scale applications.

Appendix A Appendix

Proof of (18)$\Rightarrow$(17) Let $R^{\star}_{t}$ be the $r\times r$ orthogonal matrix such that ${\rm{\textsc{Dist}}}(U_{t},V_{t};X^{\star}_{r})=\left\|{W_{t}-W^{\star}R^{\star}_{t}}\right\|_{F}$. By the triangle inequality, we have

where $(i)$ is due to triangle inequality, $(ii)$ is due to Assumption A.1 $(iii)$ is due to the fact that $\sqrt{2}\cdot\sigma_{r}(X^{\star}_{r})^{\nicefrac{{1}}{{2}}}=\sigma_{r}(W^{\star})$ and $\kappa\geq 1$. The above bound holds for every $t=0,1,\ldots$.

where $(i)$ is due to the fact that $f$ is $L$-smooth and, $(ii)$ holds by adding and subtracting $U^{\star}V^{\star\top}$ and then applying triangle inequality. To bound the last two terms on the right hand side, we observe:

where $(i)$ is due to the triangle and Cauchy-Schwartz inequalities, $(ii)$ is by Assumption A1 and (A). Similarly, one can show that $\|U_{0}V_{0}^{\top}-U^{\star}{V^{\star}}^{\top}\|_{F}\leq\frac{7}{10}\cdot\left\|{W_{0}}\right\|_{2}^{2}$. Thus, (30) becomes:

Applying (A), (A), and the above bound, we obtain the desired result. $\blacksquare$

Appendix B Appendix

In the above formulations, we use as regularizer of $g$ function $\lambda=\tfrac{1}{2}$.

Our discussion below is based on the Assumption A.1, where:

holds for the current iterate. The last equality is due to the fact that $\sigma_{r}(W^{\star})=\sqrt{2}\cdot\sigma_{r}(X^{\star}_{r})^{1/2}$, for $(U^{\star},V^{\star})$ with “equal footing”. For the initial point $(U_{0},V_{0})$, (32) holds by the assumption of the theorem. Since the right hand side is fixed, (32) holds for every iterate, as long as ${\rm{\textsc{Dist}}}(U,V;X^{\star}_{r})$ decreases.

To show this, let $R\in O_{r}$ be the minimizing orthogonal matrix such that ${\rm{\textsc{Dist}}}(U,V;X^{\star}_{r})=\left\|{W-W^{\star}R}\right\|_{F}$; here, $O_{r}$ denotes the set of $r\times r$ orthogonal matrices such that $R^{\top}R=I$. Then, the decrease in distance can be lower bounded by

where the last equality is obtaining by substituting $W^{+}$, according to its definition above. To bound the first term on the right hand side, we use the following lemma; the proof is provided in Section B.

Suppose (32) holds for $W$. Let $\mu_{\min}=\min\left\{\mu,~{}\mu_{g}\right\}$ and $L_{\max}=\max\left\{L,~{}L_{g}\right\}$ for $(\mu,L)$ and $(\mu_{g},L_{g})$ the strong convexity and smoothness parameters pairs for $f$ and $g$, respectively. Then, the following inequality holds:

For the second term on the right hand side of (B), we obtain the following upper bound:

where $(a)$ follows from the fact $\left\|{A+B}\right\|_{F}^{2}\leq 2\left\|{A}\right\|_{F}^{2}+2\left\|{B}\right\|_{F}^{2}$, $(b)$ is due to the fact $\left\|{AB}\right\|_{F}\leq\left\|{A}\right\|_{F}\cdot\left\|{B}\right\|_{2}$, and $(c)$ follows from the observation that $\left\|{U}\right\|_{2},\left\|{V}\right\|_{2}\leq\left\|{W}\right\|_{2}$.

where we use the fact that $\sigma_{r}(W^{\star})=\sqrt{2}\cdot\sigma_{r}(X^{\star}_{r})^{1/2}$.

The above lead to the following recursion:

where $\gamma_{t}=1-\tfrac{\widehat{\eta}\cdot\mu_{\min}\cdot\sigma_{r}(X^{\star}_{r})}{5}$. By the definition of $\widehat{\eta}$ in (18), we further have:

where $(i)$ is by using (A) that connects $\|W\|_{2}$ with $\|W^{\star}\|_{2}$ as $\|W\|_{2}\geq\tfrac{9}{10}\|W^{\star}\|_{2}$, and $(ii)$ is due to the fact $\|W^{\star}\|_{2}=\sqrt{2}\cdot\sigma_{1}(X^{\star}_{r})^{1/2}$. $\blacksquare$

Before we step into the proof, we require a few more notations for simpler presentation of our ideas. We use another set of stacked matrices $Y=\begin{bmatrix}U\\ -V\end{bmatrix},~{}Y^{\star}=\begin{bmatrix}U^{\star}\\ -V^{\star}\end{bmatrix}.$ The error of the current estimate from the closest optimal point is denoted by the following $\Delta_{\times}$ matrix structures:

where, for the second term in (36), we use the fact that

the Cauchy-Schwarz inequality and the fact that $\|\Delta_{U}\|_{F},\|\Delta_{V}\|_{F}\leq\|\Delta_{W}\|_{F}$; the first term in (36) follows from:

where $(i)$ is due to the $\mu$-strong convexity of $f$, $(ii)$ is by adding and subtracting $f(X^{\star})$; observe that $f(X^{\star})=f(U^{\star}{V^{\star}}^{\top})$ if and only if $\text{rank}(X^{\star})=r$, and $(iii)$ is due to the $L$-smoothness of $f$ and the fact that $\nabla f(X^{\star})=0$ (for the middle term), and due to the inequality [81, eq. (2.1.7)] (for the first term):

where $(a)$ follows from the “balance” assumption in $\mathcal{X}_{r}^{\star}$:

for the first term, and the fact that $\nabla g$ is symmetric, and therefore

for the second term; $(b)$ follows from the fact

and the Cauchy-Schwarz inequality on the second term in (38), and

where $(i)$ follow from the strong convexity, $(ii)$ is due to (37), and $(iii)$ is by construction of $g$ where $\nabla g(0)=0$. Furthermore, (B1) can be bounded below as follows:

and the first inequality holds by the fact that the inner product of two PSD matrices is non-negative.

At this point, we have all the required components to compute the desired lower bound. Combining (A1) and (B1), we get

where, in order to obtain the last inequality, we borrow the following Lemma by :

For convenience, we further lower bound the right hand side of this lemma by: $2\cdot\left(\sqrt{2}-1\right)\cdot\sigma_{r}(W^{\star})^{2}\cdot\left\|{\Delta_{W}}\right\|_{F}^{2}\geq\frac{4\sigma_{r}(W^{\star})^{2}}{5}\left\|{\Delta_{W}}\right\|_{F}^{2}.$

Given the definitions of $\mu_{\min}$ and $L_{\max}$, we have:

where in $(i)$ we used the definitions of $\mu_{\min}$ and $L_{\max}$. Note that we have not used the condition (32). It follows from (32) that

where we use the AM-GM inequality. Plugging (40) and (41) in (B), it is easy to obtain:

Appendix C Appendix

The proof follows the same framework of the sublinear convergence proof in . We use the following general lemma to prove the sublinear converegence.

Suppose that a sequence of iterates $\{W_{t}\}_{t=0}^{T}$ satisfies the following conditions

for all $t=0,\ldots,T-1$ and some values $\alpha,\beta>0$ independent of the iterates. Then it is guaranteed that

Define $\delta_{t}=f(W_{t}W_{t}^{\top})-f(W^{\star}{W^{\star}}^{\top})$. If we get $\delta_{T_{0}}\leq 0$ at some $T_{0}<T$, the desired inequality holds because the first hypothesis guarantees $\{\delta_{t}\}_{t=0}^{T}$ to be non-increasing. Hence, we can only consider the time $t$ where $\delta_{t},\delta_{t+1}\geq 0$. We have

where (a) follows from the first hypothesis, (b) follows from the second hypothesis, (c) follows from that $\delta_{t+1}\leq\delta_{t}$ by the first hypothesis. Dividing by $\delta_{t}\cdot\delta_{t+1}$, we obtain

Then we obtain the desired result by telescoping the above inequality. ∎

Now it suffices to show BFGD provides a sequence $\{W_{t}\}_{t=0}^{T}$ satisfies the hypotheses of Lemma C.1.

Obtaining (42) Although $f$ is non-convex over the factor space, it is reasonable to obtain a new estimate (with a carefully chosen steplength) which is no worse than the current one, because the algorithm takes a gradient step.

Let $f$ be a $L$-smooth convex function. Moreover, consider the recursion in Let $X=WW^{\top}$ and ${X}^{\scalebox{0.7}{+}}={W}^{\scalebox{0.7}{+}}{{W}^{\scalebox{0.7}{+}}}^{\top}$ be two consecutive estimates of BFGD. Then

Since we can fix the steplength $\eta$ based on the initial solution so that it is independent of the following iterates, we have obtained the first hypothesis of Lemma C.1.

Obtaining (43) Consider the following assumption.

Trivially (A) holds for $U_{0}$ and $V_{0}$. Now we provide key lemmas, and then the convergence proof will be presented.

Assume that (A) holds for $W$. Then we have

Combining the above two lemmas, we obtain

Plugging (44) and (45) in Lemma C.1, we obtain the desired result. $\blacksquare$

Using the Cauchy-Schwarz inequality, the second term can be bounded as follows.

To bound the third term of (C.1), we have

Plugging (C.1), (C.1), and (C.1) to (C.1), we obtain

where the last inequality follows from the condition of the steplength $\eta$. This completes the proof. $\blacksquare$

C.2 Proof of Lemma C.3

(a) follows from the convexity of $f$, (b) follows from the Cauchy-Schwarz inequality, and (c) follows from Lemma C.5. $\blacksquare$

C.3 Proof of Lemma C.5

Define $Q_{W}$, $Q_{W^{\star}}$, and $Q_{\Delta_{W}}$ as the projection matrices of the column spaces of $W$, $W^{\star}$, and $\Delta_{W}=W-W^{\star}R$, respectively. We have

where (a) follows from the Cauchy-Schwarz inequality and the fact $\left\|{AB}\right\|_{F}\leq\left\|{A}\right\|_{2}\cdot\left\|{B}\right\|_{F}$, and (b) follows from that $W-W^{\star}$ lies on the column space spanned by $W$ and $W^{\star}$. To bound the terms in (LABEL:eqn:step2_errorbound), we obtain

where $W^{\dagger}$ and ${W^{\star}}^{\dagger}$ are the Moore-Penrose pseudoinverses of $W$ and $W^{\star}$. Plugging the above into (LABEL:eqn:step2_errorbound), we get

where (a) follows from (A). $\blacksquare$

C.4 Proof of Lemma C.4

For this proof, we borrow a lemma from . Although the assumption for the lemma is stronger than Assumption (A), but a slight modification of the proof leads to the following lemma from Assumption (A).

Let Assumption (A) hold and $f(W^{+}{W^{+}}^{\top})\geq f(W^{\star}{W^{\star}}^{\top})$. Then the following lower bound holds:

where (a) follows from Lemma C.6, (b) follows from the convexity of $f$, the hypothesis of the lemma, and Lemma C.2 as follows.

Appendix D Appendix

Let us first obtain an upper bound on the first term. We have

where (a) follows from the triangle inequality, and (b) is due to Mirsky’s theorem . Plugging this bound into (52), we get

Now we bound the first term of (53). We have

where (a) follows from the $L$-smoothness, (b) and (c) follow from the $\mu$-strong convexity. Then it follows that

Applying this inequality to (53), we get the desired inequality. $\blacksquare$

References