Chen P Liu
email: michel.lagarde@insa-lyon.fr

Functional lipidomics of oxidized products from polyunsaturated fatty acids

Keywords: Lipid peroxidation, Oxygenases, Lipid mediators, Lipid analysis

Because of their high degree of unsaturation, polyunsaturated fatty acids (PUFA) in mammals, with mainly 18, 20 and 22 carbons, can easily be autooxidized, and converted into many oxidized derivatives and degradation products. This short review reports on some of those relevant to the evaluation of oxidative stress in situ. In addition, the enzyme-dependent oxygenation by both dioxygenases and monooxygenases is briefly reviewed by functional and/or metabolic categories, pointing out the structure variety and the analytical approaches.

Introduction

Polyunsaturated fatty acids (PUFA) of nutritional interest are highly susceptible to oxidative stress, which leads to lipid peroxidation with degradation products from primary hydroperoxides. Besides non enzymatic lipid peroxidation, enzyme-induced peroxidation, e.g. cyclooxygenases and lipoxygenases, is a process that leads to hydroperoxides specifically transformed into a series of more stable metabolites called eicosanoids, docosanoids and even octadecanoids from C20, C22 and C18 PUFA, respectively. Functional lipidomics may be seen as targeted lipidomics addressed to classes of lipids associated with specific types of function or biological system, e.g. eicosanoids in the vascular bed. This contrasts with the measurement of a limited number of markers associated with a physiological or pathophysiological situation. In this short review, the analysis of oxidized products from PUFA will be discussed in relation with their biological relevance.

Autooxidation of PUFA

It is well known that the autooxidation of PUFA leads to many metabolites that have been used as markers of the process (Catala, 2009). The mostly used one is malondialdehyde (MDA), which can also be produced through the cleavage of prostaglandin (PG) H 2 , releasing equimolecular amount of 12-hydroxy-heptadecatrienoic acid (HHT) [START_REF] Hecker | Products, kinetics, and substrate specificity of homogeneous thromboxane synthase from human platelets: development of a novel enzyme assay[END_REF], making MDA an index of cyclooxygenase pathway as well. In addition to be a global index of lipid peroxidation, it can derive from non-lipid precursors such as carbohydrates, amino acids and DNA [START_REF] Janero | Malondialdehyde and thiobarbituric acid-reactivity as diagnostic indices of lipid peroxidation and peroxidative tissue injury[END_REF]. Despite this lack of specificity, the advantage of its measurement in urine for instance is to provide an integrative and non invasive assessment of overall oxidative stress in the body [START_REF] Guichardant | Malondialdehyde measurement in urine[END_REF]. In contrast to this lack of specificity, products such as isoprostanes and neuroprostanes are quite specific, perhaps too much as they represent the peroxidation of only two precursors, i.e. arachidonic acid (ARA) and docosahexaenoic acid (DHA), respectively (Roberts & Fessel 2004). In addition, they are quite numerous with 64 different isoprostane isomers [START_REF] Morrow | The isoprostanes: unique bioactive products of lipid peroxidation[END_REF]. Among them, at least two isoprostanes exert biological activities. 8-iso-PGF 2 and 8-Iso-PGE 2 exhibit vasoconstricting effects [START_REF] Fukunaga | Evidence for the existence of F2-isoprostane receptors on rat vascular smooth muscle cells[END_REF], with the latter also affecting platelet aggregation [START_REF] Lahaie | A novel mechanism for vasoconstrictor action of 8-isoprostaglandin F2 alpha on retinal vessels[END_REF][START_REF] Longmire | Actions of the E2-isoprostane, 8-ISO-PGE2, on the platelet thromboxane/endoperoxide receptor in humans and rats: additional evidence for the existence of a unique isoprostane receptor[END_REF]. On the other hand, neuroprostanes could be sensitive biomarkers of brain injury in response to oxidative stress, and then be relevant to neurodegenerative diseases (Montime et al. 2004).

We have been interested in hydroxy-alkenals that may be issued from the whole series of omega-3 or omega-6 PUFA, in contrast to isoprostanes and neuroprostanes that are products of ARA and DHA only. In that case, 4-hydroxy-hexenal (4-HHE) and 4-hydroxy-nonenal (4-HNE) are indices of omega-3 and omega-6 peroxidation (or distal peroxidation in the esterified fatty chain), respectively [START_REF] Guichardant | Hydroxy-alkenals from the peroxidation of n-3 and n-6 fatty acids and urinary metabolites[END_REF][START_REF] Bacot | Evidence for in situ ethanolamine phospholipid adducts with hydroxyalkenals[END_REF]) (Figure 1). It can however be argued that hydroxyalkenals are reactive enough to make covalent adducts with protein residues (on thiol and amine groups) as well as with amino-phospholipids and other bio-amines [START_REF] Jürgens | Detection of new epitopes formed upon oxidation of lowdensity lipoprotein, lipoprotein (a) and very-low-density lipoprotein. Use of an antiserum against 4hydroxynonenal-modified low-density lipoprotein[END_REF][START_REF] Guichardant | Covalent modifications of aminophospholipids by 4-hydroxynonenal[END_REF]. Therefore, the measurement of 4-HHE and 4-HNE only represents the remaining free molecules. Yet, this argument is also valid for chemically stable isoprostanes, and neuroprostanes, that are very likely (although not reported) to be beta-oxidized as described for prostaglandins, and then found partly in urine as dinor and tetranor derivatives [START_REF] Diczfalusy | Beta-oxidation of eicosanoids[END_REF]. The advantage of measuring hydroxy-alkenals is the possibility to include in the same run of analysis the measurement of other homologues such as 4-hydroxy-dodecadienal (Figure 2) derived from the 12-lipoxygenase product of ARA, 12-HpETE [START_REF] Bacot | Evidence for in situ ethanolamine phospholipid adducts with hydroxyalkenals[END_REF]). These molecules have been measured in blood plasma [START_REF] Calzada | Subgram daily supplementation with docosahexaenoic acid protects low-density lipoproteins from oxidation in healthy men[END_REF], and a non invasive measurement of their acidic metabolites (Figure 1), 4hydroxy-hexaenoic/nonenoic acids in urine is a possible alternative [START_REF] Guichardant | Hydroxy-alkenals from the peroxidation of n-3 and n-6 fatty acids and urinary metabolites[END_REF]. In terms of biological effects, recent results show their cytotoxic activity [START_REF] Pillon | Quantitative structure-activity relationship for 4-hydroxy-2-alkenal induced cytotoxicity in L6 muscle cells[END_REF]. Hydroperoxy derivatives of PUFA can also be reduced into their stable counterparts by glutathionedependent peroxidase [START_REF] Foster | Selenium in health and disease: a review[END_REF], namely hydroxy derivatives, and then escape the spontaneous cleavage into hydroxy-alkenals. In that case, the various hydroxy-eicosatetraenoic acids (HETEs) from ARA and hydroxy-octadecadienoic acids (HODEs) from linoleic acid can easily be measured (Figure 3). As an example of clinical investigation, 9-and 13-HODE have been found as the main hydroxy derivatives present in LDL, and are significantly increased in LDL from type 2 diabetic patients compared to control LDL (Colas, 2010). However, it should be noted that HETEs and HODEs may also derive from the action of lipoxygenases/glutathione peroxidase, making these markers as indices of both autooxidation and enzyme-dependent oxygenation. The only difference between autooxidation product and lipoxygenase-derived products is the R/S racemic configuration in the former and the S configuration in the latter. The biological function of HETEs and HODEs has been thoroughly reviewed [START_REF] Spector | Hydroxyeicosatetraenoic acids (HETEs)[END_REF]).

Dioxygenase-dependent oxygenated products

The main dioxygenases involved in oxygenation of PUFA are cyclooxygenases and lipoxygenases. The former enzymes lead to the unstable cyclic prostaglandin endoperoxides, PGH, from dihomogamma-linolenic, arachidonic and eicosapentaenoic acids, namely PGH 1 , H 2 and H 3 , because the enzyme activity is associated with hydroperoxidase in the same protein (PGH synthase) [START_REF] Smith | Nutritionally essential fatty acids and biologically indispensable cyclooxygenases[END_REF]. A series of PGD, E and F  synthases convert PGHs into primary and ubiquitous PGs. More tissuespecific isomerases such as thromboxane (Tx) and prostacyclin (PGI) synthases convert PGHs into TxAs and PGIs, further hydrolyzed into TxBs and 6-oxo-PGFs (Figure 4) due to their instability. Their biological activities are diverse and tissue-specific [START_REF] Smith | Cyclooxygenases: structural, cellular, and molecular biology[END_REF]. It should be however noticed that PGI 1 cannot be formed [START_REF] Bunting | An antiserum to 5,6-dihydro prostacyclin (PGI1) which also binds prostacyclin[END_REF]) and TxA 3 is hardly produced [START_REF] Boukhchache | Interactions between prostaglandin precursors during their oxygenation by human platelets[END_REF]. Dihomo counterparts of PG 2 s can be produced from docosatetraenoic or adrenic acid [START_REF] Sprecher | Dihomo-prostaglandins andthromboxane. A prostaglandin family from adrenic acid that may be preferentially synthesized in the kidney[END_REF]. It might be worth to measure all those products by GC-MS or LC-MS (Samuelsson et al. 1975;[START_REF] Deems | Detection and quantitation of eicosanoids via high performance liquid chromatography-electrospray ionizationmass spectrometry[END_REF] to get an overview of the cyclooxygenation of PG precursors. Lipoxygenases are less substrate-specific than cyclooxygenases as they simply require a 1,4-cis,cispentadiene structure (Ewans & Sprecher 1985). The most classical ones are depicted as 5-, 12-and 15-lipoxygenases, based on the number of carbon to be oxygenated in arachidonic acid as a substrate [START_REF] Yamamoto | Mammalian lipoxygenases: molecular and catalytic properties[END_REF]. To be extended at many other possible PUFA substrates, 12-and 15lipoxygenases could better be named omega-9 and omega-6 lipoxygenases, respectively. The primary products are hydroperoxy derivatives that are easily reduced into hydroxyl counterparts if the biological system of interest is well endowed with glutathione peroxidase activity (Figure 5). The 5lipoxygenase products of arachidonic acid, mainly leukotrienes, have been involved in many pathophysiological states, especially inflammation [START_REF] Samuelsson | Leukotrienes and other lipoxygenase products[END_REF]. The scheme may however be much more complicated if we consider that primary products of PUFA may be further oxygenated if they contain the required additional double bonds [START_REF] Serhan | Lipoxins: a new series of eicosanoids (biosynthesis, stereochemistry, and biological activities)[END_REF]. Because of substantial absorption molecular coefficient of conjugated dienes, trienes, even tetraenes which are always present in those products, LC separation with UV detection may sometimes be sensitive enough. When coupled with mass spectrometry, the LC-MS combination will increase the sensitivity and specificity [START_REF] Poulsen | Identification of inflammatory and proresolving lipid mediators in bone marrow and their lipidomic profiles with ovariectomy and omega-3 intake[END_REF]). More recently, numerous bioactive hydroxyl derivatives of DHA and eicosapentaenoic acid have been described by [START_REF] Serhan | Novel lipid mediators and resolution mechanisms in acute inflammation: to resolve or not?[END_REF], and called maresins, protectins and resolvins for their various activities in the anti-inflammatory processes. In addition, a double lipoxygenation of PUFA has been pointed out with a family of compounds that share platelet anti-aggregatory properties because of a conjugated trans,cis,trans (E,Z,E) double bond geometry resulting from the double oxygenation [START_REF] Chen | Full characterization of PDX, a neuroprotectin/protectin D1 isomer, which inhibits blood platelet aggregation[END_REF][START_REF] Chen | Poxytrins, a class of oxygenated products from polyunsaturated fatty acids, potently inhibit blood platelet aggregation[END_REF]. These compounds, issued from various PUFA owning at least three double bonds at the omega-6, 9 and -12 positions, have been named poxytrins. They are dihydroxylated derivatives obtained by double lipoxygenation of PUFA, followed by reduction of the intermediate hydroperoxides. Although their biological relevance remains to be ascertained as endogenous products, their LC separation for a putative measurement in biological samples is to be considered in the same run together with the intermediates from single lipoxygenation, also reduced into mono-hydroxy derivatives. An example of separation of poxytrins and mono-hydroxy intermediates is given with a discrimination based on their max (Figure 6). Although most of the direct enzyme-dependent oxygenation leads to the S stereo-configuration, some R alcoholic carbons can be found, and the chromatographic step in the separation and evaluation is crucial, sometimes requiring chiral stationary phases.

