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

[Link]
Submitted on 22 May 2013

HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
 1

Product quantization for nearest neighbor search

Hervé Jégou, Matthijs Douze, Cordelia Schmid
INRIA

Abstract— This paper introduces a product quantization based with limited assumptions. It has been successfully used for
approach for approximate nearest neighbor search. The idea local descriptors [12] and 3D object indexing [13], [11].
is to decomposes the space into a Cartesian product of low However, for real data, LSH is outperformed by heuristic
dimensional subspaces and to quantize each subspace separately.
A vector is represented by a short code composed of its subspace methods, which exploit the distribution of the vectors. These
quantization indices. The Euclidean distance between two vectors methods include randomized KD-trees [14] and hierarchical
can be efficiently estimated from their codes. An asymmetric k-means [15], both of which are implemented in the FLANN
version increases precision, as it computes the approximate selection algorithm [9].
distance between a vector and a code. ANN algorithms are typically compared based on the trade-
Experimental results show that our approach searches for
nearest neighbors efficiently, in particular in combination with an off between search quality and efficiency. However, this trade-
inverted file system. Results for SIFT and GIST image descriptors off does not take into account the memory requirements of
show excellent search accuracy outperforming three state-of-the- the indexing structure. In the case of E2LSH, the memory
art approaches. The scalability of our approach is validated on usage may even be higher than that of the original vectors.
a dataset of two billion vectors. Moreover, both E2LSH and FLANN need to perform a final
Index Terms— High-dimensional indexing, image indexing, re-ranking step based on exact L2 distances, which requires the
very large databases, approximate search. indexed vectors to be stored in main memory if access speed
is important. This constraint seriously limits the number of
I. I NTRODUCTION vectors that can be handled by these algorithms. Only recently,
researchers came up with methods limiting the memory usage.
OMPUTING Euclidean distances between high dimen-
C sional vectors is a fundamental requirement in many
applications. It is used, in particular, for nearest neighbor
This is a key criterion for problems involving large amounts
of data [16], i.e., in large-scale scene recognition [17], where
millions to billions of images have to be indexed. In [17],
(NN) search. Nearest neighbor search is inherently expensive Torralba et al. represent an image by a single global GIST
due to the curse of dimensionality [3], [4]. Focusing on the descriptor [18] which is mapped to a short binary code.
D-dimensional Euclidean space RD , the problem is to find When no supervision is used, this mapping is learned such
the element NN(x), in a finite set Y ⊂ RD of n vectors, that the neighborhood in the embedded space defined by the
minimizing the distance to the query vector x ∈ RD : Hamming distance reflects the neighborhood in the Euclidean
NN(x) = arg min d(x, y). (1) space of the original features. The search of the Euclidean
y∈Y nearest neighbors is then approximated by the search of the
Several multi-dimensional indexing methods, such as the nearest neighbors in terms of Hamming distances between
popular KD-tree [5] or other branch and bound techniques, codes. In [19], spectral hashing (SH) is shown to outper-
have been proposed to reduce the search time. However, for form the binary codes generated by the restricted Boltzmann
high dimensions it turns out [6] that such approaches are machine [17], boosting and LSH. Similarly, the Hamming
not more efficient than the brute-force exhaustive distance embedding method of Jegou et al. [20], [21] uses a binary
calculation, whose complexity is O(nD). signature to refine quantized SIFT or GIST descriptors in a
There is a large body of literature [7], [8], [9] on algorithms bag-of-features image search framework.
that overcome this issue by performing approximate nearest In this paper, we construct short codes using quantization.
neighbor (ANN) search. The key idea shared by these algo- The goal is to estimate distances using vector-to-centroid
rithms is to find the NN with high probability “only”, instead distances, i.e., the query vector is not quantized, codes are
of probability 1. Most of the effort has been devoted to the assigned to the database vectors only. This reduces the quan-
Euclidean distance, though recent generalizations have been tization noise and subsequently improves the search quality.
proposed for other metrics [10]. In this paper, we consider the To obtain precise distances, the quantization error must be
Euclidean distance, which is relevant for many applications. limited. Therefore, the total number k of centroids should
In this case, one of the most popular ANN algorithms is be sufficiently large, e.g., k = 264 for 64-bit codes. This
the Euclidean Locality-Sensitive Hashing (E2LSH) [7], [11], raises several issues on how to learn the codebook and assign
which provides theoretical guarantees on the search quality a vector. First, the number of samples required to learn the
quantizer is huge, i.e., several times k. Second, the complexity
This work was partly realized as part of the Quaero Programme, funded of the algorithm itself is prohibitive. Finally, the amount of
by OSEO, French State agency for innovation. It was originally published as computer memory available on Earth is not sufficient to store
a technical report [1] in August 2009. It is also related to the work [2] on
source coding for nearest neighbor search. This particular version is a post- the floating point values representing the centroids.
print correcting mistakes reported by Jan Stria and Ondrej Chum (CVUT). The hierarchical k-means see (HKM) improves the effi-
 2

ciency of the learning stage and of the corresponding assign- reduce the cardinality of the representation space, in particular
ment procedure [15]. However, the aforementioned limitations when the input data is real-valued. Formally, a quantizer is a
still apply, in particular with respect to memory usage and size function q mapping a D-dimensional vector x ∈ RD to a
of the learning set. Another possibility are scalar quantizers, vector q(x) ∈ C = {ci ; i ∈ I}, where the index set I is from
but they offer poor quantization error properties in terms now on assumed to be finite: I = 0 . . . k −1. The reproduction
of the trade-off between memory and reconstruction error. values ci are called centroids. The set of reproduction values
Lattice quantizers offer better quantization properties for uni- C is the codebook of size k.
form vector distributions, but this condition is rarely satisfied The set Vi of vectors mapped to a given index i is referred
by real world vectors. In practice, these quantizers perform to as a (Voronoi) cell, and defined as
significantly worse than k-means in indexing tasks [22]. In
this paper, we focus on product quantizers. To our knowledge, Vi , {x ∈ RD : q(x) = ci }. (2)
such a semi-structured quantizer has never been considered in
The k cells of a quantizer form a partition of RD . By
any nearest neighbor search method.
definition, all the vectors lying in the same cell Vi are recon-
The advantages of our method are twofold. First, the number
structed by the same centroid ci . The quality of a quantizer
of possible distances is significantly higher than for competing
is usually measured by the mean squared error between the
Hamming embedding methods [20], [17], [19], as the Ham-
input vector x and its reproduction value q(x):
ming space used in these techniques allows for a few distinct Z
distances only. Second, as a byproduct of the method, we
2
MSE(q) = EX d(q(x), x)2 = p(x) d q(x), x dx, (3)
 
