Product Quantization for Nearest Neighbor Search
Hervé Jégou, Matthijs Douze, Cordelia Schmid

To cite this version:
Hervé Jégou, Matthijs Douze, Cordelia Schmid. Product Quantization for Nearest Neighbor
Search. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2011, 33 (1), pp.117-128.
�10.1109/TPAMI.2010.57�. �inria-00514462v2�

HAL Id: inria-00514462
https://inria.hal.science/inria-00514462v2
Submitted on 22 May 2013

HAL is a multi-disciplinary open accessarchive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come fromteaching and research institutions in France orabroad, or from public or private research centers.

L’archive ouverte pluridisciplinaire HAL, estdestinée au dépôt et à la diffusion de documentsscientifiques de niveau recherche, publiés ou non,émanant des établissements d’enseignement et derecherche français ou étrangers, des laboratoires
publics ou privés.

1

Product quantization for nearest neighbor search
Hervé Jégou, Matthijs Douze, Cordelia Schmid
INRIA

Abstract— This paper introduces a product quantization based
approach for approximate nearest neighbor search. The idea
is to decomposes the space into a Cartesian product of low
dimensional subspaces and to quantize each subspace separately.
A vector is represented by a short code composed of its subspacequantization indices. The Euclidean distance between two vectors
can be efficiently estimated from their codes. An asymmetricversion increases precision, as it computes the approximate
distance between a vector and a code.
Experimental results show that our approach searches fornearest neighbors efficiently, in particular in combination with aninverted file system. Results for SIFT and GIST image descriptors
show excellent search accuracy outperforming three state-of-theart approaches. The scalability of our approach is validated ona dataset of two billion vectors.
Index Terms— High-dimensional indexing, image indexing,
very large databases, approximate search.

I. I NTRODUCTIONOMPUTING Euclidean distances between high dimensional vectors is a fundamental requirement in many
applications. It is used, in particular, for nearest neighbor(NN) search. Nearest neighbor search is inherently expensive
due to the curse of dimensionality [3], [4]. Focusing on the
D-dimensional Euclidean space RD , the problem is to findthe element NN(x), in a finite set Y ⊂ RD of n vectors,
minimizing the distance to the query vector x ∈ RD :C

NN(x) = arg min d(x, y).
y∈Y

(1)

Several multi-dimensional indexing methods, such as the
popular KD-tree [5] or other branch and bound techniques,
have been proposed to reduce the search time. However, forhigh dimensions it turns out [6] that such approaches are
not more efficient than the brute-force exhaustive distancecalculation, whose complexity is O(nD).
There is a large body of literature [7], [8], [9] on algorithms
that overcome this issue by performing approximate nearest
neighbor (ANN) search. The key idea shared by these algorithms is to find the NN with high probability “only”, insteadof probability 1. Most of the effort has been devoted to the
Euclidean distance, though recent generalizations have been
proposed for other metrics [10]. In this paper, we consider the
Euclidean distance, which is relevant for many applications.
In this case, one of the most popular ANN algorithms isthe Euclidean Locality-Sensitive Hashing (E2LSH) [7], [11],which provides theoretical guarantees on the search quality
This work was partly realized as part of the Quaero Programme, fundedby OSEO, French State agency for innovation. It was originally published asa technical report [1] in August 2009. It is also related to the work [2] onsource coding for nearest neighbor search. This particular version is a postprint correcting mistakes reported by Jan Stria and Ondrej Chum (CVUT).

with limited assumptions. It has been successfully used for
local descriptors [12] and 3D object indexing [13], [11].
However, for real data, LSH is outperformed by heuristic
methods, which exploit the distribution of the vectors. Thesemethods include randomized KD-trees [14] and hierarchicalk-means [15], both of which are implemented in the FLANN
selection algorithm [9].
ANN algorithms are typically compared based on the tradeoff between search quality and efficiency. However, this tradeoff does not take into account the memory requirements of
the indexing structure. In the case of E2LSH, the memory
usage may even be higher than that of the original vectors.
Moreover, both E2LSH and FLANN need to perform a final
re-ranking step based on exact L2 distances, which requires theindexed vectors to be stored in main memory if access speed
is important. This constraint seriously limits the number ofvectors that can be handled by these algorithms. Only recently,
researchers came up with methods limiting the memory usage.
This is a key criterion for problems involving large amountsof data [16], i.e., in large-scale scene recognition [17], where
millions to billions of images have to be indexed. In [17],Torralba et al. represent an image by a single global GIST
descriptor [18] which is mapped to a short binary code.
When no supervision is used, this mapping is learned suchthat the neighborhood in the embedded space defined by the
Hamming distance reflects the neighborhood in the Euclidean
space of the original features. The search of the Euclidean
nearest neighbors is then approximated by the search of thenearest neighbors in terms of Hamming distances between
codes. In [19], spectral hashing (SH) is shown to outperform the binary codes generated by the restricted Boltzmann
machine [17], boosting and LSH. Similarly, the Hamming
embedding method of Jegou et al. [20], [21] uses a binary
signature to refine quantized SIFT or GIST descriptors in abag-of-features image search framework.
In this paper, we construct short codes using quantization.
The goal is to estimate distances using vector-to-centroid
distances, i.e., the query vector is not quantized, codes are
assigned to the database vectors only. This reduces the quantization noise and subsequently improves the search quality.
To obtain precise distances, the quantization error must be
limited. Therefore, the total number k of centroids should
be sufficiently large, e.g., k = 264 for 64-bit codes. This
raises several issues on how to learn the codebook and assign
a vector. First, the number of samples required to learn thequantizer is huge, i.e., several times k. Second, the complexityof the algorithm itself is prohibitive. Finally, the amount ofcomputer memory available on Earth is not sufficient to store
the floating point values representing the centroids.
The hierarchical k-means see (HKM) improves the effi-

