Latent Structure Preserving Hashing

Aiming at efficient similarity search, hash functions are designed to embed high-dimensional feature descriptors to low-dimensional binary codes such that similar descriptors will lead to binary codes with a short distance in the Hamming space. It is critical to effectively maintain the intrinsic structure and preserve the original information of data in a hashing algorithm. In this paper, we propose a novel hashing algorithm called Latent Structure Preserving Hashing (LSPH), with the target of finding a well-structured low-dimensional data representation from the original high-dimensional data through a novel objective function based on Nonnegative Matrix Factorization (NMF) with their corresponding Kullback-Leibler divergence of data distribution as the regularization term. Via exploiting the joint probabilistic distribution of data, LSPH can automatically learn the latent information and successfully preserve the structure of high-dimensional data. To further achieve robust performance with complex and nonlinear data, in this paper, we also contribute a more generalized multi-layer LSPH (ML-LSPH) framework, in which hierarchical representations can be effectively learned by a multiplicative up-propagation algorithm. Once obtaining the latent representations, the hash functions can be easily acquired through multi-variable logistic regression. Experimental results on three large-scale retrieval datasets, i.e., SIFT 1M, GIST 1M and 500 K TinyImage, show that ML-LSPH can achieve better performance than the single-layer LSPH and both of them outperform existing hashing techniques on large-scale data.


