High-Order Co-Clustering via Strictly Orthogonal and Symmetric L1-Norm Nonnegative Matrix Tri-Factorization

1 Introduction

Recent advancements of Internet and other computational technologies have brought data with much richer structures, which often convey useful information for building more effective learning models. For example, one can use different types of information to recommend movies to a user, such as the user's watching history, comments on watched movies, habits, job, age and social networks, to name a few. These different types of information are represented as different types of data that usually interrelate with each other by various means. For example, user preferences, comments and watching history can be correlated via weighted co-occurrence matrices - some movies are more favored by a user, if they appear more frequently in the user's watching history; while, if a user gave negative comments to a movie, this user is less likely to watch other movies of the same genre. Such highly heterogeneous data with intertype relationships are called multi-type relational data [Long et al., 2006]. The learning problems to deal with multitype relational data under the unsupervised setting is called as High-Order Co-Clustering (HOCC) [Wang et al., 2011a; Wang et al., 2011b], which has attracted more and more attention in recent years.

∗Corresponding author. This work was partially supported by NSF-IIS 1423591 and NSF-IIS 1652943.

Spectral clustering [Long et al., 2006] was first proposed to cluster multi-type relational data, though only inter-type information was utilized. Recently, Wang et al. [2011b] devised an indicator Nonnegative Matrix Tri-Factorization (NMTF) method to avoid soft clustering and improve the computational speed. But since this method searches the binary space that is small, it may not always find a good optimum. It was then further developed by performing optimization in a continuous space [Wang et al., 2011a]. However, the clustering performance was not improved, because the orthogonal constraints were not imposed on the factor matrices in the new objective, which, though, was very important to avoid degenerate solutions [Ding et al., 2006]. Besides the limitations mentioned above, all these previous studies routinely formulated their NMTF objectives as least-square error functions, which are notoriously known to be sensitive to outliers [Kong et al., 2011]. However, at the era of big data, noise and outliers are inevitable by nature due to the ever increasing data sizes. To tackle this problem, Liu and Wang [2015] improved these prior studies by using the ℓ1-norm distances. However, because the new objective in [Liu and Wang, 2015] computes the symmetric NMTF that involves the fourth-order matrix polynomials, the factor matrices in the solution were not orthogonally constrained due to the mathematical difficulty. Moreover, same as many existing NMTF methods that derive the solutions using Multiplicative Updating Algorithm (MUA) [Lee and Seung, 2001], this method also suffers from suboptimal solutions because it is easily trapped into local minimum [Lin, 2007].

To address all above difficulties, in this paper we propose a novel HOCC method via a strictly orthogonal and symmetric ℓ1-norm NMTF method. We derive the solution algorithm using the Alternating Direction Method of Multipliers (ADMM) [Boyd et al., 2011], such that the orthogonal constraints are strictly enforced on the factor matrices. Extensive experiments have been conducted with promising results that validate the effectiveness of our new method on clustering multi-type relational data.

2 High-Order Co-Clustering via Graph Regularized Symmetric ℓ1-Norm NMTF

Throughout this paper, we use A(ij) to denote the entry at the i-th row and j-th column of a matrix A. A K-type relational data set can be denoted as χ = {χ1, χ2, . . . , χK}, where χk = {xk1 , xk2 , . . . , xknk } represents the data of the kth type. Suppose that we are given a set of relationship matrices {Rkl ∈ ℜnk×nl }(1≤k≤K,1≤l≤K) between different types of data objects, then we assume Rkl = RT lk . Our goal is to simultaneously partition the data objects in χ1, χ2, . . . , χK into c1, c2, . . . , cK disjoint clusters respectively.

2.1 Co-clustering Two-Type Relational Data Using Symmetric ℓ1-Norm NMTF

