Hardness results for the Discrete Gauss Transform
Kyriakos Axiotis, Aleksandar Makelov, Chenyang Yuan

1

Introduction

The Discrete Gauss Transform (DGT) is an important computational primitive that arises
in various areas of mathematics and its applications. Given a set of n “source” points A ∈ Rd
that defines a Gaussian kernel and a set of n “query” points B ∈ Rd , it asks to evaluate
the Gaussian kernel at every query point. As an example, in the study of heat equations
and fluid dynamics DGT arises as a discretization of certain families of differential equation
solutions. It is thus widely used in physical simulations and computer modeling. Furthermore,evaluating and optimizing over Gaussian kernels is important for Machine Learning.
The trivial algorithm for this problem runs in time O(n2 d) and comprises of explicitly
evaluating the Gaussian kernel at each query point. Since there are Θ(n2 ) pairs of points intotal and it takes O(d) time to compute the distance between two points, the bound follows.
This algorithm, however, is somewhat impractical because of the quadratic dependence on
n. In practice, n will correspond to the number of particles or bodies in a physical system,
and it can be expected to be large. On the other hand, the dependence on d might not be
extremely crucial, at least for physics applications, in which the space can be expected to be
intrinsically low-dimensional. (Note, by the way, that this is not the case in most Machine
Learning applications, where one expects the dimension to be very large).
In light of the above, it is natural to ask whether there exists an algorithm with much better (ideally linear) dependence on n. This question was answered affirmatively by [GS91], who
proposed the Fast Gauss
Transform,
 an algorithm that computes an additive ε-approximation
1 O(d). While the dependence on n is optimal, the dependence on
to DGT in time O n log εd renders the algorithm practical only for applications with very low dimension. As such, a lotof research has been done on devising various heuristics that seem to work somewhat better inpractice, but do not enjoy theoretically better guarantees (see e.g. [GW13, RYDG05, TW09]).
As we will see in this exposition, the apparent barrier in the improvement of the runtime
of Gauss transform algorithms is inherent, as fast algorithms for this problem would imply
faster algorithms for SAT, refuting the well-known Strong Exponential Time Hypothesis.

2Preliminaries

We assume familiarity with the basics of computational complexity theory. In this section,
we describe the computational problems relevant to our reductions. Key to our approach1

