Algorithms and Hardness for Linear Algebra on Geometric Graphs

Introduction

Linear algebra has a myriad of applications throughout computer science and physics. Consider the following seemingly unrelated tasks:

Each of these tasks has seen much work throughout numerical analysis, theoretical computer science, and machine learning. The first task is a celebrated application of the fast multipole method of Greengard and Rokhlin [GR87, GR88, GR89], voted one of the top ten algorithms of the twentieth century by the editors of Computing in Science and Engineering [DS00]. The second task is spectral clustering [NJW02, LWDH13], a popular algorithm for clustering data. The third task is to label a full set of data given only a small set of partial labels [Zhu05b, CSZ09, ZL05], which has seen increasing use in machine learning. One notable method for performing semi-supervised learning is the graph-based Laplacian regularizer method [LSZ+19b, ZL05, BNS06, Zhu05a].

$n$-body simulation (one step): For each $i\in\{1,2,\ldots,d\}$, make a weighted graph $G_{i}$ on the points in $X$, in which the weight of the edge between the points $u,v\in X$ in $G_{i}$ is $\mathsf{K}_{i}(u,v):=(\frac{G_{\text{grav}}\cdot m_{u}\cdot m_{v}}{\|u-v\|_{2}^{2}})(\frac{v_{i}-u_{i}}{\|u-v\|_{2}})$, where $G_{\text{grav}}$ is the gravitational constant and $m_{x}$ is the mass of the point $x\in X$. Let $A_{i}$ denote the weighted adjacency matrix of $G_{i}$. Then $A_{i}\textbf{1}$ is the vector of $i$th coordinates of force vectors. In particular, gravitational force can be computed by doing $O(d)$ adjacency matrix-vector multiplications, where each adjacency matrix is that of the $\mathsf{K}_{i}$-graph on $X$ for some $i$.

Spectral clustering: Make a $\mathsf{K}$ graph $G$ on $X$. In applications, $\mathsf{K}(u,v)=f(\|u-v\|_{2}^{2})$, where $f$ is often chosen to be $f(z)=e^{-z}$ [vL07, NJW02]. Instead of directly running a spectral clustering algorithm on $L_{G}$, one popular method is to construct a sparse matrix $M$ approximating $L_{G}$ and run spectral clustering on $M$ instead [CFH16, CSB+11, KMT12]. Standard sparsification methods in the literature are heuristical, and include the widely used Nystrom method which uniformly samples rows and columns from the original matrix [CJK+13].

If $H$ is a spectral sparsifier of $G$, it has been suggested that spectral clustering with the top $k$ eigenvectors of $L_{H}$ performs just as well in practice as spectral clustering with the top $k$ eigenvectors of $L_{G}$ [CFH16]. One justification is that since $H$ is a spectral sparsifier of $G$, the eigenvalues of $L_{H}$ are at most a constant factor larger than those of $L_{G}$, so cuts with similar conductance guarantees are produced. Moreover, spectral clustering using sparse matrices like $L_{H}$ is known to be faster than spectral clustering on dense matrices like $L_{G}$ [CFH16, CJK+13, KMT12].

In the first, second, and third tasks above, a small number of calls to matrix-vector multiplication, spectral sparsification, and Laplacian system solving, respectively, were made on geometric graphs. One could solve these problems by first explicitly writing down the graph $G$ and then using near-linear time algorithms [SS11, CKM+14] to multiply, sparsify, and solve systems. However, this requires a minimum of $\Omega(n^{2})$ time, as $G$ is a dense graph.

In this paper, we initiate a theoretical study of the geometric graphs for which efficient spectral graph theory is possible. In particular, we attempt to determine for which (a) functions $\mathsf{K}$ and (b) dimensions $d$ there is a much faster, $n^{1+o(1)}$-time algorithm for each of (c) multiplication, sparsification, and Laplacian solving. Before describing our results, we elaborate on the choices of (a), (b), and (c) that we consider in this work.

We would also like to emphasize that many kernel functions which do not at first appear to be of the form $f(\|u-v\|_{2}^{2})$ can be rearranged appropriately to be of this form. For instance, in Section 10 below we show that the recently popular Neural Tangent Kernel is of this form, so our results apply to it as well. That said, to emphasize that our results are very general, we will mention later how they also apply to some functions of the form $\mathsf{K}(u,v)=f(\langle u,v\rangle)$, including $\mathsf{K}(u,v)=|\langle u,v\rangle|$.

Matrix-vector multiplication, spectral sparsification, and Laplacian system solving are very natural linear algebraic problems in this setting, and have many applications beyond the three we have focused on ($n$-body simulation, spectral clustering, and semi-supervised learning). See Section 1.5 below where we expand on more applications.

Finally, we discuss dependencies on the dimension $d$ and the accuracy $\varepsilon$ for which $n^{1+o(1)}$ algorithms are possible (part (b)). Define $\alpha$, a measure of the ‘diameter’ of the point set and $f$, as

It is helpful to have the following two questions in mind when reading our results:

(Low-dimensional algorithms, e.g. $d=o(\log n)$) Is there an algorithm which runs in time $(\log(n\alpha/\varepsilon))^{O(d)}n^{1+o(1)}$ for multiplication and Laplacian solving? Is there a sparsification algorithm which runs in time $(\log(n\alpha))^{O(d)}n^{1+o(1)}$ when $\varepsilon=1/2$?

We will see that there are many important functions $\mathsf{K}$ for which there are such efficient low-dimensional algorithms, but no such efficient high-dimensional algorithms. In other words, these functions $\mathsf{K}$ suffer from the classic ‘curse of dimensionality.’ At the same time, other functions $\mathsf{K}$ will allow for efficient low-dimensional and high-dimensional algorithms, while others won’t allow for either.

We now state our results. We will give very general classifications of functions $\mathsf{K}$ for which our results hold, but afterwards in Section 1.4 we summarize the results for a few particular functions $\mathsf{K}$ of interest. The main goal of our results is as follows:

Goal: For each problem of interest (part (c)) and dimension $d$ (part (b)), find a natural parameter $p_{f}>0$ associated with the function $f$ for which the following dichotomy holds:

If $p_{f}$ is high, then the problem cannot be solved in subquadratic time assuming $\mathsf{SETH}$ on points in dimension $d$.

If $p_{f}$ is low, then the problem of interest can be solved in almost-linear time ($n^{1+o(1)}$ time) on points in dimension $d$.

As we will see shortly, the two parameters $p_{f}$ which will characterize the difficulties of our problems of interest in most settings are the approximate degree of $f$, and a parameter related to how multiplicatively Lipschitz $f$ is. We define both of these in the next section.

If $f$ can be $\varepsilon$-additively-approximated by a polynomial of degree at most $o(\log n)$, then the problem can be solved in $n^{1+o(1)}$ time.

Otherwise, assuming $\mathsf{SETH}$, the problem requires time $n^{2-o(1)}$.

