Mohamed Riadh Trad

Alexis Joly
email: alexis.joly@inria.fr

Nozha Boujemaa
email: nozha.boujemaa@inria.fr

Large scale visual-based event matching

Organizing media according to real-life events is attracting interest in the multimedia community. Event-centric indexing approaches are very promising for discovering more complex relationships between data. In this paper we introduce a new visual-based method for retrieving events in photo collections, typically in the context of User Generated Contents. Given a query event record, represented by a set of photos, our method aims to retrieve other records of the same event, typically generated by distinct users. Similarly to what is done in state-of-the-art object retrieval systems, we propose a two-stage strategy combining an efficient visual indexing model with a spatiotemporal verification re-ranking stage to improve query performance. For efficiency and scalability concerns, we implemented the proposed method according to the MapReduce programming model using Multi-Probe Locality Sensitive Hashing. Experiments were conducted on LastFM-Flickr dataset for distinct scenarios, including event retrieval, automatic annotation and tags suggestion. As one result, our method is able to suggest the correct event tag over 5 suggestions with a 72% success rate.

INTRODUCTION

An event can be described as an action that occurs at a specific time in a specific place. This notion is potentially useful for connecting individual facts and discovering complex relationships. It is worth noting that photos in User Generated Content (UGC) websites, as well as in personal collections, are often organized into events. Indeed, users are usually more likely to upload or gather pictures related to the same event, such as a given holiday trip, a music concert, a wedding, etc. This applies to professional contents such as journalism or historical data that are even more systematically organized according to hierarchies of events. Defining new methods for organizing, searching and browsing media according to real-life events is therefore gaining interest in the multimedia community [START_REF] Troncy | Linking events with media[END_REF][START_REF] Gkalelis | A joint content-event model for event-centric multimedia indexing[END_REF]. In this paper, we address the problem of matching distinct records of the same event in picture datasets, typically in UGC's photo collections. Given a query event record represented by a set of photos, our method aims to retrieve other records of the same event, notably those generated by other actors or witnesses of the same real-world event. An illustration of two matching event records is presented in Figure 1. It shows how a small subset of visually similar and temporally coherent pictures might be used to match the two records, even if they include other distinct pictures covering different aspects of the event. Application scenarios related to such a retrieval paradigm are numerous. By simply uploading their own record of an event users might, for example, access to the community of other participants. They can then revive the event by browsing or collecting new data complementary to their own view of the event. If some previous event's records had already been uploaded and annotated, the system might also automatically annotate a new record or suggest some relevant tags. The proposed method might also have nice applications in the context of citizen journalism. Automatically detecting the fact that a large number of amateur users did indeed record data about the same event would be very helpful for professional journalists in order to cover breaking news. Finally, tracking events across different media has a big potential for historians, sociologists, politicians, etc. Of course, in such scenarios, time and geographic information provided with the contents have a major role to play. Our claim is that using visual content as complementary information might overcome several limitations of approaches that rely only on metadata. First of all, distinct records of the same event are not necessarily located at the same place or can be recorded at different times. Some events might, for example, have wide spatial and temporal coverage such as a volcano eruption or an eclipse, so that geocoordinates and time stamps might be not sufficiently discriminant. This lack of discrimination can be problematic even for precisely located events, typically in crowded environments such as train stations, malls or tourist locations. In such environments, many records might be produced at the same time and place while being related to very distinct realworld events. Furthermore, in a wider meaning of the event concept, several instances of an event might be recorded at different times, e.g. periodical events or events such as "a trip to Egypt" illustrated in Figure 2. Finally, location and time information is not always available or might be noisy. The Flickr dataset used in the experiments reported in this paper notably does not contain any geographic information and contains noisy time information (as discussed in section

5).

The retrieval method we propose makes use of both visual content and contextual meta-data. Visual content is used in the first stage to detect potential matches whereas geotemporal metadata are used in the second stage to re-rank the results and therefore estimate the spatiotemporal offset between records. It is important to notice that our method allows spatiotemporally coherent records to be retrieved even if they were not produced at the same time and place (such as the examples discussed above).

RELATED WORKS

To the best of our knowledge, this work is the first to deal with visual based event matching in picture collections. Some recent works are more generally concerned with event models and ontology-based event retrieval, such as [START_REF] Troncy | Linking events with media[END_REF] or [START_REF] Gkalelis | A joint content-event model for event-centric multimedia indexing[END_REF]. In the latter, a joint content-event model is proposed to facilitate the automatic enrichment of event elements with information extracted by automatic analysis of multimedia content segments. The automatic analysis of an event's visual contents is not however addressed in this work. In a survey paper on event mining in multimedia streams, Xie et al. [START_REF] Xie | Event Mining in Multimedia Streams[END_REF] give more insights on how visual features might be used to contribute to multi-cue events detection. But here again, they just mention some high level principles and do not experiment any realization. The only work we found that really uses and experiments the use of visual content in the context of event based indexing is that by [START_REF] Quack | World-scale mining of objects and events from community photo collections[END_REF]. In this work the authors present a multimedia mining framework allowing the discovery of picture clusters from multiple cues including visual content, text content and metadata. The produced clusters can then be classified in either object or event type and further annotated by linking their content to wikipedia pages. Retrieving different instances of an event is thus not really addressed.

We also mention that a large number of studies focused on detecting or recognizing events in videos [START_REF] Schuldt | Recognizing human actions: A local svm approach[END_REF][START_REF] Xu | Live sports event detection based on broadcast video and web-casting text[END_REF][START_REF] Liu | Recognizing realistic actions from videos[END_REF][START_REF] Young-Kee | Content based event retrieval using semantic scene interpretation for automated traffic surveillances[END_REF] notably human actions [START_REF] Schuldt | Recognizing human actions: A local svm approach[END_REF], sports events [START_REF] Xu | Live sports event detection based on broadcast video and web-casting text[END_REF] or video-surveillance events [START_REF] Liu | Recognizing realistic actions from videos[END_REF][START_REF] Young-Kee | Content based event retrieval using semantic scene interpretation for automated traffic surveillances[END_REF]. It might be interesting to build upon some aspects of these methods in the context of imagebased event records but they are clearly not directly applicable to our context. Most of them involve video specific algorithms such as tracking, space-time visual features, etc.

Finally, our work is, to some extent, related to object retrieval in picture collections. Our method is indeed very similar to state-of-the-art large-scale object retrieval methods combining efficient bag-of-words or indexing models with a spatial verification re-ranking stage to improve query performance [START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF][START_REF] Joly | A Posteriori Multi-Probe Locality Sensitive Hashing[END_REF]. We might give the following analogy: images are replaced by event records (picture sets), local visual features are replaced by global visual features describing each picture of a record, spatial positions of the local features are replaced by the geo-coordinates and time stamps of the pictures. Matching spatially and temporally coherent event records is finally equivalent to retrieving geometrically consistent visual objects.

VISUAL BASED EVENT MATCHING

We first describe the proposed method in the general context of event records composed of a set of geo-tagged and time coded pictures. We further restrict ourselves to time coded only pictures since our experimental dataset did not include geo-tags. We consider a set of N event records Ei, each record being composed of Ni pictures I i j captured from the same real-world event. Each picture is associated with a geocoordinate x i j and a time stamp t i j resulting in a final geotemporal coordinate vector P i j = (x i j , t i j ). The visual content of each image I i j is described by a visual feature vector F i j ∈ R d associated with a metric d : R d × R d → R. Now let Eq be a query event record represented by Nq pictures, with associated visual features F q j and geo-temporal metadata P q j . Our retrieval method works as follows:

STEP 1 -Visual Matching: Each query image feature F q j is matched to the full features dataset thanks to an efficient similarity search technique (see section 4). It typically returns the approximate K-nearest neighbors according to the used metric d (i.e the K most similar pictures). When multiple matches occur for a given query image feature and a given retrieved record, we only keep the best match according to the feature distance. The visual matching step finally returns a set of candidate event records Ei, each being associated with M q i picture matches of the form (I q m , I i m ).

STEP 2 -Stop List:

Only the retrieved records with at least two image matches are kept for the next step, i.e {Ei | M q i ≥ 2} 1≤i≤N STEP 3 -Geo-temporal consistency: For each remaining record, we compute a geo-temporal consistency score by estimating a translation model between the query record and the retrieved ones. The resulting scores Sq(Ei) are used to produce the final records ranking returned for query Eq. The translation model estimation is based on a robust re-