The simplest multi-type relational data involve only two types of data objects, which widely appear in many real world applications, like words and documents in document analysis, users and items in collaborative filtering, etc. Simultaneously clustering two-type relational data is often called coclustering or bi-clustering. Ding et al. [2006] proposed to use NMTF to simultaneously cluster the rows and columns of an input nonnegative relationship matrix Rij by decomposing it into three nonnegative factor matrices, which minimizes:

where ‖ · ‖2F denotes the squared Frobenius norm of a matrix, G1 ∈ ℜn1×c1 + and G2 ∈ ℜn2×c2 + are the cluster indicator matrices for χ1 and χ2 respectively, and S12 ∈ ℜc1×c2 + absorbs the different scales of R12, G1 and G2. Simultaneous clustering of χ1 and χ2 is then achieved by solving Eq. (1). Because the rows of resulted Gk (k ∈ {1, 2}) (with normalization) can be interpreted as the posterior probability for clustering on χk, the cluster label of xki is obtained by:

As can be seen, the data points enter the objective function in Eq. (1) as squared residual errors. Thus, outlying data samples can easily dominate the objective function because of the squared errors. To solve the problem, Liu and Wang [2015] proposed to replace the traditional squared Frobenius norm into ℓ1-norm and optimize the following objective:

where ‖X‖1 is defined as ∑ni=1 ∑mj=1 |X(ij) | for a matrix X ∈ ℜn×m.

2.2 Simultaneously ClusteringMulti-Type Relational Data Using Inter-Type Relationships

A natural way to generalize the co-clustering objective in Eq. (1) to simultaneously cluster multi-type relational data is to solve the following optimization problem:

Although Eq. (4) generalizes Eq. (1) to deal with multitype relational data, it is not straightforward to solve Eq. (4) by generalizing existing NMTF algorithms. As a contribution of this paper, we will derive the solution to Eq. (4).

We first introduce the following useful lemma.

Lemma 1 The optimization problem in Eq. (3) can be equivalently solved by the following problem:

in which

where the superscripts denote the matrix sizes, and R12 = RT 21 , S21 = ST 12 . 0n1×n1 is a matrix with all zero entries of size n1 × n1.

Proof 1 Following the definitions of R, G and S, we can derive

where P = G1S12GT2 , which proves the lemma.

Based upon Lemma 1, we have the following theorem.

Theorem 1 It is equivalent to solve Eq. (4) and to solve

in which

where Rji = RT ij and Sij = ST ji .

The proof of Theorem 1 can be easily obtained by generalizing the proof of Lemma 1 to multi-type relational data.

Theorem 1 presents a general framework via ℓ1-norm symmetric NMTF (S-NMTF) for simultaneously cluster multitype relational data using the mutual relationship matrices.

2.3 Incorporating Intra-Type Information via Graph Regularization

The optimization objectives in Eq. (4) and Eq. (7) only use the inter-type relationships of a multi-type relational data set, whereas the intra-type information, though often available, is not used. We can incorporate the intra-type relationship information through Laplacian regularization. For a multitype relational data set, given the intra-type information in the form of the pairwise affinity matrices W1,W2, . . . ,WK respectively, we can use them as following:

where Lk = Dk -Wk is the Laplacian matrix with D(ii) k = ∑j W (ij) k . Because Lk is the discrete approximation of the Laplace-Beltrami operator on the underlying data manifold, the regularization term reflects the label smoothness of the two types of data points. The smoother the data labels are with respect to the underlying data manifolds, the smaller their values will be.

Using R, S and G defined in Eq. (8), denote

and L = D-W where D(ii) = ∑j W (ij) , it is easy to prove

Combining Theorem 1, we approach simultaneous clustering of multi-type relational data with manifold in Eq. (9) by solving the following optimization problem:

It can be verified that with constraint GT G = I, then GT(i)G(i) = I, ∀ 0 < i ≤ K, which means the orthogonality of G will make membership matrices in each relational data orthogonal. The benefits of orthogonality constraint on G are two-fold: 1) get unique solution; 2) avoid soft clustering.

