Fast k-NN Search
\IEEEauthorblockN{ Ville Hyvönen\IEEEauthorrefmark{1
\IEEEauthorblockN{Sotiris Tasoulis
Efficient index structures for fast approximate nearest neighbor queries are required in many applications such as recommendation systems. In high-dimensional spaces, many conventional methods suffer from excessive usage of memory and slow response times. We propose a method where multiple random projection trees are combined by a novel voting scheme. The key idea is to exploit the redundancy in a large number of candidate sets obtained by independently generated random projections in order to reduce the number of expensive exact distance evaluations. The method is straightforward to implement using sparse projections which leads to a reduced memory footprint and fast index construction. Furthermore, it enables grouping of the required computations into big matrix multiplications, which leads to additional savings due to cache effects and low-level parallelization. We demonstrate by extensive experiments on a wide variety of data sets that the method is faster than existing partitioning tree or hashing based approaches, making it the fastest available technique on high accuracy levels. %Random projection trees have proven to be effective for approximate nearest neighbor searches in high-dimensional spaces where conventional methods are not applicable due to excessive usage of memory and computational time. We show that growing multiple trees on the same data can improve the performance even further, without increasing the total computational cost of queries when the tree parameters are controlled in a suitable way. Our experiments identify optimal parameter values to achieve accurate searches with extremely fast query times. In empirical experiments, the proposed algorithm achieves a better tradeoff between speed and accuracy than existing state-of-the-art methods. %\liang{I think we need to emphasise the voting in the abstract since it is an important contribution imo. Something like "We also propose a voting scheme to generalise the standard MPRT algorithm, which not only reduces computation complexity but also improves accuracy ..."
RP_trees_arxiv.tex
%% File list of work: IEEEtran.cls, IEEEtran_HOWTO.pdf, bare_adv.tex,
%%                    bare_conf.tex, bare_jrnl.tex, bare_jrnl_compsoc.tex
%%*************************************************************************


\documentclass[10pt, conference]{IEEEtran}



% *** CITATION PACKAGES ***
%
\usepackage{cite}
% cite.sty was written by Donald Arseneau
% V1.6 and later of IEEEtran pre-defines the format of the cite.sty package
% \cite{} output to follow that of IEEE. Loading the cite package will
% result in citation numbers being automatically sorted and properly
% "compressed/ranged". e.g., [1], [9], [2], [7], [5], [6] without using
% cite.sty will become [1], [2], [5]--[7], [9] using cite.sty. cite.sty's
% \cite will automatically add leading space, if needed. Use cite.sty's
% noadjust option (cite.sty V3.8 and later) if you want to turn this off.
% cite.sty is already installed on most LaTeX systems. Be sure and use
% version 4.0 (2003-05-27) and later if using hyperref.sty. cite.sty does
% not currently provide for hyperlinked citations.
% The latest version can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/cite/
% The documentation is contained in the cite.sty file itself.


% *** GRAPHICS RELATED PACKAGES ***
%
\ifCLASSINFOpdf
  \usepackage[pdftex]{graphicx}
  % declare the path(s) where your graphic files are
  % \graphicspath{{../pdf/}{../jpeg/}}
  % and their extensions so you won't have to specify these with
  % every instance of \includegraphics
  % \DeclareGraphicsExtensions{.pdf,.jpeg,.png}
\else
  % or other class option (dvipsone, dvipdf, if not using dvips). graphicx
  % will default to the driver specified in the system graphics.cfg if no
  % driver is specified.
  % \usepackage[dvips]{graphicx}
  % declare the path(s) where your graphic files are
  % \graphicspath{{../eps/}}
  % and their extensions so you won't have to specify these with
  % every instance of \includegraphics
  % \DeclareGraphicsExtensions{.eps}
\fi


\usepackage[cmex10]{amsmath}
\interdisplaylinepenalty=2500
\usepackage{array}
\usepackage{url}
\usepackage[utf8]{inputenc}
\usepackage{multicol}
\usepackage{multirow}
\usepackage{amssymb}
\usepackage{color}
\usepackage{balance}
\usepackage{booktabs}
\usepackage{dsfont}
\usepackage{subfig}
\usepackage{algorithm}
\usepackage[noend]{algpseudocode}




% *** SPECIALIZED LIST PACKAGES ***
\usepackage{algorithm}
\usepackage[noend]{algpseudocode}

% algorithmic.sty was written by Peter Williams and Rogerio Brito.
% This package provides an algorithmic environment fo describing algorithms.
% You can use the algorithmic environment in-text or within a figure
% environment to provide for a floating algorithm. Do NOT use the algorithm
% floating environment provided by algorithm.sty (by the same authors) or
% algorithm2e.sty (by Christophe Fiorio) as IEEE does not use dedicated
% algorithm float types and packages that provide these will not provide
% correct IEEE style captions. The latest version and documentation of
% algorithmic.sty can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/algorithms/
% There is also a support site at:
% http://algorithms.berlios.de/index.html
% Also of interest may be the (relatively newer and more customizable)
% algorithmicx.sty package by Szasz Janos:
% http://www.ctan.org/tex-archive/macros/latex/contrib/algorithmicx/




% correct bad hyphenation here
\hyphenation{op-tical net-works semi-conduc-tor}


\newcommand{\ve}[1]{\mathbf{#1}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\bigO}[1]{\mathcal{O}\left(#1\right)}
\def\independent{\perp\!\!\!\perp}


\newcommand{\teemu}{\textcolor[rgb]{1,0.5,0}}
\newcommand{\VH}{\textcolor[rgb]{.4,.2,.6}}
\newcommand{\sot}{\textcolor{green}}

\newcounter{LINumberOfComments}
\stepcounter{LINumberOfComments}

\newcommand{\liang}[1]{\textcolor{red}{\small \bf [LW\#\arabic{LINumberOfComments}\stepcounter{LINumberOfComments}: #1]}}


\begin{document}
%
% paper title
% can use linebreaks \\ within to get better formatting as desired
\title{Fast k-NN Search}



\author{\IEEEauthorblockN{
Ville Hyvönen\IEEEauthorrefmark{1}, 
Teemu Pitkänen\IEEEauthorrefmark{1},
Sotiris Tasoulis\IEEEauthorrefmark{2},
Elias Jääsaari\IEEEauthorrefmark{1},
Risto Tuomainen\IEEEauthorrefmark{1},\\
Liang Wang\IEEEauthorrefmark{3},
Jukka Corander\IEEEauthorrefmark{1}\IEEEauthorrefmark{4}\IEEEauthorrefmark{5} 
and Teemu Roos\IEEEauthorrefmark{1}
}
\IEEEauthorblockA{\IEEEauthorrefmark{1}Helsinki Institute for Information Technology HIIT, 
University of Helsinki, Finland. Email: \{firstname.lastname\}@cs.helsinki.fi}
\IEEEauthorblockA{\IEEEauthorrefmark{2}Liverpool John Moores University, UK. Email: S.Tasoulis@ljmu.ac.uk}
\IEEEauthorblockA{\IEEEauthorrefmark{3}Computer Laboratory, University of Cambridge,
UK. Email: liang.wang@cl.cam.ac.uk}
\IEEEauthorblockA{\IEEEauthorrefmark{4}Pathogen Genomics, Wellcome Trust Sanger Institute, Cambridge, UK}
\IEEEauthorblockA{\IEEEauthorrefmark{5}Department of Biostatistics, University of Oslo,  Norway. Email: jukka.corander@helsinki.fi}
}


%\author{\IEEEauthorblockN{Sotiris Tasoulis}
%\IEEEauthorblockA{Helsinki Institute for Information Technology HIIT\\
%Department of Computer Science \\
%University of Helsinki, Finland\\
%tasoulis@cs.helsinki.fi}
%
%\and \IEEEauthorblockN{Teemu Pitkänen}
%\IEEEauthorblockA{HIIT\\
%Department of Computer Science \\
%University of Helsinki, Finland\\
%teemu.pitkanen@helsinki.fi}
%
%\and \IEEEauthorblockN{Ville Hyvönen}
%\IEEEauthorblockA{HIIT\\
%Department of Computer Science \\
%University of Helsinki, Finland\\
%ville.o.hyvonen@helsinki.fi}
%}





% use for special paper notices
%\IEEEspecialpapernotice{(Invited Paper)}




% make the title area
\maketitle


\begin{abstract}
Efficient index structures for fast approximate nearest neighbor queries are required in many applications such as recommendation systems. In high-dimensional spaces, many conventional methods suffer from excessive usage of memory and slow response times. We propose a method where multiple random projection trees are combined by a novel voting scheme. The key idea is to exploit the redundancy in a large number of candidate sets obtained by independently generated random projections in order to reduce the number of expensive exact distance evaluations. The method is straightforward to implement using sparse projections which leads to a reduced memory footprint and fast index construction. Furthermore, it enables grouping of the required computations into big matrix multiplications, which leads to additional savings due to cache effects and low-level parallelization. We demonstrate by extensive experiments on a wide variety of data sets that the method is faster than existing partitioning tree or hashing based approaches, making it the fastest available technique on high accuracy levels.

%Random projection trees have proven to be effective for approximate nearest neighbor searches in high-dimensional spaces where conventional methods are not applicable due to excessive usage of memory and computational time. We show that growing multiple trees on the same data can improve the performance even further, without increasing the total computational cost of queries when the tree parameters are controlled in a suitable way. Our experiments identify optimal parameter values to achieve accurate searches with extremely fast query times. In empirical experiments, the proposed algorithm achieves a better tradeoff between speed and accuracy than existing state-of-the-art methods.
%\liang{I think we need to emphasise the voting in the abstract since it is an important contribution imo. Something like "We also propose a voting scheme to generalise the standard MPRT algorithm, which not only reduces computation complexity but also improves accuracy ..."}

\end{abstract}

\begin{IEEEkeywords}
Nearest Neighbor Search; Random Projections; High Dimensionality; Approximation Algorithms
\end{IEEEkeywords}


% For peer review papers, you can put extra information on the cover
% page as needed:
% \ifCLASSOPTIONpeerreview
% \begin{center} \bfseries EDICS Category: 3-BBND \end{center}
% \fi
%
% For peerreview papers, this IEEEtran command inserts a page break and
% creates the second title. It will be ignored for other modes.
\IEEEpeerreviewmaketitle

\section{Introduction}
\label{sec:intro}

Nearest neighbor search is an essential part of many machine learning algorithms. Often, it is also the most time-consuming stage of the procedure, see~\cite{Muja2014}. In applications areas, such as in recommendation systems, robotics and computer vision, where fast response times are critical, using brute force linear search is often not feasible. This problem is further magnified by the availability of inreasingly complex, high-dimensional data sources. 

Applications that require frequent $k$-NN queries from large data sets are for example object recognition \cite{Nister2006, Lowe2004}, shape recognition \cite{Amit1997} and image completion \cite{Hays2007} using large databases of image descriptors,  and content-based web recommendation systems \cite{Wang2016}.

Consequently, there is a vast body of literature on the algorithms for fast nearest neighbor search. These algorithms can be divided into exact and approximate nearest neighbor search. While  exact nearest neighbor search algorithms return the true nearest neighbors of the query point, they suffer from the curse of dimensionality: their performance degrades when the dimension of the data increases, rendering them no better than brute force linear search in the high-dimensional regime. 

In approximate nearest neighbor search, the main interest is the tradeoff between query time and accuracy, which can be measured either in terms of distance ratios or the probability of finding the true nearest neighbors. Several different approaches have been proposed and their implementations are commonly used in practical applications.

In this section, we first define the approximate nearest neighbor search problem, then review the existing approaches to it, and finally outline a new method that avoids the main drawbacks of existing methods.

%Finding the $k$ nearest neighbors (k-NN) of a query data object is an extremely frequent task in data mining. In the common formulation, we have a data set of $n$ points in the $d$-dimensional Euclidean space, and we look for the $k$ objects in this set with the smallest distance to a $d$-dimensional query vector. The information on neighbors can be used for numerous tasks, including classification, pattern recognition, information retrieval, and etc. Even though the exact nearest neighbors can be trivially found by \liang{"by a linear search, i.e. computing and sorting ..." because you mentioned linear search in the next paragraph ... } computing and sorting the distances between a query point to all $n$ data points, such operation takes $\bigO{dn + n\log n}$ time.

%In many modern applications, for example in the fields of computer vision, pattern recognition and natural language processing, both the dimensionality and the sample size of the data may be too large for a linear search to be feasible, especially in real-time applications where fast response times are critical. Examples of these high-dimensional data sets are image data where each pixel of the original picture is presented as one variable, and document-term matrices, even in a dimensionality-reduced form. Because exact nearest neighbor algorithms are not efficient for truly high-dimensional data sets, most of the applications where fast query times are required for high-dimensional data use approximate nearest neighbor search algorithms.

%sotiris: 
%In most modern applications, for example in the fields of computer vision, pattern recognition and natural language
%processing, both the dimensionality and the sample size of the data are too large for a linear search to be feasible,
%especially in real-time applications where fast response times are critical. 
%Examples of such high-dimensional data sets are image data where each pixel of the original picture
%is presented as one variable, and document-term matrices, even in a dimensionality-reduced form. In addition,
%exact nearest neighbor algorithms are not efficient for truly high-dimensional data sets, thus 
%in most applications where fast query times are critical the use of approximate nearest neighbor search algorithms
%is essential. 
% 
 
%In high dimensional spaces or data sets with a myriad of points such complexity is often unbearable. A recent trend is to deploy projections on random vectors and pre-built index structures in the search algorithm to improve time consumption. These methods do not always return the exact nearest neighbors, but objects that are nevertheless often close enough to the query object for practical applications.

%One of the most promising approaches to approximate nearest neighbor search is using random projection trees (RP trees). In RP trees the data points in the search set are repeatedly projected on random vectors and split into two smaller spaces, until some stopping criterion is met. These binary trees can be exploited in k-NN searches by placing the query object to a leaf using the same splitting criteria as while growing the tree, and then performing the linear time search only in the substantially smaller search set which consists of the data points at the same leaf. Previous studies have shown, that RP tree based k-NN approximations yield good results with significantly improved query times.



%sotiris:
%One of the most promising data structures for approximate nearest neighbor search in high-dimensional spaces are random projection trees (RP trees). They were first proposed for vector quantization in \cite{Dasgupta2009}, and later for $k$-NN search in \cite{Dasgupta2013}. Here we present our version of an RP tree, that especially optimized for approximate nearest neighbor search in high-dimensional space. By using the same random vector at each of the sibling node of the tree, as described in \ref{sec:growing}, we reduce the tree-growing time by one or two orders of magnitude, and also reduce the query times, while keeping the accuracy level constant.
%
%When growing an RP tree, the data set is recursively projected onto random vectors generated from the standard normal distribution and then split into two smaller sets, until a predefined maximum leaf size is met. These binary trees can be utilized in answering k-NN queries by routing a query point into a leaf by projecting it onto the same random vectors that where used when growing the tree, and then performing brute-force search only in the substantially smaller set consisting of the data points in a leaf the query point is routed into. 
%
%In high dimensions, the randomness in tree construction causes variability in the neighbours that are returned by several RP trees for a given query point. This can be taken advantage of by growing multiple RP trees and searching for the nearest neighbours in the combined search set. The trees can be constructed independently and in parallel exploiting parallel computing environments. The memory consumption is trivial since the optimized tree structures are very lightweight: at each node only a split point, a single floating-point value, is stored, and hence, a modest amount of memory is required even for hundreds of trees.

%In this paper we investigate the gains obtained by implementing multiple parallel RP trees. We provide both analytical and experimental results for choosing optimal parameter values such as the leaf size of the trees and the split criterion used in tree construction. The most important finding of our study is that a significant improvement in the accuracy of nearest neighbor queries can be achieved by using multiple RP trees \emph{even when the total computational cost,
%summed up over all parallel computation nodes, is kept fixed.} In other words, when the workload is distributed over
%multiple RP trees, the accuracy of the nearest neighbor search is actually improved compared to using a single RP tree.

%\liang{Again, we shouldn't be so modest when it comes to contribution. The current version reads like we only did some empirical evaluation which is not true. Be very specific about our contributions, e.g., we performed a theoretical/complexity analysis of mrpt; we generalise the mprt using voting scheme; we show how voting scheme is able to exploit the frequency information contained in trees then further reduce the complexity; we further increase the parallelism by using one projection per level ... }

%\liang{Do not talk like we only make marginal contribution to an existing data structure or algorithm, talk like it is a new thing ;-) so, it is better to end intro section with some explicit numbers like "the proposed mrpt algorithm improves the accuracy by ???, reduces the complexity by ???, increases the parallelism by ???. You can of course merge this into contribution part." }

%\liang{IF there is space, add "the paper is organised as follows ..." in order to help reviewers grasp the bird-view of what we do in the paper ... of course given we are not running out of space.}

%\section{Background and related work}
%\label{sec:related}

%There are many efficient algorithms that first build an index data structure and then utilize it to make fast queries for nearest neighbors of a query point. For many of these methods the time complexity for structure-building has an exponential dependency on the dimension $d$; for example the time complexity of Bentley's divide and conquer-method \cite{Bentley1980} is $\bigO{n(\log n)^{d-1}}$ and in \cite{Clarkson1983} an algorithm which takes $\bigO{c^d \log n}$ time for some constant $c$ is presented. The exponential dependency on the dimension means that these methods are not feasible for even moderately high-dimensional data sets. In moderately high-dimensional data, methods such as  metric trees \cite{Uhlmann1991} and ball-trees \cite{Omohundro1987}, are commonly used but even they tend to fail in large, very high-dimensional data sets.

%However there are also methods which are intented specifically for nearest neighbor search in high-dimensional spaces.  For many problems it is not necessary to return the true nearest neighbors of a query point, but it is sufficient that the points returned are close enough to the query point, or at least close to the query point with a high enough probability. For example in Kvasir \cite{Wang2015} it is not necessary to return the web pages that are the most similar to the web page just visited, but the user is most likely satisfied with results that are close enough to the subject matter. 

%sotiris:
%However, there are also methods which are intended specifically for nearest neighbor search in high-dimensional spaces.  For many problems it is not necessary to return the true nearest neighbors of a query point, but it is sufficient that the points returned are close enough to the query point, or at least close to the query point with a high enough probability. For example, in recommendation systems such as Kvasir \cite{Wang2016} it is not necessary to return the web pages that are exactly most similar to the viewed web page, but the user is often satisfied with pages whose topic is close enough.

%The tradeoff between the search efficiency and search accuracy explains why most of the modern methods intended for fast nearest neighbor search in high-dimensional spaces do \textit{approximate} instead of \textit{exact} nearest neighbor search. One formulation is $(1 + \epsilon)$--approximate k-NN search \cite{Arya1998}, in which the distance of the points returned by the algorithm is at most $(1 + \epsilon)$ times the distance of the true $k$th nearest neighbor of the query point. However, herein we consider a slightly different version of the problem, where $k$ points are returned and it is required that those points are among the $k$ nearest neighbors with a high probability, or in other words, that a large fraction of the $k$ nearest neghbors are successfully returned. \liang{I also recommend to give another couple of examples where approximate knn applies, otherwise people will think this is an extension of kvasir paper. E.g., Spotify is another example wherein approx. knn is used to search similar songs, they use a data set of about 300 dimensions I remember.}

\subsection{Approximate nearest neighbor search}

%The problem we consider here is approximate nearest neighbor search. One formulation is $(1 + \epsilon)$--approximate $k$-NN search \cite{Arya1998}, in which the distance of the points returned by the algorithm is at most $(1 + \epsilon)$ times the distance of the true $k$th nearest neighbor of the query point. However,  in this article we consider a slightly different version of the problem: we do not set any error bounds, but assess the efficiency of the algorithms empirically by examining the $k$-NN query time required to reach a given accuracy level.


Consider a metric space $\mathcal{M}$, a data set $\ve{X} = (\ve{x}_1, \dots, \ve{x}_n) \subseteq \mathcal{M}$, a query point $\ve{q} \in \mathcal{M}$, and a distance metric $m: \mathcal{M}^2 \rightarrow \R$. The $k$-nearest neighbor ($k$-NN) search is the task of finding the $k$ closest (w.r.t. $m$) points to $\ve{q}$ from the data set $\ve{X}$, i.e., find a set $\ve{K} \subseteq \ve{X}$ for which it holds that $|\ve{K}| = k$ and
\[
m(\ve{q},\ve{x}) \leq m(\ve{q},\ve{y})
\]
for all $\ve{x} \in \ve{K}$, $\ve{y} \in \ve{X} \setminus \ve{K}$.

In this work, we consider applications where $\mathcal{M}$ is the $d$-dimensional Euclidean space $\R^d$, and $m$ is Euclidean distance 
\[
m(\ve{x},\ve{y}) = \sqrt{\sum_{i=1}^d (x_i-y_i)^2}.
\]

An approximate nearest neighbor search algorithm can be divided into two separate phases: an offline phase, and an online phase. In the offline phase an index is built, and the online phase refers to completing fast nearest neighbor queries using the constructed index. In practical applications, the most critical considerations are the accuracy of the approximation (the definition of the accuracy and the required accuracy level depending on the application), and the query time needed to reach it.

As a measure of accuracy, we use recall, which is relevant especially for information retrieval and recommendation settings. It is defined as the proportion of true $k$-nearest neighbors of the query point returned by the algorithm:
\[
\mathrm{Recall} = \frac{|\ve{A} \cap \ve{K}|}{k}.
\]
Here $\ve{A}$ is a set of $k$ approximate nearest neighbors returned by the algorithm, and $\ve{K}$ is the set of true $k$ nearest neighbors of the query point.

Of secondary importance are the index construction time in the offline stage and the memory requirement, which are usually not visible to the user but must be feasible. Note that in many applications the index must be updated regularly when we want to add new data points into the set from where the nearest neighbors are searched from. 

\subsection{Related work}

Most effective methods for approximate nearest neighbor search in high-dimensional spaces can be classified into either hashing, graph, or space-partitioning tree based strategies.

\subsubsection{Hashing based algorithms}

The most well-known and effective hashing based algorithms are variants of locality-sensitive hashing (LSH) \cite{Indyk1998, Gionis1999, Andoni2006}. In LSH, several hash functions of the same type are used. These hash functions are \emph{locality-sensitive}, which means that nearby points fall into the same hash bucket with higher probability than the points that are far away from each other.  When answering a $k$-NN query, the query point is hashed with all the hash functions, and then a linear search is performed in the set of data points in the buckets where the query point falls into. Multi-probe LSH \cite{Lv2007, Dong2008, Andoni2015} improves the efficiency of LSH by searching also the nearby hash buckets of the bucket a query point is hashed into. Thus a smaller amount of hash tables is required to obtain a certain accuracy level. The choice of the hash function affects the performance and needs to be done carefully in each case. 

\subsubsection{Graph based algorithms}

Approximate graph based methods such as \cite{Hajebi2011, Dong2011, Wang2012} build a $k$-NN-graph of the data set in an offline phase, and then utilize it to do fast $k$-NN queries. For example in \cite{Wang2012} several graphs of random subsets of the data are used to approximate the exact $k$-NN-graph of the whole data set. The data set is divided hierarchically and randomly into subsets, and the $k$-NN-graphs of these subsets are built. This process is then repeated several times, and the resulting graphs are merged and the final graph is utilized to answer nearest neighbor queries. Graph based methods are in general efficient, but the index construction is slow for larger data sets.

\subsubsection{Space-partitioning trees}

The first space partitioning-tree based strategy proposed for nearest neighbor search was the $k$-d tree \cite{Bentley1975}, which divides the data set hierarchically into cells that are aligned with the coordinate axes. Nearest neighbors are then searched by performing a backtracking or priority search in this tree. $k$-d trees and other space-partitioning trees can be utilized in approximate nearest neighbor search by terminating the search when a predefined number of data points is checked, or all the points within some predefined error bound from the query point are checked \cite{Arya1998}. 

Instead of the hyperplanes, clustering algorithms can be utilized to partition the space hierarchically; for example the $k$-means algorithm is used in $k$-means trees \cite{Fukunaga1975, Muja2014}. Other examples of clustering based partitioning trees are cover trees \cite{Beygelzimer2006}, VP trees \cite{Yianilos1993} and ball trees \cite{Leibe2006}. Similarly to graph-based algorithms, clustering based variants have the drawback of long index construction times due to hierarchical clustering on large data sets.

An efficient scheme to increase the accuracy of space-partitioning trees is to utilize several parallel randomized space-partitioning trees. Randomized $k$-d trees \cite{Silpa2008, Muja2009} are grown by choosing a split direction at random from dimensions of the data in which it has the highest variance. Nearest neighbor queries are answered using priority search; a single priority queue, which is ordered by distance of the query point from splitting point in the splitting direction, is maintained for all the trees. 

Another variant of randomized space-partitioning tree is the random projection tree (RP tree), which was first proposed for vector quantization in \cite{Dasgupta2009}, and later modified for approximate nearest neighbor search in \cite{Dasgupta2013}. In random projection trees the splitting hyperplanes are aligned with the random directions sampled from the unit sphere instead of the coordinate axes. Nearest neighbor queries are answered using \emph{defeatist search}: first the query point is routed down in several trees, and then a brute force linear search is performed in the union of the points of all the leaves the query point fell into.

%However, we introduce here both a more efficient way to randomize space-partitioning trees, and a more efficient way utilize them in a nearest neighbor search.

%A promising approach to the nearest neighbor search is the use of random projection trees (RP trees). They were first proposed for vector quantization in \cite{Dasgupta2009}, and later modified for approximate nearest neighbor search in \cite{Dasgupta2013}. Random projection tree is a form of space-partitioning tree resembling the classic $k$-d tree. In random projection trees the splitting hyperplanes are aligned with the random directions sampled from the unit sphere instead of the coordinate axes. 

%The intuition behind using random projections is that projecting the data onto a lower dimensional subspace by multiplying it by on orthogonal random matrix can be used for dimensionality reduction as effectively as PCA, but with much lower computational cost \cite{Bingham2001, Achlioptas2001}. On the other hand the randomness inherent in the tree growing process enables growing many independent trees with different sets of random vectors and combining the results; this can be utilized to reach very high accuracy when a sufficiently high number of trees are used. A strong point of random projection trees is also their relative simplicity: they have only two tuning parameters, number and depth of the trees used, and the method is quite robust with respect to these parameters.

\subsection{MRPT algorithm}

One of the strengths of randomized space-partitioning trees is in their relative simplicity compared to hashing and graph based methods: they have only few tunable parameters, and their performance is quite robust with respect to them.
They are also perfectly suitable for parallel implementation because the trees are independent, and thus can be easily parallelized either locally or over the network with minimal communication overhead. 

However, both of the aforementioned variants of randomized space-partitioning trees have a few weak points which limit their efficiency and scalability.  First, randomization used in randomized $k$-d trees does not make the trees sufficiently decorrelated to reach high levels of accuracy efficiently. Randomization used in RP trees is sufficient, but it comes at high computational cost: computation of random projections is time-consuming for high-dimensional data sets because the random vectors have the same dimension as the data. 

Second, both the priority queue search used in \cite{Silpa2008, Muja2009} and the defeatist search used in \cite{Dasgupta2013} require checking a large proportion of the data points to reach high accuracy levels; this leads to slow query times. Third, because each node of an RP tree has its own random vector, the number of random vectors required is exponential with respect to tree depth. In the high-dimensional case this on the one hand increases the memory required by the index unnecessarily, and on the other hand slows down the index construction process.

Because of the first two reasons the query times of the randomized space-partitioning trees are not competitive with the fastest graph based methods (cf. experimental results in  Section \ref{sec:results}).

In this article we propose a method that uses multiple random projection trees (MRPT); it incorporates two additional features that lead to fast query times and accurate results. The algorithm has the aforementioned strengths of randomized space-partitioning trees but avoids their drawbacks. 

More specifically, our contributions are:

\begin{enumerate}
\item We show that sparse random projections can be used to obtain a fast variant of randomized space-partitioning trees.

\item We propose the MRPT algorithm that combines multiple trees by \emph{voting search} as a new and more efficient method to utilize randomized space-partitioning trees for approximate nearest neighbor search.

\item We present a time and space complexity analysis of our algorithm.

\item We demonstrate experimentally  that the MRPT algorithm with sparse projections and voting search outperforms the state-of-the-art methods using several real-world data sets across a wide range of sample size and dimensionality.

%\item We present the multiple random projection tree (MRPT) algorithm that is optimized for high dimensionality and large sample sizes, and which can be easily and efficiently parallelized.
%\item We show experimentally how increasing the number of trees increases the number of true nearest neighbors found, even while the size of the set where a final linear search is performed stays constant or even decreases.
%\item We demonstrate the efficiency of the algorithm for approximate nearest neighbor search by comparing its performance to the state-of-the-art methods. 
\end{enumerate} 




%An efficient space partitioning-tree based scheme of Muja and Lowe \cite{Muja2009, Muja2014} utilizes randomised k-d trees, which are grown by choosing a split direction at random among the five dimensions along which the data set has a highest variance. Nearest neighbor queries are answered using priority search; a single priority queue, that is ordered by distance of the query point from splitting point in the splitting direction, is maintained for all the trees. 
%{\bf what are the drawbacks with this approach?}


%Each of the existing approaches has limits in terms of the tradeoff between search accuracy and query time and other performance characteristics. High-dimensional data, high accuracy, and fast response times are still critical for many practical applications. In this paper, we propose methods that consistently achieve better performance than the aforementioned approaches especially in the high dimensional, high accuracy regime.

%\section{Random projection trees}
%\label{sec:RP_tree}

In the following sections we describe the proposed method by first revisiting the classic RP tree method~\cite{Dasgupta2013}, and then describing the novel features that lead to increased accuracy and reduced memory footprint and query time.

\section{Index construction}
\label{sec:growing}

\subsection{Classic random projection trees}

% An RP tree is a form of space-partitioning tree. It resembles a k-d tree, but the splitting hyperplanes are aligned with random directions instead of the coordinate axes. 

%We construct the index for the MRPT algorithm by growing $T$ different RP trees as follows.

At the root node of the tree, the points are projected onto a random vector $\ve{r}$ generated from the $d$-dimensional standard normal distribution $N_d(\ve{0}, \ve{I})$. The data set is then split into two child nodes using the median of the points in the projected space: points whose projected values are less or equal to the median in the projected space are routed into the left child node, and points with projected values greater than the median are routed into the right child node.  So the points 
\[
\ve{x}^{(1)}, \dots , \ve{x}^{(\lceil n/2 \rceil)}
\]
fall into the left branch, and the points
\[ 
\ve{x}^{(\lceil n/2 \rceil + 1)} , \dots , \ve{x}^{(n)} 
\]
fall into the right branch. Here $\ve{x}^{(1)}, \dots , \ve{x}^{(n)}$ denotes the ordering of the points $\ve{x}_1, \dots , \ve{x}_n$  with respect to the values of their projections $\ve{x}_1^T \ve{r}, \dots , \ve{x}_n^T \ve{r}$. This process is then repeated recursively for the subsets of the data in both of the child nodes, until a predefined depth $\ell$ is met. 

%{\bf TR: strictly speaking, we claim later that the number of points in different leaf nodes differs at most by one --- however, the above strategy where the left child always gets at least as many points as the right child, leads to a bit more imbalanced distribution of the points. if we do it as above in the implementation, we should in principle modify the computational complexity analysis... however, it may be just be better to just forget this as it is really a minor issue.} 

For high-dimensional data, the computation of random projections is slow, and the memory requirement for storing the random vectors is high. To reduce both the computation time and the memory footprint, we propose two improvements: 
\begin{itemize}
\item Instead of using dense vectors sampled from the $d$-dimensional normal distribution, we use sparse vectors to reduce both storage and computation complexity.
\item Instead of using different projection vectors for each intermediate node of a tree, we use one projection vector for all the nodes on the same level of a tree, so that we can further reduce the storage requirement and maximize low-level parallelism through vectorization.
\end{itemize}

In addition, instead of splitting at a fractile point chosen at random from the interval $[\frac{1}{4}, \frac{3}{4}]$, as suggested in~\cite{Dasgupta2013}, we split at the median to make the performance more predictable and to enable saving the trees in a more compact format. 

In the following subsections we briefly discuss the details of the above improvements. 
Pseudocode for constructing the index using sparse RP trees is given in detail in Algorithms \ref{algo:1}--\ref{algo:2} below.
The proposed Algorithm \ref{algo:1} consists of three embedded \textbf{for} loops. The outermost loop (line 4 - 15) builds $T$ RP trees by continuously calling the \texttt{GROW\_TREE} function in Algorithm \ref{algo:2}. For each individual RP tree, the two inner \textbf{for} loops (line 6 - 12) prepare a sparse random matrix $\ve{R}$ which will be used for projecting the data set $\ve{X}$ to $\ve{P}$ using $\ell$ random vectors (at line 13). Specifically, the innermost loop is for constructing a $d$-dimensional sparse vector for a given tree level.
In Algorithm \ref{algo:2}, the \texttt{GROW\_TREE} function constructs a binary search tree by recursively splitting the high-dimensional space $\ve{X}$ into sub-spaces using the previous projection $\ve{P}$.


\begin{algorithm}
\begin{algorithmic}[1]
\Function{grow\_trees}{$\ve{X}$, $T$, $\ell$, $a$}
  \State n $\leftarrow$ $\ve{X}$.nrows
  \State let trees[$1 \dots T$] be a new array
  \For{$t$ in $1, \dots, T$}
    %\State generate $\ve{r}_{t,1}, \dots, \ve{r}_{t,\ell}$ from $N_d(\ve{0}, \ve{I})$
    \State let $\ve{R}$ be a new $d \times \ell$ matrix
	\For{level in $1, \dots , \ell$}
		\For{$i$ in $1, \dots , d$}
			\State generate $z$ from Bernoulli($a$)
			\If{$z=1$}
            \State generate $\ve{R}$[$i$, level] from $N(0,1)$
			\Else
            \State  $\ve{R}$[$i$, level] $\leftarrow 0$ 
			\EndIf
		\EndFor
	\EndFor    
    \State $\ve{P}$ $\leftarrow$ $\ve{X} \ve{R}$
    \State trees[$t$].root $\leftarrow$ \Call{grow\_tree}{$\ve{X}$, [$1 \dots n$], 0, $\ve{P}$}
    \State trees[$t$].random\_matrix $\leftarrow$ $\ve{R}$
  \EndFor
  \State \textbf{return} trees
\EndFunction
\end{algorithmic}
\caption{Grow $T$ RP trees.}
\label{algo:1}
\end{algorithm}


\begin{algorithm}
\begin{algorithmic}[1]
    \Function{grow\_tree}{$\ve{X}$, indices, level, $\ve{P}$}
  \If{level $= \ell$}  
  \State return $\ve{X}$[indices, ] as leaf
  \EndIf
  \State proj $\leftarrow$ $\ve{P}$[indices, level]
  \State split $\leftarrow$ median(proj)
  \State left\_indices $\leftarrow$ indices[proj $\leq$ split]
  \State right\_indices $\leftarrow$ indices[proj $>$ split]
  \State left $\leftarrow$ \Call{grow\_tree}{$\ve{X}$, left\_indices, level + 1, $\ve{P}$}
  \State right $\leftarrow$ \Call{grow\_tree}{$\ve{X}$, right\_indices, level + 1, $\ve{P}$}
  \State \textbf{return} split, left, right
\EndFunction
\end{algorithmic}
\caption{Grow a single RP tree.}
\label{algo:2}
\end{algorithm}



% \begin{algorithm}
% \begin{algorithmic}[1]
% \Function{grow\_trees}{$\ve{X}$, $T$, $\ell$, $a$}
%   \State let trees[$1 \dots T$] be a new array
%   \For{$t$ in $1, \dots, T$}
%     %\State generate $\ve{r}_{t,1}, \dots, \ve{r}_{t,\ell}$ from $N_d(\ve{0}, \ve{I})$
% 	\For{level in $1, \dots , \ell$}
% 		\For{$i$ in $1, \dots , d$}
% 			\State generate $z$ from Bernoulli($a$)
% 			\If{$z=1$}
% 				\State generate $\ve{r}_{t,\mathrm{level},i}$ from $N(0,1)$
% 			\Else
% 				\State  $\ve{r}_{t,\mathrm{level},i} \leftarrow 0$ 
% 			\EndIf
% 		\EndFor
% 	\EndFor    
%     \State trees[$t$] $\leftarrow$ \Call{grow\_tree}{$\ve{X}$, 0, $\ve{r}_t$}
%   \EndFor
%   \State \textbf{return} trees
% \EndFunction
% \end{algorithmic}
% \caption{Grow $T$ RP trees.}
% \label{algo:1}
% \end{algorithm}


% \begin{algorithm}
% \begin{algorithmic}[1]
% \Function{grow\_tree}{$\ve{X}$, level, $\ve{r}_t$}
%   \If{level $= \ell$}  
%     \State return $\ve{X}$ as leaf
%   \EndIf
%   \State proj $\leftarrow \ve{X} \ve{r}_{t,\mathrm{level}}$
%   \State split $\leftarrow$ median(proj)
%   \State left $\leftarrow$ \Call{grow\_tree}{$\ve{X}$[proj $\leq$ split, ], level + 1, $\ve{r}_t$}
%   \State right $\leftarrow$ \Call{grow\_tree}{$\ve{X}$[proj $>$ split, ], level + 1, $\ve{r}_t$}
%   \State \textbf{return} split, left, right
% \EndFunction
% \end{algorithmic}
% \caption{Grow a single RP tree.}
% \label{algo:2}
% \end{algorithm}


\subsection{Sparse random projections}

With high-dimensional data sets, computation of the random vectors easily becomes a bottleneck on the performance of the algorithm. However, it is not necessary to use random vectors sampled from the $d$-dimensional standard normal distribution to approximately preserve the pairwise distances between the data points. Achlioptas  \cite{Achlioptas2001} shows that the approximately distance-preserving low-dimensional embedding of Johnson-Lindenstrauss-lemma is obtained also with sparse random vectors with components sampled from $\{-1, 0, 1\}$ with respective probabilities $\{\frac{1}{6}, \frac{2}{3}, \frac{1}{6}\}$. Li \emph{et al.} \cite{Li2006} prove that the same components with respective probabilities $\{\frac{1}{2\sqrt{d}}, 1-\frac{1}{\sqrt{d}}, \frac{1}{2\sqrt{d}}\}$, where $d$ is the dimension of the data, can be used to obtain a $\sqrt{d}$-fold speed-up without significant loss in accuracy compared to using normally distributed random vectors. 

We use sparse random vectors $\ve{r} = (r_1, \dots , r_d)$, whose components are sampled from the standard normal distribution with probability $a$, and are zeros with probability $1-a$ :
\[
r_i = 
\begin{cases}
N(0,1) \quad &\text{with probability}\,\, a \\
0 \quad &\text{with probability}\,\, 1-a.
\end{cases}
\]

The sparsity parameter $a$ can be tuned to optimize performance but we have observed that $a = \frac{1}{\sqrt{d}}$ recommended in \cite{Li2006} tends to give near-optimal results in all the data sets we tested, which suggests that further fine-tuning of this parameter is unnecessary. This proportion is small enough to provide significant computational savings through the use of sparse matrix libraries. 

\subsection{Compactness and speed with fewer vectors}

In classic RP trees, a different random vector is used at each inner node of a tree, whereas we use the same random vector for all the sibling nodes of a tree. This choice does not affect the accuracy at all, because a query point is routed down each of the trees only once; hence, the query point is projected onto a random vector $\ve{r}_i$ sampled from the same distribution at each level of a tree. This means that the query point is projected onto i.i.d.\ random vectors $\ve{r_1, \dots , r_\ell}$ in both scenarios. 

An RP tree has $2^\ell - 1$ inner nodes; therefore, if each node of a tree had a different random vector as in classic RP trees, $2^\ell - 1$ different random vectors would be required for one tree. However, when a single vector is used on each level, only $\ell$ vectors are required. This reduces the amount of memory required by the random vectors from exponential to linear with respect to the depth of the trees.

Having only $\ell$ random vectors in one tree also speeds up the index construction significantly. While some of the observed speed-up is explained by a decreased amount of the random vectors that have to be generated, mostly it is due to enabling the computation of all the projections of the tree in one matrix multiplication: the projected data set $\ve{P} \in \R^{n \times \ell}$ can be computed from the data set $\ve{X} \in \R^{n \times d}$ and a random matrix $\ve{R} \in \R^{d \times \ell} $ as 
\[
\ve{P} = \ve{X} \ve{R}.
\]
Although the total amount of computation stays the same, in practice this speeds up the index construction significantly due to the cache effects and low-level parallelization through vectorization.

%\liang{This section is very clear to me, but it is still good to use figures to illustrate the mechanisms you proposed. Actually, it would be nice to use one figure for each smart tweak you made. E.g., refer to p.11 and p.12 in my presentation slides at http://www.cl.cam.ac.uk/~lw525/publications/high-dimensional-bigdata-cambridge-workshop-slides.pdf ... feel free to use any figures you need. It is already the fourth page of the paper, we need some figures to "entertain" our reviewers :)}


%Using the same random vector for all the sibling nodes speeds up the index construction significantly: when using a different random vector for each of the nodes at one level of a tree, as many random vectors as there are internal nodes in the tree has to be generated. However, when using the same random vector for all the sibling nodes, only as many random vectors as there are levels in the tree has to be generated. When a median split is used, an RP tree is an almost complete binary tree, so for the tree of depth $\ell = \lceil \log \frac{n}{n_0} \rceil$ the amount $l$ of random vectors needed is 
%\[
%\frac{n}{2n_0} \leq 2^{\ell-1} \leq l \leq 2^{\ell} - 1 \leq \frac{2n}{n_0}
%\]       
%in the first case. In the second case only
%\[
%l = \ell = \lceil \log \frac{n}{n_0} \rceil
%\]
%random vectors are needed, so the dependence on the size of the data set is logarithmic instead of linear.


%When the dimension of the data set is not ultra high, using the same vectors also makes it possible to save random vectors instead of just the random seeds, because a tree contains only $\ell \leq \log n $ random vectors. This saves the overhead of generating random vectors from the seed when answering the query, and reduces the query time at a the price of an insignificant memory overhead.

%\subsection{Split criterion}

%A third difference to the RP trees of Dasgupta and Sinha \cite{Dasgupta2013} is a choice of split criterion. They split at the fractile point chosen uniformly at random from the range $[\frac{1}{4}, \frac{3}{4}]$ to introduce extra randomness into the tree growing process. We also tested other split criteria, such as splitting at the mean of the projected data, a completely random split, and splitting at the longest interval between two data points at the midrange of the projected data. However, since differences between the performance of the different split criteria were not significant, we use the median split because it makes the trees complete binary trees. This enables saving trees in a compact form using a preallocated array, because the size of the trees is known when a predefined depth $\ell$ is fixed. The query time is also more predictable, because each leaf of each tree lies at the same depth, and contains a same ($\pm 1$) number of data points. 

\subsection{Time and space complexity: Index construction}

At each level of each tree the whole data set is projected onto a $d$-dimensional random vector that has on average $ad$ non-zero components, so the expected\footnote{The following complexity results hold exactly, if the algorithm is modified so that each random vector has exactly $\lceil ad \rceil$ non-zero components instead of the expected number of non-zero components being $ad$.} index construction time is $\Theta(T\ell n ad)$. For classic RP trees ($a=1$) this is $\Theta(T\ell n d)$, but for sparse trees ($a = \frac{1}{\sqrt{d}}$) this is only $\Theta(T\ell n \sqrt{d})$. 

For each tree, we need to store the points allocated to each leaf node, so the memory required by the index is at least $\bigO{Tn}$. At each node only a split point is saved; this does not increase the space complexity because there are only $2^\ell - 1 < n$ inner nodes in one tree.

The expected amount of memory required  by one random vector is $\Theta(ad)$ when random vectors are saved in sparse matrix form. This is $\Theta(d)$ for dense RP trees, and $\Theta(\sqrt{d})$ for sparse RP trees with the sparsity parameter fixed to $a=\frac{1}{\sqrt{d}}$. Because an RP tree has $2^\ell - 1$ inner nodes, the memory requirement for $T$ classic RP trees, which have a different random vector for each node of a tree, is $\bigO{Tdn}$. However, in our version, in which a single vector is used on each level, there are only $\ell$ vectors; hence, the memory requirement for $T$ sparse RP trees is $\bigO{T(\sqrt{d}\log n + n)}$.  

%\subsection{Tuning the parameters}
%The performance of the MRPT algorithm is controlled by three parameters: the number of trees $T$, the depth of the trees $\ell$ and the vote threshold $v$. Using more trees leads to better recall at a cost of higher query time, as does decreasing the depth of the trees or the vote threshold. 


%While the number of trees is the primary way to tune the performance, the choice of the depth and the vote threshold are also important for obtaining optimal performance. 

% Especially in the case where memory usage is of concern, decreasing the depth is preferable to increasing the number of trees.

%Generally, several hundred trees are needed to obtain high recall rates. It is a good idea to start with a relatively low vote threshold, for example to set it as 1\% of the number of trees. 

%The vote count is sensitive to the depth of the trees and should be decreased when the depth of the trees is increased. It should also be slightly increased when the number of trees is higher. To keep accuracy intact as the number of data points grows, in general the number of trees and the vote limit should be increased and the depth of the trees decreased.

%As an example, on a data set with $n \approx 10^5$ and $d = 10^3-10^4$, to get good performance on 90\% recall, a common choice of parameters is to use around 500-750 trees with 7-10 depth and 5-15 as the vote threshold. If the dimension is lower, that is, $d = 10^2 - 10^3$, a good starting point is around 250-500 trees, 9-12 depth and 3-5 as the vote threshold.

%On harder data sets where the variance in distance between points is small (that is, the data set has a high intrinsic dimensionality ~\cite{Amsaleg2015}), it is a good idea to decrease the depth of the trees while increasing the vote limit. This is because on harder data sets the true nearest neighbors are less likely to occupy the same leaf as the query point if the leaves are small.

%To tune the parameters automatically, a parameter space exploration scheme similar to the one proposed by Muja et al. \cite{Muja2009} could be adopted. We leave this as future work.


% When the sparsity parameter is fixed to $a=\frac{1}{\sqrt{d}}$, this is $\bigO{T(\sqrt{d}\log n + n)}.$


%The expected total number of non-zero elements in the projection vectors in all the trees is $Tad\ell$. Since for sparse RP trees $\ell = \lceil \log(n) \rceil$ and $a = \frac{1}{\sqrt{d}}$, this equals $T\sqrt{d} \lceil \log(n) \rceil$. 


% Unless the data set has an ultra-high dimensionality (or has a small sample size, in which case the memory required is insignificant anyway), this is in practice almost always $\bigO{Tn}$, which means that the storage requirement for the projection vectors is dominated by the storage required by the allocation of the points to the leaf nodes. We mention that the storage of the  projection vectors can be avoided altogether by saving only the random seed, from which the random vectors used in the index construction can be generated, in the query phase. However, this tends to slightly increase the query times; hence we save the random vectors as default.  


%\subsection{Splitting criteria}
%
%\VH{Remove this chapter (Splitting criteria) and add more experiments / something else instead?}
%\sot{I agree. This part need an extensive analysis to be consistent and is out of scope at the moment ( try some ideas on the density approach for example). We can still include a sort reference if it fits somewhere maybe as future work that could further improve the method.}
%
%Another difference between our version of RP tree and RP tree of Dasgupta and Sinha is the choice of split criterion. Dasgupta and Sinha \cite{Dasgupta2013} suggest to introduce more randomness into a tree growing process by splitting from a fractile point chosen uniformly at random from the range $[\frac{1}{4}, \frac{3}{4}]$ instead of splitting at the median. Wang \emph{et al.} \cite{Wang2015} use the mean instead of the median. In addition to these, we introduce two new split criteria to test how the choice of the split criterion affects the accuracy of the results. We tested the following split criteria:
%
%\begin{itemize}
%\item Mean split:  The split point is determined by the mean of the projections.
%\item Median split: The split point is determined by the median of the projections.
%\item Random split: The split is performed at a random point chosen uniformly from the interval between the minimum and the maximum of the projections.
%\item Bounded random split:  Same as random split, except the split point is bounded to be chosen from the middle half of the interval, i.e. both sides of the split correspond to at least 25\% of the whole interval.
%\item Longest interval split:  The projections are first sorted and the split is then performed at the longest interval between two consecutive projections. To avoid extremely uneven splits, we only consider the intervals between the middle 50\% of the projections, thus both sides of the split have at least 25\% of the data objects.
%\end{itemize}
%
%The longest interval split is used as a very rough approximation of a minimum density split used in random projection based hierarchical clustering \cite{Tasoulis2013, Tasoulis2014}. It is done in a middle half, because otherwise a split point would almost always be at the tail of the distribution (of the projected data). A motivation is to find a low density area between two high density areas, and thereby split with higher probability between the natural clusters of the data set. Unlike in the case of clustering, when doing nearest neighbor search it is not clear whether splitting should be done at the high density area (mean, median) or at the low density area (longest interval split); the longest interval split is introduced to test this empirically. A completely random split is implemented to test if introducing even more randomness into a tree growing process improves the accuracy by making the trees as uncorrelated as possible. 
%
%We compared the performance of the different split criteria experimentally. However, because there emerged no clear winner, and it seems that the differences between the different split criteria were explained mostly by the random variation, we omit the results. The only expection was the median split, which seems to give consistently slightly worse results than the other split criteria, especially with the higher number of trees.
%
%However, the median split is a good choice because of its other desirable properties: the trees are always balanced when using the median split, so the memory for them can be preallocated before growing the tree, and all the matrix multiplications required for the random projections can be computed beforehand. On the other hand for example the totally random split, which seemed to give quite a nice accuracy, can lead to very imbalanced trees, and thus very unpredictable performance.
%
%Because the differences on the accuracy of the different split criteria seem negligible, the median split remains our choice out of the convenience of implementation (which leads to the reduced growing and query times in practice), and the more predictable performance.

\section{Query phase}

In many approximate nearest neighbor search algorithms the query phase is further divided into two steps: a candidate generation step, and an exact search step. 

In the candidate generation step, a candidate set $S$, for which usually $|S| \ll n$, is retrieved from the whole data set, and then in the exact search step $k$ approximate nearest neighbors of a query point are retrieved by performing a brute force linear search in the candidate set. In the MRPT algorithm, the candidate generation step consists of traversal of $T$ trees grown in the index construction phase.

The leaf to which a query point $\ve{q}$ belongs to is retrieved by first projecting $\ve{q}$ at the root node of the tree onto the same random vector as the data points, and then assigning it into the left or right branch depending on the value of the projection. If it is smaller than or equal to the cutpoint $s$ (median of the data points belonging to that node in the projected space) saved at that node, i.e. 
\[
\ve{q}^T \ve{r}_i \leq s, 
\] 
the query point is routed into the left child node, and otherwise into the right child node. This process is then repeated recursively until a leaf is met. 

The query point is routed down into a leaf in all the $T$ trees obtained in the index construction phase. The query process is thus far similar to the one described in \cite{Dasgupta2013}. The principal difference is the candidate set generation: in classic RP trees, the candidate set $S$ includes all the points that belong to the same leaf with the query point in at least one of the trees.  
%For a given query point $\ve{q}$ and its $i$:th nearest neighbor $x_{(i)}$, the events
%\[
%\{f_t(x_{(i)}; \ve{q})=0\}, \, t \in \{1, \dots, T\}
%\]
%that they end up into different leaves in $t$:th tree, are independent because the probability is over the random vectors, which are i.i.d. Hence, the probability that the algorithm fails to return $x_{(i)}$ is 
%\[
%P\{f(x_{(i)}; \ve{q})=0\} = \prod_{i=1}^n P\{f_t(x_{(i)}; \ve{q})=0\}. 
%\] 
A problem with this approach is that when a high number of trees are used in the tree traversal step, the size of the candidate set $|S|$ becomes excessively large. 

In the following, we show how the frequency information (i.e., how frequently a point falls into the same cell as query point) can be utilized to improve both query performance and accuracy.

\subsection{Voting search}

Assume that we have constructed $T$ RP trees of depth $\ell$. Each of them partitions $\R^d$ into $2^\ell$ cells (leaves) $L_1, \dots L_{2^\ell}$, all of which contain $\lceil \frac{n}{2^\ell} \rceil$ or $\lfloor \frac{n}{2^\ell} \rfloor$ data points. For $1\leq t\leq T$, let $f_t$ be an indicator function of data point $\ve{x} \in \ve{X}$ and the query $\ve{q}$, which returns $1$, if $\ve{x}$ and $\ve{q}$ reside in the same cell in tree $t$, and $0$ otherwise:
\[
f_t(\ve{x};\ve{q}) = \sum_{m=1}^{2^\ell} \mathds{1}\{\ve{x} \in L_m, \ve{q} \in L_m\}.
\]
Further, let $F$ be a count function of data point $\ve{x}$, which returns the number of trees in which $\ve{x}$ and $\ve{q}$ belong to the same leaf:
\[
F(\ve{x};\ve{q}) = \sum_{t=1}^T f_t(\ve{x};\ve{q}).
\]

We propose a simple but effective \emph{voting search} where we choose into the candidate set only data points residing in the same leaf as the query point in at least $v$ trees:
\[
S = \{\ve{x} \in \ve{X}: F(\ve{x};\ve{q}) \geq v \}.
\]
The vote threshold $v$ is a tuning parameter. A lower threshold value yields higher accuracy at the expense of increased query times.

This further pruning of the candidate set utilizes the intuitive notion that the closer the data point is to the query point, the more probably an RP tree divides the space so that they both belong to the same leaf. %We do not yet have any theoretical proof for this, but our experimental results show the superiority of this approach.
We emphasize that our voting scheme is not restricted to RP trees, but it can be used in combination with any form of space-partitioning algorithms, given that there is enough randomness involved in the process to render the partitions sufficiently independent.  

Pseudocode for the online stage of the MRPT algorithm is given in detail in Algorithms \ref{algo:3}--\ref{algo:4} below (the knn($\ve{q},k,S$) function is a regular $k$-NN search which returns $k$ nearest neighbors for the point $\ve{q}$ from the set $S$).

\begin{algorithm}
\begin{algorithmic}[1]
\Function{tree\_query}{$\ve{q}$, tree}
  \State $\ve{R}$ $\leftarrow$ tree.random\_matrix
  \State $\ve{p}$ $\leftarrow$ $\ve{q}^T \ve{R}$
  \State root $\leftarrow$ tree.root

  \For{level in $1, \dots, \ell$}
    \If{$\ve{p}$[level] $\leq$ root.split}
      \State root $\leftarrow$ root.left
    \Else
      \State root $\leftarrow$ root.right
    \EndIf
  \EndFor
  \State \textbf{return} data points in root
\EndFunction
\end{algorithmic}
\caption{Route a query point into a leaf in an RP tree.}
\label{algo:3}
\end{algorithm}

\begin{algorithm}
\begin{algorithmic}[1]
\Function{approximate\_knn}{$\ve{q}$, $k$, trees, $v$}
  \State $S \leftarrow \emptyset$
  \State let votes[$1 \dots \ve{X}\text{.nrows}$] be a new array
  \For{tree in trees}
    \For{point in \Call{tree\_query}{$\ve{q}$, tree}}
      \State votes[point] $\leftarrow$ votes[point] + 1
      \If {votes[point] $= v$}
        \State $S \leftarrow S \,\, \cup \{$point$\}$
      \EndIf
    \EndFor
    % \State $\ve{X} \leftarrow$ $\ve{X} \,\, \cup$ tree\_query($\ve{q}$, $t$, $\ell$)
  \EndFor
  \State \textbf{return} knn($\ve{q}, k, S$)
\EndFunction
\end{algorithmic}
\caption{Approximate $k$-NN search using multiple RP trees.}
\label{algo:4}
\end{algorithm}

% \begin{algorithm}
% \begin{algorithmic}[1]
% \Function{tree\_query}{$\ve{q}$, $t$, $\ell$}
%   \State tree $\leftarrow$ trees[$t$]  
%   \For{level in $1, \dots, \ell$}
%     \State projection $\leftarrow \ve{q}^T \ve{r}_{t,\mathrm{level}}$
%     \If{projection $\leq$ tree.split}
%       \State tree $\leftarrow$ tree.left
%     \Else
%       \State tree $\leftarrow$ tree.right
%     \EndIf
%   \EndFor
%   \State \textbf{return} data points in a leaf of tree
% \EndFunction
% \end{algorithmic}
% \caption{Route a query point into a leaf in an RP tree.}
% \label{algo:3}
% \end{algorithm}

% \begin{algorithm}
% \begin{algorithmic}[1]
% \Function{approximate\_knn}{$\ve{q}$, $k$, trees, $v$}
%   \State $S \leftarrow \emptyset$
%   \State let votes[$1 \dots \ve{X}\text{.nrows}$] be a new array
%   \For{$t$ in $1, \dots , T$}
%     \For{point in \Call{tree\_query}{$\ve{q}$, $t$, $\ell$}}
%       \State votes[point] $\leftarrow$ votes[point] + 1
%       \If {votes[point] $= v$}
%         \State $S \leftarrow S \,\, \cup \{$point$\}$
%       \EndIf
%     \EndFor
%     % \State $\ve{X} \leftarrow$ $\ve{X} \,\, \cup$ tree\_query($\ve{q}$, $t$, $\ell$)
%   \EndFor
%   \State \textbf{return} knn($\ve{q}, k, S$)
% \EndFunction
% \end{algorithmic}
% \caption{Approximate $k$-NN search using multiple RP trees.}
% \label{algo:4}
% \end{algorithm}

%{\bf TR: We should explain 1) what $knn(q,k,S)$ does; 2) how the union operation in Algo4 is implemented so that the list of points in the knn operation is optimized.}

% \begin{algorithm}
% \begin{algorithmic}[1]
% \Function{approximate\_knn}{$\ve{q}$, $k$, trees}
%   \For{$t$ in $1, \dots , T$}
%     \State $\ve{X} \leftarrow$ $\ve{X} \,\, \cup$ tree\_query($\ve{q}$, $t$, $\ell$)  
%   \EndFor
%   \State \textbf{return} knn($\ve{q}, k, \ve{X}$)
% \EndFunction
% \end{algorithmic}
% \caption{Approximate $k$-NN search using multiple RP trees.}
% \label{algo:4}
% \end{algorithm}

% \liang{if we have extra space, add some discussions on future work: e.g., look for theoretical foundation of voting scheme; can we use voting to completely get rid of original data set but only use mprt for query?}

\subsection{Time and space complexity: Query execution}

%\liang{I guess we can still call this section "Time and space complexity". Although voting does not help in reducing space complexity, we can still ad one sentence to clarify this.}
When the trees are grown into some predetermined depth $\ell$ using the median split, the expected running time of the tree traversal step is $\Theta(T\ell ad)$\footnote{Adding the leaf points  into the candidate set $S$ for each $T$ trees adds a term $\bigO{T \frac{n}{2^\ell}}$ into the time complexity of the tree traversal phase for defeatist search. Because checking the vote count is a constant time operation this term is the same for voting search.  However, in both cases these terms are dominated by the exact search in the search set, which is an $\bigO{T \frac{n}{2^\ell} d}$ operation, so they are not taken into account here.}, because at each level of each RP tree the query point is projected onto a $d$-dimensional random vector that has on average $ad$ non-zero components. When using random vectors sampled from the $d$-dimensional standard normal distribution, this is $\Theta(T\ell d)$, but when using value $a = \frac{1}{\sqrt{d}}$ for the sparsity parameter, it is only $\Theta(T\ell\sqrt{d})$. 

When the median split is used, each leaf has either $\lceil \frac{n}{2^\ell} \rceil$ or $\lfloor \frac{n}{2^\ell} \rfloor$ data points. Therefore, the size of the final search set satisfies $|S| \leq T \lceil \frac{n}{2^\ell} \rceil$ for both defeatist and voting search. Increasing the vote threshold $v$ decreases the size of the search set $S$ in practice, which explains the significantly reduced query times when using voting search (cf. experimental results in Section \ref{sec:results}). However, because the effect of $v$ on $|S|$ depends on the distribution of the data, a detailed analysis of the time complexity requires further assumptions on the data source and is not addressed here. Without any assumptions, both defeatist search and voting search have the same theoretical upper bound for $|S|$. 

Because of this upper bound for the size of the candidate set, the running time of the exact search in the candidate set is $\bigO{T \frac{n}{2^\ell} d}$ for both defeatist search and voting search. Thus, the total query time is $\bigO{Td(\ell + \frac{n}{2^\ell})}$ for classic RP trees, and $\bigO{Td(\frac{\ell}{\sqrt{d}} + \frac{n}{2^\ell})}$ for the MRPT algorithm. 

%which depends on two tuning parameters: number of trees ($T$), and depth of one tree ($\ell$). 

%This upper bound is necessarily not tight, because there are duplicates among the points returned by the different trees. 

% In the extreme case when all the points returned by $T$ trees are the same, $|S| = \lceil \frac{n}{2^\ell} \rceil$. In practice, however, we have noticed that in high-dimensional spaces, as long as sensible parameter values are used, most of the points returned by the different trees are unique, and thus $|S|$ is quite close to this upper bound.  {\bf TR: doesn't this suggest that voting is meaningless!?} 

%Additionally, using the voting scheme to prune the search set also loosens the upper bound, because only data points that get at least $v$ votes are chosen into the final search set. {\bf TR: We should really provide more precise comments on the effect of voting and parameter $v$ to the time complexity. Now our complexity analysis doesn't address the voting stage at all!}


 

%The obvious problem with this approach is that the nearest neighbors of the query point may very well end up into a different leaf than the query point itself. There are several strategies to mitigate this problem. In spill trees and virtual spill trees an overlap parameter is used, so that the points falling close to the split point in the projected space are not separated from the query point. However, our approach is completely different: we utilize multiple RP trees simultaneously, and combine their results. This strategy will be presented at Section \ref{sec:MRPT}.



%\section{MRPT algorithm}
%\label{sec:MRPT}
%
%In the multiple random projection trees algorithm (MRPT algorithm) $T$ random projection trees are grown from the same data set with a different set of random vectors. When answering a $k$-NN query, a query point is routed into a leaf at each of $T$ trees, and then a basic linear search is performed in the union of all the points in these leaves.
%
%If a depth of trees $\ell$ is kept constant, then growing more trees trivially increases the number of true nearest neighbors found because the size of the set $S$ where the final linear search is performed increases. When $T$ trees, each with depth $\ell$, are grown, for the search set size holds $|S| \leq T \lceil \frac{n}{2^\ell} \rceil$. 
%
%However, when increasing the number of trees used, it is possible to increase the depth, which decreases the leaf size, proportionally, so that an upper bound for the final search set size $|S|$ stays constant. This way, when a moderate amount of trees is used, the query time stays approximately constant.  
%
%There are some details in the algorithm that make it suitable for high-dimensional data. If the dimension of the data is very high, only the initial random seed is saved. When answering a $k$-NN query, all the same random vectors can be generated by just setting this initial seed.  When the projection vectors are not stored, but instead re-generated from the random seed on demand, only the split points are saved at each of the nodes. Thus the space required by the trees is independent of the dimension $d$, and depends only on the size $n$ of the data set. This keeps the size of the trees modest compared to the size of the data set in high-dimensional scenarios. 

%In case the data is only moderately high-dimensional, the random vectors can be stored to avoid the minor computational overhead of re-generating them when processing queries. We use the same random vector at all nodes on the same level of the tree. Thus, only $T\ell$ random vectors need to be stored, where $T$ is the number of trees, and $\ell$ is the depth of the trees. 


%\pagebreak


%\section{Time and space complexity}
%\label{sec:complexity}

%\subsection{Non-parallel implementation}


%The running time of the algorithm can be considered in two different ways. The first, and also the most relevant consideration here is how long it takes to search the $k$ (approximate) nearest neighbours for a new query point $q$. Approximate nearest neighbor search is very useful in applications where a fast query time is so critical that a linear search is not feasible. For example in online applications with massive databases, such as \cite{Wang2015}, a query time of even 5 seconds would be too much.

%The secondary consideration is how long it takes to grow the trees that are used for queries. Although the time to build the index is not so critical, it must be feasible, so that the indexing could be updated when new data is added; for example in the Kvasir application \cite{Wang2015} the indexing is rebuilt a few times a day and it takes a couple of hours to update it. 

%We first analyze the time complexity of a serial version of the algorithm: both the query time of one point, and the combined tree growing and query time for $n_q$ query points. Then we show that the same analysis holds for the parallel version of the algorithm.  

%\subsection{Query phase}

  

%An approximate k-NN query using the MRPT algorithm consists of two parts: queries in all $T$ trees, and a linear search in the union of all the points in the leaves a query point falls into. We call the set where a final linear search is performed the final search set, and it by $S$ and its size by $|S|$. If a stopping criterion for tree growing is set so that each leaf contains at most $n_0$ points, then  always $|S| \leq |S_{\max}| := Tn_0$, where $T$ is the number of trees used. Our experimental results show that when the dimensionality of the data is high, most of the points at the leaves are unique even when the number of trees $T$ is large, as long as $Tn_0$ stays relatively small compared to a number of data points $n$; so that an upper bound on the time complexity is obtained by using the bound $|S| \leq Tn_0$. 

%The average length of a query path in a tree is the average depth\footnote{For the median split this is the exact depth of the tree. For the other split methods trees are not perfectly balanced. For the bounded split the maximum depth is still logarithmic with respect to $n$. Time complexities presented here are worst case complexities for the median and the bounded split.} of the tree $\ell := \log \frac{n}{n_0}$, and at each level the query point is projected onto a $d$-dimensional random vector, so the time complexity of a query in a tree is $\bigO{d\log\frac{n}{n_0}}$ (all logarithms are to base 2, unless otherwise noted). When using $T$ trees, the total time complexity of the queries for one query point is $\bigO{Td\log\frac{n}{n_0}}$, if all the queries are ran serially. The time complexity of k-NN search in the final search set of size $|S| \leq Tn_0$ is $\bigO{Tdn_0}$. 



%Our experimental results show that increasing the number of trees used, while simultaneously increasing the depth of the trees proportionately, so that the upper bound for the search set size $T\lceil \frac{n}{2^\ell} \rceil$ stays constant, increases the quality of the nearest neighbor candidate set, and thus also the accuracy of the results. This means that as long as the final exact search dominates the running time, it is beneficial to increase the time spent on the tree traversal time by growing more proportionately deeper trees\footnote{When the amount of trees $T$ is doubled, depth of one tree $\ell$ has to be increased by one to keep the upper bound for the search set size $T\lceil \frac{n}{2^\ell} \rceil$ constant.}. Thus for the optimal parameter values holds
%\[
%\begin{split}
%T\ell d &= \bigO{T \frac{n}{2^\ell}} \\
%\ell 2^\ell &= \bigO{n}.
%\end{split}
%\]
%
%Because $f(x) = x 2^x$ is an increasing function, there is an unique solution for the equation 
%\begin{equation}
%\label{eq:optimal}
%\ell 2^\ell = n.
%\end{equation}
%Although asymptotic notation is used here, it is worth noting that the constant terms involved in the operations are really similar for both the tree traversal and the exact search when the dimension of the data $d$ is high, because both involve mostly matrix computations. Our experimental results show that the optimal value for $\ell$ indeed usually lies just below the solution for Equation \ref{eq:optimal}. Another remarkable fact is that the optimal value for the depth of the trees $\ell$ does not at least asymptotically depend on the number of trees $T$. This means that $T$ could used to set the desired accuracy level, while $\ell$ is set to its optimal value, which is independent of $T$.




%In order to avoid this overhead due to applying multiple trees becoming greater than the linear search in the final search set, we should make sure that time complexity of processing the queries in each of the trees, $\bigO{Td\log\frac{n}{n_0}}$, is not greater than the final search complexity, $\bigO{Tdn_0}$. This is achieved when $n_0$ is of the order $\log n$. The optimal values of $T$ and $n_0$ depend of course on the actual implementation of the algorithm.

%\subsection{Index construction phase}



%The secondary resource cost is the time it takes to grow the trees and the space required by them. When growing a tree, each data point is projected into $d$-dimensional random vector at each level of the tree; thus the time complexity of growing $T$ trees is $\bigO{Tnd \log \frac{n}{n_0}}$. As discussed above, it is only feasible to use a maximum leaf size $n_0$ that is of order $\bigO{\log n}$ or smaller, so when $\bigO{n_0} = \bigO{\log n}$ and the size of the search set $|S| \approx Tn_0$ is kept constant, the time complexity of growing the trees is:
%\[
%\begin{split}
%\bigO{Tnd \log \frac{n}{n_0}} &= \bigO{Tnd (n_0 - \log n_0) } \\
%							  &= \bigO{Tndn_0} \\
%							  &= \bigO{|S|nd}.
%\end{split}
%\]

%Besides online applications where the query times are critical, the algorithm can also be useful in applications where both the dimensionality $d$ and the size $n$ of the data are high, and the number of query points $n_q$ is several magnitudes larger than the number of data points ($n_q \gg n$). Then exact k-NN search for all the query points may be very slow, and approximate nearest neighbor search may save considerable amount of time, even when the index construction time is accounted for, because the trees have to be grown only once. With these assumptions the total time complexity of the algorithm is 
%\[
%\bigO{n|S|d + n_q|S|d} = \bigO{n_q|S|d}.						
%\]   
%Because the time complexity of exact k-NN search using basic linear search for $n_q$ query points is $\bigO{n_q nd}$, the multiple random projection trees algorithm is faster by a factor $\frac{|S|}{n} = \frac{T}{2^\ell} \leq \frac{T}{n}$. This means that the algorithm can be several magnitudes faster than linear search, even when the time to grow trees is accounted for. 

%The space complexity of the algorithm depends on the amount of memory required to store the trees and the corresponding leaf labels for each data point. Because for each node only the split point in the projected space is stored, the space required by one node is constant. So the space requirement of $T$ trees is 
%\[
%\bigO{T \left(\frac{2n}{n_0} - 1 \right) + n} = \bigO{Tn}.
%\] 

\begin{table}[tb!]
%\begin{footnotesize}
\begin{center}
\caption{Time and space complexity of the algorithm compared to classic RP trees.}
\label{table:complexity}
\begin{tabular}{lll}
\toprule
  & RP trees & MRPT $\left(a = \frac{1}{\sqrt{d}}\right)$\\
\midrule
%Candidate generation time & $\Theta(T\ell d)$ & $\Theta(T\ell\sqrt{d})$ \\ 
%Final search time& $\bigO{T \frac{n}{2^\ell} d}$ & $\bigO{T \frac{n}{2^\ell} d}$ \\
Query time & $\bigO{Td(\ell + \frac{n}{2^\ell})}$ & $\bigO{Td(\frac{\ell}{\sqrt{d}} + \frac{n}{2^\ell})}$ \\
Index construction time & $\Theta(T\ell nd)$ & $\Theta(T\ell n\sqrt{d})$ \\
Index memory & $\bigO{Tdn}$ & $\bigO{T(\sqrt{d}\log n + n)}$ \\
\bottomrule 
\end{tabular}
\end{center}
%\end{footnotesize}
\end{table}


%\section{Parallel implementation}
%
%\VH{Scrap or rewrite this section.}
%\liang{OK, I will comment after you finish rewriting. But note: when you talk about parallelism, discuss in both indexing and query phases. Especially in indexing phase, note the two levels of parallelism: 1) parallelism between independent trees which can be trivially achieved; 2) parallelism within a tree, which is not trivially achieved but your smartly solved in Section III.C; 3) then talk about traffic overhead and other practical challenges ...}
%
%Both the growing of the trees and processing of queries in the trees are independent operations for each of the $T$ trees, so the most obvious and efficient way to parallelize the algorithm is to divide the growing of trees between all the nodes. 
%
%The algorithm is easy to implement in parallel with very little communication overhead. If the trees are grown and stored in separate nodes, then the queries in the trees and the computing of distances in the final search set can be done completely in parallel, and only a final step, where the results are combined, has to be done serially. Because the distances from the query point to all points in the final search set are already computed in parallel, this final step is just a linear scan of the final search set, and so its time complexity is 
%\[
%\bigO{Tn_0} = \bigO{|S|},
%\]
%which is negligible.
%
%Because the trees grown in separate nodes can be utilized in the same nodes where they are constructed, the only communication overhead in the parallel algorithm comes from the distribution of the data set and the query points to the nodes, and then collecting the final results. Distribution of the data set into nodes has to be done only in the tree growing phase, so it does not have any effect on the query times. Because the final results communicated to the master node are just distances from the query point to all points in the search set, their space complexity is $\bigO{S}$, and the $d$-dimensional query point is transferred into each of the worker nodes, so the amount of data that has to transferred for each query point is only of order $\bigO{S+pd}$, where $p$ is the number of nodes. 
%
%If $p$ worker nodes are available for parallel computing, the total number of trees $T$ is chosen so that $T = mp$ where
%$m$ is an integer, because otherwise the workload would be unevenly distributed between the nodes. Specifically always $T \geq p$. 
%
%Now for each node $i = 1, \dots , p$ the analysis presented above holds. The total number of points in the leaves for which the distances to the query point have to be computed is $\frac{T}{p}n_0 = mn_0 = Sp^{-1}$, and the time complexity of the distance calculations in the final search set is 
%\[
%\bigO{mn_0d} = \bigO{Sp^{-1}d}
%\]
%for each of the nodes. With the same notation a query in all $m$ trees in one node takes $\bigO{m d \log \frac{n}{n_0}}$ time, and if the maximum leaf size $n_0$ is again set so that $\bigO{n_0} = \bigO{\log n}$, this becomes
%\[
%\bigO{m d \log \frac{n}{n_0}} = \bigO{mn_0 d} = \bigO{Sp^{-1}d}.
%\] 
%
%With similar calculations it can be shown that the time complexity of growing all $m$ trees is 
%\[
%\bigO{Sp^{-1}nd}
%\] 
%for each node, when $n_0$ is chosen so that $\bigO{n_0} = \bigO{\log n}$.
%
%So also for the parallel algorithm the choice of the number of trees $T$ and maximum leaf size $n_0$ so that the magnitude of the $n_0$ is $\bigO{\log n}$ gives optimal performance for a given $S \approx Tn_0$. An additional consideration for the parallel algorithm is that the number of trees should be a multiple of the number of worker nodes utilized.
% 
%One possible issue arises when the data set is so big that the whole data cannot be transferred into each of the worker nodes. In this case the algorithm can be parallelized by dividing the data set evenly between the nodes, executing MRPT algorithm at each node, and in the end searching the nearest ones from the potential nearest neighbors returned by the nodes. However this approach is not as efficient as distributing the trees between the nodes as described above: if the same values of $T$ and $n_0$ are used in both approaches, query times are almost identical because $|S_{\max}|$ is identical, and the trees used in queries are not significantly shallower in the second approach, unless the number of nodes $p$ is really large. This is because the depth of each tree is $\log \frac{n}{pn_0}$ in the second approach where the data set is divided between the nodes, and $\log \frac{n}{n_0}$ in the first approach where all trees are grown using the whole data set. So query times are practically the same, but in the second approach the partial data sets of size $\frac{n}{p}$ are used to grow each of the trees, whereas in the first approach the whole data set is used for each of the trees. Clearly the first approach is more efficient (if only query times and accuracy are considered, of course trees are grown faster in the second approach) but for really large data sets the second approach has to be taken. 
%
%%\subsection{Parallel MPRT}
%%
%%\liang{[Chaps, please have a closer look into this algorithm. If you notice any bugs, let's fix it. There are another two things need your attention: (1) I only wrote down some crude explanations to help you understand the algorithm, please polish the text to make the style consistent with your previous text; (2) feel free to tune the notations and variable to make it read better. In case we have space, please expand the "if ... then ..." logic. ]} 
%% 
%% \begin{algorithm}
%%\begin{algorithmic}[1]
%%\Function{build\_trees\_parallel}{ }
%%  \While{ $Q \neq \varnothing$ }
%%    \State $S_1, s, j \leftarrow Q.$dequeue( )
%%    \State Generate $\ve{r}$ from $N_d(\ve{0},\ve{I})$ using $s$.
%%	\State $S_2 \leftarrow S_1.$map \{ lambda $i,\ve{x}$: $(i, \ve{x}^T \ve{r})$ \}
%%	\State $T[j] \leftarrow$ split ($S_2$)
%%	\State $S_3 \leftarrow S_2.$map \{
%%	\State \quad lambda $i, x$:
%%	\State \quad \quad \textbf{if} $x < T[j]$
%%		\textbf{then} $(0, i)$ \textbf{else} $(1, i)$ 
%%	\State \}
%%	\State $L, R$ $\leftarrow S_3$.reduce\_by\_key ( )
%%
%%	\State $m \leftarrow 2j+1$; $n \leftarrow 2j+2$
%%	\If{$|L| > n_0$} 
%%		\State $Q$.enqueue$(L,s+a,m)$
%%	\Else 
%%		\State $T[m] \leftarrow L$
%%	\EndIf	
%%	\If {$|R| > n_0$} 
%%		\State $Q$.enqueue$(R,s+b,n)$
%%	\Else 
%%		\State $T[n] \leftarrow R$
%%	\EndIf
%%  \EndWhile
%%  \State \textbf{return} $T$
%%\EndFunction
%%\end{algorithmic}
%%\end{algorithm}
%%
%%\liang{[argue that the recursive form is slow and may cause stack overflow especially when the tree is highly imbalanced ... ]}
%%The algorithm above provides a parallel version of the previous recursive algorithm. Before we explain the algorithmic details, we introduce the notations and variables used in the algorithm. 
%%
%%$T$ is an array representing a RP tree. Each element in $T$ represents a node in the tree. Note that our choice of using an array to present the algorithm is for the purpose of simplicity. An array may become inefficient if the tree is highly imbalanced, since there will be many unused elements. However, we can easily and efficiently re-construct the same tree using more compact data structure from $T$. We can also directly use other data structures which may complicate the presentation of the algorithm, however, the logic of the parallel algorithm remains the same.
%%
%%$Q$ is a thread-safe \textit{First-In, First-Out} (FIFO) queue wherein each element is a tuple of $(S, s, j)$. $S$ is a set which contains the indices of vectors in the original data space. The algorithm can use the indices in $S$ to access the original vectors in the data set. $s$ is the random seed which is used for generating a random vector for projection. $j$ is the position index of a node in $T$. Given a binary tree and its array representation, we can easily calculate the indices of its left- and right-child as $2j+1$ and $2j+2$ respectively.
%%
%%Initially, $Q$ has only one element $(X, s, 0)$. $X$ contains all the indices in the data set, and $0$ indicates that we are at the root of the RP tree (i.e., node index is zero). We start the parallel tree building by spawning multiple worker processes each of which executes the exact algorithm above. A worker process first dequeues an element by calling $Q.$dequeue(), then extracts the three variables from the tuple (line 3). Then a random vector $r$ is generated using the random seed $s$.
%%
%%Line 5 utilises a map function to project the vectors in $S_1$ to $r$. For each vector in $S_1$, the actual projection operation is done by calling an anonymous function lambda which takes $(i, x)$ as its input and outputs $(i, x^T r)$ where $i$ denotes the vector index and $x$ denotes the corresponding vector value.
%%
%%Line 6 then calculates the the split point from the project space $S_2$, then saves it to the corresponding tree node in $T[j]$.
%%
%%Then the algorithm uses another map function (line 7 -- 10) to split the plane, or more precisely, to check which branch the point falls in. Similar to the previous map function in line 5, the actual operation is done by creating another anonymous function which outputs either $(0,i)$ or $(1,i)$ where $0$ represents left child and $1$ represents right child.
%%
%%Reduce by key function in line 11 creates two sets $L$ and $R$ which contain the indices of the left-space points and the right-space points respectively. Line 13 -- 20 further check $L$ and $R$, if the size (i.e., cardinality of the set) is smaller than the threshold $n_0$, then the algorithm saves the cluster to the leaf node. Otherwise, the algorithm makes a new tuple and queues it into $Q$. The algorithm stops when $Q$ becomes empty. \liang{[explain a and b, chaps, or simply refer to the recursive form.]}
%%
%%\liang{[Note this, expand a bit if you have time.]}
%%There are two levels of parallelism in the algorithm. First level is achieved by using a queue and multiple worker processes. As we are going down from the root to the leaves, the parallelism will increase as the number of leaves increases. As the problem space is being continuously split into half, the "task size" also becomes smaller and smaller accordingly.
%%
%%As we already mentioned, the initial runs in the recursive tree building algorithm are slow due to the large problem space even with first level parallelism. In other words, there are many points given a large space. This issue is solved by the second level of parallelism which utilises the \textit{MapReduce} programming paradigm. The projection operation of each vector is independent hence it fits very well into data-parallel computing paradigm. The algorithm can be implemented on multi-core computers in a straightforward way or easily ported onto Apache Spark.
%%
%%
%%\subsection{On Very Big Data Sets}
%%
%%\liang{[Maybe add some discussions on the super big data sets as I mentioned in the previous emails. E.g., when the data set has to be saved across the whole cluster, we can probably take each chunk as input and apply MPRT on top of that. When we query, the query is distributed over the cluster then the results are aggregated .... blah blah blah ... chaps, use our imagination here. :) ]}
%%

 
\section{Experimental results}
\label{sec:results}

We assess the efficiency of our algorithm by comparing its against the state-of-the-art approximate nearest neighbor search algorithms. To make the comparison as practically relevant as possible, we chose widely used methods that are available in optimized libraries. % and are known to perform well in benchmarks. {\bf TR: cite for example the Spotify benchmarks and any others that show that these are the best competitors.} 

All of the compared libraries, including ours, are implemented in C++ and compiled with similar optimizations. We make our comparison and the implementation of our algorithm available as open-source\footnote{\url{https://github.com/ejaasaari/mrpt-comparison}}.

All of the experiments were performed on a single computer with two Intel Xeon E5540 2.53GHz CPUs and 32GB of RAM. No parallelization beyond that achieved by the use of linear algebra libraries, such as caching and vectorization of matrix operations, was used in any of the experiments.

\subsection{Data sets}
We made the comparisons using several real-world data sets (see Table \ref{table:datasets}) across a wide range of sample size and dimensionality. The data sets used are typical in applications in which fast nearest neighbor search is required. 
For all data sets, we use a disjoint random sample of 100 test points to evaluate the query time and recall.

The BIGANN and GIST data sets\footnote{\url{http://corpus-texmex.irisa.fr/}} consist of local 128-dimensional SIFT descriptors and global 960-dimensional color GIST descriptors extracted from images \cite{Jegou2011A, Jegou2011B}, respectively. As the SIFT data set, we use a random sample of $n = 2500000$ data points from the BIGANN data set.

In addition to feature descriptor data sets, we consider three image data sets for which no feature extraction has been performed: 1) the MNIST data set\footnote{\url{http://yann.lecun.com/exdb/mnist/}} consists of $28 \times 28$ pixel grayscale images each of which is represented as a vector of dimension $d=784$; 2) the Trevi data set\footnote{\url{http://phototour.cs.washington.edu/patches/default.htm}} consists of $64 \times 64$ pixel grayscale image patches represented as a vector of dimension $d=4096$ \cite{Winder2007}; and the STL-10 data set\footnote{\url{http://cs.stanford.edu/~acoates/stl10}}, which is an image recognition data set containing $96 \times 96$ pixel images from 10 different classes of objects \cite{Coates2011}.

%The audio data set is a sample of features extracted from English sentences spoken by different speakers. It was used as a test data set in \cite{Dong2008}, and further information about the preprocessing steps of the data set can be found there. The data set is available at \url{http://www.cs.princeton.edu/cass/audio.tar.gz}. 

%For word vectors, distance between vectors can be used to measure semantic similarity of the corresponding words. The GloVe data sets\footnote{\url{http://nlp.stanford.edu/projects/glove/}} consist of word vectors trained from different corpora \cite{Pennington2014}. We used the data set pre-trained from tweets consisting of 200-dimensional vectors.

We also use a news data set that contains web pages from different news feeds converted into a term frequency-inverse document frequency (TF-IDF) representation. Dimensionality of the data is reduced to $d=1000$ by applying latent semantic analysis (LSA) to the TF-IDF data. More elaborate description of the preprocessing of the news data set is found in \cite{Wang2016}. 

Finally, we include results for a synthetic data set comprised of $n = 50000$ random vectors sampled from the 4096-dimensional standard normal distribution and normalized to unit length.

\begin{table}[tb!]
%\begin{footnotesize}
\begin{center}
\caption{Data sets used in the experiments}
\label{table:datasets}
\begin{tabular}{l l l l}
\toprule
Data set & $n$ & $d$ & type \\
\midrule
%GloVe & 1193514 & 200 \\
GIST & 1000000 & 960 & image descriptors \\
SIFT & 2500000 & 128 & image descriptors\\
MNIST & 60000 & 784 & image \\
Trevi & 101120 & 4096 & image\\
STL-10 & 100000 & 9216 & image\\
News & 262144 & 1000 & text (TF-IDF + LSA)\\
Random & 50000 & 4096 & synthetic\\
\bottomrule 
\end{tabular}
\end{center}
%\end{footnotesize}
\end{table}

% We provide the results for sample sizes that are exact powers of two to better illustrate the effect of different leaf sizes (caused by uneven divisions) on the results when using the median split. All the values of the precision and the query time represented are averages over the 100 query points.

%\subsection{Effect of the number of trees}
%\label{sub:comp_number}
%
%\begin{figure*}[!ht]
%\includegraphics[width=0.5\textwidth]{images/news_262144_to_16384.pdf} \includegraphics[width=0.5\textwidth]{images/news_262144_to_16384_times.pdf}
%\caption{The average proportion of true 10 nearest neighbors found as a function of the number of trees (\emph{left})
%and query time (\emph{right}). Different curves correspond to different values of the upper bound for the total search set size, $|S_{\max}|$, while different points on a given curve correspond to different number of trees. (news data set, $n = 262144$, $d=1000$).}
%\label{fig:news8}
%\label{fig:news8_times}
%\end{figure*}
%
%\begin{table*}[ht!]
%\begin{scriptsize}
%
%\begin{center}
%\caption{Recall (R) and query time per point (t) for the news data set ($n=262144$, $d=1000$) with upper bounds for the search set size $|S_{\max}|$ = 2048, 4096, 8192, and 16384, respectively. (Exhaustive search takes 154.9 ms.)}
%\label{table:news1024}
%\begin{tabular}{ccccccccc}
%\toprule
%$k$ & $(T, n_0)$ & (1,2048) & (2,1024) & (4,512) & (8,256) & (16,128) & (32,64) & (64,32)  \\
%\midrule
%& t (ms)  &  1.5 &  1.5 &  1.8 &  2.4 &  3.3 &  5.4 &  9.4 \\ 
%10 & R (\%)  &  14.7 &  16.2 &  32.1 &  45.2 &  48.7 &  61.6 &  70.8 \\ 
%50 & R (\%) &  11.2 &  11.7 &  22.6 &  31.7 &  39.0 &  47.7 &  57.2 \\ 
%100 & R (\%) &  9.7 &  9.9 &  19.4 &  26.2 &  32.8 &  41.1 &  49.4 \\  
%\bottomrule
%\end{tabular}
%
%\vspace{15pt}
%
%\begin{tabular}{cccccccccc}
%\toprule
%$k$ & $(T, n_0)$ & (1,4096) & (2,2048) & (4,1024) & (8,512) & (16,256) & (32,128) & (64,64) & (128,32) \\
%\midrule
%& t (ms) &  2.8 &  2.9 &  3.2 &  3.7 &  4.7 &  6.7 &  10.8 &  18.8 \\
%10 & R (\%) &  14.2 &  26.1 &  43.3 &  50.1 &  63.4 &  71.0 &  81.4 &  89.5 \\ 
%50 & R (\%) &  10.4 &  19.7 &  31.5 &  40.3 &  51.2 &  59.4 &  69.8 &  79.4 \\  
%100 & R (\%) &  10.0 &  16.8 &  27.3 &  34.2 &  45.5 &  52.7 &  62.1 &  71.8 \\ 
%\bottomrule 
%\end{tabular}
%
%\vspace{15pt}
%
%\begin{tabular}{ccccccccccc}
%\toprule
%$k$ & $(T, n_0)$ & (1,8192) & (2,4096) & (4,2048) & (8,1024) & (16,512) & (32,256) & (64,128) & (128,64) & (256, 32) \\
%\midrule
%& t (ms)  &  5.6 &  5.7 &  6.0 &  6.5 &  7.6 &  9.6 &  13.6 &  21.6 &  37.7 \\ 
%10 & R (\%)  &  21.2 &  35.4 &  46.1 &  62.6 &  76.3 &  84.3 &  90.1 &  93.7 &  96.5 \\
%50 & R (\%)  &  17.6 &  29.2 &  38.9 &  52.4 &  64.9 &  73.1 &  82.4 &  88.3 &  93.4 \\ 
%100 & R (\%)  &  15.4 &  25.1 &  34.5 &  46.8 &  57.0 &  66.7 &  76.3 &  82.8 &  89.7 \\ 
%\bottomrule
%\end{tabular}
%
%\vspace{15pt}
%
%\begin{tabular}{cccccccccccc}
%\toprule
%$k$ & $(T, n_0)$ & (1,16384) & (2,8192) & (4,4096) & (8,2048) & (16,1024) & (32,512) & (64,256) & (128,128) & (256, 64) & (512, 32) \\
%\midrule
%& t (ms)  &  11.1 &  11.2 &  11.5 &  12.1 &  13.1 &  15.2 &  19.3 &  27.4 &  43.8 &  76.4 \\
%10 & R (\%) &  26.6 &  42.9 &  57.7 &  76.2 &  84.7 &  92.0 &  95.1 &  96.7 &  98.2 &  99.7 \\ 
%50 & R (\%)  &  24.3 &  39.0 &  49.4 &  65.3 &  73.5 &  85.9 &  91.2 &  95.6 &  97.4 &  98.4 \\  
%100 & R (\%)  &  21.7 &  35.8 &  45.9 &  59.0 &  68.2 &  81.0 &  87.5 &  93.6 &  96.7 &  98.3 \\ 
%\bottomrule
%\end{tabular}
%
%
%\end{center}
%\end{scriptsize}
%\end{table*}
%
%
%First, we study the effect of the number of trees on the accuracy when the total search set size $|S| = Tn_0$ is kept constant. The accuracy of the algorithm is measured by an average of recall, which is the proportion of true $k$ nearest neighbors of the query returned by the algorithm:
%\[
%\mathrm{Recall} = \frac{\textrm{\# of true }k \textrm{-NN returned}}{k}.
%\]
%In the left panel of Fig.~\ref{fig:news8}, an upper bound for the search set size, $|S_{\max}|$, is kept constant for each of the curves. For example the first point of the red curve represents the situation where one tree with a maximum leaf size $n_0=2048$ is grown. At the second point of the red curve two trees with maximum leaf size $n_0 = 1024$ are grown, at the third point four trees with $n_0=512$ are grown and so on. It can clearly be seen that increasing the number of trees really increases the accuracy significantly, even when keeping the size of the final search set constant.
%
%According to the analysis of Section \ref{sec:complexity}, increasing the number of trees while keeping the search set size constant should not significantly increase the query times, as long as the final linear search dominates the total running time of the algorithm. The right panel of Fig.~\ref{fig:news8_times}, where recall is given as a function of the actual query times for the same parameter values as in the left panel of the same figure, it can be seen that this is indeed the case. At the beginning each of the curves rises almost vertically, which means that the time taken by the traversal of trees is trivial compared to a time taken by the final search. However, there is a quite sharp point of diminishing returns as the traversal of trees starts to dominate the running time. If more accuracy is required after this point, it is more effective to increase the search set size than number of trees. 
%
%The optimal  parameter combinations can be read directly from the figure. The curve that is highest at a given point gives the highest recall for the given query time, or the shortest query time for the given recall level. If for example the best average recall in 0.5 milliseconds is required, the best option from the ones tested is to grow $T=16$ trees with maximum leaf size $n_0=32$ (the fifth point of the blue curve). This parameter combination yields the mean recall of roughly 70\%. These values can also be read from Table~\ref{table:news1024}.
%
%Another main point of interest is a comparison to exhaustive brute force search. For the mnist data set a basic exhaustive search takes 15 milliseconds for one query point. With the MRPT algorithm, for example a $k=10$ nearest neighbor query with 85~\% mean accuracy can be achieved in 0.74 ms, and 95~\% mean accuracy in 1.3 ms. For the news data set and $k=10$, exhaustive search takes 155 milliseconds, whereas with the MRPT algorithm 85~\% mean accuracy is obtained in 9.6 ms, and over 95~\% accuracy is obtained in 19.3 ms.  
%
%Although the time taken by an exhaustive search for the largest data set tested is only 155 ms, these results are  suggestive of a more general trend. For 10-nn queries, mean accuracy of 90~\% is achieved with a 20 fold speed-up compared to the exact exhaustive search for the mnist data set. In the news data set, the same accuracy is obtained with a 10 fold speed-up. 
%The main applications of the algorithm will deal with ultra-high dimensional data sets, and because the algorithm is scalable and easily parallelizable, it is an attractive approach for real-time applications.

\subsection{Compared algorithms}

% We consider only recall levels over 0.8 since these are the most relevant for the applications of approximate nearest neighbor search.

We compare the MRPT algorithm to  representatives from all three major groups of approximate nearest neighbor search algorithms: tree-, graph- and hashing-based methods.

To test the efficiency of our proposed improvements to classic RP trees, we include our own implementation\footnote{The version of RP trees used in the comparisons is already a slightly improved version: we used RP trees with the same random vector at every level and median split (cf. section \ref{sec:growing}). If we had used a different random vectors at each node of the tree we would have ran out of memory on the larger data sets. However, we emphasize that using the same random vector at every level actually slightly decreases the query times while keeping the accuracy intact; hence, the real performance gap between RP trees and the MRPT algorithm is actually even wider than shown in the results.} of classic RP trees suggested in \cite{Dasgupta2013}, and an MRPT-algorithm with vote threshold $v=1$, which amounts to using sparse RP trees without the voting scheme.

Of the tree-based methods we include randomized $k$-d trees \cite{Silpa2008} and hierarchical $k$-means trees \cite{Muja2009}. Authors of \cite{Muja2014} suggest that these two algorithms give the best performance on data sets with a wide range of dimensionality. Both of the algorithms are implemented in FLANN\footnote{\url{http://www.cs.ubc.ca/research/flann/}}.

%For k-d trees, we use the grid method to search for the optimal combinations of points to check ($c$) and number of trees to use ($T$). For hierarchical k-means trees, we similarly search for the optimal $c$ and branching factor $b$. We set the number of k-means iterations as 15.

We also use as a benchmark the classic $k$-d tree algorithm~\cite{Bentley1975} modified for approximate nearest neighbor search as described in~\cite{Arya1998}; it is implemented in the ANN library\footnote{\url{http://www.cs.umd.edu/~mount/ANN/}}.

% We use standard search with different values for the error bound $\epsilon$, and the sliding midpoint split rule as suggested by the authors.

Among the graph-based methods, the NN-descent algorithm suggested in~\cite{Dong2011} is generally thought to be a leading method. It first constructs a $k$-NN graph using local search in the offline phase, and then utilizes this graph for fast $k$-NN queries in the online phase. The algorithm is implemented in the KGraph library\footnote{\url{http://www.kgraph.org}}.

%We search for the optimal parameters by varying the cost parameter $P$ and the number of neighbors for each node in the constructed graph ($K$).

As a representative of hashing-based solutions we include a multi-probe version of cross-polytope LSH suggested in~\cite{Andoni2015} as implemented in the FALCONN library\footnote{\url{http://www.falconn-lib.org}}. All these algorithms and implementations were selected because they represent the state-of-the-art in their respective categories.


% To tune the algorithm, we search for the optimal number of hash tables and hash functions to use.

For the most important tuning parameters (as stated by the authors), we use grid search on the appropriate ranges to find the optimal parameter combinations in the experiment, or the default values whenever provided by the authors. 
The performance of the MRPT algorithm is controlled by the number of trees, the depth of the trees, the sparsity parameter $a$, and the vote threshold. We fix the sparsity parameter as $a = 1/\sqrt{d}$, and optimize the other parameters using the same grid search procedure as for the other methods. For ANN, we vary the error bound $\epsilon$, and for FLANN, we control the accuracy primarily by tuning the maximum number of leafs to search, but also vary the number of trees for $k$-d trees and the branching factor for hierarchical $k$-means trees. KGraph is primarily controlled by the search cost parameter, but we also build the graph with varying number of nearest neighbors for each node. For FALCONN, we search for the optimal number of hash tables and hash functions. The supplementary material on GitHub provides the exact parameter ranges of grid search as well as interactive graphs that show the optimal values for each obtained result.

% For the most important tuning parameters (as stated by authors), we use grid search on the appropriate ranges to find the optimal parameter combinations. For the other parameters, we use values recommended by the authors.  The exact parameter values used can be found on our comparison code available on github. {\bf TR: do we provide any sensitivity analyses on github in the form of results for other than optimal parameter settings? If yes, this should be mentioned here.}

% Multi-probe LSH implemented in the LSHKIT C++ library,\footnote{\url{http://lshkit.sourceforge.net/}} is an efficient version of LSH~\cite{Dong2008}. We use automatic tuning system provided with the library to tune the number of hash functions $M$ and width of the projection $W$. We set the number of hash tables to $L=10$, and the number of points to check per hash table to $T=100$ which seemed to provide consistently efficient performance with different data sets and values of $k$.

\subsection{Results}

In the experiments we used three different values of $k$ that cover the range used in typical applications: $k=1, 10$ and $100$. However, due to space restrictions, we present results for multiple values of $k$ only for MNIST, which is representative of the general trend. For the other six data sets we only provide results for $k=10$. (The results for $k=1,100$ are available in the supplementary material on GitHub).

Figure~\ref{fig:comparison} shows results for $k = 10$ on six data sets. The times required to reach a given recall level are shown in Table \ref{table:results}. The MRPT algorithm is significantly faster than other tree-based methods and LSH on all of the data sets, expect SIFT. It is worth noting that SIFT has a much lower dimension ($d=128$) than the other data sets. These results suggest that the proposed algorithm is more suitable for high-dimensional data sets.

For the MNIST data set, results for all three values of $k$ (and additionally for $k=25$) are shown in Figure~\ref{fig:diff_k}.   For the other data sets the results for $k=1 $ and $k=100$ follow a similar pattern as the results for the MNIST data set: the choice of $k$ has no significant effect to the relative performance of the other methods, but the relative performance of KGraph seems to degrade significantly for small values of $k$. 

KGraph is the fastest method on low recall levels, but its performance degrades rapidly on high recall levels (over 90\%), thus making the MRPT algorithm the fastest method on very high recall levels (over 95\%) on all the data sets except SIFT. Our hypothesis is that the flexibility obtained on the one hand by using completely random split directions, and on the other hand by using a high number of trees, which is enabled by sparsity and voting, allows the MRPT algorithm to obtain high recall levels for high-dimensional data sets very efficiently.  

The performance of the MRPT algorithm is also superior to classic RP trees on all of the data sets. In addition, from our two main contributions, voting seems the be more important than sparsity with respect to query times: sparse RP trees are only slightly faster than dense RP trees, the gap being somewhat larger on the higher-dimensional data sets, but the voting search provides a marked improvement on all data sets. This shows that the voting search, especially when combined with sparsity, is an efficient way to reduce query time without compromising accuracy.

The numerical values in Table \ref{table:results} indicate that for recall levels $\geq$ 90\%, the MRPT method is the fastest in 14 out of 21 instances, while the KGraph method is fastest in 7 out 21 cases. Compared to brute force search, MRPT is roughly 25--100 times faster on all six real-world data sets  at 90\% recall level, and roughly 10--40 times faster even at the highest 99\% recall level.

%{\bf TR: Mention Table 1 (and add a caption to it).}

%We observe that the effect of sparsity is significant mainly on the data sets with the highest dimensionalities, Trevi and STL-4096. As matrix multiplication with dense matrices can make effective use of caches and vectorization much better than with sparse matrices, the dimensionality needs to be high before the speed-up gained by using sparse vectors is significant. On the other hand, voting can be seen to have a significant impact on all data sets. Thus, voting proves to be an efficient way to reduce the size of the search set without sacrificing accuracy.

%The MRPT algorithm is significantly faster than the other tree-based methods and multi-probe LSH on almost all data sets. Only on the data set with the lowest dimensionality, SIFT, hierarhical k-means trees seem to be more efficient on lower recall rates. Our method is also faster than KGraph on high recall rates, especially on high-dimensional data.

% It can be seen that the MRPT algorithm is significantly faster for both data sets and values of $k$ for the given recall level than multi-probe LSH. For the mnist data set MRPT is significantly faster than k-d tree (ANN). For the news data set the margin is narrower. Randomized k-d trees (FLANN) seem to be more efficient with lower levels of recall, but especially with the larger news data set their performance seems to drop rather quickly for high levels of recall. It seems that the relative performance of MRPT algorithm is highest when the number data points $n$, desired recall level, and the number of neighbors searched $k$ are all high: in particular, the steep drop of performance that is seen for randomized k-d trees does not seem to occur with MRPT until very near the 100~\% recall level. Our hypothesis is that although computing random projections slows the algorithm for  low recall levels compared to k-d trees which require no projections, the benefits of the flexibility allowed by arbitrary projections guarantee high efficiency even at  very high recall levels.


%\begin{figure*}[!ht]
%\includegraphics[width=\textwidth]{images/compare_new_version.pdf}
%\caption{Comparison of query times with mnist ($n = 32678, d = 784$) and news ($n = 262144, d = 1000)$ data sets. The methods are: multi-probe LSH (LSHKIT), k-d tree (ANN), randomized k-d tree (FLANN), and our method (MRPT).}
%\label{fig:comparison}
%\end{figure*}


% An example of a floating figure using the graphicx package.
% Note that \label must occur AFTER (or within) \caption.
% For figures, \caption should occur after the \includegraphics.
% Note that IEEEtran v1.7 and later has special internal code that
% is designed to preserve the operation of \label within \caption
% even when the captionsoff option is in effect. However, because
% of issues like this, it may be the safest practice to put all your
% \label just after \caption rather than within \caption{}.
%
% Reminder: the "draftcls" or "draftclsnofoot", not "draft", class
% option should be used if it is desired that the figures are to be
% displayed while in draft mode.
%
%\begin{figure}[!t]
%\centering
%\includegraphics[width=2.5in]{myfigure}
% where an .eps filename suffix will be assumed under latex, 
% and a .pdf suffix will be assumed for pdflatex; or what has been declared
% via \DeclareGraphicsExtensions.
%\caption{Simulation Results}
%\label{fig_sim}
%\end{figure}

% Note that IEEE typically puts floats only at the top, even when this
% results in a large percentage of a column being occupied by floats.


% An example of a double column floating figure using two subfigures.
% (The subfig.sty package must be loaded for this to work.)
% The subfigure \label commands are set within each subfloat command, the
% \label for the overall figure must come after \caption.
% \hfil must be used as a separator to get equal spacing.
% The subfigure.sty package works much the same way, except \subfigure is
% used instead of \subfloat.
%
%\begin{figure*}[!t]
%\centerline{\subfloat[Case I]\includegraphics[width=2.5in]{subfigcase1}%
%\label{fig_first_case}}
%\hfil
%\subfloat[Case II]{\includegraphics[width=2.5in]{subfigcase2}%
%\label{fig_second_case}}}
%\caption{Simulation results}
%\label{fig_sim}
%\end{figure*}
%
% Note that often IEEE papers with subfigures do not employ subfigure
% captions (using the optional argument to \subfloat), but instead will
% reference/describe all of them (a), (b), etc., within the main caption.


% An example of a floating table. Note that, for IEEE style tables, the 
% \caption command should come BEFORE the table. Table text will default to
% \footnotesize as IEEE normally uses this smaller font for tables.
% The \label must come after \caption as always.
%
%\begin{table}[!t]
%% increase table row spacing, adjust to taste
%\renewcommand{\arraystretch}{1.3}
% if using array.sty, it might be a good idea to tweak the value of
% \extrarowheight as needed to properly center the text within the cells
%\caption{An Example of a Table}
%\label{table_example}
%\centering
%% Some packages, such as MDW tools, offer better commands for making tables
%% than the plain LaTeX2e tabular which is used here.
%\begin{tabular}{|c||c|}
%\hline
%One & Two\\
%\hline
%Three & Four\\
%\hline
%\end{tabular}
%\end{table}


% Note that IEEE does not put floats in the very first column - or typically
% anywhere on the first page for that matter. Also, in-text middle ("here")
% positioning is not used. Most IEEE journals/conferences use top floats
% exclusively. Note that, LaTeX2e, unlike IEEE journals/conferences, places
% footnotes above bottom floats. This can be corrected via the \fnbelowfloat
% command of the stfloats package.

\section{Conclusion}

We propose the multiple random projection tree (MRPT) algorithm for approximate $k$-nearest neighbor search in high dimensions. The method is based on combining multiple sparse random projection trees using a novel voting scheme where the final search is focused to points occurring most frequently among points retrieved by multiple trees. The algorithm is straightforward to implement and exploits standard fast linear algebra libraries by combining calculations into large matrix--matrix operations. 

We demonstrate through extensive experiments on both real and simulated data that the proposed method is faster than state-of-the-art space-partitioning tree and hashing based algorithms on a variety of accuracy levels. It is also faster than a leading graph based algorithm (KGraph) on high accuracy levels, while being slightly slower on low accuracy levels. The good performance of MRPT is especially pronounced for high-dimensional data sets.

Due to its very competitive and consistent performance, and simple and efficient index construction stage --- especially compared to graph-based algorithms --- the proposed MRPT method is an ideal method for a wide variety of applications where high-dimensional large data sets are involved.

\begin{figure*}[hbtp]
\captionsetup[subfigure]{labelformat=empty}
\centering
\includegraphics[width=\textwidth]{comparison.png}

\caption{Query time (log scale) for 100 queries versus recall on six different data sets with $k = 10$. The compared methods are: k-d tree (ann), multi-probe LSH (falconn), randomized k-d tree (flann-kd), hierarchical k-means tree (flann-kmeans), $k$-NN graph (kgraph) and our method (mrpt). We also include results for our implementation of classic RP-trees (rp trees) and RP-trees with sparsity (sparse rp trees).}
\label{fig:comparison}
\end{figure*}


\begin{figure*}[!ht]
\captionsetup[subfigure]{labelformat=empty}
\centering
\includegraphics[width=\textwidth]{different_ks.png}
\caption{Results for $k \in \{1,10,25,100\}$ on the MNIST data set.}
\label{fig:diff_k}
\end{figure*}

\begin{table*}[t]
   \centering
   \caption{Query times required to reach at least a given recall level for recall R = 80\%, 90\%, 95\%, 99\%, and for $k=10$. Fastest times for each recall level are emphasized with bold font. MRPT consistently outperforms other algorithms at high recall levels. For example, for R $\ge$ 95\%, 71.4\% of times MPRT achieves the fastest query time (i.e., 10 out of 14) whereas it is only 28.6\% for KGraph (i.e., 4 out of 14).
   ${}^{*)}$ R = 100\% by definition for brute force search.}
   \label{table:results}
\begin{tabular}%{|c|c|c|c|c|c|c|c|c|c|c|}
{ccccccccccc}
    \toprule
    & & 
    \multicolumn{9}{c}{time (s), 100 queries} \\
    \midrule
    data set & R (\%) & ANN & FALCONN & FLANN-kd & FLANN-kmeans & KGraph & RP-trees & MRPT($v=1$) & MRPT & brute force${}^{*)}$ \\
    \midrule
    \multirow{4}{*}{MNIST}
    & 80 & 0.89 & 0.21 & {\bf 0.02} & {\bf 0.02} & {\bf 0.02} & 0.09 & 0.05 & {\bf 0.02} & \multirow{4}{*}{2.59} \\
    & 90 & 1.57 & 0.36 & 0.04 & 0.04 & 0.04 & 0.16 & 0.08 & {\bf 0.03} & \multicolumn{1}{c}{} \\
    & 95 & 2.29 & 0.55 & 0.06 & 0.05 & 0.06 & 0.2 & 0.14 & {\bf 0.04} & \multicolumn{1}{c}{} \\
    & 99 & 3.46 & 1.23 & 0.15 & 0.1 & 0.44 & 0.39 & 0.22 & {\bf 0.07} & \multicolumn{1}{c}{} \\
    \midrule
    \multirow{4}{*}{News}
    & 80 & 2.15 & 0.39 & 0.56 & 0.26 & {\bf 0.05} & 0.71 & 0.64 & 0.12 & \multirow{4}{*}{23.09} \\
    & 90 & 4.42 & 0.88 & 2.92 & 0.45 & {\bf 0.11} & 1.37 & 1.05 & 0.2 & \multicolumn{1}{c}{} \\
    & 95 & 7.35 & 1.23 & 9.55 & 1.39 & {\bf 0.25} & 1.98 & 1.62 & 0.28 & \multicolumn{1}{c}{} \\
    & 99 & 11.83 & 2.58 & 38.81 & 7.5 & 3.24 & 4.3 & 3.87 & {\bf 0.5} & \multicolumn{1}{c}{} \\
    \midrule
    \multirow{4}{*}{GIST}
    & 80 & 20.54 & 21.21 & 6.03 & 1.86 & {\bf 0.59} & 4.77 & 4.19 & 1.03 & \multirow{4}{*}{52.98} \\
    & 90 & 36.18 & 29.19 & 13.35 & 3.63 & 1.87 & 7.69 & 7.54 & {\bf 1.83} & \multicolumn{1}{c}{} \\
    & 95 & 48.13 & 37.92 & 24.41 & 6.22 & 7.94 & 13.6 & 11.89 & {\bf 2.91} & \multicolumn{1}{c}{} \\
    & 99 & 76.57 & 50.29 & 59.12 & 14.14 & - & 24.11 & 20.33 & {\bf 6.1} & \multicolumn{1}{c}{} \\
    \midrule
    \multirow{4}{*}{SIFT}
    & 80 & 0.81 & 0.93 & 0.21 & 0.09 & {\bf 0.05} & 0.4 & 0.38 & 0.24 & \multirow{4}{*}{21.32} \\
    & 90 & 1.67 & 1.47 & 0.41 & 0.18 & {\bf 0.11} & 0.71 & 0.59 & 0.31 & \multicolumn{1}{c}{} \\
    & 95 & 3.1 & 2.06 & 0.84 & 0.31 & {\bf 0.18} & 0.92 & 0.97 & 0.37 & \multicolumn{1}{c}{} \\
    & 99 & 7.88 & 3.2 & 2.81 & 0.86 & {\bf 0.35} & 3.22 & 2.16 & 0.62 & \multicolumn{1}{c}{} \\
    \midrule
    \multirow{4}{*}{Trevi}
    & 80 & 5.5 & 14.62 & 0.97 & 0.7 & {\bf 0.16} & 1.48 & 0.94 & 0.27 & \multirow{4}{*}{22.15} \\
    & 90 & 12.43 & 17.88 & 1.81 & 1.01 & 0.75 & 2.44 & 1.46 & {\bf 0.45} & \multicolumn{1}{c}{} \\
    & 95 & 18.64 & 18.8 & 3.02 & 1.69 & 18.68 & 3.74 & 2.15 & {\bf 0.67} & \multicolumn{1}{c}{} \\
    & 99 & - & 25.54 & 10.54 & 7.61 & - & 8.13 & 5.25 & {\bf 2.07} & \multicolumn{1}{c}{} \\
    \midrule
    \multirow{4}{*}{STL-10}
    & 80 & 31.31 & 18.62 & 4.29 & 2.04 & {\bf 0.18} & 5.46 & 3.06 & 0.88 & \multirow{4}{*}{49.43} \\
    & 90 & 39.38 & 28.62 & 8.76 & 2.72 & {\bf 0.39} & 9.45 & 6.32 & 1.52 & \multicolumn{1}{c}{} \\
    & 95 & 45.0 & 31.46 & 13.39 & 3.53 & {\bf 1.13} & 11.62 & 8.58 & 2.43 & \multicolumn{1}{c}{} \\
    & 99 & 49.67 & 45.42 & 29.96 & 6.1 & 32.12 & 18.56 & 16.89 & {\bf 4.28} & \multicolumn{1}{c}{} \\
    \midrule
    \multirow{4}{*}{Random}
    & 80 & 23.86 & 10.1 & 12.55 & 9.73 & 7.59 & 8.55 & 8.2 & {\bf 4.55} & \multirow{4}{*}{10.9} \\
    & 90 & 27.33 & 12.27 & 13.91 & 11.48 & 9.8 & 10.2 & 9.8 & {\bf 6.9} & \multicolumn{1}{c}{} \\
    & 95 & 29.07 & 14.28 & 14.5 & 12.55 & 10.93 & 10.88 & 10.71 & {\bf 8.6} & \multicolumn{1}{c}{}\\
    & 99 & 30.47 & 17.03 & 15.37 & 13.7 & 12.02 & 11.59 & 11.42 & {\bf 10.36} & \multicolumn{1}{c}{} \\
    \bottomrule
\end{tabular}
\end{table*}

% trigger a \newpage just before the given reference
% number - used to balance the columns on the last page
% adjust value as needed - may need to be readjusted if
% the document is modified later
%\IEEEtriggeratref{8}
% The "triggered" command can be changed if desired:
%\IEEEtriggercmd{\enlargethispage{-5in}}

% references section

% can use a bibliography generated by BibTeX as a .bbl file
% BibTeX documentation can be easily obtained at:
% http://www.ctan.org/tex-archive/biblio/bibtex/contrib/doc/
% The IEEEtran BibTeX style support page is at:
% http://www.michaelshell.org/tex/ieeetran/bibtex/
%\bibliographystyle{IEEEtran}
% argument is your BibTeX string definitions and bibliography database(s)
%\bibliography{IEEEabrv,../bib/paper}
%
% <OR> manually copy in the resultant .bbl file
% set second argument of \begin to the number of references
% (used to reserve space for the reference number labels box)

\balance

\bibliographystyle{IEEEtran}
\bibliography{RP_trees_BigData_refs}  

\end{document}



RP_trees_arxiv.bbl
% Generated by IEEEtran.bst, version: 1.13 (2008/09/30)
\begin{thebibliography}{10}
\providecommand{\url}[1]{#1}
\csname url@samestyle\endcsname
\providecommand{\newblock}{\relax}
\providecommand{\bibinfo}[2]{#2}
\providecommand{\BIBentrySTDinterwordspacing}{\spaceskip=0pt\relax}
\providecommand{\BIBentryALTinterwordstretchfactor}{4}
\providecommand{\BIBentryALTinterwordspacing}{\spaceskip=\fontdimen2\font plus
\BIBentryALTinterwordstretchfactor\fontdimen3\font minus
  \fontdimen4\font\relax}
\providecommand{\BIBforeignlanguage}[2]{{%
\expandafter\ifx\csname l@#1\endcsname\relax
\typeout{** WARNING: IEEEtran.bst: No hyphenation pattern has been}%
\typeout{** loaded for the language `#1'. Using the pattern for}%
\typeout{** the default language instead.}%
\else
\language=\csname l@#1\endcsname
\fi
#2}}
\providecommand{\BIBdecl}{\relax}
\BIBdecl

\bibitem{Muja2014}
M.~Muja and D.~G. Lowe, ``Scalable nearest neighbor algorithms for high
  dimensional data,'' \emph{Pattern Analysis and Machine Intelligence, IEEE
  Transactions on}, vol.~36, no.~11, pp. 2227--2240, 2014.

\bibitem{Nister2006}
D.~Nister and H.~Stewenius, ``Scalable recognition with a vocabulary tree,'' in
  \emph{2006 IEEE Computer Society Conference on Computer Vision and Pattern
  Recognition (CVPR'06)}, vol.~2.\hskip 1em plus 0.5em minus 0.4em\relax IEEE,
  2006, pp. 2161--2168.

\bibitem{Lowe2004}
D.~G. Lowe, ``Distinctive image features from scale-invariant keypoints,''
  \emph{International journal of computer vision}, vol.~60, no.~2, pp. 91--110,
  2004.

\bibitem{Amit1997}
Y.~Amit and D.~Geman, ``Shape quantization and recognition with randomized
  trees,'' \emph{Neural computation}, vol.~9, no.~7, pp. 1545--1588, 1997.

\bibitem{Hays2007}
J.~Hays and A.~A. Efros, ``Scene completion using millions of photographs,'' in
  \emph{ACM Transactions on Graphics (TOG)}, vol.~26, no.~3.\hskip 1em plus
  0.5em minus 0.4em\relax ACM, 2007, p.~4.

\bibitem{Wang2016}
L.~Wang, S.~Tasoulis, T.~Roos, and J.~Kangasharju, ``Kvasir: Scalable provision
  of semantically relevant web content on big data framework,'' \emph{IEEE
  Transactions on Big Data}, vol. Advance online publication.
  DOI:10.1109/TBDATA.2016.2557348, 2016.

\bibitem{Indyk1998}
P.~Indyk and R.~Motwani, ``Approximate nearest neighbors: towards removing the
  curse of dimensionality,'' in \emph{Proceedings of the thirtieth annual ACM
  symposium on Theory of computing}.\hskip 1em plus 0.5em minus 0.4em\relax
  ACM, 1998, pp. 604--613.

\bibitem{Gionis1999}
A.~Gionis, P.~Indyk, R.~Motwani \emph{et~al.}, ``Similarity search in high
  dimensions via hashing,'' in \emph{VLDB}, vol.~99, no.~6, 1999, pp. 518--529.

\bibitem{Andoni2006}
A.~Andoni and P.~Indyk, ``Near-optimal hashing algorithms for approximate
  nearest neighbor in high dimensions,'' in \emph{Foundations of Computer
  Science, 2006. FOCS'06. 47th Annual IEEE Symposium on}.\hskip 1em plus 0.5em
  minus 0.4em\relax IEEE, 2006, pp. 459--468.

\bibitem{Lv2007}
Q.~Lv, W.~Josephson, Z.~Wang, M.~Charikar, and K.~Li, ``Multi-probe lsh:
  efficient indexing for high-dimensional similarity search,'' in
  \emph{Proceedings of the 33rd international conference on Very large data
  bases}.\hskip 1em plus 0.5em minus 0.4em\relax VLDB Endowment, 2007, pp.
  950--961.

\bibitem{Dong2008}
W.~Dong, Z.~Wang, W.~Josephson, M.~Charikar, and K.~Li, ``Modeling lsh for
  performance tuning,'' in \emph{Proceedings of the 17th ACM conference on
  Information and knowledge management}.\hskip 1em plus 0.5em minus 0.4em\relax
  ACM, 2008, pp. 669--678.

\bibitem{Andoni2015}
A.~Andoni, P.~Indyk, T.~Laarhoven, I.~Razenshteyn, and L.~Schmidt, ``Practical
  and optimal lsh for angular distance,'' in \emph{Advances in Neural
  Information Processing Systems 28}.\hskip 1em plus 0.5em minus 0.4em\relax
  Curran Associates, Inc., 2015, pp. 1225--1233.

\bibitem{Hajebi2011}
K.~Hajebi, Y.~Abbasi-Yadkori, H.~Shahbazi, and H.~Zhang, ``Fast approximate
  nearest-neighbor search with k-nearest neighbor graph,'' in \emph{IJCAI
  Proceedings-International Joint Conference on Artificial Intelligence},
  vol.~22, no.~1, 2011, p. 1312.

\bibitem{Dong2011}
W.~Dong, C.~Moses, and K.~Li, ``Efficient k-nearest neighbor graph construction
  for generic similarity measures,'' in \emph{Proceedings of the 20th
  international conference on World wide web}.\hskip 1em plus 0.5em minus
  0.4em\relax ACM, 2011, pp. 577--586.

\bibitem{Wang2012}
J.~Wang, J.~Wang, G.~Zeng, Z.~Tu, R.~Gan, and S.~Li, ``Scalable k-nn graph
  construction for visual descriptors,'' in \emph{Computer Vision and Pattern
  Recognition (CVPR), 2012 IEEE Conference on}.\hskip 1em plus 0.5em minus
  0.4em\relax IEEE, 2012, pp. 1106--1113.

\bibitem{Bentley1975}
J.~L. Bentley, ``Multidimensional binary search trees used for associative
  searching,'' \emph{Communications of the ACM}, vol.~18, no.~9, pp. 509--517,
  1975.

\bibitem{Arya1998}
S.~Arya, D.~M. Mount, N.~S. Netanyahu, R.~Silverman, and A.~Y. Wu, ``An optimal
  algorithm for approximate nearest neighbor searching fixed dimensions,''
  \emph{Journal of the ACM (JACM)}, vol.~45, no.~6, pp. 891--923, 1998.

\bibitem{Fukunaga1975}
K.~Fukunaga and P.~M. Narendra, ``A branch and bound algorithm for computing
  k-nearest neighbors,'' \emph{IEEE transactions on computers}, vol. 100,
  no.~7, pp. 750--753, 1975.

\bibitem{Beygelzimer2006}
A.~Beygelzimer, S.~Kakade, and J.~Langford, ``Cover trees for nearest
  neighbor,'' in \emph{Proceedings of the 23rd international conference on
  Machine learning}.\hskip 1em plus 0.5em minus 0.4em\relax ACM, 2006, pp.
  97--104.

\bibitem{Yianilos1993}
P.~N. Yianilos, ``Data structures and algorithms for nearest neighbor search in
  general metric spaces,'' in \emph{SODA}, vol.~93, no. 194, 1993, pp. 311--21.

\bibitem{Leibe2006}
B.~Leibe, K.~Mikolajczyk, and B.~Schiele, ``Efficient clustering and matching
  for object class recognition.'' in \emph{BMVC}, 2006, pp. 789--798.

\bibitem{Silpa2008}
C.~Silpa-Anan and R.~Hartley, ``Optimised kd-trees for fast image descriptor
  matching,'' in \emph{Computer Vision and Pattern Recognition, 2008. CVPR
  2008. IEEE Conference on}.\hskip 1em plus 0.5em minus 0.4em\relax IEEE, 2008,
  pp. 1--8.

\bibitem{Muja2009}
M.~Muja and D.~G. Lowe, ``Fast approximate nearest neighbors with automatic
  algorithm configuration.'' \emph{VISAPP (1)}, vol.~2, pp. 331--340, 2009.

\bibitem{Dasgupta2009}
S.~Dasgupta and Y.~Freund, ``Random projection trees for vector quantization,''
  \emph{IEEE Transactions on Information Theory}, vol.~55, no.~7, pp.
  3229--3242, 2009.

\bibitem{Dasgupta2013}
S.~Dasgupta and K.~Sinha, ``Randomized partition trees for nearest neighbor
  search,'' \emph{Algorithmica}, vol.~72, no.~1, pp. 237--263, 2015.

\bibitem{Achlioptas2001}
D.~Achlioptas, ``Database-friendly random projections,'' in \emph{Proceedings
  of the twentieth ACM SIGMOD-SIGACT-SIGART symposium on Principles of database
  systems}.\hskip 1em plus 0.5em minus 0.4em\relax ACM, 2001, pp. 274--281.

\bibitem{Li2006}
P.~Li, T.~J. Hastie, and K.~W. Church, ``Very sparse random projections,'' in
  \emph{Proceedings of the 12th ACM SIGKDD international conference on
  Knowledge discovery and data mining}.\hskip 1em plus 0.5em minus 0.4em\relax
  ACM, 2006, pp. 287--296.

\bibitem{Jegou2011A}
H.~J{\'e}gou, M.~Douze, and C.~Schmid, ``Product quantization for nearest
  neighbor search,'' \emph{IEEE Transactions on Pattern Analysis and Machine
  Intelligence}, vol.~33, no.~1, pp. 117--128, 2011.

\bibitem{Jegou2011B}
H.~Jegou, R.~Tavenard, M.~Douze, and L.~Amsaleg, ``Searching in one billion
  vectors: Re-rank with source coding,'' in \emph{Proceedings of the {IEEE}
  International Conference on Acoustics, Speech, and Signal Processing,
  {ICASSP} 2011}, 2011, pp. 861--864.

\bibitem{Winder2007}
S.~A.~J. Winder and M.~A. Brown, ``Learning local image descriptors,'' in
  \emph{2007 {IEEE} Computer Society Conference on Computer Vision and Pattern
  Recognition ({CVPR} 2007)}, 2007.

\bibitem{Coates2011}
A.~Coates, A.~Y. Ng, and H.~Lee, ``An analysis of single-layer networks in
  unsupervised feature learning,'' in \emph{Proceedings of the Fourteenth
  International Conference on Artificial Intelligence and Statistics, {AISTATS}
  2011}, 2011, pp. 215--223.

\end{thebibliography}
