IEEE TRANSACTIONS ON JOURNAL NAME,  MANUSCRIPT ID 1 
 
Efficient and robust approximate nearest 
neighbor search using Hierarchical Navigable 
Small World graphs  
Yu. A. Malkov, D. A. Yashunin 
Abstract — We present a new approach for the approximate K -nearest neighbor search based on navigable small world 
graphs with controllable hierarchy (Hierarchical NSW, HNSW). The proposed solution is fully graph-based, without any need for 
additional search structures, which are typically used at the coarse search stage of the most proximity graph techniques. 
Hierarchical NSW incrementally builds a multi-layer structure consisting from hierarchical set of proximity graphs (layers) for 
nested subsets of the stored elements. The maximum layer in which an element is present is selected randomly with an 
exponentially decaying probability distribution. This allows producing graphs similar to the previously studied Navigable Sma ll 
World (NSW) struc tures while additionally having the links separated by their characteristic distance scales. Starting search 
from the upper layer together with utilizing the scale separation boosts the performance compared to NSW and allows a 
logarithmic complexity scalin g. Additional employment of a heuristic for selecting proximity graph neighbors significantly 
increases performance at high recall and in case of highly clustered data. Performance evaluation has demonstrated that the 
proposed general metric space search i ndex is able to strongly outperform previous  opensource state-of-the-art vector-only 
approaches. Similarity of the algorithm to the skip list structure allows straightforward balanced distributed implementation. 
Index Terms — Graph and tree search strategies, Artificial Intelligence, Information Search and Retrieval, Information Storage 
and Retrieval, Information Technology and Systems, Search process, Graphs and networks, Data Structures, Nearest neighbor 
search, Big data, Approximate search, Similarity search 
——————————      —————————— 
1 INTRODUCTION
onstantly growing amount of the available info r-
mation resources has led to high demand in scalable 
and efficient similarity search data structures. One of the 
generally used approaches for information search is the 
K-Nearest Neighbor Search (K -NNS). The K -NNS a s-
sumes you have a defined distance function between the 
data elements and aims at finding the K elements from 
the dataset which minimize the distance to a given query. 
Such algorithms are used in many applications, such as 
non-parametric machine learn ing algorithms, image fe a-
tures matching in large scale databases [1] and semantic 
document retrieval [2]. A naïve approach to K -NNS is to 
compute the distances between the query and every el e-
ment in the dataset and select the elements with minimal 
distance. Unfortunately, the complexity of the naïve a p-
proach scales linearly with the number of stored elements  
making it infeasible for large-scale datasets. This has led 
to a high interest in development of fast and scalable K -
NNS algorithms.  
Exact solutions for K-NNS [3-5] may offer a substantial 
search speedup only in case of relatively low dimensional 
data du e to “curse of dimensionality”. To overcome this 
problem a concept of Approximate Nearest Neighbors 
Search (K-ANNS) was proposed, which relaxes the cond i-
tion of the exact search by allowing a small number of 
errors. The quality of an inexact search (the re call) is d e-
fined as the ratio between the number of found true nea r-
est neighbors and K. The m ost popular K -ANNS sol u-
tions are based on approximated versions of tree alg o-
rithms [6, 7], locality -sensitive hashing (LSH) [8, 9]  and 
product quantization  (PQ) [10-17]. Proximity graph K -
ANNS algorithms [10, 18-26] have recently gained pop u-
larity offering a better performance on high dimensional 
datasets. However, the power -law scaling of the proxim i-
ty graph routing causes extreme performance degrad a-
tion in case of low dimensional or clustered data.  
In this paper we propose the Hierarchical Navigable 
Small World (Hierarchical NSW, HNSW), a new fully 
graph based incremental K -ANNS structure, which can 
offer a much better logarithmic complexity scaling. The 
main contributions are: explicit selection of the graph’s 
enter-point node, separation of links by  different scales 
and use of an advanced heuristic to select the neighbors. 
Alternatively, Hierarchical NSW algorithm can be seen as 
an extension of the probabilistic skip list structure  [27] 
with proximity graphs instead of the linked lists. Perfo r-
mance e valuation has demonstrated that the proposed 
general metric space method is able to strongly outpe r-
form previous opensource state-of-the-art approaches 
suitable only for vector spaces. 
———————————————— 
 Y. Malkov is with the Federal state budgetary institution of science 
Institute of Applied Physics of the Russian Academy of Sciences, 46 
Ul'yanov Street, 603950 Nizhny Novgorod, Russia. E-mail:  
yurymalkov@mail.ru. 
 D. Yashunin. Addres: 31-33 ul. Krasnozvezdnaya, 603104 Nizhny 
Novgorod, Russia. E-mail: yashuninda@yandex.ru 
C 
This work has been submitted to the IEEE for possible publication. Copyright 
may be transferred without notice, after which this version may no longer be 
accessible.  
2 IEEE TRANSACTIONS ON JOURNAL NAME,  MANUSCRIPT ID 
 
2 RELATED WORKS 
2.1 Proximity graph techniques 
In the vast majority of studied graph algorithms searc h-
ing takes a form of greedy routing in k -Nearest Neighbor 
(k-NN) graphs [10, 18 -26]. For a given proximity graph, 
we start the search at some enter point (it can be random 
or supplied by a separate algorithm) and iteratively t rav-
erse the graph. At each step of the traversal the algorithm 
examines the distances from a query to the neighbors of a 
current base node and then selects as the next base node 
the adjacent node that minimizes the distance, while co n-
stantly keeping track of the best discovered neighbors. 
The search is terminated when some stopping condition 
is met (e.g. the number of distance calculations).  Links to 
the closest neighbors in a k -NN graph serve as a simple 
approximation of the Delaunay graph [25, 26] ( a graph 
which guranties that the result of a basic greedy graph 
traversal is always the nearest neighbor). Unfortunately, 
Delaunay graph cannot be efficiently constructed without 
prior information about the structure of a space [4], but its 
approximation by the nearest neighbors can be done by 
using only distances between the stored elements . It was 
shown that proximity graph approaches with such  
approximation perform competitive to other k-ANNS 
thechniques, such as kd-trees or LSH [18-26]. 
The main drawbacks  of the k-NN graph approaches 
are: 1) the power law scaling of the number of steps with 
the dataset size during the routing process [28, 29] ; 2)  a 
possible loss of global connectivity which leads to poor 
search results  on clusetered data . To overcome these  
problems many hybrid approaches have been proposed 
that use auxiliary algorithms applicable only for vector 
data (such as kd -trees [18, 19] and product 
quantization [10]) to find better candidates for the enter 
nodes by doing a coarse search.  
In [25, 26, 30]  authors proposed a proximity graph 
K-ANNS algorithm called Navigable Small World (NSW, 
also known as Metricized Small World, MSW), which 
utilized navigable graphs, i.e. graphs with logarithmic or 
polylogarithmic scaling of the number of hops during  the 
greedy traversal with the respect of the network 
size [31, 32]. The NSW graph is constructed via consec u-
tive insertion of elements in random order by bidirectio n-
ally connecting them to  the M closest neighbors from the 
previously inserted elements. The M closest neighbors are 
found using the structure’s search procedure (a variant of 
a greedy search from multipl e random enter nodes ). 
Links to the closest neighbors of the elements inserted in 
the beginning of the construction later become bridges  
between the network hubs  that keep the overall graph 
connectivity and allow the logarithmic scaling  of the 
number of hops during greedy routing.  
Construction phase of the NSW  structure can be eff i-
ciently parallelized without global synchronization and 
without mesuarable effect on accuracy [26] , being a good 
choice for distributed search systems. The NSW  approach 
delivered the state -of-the-art performance on some d a-
tasets [33, 34] , h owever, due to  the overall polylogarit h-
mic complexity scaling , the algorithm  was still prone to 
severe performance degradation  on low dimensional d a-
tasets (on which NSW could lose to tree-based algorithms 
by several orders of magnitude [34]). 
2.2 Navigable small world models 
Networks with logarithmic or polylogarithmic scaling of 
the greedy graph r outing are known as the navigable 
small world networks  [31, 32] . Such networks are an i m-
portant topic of complex network theory aiming at u n-
derstanding of underlying mechanisms of real-life ne t-
works formation in order to apply them for appli cations 
of scalable routing [32, 35, 36]  and distributed similarity 
search [25, 26, 30, 37-40].  
The first works to consider spatial models of navigable 
networks were done by J.  Kleinberg [31, 41] as social ne t-
work models for the famous Milgram experiment  [42]. 
Kleinberg studied a variant of random Watts -Strogatz 
networks [43], using a regular lattice graph in d -
dimensional vector space together with augmentation of 
long-range links following a specific long link length di s-
tribution r-. For =d the number of hops to get to the 
target by greedy routing scales polylogarithmically (i n-
stead of a power law for any other value of ). This idea 
has inspired development of many K -NNS and K -ANNS 
algorithms based on the navigation effect [37-40]. But 
even though the  Kleinberg’s navigability criterion in 
principle can be extended for more general spaces, in o r-
der to build such a navigable network one has to know 
the data distribution beforehand . In addition, g reedy 
routing in Kleinberg’s graphs  suffers fr om polylogarith-
mic complexity scalability at best. 
Another well -known class of navigable networks are 
the scale -free models  [32, 35, 36] , which can reproduce 
several features of real -life networks and advertised for 
routing applications  [35]. However, n etworks produced 
by such models have even worse power law complexity 
scaling of the greedy search [44] and, just like the Klei n-
berg’s model, scale-free models require global knowledge 
of the data distribution, making them unusable for search 
applications. 
The above-described NSW algorithm uses a simpler, 
previously unknown model of navigable networks, a l-
lowing decentralized graph construction and suitable for 
data in arbitrary spaces. It was suggested [44] that the 
NSW network formation mechanism may be respon sible 
for navigability of large -scale biological neural networks 
(presence of which is disputable): similar models were 
able to describe growth of small brain networks, while the 
model predicts several high-level features observed in 
large scale neural net works. However, the NSW model 
also suffers from the polylogarithmic search complexity 
of the routing process.  
3 MOTIVATION 
The ways of improving the NSW search complexity can 
be identified through the analysis of the routing process, 
which was studied in detail in [32, 44]. The routing can be 
divided into two phases: “zoom -out” and “zoom-in” [32]. 
The greedy algorithm starts in the “zoom -out” phase 
AUTHOR ET AL.:  TITLE 3 
 
