Consistency between ordering and clustering methods for graphs

We introduce a measure to quantify the consistency between a sequence $\pi$ and partition $\sigma$. We define the sequence $\pi$ as consistent with $\sigma$ if vertices with the same group label are maximally adjacent to each other in the sequence $\pi$. For instance, if the original indices $I$ are consistent with group labels $\sigma$,

is maximized, where $\delta (a,b)$ represents the Kronecker delta; Fig. 3(a) presents an example. To evaluate the consistency between $\pi$ and $\sigma$, we introduce a measure that we refer to as the label continuity, defined by

where ${\pi }^{-1}$ is the inverse mapping of $\pi$ and $i^{\prime }$ is the index label after the permutation; that is, ${\pi }^{-1}\left(i^{\prime }\right)$ is the label $i$ in the original indices satisfying $\pi \left(i\right)=i^{\prime }$. The number of times that the group labels are flipped when following the vertices in the order of $\pi$ is expressed as $(N-1)(1-C(\text{boldmath}\pi ,\text{boldmath}\sigma ))$, for which the group labels must be flipped at least $K-1$ times. Considering this feature, we define the label continuity error (LCE) as

Hereafter, we abbreviate $C(\text{boldmath}\pi ,\text{boldmath}\sigma )$ and $\Delta (\text{boldmath}\pi ,\text{boldmath}\sigma )$ as $C$ and $\Delta$, respectively, as long as there is no possibility of confusion.

For a given partition $\sigma$ and different vertex sequences, we can evaluate which vertex sequence is more consistent with $\sigma$ using the LCE (e.g., Figs. 3(a) and 3(b)). Similarly, for a given sequence $\pi$ and different partitions, we can also evaluate which partition is more consistent with $\pi$, keeping the group sizes $\left\{N_k\right\}$ fixed (e.g., Figs. 3(a) and 3(c)).

The LCE can take only small values when the number of groups $K$ is very small or large. For example, it is obvious that $\Delta$ is zero when $K=1$ or $K=N$. In other words, the resolution of the LCE is low in such regions. Moreover, this property would depend on the distribution of the group sizes $\left\{N_k\right\}$. In this section, we quantify these intuitions.

The minimum value of $\Delta$ is zero by construction. The maximum value of $\Delta$ is obtained when labels are flipped the maximum number of times. The maximization of $\Delta$ by optimizing the sequence $\pi$, given an arbitrary partition $\sigma$ with $\left\{N_k\right\}$, is equivalent to maximizing $\Delta$ by optimizing partition $\sigma$ (constrained to $\left\{N_k\right\}$) for a given sequence $\pi$. We denote the maximum by $max\Delta$ as

As derived in Appendix A, we have

where $\lceil \cdot \rceil$ denotes the ceiling function.

We next investigate statistical properties of the LCE. First, we calculate the probability $P\left(m\right)$ at which the number of times that two consecutive vertices in a sequence have the same group label is $m$, where $m=(N-1)C(\text{boldmath}\pi ,\text{boldmath}\sigma )$. When any sequence realizes at random, we have

where $\sigma$ is an arbitrary partition with group sizes $\left\{N_k\right\}$ and the sum is over all possible sequences ($|\left\{{\text{boldmath}\pi }^{\prime }\right\}|=N!$). Note that Eq. (6) is also a distribution in which each distinguishable partition realizes at random. This equivalence might sound peculiar because there are only $N!/{\prod }_{k=1}^K{N}_k!$ distinguishable partitions, whereas there are $N!$ possible sequences. However, because every distinguishable partition is overcounted exactly ${\prod }_{k=1}^K{N}_k!$ times in the summation of Eq. (6), the distribution $P\left(m\right)$ is identical for both random sequences and random partitions.

Although Eq. (6) is a straightforward expression, a strict constraint on $\left\{N_k\right\}$ makes analytical calculations complicated. Therefore, we instead calculate the distribution of bootstrapped group labels ${\text{boldmath}\sigma }^{\ast }$ as an approximation. That is, we generate a random group assignment ${\text{boldmath}\sigma }^{\ast }$ by sampling independently from the empirical distribution $Prob\left[k\right]=N_k/N$ ($k\in \{1,\dots ,K\}$); in other words, we randomly resample group labels from $\sigma$ with replacement. The distribution of group labels ${\text{boldmath}\sigma }^{\ast }$ is

This approximation for random group labels is expected to be accurate if each element in $\left\{N_k\right\}$ is sufficiently large.

Using the bootstrapped group labels, the mean value of $C$ is obtained as

Therefore, the mean value of LCE under random partitioning is