2ciency of the learning stage and of the corresponding assignment procedure [15]. However, the aforementioned limitations
still apply, in particular with respect to memory usage and sizeof the learning set. Another possibility are scalar quantizers,
but they offer poor quantization error properties in termsof the trade-off between memory and reconstruction error.
Lattice quantizers offer better quantization properties for uniform vector distributions, but this condition is rarely satisfied
by real world vectors. In practice, these quantizers perform
significantly worse than k-means in indexing tasks [22]. In
this paper, we focus on product quantizers. To our knowledge,
such a semi-structured quantizer has never been considered inany nearest neighbor search method.
The advantages of our method are twofold. First, the numberof possible distances is significantly higher than for competing
Hamming embedding methods [20], [17], [19], as the Hamming space used in these techniques allows for a few distinct
distances only. Second, as a byproduct of the method, we
get an estimation of the expected squared distance, which is
required for ε-radius search or for using Lowe’s distance ratio
criterion [23]. The motivation of using the Hamming spacein [20], [17], [19] is to compute distances efficiently. Note,however, that one of the fastest ways to compute Hamming
distances consists in using table lookups. Our method uses
a similar number of table lookups, resulting in comparable
efficiency.
An exhaustive comparison of the query vector with all codesis prohibitive for very large datasets. We, therefore, introduce
a modified inverted file structure to rapidly access the most
relevant vectors. A coarse quantizer is used to implement thisinverted file structure, where vectors corresponding to a cluster(index) are stored in the associated list. The vectors in thelist are represented by short codes, computed by our productquantizer, which is used here to encode the residual vector
with respect to the cluster center.
The interest of our method is validated on two kindsof vectors, namely local SIFT [23] and global GIST [18]descriptors. A comparison with the state of the art showsthat our approach outperforms existing techniques, in particular spectral hashing [19], Hamming embedding [20] andFLANN [9].
Our paper is organized as follows. Section II introduces
the notations for quantization as well as the product quantizer
used by our method. Section III presents our approach for
NN search and Section IV introduces the structure used to
avoid exhaustive search. An evaluation of the parameters of
our approach and a comparison with the state of the art isgiven in Section V.
II. BACKGROUND : QUANTIZATION , PRODUCT QUANTIZER
A large body of literature is available on vector quantization, see [24] for a survey. In this section, we restrict ourpresentation to the notations and concepts used in the rest of
the paper.
A. Vector quantization
Quantization is a destructive process which has been extensively studied in information theory [24]. Its purpose is to

reduce the cardinality of the representation space, in particularwhen the input data is real-valued. Formally, a quantizer is afunction q mapping a D-dimensional vector x ∈ RD to avector q(x) ∈ C = {ci ; i ∈ I}, where the index set I is fromnow on assumed to be finite: I = 0 . . . k −1. The reproduction
values ci are called centroids. The set of reproduction values
C is the codebook of size k.
The set Vi of vectors mapped to a given index i is referred
to as a (Voronoi) cell, and defined as
Vi , {x ∈ RD : q(x) = ci }.

(2)

The k cells of a quantizer form a partition of RD . By
definition, all the vectors lying in the same cell Vi are reconstructed by the same centroid ci . The quality of a quantizer
is usually measured by the mean squared error between the
input vector x and its reproduction value q(x):Z
2
MSE(q) = EX d(q(x), x)2 = p(x) d q(x), x dx, (3)where d(x, y) = ||x − y|| is the Euclidean distance between
x and y, and where p(x) is the probability distribution function corresponding the random variable X. For an arbitrary
probability distribution function, Equation 3 is numericallycomputed using Monte-Carlo sampling, as the average of||q(x) − x||2 on a large set of samples.
In order for the quantizer to be optimal, it has to satisfy two
properties known as the Lloyd optimality conditions. First, avector x must be quantized to its nearest codebook centroid,in terms of the Euclidean distance:q(x) = arg min d(x, ci ).
ci ∈C

(4)

As a result, the cells are delimited by hyperplanes. The
second Lloyd condition is that the reconstruction value must
be the expectation of the vectors lying in the Voronoi cell:Z 
ci = EX x|i =p(x|x ∈ Vi ) x dx.
(5)
Vi

The Lloyd quantizer, which corresponds to the k-meansclustering algorithm, finds a near-optimal codebook by iteratively assigning the vectors of a training set to centroids andre-estimating these centroids from the assigned vectors. In thefollowing, we assume that the two Lloyd conditions hold, as
we learn the quantizer using k-means. Note, however, that kmeans does only find a local optimum in terms of quantization
error.
Another quantity that will be used in the following is themean squared distortion ξ(q, ci ) obtained when reconstructing
a vector of a cell Vi by the corresponding centroid ci . Denoting
by pi = P q(x) = ci the probability that a vector is assigned
to the centroid ci , it is computed asZ2
1
d x, q(x) p(x) dx.
(6)
ξ(q, ci ) =pi ViNote that the MSE can be obtained from these quantities as
X
MSE(q) =pi ξ(q, ci ).
(7)
i∈I

3

The memory cost of storing the index value, without any
further processing (entropy coding), is dlog2 ke bits. Therefore,
it is convenient to use a power of two for k, as the code
produced by the quantizer is stored in a binary memory.

