To Project More or to Quantize More: Minimize Reconstruction Bias for Learning Compact Binary Codes

1 Introduction

Approximate nearest neighbour (ANN) search [Gionis et al., 1999a] plays an important role in machine learning, computer vision and information retrieval. Using similarity preserving binary codes to represent original data points is of particular interest for ANN search [Weiss et al., 2008; Norouzi and Fleet, 2011]. The binary codes can bring about low memory cost as well as fast similarity distance computing speed. This is particular useful when dealing with large scale database [Torralba et al., 2008; Gong and Lazebnik, 2011; Weiss et al., 2008; Duan et al., 2016].

A common binary coding approach, often called Hashing, is to develop similarity preserving hashing functions for mapping data points into a Hamming space. As it is NPhard to directly learn the optimal binary codes [Weiss et al., 2008], hashing methods typically work on a two-stage strategy: projection and quantization [Kong and L, 2012; Kong et al., 2012]. Specifically, given a data point x 2 Rd , they first project x into a low dimensional vector

where real-valued functions ¶fi(·)♦ki=1 are called projection functions. Then they utilize Single Bit Quantization (SBQ) to quantize each projection element fi(x) into a single bit by thresholding [Kong and L, 2012; Wang et al., 2016].

Lots of research efforts have been devoted to the first stage, with an aim to learn powerful projections to maintain the similarity structure of the original data points. Local Sensitive Hashing (LSH) [Andoni and Indyk, 2006] adopts a random projection which is independent of training data. Similarly, Shift Invariant Kernel Hashing (SIKH) [Raginsky and Lazebnik, 2009] chooses random projection and applies shifted cosine function to generate binary codes. Both LSH and SIKH are data independent and flexible since they do not rely on any training data. However, long codes are often required to achieve satisfactory performance [Gong and Lazebnik, 2011; Raginsky and Lazebnik, 2009].

To build up more effective projection, many promising data dependent methods have been proposed. Through learning the projection functions over training data, data dependent methods usually outperform data independent methods at relatively shorter codes [Liu et al., 2010]. Representative methods include Spectral Hashing [Weiss et al., 2008], Binary Reconstructive Embedding Hashing [Kulis and Darrell, 2009], Semi-Supervised Hashing [Wang et al., 2010], Anchor Graph Hashing [Liu et al., 2010], Iterative Quantization [Gong and Lazebnik, 2011], Minimal Loss Hashing [Norouzi and Fleet, 2011], Kernal Supervised Hashing [Liu et al., 2012], Isotropic Hashing [Kong and Li, 2012], K-means Hashing [He et al., 2013], Inductive Hashing on Manifolds [Shen et al., 2013], Harmonious Hashing [Xu et al., 2013], Discrete Graph Hashing [Liu et al., 2014], Sparse Projection Hashing [Xia et al., 2015], etc.

Moreover, recent works have reported the impact of quantization on Hashing performance. Single Bit Quantization (SBQ) in most hashing methods incurs lots of quantization errors, which would seriously degrade the performance [Kong and L, 2012; Kong et al., 2012]. Thus, promising Multiple Bits Quantization (MBQ) methods have been proposed. Double Bits Quantization [Kong and L, 2012] divides each projection dimension into three regions and uses double bits code to represent each element region. Manhattan Quantization [Kong et al., 2012] proposes natural binary code (NBC) and adopts Manhattan distance to compute the distance between NBC codes. Hamming Compatible Quantization [Wang et al., 2015] aims to minimize the distance error function to preserve the capability of similarity metric between Euclidean space and Hamming space. Overall, MBQ methods do facilitate the reduction of information loss in quantization. Experiment results have demonstrated the functionality of high quality quantization in improving Hashing performance [Kong and L, 2012; Kong et al., 2012; Wang et al., 2015].

Hence, both projection and quantization are important. Optimal binary codes rely on the joint optimization of projection and quantization. However, how to elegantly connect optimizing projection with optimizing quantization, and to maximize the complementary effects of two stages, remains a challenging problem in practice. Given a specified binary code length k, shall we project more or quantize more?

