Compact Hashing with Joint Optimization of Search Accuracy and Time
Junfeng He
Columbia University
New York, New York, 10027, USA

Regunathan Radhakrishnan
Dolby Laboratories
San Francisco, CA, 94103, USA

jh2700@columbia.edu

regu.r@dolby.com

Shih-Fu Chang
Columbia University
New York, New York, 10027, USA

Claus Bauer
Dolby Laboratories
San Francisco, CA, 94103, USA

sfchang@ee.columbia.edu

cb@dolby.com

Abstract

the data dimension is high, sometimes even worse than linear scan. Recently, many hash coding based algorithms
have been proposed to handle similarity search of high dimensional data [3, 5, 6, 11, 10, 8, 13, 14, 12, 17]. Up tonow, there are two main groups of hashing methods: hashing based on random projections and hashing for compact
codes.
Among the popular randomized hash based techniques
is the locality-sensitive hashing (LSH) [3]. In LSH, random
vectors (with some specific distribution) are used as the projection bases. As an important property of LSH, points withhigh feature similarity are proven with a high probability toobtain the same hash code, which guarantees a asymptotic
theoretical property of sublinear search time. Variations of
LSH algorithms have been proposed in recent literatures toexpand its usage to the cases of inner products [4], 𝐿𝑝 norms[5], and learned metrics [7], etc.
Despite of its success, one arguable disadvantage of LSH
is the inefficiency of the hash codes. Since the hash functions in LSH are randomly generated and independent of
the data, it is not very efficient. Usually, it needs long codesin each hash table to achieve an acceptable accuracy. However, since the collision probability decreases very fast (exponentially) with the length of the code, when long codes
are used in one hash table, then many hash tables are often
needed to get a good recall. This would heavily increase
the requirement of storage, causing problems for very large
scale applications. As demonstrated in [11, 10], if we wantto scale up to Internet size databases (80 millions samples
in their experiments), fitting those large scale data into adesktop memory, e.g., a few Gigabytes, means we can only
use a small number of bits per sample. In this case, LSH
with short codes usually has a poor performance in terms ofprecision [11, 10]. So, many recent research works focus on
how to generate short compact hash codes [8, 11, 10, 13, 17]

Similarity search, namely, finding approximate nearest
neighborhoods, is the core of many large scale machinelearning or vision applications. Recently, many research
results demonstrate that hashing with compact codes canachieve promising performance for large scale similarity
search. However, most of the previous hashing methods
with compact codes only model and optimize the search accuracy. Search time, which is an important factor for hashing in practice, is usually not addressed explicitly. In thispaper, we develop a new scalable hashing algorithm withjoint optimization of search accuracy and search time simultaneously. Our method generates compact hash codesfor data of general formats with any similarity function. We
evaluate our method using diverse data sets up to 1 million samples (e.g., web images). Our comprehensive results
show the proposed method significantly outperforms several
state-of-the-art hashing approaches.

1. Introduction and related work
The explosive growth of data on the Internet has givenrise to many interesting applications, for many of which acore technical issue involves fast nearest neighbor searchover a large data collection in response to a query. Bruteforce solutions using linear scan obviously do not scale up
when the data size grows due to the excessive cost in storage
and processing. Recent works [1, 2, 3, 5, 6, 11, 10, 8, 13,
14, 12, 17], have reported interesting results in finding approximate solutions over large data sets. Early works in thisfield focused on finding appropriate spatial partitions of thefeature space via various tree structures[1, 2]. With good
performance for low dimensional data, such tree-based approaches however, are known to degrade significantly when753

2. Modeling search accuracy and search timefor each sample with only one hash table or few hash tables, which can still keep good precision. A well-known
method for hashing with compact codes is spectral hashing (SH)[11], which has been shown to achieve much better
performance than LSH in applications with compact hashcodes[11].
In this paper, we also aim at developing a hashingmethod to generate compact hash codes. However, previous
hashing methods with compact codes, e.g., [11, 10, 8], only
consider and optimize the search accuracy. Search time,
which is also an important factor, is usually not discussed
explicitly. In this paper, we develop a new hashing algorithm with a joint framework to explicitly optimize both
search accuracy and search time together. To the best to
our knowledge, our paper presents the first result on such
explicit joint optimization.
Moreover, most previous hashing methods assume data
is represented in a vector format and can be embedded as
points in a vector space. However, in applications involving
computer vision, multimedia, biology, or Web, etc., structured data are quite common [27, 25]. Though several recent previous kernel hashing methods, e.g.,[13, 14, 12] havediscussed hashing for non-vectors, but they all need a Mercer kernel function. In our proposed method, any similarity/proximity function (including kernel function of course)
is applicable, and thus broadens the possible applications ofour method.
The main contributions of our work are

