Streaming Binary Sketching based on Subspace Tracking and Diagonal Uniformization

TL;DR

This paper introduces an online method for generating binary sketches from high-dimensional streaming data, enabling efficient similarity search without storing the entire dataset.

Contribution

The proposed algorithm provides a fully online, memory-efficient way to produce binary embeddings with convergence guarantees for streaming high-dimensional data.

Findings

Effective binary sketches for high-dimensional streams.

Low time complexity and no need for dataset storage.

Successful experiments on real data for nearest neighbor search.

Abstract

In this paper, we address the problem of learning compact similarity-preserving embeddings for massive high-dimensional streams of data in order to perform efficient similarity search. We present a new online method for computing binary compressed representations -sketches- of high-dimensional real feature vectors. Given an expected code length $c$ and high-dimensional input data points, our algorithm provides a $c$-bits binary code for preserving the distance between the points from the original high-dimensional space. Our algorithm does not require neither the storage of the whole dataset nor a chunk, thus it is fully adaptable to the streaming setting. It also provides low time complexity and convergence guarantees. We demonstrate the quality of our binary sketches through experiments on real data for the nearest neighbors search task in the online setting.

Equations

$${\Sigma }_V$$

$$R$$

$$\left(\begin{array}{cc}{a}^{\prime }& b^{\prime }\\ b^{\prime }& d^{\prime }\end{array}\right)\stackrel{def}{=}\left(\begin{array}{cc}c& -s\\ s& c\end{array}\right)\left(\begin{array}{cc}a& b\\ b& d\end{array}\right)\left(\begin{array}{cc}c& s\\ -s& c\end{array}\right).$$

Full text

**Index Terms— ** Sketching, streaming, subspace tracking, Givens rotations ††A. Morvan is partly supported by the DGA (French Ministry of Defense). ††Copyright 2018 IEEE. Published in the IEEE 2018 International Conference on Acoustics, Speech, and Signal Processing (ICASSP 2018), scheduled for 15-20 April 2018 in Calgary, Alberta, Canada. Personal use of this material is permitted. However, permission to reprint/republish this material for advertising or promotional purposes or for creating new collective works for resale or redistribution to servers or lists, or to reuse any copyrighted component of this work in other works, must be obtained from the IEEE. Contact: Manager, Copyrights and Permissions / IEEE Service Center / 445 Hoes Lane / P.O. Box 1331 / Piscataway, NJ 08855-1331, USA. Telephone: + Intl. 908-562-3966.

1 Introduction

In large-scale machine learning applications such as computer vision or metagenomics, learning similarity-preserving binary codes is critical to perform efficient indexing of large-scale high-dimensional data. Storage requirements can be reduced and similarity search sped up by embedding high-dimensional data into a compact binary code. A classical method is Locality-Sensitive Hashing (LSH) [1] for nearest neighbors search: 1) input data with high dimension $d$ are projected onto a lower $c$-dimensional space through a $c\times d$ random projection with i.i.d. Gaussian entries, 2) then a hashing function is applied to the resulting projected vectors to obtain the final binary codes. Two examples for the hashing function are cross-polytope LSH [2] which returns the closest vector from the set $\{\pm 1e_{i}\}_{1\leq i\leq c}$ where $\{e_{i}\}_{1\leq i\leq c}$ stands for the canonical basis, and hyperplane LSH [3], the $\operatorname{sign}$ function applied pointwise. To reduce the storage cost of the projection matrix and the matrix-vector products computation times ($O(c\times d$)), a structured pseudo-random matrix can be used instead [3, 4] with a reduced time complexity of $O(d\log c)$ thanks to fast Hadamard and Fourier transforms.