The projection dimensionality in existing hashing methods or quantization methods are all fixed, which is not flexible. An intuitive thought is, when the original data points inherently lie in a low dimensional space, we may project data points to a lower-dimensional space (project more), while performing finer quantization for each element by using multiple bits instead of a single bit (quantization less). For example, we may project original data points into space R k2 and assign 2 bits to quantize the values of each projection element, while the target code length remains as k. Through adaptively adjusting the projection dimensionality, we may minimize the information loss over the whole course of Hashing, in which the balance between projection and quantization can be addressed by the optimal projection dimensionality.

In this paper, we propose a novel approach Minimal Reconstruction Bias Hashing (MRH) to tackle the joint optimization of projection and quantization. We summarize the main contributions of this paper as follows:

• We present a novel approach to learn similarity preserving binary codes which jointly optimizes both projection and quantization stages with adjustable projection dimensionality. To the best of our knowledge, this is the first work that can adaptively adjust the projection dimensionality to balance the information loss between projection and quantization. Our practice of jointly optimizing projection dimensionality, projection matrix, as well as quantization functions, has achieved the state-ofthe-art performance over several benchmarks.
• We reinterpret the problem of maximizing the similarity preserving in Hashing from the perspective of minimal reconstruction bias of signals. By introducing a lower bound analysis, we establish the relationship between the information loss from projection and quantization, and the Hamming approximation errors, which justifies the learning objective of our proposed MRH method.
• By analyzing the unimodal characteristics of the MRH objective function with respect to projection dimensionality, we propose an effective solution to resolve the joint optimization problem of MRH. In particular, we have reduced the complexity of searching optimal projection dimensionality from O(N) to O(log(N)).

2 Preliminaries

We firstly introduce the basic notations. Let matrix X = [x1↪ x2↪ · · · ↪ xn] 2 Rd⇥n denote the samples of data points, and k denote the length of target codes. The goal is to learn a binary string bi 2 ¶0↪ 1♦k for each data point xi 2 Rd with the maximized similarity preservation in Hamming space.

The proposed binary code learning approach involves both projection and quantization stages. At the projection stage, we adopt linear projection to transform xi 2 Rd into a subspace

where T(xi) = Rxi and R 2 R kc ⇥d . The projection matrix R is required to be orthogonal, i.e., RR> = I, to make the projections of a vector in different dimensions independent of each other. At the quantization stage, we quantize projection vector yi into

byi = Q(yi) 2 H kc ↪ where H is the set of quantization centroid(s), and each element in yi is quantized to a value in H. Finally, we use c bits to encode each element of byi to obtain a binary string bi = B(byi) 2 ¶0↪ 1♦k ▷ Note that we introduce variable c to adjust the projection dimensionality. If the given code length k is indivisible by c, the target code length will round down to b kc c⇥c bits. In this paper, we will figure out an optimal c value to make a joint optimization of projection T(·) and quantization Q(·).

The range of Hamming distance is limited to the length of a binary code. The maximum Hamming distance of c bits codes is only c. When we use c bits to encode 2c values, the distance consistency in the Hamming space cannot maintained. Let us take the example of c = 2. We quantize the projection values into 22 = 4 centroid(s) with σ0 < σ1 < σ2 < σ3, namely, (00)2, (01)2, (10)2 and (11)2. We have kσ1-σ2k < kσ1-σ3k, but dH(01↪ 10) > dH(01↪ 11) where dH(▷) denotes the Hamming distance.

To address the issue of inconsistent measurements, we may resort to other distance measurer like Manhattan distance [Kong et al., 2012]. But this would seriously degrade the retrieval efficiency [Wang et al., 2015]. By contrast, Hamming distance measurement is extremely fast, and more than 109 operations can be done per second [He et al., 2013; Weiss et al., 2008], so that Hamming distance computing is still the priority of effective and efficient ANN search. In this work, we employ an incomplete encoding strategy to keep the distance consistency in Hamming space, in which we only quantize projection values into c + 1 equidistant centroid(s) H = ¶σi♦ci=0 , where σi - σi-1 = ∆ for 1  i  c. Then, we apply unary representation [Gionis et al., 1999b] to encode each centroid B(σi) = Uc(i) for σi 2 H, where unary representation Uc(i) is defined as a c bits binary string with i ones followed by c - i zeros, e.g. U2(1) = 10, U3(0) = 000, U4(2) = 1100. With the unary representation, the Hamming distance between binary codes is proportional to the distance of the centroid(s),

