Regularity-based spectral clustering and mapping the Fiedler-carpet

1 Introduction

The spectral graph theory started developing about 50 years ago (see, e.g., Hoffman [1], Fiedler [2], Cvetković et al. [3], and Chung [4]), with the goal of characterizing certain structural properties of a graph by means of the eigenvalues of its adjacency or Laplacian matrix. In the last decades of the 20th century, the eigenvectors corresponding to the smallest few eigenvalues of the Laplacian matrix were used to obtain clusterings of the vertices into disjoint parts so that the inter-cluster relations are negligible compared to the intra-cluster ones. The motivation for finding such clusterings comes from machine learning, where the classification of data into a small number of highly similar clusters is a common task. In this setup, the famous Fiedler-vector, the eigenvector, corresponding to the smallest positive Laplacian eigenvalue was used to classify the vertices into two parts. To find a clustering with k≥2 parts, one can instead use the k−1 eigenvectors corresponding to the k−1 smallest positive eigenvalues. It was also noted in [5] that it is better to have the more eigenvectors but no exact estimates between the spectral gap and the quality of the clustering were available, except the k=1 and k=2 cases; the former case is related to the isoperimetric number and expander graphs (see, e.g., [4,6,7]), while the latter one is related to the sum of the inner variances of 2-clusterings (see [8]).

Since then, the problem has been generalized in several ways: to edge-weighted graphs and rectangular arrays of nonnegative entries (e.g., microarrays in biological genetics and forensic science [9,10, 11,12]), and to degree-corrected adjacency and Laplacian matrices [13,14]. Starting in the 21st century, physicists and social scientists introduced modularity matrices and investigated the so-called anticommunity structures (where intracluster relations are negligible compared to the intercluster ones) in contrast to the former community structures [15]. By uniting these two approaches, so-called regular cluster pairs with small discrepancy can be defined, where one looks for homogeneous clusters (e.g., in microarrays, one looks for groups of genes that similarly influence the same groups of conditions), see [10]. The existence of such a regular structure is theoretically guaranteed by the Abel-prize winner Szemeredi’s regularity lemma [16], which, for any small positive ε, guarantees a number k=k(ε) of clusters (independent of the numbers of vertices) such that by partitioning the vertices into k parts (sometimes requiring possibly a small exceptional part), the pairs of clusters have discrepancy less than ε. However, this k guaranteed by the regularity lemma can be enormously large and not applicable for practical purposes. Our goal is to give a moderate k, where the sum of the inner variances of 2k−1 clusters is estimated as mentioned earlier by the spectral gap between the (k−1)st and kth smallest positive normalized Laplacian eigenvalues, even in the worst-case scenario. This is the generalization of a theorem of [8] that was applicable only to the k=2 situation. In that situation, 2k−1=k, and so the Fiedler vector was used for clustering into k=2 parts. If k>2, then we use k−1 nontrivial eigenvectors, the ensemble of which form what we call Fiedler-carpet. In special cases, e.g., in the case of generalized multiclass random or quasirandom graphs, both the objective function of the k-means algorithm and the k-way discrepancy dramatically decrease compared to the values for k−1 [17], where the number of clusters is one more than the number of eigenvectors used in the representation. However, in the generic case, the number of eigenvectors used for the classification is much smaller than the number of clusters, which is promising from the point of view of computational complexity.

The organization of this article is as follows. In Section 2, we define the most important notions and facts concerning normalized Laplacian spectra and multiway discrepancy of rectangular arrays with nonnegative entries [10,11,13,17,18]. In Section 3, the main theorem is stated and proved. The proof is based on the energy minimizing representation, and we analyze the structure of the vertex representatives and the underlying spectral subspaces that are mapped in a convenient way. In Appendix A, a few technical details, simulation results, and plots of the so-obtained Fiedler-carpet are presented. In Section 4, an application to the directed graph of a migration dataset is discussed. This discussion is supplemented with more figures in Appendix C. Appendix B contains pseudocode for the algorithms, which are discussed. In Section 5, the main contributions are summarized and conclusions are drawn.

2 Minimal and regular cuts versus spectra

Let W be the n×n edge-weight matrix of a graph G on n vertices. It is symmetric and has 0 diagonal and nonnegative entries. In the case of a simple graph, W=(wij) is the usual adjacency matrix. The generalized degrees are di=∑j=1nwij, for i=1,…,n. Assume that dis are all positive and the diagonal degree-matrix D contains them in its main diagonal. The Laplacian of G is L=D−W, while its normalized Laplacian is

