Fixed Point and Bregman Iterative Methods for Matrix Rank Minimization

Introduction

The matrix rank minimization problem can be written as

In this paper, we are interested in methods for solving the affinely constrained matrix rank minimization problem

is a special case of (1), where $X$ and $M$ are both $m\times n$ matrices and $\Omega$ is a subset of index pairs $(i,j).$ The so called collaborative filtering problem (Rennie-Srebro-2005; Srebro-thesis-2004) can be cast as a matrix completion problem. Suppose users in an online survey provide ratings of some movies. This yields a matrix $M$ with users as rows and movies as columns whose $(i,j)$-th entry $M_{ij}$ is the rating given by the $i$-th user to the $j$-th movie. Since most users rate only a small portion of the movies, we typically only know a small subset $\{M_{ij}|(i,j)\in\Omega\}$ of the entries. Based on the known ratings of a user, we want to predict the user’s ratings of the movies that the user did not rate; i.e., we want to fill in the missing entries of the matrix. It is commonly believed that only a few factors contribute to an individual’s tastes or preferences for movies. Thus the rating matrix $M$ is likely to be of numerical low rank in the sense that relatively few of the top singular values account for most of the sum of all of the singular values. Finding such a low-rank matrix $M$ corresponds to solving the matrix completion problem (1).

When the matrix $X$ is diagonal, problem (1) reduces to the cardinality minimization problem

The basis pursuit problem has received an increasing amount of attention since the emergence of the field of compressed sensing (CS) (Candes-Romberg-Tao-2006; Donoho-2006). Compressed sensing theories connect the NP-hard problem (1.1) to the convex and computationally tractable problem (1.1) and provide guarantees for when an optimal solution to (1.1) gives an optimal solution to (1.1). In the cardinality minimization and basis pursuit problems (1.1) and (1.1), $b$ is a vector of measurements of the signal $x$ obtained using the sampling matrix $A$. The main result of compressed sensing is that when the signal $x$ is sparse, i.e., $k:=\|x\|_{0}\ll n,$ we can recover the signal by solving (1.1) with a very limited number of measurements, i.e., $m\ll n$, when $A$ is a Gaussian random matrix or when it corresponds to a partial Fourier transformation. Note that if $b$ is contaminated by noise, the constraint $Ax=b$ in (1.1) must be relaxed, resulting in either the problem

where $\theta$ and $\mu$ are parameters and $\|x\|_{2}$ denotes the Euclidean norm of a vector $x$.. Algorithms for solving (1.1) and its variants (1.1) and (1.16) have been widely investigated and many algorithms have been suggested including convex optimization methods ((Candes-Romberg-2005-l1-magic; Figueiredo-Nowak-Wright-2007; Hale-Yin-Zhang-2007; Kim-Koh-Lustig-Boyd-Gorinevsky-2007; vandenBerg-Friedlander-2008)) and heuristic methods ((Tibshirani-1996; Donoho-Tsaig-Drori-Starck-2006; Tropp-2006; Donoho-Tsaig-2006; Dai-Milenkovie-2008)).

2 Nuclear norm minimization

Nuclear norm and Operator norm. Assume that the matrix $X$ has $r$ positive singular values of $\sigma_{1}\geq\sigma_{2}\geq\ldots\geq\sigma_{r}>0$. The nuclear norm of $X$ is defined as the sum of its singular values, i.e.,

The operator norm of matrix $X$ is defined as the largest singular value of $X$, i.e.,

The nuclear norm is also known as Schatten 1-norm or Ky Fan norm. Using it as an approximation to $\mbox{rank}(X)$ in (1) yields the nuclear norm minimization problem

As in the basis pursuit problem, if $b$ is contaminated by noise, the constraint $\mathcal{A}(X)=b$ must be relaxed, resulting in either the problem

For the matrix completion problem (1), the corresponding nuclear norm minimization problem is

Candès et al.(Candes-Recht-2008) proved the following result.

Let $M$ be an $n_{1}\times n_{2}$ matrix of rank $r$ with SVD

where the family $\{u_{k}\}_{1\leq k\leq r}$ is selected uniformly at random among all families of $r$ orthonormal vectors, and similarly for the family $\{v_{k}\}_{1\leq k\leq r}$. Let $n=\max(n_{1},n_{2})$. Suppose we observe $m$ entries of $M$ with locations sampled uniformly at random. Then there are constants $C$ and $c$ such that if

the minimizer to the problem (1.2) is unique and equal to $M$ with probability at least $1-cn^{-3}$. In addition, if $r\leq n^{1/5}$, then the recovery is exact with probability at least $1-cn^{-3}$ provided that

This theorem states that a surprisingly small number of entries are sufficient to complete a low-rank matrix with high probability.

Recently, this result was strengthened by Candès and Tao in (Candes-Tao-2009), where it is proved that under certain incoherence conditions, the number of samples $m$ that are required is only $O(nr\log n).$

The dual problem corresponding to the nuclear norm minimization problem (1.2) is

where $\mathcal{A}^{*}$ is the adjoint operator of $\mathcal{A}$. Both (1.2) and (1.2) can be rewritten as equivalent semidefinite programming (SDP) problems. The SDP formulation of (1.2) is:

where $\mbox{Tr}(X)$ denotes the trace of the square matrix $X$. The SDP formulation of (1.2) is:

Thus to solve (1.2) and (1.2), we can use SDP solvers such as SeDuMi (Sturm-1999) and SDPT3 (Tutuncu-Toh-Todd-2003) to solve (1.2) and (1.2). Note that the number of variables in (1.2) is $\frac{1}{2}(m+n)(m+n+1)$. SDP solvers cannot usually solve a problem when $m$ and $n$ are both much larger than $100.$

Recently, Liu and Vandenberghe (Liu-Vandenberghe-2008) proposed an interior-point method for another nuclear norm approximation problem

Liu and Vandenberghe (Liu-Vandenberghe-2008) proposed a customized method for computing the scaling direction in an interior point method for solving the SDP (1.2). The complexity of each iteration in their method was reduced from $O(p^{6})$ to $O(p^{4})$ when $m=O(p)$ and $n=O(p)$; thus they were able to solve problems up to dimension $m=n=350.$

where $\|X\|_{F}$ denotes the Frobenius norm of the matrix $X$:

It is known that as long as $r$ is chosen to be sufficiently larger than the rank of the optimal solution matrix of the nuclear norm problem (1.2), this low-rank factorization problem is equivalent to the nuclear norm problem (1.2) (see e.g., (Recht-Fazel-Parrilo-2007)). The advantage of this low-rank factorization formulation is that both the objective function and the constraints are differentiable. Thus gradient-based optimization algorithms such as conjugate gradient algorithms and augmented Lagrangian algorithms can be used to solve this problem. However, the constraints in this problem are nonconvex, so one can only be assured of obtaining a local minimizer. Also, how to choose $r$ is still an open question.

Our algorithms have some similarity with the SVT algorithm in that they make use of matrix shrinkage (see Section 2). However, other than that, they are greatly different. All of our methods are based on a fixed point continuation (FPC) algorithm which uses an operator splitting technique for solving (1.23). By adopting a Monte Carlo approximate SVD in the FPC, we get an algorithm, which we call FPCA (Fixed Point Continuation with Approximate SVD), that usually gets the optimal solution to (1) even if the condition of Theorem 1.1, or those for the affine constrained case, are violated. Moreover, our algorithm is much faster than state-of-the-art SDP solvers such as SDPT3 applied to (1.2). Also, FPCA can recover matrices of moderate rank that cannot be recovered by SDPT3, SVT, etc. with the same amount of samples. For example, for matrices of size $1000\times 1000$ and rank 50, FPCA can recover them with a relative error of $10^{-5}$ in about 3 minutes by sampling only 20 percent of the matrix elements. As far as we know, there is no other method that has as good a recoverability property.

3 Outline and Notation

Fixed point iterative algorithm

Our fixed point iterative algorithm for solving (1.23) is the following simple two-line algorithm:

where $S_{\nu}(\cdot)$ is the matrix shrinkage operator which will be defined later.

where $g^{*}=A^{\top}(Ax^{*}-b)$. For any $\tau>0$, (2.4) is equivalent to

Note that the operator $T(\cdot):=\tau\mu\mbox{SGN}(\cdot)+\tau g(\cdot)$ on the right hand side of (2.5) can be split into two parts: $T(\cdot)=T_{1}(\cdot)-T_{2}(\cdot),$ where $T_{1}(\cdot)=\tau\mu\mbox{SGN}(\cdot)+I(\cdot)$ and $T_{2}(\cdot)=I(\cdot)-\tau g(\cdot)$.

Letting $y=T_{2}(x^{*})=x^{*}-\tau A^{\top}(Ax^{*}-b)$, (2.5) is equivalent to

Note that (2.6) is actually the optimality conditions for the following convex problem

This problem has a closed form optimal solution given by the so called shrinkage operator:

Thus, the fixed point iterative algorithm is given by

Motivated by this work, we develop a fixed point iterative algorithm for (1.23). Since the objective function in (1.23) is convex, $X^{*}$ is the optimal solution to (1.23) if and only if

Hence, we get the following optimality conditions for (1.23):

Now based on the optimality conditions (2.10), we can develop a fixed point iterative scheme for solving (1.23) by adopting the operator splitting technique described at the beginning of this section. Note that (2.10) is equivalent to

In the following we will prove that the matrix shrinkage operator applied to $Y^{*}$ gives the optimal solution to (2.14). First, we need the following definitions.

To verify that (2.15) satisfies (2.18), consider the following two cases:

Case 1: $\gamma_{1}\geq\gamma_{2}\geq\ldots\geq\gamma_{t}>\nu.$ In this case, choosing $X$ as above, with $r=t,U=U_{Y},V=V_{Y}$ and $\sigma=s_{\nu}(\gamma)=\gamma-\nu e$, where $e$ is a vector of $r$ ones, and choosing $\bar{\sigma}=0$ (i.e., $W=\mathbf{0}$) satisfies (2.18).