Clearly, unary representation is an incomplete encoding method, as a c bits code can represent 2c states. To make the Hamming distance consistent with the distance of quantization centroid(s), we have to discard parts of code space.

3 Minimal Reconstruction Bias Hashing

We formulate the problem of optimal binary coding (i.e., optimal hashing) from the perspective of minimizing the reconstruction bias of signals [Allen and Gray, 2012]. The relationship between minimal reconstruction bias and Hamming approximation errors will be studied as well.

3.1 Reconstruction Bias

The reconstructed data points are recovered from the compressed codes. Given the hashing code bi of data point xi, we obtain the reconstructed data by first decoding bi to quantization centroid(s) and then transforming the quantization vector back to the original space (see Fig.1). Specifically, we first decode bi and get byi = B-1 (bi). Quantization function Q(▷) is not invertible, so yi can't be recovered. Then, we directly apply the inverse projection transformation to byi and get bxi = T-1 (B-1 (bi)) = R> byi▷ (1) Here, bxi is called the reconstructed data of xi. The recon-

d(xi↪ bxi) = kxi - bxik2 = kxi - R> byik2 where d(·) denotes the Euclidean distance, k · k2 denotes the L2 norm of a vector.

3.2 Learning Objective

The reconstruction bias indicates the information loss incurred by mapping the data points into Hamming space. To preserve the similarity structure of original data points, we aim to minimize the reconstruction bias. Directly optimizing the objective function in Eqn.2 is intractable due to a large number of free parameters in byi and the orthogonal constraint

kxi -R> byik2  kxi -R>yik2 +kR>(yi - byi)k2▷ (3) The first term indicates the distortions by transforming yi back to xi. When R is orthogonal, for any vector a, we have 1 kR>ak2 = kak2. So the second term actually indicates the mean square error (MSE) of byi. Let Y = [y1↪ y2↪ · · · ↪ yn] denotes the projection matrix and bY = Q(Y). To resolve optimal binary coding, we formulate the problem of jointly minimizing the projection distortions and the quantization errors, in which the projection dimensionality may be variable

as well. Specifically, the learning objective is formulated as

Y =RX↪ bY 2 H kc ⇥n ↪ kHk = c +1▷ where k · kF denotes the Frobenius norm, kX-R>Yk2F denotes the sum of projection distortions, indicating the information loss in the projection stage, and kY - bYk2F denotes the sum of mean square error (MSE), indicating the information loss in the quantization stage. In particular, variable c is to adjust the projection dimensionality to adaptively balance the information loss between the projection and quantization stages in a joint optimization.

3.3 Relationship to Hamming Approximation

Similarity preserving hashing methods aim to map close data points to near binary codes [Gionis et al., 1999a; Andoni and Indyk, 2006]. Conversely, if two data points are far away in the original space, their binary codes should produce a large Hamming distance. We will show that the distance approximation error between the original distance and the root mean square Hamming distance is a lower bound of the learning objective in Eqn.4. Since the Hamming approximation quality, as a critical indicator, can significantly impact the performance of ANN search [Kulis and Darrell, 2009; He et al., 2013], this lower bound analysis may justify the rationale of the proposed learning objective.

In hashing methods, the similarity of two data points xi and xj is defined by the Hamming distance of their hashing codes, dH(bki ↪ bkj ), where bki = B(byk i ), bki and byk i denote the k-th element in vector bi and byi, respectively. Let l denote the projection dimensionality, sk the Hamming distance of the k-th hashing codes where sk = dH(bki ↪ bkj ). The root mean square Hamming distance f(bi↪ bj) of two binary strings bi and bj is defined as f(bi↪ bj) = (Pli=1 s2i ◁l) 12 . Consider the distance approximation error between the original distance d(xi↪ xj) and the root mean squared Hamming distance f(bi↪ bj). We have the theorem

