Accelerated Inexact Soft-Impute for Fast Large-Scale Matrix Completion

1 Introduction

In many applications, the data are stored in a matrix. Given a partially observed matrix O, matrix factorization attempts to recover a low-rank matrix X that best approximates O on the observed entries [Cand`es and Recht, 2009]. Matrix factorization has been widely used in a variety of domains. Examples include collaborative filtering for recommender systems [Mazumder et al., 2010], link prediction for social networks [Kim and Leskovec, 2011], image inpainting [Liu et al., 2013], and multilabel learning [Cabral et al., 2011].

However, direct minimization of the matrix rank is NPhard, and the nuclear norm (which is the sum of singular values) is often used instead. Mathematically, let O ∈ Rm×n (where, without loss of generality, we assume that m ≥ n), and the positions of the observed entries be indicated by Ω ∈ {0, 1}m×n , such that Ωij = 1 if Oij is observed, and 0 otherwise. Matrix factorization is formulated as the following optimization problem

where [PΩ(A)]ij = Aij if Ωij = 1, and 0 otherwise; and ‖X‖∗ is the nuclear norm of X. It is known that the nuclear norm is the tightest convex lower bound of the rank [Recht et al., 2010]. Besides, though the nuclear norm is only a surrogate, there are theoretical guarantees that the underlying matrix can be recovered [Cand`es and Recht, 2009].