As the LCE does not practically become worse than $\overline{\Delta }\left(\right\{N_k\left\}\right)$, this mean value is a more meaningful reference value than Eq. (5) as the upper bound.

We can also derive the variance $Var\left[\Delta \right]$ (the derivation is shown in Appendix B) as

showing that $\Delta$ converges to $E\left[\Delta \right]$ by the law of large numbers. Furthermore, in Appendix C, we show that the probability distribution is asymptotically normal when the group sizes are equal and $K=O\left(1\right)$, implying that higher-order moments will vanish.

Let us summarize the results we obtained in this section. As the number of groups $K(>1)$ increases, the upper bound of the LCE ($max\Delta$ in Eq. (5)) decreases monotonically as long as the partitions are not highly skewed, i.e., ${max}_k{N}_k<\lceil N/2\rceil$. However, as illustrated in Fig. 4, the LCE for a random sequence ($\overline{\Delta }$ in Eq. (9)) is a convex function with respect to $K$. When $K$ is small, the LCE increases because the chance for label flips increases, while the LCE decreases owing to the increase in the minimum number of label flips. For equipartitioning (Fig. 4(a)), the mean LCE $\overline{\Delta }$ is peaked at an integer of approxmately $K=\sqrt{N-1}$. As a partition becomes more skewed (Fig. 4(b)), $max\Delta$ and $\overline{\Delta }$ are peaked at larger values of $K$. Therefore, when evaluating the LCE, we must implement appropriate normalizations.

In this study, we focus on comparing partitions with the same number of groups $K$. In Appendix D, however, we discuss nested partitions (subpartitions of another partition) as an example in which different partitions have different numbers of groups.

Spectral ordering is derived as a continuous relaxation of the discrete optimization problem called envelope reduction [14]. This problem optimizes the vertex sequence $\pi$ such that each connected pair of vertices is located close to each other in the sequence. To this end, the following objective function is considered:

which is the sum of squared distances ${({\pi }_i-{\pi }_j)}^2$ with respect to the set of connected vertices. The sequence that minimizes this function is the solution to envelope reduction.

As the minimization of Eq. (11) is not computationally feasible, we consider its continuous relaxation. That is, we represent $\pi$ using a continuous vector $\text{boldmath}x\in R^N$. However, if we simply replace $\pi$ with $x$, $\text{boldmath}x=\text{boldmath}0$ would be the trivial minimizer of $H_2(\text{boldmath}x;A)$. Thus, we constrain $x$ such that ${\sum }_{i=1}^N{x}_i^2$ is a positive constant (i.e., the spherical constraint) to reflect the fact that ${\sum }_{i=1}^N{\pi }_i^2$ is positive regardless of the choice of sequence. Therefore, we consider the minimization of the following function:

where $\lambda$ is the Lagrange multiplier. The extremum condition in Eq. (12) yields the following eigenvalue equation with an eigenvector $\nu$:

Here, $L\equiv D-A$ is the unnormalized (or combinatorial) Laplacian, where $D=diag(d_1,\dots ,d_N)$ is the degree matrix ($d_i={\sum }_{j=1}^N{A}_{ij}$). Although we would have a vector proportional to $1$ (a vector of ones) as the minimizer of Eq. (12), which is also the eigenvector associated with the smallest eigenvalue of $L$, we cannot infer the optimal sequence from $1$ because all the elements are identical. Therefore, we exclude vectors proportional to $1$, which is equivalent to imposing a perpendicular constraint to $1$ in Eq. (12), i.e., ${\sum }_{i=1}^N{x}_i=0$. Then, the minimizer of the objective function is the eigenvector ${\text{boldmath}\nu }_2$ of $L$ associated with the second-smallest eigenvalue.

The estimate of the optimal sequence $\widehat{\text{boldmath}\pi }$ using the spectral ordering method is

where ${\nu }_{2i}$ is the $i$th element of ${\text{boldmath}\nu }_2$, and $rank\left({\nu }_{2i}\right)$ is the index of ${\nu }_{2i}$ in an array in which the vector elements of ${\text{boldmath}\nu }_2$ are sorted in the ascending or descending order.

A spectral ordering method with the normalized Laplacian was derived in [15]. Note that the objective function (11) does not have a periodic boundary condition. Therefore, while the distance from one vertex at the end of the sequence to another vertex ranges from $1$ to $N-1$, the distance from the vertex at the middle of the sequence ranges from $1$ to $\lfloor N/2\rfloor$, where $\lfloor \cdot \rfloor$ is the floor function. This implies that when a graph has a vertex with a considerably large degree (i.e., a hub), it is typically more beneficial for the minimization objective to assign such a vertex near the middle of the sequence. To incorporate this feature, we replace the spherical constraint in Eq. (12) with the following ellipsoidal constraint:

which tends to restrict $x_i$ with a large $d_i$ to be relatively small (recall that a variable with a large coefficient typically has relatively small values on an ellipsoid). Therefore, Eq. (15) constrains $x$ such that $x_i$ of a hub vertex $i$ is near the origin, and when $x$ is discretized, the hub vertices are likely to be located near the middle of the sequence. Note also that the mean of $\left\{x_i\right\}$ is located at the origin because of the perpendicular constraint ${\sum }_{i=1}^N{x}_i=0$.

Consequently, Eq. (13) is replaced with the following generalized eigenvalue equation with respect to its second-smallest eigenvalue ${\lambda }_2$:

This is equivalent to

where $L\equiv D^{-\frac{1}{2}}L{D}^{-\frac{1}{2}}$ is the normalized Laplacian and ${\text{boldmath}z}_2\equiv D^{\frac{1}{2}}{\text{boldmath}\nu }_2$. As ${\text{boldmath}\nu }_2$ is a continuous relaxation of the sequence $\pi$, we estimate the optimal sequence $\widehat{\text{boldmath}\pi }$ as

The modularity matrix $Q$ appears in the spectral clustering method for modularity maximization in community detection [21]. The matrix element is commonly defined as

where $M$ is the total number of edges in the graph.

To formulate the spectral ordering problem with the modularity matrix, we again consider the objective function $H_2(\text{boldmath}\pi ;A)$ in the envelope reduction problem and its continuous relaxation with the spherical constraint ${\sum }_{i=1}^N{x}_i^2=1$. Herein, we add the following penalty terms to the objective function:

The penalty terms ensure that $\left\{x_i\right\}$ are “balanced” around the origin. The first term prohibits $\left\{x_i\right\}$ for hub vertices from being located only on the positive or negative side of the real interval $[-1,1]$. Owing to the second term, $\left\{x_i\right\}$ associated with hub vertices also tend to be away from the origin. Therefore, the penalty term Eq. (20) decreases when $\left\{x_i\right\}$ are more symmetrically distributed around the origin.

Using Lagrange multipliers, the objective function to be minimized is then

Here, we do not impose the perpendicular constraint in Eq. (21), because a vector proportional to $1$ is not a trivial minimizer. The extremum conditions in Eq. (21) yield

where ${\lambda }_1$ represents the largest eigenvalue of $Q$ and ${\text{boldmath}\nu }_1$ is the associated eigenvector. ${\text{boldmath}\nu }_1$ is the minimizer in Eq. (21) provided that it is not a vector proportional to $1$. Analogously to Eq. (14), we estimate the optimal sequence $\widehat{\text{boldmath}\pi }$ as

The ellipsoidal constraint enforces $\left\{x_i\right\}$ for hub vertices to be concentrated around the origin, whereas the penalty terms (20) enforce them to be evenly distributed at both the positive and negative ends of the real line. Therefore, the results of the spectral ordering methods using the normalized Laplacian and modularity matrix are expected to be quite distinct for graphs with heterogeneous degree distributions.

Bethe Hessian is also a matrix that is originally formulated to perform spectral clustering [23, 24]. This method is inspired by the statistical inference of the stochastic block model, which will be explained in Sec. IV.1. This section considers a spectral ordering method using the Bethe Hessian.

The derivation of spectral ordering with the Bethe Hessian is analogous to that with the normalized Laplacian. However, instead of imposing an ellipsoidal constraint (15), we introduce ${\sum }_{i=1}^N{d}_i{x}_i^2$ as a penalty term. Thus, we consider the following objective function:

where $\tau$ is an arbitrary constant (hyperparameter) that can be either positive or negative. To avoid the trivial minimizer $\text{boldmath}x=\text{boldmath}0$, we impose the spherical constraint ${\sum }_{i=1}^N{x}_i^2=1$.

Using Lagrange multipliers, the objective function to be minimized is then

where

is a matrix element of Bethe Hessian. The extremum conditions in Eq. (25) yield an eigenvalue equation with respect to $B$.

We estimate the optimal sequence $\widehat{\text{boldmath}\pi }$ as follows:

where ${\text{boldmath}\nu }_2$ is the eigenvector associated with the second-smallest eigenvalue ${\lambda }_2$, i.e.,