In the context of nearest neighbors search or classification, the accuracy of the data sketches can be improved by learning this projection from data [5, 6, 7, 8, 9]. As Principal Component Analysis (PCA) is a common tool for reducing data dimensionality, data are often projected onto the first $c$ principal components. But PCA alone is not sufficient. Indeed, the $c$ first principal components are chosen with a decreasing order of explained variance: principal directions with higher variance carry more information. Thus, associating each of the $c$ directions to one of the $c$ bits is equivalent to giving more weights to less informative directions and will lead to poor performance of the obtained sketches. To remedy this problem, after data have been projected on the first principal components of the covariance matrix, a solution consists in applying a suitable rotation on the projected data before performing the hashing function: it balances variance over the principal components. In work from [10], a random rotation is used giving quite good results. In ITerative Quantization (ITQ) [7] or in Isotropic Hashing (IsoHash) [8], the rotation is rather learned. In ITQ, the rotation is iteratively computed by solving an orthogonal Procustes problem: the alternating minimization of the quantization error of mapping data to the vertices of the $2^{c}$ hypercube. This technique is currently the state-of-the-art for computing similarity-preserving binary codes based on PCA. ITQ and IsoHash come with a major drawback, though. They are completely offline since the whole dataset needs to be stored for computing the $c$ principal components. This can be prohibitive when dealing with lots of high-dimensional data.

Contributions: We introduce a streaming algorithm with convergence guarantees where data is seen only once and principal subspace plus the balancing rotation are updated as new data is seen. To obtain the principal subspace and the rotation, this requires additionally only the storage of two $c\times c$ matrices, instead of the whole initial and projected datasets as for ITQ or IsoHash. Our algorithm outperforms the known state-of-the-art online unsupervised method Online Sketching Hashing (OSH) [11] while being far less computationally demanding.

2 Related work

Two paradigms exist to build hash functions [12]: data-independent [1, 13, 14] and data-dependent methods. The latter ones learn the hash codes from a training set and perform better. The learning can be unsupervised [5, 15, 7, 8, 16, 17, 9, 18] aiming at preserving distances in the original space or (semi-)supervised which also tries to preserve label similarity [6, 19]. Some recent hashing functions involve deep learning [20, 21]. When the dataset is too large to be loaded into memory, distributed [22] and online hashing techniques [23, 11, 24] have been developed. Online Hashing (OKH) [23] learns the hash functions from a stream of similarity-labeled pair of data with a ”Passive-Aggressive” method. Recent approach from [24] relies on Mutual Information. Although claimed to be unsupervised, it is supervised as similarity labels are also needed between pairs of data. In Online Sketching Hashing (OSH) [11], the binary embeddings are learned from a maintained sketch of the dataset with a smaller size but preserving the property of interest. The proposed algorithm belongs to this latter category of online unsupervised hyperplanes-based hashing methods.

3 The proposed online unsupervised model for binary quantization

3.1 Notations and problem statement

We have a stream of $n$ data points $\{x_{t}\in\mathbb{R}^{d}\}_{1\leq t\leq n}$ supposed to be zero-centered. The goal is to have $b_{t}=\operatorname{sign}(\tilde{W}_{t}x_{t})\in\{-1,1\}^{c}$ for $t=1\ldots n$ where $c$ denotes the code length, $c\ll d$, s.t. for each bit $k=1\ldots c$, the binary encoding function is defined by $h_{k}(x_{t})=\operatorname{sign}(\tilde{w}^{T}_{k,t}x_{t})$ where $\tilde{w}_{k,t}$ are column vectors of hyperplane coefficients and $\operatorname{sign}(x)=1$ if $x\geq 0$ and $-1$ otherwise which is applied component-wise on coefficients of vectors. So $\tilde{w}_{k,t}^{T}$ is a row of $\tilde{W}_{t}$ for each $k$, with $\tilde{W}_{t}\in\mathbb{R}^{c\times d}$. $b_{t}$ is computed and returned before $x_{t+1}$ is seen by using $x_{t}$ and $\tilde{W}_{t-1}$ solely. We consider in this paper the family of hash functions $\tilde{W}_{t}$ s.t. $\tilde{W}_{t}=R_{t}W_{t}$ where $W_{t}$ is the linear dimension reduction embedding applied to data and $R_{t}$ is a suitable $c\times c$ orthogonal matrix. Here, we take $W_{t}$ as the matrix whose row vectors $w_{k,t}^{T}$ are the $c$ first principal components of the covariance matrix $\Sigma_{X,t}=X_{t}X_{t}^{T}$ where $X_{t}$ denotes the dataset seen until data $t$. So the challenge is in tracking the principal subspace in an online fashion as new data is seen and in defining an appropriate orthogonal matrix111In the sequel, we use equally the term orthogonal matrix or rotation. to rotate the projected data onto this principal subspace.