2.1. Background and notations
We first assume our data are vectors. The extension to
non-vector data will be discussed in section 3.5. Suppose
feature matrix for training data 𝑋 is a 𝑑 by 𝑁 matrix. 𝑑 isthe dimension of feature, and 𝑁 is the number of training
data. Denote 𝑋𝑝 as the 𝑝-th column of 𝑋, which is thefeature vector for 𝑝-th data in the training set. Like in mostmachine learning problems, the training samples 𝑋𝑝 , 𝑝 =1, ..., 𝑁 are assumed to be i.i.d. sampled from a random
feature vector 𝑥.
Suppose 𝑀 is the number of hash bits, and𝐻 : 𝑅𝑑 × {−1, 1}𝑀represents the hash function we are looking for, which maps
a d-dimension vector to M binary bits. Intuitively, 𝐻 consists of 𝑀 functions 𝐻𝑗 : 𝑅𝑑 × {−1, 1}, 𝑗 = 1, ..., 𝑀 , eachof which maps a 𝑑-dimension vector to 1 binary bit.
Denote 𝑌𝑝 = 𝐻(𝑋𝑝 ), the hash bits for 𝑋𝑝 , which isa 𝑀 -dimensional vector. Since {𝑋𝑝 , 𝑝 = 1, ..., 𝑁 } arei.i.d. sampled from a random vector 𝑥, {𝑌𝑝 = 𝐻(𝑋𝑝 ), 𝑝 =1, ..., 𝑁 } are i.i.d. samples from 𝑦 = 𝐻(𝑥). Here 𝑦 is a
𝑀 -dimension random binary vector. 𝑦𝑘 is denoted as 𝑘-th
dimension of 𝑦. There are 𝑀 bits, so 𝑦 has 2𝑀 outcomes:𝑎𝑖 ∈ {−1, 1}𝑀 , 𝑖 = 1, ..., 2𝑀

1. We analyze and model the search time, and prove thatit can be minimized if the hashing buckets are perfectly balanced, i.e., contain equal number of samples. A maximum entropy criterion and equivalently
a minimum mutual information criterion are proposedto achieve the bucket balancing.

Each 𝑎𝑖 is associated with a unique 𝑀 -bit hash code andthus represents a distinct hash bucket. There are 2𝑀 buckets
in total. Denote 𝑁𝑖 as the number of training samples inbucket 𝑖, for 𝑖 = 1, ..., 2𝑀 .

2.2. Similarity preserving for search accuracy
2. A joint optimization framework, which considers boththe search accuracy and the search time, is provided.

For two data samples 𝑋𝑝 and 𝑋𝑞 in the training set, suppose 𝑊𝑝𝑞 is the similarity between 𝑋𝑝 and 𝑋𝑞 . Similarity
𝑊𝑝𝑞 can come from feature similarity or label consistency,
etc., depending on the applications. The only requirement
for 𝑊 is the symmetry.
To obtain a good search accuracy, hashing methods usually have a common property: similarity preserving, i.e.,
similar points are supposed to have similar hash codes withsmall hamming distance. One popular approach to preserver similarity is to minimize the following [11]:3. Relaxation and implementation techniques are proposed such that the joint optimization problem can besolved efficiently, which outperforms the state-of-theart hashing methods with compact codes.
4. The proposed method is generalized to handle data ofgeneral formats including graphs, sequences, sets, etc.
besides vectors, with any kernel function or similarity/proximity function defined.

𝐷(𝑌 ) =We have evaluated our algorithm on several databaseswith a data size ranging from several thousand to one million. Our proposed method significantly outperforms several stat-of-the-art approaches for all the data sets/tasks, including implementations with both single hash table andmultiple hash tables.

∑∑

𝑊𝑝𝑞 ∣∣𝑌𝑝 − 𝑌𝑞 ∣∣2(1)

𝑝=1,...,𝑁 𝑞=1,...,𝑁

where 𝑌 is the set of all 𝑌𝑝 . With this criterion, samples
with high similarity, i.e., larger 𝑊𝑝𝑞 , are supposed to have2
similar hash codes, i.e., smaller ∣∣𝑌𝑝 − 𝑌𝑞 ∣∣ .
754

point falls into bucket 𝑖 is approximated by the frequency𝑖𝑝𝑖 ≈ 𝑁𝑁 , and moreover the processing time for all samples in bucket 𝑖 is 𝑆𝑖 = 𝑐𝑁𝑖 , where 𝑐 is a constant that is
related to the time of processing (e.g., loading or displaying) one individual sample. So the expected time to search
𝑀
𝑀
2∑
2∑
𝑝𝑖 𝑆𝑖 ≈ 𝑁𝑐
𝑁𝑖2 which only depends
one bucket is:

2.3. Minimal mutual information criterion forsearch time
Generally speaking, the search process for hash based
methods would consist of 3 steps:
1. Compute the hash bits for the query.
2. Find the bucket with the same hash bits as the query.
3. Load samples into memory which are in the selectedbucket. Rerank them if necessary.

on

𝑀
2
∑

𝑖=1

𝑖=1

𝑁𝑖2 .

𝑖=1

To minimize the average time complexity instep 3 equals to minimizing

It is easy to see that similarity preserving alone does
not guarantee a good hash/serach method. For example,
when the similarity matrix 𝑊 contains only non-negative
elements, as an extreme case, one can always assign thesame bits for every data, and hence 𝐷(𝑌 ) would be 0 andhence minimized. However, in this case, for every query,
all the data would be returned, which is the worst case of
search time, and actually equals linear scan. So the objective of hashing should not focus on similarity preserving
only, but aim at a good tradeoff between search accuracy
and search time. We will first model the search time in thissection, and then provide the joint optimization frameworkto achieve tradeoff between search accuracy and time in section 3.1.
Note that the main difference of search time among hashmethods is determined by the difference of time complexity
for step 3, since the time complexity of step 1 and 2 is independent of the database size. Moreover, for large scale applications, usually, samples can not be pre-stored into memory, so step 3 require harddisk operations, while step 1 and 2
only involve memory operations. So the time complexity ofstep 3 is the main factor determining the search time, which
depends on the number of samples in the selected buckets.
The following proposition shows that the time complexity of step 3 is minimized if all buckets contain the same
number of samples, which leads to a minimum mutual information criterion.
Proposition 1: Both the worst case and the averagecase of time complexity in step 3 can be minimized if thebuckets are perfectly balanced, i.e., every bucket contains the same number of samples, which happens when
𝐸𝑛𝑡𝑟𝑜𝑝𝑦(𝑦) is maximized, or equivalently, mathematical expectation 𝐸(𝑦𝑘 ) = 0 for 𝑘 = 1, ..., 𝑀 and the mutual information 𝐼(𝑦1 , ..., 𝑦𝑘 , ..., 𝑦𝑀 ) is minimized.
Proof of Proposition 1: Recall that 𝑁𝑖 is the number oftraining samples in bucket 𝑖, for 𝑖 = 1, ..., 2𝑀 . The worst
case of processing time in step 3 happens when the selectedbucket is the largest bucket, i.e., the bucket with largest 𝑁𝑖 .
So the worst case can be minimized if the buckets are perfectly balanced, i.e., every bucket contains equal number ofsamples, in other words, 𝑁𝑖 = 2𝑁𝑀 , for 𝑖 = 1, ..., 2𝑀 .
Further more, since the offline training data and newquery are i.i.d. sampled, then the probability of the query