Because of the normalization, LD is not affected by the scaling of the edge-weights. Therefore, we can assume that ∑i=1ndi=1. LD is positive semidefinite, and if W is irreducible (G is connected), then its eigenvalues are

with unit-norm pairwise orthogonal eigenvectors u0,u1,…,un−1. In particular, u0=(d1,…,dn)T≕dT.

For 1<d<n, the row vectors of the n×d matrix X∗=(D−1/2u1,…,D−1/2ud) are optimal d-dimensional representatives r1∗,…,rn∗ of the vertices. More precisely, X∗ minimizes the energy function:

where the general vertex representatives r1,…,rn∈Rd are row vectors of the n×d matrix X, and they are required to satisfy the constraints XTDX=Id and ∑i=1ndiri=0. Intuitively, minimizing the energy Qd(X) forces representatives of vertices connected with large edge-weights to be close to each other. The minimum of Qd(X), attained at X∗, is the sum of the d smallest positive eigenvalues of LD.

The weighted k-variance of these representatives is defined as follows:

where Vol(U)=∑j∈Udj, ci=1Vol(Vi)∑j∈Vidjrj is the weighted center of the cluster Vi, and the minimization is over the set Pk of proper k-partitions Pk=(V1,…,Vk) of the vertex set.

We will use the weighted k-means algorithm to approach this minimum. In [19], it is stated that if the data satisfy the k-clusterable criterion (Sk2≤ε2Sk−12 with a small enough ε), then there is a PTAS (polynomial time approximation scheme) for the k-means problem. This is the situation we usually encounter.

It is well known that the sum of the k bottom eigenvalues of the normalized Laplacian matrix estimates from below the k-way normalized cut of G, which is fk(G)=minPk∈Pkf(Pk,G), where

Here, w(Va,Vb)=∑i∈Va∑j∈Vbwij is the weighted cut between the cluster pairs. Since ∑i=1k−1λi is the overall minimum of Qk (on the orthogonality constraints) and fk(G) is the minimum over partition vectors (having stepwise constant coordinates over the parts of Pk), the relation

is easy to prove. This estimate is sharper if the eigensubspace spanned by the corresponding eigenvectors is closer to that of the partition vectors in the convenient k-partition of the vertices, the one produced by the weighted k-means algorithm. Therefore, Sk2(Xk−1∗) indicates the quality of the k-clustering based on the k−1 bottom eigenvectors (excluding the trivial one). Later, it will be used that neither Qk nor Sk2(Xk−1∗) is affected by the orientation of the orthonormal eigenvectors.

In [8], it is proved that S22(X1∗)≤λ1λ2, so the larger the gap after the first positive eigenvalue of LD, the sharper the estimate in (3) is. Here, this statement is generalized to the gap between λk−1 and λk, but instead of k clusters, we must consider 2k−1 ones, which are based on (k−1)-dimensional vertex representatives. This contrasts with the typical situation in the literature (see, e.g., [20,21,22]), where the number of clusters is the same as the number of eigenvectors used for the classification. The message of Theorem 1 of Section 3 is that the number of clusters is much higher in the generic case than the dimension of the representatives, at least in the minimum multiway cut problems. With the discrepancy objective, the famous Szemerédi’s regularity lemma [16] also suggests this phenomenon.

Now we consider the discrepancy view. The theory can be better illustrated by considering rectangular arrays of nonnegative entries; adjacency matrices of simple, edge-weighted, and directed graphs are special cases of such matrices.

In many applications, for example, when microarrays are analyzed, the data are collected in the form of an m×n rectangular matrix C=(cij) of nonnegative real entries. (If the entries are integer frequency counts, then the array is called contingency table in statistics.) We assume that C is nondegenerate, i.e., CCT (when m≤n) or CTC (when m>n) is irreducible. Consequently, the row-sums drow,i=∑j=1ncij and column-sums dcol,j=∑i=1mcij of C are strictly positive, and the diagonal matrices Drow=diag(drow,1,…,drow,m) and Dcol=diag(dcol,1,…,dcol,n) are regular. Without the loss of generality, we also assume that ∑i=1n∑j=1mcij=1, since neither our main object, the normalized contingency table

nor the multiway discrepancies to be introduced are affected by any uniform scaling of the entries of C. It is known that the singular values of CD are in the [0,1] interval. The positive singular values, enumerated in nonincreasing order, are the real numbers denoted by