get an estimation of the expected squared distance, which is
required for ε-radius search or for using Lowe’s distance ratio
where d(x, y) = ||x − y|| is the Euclidean distance between
criterion [23]. The motivation of using the Hamming space
x and y, and where p(x) is the probability distribution func-
in [20], [17], [19] is to compute distances efficiently. Note,
tion corresponding the random variable X. For an arbitrary
however, that one of the fastest ways to compute Hamming
probability distribution function, Equation 3 is numerically
distances consists in using table lookups. Our method uses
computed using Monte-Carlo sampling, as the average of
a similar number of table lookups, resulting in comparable
||q(x) − x||2 on a large set of samples.
efficiency.
In order for the quantizer to be optimal, it has to satisfy two
An exhaustive comparison of the query vector with all codes
properties known as the Lloyd optimality conditions. First, a
is prohibitive for very large datasets. We, therefore, introduce
vector x must be quantized to its nearest codebook centroid,
a modified inverted file structure to rapidly access the most
in terms of the Euclidean distance:
relevant vectors. A coarse quantizer is used to implement this
inverted file structure, where vectors corresponding to a cluster q(x) = arg min d(x, ci ). (4)
(index) are stored in the associated list. The vectors in the ci ∈C

list are represented by short codes, computed by our product As a result, the cells are delimited by hyperplanes. The
quantizer, which is used here to encode the residual vector second Lloyd condition is that the reconstruction value must
with respect to the cluster center. be the expectation of the vectors lying in the Voronoi cell:
The interest of our method is validated on two kinds Z
of vectors, namely local SIFT [23] and global GIST [18]  
ci = EX x|i = p(x|x ∈ Vi ) x dx. (5)
descriptors. A comparison with the state of the art shows Vi
that our approach outperforms existing techniques, in par- The Lloyd quantizer, which corresponds to the k-means
ticular spectral hashing [19], Hamming embedding [20] and clustering algorithm, finds a near-optimal codebook by iter-
FLANN [9]. atively assigning the vectors of a training set to centroids and
Our paper is organized as follows. Section II introduces re-estimating these centroids from the assigned vectors. In the
the notations for quantization as well as the product quantizer following, we assume that the two Lloyd conditions hold, as
used by our method. Section III presents our approach for we learn the quantizer using k-means. Note, however, that k-
NN search and Section IV introduces the structure used to means does only find a local optimum in terms of quantization
avoid exhaustive search. An evaluation of the parameters of error.
our approach and a comparison with the state of the art is
Another quantity that will be used in the following is the
given in Section V.
mean squared distortion ξ(q, ci ) obtained when reconstructing
II. BACKGROUND : QUANTIZATION , PRODUCT QUANTIZER 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
A large body of literature is available on vector quantiza-
to the centroid ci , it is computed as
tion, see [24] for a survey. In this section, we restrict our Z
presentation to the notations and concepts used in the rest of 1
2
ξ(q, ci ) = d x, q(x) p(x) dx. (6)
the paper. pi Vi

A. Vector quantization Note that the MSE can be obtained from these quantities as
X
Quantization is a destructive process which has been ex- MSE(q) = pi ξ(q, ci ). (7)
tensively studied in information theory [24]. Its purpose is to i∈I
 3

The memory cost of storing the index value, without any memory usage assignment complexity
k-means kD kD
further processing (entropy coding), is dlog2 ke bits. Therefore, bf
HKM (k − 1) D lD
it is convenient to use a power of two for k, as the code b −1
f
product k-means m k∗ D∗ = k1/m D m k∗ D∗ = k1/m D
produced by the quantizer is stored in a binary memory.
TABLE I
B. Product quantizers M EMORY USAGE OF THE CODEBOOK AND ASSIGNMENT COMPLEXITY FOR
Let us consider a 128-dimensional vector, for example the DIFFERENT QUANTIZERS . HKM IS PARAMETRIZED BY TREE HEIGHT l

SIFT descriptor [23]. A quantizer producing 64-bits codes, AND THE BRANCHING FACTOR bf .

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 complexity 0.3
m=1
of learning the quantizer are several times k. It is even m=2
0.25 m=4
impossible to store the D×k floating point values representing m=8
the k centroids. m=16

square distortion D(q)

Product quantization is an efficient solution to address these 0.2 k*=16

issues. It is a common technique in source coding, which

allows to choose the number of components to be quantized 0.15
64
jointly (for instance, groups of 24 components can be quan-
tized using the powerful Leech lattice). The input vector x is 0.1 256

split into m distinct subvectors uj , 1 ≤ j ≤ m of dimension 1024

D∗ = D/m, where D is a multiple of m. The subvectors 0.05
are quantized separately using m distinct quantizers. A given
vector x is therefore mapped as follows: 0
0 16 32 64 96 128 160


