Improved Bounded Matrix Completion for Large-Scale Recommender Systems

1 Introduction

Matrix factorization and matrix completion [Koren et al., 2009a] are widely used for personalized recommender systems. Given partially observed ratings in the user-item matrix, the goal is to predict missing ratings by solving one of the following optimization problems:

where (1.2) is the convex form, and the nonconvex form (1.1) is equivalent to the convex form if k is chosen large enough. In the above problems, Y ∈ Rm×n is our observed rating matrix, Ω is the observed set, usually |Ω| ≪ mn. W ∈ Rm×k and H ∈ Rk×n are low rank matrices with rank k ≪ min(m, n) and PΩ(·) is the indicator function with the indicate set Ω. ‖ · ‖∗ denote the nuclear norm which is commonly used as a convex relaxation technique for the nonconvex rank penalty.

Despite the excellent performance on many real world problems, traditional matrix factorization formulations do not enforce any constraint to the prediction values. However, in many real world problems, the ratings are bounded within a certain region. For example, the ratings of MovieLens and Netflix datasets are bounded in the range of [0.5, 5], but the prediction of Eq (1.1) and (1.2) often output values out of this range. In those cases, it is reasonable to impose bounded constraints to the original optimization problems, leading to the following Bounded Matrix Completion (BMC) problem:

Note that we use rmin ≤ X ≤ rmax to denote the elementwise constraints rmin ≤ Xij ≤ rmax for all i, j throughout the paper. This is a hard optimization problem because there are totally mn constraints, so a naive implementation will easily require O(mn) memory, which is prohibitive for problems with more than 10,000 users/movies.

Recently, [R. Kannan and Park., 2012] proposed an efficient block coordinate descent algorithm for solving the nonconvex form (1.3) of the BMC problem. The resulting algorithm, named BMA, often outperforms traditional matrix completion algorithms on real datasets with bounded ratings. However, the BMA algorithm does not always converge to stationary points. Instead, it can easily stuck in non-stationary points, as we will discuss in Section 2. This unstable convergence behavior often leads to a performance drop in practice.

In this paper, we propose a new algorithm for minimizing the BMC objective function. Our algorithm is based on the ADMM framework [Boyd et al., 2011] for solving the convex form (1.4), so it is guaranteed to converge to the optimal solution. Moreover, with carefully-designed update rules, our algorithm does not encounter O(mn) space complexity, and can scale to large datasets with 10 million ratings. Experimental results on real world datasets show that our algorithm can reach a better solution and is also much faster than existing bounded matrix completion approaches. More interestingly, our algorithm often outperforms state-of-the-art matrix factorization in terms of prediction accuracy.

The rest of the paper is outlined as follows. We present related work and give an example to show that coordinate descent does not work for bounded matrix completion in Section 2. Our main algorithm is proposed and analyzed in Section 3, and the experimental results on real datasets are presented in Section 4.

2 RelatedWork

There are many existing methods for solving (1.1) and (1.2), [Funk, 2006] first proposed matrix factorization for recommender systems and used SGD to solve the optimization problem. [Koren, 2008] proposed a baseline estimate and incorporate it with matrix factorization to promote predict accuracy. See also a useful survey in [Koren et al., 2009b]. Later, several coordinate descent and SGD algorithms have been proposed to scale matrix factorization to very large datasets [Yu et al., 2012; Yun et al., 2014; Gemulla et al., 2011]; unfortunately, none of the above algorithms can handle bounded constraints. There are several other recent trend of matrix completion focusing on handling side information [Rao et al., 2015; Chiang et al., 2015] and noisy observations [Hsieh et al., 2015], and our algorithm can also be potentially used in those cases.