from a low degree node and traverses the graph  simulta-
neously increasing the node’s degree until the characteris-
tic radius of the node links length reaches the scale of the 
distance to the query. Before the latter happens, the ave r-
age degree of a node can stay relatively small, which 
leads to an increased probabilit y of being stuck in a di s-
tant false local minimum.  
One can avoid the described problem in NSW by star t-
ing the search from a node with the maximum degree 
(good candidates are the first nodes inserted in the NSW 
structure [44]), directly going to the “zoom -in” phase of 
the search. Tests show that setting hubs as starting points 
substantially increases probability of successful routing in 
the structure and provides significantly better perfo r-
mance at low dimensional data. However, it still has only 
a polyloga rithmic complexity scalability of a single 
greedy search at best , and performs worse on high d i-
mensional data compared to Hierarchical NSW.  
The reason for the polylogarithmic complexity scaling 
of a single greedy search in NSW is that the overall nu m-
ber of distance computations is roughly proportional to a 
product of the average number of greedy algorithm hops 
by the average degree of the nodes on the greedy path. 
The average number of hops scales logarithmically  
[26, 44], while the average degree of the nodes on the 
greedy path also scales logarithmically due to the facts 
that: 1) the greedy search tends to go through the same 
hubs as the network grows [32, 44]; 2) the  average num-
ber of hub connections grow s logarithmically with an 
increase of the network  size. Thus we get an overall po l-
ylogarithmic dependence of the resulting complexity.  
The idea of Hierarchical NSW algorithm is to separate 
the links according to their length scale  into different lay-
ers and then search in a multilayer graph. In this case we 
can evaluate only a needed  fixed portion of the conne c-
tions for each element independently of the networks 
size, thus allowing a logarithmic scalability. In such struc-
ture the search starts from the upper layer which has only 
the longest links (the “zoom -in” phase) . The algorithm  
greedily traverses through the  elements from the  upper 
layer until a local minimum is reached (see Fig. 1 for illus-
tration). After that, the search switches to the lower layer 
(which has shorter links), restarts from the element which 
was the local minimum in the previous layer  and the pro-
cess repeats . The maximum number of connections per 
element in all layers can be made constant , thus allowing 
a logarithmic complexity scaling of routing in a navigable 
small world network. 
One way to form such a layered structure is to explici t-
ly set links with different length scales by introducing 
layers. For every element we select an integer  level l 
which defines the maximum layer for which the element 
belongs to. For all elements in a layer a proximity graph 
(i.e. graph containing only “short” links that approximate 
Delaunay graph) is built incrementally. If we set an exp o-
nentially decaying probabil ity of l (i.e. following a ge o-
metric distribution) we get a logarithmic scaling of the 
expected number of layers in the structure. The search 
procedure is an iterative greedy search starting from the 
top layer and finishing at the zero layer.  
In case we merge connections from all layers, the struc-
ture becomes similar to the NSW graph (in this case the l 
can be put in correspondence to the node degree in 
NSW). In contrast to NSW, Hierarchical NSW constru c-
tion algorithm does not require the elements to be  shuf-
fled before the insertion - the stochasticity is achieved by 
using level randomization, thus allowing truly increme n-
tal indexing even in case of temporaril y alterating data 
distribution ( though changing the order of the inser tion 
slightly alters the performace due to only partially de-
termenistic construction procedure). 
The Hierarchical NSW idea is also very similar to a 
well-known 1D probabilistic skip list structure  [27] and 
can be described using its terms. The major difference to 
skip list is that we generalize the structure by replacing 
the linked list with proximity graphs. The Hierarchical 
Layer=0
Layer=1
Layer=2
Decreasing characteristic radius
 
Cluster 2
Inserted 
element
e2
Cluster 1  
Fig. 1. Illustration of the Hierarchical NSW idea. The search starts 
from an element from the top layer (shown red). Red arrows show 
direction of the greedy algorithm from the entry point to the query 
(shown green). 
Fig. 2. Illustration of the heuristic used to select the graph neig h-
bors for two isolated clusters. A new ele ment is inserted on the 
boundary of Cluster 1. All of the closest neighbors of the element 
belong to the Cluster 1, thus missing the edges of Delaunay graph 
between the clusters. The heuristic, however, selects element e2 
from Cluster 2, thus, maintaining the global connectivity in case the 
inserted element is the closest to e2 compared to any other element 
from Cluster 1. 
 
4 IEEE TRANSACTIONS ON JOURNAL NAME,  MANUSCRIPT ID 
 
NSW approach thus can utilize the same methods for 
making the distributed approximate search/overlay 
structures [45]. 
For the selection of the proximity graph con nections 
during the element insertion we utilize a heuristic that 
takes into account the distances between the candidate 
elements to create diverse connections (a similar alg o-
rithm was utilized in the spatial approximation tree [4] to 
select the tree children) instead of just selecting the closest 
neighbors. The heuristic examines the candidates starting 
from the nearest (with respect to the inserted element ) 
and creates a connection to a candidate only if it is closer 
to the base  (inserted) element compared to any of the a l-
ready connected candidates (see Section 4 for the details).  
When the number of candidates is large enough the 
heuristic allows getting the exact relative neighborhood 
graph [46] as a subgraph , a minimal subgraph of the D e-
launay graph deducible by using only  the distances b e-
tween the nodes. The relative neighborhood graph allows 
easily keeping the global connected component, even in 
case of highly clustered data (see Fig. 2 for illustration). 
Note that the heuristic creates  extra edges compared to 
the exact relative neighborhood graph s, allowing control-
ling the number of the connections which is important for 
search performance. For the case of 1D data the heuristic 
allows getting the exact Delaunay subgraph (which in 
this case coincides with the relative neighborhood graph) 
by using only information about the distances between 
the elements, thus making a direct transition from Hiera r-
chical NSW to the 1D probabilistic skip list algorithm.  
Base variant of the Hierarchical NSW proximity 
graphs was also used in ref. [18] (called ‘sparse neighbor-
hood graphs’) for proximity graph searching. Similar 
heuristic was also a focus of the FANNG algorithm  [47] 
(published shortly after the first versions of the current 
manuscript were posted online) with a slightly different 
interpretation, based on the sparse neighborhood graph’s 
property of the exact routing [18]. 
4 ALGORITHM DESCRIPTION 
Network construction algorithm (alg. 1)  is organized via 
consecutive insertions of the stored elements into the 
graph structure. For every inserted element an integer  
maximum layer l is randomly selected with an expone n-
tially decaying probability distribution (normalized by 
the mL parameter, see line 4 in alg. 1). 
The first phase of the insertion process starts from the 
top layer by greedily traversing the graph in order to find 
the ef closest neighbors to the inserted element q in the 
layer. After that, the algorithm continues  the search from 
the next layer using the found closest neighbors from the 
previous layer as enter points, and the process repeats.  
Closest neighbors at each layer are found by a variant of 
the greedy search algorithm described in alg. 2, which is 
an updated version of the algorithm from [26]. To obtain 
the approximate ef nearest neighbors in some layer lс, a 
dynamic list W of ef closest found elements (initially 
filled with enter points) is kept during the search. The list 
is updated at each step by evaluating the neighborhood of 
the closest previously non -evaluated element in the list 
until the neighborhood of every element from the l ist is 
evaluated. Compared to limiting the number of distance 
calculations, Hierarchical NSW stop condition has a n ad-
vantage - it allows discarding candidates for evalution 
that are further from the query than the furthest element 
in the list, thus avoiding bloating of search structures. As 
in NSW, the list is emulated via two priority queues for 
better performance. The distinctions from NSW (along 
with some queue optimizations) are: 1) the enter point is a 
fixed parameter; 2)  instead of changing the numb er of 
multi-searches, the quality of the search is controlled by a 
different parameter ef (which was set to K in NSW [26]). 
Algorithm 1 
INSERT(hnsw, q, M, Mmax, efConstruction, mL) 
Input: multilayer graph hnsw, new element q, number of established 
connections M, maximum number of connections for each element 
per layer Mmax, size of the dynamic candidate list efConstruction, nor-
malization factor for level generation mL 
Output: update hnsw inserting element q 
1  W ← ∅    // list for the currently found nearest elements 
2  ep ← get enter point for hnsw 
3  L ← level of ep    // top layer for hnsw 
4  l ← ⌊-ln(unif(0..1))∙mL⌋  // new element’s level 
5  for lc ← L … l+1 
6     W ← SEARCH-LA YER(q, ep, ef=1, lc) 
7     ep ← get the nearest element from W to q 
8  for lc ← min(L, l) … 0 
9      W ← SEARCH-LAYER(q, ep, efConstruction, lc) 
10    neighbors ← SELECT-NEIGHBORS(q, W, M, lc) // alg. 3 or alg. 4 
11    add bidirectionall connectionts from neighbors to q at layer lc 
12    for each e ∈ neighbors   // shrink connections if needed 
13       eConn ← neighbourhood(e) at layer lc    
14       if │eConn│ > Mmax // shrink connections of e 
                                            // if lc = 0 then Mmax = Mmax0 