x , ..., x ∗ , ..., xD−D∗ +1 , ..., xD → q1 u1 (x)), ..., qm (um (x) , code length (bits)
| 1 {z D } | {z }
u1 (x) um (x) Fig. 1. SIFT: quantization error associated with the parameters m and k∗ .
(8)
where qj is a low-complexity quantizer associated with the
j th subvector. With the subquantizer qj we associate the index m D∗ k ∗ = k ∗ D floating points values. Quantizing an element
set Ij , the codebook Cj and the corresponding reproduction requires k ∗ D floating point operations. Table I summarizes
values cj,i . the resource requirements associated with k-means, HKM and
A reproduction value of the product quantizer is identified product k-means. The product quantizer is clearly the the only
by an element of the product index set I = I1 × . . . × Im . one that can be indexed in memory for large values of k.
The codebook is therefore defined as the Cartesian product In order to provide good quantization properties when
C = C1 × . . . × C m , (9) choosing a constant value of k ∗ , each subvector should have,
on average, a comparable energy. One way to ensure this prop-
and a centroid of this set is the concatenation of centroids erty is to multiply the vector by a random orthogonal matrix
of the m subquantizers. From now on, we assume that all prior to quantization. However, for most vector types this is not
subquantizers have the same finite number k ∗ of reproduction required and not recommended, as consecutive components
values. In that case, the total number of centroids is given by are often correlated by construction and are better quantized
together with the same subquantizer. As the subspaces are
k = (k ∗ )m . (10)
orthogonal, the squared distortion associated with the product
Note that in the extremal case where m = D, the com- quantizer is X
ponents of a vector x are all quantized separately. Then the MSE(q) = MSE(qj ), (11)
product quantizer turns out to be a scalar quantizer, where the j
quantization function associated with each component may be
where MSE(qj ) is the distortion associated with quantizer qj .
different.
Figure 1 shows the MSE as a function of the code length
The strength of a product quantizer is to produce a large for different (m,k ∗ ) tuples, where the code length is l =
set of centroids from several small sets of centroids: those m log2 k ∗ , if k ∗ is a power of two. The curves are obtained
associated with the subquantizers. When learning the subquan- for a set of 128-dimensional SIFT descriptors, see section V
tizers using Lloyd’s algorithm, a limited number of vectors is for details. One can observe that for a fixed number of bits,
used, but the codebook is, to some extent, still adapted to the it is better to use a small number of subquantizers with many
data distribution to represent. The complexity of learning the centroids than having many subquantizers with few bits. At
quantizer is m times the complexity of performing k-means the extreme when m = 1, the product quantizer becomes a
clustering with k ∗ centroids of dimension D∗ . regular k-means codebook.
Storing the codebook C explicitly is not efficient. Instead, High values of k ∗ increase the computational cost of the
we store the m × k ∗ centroids of all the subquantizers, i.e., quantizer, as shown by Table I. They also increase the memory
 4

SDC ADC
encoding x k∗ D 0
k∗ D


compute d uj (x), cj,i 0
ˆ y) or d(x,
for y ∈ Y, compute d(x, ˜ y) nm nm
find the k smallest distances n + k log k log log n

TABLE II
A LGORITHM AND COMPUTATIONAL COSTS ASSOCIATED WITH SEARCHING
THE k NEAREST NEIGHBORS USING THE PRODUCT QUANTIZER FOR
SYMMETRIC AND ASYMMETRIC DISTANCE COMPUTATIONS (SDC, ADC).
symmetric case asymmetric case
Fig. 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 is vector y is represented by q(y), but the query x is not encoded.
on average bounded by the quantization error. ˜

d(x, y) is approximated by the distance d(x, y) ,
The distance
d x, q(y) , which is computed using the decomposition
usage of storing the centroids (k ∗ × D floating point values), sX

2
which further reduces the efficiency if the centroid look-up ˜ y) = d x, q(y) =


d(x, d uj (x), qj (uj (y)) , (13)
table does no longer fit in cache memory. In the case where j
m = 1, we can not afford using more than 16 bits to keep
2
this cost tractable. Using k ∗ = 256 and m = 8 is often a where the squared distances d uj (x), cj,i : j = 1 . . . m, i =
reasonable choice. 1 . . . k ∗ , are computed prior to the search.
For nearest neighbors search, we do not compute the square
III. S EARCHING WITH QUANTIZATION roots in practice: the square root function is monotonically
increasing and the squared distances produces the same vector
Nearest neighbor search depends on the distances between
ranking.
the query vector and the database vectors, or equivalently
the squared distances. The method introduced in this section Table II summarizes the complexity of the different steps
compares the vectors based on their quantization indices, involved in searching the k nearest neighbors of a vector x in a
in the spirit of source coding techniques. We first explain dataset Y of n = |Y| vectors. One can see that SDC and ADC
how the product quantizer properties are used to compute have the same query preparation cost, which does not depend
the distances. Then we provide a statistical bound on the on the dataset size n. When n is large (n > k ∗ D∗ ), the most
distance estimation error, and propose a refined estimator for consuming operations are the summations in Equations 12
the squared Euclidean distance. and 13. The complexity given in this table for searching the
k smallest elements is the average complexity for n  k and
A. Computing distances using quantized codes when the elements are arbitrarily ordered ([25], Section 5.3.3,
Equation 17).
Let us consider the query vector x and a database vector
The only advantage of SDC over ADC is to limit the
y. We propose two methods to compute an approximate
memory usage associated with the queries, as the query vector
Euclidean distance between these vectors, a symmetric and
is defined by a code. This is most cases not relevant and one
a asymmetric one. See Figure 2 for an illustration.
should then use the asymmetric version, which obtains a lower
distance distortion for a similar complexity. We will focus on
Symmetric distance computation (SDC): both the vectors x
ADC in the rest of this section.
and y are represented by their respective centroids q(x) and
q(y). The distance d(x,

y) is approximated by the distance
ˆ y) , d q(x), q(y) which is efficiently obtained using a B. Analysis of the distance error
d(x,
product quantizer In this subsection, we analyze the distance error when using
sX ˜ y) instead of d(x, y). This analysis does not depend
d(x,

2
ˆ


d(x, y) = d q(x), q(y) = d qj (x), qj (y) , (12) on the use of a product quantizer and is valid for any
j quantizer satisfying Lloyd’s optimality conditions defined by

2 Equations 4 and 5 in Section II.
where the distance d cj,i , cj,i0 is read from a look-up In the spirit of the mean squared error criterion used for
table associated with the j th subquantizer. Each look-up table reconstruction, the distance distortion is measured by the mean
contains all the squared distances between pairs of centroids squared distance error (MSDE) on the distances:
(i, i0 ) of the subquantizer, or (k ∗ )2 squared distances1 . ZZ
˜ y) 2 p(x)dx p(y)dy. (14)


MSDE(q) , d(x, y) − d(x,
Asymmetric distance computation (ADC): the database
The triangular inequality gives
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 d x, q(y) − d y, q(y) ≤ d(x, y) ≤ d x, q(y) + d y, q(y) ,
are zeros. (15)
 5

between x, which is fully known in our asymmetric distance

computation method, and a random variable Y , subject to
1.2 q(Y ) = q(y) = ci , which represents all the hypothesis on
y knowing its quantization index.
estimated distance

ẽ(x, y) , EY (x − Y )2 |q(Y ) = ci
 
1 (19)
Z
= (x − y)2 p(y|i) dy, (20)
0.8 Vi
Z
1
= (x − ci + ci − y)2 p(y) dy. (21)
0.6 pi Vi
symmetric Developing the squared expression and observing, using
asymmetric Lloyd’s condition of Equation 5, that
0.4
0.4 0.6 0.8 1 1.2
Z
true distances (y − ci ) p(y) dy = 0, (22)
Vi

Equation 21 simplifies to
Fig. 3. Typical query of a SIFT vector in a set of 1000 vectors: comparison Z
of the distance d(x, y) obtained with the SDC and ADC estimators. We have
2
2

used m = 8 and k∗ = 256, i.e., 64-bit code vectors. Best viewed in color. ẽ(x, y) = x − q(y) + q(y) − y p y|y ∈ Vi dy
Vi
(23)
˜ y)2 + ξ q, q(y)


