Greedy Learning of Generalized Low-Rank Models

1 Introduction

In many machine learning problems, the data samples can be naturally represented as low-rank matrices. For example, in recommender systems, the ratings matrix is low-rank as users (and items) tend to form groups. The prediction of unknown ratings is then a low-rank matrix completion problem [Cand`es and Recht, 2009]. In social network analysis, the network can be represented by a matrix with entries representing similarities between node pairs. Unknown links are treated as missing values and predicted as in matrix completion [Chiang et al., 2014]. Low-rank matrix learning also have applications in image and video processing [Cand`es et al., 2011], multitask learning [Argyriou et al., 2006], multilabel learning [Tai and Lin, 2012], and robust matrix factorization [Eriksson and van den Hengel, 2012].

The low-rank matrix optimization problem is NP-hard [Recht et al., 2010], and direct minimization is difficult. To alleviate this problem, one common approach is to factorize the target m⇥ n matrix X as a product of two low-rank matrices U and V , where U 2 Rm⇥r and V 2 Rn⇥r with r ⌧ min¶m↪ n♦. Gradient descent and alternating minimization are often used for optimization [Srebro et al., 2004; Eriksson and van den Hengel, 2012; Wen et al., 2012]. However, as the objective is not jointly convex in U and V , this approach can suffer from slow convergence [Hsieh and Olsen, 2014].