The same holds for $L_{G}$, the Laplacian matrix of $G$, replaced by $A_{G}$, the adjacency matrix of $G$.

While Theorem 1.1 yields a parameter $p_{f}$ that characterizes hardness of multiplication in high dimensions, it is somewhat cumbersome to use, as it can be challenging to show that a function is far from a polynomial. We also show Theorem 5.15, which shows that if $f$ has a single point with large $\Theta(\log n)$-th derivative, then the problem requires time $n^{2-o(1)}$ assuming $\mathsf{SETH}$. The Strong Exponential Time Hypothesis ($\mathsf{SETH}$) is a common assumption in fine-grained complexity regarding the difficulty of solving the Boolean satisfiability problem; see section 3.6 for more details. Theorem 1.1 informally says that assuming $\mathsf{SETH}$, the curse of dimensionality is inherent in performing adjacency matrix-vector multiplication. In particular, we directly apply this result to the $n$-body problem discussed at the beginning:

Assuming $\mathsf{SETH}$, in dimension $d=\Theta(\log n)$ one step of the $n$-body problem requires time $n^{2-o(1)}$.

The fast multipole method of Greengard and Rokhlin [GR87, GR89] solves one step of this $n$-body problem in time $(\log(n/\varepsilon))^{O(d)}n^{1+o(1)}$. Our Corollary 1.2 shows that assuming $\mathsf{SETH}$, such an exponential dependence on $d$ is required and cannot be improved. To the best of our knowledge, this is the first time such a formal limitation on fast multipole methods has been proved. This hardness result also applies to fast multipole methods for other popular kernels, like the Gaussian kernel $\mathsf{K}(u,v)=\exp(-\|u-v\|_{2}^{2})$, as well.

1.2 Sparsification

This Theorem applies even when $d=\Omega(\log n)$. When $L$ is constant, the running time simplifies to $O(nd\sqrt{\log n}+n^{1+o(1)}\log\alpha/\varepsilon^{2})$. This covers the case when $f(x)$ is any rational function with non-negative coefficients, like $f(z)=z^{L}$ or $f(z)=z^{-L}$.

It may seem more natural to instead define $L$-multiplicatively Lipschitz functions, without the parameter $C$, as functions with $\rho^{-L}f(x)\leq f(\rho x)\leq\rho^{L}f(x)$ for all $\rho$ and $x$. Indeed, an $L$-multiplicatively Lipschitz function is also $(C,L)$-multiplicative Lipschitz for any $C>1$, so our results show that efficient sparsification is possible for such functions. However, the parameter $C$ is necessary to characterize when efficient sparsification is possible. Indeed, as in Theorem 1.3 above, it is sufficient for $f$ to be $(C,L)$-multiplicative Lipschitz for a $C$ that is bounded away from 1. To complement this result, we also show a lower bound for sparsification for any function $f$ which is not $(C,L)$-multiplicatively Lipschitz for any $L$ and sufficiently large $C$:

For example, when $L=\Theta(\log^{2+\delta}n)$ for some constant $\delta>0$, Theorem 1.4 shows that there is a $C$ for which, whenever $f$ is not $(C,L)$-multiplicatively Lipschitz, the sparsification problem cannot be solved in time $n^{1+o(1)}$ assuming $\mathsf{SETH}$.

Bounding $C$ in terms of $L$ above is important. For example, if $C$ is small enough that $C^{L}=2$, then $f$ could be close to constant. Such $\mathsf{K}$-graphs are easy to sparsify by uniformly sampling edges, so one cannot hope to show hardness for such functions.

Theorem 1.4 shows that geometric graphs for threshold functions, the exponential function, and the Gaussian kernel do not have efficient sparsification algorithms. Furthermore, this hardness result essentially completes the story of which decreasing functions can be sparsified in high dimensions, modulo a gap of $L^{.48}$ versus $L$ in the exponent. The tractability landscape is likely much more complicated for non-decreasing functions. That said, many of the kernels used in practice, like the Gaussian kernel, are decreasing functions of distance, so our dichotomy applies to them.

1.3 Laplacian solving

Laplacian system solving has a similar tractability landscape to that of adjacency matrix multiplication. We prove the following algorithmic result for solving Laplacian systems:

In Section 4, we will see, using an iterative refinement approach, that if a $\mathsf{K}$ graph can be efficiently sparsified, then there is an efficient Laplacian multiplier for $\mathsf{K}$ graphs if and only if there is an efficient Laplacian system solver for $\mathsf{K}$ graphs. Corollary 1.7 then follows using this connection: it describes functions which we have shown have efficient sparsifiers but not efficient multipliers.

Corollary 1.7, which is the first of our two hardness results in this setting, applies to slowly-growing functions that do not have low-degree polynomial approximations, like $f(z)=1/(1+z)$. Next, we state our second hardness result:

This yields a quadratic time hardness result when $L=\Omega(\log^{2}n)$. By comparison, the first hardness result, Corollary 1.7, only applied for $L=o(\log n)$. In particular, this shows that for non-Lipschitz functions like the Gaussian kernel, the problem of solving Laplacian systems and, in particular, doing semi-supervised learning, cannot be done in almost-linear time assuming $\mathsf{SETH}$.

2 Our Techniques

We begin by showing a simple, generic equivalence between $\mathsf{KLapE}$ and $\mathsf{KAdjE}$ for any $\mathsf{K}$: an algorithm for either one can be used as a black box to design an algorithm for the other with only negligible blowups to the running time and error. It thus suffices to design algorithms and prove lower bounds for $\mathsf{KAdjE}$.

We use two primary algorithmic tools: the Fast Multipole Method (FMM), and a ‘kernel method’ for approximating $A_{G}$ by a low-rank matrix.

FMM is an algorithmic technique for computing aggregate interactions between $n$ bodies which has applications in many different areas of science. Indeed, when the interactions between bodies is described by our function $\mathsf{K}$, then the problem solved by FMM coincides with our $\mathsf{KAdjE}$ problem.

Most past work on FMM either considers the low-dimensional case, in which $d$ is a small constant, or else the low-error case, in which $\varepsilon$ is a constant. Thus, much of the literature does not consider the simultaneous running time dependence of FMM on $\varepsilon$ and $d$. In order to solve $\mathsf{KAdjE}$, we need to consider the high-dimensional, high-error case. We thus give a clean mathematical overview and detailed analysis of the running time of FMM in Section 9, following the seminal work of Greengard and Strain [GS91], which may be of independent interest.

As discussed in section 1.1 above, the running time of FMM depends exponentially on $d$, and so it is most useful in the low-dimensional setting. Our main algorithmic tool in high dimensions is a low-rank approximation technique: we show that when $f(x)$ can be approximated by a sufficiently low-degree polynomial (e.g. any degree $o(\log n)$ suffices in dimension $\Theta(\log n)$), then we can quickly find a low-rank approximation of the adjacency matrix $A_{G}$, and use this to efficiently multiply by a vector. Although this seems fairly simple, in Theorem 1.1 we show it is optimal: when such a low-rank approximation is not possible in high dimensions, then $\mathsf{SETH}$ implies that $n^{2-o(1)}$ time is required for $\mathsf{KAdjE}$.

