Local Search Methods for Fast Near Neighbor Search
E.S. Telleza,b , G. Ruizc , E. Chavez∗,d , M. Graffa,b

arXiv:1705.10351v2 [cs.DS] 31 May 2017

a CONACyT

Consejo Nacional de Ciencia y Tecnologı́a, Dirección de Cátedras, Insurgentes
Sur 1582, Crédito Constructor 03940, Ciudad de México, México.
b INFOTEC Centro de Investigación e Innovación en Tecnologı́as de la Información y
Comunicación, Circuito Tecnopolo Sur No 112, Fracc. Tecnopolo Pocitos II, Aguascalientes
20313, México.
c Universidad Michoacana de San Nicolás de Hidalgo, México
d Centro de Investigación Cientı́fica y de Educación Superior de Ensenada, Baja California,
México.

Abstract
Near neighbor search (NNS) has been traditionally addressed from an algorithmic perspective, that is, given a dataset and a distance function, the goal
is to build a data structure along with discarding rules that quickly find the
near neighbor of a query under the given distance function. In this manuscript,
we take another approach. We restate NNS as a combinatorial optimization
problem, being the goal to minimize the distance from the query to the data
set. This cost function is well defined, and the minimum is realized precisely
with the nearest object to the query. Adopting this new view allows the use of
a rich collection of optimization algorithms with a long tradition.
We present three local search algorithms that solve the NNS problem from
the combinatorial optimization point of view. Two of the algorithms can be
described as the hyper-parameter optimization of two known indexing methods,
simplifying their adoption for final users. The third method achieves excellent
search tradeoffs among speed, memory, and quality of the approximation. For
instance, for moderately large dimensions, our indexes can reach above 0.99
recall reviewing less than 1% of the database.
We conducted extensive experimentation with five real-world standard benchmark and five synthetic datasets to prove our claims.
Key words: Near neighbor search, Combinatorial optimization, Approximate
algorithms.

1. Introduction
Nearest neighbor search is a pervasive problem in computer science. The
applications range from pattern recognition to textual and multimedia infor∗ Corresponding

author: E. Chavez, email: elchavez@cicese.mx.
Email addresses: eric.tellez@infotec.mx (E.S. Tellez), gruiz@dep.fie.umich.mx (G.
Ruiz), elchavez@cicese.mx (E. Chavez), mario.graff@infotec.mx (M. Graff)

Preprint submitted to Elsevier

June 1, 2017

mation retrieval, machine learning, compression, biometric identification and
authentication, and bioinformatics [37].
Given a query object q, the nearest neighbor search consists in retrieving
the set k nn(q) of k closer items to the query from a finite set S ⊂ U from
a possibly infinite set U , let n = |S|. The result set of k nn(q) is a ball of
radius rk , B(q, rk ) ⊆ S, containing those items in S being close to q; ties are
arbitrarily broken. Closeness between objects is measured with some criterion,
say a distance function d : U × U → R+ which is usually required to obey the
metric properties: d(u, v) = 0 ⇐⇒ u = v, d(u, v) = d(v, u), and the triangle
inequality |d(u, v) − d(v, w)| ≤ d(u, w) ≤ d(u, v) + d(v, w). A search can be
trivially solved just evaluating all distances from q to each object in S; however,
this approach is not of use when S is large, or whenever the cost of evaluating
the distance function makes this approach prohibitive.
A vector space V of dimension δ is a set of tuples of δ numbers (a subset
of Rδ ). There are many possible distances for vectors, being the Euclidean or
Manhattan distances the most popular. Both are Minkowski distances (Lp )
defined in general as follows:

Lp (u, v) =

δ
X

!1/p
|ui − vi |

p

(1)

i=1

The Euclidean distance is obtained for p = 2, and the Manhattan distance for
p = 1. In the limit where p → ∞, we have the Chebyshev distance L∞ (u, v) =
max1≤i≤δ |ui − vi |.
In other domains, for example in Information Retrieval and Natural Language Processing, TFIDF vectors are traditional representations of documents.
TFIDF vectors use the cosine similarity,
Pδ
simC (u, v) = qP

δ
i=1

i=1

u2i ·

ui vi
qP
δ

i=1

(2)
vi2

Metric methods should use ∠(u, v) = cos−1 simC (u, v), that is, the angle between
vectors. We use this distance in our Documents benchmark, see §4.
Another distance of interest in this manuscript is the Levenshtein or edit
distance, defined between strings of symbols. It computes the number of edit
operations needed to convert one string into another. Edit operations are insertion, deletion, and substitution. For strings u, v of sizes i, j respectively it is
defined as follows:


max(i, j) if min(i, j) = 0





ed(ui−1 , vj )
+1






ed(ui , vj−1 )
+ 1(
ed(ui , vj ) =

min



0 if last(ui ) = last(vj )





ed(ui−1 , vj−1 ) +

1 if last(ui ) 6= last(vj )

2

where last(·) is a function that returns the last symbol of a string. The edit
distance is used for the Wiktionary benchmark in §4. A related distance, the
Longest Common Subsequence (LCS) distance can be defined in terms of ed
just setting to 2 the cost when last(ui ) 6= last(vj ). The effect is to allow only
insertions and deletions.
Indexes for multidimensional data like quad-trees, oct-trees, or kd-trees,
make use of the coordinates to partition either the data or the space in such
a way that items are pruned by checking numeric values per dimension. This
approach is convenient whenever δ remains low since it can achieve logarithmic
times on n when δ is considered a constant value [34]. Partition by dimension
does not work in high dimensions because the algorithms have an exponential
dependency on δ, this is folklore and is often called the curse of dimensionality
in the literature.
On the other hand, metric access methods claim to be independent of the
representation of the objects. They use the distance function to partition the
data. The metric approach is useful, in principle, on datasets of high explicit
dimension and small intrinsic dimension. In high dimensional metric spaces the
distribution of the distance among items is concentrated around the mean, and
hence, the partitions have a lot of overlap. Therefore, the curse of dimensionality appears again, in disguise. This effect has been studied both in theory
and in practice and has been reported in the literature. The net effect of the
curse of dimensionality, either in vector or metric spaces, is a slowdown in the
performance at query time. Even for highly selective queries, it is necessary to
check virtually the entire database.
1.1. Approximate and Probabilistic Indexes
On intrinsically high dimensional datasets, when the curse of dimensionality
prevents the use of exact algorithms, a more promising approach to have a
fast search is the use of approximate and probabilistic algorithms. Trading
accuracy for speed gives roominess in the design of algorithms. For example,
instead requesting exactly the k nearest neighbors, we allow two types of error,
i) Not reporting some relevant results, or ii) Reporting objects not on the query
outcome. Since the results will contain errors, we need some way to measure the
quality of the result set. For this purpose, we define the recall of an approximate
result k nn∗ (q) as
recall(k nn∗ (q)) =

|k nn∗ (q) ∩ k nn(q)|
|k nn(q)|

It is not hard to modify the definition of recall to take into account distance
accuracy (how far the result is from the actual answer). This alternate definition
is of use when comparing results with a large number of ties, e.g. when the
distance is discrete.
There exist a number of generic methods to turn an exact metric index into
an approximate and/or probabilistic index [12, 43, 8, 36]. However, indexes
designed to be approximate scale better with both the intrinsic dimension and
3

the database size [15, 10, 20, 33, 39, 2, 16, 38, 40, 41, 42]. Below, we describe
some relevant examples.
Locality Sensitive Hashing (LSH) is a set of fast approximate searching techniques giving probabilistic guarantees on the quality of the result set, Gionis et
al. [20]. The general idea of an LSH index is to find hashing functions such that
close items share the same bucket with high probability, while distant objects
share the same bucket with low probability.
More formally, a family of hashing functions H = {g1 , g2 , . . . , gh }, gi : U →
{0, 1} is called (p1 , p2 , r1 , r2 )-sensitive if, for any p, q ∈ U ,
— if d(p, q) < r1 then P r[hash(p) = hash(q)] > p1 ,
— if d(p, q) > r2 then P r[hash(p) = hash(q)] < p2 ,
where hash(u) is the concatenation of the output of individual hashing functions
gi following a fixed order, that is, hash(u) = g1 (u)g2 (u) · · · gh (u).
Notice that LSH gives probabilistic guarantees about the maximum distance
between the actual answer and the answer provided by the algorithm, and it
cannot guarantee the quality regarding the recall. A large number of empty
buckets will appear to achieve the desired p1 and p2 values (for potential queries
of U that are far from any element in S). This behavior stress the memory
footprint of the algorithm, because to solve it several hashing functions (a so
called hashing family) are needed. Still, high-quality indexes usually need high
amounts of memory.
LSH is well defined for general metric spaces, but the process of finding
suitable hashing functions gi is not trivial and has to be engineered for each case.
There exist simple LSH functions for the Euclidean (L2 ), Hamming, Jaccard,
and Cosine distances [3, 20]. Also, there exist several approaches extending the
concept of LSH to be data sensitive. In general, the problem consist in selecting
hashing functions suitable for the data and the query set [19].
The Fast Library for Approximate Nearest Neighbors, FLANN [28], is a popular library that implements different indexing strategies and selects the more
suitable one, based on parameters like the search time, the construction time,
and the memory available. With those parameters, together with the parameters of the indexes, a cost function is defined. FLANN builds and tries different
instances of the indexes using optimization techniques to find a local minimum
to the cost function. The instance corresponding to the local minimum will be
the one selected.
The indexing methods used by the FLANN are: Randomized k-d forest,
priority search k-means tree, and the sequential search. The randomized k-d
forest consists of multiple k-d trees searched in parallel [35] where the number
of visited leafs on each tree is limited to control the search cost. The k-means
tree [28] is a data structure computing k-means clustering on the dataset producing k compact regions. At each region, the procedure is repeated recursively,
producing a tree of arity k. The recursion is stopped when regions have less than
k items. The search procedure has a best first strategy, and the algorithm uses
a priority queue to order the paths. The performance is adjusted controlling the
4

expansion factor (number of visited regions). As commented, their parameters
are optimized, and the user only needs to specify the desired recall.
Recently, a set of small and fast indexes have appeared, among the remarkable ones we found CNAPP [42], PP-Index [17], MIF [1], and the quantized
permutations [26, 27]. In these indexes, a similar structure arises: every node is
associated with a set of k nearest references (KNR), where the set of references
is a sample of the database. The similarity between items is hinted by a similarity function over the shared neighbors. In this approach, there is no navigation
and proximity queries are solved using an inverted index. This structural similarity was systematically explored in [11], adding several new indexes to the
list. In the same paper, the authors include experimental comparisons showing
the performance of KNR indexes as compared against near neighbor techniques
based on hashing, in particular, LSH [23, 4], and DSH [19]. The best indexes
in the KNR framework are KNR-Cos and KNR-LCS. Those indexes use the
cosine similarity and the longest common sub-sequence as similarity functions
as proximity predictors, respectively. The experimental evidence suggests that
KNR-Cos surpasses the majority of the state of the art KNR indexes, so we will
use it as the canonical member of the KNR family in our experimental section.
In [21] near neighbors are found using a random walk over a precomputed
table of neighbors. The required preprocessing time and memory for the table
are quadratic, preventing the scalability of the approach for large sets. In [7] the
Cover Trees (CT) are introduced. A CT is composed of levels, where nodes at
some level cover lower levels using a notion of proximity. The SASH data structure [7] is constructed by sampling the database, connecting the non-sampled
elements to their approximate nearest neighbors and repeating in the next level.
In [22], the authors introduced the Rank Cover Trees (RCT) incorporating
ideas from the above algorithms. A tree is built using ordered rank for pruning,
instead of rules derived from distances and the triangle inequality. Node descendants in the tree are obtained using a rank order. Since only rank information is
used for navigation, the number of nodes visited is known a priori. The authors
include proofs, or sketches, for complexity and accuracy bounds.
Malkov et al. [24] introduced the Approximate Proximity Graph (APG), an
index with excellent searching performance. The authors suggest that the speed
and accuracy improvement in APG is the consequence of the incremental construction. They claim it produces a graph with small world properties. APG
can trade speed and accuracy for memory usage. If memory is scarce, speed
and accuracy will also suffer.
Construction and searching are entangled in APG. The construction is incremental and consists of a simple rule: To insert the j-th element, simply
find the (approximate) t-nearest neighbors among the j − 1 elements already
indexed; then, the new item is linked (in both directions) to its t near neighbors. The authors encourage the preservation of links to distant items since
they boost the search performance. Notice that these links naturally arise using
the incremental construction.
Search is also greedy. A starting point is selected and greedily follows the
neighbor minimizing the distance to the query, and repeat until no further
5

Algorithm 1 The search algorithm for APG as described in [25]. We use our
notation N (u) to describe the connected vertices of u, see Section 2 for more
details.
Name: Search algorithm APG
Input: A transition function N , the database S, and the query q, the number
of restarts m
Output: The set res of near neighbors of q.
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:

Let res be an empty min-queue of fixed size k
Let candidates be an empty min-queue
Let visited be an empty set of object identifiers
for i = 1 to m do
Randomly select c ∈ S \ visited
Append c into visited and (d(q, c), c) into candidates and res
loop
Let (rb , best) be the nearest pair in candidates
Remove best from candidates
Let r be the biggest radius in res, empty res defines r = ∞
if rb > r then
break loop
end if
for u ∈ N (best) do
if u 6∈ visited then
Add u to visited and (d(q, u), u) to candidates and res
end if
end for
end loop
end for

improvement is possible. This simple procedure performs poorly. However, it is
likely to boost the recall by performing m searches from random starting points.
Many restarts slow down the overall searching time, yet the authors claim that
using m = O(log n) restarts achieve recall close to one.
The parameters governing accuracy and speed are the number of starting
points m and t, the number of direct neighbors in the graph. Since it is necessary
to restart many times to achieve high accuracy, then the speed is penalized.
Increasing t could also boost the speed at the expense of more memory.
In a follow-up paper, in [25], the searching procedure is updated as follows.
As the previous algorithm, the search is repeated m times, but now the restarts
are persistent. In other words, the closest open node is selected as the new
restarting point. In other words, a set of visited items and a set of candidates
are kept along the entire search process. At the beginning of each restart, the
algorithm appends a random item to the list of candidates to add diversity to
the search process. Algorithm 1 formally describes this algorithm.

6

2. Proximity Search as an Optimization Problem
Combinatorial optimization can be stated as follows. Given a set of states
Ω and a cost function f : Ω → R, find ω ∈ Ω such that for , f (ω) is optimal.
An optimal state is one that reaches either the minimum [or the maximum]. In
other words, the goal is to obtain ω such that
f (ω) = min f (x) [or, max f (x) in the case of maximization].
x∈Ω

x∈Ω

For our indexing problem we set Ω = S, and for a fixed query q ∈ U we set
f (ω) = d(w, q), and we will seek for the minimum of this function.
It should be clear that f (nn(q)) is minimal, hence optimizing f (·) implies
finding the nearest neihbor of q.
For the optimization problem, we will preserve the following restrictions,
derived from the metric space approach.
• The distance function is the only way to compare any pair of items.
• Since objects are seen as a black box, items (states) are given, and cannot
be generated (for example from averaging).
• The gradient is not defined, i.e., the notion of direction is not well-defined.
• There are no assumptions about the distance function d. e.g. uniqueness
nor convexity are assumed.
• The size of the database can be large, and the intrinsic dimension is high.
The above restrictions narrow the alternatives in the selection of optimization techniques. We cannot use gradient or steepest descent. There are, however, many options based on local search. Those techniques require some binding
to navigate between states (which are objects of the database in our case). In
other words, we need a graph with vertex the objects of the database and with a
suitable definition for the edges. Edge definition is critical for the performance
of the optimization algorithms. The three main components of our approach are
the transition function N , the function to be minimized f , and the minimizing
algorithm A.
Definition 1 (State transitions N ). N : Ω → Ω+ is a function associating
each state u with a set of possible transitions.
Definition 2 (The function to be optimized dq ). Each query q defines a
function to minimize, dq (u) = d(q, u) for u ∈ Ω.
Definition 3 (The minimizing algorithm A). A is an algorithm to minimize dq . It explores Ω using the transitions N to find a state minimizing dq for
a given q.

7

2.1. About A
Combinatorial optimization is a well-established area of computer science.
There are many alternatives to choose from. Many meta-heuristics are discussed for example in citeBurkeKendall2014. We selected local search due to
the properties of our problem statement and restrictions.

(u)
d(u, q)
q

u

d(u’, q) = min {d(v, q) | v ∈

(u)}

u’

Figure 1: An object u and its neighborhood N (u), u0 is the item closer to q in the neighborhood.

May we remind that in the problem (S, d) is a metric database, and N
describes the edges of a graph with the vertices in S. With a proper graph, it is
possible to minimize dq starting from any random vertex u ∈ S. The property
of the edges is to connect a vertex with neighbors locally minimizing the target
cost function at each step. There should be at least one v ∈ N (u) minimizing
dq . Figure 1 illustrates a step of the minimization procedure. At each point u,
we must explore its neighborhood N (u), and try to minimize dq . For example,
let u0 ∈ N (u) be an object such that d(q, u0 ) = min{d(q, v) | v ∈ N (u)∪{u}}, we
move to u0 and repeat the process, until we reach a local minimum where there is
no object closer to q than the current node. While typical metric access methods
will try to discard objects using the triangle inequality, see §1, the optimization
approach will seek to minimize dq . In this process, we will find many local
minima before we reach the global minimum or a good approximation.
We review some prominent meta-heuristics that allow us to use them as part
of our algorithms.
2.1.1. Greedy Search.
Perhaps the simpler working algorithm for local search is the Greedy search
algorithm, it works as follows to solve a nearest neighbor query (nn(q)):
i) Start in random state from Ω.

8

ii) Let s be the current best answer to nn(q). Compute d(q, u) for each u ∈
N (s).
iii) If there exists some u that decreases dq then u becomes the best answer
and repeat step ii). Otherwise stop the search.
iv) Our best guess of nn(q) is the item producing the minimum known value
of dq .
The greedy search will retrieve a local minimum for dq . The quality of the
answer can be improved by using multiple restarts of greedy search; this is a
standard amplification method.

m
800
600
400
200
0

P
P
P
P

m
40
35
30
25
20
15
10
5

=0.99
=0.95
=0.90
=0.80

0.02 0.04 0.06 0.08 0.1

P

P
P
P
P
0.1

(a) P ranges from 0.005 to 0.1

0.2

0.3

=0.99
=0.95
=0.90
=0.80
0.4

0.5

P

(b) P ranges from 0.1 to 0.5

Figure 2: The number of steps m for several desired recall (P ? ) as a function of the expected
recall per step P . Observe how low P values need extremely large m values. This behavior is
consequent with the curse of dimension.

If each search has probability P of finding the exact nn(q), then we can
achieve the desired probability P ? with m independent searches with
m=

log (1 − P ? )
;
log (1 − P )

(3)

this expression supposes independence of N (·), which can be hard to achieve.
Notice that this formula does not predict the closeness of the retrieved item at
any given m; however, it is faster than an exact method solving nn(a).
From the optimization perspective, a low P implies a high number of local
minimums. This issue is probably yet another incarnation of the curse of dimension in the proximity search problem. Figure 2 illustrates the effect. Notice
that our method is practical when
cost = m × hops × |N (·)| + 1 < n,

(4)

where hops is the expected length of the path (number of transitions) from a
random starting point to the local minimum found, and |N (·)| the expected
number of transitions per item. We show in our experimental section that the
number of transitions is low in practice.
9

2.1.2. Tabu Search.
Tabu search is a meta-heuristic allowing to take non-greedy decisions at each
step, i.e., not always decreasing dq at each transition. This strategy recognizes
the existence of local minima and plateaus; the idea is to avoid loops and to
follow the same path twice. Tabu search solves the problem adding memory to
mark already visited states as tabu, and then ignoring them on later stages of the
search process. We will store visited states in short-term, medium-term, or longterm memory. The desired behavior, along with the problem, and the available
resources, determine the kind of memory used in a given implementation of tabu
search.
2.1.3. Beam Search.
Beam search is a strategy that optimizes dq using a beam of states B ⊂ Ω of
size b. At each step, beam search computes B 0 = ∪s∈B N (s) and selects the best
b from B and B 0 . The search stops whenever is impossible to improve B. Tabu
search can be used along with beam search to avoid non-essential computations.
At the beginning of the search process, B is populated randomly from Ω. Notice
that the low memory requirement is the main advantage especially for large
datasets or systems with limited memory.
2.2. About N
A fundamental component of the local search meta-heuristic is the transition
function N , if properly defined, it can ensure the convergence to the global minimum or at least to a suitable local minimum. To give an idea of the complexity
of setting a good N , lets think in A be the Greedy Search algorithm. The
function N defines a graph where the vertices are the objects of the database,
and a vertex u is connected with the elements of N (u). Note that the greedy
search starting at a point u will always stay in the path connected component
of u. Therefore, a natural requirement for the graph is to be path connected.
Since our approach is approximate, then our transition function does not
need to ensure that A finds the global optimum dq . From the similarity search
perspective, we allow to lose some results, and also, some retrieved items will
not be part of the exact query. Nonetheless, in practice, the error rate will be
low, and the performance gain is significant. Those claims are supported with
experimentally.
The proposal of Malkov et al. [24, 25] introduce a working N construction
for its APG, already explained in §1. The underlying graph N can be used with
improved versions of greedy search and best-first search. Our contribution is
more about A. As the APG, we will use the undirected version of the graph.
3. Direct Improvements on APG
As detailed on Section 1, search in APG is repeated m times to improve
the expected recall. Authors suggest in [25], that m must be adjusted for each
dataset. This requirement supposes an additional parameter selection, which

10

Algorithm 2 The search algorithm for APG*. The main difference is that m
is not longer necessary, yet a small value σ is needed to guess a good number
of restarts. The main idea is to apply the search procedure repeatedly while
the covering radius of res improves. In practice, it produces a competitive local
minimum.
Name: Search algorithm APG*
Input: A transition function N , the database Ω = S, the query k − nn(q), a
minimum number of repeats σ (pretty small, e.g., 4)
Output: The set res of near neighbor of q.
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:

Let res be an empty min-queue of fixed size k
Let candidates be an empty min-queue
Let visited be an empty set of object identifiers
repeat
Let cov∗ ← cov(res)
for i = 1 to σ do
{This for loop is identical to that of APG (Algorithm 1)}
Randomly select c ∈ S \ visited
Append c into visited and (d(q, c), c) into candidates and res
loop
Let (rb , best) be the nearest pair in candidates
Remove best from candidates
if rb > cov(res) then
break loop
end if
for u ∈ N (best) do
if u 6∈ visited then
Add u to visited and (d(q, u), u) to candidates and res
end if
end for
end loop
end for
until cov∗ = cov(res) {Stops when there is no improvement over dq , i.e.,
cov(res)}

can be difficult because the best m ideally should be adapted as n grows. A
significant overestimation of m impact the construction time because it is incremental, based on searches. Searching is also negatively affected because of the
number of restarts.
Before we introduce our new algorithms, we will define an important function. Let, result be a min priority queue (min-queue) of fixed size k. The
cov(result) function is defined as follows: if result is under its full capacity
then it evaluates to the maximum possible distance value; if it is full, then
cov(result) is the radius of the furthest item.

11

Our first index contribution is Algorithm 2 (APG*). The objective is to
remove parameter m, checking the changes over the global result set, via the
covering radius cov(res). In some sense, we are minimizing dq , and stopping
the algorithm when our guess of dq cannot be improved after σ tries. All our
new algorithms will use the same technique. Notice that σ should be small (in
the range of 2 to 4). The idea is to relax the stopping condition, avoiding the
introduction of an extra parameter. Larger values of σ should be avoided in
favor of other parameters.
Algorithm 3 The search algorithm for APG*-R. As APG*, it tries to adjust
the number of restarts automatically.
Name: Search algorithm APG*-R
Input: A transition function N , the database Ω = S, the query k − nn(q), a
minimum number of repeats σ (pretty small, e.g., 4),
Output: The set res of near neighbor of q.
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:

Let res be an empty min-queue of fixed size k
Let visited ← ∅
repeat
Let cov∗ ← cov(res)
for i = 1 to σ do
Let res∗ be an empty min-queue of fixed size k
Randomly select sp ∈ S \ visited
Append sp into visited and (d(q, sp), sp) into res∗
Let evaluated ← ∅
repeat
for v ∈ N (sp) do
if v 6∈ evaluated then
Add (d(q, v), v) into res∗
end if
end for
Update sp to be the closer item in res∗ not in visited. If there is no
item satisfying the condition, then sp is set as undefined.
Add sp into visited
until sp is undefined
res ← res ∪ res∗ {The properties of the min-queue of fixed size should
be preserved}
end for
until cov∗ = cov(res) {Stops when dq does not improve, i.e., cov(res) stays
unchanged}

Our second contributed index is named APG*-R (Algorithm 3). The idea is
that APG and APG* can get stuck on a local minimum since they are limited
to use always the best-known candidate. APG*-R share structural similarities
with APG and APG*; however, there is a conceptual difference. On the one
12

hand, APG*-R uses random starting points, not the best known in APG, and
also, res∗ is local to each step. We also added a local memory evaluated to
avoid additional evaluations, this scheme is different to that found whenever
the candidates set has a global scope, and so, the improvement policy is also
global.
Algorithm 4 Beam Search over Ω and N .
Name: Beam Search
Input: A transition function N , the database S, the query q, the size of the
beam b
Output: The set res of near neighbor of q.
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:

Let res be an empty min-queue of fixed size k
Let beam be an empty min-queue of fixed size b
Set visited ← {}
Set beam ← {}
for i = 1 to b do
Randomly select u from S
Add u to visited and (d(u, q), u) to res and beam
end for
repeat
Let cov∗ ← cov(res)
for i = 1 to σ do
Let beam∗ ← {} {Fixed sized priority min-queue of size b as beam}
for c ∈ beam do
for u ∈ N (c) do
if u 6∈ visited then
Add u to visited and (d(q, u), u) to res and beam∗
end if
end for
end for
beam ← beam∗
end for
until cov∗ = cov(res) {Stops when there is no improvement over dq , i.e.,
cov(res)}

Our final index is BS (Algorithm 4), an index using beam search over S
and N . As in our previous approaches, we ensure at least σ steps to take the
decision to stop the search. So, if after some minimization steps, dq does not
change, then we achieved a (good) local minimum. The critical parameter is
the size of the beam b, which can be adjusted at any time. Notice that beam
search does not need to restart; however, it needs to know b. As we will see
in our experimental section, the use of beam search as an alternative to other
schemes produce significant improvements in both recall and search speed.

13

dimension

3 × n5

106

3 × 106

16
32
64

0.035s
0.050s
0.076s

0.137s
0.202s
0.342s

0.320s
0.438s
0.708s

Table 1: Time spent on a sequential scan.

4. Experimental Results
In this section, we compare the experimental performance of our technique
against APG, KNR, and the RCT. We fixed our attention in the recall and
search speed. All experiments were performed on a 24-core Intel Xeon 2.60
GHz workstation with 256GB of RAM, running CentOS 7. We do not use
the multiprocessing capabilities in the search process. Both, the index and the
database were maintained in memory.
Our benchmarks consist of both synthetic and real databases. In particular,
we selected the following databases:
— Nasa. This database is a collection of 40, 150 vectors of 20 coordinates obtained from the SISAP project (http://www.sisap.org). It uses L2 as
distance function. A sequential search completes on 0.005 seconds on our
testing machine.
— Colors. The second benchmark is a set of 112, 682 color histograms (112dimensional vectors) from SISAP, under the L2 distance. A sequential
search uses 0.053 seconds.
— Wiktionary. A third benchmark is an English dictionary with 736, 639 entries
with Levenshtein’s distance as metric.1 A sequential search completes in
0.700 seconds.
— Documents. This database is a collection of 25,157 short news articles from
the TREC-3 collection of the Wall Street Journal 1987-1989. We use their
tf-idf vectors, taken from the SISAP project [18], www.sisap.org. We
use the angle between the vectors as the distance measure [5]. We remove
100 random documents from the collection and use them as queries (thus
these 100 documents are not indexed in the database). The objects are
vectors of hundred thousands coordinates. This dataset has a very high
intrinsic dimension in the sense of Chavez et al. [13]. Most exact metric
indexes will need to evaluate the entire dataset to solve nearest neighbor
queries. As a reference, a sequential scan needs 0.185 seconds.
1 Available

at https://www.wiktionary.org/

14

— BigANN-1M. A one million subset of the TEXMEX corpus.2 The original
dataset contains a billion SIFT points (128 coordinates). We also extract
256 queries from the standard query set (not in the dataset). In average,
a sequential search of the nn under L2 needs 0.441 seconds to be solved.
— RVEC. On the one hand, for measuring the effect of varying dimension, we
used datasets of 16, 32, 64, 128 and 256 dimensions, each one with 105
randomly generated items. A query is solved by exhaustive search in
0.015, 0.014, 0.024, 0.044, and 0.081 seconds, respectively. On the other
hand, the performance as the size increases is measured using 16, 32 and
64-dimensional datasets of size 3×105 , 106 , and 3×106 items each. Table 4
shows the search time spent by an exhaustive evaluation of the dataset. All
these datasets use the L2 distance. The importance of synthetic data with
a uniform distribution is to study the the performance of the algorithms
for a fixed size and dimensionality.
To compare against LSH we used the E2 LSH tool3 by Alex Andoni. It automatically optimizes an LSH index for the available memory. We fix the index
memory to be 1, 3, and 5GB, which is many times higher than the necessary
by other indexes. To be fair, since LSH offers distance-based rather than recallbased guarantees, we used it to search by radius (i.e., return any element within
distance r to the query), setting r as the average radius of the 30-th nearest
neighbor. To avoid extreme or uninteresting cases, we removed queries yielding empty answers and those with more than 1000 answers. For the RCT [22]
we used the C++ implementation, kindly shared by the authors. Notice that
LSH, RCT, and FLANN do not offer an implementation for Levenshtein and
the cosine similarity, so they only are compared on datasets with the L2 metric.
Therefore, we include LSH and RCT only on experiments over vector datasets
fulfilling this requirement. Also, LSH is not present in RVEC benchmarks since
its performance is not preserved under this type of datasets [23, 4].
While our implementations are written in C# and run under the Mono
virtual machine, E2 LSH, RCT, and FLANN are written in C++, and run in
native form. For example, their sequential search has an advantage of around
2 to 3X in the case of E2 LSH and RCT, and 5X times in the case of FLANN
due to their code’s optimizations. To give a full picture of the performance,
we report both the absolute times and speedup. The latter is computed on the
sequential scan in its implementation. We also give the actual times to give the
full context to the reader.
4.1. Performance per Database
In the first batch of experiments, we present the indexes APG*, APG*-R,
BS, with different values for t (the number of neighbors in the construction) for
the Colors, Nasa, Wikionary, Documents, and BigANN-1M databases.
2 Available

at http://http://corpus-texmex.irisa.fr/
~andoni/LSH/

3 http://www.mit.edu/

15

The APG, and BS are labeled with t = 8, 16, and 32, for each value of t,
the parameter m for APG and b for BS varies from 8, 16, and 32. The APG*
and APG*-R takes t = 8, 16, and 32. Each KNR uses k = 7, 12 near references
among 2048 total references (larger values of k will produce better recalls but
at the cost of higher memory consumption). For each k, the index is allowed
to review at most 3% and 10% of the database (these will be the two points in
each curve). For the RCT we used a covering factor of 32, 64, and 128 (larger
values will produce better recalls but slower searches).
On Figure 3(a) we show the recall and speedup of the methods for the Colors
database. The APG’s for t = 8 are the fastest but with less recall, as we increase
the number of neighbors in the construction, the recall gets better, but the speed
is affected. When we compare this with the line of the APG*, the APG* is very
close to the corresponding APG. Now, the APG*-R achieves a higher recall
but at the cost of losing speed. Next, we have the BS, it has the best overall
performance, that is, for a given recall, it has the best speed. Note how it gets
the most recall. The RCT, KNR, and LSH are left behind.
The results of the Nasa database are shown on Figure 3(b). This figure
corresponds to the smallest database among our benchmarks, and all the indexes
reach competitive recall. The results for the APG are similar to the Colors
database, and the recall increases with the amount of memory; however, we can
find faster alternatives. The APG* is still very close the corresponding APG.
Note that APG*-R achieve a significant improvement as compared with APG
and APG*. All the BS perform with high recall at a competitive search speed.
The RCT has more appeal here than on Colors being the faster index. While
KNR is quite fast, it gets stuck and has difficulties to improve its recall just
adjusting the k parameter; this behavior is shared with LSH since their recall
values are below others.
The Wiktionary works with a costly distance function, in this type of datasets
it is critical to lower the number of computed distances. On Figure 3(c) almost
every index is hundreds of times faster than a sequential scan. The performance
of the APG is again below the others and the APG* is again essentially the
same. On the other hand, the APG*-R is above all others. For the correct
amount of neighbors, the APG*-R is an absolute improvement of the APG*.
The BS excels in performance, and it reaches the best recall. KNR achieve its
maximum performance rapidly, and it is left behind for all local search methods.
The Figure 3(d) shows the Documents database. Note how most of the
indexes have the same curve of tendency. Remember that this database has
high intrinsic dimension, so, we can have a high recall or a high search speed,
but it is quite difficult to be good on both. The APG* is practically the same as
the APG, however without hassle about the m parameter. KNR indexes exhibit
a good trade between recall and search speed, but APG*-R reaches better recall
values. As before, BS shows the best performances among the alternatives.
The results for the BigANN-1M are shown on Figure 3(e). The results are
striking similar to the Documents experiment. The principal difference is the
higher values for the speedup because BigANN-1M is much larger than Documents. Apart from that, the APG are the same as the APG*; the APG*-R
16

have higher recalls values; the BS instances have higher recalls at the cost of
the speedup; the KNR, RCT and the LSH have below the average speed.
In summary, APG*-R is better in general than the APG* which is similar
to APG but computing m online. Also, all the BS have good recall and speed,
making it the rule of thumb for application needing high quality results.
4.2. The effect of the dimension on the Search Performance
We compared indexes to see how dimensionality of datasets affects them.
Here, we selected RVEC datasets with a wide range of dimensions, specifically,
16, 32, 64, 128, and 256. We measured recall, the search cost, and the speed up.
Each line of points represents the results for each dimension.
The Figure 4 shows the cost of the indexes. We compared the recall with the
cost of the searches. Figure 4(a) shows the indexes using 1 to 10 integers per
element of the database, that is, on the APG the number of neighbors is small.
The APG and the APG* are similar and produce low recall. This performance
contrasts with the APG*-R that shows a significant improvement on it. These
results demonstrate the effectiveness of this variant. Moreover, the BS has a
broad range of action, note how the recall is much better than the APG at
the cost of computing more distances. Here, the KNR has a good performance
with one of the configurations reaching the higher recall and showing an almost
constant cost on the dimension. The RCT, like the BS, is all spread over but
with more distances computed.
Figure 5(a) shows the speedup. The APG and APG* are very fast because
of the small number of distances they have to compute but the low recall makes
them useless. The APG*-R appears in the middle of the other indexes. The
different instances of the BS show how it gains recall but loses speed when we
increase the beam. The KNR gets good recalls but is slower than the average.
Finally, the RCT is fast and with good recall for these smaller indexes. Note
how the recall for the APG’s and BS increases when adding more memory.
For Figure 4(b) we let the methods take more neighbors for the elements;
they compute more distances than smaller indexes, but they reach higher recall
values. The APG and APG* cannot get a good recall but the improve of the
APG*-R makes them comparable to the BS. Figure 5(b) contains the speeds
of those indexes. The PG and APG* are very fast, although with a low recall.
The APG*-R and BS with b = 8 are very similar, one a little faster and the
other with better recall. As expected, the other configurations of BS have higher
recalls.
Figures 4(c) and 4(d) are very similar; they show how the APG and APG*
become to get a good recall for the first dimensions but continue to have a low
value for the higher ones. With the APG*-R and BS, we see that they can have
very high recalls even for the big dimensions.
Figures 5(c) and 5(d) are again very similar. APG and APG* are not very
fast. APG*-R is faster for dimensions 16, 32, and 64 but the BS with b = 8
becomes faster for the larger dimensions.

17

90.0

60.0
50.0
40.0
30.0

16.0
14.0
12.0
10.0
8.0
6.0

20.0

4.0

10.0
0.0

APG t=8
APG t=16
APG t=32
APG*
APG*-R
BS t=8
BS t=16
BS t=32
KNR-Cos k=7
KNR-Cos k=12
RCT
LSH
FLANN

18.0

speedup

70.0

speedup

20.0

APG t=8
APG t=16
APG t=32
APG*
APG*-R
BS t=8
BS t=16
BS t=32
KNR-Cos k=7
KNR-Cos k=12
RCT
LSH
FLANN

80.0

2.0

0.60.7

0.8

0.9

0.95

0.0

0.96

0.60.7

0.8

0.9

0.95 0.96 0.97

recall

(a) Colors

0.99

35.0

APG t=8
APG t=16
APG t=32
APG*
APG*-R
BS t=8
BS t=16
BS t=32
KNR-Cos k=7
KNR-Cos k=12

300.0
250.0
200.0

APG t=8
APG t=16
APG t=32
APG*
APG*-R
BS t=8
BS t=16
BS t=32
KNR-Cos k=7
KNR-Cos k=12

30.0
25.0

speedup

350.0

20.0
15.0

150.0
10.0

100.0

5.0

50.0

0.9

0.95 0.96

0.97

0.98

0.0

0.99

0.60.7

0.8

0.9

0.95 0.96 0.97

recall

0.98

recall

(c) Wiktionary

(d) Documents

250.0

APG t=8
APG t=16
APG t=32
APG*
APG*-R
BS t=8
BS t=16
BS t=32
KNR-Cos k=7
KNR-Cos k=12
RCT
LSH
FLANN

200.0

speedup

0.0

0.999

(b) Nasa

400.0

speedup

0.98

recall

150.0

100.0

50.0

0.0

0.60.7

0.8

0.9

0.95 0.96 0.97

0.98

0.99

0.999

recall

(e) BigANN-1M

Figure 3: Recall and speedup performance for our real world datasets.

18

0.99

0.18

0.12
0.10
0.08
0.06

APG m=8
APG m=16
APG m=32
APG*
APG-R*
BS b=8
BS b=16
BS b=32

0.25

0.20
# distances / n

0.14

# distances / n

0.30

APG m=8
APG m=16
APG m=32
APG*
APG-R*
BS b=8
BS b=16
BS b=32
KNR-Cos k=7
KNR-Cos k=12
RCT c=32
RCT c=64
RCT c=128

0.16

0.15

0.10

0.04
0.05
0.02
0.00

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.00

1

0.6 0.7

0.8

recall

(a) 1 to 10 integers

0.95

0.97

0.99

(b) 10 to 30 integers

0.50

0.60

APG m=8
APG m=16
APG m=32
APG*
APG-R*
BS b=8
BS b=16
BS b=32

0.40
0.35
0.30

APG m=8
APG m=16
APG m=32
APG*
APG-R*
BS b=8
BS b=16
BS b=32

0.50

0.40
# distances / n

0.45

# distances / n

0.9
recall

0.25
0.20

0.30

0.20

0.15
0.10

0.10

0.05
0.00

0.6 0.7

0.8

0.9

0.95

0.97

0.00

0.990.999

recall

0.6 0.7

0.8

0.9

0.95

0.97

0.990.999

recall

(c) 30 to 100

(d) 100 to 300 integers

Figure 4: Performance comparison (recall vs search cost) among our indexes and several state
of the art techniques for fixed n = 105 over several dimensions. Each curve corresponds
contains points that characterize the performance in a dimension, e.g., 16, 32, 64, 128, and
256; as a hint to review figures, large dimensions correspond to smaller recalls and scale is
exponential.

19

30

20

15

10

10
8
6
4

5

0

APG m=8
APG m=16
APG m=32
APG*
APG-R*
BS b=8
BS b=16
BS b=32

12

speedup

25

speedup

14

APG m=8
APG m=16
APG m=32
APG*
APG-R*
BS b=8
BS b=16
BS b=32
KNR-Cos k=7
KNR-Cos k=12
RCT c=32
RCT c=64
RCT c=128

2

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0

1

0.6 0.7

0.8

recall

(a) 1 to 10 integers
9

0.95

0.97

0.99

(b) 10 to 30 integers
7

APG m=8
APG m=16
APG m=32
APG*
APG-R*
BS b=8
BS b=16
BS b=32

7
6

APG m=8
APG m=16
APG m=32
APG*
APG-R*
BS b=8
BS b=16
BS b=32

6
5

speedup

8

speedup

0.9
recall

5

4
3

4
2

3

1

2
1

0.6 0.7

0.8

0.9

0.95

0.97

0

0.990.999

recall

0.6

0.7

0.8

0.9

0.95

0.97

0.99 0.999

recall

(c) 30 to 100 integers

(d) 100 to 300 integers

Figure 5: Recall and speedup comparison in four different classes of memory among our
indexes and several state of the art techniques for fixed n = 105 over several dimensions.
Each point in curves corresponds to a different dimension 16, 32, 64, 128, and 256; as a hint to
review figures, large dimensions correspond to smaller recalls and scale is exponential.

20

4.3. Scalability
A central property of any search index is its capacity to scale with the size of
the database. Even a static index can be dynamized using the Bentley and Saxe
[6] approach, as it has been shown for metric indexes in [29]. This dynamization
comes with a small penalization at search time, in addition to a constant factor introduced in the construction. This technique is beyond the scope of this
document, and static indexes are a single large index without additional search
overheads. Notice that the APG’s and our BS has an incremental construction.
Hence they can handle a growing database naturally. In any case, most indexes
will be affected by database growth, and we expect the performance to be different in this setup. We used the randomly generated database of dimensions
16, 32, and 64, each with sizes 3 × 105 , 106 , and 3 × 106 . Each point in a curve
is produced by a different database size; speedups are larger as n increases. As
the database grows, so does the size of the indexes and it could not fit in main
memory. Because of that, we set the number of initial neighbors to 8 for the
APG’s and BS.
Figure 6(a) shows the recall and the search speed for the different sizes of
the database for 16-dimension. The APG and APG* scored a recall below 0.7
on all the sizes, even when using 32 searches for the query. As shown in the
previous experiments, the APG*-R gets a boost on the recall compared to the
APG*. Note how the lines of the BS are practically vertical meaning that the
index gains in speed much more than it loses in the recall. The RCT does not
get a recall as good as the BS and is considerably slower. On this and the other
dimensions, the sizes corresponding to the lines of each method are 3 × 105 , 106 ,
and 3 × 106 from right to left except for the KNR. For this method, the recall
gets better as the size of the database increases. In Figure 6(b) we see the recall
and the number of distances for the queries. This general behavior is the same
in the next figures.
Figure 6(c) shows the same tendency: the APG*-R has been an improvement
over the APG’s, the BS having very competitive recalls and speeds, the KNR
been slow, and the RCT behind the BS. The Figure 6(d) shows the same but
with the distances computed. The same can be said to the figures 6(e) and
6(f), but the interesting part is to note how the change of dimension affects
more the indexes than the size of the database for reaching high recalls. For
16 dimensions, high recall values can be achieved for all the sizes; for 32 and
64-dimension, the BS always gets better recalls.
As a summary, we arrange the results for the methods with best recall values
on Tables 2, 3, and 4 for 16, 32, and 64 dimensions respectively.
5. Tuning APG*, APG*-R, and BS
The flexibility of the combinatorial optimization framework for proximity
search is one of the sources of their great potential, but also may be problematic for an end-user. In this section, we provide a simple guide for finding a
competitive tuning, although it should be noticed that the optimal selection of
the parameters cannot be guaranteed.
21

Table 2: The achieved performance for several sizes of the random database of dimension 16.
Index
APG*-R
BS b=8
BS b=16
BS b=32
KNR-Cos k=7
KNR-Cos k=12
RCT

105
Cost Recall
1868.84
0.95
2063.57
0.96
3137.97
0.98
4937.93
0.99
4553.08
0.85
4950.76
0.92

3 × 105
Cost Recall
3561.40
0.97
2297.48
0.95
3566.05
0.98
5621.50
0.99
9998.99
0.92
10948.57
0.96
15836.4
0.95

106
Cost
4002.21
2473.05
3925.69
6200.69
28686.36
31941.04
26863.9

Recall
0.96
0.95
0.98
0.99
0.95
0.98
0.94

3 × 106
Cost Recall
4208.77
0.96
2705.92
0.94
4217.09
0.97
6637.78
0.98
83004.65
0.97
91945.43
0.99
42947
0.94

Table 3: The achieved performance for several sizes of the random database of dimension 32.
Index
APG*-R
BS b=8
BS b=16
BS b=32
KNR-Cos k=12
RCT

105
Cost Recall
4484.40
0.86
4405.64
0.87
6565.40
0.94
9450.75
0.97
5105.69
0.65

3 × 105
Cost Recall
9276.67
0.88
5529.78
0.83
8105.05
0.91
12122.14
0.95
11100.57
0.71
19229.8
0.71

106
Cost
12865.79
6480.63
10078.76
15067.28
32088.94
34085.0

Recall
0.82
0.78
0.88
0.94
0.78
0.67

3 × 106
Cost Recall
14168.41
0.82
7309.23
0.73
11446.03
0.84
17256.64
0.91
92088.18
0.83
52872.1
0.63

Table 4: The achieved performance for several sizes of the random database of dimension 64.

Index
APG*-R
BS b=8
BS b=16
BS b=32

105
Cost
10357.14
8181.72
11692.90
16510.76

Recall
0.69
0.71
0.84
0.91

3 × 105
Cost Recall
21958.30
0.74
11582.28
0.62
16907.20
0.76
24060.58
0.86

22

106
Cost
35579.05
16384.80
24275.12
35590.51

Recall
0.61
0.50
0.66
0.79

3 × 106
Cost Recall
51994.29
0.57
22398.66
0.40
32656.85
0.56
49443.82
0.71

180

120
100
80

APG m=8
APG m=16
APG m=32
APG*
APG* R
BS b=8
BS b=16
BS b=32
KNR-Cos k=7
KNR-Cos k=12
RCT c=128

100000

# distances / n

140

speedup

1e+06

APG m=8
APG m=16
APG m=32
APG*
APG* R
BS b=8
BS b=16
BS b=32
KNR-Cos k=7
KNR-Cos k=12
RCT c=128

160

10000

60
1000

40
20
0

0.60.7

0.8

0.9

0.95

0.97

0.98

100

0.99

0.6
0.7

0.8

0.9

recall

(a) Search speedup / 16 dimensions
120

0.97

0.98

0.99

(b) Search cost / 16 dimensions
1e+06

APG m=8
APG m=16
APG m=32
APG*
APG* R
BS b=8
BS b=16
BS b=32
KNR-Cos k=7
KNR-Cos k=12
RCT c=128

80

60

APG m=8
APG m=16
APG m=32
APG*
APG* R
BS b=8
BS b=16
BS b=32
KNR-Cos k=7
KNR-Cos k=12
RCT c=128

100000
# distances / n

100

speedup

0.95
recall

40

10000

20

0

0.60.7

0.8

0.9

1000

0.95

0.60.7

0.8

recall

(c) Search speedup / 32 dimensions
70

1e+06

40
30

APG m=8
APG m=16
APG m=32
APG*
APG* R
BS b=8
BS b=16
BS b=32
KNR-Cos k=7
KNR-Cos k=12
RCT c=128

100000
# distances / n

speedup

50

0.95

(d) Search cost / 32 dimensions

APG m=8
APG m=16
APG m=32
APG*
APG* R
BS b=8
BS b=16
BS b=32
KNR-Cos k=7
KNR-Cos k=12
RCT c=128

60

0.9
recall

10000

20
10
0

0.6

0.7

1000

0.8
recall

0.6

0.7

0.8
recall

(e) Search speedup / 64 dimensions

(f) Search cost / 64 dimensions

Figure 6: Speedup and search cost as a function of the expected recall of our indexes and
state of the art. We use datasets with dimensions 16, 32, and 64. Each dataset has three
instances with sizes of 3 × 105 , 106 , and 3 × 106 . We also limited memory to be under eight
direct neighbors per item, i.e. 16 undirected neighbors in average.

23

First of all, a guide for index selection. We first notice that APG* performs
as good as APG; however, APG*-R, and BS perform better. Both APG*-R
and BS need to know the number of direct near-neighbors t. APG*-R has no
additional parameters since it adjusts the number of repetitions based on a
simple iterative optimization; then the robustness is entirely based on t. For
this reason, t should be relatively large; our suggestion is to use values greater
than 16 but smaller than 64. Since the underlying graph is undirected, the
memory usage will be n(2t + 1) integers. Thus, the right selection of t should
consider this amount of main memory for the adequate performance. For the
case of BS, t can be smaller, i.e., 8 or 16; larger values of t should be used for
very high dimensions, as described in our experimental section. After fixing t,
the free parameter b (the size of the beam) should be adjusted to obtain the
desired performance. A significant b will produce high-quality search results
while increasing search time. In contrast, faster searches will be achieved for
small b values, yet the quality of the results could decrease. As a rule of thumb,
we recommend fixing b to the largest value producing a search time below the
user threshold. Nonetheless, good recall values are achieved for most of our
datasets for b values between 16 and 64. A more precise setup of b could be
accomplished with the help of a query set Q. The idea is to adjust b to obtain
a minimum recall on Q; this fitting should be applied after a fixed number of
insertions.
6. Summary and Perspectives
In this manuscript, we approach proximity search as a combinatorial optimization problem. To this end, we define the necessary framework to allow
the use of combinatorial optimization algorithms studied in the literature. In
particular, this contribution focuses on local search algorithms.
Using our new perspective, we introduced three near-neighbor indexes, namely,
APG*, APG*-R, and BS. Those indexes use an underlying graph and several
variants of local search to navigate them. As part of our contribution, we compare our techniques against state of the art indexes. We found that our indexes
significantly surpass the performance of the majority of the alternatives, in almost all of our benchmarks. In particular, we believe that our indexes are of
practical use due to their excellent trade-off among speed, recall, and index size.
In addition to our core contributions, this manuscript also gives an extensive
comparison of state of the art on the nearest neighbor problem. We tested
algorithms working in both general metrics, and those specialized on vector
spaces, in particular for the Euclidean distance. We also include E2 LSH, a well
known and used near neighbor index based on hashing. Our intention was to
provide a better and wider picture of the alternatives for near neighbor search,
disregarding the underlying theoretical basis with a comprehensive guide to
real-world applications in mind.
In the present work, we fixed the construction of N to be that defined by
APG while we focus on finding better algorithms for navigating the graph, i.e.,
A. However, new N definitions with improved or additional properties should
24

be explored. The indexes introduced in this manuscript have some parameters
to be determined experimentally. This situation is not uncommon; however,
this complicates the usage of the new indexes in real world tasks. Even when
we give guidelines for tuning our indexes, it is necessary to look forward to
construction algorithms that determine their parameters automatically using,
for example, the topology and geometry of the given dataset.
To be effectively used in production systems, our indexes should support
some degree of persistence. A secondary memory implementation could be a
solution; however, the high number of random seeks arising on graph-based
algorithms is a problem in disk-based storage. A possible solution is to use
exclusively solid state storage as secondary memory; however, I/O operations
are still far slower than main memory. An alternative scheme is a pure memory
index having a secondary memory mirror, such that any modification is performed in both memory and disk storage. Therefore, the idea is that the data
persist and the index become available as soon the process is launched in the
system. This scheme is quite attractive as it is implemented in high-performance
databases like Redis [9].4
Another important task is to support dynamic operations. This feature is
useful for applications using metric databases to represent knowledge, like those
using incremental learning. Since our indexes are incrementally constructed, the
insertion operation is native. However, the deletion algorithm is not yet studied.
The literature in metric indexes used to solve the deletion operation marking
those items as unavailable for most operations [13]. While this is enough for
workloads with just a few deletions, a high rate of removals could yield to a
significant unnecessary extra memory. In other approaches, the metric indexes
can afford real deletions. However, to maintain the invariants of the indexes, a
partial or even a complete reconstruction of the index is needed ([31, 32, 14]).
A recent approach uses the strategy of Bentley and Saxe ([6, 30]) to produce
dynamic structures using a list of log n static structures. This approach gives
the possibility of using any static indexes at the expense of a log n factor in
the search cost. Also, since only amortized times are ensured, incremental
constructions cannot be used (i.e. a single insertion or deletion operation can
be lengthy). We believe that our combinatorial optimization approach is robust
enough to support deletions, and the primary task becomes on determining
optimal replacements of the removed items. This operation deserves attention
in future research on the area.
The current construction of APG and our indexes is a sequential process,
i.e., the i-th item is connected with its neighbors among the previous i − 1
items. In the case of APG and APG*, each search is composed of several greedy
searches. Each greedy search is not independent of others since they use the
best-known candidate at each step. A possible point of parallelization is at each
hop since 2t distances are computed. In the case of BS, at each hop 2tb distance
4 Redis is an in-memory data structure store, used as a database, cache and message broker,
https://redis.io/

25

computations are performed. In any case, the necessary distance calculations
can be executed in parallel. However, the intrinsic cost of thread creation and
management can be too high to take advantage of parallelism. APG*-R, on
the other hand, can be improved with independent searches. Nonetheless, there
exists space for clever and better parallel construction algorithms.
Finally, we also should mention the lack of theoretical guarantees for both the
search time and the expected recall. This task is quite relevant since theoretical
assurances allow the usage of our indexes in those applications where a minimum
guarantee is required.
References
[1] Amato, G., Gennaro, C., Savino, P., 2014. Mi-file: using inverted files for
scalable approximate similarity search. Multimedia Tools and Applications
71 (3), 1333–1362.
[2] Amato, G., Savino, P., 2008. Approximate similarity search in metric spaces
using inverted files. In: InfoScale ’08: Proceedings of the 3rd international
conference on Scalable information systems. ICST (Institute for Computer
Sciences, Social-Informatics and Telecommunications Engineering), ICST,
Brussels, Belgium, Belgium, pp. 1–10.
[3] Andoni, A., Indyk, P., January 2008. Near-optimal hashing algorithms for
approximate nearest neighbor in high dimensions. Communications ACM
51, 117–122.
URL http://doi.acm.org/10.1145/1327452.1327494
[4] Andoni, A., Indyk, P., January 2008. Near-optimal hashing algorithms for
approximate nearest neighbor in high dimensions. Communications of the
ACM 51, 117–122.
[5] Baeza-Yates, R. A., Ribeiro-Neto, B. A., 1999. Modern Information Retrieval. ACM Press / Addison-Wesley.
URL citeseer.ist.psu.edu/baeza-yates99modern.html
[6] Bentley, J. L., Saxe, J. B., 1980. Decomposable searching problems i. staticto-dynamic transformation. Journal of Algorithms 1 (4), 301–358.
[7] Beygelzimer, A., Kakade, S., Langford, J., 2006. Cover trees for nearest
neighbor. In: Proceedings of the 23rd international conference on Machine
learning. ACM, pp. 97–104.
[8] Bustos, B., Navarro, G., 2004. Probabilistic proximity search algorithms
based on compact partitions. Journal of Discrete Algorithms 2 (1), 115–
134.
[9] Carlson, J. L., 2013. Redis in Action. Manning Publications Co., Greenwich, CT, USA.

26

[10] Chavez, E., Figueroa, K., Navarro, G., Sep. 2008. Effective proximity retrieval by ordering permutations. IEEE Transactions on Pattern Analysis
and Machine Intelligence 30 (9), 1647–1658.
[11] Chávez, E., Graff, M., Navarro, G., Téllez, E., 2015. Near neighbor
searching with K nearest references . Information Systems 51 (0), 43–61.
URL
http://www.sciencedirect.com/science/article/pii/
S0306437915000241
[12] Chávez, E., Navarro, G., 2003. Probabilistic proximity search: Fighting the
curse of dimensionality in metric spaces. Information Processing Letters 85,
39–46.
[13] Chávez, E., Navarro, G., Baeza-Yates, R., Marroquı́n, J. L., 2001. Searching
in metric spaces. ACM Comput. Surv. 33 (3), 273–321.
[14] Ciaccia, P., Patella, M., Zezula, P., 1997. M-tree: An efficient access
method for similarity search in metric spaces. In: Proceedings of the 23rd
International Conference on Very Large Data Bases. VLDB ’97. Morgan
Kaufmann Publishers Inc., San Francisco, CA, USA, pp. 426–435.
URL http://dl.acm.org/citation.cfm?id=645923.671005
[15] Clarkson, K. L., 1999. Nearest neighbor queries in metric spaces. Discrete
and Computational Geometry 22, 63–93.
[16] Esuli, A., 2009. Pp-index: Using permutation prefixes for efficient and scalable approximate similarity search. In: Proceedings of the 7th Workshop on
Large-Scale Distributed Systems for Information Retrieval (LSDS-IR’09).
Boston, USA, pp. 17–24.
[17] Esuli, A., 2012. Use of permutation prefixes for efficient and scalable approximate similarity search. Information Processing & Management 48 (5),
889–902.
[18] Figueroa, K., Chavez, E., Navarro, G., 2009. Technical description of the
metric spaces library. Technical Report. SISAP Project.
[19] Gao, J., Jagadish, H. V., Lu, W., Ooi, B. C., 2014. DSH: Data sensitive hashing for high-dimensional k-nnsearch. In: Proc. ACM International
Conference on Management of Data (SIGMOD). pp. 1127–1138.
[20] Gionis, A., Indyk, P., Motwani, R., 1999. Similarity search in high dimensions via hashing. In: Proceedings of the 25th International Conference on
Very Large Data Bases. VLDB ’99. Morgan Kaufmann Publishers Inc., San
Francisco, CA, USA, pp. 518–529.
URL http://portal.acm.org/citation.cfm?id=645925.671516
[21] Goyal, N., Lifshits, Y., Schütze, H., 2008. Disorder inequality: a combinatorial approach to nearest neighbor search. In: Proceedings of the 2008
International Conference on Web Search and Data Mining. ACM, pp. 25–
32.
27

[22] Houle, M. E., Nett, M., 2013. Rank cover trees for nearest neighbor search.
In: Similarity Search and Applications. Springer, pp. 16–29.
[23] Indyk, P., 2004. Handbook of Discrete and Computational Geometry., 2nd
Edition. CRC press, chapter 39.
[24] Malkov, Y., Ponomarenko, A., Logvinov, A., Krylov, V., 2012. Scalable
distributed algorithm for approximate nearest neighbor search problem in
high dimensional general metric spaces. In: Proc. 5th International Conference on Similarity Search and Applications (SISAP). pp. 132–147.
[25] Malkov, Y., Ponomarenko, A., Logvinov, A., Krylov, V., 2014. Approximate nearest neighbor algorithm based on navigable small world graphs.
Information Systems 45, 61–68.
[26] Mohamed, H., 2014. Scalable approximate k-nn in multidimensional big
data. phd thesis. Ph.D. thesis, University of Geneva.
[27] Mohamed, H., Marchand-Maillet, S., 2013. Quantized ranking for
permutation-based indexing. In: Brisaboa, N., Pedreira, O., Zezula, P.
(Eds.), Similarity Search and Applications. Vol. 8199 of Lecture Notes in
Computer Science. Springer Berlin Heidelberg, pp. 103–114.
URL http://dx.doi.org/10.1007/978-3-642-41062-8_11
[28] Muja, M., Lowe, D., Nov 2014. Scalable nearest neighbor algorithms for
high dimensional data. Pattern Analysis and Machine Intelligence, IEEE
Transactions on 36 (11), 2227–2240.
[29] Naidan, B., Hetland, M. L., 2012. Static-to-dynamic transformation for
metric indexing structures. In: Navarro, G., Pestov, V. (Eds.), Similarity Search and Applications - 5th International Conference, SISAP 2012,
Toronto, ON, Canada, August 9-10, 2012. Proceedings. Vol. 7404 of Lecture Notes in Computer Science. Springer, pp. 101–115.
URL http://dx.doi.org/10.1007/978-3-642-32153-5_8
[30] Naidan, B., Hetland, M. L., 2013. Static-to-dynamic transformation for
metric indexing structures (extended version). Information Systems.
[31] Navarro, G., Reyes, N., 2008. Dynamic spatial approximation trees. ACM
Journal of Experimental Algorithmics 12, 1.5:1–1.5:68.
URL http://doi.acm.org/10.1145/1227161.1322337
[32] Navarro, G., Reyes, N., 2009. Dynamic spatial approximation trees for
massive data. In: Second International Workshop on Similarity Search and
Applications, SISAP 2009, 29-30 August 2009, Prague, Czech Republic.
pp. 81–88.
URL http://doi.ieeecomputersociety.org/10.1109/SISAP.2009.28
[33] Patella, M., Ciaccia, P., 2009. Approximate similarity search: A multifaceted problem. Journal of Discrete Algorithms 7 (1), 36–48.
28

[34] Samet, H., 2006. Foundations of Multidimensional and Metric Data Structures, 1st Edition. The morgan Kaufman Series in Computer Graphics and
Geometic Modeling. Morgan Kaufmann Publishers, University of Maryland
at College Park.
[35] Silpa-Anan, C., Hartley, R., June 2008. Optimised kd-trees for fast image
descriptor matching. In: Computer Vision and Pattern Recognition, 2008.
CVPR 2008. IEEE Conference on. pp. 1–8.
[36] Skopal, T., Nov. 2007. Unified framework for fast exact and approximate
search in dissimilarity spaces. ACM Transactions on Database Systems
32 (4).
URL http://doi.acm.org/10.1145/1292609.1292619
[37] Skopal, T., 2010. Where are you heading, metric access methods?: a
provocative survey. In: Proceedings of the Third International Conference
on SImilarity Search and APplications. SISAP ’10. ACM, New York, NY,
USA, pp. 13–21.
URL http://doi.acm.org/10.1145/1862344.1862347
[38] Tellez, E. S., Chavez, E., 2010. On locality sensitive hashing in metric
spaces. In: Proceedings of the Third International Conference on SImilarity
Search and APplications. SISAP 2010. ACM, New York, NY, USA, pp. 67–
74.
URL http://doi.acm.org/10.1145/1862344.1862355
[39] Tellez, E. S., Chavez, E., Camarena-Ibarrola, A., November 2009. A brief
index for proximity searching. In: Proceedings of 14th Iberoamerican
Congress on Pattern Recognition CIARP 2009. Lecture Notes in Computer
Science. Springer Verlag, Berlin, Heidelberg, pp. 529–536.
[40] Tellez, E. S., Chavez, E., Graff, M., 2011. Scalable pattern search analysis.
In: 3rd Mexican Congress on Pattern Recognition, MCPR 2011. Springer
Verlag, Lecture Notes in Computer Science.
[41] Tellez, E. S., Chávez, E., Navarro, G., 2011. Succinct nearest neighbor
search. In: Proc. 4th International Conference on Similarity Search and
Applications, SISAP 2011. ACM Press.
[42] Tellez, E. S., Chávez, E., Navarro, G., 2012. Succinct nearest neighbor
search. Information Systems.
[43] Zezula, P., Amato, G., Dohnal, V., Batko, M., 2006. Similarity Search - The
Metric Space Approach, 1st Edition. Vol. 32. Springer, Series: Advances in
Database Systems.

29