Works from [7, 10] argued that defining an orthogonal transformation applied to the PCA-projected data which tends to balance the variances over the PCA-directions improves the quality of the hashing but $R$ was first learned to uniformize the variances for different directions only in work from [8] with Isotropic Hashing (IsoHash). In the latter, the problem is explicitly described as learning a rotation which produces isotropic variances of the PCA-projected dimensions. We propose here a new simpler online method UnifDiag for learning a rotation to uniformize the diagonal of the covariance of projected data after tracking the principal subspace $W_{t}$ from the data stream with Fast Orthonormal PAST (Projection Approximation and Subspace Tracking) [25], also named OPAST. At each iteration $t$, OPAST guarantees the orthonormality of $W_{t}$ rows and costs only $4dc+O(c^{2})$ flops while storing only $W_{t}$ and a $c\times c$ matrix. The rotation learning is completely independent from the principal subspace tracking. Indeed, any other method for online PCA [26, 27] or subspace tracking [25] can be plugged before UnifDiag.

3.2 UnifDiag: the proposed diagonal uniformization-based method for learning a suitable rotation

Let the $c\times c$ symmetric matrix $\Sigma_{V,t}=V_{t}V_{t}^{T}$ be the covariance matrix of projected data seen until data $t$ $V_{t}=W_{t}X_{t}$. $R_{t}$ is learned for each $t$ to balance the variance over the $c$ directions given by the $c$ principal components of $\Sigma_{V,t}$. $\Sigma_{V,t}$ is easy to update dynamically and we adapt OPAST algorithm to perform this while computing $W_{t}$. In the sequel, for clarity we drop the subscript $t$. Let us consider the $c$ diagonal coefficients of $\Sigma_{V}$: $\sigma^{2}_{1},...,\sigma^{2}_{c}$ s.t. $\sigma^{2}_{1}\geq...\geq\sigma^{2}_{c}$. As $\Sigma_{V}$ is symmetric, $\operatorname{Tr}(\Sigma_{V})=\sum_{i=1}^{c}\sigma^{2}_{i}=\sum_{i=1}^{c}\lambda_{i}$ where $\operatorname{Tr}$ stands for the Trace application and $\lambda_{1},...,\lambda_{c}$ are the $c$ first eigenvalues of $\Sigma_{X}$ s.t. $\lambda_{1}\geq...\geq\lambda_{c}$222If $W$ is exactly the $c$ first eigenvectors of $\Sigma_{X}$ – for instance, if $W$ is obtained through PCA –, then $\forall i\in\{1,...,c\}$, $\sigma_{i}^{2}=\lambda_{i}$.. Balancing variance over the $c$ directions can be seen as equalizing the diagonal coefficients of $\Sigma_{V}$ s.t. $\sigma^{2}_{1}=...=\sigma^{2}_{c}\overset{def}{=}\tau$. Since $\Sigma_{V}$ is symmetric, $\operatorname{Tr}(R\ \Sigma_{V}R^{T})=\operatorname{Tr}(\Sigma_{V})$. So in order to have $R\ \Sigma_{V}R^{T}$ with equal diagonal coefficients, we should set $\tau=\operatorname{Tr}(\Sigma_{V})/c$. So, similarly to IsoHash, we formulate the problem of finding $R$ as the problem of equalizing the diagonal coefficients of $\Sigma_{V}$ to the value $\tau=\operatorname{Tr}(\Sigma_{V})/c$. Our proposed optimal orthogonal matrix $R$ is built as a product of $c-1$ Givens rotations $G(i,j,\theta)$ described by Def. 3.1.

Definition 3.1

A Givens rotation $G(i,j,\theta)$ is a matrix of the form:

$G(i,j,\theta)=\begin{bmatrix}1&\cdots&0&\cdots&0&\cdots&0\\ \vdots&\ddots&\vdots&&\vdots&&\vdots\\ 0&\cdots&c&\cdots&-s&\cdots&0\\ \vdots&&\vdots&\ddots&\vdots&&\vdots\\ 0&\cdots&s&\cdots&c&\cdots&0\\ \vdots&&\vdots&&\vdots&\ddots&\vdots\\ 0&\cdots&0&\cdots&0&\cdots&1\end{bmatrix}$**