3 Difficulty in Solving the Problemwith Traditional Method

Traditional methods to solve NMTF problem are based on MUA, where an auxiliary function is introduced to give an upper bound for the objective during every iteration. However, our proposed ℓ1-norm objective function, different from traditional squared Frobenius norm which has nice mathematical properties, is non-differentiable. As a result, we cannot apply the traditional method to solve Eq. (12).

In our previous work [Liu and Wang, 2015], we attempted to solve the ℓ1 norm by transforming it into a trace problem as:

where D is a diagonal matrix defined as:

By taking the derivative of J with respect to G and S, they get the solution following the traditional method. However, the derivation is not rigorous since D is a dependent variable on G and S, while they are suppose to be independent when we take derivatives.

As a contribution of our paper, we propose an algorithm which can give the optimal solution with rigorous mathematical guarantee. Moreover, a strictly orthogonal membership matrix will be given to avoid soft clustering while the objective loss is usually lower than MUA based methods.

4 The Solution Algorithm

Before giving our algorithm, we will first introduce the ADMM, which was proposed in [Bertsekas, 1996; Boyd et al., 2011] to solve convex optimization problems by breaking them into smaller pieces that are easier to handle. Specifically, given the following objective with the equality constraint:

Algorithm 1 solves the problem by decoupling it into subproblems and optimizing each variable while fixing others [Bertsekas, 1996; Boyd et al., 2011], where y is the Lagrangian multiplier to the constraint h. It is worth noting that Algorithm 1 was proved to converge Q-linearly to the optimal solution [Bertsekas, 1996].

Algorithm 1: The ADMM algorithm.

• 1 Set 1 < ρ < 2 and initialize µ > 0 and y; 2 while not converge do 3 1. Update x by solving xk+1 = arg minx(f (x) + µ2 ‖h(x, zk ) + yk µ ‖2 ); 4 2. Update z by solving zk+1 = arg minz(g(z) + µ2 ‖h(xk+1 , z) + yk µ ‖2 ); 5 3. Update y by yk+1 = yk +µh(xk+1 , zk+1 ); 6 4. Update µ by µ = ρµ. 7 end

Huang et al. [2014] used the ADMM to solve the manifold regularized NMF problem. Following the same idea, we derive the solution algorithm to our objective.

Because there are two constraints (nonnegative and orthogonal) on G in our objective in Eq. (12), which is difficult to optimize, we introduce three auxiliary variables:

F = G; P = F;H = G. In addition, considering that the ℓ1-norm objective is non-differentiable, we introduce one additional variable of E = R- FSGT to solve the following new objective function:

The constraints in Eq. (12) now become: GT G = I, F T F = I , H ≥ 0, P ≥ 0, where each variable has only one constraint which is easier to be optimized.

Eq. (16) can be reduced to several manageable subproblems, where each subproblem yields a closed-form solution as following. We repeat the optimizing procedure until convergence, which is summarized in Algorithm 2.

Step 1. Solving S when fixing other variables:

where we write ˜R = R-E - 1µ Λ for brevity. Taking derivative of Eq. (17) with respect to S and setting it as 0, we can easily obtain the solution to this optimization problem as S = (F T F )-1 F T ˜RG (GT G)-1 .

Step 2. Solving E when fixing other variables:

where we write Z = R - FSGT - 1µ λ for brevity.

The optimization problem in Eq. (22) can be decoupled to solve the following problem for every entry of E:

Taking derivative of Eq. (19) with respect to eij and setting it as 0, we can solve Eq. (19) as follows:

Step 3. Solving G when fixing other variables:

By mathematical derivation, the problem in Eq. (21) can be rewritten as follows:

According to [Wang et al., 2013], the problem in Eq. (22) can be solved by computing the SVD of M: if svd (M) = UAV T , the solution of Eq. (22) is given by UV T .