𝑀
2
∑

𝑖=1

𝑁𝑖2 , which can be achievedby 𝑁𝑖 = 2𝑁𝑀 for 𝑖 = 1, ..., 2𝑀 , considering

𝑀
2
∑

𝑖=1

𝑁𝑖 = 𝑁 .

In sum, both the worst case and the average case of timecomplexity in step 3 can be minimized if 𝑁𝑖 = 2𝑁𝑀 for𝑖 = 1, ..., 2𝑀 .
𝑀
2∑𝑃 (𝑦 = 𝑎𝑖 ) log 𝑃 (𝑦 =Note that 𝐸𝑛𝑡𝑟𝑜𝑝𝑦(𝑦) = {−𝑎𝑖 )} = {−
𝑁
2𝑀

𝑀
2
∑

𝑖=1

𝑖=1

𝑁𝑖
𝑁𝑖
𝑁 log( 𝑁 )}.

It is easy to see that 𝑁𝑖 =𝑀

for 𝑖 = 1, ..., 2 , if and only if 𝐸𝑛𝑡𝑟𝑜𝑝𝑦(𝑦) gets its
maximum value.
As shown in [19],𝐸𝑛𝑡𝑟𝑜𝑝𝑦(𝑦) =𝑀∑

𝐸𝑛𝑡𝑟𝑜𝑝𝑦(𝑦𝑘 ) − 𝐼(𝑦1 , ..., 𝑦𝑘 , ..., 𝑦𝑀 )

𝑘=1

(2)
where 𝐼() is the mutual information.
So 𝐸𝑛𝑡𝑟𝑜𝑝𝑦(𝑦) would be maximized,
if𝑀∑
𝐸𝑛𝑡𝑟𝑜𝑝𝑦(𝑦𝑘 ) is maximized and 𝐼(𝑦1 , ..., 𝑦𝑘 , ..., 𝑦𝑀 )𝑘=1

is minimized. Moreover, note that 𝑦𝑘 is a binary randomvariable. If the mathematical expectation 𝐸(𝑦𝑘 ) = 0, half
samples would have bit +1 and the other half would havebit −1 for 𝑦𝑘 , which means 𝐸𝑛𝑡𝑟𝑜𝑝𝑦(𝑦𝑘 ) = 1, and ismaximized.
In conclusion, if 𝐸(𝑦𝑘 ) = 0, 𝑘 = 1, ..., 𝑀 and𝐼(𝑦1 , ..., 𝑦𝑘 , ..., 𝑦𝑀 ) is minimized, 𝐸𝑛𝑡𝑟𝑜𝑝𝑦(𝑦) would be
maximized, and the search time would be minimized. This
completes the proof of Proposition 1.
In practice, many hashing methods are applied with a”multi-probe” strategy, i.e., finding multiple buckets withina small hamming distance or sharing the most reliable hashbits of the query, instead of finding only one bucket with thesame hash bits. It is easy to see that the analysis of the worst
case search time is still valid in this case.
Moreover, note that if 𝐼(𝑦1 , ..., 𝑦𝑘 , ..., 𝑦𝑀 ) is minimized,
it means 𝑦1 , ..., 𝑦𝑘 , ..., 𝑦𝑀 are independent 1 . So minimiz-

1 Independence among hash bits are mentioned in spectral hashing[11],but it does not relate independence to search time, and moreover, there isno actual formulation, derivation or algorithm to achieve the independence.

755

3.3. Relaxation for minimizing 𝐼(𝑦1 , ..., 𝑦𝑘 , ..., 𝑦𝑀 )ing mutual information criterion is also to provide the mostcompact and least redundant hash codes.