Δ(Eq, Ei) = arg min Δ M q i X m=1 ρ θ " P q m -(P i m + Δ) " (1)
where P q m and P i m are the geo-temporal coordinates of the m-th match (I q m , I i m ). The cost function ρ θ is typically a robust M -estimator allowing outliers to be rejected with a tolerance θ (in our experiments we used Tukey's robust estimator). The estimated translation parameter Δ should be understood as the spatial and temporal offset required to register the query event record Eq with the retrieved event record Ei. Once this parameter has been estimated, the final score of an event Ei is finally computed by counting the number of inliers, i.e the number of visual matches that respect the estimated translation model:

Sq(Ei) = M q i X m=1 "' ' 'P q m -(P i m + b Δ) ' ' ' ≤ θ " ( 2
)
where θ is a tolerance error parameter, typically the same as the one used during the estimation phase. In practice, we use a smooth counting operator to get a better dynamic on resulting scores. When we restrict ourselves to temporal metadata (as was done in the experiments), equation 1 can be simplified to:

δ(Eq, Ei) = arg min δ M q i X m=1 ρ θ " t q m -(t i m + δ) " (3)
where δ represents the estimated temporal offset between Eq and Ei and θ is now a temporal tolerance error whose value is discussed in the experiments. Since δ is a single mono-dimensional parameter to be estimated, equation 3 can be resolved efficiently by a brut force approach testing all possible solutions δ.

Final scores then become:

Sq(Ei) = M q i X m=1 "˛t q m -(t i m + δ) ˛≤ θ " (4)
STEP 4 -Prior constraints: Depending on the application context, strong effectiveness improvements might be obtained by adding prior constraints on the tolerated values for b Δ. Rejecting events with too large a spatial and/or temporal offset from the query record is indeed a good way to reduce the probability of false alarms. In our experiments we study the impact of such a constraint on the estimated temporal offsets. Concretely, we reject from the result list all retrieved event records having an estimated offset above a given threshold δmax (regardless the matching score Sq(Ei)).

ENABLING SCALABILITY

To allow fast visual matching in large picture datasets, we implemented a distributed similarity search framework based on Multi-Probe Locality Sensitive Hashing [START_REF] Lv | Multi-probe lsh: Efficient indexing for high-dimensional similarity search[END_REF][START_REF] Joly | A Posteriori Multi-Probe Locality Sensitive Hashing[END_REF] and the MapReduce [START_REF] Dean | Mapreduce: simplified data processing on large clusters[END_REF] programming model.

Multi-Probe LSH

To process Nearest Neighbors search efficiently, we use an approximate similarity search structure, namely Multi-Probe Locality Sensitive Hashing (MP-LSH) [START_REF] Lv | Multi-probe lsh: Efficient indexing for high-dimensional similarity search[END_REF][START_REF] Joly | A Posteriori Multi-Probe Locality Sensitive Hashing[END_REF]. MP-LSH methods are built on the well-known LSH technique [START_REF] Datar | Locality-sensitive hashing scheme based on p-stable distributions[END_REF], but they intelligently probe multiple buckets that are likely to contain results. Such techniques have been proved to overcome the over-linear space cost drawback of common LSH while preserving a similar sub-linear time cost (with complexity O(N λ )). Now, let F be the dataset of all visual features F ∈ R d (i.e. the one extracted from the pictures of the N event records Ei). Each feature F is hashed with a hash function g : R d → Z k such that:

g(F) = (h1(F), ..., h k (F)) (5)
where individual hash functions hj are drawn from a given locality sensitive hashing function family. In this work we used the following binary hash function family which is known to be sensitive to the inner product:

h(F) = sgn (W.F) ( 6
)
where W is a random variable distributed according to N (0, I). The produced hash codes g i = g(Fi) are thus binary hash codes of size k. At indexing time, each feature Fi is mapped into a single hash table T according to its hash code value g i . As a result, we obtain a hash table of N b buckets where

N b ≤ 2 k .
At query time, the query vector Fq is also mapped onto the hash table T according to its hash code value g q . The multi-probe algorithm then selects a set of Np buckets {(bj)}j=1..N p as candidates that may contain objects similar to the query according to : d h (g q , bj) < δMP [START_REF] Joly | A Posteriori Multi-Probe Locality Sensitive Hashing[END_REF] where d h is the hamming distance between two binary hash codes and δMP is the multi-probe parameter (i.e. a radius of hamming space). A final step is then performed to filter the features contained in the selected buckets by computing their distance to the query and keeping the K Nearest Neighbors.

The MapReduce framework

MapReduce is a programming model introduced by Google to support distributed batch processing on large data sets. A MapReduce job splits the input dataset into independent chunks which are processed by the map tasks in a parallel manner. The framework sorts the outputs of the maps, which are then input to the reduce tasks. Chunks are processed based on key/value pairs. The map function computes a set of intermediate key/value pairs and, for each intermediate key, the reduce function iterates through the values that are associated with that key and outputs 0 or more values. The map and Reduce tasks scheduling is performed by the framework. In a distributed configuration, the framework assigns jobs to the nodes as slots become available. The number of map and reduce slots as well as chunk size can be specified for each job, depending on the cluster size. With such a granularity, large data sets processing can be distributed efficiently on commodity clusters.

Multi-Probe LSH in the MapReduce framework

The hash table T in the MapReduce framework is stored in a text file where each line corresponds to one single bucket. Each bucket is represented by a < key, value > pair: < b, ((id(F1), F1), (id(F2), F2), . . . ) > [START_REF] Liu | Recognizing realistic actions from videos[END_REF] where b is the hash code of the bucket. In order to be processed by the MapReduce framework, the table T has to be divided into a set of splits. The number of splits is deduced by the MapReduce framework according to a set of input parameters as the number of available slots and the minimal input split size which is related to the file system block size. However, in order to be entirely processed by a mapper, a bucket cannot spill over different splits.

Since MapReduce is mainly dedicated to batch processing, setting up tasks could be expensive due to process creation and data transfer. Therefore, our implementation processes multiple queries at a time, typically sets of pictures belonging to the same records.

The hash codes of all query features are computed and passed to the map instances to be executed on the different slots. The number of map instances is computed by the MapReduce framework according to the number of input splits.

Each map process iterates over its assigned input split and for each query selects the candidate buckets that are likely to contain similar features according to Equ.7. It then computes the distance to each feature within the selected buckets. For each visited feature Fi, the mapfunction outputs a < key, value > pair of the form: < id(Fq), (dist(Fq, Fi), id(Fi)) > (9) where id(F) denotes the picture identifier associated to feature F and dist(Fq, Fi) the distance between Fq and Fi.

For each query identifier id(Fq) the reduce instance sorts the set of emitted values for all map instances and filters the K-nearest neighbors. Figure 3 gives the average response time per K-NN search according to the total number of queries batched within the same MapReduce job. It shows that the MapReduce framework becomes profitable from about 50 grouped queries. The average response time becomes almost constant for more than 400 grouped queries. In our experiments, the number of images per event record range from about 5 to 200. That means that using the MapReduce framework is still reasonable for the online processing of a single event record. Finally, many MapReduce implementations materialize the entire output of each map before it can be consumed by the reducer in order to ensure that all maps successfully completed their tasks. In [START_REF] Condie | MapReduce Online[END_REF], Condel et al. propose a modified MapReduce architecture that allows data to be pipelined between operators. This extends the MapReduce programming model beyond batch processing, and can reduce completion times while improving system utilization for batch jobs as well.

EXPERIMENTS

We evaluated our method on a Flickr image dataset using last.fm tags as real-world events ground truth. It was constructed from the corpus introduced by Troncy et al. [START_REF] Troncy | Linking events with media[END_REF] for the general evaluation of event-centric indexing approaches. This corpus mainly contains events and media descriptions and was originally created from three large public event directories (last.fm, eventful and upcoming). In our case, we only used it to define a set of Flickr images labeled with last.fm tags, i.e. unique identifiers of music events such as concerts, festivals, etc. The images themselves were not provided in the data and had to be crawled resulting in some missing images. Unfortunately, in this corpus, only a small fraction had geo-tags so that we evaluated our method using only temporal metadata. We used the EXIF creation date field of the pictures to generate the time metadata used in our method. Only about 50% of the crawled images had such a valid EXIF (others had empty or null date fields). In Table 1, we report the statistics on the original, crawled and filtered dataset. To gather the pictures in relevant event records, we used both the last.fm identifier and the Flickr author field provided with each picture. An event record is then defined as the set of pictures by a given author having the same LastFM label. Our final dataset contains 41, 294 event records related to 34, 034 distinct LastFM events.

Experimental settings

We used 6 global visual features to describe a picture's visual content (including HSV Histogram [START_REF] Ferecatu | Image retrieval with active relevance feedback using both visual and keyword-based descriptors[END_REF], Hough histogram [START_REF] Ferecatu | Image retrieval with active relevance feedback using both visual and keyword-based descriptors[END_REF], Fourier histogram [START_REF] Ferecatu | Image retrieval with active relevance feedback using both visual and keyword-based descriptors[END_REF], edge orientation histogram [START_REF] Ferecatu | Image retrieval with active relevance feedback using both visual and keyword-based descriptors[END_REF]). Each feature was L2-normalized and hashed into a 1024 bits hash code using the same hash function as the one used to construct the hash table (see Equ.6). The 6 hash codes were then concatenated into a single hash code of 6144 bits. We used the Hamming distance on these hash code as visual similarity.

From the full set of 41, 294 event records in the dataset, the only queries we kept, were the records being tagged with last.fm events and having at least 7 records in the dataset. We finally get 172 query records Eq. This procedure was motivated by the fact that a very large fraction of events were represented by only one record and therefore not usable for experiments.

In all experiments, we used a leave-one-out evaluation procedure and measured performances with 2 evaluation metrics: Mean Average Precision (MAP) and Classification Rate (CR). MAP is used in most information retrieval evaluations and measures the ability of our method to retrieve all the records related to the same event as the query one. Classification rate is obtained by using our method as a nearest neighbors classifiers. The number of occurrences of retrieved events is computed from the top 10 returned records and we keep the event with the maximum score as the best prediction. It measures the ability of our method to automatically label some unknown query event record. We extend this measure to the case of multiple labels suggestion. In addition to the best retrieved event we also return the following events by decreasing scores (i.e decreasing number of occurrences found within the top-10 returned records). In this case, the success rate is measured by the percentage of query records where the correct event was retrieved among all suggested event tags. It measures the performance of our method in the context of tags suggestion rather than automatic annotation. Finally, we used the Hadoop 1 MapReduce implementation on a 5-node cluster. Nodes are equipped with Intel Xeon X5560 CPUs as well as 48Gb of RAM.

Results

1 http://hadoop.apache.org/mapreduce/

Parameters discussion

In figure 4, we report the mean average precision for varying values of the θ parameter (Eq. 3) and different numbers of K-nearest neighbors used during the visual matching step. The results show that MAP values are at their optimal for θ ∈ [300, 1800] seconds. This optimal error tolerance value is coherent with the nature of the events in the last.fm corpus. Concert's picture records indeed usually range from one to several hours. On the other hand, above 5 minutes, real-world concert scenes are too much ambiguous to be discriminated by their visual content (or at least with the global visual features used in this study). In what follows, we fix θ to 1800 as an optimal value for visual matching.

We now study the impact of adding a prior constraint δmax on the estimated temporal offsets δ. Most events in last.fm dataset being music concerts, it is unlikely that the temporal offset between two records reach high values. We therefore study the impact of rejecting all retrieved records having a temporal offset higher than δmax. Figure 5 displays the new MAP curves for varying values of δmax. It shows that the mean average precision can be consistently improved from about 0.08 without any constraint to 0.18. The optimal value for δmax is about 86,400 seconds which is exactly 1 day. That means that the records of a single realword event might have a temporal offset of up to 1 day. Our interpretation is that the EXIF creation date field might be noisy due to the different reference times of the used devices (users from different countries, etc.). It is worth noting that our method is by its very nature robust to such temporal offsets since we mainly consider temporal coherence rather than absolute time matching. On the other hand, rejecting records with temporal offsets higher than 1 day allows many visual false positives to be rejected. Figure 6 displays the results of the same experiment but for the classification rate (using a 10-NN classifier on retrieved records) rather than the mean average precision. This evaluates the ability of our method to automatically annotate a query event record rather than its ability to retrieve all records in the dataset. Here again the optimal classification rates are obtained when δmax=1 day. Furthermore, we see that the classification rate always increases with the number K of closest visual matches (returned for each query image). The interpretation is that increasing K improves recall without degrading precision too much thanks to the selectivity of our temporal consistency reranking step. We verified this from the results presented in this paper by studying the recall and the precision independently.

Event suggestion in the MapReduce framework

All the previous experiments were made using an exhaustive search for the k-NN search. In this section, we evaluate the performance of our full framework using MapReduce and the Multi-Probe LSH. As parameters, we used the optimal values discussed in the previous section (i.e. δmax=86.400, K=3000 and θ=1800).

Table 2 displays the class rates using an exhaustive search as seen in the previous section as well as class rates using a Multi-Probe LSH-based similarity search for different values of δMP .

As one might expect, all class rates values increase accordingly with the number of probes (i.e increasing δMP values) to surprisingly perform better than the exhaustive search for δMP =8.

Overall, in the best case, our method is able to suggest the correct event tag over 5 suggestions with a 72% success rate. Such performances are clearly acceptable from an application point of view. Figure 7 displays the average search time per query for both distributed and centralized search. We compare the K-NN processing time per image for a centralized setting (number of map slots = 1) to the processing time in a distributed scheme (20 map slots available on the network) for both exact and approximate similarity search. Although the multi-probe might reduce the effectiveness down to 66%, it might also reduce the search time by a factor of 13.

CONCLUSION AND FUTURE WORK

The achieved performance gain is nonetheless still less than the performance gain obtained in usual centralized settings. First, in MapReduce approaches, probing multiple buckets generates more network overhead in addition to data transfers across the network. The second reason is due to bucket occupation. In fact, imbalanced buckets generate imbalanced map chunks leading to disproportionate map execution times. Such a problem is addressed in [START_REF] Covell | Lsh banding for large-scale retrieval with memory and recall constraints[END_REF] and we plan to study such balancing methods in future work.

In this paper we presented a new visual-based method for retrieving events in photo collections, that might also be used for event tag suggestion or annotation. Our method proved to be robust to temporal offsets since we mainly rely on temporal coherence rather than absolute time matching. As one result, we are able to suggest the correct event tag with a success rate of at least 60% and even 72% if we allow multiple suggestions.

The proposed method is scalable, since it relies on efficient approximate similarity search techniques based on the MapReduce framework. We also investigated multi-probe techniques trading accuracy for efficiency, which might lead to a loss of 8.3% compared to a gain of 58.6%

Future work can focus on improving the suggestion rate as well as the efficiency of the approximate similarity search framework. The first issue can be addressed by including additional metadata during the re-ranking stage (notably geo-tags) and the use of more effective visual features. The second can be achieved through a better design of the hash functions to ensure a fair bucket occupation and therefore, balanced inputs for the map tasks.

We believe that the proposed method could have many other applications including other media event tracking or event mining in UGC's streams. Up to however we have found it difficult to collect relevant data for evaluation, and we work on that as well.

Figure 2 :

 2 Figure 2: Two records of the event "a trip in Egypt"

Figure 1 :

 1 Figure 1: Two events records of Alanis Morissette concert

Figure 3 :

 3 Figure 3: Processing time per image according to query size

Figure 4 :

 4 Figure 4: Influence of temporal error to tolerance θ

Figure 5 :Figure 6 :

 56 Figure 5: Influence of temporal offset thresholding (δmax) on MAP

Figure 7 :

 7 Figure 7: K-NN search time per image (k = 4000)

Table 1 : Test dataset Vs Original dataset

 1

		Total	Crawled Filtered
	photos 1 667 317 1637585 828902
	users	23 060	22676	10257

Table 2 : Suggestion rates

 2

	# of suggested	1	2	3	4	5	10
	events tags
	Exhaustive	0.60 0.66 0.69 0.71 0.72 0.73
	MP-Delta 0	0.39 0.48 0.50 0.51 0.52 0.54
	MP-Delta 1	0.45 0.55 0.57 0.58 0.59 0.63
	MP-Delta 2	0.48 0.59 0.61 0.65 0.66 0.69
	MP-Delta 4	0.51 0.61 0.63 0.67 0.67 0.70
	MP-Delta 8	0.61 0.67 0.70 0.72 0.72 0.74
	MP-Delta 16	0.59 0.66 0.69 0.71 0.72 0.73

http://www.glocal-project.eu/

ACKNOWLEDGMENTS

This work was funded by the EU through the Integrated Project GLOCAL 2 (contract number FP7 -248984).


M Carmona
email: mikael.carmona@cea.fr

O Michel

J-L Lacoume

B Nicolas

N Sprynski

PASSIVE IDENTIFICATION OF ACOUSTIC PROPAGATION MEDIA WITH VISCOUS DAMPING

This paper concerns passive identification of an acoustic propagation medium with viscous damping. We propose a Green correlation approach to retrieve time of arrival between two sensors and to retrieve also times of arrival of first echoes. Green correlation is introduced as the crosscorrelation of acoustic field generated by a white noise exciting the medium. This approach is experimentally validated and compared to the classical method which consists in retrieving the Green function from the Green correlation through a Ward identity. We show that this latter approach appears to be numerically inaccurate and above all unnecessary in the acoustic case with viscous damping.

1.

Passive identification of a propagation medium consists in retrieving medium parameters, by using uncontrolled noise fluctuations only [START_REF] Snieder | Unified Green function retrieval by cross-correlation; connection with energy principles[END_REF]. Such an idea has long been pursued in acoustics [START_REF] Lobkis | On the emergence of the Green function in the correlations of a diffuse field[END_REF], [START_REF] Snieder | Unified Green function retrieval by cross-correlation; connection with energy principles[END_REF] and seismology [START_REF] Campillo | Long-Range Correlations in the diffuse seismic coda[END_REF], [START_REF] Wapenaar | Tutorial on seismic interferometry. Part II: Underlying theory and new advances[END_REF] and gave rise to numerous applications and experimental validations [START_REF] Sabra | Extracting time-domain Green function estimates from ambient seismic noise[END_REF], [START_REF] Gouedard | Cross-correlation of random fields: mathematical approach and applications[END_REF].

Preceding studies [START_REF] Lobkis | On the emergence of the Green function in the correlations of a diffuse field[END_REF], [START_REF] Campillo | Long-Range Correlations in the diffuse seismic coda[END_REF], [START_REF] Sabra | Extracting time-domain Green function estimates from ambient seismic noise[END_REF], [START_REF] Gouedard | Cross-correlation of random fields: mathematical approach and applications[END_REF] rely upon the estimation of the Green function of the medium. Such estimation is made possible by exploiting the celebrated Ward identity [START_REF] Weaver | Ward identities and the retrieval of Green's functions in the correlations of a diffuse field[END_REF], which relates the noise correlation function to the Green function [START_REF] Lobkis | On the emergence of the Green function in the correlations of a diffuse field[END_REF], [START_REF] Snieder | Unified Green function retrieval by cross-correlation; connection with energy principles[END_REF], [START_REF] De Verdiere | Semi-classical analysis and passive imaging[END_REF]. The fundamental role of dissipation in Ward identity was outlined in [START_REF] Gouedard | Cross-correlation of random fields: mathematical approach and applications[END_REF], were dissipation is assumed to be constant ; however, a constant dissipation model is hardly acceptable from a physical point of view [START_REF] Landau | Kosevich Physique thorique : thorie de l'lasticit[END_REF], [START_REF] Royer | Elastic Waves in Solids[END_REF], and needs to be further discussed.

The contribution of the paper is to study and to validate experimentally a Green correlation approach in passive identification of acoustic media with viscous damping. Green correlation is introduced as the correlation of propagated white noise fields [START_REF] Lacoume | Tomographie passive: observer avec du bruit[END_REF]. The motivation for introducing Green correlation comes from classical system identification theory, where noise based identification relies strongly on the transformation of second order statistics through linear systems [START_REF] Ljung | System Identification: Theory for the User[END_REF].

To interpret acoustic Green correlation in bounded media, we derive expressions for the acoustic Green function and we relate them to the Green correlation by deriving Ward identities. Acoustic Green function is derived using Fourier theory. Low attenuation case is considered as it allows to obtain an explicit expression of the Green function and an explicit Ward identity in the temporal domain.

Those results are compared theoretically to classical acoustic results where constant damping model is considered [START_REF] Lacoume | Tomographie passive: observer avec du bruit[END_REF], [START_REF] Gouedard | Cross-correlation of random fields: mathematical approach and applications[END_REF], [START_REF] De Verdiere | Semi-classical analysis and passive imaging[END_REF]. Secondly, results are experimentally validated by showing that Green correlation allows to retrieve time of arrival between two sensors and also time of arrival of first echoes. We compared this approach to classical method which consists to retrieve the Green function from the estimated Green correlation and a Ward identity.

The organization of the paper is the following. In section 2, we introduce passive identification through a linear systems approach. We recall the definition of the Green function and its role in medium identification. Cross-correlation of random fields is recalled, and white noise notion is introduced. Then, we defined the Green correlation and we show its role in passive identification.

In section 3, acoustic waves equation with viscous damping is recalled. The acoustic Green function is computed for unbounded and bounded media. Low attenuation case is also considered.

Acoustic Green correlation is computed in section 4. Role of Green correlation in passive identification is emphasized. Ward identities are derived and compared to existing one for a constant damping model. Approximated Ward identities are derived in the low attenuation case.

In section 5, we show through an experimental validation that the Green correlation is a powerful tool in practical situation, especially when theoretical Ward identity can not be used due to signals digitalisation letting the Green function estimation unusable.

GREEN FUNCTION AND GREEN CORRELATION OF A LINEAR PROPAGATION MEDIUM : DEFINITION AND ROLE IN MEDIUM IDENTIFICATION

2.1 Green function of a linear medium.

We denote by u(t, x) the value of field u at time t and position x. A time-shift invariant linear propagation medium X satisfies the superposition theorem i.e. the value u(t, x) of the generated field can be seen as the superposition of all contributions of elementary sources f (t , x )dt dx emitted at time t during dt period in the volume centered in x and of dimensions dx , for all times t and points x . Mathematically, this can be written u(t, x) = R×X G(t -t , x, x )f (t , x )dt dx and simplified by introducing the generalised convolution ⊗ T,S as:

u = G ⊗ T,S f (1)
G is the Green function of the medium. All medium parameters are contained in its expressions. This highlights the importance of retrieving the Green function in medium identification.

Physically, G corresponds to the medium response of a spatio-temporal impulsion. According to this interpretation, G is sometimes called "impulse response" of the medium in reference to the classical impulse response of a linear system [START_REF] Ljung | System Identification: Theory for the User[END_REF]. In that system approach, f is the input and the generated field u is the output of the system.

Cross-correlation of stochastic fields.

In passive identification, source fields are not controlled. The principle relies on recording noise sources and using their statistical properties to retrieve medium parameters. With stochastic source fields, the analysis has to be performed from the cross-correlation of the generated field u defined as:

C u (t, x, x ) = lim τ→+∞ 1 τ τ 0 u(t + t , x)u(t , x )dt (2)
This formula is true for stationary and ergodic fields. This assumption is not a strong constraint in practice.

White noise.

By definition, a white noise is a field which value of a given time and position is uncorrelated to any other value taken at all other times and positions. Mathematically, the cross-correlation of a such field f is a spatio-temporal impulsion:

C f (t, x, x ) = δ (t)δ (x, x ) (3)
where δ is the Dirac distribution.

A white noise has no physical reality because it has an infinite power. However, in practice the temporal whiteness is only needed in a limited frequency band. This latter is defined by the used instrumentation. The classical approach [START_REF] Campillo | Phase and Correlation in "Random" Seismic Fields and the Reconstruction of the Green Function[END_REF] to justify that ambient noise converges to a source with a spatial whiteness consists to see the medium as a chaotic dynamical system. Then, according to equipartition theorem, it exists a time after which a coherent source snared in the medium becomes spatially white. This time, called mixing time, depends on the frequency band, medium geometry and heterogeneity. With those considerations, cross-correlation of ambient noise is stacked during a sufficient long time in order to obtain a contribution of an approximated white noise [START_REF] Lobkis | On the emergence of the Green function in the correlations of a diffuse field[END_REF], [START_REF] Campillo | Long-Range Correlations in the diffuse seismic coda[END_REF], [START_REF] Gouedard | Cross-correlation of random fields: mathematical approach and applications[END_REF].

Green correlation.

As the Green function is the field generated by a spatiotemporal impulsion source i.e. f (t, x) = δ (t,t )δ (x, x ), we define by analogy the Green correlation C by the crosscorrelation of a field generated by a white noise source. This function, introduced in [START_REF] Lacoume | Tomographie passive: observer avec du bruit[END_REF], plays by definition a fundamental role in passive identification.

We can precise Green correlation expression using definition (2) and equation [START_REF] Anderson | Optimal Filtering[END_REF]. Indeed, for every generated field u, C u can be expressed as:

C u = G ⊗ T,S G -⊗ T,S C f (4)
where G -(t, x, x ) := G(-t, x, x ). It is important to note that to establish (4), we use the property: G(t, x, x ) = G(t, x , x), true for all times t and positions x according to spatial reciprocity. Equation ( 4) is the "order two" version of equation [START_REF] Anderson | Optimal Filtering[END_REF]. When the source is a white noise, we obtain by using equation ( 4) an expression of the Green correlation:

C := G ⊗ T,S G - (5)
This shows the fundamental importance of the Green correlation in medium identification when statistical properties of ambient sources are taking into account. Equations ( 4) and [START_REF] Campillo | Phase and Correlation in "Random" Seismic Fields and the Reconstruction of the Green Function[END_REF] show that "perfect" white noise is for passive identification what "perfect" impulsion is for active identification.

ACOUSTIC PROPAGATION WITH VISCOUS

DAMPING : EQUATION AND GREEN FUNCTION.

3.1 Equation.

We consider an isotropic and homogeneous acoustic medium. Let v be the sound velocity and α be the viscous damping coefficient. Let f be a causal spatio-temporal pressure source and u the pressure field. The acoustic waves equation with viscous damping is then [START_REF] Royer | Elastic Waves in Solids[END_REF]:

∂ 2 ∂t 2 -α 2 ∂ ∂t ∆ -v 2 ∆ u = f ( 6
)
where ∆ is the Laplacian operator. To well-define this equation we have to add fields causality assumption and boundary conditions.

Green function in an unbounded medium.

For an unbounded medium, the acoustic Green function is space-shift invariant i.e. we can do the following substitution: G(t, x, x ) ↔ G(t, xx ). The case of a bounded medium is discussed in next sub-section.

We denoted by

G(ω, k) := R×R 3 G(t, x)e -i(ωt-k T x) dtdx (7)
the Fourier transform of G in the (ω, k)-domain where T transpose, and, where ω and k are the frequency variables associated with t and x, respectively.

As the Green function is the response to a spatiotemporal impulsion source i.e f := δ (t, x), we obtain from equation ( 1

): Ĝ(ω, k) = ω -2 k 2 (ω) k 2 (ω) -k 2 (8)
where:

k 2 (ω) := ω 2 v 2 (ω) (9)
with:

v 2 (ω) := v 2 + iωα 2 , k 2 := k T k. k(ω)
is the relation of dispersion of acoustic waves with viscous damping. When low attenuation is considered i.e. ωα 2 /v 2 << 1, we get:

k(ω) ≈ ω v (10)
This approximation of k(ω) will be useful to derive expression of the Green function and correlation and also Ward identities in the (t, x)-domain.

We compute now the elastic Green function in the (ω, x)domain. In that domain, fields are capped by a ∨. Applying inverse Fourier transform to equation ( 8) with respect to k, we show using residue theorem [START_REF] Arfken | Mathematical Methods for Physicists[END_REF] that:

Ǧ(ω, x) = e -i||x||k(ω) 4π||x||v 2 (ω) ≈ e -iω ||x|| v 4π||x||v 2 (11)
The approximation holds for a low viscous damping and allows to retrieve the classical Green function when there is no damping. It gives an easy interpretable expression of the Green function as a pure phase which provides an information on ||x||/v easy to extract.

Computation of a general expression for G(t, x) from equation [START_REF] Lobkis | On the emergence of the Green function in the correlations of a diffuse field[END_REF] when dissipation occurs, is difficult. In the low attenuation case, we obtain from [START_REF] Lobkis | On the emergence of the Green function in the correlations of a diffuse field[END_REF]:

G(t, x) ≈ δ (t -||x||/v) 4πv 2 ||x|| (12)
This show that retrieving [START_REF] Royer | Elastic Waves in Solids[END_REF] allows to estimate the time of arrival ||x||/v.

Green function in a bounded medium with a low viscous damping.

In a practical context, the medium is generally bounded. From equation [START_REF] Royer | Elastic Waves in Solids[END_REF] and Descartes' laws, we can decomposed the acoustic Green function in a sum of attenuated pure delays i.e.:

G(t, x, x ) = +∞ ∑ n=0 a n (x, x )δ (t -t n (x, x )) (13)
where a 0 (x, x ) := 1/(4π||xx ||), t 0 := ||xx ||/v, and, a n (x, x ) and t n (x, x ) are attenuation coefficients and echoes time of arrival, respectively. In practice, the sum appearing in ( 13) is a finite sum due to the decreasing of attenuation coefficients and because of quantification which provides only a finite set of values for fields recorded.

ACOUSTIC GREEN CORRELATION AND WARD IDENTITIES.

Unbounded case.

As the Green function is space-shift invariant in an unbounded medium, the Green correlation satisfies also this property according to [START_REF] Campillo | Phase and Correlation in "Random" Seismic Fields and the Reconstruction of the Green Function[END_REF]. Then, we can do the following substitution:

C(t, x, x ) ↔ C(t, x -x ).
We can compute the acoustic Green correlation in the (ω, k)-domain. Applying the Fourier transform to equation [START_REF] Campillo | Phase and Correlation in "Random" Seismic Fields and the Reconstruction of the Green Function[END_REF] and noting that G -= G † where † conjugate, we obtain:

Ĉ(ω, k) = Ĝ(ω, k) Ĝ(ω, k) † (14)
We can deduce a Ward identity in the (ω, k)-domain from equation ( 14) and the following easy to prove lemma:

Lemma. Let a be a non real complex. Then we have:

aa † = -Im a -1 -1 Im a.
Applying this lemma with a = Ĝ(ω, k), we obtain:

ω Ĉ(ω, k) = - 1 α 2 k 2 Im Ĝ(ω, k) (15)
In that domain, the acoustic Green correlation is proportional to the imaginary part of the Green function. The proportionally term 1/(α 2 k 2 ) is the inverse of the dispersion function of α 2 ∆ which is the dissipation operator appearing in equation ( 6). From Ward identity (15), we can easily derived exact Ward identities in (ω, x)and (t, x)domains:

ω Č(ω, x) = -α -2 ∆ -1 Im Ǧ(ω, x) (16)
∂ C(t, x) ∂t = -α -2 ∆ -1 Odd G(t, x) (17)
where Odd G := 1/2 (G -G -) is the odd part of G. Those identities extend classical ones for acoustic waves [START_REF] Lobkis | On the emergence of the Green function in the correlations of a diffuse field[END_REF], [START_REF] Gouedard | Cross-correlation of random fields: mathematical approach and applications[END_REF], [START_REF] De Verdiere | Semi-classical analysis and passive imaging[END_REF], [START_REF] Lacoume | Tomographie passive: observer avec du bruit[END_REF] in a sense that the Č is "proportional" to the imaginary part of Ǧ and the first-time derivative of C is "proportional" to the odd part of G. However, when viscous damping is considered, the proportionality term is ∆ -1 and it is not a constant operator. Then, ( 16) and ( 17) are not sufficiently explicit to interpret the relation existing between the Green function and the Green correlation. We get over this last step by considering a low viscous damping.

Unbounded and low viscous damping case.

First, we compute the acoustic Green correlation in the (ω, k)-domain. Using Ward identity [START_REF] Wapenaar | Tutorial on seismic interferometry. Part II: Underlying theory and new advances[END_REF] and expression (8), we get:

Ĉ(ω, k) = ω -4 |k(ω)| 4 |k 2 (ω) -k 2 | 2 (18)
Expression (18) proves that retrieving the Green correlation is sufficient to estimate medium parameters: v and α. Indeed, those parameters are contained in k 2 (ω) which is the pole of Ĉ(ω, k). Inverse Fourier transform of equation ( 18) with respect to k and residue theorem give:

Č(ω, x) = - Im e -i||x||k(ω) 4π||x||ω 3 α 2 ≈ sin ω ||x|| v 4π||x||ω 3 α 2 (19)
This expression shows that we can extract all the physical parameters from Č. Using relations [START_REF] Lobkis | On the emergence of the Green function in the correlations of a diffuse field[END_REF] and ( 19), we obtain the following approximated Ward identities:

ω 3 Č(ω, x) ≈ - v 2 α 2 Im Ǧ(ω, x) (20)
∂ 3 C(t, x) ∂t 3 ≈ - v 2 α 2 Odd G(t, x) (21)
(21) was deduced from (20). For a viscous damping and low attenuation framework, it is the third-time derivative which is directly proportional to the odd part of the Green function in the (t, x)-domain. The difference derivative order comes only from the form of the dissipation.

Bounded and low viscous damping case.

We show now how Green correlation can be useful in passive identification of bounded media. It is difficult to obtain an explicit expression of the Green correlation from equations ( 5) and ( 13). However, Ward identity (21) can be extended to any bounded media as:

∂ 3 C(t, x, x ) ∂t 3 ≈ - v 2 α 2 Odd G(t, x, x ) (22)
Indeed, only the propagation equation ( 6) intervenes in computation of Ward identities. Boundary conditions only intervenes in Green function and Green correlation expressions. The equation ( 22) is important as it proves that times of arrival appearing in the Green function also appear, symmetrically with respect to time origin and with an opposite amplitude, in the Green correlation. This justifies that Green correlation estimation is theoretically sufficient to retrieve times of arrival of direct and reflected waves.

EXPERIMENTAL VALIDATIONS.

We consider two microphones, denoted by m 1 and m 2 measuring pressure field in the 2 Hz-17 kHz frequency band. Every measure is sampled at 51.2 kHz and high-pass filtered at 200 Hz in order to suppress electrical noises. Sensors are disposed in an acoustic air room and spaced of 26.5 cm. The experimental protocol relies on two steps: to retrieve the Green function by active identification and to retrieve the Green correlation by passive identification. Then, results are discussed and compared.

In the first step, m 1 is seen as a source and m 2 as a receiver. Indeed, an acoustic signal is emitted from a speaker closely to sensor m 1 . Five seconds of signals recorded by m 1 and m 2 are used to retrieve the Green function of the medium between the two sensors. We use a classical recursive least square algorithm [START_REF] Anderson | Optimal Filtering[END_REF] to retrieve the finite impulsion response filter (1024 points, 20 ms) representing numerically the Green function. The convergence really appears after 500 ms of signals recorded.

In the second step, m 1 and m 2 are both seen as receivers. Speakers emitting temporal white noise, generated numerically, where displaced around the two sensors in order to "simulate" a spatio-temporal white noise. Signals received by m 1 and m 2 where used to compute cross-correlations on 20 ms. All correlations where stacked during 5 s to obtain the estimated Green correlation. The convergence really appears after 1 s.

Odd part of the estimated Green function, and, the estimated Green correlation are represented and compared in Figure 1. We normalised each estimated field by dividing by their maximum in order to not take their amplitude into account.

The estimated Green function contains its first significant pic at 0.74 ms which corresponds to a speed sound of v = 348 m/s which is an acceptable value at 22 Furthermore, we observe echoes justified by the reflections on room interfaces, this valid the model proposed by equation [START_REF] Sabra | Extracting time-domain Green function estimates from ambient seismic noise[END_REF]. We observe a similar evolution for the estimated Green correlation. There is not a perfect symmetry of the Green correlation due to the imperfect symmetry of ambient sources distribution used. We retrieve the time of arrival at ±0.74 ms with an opposite amplitude for positive and negative lags. First echoes match with the odd part of the Green function. Then, for lag t such |t| > 1 ms echoes do not really correspond between Green function ones and Green correlation ones.

• C [3].
We can conclude on this experimentation that we retrieve perfectly time of arrival of the direct path. We observe a drift on times of arrival of echoes. This is not a problem to identify sound velocity or distance sensor-sensor. However, geometrical identification of the medium from time of arrival of echoes appearing in the Green correlation is clearly less accurate. This result can also be observed in [START_REF] Lobkis | On the emergence of the Green function in the correlations of a diffuse field[END_REF].

We compare in Figure 2 the odd part of the Green function with the first-time derivative of the Green correlation, and, with the third-time derivative of the Green correlation. We see, naturally, that the more the time derivative order increases the more the signal to noise ratio decreases and times of arrival estimation becomes difficult. This proves that Ward identity is an important theoretical tool to compute and interpret the Green correlation from the Green function but it is difficult to apply numerically and above all unnecessary.

CONCLUSION.

We applied Green correlation theory to acoustic waves with viscous damping. We computed, in useful representation domains, the Green function, the Green correlation and Ward identities associated. We considered low attenuation case, which provides interesting and interpretable Ward identities in (ω, x)and (t, x)-domains. More precisely, in this latter domain, it is the third-time derivative of the elastic Green correlation which is proportional to the odd part of the elastic Green function. This is a direct consequence of viscous damping model and this result was validated experimentally and compared to classical ones in acoustic with a constant damping model. We prove that retrieving the Green correlation is sufficient to estimate the direct time of arrival and first times of arrival of echoes (before 1 ms for our experimentation). This confirms that Green correlation is an interesting practical tool in passive identification.

Figure 1 :

 1 Figure 1: Representation of the estimated Green function and the estimated Green correlation. The two functions are compared in last sub-figure.

Figure 2 :

 2 Figure 2: Comparison of the odd part of the Green function with the first-time derivative of the Green correlation, and with the third-time derivative of the Green correlation.


Ahmed Rebai
email: ahmed.rebai@inria.fr

Alexis Joly
email: alexis.joly@inria.fr

Nozha Boujemaa
email: nozha.boujemaa@inria.fr

Interpretable Visual Models for Human Perception-Based Object Retrieval

Keywords: I.4.8 [IMAGE PROCESSING AND COMPUTER VI-SION]: Scene Analysis-Object recognition Object retrieval, interpretability, feature selection, sparsity, human perception, interactivity

Understanding the results returned by automatic visual concept detectors is often a tricky task making users uncomfortable with these technologies. In this paper we attempt to build humanly interpretable visual models, allowing the user to visually understand the underlying semantic. We therefore propose a supervised multiple instance learning algorithm that selects as few as possible discriminant local features for a given object category. The method finds its roots in the lasso theory where a L1-regularization term is introduced in order to constraint the loss function, and subsequently produce sparser solutions. Efficient resolution of the lasso path is achieved through a boosting-like procedure inspired by BLasso algorithm. Quantitatively, our method achieves similar performance as current state-of-the-art, and qualitatively, it allows users to construct their own model from the original set of patches learned, thus allowing for more compound semantic queries.

INTRODUCTION

Understanding an image scene is definitely related to recognizing what it is composed of. Unfortunately, in practice, the results returned by state-of-the-art visual concept detectors are often difficult to interpret from a user point of view. The visual models produced are indeed highly depending on the used training data and might convey a different semantic than the words used to describe the originally targeted concept. This often makes users uncomfortable with these technologies since they do not get what they expected from the textual description of the trained concept. In this paper we rather attempt to build humanly interpretable visual models, allowing the user to visually understand the underlying semantic.

Our basic idea is to build concise, effective and visually interpretable models to allow for fast retrieval. Each object model consists of a concise set of visual patches representing the most discriminant image regions of a concept. Consequently, users can easily visualize the trained image regions and query the system by only choosing the visual patches that most correspond to their needs. This involves adding, subtracting and attributing weights to the visual patches displayed.

Composing a statistically and visually discriminant model provides a unique way for an interactive search. Now, the challenging question to answer is how to train the most concise set of visual patches while preserving a good retrieval effectiveness. The key to this question is provided by the image representation and the learning technique used. Using state-of-the-art bag-of-features technique is for instance not adapted because it discards the positions of the features being learned. These visual words do not provide any guaranty on whether the underlying clusters pertain to tangible parts of the objects (i.e eye, tooth, finger, etc.) or if they are just a statistical combination of some of these parts. On the other hand, using learning techniques based on such features with a SVM classifier is also not adequate for the matter since our objective is to locate a few representative image regions rather than creating an optimal separation between the positive and negative feature sets. It is worth mentioning that the method we propose is generic in the sense that it could be used with any local features or feature sets representing the content of an image region. In this work, we simply used SIFT features as the basic image primitives.

The main contribution of this paper is a novel multipleinstance learning method that allows for user interactivity when retrieving visual concepts. Applying the principle of the lasso technique to a discriminative approach for concept retrieval constitutes our second contribution. The third contribution consists in adapting the BLasso algorithm to our purpose by using weighted training data and constraining the forward steps to be only positive steps. As a last contribution, we propose an efficient implementation for our learning method to speed up the process of backward steps. The next section briefly reviews state-of-the-art feature selection methods related to our work. Section 3 explains our learning strategy in more detail. Then, we present our experiments in section 4 and set out our conclusions in section 5.

FEATURE SELECTION

Feature selection has always been an overwhelming task. The method we propose is based on a modified version of the BLasso algorithm. Our original idea is to train with a boosting-like procedure but with adding extra constraints to the loss function in order to produce sparser solutions. Added to highlighting interest features, sparsity is indeed preferable because it reduces the model complexity and subsequently enhances the interpretability of the produced models as well as it reduces the prediction time.

Constraining the empirical loss dates back to 1996 when Tibshirani observed that the ordinary least squares minimization technique is not always satisfactory since the estimates often have a low bias but a large variance. He came out with lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] which shrinks or sets some coefficients to zero. Lasso stands for Least Absolute Shrinkage and Selection Operator. The idea has two goals: first to gain more interpretation by focussing on relevant predictors and, secondly to improve the prediction accuracy by reducing the variance of the predicted values. Lasso minimizes the L2 loss function penalized by the L1 norm on the parameters. This is a quadratic programming problem with linear inequality constraints and it is intractable when the vector of parameters is very large.