Recently, [R. Kannan and Park., 2012] proposed an algorithm called BMA to solve the bounded constraint matrix factorization problem (1.3), and they demonstrated that BMA can achieve lower Test RMSE compared with SGD (without bounded constraint). BMA is built on the idea of block coordinate descent. However, block coordinate descent assumes that the feasible area of W,H can be decomposed as Cartesian product, and the bounded constraints in (1.3) violates this assumption, so BMA lacks theoretical support and may not converge to stationary points.

Here we give an example that BMA cannot converge to stationary points. Consider the problem with Y = ( 1 0 0 1 ), rmin = 0 and rmax = 1, and the initial solution is W = ( 1 -1 ), H = ( 0 0 ). When we fix W, H cannot change, and vice versa. But at this time W,H are not optimal and it is also not a stationary point. If we don't impose bounded constraints, when we fix W, H will change to (0.5, 0.5), which eventually converges to optimal solution. So we can see that the bounded constraint makes it difficult to use block coordinate descent.

Here we use a small synthetic example to illustrate the potential convergence problem of BMA. We assume there are 1,000 rows and 1,000 columns with 49,153 ratings bounded between 0 and 5. We set the low rank k = 1 and initialize W ∈ Rm×1 with standard Gaussian distribution variables and H ∈ R1×n with all zeros.

From the simulation result, BMA stucks at the initial solution and is not able to escape. However, our proposed method BMC-ADMM could converge to global optimal and get competing predict accuracy under the same initialization. Furthermore, BMA is much slower than standard matrix factorization, whose time complexity is O(k2mn) in each iteration and the space needed is O(mn/nblk+|Ω|), where nblk is the number of "blocks" used in BMA. Numerical experiments also show that the run time will increase significantly if the number of blocks is too large. The large time and space complexity imply that BMA cannot scale to large datasets very well.

3 ProposedMethod

In this section, we describe our algorithm for the bounded matrix completion problem. We begin with the mathematical formulation under the ADMM framework [Boyd et al., 2011], then describe our algorithm in detail. We also discuss the convergence property of our method and analyze the time and space complexity of our proposed algorithm.

3.1 Mathematical Formulation and Convergence Properties

To apply alternating direction method of multipliers (ADMM) framework, we first split the loss, regularization, and constraints in the original problem (1.4), leading to the following equivalent form:

This looks like a standard ADMM problem and can be solved by alternately updating X, Z, W and the Lagrange multipliers corresponding to the constraints. Unfortunately, solving (3.1) by ADMM requires O(mn) space complexity since X is a dense and full rank m by n matrix during the whole optimization procedure until converging to the optimal solution.

To resolve this issue, we propose to further decompose X into PΩ(X) + PˆΩ(X), where PΩ(X) is a sparse matrix and PˆΩ(X) is stored as the product of two low rank matrices. Here we make a transformation and write the original X as Xnew + E, where Xnew := PΩ(X) is sparse and E := PˆΩ(X). This variable transformation leads to the following equivalent problem:

For simplicity, in (3.2) we use X to denote Xnew described above, and Z,E,W ∈ Rm×n are auxiliary matrices. The corresponding augmented Lagrangian is:

This is equivalent to the following scaled form:

where U1,U2 are Lagrangian multipliers, U1,U2 ∈ Rm×n , ρ1, ρ2 are penalty parameters, U ′1 = U1/ρ1,U ′2 = U2/ρ2. Note that the inequality constraint is not included in our augmented Lagrangian equation. The motivation of this is to facilitate the calculation of auxiliary matrix W, which is discussed and named as augmented partial Lagrangian in [Du et al., 2014].

Since the objective function (3.2) is closed, proper and convex, it can be solved by finding a saddle point of the augmented Lagrangian L(X,Z,W,U1,U2). Under the GaussSeidel framework, ADMM finds the saddle point using the following iterative procedure:

Next we discuss how to solve each subproblem in detail.

