Accelerated Spectral Clustering Using Graph Filtering Of Random Signals

Introduction

Spectral clustering has become a popular clustering algorithm, due to its simplicity of implementation and high performance for many different types of datasets . Given a set of $N$ data points $(\bm{x}_{1},\bm{x}_{2},\cdots,\bm{x}_{N})$, it basically transforms them non-linearly into a $k$-dimensional space first, by computing a similarity matrix $\bm{W}$ from the data, and second by computing the $k$ first eigenvectors of its Laplacian. Calculating these eigenvectors is the computational bottleneck of spectral clustering: it becomes prohibitive when $N$ becomes large and/or when the $k$-th eigengap becomes too small . Circumventing this issue is an active area of research . Another difficulty, common to all clustering methods, is the estimation of the usually unknown number of clusters $k$ .

We propose a new method that jointly avoids to partially diagonalize the Laplacian and proposes a stability measure to estimate $k$. This method is based on the emerging field of graph signal processing , where the graph we consider here is defined by the weighted adjacency matrix $\bm{W}$. In our previous work , we proposed the use of the graph wavelet or scaling function transforms of random vectors to detect multiscale communities in networks. In this paper, we build upon this idea of using filtered random vectors as probes of the underlying graph’s structure; and further prove that ideal low-pass graph filters have deep connections with spectral clustering. Taking advantage of the fast graph filtering defined in to low-pass filter such random signals without computing the first $k$ eigenvectors of the graph’s Laplacian, we propose an accelerated spectral clustering method that has the collateral advantage of defining a stability measure to estimate the number of clusters $k$.

In Section 2, we recall the graph signal processing notations and tools we will use in the paper, as well as the classical spectral clustering algorithm. In Section 3, we prove that one can filter random signals to estimate the spectral clustering distance matrix. In Section 4, we detail our algorithm and propose a stability measure to estimate the number of clusters $k$. We finally show results obtained on a controled dataset in Section 5; before concluding in Section 6.

Background

Let $\mathcal{G}=(\mathcal{V},\mathcal{E},\mathbf{W})$ be an undirected weighted graph with $\mathcal{V}$ the set of $N$ nodes, $\mathcal{E}$ the set of edges, and $\mathbf{W}$ the weighted adjacency matrix such that $W_{ij}=W_{ji}\geq 0$ is the weight of the edge between nodes $i$ and $j$.

Consider the graph’s combinatorial LaplacianOne could use other types of Laplacians, such as the normalized Laplacian; which adds a normalization step after step 2 of the spectral clustering algorithm (see Section 2.2) and slightly changes our subsequent proofs. matrix $\bm{L}=\mathbf{S}-\mathbf{W}$ where $\mathbf{S}$ is diagonal with $S_{ii}=s_{i}=\sum_{j\neq i}W_{ij}$ the strength of node $i$. $\bm{L}$ is real symmetric, therefore diagonalizable in an orthogonal basis: its spectrum is composed of its set of eigenvalues $\{\lambda_{l}\}_{l=1\dots N}$ that we sort: $0=\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}\leq\dots\leq\lambda_{N}$; and of $\mathbf{\bm{\chi}}$ the orthonormal matrix of its eigenvectors: $\bm{\chi}=\left(\bm{\chi}_{1}|\bm{\chi}_{2}|\dots|\bm{\chi}_{N}\right)$. Considering only connected graphs, the multiplicity of eigenvalue $\lambda_{1}=0$ is one . By analogy to the continuous Laplacian operator whose eigenfunctions are the classical Fourier modes and eigenvalues their squared frequencies, the columns of $\bm{\chi}$ are considered as the graph’s Fourier modes, and $\{\sqrt{\lambda_{l}}\}_{l}$ as its set of associated “frequencies” . Other types of graph Fourier matrices have been proposed (e.g. ), but in order to exhibit the link between graph signal processing and classical spectral clustering (that partially diagonalizes the Laplacian matrix), the Laplacian-based Fourier matrix is more natural.

2 Spectral clustering

Let us recall the method of spectral clustering . The input is the set of data points $(\bm{x}_{1},\bm{x}_{2},\cdots,\bm{x}_{N})$ and $k$ the number of clusters one desires. Follow the steps:

Compute the pairwise similarities $s(\bm{x}_{i},\bm{x}_{j})$ and create a similaritySee for several choices of similarity measure $s$ as well as several ways to create $\bm{W}$ from the $s(\bm{x}_{i},\bm{x}_{j})$. graph $\bm{W}$. Compute its Laplacian $\bm{L}$.

where $\delta_{i}(j)=1$ if $j=i$ and elsewhere.