In literature, some efficient methods have been proposed to solve the exact Lasso namely the least angle regression by Efron et al. [START_REF] Efron | Least angle regression[END_REF] and the homotopy method by Osborne et al. [START_REF] Osborne | On the lasso and its dual[END_REF]. These methods were developed specifically to solve the least squares problem (i.e. using L2 loss). They work well where the number of predictors is small. However, they are not adapted to nonparametric and classification tasks. To encounter these limitations, Zhao and Yu [START_REF] Zhao | Stagewise lasso[END_REF] proposed an algorithm with the name of BLasso (Boosted Lasso). The advantage of this algorithm lies in its ability to function with an infinite number of predictors and with various loss functions. These are some characteristics of the algorithms belonging to the boosting family.

The boosting mechanism was proposed by Schapire [START_REF] Schapire | The strength of weak learnability[END_REF] in 1990. Since then, many algorithms have emerged [START_REF] Freund | A decision-theoretic generalization of on-line learning and an application to boosting[END_REF][START_REF] Friedman | Additive logistic regression: a statistical view of boosting[END_REF][START_REF] Opelt | Generic object recognition with boosting[END_REF][START_REF] Singer | Leveraged vector machines[END_REF][START_REF] Tieu | Boosting image retrieval[END_REF] and boosting has become one of the most successful machine learning techniques. The underlying idea is to combine many weak classifiers (also called hypotheses) in order to obtain one final "strong" classifier. Boosting is an additive model which builds up one hypothesis after another by re-weighting the data for the next iteration-increasing the weights of misclassified images and decreasing those of well classified ones. This concept helps to generate different hypotheses, putting emphasis on misclassified examples, typically those located near the decision boundary in the feature space. Although it is an intuitive algorithm, boosting may overfit the training data, particularly when it runs for a large number of iterations T in high dimensional and noisy data [START_REF] Grove | Boosting in the limit: Maximizing the margin of learned ensembles[END_REF][START_REF] Rätsch | Soft margins for adaboost[END_REF]. That's why, adding a regularization term is usually needed.