Let us consider a 128-dimensional vector, for example the
SIFT descriptor [23]. A quantizer producing 64-bits codes,
i.e., “only” 0.5 bit per component, contains k = 264 centroids.
Therefore, it is impossible to use Lloyd’s algorithm or even
HKM, as the number of samples required and the complexityof learning the quantizer are several times k. It is even
impossible to store the D×k floating point values representing
the k centroids.
Product quantization is an efficient solution to address theseissues. It is a common technique in source coding, whichallows to choose the number of components to be quantizedjointly (for instance, groups of 24 components can be quantized using the powerful Leech lattice). The input vector x is
split into m distinct subvectors uj , 1 ≤ j ≤ m of dimension
D∗ = D/m, where D is a multiple of m. The subvectors
are quantized separately using m distinct quantizers. A given
vector x is therefore mapped as follows:
x , ..., x ∗ , ..., xD−D∗ +1 , ..., xD → q1 u1 (x)), ..., qm (um (x) ,| 1 {z D }{z}|u1 (x)

um (x)(8)
where qj is a low-complexity quantizer associated with thej th subvector. With the subquantizer qj we associate the index
set Ij , the codebook Cj and the corresponding reproduction
values cj,i .
A reproduction value of the product quantizer is identified
by an element of the product index set I = I1 × . . . × Im .
The codebook is therefore defined as the Cartesian product
C = C1 × . . . × C m ,(9)and a centroid of this set is the concatenation of centroidsof the m subquantizers. From now on, we assume that all
subquantizers have the same finite number k ∗ of reproduction
values. In that case, the total number of centroids is given byk = (k ∗ )m .

(10)

Note that in the extremal case where m = D, the components of a vector x are all quantized separately. Then the
product quantizer turns out to be a scalar quantizer, where thequantization function associated with each component may bedifferent.
The strength of a product quantizer is to produce a large
set of centroids from several small sets of centroids: those
associated with the subquantizers. When learning the subquantizers using Lloyd’s algorithm, a limited number of vectors isused, but the codebook is, to some extent, still adapted to thedata distribution to represent. The complexity of learning thequantizer is m times the complexity of performing k-meansclustering with k ∗ centroids of dimension D∗ .
Storing the codebook C explicitly is not efficient. Instead,
we store the m × k ∗ centroids of all the subquantizers, i.e.,

assignment complexity
kD
lD

m k∗ D∗ = k1/m D

m k∗ D∗ = k1/m D

fproduct k-means

TABLE IM EMORY USAGE OF THE CODEBOOK AND ASSIGNMENT COMPLEXITY FOR
DIFFERENT QUANTIZERS . HKM IS PARAMETRIZED BY TREE HEIGHT l
AND THE BRANCHING FACTOR bf .

0.3
m=1
m=2
m=4
m=8
m=16

0.25
square distortion D(q)B. Product quantizers

memory usagekDbf(k − 1) Db −1

k-means
HKM

k*=16

0.2
0.15

64

0.1

256
1024

0.05
0
0

Fig. 1.

16

32

64
96
code length (bits)

128

160

SIFT: quantization error associated with the parameters m and k∗ .

m D∗ k ∗ = k ∗ D floating points values. Quantizing an element
requires k ∗ D floating point operations. Table I summarizes
the resource requirements associated with k-means, HKM andproduct k-means. The product quantizer is clearly the the only
one that can be indexed in memory for large values of k.
In order to provide good quantization properties whenchoosing a constant value of k ∗ , each subvector should have,on average, a comparable energy. One way to ensure this property is to multiply the vector by a random orthogonal matrix
prior to quantization. However, for most vector types this is not
required and not recommended, as consecutive components
are often correlated by construction and are better quantized
together with the same subquantizer. As the subspaces are
orthogonal, the squared distortion associated with the productquantizer isX
MSE(q) =MSE(qj ),(11)
jwhere MSE(qj ) is the distortion associated with quantizer qj .
Figure 1 shows the MSE as a function of the code lengthfor different (m,k ∗ ) tuples, where the code length is l =m log2 k ∗ , if k ∗ is a power of two. The curves are obtained
for a set of 128-dimensional SIFT descriptors, see section V
for details. One can observe that for a fixed number of bits,
it is better to use a small number of subquantizers with many
centroids than having many subquantizers with few bits. At
the extreme when m = 1, the product quantizer becomes aregular k-means codebook.
High values of k ∗ increase the computational cost of thequantizer, as shown by Table I. They also increase the memory

4

encoding x
compute d uj (x), cj,iˆ y) or d(x,˜ y)for y ∈ Y, compute d(x,find the k smallest distances

SDC

ADC

k∗ D
0

0
k∗ D

nm
nm
n + k log k log log n

TABLE IIA LGORITHM AND COMPUTATIONAL COSTS ASSOCIATED WITH SEARCHING
THE k NEAREST NEIGHBORS USING THE PRODUCT QUANTIZER FOR
SYMMETRIC AND ASYMMETRIC DISTANCE COMPUTATIONS (SDC, ADC).

symmetric caseasymmetric caseFig. 2. Illustration of the symmetric and asymmetric distance computation.
The distance d(x, y) is estimated with either the distance d(q(x), q(y)) (left)or the distance d(x, q(y)) (right). The mean squared error on the distance ison average bounded by the quantization error.

usage of storing the centroids (k ∗ × D floating point values),
which further reduces the efficiency if the centroid look-up
table does no longer fit in cache memory. In the case wherem = 1, we can not afford using more than 16 bits to keep
this cost tractable. Using k ∗ = 256 and m = 8 is often areasonable choice.
III. S EARCHING WITH QUANTIZATION
Nearest neighbor search depends on the distances between
the query vector and the database vectors, or equivalently
the squared distances. The method introduced in this section
compares the vectors based on their quantization indices,in the spirit of source coding techniques. We first explain
how the product quantizer properties are used to computethe distances. Then we provide a statistical bound on the
distance estimation error, and propose a refined estimator forthe squared Euclidean distance.
A. Computing distances using quantized codes
Let us consider the query vector x and a database vector
y. We propose two methods to compute an approximate
Euclidean distance between these vectors, a symmetric anda asymmetric one. See Figure 2 for an illustration.
Symmetric distance computation (SDC): both the vectors x
and y are represented by their respective centroids q(x) andq(y). The distance d(x, y) is approximated by the distanceˆ y) , d q(x), q(y) which is efficiently obtained using ad(x,product quantizer
sX
2

ˆ
d qj (x), qj (y) , (12)d(x, y) = d q(x), q(y) =j2where the distance d cj,i , cj,i0
is read from a look-up
table associated with the j th subquantizer. Each look-up tablecontains all the squared distances between pairs of centroids(i, i0 ) of the subquantizer, or (k ∗ )2 squared distances1 .
Asymmetric distance computation (ADC): the database
1 In fact, it is possible to store only k ∗ (k ∗ − 1)/2 pre-computed squared

distances, because this distance matrix is symmetric and the diagonal elements
are zeros.

vector y is represented by q(y), but the query x is not encoded.
˜The distance d(x, y) is approximated by the distance d(x, y) ,d x, q(y) , which is computed using the decomposition˜ y) = d x, q(y) =
d(x,

sX

2
d uj (x), qj (uj (y)) ,

(13)

j

2
where the squared distances d uj (x), cj,i : j = 1 . . . m, i =1 . . . k ∗ , are computed prior to the search.
For nearest neighbors search, we do not compute the squareroots in practice: the square root function is monotonicallyincreasing and the squared distances produces the same vector
ranking.
Table II summarizes the complexity of the different steps
involved in searching the k nearest neighbors of a vector x in adataset Y of n = |Y| vectors. One can see that SDC and ADC
have the same query preparation cost, which does not dependon the dataset size n. When n is large (n > k ∗ D∗ ), the mostconsuming operations are the summations in Equations 12
and 13. The complexity given in this table for searching thek smallest elements is the average complexity for n  k andwhen the elements are arbitrarily ordered ([25], Section 5.3.3,
Equation 17).
The only advantage of SDC over ADC is to limit the
memory usage associated with the queries, as the query vector
is defined by a code. This is most cases not relevant and one
should then use the asymmetric version, which obtains a lower
distance distortion for a similar complexity. We will focus on
ADC in the rest of this section.
B. Analysis of the distance error
In this subsection, we analyze the distance error when using
˜ y) instead of d(x, y). This analysis does not dependd(x,on the use of a product quantizer and is valid for any
quantizer satisfying Lloyd’s optimality conditions defined by
Equations 4 and 5 in Section II.
In the spirit of the mean squared error criterion used forreconstruction, the distance distortion is measured by the meansquared distance error (MSDE) on the distances:ZZ˜ y) 2 p(x)dx p(y)dy. (14)MSDE(q) ,d(x, y) − d(x,The triangular inequality givesd x, q(y) − d y, q(y) ≤ d(x, y) ≤ d x, q(y) + d y, q(y) ,
(15)

5

between x, which is fully known in our asymmetric distancecomputation method, and a random variable Y , subject toq(Y ) = q(y) = ci , which represents all the hypothesis on
y knowing its quantization index.
ẽ(x, y) , EY (x − Y )2 |q(Y ) = ci
(19)Z=(x − y)2 p(y|i) dy,
(20)
Vi
Z
1
(x − ci + ci − y)2 p(y) dy.
(21)
=
pi Vi

estimated distance

1.2

1

0.8

0.6

Developing the squared expression and observing, using
Lloyd’s condition of Equation 5, thatZ(y − ci ) p(y) dy = 0,(22)symmetric
asymmetric
0.4
0.4

0.6

0.8
1
true distances

1.2

Vi

Fig. 3. Typical query of a SIFT vector in a set of 1000 vectors: comparisonof the distance d(x, y) obtained with the SDC and ADC estimators. We haveused m = 8 and k∗ = 256, i.e., 64-bit code vectors. Best viewed in color.

Equation 21 simplifies toZ2ẽ(x, y) = x − q(y) +q(y) − y2p y|y ∈ Vi dy

Vi

(23)
and, equivalently,2
2d(x, y) − d(x, q(y)) ≤ d y, q(y) .

˜ y)2 + ξ q, q(y)= d(x,(16)

Combining this inequality with Equation 14, we obtain
Z

Z
2
MSDE(q) ≤ p(x)d y, q(y) p(y) dy dx(17)≤ MSE(q).

(18)

where MSE(q) is the mean squared error associated withquantizer q (Equation 3). This inequality, which holds forany quantizer, shows that the distance error of our method isstatistically bounded by the MSE associated with the quantizer.
For the symmetric version, a similar derivation shows that theerror is statistically bounded by 2×MSE(q). It is, therefore,
worth minimizing the quantization error, as this criterion
provides a statistical upper bound on the distance error. If an
exact distance calculation is performed on the highest ranked
vectors, as done in LSH [7], the quantization error can beused (instead of selecting an arbitrary set of elements) as acriterion to dynamically select the set of vectors on which thepost-processing should be applied.
C. Estimator of the squared distance
As shown later in this subsection, using the estimations dˆor d˜ leads to underestimate, on average, the distance between
descriptors. Figure 3 shows the distances obtained whenquerying a SIFT descriptor in a dataset of 1000 SIFT vectors.
It compares the true distance against the estimates computedwith Equations 12 and 13. One can clearly see the bias on these
distance estimators. Unsurprisingly, the symmetric version ismore sensitive to this bias.
Hereafter, we compute the expectation of the squared distance in order to cancel the bias. The approximation q(y) of agiven vector y is obtained, in the case of the product quantizer,
from the subquantizers indexes qj uj (y) , j = 1 . . . m. The
quantization index identifies the cells Vi in which y lies. We
can then compute the expected squared distance ẽ x, q(y)(24)
where we recognize the distortion ξ q, q(y) associated withthe reconstruction of y by its reproduction value.
Using the product quantizer and Equation 24, the computation of the expected squared distance between a vector x andthe vector y, for which we only know the quantization indices
qj uj (y) , consists in correcting Equation 13 asX˜ y)2 +ẽ(x, y) = d(x,ξj (y)(25)
j

where the correcting term, i.e., the average distortionξj (y) , ξ qj , qj uj (y)

(26)
th

associated with quantizing uj (y) to qj (y) using the j subquantizer, is learned and stored in a look-up table for allindexes of Ij .
Performing a similar derivation for the symmetric version,i.e., when both x and y are encoded using the productquantizer, we obtain the following corrected version of thesymmetric squared distance estimator:
XXˆ y)2 +ê(x, y) = d(x,ξj (x) +ξj 0 (y).
(27)
j

j0

Discussion: Figure 4 illustrates the probability distributionfunction of the difference between the true distance and theones estimated by Equations 13 and 25. It has been measured
on a large set of SIFT descriptors. The bias of the distanceestimation by Equation 13 is significantly reduced in thecorrected version. However, we observe that correcting thebias leads, in this case, to a higher variance of the estimator,
which is a common phenomenon in statistics. Moreover, forthe nearest neighbors, the correcting term is likely to behigher than the measure of Equation 13, which means that wepenalize the vectors with rare indexes. Note that the correcting
term is independent of the query in the asymmetric version.
In our experiments, we observe that the correction returns
inferior results on average. Therefore, we advocate the use of

empirical probability distribution function

6

In addition to the coarse quantizer, we adopt a strategy similar to that proposed in [20], i.e., the description of a vector isrefined by a code obtained with a product quantizer. However,in order to take into account the information provided by
the coarse quantizer, i.e., the centroid qc (y) associated witha vector y, the product quantizer qp is used to encode theresidual vectorr(y) = y − qc (y),(28)

biasedunbiased

0.1

0.05

0
-0.3

-0.2

-0.1
0
0.1
difference: estimator - d(x,y)

0.2

0.3

Fig. 4. PDF of the error on the distance estimation d − d˜ for the asymmetricmethod, evaluated on a set of 10000 SIFT vectors with m = 8 and k∗ = 256.
The bias (=-0.044) of the estimatord˜ is corrected (=0.002) with the errorquantization term ξ q, q(y) . However, the variance of the error increases
˜ = 0.00146.
with this correction: σ 2 (d − ẽ) = 0.00155 whereas σ 2 (d − d)Equation 13 for the nearest neighbor search. The correctedversion is useful only if we are interested in the distances
themselves.
IV. N ON EXHAUSTIVE SEARCH
Approximate nearest neighbor search with product quantizers is fast (only m additions are required per distance calculation) and reduces significantly the memory requirements forstoring the descriptors. Nevertheless, the search is exhaustive.
The method remains scalable in the context of a global image
description [17], [19]. However, if each image is described bya set of local descriptors, an exhaustive search is prohibitive,
as we need to index billions of descriptors and to perform
multiple queries [20].
To avoid exhaustive search we combine an inverted
file system [26] with the asymmetric distance computation (IVFADC). An inverted file quantizes the descriptors andstores image indices in the corresponding lists, see the step“coarse quantizer” in Figure 5. This allows rapid access to asmall fraction of image indices and was shown successful forvery large scale search [26]. Instead of storing an image index
only, we add a small code for each descriptor, as first donein [20]. Here, we encode the difference between the vector andits corresponding coarse centroid with a product quantizer, see
Figure 5. This approach significantly accelerates the search at
the cost of a few additional bits/bytes per descriptor. Furthermore, it slightly improves the search accuracy, as encoding theresidual is more precise than encoding the vector itself.
A. Coarse quantizer, locally defined product quantizer
Similar to the “Video-Google” approach [26], a codebook
is learned using k-means, producing a quantizer qc , referred toas the coarse quantizer in the following. For SIFT descriptors,
the number k 0 of centroids associated with qc typically ranges
from k 0 = 1 000 to k 0 = 1 000 000. It is therefore small
compared to that of the product quantizers used in Section III.

corresponding to the offset in the Voronoi cell. The energyof the residual vector is small compared to that of the vector
itself. The vector is approximated byÿ , qc (y) + qp y − qc (y) .
(29)
It is represented by the tuple qc (y), qp (r(y)) . By analogy
with the binary representation of a value, the coarse quantizer
provides the most significant bits, while the product quantizercode corresponds to the least significant bits.
The estimator of d(x, y), where x is the query and y the¨ y) between
database vector, is computed as the distance d(x,x and ÿ:¨ y) = d(x, ÿ) = d x − qc (y), qp y − qc (y) .
d(x,
(30)
Denoting by qp j the j th subquantizer, we use the following
decomposition to compute this estimator efficiently:
X 

 2
¨ y)2 =d(x,d uj x − qc (y) , qp j uj (y − qc (y)). (31)jSimilar to the ADC strategy, for each subquantizer qp j thedistances between the partial residual vector uj x − qc (y)and all the centroids cj,i of qp j are preliminarily computed
and stored.
The product quantizer is learned on a set of residual
vectors collected from a learning set. Although the vectors
are quantized to different indexes by the coarse quantizer, theresulting residual vectors are used to learn a unique productquantizer. We assume that the same product quantizer isaccurate when the distribution of the residual is marginalized
over all the Voronoi cells. This probably gives inferior results
to the approach consisting of learning and using a distinct
product quantizer per Voronoi cell. However, this would becomputationally expensive and would require storing k 0 product quantizer codebooks, i.e., k 0 ×d×k ∗ floating points values,
which would be memory-intractable for common values of k 0 .
B. Indexing structure
We use the coarse quantizer to implement an inverted filestructure as an array of lists L1 . . . Lk0 . If Y is the vector
dataset to index, the list Li associated with the centroid ci ofqc stores the set {y ∈ Y : qc (y) = ci }.
In inverted list Li , an entry corresponding to y contains avector identifier and the encoded residual qp (r(y)):
fieldidentifier
code

length (bits)
8–32
mdlog2 k∗ e

7Figure 5 gives an overview of how a database is indexed
and searched.
Indexing a vector y proceeds as follows:
1) quantize y to qc (y)2) compute the residual r(y) = y − qc (y)3) quantize r(y) to qp (r(y)), which, for the product quantizer, amounts to assigning uj (y) to qj (uj (y)), for j =1 . . . m.
4) add a new entry to the inverted list corresponding to
qc (y). It contains the vector (or image) identifier andthe binary code (the product quantizer’s indexes).
Searching the nearest neighbor(s) of a query x consists of1) quantize x to its w nearest neighbors in the codebook qc ;Fig. 5. Overview of the inverted file with asymmetric distance computation(IVFADC) indexing system. Top: insertion of a vector. Bottom: search.