where $c=\cos(\theta)$ and $s=\sin(\theta)$ are at the intersections of the $i$-th and $j$-th rows and columns. The nonzero elements are consequently: $\forall k\neq i,\,j,\ g_{k,k}=1$, $g_{i,i}=g_{j,j}=c$, $g_{j,i}=-s$ and $g_{i,j}=s$ for $i>j$.

The computation of $R$ follows the iterative Jacobi eigenvalue algorithm known as diagonalization process [28]:

[TABLE]

Note that left (resp. right) multiplication by $G(i,j,\theta)$ only mixes $i$-th and $j$-th rows (resp. columns). The update from Eq. 1 only modifies $i$-th and $j$-th rows and columns of $\Sigma_{V}$. The two updated diagonal coefficients $(i,i)$ and $(j,j)$ only depend on $\Sigma_{V,i,i}$, $\Sigma_{V,j,j}$, $\Sigma_{V,j,i}$ and $\theta$ which reduces the optimization of $\theta$ to a $2$-dimensional problem, a classical trick when using Givens rotations [28]. Then we have Th. 3.1 by defining: $a\overset{def}{=}\Sigma_{V,j,j}$, $d\overset{def}{=}\Sigma_{V,i,i}$, $b\overset{def}{=}\Sigma_{V,j,i}=\Sigma_{V,i,j}$,

[TABLE]

Theorem 3.1

If $\operatorname{min}(a,d)\leq\tau\leq\operatorname{max}(a,d)$ (sufficient condition)333Th. 3.1 uses only a sufficient condition. A weaker necessary and sufficient one to guarantee $|c_{2}|\leq 1$ and $s_{2}\in[0,1]$ is $\frac{a+d}{2}-\sqrt{\left(\frac{a-d}{2}\right)^{2}+b^{2}}\leq\tau\leq\frac{a+d}{2}+\sqrt{\left(\frac{a-d}{2}\right)^{2}+b^{2}}$. then there exists one $\theta\in[-\pi/2,\pi/2]$ s.t. $a^{\prime}=\tau$, $d^{\prime}=a+d-\tau$ and $b^{\prime}=-s_{2}\sqrt{\left(\frac{a-d}{2}\right)^{2}+b^{2}}$ with $\cos(\theta)=\sqrt{\frac{1+c_{1}c_{2}-s_{1}s_{2}}{2}}$ and $\sin(\theta)=-\frac{c_{1}s_{2}+c_{2}s_{1}}{2\cos\theta}$, $c_{1}=\left(\frac{a-d}{2}\right)/\sqrt{\left(\frac{a-d}{2}\right)^{2}+b^{2}}$, $s_{1}=b/\sqrt{\left(\frac{a-d}{2}\right)^{2}+b^{2}}$, $c_{2}=(\tau-\frac{a+d}{2})/\sqrt{\left(\frac{a-d}{2}\right)^{2}+b^{2}}$ and $s_{2}=\sqrt{1-c_{2}^{2}}\in[0,1]$.

Proof 3.1

*With $u(\theta)\overset{def}{=}\begin{pmatrix}\cos(2\theta)&-\sin(2\theta)\end{pmatrix}^{T}$, Eq. 3 gives: $a^{\prime}=\frac{a+d}{2}+\begin{pmatrix}\frac{a-d}{2}&b\end{pmatrix}.\ u(\theta)^{T}$,