Monooxygenation of PUFA or derivatives

Another important oxygenation pathway of PUFA is distinct from dioxygenase-dependent and enzyme-independent peroxidation by using monoxygenases that require cytochrome P 450 . Those monooxygenases may hydroxylate the substrate, making a primary (omega oxidation) or secondary alcohol, or making an epoxide at the expense of one double bond (Figure 7). Their biological relevance is less accurately stated, and this pathway for PUFA has not been so well developed compared to dioxygenase pathways. The omega hydroxylation has been sometimes described as a degradation pathway, as exemplified with leukotriene B 4 that is first transformed into 20-hydroxy-LTB 4 before being further oxidized into its 20 carboxy derivative which is entirely devoid of chemotactic activity (Shak & Goldstein 1985). Omega or omega-1 hydroxylation has also been described with PUFA, especially arachidonic acid, to provide vasoactive derivatives (Capdevilla et al. 2000). Another development has been made with the epoxidation of PUFA through cytochrome P 450 enzymes that have been called epoxygenases. The most popular products again are those issued from arachidonic acid (Capdevilla et al. 1992). The epoxy-eicosatrienoic acids (EET) formed are endowed of vasoactivities as well as their vicinal dihydroxy derivatives produced by specific epoxide hydrolases (Figure 7). The products are bioactive as well [START_REF] Spector | Action of epoxyeicosatrienoic acids on cellular function[END_REF]. Most of those products being position isomers without any conjugated double bonds, their separation can be achieved by HPLC but their measurement cannot be done in proximal UV. One way of measuring the non-derivatized compounds is to use tandem mass spectrometry coupled with LC [START_REF] Newman | The simultaneous quantification of cytochrome P450 dependent linoleate and arachidonate metabolites in urine by HPLC-MS/MS[END_REF].

Conclusion

Eicosanoids as well as docosanoids and octadecanoids are so numerous, with still new compounds of biological interest to be described, that they represent a typical example of what the functional lipidomics is made for, i.e. describe and quantify as much as possible the different lipid species to get proper structure-function relationships, and evaluate their relevance in a given physiological or pathophysiological situation. As an exhaustive analysis is rarely possible, even when targeting a limited functional field, the measurement of an array of lipid metabolites in a specific field of interest allows to get the best balance sheet of the different pathways involved and bring information on their function. with the production of PGG 1 /H 1 and PGG 3 /H 3 , respectively. The formation of the other prostanoids also occurs from the latter two PUFA, except for PGI 1 that cannot be formed because of the lack of double bond at carbon 5. The same metabolic pathways have been reported for dihomo-arachidonic or adrenic acid with the products called dihomo-prostanoids. PG: prostaglandin; Tx: thromboxane; PG-S: PG synthase; PGH-S is a bi-functional enzyme protein with two activities: cyclooxygenase producing PGG, and hydroperoxidase converting PGG into PGH.

Figure legends

Figure 1 :

 1 Figure 1: Distal peroxidation and subsequent degradation by cleavage of the cis (Z) double bond of the conjugated diene leading to a hydroxyl-alkenal. 4-hydroxy-2E-nonenal (4HNE) and 4-hydroxy-2Ehexenal (4-HHE) derives from omega-6 and omega-3 PUFA, respectively. Further oxidation of the hydroxyl-alkenal leads to corresponding carboxylic acids, 4-hydroxy-2E-nonenoic (4-HNA) and 4hydroxy-2E-hexenoic (4-HHA) acids.

Figure 2 :

 2 Figure 2: 12-lipoxygenation of an omega-6 polyunsaturated fatty acid (PUFA), such as arachidonic acid, and subsequent degradation/cleavage leading to 4-hydroxy-2E-dodecadienal (4-HDDE), further oxidized into 4-hydroxy-2E-dodecadienoic acid (4-HDDA).

Figure 3 :

 3 Figure3: Reverse phase HPLC profile of monohydroxy-derivatives of linoleic (HODEs) and arachidonic (HETEs) acids of biological relevance, with detection at 235 nm. Sensitivity can be evaluated from the less represented metabolite, 15-HETE, as the signal corresponds to around 10 pmol.

Figure 4 :

 4 Figure4: Prostanoid generation from arachidonic acid. PGG 2 and PGH 2 from arachidonic acid are shown as an example of what may occur from dihomogammalinolenic and eicosapentaenoic acids, with the production of PGG 1 /H 1 and PGG 3 /H 3 , respectively. The formation of the other prostanoids also occurs from the latter two PUFA, except for PGI 1 that cannot be formed because of the lack of double bond at carbon 5. The same metabolic pathways have been reported for dihomo-arachidonic or adrenic acid with the products called dihomo-prostanoids. PG: prostaglandin; Tx: thromboxane; PG-S: PG synthase; PGH-S is a bi-functional enzyme protein with two activities: cyclooxygenase producing PGG, and hydroperoxidase converting PGG into PGH.

Figure 5 :

 5 Figure5: 5-lipoxygenase pathway from arachidonic acid leading to leukotriene (LT) A 4 , further converted into LTB 4 or LTC 4 , the first peptido-LTs. Part of 5-HpETE is reduced into 5-HETE by glutathione peroxidase. LT 5s can be obtained from eicisapentaenoic acid, and LTB 3 can be produced from 5,8,11-eicosatrienoic acid (20:3w-9 or Mead acid). Lox: lipoxygenase; H(p)ETE: hydroxy(peroxy)-eicosatetraenoic acid; GSH: reduced glutathione.

Figure 6 :

 6 Figure6: Reverse phase HPLC profiles of monohydroxy-and dihydroxy-derivatives initiated by 15/w6lipoxygenase acting upon various PUFA having at least three Z double bonds from the omega-6 position. The dihydroxy-derivatives contain an E,Z,E conjugated motif with the E,Z,E geometry. These derivatives have been collectively named poxytrins (PUFA oxygenated trienes). The monohydroxy, and dihydroxy-derivatives have been detected at 235 and 270 nm, respectively.

Figure 7 :

 7 Figure 7: Metabolic pathways of arachidonic acid through cytochrome P 450 monooxygenases. Cyt: cytochrome; EET: epoxy-eicosatrienoic acid; DHETE: dihydroxy-eicosatrienoic acid.

Acknowledgements

The support of INSERM, Ministry of Education and Research, and LISA Carnot Institute is acknowledged.


Thi Hai

Ha Dang
email: hadang@ibisc.univ-evry.fr

Guillaume Hutzler
email: guillaume.hutzler@ibisc.univ-evry.fr

Philippe Hoppenot
email: philippe.hoppenot@ibisc.univ-evry.fr

Emotion modeling for intelligent agents -Towards a unifying framework

Existing computational models of emotions, although based on different psychological theories, share common properties and may be seen as the different facets of a common emotional process. We thus present our model GRACE -aiming at unifying existing models into a single architecture while preserving the peculiarities of each of them. We also demonstrate the generality of GRACE in emulating the behavior of these existing models.

Introduction

Research in the field of computer-human interaction has long been focusing on incorporating anthropomorphic characteristics to make applications more familiar and acceptable by humans. Together with vocal and facial abilities, emotions in particular seem to have been playing a great role in reducing the distance between humans and machines [START_REF] Forlizzi | Assistive robotics and an ecology of elders living independently in their homes[END_REF] [START_REF] Tapus | User-Robot Personality Matching and Robot Behavior Adaptation for Post-Stroke Rehabilitation Therapy[END_REF].

Since the very first model of emotions in computer science, known to be the work of Simon in 1967, several applications in computer science have been developed (see [START_REF] Gratch | Evaluating a computational model of emotion[END_REF] for a more detailed list). All of theses computational models of emotions, although based on various psychological theories, share common properties and may be seen as the different facets of a common emotional process. GRACE is thus our attempt to propose a unifying framework, encapsulating and generalizing existing computational models, while preserving the peculiarities of each of them.

Section 2 of the paper presents the three main traditions in psychology along with computational models of emotions for virtual agents. In section 3, we present GRACE, which we consider as a generic framework for modeling emotional processes. To support our claim of GRACE's generality, we develop in section 4 a projection of the four models mentioned in section 2 into our model GRACE. We end the paper with some discussion about current results.

Three main approaches in the computational modeling of emotions

In this section, we review the different psychological traditions about human emotions and the main computational models that were developed on the basis of these traditions.

Basic emotions theories

Psychologists Ekman [START_REF] Ekman | An argument for basic emotions[END_REF], Izard [START_REF] Izard | Four systems for emotion activation: Cognitive and noncognitive processes[END_REF], and Johnson-Laird and Oatley [START_REF] Johnson-Laird | Basic Emotions, Rationality, and Folk Theory[END_REF] consider human emotions as being grouped into families. In these so-called basic emotions theories, human emotions are characterized by a set of basic emotions. Each basic emotion is a family of antecedent event, appraisal, physiological symptoms, behavioral response, etc. According to Izard [START_REF] Izard | Four systems for emotion activation: Cognitive and noncognitive processes[END_REF], there are four groups of events that can trigger emotional processes. These are Neural Elicitors for those concerning brain activities, neurotransmitters; Sensorimotor Elicitors concerning sensorimotor processes (like facial expression, body gesture), Motivational Elicitors concerning all motivations that lead to emotions, like thirst, hunger; Cognitive Elicitors which can be understood as beliefs, desires, intentions, social standards, etc. This theory is also quite popular in computer science as it favors taskoriented or task-specific modeling.