Step 4. Solving F when fixing other variables, similar to optimize G, the subproblem is equal to:

where we write M = (R-E - 1µ Λ) GST + (G- 1µ Σ+P + 1µ ∆- λµ LG) for brevity. By computing the SVD of M: svd (M) = UAV T , the solution of F is given by UV T .

Step 5. Solving H when fixing other variables:

and the solution is H = max (G - 1µ Ω, 0).

Step 6. Solving P when fixing other variables, similar to optimize H, P = max (F - 1µ ∆, 0).

Step 7. Update Lagrangian Multipliers Λ, Σ, ∆, Ω.

Step 8. Update µ.

Algorithm 2: The solution algorithm to our objective.

Data: Multi-relational data: {R1 ,R2 , ...,RnK },

Number of Clusters K and set 1 ≤ ρ ≤ 2.

Result: Factor matrices: G.

1 1. Construct R and W with {R1 , R2 , ..., RnK };

2 2. Initialize G, Λ,Σ,∆,Ω;

3 while not converge do

4 3. Optimize G, S, F, E, P, H as Eq. (17-24);

5 4. Update Lagrangian Multipliers Λ,Σ,∆,Ω;

6 5. Update µ = ρµ.

7 end

8 6. Get the clustering result fromG.

We can get the corresponding solution by changing ℓ1 norm to traditional squared Frobenius norm. The only difference is to update E. When it is squared Frobenius norm, we can solve the subproblem at every entry of E as following:

Taking derivative of Eq. (25) with respect to eij and setting it as 0, we can solve Eq. (25) as follows:

Due to space limit, we skipped some details in above derivations, which will be supplied in our extended journal version of this paper.

Figure 2: Left: Objective comparison between our proposed method and MUA based algorithm. Right: Zoom in iteration 25 to 29, showing ours can eventually get a lower local minimum.

5 Experiments

Although multi-relational data can help to improve clustering, there is no such existing real-world data set which can be directly used for experiments. For our experiments, we created a new multi-relational information from ACL-IMDB data set.

5.1 Multi-Relational Data Set and Its Collection

In the ACL-IMDB [Maas et al., 2011] data set, there are 25000 reviews for highly popular movies, in which positive and negative comments come up with one half (12500) each. In our experiments, we build three relational matrices: commentword, user-comment and user-word respectively.

Comment-word matrix. Since the comments are scored from 1 to 10 points (positive from 7 and 10 while negative from 1 to 4) with detailed content (word as the unit or feature), we can build the comment-word matrix along with its corresponding labels (ground truth), which will be used for sentiment analysis (comment clustering) of our proposed algorithm. We first rule out the stopwords (stopwords, such as the, here, etc., are words which are common across documents and contribute little to the content of the document), and then get the top 700 most commonly used words such as cool, great, etc. as the features.

User-comment matrix. Comments from the same person by using the same words are more likely to belong to the same cluster (positive or negative). Since each comment in the data set has its corresponding URL, of which we can make use to identify the author. Therefore we can build the user-comment matrix.

User-word matrix. Given the comment-word and usercomment matrix respectively, the user-word matrix can be obtained by multiplying the two matrices described above.

The three matrices create the multi-relational data, which are expected help improve the clustering performance. It is worth noting that the generated comment-word matrix is noisy due to some inevitable post typos. To make our experiments more convincing, we add some noise to the three relationship matrices with a ratio up to 25% (up to 20% in amplitude). By randomly choosing 500 authors from the top 1500 who post most comments, we can generate sub-data sets to conduct our experiments.

5.2 Convergence of the Algorithm

Because our solution algorithm employs an iterative method, we first study its convergence property. From Fig. 1 we can see that, though there is an increase in the first update, it decreases sharply later. In our experiment, we initialize µ as 0.02 and keep updating it by µ = ρµ during the iterations to accelerate the decrease. we find that if µ is fixed to be very small such as 0.001, the decrease will be slow; while if it is large such as 0.5, the decrease will be insignificant, trapped into a poor local minimum.

