Tradeoffs for nearest neighbors on the sphere
Thijs Laarhoven∗
September 13, 2016

arXiv:1511.07527v2 [cs.DS] 10 Sep 2016

Abstract
We consider tradeoffs between the query and update complexities for the (approximate) nearest neighbor problem on the sphere, extending the spherical filters recently introduced by [Becker–Ducas–Gama–
Laarhoven, SODA’16] to sparse regimes and generalizing the scheme and analysis to account for different
tradeoffs. In a nutshell, for the sparse regime the tradeoff between the query complexity nρq and update complexity nρu for data sets of size n can be summarized by the following equation in terms of the
approximation factor c and the exponents ρq and ρu :
p
√
√
c2 ρq + (c2 − 1) ρu = 2c2 − 1.
For small c = 1 + ε, minimizing the time for updates leads to a linear space complexity at the cost
2
of a query time complexity of approximately n1−4ε . Balancing the query and update costs leads to op2
timal complexities of n1/(2c −1) , matching lower bounds from [Andoni–Razenshteyn, 2015] and [Dubiner,
IEEE Trans. Inf. Theory 2010] and matching the asymptotic complexities previously obtained by [Andoni–
Razenshteyn, STOC’15] and [Andoni–Indyk–Laarhoven–Razenshteyn–Schmidt, NIPS’15]. A subpolyno2
mial query time complexity no(1) can be achieved at the cost of a space complexity of the order n1/(4ε ) ,
Ω(1/ε2 )
matching the lower bound n
of [Andoni–Indyk–Pǎtraşcu, FOCS’06] and [Panigrahy–Talwar–Wieder,
FOCS’10] and improving upon results of [Indyk–Motwani, STOC’98] and [Kushilevitz–Ostrovsky–Rabani,
STOC’98] with a considerably smaller leading constant in the exponent.
2
4
For large c, minimizing the update complexity results in a query complexity of n2/c +O(1/c ) , improving
upon the related asymptotic exponent for large c of [Kapralov, PODS’15] by a factor 2, and matching
2
the lower bound nΩ(1/c ) of [Panigrahy–Talwar–Wieder, FOCS’08]. Balancing the costs leads to optimal
2
complexities of the order n1/(2c −1) , while a minimum query time complexity can be achieved with up2
4
2
4
date and space complexities of approximately n2/c +O(1/c ) and n1+2/c +O(1/c ) , also improving upon the
previous best exponents of Kapralov by a factor 2 for large n and c.
For the regime where n is exponential in the dimension, we obtain further improvements compared to
results obtained with locality-sensitive hashing. We provide explicit expressions for the query and update
complexities in terms of the approximation factor c and the chosen tradeoff, and we derive asymptotic
results for the case of the highest possible density for random data sets.

1

Introduction

Approximate nearest neighbors (ANN). A central computational problem in many areas of research,
such as machine learning, coding theory, pattern recognition, data compression, and cryptanalysis [Bis06,
Dub10, DHS00, Her15, Laa15, MO15, SDI05], is the nearest neighbor problem: given a d-dimensional data
set D ⊂ Rd of cardinality n, design a data structure and preprocess D in a way that, when later given a
query vector q ∈ Rd , we can quickly find a nearest vector to q in D. A common relaxation of this problem
is the approximate nearest neighbor problem (ANN): given that the nearest neighbor in D lies at distance
at most r from q, design an efficient algorithm that finds an element p ∈ D at distance at most c · r from
q, for a given approximation factor c > 1. We will consider the case where d scales with n; for fixed d it
is well-known that one can answer queries in timenρ with ρ = o(1) with only a polynomial increase in the
space complexity [AMN+ 94].
∗

Eindhoven University of Technology, Eindhoven, The Netherlands. E-mail: mail@thijs.com

1

