# Sub-quadratic Algorithms for Kernel Matrices
via Kernel Density Estimation

Ainesh Bakshi
MIT
<EMAIL_ADDRESS>Work done in part while at Carnegie Mellon University Piotr
Indyk
MIT
<EMAIL_ADDRESS>Praneeth Kacham
CMU
<EMAIL_ADDRESS>Sandeep Silwal
MIT
<EMAIL_ADDRESS>Samson Zhou11footnotemark: 1
UC Berkeley and Rice University
<EMAIL_ADDRESS>
###### 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 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
weighted vertex and weighted edge sampling on kernel graphs, simulating random
walks on kernel graphs, and importance sampling 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 sparsification, where
we observe a 9x decrease in the number of kernel evaluations over baselines
for LRA and a 41x reduction in the graph size for spectral sparsification.

###### Contents

  1. 1 Introduction
    1. 1.1 Our Results
  2. 2 Technical Overview
    1. 2.1 Linear Algebra Applications
    2. 2.2 Graph Applications.
  3. 3 Further Related Works
    1. 3.1 Preliminaries
  4. 4 Algorithmic Building Blocks
    1. 4.1 Multi-level KDE
    2. 4.2 Weighted Vertex Sampling
    3. 4.3 Weighted Edge Sampling and Weighted Neighbor Edge Sampling
    4. 4.4 Random Walk
  5. 5 Linear Algebra Applications
    1. 5.1 Spectral Sparsification
      1. 5.1.1 Solving Laplacian Systems Approximately
    2. 5.2 Low-rank Approximation of the Kernel Matrix
    3. 5.3 Approximating the Spectrum in EMD
    4. 5.4 First Eigenvalue and Eigenvector Approximation
  6. 6 Graph Applications
    1. 6.1 Local Clustering
    2. 6.2 Spectral Clustering
    3. 6.3 Arboricity Estimation
    4. 6.4 Computing the Total Weight of Triangles
  7. 7 Empirical Evaluation
    1. 7.1 Results
  8. 8 Auxiliary Inequalities

##  1 Introduction

For a kernel function
$k:\mathbbm{R}^{d}\times\mathbbm{R}^{d}\rightarrow\mathbbm{R}$ and a set
$X=\\{x_{1}\ldots x_{n}\\}\subset\mathbbm{R}^{d}$ of $n$ points, the entries
of the $n\times n$ kernel matrix $K$ are defined as $K_{i,j}=k(x_{i},x_{j})$.
Alternatively, one can view $X$ as the vertex set of a complete weighted graph
where the weights between points are defined by the kernel matrix $K$. Popular
choices of kernel functions $k$ include the Gaussian kernel, the Laplace
kernel, exponential kernel, etc; see [SSB+02, STC+04, HSS08] for a
comprehensive overview.

Despite their wide applicability, kernel methods suffer from drawbacks, one of
the main being efficiency – given $n$ input points in $d$ dimensions, many
kernel-based algorithms need to materialize the full $n\times n$ kernel matrix
$K$ before performing the computation. For some problems this is unavoidable,
especially if high-precision results are required [BIS17]. In this work, we
show that we can in fact break this $\Omega(n^{2})$ barrier for several
fundamental problems in numerical linear algebra and graph processing. We
obtain algorithms that run in $o(n^{2})$ time and scale inversely-proportional
to the smallest entry of the kernel matrix. This allows us to skirt several
known lower bounds, where the hard instances require the smallest kernel entry
to 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
a fixed constant [MXB15, SRB+19, BIW19, BIMW21, KAP22]. We build on recently
developed fast approximate algorithms for Kernel Density Estimation [CS17,
BCIS18, SRB+19, BIW19, CKNS20]. Specifically, these papers present fast
approximate data structures with the following functionality:

######  Definition 1.1 (Kernel Density Estimation (KDE) Queries).

For a given dataset $X\subset\mathbbm{R}^{d}$ of size $n$, kernel function
$k$, and precision parameter $\varepsilon>0$, a KDE data structure supports
the following operation: given a query $y\in\mathbbm{R}^{d}$, return a value
$\textup{KDE}_{X}(y)$ that lies in the interval
$[(1-\varepsilon)z,(1+\varepsilon)\,z]$, where $z=\sum_{x\in X}k(x,y)$,
assuming that $k(x,y)\geqslant\tau$ for all $x\in 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 a preprocessing time that is sub-quadratic (in $n$), KDE data
structures use time sublinear in $n$ to answer queries defined as above. Note
that for all of our kernels, $k(x,y)\leqslant 1$ for all inputs $x,y$.

Table 1:  Instantiations of KDE queries. The query times depend on the dimension $d$, accuracy $\varepsilon$, and lower bound $\tau$. The parameter $\beta$ is assumed to be a constant. Type | $k(x,y)$ | Preprocessing Time | Query Time | Reference
---|---|---|---|---
Gaussian | $e^{-\|x-y\|_{2}^{2}}$ | $\frac{nd}{\varepsilon^{2}\tau^{0.173+o(1)}}$ | $\frac{d}{\varepsilon^{2}\tau^{0.173+o(1)}}$ | [CKNS20]
Exponential | $e^{-\|x-y\|_{2}}$ | $\frac{nd}{\varepsilon^{2}\tau^{0.1+o(1)}}$ | $\frac{d}{\varepsilon^{2}\tau^{0.1+o(1)}}$ | [CKNS20]
Laplacian | $e^{-\|x-y\|_{1}}$ | $\frac{nd}{\varepsilon^{2}\tau^{0.5}}$ | $\frac{d}{\varepsilon^{2}\tau^{0.5}}$ | [BIW19]
Rational Quadratic | $\frac{1}{(1+\|x-y\|_{2}^{2})^{\beta}}$ | $\frac{nd}{\varepsilon^{2}}$ | $\frac{d}{\varepsilon^{2}}$ | [BCIS18]

###  1.1 Our Results

We show that given a KDE data structure as described above, it is possible to
solve a variety of matrix and graph problems in time subquadratic time
$o(n^{2})$, 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 such algorithms for problems such as eigenvalue/eigenvector estimation,
low-rank approximation, graph sparsification, local clustering, aboricity
estimation, and estimating the total weight of triangles.

Our results are obtained via the following two-pronged approach. First, we use
KDE data structures to design algorithms for the following basic primitives,
frequently used in sublinear time algorithms and property testing:

  1. 1.

sampling vertices by their (weighted) degree in $K$ (Theorems 4.7 and 4.9 and
Algorithms 4.3 / 4.6),

  2. 2.

sampling random neighbors of a given vertex by edge weights in $K$ and
sampling a random weighted edge (Theorem 4.12 and Algorithms 4.11 and 4.13),

  3. 3.

performing random walks in the graph $K$ (Theorem 4.15 and Algorithm 4.16),
and

  4. 4.

sampling the rows of the edge-vertex incident matrix and the kernel matrix
$K$, both with probability proportional to respective row norms squared
(Section 5.1, Theorem 5.3, and Section 5.2, Corollary 5.14 respectively).

In the second step, we use these primitives to implement a host of algorithms
for the aforementioned problems. We emphasize that these primitives are used
in a 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 3.1 and Section 3.1. At a high level, many of our
applications, such as spectral sparsification, are succinctly characterized by
the following parameterization.

######  Parameterization 1.2.

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 at least
$\tau$.

Our applications derived from the basic graph primitives above can be
partitioned into two overlapping 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 specific
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 $\varepsilon$ and $\tau$.

Table 2:  Summary of linear algebra and graph applications for KDE subroutines. We suppress dependence on the precision $\varepsilon$. In spectral/local clustering and low-rank approximation, $k$ denotes the number of clusters and the rank of the approximation desired, respectively. The parameter $\phi$ refers to the quality of the underlying clusters; see Section 6.1. Problem |  $\\#$ of KDE Queries | Post-processing time | Prior Work
---|---|---|---
Spectral sparsification (Thm. 1.3) | $\widetilde{O}\left(\frac{n}{\tau^{3}}\right)$ | $O\left(\frac{nd}{\tau^{3}}\right)$ |  Remark 3.1
Laplacian system solver (Thm. 1.3) | $\widetilde{O}\left(\frac{n}{\tau^{3}}\right)$ | $O\left(\frac{nd}{\tau^{3}}\right)$ |  Remark 3.1
Low-rank approx. (Thm. 1.6) | $O(n)$ | $O\left(n\cdot\mathrm{poly}\left(k\right)+nkd\right)$ |  Remark 3.3
Eigenvalue Spectrum approx. (Thm. 1.4) | $\widetilde{O}(1/\tau)$ | $O(d/\tau)$ | $\Omega(n^{2}d)$
Approximating 1st Eigenvalue (Thm. 1.5) |  Remark 3.2 | $d\cdot\mathrm{poly}(1/\tau)$ |  $\omega(n)$ (Remark 3.2)
Local clustering (Thm. 1.7) | $\widetilde{O}\left(\mathrm{poly}(k)\cdot\frac{1}{\phi^{2}}\frac{\sqrt{n}}{\tau^{1.5}}\right)$ | $\widetilde{O}\left(\mathrm{poly}(k)\cdot\frac{1}{\phi^{2}}\frac{\sqrt{n}}{\tau^{1.5}}\right)$ |  Remark 3.6
Spectral clustering (Thm. 6.12) | $\widetilde{O}\left(\frac{n}{\tau^{2}}\right)$ | $O\left(\frac{nd}{\tau^{2}}\right)+\widetilde{O}\left(nk\right)$ |  Remark 3.6
Arboricity estimation (Thm. 1.9) | $\widetilde{O}\left(\frac{n}{\tau}\right)$ | $\widetilde{O}\left(\frac{n^{2}}{\tau}\right)$ | $\widetilde{O}(n^{3})+O(n^{2}d)$
Triangle estimation (Thm. 1.10) | $\widetilde{O}\left(\frac{1}{\tau^{3}}\right)$ | $\widetilde{O}\left(\frac{1}{\tau^{3}}\right)$ | $\Omega(n^{2}d)$

One of our main results is spectral sparsification of the kernel matrix $K$
interpreted as a weighted graph. In Section 5.1, we compute a 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 $D$ is a diagonal matrix of vertex degrees. The
Laplacian matrix encodes fundamental combinatorial properties of the
underlying graph and has been well-studied for numerous applications,
including sparsification; see [Mer94, BSST13, Spi16] for a survey of the
Laplacian and its applications. Our result computes a 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.2.

######  Theorem 1.3 (Informal; see Thm. 5.3).

Let $L$ be the Laplacian matrix corresponding to the graph $K$. Then, for any
$\varepsilon\in(0,1)$, there exists an algorithm that outputs a weighted graph
$G^{\prime}$ with only $m=O({n\log n}/{(\varepsilon^{2}\tau^{3})})$ edges,
such that with probability at least $9/10$, $(1-\varepsilon)L\preceq
L_{G^{\prime}}\preceq(1+\varepsilon)L$. The algorithm makes
$\widetilde{O}(m/\tau^{3})$ KDE queries and requires
$\widetilde{O}(md/\tau^{3})$ post-processing time.

We compare our results with prior works in Remark 3.1. We also show that
Parameterization 1.2 is inherent for spectral sparsification. In particular,
we use a hardness result from [ACSS20] to show that for the Gaussian kernel,
under the strong exponential time hypothesis [IP01], any algorithm that
returns an $O(1)$-approximate spectral sparsifier with $O(n^{1.99})$ edges
requires $\Omega\left(n\cdot 2^{\log(1/\tau)^{0.32}}\right)$ time (see Theorem
5.7 for a formal statement). Obtaining the optimal dependence on $\tau$
remains an outstanding open question, even for Gaussian and Laplace kernels.
Spectral sparsification has further downstream applications in solving
Laplacian linear systems, which we present in Section 5.1.1.

Continuing the theme of the Laplacian matrix, in Section 5.3, we also obtain a
succinct 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. (2)), 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 a specific graph generative
process [CKSV18].

######  Theorem 1.4 (Informal; see Theorem 5.17).

Let $\varepsilon\in(0,1)$ be the error parameter and $L$ be the normalized
Laplacian of the kernel graph $K$. Let
$\lambda_{1}\geqslant\lambda_{2}\ldots\geqslant\lambda_{n}$ be the eigenvalues
of $L$ and let $\lambda$ be the resulting vector. Then, there exists an
algorithm that uses
$\widetilde{O}\left(\exp\left(1/\varepsilon^{2}\right)/\tau\right)$ KDE
queries and $\exp\left(1/\varepsilon^{2}\right)\cdot d/\tau$ post-processing
time and outputs a vector $\widetilde{\lambda}$ such that with probability
$99/100$,
$\textrm{EMD}\left(\lambda,\widetilde{\lambda}\right)\leqslant\varepsilon$.

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 $\Omega(n^{2}d)$. Next, we obtain truly sublinear time algorithms for
approximating the top eigenvalue and eigenvector of the kernel matrix, a
problem which was studied in [BIMW21]. Our result is the following theorem.
Our bounds, and those of prior work, depend on the parameter $p$, which refers
to the exponent of $\tau$ in the KDE query runtimes. For example for the
Gaussian kernel, $p\approx 0.173$. See Table 1 for other kernels.

######  Theorem 1.5 (Informal; see Theorem 5.22).

Given an $n\times n$ kernel matrix $K$ that admits a KDE data-structure with
query time $d/(\varepsilon^{2}\tau^{p})$ (Table 1), there exists an algorithm
that outputs a unit vector $v$ such that
$v^{T}Kv\geqslant(1-\varepsilon)\lambda_{1}(K)$ in time
$\min\left(\tilde{O}(d/(\varepsilon^{4.5}\tau^{4})),\tilde{O}(d/(\varepsilon^{9+6p}\tau^{2+2p}))\right)$,
where $\lambda_{1}(K)$ denotes the largest eigenvalue of $K$.

We discuss related works in Remark 3.2. In summary, the best prior result of
[BIMW21] had a runtime of $\Omega(n^{1+p})$ 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 Section 5.2.

######  Theorem 1.6 (Informal; see Cor. 5.14).

There exists an algorithm that outputs a rank $r$ matrix $B$ such that
$\|K-B\|_{F}^{2}\leqslant\|K-K_{r}\|_{F}^{2}+\varepsilon\|K\|_{F}^{2}$ with
probability $99\%$ where $K_{r}$ is the optimal rank-$r$ approximation of $K$.
It uses $n$ KDE queries and
$O(n\cdot\mathrm{poly}(r,1/\varepsilon)+nrd/\varepsilon)$ post-processing
time.

We give detailed comparisons between our results and prior work in Remark 3.3.
As a summary, [BCW20] obtain a relative error approximation with a running
time of $\widetilde{O}\left(nd\left(r/\varepsilon\right)^{\omega-1}\right)$,
where $\omega$ denotes the matrix multiplication constant, whereas our running
time is dominated by $O(nrd/\varepsilon)$ and we obtain only additive error
guarantees. Nevertheless, the algorithm we obtain, which builds upon the
sampling scheme of [FKV04], is a conceptually simpler algorithm than the
algorithm of [BCW20] and easier to empirically evaluate. Indeed, we implement
this algorithm in Section 7 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 a $k$-clusterable graph, formally introduced in Definition
6.4. Intuitively, it describes a graph whose vertices can be partitioned into
$k$ disjoint clusters with high-connectivity within clusters and relatively
sparse connectivity in-between clusters.

######  Theorem 1.7 (Informal; see Theorem 6.9).

Let $K$ be a $k$-clusterable kernel graph with clusters $V=\cup_{1\leqslant
i\leqslant k}V_{i}$. Let $U,W$ be one of (not necessarily distinct) clusters
$V_{i}$. Let $u,w$ be randomly chosen vertices in partitions $U$ and $W$ with
probability proportional to their degrees. There exists $c=c(\varepsilon,k)$
and an algorithm that uses
$\widetilde{O}(c(k,\varepsilon)\sqrt{n}/\tau^{1.5})$ KDE queries and post-
processing time, with the property that with probability $1-\varepsilon$, if
$U=W$ then the algorithm reports that $u$ and $w$ are in the same cluster and
if $U\neq W$, the algorithm reports that $u$ and $w$ are in different
clusters.