Unlike boosting, and in order to approximate lasso solutions, BLasso adds a backward step after each iteration of boosting. Thus, one is able to build up solutions with a coordinate descent manner and then take a look back at the consistency of these solutions regarding the model complexity. Added to that, BLasso does not suffer from the stop condition like an ordinary boosting does. In fact, and for a boosting procedure, fixing a large value of T implies a possible overfitting not to mention the long prediction time. On the other hand, setting T to a small value may lead to underfitting. Therefore, the model may be non-discriminant, inconsistent and might not cover the variability inside the category itself. A usual boosting algorithm can also be qualified as oblivious as it always functions in a forward manner aiming to minimize the empirical loss. Although the concept of re-weighting is interesting, at an iteration t + 1, we have no idea whether the t previous generated hypotheses are good enough or not versus the model complexity.

LEARNING VISUAL PATCHES

Multiple-Instance Learning with BLasso

Training a Weak Classifier

Background information sometimes plays a primordial role in recognition. Thus, we propose to train our classifiers in a context of weak learning. That is, a training image is labeled as a whole sample. It will take the label +1 if it contains the object, -1 if not. This is also known as multiple instance learning. It deals with uncertainty of instance labels. An image is viewed as a bag of multiple features which are the local visual signatures. The bag will have only one label according to whether or not it includes at least one positive instance. It follows that it is only certain for a negative bag that there are no objects. Moreover, using a multi-instance approach has the luxury of unsupervised learning. It gives more freedom to the algorithm to select background features whenever they turn out to be useful to characterize the category.