is the following classical problem, known as the Orthogonal Vectors Problem (or OV for
short):
Problem 2.1 (Orthogonal Vectors problem (monochromatic version)). Given n vectors
v1 , . . . , vn ∈ {0, 1}d , are there i, j such that hvi , vj i = 0?
More specifically, we will use the bichromatic version of OV, in which we are only concerned with the inner products between vectors from two different sets:
Problem 2.2 (Orthogonal Vectors problem (bichromatic version)). Given two sets of n
vectors A = {v1 , . . . , vn } ∈ {0, 1}d and B = {u1 , . . . , un } ∈ {0, 1}d , are there i, j such that
hvi , uj i = 0?
As we will show, the complexity of mono- and bi-chromatic OV is the same, so we will callboth problems just OV. We will also refer to the orthogonal vectors problem for n vectors (orn vectors in each set for the bichromatic version) in dimension d as OV(n, d) to emphasize
the dependence on parameters. Central to our method is the following standard conjecturefrom fine-grained complexity:
Conjecture 2.3 (Orthogonal vector conjecture). For every ε > 0, there is a c ≥ 1 such that
OV cannot be solved in n2−ε time on instances with d = c log n.
In particular, this conjecture is refuted if there exists an algorithm that solves OV in
c log n dimensions, for large enough constant c, in n2−Ω(1) time. As we will be concerned withvarious geometric problems where points can have arbitrary coordinates, a stepping stone
in our reductions is the following problem, which is a version of OV over Z (here, we only
need the bichromatic version):
Problem 2.4 (Hopcroft’s problem (bichromatic version)). Given two sets of n vectors A ={v1 , . . . , vn } ∈ Zd and B = {u1 , . . . , un } ∈ Zd , are there i, j such that hvi , uj i = 0? All input
numbers consist of D bits.
As with OV, we will also call this problem Hop(n, d, D) to emphasize the dependence onparameters. The next problem we need is the closest pair problem; here, we again only need
the bichromatic variant:
Problem 2.5 (`2 Bichromatic Closest pair). Given two sets of n vectors A = {v1 , . . . , vn } ∈Rd and B = {u1 , . . . , un } ∈ Rd , where each coordinate of each vector has ≤ D bits, find thepair (i, j) such that the `2 distance between vi and uj is the smallest.
We will call this problem BCP(n, d, D). Finally, we formally define the central problem
that we want to prove hardness of: the approximate Gauss transform. Here, we define it overZ, as this will be sufficient for proving hardness:2

Problem 2.6 (Approximate Gauss transform). Given two sets of n vectors A = {v1 , . . . , vn } ∈Zd and B = {u1 , . . . , un } ∈ Zd , and ε > 0, definen

1 X −ku−vi k22
F (u) =en i=1
The approximate Gauss transform problem is to compute Fb(u1 ), . . . , Fb(un ) such that for alli = 1, . . . , n we haveFb(ui ) − F (ui ) ≤ ε
We will call this problem GT(n, d, ε).

2.1Reductions between monochromatic and bichromatic versionsof OV

Here are the promised reductions showing that the mono- and bi-chromatic versions of OV
are equivalent for our purposes.
Lemma 2.7. An instance of monochromatic OV(n, d) reduces to an instance of bichromatic
OV(n, d + 1) in O(nd) (i.e. linear) time.
Proof. First, search for the zero vector among our input vectors and if it is present, reduce
to a trivially true instance of bichromatic OV.
Otherwise, construct the vectors A = {(v1 , 0), . . . , (vn , 0)} and B = {(v1 , 1), . . . , (vn , 1)}
in dimension d + 1 - this is our bichromatic OV instance.
Lemma 2.8. An instance of bichromatic OV(n, d) reduces to an instance of monochromatic
OV(2n, d + 2) in O(nd) (i.e. linear) time.
Proof. If our vectors are A = {u1 , . . . , un } and B = {v1 , . . . , vn }, our monochromatic OV
instance consists of the vectors (ui , 1, 0) and (vi , 0, 1) for 1 ≤ i ≤ n. Note that no two vectors
coming from A and no two vectors coming from B can now have dot product 0, so we’re ineffect searching only over the pairs of vectors in different sets.

3

Previous work

In this section, we summarize the prior work relevant to our main result about the hardnessof the approximate Gauss transform.

3.1

Williams’s reduction from OV to Hopcroft’s problem

We will use the following reduction of Williams from OV to Hop:
Lemma 3.1 (Lemma 1.1 from [Wil18]). Let ` ∈ [1, d]. There 
exists a ndO(d/`) -time reductiondfrom OV(n, d) to dO(d/`) instances of Hop n, ` + 1, O d log
.
`
3

3.2

Chen’s improved reduction from OV to Hopcroft’s problem

Lijie Chen [Che18] proves the following reduction:
Lemma 3.2 (Lemma 1.17 from [Che18]). Let ` ∈ [1, d]. There is alog∗ d d/`) poly(d) -time
O n`O(6log∗ d d/`

reduction from OV(n, d) to `O(63.3

) instances of Hop n, ` + 1, O∗
dlog `6log d
`



.

Williams’ reduction from Hopcroft’s problem to bichromatic
closest pair

Williams [Wil18] also gives the following reduction from Hop to BCP:
Lemma 3.3. There is aO nd2 (D + O(log d))2 -time
reduction from Hop(n, d, D) to BCP (n, d2 , O(D + O(log d))).
Proof. Let A = {u1 , . . . , un } and B = {v1 , . . . , vn } be the two sets of vectors. We construct
new sets
A0 = u0i = vec(ui uTi )/ kui k22 1 ≤ i ≤ nand
B 0 = vi0 = − vec(vi viT )/ kvi k221≤i≤nwhere we regard the matrices as unfolded vectors. Now note that2
vec(ui uTi ) 2 = tr(ui uTi ui uTi ) = tr uTi ui uTi ui = kui k42
and most importantlyvec(ui uTi ), − vec(vi viT ) = − tr ui uTi vi viT = − hui , vi i2
which means that2

ku0i − vi0 k2 = 2 +hui , vi i2.
kui k22 kvi k22

Thus, finding the closest pair in `2 between A0 and B 0 and checking if the distance is 2 ormore tells us the solution to Hopcroft’s problem. Note that the resulting dimension is d2 andwe have n vectors in each set.
Now it remains to deal with the bit complexity of the vectors we produce. Note that if
hui , vi i2 6= 0, then hui , vi i2 ≥ 1, which means that for such pairs we have2

ku0i − vi0 k2 ≥ 2 +41
d2 22D+1

So now suppose we truncate u0i and vi0 to D0 bits after the binary point; then ku0i − vi0 k22
d2changes by at most about 22D
0 , so we want pretty much
022D > d4 22D
or D0 ≥ D + 4 log d guarantees that the closest pair instance will be able to detect a dot
product of 0.
Finally, we scale up all our points so that they have integer coordinates.
Running time: Observe that to multiply/divide two O(B)-bit numbers naively it takesO(B 2 ) time; therefore,• computing all the ui uTi and vi viT takes O(nd2 D2 );• computing all the kui k22 , kvi k22 takes O(ndD2 );• performing the normalization takes O(nd2 (D + O(log d))2 );• scaling up the vectors to Z takes O(nd2 D)
So overall the running time is O(nd2 (D + O(log d))2 ).

4

The hardness of computing the Discrete Gauss TransformIn this section we give our main results for the inapproximability of the Gauss transform
assuming the orthogonal vectors conjecture.

4.1

From bichromatic closest pair to approximate Gauss transform

We start with a reduction from the bichromatic closest pair problem to the approximate
Gauss transform that will be used in all our results:Lemma 4.1. There is a O ndD loglog n + d2D logn -time reduction from an instance of1
BCP(n, d, D) to an instance of GT n, d, O(nO(d2.
D))
Proof. If A and B are the two sets of vectors for the closest pair instance, we run our approximate Gauss transform algorithm on λ-scaled versions A0 , B 0 of these sets with approximation
ε for λ, ε to be chosen later.
Let d1 be the closest distance between A and B, and d2 the second-closest. Note thatd22 ≥ d21 + 1 as we are working over Z. Now, for any u∗ ∈ A such that u∗ realizes the shortest
distance with some vector in B, we have
F (u∗ ) ≥5

1 −λ2 d21enand for any u that doesn’t, we have
2 2

F (u) ≤ e−λ d2 .
Our strategy now is to choose λ so that F (u∗ ) − ε > F (u) + ε for any u ∈ A which
doesn’t realize the shortest distance to B. Note that this will ensure that when we run theapproximate Gauss transform and find the smallest Fb(u∗ ) among all u∗ ∈ A, it will be fora u∗ which realizes the shortest distance to B. Then, for such a u∗ , we just enumerate overv ∈ B and find the shortest distance; this will be the solution to the bichromatic closest pairinstance.
So, we want to have
2 2

2 2

e−λ d1 − ne−λ d2 ≥ 2nε√Choosing λ = dc log ne for c ≥ 2, we have for n ≥ 2 that1
−λ2 d21−λ2 d22−λ2 d21e− ne≥e
1 − c−1n1−c2 d21 log n=e
1 − c−1n
1 1.
≥
2 nc2 d21
1
Since d21 = O(d2D ) and d2D ≥ 1, it suffices to choose ε = O( nO(d2D ) ). This finishes the proof.
√
Running time: We just need to scale the sets A and B by the integer dc log ne which
has O(log log n) bits, so the scaling takes O(ndD log log n) time, and we need to write downε, which takes O(d2D log n) time.

4.2

Reduction from OV to Gauss Transform in log n Dimensions

For instances of the Gauss transform where d = c log n, we can give a simple reduction
from the bichromatic version of the orthogonal vector conjecture to the approximate Gausstransform problem. Note that in this regime the simple exact algorithm that checks all Θ(n2 )pairs beats the FGT algorithm.
The idea is to reduce from OV to BCP directly, without passing through Hop, whichallows us to avoid the more complex reductions by Williams, and also results in a better
dependence on the dimension. Namely, we can avoid the squaring of the dimension from
Williams’ reduction from Hop to BCP.
Lemma 4.2. There is an Õ(n)-time reduction from an instance of OV(n, c log n) to an
instance of GT(n, 3d, ε) wherelog(1/ε) = O(log2 n).
6

Proof. Given the two sets of n vectors A = {u1 , . . . , un } ⊆ {0, 1}d and B = {v1 , . . . , vn } ⊆{0, 1}d in the OV instance, we construct new sets A0 , B 0 ⊆ {0, 1}3d such that for each u ∈ A
and v ∈ B we include u0 in A0 and v 0 in B 0 , whereu0 = [ u1v 0 = [ −v1

···
···

ud
−vd

1 − u10······

1 − ud00
1 − v1······0
1 − vd]]Note that hu0 , u0 i = hv 0 , v 0 i = d and hu0 , v 0 i = − hu, vi by construction. Therefore ku0 − v 0 k2 =2d+hu, vi. Moreover, hu, vi ≥ 0 for all u, v since they are in {0, 1}d . This gives us a reduction
from OV(n, d) to BCP(n, 3d, O(1)), since a dot product of zero between some u ∈ A andv ∈ B is equivalent to to the shortest distance between u0 ∈ A0 and v 0 ∈ B 0 being 2d. This
reduction runs in O(nd) = Õ(n)-time.
Now we can apply this reduction with d = c log n, and compose it with our reduction
1
from BCP(n, 3d, O(1)) to GT(n, 3d, nO(logn) ), which will work in timeeO(nd log log n + d log n) = O(n log n log log n + log2 n) = O(n)as we wanted to show, and moreover we have log(1/ε) = O(d log n) = O(log2 n).
In particular, for this lemma we get our first hardness result for the Gauss transform as
a direct corollary:
Corollary 4.3. Assuming the OV conjecture, for any η, δ > 0 and d ≥ c log n, there is noalgorithm for GT(n, d, ε) that runs in time
1−ηO n2−δ log(1/ε)d
Proof. Assuming such an algorithm existed, applying the above reduction we would get an0O(n2−δ )-time algorithm for OV(n, d0 ) for some 0 < δ 0 < δ with d0 ≥ c log n/3, refuting the
OV conjecture.

4.3

Using Williams’ reduction to get a better lower boundIn this subsection, we show how we can use Williams’ reduction to show a hardness result
for GT which has worse dependence on the degree of log(1/ε), but applies to the moreinteresting range of dimensions (ω((log log n)2 ), o(log n))
Lemma 4.4. Let ` ∈ [1, d]. Then there is a
 2 2 d log dO(d/`)
2nd+ n(` + 1) O
dO(d/`) + dO(d/`) · O n(` + 1)2 log log n + (` + 1)2 2O(d log d/`) log n
2`time reduction from OV(n, d) to dO(d/`) instances of GT(n, (` + 1)2 , ε) where−(`+1)2 2O(d log d/`)ε=O n7

Proof. This follows by combining the reduction from OV to Hop in lemma 3.1 with the
reduction from Hop to BCP in lemma 3.3 and the reduction from BCP to GT in lemma
4.1.
In particular, we have the following corollary:
Corollary 4.5. Assuming the OV conjecture, there is no algorithm for GT(n, d0 , ε) ind0 1/2−η
O(n2−δ log 1ε) time for any δ, η > 0 and d0 = (log log n)ω(1) ∩ O(log n).
Proof. Assume the contrary. In the reduction√from OV to GT, choose d = c log n and(` + 1)2 = d0 so that ` = Ω((log log n)ω(1) ) ∩ O( log n); then note that

`1−η
1
1−η
−η
−η
log
= (` + 1)O(` ) 2O(d` log d) = no(1) nO(` log log n) = no(1)εby the choice of the lower bound of `.
Running time: For this choice of parameters, the total running time of our reduction is
1+o(1)n. Indeed, one sees that dO(d/`) = (c log n)c log n/` = (c log n)o(c log n/ log log n) which is no(1) ,and it’s easy to see that all the other terms are no(1) too.
d0 1/2−η
) approximate Gauss transform algorithm toApplying our supposed O(n2−δ log 1ε2−δ o(1)this gives us a n nalgorithm for the OV instance, which contradicts the OV conjecture.

4.4

Improved bichromatic closest pair hardness and Gauss transform implications

In this section, we show that combining the improved reduction by Chen [Che18] from OV to
Hop with our reductions gives a lower bound for GT for an even better range of dimensions.
Firstly, by combining Chen’s reduction with our above reductions we obtain the following:
Lemma 4.6. There is alog∗ d d/`) poly(d)O n`O(6+ `O(∗
6log d d/`log∗ d d/`+ `O(6) O n(` + 1)2dlog ` · 6log ∗d
`

