Ahmet Iscen
email: ahmet.iscen@inria.fr

Teddy Furon
email: teddy.furon@inria.fr

Vincent Gripon
email: vincent.gripon@telecom-bretagne.eu

Michael Rabbat
email: michael.rabbat@mcgill.ca

Hervé

Memory vectors for similarity search in high-dimensional spaces

Keywords: High-dimensional indexing, image indexing, image retrieval

come

INTRODUCTION

W E consider the problem of searching for vectors similar to a query vector in a large database. The typical applications are the search and exploration in Big Media Data where documents are represented by feature vectors [START_REF] Wang | Special issue on big media data: Understanding, search, and mining[END_REF]. In this context, many papers report how the curse of dimensionality (due to the size of the vectors) makes indexing techniques ineffective [START_REF] Weber | A quantitative analysis and performance study for similarity-search methods in high-dimensional spaces[END_REF], [START_REF] Muja | Scalable nearest neighbor algorithms for high dimensional data[END_REF]. The recent paper [START_REF] Muja | Scalable nearest neighbor algorithms for high dimensional data[END_REF] describing and analyzing the popular FLANN method experimentally observes that even this state-of-the-art method performs poorly on synthetic high-dimensional vectors, and the authors conclude that "random datasets are one of the most difficult problems for nearest neighbor search". Some strategies have been proposed to (partly) overcome this problem. For instance, the vector approximation file [START_REF] Weber | A quantitative analysis and performance study for similarity-search methods in high-dimensional spaces[END_REF] first relies on exhaustive search with approximate measurements and then computes the exact similarities only for a subset of vectors deemed of interest. The cosine sketch [START_REF] Charikar | Similarity estimation techniques from rounding algorithms[END_REF] approximates cosine similarity with faster Hamming distance. Other works like spectral hashing [START_REF] Weiss | Spectral hashing[END_REF], Euclidean sketches [START_REF] Dong | Asymmetric distance estimation with sketches for similarity search in high-dimensional spaces[END_REF], product quantization [START_REF] Jégou | Product quantization for nearest neighbor search[END_REF] and inverted multi-index [START_REF] Babenko | The inverted multi-index[END_REF] also rely on compact codes to speed up neighbor search while compressing the data. An interesting strategy is the Set Compression Tree [START_REF] Arandjelović | Extremely low bit-rate nearest neighbor search using a set compression tree[END_REF], which uses a structure similar to a k-d tree to compress a set of vectors into an extremely compact representation. Again, this method is dedicated to small dimensional vectors (its authors recommend the dimension be smaller than log 2 (N ) where N is the size of the database) so that it is used in conjunction with a drastic dimension reduction via PCA to work with classical computer vision descriptors.

The main contribution of the paper is a similarity search approach specifically adapted to high-dimensional vectors such as those recently introduced in computer vision to represent images [START_REF] Perronnin | Fisher kernels on visual vocabularies for image categorization[END_REF], [START_REF] Jégou | Aggregating local descriptors into a compact image representation[END_REF]. The proposed indexing architecture consists of memory units, each of which is associated with several database vectors. A representative, called a memory vector, is produced for each memory unit and defined such that one can quickly and reliably infer whether or not at least one similar vector is stored in this unit by computing a single inner product with a new query.

Our approach is similar to the descriptor pooling problem in computer vision, but at a higher level. Many successful descriptors, such as BOV [START_REF] Csurka | Visual categorization with bags of keypoints[END_REF], [START_REF] Sivic | Video Google: A text retrieval approach to object matching in videos[END_REF], VLAD [START_REF] Jégou | Aggregating local descriptors into a compact image representation[END_REF], FV [START_REF] Perronnin | Fisher kernels on visual vocabularies for image categorization[END_REF], and EMK [START_REF] Bo | Efficient match kernels between sets of features for visual recognition[END_REF], encode and aggregate a set of local features into global representations. Yet, the new representation has a larger dimension than the local features. We use a similar approach at a higher-level: instead of aggregating local features into one global image representation, we aggregate global representations into group representations, so called memory vectors. These have no semantical meaning. Their purpose is to allow efficient search. Another difference is that we keep the same ambient dimension.

A second contribution of the paper is to study similarity search from the perspective of statistical signal processing and decision theory. Our analysis provides insight into when and why the proposed approach provides low complexity (fast response time) without sacrificing accuracy. A third contribution is the experimental work using computer vision Big Media datasets as large as 100 million images. Our results suggest that the approach we propose can achieve accuracy comparable to state-of-the-art while providing results 5-10× faster.

The paper is organized is as follows. Section 2 gives a formal problem statement. Section 3 focuses on the design of a single memory vector. We formalize the similarity of a query with the vectors of one memory unit as a hypothesis test. We derive the optimal representative vector under some design constraints and show how to compute it in an online manner. Section 4 proposes and analyzes two different ways to assign database vectors to memory units: random assignment and weakly supervised assignment which packs similar vectors into memory units. We provide a theoretical and experimental analysis of the different design and assignment strategies. Section 5 presents the results of experiments that evaluate our approach on standard benchmarks for image search. We use descriptors (vectors describing images) extracted with the most recent state-of-the-art algorithm in computer vision [START_REF] Jégou | Triangulation embedding and democratic kernels for image search[END_REF]. Our results show the potential of our approach for this application. Finally, we present the complete analysis of our method in Section 6.

PROBLEM FORMULATION

Let X = {x 1 , . . . , x n } denote a collection of d-dimensional vectors. We assume that all vectors are normalized so that x i 2 2 = 1 for all i = 1, . . . , n.

Given a query vector y ∈ R d with y 2 2 = 1, our objective is to determine which vectors x i ∈ X are most similar to y. Determining which images would be "most similar" to a human is subjective and not easily quantified. The real image datasets used in the experiments we report in Section 5 include a set of queries and the corresponding human-determined response sets, which are treated as ground truth for those experiments. While humans can provide result sets for individual queries, generalizing to produce result sets for never-seen-before queries is still an open problem.

More generally, when such a ground truth is not given to measure performance, since all vectors lie on the unit sphere, we measure similarity using the inner product x i ⊤ y. Then the problem of determining the most similar vectors to the query can be stated more precisely as: given α 0 > 0, find all vectors x i ∈ X such that x i ⊤ y ≥ α 0 . The accuracy of a technique can be measured in terms of its precision and recall.

Clearly a naive baseline approach to this problem would compute all n inner-products x i ⊤ y. Although it may provide perfect accuracy, this approach would have computational complexity (strongly related to run-time) of O(nd) operations. This is generally unacceptable when n and/or d are large, and we aim to obtain high accuracy while reducing the complexity.

MEMORY VECTORS

Given a memory unit X = {x 1 , . . . , x n } ⊂ R d , our objective in this section is to produce a representative, a so-called memory vector, such that, given a query vector Y regarded as a random variable, we can efficiently perform a "similarity" test answering: is Y a quasi-copy of, or similar to, at least one of the vectors of the memory unit?

For the sake of analysis, this section assumes that the vectors x i are independent and identically distributed samples from a uniform distribution on the d-dimensional unit hypersphere. We model the query as a random vector Y distributed according to one of the two laws:

• Hypothesis H 0 : Y is not related to any vector in X . Y is then uniformly distributed on the unit hypersphere.

• Hypothesis H 1 : Y is related to one vector in X , say x 1 without loss of generality. We write this relationship as Y = αx 1 + βZ, where Z is a random vector orthogonal to x 1 and Z = 1. This means that Y is more similar to x 1 as α gets closer to 1. We have

α 2 + β 2 = 1 because Y = 1.
We look for a representation scheme satisfying the following design constraints. First, the set of vectors X is summarized by a single vector of the same dimension, m(X ) ∈ R d , called the memory vector and denoted by m when not ambiguous. Second, the potential membership of query Y to X is tested by thresholding the inner product m ⊤ Y.

Sum-memory vector: Analysis

A very simple way to define the memory vector is

m(X ) = x∈X x, (1)
where we assume that X is composed of n different vectors. Albeit naive, this strategy offers some insights when considering highdimensional spaces. Section 6.1 derives the pdf of the score m ⊤ Y under H 0 when m is a known vector. This score has expectation 0 and variance m 2 /d, and it is asymptotically distributed as N (0, m 2 /d) as d → ∞. This gives an approximate pdf of m ⊤ Y under H 0 . In contrast, under H 1 , the inner product equals

m ⊤ Y = α + αm(X ′ ) ⊤ x 1 + βm(X ′ ) ⊤ Z, (2)
with X ′ = X -{x 1 }. This shows two sources of randomness: the interference of x 1 with the other vectors in X and with the noise vector Z. Assuming that Y is statistically independent of the vectors in X ′ (this implies that the vectors of X are mutually independent), we have

E Y [m ⊤ Y|H 1 ] = α, V[m ⊤ Y|H 1 ] = m(X ′ ) 2 /d. (3)
Assuming that X is composed of n < d statistically independent vectors on the unit hypersphere also gives E

X [ m(X ) 2 ] = n and E X ′ [ m(X ′ ) 2 ] = n -1.
To summarize, for large d, we expect the following distributions:

H 0 : m ⊤ Y ∼ N (0, n/d), (4)
H 1 : m ⊤ Y ∼ N (α, (n -1)/d). (5)
Making a hard decision by comparing the inner product to a threshold τ , the error probabilities (false positive and false negative rates) are given by

P fp ≈1 -Φ τ d/n (6)
P fn ≈Φ (τ -α) d/(n -1) , (7)
where

Φ(x) = 1 √ 2π x -∞ e -t 2 /2 dt.
The number of elements one can store in a sum-memory vector is linear with the dimension of the space when vectors are drawn uniformly on the unit hypersphere. This construction is therefore useful for high-dimensional vectors only, as opposed to traditional indexing techniques that work best in low-dimensional spaces.

If the vectors were pair-wise orthogonal, the dominant source of randomness (the interference between x 1 and vectors of X ′ ) is cancelled in [START_REF] Weber | A quantitative analysis and performance study for similarity-search methods in high-dimensional spaces[END_REF]. The variance under H 1 reduces to β 2 (n -1)/d. This prevents any false negative if β → 0. We further exploit this intuition that orthogonality helps in the next section.

Optimization of the hypothesis test per unit

We next consider optimizing the construction of the memory vector of a given set X . Denote the d × n matrix X = [x 1 , . . . , x n ]. We impose that, for all i, x i ⊤ m(X ) = 1 exactly and not only in expectation, as assumed above. In other words, X ⊤ m = 1 n where 1 n is the length-n vector with all entries equal to 1. Achieving this, when Y = x 1 , we eliminate the interference with the remaining vectors in X ′ which was previously the dominant source of noise. In other words, under H 1 , Eq. ( 2) becomes

m ⊤ Y = α + βm ⊤ Z. (8)
Under H 0 , the variance of the score remains m 2 /d. Therefore, the norm of the memory vector is the key quantity determining the false positive probability.

We thus seek the representation m minimizing the energy m 2 subject to the constraint that X ⊤ m = 1 n . If multiple solutions exist, the minimal norm solution is given by the Moore-Penrose pseudo-inverse [START_REF] Rao | Generalized inverse of a matrix and its applications[END_REF]:

m ⋆ = (X + ) ⊤ 1 n . (9)
Since n < d, m ⋆ = X(X ⊤ X) -1 1 n . If no solution exists, m ⋆ is a minimizer of X ⊤ m -1 n . This formulation amounts to treating the representation of a memory unit as a linear regression problem [START_REF] Bishop | Pattern recognition and machine learning[END_REF] with the objective of minimizing over m the quantity

X ⊤ m -1 n 2 .
Taking the gradient, setting it equal to zero, and solving for m gives back m ⋆ . When possible, and for large d, this new construction leads to the distributions:

H 0 : m ⋆⊤ Y ∼ N (0, m ⋆ 2 /d) (10) H 1 : m ⋆⊤ Y ∼ N (α, β 2 m ⋆ 2 /d). (11)
The major improvement comes from the reduction of the variance under H 1 for small values of β 2 , i.e., α 1. Section 6.2 shows that if the vectors of X are uniformly distributed then m ⋆ 2 is larger in expectation than the square norm of the naive sum representation from Section 3.1. The reduction of the variance under H 1 comes at the price of an increase of the variance under H 0 . However, this increase is small if n/d remains small. For large d, we have

P fp ≈1 -Φ τ d n -1 , (12)
P fn ≈Φ τ -α β d n -1 . ( 13
)
Note that β = √ 1α 2 is a decreasing function of α. Therefore, if τ < α, P fn is a decreasing function of α. In particular P fn → 0 when α → 1 as claimed above. In contrast to the naive sum approach from Section 3.1, there is no longer false negative when the query Y is exactly one of the vectors in X . This holds for any value of τ < 1 when α = 1, so that the false positive rate can be as low as 1 -Φ( d/n -1). Remark. This solution is identical (up to a regularization) to the "generalized max-pooling" method introduced to aggregate local image descriptors [START_REF] Murray | Generalized max-pooling[END_REF]. However in our case the aggregation is performed on the database side only. Our solution is moreover theoretically grounded by a hypothesis test interpretation. . The different curves correspond to values of α 0 ∈ {0.5, 0.6, 0.7, 0.8, 0.9}. Red and blue lines correspond to sum and pinv constructions respectively. Darker shades correspond to higher α 0 .

Weakly supervised assignment

We now analyze a scenario where the vectors packed in the same memory unit are random but no longer uniformly distributed over the hypersphere: There is some correlation among them. The vectors in a memory unit are now uniformly distributed over a spherical cap (see Sections 6.4 and 6.5). This models the effect of a pre-processing which analyses the database vectors in order to assign similar vectors to memory units. For instance, Section 4.1 uses the k-means algorithm to process batches of database vectors.

We derive the same analysis as in the previous subsection with expectations and variances which now depend on the angle of the spherical cap. These expressions are complex and their derivation is detailed in Section 6. In summary, the Kullback-Leibler distance between the distributions of m ⊤ Y under both hypotheses increases as the spherical cap gets narrow. In other words, identifying the positive memory units becomes easier when we assign correlated vectors to the same memory unit. Interestingly, this mechanism helps the sum construction more than pinv, so that when the vectors are very correlated, both constructions indeed perform equivalently.

EXPERIMENTAL INVESTIGATION

We next consider application scenarios where we need to store a large number N of vectors and perform similarity search. One memory vector is not sufficient to achieve a reliable test. We therefore consider an architecture that consists of M memory vectors. The search strategy is as follows. A given query vector is compared with all the memory vectors. Then we compare the query with the vectors stored in the memory units associated with the high responses, i.e., those likely to contain a similar vector.

Experimental Setup

Our experimental investigations are carried out using synthetic data (vectors uniformly distributed over the hypersphere) as well as real data, which are described in this section. Datasets. We use the Inria Holidays [START_REF] Jégou | Hamming embedding and weak geometric consistency for large scale image search[END_REF], Oxford5k [START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF], and UKB [START_REF] Nistér | Scalable recognition with a vocabulary tree[END_REF] image datasets in our experiments. Additionally, we conduct large scale experiments in Holidays+Flickr1M, which is created by adding images from the Flickr1M [START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF] dataset to the Holidays dataset. Similarly we use the recently introduced Yahoo100M dataset [START_REF] Thomee | Yfcc100m: The new data in multimedia research[END_REF] to increase the dataset size.

Descriptors. We use the state-of-the-art triangular embedding descriptor [START_REF] Jégou | Triangulation embedding and democratic kernels for image search[END_REF], denoted by φ ∆ . We use the off-the-shelf reference implementation provided by the authors, which can be found online. 1 Each image is represented by a feature vector. The only difference is that we do not apply the "powerlaw normalization" to better illustrate the benefit of the pinv technique for the memory vector construction compared to the sum (when applying the powerlaw normalization, both designs perform equally well since the vectors are nearly orthogonal). Ultimately, we have d = 8064 (or d = 1920) dimensional feature vectors for each image, obtained by using a vocabulary of size 64 (resp. 16). For large experiments in Yahoo100M, we use d = 1024 VLAD descriptors as extracted by [START_REF] Spyromitros-Xioufis | A comprehensive study over VLAD and product quantization in large-scale image retrieval[END_REF].

We also experiment using deep learning features (d = 4096) provided by Babenko et al. [START_REF] Babenko | Neural codes for image retrieval[END_REF]. As explained in their paper, the performance for the UKB dataset drops with adapted features trained on the Landmarks dataset. Therefore, we use the original neural codes trained on ILSVRC for the UKB dataset, and the adapted features for Holidays and Oxford5k.

Random assignment

We suppose that the N vectors in the database are randomly grouped into M units of n vectors: N = nM . We aim at finding the best value for n. When the query is related to the database (i.e., under H 1 ), we make the following assumption: α 0 < α < 1, and we fix the following requirement: P fn < ǫ < 1/2. Since P fn is a decreasing function of α, we need to ensure that P fn (α 0 ) = ǫ. This gives us the threshold τ :

τ = µ H1 + σ H1 Φ -1 (ǫ), (14)
with µ H1 and σ H1 being the expectation and the standard de- viation of m j ⊤ Y under H 1 . Note that Φ -1 (ǫ) < 0 because ǫ < 1/2 so that τ < α. The probability of false positive equals 1 -Φ(τ /σ H0 ) which depends on n, denoted by P fp (n). This is indeed an increasing function for both memory vector constructions. Now, we decide to minimize the expectation of the total computational cost C H0 when the query is not related. We need to compute one inner product m j ⊤ y per unit, and then to compute n inner products x i ⊤ y for the units giving a positive detection. In expectation, there are M • P fp (n) such units, and so

C H0 = M + M • P fp (n) • n = N (n -1 + P fp (n)). (15)
The total cost is the sum of a decreasing function (n -1 ) and an increasing function (P fp (n)).

For the random assignment strategy, there is a tradeoff between having a few big units (n large) and many small units (n small). Fig. 1 illustrates this tradeoff for different values of α 0 with synthetic data. It is not possible to find a closed form expression for the cost minimizer n ⋆ . When α 0 is close to 1, the threshold is set to a high value, producing reliable tests, and we can pack many vectors into each unit: n ⋆ is large allowing a huge reduction in complexity. Even when α 0 is as small as 0.5, n ⋆ is small but the improvement remains significant. In the setup of Fig. 1, the proposed approach has a complexity that is less than one tenth of that of searching through all database vectors (equivalent to n = 1). However, in order to increase the efficiency, we introduce an additional O(M d) = O(dN/n) memory overhead for storing memory vectors.

Figure 2 depicts the theoretical and empirical Receiver Operating Characteristic (ROC) curves for different values of α. For 1. http://www.tinyurl.com/democratic-kernel the synthetic data, P fn and P fp are evaluated using Eq. ( 6) and [START_REF] Csurka | Visual categorization with bags of keypoints[END_REF]. As expected the test performs better when α is closer to 1 and when n ≪ d. For the real data, we use cosine similaritybased ground truth, since it is directly related with the model we considered theoretically. For each query vector, we deem a database vector as relevant if their cosine similarity is greater than α 0 . To have enough ground-truth vectors, we look for these matching vectors on the Holidays+Flickr1M dataset using various α 0 values. This experiment using real data confirms the findings of the theoretical analysis. The pinv construction performs better than the sum as long as α 0 is big as explained in Sect. 3.2. However, the theoretical analysis is unable to predict performance levels on real dataset. It seems that the vectors of this real dataset have a much lower intrinsic dimensionality than their representational dimension d = 1920.

Weakly supervised assignment

A well-known technique in the approximate search literature is to partition the space R d by clustering the database vectors. This assigns similar database vectors to the same cell [START_REF] Muja | Scalable nearest neighbor algorithms for high dimensional data[END_REF]. In Section 3.3, we explain the advantage of a weakly supervised assignment by showing that the distance between the distributions of positive and negative memory units similarities increases. To show this point experimentally, we modify the spherical k-means clustering [START_REF] Dhillon | Concept decompositions for large sparse text data using clustering[END_REF], so that the clusters are represented using pinv (or sum) in the update stage 2 . Better hypothesis test. Figure 4 shows that highly ranked memory units are very likely true positives containing at least one matching vector. On the contrary, with the random assignment, a positive memory unit may have a low rank. This means that we now can analyse the database vectors of a shorter list of memory units to find most of the matching vectors. More than one match. Another byproduct of weakly supervised assignment is that positive memory units are likely to contain more than one match, since matching vectors usually have high cosine similarity with each other. This helps the search efficiency by returning most of the matching vectors by only scanning a few positive memory units. This is also experimentally shown in Figure 5. With the random assignment, we have almost surely at most one matching vector in each positive memory unit. Imbalance factor. We now analyze the cost of search with weakly supervised assignment. In Eq. ( 15), we assume that each unit contains n vectors. Up to this approximation, Fig. 3 shows that both constructions sum and pinv perform better as the inner correlation increases, but more surprisingly, they perform equivalently. It is also shown that it is possible to pack many vectors into the same memory unit with weakly supervised assignment, and still obtain a low search cost.

In practical applications, assuming that each unit contains a constant number of vectors is no longer true with weakly supervised assignment. This makes the analysis of the complexity more involved than [START_REF] Jégou | Triangulation embedding and democratic kernels for image search[END_REF]. Moreover, this is potentially problematic in some applications: the complexity and thus the runtime can change dramatically from one query to another. Imbalance factor is a metric to measure the impact of unbalanced clusters [START_REF] Jégou | Improving bag-of-features for large scale image search[END_REF]. It is defined as

δ = M M i=1 p i 2 , ( 16
)
where M is the number of clusters, and p i is the empirical probability that a database vector belongs to the i-th cluster. This is measured as frequency p i = n i /N , where n i is the cardinality of the i-th cluster. Simple derivations give the following expectation and variance:

E(n i ) = N/M and V(n i ) = (δ -1)N 2 /M 2 .
This shows that higher imbalance factor corresponds to clusters with varying sizes. This gives birth to a wide variability of the complexity from one query to another. Table 1 shows the imbalance factor for different weakly supervised assignments. It is shown that pinv variants have more balanced cells compared to traditional sum, making the search process more effective. The negative effect of high imbalance factor in practice is better observed in Figure 6. In this figure, the algorithm visits a fixed number of positive memory units: 7 (Holidays), 30 (Oxford5k), or 60 (UKB). This roughly gives us a complexity ratio of C H0 (τ ) ≈ 0.2 on average. We then show the complexity ratio per query in a histogram. It is clearly seen that the distribution for pinv has smaller standard deviation compared to sum, even though their means are almost the same. This makes pinv variant of spherical k-means a better alternative for weakly supervised assignment.

APPLICATION TO IMAGE SEARCH

This section shows that memory vectors perform extremely well on typical computer vision benchmarks. We assume two scenarios: closed-datasets in Section 5. Whereas the datasets were already introduced in Sect. 4.1, let us describe the measure of performances. We follow the standard image retrieval protocol where each image is represented by a feature vector and the ground truth is now based on the visual similarity. The goal is to return visually similar images for a given query image. The similarity of two images is measured by the cosine of their descriptor vectors, and the images are ordered accordingly. We adopt the performance measure defined for each benchmark: mAP (mean average precision) which measures the area under the precision-recall curve [START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF] or 4-recall@4, which is the average number of correct images in the top-4 positions.

As for the complexity, we first measure the similarities between the query and M memory vectors. We compare these similarities with a given threshold τ . Then, we re-rank all the vectors in positive memory units according to their similarities with the query vector. We characterize the complexity of the search per database vector by:

C H1 (τ ) = M + i:y ⊤ mi>τ n i , (17)
where n i is the number of database vectors in the i-th memory unit. We measure the complexity ratio C H1 (τ )/N and the retrieval performance for different values of the threshold τ . For large τ , no memory unit is positive, resulting in C H1 (τ )/N = M and no candidate is returned. As τ decreases, more memory units trigger reranking.

Closed dataset

Recall from Section 4 that weakly supervised assignment provides better approximate search than a random assignment. This is confirmed for the image search benchmark in Figure 7. Additionally, we show that it is possible to pack more vectors in a memory unit using weakly supervised assignment in Figure 8. We use this approach (spherical k-means with pinv) for the rest of our experiments.

The dimensionality of the descriptor linearly impacts the efficiency of any system. Dimensionality reduction with PCA is one way to improve this point. Our method is compatible with dimensionality reduction as shown in Fig. 9, where we reduce . Probability that a memory unit contains at least one match with respect to their rank. Green, red and blue lines correspond to random, sum spherical k-means and pinv spherical k-means respectively. There is a high probability of retrieving a match in highly ranked memory units, but it decreases faster with a weakly supervised assignment. 5. The average number of matches given that the memory unit is positive. Green, red and blue lines correspond to random, sum spherical k-means and pinv spherical k-means respectively. Using random assignment, we have about 1 matching vector per memory unit. The weakly supervised assignment improves this especially for higher-ranked units. the vectors to d ′ = 1024 components. The search performance is comparable to the baseline with less computational complexity. We also apply our method to features learned with deep learning (d = 4096). Fig. 10 shows that the reduction in complexity also applies when using high performance deep learning features.

Compact codes are another way to increase efficiency. We reduce the dimensionality of the vectors to d ′ =1024 and binarize them by taking the sign of each component, in the spirit of cosine sketches [START_REF] Charikar | Similarity estimation techniques from rounding algorithms[END_REF]. In the asymmetric case [START_REF] Gordo | Asymmetric distances for binary embeddings[END_REF], [START_REF] Jain | Asymmetric hamming embedding[END_REF], only memory vectors and dataset vectors are binarized, whereas in the symmetric case query vectors are also binarized during the query time.

Figures 11 and12 show the performance when using compact binary codes. For the symmetric case, the sum method seems to perform better than pinv on the Holidays and Oxford5k datasets. In the asymmetric case, both methods perform similarly. In all cases, we achieve convergence to the baseline with a complexity ratio well below 1. Implementation efficiency is further improved in the symmetric case by using the Hamming distance calculation instead of dot product. [START_REF] Muja | Scalable nearest neighbor algorithms for high dimensional data[END_REF]. Running the FLANN algorithm on the Holidays+Flickr1M dataset reveals that the convergence to the baseline is achieved with a speedup of 1.25, which translates Fig. 8. Search performance using random assignment (left) and weakly supervised assignment (right). This last option uses fewer memory units and still obtains a good search performance, which is not possible with random assignment. to a complexity ratio of around 0.8. We can achieve similar performance with a complexity ratio of only 0.3. This confirms that FLANN is not effective for high-dimensional vectors. In this experiment, we use the autotune option of the FLANN library, and set target precision = 0.95, build weight = 0.01, and memory weight = 0. Execution time. We have shown that we get close to baseline performance while executing significantly fewer operations. We now measure the difference in execution time under a simple setup: d = 1024 and N = 10 6 dataset vectors. An average dot product calculation between the query and all dataset vectors is 0.2728s. With N/10 memory vectors and ≈ 100k vectors in positive memory vectors, the execution time decreases to 0.0544s. We improve the efficiency even further with symmetric compact codes and Hamming Distance computation: the execution time becomes 0.0026s. Our method is parallelized for further improvement.

Comparison with FLANN

Large scale and streaming data

We conduct large scale experiments on Holidays+Flickr1M and Yahoo100M datasets. The main advantage of our approach is its compatibility with large scale and streaming data, where preclustering the data may not be possible. More specifically, we assume that we have streaming images which we would like to index. As the size of the data keeps growing continuously, it is not possible to apply traditional k-means in such a scenario. We investigate two different approaches: random assignment and weakly supervised assignment over mini-batches. Online indexing assumes that we would like to index items in streaming data as they become available. In such case, the random assignment is applicable provided that the successive vectors in the stream are independent.

Figure 13 shows the image retrieval performance based on random assignment with different group sizes n. With n = 10, the performance is close to the baseline while performing roughly three times fewer vector operations than exhaustive search. On the other hand, larger groups make it possible to have smaller complexity ratio with a degrading effect on the quality of search, since the scores obtained from memory vectors are noisier (see [START_REF] Dong | Asymmetric distance estimation with sketches for similarity search in high-dimensional spaces[END_REF] and ( 12)). The pinv construction performs better than sum in all cases except for a very large memory units of size n = 500, where the quality of search is low in general.

Batch assignment. The alternative approach runs weakly super- vised assignment on small batches. A batch spherical k-means is the same as the regular spherical k-means discussed in previous sections. We do not cluster the whole dataset at once because this would not be tractable with large collections. Instead, we randomly divide the dataset into batches of the same size and run a weakly supervised assignment separately for each batch. Fig. 14 shows that this strategy improves the performance while keeping the complexity of the clustering algorithm manageable. We compare our approach to a well-known mini-batch kmeans algorithm [START_REF] Sculley | Web-scale k-means clustering[END_REF] (referred to as mbk) as implemented by [START_REF] Pedregosa | Scikit-learn: Machine learning in python[END_REF]. We compare both strategies using a batch size of 10k, and show the image retrieval performance for different complexity ratio in Figure 15.

The first observation is that the plot of mbk is not smooth. This is due to the clusters being unbalanced when using such mini-batch approaches. In fact, when we measure the imbalance factor [START_REF] Rao | Generalized inverse of a matrix and its applications[END_REF], we obtain δ = 183 ± 8. As a result, few clusters contain a large number of dataset vectors, and when a cluster is accessed, the cost of the verification step becomes expensive. On the contrary, the imbalance factor observed using our batch spherical k-means is only 2.47±0.01, resulting in a more efficient verification step for positive memory units.

The complexity ratio per query for the two methods can be seen in Figure 16. When we set the number of positive memory units to 3500, the mean complexity ratio for our batch spherical k-means approach is 0.17, with a standard deviation of 0.01. On the other hand, when mbk is used, the mean increases to 0.26 with a standard deviation of 0.51.

Finally, we apply our batch assignment strategy to the Ya-hoo100M, which consists of 97.6 million vectors. We divide the dataset into batches of 100k, and run three different indexing strategies for each batch: random assignment with pinv, pinv-based spherical k-means and sum-based spherical k-means. This dataset does not have manually annotated ground truth or designated queries, therefore We use an existing evaluation protocol [START_REF] Iscen | Efficient large-scale similarity search using matrix factorization[END_REF]. Dataset vectors are considered a match if they have a similarity of 0.5 with the query vector. 1000 query vectors are randomly chosen and those that have 0 or more than 1000 matches are filtered out. In the end, we have 112 queries, with each query having 11.4 matches on average.

We present our results in Figure 17. Since we have N/10 memory vectors for each batch, a lower bound for the complexity ratio is 0.1. We see that pinv-based methods converge to the baseline much faster than sum. Using our pinv-based k-means variant, we achieve the same performance as baseline with 0.12 complexity ratio. This saves us almost 90 million vector operations at query time.

ANALYSIS

In this section, we provide the theoretical analysis of the methods presented in this paper.

Distribution of a scalar product

Let Y be a random vector uniformly distributed on the unit hypersphere in R d ( Y = 1), and m a fixed vector. This section studies the distribution of S = Y ⊤ m. To generate Y, we can first generate a multivariate Gaussian vector G = (G 1 , . . . , G d ) ⊤ ∼ N (0, I d ), where I d is d×d identity matrix, and set Y = G/ G . This means that G i are i.i.d. Gaussian distributed: G i ∼ N (0, 1). Without loss of generality, using symmetry of the Euclidean norm, assume that m = ( m , 0, . . . , 0). This simplifies into

S = Y ⊤ m = m G 1 d i=1 G 2 i . (18)
Obviously,m ≤ S ≤ m so that its cdf F S (s) equals 0 if s ≤m , and 1 if s ≥ m . For all s ∈ [0, m ], P(S ≥ s) is the probability that vector G belongs to the convex cone pointing in the direction m and with angle cos -1 (s/ m ) By a symmetry argument (replacing m by -m), P(S ≤ -s) = P(S ≥ s), giving It also implies that

F S (-s) = 1 -F S (s), ∀0 ≤ s ≤ m . ( 19
P(S 2 ≥ s 2 ) = P(S ≥ s) + P(S ≤ -s) = 2(1 -F S (s)). (20)
Going back to G, we can write

P(S 2 ≥ s 2 ) = P G 2 1 d i=1 G 2 i ≥ τ 2 = P G 2 1 d i=2 G 2 i ≥ τ 2 1 -τ 2 , (21)
with τ = s/ m . By definition, the random variable

U = (d - 1)G 2 1 / d i=2 G 2 i has an F-distribution F (1, d -1). It follows that P(S 2 ≥ τ 2 m 2 ) = 1 -I τ 2 (1/2, (d -1)/2), (22)
where I x (a, b) is the regularized incomplete beta function. In the end, ∀s ∈ [0, m ], we have

F S (s) = 1 + I s 2 m 2 1 2 , d -1 2 /2. ( 23
)
By taking the derivative of this expression and accounting for symmetry, the pdf is found to be:

f S (s) = (1 -s 2 / m 2 ) d-3 2 m B(1/2, (d -1)/2) , ∀s, -m ≤ s ≤ m
where B(a, b) is the beta function. This implies that E[S] = 0 and

V[S] = m 2 B 1 2 , d-1 2 1 0 a 2 (1 -a 2 ) d-3 2 da = m 2 B 1 2 , d-3 2 B 1 2 , d-1 2 = m 2 d .
When d → ∞, using the expansion [32, Eq. ( 26)] of the regularized incomplete beta function in [START_REF] Spyromitros-Xioufis | A comprehensive study over VLAD and product quantization in large-scale image retrieval[END_REF] yields

F S (s) ≈ Φ d -1 m 2 2s 1 + 1 -s 2 / m 2 , (24)
which is approximately Φ( d/ m 2 s) for small s, i.e., the cdf of a centered Gaussian r.v. with variance m 2 /d.

Expected value of the memory unit

We assume that the vectors x i are not related to each other so that, for n < d, X has n linearly independent columns. Then

X + = (X ⊤ X) -1 X ⊤ and m ⋆ 2 = 1 ⊤ n (X ⊤ X) -1 1 n .
The n × n Gram matrix X ⊤ X is real, symmetric, and positive semidefinite, and so it has an eigendecomposition UΛU ⊤ , where Λ is a diagonal matrix with non-negative coefficients {λ i } n i=1 and

UU ⊤ = U ⊤ U = I n . Moreover, tr(X ⊤ X) = n i=1 λ i = n because x ⊤ i x i = 1, for all i. The first construction m = X1 n has a square norm m 2 = 1 ⊤ n X ⊤ X1 n = 1 ⊤ n UΛU ⊤ 1 n . The second construction m ⋆ = X(X ⊤ X) -1 1 n has a norm m ⋆ 2 = 1 ⊤ n (X ⊤ X) -1 1 n = 1 ⊤ n UΛ -1 U ⊤ 1 n .
It is difficult to say anything more for a given X. However, if we consider X ⊤ X as a random matrix, then

E[(X ⊤ X) -1 ] -(E[X ⊤ X]) -1 = E[(X ⊤ X) -1 ] -I n is positive semi-definite [33]. This shows that E[ M ⋆ 2 ] ≥ E[ M 2 ] = n.
This means that the second construction increases the variance of the score under H 0 .

To make further progress, we resort to the asymptotic theory of random matrices, and especially to the Marcenko-Pastur distribution. Suppose that both d and n tend to infinity while remaining proportional: n = cd with c < 1. The matrix X ⊤ X can be thought of as an empirical covariance matrix of the d vectors which are the rows of X. These vectors are i.i.d. with bounded support components. Therefore, the marginal distribution of the n eigenvalues {λ i } n i=1 of X ⊤ X asymptotically follows the Marcenko-Pastur distribution: for all λ such that

(1 - √ c) 2 ≤ λ ≤ (1 + √ c) 2 , f MP (λ) = (λ -(1 - √ c) 2 )((1 + √ c) 2 -λ) 2cπλ . (25)
Moreover, for any function ψ bounded over the interval

[(1 - √ c) 2 , (1 + √ c) 2 ]: 1 n n i=1 ψ(λ i ) -ψ(λ)f MP (λ)dλ → 0. (26)
Using the eigendecomposition, X ⊤ X = UΛU ⊤ , we get (X ⊤ X) -1 = UΛ -1 U ⊤ , where the columns of U are the random eigenvectors. Then, asymptotically as d → ∞:

n -1 E[1 ⊤ n (X ⊤ X) -1 1 n ] -λ -1 f MP (λ)dλ → 0. (27)
This shows that E[ M ⋆ 2 ]/n converges to 1 1-c asymptotically (the above integral is the Stieljes transform of the Marcenko-Pastur distribution evaluated at z = 0), whereas E[ M 2 ]/n converges to E[λ] = 1. In expectation, M ⋆ has a higher variance, but only by a factor of 1/(1c) which remains acceptable if c = n/d is small.

Y uniformly drawn over a spherical cap

Assume that Y is uniformly distributed over the spherical cap C u,γ , which is the intersection of the unit hypersphere and the single hypercone of axis u, u = 1 and angle γ. In other words, Y = 1 and S ′ = Y ⊤ u > cos(γ). Denote η = cos(γ) and η = sign(η). The probability distribution function of S ′ is

f S ′ (s ′ ) = f S (s ′ ) 1 -F S (s ′ ) 1 [s ′ >η] (s ′ ). (28)
This stems into

E[S ′ ] = 2(1 -η 2 ) d-1 2 (d -1)B( 1 2 , d-1 2 ) 1 -ηI η 2 ( 1 2 , d- 1 2 )
.

Note that:

• η = -1: E[S ′ ] = 0.
The cap is the full hypersphere.

• η → 1: E[S ′ ] → 1, thanks to De l'Hospital's rule. The cap reduces to {u}.

In the same way:

E[S ′2 ] = 1 d 1 -ηI η 2 ( 3 2 , d-1 2 ) 1 -ηI η 2 ( 1 2 , d-1 2 ) , (30)
from which we can deduce V(S ′ ) with the König-Huygens formula. Note that:

• η = -1: V[S ′ ] = 1/d. • η → 1: V[S ′ ] → 0, thanks to De l'Hospital's rule.
From now on, we define

µ κ (η, d) = E[(S ′ ) κ ], (31)
with S ′ = Y ⊤ u and Y ∼ U Cu,γ .

Modeling k-means for the sum construction

We assume that k-means has packed together in a memory unit independent vectors uniformly distributed over the spherical cap:

x i ∼ U Cu,γ . Vector x i can be modeled as x i = S ′ i u + N i , where {S ′ i } n i=1 are i.i.d. according to the pdf described in Section 6.3,

N ⊤ i u = 0 and N i 2 = 1 -(S ′ i ) 2 .
Their memory vector is the sum m = n i=1 x i .

Hypothesis H 0

Y is independent from m. It follows that:

E[Y ⊤ m] = 0, V[Y ⊤ m] = E[ m 2 ]/d (32)
with

E[ m 2 ] = E[ n i=1 x i 2 + 2 i<j S ′ i S ′ j + N ⊤ i N j ] (33) = n + n(n -1)µ 1 (η, d) 2 , (34)
where µ 1 (η, d) is given in [START_REF] Sculley | Web-scale k-means clustering[END_REF].

To sum up: while E[Y ⊤ m] = 0, the variance is increasing with η.

• η = -1: It equals n/d as already shown by (4).

• η → 1: It converges to n 2 /d.

Hypothesis H 1

We now single out the role of the matching vector X 1 : m = x 1 + m ′ . The query is modeled as Y = αx 1 + βZ with Z = 1 and Z ⊤ x 1 = 0. With the same notation as above:

Y ⊤ m = α + α n i=2 S 1 S i + αN ⊤ 1 m ′ + βZ ⊤ m ′ . (35)
The expectation easily comes as

E[Y ⊤ m] = α(1 + (n -1)µ 1 (η, d) 2 ). (36)
The variance of the second summand is given by the law of total variance:

V[α n i=2 S 1 S i ] = α 2 (n -1) µ 2 (η, d) 2 -µ 1 (η, d) 4 . ( 37
)
The variance of the third term is

V[α n i=2 N ⊤ 1 N i ] = α 2 (n -1) d -1 (1 -µ 2 (η, d)) 2 . ( 38
)
The variance of the last summand is more complex to analyze. We need to decompose Z into its projection on u and on the complementary space.

V[β(Z ⊤ u)u ⊤ m ′ ] = (1 -α 2 ) n -1 d -1 (1 -µ 2 (η, d))µ 2 (η, d) V[βZ ⊤ ⊥ m ′ ] = (1 -α 2 ) n -1 d -1 (1 -µ 2 (η, d)) ×(1 + (n -2)µ 1 (η, d) 2 ) To sum up: E[Y ⊤ m] increases while V[Y ⊤ m] decreases with η. • η = -1: E[Y ⊤ m] = α while V[Y ⊤ m] = (n -1)/d as already shown by (5). • η → 1: E[Y ⊤ m] → nα whereas V[Y ⊤ m] → 0.
In the end, under the Gaussian assumption, the Kullback-Leibler distance between both distributions increases which proves that identifying the positive memory units becomes easier as η increases (see Fig. 18).

Modeling k-means for the pinv construction

Hypothesis H 0

We make the same assumption as in Section 6.4. We write X as X = uS ′⊤ + N with S ′ a n × 1 vector storing the correlations u ⊤ x i , 1 ≤ i ≤ n and N a d × n whose columns are random vectors orthogonal to u and of norm 1 -S ′2 i . The memory unit is now given by ( 9): m ⋆ = X(X ⊤ X) -1 1 n . Eq. ( 32) holds but with a new expression for the norm of m ⋆ :

E[ m ⋆ 2 ] = E[1 ⊤ n (X ⊤ X) -1 1 n ] (39) ≥ 1 ⊤ n (E[X ⊤ X]) -1 1 n . (40)
We write matrix

X ⊤ X = X ⊤ X = S ′ S ′⊤ + NN ⊤ s.t.: E[X ⊤ X] = (1 -E[S ′ ] 2 )I n + E[S ′ ] 2 1 n 1 ⊤ n , (41)
whose inverse is given by the Sherman-Morrison formula: leading to

E[X ⊤ X] -1 = 1 1 -E[S ′ ] 2 I n + E[S ′ ] 2 1 + (n -1)E[S ′ ] 2 1 n 1 ⊤ n , -
E[ m ⋆ 2 ] ≥ n 1 + (n -1)E[S ′ ] 2 (42)
To sum up: While E[Y ⊤ m ⋆ ] remains constant, the lower bound of the variance is decreasing with η. • η → 1: The lower bounds converges to 1/d and is also tight since m ⋆ → x 1 as the vectors are converging to x 1 .

Hypothesis H 1

We single out the matching vector x 1 by writing X = (x 1 , X) with X = uS ′⊤ + N, S ′ being now a (n -1) × 1 vector storing the correlations u ⊤ x i , 2 ≤ i ≤ n and N a d × (n -1) matrix whose columns are random vectors orthogonal to u and of norm 1 -S ′2 i . This makes

X ⊤ X = 1 x ⊤ 1 X X⊤ x 1 X⊤ X , (43)
whose inverse is

(X ⊤ X) -1 = 1 + x ⊤ 1 XD X⊤ x 1 -x ⊤ 1 XD -D X⊤ x 1 D . ( 44
)
with D = ( X⊤ (I -X 1 X ⊤ 1 ) X) -1 . This makes the following memory vector:

m ⋆ = x 1 + m ⋆ ⊥ (45) m ⋆ ⊥ = (I -x 1 x ⊤ 1 ) XD(1 n-1 -Xx 1 ) (46)
Note that m ⋆ 2 = 1 + m ⋆ ⊥ 2 because x ⊤ 1 m ⋆ ⊥ = 0. The query vector being defined as in Sect. 6.4.2, its correlation with the memory vector is whose expectation and variance are given by

Y ⊤ m ⋆ = α + 1 -α 2 Z ⊤ m ⋆ ⊥ , (47)
E[Y ⊤ m ⋆ ] = α, V[Y ⊤ m ⋆ ] = (1 -α 2 ) E[ m ⋆ ⊥ 2 ] d -1 (48)
with

E[ m ⋆ ⊥ 2 ] = E[ m ⋆ 2 ] -1 (49)
≥ (n -1)(1 -E[S ′ ] 2 ) 1 + (n -1)E[S ′ ] 2 . ( 50
)
To sum up: While E[Y ⊤ m ⋆ ] remains constant, the lower bound of the variance is decreasing with η.

• η = -1: The lower bound equals (1α 2 ) n-1 d-1 which is tight w.r.t. to the previous result in Sect. 6.2 (i.e. (1α 2 ) n d-n ) for n much smaller than d.

• η → 1: The lower bounds converges to 0 and is also tight since m ⋆ → x 1 as the vectors are converging to x 1 .

In the end, under the Gaussian assumption, the Kullback-Leibler distance between both distributions increases which proves that identifying the positive memory units becomes easier as η increases (see Fig. 19).

CONCLUSION

In this paper, we take a statistical signal processing point of view for image indexing, instead of traditional geometrical approaches in the literature. This shift of paradigm allows us to bring theoretical justifications for representing a set of vectors. We have presented and analyzed two strategies for designing memory vectors, enabling efficient membership tests for real-valued vectors. We have also showed two possible assignment strategies and analyzed their performance theoretically and experimentally. For random assignment, the optimized pinv construction gives better results than the simple sum aggregation. On the other hand, when the vectors in the same memory unit share some correlation, sum is on par with the pinv construction as for the quality of the hypothesis test. Yet, the pinv construction when used in the weakly supervised assignment offers a lower imbalance factor. This yields less variability of the search runtime from one query to another. This procedure is done offline and its complexity is often ignored in the image search literature. On the contrary, we did pay attention to this bottleneck: we proposed to run the weakly supervised assignment by batch and showed that it does not spoil the overall performance of the image search.

Fig. 1 .

 1 Fig.1. Ratio of the global cost by the cost of the exhaustive search C H 0 /N as a function of n using random assignment on synthetic data. Setup: d = 1000, ǫ = 10 -2 . The different curves correspond to values of α 0 ∈ {0.5, 0.6, 0.7, 0.8, 0.9}. Red and blue lines correspond to sum and pinv constructions respectively. Darker shades correspond to higher α 0 .

Fig. 2 .Fig. 3 .

 23 Fig. 2. ROC curves (1 -P fn as a function of P fp ) for sum and pinv constructions evaluated using (left) synthetic data with d = 100 and n = 10, (right) Holidays+Flickr1M real dataset and φ ∆ features with d = 1920 and n = 10.

Fig. 4

 4 Fig.[START_REF] Charikar | Similarity estimation techniques from rounding algorithms[END_REF]. Probability that a memory unit contains at least one match with respect to their rank. Green, red and blue lines correspond to random, sum spherical k-means and pinv spherical k-means respectively. There is a high probability of retrieving a match in highly ranked memory units, but it decreases faster with a weakly supervised assignment.

  Fig.5. The average number of matches given that the memory unit is positive. Green, red and blue lines correspond to random, sum spherical k-means and pinv spherical k-means respectively. Using random assignment, we have about 1 matching vector per memory unit. The weakly supervised assignment improves this especially for higher-ranked units.

Fig. 6 .Fig. 7 .

 67 Fig.6. Complexity per query for pinv k-means (blue) and sum k-means (red). Although the averaged complexity ratios over all queries are similar, pinv has less variance thanks to lower imbalance factor.

Fig. 9 .

 9 Fig. 9. The performance of memory vectors after PCA dimensionality reduction with d ′ = 1024.

Fig. 10 .

 10 Fig.10. Image retrieval performance with deep learning features (d = 4096), as trained by Babenko et al.[START_REF] Babenko | Neural codes for image retrieval[END_REF]. Features for Holidays and Oxford5K are retrained on the Landmarks dataset, whereas the ones for UKB are trained on ILSVRC.

Fig. 12 .

 12 Fig. 12. Experiments with binary codes after PCA reduction to d ′ = 1024. The quantization is asymmetric: real query vectors are compared to binary memory vectors.

Fig. 13 .Fig. 14 .Fig. 15 .

 131415 Fig.[START_REF] Sivic | Video Google: A text retrieval approach to object matching in videos[END_REF]. Image retrieval performance in Holidays+Flickr1M with different memory unit size n. The data is streamed in this scenario, and the memory untis are created randomly.

• η = - 1 :

 1 The lower bound equals n d which is tight w.r.t. to the previous result in Sect. 6.2 (i.e. n d (1 + n (d-n) )) for n much smaller than d.

  1, and large-scale and streaming data in Section 5.2.

		sum sum + norm. pinv pinv + norm.
	Holidays	2.08	2.17 1.90	2.03
	Oxford5k 2.76	2.69 2.27	2.23
	UKB	2.58	2.56 2.06	2.09

TABLE 1

 1 Imbalance factor for different datasets using sum, pinv and their normalized variants of k-means. Each dataset is clustered into M = N/10.

  Fig. 11. Experiments with binary codes after PCA reduction to d ′ = 1024. The quantization is symmetric: real query vectors are binarized and then compared to binary memory vectors.

		0.7		Holidays
		0.6
	mAP	0.1 0.2 0.3 0.4 0.5 0 0.6	0.2	0.4 Oxford5k 0.6	baseline bin. sum bin. pinv 0.8	1
		0.5
	mAP Score	0.1 0.2 0.3 0.4 0 0 0.5 1 1.5 2 2.5 4 3.5 3	0.2	0.4	UKB	0.6	baseline bin. sum bin. pinv 0.8 baseline bin. sum bin. pinv	1
			0.2	0.4		0.6	0.8	1
				Complexity Ratio
	mAP	0.1 0.2 0.3 0.4 0.8 0.7 0.6 0.5 0 0.6	0.2	Holidays 0.4 0.6 Oxford5k	baseline bin. sum bin. pinv 0.8	1
		0.5
	mAP Score	0.1 0.2 0.3 0.4 0 0 0.5 1 1.5 2 2.5 4 3.5 3	0.2	0.4	UKB	0.6	baseline bin. sum bin. pinv 0.8 baseline bin. sum bin. pinv	1
			0.2	0.4		0.6	0.8	1
				Complexity Ratio

  Complexity of each query for Holidays+Flickr1M using our batch spherical k-means (blue) and mbk (red). Our scheme has less complexity on average and less variance. Fig. 17. Comparison of pinv and sum based methods in Yahoo100M. Our pinv -based k-means variant achieves baseline performance with only 0.12 complexity ratio.

		1
		0.8
	mAP	0.6 0.4				kmeans-pinv random-pinv
	) Fig. 16. 0.2 0.15 0.20 0.25 0.30 0.35 0.40 Complexity Holidays_Flickr1M	0.1	0.12	0.14	0.16 Complexity Ratio 0.18 kmeans-sum 0.2	0.22	0.24

  Fig. 18. Expectations and Variances under both hypothesis and Kullback Leibler distance as functions of η for the sum construction. α = 0.5, d = 512.

	6					Expectation
			H 0
	4		H 1
	2
	0	1	-0.8	-0.6	-0.4	-0.2	0	0.2	0.4	0.6	0.8
	0.2					Variance
	0.15
	0.1
	0.05
	-1 0	-0.8	-0.6	-0.4	-0.2	0	0.2	0.4	0.6	0.8
	600					Kullback Leibler Distance
	400
	200
	-1 0	-0.8	-0.6	-0.4	-0.2	0	0.2	0.4	0.6	0.8
						η

  Fig. 19. Expectations and Variances under both hypothesis and Kullback Leibler distance as functions of η for the pinv construction. α = 0.5, d = 512.

	0.8				Expectation
	0.6	H 0 H 1
	0.4
	0.2
	-1 0	-0.8	-0.6	-0.4	-0.2	0	0.2	0.4	0.6	0.8
	0.02				Variance
	0.015
	0.01
	0.005
	-1 0	-0.8	-0.6	-0.4	-0.2	0	0.2	0.4	0.6	0.8
	600				Kullback Leibler Distance
	400
	200
	-1 0	-0.8	-0.6	-0.4	-0.2	0	0.2	0.4	0.6	0.8
					η

Note that, when we use such assignment techniques in spherical kmeans, the number of vectors per clusters is not evenly distributed. The dot product may be dominated by long cluster representation vectors. Hence, we also propose a normalized version of the assignment, where the cluster representation vectors obtained are normalized to the unit norm.

ACKNOWLEDGMENTS

This work was supported, in part, by the Natural Sciences and Engineering Research Council of Canada through grants RGPAS 429296-12 and RGPIN/341596-12, by the CHIST-ERA project ID IOT (ANR-16-CHR2-0005), and by the European Research Council under the European Union's Seventh Framework Programme ( FP7 / 2007 -2013 ) / ERC grant agreement 290901.