Theorem 1. The distance approximation error between the original distance and the root mean squared Hamming distance is a lower bound of objective function G,

where parameter λ and θ are two constant factors, which λ = ∆pl and θ = 32n.

Proof. According to the triangle inequality, we have

Relax the right side of the above inequality, according to CauchySchwarz Inequality, we have Pi d(xi↪ bxi)

By substituting d(bxi↪ bxj) = ∆plf(bi↪ bj) to (a1) and squaring both sides of the inequality, we obtain Theorem 1.

4 Optimization

The goal is to minimize the objective of overall information loss in Eqn.4 with variables c, R and bY.

Update bY When updating bY, the learning objective reduces to the second term kY - bYk2F . We simply assume that the data distribution is zero-centered. If c is an odd number where c = 2t + 1, the quantization centroid set H is represented as H = ¶0↪◦∆↪ ▷▷▷↪◦t∆♦. Without loss of generality, we just show the odd case. We quantize each element y 2 Y to the nearest value in HQ(y) = arg min

4.1 Update bY and R To resolve R and bY, we first fix the variable c, which is meant to fix the projection dimensionality over the course of alternating optimization of R and bY. The alternating fashion works by updating R or bY with the other fixed.

The only variable in quantization function Q(▷) is ∆. We

The above optimization problem is a multi chromatic pair problem in computational geometry [Berg et al., 2000], which can be solved by Expectation Maximization (EM) algorithm [Moon and K, 1996].

Update R

Updating R is a typical orthogonality constraint optimization problem. We apply the optimization procedure in [Z. and W., 2013] to update R. Let U be the partial derivative of the objective function with respect to R. We have U = -2 bYX>. We first define the skew-symmetric matrix [Z. and W., 2013]

Then, we adopt Crank Nicolson like [Smith, 1965] scheme to update the orthogonal matrix R

where ⌧ denotes the step size, we empirically set ⌧ = 0▷5.

Convergence

We alternatively update bY and R in several iterations until convergence. In practice, we have found the algorithm converges in about 50-100 iterations. A typical behavior of the information loss (4) is shown in Fig.5.

4.2 Update c to Balance the Information Loss

Given a specified target code length, setting c as a large value can improve quantization quality but would degrade projection quality, and vice versa. There exists a trade-off between projection and quantization. To balance the information loss between projection and quantization, we aim to optimize c to minimize the objective of overall reconstruction bias. The value of c ranges from 1 to the target code length k (say hundreds or thousands). The brute-force enumeration method would be costly.

Rather than exhaustive search, we propose a fast approach to search for the optimal c, as our empirical findings have shown that the objective function with respect to c is unimodal. To explain this important findings, we have derived the theorem 2 by assuming a moderate distribution function (i.e. uniform or Gaussian) of projection values.

Theorem 2. Function G(c) can be well-approximated by an unimodal function, which only has a single local minimum point c⇤ . G(c) is monotonically decreasing for c  c⇤ and monotonically increasing for c > c⇤ .

Proof. According to the orthogonal constraint in R, we have

Let matrix E = Y - ˆY. G(c) can be simplified as

The first term is independent of c. Each term of kYk2F or kEk2F contains n⇥ b dc c elements. The expression of G(c) depends on the distribution of projection values. We adopt statistical expectation for sample estimation. Without loss of generalization, assuming d is divisible by c, we have

When projection values Y are subject to an uniform distribution, the probability density function is given by

Accordingly, we have [Allen and Gray, 2012]

In our method, step size ∆ = p2-p1 c+1 . Then,

and

Let λ = nd(p21 + p22 + p21p22 )◁3, θ = nd(p2 - p1)2 ◁12 and ⌘ = kXk2F , G(c) can be represented as

Take the derivative of G with respect to c. We have

and

Let H(c) denote the function of right side in above equation,