The simplest way to show that $f(x)$ can be approximated by a low-degree polynomial is by truncating its Taylor series. In fact, the FMM also requires that a truncation of the Taylor series of $f$ gives a good approximation to $f$. By comparison, the FMM puts more lenient restrictions on what degree the series must be truncated to in low dimensions, but in exchange adds other constraints on $f$, including that $f$ must be monotone. See Section 9.3 and Corollary 5.13 for more details.

We now sketch the proof of Theorem 1.1, our lower bound for $\mathsf{KAdjE}$ for many functions $\mathsf{K}$ in high enough dimensions (typically $d=\Omega(\log n)$), assuming $\mathsf{SETH}$. Although $\mathsf{SETH}$ is a hardness hypothesis about the Boolean satisfiability problem, a number of recent results [AW15, Rub18, Che18, SM19] have showed that it implies hardness for a variety of nearest neighbor search problems. Our lower bound approach is hence to show that $\mathsf{KAdjE}$ is useful for solving nearest neighbor search problems.

Similarly, for any nonnegative reals $a,b\geq 0$, we can take an appropriate affine transformation of $X$ so that an algorithm for $\mathsf{KAdjE}$ can estimate

The main tool we need for this approach is a way to pick $a_{1},\ldots,a_{d},b_{1},\ldots,b_{d}$ for a function $f$ which cannot be approximated by a low degree polynomial so that $M$ has large determinant. We do this by decomposing $\det(M)$ in terms of the derivatives of $f$ using the Cauchy-Binet formula, and then noting that if $f$ cannot be approximated by a polynomial, then many of the contributions in this sum must be large. The specifics of this construction are quite technical; see section 5.4 for the details.

Prior work (e.g. [CS17], [BIS17], [BCIS18]) has shown $\mathsf{SETH}$-based fine-grained complexity results for matrix-related computations. For instance, [BIS17] showed hardness results for exact algorithms for many machine-learning related tasks, like kernel PCA and gradient computation in training neural networks, while [CS17] and [BCIS18] showed hardness results for kernel density estimation. In all of this work, the authors are only able to show hardness for a limited set of kernels. For example, [BIS17] shows hardness for kernel PCA only for Gaussian kernels. These limitations arise from the technique used. To show hardness, [BIS17] exploits the fact that the Gaussian kernel decays rapidly to obtain a gap between the completeness and soundness cases in approximate nearest neighbors, just as we do for functions $f$ like $f(x)=(1/x)^{\Omega(\log n)}$. The hardness results of [CS17] and [BCIS18] employ a similar idea.

As discussed in Lower Bound Techniques, we circumvent these limitations by showing that applying the multiplication algorithm for one kernel a small number of times and linearly combining the results is enough to solve Hamming closest pair. This idea is enough to give a nearly tight characterization of the analytic kernels for which subquadratic-time multiplication is possible in dimension $d=\Theta(\log n)$. As a result, by combining with reductions similar to those from past work, our lower bound also applies to a variety of similar problems, including kernel PCA, for a much broader set of kernels than previously known; see Section 5.7 below for the details.

Our lower bound is also interesting when compared with the Online Matrix-Vector Multiplication (OMV) Conjecture of Henzinger et al. [HKNS15]. In the OMV problem, one is given an $n\times n$ matrix $M$ to preprocess, then afterwards one is given a stream $v_{1},\ldots,v_{n}$ of length-$n$ vectors, and for each $v_{i}$, one must output $M\times v_{i}$ before being given $v_{i+1}$. The OMV Conjecture posits that one cannot solve this problem in total time $n^{3-\Omega(1)}$ for a general matrix $M$. At first glance, our lower bound may seem to have implications for the OMV Conjecture: For some kernels $\mathsf{K}$, our lower bound shows that for an input set of points $P$ and corresponding adjacency matrix $A_{G}$, and input vector $v_{i}$, there is no algorithm running in time $n^{2-\Omega(1)}$ for multiplying $A_{G}\times v_{i}$, so perhaps multiplying by $n$ vectors cannot be done in time $n^{3-\Omega(1)}$. However, this is not necessarily the case, since the OMV problem allows $O(n^{2.99})$ time for preprocessing $A_{G}$, which our lower bound does not incorporate. More broadly, the matrices $A_{G}$ which we study, which have very concise descriptions compared to general matrices, are likely not the best candidates for proving the OMV Conjecture. That said, perhaps our results can lead to a form of the OMV Conjecture for geometric graphs with concise descriptions.

2.2 Sparsification

Our algorithm for constructing high-dimensional sparsifiers for $\mathsf{K}(x,y)=f(\|x-y\|_{2}^{2})$, when $f$ is a $(2,L)$ multiplicatively Lipschitz function, involves using three classic ideas: the Johnson Lindenstrauss lemma of random projection [JL84, IM98], the notion of well-separated pair decomposition from Callahan and Kosaraju [CK93, CK95], and spectral sparsification via oversampling [SS11, KMP10]. Combining these techniques carefully gives us the bounds in Theorem 1.3.

To overcome the ‘curse of dimensionality’, we use the Lindenstrauss lemma to project onto $\sqrt{L\log n}$ dimensions. This preserves all pairs distance, with a distortion of at most $2^{O(\sqrt{\log n/L})}$. Then, using a $1/2$-well-separated pair decomposition partitions the set of projected distances into bicliques, such that each biclique has edges that are no more than $2^{O(\sqrt{\log n/L})}$ larger than the smallest edge in the biclique. This ratio will upper bound the maximum leverage score of an edge in this biclique in the original $\mathsf{K}$-graph. Each biclique in the set of projected distances has a one-to-one correspondence to a biclique in the original $\mathsf{K}$-graph. Thus to sparsify our $\mathsf{K}$-graph, we sparsify each biclique in the $\mathsf{K}$-graph by uniform sampling, and take the union of the resulting sparsified bicliques. Due to the $(2,L)$-Lipschitz nature of our function, it is guaranteed that the longest edge in any biclique (measured using $\mathsf{K}(x,y)$) is at most $2^{O(\sqrt{L\log n})}$. This upper bounds the maximum leverage score of an edge in this biclique with respect to the $\mathsf{K}$-graph, which then can be used to upper bound the number of edges we need to sample from each biclique via uniform sampling. We take the union of these sampled edges over all bicliques, which gives our results for high-dimensional sparsification summarized in Theorem 1.3. Details can be found in the proof of Theorem 6.3 in Section 6. When $L$ is constant, we get almost linear time sparsification algorithms.