Note that there is no guarantee that ${\text{boldmath}\nu }_2$ always provides the best estimate in terms of $H_2(\widehat{\text{boldmath}\pi };A)$ among all the eigenvectors. In fact, we confirmed that the eigenvector that yields the best estimate in terms of $H_2(\widehat{\text{boldmath}\pi };A)$ (when we employ the rounding rule in Eq. (27)) depends sensitively on the value of $r$, particularly when $r$ is small (see Appendix S1 for details). However, we employ Eq. (27) because the estimate with ${\text{boldmath}\nu }_2$ offers the smallest value of $H_2(\widehat{\text{boldmath}\pi };A)$ as long as $r$ is sufficiently large.

Throughout this study, we set $r=\sqrt{{\sum }_i{d}_i^2/{\sum }_i{d}_i-1}$ ($>1$), because it is a commonly employed value in spectral clustering. The hyperparameter $\tau$ is negative when $r>1$. Thus, $\left\{x_i\right\}$ for hub vertices are aligned near the ends of the real line $[-1,1]$ so that a sequence achieves a lower value of Eq. (24) using the penalty term. By contrast, when $r<1$ ($\tau >0$), $\left\{x_i\right\}$ for hub vertices are likely to be located near the origin, implying that the resulting sequence is similar to that obtained by the spectral ordering method based on the normalized Laplacian.

During the past decade, it has been found that the performance of the Laplacian-based spectral clustering can be considerably improved by adding a constant value to every element in the adjacency matrix [26, 28] or the diagonal elements in the degree matrix [25, 27]. Although the two variants of the Laplacian are often termed differently, we collectively refer to them as the regularized Laplacian [28] for simplicity, and we denote the former version of the regularized Laplacian as $L^{\left(\tau \right)}$ and the latter version as $L$. The spectral clustering method based on $L$ can also be interpreted as a continuous relaxation of the minimization of the core cut function [29]. This section considers the spectral ordering method using a regularized Laplacian.

Similar to the formulation of the spectral ordering method with the modularity matrix, we consider the continuous relaxation of $H_2(\text{boldmath}\pi ;A)$ with a penalty term. We consider the following objective function:

which is minimized with respect to the continuous vector $x$. $\tau$ is an arbitrary positive constant (hyperparameter) and

is the variance with respect to the elements in $x$. To ensure that $x$ is not a vector of zeros, we impose the following ellipsoidal constraint:

By incorporating this constraint, the objective function to be minimized is given by

Because a vector proportional to $1$ is a trivial minimizer of Eq. (32), we also impose the constraint that $x$ is perpendicular to $1$, i.e., ${\sum }_{i=1}^N{x}_i=0$. Then, the extremum conditions in Eq. (32) yield

where $I$ is the identity matrix. ${\lambda }_2$ is the second-smallest eigenvalue of the generalized eigenvalue equation and ${\text{boldmath}\nu }_2$ is the associated generalized eigenvector. Equation (33) is equivalent to

where

is the regularized Laplacian and ${\text{boldmath}z}_2={(D+\tau I)}^{1/2}{\text{boldmath}\nu }_2$. Similar to the spectral ordering method with the normalized Laplacian, we estimate the optimal sequence $\widehat{\text{boldmath}\pi }$ as

The contribution of hub vertices is more complicated than that of other methods. As shown in Fig. 5, whereas $\left\{x_i\right\}$ for the hub vertices tend to be relatively small because of the ellipsoidal constraint, the variance $Var\left[\text{boldmath}x\right]$ is minimized when all $\left\{x_i\right\}$ have the same value. Therefore, when the hyperparameter $\tau$ is small, the result is similar to that obtained by using the spectral ordering method based on the normalized Laplacian. As $\tau$ increases, the hub vertices are less likely to be located in the middle of the sequence because of the penalty term.

As mentioned above, we also consider

as the definition of a regularized Laplacian. Unlike $L^{\left(\tau \right)}$, only the degree matrix is perturbed by a constant value in $L$. If we consider

as the objective function to be minimized, and impose the ellipsoidal constraint (31), we obtain the eigenvalue equation with respect to $L$ as a result of the extremum conditions.

If we also impose the constraint ${\sum }_{i=1}^N{x}_i=0$ in Eq. (38), this objective function becomes equivalent to Eq. (29). We do not have such a constraint because $L$ does not have $1$ as a trivial eigenvector unlike $L^{\left(\tau \right)}$. The eigenvectors of $L^{\left(\tau \right)}$ and $L$ are therefore distinct. As a spectral ordering method with the regularized Laplacian $L$, we replace ${\text{boldmath}z}_2$ in Eq. (36) with the eigenvector associated with the second-smallest eigenvalue of $L$. Here, we use the second-smallest value because $L$ approaches $L$ as $\tau \to 0$, and $L$ has ${\text{boldmath}z}_1\propto \text{boldmath}1$. Hereafter, when we refer to the spectral ordering method with the regularized Laplacian, we employ $L$ because it is more computationally efficient. Throughout this study, we set $\tau$ as the average degree of the graph, as it is a commonly employed value [27].