$d^{\prime}=\frac{a+d}{2}-\begin{pmatrix}\frac{a-d}{2}&b\end{pmatrix}.\ u(\theta)^{T}$ and $b^{\prime}=\begin{pmatrix}\frac{a-d}{2}&b\end{pmatrix}.\begin{pmatrix}\sin(2\theta)\\ \cos(2\theta)\end{pmatrix}$. As the Givens angle $\theta$ should be parameterized s.t. all diagonal coefficients are set to $\tau$, $\begin{pmatrix}\frac{a-d}{2}&b\end{pmatrix}.\ u(\theta)^{T}=\tau-\frac{a+d}{2}$. Previously defined $c_{1}$, $c_{2}$, $s_{1}$, $s_{2}$ and the condition $\operatorname{min}(a,d)<\tau<\operatorname{max}(a,d)$ gives $|c_{2}|\leq 1$ and $s_{2}\in[0,1]$. Then the last equation becomes: $\begin{pmatrix}c_{1}&s_{1}\end{pmatrix}.\ u(\theta)^{T}=c_{2}$. A clear solution is: $\ u(\theta)^{T}=\begin{pmatrix}c_{1}c_{2}-s_{1}s_{2}\\ c_{1}s_{2}+c_{2}s_{1}\end{pmatrix}.$ Thus, one can take: $\cos(\theta)=\sqrt{\frac{1+\cos(2\theta)}{2}}=\sqrt{\frac{1+c_{1}c_{2}-s_{1}s_{2}}{2}}$, $\sin(\theta)=\frac{\sin(2\theta)}{2\cos(\theta)}=-\frac{c_{1}s_{2}+c_{2}s_{1}}{2\cos\theta}$ and the associated Givens rotation gives: $a^{\prime}=\tau$,

$d^{\prime}=a+d-\tau$,   $b^{\prime}=\begin{pmatrix}\frac{a-d}{2}&b\end{pmatrix}\begin{pmatrix}\sin(2\theta)\\ \cos(2\theta)\end{pmatrix}$

$=\sqrt{\left(\frac{a-d}{2}\right)^{2}+b^{2}}\begin{pmatrix}c_{1}&s_{1}\end{pmatrix}\begin{pmatrix}-c_{1}s_{2}-c_{2}s_{1}\\ c_{1}c_{2}-s_{1}s_{2}\end{pmatrix}$

$=-s_{2}\sqrt{\left(\frac{a-d}{2}\right)^{2}+b^{2}}$. *

Note that there is no need to compute explicitly $\theta$, $\theta_{1}$ or $\theta_{2}$. We now briefly describe the underlying diagonal uniformization algorithm. The mean of $\Sigma_{V}$ diagonal coefficients being equal to $\tau$, these indices sets are not empty: $iInf\overset{def}{=}\{\ l\in\{1,\ ...,\ c\}\ |\ \Sigma_{V,l,l}<\tau\}$ and $iSup\overset{def}{=}\{\ l\in\{1,\ ...,\ c\}\ |\ \Sigma_{V,l,l}>\tau\}$. Taking one index $j$ from $iInf$ and the other one $i$ from $iSup$ guarantees the condition of Th. 3.1, which allows to set $\Sigma_{V,j,j}$ to the value $\tau$. The index $j$ can then be removed from $iInf$ and as $\Sigma_{V,i,i}$ is set to $a+d-\tau$, the index $i$ reassigned to $iInf$ if $\tau>\frac{a+d}{2}$ or in $iSup$ if $\tau<\frac{a+d}{2}$. The number of diagonal coefficients of $\Sigma_{V}$ different from $\tau$ has been decreased by one. Finally, the necessary number of iterations to completely empty $iInf$ and $iSup$, i.e. uniformizing $\Sigma_{V}$ diagonal, is bounded by $c-1$. The method is summarized in Algorithm 1 where $\operatorname{pop}(list)$ and $\operatorname{add}(list,e)$ are subroutines to delete and return the first element of $list$, resp. to add $e$ in $list$.

4 Complexity analysis of existing works

Our algorithm requires the storage of two $c\times c$ matrices besides $W_{t}$ and $R_{t}$ obviously: one with OPAST to obtain $W_{t}$ and $\Sigma_{V,t}$ for $R_{t}$. One update with OPAST for $W_{t}$ and $\Sigma_{V,t}$ costs $4dc+O(c^{2})$. Then, to compute $R_{t}$, at most $c-1$ Givens rotations are needed, each implying four column or row multiplications i.e. $4c$ flops. So the final time complexity of our algorithm is $4dc+O(c^{2})$.