Introduction
Similarity search Gionis et al. 1999;Qin et al. 2015;Yu et al. 2015;Gao et al. 2015;Zhang et al. 2010;Wang et al. 2014;Bian and Tao 2010) is one of the most critical problems in information retrieval as well as in pattern recognition, data mining and machine learning. Generally speaking, effective similarity search approaches try to construct the index structure in the metric space. However, with the increase of the dimensionality of the data, how to implement the similarity search efficiently and effectively has become an significant topic. To improve retrieval efficiency, hashing algorithms are deployed to find a hash function from Euclidean space to Hamming space. The hashing algorithms with binary coding techniques mainly have two advantages: (1) binary hash codes save storage space; (2) it is efficient to compute the Hamming distance (X O R operation) between the training data and the new coming data in the retrieval procedure of similarity search. The time complexity of searching the hashing table is near O(1).
Current hashing algorithms can be roughly divided into two groups: random projection based hashing and learning based hashing. For the random projection based hashing techniques, the most well-known hashing technique that preserves similarity information is probably Locality-Sensitive Hashing (LSH) (Gionis et al. 1999). LSH simply employs random linear projections (followed by random thresholding) to map data points close in a Euclidean space to similar codes. It is theoretically guaranteed that as the code length increases, the Hamming distance between two codes will asymptotically approach the Euclidean distance between their corresponding data points. Furthermore, kernelized localitysensitive hashing (KLSH) (Kulis and Grauman 2009) has also been successfully proposed and utilized for large-scale image retrieval and classification. However, in realistic applications, LSH-related methods usually require long codes to achieve good precision, which result in low recall since the collision probability that two codes fall into the same hash bucket decreases exponentially as the code length increases.
However, the random projection based hash functions are effective only when the binary hash code is long enough. To generate more effective but compact hash codes, a number of methods such as projection learning for hashing have been introduced. Through mining the structure of the data, then being represented on the objective function, a projection learning based hashing algorithm can obtain the hash function by solving an optimization problem associated with the objective function. Spectral Hashing (SpH) (Weiss et al. 2009) is a representative unsupervised hashing method, which can learn compact binary codes that preserve the similarity between documents by forcing the balanced and uncorrelated constraints into the learned codes. Furthermore, principled linear projections like PCA Hashing (PCAH) ) has been suggested for better quantization rather than random projection hashing. Moreover, Semantic Hashing (SH), which is based a stack of Restricted Boltzmann Machines (RBM) (Salakhutdinov and Hinton 2007), was proposed in Salakhutdinov and Hinton (2009). In particular, SH involves two steps: pre-training and fine-tuning. During these two steps, a deep generative model is greedily learned, in which the lowest layer represents the high-dimensional data vector and the highest layer represents the learned binary code for that data. Liu et al. (2011) proposed an Anchor Graph-based Hashing method (AGH), which automatically discovers the neighborhood structure inherent in the data to learn appropriate compact codes. To further make such an approach computationally feasible, the Anchor Graphs used in Liu et al. (2011) were defined with tractable low-rank adjacency matrices. In this way, AGH can allow constant time hashing of a new data point by extrapolating graph Laplacian eigenvectors to eigenfunctions. More recently, Spherical Hashing (SpherH) (Heo et al. 2012) was proposed to map more spatially coherent data points into a binary code compared to hyperplane-based hashing functions. Meanwhile, the authors also developed a new distance function for binary codes, spherical Hamming distance, to achieve final retrieval tasks. Iterative Quanti-zation (ITQ) (Gong et al. 2013) was developed for more compact and effective binary coding. Particularly, a simple and efficient alternating minimization scheme for finding a orthogonal rotation of zero-centered data so as to minimize the quantization error of mapping this data and the vertices of a zero-centered binary hypercube. Additionally, Boosted Similarity Sensitive Coding (BSSC) (Shakhnarovich 2005) was designed to learn a compact and weighted Hamming embedding for task specific similarity search. Boosted binary regression stumps were used as hashing functions to map the input vectors into binary codes. A similar idea as BSSC is also applied to Evolutionary Compact Embedding (ECE) , which combines Genetic Programming with the boosting scheme to generate high-quality binary codes for large-scale data classification tasks. Besides, Self-taught Hashing (STH) (Zhang et al. 2010), in which a two-step scheme is effectively applied to learn hash functions, was also successfully utilized for visual retrieval. More hashing techniques can also be seen in , Cao et al. (2012), Song et al. (2014), Song et al. (2013), , , Lin et al. (2013).
Nevertheless, the above mentioned hashing methods have their limitations. Although the random projection based hashing methods, such as LSH, KLSH and SKLSH (Raginsky and Lazebnik 2009), can produce relatively effective codes, such simple linear hash functions cannot reflect the underlying relationship between the data points. Meanwhile, since the long codes are required for acceptable retrieval results via random projection based hashing, the storage space and the cost of computing the Hamming distance will be expensive. On the other hand, in terms of learning based hashing algorithms, most of them, e.g., Shakhnarovich (2005), Weiss et al. (2009), Liu et al. (2011, only focus on the relationship between data or sets rather than considering the combination of the intra-latent structure 1 of data and the inter-probability distribution between the high-dimensional Euclidean space and the low-dimensional Hamming space. To overcome these limitations above, in this paper, we propose a novel NMF-based approach called Latent Structure Preserving Hashing (LSPH) which can effectively preserve data probabilistic distribution and capture much of the locality structure from the high-dimensional data. In particular, the nonnegative matrix factorization can automatically learn the intra-latent information and the part-based representations of data, while the data probabilistic distribution preserving aims Fig. 1 The outline of the proposed method. Part-based latent information is learned from NMF with the regularization of data distribution. We propose two different versions of our algorithm, i.e., single layer LSPH and multi-layer LSPH. Specifically, ML-LSPH generates deep data representations which can theoretically lead to better performance for retrieval tasks with more complex data to maintain the similarity between high-dimensional data and low-dimensional codes. Moreover, incorporated with the representation of binary codes, the part-based latent information obtained by NMF based hashing, i.e., LSPH, could be regarded as independent latent attributes of samples. In other words, the binary codes determine whether the highdimensional data hits the corresponding latent attributes or not. Given an image, this kind of data-driven attributes can allow us to well describe an image and also may benefit zeroshot learning (Jayaraman and Grauman 2014;Lampert et al. 2014;Yu et al. 2013) for unseen image classification/retrieval in future work.
Specifically, because of the limitation of NMF, which cannot completely discover the latent structure of the original high-dimensional data, we provide a new objective function to preserve as much of the probabilistic distribution structure of the high-dimensional data as possible to the low-dimensional map. Meanwhile, we propose an optimization framework for the objective function and show the updating rules. Besides, to implement the optimization process, the training data are relaxed into a real-valued range. Then, we convert the real-valued representations into binary codes. Finally, we analyze the experimental results and compare them with several existing hashing algorithms. The outline of the proposed LSPH approach is depicted in Fig. 1.
LSPH is a linear hashing technique with a single-layer generative network and data distribution preserving constraints. Although it is an efficient binary coding method for large-scale retrieval tasks, such a single-layer generative network may lead to several limitations in the following cases as mentioned in [1]: (1) when it learns data which lie on or near a nonlinear manifold; (2) when it learns syntactic rela-tionships of given data; and (3) when it learns hierarchically generated data. The single-layer LSPH is apparently not fit for such cases. For instance, LSPH with a single-layer network can well tackle data with small intra-variations such as face images. However, for more complex data with extremely different viewpoints, additional degrees of freedom of the data will be required. In terms of large-scale image retrieval tasks, the sources of data can be very variant and even samples belonging to the same category can differ significantly. Naturally, the single-layer LSPH is not competent for similarity search on such heterogeneous databases. Therefore, in this paper, we also propose an extension of LSPH called multi-layer LSPH (ML-LSPH) with the multi-layer generative network [1] and distribution preserving constraints. ML-LSPH can deeply learn part-based latent information of data and preserve the joint probabilistic distribution for deep data representations. Applying the sigmoid function to each layer, ML-LSPH is a nonlinear architecture. Similar to recent deep neural networks (Hinton et al. 2006;Masci et al. 2011;Krizhevsky et al. 2012), ML-LSPH generates deep data representations which can theoretically lead to better performance than single layer LSPH for retrieval tasks with more complex data 2 in realistic scenarios. However, ML-LSPH is computationally more expensive during training and test phases compared to single layer LSPH. Thus, there exists a trade-off between ML-LSPH and LSPH in terms of performance and computational complexity, and the choice between these two versions depends on the requirement of the application. Besides, as ML-LSPH is a generalized framework of LSPH, it can easily shrink to LSPH if the number of the layers is set to 1. We evaluate our LSPH and ML-LSPH on three large-scale datasets: SIFT 1M, GIST 1M and TinyImage, and the results show that our methods significantly outperform the state-of-the-art hashing techniques. It is worthwhile to highlight several contributions of the proposed methods: -LSPH can learn compact binary codes uncovering the latent semantics and simultaneously preserving the joint probability distribution of data. -We utilize multivariable logistic regression to generate the hashing function and achieve the out-of-sample extension. -To tackle the data with more complex distribution, a multi-layer extension of LSPH (i.e., ML-LSPH) has been proposed for large-scale retrieval as well.
The rest of this paper is organized as follows. In Sect. 2, we give a brief review of NMF. The details of LSPH and ML-LSPH are described in Sects. 3 and 4, respectively. Section 5 reports the experimental results. Finally, we conclude this paper and discuss the future work in Sect. 6.

A Brief Review of NMF
In this section, we mainly review some related algorithms, focusing on Nonnegative Matrix Factorization (NMF) and its variants. NMF is proposed to learn the nonnegative parts of objects. Given a nonnegative data matrix X = [x 1 , · · · , x N ] ∈ R M×N ≥0 , each column of X is a sample data. NMF aims to find two nonnegative matrices U ∈ R M×D ≥0 and V ∈ R D×N ≥0 with full rank whose product can approximately represent the original matrix X, i.e., X ≈ U V . In practice, we always set D < min(M, N ). The target of NMF is to minimize the following objective function where · is the Frobenius norm. To optimize the above objective function, an iterative updating procedure was developed in Lee and Seung (1999) as follows: and normalization It has been proved that the above updating procedure can find the local minimum of L N M F . The matrix V obtained in NMF is always regarded as the low-dimensional representation while the matrix U denotes the basis matrix. Furthermore, there also exists some variants of NMF. Local NMF (LNMF) (Li et al. 2001) imposes a spatial localized constraint on the bases. In Hoyer (2004), sparse NMF was proposed and later, NMF constrained with neighborhood preserving regularization (NPNMF) (Gu and Zhou 2009) was developed. Besides, researchers also proposed graph regularized NMF (GNMF) (Cai et al. 2011), which effectively preserves the locality structure of data. Beyond these methods, Zhang et al. (2006) extends the original NMF with the kernel trick as kernelized NMF (KNMF), which could extract more useful features hidden in the original data through some kernel-induced nonlinear mappings. Additionally, a hashing method based on multiple kernels NMF was proposed in , where an alternate optimization scheme is applied to determine the combination of different kernels.
In this paper, we present a Latent Structure Preserving NMF framework for hashing (i.e., LSPH), which can effectively preserve the data intrinsic probability distribution and simultaneously reduce the redundancy of low-dimensional representations. Specifically, since the solution of standard NMF only focuses on optimizing matrix factorization to minimize Eq. (1), the obtained low-dimensional representation V lacks the data relationship information. In fact, most of previous NMF extensions are based on keeping the locality regularization to guarantee that, if the high-dimensional data points are close, the low-dimensional representations from NMF can still be close. However, this kind of regularization may lead to a low-quality factorization, since it ignores preserving the whole data distribution but only focuses on locality information. For a realistic scenario with noisy data, locality preserving regularization would even produce worse performance. Rather than locality-based graph regularization, we measure the joint probability of data by Kullback-Leibler divergence, which is defined over all of the potential neighbors and has been proved to effectively resist data noise (Maaten and Hinton 2008). This kind of measurement reveals the global structure such as the presence of clusters at several scales. To make LSPH more capable on data with more complex distributions, the multilayer LSPH (ML-LSPH) is also proposed, in which more discriminative, high-level representations can be learned from a multi-layer network with the distribution preserving regularization term. To the best of our knowledge, this is the first time that multi-layer NMF based hashing is successfully applied to feature embedding for large-scale similarity search. A preliminary version of our LSPH has been presented in Cai et al. (2015). In this paper, we include more details and experimental results and extend LSPH to ML-LSPH for more complex data in realistic retrieval applications.

Latent Structure Preserving Hashing
In this section, we mainly elaborate the proposed Latent Structure Preserving Hashing algorithm.

Preserving Data Structure with NMF
NMF is an unsupervised learning algorithm which can learn a parts-based representation. Theoretically, it is expected that the low-dimensional data V given by NMF can obtain locality structure from the high-dimensional data X . However, in real-world applications, NMF cannot discover the intrinsic geometrical and discriminating structure of the data space. Therefore, to preserve as much of the significant structure of the high-dimensional data as possible, we propose to minimize the Kullback-Leibler divergence (Xie et al. 2011) between the joint probability distribution in the highdimensional space and the joint probability distribution that is heavy-tailed in the low-dimensional space: In Eq. (4), P is the joint probability distribution in the high-dimensional space which can also be denoted as p i j . Q is the joint probability distribution in the low-dimensional space that can be represented as q i j . λ is the control of the smoothness of the new representation. The conditional probability p i j means the similarity between data points x i and x j , where x j is picked in proportion to their probability density under a Gaussian centered at x i . Since only significant points are needed to model pairwise similarities, we set p ii and q ii to zero. Meanwhile, it has the characteristics that p i j = p ji and q i j = q ji for ∀i, j. The pairwise similarities in the high-dimensional space p i j are defined as: where σ i is the variance of the Gaussian distribution which is centered on data point x i . Each data point x i makes a significant contribution to the cost function. In the low-dimensional map, using the probability distribution that is heavy tailed, the joint probabilities q i j can be defined as: This definition is an infinite mixture of Gaussians, which is much faster to evaluate the density of a point than the single Gaussian, since it does not have an exponential. This representation also makes the mapped points invariant to the changes in the scale for the embedded points that are far apart.
Thus, the cost function based on Kullback-Leibler divergence can effectively measure the significance of the data distribution . q i j models p i j is given by For simplicity, we define two auxiliary variables d i j and Z for making the derivation clearer as follows: Therefore, the gradient of function G with respect to v i can be given by Then ∂G ∂d i j can be calculated by Kullback-Leibler divergence in Eq. (7): Since is nonzero if and only if k = i and l = j, and k =l p kl = 1, the gradient function can be expressed as Eq. (11) can be substituted into Eq. (9). Therefore, the gradient of the Kullback-Leibler divergence between P and Q is Therefore, through combining the data structure preserving part in Eq. (4) and the NMF technique, we can obtain the following new objective function: where V ∈ {0, 1} D×N , X, U, V 0, U ∈ R M×D , X ∈ R M×N , and λ controls the smoothness of the new representation.
In most of the circumstances, the low-dimensional data only from NMF is not effective and meaningful for realistic applications. Thus, we introduce λK L(P Q) to preserve the structure of the original data which can obtain better results in information retrieval.

Relaxation and Optimization
Since the discreteness condition V ∈ {0, 1} D×N in Eq. (22) cannot be calculated directly in the optimization procedure, motivated by Weiss et al. (2009), we first relax the data V ∈ {0, 1} D×N to the range V ∈ R D×N for obtaining real-values. Then let the Lagrangian of our problem be: where matrices Φ and Ψ are two Lagrangian multiplier matrices. Since we have the gradient of C = λG: we make the gradients of L be zeros to minimize O f : In addition, we also have KKT conditions: Φ i j U i j = 0 and Ψ i j V i j = 0, ∀i, j. Then multiplying V i j and U i j in the corresponding positions on both sides of Eqs. (16) and (17) respectively, we obtain Note that Therefore, we have the following update rules for any i, j: All the elements in U and V can be guaranteed that they are nonnegative from the allocation. In Lee and Seung (2000), it has been proved that the objective function is monotonically non-increasing after each update of U or V . The proof of convergence about U and V is similar to the ones in Zheng et al. (2011), Cai et al. (2011. Once the above algorithm is converged, we can obtain the real-valued low-dimensional representation by a linear projection matrix. Since our algorithm is based on general NMF rather than Projective NMF (PNMF) (Yuan and Oja 2005;Guan et al. 2013), a direct projection does not exist for data embedding. Thus, in this paper, inspired by Cai et al. (2007), we consider using linear regression to compute our projection matrix. Particularly, we make the projection orthogonal by solving the Orthogonal Procrustes problem (Schönemann 1966) as follows: where P is the orthogonal projection. The optimal solution can be obtained by the following procedure: 1. use the singular value decomposition algorithm to decompose the matrix X T V = AΣ B T ; 2. calculate P = BΩ A T , where, Ω is a connection matrix as Ω = [I, 0] ∈ R D×M and 0 indicates all zeros matrix. Given data x ∈ R M×1 , its low-dimensional representation is v = Px. There are three advantages on using orthogonal projection according to Zhang et al. (2015): Firstly, the orthogonal projection can preserve the Euclidean distance between two points; Secondly, the orthogonal projection can distribute the variance more evenly across the dimensions; Thirdly, the orthogonal projection can learn maximally uncorrelated dimensions, which leads to more compact representations.

Hash Function Generation
The low-dimensional representations V ∈ R D×N and the bases U ∈ R M×D , where D M, can be obtained from Eq. (20) and Eq. (21), respectively. Then we need to convert the low-dimensional real-valued representations from V = [v 1 , · · · , v N ] into binary codes via thresholding: if the d-th element in v n is larger than a specified threshold, this real value will be represented as 1; otherwise it will be 0, where d = 1, · · · , D and n = 1, · · · , N .
In addition, a well-designed semantic hashing should also be entropy maximizing to ensure its efficiency (Baluja and Covell 2008). Meanwhile, from the information theory, through having a uniform probability distribution, the source alphabet can reach a maximal entropy. Specifically, if the entropy of codes over the corpus is small, the documents will be mapped to a small number of codes (hash bins). In this paper, the threshold of the elements in v n can be set to the median value of v n , which can satisfy entropy maximization. Therefore, half of the bit-strings will be 1 and the other half will be 0. In this way, the real-value code can be calculated into a binary code (Yu et al. 2014).
However, from the above procedure, we can only obtain the binary codes of the data in the training set. Therefore, given a new sample, it cannot directly find a hash function. To solve such an "out-of-sample" problem, in our approach, we are inspired by the "self-taught" binary coding scheme (Zhang et al. 2010) to use the logistic regression (Hosmer and Lemeshow 2004) which can be treated as a type of probabilistic statistical classification model to compute the hash code for unseen test data. Specifically, we learn a square projection matrix via logistic regression, which can be regarded as a rotation of V . This kind of transformation can make the codes more balanced (Gong et al. 2013;) and lead to better performance compared with directly binarizing V with the median value calculated from training data. To make it more convincing, we also show the performance difference in the later section. Before obtaining the logistic regression cost function, we define that the binary code is represented asV = [v 1 , · · · ,v N ], wherê v n ∈ {0, 1} D and n = 1, · · · , N . Therefore, the training set can be considered as The vector-valued regression function which is based on the corresponding regression matrix Θ ∈ R D×D can be represented as Therefore, with the maximum log-likelihood criterion for the Bernoulli-distributed data, our cost function for the corresponding regression matrix can be defined as: where log(·) is the element-wise logarithm function and 1 is an D × 1 all ones matrix. We use δ Θ 2 as the regularization term in logistic regression to avoid overfitting.
To find the matrix Θ which aims to minimize J(Θ), we use gradient descent and repeatedly update each parameter using a learning rate α. The updating equation is shown as follows: The updating equation stops when the norm of difference between Θ (t+1) and Θ (t) , i.e., ||Θ (t+1) − Θ (t) || 2 , is smaller than a small value. Then we can obtain the regression matrix Θ. For a new coming test data X new ∈ R M×1 , then its lowdimensional representation is V new = P X new . Note that each entry of h Θ is a sigmoid function, the hash codes for a new coming sample X new ∈ R M×1 can be represented as: where · means the nearest integer function for each entry of h Θ . Specifically, since the output of logistic regression i.e., h Θ (P X new ), indicates the probability of "1" for each entry, · is equivalent to binarizing each bit by probability 0.5. Thus, if the probability of a bit from h Θ (P X new ) is larger than 0.5, it will be represented as 1, otherwise 0. For example, through Eq. (26), vector h Θ (P X new ) = [0.17, 0.37, 0.42, 0.79, 0.03, 0.92, · · · ] can be expressed as [0, 0, 0, 1, 0, 1, · · · ]. Up to now, we can obtain the Latent Structure Preserving Hashing codes for both training and test data. The procedure of LSPH is summarized in Algorithm 1.

Input:
The training matrix X ∈ R M×N ; the objective dimension (code length) D of hash codes; the learning rate α for logistic regression; the regularization parameters {δ, λ}.

Output:
The basis matrix U , the orthogonal projection P and the regression matrix Θ. 1: Initialize U and V with uniformly distributed random values between 0 and 1. 2: repeat 3: Compute the low-dimensional representation matrix V and the basis matrix U via Eqs. (20) and (21), respectively; 4: until convergence 5: SVD decompose the matrix X T V to obtain AΣ B T and calculate P = BΩ A T ; 6: Obtain the regression matrix Θ through Eq. (25) and the final LSPH encoding for each sample is defined in Eq. (26). LSPH (ML-LSPH). ML-LSPH aims to generate more informative high-level representations compared with single-layer LSPH for data with complex distributions. Once the representation by ML-LSPH is computed, the final hashing functions are also obtained through multivariable logistic regression, same as LSPH mentioned above. Given data matrix X ∈ R M×N , inspired by recent deep learning algorithms and multi-layer NMF [1], Trigeorgis et al. (2014), we can extract latent data attributes by incorporating our LSPH algorithm with a multi-layer structure as illustrated in Fig. 2. Similar to the learning of the above representation matrix V , a matrix sequence V 1 , · · · , V n can be obtained by solving the following optimization problems: . . .
In this way, V i , i = 1, · · · , n, are the hidden factors of each layer. By introducing the nonlinear function g(·) into the network, these hidden factors are generated by the following rules: Then our task is to minimize the following objective function: Let us denote G i = K L(P||Q (i) ) and is the Hadamard product (element-wise product). By taking the derivatives of F with respect to U i and V i , we have: where matrices N i , R i ∈ R D i−1 ×N are calculated by the following rules: for i = 1, · · · , n − 1, with the initialization: Besides, the derivatives of G with respect to U i and V i are calculated by Eqs. (27) and (28). With the derivation in Sect. 3, we have the derivative where v j k is the k-th column of V j , k = 1, · · · , N and j = 1, · · · , n. To ensure that every element in U i and V i is nonnegative, we use the following symbols to split the above derivatives as: where Then we can define two matrix sequences S l and M l as follows: where In this way, the derivatives of G j with respect to U i and V i , i.e., Eqs.
(27) and (28), will be: Substitute the above equations into Eqs. (33) and (34), we obtain: Finally, similar to the procedure in Sect. 3, the update rules for multi-layer LSPH (ML-LSPH) are shown in Eqs. (29) and (30), where 0 < γ < 1 is the learning rate and i = 1, · · · , n. The convergence property of the above iteration is similar to the one in [1]. Besides, for better understanding our ML-LSPH, we aim to unify the LSPH and ML-LSPH under a same framework. Thus, in our implementation, the function g(·) applied on U 1 and V 1 is regarded as identity function f (x) = x. The function g(·) for U i and V i , i = 2, · · · , n is played by nonlinear sigmoid function f (x) = 1 1+e x . In this way, when the number of layers n = 1, the ML-LSPH will shrink to the ordinary singlelayer LSPH. It is noteworthy that we could theoretically formulate our ML-LSPH to an arbitrary number of layers according the above algorithms. However, for realistic applications with complex data distributions, the number of layers is always less than 3, since when the number of layers increases, the accumulative reconstruction error will cause the non-convergence of the proposed model (Trigeorgis et al. 2014).
For the hash code generating phase, it is similar to LSPH. In particular, acquired the low-dimensional representation V n in the n-th layer, we first solve the Orthogonal Procrustes problem min P T P=I P X − V n to achieve the orthogonal projection P. The optimal solution can be obtained by the following procedure: 1. use the singular value decomposition algorithm to decompose the matrix X T V n = AΣ B T ; 2. calculate P = BΩ A T , where, Ω is a connection matrix as Ω = [I, 0] ∈ R D×M and 0 indicates all zeros matrix. For a new coming test data x new ∈ R M×1 , the low-dimensional representation in the n-th layer is v new n = Px new and the binary codes are calculated byv new where Θ is obtained by the similar multi-variable regression scheme. The procedure of ML-LSPH is summarized in Algorithm 2.

Input:
The training matrix X ∈ R M×N ; the dimensions for n layers D 1 , · · · , D n ; the learning rate γ for the multi-layer NMF structure; the learning rate α for logistic regression; the regularization parameters {δ, λ}.

Output:
The n-th layer representation V n , the orthogonal projection P and the regression matrix Θ. 1: Initialize U i and V i with uniformly distributed random values between 0 and 1, i = 1, · · · , n. 2: repeat 3: Compute the basis matrix U i and the low-dimensional representation matrix V i via Eqs. (29) and (30) respectively for each layer; 4: until convergence 5: SVD decomposes the matrix X T V n to obtain AΣ B T and calculate P = BΩ A T ; 6: Obtain the regression matrix Θ through Eq. (25) and the final ML-LSPH encoding for each sample is defined by h Θ (P X ) .

Batch-Based Learning Scheme
With the number of the layers increasing, the computational costs will inevitably increase as well in the current multi-layer network architecture of ML-LSPH. In order to effectively reduce the computational complexity on large-scale data, we adopt a random batch-based learning strategy (RBLS) in the iteration opti- as mentioned above, which is still regarded as linear complexity in terms of N and not very demanding for large-scale data processing. However, the real bottleneck of the optimization procedure is the calculation of the KL divergence for each layer, specifically, the similarity matrices P and Q (m) , due to the complexity of O (N 2 D). Therefore, in our implementation, we adopt RBLS to effectively reduce the complexity for computing P and Q (m) in ML-LSPH. In detail, for each step of updating P and Q (m) , we randomly select a small subset of the whole training set. Then we only need to compute the pairwise similarity of this subset and the rest of the elements of P and Q (m) are replaced by zeros: where m indicates the layer index. The illustration of proposed RBLS are also shown in Fig. 3. If we assume the size of the small subset is , the complexity of our RBLS will be reduced from O(N 2 D) to O( 2 D). Usually, the can be set as = 1/100N . It is noteworthy that only the computation of P and Q (m) are replaced by the above RBLS trick and other parts of the algorithm in ML-LSPH are still the same as mentioned before. In this way, our multi-layer LSPH becomes scalable for large-scale data.

Computational Complexity Analysis
In this section, we will discuss the computational complexity of LSPH and ML-LSPH. The computational complexity of LSPH consists of three parts. The first part is for computing NMF, the complexity of which is O(N M D) , where N is the size of the dataset, M and D are the dimensionalities of the high-dimensional data and the low-dimensional data respectively. The second part is to compute the cost function Eq. (7) in the objective function which has the complexity O (N 2 D). The last part is the logistic regression procedure whose complexity is O (N D 2 ). Therefore, the total computational complexity of LSPH is: where t is the number of iterations. It is obvious that single-layer ML-LSPH is actually LSPH. The computational complexity of ML-LSPH with multiple layers consists of several NMF optimizations, the computation of matrices P and Q ( j) , and the logistic regression procedure. With the above discussion, the com- where is the batch size.

Experiments and Results
In this section, we systematically evaluate the proposed LSPH and multi-layer LSPH (ML-LSPH) on three largescale datasets. The relevant experimental results and data visualization will be included in the rest of this section. All the experiments are completed using MatLab2014a on a work-    station with a 12-core 3.2GHz CPU and 120GB of RAM running the Linux OS.

Evaluation on LSPH
In this subsection, we first apply the proposed single-layer LSPH algorithm method for large-scale similarity search tasks. Two realistic datasets are used to evaluate all methods: SIFT 1M: it contains one million 128-dim local SIFT feature vectors (Lowe 2004). GIST 1M: it contains one million 960-dim global GIST feature vectors (Oliva and Torralba 2001). The two datasets are publicly available. 3 3 http://corpus-texmex.irisa.fr.
For both SIFT 1M and GIST 1M, we respectively take 10K images as the queries by random selection and use the remaining to form the gallery database. To construct the training set, 200, 000 samples from the gallery database are randomly selected for all of the compared methods. Additionally, we also randomly choose another 50, 000 data samples as a cross-validation set for parameter tuning. In the querying phase, the returned points are regarded as true neighbors if they lie in the top 2 percentile points closest to a query for both two datasets. Since the Hamming lookup table is fast with hash codes, we will use the Hamming lookup table to measure our retrieval tasks. We evaluate the retrieval results in terms of the Mean Average Precision 4 (MAP) and the precisionrecall curves. Additionally, we also report the training time and the test time (the average searching time used for each query) for all methods.

The Selected Methods and Settings
In this experiment, we compare LSPH with the 10 selected popular hashing methods including LSH (Gionis et al. 1999), BSSC (Shakhnarovich 2005), RBM (Salakhutdinov and Hinton (2007)), SpH (Weiss et al. 2009), STH (Zhang et al. 2010), AGH (Liu et al. 2011), ITQ (Gong et al. 2013), KLSH (Kulis and Grauman 2009), PCAH ) and CH (Lin et al. 2013). In these methods, for BSSC, through the labeled pairs scheme in the boosting framework, it can obtain weights and thresholds for every hash function. RBM will be trained with several 100-100 hidden layers without fine-tuning. According to KLSH, 500 training samples and the RBF-kernel are used to output the empirical kernel map, in which we always set the scalar parameter σ to an appropriate value on each dataset. For ITQ, the number of the iterations is set as 50. In AGH with two-layer, we consider the number of the anchor points k as 200 and the number of the nearest anchors s in sparse coding as 50. CH has the same anchor-based sparse coding setting with AGH. All of the 10 methods are evaluated on different lengths of the codes, e.g., 16, 32, 48, 64, 80 and 96. Under the same experimental setting, all the parameters used in the compared methods have been strictly chosen according to their original papers. 4 The ground-truth is defined by top 2 percentile nearest neighbors of a query via linear scan based on the Euclidean distance.
In the experiments of our LSPH method, all the data are first normalized into the range of [0, 1], since the nonnegative constraint is required in our framework. We also use the validation set to tune our hyper-parameters. Particularly, for each dataset, we select one value from {0.01, 0.02, · · · , 0.10} as the optimal learning rate α of multi-variable logistic regression through 10-fold cross-validation on the validation set. The regularization parameter λ is determined from one of {10 −2 , 10 −1 , 1, 10 1 , 10 2 , 10 3 }, which yields the best performance via 10-fold cross-validation on the validation set. The regularization parameter δ in the hash function generation is fixed as δ = 0.35. We limit the maximum number of iterations with 1000 in LSPH learning phase, as well.

Results Comparison
We demonstrate MAP curves on the SIFT 1M and GIST 1M datasets compared with different methods in Fig. 4. From the general tendency, the accuracies on the GIST 1M dataset are lower than those on the SIFT 1M dataset, since features in the GIST 1M has higher dimensions with larger variations than those on SIFT 1M. In detail of Fig. 4, ITQ always achieves higher Mean Average Precision (MAP) and gets a consistent increasing condition with the change of the code length on both datasets. Furthermore, MAP of CH also proves to be competitive but a little lower than ITQ. Interestingly, on both the SIFT 1M and GIST 1M datasets, the MAP of SpH and STH are always "rise-then-fall" when the length of code increases. Due to the use of random projection, LSH and KLSH have a low MAP when the code length is short. Moreover, PCAH always gets decreasing accuracies when the code length increases. For our method LSPH, it achieves the highest performance among all the compared methods on both datasets. The proposed LSPH algorithm can automatically exploit the latent structure of the original data and simultaneously preserve the consistency of distribution between the original data and the reduced representations. The above properties of LSPH allow it to achieve better performance in large-scale visual retrieval tasks. Fig. 5 also shows the precision-recall curves with the code length of 48 bits on both SIFT 1M and GIST 1M datasets with the top 2 percentile nearest neighbors as the ground truth. From all these figures, we can further discover that, for both two datasets, LSPH achieves apparently better performance than other hashing methods by comparing the Mean Average Precision (MAP) and Area Under the Curve (AUC). Additionally, to further illustrate the necessity of using logistic regression binarization rather than direct median value binarization as mentioned in Sect. 3.3, we carry out comparison experiments on both SIFT 1M and GIST 1M datasets. In Table 1, it is easy to observe that logistic regression binarization can achieve consistently higher MAP across all code lengthes than the direct median value binarization scheme.   For ML-LSPH: d loss1 = X − g(U 1 g(U 2 · · · g(U n V n ))) 2 versus number of iterations. e loss2 = n i=1 K L(P||Q (i) ) versus number of iterations. f T otalloss = X −g(U 1 g(U 2 · · · g(U n V n ))) 2 +λ n i=1 K L(P||Q (i) ) versus number of iterations. Zoom in for better viewing Meanwhile, the training and test time for all the methods are listed in Tables 2 and 3, as well. Considering the training time, the random projection based algorithms are relatively more efficient, especially the LSH. While, RBM takes the most time for training, since it is based on a time-consuming deep learning method. Our method LSPH is significantly more efficient than STH, BSSC and RBM, but slightly slower than ITQ, AGH and SpH. It is noteworthy that once the optimal hashing functions of our method are obtained from the training phase, the optimized hashing functions will be fixed and directly used for new data. In addition, with the rapid development of silicon technologies, future computers will be much faster and even the training will become less a problem. In terms of the test phase, LSH is the most efficient methods as well, since only a simple matrix multiplication and a thresholding are needed to obtain the binary codes. AGH and SpH always take more time for the test phase. Our LSPH has the competitive efficiency with STH. Therefore, in general, it can be concluded that LSPH is an effective and relatively efficient method for the large-scale retrieval tasks.

Evaluation on ML-LSPH
In this subsection, the multi-layer LSPH is evaluated on the TinyImage dataset, which is a subset containing 500,000 images 5 from ***80 Million Tiny Images (Torralba et al. 2008a, b). Some example images from the TinyImage dataset are illustrated in Fig. 6. We further take 1K images as the queries by random selection and use the remaining to form the gallery database. Considering the cost of computation in multi-layer networks, in this experiment, only 100, 000 randomly selected samples from the gallery database form the training set. Similar to the experiments of LSPH, another 50, 000 data samples are also randomly chosen as a cross-validation set for parameter tuning. For image searching tasks, given an image, we describe it with 512dimensional GIST descriptors (Oliva and Torralba 2001) in this experiment and then learn to hash these descriptors with all compared methods. In the querying phase, a returned point is regarded as a neighbor if it lies in the top ranked 500 points for the TinyImage dataset. We evaluate the retrieval results through Hamming distance ranking and report the Mean Average Precision (MAP) and the precision-recall curves by changing the number of top ranked points. Additionally, we also report the parameter sensitive analysis and visualize some retrieved images of compared methods on this dataset.
To avoid confusion, in this experiment, we only compare with LSH, PCAH, ITQ, AGH, RBM and SpH, which have shown distinctive performance according to the results in the previous comparison. Besides, we also add a new hashing technique SKLSH in this experiment. Note that, RBM here has 100-100-100 three hidden layers without fine-tuning. All of the above methods are then evaluated on six different lengths of codes (32,48,64,80,96,128). Under the same experimental setting, all the parameters used in the compared methods have been strictly chosen according to their original papers.
For our ML-LSPH, the data zero-one normalization and hyper-parameter selection follow the same scheme as those in LSPH. Besides, to better reach the local minimum loss in ML-LSPH, the learning rate γ for iterative optimization is considered. In this experiment, we fix γ = 0.5 for ML-LSPH. More importantly, for the reason mentioned above that when the number of the layers increases, the accumulative reconstruction error will cause the non-convergence problem of the proposed model (Trigeorgis et al. 2014), we evaluate ML-LSPH with n = 2 layers, which is the same setting as in Trigeorgis et al. (2014). We further set the dimension of the middle layer 6 (i.e., V 1 ) to 256 on this dataset. 6 After attempting various network architectures with different dimensions of middle layers, the current ML-LSPH network is the best option for our task. Similar to the deep learning model (Krizhevsky et al. 2012;Szegedy et al. 2014). It is noteworthy that the dimensions of middle layers are quite sensitive to the data distribution (i.e., dataset) and there is no particular proof to explain what length of dimension for middle layers is better. Thus, we recommend to try different structure settings in order to determine what kind of structure can achieve the best performance.   7 illustrates the MAP curves of all compared algorithms on the TinyImage dataset. Our ML-LSPH algorithm achieves slightly lower MAP than ITQ when the code length is less than 48 bits but consistently outperforms all other compared methods in every length of code. RBM with 3 layers can also produce competitive search accuracies on this dataset. Different to other hashing techniques, the performance of PCAH decreases with the increase of the code length. The similar phenomenon has appeared in the previous evaluation on SIFT 1M and GIST 1M datasets. The performance of AGH, SpH and LSH is consistent with that in the previous experiments.
Besides, Fig. 8 shows a series of precision-recall curves with different code lengths on the TinyImage dataset with the 500 nearest neighbors as the ground truth. By comparing the Area Under the Curve (AUC), our ML-LSPH achieves apparently better performance than other methods on relatively long bits (code length ≥ 48 bits) and the performance slightly goes down when short hash codes are adopted. Moreover, in Table 4 and Fig. 9, we also compare the performance between ML-LSPH and LSPH in terms of the MAP and AUC on all three datasets. The ML-LSPH can achieve consistently better results than LSPH, since large intra-class variations in Tiny-Image cause complex and noisy data distributions, which are more difficult to handle by LSPH. Besides, Fig. 10 illustrates the convergence of the proposed LSPH and ML-LSPH on TinyImage with the code length of 64. We can clearly observe that, for LSPH, the loss of X −U 1 V 1 2 dramatically drops down when the number of iterations increases. While, for ML-LSPH, the loss of X − g(U 1 g(U 2 · · · g(U n V n ))) 2 first climbs up when the number of iterations is less than 50 and then goes down. With the batch-based learning scheme, the total losses of both LSPH and ML-LSPH can steadily decrease with little fluctuation.
Additionally, in Fig. 11, we also compare the retrieval performance of ML-LSPH with respect to the regularization parameter λ along different code lengths via cross-validation. When tuning the parameter λ from 0.01 to 1000 with a scale factor of 10, the MAP curves of ML-LSPH appear to be relatively stable and insensitive to λ. For code lengths equal to 64 bits and 80 bits, the best performance occurs when λ = 10. However, for the rest of the code lengths, ML-LSPH can achieve the highest retrieval accuracies with λ = 100. Specifically, to illustrate the effectiveness of the data distribution preserving (KL divergence) regularization term in ML-LSPH, we also compare the algorithm without using the KL divergence term in Table 5. The results indicate that combining the multi-layer NMF network with the data distribution preserving term could always gain better performance. Finally, the top-ranked retrieval results using compared methods on some typical queries are illustrated in Fig. 12. It can be observed that the retrieved images via ML-LSPH have more semantic consistency with the query images.

Conclusion and Future Work
In this paper, we have proposed the Latent Structure Preserving Hashing (LSPH) algorithm, which can find a wellstructured low-dimensional data representation through the Nonnegative Matrix Factorization (NMF) with the probabilistic structure preserving regularization part, and then the multi-variable logistic regression is effectively applied to generate the final hash codes. To better tackle the data with complex and noisy data distributions, a multi-layer LSPH (ML-LSPH) extension has also been developed in this paper. Extensive experiments on three large-scale datasets have demonstrated that our algorithms can lead to competitive hashing performance for large-scale retrieval tasks.
As we mentioned in the introduction of the paper, the hash codes learned from both LSPH and ML-LSPH can be regarded as independent latent attributes due to the property of NMF. In future work, this kind of learned data-driven attributes will be further explored with zero-shot learning for unseen image classification, retrieval and annotation tasks.