and, equivalently, = d(x, (24)
 2
2

d(x, y) − d(x, q(y)) ≤ d y, q(y) . (16) where we recognize the distortion ξ q, q(y) associated with
the reconstruction of y by its reproduction value.
Combining this inequality with Equation 14, we obtain Using the product quantizer and Equation 24, the computa-
tion of the expected squared distance between a vector x and
Z Z 

2
MSDE(q) ≤ p(x) d y, q(y) p(y) dy dx (17) the vector
y, for which we only know the quantization indices
qj uj (y) , consists in correcting Equation 13 as
≤ MSE(q). (18)
X
˜ y)2 +
ẽ(x, y) = d(x, ξj (y) (25)
where MSE(q) is the mean squared error associated with
j
quantizer q (Equation 3). This inequality, which holds for
any quantizer, shows that the distance error of our method is where the correcting term, i.e., the average distortion
statistically bounded by the MSE associated with the quantizer.


ξj (y) , ξ qj , qj uj (y) (26)
For the symmetric version, a similar derivation shows that the
error is statistically bounded by 2×MSE(q). It is, therefore, th
associated with quantizing uj (y) to qj (y) using the j sub-
worth minimizing the quantization error, as this criterion quantizer, is learned and stored in a look-up table for all
provides a statistical upper bound on the distance error. If an indexes of Ij .
exact distance calculation is performed on the highest ranked Performing a similar derivation for the symmetric version,
vectors, as done in LSH [7], the quantization error can be i.e., when both x and y are encoded using the product
used (instead of selecting an arbitrary set of elements) as a quantizer, we obtain the following corrected version of the
criterion to dynamically select the set of vectors on which the symmetric squared distance estimator:
post-processing should be applied. X X
ˆ y)2 +
ê(x, y) = d(x, ξj (x) + ξj 0 (y). (27)
j j0
C. Estimator of the squared distance
As shown later in this subsection, using the estimations dˆ Discussion: Figure 4 illustrates the probability distribution
or d˜ leads to underestimate, on average, the distance between function of the difference between the true distance and the
descriptors. Figure 3 shows the distances obtained when ones estimated by Equations 13 and 25. It has been measured
querying a SIFT descriptor in a dataset of 1000 SIFT vectors. on a large set of SIFT descriptors. The bias of the distance
It compares the true distance against the estimates computed estimation by Equation 13 is significantly reduced in the
with Equations 12 and 13. One can clearly see the bias on these corrected version. However, we observe that correcting the
distance estimators. Unsurprisingly, the symmetric version is bias leads, in this case, to a higher variance of the estimator,
more sensitive to this bias. which is a common phenomenon in statistics. Moreover, for
Hereafter, we compute the expectation of the squared dis- the nearest neighbors, the correcting term is likely to be
tance in order to cancel the bias. The approximation q(y) of a higher than the measure of Equation 13, which means that we
given vector y is obtained, in the case of the

product quantizer, penalize the vectors with rare indexes. Note that the correcting
from the subquantizers indexes qj uj (y) , j = 1 . . . m. The term is independent of the query in the asymmetric version.
quantization index identifies the cells Vi in which y lies. We
In our experiments, we observe that the correction returns
can then compute the expected squared distance ẽ x, q(y) inferior results on average. Therefore, we advocate the use of
 6

biased
In addition to the coarse quantizer, we adopt a strategy sim-
empirical probability distribution function

unbiased ilar to that proposed in [20], i.e., the description of a vector is

0.1 refined 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 with
a vector y, the product quantizer qp is used to encode the
residual vector
0.05 r(y) = y − qc (y), (28)
corresponding to the offset in the Voronoi cell. The energy
of the residual vector is small compared to that of the vector
itself. The vector is approximated by
0

-0.3 -0.2 -0.1 0 0.1 0.2 0.3 ÿ , qc (y) + qp y − qc (y) . (29)
difference: estimator - d(x,y)

It is represented by the tuple qc (y), qp (r(y)) . By analogy
Fig. 4. PDF of the error on the distance estimation d − d˜ for the asymmetric with the binary representation of a value, the coarse quantizer
method, evaluated on a set of 10000 SIFT vectors with m = 8 and k∗ = 256. provides the most significant bits, while the product quantizer
The bias (=-0.044) of the estimator d˜ is corrected (=0.002) with the error


quantization term ξ q, q(y) . However, the variance of the error increases code corresponds to the least significant bits.
with this correction: σ 2 (d − ẽ) = 0.00155 whereas σ 2 (d − d)˜ = 0.00146. 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 ÿ:
Equation 13 for the nearest neighbor search. The corrected 

¨ y) = d(x, ÿ) = d x − qc (y), qp y − qc (y) .
d(x, (30)
version is useful only if we are interested in the distances
themselves.
Denoting by qp j the j th subquantizer, we use the following
IV. N ON EXHAUSTIVE SEARCH decomposition to compute this estimator efficiently:
Approximate nearest neighbor search with product quantiz-
X 
 2
¨ y)2 =


d(x, d uj x − qc (y) , qp j uj (y − qc (y)) . (31)
ers is fast (only m additions are required per distance calcu- j
lation) and reduces significantly the memory requirements for
storing the descriptors. Nevertheless, the search is exhaustive. Similar to the ADC strategy, for each subquantizer qp j the


The method remains scalable in the context of a global image distances between the partial residual vector uj x − qc (y)
description [17], [19]. However, if each image is described by and all the centroids cj,i of qp j are preliminarily computed
a set of local descriptors, an exhaustive search is prohibitive, and stored.
as we need to index billions of descriptors and to perform The product quantizer is learned on a set of residual
multiple queries [20]. vectors collected from a learning set. Although the vectors
To avoid exhaustive search we combine an inverted are quantized to different indexes by the coarse quantizer, the
file system [26] with the asymmetric distance computa- resulting residual vectors are used to learn a unique product
tion (IVFADC). An inverted file quantizes the descriptors and quantizer. We assume that the same product quantizer is
stores image indices in the corresponding lists, see the step accurate when the distribution of the residual is marginalized
“coarse quantizer” in Figure 5. This allows rapid access to a over all the Voronoi cells. This probably gives inferior results
small fraction of image indices and was shown successful for to the approach consisting of learning and using a distinct
very large scale search [26]. Instead of storing an image index product quantizer per Voronoi cell. However, this would be
only, we add a small code for each descriptor, as first done computationally expensive and would require storing k 0 prod-
in [20]. Here, we encode the difference between the vector and uct quantizer codebooks, i.e., k 0 ×d×k ∗ floating points values,
its corresponding coarse centroid with a product quantizer, see which would be memory-intractable for common values of k 0 .
Figure 5. This approach significantly accelerates the search at
the cost of a few additional bits/bytes per descriptor. Further- B. Indexing structure
more, it slightly improves the search accuracy, as encoding the
residual is more precise than encoding the vector itself. We use the coarse quantizer to implement an inverted file
structure as an array of lists L1 . . . Lk0 . If Y is the vector
dataset to index, the list Li associated with the centroid ci of
A. Coarse quantizer, locally defined product quantizer qc stores the set {y ∈ Y : qc (y) = ci }.
Similar to the “Video-Google” approach [26], a codebook
is learned using k-means, producing a quantizer qc , referred to In inverted list Li , an entry corresponding to y contains a
as the coarse quantizer in the following. For SIFT descriptors, vector identifier and the encoded residual qp (r(y)):
the number k 0 of centroids associated with qc typically ranges field length (bits)
from k 0 = 1 000 to k 0 = 1 000 000. It is therefore small identifier 8–32
compared to that of the product quantizers used in Section III. code mdlog2 k∗ e
 7

Figure 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 quan-
tizer, 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 and
the binary code (the product quantizer’s indexes).

Searching the nearest neighbor(s) of a query x consists of

1) quantize x to its w nearest neighbors in the codebook qc ;
For the sake of presentation, in the two next steps we
simply denote by r(x) the residuals associated with
these w assignments. The two steps are applied to all
w assignments.