ANN on the sphere. Depending on the notion of distance, different solutions have been proposed in
the literature (e.g. [AHL01, DIIM04, GIM99, LRU14, WSSJ14]). In this work we will restrict out attention
to the angular distance, where two vectors are considered nearby iff their common angle is small [AIL+ 15,
Cha02, Laa15, SSLM14, STS+ 13]. This equivalently corresponds to spherical ANN under the `2 -norm,
where the entire data set is assumed to lie on the Euclidean unit sphere. Recent work of Andoni and
Razenshteyn [AR15a] showed how to reduce ANN in the entire Euclidean space to ANN on the sphere,
which further motivates why finding optimal solutions for the spherical case is relevant. For spherical, lowdensity settings (n = 2o(d) ), we will further focus
√ on the random setting described in [AR15a], where nearby
corresponds to a Euclidean distance of r = 1c 2 on the unit sphere (i.e. an angle θ = arccos(1 − 1/c2 )) and
√
far away corresponds to a distance c · r = 2 (angle ψ = 21 π).
Spherical locality-sensitive hashing. A well-known method for solving ANN in high dimensions is
locality-sensitive hashing (LSH) [IM98]. Using locality-sensitive hash functions, with the property that
nearby vectors are more likely to be mapped to the same output value than distant pairs of vectors, one
builds several hash tables with buckets containing sets of vectors with the same hash value. To answer a
query q, one computes q’s hash values and checks the corresponding buckets in each of the hash tables for
potential near neighbors. For spherical ANN in the random setting, two recent works [AR15a, AIL+ 15] have
shown how to solve ANN with query time Õ(nρ ) and space Õ(n1+ρ ) with ρ = 2c21−1 + o(1) where o(1) → 0
as n → ∞. For large c and n, this improved upon e.g. hyperplane LSH [Cha02] and Euclidean LSH [AI06].
Within the class of LSH algorithms, these results are known to be essentially optimal [AR15b, Dub10].
Spherical locality-sensitive filters. Recently Becker–Ducas–Gama–Laarhoven [BDGL16] introduced
spherical filters, which map the data set D to a subset D0 ⊆ D consisting of all points lying in a certain
spherical cap. Filtering could be considered a relaxation of locality-sensitive hashing: for LSH a hash
function is required to partition the space in regions, while for LSF this is not necessarily the case. Similar
filtering constructions were previously proposed in [Dub10, MO15]. For dense data sets of size n = 2Θ(d) ,
the approach of [BDGL16] led to a query exponent ρ < 1/(2c2 − 1) for the random setting.
Asymmetric ANN. The exponents ρ described so far are all for balanced or symmetric ANN: both the
time to answer a query and the time to insert/delete vectors from the data structure are then equal to
Õ(nρ ), and the time complexity for preprocessing the data and the total space complexity are both equal to
Õ(n1+ρ ). Depending on the application however, it may be desirable to obtain a different tradeoff between
these costs. In some cases it may be beneficial to use even more space and more time for the preprocessing,
so that queries can be answered even faster. In other cases, memory constraints might rule out the use of
balanced parameters, in which case one has to settle for a lower space and update complexity and it would
be interesting to know the best time complexity that can be achieved for a given space complexity. Finding
optimal tradeoffs between the different costs of ANN is therefore essential for achieving the best performance
in different contexts.
Smooth tradeoffs for asymmetric ANN. Various works have analyzed tradeoffs for ANN, among
others using multi-probing in LSH to reduce the memory complexity at the cost of a higher query complexity [AIL+ 15, AMM09, LJW+ 07, Pan06]. However, most existing techniques either describe one particular
tradeoff between the costs, or do not offer provable asymptotic bounds on the query and update exponent as
the parameters increase. Recently Kapralov [Kap15] showed how to obtain smooth and provable asymptotic
tradeoffs for Euclidean ANN, but as the exponents for the balanced setting are a factor 2 above the lower
bound ρ ≥ 1/(2c2 −1) for large c, it may be possible to improve upon these techniques not only for symmetric
but also for asymmetric ANN.

2

1.1

Contributions.

In this work we extend the symmetric ANN technique of [BDGL16] for dense regimes to asymmetric ANN
on the sphere for both sparse and dense regimes, showing how to obtain smooth and significantly improved
tradeoffs between the query and update complexities compared to e.g. [Kap15, Pan06] in both the small
c and large c regimes. For sparse settings, the tradeoff between the query complexity nρq and the update
complexity nρu can essentially be summarized by the non-negative solution pairs (ρu , ρq ) ∈ R2 to the following
equation, which can be expressed either in terms of θ (left) or c (right) by substituting cos θ = 1 − 1/c2 .
p
√
√
√
√
ρq + (cos θ) ρu = sin θ,
c2 ρq + (c2 − 1) ρu = 2c2 − 1.
(1)
The space complexity for the preprocessed data, as well as the time for preprocessing the data, are both
Õ(n1+ρu ). The resulting tradeoffs for the random case for small and large c are illustrated in Table 1 and
Figure 2, and can be derived from (1) by substituting ρu = 0 (minimize the space), ρq = ρu (balanced ANN),
or ρq = 0 (minimize the time), and computing Taylor expansions around c ∈ {1, ∞}.
Small approximation factors. In the regime of small c = 1 + ε, as described in Table 1 we obtain an
2
3
update complexity no(1) and space complexity n1+o(1) with query complexity n1−4ε +O(ε ) . This improves
upon results of Kapralov, where a sublinear query complexity and a quasi-linear space complexity could only
be achieved for approximation factors c > 1.73 [Kap15]. Balancing the complexities leads to asymptotic
exponents ρq = ρu = 1/(2c2 − 1), which means that both exponents scale as 1 − O(ε) for small c > 1.
These exponents match the asymptotic complexities previously obtained by [AR15a, AIL+ 15] and the lower
bounds from [AR15b, Dub10]. A query complexity no(1) can be achieved for arbitrary c with an update
2
complexity n4/ε +O(1/ε) , matching the asymptotic lower bounds of [AIP06, PTW10]1 and the constructions
of [IM98, KOR00] with a smaller leading constant in the exponent2 . This also improves upon [Kap15],
achieving a query complexity no(1) only for c > 1.73.
Large approximation factors. For large c, both ρq and ρu are proportional to 1/c2 , with leading constants 1/2 in the balanced regime, and leading constant 2 if the other complexity is minimized. This improves
upon results from [Kap15], whose leading constants are a factor 2 higher in all cases, and matches the lower
2
bound on the space complexity of nΩ(1/c ) of [PTW08] for query complexities no(1) .
High-density regime. Finally, for data sets of size n = 2Θ(d) , we obtain improved tradeoffs between
the query and update complexities compared to results obtained using locality-sensitive hashing, even for
balanced settings. We show that also for this harder problem we obtain query exponents less than 1 regardless
of the tradeoff, while a query exponent 0 is impossible to achieve with our methods.

1.2

Outline.

In Section 2 we describe preliminary notation and results regarding spherical filters. Section 3 describes
asymptotic tradeoffs for the sparse regime of n = 2o(d) , and Section 4 then discusses the application of these
techniques to the dense regime of n = 2Θ(d) . Section 5 concludes with a discussion on extending our methods
to slightly different problems, and open problems for future work.

2

Preliminaries

2.1

Subsets of the unit sphere.

We first recall some preliminary notation and results on geometric objects on the unit sphere, similar
to [BDGL16]. Let µ denote the canonical Lebesgue measure over Rd , and let us write h·, ·i for the standard
1
2

The constant 1/4 in the exponent even matches the lower bound for single-probe schemes of [PTW10, Theorem 1.5].
Explicit leading constants are not stated in [IM98, KOR00] but appear to be 9 and 9/2 respectively, compared to our 1/4.

3

General expressions

Small c = 1 + ε

Large c → ∞

Minimize space
(β = cos θ)

ρq = (2c2 − 1)/c4
ρu = 0

ρq = 1 − 4ε2 + O(ε3 )
ρu = 0

ρq = 2/c2 + O(1/c4 )
ρu = 0

Balance costs
(β = 1)

ρq = 1/(2c2 − 1)
ρu = 1/(2c2 − 1)

ρq = 1 − 4ε + O(ε2 )
ρu = 1 − 4ε + O(ε2 )

ρq = 1/(2c2 ) + O(1/c4 )
ρu = 1/(2c2 ) + O(1/c4 )

Minimize time
(β = 1/ cos θ)

ρq = 0
ρu = (2c2 − 1)/(c2 − 1)2

ρq = 0
ρu = 1/(4ε2 ) + O(1/ε)

ρq = 0
ρu = 2/c2 + O(1/c4 )

Table 1: The extreme points of our asymptotic tradeoffs. Answering a query takes time Õ(nρq ), updates
take Õ(nρu ) operations, and the space/preprocessing complexities are Õ(n1+ρu ). Lower order terms which
tend to 0 as d, n → ∞ are omitted for clarity. The colors match those used in Figure 2.
Euclidean inner product. We denote the unit sphere in Rd by S d−1 = {x ∈ Rd : kxk = 1} and half-spaces
by Hu,α := {x ∈ Rd : hu, xi ≥ α}. For constants α, α1 , α2 ∈ (0, 1) and vectors u, u1 , u2 ∈ S d−1 we denote
spherical caps and wedges by Cu,α := S d−1 ∩ Hu,α and Wu1 ,α1 ,u2 ,α2 := S d−1 ∩ Hu1 ,α1 ∩ Hu2 ,α2 .
For analyzing the performance of spherical filters, we would like to know the volumes of these objects
in high dimensions. The following asymptotic estimates can be found in [BDGL16, MV10], where γ =
γ(α1 , α2 , θ) satisfies γ 2 = (α12 + α22 − 2α1 α2 cos θ)/ sin2 θ and θ denotes the angle φ(u1 , u2 ) := arccos hu1 , u2 i
between u1 , u2 .
C(α) :=

2.2

d
p
µ(Cu,α )
Θ(1)
2
1
−
α
=
d
,
µ(S d−1 )

W(α1 , α2 , θ) :=

d
p
µ(Wu1 ,α1 ,u2 ,α2 )
Θ(1)
2
1
−
γ
=
d
.
µ(S d−1 )

(2)

Symmetric spherical filters in the dense regime.

We now continue with a brief description of the algorithm of Becker–Ducas–Gama–Laarhoven [BDGL16] for
solving dense ANN on the sphere.
Initialization. Let m = Θ(log d) and suppose that m|d. We partition the d coordinates into m blocks
of size d/m, and for each of these m blocks of coordinates we randomly sample t1/m code words from
S d/m−1 . This results in m subcodes C1 , . . . , Cm ⊂ S d/m−1 . Combining one code word from each subcode,
we obtain (t1/m )m = t different vectors √1m (c1 , . . . , cm ) ∈ S d−1 with ci ∈ Ci . We denote the resulting set
of vectors by the code C. The ideas behind this construction are that (1) this code C behaves as a set of
t random unit vectors in S d−1 , where the difference with a completely random code is negligible for large
parameters [BDGL16, Theorem 5.1]; and (2) the additional structure hidden in C allows us to decode faster
than with a linear search. The parameter t will be specified later.
Preprocessing. Next, given D ⊂ S d−1 , we consider each point p ∈ D one by one and compute its relevant
filters Update(p) := {c ∈ C : hp, ci ≥ α}. Naively finding these filters by a linear search over all filters would
cost time Õ(t), but as described in [BDGL16] this can be done in time O(|Update(pi )|) due to the hidden
additional structure in the code3 . Finally, we store all vectors in the respective filter buckets B1 , . . . , Bt ,
where p is stored in Bj iff cj ∈ Update(p). The parameter α will be specified later.
Answering a query. To find neighbors for a query vector q, we compute its relevant filters Query(q) :=
{c ∈ C : hq, ci ≥ α} in time proportional to the size of this set. Then, we visit all these buckets in our
data structure, and compare q to all vectors p in these buckets. The cost of this step is proportional to the
3

Note that the overhead of the enumeration-based decoding algorithm [BDGL16, Algorithm 1] mainly consists of computing
and sorting all blockwise inner products hpi , ci i, which can be done in time Õ(t1/m ) = no(1) .

4

number of vectors colliding with q in these filters, and the success probability of answering a query depends
on the probability that two nearby vectors are found due to a collision.
Updating the data structure (optional). In certain applications, it may further be important that
one can efficiently update the data structure when D is changed. Inserting or removing a vector p from
the buckets is done by computing Update(p) in time proportional to |Update(p)|, and inserting/removing
the vector from the corresponding buckets. Note that by e.g. keeping buckets sorted in lexicographic order,
updates in one bucket can be done in time Õ(log n) = dO(1) .
Correctness. To prove that this filtering construction works for certain parameters α and t, two properties
are crucial: the code C needs to be efficiently decodable, and C must be sufficiently smooth on S d−1 in the
sense that collision probabilities are essentially equal to those of uniformly random codes C ⊂ S d−1 . These
two properties were proved in [BDGL16, Lemma 5.1 and Theorem 5.1] respectively.

3

Asymmetric spherical filters for sparse regimes

To convert the spherical filter construction described above to the low-density regime, we need to make sure
that the overhead remains negligible in n. Note that costs t1/m = t1/ log d are considered no(1) in [BDGL16]
as t = 2Θ(d) and n = 2Θ(d) . In the sparse setting of n = 2Θ(d/ log d) , this may no longer be the case4 . To
overcome this potential issue, we set m = O(log2 d), so that t1/m = no(1) even if t = 2Θ(d) . Increasing m
means that the code C becomes less smooth, but a detailed inspection of the proof of [BDGL16, Thm. 5.1]
shows that also for m = logO(1) d the code is sufficiently smooth on S d−1 .
To allow for tradeoffs between the query/update complexities, we introduce two parameters αq and αu
for querying and updating the database. This means that we redefine Query(q) := {c ∈ C : hq, ci ≥ αq } and
Update(p) := {c ∈ C : hp, ci ≥ αu } where αq , αu ∈ (0, 1) are to be chosen later. Smaller parameters mean
that more filters are contained in these sets, so intuitively αq < αu means more time is spent on queries,
while αq > αu means more time is spent on updates and less time on queries.

3.1

Random sparse instances.

For studying the sparse regime of n = 2o(d) , we will consider the random model of [AR15a, AIL+ 15] defined
as follows.
Definition 1 (Random θ-ANN in the sparse regime). Given an angle θ ∈ (0, 21 π), a query q ∈ Rd , and a
data set D of size n = 2o(d) , the θ-ANN problem is defined as the problem of either finding a point p ∈ D
with φ(p, q) ≤ 21 π, or concluding that w.h.p. no vector p ∈ D exists with φ(p, q) ≤ θ.
Note that for angles ψ = 21 π − δ where δ > 0 is fixed independently of d and n, a random point u on
the sphere covers a fraction (1 − O(δ 2 ))d = 2−Θ(d) of the sphere with points at angle at most ψ from u.
Together, n = 2o(d) points therefore cover a fraction of at most t · 2−Θ(d) = 2−Θ(d) of the sphere. For a query
q sampled uniformly at random from the sphere, with high probability it is far away from all n points, i.e.
at angle at least ψ from D. In other words, we expect that there is a significant gap between the angle with
the (planted) nearest neighbor, and the angle with all other vectors, in which case solving ANN with small
approximation factor is actually equivalent to solving the exact NN problem.
If we were to define the notion of “far away” as being at some angle ψ < 12 π from q, and we somehow
expect that a significant part of the data set lies at angle at most ψ from q, then the data set and queries are
apparently concentrated on one part of the sphere and their distribution is not spherically symmetric. If this
is indeed the case, then this knowledge may be exploited using data-dependent ANN techniques [WSSJ14],
and such techniques may be preferred over data-independent filters.
4

For even sparser data sets, we can always first apply a dimension reduction using the Johnson-Lindenstrauss transform [JL84]
0
0
to transform the points to d0 -dimensional vectors with n = 2Ω(d / log d ) , and without significantly distorting inter-point distances.

5

p

•

Pru∼S d−1 [u ∈ Update(p)]
∝ C(αu )

αu

φ
0

Pru∼S d−1 [u ∈ Update(p) ∩ Query(q)]
∝ W(αq , αu , φ)

•q

•

αq

S d−1

Pru∼S d−1 [u ∈ Query(q)]
∝ C(αq )

Figure 1: The geometry of spherical filters. A vector p is inserted into/deleted from a filter u with probability
proportional to C(αu ), over the randomness of sampling u at random from S d−1 ; a filter u is queried for
nearest neighbors for q with probability C(αq ); and a vector p at angle φ from q is found as a candidate
nearest neighbor in one of the filters with probability proportional to W(αq , αu , φ).

3.2

Main result.

Before describing the analysis that leads to the optimized parameter choices, we state the main result for
the random, sparse setting described above, in terms of the nearby angle θ.
p
Theorem 1. Let θ ∈ (0, 21 π) and let β ∈ [cos θ, 1/ cos θ]. Then using parameters αq = β (2 log n)/d and
p
αu = (2 log n)/d we can solve the θ-ANN problem on the sphere with query/update exponents:

ρq =

1 − β cos θ
sin θ

2


+O

1
log d




,

ρu =

β − cos θ
sin θ

2


+O

1
log d


.

(3)

The resulting algorithm has a query time complexity Õ(nρq ), an update time complexity Õ(nρu ), a preprocessing time complexity Õ(n1+ρq ), and a total space complexity of Õ(n1+ρq ).
This result can equivalently be expressed in terms of c, by replacing θ = arccos(1 − 1/c2 ). In that case,
2
c2
θ ∈ (0, 21 π) translates to c ∈ (1, ∞), the interval for β becomes β ∈ [ c c−1
2 , c2 −1 ], and we get
(β(1 − c2 ) + c2 )2
ρq =
+O
2c2 − 1



1
log d

(1 − c2 + βc2 )2
ρu =
+O
2c2 − 1


,



1
log d


.

(4)

Due to the simple dependence of these expressions on β, we can easily compute β as a function of ρu and θ
(or c), and substitute this expression for β into ρq to express ρq in terms of ρu and θ (or c):





1
1 p 2
1
√
√
√
2
ρq = sin θ − ρu · cos θ + O
= 2
2c − 1 − ρu · (c − 1) + O
.
(5)
log d
c
log d
From these expressions we can derive both Table 1 and Figure 2 by substituting appropriate values for
ρq , ρu , β and computing Taylor expansions around c = 1 and c = ∞.

3.3

Cost analysis.

We now proceed with a proof of Theorem 1, by analyzing the costs of different steps of the filtering process
in terms of the spherical cap heights αq and αu and the angle of nearby vectors θ, and optimizing the
parameters accordingly. The analysis will be done in terms of θ.

6

Updating the data structure. The probability that a filter is considered for updates is equal to the
probability that hp, ci ≥ αu for random c, which is proportional to C(αu ) (cf. Figure 1). The size of
Update(p) ⊆ C and the time required to compute this set with efficient decoding [BDGL16, Algorithm 1]
are of the order t · C(αu ). The total preprocessing time comes down to repeating this procedure n times,
and the total space complexity is also equal to n · t · C(αu ). (We only store non-empty buckets.)
Answering a query. The probability that a filter is considered for query q is of the order C(αq ) (cf.
Figure 1), and the size of Query(p) is of the order t · C(αq ). After finding the relevant buckets, we go through
a number of collisions with distant vectors before (potentially) finding a near neighbor. The probability that
distant vectors collide in a filter is proportional to W(αq , αu , 12 π) (cf. Figure 1), so the number of comparisons
for all t filters and all n distant vectors is Õ(n · t · W(αq , αu , 21 π)).
Choosing the number of filters. Note that the probability that a nearby vector at angle at most θ
from a query q collides with q in a random filter is proportional to W(αq , αu , θ). By the union bound, the
probability that two nearby vectors collide in at least one filter is at least t · W(αq , αu , θ). To make sure that
nearby vectors are found with constant probability (say 90%), we set t ∝ 1/W(αq , αu , θ). With the choice
of t fixed in terms of αq , αu , θ and the above cost analysis in mind, the following table gives an overview of
the asymptotic costs of spherical filtering in random, sparse settings.

3.4

Quantity

Costs for general αq , αu , θ

Time: Finding relevant filters for a query
Time: Comparing a query with colliding vectors
Time: Finding relevant filters for an update
Time: Preprocessing the data
Space: Storing all filter entries

C(αq ) / W(αq , αu , θ)
n · W(αq , αu , 12 π) / W(αq , αu , θ)
C(αu ) / W(αq , αu , θ)
n · C(αu ) / W(αq , αu , θ)
n · C(αu ) / W(αq , αu , θ)

Balancing the query costs.

Next, note that answering a query consists of two steps: find the αq -relevant filters, and go through all
candidate near neighbors in these buckets. To obtain an optimal balance between these costs (with only a
polynomial difference in d in the time complexities), we must have C(αq ) = dO(1) · n · W(αq , αu , 12 π). Raising
both sides to the power 2/d, this is equivalent to 1 − αu2 = dO(1/d) n2/d (1 − αq2 − αu2 ). Isolating αu and noting
that n1/d = exp O( log1 d ), this leads to:
s
r



n2/d − 1
2 log n
1
O(1/d)
=
αu = d
1
+
O
.
(6)
d
log d
n2/d
This choice of αu guarantees that the query costs are balanced. As αq will have a similar scaling to αu , we
set αq = β · αu for β to be chosen later. Note that αq , αu = o(1) implies that the corresponding spherical
caps (cf. Figure 1) are almost-hemispheres,
similar to spherical LSH
[AR15a]. However, in our case the
√
√
4
parameters scale as αq , αu = O(1/ log d), compared to α = O(1/ d) in [AR15a].
Explicit costs. With αu fixed and the relation between αq and αu expressed in terms of β, we now
evaluate the costs for large d and n = exp O( logd d ) in terms of β and θ. Using Taylor expansions we get:




log C(αq )
log W(αq , αu , θ)
1
1 + β 2 − 2β cos θ
1
2
= −β + O
,
=−
+O
,
(7)
log n
log d
log n
log d
sin2 θ




log W(αq , αu , 12 π)
1
1
log C(αu )
2
= −1 + O
,
= −1 − β + O
.
(8)
log n
log d
log n
log d
Combining these expressions, we can derive asymptotics for all of the costs related to the filtering algorithm.
For the query/update exponents we then obtain the expressions given in Theorem 1.
7

3.5

Optimal parameter range.

Note that the best complexities are obtained by choosing β ∈ [cos θ, 1/ cos θ]; beyond this range, complexities
are strictly worse. Taking inverses in this range, we get:




√
1 − ρq · sin θ
1
1
√
= cos θ + ρu · sin θ + O
.
(9)
β=
+O
cos θ
log d
log d
Isolating

4

√

ρq then leads to (1), while (9) also shows how to choose β to achieve given complexities.

Asymmetric spherical filters for dense regimes

We now revisit the dense regime of data sets of size n = 2Θ(d) , as previously analyzed in [BDGL16] for
symmetric ANN. We will again use two parameters αq , αu ∈ [0, 1] where the optimization now leads to a
slightly different, more refined result, depending on the chosen tradeoff.

4.1

Random dense instances.

To study the dense regime, we consider the following model.
Definition 2 (Random θ-ANN in the dense regime). Given an angle θ ∈ (0, 21 π), a query q ∈ Rd , and a data
set D of n = 2Θ(d) points sampled uniformly at random from S d−1 , the random θ-ANN problem is defined
as the problem of either finding a point p ∈ D with φ(p, q) ≤ θ, or concluding that with high probability no
vector p ∈ D exists with φ(p, q) ≤ θ.
At first sight, the above definition does not seem to correspond to an approximate, but to an exact nearest
neighbor instance. However, we made a critical assumption on D here: we assume that these points are
sampled uniformly at random from the sphere. This seems to be a natural assumption in various applications
(see e.g. [Her15, Laa15, MO15]), and this implies that in fact many of the points in the data set lie at angle
approximately 21 π from q. As a result the problem is significantly easier than e.g. the worst-case ANN
setting with c ≈ 1, where the entire data set might lie at angle θ + δ from q.
For comparing this problem with the sparse setting of Section 3, observe that this dense problem is
harder : uniformly random points on the sphere (of which roughly half has angle less than 21 π with q) are
more likely to cause collisions than orthogonal points to q. The number of collisions with distant vectors
will therefore increase, and we expect the query and update exponents to be larger. This was also observed
in e.g. [BDGL16, BL15, Laa15, LdW15], where the exponent for lattice sieving with other ANN techniques
would have been smaller if one could only assume that far away means orthogonal.
Note that we could also extend the analysis of Section 3 to the dense regime simply by fixing the distant
angle at ψ = π2 . In that case, similar to [BDGL16], both the query and update exponents will become
smaller compared to the low-density regime as the problem becomes easier. However, such an assumption
would imply that the data set is not spherically symmetric and is concentrated on one part of the sphere,
in which case data-dependent methods may be preferred [LRU14, WSSJ14].

4.2

Density reduction and the critical density.

As D is sampled at random from S d−1 , a point p ∈ D is close to q with probability proportional to (sin θ)d .
With n points, we expect to find approximately n · (sin θ)d nearby neighbors p ∈ D. For n  (sin θ)−d ,
nearby vectors are rare, and we are essentially solving the exact (decisional) nearest neighbor problem with
high probability. On the other hand, if n  (sin θ)−d , then we expect there to be many (nO(1) ) solutions
p ∈ D with φ(p, q) ≤ θ.
In our analysis we will focus on the case where n = Õ((sin θ)−d ): there might not be any near neighbors
in D at all, but if there is one, we want to find it. For the regime n  (sin θ)−d , we can reduce this problem
to a regime with a lower density n0 = Õ((sin θ)−d ) through a simple transformation:
8

• Randomly select a subset D0 ⊂ D of size n0 = Õ((sin θ)−d ).
• Run the (approximate) nearest neighbor algorithm on this subset D0 .
By choosing the hidden factor inside n0 sufficiently large, with high probability there will be a solution in
this smaller subset D0 as well, which our algorithm will find with high probability. This means that in our
cost analysis, we can then simply replace n by n0 to obtain the asymptotic complexities after this density
reduction step. We denote the regime n ∝ (sin θ)−d as the critical density.
Note that if we are given a random data set of size n, then we expect the nearest neighbor to a random
query q ∈ S d−1 to lie at angle θ ≈ arcsin(n−1/d ) from q; for larger angles we will find many solutions, while
for smaller angles w.h.p. there are no solutions at all (except for planted near neighbor instances (outliers)
which are not random on the sphere). This further motivates why the critical density is important, as
θ ≈ arcsin(n−1/d ) commonly corresponds to solving the exact nearest neighbor problem for random data
sets. Setting θ  arcsin(n−1/d ) then corresponds to searching for outliers.

4.3

Main result.

We first state the main result for the random, dense setting described above without making any assumptions
on the density. A derivation of Theorem 2 can be found in Appendix B.
p
Theorem 2. Let θ ∈ (0, 21 π) and let β ∈ [cos θ, 1/ cos θ]. Then using parameters αq = β 1 − n−2/d and
p
αu = 1 − n−2/d we can solve the dense θ-ANN problem on the sphere with exponents:



 i

h

2
d
−d
−2/d 1 + β − 2β cos θ
−2/d
+
(10)
ρq =
β2 ,
log 1 − 1 − n
log
1
−
1
−
n
2 log n
2 log n
sin2 θ




2
−d
−2/d 1 + β − 2β cos θ
ρu =
log 1 − 1 − n
− 1.
(11)
2 log n
sin2 θ
Note that in the limit of n1/d → 1, corresponding to sparse data sets, we obtain the expressions from
Theorem 1. This matches our intuition that taking n points uniformly at random from the sphere with
n = 2o(d) roughly means that all points have angle ψ = 12 π with a query q, as described in Section 3.1.

4.4

Critical densities.

As a special case of the above result, we focus on the regime of n ∝ (sin θ)−d . The following result shows
the complexities obtained after substituting this density into Theorem 2.
Corollary 1. Let θ ∈ (0, 12 π), let β ∈ [cos θ, 1/ cos θ], and let n = (1/ sin θ)d . Then using parameters
αq = β cos θ and αu = cos θ, the complexities in Theorem 2 reduce to:
n

ρq


=

sin2 θ (β cos θ + 1)
β cos θ − cos 2θ

d/2
,

n

ρu


=

sin2 θ
1 − cot2 θ (β 2 − 2β cos θ + 1)

d/2
.

(12)

To obtain further intuition into the above results, let us consider the limits obtained by setting β ∈
{cos θ, 1, 1/ cos θ}. For β = cos θ, we obtain exponents of the order:
n ∝ (sin θ)−d ,

β = cos θ,

=⇒

ρu = 0,

ρq =

log 2 − log(3 + 2 cos(2 θ))
.
2 log(sin θ)

(13)

Balancing the complexities is done by setting β = 1, in which case we obtain the asymptotic expressions:
n ∝ (sin θ)−d ,

β = 1,

=⇒

ρq = ρu =

2 log(tan 2θ ) + log(2 cos θ + 1)
− 1.
2 log(sin θ)

(14)

To minimize the query complexity, we would ideally set β = 1/ cos θ, as then the query exponent approaches
0 for large n. However, one can easily verify that substituting β = 1/ cos θ into (12) leads to a denominator
9

of 0, i.e. the update costs and the space complexities blow up as β approaches 1/ cos θ. To get an idea of
how the update complexity scales in terms of the query complexity, we set ρq = δ and we compute a Taylor
expansion around δ = 0 for ρu to obtain:

log 8 δ log(1/ sin θ) tan2 θ
−d
n ∝ (sin θ) , β = 1
=⇒
ρq = δ, ρu =
− 1 + O(δ).
(15)
2 log(sin θ)
In other words, for fixed angles θ, to achieve ρq = δ the parameter ρu scales as log(1/δ). Note that except
for the latter result, we can substitute cos θ = 1 − 1/c2 and c = 1 + ε, and compute a Taylor series expansion
of these expressions around ε = 0 to obtain the expressions in Table 1 for small c = 1 + ε. This matches our
intuition that θ → 21 π for random data sets corresponds to θ ≈ 21 π for sparse settings.
Finally, we observe that substituting θ = 31 π, for a minimum space complexity we obtain ρq = log( 45 )/ log( 43 ),
while balancing both costs leads to ρq = ρu = log( 98 )/ log( 34 ). These results match those derived in [BDGL16]
for the application to sieving for finding shortest lattice vectors.

5

Discussion and open problems

We conclude this work with a brief discussion on how the described methods can possibly be extended and
modified to solve other problems and to obtain a better performance in practical applications.
Probing sequences. In LSH, a common technique to reduce the space complexity, or to reduce the number
of false negatives, is to use probing [AIL+ 15, LJW+ 07, Pan06]: one does not only check exact matches in
the hash tables for reductions, but also approximate matches which are still more likely to contain near
neighbors than random buckets. Efficiently being able to define a probing sequence of all buckets, in order
of likelihood of containing near neighbors, can be useful both in theory and in practice.
For spherical filters, an optimal probing sequence is obtained by sorting the filters (code words) according
to their inner products with the target vector. Due to the large number of buckets t, computing and sorting
all filter buckets is too costly, but for practical applications we can do the following. We first choose a
sequence 1 = α0 > α1 > · · · > αT , and then given a target t we apply our decoding algorithm to find all
code words c ∈ C with hc, ti ∈ (α1 , α0 ]. The corresponding buckets are the most likely to contain nearby
vectors. If this does not result in a nearest neighbor, we apply our decoding algorithm to find code words
c ∈ C with hc, ti ∈ (α2 , α1 ], and we repeat the above procedure until e.g. we are convinced that no solution
exists. For constant T , the overhead of this repeated decoding is small.
To implement the search for finding code words c ∈ C with hc, ti ∈ (αlow , αu ] efficiently, we can use
Algorithm 1 in Appendix C. The most costly step of this decoding algorithm is computing and sorting all
blockwise inner products hck,j , tk i, but note that these computations have to be performed only once; later
calls to this function with different intervals (αi+1 , αi ] can reuse these sorted lists.
Ruling out false negatives. An alternative to using probing sequences to make sure that there are no
false negatives, is to construct a scheme which guarantees that there are never any false negatives at all (see
e.g. [Pag16]). In the filtering framework, this corresponds to using codes C such that it is guaranteed that
nearby vectors always collide in one of the filters. In other words, for each pair of points p, q on the sphere
at angle θ, the corresponding wedge Wp,αu ,q,αq must contain a code word c ∈ C. Note that with our random
construction we can only show that with high probability, this is the case.
For spherical filters, codes guaranteeing this property correspond to spherical codes such that all possible
wedges Wp,αu ,q,αq for p, q ∈ S d−1 contain at least one code word c ∈ C. For αq = αu = α, note that at
the middle of such a wedge lies a point y = (p + q)/kp + qk at angle θ/2 from both p and q. If a code is
not covering and allows for false negatives, then there are no code words at angle 12 θ − arccos α from y. In
particular, the covering radius of the code (the smallest angle ψ such that spheres of angle ψ around all code
words cover the entire sphere) is therefore larger than 12 θ − arccos α. Equivalently, being able to construct

10

spherical codes of low cardinality with covering radius at most 12 θ − arccos α implies being able to construct
a spherical filtering scheme without false negatives.
As we make crucial use of concatenated codes C = C1 × · · · × Cm to allow for efficient decoding, covering
codes without efficient decoding algorithms cannot be used for C. Instead, one might aim at using such
α
covering codes for the subcodes: if all subcodes Ci have covering radius at most 12 θ − arccos m
, then the
1
concatenated code C = C1 ×· · ·×Cm has a covering radius of at most 2 θ −arccos α. Finding tight bounds on
α
the size of a spherical code with covering radius 12 θ −arccos m
(where θ is defined by the problem setting, and
α may be chosen later) would directly lead to an upper bound on the number of filters needed to guarantee
that there are no false positives.
Sparse codes for efficiency. As described in e.g. [Ach01, BDGL16, LHC06], it is sometimes possible
to use sparse, rather than fully random codes without losing on the performance of a nearest neighbor
scheme. Using sparse subcodes Ci might further reduce the overhead of decoding (computing blockwise
inner products). For this we could either use randomly sampled sparse subcodes, but one might also consider
using codes which are guaranteed to be “smooth” on the sphere and have a small covering radius. Similar
to Leech lattice LSH [AI06], one might consider using vectors from the Leech lattice in 24 dimensions to
define the subcodes. Asymptotically in our construction the block size d/m needs to scale with d, and fixing
d/m = 24 would invalidate the proof of smoothness of [BDGL16, Theorem 5.1], but in practice both n and
d are fixed, and only practical assessments can show whether predetermined spherical codes can be used to
obtain an even better performance.
Optimality of the tradeoff. An interesting open problem for future work is determining precise bounds
on the best tradeoff that one can possibly hope to achieve in the sparse regime. Since our tradeoff matches
various known bounds in the regimes of small and large approximation factors c [AIP06, And09, KOR00,
PTW08, PTW10] and no LSH scheme can improve upon our results in the balanced setting [AR15b, Dub10],
and since the tradeoff can be described through a remarkably simple relation (especially when described in
terms of θ), we conjecture that this tradeoff is optimal.
√
√
Conjecture 1. Any algorithm for sparse θ-ANN on the sphere must satisfy ρq + ρu · cos θ ≥ sin θ.
As a first step, one might try to (dis)prove this conjecture within the LSH framework, similar to various
other works focusing on lower bounds for schemes that fall in this category [MNP07, OWZ14].
Extension to Euclidean spaces. As mentioned in the introduction, Andoni and Razenshteyn showed
how to reduce ANN in Euclidean spaces to (sparse) random ANN on the sphere in the symmetric case,
using LSH techniques [AR15a]. An important open problem for future work is to see whether the techniques
and the reduction described in [AR15a] are compatible with locality-sensitive filters, and with asymmetric
nearest neighbor techniques such as those presented in this paper. If this is possible, then our results may
also be applicable to all of `d2 , rather than only to the angular distance on Rd or to Euclidean distances on
the unit sphere S d−1 .
Combination with cross-polytope LSH. Finally, the recent paper [AIL+ 15] showed how cross-polytope
hashing (previously introduced by Terasawa and Tanaka [TT07]) is asymptotically equally suitable for
solving Euclidean nearest neighbor problems on the sphere (and for the angular distance) as the approach of
spherical LSH of using large, completely random codes on the sphere [AR15a]. Advantages of cross-polytope
LSH over spherical LSH are that they have a much smaller size (allowing for faster decoding), and that
cross-polytope hash functions can be efficiently rerandomized using sparse and fast random projections such
as Fast Hadamard Transforms [AIL+ 15]. In that sense, cross-polytope LSH offers a significant practical
improvement over spherical LSH.
The approach of using spherical filters is very similar to spherical LSH: large, random (sub)codes are used
to define regions on the sphere. A natural question is therefore whether ideas analogous to cross-polytope
11