15          eNewConn ← SELECT-NEIGHBORS(e, eConn,  Mmax, lc) 
                                                                                         // alg. 3 or alg. 4 
16          set neighbourhood(e) at layer lc to eNewConn 
17    ep ← W 
18 if l > L 
19    set enter point for hnsw to q 
Algorithm 2 
SEARCH-LAYER(q, ep, ef, lc) 
Input: query element q, enter points ep, number of nearest to q ele-
ments to return ef, layer number lc 
Output: ef closest neighbors to q 
1  v ← ep      // set of visited elements 
2  C ← ep     // set of candidates  
3  W ← ep    // dynamic list of found nearest neighbors 
4  while │C│ > 0 
5     c ← extract nearest element from C to q 
6     f ← get furthest element from W to q 
7     if distance(c, q) > distance(f, q) 
8        break   // all elements in W are evaluated 
9     for each e ∈ neighbourhood(c) at layer lc   // update C and W 
10        if e ∉ v 
11          v ← v ⋃ e 
12          f ← get furthest element from W to q 
13          if distance(e, q) < distance(f, q) or │W│ < ef 
14             C ← C ⋃ e 
15             W ← W ⋃ e  
16             if │W│ > ef 
17                remove furthest element from W to q  
18 return W 
AUTHOR ET AL.:  TITLE 5 
 
Algorithm 3 
SELECT-NEIGHBORS-SIMPLE(q, C, M) 
Input: base element q, candidate elements C, number of neighbors to 
return M 
Output: M nearest elements to q   
return M nearest elements from C to q 
 
Algorithm 5 
K-NN-SEARCH(hnsw, q, K, ef) 
Input: multilayer graph hnsw, query element q, number of nearest 
neighbors to return K, size of the dynamic candidate list ef 
Output: K nearest elements to q 
1  W ← ∅   // set for the current nearest elements 
2  ep ← get enter point for hnsw 
3  L ← level of ep    // top layer for hnsw 
4  for lc ← L … 1 
5     W ← SEARCH-LA YER(q, ep, ef=1, lc) 
6     ep ← get nearest element from W to q 
7  W ← SEARCH-LA YER(q, ep, ef, lc =0) 
8 return K nearest elements from W to q  
During the first phase of the search  the ef parameter is set 
to 1 (simple greedy search)  to avoid introduction of add i-
tional parameters. 
When the search reaches the layer that is equal or less 
than l, the second phase of the construction algorithm is 
initiated. The second phase differs in two points: 1) the ef 
parameter is increased from 1 to  efConstruction in order 
to control the recall of the greedy search procedure; 2) the 
found closest neighbors on each layer are also used as 
candidates for the connections of the inserted element.  
Two methods for the selection of M neighbors from the 
candidates were considered: simple connection to the 
closest elements (alg. 3) and the heuristic that accounts for 
the distances between the candidate elements to create 
connections in diverse direction s (alg. 4), described in the 
Section 3. The heuristic has two additional parameters: 
extendCandidates (set to false by default) which extends 
the candidate set and useful only for extremely clustered 
data, and keepPrunedConnections which allows getting 
fixed number of connection per element . The maximum 
number of connections that an element can have per layer 
is defined by the parameter Mmax for every layer higher 
than zero (a special parameter Mmax0 is used for the 
ground layer sep arately). If a node is already full at the 
moment of making of a new connection, then its extended 
connection list gets shrunk by the same algorithm that 
used for the neighbors selection (algs. 3 or 4). 
The insertion procedure terminates when the conne c-
tions of the inserted elements are established on the zero 
layer.  
The K -ANNS search algorithm used in Hierarchical 
NSW is presented in alg. 5. It is roughly equivalent to the 
insertion algorithm for an item with  layer l=0. The differ-
ence is that the closest neighbors found at the ground 
layer which are used as candidates for the connections are 
now returned as the search result. The quality of the 
search is controlled by the ef parameter (corresponding to 
efConstruction in the construction algorithm).  
4.1 Influence of the construction parameters 
Algorithm construction parameters mL and Mmax0 are r e-
sponsible for maintaining the small world navigability in 
the constructed graphs. Setting mL to zero (this corr e-
sponds to a single layer in  the graph) and Mmax0 to M 
leads to production of directed k -NN graphs with a po w-
er-law search complexity well studied before [21, 29] (as-
suming using the alg. 3 for neighbor selection). Setting mL 
to zero and Mmax0 to infinity leads to production of NSW 
graphs with polylogarithmi c complexity [25, 26] . Finally, 
setting mL to some non -zero value leads to emergence of 
controllable hierarchy graphs w hich allow  logarithmic 
search complexity by introduction of layers (see the Se c-
tion 3).  
To achieve the optimum performance advantage of the 
controllable hierarchy, the overlap between neighbors on 
different layers (i.e. percent of element neighbors that are 
also belong to other layers) has to be small. In order to 
decrease the overlap we need to decrease the mL. Howev-
er, at the same time, decreasing mL leads to an increase of 
average hop number during a  greedy search on each la y-
er, which negatively affects the performance. This leads to 
existence of the optimal value for the mL parameter.  
A simple choice for the optimal mL is 1/ln(M), this cor-
responds to the skip list parameter p=1/M with an ave r-
age single element overlap between the layers. Simul a-
tions done on an Intel Core  i7 5930K  CPU show that the 
proposed selection of mL is a reasonable choice (see Fig. 3 
for data on 10M random d =4 vectors). In addition, the 
plot demonstrates a massive speedup on low dimensional 
data when increasing the mL from zero  and the effect of 
using the heurist ic for selection of the graph connections . 
It is hard to expect the same behavior for high dime n-
sional data since in this case the k -NN graph already has 
Algorithm 4 
SELECT-NEIGHBORS-HEURISTIC(q, C, M, lc, extendCandidates, keep-
PrunedConnections) 
Input: base element q, candidate elements C, number of neighbors to 
return M, layer number lc, flag indicating whether or not to extend 
candidate list extendCandidates, flag indicating whether or not to add 
discarded elements keepPrunedConnections 
Output: M elements selected by the heuristic 
1  R ← ∅ 
2  W ← C   // working queue for the candidates 
3  if extendCandidates   // extend candidates by their neighbors 
4    for each e ∈ C 
5       for each eadj ∈ neighbourhood(e) at layer lc 
6          if eadj ∉ W 
7             W ← W ⋃ eadj 
8  Wd ← ∅   // queue for the discarded candidates 
9  while │W│ > 0 and │R│< M 
10    e ← extract nearest element from W to q  
11    if e is closer to q compared to any element from R 
12       R ← R ⋃ e 
13    else 
14       Wd ← Wd ⋃ e 
15  if keepPrunedConnections  // add some of the discarded 
                                                   // connections from Wd 
16    while │Wd│> 0 and │R│< M 
17       R ← R ⋃ extract nearest element from Wd to q  
18 return R 
6 IEEE TRANSACTIONS ON JOURNAL NAME,  MANUSCRIPT ID 
 