2
2) compute the squared distance d uj (r(x)), cj,i for each
subquantizer j and each of its centroids cj,i ;
Fig. 5. Overview of the inverted file with asymmetric distance computation 3) compute the squared distance between r(x) and all the
(IVFADC) indexing system. Top: insertion of a vector. Bottom: search.
indexed vectors of the inverted list. Using the subvector-
to-centroid distances computed in the previous step,
this consists in summing up m looked-up values, see
The identifier field is the overhead due to the inverted file
Equation 31;
structure. Depending on the nature of the vectors to be
stored, the identifier is not necessarily unique. For instance, 4) select the K nearest neighbors of x based on the
to describe images by local descriptors, image identifiers can estimated distances. This is implemented efficiently by
replace vector identifiers, i.e., all vectors of the same image maintaining a Maxheap structure of fixed capacity, that
have the same identifier. Therefore, a 20-bit field is sufficient to stores the K smallest values seen so far. After each
identify an image from a dataset of one million. This memory distance calculation, the point identifier is added to the
cost can be reduced further using index compression [27], [28], structure only if its distance is below the largest distance
which may reduce the average cost of storing the identifier in the Maxheap.
to about 8 bits, depending on parameters2 . Note that some Only Step 3 depends on the database size. Compared with
geometrical information can also be inserted in this entry, as ADC, the additional step of quantizing x to qc (x) consists
proposed in [20] and [27]. in 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
C. Search algorithm faster than ADC, as shown in the next section.
The inverted file is the key to the non-exhaustive version
V. E VALUATION OF NN SEARCH
of our method. When searching the nearest neighbors of a
vector x, the inverted file provides a subset of Y for which In this section, we first present the datasets used for the
distances are estimated: only the inverted list Li corresponding evaluation3 . We then analyze the impact of the parameters
to qc (x) is scanned. for SDC, ADC and IVFADC. Our approach is compared to
However, x and its nearest neighbor are often not quantized three state-of-the-art methods: spectral hashing [19], Hamming
to the same centroid, but to nearby ones. To address this embedding [20] and FLANN [9]. Finally, we evaluate the
problem, we use the multiple assignment strategy of [29]. The complexity and speed of our approach.
query x is assigned to w indexes instead of only one, which
correspond to the w nearest neighbors of x in the codebook A. Datasets
of qc . All the corresponding inverted lists are scanned. Multiple We perform our experiments on two datasets, one with local
assignment is not applied to database vectors, as this would SIFT descriptors [23] and the other with global color GIST
increase the memory usage. descriptors [18]. We have three vector subsets per dataset:
2 An average cost of 11 bits is reported in [27] using delta encoding and 3 Both the software and the data used in these experiments are available at
Huffman codes. [Link]
 8

learning, database and query. Both datasets were constructed

using publicly available data and software. For the SIFT SDC
descriptors, the learning set is extracted from Flickr images 1 4096

and the database and query descriptors are from the INRIA 1024

Holidays images [20]. For GIST, the learning set consists 0.8 256
of the first 100k images extracted from the tiny image set
of [16]. The database set is the Holidays image set combined

recall@100
0.6
with Flickr1M used in [20]. The query vectors are from the 64

Holidays image queries. Table III summarizes the number of

0.4
descriptors extracted for the two datasets.
k*=16 m=1
0.2 m=2
vector dataset: SIFT GIST m=4
descriptor dimensionality d 128 960 m=8
learning set size 100,000 100,000 m=16
0
database set size 1,000,000 1,000,991 0 16 32 64 96 128 160
queries set size 10,000 500 code length (bits)
TABLE III
S UMMARY OF THE SIFT AND GIST DATASETS .
ADC
1024 4096
1
256
The search quality is measured with recall@R, i.e., the
proportion of query vectors for which the nearest neighbor 0.8 64

is ranked in the first R positions. This measure indicates

the fraction of queries for which the nearest neighbor is
recall@100
0.6
retrieved correctly, if a short-list of R vectors is verified
using Euclidean distances. Furthermore, the curve obtained by 0.4
k*=16
varying R corresponds to the distribution function of the ranks,
m=1
and the point R=1 corresponds to the “precision” measure used m=2
0.2
in [9] to evaluate ANN methods. m=4
m=8
In practice, we are often interested in retrieving the K m=16
nearest neighbors (K > 1) and not only the nearest neighbor. 0
0 16 32 64 96 128 160
We do not include these measures in the paper, as we observed code length (bits)
that the conclusions for K=1 remain valid for K > 1.
Fig. 6. SDC and ADC estimators evaluated on the SIFT dataset: recall@100
B. Memory vs search accuracy: trade-offs 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
The product quantizer is parametrized by the number of (m=16,k∗ =4096) gives recall@100=1 for both SDC and ADC.
subvectors m and the number of quantizers per subvector k ∗ ,
producing a code of length m × log2 k ∗ . Figure 6 shows the
trade-off between code length and search quality for our SIFT
descriptor dataset. The quality is measured for recall@100 for IVFADC
the ADC and SDC estimators, for m ∈ {1, 2, 4, 8, 16} and 1
k ∗ ∈ {24 , 26 , 28 , 210 , 212 }. As for the quantizer distortion in
Figure 1, one can observe that for a fixed number of bits, it 0.8
is better to use a small number of subquantizers with many
centroids than to have many subquantizers with few bits. How-
recall@100