Run $k$-means (or any clustering algorithm) with the Euclidean distance $D_{ij}=||\bm{f}_{i}-\bm{f}_{j}||$ to obtain $k$ clusters.

3 Graph filtering

and $\bm{H}_{\lambda_{c}}$ its associated graph filter operator.

4 Fast graph filtering

In order to filter a signal by $h$ without diagonalizing the Laplacian matrix, one may rely on a polynomial approximation of order $m$ of $h$ on $[0,\lambda_{N}]$: $\exists\{\alpha_{l}\}_{l\in[0,m]}\mbox{ s.t. }h(\lambda)\simeq\sum_{l=0}^{m}\alpha_{l}\lambda^{l}$. This enables us to approximate $\bm{Hx}$:

where we note that the approximation only requires matrix-vector multiplication and effectively bypasses the diagonalisation of the Laplacian; and where we introduce the notation $\mathcal{F}^{m}_{h}$ to design the fast filtering operator of order $m$ associated to $h$. This fast filtering method has a total complexity of $O(m(|\mathcal{E}|+N))$ , which in the case of sparse graph where $|\mathcal{E}|\sim N$, ends up being $O(mN)$; compared to the $O(N^{3})$ complexity needed to diagonalize the Laplacian matrix.

In this study, we only consider ideal low-pass $h_{\lambda_{c}}$ defined on $[O,\lambda_{N}]$. Notice that $h_{\lambda_{c}/\lambda_{N}}=h_{\lambda_{c}}(\lambda\times\lambda_{N})$ is always defined on $$and does not depend on any graph anymore. One can therefore beforehand tabulate$m^{*}$as a function of its sole parameter$\lambda_{c}/{\lambda_{N}}$(and$\delta$). Then, given a real filter$h_{\lambda_{c}}$defined on$[O,\lambda_{N}]$to approximate, one only needs to refer to this table and choose$m=m^{*}(\lambda_{c}/\lambda_{N};\delta)$. In the following, we fix$\delta=0.1$.

Accelerating spectral clustering

We show that we can estimate the distance $D_{ij}=||\bm{f}_{i}-\bm{f}_{j}||$ by filtering a few random signals with the ideal low-pass $h_{\lambda_{k}}$ (as defined in Eq. (2) with $\lambda_{c}=\lambda_{k}$). First of all, consider the matrix operator $\bm{H}_{\lambda_{k}}$ associated to $h_{\lambda_{k}}$. One may write:

where $\bm{I}_{k}$ is the identity of size $k$, and $\bm{0}$ null block matrices.

Proposition: Let $\epsilon,\beta>0$ be given. If $\eta$ is larger than:

then with probability at least $1-N^{-\beta}$, we have: $\forall(i,j)\in[1,N]^{2}$

As i) the columns of $\bm{U}$ are normalized to 1 and orthogonal to each other, and ii) $\bm{R}$ is a random Gaussian matrix with mean zero and variance $1/\eta$, then $\bm{R}^{\prime}=\bm{R}^{\top}\bm{U}$ is also Gaussian with same mean and variance; and Equation (9) reads:

This enables us to apply Theorem 1.1 of (an instance of the Johnson-Lindenstrauss lemma) and finish the proof. ∎

Consequence: Setting $\beta$ to 1, and therefore the failure probability of Equation (8) to $1/N$, we only need to filter $\eta\gtrsim\frac{12}{\epsilon^{2}}\log{N}$ random signals to estimate (up to an error $\epsilon$) the spectral clustering distance matrix. How this error $\epsilon$ on the distance estimation theoretically affects the performance of the spectral clustering algorithm is still, to our knowledge, an open question. We observe experimentally (see Sec. 5) that using a number $\eta\gtrsim k$ (i.e. allowing a relatively high error $\epsilon^{2}\simeq 12\frac{\log{N}}{k}$) is usually enough for satisfying performance.

2 Fast filtering of random signals

In practice, we do not exactly filter these random signals by $h_{\lambda_{k}}$ as the computation of $\bm{H}_{\lambda_{k}}$ requires the diagonalisation of the Laplacian, which is precisely what we are trying to avoid. Instead, we take advantage of the fast filtering scheme recalled in Section 2.4. Still, one question remains: the fast filtering is based on the polynomial approximation of $h_{\lambda_{k}}$, which is itself parametrized by $\lambda_{k}$. Unless we compute the first $k$ eigenvectors of $\bm{L}$, thereby loosing our efficiency edge on other methods, we cannot know exactly the value of $\lambda_{k}$.

Accelerated spectral clustering

