Jointly Clustering Rows and Columns of Binary Matrices: Algorithms and Trade-offs

I Introduction

Data matrices exhibiting both row and column cluster structure, arise in many applications, such as collaborative filtering, gene expression analysis, and text mining. For example, in recommender systems, a rating matrix can be formed with rows corresponding to users and columns corresponding to items, and similar users and items form clusters. In DNA microarrays, a gene expression matrix can be formed with rows corresponding to patients and columns corresponding to genes, and similar patients and genes form clusters. Such row and column cluster structure is of great scientific interest and practical importance. For instance, the user and movie cluster structure is crucial for predicting user preferences and making accurate item recommendations . The patient and gene cluster structure reveals functional relations among genes and helps disease detection . In practice, we usually only observe a very small fraction of entries in these data matrices, possibly contaminated with noise, which obscures the intrinsic cluster structure. For example, in Netflix movie dataset, about $99\%$ of movie ratings are missing and the observed ratings are noisy .

In this paper, we study the problem of inferring hidden row and column cluster structure in binary data matrices from a few noisy observations. We consider a simple model introduced in for generating binary data matrix from underlying row and column clusters. In the context of movie recommender systems, our model assumes that users and movies each form equal-sized clusters. Users in the same cluster give the same rating to movies in the same cluster, where ratings are either $+1$ or $-1$ with $+1$ being “like” and $-1$ being “dislike”. Each rating is flipped independently with a fixed flipping probability less than $1/2$, modeling the noisy user behavior and the fact that users (movies) in the same cluster do not necessarily give (receive) identical ratings. Each rating is further erased independently with a erasure probability, modeling the fact that some ratings are not observed. Then, from the observed noisy ratings, we aim to exactly recover the underlying user and movie clusters, i.e., jointly cluster the rows and columns of the observed rating matrix.

The binary assumption on data matrices is of practical interest. Firstly, in many real datasets like Netflix dataset and DNA microarrays, estimation of entry values appears to be very unreliable, but the task of determining whether an entry is $+1$ or $-1$ can be done more reliably . Secondly, in recommender systems like rating music on Pandora or rating posts on sites such as Facebook and MathOverflow, the user ratings are indeed binary . The equal-sized assumption on cluster size is just for ease of presentation and can be relaxed to allow for different cluster sizes.

The hardness of our cluster recovery problem is governed by the erasure probability and cluster size. Intuitively, recovery becomes harder when the erasure probability increases, meaning fewer observations, and the cluster size decreases, meaning that clusters are harder to detect. The first goal of this paper is to understand when exact cluster recovery is possible or fundamentally impossible. Furthermore, our cluster recovery problem poses a computational challenge: An algorithm exhaustively searching over all the possible row and column cluster structures would have a time complexity exponentially increasing with the matrix dimension. The second goal of this paper is to understand how the computational complexity of our cluster recovery problem changes when increasingly more observations are available.

In this paper, our contributions are as follows. We first derive a lower bound on the minimum number of observations needed for exact cluster recovery as a function of matrix dimension and cluster size. Then we propose three algorithms with different runtimes and compare the number of observations needed by them for successful cluster recovery.

The first algorithm directly searches for the optimal clustering of rows and columns separately; it is combinatorial in nature and takes exponential-time but achieves the best statistical performance among the three algorithms in the noiseless setting.

By noticing that the underlying true rating matrix is a specific type of low-rank matrix, the second algorithm recovers the clusters by solving a nuclear norm regularized convex optimization problem, which is a popular heuristic for low rank matrix completion problems; it takes polynomial-time but has less powerful statistical performance than the first algorithm.

The third algorithm applies spectral clustering to the rows and columns separately and then performs a joint clean-up step; it has lower computational complexity than the previous two algorithms, but less powerful statistical performance. We believe that this is the first such performance guarantee for exact cluster recovery, with a growing number of clusters, using spectral clustering.

These algorithms are then compared with a simple nearest-neighbor clustering algorithm proposed in . Our analytical results show smooth time-data trade-offs: when increasingly more observations are available, one can gradually reduce the computational complexity by applying simpler algorithms while still achieving the desired performance. Such time-data trade-offs is of great practical interest for statistical learning problems involving large datasets .

The rest of the paper is organized as follows. In Section II, we discuss related work. In Section III, we formally introduce our model and main results. The lower bound is presented in Section IV. The combinatorial method, convex method, spectral method are studied in Section V, Section VI and Section VII, respectively. The proofs are given in Section VIII. The simulation results are presented in Section IX. Section X concludes the paper with remarks.

II Related work

In this section, we point out some connections of our model and results to prior work. There is a vast literature on clustering and we only focus on theoretical works with rigorous performance analysis. More detailed comparisons are provided after we present the theorems.

