Fast Approximation of Similarity Graphs
with Kernel Density Estimation

Given any set $X=\{x_{1},\ldots,x_{n}\}\subset\mathbb{R}^{d}$ and a kernel function $k:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}_{\geq 0}$, a fully connected similarity graph $\mathsf{K}=(V,E,w)$ of $X$ consists of $n$ vertices, and every $v_{i}\in V$ corresponds to $x_{i}\in X$; we set $w(v_{i},v_{j})\triangleq k(x_{i},x_{j})$ for any different $v_{i}$ and $v_{j}$. We introduce an efficient algorithm that constructs a sparse graph $\mathsf{G}$ directly from $X$, such that $\mathsf{K}$ and $\mathsf{G}$ share the same cluster-structure, and the graph matrices for $\mathsf{K}$ and $\mathsf{G}$ have approximately the same eigen-gap. This ensures that spectral clustering from $\mathsf{G}$ and $\mathsf{K}$ return approximately the same result.

The design of our algorithm is based on a novel reduction from the approximate construction of similarity graphs to the problem of Kernel Density Estimation (KDE). This reduction shows that any algorithm for the KDE can be employed to construct a sparse representation of a fully connected similarity graph, while preserving the cluster-structure of the input data points. This is summarised as follows:

This result builds a novel connection between the KDE and the fast construction of similarity graphs, and further represents a state-of-the-art algorithm for constructing similarity graphs. For instance, when employing the fast Gauss transform [10] as the KDE solver, Theorem 1 shows that a sparse representation of the fully connected similarity graph with the Gaussian kernel can be constructed in $\widetilde{O}(n)$ time when $d$ is constant. As such, in the case of low dimensions, spectral clustering runs as fast (up to a poly-logarithmic factor) as the time needed to read the input data points. Moreover, any improved algorithm for the KDE would result in a faster construction of approximate similarity graphs.

To demonstrate the significance of this work in practice, we compare the performance of our algorithm with five competing algorithms from the well-known scikit-learn library [20] and FAISS library [11]: the algorithm that constructs a fully connected Gaussian kernel graph, and four algorithms that construct different variants of $k$-nearest neighbour graphs. We apply spectral clustering on the six constructed similarity graphs, and compare the quality of the resulting clustering. For a typical input dataset of 15,000 points in $\mathbb{R}^{2}$, our algorithm runs in 4.7 seconds, in comparison with between 16.1 – 128.9 seconds for the five competing algorithms from scikit-learn and FAISS libraries. As shown in Figure 1, all the six algorithms return reasonable output.

We further compare the quality of the six algorithms on the BSDS image segmentation dataset [2], and our algorithm presents clear improvement over the other five algorithms based on the output’s average Rand Index. In particular, due to its quadratic memory requirement in the input size, one would need to reduce the resolution of every image down to 20,000 pixels in order to construct the fully connected similarity graph with scikit-learn on a typical laptop. In contrast, our algorithm is able to segment the full-resolution image with over 150,000 pixels. Our experimental result on the BSDS dataset is showcased in Figure 2 and demonstrates that, in comparison with SKLearn GK, our algorithm identifies a more detailed pattern on the butterfly’s wing. In contrast, none of the $k$-nearest neighbour based algorithms from the two libraries is able to identify the wings of the butterfly.

There is a number of work on efficient constructions of $\varepsilon$-neighbourhood graphs and $k$-NN graphs. For instance, Dong et al. [9] presents an algorithm for approximate $k$-NN graph construction, and their algorithm is based on local search. Wang et al. [32] presents an LSH-based algorithm for constructing an approximate $k$-NN graph, and employs several sampling and hashing techniques to reduce the computational and parallelisation cost. These two algorithms [9, 32] have shown to work very well in practice, but lack a theoretical guarantee on the performance.

Our work also relates to a large and growing number of KDE algorithms. Charikar and Siminelakis [7] study the KDE problem through LSH, and present a class of unbiased estimators for kernel density in high dimensions for a variety of commonly used kernels. Their work has been improved through the sketching technique [24], and a revised description and analysis of the original algorithm [4]. Charikar et al. [6] presents a data structure for the KDE problem, and their result essentially matches the query time and space complexity for most studied kernels in the literature. In addition, there are studies on designing efficient KDE algorithms based on interpolation of kernel density estimators [29], and coresets [12].

Our work further relates to efficient constructions of spectral sparsifiers for kernel graphs. Quanrud [22] studies smooth kernel functions, and shows that an explicit $(1+\varepsilon)$-approximate spectral approximation of the geometric graph with $\widetilde{O}(n/\varepsilon^{2})$ edges can be computed in $\widetilde{O}(n/\varepsilon^{2})$ time. Bakshi et al. [5] proves that, under the strong exponential time hypothesis, constructing an $O(1)$-approximate spectral sparsifier with $O(n^{2-\delta})$ edges for the Gaussian kernel graph requires $\Omega\left(n\cdot 2^{\log(1/\tau)^{0.32}}\right)$ time, where $\delta<0.01$ is a fixed universal constant and $\tau$ is the minimum entry of the kernel matrix. Compared with their results, we show that, when the similarity graph with the Gaussian kernel presents a well-defined structure of clusters, an approximate construction of this similarity graph can be constructed in nearly-linear time.