Visual patches are viewed as weak classifiers. A weak classifier h k represents a coordinate of base learners. Its weight is strictly positive if it was chosen at least once during the boosting process and remains zero if not. In the context of object categorization, the weak hypothesis h k is a local image signature F k with an optimal radius r k (Opelt et al. [START_REF] Opelt | Generic object recognition with boosting[END_REF]). In other words, h k corresponds to a hypersphere centered on a local image feature F k . For a test image x, h k will output +1 if the distance between x and F k is less than r k and -1 otherwise. Discriminant radii of weak classifiers are optimal in the sense that the classification error is as minimal as possible. They are determined through a distance matrix. Each entry of the matrix represents a local feature with a corresponding ranked list of the training images. The distance between any feature F ki belonging to the image Ii and any image Ij of the training set is defined by:

d(F ki , Ij) = min 1≤k ≤M j d(F ki , F k j )
Let Si = (Ii, li) be the couple representing the training image Ii labeled with li ∈ {-1, 1} and wi be the weight associated to Ii. S = {S1, • • • , SN } represents the set of all the training data. We denote by d k,η the distance separating the feature F k from the image whose rank is η. The optimal radius r k is determined after two steps. First we compute the curve corresponding to the sums of the weighted image labels and take the index where the maximum is achieved:

η = arg max 1≤η≤N η i=1 wi • li (1)
Then, r k is given by:

r k = d k,η + d k,η+1 2 (2)
Training a weak classifier is a simple boosting iteration. At an iteration t, we compute the score s k corresponding to the feature F k based on the image weights

s k = max m m j=1 wj • lj (3)
Then, we select the feature which obtains the highest score. Note that Eq. ( 1) is in accordance with Eq. ( 3) in the sense that Eq. ( 1) looks for the index at which the score curve is maximum.

Using AdaBoost algorithm, Opelt et al. [START_REF] Opelt | Generic object recognition with boosting[END_REF] considered an infinity of base learners. At the iteration t, and after selecting the best feature F k , the base learner associated to F k is defined by the couple (F k , r k ) and is attributed a weight according to the training error. (r k is given by eq. 2.) This seems plausible because AdaBoost uses a steepest descent to converge.

In our work, and for the sake of interpretability, each feature F k represents only one base learner and should at most be selected once. It has to have its weight determined after several boosting iterations (forward stagewise). However, whenever selected, this base learner sees its radius r k changes according to the weight updates, thus the final r k attributed is the average of all the computed radii.

The Algorithm

Our algorithm is viewed as a boosting method constrained by a L1-regularization term. L1-regularization is equivalent to lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. Let β = (β1, • • • , β k , • • • ) T be the vector of parameters to estimate (the weights of the weak hypotheses). β is initially zero. The lasso loss function Γ can be written as:

Γ(β, λ) = N n=1 L2(Sn, β) + λ • ||β||1 (4)
where ||β||1 = j |βj | denotes the L1 norm of the vector β and λ ≥ 0 is the parameter controlling the amount of regularization applied to the estimate. In order to get sparse solutions with an efficient shrinkage tradeoff, λ usually takes a moderate value since a large value may set these coefficients to exactly zero, leading to the null model. On the other hand, setting λ to zero reverses the lasso problem to minimizing the unregularized empirical loss. In our case, for classification, we replace L2-loss function by the exponential loss function.

Since the exact lasso minimization is not tractable, BLasso [START_REF] Zhao | Stagewise lasso[END_REF] tries to find the same solutions as lasso with more cautious steps. Indeed it works with both forward and backward steps. Forward steps are used to minimize the empirical loss. On the other hand, backward steps minimize the regularization. In fact, at each iteration, a coordinate βj is selected and updated by a small step size ε > 0. It has been shown that it is preferable to choose a very small step size so that BLasso can approximate the lasso path perfectly. In practice, ε should always be less than 0.1. In the original BLasso algorithm, forward steps can be either positive or negative (update by ±ε). However, we chose not to. Our forward steps are always positive. We justify this choice by the fact that a selected classifier is a visual patch that contributes to build the object category model, thus it has to be positive. Algorithm 1 BLasso 1. Initialization: β = 0 Make a forward step and initialize λ 2. Backward and forward steps:

Find the backward step that leads to the minimal empirical loss. if the step decreases the lasso loss then take it. else make a forward step and relax λ if necessary 3. Repeat step 2 until λ ≤ 0.

Algorithm 1 gives a general overview on the BLasso mechanism. In our implementation, and in order to minimize the empirical loss, we used a weighting scheme as in AdaBoost. Adopting this strategy is appealing because it gives more attention to the misclassified observations by increasing their respective weights. Note that, at initialization, all weights are equal to 1/N .

At the iteration t, a forward step can be summarized in the next five points.

Update weights

w (t+1) n = w (t) n • exp -ε • ln • h (t) κ I(n) τ (5)
τ is a normalization constant s.t. A ∪ {κ} where IA is the active index set. Note that equation 5 is in accordance with the weight update mechanism used in AdaBoost. In fact, we just replaced the estimated hypothesis weight by the forward step ε.

The variable ξ, used when updating λ, is a tolerance parameter that is strictly positive. This parameter is added to gain more stability and it should be set to a very small value. Moreover, the initial value of λ can be obtained by the same formula (eq. 6) but with omitting ξ (as if ξ = 0). ξ is also used to decide whether to accept or reject a backward step. A backward step consists in finding the step that leads to the minimal empirical loss. It is given by: ĵ = arg min

j∈I (t) A N n=1 L(Sn, β (t) -ε • 1j) ( 7
)
This step is taken if and only if it decreases the lasso loss.

That is, if Γ(β (t) -ε • 1 ĵ , λt) -Γ(β (t) , λt) ≤ -ξ then β (t+1) = β (t) -ε • 1 ĵ ; λt+1 = λt

Efficient Implementation

The inescapable processing burden for training is the computation of the distance matrix. To accelerate the process, one solution consists in trading accuracy for time by relying on approximate rather than exhaustive search. This could be achieved through range queries with an appropriate index structure like LSH for example. However, extreme care should be taken because it may drastically lower the quality of weak classifiers. This paper doesn't cover the gains obtained with such optimization.

The next processing burden and greediest block during training is computing the loss function. In fact, we need to compute both of the empirical and lasso losses many times during each iteration. The complexity increases every time a new coordinate is chosen. Luckily, the lasso loss can be deduced from the empirical loss. It follows that, in order to compute the current lasso loss and the backward lasso loss, we just need to compute their respective empirical losses. On the other hand, and at each iteration, only one coordinate is modified while all the other coordinates remain unchanged. Therefore, assuming that memory is cheaper than processing time, the computation of the empirical loss can be accelerated-at an iteration t-by preserving in memory the classification values of the previous state for each coordinate βj and for each image In. Since our solutions are sparse, even when using a large number of training images, this method is still applicable. Each image In will have a storage vector ζn. To simplify the notation, we will consider that the j th entry of the vector ζn (i.e. ζn(j)) refers to the classification value of the variable βj

ζn(j) = βj • hj (In) 1≤ n ≤ N ; j ∈ IA (8)
Thus, when changing a hypothesis h k in the next iteration, we only need to update (if β k has already been selected before) or create new (if the index k is selected for the first time) N classification values. That is, we have to compute ζn(k) for all 1 ≤ n ≤ N

The storage of the classification values is more important when searching for the backward step. The empirical loss is not computed just once but at least card{IA} times where card{IA} is the cardinal of the index set (It is computed card{IA} + 1 times if the coordinate is selected for the first time.) Given an image In, for each coordinate j, we need to compute the empirical loss based on the coordinates β -ε1j . Except for the coordinate βj , the other classification values have already been computed and stored in ζn(k) with k = j. Moreover, the absolute difference between ζn(j) (already computed if the hypothesis is old) and the classification value δnj that we need to compute is ε

if ζn(j) < 0 ⇒ δnj = ζn(j) + ε ζn(j) > 0 ⇒ δnj = ζn(j) -ε (9)
It follows that the classification value of the image In (i.e.

k ζn(k)) will only change by an absolute difference of ε. We denote by ζ n (j) = k =j ζn(k) + δnj , that is:

ζ n (j) = k ζn(k) + ε if ζn(j) < 0 k ζn(k) -ε if ζn(j) > 0 (10)
This formulation helps to locate the backward step quickly.

In fact, it only takes a linear time according to the number of the selected hypotheses. After subtracting the value ε from the coordinate βj , the empirical loss Ej is computed as follows:

Ej = N n=1 exp(-ln • ζ n (j)) (11)
The backward step ĵ is then defined by ĵ = arg min j Ej [START_REF] Zhao | Stagewise lasso[END_REF] When the algorithm proceeds and selects a coordinate g at the iteration t, the stored values will be altered as follows. If g is selected for the first time, then

∀n ζ n (g) = k =g ζn(k) ; ζn(g) = ε • hg(In) (13)
and

ζ (t) n (j) = ζ (t-1) n (j) + ε • hg(In) ∀j = g (14)
Now, if g already belongs to the active index IA, then

ζ (t) n (g) = ζ (t-1) n (g) + sign(ζ (t-1) n (g)) • ε (15)
In this case, equation ( 14) is valid for all j ∈ IA. Note that when j = g, this equation automatically takes into account the backward step (i.e the term -ε • 1g) because it was not added in the first place.

Prediction

For classification purpose, each image is predicted separately by looping through all weak classifiers. However, to be more efficient in retrieval, we predict by means of range queries. The final scores of all images are computed at the same time. First, these scores are initialized to zero. Then, and for each weak classifier, we query the search engine to retrieve all the images that fall into its hypersphere. Consequently, these images will see their scores incremented each by the weight of this weak classifier. On the other hand, images that are not returned by the system will have their scores decremented by the same weight of this classifier. After looping through all weak classifiers, each image ends up with a score judging its relevance to the concept.

The process described earlier can be relatively slow depending on the database size and the number of hypotheses constituting the model. For applications that have a fixed image database (not updated online), we can compute a priori the distances separating each weak classifier from each image and load them when the search engine starts. This could be achieved with no worries because the models used are concise and so won't consume too much computer memory. Moreover, since distances are computed offline, the process can be done either exhaustively or using an index structure. Unlike retrieval with range queries, to answer a given query, images are processed one by one. For each image, the distances to the model defined are obtained through a lookup table then the score is computed by comparing these distances to the corresponding classifier radii.

