Data Clustering: Principal Components, Hopfield and Self-Aggregation Networks

1 Introduction

Principal component analysis (PCA) is widely adopted as an effective unsupervised dimension reduction method. PCA is extended in many different directions [Hastie and Stuetzle, 1989; Kramer, 1991; Lee and Seung, 1999; Scholkopf et al., 1998; Collins et al., 2001].

The main justification is that PCA uses singular value de­ composition (SVD) which is the best low rank approxima­ tion in £2 norm to original data due to Eckart-Young theo­ rem. However, this results alone is inadequate to explain the effectiveness of PCA. Here, we provide a new derivation of PCA based on optimizing suitable clustering objective func­ tions and show that principal components are actually cluster indicator vectors in clustering.

Hopfield network[Hopfield, 1982] provide a convenient framework of our study. In particular, the self-aggregation network proposed in this work uses Hebb rule to encode pat­ tern vectors. One feature of Hopfield associative-memory networks is that it can be adopted to solve hard combinato­ rial problems[Haykin, 1998 2nd ed].

Another thread of this work is the spectral graph partition­ ing [Fiedler, 1973; Pothen et al., 1990; Hagen and Kahng, 1992; Shi and Malik, 2000; Ding et al., 2001a; 2001b; Ng et at., 2001; Meila and Shi, 2001], which uses Laplacian matrix of a graph. This arises naturally for balancing the clus­ ters (see §2.4). Our approach differs from others mainly in well-motivated clustering objective functions. A further point is the recognition that spectral graph clustering is embedded in the scaled PCA thus leading to self-aggregation networks.

*Work supported by Department of Energy (Office of Science, through a LBNL LDRD) under contract DE-AC03-76SF00098.

This paper combines above three threads and develop a coherent framework for clustering. We begin with two-way clustering in §2 and generalize to multi-way clustering in §3.

2 1\vo-way clustering

We start by formulate two-way clustering as a Hopfield net­ work and derive PCA as cluster indicator vectors.

2.1 Hopfield Network and PCA

Given n data points and properly defined similarity or associ­ ation w;1 = w1; 2: 0 between points i, j, we form a network (a weighted graph) G with w; 1 as the connection between nodes i,j. We wish to partition it into two clusters c1,c2. The result of clustering can be represented by an indicator vector q, where

Consider the clustering objective,

where s(c1,c2) = LiECJ,jEC 2 w;1 is the overlap between c1, c2; s( ch ci) = Li,JEC 1 w;j is the within-cluster similar­ ity of CJ, and analogously for s( c2, c2).

Now we propose a min-max clustering principle: data points are grouped into clusters such that the overlaps( c1, c2) between different clusters are minimized while within-cluster similarities s(c~, c1), s(c2, c2) are maximized

These conditions can be simultaneously satisfied by maximiz­ ing the energy (objective function) J1 . Using Hopfield model [Hopfield, 1982], the solution is obtained by the update rule

or in vector form q(t+l) = sgn[Wq{tl]. If one relaxes q(i) from discrete indicators to continuous values in ( -1, 1 ), the solution q is given by

Since the matrix entries in Ware non-negative, the first prin­ cipal eigenvector q 1 has all positive (or all negative) entries. Thus the desired solution is q2 ·

2.2 Hopfield network for bipartite graph

Here we extend the Hopfield networks for clustering bipar­ tite graph. An example of bipartite graph is a m x n term­ document association matrix B = (bij ), where each row rep­ resents a word, each column represents a document, and b;j the counts of co-occurrence of row r; and column Cj. We show that the solution for clustering indicators is precisely the Latent Semantic Indexing[Deerwester et at., 1990].

We wish to partition the r-type nodes of R into two parts R1, R2 and simultaneously partition the c-type nodes of C into two parts c1, c2, based on the clustering principle of mini­ mizing between-cluster association and maximizing within­ cluster association. We use indicator vector f to determine how to split R into R1, R2: f( i) = 1, or -1 depending on r; E R1 orr; E R2. We use g to determine how to split C into c1, c2: g(i) = 1, or -1 depending on c; E c1 or c; E c2.

(For presentation purpose, we index the nodes such that nodes within same cluster are indexed contiguously. The clustering algorithms presented are independent to this as­ sumption. Bold face lower case letters are vectors. Matri­ ces are denoted by upper case letters.) Thus we may write

It is convenient to convert the bipartite graph into an undi­ rected graph. We follow standard procedure and combine the two types nodes to one by setting

This induces an undirected graph G, whose adjacency matrix is the symmetric weight matrix W.

Consider the following objective function,

where s(R1,c2) = L::r·ER c·EC b;j, the sum of association " 1' 1 2 between R1 and c2, and s(R" c2) is the sum of association between R1 and c2. s(R1, c2) and s(R2, c1) should be mini­ mized. s(R1, c1) is the sum of association within cluster 1 (see Fig.1), and s(R2 , c2 ) is the sum of association within cluster 2. s(R1, c1) and s(R2, c2) should be maximized. These condi­ tions are simultaneously satisfied by maximizing J2 .