Cathexis [START_REF] Velásquez | Modeling emotions and other motivations in synthetic agents[END_REF][10] is an attempt at modeling emotions in terms of proto-specialists introduced by Marvin Minsky [START_REF] Minsky | The Society of Mind[END_REF] and implementing the four types of emotion elicitors proposed by Izard. This is one of the first realizations of basic emotions theories as a process. Beside Cathexis, other works also used basic emotions theories, such as the robot Kismet of Breazeal [START_REF] Oatley | Towards sociable robots[END_REF], the agent Max of Becker, Kopp, and Wachsmuth [START_REF] Becker | Simulating the Emotion Dynamics of a Multimodal Conversational Agent[END_REF], and many others.

Appraisal process -OCC's model of emotions

Psychologist Ortony and colleagues proposed in 1988 a model of emotions aiming at being implemented in computational models, well known under the name OCC (for the names of the three inventors Ortony, Clore, and Collins) [START_REF] Ortony | The Cognitive Structure of Emotions[END_REF]. To calculate an emotional response, this model considers three aspects of the perceived event, which are the consequences of the event, the action of the agent, and the aspects of objects. Thanks to its simplicity and its clarity, OCC's theory became a really popular reference for the design of agents.

One of the first applications of OCC is Affective Reasoner [START_REF] Elliott | Using the Affective Reasoner to support social simulations[END_REF] -a computational model of emotions for virtual agents. The emotional process functions on the basis of a set of predefined rules (called Emotion Eliciting Condition Relations -EECRs) to do the emotion mapping. These EECRs were constructed based on the model OCC. Beside EECRs, the emotion mapping process also depends on the personality type of the agent, which means its goals, standards, and preferences.

There is also Agent Greta [START_REF] De Rosis | From Greta's Mind to Her Face: Modelling the Dynamics of Affective States in a Conversational Embodied Agent[END_REF] -an emotional BDI agent, developed on the BDI model [START_REF] Rao | Modeling rational agents within a BDI-architecture[END_REF], which considers the influence of emotions on its actions. The affective state of Greta corresponds to the evaluation of the event using the OCC's theory. Other works also used OCC's theory as a reference model of emotions for agents, such as FLAME of El-Nars et al [START_REF] El-Nasr | FLAME -A Fuzzy Logic Adaptive Model of Emotions[END_REF], the Oz project of the MIT [START_REF] Bates | An Architecture for Action, Emotion, and Social Behaviour[END_REF].

Coping mechanism for emotional agents in stressful situation

R. S. Lazarus suggested that the emotional response of humans is produced through two stages of evaluation -cognitive appraisal and coping [START_REF] Lazarus | Emotion and adaptation[END_REF]. The cognitive appraisal is in charge of evaluating the person-environment relationship, giving the former an adapted view about the current situation. At the stage of coping, there will be problem-focused coping or emotion-focused coping. In problem-focused coping, one tries to establish a plan to solve the problem that makes the situation stressful. In emotion-focused coping, one aims at reappraising the situation by changing personal goals and intentions so as to make the situation less stressful and more acceptable.

This theory of coping has inspired the model EMA by Gratch and Marsella [START_REF] Gratch | A domain-independent Framework for Modeling Emotion[END_REF], which aims at modeling human reactions in a simulation of stressful situations. The interesting issue of this model, in our opinion, is the implementation of the coping process into a model of emotions. Emotions in EMA can be not only a factor in the planning of the agent but also a factor that urges the agent to review/change its own beliefs and desires.

Incompleteness of existing models

Physiological reactions, which are surviving mechanisms that help living organisms get out of danger, appear to be an important aspect of emotional process. However, these reactions are taken into account in Cathexis (via proto-specialists of Fear and Anger) but not present in Greta Affective Reasoner, and EMA.

We also agree with the authors of EMA that the coping mechanism plays an important role in the emotional process. Emotions are not only involved in the planning process but are also an important factor that modifies our way of thinking about reality. But again, only EMA implements this coping ability.

These lacks lead us to propose a generic model of emotions named 'GRACE' that incorporates all necessary elements in order to better reflect the emotional process in humans and other species.

GRACE -Generic Architecture to Create Emotions in software agent

According to Scherer, emotional processes of human-beings have in general five important features [START_REF] Scherer | Emotions are emergent processes: They require a dynamic computational architecture[END_REF]. One interesting aspect of the work of Scherer is that these five general features can account for the different processes and strategies proposed by the other theories, reflecting his care to take into account existing theories into his proposal. Another aspect that draws our attention to the work of Scherer is his willingness to propose a computational model of emotions. Scherer's approach is thus consistent with our objective to have a general model of emotions that can be a blueprint for further implementations of complete emotional processes in computerized applications.

Figure 1 Architecture of GRACE

In the aim of incorporating the five features proposed by Scherer, we proposed GRACE [START_REF] Dang | GRACE -Generic Robotic Architecture to Create Emotions[END_REF][23] [START_REF] Dang | How to gain emotional rewards during human-robot interaction using music? Formulation and propositions[END_REF] as a general framework of emotions for computerized applications (figure 1).

In this model, we can simulate reflex emotional reactions, like surprise or fear by transferring the event from Sensation to Body directly through Physiological Interpretation. The cognitive analysis is achieved by transferring the event through Cognitive Interpretation. Information from the two interpretation modules then goes through the module Behavior. This is the place where the model considers the current effect of the event (interpreted by Cognitive Interpretation and Physiological Interpretation), re-evaluates that perception with the current internal state and motivation to produce a more informed emotional output. This emotional output is then expressed by the Body and triggers an update in the module Knowledge (corresponding to Experience), which further entails possible modifications of the internal state.

GRACE vs. other models of emotions

This section explains how GRACE can emulate some computational models listed in the section 2.

Basic emotions theory -Cathexis

As we already underlined, the emotional process in Cathexis is triggered by four kinds of Emotion Elicitors (Neural Elicitors, Sensorimotor Elicitors, Motivational Elicitors, and Cognitive Elicitors).

Regarding the functionality of each kind of Elicitors, we propose to implement them into GRACE as following (see figure 2):

Figure 2 GRACE emulating Cathexis

As compared to Cathexis, GRACE more clearly relates the different elicitors to one another. Moreover, with GRACE, it is possible to simulate different personalities by simply changing the way each elicitor affects the emotion intensity computation in Behavior. This would significantly facilitate the design of empathetic agents with different personalities.

Appraisal process -Affective Reasoner

Affective Reasoner, as a model of mapping rules, can be easily transferred to GRACE (see figure 3). While the rule set EECRs of Affective Reasoner embraces the agent's personality, GRACE allows separating the two elements, which makes it more flexible and scalable. GRACE, in addition, can be more flexible as it allows taking into account, among other things, the current internal state or immediate reactions, which influence the final emotional output.

Coping strategy -EMA

To emulate EMA in GRACE, we propose to implement its cognitive analysis into Cognitive Interpretation. Its appraisal values represent the understanding of the agent about the perceived event, which helps the agent in establishing a plan (either emotion-focused of problem-focused plan) in response to the event. Coping process is implemented in Behavior. In GRACE, coping can account for personality defined in the module Internal State. GRACE also provides to the design of EMA the possibility, in stressful situations, to have emotional reflexes without going through the coping process. Also, re-appraisal can be seen as a kind of neural elicitation, or brain activities, where agent can generate some imaginary events to overcome the stress.

Figure 5 GRACE emulating EMA 5 Conclusion

This paper is an attempt to put together the most important features of emotional processes that should, according to us, be implemented into virtual agents. Our aim, when proposing GRACE, is not to criticize the efficiency of existing models of emotions, but to propose a unifying framework for emotion modeling. The model is not meant, however, to be implemented in its full details but was rather thought as a general and generic framework that could be instantiated partially depending on which facets of emotions a designer wishes to incorporate into his/her agents.

Figure 3 GRACE emulating Affective Reasoner 4 . 3

 343 Figure 3 GRACE emulating Affective Reasoner 4.3 Appraisal process -Greta In GRACE, the mechanism of Greta's Mind can be implemented in Cognitive Interpretation. The output of Cognitive Interpretation is the same as the output of Mind, i.e. the appropriate emotion in OCC. That output is then transferred to Behavior, where the intensity of the emotion is determined upon agent's personality. The output of Behavior is finally expressed using Body expressive features. Emotion decay occurs in Internal State component. In addition, GRACE can enrich Greta with other properties of emotions, such as physiological symptoms, imaginary events (through Intuition), etc.

Figure 4 -

 4 Figure 4 -GRACE emulating Greta


Thi-Hai-Ha Dang
email: hadang@ibisc.univ-evry.fr

Guillaume Hutzler
email: guillaume.hutzler@ibisc.univ-evry.fr

Philippe Hoppenot
email: philippe.hoppenot@ibisc.univ-evry.fr

How to gain emotional rewards during human-robot interaction using music? Formulation and propositions

In this paper, we present arguments for the need of emotion modelling and we define elements for a study in Human-Robot Interaction (HRI) using music. We also propose an adaptation of our generic model of emotions (GRACE) to give a precise idea of how to design emotional intelligence for a robot with music-related abilities.

Introduction

Thanks to the rapid development of artificial intelligence (including pattern recognition, speech processing, motion planning, machine learning, etc.), robots nowadays have gained various abilities to do repetitive, meticulous, dangerous tasks with high precision at great speed. Robots are thus greatly appreciated for replacing human in everyday/industrial/outer space tasks.

Besides, robotic applications are also present in assistive tasks, where comfortable experience is in the first rank of exigencies. This entails demands on soft skills for robots in order to deal with human mental states, to mimic human social interaction and to gain cooperation of human. Experimentations of robotic applications have also shown the emotional rewards to human in interacting with intelligent robots.

Moreover, various researches interested in music and emotions have claimed that the emotional reward of music remains an important reason for the universal appeal of music [START_REF] Gebhardt | von Georgi: Music, mental disorder and emotional reception behavior[END_REF][START_REF] Zentner | Emotions Evoked by the Sound Music: Differentiation, Classification, and Measurement[END_REF]. So, can music be involved in human-robot interaction to gain emotional rewards to human? This paper is an attempt to formulate such a case study.