hashing can be used in combination with spherical filters, to reduce the subexponential overhead in d for
decoding to an overhead which is only polynomial in d. This is also left as an open problem for further
research.

References
[Ach01]

Dimitris Achlioptas. Database-friendly random projections. In PODS, pages 274–281, 2001.

[AHL01]

Helmut Alt and Laura Heinrich-Litan. Exact L∞ nearest neighbor search in high dimensions.
In SOCG, pages 157–163, 2001.

[AI06]

Alexandr Andoni and Piotr Indyk. Near-optimal hashing algorithms for approximate nearest
neighbor in high dimensions. In FOCS, pages 459–468, 2006.

[AIL+ 15]

Alexandr Andoni, Piotr Indyk, Thijs Laarhoven, Ilya Razenshteyn, and Ludwig Schmidt. Practical and optimal LSH for angular distance. In NIPS, 2015.

[AIP06]

Alexandr Andoni, Piotr Indyk, and Mihai Pǎtraşcu. On the optimality of the dimensionality
reduction method. In FOCS, pages 449–458, 2006.

[AMM09] Sunil Arya, Theocharis Malamatos, and David M. Mount. Space-time tradeoffs for approximate
nearest neighbor searching. Journal of the ACM, 57(1):1:1–1:54, November 2009.
[AMN+ 94] Sunil Arya, David M. Mount, Nathan S. Netanyahu, Ruth Silverman, and Angela Y. Wu. An
optimal algorithm for approximate nearest neighbor searching in fixed dimensions. In SODA,
pages 573–582, 1994.
[And09]