• Eq (3.3) can be solved in closed form by the elementwise soft-thresholding operator: Xk+1 i,j = 1 1 + ρ1 (Yi,j + ρ1Z k i,j -ρ1Ek i,j -ρ1U ′k 1i,j when (i, j) ∈ Ω and Xk+1 i,j = 0 when (i, j) /∈ Ω.
• To solve (3.4), we first notice that this subproblem is equivalent to

where A = ρ1 ρ1+ρ2 (Ek + Xk+1 + U ′k 1 ) + ρ2 ρ1+ρ2 (W k - U ′k 2 ) and ρ = ρ1 + ρ2. Therefore, the original problem has a closed form solution by soft-thresholding singular values: Z k+1 = Usoftλ/ρ(Λ)V T , where U,Λ, V is the SVD decomposition of A, softλ/ρ(Λ) = diag[(λ1 - λ/ρ)+, (λ2-λ/ρ)+, ..., (λn-λ/ρ)+]. Computing SVD is costly, but here we only need to compute singular values that are larger than the threshold λ/ρ, which can be solved iteratively by PROPACK [Larsen, 1998] or power iterations [Halko et al., 2011]. In this paper, we apply PROPACK to compute partial SVD decomposition efficiently. To use PROPACK, we have to specify the number of singular values greater than the threshold. We can either give a fixed number or dynamically predict this number in each iteration, which is presented in [Mazumder et al., 2010; Hsieh and Olsen, 2014]. In this paper, we use fixed rank k to compute partial SVD decomposition for simplicity.

• For eq (3.5), the subproblem has a simple closed form solution by setting Ek+1 i,j = 0 when (i, j) ∈ Ω, and Ek+1 i,j = Z k+1 i,j -Xk+1 i,j - U ′k 1i,j for (i, j) /∈ Ω. Since PˆΩ(Xk+1 ) = 0 and PˆΩ(U ′k 1 ) = 0, so E = PˆΩ(Z k+1 )
• For eq (3.6), W k+1 = ΠC(Z k+1 + U ′k 2 ), where ΠC is the projection into [rmin, rmax].

The resulting algorithm is summarized in Algorithm 1.

3.2 Complexity and Implementation Details

Time complexity

The most time consuming steps are (3.4) and (3.6). For (3.4), PROPACK uses iterative method to get truncated k SVD decomposition. For (3.6), the complexity is O(kmn) for checking all the elements in W k+1 and project back to the bounded constraints. Numerical experiments (refer to section 4.4) show that (3.6) is the slowest step and dominates the run time.

Space complexity

Here we use the same trick as [Lin et al., 2011] used to save memory usage, where Y,X,U1 are sparse matrices(O(|Ω|) complexity), Z is stored as the product of two low rank matrices(O(k(m + n)) complexity), Z = A1AT2 , A1 ∈ Rm×k , A2 ∈ Rn×k . Since E = PˆΩ(Z) and PΩ(E) = 0, we can express E by A1 and A2. U2 is also represented as a sparse matrix where its number of non-zero elements equals to the number of elements of W that are on the boundary. For W, W = ΠC(Z + U ′2 ), we need to store both A1,A2 and record the elements that are out of bound in order to represent W. Numerical experiments show that the number of elements out of bound is far less than |Ω|, for example, on movielens10m training set, |Ω| = 8 million, but mostly only 96, 073 elements are out of bound(nnz(U ′2 )) during the computation.

Algorithm 1 BMC-ADMM for Bounded Matrix Completion

Input Y, λ, k, rmin, rmaxρ1 > 0, ρ2 > 0, maxIter Initialize X, W, Z with baseline initialization Initialize U1, U2 with zeros matrices ∈ Rm×n for i = 1,2,...,maxIter do # Solve Eq (3.3) PΩ(Xi+1 ) = 1 1 + ρ1 PΩ(Y + ρ1(Z i - Ei - U i1 )) PˆΩ(Xi+1 ) = 0 # Solve Eq (3.4) ρ = ρ1 + ρ2 A = ρ1 ρ1 + ρ2 (Ei +Xi+1 +U i1 ) + ρ2 ρ1 + ρ2 (W i - U i2 ) (U, S, V ) = Partial svd(A, k) D = diag(diag(S(1 : k, 1 : k) - λ/ρ)+) Z i+1 = UDV T # Solve Eq (3.5) PˆΩ(Ei+1 ) = PˆΩ(Xi+1 - Z i+1 + U i1 ) PΩ(Ei+1 ) = 0 # Solve Eq (3.6) W i+1 = Z i+1 + U i2 W i+1 (W i+1 > rmax) = rmax W i+1 (W i+1 < rmin) = rmin # Update Scaled Lagrangian multipliers U i+1 1 = U i1 +Xi+1 -Z i+1 +Ei+1 U i+1 2 = U i2 +Z i+1 -W i+1 if stopping criterion is met then break end end Ouput rating matrixW

In summary, the total space complexity is O((m+ n)k + |Ω|) which is the same as standard matrix factorization methods.

Initialization

Initialization for BMC is non-trivial since W0H0 need to satisfy the bounded constraints. Here we discuss two initialization techniques.

The first one is random initialization: we generate two low rank matrices with random Gaussian variables W0 and H0, and then multiply them with a scalar α and adjust their first column and first row such that min (W0H0) = rmin and max (W0H0) = rmax. The detailed procedure is in Algorithm 2.

In Algorithm 2, the construction of AdjustTerm and α is based on the following equations:

{AdjustTerm + αvmin = rmin AdjustTerm + αvmax = rmax.

In this way, our initialized W0H0 is bounded in [rmin, rmax]. This is an improved version of the random initialization proposed in [R. Kannan and Park., 2012] since the random initialization implemented in [R. Kannan and Park., 2012] can only generate non-negative W0 and H0 and guarantee that the initialized W0H0 ∈ [0, rmax].

Algorithm 2 Random initialization

Input Y, k, rmin, rmax m, n = size(Y ) Initialized 2 random matrices W0 ∈ Rm×k and H0 ∈ Rk×n vmin = min(W0(:, 2 : end) ∗ H0(2 : end, :)) vmax = max(W0(:, 2 : end) ∗ H0(2 : end, :)) α = rmax - rmin vmax - vmin AdjustTerm = rminvmax - rmaxvmin vmax - vmin W0 = √αW0, H0 = √αH0 W0(:, 1) = AdjustTerm, H0(1, :) = 1 Ouput initialized matrices W0 and H0

Another initialization approach is called "baseline estimate". [Koren, 2008] presents that baseline estimate is a reasonable initialization for recommender systems, where the missing rating µi,j = µ+gi +hj , where µ is the global mean of ratings and gi and hj are corresponding bias terms for user i and movie j.

Note that these two initializations are also used by BMA, and we will use random initialization to compare the performance of BMA, CCD++ and our method in the following sections.

4 Experiments

In this section, we compare ADMM with existing matrix factorization methods on real datasets, in terms of efficiency, quality of the solution, and prediction accuracy. We include the following algorithms/implementations in our comparisons:

• BMC-ADMM: Our proposed method for solving the bounded matrix completion problem (Algorithm 1).
• BMA [R. Kannan and Park., 2012]: The coordinate descent algorithm for solving the BMC problem proposed in [R. Kannan and Park., 2012].
• CCD++ [Yu et al., 2012]: The coordinate descent solver for the traditional matrix completion problem (1.1).

We use five real world datasets to test the performance of the above matrix completion algorithms, and the detailed data statistics are listed in Table 1. We randomly split 80% as training data and 20% as testing data. For each algorithm and dataset, we choose the best regularization parameters from λ ∈ {0, 0.01, 0.1, 1, 10, 100} based on validation set. All the following experiments are conducted on an Intel i5-4590 3.30 GHz CPU with 16GB RAM.

Note that BMA is slow and requires large memory when data is large unless we let nblk to be very large to save memory, but this will make the BMA even slower. Therefore, we use small datasets or a subset of large datasets to compare BMC-ADMM, BMA and CCD++, and use full large dataset to compare BMC-ADMM and CCD++.

4.1 Comparison of BoundedMatrix Completion (BMC) Algorithms

In the first set of experiments, we compare BMC-ADMM with the state-of-the-art bounded matrix completion algorithm, BMA [R. Kannan and Park., 2012]. Although BMA solves the non-convex form (1.3) while we solve the convex form (1.4), it is well known that these two objective functions are equivalent since ‖X‖∗ = minU∈Rm×k ,V∈Rn×k :X=UV T 12 (‖U‖2F + ‖V ‖2F ) for sufficiently large k. Therefore, we first compare the objective function of BMC-ADMM with BMA in Figure 2.

We set ρ1 = ρ2 = 1 and try λ = 0.1 and 10, and compare the performance on all the datasets. In Figure 2, we observe that our algorithm BMC-ADMM can achieve lower objective function value than BMA, and the main reason is that BMA can stuck at non-stationary points, while our algorithm always converges to stationary points. In Figure 2, we also compare BMC-ADMM and BMA under different regularization parameter λ. Results show that BMC-ADMM can always find a solution with better objective function value.

To compare the prediction performance of these two methods, we list the test RMSE for all the 5 datasets in Table 2 using their own best parameter chosen by validation sets. The results indicate that BMC-ADMM can achieve a better test RMSE than BMA in most of the datasets except Flixster subset. Unfortunately, due to the large time complexity and huge memory requirement, BMA is not able to scale to the full Flixster dataset. In conclusion, we observe that BMC-ADMM is faster, more accurate, and more scalable than BMA.

4.2 Comparison with StandardMatrix Completion

Next we compare our algorithm, BMC-ADMM, with matrix completion algorithm without bounded constraints on large datasets. We choose the coordinate descent solver for matrix completion developed in [Yu et al., 2012] to compare with, and list the results on all the datasets in Table 2. To have a fair comparison, we tried different settings of rank k = 5, 10, 30 for all the algorithms, and select the best regularization parameter for each of them using the random validation set. The results in Table 2 indicates that BMC-ADMM outperforms other methods in most cases, and this means adding bounded constraints is really useful in practice if we want to achieve better prediction accuracy.

4.3 Scalability

To test the scalability of our BMC-ADMM, we run our algorithm on larger datasets, i.e. movielens10m, flixster and Jester. Result is presented in Figure 4. A larger rank k usually leads to a better test RMSE but slower convergence (see Figure 4 b,c), but sometimes it also leads to over-fitting, especially when the dataset has few users or movies (see Figure 4 a).

Experiments show that our algorithm is slower than standard matrix factorization method (without bounded constraint), but faster than BMA.

4.4 Time Spent on Each Part of BMC-ADMM

Table 3 present the detailed run time information of each step for our implementation. The update of W dominates our run time since its complexity is O(kmn). Note that for Eq. (3.5), E = PˆΩ(Z), so we don't need explicit computation for E once we get the update of Z. Therefore, there is no runtime for Eq. (3.5) in Table 3. The result presented here verified our computational complexity analysis, showing that computing the projection in (3.6) is the bottleneck of our algorithm.

5 Discussions and Conclusions

In this paper, we consider the bounded matrix factorization problem. We point out the convergence problem of the existing algorithm for Bounded Matrix Factorization BMA and propose a novel algorithm based on the ADMM framework. Experimental results showed that our approach has better convergence properties than BMA and can achieve lower Test RMSE values in most cases. In terms of run time and space complexity, our algorithm also outperforms BMA but takes more time than standard matrix factorization (without bounded constraint). To further speedup our algorithm, we can use parallel computing to solve subproblem (3.6) since it is the bottleneck of our algorithm. Since this step is a simple element-wise projection, we expect a linear speedup by parallelizing it.