Hervé Jégou

Teddy Furon

Jean-Jacques Fuchs

ANTI-SPARSE CODING FOR APPROXIMATE NEAREST NEIGHBOR SEARCH

Keywords: sparse coding, spread representations, approximate neighbors search, Hamming embedding

This paper proposes a binarization scheme for vectors of high dimension based on the recent concept of anti-sparse coding, and shows its excellent performance for approximate nearest neighbor search. Unlike other binarization schemes, this framework allows, up to a scaling factor, the explicit reconstruction from the binary representation of the original vector. The paper also shows that random projections which are used in Locality Sensitive Hashing algorithms, are significantly outperformed by regular frames for both synthetic and real data if the number of bits exceeds the vector dimensionality, i.e., when high precision is required.

INTRODUCTION

This paper addresses the problem of approximate nearest neighbor (ANN) search in high dimensional spaces. Given a query vector, the objective is to find, in a collection of vectors, those which are the closest to the query with respect to a given distance function. We focus on the Euclidean distance in this paper. This problem has a very high practical interest, since matching the descriptors representing the media is the most consuming operation of most state-of-the-art audio [START_REF] Casey | Content-based music information retrieval: Current directions and future challenges[END_REF], image [START_REF] Jégou | Improving bag-of-features for large scale image search[END_REF] and video [START_REF] Law-To | Video copy detection: a comparative study[END_REF] indexing techniques. There is a large body of literature on techniques whose aim is the optimization of the trade-off between retrieval time and complexity.

We are interested by the techniques that regard the memory usage of the index as a major criterion. This is compulsory when considering large datasets including dozen millions to billions of vectors [START_REF] Torralba | Small codes and large databases for recognition[END_REF][START_REF] Weiss | Spectral hashing[END_REF][START_REF] Jégou | Improving bag-of-features for large scale image search[END_REF][START_REF] Jégou | Searching in one billion vectors: re-rank with source coding[END_REF], because the indexed representation must fit in memory to avoid costly hard-drive accesses. One popular way is to use a Hamming Embedding function that maps the real vectors into binary vectors [START_REF] Torralba | Small codes and large databases for recognition[END_REF][START_REF] Weiss | Spectral hashing[END_REF][START_REF] Jégou | Improving bag-of-features for large scale image search[END_REF]: Binary vectors are compact, and searching the Hamming space is efficient (XOR operation and bit count) even if the comparison is exhaustive between the binary query and the database vectors. An extension to these techniques is the asymmetric scheme [START_REF] Dong | Asymmetric distance estimation with sketches for similarity search in high-dimensional spaces[END_REF][START_REF] Jain | Asymmetric hamming embedding[END_REF] which limits the approximation done on the query, leading to better results for a slightly higher complexity.

We propose to address the ANN search problem with an antisparse solution based on the design of spread representations recently proposed by Fuchs [START_REF] Fuchs | Spread representations[END_REF]. Sparse coding has received in the last decade a huge attention from both theoretical and practical points of view. Its objective is to represent a vector in a higher dimensional space with a very limited number of non-zeros components. Anti-sparse coding has the opposite properties. It offers a robust representation of a vector in a higher dimensional space with all the components sharing evenly the information.

Sparse and anti-sparse coding admits a common formulation. The algorithm proposed by Fuchs [START_REF] Fuchs | Spread representations[END_REF] is indeed similar to pathfollowing methods based on continuation techniques like [START_REF] Efron | Least angle regression[END_REF]. The anti-sparse problem considers a ℓ∞ penalization term where the sparse problem usually considers the ℓ1 norm. The penalization in x ∞ limits the range of the coefficients which in turn tend to 'stick' their value to ± x ∞ [START_REF] Fuchs | Spread representations[END_REF]. As a result, the anti-sparse approximation offers a natural binarization method.

Most importantly and in contrast to other Hamming Embedding techniques, the binarized vector allows an explicit and reliable reconstruction of the original database vector. This reconstruction is very useful to refine the search. First, the comparison of the Hamming distances between the binary representations identifies some potential nearest neighbors. Second, this list is refined by computing the Euclidean distances between the query and the reconstructions of the database vectors.

We also provide a Matlab package to reproduce the analysis comparisons reported in this paper 1 . The paper is organized as follows. Section 2 introduces the anti-sparse coding framework. Section 3 describes the corresponding ANN search method which is evaluated in Section 4 on both synthetic and real data.

SPREAD REPRESENTATIONS

This section briefly describes the anti-sparse coding of [START_REF] Fuchs | Spread representations[END_REF]. We first introduce the objective function and provide the guidelines of the algorithm giving the spread representation of a given input real vector.

Let A = [a1| . . . |am] be a d × m (d < m) full rank matrix. For any y ∈ R d , the system Ax = y admits an infinite number of solutions. To single out a unique solution, one add a constraint as for instance seeking a minimal norm solution. Whereas the case of the Euclidean norm is trivial, and the case of the ℓ1-norm stems in the vast literature of sparse representation, Fuchs recently studied the case of the ℓ∞-norm. Formally, the problem is:

x ⋆ = min x: Ax=y x ∞, (1)
with x ∞ = max i∈{1,...,m} |xi|. Interestingly, by minimizing the range of the components, md + 1 of them are stuck to the limit, ie. xi = ± x ∞ (a proof will be presented in a journal version). Fuchs also exhibits an efficient way to solve [START_REF] Casey | Content-based music information retrieval: Current directions and future challenges[END_REF]. He proposes to solve the series of simpler problems

x ⋆ h = min x∈R m J h (x) with J h (x) = Ax -y 2 2 /2 + h x ∞
for some decreasing values of h. As h → 0, x ⋆ h → x ⋆ . The reader may refer to the simple Matlab implementation we provide for an algorithmic presentation of the method.

The sub-differential set

For a fixed h, J h is not differentiable due to . ∞. Therefore, we need to work with sub-differential sets. The sub-differential set ∂f

(x) of function f at x is the set of gradients v s.t. f (x ′ )-f (x) ≥ v ⊤ (x ′ -x), ∀x ′ ∈ R m . For f ≡ . ∞, we have: ∂f (0) = {v ∈ R m : v 1 ≤ 1}, (2) ∂f (x) = {v ∈ R m : v 1 = 1, (3) vixi ≥ 0 if |xi| = x ∞, vi = 0 else} , for x = 0 Since J h is convex, x ⋆ h is solution iff 0 belongs to the sub-differential set ∂J h (x ⋆ h ), i.e. iff there exist v ∈ ∂f (x ⋆ h ) s.t. A ⊤ (Ax ⋆ h -y) + hv = 0 (4)

Initialization and first iteration

For h0 large enough, J h 0 (x) is dominated by x ∞, and the solution writes

x ⋆ h 0 = 0 and v = h -1 0 A ⊤ y ∈ ∂f (0). (2)
shows that this solution no longer holds for h < h1 with h1 = A ⊤ y 1.

For x ∞ small enough compared to the norm of 4) by sign(v) ⊤ on the left and substituting x * h , we have

A ⊤ y, J h (x) is dominated by y 2 -x ⊤ A ⊤ y + h x ∞ whose minimizer is x ⋆ h = x ∞sign(A ⊤ y). In this case, ∂f (x) is the set of vectors v s.t. sign(v) = sign(x) and v 1 = 1. Multiplying (
h = h1 -Asign(A ⊤ y) 2 x ∞. (5)
This shows that i) x ⋆ h can be a solution for h < h1, and ii) x ∞ increases as h decreases. Yet, Eq. ( 4) also imposes that v = ν1µ 1 x ∞, with ν1 h -1 A ⊤ y and µ 1 h -1 A ⊤ Asign(A ⊤ y). But, the condition sign(v) = sign(x) from (3) must hold, ie. no change of sign is allowed in v. This limits x ∞ by ρi 2 where ρi = ν1,i/µ1,i and i2 = arg mini:ρ i >0(ρi), which in turn translates to a lower bound h2 on h via (5).

Index partition

For the sake of simplicity, we introduce I {1, . . . , m}, and the index partition Ī {i : |xi| = x ∞} and Ȋ I \ Ī. The restriction of vectors and matrices to Ī (resp. Ȋ) are denoted alike x (resp. x). For instance, Eq. (3) translates in sign(v) = sign(x), v 1 = 1 and v = 0. The index partition splits (4) into two parts:

Ȃ⊤ Ȃx + Āsign(v) x ∞ = Ȃ⊤ y (6) Ā⊤ Ȃx + Āsign(v) x ∞ -y = -hv (7)
For h2 ≤ h < h1, we've seen that x = x, v = v, and Ā = A. Their 'tilde' versions are empty. For h < h2, the index partition Ī = I and Ȋ = ∅ can no longer hold. Indeed, when v1,i 2 is null at h = h2, the i2-th column of A moves from Ā to Ȃ s.t. now, Ȃ = [ai 2 ].

General iteration

The general iteration consists i) in determining on which interval [h k+1 , h k ] an index partition holds, ii) giving the expression of the solution x ⋆ h on this interval, iii) and proposing a new index partition to the next iteration. Provided Ȃ is full rank, (6) gives

x = ξ k + ζ k x ∞, with (8)
ξ k = ( Ȃ⊤ Ȃ) -1 Ȃ⊤ y and ζ k = -( Ȃ⊤ Ȃ) -1 Āsign(v). (7) gives: v = ν k -µ k x ∞, with (9)
µ k = Ā⊤ (I -Ā⊤ Ȃ( Ȃ⊤ Ȃ) -1 ) Āsign(v)/h and ν k = ( Ȃ⊤ yξ k )/h. Left multiplying [START_REF] Dong | Asymmetric distance estimation with sketches for similarity search in high-dimensional spaces[END_REF] by sign(v), we get:

h = η k -υ k x ∞ (10)
with

υ k = ( Āsign(v)) ⊤ I -Ȃ( Ȃ⊤ Ȃ) -1 Ȃ⊤ Āsign(v), and
η k = -sign(v) ⊤ Ā⊤ ( Ȃx -y). Note that υ k > 0 so that x ∞ increases when h decreases.
These equations extend a solution x ⋆ h to the neighborhood of h. However, we must check that this index partition is still valid as we decrease h and x ∞ increases. Two events can break the validity:

• Like in the first iteration, a component of v given in (9) becomes null, which breaks sign(v) = sign(x). This index moves from Ī to Ȋ. • A component of x given in (8) sees its amplitude equalling ± x ∞. This index moves from Ȋ to Ī, and the sign of this component will be the sign of the new component of x.

The value of x ∞ for which one of these two events first happens is translated in h k+1 thanks to (10).

Stopping condition and output

If the goal is to minimize J h t (x) for a specific target ht, then the algorithm stops when h k+1 < ht. The real value of x ⋆ h t ∞ is given by [START_REF] Efron | Least angle regression[END_REF], and the components not stuck to ± x ⋆ h t ∞ by [START_REF] Jain | Asymmetric hamming embedding[END_REF]. We obtain the spread representation x of the input vector y. The vector x has many of its components equal to ± x ∞ (md + 1 indeed if ht → 0). An approximation of the original vector y is

ŷ = Ax. ( 11
)
Remark: Alternately, the termination rule based on the pre-defined target value ht may be replaced by another one, such as using a fixed number of iterations. Empirically, this choice leads to comparable results with respect to the approximate search problem.