We can write down the Hop field network update rule for q. If one relaxes q( i) from discrete indicators to continuous val­ ues, the solution q satisfies Eq.(4). Now utilizing the explicit structures of W and q, we have

The solutions to this equation are the singular value decom­ position (SVD) of B: { fi} are left singular vectors and {gi} are right singular vectors of the SVD of B,

Again, since the matrix entries in W are non-negative, the first principal components u 1 , v 1 have entries of same sign; thus the desired solutions are u 2 , v2 . Note that this is also the SVD employed in LSI. We summarize these results in

Theorem 1. Principal components are solutions for clus­ tering indicators for clustering undirect graphs and bipartite graphs under appropriate Hopfield network models.

2.3 Principal component clustering

In the objective J1, we may explicitly enforce a balance of clusters c1, c2. From Eq.1, we need to minimize I I:;; q(i)l. Thus we consider the clustering objective,

where J-L > 0 is a parameter and e = (1, · · ·, 1)T. We adjust J-L to control the levele of balance between c1 , c2 . The principal eigenvector q 1 of W = W - J-Leer is the desired indicator vector.

What is the reasonable choice of J-L? A natural choice is to set J-L = w,./ n 2 , w.. = L:;ij W;j, the average of W;j. With this choice, the modified weight matrix satisfies the sum-to-zero condition:

This condition can be further refined by centering each col­ umn and each rows, such that

where

(here column and row sums, wi. = L:;j W;j, W.j = I:;; W;j and w. . are standard notations in statistics.) Now the desired cluster indicator vector is the principal eigenvector of

The full~entered W has a useful property that all eigen­ vectors of W with nonzero eigenvalues have the sum-to-zero property: I:;; q(i) = 0. This is because (1) <io = e is an eigenvector of W with >.0 = 0; (2) all other eigenvectors are orthogonal to <i.o, i.e, qr e = 0.

The sum-to-zero condition I:;; q( i) = 0 does not necessar­ ily imply that the sizes of the two cluster, n 1 , n2, should be equal. In fact, the cluster indicator vector should be refined to

(in this paper, all vectors are implicitly normalized to 1 using £2 norm). Correspondingly, the clustering objective becomes

Clearly, the first two terms represent the average within­ cluster similarities which are maximized, and the 3rd term represent average between-cluster similarities which are min­ imized. The factor n1n2 encourages cluster balance, since n1 + n2 = n, and max(n1n2) is reached when n1 = n2 = n/2. We summarize these results in

Theorem 2. Solution for clustering objective 11 is given by q1; and solution for clustering objective 11 is q1.

We call these clustering schemes with objective functions 11 , 11 , 11 as principal component clustering methods because of the use of principal components.

All above discussions apply to bipartite graphs. For exam­ ple, 11 becomes

where nR 1 , nR 2 are sizes of row clusters (nR 1 + nR 2 = total row size); and nc 1 , nc 2 are sizes of column clusters (nc 1 + nc 2 =total column size). Again, the first two terms represent the average within-cluster similarities which are maximized, and the last two terms represent average between-cluster sim­ ilarities which are minimized.

