Théophile Dalens

Josef Sivic

Ivan Laptev

Marine Campedel

Painting recognition from wearable cameras Project Report

Are smart glasses the new high-tech device that will guide you through a museum? In this report, we describe a system that runs "on device" on Google Glass and retrieves the painting you're looking at among a set of paintings in a database. We perform an experimental comparison of the accuracy and speed of different feature detectors and descriptors on a realistic dataset of paintings from Musée du Louvre. Based on this analysis we design an algorithm for fast and accurate image matching using bags of binary features. We develop an application that runs directly on Google Glass without sending images to external servers for processing and recognizes a query painting in a database of 100 paintings in one second. The code is available online at [1].

) uses our application in a museum to identify a painting (middle). Our application identifies the painting from the database of the museum and displays information such as the name of the painter and the title of the painting (right).

Introduction 1.Objective

Google Glass is a newly released wearable computer, providing to its user a camera and a display screen, both directly at the level of its eye. Computer vision has many 1 picture by Antonio Zugaldia, http://www.flickr.com/photos/azugaldia/7457645618 applications that would benefit the device and its user, including object detection, categorization and retrieval. For example, the device could help its user to identify and describe what he is currently seeing. However, wearable devices have limited processing power and available rapid access memory. Therefore, any computer vision-related application must be time and space efficient. In this work, we aim to develop a method for object retrieval that works in real-time on Google Glass. Our task is to make a painting recognition algorithm, that can match a painting acquired by the device in a museum and match it with the corresponding painting in the database of the museum. The whole application should work in realtime and display information about the painting, like the name of the painter. For this, we detail the general framework of such an algorithm (section 2), consisting of keypoint extraction and descriptors computation, followed if needed by grouping the descriptors into a single image descriptor, and/or geometric verification. In section 3, we compare different approaches for this task, in terms of efficiency and time performance. Then, in section 4, we detail our implementation and give qualitative results.

Related work

Using the Google Glass for computer vision has been studied in [START_REF] Altwaijry | Recognizing locations with google glass: A case study[END_REF], in which Altwaijry et al. design a location recognition on Google Glass. In this work, they have a database of images of different buildings and scenes of a city and aim to identify what the user is looking at. They compare different approaches for image description and matching, including the bag of words approach (as described in [START_REF] Sivic | Video google: A text retrieval approach to object matching in videos[END_REF]) with SIFT and SURF descriptors, which performed successfully. However, their approach is different than the one we try to develop in this work, for two reasons. First, they integrate GPS information, which we don't -our system is designed to work in a very restricted area such as a museum aisle. Secondly, data from the Google Glass is sent to a server for computation. While an online server may have more computational power than the Google Glass, this prevents the system from working without connection to the server and complicates the deployment of the system. Our system is designed to be easy to deploy and work offline, and must face the computational limitations of the wearable device.

Overview of the retrieval system

Our system is composed of two distinct parts. First, in the offline part, we compute image descriptors of all the images in our database. This part is done on a computer, and the results are then stored on the Glass. The on device part is the wearable device application. When the application starts, descriptors are loaded into the memory of the device. At the user's command, the camera on the device captures the query image, and computes the descriptors with the same algorithm as the one used for the offline part. Then, descriptors of the query image are matched to the ones in the database. A match between a descriptor of the query image and a descriptor of the database image is considered good if there is no other descriptor in the database image that matches nearly as well as the first one. The images in the database are then ranked decreasingly according to the number of good matches they have with the query image. Alternatively, descriptors are matched to visual words to compute a single image descriptor for each image, and we just compare the image descriptors. At last, a geometric verification can be performed.

Keypoint detection and descriptors

A capture of a painting by the visitor of the museum, and the painting itself, are different at various aspects. A geometric transformation is applied (translation and homothetic transformation, with also small rotation and perspective transformation), lightning conditions change and partial occlusions can occur. Hence, it is necessary for our task that image description is as robust as possible to those changes. We consider several possibilities. While SIFT is widely used and is still considered as one of the most robust descriptors, its speed is a concern for real-time applications, leading to the development of descriptors that try to match their performance while being more efficient to compute and to compare. Figure 2 illustrates the various keypoints that can be detected in a painting.

SIFT

Scale Invariant Feauture Transform ( [START_REF] David G Lowe | Object recognition from local scale-invariant features[END_REF]) is a classic algorithm to detect keypoints and compute associated descriptors. First, the input image is resampled and smoothed at different scales to create a pyramid of images. Then, keypoints are selected as extrema of differences of Gaussian. Descriptors are then a 128-bin histogram of the differences of Gaussian at different spaces and orientation.

SURF

Speeded-Up Robust Features ( [START_REF] Bay | Speeded-up robust features (surf)[END_REF]), designed as a faster alternative to SIFT, computes an approximation of the Hessian matrix in each pixel of the image at different scales.

Keypoints are selected as the maxima of the Hessian matrix. Descriptors, along with an orientation assignment, are computed as the sums of the responses of patches to wavelets, which result in 64D vectors.

ORB

ORB ( [START_REF] Rublee | Orb: an efficient alternative to sift or surf[END_REF]) is short for Oriented FAST and Rotated BRIEF. The keypoint detector uses FAST ( [START_REF] Rosten | Faster and better: A machine learning approach to corner detection[END_REF]), a corner detector, combined with a Harris corner detector ( [START_REF] Harris | A combined corner and edge detector[END_REF]), at multiple scales. It also measures a corner orientation. Then, BRIEF ( [START_REF] Calonder | Brief: Binary robust independent elementary features[END_REF]) descriptors are computed. A BRIEF descriptors is a bit string of length n, where each bit i equals the result of the test I(x i ) < I(y i ), I(a) being the intensity of the image I at pixel location a. In practice, the image is smoothed during this operation: locations x i and y i are chosen as 5 × 5 patches in a 31 × 31 window.

BRISK

Binray Robust Invariant Scalable Keypoints ( [START_REF] Leutenegger | Brisk: Binary robust invariant scalable keypoints[END_REF]) is another descriptor that uses FAST keypoints and descriptors that are similar to BRIEF. The main difference is in the sampling part that is used to generate the locations x i and y i .

Matching

To compare SIFT or SURF descriptors, the L 2 norm is used. Matching binary descriptors is faster, since it uses Hamming distance, which is the number of bits that differ between the two descriptors. To search for the nearest neighbour of a "query" descriptor in a set of "database" descriptors, we can use brute-force matching, i.e. matching the query descriptor with all the database descriptors and see which gives the highest L 2 norm or Hamming distance. However, this method takes time and can be the bottleneck of our system is the set of database descriptors is large. In [START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF] and [START_REF] Muja | Fast approximate nearest neighbors with automatic algorithm configuration[END_REF], the authors develop forests of randomized k-d trees in order to retrieve the approximate nearest neighbour of a SIFT descriptor, reducing the complexity of finding the nearest neighbour in a database of size K from O(K) to O(log(K)). For binary descriptors, a classic way to speed up matching is to use locality sensitive hashing. Descriptors are hashed into smaller hashes, which are faster to compare. Then, the real descriptor comparison is only done between the ones that have similar hashes.

Compact image descriptor

Matching the descriptors of the query image to all the descriptors of the database image is not a scalable strategy. Indeed, its complexity is linear with the number of paintings in the database. To solve this issue, the authors of [START_REF] Sivic | Video google: A text retrieval approach to object matching in videos[END_REF] propose a "bag of features", or "bag of words", approach. It consists in building a vocabulary of features, then in representing each image by a histogram. For each image, descriptors are converted into words by taking the approximate nearest neighbour (which is fast with the previous subsection), then the histogram counts the number of occurrences of each word. Histograms are then normalized and re-weighted using tf-idf so that rare words are considered more significant then frequent ones. To match images, we need only match one histogram instead of hundreds of descriptors.

To build the dictionary, the authors of [START_REF] Sivic | Video google: A text retrieval approach to object matching in videos[END_REF] propose to run k-means algorithm to features extracted from various images, if possible including the images of the database. The k centers are then the k words of the vocabulary. For binary descriptors, the authors of [START_REF] Galvez-Lopez | Bags of binary words for fast place recognition in image sequences[END_REF] replace k-means by k-medians, after "k-means ++" initialization. k-means ++ ( [START_REF] Arthur | K-means++: The advantages of careful seeding[END_REF]) is an algorithm for seeding k-means that select iteratively k initial cluster centers. The first center is chosen according to a uniform distribution. After that, at each step, the distance between each point and its nearest center is computed and squared; the next center is taken randomly according to this distribution. As a result, cluster centers chosen by this algorithm are close to the expected result of the k-means algorithm.

Geometric verification

Once we have potential good matches, we can make use of the spatial repartition of the descriptors in order to retrieve an image that is geometrically consistent. If a query image contains an image of our database, then there exists a geometric transformation, a homography, that matches part of the descriptors of the query image to the database image. To find such a homography, we can use the RANSAC algorithm ( [START_REF] Martin | Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography[END_REF]), which gives us both the best homography to fit the model and the number of inliers (descriptors that fit the homography). The geometric verification consists in applying this algorithm between the query image and all the best results in our database, then stopping at the first one that gives an important number of inliers. However, each geometric verification takes time to perform.

Evaluation of feature extraction and matching

Database and set up

