Non-Graph Data Clustering via $\mathcal{O}(n)$ Bipartite Graph Convolution

Let $\mathcal{G}=(\mathcal{V},\mathcal{E},\mathcal{W})$ be a graph where $\mathcal{V}$, $\mathcal{E}$, and $\mathcal{W}$ represent nodes, edges, and weights respectively. For an unweighted matrix, $\mathcal{W}$ can be ignored since the weights of edges are all regarded as 1. All matrices and vectors are represented by upper-case letters and lower-case letters in boldface, respectively. $|\cdot|$ represents the size of a set. $(\cdot)_{+}=\max(\cdot,0)$. For a real number $x$, $\lfloor x\rfloor$ represents the maximum integer that is not larger than $x$. $\textbf{1}_{n}$ denotes an $n$-dimension vector whose all entries equal 1. In computer science, the adjacency matrix, which is denoted by $\bm{A}$ in this paper, is a common way to represent the edges of a graph. Specifically, for an unweighted graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$, $A_{ij}=1$ if $\langle v_{i},v_{j}\rangle\in\mathcal{E}$ and $A_{ij}=0$ otherwise. A popular model of GNN is the graph convolution network (GCN) [28] where the graph convolution layer is accordingly formulated as [28]

where $\varphi(\cdot)$ is an activation function, $\bm{X}\in\mathbb{R}^{n\times d}$ is the representation of data points, $\bm{W}\in\mathbb{R}^{d\times d^{\prime}}$ denotes the coefficients to learn, $\bm{L}=\bm{D}^{-\frac{1}{2}}(\bm{A}+\bm{I})\bm{D}^{-\frac{1}{2}}$ is often viewed as a variant of the normalized Laplacian matrix, and $\bm{D}$ is the degree matrix of $(\bm{A}+\bm{I})$. Specifically, $\bm{D}$ is a diagonal matrix where $D_{ii}=\sum_{j=1}^{n}(A_{ij}+1)$.

In recent years, how to effectively extract features of graphs has attracted a lot of attention [28, 31]. Inspired by CNN, a few works [27, 38, 28] attempt to investigate convolution operation on graph-type data. With the rise of the attention mechanism, GAT [29, 39, 40] introduces self-attention into graph neural networks. Meanwhile, GNN can be applied to sequential tasks [41, 30], which is known as spatial-temporal graph neural network (STGNN). Some theoretical works [42, 43, 44, 45, 46, 47] are also investigated in recent years. Some works [46] argue the impact of over-smoothing while literature [47] proposes the over-squashing. To extend GCN [28] into unsupervised learning, graph auto-encoders [48, 49, 50, 51, 52], based on auto-encoder [53, 54], are developed. Different from traditional auto-encoders, GCN [28] is used as an encoder and the decoder calculates inner-products of embedding to reconstruct graphs. To build a symmetric architecture, GALA [50] designs a sharpening operator since the graph convolution is equivalent to a smooth operator. SDCN [37] has tried to utilize GNN for clustering, while it just constructs a KNN graph as the preprocessing and suffers from the inefficiency problem as well. Although GNN has achieved outstanding performance [28, 55], it requires too much memory and time. Some methods [33, 34, 36, 35, 56] for accelerating training and reducing consumption of memory are proposed. For instance, SGC [35] attempts to accelerate by removing non-linear activations in GCN, which reduces the number of parameters. More importantly, it compresses multiple layers into one layer such that stochastic gradient descent can be used. GraphSAGE [33] attempts to sample part of nodes for each node to reduce the computational complexity, while FastGCN [32] samples nodes independently. ClusterGCN [34] tries to find an approximate graph with multiple connected components via METIS [57] to employ batch gradient descent.