We compare here the spatial and time costs of our method with, to the best of our knowledge, the only online unsupervised method: Online Sketching Hashing (OSH) [11] which is the most similar to ours, i.e. unsupervised, hyperplanes-based and reading one data point at a time. Despite what is announced, OSH is fundamentally mini-batch: the stream is divided into chunks of data for which a matrix $S\in\mathbb{R}^{d\times l}$ as a sketch of the whole dataset $X\in\mathbb{R}^{d\times n}$ is maintained. Then the principal components are computed from the updated sketch $S$. The projection of data followed by the random rotation can be applied only after this step. Therefore there are actually two passes over the data by reading twice data of each chunk. Without counting the projection matrix and the rotation, OSH needs spatially to maintain the sketch $S$ which costs $O(d\times l)$ with $c\ll l\ll d$. The SVD decomposition then needs $O(dl+l^{2})$ space. In comparison, we only need $O(c^{2})$. For each round, OSH takes $O(dl^{2}+l^{3})$ time to learn the principal components, i.e. $O(dl+l^{2})$ for each new data seen. We also compare with IsoHash [8]. Although it counts as an offline method because no technique is proposed to approximatively estimate the principal subspace, IsoHash rotation can be applied after for instance OPAST. IsoHash rotation computation involves an integration of a differential equation using Adams-Bashforth-Moulton PECE solver which costs $O(c^{3})$ time. Even if $c$ is small in comparison to $d$ and the complexities do not either depend on $n$, our model has the advantages to have a lower time cost and to be much more simple than IsoHash. Thus, our method shows advantages in terms of spatial and time complexities over OSH and IsoHash. Moreover, binary hash codes can be directly computed as new data is seen, while OSH, as a mini-batch method, has a delay.

5 Experiments

Experiments are made on CIFAR-10 (http://www.cs.toronto.edu/~kriz/cifar.html) and GIST1M sets (http://corpus-texmex.irisa.fr/). CIFAR-10 contains $60000$ $32\times 32$ color images equally divided into 10 classes. 960-D GIST descriptors were extracted from data. GIST1M contains $1$ million $960$-D GIST descriptors, from which $60000$ instances were randomly chosen from the first half of the learning set. Quality of hashing has been assessed on the nearest neighbor (NN) search task performed on the binary codes instead of the initial descriptors. A nominal threshold of the average distance to the $50$th nearest neighbor is computed and determines the sets of neighbors and non-neighbors called Euclidean ground truth. $1000$ queries were randomly sampled and the remaining data are used as training set. We compared our method to three online baseline methods that follow the basic hashing scheme $\Phi(x_{t})=\operatorname{sgn}(\tilde{W}_{t}x_{t})$, where the projection matrix $\tilde{W}_{t}\in\mathbb{R}^{c\times d}$ is determined according to the chosen method:

1. OSH 2) RandRot-OPAST: $W_{t}$ is the PCA matrix obtained with OPAST and $R_{t}$ a constant random rotation.
2. IsoHash-OPAST: $R_{t}$ is obtained with IsoHash.
3. UnifDiag-OPAST: $R_{t}$ is obtained with UnifDiag. For OSH, the number of chunks is set to $200$ and $l=50$. Fig. 1 and 2 (best viewed in color) show the Mean Average Precision (mAP) [29] for both datasets for $c=32$ (similar results are obtained for $c\in[8,16,64]$) averaged over $5$ random training/test partitions. Our algorithm outperforms all the compared methods. Code on GitHub: annemorvan/UnifDiagStreamBinSketching/.

6 Conclusion

We introduced a novel method for learning distance-preserving binary embeddings of high-dimensional data streams with convergence guarantees. Unlike classical state-of-the-art methods, our algorithm does not need to store the whole dataset and enables to obtain without delay a binary code as a new data point is seen. Our approach shows promising results as evidenced by the experiments. It can achieve better accuracy than state-of-the-art online unsupervised methods while saving considerable computation time and spatial requirements. Besides, the Givens rotations, that are a classical tool for QR factorization, singular and eigendecomposition or joint diagonalization, can also be used for uniformizing the diagonal of a symmetric matrix via an original Givens angle tuning technique. Further work would be to investigate whether another rotation, not uniformizing the diagonal of the covariance matrix of the projected data, could be more optimal. Another interesting perspective is to evaluate the performance of the compact binary codes in other machine learning applications: instead of using the original data, one could use directly these binary embeddings to perform unsupervised or supervised learning while preserving the accuracy.