Alexandr Andoni. Nearest Neighbor Search: the Old, the New, and the Impossible. PhD thesis,
Massachusetts Institute of Technology, 2009.

[AR15a]

Alexandr Andoni and Ilya Razenshteyn. Optimal data-dependent hashing for approximate near
neighbors. In STOC, pages 793–801, 2015.

[AR15b]

Alexandr Andoni and Ilya Razenshteyn. Tight lower bounds for data-dependent locality-sensitive
hashing. Manuscript, pages 1–15, 2015.

[BDGL16] Anja Becker, Léo Ducas, Nicolas Gama, and Thijs Laarhoven. New directions in nearest neighbor
searching with applications to lattice sieving. In SODA, 2016.
[Bis06]

Christopher M. Bishop. Pattern Recognition and Machine Learning (Information Science and
Statistics). Springer-Verlag, 2006.

[BL15]

Anja Becker and Thijs Laarhoven. Efficient (ideal) lattice sieving using cross-polytope LSH.
Cryptology ePrint Archive, Report 2015/823, pages 1–25, 2015.

[Cha02]

Moses S. Charikar. Similarity estimation techniques from rounding algorithms. In STOC, pages
380–388, 2002.

[DHS00]

Richard O. Duda, Peter E. Hart, and David G. Stork. Pattern Classification (2nd Edition).
Wiley, 2000.