INDEXING AND SEARCH MECHANISMS

This section describes how Hamming Embedding functions are used for approximate search, and in particular how the anti-sparse coding framework described in Section 2 is exploited.

Problem statement

Let Y be a dataset of n real vectors, Y = {y1, . . . , yn}, where yi ∈ R d , and consider a query vector q ∈ R d . We aim at finding the k vectors in Y that are closest to the query, with respect to the Euclidean distance. For the sake of exposure, we consider without loss of generality the nearest neighbor problem, i.e., the case k = 1.

The nearest neighbor of q in Y is defined as

NN(q) = arg min y∈Y q -y 2 . ( 12
)
The goal of approximate search is to find this nearest neighbor with high probability and using as less resources as possible. The performance criteria are the following:

• The quality of the search, i.e., to which extent the algorithm is able to return the true nearest neighbor ;

• The search efficiency, typically measured by the query time ;

• The memory usage, i.e., the number of bytes used to index a vector yi of the database.

In our paper, we assess the search quality by the recall@R measure: over a set of queries, we compute the proportion for which the system returns the true nearest neighbor in the first R positions.

Approximate search with binary embeddings

A class of ANN methods is based on embedding [START_REF] Torralba | Small codes and large databases for recognition[END_REF][START_REF] Weiss | Spectral hashing[END_REF][START_REF] Jégou | Improving bag-of-features for large scale image search[END_REF]. The idea is to map the input vectors to a space where the representation is compact and the comparison is efficient. The Hamming space offers these two desirable properties. The key problem is the design of the embedding function e : R d → B m mapping the input vector y to b = e(y) in the m-dimensional Hamming space B m , here defined as {-1, 1} m for the sake of exposure. Once this function is defined, all the database vectors in Y are mapped to B m , and so is the query vector q. The search problem is translated into the Hamming space based on the Hamming distance, or, equivalently:

NN b (e(q)) = arg max y∈Y e(q) ⊤ e(y).

NN b (e(q)) is returned as the approximate NN(q).

Binarization with anti-sparse coding. Given an input vector y, the anti-sparse coding of Sect. 2 (with fixed parameters A and ht) produces x with many components equal to ±||x||∞. We consider a "pre-binarized" version ẋ(y) = x/ x ∞, and the binarized version e(y) = sign(x).

Hash function design

The locality sensitive hashing (LSH) algorithm is mainly based on random projection, though different kinds of hash functions have been proposed for the Euclidean space [START_REF] Andoni | Near-optimal hashing algorithms for near neighbor problem in high dimensions[END_REF].

Let A = [a1| . . . |am] be a d × m matrix storing the m projection vectors. The most simple way is to take the sign of the projections: b = sign(A ⊤ y). Note that this corresponds to the first iteration of our algorithm (see Sect. 2.2). We also try A as an uniform frame. A possible construction of such a frame consists in performing a QR decomposition on a m×m matrix. The matrix A is then composed of the d first rows of the Q matrix, ensuring that A × A ⊤ = I d . Sect. 4 shows that such frames significantly improve the results compared with random projections, for both LSH and anti-sparse coding embedding methods.

Asymmetric schemes

As recently suggested in the literature, a better search quality is obtained by avoiding the binarization of the query vector. Several variants are possible. We consider the simplest one derived from (13), where the query is not binarized in the inner product. For our antisparse coding scheme, this amounts to performing the search based on the following maximization:

NNa (e(q)) = arg max y∈Y ẋ(q) ⊤ e(y).

(

) 14
The estimate NNa is better than NN b . The memory usage is the same because the vectors in the database {e(yi)} are all binarized. For better efficiency, the search ( 14) is done using look-up tables computed for the query and prior to the comparisons. Doing so, the asymmetric scheme remains a bit slower than the pure bit-based comparison: comparing two optimized implementations, Jain et al. report [START_REF] Jain | Asymmetric hamming embedding[END_REF] that the comparison of binary vectors is 1.7 times faster than the corresponding asymmetric scheme on a typical computer architecture. This is slightly slower than computing the Hamming distances in (13). This asymmetric scheme is interesting for any binarization scheme (LSH or anti-sparse coding) and any definition of A (either random projections or a frame).

Explicit reconstruction

The anti-sparse binarization scheme explicitly minimizes the reconstruction error, which is traded in (1) with the ℓ∞ regularization term. Eq. ( 11) gives an explicit approximation of the database vector y up to a scaling factor: ŷ ∝ Ae(y) ||Ae(y)|| 2 . The approximate nearest neighbors NNe are obtained by computing the exact Euclidean distances ||q -ŷi||2. This is slow compared to the Hamming distance computation. That is why, it is used to operate, like in [START_REF] Jégou | Searching in one billion vectors: re-rank with source coding[END_REF], a re-ranking of the first hypotheses returned based on the Hamming distance (on the asymmetric scheme described in Sect. 3.3). The main difference with [START_REF] Jégou | Searching in one billion vectors: re-rank with source coding[END_REF] is that no extra-code has to be retrieved: the reconstruction ŷ solely relies on e(y).

SIMULATIONS AND EXPERIMENTS

This section evaluates the search quality on synthetic and real data. In particular, we measure the impact of:

• The Hamming embedding technique: LSH and binarization based on anti-sparse coding. We also compare to the spectral hashing method of [START_REF] Weiss | Spectral hashing[END_REF], using the code available online. • The choice of matrix A: random projections or frame for LSH. For the anti-sparse coding, we always assume a frame. • The search method: 1) NN b of (13) 2) NNa of ( 14) and 3) NNe as described in Sect. 3.4.

Our comparison focuses on the case m ≥ d. In the anti-sparse coding method, the regularization term h controls the trade-off between the robustness of the Hamming embedding and the quality of the reconstruction. Small values of h favors the quality of the reconstruction (without any binarization). Bigger values of h gives more components stuck to x ∞, which improves the approximation search with binary embedding. Optimally, this parameter should be adjusted to give a reasonable trade-off between the efficiency of the first stage (methods NN b or NNa) and the re-ranking stage (NNe). Note however that, thanks to the algorithm described in Section 2, the parameter is stable, i.e., a slight modification of this parameter only affects a few components. We set h = 1 in all our experiments. Two datasets are considered for the evaluation:

• A database of 10,000 16-dimensional vectors uniformly drawn on the Euclidean unit sphere (normalized Gaussian vectors) and a set of 1,000 query vectors. • A database of SIFT [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] descriptors available online2 , comprising 1 million database and 10,000 query vectors of dimensionality 128. Similar to [START_REF] Weiss | Spectral hashing[END_REF], we first reduce the vector dimensionality to 48 components using principal component analysis (PCA). The vectors are not normalized after PCA.

The comparison of LSH and anti-sparse. Figures 1 and2 show the performance of Hamming embeddings for synthetic data. On Fig. 1, the quality measure is the recall@10 (proportion of true NN ranked in first 10 positions) plotted as a function of the number of bits m. For LSH, observe the much better performance obtained by the proposed frame construction compared with random projections. The same conclusion holds for anti-sparse binarization. The anti-sparse coding offers similar search quality as LSH for m = d when the comparison is performed using NN b of (13). The improvement gets significant as m increases. The spectral hashing technique [START_REF] Weiss | Spectral hashing[END_REF] exhibits poor performance on this synthetic dataset.

The asymmetric comparison NNa leads a significant improvement, as already observed in [START_REF] Dong | Asymmetric distance estimation with sketches for similarity search in high-dimensional spaces[END_REF][START_REF] Jain | Asymmetric hamming embedding[END_REF]. The interest of anti-sparse coding becomes obvious by considering the performance of the comparison NNe based on the explicit reconstruction of the database vectors from their binary-coded representations. For a fixed number of bits, the improvement is huge compared to LSH. It is worth using this technique to re-rank the first hypotheses obtained by NN b or NNa.

Experiments on SIFT descriptors. As shown by Figure 3, LSH is slightly better than anti-sparse on real data when using the binary representation only (here m = 128), which might solved by tuning h, since the first iteration of antisparse leads the binarization as LSH. However, the interest of the explicit reconstruction offered by NNe is again obvious. The final search quality is significantly better than that obtained by spectral hashing [START_REF] Weiss | Spectral hashing[END_REF]. Since we do not specifically handle the fact that our descriptor are not normalized after PCA, our results could probably be improved by taking care of the ℓ2 norm.

CONCLUSION AND OPEN ISSUES

In this paper, we have proposed anti-sparse coding as an effective Hamming embedding, which, unlike concurrent techniques, offers an explicit reconstruction of the database vectors. To our knowledge, it outperforms all other search techniques based on binarization. There are still two open issues to take the best of the method. First, the computational cost is still a bit high for high dimensional vectors. Second, if the proposed codebook construction is better than random projections, it is not yet specifically adapted to real data.

Fig. 1 .Fig. 2 .

 12 Fig.1. Anti-sparse coding vs LSH on synthetic data. Search quality (recall@10) as a function of the number of bits of the representation.

Fig. 3 .

 3 Fig. 3. Approximate search in a SIFT vector set of 1 million vectors.

This work was realized as part of the Quaero Project, funded by OSEO, French State agency for innovation.

http://gforge.inria.fr/projects/antisparse

http://corpus-texmex.irisa.fr


Lionel Bombrun
email: lionel.bombrun@ims-bordeaux.fr

Nour- Eddine Lasmar
email: nour-eddine.lasmar@ims-bordeaux.fr

Yannick Berthoumieu
email: yannick.berthoumieu@ims-bordeaux.fr

Geert Verdoolaege
email: geert.verdoolaege@ugent.be

MULTIVARIATE TEXTURE RETRIEVAL USING THE SIRV REPRESENTATION AND THE GEODESIC DISTANCE

Keywords: Texture, Multiscale analysis, Kullback-Leibler divergence, Geodesic distance, Spherically Invariant Random Vector

This paper presents a new wavelet based retrieval approach based on Spherically Invariant Random Vector (SIRV) modeling of wavelet subbands. Under this multivariate model, wavelet coefficients are considered as a realization of a random vector which is a product of the square root of a scalar random variable (called multiplier) with an independent Gaussian vector. We propose to work on the joint distribution of the scalar multiplier and the multivariate Gaussian process. For measuring similarity between two texture images, the geodesic distance is provided for various multiplier priors. A comparative study between the proposed method and conventional models on the VisTex image database is conducted and indicates that SIRV modeling combined with geodesic distance achieves higher recognition rates than classical approaches.

INTRODUCTION

Several works in textured image indexing have shown that the wavelet representation is a well adapted domain to characterize the texture [START_REF] Do | Wavelet-Based Texture Retrieval Using Generalized Gaussian Density and Kullback-Leibler Distance[END_REF] [2] [3] [4]. This leads to a multiscale analysis scheme which consists in modeling each wavelet subband by a given probability density function (PDF). The distribution parameters are then estimated and compose the signature of the texture while a probabilistic metric is used to measure similarity. Previous works [START_REF] Do | Wavelet-Based Texture Retrieval Using Generalized Gaussian Density and Kullback-Leibler Distance[END_REF] [START_REF] Kwitt | Lightweight Probabilistic Texture Retrieval[END_REF] show that the symmetrized Kullback-Leibler divergence (SKLD, called also Jeffrey divergence), increases significantly the retrieval rate in the framework of stochastic models. Although the SKLD is a popular similarity measure, it is not a distance. The geodesic distance (GD) [START_REF] Berkane | Geodesic Estimation in Elliptical Distributions[END_REF] which relies on the information metric should increase the retrieval rate [START_REF] Verdoolaege | Multiscale Colour Texture Retrieval using the Geodesic Distance between Multivariate Generalized Gaussian Models[END_REF].

