1 Introduction
How to reconstruct signals exactly from very few measurements? This central question in signal processing has been extensively studied in the last few years and has triggered a fast emerging field of research, namely compressed sensing. Exact recovery from few measurements is actually possible if the signal is sparse in some representation domain. A related essential question has been recently considered for matrices: is it possible to reconstruct matrices exactly from very few observations? It appears that exact recovery is also possible in this setting if the matrix is lowrank. The problem of lowrank recovery from sparse observations is referred as the matrix completion problem. Several important realworld problems can be cast as a matrix completion problem, including remote sensing [26], system identification [17] and recommendation systems [27]. Throughout the paper, we will consider the recommendation system problem as an illustration of the matrix completion problem. Such systems have indeed become very common in many applications such as movie or product recommendation (e.g. Netflix, Amazon, Facebook, and Apple). The Netflix recommendation system tries to predict ratings of movies never seen by users. Collaborative filtering is widely used today to solve this problem [6], inferring recommendations by finding similar rating patterns and using them to complete missing values. This is typically a matrix completion problem where the unknown values of the matrix are computed by finding a lowrank matrix that fits the given entries.
How much information is needed for the exact recovery of lowrank matrices? In the case of random uniformly sampled entries without noise, Candès and Recht showed in [9] that, to guarantee perfect recovery, the number of observed entries must be larger than for matrices of rank (this bound has been refined more recently, see [24] and references therein). The case of noisy observations was studied in [8, 21], while a nonuniform sampling setting was considered in [25]. In this work, we propose to use additional information about rows and columns of the matrix to further improve the matrix completion solution.
In the standard matrix completion problem, rows and columns are assumed to be completely unorganized. However, in many realworld problems like the Netflix problem, there exist relationships between users (such as their age, gender, hobbies, education, etc) and movies (such as their genre, release year, actors, origin country, etc). This information can be taken advantage of, since people sharing the same tastes for a class of movies are likely to rate them similarly. We make use of graphs to encode relationships between users and movies and we introduce a new reconstruction model called matrix completion on graphs. Our main goal is to find a lowrank matrix that is structured by the proximities between users and movies. Introducing structure in sparse recovery problems is not new in the literature of compressed sensing [14, 1, 16], while similar structure inducing regularization has been proposed for factorized models for matrix completion [18]. Yet, introducing structures via graphs in the convex lowrank matrix recovery setting is novel to the best of our knowledge. We note that a large class of recommendation systems, called contentbased filtering, use graphs and clustering techniques to make predictions [15]. Along this line, our proposed methodology can be seen as a hybrid recommendation system that combines collaborative filtering (lowrank property) and contentbased filtering (graphs of users and movies).
We borrow ideas from the field of manifold learning [3, 4] and force the solution to be smooth on the manifolds of users and movies. We make the standard assumption that the graphs of users and movies are (nonuniform) discretizations of the corresponding manifolds. Forcing smoothness can be achieved through different regularizations. We use the popular Dirichlet/Laplacian energy, whose minimizing flow is the wellknown linear heat diffusion on manifold/graph, and show that the proposed model leads to a convex nonsmooth optimization problem. Convexity is a desired property for uniqueness of the solution. Nonsmoothness can be highly challenging, however, our problem belongs to the class of type optimization problems, for which several recently proposed efficient solvers exist [5, 11, 22]. The corresponding algorithm is derived in Section 4. It is tested on synthetic and real data [20] in Section 5.2.
2 Original matrix completion problem
The problem of matrix completion is to find the values of an matrix given a sparse set of observations . Problems of this kind are often encountered in collaborative filtering or recommender system applications, the most famous of which is the Netflix problem, in which one tries to predict the rating that users (columns of ) would give to films (rows of ), given only a few ratings provided by each user. A particularly popular model is to assume that the ratings are affected by a few factors, resulting in a lowrank matrix. This leads to the rank minimization problem
(1) 
where denotes the observed elements of . Problem (1) is NPhard. However, replacing with its convex surrogate known as the nuclear or trace norm [27] , where
are singular values of
, one obtains a semidefinite program(2) 
Under the assumption that is sufficiently incoherent, if the indices
are uniformly distributed and
is sufficiently large, the minimizer of (2) is unique and coincides with the minimizer of (1) [9, 24]. If in addition the observations are contaminated by noise, one can reformulate problem (2) as(3) 
where the data term in general depends on the type of noise assumed. If is the squared Frobenius norm ( here is the observations mask matrix, the Hadamard product), the distance between the solution of (3) and can be bounded by the norm of the noise [8].
One notable disadvantage of problems (23) is the assumption of a “good” distribution of the observed elements , which implies, in the movie rating example, that on average each user rates an equal number of movies, and each movie is rated by an equal number of users. In practice, this uniformity assumption is far from being realistic: for instance, in the Netflix dataset, the number of movie ratings of different users varies from to . When the sampling is nonuniform, the quality of the lower bound on deteriorates dramatically [25], from approximately constant number of observations per row in the former case, to an order of in the latter. In such settings, Salakhutdinov and Srebro [25] suggest using the weighted nuclear norm , where and are  and  dimensional row and columnmarginals of the distribution of observations, showing a significant performance improvement over the unweighted nuclear norm. Pathologically nonuniform sampling patterns, such as an entire row or column of missing, cannot be handled. Furthermore, in many situations the number of observations might be significantly smaller than the lower bounds.
3 Matrix completion on graphs
Low rank implies the linear dependence of rows/columns of . However, this dependence is unstructured. In many situations, the rows/columns of matrix possess additional structure that can be incorporated into the completion problem in the form of a regularization. In this paper, we assume that rows/columns of are given on vertices of graphs. In the Netflix example, the users (columns of ) are the vertices of a “social graph” whose edges represent e.g. friendship or similar tastes relations. Thus, it is reasonable to assume that connected users would give similar movie ratings, i.e., interpreting the ratings as an
dimensional vectorvalued function on the
vertices of the social graph, such a function would be smooth.More formally, let us be given the undirected weighted row graph with vertices and edges weighted with nonnegative weights represented by the matrix ; and respectively the column graph defined in the same way. Let be a matrix, which we will regard as a collection of dimensional column vectors denoted with subscripts , or of dimensional row vectors denoted with superscripts . Regarding the columns as a vectorvalued function defined on the vertices , the smoothness assumption implies that if . Stated differently, we want
(4) 
to be small, where , is the Laplacian of the column graph , and is the graph Dirichlet seminorm for columns. Similarly, for the rows we get a corresponding expression with the Laplacian of the row graph . These smoothness terms are added to the matrix completion problem as regularization terms (in the sequel, we treat the case where is the squared Frobenius norm),
(5) 
Relation to simultaneous sparsity models.
Low rank promoted by the nuclear norm implies sparsity in the space of outer products of the singular vectors, i.e., in the singular value decomposition
only a few coefficients are nonzero. Recent works [23] proposed imposing additional structure constraints, considering matrices that are simultaneously lowrank (i.e., sparse in the space of singular vectors outer products) and sparse (in the original representation). Our regularization can also be considered as a kind of simultaneously structured model. The column smoothness prior (4) makes the rows ofbe close to the eigenvectors of the column graph Laplacian
, i.e., each row of can be expressed as a linear combination of a few eigenvectors of . This can be interpreted as rowwise sparsity of in the column graph Laplacian eigenbasis. Similarly, the row smoothness prior results in columnwise sparsity of in the row graph Laplacian eigenbasis. Overall, the whole model (5) promotes simultaneous sparsity ofin the singular vectors outer product space, and row/columnwise sparsity in the respective Laplacian eigenspaces.
4 Optimization
Algorithm. Problems like (5) containing nondifferential terms cannot be tackled efficiently with a direct approach, while proximal based methods can be applied. We use the Alternating Direction Method of Multipliers (ADMM) that has seen great success recently [5] (other choices could include f.e. fast iterative soft thresholding [2]) by first introducing the equivalent splitting version of (5)
(6) 
This splitting step followed by an augmented Lagrangian method to handle the linear equality constraint is what constitutes ADMM. The success of ADMM for problems is mainly due to the fact that it does not require an exact solution for the iterative suboptimization problems, but rather an approximate solution. The augmented Lagrangian of (6) is . Both and are closed, proper and convex, and since we have no inequality constraints Slater’s conditions and therefore strong duality hold. Then and are primaldual optimal if is a saddle point of the augmented Lagrangian , i.e. or . ADMM finds a saddle point with the following iterative scheme
(7) 
(8) 
(9) 
Eventually the convergence of the proposed ADMM algorithm (7)(9) can be studied (and likely proved) with different mathematical approaches, for example [7].
Solving suboptimization problems. ADMM algorithms can be very fast as long as we can compute fast approximate solutions to the suboptimization problems, here (7) and (8).
Problem (7) requires finding that minimizes , i.e. , where is the proximal operator of defined as . In the case of the nuclear norm, there exists a closedform solution: , where are respectively the singular value decomposition (SVD) of , i.e. , and is the softthresholding operator defined for [7].
Problem (8) requires finding that minimizes , i.e. . Unlike Problem (7), there is no closedform solution to compute the proximal of as the solution is given by solving a linear system of equations. Precisely, the optimality condition of (8) is , which can be rewritten as as follows: , where is the columnstack vectorization operator, is the Kronecker product, and we use the formula . Also note that is symmetric positive semidefinite (s.p.s.d) as the Kronecker product of two s.p.s.d matrices, thus the conjugate gradient (CG) algorithm can be applied to compute a fast approximate solution of (8).
Computational complexity. The overall complexity of the algorithm is dominated by the computation of the nuclear proximal solutions by SVD, whose complexity is per iteration for [12]. The computational complexity of the CG algorithm is for NN graphs.
5 Numerical experiments
5.1 Synthetic ‘Netflix’ dataset
We start the evaluation of our matrix recovery model with a synthetic Netflixlike dataset, to study the behavior of model under controlled conditions. The artificial dataset is generated such that if fulfills two assumptions: (1) is lowrank and (2) its columns and rows are respectively smooth w.r.t. the column graph and the row graph . Figure 0(a) shows our synthetic dataset. It is inspired by the problem of movie recommendations as elements of are chosen to be integers from like in the Netflix prize problem. The matrix in Fig. 0(a) is noiseless, showing the ideal ratings for each pair of user and movie groups.
The row graph of the matrix is constructed as follows. The rows of are grouped into communities of different sizes. We connect nodes within a community using a
nearest neighbors graph and then add different amounts of erroneous edges, that is, edges between vertices belonging to different communities. The erroneous edges form a standard ErdősRényi graph with variable probability. We follow the same construction process for the column graph
that contains communities. For both graphs, binary edge weights are used. The intuition behind this choice of graphs is that users form communities of people with similar taste. Likewise, movies can be grouped according to their type, so that movies of the same group obtain similar ratings. The users graph is depicted in Fig 0(b), where nodes of the same community are clustered together. Note that matrix in Fig. 0(a) has rank equal to the minimum of user communities and the movie communities, in this case .5.1.1 Recovery quality versus number of observations
Two standard assumptions often made in the literature on matrix completion are that the observed elements of the matrix are sampled uniformly at random, and that the reconstructed matrix is perfectly lowrank (the case that we call noiseless).
Noiseless case. We test the performance of our method in this setting, comparing it to the standard nuclear normbased matrix completion (a particular case of our problem with ) and to a method that uses only the graphs (). We reconstruct the matrix using different levels of observed values and report the reconstruction root mean squared error (RMSE) on a fixed set of of the elements that was not observed. The result is depicted in Fig 1(a).
We use graphs with , , and of erroneous edges.
Noisy graphs alone (green lines) perform poorly compared to the nuclear norm reconstruction (blue line). However, when we use both graphs and nuclear norm (red lines), we obtain results that are better than any of the two alone.
Noisy case. We add noise to using a discretized Laplacian distribution. This type of noise models the human tendency to impulsively over or underrate a movie. In this case, the matrix that we try to reconstruct is close to lowrank, and the nuclear norm is still expected to perform well. As we see in Fig. 1(b) though, if we have highquality graphs (green line with erroneous edges), we can expect the same reconstruction quality of the nuclear norm regularization by using just half of the number of observations and only with the graph smoothness terms (green dashed line), that computationally is much cheaper to run. In this figure, the dashed black line designates the level of added noise in the data.
Note also that even if we use connectivity information of relatively bad quality (green dashed line with wrong edges), we can still benefit by combining the smoothness and the lowrank regularization terms (solid red line). Therefore the combination of nuclear and graph is robust to graph construction errors for low levels of observation. However, when the observation level is high enough ( for this specific size of matrices  note that this number may vary significantly depending on the matrix size), this benefit is lost (solid red line) and the nuclear norm regularizer (blue line) works better without the graph smoothness terms.
Nonuniform sampling. As noted in [25], the pattern of the observed values in real datasets does not usually follow a uniform distribution. In fact, the observations are such that the rating frequencies of users and movies closely follow power law distributions. In our experiment, we assume a simple generative process where users and movies are independently sampled from a power law, that is . This is a very sparse distribution with fixed expected number of observations, so we repeat this process identically times in order to control the overall density. Our final sampling is the logical OR operator of all these
‘epochs’, that follows the distribution
We find that this simple sampling scheme gives results close to the actual ratings of real datasets such as the MovieLens 10M that we use in the following.The results of our experiments for this setting are summarized in Fig. 2(a). Not surprisingly, all methods suffer from the nonuniformity of the sampling distribution. Still, the nuclear norm (blue line) crosses the line of a highquality graph ( green line) only after observations, while in the uniform case, Fig. 1(a), this happened for less than observed values. A similar behavior is exhibited for the noisy case, Fig. 2(b). There, the nuclear norm regularization quality is better than the mediumquality graph ( green line) only for more than observations, while in the uniform case, Fig. 1(b), the corresponding percentage was .
5.2 Movielens dataset
In this section, we report experiments on real data, which appear consistent with the results on the aforementioned artificial data. We work with the widely used MovieLens 10M dataset [20], containing ratings (‘stars’) from to (increments of ) given by 71,567 users for 10,677 movies. The density of the observations is . In our experiments, we use a subset of the original matrix for the reconstruction evaluation. This serves two purposes: firstly, we can choose an arbitrary density of the submatrix, and secondly, we can use ratings outside of it as features for the construction of the column and row graphs, as detailed below (see Figure 3(a)). Furthermore, the effect of nonuniformity is weaker.
The density of the observations is selected as follows. We sort the rows (users) and columns (movies) by order of increasing sampling frequency (Figure 3(b)). Then, users and movies are chosen to be close to the th and th percentile of their corresponding distributions.^{1}^{1}1Since the number of movies in the full matrix is much smaller than the number of users, we keep more frequently rating users in order to have a dense features matrix when we create the users graph. The resulting matrix has observed values that correspond to the ratings that a user has given to a movie.
After a row and column permutation, the original MovieLens 10M matrix is partitioned in blocks , where is the matrix that we use for our experiments (Figure 3(a)). We treat as the users feature matrix, as the movies feature matrix and discard the remaining matrix .
Graph construction.
Quality of graphs obviously plays an important role to our matrix recovery algorithm. Since a detailed analysis of how to construct good graphs is beyond the scope of this paper, we will resort to a simple, yet natural way of constructing the graphs for our setting, using the feature matrices and . We adapt the basic algorithm of [4] to our setting that contains missing values.
The distance we use between two users is the RMS distance between their commonly rated movies
where is the set of observed movie ratings for user (row) in and is the number of movies in that both users and have rated. We do the same to construct the movie distances from , that is, for each pair of movies we only take into account the ratings from users that have rated both. Note that distances between movies or between users, that take values from stars, share the same scale with the ratings and with the reconstruction error. Since the distances are all Euclidean, choosing the parameters of the graphs becomes more natural. The first choice we make is to use an neighborhood graph instead of a NN graph. To give weights to the edges, we use a Gaussian kernel, that is, . In the latter, denotes the minimum distance among all pairs of users and controls how fast the weights decay as distances increase. The transfer function used for the movies graph is plotted in Fig. 4(a), while the one for users is nearly identical.
We give weight values equal to for distances close to the minimum one (around stars), while the weights decay fast as the distance increases. We choose star, while is chosen so that the transfer function is already very close to for . This means that our model is equivalent to an infNN graph with the same exponential kernel. Note that the final reconstruction error is better than star in RMS, which justifies that distances that are smaller than that are trusted. We found that the results are indeed much better when a NN graph is not used. A possible explanation for this is that in that case a user that deviates a lot from the habits of other users would still have connections. These connections would not contribute positively in the recommendations quality regarding this user.
All this being said, it is here essential to emphasize that the graphs constructed for these experiments are not optimal. We foresee that the results presented in this paper can be further improved if one has access to detailed profile information about users and movies/products. This information is available to typical business companies that sell products to users.
Results. We apply a standard crossvalidation technique to evaluate the quality of our completion algorithm. For this purpose, the observations of the matrix are split into a fixed test set () and a varying size training set (from to ). We perform fold cross validation to select the parameters , and of our model (5) and only use the test set to evaluate the performance of the final models.
The recommendation error results are plotted in Fig. 4(b). The behavior of the algorithms is similar to the one exhibited by the noisy artificial data above (medium quality of graphs). For most observation levels our method combining nuclear norm and graph regularization (red line) clearly outperforms the rest. There are however two boundary phases that are noteworthy. When very few observations are available () there seems to be no benefit in adding the expensive nuclear norm term in the optimization problem, as the graph regularization alone (green line) performs best. On the other hand, for very dense observation levels () the nuclear norm (blue line) reaches the performance of the combined model. In general our combined model is very robust to observation sparsity, while the standard nuclear norm model performs worse even than the much cheaper graphsonly model for up to observations.
6 Conclusion
The main message of this work is that the standard lowrank matrix recovery problem can be further improved using similarity information about rows and columns. We solve an optimization problem seeking a lowrank solution that is structured by the proximity between rows and columns that form communities. As an application, our matrix completion model offers a new recommendation algorithm that combines the traditional collaborative filtering and contentbased filtering tasks into one unified model. The associated convex nonsmooth optimization problem is solved with a wellposed iterative ADMM scheme, which alternates between nuclear proximal operators and approximate solutions of linear systems. Artificial and real data experiments are conducted to study and validate the proposed matrix recovery model, suggesting that in reallife applications where the number of available matrix entries (ratings) is usually low and information about products and people taste is available, our model would outperform the standard matrix completion approaches.
Specifically, our model is robust to graph construction and to nonuniformly sampling of observations. Furthermore, it significantly outperforms the standard matrix completion when the number of observations is small.
The proposed matrix recovery algorithm can be improved in several ways. The effect of the nonuniformity of sampling matrix entries, as discussed in Sections 3 and 5.2, can be partially alleviated using a special weighting of the nuclear norm [25]. The nonuniform sampling of user data points and movie data points from the corresponding manifolds, which influences the quality of graph Laplacians, can also be corrected using special graph normalizations [10]. Furthermore, the optimization algorithm can be improved, firstly in terms of speed by using enhanced iterative schemes like [22]. Secondly in terms of scalability, either by using distributed schemes like [19] or by carrying out techniques from the recent work [13], which deals with nuclear norm for matrices with sizes much bigger than the Netflix dataset.
References
 [1] R.G. Baraniuk, V. Cevher, M.F. Duarte, and C. Hegde. Modelbased Compressive Sensing. IEEE Trans. Information Theory, 56(4):1982–2001, 2010.
 [2] Amir Beck and Marc Teboulle. A fast iterative shrinkagethresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences, 2(1):183–202, 2009.
 [3] M. Belkin and P. Niyogi. Laplacian eigenmaps and spectral techniques for embedding and clustering. In Proc. NIPS, 2001.
 [4] M. Belkin and P. Niyogi. Laplacian eigenmaps for dimensionality reduction and data representation. Neural Computation, 15(6):1373–1396, 2003.