Our definitions for the local clustering result are adopted from prior
literature in property testing; see Remark 3.6 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 6.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.8 (Informal; see Theorem 6.13).

Let $L$ be the Laplacian matrix of the spectral sparsifier. There exists an
algorithm that can compute $(1+\varepsilon)$-approximations of the first $k$
eigenvectors of $L$ in time
$\widetilde{O}\left(kn/(\tau^{2}\varepsilon^{2.5})\right)$.

We also give algorithms for approximating the arboricity of a graph, which is
the density of the densest subgraph of the kernel graph (see exact definition
in Section 6.3).

######  Theorem 1.9 (Informal; see Theorem 6.15).

There exists an algorithm that uses $m=\widetilde{O}(n/(\varepsilon^{2}\tau))$
KDE queries and $O(mn)$ post-processing time and outputs a sparse subgraph
$G^{\prime}$ of the kernel graph such that with high probability,
$(1-\varepsilon)\alpha_{G}\leqslant\alpha_{G^{\prime}}\leqslant(1+\varepsilon)\alpha_{G}$,
where $\alpha_{G}$ 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 $\widetilde{O}(nm)$ where $m$ is the number of edges
leading to a runtime of $\widetilde{O}(n^{3})+O(n^{2}d)$ [GGT89]. In Section
6.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 a
triangle as the product of its edge weights. This is a natural definition if
weighted edges are interpreted as parallel unweighted edges, in addition to
having applications in defining cluster coefficients of weighted graphs [KH06,
LLL07, AT08]. Our bound is similar in spirit to the bound of the unweighted
case given in [ELRS17], under a different computation model. We refer to
Remark 3.7 for discussions on related works.

######  Theorem 1.10 (Informal; see Theorem 6.17).

There exists an algorithm that makes $\widetilde{O}(1/\tau^{3})$ KDE queries
and the same bound for post-processing time and with probability at least
$\frac{2}{3}$, outputs a $(1\pm\varepsilon)$-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 [PT18, PT20a,
Tai22]. We view this direction of work as orthogonal to our line of study.
Lastly, the work [BIMW21] 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, [BIMW21] also study the problem of
estimating the sum of all entries in the kernel matrix and obtain tight bounds
for the latter.

##  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 a kernel
matrix $K$ without computing each of its entries $k(x_{i},x_{j})$.

##### Algorithmic Building Blocks

. We first describe the “multi-level” KDE data structure, which constructs 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 $O(\log n)$ levels, the overall runtime to
initialize our multi-level KDE data structure incurs only a logarithmic
overhead (see Figure 1 for an illustration).

##### Weighted vertex sampling.

We describe how to sample vertices approximately proportional to their
weighted degree, where the weighted degree of a vertex $x_{i}$ with $i\in[n]$
is $w_{i}=\sum_{j\neq i}k(x_{i},x_{j})$. We observe that performing $n$ KDE
queries suffices to get an 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 no query cost. Moreover, this preprocessing step
of taking $n$ queries only needs to be performed once. Further, we can then
perform weighted vertex sampling from a distribution that is
$\varepsilon$-close in total variation to the true distribution (see Theorem
4.9 for details). Here, we use a multi-level tree structure to iteratively
choose a 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 $O(\log n)$ runtime, which is significantly better than the naïve
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
$x$. The goal of weighted neighbor edge sampling is to efficiently output a
vertex $v$ such that
$\mathbf{Pr}[v=x_{k}]=\frac{(1\pm\varepsilon)k(x,x_{k})}{\sum_{j\in[n],x_{j}\neq
x}k(x,x_{j})}$ for all $k\in[n]$. Unlike the degree case, edge sampling 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 a KDE query provides. However, we can utilize a similar tree
procedure as in Figure 1 in conjunction with KDE queries.

In particular, consider the tree in Figure 1 where each internal node
corresponds to a subset of neighbors of $x$. The two children of a 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 as in the vertex sampling case: at every
node, we 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 $O(\log n)$ KDE queries and $O(\log n)$ time (see
theorems 4.12 and 4.14 for details).

##### 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 $v$
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(T\log n)$ 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 to the importance
sampling score / leverage score (see Definition 5.4). We remark that
approximately sampling proportional to the leverage score distribution for $H$
is a fundamental algorithmic primitive in spectral graph theory and numerical
linear algebra. We note that apriori, such a task seems impossible to perform
in $o(n^{2})$ time, even if the leverage scores are precomputed for us, since
the support of the distribution has size $\Theta(n^{2})$. 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 to approximate
leverage score sampling (see proof of Theorem 5.3). Further, we show that $H$
has a condition number (Lemma 5.6) that is bounded by $\mathrm{poly}(1/\tau)$.
Recall, the edge-vertex incident matrix is defined as the $\binom{n}{2}\times
n$ matrix with the rows indexed by all possible edges and the columns indexed
by vertices. For each $e=\\{i,j\\}$, we have
$H_{\\{i,j\\},i}=\sqrt{k(x_{i},x_{j})}$ and
$H_{\\{i,j\\},j}=-\sqrt{k(x_{i},x_{j})}$. 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 $H^{T}H=L_{G}$ where $L_{G}$ is the Laplacian
matrix of the graph corresponding to the kernel matrix $K$. For both $H$ and
$K$, we wish to sample the rows with probability proportional to row normed
squared. For example, the row $r_{e}$ corresponding to edge $e=(x_{i},x_{j})$
in $H$ satisfies $\|r_{e}\|_{2}^{2}=2k(x_{i},x_{j})$. Since the squared norm
of 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 to 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$. Naïvely performing this sampling would require us to implicitly
compute the entire kernel matrix, which as mentioned previously, is
prohibitive. However, if there exists a constant $c$ such that the kernel
function $k$ that defines the matrix $K$ satisfies $k(x,y)^{2}=k(cx,cy)$ for
all inputs $x,y$, then the $\ell_{2}^{2}$ norm of each row can be approximated
via a KDE query on the transformed dataset $X^{\prime}=cX$. In particular, the
$\ell_{2}^{2}$ row norms of $K$ are the vertex degrees of the kernel graph for
$X^{\prime}$. The property that $k(x,y)^{2}=k(cx,cy)$ 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.

###  2.1 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 $\varepsilon$ spectral
sparsifier for the kernel density graph $G$ can be computed i.e., we compute a
graph $G^{\prime}$ such that for all vectors $x$,
$(1-\varepsilon)x^{T}L_{G}x\leqslant
x^{T}L_{G^{\prime}}x\leqslant(1+\varepsilon)x^{T}L_{G^{\prime}}x$, where
$L_{G}$ and $L_{G^{\prime}}$ denote the Laplacian matrices of the graphs $G$
and $G^{\prime}$. Recall that $H$ is the $\binom{n}{2}\times n$ matrix such
that $H_{\\{i,j\\},i}=\sqrt{k(x_{i},x_{j})}$ and
$H_{\\{i,j\\},j}=-\sqrt{k(x_{i},x_{j})}$. 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^{T}H=L_{G}$. It
is known that sampling $t=O(n\log(n)/\varepsilon^{2})$ rows of the matrix $H$
by using the so-called leverage scores gives a $t\times\binom{n}{2}$
selecting-and-scaling matrix $S$ such that with probability at least $9/10$,

$\displaystyle(1-\varepsilon)L_{G}=(1-\varepsilon)H^{T}H\preceq
H^{T}S^{T}SH\preceq(1+\varepsilon)H^{T}H=(1+\varepsilon)L_{G}.$ (1)

Thus the matrix $SH$ directly corresponds to a graph $G^{\prime}$, which is an
$\varepsilon$ 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 perform leverage score sampling of $H$. On the other hand,
observe that the squared norm of row $\\{i,j\\}$ of $H$ is $2k(x_{i},x_{j})$
and with an 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 to its squared norm. It is
a standard result that sampling from squared length distribution gives a
selecting-and-scaling matrix $S$ that satisfies (1), although we have to
sample $t=O(\kappa^{2}n\log(n)/\varepsilon^{2})$ rows from this distribution,
where $\kappa=\sigma_{\max}(H)/\sigma_{\min}(H)$ denotes the condition number
of $H$ ($\sigma_{\max}(H)$/$\sigma_{\min}(H)$ denote the largest/smallest
positive singular values).

With the parameterization that for all $i\neq j$,
$k(x_{i},x_{j})\geqslant\tau$, we are able to show that $\kappa\leqslant
O(1/\tau^{1.5})$. 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 a Cheeger-type inequality for
weighted graphs. Note that
$\sigma_{\min}(H)\geqslant\sqrt{\lambda_{2}(H^{T}H)}=\sqrt{\lambda_{2}(L_{G})}$,
where we use $\lambda_{2}(M)$ to denote the second smallest eigenvalue of a
positive semidefinite matrix. Cheeger’s inequality lower bounds exactly the
quantity $\lambda_{2}(L_{G})$ in terms of graph conductance. A lower bound of
$\tau$ on every kernel value implies that every node in the Kernel Graph has a
high weighted degree and this lets us lower bound $\lambda_{2}(G)$ in terms of
$\tau$ using a Cheeger-type inequality from [FN02] and shows that
$O(n\log(n)/\tau^{3}\varepsilon^{2})$ samples from the approximate squared
length sampling distribution gives an $\varepsilon$ spectral sparsifier for
the graph $G$.

##### First eigenvalue and eigenvector approximation.

Our goal is to compute a $1-\varepsilon$ approximation to $\lambda$, 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 to at least $n\tau$, a condition
used in prior works [BIMW21], the first eigenvalue must be at least $n\tau$ by
looking at the quadratic form associated with the all-ones vector.

Now we combine two disparate families of algorithms: first the guarantees of
[BMR21, BCJ20] show that sub-sampling a $t\times t$ principal submatrix of a
PSD matrix preserves the eigenvalues of the matrix up to an additive
$O(n/\sqrt{t})$ factor. Since we’ve shown the first eigenvalue of $K$ is at
least $n\tau$, we can set $t$ roughly $O(1/(\varepsilon^{2}\tau^{2}))$ with
the guarantee that the top eigenvalue of the sub-sampled matrix is at lest
$(1-\varepsilon)\lambda$. Now we can either run the standard Krylov method
algorithm [MM15] to compute the top eigenvalue of the sampled matrix or
alternatively, we can instead use the algorithm of [BIMW21], the prior state
of the art, to compute the eigenvalues of the sampled matrix. At a 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 $n$, the size of
the dataset, whereas the prior state of the art result of [BIMW21] have a
$\omega(n)$ runtime.

##### Kernel matrix low-rank approximation.

In this setting, our goal here is to output a matrix $B$ such that

$\|K-B\|_{F}^{2}\leqslant\|K-K_{t}\|_{F}^{2}+\varepsilon\|K\|_{F}^{2}$

where $K_{t}$ is the best rank-$t$ approximation to the kernel matrix $K$. The
efficient algorithm of [FKV04] is able to achieve this guarantee if one can
sample the $i$th row $r_{i}$ of $K$ with probability
$p_{i}\geqslant\Omega(1)\cdot\|r_{i}\|_{2}^{2}/\|K\|_{F}^{2}$. We can perform
such an action using our primitive, which is capable of sampling the rows of
$K$ with probability proportional to the squared row norms for the Laplacian,
exponential, and Gaussian kernels. Thus for these kernels, we can immediately
obtain efficient algorithms for computing a low-rank approximation.

##### Spectrum approximation.

For this problem, the goal is to compute approximations of 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 [CKSV18] 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 [CKSV18] on
kernel graphs in a black-box manner.

##### Spectral clustering.

Given our spectral sparsification result, we can immediately obtain a fast
version of a heuristic algorithm used in practice for graph clustering: we
embed each vertex into $\mathbbm{R}^{k}$ using $k$ eigenvectors of the
Laplacian matrix and run $k$-means clustering. Clearly if we have a sparse
graph, the eigenvector computation is faster. Theoretically, we can show that
spectral sparsification preserves a notion of clusterability which is weaker
than the definition used in the local clustering section and we additionally
give empirical evidence of the validity of this procedure.

###  2.2 Graph Applications.

We now discuss our graph applications.

##### Local clustering.

The random walks primitive allow us to run a well-studied local clustering
algorithm on the kernel graph. The algorithm is quite standard in the property
testing literature (see [CPS15] and [Pen20]) 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 of the algorithm is to
determine if two input vertices $u$ and $v$ belong to the same cluster of the
kernel graph if the graph has a natural cluster structure (see Definition 6.4
for the formal definition). The well-studied algorithm in literature performs
approximately $O(\sqrt{n})$ random walks from $u$ and $v$ 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 $\ell_{2}$ distance which can be detected via a
standard distribution tester of [CDVV14]. The guarantees of the overall local
clustering algorithm of [CPS15] follow for kernel graphs since we only need to
access the graph via random walks.

##### Arboricity estimation.

The arboricity of a weighted graph $G=(V,E,w)$ is defined as
$\alpha:=\max_{U\subseteq V}\frac{w(E(G_{U}))}{|U|}$. Informally, the
arboricity of a (weighted) graph represents the maximum (weighted) density of
a subgraph of $G$. To approximate the weighted arboricity, we adapt a result
of [MTVV15], who observed that to estimate the arboricity on unweighted
graphs, it suffices to sample a set of
$\tilde{O}\left(\frac{|V|}{\varepsilon^{2}}\right)$ 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 [MTVV15] 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.e., importance sampling. In fact, a similar result still holds if we
only have upper bounds on the weight of each edge, provided that we increase
the number of fixed edges that we sample by the gap between the upper bound
and the actual weight of the edge. Thus, our arboricity algorithm requires
sampling a fixed number of edges, where each edge is sampled with probability
proportional 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 naïve approach of querying each edge in
the kernel graph by using our weighted edge sampling subroutine to sample a
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.

##### Weighted triangle estimation.

We define the weight of a triangle as the product of its edges, generalizing
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 [ELRS17], 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 $R$ involves sampling neighbors of the
vertices of $e$, which we can efficiently implement using our weighted
neighbor edge sampling subroutine.

##  3 Further Related Works

######  Remark 3.1.

Spectral sparsification for kernel graphs has also been studied in prior
works, notably in [ACSS20] and [Qua21]. We first compare to [ACSS20], who
obtain a spectral sparsification using an entirely different approach. They
obtain an almost linear time sparsifier ($n^{1+o(1)}$) when the kernel is
multiplicativily Lipschitz (see Section 1.1.2 in [ACSS20] for definition) and
show hardness for constructing such a sparsifier when it is not. Focusing on
the Gaussian kernel, under Parameterization 1.2, [ACSS20] obtain an algorithm
that runs in time $O\left(nd+n2^{\sqrt{\log(1/\tau)\log n}\log(\log
n)}/\varepsilon^{2}\right)$, whereas our algorithm runs in
$O\left(nd\log^{2}(n)/(\varepsilon^{2}\tau^{2.0173+o(1)})\right)$ time. We
also note that the dimension $d$ can be upper bounded by $O(\log
n/\varepsilon^{2})$ by applying Johnson-Lindenstrauss to the initial dataset.
Therefore, [ACSS20] obtain a better dependence on $1/\tau$, whereas we obtain
a better dependence on $n$. A similar comparison can be established for other
kernels as well. In practice, $\tau$ is set to be a small fixed constant,
whereas $n$ can be arbitrarily large. Indeed in practice, a common setting of
$\tau$ is $0.01$ or $0.001$, irrespective of the size of the dataset [MXB15,
SRB+19, BIW19, BIMW21, KAP22].