H(c) is monotonically increasing for c ≥ 1. Assume c⇤ is the minimal point where G0 (c⇤ ) = 0 and H(c) = λθ . Then, for c > c⇤ we have H(c) > λθ i.e. G0 (c⇤ ) > 0, and for c < c⇤ we have H(c) < λθ i.e. G0 (c⇤ ) < 0. Thus, G(c) is unimodal.

Likewise, for Gaussian distribution f(y) = 1 σp2⇡ e- (y-θ)2 2δ2 , y 2 [θ-p↪ θ+p], we can still derive the expression of G by computing the expectation of y2 and (y - (Q(y)))2 . We employ the second order Taylors expansion to represent f(y) and calculate E(y2 ) and E((y - (Q(y)))2 ) by computing the integral of the Taylors expansion, followed by the derivative and monotonic analysis for the proof.

G(c) is defined in discrete domain c 2 ¶1↪ 2↪ ▷▷▷↪ k♦. According to the unimodal property in Theorem.2, we can apply ternary search to find out the optimal c⇤ . We have Theorem 3. Let c⇤ denote the minimum point of objective function G(c). Assume we have already known c⇤ 2 [l↪ r]. Let s = r-l 3 , m1 = bl + 13 sc and m2 = bl + 23 sc. We have

• if G(m1) > G(m2), then m1  c⇤  r.
• if G(m1)  G(m2), then l  c⇤ m2.

Proof. Consider m1 6= m2. For G(m1) < G(m2), we have c⇤  m2. If not, then m1 < m2 < c⇤ . As function G(c) is monotonically decreasing for c < c⇤ , then G(m1) > G(m2). This is contradictory with the assumption. Thus, c⇤  m2 and l  c⇤  m2. For G(m1) > G(m2), the same procedure can be adapted to obtain m1  c⇤  r. If m1 = m2, we have l = r. Obviously, Theorem 3 still holds.

The ternary search algorithm works as follows. We initialize l = 1 and r = k. For each iteration, we set c = m1 and c = m2 respectively, and solve G(m1) and G(m2) by alternatively updating R and bY. If G(m1)  G(m2), we update r = m2. Otherwise, we update l = m1. The algorithm terminates when l = r. As we cut out 1◁3 search scope after each iteration, the run time order is

The reduced complexity benefits the fast search of optimal c, especially when learning long binary codes. Figure 2 shows an example of ternary search.

Algorithm 1 The algorithm to optimize Equation (4).

Input: original data points ¶xi♦li=1 and target code length k. Output: a binary string bi for each data point xi.

Initialize l = 1 and r = k.

while l < r do

Let s = r-l

Otherwise, if G(m1) > G(m2), set l =m1.

3 , m1 = bl + 13 sc and m2 = bl + 23 sc. Set c = m1, c = m2, and calculate G(m1), G(m2) by alternatively updating R and bY. If G(m1)  G(m2), then set r = m2.

end while

Set c = l. Project each xi into yi 2 R kc and quantize each projection value into c bits to obtain a binary string bi with b kc c ⇥ c bits.

1

1.5

2

4.3 Complexity Analysis

We need O(log k) recursions to find out the optimal c. In each recursion, we iteratively update bY and R. Let t denote the number of iterations for the alternatively updating, t1 the iteration number in EM algorithm and t2 the iteration number in Crank Nicolson like scheme [Smith, 1965]. It takes O(t1c2n) to update bY and O(nl + d2 + t2ld) to update R. The overall time complexity is O(log k(t1c2n + ln + d2 + t2ld)). We will show the running time cost in the next section. Algorithm 1 shows the pseudo-code of our MRH algorithm.

5 Experiments

We evaluate and compare the approaches over five benchmarks SIFT1M [Jegou et al., 2011], GIST1M [Jegou et al., 2011], CIFAR10 [Krizhevsky, 2009], LableMe22k [Torralba et al., 2008] and ImageNet1M [Deng et al., 2009]. SIFT1M [Jegou et al., 2011] and GIST1M [Jegou et al., 2011] are the feature datasets containing 1 million 128-D SIFT [Lowe, 2004] and 960-D GIST [Aude and Torralba, 2001] features, respectively. The CIFAR10 dataset contains 60,000 images. The LabelMe22K dataset contains 22,019 images. Each image in CIFAR-10 and LabelMe22K is represented by a 512 dimensional GIST feature [Aude and Torralba, 2001]. ImageNet1M is another large-scale benchmark with 1 million images. For each image, we extract a 4096D fisher vector [Perronnin and Dance, 2006] to evaluate the performance in a high dimensional space.