Much of the prior work on graph clustering, as surveyed in , focuses on graphs with a single node type, where nodes in the same cluster are more likely to have edges among them. A low-rank plus sparse matrix decomposition approach is proved to exactly recover the clusters with the best known performance guarantee in . The same approach is used to recover the clusters from a partially observed graph in . A spectral method for exact cluster recovery is proposed and analyzed in with the number of clusters fixed. More recently, proved an upper bound on the number of nodes “mis-clustered” by a spectral clustering algorithm in the high dimensional setting with a growing number of clusters.

In contrast to the above works, in our model, we have a labeled bipartite graph with two types of nodes (rows and columns). Notice that there are no edges among nodes of the same type and cluster structure is defined for the two types separately. In this sense, our cluster recovery problem can be viewed as a natural generalization of graph clustering problem to labeled bipartite graphs. In fact, our second algorithm via convex programming is inspired by the work .

A model similar to ours but with a fixed number of clusters has been considered in , where the spectral method plus majority voting is shown to approximately predict the rating matrix. However, our third algorithm via spectral method is shown to achieve exact cluster and rating matrix recovery with a growing number of clusters. This is the first theoretical result on spectral method for exact cluster recovery in with a growing number of clusters to our knowledge.

II-B Biclustering

Biclustering tries to find sub-matrices (which may overlap) with particular patterns in a data matrix. Many of the proposed algorithms are based on heuristic searches without provable performance guarantees. Our cluster recovery problem can be viewed as a special case where the data matrix consists of non-overlapping sub-matrices with constant binary entries, and our paper provides a thorough study of this special biclustering problem. Recently, there is a line of work studying another special case of biclustering problem, which tries to detect a single small submatrix with elevated mean in a large fully observed noisy matrix . Interesting statistical and computational trade-offs are summarized in .

II-C Low-rank matrix completion

Under our model, the underlying true data matrix is a specific type of low-rank matrix. If we recover the true data matrix, we immediately get the user (or movie) clusters by assigning the identical rows (or columns) of the matrix to the same cluster. In the noiseless setting with no flipping, the nuclear norm minimization approach can be directly applied to recover the true data matrix and further recover the row and column clusters. Alternate minimization is another popular and empirically successful approach for low-matrix completion . However, it is harder to analyze and the performance guarantee is weaker than nuclear norm minimization . In the low noise setting with the flipping probability restricting to be a small constant, the low-rank plus sparse matrix decomposition approach can be applied to exactly recover data matrix and further recover the row and column clusters.

The performance guarantee for our second algorithm via convex programming is better than these previous approaches and it allows the flipping probability to be any constant less than $1/2$. Moreover, our proof turns out to be much simpler. The recovery of our true data matrix can also be viewed as a specific type of one-bit matrix completion problem recently studied in . However, only focuses on approximate matrix recovery and the results there cannot be used to recover row and column clusters.

III Model and Main Results

In this section, we formally state our model and main results.

Our model is described in the context of movie recommender systems, but it is applicable to other systems with binary data matrices having row and column cluster structure. Consider a movie recommender system with $n$ users and $n$ movies. Let $R$ be the rating matrix of size $n\times n$ where $R_{ij}$ is the rating user $i$ gives to movie $j$. Assume both users and movies form $r$ clusters of size $K=n/r$. Users in the same cluster give the same rating to movies in the same cluster. The set of ratings corresponding to a user cluster and a movie cluster is called a block. Let $B$ be the block rating matrix of size $r\times r$ where $B_{kl}$ is the block rating user cluster $k$ gives to movie cluster $l$. Then the rating $R_{ij}=B_{kl}$ if user $i$ is in user cluster $k$ and movie $j$ is in movie cluster $l$. Further assume that entries of $B$ are independent random variables which are $+1$ or $-1$ with equal probability. Thus, we can imagine the rating matrix as a block-constant matrix with all the entries in each block being either $+1$ or $-1$. Observe that if $r$ is a fixed constant, then users from two different clusters have the same ratings for all movies with some positive probability, in which case it is impossible to differentiate between these two clusters. To avoid such situations, assume $r$ is at least $\Omega(\log n)$.

Suppose each entry of $R$ goes through an independent binary symmetric channel with flipping probability $p<1/2$, representing noisy user behavior, and an independent erasure channel with erasure probability $\epsilon$, modeling the fact that some entries are not observed. The expected number of observed ratings is $m=n^{2}(1-\epsilon)$. We assume that $p$ is a constant throughout the paper and $\epsilon$ could converge to $1$ as $n\to\infty$. Let $R^{\prime}$ denote the output of the binary symmetric channel and $\Omega$ denote the set of non-erased entries. Let $\widehat{R}_{ij}=R^{\prime}_{ij}$ if $(i,j)\in\Omega$ and $\widehat{R}_{ij}=0$ otherwise. The goal is to exactly recover the row and column clusters from the observation $\widehat{R}$.

III-B Main Results

The main results are summarized in Table I. Note that these results do not explicitly depend on $p$. In fact, as $p$ is assumed to be a constant strictly less than $1/2$, it affects the results by constant factors.