Many univariate models such as the Generalized Gaussian (GG) [START_REF] Do | Wavelet-Based Texture Retrieval Using Generalized Gaussian Density and Kullback-Leibler Distance[END_REF] and the Weibull (Wbl) [START_REF] Kwitt | Lightweight Probabilistic Texture Retrieval[END_REF] distributions have been introduced to successfully characterize the wavelet subbands.

Those univariate approaches are quite simple but they do not exploit the spatial dependency inside wavelet subbands. Multivariate models have consequently been studied. The multivariate Gaussian mixture (MGmix) was introduced for texture classification and segmentation [START_REF] Kim | Texture Classification and Segmentation using Wavelet Packet Frame and Gaussian Mixture Model[END_REF] and the SKLD on MGmix models can be approximated using the works of [START_REF] Hershey | Approximating the Kullback Leibler Divergence Between Gaussian Mixture Models[END_REF]. In [START_REF] Lasmar | Multivariate Statistical Modeling for Texture Retrieval Using Wavelet Transforms[END_REF], a SIRV based modeling has been proposed and a closed form of SKLD has been derived. Recently, the GD based on the multivariate Generalized Gaussian (MGG) model [START_REF] Verdoolaege | Multiscale Colour Texture Retrieval using the Geodesic Distance between Multivariate Generalized Gaussian Models[END_REF] has been proposed. It leads to higher performance recognition rates compared to the SKLD, but a GD closed form is given only for fixed MGG shape parameters.

In this paper, SIRV modeling is proposed to represent wavelet subbands and a closed form of GD is derived. Retrieval performances are compared to those of SKLD. The remainder of this paper is as follows. We introduce in section 2 the SIRV based multivariate modeling. In section 3, expressions of SKLD and GD on joint probability resulting from SIRV representation are developed. Experimental results are given in section 4 to evaluate retrieval performance. Finally, section 5 concludes the paper and suggests an outlook on future works.

SIRV MODELING

SIRV modeling has been successfully used for characterizing non-Gaussian stochastic processes such as radar clutter returns, radio fading analysis, or sonar interferences. Recently, in image processing field, SIRV representation has been employed for SAR image segmentation [START_REF] Bombrun | Hierarchical Segmentation of Polarimetric SAR Images using Heterogeneous Clutter Models[END_REF] and texture retrieval [START_REF] Lasmar | Multivariate Statistical Modeling for Texture Retrieval Using Wavelet Transforms[END_REF].

Let x s,o (n, m) be the wavelet coefficient at scale s, orientation o and centered at spatial location (2 s n, 2 s m). In the multivariate modeling context, a wavelet subband is considered as a realization of a field characterized by a spatial dependency following a given neighborhood around the current spatial location and clustered in a random vector x.

In this paper, the vector x is modeled as a SIRV which is defined as:

x = √ τ g ( 1
)
where τ is a scalar random variable called the multiplier (τ ∈ R + ) and g a real Gaussian vector with zero mean and covariance matrix [M] = E{ g g T }. The superscript T denotes the transpose operator and E{•} the expectation. The processes τ and g are supposed independant. We denote d = p × q the size of the considered neighborhood, and p τ (τ |Θ) the PDF of the scalar multiplier. The PDF of the multivariate vector x = √ τ g is obtained using Bayes' theorem:

pX ( x) = 1 (2π) d 2 |[M]| 1 2 Z +∞ 0 1 τ d 2 exp " - x T [M] -1 x 2τ " pτ (τ |Θ) dτ = 1 |[M]| 1 2 h d " x T [M] -1 x " (2)
where h d (•) is the density generator. An example of the SIRV decomposition is given in Fig. 1. The multiplier τ and the Gaussian vector g are extracted from the wavelet subband.

In literature, various model are used to represent the multiplier τ . The Gamma distribution is at the basis of the multiplier models used in [START_REF] Quegan | Statistical Models for Polarimetric SAR Data[END_REF]. Other a priori could be used, such as Inverse Gamma [START_REF] Freitas | The Polarimetric G Distribution for SAR Data Analysis[END_REF], Weibull [START_REF] Lasmar | Multivariate Statistical Modeling for Texture Retrieval Using Wavelet Transforms[END_REF] or Fisher distributions [START_REF] Bombrun | Hierarchical Segmentation of Polarimetric SAR Images using Heterogeneous Clutter Models[END_REF]. No closed form exists for the PDF of the multivariate vector x for the case of Weibull a priori but for the three other models (Gamma, Inverse Gamma and Fisher) the multivariate distribution can be derived. They are respectively the K [9], G 0 [START_REF] Freitas | The Polarimetric G Distribution for SAR Data Analysis[END_REF] and KummerU [START_REF] Bombrun | Hierarchical Segmentation of Polarimetric SAR Images using Heterogeneous Clutter Models[END_REF] PDFs. Their corresponding density generators are given by: Unfortunately, in the best of our knowledge there is no analytic expression neither for KLD or GD between these multivariate PDFs. However, characterizing texture content can be provided considering vector y = (τ, g) resulting from the SIRV representation. In this case and using the independence between τ and g, the joint PDF of vector y is

h K d (u) = 1 (2π) d 2 2 Γ(L) L m L+ d 2 2 u 2 L-d 2 2 K d 2 -L 2 Lu 2m h G 0 d (u) = 1 (2π) d 2 (Mm) M Γ(M) Γ d 2 + M u 2 + Mm -(M+ d 2 ) h KummerU d (u) = 1 (2π) d 2 1 B(L, M) L Mm d 2 Γ d 2 + M × U d 2 + M; 1 + d 2 -L; L 2Mm u
p Y ( y) = p τ (τ |Θ)p G ( g) (3)
Let x i , i = 1 . . . N , be the realizations of the vector x. The Approximated Maximum Likelihood (AML) estimator of the normalized covariance matrix is the solution of the following recursive equation:

[ M] FP = f ([ M] FP ) = d N N i=1 x i x T i x T i [ M] -1 FP x i (4)
The existence and the uniqueness, up to a scalar factor, of the Fixed Point estimator of the normalized covariance matrix have been established, as well as the convergence of the recursive algorithm, irrespective of the initialization [START_REF] Pascal | Covariance Structure Maximum-Likelihood Estimates in Compound Gaussian Noise : Existence and Algorithm Analysis[END_REF]. In this paper, the trace of the covariance matrix is normalized to d, the dimension of vector x.

For a given covariance matrix [M], the ML estimator of the multiplier for the pixel i is given by:

τi = x T i [M] -1 x i d . (5)
To complete the hyperparameters estimation, the prior on τ is used and all the expressions of ML estimators for Weibull, Gamma, Inverse Gamma and Fisher can be found in [START_REF] Johnson | Continuous Univariate Distributions[END_REF].

SIMILARITY MEASURES

Kullback-Leibler divergence

By working on the vector y = (τ, g), the Kullback-Leibler divergence (KLD) is obtained using the chain rule, since the multiplier parameter τ and the Gaussian process g are independent in the SIRV model. This yields [START_REF] Lasmar | Multivariate Statistical Modeling for Texture Retrieval Using Wavelet Transforms[END_REF]:

KLD " pY ( y)||qY ( y) " = KLD " pτ (τ )||qτ (τ ) " + KLD " pG( g)||qG( g) " (6)
The first term in [START_REF] Hershey | Approximating the Kullback Leibler Divergence Between Gaussian Mixture Models[END_REF] corresponds to the KLD between the multiplier models, while the second term corresponds to the KLD for the multivariate Gaussian process:

KLD p G ( g)||q G ( g) = 1 2 tr [M] -1 2 [M] 1 +ln |[M] 2 | |[M] 1 | -d (7)
The KLDs between two models depend to the multiplier τ prior and are given for Weibull and Gamma priors respectively in [START_REF] Lasmar | Multivariate Statistical Modeling for Texture Retrieval Using Wavelet Transforms[END_REF] and [START_REF] Mathiassen | Texture Similarity Measure using Kullback-Leibler Divergence Between Gamma Distributions[END_REF]. For a Fisher distributed random variable, the KLD has an analytical expression:

KLD " p(τ |m1, L1, M1)||q(τ |m2, L2, M2) " = ln B(L2, M2) B(L1, M1) ! + L1 ln " L1 M1m1 " -L2 ln " L2 M2m2 " +(L1 -L2) h ln m1 + (Ψ(L1) -ln L1) -(Ψ(M1) -ln M1) i -(L1 + M1) h Ψ(L1 + M1) -Ψ(M1) i + (L2 + M2) B(L1, M1) ∞ Z 0 x L 1 -1 (1 + x) L 1 +M 1 ln " 1 + L2M1m1 M2m2L1 x « dx (8)
Similarly, the KLD for an Inverse Gamma PDF is given by:

KLD p(τ |m 1 , M 1 )||q(τ |m 2 , M 2 ) = ln Γ(M 2 ) Γ(M 1 ) + M 1 ln(M 1 m 1 ) -M 2 ln(M 2 m 2 ) -M 1 + M 2 m 2 m 1 -(M 1 -M 2 ) ln m 1 -(Ψ(M 1 ) -ln(M 1 )) (9)
In the following, the Jeffrey divergence will be considered, i.e. half of the double-sided KLD:

SKLD = KLD p Y ( y)||q Y ( y) + KLD q Y ( y)||p Y ( y) 2 (10)

Rao geodesic distance

The geodesic distance (GD), which is based on the Fisher information metric, can be used as an alternative of the KLD. As the multiplier τ and the Gaussian process g are independent in the SIRV model, the GD of the joint distribution p Y ( y) is given by: GD " pY ( y)||qY ( y)

" = GD " pτ (τ )||qτ (τ ) " + GD " pG( g)||qG( g) " (11)
The second term of ( 11) corresponds to the GD on the Gaussian part. Its expression is given by [START_REF] Berkane | Geodesic Estimation in Elliptical Distributions[END_REF].

GD p G ( g|[M] 1 )||p G ( g|[M] 2 ) = 1 2 d i=1 (ln λ i ) 2 1 2 (12)
where

λ i (i = 1, . . . , d) are the eigenvalues of [M] -1 1 [M] 2 .
The first term of [START_REF] Pascal | Covariance Structure Maximum-Likelihood Estimates in Compound Gaussian Noise : Existence and Algorithm Analysis[END_REF] corresponds to the GD on the multiplier variable. For various stochastic models (Gamma, Weibull, . . .), the GD can be computed numerically. Nevertheless, for distributions that lie infinitesimally close on the probabilistic manifold, it can be proved that the KLD equals half of the squared GD between the distributions. It follows that locally [START_REF] Kullback | Information Theory and Statistics[END_REF]:

GD p(τ )||q(τ ) ≈ 2 SKLD p(τ )||q(τ ) (13)

RESULTS

Several experiments are conducted to evaluate the performance of the proposed algorithms based on the SKLD and the GD between multivariate random variables. We work with the 40 texture classes from the MIT Vision Texture database.

From each of these texture images of size 512 × 512 pixels, 16 subimages (128 × 128) are created. A dataset containing 640 texture images is finally obtained [START_REF] Do | Wavelet-Based Texture Retrieval Using Generalized Gaussian Density and Kullback-Leibler Distance[END_REF] [2].

Table 1 shows the average retrieval rate for the joint distribution p Y ( y). This analysis is carried out for the SKLD [START_REF] Hershey | Approximating the Kullback Leibler Divergence Between Gaussian Mixture Models[END_REF] and the GD [START_REF] Pascal | Covariance Structure Maximum-Likelihood Estimates in Compound Gaussian Noise : Existence and Algorithm Analysis[END_REF] with four multiplier models, i.e. Weibull, Gamma, Fisher and Inverse Gamma. Three wavelet decompositions are analyzed: the steerable pyramid with 2 scales (Nsc) and 4 orientations (Nor), the orthogonal wavelet transform (OWT, Nsc = 2) with Daubechies' filter db2 and the dual-tree complex wavelet transform (DT-CWT, Nsc = 2). First, the GD outperforms the SKLD in all our experiments (increases average retrieval rate by about 2%). Moreover, the choice of the PDF to model the multiplier does not seem critical. All multiplier PDFs yield similar performance (GD ≈ 82.5% for the steerable pyramid), except the Inverse Gamma distribution (78.33%) which does not seem well-adapted for the modeling of the multiplier. The use of the Dual-tree complex wavelet transform with two scales outperforms the two other wavelet decompositions. The best retrieval rates are obtained for the GD with the Weibull and Gamma models (GD ≈ 86.35%).

SIRV based model

Fig. 2 shows the performances of the proposed method using the joint distribution p Y ( y) with the GD for a Gamma distributed multiplier (SIRV-Gamma+GD). Comparisons are done with classical approaches based on the SKLD. Moreover, ROC curves provide an additional measure of the improved retrieval performance using the joint distribution p Y ( y) with the GD. For example, we need to use queries of size 24 to reach 90% of relevant samples while we must retrieve 35, 48 and 49 images to reach the same percentage with multivariate Gaussian mixture (GMix + SKLD), the univariate Weibull (Wbl + SKLD) and Generalized Gaussian (GG + SKLD) models respectively.

CONCLUSION

In this paper, we have proposed a multiscale texture retrieval algorithm using the GD between multivariate PDFs. The SIRV modeling has been used to characterize wavelet coefficients. The joint distribution of the multiplier and the Gaussian part resulting from SIRV representation is employed, and the GD has been derived. The proposed method outperforms the widely used univariate models and the multivariate Gaussian mixture model. Further works will concern the derivation of the GD based on the SIRV density generators, which should improve the retrieval performance.

Fig. 1 .

 1 Fig. 1. Example of the SIRV decomposition on a wavelet subband.

  where Γ(•) and B(•, •) are the Gamma and Beta functions. K(•) and U(•; •; •) are the modified Bessel function, resp. the confluent hypergeometric function of the second kind.

Fig. 2 .

 2 (a) plots the recall/precision curves for the four different models. It shows that the proposed algorithm based on the joint distribution p Y ( y) significantly improves the recognition rate.

Fig. 2 .

 2 Fig. 2. Comparison between the proposed method and conventional algorithms.

Table 1 .

 1 Average retrieval rate (in %) by using the joint distribution p Y ( y).

	Wavelet decomposition	Similarity measure	Weibull	Gamma	Fisher	Inverse Gamma
	Steerable pyramid	SKLD	79.19	79.39	78.69	74.27
	Nsc = 2, Nor = 4	GD	82.73	82.97	82.47	78.33
	OWT	SKLD	80.39	80.41	79.96	76.40
	db2, Nsc = 2	GD	83.89	83.90	83.69	79.02
	DT-CWT	SKLD	82.08	82.20	81.52	71.65
	Nsc = 2	GD	86.38	86.34	85.97	82.16


Lionel Bombrun
email: lionel.bombrun@ims-bordeaux.fr

Yannick Berthoumieu
email: yannick.berthoumieu@ims-bordeaux.fr

Nour-Eddine Lasmar
email: nour-eddine.lasmar@ims-bordeaux.fr

Geert Verdoolaege
email: geert.verdoolaege@ugent.be

MULTIVARIATE TEXTURE RETRIEVAL USING THE GEODESIC DISTANCE BETWEEN ELLIPTICALLY DISTRIBUTED RANDOM VARIABLES

Keywords: Texture, Geodesic distance, Multivariate elliptical distribution

published or not. The documents may come

Multivariate texture retrieval using the geodesic distance between elliptically distributed random variables

Lionel Bombrun, Yannick Berthoumieu, Nour-Eddine Lasmar, Geert Verdoolaege

INTRODUCTION

In the context of texture image recognition, the wavelet representation domain has been widely considered to characterize the texture. The multiscale recognition algorithm is composed by two main steps. One is the feature extraction which consists in modeling each wavelet subband by a given probability density function (pdf). Its estimated parameters compose the signature of the texture. Then, a similarity measure based on a probabilistic metric is computed between the query image and all the other images of the database.

Many univariate models such as the Generalized Gaussian [START_REF] Do | Wavelet-Based Texture Retrieval Using Generalized Gaussian Density and Kullback Leibler Distance[END_REF] and the Weibull [START_REF] Kwitt | Lightweight Probabilistic Texture Retrieval[END_REF] distributions have been proposed to characterize the wavelet subbands. Those approaches are quite simple but they do not fully exploit the texture information. Multivariate distributions have recently been introduced to model the spatial and/or color dependencies of the wavelet coefficients such as the Gaussian [START_REF] Stitou | Copulas based Multivariate Gamma Modeling for Texture Classification[END_REF] and Student-t [START_REF] Kwitt | A Joint Model of Complex Wavelet Coefficients for Texture Retrieval[END_REF] copulas, and the multivariate elliptical distributions.

In this paper, two elliptical distributions are considered to model the color and spatial dependencies of the wavelet coefficients: the multivariate Laplace [START_REF] Verdoolaege | Multiscale Colour Texture Retrieval using the Geodesic Distance between Multivariate Generalized Gaussian Models[END_REF] and G 0 distributions. For the G 0 distribution, a closed form of the geodesic distance (GD) [START_REF] Berkane | Geodesic Estimation in Elliptical Distributions[END_REF] is established for fixed shape parameters and an approximation of the GD is given for the general case by assuming the geodesic coordinate functions as straight lines. A multi-model approach is then proposed for the classification of texture images. Next, some experimental results on the VisTex database are shown in section 3 to evaluate the retrieval performance. Finally, some conclusions of this work are given.

MULTIVARIATE TEXTURE MODELING

To model the spatial or color dependencies of the wavelet coefficient, various multivariate models can be considered such as copulas [3] [4]. Here, two elliptical distributions are considered: the multivariate Laplace and multivariate G 0 distributions.

The multivariate Laplace distribution

The Multivariate Generalized Gaussian Distribution (MGGD) has been introduced in [START_REF] Verdoolaege | Multiscale Colour Texture Retrieval using the Geodesic Distance between Multivariate Generalized Gaussian Models[END_REF] to model the wavelet coefficients of color images. The Multivariate Laplace, a particular case of the MGGD, is generally considered as a reference model for color texture indexing. Its pdf is given by:

p k (k|[M ]) = 1 |[M ]| 1 2 Γ p 2 π p 2 Γ (p) 1 2 p+1 exp - 1 2 k T [M ] -1 k 1 2 , (1)
where the superscript T denotes the transpose operator and |[M ]| is the determinant of the covariance matrix [M ]. p is the dimension of the vector k. In the following, the covariance matrix [M ] is estimated using the method of moments [START_REF] Verdoolaege | Multiscale Colour Texture Retrieval using the Geodesic Distance between Multivariate Generalized Gaussian Models[END_REF].

The multivariate G 0 distribution

To model the wavelet coefficients, other multivariate models can be considered such as the Spherically Invariant Random Vectors (SIRV). This class of models has been introduced to take into account the non-Gaussianity of the signal. In such a case, the observed vector k is decomposed as:

k = √ τ z, (2)
where τ is a scalar random variable called the multiplier (τ ∈ R + ) and z a real Gaussian vector with zero mean and covariance matrix [M ] = E{zz T }. The processes τ and z are assumed independent.

In the literature, various models issued from the Pearson system have been introduced to represent the multiplier τ . For Gamma, Inverse Gamma, Fisher, Beta and Inverse Beta distributed multiplier, the SIRV vector k follows respectively a K, G 0 , KummerU, W and M distributions. An overview of those multivariate models is given in [START_REF] Bombrun | A Complete Coverage of Log-cumulant Space in Terms of Distributions for Polarimetric SAR Data[END_REF].

Let τ be an Inverse Gamma distributed random variable. Its pdf p τ (τ ) is given by:

p τ (τ ) = GI(τ |m, M) = 1 Γ(M) 1 Mm Mm τ M+1 e -Mm τ , (3)
where M ≥ 0, m > 0 and τ ∈ 0, +∞ . M is a shape parameter while m is a scale parameter. For an Inverse Gamma distributed multiplier, the vector k follows the G 0 distribution (or Student-t), its pdf is given by [START_REF] Freitas | The Polarimetric G Distribution for SAR Data Analysis[END_REF]:

p k (k|[M ], m, M) = 1 (2π) p 2 |[M ]| 1 2 (Mm) M Γ(M) Γ p 2 + M × k T [M ] -1 k 2 + Mm -(M+ p 2 )
.

Let (k 1 , . . . , k N ), be N realizations of the vector k. The Approximated Maximum Likelihood (AML) estimator of the normalized covariance matrix is the solution of the following recursive equation:

[ M ] FP = f ([ M ] FP ) = p N N i=1 z i z T i z T i [ M ] -1 FP z i . ( 5
)
The existence and the uniqueness, up to a scalar factor, of the Fixed Point estimator of the normalized covariance matrix have been established, as well as the convergence of the recursive algorithm, irrespective of the initialization [START_REF] Pascal | Covariance Structure Maximum-Likelihood Estimates in Compound Gaussian Noise : Existence and Algorithm Analysis[END_REF]. In this paper, the trace of the covariance matrix is normalized to p. For a given covariance matrix [M ], the ML estimator of the multiplier for the pixel i is given by:

τi = k T i [M ] -1 k i p . (6)
Once the N multipliers τ i are computed, the scale and shape parameters of the Inverse Gamma distribution (m and M) are estimated by using the ML method [START_REF] Johnson | Continuous Univariate Distributions[END_REF].

Rao geodesic distance

The Kullback-Leibler Divergence (KLD) has been widely used in the past to compute a similarity measure between two populations. Unfortunately, up to our knowledge, for elliptically distributed random variables, no analytical expression of the KLD exists except for the degenerate Gaussian case. But, the Rao Geodesic Distance (GD) can be used instead.