EXPERIMENTS

Performance Evaluation

Our experiments aims to prove that sparse models are still able to achieve good retrieval results. We first evaluated the retrieval effectiveness on 10 object categories belonging to the ImageNet database 1 . It consists on comparing the performance between BLasso and AdaBoost. The name of the categories as well as the composition of the database used for training are given in table 1. Note that for each object category, we built the counter-class (negative samples) by randomly collecting images from the other categories while keeping the same number. We relied on SIFT descriptors as the basic image primitives. Moreover, each image contained at most 500 features. The original image set was divided equally between training and test. Therefore, the test database has the same composition shown in table 1 (a total of 6830 images). For all our experiments, the parameters of BLasso were tuned as follows ξ = 10 -6 and ε = 1 80

This choice is made based on the experiments in [START_REF] Zhao | Stagewise lasso[END_REF]. For AdaBoost, we fixed the number of iterations T = 100.

For evaluation, we used the average precision. It is a reliable measure for retrieval and is computed as the area under the precision-recall curve. Precision is defined as the ratio of the number of correct answers to the number of the documents retrieved. Recall is defined as the ratio of the number of correct answers to the number of all the relevant documents in the database. Results are given in table BLasso outperformed AdaBoost in 6 categories (camel, laptop, penguin, sunflower, tomato and zebra). However, on average, AdaBoost was better approximately by 1.6%. This result is due to the clear domination of AdaBoost in the categories revolver and tennis_racket where the performance was respectively higher by 11.3% and 8.7%. On the other hand, it is noticeable that the feature selection performed by BLasso is better. In fact, the number of the features composing the model is fewer. Knowing that the prediction time is linear with the model size, the time needed by AdaBoost is 1.5 times the time needed by BLasso. Furthermore, and in general, AdaBoost get stuck at the stop condition defined (T = 100) and didn't stop earlier except for revolver. This is problematic for AdaBoost because it proves that the solutions found are not statistically optimal, they might be visually correlated due to an excessive selection or they might need other features to be added in order to well cover the intra-class variability. Precision-recall curves are given in figure 1. They were averaged over the 10 categories. As proven numerically, the average performance of AdaBoost is slightly higher. However, we notice that for the first results returned (recall of 5%), the precision given by BLasso is similar or slightly better.

The second experiment aims to measure to which extent the models can be generic. In fact, the selection of Ima-geNet categories was done in purpose to be able to predict with the models generated on the same concepts belonging to Caltech2562 database. Our test database was constructed from 15, 360 images from all the 256 Caltech categories with only 586 images relevant to the 10 predicted categories. Details on the composition of these images are given in table 3. In this experiment, we compare the performance of the models trained with an image set collected from Caltech to the models trained with an image set collected from ImageNet. Note that the number of the training images used in the Caltech model is the same used in prediction (cf. table 3). Results of this experiment are presented in table 4. We notice that the models belonging to ImageNet outperformed the models of Caltech in 6 categories and that they also did better, on average. This proves that the models generated are generic in terms that, statistically, they don't overfit. Therefore, they can be applied to various test databases.

Interactive Retrieval

How It Works

Interactive search is based on the user specialization of the models. Thanks to interpretability, the trained patches allow users to retrieve images, with not just the object inside, but also within a specific context, or with a desired scale or pose. For example, one may need to search for images containing cars in sand roads or airplanes in airports, etc.

Figure 2 is a snapshot of the user interface we developed. It gives an idea about what our interactive retrieval looks like. On the left side, users can choose a category from a list of predefined visual categories. The corresponding model is then loaded. Consequently, users can browse the visual patches displayed and start to form their own concept. In figure 2, the model of the sunflower category is chosen. In the left panel, we can browse all the 101 patches selected during training (cf. table 2). Note that the patches are displayed with their relative weights learned during training. However, when building their models, users have the possibility to change these confidence measures according to their own perception. That is, the more relevant a feature is, the higher the weight it will be assigned. Each visual patch chosen as well as its weight are displayed in the middle panel (11 patches in our snapshot). The visual query is now ready and the search engine can be interrogated. The retrieved images are displayed in the right panel.

Interactively Building Concepts

In this section, we will give visual illustrations on how interactivity is beneficial. The first example is given in figure 3. The user-specialized model used for retrieval (shown at the top) comprises six patches. The choice to make for displaying patches in a way that truly takes into account the description involved is still challenging. Here, we displayed the patches in gray-scale just to point out that the description is not color-related. They were taken from the sunflower model and they represent a global view of sunflowers within their context. Next, in fig. (3-2), we show the first page search results. As we can notice, the retrieved images match the query. They mostly contain sunflowers within a context (field, vegetation, sky). This result is due to the scale information brought by the patches. On the other hand, notice the presence of an image containing tomatoes. Even though it is not appreciated, such result explains the descriptor limitations. The image indeed looks like a sunflower. Now, have a look to the example shown in figure 2 (3-1) A specialized sunflower model used for retrieval. and notice the difference between first ranked results. The model used is displayed in the middle panel. Unlike the previous example, the images retrieved here tend to occupy all the image area.

As a second example, we present two different queries on the zebra category and their respective results (cf. fig. 4 and5). Both queries are related to the zebra category. The first query focus on the zebra upper-front part (head) while the second query is rather a general and global view of a zebra.

Our last experiment for interactivity is an attempt to build a hierarchical search engine that can semantically provide models for top level concepts and their sub-concepts. The idea is to draw together the categories camel, penguin, snail and zebra under the concept animal, then the categories laptop, revolver, tennis_racket and watch under the concept man-made and finally the categories sunflower and tomato under the concept vegetation. For training purpose, we used 30 images for each sub-concept. In summary, there were 120 images used in animal and man-made concepts and only 60 images used for vegetation. For prediction, we used a total of 13, 363 images. AP results are given in table 5. They are rather promising and demonstrate the feasibility of a hierarchical search engine with our method. (5-1) Zebra specialized model used for retrieval.

Category

(5-2) Search results.

Limitations

In our experiments, we noticed that it is possible for a single patch not to perform well if used by itself even if it is interpretable enough for a human. Two explanations can justify this limitation. The first possible cause is that our visual representation of the patch (gray-scale image corresponding to the description window of an interest point) doesn't faithfully reproduce the information the descriptor characterized. That is, the information we perceive from the patch is more complete than the semantic information coded by SIFT. The second explanation is that, statistically speaking, the database may be somehow skewed so that a single feature cannot be very discriminative. This problem is illustrated in figure 6 by the revolver category. The full model only comprises two patches that are shown at the top of the figure and their corresponding precision-recall curves are given below. Compare precision for the first ranked results (say 5% recall), we see that the full model behave quite well while each feature, alone, doesn't.

CONCLUSIONS

We introduced a new way to search for concepts in heterogeneous image databases through sparse and interpretable models. The models are a weighted set of visual patches that are directly mapped to image local features. For a given category, they represent the best features that summarize the intra-class variability. Consequently, we use them to provide a start page allowing users to build specialized models that, somehow, represent the concepts in their mind. Furthermore, we gave an efficient way to implement our learning strategy with the BLasso algorithm. Experiments revealed that our method gives equivalent performance than AdaBoost while consistently reducing the model complexity, giving better interpretability and saving the prediction time. Qualitatively, we showed how the interaction with the model can improve retrieved results.

In our future work, we plan to quantitatively evaluate user interactivity by generating a ground truth for more complex semantic queries. Moreover, we will think about experiment-ing new types of constraints in order to emphasize on some local aspects in images and see the subsequent effects on the model composition. We also intend to enrich patch description by using more coherent and compound representation through feature sets. Another remaining challenge is how to visually display the patches in a way to faithfully reproduce the underlaid description.

2 .

 2 Train the classifier and get the best hypothesis h (t+1) κ 3. β(t+1) = β(t) + ε • 1κ

1Figure 1 :

 1 Figure 1: Precision-recall curves for BLasso and Ad-aBoost.

Table 3 :

 3 Composition of the training/test set.

( 3 - 2 )

 32 Search results.

Figure 3 :

 3 Figure 3: Retrieving sunflowers.

Figure 2 :( 4 - 2 )

 242 Figure 2: Example of an interactive search

Figure 4 :

 4 Figure 4: Example of patches to focus on the upperfront part.

Figure 5 :

 5 Figure 5: Example of patches to look for a global view of zebras.

Figure 6 :

 6 Figure 6: Performance of the revolver category.

Table 1 : Composition of the training set.

 1

	2.

Table 2 : Performance evaluation of BLasso (BL) and AdaBoost (AB) in 10 categories of ImageNet.

 2

Table 4 : Illustration of prediction on a different database.

 4

Table 5 : Results of top semantical concepts.

 5

		# features	AP
	animal	36	0.5516
	man-made	22	0.5536
	vegetation	33	0.4224

http://www.vision.caltech.edu/Image Datasets/Caltech256/

ACKNOWLEDGMENTS

This work was funded by the EU through the European Integrated Project GLOCAL http://www.glocal-project.eu/ and the Agropolis Foundation through the project Pl@ntNet http://www.plantnet-project.org/papyrus.php?langue=en.