2 !



∗
) O n(` + 1)2 log log n + `O(6log d d/`)log∗ d d/`time reduction from OV(n, d) to `O(6) instances of GT(n, (` + 1)2 , ε) wherelog∗ d +log `d2 O ` log `·6

ε = n−(`+1) 2



(

8

)

log∗ dd
2 O ` log `·6= n−(`+1) 2

(

)

.

Proof. Follows by composing Chen’s reduction from OV to Hop in lemma 3.2 with the
reduction from Hop to BCP in lemma 3.3 and the reduction from BCP to GT in lemma
4.1.
Using this reduction, we get a lower bound for GT in an improved range of dimensions:
Corollary 4.7. Assuming the OV conjecture, there is no algorithm for GT(n, d0 , ε) ind0 1/2−η
∗O(n2−δ log 1ε) time for any δ, η > 0 and d0 = 6ω(log n) ∩ o(log n).
Proof. We apply Lemma 4.6 with d0 = (` + 1)2 . We want to geta lower bound for themrunning time of an algorithm for GT(n, d0 , ε) of the form n2−δ log 1ε . The idea is to apply
our reduction with d = c log n and an m = `1−δ (where δ > 0) chosen so that
 m
1
log
= no(1)εFor our specified range of d0 and therefore `. From our reduction with d = c log n, we have
log