0.6
ever, the comparison also reveals that MSE underestimates, for
4 8 16
a fixed number of bits, the quality obtained for a large number 0.4
2
of subquantizers against using more centroids per quantizer. m=1

As expected, the asymmetric estimator ADC significantly k’=1024, w=1

0.2
outperforms SDC. For m=8 we obtain the same accuracy for k’=1024, w=8
k’=8192, w=8
ADC and k ∗ =64 as for SDC and k ∗ =256. Given that the k’=8192, w=64
efficiency of the two approaches is equivalent, we advocate 0
0 16 32 64 96 128
not to quantize the query when possible, but only the database code length (bits)
elements.
Figure 7 evaluates the impact of the parameters for the Fig. 7. SIFT dataset: recall@100 for the IVFADC approach as a function of
IVFADC method introduced in Section IV. For this approach, the memory usage for k∗ =256 and varying values of m = {1, 2, 4, 8, 16},
we have, in addition, to set the codebook size k 0 and the k0 = {1024, 8192} and w = {1, 8, 64}.
number of neighboring cells w visited during the multiple
 9

assignment. We observe that the recall@100 strongly depends the same orientation. On SIFT descriptors, this is a slightly
on these parameters, and that increasing the code length is less efficient structure, probably because the natural order
useless if w is not big enough, as the nearest neighbors which corresponds to spatially related components. On GIST, this
are not assigned to one of the w centroids associated with the choice significantly improves the performance. Therefore, we
query are definitely lost. use this ordering in the following experiments.
This approach is significantly more efficient than SDC and
ADC on large datasets, as it only compares the query to a small Discussion: A method that automatically groups the compo-
fraction of the database vectors. The proportion of the dataset nents could further improve the results. This seems particularly
to visit is roughly linear in w/k 0 . For a fixed proportion, it is important if we have no prior knowledge about the relationship
worth using higher values of k 0 , as this increases the accuracy, between the components as in the case of bag-of-features.
as shown by comparing, for the tuple (m, w), the parameters A possible solution is the minimum sum-squared residue co-
(1024, 1) against (8192, 8) and (1024, 8) against (8192, 64). clustering [30] algorithm.

D. Comparison with the state of the art

C. Impact of the component grouping
The product quantizer defined in Section II creates the
Comparison with Hamming embedding methods: We com-
subvectors by splitting the input vector according to the
pare our approach to spectral hashing (SH) [19], which maps
order of the components. However, vectors such as SIFT
vectors to binary signatures. The search consists in comparing
and GIST descriptors are structured because they are built
the Hamming distances between the database signatures and
as concatenated orientation histograms. Each histogram is
the query vector signature. This approach was shown to
computed on grid cells of an image patch. Using a product
outperform the restricted Boltzmann machine of [17]. We have
quantizer, the bins of a histogram may end up in different
used the publicly available code. We also compare to the
quantization groups.
Hamming embedding (HE) method of [20], which also maps
The natural order corresponds to grouping consecutive vectors to binary signatures. Similar to IVFADC, HE uses
components, as proposed in Equation 8. For the SIFT descrip- an inverted file, which avoids comparing to all the database
tor, this means that histograms of neighboring grid cells are elements.
quantized together. GIST descriptors are composed of three Figures 8 and 9 show, respectively for the SIFT and the
320-dimension blocks, one per color channel. The product GIST datasets, the rank repartition of the nearest neighbors
quantizer splits these blocks into parts. when using a signature of size 64 bits. For our product
quantizer we have used m = 8 and k ∗ = 256, which
SIFT GIST give similar results in terms of run time. All our approaches
m 4 8 8 significantly outperform spectral hashing4 on the two datasets.
natural 0.593 0.921 0.338 To achieve the same recall as spectral hashing, ADC returns
random 0.501 0.859 0.286
structured 0.640 0.905 0.652
an order of magnitude less vectors.
The best results are obtained by IVFADC, which for low
TABLE IV ranks provides an improvement over ADC, and significantly
I MPACT OF THE DIMENSION GROUPING ON THE RETRIEVAL outperforms spectral hashing. This strategy avoids the exhaus-
PERFORMANCE OF ADC ( RECALL @100, k ∗ =256). tive search and is therefore much faster, as discussed in the
next subsection. This partial scan explains why the IVFADC
and HE curves stop at some point, as only a fraction of the
database vectors are ranked. Comparing these two approaches,
To evaluate the influence of the grouping, we modify HE is significantly outperformed by IVFADC. The results of
the uj operators in Equation 8, and measure the impact of HE are similar to spectral hashing, but HE is more efficient.
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 parameters Comparison with FLANN: The approximate nearest-
k ∗ =256 and m ∈ {4, 8}. neighbor search technique of Muja & Lowe [9] is based
Overall, the choice of the components appears to have a on hierarchical structures (KD-trees and hierarchical k-means
significant impact of the results. Using a random order instead trees). The software package FLANN automatically selects the
of the natural order leads to poor results. This is true even for best algorithm and parameters for a given dataset. In contrast
GIST, for which the natural order is somewhat arbitrary. with our method and spectral hashing, all vectors need to
The “structured” order consists in grouping together dimen- remain in RAM as the method includes a re-ranking stage
sions that are related. For the m = 4 SIFT quantizer, this that computes the real distances for the candidate nearest
means that the 4×4 patch cells that make up the descriptor [23] neighbors.
are grouped into 4 2 × 2 blocks. For the other two, it groups The evaluation is performed on the SIFT dataset by measur-
together dimensions that have have the same index modulo 8. ing the 1-recall@1, that is, the average proportion of true NNs
The orientation histograms of SIFT and most of GIST’s have 8 4 In defense of [17], [19], which can be learned for arbitrary distance
bins, so this ordering quantizes together bins corresponding to measures, our approach is adapted to the Euclidean distance only.
 10

SIFT, 64-bit codes

100 FLANN
1 IVFADC, R=10
IVFADC, R=100
IVFADC, R=1000 16/128
0.8
16/1024

search time (s)

16/128

0.6
recall@R

16/1024
10
4/128
0.4 SDC 4/1024
ADC 4/1024
IVFADC w=1
0.2 IVFADC w=16
HE w=1
HE w=16
spectral hashing
0 1
1 10 100 1k 10k 100k 1M 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1
R 1-recall at 1