To this aim, two cases should be considered: when the shape parameters are fixed and when they aren't.

Geodesic distance for fixed shape parameters

For fixed shape parameters, the GD in the context of elliptical distributions is given by [START_REF] Berkane | Geodesic Estimation in Elliptical Distributions[END_REF]:

GD (Θ, [M ] 1 ||Θ, [M ] 2 ) = 3b h - 1 4 p i=1 (r i ) 2 + 2 b h - 1 4 p i<j r i r j 1 2 , (7)
where

r i (i = 1, . . . , p) are the eigenvalues of ln [M ] -1 1 [M ] 2 . Moreover, b h = -E{u 2 w 2 }/(p(p+2)) and w = (∂ ln h p )/(∂u).
The pdf of U = k T [M ] -1 k is given by [START_REF] Berkane | Geodesic Estimation in Elliptical Distributions[END_REF]:

f U (u) =    π p 2 Γ(p/2) u p/2-1 h p (u) if u > 0. 0 otherwise. ( 8
)
where h p (•) is the density generator:

p k (k) = 1 |[M ]| 1 2 h p k T [M ] -1 k . (9)
For the multivariate Laplace distribution, it has been shown in [START_REF] Verdoolaege | Multiscale Colour Texture Retrieval using the Geodesic Distance between Multivariate Generalized Gaussian Models[END_REF]

that b h = 1 4 p+1
p+2 . An analytical expression of the GD between two multivariate Laplace distributions is obtained as they do not involve shape parameters.

For the multivariate G 0 distribution, one can prove that

b h = 1 4 M+ p 2 M+ p 2 +1
. Consequently, for fixed shape parameters, the GD between G 0 distributions can be computed according to [START_REF] Bombrun | A Complete Coverage of Log-cumulant Space in Terms of Distributions for Polarimetric SAR Data[END_REF]. Note that when M tends toward infinity, the G 0 pdf converges toward the Gaussian pdf and the coefficient b h = 1/4 is retrieved [START_REF] Berkane | Geodesic Estimation in Elliptical Distributions[END_REF]. Note also that it does not depend on the scale parameter m.

Geodesic distance for variable m and M

The GD between two multivariate G 0 distributions with parameters (m 1 , M 1 , [M ] 1 ) and (m 2 , M 2 , [M ] 2 ) may be com-puted as follows:

GD (m 1 , M 1 , [M ] 1 ||m 2 , M 2 , [M ] 2 ) = (m2,M2,[M ]2) (m1,M1,[M ]1) ds = 1 0 µ,ν g µν θµ θν dt, (10)
where θ = m, M, r 1 , . . . , r p are the hyper-parameters of the G 0 distribution and g µν the elements of the (p+2)×(p+2) Fisher information matrix defined by:

G =           g mm g mM g mi • • • • • • g mi g mM g MM g Mi • • • • • • g Mi g mi g Mi g ii g ij • • • g ij . . . . . . g ij . . . . . . . . . . . . . . . . . . . . . . . . g ij g mi g Mi g ij • • • g ij g ii           . (11)
After some cumbersome computations, one can express the Fisher information matrix only as a function of the two shape and scale parameters:

g mm = M m 2 p 2 p 2 + M + 1 (12) g mM = - 1 m + 2M + p 2 m p 2 + M - (M + 1) m p 2 + M + 1 (13)
g MM = - 1 M + Ψ(1, M) -Ψ 1, p 2 + M + 2 p 2 + M - (M + 1) M p 2 + M + 1 ( 14
)
g mi = - (M + 2) (M + 1) M 2 4p (Mm) 3 p 2 + M + 1 p 2 -1 (15)
g Mi = - 1 2 1 M + p 2 + (M + 2) (M + 1) 4p (Mm) 2 p 2 + M + 1 p 2 -1 (16)
g ii = 1 4 3 M + p 2 M + p 2 + 1 -1 (17)
g ij = 1 4 M + p 2 M + p 2 + 1 -1 = - 1 4 1 M + p 2 + 1 , i = j (18)
where Ψ(1, •) is the trigamma function. The GD is finally computed by using [START_REF] Johnson | Continuous Univariate Distributions[END_REF] and by assuming the geodesic coordinate functions as straight lines, i.e.

   m(t) = m 1 (1 -t) + m 2 t. M(t) = M 1 (1 -t) + M 2 t. r i (t) = r i t. (19)
This is only an approximation since the geodesic itself lives on the curved manifold.

A multi-model approach

When various models are available to represent the texture images, a multi-model approach can be considered for the classification. An image x t is assigned to the label ω which maximizes the probability p(x t |ω):

ω = Argmax ω∈Ωm p(x t |ω), (20)
where Ω m = (1, . . . , m), m is the number of texture classes.

If K models are available, let:

p(x t |ω) = K i=1 p(x t |M i , ω)p(M i |ω), (21)
where M i refers to the i st model. The second term of ( 21) is the prior probability. It characterizes the probability associated to the model M i for the class ω. An approximation is:

p(M i |ω) ≡ d Mi (ω) K j=1 d Mj (ω) , (22)
where d Mi (ω) is the median distance between the images of the ω class and all other images for the M i criterion, i.e.:

d Mi (ω) = median xu∈ω,xv / ∈ω d Mi (x u , x v ). (23)
d Mi (x u , x v ) denotes the distance computed with the M i criterion between two populations x u and x v .

Similarly, we define the first term of (21) as:

p(x t |M i , ω) ≡ 1 1+d M i (xt,ω) j∈Ωm 1 1+d M i (xt,j) (24)
where d Mi (x t , j) = min xu∈j d Mi (x t , x u ).

RESULTS

Several experiments are conducted to evaluate the performances of the proposed algorithms based on the GD between multivariate random variables. We work with the 40 texture classes from the MIT Vision Texture database. From each of these texture images of size 512 × 512 pixels, 16 subimages (128 × 128) are created. A dataset containing 640 texture images is finally obtained [START_REF] Do | Wavelet-Based Texture Retrieval Using Generalized Gaussian Density and Kullback Leibler Distance[END_REF] [START_REF] Kwitt | Lightweight Probabilistic Texture Retrieval[END_REF].

In those experiments, three types of dependencies are considered (spatial, color and spatial/color). For the spatial dependency, the wavelet coefficients located around the neighborhood of the current spatial position are clustered in the random vector k spatial . For color images, the observed vector k color contains the realizations of the wavelet coefficients for each channel of the RGB image. When dealing with both spatial and color dependency, the random vector

is k spatial/color = k spatial R ; k spatial G ; k spatial B T
. In this study, the orthogonal wavelet decomposition (Nsc = 2) with Daubechies' filter db4 is considered, and a 3×3 neighborhood is used to take into account the spatial dependency.

Table 1 shows the average retrieval rate for the three considered types of dependencies and for various multivariate models (copula, Gaussian, G 0 and Laplace). For the Gaussian copula with Gamma distributed margins [START_REF] Stitou | Copulas based Multivariate Gamma Modeling for Texture Classification[END_REF], the symmetrized version of the Kullback-Leibler divergence has been considered. For the multivariate Gaussian model, the GD has been computed with both the sample covariance matrix (SCM) and fixed point covariance matrix (FP) estimators. Table 1 shows that indexing performances are better with the color models than with those characterizing the spatial dependency. A gain of about 10% is observed. By working on k spatial/color , a slight gain (around 1%) is observed compared to the retrieval rates with the color models. To model the spatial dependency, the G 0 distribution with the geodesic paths approximated as straight lines (GD linear) gives the best performances (79.38%). For color texture indexing, the best results are achieved for the GD with the G 0 model with fixed shape parameters and the multivariate Laplace distribution (≈ 87.7%).

Fig. 1 shows the average classification rate as a function of the number of training samples. For four learned images, a gain of 2% is observed by merging the two multivariate models which characterize the best the spatial (G 0 ) and color (Laplace) dependencies (in red, 97.97%) compared to the multivariate Laplace model which represents simultaneously the color and spatial dependencies (in pink, 95.94%).

CONCLUSION

In this paper, a multiscale texture retrieval algorithm using the GD between elliptically distributed random variables has been proposed. The multivariate G 0 and Laplace distributions have been introduced to model respectively the spatial and color dependencies of the wavelet coefficients. An analytical expression of the GD for fixed shape texture parameters and an approximation of the GD when assuming the geodesic coordinate functions as straight lines have been given for the G 0 distribution. Some experiments on the VisTex database have shown that the combination of those two criteria leads to the best performances.

Fig. 1 .

 1 Fig. 1. Average classification rate as a function of the number of training samples.

Table 1 .

 1 Average retrieval rate (in %) on the VisTex database.

		Gaussian			GD
		copula Gamma	Gaussian SCM	Gaussian FP	G 0 fixed shape	G 0 linear	Laplace
	spatial	76.82	71.51	76.99	77.44	79.38	73.56
	color	85.02	86.10	86.96	87.78	86.69	87.55
	spatial & color	85.83	86.12	88.23	89.26	89.65	88.09


Benjamin Lesage

Isabelle Puaut
email: isabelle.puaut@irisa.fr

André Seznec
email: andre.seznec@inria.fr

PRETI: Partitioned REal-TIme shared cache for mixed-criticality real-time systems

Multithreaded processors raise new opportunities as well as new issues in all application domains. In the context of real-time applications, it has created one major opportunity and one major difficulty. On the one hand, the concurrent execution of multiple threads creates the opportunity to mix, on the same hardware platform, the execution of a complex real-time workload and the execution of non-critical applications. On the other hand, for real-time tasks, timing deadlines must be met and enforced. Hardware resource sharing, inherent to multithreading, hinders the timing analysis of concurrent tasks. Two different objectives are pursued in this work: enforcing timing deadlines for real-time tasks, and achieving the highest possible performance for the non-critical workload.

In this paper, we present the PRETI, Partitioned REal-TIme shared cache scheme, a flexible hardwarebased cache partitioning scheme that aims at pursuing these two objectives at the same time. Plainly considering inter-task conflicts on shared caches for real-time tasks yields very pessimistic timing estimates.

We remove this pessimism by allocating private cache space for real-time tasks. During the execution of a real-time task, our scheme reserves a fixed number of cache lines per set for the task. Therefore, uniprocessor, i.e. unithread, worst-case execution time (WCET) estimation techniques can be used, resulting in tight WCET estimates. Apart from the private spaces reserved for the real-time tasks currently running, the remaining cache space is shared by all tasks running on the processor, in particular non-critical tasks, enabling high performance for these tasks. While the PRETI cache scheme is of great help to achieving our both objectives, its hardware implementation consists only in a slight modification of the LRU cache replacement policy.

Experiments are presented to show that the PRETI cache scheme allows for both guaranteeing the schedulability of a set of real-time tasks with tight timing constraints, and enabling high performance for the non-critical tasks.

Introduction

Need for performance, as well as the economic pressure to use off-the-shelf components or standard IPs, lead to the use of standard microarchitectures in a large range of application domains. Most processors are designed to achieve the best average performance as possible on as large as possible application domains.

The use of multithreaded processors in the real-time systems domain creates the opportunity to mix on the same platform the execution of critical real-time tasks with a non-critical workload. However, using multithreaded processors to execute critical real-time tasks creates at the same time a major difficulty, since it must be guaranteed that all critical tasks meet their timing constraints in all cases, including the worst-case.