In section two, we highlight the main characteristics of HRI to be taken into account in robotic application design. We present also in this section the findings of emotional rewards during human-robot interaction and recent works on emotion modelling in computer science. In section three, the first subsection lists important elements to be focused on in studying the role of emotional rewards in HRI using music. In the second subsection, we propose an adaptation of GRACE, our generic model of emotions, in order to take into account these elements. Finally, in section four, we give some perspectives for our study on emotional rewards in HRI using music.

Robotics, since the use of diverse robotic applications in human life, has been opening large perspectives of multi-disciplinary works, e.g. machinery design, control mechanism, beautiful appearance, intelligent reaction, good performance of task accomplishment and also efficient interface of interaction (i.e. humanrobot interaction -HRI). This HRI itself is also a field of study where several researches have been focused. The number of research projects on communication channels between human and robot (keyboard, mouse, voice, speech, gesture, etc.) increases with time. Many of them have gained successes in improving communication efficiency. Soon afterward, importance of the emotional aspect of communication began to draw attention. Researchers in HRI started to study/apply knowledge in emotion theories into their applications. In this section, we will present advances in HRI research and the need of emotion modelling in robotic applications.

Human-Robot Interaction -Need of emotion modelling

Natural and affective interaction -this could be considered as the principle goal of all works on human-machine interaction. Being a specific branch of research in human-machine interaction, research in HRI, specially HRI for robotic assistance, is focused on such a kind of interaction. Kiesler and Hinds, in the introduction of Special Issues in Human-Robot Interaction [START_REF][END_REF], mention some important points to be considered in HRI research, particularly in robotic assistance, which are:

-Anthropomorphic model : this relates to the need of using natural channels of communication to facilitate the interaction of humans with the robot.

The more natural human-robot communication is, the better the robot gains cooperation from human. Several channels could be considered such as voice, speech, gesture, etc. -Physically close to other robots, people, objects: usually, robotic assistance involves human participation in a close distance. This characteristic has to be seen as of major importance to study HRI. This importance appears in the necessity of robots having abilities to work in a dynamic and challenging environment. This interaction also entails the question of safety issues for both human and robot during interaction. -Sufficient knowledge of the context of use: this is for the design of an effective interaction scheme to provide appropriate HRI applications that suit the social/ethnical standards of the user. This also includes the robot's abilities to automatically learn about themselves and their world to well adapt to changing situations.

Experimentations with intelligent robots have repeatedly demonstrated the importance of emotional rewards in robotic applications. Since 2003, several experiments with Paro [START_REF]Official site for Mental Commitment Robot[END_REF], a seal-like robot, held in Japan, in United States of America, in France, etc., have shown that people (elderly, children) experienced some kind of comfortable feeling while interacting and living with robot Paro. In a survey of [START_REF] Forlizzi | Assistive robotics and an ecology of elders living independently in their homes[END_REF] [START_REF] Forlizzi | Assistive robotics and an ecology of elders living independently in their homes[END_REF], emotional rewards were mentioned many times when people talked about the capacity of assistive materials to permanently keep elderly in connection with others (their relatives, their friends, their neighbours, etc.).

Being encouraged from works of psychologists [START_REF] Darwin | The expression of the emotion in man and animals[END_REF][START_REF] Gebhardt | von Georgi: Music, mental disorder and emotional reception behavior[END_REF][START_REF] Scherer | Emotion, the psychological structure of[END_REF]5] on the importance of emotional reactions in the adaptation abilities of human, research in artificial intelligence has begun studying the emotional intelligence of human and trying to implement it in robots to render more natural behaviour and so to gain success in interaction with human.

Researches in Modelling Emotions For Computerized Applications

Advances in robotics and artificial intelligence have strongly changed the way people experience computerized applications. With the discovery of emotional rewards for human during interaction with his robotic equipments, research on modelling emotions has come to life. This research direction requires to be based on studies on emotion theories in psychology. Among several theories, some are often used as basis by almost all works. These theories are the theory of Ortony et al [START_REF] Ortony | The cognitive structure of emotions[END_REF] on event appraisal, the theory of Lazarus [START_REF] Lazarus | Emotion and Adaptation[END_REF] about appraisal and coping, and the theory of Scherer [START_REF] Scherer | Emotion, the psychological structure of[END_REF][START_REF] Scherer | On the nature and function of emotion: A component process approach[END_REF] about emotional processes.

In the aim of modelling and implementing emotions into computerized application, several approaches have been explored. El-Nasr et al. [START_REF] El-Nasr | Ioerger: FLAME -A Fuzzy Logic Adaptive Model of Emotions[END_REF] explored the use of fuzzy logic to calculate emotional intensities. This model was implemented in a simulation named PETEEI -a Pet with Evolving Emotional Intelligence. Embodied Conversational Agents, such as ParleE (by [START_REF] Bui | ParleE: An Adaptive Plan Based Event Appraisal Model of Emotions[END_REF] [START_REF] Bui | ParleE: An Adaptive Plan Based Event Appraisal Model of Emotions[END_REF], Greta (Pelachaud, 2003) [START_REF] De Rosis | From Gretas Mind to her Face: modelling the Dynamics of Affective States in a Conversational Embodied Agent[END_REF], GALAAD (by C. Adam, 2005) [START_REF] Adam | Emotions: from psychological theories to logical formalization and implementation in a BDI agent[END_REF], have been developed, taking into account the appraisal process of emotions and even personality traits. There is also a framework for training applications called Mission Rehearsal Exercise, developed by Gratch and Marcella in 2004 [START_REF] Gratch | The Architectural Role of Emotion in Cognitive Systems[END_REF], that provides virtual scenarios for health intervention, marketing and entertainment. The emotional aspect of this proposition is the adaptation of event appraisal and coping (two important aspects of emotions, according to Lazarus [START_REF] Lazarus | Emotion and Adaptation[END_REF]). Now come the applications in the domain of HRI. In 2002, C. Breazeal at MIT lab attempted to establish an interaction between a robot and a human inspired by the relation between a baby and his/her parents [START_REF]Emotion and sociable humanoid robots[END_REF]. Emotions experimented are anger, distaste, fear, sadness, and happinesss. Recently, MIT lab has also announced their advance in modelling emotions into another robot named Nexi -a Mobile Dexterous Social Robot with more human-centric communication and interaction abilities [START_REF]Project MDS: Personal Robots Group[END_REF]. In 2008, we proposed a model of emotions for robotic applications, called GRACE [START_REF] Dang | Comparison of recent architectures of emotions[END_REF][START_REF] Dang | Experimentation with GRACE, the Generic Model of Emotions For Computational Applications[END_REF]. GRACE is based on the psychological theories that we mentioned previously. This model is considered generic as it can instantiate anterior models of emotions in the domain of computer science by using some of its components. A comparison of these models of emotions can be found in [START_REF] Dang | Comparison of recent architectures of emotions[END_REF].

Proposition for the Study on Emotional Rewards of HRI Using Music

The universal appeal of music mostly comes from the emotional reward that music offers to its listeners. According to Zentner in [START_REF] Zentner | Emotions Evoked by the Sound Music: Differentiation, Classification, and Measurement[END_REF], emotional impacts of music explain its prominent role in people's everyday lives. In [START_REF] Gebhardt | von Georgi: Music, mental disorder and emotional reception behavior[END_REF], music is ascertained to have direct effects to treatment of emotion disorders. Now imagine that a robot can detect the emotional state of its partner when he listens to the music (or he plays the music) and then react in an appropriate manner, can this robot enhance emotional rewards for its partner during/after this interaction? This is the question we try to answer in our study.

To this end, beside the study of emotional impacts of music to human, we need a model of emotions to implement the robot's behaviour so as to enable an emotionally rich interaction with humans. With GRACE, we already described and modelled the emotional process for robots. However, the use of music during human-robot interaction requires some specific adaptations. In this section we study these requirements and propose modifications to suit the case of using music.

Important Elements in Studying HRI Using Music

Given that the objective of our work is to study the mental rewards during human-robot interaction using music, the interaction scenario will be between a musician and a robot with music-related abilities. Thus, the desired robot should be equipped with a cognitive process that helps it to perceive its partner's emotions via his music (and maybe his visible behaviour) and then to react in an appropriate way. Important information used for the robot's intelligent behaviour in this study should be:

-The emotional interpretation of the music played by a musician: this represents the interpretation done by the musician. The robot has to know how its partner feels via his music. As claimed in [START_REF] Zentner | Emotions Evoked by the Sound Music: Differentiation, Classification, and Measurement[END_REF]5], it is feasible to construct a learning module that takes a piece of music as input and reproduces the emotions probably experienced by a specific individual. -The robot's mental state (specifically personality trait of robot: introversive, extroversive, aggressive, curious, etc.). This mental state has an important impact on the selection of action alternatives. So, with the same emotional state of musician (interpreted by the module of interpretation), the robot can change its preferred tendencies of action based on its own mental state. For example, if the detected emotion of the musician is sadness, then the robot, if it is extroversive, can have action tendencies to show empathy to its partner; but if the robot is aggressive, it could have tendency to disturb the musician (like making noises, performing exciting movements, etc.). -The robot's current goals (excite its partner, get rid of him, disturbing ongoing situation, etc.). This current goal gives hints to select preferred tendencies. When the robot wants to excite its partner, it should prefer actions that decrease the negative emotions of its partner and increase positive impact (such as performing funny dances, showing empathy by dancing coherently with the music played by its partner). -Ability to guess future actions of its partner (the musician). This ability allows the robot to predict upcoming events to better adapt its behaviour. This could be done by deploying a prediction module based on a predefined scenario or applying pattern recognition on the set of events over time.

To take into account important elements presented above, the model GRACE, which remains abstract for general use case, will be adjusted and concretized into an adaptation called MACE-GRACE (Music-adapted Architecture to Create Emotions -GRACE).

MACE-GRACE -Adaptation of GRACE in the Context of HRI Using Music

Being constructed by merging psychological theories on emotions and computational models of emotions, GRACE aims at simulating an entire emotional process of human. The emotional process of GRACE is described in Figure 1: For our study on emotional rewards in HRI application, some details must be precisely pointed out:

!"#$%"&'()*+"!#!,$(-#&#$( -$
-Sensation: this GRACE module does the collection of information necessary for emotional processing. In our case, information collected will be musical stimuli (let's call them 'musical events') and environment information required for robot safety verification. -Physiological interpretation: this GRACE module simulates the physiological symptoms of the emotional process (like temperature sensations, respiratory and cardiovascular accelerations and decelerations, trembling and muscle spasms) to show reflex reactions. In our case, only the safety of the human and the robot will be verified in this module to ensure a good performance of the robot; other physiological symptoms will not be taken into account. -Cognitive interpretation: in GRACE, this module is in charge of analyzing the emotion-eliciting events to have a global view of current situation (more specifically, the emotion-related aspects of the current situation). Result of this analysis is used to decide action tendencies to respond to events. In our case, what we want in output of this module are the current emotions expressed in the musical events played by the robot's partner. We tend to use a learning module here to analyze the musician's emotions via his music. The result of this module can help the robot to choose appropriate actions to respond to the musical events. -Mood : this component contains the robot's current mental state, its current stand faced to current situation, its current goals/needs in term of assigned tasks and even its habits(e.g. its preferred action tendencies, preferred situation/events). To focus on the goal of our study (mental rewards for humans during the interaction with a robot having music-related abilities), we will code in this module only the robot's goal/need, and a simplified mental state. -Feeling: this module plays the role of memorization of the running process. It captures the reaction of all components to events and is used as a database for adaptation and memory recall abilities. In our study, this part would store the musician's attitude (history/personality of musician). This module can help the robot to decide its preferred strategy of action in regard to the attitude of musician. -Intuition: this module of GRACE is supposed to do the anticipation of future events, to generate imaginary events (i.e. internal events). We propose in our study to relate it to some memory recalls or anticipations of future reaction of the musician so that the robot can react rapidly to new events. -Behaviour : in GRACE, this module is for the selection of an action to respond to an event. This selection takes into account information of both the external situation and the robot's internal state. In our case study, this would be the selection/regulation of robot's movements (accelerate or decelerate the movement of the hands, turn the body around, make some noise, etc.) according to the musician's emotion expressed in his music. -Body: this module executes the robot's motor expression in response to input events. Generally, this includes reflex reactions from Physiological Interpretation and reactions selected by Behaviour module. In our study, this part will execute the robot's music-related reactions to show that the robot understands the emotional expression in the playing music and shows its response in relation with the current situation and its own personality.

To take into account these issues, we propose MACE-GRACE, an adaptation of GRACE for HRI using music, described in Figure 2. The analysis process of robot's intelligence when it captures musical events will be as follows:

The physiological interpretation is in charge of verifying the physical impact of events to robot to ensure the robot's safety based on its position, its battery level, etc. The simulation of human-like physiological stimuli (such as body temperature, heart beat, respiration, or even reflex reactions) of the robot is not taken into account in this study. The robot's cognitive interpretation will measure its partner's emotions based on the music that he plays. This measured value will be sent to the action selection to determine an appropriate response to the current emotional state (measured by the cognitive interpretation). This phase of selection is also influenced by the robot's current goal/need (calm down Fig. 2. Adaptation of GRACE in the context of HRI using the musician or excite the current ambiance, amplify the current state of the musician) and with regard to the musician's attitude. The execution of the selected action will be stored in the robot's database along with the behaviour of the musician at that moment. This memory allows the robot to make predictions of the future behaviour/music of its partner so that it could have an adaptive behaviour in the whole conversation.

With this adaptation, we tend to integrate in our robot an artificial emotional intelligence to interact with a musician. This emotional intelligence can help the robot not only perceive the emotional state of its partner but also react to diverse situations in an adaptive manner in term of emotion-based strategy.

Conclusion

Throughout the paper we have presented major characteristics of robotic applications and then mentioned the role of emotional rewards during human-robot interaction. We also proposed an evolution of our model of emotions previously proposed for computerized applications. This evolution is in the aim of studying the emotional rewards in human-robot interaction using music. We consider the study as an attempt to implement emotional intelligence into a concrete computerized application, more specifically, an assistive/entertaining robot. This could have interesting applications in such fields as music-based therapy, emotionbased movements, preferences of robot abilities, etc.

The realization of this adaptation is a long way to go. It consists in implementing the model, doing the experimentation with the robot and then analyzing the results. Furthermore, the implementation of the model could also concern works in database construction (including choices in music style, user preferences, music descriptors, definition of musical event, etc.), system specification, programming tasks, etc. Yet, we demonstrated in this study that GRACE, our generic model of emotions, was versatile enough to be easily adapted to the specific needs of emotional interaction using music.

Fig. 1 .

 1 Fig. 1. Model GRACE -Generic Architecture to Create Emotions


Quentin Mérigot

A multiscale approach to optimal transport

In this paper, we propose an improvement of an algorithm of Aurenhammer, Hoffmann and Aronov to find a least square matching between a probability density and finite set of sites with mass constraints, in the Euclidean plane. Our algorithm exploits the multiscale nature of this optimal transport problem. We iteratively simplify the target using Lloyd's algorithm, and use the solution of the simplified problem as a rough initial solution to the more complex one. This approach allows for fast estimation of distances between measures related to optimal transport (known as Earth-mover or Wasserstein distances). We also discuss the implementation of these algorithms, and compare the original one to its multiscale counterpart.

Introduction

Engineer and mathematician Gaspard Monge proposed the following problem [START_REF] Monge | Mémoire sur la théorie des déblais et de remblais[END_REF]: what is the cheapest way to transport a pile of sand into a hole with minimum cost, knowing that moving an individual particle from a position x to another position y has a cost c(x, y) ? This problem gave birth to the field of optimal transport, which has been very vivid in the past twenty years, with applications in geometry, probability and PDEs (see e.g. [START_REF] Villani | Optimal transport: old and new[END_REF]).

However, Monge's problem was an engineering problem, and it is not very surprising that various form of optimal transport appeared in many applied fields. In computer vision, distances defined using optimal transport have been used as a way to compare the color histograms of images in [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF] under the name of Earth mover distance. Optimal transport on the circle has been used for transferring the hue of an image to another [START_REF] Delon | Fast Transport Optimization for Monge Costs on the Circle[END_REF]. In combination with Lloyd's algorithm, optimal transport seems a interesting tool for optimal quantization [START_REF] Balzer | Capacity-constrained point distributions: a variant of Lloyd's method[END_REF]. More recently, methods inspired by (or relying on) optimal transport have been proposed as a tool in several parts of geometry processing: surface comparison [START_REF] Lipman | Conformal Wasserstein distances: comparing surfaces in polynomial time[END_REF], surface reconstruction from data corrupted with outliers [CCSM10, MDGD * 10], construction of optimized primal-dual triangulations [START_REF] Mullen | Hodge-Optimized Triangulations[END_REF], reconstruction with sharp corners and edges [START_REF] De Goes | An optimal transport approach to robust reconstruction and simplification of 2D shapes[END_REF].

Yet, the lack of a practical method to compute optimal transports maps except in 1D has hindered the development of many of these applications. Even in the simplest planar situation, namely with c(x, y) = x -y 2 , there is a lack of reasonably fast and widely usable method able to efficiently compute optimal transport maps.

1.1. L 2 optimal transport. In the remaining of the paper, we will deal with the L 2 optimal transport, i.e. where the cost for moving a certain amount of mass from a point x to a point y is proportional to square of the Euclidean distance x -y 2 . This model is well-understood theoretically, and has several nice properties (such as uniqueness of the solution) that follow from strict concavity of the cost. Numerical schemes have been proposed by Brenier-Benamou [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF], Loeper [START_REF] Loeper | Numerical solution of the Monge-Ampère equation by a Newton's algorithm[END_REF] and Angenent-Haker-Tannenbaum [START_REF] Angenent | Minimizing flows for the monge-kantorovich problem[END_REF] to solve L 2 optimal transport. However, numerical instabilities makes them difficult to use for general problems: for instance, [START_REF] Loeper | Numerical solution of the Monge-Ampère equation by a Newton's algorithm[END_REF] requires a lower bound on the density of the source measure, while the gradient descent algorithm of [START_REF] Angenent | Minimizing flows for the monge-kantorovich problem[END_REF] suffers from a drift effect which produces optimal maps that are not transport plans. Another possibility to find optimal transport plans is by discretizing the source and/or target measure. These discrete approach include linear programming, the Hungarian method, and a variant known as Bertsekas' auction algorithm [START_REF] Bertsekas | The auction algorithm: A distributed relaxation method for the assignment problem[END_REF]. These methods work for general cost functions, and are often unable to take advantage of the geometric simplifications that occur when working specifically with the squared Euclidean metric.

A promising approach, both from the practical and the theoretical point of view has been proposed in [START_REF] Aurenhammer | Minkowski-type theorems and least-squares clustering[END_REF]. In this approach, the source measure has a density ρ while the target measure is given by a sum of Dirac masses supported on a finite set S. The fact that the source measure has a density ensures the existence and uniqueness of the optimal transport map. In this case, solving the optimal transport problem amounts to finding a weight vector (w p ) p∈S such that the power diagram of (S, w) has the following property: for every point s in S, the proportion of the mass of ρ contained in the corresponding power cell should be equal to the mass of the Dirac mass at s (see Section 2.2). In the same article, the authors also introduced a convex function whose minimum is attained at this optimal weight vector, thus transforming the optimal transport problem into an unconstrained convex optimization problem on R N where N is the number of points in the set S.

1.2. Contributions. We revisit the approach of [START_REF] Aurenhammer | Minkowski-type theorems and least-squares clustering[END_REF] with implementation in mind. Our main contribution is a multiscale approach for solving the unconstrained convex minimization problem introduced in [AHA98], and thus to solve L 2 optimal transport.

Let us sketch briefly the main idea. In order to solve the optimal transport problem between a measure with density µ such as a grayscale image and a discrete measure ν, we build a sequence ν 0 = ν, . . . , ν L of simplifications of ν using Lloyd's algorithm. We start by solving the much easier transport problem between µ and the roughest measure ν L using standard convex optimization techniques. Then, we use the solution of this problem to build an initial guess for the optimal transport between µ and ν L-1 . We then proceed to convex optimization starting from this guess to solve the optimal transport between µ and ν L-1 . This is repeated until we have obtained a solution to the original problem. If the original target measure has a density, we use Lloyd's algorithm to obtain the first discretization ν. Note that this idea was independently proposed in [START_REF] Bosc | Numerical approximation of optimal transport maps[END_REF].

This procedure provides a significant speedup, up to an order of magnitude, for computing optimal transport plans. Moreover, at every step of the algorithm it is possible to obtain a lower and upper bound on the Wasserstein distance (also known as Earth-mover distances [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF]) between the source measure µ and the original target measure ν. Using this approach, one can obtain rough estimates of Wasserstein distances between two images with a speedup of up to two order of magnitude over the simple convex optimization approach.

Background

We briefly recap the few concepts of measure theory and optimal transport that we use, before explaining the relation between the L 2 optimal transport and power diagrams. We also recall how optimal transport can be turned into an unconstrained convex optimization problem.

2.1. Measure theory and Optimal transport. A non-negative measure µ on the space R d is a map from (measurable) subsets of R d to a non-negative numbers, which is additive in the sense that µ (∪ i∈N B i ) = i µ(B i ) whenever (B i ) is a countable family of disjoint (measurable) subsets. The total mass of a measure µ is mass(µ) := µ(R d ). A measure µ with unit total mass is called a probability measure. The support of a measure µ, denoted by spt(µ) is the smallest closed set whose complement has zero measure.

The optimal transport problem involves two probability measures: a source measure µ, and a target measure ν. We will always suppose that the source measure µ has a density, i.e. there exists a non-negative function ρ on R d such that for every (measurable) subset

B of R d , µ(B) := B ρ(x)dx.
On the other hand, we will assume that the target measure ν is discrete, supported on a finite set S of R d . This means that there exists a family of positive coefficients (λ p ) p∈S such that for every subset B,

ν(B) = p∈S∩B λ p .
The above formula is equivalent to writing ν as the sum p∈S λ p δ p , where δ p is the unit Dirac mass at p. The integral of a continuous function φ with respect to these two measures is R d φ(x)dµ(x) := R d φ(x)ρ(x)dx, and

R d φ(x)dν(x) := p∈S λ p φ(p).
In Section 3.1 we will see how to adapt the proposed method to the case where the source and target measures both have density.

Transport map. The pushforward of a measure µ by a map T : R d → R d is another measure T # µ defined by the equation T # µ(B) := µ(T -1 (B)) for every subset B of R d . A map T is called a transport map between µ and ν if the pushforward of µ by T is ν. We denote by Π(µ, ν) the set of transport maps between µ and ν.

For instance, a map T is a transport map between the source measure µ and the target measure ν described at the beginning of this paragraph if and only if for every point p in the support S of ν,

λ p = µ(T -1 ({p})) = T -1 ({p}) ρ(x)dx .
Optimal transport maps. The cost of a transport map T between the source measure µ with density ρ and the target measure ν is defined by:

c(T ) := R d x -T (x) 2 dµ(x) = R d x -T (x) 2 ρ(x)dx
The problem of optimal transport, also called Monge's problem, is to find a transport map T opt whose cost is minimal among all transport maps between µ and ν, i.e.

T opt := arg min{c(T ); T ∈ Π(µ, ν)}

The non-convexity of both the cost function c and of the set of transport plans makes it difficult, in general, to prove the existence of a minimum. However, in this specific case where µ has a density and the cost is the squared Euclidean norm the existence follows from [START_REF] Brenier | Polar factorization and monotone rearrangement of vectorvalued functions[END_REF].

Wasserstein distance. The Wasserstein distance between two probability measures µ and ν, µ having a density, is the square root of the minimal transport cost. We denote it by Wass 2 (µ, ν). Intuitively, the Wasserstein distance measures the minimum global cost of transporting every bit of µ onto ν, supposing that moving an infinitesimal amount dµ(x) from x to y is equal to x -y 2 dµ(x). Note that our definition above is not symmetric, as it requires the source measure to have a density. However, this restriction can be leveraged using the notion of transport plan instead of transport map, leading to a definition of Wasserstein distance between any pair of probability measures [Vil09, Ch. 6]).