[5]
S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein.
Distributed optimization and statistical learning via the alternating
direction method of multipliers.
Foundations and Trends in Machine Learning
, 3(1):1–122, 2011. 
[6]
J. S. Breese, D. Heckerman, and C. Kadie.
Empirical analysis of predictive algorithms for collaborative
filtering.
In
Proc. Conf. Uncertainty in Artificial Intelligence
, 1998.  [7] J.F. Cai, E. Candès, and Z. Shen. A singular value thresholding algorithm for matrix completion. SIAM J. Optimization, 20(4):1956–1982, 2010.
 [8] E. Candès and Y. Plan. Matrix completion with noise. Proc. IEEE, 98(6):925–936, 2010.
 [9] E. Candès and B. Recht. Exact matrix completion via convex optimization. FCM, 9(6):717–772, 2009.
 [10] R.R. Coifman and S. Lafon. Diffusion Maps. Applied and Computational Harmonic Analysis, 21(1):5–30, 2006.
 [11] P. L. Combettes and J.C. Pesquet. Proximal splitting methods in signal processing. In Fixedpoint algorithms for inverse problems in science and engineering, pages 185–212. Springer, 2011.
 [12] G. Golub and C. Van Loan. Matrix computations, volume 3. JHU Press, 2012.
 [13] C.J. Hsieh and P. Olsen. Nuclear norm minimization via active subspace selection. In Proc. ICML, pages 575–583, 2014.
 [14] J. Huang, Z. Zhang, and D. Metaxas. Learning with structured sparsity. In Proc. ICML, 2009.
 [15] Z. Huang, W. Chung, T.H. Ong, and H. Chen. A graphbased recommender system for digital library. In Proc. ACM/IEEECS joint conference on Digital libraries, 2002.
 [16] R. Jenatton, J.Y. Audibert, and F. Bach. Structured Variable Selection with SparsityInducing Norms. JMLR, 12:2777–2824, 2011.
 [17] Z. Liu and L. Vandenberghe. InteriorPoint Method for Nuclear Norm Approximation with Application to System Identification. SIAM Journal on Matrix Analysis and Applications, 31(3):1235–1256, 2009.
 [18] Hao Ma, Dengyong Zhou, Chao Liu, Michael R Lyu, and Irwin King. Recommender systems with social regularization. In Proceedings of the fourth ACM international conference on Web search and data mining, pages 287–296. ACM, 2011.
 [19] M. Mardani, G. Mateos, and G. Giannakis. Distributed nuclear norm minimization for matrix completion. In SPAWC, pages 354–358, 2012.
 [20] B. Miller, I. Albert, S. Lam, J. Konstan, and J. Riedl. MovieLens unplugged: Experiences with an occasionally connected recommender system. In Proc. Conf. Intelligent User Interfaces, 2003.
 [21] S. Negahban and M. J. Wainwright. Restricted strong convexity and weighted matrix completion: Optimal bounds with noise. JMLR, 13(1):1665–1697, 2012.
 [22] Y. Nesterov and A. Nemirovski. On firstorder algorithms for l 1/nuclear norm minimization. Acta Numerica, 22:509–575, 2013.
 [23] S. Oymak, A. Jalali, M. Fazel, Y. C Eldar, and B. Hassibi. Simultaneously structured models with application to sparse and lowrank matrices. arXiv:1212.3753, 2012.
 [24] B. Recht. A simpler approach to matrix completion. JMLR, 12:3413–3430, 2011.
 [25] R. Salakhutdinov and N. Srebro. Collaborative filtering in a nonuniform world: Learning with the weighted trace norm. In Proc. NIPS, 2010.

[26]
R. Schmidt.
Multiple emitter location and signal parameter estimation.
IEEE Trans. on Antennas and Propagation, 34(3):276–280, 1986.  [27] N. Srebro, J. Rennie, and T. Jaakkola. Maximummargin matrix factorization. In Proc. NIPS, volume 17, pages 1329–1336, 2004.
Comments
There are no comments yet.