This safety property of the worst-case execution time (WCET) estimate is of particular importance for hard real-time systems, e.g. control of airplanes: one has to guarantee that the effective execution time will always be smaller than the WCET estimation. Although research on WCET estimation on uniprocessor featuring instruction level parallelism and complex memory hierarchy has made important progress in the last two decades [START_REF] Wilhelm | The worst-case execution-time problem : overview of methods and survey of tools[END_REF], research in the domain is still going to face the ever increasing complexity of hardware. In general, the more complex is the processor microarchitecture, the less precise are estimated WCETs. In the context of multithreaded processors, the major source of pessimism in WCET estimation comes from the sharing of the cache between hardware threads. While meeting deadlines is the objective for critical real-time tasks, achieving high performance is the goal for non-critical workloads. Ideally, one would like to address these two, apparently opposite, objectives with a standard microprocessor. Moreover the proposal of a new microarchitectural feature will not be considered by the industry if it only targets real-time systems unless it induces very marginal modification of the microarchitecture. In this paper, we propose a very simple modification of the replacement policy in the shared cache of a multithreaded processor that both addresses meeting strict deadlines in real-time systems and achieving high performance on non-critical workloads.

On a multithreaded processor, the threads share the hardware core resources and the memory hierarchy, in particular the first level cache. Functional units sharing can create difficulties to compute a very tight (within a few cycles) WCET estimates for each basic block on superscalar processors [START_REF] Li | Modeling out-of-order processors for WCET analysis[END_REF]. However the actual estimate is much more sensitive to the worst-case cache behavior, where a hit or a miss may result in hundreds of cycles difference on the execution time of a short basic block. Therefore, in this paper we will assume a very simple and predictable model for the execution core, whereas we will concentrate our efforts on (i) guaranteeing deadlines through worst-case analysis of the cache behavior for critical tasks and (ii) improving cache behavior for non-critical tasks.

In this paper, we present the PRETI, Partitioned REal-TIme shared cache scheme, a flexible hardwarebased cache partitioning scheme that mitigates the rigidity of classic cache partitioning schemes. In order to avoid very pessimistic timing estimates associated with inter-task conflicts on real-time critical tasks, we create a private cache space for these tasks. All along the execution of a real-time task, our scheme reserves a fixed number of cache lines per set for the task. When a task is allocated N lines per set, the PRETI scheme does not allocate a fixed subset of N ways for the task, but instead guarantees that for each set, the N most recently used blocks from the task, mapped in this set, are present in the cache, regardless of their actual assigned ways in the cache. This guarantees that the analysis of the worst-case cache behavior of the task, assuming for a N-way set-associative cache, is safe. Apart from the private spaces reserved for the real-time tasks currently running, the cache space is shared by all tasks running on the processor, i.e. non-critical tasks, but also the real-time tasks for their least recently used memory blocks. This allows for enhancing the performance of the non-critical threads and also for reducing the effective execution time of the real-time tasks. In practice, the PRETI cache scheme is implemented as a slight modification of a conventional LRU replacement policy.

In summary the PRETI cache scheme:

• allows providing predictability guarantees for critical real-time tasks, through reserved space in the shared cache for these tasks • allows as good as possible performance for non-critical tasks, through the dynamic allocation of cache lines not currently used by critical tasks,

• is a marginal modification of conventional LRU replacement, and as such could be considered by the industry as an architectural support for real-time applications In this paper, hard real-time systems, in which safety is critical, are considered to illustrate the PRETI cache. The PRETI cache scheme exhibits the same benefits in the context of soft real-time systems where the safety of WCET estimates can be relaxed.

The remainder of this paper is organized as follows. Section 2 surveys the related works on using caches in real-time systems. The PRETI cache scheme is detailed Section 3 including considering hardware implementation and WCET analysis aspects. Section 4 gives the results of experimentations conducted to illustrate the behavior of our proposal. Finally, Section 5 summarizes the present study and gives directions for future works.

Related work

Sharing caches on real-time systems has essentially been considered for multi-cores or multiprocessors, or in the context of time shared-multiprocessors. To the best of our knowledge, no specific studies has addressed sharing caches in multithreaded processors on real time caches. However, most of the research done in the context of multiprocessors or multicores is directly relevant.

For multi-core architectures, state-of-the-art shared cache analyzes [START_REF] Li | Timing analysis of concurrent programs running on shared cache multi-cores[END_REF][START_REF] Hardy | Using bypass to tighten WCET estimates for multi-core processors with shared instruction caches[END_REF], stemming from abstract interpretation based analyzes [START_REF] Ferdinand | Efficient and precise cache behavior prediction for real-timesystems[END_REF], have raised the problems of inter-task conflicts and their negative impact on tasks' predictability. The use of bypass [START_REF] Hardy | Using bypass to tighten WCET estimates for multi-core processors with shared instruction caches[END_REF], information about tasks lifetimes [START_REF] Li | Timing analysis of concurrent programs running on shared cache multi-cores[END_REF], or the conjoint use of partitioning and locking [START_REF] Suhendra | Exploring locking & partitioning for predictable shared caches on multi-cores[END_REF] have been recommended to reduce the impact of these conflicts. Similarly to these works, our proposal aims at containing inter-task conflicts, hence improving systems predictability.

Cache partitioning, in the context of real-time systems, has been mostly explored for preemptive systems, to eliminate cache related preemption delays [START_REF] Tan | A prioritized cache for multi-tasking real-time systems[END_REF][START_REF] Wolfe | Software-based cache partitioning for real-time applications[END_REF][START_REF] Mueller | Compiler support for software-based cache partitioning[END_REF]. The advent of multi-processors and multi-core architectures has led to the use of partitioning for shared caches, to preclude inter-task conflicts in such concurrent systems [START_REF] Kirk | Smart (strategic memory allocation for real-time) cache design[END_REF][START_REF] Suhendra | Exploring locking & partitioning for predictable shared caches on multi-cores[END_REF].

The partitioning schemes introduced, or assumed, in related work [START_REF] Wolfe | Software-based cache partitioning for real-time applications[END_REF][START_REF] Mueller | Compiler support for software-based cache partitioning[END_REF][START_REF] Kirk | Smart (strategic memory allocation for real-time) cache design[END_REF][START_REF] Suhendra | Exploring locking & partitioning for predictable shared caches on multi-cores[END_REF] tend to be shackled by their inherent strictness. Tasks' memory blocks are kept within the limits of their allocated partitions, and conversely, a partition in the cache can only be modified by its owner. Our proposal alleviates some of these restrictions, to maximize cache space utilization, and allows for reasonable average-case performances.

Prioritized caches [START_REF] Tan | A prioritized cache for multi-tasking real-time systems[END_REF] were proposed as a mean to advantage the highest priority tasks in the context of preemptive systems. A cache line can only be redeemed by a task of priority higher than the current owner of the line. Hence, the higher the task priority, the less it suffers from inter-task conflicts. Furthermore, resources are allocated to tasks on demand. This proposal allows to mix critical and no-critical tasks. We use a similar on-demand cache line allocation mechanism. Unlike prioritized caches, using the proposed PRETI scheme, partitions' sizes are upper-bounded to prevent the trashing of the whole cache by high priority tasks.

A software-based per-set partitioning has also been proposed for preemptive real-time systems [START_REF] Wolfe | Software-based cache partitioning for real-time applications[END_REF][START_REF] Mueller | Compiler support for software-based cache partitioning[END_REF].

Using the compiler, the memory mapping of the tasks of the system is altered to ensure their spatial isolation in the cache. However, altering tasks' mapping in the memory is far from trivial; significant modifications of the compilation tool chain are required to dispatch tasks' memory blocks to appropriate addresses in the memory, depending on the cache sets where they belong. Instead, we propose a simple hardware-based partitioning scheme. S.M.A.R.T. caching [START_REF] Kirk | Smart (strategic memory allocation for real-time) cache design[END_REF] addresses the additional issue of shared data structures. Using this hybrid hardware and software scheme, cache sets are grouped into partitions and a special, statically-sized partition, the shared pool, is used to cache data shared between multiple tasks, and bear less critical tasks. The partitioning scheme proposed in this paper uses a similar shared space to accommodate for additional tasks. However, our shared space grows and shrinks dynamically as resources are freed or allocated.

The per-set division of the cache, shared by compile-time methods [START_REF] Wolfe | Software-based cache partitioning for real-time applications[END_REF][START_REF] Mueller | Compiler support for software-based cache partitioning[END_REF] and S.M.A.R.T. caching [START_REF] Kirk | Smart (strategic memory allocation for real-time) cache design[END_REF], allows for fine-grained partitions. However, disused resources are beyond reclaim; merging set-based partitions, without flushing a part of the cache contents, could induce the duplication of memory blocks in different sets (their previous set and the post-merge one). With the PRETI scheme, we rely on a "per-way" division of the cache instead. However, this "per-way" division is virtual in the sense that the PRETI scheme guarantees that for a given task and for each set, at least the N most recently used blocks reside in the cache, but it does not impose that these blocks are stored in fixed ways.

The PRETI cache scheme has been defined for real-time systems mixing real-time tasks and non-critical tasks. It benefits from the predictability increase that can be expected from the use of partitioning for shared caches in the context of multi-core architectures [START_REF] Kirk | Smart (strategic memory allocation for real-time) cache design[END_REF], while preserving flexibility, in the sense that unclaimed cache space can be used by any task, be it real-time or not.

The PRETI cache partitioning scheme

The PRETI (Partitioned REal-TIme) shared cache is a per-way partitioning scheme intended for shared cache levels of the memory hierarchy. Unlike existing strict partitioning schemes, by per-way, we mean that, for a given task featuring a N-way allocated partition, and for each cache set, at least the N most recently used blocks are present in the cache. In the context of multithreaded architectures and real time systems, the first objective of the PRETI cache is to ensure the real-time tasks' timing predictability, through guarantees in terms of cache space available for these the critical real-time tasks. Secondly, PRETI also aims at maximizing the overall performance of the system, particularly non-critical tasks.

In the following, the behavior of the PRETI cache is first introduced ( § 3.1). Based on this description of the mechanism, possible implementations ( § 3.2), then static cache analyzes targeting the PRETI cache are discussed ( § 3.3).

Principles

The principles introduced subsequently regulate the behavior of PRETI caches, comprised of the notions which underlie the partitioning of the cache lines and the cache replacement policy.

Private and shared partitions divide the cache in distinct spaces. At creation, a task (typically a realtime task) can be allocated a private space consisting in a fixed number of ways in the set-associative cache.

Each private space is used by a single task, its owner, to store its most recently accessed memory blocks. This private space is virtual in the sense that to avoid moving blocks in the cache, the private space of N ways consists in the N most recently used blocks accessed by the task in each cache set regardless of their physical location; their effective storage place may vary during the task execution. Upon task termination, the associated private space is released.