2.2. Power diagrams. Let S be a finite set of points in R d , and w : S → R be a given weight vector. The power diagram or weighted Voronoi diagram of (S, w) is a decomposition of the ambient space in a finite number of cells, one for each point in S, defined by the property that a point x belongs to Vor w S (p) iff for every q in S, one has x -p 2 -w(p) ≤ x -q 2 -w(q). Note that if the weights are all zero, this coincides with the usual Voronoi diagram.

Given such a decomposition, we will consider the application T w S which maps every point x in the power cell Vor w S (p) to its center p. This map is well-defined everywhere except on the boundary of the power cells, but since this set has zero Lebesgue measure this has no consequence for us. The pushforward of the source measure µ by the map T w S is a sum of Dirac masses centered at every point p in the set S, whose mass is the µ-mass of the corresponding power cell:

T w S | # µ = p∈S µ(Vor w S (p))δ p ,
In [START_REF] Aurenhammer | Minkowski-type theorems and least-squares clustering[END_REF], the following theorem was proven. Note that this result, as well as Theorem 2, can be also obtained as a consequence of Brenier theorem [START_REF] Brenier | Polar factorization and monotone rearrangement of vectorvalued functions[END_REF].

Theorem 1. For any probability measure µ with density, and a weighted set of points (S, w), the map T w S is an optimal transport map between the measure µ and the pushforward measure T w S | # µ. Consequently,

Wass 2 µ, T w S | # µ =   p∈S Vor w S (p) x -p 2 ρ(x)dx   1/2 .
This theorem gives examples of optimal transport plans and the discrete probability measures supported on S that can be written as the pushforward T w S | # µ, where w is a set of weights on S. As we will see in the next section (Theorem 2), it turns out that every probability measure supported on S can be written in this way. Said otherwise, the optimal transport maps between µ and a measure supported on S can always be written as T S,w for some weight vector w.

2.

3. An unconstrained convex optimization problem. Adapted weight vectors. Let µ be a probability measure with density ρ, and ν a discrete probability measure supported on a finite set S, with ν = p∈S λ p δ p . A weight vector w : S → R on S is called adapted to the couple of measures (µ, ν) if for every site p in S, one has

λ p = µ(Vor w S (p)) = Vor w S (p) ρ(x)dx . (1)
By Theorem 1, if the weigh vector w is adapted to the couple (µ, ν) then the optimal transport between µ and the discrete measure ν is given by the map T S,w . Moreover, Theorem 2 below asserts that finding a weight vector adapted to the couple (µ, ν) amounts to finding a global minimum of the function Φ below, thus turning the very constrained original problem (minimization among convex maps) into an unconstrained convex optimization problem. Note that this theorem follows from the discussion in Section 5 of [START_REF] Aurenhammer | Minkowski-type theorems and least-squares clustering[END_REF].

Φ(w) := p∈S λ p w(p) - Vor w S (s) ( x -s 2 -w(p))dµ(x) (2)
Theorem 2. Given a measure µ with density ρ on R d , and ν = s∈S λ p δ p , the following three statements are equivalent: (i) the power map T w S realizes an optimal transport between the measures µ and ν ; (ii) w is adapted to the couple (µ, ν); (iii) w is a global minimizer of the convex function Φ.

We will recall briefly how the value of the gradient of the function Φ can be computed at a weight vector w. Incidentally, the steps necessary to this computation almost sketch a complete proof of the theorem. Consider the map

Ψ(w) := p∈S Vor w S (p) ( x -p 2 -w(p))dµ(x)
which is related to the our function by the equation Φ(w) = p∈S λ p w(p) -Ψ(w). The map Ψ is concave, as we will show by writing it as an infimum of linear functions. Consider any map T from R d to the finite set S. By definition of the power cell, for every point x in Vor w S (p) one has

x -p 2 -w(p) ≤ x -T (x) 2 -w(T (x))
As a consequence, the function Ψ can be rewritten as

Ψ(w) = inf T Ψ T (w), where Ψ T (w) := R 2 ( x -T (x) 2 -w(T (x)))dµ(x).
Since the functions Ψ T all depend linearly on w, the function Ψ is concave, and Φ is convex. Moreover, it is easy to check that the values of the functions Ψ and Ψ T w S coincide for the weight vector w. Consequently, their gradient coincide to at that point, and a simple computation shows :

∂Ψ ∂w(p) (w) = Vor w S (p) ρ(x)dx (3) ∂Φ ∂w(p) (w) = λ p - Vor w S (p) ρ(x)dx (4)
In particular, the second equation shows that the gradient ∇Φ vanishes at a weight vector w if and only if w satisfies Eq. (1).

Multiscale Approach for Minimizing Φ

The efficiency of an optimization technique relies on two important choices. The most important one is the choice of descent algorithm, as it is well-known that the difference in efficiency between (for instance) the first order simple gradient descent algorithm and the second order Newton methods can be tremendous [START_REF] Fletcher | Practical methods of optimization[END_REF].

The second one is the choice of the position from where the optimization is started. Its importance shouldn't be disregarded, even for convex optimization, as the second-order convergence in Newton's descent does only happen in a small basin around the global minimizer.

In this section, we introduce our multiscale algorithm for finding a global minimum of Φ. We start by building a decomposition of the target measure ν, i.e. a sequence of discrete measures ν 0 := ν, ν 1 , . . . , ν L that are simpler and simpler as L increases. The level ℓ of the decomposition is then used to construct a good initial weight vector for the optimal transport problem (µ, ν ℓ-1 ), in a hierarchical way.

3.1. Decomposition of the target measure. A decomposition of the target measure ν = p∈S λ p δ p is a sequence of discrete probability measures (ν ℓ ) ℓ≥0 such that ν 0 = ν and ν ℓ is supported on a set S ℓ :

ν ℓ = p∈S ℓ λ p,ℓ δ p .
Moreover, for every level ℓ we are given a transport map between ν ℓ and ν ℓ+1 , that is a map π ℓ from the support S ℓ of ν ℓ to the (smaller) support S ℓ+1 of ν ℓ+1 with the additional property that for every point p in S ℓ+1 ,

λ p,ℓ+1 = q∈π -1 ℓ (p) λ q,ℓ (5)
The decomposition that we will consider in practice are constructed using Lloyd's algorithm, as explained in Section 4.2. This means that the transport map π ℓ maps every point p in S ℓ to its nearest neighbor in S ℓ+1 . We remark that having access to such a transport map between ν ℓ and ν ℓ+1 allows to bound the Wasserstein distance between these two measures. By considering the composition of transport maps, it is also possible to bound the distance between e.g. ν and ν L . Letting π = π L-1 • . . . • π 0 one has:

Wass 2 (ν, ν L ) ≤   p∈S λ p p -π(p) 2   1/2 (6)
3.2. Algorithm. We are given the source probability measure µ, and a decomposition (ν ℓ ) 0≤ℓ≤L with L levels of the target measure. The goal is to use the solution of the optimal transport problem from µ to ν ℓ+1 at level ℓ + 1 to construct the initial weight vector for the optimal transport problem between µ and ν ℓ at level ℓ. As before, we will consider the weight vectors at level ℓ as functions from the support S ℓ of µ ℓ to R. The function Φ ℓ that we optimize at step ℓ is given by the same formula as in Eq. (2).

The multiscale descent algorithm is summarized in Algorithm 1. Note that the algorithm does not depend on the choice of the convex optimization scheme (L-BFGS), which we will discuss later in Section 4.

Algorithm 1 Multiscale minimization of Φ := Φ 0 w L := 0 for ℓ = L -1 to 0 do set w ℓ,0 (p) := w ℓ+1 (π ℓ (p)) for every p ∈ S ℓ k := 0 repeat compute w ℓ,k+1 from w ℓ,k using L-BFGS on Φ ℓ set v k+1 := ∇Φ ℓ (w ℓ,k+1 ), k := k + 1 until v k q > ε set w ℓ := w ℓ,k end for
In the stopping condition . q denotes the usual L q -norm where q > 1 or q = +∞. In particular,