where r=rank(CD)=rank(C). Provided C is nondegenerate, 1 is a single singular value; it is called the trivial singular value and is denoted by s0 since it corresponds to the trivial pair of singular vectors, which are disregarded in clustering problems. This is a well-known fact of correspondence analysis; for further explanation, see [13,23] and the subsequent paragraph.

For a given integer 1≤k≤min{m,n}, we are looking for k-dimensional representatives r1,…,rm∈Rk of the rows and q1,…,qn∈Rk of the columns of C such that they minimize the energy function:

subject to

When minimized, the objective function Qk favors k-dimensional placement of the rows and columns such that representatives of columns and rows with large association (cij) are close to each other. It is easy to prove that the minimum is obtained by the singular value decomposition (SVD):

where r≤min{n,m} is the rank of CD. The constrained minimum of Qk is 2k−∑i=0k−1si, and it is attained with row- and column-representatives that are row vectors of the matrices Drow−1/2(v0,v1,…,vk−1) and Dcol−1/2(u0,u1,…,uk−1), respectively.

Note that if the entries of C are frequency counts and their sum (N) is the sample size, then the χ2 statistic, which measures the deviation from independence, is

If the χ2 test based on this statistic indicates significant deviance from independence (i.e., from the rank 1 approximation of C), then one may look for the optimal rank k approximation (1<k<r=rank(C)), which is constructed using the first k singular vector pairs.

When bi-clustering the rows and columns of C one may also look for subtables close to independent ones. This is measured by the discrepancy. In [10], the multiway discrepancy of the rectangular matrix C of nonnegative entries in the proper k-partition R1,…,Rk of its rows and C1,…,Ck of its columns is defined as follows:

where c(X,Y)=∑i∈X∑j∈Ycij is the cut between X⊂Ra and Y⊂Cb, Vol(X)=∑i∈Xdrow,i is the volume of the row-subset X, and Vol(Y)=∑j∈Ydcol,j is the volume of the column-subset Y, whereas ρ(Ra,Cb)=c(Ra,Cb)Vol(Ra)Vol(Cb) denotes the relative density between Ra and Cb. The minimum k-way discrepancy of C itself is expressed as follows:

In [10], the following is proved:

provided 0<mdk(C)<1. In the forward direction, the following is established in [13]. Given the m×n contingency table C, consider the spectral clusters R1,…,Rk of its rows and C1,…,Ck of its columns, obtained by applying the weighted k-means algorithm to the (k−1)-dimensional row- and column-representatives. Let Sk,row2 and Sk,col2 denote the minima of the weighted k-means algorithm, applied to the rows and columns, respectively. Then, under some balancing conditions for the margins and for the cluster sizes, mdk(C)≤B(2k(Sk,row+Sk,col)+sk), where B is some constant that depends only on the balancing conditions, and does not depend on m and n. Roughly speaking, the two directions together imply that if sk is “small” and “much smaller” than sk−1, then one may expect a simultaneous k-clustering of the rows and columns of C with small k-way discrepancy. This is the case for generalized random and quasirandom graphs [17].

This notion can be extended to an edge-weighted graph G and denoted by mdk(G). In that setup, C plays the role of the weighted adjacency matrix. Here, the singular values of the normalized adjacency matrix are the absolute values of the eigenvalues, which are encountered in the decreasing order.

At this point, we introduce some new matrices, originally defined by physicists (see [15]). The modularity matrix of an edge-weighted graph G is defined as M=W−ddT, where the entries of W sum to 1. The normalized modularity matrix of G (see [24]) is expressed as follows:

The normalized modularity matrix is the normalized edge-weight matrix deprived of the trivial dyad. Obviously, LD=I−WD=I−MD−ddT, see [13].

Therefore, the k−1 largest singular values of MD are the absolute values of the k−1 largest absolute value eigenvalues of WD (except the trivial 1). These, in turn, are 1 minus the k−1 positive eigenvalues of LD, which are farthest from 1. If those are all less than 1, then these are 1−λ1,…,1−λk−1. In this case, the regularity-based spectral clustering boils down to the minimum cut objective.

Otherwise, for a 1<k<n integer fixed, in the modularity-based spectral clustering, we look for the proper k-partition V1,…,Vk of the vertices such that the within- and between-cluster discrepancies are minimized. To motivate the introduction of the exact discrepancy measure, observe that the ij entry of M is wij−didj, which is the difference between the actual edge-weight between the vertices i and j and the edge-weight that is expected under independent attachment of them with probabilities di and dj, respectively. Consequently, the difference between the actual and the expected total edge-weight between the subsets X,Y⊂V is expressed as follows:

A directed edge-weighted graph G=(V,W) is described by its quadratic, but usually not symmetric weighted adjacency matrix W=(wij) of zero diagonal, where wij is the nonnegative weight of the j→i edge (i≠j). The row-sums din,i=∑j=1nwij and column-sums dout,j=∑i=1nwij of W are the in- and out-degrees, while Din=diag(din,1,…,din,n) and Dout=diag(dout,1,…,dout,n) are the diagonal in- and out-degree matrices. The multiway discrepancy of the directed, edge-weighted graph G=(V,W) in the in-clustering Vin,1,…,Vin,k and out-clustering Vout,1,…,Vout,k of its vertices is expressed as follows:

where w(X,Y) is the sum of the weights of the X→Y edges, whereas Volin(X)=∑i∈Xdin,i and Volout(Y)=∑j∈Ydout,j are the out- and in-volumes, respectively. The minimum k-way discrepancy of the directed edge-weighted graph G=(V,W) is expressed as follows:

In [10], it is proved that

where sk is the kth largest nontrivial singular value of the normalized weighted adjacency matrix WD=Din−1/2WDout−1/2. In Section 4, we apply the SVD of WD to find migration patterns, i.e., emigration and immigration trait clusters.

3 Mapping the Fiedler-carpet: more clusters than eigenvectors

Note that only if λk−1<λk, uk and xk are not in the subspace spanned by u1,…,uk−1. Theorem 1 indicates the following clustering property of the (k−1)st and kth smallest normalized Laplacian eigenvalues: the greater the gap between them, the better the optimal k-dimensional representatives of the vertices can be classified into 2k−1 clusters.

Figure 1 shows a graph with three well-separated clusters, as well as the image of the corresponding map f:R2→R2 of the Fiedler-carpet as used in the proof of Theorem 1.

Figure 1

A graph with three well-separated clusters, and the image of its Fiedler-carpet.

The image of f contains the origin, and we can find by inspection a pair (a1,a2) for which f(a1,a2) is approximately zero; namely, choosing a1=−0.19099 and a2=−0.35688 gives f(a1,a2)≈(−0.000002,0.00000001). The two-dimensional representative of vertex i is the point (x1i,x2i) for i=1,…n; these are plotted in Figure 2, where the colors indicate the cluster memberships and the black point denotes the approximate root of f.

Figure 2

Two-dimensional vertex representatives of the graph of Figure 1. The black point denotes the approximate root of f.

Some remarks are in order:

1. Theorem 1 is the generalization of Theorem 2.2.3 of [8]. In the k=2 case, 2k−1=k, but in general, the number of clusters is much larger than that of the relevant eigenvectors.

2. The statement of the theorem has relevance, since for any k>0, the relation ∑j=1k−1λj<(k−1)λk holds; but with analysis of variance considerations, S2k−12≤k−1 also holds. In particular, when k=2, only one eigenvector is used for the representation. The total inertia of the coordinates is 1, and it can be divided into the sum of nonnegative within- and between-cluster inertias. The within-cluster inertia is the sum of the inner variances of the two clusters, which is S22, and it is at most 1. (Here, the variances are calculated with respect to the discrete distribution d1,…,dn.) When (k−1)-dimensional representatives are used, then the total inertia is tr(Xk−1TDXk−1)=tr(Ik−1)=k−1, and the sum of the inner variances is again at most k−1, but it is further bounded with ∑j=1k−1λjλk by Theorem 1.

3. The vector u1 is called Fiedler-vector of the (nonnormalized) Laplacian. Here, we use the ensemble of the first k−1 transformed eigenvectors of the normalized Laplacian together, the so-called Fiedler-carpet.

4 Applications

We investigated the international migrant stock by the country of origin and destination in the years 2015 and 2019. The focus is on 41 European countries plus the United States of America and Canada. The data^[1] are based on the official registered migrants numbers, where columns and rows correspond to the country of origin and the country of destination, respectively.

Here, the quadratic, but not symmetric edge-weight matrix contains weights of bidirected edges (the diagonal is zero): the i,j entry is the number of persons going j→i. Via SVD of the normalized table, we can find emigration (column) and immigration (row) clusters, between which the migration is the best homogeneous (in terms of discrepancy).

Both for the 2015 and 2019 data, there was a gap after four nontrivial singular values, and therefore, the corresponding four singular vector pairs were used to find five emigration and immigration trait clusters.