For low dimensional sparsification, we skip the Johnson Lindenstrauss step, and use a $(1+1/L)$-well separated pair decomposition. This gives us a nearly linear time algorithm for sparsifying $(C,L)$ multiplicative Lipschitz functions, when $(2L)^{O(d)}$ is small, which covers the case when $d$ is constant and $L=n^{o(1)}$. See Theorem 6.9 for details.

To prove lower bounds on sparsification for decreasing functions that are not $(C_{L},L)$-multiplicatively Lipschitz, we reduce from exact bichromatic nearest neighbors on two sets of points $A$ and $B$. In high dimensions, nearest neighbors is hard even for Hamming distance [Rub18], so we may assume that $A,B\subseteq\{0,1\}^{d}$. In low dimensions, we may assume that the coordinates of points in $A$ and $B$ consist of integers on at most $O(\log n)$ bits. In both cases, the set of possible distances between points in $A$ and $B$ is discrete. We take advantage of the discrete nature of these distance sets to prove a lower bound. In particular, $C_{L}$ is set so that $C_{L}$ is the smallest ratio between any two possible distances between points in $A$ and $B$. To see this in more detail, see Lemma 8.4.

Let $x_{f}$ be a point at which the function $f$ is not $(C_{L},L)$-multiplicatively Lipschitz and suppose that we want to solve the decision problem of determining whether or not $\min_{a\in A,b\in B}\|a-b\|_{2}\leq k$. We can do this using sparsification by scaling the points in $A$ and $B$ by a factor of $k/x_{f}$, sparsifying the $\mathsf{K}$-graph on the resulting points, and thresholding based on the total weight of the resulting $A$-$B$ cut. If there is a pair with distance at most $k$, there is an edge crossing the cut with weight at least $f(x_{f})$ because $f$ is a decreasing function. Therefore, the sparsifier has total weight at least $f(x_{f})/(1+\varepsilon)=f(x_{f})/2$ crossing the $A$-$B$ cut by the cut sparsification approximation guarantee. If there is not a pair with distance at most $k$, no edges crossing the cut with weight larger than $f(C_{L}x_{f})\leq C_{L}^{-L}f(x_{f})\leq(1/n^{10})\cdot f(x_{f})$ by choice of $C_{L}$. Therefore, the total weight of the $A$-$B$ cut is at most $(1/n^{8})\cdot f(x_{f})$, which means that it is at most $((1+\varepsilon)/n^{8})\cdot f(x_{f})<f(x_{f})/4$ in the sparsifier. In particular, thresholding correctly solves the decision problem and one sparsification is enough to solve bichromatic nearest neighbors.

Before proceeding to the body of the paper, we summarize our results. Recall that we consider three linear-algebraic problems in this paper, along with two different dimension settings (low and high). This gives six different settings to consider. We now define $p_{f}$ in each of these settings. In all high-dimensional settings, we have found a definition of $p_{f}$ that characterizes the complexity of the problem. In some low-dimensional settings, we do not know of a suitable definition for $p_{f}$ and leave this as an open problem. For simplicity, we focus here only on decreasing functions $f$, although all of our algorithms, and most of our hardness results, hold for more general functions as well.

High dimensions: $p_{f}$ is the minimum degree of any polynomial that $1/2^{\text{poly}(\log n)}$-additively approximates $f$. $p_{f}>\Omega(\log n)$ implies subquadratic-time hardness (Theorem 1.1 part 2), while $p_{f}<o(\log n)$ implies an almost-linear time algorithm (Theorem 1.1, part 1).

Low dimensions: Not completely understood. The fast multipole method yields an almost-linear time algorithm for some functions, like the Gaussian kernel (Theorem 2.1), but functions exist that are hard in low dimensions (Proposition 5.31).

Sparsification. In both settings, $p_{f}$ is the minimum value for which $f$ is $(C,p_{f})$-multiplicatively Lipschitz, where $C=1+1/p_{f}^{c}$ for some constant $c>0$ independent of $f$.

High dimensions: If $p_{f}>\Omega(\log^{2}n)$ and $f$ is nonincreasing, then no subquadratic time algorithm exists (Theorem 1.4). If $p_{f}<o(\log n)$, then an almost-linear time algorithm for sparsification exists (Theorem 1.3).

Low dimensions: There is some constant $t>1$ such that if $p_{f}>\Omega(n^{t})$ and $f$ is nonincreasing, then no subquadratic time algorithm exists (Theorem 2.5).If $p_{f}<n^{o(1/d)}$, then there is a subquadratic time algorithm (Theorem 2.4).

High dimensions: $p_{f}$ is the maximum of the $p_{f}$ values in the Adjacency matrix-vector multiplication and Sparsification settings, with hardness occurring for decreasing functions $f$ if $p_{f}>\Omega(\log^{2}n)$ (Corollary 1.7 combined with Theorem 1.8) and an algorithm existing when $p_{f}<o(\log n)$ (Theorem 1.6).

Low dimensions: Not completely understood, as in the low-dimensional multiplication setting. As in the sparsification setting, we are able to show that there is a constant $t$ such that if $f$ is nonincreasing and $p_{f}>\Omega(n^{t})$ where $p_{f}$ is defined as in the Sparsifiction setting, then no subquadratic time algorithm exists (Theorem 8.6).

Many of our results are not tight for two reasons: (a) some of the hardness results only apply to decreasing functions, and (b) there are gaps in $p_{f}$ values between the upper and lower bounds. However, neither of these concerns are important in most applications, as (a) weight often decreases as a function of distance and (b) $p_{f}$ values for natural functions are often either very low or very high. For example, $p_{f}>\Omega(\text{polylog}(n))$ for all problems for the Gaussian kernel ($f(x)=e^{-x}$), while $p_{f}=O(1)$ for sparsification and $p_{f}>\Omega(\text{polylog}(n))$ for multiplication for the gravitational potential ($f(x)=1/x$). Resolving the gap may also be difficult, as for intermediate values of $p_{f}$, the true best running time is likely an intermediate running time of $n^{1+c+o(1)}$ for some constant $0<c<1$. Nailing down and proving such a lower bound seems beyond the current techniques in fine-grained complexity.

4 Summary of our Results on Examples

To understand our results better, we illustrate how they apply to some examples. For each of the functions $f_{i}$ given below, make the $\mathsf{K}$-graph, where $\mathsf{K}_{i}(u,v)=f_{i}(\|u-v\|_{2}^{2})$:

$f_{1}(z)=z^{k}$ for a positive integer constant $k$.

$f_{2}(z)=z^{c}$ for a negative constant or a positive non-integer constant $c$.

$f_{4}(z)=1$ if $z\leq\theta$ and $f_{4}(z)=0$ if $z>\theta$ for some parameter $\theta>0$ (the threshold kernel).

5 Other Related Work