[DIIM04]

Mayur Datar, Nicole Immorlica, Piotr Indyk, and Vahab S. Mirrokni. Locality-sensitive hashing
scheme based on p-stable distributions. In SOCG, pages 253–262, 2004.

[Dub10]

Moshe Dubiner. Bucketing coding and information theory for the statistical high-dimensional
nearest-neighbor problem. IEEE Transactions on Information Theory, 56(8):4166–4179, Aug
2010.
12

[GIM99]

Aristides Gionis, Piotr Indyk, and Rajeev Motwani. Similarity search in high dimensions via
hashing. In VLDB, pages 518–529, 1999.

[Her15]

Gottfried Herold. Applications of nearest neighbor search techniques to the BKW algorithm
(draft). 2015.

[IM98]

Piotr Indyk and Rajeev Motwani. Approximate nearest neighbors: Towards removing the curse
of dimensionality. In STOC, pages 604–613, 1998.

[JL84]

William B. Johnson and Joram Lindenstrauss. Extensions of Lipschitz mappings into a Hilbert
space. Contemporary Mathematics, 26(1):189–206, 1984.

[Kap15]

Michael Kapralov. Smooth tradeoffs between insert and query complexity in nearest neighbor
search. In PODS, pages 329–342, 2015.

[KOR00]