Computationally, though the nuclear norm is nonsmooth, problem (1) can be solved by various optimization tools. An early attempt is based on reformulating (1) as a semidefinite program (SDP) [Cand`es and Recht, 2009]. However, SDP solvers have large time and space complexities, and are only suitable for small data sets. For large-scale matrix completion, Cai et al. [2010] pioneered the use of firstorder methods and proposed the singular value thresholding (SVT) algorithm. However, a singular value decomposition (SVD) is required in each SVT iteration. This takes O(mn2 ) time and can be computationally expensive. In [Toh and Yun, 2010], this is reduced to a partial SVD by computing only the leading singular values/vectors using PROPACK (a variant of the Lanczos algorithm) [Larsen, 1998]. Another major breakthrough is made by the Soft-Impute algorithm [Mazumder et al., 2010], which utilizes a special "sparse plus low-rank" structure associated with the SVT to efficiently compute the SVD. Empirically, this allows Soft-Impute to perform matrix completion on the entire Netflix data set with a reasonable time.

The SVT algorithm can also be viewed as a proximal gradient algorithm [Tibshirani, 2010]. Hence, it converges with a O(1/T) rate, where T is the number of iterations [Beck and Teboulle, 2009; Nesterov, 2013]. Later, this is further "accelerated", and the convergence rate is improved to O(1/T 2 ) [Ji and Ye, 2009; Toh and Yun, 2010]. However, Tibshirani [2010] suggested that this is not useful, as the special "sparse plus low-rank" structure crucial to the efficiency of Soft-Impute no longer exist. In other words, the gain in convergence rate is more than compensated by the increase in per-iteration time complexity.

In this paper, we show that accelerating Soft-Impute is indeed possible while simultaneously preserving the "sparse plus low-rank" structure. Moreover, instead of using PROPACK to compute the (exact) partial SVD as in [Toh and Yun, 2010], we propose to use the power method [Halko et al., 2011] to obtain only an approximation of the dominant singular subspace. This is more efficient than PROPACK

and also allows warm-start, which is particularly useful because of the iterative nature of Soft-Impute. Since the SVT obtained is then only approximate, we use recent results in proximal gradient algorithms [Schmidt et al., 2011] to show that convergence can still be as fast as performing exact SVT in each iteration. Hence, the resultant algorithm has low per-iteration complexity and fast O(1/T 2 ) convergence rate.

The rest of the paper is organized as follows. Section 2 provides a brief review on the related work. The proposed accelerated inexact Soft-Impute algorithm is described in Section 3. Experimental results are presented in Section 4, and the last section gives some concluding remarks.

Notation

In the sequel, the transpose of vector/ matrix is denoted by the superscript ⊤. For a vector x, ‖x‖1 = ∑i |xi| is its ℓ1norm, and ‖x‖ = √∑i x2i its ℓ2-norm. For a matrix X, σ1 ≥ σ2 ≥ . . . are its singular values, tr(X) = ∑i Xii is its trace, ‖X‖1 = ∑i,j |Xij|, ‖X‖F = tr(X⊤X) is the Frobenius norm, and ‖X‖∗ = ∑i σi the nuclear norm. Moreover, I denotes the identity matrix. For a function f, we let f ′ be a subgradient in the subdifferential ∂f(x) = {g | f(y) ≥ f(x)+gT (y-x), ∀y}. When f is differentiable, we use∇f for its gradient. Finally, we use Ω⊥ to denote the complement

of a set Ω. For a matrix A, [PΩ(A)]ij = { Aij Ωij = 1

2 RelatedWork

2.1 Proximal Gradient Algorithms

Consider minimizing composite functions of the form:

where f, g are convex, f is smooth but g is possibly nonsmooth. The proximal gradient algorithm [Beck and Teboulle, 2009; Nesterov, 2013; Parikh and Boyd, 2014] generates a sequence of estimates {xt} as

where

and proxµg (·) is called the proximal operator. When f has ρLipschitz continuous gradient (i.e., ‖∇f(x1) -∇f(x2)‖ ≤ ρ‖x1 - x2‖) and a fixed stepsize µ ≤ 1/ρ is used, this algorithm converges with rate O(1/T), where T is the number of iterations [Parikh and Boyd, 2014]. Moreover, it can be accelerated by replacing (4) with

Here, several choices for θt can be used. The resultant accelerated proximal gradient (APG) algorithm converges with the optimal O(1/T 2 ) rate [Nesterov, 2013].

In the sequel, as our focus is on matrix completion, the variable x in (2) is a matrix X.

2.2 Soft-Impute

Soft-Impute is a state-of-the-art algorithm for large-scale matrix completion [Mazumder et al., 2010]. At iteration t, the missing values in O are filled in as

where Xt is the current iterate. Xt+1 is then generated as

where SVTλ(·) is the singular value thresholding (SVT) operator defined as follows.

Lemma 2.1 ([Cai et al., 2010]). Let the SVD of a matrix Z be UΣV ⊤. The solution to the optimization problem

is a low-rank matrix given by

where [(A)+]ij = max(Aij , 0).

To compute Xt+1 in (7), we thus need to first perform SVD on Zt. In general, obtaining the rank-k SVD of a m × n matrix Z takes O(mnk) time. Its most expensive steps are on computing matrix-vector multiplications of the form Zu and v⊤Z, where u ∈ Rn and v ∈ Rm.

To make Soft-Impute efficient, an important observation by Mazumder et al. is that Zt in (6) has a special "sparse plus low-rank" structure, namely that PΩ(O -Xt) is sparse and Xt is low-rank. Let the rank of Xt be r, and its SVD be UtΣtVt⊤. For any u ∈ Rn , Ztu can then be computed as

Constructing PΩ(O - Xt) takes O(r‖Ω‖1) time while computing PΩ(O - Xt)u takes O(‖Ω‖1) time, and computing UtΣt(Vt⊤u) takes O((m + n)r) time. Thus, to obtain the rank-k SVD of Zt, Soft-Impute needs only O((r +k)‖Ω‖1 + rk(m+n)) time. Moreover, its convergence rate is O(1/T).

In (1), take f(X) ≡ 12 ‖PΩ(X - O)‖2F and g(X) ≡ λ‖X‖∗. Note that ∇f(X) = PΩ(X -O), and ‖∇f(X1) - ∇f(X2)‖2F = ‖PΩ(X1-X2)‖2F ≤ ‖X1-X2‖2F . Hence, f is convex and has 1-Lipschitz-continuous gradient (i.e., ρ = 1); while g is also convex but nonsmooth. It is obvious that Zt in (4) is the same as that in (6). Using µ = 1/ρ = 1, it can be easily shown that proxµg (Zt) = proxg (Zt) in (3) is indeed equal to SVTλ(Zt). Hence, interestingly, Soft-Impute is also a proximal gradient algorithm [Tibshirani, 2010].

3 Accelerated Inexact Soft-Impute

Since Soft-Impute is a proximal gradient algorithm, it is natural to consider accelerating it with the schemes in Section 2.1. However, Tibshirani [2010] suggested that this is not useful, as the special "sparse plus low-rank" structure crucial to the efficiency of Soft-Impute will be lost. In this Section, we show that this structure can indeed be preserved, and thus acceleration is possible.

3.1 Accelerating Soft-Impute

To accelerate Soft-Impute, recall from Sections 2.1 and 2.2 that we have to compute proxg ( ˜ Zt) = SVTλ( ˜ Zt), where

The efficiency of Soft-Impute hinges on the "sparse plus low-rank" structure of Zt in (6), which allows matrix-vector multiplications of the form Ztu and v⊤Zt to be computed inexpensively. In the following, we show that ˜ Zt also has a similar structure.

Assume that Xt and Xt-1 have ranks rt and rt-1, and their SVDs are UtΣtV ⊤ t and Ut-1Σt-1V ⊤ t-1 , respectively. Using these SVDs, the sparse PΩ(O - Yt) can be constructed in O((rt + rt-1)‖Ω‖1) time. Using (10), for any u ∈ Rn , PΩ(O - Yt)u can be obtained in O(‖Ω‖1) time; and (1 + θt)Xtu - θtXt-1u = (1 + θt)UtΣt(Vt⊤u) - θtUt-1Σt-1(Vt-1⊤u) in O((m+ n)(rt-1 + rt)) time. The same applies to the computation of v⊤ ˜ Zt for any v ∈ Rm. Thus, the rank-k SVDof ˜ Zt can be obtained in O((rt+rt-1+ k)‖Ω‖1 + (rt-1 + rt)k(m + n)) time, which is only about twice that of Soft-Impute.

3.2 Approximating the SVT

Note from (8) that SVTλ( ˜ Zt) only involves the leading singular vectors of ˜ Zt (whose singular values are greater than or equal to λ), and thus only a partial SVD is needed. Let ˆ k be the number of such singular values. The power method (Algorithm 1) [Halko et al., 2011] is a fast and accurate algorithm for obtaining an approximation Q ∈ Rm×k (where k ≥ ˆ k) of such a subspace. Besides the power method, algorithms such as PROPACK [Larsen, 1998] have also been used [Toh and Yun, 2010]. However, the power method is more efficient than PROPACK [Wu and Simon, 2000], and it also allows warm-start, which is particularly useful because of the iterative nature of Soft-Impute.

Algorithm 1 PowerMethod( ˜ Z, R, ˜ ϵ). In step 3, QR(·) is the QR factorization.

Require: ˜ Z ∈ Rm×n , initial R ∈ Rk×n for warm-start, tolerance ˜ ϵ;

• 1: Initialize Y1 ← ˜ ZR;
• 2: for τ = 1, 2, . . . , do
• 3: Qτ+1 = QR(Yτ);
• 4: Yτ+1 = ˜ Z( ˜ Z ⊤Qτ+1);
• 5:
• if ‖Qτ+1Qτ+1⊤ -QτQτ ⊤‖F ≤ ˜ ϵ then
• 6: break;
• 7: end if
• 8: end for
• 9: return Qτ+1;

With this Q, SVT( ˜ Zt) can be obtained from a much smaller SVT as follows.

Proposition 3.1. Let Q ∈ Rm×k where k ≥ ˆ k, be orthogonal and contains the subspace spanned by the top ˆ k left singular vectors of ˜ Zt. Then, SVTλ( ˜ Zt) = Q SVTλ(Q⊤ ˜ Zt).

Algorithm 2 shows how to approximate SVTλ( ˜ Zt). Let the target (exact) rank-k SVD of ˜ Zt be ˜Uk ˜Σk ˜ V ⊤ k . Step 1 first approximates ˜Uk by the power method. In steps 2 to 5, a much smaller, and thus less expensive, (exact) SVTλ(Q⊤ ˜ Zt) is obtained from (8). Finally, SVTλ( ˜ Zt) is recovered using Proposition 3.1.

Algorithm 2 Algorithm to approximate the SVT of ˜ Zt (approx-SVT( ˜ Zt, R, λ, ˜ ϵ)).

Require: ˜ Zt ∈ Rm×n , R ∈ Rn×k , thresholds λ and ˜ ϵ;

• 1: Q = PowerMethod( ˜ Zt, R, ˜ ϵ);
• 2: [U,Σ, V ] = SVD(Q⊤ ˜ Zt);
• 3: U = {ui | σi > λ};
• 4: V = {vi | σi > λ};
• 5: Σ = (Σ-λI)+;
• 6: return QU,Σ and V .

The main computation is on the power method (step 1), whose most expensive operations, in turn, are on matrix multiplications (steps 1 and 4 of Algorithm 1). By using the special structure of ˜ Zt in (10) as discussed in Section 3.1, the time complexity of Algorithm 2 can be reduced from O(mnk) to O((rt + rt-1 + k)‖Ω‖1 + (rt-1 + rt)k(m + n)).

The following Proposition guarantees the quality of ˆXt = (QU)ΣV ⊤ returned from Algorithm 2.

Proposition 3.2. Let σt,k be the kth singular value of ˜ Zt , ηt = σt,k+1/σt,k < 1, αt, βt, γt be constants depending on ˜ Zt, and that the power method terminates after j iterations. Assume that k ≥ ˆ k, and ˜ ϵ ≥ αtηj t √1 + η2 t . Then,

where hλ‖·‖∗ (X; ˜ Zt) ≡ 1 2 ‖X- ˜ Zt‖2+λ‖X‖∗ is the objective in the proximal step.

3.3 Proposed Algorithm

The whole procedure is shown in Algorithm 3. The core steps are 8-11, where an approximate SVT is performed. Similar to [Lin et al., 2010; Toh and Yun, 2010], steps 9-10 use the column space of the last two iterations (Vt and Vt-1) to warmstart the power method. For further speedup, at step 5, we employ a continuation strategy as in [Toh and Yun, 2010; Mazumder et al., 2010], in which λt is initialized to a large value and then decreased gradually. Moreover, as in [O'Donoghue and Candes, 2012; Nesterov, 2013], the algorithm is restarted (at step 12) if F(X) starts to increase.

3.4 Convergence of the Inexact APGAlgorithm

Since ˆXt only approximates SVT( ˜ Zt), the proximal objective hλ‖·‖∗ (·; ˜ Zt) is not exactly minimized. The convergence of such inexact APG algorithms has been recently studied in [Schmidt et al., 2011]. In general, for problem (2), an inexact APG algorithm may commit two types of errors: (i) an error et in the calculation of ∇f(·), and (ii) an error εt in the proximal objective achieved by xt+1, i.e.,

Algorithm 3 Accelerated Inexact Soft-Impute (AIS-Impute). Require: partially observed matrix O, parameter λ, decay

parameter ν ∈ (0, 1), threshold ϵ; 1: [U0, λ0, V0] = rank-1 SVD(PΩ(O)); 2: initialize c = 1, ˜ ϵ0 = ‖PΩ(O)‖F , X0 = X1 = λ0U0V ⊤ 0 ; 3: for t = 1, 2, . . . do 4: ˜ ϵt = ν t ˜ ϵ0; 5: λt = ν t (λ0 -λ) +λ; 6: θt = (c -1)/(c +2); 7: Yt =Xt + θt(Xt -Xt-1); 8: ˜ Zt = Yt + PΩ(O - Yt); 9: Vt-1 = Vt-1 -Vt(Vt⊤Vt-1), remove zero columns; 10: Rt = QR([Vt, Vt-1]); 11: [Ut+1, Σt+1, Vt+1] = approx-SVT( ˜ Zt, Rt, λt, ˜ ϵt); 12: if F (Ut+1Σt+1V ⊤ t+1 ) > F (UtΣtV ⊤ t ) then 13: c = 1; 14: else 15: c = c + 1; 16: end if 17: if |F(Ut+1Σt+1V ⊤ t+1 ) - F (UtΣtV ⊤ t )| ≤ ϵ then 18: break; 19: end if 20: end for 21: return Xt+1 = Ut+1Σt+1V ⊤ t+1 .

The following Theorem shows that the convergence rate remains unchanged if the errors decrease at appropriate rates. Theorem 3.3 ([Schmidt et al., 2011]). Assume that (i) f is convex and has Lipschitz-continuous gradient; and (ii) g is convex and possibly nonsmooth. If ‖et‖ and √εt decrease as O(1/t2+δ ) for some δ > 0, the inexact APG algorithm still converges with rate O(1/T 2 ).

In the proposed Algorithm 3, there is no error in ∇f(·) and so et = 0. As for εt, we have from Proposition 3.2 that εt ≤ ηt 1-ηt ˜ ϵtβtγt. From step 4, ˜ ϵt = O(ν t ) where ν ∈ (0, 1). Thus, εt converges to zero at a linear rate and is faster than the O(1/t2+δ ) rate required in Theorem 3.3. This is formally stated in the following Theorem.

Theorem 3.4. Assume the conditions in Proposition 3.2, Algorithm 3 converges to the optimal solution of (1) with a rate of O(1/T 2 ).

4 Experiments

In this section, we perform experiments on both synthetic data (Section 4.1) and the entire MovieLens and Netflix recommendation data sets (Section 4.2).

4.1 Synthetic Data

We generate a m × m data matrix O = UV + G, where the elements of U ∈ Rm×5 , V ∈ R5×m are sampled i.i.d. from the normal distribution N(0, 1), and elements of G sampled from N(0, 0.05). A total of ‖Ω‖1 = 15mlog(m) random elements in O are observed. Half of them are used for training, and the other half as validation set for parameter tuning. Testing is performed on the unobserved (missing) elements.

The following variants of proximal gradient algorithms are compared:

• accelerated proximal gradient algorithm (denoted "APG") 1 [Ji and Ye, 2009; Toh and Yun, 2010];
• Soft-Impute 2 [Mazumder et al., 2010];
• the proposed accelerated inexact Soft-Impute algorithm (denoted "AIS-Impute") in Algorithm 3.

APG and AIS-Impute are in Matlab, while Soft-Impute is in R (and is called from Matlab). For performance evaluation, we use (i) the normalized mean squared error NMSE = √‖P ⊥ Ω (X - UV )‖F /‖P ⊥ Ω (UV )‖F , where X is the recovered matrix; and (ii) the rank of X. We vary m in {500, 1000, 1500, 2000}. Each experiment is repeated five times. Experiments are performed on a PC with Intel i7 CPU and 16GB RAM.

Results are shown in Table 1. As can be seen, all algorithms are equally good at recovering the missing matrix entries, but AIS-Impute is much faster. A more detailed timing comparison is in Figure 1. APG and AIS-Impute converge much faster than Soft-Impute w.r.t. the number of iterations (Figure 1(a)), as their convergence rates are O(1/T 2 ) rather than O(1/T). Because of its inexact proximal step, AISImpute has a slightly higher objective than APG. However, when measured against time (Figure 1(b)), APG is the slowest as it does not utilize the "sparse plus low-rank" structure. Overall, AIS-Impute is the fastest, as it has both low periteration complexity and fast O(1/T 2 ) convergence rate.

• (a) objective vs #iterations.
• (b) objective vs time.

Figure 1: Performance on synthetic data (with m = 2000).

4.2 Recommendation Data

MovieLens Data

The MovieLens data set 3 (Table 2) contains ratings (from 1 to 5) of different users on movies, and has been commonly used in matrix completion experiments [Mazumder et al., 2010; Hsieh and Olsen, 2014]. Following [Wang et al., 2014], we use 50% of the observed ratings for training, 25% for validation and the rest for testing. Besides the proximal gradient algorithms in Section 4.1, we also compare with other state-of-the-art non-proximal-gradient-based matrix completion algorithms, including

1 http://www.math.nus.edu.sg/∼mattohkc/NNLS.html

2 http://cran.r-project.org/web/packages/softImpute/index.html

3 http://grouplens.org/datasets/movielens/

• active subspace selection (denoted "active") 4 algorithm [Hsieh and Olsen, 2014]: In each iteration, this algorithm uses the power method to identify the active row and column subspaces, and then reduces the nuclear norm optimization problem to a smaller problem;
• a Frank-Wolfe based algorithm with local acceleration (denoted "boost") 5 [Zhang et al., 2012];
• a recent variant of Soft-Impute, which replaces the SVD in the soft-thresholding step by alternating least squares (denoted as "ALT-Impute") [Hastie et al., 2014] and
• a second-order trust-region algorithm 6 (denoted "TR") that alternates between fixed-rank optimization and rank-one updates [Mishra et al., 2013].

For performance evaluation, as in [Hsieh and Olsen, 2014; Mazumder et al., 2010], we use (i) the root mean squared error RMSE = √‖PΩ(X -O)‖2F /‖Ω‖1; and (ii) the rank of X. The experiment is repeated five times.

Results are shown in Table 3, and convergence of the objective and testing RMSE are shown in Figure 2. Again, all algorithms are equally good at recovering the missing matrix entries. In terms of speed, AIS-Impute is the fastest. ALTImpute has the same convergence rate as Soft-Impute, but is faster (than Soft-Impute) as it does not require performing SVD. As for boost, it only needs to perform a sparse rankone SVD in each iteration. However, much time is spent on maintaining the recovery matrix in factorized form and also local acceleration. TR is the slowest, as it has to solve many fixed-rank optimization problems.

In Proposition 3.2, we assume that the rank of Rt is not smaller than the rank of Ut+1 (step 11 of Algorithm 3). Figure 3 compares their values throughout the iterations. As can be seen, this assumption always holds and the two ranks gradually converge to the final recovered rank.

Netflix Data

In this Section, we demonstrate the speedup of AIS-Impute over ALT-Impute and Soft-Impute on the Netflix data set, which contains ratings of 480,189 users on 17,770 movies. 1% of the ratings matrix are observed. We randomly sample 50% of the observed ratings for training, and the rest for testing. As in [Mazumder et al. 2010], several choices of λ are used.

4 http://www.cs.utexas.edu/∼cjhsieh/nuclear active 1.1.zip

5 http://users.cecs.anu.edu.au/∼xzhang/GCG/

6 https://sites.google.com/site/bamdevm/codes/tracenorm

Results are shown in Table 4, and convergence of the objective and testing RMSEw.r.t. the running time are shown in Figure 4. As can be seen, all algorithms can recover the same RMSE and rank. However, AIS-Impute does not involve with SVD and has a better O(1/T 2 ) rate. Thus, it is again the fastest.

5 Conclusion

In this paper, we show that Soft-Impute, as a proximal gradient algorithm, can be accelerated without losing the "sparse plus low-rank" structure crucial to the efficiency of Soft-Impute. To further reduce the per-iteration time complexity, we proposed an approximate-SVT scheme based on the power method. Theoretical analysis shows that the proposed algorithm still enjoys the fast O(1/T 2 ) convergence rate. Extensive experiments on both synthetic and large recommendation data sets show that the proposed algorithm is much faster than the state-of-the-art.

Acknowledgment

This research was supported in part by the Research Grants Council of the Hong Kong Special Administrative Region (Grant No. 614012).