Polytopes, lattices, and spherical codes for the nearest neighbor problem

TL;DR

This paper develops a framework for analyzing spherical code-based locality-sensitive hashing methods for the nearest neighbor problem, providing exact performance results and identifying optimal codes in low dimensions.

Contribution

It introduces a unified analysis framework for spherical code hashing, deriving exact collision probabilities, and identifies optimal codes in 2-5 dimensions with practical implications.

Findings

Equilateral triangle is optimal in 2D, outperforming hyperplane and cross-polytope hashing.

In 3D and 4D, tetrahedron, 5-cell, and 16-cell are best performers.

In higher dimensions, orthoplices and root lattices like D_k outperform other codes.

Abstract

We study locality-sensitive hash methods for the nearest neighbor problem for the angular distance, focusing on the approach of first projecting down onto a low-dimensional subspace, and then partitioning the projected vectors according to Voronoi cells induced by a suitable spherical code. This approach generalizes and interpolates between the fast but suboptimal hyperplane hashing of Charikar [STOC'02] and the asymptotically optimal but practically often slower hash families of Andoni-Indyk [FOCS'06], Andoni-Indyk-Nguyen-Razenshteyn [SODA'14] and Andoni-Indyk-Laarhoven-Razenshteyn-Schmidt [NIPS'15]. We set up a framework for analyzing the performance of any spherical code in this context, and we provide results for various codes from the literature, such as those related to regular polytopes and root lattices. Similar to hyperplane hashing, and unlike cross-polytope hashing, our analysis of collision probabilities and query exponents is exact and does not hide order terms which vanish only for large $d$, facilitating an easy parameter selection. For the two-dimensional case, we derive closed-form expressions for arbitrary spherical codes, and we show that the equilateral triangle is optimal, achieving a better performance than the two-dimensional analogues of hyperplane and cross-polytope hashing. In three and four dimensions, we numerically find that the tetrahedron, $5$-cell, and $16$-cell achieve the best query exponents, while in five or more dimensions orthoplices appear to outperform regular simplices, as well as the root lattice families $A_k$ and $D_k$. We argue that in higher dimensions, larger spherical codes will likely exist which will outperform orthoplices in theory, and we argue why using the $D_k$ root lattices will likely lead to better results in practice, due to a better trade-off between the asymptotic query exponent and the concrete costs of hashing.

Tables

Table 1. Table 1 : Parameters ρ 𝜌 \rho for various spherical codes. Sphere packings are from  [ 45 ] . Rows labeled “spherical caps” are lower bounds (Theorem  1.4 ). Superscripts refer to the number of caps minimizing ρ θ subscript 𝜌 𝜃 \rho_{\theta} . Bold values indicate the best values ρ 𝜌 \rho encountered up to this dimension. Results for k = 1 , 2 , d 𝑘 1 2 𝑑 k=1,2,d , were obtained analytically; for k = 3 , 4 , 5 , 6 𝑘 3 4 5 6 k=3,4,5,6 most results were obtained through numerical integration and Monte Carlo simulation.

$k$	spherical code	$c$	$\rho \left(\frac{\pi }{12}\right)$	$\rho \left(\frac{\pi }{6}\right)$	$\rho \left(\frac{\pi }{4}\right)$	$\rho \left(\frac{\pi }{3}\right)$	$\rho \left(\frac{5\pi }{12}\right)$
$1$	spherical caps		${\text{<em>0.1255</em>}}^{\left(2\right)}$	${\text{<em>0.2630</em>}}^{\left(2\right)}$	${\text{<em>0.4150</em>}}^{\left(2\right)}$	${\text{<em>0.5850</em>}}^{\left(2\right)}$	${\text{<em>0.7776</em>}}^{\left(2\right)}$
	hyperplane	$2$	$0.1255$	$0.2630$	$0.4150$	$0.5850$	$0.7776$
$2$	spherical caps		${\text{<em>0.1194</em>}}^{\left(3\right)}$	${\text{<em>0.2518</em>}}^{\left(3\right)}$	${\text{<em>0.4005</em>}}^{\left(3\right)}$	${\text{<em>0.5700</em>}}^{\left(3\right)}$	${\text{<em>0.7666</em>}}^{\left(3\right)}$
	triangle ($S_2$)	$3$	$0.1194$	$0.2518$	$0.4005$	$0.5700$	$0.7666$
	square ($O_2$, $C_2$, $D_2$)	$4$	$0.1255$	$0.2630$	$0.4150$	$0.5850$	$0.7776$
	pentagon	$5$	$0.1343$	$0.2788$	$0.4346$	$0.6040$	$0.7905$
	hexagon ($A_2$)	$6$	$0.1438$	$0.2954$	$0.4544$	$0.6222$	$0.8022$
$3$	spherical caps		${\text{<em>0.1117</em>}}^{\left(5\right)}$	${\text{<em>0.2389</em>}}^{\left(4\right)}$	${\text{<em>0.3846</em>}}^{\left(4\right)}$	${\text{<em>0.5548</em>}}^{\left(4\right)}$	${\text{<em>0.7561</em>}}^{\left(4\right)}$
	tetrahedron ($S_3$)	$4$	$0.1155$	$0.2445$	$0.3910$	$0.5600$	$0.7592$
	sphere packing	$5$	$0.1170$	$0.2481$	$0.3965$	$0.5664$	$0.7644$
	octahedron ($O_3$)	$6$	$0.1159$	$0.2465$	$0.3952$	$0.5661$	$0.7649$
	sphere packing	$7$	$0.1207$	$0.2554$	$0.4065$	$0.5772$	$0.7725$
	cube ($C_3$)	$8$	$0.1255$	$0.2630$	$0.4150$	$0.5850$	$0.7776$
	sphere packing	$9$	$0.1217$	$0.2591$	$0.4129$	$0.5850$	$0.7786$
	icosahedron	$12$	$0.1255$	$0.2678$	$0.4255$	$0.5983$	$0.7883$
	cuboctahedron ($A_3$, $D_3$)	$12$	$0.1294$	$0.2728$	$0.4301$	$0.6017$	$0.7900$
	dodecahedron	$20$	$0.1509$	$0.3077$	$0.4692$	$0.6360$	$0.8112$
$4$	spherical caps		${\text{<em>0.1050</em>}}^{\left(7\right)}$	${\text{<em>0.2285</em>}}^{\left(6\right)}$	${\text{<em>0.3730</em>}}^{\left(6\right)}$	${\text{<em>0.5433</em>}}^{\left(5\right)}$	${\text{<em>0.7466</em>}}^{\left(5\right)}$
	$5$-cell ($S_4$)	$5$	$0.1126$	$0.2392$	$0.3840$	$0.5527$	$0.7537$
	sphere packing	$6$	$0.1128$	$0.2401$	$0.3861$	$0.5555$	$0.7564$
	sphere packing	$7$	$0.1120$	$0.2392$	$0.3852$	$0.5553$	$0.7567$
	$16$-cell ($O_4$)	$8$	$0.1107$	$0.2368$	$0.3822$	$0.5528$	$0.7553$
	sphere packing	$10$	$0.1136$	$0.2426$	$0.3909$	$0.5623$	$0.7622$
	sphere packing	$13$	$0.1133$	$0.2439$	$0.3939$	$0.5666$	$0.7663$
	tesseract ($C_4$)	$16$	$0.1255$	$0.2630$	$0.4150$	$0.5850$	$0.7776$
	runcinated $5$-cell ($A_4$)	$20$	$0.1216$	$0.2586$	$0.4128$	$0.5855$	$0.7795$
	octacube ($D_4$)	$24$	$0.1202$	$0.2577$	$0.4140$	$0.5877$	$0.7823$
$5$	spherical caps		${\text{<em>0.0997</em>}}^{\left(13\right)}$	${\text{<em>0.2203</em>}}^{\left(10\right)}$	${\text{<em>0.3630</em>}}^{\left(8\right)}$	${\text{<em>0.5342</em>}}^{\left(7\right)}$	${\text{<em>0.7415</em>}}^{\left(6\right)}$
	$5$-simplex ($S_5$)	$6$	$0.1105$	$0.2354$	$0.3785$	$0.5469$	$0.7493$
	$5$-orthoplex ($O_5$)	$10$	$0.1076$	$0.2299$	$0.3733$	$0.5433$	$0.7483$
	$1_{21}$-polytope	$16$	$0.1080$	$0.2330$	$0.3794$	$0.5516$	$0.7554$
	expanded $5$-simplex ($A_5$)	$30$	$0.1167$	$0.2494$	$0.4007$	$0.5735$	$0.7713$
	$5$-hypercube ($C_5$)	$32$	$0.1255$	$0.2630$	$0.4150$	$0.5850$	$0.7776$
	rectified $5$-orthoplex ($D_5$)	$40$	$0.1139$	$0.2471$	$0.4009$	$0.5757$	$0.7739$
$6$	spherical caps		${\text{<em>0.0955</em>}}^{\left(18\right)}$	${\text{<em>0.2135</em>}}^{\left(15\right)}$	${\text{<em>0.3552</em>}}^{\left(11\right)}$	${\text{<em>0.5263</em>}}^{\left(9\right)}$	${\text{<em>0.7357</em>}}^{\left(8\right)}$
	$6$-simplex ($S_6$)	$7$	$0.1089$	$0.2319$	$0.3742$	$0.5422$	$0.7458$
	$6$-orthoplex ($O_6$)	$12$	$0.1065$	$0.2260$	$0.3670$	$0.5361$	$0.7431$
	$2_{21}$-polytope	$27$	$0.1038$	$0.2258$	$0.3712$	$0.5442$	$0.7510$
	$1_{31}$-polytope	$32$	$0.1062$	$0.2314$	$0.3788$	$0.5520$	$0.7564$
	expanded $6$-simplex ($A_6$)	$42$	$0.1133$	$0.2427$	$0.3917$	$0.5642$	$0.7647$
	rectified $6$-orthoplex ($D_6$)	$60$	$0.1107$	$0.2404$	$0.3915$	$0.5661$	$0.7673$
	hypercube ($C_6$)	$64$	$0.1255$	$0.2630$	$0.4150$	$0.5850$	$0.7776$
$d$	lower bound		0.0173	0.0718	0.1716	0.3333	0.5888
	simplex ($S_d$)	$d+1$	$0.0173$	$0.0718$	$0.1716$	$0.3333$	$0.5888$
	orthoplex ($O_d$)	$2d$	$0.0173$	$0.0718$	$0.1716$	$0.3333$	$0.5888$
	expanded simplex ($A_d$)	$d\left(d+1\right)$	$0.0173$	$0.0718$	$0.1716$	$0.3333$	$0.5888$
	rectified orthoplex ($D_d$)	$2d\left(d-1\right)$	$0.0173$	$0.0718$	$0.1716$	$0.3333$	$0.5888$
	hypercube ($C_d$; $𝐀$ orth.)	$2^d$	$0.0799$	$0.1800$	$0.3151$	$0.5201$	$0.7686$
	hypercube ($C_d$; $𝐀$ Gaussian)	$2^d$	$0.1255$	$0.2630$	$0.4150$	$0.5850$	$0.7776$

Equations

$$p_1$$

$$\displaystyle\rho=\ln\left[\frac{1}{3}+3\left(\frac{\pi-\theta}{2\pi}\right)^{2}-3\left(\frac{\arccos(\frac{1}{2}\cos\theta)}{2\pi}\right)^{2}\right]\Big{/}\ln\left[\frac{1}{3}\right].$$

$$\rho (\theta )\ge {\rho }_k(\theta ):=\underset{c\ge 2}{\min }{\rho }_k(c,\theta ),$$

$$\displaystyle\rho_{k}(c,\theta):=\log\left\{\Pr_{\boldsymbol{x},\boldsymbol{y}\sim\mathcal{N}(0,1)^{k}}\left(\frac{x_{1}}{\|\boldsymbol{x}\|}\geq\alpha_{k}(c),\frac{x_{1}\cos\theta+y_{1}\sin\theta}{\|\boldsymbol{x}\cos\theta+\boldsymbol{y}\sin\theta\|}\geq\alpha_{k}(c)\right)\right\}\bigg{/}\log\left(\frac{1}{c}\right).$$

$$\varphi (\mathbf{x},\mathbf{y})=\arccos ⟨\frac{\mathbf{x}}{\Vert \mathbf{x}\Vert },\frac{\mathbf{y}}{\Vert \mathbf{y}\Vert }⟩.$$

$$\frac{\mu ({\mathcal{C}}_{\alpha })}{\mu ({\mathcal{S}}^{k-1})}=\frac{1}{2}{I}_{1-{\alpha }^2}\left(\frac{k-1}{2},\frac{1}{2}\right).$$

$${\mathcal{V}}_i={\mathcal{V}}_i(\mathcal{C})=\left\{\mathbf{x}\in {\mathbb{R}}^k:\Vert \mathbf{x}-{\mathbf{c}}_i\Vert \le \Vert \mathbf{x}-{\mathbf{c}}_j\Vert ,\forall j=1,\dots ,c\right\}.$$

$$\text{dist}(\mathbf{x},\mathbf{y})=\varphi (\mathbf{x},\mathbf{y})=\arccos ⟨\frac{\mathbf{x}}{\Vert \mathbf{x}\Vert },\frac{\mathbf{y}}{\Vert \mathbf{y}\Vert }⟩.$$

$$\mathbf{x},\mathbf{y}\in {\mathcal{S}}^{d-1}⟹\frac{1}{2}\Vert \mathbf{x}-\mathbf{y}{\Vert }^2=1-\cos \varphi (\mathbf{x},\mathbf{y}).$$

$$\underset{\begin{array}{c}h\sim \mathcal{H}\\ \mathbf{x},\mathbf{y}\sim {\mathcal{S}}^{d-1}\end{array}}{\mathrm{Pr}}(h(\mathbf{x})=h(\mathbf{y})\mid \varphi (\mathbf{x},\mathbf{y})\le \theta )\ge p_1,\underset{\begin{array}{c}h\sim \mathcal{H}\\ \mathbf{x},\mathbf{y}\sim {\mathcal{S}}^{d-1}\end{array}}{\mathrm{Pr}}(h(\mathbf{x})=h(\mathbf{y}))\le p_2.$$

$$\rho =\frac{\log 1/p_1}{\log 1/p_2},\ell =\frac{\log n}{\log 1/p_2},t=O(n^{\rho }).$$

$$h_{\mathrm{A}}(\mathbf{x})=\mathrm{dec}(\mathrm{A}\mathbf{x}).$$

$$p_1$$

$$=\sum _{i=1}^c\underset{{\mathbf{a}}_1,{\mathbf{a}}_2\sim \mathcal{N}(0,{\mathrm{I}}_k)}{\mathrm{Pr}}(\mathrm{dec}({\mathbf{a}}_1)=\mathrm{dec}({\mathbf{a}}_1\cos \theta +{\mathbf{a}}_2\sin \theta )={\mathbf{c}}_i).$$

$$p_2$$

$$\left[\begin{array}{c}\text{— }({\mathbf{c}}_i-{\mathbf{c}}_1)\text{ —}\\ \text{— }({\mathbf{c}}_i-{\mathbf{c}}_2)\text{ —}\\ \dots \\ \text{— }({\mathbf{c}}_i-{\mathbf{c}}_c)\text{ —}\end{array}\right]\cdot \left[\begin{array}{c}|\\ \mathbf{x}\\ |\end{array}\right]\ge \left[\begin{array}{c}|\\ 0\\ |\end{array}\right].$$

$$p_1$$

$$=\sum _{i=1}^c\underset{\mathbf{z}\sim \mathcal{N}(0,{\mathrm{I}}_{2k})}{\mathrm{Pr}}\left(\left\{\left[\begin{array}{cc}1& 0\\ \cos \theta & \sin \theta \end{array}\right]\otimes {\mathrm{C}}_i\right\}\mathbf{z}\ge 0\right).$$

$$p_1$$

$$p_2$$

$$=\sum _{i=1}^c\underset{\mathbf{z}\sim \mathcal{N}(0,{\mathrm{\Theta }}_i)}{\mathrm{Pr}}{\left(\mathbf{z}\ge 0\right)}^2,{\mathrm{\Theta }}_i={\mathrm{C}}_i{\mathrm{C}}_i^T,$$

$$h(\mathbf{x})=\mathrm{dec}(\mathrm{A}\mathbf{x})\in \{1,\dots ,c\},\mathrm{A}\sim \mathcal{N}(0,1{)}^{k\times d}.$$

$$p_1$$

$$\mathrm{dec}(\mathbf{x})=i⟺\forall {\mathbf{c}}_j\in \mathcal{C}:\Vert \mathbf{x}-{\mathbf{c}}_i\Vert \le \Vert \mathbf{x}-{\mathbf{c}}_j\Vert ⟺\forall {\mathbf{c}}_j\in \mathcal{C}:⟨\mathbf{x},{\mathbf{c}}_i-{\mathbf{c}}_j⟩\ge 0.$$

$$\mathrm{dec}(\mathbf{x})=i⟺\forall {\mathbf{c}}_j\in {\mathcal{R}}_i:\Vert \mathbf{x}-{\mathbf{c}}_i\Vert \le \Vert \mathbf{x}-{\mathbf{c}}_j\Vert ⟺\forall {\mathbf{c}}_j\in {\mathcal{R}}_i:⟨\mathbf{x},{\mathbf{c}}_i-{\mathbf{c}}_j⟩\ge 0.$$

$$p_1$$

$$\mathrm{\Sigma }=\left[\begin{array}{cc}1& a\\ a& 1\end{array}\right]⟹\underset{(z_1,z_2)\sim \mathcal{N}(0,\mathrm{\Sigma })}{\mathrm{Pr}}\left(z_1\ge 0,z_2\ge 0\right)=\frac{\pi -\arccos a}{2\pi }.$$

$$p_1=1-\frac{\theta }{\pi },p_2=\frac{1}{2}.$$

$$p_1$$

$${\mathrm{\Sigma }}_1=2(1-\cos \frac{2\pi }{c})\cdot \left[\begin{array}{cc}1& \cos \theta \\ \cos \theta & 1\end{array}\right]\otimes \left[\begin{array}{cc}1& -\cos \frac{2\pi }{c}\\ -\cos \frac{2\pi }{c}& 1\end{array}\right].$$

Full text

For the two-dimensional case, we analytically derive closed-form expressions for arbitrary spherical codes, and we show that the equilateral triangle is optimal, achieving a better performance than the two-dimensional analogues of hyperplane and cross-polytope hashing. In three and four dimensions, we numerically find that the tetrahedron and $5$-cell (the $3$-simplex and $4$-simplex) and the $16$-cell (the $4$-orthoplex) achieve the best query exponents, while in five or more dimensions orthoplices appear to outperform regular simplices, as well as the root lattice families $A_{k}$ and $D_{k}$ in terms of minimizing the query exponent. We provide lower bounds based on spherical caps, and we predict that in higher dimensions, larger spherical codes exist which outperform orthoplices in terms of the query exponent, and we argue why using the $D_{k}$ root lattices will likely lead to better results in practice as well (compared to using cross-polytopes), due to a better trade-off between the asymptotic query exponent and the concrete costs of hashing.

1 Introduction

Given a large database of high-dimensional vectors, together with a target data point which does not lie in this database, a natural question to ask is: which item in the database is the most similar to the query? And can we somehow preprocess and store the database in a data structure that allows such queries to be answered faster? These and related questions have long been studied in various contexts, such as machine learning, coding theory, pattern recognition, and cryptography [25, 42, 16, 24, 31, 38], under the headers of similarity search and nearest neighbor searching. Observe that a naive solution might consist of simply storing the data set in a big list, and to search this list in linear time to find the nearest neighbor to a query point. This solution requires an amount of time and space scaling linearly in the size of the data set, and solutions we are interested in commonly require more preprocessed space (and time), but achieve a sublinear query time complexity to find the nearest neighbor.

Depending on the context of the problem, different solutions for these problems have been proposed and studied. For the case where the dimensionality of the original problem is constant, efficient solutions are known to exist [9]. Throughout the remainder of the paper, we will therefore assume that the dimensionality of the data set is superconstant. In the late 1990s, Indyk–Motwani [28] proposed the locality-sensitive hashing framework, and until today this method remains one of the most prominent and popular methods for nearest neighbor searching in high-dimensional vector spaces, both due to its asymptotic performance when using theoretically optimal hash families [5, 4], and due to its practical performance when instantiated with truly efficient locality-sensitive hash functions [1, 19, 15]. And whereas many other methods scale poorly as the dimensionality of the problem increases, locality-sensitive hashing remains competitive even in high-dimensional settings.

Although solutions for both asymptotic and concrete settings have been studied over the years, there is often a separation between both worlds: some methods work well in practice but do not scale optimally when the parameters increase (e.g. Charikar’s hash family [19]); while some other methods are known to achieve a superior performance for sufficiently large parameter sizes, but may not be quite as practical in some applications due to large hidden order terms in the asymptotic analysis (e.g. hash families studied in [3, 5, 4] and filter families from [13, 6, 22]). Moreover, the latter methods are often not as easy to deploy in practice due to these unspecified hidden order terms, making it hard to choose the scheme parameters that optimize the performance. A key problem in this area thus remains to close the gap between theory and practice, and to offer solutions that interpolate between quick–and–dirty simple approaches that might not scale well with the problem size, and more sophisticated methods that only start outperforming these simpler methods as the parameters are sufficiently large.

1.1 Related work

In this paper we will focus on methods for solving the nearest neighbor problem for the angular distance, which is closely related to the nearest neighbor problem in various $\ell_{p}$-norms: as shown by Andoni and Razenshteyn [7], a solution for the nearest neighbor problem for the angular distance (or the Euclidean distance on the sphere) can be optimally extended to a solution for the $\ell_{2}$-norm for all of $\mathbb{R}^{d}$. Solutions for the $\ell_{2}$-norm can further be translated to e.g. solutions for the $\ell_{1}$-norm via appropriate embeddings.

For the angular distance, perhaps the most well-known and widely deployed approach for finding similar items in a large database is to use the hyperplane hash family of Charikar [19]. For spherically-symmetric, random data sets on the sphere of size $n$, it can find a nearest neighbor at angle at most e.g. $\theta=\frac{\pi}{3}$ in sublinear time $\tilde{O}(n^{\rho})$ and space $\tilde{O}(n^{1+\rho})$, with $\rho\approx 0.5850$. Various improvements later showed that smaller query exponents $\rho$ can be achieved in sufficiently high dimensions [3, 5], the most practical of which is based on cross-polytopes or orthoplices [47, 4, 30]: for large $d$ and for the same target angle $\theta=\frac{\pi}{3}$, the query time complexity scales as $n^{\rho+o(1)}$ with $\rho=1/3$. Note that the convergence to the limit is rather slow [4, Theorem 1] and depends on both $d$ and $n$ being sufficiently large – for fixed $d$ and large $n$, the exponent is still larger than $1/3$. In certain practical applications with moderate-sized data sets, hyperplane hashing still appears to outperform cross-polytope hashing and other advanced hashing and filtering schemes [36, 37, 14, 2, 13] due to its low cost for hashing, and the absence of lower order terms in the exponent.

Related to the topic of this paper, various other works have further observed the relation between finding good locality-sensitive hash families for the angular distance and finding “nice” spherical codes that partition the sphere well, and allow for efficient decoding [3, 47, 48, 4, 13, 18]. The requirements on a spherical code to be a good spherical code are somewhat intricate to state, and to date it is an open problem to exactly quantify which spherical codes are the most suited for nearest neighbor searching. It may well be possible to improve upon the state-of-the-art orthoplex (cross-polytope) locality-sensitive hash family [4, 40] in practice with a method achieving the same asymptotic scaling, but with a faster convergence to the limit, and thus potentially a better performance in practice for large problem instances.

1.2 A framework for evaluating spherical codes

As a first contribution of this paper, we provide a framework for analyzing the performance of arbitrary spherical codes in the context of locality-sensitive hashing for the angular distance, where we focus on the approach of (1) projecting down onto a random low-dimensional subspace, and (2) partitioning the resulting subspace according to the Voronoi cells induced by a spherical code. More specifically, we relate the collision probabilities appearing in the analysis of these hash functions to so-called orthant probabilities of multivariate normal distributions with non-trivial correlation matrices. Below we informally state this relation for general spherical codes, which provides us a recipe for computing the collision probabilities (and the query exponent) as a function of the set of vertices $\mathcal{C}$ of the corresponding spherical code. Here a hash family being $(\theta,p_{1},p_{2})$-sensitive means that uniformly random vectors on the sphere collide with probability at most $p_{2}$, and target nearest neighbors at angle at most $\theta$ from a random query vector collide in a random hash function with probability at least $p_{1}$.

Theorem 1.1** (Spherical code locality-sensitive hashing).**

Let $\mathcal{C}\subset\mathcal{S}^{k-1}$ be a $k$-dimensional spherical code, and consider a hash family where:

• •

We first project onto a $k$-dimensional subspace using a random matrix $\boldsymbol{\mathord{\mathrm{A}}}\sim\mathcal{N}(0,1)^{k\times d}$;

• •

We then assign hash values based on which $\boldsymbol{c}\in\mathcal{C}$ is nearest to the projected vector.

Then for any $\theta\in(0,\frac{\pi}{2})$, this family is $(\theta,p_{1},p_{2})$-sensitive, where $p_{2}$ can be expressed in terms of the relative volumes of the Voronoi cells induced by $\mathcal{C}$, and $p_{1}$ can be expressed as a sum of orthant probabilities $\Pr_{\boldsymbol{z}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\mathord{\mathrm{\Sigma}}}_{i})}\left(\boldsymbol{z}\geq\boldsymbol{0}\right)$, where each $\boldsymbol{\mathord{\mathrm{\Sigma}}}_{i}$ has size $(2k)\times(2k)$.

The locality-sensitive hashing exponent $\rho=\log p_{1}/\log p_{2}$ describes the distinguishing power of a hash family, as for large $n$ this tells us that we can solve the (average-case) nearest neighbor problem with target angle $\theta$ in query time $\tilde{O}(n^{\rho})$ with space $\tilde{O}(n^{1+\rho})$. The theorem above essentially provides us a recipe which, given any spherical code $\mathcal{C}$ and target angle $\theta$ as input, tells us how to compute these probabilities $p_{1}$ and $p_{2}$ (and thus $\rho$) exactly.

While the above theorem is somewhat abstract, we show how to explicitly construct the correlation matrices $\boldsymbol{\mathord{\mathrm{\Sigma}}}_{i}$ appearing in the theorem, allowing us to always at least obtain numerical estimates on the performance of different spherical codes in the context of nearest neighbor searching. We further study how to reduce the dimensionality of the problem (and in particular, the sizes of the matrices $\boldsymbol{\mathord{\mathrm{\Sigma}}}_{i}$) when the spherical code exhibits many symmetries, such as being isogonal, and we show how using only the relevant vectors of each code word can further simplify the computations. In most cases the resulting orthant probabilities will still remain too complex to evaluate analytically, but in some cases we can obtain exact expressions this way.

1.3 A survey of known spherical codes

One and two dimensions.

Using this new framework, we then apply it to various spherical codes, starting in very low dimensions. For the one-dimensional case we rederive the celebrated result of Charikar [19] for one-dimensional spherical codes, noting that the collision probability analysis in fact dates back to an old result from the late 1800s [44]. Applying the same framework to two-dimensional codes, among others we establish a general formula for collision probabilities and query exponents for arbitrary polygons, as well as a more compact description of the collision probabilities for regular polygons.

Theorem 1.2** (Collision probabilities for regular polygons).**

Let $\mathcal{C}\subset\mathcal{S}^{1}$ consist of the vertices of the regular $c$-gon, for $c\geq 2$, and let $\theta\in(0,\frac{\pi}{2})$. Then the corresponding project–and–partition hash family $\mathcal{H}$ is $(\theta,p_{1},p_{2})$-sensitive, with:

[TABLE]

As a direct corollary of the above theorem, we establish that triangular locality-sensitive hashing achieves a superior asymptotic performance to hyperplane hashing, and achieves the lowest query exponents $\rho$ among all regular polygons.

Corollary 1.3** (Triangular locality-sensitive hashing).**

Consider the hash family where we first project onto a random plane using a random projection matrix $\boldsymbol{\mathord{\mathrm{A}}}\sim\mathcal{N}(0,1)^{2\times d}$, and then decode to the nearest corner of a fixed equilateral triangle centered at $(0,0)$. Then, for target angles $\theta\in(0,\frac{\pi}{2})$, this hash family achieves query exponents $\rho$ of the form:

[TABLE]

Among all such project–and–partition hash families based on regular $k$-gons, this family achieves the lowest values $\rho$ for any target angle $\theta\in(0,\frac{\pi}{2})$.

Note again that the above statement of $\rho$ is exact, and does not hide any order terms in $d$ or $n$ – in fact, the collision probabilities do not depend on $d$ at all, due to our choice of $\boldsymbol{\mathord{\mathrm{A}}}$ being Gaussian. Further note that the $2$-gon in the plane corresponds to the one-dimensional antipodal code of Charikar, and triangular hashing therefore strictly improves upon hyperplane hashing for any $\theta$ in terms of the query exponent $\rho$. For instance, we can find neighbors at angle at most $\frac{\pi}{3}$ in time $\tilde{O}(n^{\rho})$ with $\rho\approx 0.56996$, offering a concrete but minor improvement over the hyperplane hashing approach of Charikar [19] with query exponent $\rho\approx 0.5850$. This is the best we can do with any two-dimensional isogonal spherical code, and we numerically predict that this code achieves the lowest values $\rho$ among all (not necessarily isogonal) two-dimensional spherical codes as well.

Three and four dimensions.

For three-dimensional codes, through numerical integration of the resulting orthant probabilities we conclude that the tetrahedron appears to minimize the query exponent $\rho$ out of all three-dimensional codes, beating the three-dimensional analogues of e.g. the cross-polytope and hypercube, as well as various sphere packings and other regular polytopes appearing in the literature, such as the Platonic and Archimedean solids. In four dimensions an interesting phenomenon occurs: the so-called $5$-cell ($4$-simplex) and $16$-cell ($4$-orthoplex) are optimal in different regimes, with the $5$-cell inducing a more coarse-grained partition of the space and achieving a better performance when the nearest neighbor lies relatively far away from the data set, and the $16$-cell inducing a more fine-grained partition, and working better when the nearest neighbor lies relatively close to the target vector. This crossover effect is also visualized in Figure 1, and it strengthens our intuition that as the dimensionality goes up, or as the nearest neighbor lies closer to the query, more fine-grained partitions are necessary to obtain the best performance. An overview of some of the exponents $\rho$ for various low-dimensional spherical codes, for different target angles $\theta\in\frac{\pi}{12}\{1,2,3,4,5\}$, is given in Table 1. The best exponents $\rho$ are highlighted in bold.

Five and more dimensions.

For higher-dimensional spherical codes, we obtain further improvements in the query exponents $\rho$ through the use of suitable spherical codes, as shown in Table 1 and Figures 1–4. For dimensions $5$ and $6$, the corresponding orthoplices achieve a better performance than the simplices, and the exotic polytopes $1_{21}$ and $2_{21}$, related to the root lattice $E_{6}$, seem useful in the context of nearest neighbor searching as well.

We further study the performance of the following five non-trivial infinite families of spherical codes. Color codes below correspond to the same colors used in Table 1 and Figures 1–4 to differentiate these families of codes.

• •

The simplices $S_{k}$ on $k+1$ vertices;

• •

The orthoplices $O_{k}$ on $2k$ vertices (also known as cross-polytopes [48, 4]);

• •

The hypercubes $C_{k}$ on $2^{k}$ vertices (as studied in [33]);

• •

The expanded simplices $A_{k}$ on $k(k+1)$ vertices;

• •

The rectified orthoplices $D_{k}$ on $2k(k-1)$ vertices.

The latter two families are connected to the root lattices $A_{k}$ and $D_{k}$, while the first three are related to the lattice $\mathbb{Z}^{k}$. We conjecture that, should other nice families of dense lattices be found (e.g. the recent [49]), these may give rise to suitable spherical codes in nearest neighbor applications as well. For each of the above five families we give closed-form expressions on the correlation matrices $\boldsymbol{\mathord{\mathrm{\Sigma}}}$, but the resulting orthant probabilities that need to be evaluated for computing $\rho$ do not appear to admit simple closed-form expressions. Apart from the family of hypercubes, these are all asymptotically optimal (with $\rho\to(1-\cos\theta)/(1+\cos\theta)$ as $k\to\infty$), although the convergence to the limit may differ for each family.

1.4 Lower bounds via spherical caps

Although this work tries to be exhaustive in covering as many (families of) spherical codes as possible, better spherical codes may exist, achieving even lower query exponents $\rho$. As these codes may be hard to find, and as it may be difficult to rule out the existence of other, better spherical codes, the next best thing one might hope for is a somewhat tight lower bound on the performance of any $k$-dimensional spherical code in our framework, which hopefully comes close to the performance of the spherical codes we have considered in this survey.

As has been established in several previous works on nearest neighbor searching on the sphere [3, 5, 13, 32, 6, 22], ideally we would like the hash regions induced by the partitions to take the shape of a spherical cap. Such a shape minimizes the angular radius, given that the region has a fixed volume, and in a sense it is the most natural and smoothest shape that a region on a sphere can take. So in a utopian world, one might hope that a spherical code partitions the sphere into $c$ regions, and each region corresponds exactly to a spherical cap of volume $\operatorname{Vol}(\mathcal{S}^{k-1})/c$. Clearly such spherical codes do not exist for $k>2$ and $c>2$, but such an extremal example does give us an indication on the limits of what might be achievable in dimension $k$, and how the optimal $\rho$ decreases with $k$. Note that random spherical codes might approach this utopian setting in high dimensions.

Following the above reasoning, and using a classic result of Baernstein–Taylor [12], we state a formal lower bound on the performance parameter $\rho$ of any $k$-dimensional spherical code of size $c$, and by minimizing over $c$ for given $k$, one obtains a lower bound for any spherical code living in $k$ dimensions. Here $I_{x}(a,b)$ denotes the regularized incomplete beta function, which comes from computing the volume of a sphere in $k$ dimensions, while $\|\cdot\|$ denotes the Euclidean norm.

Theorem 1.4** (Spherical cap lower bounds).**

Let $\mathcal{C}\subset\mathcal{S}^{k-1}$ be a spherical code, and let $\mathcal{H}$ be the associated project–and–partition hash family. Then the parameter $\rho$ for $\mathcal{C}$, for target angle $\theta\in(0,\frac{\pi}{2})$, must satisfy:

[TABLE]

where, with $\alpha_{k}(c)$ denoting the solution $\alpha$ to $\frac{1}{2}\cdot I_{1-\alpha^{2}}(\frac{k-1}{2},\frac{1}{2})=\frac{1}{c}$, $\rho_{k}(c,\theta)$ is given by:

[TABLE]

The above minimization over $c\geq 2$ concerns the possible code sizes, and $\rho_{k}(c,\theta)$ describes the parameter $\rho$ one would obtain when indeed, the code consisted of $c$ equivalent regions of equal volume, and each region was shaped like a spherical cap. Equivalently, one could state that any spherical code of size exactly $c$ must satisfy $\rho(\theta)\geq\rho_{k}(c,\theta)$.

Numerical evaluation of these expressions $\rho_{k}(c,\theta)$, and the resulting minimization over $c$, leads to the values in Table 1 in the rows indicated by spherical caps. The superscripts denote the values $c$ that numerically solve the minimization problems. These results in low dimensions suggest that, especially for small angles, the optimal code size should increase superlinearly with the dimension. This inspires the following conjecture, stating that the use of orthoplices is likely not optimal in higher dimensions.

Conjecture 1.5** (Orthoplices are suboptimal for large $k$).**

For arbitrary $\theta\in(0,\frac{\pi}{2})$, there exists a dimension $k_{0}\in\mathbb{N}$ such that, for all dimensions $k\geq k_{0}$, there exist spherical codes $\mathcal{C}\subset\mathcal{S}^{k-1}$ whose query exponents $\rho$ in the project–and–partition framework are smaller than the exponents $\rho$ of the $k$-orthoplex.

Actually finding such spherical codes, or finding families of spherical codes that outperform orthoplices may again be closely related to the problem of finding nice families of dense lattices with efficient decoding algorithms. We informally conjecture that in high dimensions, and for sufficiently large code sizes $c$, random spherical codes may be close to optimal. If we care only about minimizing $\rho$, then a further study might focus on (1) getting a better grip of the optimal scaling of $c=c(k)$ with $k$, and (2) estimating the asymptotic performance of using random spherical codes of size $c(k)$, for large $k$. We predict this will lead to better values $\rho$ than those obtained with cross-polytope hashing.

1.5 Selecting spherical codes in practice

Although the exponent $\rho$, and therefore the probabilities $p_{1}$ and $p_{2}$, directly imply the main performance parameters to assess the asymptotic performance of a locality-sensitive hash family, in practice we are always dealing with concrete, non-asymptotic values $n$, $d$, and $\theta$ – if the convergence to the optimal asymptotic scaling is slow, or if the hash functions are too expensive to evaluate in practice, then hash families with lower query exponents $\rho$ may actually be less practical for small, concrete values of $n$ and $d$ than fast hash families with larger $\rho$. For example, the hash family from [3] may be “optimal”, but appears to be only of theoretical interest. In practice one needs to find a balance between decreasing $\rho$ and using hash functions which are fast to evaluate.

As a more concrete example, consider how hyperplane hashing [19] allows us to partition a sphere in $2^{k}$ regions with a single random projection matrix $\boldsymbol{\mathord{\mathrm{A}}}\in\mathbb{R}^{k\times d}$ (or equivalently $k$ matrices $\boldsymbol{\mathord{\mathrm{A}}}_{1},\dots,\boldsymbol{\mathord{\mathrm{A}}}_{k}\in\mathbb{R}^{1\times d}$ merged into one large matrix). For cross-polytope hashing [4] a $k$-dimensional projection would only divide the $k$-dimensional sphere into $2k$ regions. As the total required number of hash buckets in the data structure is often roughly the same, regardless of the chosen hash family,111While a hash family with lower query exponents $\rho$ does require a smaller number of hash buckets per hash table, this effect is marginal compared to the increase in the number of projections/rotations required by e.g. cross-polytope hashing compared to hyperplane hashing. this means that if we wish to divide the sphere into $2^{k}$ hash regions with cross-polytope hashing, we would need $m=k/\log_{2}(2k)$ independent random projection matrices $\boldsymbol{\mathord{\mathrm{A}}}_{1},\dots,\boldsymbol{\mathord{\mathrm{A}}}_{m}\in\mathbb{R}^{k\times d}$ to hash a vector to one of $2^{k}$ buckets. So even though the resulting partitions for cross-polytope hashing generate smaller exponents $\rho$, in practice this improvement may not offset the additional costs of computing the projections/rotations, which is almost a linear factor $O(d)$ more than for hyperplane hashing. So especially for data sets of small/moderate sizes, hyperplane hashing may be more practical than cross-polytope hashing, as also observed in e.g. [31, 14, 36, 37].

We explicitly quantify the trade-off between the complexity of the hash functions and the asymptotic query exponents $\rho$ in Figure 3, where on the vertical axis we plotted the query exponents $\rho$ as computed in this paper, and on the horizontal axis we plotted the number of bits extracted per row of a projection matrix; for hyperplane hashing we extract $1$ bit per projection, while e.g. for cross-polytope hashing in dimension $k$ we only extract $\log(2k)/k$ bits per projection. Depending on the application, it may be desirable to choose hash families with slightly larger values $\rho$, if this means the cost of the hashing becomes less. Figure 3 makes a case for using hashes induced by the root lattices $D_{k}$, as well as those induced by exotic polytopes such as the $2_{21}$-polytope; compared to e.g. the $5$-orthoplex, the $D_{10}$ lattice achieves lower values $\rho$ and extracts more bits per projection, while the $2_{21}$-polytope further improves upon $D_{10}$ by extracting more hash bits per projection.

To further illustrate how these trade-offs might look in a real-world setting, Figure 4 describes a case study for average-case nearest neighbor searching with different sizes $n$ for the data sets. For small data sets, hyperplane hashing is quite competitive, and the best other spherical codes are those induced by the $D_{k}$ lattices (for small $k$), and spherical codes generated by the polytopes $1_{21},1_{31}$ and $2_{21}$. This case study matches the recommendations from Figure 3. For large $n$, the role of $\rho$ becomes more prominent, and higher-dimensional orthoplices and $D_{k}$ lattices attain the best performance. We further describe how one might choose the best codes in practice.

1.6 Summary and open problems

With the project–and–partition framework outlined in this paper, we can now easily formalize and analyze the performance of any spherical code in this framework, and analyze the collision probabilities and query exponents either analytically (by attempting to simplify the corresponding multivariate orthant probabilities) or numerically (by evaluating these multivariate integrals with mathematical software, such as the mvtnorm package [27]). We observed that already in two dimensions, it is possible to improve upon hyperplane hashing with a smaller exponent $\rho$ by using triangular partitions, and we described closed-form formulas for the resulting parameters $\rho$. In three and four dimensions the improvements in the exponent $\rho$ become more significant, with e.g. the tetrahedron, the $5$-cell, and the $16$-cell achieving the best exponents $\rho$ in these dimensions. For higher dimensions the simplices and orthoplices seem to achieve the best performance in theory, thus making a strong case for the use of these partitions as advertised in [47, 48, 4, 30], and as observed in nearest neighbor benchmarks [15, 10, 40]. The family of $D_{k}$ lattices and the generalized $m$-max hash functions however seem to offer better practical trade-offs between the asymptotic performance in terms of $\rho$ and the costs of computing hashes, as discussed later on.

Finding better spherical codes.

An important open problem, both for our project–and–partition framework and for arbitrary hash families for the angular distance, is to find (if they exist) other nice families of spherical codes which achieve an even better performance than the family of orthoplices. Numerics suggest that as $k$ increases, the optimal code size $c$ should increase faster than the linear scaling of $c=2k$ offered by orthoplices. Finding the optimal scaling of $c$ as a function of $k$, and finding nicer spherical codes closely matching this appropriate scaling of $c$ is left for future work.

Faster pseudorandom projections.

To make e.g. hyperplane and cross-polytope hashing more practical, various ideas were proposed in e.g. [1, 4] to work with sparse projections and pseudorandom rotation matrices. Similar ideas can be used for any hash family, to reduce the cost of multiplication by a random Gaussian matrix $\boldsymbol{\mathord{\mathrm{A}}}$. Concretely, one can often replace it with a sparse, “sufficiently random” pseudorandom matrix $\boldsymbol{\mathord{\mathrm{A}}}$ while still achieving good query exponents $\rho$ in practice. Depending on how these pseudorandom projections are instantiated, this may lead to a different practical evaluation than what we described in our practical case analyses. In particular, as this reduces the relative cost of the hashing, this will further favor schemes which reduce $\rho$ at the cost of increasing the (naive) complexity of hashing.

Using orthogonal projection matrices.

In our framework, we focused on projecting down onto a low-dimensional subspace using Gaussian matrices $\boldsymbol{\mathord{\mathrm{A}}}\sim\mathcal{N}(0,1)^{k\times d}$. For $k\ll d$, using Gaussian matrices or orthogonal matrices (i.e. $\boldsymbol{\mathord{\mathrm{A}}}\boldsymbol{\mathord{\mathrm{A}}}^{T}=I_{k}$) is essentially equivalent [23, 29], but for large $k$ it may be beneficial to use proper rotations induced by orthogonal matrices. For instance, recent work [33] showed that for hypercube hashing in the ambient space, there is a clear gap between using random or orthogonal matrices; using orthogonal matrices generally works better than using random Gaussian projection matrices.

The main issue when analyzing the same framework with orthogonal matrices is that the dependence on $d$ then does not disappear; the distribution of $\boldsymbol{\mathord{\mathrm{A}}}\boldsymbol{x}$ then relies on both $k$ and $d$, rather than only on $\boldsymbol{k}$. Our framework is in a sense dimensionless, as the performance can be computed without knowledge of $d$, and for $k=2$ this even allowed us to obtain explicit analytic expressions for the collision probabilities for arbitrary $k$-gons. Using orthogonal matrices, the collision probabilities will likely become complicated functions of $d$, and comparing different spherical codes may then become a more difficult task. We leave it as an open problem to adjust the above framework to the setting where $\boldsymbol{\mathord{\mathrm{A}}}$ is orthogonal, and to see how much can be gained by using orthogonal rather than Gaussian matrices222Note that one obtains different asymptotics when analyzing the hypercube for Gaussian [19] and orthogonal projection/rotation matrices [33]. For constant $k=O(1)$ and large $d$ both approaches are asymptotically equivalent, but for large $k=O(d)$ we expect orthogonal matrices to give better results..

Outline

The remainder of the paper is structured as follows. In Section 2 we first describe preliminary results and notation. Section 3 defines the model for the nearest neighbor problem, and the framework for analyzing the project–and–partition approach for arbitrary spherical codes. Sections 4–7 continue with applying this framework to spherical codes in dimensions 1 through 4, and Section 8 describes results for spherical codes in higher dimensions. Section 9 continues with (conjectured) lower bounds on the parameter $\rho$ in this framework. Section 10 concludes with a practical analysis, weighing the costs of hashing against the exponent $\rho$.

2 Preliminaries

Basic notation.

Throughout the paper, we denote vectors in lower-case boldface (e.g. $\boldsymbol{x}$), and matrices in upright, upper-case boldface (e.g. $\boldsymbol{\mathord{\mathrm{A}}}$). Coordinates of vectors and matrices are written without boldface (e.g. $x_{i}$ and $\mathord{\mathrm{A}}_{i,j}$). We write $\boldsymbol{x}\geq\boldsymbol{y}$ to denote the event that $x_{i}\geq y_{i}$ for all $i$. We write $\boldsymbol{e}_{i}$ for the $i$th unit vector with a $1$ on position $i$, and a [math] in all other coordinates. We write $\boldsymbol{\mathord{\mathrm{I}}}_{k}$ for the $k\times k$ identity matrix, and $\boldsymbol{\mathord{\mathrm{J}}}_{k}$ for the $k\times k$ matrix with all entries equal to $1$. We write $\boldsymbol{\mathord{\mathrm{K}}}=\boldsymbol{\mathord{\mathrm{A}}}\otimes\boldsymbol{\mathord{\mathrm{B}}}$ for the standard Kronecker (tensor) product of matrices; if $\boldsymbol{\mathord{\mathrm{A}}}\in\mathbb{R}^{r_{a}\times c_{a}}$ and $\boldsymbol{\mathord{\mathrm{B}}}\in\mathbb{R}^{r_{b}\times c_{b}}$, then $\boldsymbol{\mathord{\mathrm{K}}}\in\mathbb{R}^{(r_{a}r_{b})\times(c_{a}c_{b})}$ is a matrix with entries $\mathord{\mathrm{A}}_{i,j}\mathord{\mathrm{B}}_{k,l}$. We write $\|\boldsymbol{x}\|=(\sum_{i}x_{i}^{2})^{1/2}$ for the Euclidean norm of $\boldsymbol{x}$, and we denote the standard dot product in $\mathbb{R}^{k}$ by $\langle\boldsymbol{x},\boldsymbol{y}\rangle=\sum_{i}x_{i}y_{i}$. Given two vectors $\boldsymbol{x},\boldsymbol{y}\in\mathbb{R}^{k}$, we further denote their mutual angle by

[TABLE]

Angles between vectors are invariant under scalar multiplication of either vector (i.e. $\phi(\lambda\boldsymbol{x},\mu\boldsymbol{y})=\phi(\boldsymbol{x},\boldsymbol{y})$ for all $\lambda,\mu>0$), and if $\boldsymbol{\mathord{\mathrm{A}}}$ is a rotation matrix ($\boldsymbol{\mathord{\mathrm{A}}}$ is orthogonal), then $\phi(\boldsymbol{\mathord{\mathrm{A}}}\boldsymbol{x},\boldsymbol{\mathord{\mathrm{A}}}\boldsymbol{y})=\phi(\boldsymbol{x},\boldsymbol{y})$.

Geometry on the sphere.

We write $\mathcal{S}^{k-1}=\{\boldsymbol{x}\in\mathbb{R}^{k}:\|\boldsymbol{x}\|=1\}$ for the Euclidean unit sphere in $k$ dimensions. We write $\mu$ for the standard Haar measure, so that the relative volume of $\mathcal{A}\subseteq\mathcal{S}^{k-1}$ on the sphere can be expressed as $\mu(\mathcal{A})/\mu(\mathcal{S}^{k-1})$. Let $\mathcal{C}_{\alpha}=\{\boldsymbol{x}\in\mathcal{S}^{k-1}:x_{1}>\alpha\}$ denote a spherical cap of height $1-\alpha$. Writing $I_{x}(a,b)$ for the regularized incomplete beta function, we can compute the volume of $\mathcal{C}_{\alpha}$ as [35]:

[TABLE]

Asymptotically (for large $k$), the above ratio scales as $k^{\Theta(1)}\cdot(1-\alpha^{2})^{k/2}$ [43, 13, 6].

Spherical codes and Voronoi cells.

Given a set of points (code words) $\mathcal{C}=\{\boldsymbol{c}_{1},\dots,\boldsymbol{c}_{c}\}\subset\mathbb{R}^{k}$, the Voronoi cells of $\mathcal{C}$ are subsets of $\mathbb{R}^{k}$, associating to each point $\boldsymbol{c}_{i}\in\mathcal{C}$ a Voronoi region

[TABLE]

Note that $\bigcup_{i}\mathcal{V}_{i}=\mathbb{R}^{k}$, while $\sum_{i\neq j}\mu(\mathcal{V}_{i}\cap\mathcal{V}_{j})/\mu(\mathcal{S}^{k-1})=0$; except for the boundaries between cells, which together have measure [math], this forms a proper partition of $\mathbb{R}^{k}$. For a partition of a space into Voronoi cells induced by $\mathcal{C}$, we write $\mathcal{R}_{i}$ for the set of relevant vectors of $\boldsymbol{c}_{i}\in\mathcal{C}$, i.e. the vectors $\boldsymbol{c}_{j}\in\mathcal{C}$ such that $\mathcal{V}_{i},\mathcal{V}_{j}$ share a non-trivial boundary. If $\mathcal{C}\subset\mathcal{S}^{k-1}$, we call $\mathcal{C}$ a spherical code, and in that case the Voronoi cells take conical shapes, and are invariant under (positive) scalar multiplication; $\boldsymbol{x}\in\mathcal{V}_{i}$ iff $\boldsymbol{x}/\|\boldsymbol{x}\|\in\mathcal{V}_{i}$. Each pair of cells $\mathcal{V}_{i},\mathcal{V}_{j}$ then shares a common point $\boldsymbol{0}$; for $\boldsymbol{c}_{j}$ to be considered a relevant vector of $\boldsymbol{c}_{i}$, the intersection between $\mathcal{V}_{i},\mathcal{V}_{j}$ must contain more than only $\boldsymbol{0}$. We further say a spherical code $\mathcal{C}$ is uniform if the associated Voronoi partition satisfies $\mu(\mathcal{V}_{i})=\mu(\mathcal{V}_{j})$ for all $i,j$, and we say $\mathcal{C}$ is isogonal (vertex-transitive) if, for each $\boldsymbol{c}_{i},\boldsymbol{c}_{j}\in\mathcal{C}$, there exists an isometry $f$ with $f(\boldsymbol{c}_{i})=f(\boldsymbol{c}_{j})$ and $f(\mathcal{C})=\mathcal{C}$.

Lattices and root systems.

Given a set $\boldsymbol{\mathord{\mathrm{B}}}=\{\boldsymbol{b}_{1},\dots,\boldsymbol{b}_{k}\}$ of linearly independent vectors, we define the lattice $\mathcal{L}=\mathcal{L}(\boldsymbol{\mathord{\mathrm{B}}})=\{\sum_{i=1}^{k}\lambda_{i}\boldsymbol{b}_{i}:\boldsymbol{\lambda}\in\mathbb{Z}^{k}\}$. A lattice $\mathcal{L}$ is an integral lattice if for all $\boldsymbol{v},\boldsymbol{w}\in\mathcal{L}$ we have $\langle\boldsymbol{v},\boldsymbol{w}\rangle\in\mathbb{Z}$. A lattice $\mathcal{L}$ is a root lattice if it is integral and generated by a set of roots $\boldsymbol{a}\in\mathcal{L}$ satisfying $\|\boldsymbol{a}\|^{2}=2$. Besides the three exceptional lattices $E_{6},E_{7}$ and $E_{8}$, the only irreducible root lattices are those from the following two families [17]:

• •

$A_{k}$: Described as a subset of $\mathbb{R}^{k+1}$, this lattice has $k(k+1)$ roots $\boldsymbol{e}_{i}-\boldsymbol{e}_{j}$ for $i\neq j$.333As all vertices lie on a hyperplane, the vertices can be projected onto $\mathbb{R}^{k}$ preserving mutual distances.

• •

$D_{k}$: This lattice has roots $\pm\boldsymbol{e}_{i}\pm\boldsymbol{e}_{j}$, for $i\neq j$. There are $2k(k-1)$ such roots.

Note that the roots of the root systems $B_{k}$ and $C_{k}$ are of different norms, and these sets of roots therefore do not define proper spherical codes.

Regular polytopes.

Throughout the paper we will study spherical codes generated by regular polytopes, i.e. polytopes with regular facets and regular vertex figures. Besides several regular polytopes in low dimensions, there are only three families of regular polytopes in high dimensions.

• •

$S_{k}$: Described as a subset of $\mathbb{R}^{k+1}$, the $k$-simplex has vertices $\boldsymbol{e}_{i}$ for $i=1,\dots,k+1$. 11footnotemark: 1

• •

$O_{k}$: The $k$-orthoplex or $k$-cross polytope has $2k$ vertices of the form $\pm\boldsymbol{e}_{i}$, for $i=1,\dots,k$.

• •

$C_{k}$: The $k$-hypercube has $2^{k}$ vertices $\pm\boldsymbol{e}_{1}\pm\dots\pm\boldsymbol{e}_{k}$, taking all possible signs.

Besides these infinite families, there are several regular polytopes in dimension up to four, such as the Platonic and Archimedian solids.

Probabilities.

Given a set $\mathcal{A}$, we write $a\sim\mathcal{A}$ to denote the process of drawing $a$ uniformly at random from $\mathcal{A}$. With slight abuse of notation, we similarly write $a\sim\xi$ to denote that $a$ is drawn from a probability distribution $\xi$; from the context it will be clear which case we are referring to. We write $\mathcal{N}(\mu,\sigma^{2})$ for the (univariate) normal distribution with mean $\mu$ and variance $\sigma^{2}$, and we write $\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\mathord{\mathrm{\Sigma}}})$ for the multivariate normal distribution with mean $\boldsymbol{\mu}$ and correlation matrix $\boldsymbol{\mathord{\mathrm{\Sigma}}}$. We refer to $\mathcal{N}(0,1)$ and $\mathcal{N}(\boldsymbol{0},\boldsymbol{\mathord{\mathrm{I}}}_{k})$ as the univariate and multivariate standard normal distributions.

3 A framework for nearest neighbor searching with spherical codes

In this section we will describe a framework for solving a natural average-case version of the nearest neighbor problem with well-chosen spherical codes, and we will show how the performance of these codes can be estimated and compared, to find the best spherical codes. In subsequent sections we will then apply this framework to different spherical codes.

3.1 The nearest neighbor problem

Let us start with a broad definition of the nearest neighbor problem, the topic of this paper.

Definition 3.1** (Nearest neighbor problem).**

Given a data set $\mathcal{D}=\{\boldsymbol{p}_{1},\dots,\boldsymbol{p}_{n}\}\subset\mathbb{R}^{d}$, the nearest neighbor problem asks to index $\mathcal{D}$ in a data structure such that, when later given a query vector $\boldsymbol{q}\in\mathbb{R}^{d}$, one can quickly find a nearest neighbor $\boldsymbol{p}^{*}\in\mathcal{D}$ to $\boldsymbol{q}$.

Depending on the similarity measure, which defines when points are considered nearby, and depending on how the above problem is modeled, different solutions exist. In this paper we will be concerned with the nearest neighbor problem for the angular distance, where the similarity between two vectors is measured by their angle:

[TABLE]

With this notion of similarity, two vectors are considered nearby if they have a small angle, and so for the nearest neighbor problem we are looking for the point $\boldsymbol{p}^{*}\in\mathcal{D}$ which has the smallest angle with the query vector $\boldsymbol{q}$. As the angle between two vectors is independent of their norms, this can equivalently be viewed as the nearest neighbor problem for the Euclidean distance on the unit sphere ($\mathcal{D}\subset\mathcal{S}^{d-1}$) by scaling all vectors to have norm $1$: a small angle then corresponds to a small distance on the unit sphere through the following one-to-one relation between Euclidean distances and angles:

[TABLE]

As described in e.g. [7, 4, 8, 6], the angular nearest neighbor problem (or the Euclidean nearest neighbor problem on the unit sphere) is in a sense the canonical hard problem; through various reductions, a solution to the nearest neighbor problem for this problem would immediately yield an (efficient) algorithm for the nearest neighbor problem in the $\ell_{1}$ and $\ell_{2}$-norms in all of $\mathbb{R}^{d}$. Studying this problem therefore may also be of interest when considering the nearest neighbor problem for other metrics.

Since it has long been known that solving this problem exactly for worst-case instances cannot be done in sublinear time in the size of the data set $n=|\mathcal{D}|$ [28, 8], a large body of work has focused on approximation versions of this problem under worst-case assumptions, where an approximate solution $\hat{\boldsymbol{p}}\in\mathcal{D}$ with $\text{dist}(\hat{\boldsymbol{p}},\boldsymbol{q})\leq c\cdot\text{dist}(\boldsymbol{p}^{*},\boldsymbol{q})$ suffices for solving the problem, given some approximation factor $c>1$ and the exact nearest neighbor $\boldsymbol{p}^{*}\in\mathcal{D}$ to $\boldsymbol{q}$ in $\mathcal{D}$. Here instead we will focus on methods for solving this problem exactly in sublinear time for average-case instances, by making suitable, natural assumptions on the distribution of the data points $\boldsymbol{p}_{i}$ and the query vector $\boldsymbol{q}$.

Definition 3.2** (Nearest neighbor problem on the sphere).**

Given a data set $\mathcal{D}$ of size $n$ drawn uniformly at random from $\mathcal{S}^{d-1}$, and a parameter $\theta\in(0,\frac{\pi}{2})$, the nearest neighbor problem on the sphere asks to index $\mathcal{D}$ in a data structure such that, when later given a uniformly random query $\boldsymbol{q}\in\mathcal{S}^{d-1}$, we can quickly find a neighbor $\boldsymbol{p}^{*}\in\mathcal{D}$ with $\phi(\boldsymbol{p}^{*},\boldsymbol{q})\leq\theta$, or conclude that with high probability no such vector $\boldsymbol{p}^{*}\in\mathcal{D}$ exists.

This description covers many applications where the data set and the query are indeed (close to) uniformly random, and we wish to detect highly similar objects in the database, if they exist. Note that there is often a one-to-one correspondence between such average-case, exact solutions and worst-case, approximate solutions – we can often substitute $\mathcal{D}$ with a (worst-case) data set $\mathcal{D}^{\prime}$ with the guarantee that no vectors in $\mathcal{D}^{\prime}\setminus\{\boldsymbol{p}^{*}\}$ have angle smaller than $\theta^{\prime}>\theta$ with the query vector, and use the same algorithms in both cases.

3.2 Locality-sensitive hashing

To solve the nearest neighbor problem on the unit sphere, we will follow the celebrated locality-sensitive hashing framework, first introduced in [28]. Let us first restate the folklore definition of locality-sensitive hash function, slightly adjusted to fit our model of exact nearest neighbor searching for the angular distance for random data points.

Definition 3.3** (Locality-sensitive hash functions).**

Let $c\in\mathbb{N}$. A set $\mathcal{H}$ of hash functions $h:\mathbb{R}^{d}\to\{1,\dots,c\}$ is called $(\theta,p_{1},p_{2})$-sensitive if (1) two vectors at an angle at most $\theta$ collide with probability at least $p_{1}$; and (2) two random vectors collide with probability at most $p_{2}$:

[TABLE]

With slight abuse of notation, we will sometimes use $p_{1}$ and $p_{2}$ to refer to the actual probabilities above (rather than their lower and upper bounds, respectively). Now, assuming such hash functions exist, and they are efficiently evaluable444Being able to compute $h(\boldsymbol{x})$ given $\boldsymbol{x}$ in $n^{o(1)}$ time may be sufficient to guarantee the stated asymptotic time complexities, but in practice a superpolynomial decoding cost $d^{\omega(1)}$ may already render the method impractical., we can build efficient nearest neighbor data structures as follows. After selecting suitable parameters $\ell$ and $t$, we construct $t$ hash tables $T_{1},\dots,T_{t}$, each consisting of $c^{\ell}$ hash buckets, corresponding to the concatenated outputs of $\ell$ randomly chosen hash functions $h_{1},\dots,h_{k}\sim\mathcal{H}$ from our hash family. For each vector $\boldsymbol{p}\in\mathcal{D}$ and each hash table $T_{i}$, we compute its combined hash value, and insert $\boldsymbol{p}$ in the corresponding hash bucket with label $(h_{i,1}(\boldsymbol{p}),\dots,h_{i,\ell}(\boldsymbol{p}))$. Then, given a query $\boldsymbol{q}$, for each hash table $T_{i}$ we compute its corresponding hash bucket $(h_{i,1}(\boldsymbol{q}),\dots,h_{i,\ell}(\boldsymbol{q}))$, look up vectors stored in this hash bucket, and check if any of them have angle less than $\theta$ with $\boldsymbol{q}$. After going through all the hash tables, we hope to have encountered the nearest neighbor $\boldsymbol{p}^{*}\in\mathcal{D}$ in one of these hash tables, while hopefully the total amount of work done at this point (computing hashes, doing hash table look-ups, going through potential nearest neighbors in the hash buckets) is still sublinear in $n$.

Besides the initialization costs of locality-sensitive hashing (e.g. setting up the hash functions, and initializing and populating the hash tables), the query costs can be divided in two main costs: (1) computing $t\cdot\ell$ hash values and doing $t$ lookups in memory to locate the right hash buckets; and (2) going through all potential nearest neighbors encountered in these hash buckets, to find a solution. By choosing $\ell$ and $t$ carefully, as a function of $n,\theta,p_{1},p_{2}$, we can balance these costs and minimize the query costs as stated in the folklore result below.

Proposition 3.4** (Locality-sensitive hashing for the nearest neighbor problem).**

*Let $\mathcal{H}$ be

an efficiently evaluable $(\theta,p_{1},p_{2})$-sensitive hash family, and let*

[TABLE]

Then, using locality-sensitive hashing with these parameters, we can solve the nearest neighbor problem with the following costs:

Initializing and populating the data structure can be done in time $\tilde{O}(n^{1+\rho})$; 2. 2.

The data structure requires $\tilde{O}(n^{1+\rho})$ memory to store; 3. 3.

With prob. $\geq 0.99$, within time $\tilde{O}(n^{\rho})$ the algorithm finds a neighbor at angle $\leq\theta$, if it exists.

Here $\tilde{O}$ hides multiplicative factors which are polynomial in $d$ and $\log n$, while the success probability can be fine-tuned through the leading constant of $t$; a larger $t$ leads to a higher success probability.

Note that the memory and query costs scale as $n^{1+\rho}$ and $n^{\rho}$ respectively, and so the goal is to find a locality-sensitive hash family $\mathcal{H}$ minimizing the associated parameter $\rho$, for the parameter $\theta$ of interest, preferably with an efficient decoding algorithm. Beyond this point, we will mostly ignore all other aspects of locality-sensitive hashing, and focus on finding families $\mathcal{H}$ minimizing $\rho$; in Section 10 we will come back to these practical aspects of locality-sensitive hashing.

3.3 Space partitions

Various locality-sensitive hash families achieving sublinear query times have been proposed over time. Although some are strictly more efficient than others, often there is a trade-off between e.g. the different performances for different angles $\theta$ (approximation factors $c$); the asymptotic efficiency for large $d$ and $n$; and the practical performance for small/moderate values $d$ and $n$. Examples of hash families for the angular distance from the literature include:

• •

The fast and practical hyperplane-based approach of Charikar [19];

• •

The asymptotically optimal cross-polytope hash family [47, 26, 48, 4, 30];

• •

The closely related (asymptotically optimal) simplex hash family [47, 48];

• •

The hypercube hash family [47, 48, 33];

• •

Hash families based on spherical caps [3, 5, 13, 6];

• •

Approaches based on using suitable integer lattices [3, 18].

More recently, it has been shown that deviating from the locality-sensitive hashing framework (by relaxing the condition that $h$ induces a proper partition of the space) may lead to better asymptotic results in certain applications [14, 13, 6, 22]. However, for data sets of subexponential size ($n=2^{o(d)}$), the best locality-sensitive hash methods achieve the same asymptotic performance as these locality-sensitive filters, and appear to perform better in practice [4, 40, 15, 11, 10]. Here we restrict our attention to locality-sensitive hash families generated by proper partitions of the unit sphere.

Out of the above methods, the hyperplane-based approach of Charikar [19] is often (one of) the most efficient in practice, due to its simple decoding algorithm which can be made extremely fast with heuristic tweaks [1]. For sufficiently large $d$ and $n$, the cross-polytope hash family will perform better, but e.g. in applications in cryptography, the cross-over point between hyperplane hashing and “optimal” hash families was found to be quite high [31, 14, 13, 37]. Rather than only focusing on asymptotic performance, we will search for hash families with similar properties to the hyperplane-based approach, ideally allowing us to interpolate between the fast hyperplane LSH and the asymptotically superior cross-polytope LSH.

3.4 Project–and–partition

A common approach for locality-sensitive hashing encompasses two steps:

Project.

Given vectors $\boldsymbol{x}\in\mathbb{R}^{d}$, we first project down onto a $k$-dimensional subspace through a random projection matrix $\boldsymbol{\mathord{\mathrm{A}}}\in\mathbb{R}^{k\times d}$, to obtain $\boldsymbol{\mathord{\mathrm{A}}}\boldsymbol{x}\in\mathbb{R}^{k}$. This projection is not lossless, and mutual distances/angles generally cannot be perfectly preserved through this transformation, but ideally these distances are preserved as best as possible. For this, here we generate $\boldsymbol{\mathord{\mathrm{A}}}$ by choosing each entry independently as a standard normal random variable. Note that with high probability, for constant $k=O(1)$ and large $d$ the $k$ rows of $\boldsymbol{\mathord{\mathrm{A}}}$ will be approximately orthogonal and have approximately the same norms [29]. Through the above definition of $\boldsymbol{\mathord{\mathrm{A}}}$ however, we eliminate the dependence on $d$ in the analysis of the collision probabilities $p_{1}$ and $p_{2}$ – the distribution of any particular column of $\boldsymbol{\mathord{\mathrm{A}}}$ is independent of $d$.

Partition.

Using a partition of $\mathbb{R}^{k}$ into $c$ regions, we assign a hash value to $\boldsymbol{x}$ according to the region $\boldsymbol{\mathord{\mathrm{A}}}\boldsymbol{x}$ landed in, using a decoding function $\mathrm{dec}:\mathbb{R}^{k}\to\{1,\dots,c\}$. A natural choice for the partition is to use the Voronoi cells induced by a suitable code $\mathcal{C}\subset\mathbb{R}^{k}$ of size $|\mathcal{C}|=c$. When all code words $\boldsymbol{c}\in\mathcal{C}$ have the same (unit) length, $\mathcal{C}$ forms a spherical code, and in the remainder we will exclusively focus on such codes $\mathcal{C}\subset\mathcal{S}^{k-1}$.

Given a spherical code $\mathcal{C}$, we thus obtain the following hash function family $\mathcal{H}$.

Definition 3.5** (Spherical code hash families).**

Let $\mathcal{C}\subset\mathcal{S}^{k-1}$ be a spherical code of size $|\mathcal{C}|=c$, and let $\operatorname{dec}:\mathbb{R}^{k}\to\{1,\dots,c\}$ be a decoding algorithm for the Voronoi cells of $\mathcal{C}$. Then we associate to $\mathcal{C}$ hash functions $h=h_{\boldsymbol{\mathord{\mathrm{A}}}}$ as follows:

[TABLE]

For the corresponding hash family $\mathcal{H}$, sampling $h\sim\mathcal{H}$ corresponds to sampling $\boldsymbol{\mathord{\mathrm{A}}}\sim\mathcal{N}(0,1)^{k\times d}$.

Observe that in all cases the decoding step can be done in time $O(ck)$, by computing all distances between $\boldsymbol{\mathord{\mathrm{A}}}\boldsymbol{x}$ and code words $\boldsymbol{c}\in\mathcal{C}$, but for some codes with additional structure, faster decoding algorithms may exist. Projections can naively be computed in time $O(dk)$, for a total hash time complexity of $O(dk+ck)$.

As an example, for hyperplane LSH [19] we have $k=1$ (i.e. we project onto a line), and the matrix $\boldsymbol{\mathord{\mathrm{A}}}\in\mathbb{R}^{1\times d}$ is a vector drawn from a spherically-symmetric distribution (e.g. a Gaussian distribution). The partition is then defined by the Voronoi cells of the one-dimensional spherical code $\mathcal{C}=\{-1,1\}$; the sign of the inner product $\boldsymbol{\mathord{\mathrm{A}}}\boldsymbol{x}\in\mathbb{R}$ determines whether the hash value equals [math] or $1$.

3.5 Orthant probabilities

Without loss of generality, let $\boldsymbol{x}=(1,0,\dots,0)$ and, in case of a predetermined angle $\theta$ with $\boldsymbol{x}$, let $\boldsymbol{y}=(\cos\theta,\sin\theta,0,\dots,0)$. Denoting the $k$-dimensional columns of $\boldsymbol{\mathord{\mathrm{A}}}$ by $\boldsymbol{a}_{i}$, observe that $\boldsymbol{\mathord{\mathrm{A}}}\boldsymbol{x}=\boldsymbol{a}_{1}$ and $\boldsymbol{\mathord{\mathrm{A}}}\boldsymbol{y}=(\cos\theta)\boldsymbol{a}_{1}+(\sin\theta)\boldsymbol{a}_{2}$. Note that while $\boldsymbol{x}$ and $\boldsymbol{y}$ have angle exactly $\theta$, $\boldsymbol{\mathord{\mathrm{A}}}\boldsymbol{x}$ and $\boldsymbol{\mathord{\mathrm{A}}}\boldsymbol{y}$ only have angle approximately $\theta$, depending on how orthogonal the two random vectors $\boldsymbol{a}_{1},\boldsymbol{a}_{2}\sim\mathcal{N}(0,1)^{k}$ are. This means that we can compute $p_{1}$ as follows:

[TABLE]

For two independently sampled vectors $\boldsymbol{x},\boldsymbol{y}\sim\mathcal{S}^{d-1}$, also $\boldsymbol{\mathord{\mathrm{A}}}\boldsymbol{x}$ and $\boldsymbol{\mathord{\mathrm{A}}}\boldsymbol{y}$ are independent. Given any spherical code $\mathcal{C}=\{\boldsymbol{c}_{1},\dots,\boldsymbol{c}_{c}\}\in\mathcal{S}^{k-1}$, we can thus describe the probability $p_{2}$ as follows:

[TABLE]

Here $\mathcal{V}_{i}\subset\mathbb{R}^{k}$ denotes the Voronoi region associated to $\boldsymbol{c}_{i}$, and $\mu(\mathcal{V}_{i})$ its measure. Note that for a size-$c$ partition of $\mathbb{R}^{k}$, $p_{2}$ is smallest if all regions have the same size, in which case $p_{2}=1/c$.

To further study $p_{1}$, for general codes $\mathcal{C}$ by definition we have $\operatorname{dec}(\boldsymbol{x})=\boldsymbol{c}_{i}$ if and only if $\|\boldsymbol{x}-\boldsymbol{c}_{i}\|\leq\|\boldsymbol{x}-\boldsymbol{c}_{j}\|$ for all $j=1,\dots,c$.555For simplicity we ignore cases where $\boldsymbol{x}$ lies on the boundary between two Voronoi cells and decoding is not unique, as the boundaries together have measure [math]. Since we assume that $\mathcal{C}$ is a spherical code, this is equivalent to $\langle\boldsymbol{x},\boldsymbol{c}_{i}\rangle\geq\langle\boldsymbol{x},\boldsymbol{c}_{j}\rangle$ for all $j=1,\dots,c$, or equivalently $\langle\boldsymbol{x},\boldsymbol{c}_{i}-\boldsymbol{c}_{j}\rangle\geq 0$ for all $j=1,\dots,c$. The event that $\boldsymbol{x}\in\mathbb{R}^{k}$ is mapped to the hash region corresponding to $\boldsymbol{c}_{i}$ (i.e. $\boldsymbol{x}\in\mathcal{V}_{i}$) can thus equivalently be described by the following matrix inequality:

[TABLE]

Let $\boldsymbol{\mathord{\mathrm{C}}}_{i}$ denote the matrix on the left hand side above, with rows $\boldsymbol{c}_{i}-\boldsymbol{c}_{j}$ for $j=1,\dots,c$. Recall that for the event of two $\theta$-correlated random vectors to both be mapped to $\boldsymbol{c}_{i}$, we need $\boldsymbol{\mathord{\mathrm{C}}}_{i}\boldsymbol{x}\geq\boldsymbol{0}$ to hold for both $\boldsymbol{x}=\boldsymbol{a}_{1}$ and $\boldsymbol{x}=(\cos\theta)\boldsymbol{a}_{1}+(\sin\theta)\boldsymbol{a}_{2}$, where $\boldsymbol{a}_{1},\boldsymbol{a}_{2}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\mathord{\mathrm{I}}}_{k})$. The probability $p_{1}$ can thus be described by a $2k$-dimensional system of inequalities:

[TABLE]

Here $\otimes$ denotes the standard Kronecker product for matrices. Note that if $\boldsymbol{z}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\mathord{\mathrm{I}}}_{2k})$, then $\boldsymbol{z}^{\prime}=\boldsymbol{\mathord{\mathrm{B}}}\boldsymbol{z}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\mathord{\mathrm{B}}}\boldsymbol{\mathord{\mathrm{B}}}^{T})$. Furthermore, using $(\boldsymbol{\mathord{\mathrm{A}}}\otimes\boldsymbol{\mathord{\mathrm{B}}})\cdot(\boldsymbol{\mathord{\mathrm{A}}}\otimes\boldsymbol{\mathord{\mathrm{B}}})^{T}=(\boldsymbol{\mathord{\mathrm{A}}}\boldsymbol{\mathord{\mathrm{A}}}^{T})\otimes(\boldsymbol{\mathord{\mathrm{B}}}\boldsymbol{\mathord{\mathrm{B}}}^{T})$ for arbitrary $\boldsymbol{\mathord{\mathrm{A}}}$ and $\boldsymbol{\mathord{\mathrm{B}}}$, we obtain:

[TABLE]

The probabilities above are also known as orthant probabilities for the multivariate normal distribution with covariance matrix $\boldsymbol{\mathord{\mathrm{\Sigma}}}_{i}$. For completeness, note that $p_{2}$ can equivalently be described as an orthant probability, as:

[TABLE]

This description is computationally more difficult to work with than the approach using volumes of the corresponding Voronoi cells, as in Equation (14). To summarize, the above derivation can be seen as a proof of the following theorem. Here $\operatorname{dec}(\cdot)=\operatorname{dec}_{\mathcal{C}}(\cdot)$ denotes a decoding function for the corresponding spherical code, in terms of its Voronoi cells $\mathcal{V}_{i}$: $\operatorname{dec}(\boldsymbol{x})=i$ if $\boldsymbol{x}\in\mathcal{V}_{i}$ lies closer to $\boldsymbol{c}_{i}$ than to all other code words $\boldsymbol{c}_{j}\in\mathcal{C}$.

Theorem 3.6** (Spherical code locality-sensitive hashing).**

Let $\mathcal{C}=\{\boldsymbol{c}_{1},\dots,\boldsymbol{c}_{c}\}\subset\mathcal{S}^{k-1}$ be a $k$-dimensional spherical code, and for $i=1,\dots,c$, let $\boldsymbol{\mathord{\mathrm{C}}}_{i}\in\mathbb{R}^{c\times k}$ denote the matrix whose $j^{\text{th}}$ row is the vector $\boldsymbol{c}_{i}-\boldsymbol{c}_{j}$. Let $\mathcal{H}$ consist of functions $h=h_{\boldsymbol{\mathord{\mathrm{A}}}}$ of the form:

[TABLE]

Then for any $\theta\in(0,\frac{\pi}{2})$, $\mathcal{H}$ is a $(\theta,p_{1},p_{2})$-sensitive hash family, with:

[TABLE]

Note again that, by our choice of $\boldsymbol{\mathord{\mathrm{A}}}$, the parameters $p_{1},p_{2}$ are independent of $d$. For large $d\gg k$, choosing $\boldsymbol{\mathord{\mathrm{A}}}\sim\mathcal{N}(0,1)^{k\times d}$ or choosing $\boldsymbol{\mathord{\mathrm{A}}}$ to have orthogonal, unit-norm rows is essentially equivalent [23, 29]. For smaller $d\approx k$ however, these definitions are not the same, and it may be worthwhile to use orthogonal rows instead, to decrease the distortion of mutual distances by the projection $\boldsymbol{\mathord{\mathrm{A}}}$. For instance, when using dense $d$-dimensional spherical codes for partitions, it may be impossible for two almost-orthogonal vectors to be in the same hash region, and $p_{2}$ may quickly approach [math] while $p_{1}$ remains large. Although our study may lead to advances and insights in the regime of $k\approx d$ as well, the main topic of interest here however is to first project down onto a $k$-dimensional space with $k\ll d$, and then use a suitable spherical code there for partitions.

3.6 Special spherical codes

For special classes of spherical codes, the expressions of Theorem 3.6 can be further simplified. We state two specific cases of interest below. First, if $\mathcal{C}$ is such that the Voronoi cells of the vertices each cover an equal part of the sphere, i.e. $\mu(\mathcal{V}_{i})=\mu(\mathcal{V}_{j})$ for all $i$ and $j$, then we get the following slightly simplified corollary. Note that given a fixed code size $c$, the value $p_{2}$ is minimized exactly when all Voronoi cells are of equal size. As we want to minimize $\rho=\rho(p_{1},p_{2})$ which is an increasing function of $p_{2}$, this suggests that using such codes may well lead to the best results.666The stated intuition does not provide a formal proof that $\rho$ is minimized when all Voronoi cells have the same volume, as there is also an intricate dependence of $p_{1}$ (and $\rho$) on the shape of the Voronoi cells. If all Voronoi cells are of equal size, but their shapes are “bad”, then the resulting code may not be useful in our framework.

Corollary 3.7** (Collision probabilities for uniform spherical codes).**

Let $\mathcal{C}=\{\boldsymbol{c}_{1},\dots,\boldsymbol{c}_{c}\}\subset\mathcal{S}^{k-1}$ be a $k$-dimensional uniform spherical code, in the sense that $\mu(\mathcal{V}_{i})=\mu(\mathcal{V}_{j})$ for all $i,j$. Let $\boldsymbol{\mathord{\mathrm{C}}}_{i}$ and $\mathcal{H}$ be as in Theorem 3.6. Then for any $\theta\in(0,\frac{\pi}{2})$, $\mathcal{H}$ is a $(\theta,p_{1},p_{2})$-sensitive hash family, with

[TABLE]

In case $\mathcal{C}$ has many symmetries, and in particular is isogonal (vertex-transitive), then Theorem 3.6 can be further simplified, as in that case all hash regions are equivalent, and only one of the orthant probabilities needs to be evaluated to determine $p_{1},p_{2}$ and $\rho$. Note that by definition, isogonal codes are a subset of the uniform spherical codes above.

Corollary 3.8** (Collision probabilities for isogonal spherical codes).**

Let $\mathcal{C}=\{\boldsymbol{c}_{1},\dots,\boldsymbol{c}_{c}\}\subset\mathcal{S}^{k-1}$ be a $k$-dimensional isogonal spherical code, in the sense that for each $\boldsymbol{c}_{i},\boldsymbol{c}_{j}\in\mathcal{C}$, there exists an isometry $f$ with $f(\boldsymbol{c}_{i})=\boldsymbol{c}_{j}$ and $f(\mathcal{C})=\mathcal{C}$. Let $\boldsymbol{\mathord{\mathrm{C}}}_{i}$ and $\mathcal{H}$ be as in Theorem 3.6. Then for any $\theta\in(0,\frac{\pi}{2})$, $\mathcal{H}$ is a $(\theta,p_{1},p_{2})$-sensitive hash family, with:

[TABLE]

For isogonal codes, computing the query exponent $\rho$ comes down to computing a single orthant probability.

3.7 Relevant vectors

Finally, the computations of the collision probabilities $p_{1}$ and $p_{2}$ for a given spherical code can sometimes be simplified by reducing the number of conditions that define when a point is decoded to a given code word. Recall that the matrix $\boldsymbol{\mathord{\mathrm{C}}}_{i}$, containing difference vectors $\boldsymbol{c}_{i}-\boldsymbol{c}_{j}$ as its rows, originates from the equivalence below:

[TABLE]

In many cases however, we can already conclude that $\operatorname{dec}(\boldsymbol{x})=i$ by only computing some of the distances $\|\boldsymbol{x}-\boldsymbol{c}_{j}\|$ – if $\boldsymbol{x}$ is closer to $\boldsymbol{c}_{i}$ than to any of the closest neighbors of $\boldsymbol{c}_{i}$ in $\mathcal{C}$, then it must be closest to $\boldsymbol{c}_{i}$. More generally, we may define $\mathcal{R}_{i}\subseteq\mathcal{C}$ as the set of relevant vectors for $\boldsymbol{c}_{i}$, i.e. those vectors $\boldsymbol{c}_{j}\in\mathcal{C}$ whose Voronoi cells $\mathcal{V}_{j}$ share a non-trivial boundary with $\mathcal{V}_{i}$. The conditions to verify that decoding $\boldsymbol{x}$ results in $i$ may then be simplified:

[TABLE]

Computationally, for evaluating the orthant probabilities in $p_{1}$ in Theorem 3.6, this allows us to replace the matrix $\boldsymbol{\mathord{\mathrm{C}}}_{i}$ (with $c$ rows $\boldsymbol{c}_{i}-\boldsymbol{c}_{j}$ for $\boldsymbol{c}_{j}\in\mathcal{C}$), with a smaller matrix $\hat{\boldsymbol{\mathord{\mathrm{C}}}}_{i}$ (consisting of $|\mathcal{R}_{i}|$ rows $\boldsymbol{c}_{i}-\boldsymbol{c}_{j}$ for $\boldsymbol{c}_{j}\in\mathcal{R}_{i}$). This can significantly reduce the dimensionality of the resulting orthant probability problem, and this may assist in evaluating $p_{1}$ more accurately via numerical integration when no closed form expression is available.

Proposition 3.9** (Computing collision probabilities using relevant vectors).**

*Theorem 3.6

and Corollaries 3.7 and 3.8 remain valid when we replace the matrices $\boldsymbol{\mathord{\mathrm{C}}}_{i}\in\mathbb{R}^{c\times k}$ by $\hat{\boldsymbol{\mathord{\mathrm{C}}}}_{i}\in\mathbb{R}^{|\mathcal{R}_{i}|\times k}$, consisting of the rows $\boldsymbol{c}_{i}-\boldsymbol{c}_{j}$ with $\boldsymbol{c}_{j}\in\mathcal{R}_{i}$.*

With this general framework in hand, we can now compute collision probabilities for arbitrary spherical codes, and compute the resulting nearest neighbor exponents $\rho$, given $\mathcal{C}$ and $\theta$.

4 One-dimensional codes

Let us first illustrate the above framework by rederiving the hyperplane locality-sensitive hashing results of Charikar [19]. For the $1$-dimensional unit sphere $\mathcal{S}^{0}=\{(-1),(1)\}$, there is only one choices for a spherical code $\mathcal{C}\subseteq\mathcal{S}^{0}$ of size at least two, namely $\mathcal{C}=\{(-1),(1)\}$. This code is an isogonal code, and $\mathcal{R}_{(1)}=\{(-1)\}$ and $\mathcal{R}_{(-1)}=\{(1)\}$ are the relevant vectors. As in Proposition 3.9, $\hat{\boldsymbol{\mathord{\mathrm{C}}}}_{1}$ consists of one row, containing the vector $(1)-(-1)=(2)$. By Corollary 3.8 and Proposition 3.9, we can compute the collision probabilities, after projecting down using a random vector $\boldsymbol{\mathord{\mathrm{A}}}\sim\mathcal{N}(0,1)^{1\times d}$, as:

[TABLE]

For $p_{1}$, observe that $\Pr_{\boldsymbol{z}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\mathord{\mathrm{\Sigma}}}_{1})}\left(\boldsymbol{z}\geq\boldsymbol{0}\right)=\Pr_{\boldsymbol{z}\sim\mathcal{N}(\boldsymbol{0},\mu\cdot\boldsymbol{\mathord{\mathrm{\Sigma}}}_{1})}\left(\boldsymbol{z}\geq\boldsymbol{0}\right)$ for all scalars $\mu$. We can thus eliminate the factor $4$, and we are left with a standard orthant probability for the bivariate normal distribution with covariance $\cos\theta$. The derivation of a closed-form formula for the orthant probability of the bivariate normal distribution is often attributed to Sheppard in the late 1800s [44], who showed that the orthant probability for the bivariate normal distribution with Pearson correlation coefficient $a$ is:

[TABLE]

For our application, this implies the well-known result of Charikar [19], stated below.

Proposition 4.1** (Collision probabilities for antipodal codes).**

Let $\mathcal{C}$ consist of two antipodal points on the unit sphere. Then the corresponding project–and–partition hash family $\mathcal{H}$ is $(\theta,p_{1},p_{2})$-sensitive, with:

[TABLE]

For the nearest neighbor query exponent $\rho=\log p_{1}/\log p_{2}$, a simple computation leads to e.g. $\rho=\log_{2}(3/2)\approx 0.5850$ for $\theta=\frac{\pi}{3}$. Table 1 lists several query exponents $\rho$ for this hash family, for different parameters $\theta\in\frac{\pi}{12}\{1,2,3,4,5\}$. Note that for $\theta\to 0$ we have $p_{1}\to 1$ and $\rho\to 0$, while for $\theta\to\frac{\pi}{2}$ we have $p_{1}\to p_{2}=\frac{1}{2}$ and $\rho\to 1$, as expected.

5 Two-dimensional codes

The simplest previously unstudied case is to project down onto a two-dimensional subspace through a random projection matrix $\boldsymbol{\mathord{\mathrm{A}}}\sim\mathcal{N}(0,1)^{2\times d}$, and to use a suitable spherical code on the two-dimensional sphere (circle) to divide these projected vectors into hash buckets. Using two (antipodal) points on the circle would be equivalent to using the one-dimensional hyperplane hashing described above, so the non-trivial cases start at spherical codes of size at least three.

5.1 Two-dimensional isogonal spherical codes

The following theorem summarizes the collision probabilities for what is arguably the most natural choice for a spherical code on the circle, which is to use any number of equidistributed points on the circle (i.e. to use the vertices of a regular polygon to partition the sphere).

Theorem 5.1** (Collision probabilities for regular polygons).**

Let $\mathcal{C}\subset\mathcal{S}^{1}$ consist of the vertices of the regular $c$-gon, for $c\geq 2$, and let $\theta\in(0,\frac{\pi}{2})$. Then the corresponding project–and–partition hash family $\mathcal{H}$ is $(\theta,p_{1},p_{2})$-sensitive, with:

[TABLE]

To prove the above result, we will first derive the shape of the correlation matrix below.

Lemma 5.2** (The correlation matrix for regular polygons).**

Let $\mathcal{C}\subset\mathcal{S}^{1}$ consist of the vertices of the regular $c$-gon, for $c\geq 2$. Then the correlation matrix $\boldsymbol{\mathord{\mathrm{\Sigma}}}_{1}$ in Theorem 3.8 (Prop. 3.9) satisfies:

[TABLE]

Proof 5.3**.**

Without loss of generality777Note that although the matrices $\boldsymbol{\mathord{\mathrm{C}}}_{i}$ depend on the absolute positioning of the code words, the dot products $\boldsymbol{\mathord{\mathrm{C}}}_{i}\boldsymbol{\mathord{\mathrm{C}}}_{i}^{T}$ only depends on the relative positions between code words, and is invariant under rotations of $\mathcal{C}$., let us fix one vertex of $\mathcal{C}$ at $(1,0)$. For ease of notation, let us write $c_{j}=\cos\frac{2j\pi}{c}$ and $s_{j}=\sin\frac{2j\pi}{c}$, so that $\mathcal{C}=\{(c_{j},s_{j}):j=1,\dots,c\}$. Now, a point in $\mathbb{R}^{2}$ is closest to $\boldsymbol{c}_{c}=(1,0)$ out of all code words in $\mathcal{C}$ iff it lies closer to $(1,0)$ than to the relevant vectors of $(1,0)$, which are its left and right neighbors on the circle, $\mathcal{R}_{(1,0)}=\{(c_{1},s_{1}),(c_{c-1},s_{c-1})\}=\{(c_{1},s_{1}),(c_{1},-s_{1})\}$. Let $\hat{\boldsymbol{\mathord{\mathrm{C}}}}_{i}$ as in Proposition 3.9 denote the matrix $\boldsymbol{\mathord{\mathrm{C}}}_{i}$ reduced to its relevant rows. For the matrix $\hat{\boldsymbol{\mathord{\mathrm{D}}}}_{c}=\hat{\boldsymbol{\mathord{\mathrm{C}}}}_{c}\hat{\boldsymbol{\mathord{\mathrm{C}}}}_{c}^{T}\in\mathbb{R}^{2\times 2}$, we therefore obtain:

[TABLE]

Due to the spherical/circular symmetry of $\mathcal{C}$, all matrices $\hat{\boldsymbol{\mathord{\mathrm{D}}}}_{i}=\hat{\boldsymbol{\mathord{\mathrm{C}}}}_{i}\hat{\boldsymbol{\mathord{\mathrm{C}}}}_{i}^{T}$ are equal to $\hat{\boldsymbol{\mathord{\mathrm{D}}}}_{c}$. Together with Theorem 3.8, this leads to the given expression for $\boldsymbol{\mathord{\mathrm{\Sigma}}}=\boldsymbol{\mathord{\mathrm{\Sigma}}}_{1}$. ∎

Orthant probabilities for general quadrivariate normal distributions are difficult to solve, but for certain specific forms of $\boldsymbol{\mathord{\mathrm{\Sigma}}}$, explicit formulas for the resulting orthant probabilities in terms of the off-diagonal entries of $\boldsymbol{\mathord{\mathrm{\Sigma}}}$ are known. In the 1960s, Cheng [20, 21] used a clever path integration technique, together with Plackett’s reduction formula [39], to ultimately end up with the following expressions for the orthant probabilities for matrices of the above form.

Lemma 5.4** (Orthant probabilities for a class of quadrivariate normal distributions).**

[20, Appendix]** Let $a,b\in\mathbb{R}$, and let the correlation matrix $\boldsymbol{\mathord{\mathrm{\Sigma}}}=\boldsymbol{\mathord{\mathrm{\Sigma}}}_{a,b}$ be of the following form:

[TABLE]

Then the orthant probability $\Phi(a,b)=\Pr_{\boldsymbol{z}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\mathord{\mathrm{\Sigma}}})}\left(\boldsymbol{z}\geq\boldsymbol{0}\right)$ for the correlation matrix $\boldsymbol{\mathord{\mathrm{\Sigma}}}$ satisfies:

[TABLE]

The orthant probability $\Phi(a,b)$ above can be equivalently written in the following form:

[TABLE]

Note that $a,b,ab$ are the off-diagonal entries of $\boldsymbol{\mathord{\mathrm{\Sigma}}}$, or equivalently, the arguments of the arccosines $-a,-b,ab$ are the off-diagonal entries of $\boldsymbol{\mathord{\mathrm{\Sigma}}}^{-1}$, thus showing similarities with the orthant probabilities for $2\times 2$ matrices related to antipodal spherical codes, where $\Phi(a)=\arccos(-a)/(2\pi)$. An interesting open problem, as we will see later, is to establish a pattern in the orthant probabilities for the matrices $\boldsymbol{\mathord{\mathrm{\Sigma}}}$ appearing in Theorem 3.6 and Corollary 3.8.

Now, with the above lemma in hand, we are ready to prove Theorem 5.1.

Proof 5.5** (Proof of Theorem 5.1).**

From Theorem 3.8, and the observation that regular polygons form isogonal codes on the sphere, we obtain the correlation matrix $\boldsymbol{\mathord{\mathrm{\Sigma}}}$ as described by Lemma 5.2. Recall that orthant probabilities are invariant under scalar multiplication of the random vector $\boldsymbol{z}$; removing the leading scalar factor of $\boldsymbol{\mathord{\mathrm{\Sigma}}}$ in Lemma 5.2, the collision probability $p_{1}$ remains the same. After removing this factor, we obtain a matrix $\boldsymbol{\mathord{\mathrm{\Sigma}}}$ of the form of Lemma 5.4, where $a=\cos\theta$ and $b=-\cos\frac{2\pi}{c}$. Through some elementary trigonometric reductions, we then obtain the stated result, where we obtain an additional leading factor $c$ as $\Pr(h(\boldsymbol{x})=h(\boldsymbol{y}))=c\cdot\Pr(h(\boldsymbol{x})=h(\boldsymbol{y})=1)$. ∎

Let us draw some observations from Theorem 5.1. First, observe that if we substitute $c=2$, i.e. using a regular $2$-gon, we obtain $p_{1}=1-\theta/\pi$ and $p_{2}=1/2$, as we expect from Proposition 4.1. Similarly, substituting $c=4$ into Theorem 5.1 we obtain $p_{1}=(1-\theta/\pi)^{2}$ and $p_{2}=(1/2)^{2}$, matching the combined collision probability of two independent antipodal codes. As the nearest neighbor exponents $\rho=\log p_{1}/\log p_{2}=\log p_{1}^{2}/\log p_{2}^{2}$ are invariant under applying the same exponentiation to both $p_{1}$ and $p_{2}$, both the regular $2$-gon and $4$-gon (square) achieve the same performance when using Gaussian projection matrices $\boldsymbol{\mathord{\mathrm{A}}}$.

For other values of $c\geq 2$, the collision probabilities do not simplify much further than the given expression in Theorem 5.1 with $c$ replaced by the corresponding number. Looking at the nearest neighbor exponents $\rho=\log p_{1}/\log p_{2}$ for e.g. $\theta=\frac{\pi}{3}$ and $\theta=\frac{\pi}{4}$, we obtain the two graphs in Figure 5. As we can see, the query exponents for $c=2$ and $c=4$ are the same, and using very large polygons does not improve the performance; any polygon with more than $4$ vertices achieves worse query exponents $\rho$ than antipodal locality-sensitive hashing. For $c=3$ however we observe a non-trivial minimum: using equilateral triangles, the query time complexity for nearest neighbor searching is lower than for $c=2,4$.

Theorem 5.6** (In two dimensions, the triangle is optimal).**

Let $\theta\in(0,\frac{\pi}{2})$. Then, among all project–and–partition hash families induced by two-dimensional isogonal spherical codes, the query exponent $\rho=\log p_{1}/\log p_{2}$ is minimized by the equilateral triangle.

Proof 5.7**.**

By inspecting the function $f(c,\theta)=\log p_{1}/\log p_{2}$ with $p_{1}$ and $p_{2}$ as in (1), we see that on the positive reals, its derivative with respect to $c$ always has a single root $c_{0}=c_{0}(\theta)\in(2,3)$. For $c\leq c_{0}$ the derivative is negative and for $c\geq c_{0}$ the derivative is positive, hence on the positive reals, $f(c,\theta)$ is minimized at $c=c_{0}$, where the exact value $c_{0}$ depends on $\theta$. On the natural numbers, $f(c,\theta)$ is therefore minimized either at $c=2$ or at $c=3$. Since $f(2,\theta)=f(4,\theta)>f(3,\theta)$, the value at $c=3$ is the unique global minimum, and therefore for all target nearest neighbor angles $\theta\in(0,\frac{\pi}{2})$, the equilateral triangle achieves the lowest values $\rho$. ∎

For example, when $\theta=\frac{\pi}{3}$ we obtain a query exponent $\rho=0.56996\dots$ for the equilateral triangle, compared to $\rho=0.58496\dots$ for the square or for two antipodal points. Recall again that this improvement holds unconditionally, for all $d$; these collision probabilities are exact, and do not hide order terms in $d$. A further overview of query exponents $\rho$ for regular polygons is given in Table 1, in the rows labeled $k=2$.

5.2 Arbitrary two-dimensional spherical codes

For arbitrary (non-isogonal) two-dimensional spherical codes, where code word $\boldsymbol{c}_{i}$ covers a fraction $f_{i}$ of the circle, the probability that two vectors collide in one of the $c$ hash regions, given their mutual angle $\theta$, can be computed in terms of the vector $\boldsymbol{f}$ as:

[TABLE]

The set of query exponents $\rho$ that can be achieved by two-dimensional codes $\mathcal{C}\subset\mathcal{S}^{1}$ in the project–and–partition framework, given a target angle $\theta$, is then equal to:

[TABLE]

Here $\|\boldsymbol{f}\|_{2}$ and $\|\boldsymbol{f}\|_{1}$ denote the $\ell_{2}$-norm and $\ell_{1}$-norm of $\boldsymbol{f}$. Although a numerical evaluation indicates that, for arbitrary $\theta\in(0,\frac{\pi}{2})$, the minimum over $R_{\theta}$ is obtained for $\boldsymbol{f}=(\frac{1}{3},\frac{1}{3},\frac{1}{3})$, we leave a formal proof of the optimality of the equilateral triangle over all two-dimensional codes as an open problem.

6 Three-dimensional codes

As the dimensionality increases, the dimensions of the matrices $\boldsymbol{\mathord{\mathrm{C}}}_{i}$ in Theorem 3.6 and Corollary 3.8 generally also increase, even when using the relevant vectors optimization of Proposition 3.9. In particular, as each Voronoi region in three dimensions is defined by at least three relevant vectors, the matrix $\boldsymbol{\mathord{\mathrm{\Sigma}}}$ will be at least six-dimensional with non-trivial off-diagonal entries. The bad news is that for none of the three-dimensional codes we considered, we could derive analytic, closed-form expressions for the collision probabilities; numerical integration seems to be necessary to evaluate their performance. The good news however is that Section 3 provides us with a computational framework to evaluate the suitability of a spherical code for nearest neighbor searching, given the vertices defining the spherical code.

Table 1 summarizes the results of a brief survey of various spherical codes appearing in the literature in different contexts. The three-dimensional codes are given in the rows labeled $k=3$. The codes we considered and evaluated numerically are:

• •

The tetrahedron, or the $3$-simplex, consisting of $4$ equidistributed points on the sphere;

• •

The octahedron, or the $3$-orthoplex (cross-polytope), consisting of the $6$ vectors $\pm\boldsymbol{e}_{i}$;

• •

The cube, consisting of $8$ vertices of the form $\frac{1}{3}\sqrt{3}\cdot(\pm 1,\pm 1,\pm 1)$;

• •

The icosahedron, one of the five Platonic solids, with $12$ vertices;

• •

The cuboctahedron, the vertex figure of the root lattices $A_{3}$ and $D_{3}$, with $12$ vertices;

• •

The dodecahedron, the largest of the Platonic solids, with $20$ vertices;

• •

Three sphere packings of $5,7,9$ vertices from the sphere packing database of Sloane [45].

Numerically evaluating the resulting query exponents $\rho$, it seems that for all $\theta\in(0,\frac{\pi}{2})$, the tetrahedron gives the best query exponents $\rho$. Similar to the two-dimensional case, using large, dense spherical codes does not seem to help; a fine-grained partition may make $p_{2}$ smaller, but at the same time $p_{1}$ then also becomes so much smaller that $\rho$ still increases. Note that, similar to the triangle, the tetrahedron belongs to the class of simplices.

Conjecture 6.1** (In three dimensions, the tetrahedron is optimal).**

*Let $\theta\in(0,\frac{\pi}{2})$. Then,

among all project–and–partition hash families induced by three-dimensional (isogonal) spherical codes, the query exponent $\rho=\log p_{1}/\log p_{2}$ is minimized by the regular tetrahedron.*

For the tetrahedron, the collision probabilities are $p_{2}=1/4$ and:

[TABLE]

A further simplification of $p_{1}$, i.e. finding a closed-form expression for this six-dimensional orthant probability is left as an open problem [34].

For numerically computing the values $\rho$ in Table 1 for three-dimensional codes, we used the statistical software R together with the mvtnorm-package [27], which specifically implements fast and precise evaluation of orthant probabilities. For instance, for the tetrahedron with $\theta=\frac{\pi}{3}$, we numerically compute $\rho=0.559998$. Using an independent C implementation, we performed Monte Carlo experiments, which for $\theta=\frac{\pi}{3}$ led to the estimate $\rho=0.559997$ based on $10^{8}$ random trials.

Note that the (numerical) integration task becomes more complex, and likely less precise, as the dimensionality of the problem increases (i.e. the size of $\Sigma$ increases). For isogonal codes, the dimensions of $\Sigma$ are not determined by the size of the code but by the number of relevant vectors for each code word. The entries in Table 1 may therefore be less precise for e.g. the icosahedron, with $|\mathcal{R}|=5$, and for large irregular sphere packings from [45].

7 Four-dimensional codes

As for the three-dimensional case, the orthant probabilities in the four-dimensional case generally appear too hard to evaluate analytically. Instead, we numerically computed parameters $\rho$ for different angles $\theta$ for the following list of four-dimensional spherical codes, by estimating the corresponding multivariate integrals with the R-library [27]:

• •

The pentatope or $\boldsymbol{5}$-cell, also known as the $4$-simplex, with $5$ vertices;

• •

The $\boldsymbol{16}$-cell, the four-dimensional version of the orthoplex, with $8$ vertices;

• •

The tesseract, or the four-dimensional hypercube, with $16$ vertices;

• •

The octacube or $\boldsymbol{24}$-cell, the vertex figure of the $D_{4}$ root lattice, with $24$ vertices;

• •

The runcinated $5$-cell, corresponding to the $20$ roots of the $A_{4}$-lattice;

• •

Four sphere packings of $6,7,10,13$ vertices from the sphere packing database of Sloane [45].

The only other regular convex polytopes not mentioned above are the $120$-cell and $600$-cell which, as their names suggest, are rather large; not only will the numerical evaluation (or Monte Carlo estimation) of the corresponding orthant probabilities likely be imprecise, but we also expect that such large codes are unlikely to offer an improvement, if we extrapolate the intuition obtained from the results in two and three dimensions to the four-dimensional setting.

As shown in Table 1, an interesting phenomenon occurs in dimension $4$. For two- and three-dimensional spherical codes we saw that the simplices give the best performance, beating other codes such as the orthoplices. In four dimensions however, for small angles $\theta$, the more fine-grained partitioning offered by the orthoplex gives a better performance for the associated locality-sensitive hashing scheme. Note that the best parameters $\rho$ are achieved by the simplex and orthoplex; these polytopes seem well-suited for locality-sensitive hashing, and even the nicest other spherical codes cannot compete with these highly regular and symmetric shapes.

Conjecture 7.1** (In four dimensions, the $5$-cell and $16$-cell are optimal).**

*Let $\theta\in(0,\frac{\pi}{2})$.

There exists a parameter $\theta_{0}\in(0,\frac{\pi}{2})$ such that, among all project–and–partition hash families induced by four-dimensional (isogonal) spherical codes, the query exponent $\rho$ is minimized by:*

• •

the $16$-cell, if $\theta\in(0,\theta_{0})$;

• •

the $5$-cell, if $\theta\in(\theta_{0},\frac{\pi}{2})$.

Numerically, we estimate the cross-over point between the $5$-cell and $16$-cell in terms of the query exponent $\rho$ to lie at $\theta_{0}\approx 0.33\pi$, slightly below $\theta=\frac{\pi}{3}$. See also Figure 1, showing the improvements over hyperplane hashing for both the $5$-cell and the $16$-cell, and showing that the cross-over point between both four-dimensional codes lies close to $\theta=\frac{\pi}{3}$. From a practical point of view however, as discussed in Section 10, using larger spherical codes (in a fixed dimension) is generally better, as the hash computations are then cheaper. In that sense, using the $16$-cell may be a better choice than using a $5$-cell for partitions even for angles $\theta>\theta_{0}$.

For completeness, for the $5$-cell and $16$-cell we performed independent Monte Carlo estimations of the collision probabilities using $10^{8}$ trials for each value $\theta$ in Table 1, to verify that these numbers are accurate. The first three decimals in Table 1 were always an exact match, and small deviations in the fourth decimal were occasionally observed for small $\theta$.

Again, we stress that the list of spherical codes given in Table 1 is by no means exhaustive, and the framework provided in Section 3 may be useful for evaluating any other spherical codes, if one suspects better spherical codes exist for partitioning the space.

8 Higher-dimensional codes

For high dimensions $k$, the orthant probabilities become increasingly difficult to evaluate, and the number of known, nice spherical codes quickly decreases. A well-known result is that in five or more dimensions, there exist only three classes of regular convex polytopes, which have also been studied before in the context of locality-sensitive hashing [47, 48, 4, 33, 30]. These are the simplices, orthoplices, and hypercubes. Besides these families of polytopes, we will also study the families of expanded simplices and rectified orthoplices, which correspond to the vertex figures of the $A_{k}$ and $D_{k}$ root lattices. Together with these families, we further study three more special polytopes in five and six dimensions, listed in Table 1:

• •

The $\mathbf{1_{21}}$-polytope or the $5$-demicube, obtained by taking $16$ vertices of the $5$-cube;888Viewing the hypercube as a bipartite graph, the vertices of the demicube are all vertices from one of the two sets.

• •

The $\mathbf{1_{31}}$-polytope, or the $6$-dimensional demicube on $32$ vertices;

• •

The $\mathbf{2_{21}}$-polytope or Schäfli polytope on $27$ vertices, with symmetries of the $E_{6}$ lattice.

The other polytopes in Table 1 for dimensions higher than four are all from the families $S_{k}$, $O_{k}$, $C_{k}$, $A_{k}$ and $D_{k}$ which we further discuss separately below.

8.1 Simplices ($S_{k}$)

Recall that simplices can most easily be described as the subset in $\mathbb{R}^{k+1}$ containing all standard unit vectors $\boldsymbol{e}_{1},\dots,\boldsymbol{e}_{k+1}$. We are unable to obtain closed-form expressions for $p_{1}$ and $\rho$ for $k\geq 3$, beyond the multivariate integral representation from Corollary 3.8, and we only obtain a closed-form expression of the correlation matrices appearing in the computation of $p_{1}$ below. In this case the matrix $\boldsymbol{\mathord{\mathrm{D}}}\in\mathbb{R}^{k\times k}$ has the simple shape of having ones on the diagonal and all off-diagonal entries being equal to $1/2$. (Observe that orthant probabilities for such matrices, with all off-diagonal entries equal to the same value, were studied in e.g. [46].) Yet, due to the Kronecker product, an exact evaluation of the resulting orthant probabilities for $k\geq 3$ seems difficult. In the following proposition, recall that $\boldsymbol{\mathord{\mathrm{I}}}_{k}$ denotes the identity matrix, and $\boldsymbol{\mathord{\mathrm{J}}}_{k}$ the all-$1$ matrix.

Proposition 8.1** (Collision probabilities for simplices).**

Let $\mathcal{C}\subset\mathcal{S}^{k-1}$ consist of the $k+1$ vertices of the simplex, and let $\theta\in(0,\frac{\pi}{2})$. Then the corresponding project–and–partition hash family $\mathcal{H}$ is $(\theta,p_{1},p_{2})$-sensitive, with $p_{2}=1/(k+1)$ and:

[TABLE]

Proof 8.2**.**

For ease of computations, let $\mathcal{C}=\{\boldsymbol{e}_{i}:i=1,\dots,k+1\}$. Although this code is not centered around the origin, and lies in $\mathbb{R}^{k+1}$ rather than $\mathbb{R}^{k}$, none of these things matter when computing $\boldsymbol{\mathord{\mathrm{D}}}=\hat{\boldsymbol{\mathord{\mathrm{C}}}}_{i}\hat{\boldsymbol{\mathord{\mathrm{C}}}}_{i}^{T}$, which is a function only of the difference vectors $\boldsymbol{c}_{i}-\boldsymbol{c}_{j}$. Now, a vector lies in the Voronoi cell of the vertex $\boldsymbol{e}_{1}$ if and only if it lies closer to $\boldsymbol{e}_{1}$ than to all its neighbors, which in this case are all other vectors in the code. The difference vectors $\boldsymbol{e}_{1}-\boldsymbol{e}_{j}$ all have norm $\sqrt{2}$ and pairwise inner products $1$, hence after dividing $\boldsymbol{\mathord{\mathrm{D}}}$ by $2$ we obtain a matrix with ones on the diagonal and $1/2$ everywhere off the diagonal. The result then follows from Theorem 3.8. ∎

Asymptotically, as $d\to\infty$, using simplices is equivalent to using cross-polytopes, and the last rows of Table 1 illustrate this asymptotic scaling as $k$ increases.

8.2 Orthoplices ($O_{k}$)

For cross-polytopes or orthoplices, consisting of the $2k$ code words $\pm\boldsymbol{e}_{1},\dots,\pm\boldsymbol{e}_{k}$, unfortunately we also cannot easily derive an analytic expression for the collision probabilities for $k\geq 3$. We state the reduced form of the correlation matrix below, where the matrix $\boldsymbol{\mathord{\mathrm{D}}}\in\mathbb{R}^{(2k-2)\times(2k-2)}$ is a Toeplitz matrix with ones on the diagonal, $1/2$ on most of the other diagonals, and one diagonal strip above and below the diagonal with all zeros.

Proposition 8.3** (Collision probabilities for orthoplices).**

Let $\mathcal{C}\subset\mathcal{S}^{k-1}$ consist of the $2k$ vertices of the orthoplex, and let $\theta\in(0,\frac{\pi}{2})$. Then the corresponding project–and–partition hash family $\mathcal{H}$ is $(\theta,p_{1},p_{2})$-sensitive, with $p_{2}=1/(2k)$ and:

[TABLE]

Proof 8.4**.**

Without loss of generality, let $\mathcal{C}=\{\pm\boldsymbol{e}_{j}:j=1,\dots,k\}$. For the vector $\boldsymbol{c}_{1}=\boldsymbol{e}_{1}$, the relevant vectors are $\mathcal{R}_{1}=\mathcal{C}\setminus\{\pm\boldsymbol{e}_{1}\}$, and taking differences with the relevant vectors, we observe that the matrix $\hat{\boldsymbol{\mathord{\mathrm{C}}}}_{1}$ has $2k-2$ rows $\boldsymbol{e}_{1}\pm\boldsymbol{e}_{j}$, for $j=2,\dots,k$. For the product $\boldsymbol{\mathord{\mathrm{D}}}=\hat{\boldsymbol{\mathord{\mathrm{C}}}}_{1}\hat{\boldsymbol{\mathord{\mathrm{C}}}}_{1}^{T}$ we get twos on the diagonal, ones for most off-diagonal entries, and each row is paired up with exactly one other row whose joint correlation is [math]. Through a suitable reordering of the rows/columns of $\boldsymbol{\mathord{\mathrm{D}}}$, we thus get the claimed shape of $\boldsymbol{\mathord{\mathrm{D}}}$. Finally, with $\mathcal{C}$ being isogonal and of size $2k$, using Theorem 3.8 we get the claimed expressions for $p_{1}$ and $p_{2}$. ∎

As $k\to\infty$, from [4] we know that the parameters $\rho$ scale optimally as $\rho\to 1/(2c^{2}-1)$, where $c=\sqrt{1/(1-\cos\theta)}$ is the approximation factor for worst-case approximate nearest neighbor searching. Substituting values $\theta\in\frac{\pi}{12}\{1,2,3,4,5\}$, we obtain the values in the bottom rows of Table 1.

8.3 Hypercubes ($C_{k}$)

The easiest of these three families of convex regular polytopes to analyze is the family of hypercubes. For arbitrary $k\geq 1$, these codes consist of $2^{k}$ vertices of the form $\frac{1}{\sqrt{k}}(\pm\boldsymbol{e}_{1}\pm\dots\pm\boldsymbol{e}_{k})$. When using normal matrices $\boldsymbol{\mathord{\mathrm{A}}}\sim\mathcal{N}(0,1)^{k\times d}$ to project down onto a $k$-dimensional subspace, the following proposition shows that the query exponent does not change at all as $k$ increases.

Proposition 8.5** (Collision probabilities for hypercubes).**

Let $\mathcal{C}\subset\mathcal{S}^{k-1}$ consist of the $2^{k}$ vertices of the hypercube, and let $\theta\in(0,\frac{\pi}{2})$. Then the corresponding project–and–partition hash family $\mathcal{H}$ is $(\theta,p_{1},p_{2})$-sensitive, with:

[TABLE]

As a result, when using Gaussian projection matrices $\boldsymbol{\mathord{\mathrm{A}}}$ the query exponents $\rho$ for hypercubes are equal to those for the antipodal spherical code of Proposition 4.1.

Proof 8.6**.**

As before, scaling $\mathcal{C}$ by a constant factor does not affect collision probabilities, so w.l.o.g. we may let $\mathcal{C}$ consist of the $2^{k}$ vectors $(\pm 1,\pm 1,\dots,\pm 1)$. A vector $\boldsymbol{x}$ is mapped to the code word $\boldsymbol{c}_{0}=(1,1,\dots,1)$ if and only if it is closer to this code word than to its relevant vectors $\mathcal{R}_{0}$, which in this case are the vectors with $(k-1)$ entries equal to $+1$ and one entry equal to $-1$. For $j=1,\dots,c$ let us denote by $\boldsymbol{c}_{j}$ the vector with ones everywhere, except at position $j$ where we have a $-1$. For the difference vectors $\boldsymbol{c}_{0}-\boldsymbol{c}_{j}$ with $\boldsymbol{c}_{j}\in\mathcal{R}_{1}$ we get $\boldsymbol{c}_{0}-\boldsymbol{c}_{j}=2\boldsymbol{e}_{j}$, and as a result we obtain $\hat{\boldsymbol{\mathord{\mathrm{C}}}}_{0}\hat{\boldsymbol{\mathord{\mathrm{C}}}}_{0}^{T}=4\cdot\boldsymbol{\mathord{\mathrm{I}}}_{k}$, which is a scalar multiple of the identity matrix. Now, note that a correlation matrix $\boldsymbol{\mathord{\mathrm{I}}}_{k}$ essentially means that the entries of the random vector are independent. It is therefore an easy exercise to see that, for all matrices $\boldsymbol{\mathord{\mathrm{A}}}$:

[TABLE]

Together with the derivation of the orthant probability for the one-dimensional antipodal spherical code, this leads to $p_{1}=2^{k}\cdot(\pi-\theta)^{k}/(2\pi)^{k}=(1-\theta/\pi)^{k}$. As $\mathcal{C}$ is isogonal, we further obtain $p_{2}=1/2^{k}$ as claimed. ∎

Note that the paper [33] showed that, if $\boldsymbol{\mathord{\mathrm{A}}}\in\mathbb{R}^{d\times d}$ is chosen as a random orthogonal matrix, then in fact the performance of hypercube hashing gets better as $k$ increases. Here however we assumed $\boldsymbol{\mathord{\mathrm{A}}}\sim\mathcal{N}(0,1)^{k\times d}$, which for small $k$ is essentially equivalent to sampling $\boldsymbol{\mathord{\mathrm{A}}}$ as an orthogonal matrix, but is not quite the same for large $k$; see e.g. [23, 29] for more details on how the gap between orthogonal and Gaussian matrices changes with the relation between $k$ and $d$. This leads to two different asymptotic scalings, as indicated in the bottom rows of Table 1.

8.4 Expanded simplices ($A_{k}$)

Similar to simplices, so-called expanded simplices (the vertex figures of the $A_{k}$ root lattices) can most easily be described as a subset in $\mathbb{R}^{k+1}$. These spherical codes contain all differences between unit vectors, $\boldsymbol{e}_{i}-\boldsymbol{e}_{j}$ for $i\neq j$. These vertices all lie on a $k$-dimensional hyperplane, and can thus be viewed as a $k$-dimensional spherical code. For these codes the matrix $\boldsymbol{\mathord{\mathrm{D}}}\in\mathbb{R}^{(2k-1)\times(2k-1)}$ has a block structure described below.

Proposition 8.7** (Collision probabilities for expanded simplices).**

Let $\mathcal{C}\subset\mathcal{S}^{k-1}$ consist of the $k(k+1)$ vertices of the expanded simplex, and let $\theta\in(0,\frac{\pi}{2})$. Then the corresponding project–and–partition hash family $\mathcal{H}$ is $(\theta,p_{1},p_{2})$-sensitive, with $p_{2}=1/k(k+1)$ and:

[TABLE]

We omit the proof, which is analogous to the proofs of Propositions 8.1 and 8.3.

Observe that, as hashing in the expanded simplex corresponds to mapping to the two largest coordinates of a vector, for large $d$ up to order terms this is essentially equivalent to using two independent simplex hashes. For expanded simplices we therefore again get the optimal scaling for large $k$ also obtained with orthoplices and simplices, as illustrated in Table 1.

8.5 Rectified orthoplices ($D_{k}$)

Rectified orthoplices, related to the root lattices $D_{k}$, contain as vertices all pairwise combinations of unit vectors; in contrast to $A_{k}$, where the signs must be opposite, here any combination of signs is included. As a subset of $\mathbb{R}^{k}$, this code thus contains the $2k(k-1)$ vertices of the form $\pm\boldsymbol{e}_{i}\pm\boldsymbol{e}_{j}$ for $i\neq j$. For these codes the matrix $\boldsymbol{\mathord{\mathrm{D}}}\in\mathbb{R}^{(2k-1)\times(2k-1)}$ has a $2\times 2$ block structure described below.

Proposition 8.8** (Collision probabilities for rectified orthoplices).**

Let $\mathcal{C}\subset\mathcal{S}^{k-1}$ consist of the $2k(k-1)$ vertices of the rectified orthoplex, and let $\theta\in(0,\frac{\pi}{2})$. Then the corresponding project–and–partition hash family $\mathcal{H}$ is $(\theta,p_{1},p_{2})$-sensitive, with $p_{2}=1/2k(k-1)$ and:

[TABLE]

All matrices $\boldsymbol{\mathord{\mathrm{I}}}=\boldsymbol{\mathord{\mathrm{I}}}_{k-2}$ and $\boldsymbol{\mathord{\mathrm{J}}}=\boldsymbol{\mathord{\mathrm{J}}}_{k-2}$ above are square matrices of size $(k-2)\times(k-2)$.

We omit the proof, which is analogous to the proofs of Propositions 8.1 and 8.3. Note again that, as we are mapping to the two largest coordinates after projection, the hash function for the rectified orthoplex has the same asymptotic scaling as the simplex and orthoplex.

8.6 Comparison

To compare these families of spherical codes in terms of their performance for nearest neighbor searching, Figure 2 depicts the resulting parameters $\rho$ for different values $k$, focusing on the two cases $\theta=\frac{\pi}{3}$ and $\theta=\frac{\pi}{4}$. As observed previously there is a cross-over between simplices and orthoplices around $k=4$, after which the orthoplex consistently performs better than the simplex. The hypercube, as shown in Proposition 8.5, remains stable at the same performance as the antipodal code, regardless of the choice of dimension $k$.

9 Utopian spherical cap estimates

To get an idea how good the spherical codes considered above are, ideally we would find a matching lower bound, stating that any spherical code of dimension at most $k$ must satisfy $\rho(\theta)\geq\rho_{k}(\theta)$, for some monotonically increasing function $\rho_{k}:(0,\frac{\pi}{2})\to(0,1)$. However, ruling out the existence of exotic spherical codes which happen to perform better than the ones we considered seems difficult.

As has been observed several times before [5, 13, 6], using spherical caps to partition the sphere into regions seems to work best, so that the regions have nice, symmetric shapes that ultimately seem to minimize the exponent $\rho$. We can therefore get lower bounds on the parameter $\rho$ by assuming that the sphere can be perfectly divided into an integer number of spherical caps, and estimating the resulting parameter $\rho$ for these utopian999Note that for $k\geq 3$ it is impossible to partition the sphere into spherical caps, as regions either overlap, or part of the sphere remains uncovered. These conjectured lower bounds therefore do not correspond to actual, achievable bounds, and are likely not tight. partitions of the sphere. This leads to the following result, where we assume such perfect spherical cap partitions exist, and where we compute collision probabilities $p_{1},p_{2}$ and the query exponents $\rho$ accordingly.

Theorem 9.1** (Spherical cap lower bounds).**

Let $\mathcal{C}\subset\mathcal{S}^{k-1}$ be a spherical code, and let $\mathcal{H}$ be the associated project–and–partition hash family. Then the parameter $\rho$ for $\mathcal{C}$, for angle $\theta\in(0,\frac{\pi}{2})$, must satisfy:

[TABLE]

where $\rho_{k}(c,\theta)$ is given by:

[TABLE]

and $\alpha_{k}(c)$ is the solution to $\frac{1}{2}\cdot I_{1-\alpha^{2}}(\frac{k-1}{2},\frac{1}{2})=\frac{1}{c}$.

Proof 9.2**.**

Suppose we use a spherical code in $k$ dimensions, partitioning the $k$-sphere into $c$ regions $\mathcal{V}_{1},\dots,\mathcal{V}_{c}$. First we will argue that the performance parameter $\rho$ of any such spherical code is at most the corresponding parameter where each region $\mathcal{V}_{i}$ was replaced by a spherical cap of the same volume. Then we will argue that if indeed we used spherical caps to form our code, it is optimal to use spherical caps of equal size. Finally we will tie up the loose ends to derive the stated lower bounds.

Proof, part 1: Spherical caps are optimal.

First, suppose we fix the size of the code $c$ as well as the volumes of each Voronoi region. It is clear that the collision probability $p_{2}=\sum_{i=1}^{c}\operatorname{Vol}(\mathcal{V}_{i})^{2}/\operatorname{Vol}(\mathcal{S}^{k-1})^{2}$ remains invariant under any changes we make to $\mathcal{C}$ adhering to these constraints, and so the only question is whether $p_{1}$ is largest when each region is a spherical cap.

To prove this we use a classical result of Baernstein–Taylor **[12]**, which extends Riesz’ inequality **[41]** to the Euclidean sphere. Suppose we have functions $f,g:\mathcal{S}^{d-1}\to\mathbb{R}$ and $h:[-1,1]\to\mathbb{R}$, and let the surface measure $\sigma$ on $\mathcal{S}^{d-1}$ be normalized so that $\sigma(\mathcal{S}^{d-1})=1$. Let $h$ further be non-decreasing, bounded, and measurable, and let $f,g\in L_{1}(\mathcal{S}^{d-1})$. Let:

[TABLE]

Then **[12, Theorem 2]** states:

[TABLE]

where $f^{*}=f^{*}(x_{1})$ is a non-decreasing function which only depends on the first coordinate $x_{1}$, $g^{*}=g^{*}(y_{1})$ is a non-decreasing function which only depends on $y_{1}$, and $\sigma(\{f^{*}>\lambda\})=\sigma(\{f>\lambda\})$ and $\sigma(\{g^{*}>\lambda\})=\sigma(\{g>\lambda\})$ for all $\lambda\in\mathbb{R}$.

To instantiate it, for each Voronoi region $i=1,\dots,c$, we apply **[12, Theorem 2]** with $f(\boldsymbol{x}):=\mathds{1}\{\boldsymbol{x}\in\mathcal{V}_{i}\}$, $g(\boldsymbol{x}):=\mathds{1}\{\boldsymbol{y}\in\mathcal{V}_{i}\}$, and $h(\boldsymbol{x}\cdot\boldsymbol{y})=\mathds{1}\{\boldsymbol{x}\cdot\boldsymbol{y}\geq\cos\theta\}$. Note that for these functions $f$ and $g$ we have $f^{*}(\boldsymbol{x})=\mathds{1}(\boldsymbol{x}\cdot\boldsymbol{e}_{1}\geq\alpha_{i}\}$ and $g^{*}(\boldsymbol{x})=\mathds{1}(\boldsymbol{x}\cdot\boldsymbol{e}_{1}\geq\alpha_{i}\}$ with $\alpha_{i}$ depending on the volume of $\mathcal{V}_{i}$: the functions $f^{*}$ and $g^{*}$ correspond to spherical caps $\mathcal{C}_{i}$.

[TABLE]

On the other hand, for $\mathcal{T}(f^{*},g^{*},h)$ we obtain:

[TABLE]

By the inequality $\mathcal{T}(f,g,h)\leq\mathcal{T}(f^{*},g^{*},h)$ it follows that each such probability is upper bounded by the corresponding spherical cap probability. Thus, if our original hash function $h$ uses a spherical code $\mathcal{C}$ with Voronoi regions $\mathcal{V}_{i}$, and a utopian spherical code $\mathcal{C}^{*}$ were to consist of induced Voronoi regions in the shape of spherical caps $\mathcal{C}_{i}$ and generated a hash function $h^{*}$, then the previous results would show that:

[TABLE]

In other words, in a utopian world we would always replace hash regions with spherical caps of the same volume if we could, as the resulting parameter $\rho$ can only get smaller.

Proof, part 2: Equal spherical caps are optimal.

The previous derivation shows that the parameter $\rho$ can be lower-bounded by replacing each region with a spherical cap of the same volume. The next question is: if we do use such spherical cap-shaped hash regions, does it ever make sense to use regions of different volumes?

In the above derivation, where each spherical region is independent of the others (we assume in our utopian setting that the spherical caps are perfectly disjoint, and the only boundary conditions for regions $\mathcal{V}_{i}$ are defined by $\alpha_{i}$), the performance parameter $\rho$ of the entire code is a weighted average of the individual performance parameters $\rho$ of each of these spherical caps on parts of the sphere. This weighted average is minimized when we choose all the regions to be of the size that minimizes $\rho$, i.e. we find the spherical cap size that offers the best balance between the collision probabilities $p_{1}$ and $p_{2}$ and use that cap size for all our spherical caps.101010Note that the resulting optimal code size may technically not be integral, and so the optimal integral-sized code may not consist of equal regions. However, minimizing over real-valued $c$ will always provide a lower bound on the parameter $\rho$ of any $k$-dimensional spherical code.

Proof, part 3: Obtaining the final expression.

*Finally, observe that w.l.o.g. we may model two vectors $\boldsymbol{v},\boldsymbol{w}$ at angle $\theta$ in $d$-dimensional space by the two vectors $\boldsymbol{v}=\boldsymbol{e}_{1}$ and $\boldsymbol{w}=\boldsymbol{e}_{1}\cos\theta+\boldsymbol{e}_{2}\sin\theta$. A projection via a random Gaussian matrix $A\sim\mathcal{N}(0,1)^{d\times k}$, whose first two columns are denoted $\boldsymbol{x}$ and $\boldsymbol{y}$, then projects $\boldsymbol{v}$ onto $A\boldsymbol{v}=\boldsymbol{x}$ and $\boldsymbol{w}$ onto $A\boldsymbol{w}=\boldsymbol{x}\cos\theta+\boldsymbol{y}\sin\theta$, with $\boldsymbol{x},\boldsymbol{y}\sim\mathcal{N}(0,1)^{k}$. The condition of being contained in a spherical cap then requires an additional normalization, as these vectors $\boldsymbol{x}$ and $\boldsymbol{y}$ are not necessarily normalized. This finally leads to the given expression for $\rho_{k}(c,\theta)$, which denotes the parameter $\rho$ for a utopian code consisting of $c$ equal-sized spherical caps. Note that the given parameter $\alpha_{k}(c)$ is such that $\Pr_{\boldsymbol{z}\in\mathcal{S}^{k-1}}(z_{1}\geq\alpha_{k}(c))=1/c$. *

For analyzing actual spherical codes, we obtained conditions of the form $\|\boldsymbol{x}-\boldsymbol{c}_{i}\|\geq\|\boldsymbol{x}-\boldsymbol{c}_{j}\|$, which then translated to conditions $\langle\boldsymbol{x},\boldsymbol{c}_{i}-\boldsymbol{c}_{j}\rangle\geq 0$, leading to orthant probabilities. Here however we have different conditions of the form $\|\boldsymbol{x}/\|\boldsymbol{x}\|-\boldsymbol{c}_{i}\|\leq u$ for $\boldsymbol{x}$ to be contained in the Voronoi cell (spherical cap) defined by code word $\boldsymbol{c}_{i}$. Assuming w.l.o.g. $\boldsymbol{c}=\boldsymbol{e}_{1}$, the formulas reduce to those given in (46), but these unfortunately cannot be further simplified or expressed in terms of orthant probabilities.

For small values $k=1,\dots,6$, and for the angles $\theta\in\frac{\pi}{12}\{1,\dots,5\}$, we numerically computed values $\rho_{k}(c,\theta)$, and taking the minimum over $c$, we obtained the values $\rho_{k}(\theta)$ given in Table 1. The superscripts in Table 1 indicate the value $c$ leading to the minimum value $\rho$, i.e. the ideal code size for these utopian spherical cap codes. Although these codes are not actually instantiable, these lower bounds may be useful as a guideline for choosing the code size $c$, when trying to come up with the best spherical codes in this framework.

From Table 1, observe that as the angle $\theta$ decreases, and the nearest neighbor is expected to lie very close to the query vector, the optimal code size $c$ increases; for problem instances where the neighbor lies very close to the query vector, a fine-grained partition will likely work better. Also note that as $k$ increases, the optimal code size $c$ for these utopian codes increases as well, and whereas for small $k$ the optimal code size is always less than the number of vertices of the simplex and orthoplex, for larger $k$ the optimal code size increases beyond that of the simplex, and also beyond that of the orthoplex. This is in line with our observation that as $k$ increases, orthoplices start outperforming simplices, and this also hints at the idea that as $k$ increases further, even denser codes (likely with a superlinear number of vertices) will exhibit a better performance (smaller query exponents $\rho$) than partitions based on the cross-polytopes.

Conjecture 9.3** (Orthoplices are suboptimal for large $k$).**

There exists a dimension $k_{0}\in\mathbb{N}$ such that, for all dimensions $k\geq k_{0}$, there exist spherical codes $\mathcal{C}\subset\mathcal{S}^{k-1}$ whose query exponents $\rho$ in the project–and–partition framework are smaller than the exponents $\rho$ of the $k$-orthoplex.

An interesting open question, which may well be closely related to other sphere packing and covering problems, is to find higher-dimensional spherical codes that are more suitable for nearest neighbor tasks than using the vertices of simplices or orthoplices.

10 Choosing a code in practice

With the above framework in mind, and looking at the results from Table 1 and Figures 1 and 2, one might ask: how do I choose the best spherical code for my application? Is this purely a matter of choosing the code with the smallest exponent $\rho$? Or is the practical performance affected by other aspects of the code as well?

To look for an answer to this question, recall that in locality-sensitive hashing (as in Equation (10)) there are several parameters that play a role:

[TABLE]

Here the constant in $t$ determines the success probability of the scheme. Now, the concrete costs for finding nearest neighbors with locality-sensitive hashing can be summarized as follows:

Projections:

For a $k$-dimensional code, we use Gaussian projection matrices $\boldsymbol{\mathord{\mathrm{A}}}\sim\mathcal{N}(0,1)^{k\times d}$, and without further optimizations to $\boldsymbol{\mathord{\mathrm{A}}}$, computing a product $\boldsymbol{\mathord{\mathrm{A}}}\boldsymbol{x}$ requires $k\cdot d$ multiplications. Generating such a matrix $\boldsymbol{\mathord{\mathrm{A}}}$ can also be done in time proportional to $k\cdot d$.

Decoding:

Although for some spherical codes faster decoding algorithms may exist, in general the decoding cost for a spherical code of size $c$ in dimension $k$ is proportional to $k\cdot c$ multiplications.

Hash tables:

We need $t=O(n^{\rho})$ hash tables for our data structure.

Hash functions:

For each of the $t$ tables, we need to initialize $\ell$ random hash functions, for a cost of $\ell\cdot t$ hash functions. Each of these requires initializing a random projection matrix $\boldsymbol{\mathord{\mathrm{A}}}\sim\mathcal{N}(0,1)^{k\times d}$, which can all be generated in time proportional to $\ell\cdot t\cdot k\cdot d$ multiplications.

Preprocessing time:

For each of $n$ vectors, and each of $t$ hash tables, we need to compute the $k$-concatenated hash value, and insert the vector in the corresponding bucket. Timewise, computing all hash values for all vectors takes time proportional to $\ell\cdot t\cdot k\cdot d\cdot n$.

Memory requirement:

We need to store $t$ hash tables, each of which contains all $n$ vectors, for a total cost of $t\cdot n$ entries in the hash tables.111111Observe that we can keep a single copy of the entire list in memory, and only store pointers of essentially constant space to these vectors in the hash tables.

Query time – Hashing:

Given a vector, we need to compute its hash values in all $t$ hash tables, for a cost of $\ell\cdot t\cdot k\cdot d$ multiplications.

Query time – Comparisons:

In total we expect that, apart from the nearest neighbor, approximately $t$ vectors will collide with a query vector in the $t$ hash tables. Going through these vectors to find the actual nearest neighbor can be done in time proportional to $t\cdot d$.

Here, an important observation is that the projections generally become more costly as the dimension $k$ of the spherical code increases; in higher dimensions, spherical codes with better parameters $\rho$ exist, but at the same time we pay for these more complex spherical codes in the decoding costs. The highest time and space complexities for the different parts of the nearest neighbor data structure can be summarized as follows:

[TABLE]

For the query time complexity $T_{\text{query}}$, the factor $d$ is constant and does not depend on the choice of the spherical code $\mathcal{C}$. Similarly, the term $\ell$ contains a factor $\log n$ which does not depend on the choice of $\mathcal{C}$. Dividing out a factor $n\log_{2}n$, the product $t\cdot\ell\cdot k/\log_{2}n$ has the following form for isogonal codes $\mathcal{C}$:

[TABLE]

Now, for large $n$ this product is clearly minimized especially when $\rho$ is smallest. For concrete values of $n$ however, the other main factor that determines the query complexity is the ratio $\log_{2}(c)/k$. Here $\log_{2}(c)$ can be seen as the number of hash bits that are extracted from a single hash value $h(\boldsymbol{x})\in\{1,\dots,c\}$, while the factor $k$ comes from the cost of computing a projection onto a $k$-dimensional subspace. In general, codes with a large ratio $\log_{2}(c)/k$ require fewer random projections to extract the necessary number of hash bits to determine the hash bucket, while small codes in high dimensions are in that sense rather costly. This quantitatively explains why e.g. hyperplane hashing is fast in practice, as $\log_{2}(c)/k=1$, but cross-polytope hashing has a high cost for computing projections, as $\log_{2}(c)/k=O(\log(k)/k)$ is almost a linear factor smaller.

To compare all codes in terms of their extracted number of bits per projection, as well as their query exponents $\rho$, Figure 3 shows the resulting trade-offs for the spherical codes considered in this paper. As we can see, there is often a concrete trade-off between the query exponent $\rho$ and the ratio of extracted hash bits $\log_{2}(c)/k$; codes with small $\rho$ generally have a high hash cost, and codes which partition a low-dimensional space in many regions pay the price with a higher value $\rho$. Figure 3 also paints a different picture as to which codes are best suited in practice, compared to Table 1; the lattice-based codes $D_{k}$ may be better choices in practice than the pure orthoplex hashing of [4], achieving better values $\rho$ while keeping the extracted hash bits per projection vector larger than for orthoplex hashing.

For an even more concrete comparison of different codes in practice, Figure 4 describes concrete cost estimates for processing queries for different spherical codes, for different parameters $n$. Here the query cost estimate is simplified to $t\cdot\ell\cdot k\cdot d$, without taking into account further multiplicative order terms (which are likely to be roughly the same for different spherical codes). As the figure shows, for large $n$ the query exponent $\rho$ dominates, and codes with small values $\rho$ achieve the best query cost estimates. For smaller $n$ (e.g. $n=10^{5}$, which might be the size of the database in applications in cryptanalysis [31, 13]), the cross-over point with hyperplane hashing quickly goes down, and using full-size cross-polytopes may not be as practical as using hyperplanes to partition the space; the smaller $\rho$ of cross-polytope hashing does not outweight the higher cost of the many random projections. This matches practical results from e.g. [31, 14, 13, 36, 37], which established that in applications with rather small $n$, hyperplane hashing dominates. For other applications, where $n$ is larger, using $D_{k}$ lattices or e.g. the $2_{21}$-polytope seems to give the best performance.

As a generalization of the (rectified) orthoplex families of hash functions $O_{k}$ and $D_{k}$, we finally propose the following hash function families, mapping a vector to the $m$ absolute largest coordinates.

Definition 10.1** ($m$-max locality-sensitive hash functions).**

Let $1\leq m\leq k\leq d$. Then we define the $m$-max hash function family $\mathcal{H}$ as the project–and–partition hash family with code:

[TABLE]

For arbitrary $k$ and $m$, this spherical code has size $c=2^{m}\binom{k}{m}$.

As special cases, $m=1$ corresponds to orthoplex hashing ($O_{k}$), $m=2$ corresponds to the rectified orthoplices ($D_{k}$), and $m=k$ corresponds to hypercube hashing ($C_{k}$). As Figure 3 suggests that larger codes ($m=2$) may work better in practice than smaller codes ($m=1$), we predict that for large $k$, larger values $m$ will outperform both orthoplex hashing and rectified orthoplex hashing. An interesting open question would be to estimate how $m$ should ideally scale with $k$ to obtain the best results in practice.