We note that ] 1 is partially motived by the classic scal­ ing (also called principal coordinate analysis) in statistics [Borg and Gronen, 1997]. Suppose we are given pairwise distances diJ for a dataset, define the pairwise similarity as Wij = - ~dTJ, follow through the same centering procedures and solve for the eigenvectors of W. The classic scaling the­ orem proves that if dij are the Euclidean distances, the first K eigenvectors of W will recover the original coordinates. This theorem justifies the use of the K principal eigenvec­ tors as the coordinates in multidimensional scaling. However, our clustering approach does not emphasize the recovery of original coordinates, and the objective function ] 1 is well­ motivated for clustering only.

2.4 Cluster balance and MinMaxCut

Principal component clustering of §2.3 is motivated for bal­ ancing clusters by adding a panelity term to 11 . Here we balance clusters through a divisive approach.

Maximizing 11 of Eq.2 is equivalent to 1

Although s( c1, ci) + s( c2, c2) are maximized, it often oc­ curs that one cluster is much larger than another, s( c1, c1) » s( c2, c2) or vice versa.

1 Since exp(x) is monotonic, maxJ1 ~ maxexp(J1), which is minexp [s(c1 ,c2)]/ exp [s(CJ ,cJ) + s(c2,c2)] . This is approx­ imately Eq.l7.

To overcome this problem, we seek to prevent either s( c1, c1) or s( c2, c2) become very small. We optimize[Ding etal.,2001b]

One can show that

subject to qTW e = qT De = 0, where D = ( wi.) is a diago­ nal matrix. We relax q( i) from discrete indicators to real val­ ues in ( -1, 1). The solution of q for minimizing the Rayleigh quotient of Eq.(19) is given by

which can be written as

Again, the desired solutions is the eigenvector associated with the second largest eigenvalue. Note that by comparing Eq.(21) with Eq.(4), the net effect for cluster balance is di­ agonal scaling. This diagonal scaling, however, leads to an important feature of self-aggregation (see §3.2).

2.5 Cluster balance for bipartite graphs

For balanced clustering of bipartite graphs, object function of Eq.(7) becomes

Using the representation Wand q ofEq.6, and similar deriva­ tion [Ding, 2003; Zha et al., 2001], the solution for optimiza­ tion of 14 is also given by Eq.21. Let Dr = diag(Be) and De = diag(BT e), we have

Substituting into Eq.(21), we have

where

The solutions to Eq.(24) are SVD of B: B = L:k uk>-.kv"k; The clustering indicator vectors are fk = D:; 112 uk for row objects and gk = D-; 112 vk for column objects. Note that Eq.(24) is identical Eq.(8), with the corresponding relation-

clustering balancing using MinMaxCut objective of Eq.(22) over the simple objective Eq.(7) is the scaling of the associ­ ation matrix B in Eq.(25). However, with this scaling, the self-aggregation property emerges (see §3.4).

3 K-way clustering

The above mainly focus on 2-way clustering. Below we gen­ eralize to K -way clustering, K ;::: 2. The generalization is based on the key observation that the solution for cluster in­ dicator vectors (see §2.4 and §2.5) are scaled principal com­ ponents, as we discuss next.

3.1 Scaled principal components

Associations among data objects are mostly quantified by a similarity metric. The scaled principal component approach starts with a nonlinear (non-uniform) scaling of W . Noting that W = D112(D-112WD-112)D112, we apply spectral decomposition on the scaled matrix

instead of on W. This leads to

Here we call qk = D-1 1 2zk the scaled principal compo­ nents. they are obtained by solving the eigenvalue system

Note Eq.(27) is identical to Eq.(20), or Eq.(21), thus scaled principal components are cluster indicator vectors in Min­ MaxCut (see §2.3).

3.2 Self-aggregation network

In Hopfield networks, a pattern q1 is encoded into the net­ work as q 1 qf (the Hebb rule); multiple patterns are encoded additively. In our problem, a pattern is a cluster partitioning indicator vector. We define a self-aggregation network with the connection weights