5.1 Baseline Methods

We perform extensive comparison with 7 state-of-the-art methods: Local sensitive hashing (LSH) [Andoni and Indyk, 2006], Iterative quantization (ITQ) [Gong and Lazebnik, 2011], Shift invariant kernels hashing (SIKH) [Raginsky and Lazebnik, 2009], Principal component analysis hashing (PCAH) [Wang et al., 2006], Spectral hashing (SH) [Weiss et al., 2008], K-means hashing (KMH) [He et al., 2013] and Sparse Projection Hashing (SP) [Xia et al., 2015]. All the methods are run with released source codes in default settings.

5.2 Setting and Configuration

We follow most of previous hashing works to adopt the Hamming distance ranking for ANN search. Recall rate and mean average precision (mAP) are used for performance comparison. For each benchmark, we select 1000 data points as queries and leave the rest as database. For each query, the top 100 nearest data points in Euclidean distance are used as the ground truth. We random select 10K data points from each dataset for training. We empirically set t = 100, t1 = 50 and t2 = 50. Figure 5(b) shows the training time. Except for KMH [He et al., 2013], other baseline methods are run in a single thread. All experiments are carried out over a server with Intel i5 CPU at 3.20GHz and 64Gb memory cache.

5.3 Performance

Figure 4 and Figure 3 shows the results of recall rate of baseline methods over datasets GIST1M and SIFT1M. The proposed MRH consistently outperforms the state-of-the-art methods. Over dataset GIST1M, the performance gains of MRH are with 9.1%, 13.4% and 17.8% recall rate when ranking 20,000 in code length 32 bits, 64 bits and 128 bits, respectively, and with 7.1%, 5.4% and 3.2% recall rate on SIFT1M.

Figure 6 shows the mAP results of baseline methods. We report the results over LabelMe, CIFAR and ImageNet1M. MRH outperforms the baselines on all settings even at small codes. As the code length increases, the performance gap becomes more significant. For example, MRH outperforms the competitive method SP by 2.1%, 2.9%, 4.2%, 7.4%, 8.5% and 7.1% at code length 8 bits to 256 bits over LabelMe, respectively. Considerable improvements are also obtained over CIFAR10 as well as ImageNet1M.

Figure 7 shows the impact of projection dimensionality for learning 512 and 1024 bits codes over ImageNet1M. Variable c produces balancing effects on projection and quantization.

Increasing c reduces the quantization error but it incurs more projection distortions, and vice versa. There exists a trade-off between projection and quantization. From figure 7, the best setting is c = 4 for 512 bits and c = 6 for 1024 bits. We notice that MRH tends to set c to a large value for long size codes, which means that more bits are allowed for quantizing the values of each projection element. By adaptively adjusting the projection dimensionality, MRH obtains discriminative codes with overall minimal information loss. Table 1 lists the optimal settings of c on three datasets.

6 Conclusion

The joint optimization of projection and quantization impacts the similarity preserving of binary codes, which is important for generating discriminative binary codes. To optimize the projection dimensionality has tackled the problem of balancing the information loss between the projection and quantization stages. The practice of jointly optimizing projection dimensionality, projection matrix, as well as quantization functions is expected to facilitate the state-of-the-art Hashing methods.

Acknowledgments

This work was supported by the National Hightech R&D Program of China (863 Program): 2015AA016302, and Chinese Natural Science Foundation: 61271311, 61390515, 61421062. Part of this work was carried out at the RapidRich Object Search (ROSE) Lab at the Nanyang Technological University, Singapore. The ROSE Lab is supported by the National Research Foundation, Singapore, under its Interactive Digital Media (IDM) Strategic Research Programme. Ling-Yu Duan is the corresponding author.