very short greedy algorithm paths  [28]. Surprisingly, i n-
creasing the mL from zero leads to a measurable increase 
in speed on very high dimensional data (100k dense ra n-
dom d=1024 vectors, see plot in F ig. 4), and does not i n-
troduce any penalty for the Hierarchical NSW approach.  
For real data such as SIFT vectors [1] (which have co m-
plex mixed structure ), the performance improvement by 
increasing the mL is higher, but  less prominent at current 
settings compared to improvement from the heuristic (see 
Fig. 5 for 1 -NN search performance on 5 million 128 -
dimensional SIFT vectors from the learning set of BI G-
ANN [13]).  
Selection of the Mmax0 (the maximum number of co n-
nections that an element can have in the zero layer) also 
has a strong in fluence on the search performance, esp e-
cially in case of high quality (high recall) search. Simul a-
tions show that setting Mmax0 to M (this corresponds to k -
NN graphs on each layer if the neighbors selection heuris-
tic is not used) leads to a very strong performance penalty 
at high recall. Simulations also suggest that 2∙ M is a good 
choice for Mmax0: setting the parameter higher leads to 
performance degradation and excessive memory usage . 
In Fig. 6 there are presented results of search performance 
for the 5M SIFT learn dataset depending on the Mmax0 pa-
rameter (done on an Intel Core i5 2400 CPU) . The su g-
gested value gives performance close to optimal at diffe r-
ent recalls. 
In all of the considered cases, use of the heuristic for 
proximity graph neighbors selection (alg. 4) leads to a 
higher or similar search performance compared to the 
naïve connection to the nearest neighbors (alg. 3). The 
effect is the most prominent for low dimensional data, at 
high recall for mid -dimensional data and for the case of 
highly clust ered data (ideologically discontinuity can be 
regarded as a local low dimensional feature), see the 
comparison in Fig.  7 (Core i5 2400 CPU) . When using the 
closest neighbors as connections for the proximity graph, 
the Hierarchical NSW algorithm fails to ac hieve a high 
recall for clustered data because the search stucks at the 
clusters boundaries. Contrary, when the heuristic is used  
(together with candidates ’ extension, line 3 in Alg. 4) , 
clustering leads to even higher performance. For uniform 
and very high dimensional data there is a little difference 
between the neighbors selecting methods (see Fig. 4), pos-
sibly due to the fact that in this case almost all of the 
nearest neighbors are selected by the heuristic. 
The only meaningful construction parameter left for 
the user is M. A reasonable range of M is from 5 to 48. 
Simulations show that smaller M generally produces be t-
ter results for lower recalls and/or lower dimensional 
data, while b igger M is better for high recall and/or high 
dimensional data (see Fig. 8 for illustration , Core i5 2400 
CPU). The parameter also defines the memory consum p-
tion of the algorithm (which is proportional to M), so it 
should be selected with care. 
Selection o f the efConstruction parameter is straigh t-
forward. As it was suggested in [26] it has to be large 
enough to produce K -ANNS recall close to unity during 
the construction process (0.95 is enough for the most use -
cases). And just like in [26], this parameter can possibly 
0,0 0,5 1,0 1,5 2,0
0,01
0,02
0,03
0,04Query time, ms
mL
10M d=4 random vectors,
M=6, Mmax0=12
Recall 0.9, 1-NN
 Simple neighbors
 Heuristic
Autoselect
 
0,0 0,5 1,0 1,5 2,0
14,1
14,2
14,3
14,4
14,5Query time, ms
mL
100k random vectors, d=1024
M=20, Mmax0=40,
Recall=0.9, 1-NN
 Simple neighbors
 Heuristic
Autoselect  
0,0 0,5 1,0 1,5 2,0
0,09
0,10
0,11
0,12
0,13
0,14
0,15Query time, ms
mL
5M SIFT, d=128,
M=16, Mmax0=32,
Recall=0.9, 1-NN
 Simple Neighbors
 Heuristic
Autoselect  
Fig. 3. Plots for query time vs mL parameter 
for 10M random vectors with d=4. The a u-
toselected value 1/ln( M) for mL is shown by 
an arrow. 
Fig. 4. Plots for query time vs mL parame-
ter for 100k random vectors with d=1024. 
The autoselected value 1/ln( M) for mL is 
shown by an arrow. 
Fig. 5. Plots for query time vs mL parameter 
for 5M SIFT learn dataset. The autoselected 
value 1/ln(M) for mL is shown by an arrow. 
 
0 20 40 60 80 100 120 140 160 180 200
0.1
1
Query time, ms
Mmax0
5M SIFT, d=128, M=20,
mL=0.33, 10-NN
 Recall=0.4
 Recall=0.8
 Recall=0.94
Autoselect
 
0,2
0,4
0,6
0,8
1,0
0,01 0,1
Query time, ms
Recall
10M random vectors, d=10 
M=16, 10-NN
 baseline - no clusters
 heuristic - no clusters
 baseline - 100 clusters
 heuristic - 100 clusters  
10
-3
10
-2
10
-1
10
0
0,01 0,1 1
Query time, ms
Recall error (1-recall)
5M SIFT, d=128,
10-NN
 M=2
 M=3
 M=6
 M=12
 M=20
 M=40  
Fig. 6. Plots for query time vs Mmax0 pa-
rameter for 5M SIFT learn dataset. The 
autoselected value 2∙M for Mmax0 is shown 
by an arrow. 
Fig. 7. Effect of the method of neighbor 
selections (baseline corresponds to alg. 3, 
heuristic to alg. 4) on clustered (100 ra n-
dom isolated clusters) and non-clustered 
d=10 random vector data. 
Fig. 8. Plots for recall error vs query time 
for different parameters of M for Hiera r-
chical NSW on 5M SIFT learn dataset. 
 
AUTHOR ET AL.:  TITLE 7 
 
be auto-configured by using sample data. 
The construction process can be easily and efficiently 
parallelized with only few synchronization points (as 
demonstrated in Fig. 9) and no measurable effect on index 
quality. Construction speed/index quality tradeoff is con-
trolled via  the efConstruction parameter. The tradeoff 
between the search time and the index construction time 
is presented in Fig. 10 for a 10M SIFT dataset and shows 
that a reasonable quality index can be constructed for 
efConstruction=100 on a 4X 2.4 GHz  10-core Xeon E5 -
4650 v2 CPU server in just 3 minutes. Further increase of 
the efConstruction leads to little extra performance but in 
exchange of significantly longer construction time.  
 
4.2 Complexity analysis 
4.2.1 Search complexity 
The complexity scaling of a single search can be strictly 
analyzed under the assumption that we build exact D e-
launay graphs instead of the approximate ones. Suppose 
we have found the closest element on some layer (this is 
guaranteed by having the Delaunay graph) and then d e-
scended to the next layer. One can show that the average 
number of steps before we find the closest element in the 
layer is bounded by a constant.  
Indeed, the layers are not correlated with the  spatial 
positions of the data elements and , thus, when we tra v-
erse the graph there is a fixed probability p=exp(-mL) that 
the next node belongs to the upper layer. However, the 
search on the layer always terminates before it reaches the 
element which belongs to the higher layer (otherwise the 
search on the upper layer would have stopped on a dif-
ferent element), so the probability of not reaching the ta r-
get on  s-th step is bounded by  exp(-s· mL). Thus the ex-
pected number of steps in a layer is bounded by a sum of 
geometric progression S =1/(1-exp(-mL)), which is inde-
pendent of the dataset size. 
If we assume that the average degree of a node in the 
Delaunay graph is capped  by a constant C in the limit of 
the large dataset (this is the case for random Euclid d a-
ta [48], but can be in principle violated in exotic spaces), 
then the overall average number of distance evaluations 
in a layer is bounded by a constant C· S, independently of 
the dataset size.  
And since the expectation of the maximum layer index 
by the construction  scales as O(log(N)), the overall co m-
plexity scaling is O(log(N)), in agreement with  the simu-
lations on low dimensional datasets.  
The inital assumption of having the  exact Delaunay 
graph violates in Hierarchical NSW due to usage of a p-
proximate edge selection heuristic with a fixed number of 
neighbors per element. Thus, to avoid stucking into a l o-
cal minimum the greedy search algo rithm em ploys a 
backtracking procedure on the zero layer.  Simulations 
show that at least for  low dimensional data  (Fig. 11, d=4) 
the dependence of the  required ef parameter (which de-
termines the complexity via the minimal number of hops 
during the backtracking) to get a fixed recall  saturates 
with the rise of the dataset size . The b acktracking co m-
plexity is an additive term in respect to the final comple x-
ity, thus, as follows from the empirical data,  inaccuracies 
of the Delauna y graph approximation do not alter the 
scaling.  
Such empirical investigation  of the  Delaunay graph 
approximation resilience requires having the average 
number of Delaunay graph edges independent of the d a-
taset to evidence how well the edges  are approximated 
with a constant number of connections in Hierarchical 
NSW. However, the average degree of Delaunay graph 
scales exponentially with the dimensionality  [39]), thus 
for high dimensional data (e.g. d=128) the aforeme n-
tioned condition requires having extremely larg e d a-
tasets, making such empricial investigation unfeasible. 
Further analitical evidence is required to confirm whether 
the resilience of Delaunay graph aproximations genera l-
izes to higher dimensional spaces.  
4.2.2 Construction complexity 
The construction is done by iterative insertions of all el e-
ments, while the  insertion of an element is  merely a s e-
quence of K-ANN-searches at different layers  with a sub-
sequent use of heuristic (which has fixed complexity at 
fixed efConstruction). The average number of lay ers for 
an element to be added  in is a constant that depends on 
mL: 
 
   1 ln( (0,1)) 1 1 LLE l E unif m m     
   (1) 
