Paolo Piro
email: piro@i3s.unice.fr

Sandrine Anthoine
email: anthoine@i3s.unice.fr

Eric Debreuve
email: debreuve@i3s.unice.fr

Michel Barlaud
email: barlaud@i3s.unice.fr

IMAGE RETRIEVAL VIA KULLBACK-LEIBLER DIVERGENCE OF PATCHES OF MULTISCALE COEFFICIENTS IN THE KNN FRAMEWORK

Keywords: Image retrieval, sparse wavelet description, intra/interscale dependency, Kullback-Leibler divergence, k-th nearest neighbors

In this paper, we define a similarity measure between images in the context of (indexing and) retrieval. We use the Kullback-Leibler (KL) divergence to compare sparse multiscale image representations. The KL divergence between parameterized marginal distributions of wavelet coefficients has already been used as a similarity measure between images.

Here we use the Laplacian pyramid and consider the dependencies between coefficients by means of non parametric distributions of mixed intra/interscale and interchannel patches. To cope with the high-dimensionality of the resulting description space, we estimate the KL divergences in the k-th nearest neighbor (kNN) framework (instead of classical fixed size kernel methods). Query-by-example experiments show the accuracy and robustness of the method.

INTRODUCTION

A central question in content-based image indexing is to define a similarity measure between images that matches -or at least is close enough to -our perception of their similarity. Then, database images can be simply ranked in increasing order of their similarity to the reference (or example) image for a query-by-example task. Understanding how human perceive the similarity between images via perceptual studies is still a topic of active research. Thus, content-based image indexing systems relying on such studies may be subjective and very hard to implement. Here, we focus on developing an objective and mathematically defined measure that will be easily implementable.

In this context, used measures often rely on global descriptions such as dominant colors or color distribution, or on extracted information such as salient points/regions together with local features or segmentation along with region arrangement [START_REF] Mezaris | Object-based mpeg-2 video indexing and retrieval in a collaborative environment[END_REF][START_REF] Zhang | Optimizing metrics combining low-level visual descriptors for image annotation and retrieval[END_REF]. The philosophy here is to use a sparse multiscale image description. The Kullback-Leibler (KL) divergence has This work is supported in part by the ANR "ICOS-HD". already been used as a similarity measure between parameterized marginal distributions of wavelet coefficients at different scales [START_REF] Do | Waveletbased texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF][START_REF] Wang | Quality-aware images[END_REF]. Nevertheless, independence between the coefficients was assumed, preventing from taking into account local image structures such as texture. In contrast, we propose to consider dependency by means of distributions of mixed intra/interscale patches of the Laplacian pyramid coefficients. In addition, for the case of color images, we take into account the statistical dependencies amongst the three color channels; hence patches of coefficients are also interchannel. This approach implies to deal with a high-dimensional statistical description space. The number of samples being too small to reasonably fill this space, fixed size kernel options to estimate distributions or divergences fail. Alternatively, we propose to estimate the KL divergence in the k-th nearest neighbor (kNN) framework [START_REF] Boltz | High-dimensional kullback-leibler distance for region-of-interest tracking: Application to combining a soft geometric constraint with radiometry[END_REF], i.e., adapting to the local sample density and directly from the samples.

SIMILARITY BETWEEN IMAGES

In this section we define our similarity measure between images. It is a combination of Kullback-Leibler divergences in a multiscale feature space. This feature space consists of inter/intrascale and interchannel patches of Laplacian pyramid coefficients for color images. We first describe the feature space and then the use and estimation of the KL divergences.

Patches of Laplacian pyramid coefficients

Let us denote by w(I) j,k the coefficient for image I at scale j and location in space k for a general multiresolution decomposition.

The concept of patches of multiresolution coefficients was introduced by [START_REF] Portilla | Image denoising using a scale mixture of Gaussians in the wavelet domain[END_REF] for the wavelet decomposition under the name "neighborhoods of wavelet coefficients". It stems from the following ideas.

The wavelet domain provides a sparse representation of images, meaning that it concentrates the information into a few coefficients of large amplitude (the rest of the coefficients being small) and enjoys a fast transform. This is what makes wavelet thresholding methods such powerful tools in image processing. Initial thresholding wavelet methods treat each Preprint Accepted in CBMI, 18-20/06/2008, London, UK -Copyright IEEE coefficient separately relying on the fact that these coefficients are decorrelated. However, they are not independent and these dependencies are the signature of structures present in the image. For example, a discontinuity between smooth regions at point k 0 will give large coefficients at this point at all scales j (w(I) j,k0 large for all j). The most significant dependencies are seen between a wavelet coefficient w(I) j,k and its closest neighbors in scale: w(I) j-1,k , or in space: w(I) j,k±(0,1) , w(I) j,k± [START_REF] Mezaris | Object-based mpeg-2 video indexing and retrieval in a collaborative environment[END_REF]0) . Several models using these dependencies have been proposed and used in image enhancement [START_REF] Portilla | Image denoising using a scale mixture of Gaussians in the wavelet domain[END_REF][START_REF] Romberg | Bayesian tree-structured image modeling using wavelet-domain hidden markov models[END_REF]. The concept of patches (or neighborhoods) of wavelet coefficients, was introduced in [START_REF] Portilla | Image denoising using a scale mixture of Gaussians in the wavelet domain[END_REF]. These are vectors of the form: w(I) j,k = w(I) j,k , w(I) j,k±(1,0) , w(I) j,k±(0,1) , w(I) j-1,k

(1) for a greyscale image.

The probability density function (pdf) of such patches was shown to characterize and estimate fine spatial structures in greyscale images [START_REF] Portilla | Image denoising using a scale mixture of Gaussians in the wavelet domain[END_REF][START_REF] Pierpaoli | Reconstructing sunyaev-zeldovich clusters in future cmb experiments[END_REF]. Hence such patches are expected to be relevant features to represent the image content.

Critically sampled tensor wavelet transforms lack of rotation and translation invariance and so would the neighborhoods made of such coefficients. We prefer to use the Laplacian pyramid [START_REF] Burt | The laplacian pyramid as a compact image code[END_REF] which shares with the wavelet transform the sparsity and inter/intrascale dependencies properties and is more robust to rotations. Hence we will define our feature space from the set of patches of Laplacian pyramid coefficients of the form of [START_REF] Mezaris | Object-based mpeg-2 video indexing and retrieval in a collaborative environment[END_REF].

Here, we consider colored images in the luminance -chrominances space: I = (I Y , I U , I V ). Since the pyramid coefficients are correlated through channels, we aggregate in the patches the coefficients of the three channels (still noting the patches w):

w(I) j,k = w(I Y ) j,k , w(I U ) j,k , w(I V ) j,k (2)
with

w(I c ) j,k = w(I c ) j,k , w(I c ) j,k±(1,0) , w(I c ) j,k±(0,1) , w(I c ) j-1,k . (3)
Hence our feature space is the set of the patches as in Eq. ( 2) for all scales j and locations k. Fig. 1 illustrates how to build one patch from two near subbands of the Laplacian pyramid in the case of a single color component image (as in Eq. ( 1)). The low-frequency subband is also considered in a similar way. Namely, intrascale and interchannel patches (of dimension 27) are built by joining the spatial 3 × 3 neighborhoods in the three channels. As a result, our feature space is the set of Laplacian pyramid coefficients defined in Eq. ( 2) (of dimension 18), for all scales j and locations k in addition to the low-frequency patches.

Let us now turn to the measure of similarity on this space.

Similarity measure between images

Since geometrically modified or slightly degraded versions of the same image as well as images containing similar objects should be close, one cannot define a measure comparing directly the neighborhoods one by one, but rather their probability distributions. More specifically, we consider the pdf of the neighborhoods of Eq. ( 2) for each scale, i.e. we consider the pdf p wj (I) of the set neighborhoods {w(I) j,k } k for each fixed j.

The considered pdf are those of coefficients that carry the informational content on the signal when they are large. The natural way to compare such pdf is to use measures derived from information theory. Here, we use the KL divergence between pdf, an approach that has also been successfully taken for other applications [START_REF] Boltz | High-dimensional kullback-leibler distance for region-of-interest tracking: Application to combining a soft geometric constraint with radiometry[END_REF][START_REF] Angelino | A nonparametric minimum entropy image deblurring algorithm[END_REF]. This was also done in [START_REF] Do | Waveletbased texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF][START_REF] Wang | Quality-aware images[END_REF] to evaluate the similarity between images using the marginal pdf of the wavelet coefficients. We propose to use this measure on the multidimensional pdf of the neighborhoods of coefficients: the similarity between images I 1 and I 2 is a weighted sum over scales of the KL divergences between the Preprint Accepted in CBMI, 18-20/06/2008, London, UK -Copyright IEEE pdf p wj (I1) and p wj (I2) :

S(I 1 , I 2 ) = j α j D kl (p wj (I1) ||p wj (I2) ) (4)
where p wj (Ii) is the pdf of the pyramid patches (Eq.( 2)) of image I i at scale j and α j > 0 are weights (chosen according to the redundancy of the transform in each scale). Previous works on neighborhoods of wavelet coefficients or indexation using marginal pdf of the wavelet coefficients all assumed a parametric model for the pdf involved. In the marginal case, efficient models (e.g. generalized Gaussian [START_REF] Do | Waveletbased texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF][START_REF] Wang | Quality-aware images[END_REF]) lead to an analytic expression of the KL divergence as a function of the model parameters; but they are not easily generalizable to the multidimensional correlated case of multiscale patches. On the other hand, efficients multidimensional models including correlations (e.g. Gaussian mixtures [START_REF] Pierpaoli | Reconstructing sunyaev-zeldovich clusters in future cmb experiments[END_REF]) fit a wide variety of multidimensional pdf but impose to estimate the KL divergence after estimating model parameters. Besides the heavy computational cost of the consecutive estimations, the stability of such cascading estimates is likely to be difficult to obtain numerically. We prefer to make no hypothesis on the pdf at hand, hence sparing the cost of fitting the model parameters but with the need of estimating the KL divergences in this non-parametric case.

Estimation of the Kullback-Leibler (KL) divergence

Let us first remind that the expression of the KL divergence between two continuous pdf p 1 and p 2 is:

D kl (p 1 ||p 2 ) = p 1 (x) log p 1 (x) p 2 (x) dx = H x (p 1 , p 2 ) -H(p 1 )
(5) with H the differential entropy and H x the cross entropy.

The estimation of statistical measures in the multidimensional case is hard. In particular, kernel-based methods such as Parzen estimates become unadapted due to the sparsity of samples in high dimension (curse of dimensionality): the tradeoff between a kernel with a large bandwidth to perform well in low local sample density (which oversmoothes the estimator) and a kernel with a smaller bandwidth to preserve local statistical variability (which results in an unstable estimator) cannot always be achieved. We use instead the k-th nearest neighbor (kNN) framework [START_REF] Scott | Multivariate Density Estimation: Theory, Practice, and Visualization[END_REF] to compute the KL divergence. Indeed it follows the dual approach to the above fixed size kernel: the bandwidth adapts to the local sample density by letting the kernel contain exactly k neighbors of a given sample. Moreover it allows direct estimation of the divergence without explicitly estimating the pdf.

Fix ν, a set of N ν samples w 1 , w 2 ,..,w Nν of pdf p ν and k, a non-negative integer. kNN balloon estimates are based on the principle that p ν (s) is inversely proportional to the volume of the sphere containing the k nearest neighbors of s in ν [START_REF] Scott | Multivariate Density Estimation: Theory, Practice, and Visualization[END_REF]:

p ν (s) ∼ k v d ρ d k,ν (s) (6)
with v d the volume of the unit sphere in R d and ρ k,ν (s) the distance of s to its k-th nearest neighbor in ν -{s}.

Plugging this density estimator in Ahmad-Lin [START_REF] Ahmad | A nonparametric estimation of the entropy absolutely continuous distributions[END_REF] entropy estimator:

H al (p ν ) = - 1 N ν Nν n=1 log p ν (w n ) (7)
we get a biased entropy estimate. This bias can be rid of by replacing log(k) by the digamma function ψ(k), giving the entropy estimate [START_REF] Goria | A new class of random vector entropy estimators and its applications in testing statistical hypotheses[END_REF]:

H(p ν ) = log (N ν -1)v d -ψ(k) + d N ν Nν n=1 log ρ k,ν (w n ) (8)
The cross entropy estimate is then [START_REF] Boltz | High-dimensional kullback-leibler distance for region-of-interest tracking: Application to combining a soft geometric constraint with radiometry[END_REF]:

H x (p ν1 , p ν2 ) = log N ν2 v d -ψ(k)+ d N ν2 Nν 2 n=1 log ρ k,ν1 (w 2 n ) (9)
and the KL divergence estimate is:

D kl (p ν1 ||p ν2 ) = log Nν 2 Nν 1 -1 + d N ν2 Nν 2 n=1 log ρ k,ν1 (w 2 n ) - d N ν1 Nν 1 n=1 log ρ k,ν1 (w 1 n ) (10)
This expression is valid in any dimension and it is robust to the choice of k.

NUMERICAL EXPERIMENTS

Settings

We used for our numerical experiments a database containing images from the Recognition Benchmark collection, which is a ground-truth dataset already used in [START_REF] Nistér | Scalable recognition with a vocabulary tree[END_REF]. The database consists in 640x480 images grouped by sets of four images of the same scene. Hence, for any query image, exactly three other relevant images are to be retrieved. A subset of 100 images from this dataset was used to adjust some parameters of our algorithm and compare its performance with a reference method (section 3.2 and 3.5). The method we propose was then validated on a larger dataset, containing 1,000 images from the same collection (section 3.3 and 3.4).

In all experiments a Laplacian pyramid was computed for each channel of the images (in the YUV color space) with a 5-points binomial filter. The first three high-frequency subbands and the low-frequency image approximation were used to build patches.

To perform the KL divergence estimations, we used a fast kNN search algorithm of complexity O(N log N ), N being the number of sample points [START_REF] Arya | An optimal algorithm for approximate nearest neighbor searching fixed dimensions[END_REF]. This algorithm was run on In order to speed up the computing time for queries on a large dataset, we used a parallel implementation of the kNN searching on Graphic Processing Unit (GPU) [START_REF] Garcia | Fast k Nearest Neighbor Search usign GPU[END_REF]. This is based on a brute-force approach and was written in the CUDA development environment. It was run on a NVIDIA GeForce 8800 GTX graphic card and allowed to compute a single similarity measure in less than 0.2 s on average.

The sparsity property of the multiresolution transform allowed us to retain not all feature vectors but only a small proportion of them in the high-frequency subbands, thus reducing the number of operations required for kNN estimations. For this purpose, only a certain proportion of patches was taken into account, according to a significance criterion. Namely, all subbands in the luminance channel were thresholded based on either the energy (quadratic norm) of the whole patch or the amplitude of the patch center, so that all coefficients being under the threshold do not contribute to entropy estimation. Then, each significant patch in the luminance channel and the corresponding ones in the two chrominance channels were joined together as described in section 2.1.

With respect to the kNN estimation of the KL divergence, we fixed k to ten and the weights α j of Eq.( 4) to one, so that all subbands contribute equally to the similarity measure.

Robustness to sparsity of multiscale features

The similarity measure we propose is robust with respect to the sparsity of the subband image decomposition. Indeed, retrieval results are consistent with ground-truth, as we decrease the number of subband patches contributing to distance esti-mations. In order to select the most significant patches, (i.e. those containing the most important information of the image), we selected the largest ones (either in the sense of overall energy or in the sense of center amplitude). Given a subset of the database containing both relevant and non relevant images for a given query, we computed the similarity measures using different proportions of significant patches. Results are generally consistent as sparsity increases (i.e. as the number of selected patches decreases), in the sense that relevant images remain significantly closer to the query than others and the rank order of retrieved images does not change. An example is given in Fig. 2, where similarity between one query and the first four retrieved images is shown for different proportions of selected subband coefficients. This trend is general for all pictures and confirms the suitability of the selection criterion: it allows to reduce the computing time for the distance estimation, while preserving its accuracy and sensitivity.

Retrieval results

Results for five queries in the database containing 1,000 images are displayed in Fig. 3. In this figure, each row displays the retrieval result for the query image shown on the leftmost column. From the second column on, one can see the first 3 retrieved images ranked in increasing order of their similarity measure from the query. Hence the second leftmost image is the most similar, excluding the reference image which is always ranked first with a distance of zero. The criterion we adopted to select feature vectors was based on patch energy. Proportions of selected patches in the three considered highfrequency subbands were respectively: 1/64, 1/32, 1/16.

As shown by experimental results, the method we propose seems to perform very well. In particular, images being relevant for a given query are generally ranked first; this holds in spite of the fact that images belonging to the same group have been often subjected to different geometric transformations. Hence our method is robust to such transformations, e.g. rotations (such as the third row in Fig. 3) or changes in viewpoint (rows 2 and 5) and zoom (rows 1 and 4).

Evaluation of image retrieval performances

A wide variety of measures have been proposed to evaluate image retrieval performance [START_REF] Smith | Image retrieval evaluation[END_REF]. We adopt a standard criterion in information retrieval, which has been also used in the context of image retrieval [START_REF] Mikolajczyk | A performance evaluation of local descriptors[END_REF]. It is based on the number R of expected results for a given query (relevant images), the number D of correct results (detected images) and the number W of wrong results (false positive images), with respect to the ground-truth relevance. Naturally, these measures depend on the number C of retrieved images (cut-off ). Together, they enable to define the standard measures of precision and recall in the following way: Varying the cut-off value, one obtains the ROC (Receiver Operating Characteristic) curve, which gives the performance of the retrieval system for each request. This curve is generally shown as recall versus 1 -precision, that represent, respectively, the detection rate and the false positive rate. Hence the larger are precision and recall values, the better is retrieval performance.

precision = D C = D D+W , recall = D R . (11)
To evaluate the overall performance of our method, we used the curves obtained by querying the whole database with each of the 1,000 images. The significance criterion for patch selection was that based on patch energy, that showed the best performance.

By averaging over all queries the individual values of precision and recall, we obtained one average curve. It is displayed in Fig. 4. This graph shows that the best trade-off between precision and recall was reached when we retrieved three images; in this case, cut-off value matches exactly the number of relevant images or, in other words, there is a high probability that retrieved images are all and only the relevant ones.

Comparison with SIFT-based retrieval

Retrieval experiments were done by using a state-of-the-art method as well; this is based on SIFT descriptors, that represent the gradient orientations at interest points. It uses a matching criterion between the interest points of images to measure their similarity. We used a Matlab implementation of the SIFT algorithm [START_REF] Lowe | Sift keypoint detector[END_REF]. Since it took about 4.8 s on average to compute the similarity between two images, experiments on the whole database could not be performed in a reasonable time. As a result, we compared the two methods by querying a subset containing 100 images. The results are shown in Fig. 5 in terms of the ROC curves. Namely, the KL divergence between patch distributions and the number of matched SIFT keypoints between images are compared as similarity measures for image retrieval.

The performances of our method are shown to be very close to that of the SIFT-based method, which is a reference method for the task of content-based image retrieval.

CONCLUSION

In this paper, we proposed a new image retrieval method based on multidimensional probability distributions of multiscale coefficients. These are grouped in coherent patches, that are selected by a significance criterion in order to effectively represent image features. The patches are build by taking into account intrascale, interscale and interchannel dependencies of subband coefficients for color images. The similarity measure we used is the sum over scales of the KL divergences between probability distributions of the image features, and it is estimated non-parametrically via a kNN approach. It was proven to be consistent with respect to the sparsity of the subband transform. Query-by-example experiments on real images confirm the suitability of the proposed measure in the context of image retrieval. It is in particular robust to different geometric transformations, such as change in viewpoint, rotation and zoom. Moreover, retrieval performances are comparable to those of a reference algorithm, based on local SIFT descriptors.

Further investigations are planned with larger databases and others benchmark image collections. Future works will also include a more extensive comparison with state-of-theart CBIR systems and the extension of the proposed method to the retrieval of video sequences.

Fig. 1 .

 1 Fig. 1. How to build a patch by grouping multiscale coefficients, for a single color channel.

Fig. 2 .

 2 Fig. 2. Evolution of the similarity measure with the proportion of patches selected. Retrieval results shown for image query 1 and its four most similar images (images 2, 3 and 4 are relevant for this query).

Fig. 3 .

 3 Fig. 3. Retrieval results for 5 images of the database. For each row, left to right: reference image; first 3 ranked images of the database (excluding reference image). For each retrieved image our distance to the query is also shown.

Fig. 4 .

 4 Fig. 4. Receiver Operating Characteristic (ROC) from cut-off variation.

Fig. 5 .

 5 Fig. 5. Performance of image retrieval for our method and SIFT-based algorithm in terms of the ROC curve.


Frank Nielsen

Paolo Piro

Michel Barlaud

Tailored Bregman Ball Trees for Effective Nearest Neighbors *

Nearest Neighbor (NN) search is a crucial tool that remains critical in many challenging applications of computational geometry (e.g., surface reconstruction, clustering) and computer vision (e.g., image and information retrieval, classification, data mining). We present an effective Bregman ball tree [5] (Bb-tree) construction algorithm that adapts locally its internal node degrees to the inner geometric characteristics of the data-sets. Since symmetric measures are usually preferred for applications in content-based information retrieval, we furthermore extend the Bb-tree to the case of symmetrized Bregman divergences. Exact and approximate NN search experiments using high-dimensional real-world data-sets illustrate that our method improves significantly over the state of the art [5], sometimes by an order of magnitude.

Introduction and prior work

Finding nearest neighbors (NNs) is a very common task occurring in many applications ranging from computer vision to machine learning and data mining. Let S = {p 1 , ..., p n } be a set of n d-dimensional points (with d typically ranging up to a few thousand dimensions). Given a query point q, the nearest neighbor NN(q) is defined to be the "closest" point of S with respect to a dissimilarity measure D: NN(q) = arg min i D(q, p i ). Instead of considering the closest neighbor, queries can be enlarged to report the first k "closest" points (useful for k-NN classification rule in machine learning). It is usually enough to get a good neighbor [START_REF] Chazelle | Technical perspective: finding a good neighbor, near and fast[END_REF] by relaxing the exact search to get near and fast neighbors. Besides the theoretical puzzling questions related to the dreaded curse of dimensionality [START_REF] Chazelle | Technical perspective: finding a good neighbor, near and fast[END_REF], practitioners make every endeavor to speed up applications by designing tailored schemes for improving over the naïve linear-time O(dn) brute-force method. Among the flourishing literature of nearest neighbor techniques, we distinguish two main sets of methods: (1) those relying on tree-like space partitions with branch-and-bound * Financially supported by ANR-07-BLAN-0328-01 GAIA (Computational Information Geometry and Applications) and DIGITEO GAS 2008-16D (Geometric Algorithms & Statistics)

† École Polytechnique, LIX, Palaiseau, France, nielsen@lix.polytechnique.fr ‡ CNRS I3S, Sophia-Antipolis, France, {piro,barlaud}@i3s.unice.fr queries, such as kD-trees, metric ball and vantage point trees [START_REF] Cayton | Fast nearest neighbor retrieval for Bregman divergences[END_REF][START_REF] Kumar | What is a good nearest neighbors algorithm for finding similar patches in images[END_REF], and (2) those based on mapping techniques [START_REF] Athitsos | Nearest neighbor retrieval using distancebased hashing[END_REF] (e.g., locality-sensitive hashing, random projections). The former tree-based methods improve over the brute force algorithm by pruning sub-trees whenever traversing the trees with queries (exact NN) or stopping after visiting a given "budget" of leaves (approximate NN). The latter methods concentrate on reducing dimensions while preserving as much as possible distances and controlling geometrically collisions of hash functions (only work for approximate NN). Since the Euclidean distance is often inappropriate for meaningfully measuring the proximity of feature points arising from, say, image applications, a wide range of distortion measures has been proposed and NN data-structures have been first extended to arbitrary metrics (e.g., vantage point trees [START_REF] Kumar | What is a good nearest neighbors algorithm for finding similar patches in images[END_REF]). However, these NN search methods rely fundamentally on the triangle inequality property, which is not satisfied by informationtheoretic statistical distances such as the Kullback-Leibler divergence. Cayton [START_REF] Cayton | Fast nearest neighbor retrieval for Bregman divergences[END_REF] recently proposed the analogous of metric ball trees for the broad class of Bregman divergences. Bregman divergences D F on vector set X ⊂ R d are defined for a strictly convex and differentiable generator F (x) :

X ⊂ R d → R as D F (p||q) = F (p) -F (q) -(p -q) T ∇F (q)
, where ∇F denotes the gradient of F . Those divergences (parameterized by a generator F ) include all quadratic distances (also known as Mahalanobis squared distances, which are the only symmetric Bregman divergences and are obtained for F (x) = Σ -1 x, where Σ ≻ 0 is the positive-definite variance-covariance matrix) and the asymmetric KL divergence (F (x) = d j=1 x j log x j , the negative Shannon entropy), for which D F (p||q) = D F (q||p). Many fundamental algorithms (e.g., kmeans [START_REF] Banerjee | Clustering with Bregman divergences[END_REF], PCA) and data-structures (e.g., Voronoi diagrams [START_REF] Nielsen | On Bregman Voronoi diagrams[END_REF]) have been generalized to that class of divergences, thus offering meta-algorithms that can work for any Bregman divergence. For example, Banerjee et al. [START_REF] Banerjee | Clustering with Bregman divergences[END_REF] showed that the celebrated Lloyd k-means algorithm extends to (and only to) the class of Bregman divergences unifying various algorithms. Since Bregman divergences are typically non-metric asymmetric distortion measures, we consider both left/right-sided and symmetrized NN searches defined as follows: NN r F (q) = arg min i D F (p i ||q) (rightsided), NN l F (q) = arg min i D F (q||p i ) (left-sided) and NN F (q) = arg min i (D F (p i ||q) + D F (q||p i ))/2 (sym-metrized).

Bregman ball trees (Bb-trees)

Without loss of generality, we consider only rightsided NN queries, as left-sided NN queries can be handled similarly by considering the dual divergence D F * (∇F (q)||∇F (p)) = D F (p||q) arising from the Legendre-Fenchel conjugate F * of F (with ∇F * the functional inverse of ∇F -1 ). See [START_REF] Nielsen | On Bregman Voronoi diagrams[END_REF][START_REF] Cayton | Fast nearest neighbor retrieval for Bregman divergences[END_REF] for details. Bregman generators come by pairs (F, F * ), e.g. the dual Legendre conjugates F (x) = x log x and F * (y) = exp y (with F ′ (x) = (F * ′ (y)) -1 ), see [START_REF] Nielsen | On Bregman Voronoi diagrams[END_REF].

Outline of Bregman ball trees (Bb-trees) [5]

Similar to metric ball trees, a Bb-tree is built in a topdown fashion by applying recursively a partitioning scheme. First, the root is created to handle the source data-set S. A 2-means (k-means [START_REF] Banerjee | Clustering with Bregman divergences[END_REF] with k = 2) is computed, splitting S according to the two centroid points, say c l and c r , with respect to D F . These two centroids defines a partition S = S l ∪ S r that can be covered geometrically with corresponding Bregman balls B(c l , R l ) and B(c r , R r ) (possibly overlapping). This hierarchical decomposition of S is carried out recursively on S l and S r until a stopping criterion is eventually met. Two such typical termination criteria are: (1) a predefined maximum number of points l 0 (stored at leaves), or (2) a prescribed maximum radius (r 0 , comes in handy for approximate NN). Note that the source points are stored only at leaves of the Bbtree, and the sub-sets S l and S r may be theoretically unbalanced. Internal nodes store only two left/right Bregman balls covering the point sets stored at their left/right sub-trees.

In both exact and approximate NN queries, we perform a branch-and-bound search to answer queries. For a given query q, the tree is traversed in depthfirst-search order from the root to the leaves. At an internal node, we choose to branch first on the sub-tree whose corresponding ball is "closer" to the query q (the sibling is temporarily ignored). Once a leaf node is reached, the closest point to q among the points stored at the leaf is computed using the brute-force method. This first visited leaf yields the very first NN point candidate p ′ , thus giving an upper bound D F (p ′ ||q) to the NN distance. In exact search, the tree traversal goes on through all formerly ignored subtrees. To decide whether a subtree must be explored or not, we check whether D F (p ′ ||q) > min x∈B(c,R) D F (x||q), where p ′ is the current NN candidate. This test is performed by projecting q onto the Bregman ball B(c, R). The projection of a point q onto a Bregman ball B is the unique minimizer q B = arg min x∈B D F (x||q) [START_REF] Nielsen | On Bregman Voronoi diagrams[END_REF]. Instead of projecting exactly the query point onto the ball, we rather make

use of a Bregman annulus B(c, R, R ′ ) = {x | R ≤ D F (x||c) ≤ R ′ } that encloses the exact projection by construction: q B ∈ B(c, R, R ′ ). If D F (p ′ ||q) < R then the node can be pruned; If D F (p ′ ||q) > R ′
then the node must be explored. These lower/upper bounds are computed during the geodesic bisection search of the projection [START_REF] Cayton | Fast nearest neighbor retrieval for Bregman divergences[END_REF] (see Section 2.4). Overall, observe that the less the number of pairwise intersecting balls stored at nodes, the better the Bb-tree performance. Although exact NN retrieval on Bb-trees can often be achieved with much smaller computational cost than the brute-force search, the practical interest of Bb-trees is to get significant speed-up search when performing approximate NN queries. The approximate search allows one for large speed-up as it does not perform exhaustive branch-and-bound search. A common criterion to perform approximate search is to stop the branch-and-bound algorithm after exploring a prescribed number of leaves [START_REF] Cayton | Fast nearest neighbor retrieval for Bregman divergences[END_REF].

Speeding up construction time: Bb-tree++

In order to preprocess datasets efficiently by means of Bregman ball trees, instead of running the full regular Bregman k-means algorithm [START_REF] Banerjee | Clustering with Bregman divergences[END_REF], we just perform a careful light initialization of the two cluster centers (seeds). Initialization turns out to be the crucial stage of k-means. That is, k-means locally optimizes the potential P (C) = p∈S min c∈C D F (p||c), where C denotes the set of k centers. Each round assignment/cluster adjustment decreases this potential function so that monotonous convergence is guaranteed. It is striking to know that the worst-case running-time of k-means is theoretically exponential with the dimension although a recent smooth analysis yields polynomial amortized time [START_REF] Manthey | Improved Smoothed Analysis of the k-Means Method[END_REF]. A bad initialization may trap k-means into a local optimum. Initialization is thus all the more important, as k-means is called at each internal node of the Bb-tree to partition the data-sets into two sub-sets. Therefore, the idea is to replace the k-means local iterative algorithm by a well-chosen initialization that provides a guaranteed upper-bound on the optimal partition. This quantum leap for k-means was discovered by Arthur and Vassilvitskii [START_REF] Arthur | k-means++: the advantages of careful seeding[END_REF] and later extended to Bregman divergences by Nock et al. [START_REF] Nock | Mixed Bregman clustering with approximation guarantees[END_REF]. We describe the initialization procedure of Bregman 2-means++ for a set S. First, draw uniformly at random the seed c l (cluster "center") among the points of the data-set. Then compute the Bregman divergence of this point to all points of S, and draw the second seed c r according to the divergence distribution: π i = DF (pi||c l ) p j ∈S DF (pj ||c l ) . Thus c r can never be c l , since D F (c l ||c l ) = 0 (the probability of drawing c l is zero). This careful initialization yields an extremely fast tree construction, which statistically provides nice splitting and tends to balance sub-sets stored at leaves. Arthur and Vassilvitskii [START_REF] Arthur | k-means++: the advantages of careful seeding[END_REF] proved that this initialization yields a k-means score P (C) at most 8(2+log k) of the optimal value. Similar bounds in O(log k) for generic Bregman divergences were later reported by Nock et al. [START_REF] Nock | Mixed Bregman clustering with approximation guarantees[END_REF]. Note that the two left/right Bregman balls stored at internal nodes tend to minimize the Bregman information [START_REF] Banerjee | Clustering with Bregman divergences[END_REF] (i.e., variance for the square potential F (x) = x 2 , mutual information for the negative Shannon entropy F (x) = x log x, etc). Next, we improve the partition at each internal node by learning its degree (number of siblings) from the local data-sets.

Learning the tree branching factor

Answering queries may range from optimal logarithmic time (i.e., the shortest path to a leaf) up to linear time for a complete tree traversal. Therefore, it is important to partition the source data-sets into as many as possible non-overlapping Bregman balls. However, consider a data-set consisting of three separated Gaussian point samples. Forcing this set to be split into two will likely create overlapping balls. Thus we better learn the number of sub-sets when partitioning the data so that the induced Bregman balls better fit the intrinsic geometric characteristics of the set. 1 We adapt the branching factor bf i (up to a maximum branching factor BF ) of each internal node of the Bbtree to the underlying distribution of the points by using the G-means strategy [START_REF] Hamerly | Learning the k in k-means[END_REF]. The underlying idea is to assume Gaussian distribution of each group of points (hence the name G-means). 2 Our use Gmeans [START_REF] Hamerly | Learning the k in k-means[END_REF] algorithm starts by setting k = 2 and then test for Gaussian distribution of the points using the Anderson-Darling statistical test. Given a confidence level α, if the Anderson-Darling normality test returns true, we keep the center, otherwise we split it into two. Between two rounds, we simply run Bregman k-means++ initialization on the data-set and get all the new centers that hopefully refine the partitioning. We enforce a maximum degree to each internal node in order to strike a balance between the average tree depth and the overall ball tree shape. The initialization to k clusters follows the same principle: Namely, we draw the l-th seed from the data-set uniformly according to the distribution π i = DF (pi||C l-1 ) p j ∈S DF (pj ||C l-1 ) , where C l-1 denotes the formerly chosen (l -1)-th seed and D F (p||S) = min x∈S D F (p||x).

When visiting the Bb-tree for answering nearest neighbor queries, we use a priority queue. To relax the exact NN to approximate NN queries, the criterion of stopping the search once a few leaves have been visited was proposed and successfully demonstrated by Cayton [START_REF] Cayton | Fast nearest neighbor retrieval for Bregman divergences[END_REF]. However, there was no guarantee to 1 Refer to Figure 1 of [START_REF] Hamerly | Learning the k in k-means[END_REF] for examples. 2 This is not restrictive as any smooth density function may be arbitrarily well approximated using a Gaussian mixture model.

get a good approximation to the exact NN because leaves containing sub-sets that are close to the exact NN are not necessarily close in the tree. We improve this point by a careful non-recursive implementation of the tree traversal, which allows us to order the nodes to be explored by their divergence to the query point (i.e., the divergence of the query to the ball centers). Using a priority queue guarantees that nearest nodes are always explored first when traversing back the tree, by jumping appropriately to the most likely not yet visited sub-tree.

Bregman projection onto balls

The geodesic Γ pq linking p to q is defined as Γ pq = {LERP(λ, p, q) | λ ∈ R}, with LERP(λ, p, q) = ∇F -1 ((1 -λ)∇F (p) + λ∇F (q)); see [START_REF] Nielsen | On Bregman Voronoi diagrams[END_REF][START_REF] Cayton | Fast nearest neighbor retrieval for Bregman divergences[END_REF][START_REF] Nielsen | Bregman sided and symmetrized centroids[END_REF]. To find the Bregman projection q B = arg min x∈B(c,R) D F (x||c) of a query point q onto a Bregman ball B(c, R), we first check whether q is outside the ball or not: D F (q||c) > R. If not, D F (q||c) ≤ R and the projection is simply the point itself: q B = q. Now consider the geodesic segment Γ cq , the projection point q B belongs necessarily to the geodesic Γ cq : q B = LERP(λ q , c, q) for some λ q ∈ [0, 1]. The value λ q is approximated by a bisection search as follows. Initially, let λ

M ] = [λ (0) , λ (1)
M ] and so on until the range size λ

(i) M -λ (i)
m goes below a threshold ǫ (typically ǫ ∈ [10 -10 , 10 -5 ]). This yields a fine approximation of q B ∼ LERP(λ (i) , c, q) by splitting the "λ" intervals a dozen of times. (Recall that the algorithm keeps a Bregman annulus, and refine the inner/outer radii to decide whether to prune or explore a sub-tree.)

Handling symmetrized Bregman divergences

Many content-based information retrieval systems (CBIRs) need symmetrical distortions measures.

Except for generalized quadratic distances (Mahalanobis), the symmetrized Bregman divergence is technically not a Bregman divergence [START_REF] Nielsen | On Bregman Voronoi diagrams[END_REF].

The symmetric KL divergence (SKL) goes by the name of Jensen-Shannon JS(p; q) = 1 2 KL(p|| p+q 2 ) + 1 2 KL(q|| p+q 2 ). It turns out that these symmetrized Bregman divergences [START_REF] Nielsen | Bregman sided and symmetrized centroids[END_REF] are all generalizations of the Jensen remainder obtained for convex generators F : JS F (p; q) = 1 2 (F (p) + F (q)) -F ( p+q 2 ). (Besides the Jensen-Shannon divergence for F (x) = x log x, these so-called Burbea-Rao divergences also includes the important COSH distance used in speech/sound processing for the Burg entropy F (x) = -log x.) The symmetrized Bregman centroid can be explicitly computed following the geodesic-walk algorithm of [START_REF] Nielsen | Bregman sided and symmetrized centroids[END_REF].

Experiments

Table 1 compares different partitioning methods when building a Bb-tree on the SIFT dataset [START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF] (10, 000 visual feature points extracted from images, encoded in dimension 1111). Bregman 2-means++ with iter = 0 means that only the initialization step is carried out. The speed-up is by comparison with the most computationally expensive k-means method. We used the maximum number of explored leaves as a parameter for stopping the branch-and-bound search. In each experiment of approximate search, we fixed a value of this parameter (from near-exact search to visiting only a single leaf), then we evaluated error rate and computational cost. Namely, the error rate is expressed as the number of closer points to the approximated NN (a quantity called Number Closer), while the speed-up is given by the ratio between the number of divergence computations in Bbtree over the brute-force search. Figure 1 displays a comparison of Cayton's Bb-tree with our Bb-tree++ in approximate search on SIFT dataset (log-log plot).

Figure 2 presents the results of approximate search queries on SIFT dataset wrt. SKL divergence (symmetrized Bregman). The marked point shows a speedup of 60.3 (10 1.78 ) when the Number Closer equals 14 (10 1.15 ), i.e. 0.14% of database points (log-log plot).

2 .

 2 midpoint m = LERP(λ (0) , c, q) obtained for λ (0) = λ If D F (m||c) > R then recurse on [λ

Figure 1 :

 1 Figure 1: BB-tree versus BB-tree++ (log-log plot scale). The axis denotes the error rate expressed as the number of closer points to the approximated NN.

Figure 2 :

 2 Figure 2: NN queries with respect to Jensen-Shannon divergence (symmetric Kullback-Leibler, SKL).

Table 1 :

 1 Construction time of Bregman ball trees.

		Bb-tree construction (bs = 50)
	method	iter depth depthavg Leav. speed-up
	2-m.	10	53	28.57	594	1
	2-m.++	10	58.33	31.18	647	1.03
	2-m.++	0	20	10.76	362	19.71

Acknowledgments

We thank L. Cayton for sharing with us his code and SIFT high-dimensional datasets.