In traditional clustering, $k$-means [58] works efficiently on large scale data but it can only deal with spherical data, which limits the application of $k$-means. By virtue of deep learning, deep clustering [22, 17, 23, 18, 19, 59] often aims to learn expressive representations and then feed the representations into the simple clustering module. The whole optimization is based on batch gradient descent so that they are still available on large scale data. However, they frequently suffers from the over-fitting and the performance heavily relies on the tuning of neural networks. For graph-based clustering methods [3, 2, 1], both computational complexity and space complexity of graph construction are $\mathcal{O}(n^{2})$. In particular, solving spectral clustering needs eigenvalue decomposition, which is also time-consuming. To apply spectral clustering on large scale data, a lot of variants [12, 9, 10, 11] have been developed. Anchor-based methods [12, 13, 14] are important extensions of spectral clustering. The core idea is to construct a graph implicitly via some representative points, which are usually named as anchors. By constructing a bipartite graph between samples and anchors, they only need to perform singular value decomposition on the Laplacian of the bipartite graph [13]. The anchor-based spectral clustering provides an elegant scheme to improve the efficiency of graph-based models.

Before elaborating on how to efficiently improve the capacity of graph-based clustering via GNNs on non-graph data, we introduce the major hurdles respectively.

No Prior Graph: In the scenarios of clustering, only features are given as the prior information. The graph is usually constructed as the preprocessing in graph-based clustering methods. The key of performance is the quality of the constructed graph. When the scattering of data is too complicated to precisely capture the key information based on the original data (which is the main bottleneck of the most existing graph-based clustering methods), it is a rational scheme to map the data into a latent space and then build a more precise graph in this space. As we focus on graph-based models, GNNs are more rational schemes compared with classical neural networks. Moreover, the utilization of the estimated graph (via GNNs) to update graph leads to a self-supervised procedure. The graph construction usually results in $\mathcal{O}(n^{2})$ costs, which also causes the inefficiency of graph-based clustering methods. The scheme to solve this problem will be elaborated in Section 3.2.

Inefficiency of GNNs: In the forward propagation of GNNs, it usually require neighbored nodes to compute the hidden feature of one node. This computation results in the inefficiency of GNNs on large graph [60], especially with the increase of depth. In this paper, we simply focus on GCN, the widely used variant of GNN, since the theoretical analysis does not depend on specific GNN modules. For instance, it requires $\mathcal{O}(|\mathcal{E}|d)$ time to compute $\bm{L}\bm{X}$ in the forward propagation. If $\mathcal{G}$ is sparse enough, the computational complexity is acceptable. If there are plenty of edges in $\mathcal{G}$, the computational complexity can be roughly regarded as $\mathcal{O}(n^{2}d)$, so that GCNs can not be applied in this case. We will show how to alleviate it into $\mathcal{O}(n)$ in Sections 3.3 and 3.4.

Inefficiency of Graph-Based Clustering Methods: Even after the $\mathcal{O}(n^{2})$ graph construction, the computation complexity reaches $\mathcal{O}(n^{3})$ due to the eigenvalue decomposition. The problem is easy to address after introducing the $\mathcal{O}(n)$ bipartite graph convolution and the details can be found in Section 3.6.

By virtue of the generative perspective introduced in this section, the probability transition matrix is accordingly built so that the Markov process can be used to implement the efficient bipartite graph convolution, which is elaborated in Section 3.3. The generative perspective also provides a training goal, reconstructing the distribution, for AnchorGAE.

As the goal is to capture the structure information, a multi-layer GCN is a rational scheme to implement $f(\bm{x}_{i})$ which is defined in Eq. (2). In this subsection, we will show why the bipartite graph convolution accelerates GCN.

We utilize the bipartite graph, defined in Eq. (12), to accelerate the graph convolution. In virtue of the Markov process [61], $p(v_{j}|v_{i})$ can be obtained by the one-step transition probability, i.e.,

Similarly, $p(u_{j}|u_{i})=\sum_{l=1}^{n}p(u_{j}|v_{l})p(v_{l}|u_{i})$. Formally, the constructed graphs of $\mathcal{V}$ and $\mathcal{U}$ are defined as

Accordingly, $\bm{A}$ is the weighted adjacency matrix, constructed by anchors, of $n$ samples, and $\bm{A}_{t}$ is the weighted adjacency of $m$ anchors.

Therefore, $\bm{A}$ can be directly used as a part of the graph convolution operator. Since

the degree matrix of $\bm{A}$ is the identity matrix, i.e., $\bm{L}=\bm{A}$. Therefore, the output of GCN is formulated as

