\section{Introduction}
For a kernel function $k: \R^d \times \R^d \rightarrow \R$ and a set $X=\{x_1 \ldots x_n\} \subset \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 \cite{ scholkopf2002learning,shawe2004kernel, hofmann2008kernel} 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~\cite{BackursIS17}.  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~\cite{march2015askit, siminelakis2019rehashing, backurs2019space, backurs2021,  karppa2021deann}. We build on recently developed fast approximate algorithms for Kernel Density Estimation~\cite{charikar2017hashing,backurs2018, siminelakis2019rehashing, backurs2019space,charikar2020}. 
Specifically, these papers present fast approximate data structures with the following functionality:


\begin{definition}[Kernel Density Estimation (KDE) Queries]\label{def:kde_query}
For a given dataset $X \subset \R^d$ of size $n$, kernel function $k$, and precision parameter $\eps>0$, a KDE data structure supports the following operation: given a query $y \in \R^d$, return a value $\kde_X(y)$ that lies in the interval $[(1-\eps)z, (1+\eps) \, z]$, where 
$z = \sum_{x \in X}k(x,y)$, assuming that $k(x,y) \ge \tau$ for all $x \in X$.
\end{definition}
The performance of the state of the art algorithms for KDE also scales proportional to the smallest kernel value of the dataset (see Table \ref{table:kde_instantiation}). 
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) \le 1$ for all inputs $x,y$. 

\begin{table*}[!h]
\centering
\caption{\label{table:kde_instantiation} Instantiations of KDE queries. The query times depend on the dimension $d$, accuracy $\eps$, and lower bound $\tau$. The parameter $\beta$ is assumed to be a constant.}
\begin{tabular}{c|c|c|c|c}
Type               & $k(x,y)$                                                                                                                                 & Preprocessing Time & Query Time & Reference \\ \hline
Gaussian    & $e^{-\|x-y\|_2^2}$     &    $\frac{nd}{\eps^2 \tau^{0.173 + o(1)}}$                & $\frac{d}{\eps^2 \tau^{0.173 + o(1)}}$ &  \cite{charikar2020}         \\
Exponential        & $e^{-\|x-y\|_2}$                                                                                                                          &    $\frac{nd}{\eps^2 \tau^{0.1 + o(1)}}$                & $\frac{d}{\eps^2 \tau^{0.1 + o(1)}}$  &   \cite{charikar2020}        \\
Laplacian & $e^{-\|x-y\|_1}$ & $\frac{nd}{\eps^2 \tau^{0.5}}$ & $\frac{d}{\eps^2 \tau^{0.5}}$ & \cite{backurs2019space} \\
Rational Quadratic & $\frac{1}{(1+\|x-y\|_2^2)^\beta}$ &   $\frac{nd}{\eps^2}$                 & $\frac{d}{\eps^2}$                    &   \cite{backurs2018}       
\end{tabular}

\end{table*}

\subsection{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 \emph{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:
\begin{enumerate}
    \item sampling vertices by their (weighted) degree in $K$ (Theorems \ref{thm:vertex_sampling_1} and \ref{thm:vertex_sampling_2} and Algorithms \ref{alg:vertex_sampling} / \ref{alg:vertex_sampling_complete}),
    \item sampling random neighbors of a given vertex by edge weights in $K$ and sampling a random weighted edge (Theorem \ref{thm:neighbor_sampling} and Algorithms \ref{alg:neighbor_sampling} and \ref{alg:edge_sampling}),
    \item performing random walks in the graph $K$ (Theorem \ref{thm:random_walk} and Algorithm \ref{alg:random_walk}), and
    \item sampling the rows of the edge-vertex incident matrix and the kernel matrix $K$, both with probability proportional to respective row norms squared (Section \ref{sec:spectral_sparsification}, Theorem \ref{thm:spec-spar-kde-matrix}, and Section \ref{sec:lra}, Corollary \ref{cor:additive_lra} respectively).
\end{enumerate}

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 \ref{remark:sparsifier} and Section \ref{sec:initialization}. At a high level, many of our applications, such as spectral sparsification, are succinctly characterized by the following parameterization. 

\begin{parameterization}\label{parameterization:smallest_edge}
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$.
\end{parameterization}

Our applications derived from the basic graph primitives above can be partitioned into two overlapping classes, linear-algebraic and graph theoretic results. Table \ref{table:summary} 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 \emph{all of our theorems below}, we assume access to a KDE data structure of Definition \ref{def:kde_query} with parameters $\eps$ and $\tau$.



\begin{table*}[t]
\caption{\label{table:summary} Summary of linear algebra and graph applications for KDE subroutines. We suppress dependence on the precision $\eps$. 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 \ref{sec:local_clustering}.} 
\centering
\scriptsize


\begin{tabular}{c|c|c|c}
Problem  & $\#$ of KDE Queries & Post-processing time & Prior Work \\\midrule
Spectral sparsification  (Thm. \ref{thm:sparsification_informal})& $\widetilde{O} \left( \frac{n }{\tau^3} \right)$ & $O \left( \frac{nd  }{ \tau^3} \right)$ & Remark \ref{remark:sparsifier} \\
Laplacian system solver (Thm. \ref{thm:sparsification_informal}) & $\widetilde{O} \left( \frac{n  }{ \tau^3} \right)$ & $O \left( \frac{nd }{ \tau^3} \right)$ & Remark \ref{remark:sparsifier}\\
Low-rank approx. (Thm. \ref{thm:lra_informal}) & $O(n)$ & $O\left(n \cdot \poly\left(k\right) + nkd\right)$ & Remark \ref{remark:lra}\\
Eigenvalue Spectrum approx. (Thm. \ref{thm:emd_informal}) & $\widetilde{O}(1/\tau)$ & $O(d/\tau)$ & $\Omega(n^2d)$\\
Approximating 1st Eigenvalue (Thm. \ref{thm:first_eig_informal}) & Remark \ref{remark:first_eig_comparision} & $d \cdot \poly(1/\tau)$ & $\omega(n)$ (Remark \ref{remark:first_eig_comparision}) \\
Local clustering (Thm. \ref{thm:local_clustering_informal})& $\widetilde{O}\left(\poly(k) \cdot \frac{1}{\phi^2} \frac{\sqrt{n}}{\tau^{1.5}}\right)$ & $\widetilde{O}\left(\poly(k) \cdot \frac{1}{\phi^2} \frac{\sqrt{n}}{\tau^{1.5}} \right)$ & Remark \ref{rem:local_clustering}\\
Spectral clustering (Thm. \ref{thm:spectral_clustering}) & $\widetilde{O} \left( \frac{n  }{ \tau^2} \right)$ & $O \left( \frac{nd}{ \tau^2} \right) + \widetilde{O}\left( nk \right)$  & Remark \ref{rem:local_clustering}\\
Arboricity estimation (Thm. \ref{thm:arboricity_informal}) & $\widetilde{O}\left(\frac{n}{\tau}\right)$ & $\widetilde{O}\left(\frac{n^2}{\tau}\right)$  & $\widetilde{O}(n^3) + O(n^2d)$\\
Triangle estimation (Thm. \ref{thm:triangles_informal}) & $\widetilde{O}\left(\frac{1}{\tau^3}\right)$ & $\widetilde{O}\left(\frac{1}{\tau^3}\right)$ & $\Omega(n^2d)$
\end{tabular}
\end{table*}


One of our main results is spectral sparsification of the kernel matrix $K$ interpreted as a weighted graph. In Section \ref{sec:spectral_sparsification}, 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 \cite{merris1994Laplacian, batson2013spectral, spielman2016Laplacian} 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 \ref{parameterization:smallest_edge}.

\noindent
\begin{theorem}[Informal; see Thm. \ref{thm:spec-spar-kde-matrix}]\label{thm:sparsification_informal}
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'$ 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'} \preceq (1+\varepsilon) L$. The algorithm makes $\widetilde{O}(m/\tau^3)$ KDE queries and requires $\widetilde{O}(md/\tau^3)$ post-processing time.
\end{theorem}

We compare our results with prior works in Remark \ref{remark:sparsifier}. We also show that Parameterization \ref{parameterization:smallest_edge} is \textbf{inherent} for spectral sparsification. 
In particular, we use a hardness result from~\cite{alman2020algorithms} to show that for the Gaussian kernel, under the strong exponential time hypothesis~\cite{impagliazzo2001complexity}, 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~\ref{thm:lowerbound-sparsification} 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 \ref{sec:laplacian_system}. 











Continuing the theme of the Laplacian matrix, in Section \ref{sec:EMD}, 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.\ \eqref{eq:emd}), 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 \cite{Cohen-SteinerKS18}.

\begin{theorem}[Informal; see Theorem \ref{thm:emd}]\label{thm:emd_informal}
Let $\eps \in (0, 1)$ be the error parameter and $L$ be the normalized Laplacian of the kernel graph $K$.  Let $\lambda_1 \geq \lambda_2 \ldots \geq \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) \leq \varepsilon$.
\end{theorem}

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^2d)$. Next, we obtain truly sublinear time algorithms for approximating the top eigenvalue and eigenvector of the kernel matrix, a problem which was studied in \cite{backurs2021}. 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 \ref{table:kde_instantiation} for other kernels.

\begin{theorem}[Informal; see Theorem \ref{thm:first_eig}]\label{thm:first_eig_informal}
Given an $n\times n$ kernel matrix $K$ that admits a KDE data-structure with query time $d/(\eps^2 \tau^p)$ (Table \ref{table:kde_instantiation}), there exists an algorithm that outputs a unit vector $v$ such that $v^TKv \ge (1-\eps)\lambda_1(K)$ in time $\min\Paren{\tilde{O}(d/(\eps^{4.5}\tau^4)) , \tilde{O}(d/(\eps^{9+6p} \tau^{2+2p})) }$, where $\lambda_1(K)$ denotes the largest eigenvalue of $K$.
\end{theorem}
We discuss related works in Remark \ref{remark:first_eig_comparision}. In summary, the best prior result of \cite{backurs2021} 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 \ref{sec:lra}.

\begin{theorem}[Informal; see Cor. \ref{cor:additive_lra}]\label{thm:lra_informal}
There exists an algorithm that outputs a rank $r$ matrix $B$ such that $\|K-B\|_F^2 \le \|K-K_r\|_F^2 + \eps \|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 \poly (r,1/\eps) + nrd/\eps)$ post-processing time.
\end{theorem}

We give detailed comparisons between our results and prior work in Remark \ref{remark:lra}. As a summary, \cite{bakshi2020robust} obtain a \emph{relative} error approximation with a running time of $\widetilde{O}\left(n d \left(r/\varepsilon\right)^{\omega-1}  \right)$, where $\omega$ denotes the matrix multiplication constant, whereas our running time is dominated by $O(nrd/\eps)$ and we obtain only additive error guarantees.  Nevertheless, the algorithm we obtain, which builds upon the sampling scheme of \cite{frieze2004fast}, is a conceptually simpler algorithm than the algorithm of \cite{bakshi2020robust} and easier to empirically evaluate. Indeed, we implement this algorithm in Section \ref{sec:experiments} 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~\ref{def:k_clusterable_graph2}. 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.

\begin{theorem}[Informal; see Theorem \ref{thm:localalg_correctness}]\label{thm:local_clustering_informal}
Let $K$ be a $k$-clusterable kernel graph with clusters $V = \cup_{1 \le i \le 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(\eps, k)$ and an algorithm that uses $\widetilde{O}(c(k,\eps) \sqrt{n}/\tau^{1.5})$ KDE queries and post-processing time, with the property that with probability $1-\eps$, if $U = W$ then the algorithm reports that $u$ and $w$ are in the same cluster and if $U \ne W$, the algorithm reports that $u$ and $w$ are in different clusters. 
\end{theorem}

Our definitions for the local clustering result are adopted from prior literature in property testing; see Remark \ref{rem:local_clustering} 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 \ref{sec:spectral_clustering}. 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.

\begin{theorem}[Informal; see Theorem \ref{thm:laplacian_eigv}]
Let $L$ be the Laplacian matrix of the spectral sparsifier. 
There exists an algorithm that can compute $(1+\eps)$-approximations of the first $k$ eigenvectors of $L$ in time $\widetilde{O}\left( kn/(\tau^2 \eps^{2.5})\right)$.
\end{theorem}

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 \ref{sec:arboricity}).

\begin{theorem}[Informal; see Theorem \ref{thm:arboricity}]\label{thm:arboricity_informal}
There exists an algorithm that uses $m = \widetilde{O}(n/(\eps^2 \tau))$ KDE queries and $O(mn)$ post-processing time and outputs a sparse subgraph $G'$ of the kernel graph such that with high probability, $(1 - \eps)\alpha_G \le \alpha_{G'} \le (1+\eps)\alpha_G$, where $\alpha_G$ is the arboricity of $G$.
\end{theorem}

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^2d)$ \cite{GalloGT89}. In Section \ref{sec:triangles}, 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 \cite{kalna2006clustering, li2007structure, antoniou2008statistical}. Our bound is similar in spirit to the bound of the unweighted case given in \cite{EdenLRS17}, under a different computation model. We refer to Remark \ref{remark:triangles} for discussions on related works.



\begin{theorem}[Informal; see Theorem \ref{thm:triangles}]\label{thm:triangles_informal}
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\eps)$-approximation to the total weight of the triangles in the kernel graph. 
\end{theorem}