Thus, the insertion complexiy scaling is the same as the 
one for the search, meaning that at least  for relatively low 
dimensional datasets the construction time scales as 
O(N∙log(N)). 
0 5 10 15 20 25 30 35 40
0
10
20
30
40
50
60
70
80
90Build time, minutes
Thread count
10M SIFT, d=128, M=16, 
efConstruction=100
 4X Xeon E5-4650 v2 (4x10 cores)
 Core i7-6850K (6 cores+HT)
 
0 2 4 6 8 10 12
0.2
0.3
0.4
0.5
0.6
0.7Query time, ms
Build time, minutes
10M SIFT, d=128, M=16,
Recall 0.9, 1-NN  
10
0
10
1
10
2
10
3
10
4
10
5
10
6
10
7
10
8
10
9
0
2
4
6
8
10
12Average ef to reach target recall
Dataset size
Random vectors, d=4
M=6, Mmax0=12, 1-NN
 Recall=0.9
 Recall=0.99
 Recall=0.999  
Fig. 9. Construction time for Hierarchical 
NSW on 10M SIFT dataset for different 
numbers of threads on two CPUs. 
Fig. 10. Plots of t he query time vs constru c-
tion time tradeoff for Hierarchical NSW on 
10M SIFT dataset. 
Fig. 11. Plots of the ef parameter required 
to get fixed accuracies vs the dataset size 
for d=4 random vector data. 
 
8 IEEE TRANSACTIONS ON JOURNAL NAME,  MANUSCRIPT ID 
 
10
-5
10
-4
10
-3
10
-2
10
-1
10
0
100
1000
10000
100000Distance computations
Recall error (1-recall)
(a) 10M random vectors, d=4
M=6, 10-NN
 NSW
 Hierarchical NSW
10
2
10
3
10
4
10
5
10
6
10
7
10
8
0
200
400
600
800
1000
1200Distance computations
Dataset size
(b) Random vectors, d=8
M=6, Recall=0.95, 10-NN
 NSW
 Hierarchical NSW
10
2
10
3
10
4
10
5
10
6
10
7
10
8
0,0
0,5
1,0
1,5
2,0
2,5Query time, ms
Dataset size
(c) Random vectors, d=8
M=6, Recall=0.95, 10-NN
 NSW
 Hierarchical NSW
 
Fig. 12. Comparison between NSW and Hierarchical NSW: (a) distance calculation number vs accuracy tradeoff for a 10 million 4 -
dimensional random vectors dataset; (b-c) performance scaling in terms of number of distance calculations (b) and raw query(c) time on 
a 8-dimensional random vectors dataset. 
 
4.2.3 Memory cost 
The memory consumption of the Hierarchical NSW is 
mostly defined by the  storage of  graph connections. The 
number of connection s per element is Mmax0 for the zero  
layer and Mmax for all other l ayers. Thus, the average 
memory consumption per element is  
(Mmax0+mL ∙Mmax)∙bytes_per_link. If we limit the maximum 
total number of elements by approximately four billions, 
we can use four -byte unsigned integers to store the co n-
nections. Tests suggest that typical close to optimal M 
values usually lie in a range between 6 and 48. This 
means that the typical memory requirements for the i n-
dex (excluding the size of the data) are about 60-450 bytes 
per object, which is in a good agreement with the simul a-
tions.  
5 PERFORMANCE EVALUATION 
The Hierarchical NSW algorithm was implemented in 
C++ on top of the Non Metric Space Library 
(nmslib) [49]1, which already had a functional NSW i m-
plementation (under name “sw -graph”). D ue to several 
limitations posed by the library, to achieve a better pe r-
formance, the Hierarchical NSW implementation uses 
custom distance functions together with C -style memory 
management, which avoids unnecessary implicit address-
ing and allows efficient hardware and software prefetc h-
ing during the graph traversal. 
Comparing the performance of K -ANNS algorithms is 
a nontrivial task since the state -of-the-art is constantly 
changing as new algorithms and implementations are 
emerging. In this work we concentr ated on comparison 
with the best algorithms in Euclid spaces that have open 
source implementations. An implementation of the Hie r-
archical NSW algorithm presented in this paper is also 
distributed as a part of the open source nmslib library 1 
together with a n external C++ memory-efficient header-
only version with support for incremental index construc-
tion2.  
The comparison section consists of four parts: compar-
ison to the baseline NSW (5.1), comparison to the state -of-
the-art algorithms in Euclid spaces (5.2) , rerun of the su b-
 
1 https://github.com/searchivarius/nmslib 
2 https://github.com/nmslib/hnsw 
set of tests [34] in general metric spaces in which NSW 
failed (5.3) and comparison to state -of-the-art PQ -
algorithms on a large 200M SIFT dataset (5.4). 
5.1 Comparison with baseline NSW 
For the baseline NSW algorithm implementation, we used 
the “sw-graph” from nmslib  1.1 (which is slightly updat-
ed compared to the implementation tested in [33, 34] ) to 
demonstrate the improvements in speed and algorithmic 
complexity (measured by the number of distance comp u-
tations). 
Fig. 12(a) presents a co mparison of Hierarchical NSW 
to the basic NSW algorithm for d=4 random hypercube 
data made on a Core i5 2400 CPU (10-NN search). Hierar-
chical NSW uses much less distance computations during 
a search on the dataset, especially at high recalls. 
The scalings of the algorithms on a d=8 random hype r-
cube dataset for a 10 -NN search with a fixed recall of 0.95 
are presented in Fig. 12 (b). It clearly demostrates that Hi-
erarchical NSW has a complexity scaling for this setting 
not worse than logarithmic and outperforms NSW at any 
dataset size. The performance advantage in absolute time 
(Fig. 12(c)) is even higher due to improved algorithm i m-
plementaion. 
5.2 Comparison in Euclid spaces 
The main part of the comparison was carried out on ve c-
tor datasets with use  of the popular K -ANNS benchmark 
ann-benchmark3 as a testbed. The testing system utilize s 
python bindings of the algorithms  – it consequentially 
runs the K -ANN search for one thousand queries (ra n-
domly extracted from the initial dataset) with preset alg o-
rithm parameters producing an output containing recall 
and average time of a single search. The considered alg o-
rithms are:  
1. Baseline NSW algorithm from nmslib 1.1 
(“sw-graph”). 
2. FLANN 1.8.4 [6]. A popular library 4 containing 
several algorithms, built-in in OpenCV5. We used 
the available auto-tuning procedure with several 
reruns to infer the best parameters. 
3. Annoy6, 02.02.2016 build. A popular algorithm 
 
3 https://github.com/erikbern/ann-benchmarks 
4 https://github.com/mariusmuja/flann 
5 https://github.com/opencv/opencv 
6 https://github.com/spotify/annoy 
AUTHOR ET AL.:  TITLE 9 
 
based on random projection tree forest.  
4. VP-tree. A general metric space algorithm with 
metric pru ning [50] implemented as a part of 
nmslib 1.1. 
5. FALCONN7, version 1.2. A new efficient LSH a l-
gorithm for cosine similarity data [51]. 
The comparison was done on a 4X Xeon E5 -4650 v2 
Debian OS system with 128 Gb of RAM. For every alg o-
rithm we carefully chose  the best results at every recall 
range to evaluate the best possible performance (with 
initial values from the testbed defaults). All tests were 
done in a single thread regime. Hierarchical NSW was 
compiled using the GCC 5.3 with -Ofast optimization 
flag. 
The parameters and description of the used datasets 
are outlined in Table 1. For all of the datasets except 
GloVe we used the L 2 distance. For GloVe we used the 
cosine similarity which is equivalent to L 2 after vector 
normalization. The brute -force (BF) t ime is measured by 
the nmslib library.  
Results for the vector data are presented in Fig. 13. For 
SIFT, GloVE, DEEP and CoPhIR datasets Hierarchical 
NSW clearly outperforms the rivals by a large margin. 
For low dimensional data (d=4) Hiera rchical NSW is 
 
7 https://github.com/FALCONN-LIB/FALCONN 
slightly faster at high recall compared to the Annoy while 
strongly outperforms the other algorithms.  
5.3 Comparison in general spaces 
A recent comparison of algorithms [34] in general spaces 
(i.e. non-symmetric or with violation of triangle inequal i-
ty) showed that the baseline NSW algorithm has severe 
problems on low dimensional datasets. To test the pe r-
formance of the Hierarchical NSW algorithm we have 
repeated a subset of tests from [34] on which NSW pe r-
formed poorly or suboptimal. For that purpose we used a 
built-in nmslib testing system which had scripts to run 
tests from [34]. The evaluated algorithms included the 
VP-tree, permutation techniques (NAPP and bruteforce 
filtering) [49, 55 -57], the basic NSW algorithm and 
NNDescent-produced proximity graphs  [29] (both in pair 
with the NSW graph search algorithm). As in the original 
tests, for every dataset the test includes the results of e i-
ther NSW or NNDescent, depending on which structure 
performed better. No custom distance functions or special 
TABLE 1 
 Parameters of the used datasets on vector spaces 