Case 2: $\gamma_{1}\geq\gamma_{2}\geq\ldots\geq\gamma_{k}>\nu\geq\gamma_{k+1}\geq\ldots\geq\gamma_{t}.$ In this case, by choosing $r=k,\hat{U}(1:t)=U_{Y},\hat{V}(1:t)=V_{Y},\sigma=s_{\nu}((\gamma_{1},\ldots,\gamma_{k}))$ and $\bar{\sigma}_{1}=\gamma_{k+1}/\nu,\ldots,\bar{\sigma}_{t-k}=\gamma_{t}/\nu,\bar{\sigma}_{t-k+1}=\ldots=\bar{\sigma}_{m-r}=0,$ $X$ and $W$ satisfy (2.18).

Note that in both cases, $X$ can be written as the form in (2.15) based on the way we construct $X$. ∎

Based on the above we obtain the fixed point iterative scheme (2) stated at the beginning of this section for solving problem (1.23).

Moreover, from the discussion following Theorem 2.1 we have

$X^{*}$ is an optimal solution to problem (1.23) if and only if $X^{*}=S_{\tau\mu}(h(X^{*}))$, where $h(\cdot)=I(\cdot)-\tau g(\cdot).$

Convergence results

In this section, we analyze the convergence properties of the fixed point iterative scheme (2). Before we prove the main convergence result, we need some lemmas.

Without loss of generality, we assume $m\leq n$. Assume SVDs of $Y_{1}$ and $Y_{2}$ are $Y_{1}=U_{1}\Sigma V_{1}^{\top}$ and $Y_{2}=U_{2}\Gamma V_{2}^{\top}$, respectively, where

$\bar{\sigma}=(\sigma_{1}-\nu,\ldots,\sigma_{k}-\nu)$ and $\bar{\gamma}=(\gamma_{1}-\nu,\ldots,\gamma_{l}-\nu)$. Thus,

where $U=U_{1}^{\top}U_{2},V=V_{1}^{\top}V_{2}$ are clearly orthogonal matrices. Now let us derive an upper bound for $\mbox{Tr}(Y_{1}^{\top}Y_{2}-\bar{Y}_{1}^{\top}\bar{Y}_{2})$. It is known that an orthogonal matrix $U$ is a maximizing matrix for the problem

if and only if $AU$ is positive semidefinite matrix (see 7.4.9 in (Horn-Johnson-book-1985)). It is also known that when $AB$ is positive semidefinite,

Thus, $\mbox{Tr}((\Sigma-\bar{\Sigma})^{\top}U(\Gamma-\bar{\Gamma})V^{\top})$, $\mbox{Tr}((\Sigma-\bar{\Sigma})^{\top}U\bar{\Gamma}V^{\top})$ and $\mbox{Tr}(\bar{\Sigma}U(\Gamma-\bar{\Gamma})V^{\top})$ achieve their maximum, if and only if $(\Sigma-\bar{\Sigma})^{\top}U(\Gamma-\bar{\Gamma})V^{\top}$, $(\Sigma-\bar{\Sigma})^{\top}U\bar{\Gamma}V^{\top}$ and $\bar{\Sigma}U(\Gamma-\bar{\Gamma})V^{\top}$ are all positive semidefinite. Applying (3.4) to these three terms, we get $\mbox{Tr}((\Sigma-\bar{\Sigma})^{\top}U(\Gamma-\bar{\Gamma})V^{\top})\leq\sum_{i}\sigma_{i}(\Sigma-\bar{\Sigma})\sigma_{i}(\Gamma-\bar{\Gamma})$, $\mbox{Tr}((\Sigma-\bar{\Sigma})^{\top}U\bar{\Gamma}V^{\top})\leq\sum_{i}\sigma_{i}(\Sigma-\bar{\Sigma})\sigma_{i}(\bar{\Gamma})$ and $\mbox{Tr}(\bar{\Sigma}U(\Gamma-\bar{\Gamma})V^{\top})\leq\sum_{i}\sigma_{i}(\bar{\Sigma})\sigma_{i}(\Gamma-\bar{\Gamma}).$ Thus, without loss of generality, assuming $k\leq l\leq s\leq t$, we have,

since $t\geq s$ and $\sigma_{i}^{2}+\gamma_{i}^{2}-2\sigma_{i}\gamma_{i}\geq 0$. Also, since the function $f(x):=2\gamma_{i}x-x^{2}$ is monotonely increasing in $(-\infty,\gamma_{i}]$ and $\sigma_{i}<\nu\leq\gamma_{i},i=k+1,\ldots,l$,

Also, $D(Y_{1},Y_{2})$ achieves its minimum value if and only if $\mbox{Tr}((\Sigma-\bar{\Sigma})^{\top}U(\Gamma-\bar{\Gamma})V^{\top})$, $\mbox{Tr}((\Sigma-\bar{\Sigma})^{\top}U\bar{\Gamma}V^{\top})$ and $\mbox{Tr}(\bar{\Sigma}U(\Gamma-\bar{\Gamma})V^{\top})$ achieve their maximum values simultaneously.

Furthermore, if equality in (3.1) holds, i.e., $D(Y_{1},Y_{2})$ achieves its minimum, and its minimum is zero, then $k=l$, $s=t$, and $\sigma_{i}=\gamma_{i},i=k+1,\ldots,s$, which further implies $\Sigma-\bar{\Sigma}=\Gamma-\bar{\Gamma}$ and $\mbox{Tr}((\Sigma-\bar{\Sigma})^{\top}U(\Gamma-\bar{\Gamma})V^{\top})$ achieves its maximum. By applying the result 7.4.13 in (Horn-Johnson-book-1985), we get