Q K = ( q1, · · · , qK). Here we highlights several important properties of WsA and provides several example applica­ tions. (Self-aggregation is first studied in [Ding et al., 2002].)

WsA has an interesting self-aggregation property enforced by within-cluster association (connectivity). To prove, we ap­ ply perturbation analysis by writing w = w<o) + w<l)' where w<o) is the similarity matrix when clusters are well­ separated (zero-overlap) and W( 1 ) accounts for the overlap among clusters and is treated as a perturbation[Ding et al., 2001a]. This perturbation approach is standard in quantum physics[Mathews and Walker, 1971].

Well-separated clusters

First, we consider the case where clusters are well sepa­ rated, i.e., no overlap (no connectivity) exist between dif­ ferent clusters. Let W = (Wpq), p, q = 1, · · ·, K. be the subblock decomposition of W. No overlap is represented by the fact that off-diagonal blocks are zero: Wpq = O,p =/= q. Now, W - D = diag(Wu - Du, · · ·, WKK - DKK) The eigenvalue problem decouples into K independent system (Wkk - Dkk)q = >..Dkkq. Clearly, ek = (1, · · ·, 1Y with the size of cluster Gk, is the eigenvector for Gk . Let 8pq =

LiEGp LjEG· W;j, and C = (eke) be an arbitrary orthonor­ mal matrix. Then

are the K eigenvectors with eigenvalue (k = 1. qk is a lin­ ear combination of K step functions, i.e., piece-wise constant function. Clearly, all data objects within the same cluster have identical elements in q. The coordinate of object i in the K­ dim SPCA space is r; = (q1(i), · · ·, qK(i)Y. Thus objects within a cluster are located at (self-aggregate into) the same point in SPAC space.

Scaled principal components are not unique, since C could be any orthogonal matrix. However,

is unique. Thus, self-aggregation of cluster member is equiv­ alent to the fact that QKQ~ has a block diagonal structure, where elements within the same diagonal block all have the same value.

It happens that two objects i , j within the same cluster are not initially connected, ie., w;1 = 0. However, (QKQ'£.);1 gets the same connection strength as any other pairs within the same cluster, i.e., SA network performs a transitive clo­ sure on each connected component. Thus, in SA network, the within-cluster connections are enhanced.

Clusters overlap

Consider the case when overlaps among different clusters exist. The overlaps are treated as a perturbation. Theorem 3. At the first order, the K scaled principal components and their eigenvalues have the form

where y and ( satisfy the eigensystem ry = >..y. The matrix r has the form r = n- 1 1 2 r n- 1 1 2 ' where

hkk = Lp,p# 8kp and n = diag(8u , ... ' 8K)· This analysis is accurate to order IIW(l)II 2 /IIW(o)ll 2 for eigenvalues and to order IIW( 1 ) 11/IIW(O)II for eigenvectors. D

Several features of SPCA can be obtained from Theorem 3: Corollary 3.1. SA network QKQk has the same block di­ agonal form of Eq.(30), within the accuracy of Theorem 2. This assures that objects within the same cluster will self­ aggregate as in Theorem 1.

Corollary 3.2. The first scaled principal component is q1 = D- 1 1 2 XKY1 = (1, · · ·, 1Y with >..1 = 0, (1 = 1. >..1 and q1 are also the exact solutions to the original Eq.(20).

Corollary 3.3. When K = 2, the second principal component is

Its eigenvalue is

Note that this is precisely J3 ofEq.l8, the clustering objective we started with: our clustering framework is consistent.

Example 1. A dataset of 3 clusters with substantial ran­ dom overlap between the clusters. All edge weights are 1. The similarity matrix and results are shown in Fig.l, where nonzero matrix elements are shown as dots. The exact >.2 and approximate .\2 from Theorem 3 are close:

WsA is much sharper than the original weight matrix W clearly due to self-aggregation: connections between differ­ ent clusters are substantially suppressed while connections within same clusters are substantially enhanced.

Application 1. In DNA micro-array gene expression pro­ filing, responses of thousands of genes from tumor tissues are simultaneously measured. We apply SPCA framework to gene expression profiles of lymphoma cancer samples[Al­ izadeh et al., 2000]. Three cancer and three normal subtypes are shown in Fig.2. This is a difficult case due to the large variations of cluster sizes (the number of samples in each sub­ type are shown in parentheses in Fig.2B). Self-aggregation is evident in Figure 2B and 2C. The computed clusters corre­ spond quite well to the normal and cancer subtypes identified by human experts.

3.3 Dynamic aggregation

The self-aggregation process can be repeated to obtain sharper clusters. WK is the low-dimensional projection that contains the essential cluster structure. Combining this struc­ ture with the original similarity matrix, we obtain a new sim­ ilarity matrix containing sharpened cluster information:

where Wkt) = DWsAD for wCt), the weight matrix at t-th iteration, Setting Ak = 1 is crucial for enforcing the cluster structure. w(l) = w and 0: = 0.5.

Applying SA net on W( 2 ) leads to further aggregation (see Figure 2C). The eigenvalues of the 1st and 2nd SA net are shown in the insert in Figure 1 C. As iteration proceeds, a clear gap is developed, indicating that clusters becoming more sep­ arated.

Noise reduction

SA net has noises. For example, WsA has sometimes nega­ tive weights (WsA)ij whereas we expect them to be nonnega­ tive. However, by Corollaries 2.1 and 3.1, WsA has a diagonal block structure and every elements in the block are identical (Eq.30) even when overlaps exit. This property allows us to interpret WsA = QQT as the probability that two objects i, j belong to the same cluster:

To reduce noise in the dynamic aggregation, we set

where 0 < {J < 1 and we chose {J = 0.8. Noise reduction is an integral part of SA net. In our experiments, final re­ sults are insensitive to o:, {J. The above dynamic aggregation repeats self-aggregation process and forces data objects move towards the attractors, which are the desired clusters and their principal eigenvalues approach 1 (see insert in Fig.2C). Usu­ ally, after one or two iterations the cluster structure becomes evident.

3.4 Bipartite graphs

SPCA applies to bipartite graph problems as well . The non­ linear scaling factors are Dr and De (cf. Eq.23) . Let B = D; 1 \B)D~/ 2 , where B is given in Eq.(25). Applying SVD on B, we obtain

Scaled principal components are fk = D;: 1 1 2 uk for row objects and gk = D-; 112 vk for column objects. They have the same self-aggregation and related properties. We note that SPCA on bipartite graph leads to correspondence analysis[Greenacre, 1984] from a new perspective.

The structure of self-aggregation network QKQ~ can be meaningfully further decomposed:

where FK = (f1, · · ·, fK) and GK = (g1, · · ·, gK) · FKF[ provides the cluster structure for row objects, while �K�~ provides the cluster structure for column objects. The off­ diagonal block matrix FKGk provides the sharpened associa­ tion between row and column objects[Ding, 2003].

4 Discussions

In this paper, we present a data clustering framework based on properly motivated Hopfield network, and show the clus­ tering indicators are principal components. Further motivated by cluster balance, we extend the framework to MinMaxCut that utilized the Laplacian matrix of a graph. The frame­ work is generalized to multi-way clustering, by using scaled principal components, and self-aggregation networks are con­ structed. We prove the cluster member self-aggregation prop­ erty of the network. This framework extends naturally to bi­ partite graphs which leads to row-row, column-column and row-column SA nets that simultaneously cluster the row and column objects.

In self-aggregation, data objects move towards each other guided by connectivity. This is similar to the self-organizing map [Kohonen, 1989], where feature vectors self-organize into a 2D feature map while data objects remain fixed. All these have a connection to recurrent networks [Hopfield, 1982; Haykin, 1998 2nd ed]. In Hopfield network, features are stored as associative memories. In more complicated net­ works, connection weights are dynamically adjusted to learn or discover the patterns. The self-aggregation network pro­ vides a new mechanism to realize this unsupervised learning.