We now compare our guarantees to that of [Qua21]. The author studies spectral
sparsification resurrected to smooth kernels (for example kernels of the form
$1/\|x-y\|_{2}^{t}$ which have a polynomial decay; see [Qua21] for a formal
definition). This family does not include Gaussian, Laplacian, or exponential
kernels. For smooth kernels, [Qua21] obtained a sparsifier with a nearly
optimal $\tilde{O}(n/\varepsilon^{2})$ number of edges in time
$\tilde{O}(nd/\varepsilon^{2})$. Our algorithm obtains a similar dependence in
$n,d,\varepsilon$ but includes an additional $1/\tau^{3}$ factor. However, it
generalizes for any kernel supporting a KDE data structure, which includes
smooth kernels [BCIS18] (see Table 1 for a summary of kernels where our
results apply). Our techniques are also different: [Qua21] 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 sparsifier, another interesting
feature of [Qua21] is that it enriches the connections between spectral
sparsification of kernel graphs and KDE data structures. Indeed, the data
structures used in [Qua21] are inspired by and were used in the prior work of
[BCIS18] 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 [Qua21] for details.

######  Remark 3.2.

Our algorithm returns a sparse vector $v$ supported on roughly
$O(1/(\varepsilon^{2}\tau^{2}))$ coordinates. The best prior result is that of
[BIMW21] which presented an algorithm with total runtime
$O\left(\frac{dn^{1+p}\log(n/\varepsilon)^{2+p}}{\varepsilon^{7+4p}}\right)$.

In comparison, our bound has no dependence on $n$ and is thus a truly
sublinear runtime. Note that the bound of [BIMW21] does not depend on $\tau$.
We do not state the number of KDE queries used explicitly in Table 2 since our
algorithm uses KDE queries on a subsampled dataset and in addition, only uses
them by calling the algorithm of [BIMW21] as a subroutine (on the subsampled
dataset). The algorithm of [BIMW21] uses $\tilde{O}(1/\varepsilon)$ KDE
queries but with various different initialization of $\tau$ so it is not
meaningful to state “one” bound for the number of KDE queries used and thus
the final runtime is a more meaningful quantity to state. Lastly, the authors
in [BIMW21] present a lower bound of $\Omega(nd)$ for estimating the top
eigenvalue $\lambda_{1}$, which ostensibly seems at odds with our stated bound
which has no dependence on $n$. However, the lower bound presented in [BIMW21]
essentially sets $\tau=1/\mathrm{poly}(n)$ for a large polynomial factor
depending on $n$ (we estimate this factor to be $\Omega(n^{2})$). Since we
parameterize our dependence via $\tau$, which in practice is often set to a
fixed constant, we can bypass the lower bound.

######  Remark 3.3.

We now compare our low-rank approximation result with a recent work of [MW17,
BCW20]. They showed the following theorem:

######  Theorem 3.4 (Theorem 4.2, [BCW20]).

Given a $n\times n$ PSD matrix $A$, target rank $r\in[n]$ and accuracy
parameter $\varepsilon\in(0,1)$, there exists an algorithm that queries
$\widetilde{O}\left(nr/\varepsilon\right)$ entries in $A$ and with probability
at least $9/10$, outputs a rank-$r$ matrix $B$ such that

$\left\lVert A-B\right\rVert_{F}^{2}\leqslant(1+\varepsilon)\left\lVert
A-A_{r}\right\rVert_{F}^{2},$

where $A_{r}$ is the best rank-$r$ approximation to $A$. Further, the running
time is $\widetilde{O}\left(n\left(r/\varepsilon\right)^{\omega-1}\right)$,
where $\omega$ is the matrix multiplication constant.

We note that their result applies to kernel matrices as well via the following
fact.

######  Fact 3.5 (Kernel Matrices are PSD, [SSB+02]).

Let $k$ be a reproducing kernel and $X$ be $n$ data points in
$\mathbb{R}^{d}$. Let $K$ be the associated $n\times n$ kernel matrix such
that $K_{i,j}=k(x_{i},x_{j})$. Then, $K\succ 0$.

Here, the family of reproducing kernels is quite broad and includes polynomial
kernels, 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 [BCW20] have comparable
runtimes. While [BCW20] obtain relative-error guarantees, we only obtain
additive-error guarantees.

However, reading each entry of the kernel matrix require $O(d)$ time and thus
[BCW20] obtain an running time of
$\widetilde{O}\left(nd\left(r/\varepsilon\right)^{\omega-1}\right)$, whereas
our running time is dominated by $O(nrd/\varepsilon)$. We note that similar
ideas as our algorithm for additive error LRA were previously used to design
subquadratic algorithms running in time $o(n^{2})$ for low-rank approximation
of distance matrices [BW18, IVWW19].

######  Remark 3.6.