To conclude, clearly $\|S_{\nu}(Y_{1})-S_{\nu}(Y_{2})\|_{F}=\|Y_{1}-Y_{2}\|_{F}$ if (3.9) holds. ∎

The following two lemmas and theorem and their proofs are analogous to results and their proofs in Hale et al.(Hale-Yin-Zhang-2007).

Let $\mathcal{A}X=A\mbox{vec}(X)$ and assume that $\tau\in(0,2/\lambda_{max}(A^{\top}A))$. Then the operator $h(\cdot)=I(\cdot)-\tau g(\cdot)$ is non-expansive, i.e., $\|h(X)-h(X^{\prime})\|_{F}\leq\|X-X^{\prime}\|_{F}$. Moreover, $h(X)-h(X^{\prime})=X-X^{\prime}$ if and only if $\|h(X)-h(X^{\prime})\|_{F}=\|X-X^{\prime}\|_{F}$.

First, we note that since $\tau\in(0,2/\lambda_{max}(A^{\top}A))$, $-1<\lambda_{i}(I-\tau A^{\top}A)\leq 1,\forall i,$ where $\lambda_{i}(I-\tau A^{\top}A)$ is the $i$-th eigenvalue of $I-\tau A^{\top}A$. Hence,

Moreover, $\|h(X)-h(X^{\prime})\|_{F}=\|X-X^{\prime}\|_{F}$ if and only if the inequalities above are equalities, which happens if and only if

i.e., if and only if $h(X)-h(X^{\prime})=X-X^{\prime}.$∎

Let $X^{*}$ be an optimal solution to problem (1.23), $\tau\in(0,2/\lambda_{max}(A^{\top}A))$ and $\nu=\tau\mu$. Then $X$ is also an optimal solution to problem (1.23) if and only if

The “only if” part is an immediate consequence of Corollary 1. For the “if” part, from Lemmas 1 and 2,

Hence, both inequalities hold with equality. Therefore, first using Lemma 1 and then Lemma 2 we obtain

which implies $S_{\nu}(h(X))=X$ since $S_{\nu}(h(X^{*}))=X^{*}$. It then follows from Corollary 1 that $X$ is an optimal solution to problem (1.23). ∎

We now claim that the fixed-point iterations (2) converge to an optimal solution of problem (1.23).

The sequence $\{X^{k}\}$ generated by the fixed point iterations with $\tau\in(0,2/\lambda_{max}(A^{\top}A))$ converges to some $X^{*}\in\mathcal{X}^{*},$ where $\mathcal{X}^{*}$ is the set of optimal solutions of problem (1.23).

Since both $S_{\nu}(\cdot)$ and $h(\cdot)$ are non-expansive, $S_{\nu}(h(\cdot))$ is also non-expansive. Therefore, $\{X^{k}\}$ lies in a compact set and must have a limit point, say $\bar{X}=\lim_{j\rightarrow\infty}X^{k_{j}}.$ Also, for any $X^{*}\in\mathcal{X}^{*}$,

which means that the sequence $\{\|X^{k}-X^{*}\|_{F}\}$ is monotonically non-increasing. Therefore,

where $\bar{X}$ can be any limit point of $\{X^{k}\}$. By the continuity of $S_{\nu}(h(\cdot))$, the image of $\bar{X}$,

is also a limit point of $\{X^{k}\}$. Therefore, we have

which allows us to apply Lemma 3 to get that $\bar{X}$ is an optimal solution to problem (1.23).

Finally, by setting $X^{*}=\bar{X}\in\mathcal{X}^{*}$ in (3.11), we get that

i.e., $\{X^{k}\}$ converges to its unique limit point $\bar{X}.$ ∎

Fixed point continuation

In this section, we discuss a continuation technique (i.e., homotopy approach) for accelerating the convergence of the fixed point iterative algorithm (2).

Inspired by the work of Hale et al.(Hale-Yin-Zhang-2007), we first describe a continuation technique to accelerate the convergence of the fixed point iteration (2). Our fixed point continuation (FPC) iterative scheme for solving (1.23) is outlined below.

The parameter $\eta_{\mu}$ determines the rate of reduction of the consecutive $\mu_{k}$, i.e.,

2 Stopping criteria for inner iterations

Note that in the fixed point continuation algorithm, in the $k$-th inner iteration we solve problem (1.23) for a fixed $\mu=\mu_{k}$. There are several ways to determine when to stop this inner iteration, decrease $\mu$ and go to the next inner iteration. The optimality conditions for (1.23) is given by (2.11a) and (2.11b). Thus we can use the following condition as a stopping criterion:

where $gtol$ is a small positive parameter. However, the expense of computing the largest singular value of a large matrix greatly decreases the speed of the algorithm. Hence, we do not use this criterion as a stopping rule for large matrices. Instead, we use the criterion

where $xtol$ is a small positive number, since when $X^{k}$ gets close to an optimal solution $X^{*}$, the distance between $X^{k}$ and $X^{k+1}$ should become very small.