benchmark. 
Dataset Description Size d BF time Space 
SIFT Image feature vectors [13] 1M 128 94 ms L2 
GloVe Word embeddings trained on tweets [52] 1.2M 100 95 ms cosine 
CoPhIR MPEG-7 features extracted from the images [53] 2M 272 370 ms L2 
Random vectors Random vectors in hypercube 30M 4 590 ms L2 
DEEP One million subset of the billion deep image 
features dataset [14] 
1M 96 60 ms L2 
MNIST Handwritten digit images [54] 60k 784 22 ms L2 
 
 
TABLE 2. 
 Used datasets for repetition of the Non-Metric data tests 
subset. 
Dataset Description Size d BF time Distance 
Wiki-sparse TF-IDF (term frequency –inverse document 
frequency) vectors (created via GENSIM [58]) 
4M 105 5.9 s Sparse cosine  
 
Wiki-8 Topic histograms created from sparse TF -IDF 
vectors of the wiki -sparse dataset (created via 
GENSIM [58]) 
2M 8 - Jensen–
Shannon (JS) 
divergence 
Wiki-128 Topic histograms created from sparse TF -IDF 
vectors of the wiki -sparse dataset (created via 
GENSIM [58]) 
2M 128 1.17 s Jensen–
Shannon (JS) 
divergence 
ImageNet Signatures extracted from LSVRC -2014 with 
SQFD (signature quadratic form) distance [59] 
1M 272 18.3 s SQFD 
DNA DNA (deoxyribonucleic acid) dataset sampled 
from the Human Genome 5 [34]. 
1M - 2.4 s Levenshtein 
 
 
0,0
0,2
0,4
0,6
0,8
1,0
0,01 0,1 1 10 100
Query time, ms
Recall
1M SIFT, 10-NN
BF: 94 ms
 Hierarchical NSW
 NSW
 Annoy
 VP-tree
 FLANN
0,0
0,2
0,4
0,6
0,8
1,0
0,01 0,1 1 10 100
Query time, ms
Recall
1.2M GloVe, 10-NN
BF: 95 ms
 Hierarchical NSW
 NSW
 Annoy
 VP-tree
 FLANN
 FALCONN
0,01 0,1 1 10 100
0,0
0,2
0,4
0,6
0,8
1,0Recall
Query time, ms
2M CoPhIR, 10-NN
BF: 370 ms
 Hierarchical NSW
 NSW
 Annoy
 VP-tree
 FLANN
0,01 0,1 1
0,2
0,4
0,6
0,8
1,0Recall
Query time, ms
30M random d=4,10-NN
BF: 590 ms
 Hierachical NSW
 NSW
 Annoy
 VP-tree
 FLANN
0,1 1 10
0,0
0,2
0,4
0,6
0,8
1,0Recall
Query time, ms
60k MNIST, 10-NN
BF: 22 ms
 Hierarchical NSW
 NSW
 Annoy
 VP-tree
 FLANN
0,01 0,1 1 10 100
0,0
0,2
0,4
0,6
0,8
1,0Recall
Query time, ms
1M Deep,10-NN
BF: 60 ms
 Hierarchical NSW
 NSW
 Annoy
 VP-tree
 FLANN
 
Fig. 13. Results of the comparison of Hierarchical NSW with open source implementations of K-ANNS algorithms on five datasets for 10-
NN searches. The time of a brute-force search is denoted as the BF. 
 
10 IEEE TRANSACTIONS ON JOURNAL NAME,  MANUSCRIPT ID 
 
memory management were used in this case for Hiera r-
chical NSW leading to some performance loss.  
The datasets are summarized in Table 2. Further d e-
tails of the datasets , spaces and a lgorithm parameter s e-
lection can be found in the original work  [34]. The brute -
force (BF) time is measured by the nmslib library. 
The results are presented in Fig. 14. Hierarchical NSW 
significantly improves the performance of NSW and is a 
leader for any of the tested datasets. The strongest e n-
hancement over NSW, almost by 3 orders of magnitude is 
observed for the dataset with the lo west dimensionality, 
the wiki-8 with JS -divergence. This is an important result 
that demonstrates the robustness of Hierarchical NSW , as 
for the original NSW this dataset was a stumbling block. 
Note that for the wiki -8 to nullify the effect of impleme n-
tation results are presented for the distance computations 
number instead of the CPU time. 
5.4 Comparison with product quantization based 
algorithms. 
Product quantization K -ANNS algorithms [10-17] are 
considered as the state -of-the-art on billion scale dataset s 
since they can efficiently compress stored data , allowing 
modest RAM usage while achieving millisecond search 
times on modern CPUs.  
To compare the performance of Hierarchical NSW 
against PQ algorithms we used the facebook Faiss library 8 
as the baseline (a new library with state -of-the-art PQ 
algorithms [12, 15]  implementations, released after  the 
current manuscript was submitted) compiled with the 
OpenBLAS backend. The tests where done for a 200M 
subset of 1B  SIFT dataset  [13] on a 4X Xeon E5 -4650 v2 
server with 128Gb of RAM. The ann-benchmark testbed 
was not feasible for these experiments because of its reli-
ance on 32 -bit floating point format (requiring more than 
100 Gb just to store the data). To get the results for Faiss 
PQ algorithms we have utilized  built-in scripts with the 
parameters from Faiss wiki 9. For the Hierarchical NSW 
algorithm we used a special build outside of the nmslib 
with a small memory footprint, simple non -vectorized 
 
8 https://github.com/facebookresearch/faiss 2017 May build. From 2018 
Faiss library has its own implementation of Hierarchical NSW. 
9 https://github.com/facebookresearch/faiss/wiki/Indexing-1G-vectors 
TABLE 3. 
 Parameters for comparison between Hierarchical NSW 
and Faiss on a 200M subset of 1B SIFT dataset. 
Algorithm Build time Peak memory (runtime) Parameters 
Hierarchical NSW 5.6 hours 64 Gb M=16, efConstruction=500 (1) 
Hierarchical NSW 42 minutes 64 Gb M=16, efConstruction=40 (2) 
Faiss 12 hours 30 Gb OPQ64, IMI2x14, PQ64 (1) 
Faiss 11 hours 23.5 Gb OPQ32, IMI2x14, PQ32 (2) 
 
 
0,0
0,2
0,4
0,6
0,8
1,0
0,1 1 10 100
Query time, ms
Recall
200M SIFT, 1-NN
 Faiss 1
 Faiss 2
 HNSW 1
 HNSW 2
1 10 100
0,1
0,2
0,3
0,4
0,5
0,6Query time, ms
Dataset size, millions
 
Fig. 15 Results of comparison with Faiss library on the 200M  
SIFT dataset from [13]. The inset shows the scaling of the query 
time vs the dataset size for Hierarchical NSW. 
 
0,2
0,4
0,6
0,8
1,0
0,1 1 10 100
Query time, ms
Recall
2M Wiki-128 (JS-div), 10-NN
BF: 1.17 s
 Hierarchical NSW
 NSW
 VP-tree
 NAPP
0,80
0,85
0,90
0,95
1,00
10 100 1000
Query time, ms
Recall
1M ImageNet, 10-NN
BF: 18.3 s
 Hierarchical NSW
 NSW
 VP-tree
 NAPP
 bruto-force filtering
0,6
0,7
0,8
0,9
1,0
100 1000 10000 100000
Number of distance computations
Recall
2M Wiki-8 (JS-div), 10-NN
 Hierarchical NSW
 Old kNN(NNDescent)
 VP-tree
 NAPP
0,6
0,8
1,0
10 100 1000
Query time, ms
Recall
1M DNA (edit dist), 10-NN
BF: 2.4 s 
 Hierarchical NSW
 Old kNN (NNDescent)
 VP-tree
 NAPP
 brute-force filt. bin.
0,0
0,2
0,4
0,6
0,8
1,0
1 10 100 1000
Query time, ms
Recall
4M Wiki-sparse, 10-NN
BF: 5.9 s
 Hierarchical NSW
 NSW
 NAPP
 
Fig. 14. Results of the comparison of Hierarchical NSW with general space K -ANNS algorithms from the Non Metric Space Library on 
five datasets for 10-NN searches. The time of a brute-force search is denoted as the BF. 
 
AUTHOR ET AL.:  TITLE 11 
 