log n
1
log∗ log n) log n= (` + 1)2 2O( ` log `·6εand we want this to the m-th power to be no(1) , for which it suffices that each term in the
product raised to m be no(1) . We see that to have(log n)m = no(1)
we need m = o logloglogn n , which is true if d0 = (` + 1)2 = o(log n). For the (` + 1)2 term it iseasy to see that for our allowed values of `:m log (` + 1)
`1−δ log (` + 1)=
= o(1)
log n
log n
But the middle term is most troublesome. Since we choose m = `1−δ ,
∗

m

∗

log n
log log n) = n`−δ log `·6log log n2( ` log `·6∗∗

Which is no(1) if ` ≥ 6ω(log n) which also implies d0 ≥ 6ω(log n) . This means that, for example,
the lower bound works for dimensions d0 ≥ log · · · log n, for any finite number of logs.
log∗ d d/`∗) = n6log∗ n log `/` = no(1) , and it
Running time: For our choices of ` ≥ 6ω(log n) , `O(6is easy to see that the other terms are no(1) too, showing that the reduction takes O(n2−ξ )time in total for some ξ > 0

9

4.5

Barriers to Improving the Lower Bound

As we can see from the previous sections, the quality of dimension reduction from OV to
Hop determine the lower range of allowed d0 in our lower bounds for the Gauss transform.
∗
Although 6ω(log n) is already very small, it is conjectured in [Che18] that this can be improvedto ω(1), which would imply that our reduction will work for d ≥ ω(1) dimensions.
In our lower bounds for small dimensions, the exponent on log(1/) is d0 1/2−η , and thiscomes from the reduction from Hop(n, d, D) to BCP (n, d2 , O(D)) in lemma 3.3. If we havea reduction to BCP (n, O(d), O(D)) instead, then we can get a tight exponent of d0 1−η . One
way to get such a reduction is to solve the following problem:
Problem 4.8. Given two sets of n vectors each, A, B ⊂ Zd , can we find a pair of mappings
f , g : Zd → ZO(d) , so that the norms of f (u) are the same for all u ∈ A, the norms of
all g(v) are the same for all v ∈ B, and moreover maxu∈A,v∈B hf (u), g(v)i is achieved iff
hu, vi = 0? The mappings f , g should also not increase the bit-complexity of u, v by morethan a polynomial factor.
2

The reduction in lemma 3.3 can be seen as using the maps f , g : Zd → Zd , wheref (x) = xxT /kxk2 and g(x) = −xxT /kxk2 . While we have spent a long time unsuccessfully
trying to improve this reduction, it is unclear to us whether this is possible, or if there is afundamental barrier preventing us from improving it.
In particular, if we look for a reduction with the property that hf (u), g(v)i ≤ 0 withequality iff hu, vi = 0 (which is the case for the reduction in lemma 3.3), one can show thatin such a reduction the functions f and g cannot depend on their inputs coordinate-wise, i.e.
they cannot simply map each coordinate of u to a number of coordinates of f (u) withoutlooking at the other coordinates of u. This gives some evidence that we may need to lookat all d2 combinations of coordinates (and in fact, f (u) and g(v) from lemma 3.3 do havea coordinate for each pair of coordinates of u and v). However, it does not rule out theexistence of a smarter reduction.

5

Algorithms for the Gauss TransformIn this section we discuss different algorithms for computing the Gauss Transform. We arepresenting three algorithms, each one corresponding to one parameter regime, depending on
the desired runtime dependence on the approximation parameter ε and the dimension d. In
Section 5.1, we describe the naive, exact algorithm, that however is quadratic in the numberof points. In Section 5.2, we present the Fast Gauss Transform, which has a linear dependence
on n, and works well for very low dimensions and good approximation guarantees. Finally,by using Random Fourier Features, we show a simple way to get an algorithm that is linear
in the number of points and the dimension, albeit with a polynomial dependence on theerror.

10

5.1

Naive algorithm

Naively, one can compute the Gauss Transform by iterating over all query points, and evaluating the Gaussian kernel in linear time. This can be done by computing the distance
between any two points. The total runtime of this algorithm is O(n2 d), where d comes fromthe time needed to compute the distance between two points.
Although trivial, any exact algorithm that has subquadratic dependence on n and polynomial dependence in d would refute SETH, as we have seen in the previous hardness results.

5.2

Fast Gauss TransformIn this section we will outline the Fast Gauss Transform [GS91]. The intuition behind thealgorithm is that the Gaussian kernel decays exponentially in the distance squared, and so
if one discretizes the space by subdividing into boxes, the value of a query point will only
depend on points that are in boxes close to the query point. Furthermore, one does not needto evaluate the kernel between the query point and every point in a box: by computing
a Taylor approximation around the center of the box, one only needs to evaluate just one
polynomial for each box B.
More specifically, for any t, s, sB ∈ Rd where sB is the center of box B, the Taylor expan2
sion of e−||t−s||2 for some s ∈ B around the center sB is
X (t − sB )α1 +···+αd2hα (s − sB )e−||t−s||2 =
α
!
.
.
.
α
!
1
d
α ,...,α ≥0
1

d

where

∂ αd −t21 −···−t2d
∂ α1
.
.
.
e
∂tα1 1∂tαd dWe now sum these Taylor series over all s ∈ B and switch the order of summation to
obtain a single Hermite expansion
X
G(t) =Aα hα (t − sB )hα (t) = (−1)α1 +···+αdα1 ,...,αd ≥0

where
Aα =

X1
1
(s1 − sB,1 )α1 . . . (sd − sB,d )αdα1 ! . . . αn ! s∈B n

Now we truncate this series expansion to include only those monomials where all αi are ≤ p.
By using Cramer’s inequality, i.e. the fact that1122|hα (t)| ≤ Ke−(t1 +···+td )/2 2(a1 +···+ad )/2 √α!
a1 ! . . . ad !
where K < (1.09)d , the error due to this truncation is at mostp−Ω(pd) rO(pd)11where r is the side length of the box B. Setting r = 12 and p = O(log 1ε ) we get that the error
is at most ε.
Given any two points p, q with ||p − q||∞ ≥ Ω(log 1ε ), note that22e−||p−q||2 ≤ e−||p−q||∞ ≤ εΩ(1)and so we only need to care about boxes within an `∞ distance of O log 1ε . The number ofO(d).
such boxes is log 1εFrom the above, we deduce that both the time needed to compute the polynomials for allO(d)boxes and the time to evaluate the Gaussian kernel at all query points is n log 1ε. Thisgives the total runtime of the algorithm.

5.3

Random Fourier Features

Another approach for approximately computing the Fast Gauss Transform is by using Random Fourier Features, as introduced by [RR08]. Given a set S ⊆ [−M , M ]d and an error
parameter ε > 0, there exists a mapping f : Rd → RD with D = O( εd2 log Mε ), such that22∀x, y ∈ S e−||x−y||2 − ε ≤ hf (x), f (y)i ≤ e−||x−y||2 + ε
Furthermore, such a mapping can be computed with high probability and f can be appliedin O(D) time. Now this implies that also for all x
1 X −||x−yi ||22
1X
1 X −||x−yi ||22e
−ε≤hx, yi i ≤
e
+ε
n
n
n
i∈[n]

i∈[n]

i∈[n]

however, by linearity, the middle quantity is also equal to
hx,

1X
yi i
n
i∈[n]

so if we just precompute Y = n1

P
i∈[n]

yi , then n1P2e−||x−yi ||2 can be ε-approximated in O(D)

i∈[n]

time.
The time to compute the mappings of all our points is O(nD), and then the time tocompute the inner products is also O(nD), so the final runtime is


nd
M
O
log
ε2
ε

12

References[Che18]

Lijie Chen. On the hardness of approximate and exact (bichromatic) maximum
inner product. arXiv preprint arXiv:1802.02325, 2018.

[GS91]

Leslie Greengard and John Strain. The fast gauss transform. SIAM Journal on
Scientific and Statistical Computing, 12(1):79–94, 1991.

[GW13]

Michael Griebel and Daniel Wissel. Fast approximation of the discrete gausstransform in higher dimensions. Journal of Scientific Computing, 55(1):149–172,
2013.

[RR08]Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In Advances in neural information processing systems, pages 1177–1184,
2008.

[RYDG05] Vikas Chandrakant Raykar, Changjaing Yang, Ramani Duraiswami, and Nail
Gumerov. Fast computation of sums of gaussians in high dimensions. Technicalreport, 2005.
[TW09]

Johannes Tausch and Alexander Weckiewicz. Multidimensional fast gausstransforms by chebyshev expansions. SIAM Journal on Scientific Computing,31(5):3547–3565, 2009.

[Wil18]Ryan Williams. On the difference between closest, furthest, and orthogonal pairs:
Nearly-linear vs barely-subquadratic complexity. In Proceedings of the TwentyNinth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 1207–1215.
SIAM, 2018.

13