3 Debiasing

Debiasing is another technique that can improve the performance of FPC. Debiasing has been used in compressed sensing algorithms for solving (1.1) and its variants, where debiasing is performed after a support set $\mathcal{I}$ has been tentatively identified. Debiasing is the process of solving a least squares problem restricted to the support set $\mathcal{I}$, i.e., we solve

where $A_{\mathcal{I}}$ is a submatrix of $A$ whose columns correspond to the support index set $\mathcal{I}$, and $x_{\mathcal{I}}$ is a subvector of $x$ corresponding to $\mathcal{I}$.

where $r$ is the rank of current matrix $X^{k}$. Because debiasing can be costly, we use a rule proposed in (Wen-Yin-Goldfarb-Zhang-2009) to decide when to do it. In the continuation framework, we know that in each subproblem with a fixed $\mu$, $\|X_{k+1}-X_{k}\|_{F}$ converges to zero, and $\|g\|_{2}$ converges to $\mu$ when $X_{k}$ converges to the optimal solution of the subproblem. We therefore choose to do debiasing when $\|g\|_{2}/\|X_{k+1}-X_{k}\|_{F}$ becomes large because this indicates that the change between two consecutive iterates is relatively small. Specifically, we call for debiasing in the solver FPC3 (see Section 7) when $\|g\|_{2}/\|X_{k+1}-X_{k}\|_{F}>10.$

Bregman iterative algorithm

Algorithm FPC is designed to solve (1.23), an optimal solution of which approaches an optimal solution of the nuclear norm minimization problem (1.2) as $\mu$ goes to zero. However, by incorporating FPC into a Bregman iterative technique, we can solve (1.2) by solving a limited number of instances of (1.23), each corresponding to a different $b$.

Given a convex function $J(\cdot)$, the Bregman distance (Bregman-1967) of the point $u$ from the point $v$ is defined as

where $p\in\partial J(v)$ is some subgradient in the subdifferential of $J$ at the point $v$.

Bregman iterative regularization was introduced by Osher et al.in the context of image processing (Osher-Burger-Goldfarb-Xu-Yin-2005). Specifically, in (Osher-Burger-Goldfarb-Xu-Yin-2005), the Rudin-Osher-Fatemi (Rudin-Osher-Fatemi-1992) model

was extended to an iterative regularization model by replacing the total variation functional

by the Bregman distance with respect to $J(u)$. This Bregman iterative regularization procedure recursively solves

for $k=0,1,\ldots$ starting with $u^{0}=\mathbf{0}$ and $p^{0}=\mathbf{0}$. Since (5.3) is a convex programming problem, the optimality conditions are given by $\mathbf{0}\in\partial J(u^{k+1})-p^{k}+u^{k+1}-b,$ from which we get the update formula for $p^{k+1}:$

Therefore, the Bregman iterative scheme is given by

Interestingly, this turns out to be equivalent to the iterative process

which can be easily implemented using existing algorithms for (5.2) with different inputs $b$.

Subsequently, Yin et al.(Yin-Osher-Goldfarb-Darbon-2008) proposed solving the basis pursuit problem (1.1) by applying the Bregman iterative regularization algorithm to

for $J(x)=\mu\|x\|_{1},$ and obtained the following two equivalent iterative schemes analogous to (5) and (5), respectively:

$x^{0}\leftarrow\mathbf{0},p^{0}\leftarrow\mathbf{0},$

$x^{k+1}\leftarrow\mbox{argmin}_{x}D_{J}^{p^{k}}(x,x^{k})+\displaystyle{\frac{1}{2}}\|Ax-b\|_{2}^{2}$

$p^{k+1}\leftarrow p^{k}-A^{\top}(Ax^{k+1}-b)$

$b^{0}\leftarrow\mathbf{0},x^{0}\leftarrow\mathbf{0},$

$x^{k+1}\leftarrow\mbox{argmin}_{x}J(x)+\displaystyle{\frac{1}{2}}\|Ax-b^{k+1}\|_{2}^{2}.$

One can also use the Bregman iterative regularization algorithm applied to the unconstrained problem (1.23) to solve the nuclear norm minimization problem (1.2). That is, one iteratively solves (1.23) by

Equivalently, one can also use the following iterative scheme:

Thus, our Bregman iterative algorithm for nuclear norm minimization (1.2) can be outlined as follows.

The last step can be solved by Algorithm FPC.

An approximate SVD based FPC algorithm: FPCA

The outputs $\sigma_{t}(C),t=1,\ldots,k_{s}$ are approximations to the largest $k_{s}$ singular values and $H_{k_{s}}^{(t)},t=1,\ldots,k$ are approximations to the corresponding left singular vectors of the matrix $A$. Thus, the SVD of $A$ is approximated by

Drineas et al.(Drineas-Kannan-Mahoney-2006) prove that with high probability, the following estimate holds for both $\xi=2$ and $\xi=F$:

There are many ways to choose the probabilities $p_{i}$. In our numerical experiments in Section 7, we used the simplest one, i.e., we set all $p_{i}$ equal to $1/n$. For other choices of $p_{i}$, see (Drineas-Kannan-Mahoney-2006) and the references therein.

In our numerical experiments, we set $k_{s}$ using the following procedure. In the $k$-th iteration, when computing the approximate SVD of $Y^{k}=X^{k}-\tau g^{k}$, we set $k_{s}$ equal to the number of components in $\bar{s}_{k-1}$ that are no less than $\epsilon_{k_{s}}\max\{\bar{s}_{k-1}\},$ where $\epsilon_{k_{s}}$ is a small positive number and $\max\{\bar{s}_{k-1}\}$ is the largest component in the vector $\bar{s}_{k-1}$ used to form $X^{k}=U^{k-1}\mbox{Diag}(\bar{s}_{k-1}){V^{k-1}}^{\top}$. Note that $k_{s}$ is non-increasing in this procedure. However, if $k_{s}$ is too small at some iteration, the non-expansive property (3.1) of the shrinkage operator $S_{\nu}$ may be violated since the approximate SVD is not a valid approximation when $k_{s}$ is too small. Thus, in algorithm FPCA, if (3.1) is violated 10 times, we increase $k_{s}$ by $1$. Our numerical experience indicates that this technique makes our algorithm very robust.

Our numerical results in Section 7 show that this approximate SVD based FPC algorithm: FPCA, is very fast, robust, and significantly outperforms other solvers (such as SDPT3) in recovering low-rank matrices. This result is not surprising. One reason for this is that in the approximate SVD algorithm, we compute a low-rank approximation to the original matrix. Hence, the iterative matrices produced by our algorithm are more likely to be of low-rank than an exact solution to the nuclear norm minimization problem (1.2), or equivalently, to the SDP (1.2), which is exactly what we want. Some convergence/recoverability properties of a variant of FPCA, which uses a truncated SVD rather than a randomized SVD at each step, are discussed in (Goldfarb-Ma-2009).

Numerical results

In this section, we report on the application of our FPC, FPCA and Bregman iterative algorithms to a series of matrix completion problems of the form (1) to demonstrate the ability of these algorithms to efficiently recover low-rank matrices.

To illustrate the performance of our algorithmic approach combined with exact and approximate SVD algorithms, different stopping rules, and with or without debiasing, we tested the following solvers.

FPC1. Exact SVD, no debiasing, stopping rule: (4.2).

FPC2. Exact SVD, no debiasing, stopping rule: (4.1) and (4.2).

FPC3. Exact SVD with debiasing, stopping rule: (4.2).

FPCA. Approximate SVD, no debiasing, stopping rule: (4.2).

Bregman. Bregman iterative method using FPC2 to solve the subproblems.

to estimate the closeness of $X_{opt}$ to $M$, where $X_{opt}$ is the “optimal” solution to (1.2) produced by our algorithms. We declared $M$ to be recovered if the relative error was less than $10^{-3}$, which is the criterion used in (Recht-Fazel-Parrilo-2007) and (Candes-Recht-2008). We use $RA,RU,RL$ to denote the average, largest and smallest relative error of the successfully recovered matrices, respectively.

We summarize the parameter settings used by the algorithms in Table 1. We use $I_{m}$ to denote the maximum number of iterations allowed for solving each subproblem in FPC, i.e., if the stopping rules (4.2) (and (4.1)) are not satisfied after $I_{m}$ iterations, we terminate the subproblem and decrease $\mu$ to start the next subproblem.

All numerical experiments were run in MATLAB 7.3.0 on a Dell Precision 670 workstation with an Intel Xeon(TM) 3.4GHZ CPU and 6GB of RAM.

The comparisons between FPC1, FPC2, FPC3 and SDPT3 for small matrix completion problems are presented in Table 2. From Table 2 we can see that FPC1 and FPC2 achieve almost the same recoverability and relative error, which means that as long as we set $xtol$ to be very small (like $10^{-10}$ ), we only need to use (4.2) as the stopping rule for the inner iterations. That is, use of stopping rule (4.1) does not affect the performance of the algorithm. Of course FPC2 costs more time than FPC1 since more iterations are sometimes needed to satisfy the stopping rules in FPC2. While FPC3 can improve the recoverability, it costs more time for performing debiasing. SDPT3 seems to obtain more accurate solutions than FPC1, FPC2 or FPC3.

To illustrate the performance of our Bregman iterative algorithm, we compare the results of using it versus using FPC2 in Table 3. From our numerical experience, for those problems for which the Bregman iterative algorithm greatly improves the recoverability, the Bregman iterative algorithm usually takes 2 to 3 iterations. Thus, in our numerical tests, we fixed the number of subproblems solved by our Bregman algorithm to 3. Since our Bregman algorithm achieves as good a relative error as the FPC algorithm, we only report how many of the examples that are successfully recovered by FPC, are improved greatly by using our Bregman iterative algorithm. In Table 3, NIM is the number of examples that the Bregman iterative algorithm outperformed FPC2 greatly (the relative errors obtained from FPC2 were at least $10^{4}$ times larger than those obtained by the Bregman algorithm). From Table 3 we can see that for more than half of the examples successfully recovered by FPC2, the Bregman iterative algorithm improved the relative errors greatly (from [$10^{-10}$, $10^{-9}$] to [$10^{-16}$, $10^{-15}$]). Of course the run times for the Bregman iterative algorithm were about three times that for algorithm FPC2, since the former calls the latter three times to solve the subproblems.

In the following, we discuss the numerical results obtained by our approximate SVD based FPC algorithm (FPCA). We will see from these numerical results that FPCA achieves much better recoverability and is much faster than any of the solvers FPC1, FPC2, FPC3 or SDPT3.

We present the numerical results of FPCA for small (m=n=40) and medium (m=n=100) problems in Tables 4, and 5 respectively. Since we found that $xtol=10^{-6}$ is small enough to guarantee very good recoverability, we set $xtol=10^{-6}$ in algorithm FPCA and used only (4.2) as stopping rule for the inner iterations. From these tables, we can see that our FPCA algorithm is much more powerful than SDPT3 for randomly created matrix completion problems. When $m=n=40$ and $p=800$, and the rank $r$ was less than or equal to 8, FPCA recovered the matrices in all 50 examples. When rank $r=9$, it failed on only one example. Even for rank $r=10$, which is almost the largest rank that satisfies $FR\leq 1$, FPCA still recovered the solution in more than $60\%$ of the examples. However, SDPT3 started to fail to recover the matrices when the rank $r=2$. When $r=6$, there was only one example out of 50 where the correct solution matrix was recovered. When $r\geq 7$, none of the 50 examples could be recovered. For the medium sized matrices $(m=n=100)$ we used $p=2000$, which is only a $20\%$ measurement rate, FPCA recovered the matrices in all 50 examples when $r\leq 6$. For $r=7,$ FPCA recovered the matrices in most of the examples (49 out of 50). When $r=8,$ more than $60\%$ of the matrices were recovered successfully by FPCA. Even when $r=9,$ FPCA still recovered 1 matrices. However, SDPT3 could not recover all of the matrices even when the rank $r=1$ and none of the matrices were recovered when $r\geq 4.$ When we increased the number of measurements to $3000$, we recovered the matrices in all 50 examples up to rank $r=12.$ When $r=13,14,$ we still recovered most of them. However, SDPT3 started to fail for some matrices when $r=3.$ When $r\geq 8$, SDPT3 failed to recover any of the matrices. We can also see that for the medium sized problems, FPCA was much faster than SDPT3.

2 Comparison of FPCA and SVT

In this subsection we compare our FPCA algorithm against the SVT algorithm proposed in (Cai-Candes-Shen-2008). The SVT code is downloaded from http://svt.caltech.edu. We constructed two sets of test problems. One set contained “easy” problems. These problems are “easy” because the matrices are of very low-rank compared to the matrix size and the number of samples, and hence they are easy to recover. For all problems in this set, $FR$ was less than 0.34. The other set contained “hard” problems, i.e., problems that are very challenging. These problems involved matrices that are not of very low rank and for which sampled a very limited number of entries. For this set of problems, $FR$ ranged from 0.40 to 0.87. The maximum iteration number in SVT was set to be 1000. All other parameters were set to their default values in SVT. The parameters of FPCA were set somewhat loosely for easy problems. Specifically, we set $\bar{\mu}=10^{-4},xtol=10^{-4},\tau=2,I_{m}=10$ and all other parameters were set to the values given in Table 1. Relative errors and times were averaged over 5 runs. In this subsection, all test matrices were square, i.e., $m=n.$

From Table 6, we can see that for the easy problems except for one problem which is exceptionally sparse as well as having low rank, FPCA was much faster and usually provided more accurate solutions than SVT.

For hard problems, all parameters of FPCA were set to the values given in Table 1, except that we set $xtol=10^{-6}$ since this value is small enough to guarantee very good recoverability. Also, for small problems ( i.e., $\max\{m,n\}<1000$ ), we set $I_{m}=500$; and for large problems ( i.e., $\max\{m,n\}\geq 1000$ ), we set $I_{m}=20.$ We use “—” to indicate that the algorithm either diverges or does not terminate in one hour. Relative errors and times were averaged over 5 runs.

From Table 7, we can see that for the hard problems, SVT either diverged or did not solve the problems in less than one hour, or it yielded a very inaccurate solution. In contrast, FPCA always provided a very good solution efficiently.

We can also see that FPCA was able to efficiently solve large problems ($m=n=1000$) that could not be solved by SDPT3 due to the large size of the matrices and the large number of constraints.

3 Results for real data matrices

In this section, we consider matrix completion problems based on two real data sets: the Jester joke data set (Goldberg-Roeder-Gupta-Perkins-2001) and the DNA data set (Spellman-1998). The Jester joke data set contains 4.1 million ratings for 100 jokes from 73,421 users and is available on the website http://www.ieor.berkeley.edu/˜Egoldberg/jester-data/. Since the number of jokes is only 100, but the number of users is quite large, we randomly selected $n_{u}$ users to get a modestly sized matrix for testing purpose. As in (Srebro-Jaakkola-2003), we randomly held out two ratings for each user. Since some entries in the matrix are missing, we cannot compute the relative error as we did for the randomly created matrices. Instead, we computed the Normalized Mean Absolute Error (NMAE) as in (Goldberg-Roeder-Gupta-Perkins-2001) and (Srebro-Jaakkola-2003). The Mean Absolute Error (MAE) is defined as