Consider a set of data points $(\bm{x_{1}},\bm{x_{2}},\cdots,\bm{x_{N}})$ and $k$ the number of desired clusters. The first step of the algorithm does not change as compared to Section 2.2: compute the pairwise similarities $s_{ij}$, create a similarity graph $\bm{W}$, and compute its Laplacian $\bm{L}$. Then:

Estimate $\bm{L}$’s largest eigenvalue $\lambda_{N}$ (necessary for steps 2 and 3).

2 Complexity considerations

We compare the time complexity of this algorithm and the classical spectral clustering algorithm. Let us separate our algorithm (resp. the classical algorithm) into two parts: the spectral estimation part consisting of steps 1 and 2 (resp. step 2) and the clustering part consisting in steps 4 to 6 (resp. steps 3 and 4). According to , $k$-means has complexity $O(\eta N)$ where $\eta$ is the dimension of the $N$ considered vectors. Considering the discussion of Section 2.4, the time complexity of the clustering part is thus $O(\eta(m+1)N)$ (resp. $O(kN)$). The time complexity of the spectral estimation part is difficult to estimate, as it depends on the graph-dependent eigengap $\lambda_{k+1}-\lambda_{k}$. In both algorithms, the larger is this eigengap, the faster is the convergence. Nevertheless, empirical observations show that, even though the clustering part of the classical algorithm is faster than ours; our algorithm makes up to that difference by computing even faster (especially for large $N$) its spectral estimation part than the classical one.

3 Estimate the number of clusters k𝑘k

In real data, the number of clusters $k$ to find is usually unknown. Instead, one has access to a (possibly large) interval of values $[k_{min},k_{max}]$. Notions of stability to estimate $k$ have been used in various contexts ; and we propose here a new one that naturally comes from the random vectors’ stochasticity. For all $k\in[k_{min},k_{max}]$, we perform our algorithm $J$ times using $J$ different realisations of the $\eta$ random signals, to obtain $J$ different clusterings $\{\mathcal{C}^{j}_{k}\}_{j\in[1,J]}$. We define a stability measure as the mean of the Adjusted Rand Index similarity between all pairs of clusterings :

Denote $k^{*}$ the value of $k$ for which $\gamma$ reaches its global maximum: we consider it as the relevant number of clusters.

Results

From these N points, we create a similarity matrix by building a $K$ nearest neighbours graph with $K~{}\sim\log{N}$ as suggested in as a classical way of generating sparse similarity graphs. Other wiring possibilities exist to create such a graph (see ) but we only show this particular one as the choice of construction does not affect our results (but still considering graphs with same sparsity).

We use Matlab’s $k$-means function (with 20 replicates) and the GSP toolbox for steps 2 and 3 of our algorithm.

First, we remark in this dataset that values of $\lambda_{k=10}/\lambda_{N}$ are small (of the order of $10^{-4}$), which necessitates a large $m\sim 200$ for a correct approximation of the ideal low-pass on $[0,\lambda_{N}]$ (see Section 2.4).

Let us now illustrate in Fig. 2 the stability measure $\gamma$ obtained using our method with $N=5000$, $\eta=2k$ and $J=20$: the global maximum correctly detects $k=10$ and Fig. 1 (right) shows one recovered labeling for $k=10$.

We compare our method vs. classical spectral clustering in Figure 3 in terms of performance and time of computation. We note that for $\eta=2k$ and $\eta=3k$, our method performs as well as the classical algorithm; while for $\eta=k$, on the other hand, the number of random signals becomes insufficient as we observe the recovery starting to fail. Up to $N\simeq 4.10^{4}$, the computation time is slightly faster with the classical algorithm. But as $N$ increases, the classical algorithm’s computing time increases significantly faster than our proposition’s: for $N=10^{5}$ for instance, computation time is 2 (resp. 2.5, 3) times faster when one uses $\eta=3k$ (resp. $2k$, $k$) random signals than the classical algorithm.

Conclusion

We propose a new method that paves the way to alternative spectral clustering methods bypassing the usual computational bottleneck of extracting the Laplacian’s first $k$ eigenvectors. We take advantage of the fast graph low-pass graph filtering of a few random vectors to estimate the spectral clustering distance. The use of random vectors makes our algorithm stochastic, which in turn enables us to define a stability measure $\gamma$ for any $k$: scales of interest maximize $\gamma$. Results on synthetic data show that our method is scalable and for $\eta\gtrsim k$, one has the same performance as with the classical spectral algorithm while reducing the time complexity by a few factors.

We prooved that the error on the estimation of the distance $D_{ij}$ is well controlled, but the question of how such an error propagates on the estimation of the clusters themselves is open. Moreover, the impact of the error $\delta$ of the polynomial approximation on the rest of the algorithm is still largely unknown and matter of future work.

References