Linear Program is a fundamental problem in convex optimization. There is a long list of work focused on designing fast algorithms for linear program [Dan47, Kha80, Kar84, Vai87, Vai89, LS14, LS15, CLS19, LSZ19a, Son19, Bra20, BLSS20, SY20, JSWZ20]. For the dense input matrix, the state-of-the-art algorithm [JSWZ20] takes $n^{\max\{\omega,2+1/18\}}\log(1/\varepsilon)$ time, $\omega$ is the exponent of matrix multiplication [AW21]. The solver can run faster when matrix $A$ has some structures, e.g. Laplacian matrix.

A recent line of work by Charikar et al. [CS17, BCIS18] also studies the algorithmic KDE problem. They show, among other things, that kernel density functions for “smooth” kernels $\mathsf{K}$ can be estimated in time which depends only polynomially on the dimension $d$, but which depends polynomially on the error $\varepsilon$. We are unfortunately unable to use their algorithms in our setting, where we need to solve $\mathsf{KAdjE}$ with $\varepsilon=n^{-\Omega(1)}$, and the algorithms of Charikar et al. do not run in subquadratic time. We instead design and make use of algorithms whose running times have only polylogarithmic dependences on $\varepsilon$, but often have exponential dependences on $d$.

Kernel functions are useful functions in data analysis, with applications in physics, machine learning, and computational biology [Sou10]. There are many kernels studied and applied in the literature; we list here most of the popular examples.

The following kernels are of the form $\mathsf{K}(x,y)=f(\|x-y\|_{2}^{2})$, which we study in this paper: the Gaussian kernel [NJW02, RR08], exponential kernel, Laplace kernel [RR08], rational quadratic kernel, multiquadric kernel [BG97], inverse multiquadric kernel [Mic84, Mar12], circular kernel [BTFB05], spherical kernel, power kernel [FS03], log kernel [BG97, Mar12], Cauchy kernel [RR08], and generalized T-Student kernel [BTF04].

For these next kernels, it is straightforward that their corresponding graphs have low-rank adjacency matrices, and so efficient linear algebra is possible using the Woodbury Identity (see Section 3.4 below): the linear kernel [SSM98, MSS+99, Hof07, Shl14] and the polynomial kernel [CV95, GE08, BHOS+08, CHC+10].

Finally, the following relatively popular kernels are not of the form we directly study in this paper, and we leave extending our results to them as an important open problem: the Hyperbolic tangent (Sigmoid) kernel [HS97, BTB05a, JKL09, KSH12, ZSJ+17, ZSD17, SSB+17, ZSJD19], spline kernel [Gun98, Uns99], B-spline kernel [Ham04, Mas10], Chi-Square kernel [VZ12], and the histogram intersection kernel and generalized histogram intersection [BTB05b]. More interestingly, our result also can be applied to Neural Tangent Kernel [JGH18], which plays a crucial role in the recent work about convergence of neural network training [LL18, DZPS19, AZLS19b, AZLS19a, SY19, BPSW21, LSS+20, JMSZ20]. For more details, we refer the readers to Section 10.

The authors would like to thank Lijie Chen for helpful suggestions in the hardness section and explanation of his papers. The authors would like to thank Sanjeev Arora, Simon Du, and Jason Lee for useful discussions about the neural tangent kernel.

Summary of Low Dimensional Results

In the results we’ve discussed so far, we show that in high-dimensional settings, the curse of dimensionality applies to a wide variety of functions that are relevant in applications, including the Gaussian kernel and inverse polynomial kernels. Luckily, in many settings, the points supplied as input are very low-dimensional. In the classic $n$-body problem, for example, the input points are 3-dimensional. In this subsection, we discuss our results pertaining to whether algorithms with runtimes exponential in $d$ exist; such algorithms can still be efficient in low dimensions $d=o(\log n)$.

The prior work on the fast multipole method [GR87, GR88, GR89] yields algorithms with runtime $(\log(n/\varepsilon))^{O(d)}n^{1+o(1)}$ for $\varepsilon$-approximate adjacency matrix-vector multiplication for a number of functions $\mathsf{K}$, including when $\mathsf{K}(u,v)=\frac{1}{\|u-v\|_{2}^{c}}$ for a constant $c$ and when $\mathsf{K}(u,v)=e^{-\|u-v\|_{2}^{2}}$. In order to explain what functions $\mathsf{K}$ the fast multipole methods work well for, and to clarify dependencies on $d$ in the literature, we give a complete exposition of how the fast multipole method of [GS91] works on the Gaussian kernel:

The fast multipole method is fairly general, and so similar algorithms also exist for a number of other functions $\mathsf{K}$; see Section 9.3 for further discussion. Unlike in the high-dimensional case, we do not have a characterization of the functions for which almost-linear time algorithms exist in near-constant dimension. We leave this as an open problem. Nonetheless, we are able to show lower bounds, even in barely super-constant dimension $d=\exp(\log^{*}(n))$Here, $\log^{*}(n)$ denotes the very slowly growing iterated logarithm of $n$., on adjacency matrix-vector multiplication for kernels that are not multiplicatively Lipschitz:

This includes threshold functions, but does not include piecewise exponential functions. Piecewise exponential functions do have efficient adjacency multiplication algorithms by Theorem 9.13.

To illustrate the complexity of the adjacency matrix-vector multiplication problem in low dimensions, we are also able to show hardness for the function $\mathsf{K}(u,v)=|\langle u,v\rangle|$ in nearly constant dimensions. By comparison, we are able to sparsify for this function $\mathsf{K}$, even in very high $d=n^{o(1)}$ dimensions (in Theorem 1.5 above).

2 Sparsification

We are able to give a characterization of the decreasing functions for which sparsification is possible in near-constant dimension. We show that a polynomial dependence on the multiplicative Lipschitz constant is allowed, unlike in the high-dimensional setting:

Let $f$ be a $(1+1/L,L)$-multiplicatively Lipschitz function and let $\mathsf{K}(u,v)=f(\|u-v\|_{2}^{2})$. Then an $(1\pm\varepsilon)$-spectral sparsifier for the $\mathsf{K}$-graph on $n$ points can be found in $n^{1+o(1)}L^{O(d)}(\log\alpha)/\varepsilon^{2}$ time.

Thus, geometric graphs for piecewise exponential functions with $L=n^{o(1)}$ can be sparsified in almost-linear time when $d$ is constant, unlike in the case when $d=\Omega(\log n)$. In particular, spectral clustering can be done in $O(kn^{1+o(1)})$ time for $k$ clusters in low dimensions. Unfortunately, not all geometric graphs can be sparsified, even in nearly constant dimensions:

There are constants $c^{\prime}\in(0,1),c>1$ and a value $C_{L}$ given $L>1$ for which any decreasing function $f$ that is not $(C_{L},L)$-multiplicatively Lipschitz does not have an $O(nL^{c^{\prime}})$ time sparsification algorithm for $\mathsf{K}$-graphs on $c^{\log^{*}n}$ dimensional points, where $\mathsf{K}(u,v)=f(\|u-v\|_{2}^{2})$.