The parameter regime where exact cluster recovery is fundamentally impossible for any algorithm is proved in Section IV. The combinatorial method, convex method and spectral method are studied in Section V, Section VI and Section VII, respectively. We only analyze the combinatorial method in the noiseless case where $p=0$, but we believe similar result is true for the noisy case as well. The parameter regime in which the convex method succeeds is obtained by assuming that a technical conjecture holds, which is justified through extensive simulation. The parameter regime in which the spectral method succeeds is obtained for the first time for exact cluster recovery with a growing number of clusters. The nearest-neighbor clustering algorithm was proposed in . It clusters the users by finding the $K-1$ most similar neighbors for each user. The similarity between user $i$ and $i^{\prime}$ is measured by the number of movies with the same observed rating, i.e.,

The number of observations needed for successful cluster recovery can be derived from the corresponding parameter regime using the identity $m=n^{2}(1-\epsilon)$ as shown in Table I. For better illustration, we visualize our results in Figure 1. In particular, we take $\log(m/n)$ as $x$-axis and $\log K$ as $y$-axis and normalize both axes by $\log n$. Since exact cluster recovery becomes easy when the number of observations $m$ and cluster size $K$ increase, we expect that exact cluster recovery is easy near $(1,1)$ and hard near $(0,0)$.

From Figure 1, we can observe interesting trade-offs between algorithmic runtime and statistical performance. In terms of the runtime, the combinatorial method is exponential, while the other three algorithms are polynomial. In particular, the convex method can be casted as a semidefinite programming and solved in polynomial-time. For the spectral method, the most computationally expensive step is the singular value decomposition of the observed data matrix which can always be done in time $O(n^{3})$ and more efficiently when the observed data matrix is sparse. It is not hard to see that the time complexity for the nearest-neighbor clustering algorithm is $O(n^{2}r)$ and more careful analysis reveals that its time complexity is $O(mr)$. On the other hand, in terms of statistical performance, the combinatorial method needs strictly fewer observations than the other three algorithms when there is no noise, and the convex method always needs fewer observations than the spectral method. It is somewhat surprising to see that the simple nearest-neighbor clustering algorithm needs fewer observations than the more sophisticated convex method when the cluster size $K$ is $O(\sqrt{n})$.

In summary, we see that when more observations available, one can apply algorithms with less runtime while still achieving exact cluster recovery. For example, consider the noiseless case with cluster size $K=n^{0.8}$, the number of observations per user required for cluster recovery by the combinatorial method, convex method, spectral method and nearest-neighbor clustering algorithm are $\Omega(n^{0.1})$, $\Omega(n^{0.2})$, $\Omega(n^{0.4})$ and $\Omega(n^{0.5})$, respectively. Therefore, when the number of observations per user increases from $\Omega(n^{0.1})$ to $\Omega(n^{0.5})$, one can gradually reduces the computational complexity from exponential-time to polynomial-time as low as $O(n^{1.7})$.

The main results in this paper can be easily extended to the more general case with $n_{1}$ rows and $n_{2}=\Theta(n_{1})$ columns and $r_{1}$ row clusters and $r_{2}=\Theta(r_{1})$ column clusters. The sizes of different clusters could vary as long as they are of the same order. Likewise, the flipping probability $p$ and the erasure probability $\epsilon$ could also vary for different entries of the data matrix as long as they are of the same order. Due to space constraints, such generalizations are omitted in this paper.

III-C Notations

Throughout the paper, we say that an event occurs “a.a.s.” or “asymptotically almost surely” when it occurs with a probability which tends to one as $n$ goes to infinity.

IV Lower Bound

Fix $0<\delta<1$. If $nK^{2}(1-\epsilon)^{2}<\delta$, then with probability at least $1-\delta$, it is impossible for any algorithms to recover the user clusters or movie clusters.

Intuitively, Theorem 1 says that when the erasure probability is high and the cluster size is small that $nK^{2}(1-\epsilon)^{2}=O(1)$, the observed rating matrix $\widehat{R}$ does not carry enough information to distinguish between different possible cluster structures.

V Combinatorial Method

In this section, we study a combinatorial method which clusters users or movies by searching for a partition with the least total number of “disagreements”. We describe the method in Algorithm 1 for clustering users only. Movies are clustered similarly. The number of disagreements $D_{ii^{\prime}}$ between a pair of users $i,i^{\prime}$ is defined as the number of movies satisfying that: The two ratings given by users $i,i^{\prime}$ are both observed and the observed two ratings are different. In particular, if for every movie, the two ratings given by users $i,i^{\prime}$ are not observed simultaneously, then $D_{ii^{\prime}}=0$.

The idea of Algorithm 1 is to reduce the problem of clustering both users and movies to a standard user clustering problem without movie cluster structure. In fact, this algorithm looks for the optimal partition of the users which has the minimum total in-cluster distance, where the distance between two users is measured by the number of disagreements between them. The following theorem shows that such simple reduction does not achieve the lower bound given in Theorem 1. The optimal algorithm for our cluster recovery problem might need to explicitly make use of both user and movie cluster structures.