Our definitions for the local clustering result are adopted from prior
literature in property testing; see [KS08, CS10, GR11, CPS15, CKK+18, DPRS19,
Pen20, GKL+21] and the references within. Our algorithmic details for the
local cluster section are also derived from prior works, such as the works of
[CPS15] and [Pen20]; indeed, many of the lemmas of the local clustering
section follow in a straightforward fashion from [CPS15] and [Pen20]. 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 on local clustering requires having access to
the entire neighbor of a vertex when performing random walks, thereby implying
the runtime of $\Omega(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 a straightforward manner via our queries can be viewed a major
strength of our work.

######  Remark 3.7.

Our general bound of the number of KDE qeuries required to approximate the
total weight of triangles in Theorem 6.17 is
$\widetilde{O}(m\sqrt{w_{G}}/w_{T})$, where $w_{G}$ is the sum of all entries
of $K$ and $w_{T}$ is the total weight of triangles we wish to approximate.
This bound is a natural generalization of the result of [ELRS17]. 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 a 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 4 play a similar role to the queries used in [ELRS17].
There the authors give a bound of $\widetilde{O}(m^{3/2}/T)$ queries where $T$
is the total number of triangles. In our case, we indeed get a bound of the
order of $m^{3/2}$ in the numerator as $w_{G}\leqslant mw_{\max}$ and $w_{T}$
is the natural analogue of $T$ in [ELRS17]. Finally note that under our
parameterization of every edge in the kernel graph possessing weight at most
$1$ and at least $\tau$, our bound reduces to at most
$\widetilde{O}(1/\tau^{3})$ KDE queries.

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 $\Omega(n^{2}d)$ in our
setting.

We also note that our paper is closely related to the field of (graph)
property testing. In graph property testing, it is customary to assume query
access to an unknown graph via vertex and edge queries [Gol17]. 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^{2}$, the size of the graph. We can interpret the graph primitives we
construct as a realization of the property testing model where queries are
explicitly constructed.

###  3.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 polynomial in the
dimension $d$. Those data structures have construction times of the form
$O(dn/(\tau^{p}\varepsilon^{2}))$ and answer KDE queries in time
$O(d/(\tau^{p}\varepsilon^{2}))$, under the condition that for all queries $y$
we have $\frac{1}{n}\sum_{x\in X}k(x,y)\geqslant\tau$ (which clearly holds
under our Parameterization 1.2). The algorithms are randomized, 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 a random subset $R\subset X$ of size
$O(1/(\tau\varepsilon^{2}))$ and reports $\frac{n}{|R|}\sum_{x\in R}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 $\tau$, 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 $\sigma$ (for example, the
exponential kernel is of the form $k(x,y)=\exp(-\|x-y\|_{2}/\sigma)$). In
practice, a common choice of $\sigma$ follows the so called ‘median rule’
where $\sigma$ is set to be the median distance among all pairs of points in
$X$. Thus, according to the median rule, the ‘typical’ kernel values in the
graph $K$ are $\Omega(1)$. 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 $\tau$ is a small constant,
independent of the size of the dataset [KAP22].

We note that, in addition to the aforementioned algorithms with theoretical
guarantees, there are other practical algorithms based on random sampling,
space partition trees [GM01, GM03, LMG06, LG08, MSR+08, RLMG09, MXB15],
coresets [Phi13, ZJPL13, PT20b], or combinations of these methods [KAP22],
which support queries needed in Definition 1.1; see [KAP22] 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.

##  4 Algorithmic Building Blocks

$A_{1,n/4}$$\ldots$$\ldots$$\ldots$$A_{1,n/2}$$A_{n/2+1,n}$$A_{1,n}$ Figure 1:
Multi-level Kernel Density Estimation Data Structure.

###  4.1 Multi-level KDE

We first describe the “multi-level” KDE data structure, which is required in
our algorithms. The data structure recursively constructs a KDE data structure
on the entire dataset $X$, and then recursively partitions $X$ into two
halves, building a KDE data structure on each half. See Algorithm 4.1 for more
details.

######  Algorithm 4.1 (Multi-level KDE Construction).

Input: Dataset $X\subset\mathbbm{R}^{d}$, precision $\varepsilon>0$.
Operation: 1. Let $T=X$. 2. While $|T|>1$, (a) Construct $\textup{KDE}_{X}$
queries (see Definition 1.1). (b) Recursively apply Multi-level KDE
Construction to $T[1:\lfloor m/2\rfloor]$ and $T[\lfloor m/2\rfloor+1:m]$
Output: All the data structures associated with the KDE query constructions

######  Lemma 4.2.

Given a dataset $X\subset\mathbbm{R}^{d}$, suppose the initialization of the
KDE data structure defined in Definition 1.1 uses runtime $f(n,\varepsilon)$
for some function linear in $n$. Then the total construction time of Algorithm
4.1 is $f(n\log n,\varepsilon)$.

###### Proof.

The proof follows from the fact that at each recursive level, we do
$O(f(n,\varepsilon))$ total work since $f$ is linear in $n$ and there are
$O(\log n)$ levels. ∎

###  4.2 Weighted Vertex Sampling

We now discuss our fundamental primitives. The first one computes approximate
weighted degrees for all vertices. Algorithm 4.6 then performs vertex sampling
by their (weighted) degree.

######  Algorithm 4.3 (Computing Approximate (Weighted) Degrees).

Input: Dataset $X\subset\mathbbm{R}^{d}$, precision $\varepsilon>0$.
Operation: 1. For $i\in[1,n]$ : (a)
$p_{i}\leftarrow\textup{KDE}_{X}(x_{i})-(1-\varepsilon)\,k(x_{i},x_{i})$
Output: Reals $\\{p_{i}\\}_{i=1}^{n}$ such that
$(1-\varepsilon)\textup{deg}(x_{i})\leqslant
p_{i}\leqslant(1+\varepsilon)\textup{deg}(x_{i})$ for all $1\leqslant
i\leqslant n$

######  Definition 4.4 (Weighted Vertex Sampling).

The weighted degree of a vertex $x_{i}$ with $i\in[n]$ is $w_{i}=\sum_{j\neq
i}k(x_{i},x_{j})$. The goal of weighted vertex sampling is to output a vertex
$v$ such that
$\mathbf{Pr}[v=x_{i}]=\frac{(1\pm\varepsilon)w_{i}}{\sum_{j\in[n]}w_{j}}$ for
all $i\in[n]$.

This is a straightforward application of using $n$ 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 a preprocessing step that uses $O(n)$ queries upfront and then
allows for arbitrary access at any point in the future with no query cost.

Once we acquire $\\{p_{i}\\}_{i=1}^{n}$, we can perform a fast sampling
procedure through the following algorithm, which we state in slightly more
general terms.

######  Algorithm 4.5 (Sample from Positive Array).

Input: Array $A=[a_{1},\cdots,a_{n}]$ with $a_{i}>0$ for all $i$. Access to
queries $A_{i,j}=\sum_{i\in[t,j]}a_{t}$ for $1\leqslant i\leqslant j\leqslant
n$. Operation: 1. Let $T=A$. While $|T|>1$: (a) Let $m=\text{len}(T)$. (b) Let
$a\leftarrow\sum(T[1:\lfloor m/2\rfloor])$ //Can be simulated using an
$A_{i,j}$ query (c) Let $b\leftarrow\sum(T[\lfloor m/2\rfloor+1:m])$. (d) If
Unif$[0,1]\leqslant a/(a+b)$, $T\leftarrow T[1:\lfloor m/2\rfloor]$ (e) Else
$T\leftarrow T[\lfloor m/2\rfloor+1:m]$. Output: The single remaining element
in $T$

Combining Algorithms 4.3 and 4.5, we can sample from the degree distribution
of the graph $K$.

######  Algorithm 4.6 (Degree Sampling).

Input: Dataset $X\subset\mathbbm{R}^{d}$, precision $\varepsilon>0$.
Operation: 1. Use Algorithm 4.3 to compute reals $\\{p_{i}\\}_{i=1}^{n}$ such
that $(1-\varepsilon)\textup{deg}(x_{i})\leqslant
p_{i}\leqslant(1+\varepsilon)\textup{deg}(x_{i})$ for all $1\leqslant
i\leqslant n$ (only needs to be done once). 2. $i\leftarrow$ index in $[n]$,
which is the output of running Algorithm 4.5 on the array
$\\{p_{i}\\}_{i=1}^{n}$. Output: $x_{i}$ with probability
$p_{i}/\sum_{j\in[n]}p_{j}$.

We now analyze the correctness and the runtimes of the algorithms proposed in
Section 4. First, we give guarantees on Algorithm 4.3.

######  Theorem 4.7.

Algorithm 4.3 returns $\\{p_{i}\\}_{i=1}^{n}$ such that
$(1-\varepsilon)\textup{deg}(x_{i})\leqslant
p_{i}\leqslant(1+\varepsilon)\textup{deg}(x_{i})$ for all $1\leqslant
i\leqslant n$.

###### Proof.

The proof follows by the Definition of a KDE query, Definition 1.1. ∎

We now analyze Algorithm 4.5, which samples from an array based on a 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 4.

######  Lemma 4.8.

Algorithm 4.5 samples an index $i\in[n]$ proportional to $a_{i}$ in $O(\log
n)$ time with $O(\log n)$ queries.

###### Proof.

Consider the sampling diagram given in Figure 1. Algorithm 4.5 does the
following: it first queries the root node $A_{1,n}$ and then its two children
$A_{1,m},A_{m+1,n}$ where $m=\lfloor n/2\rfloor$. Note that
$A_{1,n}=A_{1,m}+A_{m+1,n}$. It then picks the tree rooted at $A_{1,m}$ with
probability $\frac{\sum_{i\in[m]}a_{i}}{\sum_{i\in[n]}a_{i}}$ and otherwise,
picks the tree rooted at $A_{m+1,n}$. The procedure recursively continues by
querying the root node, its two children, and picking one of its children to
be the new root node with probability proportional to the child’s weight given
by an appropriate query access. This is done until we reach a leaf node that
corresponds to an index $i\in[n]$.

We now prove correctness. Note that each node of the tree in Figure 1
corresponds to a subset $S\subseteq[n]$. We prove inductively that the
probability of landing on the vertex is equal to $\sum_{i\in S}a_{i}$. 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_{1},S_{2}$. We know that we are at node $S$ with probability $\sum_{i\in
S}a_{i}/\sum_{j}a_{j}$. Furthermore, we transition to $S_{1}$ with probability
$\sum_{i\in S_{1}}a_{i}/\sum_{j\in S}a_{j}$. Therefore, the probability of
being at $S_{1}$ is equal to

$\frac{\sum_{i\in S_{1}}a_{i}}{\sum_{j\in S}a_{j}}\cdot\frac{\sum_{i\in
S}a_{i}}{\sum_{j}a_{j}}=\frac{\sum_{i\in S_{1}}a_{i}}{\sum_{j}a_{j}}.$

Since there is only one path from the root node to any vertex of a tree, this
completes the induction.

The runtime and the number of queries taken follows from the fact that the
sampling procedure descends on a tree with $O(\log n)$ height. ∎

Combining Algorithms 4.3 and 4.5 allows us to sample from the degree
distribution of the graph $K$ up to low error in total variation (TV)
distance.

######  Theorem 4.9.

Algorithm 4.6 samples from the degree distribution of $K$ up to TV error
$O(\varepsilon)$ using a fixed overhead of $n$ KDE queries and runtime $O(\log
n)$.

###### Proof.

Since $p_{i}$ is with a $1\pm\varepsilon$ factor of $\textup{deg}(x_{i})$ for
all $i$, then $\\{p_{i}\\}_{i=1}^{n}$ is $O(\varepsilon)$ close in total
variation distance from the true degree distribution. Moreover, Algorithm 4.5
perfectly samples from the array $\\{p_{i}\\}_{i=1}^{n}$, which proves the
first part of the theorem.

For the second part, note that acquiring $\\{p_{i}\\}_{i=1}^{n}$ requires $n$
KDE queries. We can then construct the data structure for Algorithm 4.5 by
computing all the partial prefix sums in $O(n)$ time. Now the query access
required by Algorithm 4.5 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 4.5.
It follows from Lemma 4.8 that Algorithm 4.6 takes $O(\log n)$ time. ∎

###  4.3 Weighted Edge Sampling and Weighted Neighbor Edge Sampling

We describe how to perform weighted neighbor edge sampling.

######  Definition 4.10 (Weighted Neighbor Edge Sampling).

Given a vertex $x_{i}$, the goal of weighted neighbor edge sampling is to
output a vertex $v$ such that
$\mathbf{Pr}[v=x_{k}]=\frac{(1\pm\varepsilon)k(x_{i},x_{k})}{\sum_{j\in
n,j\neq i}k(x_{i},x_{j})}$ for all $i\in[n]$.

######  Algorithm 4.11 (Sample Random Neighbor).

Input: Dataset $X\subset\mathbbm{R}^{d}$, precision $\varepsilon>0$, input
vertex $x_{i}\in X$. Operation: 1. Let $\varepsilon^{\prime}=\varepsilon/\log
n$ and $T\leftarrow X\setminus\\{x_{i}\\}$. 2. While $|T|>1$ (a) Let
$m\leftarrow|T|$. (b) Compute
$a\leftarrow\textup{KDE}_{T[1:m/2],\varepsilon^{\prime}}(x_{i})$ and
$b\leftarrow\textup{KDE}_{T[m/2+1:m],\varepsilon^{\prime}}(x_{i})$. (c) If
$x_{i}\in T[1:m/2]$, set $a\leftarrow
a-(1-\varepsilon^{\prime})k(x_{i},x_{i})$. (d) If $x_{i}\in T[m/2+1:m]$, set
$b\leftarrow b-(1-\varepsilon^{\prime})k(x_{i},x_{i})$. (e) If
Unif$[0,1]\leqslant a/(a+b)$, let $T\leftarrow T[1:m/2]$. Else, let
$T\leftarrow T[m/2+1:m]$. Output: Return the last element $x\in T$ such that
$x\in X\setminus\\{x_{i}\\}$ and the probability of selecting $x$ is
proportional to $k(x_{i},x)$.

We now prove the correctness of Algorithm 4.11 based on the ideas in Lemma
4.8. Note that Algorithm 4.11 takes in input a precision level $\varepsilon$,
which can be adjusted and impacts the accuracy of KDE queries. We will discuss
the cost of initializing KDE queries with various precisions in Section 3.1.

######  Theorem 4.12.

Let $x_{i}\in X$ be an input vertex. Consider the distribution $\mathcal{D}$
over $X\setminus\\{x_{i}\\}$, the neighbors of $x_{i}$ in the graph $K$,
induced by the edge weights in $K$. Algorithm 4.11 samples a neighbor from a
distribution that is within TV distance $O(\varepsilon)$ from $\mathcal{D}$
using $O(\log n)$ KDE queries and $O(\log n)$ time. In addition, we can
perfectly sample from $\mathcal{D}$ using $O(\log n/\tau)$ additional kernel
evaluations in expectation.

###### Proof.

The proof idea is similar to that of Lemma 4.8. Given a vertex $x_{i}$, 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 a straightforward KDE query as an edge only cares about the
kernel value between two points, rather than the sum of kernel values that a
KDE query provides. Nevertheless, we can utilize the tree procedure outline in
the proof of Lemma 4.8 in conjunction with KDE queries with over various
subsets of $X$.

Imagine the same tree as in Figure 1 where each subset corresponds to a subset
of neighbors of $x_{i}$ (note that $x_{i}$ cannot be its own neighbor and
hence we subtract $k(x_{i},x_{i}$) in line $7$ or line $10$). Algorithm 4.11
descends down the tree using the same probabilistic procedure as in the proof
of Lemma 4.8: 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 is the sum of the weights of the edges connecting to the
corresponding neighbors of $x_{i}$.

Now compare the telescoping product of probabilities that lands us in some
leaf node $a_{j}$ to the ideal telescoping product if we knew the exact array
of edge weights as in the proof of Lemma 4.13. Suppose the tree has height
$\ell$. 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 a factor of $1+\varepsilon^{\prime}$ or
$1-\varepsilon^{\prime}$ factor respectively.

Therefore, we land in the correct leaf node with the same probability as in
the ideal telescoping product, except our probability can be off by a
multiplicative $(1\pm\varepsilon^{\prime})^{\ell}$ factor. However, since
$\varepsilon^{\prime}=\varepsilon/\log n$ and $\ell\leqslant\log n$, this
factor is within $1\pm\varepsilon$. Thus, we sample from the correct
distribution over the leaves of the trees in Figure 1 up to TV distance
$O(\varepsilon)$. Now by doing $O(1/\tau)$ steps of rejection sampling, we can
actually get a prefect sample of the edge. This is because the denominator of
the fraction for $\mathbf{Pr}[v=x_{k}]$ is at least $\Omega(n\tau)$ and at
most $n$ so we can estimate the proportionality constant in the denominator by
$n$ which is only at most $O(1/\tau)$ multiplicative factor larger. Hence by
standard guarantees of rejection sampling, we only need repeat the sampling
procedure $O(1/\tau)$ additional times. ∎

######  Algorithm 4.13 (Sample Random Edge by Weight).

Input: Dataset $X\subset\mathbbm{R}^{d}$, precision $\varepsilon>0$.
Operation: 1. Compute $x_{i}\leftarrow$ random vertex by using Algorithm 4.6.
2. Compute $x_{j}\leftarrow$ random Neighbor of $x_{i}$ using Algorithm 4.11.
Output: Edge $(x_{i},x_{j})$ such that $(x_{i},x_{j})$ is sampled with
probability at least $(1-\varepsilon)k(x_{i},x_{j})$.

######  Theorem 4.14 (Weighted Edge Sampling).

Algorithm 4.13 returns a random edge of $K$ with probability proportional to
at least $(1-\varepsilon)$ its weight using $1$ call to Algorithm 4.11.

###### Proof.

Consider an edge $(u,v)$. Vertex $u$ is sampled with probability at least
$(1-2\varepsilon)\frac{\textup{deg}(u)}{\sum_{x\in X}\textup{deg}(x)}$. Given
this, $v$ is then sampled with probability at least
$(1-2\varepsilon)\frac{k(u,v)}{\sum_{x\in
X\setminus{u}}k(u,x)}=(1-2\varepsilon)\frac{k(u,v)}{\textup{deg}(u)}.$ Using
the same analysis for sampling $v$ and then $u$, we have that any edge $(u,v)$
is sampled with probability at least $1-2\varepsilon$ times $k(u,v)/\sum_{e\in
K}w(e)$. Note that the same rejection sampling remark as in the proof of
Theorem 4.12 applies and we can perfectly sample an edge proportional to its
weight with an addition $O(1/\tau)$ rejection sampling steps. ∎

###  4.4 Random Walk

######  Theorem 4.15.

Algorithm 4.16 outputs a vertex from a vertex within $O(T\varepsilon)$ total
variation distance from the true random walk distribution. Each step of the
walk requires $1$ call to Algorithm 4.11 .

###### Proof.

The proof follows from the correctness of Algorithm 4.11 given in Theorem
4.12. Lastly we again note that by performing an additional $O(1/\tau)$ rounds
of rejection sampling steps (as outlined in the end of the proof of Theorem
4.12), we can make sure that we are sampling from the true random walk
distribution at each step of the walk. ∎

######  Algorithm 4.16 (Perform Random Walk).

Input: Dataset $X\subset\mathbbm{R}^{d}$, vertex $x_{i}\in X$, length of walk
$T\geqslant 1$. Operation: 1. Start at vertex $v\leftarrow x_{i}$. 2. For
$j=1$ to $T$: (a) Sample a random neighbor of $v$ using Algorithm 4.11. Let
$w$ be the resulting output. (b) Set $v\leftarrow w$. Output: Data point $v$.

##  5 Linear Algebra Applications

We now present a wide array of applications of the algorithmic building blocks
constructed in Section 4. Altogether, these applications allow us to
understand or approximate fundamental and properties of the kernel matrix and
the graph $K$. In this section we present the linear algebra applications and
the graph applications are given in Section 6.

###  5.1 Spectral Sparsification

######  Algorithm 5.1 (Spectral Sparsification of the Kernel Graph).

Input: Dataset $X\subset\mathbbm{R}^{d}$, accuracy parameter $\varepsilon$.
Operation: 1. Let $t=O(n\log(n)/\varepsilon^{2}\tau^{3})$ be the number of
edges that are to be sampled 2. Let $\hat{p}$ denote the distribution returned
by Algorithm 4.3 for a small enough constant $\varepsilon$. 3. Initialize
$G^{\prime}=\emptyset$. For $i=1,\ldots,t$: (a) Sample a vertex $u$ from the
distribution $\hat{p}$. (b) Sample a neighbor $v$ of $u$ using Algorithm 4.11
with constant $\varepsilon$. (c) Compute $\hat{q}_{uv}$, the probability that
Algorithm 4.11 samples $v$ given $u$ as input. (d) Similarly define and
compute $\hat{q}_{vu}$. Let
$w_{uv}=1/(t(\hat{p}_{u}\hat{q}_{uv}+\hat{p}_{v}\hat{q}_{vu}))$. (e) Add the
weighted edge $(\\{u,v\\},w_{uv})$ to the graph $G^{\prime}$. Output:

Given a set $X$, $|X|=n$, and a kernel $k:X\times
X\rightarrow\mathbbm{R}^{+}$, we describe how to construct a spectral
sparsifier for the weighted complete graph on $X$ where weight of the edge
$\\{x_{i},x_{j}\\}$ is given by $k(x_{i},x_{j})$.

######  Definition 5.2 (Graph Laplacian).

Given a weighted graph $G=(V,E,w)$, the Laplacian of $G$, denoted by
$L_{G}=D-A$, where $A$ is the adjacency matrix of $G$ with
$A_{i,j}=w(\\{i,j\\})$ and $D$ is a diagonal matrix such that for all
$i\in[n]$, $D_{i,i}=\sum_{j\neq i}A_{i,j}$.

######  Theorem 5.3 (Spectral Sparsification of Kernel Density Graphs).

Given a dataset $X$ of $n$ points in $\mathbb{R}^{d}$, and a kernel $k:X\times
X\to\mathbb{R}^{+}$, let $G=(X,{\binom{X}{2}},w)$ be the weighted complete
graph on $X$ with the weights $w(\\{x_{i},x_{j}\\})=k(x_{i},x_{j})$. Further,
for all $x_{i},x_{j}\in X$, let $k(x_{i},x_{j})\geqslant\tau$, for some
$\tau\in(0,1)$. Let $L_{G}$ be the Laplacian matrix corresponding to the graph
$G$. Then, for any $\varepsilon\in(0,1)$, Algorithm 5.1 outputs a graph
$G^{\prime}$ with only $m=O({n\log n}/{(\varepsilon^{2}\tau^{3})})$ edges,
such that with probability at least $9/10$,

$(1-\varepsilon)L_{G}\preceq L_{G^{\prime}}\preceq(1+\varepsilon)L_{G}.$

The algorithm makes $\widetilde{O}(m/\tau^{3})$ KDE queries and requires
$\tilde{O}(md/\tau^{3})$ post-processing time.

Let $G_{d}$ be the weighted directed graph obtained by arbitrarily orienting
the edges of the graph $G$ and let $H$ be an edge-vertex incidence matrix
defined as follows : for each $e=(x_{i},x_{j})$ in graph $G_{d}$, let
$H_{e,x_{i}}=\sqrt{k(x_{i},x_{j})}$ and $H_{e,x_{j}}=-\sqrt{k(x_{i},x_{j})}$.
Note that $H^{\top}H=L_{G}$. Our idea to construct spectral sparsifier is to
compute a sampling-and-reweighting matrix $S$, i.e., a matrix that has at most
one nonzero entry in each row, that with probability $\geqslant 9/10$,
satisfies

$(1-\varepsilon)L_{G}=(1-\varepsilon)H^{\top}H\preceq
H^{\top}S^{\top}SH\preceq(1+\varepsilon)H^{\top}H=(1+\varepsilon)L_{G}.$

The edges sampled by $S$ form the edges of the graph $G^{\prime}$. We
construct this matrix $S$ by sampling rows of the matrix $H$ from a
distribution close to the distribution that samples a row of $H$ with a
probability proportional to its squared norm. We show that this gives a
spectral sparsifier by showing that such a distribution approximates the
“leverage score sampling” distribution.

######  Definition 5.4 (Leverage Scores).

Let $M$ be a $n\times d$ matrix and $m_{i}$ denote the $i$-th row of $M$.
Then, for all $i\in[n]$, $\tau_{i}$, the $i$-th leverage of $M$ is defined as
follows:

$\tau_{i}=m_{i}(M^{\top}M)^{+}m_{i}^{\top},$

where $X^{+}$ is the Moore-Penrose pseudoinverse for a matrix $X$.

We introduce the following intermediate lemmas. We begin by recalling that
sampling edges proportional to leverage scores (effective resistances on a
graph) suffices to obtain spectral sparsification [SS11, Woo14].

######  Lemma 5.5 (Leverage Score Sampling implies Sparsification).

Given an $n\times d$ matrix $M$ and $\varepsilon\in(0,1)$, for all $i\in[t]$,
let $\tau_{i}$ be the $i$-th leverage score of $M$. Let
$p=\\{p_{1},p_{2},\ldots,p_{n}\\}$ be a distribution over the rows of $M$ such
that $p_{i}=\tau_{i}/\sum_{j\in[n]}\tau_{j}$. Further, for some
$\phi\in(0,1)$, let $\hat{p}=\\{\hat{p}_{1},\hat{p}_{2},\ldots,\hat{p}_{n}\\}$
be a distribution such that $\hat{p}_{i}\geqslant\phi p_{i}$ and let
$t=O\left(\frac{d\log(d)}{\varepsilon^{2}\phi}\right)$. Let
$S\in\mathbbm{R}^{t\times n}$ be a random matrix where for all $j\in[t]$, the
$j$-th row is independently chosen as $(1/\sqrt{t\hat{p}_{i}})e_{i}^{\top}$
with probability $\hat{p}_{i}$. Then, with probability at least $99/100$,

$(1-\varepsilon)M^{\top}M\preceq
M^{\top}S^{\top}SM\preceq(1+\varepsilon)M^{\top}M.$

Next, we show that the matrix $H$ is well-conditioned, in fact the condition
number is 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 $\varepsilon$ spectral sparsifier.

######  Lemma 5.6 (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 $1$ and at least $\sqrt{\tau}$. Let $\sigma_{\max}(H)$ be the maximum
singular value of $H$ and $\sigma_{\min}(H)$ be the minimum nonzero singular
value of $H$. Then $\sigma_{\max}(H)/\sigma_{\min}(H)\leqslant
4\sqrt{2}/\tau^{1.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$:

$\|A\|_{2}\leqslant\sqrt{\left(\max_{i}\sum_{j}|A_{i,j}|\right)\left(\max_{j}\sum_{i}|A_{i,j}|\right)}.$

For the matrix $H$, as each column has at most $n$ nonzero entries and each
row has at most $2$ non-zero entries and from the assumption that all the
entries have magnitude at most $1$, we obtain that
$\|H\|_{2}\leqslant\sqrt{2n}$. To obtain lower bounds on $\sigma_{\min}(H)$,
we appeal to a Cheeger-type inequality for weighted graphs from [Fri92, FN02].
First, we note that
$\sigma_{\min}(H)=\sqrt{\sigma_{\min}(H^{\top}H)}=\sqrt{\sigma_{\min}(L_{G})}$
where $G$ is the kernel graph that we are considering with each edge having a
weight of at least $\tau$. Let
$0=\lambda_{1}\leqslant\lambda_{2}\leqslant\cdots\leqslant\lambda_{n}$ be the
eigenvalues of the positive semi-definite matrix $L_{G}$. Now we have that

$\sigma_{\min}(L_{G})=\lambda_{2}(L_{G})\geqslant\min_{i}(\delta_{i}/2)\varepsilon(G)^{2}$

where $\delta_{i}=\sum_{j\neq i}k(x_{i},x_{j})$ i.e., the weighted degree of
vertex $x_{i}$ in graph $G$ and

$\varepsilon(G)=\min_{\phi\neq U\subset V,|U|\leqslant
n/2}\frac{|E(U,\bar{U})|}{|E(U)|}$

where $|E(U)|$ denotes the sum of weighted degrees of vertices in $U$ and
$|E(U,\bar{U})|$ denotes the total weight of edges with one end point in $U$
and the other outside $U$. Using the fact that $G$ is a complete graph with
each edge having a weight of at least $\tau$ and at most $1$, we obtain
$|E(U,\bar{U})|\geqslant\tau|U||\bar{U}|$ and $|E(U)|\leqslant n|U|$, which
implies that $\varepsilon(G)\geqslant\min_{\phi\neq U\subset V,|U|\leqslant
n/2}\tau|\bar{U}|/n\geqslant\tau/2$. We also similarly have that
$\min_{i}\delta_{i}\geqslant(n-1)\tau$, which overall implies that
$\lambda_{2}(L_{G})\geqslant n\tau^{3}/16$ and that
$\sigma_{\min}(H)\geqslant\sqrt{n}\tau^{1.5}/4$. Thus, we obtain that
$\sigma_{\max}(H)/\sigma_{\min}(H)\leqslant 4\sqrt{2}/\tau^{1.5}$. ∎

We are now ready to complete the proof of our main theorem:

###### Proof of Theorem 5.3.

Let $q=\\{q_{1},q_{2},\ldots,q_{{\binom{n}{2}}}\\}$ be a distribution over the
rows of $H$ such that for all edges $e=\\{i,j\\}$, $q_{e}\geqslant
c\frac{\left\lVert H_{e,*}\right\rVert_{2}^{2}}{\left\lVert
H\right\rVert_{F}^{2}}=\frac{k(x_{i},x_{j})}{\sum_{e^{\prime}=\\{i^{\prime},j^{\prime}\\}}k(x_{i^{\prime}},x_{j^{\prime}})}$,
for a fixed universal constant $c$.

Next, we show that this distribution is $\Theta(1/\kappa^{2})$ approximation
to the leverage score distribution for $H$. Let $H=U\Sigma V^{\top}$ be the
“thin” singular value decomposition of $H$ and therefore all the diagonal
entries of $\Sigma$ are nonzero. By definition $\tau_{i}=\|U_{i*}\|_{2}^{2}$.
We have

$\displaystyle\|h_{i}\|_{2}^{2}=\|U_{i*}\Sigma
V^{\top}\|_{2}^{2}=\|U_{i*}\Sigma\|_{2}^{2}$

where the equality follows from the fact that $V^{\top}$ has orthonormal rows.
Now, $\|U_{i*}\Sigma\|_{2}^{2}\geqslant\|U_{i*}\|_{2}^{2}\sigma_{\min}^{2}$
and $\|U_{i*}\Sigma\|_{2}^{2}\leqslant\|U_{i*}\|_{2}^{2}\sigma_{\max}^{2}$.
Therefore, for all $i\in{\binom{n}{2}}$, defining
$\kappa=\sigma_{\min}/\sigma_{\max}$, we have

$\frac{\tau_{i}}{\sum_{j}\tau_{j}}=\frac{\|U_{i*}\|_{2}^{2}}{\sum_{j}\|U_{j*}\|_{2}^{2}}\geqslant\frac{\|h_{i}\|_{2}^{2}/\sigma_{\max}^{2}}{\sum_{j}\|h_{j}\|_{2}^{2}/\sigma_{\min}^{2}}=\frac{1}{\kappa^{2}}\frac{\|h_{i}\|_{2}^{2}}{\|H\|_{F}^{2}}.$

Then, we invoke Lemma 5.5 with $\phi=\Omega(1/\kappa^{2})$ and conclude that
sampling $t=O\left(\frac{n\log n}{\varepsilon^{2}\kappa^{2}}\right)$ rows of
$H$ results in a sparse graph $G^{\prime}$ with corresponding Laplacian
$L_{G^{\prime}}$ such that with probability at least $99/100$,

$(1-\varepsilon/2)L_{G}\preceq L_{G^{\prime}}\preceq(1+\varepsilon/2)L_{G}.$

Further, by Lemma 5.6, we can conclude $\kappa^{2}\leqslant 32/\tau^{3}$ and
thus sampling $t=O\left(\frac{n\log n}{\varepsilon^{2}\tau^{3}}\right)$ edges
suffices.

We do not use Algorithm 4.13 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, we first use
Algorithm 4.6 with constant $\varepsilon$ (say 1/2) to sample a vertex $u$ and
Algorithm 4.11 with constant $\varepsilon$ (say 1/2) to sample a neighbor $v$
of $u$. Note that Algorithms 4.6 and Algorithms 4.11 can be modified to also
return the probabilities $\hat{p}_{u}$ and $\hat{q}_{vu}$ with which the
vertex $i$ and the neighbor $j$ of $i$ are sampled. We can further query the
algorithms to return $\hat{p}_{v}$ and $\hat{q}_{uv}$. Now,
$q_{\\{u,v\\}}=\hat{p}_{u}\hat{q}_{vu}+\hat{p}_{v}\hat{q}_{uv}$ is the
probability with which this sampling process samples the edge $\\{u,v\\}$ and
we have that $\hat{p}_{u}\hat{q}_{vu}+\hat{p}_{v}\hat{q}_{uv}\geqslant
c\frac{k(x_{u},x_{v})}{\sum_{i\neq j}k(x_{i},x_{j})}$ and we use this
distribution $q$ to implement spectral sparsification as described above. As
already seen (Theorem 4.12), to compute vertex sampling distribution
$\hat{p}$, we use $n$ KDE queries and for each neighbor sampling step, we use
$O(\log n)$ KDE queries. Thus, we overall use
$O(n\log^{2}n/(\varepsilon^{2}\tau^{3}))$ constant approximate KDE queries to
obtain an $\varepsilon$ spectral sparsifier. ∎

We can further compute another graph $G^{\prime\prime}$ with only
$O(n/\varepsilon^{2})$ edges by computing an $\varepsilon/2$ spectral
sparsifier for $G^{\prime}$ using the spectral sparsification algorithm of Lee
an Sun [LS18] (see Theorem 1.1). This procedure doesn’t require any KDE
queries and solely operates on the weighted graph $G^{\prime}$. The overall
running time is
$O\left(\frac{n^{1+1/c}\log(n)}{\epsilon^{6+1/c}\tau^{3}}\right)$, for a large
fixed constant $c$.

##### Hardness for spectral sparsification.

We observe that we can use the lower bound from Alman et. al. to establish
hardness in terms of $\tau$ from Parameterization 1.2. The lower bound we
obtain is as follows:

######  Theorem 5.7 (Lower Bound for Spectral Sparsification under
Parameterization 1.2).

Let $k$ be the Gaussian kernel and let $X$ be dataset such that $\min_{x,y\in
X}k(x,y)=\tau$, for some $1>\tau>0$. Then, any algorithm that with probability
$9/10$ outputs an $O(1)$-approximate spectral sparsifier for the kernel graph
associated with $X$, with $O(n^{2-\delta})$ edges, where $\delta<0.01$ is a
fixed universal constant, requires $\Omega\left(n\cdot
2^{\log(1/\tau)^{0.32}}\right)$ time, assuming the strong exponential time
hypothesis.

First, we begin with the definition of a multiplicatively-Lipschitz function:

######  Definition 5.8 (Multiplicatively-Lipschitz Kernels).

A kernel $k$ over a set $X$ is $(c,L)$-multiplicatively Lipschitz if for any
$\rho\in(1/c,c)$, and for any $x,y\in X$, $c^{-L}k(x,y)\leqslant k(\rho x,\rho
y)\leqslant c^{L}k(x,y)$.

We will require the following theorem showing hardness for sparsification when
the kernel function is not multiplicatively-Lipschitz:

######  Theorem 5.9 (Theorem 8.3 [ACSS20]).

Let $k$ be a function and $X$ be a dataset such that $k$ is not
$(c,L)$-multiplicatively-Lipschitz on $X$ for some $L>1$ and
$c=1+2\log\left(10\cdot 2^{L^{0.48}}\right)/L$. Then, there is no algorithm
that returns a sparsifier of the kernel graph associated with $X$ with
$O(n^{2-\delta})$ edges, where $\delta<0.01$ is a fixed universal constant, in
less than $O\left(n\cdot 2^{L^{0.48}}\right)$ time, assuming the strong
exponential time hypothesis.

###### Proof of Theorem 5.7 .

First, we show that for any $c>1$, if $L<\log(1/\tau)(c-1)$, then the Gaussian
kernel $k$ is not $(c,L)$-multiplicatively Lipschitz. Let $z=\left\lVert
x-y\right\rVert_{2}^{2}$ and let $f(z)=e^{-z}$. Observe, it suffices to show
that there exists a $z$ such that $f(cz)\leqslant c^{-L}f(z)$. Let $z$ be such
that $f(z)=e^{z}=\min_{x,y}k(x,y)=\tau$, i.e. $z=\log(1/\tau)$. Then,

$f(c\log(1/\tau))=e^{-c\log(1/\tau)},$

and for $L<\log(1/\tau)(c-1)$

$c^{-L}f(\log(1/\tau))>e^{-c\log(1/\tau)}.$

Then, applying Theorem 5.9 with $c=1+\frac{1}{\sqrt{L}}$, it suffices to
conclude $k$ is not $(c,L)$-multiplicatively Lipschitz when
$L<\log^{2/3}(1/\tau)$, which concludes the proof. ∎

####  5.1.1 Solving Laplacian Systems Approximately

We describe how to approximately solve the Laplacian system $L_{G}x=b$ using
the spectral sparsifier $L_{G^{\prime}}$. First, we note the following theorem
that states the running time and approximation guarantees of fast Laplacian
solvers.

######  Theorem 5.10 ([KMP11], [ST04]).

There is an algorithm that takes an input a graph Laplacian $L$ of a graph
with $m$ weighted edges, a vector $b$, and an error parameter $\alpha$ and
returns $x$ such that with probability at least $99/100$,

$\|x-L^{+}b\|_{L}\leqslant\alpha\|L^{+}b\|_{L},$

where $\|x\|_{L}=\sqrt{x^{\top}Lx}$. The algorithm runs in time
$\widetilde{O}(m\log(1/\alpha))$.

We have the following theorem that bounds the difference between solutions for
the exact Laplacian system and the spectral sparsifier Laplacian.

######  Theorem 5.11.

Let $L_{G}$ be the Laplacian of a connected graph $G$ on $n$ vertices and let
$L_{G^{\prime}}$ be the Laplacian of an $\varepsilon$-spectral sparsifier
$G^{\prime}$ of graph $G$ i.e.,

$\displaystyle(1-\varepsilon)L_{G}\preceq
L_{G^{\prime}}\preceq(1+\varepsilon)L_{G},$

for $\varepsilon<c$ for a small enough constant $c$. Then, for any vector $b$
with $1^{\top}b=0$, $\|L_{G}^{+}b-L_{G^{\prime}}^{+}b\|_{L_{G}}\leqslant
2\sqrt{\varepsilon}\|L_{G}^{+}b\|_{L_{G}}$.

###### Proof.

Note that for $\varepsilon<1$, the graph $G^{\prime}$ also has to be connected
and therefore the only eigen vectors corresponding to eigen value $0$ of the
matrices $L_{G}$ and $L_{G^{\prime}}$ are of the form $a\cdot 1$ for $a\neq 0$
and hence columns (and rows) of $L_{G}$ span all vectors orthogonal to $1$.
Therefore $L_{G}L_{G}^{+}=L_{G}^{+}L_{G}=I-(1/n)11^{\top}$. Now,

$\displaystyle\|L_{G}^{+}b-L_{G^{\prime}}^{+}b\|_{L_{G}}^{2}$
$\displaystyle=b^{\top}(L_{G}^{+}-L_{G^{\prime}}^{+})L_{G}(L_{G}^{+}-L_{G^{\prime}}^{+})b$
$\displaystyle=b^{\top}L_{G}^{+}L_{G}L_{G}^{+}b-b^{\top}L_{G^{\prime}}^{+}L_{G}L_{G}^{+}b-b^{\top}L_{G^{\prime}}^{+}L_{G}L_{G}^{+}b+b^{\top}L_{G^{\prime}}^{+}L_{G}L_{G^{\prime}}^{+}b$
$\displaystyle\leqslant
b^{\top}L_{G}^{+}b-b^{\top}L_{G^{\prime}}^{+}b-b^{\top}L_{G^{\prime}}^{+}b+\frac{1}{1-\varepsilon}b^{\top}L_{G^{\prime}}^{+}b$

where in the last inequality, we used $L_{G}L_{G}^{+}b=1$ and that for any
vector $v$,
$v^{\top}L_{G}v\leqslant\frac{1}{1-\varepsilon}v^{\top}L_{G^{\prime}}v$. As
the null spaces of both $L_{G}$ and $L_{G^{\prime}}$ are given by
$\\{a1\,|\,a\in\mathbb{R}\\}$, we also obtain that

$\displaystyle(1-\varepsilon)L_{G}^{+}\preceq
L_{G^{\prime}}^{+}\preceq(1+\varepsilon)L_{G}^{+}$

using which we further obtain that

$\displaystyle\|L_{G}^{+}b-L_{G^{\prime}}^{+}b\|_{L_{G}}^{2}$
$\displaystyle\leqslant\left(\frac{2}{1-\varepsilon}-2\right)b^{\top}L_{G^{\prime}}^{+}b\leqslant\frac{2\varepsilon(1+\varepsilon)}{1-\varepsilon}b^{\top}L_{G}^{+}b\leqslant
4\varepsilon\|L_{G}^{+}b\|_{L_{G}}^{2}.$

Thus, $\|L_{G}^{+}b-L_{G^{\prime}}^{+}b\|_{L_{G}}\leqslant
2\sqrt{\varepsilon}\|L_{G}^{+}b\|_{L_{G}}$. ∎

Therefore, if $x$ is a vector such that
$\|x-L_{G^{\prime}}b\|_{L_{G^{\prime}}}\leqslant\alpha\|L_{G^{\prime}}^{+}b\|_{L_{G^{\prime}}}$
obtained using the fast Laplacian solver, then

$\displaystyle\|x-L_{G}^{+}b\|_{L_{G}}^{2}$
$\displaystyle=\|x-L_{G^{\prime}}^{+}b+L_{G^{\prime}}^{+}b-L_{G}^{+}b\|_{L_{G}}^{2}$
$\displaystyle\leqslant
2(\|x-L_{G^{\prime}}^{+}b\|_{L_{G}}^{2}+\|L_{G^{\prime}}^{+}b-L_{G}^{+}b\|_{L_{G}}^{2})$
$\displaystyle\leqslant\frac{2}{1-\varepsilon}\|x-L_{G^{\prime}}^{+}b\|_{L_{G}^{\prime}}^{2}+4\varepsilon\|L_{G}^{+}b\|_{L_{G}}^{2}.$

Here we used the above theorem and the fact that
$L_{G}\preceq(1/(1-\varepsilon))L_{G^{\prime}}$. Now,
$\|x-L_{G^{\prime}}^{+}b\|_{L_{G^{\prime}}}^{2}\leqslant\alpha^{2}\|L_{G^{\prime}}^{+}b\|_{L_{G^{\prime}}}^{2}$
and
$\|L_{G^{\prime}}^{+}b\|_{L_{G^{\prime}}}^{2}=b^{\top}L_{G^{\prime}}^{+}L_{G^{\prime}}L_{G^{\prime}}^{+}b=b^{\top}L_{G^{\prime}}^{+}b\leqslant(1+\varepsilon)b^{\top}L_{G}^{+}b\leqslant(1+\varepsilon)\|L_{G}^{+}b\|_{L_{G}}^{2}$,
which finally implies that

$\displaystyle\|x-L_{G}^{+}b\|_{L_{G}}^{2}\leqslant\left(\frac{2(1+\varepsilon)^{2}}{1-\varepsilon}\alpha^{2}+4\varepsilon\right)\|{L_{G}^{+}b}\|_{L_{G}}^{2}.$

Thus, using a $\varepsilon$ spectral sparsifier $G^{\prime}$ with $m$ edges,
we can in time $\widetilde{O}(m\log(1/\varepsilon))$ can obtain a vector $x$
such that $\|x-L_{G}^{+}b\|_{L_{G}}\leqslant
C\sqrt{\varepsilon}\|L_{G}^{+}b\|_{L_{G}}$ for a large enough constant $C$.

###  5.2 Low-rank Approximation of the Kernel Matrix

We derive algorithms for low-rank approximations of the kernel matrix via KDE
queries. We present a algorithm for additive error approximation and compare
to prior work for relative error approximation.

We first recall the following two theorems. Let $A_{i,*}$ denote the $i$th row
of a matrix $A$.

######  Theorem 5.12 ([FKV04]).

Let $A\in\mathbbm{R}^{n\times m}$ be any matrix. Let $S$ be a sample of
$O(k/\varepsilon)$ rows according to a probability distribution
$(p_{1},\ldots,p_{n})$ that satisfies
$p_{i}\geqslant\Omega(1)\cdot\|A_{i,*}\|_{2}^{2}/\|A|_{F}^{2}$ for every
$1\leqslant i\leqslant n$. Then, in time
$O(mk/\varepsilon\cdot\mathrm{poly}(k,1/\varepsilon))$, we can compute from
$S$ a matrix $U\in\mathbbm{R}^{k\times m}$, that with probability at least
$0.99$ satisfies

$\|A-AU^{T}U\|_{F}^{2}\leqslant\|A-A_{k}\|_{F}^{2}+\varepsilon\|A\|_{F}^{2}.$

######  Theorem 5.13 ([CP17], also see [IVWW19]).

There is a randomized algorithm that given matrices $A\in\mathbbm{R}^{n\times
m}$ and $U\in\mathbbm{R}^{k\times m}$, reads only $O(k/\varepsilon)$ columns
of $A$, runs in time $O(mk)+\mathrm{poly}(k,1/\varepsilon)$, and returns
$V\in\mathbbm{R}^{n\times k}$ that with probability $0.99$ satisfies

$\|A-VU\|_{F}^{2}\leqslant(1+\varepsilon)\min_{X\in\mathbbm{R}^{n\times
k}}\|A-X\|_{F}^{2}.$

Therefore to compute the low rank approximation, we just need sample from the
distribution on rows required by Theorem 5.12. 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
$p_{i}$ in the distribution $(p_{1},\ldots,p_{n})$ is the sum of edge weights
squared for vertex $x_{i}$. From vertex queries (Algorithm 4.6), we know that
we can get the degree of each vertex, which is the sum of edge weights. We can
extend Algorithm 4.6 to sample from the sum of squared edge weights of each
vertex as follows. Consider a kernel $k$ such that there exists an absolute
constant $c$ that satisfies $k(x,y)^{2}=k(cx,cy)$ for all $x,y$. Such a $c$
exists for the most popular kernels such as the Laplacian, exponential, and
Gaussian kernels for which $c=2,2,$ and $4$ respectively. Thus give our
dataset $X$, we simply construct KDE queries for the dataset $X^{\prime}:=cX$.
Then by sampling the degrees of the vertices associated with the kernel graph
$K^{\prime}$ of $X^{\prime}$, we can sample from the distribution required by
Theorem 5.12 by invoking Algorithm 4.6 on the dataset $X^{\prime}$. In
particular, using $n$ KDE queries for $X^{\prime}$, we can get row norm
squared values for all rows of our original kernel matrix $K$. We can then
sample the rows according to Theorem 5.12 and fully construct the rows that
are sampled. Altogether, this takes $n$ KDE queries and $O(nk/\varepsilon)$
kernel function evaluations to construct a rank $k$ approximation of $K$; see
Algorithm 5.15.

######  Corollary 5.14.

Given a dataset $X$ of size $n$, there exists an algorithm that outputs a rank
$k$ matrix $B$ such that

$\|K-B\|_{F}^{2}\leqslant\|K-K_{k}\|_{F}^{2}+\varepsilon\|K\|_{F}^{2}$

with probability $99/100$, where $K$ is a kernel matrix associated with $X$
based on a Laplacian, exponential, or Gaussian kernel, and $K_{k}$ is the
optimal rank-$k$ approximation of $K$. It uses $n$ KDE queries and
$O(nk/\varepsilon\cdot\mathrm{poly}(k,1/\varepsilon)+nkd/\varepsilon)$ post-
processing time.

We remark that for the application presented in this subsection, we can we can
replace 1.2. Indeed, since we only estimate row sums, we only require that the
value of a KDE query is at least $\tau$, that is, the average value
$\frac{1}{|X|}\sum_{x\in X}k(x,y)\geqslant\tau$ for a query $y$. Note that via
Cauchy Schwartz, this automatically implies a lower bound for the average
squared sum:

$\frac{1}{|X|}\sum_{x\in
X}k(x,y)^{2}\geqslant\frac{1}{|X|^{2}}\left(\sum_{x\in
X}k(x,y)\right)^{2}\geqslant\tau^{2}.$

######  Algorithm 5.15 (Additive-error Low-rank Approximation).

Input: Kernel matrix $K\in\mathbbm{R}^{n\times n}$, data points $X\subset
R^{d}$, accuracy parameter $\varepsilon$, rank parameter $k$. Operation: 1.
Let $c$ be the constant such that $k(x,y)^{2}=k(cx,cy)$ for all inputs $x,y$.
For $i=1$ to $i=n$: (a) Compute the value
$p_{i}=\sum_{j=1}^{n}k(cx_{i},cx_{j})$ using KDE queries for the dataset $cX$.
2. Sample and construct $O(k/\varepsilon)$ rows of $K$ according to
probability proportional to $\\{p_{i}\\}_{i=1}^{n}$. 3. Compute $U$ from the
sample, using Theorem 5.12. 4. Compute $V$ from the sample, using Theorem
5.13. Output: Factors $U,V$ such that $\left\lVert
K-UV\right\rVert_{F}^{2}\leqslant\left\lVert
K-K_{k}\right\rVert_{F}^{2}+\varepsilon\left\lVert K\right\rVert_{F}^{2}$

###  5.3 Approximating the Spectrum in EMD

In this subsection, we obtain a sublinear time algorithm to approximate the
spectrum of the normalized 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
[LGT12, LRTV12, KLL+13, CPS15, GKL+21], which show that the $j$th 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 a
succinct summary of all eigenvalues, i.e. the spectrum.

Additionally, models of 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 an example from [CKSV18], “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 a 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 5.16 (Earth Mover Distance).

Given two multi-sets of $n$ points in $\mathbb{R}^{d}$, denoted by $A$ 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.

$\textrm{EMD}(A,B)=\min_{\pi:A\to B}\sum_{a\in A}\left\lVert
a-\pi(a)\right\rVert_{2},$ (2)

where $\pi$ ranges over all one-to-one mappings.

We can now invoke the algorithm ApproxSpectralMoment of [CKSV18]. The
algorithm first selects uniformly random vertices of a weighted graph $A$. It
then performs a random walk of a specified length $\ell$ starting from the
chosen vertex and then counts the number of times the walk returns back to the
original vertex. Now Theorem 4.15 allows us to perform one step of a random
walk using $O(\log n)$ KDE queries. Note that we perform an additional
$\tilde{O}(1/\tau)$ of rejection sampling in Algorithm 4.11 to perfectly
sample from the true neighbor distribution. Thus we immediately have the
following guarantee:

######  Theorem 5.17 (Corollary of Theorem 1 in [CKSV18] and Theorem 4.15).

Given a $n\times n$ kernel matrix $K$ and accuracy parameter
$\varepsilon\in(0,1)$, let $G$ be the corresponding weighted graph, and let
$L_{G}=I-D^{-1}KD^{-1}$ be the normalized Laplacian, where
$D_{i,i}=\sum_{j}K_{i,j}$. Let
$\lambda_{1}\geqslant\lambda_{2}\ldots\geqslant\lambda_{n}$ be the eigenvalues
of $L_{G}$ and let $\lambda$ be the resulting vector. Then, there exists an
algorithm that uses
$\widetilde{O}\left(\exp\left(1/\varepsilon^{2}\right)/\tau\right)$ KDE
queries and $\exp\left(1/\varepsilon^{2}\right)\cdot d/\tau$ post-processing
time and outputs a vector $\widetilde{\lambda}$ such that with probability
$99/100$,

$\textrm{EMD}\left(\lambda,\widetilde{\lambda}\right)\leqslant\varepsilon.$

We remark that the bound of $\exp\left(1/\varepsilon^{2}\right)$ is
independent of $n$, which is the size of the dataset.

###  5.4 First Eigenvalue and Eigenvector Approximation

Our goal is to approximate the top eigenvalue of the kernel matrix and find a
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 $n\times n$ kernel matrix $K$
satisfies that it’s sum is at least $n\tau$, 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 a ‘noisy’ power method of [BIMW21] on the
sampled submatrix. We note that this step employs a KDE data-structure
initialized only on the sampled indices of $K$. The algorithm and details
follow.

######  Algorithm 5.18 (First Eigenvalue and Eigenvector Approximation).

Input: Input dataset $X\subset\mathbbm{R}^{d}$ of size $|X|=n$, precision
$\varepsilon>0$. Operation: 1. Let $t\leftarrow
O(1/(\varepsilon^{2}\tau^{2})$. Let $S\leftarrow$ random subset of $[n]$ of
size $t$. Let $X_{S}$ be the samples restricted to the indices in $S$. 2. Let
$K_{S}\leftarrow$ principal submatrix of $K$ on indices in $S$ and let
$\tilde{K}\leftarrow(n/s)\cdot K_{S}$. // Just for notation; we do not
initialize $K$ or $K_{S}$ 3. Construct a KDE data structure for $X_{S}$. Run
Algorithm $1$ of [BIMW21] (Kernel Noisy Power Method) on $K_{S}$. Let
$\hat{\lambda}_{\max}$ be the resulting eigenvalue and $\hat{v}_{\max}$ be the
resulting eigenvector. Output: $\hat{\lambda}_{\max}$ and $\hat{v}_{\max}$.

We remark that the eigenvector returned by Algorithm 5.18 will be a sparse
vector supported only on the coordinates in $S$.

We first state the necessary auxiliary statements needed to prove the
guarantees of Algorithm 5.18.

######  Lemma 5.19.

If each row of $K$ satisfies that its sum is at least $n\tau$ for parameter
$\tau\in(0,1)$, then the largest eigenvalue of $K$, denoted as $\lambda_{1}$,
satisfies $\lambda_{1}\geqslant n\tau$.

###### Proof.

This follows from looking at the quadratic form $\textbf{1}^{T}K\textbf{1}$
where 1 is the vector with all entries equal to $1$:

$\lambda_{1}\geqslant\frac{\textbf{1}^{T}K\textbf{1}}{\textbf{1}^{T}\textbf{1}}\geqslant\frac{n^{2}\tau}{n}=n\tau.\qed$

We now state the guarantees of Algorithm $1$ in [BIMW21].

######  Theorem 5.20 ([BIMW21]).

Suppose the kernel function for a $m\times m$ kernel matrix $K$ has a KDE data
structure with query time $d/(\varepsilon^{2}\tau^{p})$ (see Table 1). Then
Algorithm $1$ of [BIMW21] returns $\lambda$ such that
$\lambda\geqslant(1-\varepsilon)\lambda_{1}(K)$ in time
$O\left(\frac{dm^{1+p}\log(m/\varepsilon)^{2+p}}{\varepsilon^{7+4p}}\right)$

Finally, we need the following result on eigenvalues of sampled PSD matrices,
proven in [BMR21].

######  Lemma 5.21 ([BMR21]).

Let $A\in\mathbbm{R}^{n\times n}$ be PSD with $\left\lVert
A\right\rVert_{\infty}\leqslant 1$. Let $S\subset[n]$ be a random subset of
size $t$ and let $A_{S\times S}$ be the submatrix restricted to columns and
rows in $S$ and scaled by $n/s$. Then, for all $i\in[|S|]$,
$\lambda_{i}\left(A_{S\times S}\right)=\lambda_{i}(A)\pm\frac{n}{\sqrt{t}}$.

We are now ready to prove the guarantees of Algorithm 5.18.

######  Theorem 5.22.

Given a $n\times n$ kernel matrix $K$ admitting a KDE data-structure with
query time $d/(\varepsilon^{2}\tau^{p})$, Algorithm 5.18 returns $\lambda$
such that $\lambda\geqslant(1-\varepsilon)\lambda_{1}(K)$ in total time

$\min\left(O\left(\frac{d\log(d/\varepsilon)}{\varepsilon^{4.5}\tau^{4}}\right),O\left(\frac{d}{\varepsilon^{9+6p}\tau^{2+2p}}\log\left(\frac{1}{\varepsilon\tau}\right)^{2+p}\right)\right).$

######  Remark 5.23.

Two remarks are in order. First we recall that the runtime of [BIMW21] has a
$n^{1+p}$ factor while our bound has no dependence on $n$ and is thus a truly
sublinear runtime. Second, if we skip the Kernel Noisy Power method step and
directly initialize and calculate the top eigenvalue of $K_{S}$ (using the
standard gap independent power method of [MM15]), we would get a runtime of
$\tilde{O}(d/(\varepsilon^{4.5}\tau^{4}))$ which has a polynomially better
$\varepsilon$ dependence but a worse $\tau$ dependence than the guarantees of
Algorithm 5.18.

###### Proof of Theorem 5.22.

We first prove the approximation guarantee. By our setting of $t$ and using
Lemma 5.21, we see that the additive error in approximating the first
eigenvalue of $K$ by that of $\tilde{K}$ is at most

$\frac{n}{\sqrt{t}}\leqslant\varepsilon\tau
n\leqslant\varepsilon\lambda_{1}(K),$

and thus $\lambda_{1}(\tilde{K})\geqslant(1-\varepsilon)\lambda_{1}(K)$. Then
by the guarantees of Theorem 5.20, it follows that we find a $1-\varepsilon$
multiplicative approximation to $\lambda_{1}(\tilde{K})$ and thus a
$1-O(\varepsilon)$ multiplicative approximation to that of $\lambda_{1}(K)$.

We now prove the runtime bound. It easily follows from plugging in
$m=O(1/(\varepsilon^{2}\tau^{2}))$ in Theorem 5.20. ∎

##  6 Graph Applications

In this section, we present our graph applications, including local
clustering, spectral clustering, arboricity estimation, and estimating the
total weight of triangles.

###  6.1 Local Clustering

######  Algorithm 6.1 (Local $k$-Clustering).

Input: Input dataset $X\subset\mathbbm{R}^{d}$ of size $|X|=n$, vertices
$u,w$, random walk length $t$. Operation: 1. For a given $v$, let $p_{v}^{t}$
be the endpoint distribution of a random walk of length $t$ starting at $v$.
Output: “$u,w$ are in the same cluster” if $\ell_{2}$ distribution tester (see
Theorem 6.5) outputs $\|p_{u}^{t}-p_{w}^{t}\|_{2}\leqslant 1/(7n)$. Otherwise,
output “$u,w$ are in different clusters”.

We give a local clustering algorithm on graphs. The advantage of this method
is that it is local as it 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 at 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 definition in the property testing and sublinear algorithms
community (see [KS08, CS10, GR11, CPS15, CKK+18, DPRS19, GKL+21] and the
references within).

First, we need to define the notion of conductance.

######  Definition 6.2 (Conductance).

Let $G=(V,E,w)$ be a weighted graph. The conductance of a set $S\subset V$ is
defined as

$\phi_{G}(S)=\frac{w(S,S^{c})}{\min(w(S),w(S^{c}))}$

where $w(S,S^{c})$ denotes the sum of edge weights crossing the cut
$(S,S^{c})$ and $w(S)$ denotes the sum of (weighted) degrees of vertices in
$S$. The conductance of the graph $G$ is then the minimum of $\phi_{G}(S)$
over all sets $S$:

$\phi(G)=\min_{S}\phi_{G}(S).$

######  Definition 6.3 (Inner/Outer Conductance).

For a subset $U\subseteq V$, we define $\phi(G[U])$ to be the conductance of
the induced graph on $U$. $\phi(G[U])$ is also referred to as the inner
conductance of $U$. Conversely, $\phi_{G}(U)$ is refereed to as the outer
conductance of $U$.

######  Definition 6.4 ($k$-clusterable Graph).

A graph $G$ is $(k,\phi_{in},\phi_{out})$-clusterable if the following holds:
There exists a partition of the vertex set into $h\leqslant k$ parts
$V=\cup_{1\leqslant i\leqslant h}V_{i}$ such that
$\phi(G[V_{i}])\geqslant\phi_{in}$ and $\phi_{G}(V_{i})\leqslant\phi_{out}$.

Definition 6.4 captures the intuition that one can partition the graph into
$h\leqslant k$ pieces where each piece has a strong cluster structure
(captured by $\phi_{in}$) and distinct pieces are separated by sparse cuts
(captured by $\phi_{out}$). Note that we are interested in the regime where
$\phi_{out}$ is smaller than $\phi_{in}$. We will also assume that each
$|V_{i}|\geqslant n/\mathrm{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 6.4 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 $u$ and $w$ that are in parts $V_{1}$ and
$V_{2}$ respectively of the graph (as defined in Definition 6.4), we can
efficiently test if $V_{1}=V_{2}$ or $V_{1}\neq V_{2}$. That is, we can
efficiently test if $u$ and $w$ belong to the same or distinct clusters. The
underlying idea behind the algorithm is that if $u$ 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 detected using distribution
testing. Our algorithm is given in Algorithm 6.1. The flavor of the algorithm
presented is quite standard in property testing literature, see [CPS15] and
[Pen20].

The $\ell_{2}$ distribution tester we need is a standard result in
distribution testing with the following guarantees.

######  Theorem 6.5 (Theorem $1.2$ in [CDVV14]).

Let $\delta,\xi>0$ and let $p,q$ be two discrete distributions over a set of
size $n$ with $b\geqslant\max\\{\|p\|_{2}^{2},\|q\|_{2}^{2}\\}.$ Let
$r\geqslant c\sqrt{b}\log(1/\delta)/\xi$ for an appropriate constant $c$.
There exists $\ell_{2}$ distribution tester that takes as input $r$ samples
from each distribution $p,q$ and accepts the distributions if
$\|p-q\|_{2}^{2}\leqslant\xi$, and rejects the distributions if
$\|p-q\|_{2}^{2}\geqslant 4\xi$ with probability at least $1-\delta$. The
running time of the tester is linear in its sample size.

We now prove the correctness of Algorithm 6.1. We note that many arguments
from prior works are re-derived in the proof below, rather than stating them
in a 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 is given by $M=AD^{-1}$ where $A$ is the adjacency
matrix and $D$ is the diagonal degree matrix. The normalized Laplacian matrix
$L$ is defined as $L=I-D^{-1/2}AD^{-1/2}$.

Our first result is that vertices in the same well connected cluster of $G$
have a quantitative relationship captured by the eigenvectors of $L$. This is
in similar spirit to Lemma $5.3$ of [CPS15] but we must show it holds for
weighted graphs arising from kernel matrices whereas [CPS15] is interested in
bounded degree unweighted graphs.

######  Lemma 6.6.

Let $\textbf{v}_{i}$ be the $i$th eigenvector of the normalized Laplacian of
the kernel graph $K$ and let $C$ be any subset such that
$\phi(K[C])\geqslant\phi_{in}$. Then for any $1\leqslant i\leqslant h$, the
following holds:

$\sum_{u,v\in
C}\left(\frac{\textbf{v}_{i}(u)}{\sqrt{w(u)}}-\frac{\textbf{v}_{i}(v)}{\sqrt{w(v)}}\right)^{2}\lesssim\frac{\phi_{out}n}{\phi_{in}^{2}|C|\tau^{2}}.$

###### Proof.

By Lemma $5.2$ in [CPS15] and Theorem $1.2$ in [LGT12], we have that
$\phi_{in}^{2}/h^{4}\lesssim\lambda_{h+1}$ and $\lambda_{i}\leqslant
2\phi_{out}$ for any $1\leqslant i\leqslant h$. Now by the variational
principle for eigenvalues [CG97], we have

$\lambda_{i}=\sum_{(u,v)}\left(\frac{\textbf{v}_{i}(u)}{\sqrt{w(u)}}-\frac{\textbf{v}_{i}(v)}{\sqrt{w(v)}}\right)^{2}w(u,v)\leqslant
2\phi_{out}.$

Now let $H=K[C]$. From [CG97] and our assumptions on $C$, we have that

$\text{vol}_{H}(V_{H})\cdot\frac{2\cdot\sum_{(u,v)\in
E_{H}}\left(\frac{\textbf{v}_{i}(u)}{\sqrt{w(u)}}-\frac{\textbf{v}_{i}(v)}{\sqrt{w(v)}}\right)^{2}w(u,v)}{\sum_{u,v\in
V_{H}}\left(\frac{\textbf{v}_{i}(u)}{\sqrt{w(u)}}-\frac{\textbf{v}_{i}(v)}{\sqrt{w(v)}}\right)^{2}d_{H}(u)d_{H}(v)}\geqslant\lambda_{2}(H)\geqslant\frac{\phi_{in}^{2}}{2},$

where $\text{vol}_{H}(V_{H})$ denotes the sum of the degrees of vertices in
$H$ and $d_{H}(\cdot)$ 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

$\sum_{(u,v)\in
E_{H}}\left(\frac{\textbf{v}_{i}(u)}{\sqrt{w(u)}}-\frac{\textbf{v}_{i}(v)}{\sqrt{w(v)}}\right)^{2}\leqslant\sum_{(u,v)\in
E_{K}}\left(\frac{\textbf{v}_{i}(u)}{\sqrt{w(u)}}-\frac{\textbf{v}_{i}(v)}{\sqrt{w(v)}}\right)^{2}\leqslant
2\phi_{out}.$

This implies that

$\displaystyle|C|^{2}\tau^{2}\sum_{(u,v)\in
V_{H}}\left(\frac{\textbf{v}_{i}(u)}{\sqrt{w(u)}}-\frac{\textbf{v}_{i}(v)}{\sqrt{w(v)}}\right)^{2}$
$\displaystyle\leqslant\sum_{(u,v)\in
V_{H}}\left(\frac{\textbf{v}_{i}(u)}{\sqrt{w(u)}}-\frac{\textbf{v}_{i}(v)}{\sqrt{w(v)}}\right)^{2}d_{H}(u)d_{H}(v)$
$\displaystyle\lesssim\frac{\phi_{out}\text{vol}_{H}(V_{H})}{\phi_{in}^{2}}$

where we have used the fact that all edge weights in $K$ are at least $\tau$.
Using the fact that $\text{vol}_{H}(V_{H})\leqslant|C|n$, it follows that

$\sum_{(u,v)\in
V_{H}}\left(\frac{\textbf{v}_{i}(u)}{\sqrt{w(u)}}-\frac{\textbf{v}_{i}(v)}{\sqrt{w(v)}}\right)^{2}\lesssim\frac{\phi_{out}n}{\phi_{in}^{2}|C|\tau^{2}},$

as desired. ∎

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 [CPS15] but we must show it holds for weighted graphs.

######  Lemma 6.7.

Let $0<\beta<1/2$. If graph $K$ is $(k,\phi_{in},\phi_{out})$-clusterable, and
$C\subseteq V$ is any subset such that $|C|\geqslant n/\mathrm{poly}(k)$ and
$\phi(K[C])\geqslant\phi_{in}$. There exists a constant $c=c(\beta)>0$ and
$c^{\prime}=c^{\prime}(\beta,k)$ such that for any $t\geqslant c\log
n/\phi_{in}^{2},\phi_{out}\leqslant c^{\prime}\phi_{in}^{2}$, there exists a
subset $\widetilde{C}\subseteq C$ satisfying
$\text{vol}(\widetilde{C})\geqslant(1-\beta)\text{vol}(C)$ such that for any
$u,v\in\widetilde{C}$, the following holds:

$\|p_{u}^{t}-p_{v}^{t}\|_{2}^{2}\leqslant\frac{1}{8n}.$

###### Proof.

Let $\textbf{v}_{1},\cdots,\textbf{v}_{n}$ denote the eigenvectors of $L$ with
eigenvalues $\lambda_{1},\cdots,\lambda_{n}$ in non-decreasing order. We know
that the eigenvalues of $M$ are given by $1-\lambda_{i}$ with corresponding
eigenvalues $\textbf{y}_{i}=D^{1/2}\textbf{v}_{i}$. The vector $p_{u}^{t}$ is
the vector $1_{u}$ with a one value in the $u$th coordinate applied to
$M^{t}$. Write

$1_{u}=\sum_{i}\alpha_{i}y_{i}=\sum_{i}\alpha_{i}D^{-1/2}\textbf{v}_{i}.$

Taking the innerproduct of $1_{u}$ with $D^{-1/2}\textbf{v}_{i}$ tells us that
$\alpha_{i}=\textbf{v}_{i}(u)/\sqrt{w(u)}$. Thus,

$\displaystyle p_{u}^{t}-p_{v}^{t}$
$\displaystyle=\sum_{i=1}^{n}\frac{\textbf{v}_{i}(u)}{\sqrt{w(u)}}(1-\lambda_{i})^{t}\textbf{y}_{i}-\sum_{i=1}^{n}\frac{\textbf{v}_{i}(v)}{\sqrt{w(v)}}(1-\lambda_{i})^{t}\textbf{y}_{i}$
$\displaystyle=D^{1/2}\sum_{i=1}^{n}\textbf{v}_{i}\left(\frac{\textbf{v}_{i}(u)}{\sqrt{w(u)}}-\frac{\textbf{v}_{i}(v)}{\sqrt{w(v)}}\right)(1-\lambda_{i})^{t}.$

This means that

$\|p_{u}^{t}-p_{v}^{t}\|_{2}\leqslant\|D^{1/2}\|\left\|\sum_{i=1}^{n}\textbf{v}_{i}\left(\frac{\textbf{v}_{i}(u)}{\sqrt{w(u)}}-\frac{\textbf{v}_{i}(v)}{\sqrt{w(v)}}\right)(1-\lambda_{i})^{t}\right\|_{2}.$

Since the $\textbf{v}_{i}$’s are orthogonal, we know that

$\left\|\sum_{i=1}^{n}\textbf{v}_{i}\left(\frac{\textbf{v}_{i}(u)}{\sqrt{w(u)}}-\frac{\textbf{v}_{i}(v)}{\sqrt{w(v)}}\right)(1-\lambda_{i})^{t}\right\|_{2}^{2}\leqslant\sum_{i=1}^{n}\left(\frac{\textbf{v}_{i}(u)}{\sqrt{w(u)}}-\frac{\textbf{v}_{i}(v)}{\sqrt{w(v)}}\right)^{2}(1-\lambda_{i})^{2t}.$

Now the rest of the proof follows from Lemma $4.2$ in [CPS15]. In particular,
it tells us that in the above summation, each of the terms for $1\leqslant
i\leqslant h$ can be bounded by
$\lesssim\frac{\phi_{out}n}{\beta|C|^{3}\phi_{in}^{2}}$ whereas the rest of
the sum can be bounded by $1/\mathrm{poly}(n)$ for sufficiently large
$\mathrm{poly}(n)$ by adjusting the constant in front of $t$. Our choice for
$|C|\geqslant n/\mathrm{poly}(k)$ imply that the overall sum is bounded by
$\lesssim\frac{\phi_{out}\mathrm{poly}(k)}{\beta n^{2}\phi_{in}^{2}}$. Since
$\|D^{1/2}\|^{2}\leqslant n$, and $\phi_{out}\leqslant
c^{\prime}\phi_{in}^{2}$, we have that
$\|p_{u}^{t}-p_{v}^{t}\|_{2}^{2}\leqslant 1/(8n)$, as 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 [CS10].

######  Lemma 6.8.

Let $G$ be a $(k,\phi_{in},\phi_{out})$-clusterable graph with parts
$V=\cup_{1\leqslant i\leqslant h}V_{i}$. There exists a constant $c>0$ such
that if $t\cdot\phi_{out}\leqslant c\varepsilon$, then there exists a subset
$V_{1}^{\prime}\subseteq V_{1}$ satisfying
$\text{vol}(V_{1}^{\prime})\geqslant(1-\varepsilon)\text{vol}(V_{1})$ such
that a $t$-step random walk from any vertex in $V_{1}^{\prime}$ does not leave
$V_{1}$ with probability $1-\varepsilon$.

###### Proof.

We first bound the probability that the random walks always stay in their
respective clusters. Consider a fixed partition $V_{1}$; the same arguments
apply for any partition. Let $G^{\prime}$ be the graph with the same vertex
set as $K$ but with only the following edges: edges among vertices in $V_{1}$
and edges from vertices in $V_{1}$ to $V\setminus V_{1}$. Consider a random
walk on $G^{\prime}$ of length $t$ with the initial vertex $u^{\prime}$ chosen
from the stationary distribution of $G^{\prime}$, i.e., the distribution that
chooses each vertex in $G^{\prime}$ with probability proportional to its
weight. Let $Y_{i}$ denote the indicator random variable for the event that
the $i$th vertex of the random walk is in $V\setminus V_{1}$. Since we are
simulating the stationary distribution, we have that

$\mathbf{Pr}[Y_{i}=1]=\frac{w(V_{1},V\setminus V_{1})}{w(G^{\prime})}$

where $w$ is the weight of edges in the original graph $K$. By linearity of
expectations, the number of vertices that land in $V\setminus V_{1}$ is

$\mathbf{E}\left[\sum_{i=1}^{t}Y_{i}\right]=(t+1)\frac{w(V_{1},V\setminus
V_{1})}{w(G^{\prime})}\lesssim t\phi_{out}$

due to our requirement of $\phi_{out}$. Therefore by Markov’s inequality, the
probability that any vertex in $V\setminus V_{1}$ is ever visited is $\lesssim
t\phi_{out}$.

We now move our random walk back to the original graph $K$. The preceding
calculation implies that the probability that an $t$ step random walk in $K$
starting at a vertex chosen at random from $V_{1}$ according to the stationary
distribution will remain in $V_{1}$ with probability at least $1-\phi_{out}t$.
If $\phi_{out}\lesssim\varepsilon/t$, then we know that the random walk stays
in $V_{1}$ with probability at least $1-O(\varepsilon)$ so there must be a set
of vertices $V^{\prime}\subseteq V_{1}$ of at least $1-O(\varepsilon)$
fraction of the total volume of $V_{1}$ such that a random walk starting from
a vertex in $V^{\prime}$ remains in $V_{1}$ with probability at least
$1-O(\varepsilon)$. ∎

We can now prove the correctness of Algorithm 6.1.

######  Theorem 6.9.

Let $K$ be a $(k,\phi_{in},\phi_{out})$-clusterable kernel graph with parts
$V=\cup_{1\leqslant i\leqslant h}V_{i}$. Let $U,W$ be one of (not necessarily
distinct) partitions $V_{i}$. Let $u,w$ be randomly chosen vertices in
partitions $U$ and $W$ with probability proportional to their degrees. There
exists $c=c(\varepsilon,k)$ such that if $\phi_{out}\leqslant
c\phi_{in}^{2}/\log n$, then with probability at least $1-\varepsilon$, if
$U=W$ then Algorithm 6.1 returns that $u$ and $w$ are in the same cluster and
if $U\neq W$, Algorithm 6.1 returns that $u$ and $w$ are in different
clusters. The algorithm requires
$O(\sqrt{nk/(\varepsilon\tau)}\log(1/\varepsilon))$ random walks of length
$t\geqslant c\log n/\phi_{in}^{2}$.

###### Proof.

We first consider the case that $U\neq W$. From Lemma 6.8, we know that there
are ‘non-escaping’ subsets $U^{\prime}$ and $W^{\prime}$ of $U$ and $W$
respectively such that vertices $u,w$ from $U^{\prime}$ and $W^{\prime}$
respectively don’t leave $U$ and $W$ with probability $1-\varepsilon$.
Conditioning on $u$ and $w$ being in those subsets, we have that with
probability $1-O(\varepsilon)$, $p_{u}^{t}$ and $p_{v}^{t}$ will be disjointly
supported and thus,
$\|p_{u}^{t}-p_{v}^{t}\|_{2}^{2}=\|p_{u}^{t}\|_{2}^{2}+\|p_{v}^{t}\|_{2}^{2}\geqslant
2/n$.

Now if $U=W$, we know from Lemma 6.7 that
$\|p_{u}^{t}-p_{v}^{t}\|_{2}^{2}\leqslant 1/(8n)$ if we condition on $u$ and
$v$ coming from the large volume subset of $U$.

Finally, we need one last ingredient. Lemma $4.3$ in [CPS15] readily implies
that there exists a $V^{\prime}\subseteq V$ satisfying
$\text{vol}(V^{\prime})\geqslant(1-\varepsilon)\text{vol}(V)$ such that
$\|p_{u}^{t}\|_{2}^{2}\leqslant 2k/(\varepsilon\tau^{2}n)$. Now we can set
$\xi=1/(7n)$ and $b=2k/(\varepsilon\tau^{2}n)$ in Theorem 6.5, which tells us
that $r=O(\sqrt{nk/(\varepsilon\tau)}\log(1/\varepsilon))$ samples of the
distributions $p_{u}^{t}$ and $p_{w}^{t}$ suffice to distinguish the cases
$\|p_{u}^{t}-p_{v}^{t}\|_{2}^{2}\geqslant 2/n$ or
$\|p_{u}^{t}-p_{v}^{t}\|_{2}^{2}\leqslant 1/(8n)$, i.e., $r$ samples allow us
to determine if $U=W$ or $U\neq W$, conditioned on a $1-O(\varepsilon)$
probability event. ∎

It is straightforward to translate the requirements of Theorem 6.9 in terms of
the number of KDE queries required. Note that since we only take random walks
of length $O(\log n/\phi_{in}^{2})$, 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 4.15. Alternatively, we can perform
rejection sampling as stated in the proof of Theorem 4.12.

######  Corollary 6.10.

Algorithm 6.1 and Theorem 6.9 require
$\widetilde{O}(c(k,\varepsilon)\sqrt{nk/\varepsilon}\cdot
1/(\tau^{1.5}\phi_{in}^{2}))$ KDE queries (via calls to Algorithm 4.16, which
performs random walks) as well as the same bound for post-processing time.

###  6.2 Spectral Clustering

We present applications to spectral clustering. In data science, spectral
clustering is often the following clustering procedure: (a) compute $k$
eigenvalues of the Laplacian matrix in order, (b) perform $k$-means clustering
on the Laplacian eigenvector embeddings of the vertices.

The theory behind spectral clustering relies on the fact that the Lapalacian
eigenvectors are effective in representing the cluster structure of the
underlying graph. We refer the reader to [VL07] and references within for more
information. For our application to spectral clustering, we 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 is $k$-weakly clusterable if 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 is 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 6.4.

######  Definition 6.11 (Weakly clusterable Graph).

A graph is $(k,\phi_{out})$-clusterable if the following holds: There exists a
partition of the vertex set into $h\leqslant k$ parts $V=\cup_{1\leqslant
i\leqslant h}V_{i}$ such that $\phi_{G}(V_{i})\leqslant\phi_{out}$.

We now prove the following result that says spectral sparsification preserves
cluster structure according to Definition 6.11. We first remark that the
spectral sparsifier obtained in the previous section is a cut sparsifier as
well. Recall that a cut sparsifier is a subgraph that preserves the values
across all cuts up to relative error $1\pm\varepsilon$. The implication
follows immediately by noting that cuts are induced by quadratic forms on the
Laplacian matrix using $\\{-1,1\\}^{n}$ vectors.

######  Theorem 6.12.

Let $G$ be $(k,\phi_{out})$-clusterable and let $G^{\prime}$ be a cut
sparsifier for $G$. Then $G^{\prime}$ is
$(k,(1\pm\varepsilon)\phi_{out})$-clusterable.

###### Proof.

Let $V_{i}$ be one of the $h\leqslant k$ vertex partitions of $G$. Consider
the conductance of $V_{i}$ defined in Definition 6.2. The numerator represents
the value of a cut separating $V_{i}$ and each term in the denominator is the
sum of the degrees of single vertices. Both values are appropriate cuts in the
graph. Since $G^{\prime}$ is a cut sparsifier, this implies that both the
numerator and denominator are preserved up to a $1\pm\varepsilon$ factor and
thus, the entire ratio is also preserved up to a $1\pm O(\varepsilon)$ factor.
∎

Theorem 6.12 implies that the cluster structure of the sparsified graph
$G^{\prime}$ 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^{\prime}$ as it would have on
$G$. Indeed, we verify this empirically in Section 7.

Furthermore, spectral clustering requires us to compute the first $k$
eigenvectors of the Laplacian matrix. Since our sparsifier has few edges, we
can use Theorem $1$ of [MM15], which says (a variant of) the power method can
quickly find good approximations Laplacian eigenvectors if the matrix is
sparse.

######  Theorem 6.13 (Corollary of Theorem $1$ in [MM15] and Theorem 5.3).

Let $L$ be the Laplacian matrix of the sparsifier computed in Theorem 5.3. Let
$u_{1},\cdots,u_{k}$ be the first $k$ eigenvectors of $L$. Using Theorem $1$
of [MM15], we can find $k$ vectors $v_{1},\cdots,v_{k}$ in time
$\widetilde{O}\left(\frac{kn\log n}{\tau^{2}\varepsilon^{2.5}}\right)$ such
that with probability $99/100$,

$|u_{i}^{T}Lu_{i}-v_{i}^{T}Lv_{i}|\leqslant\varepsilon\lambda_{k+1}^{2}$

for all $i\in[k]$.

###  6.3 Arboricity Estimation

######  Algorithm 6.14 (Arboricity Estimation).

Input: Input dataset $X\subset\mathbbm{R}^{d}$ of size $|X|=n$, weight
function $w(\cdot)$, accuracy parameter $\varepsilon$. Operation: 1. Let
$\Delta=\max_{e,e^{\prime}\in E}\frac{w(e)}{w(e^{\prime})}$,
$m=O\left(\frac{n\Delta\log n}{\varepsilon^{2}}\right)$, and
$G^{\prime}=\emptyset$. 2. For $i=1$ to $i=m$: (a) Sample an edge $e$ with
probability $p_{e}=\frac{\widehat{w_{e}}}{\sum\widehat{w_{e}}}$, where
$\widehat{w_{e}}\in[w_{e},2w_{e}]$ (b) Add $e$ to $G^{\prime}$ with weight
$\frac{1}{mp_{e}}$ Output: $\max_{U\subseteq V}d(G^{\prime}_{U})$.

We now apply our algorithmic building blocks to the task of arboricity
estimation. Consider a weighted graph $G=(V,E,w)$. Let $G_{U}$ be an induced
subgraph of $G$ on the subset of nodes $U$. The density of $G_{U}$ is defined
as

$d(G_{U}):=\frac{w(E(G_{U}))}{|U|}$

where $w(E(G_{U}))$ is the sum of the edge weights of $G_{U}$. The arboricity
of $G$ is defined as

$\alpha:=\max_{U\subseteq V}d(G_{U}).$

The arboricity measures the density of the densest subgraph in a graph.
Intuitively, it informs if there is a strong cluster structure among some
subset of the vertices of $G$. Therefore, it is an important primitive in the
analysis of massive graphs with applications ranging from community detection
in social networks, spam link identification, and many more; see [LRJA10] for
a survey of applications and algorithmic results.

Although polynomial time algorithms exist, we are interested in efficiently
approximating the value of $\alpha$ using the building blocks of Section 4.
Inspired by the unweighted version of the arboricity estimation algorithm from
[MTVV15], we first prove the following result.

######  Theorem 6.15.

Let $U^{\prime}=\operatorname*{arg\,max}_{U}d(G_{U}^{\prime})$ and let
$G^{\prime}$ be the output of Algorithm 6.14. Then with probability at least
$1-1/\mathrm{poly}(n)$,

$(1-\varepsilon)\alpha\leqslant
d(G^{\prime}_{U^{\prime}})\leqslant(1+\varepsilon)\alpha.$

Algorithm 6.14 uses $m=\widetilde{O}(n\log n/(\varepsilon^{2}\tau))$ KDE
queries and $O(mn)$ post-processing time.

###### Proof.

Let $U$ be an arbitrary set of $k$ nodes, let $W=\sum_{e\in E}w(e)$, and let
$W_{U}=W\cdot d(G_{U})$ Since $G$ has weight $W$, then the arboricity $\alpha$
satisfies $\alpha\geqslant\frac{W}{n}$, so that

$m\geqslant\frac{\log n}{\alpha\Delta\varepsilon^{2}}.$

Let $X_{i}$ be the random variable denoting the contribution of the $i$-th
sample to weight of the edges in $G^{\prime}_{U}$ and observe that
$\mathbb{E}[X_{i}]=\frac{W}{m}\cdot d(G_{U})$ so that $\mathbb{E}[X]=W\cdot
d(G_{U})$, for $X=\sum_{i=1}^{m}$. Similarly, we have

$\mathbb{E}[X_{i}^{2}]=\sum_{e=(u,v),u,v,\in
U}p_{e}w_{e}^{2}\cdot\frac{1}{p_{e}^{2}m^{2}}=\sum_{e=(u,v),u,v,\in
U}\frac{Ww_{e}}{m^{2}}=\frac{W\cdot W_{U}}{m^{2}}.$

Since $m=\frac{Cn\Delta\log n}{\varepsilon^{2}}$ for an absolute constant
$C>0$, then

$\mathbb{E}[X_{i}^{2}]\leqslant\frac{k\varepsilon^{2}\alpha^{2}}{Cm\log^{2}n}$

and

$\sum_{i=1}^{m}\mathbb{E}[X_{i}^{2}]\leqslant\frac{k\varepsilon^{2}\alpha^{2}}{Cm\log^{2}n}.$

We also have
$X_{i}-\mathbb{E}[X_{i}]\leqslant\frac{\alpha\varepsilon^{2}}{Cm\log^{2}n}$.
Thus by Bernstein’s inequality for sufficiently large $C$,

$\mathbf{Pr}\left[d(G^{\prime}_{U})\geqslant\frac{\alpha}{10}\right]\leqslant
n^{-10k},$

for $d(G_{U})\leqslant\frac{\alpha}{60}$ and

$\mathbf{Pr}\left[|d(G^{\prime}_{U})-d(G_{U})|\geqslant\frac{\varepsilon\alpha}{10}\right]\leqslant
2n^{-10k},$

for $d(G_{U})>\frac{\alpha}{60}$.

Since there are $\binom{n}{k}\leqslant n^{k}$ subsets of $V$ with size $k$,
then by a union bound, we have that with probability at least $1-3n^{-9k}$,
both

$d(G^{\prime}_{U})\leqslant\frac{\alpha}{10}$

for all subsets $U$ with $d(G_{U})\leqslant\frac{\alpha}{60}$ and

$(1-\varepsilon)d(G_{U})\leqslant
d(G^{\prime}_{U})\leqslant(1+\varepsilon)d(G_{U})$

for all subsets $U$ with $d(G_{U})>\frac{\alpha}{60}$.

Hence for a set $U^{*}$ such that $d(U^{*})=\alpha$, we have
$d(G_{U^{*}})\geqslant(1-\varepsilon)\alpha$ so that
$d(G_{U^{\prime}})\geqslant d(G_{U^{*}})\geqslant(1-\varepsilon)\alpha$, where
$U^{\prime}=\operatorname*{arg\,max}_{U\subseteq V}d(G^{\prime}_{U})$. Thus
with high probability, we have that

$(1-\varepsilon)\alpha\leqslant
d(G_{U^{\prime}})\leqslant(1+\varepsilon)\alpha,$

as desired. ∎

To estimate the arboricity of the input graph $G$, it then suffices Theorem
6.15 to compute the arboricity of the subsampled graph $G^{\prime}$ output by
Algorithm 6.14. This can be efficiently achieved by running an offline
algorithm such as [Cha00], which requires solving a linear program on $m$
variables, where $m$ is the number of edges of the input graph. Thus our
subsampling procedure serves as a preprocessing step that ultimately
significantly improves the overall runtime.

###  6.4 Computing the Total Weight of Triangles

We apply the tools developed in prior section to counting the number of
weighted triangles of a 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 shared-
memory and MapReduce to name a few [SPK13, AGM08, BC17, KMPT10, CEI+22].
Applications include discovering motifs in protein interaction networks
[MSOI+02], understanding social networks [FVC10], and evaluating large graph
models [LBKT08]; see the survey [AHD18] for further information.

We define the weight of a triangle as the product of its edges. This
definition is natural since it generalizes the case were the edges have
integer lengths. In this case, an 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
to the case of arbitrary real non-negative weights. This definition has also
been used in definitions of clustering-coefficient for weighted graphs [KH06,
LLL07, AT08].

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 6.16.

Let $G=(V,E,w)$ with $w:E\rightarrow\mathbbm{R}^{\geqslant 0}$ be a weighted
graph. Given a triangle $(x,y,z)\subset E$, we define its weight as

$w_{(x,y,z)}=w(x,y)\cdot w(y,z)\cdot w(x,z),$

where we abuse notation by defining $w(x,y):=w((x,y))$.

For this definition, we present the following modified algorithm from
[ELRS17], which considers the problem in unweighted graphs in a different
model of computing. See Remark 3.7 for comparison.

######  Theorem 6.17.

There exists an algorithm that makes
$\widetilde{O}\left(\frac{m\sqrt{w_{G}}w_{\max}^{3/2}}{w_{T}\cdot\varepsilon^{2}}\right)$
KDE queries and the same bound for post-processing time and with probability
at least $\frac{2}{3}$, outputs a $(1\pm\varepsilon)$-approximation to the
total weight $w_{T}$ of the triangles in the kernel graph.

###### Proof.

Given a graph $G=(V,E)$, let $|V|=n$, $|E|=m$, $\sum_{e\in E}w(e)=w_{G}$, $T$
be the number of triangle in $G$, and $w_{T}$ be the sum of the weighted
triangles in $G$, where the weight of a triangle $(x,y,z)\subset E$, is the
product of the weights of its edges

$w_{(x,y,z)}=w(x,y)\cdot w(y,z)\cdot w(x,z).$

For a vertex $v\in V$, let $w(v)=\sum_{e\in E,e=(u,v),u\in V}w(e)$, so that we