A Singular Value Thresholding Algorithm for Matrix Completion

Introduction

There is a rapidly growing interest in the recovery of an unknown low-rank or approximately low-rank matrix from very limited information. This problem occurs in many areas of engineering and applied science such as machine learning , control and computer vision, see . As a motivating example, consider the problem of recovering a data matrix from a sampling of its entries. This routinely comes up whenever one collects partially filled out surveys, and one would like to infer the many missing entries. In the area of recommender systems, users submit ratings on a subset of entries in a database, and the vendor provides recommendations based on the user’s preferences. Because users only rate a few items, one would like to infer their preference for unrated items; this is the famous Netflix problem . Recovering a rectangular matrix from a sampling of its entries is known as the matrix completion problem. The issue is of course that this problem is extraordinarily ill posed since with fewer samples than entries, we have infinitely many completions. Therefore, it is apparently impossible to identify which of these candidate solutions is indeed the “correct” one without some additional information.

In many instances, however, the matrix we wish to recover has low rank or approximately low rank. For instance, the Netflix data matrix of all user-ratings may be approximately low-rank because it is commonly believed that only a few factors contribute to anyone’s taste or preference. In computer vision, inferring scene geometry and camera motion from a sequence of images is a well-studied problem known as the structure-from-motion problem. This is an ill-conditioned problem for objects may be distant with respect to their size, or especially for “missing data” which occur because of occlusion or tracking failures. However, when properly stacked and indexed, these images form a matrix which has very low rank (e.g. rank 3 under orthography) . Other examples of low-rank matrix fitting abound; e.g. in control (system identification), machine learning (multi-class learning) and so on. Having said this, the premise that the unknown has (approximately) low rank radically changes the problem, making the search for solutions feasible since the lowest-rank solution now tends to be the right one.

provided that the number of samples obeys

for some positive numerical constant $C$.Note that an $n\times n$ matrix of rank $r$ depends upon $r(2n-r)$ degrees of freedom. In (1.1), the functional $\|\bm{X}\|_{*}$ is the nuclear norm of the matrix $\bm{M}$, which is the sum of its singular values. The optimization problem (1.1) is convex and can be recast as a semidefinite program . In some sense, this is the tightest convex relaxation of the NP-hard rank minimization problem

since the nuclear ball $\{\bm{X}:\|\bm{X}\|_{*}\leq 1\}$ is the convex hull of the set of rank-one matrices with spectral norm bounded by one. Another interpretation of Candès and Recht’s result is that under suitable conditions, the rank minimization program (1.3) and the convex program (1.1) are formally equivalent in the sense that they have exactly the same unique solution.

2 Algorithm outline

Because minimizing the nuclear norm both provably recovers the lowest-rank matrix subject to constraints (see for related results) and gives generally good empirical results in a variety of situations, it is understandably of great interest to develop numerical methods for solving (1.1). In , this optimization problem was solved using one of the most advanced semidefinite programming solvers, namely, SDPT3 . This solver and others like SeDuMi are based on interior-point methods, and are problematic when the size of the matrix is large because they need to solve huge systems of linear equations to compute the Newton direction. In fact, SDPT3 can only handle $n\times n$ matrices with $n\leq 100$. Presumably, one could resort to iterative solvers such as the method of conjugate gradients to solve for the Newton step but this is problematic as well since it is well known that the condition number of the Newton system increases rapidly as one gets closer to the solution. In addition, none of these general purpose solvers use the fact that the solution may have low rank. We refer the reader to for some recent progress on interior-point methods concerning some special nuclear norm-minimization problems.

This paper develops the singular value thresholding algorithm for approximately solving the nuclear norm minimization problem (1.1) and by extension, problems of the form

until a stopping criterion is reached. In (1.6), $\text{shrink}(\bm{Y},\tau)$ is a nonlinear function which applies a soft-thresholding rule at level $\tau$ to the singular values of the input matrix, see Section 2 for details. The key property here is that for large values of $\tau$, the sequence $\{\bm{X}^{k}\}$ converges to a solution which very nearly minimizes (1.5). Hence, at each step, one only needs to compute at most one singular value decomposition and perform a few elementary matrix additions. Two important remarks are in order:

Sparsity. For each $k\geq 0$, $\bm{Y}^{k}$ vanishes outside of $\Omega$ and is, therefore, sparse, a fact which can be used to evaluate the shrink function rapidly.

Low-rank property. The matrices $\bm{X}^{k}$ turn out to have low rank, and hence the algorithm has minimum storage requirement since we only need to keep principal factors in memory.

Our numerical experiments demonstrate that the proposed algorithm can solve problems, in Matlab, involving matrices of size $30,000\times 30,000$ having close to a billion unknowns in 17 minutes on a standard desktop computer with a 1.86 GHz CPU (dual core with Matlab’s multithreading option enabled) and 3 GB of memory. As a consequence, the singular value thresholding algorithm may become a rather powerful computational tool for large scale matrix completion.

3 General formulation

The singular value thresholding algorithm can be adapted to deal with other types of convex constraints. For instance, it may address problems of the form

where each $f_{i}$ is a Lipschitz convex function (note that one can handle linear equality constraints by considering pairs of affine functionals). In the simpler case where the $f_{i}$’s are affine functionals, the general algorithm goes through a sequence of iterations which greatly resemble (1.6). This is useful because this enables the development of numerical algorithms which are effective for recovering matrices from a small subset of sampled entries possibly contaminated with noise.

4 Contents and notations

The rest of the paper is organized as follows. In Section 2, we derive the singular value thresholding (SVT) algorithm for the matrix completion problem, and recasts it in terms of a well-known Lagrange multiplier algorithm. In Section 3, we extend the SVT algorithm and formulate a general iteration which is applicable to general convex constraints. In Section 4, we establish the convergence results for the iterations given in Sections 2 and 3. We demonstrate the performance and effectiveness of the algorithm through numerical examples in Section 5, and review additional implementation details. Finally, we conclude the paper with a short discussion in Section 6.

