Matrix completion is the task of inferring the missing entries of a matrix given a subset of known entries. Typically, this is possible because the matrix to be completed has (at least approximately) low rank . This problem has witnessed a burst of activity, see e.g. candes2009exact ; candes2010power ; kmo10 , motivated by many applications such as collaborative filtering candes2009exact , quantum tomography gross2010quantum in physics, or the analysis of a covariance matrix candes2009exact .A commonly studied model for matrix completion assumes the matrix to be exactly low rank, with the known entries chosen uniformly at random and observed without noise. The most widely considered question in this setting is how many entries need to be revealed such that the matrix can be completed exactly in a computationally efficient way candes2009exact ; kmo10 . While our present paper assumes the same model, the main questions we investigate are different.
The first question we address in this paper is detectability, i.e. how many random entries do we need to reveal in order to be able to estimate the rank reliably. This question is motivated by the more generic problem of detecting structure (in our case, low rank) hidden in partially observed data. It is reasonable to expect the existence of a region where exact completion is hard or even impossible yet the rank estimation is tractable. A second question we address is what is the minimum achievable root-mean-square error (RMSE) in estimating the unknown elements of the matrix. In practice, even if exact reconstruction is not possible, having a procedure that provides a very small RMSE might be quite sufficient.
In this paper we propose an algorithm called MaCBetH that gives the best known empirical performance for the two tasks above when the rank is small. The rank in our algorithm is estimated as the number of negative eigenvalues of an associated Bethe Hessian matrix saade2014spectral ; saade2015spectral , and the corresponding eigenvectors are used as an initial condition for the local optimization of a cost function commonly considered in matrix completion (see e.g. kmo10 ). In particular, in the random matrix setting, we show that MaCBetH detects the rank of a large matrix from entries, where is a small constant, see Fig. 2, and as . The corresponding RMSE is evaluated empirically, and in the regime close to , it compares very favorably to existing approaches, in particular to OptSpace kmo10 .
This contribution is organized as follows. First, in Sec. II we define the problem and present generally our approach in the context of existing works. In Sec. III we describe our algorithm and motivate its construction via a spectral relaxation of the Hopfield model of neural network. Next, in Sec. IV
we show how the performance of the proposed spectral method can be analyzed using, in parts, results from spin glass theory and phase transitions, and rigorous results on the spectral density of large random matrices. Finally, in Sec.V we present numerical simulations that demonstrate the efficiency of MaCBetH.
Ii Problem definition and relation to other works
Let be a rank- matrix such that
where and are two (unknown) tall matrices. We observe only a small fraction of the elements of , chosen uniformly at random. We call the subset of observed entries, and the (sparse) matrix supported on whose nonzero elements are the revealed entries of . The aim is to reconstruct the rank matrix given . An important parameter which controls the difficulty of the problem is . In the case of a square matrix , this is the average number of revealed entries per line or column.
In our numerical examples and theoretical justifications we shall generate the low rank matrix , using tall matrices and with iid Gaussian elements, we call this the random matrix setting. The MaCBetH algorithm is, however, non-parametric and does not use any prior knowledge about and . The analysis we perform applies to the limit while and .
The matrix completion problem was popularized in candes2009exact who proposed nuclear norm minimization as a convex relaxation of the problem. The algorithmic complexity of the associated semidefinite programming is, however, . A low complexity procedure to solve the problem was later proposed by cai2010singular
and is based on singular value decomposition (SVD). A considerable step towards theoretical understanding of matrix completion from few entries was made inkmo10 who proved that with the use of trimming the performance of SVD-based matrix completion can be improved and a RMSE proportional to can be achieved. The algorithm of kmo10 is referred to as OptSpace, and empirically it achieves state-of-the-art RMSE in the regime of very few revealed entries.
OptSpace proceeds in three steps kmo10 . First, one trims the observed matrix by setting to zero all rows (resp. columns) with more revealed entries than twice the average number of revealed entries per row (resp. per column). Second, a singular value decompositions is performed on the matrix and only the first components are kept. When the rank is unknown it is estimated as the index for which the ratio between two consecutive singular values has a minimum. Third, a local minimization of the discrepancy between the observed entries and a low-rank estimate is performed. The initial condition for this minimization is given by the first
left and right singular vectors from the second step.
In this work we improve upon OptSpace by replacing the first two steps by a different spectral procedure that detects the rank and provides a better initial condition for the discrepancy minimization. Our method leverages on recent progress made in the task of detecting communities in the stochastic block model krzakala2013spectral ; saade2014spectral with spectral methods. Both in community detection and matrix completion, traditional spectral methods fail in the very sparse regime due to the existence of spurious large eigenvalues (or singular values) corresponding to localized eigenvectors krzakala2013spectral ; kmo10 . The authors of krzakala2013spectral ; saade2014spectral ; blm15 showed that using the non-backtracking matrix or the closely related Bethe Hessian as a basis for the spectral method in community detection provides reliable rank estimation and better inference performance. The present paper provides an analogous improvement for the matrix completion problem. In particular, we shall analyze the algorithm using tools from spin glass theory in statistical mechanics, and show that there exists a phase transition between a phase where it is able to detect the rank, and a phase where it is unable to do so.
Iii Algorithm and motivation
iii.1 The MacBetH algorithm
A standard approach to the completion problem (see e.g. kmo10 ) is to minimize the cost function
over and . This function is non-convex, and global optimization is hard. One therefore resorts to a local optimization technique with a careful choice of the initial conditions . In our method, given the matrix , we consider a weighted bipartite undirected graph with adjacency matrix
We will refer to the graph thus defined as . We now define the Bethe Hessian matrix to be the matrix with elements
where is a parameter that we will fix to a well-defined value depending on the data, and stands for the neighbors of in the graph . The MaCBetH algorithm that is the main subject of this paper is then, given the matrix , which we assume to be centered:
Numerically solve for the value of such that
Build the Bethe Hessian following eq. (4).
Compute all its negative eigenvalues and corresponding eigenvectors . is our estimate for the rank . Set (resp. ) to be the first lines (resp. the last lines) of the matrix .
Perform local optimization of the cost function (2) with rank and initial condition .
The function , in the first step, being an increasing function of , can be found efficiently, e.g. by dichotomy. Alternatively, in step 1 can be tuned in such a way that the number of negative eigenvalues of the Bethe Hessian is the largest possible. In step 2 we could also use the non-backtracking matrix weighted by , it was shown in saade2014spectral that the spectrum of the Bethe Hessian and the non-backtracking matrix are closely related. In the next section, we will motivate and analyze this algorithm (in the setting where was generated from elements-wise random and ) and show that in this case MaCBetH is able to infer the rank whenever . Fig. 1 illustrates the spectral properties of the Bethe Hessian that justify this algorithm: the spectrum is composed of a few informative negative eigenvalues, well separated from the bulk (which remains positive). This algorithm is computationally efficient as it is based on the eigenvalue decomposition of a sparse, symmetric matrix.
iii.2 Motivation from a Hopfield model
We shall now motivate the construction of the MaCBetH algorithm from a graphical model perspective and a spectral relaxation. Given the observed matrix from the previous section, we consider the following graphical model
where the and
are binary variables, andis a parameter controlling the strength of the interactions. This model is a (generalized) Hebbian Hopfield model hopfield1982neural on a bipartite sparse graph. To study it, we can use the standard Bethe approximation which is widely believed to be exact for such problems on large random graphs yedidia2001bethe ; MM2009 . In this approximation the means
and momentsof each variable are approximated by the parameters and that minimize the so-called Bethe free energy that reads
where , and are the degrees of nodes and in the graph .
Neural networks models such as eq. (6) have been extensively studied over the last decades (see e.g. amit1985spin ; wemmenhove2003finite ; castillo2004little ; MM2009 ; zhang2015nonbacktracking and references therein) and the phenomenology, that we shall review briefly here, is well known. In particular, for small enough, the global minimum of the Bethe free energy corresponds to the so-called paramagnetic state
As we increase , above a certain value , the model enters a retrieval phase, where the free energy has local minima correlated with the factors and . There are local minima, called retrieval states indexed by such that, in the large limit,
These retrieval states are therefore well-suited as initial conditions for the local optimization of eq. (2), and we expect their number to tell us the correct rank. Increasing above a critical value the system eventually enters a spin glass phase, marked by the appearance of many spurious minima.
It would be tempting to continue the Bethe approach and to derive the belief propagation equations, but we shall here instead consider a simpler spectral relaxation of the problem, following the same strategy as used in saade2014spectral ; saade2015spectral for graph clustering. First, we use the fact that the paramagnetic state (8) is always a stationary point of the Bethe free energy, for any value of mooij2004validity ; ricci2012bethe . In order to detect the retrieval states, we thus study its stability by looking for negative eigenvalues of the Hessian of the Bethe free energy evaluated at the paramagnetic state (8). At this point, the elements of the Hessian involving one derivative with respect to vanish, while the block involving two such derivatives is a diagonal positive definite matrix saade2014spectral ; mooij2004validity . The remaining part is the matrix called Bethe Hessian in saade2014spectral , and eigenvectors corresponding to its negative eigenvalues are thus expected to give an approximation of the retrieval states (9). The picture exposed in this section is summarized in figure 1. This motivates the use of the MaCBetH algorithm.
Note that a similar approach was used in zhang2015nonbacktracking to detect the retrieval states of a Hopfield model using the weighted non-backtracking matrix krzakala2013spectral , which linearizes the belief propagation equations rather than the Bethe free energy, resulting in a larger, non-symmetric matrix. The Bethe Hessian, while mathematically closely related, is also simpler to handle in practice.
Iv Analysis of performance in detection
We now show how the performance of MaCBetH can be analyzed, and the spectral properties of the matrix characterized using both tools from statistical mechanics and rigorous arguments.
iv.1 Analysis of the phase transition
We start by investigating the phase transition above which our spectral method will detect the correct rank. Let be random vectors with the same empirical distribution as the lines of and respectively. Using the statistical mechanics correspondence between the negative eigenvalues of the Bethe Hessian and the appearance of phase transitions in model (6), we can compute the values and where instabilities towards, respectively, the retrieval states and the spurious glassy states, arise. We have repeated the computations of amit1985spin ; wemmenhove2003finite ; castillo2004little ; zhang2015nonbacktracking in the case of model (6), using the cavity method MM2009 . We refer the reader interested in the technical details of the statistical mechanics approach to neural networks to wemmenhove2003finite ; castillo2004little ; zhang2015nonbacktracking .
Following a standard computation for locating phase transitions in the Bethe approximation (see e.g. MM2009 ; zdeborova2009statistical ), the stability of the paramagnetic state (8) towards these two phases can be monitored in terms of the two following parameters:
where the expectation is over the distribution of the vectors . The parameter controls the sensitivity of the paramagnetic solution to random noise, while measures its sensitivity to a perturbation in the direction of a retrieval state. and are defined implicitly as and , i.e. the value beyond which the perturbation diverges. The existence of a retrieval phase is equivalent to the condition , so that there exists a range of values of where the retrieval states exist, but not the spurious ones. If this condition is met, by setting in our algorithm, we ensure the presence of meaningful negative eigenvalues of the Bethe Hessian.
We define the critical value of such that if and only if . In general, there is no closed-form formula for this critical value, which is defined implicitly in terms of the functions and . We thus computed numerically using a population dynamics algorithm MM2009 and the results for are presented on Figure 2. Quite remarkably, with the definition , the critical value does not depend on the ratio , only on the rank .
In the limit of large and
it is possible to obtain a simple closed-form formula. In this case the observed entries of the matrix become jointly Gaussian distributed, and uncorrelated, and therefore independent. Expression (10) then simplifies to
We stress that the MaCBetH algorithm uses an empirical estimator (5) of this quantity to compute an approximation of purely from the revealed entries.
In the large regime, both decay to , so that we can further approximate
so that we reach the simple asymptotic expression, in the large limit, that , or equivalently . It is interesting to note that this result was obtained as the detectability threshold in completion of rank matrices from entries in the Bayes optimal setting in kabashima2014phase . Notice, however, that exact completion in the setting of kabashima2014phase is only possible for : clearly detection and exact completion are different phenomena.
iv.2 Computation of the spectral density
In this section, we show how the spectral density of the Bethe Hessian can be computed analytically on tree-like graphs such as those generated by picking uniformly at random the observed entries of the matrix . The spectral density is defined as
where the ’s are the eigenvalues of the Bethe Hessian. Using again the cavity method, It can be shown rogers2008cavitySym that the spectral density (in which potential delta peaks have been removed) is given by
where the are complex variables living on the vertices of the graph , which are given by:
where is the set of neighbors of . The are the (linearly stable) solution of the following belief propagation recursion:
The ingredients to derive this formula are to turn the computation of the spectral density into a marginalization problem for a graphical model on the graph , and then write the belief propagation equations to solve it. Quite remarkably, it has been shown RSA:RSA20313 that this approach leads to an asymptotically exact (and rigorous) description of the spectral density on Erdős-Rényi random graphs, which are locally tree-like in the limit where . We can again solve equation (18) numerically using the belief propagation algorithm. The results are shown on Fig. 1: the bulk of the spectrum is always positive.
We now demonstrate that for any value of , there exists an open set around where the spectral density vanishes. This justifies independently or choice for the parameter . The proof follows saade2014spectral and begins by noticing that is a fixed point of the recursion (18) for . Since this fixed point is real, the corresponding spectral density is . Now consider a small perturbation of this solution such that . The linearized version of (18) writes . The linear operator thus defined is a weighted version of the non-backtracking matrix of krzakala2013spectral . Its spectral radius is given by , where is defined in 10. In particular, for , , so that a straightforward application saade2014spectral of the implicit function theorem allows to show that there exists a neighborhood of such that for any , there exists a real, linearly stable fixed point of (18), yielding a spectral density equal to .
V Numerical tests
The algorithm was implemented in Julia BEKS14 , using the NLopt optimization package johnson2014nlopt for the minimization of the discrepancy (2). A matlab demo using the implementation of the limited-memory BFGS algorithm of schmidt2005minfunc is also available. Both demos can be downloaded from website .
Figure 3 illustrates the ability of the Bethe Hessian to infer the rank above the critical value in the limit of large size (see section IV.1). In Figure 4, we demonstrate the suitability of the eigenvectors of the Bethe Hessian as starting point for the minimization of the cost function (2). We compare the final RMSE achieved on the reconstructed matrix with other initializations of the optimization, including the largest singular vectors of the trimmed matrix kmo10 . MaCBetH systematically outperforms all the other choices of initial conditions. Remarkably, the performance achieved by MaCBetH with the inferred rank is essentially the same as the one achieved with an oracle rank. By contrast, estimating the correct rank from the (trimmed) SVD is more challenging. We note that for the choice of parameters we consider, trimming had a negligible effect. Along the same lines, OptSpace kmo10 uses a different minimization procedure, but from our tests we could not see any difference in performance due to that.
In this paper, we have presented MaCBetH, an algorithm for matrix completion that is efficient for two distinct, complementary, tasks: (i) it has the ability to estimate a finite rank reliably from fewer random entries than other existing approaches, and (ii) it gives lower root-mean-square reconstruction errors than its competitors. The algorithm is built around the Bethe Hessian matrix and leverages both on recent progresses in the construction of efficient spectral methods for clustering of sparse networks krzakala2013spectral ; saade2014spectral ; blm15 , and on the OptSpace approach kmo10 for matrix completion. Demos in Julia and matlab are available for download website .
The method presented here offers a number of possible future directions, including replacing the minimization of the cost function by a message-passing type algorithm, the use of different neural network models, or a more theoretical direction involving the computation of information theoretically optimal transitions for detectability.
The research leading to these results has received funding from the European Research Council under the European Union’s Framework Programme (FP/2007-2013/ERC Grant Agreement 307087-SPARCS).
- (1) E. J. Candès and B. Recht, “Exact matrix completion via convex optimization,” Foundations of Computational mathematics, vol. 9, no. 6, pp. 717–772, 2009.
- (2) E. J. Candès and T. Tao, “The power of convex relaxation: Near-optimal matrix completion,” Information Theory, IEEE Transactions on, vol. 56, no. 5, pp. 2053–2080, 2010.
- (3) R. H. Keshavan, A. Montanari, and S. Oh, “Matrix completion from a few entries,” Information Theory, IEEE Transactions on, vol. 56, no. 6, pp. 2980–2998, 2010.
- (4) D. Gross, Y.-K. Liu, S. T. Flammia, S. Becker, and J. Eisert, “Quantum state tomography via compressed sensing,” Physical review letters, vol. 105, no. 15, p. 150401, 2010.
A. Saade, F. Krzakala, and L. Zdeborová, “Spectral clustering of graphs with the bethe hessian,” inAdvances in Neural Information Processing Systems, 2014, pp. 406–414.
- (6) A. Saade, F. Krzakala, M. Lelarge, and L. Zdeborová, “Spectral detection in the censored block model,” IEEE International Symposium on Information Theory (ISIT2015), to appear, 2015.
- (7) J.-F. Cai, E. J. Candès, and Z. Shen, “A singular value thresholding algorithm for matrix completion,” SIAM Journal on Optimization, vol. 20, no. 4, pp. 1956–1982, 2010.
- (8) F. Krzakala, C. Moore, E. Mossel, J. Neeman, A. Sly, L. Zdeborová, and P. Zhang, “Spectral redemption in clustering sparse networks,” Proc. Natl. Acad. Sci., vol. 110, no. 52, pp. 20 935–20 940, 2013.
- (9) C. Bordenave, M. Lelarge, and L. Massoulié, “Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs,” 2015, arXiv:1501.06087.
- (10) J. J. Hopfield, “Neural networks and physical systems with emergent collective computational abilities,” Proc. Nat. Acad. Sci., vol. 79, no. 8, pp. 2554–2558, 1982.
- (11) J. S. Yedidia, W. T. Freeman, and Y. Weiss, “Bethe free energy, kikuchi approximations, and belief propagation algorithms,” Advances in neural information processing systems, vol. 13, 2001.
- (12) M. Mezard and A. Montanari, Information, Physics, and Computation. Oxford University Press, 2009.
- (13) D. J. Amit, H. Gutfreund, and H. Sompolinsky, “Spin-glass models of neural networks,” Physical Review A, vol. 32, no. 2, p. 1007, 1985.
- (14) B. Wemmenhove and A. Coolen, “Finite connectivity attractor neural networks,” Journal of Physics A: Mathematical and General, vol. 36, no. 37, p. 9617, 2003.
- (15) I. P. Castillo and N. Skantzos, “The little–hopfield model on a sparse random graph,” Journal of Physics A: Mathematical and General, vol. 37, no. 39, p. 9087, 2004.
- (16) P. Zhang, “Nonbacktracking operator for the ising model and its applications in systems with multiple states,” Physical Review E, vol. 91, no. 4, p. 042120, 2015.
- (17) J. M. Mooij and H. J. Kappen, “Validity estimates for loopy belief propagation on binary real-world networks.” in Advances in Neural Information Processing Systems, 2004, pp. 945–952.
- (18) F. Ricci-Tersenghi, “The bethe approximation for solving the inverse ising problem: a comparison with other inference methods,” J. Stat. Mech.: Th. and Exp., p. P08015, 2012.
- (19) L. Zdeborová, “Statistical physics of hard optimization problems,” acta physica slovaca, vol. 59, no. 3, pp. 169–303, 2009.
- (20) Y. Kabashima, F. Krzakala, M. Mézard, A. Sakata, and L. Zdeborová, “Phase transitions and sample complexity in bayes-optimal matrix factorization,” 2014, arXiv:1402.1298.
- (21) T. Rogers, I. P. Castillo, R. Kühn, and K. Takeda, “Cavity approach to the spectral density of sparse symmetric random matrices,” Phys. Rev. E, vol. 78, no. 3, p. 031116, 2008.
- (22) C. Bordenave and M. Lelarge, “Resolvent of large random graphs,” Random Structures and Algorithms, vol. 37, no. 3, pp. 332–352, 2010.
- (23) J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, “Julia: A fresh approach to numerical computing,” 2014, arXiv:1411.1607.
- (24) S. G. Johnson, “The nlopt nonlinear-optimization package,” 2014.
- (25) M. Schmidt, “minfunc: unconstrained differentiable multivariate optimization in matlab,” http://www.cs.ubc.ca/s̃chmidtm/Software/minFunc.html, 2005.
- (26) The matlab (http://github.com/alaa-saade/macbeth_matlab) and the Julia (http://github.com/alaa-saade/macbeth_julia) implementaions are distributed on github. They are also linked on their author’s webpage http://alaasaade.wordpress.com/ and are listed in the SPHINX group software page as well ( http://www.lps.ens.fr/k̃rzakala/WASP.html ).
- (27) D. C. Liu and J. Nocedal, “On the limited memory bfgs method for large scale optimization,” Mathematical programming, vol. 45, no. 1-3, pp. 503–528, 1989.
- (28) R. H. Keshavan, A. Montanari, and S. Oh, “Low-rank matrix completion with noisy observations: a quantitative comparison,” in 47th Annual Allerton Conference on Communication, Control, and Computing, 2009, pp. 1216–1222.