The identifier field is the overhead due to the inverted filestructure. Depending on the nature of the vectors to be
stored, the identifier is not necessarily unique. For instance,
to describe images by local descriptors, image identifiers canreplace vector identifiers, i.e., all vectors of the same image
have the same identifier. Therefore, a 20-bit field is sufficient to
identify an image from a dataset of one million. This memory
cost can be reduced further using index compression [27], [28],which may reduce the average cost of storing the identifier
to about 8 bits, depending on parameters2 . Note that somegeometrical information can also be inserted in this entry, as
proposed in [20] and [27].

C. Search algorithm
The inverted file is the key to the non-exhaustive versionof our method. When searching the nearest neighbors of avector x, the inverted file provides a subset of Y for which
distances are estimated: only the inverted list Li corresponding
to qc (x) is scanned.
However, x and its nearest neighbor are often not quantized
to the same centroid, but to nearby ones. To address thisproblem, we use the multiple assignment strategy of [29]. The
query x is assigned to w indexes instead of only one, whichcorrespond to the w nearest neighbors of x in the codebookof qc . All the corresponding inverted lists are scanned. Multiple
assignment is not applied to database vectors, as this would
increase the memory usage.
2 An average cost of 11 bits is reported in [27] using delta encoding andHuffman codes.

For the sake of presentation, in the two next steps we
simply denote by r(x) the residuals associated withthese w assignments. The two steps are applied to all
w assignments.
22) compute the squared distance d uj (r(x)), cj,i for eachsubquantizer j and each of its centroids cj,i ;3) compute the squared distance between r(x) and all theindexed vectors of the inverted list. Using the subvectorto-centroid distances computed in the previous step,
this consists in summing up m looked-up values, see
Equation 31;
4) select the K nearest neighbors of x based on theestimated distances. This is implemented efficiently by
maintaining a Maxheap structure of fixed capacity, that
stores the K smallest values seen so far. After eachdistance calculation, the point identifier is added to thestructure only if its distance is below the largest distancein the Maxheap.
Only Step 3 depends on the database size. Compared with
ADC, the additional step of quantizing x to qc (x) consistsin computing k 0 distances between D-dimensional vectors.
Assuming that the inverted lists are balanced, about n × w/k 0
entries have to be parsed. Therefore, the search is significantly
faster than ADC, as shown in the next section.
V. E VALUATION OF NN SEARCH
In this section, we first present the datasets used for theevaluation3 . We then analyze the impact of the parameters
for SDC, ADC and IVFADC. Our approach is compared tothree state-of-the-art methods: spectral hashing [19], Hamming
embedding [20] and FLANN [9]. Finally, we evaluate thecomplexity and speed of our approach.
A. Datasets
We perform our experiments on two datasets, one with local
SIFT descriptors [23] and the other with global color GIST
descriptors [18]. We have three vector subsets per dataset:
3 Both the software and the data used in these experiments are available at
http://www.irisa.fr/texmex/people/jegou/ann.php.