Here, we explain why the above equation accelerates the operation. It should be emphasized that $\bm{A}$ is not required to be computed explicitly. The core is to rearrange the order of matrix multiplications as follows

According to the above order of computation, the computational complexity is reduced to $\mathcal{O}(nmd+m^{2}d+mdd^{\prime}+nmd^{\prime})$. Since $d$, $d^{\prime}$, and $m$ are usually much smaller than $n$, the complexity can be regarded as $\mathcal{O}(nm)$. Moreover, if the amount of anchors is small enough, the required time of each forward propagation of a GCN layer is $\mathcal{O}(n)$.

After showing the accelerated GCN mapping based on anchors and the bipartite graph, we discuss more details of the whole procedure in this section, especially the changes of architecture brought by anchors.

As the task is unsupervised, the whole architecture is based on graph auto-encoder (GAE). The encoder consists of multiple bipartite GCN layers defined in Eq. (17),

where $L$ denotes the number of layers. The embedding (denoted by $\bm{Z}$), which is produced by the encoder, is a non-linear mapping for the raw representation of $\mathcal{V}$.

However, the mentioned encoder just transforms $\mathcal{V}$ and there is no edge between $\mathcal{V}$ and $\mathcal{U}$. In other words, $\mathcal{U}$ can not be projected via the network with the adjacency $\bm{A}$. Meanwhile, GCN layers suffer from the out-of-sample problem. Accordingly, $\mathcal{U}$ can not be directly projected into the deep feature space by the encoder. To map anchors into the same feature space, a siamese [62] encoder is thus designed. The encoder for anchors should share the same parameters but has its own graph structure. Surprisingly, the one-step probability transition matrix defined in Eq. (15) gives us a graph that is only composed of anchors, $\bm{A}_{t}=\bm{\Delta}^{-1}\bm{B}^{T}\bm{B}$. As

the corresponding degree matrix, $\bm{D}_{t}$, is also the identity matrix, i.e.,

Accordingly, the embedding of anchors is formulated as

where $\bm{C}\in\mathbb{R}^{m\times d}$ represents $m$ anchors.

Loss: As the underlying connectivity distributions have been estimated, we use them, rather than the adjacency in the existing GAEs, as the target for the decoder to rebuild. According to Assumption 1, the connectivity probability should be reconstructed according to Euclidean distances between $\mathcal{V}$ and $\mathcal{U}$. Formally speaking, the decoder is

Instead of using mean-square error (MSE) to reconstruct $A$, AnchorGAE intends to reconstruct the underlying connectivity distribution and employ the cross-entropy to measure the divergence,

The GCN part of AnchorGAE is trained via minimizing $\mathcal{L}$.

So far, how to get anchors at the beginning is ignored. In deep learning, a common manner is to pre-train a neural network. Nevertheless, the pre-training of raw graph auto-encoders will result in nearly $\mathcal{O}(n^{2})$ time, which violates the motivation of our model. Therefore, we first compute anchors and transition probability on the raw features. Then anchors and transition probability will be rectified dynamically, which will be elaborated in Section 3.5.

Since the initialized anchors and transition probability is calculated from the original features, the information hidden in $\bm{A}$ may be limited especially when the data distribution is complicated. To precisely exploit the high-level information in a self-supervised way, a feasible method is to adaptively update anchors according to the learned embedding.

To achieve this goal, we first need to recompute anchors, $\mathcal{U}$. According to Eq. (7) and (10), we can obtain anchors under the deep representations, $\{f(\bm{c}_{j})\}_{j=1}^{m}$. However, the input of our model requires anchors under the raw features. In other words, it is inevitable to compute $\bm{c}_{j}=f^{-1}(f(\bm{c}_{j}))$. Unfortunately, it is hard to ensure $f(\cdot)$ an invertible mapping. A scheme is to utilize the classical auto-encoders since the decoder that tries to reconstruct the raw input from the embedding can be regarded as an approximation of the inverse mapping of the encoder. To simplify the discussion in this paper, we simply set

to estimate $\bm{c}_{j}$ under the original feature space.