Before continuing, we provide here a brief summary of the notations used throughout the paper. Matrices are bold capital, vectors are bold lowercase and scalars or entries are not bold. For instance, $\bm{X}$ is a matrix and $X_{ij}$ its $(i,j)$th entry. Likewise, $\bm{x}$ is a vector and $x_{i}$ its $i$th component. The nuclear norm of a matrix is denoted by $\|\bm{X}\|_{*}$, the Frobenius norm by \|\mbox{\boldmathX}\|_{F} and the spectral norm by \|\mbox{\boldmathX}\|_{2}; note that these are respectively the 1-norm, the 2-norm and the sup-norm of the vector of singular values. The adjoint of a matrix $\bm{X}$ is $\bm{X}^{*}$ and similarly for vectors. The notation \operatorname{diag}(\mbox{\boldmathx}), where $\bm{x}$ is a vector, stands for the diagonal matrix with $\{x_{i}\}$ as diagonal elements. We denote by $\langle\bm{X},\bm{Y}\rangle=\operatorname{trace}(\bm{X}^{*}\bm{Y})$ the standard inner product between two matrices ($\|\bm{X}\|_{F}^{2}=\langle\bm{X},\bm{X}\rangle$). The Cauchy-Schwarz inequality gives \langle\mbox{\boldmathX},\mbox{\boldmathY}\rangle\leq\|\mbox{\boldmathX}\|_{F}\|\mbox{\boldmathY}\|_{F} and it is well known that we also have \langle\mbox{\boldmathX},\mbox{\boldmathY}\rangle\leq\|\mbox{\boldmathX}\|_{*}\|\mbox{\boldmathY}\|_{2} (the spectral and nuclear norms are dual from one another), see e.g. .

The Singular Value Thresholding Algorithm

This section introduces the singular value thresholding algorithm and discusses some of its basic properties. We begin with the definition of a key building block, namely, the singular value thresholding operator.

where $\bm{U}$ and $\bm{V}$ are respectively $n_{1}\times r$ and $n_{2}\times r$ matrices with orthonormal columns, and the singular values $\sigma_{i}$ are positive (unless specified otherwise, we will always assume that the SVD of a matrix is given in the reduced form above). For each $\tau\geq 0$, we introduce the soft-thresholding operator $\mathcal{D}_{\tau}$ defined as follows:

where $t_{+}$ is the positive part of $t$, namely, $t_{+}=\max(0,t)$. In words, this operator simply applies a soft-thresholding rule to the singular values of $\bm{X}$, effectively shrinking these towards zero. This is the reason why we will also refer to this transformation as the singular value shrinkage operator. Even though the SVD may not be unique, it is easy to see that the singular value shrinkage operator is well defined and we do not elaborate further on this issue. In some sense, this shrinkage operator is a straightforward extension of the soft-thresholding rule for scalars and vectors. In particular, note that if many of the singular values of $\bm{X}$ are below the threshold $\tau$, the rank of $\mathcal{D}_{\tau}(\bm{X})$ may be considerably lower than that of $\bm{X}$, just like the soft-thresholding rule applied to vectors leads to sparser outputs whenever some entries of the input are below threshold.

The singular value thresholding operator is the proximity operator associated with the nuclear norm. Details about the proximity operator can be found in e.g. .

for all $\bm{X}$. Now $\hat{\bm{X}}$ minimizes $h_{0}$ if and only if is a subgradient of the functional $h_{0}$ at the point $\hat{\bm{X}}$, i.e.

Set $\hat{\bm{X}}:={\cal D}_{\tau}(\bm{Y})$ for short. In order to show that $\hat{\bm{X}}$ obeys (2.5), decompose the SVD of $\bm{Y}$ as

where $\bm{U}_{0}$, $\bm{V}_{0}$ (resp. $\bm{U}_{1}$, $\bm{V}_{1}$) are the singular vectors associated with singular values greater than $\tau$ (resp. smaller than or equal to $\tau$). With these notations, we have

By definition, $\bm{U}_{0}^{*}\bm{W}=0$, $\bm{W}\bm{V}_{0}=0$ and since the diagonal elements of $\bm{\Sigma}_{1}$ have magnitudes bounded by $\tau$, we also have $\|\bm{W}\|_{2}\leq 1$. Hence $\bm{Y}-\hat{\bm{X}}\in\tau\partial\|\hat{\bm{X}}\|_{*}$, which concludes the proof.

2 Shrinkage iterations

We are now in the position to introduce the singular value thresholding algorithm. Fix $\tau>0$ and a sequence $\{\delta_{k}\}$ of positive step sizes. Starting with $\bm{Y}_{0}$, inductively define for $k=1,2,\ldots$,

until a stopping criterion is reached (we postpone the discussion this stopping criterion and of the choice of step sizes). This shrinkage iteration is very simple to implement. At each step, we only need to compute an SVD and perform elementary matrix operations. With the help of a standard numerical linear algebra package, the whole algorithm can be coded in just a few lines.

Before addressing further computational issues, we would like to make explicit the relationship between this iteration and the original problem (1.1). In Section 4, we will show that the sequence $\{\bm{X}^{k}\}$ converges to the unique solution of an optimization problem closely related to (1.1), namely,

Furthermore, it is intuitive that the solution to this modified problem converges to that of (1.5) as $\tau\to\infty$ as shown in Section 3. Thus by selecting a large value of the parameter $\tau$, the sequence of iterates converges to a matrix which nearly minimizes (1.1).

As mentioned earlier, there are two crucial properties which make this algorithm ideally suited for matrix completion.

Low-rank property. A remarkable empirical fact is that the matrices in the sequence $\{\bm{X}^{k}\}$ have low rank (provided, of course, that the solution to (2.8) has low rank). We use the word “empirical” because all of our numerical experiments have produced low-rank sequences but we cannot rigorously prove that this is true in general. The reason for this phenomenon is, however, simple: because we are interested in large values of $\tau$ (as to better approximate the solution to (1.1)), the thresholding step happens to ‘kill’ most of the small singular values and produces a low-rank output. In fact, our numerical results show that the rank of $\bm{X}^{k}$ is nondecreasing with $k$, and the maximum rank is reached in the last steps of the algorithm, see Section 5.