If $nK(1-\epsilon)^{2}\leq\frac{1}{4}$, then with probability at least $3/4$, Algorithm 1 cannot recover user and movie clusters.

Next we show that the above necessary condition for the combinatorial method is also sufficient up to a logarithmic factor when there is no noise, i.e., $p=0$. We suspect that the theorem holds for the noisy setting as well, but we have not yet been able to prove this.

If $p=0$ and $nK(1-\epsilon)^{2}>C\log n$ for some constant $C$, then a.a.s. Algorithm 1 exactly recovers user and movie clusters.

This theorem is proved by considering a conceptually simpler greedy algorithm that does not need to know $K$. After computing the number of disagreements for every pair of users, we search for a largest set of users which have no disagreement between each other, and assign them to a new cluster. We then remove these users and repeat the searching process until there is no user left. In the noiseless setting, the $K$ users from the same true cluster have no disagreement between each other. Therefore, it is sufficient to show that, for any set of $K$ users consisting of users from more than one cluster, they have more than one disagreement with high probability under our assumption.

VI Convex Method

In this section, we show that the rating matrix $R$ can be exactly recovered by a convex program, which is a relaxation of the maximum likelihood (ML) estimation. When $R$ is known, we immediately get the user (or movie) clusters by assigning the identical rows (or columns) of $R$ to the same cluster.

Let $\mathcal{Y}$ denote the set of binary block-constant rating matrix with $r^{2}$ blocks of equal size. As the flipping probability $p<1/2$, Maximum Likelihood (ML) estimation of $R$ is equivalent to finding a $Y\in\mathcal{Y}$ which best matches the observation $\widehat{R}$:

Since $|\mathcal{Y}|=\Omega(e^{n})$, solving (1) via exhaustive search takes exponential-time. Observe that $Y\in\mathcal{Y}$ implies that $Y$ is of rank at most $r$. Therefore, a natural relaxation of the constraint that $Y\in\mathcal{Y}$ is to replace it with a rank constraint on $Y$, which gives the following problem:

Further by relaxing the integer constraint and replacing the rank constraint with the nuclear norm regularization, which is a standard technique for low-rank matrix completion, we get the desired convex program:

The clustering algorithm based on the above convex program is given in Algorithm 2

Denote the SVD of the block rating matrix $B$ by $B=U_{B}\Sigma_{B}V_{B}^{\top}$. The next lemma shows that

and thus $\mu$ is upper bounded by $\sqrt{r}$.

The following theorem provides a sufficient condition under which Algorithm 2 exactly recovers the rating matrix and thus the row and column clusters as well.

If $n(1-\epsilon)\geq C^{\prime}\log^{2}n$ for some constant $C^{\prime}$, and

where $C$ is a constant and $\mu$ is the incoherence parameter for $R$, then a.a.s. the rating matrix $R$ is the unique maximizer to the convex program (2) with $\lambda=3\sqrt{(1-\epsilon)n}$.

Note that Algorithm 2 is easy to implement as $\lambda$ only depends on the erasure probability $\epsilon$, which can be reliably estimated from $\widehat{R}$. Moreover, the particular choice of $\lambda$ in the theorem is just to simplify notations. It is straightforward to generalize our proof to show that the above theorem holds with $\lambda=C_{1}\sqrt{(1-\epsilon)n}$ for any constant $C_{1}\geq 3$.

Using Lemma 1, we immediately conclude from the above theorem that the convex program succeeds when $m>Cnr^{2}$ for some constant $C$. However, based on extensive simulation in Fig 2, we conjecture that the following result is true.

By (3), Conjecture 1 is equivalent to $\|U_{B}V_{B}^{\top}\|_{\infty}=\Theta(\sqrt{\frac{\log r}{r}})$. For a fixed $r$, we simulate $1000$ independent trials of $B$, pick the largest value of $\|U_{B}V_{B}^{\top}\|_{\infty}$, scale it by dividing $\sqrt{\log r/r}$, and get the plot in Fig 2.

Assuming this conjecture holds, Theorem 4 implies that

for some constant $C$ is sufficient to recover the rating matrix, which is better than the previous condition by a factor of $r$. We do not have a proof for the conjecture at this time.

Comparison to previous work In the noiseless setting with $p=0$, the nuclear norm minimization approach can be directly applied to recover data matrix and further recover the row and column clusters. It is shown in that the nuclear norm minimization approach exactly recovers the matrix with high probability if $m=\Omega(\mu^{2}nr\log^{2}n)$. The performance guarantee for Algorithm 2 given in (4) is better by at least a factor of $\log n$. Theorem 3 shows that the combinatorial method exactly recovers the row and column clusters if $m=\Omega(nr^{1/2}\log^{1/2}n)$, which is substantially better than the two previous conditions by at least a factor of $r^{1/2}$. This suggests that a large performance gap might exist between exponential-time algorithms and polynomial-time algorithms. Similar performance gap due to computational complexity constraint has also been observed in other inference problems like Sparse PCA and sparse submatrix detection .

In the low noise setting with $p$ restricting to be a small constant, the low-rank plus sparse matrix decomposition approach can be applied to exactly recover data matrix and further recover the row and column clusters. It is shown in that a weighted nuclear norm and $l_{1}$ norm minimization succeeds with high probability if $m=\Omega(\rho_{r}\mu^{2}nr\log^{6}n)$ and $p\leq\rho_{s}$ for two constants $\rho_{r}$ and $\rho_{s}$. The performance guarantee for Algorithm 2 given in (4) is better by several $\log n$ factors and we allow the fraction of noisy entries $p$ to be any constant less than $1/2$. Moreover, our proof turns out to be much simpler.

VII Spectral Method

In this section, we study a polynomial-time clustering algorithm based on the spectral projection of the observed rating matrix $\widehat{R}$. The description is given in Algorithm 3.

The following theorem shows that the spectral method exactly recovers the user and movie clusters under a condition stronger than (4). In particular, we show that Step 3 exactly recovers the block rating matrix $B$ and Step 4 cleans up clustering errors made in Step 2.

for a positive constant $C$, then Algorithm 3 with $\tau=12(1-\epsilon)^{1/2}r\log n$ a.a.s. exactly recovers user and movie clusters, and the rating matrix $R$.

Algorithm 3 is also easy to implement as $\tau$ only depends on parameters $\epsilon$ and $r$. As mentioned before, the erasure probability $\epsilon$ can be reliably estimated from $\widehat{R}$ using empirical statistics. The number of clusters $r$ can be reliably estimated by searching for the largest eigen-gap in the spectrum of $\widehat{R}$ (See Algorithm 2 and Theorem 3 in for justification). We further note that the threshold $\tau$ used in the theorem can be replaced by $C_{1}(1-\epsilon)^{1/2}r\log n$ for any constant $C_{1}\geq 12$.

Comparison to previous work Variants of spectral method are widely used for clustering nodes in a graph. Step 2 of Algorithm 3 for approximate clustering has been previously proposed and it is analyzed in . In , an adaptation of Step 1 is shown to exactly recover a fixed number of clusters under the planted partition model. More recently, proves an upper bound on the number of nodes “mis-clustered” by spectral method under the stochastic block model with a growing number of clusters.

Compared to previous work, the main novelty of Algorithm 3 is Steps 1, 3, and 4 which allow for exact cluster recovery even with a growing number of clusters. To our knowledge, Theorem 5 provides the first theoretical result on spectral method for exact cluster recovery with a growing number of clusters.

VIII Proofs

VIII-B Proof of Theorem 2

Consider a genie-aided scenario where the set of flipped entries is revealed as side information, which is equivalent to saying that we are in the noiseless setting with $p=0$. Then the true partition corresponding to the true user cluster structure has zero disagreement. Suppose users $1,3,\ldots,2K-1$ are in true cluster $1$ and users $2,4,\ldots,2K$ are in true cluster $2$. We construct a new partition different from the true partition by swapping user $1$ and user $2$. In particular, under the new partition, user $1$ forms a new cluster $\widehat{C}_{2}$ with users $2i,i=2,\ldots,K$, user $2$ forms a new cluster $\widehat{C}_{1}$ with users $2i-1,i=2,\ldots,K$. It suffices to show that for $k=1,2$, any two users in $\widehat{C}_{k}$ has zero disagreement with probability at least $3/4$, in which case the new partition has zero agreement and Algorithm 1 cannot distinguish between the true partition and the new one.

For $k=1,2$, we lower bound the probability that any two users in $\widehat{C}_{k}$ has zero disagreement.

By union bound, the probability that for $k=1,2$, any two users in $\widehat{C}_{k}$ has zero disagreement is at least $3/4$.

VIII-C Proof of Theorem 3

Consider a compatibility graph with $n$ vertices representing users. Two vertices $i,i^{\prime}$ are connected if users $i,i^{\prime}$ have zero disagreement, i.e., $D_{ii^{\prime}}=0$. In the noiseless setting, each user cluster forms a clique of size $K$ in the compatibility graph. We call a clique of size $K$ in the compatibility graph a bad clique if it is formed by users from more than one cluster. Then to prove the theorem, it suffices to show that there is no bad clique a.a.s. Since the probability that bad cliques exist increases in $\epsilon$, without loss of generality, we assume $K(1-\epsilon)<1$.

Recall that $B_{kl}$ is $+1$ or $-1$ with equal probability. Define $S_{k}=\{l:B_{kl}=+1\}$ for $k=1,\ldots,r$. As $r\to\infty$, by Chernoff bound, we get that a.a.s., for any $k_{1}\neq k_{2}$

Assume this condition holds throughout the proof.

Let $p_{n_{1}\dots n_{t}}$ be the probability that $K$ users, out of which $n_{k}$ are from cluster $k$, form a bad clique. Because $\{E_{j}\}$ are independent and there are $K$ movies in each movie cluster,

for some constant $C_{1}$. For a large enough constant $C$ in the assumption regarding $m$ in the statement of the theorem and using the fact that $m=n^{2}(1-\epsilon)$, we have

Below we show that the probability of bad cliques existing goes to zero. By the Markov inequality and linearity of expectation,

where the first term in last inequality corresponds to the case of $n_{\max}\leq K/2$ and the second term corresponds to the case of $n_{\max}>K/2$. They follows from (7), (8), (9) and the fact that $\binom{K}{n_{k}}\leq\min\{K^{n_{k}},K^{K-n_{k}}\}$.

VIII-D Proof of Theorem 4

We first introduce some notations. Let $u_{C,k}$ be the normalized characteristic vector of user cluster $k$, i.e., $u_{C,k}(i)=1/\sqrt{K}$ if user $i$ is in cluster $k$ and $u_{C,k}(i)=0$ otherwise. Thus, $||u_{C,k}||_{2}=1$. Let $U_{C}=[u_{C,1},\dots,u_{C,r}]$. Similarly, let $v_{C,l}$ be the normalized characteristic vector of movie cluster $l$ and $V_{C}=[v_{C,1},\dots,v_{C,r}]$. It is not hard to see that the rating matrix $R$ can be written as $R=KU_{C}BV_{C}^{\top}$. Denote the SVD of the block rating matrix $B$ by $B=U_{B}\Sigma_{B}V_{B}^{\top}$, then the SVD of $R$ is simply $R=UK\Sigma_{B}V^{\top}$, where $U=U_{C}U_{B}$ and $V=V_{C}V_{B}$. When $r\to\infty$, $B$ has full rank almost surely . We will assume $B$ is full rank in the following proofs, which implies that $U_{B}U_{B}^{\top}=I$ and $V_{B}V_{B}^{\top}=I$. Note that $UU^{\top}=U_{C}U_{C}^{\top},VV^{\top}=V_{C}V_{C}^{\top}$ and $UV^{\top}=U_{C}U_{B}V_{B}^{\top}V_{C}^{\top}$.

Assume user $i$ is from user cluster $k$ and movie $j$ is in movie cluster $l$, then

where the inequality follows from the Cauchy-Schwartz inequality. By definition $\mu\leq\sqrt{r}$. ∎

The following corollary applies Theorem 1.4 in to bound the spectral norm $\|\widehat{R}-\bar{R}\|$.

If $\sigma^{2}\geq C^{\prime}\log^{4}n/n$ for a constant $C^{\prime}$, then conditioned on $R$,

We adopt the trick called dilations . In particular, define $A$ as

where the last inequality holds when $n$ is sufficiently large. ∎

For any feasible $Y$ that $Y\neq R$, we have to show that $\Delta(Y)=\langle\widehat{R},R\rangle-\lambda||R||_{*}-(\langle\widehat{R},Y\rangle-\lambda||Y||_{*})>0.$ Rewrite $\Delta(Y)$ as

where the last inequality follows from definition of the incoherence parameter $\mu$. Below we bound the term $||\mathcal{P}_{T}(W)||_{\infty}$. From the definition of $\mathcal{P}_{T}$ and the fact that $U_{B}U_{B}^{\top}=I$ and $V_{B}V_{B}^{\top}=I$,

We first bound $||U_{C}U_{C}^{\top}W||_{\infty}$. To bound the term $(U_{C}U_{C}^{\top}W)_{ij}$, assume user $i$ belongs to user cluster $k$ and let $\mathcal{C}_{k}$ be the set of users in user cluster $k$. Recall that $u_{C,k}$ is the normalized characteristic vector of user cluster $k$. Then

which is the average of $K$ independent random variables. By Bernstein’s inequality (stated in the supplementary material), with probability at least $1-n^{-3}$,

Then $||U_{C}U_{C}^{\top}W||_{\infty}\leq\frac{1}{K}\left(\sqrt{\frac{2}{3r}\log n}+\frac{2\log n}{\lambda}\right)$ with probability at least $1-n^{-1}$. Similarly we bound $||WV_{C}V_{C}^{\top}||_{\infty}$ and $||U_{C}U_{C}^{\top}WV_{C}V_{C}^{\top}||_{\infty}$. Therefore, with probability at least $1-3n^{-1}$,

for some constants $C_{1}$ and $C_{2}$, where the second inequality follows from assumption (4). Substituting (17) into (16) and by assumption (4) again, we conclude that $\Delta(Y)>0$ a.a.s. ∎

VIII-E Proof of Theorem 5

The proof is divided into three parts. Recall that $x_{i}$ denotes the $i$-th row of $Pr(\widehat{R}^{(1)})$. We first show that, for most users, ${x}_{i}$ is close to the expected value conditioned on $R$. Then we show that the clusters output by Step $2$ are close to the true clusters. Finally, we show that Step $3$ exactly recovers the block rating matrix $B$ and Step $4$ exactly recovers clusters.

If $\sigma^{2}\geq C^{\prime}\log^{4}n/n$ for a constant $C^{\prime}$, then a.a.s., $|\mathcal{I}^{c}|\leq K\log^{-2}n$ and $|\mathcal{J}^{c}|\leq K\log^{-2}n$.

Let $(\sigma^{(1)})^{2}=\frac{1}{2}(1-\epsilon)$. By Corollary 1, $\|\widehat{R}^{(1)}-\bar{R}^{(1)}\|\leq 3\sigma^{(1)}\sqrt{n}$. Note that

where the second inequality follows from the definition of $P_{r}(\widehat{R}^{(1)})$ and the fact that $\bar{R}$ has rank $r$. Since both $P_{r}(\widehat{R}^{(1)})$ and $\bar{R}$ have rank $r$, the matrix $P_{r}(\widehat{R}^{(1)})-\bar{R}$ has rank at most $2r$, which implies that

As $\sum_{i=1}^{n}\|x_{i}-\bar{x}_{i}\|^{2}_{2}=\|P_{r}(\widehat{R}^{(1)})-\bar{R}\|^{2}_{F}$, we conclude that there are at most $K\log^{-2}n$ users with

Similarly we can prove the result for movies. ∎

The following proposition upper bounds the set difference between the estimated clusters and the true clusters by $K\log^{-2}n$. Let $C^{*}_{1},\ldots,C^{*}_{r}$ be the true user clusters and $\Delta$ denote the set difference.

Assume the assumption of Theorem 5 holds. Step 2 of Algorithm 3 outputs $\{\widehat{C}_{k}\}_{k=1}^{r}$ and $\{\widehat{D}_{l}\}_{1=1}^{r}$ such that, up to a permutation of cluster indices, a.a.s., $\widehat{C}_{k}\Delta C^{*}_{k}\subset\mathcal{I}^{c}$ and $\widehat{D}_{l}\Delta D^{*}_{l}\subset\mathcal{J}^{c}$ for all $k,l$. It follows that for all $k,l$,

It suffices to prove the conclusion for the user clusters. Consider two good users $i,i^{\prime}\in\mathcal{I}$. If they are from the same cluster, we have $\bar{x}_{i}=\bar{x}_{i^{\prime}}$ and

where the last inequality follows from Lemma 2. If they are from different clusters, by (6), we have a.a.s.

where $R_{i}$ denotes the $i$-th row of $R$. Thus,

where the last inequality follows from the assumption (5). Therefore, in the clustering procedure of Step $2$, if we choose a good initial user at some iteration, the corresponding estimated cluster will contain all the good users from the same cluster as the initial user and no good user from other clusters. It is not hard to see that the probability of the event that we choose a good initial user in every iteration is lower bounded by

Assume the above event holds. Under proper permutation, the initial good user in the $k$-th iteration is from cluster $C_{k}^{*}$ for all $k$. By the above argument, the set difference $\widehat{C}_{k}\Delta C^{*}_{k}\subset\mathcal{I}^{c}$. By Lemma 2, (18) follows. ∎

We first show that Step $3$ of Algorithm 3 exactly recovers the block rating matrix $B$. Let $V_{kl}$ denote the total vote that the true user cluster $k$ gives to the true movie cluster $l$, i.e.,

Then by definition of $\widehat{V}_{kl}$,

Without loss of generality, assume $B_{kl}=1$. By Bernstein inequality and assumption (5), $V_{kl}\geq\frac{1}{4}(1-\epsilon)(1-2p)K^{2}$ a.a.s. On the other hand, as $\Omega_{2}$ and $\widehat{R}^{(1)}$ are independent, $\Omega_{2}$ is independent from $\{\widehat{C}_{k}\}$ and $\{\widehat{D}_{l}\}$. It follows from (18) and the Chernoff bound that each term on the right hand side of (21) is upper bounded by $(1-\epsilon)K^{2}\log^{-2}n$ a.a.s. Hence, when assumption (5) holds for some large enough constant $C$, we have $\widehat{V}_{kl}>0$ thus $\widehat{B}_{kl}=B_{kl}$.

Next we prove that Step $4$ clusters the users and movies correctly. Without loss of generality, we only prove the correctness for users. Suppose user $i$ is from cluster $k$. Recall that $R_{i}$ denotes the $i$-th row of $R$. When $\widehat{B}=B$, we have $\mu_{kj}=R_{ij}$ for $j\in\mathcal{J}$ by definition and Proposition 1. Then

Similarly, for some user $i^{\prime}$ from cluster $k^{\prime}\neq k$,

and $\text{Var}[\langle\widehat{R}^{(2)}_{i},R_{i^{\prime}}\rangle]\leq n(\sigma^{(2)})^{2}$. Now by the Bernstein inequality and assumption (5), we have that conditioned on $R$, a.a.s. $\langle\widehat{R}^{(2)}_{i},R_{i}\rangle>7t/8$ and $\langle\widehat{R}^{(2)}_{i},R_{i^{\prime}}\rangle<5t/8$ for all $i\neq i^{\prime}$.

