Subquadratic Algorithms for Kernel Matrices via Kernel Density Estimation
Authors: Ainesh Bakshi
Piotr Indyk
Praneeth Kacham
Sandeep Silwal
Samson Zhou
Venue: ICLR 2023 notable top 25% (2023)
Keywords: kernel density estimation
sublinear time algorithms

Abstract: Kernel matrices, as well as weighted graphs represented by them, are ubiquitous objects in machine learning, statistics and other related fields. The main drawback of using kernel methods (learning and inference using kernel matrices) is efficiency -- given $n$ input points, most kernel-based algorithms need to materialize the full $n \times n$ kernel matrix before performing any subsequent computation, thus incurring $\Omega(n^2)$ runtime. Breaking this quadratic barrier for various problems has therefore, been a subject of extensive research efforts. 

We break the quadratic barrier and obtain \emph{subquadratic} time  algorithms for several fundamental linear-algebraic and graph processing primitives, including approximating the top eigenvalue and eigenvector, spectral sparsification, solving linear systems, local clustering, low-rank approximation, arboricity estimation and counting weighted triangles. We build on the recently developed Kernel Density Estimation framework, which (after preprocessing in time subquadratic in $n$) can return estimates of row/column sums of the kernel matrix. In particular, we develop efficient reductions from \emph{weighted vertex} and \emph{weighted edge sampling} on kernel graphs, \emph{simulating random walks} on kernel graphs, and \emph{importance sampling} on matrices to Kernel Density Estimation and show that we can generate samples from these distributions in \emph{sublinear} (in the support of the distribution) time. Our reductions are the central ingredient in each of our applications and we believe they may be of independent interest. We empirically demonstrate the efficacy of our algorithms on low-rank approximation (LRA) and spectral sparsification, where we observe a $\textbf{9x}$ decrease in the number of kernel evaluations over baselines for LRA and a $\textbf{41x}$ reduction in the graph size for spectral sparsification.

Published as a conference paper at ICLR 2023
SUBQUADRATIC ALGORITHMS FOR KERNEL MATRI- CES VIA KERNEL DENSITY ESTIMATION
Ainesh Bakshi MIT einesh@mit.edu
Piotr Indyk MIT indykemit.eed
Praneeth Kacham CMU pkachamêcs cmu edu
Sandeep Silwal MIT si silwal@mit edu
Samson Zhou
UC Berkeley and Rice University
s aams onouu gmai com
ABSTRACT
Kernel matrices, as well as weighted graphs represented by them, are ubiquitous objects in machine learning, statistics and other related fields. The main drawback using kernel methods (learning and inference using kernel matrices) is effi- ciency- given n input points, most kernel- based algorithms need to materialize the full n x n kernel matrix before performing any subsequent computation, thus incurring (2n?) runtime. Breaking this quadratic barrier for various problems has therefore, been a subject of extensive research efforts.
We break the quadratic barrier and obtain subquadratic time algorithms for several fundamental linear-algebraic and graph processing primitives, including approx- imating the top eigenvalue and eigenvector, spectral sparsification, solving lin- ear systems, local clustering, low-rank approximation, arboricity estimation and counting weighted triangles. We build on the recently developed Kernel Density Estimation framework, which (after preprocessing in time subquadratic in n) can return estimates of row/column sums of the kernel matrix. In particular, we de- velop efficient reductions from weighted vertex and weighted edge sampling on kernel graphs, simulating random walks on kernel graphs, and importance pling on matrices to Kernel Density Estimation and show that we can generate samples from these distributions in sublinear (in the support of the distribution) time. Our reductions are the central ingredient in each of our applications and we believe they may be of independent interest. We empirically demonstrate the efficacy of our algorithms on low-rank approximation (LRA) and spectral sparsi- fication, where we observe a 9x decrease the number of kernel evaluations over baselines for LRA and a 41x reduction in the graph size for spectral sparsification.
Introduction
For kernel function & R2d x Rd + R and a set X of the n kernel matrix K are defined K,, k,,,,,,,. Alternatively, one can view X as the vertex set of complete weighted graph where the weights between points are defined by the kernel matrix K. Popular choices of kernel functions is include the Gaussian kernel, the Laplace kernel, exponential kernel, etc; see (Schölkopf et al., 2002; Shawe-Taylor et al., 2004; Hofmann al., 2008) for a comprehensive overview.
Inh C Id of n points, the entries
Despite their wide applicability, kernel methods suffer from drawbacks, one of the main being effi- ciency given n input points d dimensions, many kernel-based algorithms need to materialize the full kernel matrix K before performing the computation. For some problems this is unavoid- able, especially if hhih-precision results are required (Backurs et al., 2017). In this work, we show that we in fact break this Q(n?) barrier for several fundamental problems in numerical linear algebra and graph processing. We obtain algorithms that run in o(n²) time and scale inversely- proportional the smallest entry of the kernel matrix. This allows to skirt several known lower

---