Denote 𝑧𝑘 = 𝑇𝑘𝑇 𝑥, where 𝑇𝑘 is the 𝑘-th column of 𝑇 .
So 𝑦𝑘 = 𝑠𝑖𝑔𝑛(𝑧𝑘 − 𝑏).
It is easy to see that if 𝑧𝑘 are indepenHence if
dent, 𝑦𝑘 would be independent too.
𝑇=
𝐼(𝑇1𝑇 𝑥, ..., 𝑇𝑘𝑇 𝑥, ..., 𝑇𝑀𝑥)𝐼(𝑧1 , ..., z𝑘 , ..., 𝑧𝑀 )is minimized, 𝐼(𝑦1 , ..., 𝑦𝑘 , ..., 𝑦𝑀 ) would be minimized. So we will minimize 𝐼(𝑧1 , ..., z𝑘 , ..., 𝑧𝑀 ) =𝑇𝑥) instead of 𝐼(𝑦1 , ..., 𝑦𝑘 , ..., 𝑦𝑀 )𝐼(𝑇1𝑇 𝑥, ..., 𝑇𝑘𝑇 𝑥, ..., 𝑇𝑀
in equation (3).
In the field of ICA, independence or mutual information are well studied for a long time. As discussed in [19],minimizing 𝐼(𝑧1 , ..., 𝑧𝑘 , ..., 𝑧𝑀 ) can be well approximated𝑀∑
as minimizing 𝐶0 −∣∣𝑔0 − 𝐸(𝐺(𝑧𝑘 ))∣∣2 , which equals

3. Formulation and relaxation
3.1. Formulation of hashing with joint optimizationof search accuracy and time
Note that 𝐸(𝑦𝑘 ) = 0, 𝑘 = 1, ..., 𝑀 means 𝐸(𝑦) = 0,𝑁∑which can further be rewritten as
𝑌𝑝 = 0 with sam𝑝=1

ples {𝑌𝑝 , 𝑝 = 1, ..., 𝑁 }. By incorporating the similarity
preserving term 𝐷(𝑌 ) for search accuracy and the mutual
information criterion for search time, the problem of joint
optimization for search accuracy and time together can now
be formulated as:
min 𝐼(𝑦1 , ..., 𝑦𝑘 , ..., 𝑦𝑀 )
𝐻

𝑠.𝑡.,

𝑁
∑
𝑝=1

𝑌𝑝 = 0

𝑀∑(3)

𝐷(𝑌 ) ≤ 𝜂
𝑌𝑝 = 𝐻(𝑋𝑝 )𝑘=1

3.4. Similarity Preserving Independent Component
Analysis (SPICA)

3.2. Relaxation for 𝐷(𝑌 )In sum, after relaxation with equation (5), (6), and (7),the problem in equation (3) can now be formulated as:First of all, recall that 𝑌𝑝 = 𝐻(𝑋𝑝 ) = 𝑠𝑖𝑔𝑛(𝑇 𝑇 𝑋𝑝 − 𝑏).
A traditional relaxation as in many algorithms (e.g., [11]),
is to remove the binary constraint by ignoring the 𝑠𝑖𝑔𝑛()function such that 𝐷(𝑌 ) is differentiable. In other words,
𝐷(𝑌 ) is relaxed as:𝑝,𝑞=1,...,𝑁

𝑀
∑

max

𝑇𝑘𝑇 𝐶𝑇𝑘

𝑘=1(5)
where 𝐶 = 𝑋𝐿𝑋 𝑇 . Here 𝐿 is the Laplacian matrix 𝐿 =𝑁∑
𝐷 − 𝑊 . 𝐷 is a diagonal matrix, 𝐷𝑝,𝑝 =𝑊𝑝,𝑞 .
Moreover,
1 𝑇
𝑁𝑇

𝑁
∑
𝑝=1

𝑁
∑
𝑝=1

𝑌𝑝 = 0 ⇒

𝑁
∑
𝑝=1

(6)

𝐸{𝑧𝑘 𝑧𝑗 } = 𝛿𝑘𝑗⇒ 𝐸{(𝑇𝑘𝑇 𝑥)(𝑇𝑗𝑇 𝑥)} = 𝑇𝑘𝑇 𝐸{𝑥𝑥𝑇 }𝑇𝑗 = 𝑇𝑘𝑇 Σ𝑇𝑗 = 𝛿𝑘𝑗(7)
for 1 ≤ 𝑘, 𝑗 ≤ 𝑀 .
Here 𝐶0 is a constant, 𝐸() means the expectation, 𝐺() is2some non-quadratic functions such as 𝐺(𝑢) = −𝑒−𝑢 /2 , or𝐺(𝑢) = log cosh(𝑢), etc., and 𝑔0 is a constant. 𝛿𝑘𝑗 = 1, if
𝑘 = 𝑗 ; 𝛿𝑘𝑗 = 0 , if 𝑘 ∕= 𝑗. Σ = 𝐸(𝑥𝑥𝑇 ).

and later on, we will provide a generalized version in section 3.5 to handle data with general format. Here 𝑇 is a
projection matrix of 𝑑 × 𝑀 and 𝑏 is a vector.
Even with the parameterizations of 𝐻, the problem in(3) is difficult to optimize (e.g., non-differential), and hencesome relaxation is needed. In the following discussion, weshow how to relax equation (3) with the parameterized H.

𝑊𝑝𝑞 ∣∣(𝑇 𝑇 𝑋𝑝 − 𝑏) − (𝑇 𝑇 𝑋𝑞 − 𝑏)∣∣2 =𝑁

1 ∑∣∣𝑔0 −
𝐺(𝑇𝑘𝑇 𝑋𝑝 )∣∣2
𝑁 𝑝=1

under the constraint of whiten condition, i.e.,

We first parameterize 𝐻, so that it can be optimized moreeasily. For simplicity, we first assume data are in vector
format, and 𝐻 is a linear function with a sign threshold,i.e.,𝐻(𝑥) = 𝑠𝑖𝑔𝑛(𝑇 𝑇 𝑥 − 𝑏)
(4)

∑

𝑘=1

maximizing𝑀∑∣∣𝑔0 − 𝑁1

𝑁∑(𝐺(𝑇𝑘𝑇 𝑋𝑝 ))∣∣2

𝑇𝑘 ,𝑘=1...𝑀 𝑘=1
𝑝=1
𝑠.𝑡., 𝑇𝑘𝑇 Σ𝑇𝑗 = 𝛿𝑘𝑗 , 1 ≤ 𝑘, 𝑗 ≤ 𝑀
𝑀∑
𝑇𝑘𝑇 𝐶𝑇𝑘 ≤ 𝜂
𝑘=1

(8)

where 𝐶 = 𝑋𝐿𝑋 𝑇 and Σ = 𝐸(𝑥𝑥𝑇 ). The hash bits for
𝑋𝑝 can be computed as 𝑌𝑝 = 𝑠𝑖𝑔𝑛(𝑇 𝑇 𝑋𝑝 − 𝑏), for 𝑝 =𝑁∑
𝑋𝑝 .
1, ..., 𝑁 , where 𝑏 = 𝑁1 𝑇 𝑇
𝑝=1

𝑞=1

Surprisingly, after relaxation steps mentioned above thesolution becomes quite intuitive. We call this method
SPICA (Similarity Preserving Independent Component
Analysis), because it incorporates a similarity preserving
term into the Fast-ICA formulation[19].

(𝑇 𝑇 𝑋𝑝 − 𝑏) = 0, so 𝑏 =

𝑋𝑝 .

756

3.5. Generalization–GSPICA

solve them efficiently, especially when the data set is verylarge. Inspired by the work in [19, 20], here we provide afast and efficient approximate method to solve the problemsof equation (8) and (10). The workflow of the optimization is described in algorithm 1 with details. The method
is shown to converge quite fast and performs well in theextensive experiments to be described later. The details ofderivation for the algorithm can be found in the supplement
material. Note that 𝛾 in algorithm 1 is equivalent to parameter 𝜂 in (8) and (10). Actually, 𝛾 = 0 means 𝜂 = ∞. Larger
𝛾 is equivalent to smaller 𝜂.

In practice, many applications involve structured data in
the forms of graphs, trees, sequences, sets, or other formats. For such general data types, usually certain kernel
functions are defined to compute the data similarities, e.g.,[27, 25]. Moreover, even if the data are stored in the vector format, many machine learning solutions benefit fromthe use of domain-specific kernels, for which the underlying
data embedding to the high-dimensional space is not known
explicitly, namely only the pair-wise kernel function is computable. We can obtain the kernelized version of SPICA todeal with data of general format by parameterizing the hashfunctions as:
𝑦 = 𝐻(𝑥) = 𝑠𝑖𝑔𝑛(𝑇 𝑇 𝐾𝑥 − 𝑏)4.2. Complexity and Scalability
In algorithm 1 for SPICA and GSPICA, the bottleneck of time complexity is the computation of 𝑋𝑊 𝑋 𝑇 or
𝐾𝑃 ×𝑁 𝑊 𝐾𝑃⊤×𝑁 in step 3. When we have a large scale dataset, what may consist of millions of samples, 𝑋𝑊 𝑋 𝑇 or
𝐾𝑃 ×𝑁 𝑊 𝐾𝑃⊤×𝑁 would be very expensive to compute, witha time complexity of 𝑂(𝑑𝑁 2 ) and 𝑂(𝑃 𝑁 2 ) respectively.
One way to overcome the computation complexity is touse a sparse 𝑊 . Generally speaking, one can always sample a small subset of training samples to compute similarity
matrix, so that 𝑊 is sparse, and hence the time complexityof 𝑋𝑊 𝑋 𝑇 or 𝐾𝑃 ×𝑁 𝑊 𝐾𝑃⊤×𝑁 is acceptable.
Another approach is by low rank representation/approximation for 𝑊 such that 𝑊 = 𝑅𝑄𝑅⊤ ,where 𝑅 and 𝑄 are low rank matrix. In this case,
𝐾𝑃 ×𝑁 𝑊 𝐾𝑃⊤×𝑁 can be computed as: 𝐾𝑃 ×𝑁 𝑊 𝐾𝑃⊤×𝑁 =(𝐾𝑃 ×𝑁 𝑅)𝑄(𝐾𝑃 ×𝑁 𝑅)⊤ which involves small matrices
only. There are several ways to obtain the low rank approximation, for example, we can choose 𝑊 as 𝑊 = 𝑋𝑋 ⊤ if
𝑋 are normalized data, or apply Nyström algorithm [23]to get a low rank approximation for 𝑊 . Moreover, when
the 𝑊 is defined by some ”shift-invariant” kernel functions
2
2
such as 𝑊 (𝑖, 𝑗) = 𝑒−∣∣𝑋𝑖 −𝑋𝑗 ∣∣ /𝜎 , we can apply the
kernel linearization technique [24] to approximate 𝑊 as
𝑊 = 𝑍𝑍 ⊤ , where 𝑍 are the random Fourier Features as in[24]. Note that we don’t need to compute or store 𝑊 , but
only compute or store the low rank matrix.
With the speed up, algorithm 1 would take about𝑂(𝑑2 𝑁 ) or 𝑂(𝑃 2 𝑁 ) for SPICA or GSPICA respectively,
which is close to state-of-the-art methods like spectral
hashing[11] (𝑂(𝑑2 𝑁 )) or OKH [17] (𝑂(𝑃 2 𝑁 )).

(9)where K𝑥 = [𝐾(𝑥, 𝑍1 ), ..., 𝐾(𝑥, 𝑍𝑖 ), ..., 𝐾(𝑥, 𝑍𝑃 )]𝑇 .
Here 𝐾 is the kernel function, and 𝑍𝑖 , 𝑖 = 1, ..., 𝑃 are somelandmark samples, which for example can be a subset chosen from the original 𝑁 samples. Usually 𝑃 ≪ 𝑁 .
Denote 𝑧𝑘 = 𝑇𝑘𝑇 𝐾𝑥 . With similar relaxation and derivation, we can get
max

𝑀
∑

∣∣𝑔0 − 𝑁1

𝑁∑(𝐺(𝑇𝑘𝑇 𝐾𝑋𝑝 ))∣∣2

𝑇𝑘 ,𝑘=1...𝑀 𝑘=1
𝑝=1
𝑀
∑
𝑇
𝑠.𝑡.,
𝑇𝑘 𝐶𝑇𝑘 ≤ 𝜂
𝑘=1
𝑇𝑘𝑇 Σ𝑇𝑗 = 𝛿𝑘𝑗 , 1 ≤ 𝑘, 𝑗 ≤ 𝑀

(10)

where
K𝑋𝑝 = [𝐾(𝑋𝑝 , 𝑍1 ), ..., 𝐾(𝑋𝑝 , 𝑍𝑃 )]𝑇
Moreover 𝑏 = 𝑁1 𝑇 𝑇

𝑁
∑
𝑝=1

𝐾𝑋𝑝 and 𝑌𝑝 = 𝑠𝑖𝑔𝑛(𝑇 𝑇 𝐾𝑋𝑝 −𝑏),for 𝑝 = 1, ..., 𝑁 . Here, 𝐶 = 𝐾𝑃 ×𝑁 (𝐷 − 𝑊 )𝐾𝑃𝑇 ×𝑁 .
𝐾𝑃 ×𝑁 is defined as(𝐾𝑃 ×𝑁 )𝑝,𝑖 = 𝐾(𝑍𝑝 , 𝑋𝑖 ), 𝑝 = 1, ⋅ ⋅ ⋅ , 𝑃, 𝑖 = 1, ⋅ ⋅ ⋅ , 𝑁.
(11)
𝑁
∑
𝑇
1
𝑇
And Σ = 𝐸{𝐾𝑥 𝐾𝑥 } = 𝑁
𝐾𝑋 𝑝 𝐾𝑋 𝑝 =1
𝑇
𝑁 𝐾𝑃 ×𝑁 𝐾𝑃 ×𝑁 .

𝑝=1

In equation(10), one can see that Mercer condition is unnecessary for function 𝐾. Actually, any similarity function
is applicable. We call the method in equation (10) Generalized SPICA (GSPICA), which can handle both vector
data and structured data with any kernel function or similarity/proximity function 𝐾 defined.

5. Experiments
5.1. Experiment setup
We compare our GSPICA algorithm with several stateof-the-art methods, including spectral hashing (SH) [11],locality sensitive hashing(LSH) and kernelized locality sensitive hashing (KLSH)[12].
All algorithms are compared using the same numberof hash bits. For a fair comparison, we always use the4. Optimization
4.1. Optimization algorithm
The optimization problem in both equation (8) and (10)is nonconvex. It is not trivial to obtain a fast algorithm to757

5.2. Evaluation metrics

Algorithm 1 Workflow for optimization of SPICA.
Input: data 𝑋𝑝 , 𝑝 = 1, ..., 𝑁 , similarity matrix 𝑊 , thenumber of required bits: 𝑀 .
(By replacing 𝑋𝑝 with 𝐾𝑋𝑝 and 𝑋 with 𝐾𝑃 ×𝑁 , we canobtain the optimization for GSPICA.)
Output: hash functions to generate 𝑀 bits foreach sample, i.e., 𝑌𝑝 = 𝐻(𝑋𝑝 ) = 𝑠𝑖𝑔𝑛(𝑇 𝑇 𝑋𝑝 − 𝑏)
Workflow:𝑁∑
1. Compute Σ = 𝑁1
𝑋𝑝 𝑋𝑝 𝑇 ; Apply SVD to Σ,

For evaluation, we will first compare the precision-recall
curve, which is one of the most papular evaluation methods
in large scale retrieval research.
We also report comparison of accuracy-time tradeoff fordifferent hashing methods. Search accuracy is represented
by recall rate, i.e., percentage of groundtruth neighbors
found. However, direct comparison of machine time foreach algorithm is not practical, since different implementation (e.g., different programming languages) may resultin varied search times of the same method. So in our experiments, search time is represented via the number of retrieved samples in the selected buckets. By this, we try toprovide an unbiased comparison of search time.

𝑝=1

Σ = ΩΛΩ𝑇 ;−1
2. 𝑄 = Ω𝑀 Λ𝑀2 , where Λ𝑀 is a diagonal matrix consisting of 𝑀 largest eigen values of Λ, Ω𝑀 is the corresponding column of Ω
3. Compute 𝐶 = 𝑋𝐿𝑋 𝑇 = 𝑋(𝐷 − 𝑊 )𝑋 𝑇 , Compute
C̃=𝑄𝑇 𝐶𝑄
4.
for 𝑘 = 1, ..., 𝑀 do
if 𝑘 = 1 then
𝐵=𝐼
else
apply QR decomposition to matrix [𝑇˜1 , ...𝑇˜𝑘−1 ] toget matrix 𝐵, such that [𝑇˜1 , ...𝑇˜𝑘−1 , 𝐵] is a full-rankorthogonal matrix.
end if˜ 𝑋
ˆ𝑝 = 𝐵 𝑇 𝑋
˜ 𝑝 = 𝐵 𝑇 𝑄𝑇 𝑋𝑝
𝐴 = 𝐵 𝑇 𝐶𝐵,
Random intialize w
repeat𝑁∑
ˆ𝑝 𝐺′ (𝑤𝑇 𝑋ˆ𝑝 )) − 𝛾𝑤𝑇 𝐴𝑤.
𝛽 = 𝑁1
(𝑤𝑇 𝑋
𝑝=1

𝑤+

=

𝛾𝐴]−1 [ 𝑁1

𝑤 − [ 𝑁1𝑁∑
𝑝=15.3. Experiment results
1 million web image data set
This is a data set consisting of 1𝑀 web images downloaded from flickr web site: www.flickr.com. 512 dimension gist features[22] are extracted for each image. RBF
kernel is used for this data set. 32 hash bits are used for all
the hashing methods.
To get the precision or recall, we need to obtaingroundtruth of the true nearest neighbors for each query
sample. Similar to the previous works [11, 12], we establish
such groundtruth by choosing the top samples (e.g.,top 100
samples) found via linear scan.
In Figure 1, we first show the comparison of precisionrecall curves. GSPICA performs significantly better thanother methods, confirming its superiority on search performance. Then we report accuracy-time comparison for different hashing methods. GSPICA also achieves a much better tradeoff between search accuracy and time. For instance,at the same search time (1 k retrieval samples), the recall of
our method is several times higher than other methods.
In Figure 2, some example query images and the top 5
search results are also provided. The results of the proposed
method are confirmed to be much better than others.
100K Photo Tourism image patch data set with multihash tables
One key property of LSH and its variations like KLSH
is the capacity to create multiple hash tables to improve therecall. Though only GSPICA with single hash table is discussed above, it can be easily extended to use multiple hashtables, for example, by using a subset of training set for eachtable.
We test our GSPICA method with multi hash tables on
Photo Tourism image patch set[26].
In our experiment, we use 100K patches, which are extracted from a collection of Notre Dame pictures. 1𝐾
patches are randomly chosen as queries, 90𝐾 are used astraining set to learn the hashing function. For each patch,
512 dimension gist features [22] are extracted. The task isˆ𝑝 ))𝐼 − 𝛽𝐼 −(𝐺′′ (𝑤𝑇 𝑋

𝑁∑ˆ𝑝 )) − 𝛽𝑤 − 𝛾𝐴𝑤].
ˆ𝑝 𝐺′ (𝑤𝑇 𝑋(𝑋

𝑝=1

𝑤 = 𝑤+ /∣∣𝑤+ ∣∣.
until converge
𝑇˜𝑘 = 𝐵𝑤
end for𝑁∑
5. 𝑇𝑘 = 𝑄𝑇˜𝑘 and 𝑏 = 1 𝑇 𝑇
𝑋𝑝
𝑁

𝑝=1

6. For any sample 𝑋𝑝 , compute its 𝑀 hash bits
𝑌𝑝 = 𝑠𝑖𝑔𝑛(𝑇 𝑇 𝑋𝑝 − 𝑏)same kernel function with the same parameters (if any),
when kernels are used in methods including GSPICA and
KLSH. The same number of landmark samples are used for
GSPICA and KLSH. And moreover, the number of landmark points are set close to the number of feature dimensions, so that GSPICA and SH would have almost the sameindexing time. The parameter 𝛾 is chosen with a validation
set, which is independent of the training set or the query set.
758

0

10

GSPICA with Linear Kernel
GSPICA with RBF Kernel
Spectral hashing
KLSH with Linear Kernel
KLSH with RBF Kernel
LSH−1Precision

10

−2

10

−3

10

−4

10

0

0.2

0.4

0.6

0.8

1

Recall(a) Precision-recall curve on 1M web image data
1

Recall0.8
0.6

GSPICA with Linear Kernel
GSPICA with RBF Kernel
Spectral hashing
KLSH with Linear Kernel
KLSH with RBF Kernel
LSH

0.4
0.2
0 −2
10

0

2

4

10
10
10
Number of retrieved samples

610(b) Accuracy-time comparison on 1M web image data
Figure 1. Search results on 1M web image data set with 32 hashbits. In (a), the comparison of precision-recall curve is provided.
In (b), comparison of accuracy-time curve is shown where recall
represents search accuracy, and the number of retrieved samples inselected buckets represents search time. Graphs are best viewed incolor.

(a) Queries(b) Top 5 search results.

Figure 2. On 1M web image data set, example query images andtop 5 search results of GSPICA, SH, KLSH, and LSH ranked byhamming distance with 32 hash bits are shown. Note that we areusing gist features, so color information is not considered.

to identify the neighbors, i.e., near-duplicate patches, in the
training set for each query patch. For each patch, its nearduplicates are used as the groundtruth. In our experiments,
we randomly sample 10,000 training samples to compute48 bits for each hash table. The same procedure are alsodone to create multi tables for spectral hashing.
The results are shown in Figure 3. By using multihash tables, the recall of our GSPICA method is improved.
Moreover, GSPICA works significantly better than LSH,
KLSH or spectral hashing, no matter with a single hash table or multi hash tables.
We also explores how the change of experiment setting,e.g., the number of landmark samples P, or the parameter 𝛾
in algorithm 1, would affect our results. As shown in Figure
4, the proposed method is quite insensitive and stable toreasonable change of 𝑃 and 𝛾. And not surprisingly, theperformance increases slowly with 𝑃 .

while keeping the data similarity at the same time, which
is highly related to problems like dimension reduction, subspace learning, etc. We will explore application of SPICA
for such tasks and other related domains in our future work.
Acknowledgments”This work has been supported in part by NSF awards#CNS-07-51078 and #CNS-07-16203.”

References[1] J. Freidman, J. Bentley, and A. Finkel. An Algorithm for Finding BestMatches in Logarithmic Expected Time. ACM Transactions on Mathematical Software, page 209-226, 1977.

6. Conclusion[2] J. Uhlmann. Satisfying general proximity / similarity queries
with metric trees. Information Processing Letters, page 175179, 1991.

In this paper, we develop a new hashing algorithm withjoint optimization of search accuracy and search time together, and support data of general formats Our methods
significantly outperforms several state-of-the-art hashing algorithms over diverse test set and data types.
Our proposed technique, SPICA and GSPICA, try tofind independent (and hence least redundant) projections[3] P. Indyk and R. Motwani. Approximate nearest neighbors:
towards removing the curse of dimensionality. In Proceedingsof STOC, 1998.
[4] M. Charikar. Similarity estimation techniques from rounding
algorithms. In Proceedings of STOC, 2002.

759

SH
1 tableGSPICA 1 table
GSPICA 5 tables
KLSH
5 tablesKLSH
1 tableSH
5 tablesLSH
1 tableLSH
5 tables1

Precision

0.8
0.6

learned metrics. In Proceedings of CVPR, 2008.
[8] R. Salakhutdinov and G. Hinton. Semantic hashing. In Proceedings of SIGIR, 2007.
[9] R. Salakhutdinov and G. Hinton. Learning a nonlinear embedding by preserving class neighbourhood structure. In Proceedings of AI and Statistics, 2007.
[10] A. Torralba, R. Fergus, and Y.Weiss. Small Codes and Large
Image Databases for Recognition. In Proceedings of CVPR,
2008.
[11] Y. Weiss, A. Torralba, and R. Fergus. Spectral hashing. In
Proceedings of NIPS, 2008.
[12] Brian Kulis and Kristen Grauman. Kernelized LocalitySensitive Hashing for Scalable Image Search. In Proceedingsof ICCV, 2009.
[13] Brian Kulis and Trevor Darrell. Learning to hash with binary
reconstructive embeddings. In Proceedings of NIPS, 2009.
[14] M. Raginsky and S. Lazebnik. Locality sensitive binary
codes from shift-invariant kernels. In Proceedings of NIPS,
2009.
[15] Ruei-Sung Lin David A. Ross Jay Yagnik. SPEC Hashing:
Similarity Preserving algorithm for Entropy-based Coding In
Proceedings of CVPR, 2010.
[16] Yadong Mu, Jialie Shen, Shuicheng Yan. Weakly-Supervised
Hashing in Kernel Space In Proceedings of CVPR, 2010.
[17] Junfeng He, Wei Liu and Shih-Fu Chang. Scalable Similarity Search with Optimized Kernel Hashing In Proceedings of
KDD, 2010.
[18] Shumeet Baluja and Michele Covell. Learning to hash: forgiving hash functions and applications. Data Mining andKnowledge Discovery, 2008.
[19] Aapo Hyvarinen and Erkki Oja. Independent Component
Analysis: Algorithms and Applications. Neural Networks,
page 411-430, 2000.
[20] Alper T. Erdogan. On the Convergence of ICA Algorithms
With Symmetric Orthogonalization. IEEE Transactions on
Signal Processing, 2009.
[21] Xiaofei He, and Partha Niyogi. Locality Preserving Projections. In Proceedings of NIPS, 2003.
[22] A. Oliva and A. Torralba. Modeling the shape of the scene:a holistic representation of the spatial envelope. International
Journal on Computer Vision, 2001.
[23] Christopher Williams , Matthias Seeger. Using the nystrom
method to speed up kernel machines. In Proceedings of NIPS,
2001.
[24] Ali Rahimi and Benjamin Recht. Random Features forLarge-Scale Kernel Machines. In Proceedings of NIPS, 2007.
[25] S. Lazebnik, C. Schmid and J. Ponce. Beyond bags of features, spatial pyramid matching for recognizing natural scene
categories. In Proceedings of CVPR, 2006.
[26] S. Winder and M. Brown. Learning Local Image Descriptors.
In Proceedings of Computer Vision and Pattern Recognition,
2007.
[27] Nino Shervashidze, Karsten M. Borgwardt. Fast subtree kernels on graphs. In Proceedings of NIPS, 2009.

0.4
0.2
0
0

0.2

0.4

0.6

0.8

1

Recall(a) Precision-recall curve on Photo Tourism data
1

Recall

0.8
SH
1 tableGSPICA 1 table
GSPICA 5 tables
KLSH
5 tablesKLSH
1 tableSH
5 tablesLSH
1 tableLSH
5 tables

0.6
0.4
0.2
0 −2
10

0

2

4

10
10
10
Average number of retrieved samples

610(b) Accuracy-time comparison on Photo Tourism data
Figure 3. Search results on 100K Photo Tourism image patch datawith 48 hash bits for each hash table. In (a), the comparison of
precision-recall curve is provided. In (b), comparison of accuracytime curve is shown where recall represents search accuracy, andthe number of samples in selected buckets represents search time.
Graphs are best viewed in color.
1γ=10, P = 300
γ=1, P = 200
γ=1, P = 400
γ=0.1, P = 300
γ=1, P = 300

Precision

0.8
0.6
0.4
0.2
0
0

0.2

0.4

0.6

0.8

1

Recall(a) Precision-recall curve for different experiment setting
Figure 4. Search results on 100K Photo Tourism image patch datawith 48 hash bits for one hash table with different number of landmark samples 𝑃 and parameter 𝛾 in algorithm 1. Note that thegreen curve with square marks are covered by other curves andcan not be seen. As shown, the proposed algorithm is quite stable and not sensitive to reasonable change of 𝑃 and parameter 𝛾.
Graphs are best viewed in color.

[5] M. Datar, N. Immorlica, P. Indyk, and V. Mirrokni. Locality sensitive hashing scheme based on p-Stable distributions.
In Proceedings of Symposium on Computational Geometry(SOCG), 2004.
[6] K. Grauman and T. Darrell. Pyramid match hashing: subLinear time indexing over partial correspondences. In Proceedings of CVPR, 2007.
[7] P. Jain, B. Kulis, and K. Grauman. Fast image search for760