8vector dataset:
descriptor dimensionality dlearning set size
database set size
queries set size

SIFT
128
100,000
1,000,000
10,000

SDC
1

4096
1024

0.8
recall@100

learning, database and query. Both datasets were constructedusing publicly available data and software. For the SIFT
descriptors, the learning set is extracted from Flickr images
and the database and query descriptors are from the INRIA
Holidays images [20]. For GIST, the learning set consistsof the first 100k images extracted from the tiny image setof [16]. The database set is the Holidays image set combinedwith Flickr1M used in [20]. The query vectors are from theHolidays image queries. Table III summarizes the number ofdescriptors extracted for the two datasets.

256

0.6

64

0.4
k*=16

m=1
m=2
m=4
m=8
m=16

0.2

GIST
960
100,000
1,000,991
500

0
0

16

32

64
96
code length (bits)

128

160

TABLE IIIS UMMARY OF THE SIFT AND GIST DATASETS .
ADC
1024

1

B. Memory vs search accuracy: trade-offs
The product quantizer is parametrized by the number ofsubvectors m and the number of quantizers per subvector k ∗ ,producing a code of length m × log2 k ∗ . Figure 6 shows thetrade-off between code length and search quality for our SIFT
descriptor dataset. The quality is measured for recall@100 for
the ADC and SDC estimators, for m ∈ {1, 2, 4, 8, 16} andk ∗ ∈ {24 , 26 , 28 , 210 , 212 }. As for the quantizer distortion in
Figure 1, one can observe that for a fixed number of bits, it
is better to use a small number of subquantizers with many
centroids than to have many subquantizers with few bits. However, the comparison also reveals that MSE underestimates, fora fixed number of bits, the quality obtained for a large numberof subquantizers against using more centroids per quantizer.
As expected, the asymmetric estimator ADC significantly
outperforms SDC. For m=8 we obtain the same accuracy for
ADC and k ∗ =64 as for SDC and k ∗ =256. Given that theefficiency of the two approaches is equivalent, we advocate
not to quantize the query when possible, but only the databaseelements.
Figure 7 evaluates the impact of the parameters for the
IVFADC method introduced in Section IV. For this approach,we have, in addition, to set the codebook size k 0 and thenumber of neighboring cells w visited during the multiple

recall@100

0.8

64

0.6
k*=16