where $r_{j}^{i}$ and $\hat{r}_{j}^{i}$ are the withheld and predicted ratings of movie $j$ by user $i$, respectively, for $j=i_{1},i_{2}.$ NMAE is defined as

where $r_{\min}$ and $r_{\max}$ are lower and upper bounds for the ratings. Since all ratings are scaled to the range $[-10,+10]$, we have $r_{\min}=-10$ and $r_{\max}=10.$

The numerical results for the Jester data set using FPC1 and FPCA are presented in Tables 8 and 9, respectively. In these two tables, $\sigma_{\max}$ and $\sigma_{\min}$ are the largest and smallest positive singular values of the recovered matrices, and $rank$ is the rank of the recovered matrices. The distributions of the singular values of the recovered matrices are shown in Figures 1 and 2. From Tables 8 and 9 we can see that by using FPC1 and FPCA to recover these matrices, we can get relatively low NMAEs, which are comparable to the results shown in (Srebro-Jaakkola-2003) and (Goldberg-Roeder-Gupta-Perkins-2001).

We also used two data sets of DNA microarrays from (Spellman-1998). These data sets are available on the website http://cellcycle-www.stanford.edu/. The first microarray data set is a matrix that represents the expression of 6178 genes in 14 experiments based on the Elutriation data set in (Spellman-1998). The second microarray data set is based on the Cdc15 data set in (Spellman-1998), and represents the expression of 6178 genes in 24 experiments. However, some entries in these two matrices are missing. For evaluating our algorithms, we created complete matrices by deleted all rows containing missing values. This is similar to how the DNA microarray data set was preprocessed in (Troyanskaya-2001). The resulting complete matrix for the Elutriation data set was $5766\times 14$. The complete matrix for the Cdc15 data set was $4381\times 24$. We must point out that these DNA microarray matrices are neither low-rank nor even approximately low-rank although such claims have been made in some papers. The distributions of the singular values of these two matrices are shown in Figure 3. From this figure we can see that in each microarray matrix, only one singular value is close to zero, while the others are far away from zero. Thus there is no way to claim that the rank of the Elutriation matrix is less than 13, or the rank of the Cdc15 matrix is less than 23. Since these matrices are not low-rank, we cannot expect our algorithms to recover these matrices by sampling only a small portion of their entries. Thus we needed to further modify the data sets to yield low-rank matrices. Specifically, we used the best rank-2 approximation to the Elutriation matrix as the new complete Elutriation matrix and the best rank-5 approximation to the Cdc15 matrix as the new complete Cdc15 matrix. The numerical results for FPCA for recovering these two matrices are presented in Table 10. In the FPCA algorithm, we set $\epsilon_{k_{s}}=10^{-2}$ and $xtol=10^{-6}$. For the Elutriation matrix, we set $c_{s}=115$ and for the Cdc15 matrix, we set $c_{s}=88$. The observed entries were randomly sampled. From Table 10 we can see that by taking 60% of the entries of the matrices, our FPCA algorithm can recover these matrices very well, yielding relative errors as low as $10^{-5}$ and $10^{-6}$, which is promising for practical use.

To graphically illustrate the effectiveness of FPCA, we applied it to image inpainting (Bertalmio-Sapiro-Caselles-Ballester-2000). Grayscale images and color images can be expressed as matrices and tensors, respectively. In grayscale image inpainting, the grayscale value of some of the pixels of the image are missing, and we want to fill in these missing values. If the image is of low-rank, or of numerical low-rank, we can solve the image inpainting problem as a matrix completion problem (1). In our test we applied SVD to the $512\times 512$ image in Figure 4(a), and truncated this decomposition to get the rank-40 image which is shown in Figure 4(b). Figure 4(c) is a masked version of the image in Figure 4(a), where one half of the pixels in Figure 4(a) were masked uniformly at random. Figure 4(d) is the image obtained from Figure 4(c) by applying FPCA. Figure 4(d) is a low-rank approximation to Figure 4(a) with a relative error of $8.41e-2.$ Figure 4(e) is a masked version of the image in Figure 4(b), where one half of the pixels in Figure 4(b) were masked uniformly at random. Figure 4(f) is the image obtained from Figure 4(e) by applying FPCA. Figure 4(f) is an approximation to Figure 4(b) with a relative error of $3.61e-2.$ Figure 4(g) is another masked image obtained from Figure 4(b), where 4 percent of the pixels were masked in a non-random fashion. Figure 4(h) is the image obtained from Figure 4(g) by applying FPCA. Figure 4(g) is an approximation to Figure 4(b) with a relative error of $1.70e-2$.

Conclusions and discussions

In this paper, we derived a fixed point continuation algorithm and a Bregman iterative algorithm for solving the linearly constrained nuclear norm minimization problem, which is a convex relaxation of the NP-hard linearly constrained matrix rank minimization problem. The convergence of the fixed point iterative scheme was established. By adopting an approximate SVD technique, we obtained a very powerful algorithm (FPCA) for the matrix rank minimization problem. On matrix completion problems, FPCA greatly outperforms SDP solvers such as SDPT3 in both speed and recoverability of low-rank matrices. Further study is needed to prove the convergence of algorithm FPCA.

References