The shared space, on the other hand, keeps blocks for all the tasks accessing the cache, whether they have a dedicated private space or not. It holds tasks' least recently accessed memory blocks. The shared space is dynamically composed from all the cache lines that do not belong to any private space: (i) cache lines occupied by blocks from tasks that do not own any private cache space, and (ii) cache lines used by a task beyond its private space allocated capacity. Figure 1 illustrates the division of a 4-way cache between shared and private spaces during the execution of two tasks, Task 1 and Task 2. Task 1 has been allocated a maximum of two cache ways, whereas Task 2 can claim the ownership of at most one way. Black (resp. grey) blocks depict blocks from the private space of Task 1 (resp. Task 2), whereas white blocks are blocks from the shared space. In Set 0, Task 1 and Task 2 fully use their private partition; the shared space of Set 0 is then reduced to one cache line, the statically unallocated one. In Set 1, Task 1 has only claimed a single line in its private space, whereas Task 2 fully uses its private space; the Set 1 shared space is composed of two cache lines, the unallocated one and one line unclaimed by Task 1 at this time instant. This example also illustrates an interesting property of the PRETI cache: real-time task are often small tasks, that sometimes do not touch every cache line before their termination. When a real-time task does not claim its overall private space, the scheme dynamically enables the other tasks (non-critical tasks or real-time tasks) to use this cache space.

The cache replacement policy is directly derived from the Least Recently Used (LRU) policy. The LRU order is maintained among all blocks present in a cache set. However, in case of a miss, the LRU block is not systematically evicted, to guarantee a minimum reserved capacity in the cache to tasks with private partitions. In case of a miss, the block to be replaced is the oldest block that:

• either belongs to the shared space,

• or belongs to the private space of the task triggering the miss The block is systematically inserted at the Most Recently Used position in the set. This replacement policy, for a task with a N-way private space, enforces the property that each set in the cache maintains at least the N most recently blocks used by the task. From the task perspective, it means that the private space acts as a N-way LRU set-associative cache. This property is of particular interest for hard real-time systems because of the timing predictability properties of LRU replacement, as demonstrated in [START_REF] Reineke | Timing predictability of cache replacement policies[END_REF].

Implementation of a PRETI cache

The hardware overhead of PRETI caches over a conventional set-associative LRU cache is only on the replacement policy. The access hit time is not modified.

On a hit, the LRU tags must be modified as on a conventional LRU cache. On a miss, the selection of the replaced block is different from the one on the LRU cache. In order to implement this selection, an extra information must be associated with the cache block, the task id of the owner task, updated on an access to the block. This is illustrated on Figure 2. The tasks having a private partition (typically real-time tasks) have distinct task ids, whereas the other tasks (typically non-critical tasks) may share the same task identifier.

Cache line data Tag

Task id The cost, in terms of storage requirements, of the proposed implementation is low: ⌈log 2 (T + 1)⌉ * ℓ bits for the task ids, where ℓ is the number of cache lines and T the maximum number of distinct private spaces which may reside in the cache at the same time. On the other hand, some extra logic is needed to perform the selection of the replaced block, but this logic is limited and is not on the critical path.

Valid bits

Analyzing PRETI cache behavior on a real-time task

As the PRETI cache is designed to provide guarantees for real-time tasks, static analysis will be used to compute WCET estimates for such tasks. WCET estimation requires to identify tasks' worst-case behavior with regard to the memory hierarchy.

In the context of a PRETI cache, the partitioning scheme enforces strong properties that enable static worst-case analysis of real-time tasks. First, by construction, private spaces, used to store tasks' most recently accessed memory blocks, can only be altered by their respective owner; they are similar to small private caches in the sense they are free of inter-task conflicts. Furthermore, each task's private space is statically known to be the very least amount of cache space available to this task during execution, and replacement policy-wise, private spaces mimic LRU-enforcing private caches.

The analysis of a task's behavior can thus be safely reduced to the analysis of said task's private space.

Such a private space analysis can be performed using state of the art private cache contents analyzes [START_REF] Theiling | Fast and precise WCET prediction by separated cache and path analyses[END_REF][START_REF] Hardy | WCET analysis of multi-level non-inclusive set-associative instruction caches[END_REF][START_REF] Lesage | WCET analysis of multi-level set-associative data caches[END_REF], while assuming a LRU cache of the same dimensions than the task's statically bounded private space.

With respect to WCET analysis, PRETI and strictly-partitioned caches behave the same way.

Nevertheless, it should be noted that sharing data among tasks causes difficulties in the worst-case analysis of the cache content. Experiments conducted by previous studies on WCET estimation for multi-core architectures with shared caches [START_REF] Li | Timing analysis of concurrent programs running on shared cache multi-cores[END_REF][START_REF] Hardy | Using bypass to tighten WCET estimates for multi-core processors with shared instruction caches[END_REF] show that data sharing yield to pessimistic WCET estimates; this is because data sharing may result in both destructive and constructive effects on the shared cache, both of them being hard to predict statically. The same difficulty applies when using the PRETI cache. A possible solution in the context of PRETI cache would be to avoid tagging shared data with the private tag of the task.

Experimental results

In this section, after a presentation of the experimental setup ( § 4.1), we illustrate the ability of the PRETI partitioning scheme, for mixed-criticality workloads, to both guarantee the schedulability of real-time task sets ( § 4.2), and to enable high performances for the non-critical tasks ( § 4.3).

Experimental setup

Simulated hardware platform For the purpose of our experiments, we will consider a very simple model of time multithreaded processor. 3 hardware threads are assumed and each thread is granted the access to the hardware every 3 cycles in a round-robin fashion.

The processor executes an instruction per cycle and the execution of a thread is stalled on any cache miss until the miss is resolved. The processor features a 8-way set-associative 4KB data cache and a same-sized 8-way set-associative instruction cache. 32 bytes cache blocks are used and LRU replacement policy is assumed for both caches. Upon a miss, the miss penalty of 150 cycles is observed by the thread, but the memory can service up to 6 misses concurrently, one data miss and one instruction miss per thread. The caches sizes correspond to the use of a small embedded processor, suited to the assumed workloads.

Software environment Each task, be it time-critical or not, is statically allocated to a single hardware thread.

Then, on each hardware thread, a non-preemptive Earliest Deadline First scheduling (NP-EDF [START_REF] Jeffay | On non-preemptive scheduling of period and sporadic tasks[END_REF]) is used for time-critical tasks. Upon each release of a time-critical task, tasks are sorted by increasing deadlines.

If no real-time task is currently running on the hardware thread, the released task with the earliest deadline is executed first, without preemption, until completion. The non-preemptive scheduling for critical tasks eliminate cache-related preemption delays from impacting the performances of time-critical tasks and allows for static WCET estimation.

Non-critical tasks execute in the background; they are allowed to execute on their assigned hardware thread only to when no time-critical task is ready, and are preempted as soon as a time-critical task is released.

Cache analysis and WCET estimation

The presented experiments have been conducted on MIPS R2000/R3000 binary code compiled with gcc 4.1 without any code optimization and using the default linker memory lay-out. WCET estimates of real-time tasks were computed using a state-of-the-art static WCET estimation tool [?] 1 using the standard Implicit Path Enumeration Technique (IPET) technique for WCET computation [START_REF] Wilhelm | The worst-case execution-time problem : overview of methods and survey of tools[END_REF]. Memory references were classified using the most recent static cache analysis techniques presented in [START_REF] Hardy | WCET analysis of multi-level non-inclusive set-associative instruction caches[END_REF][START_REF] Chattopadhyay | Unified cache modeling for WCET analysis and layout optimizations[END_REF]; every reference is classified as hit if the block is for sure in the cache, and as a miss otherwise. Focusing on caches WCET estimates only take their contribution into account.

Timing simulations

Simulations have been performed using trace driven simulations. A simple timing model is used: an instruction is fetched and executed per cycle, but the execution is stalled upon any cache miss. For the sake of simplicity, we assume that the execution overhead of the task scheduler is negligible.

Moreover, a unique memory access trace per benchmark is used.

Benchmarks Unfortunately, very few realistic real-time benchmarks are available for the research community; only benchmarks with limited execution time and data and instruction footprints are publicly available. Consequently, our experiments were conducted using the benchmark sets that are used in the context of research on hard real-time systems, in particular for the context of WCET estimation [START_REF] Gustafsson | The mälardalen WCET benchmarks -past, present and future[END_REF].

The experiments were conducted using two small but standalone real-time applications, named Debie and Papabench, as well as benchmarks from the WCET benchmark suite maintained by the Mälardalen WCET research group. The Debie software [START_REF] Holsti | Using a worst-case execution time tool for real-time verification of the debie software[END_REF], developed by Space Systems Finland, monitors the impact of space debris and micro-meteoroids using electrical and mechanical sensors. The Papabench tasks [START_REF] Nemer | PapaBench : A Free Real-Time Benchmark[END_REF] come from the Paparazzi project [1], used to build autonomous unmanned aerial vehicles.

In our experiments, the tasks from the Debie and Papabench real-time benchmarks are used as the realtime tasks of the studied workloads. The real-time tasks are released periodically, with a deadline equal to the task period. The non-critical tasks were selected among the ones maintained by the Mälardalen WCET research group. Non-critical tasks are not periodic and use all slack time left by real-time tasks, if any. The characteristics of each task are exposed in Table 1 (from left to right are the sizes in bytes of code, data, bss -uninitialized data and maximum stack sections, and execution frequency for the periodic, real-time tasks).

Task sets Three different sets of periodic time-critical tasks have been composed from the tasks presented in

Impact of cache partitioning on schedulability

In this section, we illustrate the benefits of using cache partitioning to allow for the schedulability of real-time tasks with tight timing constraints. For each real-time workload, we determine the minimum functioning frequency of the system for which the real-time workload would be safely schedulable. By safely schedulable, we mean that there exists a cache partitioning, a task allocation and a schedule such that even in the unrealistic case where the execution time of every task is equal to the WCET estimate computed by the static analysis, all the timing constraints are met.

To compute this minimum frequency, we determine the best pair (cache partitioning, task allocation)

among the threads through a systematic exploration of the possible (cache partitioning, task allocation) pairs. Assuming a non-preemptive system (i.e., no real-time task can be interrupted), when the WCETs and the timing constraints of a set of tasks are known, it has been showed in [START_REF] Jeffay | On non-preemptive scheduling of period and sporadic tasks[END_REF] that the system can be scheduled on a uniprocessor if 1) the WCET of any task is shorter than the delay between two consecutive deadlines and if 2) when a task is launched, it can complete before the deadline of any other task assigned to the same thread. This is a sufficient and necessary condition.

We determine this best pair (cache partitioning, task allocation) for the PRETI cache as well as the best task allocation for the shared cache. These best configurations are illustrated on Table 3. We also illustrate for each hardware thread the predicted number of instructions that are executed per second on its assigned real-time subset task subset. As mentioned in the introduction, the contribution of the memory hierarchy to the WCET is very large (tens of cycle per instruction) and plainly justify our simple processor model which essentially ignores the contributions of instructions dependencies and pipeline hazards. It is noticeable that assuming a fully shared cache leads to a minimum frequency for guaranteeing schedulability much higher respectively (13x, 5.8x and 5.6x ) than the minimum frequency needed for the PRETI cache.