Thus, when the rank of the solution is substantially smaller than either dimension of the matrix, the storage requirement is low since we could store each $\bm{X^{k}}$ in its SVD form (note that we only need to keep the current iterate and may discard earlier values).

Sparsity. Another important property of the SVT algorithm is that the iteration matrix $\bm{Y}^{k}$ is sparse. Since $\bm{Y}^{0}=\bm{0}$, we have by induction that $\bm{Y}^{k}$ vanishes outside of $\Omega$. The fewer entries available, the sparser $\bm{Y}^{k}$. Because the sparsity pattern $\Omega$ is fixed throughout, one can then apply sparse matrix techniques to save storage. Also, if $|\Omega|=m$, the computational cost of updating $\bm{Y}^{k}$ is of order $m$. Moreover, we can call subroutines supporting sparse matrix computations, which can further reduce computational costs.

One such subroutine is the SVD. However, note that we do not need to compute the entire SVD of $\bm{Y}^{k}$ to apply the singular value thresholding operator. Only the part corresponding to singular values greater than $\tau$ is needed. Hence, a good strategy is to apply the iterative Lanczos algorithm to compute the first few singular values and singular vectors. Because $\bm{Y}^{k}$ is sparse, $\bm{Y}^{k}$ can be applied to arbitrary vectors rapidly, and this procedure offers a considerable speedup over naive methods.

3 Relation with other works

Finally, we would like to contrast the SVT iteration (2.7) with the popular iterative soft-thresholding algorithm used in many papers in imaging processing and perhaps best known under the name of Proximal Forward-Backward Splitting method (PFBS), see for example. The constrained minimization problem (1.5) may be relaxed into

for some $\lambda>0$. Theorem 2.1 asserts that $\mathcal{D}_{\lambda}$ is the proximity operator of $\lambda\|\bm{X}\|_{*}$ and Proposition 3.1(iii) in gives that the solution to this unconstrained problem is characterized by the fixed point equation $\bm{X}=\mathcal{D}_{\lambda\delta}(\bm{X}+\delta P_{\Omega}(\bm{M}-\bm{X}))$ for each $\delta>0$. One can then apply a simplified version of the PFBS method (see (3.6) in ) to obtain iterations of the form

Introducing an intermediate matrix $\bm{Y}^{k}$, this algorithm may be expressed as

The difference with (2.7) may seem subtle at first—replacing $\bm{X}^{k}$ in (2.10) with $\bm{Y}^{k-1}$ and setting $\delta_{k}=\delta$ gives (2.7) with $\tau=\lambda\delta$—but has enormous consequences as this gives entirely different algorithms. First, they have different limits: while (2.7) converges to the solution of the constrained minimization (2.8), (2.10) converges to the solution of (2.9) provided that the sequence of step sizes is appropriately selected. Second, selecting a large $\lambda$ (or a large value of $\tau=\lambda\delta$) in (2.10) gives a low-rank sequence of iterates and a limit with small nuclear norm. The limit, however, does not fit the data and this is why one has to choose a small or moderate value of $\lambda$ (or of $\tau=\lambda\delta$). However, when $\lambda$ is not sufficiently large, the $\bm{X}^{k}$ may not have low rank even though the solution has low rank (and one may need to compute many singular vectors), and $\bm{Y}^{k}$ is not sufficiently sparse to make the algorithm computationally attractive. Moreover, the limit does not necessary have a small nuclear norm. These are reasons why (2.10) is not suitable for matrix completion.

4 Interpretation as a Lagrange multiplier method

In this section, we recast the SVT algorithm as a type of Lagrange multiplier algorithm known as Uzawa’s algorithm. An important consequence is that this will allow us to extend the SVT algorithm to other problems involving the minimization of the nuclear norm under convex constraints, see Section 3. Further, another contribution of this paper is that this framework actually recasts linear Bregman iterations as a very special form of Uzawa’s algorithm, hence providing fresh and clear insights about these iterations.

In what follows, we set $f_{\tau}(\bm{X})=\tau\|\bm{X}\|_{*}+\frac{1}{2}\|\bm{X}\|_{F}^{2}$ for some fixed $\tau>0$ and recall that we wish to solve (2.8)

The Lagrangian for this problem is given by

(The function $g_{0}(\bm{Y})=\inf_{\bm{X}}{\cal L}(\bm{X},\bm{Y})$ is called the dual function.) Uzawa’s algorithm approaches the problem of finding a saddlepoint with an iterative procedure. From $\bm{Y}_{0}=\bm{0}$, say, inductively define

where $\{\delta_{k}\}_{k\geq 1}$ is a sequence of positive step sizes. Uzawa’s algorithm is, in fact, a subgradient method applied to the dual problem, where each step moves the current iterate in the direction of the gradient or of a subgradient. Indeed, observe that

It remains to compute the minimizer of the Lagrangian (2.12), and note that

However, we know that the minimizer is given by $\mathcal{D}_{\tau}(\mathcal{P}_{\Omega}(\bm{Y}))$ and since $\bm{Y}^{k}=\mathcal{P}_{\Omega}(\bm{Y}^{k})$ for all $k\geq 0$, Uzawa’s algorithm takes the form

General Formulation

This section presents a general formulation of the SVT algorithm for approximately minimizing the nuclear norm of a matrix under convex constraints.

Set the objective functional $f_{\tau}(\bm{X})=\tau\|\bm{X}\|_{*}+\frac{1}{2}\|\bm{X}\|_{F}^{2}$ for some fixed $\tau>0$, and consider the following optimization problem:

The iteration (3.3) is of course the same as (2.7) in the case where $\mathcal{A}$ is a sampling operator extracting $m$ entries with indices in $\Omega$ out of an $n_{1}\times n_{2}$ matrix. To verify this claim, observe that in this situation, $\mathcal{A}^{*}\mathcal{A}=\mathcal{P}_{\Omega}$, and let $\bm{M}$ be any matrix obeying $\mathcal{A}(\bm{M})=\bm{b}$. Then defining $\bm{Y}^{k}=\mathcal{A}^{*}(\bm{y}^{k})$ and substituting this expression in (3.3) gives (2.7).