Another approach is to replace the matrix rank by the nuclear norm (i.e., sum of its singular values). It is known that the nuclear norm is the tightest convex envelope of the matrix rank [Cand`es and Recht, 2009]. The resulting optimization problem is convex, and popular convex optimization solvers such as the proximal gradient algorithm [Beck and Teboulle, 2009] and Frank-Wolfe algorithm [Jaggi, 2013] can be used. However, though convergence properties can be guaranteed, singular value decomposition (SVD) is required in each iteration to generate the next iterate. This can be prohibitively expensive when the target matrix is large. Moreover, nuclear norm regularization often leads to biased estimation. Compared to factorization approaches, the obtained rank can be much higher and the prediction performance is inferior [Mazumder et al., 2010].

Recently, greedy algorithms have been explored for lowrank optimization [Shalev-shwartz et al., 2011a; Wang et al., 2014]. The idea is similar to orthogonal matching pursuit (OMP) [Pati et al., 1993] in sparse coding. For example, the state-of-the-art in matrix completion is the rank-one matrix pursuit (R1MP) algorithm [Wang et al., 2014]. In each iteration, it performs an efficient rank-one SVD on a sparse matrix, and then greedily adds a rank-one matrix to the matrix estimate. Unlike other algorithms which typically require a lot of iterations, it only takes r iterations to obtain a rankr solution. Its prediction performance is also comparable or even better than others.

However, R1MP is only designed for matrix completion with the square loss. As recently discussed in [Udell et al., 2015], different loss functions may be required in different learning scenarios. For example, in link prediction, the presence or absence of a link is naturally represented by a binary variable, and the logistic loss is thus more appropriate. In robust matrix learning applications, the `1 loss or Huber loss can be used to reduce sensitivity to outliers [Cand`es et al., 2011]. While computationally R1MP can be used with these loss functions, its convergence analysis is tightly based on OMP (and thus the square loss), and cannot be easily extended.

This motivates us to develop more general greedy algorithms that can be used in a wider range of low-rank matrix learning scenarios. In particular, we consider low-rank matrix optimization problems of the form

where r is the target rank, and the objective f can be smooth or nonsmooth, (general) convex or strongly convex. The proposed algorithm is an extension of R1MP, and can be reduced to R1MP when f is the square loss. In general, when f is convex and Lipschitz-smooth, convergence is guaranteed with a rate of O(1◁T). When f is strongly convex, this is improved to a linear rate. When f is nonsmooth, we obtain a O(1◁pT ) rate for (general) convex objectives and O(log(T)◁T) for strongly convex objectives. Experiments on large-scale data sets demonstrate that the proposed algorithms are much faster than the state-of-the-art, while achieving comparable or even better prediction performance.

Notation: The transpose of vector / matrix is denoted by the superscript >. For matrix A = [Aij] 2 Rm⇥n (without loss of generality, we assume that m  n), its Frobenius norm is kAkF = qPi↪j A2ij , `1 norm is ♣♣X♣♣1 = Pi↪j ♣Xij♣, nuclear norm is kAk⇤ = Pi σi(A), where σi(A)'s are the singular values, and σmax(A) is its largest singular value. For two vectors x↪ y, the inner product hx↪ yi = Pi xiyi; whereas for two matrices A and B, hA↪Bi = Pi↪j AijBij. For a smooth function f, rf denotes its gradient. When f is convex but nonsmooth, g 2 ¶u ♣ f(y) ≥ f(x) + hu↪ x - yi♦ is its subgradient at x. Moreover, given ⌦ 2 ¶0↪ 1♦m⇥n , [P⌦(A)]ij = Aij if ⌦ij = 1, and 0 otherwise.

2 Review: Rank-One Matrix Pursuit

The rank-one matrix pursuit (R1MP) algorithm [Wang et al., 2014] is designed for matrix completion [Cand`es and Recht, 2009]. Given a partially observed m ⇥ n matrix O = [Oij ], indices of the observed entries are contained in the matrix ⌦ 2 ¶0↪ 1♦m⇥n , where ⌦ij = 1 if Oij is observed and 0 otherwise. The goal is to find a low-rank matrix that is most similar to O at the observed entries. Mathematically, this is formulated as the following optimization problem:

where r is the target rank. Note that the square loss has to be used in R1MP.

The key observation is that if X has rank r, it can be written as the sum of r rank-one matrices, i.e., X = Pri=1 ✓iuiv> i , where ✓i 2 R and kuik2 = kvik2 = 1. To solve (2), R1MP starts with an empty estimate. At the tth iteration, the (ut↪ vt) pair that is most correlated with the current residual Rt = P⌦(O - Xt-1) is greedily added. It can be easily shown that this (ut↪ vt) pair are the leading left and right singular vectors of Rt, and can be efficiently obtained from the rankone SVD of Rt. After adding this new utv> t basis matrix, all coefficients of the current basis can be updated as

Because of the use of the square loss, this is a simple leastsquares regression problem with closed-form solution.

To save computation, R1MP also has an economic variant. This only updates the combination coefficients of the current estimate and the rank-one update matrix as:

Algorithm 1 R1MP [Wang et al., 2014].

• 1: Initialize: X0 = 0;
• 2: for t = 1↪ ▷ ▷ ▷ ↪ T do
• 3: Rt = P⌦(O-Xt-1);
• 4: [ut↪ st↪ vt] = rank1SVD(Rt);
• 5: update coefficients using (3) (standard version), or (4) (economic version);
• 8: return XT .
• 6: Xt = Pti=1 ✓iuiv> i ; 7: end for

Note that each R1MP iteration is computationally inexpensive. Moreover, as the matrix's rank is increased by one in each iteration, only r iterations are needed in order to obtain a rank-r solution. It can also be shown that the residual's norm decreases at a linear rate, i.e., kRtk2F  γ t-1 kP⌦(O)k2F for some γ 2 (0↪ 1).

3 Low-RankMatrix Learning with Smooth Objectives

Though R1MP is scalable, it can only be used for matrix completion with the square loss. In this Section, we extend R1MP to problems with more general, smooth objectives. Specifically, we only assume that the objective f is convex and LLipschitz smooth. This will be further extended to nonsmooth objectives in Section 4.

Definition 1. f is L-Lipschitz smooth if f(X)  f(Y ) + hX - Y↪ rf (Y )i + L2 kX - Y k2F for any X↪ Y .

3.1 Proposed Algorithm

Let the matrix iterate at the tth iteration be Xt-1. We follow the gradient direction rf(Xt-1) of the objective f, and find the rank-one matrix M that is most correlated with rf (Xt-1):

Similar to [Wang et al., 2014], its optimal solution is given by utv> t , where (ut↪ vt) are the leading left and right singular vectors of rf(Xt-1). We then set the coefficient ✓t for this new rank-one update matrix to -st◁L, where st is the singular value corresponding to (ut↪ vt). Optionally, all the coefficients ✓1↪ ▷ ▷ ▷ ↪ ✓t can be refined as

As in R1MP, an economic variant is to update the coefficients as [θ✓1↪ ▷ ▷ ▷ ↪ θ✓t-1↪ ⇢], where θ and ⇢ are obtained as

The whole procedure, which will be called "greedy lowrank learning" (GLRL), is shown in Algorithm 2. Its economic variant will be called EGLRL. Obviously, on matrix completion problems with the square loss, Algorithm 2 reduces to R1MP.

Algorithm 2 GLRL for low-rank matrix learning with smooth convex objective f .

• 1: Initialize: X0 = 0;
• 2: for t = 1↪ ▷ ▷ ▷ ↪ T do
• 3: [ut↪ st↪ vt] = rank1SVD(rf(Xt-1));
• 5: (optional:) refine coefficients using (6) (standard version) or (7) (economic version);
• 4: Xt =Xt-1 - st L utv> t ;
• 6: end for
• 7: return XT .

Note that (6), (7) are smooth minimization problems (with  r and 2 variables, respectively). As the target matrix is low-rank, r should be small and thus (6), (7) can be solved inexpensively. In the experiments, we use the popular limited-memory BGFS (L-BFGS) solver [Nocedal and Wright, 2006]. Empirically, fewer than five L-BFGS iterations are needed. Preliminary experiments show that using more iterations does not improve performance.

Unlike R1MP, note that the coefficient refinement at step 5 is optional. Convergence results in Theorems 2 and 3 below still hold even when this step is not performed. However, as will be illustrated in Section 5.1, coefficient refinement is always beneficial in practice. It results in a larger reduction of the objective in each iteration, and thus a better rank-r model after running for r iterations.

3.2 Convergence

The analysis of R1MP is based on orthogonal matching pursuit [Pati et al., 1993]. This requires the condition 12 krf(X)k2F = f(X), which only holds when f is the square loss. In contrast, our analysis for Algorithm 2 here is novel and can be used for any Lipschitz-smooth f .

The following Proposition shows that the objective is decreasing in each iteration. Because of the lack of space, all the proofs will be omitted.

Proposition 1. If f is L-Lipschitz smooth,

where

If f is strongly convex, a linear convergence rate can be obtained.

Definition 2. f is θ-strongly convex if f(X) ≥ f(Y ) + hX - Y↪ rf (Y )i + θ2 kX - Y k2F for any X↪ Y .

Theorem 2. Let X⇤ be the optimal solution of (1). If f is θ-strongly convex and L-Lipschitz smooth, then

If f is only (general) convex, the following shows that Algorithm 2 converges at a slower O(1◁T) rate.

Theorem3. If f is (general) convex and L-Lipschitz smooth, then

The square loss in (2) is only general convex and 1Lipschitz smooth. From Theorem 3, one would expect GLRL to only have a sublinear convergence rate of O(1◁T) on matrix completion problems. However, our analysis can be refined in this special case. The following Theorem shows that a linear rate can indeed be obtained, which also agrees with Theorem 3.1 of [Wang et al., 2014].

Theorem 4. When f is the square loss, f(XT ) - f(X⇤)  (1-d21 )T [f(X0) - f(X⇤)].

3.3 Per-Iteration Time Complexity

The per-iteration time complexity of Algorithm 2 is low. Here, we take the link prediction problem in Section 5.1 as an example. With f defined only on the observed entries of the link matrix, rf is sparse, and computation of rf(Xt-1) in step 3 takes O(♣♣⌦♣♣1) time. The rank-one SVD on rf(Xt-1) can be obtained by the power method [Halko et al., 2011] in O(♣♣⌦♣♣1) time. Refining coefficients using (6) takes O(t♣♣⌦♣♣1) time for the tth iteration. Thus, the total periteration time complexity of GLRL is O(t♣♣⌦♣♣1). Similarly, the per-iteration time complexity of EGLRL is O(♣♣⌦♣♣1). In comparison, the state-of-the-art AIS-Impute algorithm [Yao and Kwok, 2015] (with a convergence rate of O(1◁T 2 )) takes O(r2 ♣♣⌦♣♣1) time in each iteration, whereas the alternating minimization approach in [Chiang et al., 2014] (whose convergence rate is unknown) takes O(r♣♣⌦♣♣1) time per iteration.

3.4 Discussion

To learn the generalized low-rank model, Udell et al. (2015) followed the common approach of factorizing the target matrix as a product of two low-rank matrices and then performing alternating minimization [Srebro et al., 2004; Eriksson and van den Hengel, 2012; Wen et al., 2012; Chiang et al., 2014; Yu et al., 2014]. However, this may not be very efficient, and is much slower than R1MP on matrix completion problems [Wang et al., 2014]. More empirical comparisons will be demonstrated in Section 5.1.

Similar to R1MP, the greedy efficient component optimization (GECO) [Shalev-shwartz et al., 2011a] is also based on greedy approximation but can be used with any smooth objective. However, GECO is even slower than R1MP [Wang et al., 2014]. Moreover, it does not have convergence guarantee.

4 Low-RankMatrix Learning with Nonsmooth Objectives

Depending on the application, different (convex) nonsmooth loss functions may be used in generalized low-rank matrix models [Udell et al., 2015]. For example, the `1 loss is useful in robust matrix factorization [Cand`es et al., 2011], the scalene loss in quantile regression [Koenker, 2005], and the hinge loss in multilabel learning [Yu et al., 2014]. In this Section, we extend the GLRL algorithm, with simple modifications, to nonsmooth objectives.

4.1 Proposed Algorithm

As the objective is nonsmooth, one has to use the subgradient gt of f(Xt-1) at the tth iteration instead of the gradient in Section 3. Moreover, refining the coefficients as in (6) or (7) will now involve nonsmooth optimization, which is much harder. Hence, we do not optimize the coefficients. To ensure convergence, a sufficient reduction in the objective in each iteration is still required. To achieve this, instead of just adding a rank-one matrix, we add a rank-k matrix (where k may be greater than 1). This matrix should be most similar to gt, which can be easily obtained as:

where ¶(u1↪ v1)↪ ▷ ▷ ▷ ↪ (uk↪ vk)♦ are the k leading left and right singular vectors of gt, and ¶s1↪ ▷ ▷ ▷ ↪ sk♦ are corresponding singular values. The proposed procedure is shown in Algorithm 3. The stepsize in step 3 is given by

where c1 ≥ 1◁θ and c2 > 0.

Algorithm 3 GLRL for low-rank matrix learning with nonsmooth objective f .

1: Initialize: X0 = 0 and choose ⌫ 2 (0↪ 1); 2: for t = 1↪ ▷ ▷ ▷ ↪ T do 3: set ⌘t as in (10); 4: compute subgradient gt of f(Xt-1), ht = 0; 5: for i = 1↪ 2↪ ▷ ▷ ▷ do 6: [ui↪ si↪ vi] = rank1SVD(gt - ht); 7: ht = ht + siuiv> i ; 8: if kgt - htk2F  ⌫kgt-1 - ht-1k2F then 9: break; 10: end if 11: end for 12: Xt =Xt-1 - ⌘tht; 13: end for 14: return XT .

4.2 Convergence

The following Theorem shows that when f is nonsmooth and strongly convex, Algorithm 3 has a convergence rate of O(log T ◁T ).

Theorem 5. Assume that f is θ-strongly convex, and kgtkF  b1 for some b1 (t = 1↪ ▷ ▷ ▷ ↪ T), then

where c1 is as defined in (10), and b2 is a constant (depending on X0, θ, ⌫ and c1).

When f is only (general) convex, the following Theorem shows that the rate is reduced to O(1◁pT).

Theorem 6. Assume that f is (general) convex, kgtkF  b1 for some b1 and kXt - X⇤kF  b3 for some b3 (t = 1↪ ▷ ▷ ▷ ↪ T), then

where c2 is as defined in (10).

For other convex nonsmooth optimization problems, the same O(log(T)◁T) rate for strongly convex objectives and and O(1◁pT) rate for general convex objectives have also been observed [Shalev-Shwartz et al., 2011b]. However, their analysis is for different problems, and cannot be readily applied to our low-rank matrix learning problem here.

4.3 Per-Iteration Time Complexity

To study the per-iteration time complexity, we take the robust matrix factorization problem in Section 5.2 as an example. The main computations are on steps 4 and 6. In step 4, since the subgradient gt is sparse (nonzero only at the observed entries), computing gt takes O(♣♣⌦♣♣1) time. At outer iteration t and inner iteration i, gt in step 6 is sparse and ht has low rank (equal to i - 1). Thus, gt - ht admits the so-called "sparse plus low-rank" structure [Mazumder et al., 2010; Yao et al., 2015]. This allows matrix-vector multiplications and subsequently rank-one SVD to be performed much more efficiently. Specifically, for any v 2 Rn , the multiplication (gt -ht)v takes only O(♣♣⌦♣♣1 +(i-1)n) time (and similarly for the multiplication u>(gt - ht) with any u 2 Rm), and rank-one SVD using the power method takes O(♣♣⌦♣♣1 + (i - 1)n) time. Assuming that it inner iterations are run at (outer) iteration t, it takes a total of O(it♣♣⌦♣♣1 + (it - 1)2n) time. Typically, it is small (empirically, usually 1 or 2).

In comparison, though the ADMM algorithm in [Lin et al., 2010] has a faster O(1◁T) convergence rate, it needs SVD and takes O(m2n) time in each iteration. As for the Wiberg algorithm [Eriksson and van den Hengel, 2012], its convergence rate is unknown and a linear program with mr + ♣♣⌦♣♣1 variables needs to be solved in each iteration. As will be seen in Section 5.2, this is much slower than GLRL.

5 Experiments

In this section, we compare the proposed algorithms with the state-of-the-art on link prediction and robust matrix factorization. Experiments are performed on a PC with Intel i7 CPU and 32GB RAM. All the codes are in Matlab.

5.1 Social Network Analysis

Given a graph with m nodes and an incomplete adjacency matrix O 2 ¶◦1♦m⇥m, link prediction aims to recover a low-rank matrix X 2 Rm⇥m such that the signs of Xij's and Oij's agree on most of the observed entries. This can be formulated as the following optimization problem [Chiang et al., 2014]:

where ⌦ contains indices of the observed entries. Note that (11) uses the logistic loss, which is more appropriate as Oij's are binary.

The objective in (11) is convex and smooth. Hence, we compare the proposed GLRL (Algorithm 2 with coefficient update step (6)) and its economic variant EGLRL (using coefficient update step (7)) with the following:

1. AIS-Impute 1 [Yao and Kwok, 2015]: This is an accelerated proximal gradient algorithm with further speedup based on approximate SVD and the special "sparse plus low-rank" matrix structure in matrix completion;
2. Alternating minimization ("AltMin") [Chiang et al., 2014]: This factorizes X as a product of two low-rank matrices, and then uses alternating gradient descent for optimization.

As a further baseline, we also compare with the GLRL variant that does not perform coefficient update. We do not compare with greedy efficient component optimization (GECO) [Shalev-shwartz et al., 2011a], matrix norm boosting [Zhang et al., 2012] and active subspace selection [Hsieh and Olsen, 2014], as they have been shown to be slower than AIS-Impute and AltMin [Yao and Kwok, 2015; Wang et al., 2014].

Experiments are performed on the Epinions and Slashdot data sets 2 [Chiang et al., 2014] (Table 1). Each row/column of the matrix O corresponds to a user (users with fewer than two observations are removed). For Epinions, Oij = 1 if user i trusts user j, and -1 otherwise. Similarly for Slashdot, Oij = 1 if user i tags user j as friend, and -1 otherwise.

Table 1: Signed network data sets used.

5▷0 ⇥ 105

As in [Wang et al., 2014], we fix the number of power method iterations to 30. Following [Chiang et al., 2014], we use 10-fold cross-validation and fix the rank r to 40. Note that AIS-Impute uses the nuclear norm regularizer and does not explicitly constrain the rank. We select its regularization parameter so that its output rank is 40. To obtain a rank-r solution, GLRL is simply run for r iterations. For AIS-Impute and AltMin, they are stopped when the relative change in the objective is smaller than 10-4 . The output predictions are binarized by thresholding at zero. As in [Chiang et al., 2014], the sign prediction accuracy is used as performance measure.

1 https://github.com/quanmingyao/AIS-impute

2 https://snap.stanford.edu/data/

Table 2 shows the sign prediction accuracy on the test set. All methods, except the GLRL variant that does not perform coefficient update, have comparable prediction performance. However, as shown in Figure 1, AltMin and AIS-Impute are much slower (as discussed in Section 3.3). EGLRL has the lowest per-iteration cost, and is also faster than GLRL.

Table 2: Testing sign prediction accuracy (%) on link prediction. The best and comparable results (according to the pairwise t-test with 95% confidence) are highlighted.

93.6◦0.1 84.4◦0.3

5.2 Robust Matrix Factorization

Instead of using the square loss, robust matrix factorization uses the `1 loss to reduce sensitivities to outliers [Cand`es et al., 2011]. This can be formulated as the optimization problem [Lin et al., 2010]:

Note that the objective is only general convex, and its subgradient is bounded ( k⌦kF ). Since there is no smooth component in the objective, AIS-Impute and AltMin cannot be used. Instead, we compare GLRL in Algorithm 3 (with ⌫ = 0▷99 and c2 = 0▷05) with the following:

1. Alternating direction method of multipliers (ADMM) 3 [Lin et al., 2010]: The rank constraint is replaced by

3 http://perception.csl.illinois.edu/matrix-rank/Files/ inexact alm rpca.zip

the nuclear norm regularizer, and ADMM [Boyd et al., 2011] is then used to solve the equivalent problem: minX↪Y P(i↪j)2⌦ ♣Xij - Oij ♣ + λkY k⇤ s.t. Y = X.

2. The Wiberg algorithm [Eriksson and van den Hengel, 2012]: It factorizes X into UV > and optimizes them by linear programming. Here, we use the linear programming solver in Matlab.

Experiments are performed on the MovieLens data sets 4 (Table 3), which have been commonly used for evaluating recommender systems [Wang et al., 2014]. They contain ratings ¶1↪ 2↪ ▷ ▷ ▷ ↪ 5♦ assigned by various users on movies. The setup is the same as in [Wang et al., 2014]. 50% of the ratings are randomly sampled for training while the rest for testing. The ranks used for the 100K, 1M, 10M data sets are 10, 10, and 20, respectively. For performance evaluation, we use the mean absolute error on the unobserved entries ⌦?:

where ˆX is the predicted matrix [Eriksson and van den Hengel, 2012]. Experiments are repeated five times with random training/testing splits.

Results are shown in Table 4. As can be seen, ADMM performs slightly better on the 100K data set, and GLRL is more accurate than Wiberg. However, ADMM is computationally expensive as SVD is required in each iteration. Thus, it cannot be run on the larger 1M and 10M data sets. Figure 2 shows the convergence of MABS with CPU time. As can be seen, GLRL is the fastest, which is then followed by Wiberg, and ADMM is the slowest.

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

Table 4: Testing MABS on the MovieLens data sets. The best results (according to the pairwise t-test with 95% confidence) are highlighted. ADMM cannot converge in 1000 seconds on the 1M and 10M data sets, and thus is not shown.

GLRL 0.724◦0.004 0.694◦0.001 0.683◦0.001

6 Conclusion

In this paper, we propose an efficient greedy algorithm for the learning of generalized low-rank models. Our algorithm is based on the state-of-art R1MP algorithm, but allows the optimization objective to be smooth or nonsmooth, general convex or strongly convex. Convergence analysis shows that the proposed algorithm has fast convergence rates, and is compatible with those obtained on other (convex) smooth/nonsmooth optimization problems. Specifically, on smooth problems, it converges with a rate of O(1◁T ) on general convex problems and a linear rate on strongly convex problems. On nonsmooth problems, it converges with a rate of O(1◁pT) on general convex problems and O(log(T)◁T) rate on strongly convex problems. Experimental results on link prediction and robust matrix factorization show that the proposed algorithm achieves comparable or better prediction performance as the state-of-the-art, but is much faster.

Acknowledgments

This research was supported in part by the Research Grants Council of the Hong Kong Special Administrative Region

(Grant 614513).