Note that the best pair (cache partitioning, task allocation) does not always allocate the 8 cache ways among the three hardware threads. This is due to the WCET analysis that sometimes does not lead to better WCET estimates when allocating more ways to a task. In the context of PRETI cache, allocating a smaller number of ways in the private cache space allows increasing the shared space and therefore may enable higher performance on the non-critical tasks.

Performances of non-critical tasks using PRETI

In this section, we illustrate the ability of the PRETI partitioning scheme, for mixed-criticality workloads, to enable high performances for the non-critical tasks while guaranteeing schedulability.

In order to illustrate this property we run simulations assuming 3 different cache configurations, the PRETI cache, a shared cache and a completely partitioned cache. We run 3 different sets of simulations on which we fix the frequency of the processor, respectively to 1, 2 and 3 times the minimum frequency necessary to schedule the real-time workload on the PRETI cache. The system is run for 1.5 billion cycles for papabench.auto and papabench.manual and for 1 second of execution i.e one complete execution period for debie (i.e. respectively 1.5 billion, 3 billion and 4.5 billions cycles).

1. For the PRETI cache, we choose the best pair (cache partitioning, task allocation) obtained above for each workload, i.e. the pair that requires the minimum frequency to get the workload schedulable.

2. The shared cache should not be used in a hard real-time system , since there is no safety guarantee for the frequency that we will simulate. This shared cache where is presented in order to illustrate the performance that could be reached on non-critical tasks if there were no safety constraint on the real-time tasks.

3. By completely partitioned cache, we denote a cache where each hardware thread will have a fixed set partition during the complete execution, that is a hardware thread cannot access the blocks from the other hardware threads. With a completely partitioned cache, assuming the same task allocation as for the PRETI cache, the safety constraints on the real-time tasks will be respected if the partition sizes are larger or equal that the ones for the PRETI cache 2 . In our experiments, we allocate all the cache ways among the hardware threads(see Table 4). As a first round of simulations, on each hardware thread, when not running one critical task, we execute repetitively a non-critical task. The respective non-critical tasks executed are lms, adpcm and fft on the three hardware contexts. We measured the performance of each non-critical task in Instruction Per Cycle (IPC) over the whole execution time. The maximum achievable performance for the thread would be 1 instruction every 3 cycles if there were no cache miss and no interruption by the critical tasks.

The simulation results are provided in Figure 3 and illustrate the benefits of sharing caches for the noncritical tasks. The results are provided for the debie, papabench.auto and papabench.manual, from top to bottom. For each workload, from left to right, the frequency of the processor was fixed to respectively 1, 2 and 3 times the minimum frequency necessary to schedule the real-time workload on the PRETI cache. The cumulated IPC is the sum of the IPC of the non-critical tasks running on each hardware thread.

Clearly always partitioning the cache is a major handicap for Thread 2 and on a less extent for Thread 1on the debie workload: the ways of the caches are mainly reserved to guarantee the safety on Thread 0. On debie, the PRETI cache is able to share the cache space for a large portion of the time (when the Thread 1 is not running the real-time workload). On Thread 1, for the three simulations, the execution time of the realtime tasks is in the range of 350 million cycles for the PRETI cache per second. That is the overall cache is completely shared for respectively for 77 %, 88 % and 92 % of the time for the three simulated frequencies.

Therefore the performance on the non-critical tasks is in the same range when using the PRETI cache and the shared cache. The same phenomenon also occurs on the two papabench workloads. On papabench.manual, the performance on the non-critical workloads are very close for the PRETI cache and the shared cache for all threads while performance on Thread1 and Thread 2 is much lower for the partitioned cache.

Notice that the lms task, running on Thread 0, has a particular behavior. It performs slightly worse when using shared or PRETI caches, than when using the partitioned cache. Its blocks tend to be evicted from the shared space by the other tasks because it exhibit lower temporal locality. On the other hand, the performance of the other tasks improve when sharing the cache space (shared and PRETI cache).

As a second round of simulations, we only execute one non-critical task concurrently with the real-time workload. This non-critical task is scheduled on the first hardware thread that becomes available. We run 3 sets of simulations using respectively lms, adpcm and fft as this single non-critical task. Simulation results in Figure 4 globally illustrates that the PRETI cache structure allows the non-critical task to benefit from shared cache space. The IPC of the non-critical task is presented by processor frequency, from top to bottom, then by concurrent real-time workload and non-critical task.

Completely partitioning the cache results in a worse cache behavior than using a fully shared cache, especially at higher frequencies, 2x and 3x. As stated earlier, using PRETI caches, the overall cache is completely shared for more than 77% of the time for each simulation. Therefore, as frequencies increases, the gap between the performance of the non-critical task using completely partitioned caches, and their performance using shared caches, increase. Conversely, the gap between the performance of PRETI caches and shared caches shrinks.

As when running 3 threads, the sensitivity of the lms benchmark to sharing the cache space with other tasks induces some particular behavior: on papabench.auto and for frequency 1x, Thread 0 is busy for 64 % of the time with the real-time tasks, this induces a slight performance degradation for lms for shared cache and PRETI cache compared with the partitioned cache. For higher frequencies, lms runs alone for a longer time and the use of shared or PRETI cache on this interval compensates for the performance loss encountered when sharing.

Conclusion

During the past decade, the need for performance and the tremendous integration technology progress have allowed inexpensive superscalar processors featuring instruction level parallelism as well as memory hierarchy. Moreover, nowadays, these processors targeting large scale volumes are also featuring hardware thread level parallelism. Economic pressure is pushing for the use of these standard processor components or IPs in every application domain. As a consequence, an application domain can only marginally influence the design of the standard processors. That is the processor manufacturers will consider adding a feature to a standard processor only if the application domain generates huge volume, or if the modification is very marginal and does not impair the rest of the design. This applies for real-time systems.

Standard processors are optimized for average case performance. Unfortunately, real-time systems must be designed to enforce timing constraints in all situations, including the worst-case scenarios. The increase of performance of standard processors has been coupled with an explosion of the complexity (instruction level parallelism, memory hierarchy and thread level parallelism). This complexity has been widening the gap between the effective average run-time execution time of real-time tasks and the WCET estimate that can be safely obtained through static analysis. Every increase in hardware complexity (memory hierarchy, shared cache levels) has resulted in a larger gap between actual performance and predicted worst-case performance [START_REF] Hardy | WCET analysis of multi-level non-inclusive set-associative instruction caches[END_REF]. This has lead some authors to propose to completely avoid cache sharing, and instead rely on only private caches or equivalently strictly partitioning schemes [START_REF] Paolieri | Ia3: An interference aware allocation algorithm for multicore hard real-time systems[END_REF]. On the other hand, the use of standard processors in real time systems has raised the possibility to execute non-critical tasks concurrently with critical tasks on the same platform. These non-critical tasks are not concerned with meeting deadlines but with maximizing their performance: for these non-critical tasks cache sharing is helping to ensure high performance.

In this paper, we have proposed a small adaptation of the cache replacement policy on a shared cache for a multithreaded processor, the PRETI cache scheme. The PRETI cache scheme is a small variation of a shared LRU cache scheme. But this small modification can be of major benefit for using the multithreaded processor in a real-time system. By allowing critical real-time tasks to grab a private cache space, the PRETI cache allows WCET analysis to produce WCET estimates that are tighter than if a shared cache was used. Through this slight modification, we become able to guarantee the schedulability of some real-time workloads that would not be schedulable using a shared cache. By releasing this private space to all the tasks upon the critical task termination, the PRETI cache enables high performance on all the tasks, i.e.

performance closer to the one obtained on a completely shared LRU cache.

While our study and experiments have been done in the context of hard real-time systems, where missing a constraint can be safety critical (e.g. controlling a car airbag), the PRETI cache scheme could also be used in soft real-time systems where missing a constraint is less critical (e.g. missing a frame in a video).

As for hard real-time systems, the PRETI cache scheme can provide the soft real-time task with a private cache space that isolates its working set from interferences from other tasks, thus providing the same cache behavior predictability as a private cache. It also allows the tasks to benefit from the overall cache space when the real-time task is finished.

Figure 1 :

 1 Figure 1: Division of a 4-way cache between shared and private spaces during execution.

Figure 2 :

 2 Figure 2: Storage requirements for PRETI caches: each cache line is extended with a task id.

  (a) Performances of the non-critical tasks alongside the debie workload. (b) Performances of the non-critical tasks alongside the papabench.auto workload.(c) Performances of the non-critical tasks alongside the papabench.manual workload.

Figure 3 :

 3 Figure 3: IPC of the non-critical tasks, allocated to different threads, with the real-time workloads.

  (a) Performances of non-critical tasks under different configurations, at frequency 1x. (b) Performances of non-critical tasks under different configurations, at frequency 2x. (c) Performances of non-critical tasks under different configurations, at frequency 3x.

Figure 4 :

 4 Figure 4: IPC of each non-critical tasks executed as the sole non-critical task, concurrently with a real-time workload, per processor frequency (1x, 2x or 3x).

Table 1 .

 1 The resulting task sets are given in Table2. The debie task set is composed of all tasks of 1 Reference omitted for blind review

	Task	Code	Data	Bss	Stack Frequency
		(bytes) (bytes) (bytes) (bytes)
	Debie
	acquisition task	5968	1280 101904	200	100Hz
	hit trigger handler	2408	647	218	104	400Hz
	monitoring task	12372	667	65772	296	1Hz
	tc execution task	12928	1280 101904	200	100Hz
	tc interrupt handler	2800	134	35648	64	1000Hz
	tm interrupt handler	872	23	168	56	1000Hz
	Papabench
	task 1	3960	74	78	32	40Hz
	task 2	696	52	242	32	40Hz
	task 3	1492	64	66	32	20Hz
	task 4	372	4532	568	32	20Hz
	task 5	1440	4257	20	32	20Hz
	task 6	8508	572	302	56	40Hz
	task 7	1352	194	158	32	20Hz
	task 8	264	212	46	8	20Hz
	task 9	10284	552	511	88	4Hz
	task 10	33700	8192	302	1160	4Hz
	task 11	268	90	0	32	4Hz
	task 12	1120	145	0	32	4Hz
	task 13	21176	284	255	0	10Hz
	Mälardalen
	adpcm	7872	1300	436	152	-
	fft	3772	88	128	296	-
	lms	2992	248	1092	184	-

Table 1 :

 1 Benchmark characteristics

	Task set	Tasks
	debie	acquisition task,	hit trigger handler,	monitoring task,
		tc execution task, tc interrupt handler, tm interrupt handler
	papabench.manual tasks 1, 2, 3, 4, 5, 6, 7, 8, and 13
	papabench.auto	tasks 3, 4, 5, 7, 8, 9, 10, 11, 12, and 13

the Debie benchmark. The papabench.manual and papabench.auto task sets correspond respectively to a remote-controlled and an automatic operating mode of a drone.

Table 2 :

 2 Task sets

Table 3 :

 3

Best configurations for the debie, papabench.auto and papabench.manual under the best partitioning or using shared caches.

Table 4 :

 4 Cache partitions allocated to the hardware threads under completely partitioned caches, for the debie, papabench.auto and papabench.manual task sets.

	Thread 0 Thread 1 Thread 2

As mentioned in the previous section, on the PRETI cache, we do not need to allocate all the cache ways to private spaces, as long as the timing constraints are respected