0.4

m=1
m=2
m=4
m=8
m=16

0.2

0
0

16

32

64
96
code length (bits)

128

160

Fig. 6. SDC and ADC estimators evaluated on the SIFT dataset: recall@100
as a function of the memory usage (code length=m × log2 k∗ ) for different
parameters (k∗ =16,64,256,...,4096 and m=1,2,4,8,16). The missing point(m=16,k∗ =4096) gives recall@100=1 for both SDC and ADC.

IVFADC1

0.8
recall@100

The search quality is measured with recall@R, i.e., theproportion of query vectors for which the nearest neighboris ranked in the first R positions. This measure indicatesthe fraction of queries for which the nearest neighbor isretrieved correctly, if a short-list of R vectors is verified
using Euclidean distances. Furthermore, the curve obtained byvarying R corresponds to the distribution function of the ranks,
and the point R=1 corresponds to the “precision” measure usedin [9] to evaluate ANN methods.
In practice, we are often interested in retrieving the K
nearest neighbors (K > 1) and not only the nearest neighbor.
We do not include these measures in the paper, as we observedthat the conclusions for K=1 remain valid for K > 1.

4096

256

0.6
2

0.4

4

16

8

m=1

k’=1024, w=1
k’=1024, w=8
k’=8192, w=8
k’=8192, w=64

0.2

0
0

16

32

64
code length (bits)

96

128

Fig. 7. SIFT dataset: recall@100 for the IVFADC approach as a function of
the memory usage for k∗ =256 and varying values of m = {1, 2, 4, 8, 16},k0 = {1024, 8192} and w = {1, 8, 64}.

9

assignment. We observe that the recall@100 strongly depends
on these parameters, and that increasing the code length isuseless if w is not big enough, as the nearest neighbors which
are not assigned to one of the w centroids associated with thequery are definitely lost.
This approach is significantly more efficient than SDC and
ADC on large datasets, as it only compares the query to a small
fraction of the database vectors. The proportion of the dataset
to visit is roughly linear in w/k 0 . For a fixed proportion, it isworth using higher values of k 0 , as this increases the accuracy,
as shown by comparing, for the tuple (m, w), the parameters(1024, 1) against (8192, 8) and (1024, 8) against (8192, 64).

The product quantizer defined in Section II creates thesubvectors by splitting the input vector according to theorder of the components. However, vectors such as SIFT
and GIST descriptors are structured because they are built
as concatenated orientation histograms. Each histogram iscomputed on grid cells of an image patch. Using a productquantizer, the bins of a histogram may end up in different
quantization groups.
The natural order corresponds to grouping consecutive
components, as proposed in Equation 8. For the SIFT descriptor, this means that histograms of neighboring grid cells arequantized together. GIST descriptors are composed of three320-dimension blocks, one per color channel. The productquantizer splits these blocks into parts.
SIFT
4
8
0.593
0.921
0.501
0.859
0.640
0.905

Discussion: A method that automatically groups the components could further improve the results. This seems particularly
important if we have no prior knowledge about the relationshipbetween the components as in the case of bag-of-features.
A possible solution is the minimum sum-squared residue coclustering [30] algorithm.
D. Comparison with the state of the art

C. Impact of the component groupingmnaturalrandomstructured

the same orientation. On SIFT descriptors, this is a slightly
less efficient structure, probably because the natural ordercorresponds to spatially related components. On GIST, thischoice significantly improves the performance. Therefore, we
use this ordering in the following experiments.

GIST
8
0.338
0.286
0.652

TABLE IVI MPACT OF THE DIMENSION GROUPING ON THE RETRIEVALPERFORMANCE OF ADC ( RECALL @100, k ∗ =256).

To evaluate the influence of the grouping, we modify
the uj operators in Equation 8, and measure the impact of
their construction on the performance of the ADC method.
Table IV shows the effect on the search quality, measured
by recall@100. The analysis is restricted to the parametersk ∗ =256 and m ∈ {4, 8}.
Overall, the choice of the components appears to have asignificant impact of the results. Using a random order insteadof the natural order leads to poor results. This is true even for
GIST, for which the natural order is somewhat arbitrary.
The “structured” order consists in grouping together dimensions that are related. For the m = 4 SIFT quantizer, thismeans that the 4×4 patch cells that make up the descriptor [23]are grouped into 4 2 × 2 blocks. For the other two, it groups
together dimensions that have have the same index modulo 8.
The orientation histograms of SIFT and most of GIST’s have 8bins, so this ordering quantizes together bins corresponding toComparison with Hamming embedding methods: We compare our approach to spectral hashing (SH) [19], which maps
vectors to binary signatures. The search consists in comparing
the Hamming distances between the database signatures andthe query vector signature. This approach was shown tooutperform the restricted Boltzmann machine of [17]. We haveused the publicly available code. We also compare to theHamming embedding (HE) method of [20], which also maps
vectors to binary signatures. Similar to IVFADC, HE uses
an inverted file, which avoids comparing to all the databaseelements.
Figures 8 and 9 show, respectively for the SIFT and the
GIST datasets, the rank repartition of the nearest neighbors
when using a signature of size 64 bits. For our productquantizer we have used m = 8 and k ∗ = 256, whichgive similar results in terms of run time. All our approaches
significantly outperform spectral hashing4 on the two datasets.
To achieve the same recall as spectral hashing, ADC returns
an order of magnitude less vectors.
The best results are obtained by IVFADC, which for low
ranks provides an improvement over ADC, and significantly
outperforms spectral hashing. This strategy avoids the exhaustive search and is therefore much faster, as discussed in thenext subsection. This partial scan explains why the IVFADC
and HE curves stop at some point, as only a fraction of thedatabase vectors are ranked. Comparing these two approaches,
HE is significantly outperformed by IVFADC. The results of
HE are similar to spectral hashing, but HE is more efficient.
Comparison with FLANN: The approximate nearestneighbor search technique of Muja & Lowe [9] is basedon hierarchical structures (KD-trees and hierarchical k-means
trees). The software package FLANN automatically selects thebest algorithm and parameters for a given dataset. In contrast
with our method and spectral hashing, all vectors need toremain in RAM as the method includes a re-ranking stage
that computes the real distances for the candidate nearest
neighbors.
The evaluation is performed on the SIFT dataset by measuring the 1-recall@1, that is, the average proportion of true NNs
4 In defense of [17], [19], which can be learned for arbitrary distancemeasures, our approach is adapted to the Euclidean distance only.

10

SIFT, 64-bit codes

100

FLANN
IVFADC, R=10
IVFADC, R=100
IVFADC, R=1000

1

16/128search time (s)recall@R

0.8

0.6

0.4

SDC
ADC
IVFADC w=1
IVFADC w=16
HE w=1
HE w=16
spectral hashing

0.2

16/102410
4/128
4/1024
4/1024

0

1
1

10

100

1k
R

10k

100k

1M