To evaluate our system, we use as a database a set of paintings from Musée du Louvre, collected from a CD-ROM guide of the museum. It includes more than 600 paintings. The resolution of the images is not high (typically 250x350). From these paintings, we have selected a representative set of 100, that include portraits, landscapes and other sceneries. The styles of the paintings are loosely consistent (mainly Italian and French paintings), as to maintain a certain difficulty and to simulate a typical set of paintings in a aisle of the museum. A set of paintings is shown figure 3. As query images, we use data collected from the Google Glass in the Musée du Louvre. We went to the room 77, also called salle Mollien, and to the aisles with the French paintings, from which we casually captured 10 images that include landscapes, portraits and other sceneries. Each picture has been rescaled using the same factors. They include the typical small occlusions and perspective deformations that can be expected when a visitor looks at paintings. They are shown in figure 4. For all descriptors, we use the OpenCV [5] implementation.

Accuracy

To measure the accuracy of the different descriptors, we match the descriptors of computed on our input images and the descriptors of the images of our database. Then, we rank the database images decreasingly according to the number of matches. We measure recall at top 1 (number of times that the correct database painting arrives in first position) and recall at top 5 (number of times that the correct painting arrives in at worst at fifth position). The recall at top 5 is interesting because if we complete the matching with geometric verification, having the correct answer in the top 5 will assure that we only need to perform this costly operation on at most 5 images. We have tested the different keypoint detectors and descriptors with various parameters and only display the ones that give the best results (default parameters for SIFT and BRISK, Hessian threshold of 475 for SURF, 1000 keypoints for ORB). We also display the results for ORB with only 500 keypoints, as decreasing the number of keypoints decreases the computational time for matching. At last, we have found that using SURF keypoints gave better results than the BRISK ones, when using BRISK descriptors. Results are shown table 1. While all the systems perform similarly for the landscapes, their performance gener-

System

Recall at top ally decrease for the portraits. ORB performs well every time when it computes enough descriptors; otherwise, the keypoints captured on the images are not dense enough.

For the ORB descriptors, we have also evaluated the performance of a bag of features approach. For these tests, we have computed dictionaries of 10,000 and 100,000 words from the extraction of ORB features from each of 600 paintings, including the 100 ones used for testing. The performances are competitive with the previous approach (table 2).

Timing

Descriptor extraction

We have measured, on the Glass, the time taken to detect keypoints and compute their descriptors. Comparisons are made on the same input image. Results are displayed in table 3. As expected, binary features are faster to compute than SURF and SIFT. Time used System Recall at top 1 Recall at top 5 Bag of ORB (10,000 words) 80% 80% Bag of ORB (100,000 words) 90% 90% to detect BRISK is unexpectedly high, and most probably comes from a flaw in its OpenCV implementation.

Descriptor matching and hierarchical bag of words

We have also timed binary feature matching. Matching 500 ORB features of a query image with 500 ORB features of a database image takes approximately 400 milliseconds.

When matching 500 descriptors to a vocabulary of 100,000 words, this matching becomes the bottleneck of our application. For this reason, we use a hierarchical bag of words, similar the hierarchical k-means described in [START_REF] Nister | Scalable recognition with a vocabulary tree[END_REF] and [START_REF] Muja | Fast approximate nearest neighbors with automatic algorithm configuration[END_REF], but applied here to binary features, and building the tree of features from leaves to roots. From the 100,000 words that have been extracted as a vocabulary, we extract a smaller number of words (typically 500). This new, small vocabulary constitutes another layer in our bag of visual words. Then, we assign each word of the large vocabulary to its k nearest neighbours in the small one (typically k = 2). To get a bag of words representation of the descriptors of our image, we match the descriptors of our database with the words of the small vocabulary. Then, we match them with the words of the large vocabulary that are assigned to matched words of the small one. As a result, the complexity in time of the matching part is in the square root of the size of the vocabulary, instead of being linear. Typically, the 500 descriptors extracted are matched one-to-one with the words of the small vocabulary (500), then with the linked words of the large one (in average 400, up to 2400), which results in less than 1 second of computation. With these parameters, we do not observe drop in accuracy. However, reducing k to 1, which would allow a faster matching, reduces the accuracy significantly.

Implementation on Google Glass

The Google Glass is a device running on Android. To implement the on device part on the algorithm, we use Java, C++ (via JNI) and OpenCV. The code described in this report can be found online at [START_REF]GlassPainting code repository[END_REF]. Inspired from the results of section 3, in which we show the efficiency of the ORB descriptors in accuracy and speed, we propose two implementations of our system.

Implementation using direct feature to feature matching

In the first version, the offline part consists in computing ORB features for all our database. These features are saved into XML files, which are then transferred to the device at the installation of the application. Then, at the start of the application, these descriptors are loaded. When the user taps the touchpad of the device, an image of what the user is currently seeing is captured. Then, the extraction of the ORB features, matching and ranking of images is performed in C++ and native OpenCV code. The result is handled by the Java part of the application, that has loaded a table with the names of the painters of each painting. Then name of the painter is then displayed.

Implementation using bag-of-visual-words matching

The second version uses bags of ORB features. In the offline part, we compute a vocabulary of 100,000 ORB words. To compute the vocabulary, we extract ORB features from the paintings of the database, and perform k-means++ seeding. We then compute the histograms of words for each painting in the database. These histograms are stored into a single sparse matrix. Using a sparse matrix is essential here, to limit both the space and time resources that will be used by the application when loading the matrix and performing comparisons of histograms. Still in the offline part, we compute a second layer of 500 visual words. These are selected via k-means++ from the large vocabulary. Then, each of the words of the large vocabulary is assigned to its 2 nearest neighbours in the small one. These assignments are also saved and stored to the device. At the start of the Glass application, the histogram, the small and the large vocabulary and the assignments are loaded. When the user touches the touchpad, 500 ORB features are extracted from the image captured by the camera, and transformed into words using brute-force matching with the small vocabulary. Then, we match each feature to a word of the large vocabulary whose nearest neighbours include the word of the small vocabulary, using the assignments matrix. The resulting histogram is compared to the ones of the database by performing a sparse matrix multiplication. The name of the painter of the image that has the most similar histogram is then retrieved and displayed at the Google Glass screen. The whole process takes approximately 1.0 second.

Results

While the first version performs successfully most of the time, it is not scalable to many images. With 10 images in the database, it runs in approximately 5 seconds, and the time is expected to be linear with the number of images in the database. Indeed, matching features from one image to another takes 0.4 seconds. The second version, on the other hand, performs in 1 second for 100 images. However, it is more slightly more likely to fail and not retrieve the correct painting. A sample of the results is shown figure 5.

Figure 5: Examples of correct and incorrect retrievals of the second method, which is expected to retrieve the correct painting nine times out of ten.

Conclusion

In this work, we have compared different descriptors in order to perform object retrieval, in the context of paintings in a museum. We have successfully implemented on Google Glass a system able to tell a user which painting he is looking at, among the different paintings of a section of a museum. It uses a hierarchical bag of ORB features, which are fast to extract and efficient in matching.

Figure 1 :

 1 Figure 1: Our application. A Google Glass user (left 1 ) uses our application in a museum to identify a painting (middle). Our application identifies the painting from the database of the museum and displays information such as the name of the painter and the title of the painting (right).

Figure 2 :

 2 Figure 2: Detected keypoints on a painting with various keypoint detectors, all tuned to detect approximately 500 keypoints: SIFT (top left), SURF (top right), ORB (botttom left) and BRISK (bottom right).

Figure 3 :

 3 Figure 3: Sample of the 100 paintings in the database used for our tests.

Figure 4 :

 4 Figure 4: Input images used for tests. These paintings are from (top to bottom, then left to right) Corot (3 paintings), David, Lebrun (2 paintings), Delacroix, Géricault, Benoist and Ingres

Table 1 :

 1 Performance measurements for the different descriptors. "SURF+BRISK" designs a combination of the SURF keypoint detector and the BRISK feature descriptor.

	1 Recall at top 5

Table 2 :

 2 Accuracy measurements for the bag of binary features (ORBs). The number of words shown is the size of the vocabulary.

	Description method Time for keypoint detection (ms) Time for computation (ms)
	SIFT	4180	3,86
	SURF	2008	2,47
	ORB	109	0,39
	BRISK	8106	0,15

Table 3 :

 3 Time to detect keypoints on a 640x480 picture on the Google Glass, and to compute a single descriptor.

Acknowledgements

This work was partially funded by the EIT-ICT labs. The authors thank Jia Li from Télécom ParisTech for sharing useful code.


Damien Garreau
email: damien.garreau@ens.fr

Rémi Lajugie
email: remi.lajugie@inria.fr

Sylvain Arlot
email: sylvain.arlot@ens.fr

Francis Bach
email: francis.bach@inria.fr

Metric Learning for Temporal Sequence Alignment

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