This theorem shows, in particular, that geometric graphs of threshold functions are not sparsifiable in subquadratic time even for low-dimensional pointsets. These two theorems together nearly classify the decreasing functions for which efficient sparsification is possible, up to the exponent on $L$.

3 Laplacian solving

As in the case of multiplication, we are unable to characterize the functions for which solving Laplacian systems can be done in almost-linear time in low dimensions. That said, we still have results for many functions $\mathsf{K}$, including most kernel functions of interest in applications. We prove most of these using the aforementioned connection from Section 4: if a $\mathsf{K}$ graph can be efficiently sparsified, then there is an efficient Laplacian multiplier for $\mathsf{K}$ graphs if and only if there is an efficient Laplacian system solver for $\mathsf{K}$ graphs.

For the kernels $\mathsf{K}(u,v)=1/\|u-v\|_{2}^{c}$ for constants $c$ and the piecewise exponential kernel, we have almost-linear time algorithms in low dimensions by Theorems 6.3 and 6.9 respectively. Furthermore, the fast multipole method yields almost-linear time algorithms for multiplication. Therefore, there are almost-linear time algorithm for solving Laplacian systems in geometric graphs for these kernels.

A similar approach also yields hardness results. Theorem 1.5 above implies that an almost-linear time algorithm for solving Laplacian systems on $\mathsf{K}$-graphs for $\mathsf{K}(u,v)=|\langle u,v\rangle|$ yields an almost-linear time algorithm for $\mathsf{K}$-adjacency multiplication. However, no such algorithm exists assuming $\mathsf{SETH}$ by Theorem 2.3 above. Therefore, $\mathsf{SETH}$ implies that no almost-linear time algorithm for solving Laplacian systems in this kernel can exist.

Preliminaries

Our results build off of algorithms and hardness results from many different areas of theoretical computer science. We begin by defining the relevant notation and describing the important past work.

For any function $f$, we write $\widetilde{O}(f)$ to denote $f\cdot\log^{O(1)}(f)$. In addition to $O(\cdot)$ notation, for two functions $f,g$, we use the shorthand $f\lesssim g$ (resp. $\gtrsim$) to indicate that $f\leq Cg$ (resp. $\geq$) for an absolute constant $C$. We use $f\eqsim g$ to mean $cf\leq g\leq Cf$ for constants $c,C$.

For a matrix $A$, we use $\|A\|_{2}$ to denote the spectral norm of $A$. Let $A^{\top}$ denote the transpose of $A$. Let $A^{\dagger}$ denote the Moore-Penrose pseudoinverse of $A$. Let $A^{-1}$ denote the inverse of a full rank square matrix.

We use $G_{\text{grav}}$ to denote the Gravitational constant.

We define $\alpha$ slightly differently in different sections. Note that both are less than the value of $\alpha$ used in Theorem 1.3:

2 Graph and Laplacian Notation

A useful notion related to Laplacian matrices is the effective resistance of a pair of nodes:

The effective resistance of a pair of vertices $u,v\in V_{G}$ is defined as

Using effective resistance, we can define leverage score

The leverage score of an edge $e=(u,v)\in E_{G}$ is defined as

We define a useful notation called electrical flow

We let $d(i)$ denote the degree of vertex $i$. For any set $S\subseteq V$, we define volume of $S$: $\mu(S)=\sum_{i\in S}d(i)$. It is obvious that $\mu(V)=2|E|$. For any two sets $S,T\subseteq V$, let $E(S,T)$ be the set of edges connecting a vertex in $S$ with a vertex in $T$. We call $\Phi(S)$ to be the conductance of a set of vertices $S$, and can be formally defined as

We define the notation conductance, which is standard in the literature of graph partitioning and graph clustering [ST04, KVV04, ACL06, AP09, LRTV11, ZLM13, OZ14].

The conductance of a graph $G$ is defined as follows:

A graph $G$ with minimum conductance $\Phi_{G}$ has the property that for every pair of vertices $u,v$,

where $c_{u}$ is the sum of the weights of edges incident with $u$. Furthermore, for every pair of vertices $u,v$,

3 Spectral Sparsification via Random Sampling

Here, we state some well known results on spectral sparsification via random sampling, from previous works. The theorems below are essential for our results on sparsifying geometric graphs quickly.

Consider a graph $G=(V,E)$ with edge weights $w_{e}>0$ and probabilities $p_{e}\in(0,1]$ assigned to each edge and parameters $\delta\in(0,1),\varepsilon\in(0,1)$. Generate a reweighted subgraph $H$ of $G$ with $q$ edges, with each edge $e$ sampled with probability $p_{e}/t$ and added to $H$ with weight $w_{e}t/(p_{e}q)$, where $t=\sum_{e\in E}p_{e}$. If

$q\geq C\cdot\varepsilon^{-2}\cdot t\log t\cdot\log(1/\delta)$, where $C>1$ is a sufficiently large constant

$p_{e}\geq w_{e}\cdot\mathtt{Reff}_{G}(u,v)$ for all edges $e=\{u,v\}$ in $G$

then $(1-\varepsilon)L_{G}\preceq L_{H}\preceq(1+\varepsilon)L_{G}$ with probability at least $1-\delta$.

There is a $\widetilde{O}(m(\log\alpha)/\varepsilon^{2})$ time algorithm which on input $\varepsilon>0$ and $G=(V,E,w)$ with $\alpha=w_{\max}/w_{\min}$ computes a $(24\log n/\varepsilon^{2})\times n$ matrix $\widetilde{Z}$ such that with probability at least $1-1/n$,

The following is an immediate corollary of Theorems 3.6 and 3.7:

There is a $\widetilde{O}(m(\log\alpha)/\varepsilon^{2})$ time algorithm which on input $\varepsilon>0$ and $G=(V,E,w)$ with $\alpha=w_{\max}/w_{\min}$, produces an $(1\pm\varepsilon)$-approximate sparsifier for $G$.

4 Woodbury Identity

where $A,U,C$ and $V$ all denote matrices of the correct (conformable) sizes: For integers $n$ and $k$, $A$ is $n\times n$, $U$ is $n\times k$, $C$ is $k\times k$ and $V$ is $k\times n$.

The Woodbury identity is useful for solving linear systems in a matrix $M$ which can be written as the sum of a diagonal matrix $A$ and a low-rank matrix $UV$ for $k\ll n$ (setting $C=I$).

5 Tail Bounds

We will use several well-known tail bounds from probability theory.

Let $X=\sum_{i=1}^{n}X_{i}$, where $X_{i}=1$ with probability $p_{i}$ and $X_{i}=0$ with probability $1-p_{i}$, and all $X_{i}$ are independent. Let $\mu=\mathbf{E}[X]=\sum_{i=1}^{n}p_{i}$. Then 1. $\Pr[X\geq(1+\delta)\mu]\leq\exp(-\delta^{2}\mu/3)$, $\forall\delta>0$ ; 2. $\Pr[X\leq(1-\delta)\mu]\leq\exp(-\delta^{2}\mu/2)$, $\forall 0<\delta<1$.