• Singular values, 2015:

  1, 0.79098, 0.71857, 0.67213, 0.56862, 0.45293, 0.40896, 0.38178, 0.36325, 0.34785, 0.32648, 0.31769, 0.2996, 0.27927, 0.26566, 0.24718, 0.22638, 0.20632, 0.18349, 0.1651, 0.14384, 0.1359, 0.12721, 0.12092, 0.11816, 0.10374, 0.09545, 0.08278, 0.0738, 0.06371, 0.05673, 0.04553, 0.03488, 0.03107, 0.02967, 0.02693, 0.01557, 0.00788, 0.00584, 0.00519, 0.00191, 0.0017, 0.00099.

• Singular values, 2019:

  1, 0.77844, 0.70989, 0.65059, 0.55122, 0.43612, 0.39512, 0.36194, 0.3558, 0.33882, 0.32174, 0.30719, 0.29601, 0.28181, 0.26865, 0.259, 0.22421, 0.1917, 0.17988, 0.1516, 0.13671, 0.13243, 0.12397, 0.11542, 0.10598, 0.09216, 0.08889, 0.07958, 0.06835, 0.06154, 0.05377, 0.04412, 0.03436, 0.03124, 0.02899, 0.02745, 0.01507, 0.00814, 0.00619, 0.0051, 0.00216, 0.00129, 0.00089.

In the discrepancy-based spectral clustering, the number of clusters is one more than the number of structural singular values (excluding the trivial 1). In our case, the number of clusters was five. The five emigration and immigration trait clusters for 2015 are shown in Tables 1 and 2, whereas those for 2019 are shown in Tables 3 and 4. Figure 3(a) and (b) illustrate the immigration–emigration cluster-pairs with the countries rearranged by their cluster memberships. The frequency counts are represented with light to dark squares.

Figure 3

Immigration–emigration cluster-pairs with the countries rearranged by their cluster memberships. The frequency counts are represented with light to dark squares. (a) 2015 and (b) 2019.

In 2015, by χ2 test, we found the smallest discrepancy between the emigration trait cluster number 2 and immigration trait cluster number 4, i.e., the subtable formed by them was close to an independent table of rank 1. The clusters were similar in the 2 years; in 2019, the smallest discrepancy was found between the emigration trait cluster number 2 and immigration trait cluster number 5, i.e., between the Balcanian and Baltic countries in both years. Correspondence analysis results with cluster memberships are found in Appendix C.

Here, we wanted to find as homogeneous cluster pairs as possible, as for migration patterns vis-a-vis the emigration–immigration. Our graph was a relatively small and sparse one, so Figure 3(a) and (b) show a large and some small clusters for both emigration and immigration. However, the emigration–immigration cluster pair with the smallest discrepancy is spotted by the χ2 test. In larger and more dense networks, e.g., in metabolic networks with bidirected edges between the vertices (enzymes), the weights representing the intensity of chemical reactions, so-called autocatalytic subnetworks could be found with our method, based on discrepancies.

5 Summary and author contributions

Spectral clustering is a collection of methods for clustering data points or vertices of a graph based on combinatorial criteria with spectral relaxation. Here, we generalize spectral clustering in several ways:

1. Instead of simple or edge-weighed graphs, we consider directed graphs and rectangular arrays of nonnegative entries. For these, the method of correspondence analysis gives low-dimensional representation of the row and column items by means of SVD of the normalized table. This is applied to real-life data.

2. Instead of the usual multiway minimum cut objective, we consider discrepancy minimization, for which a near optimal solution is given by the k-means algorithm applied to the row and column representatives, where the number of clusters (k) is concluded from the gap in the singular spectrum. There are theoretical results supporting this. Eventually, the near optimal clusters can be refined to decrease discrepancy.

3. The number of clusters can be larger than the number of eigenvectors entered in the classification. The two are the same only in case of generalized random or quasirandom graphs, see [17]. In Theorem 1, an exact estimate for the sum of the inner variances of 2k−1 clusters is given by means of the spectral gap between the (k−1)st and kth smallest normalized Laplacian eigenvalues by using the Fiedler-carpet of the corresponding eigenvectors.

4. When the number of clusters is the same as the number of eigenvectors entered in the classification, Davis-Kahan type subspace perturbation theorems are applicable, as in this case, the k-partition vectors span a k-dimensional subspace. However, if the number of parts is larger than that of the eigenvectors, then finer methods are applicable, like our construction in the proof, using a “witness.”