The problem of aligning temporal sequences is ubiquitous in applications ranging from bioinformatics [START_REF] Cuturi | A kernel for time series based on global alignments[END_REF][START_REF] Aach | Aligning gene expression time series with time warping algorithms[END_REF][START_REF] Thompson | Balibase: a benchmark alignment database for the evaluation of multiple alignment programs[END_REF] to audio processing [START_REF] Cont | Evaluation of real-time audioto-score alignment[END_REF][START_REF] Dixon | Match: A music alignment tool chest[END_REF]. The idea is to align two similar time series that have the same global structure but local temporal differences. Most alignments algorithms rely on similarity measures, and having a good metric is crucial, especially in the high-dimensional setting where some features of the signals can be irrelevant to the alignment task. The goal of this paper is to show how to learn this similarity measure from annotated examples in order to improve the precision of the alignments. For example, in the context of music information retrieval, alignment is used in two different cases: (1) audio-to-audio alignment and (2) audio-to-score alignment. In the first case, the goal is to match two audio interpretations of the same piece that are potentially different in rythm, whereas audio-to-score alignment focuses on matching an audio signal to a symbolic representation of the score. In the second case, there are some attempts to learn from annotated data a measure for performing the alignment. Joder et al. [START_REF] Joder | Learning optimal features for polyphonic audio-to-score alignment[END_REF] propose to fit a generative model in that context, and Keshet et al. [START_REF] Keshet | A large margin algorithm for speech-to-phoneme and music-to-score alignment[END_REF] learn this measure in a discriminative setting. Similarly to Keshet et al. [START_REF] Keshet | A large margin algorithm for speech-to-phoneme and music-to-score alignment[END_REF], we use a discriminative loss to learn the measure, but our work focuses on audio-to-audio alignment. In that context, the set of authorized alignments is much larger, and we explicitly cast the problem as a structured prediction task, that we solve using off-the-shelf stochastic optimization techniques [START_REF] Lacoste-Julien | Block-coordinate Frank-Wolfe optimization for structural SVMs[END_REF] but with proper and significant adjustements, in particular in terms of losses. The need for metric learning goes far beyond unsupervised partitioning problems. Weinberger and Saul [START_REF] Weinberger | Distance metric learning for large margin nearest neighbor classification[END_REF] proposed a large-margin framework for learning a metric in nearest-neighbour algorithms based on sets of must-link/must-not-link constraints. Lajugie et al. [START_REF] Lajugie | Large-margin metric learning for constrained partitioning problems[END_REF] proposed to use a large margin framework to learn a Mahalanobis metric in the context of partitioning problems. Since structured SVM have been proposed by Tsochantaridis et al. [START_REF] Tsochantaridis | Large margin methods for structured and interdependent output variables[END_REF], Taskar et al. [START_REF] Taskar | Max-margin Markov networks[END_REF], they have successfully been used to solve many learning problems, for instance to learn weights for graph matching [START_REF] Caetano | Learning graph matching[END_REF] or a metric for ranking tasks [START_REF] Mcfee | Metric learning to rank[END_REF]. They have also been used to learn graph structures using graph cuts [START_REF] Szummer | Learning crfs using graph cuts[END_REF].

Contributions. We make the following five contributions: (1) we cast the learning of a Mahalanobis metric in the context of alignment as a structured prediction problem, (2) we show that on real musical datasets this metric improves the performance of alignment algorithms using high-level features, (3) we propose to use the metric learning framework to learn combinations of basic audio features and get good alignment performances, (4) we show experimentally that the standard Hamming loss, although tractable computationnally does not permit to learn a relevant similarity measure in some real world settings, [START_REF] Cont | Evaluation of real-time audioto-score alignment[END_REF] we propose a new loss, closer to the true evaluation loss for alignments, leading to a tractable learning task, and derive an efficient Frank-Wolfe based algorithm to deal with this new loss.

Matricial formulation of alignment problems 2.1 Notations

In this paper, we consider the alignment problem between two multivariate time series sharing the same dimension p, but possibly of different lengths T A and T B , namely A ∈ R T A ×p and B ∈ R T B ×p . We refer to the rows of A as a 1 , . . . , a T A ∈ R p and those of B as b 1 , . . . , b T B ∈ R p . From now on, we denote by X the pair of signals (A, B). Let C(X) ∈ R T A ×T B be an arbitrary pairwise affinity matrix associated to the pair X, that is, C(X) i,j encodes the affinity between a i and b j . Note our framework can be extended to the case where A and B are multivariate signals of different dimensions, as long as C(X) is well-defined. The goal of the alignment task is to find two nondecreasing sequences of indices α and β of same length u ≥ max(T A , T B ) and to match each time index α(i) in time series A to the time index β(i) in the time series B, in such a way that u i=1 C(X) α(i),β(i) is maximal, and that (α, β) satisfies:

     α(1) = β(1) = 1 (matching beginning) α(u) = T A , β(u) = T B
(matching ending) ∀i, (α(i + 1), β(i + 1)) -(α(i), β(i)) ∈ {(1, 0), (0, 1), (1, 1)} (three type of moves)

(1) For a given (α, β), we define the binary matrix Y ∈ {0, 1} T A ×T B such that Y α(i),β(i) = 1 for every i ∈ {1, . . . , u} and 0 otherwise. We denote by Y(X) the set of such matrices, which is uniquely determined by T A and T B . An example is given in Fig. 1. A vertical move in the Y matrix means that the signal B is waiting for A and an horizontal one that A is waiting for B. In this sense we can say that the time reference is "warped". When C(X) is known, the alignment task can be cast as the following linear program (LP) over the set Y(X): max

Y ∈Y(X) Tr(C(X) ⊤ Y ). ( 2
)
Our goal is to learn how to form the affinity matrix: once we have learned C(X), the alignment is obtained from Eq. ( 2). The optimization problem in Eq. (2) will be referred to as the decoding of our model.

Dynamic time warping.

Given the affinity matrix C(X) associated with the pair of signals X = (A, B), finding the alignment that solves the LP of Eq. ( 2) can be done efficiently in O(T A T B ) using a dynamic programming algorithm. It is often referred to as dynamic time warping [START_REF] Cuturi | A kernel for time series based on global alignments[END_REF][START_REF] Müller | Information retrieval for music and motion[END_REF]. This algorithm is described in Alg. 1 of the supplementary material. Various additional constraints may be used in the dynamic time warping algorithm [START_REF] Müller | Information retrieval for music and motion[END_REF], which we could easily add. The cardinality of the set Y(X) is huge: it corresponds to the number of paths on a rectangular grid from the southwest (1, 1) to the northeast corner (T A , T B ) with only vertical, horizontal and diagonal moves allowed. This is the definition of the Delannoy numbers [START_REF] Banderier | Why Delannoy Journal of statistical planning and inference[END_REF]. As noted in [START_REF] Torres | An exact formula for the number of alignments between two dna sequences[END_REF], when

t = T A = T B is big, one can show that #Y t,s ∼ (3+2 √ 2) t √ πt √ 3 √ 2-4
.

The Mahalanobis metric

In many applications, for a pair X = (A, B), the affinity matrix is computed by

C(A, B) i,j = -a i,k -b j,k 2 .
In this paper we propose to learn the metric to compare Red upper triangles show the (i, j) such that Y1 i,j = 1, and the blue lower ones show the (i, j) such that Y 2 i,j = 1. The grey zone corresponds to the area loss δ abs between Y 1 and Y 2 , whereas the δ max loss corresponds to the maximum of the δ t , t ∈ {1, . . . , T 1 }. a i and b j instead of using the plain Euclidean metric. That is, C(X) is parametrized by a matrix W ∈ W ⊂ R p×p , where W ⊂ R p×p is the set of semi-definite positive matrices, and we use the corresponding Mahalanobis metric to compute the pairwise affinity between a i and b j :

C(X; W ) i,j = -(a i -b j ) ⊤ W (a i -b j ).
Note that the decoding of Eq. ( 2) is the maximization of a linear function in the parameter W :

max Y ∈Y(X) Tr(C(X; W ) ⊤ Y ) ⇔ max Y ∈Y(X) Tr(W ⊤ φ(X, Y )), (3)
if we define the joint feature map φ(X, Y

) = -i,j Y i,j (a i -b j )(a i -b j ) ⊤ ∈ R p×p .

Learning the metric

From now on, we assume that we are given n pairs of training instances 1

(X i , Y i ) = ((A i , B i ), Y i ) ∈ R T i A ×p ×R T i B ×p ×R T i A ×T i B , i = 1, . . . , n.
Our goal is to find a matrix W such that the predicted alignments are close to the groundtruth on these examples, as well as on unseen examples. We first define a loss between alignments, in order to quantify this proximity between alignments.

Losses between alignments

In our framework, the alignments are encoded by matrices in Y(X), thus we are interested in functions ℓ : Y(X) × Y(X) → R + . Let us define the Frobenius norm by

M 2 F = i,j M 2 i,j .
Hamming loss. A simple loss between matrices is the Frobenius norm of their difference, which turns out to be the unnormalized Hamming loss [START_REF] Hamming | Error detecting and error correcting codes[END_REF] for 0/1 valued matrices. For two matrices Y 1 , Y 2 ∈ Y(X), it is defined as:

ℓ H (Y 1 , Y 2 ) = Y 1 -Y 2 2 F (4) = Tr(Y ⊤ 1 Y 1 ) + Tr(Y ⊤ 2 Y 2 ) -2 Tr(Y ⊤ 1 Y 2 ) = Tr(Y 1 1 T B 1 ⊤ T A ) + Tr(Y 2 1 T B 1 ⊤ T A ) -2 Tr(Y ⊤ 1 Y 2 ),
where 1 T is the vector of R T with all coordinates equal to 1. The last line of Eq. ( 4) comes from the fact that the Y i have 0 -1 values; that makes the Hamming loss linear in Y 1 and Y 2 . This loss is often used in other structured prediction tasks [START_REF] Lacoste-Julien | Block-coordinate Frank-Wolfe optimization for structural SVMs[END_REF]. In the audio-to-score setting, Keshet et al. [START_REF] Keshet | A large margin algorithm for speech-to-phoneme and music-to-score alignment[END_REF] use a modified version of this loss, which is the average number of times the difference between the two alignments is greater than a fixed threshold. This loss is easy to optimize since it is linear in our parametrization of the alignement problem, but not optimal for audio-to-audio alignment. Indeed, a major drawback of the Hamming loss is, for alignments of fixed length, it depends only on the number of "crossings" between alignment paths: one can easily find

Y 1 , Y 2 , Y 3 such that ℓ H (Y 2 , Y 1 ) = ℓ H (Y 3 , Y 1 ) but Y 2 is intuitively much closer to Y 1 than Y 3 (see Fig. 2).
It is important to notice this is often the case as the length of the signals grows.

Area loss. A more natural loss can be computed as the mean distance beween the paths depicted by two matrices Y 1 , Y 2 ∈ Y(X). This loss is represented by the grey zone on Fig. 1 and corresponds to the area between the paths of two matrices Y . Formally, as in Fig. 1, for each t ∈ {1, . . . , T B } we define δ t as the minimum between

| min{k, Y 1 k,t = 1} -max{k, Y 2 k,t = 1}| and | max{k, Y 1 k,t = 1} -min{k, Y 2 k,t = 1}|.
Then the area loss is the mean of the δ t . In the audio literature [START_REF] Kirchhoff | Evaluation of features for audio-to-audio alignment[END_REF], this loss is sometimes called the "mean absolute deviation" loss and is noted δ abs (Y 1 , Y 2 ). Unfortunately, in the general case of alignment problem δ abs is not linear in the matrices Y . But in the context of alignment of sequences of two different nature, one of the signal is a reference and thus the index sequence α defined in Eq. ( 1) is increasing, e.g. for the audio to partition alignment problem [START_REF] Joder | Learning optimal features for polyphonic audio-to-score alignment[END_REF]. This loss is then linear in each of its arguments. More precisely, if we introduce the matrices L ∈ R T A ×T A which is lower triangular with ones, we can write the loss as

ℓ O = L(Y 1 -Y 2 ) 2 F (5) = Tr(LY 1 1 T B 1 ⊤ T A ) + Tr(LY 2 1 T B 1 ⊤ T A ) -2 Tr(LY 1 Y ⊤ 2 L ⊤ ).
We now prove that this loss corresponds to the area loss in this special case. Let Y be an alignment matrix and L ∈ R T A ×T A be the matrix such that L r,s = 1 if and only if

r ≥ s. Then it is easy see that (LY ) i,j = k L i,k Y k,j = i k=0 Y k,j .
If Y does not have vertical moves, i.e. for each j there is an unique k j such that Y k,j = 1, we have that (LY ) i,j = 1 if and only if i ≥ k j . So i,j (LY ) i,j = #{(i, j), i ≥ k j }, which is exactly the area under the curve determined by the path of Y . In all our experiments, we use δ abs for evaluation but not for training.

Approximation of the area loss: the symmetrized area loss. In many real world applications [START_REF] Kirchhoff | Evaluation of features for audio-to-audio alignment[END_REF], the best loss to assess the quality of an alignment is the area loss. As shown by our experiments, if the Hamming loss is sufficient in some simple situations and allows to learn a metric that leads to good alignment performance in terms of area loss, on more challenging datasets it does not work at all (see Sec. 5). This is due to the fact that two alignments that are very close in terms of area loss can suffer a big Hamming loss. In Fig. 2, we provide examples where the Hamming loss is not sufficient to assess performance. Thus it is natural to extend the formulation of Eq. ( 5) to matrices in Y(X). We first start by symmetrizing the formulation of Eq. ( 5) to overcome problems of overpenalization of vertical vs. horizontal moves. Let L 1 ∈ R T B ×T B be the matrix such that L r,s = 1 if and only if r ≥ s. We define, for any binary matrices Y 1 and Y 2 ,