In our expectations, the self-supervised update of anchors and transition probability would become more and more precise. However, the performance of AnchorGAE will become worse and worse after several updates of anchors, which has been proved theoretically and empirically. An interesting phenomenon is that a better reconstruction indicates a more severe collapse. The following theorem provides a rigorous description of this collapse.

The above conclusion means that $\hat{p}(\cdot|v_{i})$ is a sparse and uniform distribution. In other words, the constructed graph degenerates into an unweighted graph, which is usually unexpected in clustering. More intuitively, the degeneration is caused by projecting samples from the identical cluster into multiple tiny and cohesive clusters as shown in Figure 2. These tiny clusters have no connection with each other and they scatter randomly in the representation space, such that the clustering result is disturbed. To avoid this collapse, we propose two strategies: (1) increase $k$ which represents the sparsity of neighbors; (2) decrease the number of anchors, $m$. In this paper, we dynamically increase $k$, which is equivalent to attempt to connect these tiny groups before they scatter randomly, as follows

Suppose that there are $c$ clusters and the number of samples in the smallest cluster is represented by $n_{s}$. Then we can define the upper-bound of $k$ and increment of sparsity $\Delta k$ as

where $k_{0}$ denotes the sparsity of the initialized anchor graph and $E$ be the number of iterations to update anchors. If we have no prior information of $n_{s}$, $n_{s}$ can be simply set as $\lfloor n/c\rfloor$ or $\lfloor n/(2c)\rfloor$. The brief modification helps a lot to guarantee the stability and performance of AnchorGAE.

When AnchorGAE is trained, we can perform some fast clustering algorithms on $\bm{B}$ or $\bm{Z}$, which is elaborated in Section 3.6. The optimization of AnchorGAE is summarized in Algorithm 2.

After training AnchorGAE, we can simply run $k$-means on the learned embedding. Due to Assumption 1, Euclidean distances of similar points will be small so that it is appropriate to use $k$-means on the deep representations.

Another method is to perform spectral clustering on the constructed graph, which is the same as anchor-based spectral clustering. Let

be the normalized Laplacian matrix. It should be emphasized that $\bm{\mathcal{L}}$ is also the unnormalized Laplacian matrix. Accordingly, the objective of spectral clustering is

where $\bm{F}\in\mathbb{R}^{n\times c}$ is the soft indicator matrix. Let $\hat{\bm{B}}=\bm{B}\bm{\Delta}^{-\frac{1}{2}}$. Then the eigenvalue decomposition of $\bm{A}$ can be transformed into the singular value decomposition of $\hat{{\bm{B}}}$, which reduces the complexity from $\mathcal{O}(n^{2}c)$ to $\mathcal{O}(mnc)$.

Finally, one can perform clustering on the bipartite graph directly rather than the constructed graph. Let the bipartite graph be

The degree matrix is represented as

Then the normalized cut used spectral clustering is to optimize

where

Since $\bm{D}_{v}=\bm{I}$, the above problem can be further reformulated as

The closed-form solution can be calculated by

where $\tilde{\bm{V}}$ and $\tilde{\bm{U}}$ are $c$ leading left and right singular vectors of $\bm{B}\bm{D}_{u}^{-\frac{1}{2}}$. The details can be found in [65]. In our experiments, we report the results obtained by the above technique based on bipartite graphs, since it usually outperforms the simple $k$-means.

The performance of the AnchorGAE is evaluated on 6 non-graph datasets, including 2 UCI datasets (ISOLET [66] and SEGMENT [66]) and 4 image datasets (USPS [67], MNIST-test [68], MNIST-full [68], and Fashion-MNIST [69]). Remark that the two UCI datasets show the difference between AnchorGAE and CNN-based clustering models. Since each data point can only be represented by a vector, the CNN-based model can not be applied to these datasets. Note that Fashion-MNIST is denoted by Fashion in Table II The details of these datasets are shown in Table II.

AnchorGAE is compared with 9 methods, including 5 anchor-based methods (Nystrom [9], CSC [10], KASP [11], LSC [13], and SNC [14]), 3 GAE-based methods (SGC [35], GAE [48] and ClusterGCN [34]), 2 deep methods (DEC [17] and SpectralNet [25]), and classical K-Means [58]. All codes are downloaded from the homepages of authors.