We use $X\triangleq\{x_{1},\ldots x_{n}\}$ to represent the set of input data points, where every $x_{i}\in\mathbb{R}^{d}$. Given $X$ and some kernel function $k:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}_{\geq 0}$, we use $\mathsf{K}=(V_{\mathsf{K}},E_{\mathsf{K}},w_{\mathsf{K}})$ to represent the fully connected similarity graph from $X$, which is constructed as follows: every $v_{i}\in V_{\mathsf{K}}$ corresponds to $x_{i}\in X$, and any pair of different $v_{i}$ and $v_{j}$ is connected by an edge with weight $w_{\mathsf{K}}(v_{i},v_{j})=k(x_{i},x_{j})$. One of the most common kernels used for clustering is the Gaussian kernel, which is defined by

for some bandwidth parameter $\sigma$. Other popular kernels include the Laplacian kernel and the exponential kernel which use $\|x_{i}-x_{j}\|_{1}$ and $\|x_{i}-x_{j}\|_{2}$ in the exponent respectively.

Our work is based on algorithms for kernel density estimation (KDE), which is defined as follows. Given a kernel function $k:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}_{\geq 0}$ with $n$ source points $x_{1},\ldots,x_{n}\in\mathbb{R}^{d}$ and $m$ target points $y_{1},\ldots,y_{m}\in\mathbb{R}^{d}$, the KDE problem is to compute $g_{[1,n]}(y_{1}),\ldots g_{[1,n]}(y_{m})$, where

for $1\leq i\leq m$. While a direct computation of the $m$ values $g_{[1,n]}(y_{1}),\ldots g_{[1,n]}(y_{m})$ requires $mn$ operations, there is substantial research to develop faster algorithms approximating these $m$ quantities.

In this paper we are interested in the algorithms that approximately compute $g_{[1,n]}(y_{i})$ for all $1\leq i\leq m$ up to a $(1\pm\epsilon)$-multiplicative error, and use $T_{\mathsf{KDE}}(m,n,\epsilon)$ to denote the asymptotic complexity of such a KDE algorithm. We also require that $T_{\mathsf{KDE}}(m,n,\epsilon)$ is superadditive in $m$ and $n$; that is, for $m=m_{1}+m_{2}$ and $n=n_{1}+n_{2}$, we have

it is known that such property holds for many KDE algorithms (e.g., [1, 6, 10]).

At a very high level, our designed algorithm applies a KDE algorithm as a black-box, and constructs a cluster-preserving sparsifier by simultaneous sampling of the edges from a non-explicitly constructed fully connected graph. To explain our technique, we first claim that, for an arbitrary $x_{i}$, a random $x_{j}$ can be sampled with probability $k(x_{i},x_{j})/\deg_{\mathsf{K}}(v_{i})$ through $O(\log n)$ queries to a KDE algorithm. To see this, notice that we can apply a KDE algorithm to compute the probability that the sampled neighbour is in some set $X_{1}\subset X$, i.e.,

where we use $g_{X}(y_{i})$ to denote that the KDE is taken with respect to the set of source points $X$. Based on this, we recursively split the set of possible neighbours in half and choose between the two subsets with the correct probability. The sampling procedure is summarised as follows, and is illustrated in Figure 3. We remark that the method of sampling a random neighbour of a vertex in $\mathsf{K}$ through KDE and binary search also appears in Backurs et al. [5]

1. 1.

   Set the feasible neighbours to be $X=\{x_{1},\ldots,x_{n}\}$.

2. 2.

   While $|X|>1$:

   • •

     Split $X$ into $X_{1}$ and $X_{2}$ with $|X_{1}|=\lfloor|X|/2\rfloor$ and $|X_{2}|=\lceil|X|/2\rceil$.

   • •

     Compute $g_{X}(x_{i})$ and $g_{X_{1}}(x_{i})$; set $X\leftarrow X_{1}$ with probability $g_{X_{1}}(x_{i})/g_{X}(x_{i})$, and $X\leftarrow X_{2}$ with probability $1-g_{X_{1}}(x_{i})/g_{X}(x_{i})$.

3. 3.

   Return the remaining element in $X$ as the sampled neighbour.

Next we generalise this idea and show that, instead of sampling a neighbour of one vertex at a time, a KDE algorithm allows us to sample neighbours of every vertex in the graph “simultaneously”. Our designed sampling procedure is formalised in Algorithm 1.