integer distance functions and support for incremental 
index construction10.  
The results are presented in Fig. 15 with summariz a-
tion of the parameters in Table 3. The peak memory co n-
sumption was measured by using linux “time –v” tool in 
separate test runs after index construction for both of the 
algorithms. Even though Hierarchical NSW requires si g-
nificantly more RAM, it can achieve much higher accur a-
cy, while offering a massive advance in search speed and 
much faster index construction.  
The inset in Fig. 15 presents the scaling of the query 
time vs the dataset size for Hierarchical NSW.  Note that 
the scaling deviates from the pure logarithm , possibly 
due to relatively high dimensionality of the dataset. 
6 DISCUSSION 
By using structure decomposition of navigable small 
world graphs together with the smart neighbor selectio n 
heuristic the proposed Hierarchical NSW approach ove r-
comes several important problems of the basic NSW 
structure advancing the state-of–the-art in K-ANN search. 
Hierarchical NSW offers an excellent performance and is 
a clear leader on a large variety of the datasets, surpa s-
sing the  opensource rivals by a large margin in case of 
high dimensional data. Even for the datasets where the 
previous algorithm (NSW) has lost by orders of magn i-
tude, Hierarchical NSW was able to come first. Hiera r-
chical NSW supports continuous incremental indexing 
and can also be used as an efficient method for getting 
approximations of the k-NN and relative neighborhood 
graphs, which are byproducts of the index construction. 
Robustness of the approach is a strong feature which 
makes it very attractive for practical applications. The 
algorithm is applicable in generalized metric spaces pe r-
forming the best on any of the datasets tested in this p a-
per, and thus eliminating the need for complicated sele c-
tion of the best algorithm for a specific problem. We stress 
the importance of the algorithm’s robustness since the 
data may have a complex structure with different effe c-
tive dimensionality across the scales. For instance, a d a-
taset can consist of points lying on a curve that randomly 
fills a high dimensional cube, thus being high dimensio n-
al at  large scale and low dimensional at small scale . In 
order to perform efficient search in such datasets an a p-
proximate nearest neighbor algorithm has to work well 
for both cases of high and low dimensionality. 
There are several ways to further increase the efficie n-
cy and applicability of the Hierarchical NSW approach. 
There is still one meaningful parameter left which stron g-
ly affects the construction of the index – the number of 
added connections per la yer M. Potentially, this parame-
ter can be inferred directly by using different 
heuristics [4]. It would also be interesting to compare H i-
erarchical NSW on the full 1B SIFT and 1B DEEP 
datasets [10-14] and add support for element updates and 
removal. 
One of the apparent shortcomings of the proposed a p-
 
10 https://github.com/nmslib/hnsw 
proach compared to the basic NSW is the loss of the po s-
sibility of distributed search. The search in the Hiera r-
chical NSW structure always starts from the top layer, 
thus the structure cannot be made distributed by using 
the same techniques as described in  [26] due to cognes-
tion of  the higher layer elements. Simple workarounds 
can be used to distribute the structure, such as partitio n-
ing the data across cluster nodes studied in  [6], however 
in this case, the total parallel throughput of the system 
does not scale well with the number of computer nodes. 
Still, there are other possible known ways to make this 
particular structure distributed. Hierarchical NSW is ide-
ologically very  similar to the well -known one -
dimensional exact search probabilistic skip list st ructure, 
and thus can use the same techniques to make the stru c-
ture distributed [45]. Potentially this can lead to even be t-
ter distributed performance compared to the base NSW 
due to logarithmic scalability and ideally uniform load on 
the nodes.  
7 ACKNOWLEDGEMENTS 
We thank Leonid Boytsov for many helpful discussions, 
assistance with Non-Metric Space Library integration and 
comments on the manuscript. We thank Seth Hoffert  and 
Azat Davletshin  for the suggestions on the manuscript 
and the algorithm and fellows who contributed to the 
algorithm on the github repository. We also thank Valery 
Kalyagin for support of this work. 
The reported study was funded by RFBR, according to 
the research project No. 16-31-60104 mol_а _dk. 
8 REFERENCES 
[1] D. G. Lowe, "Distinctive image features from scale -invariant 
keypoints," International journal of computer vision, vol. 60, no. 2, 
pp. 91-110, 2004. 
[2] S. Deerwester, S. T. Dumais, T. K. Landauer, G. W. Furnas, and 
R. A . Harshman, "Indexing by Latent Semantic Analysis," J. 
Amer. Soc. Inform. Sci., vol. 41, pp. 391-407, 1990. 
[3] P. N. Yianilos, "Data structures and algorithms for nearest 
neighbor search in general metric spaces," in SODA, 1993, vol. 
93, no. 194, pp. 311-321. 
[4] G. Navarro, "Searching in metric spaces by spatial approxim a-
tion," The VLDB Journal, vol. 11, no. 1, pp. 28-46, 2002. 
[5] E. S. Tellez, G. Ruiz, and E. Chavez, "Singleton indexes for 
nearest neighbor search," Information Systems, 2016. 
[6] M. Muja and D. G. Lowe, "Scalable nearest neighbor algorithms 
for high dimensional data," Pattern Analysis and Machine Intell i-
gence, IEEE Transactions on, vol. 36, no. 11, pp. 2227-2240, 2014. 
[7] M. E. Houle and M. Nett, "Rank-based similarity search: Reduc-
ing the dimensional depen dence," Pattern Analysis and Machine 
Intelligence, IEEE Transactions on, vol. 37, no. 1, pp. 136 -150, 
2015. 
[8] A. Andoni, P. Indyk, T. Laarhoven, I. Razenshteyn, and L. 
Schmidt, "Practical and optimal LSH for angular distance," in 
Advances in Neural Informati on Processing Systems , 2015, pp. 
1225-1233. 
[9] P. Indyk and R. Motwani, "Approximate nearest neighbors: 
towards removing the curse of dimensionality," in Proceedings of 
the thirtieth annual ACM symposium on Theory of computing , 1998, 
pp. 604-613: ACM. 
12 IEEE TRANSACTIONS ON JOURNAL NAME,  MANUSCRIPT ID 
 