There are three clustering evaluation criteria used in our experiments, including the clustering accuracy (ACC), normalized mutual information (NMI), and running time. In our experiments, the number of anchor points is set from 100 to 1000 and the increasing step is 100. The best results are reported in Table I. The activation function of the last layer is set to linear, while the other layers use the ReLU function. Particularly, ClusterGCN recommends 4 layers convolution. To ensure the convolution structure of ClusterGCN the same as others, we modify the 4 layers of convolution structure (denoted by ClusterGCN-L4) to 2 layers (denoted by ClusterGCN-L2) and remain 4 layers of the convolution structure as a competitor. To apply ClusterGCN to unsupervised tasks, we use the same GAE loss and decoder architecture that is same as the other models. To test the scalability of AnchorGAE in a deeper network, we design 2 layers (denoted by AnchorGAE-L2) and 4 layers (denoted by AnchorGAE-L4) convolution structure separately. The 2 layers neurons are designed as 128 and 64 or 256 and 32. For the 4 layers neurons, we set the first three layers as 16 and set the last layer as $L$+1. $L$ is the number of labels. The maximum iterations to update network are 200, while the maximum epochs to update anchors are set as 5. The initial sparsity $k_{0}$ is set as 3.

To investigate the impact of anchors, we run AnchorGAE with different numbers of anchors and initial sparsity, i.e., $m$ and $k_{0}$. We also conduct 3 ablation experiments based on AnchorGAE-L2 to study the impact of different parts. Firstly, we fix anchors and transition probability during training, which is denoted by Ours-A. Secondly, we fix the sparsity $k$, which is denoted by Ours-B. Thirdly, we use a KNN graph, instead of a weighted graph, with the adaptive update and incremental $k$ in our model, which is denoted by Ours-C. All methods are run 10 times and the average results are reported. All experiments are conducted on a PC with an NVIDIA GeForce GTX 1660 GPU SUPER.

ACC and NMI are shown in Table I while the runtime is reported in Figure 5. For ACC and NMI, the best results are bolded and suboptimal results are underlined. Especially, some methods either spend too much memory to run on MNIST-full and Fashion-MNIST, or the codes provided by authors fail to run on them. So they are denoted by a horizontal line in Table I. Similarly, they are represented as $\infty$ in Figure 5. Besides, the visualization of SEGMENT is shown in Figure 3. Figure 4 shows the impact of the number of anchors while the effect of the initial sparsity is illustrated in Figure 6.

From the experimental results, we can conclude that:

1. 1.

   On almost all datasets, AnchorGAE-L2 obtains the best results in all metrics. On MNIST-full, AnchorGAE obtains the suboptimal NMI, while DEC, which consists of much more layers, outperforms AnchorGAE. Meanwhile, AnchorGAE, with only 2 layers, achieves better ACC than DEC. AnchorGAE needs to compute anchors such that it requires more time.

2. 2.

   On MNIST-full, normal GAE-based methods fail to run due to out-of-memory since these methods need to save the graph structure in each convolution. AnchorGAE can convolute without constructing a graph explicitly, which not only needs less memory but also reduces runtime.

3. 3.

   When sparsity $k$ is fixed or anchors and transition probability (i.e., matrix $\bm{B}$) is not updated, AnchorGAE usually gets poor results on many datasets. For instance, Ours-B shrinks about 55% in ACC and 55% in NMI on MNIST-full, while Ours-A reduces about 40% in ACC and reduces about 15% in NMI.

4. 4.

   To ensure fair competition and distinguish between anchor-based methods and GAE-based methods, we set the same number of anchors for AnchorGAE-L4 and other anchor-based methods. Note that the first two layers of AnchorGAE-L4 are the same as AnchorGAE-L2. Through observing the results of experiments, AnchorGAE-L4 obtains acceptable performance in most of the datasets but not the best results. It may be caused by the too small dimension of the final embedding.

5. 5.

   If the initial sparsity, $k_{0}$, is set as a too large number, then it results in dense connections. As a result, a great quantity of wrong information is captured.

6. 6.

   From the comparison between AnchorGAE and Ours-C, we find that the weighted graphs are still important for clustering even after introducing non-linear GNN mapping modules.