2 General convex constraints

One can also adapt the algorithm to handle general convex constraints. Suppose we wish to minimize $f_{\tau}(\bm{X})$ defined as before over a convex set $\bm{X}\in\mathcal{C}$. To simplify, we will assume that this convex set is given by

where the $f_{i}$’s are convex functionals (note that one can handle linear equality constraints by considering pairs of affine functionals). The problem of interest is then of the form

Just as before, it is intuitive that as $\tau\to\infty$, the solution to this problem converges to a minimizer of the nuclear norm under the same constraints (1.7) as shown in Theorem 3.1 at the end of this section.

Put $\mathcal{F}(\bm{X}):=(f_{1}(\bm{X}),\ldots,f_{m}(\bm{X}))$ for short. Then the Lagrangian for (3.4) is equal to

Above, $\bm{x}_{+}$ is of course the vector with entries equal to $\max(x_{i},0)$. When $\mathcal{F}$ is an affine mapping of the form $\bm{b}-\mathcal{A}(\bm{X})$ so that one solves

and thus the extension to linear inequality constraints is straightforward.

3 Example

An interesting example concerns the extension of the Dantzig selector to matrix problems. Suppose we have available linear measurements about a matrix $\bm{M}$ of interest

where $\bm{E}$ is an array of tolerances, which is adjusted to fit the noise statistics . Above, $\text{\bf vec}(\bm{A})\leq\text{\bf vec}(\bm{B})$, for any two matrices $\bm{A}$ and $\bm{B}$, means componentwise inequalities; that is, $A_{ij}\leq B_{ij}$ for all indices $i,j$. We use this notation as not to confuse the reader with the positive semidefinite ordering. In the case of the matrix completion problem where $\mathcal{A}$ extracts sampled entries indexed by $\Omega$, one can always see the data vector as the sampled entries of some matrix $\bm{B}$ obeying $\mathcal{P}_{\Omega}(\bm{B})=\mathcal{A}^{*}(\bm{b})$; the constraint is then natural for it may be expressed as

If $\bm{z}$ is white noise with standard deviation $\sigma$, one may want to use a multiple of $\sigma$ for $E_{ij}$. In words, we are looking for a matrix with minimum nuclear norm under the constraint that all of its sampled entries do not deviate too much from what has been observed.

where again $[\cdot]_{+}$ is applied componentwise.

We conclude by noting that in the matrix completion problem where $\mathcal{A}^{*}\mathcal{A}=\mathcal{P}_{\Omega}$ and one observes $\mathcal{P}_{\Omega}(\bm{B})$, one can check that this iteration simplifies to

Again, this is easy to implement and whenever the solution has low rank, the iterates $\bm{X}^{k}$ have low rank as well.

4 When the proximal problem gets close

We now show that minimizing the proximal objective $f_{\tau}(\bm{X})=\tau\|\bm{X}\|_{*}+\frac{1}{2}\|\bm{X}\|_{F}^{2}$ is the same as minimizing the nuclear norm in the limit of large $\tau$’s. The theorem below is general and covers the special case of linear equality constraints as in (2.8).

Let $\bm{X}_{\tau}^{\star}$ be the solution to (3.4) and $\bm{X}_{\infty}$ be the minimum Frobenius-norm solution to (1.7) defined as

Assume that the $f_{i}(\bm{X})$’s, $1\leq i\leq m$, are convex and lower semi-continuous. Then

Proof. It follows from the definition of $\bm{X}_{\tau}^{\star}$ and $\bm{X}_{\infty}$ that

which implies that $\|\bm{X}_{\tau}^{\star}\|_{F}^{2}$ is bounded uniformly in $\tau$. Thus, we would prove the theorem if we could establish that any convergent subsequence $\{\bm{X}^{\star}_{\tau_{k}}\}_{k\geq 1}$ must converge to $\bm{X}_{\infty}$.

Consider an arbitrary converging subsequence $\{\bm{X}^{\star}_{\tau_{k}}\}$ and set $\bm{X}_{c}:=\lim_{k\rightarrow\infty}\bm{X}^{\star}_{\tau_{k}}$. Since for each $1\leq i\leq m$, $f_{i}(\bm{X}^{\star}_{\tau_{k}})\leq 0$ and $f_{i}$ is lower semi-continuous, $\bm{X}_{c}$ obeys

Furthermore, since $\|\bm{X}_{\tau}^{\star}\|_{F}^{2}$ is bounded, (3.13) yields

An immediate consequence is $\lim_{\tau\to\infty}\|\bm{X}_{\tau}^{\star}\|_{*}=\|\bm{X}_{\infty}\|_{*}$ and, therefore, $\|\bm{X}_{c}\|_{*}=\|\bm{X}_{\infty}\|_{*}$. This shows that $\bm{X}_{c}$ is a solution to (1.1). Now it follows from the definition of $\bm{X}_{\infty}$ that $\|\bm{X}_{c}\|_{F}\geq\|\bm{X}_{\infty}\|_{F}$, while we also have $\|\bm{X}_{c}\|_{F}\leq\|\bm{X}_{\infty}\|_{F}$ because of (3.14). We conclude that $\|\bm{X}_{c}\|_{F}=\|\bm{X}_{\infty}\|_{F}$ and thus $\bm{X}_{c}=\bm{X}_{\infty}$ since $\bm{X}_{\infty}$ is unique.

Convergence Analysis

This section establishes the convergence of the SVT iterations. We begin with the simpler proof of the convergence of (2.7) in the special case of the matrix completion problem, and then present the argument for the more general constraints (3.5). We hope that this progression will make the second and more general proof more transparent.

We begin by recording a lemma which establishes the strong convexity of the objective $f_{\tau}$.