Let $X_{1},\cdots,X_{n}$ denote $n$ independent bounded variables in $[a_{i},b_{i}]$. Let $X=\sum_{i=1}^{n}X_{i}$, then we have

6 Fine-Grained Hypotheses

Impagliazzo and Paturi [IP01] introduced the Strong Exponential Time Hypothesis ($\mathsf{SETH}$) to address the complexity of CNF-SAT. Although it was originally stated only for deterministic algorithms, it is now common to extend $\mathsf{SETH}$ to randomized algorithms as well.

For every $\varepsilon>0$ there exists an integer $k\geq 3$ such that CNF-SAT on formulas with clause size at most $k$ (the so called $k$-SAT problem) and $n$ variables cannot be solved in $O(2^{(1-\varepsilon)n})$ time even by a randomized algorithm.

For every $\varepsilon>0$, there is a $c\geq 1$ such that OV cannot be solved in $n^{2-\varepsilon}$ time on instances with $d=c\log n$.

In particular, it is known that $\mathsf{SETH}$ implies OVC [Wil05]. $\mathsf{SETH}$ and OVC are the most common hardness assumption in fine-grained complexity theory, and they are known to imply tight lower bounds for a number of algorithmic problems throughout computer science. See, for instance, the survey [Wil18] for more background.

7 Dimensionality Reduction

We make use of the following binary version of the Johnson-Lindenstrauss lemma due to Achlioptas [Ach03]:

We will also use the following variant of Johnson Lindenstrauss for Euclidean space, for random projections onto $o(\log n)$ dimensions:

For $k=o(\log n)$, with high probability the maximum distortion in pairwise distance obtained from projecting $n$ points into $k$ dimensions (with appropriate scaling) is at most $n^{O(1/k)}$.

8 Nearest Neighbor Search

Our results will make use of a number of prior results, both algorithms and lower bounds, for nearest neighbor search problems.

We provide the definition of the Approximate Nearest Neighbor search problem

By comparison, the best known algorithm for $d=\Omega(\log n)$ for each of these distance measures other than edit distance runs in time about $dn+n^{2-\Omega(\varepsilon^{1/3}/\log(1/\varepsilon))}$ [ACW16].

9 Geometric Laplacian System

Building off of a long line of work on Laplacian system solving [ST04, KMP10, KMP11, KOSZ13, CKM+14], we study the problem of solving geometric Laplacian systems:

where $L_{G}^{\dagger}$ denotes the pseudo-inverse of $L_{G}$ and matrix norm is defined as $\|c\|_{A}=\sqrt{c^{\top}Ac}$.

Equivalence of Matrix-Vector Multiplication and Solving Linear Systems

In this section, we show that for linear systems with sparse preconditioners, approximately solving them is equivalent to approximate matrix multiplication. We begin by formalizing our notion of approximation.

Given a matrix $M$ and a vector $x$, we say that $b$ is an $\varepsilon$-approximate multiplication of $Mx$ if

Given a vector $d$, we say that a vector $y$ is an $\varepsilon$-approximate solution to $My=d$ if $y$ is an $\varepsilon$-approximate multiplication of $M^{{\dagger}}d$.

Before stating the desired reductions, we state a folklore fact about the Laplacian norm:

Lower Bound for $L_{G}$: Let $\lambda_{\min}$ and $\lambda_{\max}$ denote the minimum nonzero and maximum eigenvalues of $D_{G}^{-1/2}L_{G}D_{G}^{-1/2}$ respectively, where $D_{G}$ is the diagonal matrix of vertex degrees. Since $G$ is connected, all cuts have conductance at least $w_{\min}/(n^{2}w_{\max})=1/(n^{2}\alpha)$. Therefore, by Cheeger’s Inequality [Che70], $\lambda_{\min}\geq 1/(2n^{4}\alpha^{2})$. It follows that,

Upper bound for $L_{G}$: $\lambda_{\max}\leq 1$. Therefore,

Proposition 4.2 implies the following equivalent definition of $\varepsilon$-approximate multiplication:

,$b^{\top}{\bf 1}=0$, and $x^{\top}{\bf 1}=0$, then $b$ is an $2n^{3}\alpha^{2}\varepsilon$-approximate multiplication of $L_{G}x$.

Since $b^{\top}{\bf 1}=0$, $(b-L_{G}x)^{\top}{\bf 1}=0$ also. By the upper bound for $L_{G}^{\dagger}$-norms in Proposition 4.2,

Since $x^{\top}{\bf 1}=0$, $x$ has both nonnegative and nonpositive coordinates. Therefore, since $G$ is connected, there exists vertices $a,b$ in $G$ for which $\{a,b\}$ is an edge and for which $|x_{a}-x_{b}|\geq\|x\|_{\infty}/n$. Therefore,

This is the desired result by definition of $\varepsilon$-approximate multiplication. ∎

Since $b^{\top}{\bf 1}=0$, $(b-L_{G}x)^{\top}{\bf 1}=0$ as well. By the lower bound for $L_{G}^{\dagger}$-norms in Proposition 4.2 and the fact that $b$ is an approximate multiplication for $L_{G}x$,

Taking square roots gives the desired result. ∎

Consider an $n$-vertex $w$-weighted graph $G$, let $w_{\min}=\min_{e\in G}w_{e}$, $w_{\max}=\max_{e\in G}w_{e}$, $\alpha=w_{\max}/w_{\min}$, and $H$ be a known graph for which

The algorithm MultiplyG (Algorithm 2) uses standard preconditioned iterative refinement. It is applied in the opposite from the usual way. Instead of using iterative refinement to solve a linear system given matrix-vector multiplication, we use iterative refinement to do matrix-vector multiplication given access to a linear system solver.

Reduction in residual and iteration bound: First, we show that

Since $(I-L_{H}L_{G}^{{\dagger}})L_{G}(I-L_{G}^{{\dagger}}L_{H})\preceq 3(1/900)L_{G}$,

Let $x_{\text{final}}$ be the lowest element of the call stack and let $k$ be the number of recursive calls to MultiplyGAdditive (Algorithm 2). By Proposition 4.2,

Error: We start by bounding error in the $L_{G}^{{\dagger}}$ norm. Let $b$ be the output of $\textsc{MultiplyGAdditive}(x,\tau)$ (Algorithm 2). We bound the desired error recursively:

Because $0=\textsc{MultiplyGAdditive}(x_{\text{final}},\tau)$ (Algorithm 2),

By Proposition 4.2 applied to $L_{G}^{{\dagger}}$, $\|b-L_{G}x\|_{L_{G}^{{\dagger}}}\geq\|b-L_{G}x\|_{\infty}/(n\sqrt{w_{\max}})$. Therefore,