In addition, we compare clustering method using the squared Frobenius norm loss function (solved by the ADMM method by updating E using Eq. (26)) against the same objective solved by the MUA based method in [Wang et al., 2011a]. As demonstrated in Fig. 2, our proposed method can get a lower minimum that may lead to a better clustering result.

5.3 Time Consumption

We also compare the time consumption between our proposed ℓ1-norm algorithm against that in [Liu and Wang, 2015] (even though the proof is not rigorous and the loss does not decrease monotonically as expected) . We find that our algorithm is faster during each iteration, and it only takes 1/3 to 1/2 time as the counterpart does. We experiment with different sizes of matrices, and find our method is always faster as Fig. 3 shows.

5.4 Parameter Study and Robustness of ℓ1-norm Objective Function

We incorporate the intra-type data through manifold regularization. We study the influence of regularization item to clustering accuracy by changing λ from 0.1 to 2 with a step size of 0.1, and find that ℓ1-norm objective function remains robust with parameter changes, and its accuracy is constantly higher than the counterpart as Fig. 4 demonstrates. In our new method, we strictly enforce the orthogonal constraints onto the clustering membership matrix G, which plays an important role in getting rid of soft clustering (weights to different clusters are close) may lead to more distinguishable in clustering result. As we can see from Fig. 5, our solution of G satisfies the orthogonality constraint. In contrast, methods based on traditional MUA all fail to yield a strict orthogonal clustering matrix, which results in many soft clustering cases.

5.5 Clustering Result

We compare our proposed methods with some typical algorithms w.r.t. Accuracy (Acc), Normalized Mutual Information (NMI) [Xu et al., 2003] and Adjusted Rand Index [Yeung and Ruzzo, 2001]:

K-means: We conduct K-means on comment-word and user-comment matrix respectively and take the best result.

ONMTF: Ding et al. [2006] solves symmetric matrix trifactorization with orthogonal constraint while ignoring the intra-type data.

GNMF: By making use of graph, Cai et al. [2011] imposes the graph regularization to improve clustering.

HOCC: We only utilize two types of relational data (comment-word matrix) as Eq. (9) (where K = 2) for clustering, while ignoring other relational type data.

ℓ2-MUA: Wang et al. [2011a] solves ℓ2-norm HOCC with traditional MUA based method.

ℓ1-MUA: Liu et al. [2015] solves ℓ1-norm HOCC which is supposed to be robust with noise and outliers.

ℓ2-ADMM: Our proposed ℓ2-norm objective algorithm which is sensitive to outliers.

ℓ1-ADMM: Our proposed ℓ1-norm objective algorithm.

In Table 1, we change the number of users to generate symmetric matrices with different dimensions, and run the algorithms for 20 times and take the average for comparison. We also add some artificial noise to the data as mentioned in Section 5.1. In Table 2, the top section are clustering accuracies and NMI with different matrix sizes without noise while the middle and bottom sections show the clustering results before and after the noise added. We see that, in most cases, our method can get the best clustering result.

Due to space limit, we cannot report the experimental results on two-type relational data, which will be supplied in our extended journal version of this paper.

6 Conclusion

In this paper, we explored the high-order co-clustering problem and solved with the symmetric nonnegative matrix trifactorization method. Different from traditional squared error loss, we proposed a non-differentiable ℓ1-norm objective to stay robust with both sample and feature outliers. Instead of following traditional MUA method, we propose an ADMMbased algorithm and compare our method with other methods on real-world multi-relational data set we collect. All experiment results demonstrate that our methods are superior in terms of less time consumption, lower local minimum, higher clustering accuracy, NMI and Adjusted Rand Index, as well as strictly orthogonal membership matrix to get rid of soft clustering. Our method is flexible and could also be a new framework to solve non-smooth, non-convex problems with promising preliminary experimental results.