Eyal Kushilevitz, Rafail Ostrovsky, and Yuval Rabani. Efficient search for approximate nearest
neighbor in high dimensional spaces. SIAM Journal on Computing, 30(2):457–474, 2000.

[Laa15]

Thijs Laarhoven. Sieving for shortest vectors in lattices using angular locality-sensitive hashing.
In CRYPTO, pages 3–22, 2015.

[LdW15]

Thijs Laarhoven and Benne de Weger. Faster sieving for shortest lattice vectors using spherical
locality-sensitive hashing. In LATINCRYPT, pages 101–118, 2015.

[LHC06]

Ping Li, Trevor J. Hastie, and Kenneth W. Church. Very sparse random projections. In KDD,
pages 287–296, 2006.

[LJW+ 07] Qin Lv, William Josephson, Zhe Wang, Moses Charikar, and Kai Li. Multi-probe LSH: efficient
indexing for high-dimensional similarity search. In VLDB, pages 950–961, 2007.
[LRU14]

Jure Leskovec, Anand Rajaraman, and Jeffrey D. Ullman. Mining of Massive Datasets. Cambridge University Press, 2014.

[MNP07]

Rajeev Motwani, Assaf Naor, and Rina Panigrahy. Lower bounds on locality sensitive hashing.
SIAM Journal of Discrete Mathematics, 21(4):930–935, 2007.