[10] J. Wang , J. Wang, G. Zeng, R. Gan, S. Li, and B. Guo, "Fast 
neighborhood graph search using cartesian concatenation," in 
Multimedia Data Mining and Analytics : Springer, 2015, pp. 397 -
417. 
[11] M. Norouzi, A. Punjani, and D. J. Fleet, "Fast exact search in 
hamming space with multi-index hashing," Pattern Analysis and 
Machine Intelligence, IEEE Transactions on, vol. 36, no. 6, pp. 
1107-1119, 2014. 
[12] A. Babenko and V. Lempitsky, "The inverted multi -index," in 
Computer Vision and Pattern Recognition (CVPR), 2012 IEEE Con-
ference on, 2012, pp. 3069-3076: IEEE. 
[13] H. Jegou, M. Douze, and C. Schmid, "Product quantization for 
nearest neighbor search," Pattern Analysis and Machine Intell i-
gence, IEEE Transactions on, vol. 33, no. 1, pp. 117-128, 2011. 
[14] A. Babenko and V. Lempitsky , "Efficient indexing of billion -
scale datasets of deep descriptors," in Proceedings of the IEEE 
Conference on Computer Vision and Pattern Recognition , 2016, pp. 
2055-2063. 
[15] M. Douze, H. Jégou, and F. Perronnin, "Polysemous codes," in 
European Conference on  Computer Vision , 2016, pp. 785 -801: 
Springer. 
[16] Y. Kalantidis and Y. Avrithis, "Locally optimized product quan-
tization for approximate nearest neighbor search," in Proceed-
ings of the IEEE Conference on Computer Vision and Pattern Reco g-
nition, 2014, pp. 2321-2328. 
[17] P. Wieschollek, O. Wang, A. Sorkine -Hornung, and H. Lensch, 
"Efficient large -scale approximate nearest neighbor search on 
the gpu," in Proceedings of the IEEE Conference on Computer Vision 
and Pattern Recognition, 2016, pp. 2027-2035. 
[18] S. Arya and D.  M. Mount, "Approximate Nearest Neighbor 
Queries in Fixed Dimensions," in SODA, 1993, vol. 93, pp. 271 -
280. 
[19] J. Wang and S. Li, "Query -driven iterated neighborhood graph 
search for large scale indexing," in Proceedings of the 20th ACM 
international conference on Multimedia, 2012, pp. 179-188: ACM. 
[20] Z. Jiang, L. Xie, X. Deng, W. Xu, and J. Wang, "Fast Nearest 
Neighbor Search in the Hamming Space," in MultiMedia Model-
ing, 2016, pp. 325-336: Springer. 
[21] E. Chávez and E. S. Tellez, "Navigating k -nearest neighbor 
graphs to solve nearest neighbor searches," in Advances in Pa t-
tern Recognition: Springer, 2010, pp. 270-280. 
[22] K. Aoyama, K. Saito, H. Sawada, and N. Ueda, "Fast approx i-
mate similarity search based on degree -reduced neighborhood 
graphs," in Proceedings of the 17th ACM SIGKDD international 
conference on Knowledge discovery and data mining, 2011, pp. 1055-
1063: ACM. 
[23] G. Ruiz, E. Chávez, M. Graff, and E. S. Téllez, "Finding Near 
Neighbors Through Local Search," in Similarity Search and A p-
plications: Springer, 2015, pp. 103-109. 
[24] R. Paredes, "Graphs for metric space searching," PhD thesis, 
University of Chile, Chile, 2008. Dept. of Computer Science 
Tech Report TR/DCC -2008-10. Available at 
http://www.dcc.uchile.cl/~raparede/publ/08PhDthesis.pdf, 
2008. 
[25] Y. Malkov, A. Ponomarenko, A. Logvinov, and V. Krylov, 
"Scalable distributed algorithm for approximate nearest neig h-
bor search problem in high dimensional general metric spaces," 
in Similarity Search and Applications: Springer Berlin Heidelberg, 
2012, pp. 132-147. 
[26] Y. Malkov, A. Ponomarenko, A. Logvinov, and V. Krylov, "A p-
proximate nearest neighbor algorithm based on navigable small 
world graphs," Information Systems, vol. 45, pp. 61-68, 2014. 
[27] W. Pugh, "Skip lists: a probabilistic alternative to balanced 
trees," Communications of the ACM, vol. 33, no. 6, pp. 668 -676, 
1990. 
[28] C. C. Cartozo and P. De Los Rios, "Extended navigability of 
small world networks: exact results and new insights," Physical 
review letters, vol. 102, no. 23, p. 238703, 2009. 
[29] W. Dong, C. Mo ses, and K. Li, "Efficient k -nearest neighbor 
graph construction for generic similarity measures," in Proceed-
ings of the 20th international conference on World wide web , 2011, 
pp. 577-586: ACM. 
[30] A. Ponomarenko, Y. Malkov, A. Logvinov, and V. Krylov, "A p-
proximate Nearest Neighbor Search Small World Approach," in 
International Conference on Information and Communication Tec h-
nologies & Applications, Orlando, Florida, USA, 2011. 
[31] J. M. Kleinberg, "Navigation in a small world," Nature, vol. 406, 
no. 6798, pp. 845-845, 2000. 
[32] M. Boguna, D. Krioukov, and K. C. Claffy, "Navigability of 
complex networks," Nature Physics, vol. 5, no. 1, pp. 74-80, 2009. 
[33] A. Ponomarenko, N. Avrelin, B. Naidan, and L. Boytsov, 
"Comparative Analysis of Data Structures for Approximate 
Nearest Neighbor Search," In Proceedings of The Third Intern a-
tional Conference on Data Analytics, 2014. 
[34] B. Naidan, L. Boytsov, and E. Nyberg, "Permutation search 
methods are efficient, yet faster search is possible," VLDB Pro-
cedings, vol. 8, no. 12, pp. 1618-1629, 2015. 
[35] D. Krioukov, F. Papadopoulos, M. Kitsak, A. Vahdat, and M. 
Boguná, "Hyperbolic geometry of complex networks," Physical 
Review E, vol. 82, no. 3, p. 036106, 2010. 
[36] A. Gulyás, J. J. Bíró, A. K őrösi, G. Rétvári, and D. Krioukov, 
"Navigable networks as Nash equilibria of navigation games," 
Nature Communications, vol. 6, p. 7651, 2015. 
[37] Y. Lifshits and S. Zhang, "Combinatorial algorithms for nearest 
neighbors, near-duplicates and small-world design," in Proceed-
ings of the Twentieth Annual ACM -SIAM Symposium on Discrete 
Algorithms, 2009, pp. 318 -326: Society for Industrial and A p-
plied Mathematics. 
[38] A. Karbasi, S. Ioannidis, and L. Massoulie, "From Small -World 
Networks to Comparison -Based Search," Information Theory, 
IEEE Transactions on, vol. 61, no. 6, pp. 3056-3074, 2015. 
[39] O. Beaumont, A. -M. Kermarrec, and É. Rivière, "Peer to peer 
multidimensional overlays: Approximating complex stru c-
tures," in Principles of Distributed Systems : Springer, 2007, pp. 
315-328. 
[40] O. Beaumont, A. -M. Kermarrec, L. Marchal, and É. Rivière, 
"VoroNet: A scalable object network based on Voronoi tessell a-
tions," in Parallel and Distributed Processing Symposium, 2007. 
IPDPS 2007. IEEE International, 2007, pp. 1-10: IEEE. 
[41] J. Kleinberg, "The small -world phenomenon: An algorithmic 
perspective," in Proceedings of the thirty-second annual ACM sym-
posium on Theory of computing, 2000, pp. 163-170: ACM. 
[42] J. Travers and S. Milgram, "An experimental study of the small 
world problem," Sociometry, pp. 425-443, 1969. 
[43] D. J. Watts and S. H. Strogatz , "Collective dynamics of ‘small -
world’networks," Nature, vol. 393, no. 6684, pp. 440-442, 1998. 
[44] Y. A. Malkov and A. Ponomarenko, "Growing homophilic ne t-
works are natural navigable small worlds," PloS one, p. 
e0158162, 2016. 
[45] M. T. Goodrich, M. J. Nelson, a nd J. Z. Sun, "The rainbow skip 
graph: a fault -tolerant constant -degree distributed data stru c-
ture," in Proceedings of the seventeenth annual ACM -SIAM sympo-
sium on Discrete algorithm , 2006, pp. 384 -393: Society for Indu s-
trial and Applied Mathematics. 
[46] G. T. Toussaint, "The relative neighbourhood graph of a finite 
planar set," Pattern recognition, vol. 12, no. 4, pp. 261-268, 1980. 
[47] B. Harwood and T. Drummond, "FANNG: fast approximate 
nearest neighbour graphs," in Proceedings of the IEEE Conference 
on Computer Vision and Pattern Recognition, 2016, pp. 5713-5722. 
[48] R. A. Dwyer, "Higher -dimensional Voronoi diagrams in linear 
expected time," Discrete & Computational Geometry, vol. 6, no. 3, 
pp. 343-367, 1991. 
AUTHOR ET AL.:  TITLE 13 
 
[49] L. Boytsov and B. Naidan, "Engineering Efficient and Effe ctive 
Non-metric Space Library," in Similarity Search and Applications: 
Springer, 2013, pp. 280-293. 
[50] L. Boytsov and B. Naidan, "Learning to prune in metric and 
non-metric spaces," in Advances in Neural Information Processing 
Systems, 2013, pp. 1574-1582. 
[51] A. Andoni and I. Razenshteyn, "Optimal Data -Dependent 
Hashing for Approximate Near Neighbors," presented at the 
Proceedings of the Forty-Seventh Annual ACM on Symposium 
on Theory of Computing, Portland, Oregon, USA, 2015. 
[52] J. Pennington, R. Socher, and C. D. Manning, "Glove: Global 
vectors for word representation," Proceedings of the Empiricial 
Methods in Natural Language Processing (EMNLP 2014), vol. 12, 
pp. 1532-1543, 2014. 
[53] P. Bolettieri  et al. , "CoPhIR: a test collection for content -based 
image retrieval," arXiv preprint arXiv:0905.4627, 2009. 
[54] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner, "Gradient-based 
learning applied to document recognition," Proceedings of the 
IEEE, vol. 86, no. 11, pp. 2278-2324, 1998. 
[55] E. Chávez, M. Graff, G. Navarro, and E. Téllez, "Near neighbor 
searching with K nearest references," Information Systems, vol. 
51, pp. 43-61, 2015. 
[56] E. C. Gonzalez, K. Figueroa, and G. Navarro, "Effective proxim-
ity retrieval by ordering permutations," Pattern Analysis and 
Machine Intelligence, IEEE Transactions on, vol. 30, no. 9, pp. 
1647-1658, 2008. 
[57] E. S. Tellez, E. Chávez, and G. Navarro, "Succinct nearest 
neighbor search," Information Systems, vol. 38, no. 7, pp. 1019 -
1030, 2013. 
[58] P. Sojka, "Software fram ework for topic modelling with large 
corpora," in In Proceedings of the LREC 2010 Workshop on New 
Challenges for NLP Frameworks, 2010: Citeseer. 
[59] C. Beecks, "Distance -based similarity models for content -based 
multimedia retrieval," Hochschulbibliothek der R heinisch-
Westfälischen Technischen Hochschule Aachen, 2013. 
 
 
Yury A. Malkov  received a Master’s degree in 
physics from Nizhny Novgorod State University in 
2009, and a PhD degree in laser physics from the 
Institute of Applied Physics RAS in 2015. He is 
author of 20+ papers on physics and computer 
science. Yury  currently occupies a position of a  
Project Leader in Samsung AI Center in Moscow . 
His current research interests include deep learn-
ing, scalable sim ilarity search, biological and 
artificial neural networks. 
 
Dmitry A. Yashunin received a Master’s degree 
in physics from Nizhny Novgorod State University 
in 2009, and a PhD degree in laser physics from 
the Institute of Applied Physics RAS in 2015. 
From 2008 to 2012 he was working in Mera Net-
works. He is author of 10+ papers on physics. 
Dmitry currently woks at Intelli-Vision in the pos i-
tion of  a leading research engineer. His current research interests 
include scalable similarity search, computer vision and deep learn-
ing. 
 