Runtime: There is one call to SolveG and one multiplication by $L_{H}$ per call to MultiplyGAdditive (Algorithm 2). Each multiplication by $L_{H}$ takes $O(Z)$ time. As we have shown, there are only $k\leq O(\log(n\alpha/\varepsilon))$ calls to MultiplyGAdditive (Algorithm 2). Therefore, we are done. ∎

2 Matrix-Vector Multiplication Implies Solving Linear Systems

The converse is well-known to be true [ST04, KMP11]:

Consider an $n$-vertex $w$-weighted graph $G$, let $w_{\min}=\min_{e\in G}w_{e}$, $w_{\max}=\max_{e\in G}w_{e}$, $\alpha=w_{\max}/w_{\min}$, and $H$ be a known graph with at most $Z$ edges for which

3 Lower bound for high-dimensional linear system solving

We have shown in this section that if a $\mathsf{K}$ graph can be efficiently sparsified, then there is an efficient Laplacian multiplier for $\mathsf{K}$ graphs if and only if ther is an efficient Laplacian system solver for $\mathsf{K}$ graphs. Here we give one example of how this connection can be used to prove lower bounds for Laplacian system solving:

Matrix-Vector Multiplication

$A_{\mathsf{K},P}$ and $L_{\mathsf{K},P}$ are the adjacency matrix and Laplacian matrix, respectively, of the complete weighted graph on $n$ nodes where the weight between node $i$ and node $j$ is $\mathsf{K}(x_{i},x_{j})$.

In this section, we study the algorithmic problem of computing the linear transformations defined by these matrices:

We make a few important notes about these problems:

In both of the above, problems, $w_{\max}:=\max_{u,v\in P}|\mathsf{K}(u,v)|$.

Suppose the function $\mathsf{K}$ can be evaluated in time $T$ (in this paper we’ve been assuming $T=\widetilde{O}(1)$). Then, both the $\mathsf{KAdjE}$ and $\mathsf{KLapE}$ problems can be solved in $O(Tn^{2})$ time, by computing all $n^{2}$ entries of the matrix and then doing a straightforward matrix-vector multiplication. However, since the input size to the problem is only $O(nd)$ real numbers, we can hope for much faster algorithms when $d=o(n)$. In particular, we will aim for $n^{1+o(1)}$ time algorithms when $d=n^{o(1)}$.

For some functions $\mathsf{K}$, like $\mathsf{K}(x,y)=\|x-y\|_{2}^{2}$, we will show that a running time of $n^{1+o(1)}$ is possible for all $d=n^{o(1)}$. For others, like $\mathsf{K}(x,y)=1/\|x-y\|_{2}^{2}$ and $\mathsf{K}(x,y)=\exp(-\|x-y\|_{2}^{2})$, we will show that such an algorithm is only possible when $d\ll\log(n)$. More precisely, for these $\mathsf{K}$:

When $d=O(\log(n)/\log\log(n))$, we give an algorithm running in time $n^{1+o(1)}$, and

For $d=\Omega(\log n)$, we prove a conditional lower bound showing that $n^{2-o(1)}$ time is necessary.

Finally, for some functions like $\mathsf{K}(x,y)=|\langle x,y\rangle|$, we will show a conditional lower bound showing that $\Omega(n^{2-\delta})$ time is required even when $d=2^{\Omega(\log^{*}n)}$ is just barely super-constant.

In fact, assuming $\mathsf{SETH}$, we will characterize the functions $f$ for which the $\mathsf{KAdjE}$ and $\mathsf{KLapE}$ problems can be efficiently solved in high dimensions $d=\Omega(\log n)$ in terms of the approximate degree of $f$ (see subsection 5.2 below). The answer is more complicated in low dimensions $d=o(\log n)$, and for some functions $f$ we make use of the Fast Multipole Method to design efficient algorithms (in fact, we will see that the Fast Multipole Method solves a problem equivalent to our $\mathsf{KAdjE}$ problem).

Although our goal in this section is to study the $\mathsf{K}$ Laplacian Evaluation problem, it will make the details easier to instead look at the $\mathsf{K}$ Adjacency Evaluation problem. Here we show that any running time achievable for one of the two problems can also be achieved for the other (up to a $\log n$ factor), and so it will be sufficient in the rest of this section to only give algorithms and lower bounds for the $\mathsf{K}$ Adjacency Evaluation problem.

Suppose the $\mathsf{KAdjE}$ (Problem 5.1) can be solved in $\mathcal{T}(n,d,\varepsilon)$ time. Then, the $\mathsf{KLapE}$ (Problem 5.2) can be solved in $O(\mathcal{T}(n,d,\varepsilon/2))$ time.

Suppose the $\mathsf{KLapE}$ (Problem 5.2) can be solved in $\mathcal{T}(n,d,\varepsilon)$ time, and that $\mathcal{T}$ satisfies $\mathcal{T}(n_{1}+n_{2},d,\varepsilon)\geq\mathcal{T}(n_{1},d,\varepsilon)+\mathcal{T}(n_{2},d,\varepsilon)$ for all $n_{1},n_{2},d,\varepsilon$. Then, the $\mathsf{KAdjE}$ (Problem 5.1) can be solved in $O(\mathcal{T}(n\log n,d,0.5\varepsilon/\log n))$ time.

We will show that the $\mathsf{K}$ Adjacency Evaluation problem can be solved in

time, and then apply the superadditive identity for $T$ to get the final running time. For a fixed $d$, we proceed by strong induction on $n$, and assume the $\mathsf{K}$ Adjacency Evaluation problem can be solved in this running time for all smaller values of $n$.

in $O(\mathcal{T}(n,d,0.5\varepsilon/\log n))$ time.

and our two initial calls took time $O(\mathcal{T}(n,d,0.5\varepsilon/\log n))$, leading to the desired running time. Each output entry is ultimately the sum of at most $2\log n$ terms from calls to the given algorithm, and hence has error $\varepsilon$ (since we perform all recursive calls with error $0.5\varepsilon/\log n$ and the additive error guarantees in recursive calls can only be more stringent). ∎

In both Proposition 5.3 and Proposition 5.4, if the input to the $\mathsf{KLapE}$ (resp. $\mathsf{KAdjE}$) problem is a $\{0,1\}$ vector, then we only apply the given $\mathsf{KAdjE}$ ($\mathsf{KLapE}$) algorithm on $\{0,1\}$ vectors. Hence, the two problems are equivalent even in the special case where the input vector $y$ must be a $\{0,1\}$ vector.

2 Approximate Degree

The Stone-Weierstrass theorem says that, for any $\varepsilon>0$, and any continuous function $f$ which is bounded on $$, there is a positive integer$k$such that$f$is$\varepsilon$-close to a polynomial of degree$k$. That said,$k$can be quite large for some natural and important continuous functions$f$. For some examples: