Hypercube LSH for approximate near neighbors
Thijs Laarhoven

arXiv:1702.05760v1 [] 19 Feb 2017

IBM Research
Rüschlikon, Switzerland
mail@thijs.com

Abstract. A celebrated technique for finding near neighbors for the angular distance involves using a set
of random hyperplanes to partition the space into hash regions [Charikar, STOC 2002]. Experiments later
showed that using a set of orthogonal hyperplanes, thereby partitioning the space into the Voronoi regions
induced by a hypercube, leads to even better results [Terasawa and Tanaka, WADS 2007]. However, no
theoretical explanation for this improvement was ever given, and it remained unclear how the resulting
hypercube hash method scales in high dimensions.
In this work, we provide explicit asymptotics for the collision probabilities when using hypercubes to
partition the space. For instance, two near-orthogonal vectors are expected to collide with probability
( π1 )d+o(d) in dimension d, compared to ( 12 )d when using random hyperplanes. Vectors at angle π3 collide
√
with probability ( π3 )d+o(d) , compared to ( 32 )d for random hyperplanes, and near-parallel vectors collide
with similar asymptotic probabilities in both cases.
For c-approximate nearest neighbor searching, this translates to a decrease in the exponent ρ of localitysensitive hashing (LSH) methods of a factor up to log2 (π) ≈ 1.652 compared to hyperplane LSH. For c = 2,
we obtain ρ ≈ 0.302 + o(1) for hypercube LSH, improving upon the ρ ≈ 0.377 for hyperplane LSH. We
further describe how to use hypercube LSH in practice, and we consider an example application in the
area of lattice algorithms.
Keywords: (approximate) near neighbors, locality-sensitive hashing, large deviations, dimensionality reduction, lattice algorithms

1

Introduction

Finding (approximate) near neighbors. A key computational problem in various research areas, including machine learning, pattern recognition, data compression, coding theory, and cryptanalysis [SDI05, Bis06, Dub10, DHS00, MO15, Laa15], is finding near neighbors: given a data set D ⊂ Rd
of cardinality n, design a data structure and preprocess D in a way that, when given a query vector
q ∈ Rd , one can efficiently find a near point to q in D. Due to the “curse of dimensionality” [IM98] this
problem is known to be hard to solve exactly (in the worst case) in high dimensions d, so a common
relaxation of this problem is the (c, r)-approximate near neighbor problem ((c, r)-ANN): given that the
nearest neighbor lies at distance at most r from q, design an algorithm that finds an element p ∈ D
at distance at most c · r from q.
Locality-sensitive hashing (LSH) and filtering (LSF). A prominent class of algorithms for finding
near neighbors in high dimensions is formed by locality-sensitive hashing (LSH) [IM98] and localitysensitive filtering (LSF) [BDGL16]. These solutions are based on partitioning the space into regions,
in a way that nearby vectors have a higher probability of ending up in the same hash region than
distant vectors. By carefully tuning (i) the number of hash regions per hash table, and (ii) the number
of randomized hash tables, one can then guarantee that with high probability (a) nearby vectors will
collide in at least one of the hash tables, and (b) distant vectors will not collide in any of the hash
tables. For LSH, a simple lookup in all of q’s hash buckets then provides a fast way of finding near
neighbors to q, while for LSF the lookups are slightly more involved. For various metrics, LSH and
LSF currently provide the best performance in high dimensions [AR15, BDGL16, ALRW17, Chr17].
Near neighbors on the sphere. In this work we will focus on the near neighbor problem under the
angular distance, where two vectors x, y are considered nearby iff their common angle θ is small [Cha02,
STS+ 13, SSLM14, AIL+ 15]. This equivalently corresponds to near neighbor searching for the `2 -norm,
where the entire data set is assumed to lie on a sphere. A special
on the sphere,
√ case of (c, r)-ANN
√
often considered in the literature, is the random case r = 1c 2 and c · r = 2, in part due to a
reduction from near neighbor under the Euclidean metric for general data sets to (c, r)-ANN on the
sphere with these parameters [AR15].

2

1.1

Thijs Laarhoven

Related work

Upper bounds. Perhaps the most well-known and widely used solution for ANN for the angular distance
is Charikar’s hyperplane LSH [Cha02], where a set of random hyperplanes is used to partition the
space into regions. Due to its low computational complexity and the simple form of the collision
probabilities (with no hidden order terms in d), this method is easy to instantiate in practice and
commonly achieves the best performance out of all LSH methods when d is not too large. For large
d, both spherical cap LSH [AINR14, AR15] and cross-polytope LSH [TT07, ER08, AIL+ 15, KW17]
are known to perform better than hyperplane LSH. Experiments from [TT07, TT09] showed that
using orthogonal hyperplanes, partitioning the space into Voronoi regions induced by the vertices
of a hypercube, also leads to superior results compared to hyperplane LSH; however, no theoretical
guarantees for the resulting hypercube LSH method were given, and it remained unclear whether the
improvement persists in high dimensions.
Lower bounds. For the case of random data sets, lower bounds have also been found, matching the
performance of spherical cap and cross-polytope LSH for large c [MNP07, OWZ11, AIL+ 15]. These
lower bounds are commonly in a model where it is assumed that collision probabilities are “not too
small”, and in particular not exponentially small in d. Therefore it is not clear whether one can further
improve upon cross-polytope LSH when the number of hash regions is exponentially large, which would
for instance be the case for hypercube LSH. Together with the experimental results from [TT07,TT09],
this naturally begs the question: how efficient is hypercube LSH? Is it better than hyperplane LSH
and/or cross-polytope LSH? And how does hypercube LSH compare to other methods in practice?
1.2

Contributions

Hypercube LSH. By carefully analyzing the collision probabilities for hypercube LSH using results
from large deviations theory, we show that hypercube LSH is indeed different from, and superior to
hyperplane LSH for large d. The following main theorem states the asymptotic form of the collision
probabilities when using hypercube LSH, which are also visualized in Figure 1 in comparison with
hyperplane LSH.
Theorem 1 (Collision probabilities for hypercube LSH). Let X, Y ∼ N (0, 1)d , let θ ∈ [0, π]
denote the angle between X and Y , and let p(θ) denote the probability that X and Y are mapped to
the same hypercube hash region. For θ ∈ (0, arccos π2 ) (respectively θ ∈ (arccos π2 , π3 )), let β0 ∈ (1, ∞)
(resp. β1 ∈ (1, ∞)) be the unique solution to:
p
p
 
 
(β0 − cos θ) β02 − 1
(β1 + cos θ) β12 − 1
−1
1
arccos
=
,
arccos
=
.
(1)
β0
β0 (β0 cos θ − 1)
β1
β1 (β1 cos θ + 1)
Then, as d tends to infinity, p(θ) satisfies:

d+o(d)

(β0 − cos θ)2


,



πβ0 (β0 cos θ − 1) sin θ





d+o(d)



(β1 + cos θ)2


,
πβ1 (β1 cos θ + 1) sin θ
p(θ) =






1 + cos θ d+o(d)



,


π
sin
θ





0,

if θ ∈ [0, arccos π2 ];
if θ ∈ [arccos π2 , π3 ];

(2)

if θ ∈ [ π3 , π2 );
if θ ∈ [ π2 , π].

Denoting the query complexity of LSH methods by nρ+o(1) , the parameter ρ for hypercube LSH is
up to log2 (π) ≈ 1.65 times smaller than for
√ hyperplane LSH. For large d, hypercube LSH is dominated
by cross-polytope LSH (unless c · r > 2), but as the convergence to the limit is rather slow, in
practice either method might be better, depending on the exact parameter setting. For the random

Hypercube LSH for approximate near neighbors

3

1
Hyperplane LSH

→ p(θ)1/d

Hypercube LSH

ν

3
π
1
π

0

0

arccos( π2 )

π
3

π
2

π

→θ
√
Fig. 1. Asymptotics of collision probabilities for hypercube LSH, compared to hyperplane LSH. Here ν = π/(2 π 2 − 4),
and the dashed vertical lines correspond to boundary points of the piecewise parts of Theorem 1. The blue line indicates
hyperplane LSH with d random hyperplanes.

1

→ρ

0.5
0.2
0.1
0.05
1

Hyperplane LSH
Hypercube LSH
Cross-polytope LSH

2

2

2

2

4

→c
Fig. 2. Asymptotics for√the LSH exponent ρ when using hyperplane LSH, hypercube LSH, and cross-polytope LSH, for
(c, r)-ANN with c · r = 2. The curve for hyperplane LSH is exact for arbitrary d, while for the other two curves, order
terms vanishing as d → ∞ have been omitted.

setting, Figure 2 shows limiting values for ρ for hyperplane, hypercube and cross-polytope LSH. We
again remark that these are asymptotics for d → ∞, and may not accurately reflect the performance
of these methods for moderate d. We further briefly discuss how the hashing for hypercube LSH can
be made efficient.
Partial hypercube LSH. As the number of hash regions of a full-dimensional hypercube is often prohibitively large, we also consider partial hypercube LSH, where a d0 -dimensional hypercube is used
to partition a data set in dimension d. Building upon a result of Jiang [Jia06], we characterize when
hypercube and hyperplane LSH are asymptotically equivalent in terms of the relation between d0
and d, and we empirically illustrate the convergence towards either hyperplane or hypercube LSH for
larger d0 . An important open problem remains to identify how large the ratio d0 /d must be for the
asymptotics of partial hypercube LSH to be equivalent to those of full-dimensional hypercube LSH.
Application to lattice sieving. Finally, we consider a specific use case of different LSH methods, in
the context of lattice cryptanalysis. We show that the heuristic complexity of lattice sieving with

4

Thijs Laarhoven

hypercube LSH is expected to be slightly better than when using hyperplane LSH, and we discuss
how experiments have previously indicated that in this application, hypercube LSH is superior to
other dimensions up to dimensions d ≈ 80.

2

Preliminaries

Notation. We denote probabilities with P(·) and expectations with E(·). Capital letters commonly
denote random variables, and boldface letters denote vectors. We informally write P(X = x) for
continuous X to denote the density of X at x. For probability distributions D, we write X ∼ D to
denote that X is distributed according to D. For sets S, with abuse of notation we further write X ∼ S
to denote X is drawn uniformly at random from S. We write N (µ, σ 2 ) for the normal distribution
with mean µ and variance σ 2 , and H(µ, σ 2 ) for the distribution of |X| when X ∼ N (µ, σ 2 ). For µ = 0
the latter corresponds to the half-normal distribution. We write X ∼ Dd to denote a d-dimensional
qP
2
vector where each entry is independently distributed according to D. In what follows, kxk =
i xi
P
denotes the Euclidean norm, and hx, yi = i xi yi denotes the standard inner product. We denote the
angle between two vectors by φ(x, y) = arccoshx/kxk, y/kyki.
Lemma 1 (Distribution of angles between random vectors [BDGL16, Lemma 2]). Let
X, Y ∼ N (0, 1)d be two independent standard normal vectors. Then P(φ(X, Y ) = θ) = (sin θ)d+o(d) .
Locality-sensitive hashing. Locality-sensitive hash functions [IM98] are functions h mapping a ddimensional vector x to a low-dimensional sketch h(x), such that vectors which are nearby in Rd
are more likely to be mapped to the same sketch than distant vectors. For the angular distance1
φ(x, y), we quantify a set of hash functions H as follows (see [IM98]):
Definition 1. A hash family H is called (θ1 , θ2 , p1 , p2 )-sensitive if for x, y ∈ Rd we have:
– If φ(x, y) ≤ θ1 then Ph∼H (h(x) = h(y)) ≥ p1 ;
– If φ(x, y) ≥ θ2 then Ph∼H (h(x) = h(y)) ≤ p2 .
The existence of locality-sensitive hash families implies the existence of fast algorithms for (approximate) near neighbors, as the following lemma describes2 . For more details on the general principles
of LSH, we refer the reader to e.g. [IM98, And09].
Lemma 2 (Locality-sensitive hashing [IM98]). Suppose there exists a (θ1 , θ2 , p1 , p2 )-sensitive
log(p1 )
family H. Let ρ = log(p
. Then w.h.p. we can either find an element p ∈ L at angle at most θ2 from
2)
q, or conclude that no elements p ∈ L at angle at most θ1 from q exist, in time nρ+o(1) with space
and preprocessing costs n1+ρ+o(1) .
Hyperplane LSH. For the angular distance, Charikar [Cha02] introduced the hash family H = {ha :
a ∼ D} where D is any spherically symmetric distribution on Rd , and ha satisfies:
(
+1, if ha, xi ≥ 0;
ha (x) =
(3)
−1, if ha, xi < 0.
The vector a can be interpreted as the normal vector of a random hyperplane, and the hash value
depends on which side of the hyperplane x lies on. For this hash function, the probability of a collision
is directly proportional to the angle between x and y:

φ(x, y)
Ph∼H h(x) = h(y) = 1 −
.
π
For any two angles θ1 < θ2 , the above family H is (θ1 , θ2 , 1 −
1
2

θ1
π ,1

−

(4)
θ2
π )-sensitive.

Formally speaking, the angular distance is only a similarity measure, and not a metric.
Various conditions and order terms (which are commonly no(1) ) are omitted here for brevity.

Hypercube LSH for approximate near neighbors

5

Large deviations theory. Let P
{Z d }d∈N ⊂ Rk be a sequence of random vectors corresponding to an
empirical mean, i.e. Z d = d1 di=1 U i with U i i.i.d. We define the logarithmic moment generating
function Λ of Z d as:
Λ(λ) = ln EU 1 [exphλ, U 1 i] .

(5)

Define DΛ = {λ ∈ Rk : Λ(λ) < ∞}. The Fenchel-Legendre transform of Λ is defined as:
Λ∗ (z) = sup {hλ, zi − Λ(λ)} .

(6)

λ∈Rk

The following result describes that under certain conditions on {Z 0d }, the asymptotics of the probability
measure on a set F are related to the function Λ∗ .
Lemma 3 (Gärtner-Ellis theorem [DZ10, Theorem 2.3.6 and Corollary 6.1.6]). Let 0 be
contained in the interior of DΛ , and let Z d be an empirical mean. Then for arbitrary sets F ,
lim 1
d→∞ d

ln P(z ∈ F ) = − inf Λ∗ (z).
z∈F

(7)

The latter statement can be read as P(z ∈ F ) = exp(−d inf z∈F Λ∗ (z) + o(d)), and thus tells us
exactly how P(z ∈ F ) scales as d tends to infinity, up to order terms.

3

Hypercube LSH

In this section, we will analyze full-dimensional hypercube hashing, with hash family H = {hA : A ∈
SO(d)} where SO(d) ⊂ Rd×d denotes the rotation group, and hA satisfies:
(
+1, if xi ≥ 0;
hA (x) = (h1 (Ax), . . . , hd (Ax)),
hi (x) =
(8)
−1, if xi < 0.
In other words, a hypercube hash function first applies a uniformly random rotation, and then maps
the resulting vector to the orthant it lies in. This equivalently corresponds to a concatenation of d
hyperplane hash functions, where all hyperplanes are orthogonal. Collision probabilities for prescribed
angles θ between x and y are denoted by:
p(θ) = P(hA (x) = hA (y) | φ(x, y) = θ).

(9)

Above, the randomness is over hA ∼ H, with x and y arbitrary vectors at angle θ (e.g. x = e1 and
y = e1 cos θ + e2 sin θ). Alternatively, the random rotation A inside hA may be omitted, and the
probability can be computed over X, Y drawn uniformly at random from a spherically symmetric
distribution, conditioned on their common angle being θ.
3.1

Outline of the proof of Theorem 1

Although Theorem 1 is a key result, due to space restrictions we have decided to defer the full proof
(approximately 5.5 pages) to the appendix. The approach of the proof can be summarized by the
following four steps:
– Rewrite the collision probabilities in terms of (normalized) half-normal vectors X, Y ;
– Introduce dummy variables x, y for the norms of these half-normal vectors, so that the probability
can be rewritten in terms of unnormalized half-normal vectors;
– Apply the
Gärtner-Ellis
theorem (Lemma 3) to the three-dimensional vector given by
P 2 P
1 P
Z = d ( i Xi Yi , i Xi , i Yi2 ) to compute the resulting probabilities for arbitrary x, y;
– Maximize the resulting expressions over x, y > 0 to get the final result.
The majority of the technical part of the proof lies in computing Λ∗ (z), which involves a somewhat
tedious optimization of a multivariate function through a case-by-case analysis.

6

Thijs Laarhoven

A note on Gaussian approximations. From the (above outline of the) proof, and the observation
that the final optimization over x, y yields x = y = 1 as the optimum, one might wonder whether a
simpler analysis might be possible by assuming (half-)normal vectors are already normalized. Such a
computation however would only lead to an approximate solution, which is perhaps easiest to see by
computing collision probabilities for θ = 0. In the exact computation, where vectors are normalized,
hX, Y i = 1 implies X = Y . If however we do not take into account the norms of X and Y , and
do not condition on the norms being equal to 1, then hX, Y i = 1 could also mean that X, Y are
slightly longer than 1 and have a small, non-zero angle. In fact, such a computation would indeed
yield p(θ)1/d 6→ 0 as θ → 0.

3.2

Consequences of Theorem 1

From Theorem 1, we can draw
√ several conclusions. Substituting values for θ, we can find asymptotics
π 1/d
for p(θ), such as p( 3 )
= π3 + o(1) and p( π2 )1/d = π1 + o(1). We observe that the limiting function
of Theorem 1 (without the order terms) is continuous everywhere except at θ = π2 . To understand
the boundary θ = arccos π2 of the piece-wise limit function, note that two (normalized) half-normal
vectors X, Y have expected inner product EhX, Y i = π2 .
LSH exponents ρ for random settings. Using Theorem 1, we can explicitly compute LSH exponents
√
ρ for given angles θ1 and θ2 for large d. As an example, consider the random setting3 with c = 2,
corresponding to θ2 = π2 and θ1 = π3 . Substituting the collision probabilities from Theorem 1, we get
ρ → 1 − 21 logπ (3) ≈ 0.520 as d → ∞. To compare, if we had used random hyperplanes, we would have
gotten a limiting value ρ → log2 ( 32 ) ≈ 0.585. For the random case, Figure 2 compares limiting values ρ
using random and orthogonal hyperplanes, and using the asymptotically superior cross-polytope LSH.
Scaling at θ → 0 and asymptotics of ρ for large c. For θ close to 0, by Theorem 1 we are in the regime
defined by β0 . For cos θ = 1−ε with ε > 0 small, observe that β0 ≈ 1 satisfies
β0 > 1/ cos θ. Computing
√
2 2 3/2
a Taylor expansion around ε = 0, we eventually find β0 = 1 + ε + π ε + O(ε2 ). Substituting this
value β0 into p(θ) with cos θ = 1 − ε, we find:
√

p(θ) =

2√
1−
ε + O(ε)
π

!d+o(d)
.

(10)

To compare this with hyperplane LSH, recall that the collision
√ probability for d random hyperplanes
θ d
is equal to (1 − π ) . Since cos θ = 1 − ε translates to θ = 2ε(1 + O(ε)), the collision probabilities
√ √
for hyperplane hashing in this regime are also (1 − π2 ε + O(ε))d . In other words, for angles θ → 0,
the collision probabilities for hyperplane hashing and hypercube hashing are similar. This can also be
observed in Figure 1. Based on this result, we further deduce that in random settings with large c, for
hypercube LSH we have:

ln 1 −
ρ→

√

2
πc

+O

ln(1/π)

1
c2



√

2
=
+O
πc ln π



1
c2



0.393
≈
+O
c



1
c2


.

(11)

For hyperplane LSH, the numerator is the same, while the denominator is ln( 12 ) instead of ln( π1 ),
leading to values ρ which are a factor log2 π + o(1) ≈ 1.652 + o(1) larger. Both methods are inferior
to cross-polytope LSH for large d, as there ρ = O(1/c2 ) for large c [AIL+ 15].
3

√
√
Here we assume that c · r → ( 2)− ,√i.e. c · r approaches 2 from below. Alternatively, one might interpret this as
that if distant
points lie at distance 2 ± o(1), then we might expect approximately
half of them to lie at distance
√
√
less than 2, with query complexity O(n/2)ρ+o(1) = nρ+o(1) . If however c · r ≥ 2 then clearly ρ = 0, regardless of d
and c.

Hypercube LSH for approximate near neighbors

7

1◆
●
■
▲
▼
●
■
◆
▲
▼

→ p(θ)1/d

3
4

1
2

●
■

1
4

◆
▲
▼

0

●
▲
◆
■
▼

●
◆
▼ ●
■
▲
▲
■
◆
▼ ◆
●
■
▲ ●
■
▲ ●
◆
■
◆
▲ ◆
●
■
▲ ◆
●
▲ ●
■
▲ ●
■
◆
■
▲●
◆
Full hypercube LSH
■
◆
▲◆
●
■●
2-dimensional hypercube
■●
◆
■●
4-dimensional hypercube
◆
■●
◆
8-dimensional hypercube
■●
◆
■
16-dimensional hypercube
◆
●
■
32-dimensional hypercube
π
8

3π
8

π
4

π
2

→θ
Fig. 3. Empirical collision probabilities for hypercube LSH for small d. The green curve denotes the exact collision
probabilities for d = 2 from Proposition 1.

3.3

Convergence to the limit

To get an idea how hypercube LSH compares to other methods when d is not too large, we start by
giving explicit collision probabilities for the first non-trivial case, namely d = 2.
Proposition 1 (Square LSH). For d = 2, p(θ) = 1 −

2θ
π

for θ ≤

π
2

and p(θ) = 0 otherwise.

Proof. In two dimensions, two randomly rotated vectors X, Y at angle θ can be modeled as X =
(cos ψ, sin ψ) and Y = (cos(ψ + θ), sin(ψ + θ)) for ψ ∼ [0, 2π). The conditions X, Y > 0 are then
over the
equivalent to ψ ∈ (0, π2 ) ∩ (−θ, π2 − θ), which for θ < π2 occurs with probability π/2−θ
2π
randomness of ψ. As a collision can occur in any of the four quadrants, we finally multiply this
probability by 4 to obtain the stated result.
Figure 3 depicts p(θ)1/2 in green, along with hyperplane LSH (blue) and the asymptotics for
hypercube LSH (red). For larger d, computing p(θ) exactly becomes more complicated, and so instead
we performed experiments to empirically obtain estimates for p(θ) as d increases. These estimates are
also shown in Figure 3, and are based on 105 trials for each θ and d. Observe that as θ → π2 and/or d
grows larger, p(θ) decreases and the empirical estimates become less reliable. Points are omitted for
cases where no successes occurred.
Based on these estimates and our intuition, we conjecture that (1) for θ ≈ 0, the scaling of p(θ)1/d
is similar for all d, and similar to the asymptotic behavior of Theorem 1; (2) the normalized collision
probabilities for θ ≈ π2 approach their limiting value from below; and (3) p(θ) is likely to be continuous
for arbitrary d, implying that for θ → π2 , the collision probabilities tend to 0 for each d. These together
suggest that values for ρ are actually smaller when d is small than when d is large, and the asymptotic
estimate from Figure 2 might be pessimistic in practice. For the random setting, this would suggest
that ρ ≈ 0 regardless of c, as p(θ) → 0 as θ → π2 for arbitrary d.
Comparison with hyperplane/cross-polytope LSH. Finally, [TT07, Figures 1 and 2] previously illustrated that among several LSH methods, the smallest values ρ (for their parameter sets) are obtained
with hypercube LSH with d = 16, achieving smaller values ρ than e.g. cross-polytope LSH with
d = 256. An explanation for this can be found in:
– The (conjectured) convergence of ρ to its limit from below, for hypercube LSH;

8

Thijs Laarhoven

– The slow convergence of ρ to its limit (from above) for cross-polytope LSH4 .
This suggests that the actual values ρ for moderate dimensions d may well be smaller for hypercube
LSH (and hyperplane LSH) than for cross-polytope LSH. Based on the limiting cases d = 2 and
d → ∞, we further conjecture that compared to hyperplane LSH, hypercube LSH achieves smaller
values ρ for arbitrary d.
3.4

Fast hashing in practice

To further assess the practicality of hypercube LSH, recall that hashing is done as follows:
– Apply a uniformly random rotation A to x;
– Look at the signs of (Ax)i .
Theoretically, a uniformly random rotation will be rather expensive to compute, with A being a real,
dense matrix. As previously discussed in e.g. [Ach01], it may suffice to only consider a sparse subset
of all rotation matrices with a large enough amount of randomness, and as described in [AIL+ 15,
KW17] pseudo-random rotations may also be help speed up the computations in practice. As described
in [KW17], this can even be made provable, to obtain a reduced O(d log d) computational complexity
for applying a random rotation.
Finally, to compare this with cross-polytope LSH, note that cross-polytope LSH in dimension d
partitions the space in 2d regions, as opposed to 2d for hypercube hashing. To obtain a similar finegrained partition of the space with cross-polytopes, one would have to concatenate Θ(d/ log d) random
cross-polytope hashes, which corresponds to computing Θ(d/ log d) (pseudo-)random rotations, compared to only one rotation for hypercube LSH. We therefore expect hashing to be up to a factor
Θ(d/ log d) less costly.

4

Partial hypercube LSH

Since a high-dimensional hypercube partitions the space in a large number of regions, for various
applications one may only want to use hypercubes in a lower dimension d0 < d. In those cases, one would
first apply a random rotation to the data set, and then compute the hash based on the signs of the first
d0 coordinates of the rotated data set. This corresponds to the hash family H = {hA,d0 : A ∈ SO(d)},
with hA,d0 satisfying:
hA,d0 (x) = (h1 (Ax), . . . , hd0 (Ax)),

(
+1, if xi ≥ 0;
hi (x) =
−1, if xi < 0.

(12)

When “projecting” down onto the first d0 coordinates, observe that distances and angles are distorted:
the angle between the vectors formed by the first d0 coordinates of x and y may not be the same as
φ(x, y). The amount of distortion depends on the relation between d0 and d. Below, we will investigate
how the collision probabilities pd0 ,d (θ) for partial hypercube LSH scale with d0 and d, where pd0 ,d (θ) =
P(h(x) = h(y) | φ(x, y) = θ).
4.1

Convergence to hyperplane LSH

First, observe that for d0 = 1, partial hypercube LSH is equal to hyperplane LSH, i.e. p1,d (θ) = 1 − πθ .
For 1 < d0  d, we first observe that both (partial) hypercube LSH and hyperplane LSH can be
modeled by a projection onto d0 dimensions:
0

– Hyperplane LSH: x 7→ Ax with A ∼ N (0, 1)d ×d ;
0
– Hypercube LSH: x 7→ (A∗ )x with A ∼ N (0, 1)d ×d .
4

[AIL+ 15, Theorem 1] shows that the leading term in the asymptotics for ρ scales as Θ(ln d), with a first order term
scaling as O(ln ln d), i.e. a relative order term of the order O(ln ln d/ ln d).

Hypercube LSH for approximate near neighbors

9

Here A∗ denotes the matrix obtained from A after applying Gram-Schmidt orthogonalization to the
rows of A. In both cases, hashing is done after the projection by looking at the signs of the projected
vector. Therefore, the only difference lies in the projection, and one could ask: for which d0 , as a
function of d, are these projections equivalent? When is a set of random hyperplanes already (almost)
orthogonal?
This question was answered in [Jia06]: if d0 = o(d/ log d), then maxi,j |Ai,j −A∗i,j | → 0 in probability
as d → ∞ (implying A∗ = (1 + o(1))A), while for d0 = Ω(d/ log d) this maximum does not converge
to 0 in probability. In other words, for large d a set of d0 random hyperplanes in d dimensions is
(approximately) orthogonal iff d0 = o(d/ log d).
Proposition 2 (Convergence to hyperplane LSH). Let pd0 ,d (θ) denote the collision probabilities
0
for partial hypercube LSH, and let d0 = o(d/ log d). Then pd0 ,d (θ)1/d → 1 − πθ .
As d0 = Ω(d/ log d) random vectors in d dimensions are asymptotically not orthogonal, in that case
one might expect either convergence to full-dimensional hypercube LSH, or to something in between
hyperplane and hypercube LSH.
4.2

Convergence to hypercube LSH

To characterize when partial hypercube LSH is equivalent to full hypercube LSH, we first observe that
if d0 is large compared to ln n, then convergence to the hypercube LSH asymptotics follows from the
Johnson-Lindenstrauss lemma.
Proposition 3 (Sparse data sets). Let d0 = ω(ln n). Then the same asymptotics for the collision
probabilities as those of full-dimensional hypercube LSH apply.
Proof. Let θ ∈ (0, π2 ). By the Johnson-Lindenstrauss lemma [JL84], we can construct a projection
x 7→ Ax from d onto d0 dimensions, preserving all pairwise distances up to a factor 1 ± ε for ε =
Θ((ln n)/d0 ) = o(1). For fixed θ ∈ (0, π2 ), this implies the angle φ between Ax and Ay will be in the
interval θ ± o(1), and so the collision probability lies in the interval p(θ ± o(1)). For large d, this means
that the asymptotics of p(θ) are the same.
To analyze collision probabilities for partial hypercube LSH when neither of the previous two
propositions applies, note that through a series of transformations similar to those for full-dimensional
hypercube LSH, it is possible to eventually end up with the following probability to compute, where
d1 = d0 and d2 = d − d0 :
max P

x,y,u,v,φ

d1
1 X
Xi Yi = xy cos φ,
d1
i=1

d2
1 X
Ui Vi = uvf (φ, θ),
d2
i=1

d1
1 X
Xi2 = x2 ,
d1

1
d2

i=1
d2
X
i=1

Ui2 = u2 ,

d1
1 X
Yi2 = y 2 ,
d1
i=1
!
d2
1 X
Vi2 = v 2 .
d2

(13)

(14)

i=1

Here f is some function of φ and θ. The approach is comparable to how we ended up with a similar
probability to compute in the proof of Theorem 1, except that we split the summation indices I = [d]
into two sets I1 = {1, . . . , d0 } of size d1 and I2 = {d0 +1, . . . , d} of size d2 . We then substitute Ui = Xd0 +i
and Vi = Yd0 +i , and add dummy variables x, y, u, v for the norms of the four partial vectors, and a
dummy angle φ for the angle between the d1 -dimensional vectors, given the angle θ between the
d-dimensional vectors.
Although the vector Z formed by the six random variables in (14) is not an empirical mean
over a fixed number d of random vectors (the first three are over d1 terms, the last three over d2
terms), one may expect a similar large deviations result such as Lemma 3 to apply here. In that
case, the function Λ∗ (z) = Λ∗ (z1 , . . . , z6 ) would be a function of six variables, which we would like to
evaluate at (xy cos φ, x2 , y 2 , uvf (φ, θ), u2 , v 2 ). The function Λ∗ itself involves an optimization (finding
a supremum) over another six variables λ = (λ1 , . . . , λ6 ), so to compute collision probabilities for

10

Thijs Laarhoven

1◆
●
■
▲
▼
■
◆
●
▲
▼

→ p(θ)1/d

3
4

1
2

●
■
◆
▲
▼

●
■
◆
▲ ◆
▼
●
■
▲ ◆
▼
●
■ ●
▲
■
▼ ◆
▲ ◆
●
■
▼
●
▲
■
▼ ◆
●
▲
■
●
▼ ◆
■
▲
●
■
▲
▼ ◆
●
▼ ◆
■
▲
●
■◆
▼
▲◆
●●
▼◆
▲
■◆
▼ ■
■◆
▲◆
●
▲●
■◆
■◆
●
▼
■
▼▲
▲●
▼▼
Partial hypercube LSH (d = 50)
▲
▲
▼
▼
● 2-dimensional hypercube
■ 4-dimensional hypercube
◆ 8-dimensional hypercube

1
4

▲ 16-dimensional hypercube
▼ 32-dimensional hypercube

0
0

π
8

π
4

3π
8

π
2

→θ
0

Fig. 4. Experimental values of pd0 ,50 (θ)1/d , for different values d0 , compared with the asymptotics for hypercube LSH
(red) and hyperplane LSH (blue).

given d, d0 , θ exactly, using large deviations theory, one would have to compute an expression of the
following form:
)
(
min

x,y,u,v,φ

sup

Fd,d0 ,θ (x, y, u, v, φ, λ1 , λ2 , λ3 , λ4 , λ5 , λ6 ) .

(15)

λ1 ,λ2 ,λ3 ,λ4 ,λ5 ,λ6

As this is a very complex task, and the optimization will depend heavily on the parameters d, d0 , θ
defined by the problem setting, we leave this optimization as an open problem. We only mention that
intuitively, from the limiting cases of small and large d0 we expect that depending on how d0 scales
with d (or n), we obtain a curve somewhere in between the two curves depicted in Figure 1.
4.3

Empirical collision probabilities

To get an idea of how pd0 ,d (θ) scales with d0 in practice, we empirically computed several values for
fixed d = 50. For fixed θ we then applied a least-squares fit of the form ec1 d+c2 to the resulting
data, and plotted ec1 in Figure 4. These data points are again based on at least 105 experiments for
each d0 and θ. We expect that as d0 increases, the collision probabilities slowly move from hyperplane
hashing towards hypercube hashing, this can also be seen in the graph – for d0 = 2, the least-squares
fit is almost equal to the curve for hyperplane LSH, while as d0 increases the curve slowly moves
down towards the asymptotics for full hypercube LSH. Again, we stress that as d0 becomes larger, the
empirical estimates become less reliable, and so we did not consider even larger values for d0 .
Compared to full hypercube LSH and Figure 3, we observe that we now approach the limit from
above (although the fitted collision probabilities never seem to be smaller than those of hyperplane
LSH), and therefore the values ρ for partial hypercube LSH are likely to lie in between those of
hyperplane and (the asymptotics of) hypercube LSH.

5

Application: Lattice sieving for the shortest vector problem

We finally consider an explicit application for hypercube LSH, namely lattice sieving algorithms
for the
P
d
shortest vector problem. Given a basis B = {b1 , . . . , bd } ⊂ R of a lattice L(B) = { i λi bi : λi ∈ Z},

Hypercube LSH for approximate near neighbors

11

the shortest vector problem (SVP) asks to find a shortest non-zero vector in this lattice. Various
different methods for solving SVP in high dimensions are known, and currently the algorithm with the
best heuristic time complexity in high dimensions is based on lattice sieving, combined with nearest
neighbor searching [BDGL16].
In short, lattice sieving works by generating a long list L of pairwise reduced lattice vectors, where
x, y are reduced iff kx − yk ≥ min{kxk, kyk}. The previous condition is equivalent to φ(x, y) ≤ π3 ,
and so the length of L can be bounded by the kissing constant in dimension d, which is conjectured to
scale as (4/3)d/2+o(d) . Therefore, if we have a list of size n = (4/3)d/2+o(d) , any newly sampled lattice
vector can be reduced against the list many times to obtain a very short lattice vector. The time
complexity of this method is dominated by doing poly(d) · n reductions (searches for nearby vectors)
with a list of size n. A linear search trivially leads to a heuristic complexity of n2+o(1) = (4/3)d+o(d)
(with space n1+o(1) ), while nearest neighbor techniques can reduce the time complexity to n1+ρ+o(1)
for ρ < 1 (increasing the space to n1+ρ+o(1) ). For more details, see e.g. [NV08, Laa15, BDGL16].
Based on the collision probabilities for hypercube LSH, and assuming the asymptotics for partial
hypercube LSH (with d0 = O(d)) are similar to those of full-dimensional hypercube LSH, we obtain
the following result. An outline of the proof is given in the appendix.
Proposition 4 (Complexity of lattice sieving with hypercube LSH). Suppose the asymptotics
for full hypercube LSH also hold for partial hypercube LSH with d0 ≈ 0.1335d. Then lattice sieving with
hypercube LSH heuristically solves SVP in time and space 20.3222d+o(d) .
As expected, the conjectured asymptotic performance of (sieving with) hypercube LSH lies in
between those of hyperplane LSH and cross-polytope LSH.
–
–
–
–
–
–

Linear search [NV08]:
20.4150d+o(d) .
Hyperplane LSH [Laa15]:
20.3366d+o(d) .
Hypercube LSH:
20.3222d+o(d) .
Spherical cap LSH [LdW15]: 20.2972d+o(d) .
Cross-polytope LSH [BL16]: 20.2972d+o(d) .
Spherical LSF [BDGL16]:
20.2925d+o(d) .

In practice however, the picture is almost entirely reversed [SG15]. The lattice sieving method used to
solve SVP in the highest dimension to date (d = 116) used a very optimized linear search [Kle14]. The
furthest that any nearest neighbor-based sieve has been able to go to date is d = 107, using hypercube
LSH [MLB15, MB16]5 . Experiments further indicated that spherical LSF only becomes competitive
with hypercube LSH as d & 80 [BDGL16, MLB17], while sieving with cross-polytope LSH turned
out to be rather slow compared to other methods [BL16, Mar16]. Although it remains unclear which
nearest neighbor method is the “most practical” in the application of lattice sieving, hypercube LSH
is one of the main contenders.
Acknowledgments. The author is indebted to Ofer Zeitouni for his suggestion to use results from large
deviations theory, and for his many helpful comments regarding this application. The author further
thanks Brendan McKay and Carlo Beenakker for their comments. The author is supported by the
SNSF ERC Transfer Grant CRETP2-166734 FELICITY.

References
Ach01.
AIL+ 15.

Dimitris Achlioptas. Database-friendly random projections. In PODS, pages 274–281, 2001.
Alexandr Andoni, Piotr Indyk, Thijs Laarhoven, Ilya Razenshteyn, and Ludwig Schmidt. Practical and
optimal LSH for angular distance. In NIPS, pages 1225–1233, 2015.
AINR14. Alexandr Andoni, Piotr Indyk, Huy Lê Nguyên, and Ilya Razenshteyn. Beyond locality-sensitive hashing. In
SODA, pages 1018–1028, 2014.
ALRW17. Alexandr Andoni, Thijs Laarhoven, Ilya Razenshteyn, and Erik Waingarten. Optimal hashing-based timespace trade-offs for approximate near neighbors. In SODA, pages 47–66, 2017.
And09.
Alexandr Andoni. Nearest Neighbor Search: the Old, the New, and the Impossible. PhD thesis, Massachusetts
Institute of Technology, 2009.
5

Although phrased as hyperplane LSH, the implementations from [Laa15, MLB15, MB16] are using hypercube LSH.

12

Thijs Laarhoven

AR15.

Alexandr Andoni and Ilya Razenshteyn. Optimal data-dependent hashing for approximate near neighbors.
In STOC, pages 793–801, 2015.
AS72.
Milton Abramowitz and Irene A. Stegun. Handbook of Mathematical Formulas. Dover Publications, 1972.
BDGL16. Anja Becker, Léo Ducas, Nicolas Gama, and Thijs Laarhoven. New directions in nearest neighbor searching
with applications to lattice sieving. In SODA, pages 10–24, 2016.
Bis06.
Christopher M. Bishop. Pattern Recognition and Machine Learning (Information Science and Statistics).
Springer-Verlag, 2006.
BL16.
Anja Becker and Thijs Laarhoven.
Efficient (ideal) lattice sieving using cross-polytope LSH.
In
AFRICACRYPT, pages 3–23, 2016.
Cha02.
Moses S. Charikar. Similarity estimation techniques from rounding algorithms. In STOC, pages 380–388,
2002.
Chr17.
Tobias Christiani. A framework for similarity search with space-time tradeoffs using locality-sensitive filtering.
In SODA, pages 31–46, 2017.
DHS00.
Richard O. Duda, Peter E. Hart, and David G. Stork. Pattern Classification (2nd Edition). Wiley, 2000.
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.
DZ10.
Amir Dembo and Ofer Zeitouni. Large deviations techniques and applications (2nd edition). Springer, 2010.
ER08.
Kave Eshghi and Shyamsundar Rajaram. Locality sensitive hash functions based on concomitant rank order
statistics. In KDD, pages 221–229, 2008.
IM98.
Piotr Indyk and Rajeev Motwani. Approximate nearest neighbors: Towards removing the curse of dimensionality. In STOC, pages 604–613, 1998.
Jia06.
Tiefeng Jiang. How many entries of a typical orthogonal matrix can be approximated by independent normals?
The Annals of Probability, 34(4):1497–1529, 2006.
JL84.
William B. Johnson and Joram Lindenstrauss. Extensions of Lipschitz mappings into a Hilbert space. Contemporary Mathematics, 26(1):189–206, 1984.
Kle14.
Thorsten Kleinjung. Private communication, 2014.
KW17.
Christopher Kennedy and Rachel Ward. Fast cross-polytope locality-sensitive hashing. In ITCS, 2017.
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 localitysensitive hashing. In LATINCRYPT, pages 101–118, 2015.
Mar16.
Artur Mariano. Private communication., 2016.
MB16.
Artur Mariano and Christian Bischof. Enhancing the scalability and memory usage of HashSieve on multi-core
CPUs. In PDP, pages 545–552, 2016.
MLB15. Artur Mariano, Thijs Laarhoven, and Christian Bischof. Parallel (probable) lock-free HashSieve: a practical
sieving algorithm for the SVP. In ICPP, pages 590–599, 2015.
MLB17. Artur Mariano, Thijs Laarhoven, and Christian Bischof. A parallel variant of LDSieve for the SVP on lattices.
PDP, 2017.
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.
NV08.
Phong Q. Nguyên and Thomas Vidick. Sieve algorithms for the shortest vector problem are practical. Journal
of Mathematical Cryptology, 2(2):181–207, 2008.
OWZ11. Ryan O’Donnell, Yi Wu, and Yuan Zhou. Optimal lower bounds for locality sensitive hashing (except when
q is tiny). In ICS, pages 276–283, 2011.
SDI05.
Gregory Shakhnarovich, Trevor Darrell, and Piotr Indyk. Nearest-Neighbor Methods in Learning and Vision:
Theory and Practice. MIT Press, 2005.
SG15.
Michael Schneider and Nicolas Gama. SVP challenge, 2015.
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.
TT09.
Kengo Terasawa and Yuzuru Tanaka. Approximate nearest neighbor search for a dataset of normalized
vectors. In IEICE Transactions on Information and Systems, volume 92, pages 1609–1619, 2009.

A

Proof of Theorem 1

Theorem 1 will be proved through a series of lemmas, each making partial progress towards a final
solution. Reading only the claims made in the lemmas may give the reader an idea how the proof is
built up. Before starting the proof, we begin with a useful lemma regarding integrals of (exponentials
of) quadratic forms.

Hypercube LSH for approximate near neighbors

13

Lemma 4 (Integrating an exponential of a quadratic form in the positive quadrant). Let
a, b, c ∈ R with a, c < 0 and D = b2 − 4ac < 0. Then:


b
Z ∞Z ∞
π + 2 arctan √−D
√
exp(ax2 + bxy + cy 2 ) dx dy =
.
(16)
2 −D
0
0
Proof. The proof below is based on substituting y = xs (and dy = x ds) before computing the integral
over x. An integral over 1/(a + bs + cs2 ) then remains, which leads to the arctangent solution in case
b2 < 4ac.
Z ∞Z ∞
I=
exp(ax2 + bxy + cy 2 ) dx dy
(17)
y=0 0

Z ∞ Z ∞

x exp (a + bs + cs2 )x2 dx ds
(18)
=
s=0
0
 #∞
Z ∞"
exp (a + bs + cs2 )x2
=
ds
(19)
2(a + bs + cs2 )
0
x=0

Z ∞
1
=
0−
ds
(20)
2(a + bs + cs2 )
0
Z
1
−1 ∞
ds.
(21)
=
2 0 a + bs + cs2
The last equality used the assumptions a, c < 0 and b2 < 4ac so that a + bs + cs2 < 0 for all s > 0.
We then solve the last remaining integral (see e.g. [AS72, Equation (3.3.16)]) to obtain:


∞
−1
2
b + 2cs
√
√
I=
arctan
(22)
2
4ac − b2
4ac − b2 s=0



−1
b
= √
−π − 2 arctan √
.
(23)
2 4ac − b2
4ac − b2
Eliminating minus signs and substituting D = b2 − 4ac, we obtain the stated result.
Next, we begin by restating the collision probability between two vectors in terms of half-normal
vectors.
Lemma 5 (Towards three-dimensional large deviations). Let H denote the hypercube hash
family in d dimensions, and as before, let p be defined as:
p(θ) = Ph∼H (h(x) = h(y) | φ(x, y) = θ).
Let X̂, Ŷ ∼ H(0, 1)d and let the sequence {Z d }d∈N ⊂ R3 be defined as:
!
d
d
d
X
X
1 X
2
2
Zd =
X̂i Ŷi ,
X̂i ,
Ŷi .
d
i=1

i=1

(24)

(25)

i=1

Then:

p(θ) =

1
2 sin θ

d+o(d)

max P(Z d = (xy cos θ, x2 , y 2 )).

x,y>0

(26)

Proof. First, we write out the definition of the conditional probability in p, and use the fact that
each of the 2d hash regions (orthants) has the same probability mass. Here X, Y ∼ N (0, 1)d denote
random Gaussian vectors, and subscripts denoting what probabilities are computed over are omitted
when implicit.
p(θ) = Ph∼H (h(x) = h(y) | φ(x, y) = θ)
= 2d · PX,Y ∼N (0,1)d (X > 0, Y > 0 | φ(X, Y ) = θ)
=

2d

· P(X > 0, Y > 0, φ(X, Y ) = θ)
.
P(φ(X, Y ) = θ)

(27)
(28)
(29)

14

Thijs Laarhoven

By Lemma 1, the denominator is equal to (sin θ)d+o(d) . The numerator of (29) can further be rewritten
as a conditional probability on {X > 0, Y > 0}, multiplied with P(X > 0, Y > 0) = 2−2d . To
incorporate the conditionals X, Y > 0, we replace X, Y ∼ N (0, 1)d by half-normal vectors X̂, Ŷ ∼
H(0, 1)d , resulting in:
p(θ) =

PX̂,Ŷ ∼H(0,1)d (φ(X̂, Ŷ ) = θ)
(2 sin θ)d+o(d)

=

q(θ)
(2 sin θ)d+o(d)

.

(30)

To incorporate the normalization over the (half-normal)
vectors
X̂ and Ŷ , we introduce dummy
√
√
variables x, y corresponding to the norms of X̂/ d and Ŷ / d, and observe that as the probabilities
are exponential in d, the integrals will be dominated by the maximum value of the integrand in the
given range:
Z ∞Z ∞
P(hX̂, Ŷ i = x y d cos θ, kX̂k2 = x2 d, kŶ k2 = y 2 d) dx dy
(31)
q(θ) =
0
0


= 2o(d) max P hX̂, Ŷ i = x y d cos θ, kX̂k2 = x2 d, kŶ k2 = y 2 d .
(32)
x,y>0

Substituting Z d = d1 (hX̂, Ŷ i, kX̂k2 , kŶ k2 ), we obtain the claimed result.
Note that Z1 , Z2 , Z3 are pairwise but not jointly independent. To compute the density of Z d at
(xy cos θ, x2 , y 2 ) for d → ∞, we use the Gärtner-Ellis theorem stated in Lemma 3.
Lemma 6 (Applying the Gärtner-Ellis theorem to Z d ). Let {Z d }d∈N ⊂ R3 as in Lemma 5,
and let Λ and Λ∗ as in Section 2. Then 0 lies in the interior of DΛ , and therefore

P(Z d = (xy cos θ, x2 , y 2 )) = exp −Λ∗ (xy cos θ, x2 , y 2 )d + o(d) .
(33)
Essentially, all that remains now is computing Λ∗ at the appropriate point z. To continue, we first
compute the logarithmic moment generating function Λ = Λd of Z d :
Lemma 7 (Computing Λ). Let Z d as before, and let D = D(λ1 , λ2 , λ3 ) = λ21 − (1 − 2λ2 )(1 − 2λ3 ).
Then for λ ∈ DΛ = {λ ∈ R3 : λ2 , λ3 < 21 , D < 0} we have:



λ1
Λ(λ) = ln π + 2 arctan √
− ln π − 12 ln(−D).
(34)
−D
Proof. By the definition of the LMGF, we have:
h

i
Λ(λ) = ln EX̂1 ,Ŷ1 ∼H(0,1) exp λ1 X̂1 Ŷ1 + λ2 X̂12 + λ3 Ŷ12 .

(35)

We next compute the inner expectation over the random variables X̂1 , Ŷ1 , by writing out the double
integral over the product of the argument with the densities of X̂1 and Ŷ1 .


EX1 ,Y1 exp λ1 X1 Y1 + λ2 X12 + λ3 Y12
(36)
 2 Z ∞ r
 2
Z ∞r

2
x
2
y
=
exp −
dx
exp −
dy exp λ1 xy + λ2 x2 + λ3 y 2
(37)
π
2
π
2
0
0
Z ∞Z ∞

 
2
=
exp λ1 xy + λ2 − 12 x2 + λ3 − 21 y 2 dx dy .
(38)
π 0
0
Applying Lemma 4 with (a, b, c) = (λ2 − 21 , λ1 , λ3 − 12 ) yields the claimed expression for Λ, as well as
the bounds stated in DΛ which are necessary for the expectation to be finite.
We now continue with computing the Fenchel-Legendre transform of Λ, which involves a rather
complicated maximization (supremum) over λ ∈ R3 . The following lemma makes a first step towards
computing this supremum.

Hypercube LSH for approximate near neighbors

15

Lemma 8 (Computing Λ∗ (z) – General form). Let z ∈ R3 such that z2 , z3 > 0. Then the
Fenchel-Legendre transform Λ∗ of Λ at z satisfies

√
z2 z3
1
Λ∗ (z) = ln π + sup
+
+ λ1 z1 − |λ1 |β z2 z3 + ln(β 2 − 1) + ln |λ1 |
(39)
2
2
2
λ1 ,β
β>1

)

λ1
p
.
− ln π + 2 arctan
|λ1 | β 2 − 1




(40)

Proof. First, we recall the definition of Λ∗ and substitute the previous expression for Λ:
Λ∗ (z) = sup {hλ, zi − Λ(λ)}
λ∈R3




√
λ1
= ln π + sup hλ, zi + ln −D − ln π + 2 arctan √
.
−D
λ∈R3

(41)
(42)

Here as before D = λ21 − (1 − 2λ2 )(1 − 2λ3 ) < 0. Let the argument of the supremum above be denoted
by f (z, λ). We make a change of variables by setting t2 = 1 − 2λ2 > 0 and t3 = 1 − 2λ3 > 0, so that
D becomes D = λ21 − t2 t3 < 0:
z2 z3
t2 z2 t3 z3
+
+ λ 1 z1 −
−
2
2
2
2

2
1
+ 2 ln(t2 t3 − λ1 ) − ln π + 2 arctan √

f (z, λ1 , t2 , t3 ) =

(43)
λ1
t2 t3 −λ21



.

(44)

We continue by making a further change of variables u = t2 t3 > λ21 so that t2 = u/t3 . As a result the
dependence of f on t3 is only through the fourth and
p fifth terms above, from which one can
p easily
deduce that the supremum over t3 occurs at t3 = uz2 /z3 . This also implies that t2 = uz3 /z2 .
Substituting these values for t2 , t3 , we obtain:



√
f (z, λ1 , u) = z22 + z23 + λ1 z1 − uz2 z3 + 21 ln(u − λ21 ) − ln π + 2 arctan √ λ1 2 .
u−λ1

Finally, we use the substitution u = β 2 ·λ21 . From D < 0 it follows that u/λ21 = β > 1. This substitution
and some rewriting of f leads to the claimed result.
The previous simplifications were regardless of z1 , z2 , z3 , where the only assumption that was
made during the optimization of t3 was that z2 , z3 > 0. In our application, we want to compute Λ∗ at
z = (xy cos θ, x2 , y 2 ) for certain x, y > 0 and θ ∈ (0, π2 ). Substituting these values for z, the expression
from Lemma 5 becomes:

x2 y 2
1
∗
2 2
Λ (xy cos θ, x , y ) = ln π +
+
+ sup (λ1 cos θ − |λ1 |β)xy + ln(β 2 − 1)
(45)
2
2
2
λ1 ,β
β>1

)

λ1
p
+ ln |λ1 | − ln π + 2 arctan
.
|λ1 | β 2 − 1




(46)

The remaining optimization over λ1 , β now takes slightly different forms depending on whether λ1 < 0
or λ1 > 0. We will tackle these two cases separately, based on the identity:
n
o
Λ∗ (z) = max sup {hλ, zi − Λ(λ)} , sup {hλ, zi − Λ(λ)} = max{Λ∗+ (z), Λ∗− (z)}.
λ∈R3
λ1 >0

λ∈R3
λ1 <0

Lemma 9 (Computing Λ∗ (z) for positive λ1 ). Let z = (xy cos θ, x2 , y 2 ) with x, y > 0 and θ ∈
(0, π2 ). For θ ∈ (0, arccos π2 ), let β0 = β0 (θ) ∈ (1, ∞) be the unique solution to (1). Then the FenchelLegendre transform Λ∗ at z, restricted to λ1 > 0, satisfies
 

πβ0 (β0 cos θ − 1)


ln
, if θ ∈ (0, arccos π2 );
2
x2 y 2
2(β
−
cos
θ)
∗
0
Λ+ (z) =
+
− 1 − ln(xy) +
(47)

2
2

0,
if θ ∈ [arccos 2 , π ).
π

2

16

Thijs Laarhoven

Proof. Substituting λ1 > 0 into (46), we obtain:
Λ∗+ (xy cos θ, x2 , y 2 ) = ln π +

n
o
x2 y 2
+
+ sup g+ (λ1 , β) ,
2
2
λ1 >0

(48)

β>1

g+ (λ1 , β) = (cos θ − β)λ1 xy +

ln(β 2

− 1)

2




1
.
+ ln λ1 − ln π + 2 arctan p
β2 − 1

(49)

Differentiating w.r.t. λ1 gives (cos θ − β)xy + λ11 . Recall that β > 1 > cos θ. For λ1 → 0+ the derivative
is therefore positive, for λ1 → ∞ it is negative, and there is a global maximum at the only root
λ1 = 1/((β − cos θ)xy). In that case, the expression further simplifies and we can pull out more terms
that do not depend on β, to obtain:
n
o
x2 y 2
Λ∗+ (xy cos θ, x2 , y 2 ) = ln π +
+
− 1 − ln(xy) + sup g+ (β) ,
(50)
2
2
β>1


p
2
β −1

  = ln h+ (β).
g+ (β) = ln 
(51)
1
(β − cos θ) π + 2 arcsin β
p
Here we used the identity arctan(1/ β 2 − 1) = arcsin(1/β). Now, for β → 1+ we have h+ (β) → 0+ ,
while for β → ∞, we have


 
1
1
2
1
h+ (β) = +
cos θ −
+O
.
(52)
π πβ
π
β2
In other words, if cos θ ≤ π2 or θ ≥ arccos π2 , we have h+ (β) → ( π1 )− (the second order term is negative
for cos θ = π2 ), while for θ < arccos π2 we approach the same limit from above as h+ (β) → ( π1 )+ . For
θ < arccos π2 there is a non-trivial maximum at some value β = β0 ∈ (1, ∞), while for θ ≥ arccos π2 ,
we can see from the derivative h0+ (β) that h+ (β) is strictly increasing on (1, ∞), and the supremum
is attained at β → ∞. We therefore obtain two different results, depending on whether θ < arccos π2
or θ ≥ arccos π2 .
Case 1: arccos π2 ≤ θ < π2 . The supremum is attained in the limit of β → ∞, which leads to
h+ (β) → π1 and the stated expression for Λ∗+ (xy cos θ, x2 , y 2 ).
Case 2: 0 < θ < arccos π2 . In this case there is a non-trivial maximum at some value β = β0 , namely
there where the derivative h0+ (β0 ) = 0. After computing the derivative, eliminating the (positive)
denominator and rewriting, this condition is equivalent to (1).This allows us to rewrite g and Λ∗ in
terms of β0 , by substituting the given expression for arcsin β10 , which ultimately leads to the stated
formula for Λ∗+ .
Lemma 10 (Computing Λ∗ (z) for negative λ1 ). Let z = (xy cos θ, x2 , y 2 ) with x, y > 0 and
θ ∈ (0, π2 ). For θ ∈ (arccos π2 , π3 ), let β1 ∈ (1, ∞) be the unique solution to (1). Then the FenchelLegendre transform Λ∗ at z, restricted to λ1 < 0, satisfies


0,
if θ ∈ (0, arccos π2 ];







  πβ1 (β1 cos θ + 1) 
2
2
y
x
∗
, if θ ∈ (arccos π2 , π3 );
Λ− (z) =
+
− 1 − ln(xy) + ln
(53)
2
2(cos
θ
+
β
)
1

2
2

 




π


,
if θ ∈ [ π3 , π2 ).
ln
2(1 + cos θ)
Proof. We again start by substituting λ1 < 0 into (46):
n
o
x2 y 2
Λ∗− (xy cos θ, x2 , y 2 ) = ln π +
+
+ sup g− (λ1 , β) ,
2
2
λ1 <0
β>1

g− (λ1 , β) = (cos θ + β)λ1 xy +

ln(β 2
2

− 1)


 −1 
+ ln(−λ1 ) − ln π + 2 arctan p
.
β2 − 1

(54)

Hypercube LSH for approximate near neighbors

17

Differentiating w.r.t. λ1 gives (cos θ + β)xy + λ11 . For λ1 → −∞ this is positive, for λ1 → 0− this
is negative, and so the maximum is at λ1 = −1/((cos θ + β)xy). Substituting this value for λ1 , and
pulling out terms which do not depend on β yields:
Λ∗− (xy cos θ, x2 , y 2 ) = ln
g− (β) = ln

π 
2

+

n
o
x2 y 2
+
− 1 − ln(xy) + sup g− (β) ,
2
2
β>1
!

p
β2 − 1
(cos θ + β) arccos β1

(55)

= ln h− (β).

p
Above we used the identity π + 2 arctan(−1/ β 2 − 1) = 2 arccos β1 , where the factor 2 has been pulled
outside the supremum. Now, differentiating h− w.r.t. β results in:
p

2 − 1(β cos θ + 1) arccos 1 − β 2 − 1 (cos θ + β)
β
β
β
h0− (β) =
.
(56)
β (β 2 − 1) (cos θ + β)2 arccos β1
Clearly the denominator is positive, while for β → 1+ the limit is negative iff cos θ < 21 . For β → ∞
we further have h0− (β) → 0− for cos θ ≤ π2 and h0− (β) → 0+ for cos θ > π2 . We therefore analyze three
cases separately below.
Case 1: π3 ≤ θ < π2 . In this parameter range, h0− (β) is negative for all β > 1, and the supremum
1
∗
lies at β → 1+ with limiting value h− (β) → 1+cos
θ . This yields the given expression for Λ− .
Case 2: arccos π2 < θ < π3 . For θ in this range, h0− (β) is positive for β → 1+ and negative for
β → ∞, and changes sign exactly once, where it attains its maximum. After some rewriting, we find
that this is at the value β = β1 (θ) ∈ (1, ∞) satisfying the relation from (1). Substituting this expression
for arccos β11 into h− , we obtain the result for Λ∗− .
Case 3: 0 < θ ≤ arccos π2 . In this case h0− is positive for all β > 1, and the supremum lies at
β → ∞. For β → ∞ we have h− (β) → π2 (regardless of θ) and we therefore get the final claimed result.
Proof (Proof of Theorem 1). Combining the previous two results with Lemma 6 and Equation 47,
we obtain explicit asymptotics for P(Z d ≈ (xy cos θ, x2 , y 2 )). What remains is a maximization over
2
2
x, y > 0 of p, which translates to a minimization of Λ∗ . As x2 + y2 − 1 − ln(xy) attains its minimum
at x = y = 1 with value 0, we obtain Theorem 1.

B

Proof of Proposition 4

We will assume the reader is familiar with (the notation from) [Laa15]. Let t = 2ct d+o(d) denote the
number of hash tables, and n = (4/3)d/2+o(d) . Going through the proofs of [Laa15, Appendix A] and
replacing the explicit instantiation of the collision probabilities (1 − θ/π) by an arbitrary function
p(θ), we get that the optimal number of hash functions concatenated into one function for each hash
table, denoted k, satisfies
k=

ln t
ct d
√ .
=
0
− ln p(θ1 )
d log2 (π/ 3)

(57)

The latter equality follows when substituting θ1 = π/3 and substituting the collision probabilities for
partial hypercube LSH in some dimension d0 ≤ d. As we need k ≥ 1, the previous relation translates
ct √
d. As we expect the collision probabilities to be closer to those of
to a condition on d0 as d0 ≤ log (π/
3)
2

full-dimensional hypercube LSH when d0 is closer to d, we replace the above inequality by an equality,
and what remains is finding the minimum value ct satisfying the given constraints.
By carefully checking the proofs of [Laa15, Appendix A.2-A.3], the exact condition on ct to obtain
the minimum asymptotic time complexity is the following:


ct
−cn = max log2 sin θ2 + π
.
(58)
ρ( 3 , θ2 )
θ2 ∈(0,π)

18

Thijs Laarhoven

Here cn = 21 log2 ( 43 ) ≈ 0.20752, and ρ(θ1 , θ2 ) = ln p(θ1 )/ ln p(θ2 ) corresponds to the exponent ρ for
given angles θ1 , θ2 . Note that in the above equation, only ct is an unknown. Substituting the asymptotic
collision probabilities from Theorem 1, we find a solution at ct ≈ 0.11464, with maximizing angle
θ2 ≈ 0.45739π. This corresponds to a time and space complexity of (n · t)1+o(1) = 2(cn +ct )d+o(d) ≈
20.32216d+o(d) as claimed.