Published as a conference paper at ICLR 2023
bounds, where the hard instances require the smallest kernel entry be polynomially small in n. Our parameterization in terms of the smallest entry is motivated by the fact in practice, the smallest kernel value is often fixed constant (March et al., 2015; Siminelakis et al., 2019; Backurs et al., 2019; 2021; Karppa et al., 2022). We build on recently developed fast approximate algorithms for Kernel Density Estimation (Charikar & Siminelakis, 2017; Backurs et al., 2018; Siminelakis al., 2019; Backurs et al., 2019; Charikar et al., 2020). Specifically, these papers present fast approximate data structures with the following functionality:
Definition 1.1 (Kernel Density Estimation (KDE) Queries). For a given dataset R of size kernel function k, and precision parameter > O, KDE data structure supports the following operation: given a query y E Rd,, return value xDEx (y) that lies in the interval [(1 - 8)2, (1+ e) z), where Exex k(s, v), assuming that k(3, v) 2 for all E X.
The performance of the state of the art algorithms for KDE also scales proportional to the smallest kernel value of the dataset (see Table 1). In short, after preprocessing time that is sub-quadratic (in n). KDE data structures use time sublinear in n to answer queries defined above. Note that for all of our kernels, k(x,y) y) I 1 for all inputs 2;
Instantiations of KDE queries. The query times depend the dimension d, accuracy e, and lower bound т. The parameter 8 is assumed to be a constant.
Table 1:
Type Gaussian 1 kxx y)
Preprocessing Time
Query Time
Reference
v
nd (7701774-001 nd (2r0.15001
(Charikar et al., 2020)
e e² 17340(1 d
Exponential
vil2
(Charikar et al., 2020)
Laplacian
|--v|
(Backurs et al., 2019)
22p0. 2
Rational Quadratic
(Backurs et al., 2018)
1.1 Our Results
We show that given KDE data structure as described above, is possible to solve a variety of matrix and graph problems in subquadratic time o(n?), i.e., sublinear in the matrix size. We emphasize that in our applications, we only require black-box access to KDE queries. Given this, we design algorithms for problems such eigenvallueeigeenvetoor estimation, low-rank approximation, graph sparsification, local clustering, aboricity estimation, and estimating the total weight of triangles.
Our results obtained via the following two-pronged approach. First, KDE data structures to design algorithms for the following basic primitives, frequently used in sublinear time algorithms and property testing:
sampling vertices by their (weighted) degree K (Theorems C.2 and C.4 and Algorithms 2 4), sampling random neighbors of a given vertex by edge weights and sampling a random weighted edge (Theorem C.5 and Algorithms and 6),
performing random walks in the graph (Theorem C.7 and Algorithm 7), and
4. sampling the rows of the edge-vertex incident matrix and the kernel matrix K, both with proba- bility proportional to respective row norms squared (Section D.I, Theorem D.I, and Section D.2, Corollary 1.10 respectively).
In the second step, we use these primitives to implement host of algorithms for the aforementioned problems. We emphasize that these primitives are used in black-box manner, meaning that any further improvements to their running times will automatically translate into improved algorithms for the downstream problems. For our applications, we make the following parameterization, which we expand upon in Remark B. and Section B. At high level, many of our applications, such as spectral sparsification, are succinetly characterized by the following parameterization.
Parameterization 1.1. All of our algorithms are parameterized by the smallest edge weight in the kernel matrix, i.e., the smallest edge weight in the matrix K is least
2

---

Published as a conference paper at ICLR 2023
Table 2: Summary of linear algebra and graph applications for KDE subroutines. We suppress dependence on the precision e. In spectral/local clustering and low-rank approximation, & denotes the number of clusters and the rank of the approximation desired, respectively. The parameter refers to the quality of the underlying clusters; see Section E.I.
Problem
of KDE Queries
Post-procsssing time
Prior Work
Spectral sparsification (Thm. 1.2)
Romark B.1
Luplacian system solver (Thm. 1.2)
Remark B.1
Low rank approx. (Thm. 1.5)
O(n)
O (n poly (k) nkd)
Remark B.3 ((n+
Eigenvalue Spectrum appeox. (Thm. 1.3)
Õ(17))
Old/r)
Approximating 1st Higenvalue (Thm. 1.4)
Remark B.2
d (-poly11/))
w(n) (Remark B.2)
Local clustering (Thm, 1.6)
δ (pnl((A) pr p3.
v
δ (poly(A) f 21.3)
Remark B.4
Spoctral clustering (Thm. E.7)
= + +Ō ((k)) 2
Romark B.4
Arboricity estimation (Thm. 1.8)
(2)
O(x²) O(n2d)
Triangle estimation (Thm. 1.9)
((n2d)
Our applications derived from the basic graph primitives above can be partitioned into two overlap- ping classes, linear-algebraic and graph theoretic results. Table 2 lists our applications along with the number of KDE queries required in addition to any post-processing time. We refer to the spe- cific sections of the body listed below for full details. We note that in all of our theorems below, we assume access to a KDE data structure of Definition 1.1 with parameters & and T.
One of our main results is spectral sparsification of the kernel matrix K interpreted as a weighted graph. In Section D.I, we compute sparse subgraph whose associated matrix closely approximates that of the kernel matrix K. The most meaningful matrix to study for such a sparsification is the Laplacian matrix, defined as D K where /) is a diagonal matrix of vertex degrees. The Laplacian matrix encodes fundaaeental combinatorial properties of the underlying graph and has been well- studied for numerous applications, including sparsification; see (Merris, 1994; Batson et al., 2013; Spielman, 2016) for survey of the Laplacian and its applications. Our result computes sparse graph, with a number of edges that is linear in n, whose Laplacian matrix spectrally approximates the Laplacian matrix of the original graph K under Parameterization 1.1.
Theorem 1.2 (Informal; see Thm. D.I). Let L be the Laplacian matrix corresponding to the graph K. Then, for any E E (0, 1), there exists an algorithm that outputs a weighted graph G' with only m O(n log (ogn((627)) edges, such that with probability at least 9/10, 10, (1-e) - L.g 3 (I+
The algorithm makes ((m) KDE queries and requires (md) post-processing time.
We compare our results with prior works in Remark B.I. We also show that Parameterization 1.1 is inherent for spectral sparsification. In particular, we use a hardness result from (Alman al., 2020) to show that for the Gaussian kernel, under the strong exponential time hypothesis (Impagliazzo & Paturi, 2001), any algorithm that returns an apparoiite spectral sparsifier with O(n 1999
edges requires Ω nn- 2log(1/7))
time (see Theorem D.4 for a formal statement). Obtaining the
optimal dependence on remains an outstanding open question, even for Gaussian and Laplace ker- nels. Spectral sparsification has further downstream apppattaatmailittttiinnni in solving Laplacian linear systems, which we present in Section 1.
Continuing the theme of the Laplacian matrix, in Section D.3, we also obtain a suceinet summary of the entire eigenvalue spectrum of the (normalized) Laplacian matrix using a total number of KDE queries independent of n, the size of the dataset. The error of the approximation is measured in terms of the earth mover distance (see Eq. (D.1)), or EMD, between the approximation and the true set of eigenvalues. Such a result has applications in determining whether an underlying graph can be modeled from specific graph generative process
al., 2018).
Theorem 1.3 (Informal; see Theorem D.1I). Let E E (0,1) be the error parameter and be the normalized Laplacian of the kernel graph K. Let X1 2 >2 > X, be the eigenvalues of L and let X be the resulting vector. Then, there exists an algorithm that uses Ö (exp (1/82) /T) KDE queries
3

---

Published as a conference paper at ICLR 2023
and exp (1/e2) d/T post-processing time and outputs vector À such that with probability 99/199 99/100, EMD 2, S.
Again to the best of our knowledge, all prior works for approximating the spectrum in EMD require constructing the full graph beforehand, and thus have runtime ((n=d). Next, we obtain truly sub- linear time algorithms for approximating the top eigenvalue and eigenvector of the kernel matrix, problem which was studied in (Backurs et al., 2021). Our result is the following theorem. Our bounds, and those of prior work, depend on the parameter p. which refers to the exponent of in the KDE query runtimes. For example for the Gaussian kernel, p z 0.173. See Table 1 for other kernels.
Theorem 1.4 (Informal; see Theorem D.15). Given an n x n kernel matrix K that admits a KDE data-structure with query time d/(=²r") (Table 1), there exists an algorithm that outputs unit vector u such that v² Kv 2 (1 - ()K(K) in time min (Ö(d/(s*57*)), ((dd/ 2t22))) where A (K) denotes the largest eigenvalue of K.
We discuss related works in Remark B.2. In summary. the best prior result of (Backurs et al., 2021) had a runtime of 22n +7) whereas our bound has no dependence on n. Finally, our last linear- algebraic result is an additive-error low-rank approximation of the kernel matrix, presented in Sec- tion D.2.
Theorem 1.5 (Informal; see Cor. D.10). There exists an algorithm that outputs a rank- r matrix B such that || BII < ||K - K2 K with probability 99% where K2 is the optimal rank-r approximation of K. It uses n KDE queries and O( poll((r 1/e) + +rd/e) post- post-processing time.
We give detailed comparisons between our results and prior work in Remark B.3. As a summary, (Bakshi et al., 2020b) obtain a relative error approximation with a running time of (nd (r/E) where denotes the matrix multiplication constant, whereas our running time is dominated by O(nr//s) and we obtain only additive error guarantees. Nevertheless, the algorithm we obtain, which builds upon the sampling scheme of (Frieze et al., 2004), is a conceptually simpler algorithm than the algorithm of (Bakshi et al., 2020b) and easier to empirically evaluate. Indeed, we implement this algorithm in Section 2 and show that it is highly competitive to the SVD.
We now move onto graph applications. We obtain an algorithm for local clustering, where we are asked whether two vertices belong to the same or different vertex communities. The notion of a cluster structure is based on the definition of k-clusterable graph, formally introduced in Defini tion E.3. Intuitively, describes a graph whose vertices can be partitioned into & disjoint clusters with high-connectivity within clusters and relatively sparse connectivity in in-between clusters.
Theorem 1.6 (Informal; see Theorem E.5). Let K be a k-clusterable kernel graph with clusters V = UigighV V, Let U, W be one of (not necessarily distinct) clusters V,. Let u, w be randomly chosen vertices in partitions and W with probability proportional to their degrees. There exists == c(EE k) and an algorithm that uses Öcc(k, e) /n/r//.)) KDE queries and post-processe with the property that with probability 1- U = W then the algorithm reports that and are in the same cluster and if U + 4W, the algorithm reports that u and w are in different clusters.
time,
Our definitions for the local clustering result are adopted from prior literature in property testing; see Remark B.4 for an overview of related works. Our sparsification result also automatically lends itself to an application in spectral clustering, an algorithm that clusters vertices based on the eigenvectors of the Laplacian matrix, which is outlined in Section E.2. We obtain an algorithm for approximately computing the top few eigenvectors of the Laplacian matrix, which is one of the main bottlenecks in spectral clustering in practice, with subquadratic runtime. These approximate eigenvectors are used to form the clusters.
Theorem 1.7 (Informal; Theorem E.8). Let L be the Laplacian matrix of the spectral sparsifier. There exists an algorithm that can compute (1 + appprooiimations of the first k eigenvectors of L in time Ô (km/ (x22)(7222.5))
We also give algorithms for approximating the arboricity of graph, which is the density of the densest subgraph of the kernel graph (see exact definition in Section E.3).
4

---

Published as a conference paper at ICLR 2023
Theorem 1.8 (Informal; see Theorem E.9). There exists algorithm that uses (E2((E)T) KDE queries and (mm) post-procsssing time and outputs a sparse subgraph G' of kernel graph such that with high probability, (1 =))ag S ag' S +1+ where ag is the arboricity of G.
To the best of our knowledge, all prior works on computing the arboricity require the entire graph to be known beforehand. In addition, computing the arboricity requires time ((m)) where is the number of edges leading to a runtime of Ö(z²) + O(r2d) (Gallo et al., 1989). In Section E.4, we also give an algorithm for approximating the total weight of all triangles of K, again interpreted as a weighted graph. We define weight of triangle as the product of its edge weights. This is natural definition if weighted edges are interpreted as parallel unweighted edges, in addition to having applications in defining cluster coefficients of weighted graphs (Kalna & Higham, 2006; Li et al., 2007; Antoniou & Tsompa, 2008). Our bound is similar in spirit to the bound of the unweighted case given in (Eden et al., 2017), under a different computation model. We refer to Remark B.5 for discussions related works.
Theorem 1.9 (Informal; see Theorem E. There exists algorithm that makes 0(1/7²) KDE queries and the same bound for post-processing time and with probability at least z. outputs a (1 = = approximation to the total weight of the triangles in the kernel graph.
On the other hand, there is a line of work that considers dimensionality reduction for kernel density estimation e.g., through coresets (Phillips& & Tai, 2018; 2020a; Tai, 2022). We view this direction of work as orthogonal to our line of study. Lastly, the work (Backurs et al., 2021) is similar in spirit to our work as they also utilize KDE queries to speed up algorithms for kernel matrices. Besides top eigenvalue estimation mentioned before, (Backurs et al., 2021) also study the problem of estimating the sum of all entries in the kernel matrix and obtain tight bounds for the latter.
1.2 Technical Overview
We provide a high-level overview and intuition for our algorithms. We first highlight our algorithmic building blocks for fundamental tasks and then describe how these components can be used to handle a wide range of problems. We note that our building blocks use KDE data structures in a black-box way and thus we describe their performance in terms of the number of queries to a KDE oracle. We also note that a permeating theme across all subsequent applications is that we want to perform some algorithmic task on kernel matrix K without computing each of its entries k(x,,
Algorithmic Building Blocks. We first describe the "multi-level" KDE data structure, which con- structs a KDE data structure on the entire input dataset X, and then recursively partitions X into two halves, building a KDE data structure on each half. The main observation here is that if the initialization of a KDE data structure uses runtime linear in the size n of X, then at each recursive level, the initialization of the KDE data structures across all partitions remains linear. Since there are Ologn) levels, the overall runtime to initialize our multi-level KDE data structure incurs only logarithmic overhead (see Figure 1 for illustration).
Weighted vertex sampling. We describe how to sample vertices approximately proportional to their weighted degree, where the weighted degree of vertex 2, with i E [r] is w. Σj k(z,, 2z). We observe that performing n KDE queries suffices to get approximation of the weighted vertex degree of all n vertices. We can thus think of vertex sampling as a preprocessing step that uses n queries upfront and then allows for arbitrary sample access at any point in the future with query cost. Moreover, this preprocessing step of taking n queries only needs to be performed once. Further, can then perform weighted vertex sampling from a distribution that is E-close in total variation to the true distribution (see Theoorm C.4 for details). Here, we use a multi-level tree structure to itera- tively choose subset of vertices with probability proportional to its approximate sum of weighted degrees determined by the preprocessing step, until the final vertex is sampled. Hence after the initial n KDE queries, each query only uses Ollog n) runtime, which significantly better than the naive implementation that uses quadratic time to compute the entire kernel matrix.
Weighted neighbor edge sampling. We describe how to perform weighted neighbor edge sampling for a given vertex æ. The goal of weighted neighbor edge sampling is to efficiently output a vertex (1+t))((),4))
for all k [n]. Unlike the degree case, edge sampling such that Pr(v- xt] Sielllayy k(a,a;)
5

---

Published as a conference paper at ICLR 2023
is not a straightforward KDE query since the sampling probability is proportional to the kernel value between two points, rather than the sum of multiple kernel values that KDE query provides. However, can utilize a similar tree procedure as in Figure 1 in conjunction with KDE queries.
In particular, consider the tree in Figure where each internal node corresponds to a subset of neighbors of x. The two children of parent node in the tree are simply the two approximately equal subsets whose union make up the subset representing the parent node. We can descend down the tree using the same probabilistic procedure pick one of the children to descend into with probability proportional to the sum of the edge weights represented by the children. The sum of edge weights of the children can be approximated by a query to an appropriate KDE data structure in the "multi-level" KDE data structure described previously. By appropriately decreasing the error of KDE data structures at each level of the tree, the sampled neighbor satisfies the aforementioned sampling guarantee. Since the tree has height O(log n), then we can perform weighted neighbor edge sampling, up to a tunably small total variation distance, using (log n) KDE queries and Ollg n) time (see theorems C.5 and C.6 for details).
in the vertex sampling case: at every node, we
Random walks. We use our edge sampling procedure to output a random walk on the kernel graph, where at any current vertex v of the walk, the next neighbor of visited by the random walk is chosen with probability proportional to the edge weights adjacent to v. In particular, for a random walk with T steps. we can simply sequentially call our edge sampling procedure T times, with each instance corresponding to a separate step in the random walk. Thus we can perform T steps of a random walk, again up to a tunably small total variation distance, using O(Tlog KDE queries and O(T log n) additional time.
Importance Sampling for the edge-vertex incidence matrix and the kernel matrix. We now describe how to sample the rows of the edge vertex incident matrix H and the kernel matrix K with probability proportional the importance sampling score / leverage score (see Definition D.2). We remark that approximately sampling proportional to the leverage score distribution for H is a fundamental algorithmic primitive spectral graph theory and numerical linear algebra. We note that apriori, such a task seems impossible to perform in o(n?) time, even if the leverage scores are precomputed for us, since the support of the distribution has size On²²). However, note we do not need to compute (even approximately) each leverage score to perform the sampling, but rather just output an edge proportional to the right distribution.
We accomplish this by instead sampling proportional to the squared Euclidean norm of the rows of H. It is known that oversampling the rows of a matrix by a factor that depends on the condition number is sufficient we show that H has a condition number (Lemma D.3) that is bounded by poly(1/7). Recall, the edge-vertex incident matrix is defined as the (2) x n matrix with the rows indexed by all possible edges and the columns indexed by vertices. For each e = {i, we have H(JJ). = / Vk(z., 23) and H(WJ) - (kz,, zy). We pick the ordering of i and j arbitrarily. Note that this is a weighted analogue of the standard edge-vertex incident matrix and satisfies HT H = L.g where Lg is the Laplacian matrix of the graph corresponding to the kernel matrix K. For both I and K, we wish to sample the rows with probability proportional to row normed squared. For example, the row re (I;, 2y) in H satisfies ||rel12 2k(2,, Since the squared norm of
approximate leverage score sampling (see proof of Theorem D.1). Further,
corresponding to edge each row is proportional to the weight of the corresponding edge, we can perform this sampling by combining the weighted vertex sampling and weighted neighbor edge sampling primitives: we first sample a vertex with probability proportional to its degree and then sample an appropriate random neighbor. Thus our row norm sampling procedure is sufficient simulate leverage score sampling (up to a condition number factor), which implies our downstream application of spectral sparsification.
We now describe the related primitive of sampling the rows of the kernel matrix K. Naively per- forming this sampling would require us to implicitly compute the entire kernel matrix, which as mentioned previously, is prohibitive. However, if there exists a constant e such that the kernel func- tion k that defines the matrix K satisfies k(x, y)² = k(cr, ey) for all inputs x. y. then the & norm of each row can be approximated via a KDE query on the transformed dataset X' = c.X In par- ticular, the L row norms of K are the vertex degrees of the kernel graph for X'. The property that k(x, y)² = k(cr, ey) holds for the most popular kernels such as the Laplacian, exponential, and Gaussian kernels. Thus, we can sample the rows of the kernel matrix with the desired probabilities.

---

Published as a conference paper at ICLR 2023
Linear Algebra Applications. We now discuss our linear algebra applications.
Spectral sparsification. Using the previously described primitives of weighted vertex sampling and weighted neighbor edge sampling, we show that an E spectral sparsifier for the kernel density graph G can be computed i.c., we compute a graph G' such that for all vectors x, (1 - E)x² LG S x. L I (1 + E)x² L, where Lg and Lo denote the Laplacian matrices of the graphs G and G'. Recall that H is the (2) x matrix such that Hi.,),.
Vk(wi, 2y)) and H(.s)
-V)(.,,,)) Here we use subsets of [n] of size 2 to index the rows of H and the entry to be made negative in the above definition is picked arbitrarily. It can be verified that H² H = LG. It is known that sampling/ 1. = (n log(n) (n) /e²) rows of the matrix H by using the so-called leverage scores gives atx x (2) selecting-and-sent
matrix S such that with probability least 9/10,
(1 = (1 (1- e) 3 H² S SH 3 (1 +e)HH H = (1 (1+
(1.1)
Thus the matrix SH directly corresponds to graph G', which is an E spectral sparsifier for graph G. The leverage scores of rows of H are also called "effective resistances" of edges of graph G. Unfortunately, with the edge and neighbor vertex sampling primitives that we have, we cannot per- form leverage score sampling of H. On the other hand, observe that the squared norm of row {i.]] of H is 2k 22, xy) and with application of vertex sampling and edge sampling, we can sample a row of H from the length squared distribution i.e., the distribution on rows where probability of sampling a row is proportional its squared norm. It is standard result that sampling from squared length distribution gives a selecting-and-scalin
matrix S that satisfies (1.1), although we have to sample 1 = Oxx²n (og(n)/=) rows from this distribution, where K - Tmax (H) (H)/o min (H) denotes the condition number of H (amax (H) (resp. (H)) denotes the largest (resp. smallest) positive singular values).
With the parameterization that for all i * j. k(I,, xz)) 2 T, we are able to show that x S O(1/~15). Importantly, our upper bound on the condition number is independent of the data dimension and number of input points. We obtain the upper bound on condition number by using Cheeger-type inequality for weighted graphs. Note that minn (H) 2 H)
(La() where we use A2( to denote the second smallest eigenvalue of a positive semidefinite matrix. Cheeger's in- equality lower bounds exactly the quantity A₂ A²(LG) in terms of graph conductance. A lower bound of t on every kernel value implics that every node in the Kernel Graph has a high weighted degree and this lets lower bound >z in terms of using a Cheeger-typ inequality from (Friedland & Nabben, 2002) and shows that On log(n)/>² samples from the approximate squared length sampling distribution gives an spectral sparsifier for the graph G.
First eigenvalue and eigenvector approximation. Our goal is to compute a 1 - & approximation to X, the first eigenvalue of K, and an accompanying approximate eigenvector. Such a task is key in kernel PCA and related methods. We begin by noting that under the natural constraint that each row of K sums must be at least nt by looking at the quadratic form associated with the all-ones vector.
least nt, condition used in prior works (Backurs et al., 2021), the first eigenvalue
Now we combine two disparate families of algorithms: first the guarantees of (Bhattacharjee al., 2021: Bakshi et al., 2020a) show that sub-sampling t x t principal submatrix of a PSD matrix preserves the eigenvalues of the matrix up to an additive O(n/v (t) factor. Since we've ve shown the first eigenvalue of K is at least nt. we can set & roughly (L1((=77)) with the guarantee that the top eigenvalue of the sub-sampled matrix is at lest (1 =)) Now we can either run the standard Krylov method algorithm (Musco & Musco, 2015) to compute the top eigenvalue of the sampled matrix or alternatively, we can instead use the algorithm of (Backurs et al., 2021), the prior state of the art, to compute the eigenvalues of the sampled matrix. At high level, their algorithm utilizes KDE queries to approximately perform power method on the kernel graph without creating the kernel matrix. In our case, we can instead run their algorithm on the smaller sampled dataset, which represents a smaller kernel matrix. Our final runtime is independent of the size of the dataset, whereas the prior state of the art result of (Backurs et al., 2021) 1) have a wn)) runtime.
Graph Applications. We now discuss our graph applications.
Local clustering. The random walks primitive allow us to run a well-studied local clustering al- gorithm on the kernel graph. The algorithm is quite standard in the property testing literature (see (Czumaj et al., 2015) and (Peng. 2020)) so we see our main contribution here as showing how the algorithm can be initialized for kernel matrices using our building blocks. At a high level, the goal
7

---

Published as a conference paper at ICLR 2023
of the algorithm is to determine if two input vertices u and belong to the same cluster of the kernel graph if the graph has a natural cluster structure (see Definition E.3 for the formal definition). The well-studied algorithm in literature performs approximately ((vr) random walks from u and e of a logarithmic length which is sufficient to estimate the distance between the endpoint distribution of the random walks. If the vertices belong to the same cluster, the distributions are close in lg distance which can be detected via a standard distribution tester of (Chan et al., 2014). The guarantees of the overall local clustering algorithm of (Czumaj et al., 2015) follow for kernel graphs since we only need to access the graph via random walks.
Arboricity estimation. The arboricity of weighted graph G (V, w) is defined as a == w((BGo)) Informally, the arborieity of a (weighted) graph represents the maximum maxugy |0] (weighted) density of a subgraph of G. To approximate the weighted arboricity, we adapt a re- sult of (McGregor et al., 2015), who observed that to estimate the arboricity on unweighted graphs, it suffices to sample a set of ((V)/=²) edges of G and computes the arboricity of the subsampled graph, after rescaling the weight of edges inversely proportional to their sampling probabilities.
We show that a similar idea works for estimating arboricity on weighted graphs. Although (MeGre- gor et al., 2015) showed that each edge should be sampled independently without replacement, we show that it suffices to sample a fixed number of edges with replacement. Moreover, we show that each edge should be one of the weighted edges with probability proportional to the weight of the edges, i.c., importance sampling. In fact, a similar result still holds if we only have upper bounds on the weight of each edge, provided that increase the number of fixed edges that sample by the gap between the upper bound and the actual weight of the edge. Thus, our arboricity algorithm requires sampling fixed number of edges, where each edge is sampled with probability propor- tional to some known upper bound on its weight. However for kernel density graphs, this is just our weighted edge sampling subroutine. Therefore, we achieve improved runtime over the naive ap- proach of querying each edge in the kernel graph by using our weighted edge sampling subroutine to sample fixed number of edges. Finally, we compute and output the arboricity of the subsampled graph as an approximation to the arboricity of the input graph.
2 Empirical Evaluation
We present empirical evaluations for our algorithms. We chose to evaluate algorithms for low-rank approximation (LRA) and spectral sparsification (and spectral clustering as a corollary) as they are arguably two of the most well studied examples in our applications and utilize wide variety of techniques present in our other examples of Sections D and E. Our evaluations serve as proof of concept that our queries which we constructed are efficient and easy to implement in practice. For our experiments, we use the Laplacian kernel k(x, v) = exp(-I<z- /ll1/o). A fast KDF implementation of this kernel exists due to (Backurs et al., 2019), which builds upon the techniques of (Charikar & Siminelakis, 2017). Note that the focus of our work is to use KDE queries in mostly black box fashion to solve important algorithmic problems for kernel matrices. This viewpoint has the important advantage that it is flexible the choice of any particular KDE query instantiation. We chose to work with the implementation of (Backurs et al., 2019) since it possesses theoretical guarantees, has an accessible implementation', and has been used in experiments in prior works such as (Backurs et al., 2019; 2021). However, we envision other choices of KDE queries, which maybe have practical benefits but are theoretically incomparable would also work well due to our flexibility.
Datasets. We use two real and two synthetic datasets in our experiments. The datasets used in the low-rank approximation experiments are MNIST (points in IR 784) (LeCun, 1998) and Glove word embeddings (points in 2000) (Pennington et al., 2014). We use 10 points from each of the test datasets. These datasets have been used in prior experimental works on kernel density estimation (Siminelakis et al., 2019; Backurs et al., 2019). The datasets experimental results for spectral sparsification and clustering are described in detail in F.
Evaluation Metrics. For LRA, we use the additive error algorithm detailed in Corollary D.10 of Section D.2. It requires sampling the rows of the kernel matrix according to squared row norms,
wwt/ta/oooaaaaaeeee from htt.ghtc..c///aalaaeeeeeeeaiiaaa http://wwwnee/efffiiiiccoo..
8

---

Published as a conference paper at ICLR 2023
which can be done via KDE queries outlined there. Once the (small) number of rows are sampled, we explicitly construct these rows using kernel evaluations. We compare the approximation error of this method computed via the standard frobenius norm error state of the art sketching algorithm for computing low-rank approximations, which is the input-sparsity time algorithm of Clarkson and Woodruff (Clarkson & Woodruff, 2013) (IS). We also compare to an iterative SVD solver (SVD). All linear algebra subroutines rely on Numpy, Scipy, and Numba implementations when applicable.
Low-rank approximation results. Note that the algorithm in Corollary D.10 has a O(k) depen- dence on the number of rows sampled. Concretely we sample 25k rows for a rank & approximation which we fix it for all experiments. For the MNIST dataset, the rank versus approximation error is shown in Figure 2a. The performance of our algorithm labeled as KDE is given by the blue curve while the orange curve represents the IS algorithm. The green curve represents the SVD error, which is a lower bound on the error for any algorithm. Note that for SVD calculations, we do not calculate the full SVD since that is computationally prohibitive; instead, we use an iterative solver. We can see that the errors of all three methods are comparable to each other. In terms of runtime, the KDE based method took 24.7 seconds on average for the rank 50 approximation whereas IS took 71.5 seconds and iterative SVD took 74.72 seconds
average. This represents 2.9x decrease in the running time. The time measured includes the time to initialize the data structures and matrices used for the respective algorithms. In terms of the number of kernel evaluations, both IS and iterative SVD require the kernel matrix, which is 10s kernel evaluations. On the other hand for the rank 50 approximation, our method required only 1.1 10² kernel evaluations, which is a 9x decrease in the number of evaluations. In terms of space, IS and iterative SVD require 10 floating point numbers stored due to initializing the full 10 x 10t matrix whereas our method only requires 10a 25 50 floating point numbers for the rank equal to 50 case and smaller for other. This is a 8x decrease in the space required. Lastly, we verify that we are indeed sampling from the correct distribution required by Corollary 1. In Figure 2b, we plot the points (æ;, Uh) where is the row norm squared for the ith row of the kernel matrix K and y/ is the row norm squared computed in our approximation algorithm (see Algorithm 9). As shown in Figure 2b, the data points fall very close to the y line indicating that our algorithm is indeed sampling from approximately the correct idcal distribution.
The qualitatively similar results for the Glove dataset are given in Figures 2c and 2d. For the glove dataset, the average time taken by the three algorithms were 37.7s, 37.7s, and 44.2s respectively, indicating that KDE and IS were comparable in runtime whereas SVD took slightly longer. How- ever, the number of kernel evaluations required by the latter two algorithms was significantly larger: for rank equal to 10, our algorithm only required 2.6 106 kernel evaluations while the other meth- ods both required 10s due to initializing the matrix. The space required by our algorithm was also smaller by factor of 40 since we only explicitly compute 25- 10 rows for the rank 10 case. For Glove, we only perform our experiments up to rank 10 since the iterative SVD failed to converge for higher ranks. While computing the full SVD avoids the convergence issue, it's computationally pro- hibitive. For example for MNIST, computing the full SVD of the kernel matrix took 552.9s, which is approximately an order of magnitude longer than any of the other methods.
Acknowledgements Ainesh Bakshi was supported by Ankur Moitra's ONR grant. Praneeth Kacham supported by National Institute of Health (NIH) grant 5401 HG 10798-2, Simons In- vestigator Award of David P. Woodruff, and Google as part of the "Research Collabs" " program. Piotr Indyk and Sandeep Silwal were supported by the NSF TRIPODS program (award DMS-2022448), Simons Investigator Award, MIT-IBM Watson AI Lab and NSF Graduate Research Fellowship un- der Grant No. 1745302. Work done in part while Samson Zhou was at Carnegie Mellon University and supported by Simons Investigator Award of David P. Woodruff and by the National Science Foundation under Grant No.
References
Mohammad Al Hasan and Vachik S Dave. Triangle counting in large networks: a review. Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery, 8(2):e1226, 2018.

---

Published as a conference paper at ICLR 2023
Josh Alman, Timothy Chu, Aaron Schild, and Zhao Song. Algorithms and hardness for linear algebra on geometric graphs. In 2020 JEEE 61st Annual Symposium on Foundations of Computer Science (FOCS), pp. 541-552. IEEE, 2020.
Ioannis E Antoniou and ET Tsompa. Statistical analysis of weighted networks. Discrete dynamics in Nature and Society, 2008, 2008.
Albert Atserias, Martin Grohe, and Daniel Marx. Size bounds and query plans for relational joins. In 49th Annual IEEE Symposium on Foundations of Computer Science, 2008.
Arturs Backurs, Piotr Indyk, and Ludwig Schmidt. On the fine- aue--grined complexity of empirical risk minimization: Kernel methods and neural networks. In Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems, pp. 4308- 4318, 2017.
Arturs Backurs, Moses Charikar, Piotr Indyk, and Paris Siminelakis. Efficient density evaluation for smooth kernels. 2018 IEFE 59th Annual Symposium on Foundations of Computer Science (FOCS), pp. 615-626, 2018.
Arturs Backurs, Piotr Indyk, and Tal Wagner. Space and time efficient kernel density estimation in high dimensions. In Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems, NeurlPS, pp. 15773-15782, 2019.
Arturs Backurs, Piotr Indyk, Cameron Musco, and Tal Wagner. Faster kernel matrix algebra via density estimation. In Proceedings of the 38th International Conference on Machine Learning, pp. 500-510, 2021.
Ainesh Bakshi and David Woodruff. Sublinear time low-rank approximation of distance matrices. Advances in Neural Information Processing Systems, 31, 2018.
Ainesh Bakshi, Nadiia Chepurko, and Rajesh Jayaram. Testing positive semi-definiteness via ran- dom submatrices. In 2020 IEEE 6/st Annual Symposium on Foundations of Computer Science (FOCS), pp. 1191-1202. IEEE, 2020a.
Ainesh Bakshi, Nadiia Chepurko, and David P Woodruff. Robust and sample optimal algorithms for psd low rank approximation. In 2020 IEEE 6/ Annual Symposium on Foundations of Computer Science (FOCS), pp. 506-516. IEEE, 2020b.
Joshua Batson, Daniel A Spielman, Nikhil Srivastava, and Shang-Hua Teng. Spectral sparsification of graphs: theory and algorithms. Communications of the ACM, 56(8):87-94, 2013.
Suman K. Bera and Amit Chakrabarti. Towards tighter space bounds for counting triangles and other substructures in graph streams. In Symposium on Theoretical Aspects of Computer Science (STACS 2017), 2017.
Rajarshi Bhattacharjee, Cameron Musco, and Archan Ray. Sublinear time eigenvalue approxima- tion via random sampling. CoRR, abs/2109.07647, 2021. URL https: https:/:///axxi http:::://axxv..rrogg/oo /bs/ 2109. 07 647.
Siu-on Chan, Ilias Diakonikolas, Paul Valiant, and Gregory Valiant. Optimal algorithms for test- ing closeness of discrete distributions. In Proceedings of the Twenty-Fitth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA, pp. 1193-1203, 2014.
Moses Charikar. Greedy approximation algorithms for finding dense components in graph. In Approximation Algorithms for Combinatorial Optimization, Third International Workshop, AP- PROX, Proceedings, pp. 84-95, 2000.
Moses Charikar and Paris Siminelakis. mashingg-ased-esttmatte for kernel density in high di- mensions. In 58th IEEE Annual Symposium on Foundations of Computer Science, FOCS, pp. 1032-1043, 2017.
Moses Charikar, Michael Kapralov, Navid Nouri, and Paris Siminelakis. Kernel density estima- tion through density constrained near neighbor search. 2020 IEEE 61st Annual Symposium on Foundations of Computer Science (FOCS), pp. 172-183, 2020.
10

---

Published as a conference paper at ICLR 2023
Justin Y. Chen, Talya Eden, Piotr Indyk, Honghao Lin, Shyam Narayanan, Ronitt Rubinfeld, Sandeep Silwal, Tal Wagner, David Woodruff, and Michael Zhang. Triangle and four cycle counting with predictions in graph streams. CoRR, abs/2203.095 2022.
Xue Chen and Eric Price. Condition number-free query and active learning of linear families. CoRR, abs/1711,10055 2017.
Ashish Chiplunkar, Michael Kapralov, Sanjeev Khanna, Aida Mousavifar, and Yuval Peres. Testing graph clusterability: Algorithms and lower bounds. In 2018 JEEE 59th Anmual Symposiun on Foundations of Computer Science (FOCS), pp. 497-508. IEEE, 2018.
Fan RK Chung and Fan Chung Graham. Spectral graph theory. American Mathematical Soc., 1997.
Kenneth L. Clarkson and David P. Woodruff. Low rank approximation and regression in input sparsity time. Symposium on Theory of Computing Conference, STOC, pp. 81-90, 2013.
David Cohen-Steiner, Weihao Kong, Christian Sohler, and Gregory Valiant. Approximating the spectrum of a graph. In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, KDD, pp. 1263-1271, 2018.
Artur Czumaj and Christian Sohler. Testing expansion bounded-degree graphs. Combinatorics, Probability and Computing, 19(5-6):693-709, 2010.
Artur Czumaj, Pan Peng, and Christian Sohler. Testing cluster structure of graphs. In Proceedings of the Forty-Seventh Annual ACM Symposium on Theory of Computing, STOC, pp- 723-732, 2015.
Tamal K Dey, Pan Peng, Alfred Rossi, and Anastasios Sidiropoulos. Spectral concentration and greedy k-clustering. Computational Geometry, 76: 76:19-32, 2019.
Talya Eden, Amit Levi, Dana Ron, and Seshadhri. Approximately counting triangles sublinear time. SIAM J. Comput., 46(5): 160--1666, 2017.
Brooke Foucault Welles, Anne Van Devender, and Noshir Contractor. "friend" a friend? inves- tigating the structure of friendship networks in virtual worlds. In CHIl The 28th Annual CHI Conference on Human Factors in Computing Systems, Conference Proceedings and Extended Abstracts, pp. 4027-4032, 2010.
Shmuel Friedland. Lower bounds for the first eigenvalue of certain m-matrices associated with graphs. Linear Algebra and its Applications, 172:71-80 84, 1992.
Shmuel Friedland and Reinhard Nabben. On Cheeger type inequalities for weighted graphs. Journal of Graph Theory, 41(1):1-17, 2002.
Alan Frieze, Ravi Kannan, and Santosh Vempala. Fast monte-carlo algorithms for finding low-rank approximations, Journal of the ACM (JACM), 5166: 2004.
Giorgio Gallo, Michael D. Grigoriadis, and Robert Endre Tarjan. A fast parametric maximum flow algorithm and applications. SIAM J. Comput., 18(1):30-55 55, 1989.
Grzegorz Gluch, Michael Kapralov, Silvio Lattanzi. Aida Mousavifar, and Christian Sohler. Spec- tral clustering oracles in sublinear time. Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms (SODA), pp. 1598-1617. SIAM, 2021.
Oded Goldreich. Introduction to property testing. Cambridge University Press, 2017.
Oded Goldreich and Dana Ron. On testing expansion in bounded-degree graphs. In Studies in Com- plexity and Cryptography. Miscellanea on the Interplay between Randomness and Computation, pp. 68 25-75. Springer, 2011.
Alexander G Gray and Andrew W Moore. N-body problems statistical learning. Advances in neural information processing systems, pp. 521-527, 2001.

---

Published as a conference paper at ICLR 2023
Alexander G Gray and Andrew W Moore. Nonparametric density estimation: Toward computational tractability, In Proceedings of the 2003 SIAM International Conference on Data Mining, pp. 203- 211. SIAM, 2003.
Thomas Hofmann, Bernhard Schölkopf, and Alexander J Smola. Kernel methods in machine learn- ing. The annals of statistics, 16(3):1171-1220, 2008.
Russell Impagliazzo and Ramamohan Paturi. On the complexity of k-sat. Journal of Computer and System Sciences, 62(2):367-31) 2001.
Piotr Indyk, Ali Vakilian, Tal Wagner, and David P. Woodruff. Sample-optimal low-rank approxi- mation of distance matrices. Conference on Learming Theory, COLT, pp. 1723-1751, 2019.
Satyen Kale and Seshadhri. Testing expansion in bounded degree graphs. 35th ICALP, pp. 527- 538, 2008
Gabriela Kalna and Desmond J Higham. Clustering coefficients for weighted networks. In Sympo- sium on network analysis in natural sciences and engineering, pp. 45, 2006.
Matti Karppa, Martin Aumüller, and Rasmus Pagh. Deann: Speeding up kernel-density estimation using approximate nearest neighbor search. In Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, pp. 3108-3137 2022.
Mihail N. Kolountzakis, Gary L. Miller, Richard Peng, and Charalampos E. Tsourakakis. Efficient triangle counting in large graphs via degree-basee based vertex partitioning. Lecture Notes in Computer Science, pp. 15-24, 2010.
Joannis Koutis, Gary L. Miller, and Richard Peng. A nearly-m log n time solver for sdd linear systems. In IEEE 52nd Annual Symposium on Foundations of Computer Science, pp. 590-598, 2011.
Tsz Chiu Kwok, Lap Chi Lau, Yin Tat Lee, Shayan Oveis Gharan, and Luca Trevisan. Improved cheeger's inequality: analysis of spectral partitioning algorithms through higher order spectral gap. In Symposium on Theory of Computing Conference, STOC'13 pp. 11-20, 2013.
Yann LeCun. The mnist database of handwritten digits. http://yann. lecun. com/exdd/mmistt/ 1998.
Dongryeol Lee and Alexander Gray. Fast high-dimensional kernel summations using the monte carlo multipole method. Advances in Neural Information Processing Systems, 21:929 2008
Dongryeol Lee, Andrew W Moore, and Alexander G Gray. Dual-tree fast gauss transforms. In Advances in Neural Information Processing Systems, pp. 747-754, 2006.
James R. Lee, Shayan Oveis Gharan, and Luca Trevisan. Multi-way spectral partitioning and higher- order cheeger inequalities. In Proceedings of the 44th Symposium on Theory of Computing Con- ference, STOC, pp. 1117-1130, 2012.
Victor E. Lee, Ning Ruan, Ruoming Jin, and Charu C. Aggarwal. A survey of algorithms for dense subgraph discovery. In Managing and Mining Graph Data, volume 40 of Advances in Database Systems, pp. 303-336. Springer, 2010.
Jure Leskovec, Lars Backstrom, Ravi Kumar, and Andrew Tomkins. Microscopic evolution of social networks. In Proceedings of the 14th ACM SIGKDD International Conference Knowledge Discovery and Data Mining, pp. 462 (62-470, 2008
Wenyuan Li, Yongjing Lin, and Ying Liu. The structure of weighted small-world networks. Physica A: Statistical Mechanics and its Applications, 376:708-716 2007.
Anand Louis, Prasad Raghavendra. Prasad Tetali, and Santosh S. Vempala. Many sparse cuts via higher eigenvalues. In Proceedings of the 44th Symposium on Theory of Computing Conference, STOC, pp. 1131-1140, 2012.
William B March, Bo Xiao, and George Biros. Askit: Approximate skeletonization kemel- independent treecode in high dimensions. SIAM Journal on Scientific Computing. 37(2):A1089 A1110, 2015.
12

---

Jeff M Phillips and Wai Ming Tai. Near-optimal coresets of kernel density estimates. Discrete & Published as a conference paper at ICLR 2023
Andrew McGregor, David Tench, Sofya Vorotnikova, and Hoa T. Vu. Densest subgraph in dynamic graph streams. In Mathematical Foundations of Computer Science 2015 40th International Symposium, MFCS, Proceedings, Part II. pp. 472-482, 2015.
Russell Merris. Laplacian matrices of graphs: survey. Linear algebra and its applications, 197: 143-176, 1994.
R. Milo, S. Shen-Orr, S. Itzkovitz, N. Kashtan, D. Chklovskii, and U. Alon. Network motifs: Simple building blocks of complex networks. Science, 298(5594):824 827, 2002.
Vlad I Morariu, Balaji Vasan Srinivasan, Vikas C Raykar, Ramani Duraiswami, and Larry S Davis. Automatic online tuning for fast gaussian summation. In NIPS, pp. 1113-1120, 2008.
Cameron Musco and Christopher Musco. Randomized block krylov methods for stronger and faster approximate singular value decomposition. In Advances in Neural Information Processing Sys- tems 28: Annual Conference on Neural Information Processing Systems, pp. 1396-1404, 2015.
Cameron Musco and David P Woodruff. Sublinear time low-rank approximation positive semidef- inite matrices. In 2017 IEEE 58th Annual Symposium Foundations of Computer Science (FOCS), pp. 672-683. 683. IEEE, 2017.
Pan Peng. Robust clustering oracle and local reconstructor of cluster structure graphs In Pro- ceedings of the 2020 ACM-SIAM Symposium Discrete Algorithms, SODA, pp. 2953-2972. SIAM, 2020. doi: 11177/997811111599911111 URL httph::p///oott org 1137/1.
9781611975994.1
Jeffrey Pennington, Richard Socher, and Christopher 1) Manning. Glove: Global vectors for word representation. In Proceedings of the 2014 conference on empirical methods in natural language processing (EMNLP). pp. 1532-1543, 2014.
Jeff M Phillips. e-samples for kernels. In Proceedings of the twenty-fourth annual ACM-SIAM symposium on Discrete algorithms, pp. 1622-1632. SIAM, 2013.
Jeff M. Phillips and Wai Ming Tai. Improved coresets for kernel density estimates. In Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA, pp. 2718-2727, 2018.
Jeff M. Phillips and Wai Ming Tai. Near-optimal coresets of kernel density estimates. Comput. Geom., 63(4))867-887, 2020a.
Discret.
Computational Geometry, 63(4))867-887, 2020b.
Kent Quanrud. Spectral sparsification metrics and kernels. In Proceedings of 2021 ACM-SIAM Symposium on Discrete Algorithms, SODA, pp. 1445-1464, 2021.
Parikshit Ram, Dongryeol Lee, William March, and Alexander Gray. Linear-time algorithms for pairwise statistical problems. Advances in Neural Information Processing Systems, 22:1527- 1535, 2009.
Bernhard Schölkopf, Alexander J Smola, Francis Bach, al. Learning with kernels: support vector machines, regularization, optimization, and beyond. MIT press, 2002.
C. Seshadhri, Ali Pinar, and Tamara G. Kolda. Fast triangle counting through wedge sampling. In the International Conference on Data Mining (ICDM), 2013.
John Shawe-Taylor, Nello Cristianini, et al. Kernel methods for pattern analysis. Cambridge uni- versity press, 2004.
Paris Siminclakis, Kexin Rong, Peter Bailis, Moses Charikar, and Philip Alexander Levis. Rehashing kernel evaluation in high dimensions. In Proceedings of the 36th International Conference Machine Learning, ICML, pp. 5789-5798, 2019.
Daniel A Spielman. The laplacian matrices of graphs. In Plenary Talk, IEEE Intern. Symp. Inf. Theory (ISTT),, 2016.
13

---

Published as a conference paper at ICLR 2023
Daniel A Spielman and Nikhil Srivastava. Graph sparsification by effective resistances. SIAM Journal on Computing, 40(6):1913-1926, 2011.
Daniel A. Spielman and Shang-Hua Teng. Nearly-linear time algorithms for graph partitioning, graph sparsification, and solving linear systems. In Proceedings of the 36th Annual ACM Sympo- sium on Theory of Computing, pp. 81-90, 2004.
Wai Ming Tai. Optimal coreset for gaussian kernel density estimation. In 38th International Sym- posium on Computational Geometry, SoCG, pp. 63:1-63:15, 2022.
Ulrike Von Luxburg. A tutorial on spectral clustering. Statistics and computing, 17(4)339-446, 2007.
David P. Woodruff. Sketching as a tool for numerical linear algebra. Foundations and Trends@) in Theoretical Computer Science, 10(1-2):1-157, 2014.
Yan Zheng, Jeffrey Jestes, Jeff M Phillips, and Feifei Li. Quality and efficiency for kernel density estimates in large data. In Proceedings of the 2013 ACM SIGMOD International Conference on Management Data, pp. 433-444, 2013.
A Omitted Technical Overview
We provide technical overview for applications whose discussions were omitted from the main text.
Kernel matrix low-rank approximation. this setting, our goal here is to output matrix B such that
K - BB
K - +
+ =||
where K, is the best rank-t approximation to the kernel matrix K. The efficient algorithm of (Frieze et al., 2004) is able to achieve this guarantee one can sample the ith row r; of with probability 2 2(1) ||r:112/ || | 1| We can perform such action using primitive, which is capable of sampling the rows K with probability proportional to the squared row norms for the Lapla- cian, exponential, and Gaussian kernels. Thus for these kernels, algorithms for computing a low-rank approximation.
immediately obtain efficient
Spectrum approximation. For this problem, the goal to compute approximations all the eigenvalues of the normalized Laplacian matrix of the kernel graph such that the error between the approximations and the true set of eigenvalues has small error in the earth mover metric. The algorithm of (Cohen-Steiner al., 2018) achieves this guarantee in time independent in the graph size given the ability to perform random walks on uniformly sampled vertices. Surprisingly, the number of random walks and their length does not depend on the number of vertices. Thus given our random walk primitive, we can efficiently simulate the algorithm of (Cohen-Steiner et al., 2018) on kernel graphs in a black-box manner.
Spectral clustering. Given our spectral sparsification result, we can immediately obtain a fast version of heuristic algorithm used in practice for graph clustering: we embed each vertex into R* using k eigenvectors of the Laplacian matrix and run k-means clustering. Clearly we have a sparse graph, the eigenvector computation is faster. Theoretically. we show that spectral sparsification preserves a notion clusterability which is weaker than the definition used the local clustering section and additionally give empirical evidence of the validity of this procedure.
14

---

Published as a conference paper at ICLR 2023
Weighted triangle estimation. We define the weight of a triangle the product of its edges, gen- eralizing the case were the edges have integer lengths, so that an edge can be thought of as multiple parallel edges. Under this definition, we adapt an algorithm from (Eden et al., 2017), who considered the problem in unweighted graphs given query access to the underlying graph. Specifically, we show that it suffices to sample a "small" set R of edges uniformly at random and then estimate the total weight of triangles including the edges of R under some predetermined ordering. In particular, the procedure of estimating the total weight of triangles including the edges e of It involves sampling neighbors of the vertices of e, which we can efficiently implement using our weighted neighbor edge sampling subroutine.
B Further Related Works
Remark B.1. Spectral sparsification for kernel graphs has also been studied in prior works, no- tably in (Alman al., 2020) and (Quanrud, 2021). We first compare to (Alman al., 2020), who obtain spectral sparsification using an entirely different approach. They obtain an almost linear time sparsifier (n (n)+)(1)) when the kernel is multiplicativily Lipschitz (see Section 1.1.2 in (Alman al., 2020) for definition) and show hardness for constructing such a sparsifier when it is not. Focusing on the Gaussian kernel, under Parameterization 1.1, (Alman et al., 2020) obtain an al- gorithm that runs in time O (nd + (22/log(1/r) log n log(log /E²) whereas our algorithm runs in O (nd (nd)oog(n)(((++))) (e2)(((7220173+0(())
time. We also note that the dimension d can be upper bounded
by Ollog n//e)) by applying Johnson--Lindenstrauss to the initial dataset. Therefore, (Alman et al., 2020) obtain a better dependence on 1/T, whereas we obtain a better dependence on n. A similar comparison can be established for other kernels as well. In practice, is set to be small fixed con- stant, whereas n can be arbitrarily large. Indeed in practice, a common setting of T is 0.01 or 0.001, irrespective of the size of the dataset (March et al., 2015; Siminelakis et al., 2019; Backurs et al., 2019; 2021; Karppa et al., 2022).
We now compare guarantees to that of (Quanrud, 2021). The author studies spectral sparsifi- cation resurrected to smooth kernels (for example kernels of the form 1/1| - vl12 which have a polynomial decay; see (Quanrud, 2021) for formal definition). This family does not include Gaus- sian, Laplacian, or exponential kernels. For smooth kernels, (Quanrud, 2021) obtained a spansifier with a nearly optimal (n//²²) number of edges in time ((nd/=²). Our algorithm obtains a similar dependence in n2, d, e but includes an additional 1/73 factor. However, it generalizes for any kernel supporting& a KDE data structure, which includes smooth kernels (Backurs et al., 2018) (see Table for summary of kernels where our results apply). Our techniques are also different: (Quanrud, 2021) does not use KDE data structures in a black-box manner to compute the sparsification as we do. Rather, they simulate importance sampling on the edges of the kernel graph directly. In addition to the nearly linear sparisifier, another interesting feature of (Quanrud, 2021) is that it enriches the connections between spectral sparsification of kernel graphs and KDE data structures. Indeed, the data structures used in (Quanrud, 2021) are inspired by and were used in the prior work of (Backurs et al., 2018) to create KDE query data structures themselves. Furthermore, the paper demonstrates how to instantiate KDE data structures for smooth kernels using the kernel graph sparsifier itself. We refer to (Quanrud, 2021) for details.
Remark B.2. We remark that our algorithm returns a spurse vector e supported on roughly (())(=?7)) coordinates. The best prior result is that of (Backurs et al., 2021) which presented an algorithm with total runtime O dn (log(n/))3
In comparison, our bound has no dependence on n and is thus a truly sublinear runtime. Note that the bound of (Backurs et al., 2021) does not depend on т. We do not state the number of KDE queries used explicitly in Table since our algorithm uses KDE queries on subsampled dataset and in addition. only uses them by calling the algorithm of (Backurs et al., 2021) as a subroutine (on the subsampled dataset). The algorithm of (Backurs et al., 2021) uses Ö(1/s) KDE queries but with various different initialization of so it is not meaningful to state "one" bound for the number of KDE queries used and thus the final runtime is more meaningful quantity to state. Lastly, the authors in (Backurs et al., 2021) present a lower bound of Q(nd) for estimating the top eigenvalue d1 which ostensibly seems odds with our stated bound which has dependence on n. However,
15

---

Published as a conference paper at ICLR 2023
the lower bound presented in (Backurs al., 2021) essentially sets 1/ poly(n) for a large polynomial factor depending on n (we estimate this factor to be )((²)). Since we parameterize our dependence via which in practice is often set fixed constant, can bypass the lower bound.
Remark B.3. We now compare our low-rank approximation result with recent work of (Musco & Woodruff, 2017: Bakshi al., 2020b). They showed the following theorem:
Theorem B.1 (Theorem 4.2, (Bakshi et al., 2020b)). Given X PSD matrix A, target rank E [n] and accuracy parameter E (0, 1), there exists algoritim that queries Ô (nr/s) entries in A and with probability at least 9/ /10, outputs a rank-r matrix B such that
A
B Bl I (1 (1+ 1A- Ar
where Ar is the best rank-r approximation to A. Further, the running time is Ō ( (n(r//))=)) where is the matrix multiplication constant.
We note that their result applies to kernel matrices as well via the following fact. Fact B.2 (Kernel Matrices are PSD, (Schölkopf al., 2002)). Let k be a reproducing kernel and X be n data points in Rd. Let K be the associated n x n kernel matrix such that Kij, = k(T5, Then, K> 0.
Here, the family of reproducing kerness is quite broad and includes polynomial kemelss Gaussian, and Laplacian kernel, among others. Therefore, their theorem immediately implies a relative error low-rank approximation algorithm for kernel matrices. Our result and the theorem of (Bakshi et al., 2020b) have comparable runtimes. While (Bakshi et al., 2020b) obtain relative-error guarantees, we only obtain additive-error guarantees.
However, reading each entry of the kernel matrix require (dd) time and thus (Bakshi et al., 2020b) obtain an running time of Ô (nd (r/s)==t). whereas our running time is dominated by (nnd//e). We note that similar ideas as our algorithm for additive error LRA were previously used to de- sign subquadratic algorithms running in time o(n²) for low-rank approximation of distance matri- ces (Bakshi & Woodruff, 2018; Indyk et al., 2019).
Remark B.4. Our definitions for the local clustering result are adopted from prior literature in property testing; see (Kale Seshadhri, 2008: Czumaj & Sohler. 2010: Goldreich Ron, 2011; Czumaj et al., 2015; Chiplunkar et al., 2018; Dey et al., 2019; Peng. 2020; Gluch et al., 2021) and the references within. Our algorithmic details for the local cluster section are also derived from prior works, such as the works of (Czumaj et al., 2015) and (Peng, 2020): indeed, many of the lemmas of the local clustering section follow in a straightforward fashion from (Czumaj et al., 2015) and (Peng, 2020). However, the key difference between these works and our work is that they are in the property testing model where one assumes access to various graph queries in order to design sublinear graph algorithms. To the best of our knowledge, implementation of prior works local clustering requires having access to the entire neighbor a vertex when performing random walks, thereby implying the runtime of ((nd) per step of the walk. In contrast, we give efficient constructions of these commonly assumed queries for kernel graphs, rather than assuming oracle access. Indeed, the fact that one can easily take existing algorithms which hold in non kernel settings and apply them to kernel settings in straightforward manner via our queries can be viewed a major strength of our work.
Remark B.5. Our general bound of the number of KDE qeuries required to approximate the total weight of triangles in Theorem E. is Ō (mm/wa//u wr), where we is the sum of all entries of K and WT is the total weight of triangles we wish to approximate. This bound is a natural generalization of the result of (Eden et al., 2017). There, the goal is to approximate the total number of triangles in an unweighted graph given access to queries of an underlying graph in the form of random vertices and random neighbors of given vertex (assuming the entire graph is stored in memory). While their model differs from our work, we note that KDE queries constructed in Section C play a similar role to the queries used in (Eden et al., 2017). There the authors give bound of (m23/2 /2 /T) queries where T is the total number of triangles. In our case, we indeed get bound of the order of m33 in the numerator as we < mWmaa and WT is the natural analogue of T' in (Eden et al., 2017). Finally note that under our parameterization of every edge in the kernel graph possessing weight at most and at least 7, our bound reduces to at most ((1/73) KDE queries.
16

---

Published as a conference paper at ICLR 2023
We finally note that to the best of our knowledge, all prior works for approximating the number of triangles in a graph require the full graph to be instantiated, which implies a lower bound of time ((ndd) in our setting.
We also note that our paper is closely related to the field of (graph) property testing. In graph prop- erty testing, it is customary to assume query access to an unknown graph via vertex and edge queries (Goldreich, 2017). While specific details vary, common queries include access to random vertices and random neighbors of a given vertex, among others. The goal of the field is to design algorithms that require queries sublinear in n, the number of vertices, or n², the size of the graph. We can interpret the graph primitives we construct as realization of the property testing model where queries are explicitly constructed.
B.1 Preliminaries
First, we discuss the cost of constructing KDE data structure and performing the queries described in Definition 1.1. Table 1 summarizes previous work on kernel density estimation though for the sake of uniformity, we list only "high-dimensional" data structures, whose running times are poly- nomial in the dimension d. Those data structures have construction times of the form O(dn/(rPe")) and answer KDE queries in time Odd)(r?=²), under the condition that for all queries y have 1 ExEEx k(z y) 2 (which clearly holds under our Parameterization 1.1). The algorithms are ran- domized, and report correct answers with a constant probability. The values of p lie in the interval [0, 1). and depend on the kernel. For comparison, note that a simple random sampling approach, which selects random subset R C X of size ((1/(re²) and reports [R] EvER k(x, y), achieves the exponent of p = 1 for any kernel whose values lie in [0, 1].
We view our algorithms as parameterized in terms of т, the smallest edge length. We argue this is a natural parameterization. When picking a kernel function k, we also have to pick a scale term a (for example, the exponential kernel is of the form k(z, y) cxp(-II> 1|x 3/l2/or)). common choice of a follows the so called 'median rule' where a is set to be the median distance among all pairs of points in X. Thus, according to the median rule, the 'typical' kernel values the graph K are (11). While this is only true for 'typical,' and not all, edge weights in K, we believe the KDE query abstraction of Definition 1.1 still provides nontrivial and useful algorithms for working with kernel graphs. Typically in practice, the setting of T is small constant, independent of the size of the dataset (Karppa et al., 2022).
practice,
We note that, in addition to the aforementioned algorithms with theoretical guarantees, there other practical algorithms based on random sampling, space partition trees (Gray & Moore, 2001; 2003; Lee et al., 2006; Lee & Gray, 2008; Morariu al., 2008; Ram al., 2009; March et al., 2015), coresets (Phillips, 2013; Zheng et al., 2013; Phillips & Tai, 2020b), or combinations of these methods (Karppa et al., 2022), which support queries needed in Definition 1.1; see (Karppa al., 2022) for an in-depth discussion on applied works.
While these algorithms do not necessarily have as strong theoretical guarantees as the ones discussed above and in Table 1, we can nonetheless use them via black box access in our algorithms and utilize their practical benefits.
C Algorithmic Building Blocks
C.1 Multi-level KDE
We first describe the "multi-level" KDE data structure, which required in our algorithms. The data structure recursively constructs KDE data structure the entire dataset X, and then recursively partitions X into two halves, building a KDE data structure on each half. See Algorithm for more details.
17

---

Published as a conference paper at ICLR 2023
A
11/122
1n/2+1,0
A1/4
Figure I: Multi-level Kernel Density Estimation Data Structure.
Algorithm Multi-level KDE Construction
Require: Input dataset X Rd precision
1: T' X
2: while [7] > do
3: Construct KDEx queries
4: Recursively apply Multi-level KDE Construction T[1 [m/2]] and T[[m/2] + m]
5: end while
6: Return all the data structures associated with the KDE query constructions
Lemma C.1. Given a dataset X C Rdd, suppose the initialization of the KDE data structure defined Definition 1.1 uses runtime f(n, E) for some function linear in n. Then the total construction time of Algorithm is (n E).
Proof. The proof follows from the fact that at each recursive level, we do O(f(n.))) total work since f is linear in and there (loo n)) levels.
C.2 Weighted Vertex Sampling
We discuss fundamental primitives. The first one is sampling vertices by their (weighted)
degree.
Definition C.1 (Weighted Vertex Sampling). The weighted degree of vertex x, with i [n] is w; = Ejzt kz),, The goal of weighted vertex sampling to output a vertex v such that Pr(u (1+e)ws, x.] Eyetn my for all i E [n].
This is straightforward application of using KDE queries to get the (weighted) vertex degree of all n vertices. Note that this only takes n queries and only has to be done once. Therefore, we can think of vertex sampling as preprocessing step that uses (nz) queries upfront and then allows for arbitrary access at any point the future with no query cost.
Algorithm Vertex Sampling by (Weighted) Degree
Require: Precision
Ensure: Reals pt such that (1 e)deg(x)) SP. S (1 e)deg(z;) for all 1 S S
1: for = 1 to = n do
2: p; <- KDEx (z.) (1 (1-=) k(z,
3: end for
4: Return {p.)21
Once we acquire {pr } we can perform a fast sampling procedure through the following algo- rithm, which state in slightly more general terms.
18

---

Published as a conference paper at ICLR 2023
Algorithm 3 Sample from Positive Array
Require: Input array A [01,'., ,0n] with a, 0 for all i. Access to queries A.j Eisis at
for 1 Sisj S.SJm n.
1: T A
2: while [T] 1 do
3:
m IT)(T)
4:
a - [m/2]] // Can be simulated using an Ais query
5:
b
b-- 2 (T[ + +1:m)
6:
if Unif
a/ then T+T[1, T[1: 1m/2]]
7:
T
8:
else
9:
T T 7 + 1: m]
end if
11:
end while
12: Return the single element in T
Combining Algorithms 2 and 3, we can perfectly sample from the degree distribution of the graph
K.
Algorithm 4 Faster Degree Sampling
Require: Reals Py such that (1 e))eg((.) S p, S (1+ (d)g((z)) for all < S 1: < index in [n], which is the output of running Algorithm 3 on the array {pi}
2: Return T;;
We now analyze the correctness and the runtimes of the algorithms proposed in Section C. First, we give guarantees on Algorithm 2.
Theorem C.2. Algorithm 2 returns {p: {p:)- such that (1 - =)dcg(Z4) S pt S (1 (x) all
Isisn.
Proof. The proof follows by the Definition of KDE query, Definition 1.1.
We now analyze Algorithm 3, which samples from arryy based on tree data structure given access to consecutive sum queries. The analysis of this process will also greatly facilitate the analysis of other algorithms from Section C.
Lemma C.3. Algorithm 3 samples index i E [n]] proportional in O(log n) time with O(logn) n) queries.
Proof. Consider the sampling diagram given Figure Algorithm 3 does the following: first queries the root node A1, and then its two children A1, Am+1,n where m = |n/2]. Note that
A A + Am+ It then picks the tree rooted at AI.m with probability Σ a. E. 9,
and oth-
erwise, picks the tree rooted at Am+ The procedure recursively continues by querying the root node, its children, and picking one of its children be the new root node with probability pro- portional to the child's weight given by an appropriate query access. This is done until we reach leaf node that corresponds an index [n].
We now prove correctness. Note that each node of the tree in Figure 1 corresponds to subset C [n]. We prove inductively that the probability of landing on the vertex equal to E dg- This is true for the root node of the tree since the algorithm begins at the root note. Now consider transitioning from some node S to one of its children S We know that we are at node S with probability Eies a:/ >.s 2; a,. Furthermore, we transition to S₁ with probability Eies, a:/ Ejes aj.
Therefore, the probability of being at S₁ is equal to
ΣiEs, ai Eiesa az Ejesaj Ejag
Eies, a; Σ, aj
19

---

Published as a conference paper at ICLR 2023
Since there is only path from the root node to any vertex of tree, this completes the induction.
The runtime and the number of queries taken follows from the fact that the sampling procedure descends tree with (logg) height.
Combining Algorithms 2 and 3 allows us to sample from the degree distribution of the graph K up low error in total variation (TV) distance.
Theorem C.A. Algorithm samples from the degree distribution of K up to TV error O(e) using a fixed overhead of KDE queries and runtime (log n)).
Proof. Since pj is with a 1 t factor of deg(z;) for all i, then {pr is O(E) close in total variation distance from the true degree distribution. Moreover, Algorithm perfectly samples from the array {p:)-= which proves the first part of the theorem.
For the second part, note that acquiring {p:]--1 requires n KDE queries. We can then construct the data structure for Algorithm 3 by computing all the partial prefix sums in O(n) time. Now the query access required by Algorithm can be computed in O(1) time through an appropriate subtraction of two prefix sums. Note that the previous steps need to be only done once and can be utilized for all
future runs of Algorithm 3. It follows from Lemma C.3 that Algorithm 4 takes (log n) time.
C.3 Weighted Edge Sampling and Weighted Neighbor Edge Sampling
We describe how perform weighted neighbor edge sampling.
Definition C.2 (Weighted Neighbor Edge Sampling). Given a vertex v, the goal of weighted neigh- (1=))M(s,
for all i E [m]] bor edge sampling is to output vertex such that Pr(u xn] Seen k(z..)))
Algorithm Sample Random Neighbor
Require: Input vertex x; X, precision & Ensure: E XX \\ {x;} such that the probability of selecting is proportional to k(x;, z)
1: e' = // log
2: while [T] > 1 do
3: m + [T]
4:
KDEr(I.m/)).e (I;)
5: b b (DEr(mm2+11 11mmm/11mmmmm)) e (z:)
6:
if E T[1 mm2 then
7:
a 4- a - (1 - ) 2i)
8:
end if
9:
if E Tim/2 + 1 m then
10:
b b- - (1 (')=E')(xx,....
11:
end if
12:
if Unif 1] a/(a (a b) then
13:
T T[1 : m/2]
else
15: T' T T(m/22 + m]
end if
17:
end while
18:
Return single element in T
We now prove the correctness Algorithm based on the ideas in Lemma C.3. Note that Algorithm 5 takes in input precision level which can be adjusted and impacts the accuracy of KDE queries. We will discuss the cost initializing KDE queries with various precisions in Section B.I.
Theorem C.S. Let æ; E X be input vertex. Consider the distribution D over X (x;), the neigh- bors of x, in the graph K. induced by the edge weights K. Algorithm 5 samples a neighbor from distribution that is within TV distance O(E) from D using Ollogn) n) KDE queries and Ollog n) time. In addition, we can perfectly sample from D using O(logn) n/T) additional kernel evaluations in expectation.
20

---

Published as a conference paper at ICLR 2023
Proof. The proof idea is similar to that of Lemma C.3. Given vertex 24, its adjacent edges have associated weights and our goal is to sample an edge proportion to these weights. However, unlike the degree case, performing edge sampling is not straightforward KDE query cares about the kernel value between two points, rather than the sum of kernel values that a KDE query provides. Nevertheless, in conjunction with KDE queries with over various subsets of X.
edge only
utilize the tree procedure outline in the proof of Lemma C.3
Imagine the same tree as in Figure I where each subset corresponds to a subset of neighbors of I; (note that x; cannot be its own neighbor and hence we subtract k(a;, 2;;) in line 7 or line 10). Algorithm descends down the tree using the same probabilistic procedure as in the proof of Lemma C.3: at every node, it picks one of the children to descend to with probability proportional to its weight. Here, the weight of a child node in the tree in Figure 1 the sum of the weights of the edges connecting to the corresponding neighbors of Ti+-
Now compare the telescoping product of probabilities that lands us in some leaf node aj to the ideal telescoping product if we knew the exact array of edge weights as in the proof of Lemma 6. Suppose the tree has height& f. At each node in our actual path descending down the tree, we take the next step according to the ideal descent (according to the ideal telescoping product), with the same probability, except for possibly an overestimate or underestimate by factor + or 1 - factor respectively.
Therefore, we land in the correct leaf node with the same probability as in the ideal telescoping product, except our probability
be off by a multiplicative (1 + e') factor. However, since e' e/ log and & log n, this factor is within s Thus, we sample from the correct distribution over the leaves of the trees in Figure up TV distance (s)) Now by doing O(1/T) steps of rejection sampling, we can actually get a prefect sample of the edge. This is because the denominator of the fraction for Pr(u zy} is at least ((nT) and most n so we can estimate the proportionality constant in the denominator by which is only at most ((1/T) multiplicative factor larger. Hence by standard guarantees of rejection sampling, we only need repeat the sampling procedure (1/T)) additional times.
Algorithm 6 Sample Random Edge by Weight
Ensure: Sample edge (I) with probability at least (1 = (x),
I: v + random vertex by using Algorithm
2: random Neighbor of using Algorithm
3: Return: Edge (u, w).
Theorem C.6 (Weighted Edge Sampling), Algorithm 6 returns random edge of K with probability proportional to at least (1 - 8) its weight using call to Algorithm 5.
Proof. Consider an edge (u, v). The vertex is sampled with probability at least (1 2e)deg(u)/ Ezex deg(z). Given this, the vertex v is then sampled with probability at least k(u,t)
(1-2e) Ezex,, R(u,)) (1 2s) k(a,)) Using the same analysis for sampling v and then u, we have that any edge (u, v) is sampled with probability at least 2= times k(u, u)/ EeeK w(e). Note that the same rejection sampling remark as in the proof of Theorem C.5 applies and we can perfectly sample an edge proportional to its weight with an addition ((1/r) rejection sampling steps.
deg(m))
C.4 Random Walk
Theorem C.7. Algorithm outputs vertex from a vertex within O(Tε) total variation distance from the true random walk distribution. Each step of the walk requires call to Algorithm 5.
Proof. The proof follows from the correctness of Algorithm 5 given in Theorem C.5. Lastly we again note that by performing additional O(1/r) rounds of rejection sampling steps (as outlined in the end of the proof of Theorem C.5), we can make sure that we are sampling from the tre random walk distribution at each step of the walk.
21

---

Published as a conference paper at ICLR 2023
Algorithm 7 Perform Random Walk
Require: Input vertex 2; X, length of walk T 2 1.
1: Start at vertex 2;
2: - 2i
3: for j to do
4: f- Sample random neighbor of v (Algorithm 5)
5: t 4- w
6: end for
7: Return 2
D Linear Algebra Applications
We now present a wide array of applications of the algorithmic building blocks constructed Sec- tion C. Altogether, these applications allow us to understand or approximate fundamental and prop- erties the kernel matrix and the graph K. In this section we present the linear algebra applications and the graph applications are given in Section E.
D.1 Spectral Sparsification
Algorithm Spectral Sparsification of the Kernel Graph
1: Let O log(n)/=²r*) be the number of edges that are to be sampled
2: Let ş denote the distribution returned by Algorithm for small enough constant
3: Initialize G =
4: for 14
5: Sample vertex from the distribution
6: Sample a neighbor v of u using Algorithm with constant e
7: Let Jue be the probability that Algorithm 5 samples given input
8: Similarly define and compute qva
9: wwe 1/(t(fa.qwv (o))))
10: Add the weighted edge ((u, v), Was) to the graph G'
end for
12: Compute an E/2 spectral sparsifier G" of graph G' using (Batson et al., 2013)
return G"
Given set X, |X| sparsifier for the weighted complete graph on X where weight of the edge (æ,, is given by kg);
and kernel A X X - R we describe how to construct spectral
Definition D.1 (Graph Laplacian). Given a weighted graph G (V, E, w). the Laplacian of G, denoted by LG = =D - A, where A is the adjucency matrix of G with A = wi vo((i,j)) and is diagonal matrix such that for all E [n), Da,; Ejzt Ai,j-
Theorem D.I (Spectral Sparsification of Kernel Density Graphs). Given a dataset X of points in Rdd and kernel k: X X TR+, let G (X, (x), w) be the weighted complete graph on X with the weights w(fx;, x, 22 k(#,, 2z)). Further. for all 2i, X, let k( 2x)) 2 for some T E (0, 1). Let LG the Laplacian matrix corresponding to the graph G. Then, for any E (0. 1), Algorithm 8 outputs graph G" with only m O(n log (s(773) edges, such that with probability
at least 9/10,
eG To < TT
The algorithm makes Ö(m/7²) KDE queries and requires Ö(md/7²) post-processing time.
Let G. be the weighted directed graph obtained by arbitrarily orienting the edges of the graph G and let H be edge-vertex incidence matrix defined as follows for each (I,, 25) in graph G.,,
let H Haz: Vk(1,, 255) and Ha,ry /k( T.). Note that H H IG- Our idea construct
22

---

Published as a conference paper at ICLR 2023
spectral sparsifier is to compute sampling- rewieghigg matrix S, i.e., matrix that has at most one nonzero entry in each row, that with probability 2 9/10, satisfies
(1 - LG = (1- = H S
SH (1 + E)H (1+++))
The edges sampled by S form the edges of the graph G'. We construct this matrix S by sampling rows of the matrix H from distribution close to the distribution that samples a row of H with probability proportional to its squared norm. We show that this gives a spectral sparsifier by showing that such distribution approximates the "leverage score sampling" distribution.
Definition D.2 (Leverage Scores), Let M be for all [n], To the i-th leverage of M is defined as follows:
d matrix and m, denote the i-th row of M. Then,
Tq = T,m,,MM T M) al)",
where X+ is the Moore-Penrose pseudoinverse for matrix X.
We introduce the following intermediate lemmas. We begin by recalling that sampling edges propor- tional to leverage scores (effective resistances on a graph) suffices to obtain spectral sparsification (Spielman & Srivastava, 2011; Woodruff, 2014).
Lemma D.2 (Leverage Score Sampling implies Sparsification). Given an n. X d matrix M and E (0, 1), for all i E [E], let т; be the i-th leverage score of M. Let p {p1 (PI.P2.... ..... be distribution over the rows of M such that pi Ts/ T// )ee(n] Tj. Further, for some & (0, (), let
p (21,02,.... ,..] be distribution such that P. 2 óp, and let l =
Let S E B be a random matrix where for all j E [1] the j-th row is independently chosen as (1/ [D.ee with probability pa- Then, with probability at least 99/ 99/100.
[dlog(d) e²
(1 ()- M 3 M S SM = (1++)M
Next, we show that the matrix H is well-conditioned, in fact the condition number independent of the dimension and only depends on the minimum kernel value between any two points in the dataset. This lets us use our edge sampling routines to compute an spectral sparsifier.
Lemma D.3 (Bounding Condition Number). Let H be the edge-vertex incidence matrix as defined and also has the property that all nonzero entries in the matrix have an absolute value of at most and at least VT. Let max (H) be the maximum singular value of H and m (H) ) be the minimum nonzero singular value of H. Then maxx (H)/ mii (H) < <4V2/71.5
Proof. We use the following standard upper bound on the spectral norm of an arbitrary matrix A to upper bound the spectral norm of the matrix H:
(( max Σ
|All2 =
max
|Aepl
For the matrix H, as each column has at most nonzero entries and each row has most 2 non- zero entries and from the assumption that all the entries have magnitude at most 1, we obtain that ||Hl||2 < V2n. To obtain lower bounds on Jmin (H), we appeal to Cheeger-type inequality for weighted graphs from (Friedland, 1992; Friedland & Nabben, 2002). First, we note that mmin (H) Vamin (HT)) H) v amin (La) where G is the kernel graph that we considering with each edge having a weight of at least r. Let 0 = X1 < A S ... < Xn be the eigenvalues of the positive semi-definite matrix LG. Now we have that
Omin = >a(Ta) 2 m)n(()/2)e((A)
where 6, = Ejee 21)) i.e., the weighted degree of vertex o; in graph and
E(U,Ü) U)
c(G)
min SS/
|E(U)
where [E(U)] denotes the sum of weighted degrees of vertices in U and |E U) denotes the total weight of edges with one end point in U and the other outside U. Using the fact that G is
23

---

Published as a conference paper at ICLR 2023
complete graph with each edge having a weight of at least and at most 1, we obtain E(U =| and [E()) & n/U/, which implies that E(G) 2 ninestrcv//jssss2
2 7/2. We
also similarly have that min; 6, 2 (n - 1)r, which overall implies that X2 >2( 2 nr /16 mmin (H) 2 //rr1./4. Thus, we obtain that max (H) mnin(H) N 4v///11.5
916 and that
We are now ready to complete the proof of our main theorem:
Proof of T'horre D.1. Let q (m, 92, 1.2..-..911 (3)) be a distribution over the rows of H such that for
all edges [i,j], 2 |[H |M -e- k(r,,ag))> for fixed universal constant C.
A(z..z))
Next, we show that this distribution is @(1/k²) approximation to the leverage score distribution for H. Let = UEV be the "thin" singular value decomposition of and therefore all the diagonal entries of Σ are nonzero. By definition T; ||U 112- We have
|1:1122 = I0 .V 112 = IU
where the equality follows from the fact that VT has orthonormal rows. Now, 10 > || 11g and || & |0+e 12 Therefore, for all i E (2), defining miin/ m
we have
11 11 Ti Σ, Σ, 1 Σ, 1/111//22nn x² |
|/4 1hal3//omas
1 |
Then, we invoke Lemma D.2 with Ø = 2(1/x²) and conclude that sampling ( =
log n. aس
rows
of H results in a sparse graph with corresponding Laplacian Lo such that with probability at least 9000
- E/2)/ 3 Lo = e/2)LG
Further, by Lemma D.3, we can conclude x² < 32/7* and thus sampling t = O suffices.
n log 2
edges
We do not use Algorithm 6 to sample random edges from the perfect distribution to implement spectral sparsification as we cannot compute the exact sampling probability of the edge that is sampled. So, first use Algorithm 4 with constant (say 1/2) to sample vertex u and Algo- rithm 5 with constant (say 1/2) to sample a neighbor v of u. Note that Algorithms 4 and Algo- rithms 5 can be modified to also return the probabilities Pu and qva with which the vertex i and the neighbor j of i are sampled. We can further query the algorithms to return P. and Puv. Now, Hu,w) Pulvu + Pu is the probability with which this sampling process samples the edge {u, v) k(zn.za) and we have that Pugeu + Pu 2 EEes k(z.2,)) and we use this distribution q to implement spec- tral sparsification as described above. As already seen (Theorem C.5), to compute vertex sampling distribution P, queries. Thus, we overall use On log² n=)(=>?)) constant approximate KDE queries to obtain spectral sparsifier.
KDE queries and for each neighbor sampling step, use (llog KDE
We can further compute another graph G" with only O(n/=2) edges by computing an e/2 spectral sparsifier for G' using the deterministic spectral sparsification algorithm of (Batson et al., 2013). Conditioning on the graph G' being an spectral sparsifier of G, the graph G" will then be an sparsifier for G. This procedure doesn't 't require any KDE queries and solely operates on the weighted graph G'.
Hardness for spectral sparsification. We observe that we the lower bound from Alman et. al. to establish hardness in terms of 7 from Parameterization 1.1. The lower bound we obtain is follows:
Theorem D.4 (Lower Bound for Spectral Sparsification under Parameterization 1.1). Let k be the Gaussian kernel and let be dataset such that ming.ye k(x, y) = T, for some 1 > T >> O. Then, any algorithm that with probability 9/ 10 outputs O(1 -approximaa spectral sparsifier for the
24

---

Published as a conference paper at ICLR 2023
kernel graph associated with X, with O(n²-6) edges, where 6 < 0.01 is fixed universal constant, requires Q n 210R(1/+)0.33) time, assuming the strong exponential time hypothesis.
First, we begin with the definition of a muliiplicatively-Lippe Definition D.3 MMulippiicaively-LLipcchii Kernels). A kernel k over a set X is (c, L)- multiplicatively Lipschitz if for any p (1/c,c), and for any 2, y X, kk(,y)) S k(px, S k(,
function:
We will require the following theorem showing hardness for sparsification when the kernel function is not serptiicaavvll--ipcchhe
Theorem D.5 (Theorem 8.3 (Alman et al., 2020)). Let k be a functio and X be dataset such that k is not (c, mimuliipiiicaiive--iissiiiii on X for some L > and 1 + 2 log (10 22 L.. Then, there no algorithm that returns sparsifier of the kernel graph associated with X with O edges, where 8 < 0.01 is a fixed universal constant, in less than O (n 220.40) time, assuming the strong exponential time hypothesis.
Proof of Theorem D.4. First, we show that for any i, < log((17)(c-1)) then the Gaussian kernel ie is not (c, I)-multiplicatively Lipschitz. Let || y12 and let f(=) e Observe, it suffices show that there exists such that f(cz) S c- f(2). Let be such that f(z) e² min,, k(2, 8) T, i... 2 = (1/r) Then,
f(clog(1/T)) -clg((11
and for < log(1/))(c-)) 1)
c- f
Then, applying Theorem D.5 with c = 1+ Vズ it suffices to conclude k is not (c. m)mlliippiiicaiivlllyy Lipschitz when L (1gr) which concludes the proof.
D.1.1 Solving Laplacian Systems Approximately
We describe how to approximately solve the Laplacian system LGa b using the spectral sparsifier Lg.. First, we note the following theorem that states the running time and approximation guarantees of fast Laplacian solvers.
Theorem D.6 (KKoutis al., 2011), (Spielman Teng, 2004)). There is an algorithm that takes an input graph Laplacian L of graph with weighted edges, vector b, and an error parameter and returns x such that with probability at least 99/ /000,
bll S alIL+bbll.
where ||0|| = V La. The algorithm runs in time Ö(m log(1/a)).
We have the following theorem that bounds the difference between solutions for the exact Laplacian system and the spectral sparsifier Laplacian.
Theorem D.7. Let LG be the Laplacian of connected graph G on vertices and let Lo be the Laplacian of an E-spectral sparsifier G' of graph G i.e.,
L1 3 Lo 3 (1+e)La.
for E < for small enough constant c. Then, for any vector b with 1 || 2/
L S
Proof. Note that for & < 1. the graph G' also has to be connected and therefore the only eigen vectors corresponding to eigen value 0 of the matrices La and To of the form 1 for
25

---

Published as a conference paper at ICLR 2023
and hence columns (and rows) of LG span all vectors orthogonal to 1. Therefore LoLiG LoGG I (1/n)11 Now,
|L b126 =b (L0 - -G Lb
=b LoL b' LOLL - bb L bo Lo 1 b Ligh
+b
b b% ligh 1-
where in the last inequality, we used L = and that for any vector v, T L S - Lg As the null spaces of both Lg and Lo are given by {al| la E IR), we also obtain that
(1 (1- 1 L 3 (1 =)
using which we further obtain that
|L bi 2 (1-2-2)
2e(1 +) 6J L
1-6 6J Lig < 4.
Thus, - bo $2VE
Therefore, if x is vector such that || Larblleo S all obtained using the fast Laplacian solver, then
| Lobllzo 1| L + L L
$ 2(11 - blb + +
- bb
2 -ー- |> l lig 4el| Tel
Here used the above theorem and the fact that La 3 (1/(1 Lor Now, ||- b S b1 and b = 6T Lo Lay Lig = T L b S (1 (1+ Lig s (1 (1+b)
which finally implies that
2(1
1 Libli <
1-8 a² 4e |
Thus, using spectral sparsifier G' with edges, in time Ö(m log(L/=)) vector such that || Lobllr I CCVE\ for large enough constant C.
obtain
D.2 Low-rank Approximation of the Kernel Matrix
We derive algorithms for low-rank approximations of the kernel matrix via KDE queries. We present algorithm for additive error approximation and compare to prior work for relative error approxi- mation.
We first recall the following two theorems. Let Ai,, denote the ith row of a matrix A.
Theorem D.8 ((Frieze et al., 2004)). Let A E Raxm be any matrix. Let S be sample of O(k/E)
rows according to a probability distribution (pI,....Pv) that satisfies pa 2 2(1) | / A for every 1 < i S Then, in time O Omnk/s poly(k, 1/E)), we can compute from S a matrix E Rkxm, that with probability at least 0.99 satisfies
||A - AA-AUT AU S JA As 1/-- + =|[AI|?
Theorem D.9 (CChen & Price, 2017), also see (Indyk et al., 2019)). There is randomized algorithm that given matrices A E gm and U R reads only O(k/E) columns of A, runs in time ((m)) + (kl 1/E), and returns V E Bank that with probability 0.99 satisfies
||A VU 3 S (1++)
min ||A XEReXt
X
26

---

Published as a conference paper at ICLR 2023
Therefore to compute the low rank approximation, we just need sample from the distribution required by Theorem D.8. We reduce this question to evaluating KDE queries as follows: If K is the kernel matrix, each row of K is the weight of the edges of the corresponding vertex. Therefore, each pt in the distribution (p1,...,,Pa Pa) the sum of edge weights squared for vertex I;. From vertex queries (Algorithm 4), we know that of edge weights. We can extend Algorithm 4 to sample from the sum of squared edge weights of each vertex as follows. Consider kernel such that there exists an absolute constant that satisfies k(x, y)² = k(cz, cy) for all x. y, Such exists for the most popular kernels such as the Laplacian, exponential, and Gaussian kernels for which we simply construct KDE queries for the dataset X" == cX. Then by sampling the degrees of the vertices associated with the kernel graph K" of X', we can sample from the distribution required by Theorem D.8 by invoking Algorithm 4 on the dataset X". In particular, using v. KDE queries for X', we can get row norm squared values for all rows our original kernel matrix K. We can then sample the rows according to Theorem D.8 and fully construct the rows that are sampled. Altogether, this takes n. KDE queries and O(nk/s) kernel function evaluations to construct a rank & approximation of K; see Algorithm 9.
get the degree of each vertex, which is the sum
2, and respectively. Thus give our dataset X,
Corollary D.10. Given a dataset of size n, there exists an algorithm that outputs a rank matrix B such that
||K BIB S S +
with probability 99/100, where K is kernel matrix associated with X based Laplacian, exponential, or Gaussian kernel, and K₁ is the optimal rank-k approximation of K. It uses n. KDE queries and O(nk/e poly(k, 1/s) + nkd/s) post-processing time.
We remark that for the application presented in this subsection, we can we can replace 1.1. Indeed, since we only estimate row sums, we only require that the value of KDE query is at least т, that is, the average value Іx Lmxx k(æ, y) 2 t for query y. Note that via Cauchy Schwartz, this automatically implies a lower bound for the average squared sum:
Mx,z)²
k(æ, y) k(2,)) Σ XE
2
X
X
ZE
Algorithm Additive-error Low-rank Approximation
1: Let c be the constant such that k(z, R k(cz, cy) for all inputs x,y
2: for = 1 to = n do
3: Compute the value p; E k(ex,, cry) using KDE queries for the dataset c.X
4: end for
5: Sample and construct Okk/E) rows of K according to probability proportional to {p.)
6: Compute U from the sample, using Theorem D.8
7: Compute V from the sample, using Theorem D.9
8: return Return V,
D.3 Approximating the Spectrum in EMD
In this subsection, we obtain sublinear time algorithm to approximate the spectrum of the nor- malized Laplacian associated with the graph whose adjacency matrix is given by the kernel matrix K.
The eigenvalues of the Laplacian capture fundamental combinatorial properties of the graph such as community structures at varying scales. See the works (Lee et al., 2012: Louis et al.. 2012; Kwok et al., 2013; Czumaj et al., 2015; Gluch et al., 2021), which show that the jth eigenvalue of the Laplacian informs us if the graph can be partitioned into j distinct clusters. However, computing a large number of eigenvalues of the Laplacian may not be computationally feasible. Thus, it is desirable to obtain succinct summary of all eigenvalues, i.e. the spectrum.
Additionally, models random graphs that aim to describe social or biological networks often times have closed form descriptions of the spectrum for graphs drawn from the model. Borrowing
27

---

Published as a conference paper at ICLR 2023
an example from (Cohen-Steie Sreen-Stineereee et al., 2018), "if the spectrum of random power-law graphs does not closely resemble the spectrum of the Twitter graph, it suggests that a random power-law graph might be poor model for the Twitter graph." Thus, another application of computing an approximation of the spectrum of eigenvalues is to test the applicability of generative graph models.
Our notion of approximation deals with the Earth mover (EMD) distance.
Definition D.4 (Farth Mover Distance). Given two multi-sets of n points in Rd,, denoted by and B, the earth-mover distance between A and B is defined as the minimum cost of a perfect matching between the two sets, i.e.
(M, = min Σ T:A-BB
-laa- (
(D.1)
aEA
where N ranges over all one-to-one mappings.
of Stohen-Steiner al., 2018). The algorithm first selects uniformly random vertices of a weighted graph A. It then performs random walk of specified length starting from the chosen vertex and then counts number of times the walk returns back the original vertex. Now Theorem C.7 allows to perform one step of random walk using Ollog n) KDE queries. Note that perform additional ((1/T) of rejection sampling in Algorithm to perfectly sample from the true neighbor distribution. Thus we immediately have the following guarantee:
We can now invoke the algorithm ApproxSpectra
Theorem D.11 (Corollary of Theorem in (Cohen-Steiner et al., 2018) and Theorem C.7). Given X kernel matrix K and accuracy parameter E (0, 1), let G be the corresponding weighted graph, and let LG = I D-*KD"" be the normalized Laplacian, where Di, = -E; Kig. Let X1 2 >2... 2 An be the eigenvalues of LG and let be the resulting vector. Then, there exists an algorithm that uses Ô (exp (1/=²)/ /T) KDE queries and exp (1/=²) d//T post-processing time and outputs vector X such that with probability 99/100,
EMD (^.X) SE.
We remark that the bound of exp (1/6²) independent of n, which is the size of the dataset.
D.4 First Eigenvalue and Eigenvector Approximation
Our goal is to approximate the top eigenvalue of the kernel matrix and find vector witnessing this approximation. Our overall algorithm can be split into two steps: first sample a random principal submatrix of the kernel matrix. Under the condition that each row of the satisfies that it's sum is at least nT, we can easily show that it must have a large first eigenvalue and thus prior works on sampling bounds automatically imply the first eigenvalue of the sampled matrix approximates that of K. The next step is to use 'noisy' power method of (Backurs et al., 2021) 1) on the sampled submatrix. We note that this step employs KDE structtr initialized only on the sampled indices of K. The algorithm and details follow.
n kernel matrix K
Algorithm 10 First Eigenvalue and Eigenvector Approximation
Require: Input dataset X Rdd of size |X| n,, precision
1: t O(1 (11((2?7)
2: S random subset of [n] of size
3: Ks < principal submatrix of K on indices in S {Just for notation; we do not initializing K or
Ks}
4: K + (n/s). Ks
5: Return the eigenvalue and eigenvector found by Algorithm of (Backurs et al., 2021) (Kernel
Noisy Power Method) on Kg.
We remark that the eigenvector returned by Algorithm 10 will be a sparse vector supported only on the coordinates in S.
28

---

Published as a conference paper at ICLR 2023
We first state the necessary auxiliary statements needed to prove the guarantees of Algorithm 10. Lemma D.12. If each row of K satisfies that its sum is at least nr for parameter (0, 1), then the largest eigenvalue of K, denoted as 21, satisfies X 2 WT.
Proof. This follows from looking at the quadratic form 1" K1 where 1 is the vector with all entries equal to 1:
1 K1
n²T
A1
1
nT.
We now state the guarantees of Algorithm in (Backurs et al., 2021).
Theorem D.13 (Backurs et al., 2021)). Suppose the kernel function for a m X m kernel matrix K has KDE data structure with query time d/(=?r") (see Table 1). Then Algorithm of (Backurs et al., 2021) returns À such that X 2 (1-e)Xi(K) in time O dm +2 (og(m/e)
Finally, need the following result eigenvalues of sampled PSD matrices, proven in (Bhat- tacharjee et al., 2021).
Lemma D.14 (Bhattacharjje et al., 2021)). Ler A E mn be PSD with || 1. Let S C [n] be a random subset of size and let Asxs the submatrix restricted columns and rows S and scaled by n/s Then, for all i E [[S|]. A (Asxs) (x(A)
We are now ready to prove the guarantees of Algorithm 10.
Theorem D.15. Given n x n kernel matrix K admitting KDE data-structure with query time dP( Algorithm 10 returns X such that À 2 (1 - (K) total time
min (o (o(dlog(d) g.574
24t )
d
2+p
log
Remark D.1. Two remarks are in order. First we recall that the runtime of (Backurs et al., 2021) has I+ factor while bound has dependence
and is thus truly sublinear runtime.
Second, if we skip the Kernel Noisy Power method step and directly initialize and calculate the top eigenvalue of Ks (using the standard gap independent power method of (Musco Musco, 2015)), we would get a runtime of (dd((s1.5/1 which has polynomially better dependence but worse r dependence than the guarantees of Algorithm 10.
Proof of Theorem D.15. We first prove the approximation guarantee. By our setting of and using Lemma D.14, we see that the additive error in approximating the first eigenvalue of K by that of K most
VE S
and thus (K) (1-))\( Then by the guarantees of Theorem D.13, it follows that we find a multiplicative approximation to A1 \T(K) and thus 1 O(s) multiplicative approximation to that of X₁ (K).
We now prove the runtime bound. It easily follows from plugging in m 0(1/(E?73) in Theorem D.13.
E Graph Applications
In this section, we present our graph applications, including local clustering, spectral clustering, arboricity estimation, and estimating the total weight of triangles.
29

---

Published as a conference paper at ICLR 2023
E.1 Local Clustering
Algorithm 11 Local k- k-Clustering
Require: Vertices u, w, random walk length &
For given u. let pl be the endpoint distribution random walk of length starting at v
2: if f2 distribution tester outputs lpp Pilll S 1/((m) then
3: return u, w are in the same cluster
4: end if
5: return u, w are in different clusters
We give local clustering algorithm on graphs. The advantage of this method is that it is local allows us to cluster one vertex at a time. This is especially useful in the setting of local clustering where one might not wish to classify all vertices once or only a small subset of vertices are of interest.
We now present a definition for a clusterable graph that has been an extremely popular model def- inition in the property testing and sublinear algorithms community (see (Kale & Seshadhri, 2008; Czumaj Sohler, 2010; Goldreich & Ron, 2011; Czumaj et al.. 2015; Chiplunkar et al., 2018; Dey et al., 2019; Gluch et al., 2021) and the references within).
First, we need to define the notion of conductance.
Definition E.1 (Conductance). Lei G = (V, F,, w) be a weighted graph. The conductance of a set S C V is defined as
w(S.5 5
Sc( min(w((), (S))
where (S) denotes the sum of edge weights crossing the cut (S, S and w(S) denotes the sum of (weighted) degrees of vertices in S. The conductance of the graph G is then the minimum of (a(S) over all sets S:
((G) = min (S))
Definition F.2 (Inner/Outer Conductance). For a subset U C V. we define (GUU) to be the conductance of the induced graph U. (G) is also referred as the inner conductance U. Conversely, @U( is refereed to as the outer conductance Uf
Definition E.3 (k-clusterable Graph). A graph G is (k, din. dout) clusttebbeee if the following holds: There exists partition of the vertex set into h S k parts V Uisigh V4 such that (GG[Vi) din and dG (V,) < Pout-
Definition E.3 captures the intuition that one can partition the graph into < k pieces where each piece has strong cluster structure (captured by Φin) and distinct pieces are separated by sparse cuts (captured by oowx). Note that we are interested in the regime where dout is smaller than Φin. We will also assume that each [V] 2 n/ poly(k) where we allow for an arbitrary polynomial dependence on k. This means that each cluster size is not too small.
Since we are interested in clustering, through this section, we will assume our kernel graph K is k-clusterable according to Definition E.3 but we do not know what the partitions are.
The main algorithmic result of this section is that given a k-clusterable kernel graph and two vertices and w that are in parts V1 and V22 respectively of the graph (as defined in Definition E.3), we can efficiently test Vi = V/2 or V * V2. That is, we can efficiently test if and w belong to the same distinct clusters. The underlying idea behind the algorithm is that if and w belong to the same cluster, then random walks starting from these vertices will rapidly mix inside the corresponding cluster. Therefore, random walks in distinct clusters will be substantially different and can be de- tected using distribution testing. Our algorithm is given in Algorithm 11. The flavor of the algorithm presented quite standard in property testing literature, see (Czumaj et al., 2015) and (Peng, 2020).
The & distribution tester we need a standard result in distribution testing with the following guarantees.
30

---

Published as a conference paper at ICLR 2023
Theorem E.1 (Theorem 1.2 in (Chan et al., 2014)). Let 8, 0 and let p. be discrete distri- butions over set of size with b 2 max llq12). Letr coblog(L/))/ for appropriate constant c. There exists lz distribution tester that takes as input r samples from each distribution p, q and accepts the distributions [|p gll2 S and rejects the distributions if 1/ qll 2 4 with probability at least 8. The rumning time of the tester is linear in its sample size.
We now prove the correctness of Algorithm 11. We note that many arguments from prior works re-derived in the proof below, rather than stating them in black box manner, for completeness since our setting is of weighted graphs and the usual setting in literature is unweighted or regular graphs. We first need the following lemmas. Recall that the random walk matrix of an arbitrary weighted graph given by M AD where A is the adjacency matrix and D is the diagonal degree matrix. The normalized Laplacian matrix is defined as T.
D- AD-1/2
Our first result that vertices the same well connected cluster of G have a quantitative relationship captured by the eigenvectors of L. This is in similar spirit Lemma 5.3 of (Czumaj et al.. 2015) but we must show it holds for weighted graphs arising from kernel matrices whereas (Czumaj et al., 2015) is interested in bounded degree unweighted graphs.
Lemma E.2. Let v be the ith eigenvector of the normalized Laplacian of the kernel graph K and Z Din. Then for any s S h, the following holds:
let C be any subset such that & {[([)])
v:(u)
Ps(T))
Dout 8 C|
v /w(uu
Vw(v)
aEVEC
Proof. By Lemma 5.2 in (Czumaj et al., 2015) and Theorem 1.2 in (Lee et al., 2012), we have that & Laa/hi S An+1 and X S 2d for any 1 S h. Now by the variational principle for eigenvalues (Chung & Graham, 1997), we have
vu(u) Vro(u)
vi(u)
Σ (u,v)
w(u, v) S 2dout.
v Vw(w)
Now let H K[C]. From (Chung & Graham, 1997) and our assumptions on C, we have that
2
2. E(u,oeEE
ve(w)
v(le) mn(w) du (u)dn(v)
w(u, v)
Vu(u)
Φ 2 A2 2
vol (VH)
2
ΣxveVH
vu(u) Vi(u)
vi(D) Vi((v)
where vol (Vn) denotes the sum of the degrees of vertices H and d (.) denotes the degree in H. Note the last step is due to Cheeger's inequality. Combining the preceding result with our earlier derivation, we have
2
2
Vu(u)
vi(v)
Vi(u)
vsle) Vw(c)
Σ (n,
Σ
2do
(w)))
ww(v)
Vw(u)
(u,o)EEK
This implies that
2 du dH(u))H(v)
vi(u)
ve(t) V /w(lv)
vi(u)
vi(u) /w(o)
10
Σ
v w(()
Vw(u)
(W,o)VV
Dout vol (VH)
where we have used the fact that all edge weights in K are at least 7. Using the fact that vol (VH) s [C|n, it follows that
Σ
vi(u)
vi(u)
Douth
Vuluu)
/w(v)
2
CC
(o
desired.
31

---

Published as a conference paper at ICLR 2023
The second result states that vertices in the same well-connected cluster have similar random walk distributions. This is again the analogue of Lemma 4.2 in (Czumaj et al., 2015) but we must show holds for weighted graphs.
Lemma E.3. Let 0 B 1/2. If graph K is (k, Psn. Dour) clusterable, and C C V is any subset such that |C| Z n/ poly(k) and ((K[C]) 2 Oin. There exists a constant c = c(3) >> 0 and
è' c'= c(3
(5. k) such for any 2 clog n Pout S <do there exists subset C C C' satisfying CCl(C) (1 ())ol(C)) such that for any u, U E C, the following holds:
1 8n
ll 111 <
Proof. Let Vi...., Vn denote the eigenvectors of L with eigenvalues ,., ' in non-decreasing order. We know that the eigenvalues of M are given by X; with corresponding eigenvalues y; y,, The vector pt, is the vector 1, with one value in the wth coordinate applied to M Write
lu = Σ = Σ E..D-122
Taking the innerproduct of lu with D-1/2y, tells us that a, = v h00l/V/il /w()). Thus,
pu - p = Σ vz(u)
Σ vi(u)
(1
>o);y;
(1 Ar)'y:
i=1 v Vw(u)
i=1 Vw(v)
vi(u)
v.(()
= D1/2 Σ vi Vw(u)
(1 (1-Aa)
Vw(r)
This means that
vi(u)
ve(o)
lp pt,122 - |
Σv vi
(1 2,)²
Vw(u)
Vw(())
Since the ve's are orthogonal, we know that
n
ve(v:) Vw(u) Vw(()) )
vi(u)
vu(u)
ve(u) )² /w(w) (1 (1- 2s)
Σ V:
(1 (1- A:)²
Vu(u)
Now the rest of the proof follows from Lemma 4.2 in (Czumaj et al., 2015). In particular, it tells us that the above summation, each of the terms for S i S can be bounded by S 3 whereas the rest of the sum can be bounded by 1/ poly(n) n) for sufficiently large poly(n) by adjusting the constant in front of l. Our choice for |C| 2 n/ poly(k) imply that the overall sum is bounded
Caut
done (kly(k) Since | n, and Dont c' d we have that llll. pt 1 < 1/(8n). as 3n
by S desired.
Our next goal is to show that vertices from different well-connected partitions have very different random walk endpoint distributions. The argument we borrow is from (Czumaj & Sohler, 2010).
Lemma E.4. Let G be a (k, Φin, Dous) ) cclusterable graph with parts = Ui Vz. There exists constant c >> O such that if t- L-Qout S ce, then there exists subset Vi S V1 satisfying vol(V))) 2 (1
e))ol(Vi) such that a t-step random walk from any vertex in Ví does not leave V1 with probability 1-E.
Proof. We first bound the probability that the random walks always stay in their respective clusters. Consider fixed partition Vi; the same arguments apply for any partition. Let G" be the graph with the same vertex set as K but with only the following edges: edges among vertices in Vi and edges from vertices in V1 V \ Vi. Consider random walk on G' of length t with the initial vertex u' chosen from the stationary distribution of G', i... the distribution that chooses each vertex in G' with probability proportional to its weight. Let Y, denote the indicator random variable for the event that
32

---

Published as a conference paper at ICLR 2023
the ith vertex of the random walk is in V Vi. Since are simulating the stationary distribution, we have that
Pr[y, = 1| = w(V1, V Vi)
w(G")
where w is the weight of edges in the original graph K. By linearity of expectations, the number of vertices that land in V 1/2 is
t
(VV, VAIV Vi)
E
Σ 1=1
Y2
= (t
S toout
wo w=(')
due to our requirement of dout Therefore by Markov's inequality, the probability that any vertex in V V is ever visited is $ Doott-
We now move random walk back the original graph K. The preceding calculation implies that the probability that an step random walk in K starting at a vertex chosen at random from V1 according the stationary distribution will remain in V1 with probability least Doutt. t. If Pout S e/t, then we know that the random walk stays in V1 with probability at least 1 - O(e) so there must be set of vertices V' S V1 of at least O(s) fraction of the total volume of V such that random walk starting from vertex in V remains in Vi with probability at least 1 Os)). O
We can now prove the correctness of Algorithm 11.
Theorem E.S. Let K be (k, oin dous) ) cluuttrrbble kernel graph with parts V = Uisign V₁. Let U, W be one of (not necessarily distinct) partitions V4. Let u, w be randomly chosen vertices in partitions U and W with probability proportional to their degrees. There exists e d(Ek) k) such that if Dout S co log n, then with probability at least 1 E, if = W then Algorithm 11 returns that u and w are in the same cluster and if + W, Algorithm 11 returns that and w are in different clusters. The algorithm requires ol/nk//(es) log(1/e)) random walks of length & 2 clog n /n
Proof. We first consider the case that U * W. From Lemma E.4, we know that there are 'non- escaping" subsets U' and W' of U and W respectively such that vertices u, from U' and W respectively don't 't leave U and W with probability 1 - e. Conditioning on u and w being in those subsets, have that with probability 1 O(E), ph, and pt, will be disjointly supported and thus, llp' 11 = Il + 11p:112 2 2/n.
Now if U W, know from Lemma E.3 that lr - p311 1/(88) if condition on and coming from the large volume subset of U.
Finally, we need one last ingredient. Lemma 4.3 in (Czumaj et al., 2015) readily implies that there exists a V C V satisfying vol(V') > (1 (1-e) eV)o(V such that [p < 2k/((2 Now we can set & = 1/((n) and 2k/(E72n) in Theorem E. which tells us that = O0 (/n/////eeit (1/))) samples of the distributions pt, and pt, suffice to distinguish the cases ||pt - 1 2 2/n or l|pt
pt 1/(8n), i.e., r samples allow us to determine if U W U + W, conditioned on 1 - OEe) probability event.
straightforward to translate the requirements of Theorem E.5 in terms of the number of KDE queries required. Note that since we only take random walks of length O(log n/ ). we can just reduce the total variation distance from the distribution we sample our walks from and the true random walk distribution appropriately in Theorem C.7. Alternatively, we can perform rejection sampling as stated in the proof of Theorem C.5.
Corollary E.6. Algorithmm 1I and Theorem E.5 require /uk//e- (nn) KDE queries (via calls to Algorithm 7, which performs random walks) well the same bound for post- processing time.
E.2 Spectral Clustering
We present applications to spectral clustering. In data science, spectral clustering is often the follow- ing clustering procedure: (a) compute & eigenvalues of the Laplacian matrix in order, (b) perform k-means clustering the Laplacian eigenvector embeddings of the vertices.
33

---

Published as a conference paper at ICLR 2023
The theory behind spectral clustering relies on the fact that the Lapalacian eigenvectors effective in representing the cluster structure of the underlying graph. We refer the reader to (Von Luxburg, 2007) and references within for more information. For application spectral clustering, show that a spectral sparsifier, for example one computed from the prior sections, also preserves the cluster structure of the graph.
Next we define a model of a "weakly clusterable" graph. Intuitively our model says that a graph k-weakly clusterable its vertex set can be partitioned into k "well-connected" pieces separated by sparse cuts in between. Furthermore, this definition captures the notion of a well-defined cluster structure without which performing spectral clustering is meaningless. Note that this notion less stringent that the definitions of clusterable graphs commonly used in the property testing literature which additionally require each piece to be well-connected internally, see Definition E.3.
Definition E.4 (Weakly clusterable Graph). A graph (k, dout) clusterbbll if the following holds: There exists partition of the vertex set into h S k parts V = Uisigh V; such that (G(V.) & Pout-
We now prove the following result that says spectral sparsification preserves cluster structure ac- cording to Definition E.4. We first remark that the spectral sparsifier obtained in the previous section is a cut sparsifier as well. Recall that cut sparsifier is subgraph that preserves the values across all cuts up to relative error + e. The implication follows immediately by noting that cuts are induced by quadratic forms on the Laplacian matrix using {
13 vectors.
Theorem E.7. Let G be (k, (ou clusserrl and let G' be cut sparsifier for G. Then G' is (k, (1 (e==)))))) clusterblle
Proof. Let V; be of the h S k vertex partitions of G. Consider the conductance of V, defined Definition E. 1. The numerator represents the value of a cut separating V; and each term in the denominator is the sum of the degrees of single vertices. Both values appropriate cuts in the graph. Since G' is a cut sparsifier, this implies that both the numerator and denominator are preserved up factor and thus, the entire ratio is also preserved up O(e) factor.
Theorem E.7 implies that the cluster structure of the sparsified graph G' is approximately identical to that of G. Thus, we can be confident that the spectral clustering procedure described at the beginning of the section would perform equally as well on G' as it would have on G. Indeed, we verify this empirically in Section 2.
Furthermore, spectral clustering requires to compute the first k eigenvectors of the Laplacian matrix. Since our sparsifier has few edges, we can use Theorem (Musco Musco, 2015), which says (a variant of) the power method can quickly find good approximations Laplacian eigenvectors the matrix is sparse.
Theorem E.8 (Corollary Theorem (Musco & Musco, 2015) and Theorem D.1), Let L be the Laplacian matrix of the sparsifier computed in Theorem D. Let u1, ug be the first k eigenvec- tors of L. Using Theorem 1 of (Musco & Musco, 2015), we can find l vectors U1,"." Ug in time kn logn such that with probability 99/100 19/100. 7222.5
|/ Lu, of Luil S EX2211
for all i E [k).
34

---

Published as a conference paper at ICLR 2023
E.3 Arboricity Estimation
Algorithm 12 Arboricity Estimation
1: A EE w(e) m((e))
2: m < ni log n
3: for to
do
4: Sample an edge with probability
where the E [Wee 2wcel
Σ
5: Add e to G' with weight
6: end for
7: return maxugv d(G)))
We now apply our algorithmic building blocks to the task of arboricity estimation. Consider a weighted graph G = (V, E, Let Gy be induced subgraph of G on the subset of nodes U. The density of Gy is defined as
d(Gu) w(E(Gu))
[U]
where wi( w(E(Gu)) is the sum of the edge weights of Gu. The arboricity of G' is defined as
max d(Gu).
The arboricity measures the density of the densest subgraph in a graph. Intuitively, it informs if there strong cluster structure among some subset of the vertices of Therefore, an important primitive in the analysis of massive graphs with applications ranging from community detection social networks, spam link identification, and many more; see (Lee al., 2010) for survey of applications and algorithmic results.
Although polynomial time algorithms exist, we are interested in efficiently approximating the value of a using the building blocks of Section C. Inspired by the unweighted version of the arboricity estimation algorithm from (McGregor et al., 2015), we first prove the following result.
Theorem E.9. Let U' arg maxy d(G')) and let G' be the output of Algorithm 12. Then with probability at least 1/ poly(n).
e)a I d(Grr I (1+
Algorithm 12 uses m Ö(n logglog//(s (s(+)) KDE queries and ((m)) post-processing time.
Proof. Let U be arbitrary set of k nodes, let W Eeer w(c), and let Wo -W- :d((G)) Since G has weight W, then the arboricity a satisfies a > n
W so that
log n a
Let X; be the random variable denoting the contribution of the i-th sample to weight of the edges in G and observe that [[X]] m d(Gu) so that E[X] = W d(Gu), for X = mm Similarly, we have
E
Σ
2
Σ
Ww
W Wu
p2m²
m²
m2
e= M..
Since
CnA logn for absolute constant C >>> O. then
k
E|X?|
<
C'm log
and
kE²a²
> E
C' log" n
35

---

Published as a conference paper at ICLR 2023
We also have X; B[X,] I
ae
Thus by Bernstein's inequality for sufficiently large C, a
C' loga
Pr (r(d(G')) z
In
10k
10.
for d(Gtr) S 60
and
Pr [d((') d(Gu)l
EQ
2n -1kk
10
for d(Gu) 60
Since there are (*) s n² subsets of with size then by a union bound, we have that with probability least 1 3n 9k both
d(Gr)) $ 10
for all subsets U with d(Gu)
S 60
and
- ())r)
S d(G))) s (1+ (1+e))(Ge))
for all subsets U with d(Gr)
> 60
Hence for a set U such that d(U') = a, we have d(Gu-) 2 (1-2)a a so that d(G 2 d(Gr.) (1 E)a, where U arg uax d (G'y) Thus with high probability, we have that
8)a S d(Grr) S + e)a,
as desired.
To estimate the arboricity of the input graph G. it then suffices Theorem E.9 to compute the arboricity of the subsampled graph G' output by Algorithm 12. This can be efficiently achieved by running offline algorithm such as (Charikar, 2000), which requires solving a linear program on m variables, where is the number of edges of the input graph. Thus subsampling procedure serves preprocessing step that ultimately significantly improves the overall runtime.
E.4 Computing the Total Weight of Triangles
Algorithm 13 Weighted Triangle Counting
1: Let C E be a random set of O
my /wgwwal
max
edges
e²
m 2
2: + -0
Wiiin
3: For o E V, g(v) := E(u.)eg w(u, v)
4: For vu E E, E,g(u, == min(w(u), w()))
5: for 1 to
do
6: Sample E It with probability E.ER s(e)
g(e)
7: Let e (x. with x <
8: Sample neighbor of with probability
w(w)
s(e)
9:
if (x, y, 2) is a triangle assigned to u then
10:
Xs +-1
11: else
12:
X: + 0
13: end if
14:
end for
15:
return
[R Σ: Xi
We apply the tools developed in prior section counting the number of weighted triangles of kernel graph. Counting triangles is a fundamental graph algorithm task that has been explored in numerous models and settings, including streaming algorithms, fine-grained complexity, distributed
36

---

for Published as a conference paper at ICLR 2023
shared- and MapReduce to name few (Seshadhri al., 2013; Atserias al., 2008; Bera & Chakrabarti, 2017: Kolountzakis et al., 2010; Chen et al., 2022). Applications include discovering motifs in protein interaction networks (Milo et al., 2002), understanding social networks (Foucault Welles et al., 2010), and evaluating large graph models (Leskovec et al., 2008); see the survey (AI Hasan & Dave, 2018) for further information.
We define the weight of a triangle as the product of its edges. This definition natural since generalizes the case were the edges have integer lengths. In this case, edge can be thought of as multiple parallel edges. The number of triangles on any set of three vertices must account for all the parallel edge combinations. The product definition just extends this the case arbitrary real non non-negative weights. This definition has also been used in definitions of clustering-coeffiee weighted graphs (Kalna & Higham, 2006; Li al., 2007; Antoniou Tsompa, 2008).
Note that there is an alternate definition for the weight of a triangle in weighted graphs, which is just the sum of edge weights. In the case of kernel graphs, this is not an interesting definition since we can approximately compute the sum of all degrees using n. KDE queries and divide by 3 to get an accurate approximation.
Definition E.5. Let G = (V, (,Ew w) with w : R20 be weighted graph. Given a triangle (z, C E, we define its weight as
W(z,a.s) wwx w(y, 2) w(x, 2);
where we abuse notation by defining w(z, v) = w w(2 y)).
For this definition, we present the following modified algorithm from (Eden et al., 2017), which considers the problem in unweighted graphs in a different model of computing. See Remark B.S for comparison.
Theorem E.10. There exists an algorithu that makes Ō (
KDE queries and the same bound for post-procsssing time and with probability at least 3. outputs (1 + appproxiimtii to the total weight wT of the triangles in the kernel graph.
n, | = m. Eesg w(e) wa, T be the number of triangle in G, and wr be the sum of the weighted triangles in G, where the weight of a triangle (r, 2) C E, the product of the weights of its edges
Proof. Given graph G E), let [V|
W(z,0.2) = w(x y) wly, 2) w(x, z).
For vertex E V, let w(e) = set V by < e if and only if either wr(u) S wo()) or w(u) w((v) and appears in the dictionary ordering of the vertices. For each edge = (u, v). we assign to e all triangle (u, U, w) such that W. Let We denote the weight of the triangles assigned to e.
=(t)= -E (u..), w(e), that we have ordering on the vertex
Suppose, by way of contradiction, there exists E with We > VwG Since each triangle can contribute at most wi waxx weight to We, then more than / 3/2 triangles must be assigned
3/2
to e. Then there must be more than Vwgw-8 33/2 vertices with weight at least /ugwmaxx which contradicts the fact that the graph G has weight wG. Thus, have that W. Vwgwwas for all e E E. Moreover, we have that D.er We = wr so that
3/2
3/2
Wr m E [ee
and
WT W
E.Ep [
<
m By Chebyshev's incquality, it suffices to sample set R with
3/2 m m wwi Wmax
R = 0
euT
37

---

Published as a conference paper at ICLR 2023
edges uniformly at random, that
Pr
Σ LEER
W. E (1
WT 2 0.99.
For each vertex V E V, let g(v) = L(aw)EE w(u, v) and for each edge = (u, v), let g(e) =
min(w(u), w(b)) We write g Leen O(e). Now for each i [|R]], we have
We
WR g
B
Σ cE
B 1
(
Hence by Bernstein bounds, there exists constant C > 0 such that it suffices repeat the procedure
with probability least 2/3. We have that WR 2 wr [R] and Big( R)] Σ Vwgwwax 3/2 |R so that
R
times to get (1 = apprrxxiimtioo of s(R)
Wa
2m
g R
2m WGWmax
WR
WT
F
Omitted Experimental Results
MNIST
MNIST
220
KDE
350
300
SVD
NDE Approximatton
250
Frabohii
200
140
1.00
120
Rank
Normed Squared
(a) Rank versus Error for Low-rank Approximation
(b) Real vs Approximate Row Norm Squared Values
Glove
Glove
137.5
KDE Acorotimation
800
135.0
132.5
300
130.0
200
127.5
100
122.5
10
50
100
200 250 300 350 400 True Normed Sauared
150
Rank
(c) Rank versus Error for Low-rank Approximation
(d) Real Approximate Row Norm Sqaured Values
Figure 2: Figures for low rank approximation experiments.
Parameter settings. For low-rank approximation datasets, choose the bandwidth value cording to the choice made in prior experiments, in particular the values given in (Backurs et al., 2019). There, is chosen according the popular median distance rule; see their experimental sec- tion for further information. For our clustering experiments, we pick the value a, which results in spectral clustering (running the full kernel matrix) successfully clustering the input.
38

---

Published as a conference paper at ICLR 2023
Datasets for spectral sparsification and clustering. For spectral sparsification and clustering, we use construct two synthetic datasets, which are challenging for other clustering method such as k-means clustering². The first dataset denoted 'Nested' consists of 5, 000 points, equally split among the origin and a circle of radius 1. The two natural clusters are the points at the origin and the points on the circle. Since one cluster is contained in the convex hull of the other, a method like k-means clustering will not be able to separate the two clusters, but it is known that spectral clustering can. Our second dataset, labeled 'Rings', is an even more challenging clustering dataset. We consider two tori in three dimensions that pass through the interior hole of each other, i.e., they interlock. The 'small' radius of each tori is 5 while the 'large' radius 100. Our dataset consists of 2500 points uniformly distributed on the two tori; see Figure 3b. Note that our focus is not to compare the efficacy of various clustering methods, which done in other prior works (e.g., see footnote 2). Rather, we show that spectral clustering itself can be optimized in terms of runtime, space usage, and the number of kernel evaluations performed via our algorithms.
Evaluation metries for spectral clustering and sparsification. For spectral sparsification and clustering, we compare the accuracy of our method to the clustering solution when run on the full initialized kernel matrix.
Note that prior works such as (Backurs et al., 2019; 2021) have used use the number of kernel evaluations performed (i.e., how many entries of K are computed) as measure of computational cost. While this is a software and architecture independent basis of comparison, which is unaffected by access to specialized libraries or hardware (e.g., SIMD, GPU), it is of interest to go beyond this measure. Indeed, we use this measure as well as other important metrics as space usage and runtime as points of comparison.
Spectral sparsification and clustering results. Our algorithm consists of running the spectral spar- sification algorithm of Theorem D. (Algorithm 8) and computing the first two eigenvectors of the normalized Laplacian of the resulting sparse graph. We then run k-means clustering on the com- puted Laplacian embedding for k 2. As noted above, we use two datasets that pose challenges for traditional clustering methods such as k-means clustering. The Nested dataset is shown in Figure 3a. We sampled 3- 10 many edges, which is 2.5% of total edges. Figure 4a shows the Laplacian embedding of the sampled graph based on the first two eigenvectors. The colors of the red and blue points correspond to their cluster in Figure 3a as identified by running k-means clustering on the Laplacian embedding. The orange crosses are the points that the spectral clustering method failed to correctly classify. These are only 23 points, which represent a 0.5% of total points. Furthermore, Figure 4a shows that the Laplacian embedding of the sampled graph is able to embed the two clus- ters into distinct and disjoint regions. Note that the total space savings of the sampled graph over storing the entire graph is 4lx. In terms of the time taken, the iterative SVD method used to calculate the Laplacian eigenvectors took 0.18 seconds on the sparse graph whereas the same method took 0.81 seconds on the entire graph. This is a 4.5x factor reduction.
We recorded qualitatively similar results for the rings dataset. Figure 3b shows a plot of the dataset. We sampled 10 many edges for the approximation, which represents a 3.3% of total edges for form the sparse graph. The Laplacian embedding of the sparse graph is shown in Figure 4b. In this case, the embedding constructed from the sparse graph was able to separate the two rings into disjoint regions perfectly. The time taken for computing the Laplacian eigenvectors for the sparse graph was 0.08 seconds whereas it took 0.27 seconds for the full dense matrix.
G Auxiliary Inequalities
Theorem G.1 (Bernstei'' 's inequality). Let Xi...., X, be independent random variables such that E[X?] < oc and X, 2 0 for all E [n]. Let X = E; X, and 7: > 0. Then
Pr [r < E 7 $ exp
For see https::/acciittteea htt..org//talll/atoeeee..e. es/cluster/pl cl cluster iuster_compariss son, html.
39

---

Published as a conference paper at ICLR 2023
Nested Dataset
Rings Dataset
10.0
7.5
5.0
-100
-10.0
10
-200
150
-10.0
-7.
-5.0
-25
0.0
2.5
5.0
7.5
10.0
(a)
(b)
Figure 3: (a) Nested Dataset, (b) Rings Dataset
Rings Dataset
Orighn
Lapoccan pigenvociar
0.01
-0.03
0..00
3.300 0.010
0.012
3:054
5.005
0.51 0.02 0.03
0.06
005
0.06
0.07
0.0
1st Laplacian Eigenvector
(a)
(b)
Figure 4: Spectral embedding of sparsified graph for (a) Nested dataset and (b) Rings dataset, re- spectively.
If - E & A for all i, then for of = E[ - E(X,}
-72
Pr Pr 2 E + 5
2f + 2/4/33
40


Reviews:
[{"review_id": "cEckEAn7eU", "reviewer": "Reviewer_vw4c", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The primary focus of the paper is in the use of KDE queries to speed up many of the primitive operations applied to kernel matrices (Gram matrices) for the purpose of speeding up approximate executions of important linear-algebraic and graph theoretic operations/computations. The fundamental idea is that KDE queries allow for the approximation of the sums of kernel evaluations between an entire data set and a given query point in time independent of data set size and pre-processing which is linear in the data set size*.\n Algorithms are given for how these primitive operations are subsequently applied to important applications in linear algebra and computational graph theory.\n The paper concludes with an experiment illustrating the potential of the proposed method for the ubiquitous Singular Value Decomposition operation.", "summary_of_the_review": "Overall the paper is reasonably well written and includes potentially influential and important methodology in very relevant areas in machine learning. The theory is comprehensive enough to likely even allow readers to fairly easily extend the general method to operations/problems not explicitly in the scope of the paper.\n\nMy main concern relates to the fact that in many practical situations, the important \"parameter\" tau depends on sample size or likely is very small in practice, but this point seems to be almost glossed over in the paper (as well as the related literature). I think the contribution is significant enough that such a suppression (I do not intend to imply that it is deliberate) is unnecessary, and it is important to clarify/discuss where this is and is not an important factor. We see in the experiment that the performance is good FOR THE GIVEN GRAM MATRIX, however it may be that the scaling used may have oversmoothed the structure in the data to the extent that much of the local information is lost. Again, I do not mean to imply this is the case, only that as a reader one does not know.\n\nBelow I will conclude with some smaller questions/clarifications/suggestions/issues:\n- Using the term sub-linear suggests to a reader, in my opinion, sublinear in the data set size and not sublinear in the size of the Gram matrix. This is unnecessary.\n- Definition 1.1 does not seem to be correctly stated. Surely it should read that the query point y satisfies k(x, y) \\geq \\tau for all x \\in X\n- Theorem 1.3 requires setting/defining \\epsilon\n- Theorem 1.4 needs either to state that v is unit norm, or the LHS in the inequality should be divided by ||v||^2.\n- Sometimes one sees \"... with probability a/b...\" and other times \"... with probability d%...\". Why the inconsistency?\n- Below Theorem 1.5: \\omega has not been defined\n- The statement of Theorem 1.6 is rather clumsy, and could be made clearer. Also, \\epsilon features in the number of queries without being defined and also not present in any of the other terms. I presume it is somehow hidden in a probability somewhere, but the statement is \"with certainty\" and no probability is mentioned.\n- In the paragraph below Theorem 1.6, I presume the authors mean they compute the top few eigenVECTORS of the Laplacian\n- The notation of the edge-vertex incident matrix is quite confusing. Should the subscript x_i not be i? Also, although not ambiguous, subscripting with the edge is not mathematically precise. Is there a way to maybe define H as a map rather than a matrix?\n Further the expression  pi >= (1-\\epsilon)||r_i||^2 = 2(1-\\epsilon)k(x_i, x_j) is confusing since index j has not been defined. I presume first you set j and then can sample the rows with such probability?\n- Am I totally out of date by thinking the Nystrom approximation is a relevant alternative? How does the empirical performance on the experiment in Section 2 compare when using this ubiquitous approximation method?\n- The authors say that in the figure in the experiment, all the points in the true and approximate squared row norms lie close to the y=x line. This is very questionable, especially in the Glove data. Is it not important that lots of the approximate row norms are essentially zero when their actual values are far from that?", "strength_and_weaknesses": "Strengths:\n- The paper includes a fairly extensive coverage of the important operations required in popular machine learning applications.\n- The content is comprehensive enough that a reader can conceivably extend the algorithms given to any (most?) operations not covered explicitly in the paper.\n- The presentation is reasonably good overall, but I will provide some comments on this going forward.\n- The approach taken in the paper is a natural application of the KDE queries to the problems covered, which makes much of the content intuitive and easier than it would otherwise be to follow.\n\nWeaknesses:\n- There are some issues with clarity (I will provide some examples in the following)\n- I have some, what I believe to be important qualms with the way in which the content is presented. In particular, although there is no apparent incorrectness in the paper, the statements are arguably overblown in respect of their practical relevance. For example, as my closest connections to the content in the paper are from the points of view of spectral clustering and non-parametric smoothing using kernels, it is, in my experience, extremely rare that the scaling parameter (\\sigma in the paper) used in determining the kernel weights/similarities does not depend at least implicitly on the data set size (n). For example:\n-- For most non-parametric smoothing problems one has \\sigma ~ scale(X) n^(-1/(4+d)), where scale(X) is a measure of the scale of the data set. Using this, and for example the Gaussian kernel, one has \\tau ~ exp(-c*n^(2/(4+d))) for c = diam(X)^2/2scale(X)^2, where diam(X) is the diameter of the data set.\n NB: I see that the authors mention that some literature states that often the smallest element in the Gram matrix is a fixed constant, however that does not seem to be consistent with the kernels suggested in the paper unless the data distribution is compact and the scaling parameter does not tend to zero. In addition, I imagine that in general for good performance the bound on tau needs to be very small, and then the question is \"For what sort of data set size do we actually start to realise a substantial improvement in run-time while still maintaining accuracy?\"\n\n\nLet me clarify that I in no way think this invalidates the proposed approach, but I think it is a very important point which is worth discussing and, at the very least noting in the paper.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "ZpNFL0NJPY", "reviewer": "Reviewer_dkUW", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": " Using KDE queries (as black box), the paper proposes to effciently solve various gram matrix related problems such as spectral approximation/sparsification, simulating random walks, weighted sampling, etc. In each case, a theoretical bound on the approximation/quality along with no. KDE queries is presented. Simulation results on low rank approximation are provided.", "summary_of_the_review": "In view of the strong theory and nice presentation I tend to recommend an acceptance though I feel simulations section needs a heavy upgrading.\n\n--\nAfter seeing the rebuttal and other reviews, I increase my score.", "strength_and_weaknesses": "Strengths:\n1. The theoretical results in the theorems as well as the remarks are interesting. They also seem to compare well wrt. existing works.\n2. The writeup is very well organized by summarizing the theoretically results and proofs in an intuitive way. I appreciate the discussion in section 1.2\n\nWeakness:\n1. Though I understand the contribution is theoretical, the simulations seem  to be too restricitive. Only low rank based approximations are compared. Also, here, it would have been nice to see a comparison between time vs approximation, without which it is hard to evaluate the improvement.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "PXw2psncV2X", "reviewer": "Reviewer_Eddv", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This work studies broad linear algebraic and graph processing applications using the kernel density estimation (KDE) framework which provides sublinear-time algorithms in the number of vertices. A core approach is to sample vertices and edges with respect to the (kernel-based) weights using the KDE in sublinear time and construct a sparsified Laplacian, which guarantees a spectral approximation. Such a (spectrally) well-approximated Laplacian provides rigorous guarantees of downstream problems, such as low-rank approximation, graph clustering, and so on. Experimental results on low-rank approximation and graph sparsification are provided. ", "summary_of_the_review": "The results of this paper are strong and impactful. The paper is well-written and very clear.", "strength_and_weaknesses": "Strength:\n\n- The KDE is a fundamental problem in a wide range of areas and this work enhances the impact of the KDE by exploring various practical problems. Moreover, all problems preserve sublinear runtimes that are from the efficient KDE. \n- For efficient usage of the KDE, these works propose a multi-level KDE, which combines a tree structure. This incurs a logarithmic overhead hence sublinear-time advantage is preserved.\n- All the proposed methods use the KDE as a black-box approach, so further improvements on the KDE will automatically improve that of these tasks.\n- All applications have rigorous guarantees, which make the usage of the KDE reliable and much stronger.\n\nWeakness:\n\n- Due to the page limit, the main paper contains the overview only and all key contents are in the appendix. \n- Some of the applications are somewhat straightforward and trivial. So, I think it might be better if the authors reduce all applications to some specific ones, and make the main contributions thin that are not straightforward.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "Not applicable"}, "extra_text": {}}]

Author rebuttal:

Decision: Accept: notable-top-25%