Fig. 8. SIFT dataset: recall@R for varying values of R. Comparison of
the different approaches SDC, ADC, IVFADC, spectral hashing [19] andHE [20]. Parameters: m=8 and k∗ =256 for SDC/ADC, k0 =1024 for HE [20]and IVFADC.
GIST, 64-bit codes

0.8

0.6

0.4
SDC
ADC
IVFADC w=1
IVFADC w=8
IVFADC w=64
spectral hashing

0.2

0
1

10

100

1k
R

10k

100k

0.6

0.65

0.7

0.75
0.8
0.85
1-recall at 1

0.9

0.95

1

Fig. 10. IVFADC vs FLANN: trade-offs between search quality (1-recall@1)
and search time. The IVFADC method is parametrized by the shortlist size R
used for re-ranking the vector with the L2 distance, and the two parameters
w and k0 of the inverted file, which correspond to the number of assignments
and to the number of coarse centroids.

probably not optimal for all operating points.
Figure 10 shows that the results obtained by IVFADC are
better than those of FLANN for a large range of operating
points. Moreover, our method has a much smaller memory
footprint than FLANN: the indexing structure occupies lessthan 25 MB, while FLANN requires more than 250 MB ofRAM. Note, however, that both are negligible compared to thememory occupied by the vectors in the case of large datasets.
On such a scale, the re-ranking stage is not feasible and only
memory-aware approaches (HE, SH and our methods) can be
used.

1

recall@R

16/1024
16/1281M

Fig. 9. GIST dataset: recall@R for varying values of R. Comparison of
the different approaches SDC, ADC, IVFADC and spectral hashing [19].
Parameters: m=8, k∗ =256 for SDC/ADC and k0 = 1 024 for IVFADC.

ranked first in the returned vectors. This measure is referred
to as precision in [9].
For the sake of comparison with FLANN, we added averification stage to our IVFADC method: IVFADC queries
return a shortlist of R candidate nearest neighbors using the
distance estimation. The vectors in the shortlist are re-orderedusing the real distance, as done in [7], [9], and the closest
one is returned. Note that, in this experimental setup, all thevectors are stored in main memory. This requirement seriously
limits the scale on which re-ordering can be used.
The IVFADC and FLANN methods are both evaluated at
different operating points with respect to precision and searchtime. For FLANN, the different operating points are obtainedwith parameters generated automatically for various target
precisions. For IVFADC, they are obtained by varying thenumber k 0 of coarse centroids, the number w of assignments
and the short-list size R. The product quantizer is generatedusing k ∗ =256 and m=8, i.e., 64-bit codes. This choice is

E. Complexity and speed
Table V reports the search time of our methods. For
reference, we report the results obtained with the spectral
hashing algorithm of [19] on the same dataset and machine(using only one core). Since we use a separate learning set,
we use the out-of-sample evaluation of this algorithm. Notethat for SH we have re-implemented the Hamming distancecomputation in C in order to have the approaches similarly
optimized. The algorithms SDC, ADC and SH have similar
run times. IVFADC significantly improves the performance by
avoiding an exhaustive search. Higher values of k 0 yield higher
search efficiencies for large datasets, as the search benefits
from parsing a smaller fraction of the memory. However, forsmall datasets, the complexity of the coarse quantizer may be
the bottleneck if k 0 × D > n/k 0 when using a exhaustive
assignment for qc . In that case the ADC variant may be
preferred. For large datasets and using an efficient assignment
strategy for the coarse quantizer, higher values of k 0 generally
lead to better efficiency, as first shown in [15]. In this work, theauthors propose a hierarchical quantizer to efficiently assign
descriptors to one million centroids.
F. Large-scale experiments
To evaluate the search efficiency of the product quantizer
method on larger datasets we extracted about 2 billion SIFT

11

methodSDC
ADC
IVFADCparameters

k 0 = 1 024, w=1
k 0 = 1 024, w=8
k 0 = 1 024, w=64
k 0 = 8 192, w=1
k 0 = 8 192, w=8
k 0 = 8 192, w=64

SH

searchtime (ms)
16.8
17.2
1.5
8.8
65.9
3.8
10.2
65.3
22.7

average number ofcode comparisons
1 000 991
1 000 991
1 947
27 818101 158361
2 709
19 101
1 000 991

recall@100
0.446
0.652
0.308
0.682
0.744
0.240
0.516
0.610
0.132

TABLE V
GIST DATASET (500 QUERIES ): SEARCH TIMINGS FOR 64- BIT CODES AND DIFFERENT METHODS . m=8 AND k∗ =256 FOR SDC, ADC AND IVFADC.

3.5

IVF+HE 64 bits
IVFADC 64 bits (m=8)
IVF+HE 32 bits
IVFADC 32 bits (m=4)

0.75
0.7

2.5

0.65

2

mAP

search time (ms/vector)

3

0.8

HE
IVFADC

1.5

0.6
0.55

1

0.5

0.5
0

0.45
10M

100M
database size

1G

Fig. 11. Search times for SIFT descriptors in datasets of increasing sizes,
with two search methods. Both use the same 20 000-word codebook, w = 1,
and 64-bit signatures.

descriptors from one million images. Search is performed
with 30 000 query descriptors from ten images. We compared
the IVFADC and HE methods with similar parameters. In
particular, the amount of memory that is scanned for eachmethod and the cost of the coarse quantization are the same.
The query times per descriptor are shown on Figure 11.
The cost of the extra quantization step required by IVFADC
appears clearly for small database sizes. For larger scales,
the distance computation with the database vectors become
preponderant. The processing that is applied to each element of
the inverted lists is approximately as expensive in both cases.
For HE, it is a Hamming distance computation, implemented
as 8 table lookups. For IVFADC it is a distance computationthat is also performed by 8 table lookups. Interestingly, thefloating point operations involved in IVFADC are not much
more expensive than the simple binary operations of HE.
G. Image search
We have evaluated our method within an image searchsystem based on local descriptors. For this evaluation, we

0.4
1k

10k
100k
database size (number of images)

1M

Fig. 12. Comparison of IVFADC and the Hamming Embedding methodof [20]. mAP for the Holidays dataset as function of the number of distractor
images (up to 1 million).

compare our method with the HE method of [20] on the INRIA
Holidays dataset, using the pre-processed set of descriptors
available online. The comparison is focused on large scaleindexing, i.e., we do not consider the impact of a postverification step [23], [31] or geometrical information [20].
Figure 12 shows the search performance in terms of mean
average precision as a function of the size of the dataset. We
have used the same coarse quantizer (k 0 =20,000) and a single
assignment strategy (w=1) for both the approaches, and fixedk ∗ =256 for IVFADC. For a given number of bits (32 or 64),
we have selected the best choice of the Hamming thresholdfor HE. Similarly, we have adjusted the number of nearest
neighbors to be retrieved for IVFADC.
One can observe that the gain obtained by IVFADC issignificant. For example, for one million distractors, the mAPof 0.497 reported in [20] with 64-bit signatures is improved
to 0.517.

12