Fig. 10. IVFADC vs FLANN: trade-offs between search quality (1-recall@1)

Fig. 8. SIFT dataset: recall@R for varying values of R. Comparison of and search time. The IVFADC method is parametrized by the shortlist size R
the different approaches SDC, ADC, IVFADC, spectral hashing [19] and used for re-ranking the vector with the L2 distance, and the two parameters
HE [20]. Parameters: m=8 and k∗ =256 for SDC/ADC, k0 =1024 for HE [20] w and k0 of the inverted file, which correspond to the number of assignments
and IVFADC. and to the number of coarse centroids.

GIST, 64-bit codes

1 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
0.8
points. Moreover, our method has a much smaller memory
footprint than FLANN: the indexing structure occupies less
0.6
recall@R

than 25 MB, while FLANN requires more than 250 MB of

RAM. Note, however, that both are negligible compared to the
0.4 memory occupied by the vectors in the case of large datasets.
SDC
ADC On such a scale, the re-ranking stage is not feasible and only
0.2 IVFADC w=1 memory-aware approaches (HE, SH and our methods) can be
IVFADC w=8
IVFADC w=64 used.
spectral hashing
0
1 10 100 1k 10k 100k 1M E. Complexity and speed
R Table V reports the search time of our methods. For
reference, we report the results obtained with the spectral
Fig. 9. GIST dataset: recall@R for varying values of R. Comparison of hashing algorithm of [19] on the same dataset and machine
the different approaches SDC, ADC, IVFADC and spectral hashing [19].
Parameters: m=8, k∗ =256 for SDC/ADC and k0 = 1 024 for IVFADC. (using only one core). Since we use a separate learning set,
we use the out-of-sample evaluation of this algorithm. Note
that for SH we have re-implemented the Hamming distance
ranked first in the returned vectors. This measure is referred computation in C in order to have the approaches similarly
to as precision in [9]. optimized. The algorithms SDC, ADC and SH have similar
For the sake of comparison with FLANN, we added a run times. IVFADC significantly improves the performance by
verification stage to our IVFADC method: IVFADC queries avoiding an exhaustive search. Higher values of k 0 yield higher
return a shortlist of R candidate nearest neighbors using the search efficiencies for large datasets, as the search benefits
distance estimation. The vectors in the shortlist are re-ordered from parsing a smaller fraction of the memory. However, for
using the real distance, as done in [7], [9], and the closest small datasets, the complexity of the coarse quantizer may be
one is returned. Note that, in this experimental setup, all the the bottleneck if k 0 × D > n/k 0 when using a exhaustive
vectors are stored in main memory. This requirement seriously assignment for qc . In that case the ADC variant may be
limits the scale on which re-ordering can be used. preferred. For large datasets and using an efficient assignment
The IVFADC and FLANN methods are both evaluated at strategy for the coarse quantizer, higher values of k 0 generally
different operating points with respect to precision and search lead to better efficiency, as first shown in [15]. In this work, the
time. For FLANN, the different operating points are obtained authors propose a hierarchical quantizer to efficiently assign
with parameters generated automatically for various target descriptors to one million centroids.
precisions. For IVFADC, they are obtained by varying the
number k 0 of coarse centroids, the number w of assignments F. Large-scale experiments
and the short-list size R. The product quantizer is generated To evaluate the search efficiency of the product quantizer
using k ∗ =256 and m=8, i.e., 64-bit codes. This choice is method on larger datasets we extracted about 2 billion SIFT
 11

method parameters search average number of recall@100

time (ms) code comparisons
SDC 16.8 1 000 991 0.446
ADC 17.2 1 000 991 0.652
IVFADC k0 = 1 024, w=1 1.5 1 947 0.308
k0 = 1 024, w=8 8.8 27 818 0.682
k0 = 1 024, w=64 65.9 101 158 0.744
k0 = 8 192, w=1 3.8 361 0.240
k0 = 8 192, w=8 10.2 2 709 0.516
k0 = 8 192, w=64 65.3 19 101 0.610
SH 22.7 1 000 991 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 0.8
HE IVF+HE 64 bits
IVFADC IVFADC 64 bits (m=8)
3 0.75
IVF+HE 32 bits
search time (ms/vector)

IVFADC 32 bits (m=4)

0.7
2.5
0.65
2
mAP

0.6
1.5
0.55
1
0.5
0.5
0.45

0 0.4
10M 100M 1G 1k 10k 100k 1M
database size database size (number of images)

Fig. 11. Search times for SIFT descriptors in datasets of increasing sizes, Fig. 12. Comparison of IVFADC and the Hamming Embedding method
with two search methods. Both use the same 20 000-word codebook, w = 1, of [20]. mAP for the Holidays dataset as function of the number of distractor
and 64-bit signatures. images (up to 1 million).

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 compare our method with the HE method of [20] on the INRIA
particular, the amount of memory that is scanned for each Holidays dataset, using the pre-processed set of descriptors
method and the cost of the coarse quantization are the same. available online. The comparison is focused on large scale
The query times per descriptor are shown on Figure 11. indexing, i.e., we do not consider the impact of a post-
The cost of the extra quantization step required by IVFADC verification step [23], [31] or geometrical information [20].
appears clearly for small database sizes. For larger scales,
the distance computation with the database vectors become Figure 12 shows the search performance in terms of mean
preponderant. The processing that is applied to each element of average precision as a function of the size of the dataset. We
the inverted lists is approximately as expensive in both cases. have used the same coarse quantizer (k 0 =20,000) and a single
For HE, it is a Hamming distance computation, implemented assignment strategy (w=1) for both the approaches, and fixed
as 8 table lookups. For IVFADC it is a distance computation k ∗ =256 for IVFADC. For a given number of bits (32 or 64),
that is also performed by 8 table lookups. Interestingly, the we have selected the best choice of the Hamming threshold
floating point operations involved in IVFADC are not much for HE. Similarly, we have adjusted the number of nearest
more expensive than the simple binary operations of HE. neighbors to be retrieved for IVFADC.

One can observe that the gain obtained by IVFADC is

G. Image search significant. For example, for one million distractors, the mAP
We have evaluated our method within an image search of 0.497 reported in [20] with 64-bit signatures is improved
system based on local descriptors. For this evaluation, we to 0.517.
 12