The constraints and penalty terms for each method are summarized in Table 1. Compared to the classical method based on the normalized Laplacian, where hub vertices are concentrated around the middle of the sequence (“hub-centered”), the spectral ordering methods obtained with the modularity matrix, Bethe Hessian, and the regularized Laplacian may assign hub vertices at both ends of the sequence (“hub-at-the-corner”). Particularly, for the Bethe Hessian and the regularized Laplacian, we can choose the “hub-centered” or “hub-at-the-corner” alignment by tuning the hyperparameter. The next section investigates the practical performance of these spectral ordering and clustering methods using synthetic and real-world datasets. Although one might expect all the methods to work similarly when the graph is close to regular, it is not trivial to determine whether this always holds; this is investigated using synthetic datasets in Secs. IV.1 and IV.2. The effect of heterogeneous degree distribution in each method is examined using real-world datasets in Sec. IV.3.

The SBM is often used as a generative model for the inference of module structures in graphs [33, 34] and in several theoretical studies in the community detection literature [35]. In the SBM, each vertex has a “planted” (or preassigned) group assignment; we denote the corresponding partition as ${\text{boldmath}\sigma }^B$. Each vertex pair is connected by an edge, independently and randomly, based on the planted group assignments. The probabilities for the upper-right elements of the adjacency matrix are given as follows:

where $p_{k\ell }$ is the probability that a vertex in group $k$ and vertex in group $\ell$ are connected (in Eq. (39), $k={\sigma }_i^B$ and $\ell ={\sigma }_j^B$). We have $A_{ij}=A_{ji}$ for any pair of elements because we consider undirected graphs. In general, the SBM can generate graphs with complex module structures. Herein, however, we focus on the SBM with a community structure that is characterized by the following group-wise connection probability:

where $0<p_{out}\le p_{in}\le 1$, that is, vertices are more densely connected within the same planted group than between different groups. In particular, when the group sizes are equal, it is common to parametrize the model using the average degree $c$ and the fuzziness parameter $ϵ$, which are related to $p_{in}$ and $p_{out}$ as

As $ϵ$ approaches unity, the planted community structure becomes less clear. This particular case of the SBM is known as the planted partition model [36]. For a given average degree $c$, the critical value of $ϵ$ above which an algorithm cannot detect the planted block structure better than chance is called the (algorithmic) detectability limit [37, 38, 39, 40, 41, 42].

Using the SBM and spectral ordering methods, we investigate the following questions:

1. How would the reordered adjacency matrix look like? Can we visually identify the community structure through the matrix?

2. When and how would the spectral ordering methods lose their correlations with the planted partition in the SBM? Are the spectral ordering methods superior or inferior to their clustering counterparts in detecting the planted partition? Does the choice of matrix matter?

To answer these questions, we apply both spectral ordering and clustering methods to graphs generated by the SBM.

We first investigate the former question. Figure 6 shows the results of spectral ordering applied to instances of SBM. Vertices in the same planted group are indeed located closely in the inferred sequence when the community structure is strong. Even when the community structure is weak, the planted group labels and the inferred sequence are correlated. In both examples, the boundaries of the groups are ambiguous. Therefore, if we do not know the planted group labels (and without the coloring of the adjacency matrix elements), it is not clear whether the identified structure is a community structure or a banded structure from the reordered adjacency matrix. Note that, as discussed in [10], even when we generate a graph from a uniformly random graph model, one could identify a weak banded structure owing to the ordering of vertices.

Next, we address the latter question. The consistency between the sequence inferred by a spectral ordering method $\widehat{\text{boldmath}\pi }$ and planted partition ${\text{boldmath}\sigma }^B$ is measured with the normalized LCE $\Delta (\widehat{\text{boldmath}\pi },{\text{boldmath}\sigma }^B)/\overline{\Delta }\left(\right\{N_k^B\left\}\right)$. Here, $\left\{N_k^B\right\}$ is the set of group sizes in ${\text{boldmath}\sigma }^B$ and $\overline{\Delta }\left(\right\{N_k^B\left\}\right)$ is the LCE under a random sequence defined in Eq. (9). When $\Delta (\widehat{\text{boldmath}\pi },{\text{boldmath}\sigma }^B)$ saturates (i.e., the normalized LCE equals unity) as $ϵ$ increases, the spectral ordering method does not infer ${\text{boldmath}\sigma }^B$ better than random; it is deemed that the algorithm has reached the detectability limit.