On the other hand, there is a line of work that considers dimensionality reduction for kernel density estimation e.g., through coresets~\cite{PhillipsT18,PhillipsT20,Tai22}. 
We view this direction of work as orthogonal to our line of study. Lastly, the work \cite{backurs2021} 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, \cite{backurs2021} also study the problem of estimating the sum of all entries in the kernel matrix and obtain tight bounds for the latter. 

\section{Technical Overview}\label{sec: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)$.

\paragraph{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~\ref{fig:tree} for an illustration). 


\paragraph{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 $\eps$-close in total variation to the true distribution (see Theorem \ref{thm:vertex_sampling_2} 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\"{i}ve implementation that uses quadratic time to compute the entire kernel matrix. 

\paragraph{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 $\Pr[v=x_k]=\frac{(1\pm\eps)k(x ,x_k)}{\sum_{j\in [n], x_j \ne 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~\ref{fig:tree} in conjunction with KDE queries.

In particular, consider the tree in Figure \ref{fig:tree} 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 \ref{thm:neighbor_sampling}  and \ref{thm:weighted-edge-sampling} for details). 

\paragraph{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. 

\paragraph{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 \ref{def:leverage-scores}). 
We remark that approximately sampling proportional to the \emph{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 \ref{thm:spec-spar-kde-matrix}). Further, we show that $H$ has a condition number (Lemma \ref{lem:bounding-condition-number}) that is bounded by $\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^TH = 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\"{i}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' = cX$. 
In particular, the $\ell_2^2$ row norms of $K$ are the vertex degrees of the kernel graph for $X'$. 
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. 

\subsection{Linear Algebra Applications} We now discuss our linear algebra applications.

\paragraph{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'$ such that for all vectors $x$, $(1-\varepsilon)x^{T}L_G x \le x^T L_{G'}x \le (1+\varepsilon)x^T L_{G'}x$, where $L_G$ and $L_{G'}$ denote the Laplacian matrices of the graphs $G$ and $G'$. 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$,
\begin{align}
    (1-\varepsilon)L_G = (1-\varepsilon)H^T H \preceq H^TS^TSH \preceq (1+\varepsilon)H^T H = (1+\varepsilon)L_G.
    \label{eqn:spectral-sparsification}
\end{align}
Thus the matrix $SH$ directly corresponds to a graph $G'$, 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 \eqref{eqn:spectral-sparsification}, although we have to sample $t = O(\kappa^2n\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 \emph{positive} singular values).

With the parameterization that for all $i \ne j$, $k(x_i, x_j) \ge \tau$, we are able to show that $\kappa \le 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) \ge \sqrt{\lambda_2(H^TH)} = \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 \cite{friedland2002cheeger} 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$.




\paragraph{First eigenvalue and eigenvector approximation.}
Our goal is to compute a $1 - \eps$ 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 \cite{backurs2021}, 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 \cite{bhattacharjee21,bakshi2020testing} 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/(\eps^2 \tau^2))$ with the guarantee that the top eigenvalue of the sub-sampled matrix is at lest $(1-\eps) \lambda$. Now we can either run the standard Krylov method algorithm~\cite{MuscoM15} to compute the top eigenvalue of the sampled matrix or alternatively, we can instead use the algorithm of \cite{backurs2021}, 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 \cite{backurs2021} have a $\omega(n)$ runtime.



\paragraph{Kernel matrix low-rank approximation.} 
In this setting, our goal here is to output a matrix $B$ such that
\[\|K-B\|_F^2 \le \|K-K_t\|_F^2 + \eps \|K\|_F^2 \]
where $K_t$ is the best rank-$t$ approximation to the kernel matrix $K$. 
The efficient algorithm of \cite{frieze2004fast} is able to achieve this guarantee if one can sample the $i$th row $r_i$ of $K$ with probability $p_i \ge \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. 



\paragraph{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 \cite{Cohen-SteinerKS18} 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 \cite{Cohen-SteinerKS18} on kernel graphs in a black-box manner.


\paragraph{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 $\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.


\subsection{Graph Applications.} We now discuss our graph applications.

\paragraph{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 \cite{CzumajPS15} and \cite{Peng20}) 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 \ref{def:k_clusterable_graph2} 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 \cite{ChanDVV14}. The guarantees of the overall local clustering algorithm of \cite{CzumajPS15} follow for kernel graphs since we only need to access the graph via random walks.


\paragraph{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~\cite{McGregorTVV15}, who observed that to estimate the arboricity on unweighted graphs, it suffices to sample a set of $\tilde{O}\left(\frac{|V|}{\eps^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 \cite{McGregorTVV15} 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\"{i}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. 



\paragraph{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 \cite{EdenLRS17}, 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.




\section{Further Related Works}

\begin{remark}\label{remark:sparsifier}
\normalfont
Spectral sparsification for kernel graphs has also been studied in prior works, notably in \cite{alman2020algorithms} and \cite{Quanrud21}. We first compare to \cite{alman2020algorithms}, 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 \textit{multiplicativily Lipschitz} (see Section 1.1.2 in~\cite{alman2020algorithms} for definition) and show hardness for constructing such a sparsifier when it is not. Focusing on the Gaussian kernel, under Parameterization \ref{parameterization:smallest_edge}, \cite{alman2020algorithms} obtain an algorithm that runs in time $O\left(n d + n 2^{\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, \cite{alman2020algorithms} 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 ~\cite{march2015askit, siminelakis2019rehashing, backurs2019space, backurs2021,  karppa2021deann}.


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

\end{remark}

\begin{remark}\label{remark:first_eig_comparision}
\normalfont
Our algorithm returns a \emph{sparse} vector $v$ supported on roughly $O(1/(\eps^2 \tau^2))$ coordinates. The best prior result is that of \cite{backurs2021} which presented an algorithm with total runtime $O\left(\frac{dn^{1+p}\log(n/\eps)^{2+p}}{\eps^{7+4p}} \right)$.

In comparison, our bound has \emph{no dependence} on $n$ and is thus a \emph{truly sublinear} runtime. Note that the bound of \cite{backurs2021} does not depend on $\tau$. We do not state the number of KDE queries used explicitly in Table \ref{table:summary} since our algorithm uses KDE queries on a subsampled dataset and in addition, only uses them by calling the algorithm of \cite{backurs2021} as a subroutine (on the subsampled dataset). The algorithm of \cite{backurs2021} uses $\tilde{O}(1/\eps)$ 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 \cite{backurs2021} 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 \cite{backurs2021} essentially sets $\tau = 1/\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.
\end{remark}


\begin{remark}\label{remark:lra}
\normalfont
We now compare our low-rank approximation result with a recent work of \cite{musco2017sublinear, bakshi2020robust}. They showed the following theorem:

\begin{theorem}[Theorem 4.2, \cite{bakshi2020robust}]
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(n r/\varepsilon \right)$ entries in $A$ and with probability at least $9/10$, outputs a rank-$r$ matrix $B$ such that 
\begin{equation*}
    \norm{A - B}_F^2 \leq (1+\varepsilon)\norm{A - A_r}_F^2, 
\end{equation*}
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. 
\end{theorem}

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

\begin{fact}[Kernel Matrices are PSD, \cite{scholkopf2002learning}]
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$.
\end{fact}

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 \cite{bakshi2020robust} have comparable runtimes. While~\cite{bakshi2020robust} obtain relative-error guarantees,  we only obtain additive-error guarantees.

However, reading each entry of the kernel matrix require $O(d)$ time and thus \cite{bakshi2020robust} obtain an  running time of $\widetilde{O}\left(n d \left(r/\varepsilon\right)^{\omega-1}  \right)$, whereas our running time is dominated by $O(nrd/\eps)$.  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~\cite{bakshi2018sublinear,IndykVWW19}.
\end{remark}

\begin{remark}\label{rem:local_clustering} \normalfont
Our definitions for the local clustering result are adopted from prior literature in property testing; see \cite{kale2008testing, czumaj2010testing, goldreich2011testing, CzumajPS15, chiplunkar2018testing, dey2019spectral, Peng20, gluch2021spectral} and the references within. Our algorithmic details for the local cluster section are also derived from prior works, such as the works of \cite{CzumajPS15} and \cite{Peng20}; indeed, many of the lemmas of the local clustering section follow in a straightforward fashion from  \cite{CzumajPS15} and \cite{Peng20}. 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. \end{remark}

\begin{remark}\label{remark:triangles}
\normalfont
Our general bound of the number of KDE qeuries required to approximate the total weight of triangles in Theorem \ref{thm:triangles} 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 \cite{EdenLRS17}. 
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 \ref{sec:blocks} play a similar role to the queries used in \cite{EdenLRS17}. 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 \le m w_{\max}$ and $w_T$ is the natural analogue of $T$ in \cite{EdenLRS17}. 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^2d)$ 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~\cite{goldreich2017introduction}. 
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.
\end{remark}


\subsection{Preliminaries}\label{sec:initialization}
First, we discuss the cost of constructing KDE data structure and performing the queries described in Definition \ref{def:kde_query}. 
Table \ref{table:kde_instantiation} 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} \eps^2))$ and answer KDE queries in time $O(d/(\tau^{p} \eps^2))$, under the condition that for all queries $y$ we have $\frac{1}{n} \sum_{x \in X}k(x,y) \ge \tau$ (which clearly holds under our Parameterization \ref{parameterization:smallest_edge}).  
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 \eps^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]$.

\iffalse
\begin{definition}[Fast KDE Query]\label{def:fast_kde_query}
A kernel function $k$ admits a fast KDE query if given a dataset $X \subset \R^d$, and error parameter $\eps$, there exists $p \ge 0$ such that given any query $y \in \R^d$, we can output a $1+\eps$ relative error estimate of $\frac{1}{|X|} \sum_{x \in X}k(x,y)$ with probability $\ge 2/3$ assuming that $\frac{1}{|X|} \sum_{x \in X}k(x,y) \ge \tau$. The query time is $O(d/(\tau^{p} \eps^2))$.
\end{definition}
\fi



\iffalse

There is one key difference between the theoretical guarantees of Definition \ref{def:fast_kde_query} and Table \ref{table:kde_instantiation} and the queries assumed in Definition \ref{def:kde_query}. Namely, Definition \ref{def:fast_kde_query} operates under the condition that the \emph{average} kernel value $\frac{1}{|X|} \sum_{x \in X}k(x,y)$ is at least $\tau$. 
\fi

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 \emph{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 \ref{def:kde_query} 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 \cite{karppa2021deann}.

\iffalse
Such lower bound assumptions are unavoidable if we wish to understand the structure of the kernel graph $K$ without fully constructing it. Consider a clusterable graph with two well connected clusters with a sparse in between. It is clear that identifying the clusters give us valuable insight into the graph. We can arbitrarily scale the edge weights to be extremely small while still preserving the cluster structure. However, now we have little hope of accomplishing anything nontrivial since we are required to multiplicatively estimate edge weights that are arbitrarily close to $0$, which seems prohibitive.
\fi

\iffalse
Lastly, we get non trivial algorithms based on the queries in Definition \ref{def:kde_query}. For example, our algorithms run in time $\widetilde{O}(d/(\tau^p))$ whereas a naive algorithm such as rejection sampling must take time $O(d/\tau)$. For $p < 1$, which is the case the theoretical algorithms in Table \ref{table:kde_instantiation}, this represents significant savings.
\fi

We note that, in addition to the aforementioned algorithms with theoretical guarantees, there are other practical algorithms based on random sampling, space partition trees~\cite{gray2001n, gray2003nonparametric, lee2006dual, lee2008fast, morariu2008automatic, ram2009linear, march2015askit}, coresets~\cite{phillips2013varepsilon, zheng2013quality, phillips2020near}, or combinations of these methods~\cite{karppa2021deann}, which support queries needed in Definition \ref{def:kde_query}; see \cite{karppa2021deann} 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 \ref{table:kde_instantiation}, we can nonetheless use them via black box access in our algorithms and utilize their practical benefits.

\section{Algorithmic Building Blocks}\label{sec:blocks}

\begin{figure*}[!htb]
\centering
\begin{tikzpicture}[scale=0.6]

\draw (0,0) rectangle+(1,1);
\draw (1+1*0.3,0) rectangle+(1,1);
\draw (2+2*0.3,0) rectangle+(1,1);
\draw (3+3*0.3,0) rectangle+(1,1);
\draw (4+4*0.3,0) rectangle+(1,1);
\draw (5+5*0.3,0) rectangle+(1,1);
\draw (6+6*0.3,0) rectangle+(1,1);
\draw (7+7*0.3,0) rectangle+(1,1);
\draw (0.5+0*0.3,1) -- (0.5+0*0.3,1.4);
\draw (1.5+1*0.3,1) -- (1.5+1*0.3,1.4);
\draw (2.5+2*0.3,1) -- (2.5+2*0.3,1.4);
\draw (3.5+3*0.3,1) -- (3.5+3*0.3,1.4);
\draw (4.5+4*0.3,1) -- (4.5+4*0.3,1.4);
\draw (5.5+5*0.3,1) -- (5.5+5*0.3,1.4);
\draw (6.5+6*0.3,1) -- (6.5+6*0.3,1.4);
\draw (7.5+7*0.3,1) -- (7.5+7*0.3,1.4);


\draw (0+1*0.4,1.4) rectangle+(2,1);
\node at (1.3,1.9){$A_{1,n/4}$};
\draw (2+2*0.4,1.4) rectangle+(2,1);
\node at (3.7,1.9){$\ldots$};
\draw (4+3*0.4,1.4) rectangle+(2,1);
\node at (6.1,1.9){$\ldots$};
\draw (6+4*0.4,1.4) rectangle+(2,1);
\node at (8.5,1.9){$\ldots$};
\draw (1+1*0.4,1.4+1) -- (1+1*0.4,2*1.4);
\draw (3+2*0.4,1.4+1) -- (3+2*0.4,2*1.4);
\draw (5+3*0.4,1.4+1) -- (5+3*0.4,2*1.4);
\draw (7+4*0.4,1.4+1) -- (7+4*0.4,2*1.4);

\draw (0+1*0.6,2*1.4) rectangle+(4,1);
\node at (2.7,3.3){$A_{1,n/2}$};
\draw (4+2*0.6,2*1.4) rectangle+(4,1);
\node at (7.1,3.3){$A_{n/2+1,n}$};
\draw (2+1*0.6,2*1.4+1) -- (2+1*0.6,3*1.4);
\draw (6+2*0.6,2*1.4+1) -- (6+2*0.6,3*1.4);

\draw (0+1*0.8,3*1.4) rectangle+(8,1);
\node at (4.7,4.7){$A_{1,n}$};


\end{tikzpicture}
\caption{Multi-level Kernel Density Estimation Data Structure.}
\label{fig:tree}
\end{figure*}


\subsection{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~\ref{alg:kde_construction} for more details. 



\begin{mdframed}
  \begin{algorithm}[Multi-level KDE Construction]
    \label{alg:kde_construction}\mbox{}
    \begin{description}
    \item[Input:] Dataset $X \subset \R^d$, precision $\eps > 0$.
    \item[Operation:]\mbox{}
    \begin{enumerate}
    \item Let $T = X$. 
    \item While $|T|>1$,
    \begin{enumerate}
    \item Construct $\kde_{X}$ queries (see Definition \ref{def:kde_query}). 
    \item Recursively apply Multi-level KDE Construction to $T[1:\lfloor m/2 \rfloor]$ and $T[\lfloor m/2 \rfloor +1:m]$ 
    \end{enumerate}
    \end{enumerate}
     \item[Output:] All the data structures associated with the KDE query constructions
    \end{description}
  \end{algorithm}
\end{mdframed}


\begin{lemma}
Given a dataset $X \subset \R^d$, suppose the initialization of the KDE data structure defined in Definition \ref{def:kde_query} uses runtime $f(n,\eps)$ for some function linear in $n$. 
Then the total construction time of Algorithm~\ref{alg:kde_construction} is $f(n \log n, \eps)$.
\end{lemma}
\begin{proof}
The proof follows from the fact that at each recursive level, we do $O(f(n,\eps))$ total work since $f$ is linear in $n$ and there are $O(\log n)$ levels.
\end{proof}

\subsection{Weighted Vertex Sampling}
We now discuss our fundamental primitives. 
The first one computes approximate weighted degrees for all vertices. Algorithm \ref{alg:vertex_sampling_complete} then performs vertex sampling by their (weighted) degree. 


\begin{mdframed}
  \begin{algorithm}[Computing Approximate (Weighted) Degrees]
  \label{alg:vertex_sampling}\mbox{}
    \begin{description}
    \item[Input:] Dataset $X \subset \R^d$, precision $\eps > 0$.
    \item[Operation:]\mbox{}
    \begin{enumerate}
    \item For $i  \in [1,n]$ :
    \begin{enumerate}
    \item $p_i \gets \kde_X(x_i)- (1-\eps) \, k(x_i, x_i)$
    \end{enumerate}
    \end{enumerate}
     \item[Output:] Reals $\{p_i\}_{i=1}^n$ such that $(1-\eps) \dg(x_i) \le p_i \le (1+\eps) \dg(x_i)$ for all $1 \le i \le n$
    \end{description}
  \end{algorithm}
\end{mdframed}

\begin{definition}[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 $\Pr[v=x_i]=\frac{(1\pm\eps)w_i}{\sum_{j\in [n]} w_j}$ for all $i\in[n]$. 
\end{definition}

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 \emph{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.



\begin{mdframed}
  \begin{algorithm}[Sample from Positive Array]
  \label{alg:array_sampling}\mbox{}
    \begin{description}
    \item[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 \le i \le j \le n$.
    \item[Operation:]\mbox{}
    \begin{enumerate}
    \item Let $T = A$. While $|T| > 1$:
    \begin{enumerate}
    \item Let $m = \text{len}(T)$.
    \item Let $a \gets \sum(T[1:\lfloor m/2 \rfloor])$ \texttt{ //Can be simulated using an $A_{i,j}$ query}
    \item Let $b \gets \sum(T[\lfloor m/2 \rfloor +1:m])$. 
    \item If Unif$[0,1] \le a/(a+b)$, $T \gets T[1:\lfloor m/2 \rfloor]$
    \item Else $T \gets T[\lfloor m/2 \rfloor+1:m]$.
    \end{enumerate}
    \end{enumerate}
     \item[Output:] The single remaining element in $T$
    \end{description}
  \end{algorithm}
\end{mdframed}

Combining Algorithms \ref{alg:vertex_sampling} and \ref{alg:array_sampling}, we can sample from the degree distribution of the graph $K$.



\begin{mdframed}
  \begin{algorithm}[Degree Sampling]
  \label{alg:vertex_sampling_complete}\mbox{}
    \begin{description}
    \item[Input:]  Dataset $X \subset \R^d$, precision $\eps > 0$.
    \item[Operation:]\mbox{}
    \begin{enumerate}
    \item Use Algorithm \ref{alg:vertex_sampling} to compute reals $\{p_i\}_{i=1}^n$ such that $(1-\eps) \dg(x_i) \le p_i \le (1+\eps) \dg(x_i)$ for all $1 \le i \le n$ (only needs to be done once).
    \item $i \gets$ index in $[n]$, which is the output of running Algorithm \ref{alg:array_sampling} on the array $\{p_i\}_{i=1}^n$.
    \end{enumerate}
     \item[Output:] $x_i$ with probability $p_i /\sum_{j \in [n]} p_j$.
    \end{description}
  \end{algorithm}
\end{mdframed}


We now analyze the correctness and the runtimes of the algorithms proposed in Section \ref{sec:blocks}. First, we give guarantees on Algorithm \ref{alg:vertex_sampling}.

\begin{theorem}\label{thm:vertex_sampling_1}
Algorithm \ref{alg:vertex_sampling} returns $\{p_i\}_{i=1}^n$ such that $(1-\eps) \dg(x_i) \le p_i \le (1+\eps) \dg(x_i)$ for all $1 \le i \le n$.
\end{theorem}
\begin{proof}
The proof follows by the Definition of a KDE query, Definition \ref{def:kde_query}.
\end{proof}
We now analyze Algorithm \ref{alg:array_sampling}, 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 \ref{sec:blocks}.


\begin{lemma}\label{lem:array_sampling}
Algorithm \ref{alg:array_sampling} samples an index $i \in [n]$ proportional to $a_i$ in $O(\log n)$ time with $O(\log n)$ queries. 
\end{lemma}
\begin{proof}
Consider the sampling diagram given in Figure \ref{fig:tree}. Algorithm \ref{alg:array_sampling} 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 \ref{fig:tree} 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.
\end{proof}

Combining Algorithms \ref{alg:vertex_sampling} and \ref{alg:array_sampling} allows us to sample from the degree distribution of the graph $K$ up to low error in total variation (TV) distance.

\begin{theorem}\label{thm:vertex_sampling_2}
Algorithm \ref{alg:vertex_sampling_complete} samples from the degree distribution of $K$ up to TV error $O(\eps)$ using a fixed overhead of $n$ KDE queries and runtime $O(\log n)$.
\end{theorem}
\begin{proof}
Since $p_i$ is with a $1\pm \eps$ factor of $\dg(x_i)$ for all $i$, then $\{p_i\}_{i=1}^n$ is $O(\eps)$ close in total variation distance from the true degree distribution. 
Moreover, Algorithm \ref{alg:array_sampling} 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 \ref{alg:array_sampling} by computing all the partial prefix sums in $O(n)$ time. 
Now the query access required by Algorithm \ref{alg:array_sampling} 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 \emph{once} and can be utilized for all future runs of Algorithm \ref{alg:array_sampling}. 
It follows from Lemma \ref{lem:array_sampling} that Algorithm \ref{alg:vertex_sampling_complete} takes $O(\log n)$ time. 
\end{proof}

\subsection{Weighted Edge Sampling and Weighted Neighbor Edge Sampling}
We describe how to perform weighted neighbor edge sampling. 
\begin{definition}[Weighted Neighbor Edge Sampling]
Given a vertex $x_i$, the goal of weighted neighbor edge sampling is to output a vertex $v$ such that $\Pr[v=x_k]=\frac{(1\pm\eps)k(x_i,x_k)}{\sum_{j\in n, j \ne i}k(x_i,x_j)}$ for all $i\in[n]$. 
\end{definition}


\begin{mdframed}
  \begin{algorithm}[Sample Random Neighbor]
  \label{alg:neighbor_sampling}\mbox{}
    \begin{description}
    \item[Input:]  Dataset $X \subset \R^d$, precision $\eps > 0$,  input vertex $x_i \in X$.
    \item[Operation:]\mbox{}
    \begin{enumerate}
    \item Let $\eps' = \eps/\log  n$ and $T \gets X \setminus \{x_i\}$. 
    \item While $|T|>1$ 
    \begin{enumerate}
        \item Let $m \gets |T|$. 
        \item Compute $a \gets \kde_{T[1:m/2], \eps'}(x_i)$ and $b \gets \kde_{T[m/2+1:m] , \eps'}(x_i)$. 
        \item If $x_i \in T[1:m/2]$, set $a \gets a- (1-\eps')k(x_i, x_i)$.
        \item If $x_i \in T[m/2+1:m]$, set $b \gets b- (1-\eps')k(x_i, x_i)$.
        \item If  Unif$[0,1] \le a/(a+b)$, let $T \gets T[1:m/2]$. Else, let $T \gets T[m/2+1:m]$. 
    \end{enumerate}
    \end{enumerate}
     \item[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)$.
    \end{description}
  \end{algorithm}
\end{mdframed}

We now prove the correctness of Algorithm \ref{alg:neighbor_sampling} based on the ideas in Lemma \ref{lem:array_sampling}. 
Note that Algorithm \ref{alg:neighbor_sampling} takes in input a precision level $\eps$, 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 \ref{sec:initialization}.

\begin{theorem}\label{thm:neighbor_sampling}
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 \ref{alg:neighbor_sampling} samples a neighbor from a distribution that is within TV distance $O(\eps)$ 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.
\end{theorem}

\begin{proof}
The proof idea is similar to that of Lemma \ref{lem:array_sampling}. 
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 \ref{lem:array_sampling} in conjunction with KDE queries with over various subsets of $X$.

Imagine the same tree as in Figure \ref{fig:tree} 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 \ref{alg:neighbor_sampling} descends down the tree using the same probabilistic procedure as in the proof of Lemma \ref{lem:array_sampling}: 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 \ref{fig:tree} 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 \ref{alg:edge_sampling}. 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+\eps'$ or $1-\eps'$ 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\eps')^\ell$ factor. However, since $\eps' = \eps/\log n$ and $\ell \le \log n$, this factor is within $1 \pm \eps$. Thus, we sample from the correct distribution over the leaves of the trees in Figure \ref{fig:tree} up to TV distance $O(\eps)$. 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 $\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.
\end{proof}



\begin{mdframed}
  \begin{algorithm}[Sample Random Edge by Weight]
  \label{alg:edge_sampling}\mbox{}
    \begin{description}
    \item[Input:]  Dataset $X \subset \R^d$, precision $\eps > 0$.
    \item[Operation:]\mbox{}
    \begin{enumerate}
    \item Compute $x_i \gets$ random vertex by using Algorithm \ref{alg:vertex_sampling_complete}.
    \item Compute $x_j \gets$ random Neighbor of $x_i$ using Algorithm \ref{alg:neighbor_sampling}.
    \end{enumerate}
     \item[Output:] Edge $(x_i, x_j)$ such that $(x_i, x_j)$ is sampled with probability at least $(1-\eps) k(x_i, x_j)$.
    \end{description}
  \end{algorithm}
\end{mdframed}

\begin{theorem}[Weighted Edge Sampling]
\label{thm:weighted-edge-sampling}
Algorithm \ref{alg:edge_sampling} returns a random edge of $K$ with probability proportional to at least $(1-\eps)$ its weight using $1$ call to Algorithm \ref{alg:neighbor_sampling}.
\end{theorem}
\begin{proof}
Consider an edge $(u,v)$. Vertex $u$ is sampled with probability at least $(1-2 \eps)\frac{\dg(u)}{\sum_{x \in X} \dg(x)}$. Given this, $v$ is then sampled with probability at least $(1-2\eps) \frac{k(u,v)}{\sum_{x \in X \setminus {u}} k(u,x)} = (1-2\eps) \frac{k(u,v)}{\dg(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\eps$ times $k(u,v)/\sum_{e \in K} w(e)$. Note that the same rejection sampling remark as in the proof of Theorem \ref{thm:neighbor_sampling} applies and we can perfectly sample an edge proportional to its weight with an addition $O(1/\tau)$ rejection sampling steps.
\end{proof}

\subsection{Random Walk}

\begin{theorem}\label{thm:random_walk}
Algorithm \ref{alg:random_walk} outputs a vertex from a vertex within $O(T\eps)$ total variation distance from the true random walk distribution. Each step of the walk requires $1$ call to Algorithm \ref{alg:neighbor_sampling} .
\end{theorem}
\begin{proof}
The proof follows from the correctness of Algorithm \ref{alg:neighbor_sampling} given in Theorem \ref{thm:neighbor_sampling}. 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 \ref{thm:neighbor_sampling}), we can make sure that we are sampling from the \emph{true} random walk distribution at each step of the walk.
\end{proof}

\begin{mdframed}
  \begin{algorithm}[Perform Random Walk]
  \label{alg:random_walk}\mbox{}
    \begin{description}
    \item[Input:]  Dataset $X \subset \R^d$, vertex $x_i \in X$, length of walk $T \ge 1$.
    \item[Operation:]\mbox{}
    \begin{enumerate}
    \item Start at vertex $v \gets x_i$.
    \item For $j = 1$ to $T$:
    \begin{enumerate}
        \item Sample a random neighbor of $v$ using Algorithm \ref{alg:neighbor_sampling}. Let $w$ be the resulting output. 
        \item Set $v \gets w$.
    \end{enumerate}
    \end{enumerate}
     \item[Output:] Data point $v$. 
    \end{description}
  \end{algorithm}
\end{mdframed}


\section{Linear Algebra Applications}\label{sec:linalgebra_applications}
We now present a wide array of applications of the algorithmic building blocks constructed in Section \ref{sec:blocks}. 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 \ref{sec:graph_applications}.


\subsection{Spectral Sparsification}\label{sec:spectral_sparsification}





\begin{mdframed}
  \begin{algorithm}[Spectral Sparsification of the Kernel Graph]
  \label{alg:spectral_sparsification}\mbox{}
    \begin{description}
    \item[Input:]  Dataset $X \subset \R^d$, accuracy parameter $\eps$. 
    \item[Operation:]\mbox{}
    \begin{enumerate}
    \item Let $t = O(n\log(n)/\varepsilon^2\tau^3)$ be the number of edges that are to be sampled
    \item Let $\hat{p}$ denote the distribution returned by Algorithm~\ref{alg:vertex_sampling} for a small enough constant $\varepsilon$. 
    \item Initialize $G' = \emptyset$. For $i = 1,\ldots,t$:
    \begin{enumerate}
        \item Sample a vertex $u$ from the distribution $\hat{p}$.
        \item Sample a neighbor $v$ of $u$ using Algorithm~\ref{alg:neighbor_sampling} with constant $\varepsilon$.
        \item Compute $\hat{q}_{uv}$, the probability that Algorithm~\ref{alg:neighbor_sampling} samples $v$ given $u$ as input.
        \item Similarly define and compute $\hat{q}_{vu}$. Let $w_{uv} = 1/(t(\hat{p}_u\hat{q}_{uv} + \hat{p}_v\hat{q}_{vu}))$. 
        \item Add the weighted edge $(\{u,v\}, w_{uv})$ to the graph $G'$.
    \end{enumerate}
   
    \end{enumerate}
     \item[Output:] 
    \end{description}
  \end{algorithm}
\end{mdframed}


Given a set $X$, $|X| = n$, and a kernel $k : X \times X \rightarrow \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)$.

\begin{definition}[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 \ne i} A_{i,j}$.   
\end{definition}

\begin{theorem}[Spectral Sparsification of Kernel Density Graphs]\label{thm:spec-spar-kde-matrix}
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) \geq \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~\ref{alg:spectral_sparsification} outputs a graph $G'$ with only $m = O( {n \log n }/{(\varepsilon^2 \tau^3)} )$ edges, such that with probability at least $9/10$,
\begin{equation*}
    (1-\varepsilon) L_G \preceq L_{G'} \preceq (1+\varepsilon) L_G.
\end{equation*}
The algorithm makes $\widetilde{O}(m/\tau^3)$ KDE queries and requires $\tilde{O}(md/\tau^3)$ post-processing time. 
\end{theorem}


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 $\ge 9/10$, satisfies
\begin{equation*}
    (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.
\end{equation*}
The edges sampled by $S$ form the edges of the graph $G'$. We construct this matrix $S$ by sampling rows of the matrix $H$ from 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. 
\begin{definition}[Leverage Scores]
\label{def: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:
\begin{equation*}
    \tau_i = m_i ( M^\top M)^{+} m_i^\top,
\end{equation*}
where $X^+$ is the Moore-Penrose pseudoinverse for a matrix $X$. 
\end{definition}

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 \cite{spielman2011graph,dw-tool}.
\begin{lemma}[Leverage Score Sampling implies Sparsification]
\label{lem:sampling-to-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 \geq \phi p_i$ and let $t = O\left( \frac{d \log(d) }{\varepsilon^2 \phi} \right)$. Let $S \in \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$,
\begin{equation*}
    (1-\varepsilon) M^\top M \preceq M^\top S^\top S M \preceq (1+\varepsilon)   M^\top M  . 
\end{equation*}
\end{lemma}





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.

\begin{lemma}[Bounding Condition Number]
\label{lem: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) \leq 4\sqrt{2}/\tau^{1.5}$.
\end{lemma}
\begin{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$:
\begin{equation*}
    \opnorm{A} \le \sqrt{\left(\max_i \sum_j |A_{i,j}|\right)\left(\max_j \sum_i |A_{i,j}|\right)}.
\end{equation*}
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 $\opnorm{H} \le \sqrt{2n}$. To obtain lower bounds on $\sigma_{\min}(H)$, we appeal to a Cheeger-type inequality for weighted graphs from \cite{friedlandcheegeroriginal, friedland2002cheeger}. 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 \le \lambda_2 \le \cdots \le \lambda_n$ be the eigenvalues of the positive semi-definite matrix $L_G$. Now we have that
\begin{equation*}
    \sigma_{\min}(L_G) = \lambda_2(L_G) \ge \min_i (\delta_i/2)\varepsilon(G)^2
\end{equation*}
where $\delta_i = \sum_{j \ne i}k(x_i, x_j)$ i.e., the weighted degree of vertex $x_i$ in graph $G$ and
\begin{equation*}
    \varepsilon(G) = \min_{\phi \ne U \subset V, |U| \le n/2} \frac{|E(U, \bar{U})|}{|E(U)|}
\end{equation*}
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})| \ge \tau|U||\bar{U}|$ and $|E(U)| \le n|U|$, which implies that $\varepsilon(G) \ge \min_{\phi \ne U \subset V, |U| \le n/2} \tau|\bar{U}|/n \ge \tau/2$. We also similarly have that $\min_i \delta_i \ge (n-1)\tau$, which overall implies that $\lambda_2(L_G) \ge n\tau^3/16$ and that $\sigma_{\min}(H) \ge \sqrt{n}\tau^{1.5}/4$. Thus, we obtain that $\sigma_{\max}(H)/\sigma_{\min}(H) \le 4\sqrt{2}/\tau^{1.5}$.
\end{proof}


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

\begin{proof}[Proof of Theorem \ref{thm:spec-spar-kde-matrix}] 
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 \geq c \frac{ \norm{H_{e,*} }_2^2 }{ \norm{H}_F^2} = \frac{k(x_i, x_j)}{\sum_{e' = \{i',j'\}}k(x_{i'}, x_{j'})}$, 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 = \opnorm{U_{i*}}^2$. We have
\begin{align*}
    \opnorm{h_i}^2 = \opnorm{U_{i*}\Sigma V^\top}^2 = \opnorm{U_{i*}\Sigma}^2
\end{align*}
where the equality follows from the fact that $V^\top$ has orthonormal rows. Now, $\opnorm{U_{i*}\Sigma}^2 \ge \opnorm{U_{i*}}^2\sigma_{\min}^2$ and $\opnorm{U_{i*}\Sigma}^2 \le \opnorm{U_{i*}}^2\sigma_{\max}^2$. Therefore, for all $i \in {\binom{n}{2}}$, defining $\kappa = \sigma_{\min}/\sigma_{\max}$, we have
\begin{equation*}
    \frac{\tau_i}{\sum_j \tau_j} = \frac{\opnorm{U_{i*}}^2}{\sum_{j} \opnorm{U_{j*}}^2} \ge \frac{\opnorm{h_i}^2/\sigma_{\max}^2}{\sum_j \opnorm{h_j}^2/\sigma_{\min}^2} = \frac{1}{\kappa^2}\frac{\opnorm{h_i}^2}{\frnorm{H}^2}.
\end{equation*}
Then, we invoke Lemma \ref{lem:sampling-to-sparsification} 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'$ with corresponding Laplacian $L_{G'}$ such that with probability at least $99/100$, 
\begin{equation*}
    (1-\varepsilon/2) L_G \preceq L_{G'} \preceq (1+\varepsilon/2) L_G.
\end{equation*}

Further, by Lemma \ref{lem:bounding-condition-number}, we can conclude $\kappa^2 \leq 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~\ref{alg:edge_sampling} 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~\ref{alg:vertex_sampling_complete} with constant $\varepsilon$ (say 1/2) to sample a vertex $u$ and Algorithm~\ref{alg:neighbor_sampling} with constant $\varepsilon$ (say 1/2) to sample a neighbor $v$ of $u$. Note that Algorithms~\ref{alg:vertex_sampling_complete} and Algorithms~\ref{alg:neighbor_sampling} 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} \ge c \frac{k(x_u, x_v)}{\sum_{i \ne j} k(x_i, x_j)}$ and we use this distribution $q$ to implement spectral sparsification as described above. As already seen (Theorem~\ref{thm:neighbor_sampling}), 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. 
\end{proof}
We can further compute another graph $G''$ with only $O(n/\varepsilon^2)$ edges by computing an $\varepsilon/2$ spectral sparsifier for $G'$ using the  spectral sparsification algorithm of Lee an Sun~\cite{lee2018constructing} (see Theorem 1.1). 
This procedure doesn't require any KDE queries and solely operates on the weighted graph $G'$. The overall running time is $O\Paren{ \frac{n^{1+1/c} \log(n) }{\epsilon^{6 +1/c} \tau^3 } }$, for a large fixed constant $c$.  

\paragraph{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 \ref{parameterization:smallest_edge}. The lower bound we obtain is as follows:

\begin{theorem}[Lower Bound for Spectral Sparsification under Parameterization \ref{parameterization:smallest_edge}]
\label{thm:lowerbound-sparsification}
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. 
\end{theorem}


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

\begin{definition}[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) \leq k(\rho x, \rho y) \leq c^L k(x, y)$.
\end{definition}

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

\begin{theorem}[Theorem 8.3~\cite{alman2020algorithms}]
\label{thm:alman-et-al}
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.
\end{theorem}


\begin{proof}[Proof of Theorem \ref{thm:lowerbound-sparsification} ]
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 = \norm{x- y}_2^2$ and let $f(z) = e^{-z}$.   Observe, it suffices to show that there exists a $z$ such that $f(cz) \leq 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, 
\begin{equation*}
    f(c \log(1/\tau) ) = e^{-c \log(1/\tau)},
\end{equation*}
and for $L < \log(1/\tau)(c-1)$
\begin{equation*}
    c^{-L} f(\log(1/\tau)) >   e^{-c \log(1/\tau)}.
\end{equation*}
Then, applying Theorem \ref{thm:alman-et-al} 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. 
\end{proof}


\subsubsection{Solving Laplacian Systems Approximately}\label{sec:laplacian_system}
We describe how to approximately solve the Laplacian system $L_G x = b$ using the spectral sparsifier $L_{G'}$. First, we note the following theorem that states the running time and approximation guarantees of fast Laplacian solvers.
\begin{theorem}[\cite{nearlinear-kmp}, \cite{nearlinear-spielman-teng}]
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$,
\begin{equation*}
    \|x - L^{+}b\|_L \le \alpha\|L^+ b\|_L,
\end{equation*}
where $\|x\|_L = \sqrt{x^\top Lx}$. The algorithm runs in time $\widetilde{O}(m\log(1/\alpha))$.
\end{theorem}
We have the following theorem that bounds the difference between solutions for the exact Laplacian system and the spectral sparsifier Laplacian.
\begin{theorem}\label{thm:laplacian_system_solving}
Let $L_G$ be the Laplacian of a connected graph $G$ on $n$ vertices and let $L_{G'}$ be the Laplacian of an $\varepsilon$-spectral sparsifier $G'$ of graph $G$ i.e., 
\begin{align*}
    (1 - \varepsilon)L_G \preceq L_{G'} \preceq (1+\varepsilon) L_G,
\end{align*}
for $\varepsilon < c$ for a small enough constant $c$. Then, for any vector $b$ with $1^\top b = 0$, $\|L_G^+b - L_{G'}^+b\|_{L_G} \le 2\sqrt{\varepsilon}\|L_{G}^+ b\|_{L_G}$.
\end{theorem}
\begin{proof}
Note that for $\varepsilon < 1$, the graph $G'$ also has to be connected and therefore the only eigen vectors corresponding to eigen value $0$ of the matrices $L_G$ and $L_{G'}$ are of the form $a \cdot 1$ for $a \ne 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,
\begin{align*}
    \|L_G^+b - L_{G'}^+b\|_{L_G}^2 &= b^\top (L_G^+ - L_{G'}^+)L_G(L_G^+ - L_{G'}^+)b\\
    &= b^\top L_G^+ L_G L_G^+ b - b^\top L_{G'}^+ L_GL_{G}^+b - b^\top L_{G'}^+L_GL_G^+b + b^\top L_{G'} ^+ L_G L_{G'}^+ b\\
    &\le b^\top L_{G}^+ b - b^\top L_{G'}^+b - b^{\top}L_{G'}^+b + \frac{1}{1-\varepsilon}b^{\top}L_{G'}^+b
\end{align*}
where in the last inequality, we used $L_{G}L_{G}^+b = 1$ and that for any vector $v$, $v^\top L_G v \le \frac{1}{1-\varepsilon} v^\top L_{G'}v$. As the null spaces of both $L_G$ and $L_{G'}$ are given by $\{a1\,|\, a \in \mathbb{R}\}$, we also obtain that
\begin{align*}
    (1-\varepsilon)L_G^{+} \preceq L_{G'}^+ \preceq (1+\varepsilon)L_G^+
\end{align*}
using which we further obtain that
\begin{align*}
    \|L_G^+b - L_{G'}^+b\|_{L_G}^2 &\le \left(\frac{2}{1-\varepsilon} - 2\right)b^{\top}L_{G'}^+b \le \frac{2\varepsilon(1+\varepsilon)}{1-\varepsilon}b^\top L_G^+ b \le 4\varepsilon\|L_{G}^+b\|_{L_G}^2.
\end{align*}
Thus, $\|L_G^+b - L_{G'}^+b\|_{L_G} \le 2\sqrt{\varepsilon}\|L_{G}^+ b\|_{L_G}$.
\end{proof}
Therefore, if $x$ is a vector such that $\|x - L_{G'}b\|_{L_{G'}} \le \alpha\|L_{G'}^+b\|_{L_{G'}}$ obtained using the fast Laplacian solver, then 
\begin{align*}
    \|x - L_G^+b\|_{L_G}^2 &= \|x - L_{G'}^+b + L_{G'}^+b - L_{G}^+b\|_{L_G}^2\\
    &\le 2(\|x- L_{G'}^+b\|_{L_G}^2  + \|L_{G'}^+b - L_{G}^+b\|_{L_G}^2)\\
    &\le \frac{2}{1-\varepsilon}\|x - L_{G'}^+b\|_{L_G'}^2 + 4\varepsilon\|L_{G}^+b\|_{L_G}^2.
\end{align*}
Here we used the above theorem and the fact that $L_{G} \preceq (1/(1-\varepsilon))L_{G'}$. Now, $\|x - L_{G'}^+b\|_{L_{G'}}^2 \le \alpha^2\|L_{G'}^+b\|_{L_{G'}}^2$ and $\|L_{G'}^+b\|_{L_{G'}}^2 = b^{\top} L_{G'}^+L_{G'}L_{G'}^+b = b^\top L_{G'}^+ b \le (1+\varepsilon)b^\top L_G^+ b \le (1+\varepsilon)\|L_{G}^+b\|_{L_G}^2$, which finally implies that
\begin{align*}
    \|x - L_{G}^+b\|_{L_G}^2 \le \left(\frac{2(1+\varepsilon)^2}{1-\varepsilon}\alpha^2 + 4\varepsilon\right)\|{L_{G}^+b}\|_{L_G}^2.
\end{align*}
Thus, using a $\varepsilon$ spectral sparsifier $G'$ 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} \le C\sqrt{\varepsilon}\|L_{G}^+b\|_{L_G}$ for a large enough constant $C$.







\subsection{Low-rank Approximation of the Kernel Matrix}\label{sec:lra}
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$.

\begin{theorem}[\cite{frieze2004fast}]\label{thm:fkv}
Let $A \in \R^{n\times m}$ be any matrix. Let $S$ be a sample of $O(k/\eps)$ rows according to a probability distribution $(p_1, \ldots,p_n)$ that satisfies $p_i \ge \Omega(1) \cdot \|A_{i,*}\|_2^2/ \|A|_F^2$ for every $1 \le i \le n$. Then, in time $O(m k/\eps \cdot \poly(k,1/\eps))$, we can compute from $S$ a matrix $U \in \R^{k \times m}$, that with probability at least $0.99$ satisfies
\[ \|A-AU^TU \|_F^2 \le \|A-A_k\|_F^2 + \eps \|A\|_F^2. \]
\end{theorem}

\begin{theorem}[\cite{chen2017}, also see \cite{IndykVWW19}]\label{thm:chen}
There is a randomized algorithm that given matrices $A \in \R^{n \times m}$ and $U \in \R^{k \times m}$, reads only $O(k/\eps)$ columns of $A$, runs in time $O(mk) + \poly (k,1/\eps)$, and returns $V\in \R^{n \times k}$ that with probability $0.99$ satisfies 
\[ \| A - VU\|_F^2 \le (1 + \eps)  \min_{X \in \R^{n\times k}} \| A - X \|_F^2.\]
\end{theorem}

Therefore to compute the low rank approximation, we just need sample from the distribution on rows required by Theorem \ref{thm:fkv}. 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 \emph{squared} for vertex $x_i$. From vertex queries (Algorithm \ref{alg:vertex_sampling_complete}), we know that we can get the degree of each vertex, which is the sum of edge weights. We can extend Algorithm \ref{alg:vertex_sampling_complete} to sample from the sum of \emph{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' := cX$. Then by sampling the degrees of the vertices associated with the kernel graph $K'$ of $X'$, we can sample from the distribution required by Theorem \ref{thm:fkv} by invoking Algorithm \ref{alg:vertex_sampling_complete} on the dataset $X'$. In particular, using $n$ KDE queries for $X'$, 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 \ref{thm:fkv} and fully construct the rows that are sampled. Altogether, this takes $n$ KDE queries and $O(nk/\eps)$ kernel function evaluations to construct a rank $k$ approximation of $K$; see Algorithm \ref{alg:additive_lra}.

\begin{corollary}\label{cor:additive_lra}
Given a dataset $X$ of size $n$, there exists an algorithm that outputs a rank $k$ matrix $B$ such that
\[\|K-B\|_F^2 \le \|K-K_k\|_F^2 + \eps \|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/\eps \cdot \poly (k,1/\eps) + nkd/\eps)$ post-processing time.
\end{corollary}

We remark that for the application presented in this subsection, we can we can replace \ref{parameterization:smallest_edge}. 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) \ge \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 \ge \frac{1}{|X|^2} \left( \sum_{x \in X} k(x,y) \right)^2 \ge \tau^2.\]




\begin{mdframed}
  \begin{algorithm}[Additive-error Low-rank Approximation]
  \label{alg:additive_lra}\mbox{}
    \begin{description}
    \item[Input:]  Kernel matrix $K \in \R^{n \times n}$, data points $X \subset R^{d}$, accuracy parameter $\eps$, rank parameter $k$. 
    \item[Operation:]\mbox{}
    \begin{enumerate}
    \item 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$:
    \begin{enumerate}
        \item Compute the value $p_i = \sum_{j=1}^n k(cx_i, cx_j)$ using KDE queries for the dataset $cX$.
    \end{enumerate}
    \item Sample and construct $O(k/\eps)$ rows of $K$ according to probability proportional to $\{p_i\}_{i=1}^n$.
    \item Compute $U$ from the sample, using Theorem \ref{thm:fkv}.
    \item Compute $V$ from the sample, using Theorem \ref{thm:chen}. 
    \end{enumerate}
     \item[Output:] Factors $U, V$ such that $\norm{K - UV}_F^2 \leq \norm{K - K_k }_F^2 + \eps \norm{K}_F^2 $
    \end{description}
  \end{algorithm}
\end{mdframed}


\subsection{Approximating the Spectrum in EMD}\label{sec: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 \cite{LeeGT12, LouisRTV12, KwokLLGT13, CzumajPS15, gluch2021spectral}, 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 \cite{Cohen-SteinerKS18}, ``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.

\begin{definition}[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. 
\begin{equation}\label{eq:emd}
    \textrm{EMD}(A,B) = \min_{\pi : A \to B} \sum_{a \in A} \norm{a - \pi(a) }_2,
\end{equation}
where $\pi$ ranges over all one-to-one mappings. 
\end{definition}



We can now invoke the algorithm \texttt{ApproxSpectralMoment} of \cite{Cohen-SteinerKS18}. 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 \ref{thm:random_walk} 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 \ref{alg:neighbor_sampling} to perfectly sample from the true neighbor distribution. Thus we immediately have the following guarantee:


\begin{theorem}[Corollary of Theorem 1 in \cite{Cohen-SteinerKS18} and Theorem \ref{thm:random_walk}]\label{thm:emd}
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} K D^{-1}$ be the normalized Laplacian, where $D_{i,i} = \sum_{j} K_{i,j}$.  Let $\lambda_1 \geq \lambda_2 \ldots \geq \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) \leq \varepsilon.\]
\end{theorem}

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


\subsection{First Eigenvalue and Eigenvector Approximation}\label{sec:eigenvalue}

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 \cite{backurs2021} on the \emph{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.





\begin{mdframed}
  \begin{algorithm}[First Eigenvalue and Eigenvector Approximation]
  \label{alg:first_eig}\mbox{}
    \begin{description}
    \item[Input:]  Input dataset $X \subset \R^d$ of size $|X| = n$, precision $\eps > 0$. 
    \item[Operation:]\mbox{}
    \begin{enumerate}
    \item Let $t \gets O(1/(\eps^2 \tau^2)$. Let $S \gets$ random subset of $[n]$ of size $t$. Let $X_S$ be the samples restricted to the indices in $S$.  
    \item Let $K_S \gets$ principal submatrix of $K$ on indices in $S$ and let $\tilde{K} \gets (n/s) \cdot K_S$.  \texttt{// Just for notation; we do not initialize $K$ or $K_S$}
    \item Construct a KDE data structure for $X_S$. Run Algorithm $1$ of \cite{backurs2021} (Kernel Noisy Power Method) on $K_S$. Let $\hat{\lambda}_{\max}$ be the resulting 
    eigenvalue and $\hat{v}_{\max}$ be the resulting eigenvector.  
    \end{enumerate}
     \item[Output:] $\hat{\lambda}_{\max}$ and $\hat{v}_{\max}$.
    \end{description}
  \end{algorithm}
\end{mdframed}

We remark that the eigenvector returned by Algorithm \ref{alg:first_eig} will be a \emph{sparse} vector supported only on the coordinates in $S$.

We first state the necessary auxiliary statements needed to prove the guarantees of Algorithm \ref{alg:first_eig}. 

\begin{lemma}\label{lem:first_eig_large}
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 \ge n\tau$.
\end{lemma}
\begin{proof}
This follows from looking at the quadratic form $\textbf{1}^TK\textbf{1}$ where $\textbf{1}$ is the vector with all entries equal to $1$:
\[\lambda_1 \ge \frac{\textbf{1}^TK\textbf{1}}{\textbf{1}^T\textbf{1}}\ge \frac{n^2 \tau}{n} = n\tau. \qedhere \]
\end{proof}

We now state the guarantees of Algorithm $1$ in \cite{backurs2021}.

\begin{theorem}[\cite{backurs2021}]\label{thm:noisy_pm}
Suppose the kernel function for a $m \times m$ kernel matrix $K$ has a KDE data structure with query time $d/(\eps^2 \tau^p)$ (see Table \ref{table:kde_instantiation}). Then Algorithm $1$ of \cite{backurs2021} returns $\lambda$ such that $\lambda \ge (1-\eps) \lambda_1(K)$ in time $O\left( \frac{d m^{1+p}\log(m/\eps)^{2+p}}{\eps^{7+4p}}\right)$
\end{theorem}

Finally, we need the following result on eigenvalues of sampled PSD matrices, proven in \cite{bhattacharjee21}.

\begin{lemma}[\cite{bhattacharjee21}]\label{lem:psd_eig_sampling}
Let $A \in \R^{n \times n}$ be PSD with $\norm{A}_{\infty} \leq 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}\Paren{A_{S\times S} } = \lambda_i(A) \pm \frac{n}{\sqrt{t}}$. 
\end{lemma}

We are now ready to prove the guarantees of Algorithm \ref{alg:first_eig}.

\begin{theorem}\label{thm:first_eig}
Given a $n \times n$ kernel matrix $K$ admitting a KDE data-structure with query time $d/(\eps^2 \tau^p)$, Algorithm \ref{alg:first_eig} returns $\lambda$ such that $\lambda \ge (1-\eps)\lambda_1(K)$ in total time 
\[
\min \Paren{ O\left( \frac{d \log(d/\eps)}{\eps^{4.5} \tau^4} \right) ,  O\left( \frac{d}{\eps^{9+6p} \tau^{2+2p}} \log\left(\frac{1}{\eps \tau} \right)^{2+p} \right) }. \]
\end{theorem}

\begin{remark}
\normalfont
Two remarks are in order. First we recall that the runtime of \cite{backurs2021} has a $n^{1+p}$ factor while our bound \emph{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 \cite{MuscoM15}), we would get a runtime of $\tilde{O}(d/(\eps^{4.5}\tau^4))$ which has a polynomially better $\eps$ dependence but a worse $\tau$ dependence than the guarantees of Algorithm \ref{alg:first_eig}. 
\end{remark}

\begin{proof}[Proof of Theorem \ref{thm:first_eig}]
We first prove the approximation guarantee. By our setting of $t$ and using Lemma \ref{lem:psd_eig_sampling}, 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}} \leq \varepsilon \tau n \leq \varepsilon \lambda_1(K), \]
and thus $\lambda_1(\tilde{K}) \ge (1-\eps) \lambda_1(K)$. Then by the guarantees of Theorem \ref{thm:noisy_pm}, it follows that we find a $1-\eps$ multiplicative approximation to $\lambda_1(\tilde{K})$ and thus a $1-O(\eps)$ multiplicative approximation to that of $\lambda_1(K)$. 

We now prove the runtime bound. It easily follows from plugging in $m = O(1/(\eps^2 \tau^2))$ in Theorem \ref{thm:noisy_pm}.
\end{proof}

\section{Graph Applications}\label{sec:graph_applications}
In this section, we present our graph applications, including local clustering, spectral clustering, arboricity estimation, and estimating the total weight of triangles. 

\subsection{Local Clustering}\label{sec:local_clustering}



\begin{mdframed}
  \begin{algorithm}[Local $k$-Clustering]
  \label{alg:local_clustering}\mbox{}
    \begin{description}
    \item[Input:]  Input dataset $X \subset \R^d$ of size $|X| = n$, vertices $u,w$, random walk length $t$.  
    \item[Operation:]\mbox{}
    \begin{enumerate}
    \item For a given $v$, let $p_v^{t}$ be the endpoint distribution of a random walk of length $t$ starting at $v$.
    \end{enumerate}
     \item[Output:] ``$u,w$ are in the same cluster'' if $\ell_2$ distribution tester (see Theorem~\ref{thm:l2_tester}) outputs $\|p_u^{t} - p_w^{t}\|_2 \le 1/(7n) $. Otherwise, output ``$u,w$ are in different clusters''.
    \end{description}
  \end{algorithm}
\end{mdframed}

We give a local clustering algorithm on graphs. The advantage of this method is that it is \emph{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 \cite{kale2008testing, czumaj2010testing, goldreich2011testing, CzumajPS15, chiplunkar2018testing, dey2019spectral, gluch2021spectral} and the references within).

First, we need to define the notion of conductance.

\begin{definition}[Conductance]\label{def: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).\]
\end{definition}

\begin{definition}[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$.
\end{definition}

\begin{definition}[$k$-clusterable Graph]\label{def:k_clusterable_graph2}
A graph $G$ is $(k, \phi_{in}, \phi_{out})$-clusterable if the following holds: There exists a partition of the vertex set into $h \le k$ parts $V = \cup_{1 \le i \le h} V_i$ such that $\phi(G[V_i]) \ge \phi_{in}$ and $\phi_G(V_i) \le \phi_{out}$.
\end{definition}

Definition \ref{def:k_clusterable_graph2} captures the intuition that one can partition the graph into $h \le 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| \ge n/\poly(k)$ where we allow for an arbitrary polynomial dependence on $k$. This means that each cluster size is not too small.

Since we are interested in clustering, through this section, we will assume our kernel graph $K$ is $k$-clusterable according to Definition \ref{def:k_clusterable_graph2} 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 \ref{def:k_clusterable_graph2}), we can efficiently test if $V_1 = V_2$ or $V_1 \ne 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 \ref{alg:local_clustering}. The flavor of the algorithm presented is quite standard in property testing literature, see \cite{CzumajPS15} and \cite{Peng20}. 

The $\ell_2$ distribution tester we need is a standard result in distribution testing with the following guarantees.

\begin{theorem}[Theorem $1.2$ in \cite{ChanDVV14}]\label{thm:l2_tester}
Let $\delta, \xi > 0$ and let $p,q$ be two discrete distributions over a set of size $n$ with $b \ge \max\{ \|p\|_2^2, \|q\|_2^2\}.$ 
Let $r \ge 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 \le \xi$, and rejects the distributions if $\|p-q\|_2^2 \ge 4\xi$ with probability at least $1-\delta$. The running time of the tester is linear in its sample size.
\end{theorem}


We now prove the correctness of Algorithm \ref{alg:local_clustering}.
 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 \cite{CzumajPS15} but we must show it holds for weighted graphs arising from kernel matrices whereas \cite{CzumajPS15} is interested in bounded degree unweighted 
graphs.

\begin{lemma}
Let $\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]) \ge \phi_{in}$. Then for any $1 \le i \le h$, the following holds:
\[ \sum_{u,v \in C}\left(\frac{\V_i(u)}{\sqrt{w(u)}} - \frac{\V_i(v)}{\sqrt{w(v)}}\right)^2 \lesssim \frac{\phi_{out}n}{\phi_{in}^2|C| \tau^2}. \]
\end{lemma}
\begin{proof}
By Lemma $5.2$ in \cite{CzumajPS15} and Theorem $1.2$ in \cite{LeeGT12}, we have that $\phi_{in}^2/h^4 \lesssim \lambda_{h+1}$ and $\lambda_i \le 2 \phi_{out}$ for any $1 \le i \le h$. Now by the variational principle for eigenvalues \cite{chung1997spectral}, we have
\[ \lambda_i =\sum_{(u,v)} \left(\frac{\V_i(u)}{\sqrt{w(u)}} - \frac{\V_i(v)}{\sqrt{w(v)}} \right)^2 w(u,v) \le 2 \phi_{out}.\]
Now let $H = K[C]$. From \cite{chung1997spectral} 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{\V_i(u)}{\sqrt{w(u)}} - \frac{\V_i(v)}{\sqrt{w(v)}}\right)^2w(u,v)}{\sum_{u,v \in V_H} \left(\frac{\V_i(u)}{\sqrt{w(u)}} - \frac{\V_i(v)}{\sqrt{w(v)}}\right)^2d_H(u)d_H(v)} \ge \lambda_2(H) \ge \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{\V_i(u)}{\sqrt{w(u)}} - \frac{\V_i(v)}{\sqrt{w(v)}} \right)^2 \le \sum_{(u,v)\in E_K} \left(\frac{\V_i(u)}{\sqrt{w(u)}} - \frac{\V_i(v)}{\sqrt{w(v)}} \right)^2 \le 2 \phi_{out}. \]
This implies that 
\begin{align*}
     |C|^2 \tau^2 \sum_{(u,v) \in V_H} \left(\frac{\V_i(u)}{\sqrt{w(u)}} - \frac{\V_i(v)}{\sqrt{w(v)}} \right)^2  &\le \sum_{(u,v) \in V_H} \left(\frac{\V_i(u)}{\sqrt{w(u)}} - \frac{\V_i(v)}{\sqrt{w(v)}} \right)^2  d_H(u)d_H(v) \\
     &\lesssim \frac{\phi_{out} \text{vol}_H(V_H)}{\phi_{in}^2}
\end{align*}
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) \le |C|n$, it follows that 
\[ \sum_{(u,v) \in V_H} \left(\frac{\V_i(u)}{\sqrt{w(u)}} - \frac{\V_i(v)}{\sqrt{w(v)}} \right)^2   \lesssim \frac{\phi_{out}n}{\phi_{in}^2|C| \tau^2}, \]
as desired.
\end{proof}

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 \cite{CzumajPS15} but we must show it holds for weighted graphs.

\begin{lemma}\label{lem:local_same_cluster}
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| \ge n/\poly(k)$ and $\phi(K[C]) \ge \phi_{in}$. There exists a constant $c = c(\beta) > 0$ and $c' = c'(\beta, k)$ such that for any $t \ge c \log n/\phi_{in}^2, \phi_{out} \le c' \phi_{in}^2$, there exists a subset $\widetilde{C} \subseteq C$ satisfying
$\text{vol}(\widetilde{C}) \ge (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 \le \frac{1}{8n} .\]
\end{lemma}
\begin{proof}
  Let $\V_1, \cdots, \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 $\y_i = D^{1/2}\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_iy_i = \sum_i \alpha_i D^{-1/2}\V_i.\] Taking the innerproduct of $1_u$ with $D^{-1/2}\V_i$ tells us that $\alpha_i = \V_i(u)/\sqrt{w(u)}$. Thus,
\begin{align*}
    p_u^t - p_v^t &= \sum_{i=1}^n \frac{\V_i(u)}{\sqrt{w(u)}} (1-\lambda_i)^t \y_i - \sum_{i=1}^n \frac{\V_i(v)}{\sqrt{w(v)}} (1-\lambda_i)^t \y_i \\
    &= D^{1/2} \sum_{i=1}^n \V_i \left( \frac{\V_i(u)}{\sqrt{w(u)}} -\frac{\V_i(v)}{\sqrt{w(v)}} \right)(1-\lambda_i)^t.
\end{align*}
This means that 
\[ \|p_u^t - p_v^t  \|_2 \le \|D^{1/2}\| \left\| \sum_{i=1}^n \V_i \left( \frac{\V_i(u)}{\sqrt{w(u)}} -\frac{\V_i(v)}{\sqrt{w(v)}} \right)(1-\lambda_i)^t \right\|_2. \]
Since the $\V_i$'s are orthogonal, we know that 
\[ \left\| \sum_{i=1}^n \V_i \left( \frac{\V_i(u)}{\sqrt{w(u)}} -\frac{\V_i(v)}{\sqrt{w(v)}} \right)(1-\lambda_i)^t \right\|_2^2 \le \sum_{i=1}^n \left( \frac{\V_i(u)}{\sqrt{w(u)}} -\frac{\V_i(v)}{\sqrt{w(v)}} \right)^2(1-\lambda_i)^{2t}. \]
Now the rest of the proof follows from Lemma  $4.2$ in \cite{CzumajPS15}. In particular, it tells us that in the above summation, each of the terms for $1 \le i \le 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/\poly(n)$ for sufficiently large $\poly(n)$ by adjusting the constant in front of $t$. Our choice for $|C| \ge n/\poly(k)$ imply that the overall sum is bounded by $\lesssim \frac{\phi_{out}\poly(k) }{ \beta n^2 \phi_{in}^2}$. Since $\|D^{1/2}\|^2 \le n$, and $\phi_{out} \le c' \phi_{in}^2$, we have that $\|p_u^t - p_v^t \|_2^2 \le 1/(8n)$, as desired.
\end{proof}

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 \cite{czumaj2010testing}.

\begin{lemma}\label{lem:local_not_escape}
Let $G$ be a $(k, \phi_{in}, \phi_{out})$-clusterable graph with parts $V = \cup_{1 \le i \le h} V_i$. There exists a constant $c > 0$ such that if $t \cdot \phi_{out} \le c\eps$, then there exists a subset $V_1' \subseteq V_1$ satisfying $\text{vol}(V_1') \ge (1-\eps)\text{vol}(V_1)$ such that a $t$-step random walk from any vertex in $V_1'$ does not leave $V_1$ with probability $1-\eps$.
\end{lemma}
\begin{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'$ 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'$ of length $t$ with the initial vertex $u'$ chosen from the stationary distribution of $G'$, i.e., the distribution that chooses each vertex in $G'$ 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 
\[\Pr[Y_i = 1] =  \frac{w(V_1, V\setminus V_1)}{w(G')} \]
where $w$ is the weight of edges in the original graph $K$. By linearity of expectations, the number of vertices that land in $V \setminus V_1$ is 
\[\E\left[\sum_{i=1}^{t} Y_i\right] = (t+1) \frac{w(V_1, V\setminus V_1)}{w(G')}  \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 \eps/t$, then we know that the random walk stays in $V_1$ with probability at least $1-O(\eps)$ so there must be a set of vertices $V' \subseteq V_1$ of at least $1-O(\eps)$ fraction of the total volume of $V_1$ such that a random walk starting from a vertex in $V'$ remains in $V_1$ with probability at least $1-O(\eps)$.
\end{proof}

We can now prove the correctness of Algorithm \ref{alg:local_clustering}.

\begin{theorem}\label{thm:localalg_correctness}
Let $K$ be a $(k, \phi_{in}, \phi_{out})$-clusterable kernel graph with parts $V = \cup_{1 \le i \le 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(\eps, k)$ such that if $\phi_{out} \le c\phi_{in}^2/\log n$, then with probability at least $1-\eps$, if $U = W$ then Algorithm \ref{alg:local_clustering} returns that $u$ and $w$ are in the same cluster and if $U \ne W$, Algorithm \ref{alg:local_clustering} returns that $u$ and $w$ are in different clusters. The algorithm requires $O(\sqrt{nk/(\eps \tau)} \log (1/\eps))$ random walks of length $t \ge c \log n/\phi_{in}^2$. \end{theorem}

\begin{proof}
We first consider the case that $U \ne W$.
From Lemma \ref{lem:local_not_escape}, we know that there are `non-escaping' subsets $U'$ and $W'$ of $U$ and $W$ respectively such that vertices $u,w$ from $U'$ and $W'$ respectively don't leave $U$ and $W$ with probability $1-\eps$. Conditioning on $u$ and $w$ being in those subsets, we have that with probability $1-O(\eps)$, $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  \ge 2/n$.

Now if $U = W$, we know from Lemma \ref{lem:local_same_cluster} that $\|p_u^t - p_v^t\|_2^2 \le 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 \cite{CzumajPS15} readily implies that there exists a $V' \subseteq V$ satisfying $\text{vol}(V') \ge (1-\eps) \text{vol}(V)$ such that $\|p_u^t\|_2^2 \le 2k/(\eps \tau^2 n)$. Now we can set $\xi = 1/(7n)$ and $b = 2k/(\eps \tau^2 n)$ in Theorem \ref{thm:l2_tester}, which tells us that $r = O(\sqrt{nk/(\eps \tau)} \log (1/\eps))$ samples of the distributions $p_u^t$ and $p_w^t$ suffice to distinguish the cases $\|p_u^t - p_v^t\|_2^2  \ge 2/n$ or $\|p_u^t - p_v^t\|_2^2 \le 1/(8n)$, i.e., $r$ samples allow us to determine if $U = W$ or $U \ne W$, conditioned on a $1-O(\eps)$ probability event.
\end{proof}

It is straightforward to translate the requirements of Theorem \ref{thm:localalg_correctness} 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 \ref{thm:random_walk}. Alternatively, we can perform rejection sampling as stated in the proof of Theorem \ref{thm:neighbor_sampling}.

\begin{corollary}
Algorithm \ref{alg:local_clustering} and Theorem \ref{thm:localalg_correctness} require $\widetilde{O}(c(k, \eps) \sqrt{nk/\eps} \cdot 1/( \tau^{1.5} \phi_{in}^2))$ KDE queries (via calls to Algorithm \ref{alg:random_walk}, which performs random walks) as well as the same bound for post-processing time.
\end{corollary}

\subsection{Spectral Clustering}\label{sec: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 \cite{von2007tutorial} 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 \ref{def:k_clusterable_graph2}.

\begin{definition}[Weakly clusterable Graph]\label{def:k_clusterable_graph}
A graph is $(k,  \phi_{out})$-clusterable if the following holds: There exists a partition of the vertex set into $h \le k$ parts $V = \cup_{1 \le i \le h} V_i$ such that $\phi_G(V_i) \le \phi_{out}$.
\end{definition}


We now prove the following result that says spectral sparsification preserves cluster structure according to Definition \ref{def:k_clusterable_graph}. We first remark that the spectral sparsifier obtained in the previous section is a \emph{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 \eps$. The implication follows immediately by noting that cuts are induced by quadratic forms on the Laplacian matrix using $\{-1,1\}^n$ vectors.

\begin{theorem}\label{thm:spectral_clustering}
Let $G$ be $(k, \phi_{out})$-clusterable and let $G'$ be a cut sparsifier for $G$. Then $G'$ is $(k, (1\pm \eps)\phi_{out})$-clusterable.
\end{theorem}
\begin{proof}
Let $V_i$ be one of the $h \le k$ vertex partitions of $G$. Consider the conductance of $V_i$ defined in Definition \ref{def:conductance}. 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'$ is a cut sparsifier, this implies that both the numerator and denominator are preserved up to a $1\pm \eps$ factor and thus, the entire ratio is also preserved up to a $1\pm O(\eps)$ factor.
\end{proof}

Theorem \ref{thm:spectral_clustering} implies that the cluster structure of the sparsified graph $G'$ is approximately identical to that of $G$. Thus, we can be confident that the spectral clustering procedure described at the beginning of the section would perform equally as well on $G'$ as it would have on $G$. Indeed, we verify this empirically in Section \ref{sec:experiments}.

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 \cite{MuscoM15}, which says (a variant of) the power method can quickly find good approximations Laplacian eigenvectors if the matrix is sparse.

\begin{theorem}[Corollary of Theorem $1$ in \cite{MuscoM15} and Theorem \ref{thm:spec-spar-kde-matrix}]\label{thm:laplacian_eigv}
Let $L$ be the Laplacian matrix of the sparsifier computed in Theorem \ref{thm:spec-spar-kde-matrix}. Let $u_1 , \cdots, u_k$ be the first $k$ eigenvectors of $L$. Using Theorem $1$ of \cite{MuscoM15}, we can find $k$ vectors $v_1, \cdots, v_k$ in 
time $\widetilde{O}\left( \frac{kn \log n}{\tau^2 \eps^{2.5}}\right)$ such that with probability $99/100$,
\[|u_i^TLu_i - v_i^TLv_i| \le \eps\lambda_{k+1}^2 \]
for all $i \in [k]$.
\end{theorem}

\subsection{Arboricity Estimation}\label{sec:arboricity}


\begin{mdframed}
  \begin{algorithm}[Arboricity Estimation]
  \label{alg:arboricity:est}\mbox{}
    \begin{description}
    \item[Input:] Input dataset $X \subset \R^d$ of size $|X| = n$, weight function $w(\cdot)$, accuracy parameter $\eps$. 
    \item[Operation:]\mbox{}
    \begin{enumerate}
    \item Let $\Delta=\max_{e,e'\in E}\frac{w(e)}{w(e')}$, $m= O\left(\frac{n\Delta\log n}{\eps^2}\right)$, and $G'=\emptyset$.
    \item For $i=1$ to $i=m$:
    \begin{enumerate}
        \item 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]$
        \item Add $e$ to $G'$ with weight $\frac{1}{mp_e}$
    \end{enumerate}
    \end{enumerate}
     \item[Output:] $\max_{U\subseteq V} d(G'_U)$.
    \end{description}
  \end{algorithm}
\end{mdframed}

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 \cite{LeeRJA10} 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 \ref{sec:blocks}. 
Inspired by the unweighted version of the arboricity estimation algorithm from \cite{McGregorTVV15}, we first prove the following result.


\begin{theorem}\label{thm:arboricity}
Let $U' = \argmax_U d(G_U')$ and let $G'$ be the output of Algorithm~\ref{alg:arboricity:est}. 
Then with probability at least $1-1/\poly(n)$,
\[(1 - \eps)\alpha \le d(G'_{U'}) \le (1+\eps)\alpha. \]
Algorithm~\ref{alg:arboricity:est} uses $m = \widetilde{O}(n \log n/(\eps^2 \tau))$ KDE queries and $O(mn)$ post-processing time.
\end{theorem}
\begin{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\ge\frac{W}{n}$, so that 
\[m\ge\frac{\log n}{\alpha\Delta\eps^2}.\]

Let $X_i$ be the random variable denoting the contribution of the $i$-th sample to weight of the edges in $G'_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_ew_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}{\eps^2}$ for an absolute constant $C>0$, then 
\[\mathbb{E}[X_i^2]\le\frac{k\eps^2\alpha^2}{Cm\log^2 n}\]
and
\[\sum_{i=1}^m\mathbb{E}[X_i^2]\le\frac{k\eps^2\alpha^2}{Cm\log^2 n}.\]
We also have $X_i-\mathbb{E}[X_i]\le\frac{\alpha\eps^2}{Cm\log^2 n}$. 
Thus by Bernstein's inequality for sufficiently large $C$,
\[\Pr\left[d(G'_U)\ge\frac{\alpha}{10}\right]\le n^{-10k},\]
for $d(G_U)\le\frac{\alpha}{60}$ and
\[\Pr\left[|d(G'_U)-d(G_U)|\ge\frac{\eps\alpha}{10}\right]\le 2n^{-10k},\]
for $d(G_U)>\frac{\alpha}{60}$. 

Since there are $\binom{n}{k}\le 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'_U)\le\frac{\alpha}{10}\]
for all subsets $U$ with $d(G_U)\le\frac{\alpha}{60}$ and 
\[(1-\eps)d(G_U)\le d(G'_U)\le(1+\eps)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^*})\ge(1-\eps)\alpha$ so that $d(G_{U'})\ge d(G_{U^*})\ge(1-\eps)\alpha$, where $U'=\argmax_{U\subseteq V} d(G'_U)$. 
Thus with high probability, we have that
\[(1-\eps)\alpha\le d(G_{U'})\le(1+\eps)\alpha,\]
as desired.
\end{proof}

To estimate the arboricity of the input graph $G$, it then suffices Theorem~\ref{thm:arboricity} to compute the arboricity of the subsampled graph $G'$ output by Algorithm~\ref{alg:arboricity:est}. 
This can be efficiently achieved by running an offline algorithm such as~\cite{Charikar00}, 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. 



\subsection{Computing the Total Weight of Triangles}\label{sec: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 \cite{S13, AGD08, BC17, Kolountzakis_2010, ChenEILNRSWWZ22}. Applications include discovering motifs in protein interaction networks \cite{Milo824}, understanding social networks \cite{Welles10}, and evaluating large graph models \cite{Leskovec08}; see the survey \cite{al2018triangle} 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 \cite{kalna2006clustering, li2007structure, antoniou2008statistical}.

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.


\begin{definition}
Let $G = (V, E, w)$ with $w : E \rightarrow \R^{\ge 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))$.
\end{definition}

For this definition, we present the following modified algorithm from \cite{EdenLRS17}, which considers the problem in unweighted graphs in a different model of computing. See Remark \ref{remark:triangles} for comparison.

\begin{theorem}\label{thm:triangles}
There exists an algorithm that makes $\widetilde{O}\left(\frac{m\sqrt{w_G}w_{\max}^{3/2}}{w_T \cdot \eps^2}\right)$ KDE queries and the same bound for post-processing time and with probability at least $\frac{2}{3}$, outputs a $(1\pm\eps)$-approximation to the total weight $w_T$ of the triangles in the kernel graph. 
\end{theorem}
\begin{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 have an ordering on the vertex set $V$ by $u\prec v$ if and only if either $w(u)\le w(v)$ or $w(u)=w(v)$ and $u$ appears $v$ in the dictionary ordering of the vertices. 
For each edge $e=(u,v)$, we assign to $e$ all triangle $(u,v,w)$ such that $u\prec v\prec w$. 
Let $W_e$ denote the weight of the triangles assigned to $e$. 

Suppose, by way of contradiction, there exists $e\in E$ with $W_e>\sqrt{w_G}w_{\max}^{3/2}$. 
Since each triangle can contribute at most $w_{\max}^3$ weight to $W_e$, then more than $\sqrt{w_G}w^{-3/2}$ triangles must be assigned to $e$. 
Then there must be more than $\sqrt{w_G}w^{-3/2}$ vertices with weight at least $\sqrt{w_G}w_{\max}^{3/2}$, which contradicts the fact that the graph $G$ has weight $w_G$. 
Thus, we have that $W_e\le\sqrt{w_G}w_{\max}^{3/2}$ for all $e\in E$. 
Moreover, we have that $\sum_{e\in E}W_e=w_T$ so that
\[\mathbb{E}_{e\in E}[W_e]=\frac{w_T}{m}\]
and
\[\mathbb{E}_{e\in E}[W_e^2]\le\frac{w_T\sqrt{w_G}w_{\max}^{3/2}}{m}.\]
By Chebyshev's inequality, it suffices to sample a set $R$ with \[|R|=O\left(\frac{m\sqrt{w_G}w_{\max}^{3/2}}{\eps^2 w_T}\right)\]
edges uniformly at random, so that
\[\Pr\left[\sum_{e\in R}W_e\in (1\pm \eps)|R|\cdot\frac{w_T}{m}\right]\ge 0.99.\]

For each vertex $v\in V$, let $g(v)=\sum_{(u,v)\in E} w(u,v)$ and for each edge $e=(u,v)$, let $g(e)=\min(w(u),w(v))$. 
We write $g(R)=\sum_{e\in R}g(e)$. 
Now for each $i\in[|R|]$, we have
\[\mathbb{E}[\chi_i]=\sum_{e\in R}\frac{W_e}{g(R)}=\frac{W_R}{g(R)},\qquad\mathbb{E}[\chi_i^2]\le 1.\]
Hence by Bernstein bounds, there exists a constant $C>0$ such that it suffices to repeat the procedure $\frac{C}{\eps^2}\cdot\frac{\cdot g(R)}{W_R}$ times to get a $(1\pm\eps)$-approximation of $\frac{W_R}{g(R)}$ with probability at least $2/3$. 
We have that $W_R\ge\frac{w_T}{2m}\cdot|R|$ and $\mathbb{E}[g(R)]\le\sqrt{w_G}w_{\max}^{3/2}\cdot|R|$ so that 
\[\frac{g(R)}{W_R}\le\frac{2m\sqrt{w_G}w_{\max}^{3/2}}{w_T}. \qedhere\]
\end{proof}


\section{Empirical Evaluation}\label{sec:experiments}
We present empirical evaluations for our algorithms. We chose to evaluate algorithms for low-rank approximation and spectral sparsification (and spectral clustering as a corollary) as they are arguably two of the most well studied examples in our applications and utilize a wide variety of techniques present in our other examples of Sections \ref{sec:linalgebra_applications} and \ref{sec:graph_applications}. For our experiments, we use the Laplacian kernel $k(x,y) = \exp(-\|x-y\|_1/\sigma)$. A fast KDE implementation of this kernel exists due to \cite{backurs2019space}, which builds upon the techniques of \cite{charikar2017hashing}. Note that the focus of our work is to use KDE queries in a mostly black box fashion to solve important algorithmic problems for kernel matrices. This viewpoint has the important advantage that it is flexible to the choice of any particular KDE query instantiation. We chose to work with the implementation of \cite{backurs2019space} since it possesses theoretical guarantees, has an accessible implementation\footnote{from \url{https://github.com/talwagner/efficient_kde}}, and has been used in experiments in prior works such as \cite{backurs2019space, backurs2021}. However, we envision other choices of KDE queries, which maybe have practical benefits but are theoretically incomparable would also work well due to our flexibility.


\begin{figure*}[ht]
\centering
\begin{subfigure}{.5\textwidth}
  \centering
 \includegraphics[width=0.95\linewidth]{images/nested_dataset.pdf}  \caption{\label{fig:nested_dataset} }
\end{subfigure}%
\begin{subfigure}{.5\textwidth}
  \centering
 \includegraphics[width=0.95\linewidth]{images/rings_dataset.pdf} 
 \caption{\label{fig:rings_dataset} }
\end{subfigure}
\caption{(a) Nested Dataset, (b) Rings Dataset}
\label{fig:datasets}
\end{figure*}


\paragraph{Datasets.}
We use two real and two synthetic datasets in our experiments. The datasets used in the low-rank approximation experiments are MNIST (points in $\R^{784}$) \cite{lecun1998mnist} and Glove word embeddings (points in $\R^{200}$)\cite{pennington2014glove}. We use $10^4$ points from each of their test datasets. These datasets have been used in prior experimental works on kernel density estimation \cite{siminelakis2019rehashing,backurs2019space}.

For spectral sparsification and clustering, we use construct two synthetic datasets, which are challenging for other clustering method such as $k$-means clustering\footnote{For example, see \url{https://scikit-learn.org/stable/auto_examples/cluster/plot_cluster_comparison.html}.}. The first dataset denoted as `Nested' consists of $5,000$ points, equally split among the origin and a circle of radius $1$. The two natural clusters are the points at the origin and the points on the circle. Since one cluster is contained in the convex hull of the other, a method like $k$-means clustering will not be able to separate the two clusters, but it is known that spectral clustering can. Our second dataset, labeled `Rings', is an even more challenging clustering dataset. We consider two tori in three dimensions that pass through the interior hole of each other, i.e., they interlock. The `small' radius of each tori is $5$ while the `large' radius is $100$. Our dataset consists of $2500$ points uniformly distributed on the two tori; see Figure \ref{fig:rings_dataset}. Note that our focus is not to compare the efficacy of various clustering methods, which is done in other prior works (e.g., see footnote $2$). Rather, we show that spectral clustering itself can  be optimized in terms of runtime, space usage, and the number of kernel evaluations performed via our algorithms.

\paragraph{Evaluation metrics.}
For low-rank approximation, we use the additive error algorithm detailed in Corollary \ref{cor:additive_lra} of Section \ref{sec:lra}. It requires sampling the rows of the kernel matrix according to squared row norms, which can be done via KDE queries as outlined there. Once the (small) number of rows are sampled, we explicitly construct these rows using kernel evaluations. We compare the approximation error of this method computed via the standard Frobenius norm error  to a state of the art sketching algorithm for computing low-rank approximations, which is the input-sparsity time algorithm of Clarkson and Woodruff \cite{ClarksonW13} \textbf{(IS)}. We also compare to an iterative SVD solver \textbf{(SVD)}. All linear algebra subroutines rely on Numpy and Scipy implementations and Numba complication when applicable. 

Note that prior works such as \cite{backurs2019space, backurs2021} have used use the number of kernel evaluations performed (i.e., how many entries of $K$ are computed) as a measure of computational cost. While this is a software and architecture independent basis of comparison, which is unaffected by access to specialized libraries or hardware (e.g., SIMD, GPU), it is of interest to go beyond this measure. Indeed, we use this measure as well as other important metrics as space usage and runtime as points of comparison. For spectral sparsification and clustering, we compare the accuracy of our method to the clustering solution when run on the full initialized kernel matrix. 

\paragraph{Parameter settings.}
For low-rank approximation datasets, we choose the bandwidth value $\sigma$ according to the choice made in prior experiments, in particular the values given in \cite{backurs2019space}. There, $\sigma$ is chosen according to the popular median distance rule; see their experimental section for further information. For our clustering experiments, we pick the value of $\sigma$, which results in spectral clustering (running on the full kernel matrix) successfully clustering the input.


\subsection{Results}

\paragraph{Low-rank approximation.}
Note that the algorithm in Corollary \ref{cor:additive_lra} has a $O(k)$ dependence on the number of rows sampled. Concretely we sample $25k$ rows for a rank $k$ approximation which we fix it for all experiments. For the MNIST dataset, the rank versus approximation error is shown in Figure \ref{fig:lra_mnist}. The performance of our algorithm labeled as $\textbf{KDE}$ is given by the blue curve while the orange curve represents the $\textbf{IS}$ algorithm. The green curve represents the SVD error, which is a lower bound on the error for any algorithm. Note that for SVD calculations, we do not calculate the full SVD since that is computationally prohibitive; instead, we use an iterative solver. We can see that the errors of all three methods are comparable to each other.  In terms of runtime, the KDE based method took $24.7$ seconds on average for the rank $50$ approximation whereas $\textbf{IS}$ took $71.5$ seconds and iterative SVD took $74.72$ seconds on average. This represents a \textbf{2.9x} decrease in the running time. The time measured includes the time to initialize the data structures and matrices used for the respective algorithms. In terms of the number of kernel evaluations, both $\textbf{IS}$ and iterative SVD require the kernel matrix, which is $10^{8}$ kernel evaluations. On the other hand for the rank $50$ approximation, our method required only $1.1 \cdot 10^7$ kernel evaluations, which is a \textbf{9x} decrease in the number of evaluations. In terms of space, $\textbf{IS}$ and iterative SVD require $10^8$ floating point numbers stored due to initializing the full $10^4 \times 10^4$ matrix whereas our method only requires $10^4 \cdot 25 \cdot 50$ floating point numbers for the rank equal to $50$ case and smaller for other. This is a \textbf{8x} decrease in the space required. Lastly, we verify that we are indeed sampling from the correct distribution required by Corollary \ref{cor:additive_lra}. In Figure \ref{fig:mnist_dist}, we plot the points $(x_i, y_i)$ where $x_i$ is the row norm squared for the $i$th row of the kernel matrix $K$ and $y_i$ is the row norm squared computed in our approximation algorithm (see Algorithm \ref{alg:additive_lra}). As shown in Figure \ref{fig:mnist_dist}, the data points fall very close to the $y=x$ line indicating that our algorithm is indeed sampling from approximately the correct ideal distribution. 

The qualitatively similar results for the Glove dataset are given in Figures \ref{fig:lra_glove} and \ref{fig:glove_dist}. For the glove dataset, the average time taken by the three algorithms were $37.7, 37.7,$ and $44.2$ seconds respectively, indicating that \textbf{KDE} and \textbf{IS} were comparable in runtime whereas SVD took slightly longer. However, the number of kernel evaluations required by the latter two algorithms was significantly larger: for rank equal to $10$, our algorithm only required $2.6 \cdot 10^6$ kernel evaluations while the other methods both required $10^8$ due to initializing the matrix. Lastly, the space required by our algorithm was smaller by a factor of \textbf{$40$} since we only explicitly compute $25 \cdot 10$ rows for the rank $= 10$ case. For Glove, we only perform our experiments up to rank equal to $10$ since the iterative SVD failed to converge for higher ranks. While computing the full SVD would avoid the convergence issue, it would take significantly longer time in general. For example for the MNIST dataset, computing the full SVD of the kernel matrix took $552.9$ seconds, which is approximately an order of magnitude longer than any of the other methods.

\begin{figure*}[!ht] 
  \begin{subfigure}[b]{0.5\linewidth}
    \centering
    \includegraphics[width=0.75\linewidth]{images/LRA_mnist.pdf} 
    \caption{Rank versus Error for Low-rank Approximation} 
    \label{fig:lra_mnist} 
  \end{subfigure
  \begin{subfigure}[b]{0.5\linewidth}
    \centering
    \includegraphics[width=0.75\linewidth]{images/LRA_mnist_normapprox.pdf} 
    \caption{Real vs Approximate Row Norm Squared Values} 
    \label{fig:mnist_dist} 
  \end{subfigure} 
  \begin{subfigure}[b]{0.5\linewidth}
    \centering
    \includegraphics[width=0.75\linewidth]{images/LRA_glove.pdf} 
    \caption{Rank versus Error for Low-rank Approximation} 
    \label{fig:lra_glove} 
  \end{subfigure
  \begin{subfigure}[b]{0.5\linewidth}
    \centering
    \includegraphics[width=0.75\linewidth]{images/LRA_glove_normapprox.pdf} 
    \caption{Real vs Approximate Row Norm Sqaured Values} 
    \label{fig:glove_dist} 
  \end{subfigure} 
  \caption{Figures for low rank approximation experiments.}
  \label{fig:lra_experiments} 
\end{figure*}


\begin{figure*}[!ht]
\centering
\begin{subfigure}{.5\textwidth}
  \centering
\includegraphics[width=0.8\linewidth]{images/nested_spectral_approx.pdf} 
 \caption{\label{fig:nested_spectral} }
\end{subfigure}%
\begin{subfigure}{.5\textwidth}
  \centering
\includegraphics[width=0.8\linewidth]{images/rings_spectral_approx.pdf}
 \caption{\label{fig:rings_spectral} }
\end{subfigure}
\caption{Spectral embedding of sparsified graph for (a) Nested dataset and (b) Rings dataset, respectively.}
\label{fig:spectral_experiments}
\end{figure*}


\paragraph{Spectral sparsification and clustering.}
Our algorithm consists of running the spectral sparsification algorithm of Theorem \ref{thm:spec-spar-kde-matrix} (Algorithm \ref{alg:spectral_sparsification}) and computing the first two eigenvectors of the normalized Laplacian of the resulting sparse graph. We then run $k$-means clustering on the computed Laplacian embedding for $k=2$. 

As noted above, we use two datasets that pose challenges for traditional clustering methods such as $k$-means clustering. The Nested dataset is shown in Figure \ref{fig:nested_dataset}. We sampled $3 \cdot 10^5$ many edges, which is $2.5\%$ of total edges. Figure \ref{fig:nested_spectral} shows the Laplacian embedding of the sampled graph based on the first two eigenvectors. The colors of the red and blue points correspond to their cluster in Figure \ref{fig:nested_dataset} as identified by running $k$-means clustering on the Laplacian embedding. The orange crosses are the points that the spectral clustering method failed to correctly classify. These are only $23$ points, which represent a $0.5 \%$ of total points. Furthermore, Figure \ref{fig:nested_spectral} shows that the Laplacian embedding of the sampled graph is able to embed the two clusters into distinct and disjoint regions. Note that the total space savings of the sampled graph over storing the entire graph is $\textbf{41x}$. In terms of the time taken, the iterative SVD method used to calculate the Laplacian eigenvectors took $0.18$ seconds on the sparse graph whereas the same method took $0.81$ seconds on the entire graph. This is a $\textbf{4.5x}$ factor reduction.

We recorded qualitatively similar results for the rings dataset. Figure \ref{fig:rings_dataset} shows a plot of the dataset. We sampled $10^5$ many edges for the approximation, which represents a $3.3\%$ of total edges for form the sparse graph. The Laplacian embedding of the sparse graph is shown in Figure \ref{fig:rings_spectral}. In this case, the embedding constructed from the sparse graph was able to separate the two rings into disjoint regions perfectly. The time taken for computing the Laplacian eigenvectors for the sparse graph was $0.08$ seconds whereas it took $0.27$ seconds for the full dense matrix.


\section{Auxiliary Inequalities}
\begin{theorem}[Bernstein's inequality]
\label{thm:bernstein:conc}
Let $X_1,\ldots,X_n$ be independent random variables satisfying $\mathbb{E}[X_i^2]<\infty$ and $X_i\ge 0$ for all $i\in[n]$. 
Let $X=\sum_i X_i$ and $\gamma>0$. 
Then
\[\Pr\left[X\le\mathbb{E}[X]-\gamma\right]\le\exp\left(\frac{-\gamma^2}{2\sum_i\mathbb{E}[X_i^2]}\right).\]
If $X_i-\mathbb{E}[X_i]\le\Delta$ for all $i$, then for $\sigma_i^2=\mathbb{E}[X_i^2]-\mathbb{E}[X_i]^2$,
\[\Pr\left[X\ge\mathbb{E}[X]+\gamma\right]\le\exp\left(\frac{-\gamma^2}{2\sum_i\sigma_i^2+2\gamma\Delta/3}\right).\]
\end{theorem}

\section*{Acknowledgements}
Piotr Indyk and Sandeep Silwal were supported by an NSF Graduate Research Fellowship under Grant No.\ 1745302, and NSF TRIPODS program (award DMS-2022448), NSF award CCF-2006798, and Simons Investigator Award. 
Samson Zhou was supported by a Simons Investigator Award and by the National Science Foundation under Grant No. CCF-1815840. Praneeth Kacham was supported by National Institute of Health (NIH) grant 5401 HG 10798-2, a Simons Investigator Award and Google as part of the ``Research Collabs'' program.

\bibliographystyle{alpha}