Let $\bm{Z}\in\partial f_{\tau}(\bm{X})$ and $\bm{Z}^{\prime}\in\partial f_{\tau}(\bm{X}^{\prime})$. Then

Proof. An element $\bm{Z}$ of $\partial f_{\tau}(\bm{X})$ is of the form $\bm{Z}=\tau\bm{Z}_{0}+\bm{X}$, where $\bm{Z}_{0}\in\partial\|\bm{X}\|_{*}$, and similarly for $\bm{Z}^{\prime}$. This gives

and it thus suffices to show that the first term of the right-hand side is nonnegative. From (2.6), we have that any subgradient of the nuclear norm at $\bm{X}$ obeys $\|\bm{Z}_{0}\|_{2}\leq 1$ and $\langle\bm{Z}_{0},\bm{X}\rangle=\|\bm{X}\|_{*}$. In particular, this gives

This lemma is key in showing that the SVT algorithm (2.7) converges.

Suppose that the sequence of step sizes obeys $0<\inf\delta_{k}\leq\sup\delta_{k}<2$. Then the sequence $\{\bm{X}^{k}\}$ obtained via (2.7) converges to the unique solution of (2.8).

Proof. Let $(\bm{X}^{\star},\bm{Y}^{\star})$ be primal-dual optimal for the problem (2.8). The optimality conditions give

for some $\bm{Z}^{k}\in\partial f_{\tau}(\bm{X}^{k})$ and some $\bm{Z}^{\star}\in\partial f_{\tau}(\bm{X}^{\star})$. We then deduce that

and, therefore, it follows from Lemma 4.1 that

We continue and observe that because $\mathcal{P}_{\Omega}\bm{X}^{\star}=\mathcal{P}_{\Omega}\bm{M}$,

Therefore, setting $r_{k}=\|\mathcal{P}_{\Omega}(\bm{Y}^{k}-\bm{Y}^{\star})\|_{F}$,

since for any matrix $\bm{X}$, $\|\mathcal{P}_{\Omega}(\bm{X})\|_{F}\leq\|\bm{X}\|_{F}$. Under our assumptions about the size of $\delta_{k}$, we have $2\delta_{k}-\delta_{k}^{2}\geq\beta$ for all $k\geq 1$ and some $\beta>0$ and thus

The sequence $\{\|\mathcal{P}_{\Omega}(\bm{Y}^{k}-\bm{Y}^{\star})\|_{F}\}$ is nonincreasing and, therefore, converges to a limit.

As a consequence, $\|\bm{X}^{k}-\bm{X}^{\star}\|_{F}^{2}\rightarrow 0$ as $k\rightarrow\infty$.

2 General convergence theorem

Our second result is more general and establishes the convergence of the SVT iterations to the solution of (3.4) under general convex constraints. From now now, we will only assume that the function $\mathcal{F}(\bm{X})$ is Lipschitz in the sense that

We will assume to simplify that strong duality holds which is automatically true if the constraints obey constraint qualifications such as Slater’s condition .

We first establish the following preparatory lemma.

Let $(\bm{X}^{\star},\bm{y}^{\star})$ be a primal-dual optimal pair for (3.4). Then for each $\delta>0$, $\bm{y}^{\star}$ obeys

Proof. Recall that the projection $\bm{x}_{0}$ of a point $\bm{x}$ onto a convex set $\mathcal{C}$ is characterized by

Now because $\bm{y}^{\star}$ is dual optimal we have

Substituting the expression for the Lagrangian, this is equivalent to

We are now in the position to state our general convergence result.

Suppose that the sequence of step sizes obeys $0<\inf\delta_{k}\leq\sup\delta_{k}<2/\|L(\mathcal{F})\|^{2}$, where $L(\mathcal{F})$ is the Lipschitz constant in (4.6). Then assuming strong duality, the sequence $\{\bm{X}^{k}\}$ obtained via (3.5) converges to the unique solution of (3.4).

Proof. Let $(\bm{X}^{\star},\bm{y}^{\star})$ be primal-dual optimal for the problem (3.4). We claim that the optimality conditions give that for all $\bm{X}$

for some $\bm{Z}^{k}\in\partial f_{\tau}(\bm{X}^{k})$ and some $\bm{Z}^{\star}\in\partial f_{\tau}(\bm{X}^{\star})$. We justify this assertion by proving one of the two inequalities since the other is exactly similar. For the first, $\bm{X}^{k}$ minimizes $\mathcal{L}(\bm{X},\bm{y}^{k-1})$ over all $\bm{X}$ and, therefore, there exist $\bm{Z}^{k}\in\partial f_{\tau}(\bm{X}^{k})$ and $\bm{Z}_{i}^{k}\in\partial f_{i}(\bm{X}^{k})$, $1\leq i\leq m$, such that

Now write the first inequality in (4.8) for $\bm{X}^{\star}$, the second for $\bm{X}^{k}$ and sum the two inequalities. This gives

The rest of the proof is essentially the same as that of Theorem 4.5. It follows from Lemma 4.1 that

We continue and observe that because $\bm{y}^{\star}=[\bm{y}^{\star}+\delta_{k}\mathcal{F}(\bm{X})]_{+}$ by Lemma 4.3, we have

where we have put $L$ instead of $L(\mathcal{F})$ for short. Under our assumptions about the size of $\delta_{k}$, we have $2\delta_{k}-\delta_{k}^{2}L^{2}\geq\beta$ for all $k\geq 1$ and some $\beta>0$. Then

The problem (3.1) with linear constraints can be reduced to (3.4) by choosing

Suppose that the sequence of step sizes obeys $0<\inf\delta_{k}\leq\sup\delta_{k}<2/\|\mathcal{A}\|_{2}^{2}$. Then the sequence $\{\bm{X}^{k}\}$ obtained via (3.3) converges to the unique solution of (3.1).

