Simple Atom Selection Strategy for Greedy Matrix Completion

1 Introduction

Low rank matrix completion is among the most basic problems in machine learning and data analysis. It plays a key role in solving many important problems, such as collaborative filtering [Rennie and Srebro, 2005; Koren et al., 2009; Rendle et al., 2009], dimensionality reduction [Weinberger and Saul, 2006; So and Ye, 2005], clustering [Eriksson et al., 2011; Yi et al., 2012], and multi-class learning [Argyriou et al., 2008; Obozinski et al., 2010; Xu et al., 2013].

Matrix completion can be formulated as seeking the matrix with lowest rank that fits the observed data. However, directly solving such problem is NP-hard and of little practical use [Chistov and Grigor'ev, 1984] which leads to many approximation strategies. One principled approach is to adopt nuclear norm as surrogate for the rank [Cai et al., 2010; Jain et al., 2010; Lin et al., 2010; Mazumder et al., 2010; Toh and Yun, 2010]. Non-convex surrogates, usually complex pseudo norms, have also been brought forth to gain better accuracy or nearly unbiased estimation [Liu et al., 2013]. Although recovery guarantees in these contexts are established [Cand`es and Recht, 2009; Cand`es and Tao, 2010; Keshavan et al., 2010], the demand of expensive truncated Singular Value Decomposition (SVD) prevents the applications of nuclear norm based methods to large real-world problems.

∗Corresponding author

Recently, remarkable progress has been made for greedy matrix completion techniques. The idea behind them is to represent matrix A ∈ Rm×n as a sparse code over a dictionary of infinite unit rank-one matrices which are referred to as atoms [Lee and Bresler, 2010]. Such representation makes low rank matrix completion a natural extension of greedy selection for optimization with sparsity constraint [Mallat and Zhang, 1993; Pati et al., 1993; Tropp, 2004; Shalev-Shwartz et al., 2010; Zhang, 2011] to the matrix case.

Typically, greedy matrix completion algorithms, like GECO [Shalev-Shwartz et al., 2011], R1MP and ER1MP [Wang et al., 2014], BOOST [Zhang et al., 2012], and JS [Jaggi et al., 2010], proceed in two core steps in each iteration. The first step selects a locally optimal atom. The second step refines the weights of all atoms chosen up to this iteration. Since atom selection and weight refinement can be much cheaper than truncated SVD, such a two-step scheme brings us better scalability than nuclear norm based methods.

In the previous works, the second step, i.e., the weight refinement step was the research focus and almost all the existing greedy matrix algorithms differ mainly in their refinement steps[Wang et al., 2014]. For the first step, i.e., the atom selection step, only one strategy, called T1SVD in our paper, was used in current greedy matrix completion literatures, to our best knowledge. The main reason is that the T1SVD strategy corresponds to a Top-1 SVD problem which is numerically easy to solve and has plenty of efficient algorithms.

In this paper, we further explore the atom selection problem and present a simple strategy, called Optimal Atom (OA), to select the best atom. Our research is partially inspired by the work of optimization problem with sparsity constraint [Shalev-Shwartz et al., 2010; Liu et al., 2014]. We directly solve a coordinate optimization problem to find the optimal atom in each step, rather than deal with the first order approximation to the optimal atom choice as T1SVD does. In this line, our Optimal Atom based Matrix Completion algorithm (OAMC) adopts an alternating method instead of the common Top-1 SVD solver to conduct the matrix completion. We show that OA is a better strategy to construct the greedy matrix completion algorithm than T1SVD. The major contributions are summarized as follows:

• We propose a simple atom selection strategy, called OA, which finds the optimal atom in each iteration by alternating minimization with computational complexity comparable to common Top-1 SVD solvers.
• Under suitable assumptions, we construct an upper bound of the approximating error for OA strategy, which is independent of the largest singular value of underlying residual matrix. Such result is important since it provides a tighter training error bound for OAMC than for the T1SVD baselines.
• A greedy algorithm, called OAMC, is devised to solve the matrix completion problem. Several variants are designed according to different weight refinement methods. We prove that OAMC and three of its variants have the property of linear convergence.

In the experiments, we evaluate the performance of the proposed OAMC algorithm by applying it to the tasks of Recommendation and Image Recovery. We perform matrix completion on three largest publicly available recommendation datasets: MovieLens 10M, NetFlix, and Yahoo Music. In all these experiments, OAMC and its variants significantly outperform their T1SVD-based competitors in terms of both speed and accuracy. Some of these variants are 5 times faster than existing methods when using efficient random initialization for the alternating minimization scheme. Besides, we are able to process Yahoo Music dataset in about 700 seconds in PC workstation and achieve a root mean square error 24.5 on the test set.

2 Preliminaries

For a matrixM∈ Rm×n , let Ω ⊂ {1, · · · ,m}×{1, · · · , n} denote the indices of observed entries. In this paper, we always assume m ≤ n. We consider the following low rank matrix completion problem:

where K ≪ min(m, n) is a constant, ‖ · ‖F is the Frobenius norm, and the operator PΩ : Rm×n ↦→ Rm×n is defined as follows:

Usually, we call the Ω as the support of PΩ(A). In practice, problem (1) is relaxed by replacing the rank(X) with a surrogate function to fit more effective algorithms.

Actually, we can also depict problem (1) in an infinite vector form 1 . Consider a dictionary like D = U × V, in which U = {u ∈ Rm : ‖u‖ = 1}, V = {v ∈ Rn : ‖v‖ = 1}, where we let ‖ · ‖ be the 2-norm of vector. A pair (u, v) ∈ D denotes an atom in this dictionary. Over D, a vector λ ∈ R|D|

1 We follow the notation of [Shalev-Shwartz et al., 2011], in which an assumption of finite representation for real number is used to simplify the presentation.

can be used to represent an arbitrary m× n matrix. That is, given vector λ ∈ R|D| , we have a correspondent m×n matrix

where X : R|D| ↦→ Rm×n is a linear map, and λ(u,v) ∈ R denotes the value of λ in the coordinate indexed by the pair (u, v). From the SVD theorem, if rank(X(λ)) ≤ r, then there must be a λ satisfying ‖λ‖0 ≤ r. We also define a standard basis vector e(u,v) ∈ R|D| over this dictionary as

where p ∈ U, and q ∈ V. The difference between e(u,v) and standard basis vector ei in Euclidean space is that e(u,v) is indexed by pair (u, v) instead of the number i.

Thus, greedy algorithm for (1) can be developed by resorting to the following equivalent problem:

where Q(λ) ≜ L(X(λ)) = 1 2 ‖PΩ(M - X(λ))‖2F. For convenience, we also define residual function R : R|D| ↦→ Rm×n as R(λ) = PΩ(M-X(λ)). We call M-X(λ) the underlying matrix of R(λ). To solve problem (2), the current state-of-the-art greedy methods choose to find the atom using maximum gradient in each iteration. That is, the selected atom comes from following formulation

We call this atom selection method maximum gradient strategy or T1SVD since it is based on solving Top-1 SVD problem.

A few more notations will be useful in our narration. Given a rank r matrix A with SVD ∑ri=1 σiuivi⊤, σ1 ≥ · · · ≥ σr. We define S1(A) = σ1u1v⊤ 1 , S2(A) = ∑ri=2 σiuiv⊤ i , and σ1(A) = σ1(the maximum σi, for i = 1, · · · , r). We use Ai to represent the ith column of A and Ai,: to represent the ith row. O denotes the big-O notation in mathematics. 〈·, ·〉 represents the inner product of two matrices, and ∇ denotes the Del operator of a function.

3 Methodology

We start from λ(0) = 0. In the k-th iteration , suppose λ = λ(k) . The locally optimal u∗ , v∗ , and α∗ need to be estimated to approximate the residual R(λ). First of all, for a fixed pair (u, v), we should find an α that minimizes Q(λ + αe(u,v) ) which results in the optimization problem minαQ(λ + αe(u,v) ). Second, we expect that (u, v) can make the maximum progress after an increment αe(u,v) is added into λ. Thus, by combining these two goals, we have the following optimization problem for each iteration.

Algorithm 1 OAMC

Input: Ω,PΩ(M),K Output: X(K) 1: Initialize: Set X(0) = 0 and R(0) = PΩ(M) 2: for k = 1, 2, . . . ,K do 3: ¯u := INIT() 4: repeat 5: ¯v := argminv ‖PΩ(R(k-1) - ¯uv⊤)‖F 6: ¯u := argminu ‖PΩ(R(k-1) - u¯v⊤)‖F 7: until converge 8: X(k) :=X(k) + ¯u¯v⊤ 9: R(k) := PΩ(M-X(k) ) 10: end for

3.1 Optimal Atom

Basically, (4) indicates that we should find an optimal coordinate e(u∗ ,v∗ ) and a proper step α∗ . From the definition of Q(λ), we have the following objective

Problem (5) has been widely thought to be complicated to solve because the objective is jointly non-convex over u and v. However, it's easy to know that αuv⊤ is a rank 1 matrix. To such a rank invariant situation, if we fixed u, the original problem becomes convex, specifically a typical least square problem. The same can be found when v is fixed. Thus, we can alternately fix one and optimize over the other until convergence. This alternating optimization scheme is a special case of the Alternating Least Square Method [Koren et al., 2009; Jain et al., 2013; Gunasekar et al., 2013].

Practically, proper initialization of u is crucial for convergence for an alternating optimization procedure. In our method, random and prior-knowledge based initialization methods, are tested (see [Gunasekar et al., 2013] for details), and encapsulated in a macro INIT. In addition, α is redundant. Real implementation can use only two variables: ¯u ∈ Rm and ¯v ∈ Rn since it can be easily derived that u = ¯u/‖¯u‖,v = ¯v/‖¯v‖ and α = ‖¯u‖‖¯v‖. Furthermore, we use m × n matrix X and R to represent the function X(λ) and R(λ) since λ with infinite dimension is simply used in derivation and analysis. We summarize our pseudo code in Algorithm 1.

Note that, the first important difference between OA and T1SVD is that the selected atom in the latter is not derived from (4). Actually T1SVD simplifies (4) by replacing Q(λ+ αe(u,v) ) with

which is a first-order approximation of Q(λ+αe(u,v) ). Since the first item of (6) is irrelevant to (u, v), by restricting α to be finite, we have

instead of problem (4). Apparently, α will always take the reverse sign of 〈∇Q(λ), e(ˆu,ˆv) 〉, thus we have

which is the problem (3). It is reasonable to infer that the error from the linear approximation may amplify the necessity of a fully corrective procedure in current state-of-the-art methods.

The second difference between OA and T1SVD lies in that the optimization problem of OA involves only the matrix entries in Ω while T1SVD deals with the whole matrix with plenty of zeros indicating the missing entries. For T1SVD, the chain rule of partial derivative allows us to write (3) as

in which equation (7) comes from the fact that

And for equation (8), ∇L(X(λ)) = -R(λ) is easy to verify. Comparing (7) with (5), T1SVD chooses the best atom (ˆu, ˆv) to approximate the whole R(λ), while OA goes through the same procedure only within Ω which contains the indices of observed entries.

3.2 Variants for Fully Corrective Selection

Our Algorithm (1) is totally non-corrective, that is, at each iteration, we only modify the weights of the current atom. Based on OA, variants for fully corrective selection can also be developed for better accuracy, especially when sufficient computational capacity is available.

In order to compare our OA strategy with other state-ofthe-art algorithms, we design OA variants with five mainstream schemes for weight updating and encapsulate them in the form of macros to keep the algorithm succinct. The inputs of these macros are {Ω,PΩ(M),U ∈ Rm×q ,V ∈ Rm×q , Φ ∈ R(q-1)×(q-1)}, where U and V have their column Ui ∈ U and Vi ∈ V respectively, Φ ∈ R(q-1)×(q-1) is obtained from the previous iteration and q is the iteration count. The outputs are {Φ ∈ Rq×q , U, V}. U, V, and Φ are in both input and output set. Our variants for fully corrective selection are summarized in Algorithm 2. We briefly explain the weight updating schemes as follows.

ADJ-GECO follows the method in literature [ShalevShwartz et al., 2011] that solves the following regression problem

Algorithm 2 OA variants: OA-GECO,OA-R1MP, OAER1MP, OA-JS, and OA-BOOST

Let USΦSV⊤ S be the SVD of S∗ . We construct outputs as U := UUS, V := VVS and Φ := ΦS.

ADJ-R1MP is a simplification of ADJ-GECO with a constraint that S is diagonal [Wang et al., 2014].

ADJ-ER1MP further simplifies above ADJ-R1MP by setting Si,i = a1Φi,i for i ∈ {1, · · · , (q - 1)} and Sq,q = a2 [Wang et al., 2014].

ADJ-JS preprocesses the data as in [Jaggi et al., 2010] and proceeds like ADJ-ER1MP except for an additional constraint (a1 + a2 = 1) and a1 and a2 are the only variables.

ADJ-BOOST is similar to ADJ-ER1MP. It adds a regularization term β(a1 ∑q-1 i=1 Φi,i+a2) into the objective function with constraints: a1 ≥ 0 and a2 ≥ 0 [Zhang et al., 2012].

4 Analysis

In this section, we investigate how well OA approximates the residual R. The upper bound that we obtain is independent of the largest singular value of the underlying matrix of residual. We also prove that OAMC and three of its variants converge linearly.

We start our analysis with the following definition of incoherence property of matrix.

Definition 1 (µ-incoherent). Let A ∈ Rm×n be a rank k matrix with singular value decomposition UΦV⊤. A is said to be µ-incoherent if there exists a constant µ such that

Lemma 1. Let M = σuv⊤ be a µ-incoherent matrix with u ∈ U and v ∈ V and N ∈ Rm×n be a noise matrix with maxi,j Ni,j ≤ cσ n , where c is a constant. Suppose that the support Ω of R is obtained by uniformly and independently sampling from {1, · · · ,m} × {1, · · · , n} with probability p. Let G = (‖PΩ(N)‖F/√p)/σ. If G ≤ C1 ≤ √3/µ and

where C1, C2 and δ2 ≤ (1/64) are constants, then

with probability at least 1 - (1/n3 ), where (¯u, ¯v) is the solution of step 1 of Algorithm 2 initialized as in [Gunasekar et al., 2013] with O(log(1/µG)) iterations, and C is a constant.

Lemma 1 will be used to show that our algorithm needs approximately pmn ≥ C2µ4 n logn log 1 µG min{δ2,n(µG)4 } samples to fulfill the completion, whose proof is placed in the long version of this paper. Note that O(n logn) is the optimum sampling complexity to complete a rank-1 matrix according to the Coupon collector's problem. Thus our atom selection strategy is also optimum in terms of sampling complexity with respect to the matrix size.

Theorem 1. Let L ∈ Rm×n be the underlying rank-r matrix of the residual R and suppose maxi,j S2(L) ≤ cσ1(L) n , with a constant c. If the support Ω of R is obtained by uniformly and independently sampling from {1, · · · ,m} × {1, · · · , n} with probability p defined in Lemma 1 using µ and G derived from S1(L) and S2(L), then

with probability at least 1- 1 n3 , where ¯u¯v⊤ is the atom constructed by step 1 of Algorithm 2 and C is a positive constant.

Proof. By the subadditivity of Frobenius norm, we have:

where u1 and v1 are vectors from S1(L). Taking α = σ1(L),

We bound the second term of (11) by Cµ‖PΩ(L≥2)‖F using Lemma 1. Since PΩ(R) = PΩ(L), we have the result.

This bound shows that the result of Algorithm 2 approximates R with an error independent of the largest singular value of L.

The following Lemma is used to prove the linear convergence of OAMC and three of its variants.

Lemma 2. Let (¯u, ¯v) be the atom selected by the step 1 of Algorithm 2 and the step 1 is initialized by any pair (ˆu, ˆv). We have

The result of Lemma 2 can be easily obtained since the alternating minimization does not increase the value of the objective in any iteration. Now we present a simple proof to show that OAMC and three of its variants converge linearly with proper initialization. We choose OA-R1MP as an example and use notations conforming to the description of Algorithm 2.

Theorem 2. The residual R(k+1) ∈ Rm×n of OA-R1MP satisfies

with the ApproxSV defined in [Shalev-Shwartz et al., 2011] as the initialization for step 1 in each iteration, where γ ∈ [0, 1) is a constant.

Proof. From the definition of ADJ-R1MP, we have

Let Φ′ be in the form of [[(Φ(k) )⊤, 0]⊤, [0, α]⊤]. We have

Let ˆuˆv⊤ be the initial atom constructed by ApproxSV. By applying Lemma 2 to (14) we have

For (15), α has a close form solution

Substituting α∗ into (15), we have

We know ‖PΩ(ˆuˆv⊤)‖2F ≤ 1 and 〈PΩ(R(k) ), PΩ(ˆuˆv⊤)〉 = 〈PΩ(R(k) ), ˆuˆv⊤〉 = 〈R(k) , ˆuˆv⊤〉 = ˆu⊤R(k) ˆv ≥ (1 - δ)σ1(R(k) ), where δ is a constant smaller than 1. The last inequation comes from [Shalev-Shwartz et al., 2011]. Thus,

Further, we have

It is well known that σ2 1 (R(k) ) ≥ ‖R(k) ‖2F rank(R(k) ) ≥ ‖R(k) ‖2F min(m,n) , from which we have (1-δ)2 σ1(R(k) )2 ‖R(k) ‖2F < 1. Thus, there must be a constant γ < 1 that satisfies ‖R(k+1)‖F ≤ γ k ‖PΩ(M)‖F.

Note that OAMC, OA-GECO, and OA-ER1MP can also be found convergent linearly by slightly modifying (14) together with a proper initialization of the alternating minimization step.

5 Experiment

To make empirical evaluation, we conduct experiments of Recommendation and Image Recovery. OAMC and its variants are compared with five state-of-the-art T1SVD based competitors including GECO, R1MP, ER1MP, JS, and BOOST. All experiments are conducted on the same PC (Windows Server 2012 R2, Intel Xeon E5 2690v2*2 CPU, and 128G RAM).

We call the PROPACK to solve Top-1 SVD for T1SVD strategy. As for alternating minimization, we use random and the ApproxSV initialization respectively and set the maximum number of iterations to be ten. It turns out that random initialization does little harm to the convergence rate and the accuracy in our experiments. We only report the experiments with random initialization due to limited space.

For the parameter setting, we set the same maximum number of iteration for all the algorithms. And λ, the regularization parameter for BOOST, is selected by 3-fold cross validation. Additionally, JS requires a regularization parameter t, which is set to a doubled value of the nuclear norm solved by OA-ER1MP (This value is close in all OA variants).

To measure the performance, Root Mean Square Errors (RMSE) on both training set and testing set are calculated. We also record their running time (in seconds). Experimental results show that using OA strategy speeds up the convergence and reduces the test error. Additionally, noncorrective algorithm OAMC performs well along with its variants, which implies that OA strategy reduces the importance of weight refinement procedure.

5.1 Recommendation

We use three largest publicly available datasets: MovieLens10M, NetFlix, and Yahoo Music to test the matrix completion based recommendation. The statistics of these datasets are listed in Table 1. All datasets are randomly split into equal-sized training and testing parts. The maximum iteration K for them are {20, 20, 50} respectively. λ is chosen from {10i , i ∈ {-1, · · · , 3}} by 3-fold cross validation.

In Figure 1, we plot the logarithm of training error of Yahoo Music dataset in each iteration to compare the convergence rates of different algorithms. We can observe two phases in all OA based methods in the figure. The first few iterations drastically reduce the training error, which we attribute to the approximation guarantee of OA: Theorem 1 shows that, when the gap between the singular values of the underlying matrix is large, OA can remove the impact of the largest singular value on the training error. In the following iterations, the training error decreases slower, but still linearly. This can be explained by Theorem 2: convergence is linear regardless of the distribution of the singular values. Additionally, all OA based methods have smaller training error than T1SVD based methods, even just after one iteration. Further, OAMC, OA-GECO, OA-R1MP, and OA-ER1MP have similar performance. This suggests the inequality (14) is indeed tight, which we attribute to the construction of OA. As for T1SVD based methods, GECO outperforms the rest by much. Such contrast makes it reasonable to infer that OA strategy substantially diminishes the contribution of weight update procedure.

We then plot test error over the running time to show the accuracy and efficiency of our methods in Figure 2.

It shows that OA based methods achieve small test error with quite little known entries while T1SVD based methods fail in these situations. One reason is that the alternating minimization has recovery guarantee [Jain et al., 2013; Gunasekar et al., 2013] which ensures the better reconstruction of underlying matrix. We can also see that all OA based methods are more efficient than their T1SVD competitors and some OA variants are even 5 times faster. It is worth emphasizing that some of our algorithms take only 700 seconds to process Yahoo Music dataset in PC, yet achieve a root mean square error 24.5. We attribute this to the efficiency of alternating minimization procedure.

5.2 Image Recovery

In Image Recovery, we use five 512 × 512 sized gray-scale benchmark images 2 . Since images are typically high rank, we set rank K = 200. We uniformly retain 20% pixels as known entries. By conducting 10 independent trials for each image, average RMSE and running time are presented in Table 2.

We first compare the test error of two methods using the same weight refinement strategy but with different atom selection methods. The advantage of OA over T1SVD is clear, as all OA based methods outperform their T1SVD based competitor. Furthermore, we can see that even without further weight refinement, OAMC outperforms T1SVD based algorithms in most cases. Besides, OA based methods are also much faster, for example, OA-ER1MP is at least 7 times faster than ER1MP on every image. Explanation for such observation is similar to the one in Recommendation.

6 Conclusion

We propose a novel atom selection strategy for greedy matrix completion called Optimal Atom, based on which several algorithms are derived as well. Both approximation guarantee of OA and the convergence rate of our variants are established. Through two applications, Recommendation and Image Recovery, we demonstrate the superiority of our methods over existing T1SVD based algorithms. In the future work, we will further investigate the weight refinement step for OA.

2 http://www.utdallas.edu/ cxc123730/mh bcs spl.html

Acknowledgments

This work is partially supported by National Natural Science Foundation of China (Grant No: 61472347). Hui Qian is supported by China Scholar Council (Grant No: 201406325003).