On the other hand, because $\mathcal{J}$ and $\Omega_{2}$ are independent, by the Chernoff bound, a.a.s. $\sum_{j\in\mathcal{J}^{c}}|\widehat{R}^{(2)}_{ij}|$ is upper bounded by $(1-\epsilon)K\log^{-2}n<t/16$ for all $i$, when assumption (5) holds for some large enough constant $C$.

Therefore, from (22) and (23), $\langle\widehat{R}^{(2)}_{i},\mu_{k}\rangle>\langle\widehat{R}^{(2)}_{i},\mu_{k^{\prime}}\rangle$ for all $k^{\prime}\neq k$. ∎

IX Numerical Experiments

In this section, we illustrate the performance of the convex method and the spectral method using synthetic data.

The convex program (2) can be formulated as a semidefinite program (SDP) and solved using a general purpose SDP solver. However this method does not scale well for our problem when the matrix dimension $n$ is large. Instead we propose a first-order algorithm motivated by the Singular Value Thresholding algorithm proposed in . Consider the following convex program which introduces an additional term $\frac{\tau}{2}||Y||_{F}^{2}$ for $\tau>0$,

When $\tau$ is small, the solution of (24) is close to that of (2). We solve the convex program (24) using the dual gradient descent method given in Algorithm 4. Let matrix $E$ denote the matrix of all ones. Let matrices $S\geq 0$ For two matrices $X$ and $X^{\prime}$, $X\geq X^{\prime}$ means that $X_{ij}\geq X^{\prime}_{ij}$ for all entries $(i,j)$. and $T\geq 0$ denote the Lagrangian multipliers for the constraints $Y\leq E$ and $Y\geq-E$, respectively. The Lagrangian function is given by

and the dual function is given by $\min_{Y}L(Y,S,T)$. The gradients of the dual function with respect to $S$ and $T$ are $Y-E$ and $-Y-E$, respectively.

Intuitively, the soft-thresholding operator $D$ shrinks the singular values of $X$ towards zero. Applying Theorem 2.1 in , we get

Thus, the first update equation in (25) finds the minimizer of the Lagrangian function at the current estimate of the dual variables. The next two update equations in (25) move the current estimate of the dual variables in the direction of the corresponding gradients of the dual function, and then project the new estimates to the set of matrices with nonnegative entries. The parameter $\delta>0$ is the step size.

We simulate Algorithm 4 on the synthetic data. Assume $K$ and $\epsilon$ take the form given by

Theorem 4 shows that the convex program (4) recovers the rating matrix exactly when $\alpha<\beta$, assuming Conjecture 1 holds.

IX-B Spectral Method

We simulate the spectral method given in Algorithm 3 on synthetic data. Assume $K$ and $\epsilon$ take the form of (26). Theorem 5 shows that the spectral method exactly recovers the clusters when $\alpha<\frac{1}{2}(\beta+1)$.

We generate the observed data matrix according to our model with $n=1000$ and $p=0.05$, and various choices of $\beta,\alpha\in(0,1)$. We apply Algorithm 3 with slight modifications. Firstly, we do not split the observation as in Step $1$ but use all the observations for the later steps, i.e., $\Omega_{1}=\Omega_{2}=\Omega$. Secondly, in Step $2$ we use the more robust $k$-means algorithm to cluster users and movies instead of the thresholding based clustering algorithm.

To calculate the number of mis-clustered users (movies), we need to consider all possible permutations of the cluster indices, which is computationally expensive for large $r$. Thus, the clustering error is instead measured by the fraction of misclassified pairs of users and movies. In particular, we say a pair of users (movies) misclassified if they are either from the same true cluster but assigned to two different clusters or from two different true clusters but assigned to the same cluster. We say the algorithm succeeds if the clustering error is less than $5\%$.

For each $\beta$, we run the algorithm for several values of $\alpha$ and record the largest $\alpha$ for which the algorithm succeeds. The result is depicted in Fig 4. The solid blue line represents $\alpha=\frac{1}{2}(\beta+1)$, which shows the performance guarantee of the spectral method given by Theorem 5. The dotted red line represents $\alpha=\beta$, which shows the performance guarantee of the convex method given by Theorem 4. We can see that our simulation results are slightly better than the theoretical performance guarantee.

X Concluding Remarks

This paper studies the problem of inferring hidden row and column clusters of binary matrices from a few noisy observations through theoretical analysis and numerical experiments. More extensive simulation results will be presented in a longer version of the paper. Several future directions are of interest. First, proving Conjecture 1 is important to fully understand the performance of the convex method. Second, a tight performance analysis of the ML estimation (1) is needed to achieve the lower bound. Third, it is interesting to extend our analysis to block rating matrices having real-valued entries.

References