VI. C ONCLUSION [13] B. Matei, Y. Shan, H. Sawhney, Y. Tan, R. Kumar, D. Huber, and
M. Hebert, “Rapid object indexing using locality sensitive hashing and
We have introduced product quantization for approximate joint 3D-signature space estimation,” IEEE Transactions on Pattern
nearest neighbor search. Our compact coding scheme pro- Analysis and Machine Intelligence, vol. 28, pp. 1111 – 1126, July 2006.
vides an accurate approximation of the Euclidean distance. [14] C. Silpa-Anan and R. Hartley, “Optimized kd-trees for fast image de-
scriptor matching,” in Proceedings of the IEEE Conference on Computer
Moreover, it is combined with an inverted file system to Vision and Pattern Recognition, 2008.
avoid exhaustive search, resulting in high efficiency. Our [15] D. Nistér and H. Stewénius, “Scalable recognition with a vocabulary
approach significantly outperforms the state of the art in terms tree,” in Proceedings of the IEEE Conference on Computer Vision and
Pattern Recognition, pp. 2161–2168, 2006.
of the trade-off between search quality and memory usage. [16] A. Torralba, R. Fergus, and W. T. Freeman, “80 million tiny images: a
Experimental results for SIFT and GIST image descriptors are large database for non-parametric object and scene recognition,” IEEE
excellent and show that grouping the components based on Transactions on Pattern Analysis and Machine Intelligence, vol. 30,
pp. 1958–1970, November 2008.
our prior knowledge of the descriptor design further improves [17] A. Torralba, R. Fergus, and Y. Weiss, “Small codes and large databases
the results. The scalability of our approach is validated on a for recognition,” in Proceedings of the IEEE Conference on Computer
dataset of two billion vectors. 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 Com-
R EFERENCES puter Vision, vol. 42, no. 3, pp. 145–175, 2001.
[19] Y. Weiss, A. Torralba, and R. Fergus, “Spectral hashing,” in Advances
[1] H. Jégou, M. Douze, and C. Schmid, “Searching with quantization: in Neural Information Processing Systems, 2008.
approximate nearest neighbor search using short codes and distance [20] H. Jégou, M. Douze, and C. Schmid, “Hamming embedding and weak
estimators,” Tech. Rep. RR-7020, INRIA, August 2009. geometric consistency for large scale image search,” in Proceedings of
[2] H. Sandhawalia and H. Jégou, “Searching with expectations,” in IEEE the European Conference on Computer Vision, October 2008.
International Conference on Acoustics Speech and Signal Processing, [21] M. Douze, H. Jégou, H. Singh, L. Amsaleg, and C. Schmid, “Evaluation
Signal Processing, pp. 1242–1245, IEEE, March 2010. of GIST descriptors for web-scale image search,” in Proceedings of the
[3] K. Beyer, J. Goldstein, R. Ramakrishnan, and U. Shaft, “When is International Conference on Image and Video Retrieval, 2009.
”nearest neighbor” meaningful?,” in Proceedings of the International [22] J. Philbin, O. Chum, M. Isard, J. Sivic, and A. Zisserman, “Lost in
Conference on Database Theory, pp. 217–235, August 1999. quantization: Improving particular object retrieval in large scale image
[4] C. Böhm, S. Berchtold, and D. Keim, “Searching in high-dimensional databases,” in Proceedings of the IEEE Conference on Computer Vision
spaces: Index structures for improving the performance of multimedia and Pattern Recognition, 2008.
databases,” ACM Computing Surveys, vol. 33, pp. 322–373, October [23] D. Lowe, “Distinctive image features from scale-invariant keypoints,”
2001. International Journal of Computer Vision, vol. 60, no. 2, pp. 91–110,
[5] J. Friedman, J. L. Bentley, and R. A. Finkel, “An algorithm for 2004.
finding best matches in logarithmic expected time,” ACM Transaction [24] R. M. Gray and D. L. Neuhoff, “Quantization,” IEEE Transactions on
on Mathematical Software, vol. 3, no. 3, pp. 209–226, 1977. Information Theory, vol. 44, pp. 2325–2384, Oct. 1998.
[6] R. Weber, H.-J. Schek, and S. Blott, “A quantitative analysis and [25] D. E. Knuth, The Art of Computer Programming, Sorting and Searching,
performance study for similarity-search methods in high-dimensional vol. 3. Addison Wesley, 2 ed., 1998.
spaces,” in Proceedings of the International Conference on Very Large [26] J. Sivic and A. Zisserman, “Video Google: A text retrieval approach
DataBases, pp. 194–205, 1998. to object matching in videos,” in Proceedings of the International
[7] M. Datar, N. Immorlica, P. Indyk, and V. Mirrokni, “Locality-sensitive Conference on Computer Vision, pp. 1470–1477, 2003.
hashing scheme based on p-stable distributions,” in Proceedings of the [27] M. Perdoch, O. Chum, and J. Matas, “Efficient representation of local
Symposium on Computational Geometry, pp. 253–262, 2004. geometry for large scale object retrieval,” in Proceedings of the IEEE
[8] A. Gionis, P. Indyk, and R. Motwani, “Similarity search in high Conference on Computer Vision and Pattern Recognition, June 2009.
dimension via hashing,” in Proceedings of the International Conference [28] H. Jégou, M. Douze, and C. Schmid, “Packing bag-of-features,” in Pro-
on Very Large DataBases, pp. 518–529, 1999. ceedings of the International Conference on Computer Vision, September
[9] M. Muja and D. G. Lowe, “Fast approximate nearest neighbors with 2009.
automatic algorithm configuration,” in Proceedings of the International [29] H. Jégou, H. Harzallah, and C. Schmid, “A contextual dissimilarity
Conference on Computer Vision Theory and Applications, 2009. measure for accurate and efficient image search,” in Proceedings of the
[10] B. Kulis and K. Grauman, “Kernelized locality-sensitive hashing for IEEE Conference on Computer Vision and Pattern Recognition, 2007.
scalable image search,” in Proceedings of the International Conference [30] H. Cho, I. Dhillon, Y. Guan, and S. Sra, “Minimum sum-squared
on Computer Vision, October 2009. residue co-clustering of gene expression data,” in SIAM International
[11] G. Shakhnarovich, T. Darrell, and P. Indyk, Nearest-Neighbor Methods Conference on Data Mining, pp. 114–125, April 2004.
in Learning and Vision: Theory and Practice, ch. 3. MIT Press, March [31] J. Philbin, O. Chum, M. Isard, J. Sivic, and A. Zisserman, “Object re-
2006. trieval with large vocabularies and fast spatial matching,” in Proceedings
[12] Y. Ke, R. Sukthankar, and L. Huston, “Efficient near-duplicate detection of the IEEE Conference on Computer Vision and Pattern Recognition,
and sub-image retrieval,” in ACM International conference on Multime- 2007.
dia, pp. 869–876, 2004.