Let $\|\mathcal{A}\|_{2}:=\sup\{\|\mathcal{A}(\bm{X})\|_{F}:\|\bm{X}\|_{F}=1\}$. With $\mathcal{F}(\bm{X})$ given as above, we have $|L(\mathcal{F})|^{2}=2\|\mathcal{A}\|^{2}_{2}$ and thus, Theorem 4.4 guarantees convergence as long as $0<\inf\delta_{k}\leq\sup\delta_{k}<1/\|\mathcal{A}\|_{2}^{2}$. However, an argument identical to the proof of Theorem 4.2 would remove the extra factor of two. We omit the details.

Implementation and Numerical Results

This section provides implementation details of the SVT algorithm—as to make it practically effective for matrix completion—such as the numerical evaluation of the singular value thresholding operator, the selection of the step size $\delta_{k}$, the selection of a stopping criterion, and so on. This section also introduces several numerical simulation results which demonstrate the performance and effectiveness of the SVT algorithm. We show that $30,000\times 30,000$ matrices of rank 10 are recovered from just about 0.4% of their sampled entries in a matter of a few minutes on a modest desktop computer with a 1.86 GHz CPU (dual core with Matlab’s multithreading option enabled) and 3 GB of memory.

To apply the singular value tresholding operator at level $\tau$ to an input matrix, it suffices to know those singular values and corresponding singular vectors above the threshold $\tau$. In the matrix completion problem, the singular value thresholding operator is applied to sparse matrices $\{\bm{Y}^{k}\}$ since the number of sampled entries is typically much lower than the number of entries in the unknown matrix $\bm{M}$, and we are hence interested in numerical methods for computing the dominant singular values and singular vectors of large sparse matrices. The development of such methods is a relatively mature area in scientific computing and numerical linear algebra in particular. In fact, many high-quality packages are readily available. Our implementation uses PROPACK, see for documentation and availability. One reason for this choice is convenience: PROPACK comes in a Matlab and a Fortran version, and we find it convenient to use the well-documented Matlab version. More importantly, PROPACK uses the iterative Lanczos algorithm to compute the singular values and singular vectors directly, by using the Lanczos bidiagonalization algorithm with partial reorthogonalization. In particular, PROPACK does not compute the eigenvalues and eigenvectors of $(\bm{Y}^{k})^{*}\bm{Y}^{k}$ and $\bm{Y}^{k}(\bm{Y}^{k})^{*}$, or of an augmented matrix as in the Matlab built-in function ‘svds’ for example. Consequently, PROPACK is an efficient—both in terms of number of flops and storage requirement—and stable package for computing the dominant singular values and singular vectors of a large sparse matrix. For information, the available documentation reports a speedup factor of about ten over Matlab’s ‘svds’. Furthermore, the Fortran version of PROPACK is about 3–4 times faster than the Matlab version. Despite this significant speedup, we have only used the Matlab version but since the singular value shrinkage operator is by-and-large the dominant cost in the SVT algorithm, we expect that a Fortran implementation would run about 3 to 4 times faster.

1.2 Step sizes

There is a large literature on ways of selecting a step size but for simplicity, we shall use step sizes that are independent of the iteration count; that is $\delta_{k}=\delta$ for $k=1,2,\ldots$. From Theorem 4.2, convergence for the completion problem is guaranteed (2.7) provided that $0<\delta<2$. This choice is, however, too conservative and the convergence is typically slow. In our experiments, we use instead

i.e. $1.2$ times the undersampling ratio. We give a heuristic justification below.

provided that the rank of $\bm{A}$ is not too large. The probability model is that $\Omega$ is a set of sampled entries of cardinality $m$ sampled uniformly at random so that all the choices are equally likely. In (5.2), we want to think of $\epsilon$ as a small constant, e.g. smaller than 1/2. In other words, the ‘energy’ of $\bm{A}$ on $\Omega$ (the set of sampled entries) is just about proportional to the size of $\Omega$. The near isometry (5.2) is a consequence of Theorem 4.1 in , and we omit the details.

Now returning to the proof of Theorem 4.2, we see that a sufficient condition for the convergence of (2.7) is

The reason why this is not a rigorous argument is that (5.2) cannot be applied to $\bm{A}=\bm{X}^{\star}-\bm{X}^{k}$ even though this matrix difference may obey the incoherence assumption. The issue here is that $\bm{X}^{\star}-\bm{X}^{k}$ is not a fixed matrix, but rather depends on $\Omega$ since the iterates $\{\bm{X}^{k}\}$ are computed with the knowledge of the sampled set.

1.3 Initial steps

The SVT algorithm starts with $\bm{Y}^{0}=\bm{0}$, and we want to choose a large $\tau$ to make sure that the solution of (2.8) is close enough to a solution of (1.1). Define $k_{0}$ as that integer obeying

Since $\bm{Y}^{0}=\bm{0}$, it is not difficult to see that

To save work, we may simply skip the computations of $\bm{X}^{1},\ldots,\bm{X}^{k_{0}}$, and start the iteration by computing $\bm{X}^{k_{0}+1}$ from $\bm{Y}^{k_{0}}$.

This strategy is a special case of a kicking device introduced in ; the main idea of such a kicking scheme is that one can ‘jump over’ a few steps whenever possible. Just like in the aforementioned reference, we can develop similar kicking strategies here as well. Because in our numerical experiments the kicking is rarely triggered, we forgo the description of such strategies.

1.4 Stopping criteria

Here, we discuss stopping criteria for the sequence of SVT iterations (2.7), and present two possibilities.

The first is motivated by the first-order optimality conditions or KKT conditions tailored to the minimization problem (2.8). By (2.14) and letting $\partial_{\bm{Y}}g_{0}(\bm{Y})=\bm{0}$ in (2.13), we see that the solution $\bm{X}^{\star}_{\tau}$ to (2.8) must also verify

where $\bm{Y}$ is a matrix vanishing outside of $\Omega^{c}$. Therefore, to make sure that $\bm{X}^{k}$ is close to $\bm{X}^{\star}_{\tau}$, it is sufficient to check how close $(\bm{X}^{k},\bm{Y}^{k-1})$ is to obeying (5.4). By definition, the first equation in (5.4) is always true. Therefore, it is natural to stop (2.7) when the error in the second equation is below a specified tolerance. We suggest stopping the algorithm when