ℓ S (Y 1 , Y 2 ) = 1 2 L(Y 1 -Y 2 ) F + (Y 1 -Y 2 )L 1 ) 2 F (6) = 1 2 Tr(Y ⊤ 1 L ⊤ LY 1 ) + Tr(LY 2 1 T B 1 ⊤ T A ) -2 Tr(Y ⊤ 2 L ⊤ LY 1 ) + Tr(Y 1 L 1 L ⊤ 1 Y ) + Tr(Y ⊤ 2 1 T A 1 T B L 1 L ⊤ 1 Y 2 ) -2 Tr(Y 2 L 1 L ⊤ 1 Y ⊤ 1 .
We propose to use the following trick to obtain a concave loss over Y(X), the convex hull of Y(X). Let us introduce D = λ max (L ⊤ L)I T A ×T A and D 1 = λ max (L 1 L ⊤ 1 )I T B ×T B with λ max (U ) the largest eigenvalue of U . For any binary matrices Y 1 , Y 2 , we have that:

ℓ S (Y 1 , Y 2 ) = 1 2 Tr(Y ⊤ 1 (L ⊤ L -D)Y 1 ) + Tr(DY 1 1 T B 1 ⊤ T A ) + Tr(LY 2 1 T B 1 ⊤ T A ) -2 Tr(Y ⊤ 2 (L ⊤ L -D)Y 1 ) + Tr(Y 1 (L 1 L ⊤ 1 -D)Y ) + Tr(Y 1 D 1 1 T B 1 ⊤ T A ) Tr(Y ⊤ 2 L 1 L ⊤ 1 Y 2 ) -2 Tr(Y 2 L 1 L ⊤ 1 Y ⊤ 1 )
and we get a concave function over Y(X) that coincides with ℓ S on Y(X).

Empirical loss minimization

Recall that we are given n alignment examples (X i , Y i ) 1≤i≤n . For a fixed loss ℓ, our goal is now to solve the following minimization problem in W :

min W ∈W 1 n n i=1 ℓ Y i , argmax Y ∈Y T i A ,T i B Tr(C(X i ; W ) ⊤ Y ) + λΩ(W ) , (7)
where Ω = λ 2 W 2 F is a convex regularizer preventing from overfitting, with λ ≥ 0.

Large margin approach

In this section we describe a large margin approach to solve a surrogate to the problem in Eq. ( 7), which is untractable. As shown by Eq. ( 3), the decoding task is the maximum of a linear function in the parameter W and aims at predicting an output over a large and discrete space (the space of potential alignments with respect to the constraints in Eq. ( 1)). Learning W thus falls in the structured prediction framework [START_REF] Tsochantaridis | Large margin methods for structured and interdependent output variables[END_REF][START_REF] Taskar | Max-margin Markov networks[END_REF]. We define the hinge-loss, a convex surrogate to ℓ, by

L(X, Y ; W ) = max Y ′ ∈Y(X) ℓ(Y, Y ′ ) -Tr(W ⊤ [φ(X, Y ) -φ(X, Y ′ )]) . (8)
The evaluation of L is usually referred to as the "loss-augmented decoding", see [START_REF] Tsochantaridis | Large margin methods for structured and interdependent output variables[END_REF]. Among the aforementionned losses, the Hamming loss ℓ H is the only one leading directly to a tractable loss-augmented decoding problem and thus that falls directly into the structured prediction framework. Indeed, plugging Eq. ( 4) into (8) leads to a lossaugmented decoding that is a LP over the set Y(X) and that can therefore be solved using the dynamic time warping algorithm. If we define Y i as the argmax in Eq. ( 8) when (X, Y ) = (X i , Y i ), then elementary computations show that

Y i = argmin Y ∈Y T A ,T B Tr((U ⊤ -2Y i⊤ -C(X i ; W ) ⊤ )Y ),
where

U = 1 T B 1 ⊤ T B ∈ R T A ×T B .
We now aim at solving the following problem, sometimes called the margin-rescaled problem:

min W ∈W λ 2 W 2 F + 1 n n i=1 max Y ∈Y T i A ,T i B ℓ(Y, Y i ) -Tr(W ⊤ φ(X i , Y i ) -φ(X i , Y ) ) . (9)
Hamming loss case. From Eq. ( 3), one can notice that our joint feature map is linear in Y . Thus, if we take a loss that is linear in the first argument of ℓ, for instance the Hamming loss, the loss-augmented decoding is the maximization of a linear function over the spaces Y(X) that we can do efficiently using dynamic programming algorithms (see Sec. 2.1 and supplementary material). That way, plugging the Hamming loss (Eq. ( 4)) in Eq. ( 9) leads to a convex structured prediction problem. This problem can be solved using standard techniques such that cutting plane methods [START_REF] Joachims | Cutting-plane training of structural svms[END_REF], stochastic gradient descent [START_REF] Shalev-Shwartz | Pegasos: Primal estimated sub-gradient solver for SVM[END_REF], or block-coordinate Frank-Wolfe in the dual [START_REF] Lacoste-Julien | Block-coordinate Frank-Wolfe optimization for structural SVMs[END_REF]. Note that we adapted the standard unconstrained optimization methods to our setting, where W 0.

Optimization using the symmetrized area loss. In this section we propose to show that it is possible to deal The symmetrized area loss is concave in its first argument, thus the problem of Eq. 9 using it is in a min/max form and thus deriving a dual is straightforward. Details can be found in the supplementary material. If we plug the symmetrized area loss ℓ S (SAL) defined in Eq. ( 6) into our problem (9), we can show that the dual of ( 9) has the following form:

min (Z 1 ,...,Z n )∈Y 1 2λn 2 n i=1 -j,k (Y i -Z i ) j,k (a j -b k )(a j -b k ) T 2 F -1 n n i=1 ℓ S (Z, Z i ), (10)
if we denote by Y(X i ) the convex hull of the sets Y(X i ), and by Y the cartesian product over all the training examples i of such sets. Note that we recover a similar result as [START_REF] Lacoste-Julien | Block-coordinate Frank-Wolfe optimization for structural SVMs[END_REF]. Since the SAL loss is concave, the aforrementionned problem is convex. The problem ( 10) is a quadratic program over the compact set Z . Thus we can use a Frank-Wolfe [START_REF] Frank | An algorithm for quadratic programming[END_REF] algorithm. Note that it is similar to the one proposed by Lacoste-Julien et al. [START_REF] Lacoste-Julien | Block-coordinate Frank-Wolfe optimization for structural SVMs[END_REF] but with supplementary term due to the concavity of the loss.

Experiments

We applied our method to the task of learning a good similarity measure for aligning audio signals. In this field researchers have spent a lot of efforts in designing wellsuited and meaningful features [START_REF] Joder | Learning optimal features for polyphonic audio-to-score alignment[END_REF][START_REF] Cont | Evaluation of real-time audioto-score alignment[END_REF]. But the problem of combining these features for aligning temporal sequences is still challenging. [START_REF] Kirchhoff | Evaluation of features for audio-to-audio alignment[END_REF] Dataset description. First, we applied our method on the dataset of Kirchhoff and Lerch [START_REF] Kirchhoff | Evaluation of features for audio-to-audio alignment[END_REF]. In this dataset, pairs of aligned examples (A i , B i ) are artificially created by stretching an original audio signal. That way the groundtruth alignment Y i is known and thus the data falls into our setting A more precise description of the dataset can be found in [START_REF] Kirchhoff | Evaluation of features for audio-to-audio alignment[END_REF].

Dataset of Kirchhoff and Lerch

The N = 60 pairs are stretched along two different tempo curves. Each signal is made of 30s of music that are divided in frames of 46ms with a hopsize of 23ms, thus leading to a typical length of the signals of T ≈ 1300 in our setting. We keep p = 11 features simple to implement and that are known to perform well for alignment tasks [START_REF] Kirchhoff | Evaluation of features for audio-to-audio alignment[END_REF]. Those were: five MFCC [START_REF] Gold | Speech and audio signal processing: processing and perception of speech and music[END_REF] (labeled M1,. . . ,M5 in Fig. 3), the spectral flatness (SF), the spectral centroid (SC), the spectral spread (SS), the maximum of the envelope (Max), and the power level of each frame (Pow), see [START_REF] Kirchhoff | Evaluation of features for audio-to-audio alignment[END_REF] for more details on the computation of the features. We normalize each feature by subtracting the median value and dividing by the standard deviation to the median, as audio data are subject to outliers.

Experiments. We conducted the following experiment: for each individual feature, we perform alignment using dynamic time warping algorithm and evaluate the performance of this single feature in terms of losses typically used to asses performance in this setting [START_REF] Kirchhoff | Evaluation of features for audio-to-audio alignment[END_REF]. In Fig. 3, we report the results of these experiments.

Then, we plug these data into our method, using the Hamming loss to learn a linear positive combination of these features. The result is reported in Fig 3 . Thus, combining these features on this dataset yields to better performances than only considering a single feature.

Chorales dataset

Dataset. The Bach 10 dataset2 consists in ten J. S. Bach's Chorales (small quadriphonic pieces). For each Chorale, a MIDI reference file corresponding to the "score", or basically a representation of the partition. The alignments between the MIDI files and the audio file are given, thus we have converted these MIDI files into audio following what is classically done for alignment (see e.g, [START_REF] Hu | Polyphonic audio matching and alignment for music retrieval[END_REF]). That way we fall into the audio-to-audio framework in which our technique apply. Each piece of music is approximately 25s long, leading to similar signal length.

Experiments. We use the same features as in Sec. 5.1. As depicted in Fig. 4, the optimization with Hamming loss performs poorly on this dataset. In fact, the best individual feature performance is far better than the performance of the learned W . Thus metric learning with the "practical" Hamming loss performs much worse than the best single feature. Then, we conducted the same learning experiment with the symetrized area loss ℓ S . The resulting learned parameter is far better than the one learned using the Hamming loss. We get a performance that is similar to the one of the best feature. Note that these features were handcrafted and reaching their performance on this hard task with only a few training instances is already challenging. In Fig. 2, we have depicted the result, for a learned parameter W , of the loss augmented decoding performed either using the area. As it is known for structured SVM, this represents the most violated constraint [START_REF] Tsochantaridis | Large margin methods for structured and interdependent output variables[END_REF]. We can see that the most violated constraint for the Hamming loss leads to an alignment which is totally unrelated to the groundtruth alignment whereas the one for the symmetrized area loss is far closer and much more discriminative.

Feature selection

Last, we conducted experiments over the same datasets. Starting from low level features, namely the 13 leading MFCCs coefficients and their first two derivatives, we learn a linear combination of these that achieves good alignment performance in terms of the area loss.Note that very little musical prior knowledge is put into these. Moreover we either improve on the best handcrafted feature on the dataset of [START_REF] Kirchhoff | Evaluation of features for audio-to-audio alignment[END_REF] or perform similarly. On both datasets, the performance of learned combination of handcrafted features performed similarly to the combination of these 39 MFCCs coefficients

(1) (2) (3) (4) (5

Conclusion

In this paper, we have presented a structured prediction framework for learning the metric in temporal alignment problems. We were able to combine hand-crafted features, as well as building automatically new state-of-the-art features from basic low-level information with very little expert knowledge. Technically, this is made possible by considering a loss beyond the usual Hamming loss which is typically used because it is "practical" within a structured prediction framework (linear in the output representation). The present work may be extended in several ways, the main one being to consider cases where only partial information about the alignments is available. This is often the case in music [START_REF] Cont | Evaluation of real-time audioto-score alignment[END_REF] or bioinformatics applications. Note a simple alternating optimization between metric learning and constrained alignment provide a simple first solution, which could probably be improved upon. solutions. We call the maximization of the hinge over Y the loss augmented infer-Now we can write a new optimization objective:

min W ∈W λ 2 W 2 2 + max (Z1,...,Zn)∈Y 1 n n i=1 ℓ i (Z i , Y i ) -W, ψ i (Z i ) . (12)
Note that since our joint feature map φ(X

ℓ i (Y, Y i ) -W, ψ i (Z i ) }.
(13) From the above equation, we deduce the following relation linking primal variable W and dual variables (Z 1 , . . . , Z n ) ∈ Y(X 1 ) × . . . × Y(X n ): In the specific case when W is unconstrained and simply equals to R p×p , this reduces to:

W = 1 λ n i=1 ψ i (Z i ). ( 14
)
If W is the set of symmetric semidefinite positive matrices we get:

W = 1 λ n i=1 (ψ i (Z i )) + , (15)
with (ψ i (Z i )) + the projection of (ψ i (Z i )) over W.

Eventually if we consider W the set of diagonal matrices, and denote by Diag the operator associating to a matrix the matrix composed of its diagonal:

W = 1 λ n i=1 Diag(ψ i (Z i )). (16)
These relations are also known as the "representer theorems".

For what follows we consider the case of W = R p×p but dealing with the other cases is similar.

In that case the dual can be written simply as: max (Z1,...,Zn)∈Y(X 1 )×...×Y(X n )

- 1 2λn 2 n i=1 ψ i (Z i ) 2 F + 1 n n i=1 ℓ(Y i , Z i ). (17)
We recover a result similar to the ones of Lacoste-Julien et al. [START_REF] Lacoste-Julien | Block-coordinate Frank-Wolfe optimization for structural SVMs[END_REF].

Cumulated affinity matrix: T, S ← size(C), D ← zeros(T + 1, T B + 1) for i 1 to T A do D(i, 0) ← -∞ end for for j = 1 to T B do D(0, j) ← -∞ end for for i = 1 to T A do for j = 1 to T B do D(i, j) ← C(i, j) + max(D(i -1, j), D(i, j -1), D(i -1, j -1)) end for end for Backtracking: Y ← zeros(T A , T B ), i ← T A , j ← T B while i > 1 or j > 1 do Y (i, j) ← 1 if i == 1 then j ← j -1 else if j == 1 then i ← i -1 else m ← max(D(i -1, j), D(i, j -1), D(i -1, j -1)) if D(i -1, j) == m then i ← i -1 else if D(i, j -1) == m then j ← j -1 else i ← i -1, j ← j -1 end if end if end while return Y Figure 5: The dynamic time-warping algorithm that solves the LP (2), for a given similarity matrix C.

Figure 1 :

 1 Figure 1: Example of two valid alignments encoded by matrices Y 1 and Y 2 . Red upper triangles show the (i, j) such that Y 1i,j = 1, and the blue lower ones show the (i, j) such that Y 2 i,j = 1. The grey zone corresponds to the area loss δ abs between Y 1 and Y 2 , whereas the δ max loss corresponds to the maximum of the δ t , t ∈ {1, . . . , T 1 }.

Figure 2 :

 2 Figure 2: On the real world Bach chorales dataset, we have represented the most violated constrained at the end of learning, when the training loss is the Hamming one or the symmetrized area loss. Note also that, in terms of Hamming loss the most violated contraint for ℓ S and the Hamming one are the same.

W

  Pow M1 SC M4 SR SF M3 Max SS M2

Figure 3 :

 3 Figure 3: Comparison of performance between individual features and the learned metric. The left bars are mean δ abs error in frame on the dataset and the right are δ max . Error bars for the performance of the learned metric were determined with the best and the worst performance on 5 different experiments.

  Performance of our algorithms on the Chorales dataset. From left to right: (1) Best single feature, (2) Best learned combination of features using the symmetrized area loss ℓ S , (3) Best combination of MFCC and derivatives learned with ℓ S , (4) Best combination of MFCCs and derivatives learned with Hamming loss, (5) Best combination of features of[START_REF] Kirchhoff | Evaluation of features for audio-to-audio alignment[END_REF] using Hamming loss.

  i , Y ) is linear in Y , if ℓ is linear as well (for instance if ℓ is the Hamming loss), this problem is strictly equivalent to[START_REF] Hu | Polyphonic audio matching and alignment for music retrieval[END_REF] since in that case, the loss-augmented inference is a LP over Y(X i ), which has necessary a solution in Y(X i ) (see, e.g, [Prop. B.21] of In general, in order to be convex and thus tractable, the aforrementioned problem requires a loss which is concave over the convex sets Y(X i ).A.2 Dual of the structured SVM Since Prob. (11) is in saddle point form, we get the dual by switching the max and the min: max (Z1,...,Zn)∈Y(X 1 )×...×Y(X n )

	min W ∈W	λ 2	W 2 2 + {	1 n	n i=1

We will see that it is necessary to have fully labelled instances, which means that for each pair X i we need an exact alignment Y i between A i and B i . Partial alignment might be dealt with by alternating between metric learning and constrained alignment.

http://music.cs.northwestern.edu/data/Bach10.html.

Acknowledgements

We acknowledge the support of the GARGANTUA project (Mastodons program of CNRS), the grant SIERRA-23999 from the European Research Council and a PhD fellowship from the EADS Foundation.

A Derivation of the BCFW-like algorithm for the quadratic loss.

A.1 Relaxing the set for loss augmented inference.

Let us start from the global structured objective equation of the paper. Recall that we dispose of the training examples ((X 1 , Y 1 ), . . . , (X n , Y n )). In order to make the derivation easier, and following Lacoste-Julien et al. [START_REF] Lacoste-Julien | Block-coordinate Frank-Wolfe optimization for structural SVMs[END_REF], we denote the difference between the feature map associated to any Y ∈ Y(X i ) and the one associated to the true training example label Y i by: Tr

. The objective of structured prediction is thus:

The term max Y ∈Y(X i ) ℓ i (Y, Y i ) -W, ψ i (Y ) corresponds to the structural hinge loss for our problem. Let us introduce Y(X i ) the convex hull of the sets Y(X i ). We will also use Y = Y(X 1 ) × . . . × Y(X n ). From now on, we will perform the loss augmented decoding on this relaxed set. This problem has potentially non integral A.3 A Frank-Wolfe algorithm for solving Prob. [START_REF] Lajugie | Large-margin metric learning for constrained partitioning problems[END_REF] Now, we can derivate a Frank-Wolfe solving the dual problem of 17. As noted in the paper, we are able to maximize or minimize any linear form over the sets Y(X i ), thus we are able to solve LPs over the convex hulls Y(X i ) of such sets.

Plugging back the specific form of our joint feature map directly into Eq. ( 17) we get that

T and thus we can write the dual problem as:

Now, as in the paper, let us introduce

If U i is the matrix of ones of the same size as Z i , we consider the following loss:

This loss is sound for alignments problems since, when Y i and Z i are in Y, this is simply the ℓ S loss LY i -

F . Thus we get the following overall dual objective:

We recall that D is a diagonal matrix such that A ⊤ A -D 0 and thus our objective is convex. Our dynamic programming algorithm (DTW) is able to maximize any linear function over the sets. Thus we can use a Frank-Wolfe [START_REF] Frank | An algorithm for quadratic programming[END_REF] technique. At iteration t, this algorithm iteratively computes a linearization of the function at the current point (Z 1 , . . . Z n ) k , computes a linearization of the function, optimize it, get a new point (Z 1 , . . . Z n ) * k and then make a convex combination using a stepsize γ. Note that we have directly a stochastic version of such an algorithm. As noted in Lacoste-Julien et al. [START_REF] Lacoste-Julien | Block-coordinate Frank-Wolfe optimization for structural SVMs[END_REF] instead of computing a gradient for each block of variable Z i , we simply need to choose randomly one block at each timestep and make an update on these variables. The linearization simply consists in computing the matrix gradient for each of the matrix variables Z i which turns out to be:

where C is simply the affinity matrix of dynamic time warping.

B The dynamic time warping algorithm

Let us the pseudocode of the dynamic time that maximize the LP (2) of the article.In opposition to Müller [START_REF] Müller | Information retrieval for music and motion[END_REF], we give a version of the algorithm for the affinity matrix C. Intuitively, the cost matrix is the opposite of a cost matrix, thus we aim to maximize the cumulated affinity instead of minimizing the cumulated cost. This algorithm is O(T A T B ), making it very costly to compute for large time series.


Mauricio Dwek
email: mauricio.dwek@grenoble-inp.fr

Peggy Zwolinski

MODELING RECYCLABILITY: HOW CAN WE PREDICT THE EVOLUTION OF RECYCLING CHAINS?

Keywords: recycling chains, recyclability, MFA, design

Sustaining the material basis of society is becoming increasingly complex as demand grows for ever more sophisticated materials. As concerns regarding material availability arise, many studies following material stocks and flows and defining material criticality are being conducted. These assessments provide information that can be decisive for the industrial deployment of sustainable and innovative technologies. Disruption risks to the supply chains must be predicted to prevent material shortages at the corporate, national and global scale. Designers can play a major role for the preservation of material resources by considering the evolution of material availability at the material selection stage. Bearing in mind the lifecycle of any product, material recyclability will progressively become a key factor for the design process. This paper is part of an on-going research being conducted at the G-SCOP lab of the Grenoble Institute of Technology, whose aim is to provide dynamic resource scenarios and additional input to Design for Recycling guidelines, so as to assist material selection in the design process. It presents a framework to identify the parameters that determine the evolution of recycling chains, based on material flow data and historical accounts of the shifts, ascent and decline of recycling activities. This should allow designers to incorporate material criticality and recyclability to their lifecycle assessments.

INTRODUCTION

In the last 25 years, the array of materials synthesized in laboratories and present in industrial applications has expanded considerably. Not only has the economy become omnivorous in its consumption of elements from the periodic table [START_REF] Greenfield | The omnivorous diet of modern technology[END_REF], but scientists have also gone to great lengths to produce materials that are ever more sophisticated, complex and architected, down to the nano-scale. The advances in sectors such as carbonlean energies, information and communication technologies have led to a steep increase in demand for materials whose reserves seem inadequate to fulfil future scenarios [START_REF] Graedel | Will metal scarcity impede routine industrial use ?[END_REF]- [START_REF] Erdmann | Criticality of Non-Fuel Minerals : A Review of Major Approaches and Analyses[END_REF]. Yet, despite the economic and ecological need of preserving resources through recycling, by establishing what is called a closed-loop economy, most industries and material valuechains are still open-looped. This means that materials follow more complex routes and a wider range of end-of-life applications are currently in place. Today, material lifecycles look more like a spiral or a coil than a circle.

Even though research has been performed on recyclable materials and more efficient Theme: Sustainable and eco-design of processes, products and services 2/10 recycling processes, secondary lifecycles in the industry are still commonly placed beyond the range of the design process. In order to promote product and material recycling, designers have compiled Design for Recycling (DfR) guidelines for more than two decades, yet these seem to "lack a combination of concrete instructions, prioritization, and recyclability performance feedback" [START_REF] Peters | Prioritizing 'Design for Recyclability' Guidelines, Bridging the Gap between Recyclers and Product Developers[END_REF]. In many cases, designers have little to no contact with the recyclers of their products. How is the loop to be closed if both ends never meet?

Recyclers have gathered decades of recycling experiences and there is enough data today to evaluate how -and above all why -recycling chains tend to become necessary and grow, stagnate or disappear depending on the social, economic, political and technological environment of their time. This information is capital to evaluate how designers' decisions affect the shape of material flows and meet the material needs that future technology will require. This paper is part of an on-going research being conducted at the G-SCOP lab of the Grenoble Institute of Technology, whose aim is to provide dynamic resource scenarios and additional input to Design for Recycling guidelines, so as to assist material selection in the design process. It presents a framework to identify the parameters that determine the evolution of recycling chains, based on material flow data and historical accounts of the shifts, ascent and decline of recycling activities. This should allow designers to incorporate material criticality and recyclability to lifecycle assessments.

THE LINK BETWEEN AVAILABILITY, VALUE AND RECYCLABILITY

Addressing material availability in design

Designers are facing interesting times. They are expected to constantly provide technological innovations that enhance everyday life, improve production systems, create value and reduce environmental impacts. However, in their attempt to meet these expectations on a global level, a new constraint has surfaced: material scarcity. The increasing demand for manufactured goods and energy resources in developing markets in Asia, Africa and Latin America has led to an intensification of the consumption of commodities in general and some materials in particular. Many of today's information and communication technologies, as well as most carbon-lean energy systems, depend on materials whose future availability is uncertain at best. In their assessment of the risks to the supply chain associated with material availability, Alonso et al. [START_REF] Alonso | Material availability and the supply chain: risks, effects, and responses[END_REF] state that "over the long term, market forces and technology will effectively ensure that responses such as substitution and recycling will occur" in case of material scarcity or disruption. They also suggest that material selection decision-making is a means to develop a strategy to reduce vulnerability and mitigate material shortages.

Increasingly, mastering material sourcing will become a key issue for industries developing sophisticated materials for innovative products. The designer should consider the present and also future availability -and therefore value -of material supplies during the design process. Thus, it is crucial to address material availability and recyclability at the material selection stage, aiming for more efficient material uses and recovery at the end of the product's lifecycle.

Material value considerations

The value of a given material depends on several factors regarding supply and demand issues. Supply is defined by the raw materials' extraction and production processes. In the case of ceramics and metals, geological occurrence and concentration are key factors that involve geopolitical relations and may be a cause for concern in some cases, especially when few countries possess the bulk of reserves. Mining activities also depend on long-term capital-intensive investments that are usually based on feasibility studies attempting to anticipate profitability and market fluctuations. Some ceramic materials and most metals are the by-products of the extraction of major carrier metals that constitute the ore and have less efficient processing rates. Price increases of these companion materials may encourage improved recovery of these materials. In the case of polymers, the fluctuations of fossil fuel reserves affect supply to the point that investments in plastic recycling and the so-called bio-plastics may be fostered or hindered depending on the rise and fall of oil prices.

In terms of demand, Binder et al. [START_REF] Binder | Explanatory Variables for per Capita Stocks and Flows of Copper and Zinc A Comparative Statistical Analysis[END_REF] have indicated that metal use increases with income growth, a statement that could be generalized to other commodities such as ceramics and polymers. Thus, since global population and economy are growing, one could consider the consumption of some materials as ever expanding in the near future [START_REF] Graedel | Metal spectra as indicators of development[END_REF]. This may not be true in some cases, where predictions have shown it may saturate [START_REF] Müller | Exploring the engine of anthropogenic iron cycles[END_REF]. Furthermore, adding a note of unpredictability to the mix, the very own innovation process can disrupt any forecasting model of demand, as technological innovations relying on new materials can quickly transform the demand of little-known, undervalued elements of the periodic table into overnight industrial favorites [START_REF] Graedel | Will metal scarcity impede routine industrial use ?[END_REF].

Recycling issues

Recycling has therefore become a strategic issue due to the realization of the intrinsic value of materials, to the point that urban mining -as recycling is now referred to -is considered an alternative to conventional mining. For instance, data from recyclers has shown that the concentrations of precious metals such as gold, platinum, palladium and silver, along with other metals such as copper, aluminium, iron, tin, nickel, gallium and germanium, present in current electric and electronic equipment waste, are superior than the concentrations of these metals in the virgin ores [START_REF] Johnson | Dining at the Periodic Table : Metals Concentrations as They Relate to Recycling[END_REF]. However, not every available metal is being recovered from Waste Electrical and Electronic Equipment (WEEE) as flows are still being mustered.

Recycling is hampered by a number of factors related to the management of waste flows as well as to intrinsic properties of the materials themselves. Usually, profitability, and therefore viability, of material recycling is considered to be relatively poor, which, added to the volatility of recycled material prices and the competition with original raw materials, makes investing in recycling operations a high-risk venture [START_REF] Mathieux | Product ecodesign and materials: current status and future prospects[END_REF]. Recycling is not always viable depending on the material's applications. While some uses are both technically and economically feasible, others can involve recycling technology that may not be profitable yet due to lack of reserve waste volumes or waste prices. In some cases, recycling is inherently not feasible due to material dissipation during the lifecycle. Moreover, recycling operators often struggle due to lack of expertise for managing complex inputs from end-of-life streams, since the technological, economic and thermodynamic knowledge required to process waste is not the same as the one required for processing raw materials and ore [START_REF] Unep | Metal Recycling: Opportunities, Limits, Infrastructure[END_REF].

In literature, there are several definitions of recyclability. In broad terms, recyclability could be defined as the ability of a material to be managed at the end of its lifecycle, collected, treated and re-enter the value chain, in an open or closed loop. For Villalba et al., "recyclability is the ability a material has to reacquire the same properties it originally had" [START_REF] Villalba | A proposal for quantifying the recyclability of materials[END_REF]. Peters et al. assume a designer's point of view and define recyclability as "the affordance a product has for recovering as much components and materials as possible (quantity) with the highest possible purity (quality) by the least amount of effort (ease) with existing recycling technologies" [START_REF] Peters | Prioritizing 'Design for Recyclability' Guidelines, Bridging the Gap between Recyclers and Product Developers[END_REF]. The French ADEME/AFNOR BPX30-323 environmental labeling guidelines define recyclability by either the recycled content of a given material in the products it is made with, or the recycling rate at the end of life, depending on the maturity of the recycling market: in mature cases, for which most scrap is recovered and recycled, the recycling rate is considered, whereas in developing recycling markets, the reward goes to products that contain secondary material, so as to encourage these burgeoning activities [START_REF] Thomas | Methodologies to measure the sustainability of materials -focus on recycling aspects[END_REF]. Both cases assume a more practical point of view, based on actual data from the return of materials in industrial cycles. Yet the recycled content approach is still disputed for being "highly dependent of the system conditions, export share, growth and the existing open loop recycling practice" [START_REF] Rombach | Time consideration in mass flow analysis -A contradiction?[END_REF]. Mizuno et al. state that although DfR guidelines propose a theoretical calculation of the recyclability rate as "an expectation value of a rate of recyclable materials to total weight of the product at the design stage", the recyclability rate of most products depends on a number of factors such as "legislation, economic situation, corporate policy, recycling technology, recycling processes, etc" [START_REF] Mizuno | Recyclability Evaluation of LCD TVs Based on End-of-Life Scenarios[END_REF].

Several recyclability rates are used in different contexts so as to better convey the potential recycling of materials but most of them circumvent the contextual information and are more interested in the present-time, practical and static value of the recycling rates rather than in understanding their dynamics. While they are useful for policy-making purposes, these indexes provide no information as to how the secondary material is employed, whether there is potential for improving these rates and if a closed-loop is attained or attainable.

FRAMEWORK FOR STUDYING THE EVOLUTION OF RECYCLING CHAINS Scope delimitation

Materials differ in their flows, history and context so that each material has its own specificities. Nevertheless, some materials share the same qualities and can be compared to a certain degree, especially inside material classes. This comparison is interesting to obtain the factors that characterize anthropogenic cycles as a whole and those that are intrinsic to given material properties. The more chains are analyzed, the more parameters are considered, and the more general the model may be. However, it would be impracticable to provide a theoretical model of recyclability based on the evolution of availability for all materials. Therefore, the study will concentrate on metals, since they possess a larger information database and have some of the most developed, mature and well-recorded recycling chains.

Critical materials

Recent concern about material shortages has led to the emergence of studies on critical materials, i.e. materials whose supply disruption risks and economic importance are particularly disquieting. Because many of these materials are necessary for the largescale deployment of strategic technologies for the near future (carbon-lean energy systems, information and communication technologies), predicting their demand and the vulnerability of the respective systems to fluctuations in supply can be decisive to ventures and their success in the long-run. Even though critical material listings vary according to the methodologies used in their assessment, they generally include rare-earth elements, the platinoid family, as well as specialty metals such as indium, scandium, tungsten and yttrium [START_REF] Erdmann | Criticality of Non-Fuel Minerals : A Review of Major Approaches and Analyses[END_REF]. Avoiding the dissipative uses of these materials and securing their availability by fostering their recycling is the focus of many studies [START_REF] Senk | Strategies for fulfilment of critical raw materials demand in Europe[END_REF]- [START_REF] Rademaker | Recycling as a Strategy against Rare Earth Element Criticality : A Systemic Evaluation of the Potential Yield of NdFeB Magnet Recycling[END_REF]. In their proposition of a methodology for metal criticality determination, Graedel et al. [START_REF] Graedel | Methodology of Metal Criticality Determination[END_REF] suggest the relevant material characteristics for evaluating criticality in a temporal perspective, according to the organizational levels, as shown in Table 1. Clarifying the context, timeframe and assessment method is fundamental to any attempt at an assessment of criticality, as well as stakeholder validation and accounting for uncertainties through sensitivity analyses [START_REF] Erdmann | Criticality of Non-Fuel Minerals : A Review of Major Approaches and Analyses[END_REF].

Material flows into the end of life

Following the flows of materials in society is the first step in establishing a model for material availability and understanding the dynamics that drive material cycles. The main tool used in these studies is Material Flow Analysis (MFA) and, more recently, Dynamic Material Flow Analysis (DMFA). MFA is considered a fundamental industrial ecology tool and has many applications due to its very synthetic display of an element's -or substance's, in the case of Substance Flow Analysis -transformations and exchanges in society or with the environment. It provides a quantitative partitioning of a material in its different life stages for a given region or time period and may serve as a solid basis for sustainability assessments [START_REF] Dittrich | The physical dimension of international trade, part 2: Indirect global resource flows between 1962 and 2005[END_REF], [START_REF] Van Der | Environmental Impact of the use of Natural Resources and Products[END_REF], urban planning [START_REF] Müller | Stock dynamics for forecasting material flows -Case study for housing in The Netherlands[END_REF]- [START_REF] Obernosterer | Urban metal management: The Example of Lead[END_REF] or policy making [START_REF] Reimann | Evaluation of environmental life cycle approaches for policy and decision making support in micro and macro level applications[END_REF]. When performed on a local level, they offer information on the social mechanisms of management systems and the economic interactions with neighboring regions, whereas analyses performed on a global scale -usually collected from several more local studies -portray general tendencies for material consumption and scarcity [START_REF] Harper | Making Metals Count : Applications of Material Flow Analysis[END_REF]. The general model for a MFA is presented in Figure 1.

Figure 1: MFA general model [START_REF] Mao | The multilevel cycle of anthropogenic lead I . Methodology[END_REF] 6/10

To the general purpose of the research presented here, MFA is essential to quickly access an important volume of data, compiled in a comprehensive manner that quantifies primary and secondary material production, stocks in the economy according to the applications therein, different waste management strategies as well as material dissipation. Several MFAs have already been conducted, constituting an extensive albeit not exhaustive database. Many studies exist for specific materials, regions and timeframes and some attempts at an exhaustive inventory have been undertaken such as the Stocks and Flows Project led by Prof. Graedel from the Center for Industrial Ecology at the Yale School of Forestry and Environmental Studies, which has been studying the anthropogenic cycles of metals.

However, the MFA database is still incomplete and, although studies exist for hazardous chemical substances [START_REF] Long | Flow analysis of heavy metals in a pilot-scale incinerator for residues from waste electrical and electronic equipment dismantling[END_REF], plastics [START_REF] Kleijn | Dynamic substance flow analysis: the delaying mechanism of stocks, with the case of PVC in Sweden[END_REF] and building materials [START_REF] Müller | Stock dynamics for forecasting material flows -Case study for housing in The Netherlands[END_REF], it is quite concentrated on metal cycles. Most analyses are also focused on developed countries, perhaps due to information availability and reliability as the flows in these areas are generally more controlled and formally registered. Most MFA studies are static and contained within a given region and timeframe, even if dynamic models are increasingly available.

One can extract from the mass balances indexes that characterize each stage of the lifecycle, such as those presented in Table 3. These coefficients allow the comparison between MFAs in regards to different timeframes, nations or materials. They also provide information on trends in the configuration of the material flows and can be used to validate claims made by the stakeholders on the evolution of the material's lifecycle, particularly in the end of life. Table 3: Lifecycle coefficients obtained from MFA ( 1 : coefficients proposed by [START_REF] Mao | The multilevel cycle of anthropogenic lead II . Results and discussion[END_REF])

Production efficiency = 1 - + Utilization efficiency 1 = 1 - + Accumulation ratio 1 =
End of life recycling ratio 1 = -+ Secondary supply ratio 1 = Discard management ratio 1 = + But, MFA is a lifecycle assessment tool that usually presents the material's lifecycle in broad strokes and is not very adapted at a micro-level [START_REF] Reimann | Evaluation of environmental life cycle approaches for policy and decision making support in micro and macro level applications[END_REF]. Moreover, static and dynamic MFAs provide little to no contextual information and authors generally infer the circumstances that have shaped the flows and hypothesize their interpretations. To complement and confirm the inferences from the MFA, a comparison with historical data is needed. A look back to the origins of the material flows, and particularly of their end-of-life management, can clarify which contextual elements have an influence on the anthropogenic cycle and, ultimately, the availability of a given material. Stakeholders such as government agencies, eco-organisms, waste collectors and recyclers are a qualitative source of information that can provide empirical information from the industry. Associating the evolution of material availability with the history of recycling chains will point out the factors that determine recyclability. This study will use Material Flow Analysis (MFA) as a primary source of data and historical records of recycling chains from literature and stakeholders as a guideline to in-8/10 gical processes, recycling is still hampered by a lack of technical optimization [START_REF] Unep | Metal Recycling: Opportunities, Limits, Infrastructure[END_REF].

The information from the MFAs and the bibliography covered so far are not conclusive and need confirmation from historical and practical data, field experts and stakeholders. This should also provide data on the secondary cycles in case of downcycling or upcycling. Stock balances are an important source of information and provide parameters such as primary source availability, primary over secondary production ratios, as well as import and export rates. DMFAs will be used to forecast the evolution of materials because they have the advantage of presenting trends by using life-span distributions to estimate derivative data from material flows.

There is a clear distinction between theoretical recyclability and practical results from WEEE chains, of up to 35% [START_REF] Mizuno | Recyclability Evaluation of LCD TVs Based on End-of-Life Scenarios[END_REF]. A 100% recyclability rate in theory does not entail in complete recovery of material. Connecting information from actual end-of-life chains with flow data can improve DfR guidelines but the criteria must account for context and product type. The preliminary parameters identified to constitute the model can be gathered into the following categories: geopolitical, economic, social/legislative, technological and logistical.

FUTURE WORK

The research will focus primarily on metals, since the MFA database is larger for this class. A probable distinction will be made among ferrous, non-ferrous, precious and specialty metals, though the latter encompasses a very broad spectrum of materials and properties and is a more heterogeneous group.

The interactions and relative impacts of each parameter will be equated and compared to the historical and present day scenarios in order to validate the model. Since the MFA database is quite varied in terms of materials, timeframes and regions covered, a general model is still distant and most hypotheses will be confirmed or disallowed by examining the experts' interviews. This will fill an important gap in the comprehension of how recycling chains evolve, since there are very few wide-ranging compilations of data describing the history of the recycling processes and industry in literature.

The integration of MFA in the design process will be detailed, especially in regards to its quick provision of information from the economy and its potential uses as a lifecycle assessment method. This approach should constitute a practical decision-making tool for material selection in design. Because material cycles are closely linked with product lifecycles, a case around WEEE will be studied to verify the applicability of the tool developed from this project.

Table 1 :

 1 Relevant material-related characteristics for different organizational levels[START_REF] Graedel | Methodology of Metal Criticality Determination[END_REF]

		Using corporation	Using nation	Global
	Focus	Relevance to that	Relevance to national	All uses of a material,
		firm's product line	industry and population	wherever they happen
	Time scale	1-5 years	5-10 years	10-100 years
	Supply potential	Crucial	Very important	Very important
	Technological	Very important	Worth consideration	Impossible to predict
	change
	Geopolitical	Crucial	Important	Unimportant
	factors
	Social factors	Moderately	Very important	Unimportant
		important
	Environmental	Important	Important	Moderately important
	implications
	Intensity of	Crucial	Depends on national	Unimportant
	competition		industry composition

terpret the evolution of the material flows. Comparing the information from MFAs and the general context is what will allow the establishment of correlations between the parameters of the recycling chains' evolution and serve as the basis for a model of recyclability. Once it is finished, the model will be a synthesis of both quantitative and qualitative data, allowing the designer to address disruption risks to the supply chain, identify the need for substitution and conceive methods to prevent material dissipation.

PRELIMINARY REMARKS

The objective of this research is to provide a tool for designers to consider material availability when developing a product. Examining whether a material will be available in the future amounts to evaluating its criticality. Indicator-based methodologies for metal criticality determination already exist (Table 4). Nevertheless, recyclability is not a consequence of criticality exclusively. A designer considering the use of a material qualified as critical could be positioned both in the context of his corporation and his nation, as shown in Table 1, rendering concerns about supply potential, technological change, geopolitical and social factors, environmental indicators and the intensity of competition, primordial to the selection of a material. These assessments should be prepared based on projections made on a 5 to 10 years span, which corresponds to both the corporate and national timeframes of Table 1, since the designer is subject to his company's commercial interests as well as his countries social demands, economic policies and legislation. Although the global perspective is important for the long-term sustainability of the planet's biosphere and technosphere, its governance depends on the adoption of an international natural resources management policy, whose enforcement and governance conditions are beyond the scope of this study.

Corporate innovation

If a material is deemed as having a relative risk in being sourced, some alternatives are possible such as enhanced dematerialization, material substitution or a closed loop of a given strategic component. However, substitution is not always possible due to the lack or poor performance of a substitute, as is the case of magnesium, manganese, yttrium, rhodium, rhenium, thallium, lead, europium and dysprosium [START_REF] Graedel | On the materials basis of modern society[END_REF]. And though primary and secondary production of metals are based on similar hydrometallurgical or pyrometallur-