Fast Approximation of Similarity Graphs with Kernel Density Estimation
Authors: Peter Macgregor
He Sun
Venue: NeurIPS 2023 spotlight (2023)
Keywords: similarity graphs
spectral clustering

Abstract: Constructing a similarity graph from a set $X$ of data points in $ \mathbb{R}^d$ is the first step of many modern clustering algorithms. However, typical constructions of a similarity graph have high time complexity, and a quadratic space dependency with respect to $|X|$. We address this limitation and present a new algorithmic framework that constructs a sparse approximation of the fully connected similarity graph while preserving its cluster structure. Our presented algorithm is based on the kernel density estimation problem, and is applicable for arbitrary kernel functions. We compare our designed algorithm with the  well-known implementations from the scikit-learn library and the FAISS library,  and find that our method significantly outperforms the implementation from both libraries on a variety of datasets.

Fast Approximation of Similarity Graphs with Kernel Density Estimation
Peter Macgregor School of Informatics University of Edinburgh United Kingdom
He Sun School of Informatics
University of Edinburgh United Kingdom
Abstract
Constructing similarity graph from set X of data points in Rd is the first step of many modern clustering algorithms. However, typical constructions a similarity graph have high time complexity, and quadratic space dependency with respect to [X]. We address this limitation and present a new algorithmic framework that constructs a sparse approximation of the fully connected similarity graph while preserving its cluster structure. Our presented algorithm is based on the kernel density estimation problem, and is applicable for arbitrary kernel functions. We compare our designed algorithm with the well-known implementations from the scikit-learn library and the FAISS library, and find that our method significantly outperforms the implementation from both libraries on a variety of datasets.
1 Introduction
Given a set X (.1 } C Rd of data points and similarity function k: Rd x Rd - Rz0 for any pair of data points x; and zy, the objective of clustering is to partition these n data points into clusters such that similar points are in the same cluster. As fundamental data analysis technique, clustering has been extensively studied in different disciplines ranging from algorithms and machine learning, social network analysis, data science and statisties.
One prominent approach for clustering data points in Euclidean space consists of two simple steps: the first step is to construct similarity graph K (V, E, w) from X, where every vertex v; of corresponds to x; E X, and different vertices v; and vy are connected by an edge with weight w(r, 0j) if their similarity k(x1, is positive, higher than some threshold. Secondly, spectral clustering on G, and its output naturally corresponds to some clustering on [19]. Because of its out-performanc over traditional clustering algorithms like k-means, this approach has become one of the most popular modern clustering algorithms.
apply
On the other side, different constructions of similarity graphs have significant impact on the quality and time complexity of spectral clustering, which is clearly acknowledged and appropriately discussed by von Luxburg [31. Generally speaking, there are two types of similarity graphs:
the first one the h--narestt neighbour graph (k-NN graph), in which every vertex U; connects to ty if U is among the k-nearest neighbours of u; A k-NN graph is sparse by construction, but loses some of the structural information in the dataset since k is usually small and the added edges unweighted.
the second one the fully connected graph, in which different vertices i and #y are connected with weight w('; 02) k(x;, While a fully connected graph maintains most of the distance information about X, this graph is dense and storing such graphs requires quadratic memory in
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Taking the pros and cons of the two constructions into account, one would naturally ask the question:
Is it possible directly construct sparse graph that preserves the cluster structure of fully connected similarity graph?
We answer this question affirmatively, and present fast algorithm that constructs an approximation of the fully connected similarity graph. Our constructed graph consists of only O(z) edges and preserves the cluster structure of the fully connected similarity graph.
1.1 Our Result
Given any set X {01..... .......,...... In. C. I' and kernel function k R x R + R2o. fully connected similarity graph K = (V, F,, w) of X consists of n vertices, and every o2 V corresponds to E X; we set w(v;, Uz) 2 k(r; xyy for any different and vj. We introduce efficient algorithm that constructs a sparse graph directly from X, such that K and G share the same cluster-structure, and the graph matrices for K and G have approximately the same eigen-gap. This ensures that spectral clustering from G and return approximately the same result.
The design of our algorithm is based on novel reduction from the approximate construction of similarity graphs to the problem Kernel Density Estimation (KDE). This reduction shows that any algorithm for the KDE can be employed to construct a sparse representation of fully connected similarity graph, while preserving the cluster-structure of the input data points. This is summarised as follows:
Theorem 1 (Informal Statement of Theorem 2}.
Given a set of data points X = {x1,-.... ......Tnn C TRd input, there is randomised algorithm that constructs sparse graph G of X, such that it holds with probability at least 9/10 that
1. graph G has Ön) edges,
2. graph G has the same cluster structure as the fully connected similarity graph K of X.
The algorithm uses an approximate KDE algorithmm as a black-box, and has running time ((TKEE (n,n,,)) e)) for e S 1/(6 log(n)), where [KDE(n, is the running time of solving the KDE problem for data points up to a (1 te)pppproximaiinnn
This result builds a novel connection between the KDE and the fast construction similarity graphs, and further represents state-of-the-ar algorithm for constructing similarity graphs. For instance, when employing the fast Gauss transform [TO] as the KDE solver, Theorem shows that a sparse representation of the fully connected similarity graph with the Gaussian kernel can be constructed in ((n) time when d is constant. As such, in the case of low dimensions, spectral clustering runs as fast (up to poly-logarithmic factor) as the time needed to read the input data points. Moreover, any improved algorithm for the KDE would result in faster construction of approximate similarity graphs.
To demonstrate the significance of this work in practice, compare the performance of our algorithm with five competing algorithms from the well-known scikit-learn library (20] and FAISS library [II]: the algorithm that constructs fully connected Gaussian kernel graph, and four algorithms that construct different variants of k-nearest neighbour graphs. We apply spectral clustering the six constructed similarity graphs, and compare the quality of the resulting clustering. For typical input dataset of 1 15,000 points in R.? 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 I all the six algorithms return reasonable output.
We further compare the quality of the six algorithms 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-lean able to segment the full-resolution image with over 150,000 pixels. Our experimental result on
typical laptop. In contrast, our algorithm
We use Ö(n) to represent O log((n) for some constant

---

SKLEARN N-
SKLEARN GK
FAISS EXACT
FAISS HNSW
FAISS IVF
OUR ALGORITHM
95.91s
65.644
48.40s
26.76s
4.46s
128 88g
56.48s
73.34s
37.87 87
47. 10s
4.52s
133
55.085
31,20
12
16,075
Figure Output of spectral clustering with different similarity graph constructions.
the BSDS dataset is showcased Figure 2 and demonstrates that, in comparison with SKLEARN GK, our algorithm identifies a more detailed pattern on the butterfly's 's wing. In contrast, none of the k-nearest neighbour based algorithms from the libraries is able identify the wings of the butterfly.
(a) Original Image
(b) SKLEARN
(c) OUR ALGORITHM
(d) SKLEARN i--N
(e) FAISS EXACT
() FAISS HNSW
(g) FAISS IVF
Figure 2: Comparison on the performance of spectral clustering with different similarity graph constructions. Here, SKLEARN GK is based on the fully connected similarity graph construction, and (d) (g) based different k-nearest neighbour graph constructions from the two libraries.
1.2 Related Work
There are a number of works on efficient constructions of e-neighbourhood graphs and k-NN graphs. For instance, Dong et al. [9] presents an algorithm for approximate k-NN graph construction, and their algorithm based on local search. Wang et [32] presents an LSH-based algorithm for constructing approximate k-NN graph, and employs several sampling and hashing techniques to reduce the computational and parallelisation cost. These two algorithms [ 32) have shown to work very well in practice, but lack 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 class of unbiased estimators for kernel density in high dimensions for 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. Quan- rud [22] studies smooth kernel functions, and shows that an explicit (1 + aaaapooximatee spectral approximation of the geometric graph with (nn/²²) edges can be computed in (n//s²) time. Bakshi al. [5] proves that, under the strong exponential time hypothesis, constructing
(a1))pppoxiimat spectral sparsifier with O edges for the Gaussian kernel graph requires Q n 2log(1/>)0
time, where < 0.01 is a fixed universal constant and r is the minimum entry of the kernel matrix. Compared with their results, we show that, when the similarity graph with the Gaussian kernll presents a well-defined structure of clusters, an approximate construction this similarity graph can be constructed in nearly-lineas time.
2 Preliminaries
Let G = (V, E, wa) be an undirected graph with weight function WG E - R20, and n & [V]. The degree of any vertex v is defined as degc(v) _ Enoo wG(u,, v), where we write a ~ {u, v) E For any S C V. the volume of S is defined by volc(S) Eves dvg, and the conductance of S is defined by
(c(S) 0G(S)
vola(S)
.,, Aş. a k-way partition if A; + 0 for any 1 sis k. A, 0 A₂ 0 for any i + j. and U A; V. Moreover, define the k-way expansion constant by
where 0c( & Eves,ogs wc(u, v). For any 2 2, we call subsets of vertices Ai....,
PC Pc(k) A
max de (A,). purtition A Ag Isisk
min
Note that a lower value of pc(k) ensures the existence of & clusters A A As of low conductance, ie, G has at least & clusters.
For any undirected graph G, the adjacency matrix AG of G is defined by Ac(u, v) wG(u, v) v, and Ac(u, v) = 0 otherwise. We write De as the diagonal matrix defined by Da(u, dvga(v), and the normalised Laplacian of G is defined by Nç e D-D 1/2 A 1/2 For any PSD
matrix B Raxn, write the eigenvalues of B as X (B) S ... < (B).
well-knows that, while computing (k(k) exactly is NP-hard, pc(k) is closely related to A through the higher-order Cheeger inequality [3: it holds for any k that
Au (Nc) < (k(k) S O(k)) ((k)))><<Ne))))
2.1 Fully Connected Similarity Graphs
We use A [x1,....i to represent the set of input data points, where every x; Rd, Given X and some kernel function I E R20, we use K = (Vk. Ek, wk) to represent the fully connected similarity graph from X, which is constructed as follows: every E Vk corresponds to a; E X, and any pair of different v; and ey is connected by an edge with weight ux uk(u;, 02) = (.)) One of the most common kernels used for clustering is the Gaussian kernel, which is defined by
||:;
k(2,, 2y) exp
for some bandwidth parameter a. Other popular kernels include the Laplacian kernel and the exponential kernel which use || æ1 and ||, 225ll2 in the exponent respectively.

---

2.2 Kernel Density Estimation
Our work is based on algorithms for kernel density estimation (KDE), which defined as follows. Given a kernel function k R x I² + R20 with source points .i..... 2n E R and m target
points ....... Vm E RR, the KDE problem is to compute 9[1.x) (V1), 91,,) (Um), where
(1.) (16)
(1)
for 1 iiiS m. While direct computation of the m values 911, operations, there is substantial research to develop faster algorithms approximating these m quantities.
Ii.,ni (y1)... 911.. (ym) requires mn
In this paper we are interested in the algorithms that approximately compute (1s) for all 1 S < up to (1 + mulliipliicatie error, and use (m, e) to denote the asymptotic complexity of such a KDE algorithm. We also require that (mDE(m, n, e) is superadditive in m and n; that is, for my m2 and n2 n2,, have
Tm n1. <) + TKDE(m2, r22,e) e) S TKDEE(M,,,)) n1,))
it is known that such property holds for many KDE algorithms (e.g.. [ 16 10
3 Cluster-Preserving Sparsifiers
A graph sparsifier a sparse representation of an input graph that inherits certain properties of the original dense graph. The efficient construction of sparsifiers plays an important role in designing number of nearly-linear time graph algorithms. However, most algorithms for constructing sparsifiers rely the recursive decomposition of input graph [26). sampling with respect effective resistances [15] 25), or fast SDP solvers [[4]; all of these need the explicit representation of an input
graph, requiring ( time and space complexity for fully connected graph.
Sun and Zanetti 27 study a variant of graph sparsifiers that mainly preserve the cluster structure of an input graph, and introduce the notion of cluster-preservigg sparsifier defined as follows:
Definition 1 (Cluster- peeettrrpreervingg sparsifier). Let = EV wk) be any graph, and {A,) the k-way partition of K corresponding Pk(K). We call re-weighted subgraph G (V, F C E, wG) a cluster-preserving sparsifier of K if (A,) = O(k Φx (A:)) for 1 < k; and (NG) =
((Ax+1 (Nk)).
Notice that graph K has exactly & clusters (i) K has i disjoint subsets A AI...., Ag of low conduc- tance, and (ii) any (A: 1)-way partition of K would include some A C V of high conductance, which would be implied by a lower bound on Ax+1 (Nk) due to the higher-order Cheeger inequality. Together with the well-known eigen-gap heuristic [13 31 and theoretical analysis on spectral cluster-
ing [16] 21), the two properties in Definition
1 ensures that spectral clustering returns approximately
the same output from K and H
Now we present the algorithm in [27 for constructing a clluster-preserving sparsifier, and we call the SZ algorithm for simplicity. Given any input graph K (V, E, wk) with weight function wk. the algorithm computes
(v) A min {c. log wK(u, v)
1 J
and (u( e min fc log WK(U, u) 1) 1
Ae+1 degx(v) :
AR+1 degx(u)
for every edge {u, v}, where C' R+ is some constant. Then, the algorithm samples every edge e = {u, o. with probability
A Pu (v) + Pu()) pa(v) (u), Pe
and sets the weight of every sampled e = {u, v) in G as wG v) & wK(u, u)/Pe. By setting F as the set of the sampled edges, the algorithm returns G (V, wG) output. is shown in 27 that, with high probability, the constructed has ((m) edges and is dluster-preservin sparsifier of K.
'The most interesting regime for this definition is A Ō(1) and ( = 2(1), and we assume this in the rest of the paper.
5

---

We remark that a cluster- preseeving sparsifier is much weaker notion than a spectral sparsifier, which approximately preserves all the cut values and the eigenvalues of the graph Laplacian matrices. On the other side, while cluster-preserving sparsifier is sufficient for the task of graph clustering, the SZ algorithm runs (nn") time for fully connected input graph, since it's based on the computation of the vertex degrees as well as the sampling probabilities Pu (v) for every pair of vertices and v.
4 Algorithm
This section presents our algorithm that directly constructs approximation of fully connected similarity graph from X S I with |X| = n. As the main theoretical contribution, we demonstrate that neither the quadratic space complexity for directly constructing fully connected similarity graph nor the quadratic time complexity the SZ algorithm is necessary when approximately constructing fully connected similarity graph for the purpose of clustering. The performance of our algorithm is as follows:
Theorem 2 (Main Result). Given set of data points X = {xI1...., .In) C RRd as inpul, there is randomised algorithm that constructs a sparse graph G of X, such that it holds with probability at least 9/10 10 that
1. graph G has (n) edges,
2. graph G has the same cluster structure as the fully connected similarity gruph K of X; that is, if K has k well-defined clusters, then it holds that pc = O pk(()) and (NG) = ((2++1((K)))
The algorithm uses an approximate KDE algorithm as a black-box, and has running time Ô(TrDE n, e)) for S 11((log((m)
4.1 Algorithm Description
At a very high level, our designed algorithm applies a KDE algorithm as black-box, and constructs dluster-peserrving 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 2; a random æj can be sampled with probability k(2,, / (eg)((;) through Oloo n) queries to KDE algorithm. To see this, notice that we can apply a KDE algorithm to compute the probability that the sampled neighbour some set X1 X, i.e.
Σ k(z,)z)) (x, 2yEX, degx(1;)
gx (x:)
where Ox (3:) 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 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 random neighbour of a vertex in through KDE and binary search also appears in Backurs al. [
1. Set the feasible neighbours to be X {ø1.................... ).
2. While X| >> 1:
Split X into X1 and X2 with |X1 [|| and |x21= [1x]/2].
Compute gx (x;) and gx, (z+); set X < X1 with probability 9x, (x:) (x;). and 9x, (æ.)/gx (æ:).
X + X₂ with probability
3. Return the remaining element in X the sampled neighbour.
Next generalise this idea and show that, instead of sampling neighbour of one vertex at time, KDE algorithm allows us to sample neighbours of every vertex in the graph "simultaneously". Our designed sampling procedure is formalised in Algorithm T
Finally, to construct cluster-preservine sparsifier, we apply Algorithm I to sample Ollog n) neigh- bours for every vertex V; and set the weight of every sampled edge et Uj
k(0,, xT))
(U, wG WG(U,, 0z)) p(i,j)
(2)

---

Unif[0, 1]
Step 1
(1.2)/2(225)
9|n/2,n (2+))
Step 2
(0,n/4/((z))
/n/44,//2(22
Step (ng
k( mry)
Figure 3: The procedure of sampling a neighbour vj of u; with probability k(a; xr (v:). Our algorithm performs binary search to find the sampled neighbour. At each step, the value of two kernel density estimates are used to determine where the sample lies. Notice that the algorithm doesn't compute any edge weights directly until the last step.
Algorithm 1 SAMPLE
1: Input: set S of {y:}
10:
Si S2 = 0
set X of {z;)
11:
for yi E S do
2: Output:
12:
T~ Unif
E {1i,, m.) for some and j}
13:
if r s gx, (,x)/(xx,(x) + gxx(v)))
3: if |X| = then
then
4:
return S x X
14:
S S U {ri} else
5: else
15:
6:
X1
[z j
X /2]
16:
S S₂ U {3
7: X₂ - : X /
17:
end if
8:
Compute gx, (3:)
for all with KDE al-
18:
end for
gorithm
19:
return SAMPLE(S), X, XI)SSAMPLE(S, X2)
9: Compute 9x, (3:) for all with KDE al-
20:
end if
gorithm
where ((i.j & P1()+ +7, P.0) (1) is an estimate of the sampling probability of edge Uj. and
k(ae, xrn)
(()) A min 6C logn
Mp.m (x:)
for some constant C E R+; see Algorithm 2 for the formal description.
4.2 Algorithm Analysis
Now we analyse Algorithm 2 and sketch the proof of Theorem We first analyse the running time of Algorithm 2 Since it involves O(log n) executions of Algorithm total, sufficient to examine the running time of Algorithm [
(So,
X1,11
(S1,2, X1,2)
We visualise the recursion of Algorithm [ with spect to and in Figure 4 Notice that, although the number of nodes doubles each level of the recursion tree, the total number of samples S and data points X remain constant for each level of the
(2,11,
Figure 4: The recursion tree for Algorithm
tree: it holds for any that U St. So and

---

Algorithm 2 FASTSIMILLRTTYGRAPH
1: Input: data point set X {II (I1....,Zn)
8: for (U, 0y) E do
2: Output: similarity graph
9:
(.(1) min { ////
3: E == 0, L 6C (og(n))
10:
P5( min {L k( 1/
4: for E [1, do
11:
p(i,j) = p:( Pr p
5: E EEUSAMPLEXXA
12: Set wa(w;, k(2,,,;;//(L..
6: end for
13: end for
7: Compute 911,1 (xi) for each i with KDE
14:
return graph (X, E, wc)
algorithm
- Xo. Since the running time of the KDE algorithm is superadditive, the combined
running time of all nodes at level of the tree
2'
2'
2'
T = Troe
& TKDE
Σ
Sh Xislic (n, n, n)...
Σ j=1
Hence, the total running time Algorithm is ÖTKoE (n, n, as the tree has log2(n) levels. Moreover, when applying the Fast Gauss Transform (FGT) [10] as the KDE algorithm, the total running time of Algorithm 1 (n) when d O())
Finally, we prove that the output of Algorithm 2 is a cluster preserving sparsifier fully connected
similarity graph, and has (()) edges. Notice that, comparing with the sampling probabilities Pu. (v) and p. used in the SZ algorithm, Algorithm 2 samples each edge through Olog n) recursive executions of a KDE algorithm, each of which introduces (1 + --multiplicative error. We carefully examine these (1 + m)muliiplicaaive errors and prove that the actual sampling probability p(i, j) used in Al Algorithm is an
of p. for every e = (vi, )- As such the output of Algorithm is a cluster-preserving sparsifier fully connected similarity graph, and satisfies the two properties of Theorem 2 The total number of edges in G follows from the sampling scheme. We
refer the reader to Appendix A for the complete proof of Theorem
Experiments
In this section, empirically evaluate the performance of spectral clustering with our new algorithm for constructing similarity graphs. We compare our algorithm with the algorithms for similarity graph construction provided by the scikit-lee -earn library (20 and the FAISS library [IT] which are commonly used machine learning libraries for Python. In the remainder of this section, we compare the following six spectral clustering methods.
1. SKLEARN GK: spectral clustering with the fully connected Gaussian kernel similarity graph constructed with the scikit-learn library.
2. SKLEARN k-NN: spectral clustering with the k-nearest neighbour similarity graph con- structed with the scikit-learn library.
3. FAISS EXACT: spectral clustering with the exact k-nearest neighbour graph constructed with the FAISS library.
4. FAISS HNSW: spectral clustering with the approximate k-nearest neighbour graph con- structed with the "Hierarchical Navigable Small World" method [1]
5. FAISS IVF: spectral clustering with the approximate k-nearest neighbour graph constructed with the "Invertex File Index" method B
6. OUR ALGORITHM: spectral clustering with the Gaussian kernel similarity graph constructed by Algorithm 2
We implement Algorithm Yang et al. (33, and use the Python SciPy (30 and stag [17] libraries for eigenvector computation and graph operations respectively. The scikit- learn and FAISS libraries well-optimised
2 in C++, using the implementation of the Fast Gauss Transform provided by

---

and use C, C++, and FORTRAN for efficient implementation of core algorithms. We first employ classical synthetic clustering datasets to clearly compare how the running time of different algorithms is affected by the number of data points. This experiment highlights the nearly-linear time complexity of our algorithm. Next we demonstrate the applicability of our new algorithm for image segmentation on the Berkeley Image Segmentation Dataset (BSDS) (2).
For each experiment, we k 10 for the approximate nearest neighbour algorithms. For the synthetic datasets, we set the value of the Gaussian kernel used by SKLEARN GK and OuR ALGORITHM to be 0.1, and for the BSDS experiment we set σ 0.2. This choice follows the heuristic suggested by von Luxburg [31 to choose o to approximate the distance from point to its ktth nearest neighbour. All experiments are performed on an HP ZBook laptop with an 11th Gen Intel(R) Core TM) 17-11800H @ 2.30GHz processor and 32 GB RAM. The code to reproduce our results is available at battps://githhe cwm. hoom/pmacg/kde-eee http:/ddeesimmlmaat..... graph
5.1 Synthetic Dataset
In this experiment we use the scikit-learn library to generate 15,000 data points in IR² from variety of classical synthetic clustering datasets. The data is generated with the make make_circless circles, make_moons, and make_blobs methods of the scikit-leara 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-leary and FAISS.
We further compare the speedup of our algorithm against the five competitors on the two two__ow _oooon dataset with increasing number of data points, and result is reported in Figure Notice that the running time of the scikit-leary and FAISS algorithms scales roughly quadratically with the size of the input, while the running time of our algorithm 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.
80
80
SKLFARN NNNN
70
70
SKLEARN GK
60
60
FAISS EXACT
FAISS HNSW
Rumning time
50
Rouning time
50
FAISS IVF
40
OUR ALGORITHM
30
30
20
10
10
5,000
15,000
20,000
50,000
100,00
Number of data points
Number data points
Figure 5: Comparison of the running time of spectral clustering on the two moons dataset. In every case, all algorithms return the correct clusters.
5.2 BSDS Image Segmentation Dataset
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)T R where r, g. are the red, green, blue pixel values and x, 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 reported in Table [ 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).

---

Table 1: The average Rand Index of the output from the six algorithms. Due to the quadratic space complexity constraint, the SKLEARN GK cannot be applied the full resolution image.
Algorithm
Downsampled
Full Resolution
SKLEARN GK
0.681
SKLEARN k-NN
0.675
0.682
FAISS EXACT
0.675
0.680
FAISS HNSW
0,679
0,677
FAISS IVF
0,675
0.680
OUR ALGORITHM
0,680
0.702
6 Conclusion
In this paper we develop a new technique that constructs a similarity graph, and show that an approximation algorithm for the KDE can be employed to construct similarity graph with proven approximation guarantee. Applying the publicly available implementations
the KDE as a black- box, our algorithm outperforms five competing ones from the well-known scikit-leane and FAISS libraries for the classical datasets of low dimensions. We believe that our work will motivate more research faster algorithms for the KDE in higher dimensions and their efficient implementation, well as more efficient constructions of similarity graphs.
Acknowledgements
We would like to thank the anonymous reviewers for their valuable comments on the paper. This work is supported by EPSRC Early Career Fellowship (EP/T00729X/I).
References
[1] Josh Alman, Timothy Chu, Aaron Schild, and Zhao Song. Algorithms and hardness for linear algebra on geometric graphs. In 61st Annual IEEE Symposium on Foundations of Computer Science (FOCS pages 541-552 2020.
[2] Pablo Arbelaez, Michael Maire, Charless Fowlkes, and Jitendra Malik. Contour detection and hierarchical image segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 33(5):898-99) 01088-996127 2011.
[3] Artem Babenko and Victor Lempitsky. The inverted multi-index. IEEE Transactions on Pattern Analysis and Machine Intelligence, 37(6):1247-1266 2014.
[4] Arturs Backurs, Piotr Indyk, and Tal Wagner. Space and time efficient kernel density estimation in high dimensions. In 33rd Advances in Neural Information Processing Systems (NeurIPS' 19), pages 15773-15782, 2019.
[5] Ainesh Bakshi, Piotr Indyk, Pranceth Kacham, Sandeep Silwal, and Samson Zhou. Subquadratic algorithms for kernel matrices via kernel density estimation. In Ith International Conference on Learning Representations (ICLR '23), 2023.
[6] Moses Charikar, Michael Kapralov, Navid Nouri, and Paris Siminelakis. Kernel density estimation through density constrained near neighbor search. In 61 Annual IEEE Symposium on Foundations Computer Science (FOCS'20), pages 172-183, 2020.
[7] Moses Charikar and Paris Siminelakis. fashiing-based-sstimal for kernel density in high dimensions. In 58th Annual IEEE Symposium on Foundations of Computer Science (FOCS' 17), pages 1032-1043, 2017.
[8] Fan Chung and Linyuan Lu. Concentration inequalities and martingale inequalities: a survey. Internel mathematics, 3(1):79-127, 2006.
[9] Wei Dong, Moses Charikar, and Kai Li. Efficient k-nearest neighbor graph construction for generic similarity measures. 20th International Conference on World Wide Web (WWW 11)), pages 577-586, 2011.
10

---

[10] Leslie Greengard and John Strain. The fast gauss transform. SIAM Journal on Scientific & Statistical Computing, 12(1):79-99 -947 1991.
[11] Jeff Johnson, Matthijs Douze, and Hervé Jégou. Billion-scale similarity search with GPUs. IEEE Transactions on Big Data, 7(3):535-547, 2021.
[12] Zohar Karnin and Edo Liberty. Discrepancy. coresets, and sketches machine learning. In 32nd Conference Learning Theory (COLT 19), pages 1975-1993, 2019,
[13] James R. Lee, Shayan Oveis Gharan, and Luca Trevisan. Multiway spectral partitioning and higher-order Cheeger inequalities. Journal of the ACM, 61(6)37:1-37:300 2014.
[14] Yin Tat Lee and He Sun. An SDP-based algorithm for linear-sized spectral sparsification. In
49th Anmual ACM Symposium on Theory of Computing (STOC'17), pages 978-687 2017.
[15] Yin Tat Lee and He Sun. Constructing linear-sized spectral sparsification in almost-linear time.
SIAM Journal on Computing, 47(6)) --333,,,,,,,, 2018.
[16] Peter Macgregor and He Sun. A tighter analysis of spectral clustering, and beyond. In 39th International Conference on Machine Learning (ICML'22), pages 14717-14742, 2022.
[17] Peter Macgregor and He Sun. Spectral toolkit of algorithms for graphs: Technical report (1). CoRR, abs/2304.031700, 2023.
[18] Yu A Malkov and Dmitry A Yashunin. Efficient and robust approximate nearest neighbor search using hierarchical navigable small world graphs. IEEE transactions on Puttern Analysis and Machine Intelligence, 42(4):824-836, 2018.
[19] Andrew Y. Ng. Michael I. Jordan, and Yair Weiss. On spectral clustering: Analysis and an algorithm. In 14th Advances in Neural Information Processing Systems (NeurIPS'O1), pages 849-856, 2001.
[20] Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825-288 2830. 2011.
[21] Richard Peng. He Sun, and Luca Zanetti. Partitioning well-clustered graphs: Spectral clustering works! SIAM Journal on Computing, 46(2):710-743, 2017.
[22] Kent Quanrud. Spectral sparsification of metrics and kernels. 32nd Annual ACM-SIAM Symposium on Discrete Algorithms (SODA 22), pages 1445-1464, 2021.
[23] William M. Rand. Objective criteria for the evaluation of clustering methods. Journal of the American Statistical Association, 66(336)844-850, 1971.
[24] Paris Siminelakis. Kexin Rong. Peter Bailis, Moses Charikar, and Philip Levis. Rehashing kernel evaluation in high dimensions. In 36th International Conference on Machine Learning (ICML'I9), Proceedings of Machine Learning Research, pages 5789-5799. 2019.
[25] Daniel A. Spielman and Nikhil Srivastava. Graph sparsification by effective resistances. SIAM
Journal on Computing, 40(6):1913-199
2011.
[26] Daniel Spielman and Shang- Hua Teng. Spectral sparsification of graphs. SIAM Journal on Computing, 40(4):981-1025, 2011.
[27] He Sun and Luca Zanetti. Distributed graph clustering and sparsification. ACM Transactions on Parallel Computing, 63)::17:1 2019.
[28] Joel A Tropp. User-friendly tail bounds for sums of random matrices. Foundations of Computa- tional Mathematics, 12(4)3389-433 2012.
[29] Paxton Turner, Jingbo Liu, and Philippe Rigollet. Efficient interpolation of density estimators.
In 24th International Conference Artificial Intelligence and Statistics (AISTATS '21), pages 2503-2511, 2021.
[30] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, CJ Carey, ilhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E, A. Quintero. Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul
11

---

van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17:261-272, 2020.
[31] Ulrike von Luxburg. A tutorial spectral clustering. Statistics and Computing, 17(4):395-416 2007.
[32] Yiqiu Wang, Anshumali Shrivastava, Jonathan Wang, and Junghee Ryu. Randomized algorithms accelerated over CPU-GPU for ultra-high dimensional similarity search. In 2018 International Conference on Management of Data (SIGMOD '18), pages $89-903, 2018.
[33] Changjangg Yang. Ramani Duraiswami, Nail Gumerov, and Larry Davis. Improved fast Gauss transform and efficient kernel density estimation. In 9th International Conference on Computer Vision (ICCV'03), pages 674-671, 2003.
12

---

A Proof of Theorem 2
This section presents the complete proof of Theorem 2 equal to the indices of the L points sampled for æ;. Recall that by the SZ algorithm, the "ideal" sampling probability for z; from 2; is
Let Ii.l.... L be random variables which
(i( min cc. log(n)
k(rs; xz)
degg((;)) 1)
Au+1
We denote the actual sampling probability that zy is sampled from 2; under Algorithm to be
(i) A P [ E (vr, Co..I........ v...H
Finally, for each added edge, Algorithm 2 also computes an estimate of p:(pa) which we denote
(.()) min 6C log(n) k(ri,, 2tz)
Ak+1
[[1.n (24)
Similar with the definition p. in Section 3
we define
p(i,j) p.(O) (i) pz(j) (i))),
p(i.j) p.()) ()) (,0))-,(0 and
p(i.j) P.)) (1) - -A()) (())),(()
Notice that, following the convention of (27). we use (j) to refer to the probability that given edge is sampled from the vertex x; and j)i,j) is the probability that the given edge {w), Uy } is sampled at all by the algorithm. We use the same convention for ()( and (j).
We first prove sequence of lemmas showing that these probabilities are all within a constant factor of each other.
Lemma 1. For any point ou the probability that given sampled neighbour ys is equal to j is given by
k(I;,1;) 2 degg((s))
2k 2k(x,; 22))
< [3H, = 51 degk (v:)
Proof. Let X = [II....,, Znn } be the input data points for Algorithm 2. and [n] = (1..... n) be the indices of the input data points. Furthermore, let [a. b] {a....................... b} be the set of indices between and b. Then, in each recursive call to Algorithm| 1 we are given range of indices [u, b] input and assign Yı, to one half of it: either (a, [b/2]] or [[b/a] 1. By Algorithm| we have that the probability of assigning Vi, to (a, [b/2]] is
Dla. [b/2] (I;)
P Iy, E [a, [6/2]] Vaa E [a, b]
ga (I:)
By the performance guarantee of the KDE algorithm, have that (x)) E + c) degja, (v:). where we define
deg (I)))) - Σ k(),,
J=a
This gives
1 degia, 16/21 (v:) (+t)
1 1+0
10/21 (3)
deg (v)) < P [Bn. E X(a.1//2 Iu. E X(ao)] $
1
degia.b) (v:)
Next, notice that we can write for a sequence of intervals (a1, bı] C [az, be] c... c [1, n] that
P [yu. = f] [yr, = jlyi, E (ai, bi]] x IP [ya. E bilbbrl E [a2, b2]]
X x P [VA E [ax, by blb.. E [1,n]],
13

---

where each term corresponds to one level of recursion of Algorithm lnnd there are at most [log2(1) terms. Then, by (3) and the fact that the denominator and numerator of adjacent terms cancel out, have
1-6 (1))
(5) degk (2)) _ I [Vi. = =
(1 + 1-
[1gg((n) k())
5
1+e)
degk (u,)
since degrs,) (vi) ) = k(22, 2y)) and degsu,)) = degk (v:).
For the lower bound, we have that
[x)2(1)]
1 I+c
2 - (e) [082(1)) -1-3 -1-3logz(n)e 2 1/2,
where the final inequality follows by the condition e that e S 1 For the upper bound, similarly have
[log2(n]] (1 + 3c) [[og[(]]] & exp (3(log2(n)]c) & 02/3 2
1+
1-c
where the first inequality follows since e 1/(6l0g2
The next lemma shows that the probability of sampling each edge {or Uy} is approximately equal to the "ideal" sampling probability Pı(j).
Lemma 2. For every and j * i, we have
< 12 px()).
1O
Proof. Let Y = 1......... 1 be the neighbours of x; sampled by Algorithm 2 where = 6Clog(n)/Aes Then,
PÍEYY = II (1 1 11213 (1-
k( x4) )² exp (-L.
A(x;, xr2). degk (v:)
-1)=-1)1---- 2deg (w:)
The proof proceeds by case distinction.
Case 1: p((j) 9/10.. In this case, we have
1 2 p(()).
Case 2: p(()) > 9/10. In this case, we have
9
P P
20
10
which completes the proof on the lower bound of (.))).
For the upper bound, we have
PhjMss1 (1- 2k Rz5)
2k(r;, degx
log(n) k(x;, 2.,) degk(u.) (v;)'
12C
degx (vs)
Aut1
from which the statement follows.
An immediate corollary of Lemma is as follows. Corollary 1. For all different i,j [n], it holds that
9
10 p(.,)) S p(i,)) S 12 p(i.j)
and
7 p(i.)) S p(1.)) < 36
6
p(i,)).
14

---

Proof. For the upper bound of the first statement, we can assumee that p:(5) S 1/12 and p,(i) S 1/12. since otherwise we have p(i. j) - I 12-p/ 12-p(i, j) and the statement holds trivially. Then, by Lemma 2 we have
p(i.)) ()) (1)
$ 112p())+1pp,0))220 -2) 12 (p:(())++ -x()))))())))))) .p)()) p(i.j)
and
= - + 1,(0
9
(x: + 9
10 (1)
T0 p.
TO
TO 20)
9
10
9
(p:0)+p;00 - p((ips(())
10
9
p((,)).
10
For the second statement, notice that
6Clog(n) k(i.))
(i( = min
Aut1
C'lo((n) k(i,j) 1}
min
1+8 Aut1
(1))
1
C'log(n) k(i,5)
2
min
1+e
Act1
degk(1s))
1
-P:(i)
1+8
P.()).
since (x;) is (1 + s) approximation of degk (v:) and < 1/6. Similarly,
6 (.(i) - 1-E pi(j) p.())-
36
Then, for the upper bound of the second statement, we can assume that p;( S 5/36 and (i < 5/36, since otherwise ()i.j) S < (36/5) ((i.j) and the statement holds trivially. This implies that
p(i,j) = P:(5) + (1)) 7,(1
36
36
36
< (i)
36
< (v:(i) - (1)
36
5 p(i.j)
and
p(1,j) > -p,(1) (i) -p.(i) 6 (p,(0)
6
7 (p:() pz()) p(()) (z)
7 p(i,)),
which completes the proof.
15

---

We are now ready to prove Theorem 2 It is important to note that, although some of the proofs below are parallel that of (27]. our analysis needs carefully take into account the approximation ratios introduced by the approximate KDE algorithm, which makes our analysis more involved. The following concentration inequalities will be used in our proof.
Lemma (Bernstein's s Inequality () Let XI,..., ,x, be independent random variables such that Then, it holds that
[X, <
M for any i E (I.....n). Let X Ei= Xs, and R >i=1 E [X?].
P| - E 2 2cx (
t²
2(R + Mt/3)
Lemma 4 (Matrix Chernoff Bound (28) Consider finite sequence {Xi) of independent, random,
PSD matrices of dimension that satisfy || S R. Let Pmnin Amia (B 12, X:)) and Pmax Amax (E 1 X,)). Then, it holds that
Σxx, Xi
6 S(1 0) <d (1
Mmin
P Amin
8) 1-8
for 8 E [0, 1), and
e
Minxx/
P
Amax
Σx,
(1 + 8) Mmax
S
(1+8)1+ 8)143
forr 8 2 0.
Proof of Theorem 2 We first show that the degrees of all the vertices in the similarity graph K are preserved with high probability the sparsifier G. For any vertex U; let yi. Vi,L be the indices of the neighbours of sampled by Algorithm 2
For every pair of indices * and for every < l S L, define the random variable Zeja to be the weight of the sampled edge j is the neighbour sampled from at iteration 1, and O otherwise:
k(aan))) pOJ)
yu-
Zi
otherwise.
Then, fixing an arbitrary vertex æi, we can write
dego(w:)
Zush + Zza.1-
We have
[dego (m:) =
[71g,1] + E [[j.,..
1=1 8z
k(1r) p(i.5)
k( 24
TP li. =
P( = 1).
f(i,g)
By Lemmas and and Corollary
we have
k(1;;
L k.)(),))))))
L kk(a,, 2.))
[dega(v:)] 2
2 2e (v:)
2 degk (v²)
3k xz))
(p.(i) pp(())
p(i,j)
5k(x; ..)
5 degg (1:)
12
12
16

---

where the last inequality follows by the fact that p(i,5) S (36/5)-p(e) - (36/5) (p:(i) + py)) Similarly, we have
2 (2--A((Z,,
2 L (b((2))))))) degk(ug)
E [dego(
5
p(i,5)
degk(wi)
jti
12. 2g)
(p:(i) + +z)(1)
p(i,j)
jti
Σ 28 k(#;, = 28 (v:),
where the inequality follows by the fact that
6 p(i,)) 7 p(i.j) - (p:(i + p,) (;)) (p:(i (,))
In order to prove a concentration bound on this degree estimate, we would like to apply the Bernstein inequality for which need to bound
R
( +E [2
iti
[14, =] (xx,221)²
+ Iypu = i] k(z1,1,) 5(i,5)²
p(i.j)²
12 2(,;,
(p:() + p.))
p(1,))2
321
2,)
Σ
28
p(i,j)
Σ
98
k(,,,,,)²
3
(i(
98
k(x,, x,))
degk (vi) Au+1
3
C Clog(n)
98
degk (v:)² AR+1
C Clog(n)
where the third inequality follows by the fact that
6
p(i.)) 2 p(i.)) 2 7 p.()).
Then, by applying Bernstein's inequality we have for any constant C2 that
degk (v:k(c)))/(2) C3)
P (dege(rs) (v:) [(degg((:)]]
degk(1;))
< 2exp
C2
98 (v. A
7 degg (hett CCC2--lo(m)
30 log(n)
Cogg()))
$ 2 2exp
((196/3) C2 + (7/3) C2) Abt1
=0(1nn)
where we use the fact that
7 ddeg (v.) Ax+1
Zisa <
6p:(7)
6C log(n)
Therefore, by taking to be sufficiently large and by the union bound, holds with high probability that the degree of all the nodes in G are preserved up constant factor. For the remainder of the
17

---

proof, we assume that this is the case. Note in particular that this implies volG(5) = ( (8)) for any subset S C V.
Next, we prove holds for G that (G((;)) = O (k .o((s)))))))))))))))))))))))). for any S i S k, where , Sk form optimal clustering in K.
By the definition of Zr holds for any 1 S S k that
L ΣΣΣ
E [S.))
E
Σ
Σ
Za,bat + Zb
aES, bES; 1=1
( Tbb)
(pa +p(b))
Ma, b)
aeS; bgS.
O (Ok(S.))
where the last line follows by Corollary By Markov's inequality and the union bound, with constant probability it holds for all i = l...,, k that
@c(S.) O(k (k()))
Therefore, it holds with constant probability that
(k(k) & max (G())) = max O xx = O(k O(k--((k)))
Isisk
1SiSk
Finally, we prove that A (Nc) eigenspaces, and notice that Nk can be written
2(A++1 (Nk)). Let Nu be the projection of Nk on its top n
NK Σ (NNkffff
where f1,...,Ja fn are the eigenvectors of Nk. Let N 1/2 be the square root of the pseudoinverse of Nk.
We prove that the top eigenvalues of NK preserved, which implies that ( ((4++1 (Nk)). To prove this, for each data point x; and sample 1 S / < L, we define random matrix X E Rax= by
X = wG(v,, Uz)) NK 1/2 b.
1/2
where Vi,ls be A Xv: Xv, the edge indicator vector, and Dap R defined by
Xx, (a) 1 { donu(os) (v:)
if otherwise.
Notice that
Xu
Σ
wa( 0z) Nx 1/2 b.b.N
1/2
NK 1/2
N'
1=1
sampled edges (1:03)
where
N'C =
Σ
wo(u,, vy) bab
sampled edges e= {
the Laplacian matrix of G normalised with respect to the degrees of the nodes in K. We prove that, with high probability, the top eigenvectors of N' and NK are approximately the same. Then, show the same for Nç and N' which implies that A (Nc) ((k))1(Kk)..
18

---

We begin by looking at the first moment of Ei=1 El=r X and have that
Amin
E
Xay 1=1
L
1
Amin
L
P Ply. = 5 (()1)
1/2 b/ b.bINK Nx2
NK
1/2
p(i.5)
fFi (11155)
k(I,
NK b.b.NK 1/2
1/2
Amin
3p:(i)
p(i,j)
iti
Amin
5
NK
NK
122)
1/2
12
12'
where the last inequality follows by the fact that
36
36
((i,)) S
p(i,j)
<
(vx(i)
p)))).
5
5
Similarly,
Amax
E
Σ
Xia
Amax
Pin.a 5] k(a)))))
NK 1/2
12 b.b.NK
1/2
((1,5)
6= (1500))
Nz1/2 1/2 1/22
1/2
Amax
12 p((j)
NN b.bjNK
p(i,5)
(pcion)
S Amax
28 NK
1/2
N
Nk ) 28,
1/2)
where the last inequality follows by the fact that
3 p(i.i) 2 7 p(i, 7 (p:(f) (i))
Additionally, for any and j V4.1 and {us, Uy we have that
b 1/2 Nk
1/2 be
|X s wo(u,, 1zy)
k(
.b,K' be
p(i,j)
bell
p(i,j)
Ae+1
7 Au+1
1
1
+
3C 3C log(n)
Ac+1
degx(u:) (v:)
degk (1,)
Jeggtu.]
degn (v,)
7
3o')))())
Now, we apply the matrix Chernoff bound and have that
1/2
12C-log(n)
IP
Amax
ΣΣ X
42
S m
OL/nn)
( (1+ 1/2)3/2
1=1
for some constant c. The other side of the matrix Chernoff bound gives us that
ΣΣ X
PP
Amin
< 5/24
< 0(1/n9).
i=1
1m1
19

---

Combining these, with probability 1 - O(1/n²) it holds for any non-zero E R" in the space spanned by fk+1...... fn that
1/2 N 1/2 E (5/24,42)
xTx
By setting y NK 1/2
1/2/2
rewrite this
N
GV
y NNN
y NN yJy y yTNky
E (5/24,42).
VTN
K
y
y yTNKY
Since dim (f (n)) = n - k: we have shown that there exist n k orthogonal vectors whose Rayleigh quotient with respect NG is ((4++1(Nk). By the Courant-Fischar Theorem, have (NN) (Nk))
only remains to show that A6+1 (NG) a(Ax+ (Nk)). By the definition of I we have that NG DG 1/2 D N D2 Therefore, for any x E R and =D De x, holds that
(((A++((N()) which implies that AІ-1 (Ng) 1/2
x N
y N IT
y yTy
= Ω
where the final guarantee follows from the fact that the degrees G are preserved up to a constant factor. The conclusion of the theorem follows by the Courant-Fischer Theorem.
Finally, we bound running time of Algorithm 2 which dominated by the recursive calls to Algorithm[] [ We note that, although the number of nodes doubles at each level of the recursion tree (visualised in Figure the total number of samples S and data points X remain constant for each level of the tree. Then, since the running time of the KDE algorithm is superadditive, the total running time the KDE algorithms at level i the tree is
T:
TKDE Aaoee((S, (|Stsb(X.slie)
2'
2
< TKDE
Stal. Σ Naslle = TRDE (RDE((S), (X),e..
Since there are Ollog(n)) levels of the tree, the total running time of Algorithm O(TKOE([S); ()). This completes the proof.
is
B Additional Experimental Results
this section, we include Figures 6 and some additional examples of the performance of the six spectral clustering algorithms the BSDS imgge segmentation dataset. Due to its quadratic memory requirement, the SKLEARN GK algorithm cannot be used on the full-resolution image. Therefore, present its results on each image downsampled to 20,000 pixels. For every other algorithm, show the results on the full-resolution image. In every case, we find that our algorithm able to identify refined detail of the image when compared with the alternative algorithms.
20

---

(a) Original Image
(b) SKLEARN GK
(C) OUR ALGORITHM
(d) SKLEARN N--N
(c) FAISS EXACT
(f) FAISS HNSW
FAISS IVF
(h) Original Image
(i) SKLEARN GK
(j) OUR ALGORITHM
(k) SKLEARN k-NN
(I) FAISS EXACT
(m) FAISS HNSW
(n) FAISS IVF
(0) Original Image
(p) SKLEARN GK
(q) OUR ALGORRITM
(r) SKLEARN İ-NN
FAISS EXACT
(U FAISS HNSW
(u) FAISS IVF
Figure 6: More examples on the performance of the spectral clustering algorithms for image segmen- lation.
21

---

(a) Original Image
(b) SKLEARN GK
(C) OUR ALGORITHM
(d) SKLEARN k-NN
(c) FAISS EXACT
(1) FAISS HNSW
FAISS IVF
(h) Original Image
(i) SKLEARN GK
() OUR ALGORITHM
(k) SKLEARN K-NN
()) FAISS EXACT
(m) FAISS HNSW
(n) FAISS IVF
Figure 7: More examples on the performance of the spectral clustering algorithms for image segmen- tation.
22


Reviews:
[{"review_id": "C6r9DxGY7B", "reviewer": "Reviewer_awKP", "summary": "Authors proposes a fast approximation method for constructing sparse similarity graphs from data for spectral clustering. They leverage Sun and Zanetti (2019)'s sampling-based approach to find \"cluster-preserving sparsifier\" for similarity graphs. This is combined with black box KDE and binary search algorithms to generate a fast & memory efficient approximation of SZ's algorithm. The experimental results are shown to be fast and accurate compared to Scikit-learn and FAISS's implementation.", "questions": "1. Given how it was used as a major motivation in this paper, I expected to see a comparison with the output of SZ algorithm in the experiments section. Have the authors compared SZ and proposed method experimentally for clustering?\n\n2. What do the calculated probabilities and sampled sparse graphs look like? Can the authors show any examples on a small graph or using a synthetic dataset? And how does the graph actually compare to that from SZ?\n\n3. It's missing a discussion on limitations of this work.", "limitations": "Not really.", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "This paper combines existing ideas (e.g. SZ algorithm, KDE algorithm, binary search) to create a novel approach to kernel graph construction that has potential to be used by many researchers in the community. The experimental results including runtime comparison and image segmentation (e.g. Figure 2) are quite impressive. It was also well written and clearly presented.", "weaknesses": "Please see questions below.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Many thanks for your very positive report.   Here are our responses to your questions.\n\n**Question 1.**\nIn Figure 1 of the attached PDF file, we provide an updated version of the experiment on the two moons dataset to include the results of the SZ algorithm.\nThe SZ algorithm iterates over every edge weight in the fully connected similarity graph in order to construct a sparsifier, and so its running time has a quadratic dependency on the number of data points.\n\n**Question 2.**\nFigure 2 of the attached PDF file illustrates the structure of the constructed graph on a small toy dataset with 100 data points.\nWe observe that both graphs preserve the cluster structure of a fully connected similarity graph.\nWe will include an illustration demonstrating the sparse graph structure in the final version of the paper.\n\n**Question 3.** We will add more discussion on the limitations of our work in the final version of the paper."}, {"review_id": "wv8OGne30j", "reviewer": "Reviewer_eUVC", "summary": "The paper proposed an efficient approximation of the similarity graph using the Kernel Density Estimation (KDE) method as a black box. The proposed framework can preserve the clustering structure and can be run in nearly linear time. Lastly, experiments on synthetic and image segmentation datasets are provided to support these claims. \n", "questions": "1. Related to Weakness #2, in practice, people usually use exponentially decaying kernel $k$ for the fully connected graph. You can approximate the fully connected graph with $\\varepsilon$-radius graph by choosing the radius $\\varepsilon$ and the bandwidth $\\sigma$ appropriately, e.g., if using the exponential kernel as per L132, you can choose $\\varepsilon = 3\\cdot \\sigma$; that way the far enough pairs will have small kernel values, i.e., $=\\exp(-9) \\sim 10^{-4}$, thus is negligible. I am curious if your method can be extended to the $\\varepsilon$-radius graph under this scenario as well?\n1. I am curious how difficult it is to extend this framework to the scenario of higher-order (neighborhood) simplicial complex and/or k-Laplacian. There are several works trying to use this information in data analysis [A-D], as well as algorithms [E] to sparsify the higher-order simplicial complexes with higher-order Cheerger constant. The reason I am asking is because the Higher-order k-Laplacian has dependency on $O(n_k)$ where $n_k$ is the cardinality of the k-complex (e.g., size of edges is $n_1 = |E|$), so if we can successfully sparsify the simplicial complex and/or k-Laplacian it can reduce the runtime drastically for those studies. \n1. How can this framework be extended to the manifold learning or diffusion map setting? In this scenario, we care not only about the approximation of the null space but also the first few eigne-values. Can we have similar guarantees for these scenarios?\n1. [Minor language usage] L21-23 (“Thanks to its out-performance over traditional clustering algorithms like k-means, this approach has…”) might be able to be improved by something like “Due to its superior performance compared to conventional clustering algorithms such as k-means, this approach has….”\n1. [Minor notational consistency] For the equation between L194-195, you might be able to improve the clarity by changing the $z \\in X_1$ to $x_j \\in X_1$ and the $x_j$ in summation to another variable. This is because you used $x_i$ and $x_j$ in the discussion above (L192). Similar to L13 in Algorithm 1, you might want to reconsider the use of $y_i$ there to make the notation more consistent.\n\n---\n[A] Dey, Tamal K., Jian Sun, and Yusu Wang. “Approximating Loops in a Shortest Homology Basis from Point Data.” In Proceedings of the Twenty-Sixth Annual Symposium on Computational Geometry, 166–75, 2010.\n\n[B] Chazal, Frédéric, Brittany Fasy, Fabrizio Lecci, Bertrand Michel, Alessandro Rinaldo, Alessandro Rinaldo, and Larry Wasserman. “Robust Topological Inference: Distance to a Measure and Kernel Distance.” The Journal of Machine Learning Research 18, no. 1 (2017): 5845–84.\n\n[C] Chen, Yu-Chia, and Marina Meila. “The Decomposition of the Higher-Order Homology Embedding Constructed from the k-Laplacian.” Advances in Neural Information Processing Systems 34 (2021).\n\n[D] Keros, Alexandros D., Vidit Nanda, and Kartic Subr. “Dist2cycle: A Simplicial Neural Network for Homology Localization.” In Proceedings of the AAAI Conference on Artificial Intelligence, 36:7133–42, 2022.\n\n[E] Osting, Braxton, Sourabh Palande, and Bei Wang. “Spectral Sparsification of Simplicial Complexes for Clustering and Label Propagation.” ArXiv:1708.08436 [Cs], February 1, 2019. http://arxiv.org/abs/1708.08436.\n\n", "limitations": "I believe that the authors have covered all the limitations of the algorithms (e.g., it is only cluster-preserving not spectral sparsification, runtime will depend on the implementation of $T_KDE$, etc.) in the current version, but adding it and summarize it in a specific section will be beneficial.\n\nNegative social impact is not applicable in this case, as this work primarily constitutes a theoretical contribution. \n", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "1. [Originality] The author proposed a novel algorithm by combining the ideas of subset/subspace dividing (somewhat related to LHS) with the help of probability measure estimated from KDE.\n1. [Clarity] An excellent written paper that provides clear insights into motivations and formulation. Great overview of how the algorithm works Figures 3-4 and Algorithms 1-2. The fast runtime is achieved by the $O(\\log(n))$ recursive call to the KDE as well as the superaddictive property of the $T_{KDE}$. \n1. [Quality] The authors compared their proposed algorithm with publicly available nearest neighbor implementations (FAISS, sklearn) and show drastic speedup in the large sample size regime. \n1. [Significance] I believe this paper makes a significant contribution in the unsupervised learning domains. The graph sparsification/Laplacian construction usually is done in a two-way approach, i.e., calculate pairwise distances and compute the similarity graph (with an additional sparsification step). The proposed method can do the aforementioned steps in a single step, thus drastically reducing the run time. \n", "weaknesses": "1. The authors have shown theoretically that the proposed algorithm is a {\\em cluster-preserving sparsifier}. This claim can be further supported by showing a comparison between the eigenvalues (or eigengaps) of the true, graph sparsifiers (e.g., from [13, 14, 19]), and the proposed method. \n1. There is another kind of similarity graph, called the $\\varepsilon$-radius graph. The graph $G(V, E)$ is constructed with $E = { (x_i, x_j) \\in V^2 :  \\|x_i - x_j\\| \\leq \\varepsilon}$. I think it will be informative to discuss this in the Introduction (L27-34).\n1. It is quite confusing to see Theorem 1 without the full definition, or even a high-level description, of what $\\rho$, $\\lambda$, $N_K$ are in Section 1. Readers can somewhat guess what they are based on notational convention, but to improve readability, I would suggest to either include the full definition or (even better) simplify it to make it less rigorous when it is first introduced in Section 1 (then provided a detailed version later in Section 4 or Appendix if no space).\n1. For the runtime of the FGT (L148), it might be more informative to write out all the terms inside the log rather than writing $\\widetilde O(m+n)$, because the error term $\\epsilon$ is hidden in the log.\n1. It will be interesting to see experiments on higher dimensional (real or synthetic) datasets and compare the runtime correspondingly. The experiments are run on (relatively) low-dimensional and straightforward examples, e.g., synthetic datasets (d=2) and image data (d=5). Some potential high-dimensional real datasets to consider e.g., [A,B].\n\n---\n[A] Mahmoud, Eman, Ali Takey, and Amin Shoukry. “Spectral Clustering for Optical Confirmation and Redshift Estimation of X-Ray Selected Galaxy Cluster Candidates in the SDSS Stripe 82.” Astronomy and Computing 16 (2016): 174–84.\n\n[B] Chmiela, Stefan, Alexandre Tkatchenko, Huziel E. Sauceda, Igor Poltavsky, Kristof T. Schütt, and Klaus-Robert Müller. “Machine Learning of Accurate Energy-Conserving Molecular Force Fields.” Science Advances 3, no. 5 (2017): e1603015.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Many thanks for your positive and detailed report, and valuable comments. We will take these into account when preparing the final version of our paper. Here are our responses to your questions and comments.\n\n**Question 1.**\nOur method constructs a sparsifier which has the same cluster structure as the fully-connected kernel graph.\nAs you point out, for an exponentially-decaying kernel the $\\epsilon$-radius graph will have a very similar structure to the kernel graph.\nHence our constructed graph could also be viewed as a cluster-preserving sparsifier of an $\\epsilon$-radius graph.\n\n**Questions 2 and 3.**\nThank you for the questions on generalising our approach to higher-order structures, and to manifold learning.\nThese are excellent directions for future work, but we have not considered these so far. \nWe will add some discussion in the final version of the paper.\n\n**Questions 4 and 5.**\nWe will make necessary changes based on your suggestions.\n\n**Weakness 1.** Classical algorithms for constructing spectral sparsifiers [13, 14, 19] are either based on a complicated graph decomposition framework, or Laplacian solvers. Despite of their nearly-linear time running time proven in theory, implementing these algorithms is a very challenging task. In fact, as far as we know, there is no publicly available implementation of nearly-linear time algorithms for spectral sparsifiers that work for a large-scale graph, and hence comparing the eigenvalues (or eigen-gap) between our constructed graph and spectral sparsifiers isn't feasible for a large-scale graph.  \n\nTaking this into account, we compare the eigen-gap of our sparsifier with the ones of a fully connected graph, and an SZ sparsifier. We study the two moon dataset, and in Table 1 in the attached PDF we report the eigengap $\\lambda_3 / \\lambda_2$, where $\\lambda_i$ is the $i$th smallest eigenvalue of the normalised Laplacian matrix of each constructed graph. We find that both the SZ algorithm and our proposed algorithm preserve a large eigen-gap between $\\lambda_2$ and $\\lambda_3$,  which implies that the cluster structure of the graph is preserved.\n\n**Weaknesses 2 - 4.**\nWe will incorporate these editorial suggestions into the final version of the paper.\n\n**Weakness 5.**\nWe will add some discussion on applying our algorithm to high-dimensional data in the final version of the paper.\n"}, {"review_id": "DO6Ewkvr1a", "reviewer": "Reviewer_hhPL", "summary": "This work presents an algorithm framework for creating a (weighted) sparse similarity graph, which is an important object for a variety of ML problems. The algorithm samples pairs of points to include in the graph using a sampler based on kernel density estimation (KDE). The idea is to recursively split the data into subsets, compute the KDE on each subset, and select a subset proportional to the density. This process is done for every point in the dataset (node in the graph), resulting in a collection of edges in the sparse graph.\n\nA nice feature of this method is that it provably preserves the spectral clustering properties of the graph. To my knowledge, this is the only sub-quadratic (in space / time) algorithm that does this. The method is also agnostic to the choice of KDE solver, leading to interesting practical algorithms. When paired with a classical KDE method (the Fast Gauss Transform), the method outperforms a number of reasonable baselines for graph construction + spectral clustering in low dimensions.", "questions": "1. Does T_KDE include preprocessing + query time? Most of the recent work in fast KDE requires a preprocessing algorithm (typically $O(N)$) followed by a query algorithm (typically much faster than $O(N)$). It looks like this is the case based on the T_KDE reported for FGT and the LSH method, but it wasn't clear from the definition. This is important because Algorithm 1 calls the KDE subroutine on different subsets of X.\n\n2. What is the state of the art result (using e.g. neural network models) on the BSDS image segmentation dataset? Is it close to the result for k-NN / similarity graph algorithms or is it much better?\n\n**Minor items:**\nTo get this to compile under MacOS arm64 (M1 chips), I did the following.\n1. Pass an up-to-date clang compiler to cmake (via -D CMAKE_CXX_COMPILER).\n2. Add the flag \"-I/usr/local/include\" to the compile_args in setup.py.\n3. Break the function declarations in stag_lib/KMeansRex/mersenneTwister2002.c into a header.\n4. pip install stag (stag isn't in requirements.txt, for some reason)\n\nNitpick on Line 92: \"Despite extensive research on LSH-based algorithms, it remains unclear whether such techniques can be employed to approximately construct a similarity graph in nearly-linear time.\"\nThis has been done at scale (albeit without theoretical guarantees on the quality of the similarity graph), see section 4 of [the FLASH paper](https://arxiv.org/pdf/1709.01190.pdf).", "limitations": "I do not foresee any negative societal impacts of this work. ", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 2, "contribution": 3, "strengths": "**S1: Novel and interesting reduction,** from KDE to similarity graph construction. The sampling mechanism (Algorithm 1) is particularly nice, and it likely has further applications to other settings where we want to sample proportional to kernel sums. The proposed framework also benefits from and motivates the growing body of work on the KDE problem. This leads to exciting possibilities for practical algorithms, depending on which KDE method is used in the framework.\n\n**S2: Good theoretical results.** It is great to see rigorous guarantees on the quality of the subsequent partitioning of the similarity graph as well as the time / space complexity. Many studies only perform the latter, as it is difficult to prove this kind of result.\n\n**S3: Well-engineered implementation,** with reproducible experiments. I was able to fully reproduce all of the experiments in this paper. I think the implementation will be of use to the community.", "weaknesses": "**W1: Notation.** The notation in Section 2 is understandable but would benefit from some polishing. For example, $w_G(u, v)$ is the edge weight function over node inputs and is defined only for those $(u,v)$ that are in the graph. But there is a $w_G(S, V\\setminus S)$ (with a different-script $G$ and set-of-node inputs), which is the sum of outbound edge weights from within $S$ to other parts of the graph. There is also a $w(u,v)$ used in the definition of the node degree, which isn't defined anywhere. A graph is defined as $G = (V, E)$ (with a separately-specified weight function) but the similarity graph is given as $K = (V, E, \\mathrm{weights})$. The kernel density sum KDE is introduced with a subscript $g_{[a,b]}$, but the analysis doesn't use this notation and later introduces $g_X$ as the KDE sum over a set of points $X$. This is all a bit confusing.\n\nOne possible fix is to define everything in terms of a weighted graph $G = (V, E, w)$. Here, $w$ is the weighting function from $(u,v)\\to \\mathbb{R}^{+}$ where $w(u,v) = $ edge weight if $(u,v) \\in G$ and 0 otherwise. This would give a much simpler definition for the adjacency matrix ($A_{ij} = w(v_i, v_j)$), the degree ($d(v) = \\sum_{u \\in V} w(v, u)$), and the outbound edge weight sum, which can be written as a unary function $w_o(S) = \\sum_{u \\in S, v \\not \\in S} w(u,v)$ to further reduce confusion with $w(u,v)$. Also meanings of $\\rho_K$ and $\\lambda$ should also be introduced before they are used in the main theorem, even if the full definition is deferred until later.\n\n**W2: Literature positioning / context.** The flexibility of the framework is a major strength of the paper, so it would be nice to highlight the (growing) number of KDE algorithms. For example:\n- Hashing-Based Estimators (HBE): Uses LSH tables to do KDE, where we only compute the KDE over the points that collide in the LSH table. The original paper was in [[FOCS 2017]](https://arxiv.org/abs/1808.10530), with a practical version in [[ICML 2019]](http://proceedings.mlr.press/v97/siminelakis19a) and some incremental progress in [[NeurIPS 2019]](https://openreview.net/forum?id=H1xEABrgIH). (only provably works for some kernels, fast in high dimensions, good in practice).\n- Interpolation KDE: [[AISTAT 2021]](http://proceedings.mlr.press/v130/turner21a/turner21a.pdf) performs KDE via polynomial interpolation (for arbitrary smooth kernels).\n- RACE: [[WWW 2020]](https://dl.acm.org/doi/fullHtml/10.1145/3366423.3380244) performs fast, online KDE, using a few MB of space (only for a limited set of kernels, but extremely space / time efficient). \n- Discrepancy coresets: [[COLT 2019]](https://proceedings.mlr.press/v99/karnin19a.html) performs KDE using online coresets (a re-weighted sub-sample of the dataset). In practice, this requires slower preprocessing than other methods but is more general / can be more accurate.\n\nEach of these methods will present unique tradeoffs when incorporated into this framework. While a full comparison of various KDE subroutines is probably out-of-scope, it would be nice to discuss some of these tradeoffs / make recommendations. A few similarity graph construction methods have also been developed since the review paper by Luxburg in 2007, which should be mentioned:\n\n- NN Descent: [[WWW 2011]](https://dl.acm.org/doi/abs/10.1145/1963405.1963487) This is an iterative technique that progressively refines a random guess for the initial k-NN graph.\n- FLASH: [[SIGMOD 2018]](https://dl.acm.org/doi/abs/10.1145/3183713.3196925) This is an LSH-based method that constructs a similarity graph (i.e. finds all pairs $(u,v)$ such that $w(u,v)$ exceeds a threshold). This is much faster than other LSH methods because it uses clever counting tricks to avoid doing any explicit distance / kernel calculations.\n\nFinally, the method in the paper, \"Learning space partitions for nearest neighbor search\" [[ICLR 2020]](https://openreview.net/forum?id=rkenmREFDr), might be relevant. This paper provides a framework for approximate near neighbor search that starts with a similarity graph, then partitions the data via a balanced min-cut partitioning of the graph. These partitions are then used to do a partition-based (FAISS IVF-style) similarity search, with provable guarantees (first theory for this type of algorithm). From a practical perspective, the hard part of this framework is to get the initial similarity graph, so your algorithm might help.\n\n\n**W3: Experiments.** A couple of highly competitive baselines are not represented in the experimental evaluation, notably FLASH (which is $O(N \\mathrm{polylog}(N))$) and NN-Descent (which scales about $O(N^{1.15})$, based on empirical results). For example, the C++ implementation of FLASH can construct a similarity graph for the webspam dataset (N = 350k) in under 10 seconds and for the friendster dataset (N = 65 million) in 1578 seconds. This 3-10x faster than the result of extrapolating Figure 5, and these datasets are of much higher dimension than two moons.\n\nThe datasets in the experiments are also all low-dimensional: even the BSDS task is only 5-dimensional. Many clustering applications (e.g. those that handle text embeddings) have inputs in the hundreds of dimensions. It would be nice to understand how this method performs in the higher-dimensional setting.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": " Many thanks for your positive and detailed report, and valuable comments. We will take these into account when preparing the final version of our paper. Here are our responses to your questions and comments.\n\n**Question 1: KDE preprocessing and query time.**\nIn our paper, $T_{\\mathrm{KDE}}$ includes both the preprocessing and query time of the KDE algorithm.\nAs you correctly point out, our algorithm applies the KDE algorithm with  different subsets of the data points, and the preprocessing step is required each time.\nWe will clarify this in the final version of the paper.\n\n**Question 2: State-of-the-art on the BSDS dataset.** The state-of-the-art for the (specific) BSDS dataset is recorded in [this paper](https://ieeexplore.ieee.org/document/5557884), and this achieves a Rand Index of 0.85. However, to the best of our knowledge, the Segment Anything Model (SAM) available [here](https://segment-anything.com/) presents the state-of-the-art for most image segmentation datasets. SAM has been trained on a dataset of 11 million images, and produces much better segmentation results than most unsupervised approaches. We emphasise that our algorithm is a general clustering algorithm, and is not developed specifically for image segmentation.\n\n**Weakness 1.**  We agree  with your comment. We'll take this into account, and improve our use of notation in the final version of the paper.\n\n**Weakness 2.** We're pleased to see that your consider the flexibility of our developed framework as a major strength of the paper. Following your suggestion, we will spend more effort discussing  recent work on KDE algorithms, and their tradeoffs when incorporated into our developed framework.  \n\n**Weakness 3.** We'll discuss the application of our algorithm on higher-dimensional datasets, and compare it against the algorithms (FLASH and NN-Descent) that you mentioned in the final version. Conducting such experimental studies is difficult at the moment, due to the limited time of the rebuttal stage.\n\n**Minor items.**\nWe are pleased that you were able to reproduce our experimental results with the provided code and for your comments on compiling and running the code.\nWe will update our README based on your feedback."}, {"review_id": "yLfXMpMgDP", "reviewer": "Reviewer_oQAs", "summary": "Given a set of n points and a pairwise kernel k(x,y), the fully connected similarity graph K is the weighted graph that has an edge between every pair x,y with weight k(x,y). The paper presents an algorithm for computing a sparse graph G that approximates K in a certain sense that is useful for spectral clustering. The algorithm is based on black-box calls to an efficient kernel density estimation (KDE) algorithm. For some kernels such black-boxes are available and have been the focus of much work recently, and they lead to efficient graph construction in this paper. The result is given in the form of a formal theorem and empirical evaluation.\n", "questions": "1. Please clarify connection to sampling in [6] (see above)\n\n2. The paragraph about LSH in lines 90-94 (related work) sounds odd. It says it's unclear how to use LSH for approximating geometric similarity graphs. However, much of the work you cite for efficient KDE builds directly on LSH (and consequently, so do some instantiations of your algorithm, as well as similar graph algorithms from prior work mentioned earlier).", "limitations": "No major concerns, though some limitations on the generality of the results may not be sufficiently discussed, as mentioned above in weaknesses.\n", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper is generally clearly written and contains nice and convincing theoretical and empirical findings that contribute to the literature on the topic and could be useful in practice. It is a nice application of the recent progress on efficient KDE to graph clustering. Even though it may not be as novel and general and advertised (see below), I still find it above the bar for acceptance.", "weaknesses": "(1) Novelty is somewhat limited, and seems overstated. The intro (lines 61-62) claims \"a novel connection between the KDE and the fast construction of similarity graphs\". This connection is quite well-known already and has been the focus of some recent works, including [1], [6] and \"Spectral Sparsification of Metrics and Kernels\" (Quanrud, SODA 2021). Anyway this connection between efficient KDE and kernel similarity graphs is very natural and straightforward, and it is unsurprising that all these works appeared shortly after the recent progress on efficient KDE. These works instantiate the connection in possibly different ways, but the general idea is similar. The novelty in this work seems to me restricted to the specific technical application to the SZ clustering framework, and not so much a conceptual point about kernels and graphs. \n\n(2) A more specific concern about novelty is that the main sampling method at the heart of the algorithm (lines 195-204, Algorithm 1) seems to have appeared already in [6] (see their \"sample random neighbor\" primitive). The work is cited only for a certain conditional hardness result, while its very similar algorithmic ideas are not mentioned. How do these techniques relate to each other?\n\n(3) The generality of the result itself is also overstated, particularly in the abstract that the result is \"applicable to arbitrary kernels\". I suppose this alludes to the fact that the algorithm uses KDE as a black-box, this is rather misleading, since there relatively few and specialized instances of known black-boxes of the type needed for this algorithm, and for a relatively limited set of kernels. For arbitrary kernels the black-box does not exist and the result of the paper is not applicable.\n\nThis affects the writing in the intro. Lines 35-40 formulate the main question as \"constructing a sparse graph that preserves the cluster structure\", focusing on sparsity and eschewing running time, possibly in order to circumvent said limitations on generality. However, this renders the main question moot, since trivially the full graph can be constructed and then spectrally sparsified. It would probably have been preferable to promptly discuss running times with their inevitable limitations.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Many thanks for your positive and detailed report, and valuable comments. We will take these into account when preparing the final version of our paper. Here are our responses to your questions and comments.\n\n**Question 1.**\nWe agree that there is some similarity between our algorithm and the one in [6].\nWe remark that our algorithm involves sampling a random neighbour of *every* vertex in a single procedure, whereas the one in [6] samples a neighbour of one vertex.\nDue to this difference, the analysis of our algorithm requires a careful consideration of each row of the recursion tree in Figure 4 with respect to the sets $S_{i, j}$ and $X_{i, j}$.\nFurthermore, our application of these samples to construct a cluster-preserving sparsifier is novel and requires a much more involved analysis than the analysis of SZ [21] due to the sampling method based on KDE.\n\n**Question 2.**\nWe will revise this paragraph in the final version.\n\n**Weakness 1.**\nWe agree that the novelty in our work is the application of efficient KDE to the construction of a cluster-preserving sparsifier of the complete kernel similarity graph.\nIt was not our intention to claim more than this, and following your comment we will make it clearer in the final version.\n\n**Weakness 2.**\nSee the answer to Question 1.\n\n**Weakness 3.**\nBased on our black-box reduction, our algorithm is applicable to any kernel with an efficient KDE algorithm.\nWe will clarify this and remove the claim that our algorithm is applicable to arbitrary kernels.\nAs you point out, our goal is to develop an algorithm with  fast running time for constructing a sparse similarity graph, and we will make this clearer in the introduction."}]

Author rebuttal:
We thank all reviewers for their positive and detailed reviews. We respond to their specific questions individually, with reference to some additional figures and table in the attached PDF.

Meta-review: This paper proposes an efficient approach for similarity graph construction based on kernel density estimation. The reviewers appreciate the theoretical soundness and high originality of the proposed approach. Since all the reviewers are positive for the paper, I recommend accepting the paper.
Decision: Accept (spotlight)