Finally, to construct a cluster-preserving sparsifier, we apply Algorithm 1 to sample $O(\log n)$ neighbours for every vertex $v_{i}$, and set the weight of every sampled edge $v_{i}\sim v_{j}$ as

where $\widehat{p}(i,j)\triangleq\widehat{p}_{i}(j)+\widehat{p}_{j}(i)-\widehat{p}_{i}(j)\cdot\widehat{p}_{j}(i)$ is an estimate of the sampling probability of edge $v_{i}\sim v_{j}$, and

for some constant $C\in\mathbb{R}^{+}$; see Algorithm 2 for the formal description.

Now we analyse Algorithm 2, and sketch the proof of Theorem 2. We first analyse the running time of Algorithm 2. Since it involves $O(\log n)$ executions of Algorithm 1 in total, it is sufficient to examine the running time of Algorithm 1.

We visualise the recursion of Algorithm 1 with respect to $S$ and $X$ in Figure 4. Notice that, although the number of nodes doubles at each level of the recursion tree, the total number of samples $S$ and data points $X$ remain constant for each level of the tree: it holds for any $i$ that $\bigcup_{j=1}^{2^{i}}S_{i,j}=S_{0}$ and $\bigcup_{j=1}^{2^{i}}X_{i,j}=X_{0}$. Since the running time of the KDE algorithm is superadditive, the combined running time of all nodes at level $i$ of the tree is

Hence, the total running time of Algorithm 2 is $\widetilde{O}(T_{\textsf{KDE}}(n,n,\epsilon))$ as the tree has $\log_{2}(n)$ levels. Moreover, when applying the Fast Gauss Transform (FGT) [10] as the KDE algorithm, the total running time of Algorithm 1 is $\widetilde{O}(n)$ when $d=O(1)$ and $\tau=\Omega(1/\mathrm{poly}(n))$.

Finally, we prove that the output of Algorithm 2 is a cluster preserving sparsifier of a fully connected similarity graph, and has $\widetilde{O}(n)$ edges. Notice that, comparing with the sampling probabilities $p_{u}(v)$ and $p_{v}(u)$ used in the SZ algorithm, Algorithm 2 samples each edge through $O(\log n)$ recursive executions of a KDE algorithm, each of which introduces a $(1+\epsilon)$-multiplicative error. We carefully examine these $(1+\epsilon)$-multiplicative errors and prove that the actual sampling probability $\widetilde{p}(i,j)$ used in Algorithm 2 is an $O(1)$-approximation of $p_{e}$ for every $e=\{v_{i},v_{j}\}$. As such the output of Algorithm 2 is a cluster-preserving sparsifier of a fully connected similarity graph, and satisfies the two properties of Theorem 2. The total number of edges in $\mathsf{G}$ follows from the sampling scheme. We refer the reader to the supplementary material for the complete proof of Theorem 2.

In this experiment we use the scikit-learn library to generate 15,000 data points in $\mathbb{R}^{2}$ from a variety of classical synthetic clustering datasets. The data is generated with the make_circles, make_moons, and make_blobs methods of the scikit-learn library with a noise parameter of $0.05$. The linear transformation $\{\{0.6,-0.6\},\{-0.4,0.8\}\}$ is applied to the blobs data to create asymmetric clusters. As shown in Figure 1, all of the algorithms return approximately the same clustering, but our algorithm runs much faster than the ones from scikit-learn and FAISS.

We further compare the speedup of our algorithm against the five competitors on the two_moons dataset with an increasing number of data points, and our result is reported in Figure 5. Notice that the running time of the scikit-learn and FAISS algorithms scales roughly quadratically with the size of the input, while the running time of our algorithm is nearly linear. Furthermore, we note that on a laptop with 32 GB RAM, the SKLearn GK algorithm could not scale beyond 20,000 data points due to its quadratic memory requirement, while our new algorithm can cluster 1,000,000 data points in 240 seconds under the same memory constraint.

Finally we study the application of spectral clustering for image segmentation on the BSDS dataset. The dataset contains $500$ images with several ground-truth segmentations for each image. Given an input image, we consider each pixel to be a point $(r,g,b,x,y)^{\intercal}\in\mathbb{R}^{5}$ where $r$, $g$, $b$ are the red, green, blue pixel values and $x$, $y$ are the coordinates of the pixel in the image. Then, we construct a similarity graph based on these points, and apply spectral clustering, for which we set the number of clusters to be the median number of clusters in the provided ground truth segmentations. Our experimental result is reported in Table 1, where the “Downsampled” version is employed to reduce the resolution of the image to 20,000 pixels in order to run the SKLearn GK and the “Full Resolution” one is to apply the original image of over 150,000 pixels as input. This set of experiments demonstrates our algorithm produces better clustering results with repsect to the average Rand Index [23].