[MO15]

Alexander May and Ilya Ozerov. On computing nearest neighbors with applications to decoding
of binary linear codes. In EUROCRYPT, pages 203–228, 2015.

[MV10]

Daniele Micciancio and Panagiotis Voulgaris. Faster exponential time algorithms for the shortest
vector problem. In SODA, pages 1468–1480, 2010.

[OWZ14]

Ryan O’Donnell, Yi Wu, and Yuan Zhou. Optimal lower bounds for locality-sensitive hashing
(except when q is tiny). ACM Transactions on Computation Theory, 6(1):5:1–5:13, 2014.

[Pag16]

Rasmus Pagh. Locality-sensitive hashing without false negatives. In SODA, 2016.

[Pan06]

Rina Panigrahy. Entropy based nearest neighbor search in high dimensions. In SODA, pages
1186–1195, 2006.

[PTW08]

Rina Panigrahy, Kunal Talwar, and Udi Wieder. A geometric approach to lower bounds for
approximate near-neighbor search and partial match. In FOCS, pages 414–423, 2008.

[PTW10]

Rina Panigrahy, Kunal Talwar, and Udi Wieder. Lower bounds on near neighbor search via
metric expansion. In FOCS, pages 805–814, Oct 2010.

[SDI05]

Gregory Shakhnarovich, Trevor Darrell, and Piotr Indyk. Nearest-Neighbor Methods in Learning
and Vision: Theory and Practice. MIT Press, 2005.
13

1.0

)
c2
-1
(2
1/
=
ρu
=

0.4

c = 1.1

c = 1.2

c2

c = 1.3

ρq + (c 2 - 1)

ρu =

2 c2 - 1

1.6

1.8

ρq

0.6

c = 1.05

ρq = (2c2-1)/c4
ρu = 0

⟶ Query exponent ρq

0.8

c = 1.02

c = 1.5
0.2

c = 1.8

ρq = 0
ρu = (2c2 -1)/(c2 -1)2

c = 2.5
0.0
0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

2.0

⟶ Update exponent ρu
Figure 2: Tradeoffs between the query and update complexities, for various c. Informally, the x-axis represents the space and the y-axis the time (for queries). The diagonal represents symmetric ANN.
[SSLM14] Ludwig Schmidt, Matthew Sharifi, and Ignacio Lopez-Moreno. Large-scale speaker identification.
In ICASSP, pages 1650–1654, 2014.
[STS+ 13]

Narayanan Sundaram, Aizana Turmukhametova, Nadathur Satish, Todd Mostak, Piotr Indyk,
Samuel Madden, and Pradeep Dubey. Streaming similarity search over one billion tweets using
parallel locality-sensitive hashing. VLDB, 6(14):1930–1941, 2013.

[TT07]

Kengo Terasawa and Yuzuru Tanaka. Spherical LSH for approximate nearest neighbor search
on unit hypersphere. In WADS, pages 27–38, 2007.

[WSSJ14] Jingdong Wang, Heng Tao Shen, Jingkuan Song, and Jianqiu Ji. Hashing for similarity search:
A survey. arXiv:1408.2927 [cs.DS], pages 1–29, 2014.

A

Tradeoff figure in the sparse regime

Figure 2 describes asymptotic tradeoffs for different values of c. Note that the query exponent is always
smaller than 1, regardless of c > 1,pand note that the update exponent is smaller than 1 (and the query
√
exponent less than 1/2) for all c > 2 + 2 ≈ 1.85 corresponding to θ = 14 π.

B

Analysis for dense regimes

To derive the complexities for the dense regime of n = 2Θ(d) , we follow the same approach as for the sparse
regime of n = 2o(d) , but without making the assumption that e.g. n1/d = 1 + o(1).

14

B.1

General expressions.

First, we observe that the cost analysis described in Section 3 applies in the dense setting as well, with
the modification that we no longer assume that q is orthogonal to all of D (except for a potential nearest
neighbor). The update and query costs remain the same as before in terms of t and the volumes C(·) and
W(·), but to obtain the number of collisions with distant vectors, we take a different angle. First, observe
that each vector is added to Õ(t · C(αu )) filters, and that we have n vectors, leading to Õ(n · t · C(αu )) total
entries in the filters or Õ(n · t · C(αu )) entries in each filter5 . For a given vector q, we then query Õ(t · C(αq ))
buckets for nearest neighbors. In total, we therefore expect to find Õ(n · t · C(αu ) · C(αq )) colliding vectors for
a query vector. Again setting t ∝ 1/W(αq , αu , θ) to make sure that nearby vectors are found with constant
probability, we obtain the following updated table of the asymptotic costs of spherical filtering in random,
dense settings.

B.2

Quantity

Costs for αq , αu , θ

Time: Finding relevant filters for a query
Time: Comparing a query with colliding vectors
Time: Finding relevant filters for an update
Time: Preprocessing the data
Space: Storing all filter entries

C(αq ) / W(αq , αu , θ)
n · C(αq ) · C(αu ) / W(αq , αu , θ)
C(αu ) / W(αq , αu , θ)
n · C(αu ) / W(αq , αu , θ)
n · C(αu ) / W(αq , αu , θ)

Balancing the query costs.

Next, to make sure that the query costs are balanced, and not much more time is spent on looking for
relevant filters rather than actually doing comparisons, we again look for parameters such that these costs
are balanced. In this case we want to solve the asymptotic
equation C(αq ) = n · C(αq ) · C(αu ) or C(αu ) =
p
2
d/2
−2/d
(1 − αu )
= 1/n. Solving for αu leads to αu = 1 − n
leading to the parameter choice described in
Theorem 2. For now we again set αq = β · αu with β to be chosen later.

B.3

Explicit costs.

We now evaluate the costs for large d and n = 2Θ(d) , in terms of the ratio β between the two parameters αq
and αu , and the nearby angle θ. This leads to C(αu ) = 1/n and:


C(αq ) = 1 − (1 − n

−2/d

)β

2

d/2


,

W(αq , αu , θ) =

1 − (1 − n

−2/d

1 + β 2 − 2β cos θ
)
sin2 θ

d/2
.

(16)

Combining these expressions, we can then derive asymptotic estimates for all of the costs of the algorithm.
For the query and update exponents ρq = log[C(αq )/W(αq , αu , θ)]/ log n and ρu = log[C(αu )/W(αq , αu , θ)]/ log n
we then obtain:




h

 i
2
d
−d
−2/d 1 + β − 2β cos θ
−2/d
ρq =
log 1 − 1 − n
+
log
1
−
1
−
n
β2 ,
(17)
2 log n
2 log n
sin2 θ




2
−d
−2/d 1 + β − 2β cos θ
ρu =
log 1 − 1 − n
− 1.
(18)
2 log n
sin2 θ
These are also the expressions given in Theorem 2.

B.4

Optimal parameter range.

We observe that again the best exponents ρq and ρu are obtained by choosing β ∈ [cos θ, 1/ cos θ]; beyond
this range, the complexities are strictly worse. This completes the derivation of Theorem 2.
5

Here a crucial property of the filters is that each filter covers a region of exactly the same size on the sphere.

15

Algorithm 1 EfficientIntervalDecoding(C, t, αlow , αhigh )
Require: The description C1 , . . . , Cm of the code C; a target vector t ∈ Rd ; and 0 ≤ αlow < αhigh ≤ 1.
Ensure: Return all code words c ∈ C with ht, ci ∈ (αlow , αhigh ]
1: Sort each list Ck by decreasing dot-products dk,j = htk , ck,j i with tk .
Pm
2: Precompute m bounds Lk = αlow − Pi=k+1 di,t1/m .
m
3: Precompute m bounds Uk = αhigh − i=k+1 di,1 .
4: Initialize an empty output set S ← ∅.
5: Compute the lower bound `1 = min{j1 : d1,j1 > L1 }.
. do a binary search over [1, t1/m ]
6: Compute the upper bound u1 = max{j1 : d1,j1 ≤ U1 }.
. do a binary search over [1, t1/m ]
7: for each j1 ∈ {`1 , . . . , u1 } do
8:
Compute the lower bound `2 = min{j2 : d2,j2 > L2 − d1,j1 }.
9:
Compute the upper bound u2 = max{j2 : d2,j2 ≤ U2 − d1,j1 }.
10:
for each j2 ∈ {`2 , . . . , u2 } do
11:
[...]
Pm−1
12:
Compute the lower bound `m = min{jm : dm,jm > Lm − k=1 dk,jk }.
Pm−1
13:
Compute the upper bound um = max{jm : dm,jm ≤ Um − k=1 dk,jk }.
14:
for each jm ∈ {`m , . . . , um } do
15:
Add the code word c = (c1,j1 , . . . , cm,jm ) to S.
16:
end for
17:
[...]
18:
end for
19: end for
20: return S

C

Interval decoding

Algorithm 1 describes how to perform list-decoding for intervals, which may be relevant in practice for e.g.
computing probing sequences as described in Section 5. The algorithm is based on [BDGL16, Algorithm 1],
where now two sets of bounds are maintained to make sure that we only consider solutions which lie within
the given range, rather than above a threshold. The bounds Lk and Uk indicate the minimum and maximum
sum of inner products that can still be obtained in the last m
P− k sorted lists of vectors and inner products;
if in the nested for-loops, the current sum of inner products i di,ji is not in the interval (Lk , Uk ], then there
are no solutions anymore in the remaining part of the tree. Conversely, if this sum of inner products does
lie in the interval, then there must be at least one solution.

16