where $\epsilon$ is a fixed tolerance, e.g. $10^{-4}$. We provide a short heuristic argument justifying this choice below.

In the matrix completion problem, we know that under suitable assumptions

which is just (5.2) applied to the fixed matrix $\bm{M}$ (the symbol $\asymp$ here means that there is a constant $\epsilon$ as in (5.2)). Suppose we could also apply (5.2) to the matrix $\bm{X}^{k}-\bm{M}$ (which we rigorously cannot since $\bm{X}^{k}$ depends on $\Omega$), then we would have

In words, one would control the relative reconstruction error by controlling the relative error on the set of sampled locations.

A second stopping criterion comes from duality theory. Firstly, the iterates $\bm{X}^{k}$ are generally not feasible for (2.8) although they become asymptotically feasible. One can construct a feasible point from $\bm{X}^{k}$ by projecting it onto the affine space $\{\bm{X}:\mathcal{P}_{\Omega}(\bm{X})=\mathcal{P}_{\Omega}(\bm{M})\}$ as follows:

Secondly, using the notations of Section 2.4, duality theory gives that

Therefore, $b_{k}-a_{k}$ is an upper bound on the duality gap and one can stop the algorithm when this quantity falls below a given tolerance.

1.5 Algorithm

We conclude this section by summarizing the implementation details and give the SVT algorithm for matrix completion below (Algorithm 1). Of course, one would obtain a very similar structure for the more general problems of the form (3.1) and (3.4) with linear inequality constraints. For convenience, define for each nonnegative integer $s\leq\min\{n_{1},n_{2}\}$,

where $\bm{U}^{k}=[\bm{u}_{1}^{k},\ldots,\bm{u}_{s}^{k}]$ and $\bm{V}^{k}=[\bm{v}_{1}^{k},\ldots,\bm{v}_{s}^{k}]$ are the first $s$ singular vectors of the matrix $\bm{Y}^{k}$, and $\bm{\Sigma}^{k}$ is a diagonal matrix with the first $s$ singular values $\sigma^{k}_{1},\ldots,\sigma_{s}^{k}$ on the diagonal.

2 Numerical results

Our implementation is in Matlab and all the computational results we are about to report were obtained on a desktop computer with a 1.86 GHz CPU (dual core with Matlab’s multithreading option enabled) and 3 GB of memory. In our simulations, we generate $n\times n$ matrices of rank $r$ by sampling two $n\times r$ factors $\bm{M}_{L}$ and $\bm{M}_{R}$ independently, each having i.i.d. Gaussian entries, and setting $\bm{M}=\bm{M}_{L}\bm{M}_{R}^{*}$ as it is suggested in . The set of observed entries $\Omega$ is sampled uniformly at random among all sets of cardinality $m$.

The recovery is performed via the SVT algorithm (Algorithm 1), and we use

Our computational results are displayed in Table 1. There, we report the run time in seconds, the number of iterations it takes to reach convergence (5.7), and the relative error of the reconstruction

where $\bm{M}$ is the real unknown matrix. All of these quantities are averaged over five runs. The table also gives the percentage of entries that are observed, namely, $m/n^{2}$ together with a quantity that we may want to think as the information oversampling ratio. Recall that an $n\times n$ matrix of rank $r$ depends upon $d_{r}:=r(2n-r)$ degrees of freedom. Then $m/d_{r}$ is the ratio between the number of sampled entries and the ‘true dimensionality’ of an $n\times n$ matrix of rank $r$.

The first observation is that the SVT algorithm performs extremely well in these experiments. In all of our experiments, it takes fewer than 200 SVT iterations to reach convergence. As a consequence, the run times are short. As indicated in the table, we note that one recovers a $1,000\times 1,000$ matrix of rank $10$ in less than a minute. The algorithm also recovers $30,000\times 30,000$ matrices of rank $10$ from about $0.4\%$ of their sampled entries in just about 17 minutes. In addition, higher-rank matrices are also efficiently completed: for example, it takes between one and two hours to recover $10,000\times 10,000$ matrices of rank $100$ and $20,000\times 20,000$ matrices of rank $50$. We would like to stress that these numbers were obtained on a modest CPU (1.86GHz). Furthermore, a Fortran implementation is likely to cut down on these numbers by a multiplicative factor typically between three and four.

We emphasized all along an important feature of the SVT algorithm, which is that the matrices $\bm{X}^{k}$ have low rank. We demonstrate this fact empirically in Figure 1, which plots the rank of $\bm{X}^{k}$ versus the iteration count $k$, and does this for unknown matrices of size $5,000\times 5,000$ with different ranks. The plots reveal an interesting phenomenon: in our experiments, the rank of $\bm{X}^{k}$ is nondecreasing so that the maximum rank is reached in the final steps of the algorithm. In fact, the rank of the iterates quickly reaches the value $r$ of the true rank. After these few initial steps, the SVT iterations search for that matrix with rank $r$ minimizing the objective functional. As mentioned earlier, the low-rank property is crucial for making the algorithm run fast.

Finally, we demonstrate the results of the SVT algorithm for matrix completion from noisy sampled entries. Suppose we observe data from the model

where $\bm{Z}$ is a zero-mean Gaussian white noise with standard deviation $\sigma$. We run the SVT algorithm but stop early, as soon as $\bm{X}^{k}$ is consistent with the data and obeys

where $\epsilon$ is a small parameter. Our reconstruction $\hat{\bm{M}}$ is the first $\bm{X}^{k}$ obeying (5.10). The results are shown in Table 2 (the quantities are averages of 5 runs). Define the noise ratio as

and the relative error by (5.8). From Table 2, we see that the SVT algorithm works well as the relative error between the recovered and the true data matrix is just about equal to the noise ratio.

The theory of low-rank matrix recovery from noisy data is nonexistent at the moment, and is obviously beyond the scope of this paper. Having said this, we would like to conclude this section with an intuitive and nonrigorous discussion, which may explain why the observed recovery error is within the noise level. Suppose again that $\hat{\bm{M}}$ obeys (5.6), namely,