VI. C ONCLUSION
We have introduced product quantization for approximate
nearest neighbor search. Our compact coding scheme provides an accurate approximation of the Euclidean distance.
Moreover, it is combined with an inverted file system toavoid exhaustive search, resulting in high efficiency. Our
approach significantly outperforms the state of the art in termsof the trade-off between search quality and memory usage.
Experimental results for SIFT and GIST image descriptors are
excellent and show that grouping the components based on
our prior knowledge of the descriptor design further improves
the results. The scalability of our approach is validated on adataset of two billion vectors.
R EFERENCES[1] H. Jégou, M. Douze, and C. Schmid, “Searching with quantization:
approximate nearest neighbor search using short codes and distanceestimators,” Tech. Rep. RR-7020, INRIA, August 2009.
[2] H. Sandhawalia and H. Jégou, “Searching with expectations,” in IEEE
International Conference on Acoustics Speech and Signal Processing,
Signal Processing, pp. 1242–1245, IEEE, March 2010.
[3] K. Beyer, J. Goldstein, R. Ramakrishnan, and U. Shaft, “When is”nearest neighbor” meaningful?,” in Proceedings of the International
Conference on Database Theory, pp. 217–235, August 1999.
[4] C. Böhm, S. Berchtold, and D. Keim, “Searching in high-dimensional
spaces: Index structures for improving the performance of multimedia
databases,” ACM Computing Surveys, vol. 33, pp. 322–373, October
2001.
[5] J. Friedman, J. L. Bentley, and R. A. Finkel, “An algorithm forfinding best matches in logarithmic expected time,” ACM Transaction
on Mathematical Software, vol. 3, no. 3, pp. 209–226, 1977.
[6] R. Weber, H.-J. Schek, and S. Blott, “A quantitative analysis andperformance study for similarity-search methods in high-dimensional
spaces,” in Proceedings of the International Conference on Very Large
DataBases, pp. 194–205, 1998.
[7] M. Datar, N. Immorlica, P. Indyk, and V. Mirrokni, “Locality-sensitive
hashing scheme based on p-stable distributions,” in Proceedings of the
Symposium on Computational Geometry, pp. 253–262, 2004.
[8] A. Gionis, P. Indyk, and R. Motwani, “Similarity search in highdimension via hashing,” in Proceedings of the International Conference
on Very Large DataBases, pp. 518–529, 1999.
[9] M. Muja and D. G. Lowe, “Fast approximate nearest neighbors with
automatic algorithm configuration,” in Proceedings of the International
Conference on Computer Vision Theory and Applications, 2009.
[10] B. Kulis and K. Grauman, “Kernelized locality-sensitive hashing for
scalable image search,” in Proceedings of the International Conference
on Computer Vision, October 2009.
[11] G. Shakhnarovich, T. Darrell, and P. Indyk, Nearest-Neighbor Methods
in Learning and Vision: Theory and Practice, ch. 3. MIT Press, March
2006.
[12] Y. Ke, R. Sukthankar, and L. Huston, “Efficient near-duplicate detection
and sub-image retrieval,” in ACM International conference on Multimedia, pp. 869–876, 2004.

[13] B. Matei, Y. Shan, H. Sawhney, Y. Tan, R. Kumar, D. Huber, and
M. Hebert, “Rapid object indexing using locality sensitive hashing andjoint 3D-signature space estimation,” IEEE Transactions on PatternAnalysis and Machine Intelligence, vol. 28, pp. 1111 – 1126, July 2006.
[14] C. Silpa-Anan and R. Hartley, “Optimized kd-trees for fast image descriptor matching,” in Proceedings of the IEEE Conference on Computer
Vision and Pattern Recognition, 2008.
[15] D. Nistér and H. Stewénius, “Scalable recognition with a vocabularytree,” in Proceedings of the IEEE Conference on Computer Vision and
Pattern Recognition, pp. 2161–2168, 2006.
[16] A. Torralba, R. Fergus, and W. T. Freeman, “80 million tiny images: alarge database for non-parametric object and scene recognition,” IEEE
Transactions on Pattern Analysis and Machine Intelligence, vol. 30,
pp. 1958–1970, November 2008.
[17] A. Torralba, R. Fergus, and Y. Weiss, “Small codes and large databases
for recognition,” in Proceedings of the IEEE Conference on Computer
Vision and Pattern Recognition, 2008.
[18] A. Oliva and A. Torralba, “Modeling the shape of the scene: a holistic
representation of the spatial envelope,” International Journal of Computer Vision, vol. 42, no. 3, pp. 145–175, 2001.
[19] Y. Weiss, A. Torralba, and R. Fergus, “Spectral hashing,” in Advancesin Neural Information Processing Systems, 2008.
[20] H. Jégou, M. Douze, and C. Schmid, “Hamming embedding and weakgeometric consistency for large scale image search,” in Proceedings ofthe European Conference on Computer Vision, October 2008.
[21] M. Douze, H. Jégou, H. Singh, L. Amsaleg, and C. Schmid, “Evaluationof GIST descriptors for web-scale image search,” in Proceedings of theInternational Conference on Image and Video Retrieval, 2009.
[22] J. Philbin, O. Chum, M. Isard, J. Sivic, and A. Zisserman, “Lost inquantization: Improving particular object retrieval in large scale image
databases,” in Proceedings of the IEEE Conference on Computer Vision
and Pattern Recognition, 2008.
[23] D. Lowe, “Distinctive image features from scale-invariant keypoints,”
International Journal of Computer Vision, vol. 60, no. 2, pp. 91–110,
2004.
[24] R. M. Gray and D. L. Neuhoff, “Quantization,” IEEE Transactions on
Information Theory, vol. 44, pp. 2325–2384, Oct. 1998.
[25] D. E. Knuth, The Art of Computer Programming, Sorting and Searching,vol. 3. Addison Wesley, 2 ed., 1998.
[26] J. Sivic and A. Zisserman, “Video Google: A text retrieval approachto object matching in videos,” in Proceedings of the International
Conference on Computer Vision, pp. 1470–1477, 2003.
[27] M. Perdoch, O. Chum, and J. Matas, “Efficient representation of local
geometry for large scale object retrieval,” in Proceedings of the IEEE
Conference on Computer Vision and Pattern Recognition, June 2009.
[28] H. Jégou, M. Douze, and C. Schmid, “Packing bag-of-features,” in Proceedings of the International Conference on Computer Vision, September
2009.
[29] H. Jégou, H. Harzallah, and C. Schmid, “A contextual dissimilarity
measure for accurate and efficient image search,” in Proceedings of theIEEE Conference on Computer Vision and Pattern Recognition, 2007.
[30] H. Cho, I. Dhillon, Y. Guan, and S. Sra, “Minimum sum-squared
residue co-clustering of gene expression data,” in SIAM International
Conference on Data Mining, pp. 114–125, April 2004.
[31] J. Philbin, O. Chum, M. Isard, J. Sivic, and A. Zisserman, “Object retrieval with large vocabularies and fast spatial matching,” in Proceedingsof the IEEE Conference on Computer Vision and Pattern Recognition,2007.