∇Φ(w) ∞ = sup p∈S |λ p -µ(Vor w S (p))| ∇Φ(w) 1 = p∈S |λ p -µ(Vor w S (p))|
The first quantity measures the maximum error that has been made by considering the weight vector w instead of the optimal one. In particular, if ∇Φ(w) ≥ min p∈S λ p , then one is sure that all the cells in the power diagram of (S, w) intersect µ non-trivially. This is important especially for visualization purpose, as the existence of cells with zero µ-mass lead to black spots in the interpolated pictures (see Section 5.3). The choice of . 1 plays a different role which we describe in the next paragraph.

3.3. Computation of Wasserstein distances. Simple lower and upper bounds on Wasserstein distance can be obtained at every step of the multiscale algorithm, using the fact that ∇Φ(w) 1 corresponds to the twice the amount of mass that has been misaffected. This follows from the following proposition:

Proposition 1. Let w be the a weight vector on S, and consider the image measure ν := T S,w | # µ. Then,

Wass 2 (ν, ν) ≤ D × ∇Φ(w) 2/1
1 , where D is the diameter of the support S of ν. As a consequence of this proposition, stopping Algorithm 1 at level ℓ with weight vector w ℓ yields the following estimation of Wass 2 (µ, ν):

Wass 2 (µ, ν) - R d x -T S ℓ ,w ℓ (x) 2 ρ(x)dx 1/2 ≤ D × ∇Φ ℓ (w ℓ ) 1/2 1 + Wass 2 (ν, ν ℓ ) (7)
Said otherwise, if one wants to compute the Wasserstein distance between µ and ν up to a certain error ε, it is not necessary to consider the levels of the decomposition below the first level ℓ 0 such that Wass 2 (ν ℓ 0 , ν) < ε. Note that this quantity can be estimated thanks to Eq. ( 6). The effectiveness of this approach is discussed in Section 5.2.

Proof of Eq. (7). By Theorem 1, the map T S ℓ ,w ℓ is an optimal transport between the measure µ and ν = T S ℓ ,w ℓ | # µ. This means that the Wasserstein distance Wass 2 (µ, ν) is equal to c ℓ . Moreover, by the reverse triangle inequality,

| Wass 2 (µ, ν) -Wass 2 (µ, ν)| ≤ Wass 2 (ν, ν ℓ ) + Wass 2 (ν ℓ , ν).
One then concludes using the previous proposition.

Convergence of Optimal Transport Maps.

In this paragraph, we discuss the soundness of constructing an initial weight vector for the optimal transport problem (µ, ν ℓ ) from an adapted weight vector for the problem (µ, ν ℓ+1 ). The result of this section are summarized in the following theorem. Note that the definition of zero-mean convex potential is given below, and is necessary to define uniquely the adapted weight vector: without this, adapted weight vectors are defined up to an additive constant.

Theorem 3. Let ν and (ν n ) n≥1 be discrete probability measures supported on finite sets S and (S n ) n≥1 respectively, such that lim n Wass 2 (ν, ν n ) = 0. Let:

w n :S n → R be adapted to (µ, ν n ) w :S → R be adapted to (µ, ν)
Suppose that both weight vectors yield zero-mean convex potentials (see below), and that the assumptions of Proposition 2 are satisfied. Then, for every sequence of points p n ∈ S n converging to a point p in S, one has w(p) = lim w n (p n ).

Before proving this theorem, we need to introduce a few definitions and auxiliary results.

Convex potential. Let ν be a probability measure supported on a finite set S, and let w denote the weight vector adapted to the optimal transport problem between µ and ν. Set

φ w S (x) := 1 2 x 2 -min p∈S x -p 2 -w(p) (8) = max p∈S x|p + 1 2 (w(p) -p 2 ) (9)
where v|w = i v i w i denotes the usual Euclidean scalar product. From these two formulations, it is easy to see that the function φ w S is convex, and that its gradient ∇φ w S coincides with the transport map T w S . We call such a function a convex potential for the optimal transport plan. Since adding a constant to the weight vector (or to φ w S ) does not change the transport plan, we consider the zero-mean convex potential which is uniquely defined by the extra assumption that R d φ w S (x)ρ(x)dx = 0. Proposition 2. Let µ and ν be two probability measures, µ having density ρ and ν supported on a finite set S. Let (ν n ) be a sequence of probability measures supported on finite sets (S n ), s.t. lim n Wass 2 (ν n , ν) = 0. Assume that:

(i) the support of ρ is the closure of a connected open set Ω with regular (piecewise C 1 ) boundary ; (ii) there exists a positive constant m such that ρ ≥ m on Ω; (iii) the support of all the measures ν n is contained in a fixed ball B(0, L) ; Denote φ := φ w S (resp. φ n := φ wn Sn ) the zero-mean convex potential of the optimal transport between µ and ν (resp. ν n ). Then, φ n converges to φ uniformly on Ω as n grows to infinity. This proposition is similar to [Bos10, Theorem 1], but without the requirement that the source and target measure have to be supported on convex sets. It relies on the following result (cf [START_REF] Villani | Optimal transport: old and new[END_REF], Corollary 5.23): Fact. Let ν n be a sequence of measures converging to µ, and T n (resp. T ) denotes the optimal optimal transport map between µ and ν n (resp. ν). Then, for every positive ε,

lim n→+∞ µ(∆ ε (T, T n )) = 0 ( 10
)
where

∆ ε (T, T n ) := {x ∈ R d ; T (x) -T n (x) ≥ ε}.
Proof of Proposition 2. For almost every x in Ω, the gradient ∇φ n (x) = T n (x) is included in the support of ν n , hence in the ball B(0, L) by (iii). The same holds for T , so that the inequality T -T n ≤ 2L holds for almost every x in Ω. For every p ≥ 1,

T -T n L p (µ) := Ω T (x) -T n (x) p ρ(x)dx = Ω\∆ε(T,Tn) T (x) -T n (x) p ρ(x)dx + ∆ε(T,Tn) T (x) -T n (x) p ρ(x)dx ≤ ε p + (2L) p µ(∆ ε (T, T n ))
Using Eq. (10), we obtain the convergence of T n to T in the L p (µ)-sense. Thanks to the the assumptions (i) and (ii), we can apply the Poincaré inequality on the domain (Ω, µ) to the zero-mean potentials φ n and φ to get:

φ n -φ L p (µ) := R d φ n (x) -φ(x) ρ(x)dx ≤ const(Ω, ρ) T n -T L p (µ) .
In other words, φ n converges to φ in the L p (µ)-sense. Since the support of all target measures is contained in a ball of size L, ∇φ n ≤ L, and φ w S is L-Lipschitz. Hence, φ n also converges uniformly to φ on the support of µ.

Proof of Theorem 3. We begin by proving that there for every point p in S, there exists a sequence of points q n ∈ S n converging to p such that w n (q n ) also converges to w(p). Applying Eq. ( 10) with ε equal to half the minimum distance between two points in S ensures that T n converges to T on a set F with full Lebesgue measure in Ω. Choose a point x in the intersection of the cell Vor w S (x) and of F , and consider the sequence q n = T n (p). This sequence converges to p, and by definition one has :

φ n (p n ) -φ(p) = 1 2 ( x -p 2 -x -p n 2 + w(p) -w(p n ))
Using the uniform convergence of φ n to φ, one deduces that w(q n ) converges to w(p).

We now prove by contradiction that if (p n ) converges to p, then lim sup w n (p n ) is at most w(p). Suppose not: taking subsequence if necessary, the limit of w n (p n ) is larger than w(p) by a positive η. For every point x in Ω, we use the triangle inequality to get

x -p n 2 -w n (p n ) ≤ x -q n 2 -w n (q n ) + r n (11) with r n := q n -p n 2 + 2D p n -q n + w n (q n ) -w n (p n )
and D is the maximum distance between a point in Ω and a point in the ball B(0, L) defined by the assumption (iii). Using the convergence of (p n ) and (q n ) to the same point p, and the assumption on the limits of w n (p n ) and w n (q n ), we obtain lim n→+∞ r n ≤ -η. Combining this with Eq. ( 11) shows that the cell of p n in the power diagram Vor wn Sn does not intersect Ω for n large enough. This contradicts the hypothesis that p n is a Dirac with positive mass in the support of ν n .

The proof that lim inf w n (p n ) is larger than w(p) is very similar but a bit longer, as it requires the use of the zero-mean assumption for the convex potentials. These two bounds for lim inf / sup w n (p n ) conclude the proof.

Implementation

In the first paragraph, we give some details our implementation of the convex optimization method proposed in [START_REF] Aurenhammer | Minkowski-type theorems and least-squares clustering[END_REF] for a fixed target measure. Then, we explain how we compute the hierarchical decomposition of the target measure needed for the multiscale algorithm. 4.1. For a fixed target measure. Solving optimal transport between a probability measure µ with density ρ and a discrete measure ν = p∈S λ p δ p amounts to finding the minimum of the convex function Φ given in Theorem 2.(iii):

Φ(w) = p∈S λ p w(p) - Vor w S (p) ( x -p 2 -w(p))ρ(x)dx with ∂Φ ∂w(p) (w) = λ p - Vor w S (p)

ρ(x)dx

We need three ingredients to achieve this goal: an efficient and robust implementation of power diagram computation, robust numerical integration functions, and convex optimization software. In this paragraph, we discuss and motivate our choices regarding these three aspects.

Power diagrams. We use the Regular_triangulation_2 package from CGAL [cga]. It is tempting to try to avoid recomputing the whole power diagram for every evaluation of the function Φ by using the same approach that was used in [START_REF] Machado Manhães De Castro | Filtering relocations on a Delaunay triangulation[END_REF] to maintain the Delaunay triangulation. However, as shown in Figure 3(a), the topology of the power diagram keeps changing until the very last steps of the optimization, thus discarding this approach.

Numerical integration. In our C++ implementation, a measure µ with density ρ is represented by an object which can answer the following two queries. Given a convex polygon P = [a 0 , . . . , a N = a 0 ], and a function f from P to R, the class should provide a way to compute:

(1) the mass of P , i.e. P ρ(x)dx ;

(2) the integral of f over P , i.e. P f (x)ρ(x)dx.

In practice, we only use it the second query for the functions f : x → x-x 0 2 . We developped two different models of measure with density.

The first one is the uniform measure on a convex polygon R. In this case, computing the mass of a polygon P amounts to computing the area of the intersection P ∩R of two convex polygons. The integral of the squared distance function x → x -x 0 2 over the polygon P ∩ R is computed by triangulating P and summing the integral over each triangle T . The integral on T can be obtained in closed-form: if one denotes by cov(T, x 0 ) the covariance matrix of T with base point x 0 , then

T x -x 0 2 dx = cov(T, x 0 ) 11 + cov(T, x 0 ) 22
The second model corresponds to the density obtained from a grayscale image. We assume that the density ρ is constant on each square pixel p i,j = [i, i + 1) × [j, j + 1), equal to the value a i,j . We then consider:

P ρ(x)dx = i,j a i,j area(P ∩ p i,j ) (12)
P f (x)ρ(x)dx ≃ i,j a i,j area(P ∩ p i,j )f (i, j) (13)
Note that it is not possible to simply replace the area of P ∩ p i,j by zero or one depending on whether P intersects p i,j or not, thus disallowing a more efficient GPU implementation. However, since the area of P ∩ p i,j needs to be computed only for pixels containing edges or vertices of P the algorithm we use remains rather efficient. Pixels on edges are dealt with while applying Bresenham's algorithm to raster the polygon. The coverage of pixels containing vertices of P is obtained simply by computing the intersection of the polygon P with the square p i,j .

Convex optimization. We tried several approaches for the actual convex optimization. All of these methods use the following rough scheme to construct the sequence of weight vectors (w k ):

(i) Determine a descent direction d k ;

(ii) Determine a timestep s k , and set w k+1 = w k + s k d k . Methods to choose the descent direction d k include gradient methods, where d k is simply -∇Φ(w k ), Newton methods for which

d k = -[D 2 Φ(w k )] -1 ∇Φ(w k ))
and quasi-Newton methods. In quasi-Newton methods D 2 Φ(w k ) is not computed exactly, but estimated from previous evaluations of the gradients. We chose the widely used low-storage version of the BFGS scheme [START_REF] Fletcher | Practical methods of optimization[END_REF], implemented in C in libLBFGS.

The timestep s k is determined by a search along the line starting from w k with direction d k . Here again, the literature is very vast, as there is a tradeoff between finding a good step size (the best choice would be to minimize the function s → Φ(w k + sd k )) and requiring as few functions evaluations as possible -recall that in our case a function evaluation requires the construction of a complete Power diagram! Figure 2.(a) shows that gradient descent methods are outperformed by quasi-Newton ones, regardless of the choice of line search. It also shows that the choice of line search method is not as important -barring the fixed-step scheme. For all remaining experiments, we use the low-storage BFGS method with Moré-Thuente line search [START_REF] Moré | Line search algorithms with guaranteed sufficient decrease[END_REF]. 4.2. Decomposition of the target measure. Suppose for now that the measure ν is discrete; we will explain in the next paragraph how to convert an image to such a measure. From this measure, we construct a sequence of discrete probability measures (ν ℓ ), with ν ℓ = p∈S ℓ λ p,ℓ δ p such that ν 0 = ν, and that the number of points of the support of ν ℓ decreases as ℓ increases. The parameters of our algorithm are the number L of levels in the decomposition, and for each level ℓ, the number of points n(ℓ) in the support of the measure ν ℓ . In practice, we found that choosing n(ℓ) = n(0)/k ℓ with k = 5 usually provides good results.

Lloyd's algorithm. Theorem 3 suggests that if we want to be able to construct a good initial weight vector for the problem (µ, ν ℓ ) from a weight vector adapted to (µ, ν ℓ+1 ) we need to have ν ℓ+1 as close as possible to ν ℓ in the Wasserstein sense. Given the constraints that ν ℓ+1 is supported on n(ℓ + 1) points, this means

ν ℓ+1 ∈ arg min{Wass 2 (ν, ν ℓ ); | spt(ν)| ≤ n(ℓ + 1)}.
This minimization problem is equivalent to a weighted k-means problem, with k = n(ℓ + 1). Since it is hopeless to solve this problem exactly, we use the standard Lloyd's iterative algorithm to find a good local minimum.

We initialize the algorithm using a random sample S 0 ℓ+1 of n(ℓ + 1) points drawn independently from ν ℓ . We then apply Lloyd's descent step to S n ℓ+1 to obtain S n+1 ℓ+1 , stopping when the points do not move more than a given threshold between two successive steps. This procedure provides us with the support S ℓ+1 of our measure. We define π ℓ to be the application which maps a point p in S ℓ to its nearest neighbor in S ℓ+1 . The values of (λ p,ℓ+1 ) p∈S ℓ+1 are defined by Eq. (5).

Initial quantization of the target measure. Often, the target measure is not a discrete measure but a measure ν I with density σ : Ω → R (such as a grayscale image). In this case we apply Lloyd's algorithm to the measure ρ ′ in order to obtain an initial quantization ν = p∈S λ p δ p of the original measure ν I with a prescribed number of points N .

Results

We will use the following datasets in our experiments. We denote by λU the uniform probability measure on the square λS = [0, λ512] × [0, λ512]. For λ = 1 we will simply write U and S. By L, we denote is the standard grayscale picture of Lena on the square S. Given a measure with density D, we will denote by D N a quantization of this measure with N points, obtained using Lloyd's algorithm. The decomposition of measures we work with are all obtained with the same parameters: 5 levels in the decomposition (including the original one), and level ℓ being made of N/5 ℓ Dirac masses. we show the evolution of the . ∞ distance between the weight vector obtained at a given time, and the optimal one w sol . This optimal weight vector had been previously obtained by running the algorithm with a target accuracy of ∇Φ(w) ∞ < 10 -9 . The advantage of our multiscale method over the original convex optimization is especially important when the source and target measure are far from each other. Table 1 compares the running time of the original and multiresolution algorithms to compute a weight vector adapted to the problem of optimally transporting λU to U 1000 with a given accuracy ∇Φ(w) ∞ < ε. The speedup increases as λ goes to zero, i.e. as the measure λU becomes more concentrated around the lower-left corner of the original square S.

Computation of Wasserstein distances.

We use the approach described in Section 3.3 to obtain lower and upper bounds on the Wasserstein distance between µ and ν at every step of the algorithm. 1. Running time of the original and multiscale algorithm to find a weight vector such that ∇Φ(w) ∞ < 10 -6 . of the two corresponding measures, which can be computed from the L 2 optimal transport. uses optimal transport maps as a replacement for the linear interpolation µ t = (1 -t)µ + tν. Displacement interpolation can be a useful tool for the interpolation of grayscale image, when the gray value of a pixel can be interpreted as a density of some quantity (e.g. satellite views of clouds, preprocessed so that the gray level ranges from black to white depending on the thickness of the cloud). We make use of the transport map computed using the multiscale algorithm. Recall that in order to apply this algorithm to a target measure with density σ : Ω → R, we had to compute a first quantization of σ, ν = p∈S λ p δ p using Lloyd's algorithm. By construction of ν, and by definition of the optimal weight vector ω, one has for every point p in S This suggests a way to construct an interpolation between σ and ρ. Given a time t, consider the weight vector w t = tw, and the corresponding Power diagram (Vor S,wt ). Now, we define the interpolant ρ t at time t as the only piecewise-constant function ρ t on Ω obtained by spreading the mass of λ p on the intersection of the cell Vor S,wt (p) with Ω, i.e. for every point x in Vor S,wt (p), define ρ t (x) := λ p / area(Vor S,wt (p)). An example of this interpolation is presented in Figure 4.

Discussion

In this paper we have presented a simple way to increase the efficiency of the convex optimization algorithm introduced in [AHA98] to solve the optimal transport problem. We also discussed how our multiscale approach can be used to obtain fast estimation of Wasserstein distances between images.

This first step suggests that, in order to obtain faster computations of optimal transport, one has to better understand the geometry of the function Φ. For instance, it is currently not possible to obtain complexity estimates for this approach because: (i) nothing is known about the shape and size of the basin around the minimizer where Newton's method has quadratic convergence and (ii) the stability result (Theorem 3) is not quantitative. Understanding these two problems could open the way to even more efficient computations of optimal transport maps.

We also believe that this multiscale approach can be useful in the solution of more geometric problems with a similar structure. An example of such a problem is Minkowski's problem: given a set of normals n 1 , . . . , n N and a set of areas λ 1 , . . . , λ N such that i λ i n i vanishes, find a convex polytope whose facets normals are among the ( n i ), and such that the facet corresponding to n i has an area of exactly λ i . This problem has a similar multiscale structure as optimal transport, and can be also solved by minimizing a convex functional [START_REF] Lachand-Robert | Minimizing within convex bodies using a convex hull method[END_REF], and would probably benefit from a multiscale approach. A second example is the problem of designing a reflector antenna with prescribed image measure at infinity, which can also be formally cast as an optimal transport problem (Section 4.2.5 in [START_REF] Oliker | Mathematical aspects of design of beam shaping surfaces in geometrical optics[END_REF]).

Figure 1 .

 1 Figure 1. Interpolation between the standard pictures Photograph and Peppers obtained with our algorithm (see Section 5.3). The target image was quantized with 625 Dirac masses for the first row, and 15625 Dirac masses for the second row. The cells in these pictures are interpolation between power cells (on the left) and Voronoi cells (on the right), obtained by linearly varying the weights from their value on the left picture to zero. Their color vary so that the product of the area of a cell multiplied by its gray level remains constant over time.

Proof.

  By definition, both this measure ν and the target measure ν are supported on the same set S. Moreover, m := ∇Φ(w) 1 = p∈S λ p -Vor w S (p) ρ(x)dx corresponds to the amount of mass that has been mistransported. The cost of transporting this mass back at the right place in S is at most √ mD 2 , where D = diam(S).

Figure 2 .

 2 Figure 2. Speed of convergence, measured by the L ∞ distance between the weight vector at a given time/step and the optimal one. (a) Comparison of simple convex optimization algorithms: gradient descent (red) with fixed step (solid) or strong Wolfe line search (dashed), and low-storage BFGS algorithm (blue) with strong Wolfe (solid) or Moré-Thuente line-search (dashed). (b) and (c) Comparison between the original algorithm of [AHA98] (red) and the multiscale one (blue).

Figure 3 .

 3 Figure 3. (a) Percentage of points in the support of the target measure whose Power cell intersects the support of the source measure during the execution of various convex optimization algorithms (colors are the same as in Fig. 2(a)). (b) and (c) Estimation of Wasserstein distance: in red (resp. blue), the lower and upper bounds obtained by the multiscale (resp. original) algorithm as a function of time, and in green the correct value.

  Figure 3(b) and 3(c) compare the evolution of these two bounds as a function of the runtime of the original and the multiscale algorithm. 5.3. Displacement interpolation of images. The concept of displacement interpolation of two probability measures was introduced in [McC97]. It

Figure 4 .

 4 Figure 4. First and second rows: An interpolation between a picture of G. Monge and photograph of B. Riemann (with N =625 and 15k respectively). The intermediary steps are obtained using McCann's displacement interpolation [McC97] of the two corresponding measures, which can be computed from the L 2 optimal transport.

  Vor S (p)∩Ω σ(x)dx = λ p = Vor S,w (p)∩Ω ρ(x)dx.

Acknowledgements. ANR grant GIGA ANR-09-BLAN-0331-01 and Université Grenoble I MSTIC grant GEONOR.