As mentioned earlier, one condition for this to happen is that $\bm{M}$ and $\hat{\bm{M}}$ have low rank. This is the reason why it is important to stop the algorithm early as we hope to obtain a solution which is both consistent with the data and has low rank (the limit of the SVT iterations, $\lim_{k\rightarrow\infty}\bm{X}^{k}$, will not generally have low rank since there may be no low-rank matrix matching the noisy data). From

and the fact that both terms on the right-hand side are on the order of $\sqrt{m\sigma^{2}}$, we would have $p\|\hat{\bm{M}}-\bm{M}\|_{F}^{2}=O(m\sigma^{2})$ by (5.11). In particular, this would give that the relative reconstruction error is on the order of the noise ratio since $\|\mathcal{P}_{\Omega}(\bm{M})\|_{F}^{2}\asymp p\|\bm{M}\|_{F}^{2}$—as observed experimentally.

2.2 Linear inequality constraints

We now examine the speed at which one can solve similar problems with linear inequality constraints instead of linear equality constraints. We assume the model (5.9), where the matrix $\bm{M}$ of rank $r$ is sampled as before, and solve the problem (3.8) by using (3.10). We formulate the inequality constraints in (3.8) with $E_{ij}=\sigma$ so that one searches for a solution $\hat{\bm{M}}$ with minimum nuclear norm among all those matrices whose sampled entries deviate from the observed ones by at most the noise level $\sigma$.This may not be conservative enough from a statistical viewpoint but this works well in this case, and our emphasis here is on computational rather than statistical issues. In this experiment, we adjust $\sigma$ to be one tenth of a typical absolute entry of $\bm{M}$, i.e. $\sigma=0.1\,\sum_{ij\in\Omega}|M_{ij}|/m$, and the noise ratio as defined earlier is 0.780. We set $n=1,000$, $r=10$, and the number $m$ of sampled entries is five times the number of degrees of freedom, i.e. $m=5d_{r}$. Just as before, we set $\tau=5n$, and choose a constant step size $\delta=1.2p^{-1}$.

The results, reported in Figure 2, show that the algorithm behaves just as well with linear inequality constraints. To make this point, we compare our results with those obtained from noiseless data (same unknown matrix and sampled locations). In the noiseless case, it takes about 150 iterations to reach the tolerance $\epsilon=10^{-4}$ whereas in the noisy case, convergence occurs in about 200 iterations (Figure 2(a)). In addition, just as in the noiseless problem, the rank of the iterates is nondecreasing and quickly reaches the true value $r$ of the rank of the unknown matrix $\bm{M}$ we wish to recover (Figure 2(b)). As a consequence the SVT iterations take about the same amount of time as in the noiseless case (Figure 2(c)) so that the total running time of the algorithm does not appear to be substantially different from that in the noiseless case.

We close by pointing out that from a statistical point of view, the recovery of the matrix $\bm{M}$ from undersampled and noisy entries by the matrix equivalent of the Dantzig selector appears to be accurate since the relative error obeys $\|\hat{\bm{M}}-\bm{M}\|_{F}/\|\bm{M}\|_{F}=0.0769$ (recall that the noise ratio is about $0.08$).

Discussion

This paper introduced a novel algorithm, namely, the singular value thresholding algorithm for matrix completion and related nuclear norm minimization problems. This algorithm is easy to implement and surprisingly effective both in terms of computational cost and storage requirement when the minimum nuclear-norm solution is also the lowest-rank solution. We would like to close this paper by discussing a few open problems and research directions related to this work.

Our algorithm exploits the fact that the sequence of iterates $\{\bm{X}^{k}\}$ have low rank when the minimum nuclear solution has low rank. An interesting question is whether one can prove (or disprove) that in a majority of the cases, this is indeed the case.

It would be interesting to explore other ways of computing $\mathcal{D}_{\tau}(\bm{Y})$—in words, the action of the singular value shrinkage operator. Our approach uses the Lanczos bidiagonalization algorithm with partial reorthogonalization which takes advantages of sparse inputs but other approaches are possible. We mention two of them.

A series of papers have proposed the use of randomized procedures for the approximation of a matrix $\bm{Y}$ with a matrix $\bm{Z}$ of rank $r$ . When this approximation consists of the truncated SVD retaining the part of the expansion corresponding to singular values greater than $\tau$, this can be used to evaluate $\mathcal{D}_{\tau}(\bm{Y})$. Some of these algorithms are efficient when the input $\bm{Y}$ is sparse , and it would be interesting to know whether these methods are fast and accurate enough to be used in the SVT iteration (2.7).

A wide range of iterative methods for computing matrix functions of the general form $f(\bm{Y})$ are available today, see for a survey. A valuable research direction is to investigate whether some of these iterative methods, or other to be developed, would provide powerful ways for computing $\mathcal{D}_{\tau}(\bm{Y})$.

In practice, one would like to solve (2.8) for large values of $\tau$. However, a larger value of $\tau$ generally means a slower rate of convergence. A good strategy might be to start with a value of $\tau$, which is large enough so that (2.8) admits a low-rank solution, and at the same time for which the algorithm converges rapidly. One could then use a continuation method as in to increase the value of $\tau$ sequentially according to a schedule $\tau_{0},\tau_{1},\ldots$, and use the solution to the previous problem with $\tau=\tau_{i-1}$ as an initial guess for the solution to the current problem with $\tau=\tau_{i}$ (warm starting). We hope to report on this in a separate paper.

J-F. C. is supported by the Wavelets and Information Processing Programme under a grant from DSTA, Singapore. E. C. is partially supported by the Waterman Award from the National Science Foundation and by an ONR grant N00014-08-1-0749. Z. S. is supported in part by Grant R-146-000-113-112 from the National University of Singapore. E. C. would like to thank Benjamin Recht and Joel Tropp for fruitful conversations related to this project, and Stephen Becker for his help in preparing the computational results of Section 5.2.2.

References