Adrien Blanchet
email: adrien.blanchet@univ-tlse1.fr

Guillaume Carlier
email: carlier@ceremade.dauphine.fr

OPTIMAL TRANSPORT AND COURNOT-NASH EQUILIBRIA

Keywords: pdfcreator=HAL, pdfproducer=PDFLaTeX, pdfsubject=Mathematics [math]/Analysis of PDEs [math.AP], Humanities and Social Sciences/Economics Cournot-Nash equilibria, mean-field games, optimal transport, externalities, Monge-Ampère equations, convexity along generalised geodesics

published or not. The documents may come

,

INTRODUCTION

Since Aumann's seminal works [START_REF] Aumann | Existence of competitive equilibria in markets with a continuum of traders[END_REF][START_REF] Aumann | Markets with a continuum of traders[END_REF], models with a continuum of agents have occupied a distinguished position in economics and game theory. [START_REF] Schmeidler | Equilibrium points of nonatomic games[END_REF] introduced a notion of non-cooperative equilibrium in games with a continuum of agents and established several existence results. In Schmeidler's own words: Non-atomic games enable us to analyze a conflict situation where the single player has no influence on the situation but the aggregative behavior of "large" sets of players can change the payoffs. The examples are numerous: Elections, many small buyers from a few competing firms, drivers that can choose among several roads, and so on.

Following the approach of [START_REF] Kohlberg | On equilibrium allocations as distributions on the commodity space[END_REF] to Walras equilibrium analysis, Mas-Colell [1984] reformulated Schmeidler's analysis in terms of joint distributions over agents' actions and characteristics and, in particular, the concept of Cournot-Nash equilibrium distributions. Not only Mas-Colell's reformulation enabled him to obtain general existence results in an easy and elegant way but it is flexible enough to accommodate quite weak assumptions on the data (which is relevant in the framework of games with incomplete information and a continuum of players for instance). Roughly speaking, in analysing Cournot-Nash equilibria in the sense of [START_REF] Mas-Colell | On a theorem of Schmeidler[END_REF] one can take great advantage of (topological but also geometric) properties of spaces of probability measures. With this respect, it is natural to expect that optimal transport theory (which is an extremely active field in research in mathematics both from an applied and fundamental point, as illustrated by the monumental textbook [START_REF] Villani | Optimal transport: old and new[END_REF]) may be useful.

Even though there are very general existence results for Cournot-Nash equilibria (see for instance [START_REF] Kahn | On Cournot-Nash equilibrium distributions for games with a nonmetrizable action space and upper semi-continuous payoffs[END_REF]) in the literature, we are not aware of classes of problems where there is uniqueness and a full characterisation of such equilibria which is tractable enough to obtain close-form solutions or efficient numerical computation schemes. One of our goals is precisely to go one step beyond abstract existence results (in mixed or pure strategies) and to identify classes of non-atomic games where Cournot-Nash equilibria are unique and can be fully characterised or numerically computed.

Given a space of players types X endowed with a probability measure µ ∈ P(X) (which gives the exogenous distribution of the type of the agents), an action space Y and a cost Φ: X × Y × P(Y ) → R, x-type agents taking action y pay the cost Φ(x, y, ν) where ν ∈ P(Y ) represents the action distribution. The fact that this cost depends on the other agents actions only through the distribution ν means that who plays what does not matter i.e. the game is anonymous. A Cournot-Nash equilibrium is a joint probability measure γ ∈ P(X × Y ) with first marginal µ such that

(1.1) γ({(x, y) ∈ X × Y : Φ(x, y, ν) = min z∈Y Φ(x, z, ν)}) = 1
where ν represents γ's second marginal. The probability γ is naturally interpreted by saying that γ(A × B) is the probability that agents have their type in A and an action in B. The equilibrium is called pure if, in addition, γ is carried by a graph i.e. µ-a.e. the agents play in pure strategy. Condition (1.1) means that agents choose cost minimising strategies given their type and ν so that, finally, imposing that ν is the second marginal of γ is a simple self-consistency requirement.

In the sequel, we will restrict ourselves to the additively separable case where Φ(x, y, ν) = c(x, y) + V[ν](y), which seems to be a necessary limitation for the optimal transport approach we will develop. Under this separability specification, the connection with optimal transport is almost obvious: if γ is a Cournot-Nash equilibrium it necessarily minimises the average of c among probability measures having µ and ν (which is a priori unknown) as marginals i.e. it solves the optimal transport problem:

(1.2) W c (µ, ν) := inf γ∈Π(µ,ν) X×Y c(x, y) dγ(x, y)

where Π(µ, ν) is the set of joint probabilities having µ and ν as marginals. In an euclidean setting, there are well-known conditions on c and µ which guarantee that such an optimal γ necessarily is pure whatever ν is and this of course implies purity of equilibria.

If we go one step further and assume that V[ν] is the differential of some functional E (see Section 3 for a precise definition), it turns out that if ν is a minimiser of E[ν]+ W c (µ, ν) and γ solves (1.2) then it is a Cournot-Nash equilibrium. This gives a variational device to find equilibria: first find ν by minimising E[ν] + W c (µ, ν) and then find γ by solving the optimal transport problem (1.2) between µ and ν. This variational approach actually gives new existence results. To the best of our knowledge, usual general existence proofs are via fixed-point arguments and thus require a lot of regularity for the dependence of V[ν] with respect to ν, in a compact metric setting, it is typically asked that ν → V[ν] is continuous (or at least upper-semi continuous in some sense) from the set of probabilities equipped with the weak- * topology to the set of continuous functions equipped with the supremum norm. This is harmless if Y is finite but extremely restrictive in general, in particular it excludes the case of a purely local dependence which is relevant to capture congestion effects (actions that are frequently played are more costly). In contrast, the variational approach will enable us to treat such local congestion effects. If E (the primitive of V in some sense) is convex then equilibria and minimisers coincide and strict convexity gives uniqueness of the second marginal, ν, of the equilibrium. Such a convexity is quite demanding in applications but we shall prove that in an euclidean setting and for a quadratic c (and more generally strictly convex c's in dimension one), there is some hidden convexity (in the spirit of the seminal results of [START_REF] Mccann | A convexity principle for interacting cases[END_REF]) in the problem from which one can deduce uniqueness of equilibria but also a characterisation in terms of a nonlinear partial differential equation of Monge-Ampère type. This partial differential equation cannot be solved explicitly in general but, in dimension one, it is easy to solve the variational problem numerically in an efficient way as we shall illustrate on several examples. Another advantage of the variational approach is that it allows for an elementary (in-)efficiency analysis of the equilibrium and the design of a tax system to restore the efficiency of the equilibrium (see Section 5). Of course, the variational approach described above presents strong similarities with the potential games of [START_REF] Monderer | Potential games[END_REF] and our framework is very close to that of [START_REF] Konishi | Pure strategy nash equilibrium in a group formation game with positive externalities[END_REF] or [START_REF] Lebreton | Games of social interactions with local and global externalities[END_REF] in the case of a finite number of players; however we are not aware of any extension of the analysis of [START_REF] Monderer | Potential games[END_REF] to the case of a continuum of players.

Apart from our results on Cournot-Nash equilibria, another objective of the paper is to contribute to popularise the use of optimal transport in economics. Several recent papers have fruitfully used optimal transport arguments in such different fields as hedonics and matching problems [START_REF] Chiappori | Hedonic price equilibria, stable matching, and optimal transport: equivalence, topology, and uniqueness[END_REF], [START_REF] Ekeland | Existence, uniqueness and efficiency of equilibrium in hedonic markets with multidimensional types[END_REF]), multidimensional screening [START_REF] Carlier | Duality and existence for a class of mass transportation problems and economic applications[END_REF], [START_REF] Figalli | When is multidimensional screening a convex program?[END_REF]) or urban economics [START_REF] Blanchet | City equilibria via optimal transport and differential equations[END_REF], [START_REF] Carlier | Equilibrium structure of a bidimensional asymmetric city[END_REF]). We believe that cross-fertilisation between economics and optimal transport will rapidly develop. This is why we have included in Appendix A some basic results from optimal transport theory which we hope can serve as a comprehensive introduction to this vast subject to an economists readership.

The present introduction would neither be complete nor fair without an explicit reference to the mean-field games theory of Lasry andLions [2006a,b, 2007]. Indeed, our variational approach is largely inspired by the Lasry and Lions optimal control approach to mean-field games (that has some similarities with optimal transport), but also mean-field games theory enables to treat considerably richer situations than the somehow static one we treat here. Another line of research we would like to mention concerns congestion games (another example of potential games) and the literature on the cost of anarchy (see [START_REF] Roughgarden | Selfish routing and the price of anarchy[END_REF] and the references therein), indeed the variational approach we develop presents some similarities with the variational approach to Wardrop equilibria on congested networks and in both cases equilibria are socially inefficient.

The paper is organised as follows. In Section 2, we introduce the model, define equilibria and emphasise some connections with optimal transport. In Section 3, we adopt a variational approach and prove that for a large class of interactions, equilibria naturally arise as local minimisers of a certain functional. Section 4 is devoted to further uniqueness and variational characterisation of equilibria results thanks to notions of displacement convexity arising in optimal transport, we also characterise the equilibrium via a certain nonlinear partial differential equation and compute numerically the equilibrium in dimension one. Section 5 concludes. The proofs as well as well as a presentation of various results from optimal transport theory which are used throughout the paper are gathered in the Appendix.

THE EQUILIBRIUM MODEL

The model consists of a compact metric type space X equipped with a Borel probability measure µ ∈ P(X), giving the distribution of types, a compact metric action space Y , a reference 1 Borel non-negative measure m 0 , a continuous function c ∈ C(X × Y ) and interactions are captured by a map which to every action distribution ν ∈ P(Y ) ∩ L 1 (m 0 ) associates a function V[ν] defined m 0 -almosteverywhere. Given an action distribution ν, x-type agents taking action y then incur the additively separable cost c(x, y) + V[ν](y). The unknown is a probability distribution γ ∈ P(X ×Y ), with the interpretation that γ(A×B) is the probability that an agent has her type in A and takes an action in Y , such a γ induces as action distribution ν, its second marginal which we denote ν = π Y # γ. By construction, the first marginal of γ, π X # γ should be equal to µ. Since we will be interested by efficiency (or rather inefficiency) properties of equilibria, we will also impose that γ has finite social cost, where the latter is given by

SC = X×Y (c(x, y) + V[ν](y)) dγ(x, y) = X×Y c(x, y) dγ(x, y) + Y V[ν](y) dν(y)
since c is continuous the first term is finite for every γ, but the second requires the action marginal ν to belong to the domain

(2.1) D := {ν ∈ L 1 (m 0 ) : V[ν] ∈ L 1 (ν)} = {ν ∈ L 1 (m 0 ) : Y |V[ν]| dν < +∞}.
Cournot-Nash equilibria are then defined as follows.

Definition 2.1 γ ∈ P(X × Y ) is a Cournot-Nash equilibria if its first marginal is µ, its second marginal, ν, belongs to D and there exists ϕ ∈ C(X) such that

(2.2) c(x, y) + V[ν](y) ≥ ϕ(x)
for all x ∈ X and m 0 -a.e. y with equality γ-a.e..

A Cournot-Nash equilibrium γ is called pure whenever it is carried by a graph i.e. is of the form γ = (id, T ) # µ for some Borel map T : X → Y .

The previous definition is slightly different from that of [START_REF] Mas-Colell | On a theorem of Schmeidler[END_REF] because we require the action distribution to be absolutely continuous with respect to m 0 , so as to take into account congestion effects as explained in the examples below. This makes the existence of equilibria nontrivial, indeed, when ν → V[ν] is continuous from (P(Y ), w-*) to (C(Y ), . ∞ ) (as is the case for instance when V[ν](y) = Y φ(y, z) dν(z) with φ continuous) standard fixed-point arguments immediately give the existence of Cournot-Nash equilibria but here, we do not have such regularity.

Examples Holiday choice

Let us consider a population of agents whose location is distributed according to some probability distribution µ ∈ P(X) where X is some compact subset of R 2 (say). These agents have to choose their holidays destination (possibly in mixed strategy). The set of possible holiday destinations is some compact subset of the plane Y (it can be X, a finite set, ...). The commuting cost from x to y is c(x, y). In addition to the commuting cost, agents incur costs resulting from interactions with other agents, this is captured by a map ν → V[ν] that can be modelled as follows.

A natural effect that has to be taken into account is congestion, i.e. the fact that more crowded location results in more disutility for the agents. Congestion thus requires to consider local effects and actually imposes that ν is not too concentrated; a way to capture this is to impose that ν is absolutely continuous with respect to some reference probability measure m 0 . Still denoting by ν the Radon-Nikodym derivative of ν, a natural congestion cost is of the form y → f (ν(y)) with f non-decreasing. In addition to the negative externality due to congestion effect, there may be a positive externality effect due to the positive social interactions between agents which can be captured through a non-local term of the form y → Y φ(y, z) dν(z) where for instance φ(y, .) is minimal for z = y so that the previous term represents a cost for being far from the rest of the population. Finally, the presence of purely geographical factors (e.g. distance to the sea) can be reflected by a term of the form y → v(y). The total externality cost generated by the distribution ν combines the three effects of congestion, positive interactions and geographical factors and can then be taken of the form

V[ν](y) = f (ν(y)) + Y φ(y, z) dν(z) + v(y).

Technological choice

Consider now a simple model of technological choice in the presence of externalities. There is a set of consumers indexed by a type x ∈ X drawn according to the probability µ, and a set of technologies Y for a certain good (cell-phone, computer, tablet...). On the supply side, assume there is a single profit maximising profit firm with convex production cost F (y, .) producing technology y, the supply (equals demand at equilibrium) of this firm is thus determined by the marginal pricing rule p(y) = ∂ ν F (y, ν(y)). Agents aim to minimise with respect to y a total cost which is the sum of their individual purchasing cost c(x, y)+ p(y) = c(x, y)+ ∂ ν F (y, ν(y)) and an additional usage/maintenance or accessibility cost which is positively affected by the number of consumers having purchased similar technologies i.e. a term of the form Y φ(y, z) dν(z) where φ is increasing in the distance between technologies y and z.

Connection with optimal transport and purity of equilibria

For ν ∈ P(Y ), let Π(µ, ν) denote the set of probability measures on X ×Y having µ and ν as marginals and let W c (µ, ν) be the least cost of transporting µ to ν for the cost c i.e. the value of the Monge-Kantorovich optimal transport problem:

W c (µ, ν) := inf γ∈Π(µ,ν) X×Y c(x, y) dγ(x, y)
let us also denote by Π o (µ, ν) the (nonempty) set of optimal transport plans i.e.

Π o (µ, ν) := {γ ∈ Π(µ, ν) : X×Y c(x, y) dγ(x, y) = W c (µ, ν)}.
A first link between Cournot-Nash equilibria and optimal transport is based on the following straightforward observation.

Lemma 2.2 If γ is a Cournot-Nash equilibrium and ν denotes its second marginal then γ ∈ Π o (µ, ν).

Proof: Indeed, let ϕ ∈ C(X) be such that (2.2) holds and let η ∈ Π(µ, ν) then we have

X×Y c(x, y) dη(x, y) ≥ X×Y (ϕ(x) -V[ν](y)) dη(x, y) = X ϕ(x) dµ(x) - Y V[ν](y) dν(y) = X×Y c(x, y) dγ(x, y) so that γ ∈ Π o (µ, ν). Q.E.D.
The previous proof also shows that ϕ solves the dual of W c (µ, ν) (see Appendix (A.3)) i.e. maximises the functional In an euclidean setting, there are well-known conditions on c and µ which guarantee that such an optimal γ necessarily is pure whatever ν is. It is the case for instance if µ is absolutely continuous with respect to the Lebesgue measure, c(x, y) is a smooth and strictly convex function of x -y (see [START_REF] Mccann | The geometry of optimal transportation[END_REF] who extended the seminal results of [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF] in the quadratic cost case), or more generally, when it satisfies a generalised Spence-Mirrlees condition (see [START_REF] Carlier | Duality and existence for a class of mass transportation problems and economic applications[END_REF] for details):

Corollary 2.3 Assume that X = Ω
where Ω is some open connected bounded subset of R d with negligible boundary, that µ is absolutely continuous with respect to the Lebesgue measure, that c is differentiable with respect to its first argument, that ∇ x c is continuous on R d × Y and that it satisfies the generalised Spence-Mirrlees condition:

for every x ∈ X, the map y ∈ Y → ∇ x c(x, y) is injective, then for every ν ∈ P(Y ), Π o (µ, ν) consists of a single element and the latter is of the form γ = (id, T ) # µ hence every Cournot-Nash equilibrium is pure.

A VARIATIONAL APPROACH

In this section, we will see that in many relevant cases, one may obtain equilibria by the minimisation of some functional over a set of probability measures2 . The main assumption for this variational approach to be valid is that the interaction map V[ν] has the structure of a differential i.e. that V[ν] can be seen as the first variation of some function ν → E [ν]. In this case, the variational approach is based on the observation that the equilibrium condition is the first-order optimality condition for the minimisation of W c (µ, ν) + E[ν].

Interaction maps which are differentials

The main assumption for the variational approach to be valid is that ν → V[ν] is a differential in the following sense:

Definition 3.1 (Differential) Let D be defined by (2.1). The map ν ∈ D → V[ν] is a differential on D if D is convex and there exists E: D → R such that for every (ρ, ν) ∈ D 2 , V[ν] ∈ L 1 (ρ) and lim ε→0 + E[(1 -ε)ν + ερ] -E[ν] ε = Y V[ν] d(ρ -ν) i.e. V[ν] is the first variation of E which we denote V[ν] = δE δν .
Before going any further, let us consider some examples to illustrate the previous definition.

Local term

Let us consider first the case of a local dependence, again m 0 is our reference measure and for ν ∈ D := P(Y ) ∩ L 1 (m 0 ) and m 0 -a.e. y:

V[ν](y) = f (y, ν(y))
for some continuous f . Assume first that f is bounded and define, for all ν ∈ D:

(3.1) F (y, ν) = ν 0 f (y, s) ds, E[ν] = Y F (y, ν(y)) dm 0 (y) .
Then since F is Lipschitz in ν uniformly in y, it easily follows from Lebesgue's dominated convergence theorem that V[ν] is the differential of E on D. Now, rather assume that f satisfies the growth condition

(3.2) a(ν α -1) ≤ f (y, ν) ≤ b(ν α + 1)
for some a ≥ 0, b > 0, α > 0, m 0 -a.e. y and every ν ≥ 0. For p = α + 1, the corresponding energy functional E is then defined for all ν ∈ D := P(Y ) ∩ L p (m 0 ) as above by (3.1). Thanks to (3.2) and Lebesgue's dominated convergence theorem, V is the differential of E on D, V[ν] ∈ L p ′ (m 0 ) (with p ′ the conjugate exponent of p i.e. p ′ = p/(p -1)) as soon as ν ∈ D. Apart from the technical growth condition (which is useful to apply Lebesgue's theorem and guarantee that V[ν]ν is integrable) we therefore see that local V's are differentials. In Section 3.3, we will treat local V's under a different Inada-like condition on f which is more customary in economics and will ensure that ν remains positive hence simplifying the equilibrium/optimality condition.

Non-local interaction term and the role of symmetry

Let us now consider the case of (pairwise) interactions where V[ν] is defined by

V[ν](y) = Y φ(y, z) dν(z)
for some φ ∈ C(Y × Y ). It is then natural to define the quadratic functional

E[ν] = 1 2 Y ×Y φ(y, z) dν(y) dν(z) .
By expanding in ε,

E[ν + ε(ρ -ν)], its differential is immediate to compute lim ε→0 E[ν + ε(ρ -ν)] -E[ν] ε = 1 2 φ(y, z)[ dν(y) d(ρ -ν)(z) + dν(z) d(ρ -ν)(y)] = 1 2 [φ(y, z) + φ(z, y)] dν(z) d(ρ -ν)(y) . So that δE δν (y) = Y φ sym (y, z) dν(z) : φ sym (y, z) = φ(y, z) + φ(z, y) 2 .
Hence V is the differential of E on P(Y ) as soon as φ is symmetric 3 i.e. φ(y, z) = φ(z, y) (which is the case for instance if φ is the function of the distance between y and z). Note that the assumption that V is a differential requires φ to be symmetric. 4 Of course, one can combine the previous examples and consider a V which is the sum of a symmetric interaction term and a local term, such V's still have the structure of a differential.

Minimisers are equilibria

Throughout this paragraph, we assume that

(3.4) V[ν] = δE δν on D.
We then consider the variational problem

(3.5) inf ν∈D J µ [ν] where J µ [ν] := W c (µ, ν) + E[ν].
To prove that minimisers of (3.5) are equilibria, we first need to be able to differentiate the term W c (µ, ν) with respect to ν, this is possible thanks to Lemma A.1 proved in Appendix A but it requires more structure on X, µ and c: in particular X is a connected subset of R d , c is differentiable with respect to x and µ is equivalent to the Lebesgue measure. See Appendix B.1 for the proof. Let us mention however that the optimality condition for (3.5) is the following: there is a constant M such that

(3.6) ϕ c + V[ν] ≥ M ϕ c + V[ν] = M ν-a.e. ,
4 In a similar way, if we consider the case of higher-order interactions

V[ν](y) = Y m φ(y, .) dν ⊗m = Y m φ(y, z1, . . . , zm) dν(z1) . . . dν(zm)
where φ ∈ C(Y m+1 ) satisfies the symmetry relations

(3.3) φ(y, z1, . . . , zm) = φ(z1, y, . . . , zm) = • • • = φ(zm, z1, . . . , y),
for all (y, z1, . . . , zm) ∈ Y m+1 , then V is the differential of

E [ν] = 1 m + 1 Y m+1 φ dν ⊗(m+1) .
where ϕ c is the c-transform of ϕ as in (2.3).

To deduce an existence result from Theorem 3.5, assume that V[ν] is defined for

ν ∈ P(Y ) ∩ L 1 (m 0 ) by (3.7) V[ν](y) := f (y, ν(y)) + Y φ(y, z) dν(z)
where φ ∈ C(Y × Y ) is symmetric, f is continuous and non-decreasing with respect to its second argument and satisfies the growth condition (3.2) for some a > 0, b > 0 and α > 0. For p = α+1, the corresponding energy functional is then defined for all ν ∈ D := P(Y ) ∩ L p (m 0 ) by

E[ν] = Y F (y, ν(y)) dm 0 (y) + 1 2 Y 2 φ(y, z) dν(y) dν(z)
where F is defined by (3.1). The functional F (y, .) is convex and satisfies the growth condition

a(p -1 ν p -ν) ≤ F (y, ν) ≤ b(p -1 ν p + ν) .
Hence V[ν] ∈ L p ′ (m 0 ) as soon as ν ∈ D and thus, by Hölder's inequality, V[ν]ρ ∈ L 1 (m 0 ) for every ρ and ν in D.

Corollary 3.3 (Existence of equilibria by minimisation) Assume that the assumptions of Lemma A.1 hold, that ν → V[ν] is of the form (3.7) where f and φ satisfy the assumptions above, then (3.5) admits minimisers in P(Y ) ∩ L p (m 0 ) so that there exists Cournot-Nash equilibria.

The proof is given in Appendix B.2. Note that this in particular provides existence of equilibria results for the holiday and technological choice model examples above.

Remark 3.4 Under the assumptions of the previous corollary, one can prove that the minimisers are actually bounded: indeed let ν be such a minimiser either ν(y) = 0 or ν(y) > 0 and for m 0 -a.e. such points by the optimality condition (3.6) and (3.2) one should have for some constant M a(ν(y

) α -1) ≤ f (y, ν(y)) = M -ϕ c (y) - Y φ(y, z) dν(z) .
Since ϕ is a c-transform, it is continuous hence bounded on Y and the integral term is bounded since φ is. We therefore have ν ∈ L ∞ (m 0 ).

Let us now emphasise the role of convexity in the variational approach. As expected if E is convex, then finding equilibria and minimising J µ are equivalent: Proposition 3.5 (Equivalence in the convex case) Assume that the assumptions of Theorem 3.2 are satisfied. If E is convex on D then the following statements are equivalent:

• ν solves (3.5) and γ ∈ Π o (µ, ν),

• γ is an equilibrium and ν = Π Y # γ.

If moreover E is strictly convex the following uniqueness result holds:

Corollary 3.6 (Uniqueness in the strictly convex case) Assume that the assumptions of Theorem 3.2 are satisfied. If E is strictly convex then all equilibria share the same second marginal ν. If in addition, the assumptions of Corollary 2.3 are satisfied then there is at most one Cournot-Nash equilibrium.

As an application, let us observe that if the assumptions of Lemma A.1 are satisfied and if V[ν](y) = f (y, ν(y)) with an f which is increasing in ν and satisfies (3.2) then there exists a unique minimiser so the previous uniqueness result holds. This applies naturally to the technological choice equilibrium problem as well as to the holiday choice example with pure congestion or, more generally, in the case where the congestion effects dominate as explained below.

In the case where

E[ν] = Y F (y, ν(y)) dm 0 (y) + 1 2 Y 2 φ(y, z) dν(y) dν(z)
the second non-local term typically favours the concentration of ν (when φ(y, z) is increasing with the distance between y and z for instance) and it is not convex, while the congestion terms fosters dispersion and is convex. There may however be some compensation between the two terms that makes E convex. For instance, by Cauchy-Schwarz inequality, the quadratic form

Y ν 2 (y) dm 0 (y) + Y 2 φ(y, z) dν(y) dν(z)
is positive definite hence convex as soon as

Y 2 φ 2 (y, z) dm 0 (y) dm 0 (z) < 1.
Whence in this case, the uniqueness result of Corollary 3.6 applies.

The case of Inada's condition

We now consider the case where V contains a local congestion term that satisfies an Inada-like condition:

(3.8) lim ν→0 + f (ν) = -∞ and lim ν→+∞ f (ν) = +∞.
This will imply that minimisers of (3.5) are positive m 0 -a.e.. The optimality condition ϕ c + V[ν] = M , for some constant M , will therefore be satisfied m 0 -a.e. which implies the regularity of ν. More precisely, let us consider the case where the interaction are given for ν ∈ P(Y ) ∩ L 1 (m 0 ) by the map:

(3.9)

V[ν](y) = f (ν(y)) + Y φ(y, z) dν(z)
(for the sake of simplicity we have dropped the dependence in y of f ) where

• φ ∈ C(Y × Y ) is symmetric,
• f : (0, +∞) → R is continuous increasing, locally integrable on [0, +∞) and satisfies the Inada condition (3.8). We then define F by F (0) = 0 and F ′ = f so that F is strictly convex, continuous on [0, +∞) and C 1 (0, +∞), bounded from below and coercive i.e. F (ν)/ν → +∞ as ν → +∞. As before, for any ν in P(Y ) ∩ L 1 (m 0 ), we define the associated cost functional

E[ν] =      Y F (ν) dm 0 + 1 2 Y 2 φ(y, z) dν(y) dν(z) if Y F (ν) dm 0 < +∞ +∞ otherwise.
The typical example we have in mind is f (ν) = log(ν) and F (ν) = ν log ν -ν, or simply F (ν) = ν log ν since we are only dealing with probability measures. In this case, the domain of E consists of absolutely continuous measures with finite entropy. Again, we look for equilibria by solving the minimisation problem (3.5). The implication of Inada's condition on the interiority of minimisers is given by the following: Lemma 3.7 (Existence and positivity of minimisers) Under the above assumptions, the variational problem (3.5) admits solutions and if ν is such a solution ν ≥ δ m 0 -a.e. for some δ > 0 and ν ∈ L ∞ (m 0 ).

The proof (see Appendix B.4) relies on the fact that since f (0 + ) = -∞ the functional E abhors a vacuum. Now that we know that minimisers ν exist and are bounded from above and bounded away from 0, under the assumptions of Lemma A.1 it is easy to see, as in the previous paragraph, that they necessarily are equilibria and satisfy the optimality condition (3.6):

f (ν(y)) + Y φ(y, z) dν(z) + ϕ c (y) = M
for some constant M and where as usual ϕ is the Kantorovich potential between µ and ν and ϕ c is its c-transform. Note that this equality is true not only ν-a.e. but m 0 -a.e., one can then invert this relation to deduce that ν coincides m 0 with the continuous function

(3.10) ν(y) = f -1 M -ϕ c (y) - Y φ(y, z) dν(z) .
In particular there exists equilibria that have a continuous representative. 5 Relation (3.10) is however not very tractable in general since it involves the very indirect quantity ϕ c and an integral term. We will see in Section 4 how it can be simplified and reformulated as a nonlinear partial differential equation in the case of a quadratic cost.

For the moment, the Inada condition has just enabled us to prove some further regularity properties of minimisers hence of some special equilibria. Let us summarise all this by: Theorem 3.8 (Main results under the Inada condition) Let V be of the form (3.9) where f and φ satisfy the assumptions of this paragraph. If ν solves (3.5) and γ ∈ Π o (µ, ν) then γ is an equilibrium; in particular, there exists equilibria. Moreover any minimiser ν of (3.5) is bounded and bounded away from 0 and coincides m 0a.e. with the continuous function given by (3.10).

If, in addition, E is convex γ is an equilibrium if and only if γ ∈ Π o (µ, ν) where ν solves (3.5).

If, in addition, E is strictly convex there is a uniqueness of the equilibrium second marginal ν.

HIDDEN CONVEXITY AND FURTHER UNIQUENESS RESULTS

So far, our variational approach has enabled us to prove the existence of equilibria by the minimisation problem (3.5). However, the previous results are not totally satisfying since in general there might exist equilibria that are not minimisers and even if we are only interested in the special equilibria obtained by minimisation, optimality conditions like (3.10) are not tractable enough to provide a full characterisation. Under further convexity conditions that are quite stringent we have seen that equilibria necessarily are minimisers and obtained uniqueness of both. In the case where

E[ν] = Y F (y, ν(y)) dm 0 (y) + 1 2 Y ×Y φ(y, z) dν(y) dν(z)
5 Inada's condition is actually not essential to obtain a relation of the form (3.10). Indeed, in the case of a power congestion function, f (ν) = ν α , α > 0, using the positive part function, one obtains a similar relation

ν(y) = M -ϕ c (y) - Y φ(y, z) dν(z) 1/α + .
there is a competition between the convexity of the congestion term that favours dispersion and the non-convexity of the interaction term so that in general nothing can be said about the convexity of E in the usual sense. We shall see however, that some convexity structure, more adapted to optimal transport, can be used to derive new uniqueness and characterisation results. The aim of this section is precisely to exploit some hidden convexity structure in one dimension and in higher dimensions when the cost is quadratic. This goal can be achieved thanks to the very powerful notion of displacement convexity (or some slight variant of it) due to [START_REF] Mccann | A convexity principle for interacting cases[END_REF]. In recent years, these notions of convexity, intimately linked to optimal transport, have proved to be an extremely useful and flexible tool in particular in the study of nonlinear diffusions, to our knowledge, this is the first time they are used in an economic context, see also [START_REF] Blanchet | City equilibria via optimal transport and differential equations[END_REF]. We refer to Appendix A for a very short presentation and Section 4.1 for a detailed exposition in the easier one-dimensional case. Much more on this rich subject can be found in the books [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF], [START_REF] Villani | Topics in optimal transportation[END_REF][START_REF] Villani | Optimal transport: old and new[END_REF].

Hidden convexity in dimension one

Let us start with the simple one-dimensional case where the intuition is easy to understand: the functional J µ is not convex with respect to ν but it is with respect to T , the optimal transport map from µ to ν. Let us take X = Y = [0, 1], m 0 is the Lebesgue measure on [0, 1], µ is absolutely continuous with respect to the Lebesgue measure, and assume that V[ν] takes the form:

V[ν](y) = f (ν(y)) + v(y) + [0,1] φ(y, z) dν(z)
and that

• the transport cost c is of the form c(x, y) = C(x -y) where C is strictly convex and differentiable, • f is increasing, • v is convex on [0, 1] and φ is convex, symmetric, differentiable and has a locally Lipschitz gradient. As already noted the corresponding cost

E[ν] := 1 0 F (ν(y)) dy + 1 2 [0,1] 2 φ(y, z) dν(y) dν(z) + 1 0 v(y) dν(y)
(with F ′ = f ) is not convex in the usual sense in general and neither is the functional J µ = W c (µ, .) + E. However, we shall see that J µ has good convexity properties when one considers the following interpolation. Let (ρ, ν) ∈ P([0, 1]) 2 then there is a unique optimal transport map T 0 (respectively T 1 ) from µ to ν (respectively from µ to ρ) for the cost c and it is non-decreasing (see [START_REF] Villani | Topics in optimal transportation[END_REF]). For t ∈ [0, 1], let us define:

ν t := T t# µ where T t := ((1 -t)T 0 + tT 1 )
then by construction, the curve t → ν t connects ν 0 = ν to ν 1 = ρ.

Definition 4.1 A functional J : P(Y ) → R ∪ {+∞} is called displacement convex whenever t ∈ [0, 1] → J [ν t ]
is convex (for every choice of endpoints ν and ρ), it is called strictly displacement convex when, in addition

J [ν t ] < (1 -t)J [ν] + tJ [ρ] when t ∈ (0, 1) and ρ = µ.
We claim that J µ is strictly displacement convex; indeed, take (ν, ρ) two probability measures in the domain of E (which is convex by convexity of F ), define ν t as above and, let us consider the four terms in J µ separately:

• By definition of W c , ν t and the strict convexity of C we have

W c (µ, ν t ) ≤ 1 0 C(x -((1 -t)T 0 (x) + tT 1 (x))) dµ(x) ≤ (1 -t) 1 0 C(x -T 0 (x)) dµ(x) + t 1 0 C(x -T 1 (x)) dµ(x) = (1 -t)W c (µ, ν) + tW c (µ, ρ)
with a strict inequality if t ∈ (0, 1) and ν = ρ, • By construction

1 0 v dν t = 1 0 v(T t (x)) dµ(x)
which is convex with respect to t, by convexity of v,

• Similarly [0,1] 2 φ(y, z) dν t (y) dν t (z) = [0,1] 2 φ(T t (x), T t (y)) dµ(x) dµ(y)
is convex with respect to t, by convexity of φ, • The convexity of the remaining congestion term is more involved. Since ν t = T t# µ and T t is non-decreasing, at least formally6 we have ν t (T t (x))T ′ t (x) = µ(x), by the change of variables formula we also have

1 0 F (ν t (y)) dy = 1 0 F (ν t (T t (x)))T ′ t (x) dx = 1 0 F µ(x) T ′ t (x) T ′ t (x) dx
and we conclude by observing that α → F (µ(x)α -1 )α is convex and that T ′ t (x) is linear in t.

Under the assumptions above, J µ is therefore strictly displacement convex and thus admits at most one minimiser (indeed if ν and ρ were different minimisers, by strict displacement convexity, one would have

J µ [ν 1/2 ] < 1 2 J µ [ν] + 1 2 J µ [ρ]
). Actually more is true (see Appendix B.5 for details): if γ is an equilibrium then its second marginal solves (3.5) and therefore is unique. Since c satisfies the generalised Spence-Mirrlees condition (see Corollary 2.3), we deduce the following uniqueness result Theorem 4.2 (Uniqueness of an equilibrium by displacement convexity in dimension one) Under the assumptions above, we have the equivalence ν is a minimiser to (3.5) and γ ∈ Π o (µ, ν)

⇔ γ is an equilibrium and since J µ is strictly displacement convex, there is uniqueness of the equilibrium (which is actually necessarily pure).

Hidden convexity under quadratic cost

The arguments of the previous paragraph can be generalised in higher dimensions when the transport cost is quadratic. Throughout this section, we will assume the following:

• X = Y = Ω where Ω is some open bounded convex subset of R d ,
• µ is absolutely continuous with respect to the Lebesgue measure (that will be the reference measure m 0 from now on) and has a positive density on Ω, • c is quadratic i.e.

c(x, y)

:= 1 2 |x -y| 2 , (x, y) ∈ R d × R d ,
• V again takes the form

V[ν](y) = f (ν(y)) + v(y) + Y φ(y, z) dν(z)
where v is convex, f satisfies the assumptions of Section 3.3 and φ ∈ C(R d × R d ) is symmetric and C 1,1 loc (i.e. C 1 with a locally Lipschitz gradient). Again denoting by F the primitive of f that vanishes at 0, the corresponding energy reads

E[ν] = Y F (ν(y)) dy + Y v(y) dν(y) + 1 2 Y 2 φ(y, z) dν(z) dν(y).
Note that as c is quadratic, Brenier's Theorem (see Theorem A.2) implies the uniqueness and the purity of optimal plans γ between µ and an arbitrary ν. The

J µ [ν] where J µ [ν] := 1 2 W 2 2 (µ, ν) + E[ν]
with W 2 2 (µ, ν) is the squared-2-Wasserstein distance between µ and ν i.e.:

W 2 2 (µ, ν) := inf γ∈Π(µ,ν) X 2 |x -y| 2 dγ(x, y).
Two more structural assumptions are needed to guarantee the strict convexity of J µ along generalised geodesics with base µ (see Appendix A for details), namely McCann's condition:

(4.2) ν → ν d F (ν -d ) is convex non-increasing on (0, +∞)
and that φ is convex. Note that McCann's condition is satisfied by power functions with an exponent larger than 1 as well as by the entropy F (ν) = ν log(ν).

In contrast with Theorem 3.8 where minimisers are equilibria but the reverse is not always true, convexity along generalised geodesics ensures the converse reciprocal property: Theorem 4.3 (Equilibria and minimisers coincide, uniqueness and regularity under generalized convexity) Under the assumptions above, we have the equivalence ν is a minimiser to (4.1) and γ ∈ Π o (µ, ν)

⇔ γ is an equilibrium Moreover, there exists a unique equilibrium (which is actually pure) and the second marginal ν of this equilibrium has a continuous density.

A partial differential equation for the equilibrium

In the quadratic cost framework of Paragraph 4.2, our aim now is to write the optimality condition (3.10) in the form of a nonlinear and non-local equation partial differential equation of Monge-Ampère type. For computational simplicity, we take v = 0 and f (ν) = log(ν) but any convex, C 1,1 loc , symmetric v and any increasing f satisfying McCann's condition would lead to a similar partial differential equation. Let us recall that the unique minimiser/equilibrium ν satisfies (3.10) and is actually characterised by this condition. In this equation, the less explicit term is ϕ c . Thanks to Brenier's theorem (see Appendix A), this term can be made more explicit, as follows: the Brenier map T between µ and ν is the gradient of some convex function, T = ∇u with u convex, and similarly the Brenier map between ν and µ is ∇u * where u * is the Legendre transform of u. In the case of a quadratic cost, u and u * are related to the Kantorovich potential ϕ and its c-transform ϕ c through

(4.3) ϕ(x) = 1 2 |x| 2 -u(x), ϕ c (y) = 1 2 |y| 2 -u * (y), ∀(x, y) ∈ Ω × Ω .
When ∇u is smooth enough, the Monge-Ampère equation, see (A.6), reads:

µ(x) = det(D 2 u(x)) ν(∇u(x)), ∀x ∈ Ω
which has to be supplemented with the natural sort of boundary condition

(4.4) ∇u(Ω) = Ω.
We may then rewrite the optimality condition (3.10) as

(4.5) ν(y) = C exp - 1 2 |y| 2 + u * (y) - Y φ(y, z) dν(z)
where C is a normalisation constant that makes the total mass of the right hand side be 1 on Ω. Since u is defined up to an additive constant, one may actually choose

C = 1. As ∇u # µ = ν, we first have Y φ(∇u(x), z) dν(z) = Ω φ(∇u(x), ∇u(z)) dµ(z) .
On the other hand, using the well-known convex analysis identity

u * (∇u(x)) = x • ∇u(x) -u(x)
and performing the change of variable y = ∇u(x) in (4.5), the Monge-Ampère equation then becomes

(4.6) µ(x) = det(D 2 u(x)) exp - 1 2 |∇u(x)| 2 + x • ∇u(x) -u(x) × exp - Ω φ(∇u(x), ∇u(z)) dµ(z) .
The equilibrium problem is therefore equivalent to a non-local and nonlinear partial differential equation.

This kind of partial differential equation is rather complicated. However, in dimension 1, i.e. when Ω is an open interval, which we can assume to be (0, 1), the boundary condition (4.4) is u ′ (0) = 0, u ′ (1) = 1 and the Monge-Ampère equation (4.6) simplifies to

µ(x) = u ′′ (x) exp - 1 2 u ′ (x) 2 + x • u ′ (x) -u(x) - (0,1) φ(u ′ (x), u ′ (z)) dµ(z) .
Note that this differential equation automatically implies the strict convexity of u. We do not know whether this equation can be solved numerically in an iterative way but will see how the equilibrium can be computed in one dimension thanks to a convenient reformulation of (4.1) as explained below.

Numerical computations in dimension one

Let Ω = (0, 1), m 0 be the Lebesgue measure on (0, 1), µ be absolutely continuous with respect to m 0 with a positive density still denoted µ. Consider again7 the variational problem

(4.7) inf ν J µ [ν]
where

J µ [ν] := 1 2 W 2 2 (µ, ν) + 1 0 ν log ν + 1 2 [0,1] 2 φ(y, z) dν(y) dν(z)
where φ is C 1,1 loc , convex and symmetric. Looking for ν amounts to look for its rearrangement or quantile function:

G(x) := inf{λ : ν([0, λ]) ≥ x}, ∀x ∈ (0, 1).
Note that G is non-decreasing and G # m 0 = ν. We also denote by H the quantile of µ. It is then well known (see [Villani, 2003, Section 2.2] or [Ambrosio et al., 2005, Theorem 6

.0.2]) that W 2 2 (µ, ν) = 1 0 |G(x) -H(x)| 2 dx.
Moreover since G # m 0 = ν, we have

[0,1] 2 φ(y, z) dν(y) dν(z) = [0,1] 2 φ(G(x), G(θ)) dx dθ .
And since ν is regular the change of variable formula yields

1 0 ν(x) log(ν(x)) dx = 1 0 ν(G(x)) log(ν(G(x))) G ′ (x) dx = - 1 0 log(G ′ (x)) dx , as, by definition of G, ν(G(x)) G ′ (x) = 1 a.e.
Therefore reformulating (4.7) in terms of quantile consists in minimising the strictly convex functional 1 2

1 0 |G -H| 2 - 1 0 log(G ′ (x)) dx + 1 2 [0,1] 2 φ(G(x), G(θ)) dx dθ
subject to the boundary conditions G(0) = 0, G(1) = 1. The discretisation of this variational problem is easy to solve using standard gradient descent methods, see Figure 1 8 .

CONCLUDING REMARKS

We conclude the paper, by two remarks, for the sake of simplicity, we adopt exactly the same framework and notations, as in Section 4.2.

Implementation by taxes

By Theorem 4.3 the unique equilibrium is the unique minimiser of the functional J µ . It would therefore be tempting to interpret this result as a kind of welfare theorem. A simple comparison between J µ and the total social cost tells us however that the equilibrium is not efficient. Indeed, the total social cost SC[ν] is the sum of the transport cost W 2 2 (µ, ν)/2 and the additional cost Y V[ν](y) ν(y) dy i.e.

SC[ν] = 1 2 W 2 2 (µ, ν) + Y f (ν(y)) dν(y) + Y v dν + Y 2 φ(y, z) dν(y) dν(z) .
The second term represents the total congestion cost and the last one the total interaction cost. The functional J µ whose minimiser is the equilibrium has a similar form, except that in its second term f (ν)ν is replaced by F (ν) (with F ′ = f ) and the interaction term is divided by 2. The equilibrium corresponds indeed to the case where agents selfishly minimise their own cost c(x, .)

+ V[ν] = c(x, .) + f (ν(.)) + v + Y φ(y, z) dν(z) .
This individual minimisation has of course no reason to correctly estimate the marginal effect of individual behaviour on the total social cost. In other words, there is some gap between the equilibrium and the efficient (social-cost minimising) configurations, and, since we are dealing with a situation with externalities, this is actually not surprising. The computation of the equilibrium and the optimum can be done numerically in dimension 1 by using the same kind of numerical computations as explained in Section 4.4, see Figure 2. The natural way to restore efficiency of the equilibrium is the design by some social planner of a proper system of tax/subsidies which, added to V[ν], will implement the efficient configuration (or at least a stationary point of the social cost). Thanks to our variational approach, a tax system that restores the efficiency is easy to compute (up to an additive constant):

Tax[ν](y) = f (ν(y)) ν(y) -F (ν(y)) + Y φ(y, z) dν(z).
The two terms in Tax[ν] represent respectively a correction to the individual estimation of congestion cost and to the individual estimation of interaction cost. A similar inefficiency of equilibria, arises in the slightly different framework of congestion games, where it is usually referred under the name cost of anarchy, which has been extensively studied in recent years (see [START_REF] Roughgarden | Selfish routing and the price of anarchy[END_REF] and the references therein). In our Cournot-Nash context, we may similarly define the cost of anarchy as the ratio of the worst social cost of an equilibrium to the minimal social cost value:

Cost of anarchy := max{SC[ν e ] : ν e equilibrium} min ν SC[ν] .

In the previous numerical example of Figure 2, both the equilibrium and the optimum are unique and the cost of anarchy can be numerically computed as being approximately 1.8.

A dynamical perspective

Instead of minimising J µ directly, we may think that agents start with some distribution of strategies (that is not an equilibrium) and adjust it with time by a sort of gradient descent dynamics to decrease their individual cost dynamically. At least formally, a way to reach the equilibrium (or minimiser of J µ ) is then to put it into the dynamical perspective of the minimising movement scheme as follows. Fix a time step τ > 0 and start with an initial configuration of strategies ν 0 . The first step of the minimising movement scheme selects a new distribution of strategies ν 1 close to ν 0 (in W 2 ) but also decreasing J µ by

ν 1 ∈ argmin ν 1 2τ W 2 2 (ν 0 , ν) + J µ [ν] .
And then it iterates the process by choosing

(5.1)

ν k+1 ∈ argmin ν 1 2τ W 2 2 (ν k , ν) + J µ [ν] .
This sort of Wasserstein Euler scheme was first introduced in Jordan et al. [1998] for the Fokker-Planck equation. Under suitable conditions it is possible to pass the continuous limit τ → 0 + in the minimising scheme (5.1) and prove that the solution converges in some sense to the continuous evolution equation

     ∂ t ν + div -ν∇ δJ µ δν = 0, ν t=0 = ν 0
which is the gradient flow of J µ in the Wasserstein space (see [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF] for a detailed exposition of the theory). By construction J µ is a Lyapunov function of this equation and even though the equation may have non-unique solutions, by Lyapunov theory, it can be shown under appropriate conditions that its trajectories converge in large time to the unique minimiser of J µ i.e. the equilibrium. If we go back to the individual level, it can be shown that the equation above corresponds to the fact that each agent modifies her strategy according to the gradient flow of her individual cost.

We obtain a sequence of densities which converges to the equilibrium, see

Kantorovich duality

Let X and Y be two compact spaces equipped respectively with the Borel probability measures µ ∈ P(X) and ν ∈ P(Y ). For µ ∈ P(X) and T , Borel: X → Y , T # µ denotes the push forward (or image measure) of µ through T which is defined by A transport map between µ and ν is a Borel map such that T # µ = ν. Now, let c ∈ C(X × Y ) be some transport cost function, the Monge optimal transport problem for the cost c consists in finding a transport T between µ and ν that minimises the total transport cost X c(x, T (x)) dµ(x). A minimiser is then called an optimal transport. Monge problem is in general difficult to solve (it may even be the case that there is no transport map, for instance it is impossible to transport one Dirac mass to a sum of distinct Dirac masses), this is why Kantorovich relaxed Monge's formulation as

T # µ(B) = µ(T -1 (B))
(A.2) Wc(µ, ν) := inf γ∈Π(µ,ν) X×Y c(x, y) dγ(x, y)
where Π(µ, ν) is the set of transport plans between µ and ν i.e. Borel probability measures on X ×Y having µ and ν as marginals. Since Π(µ, ν) is weakly * compact and c is continuous, it is easy to see that the infimum of the linear program defining Wc(µ, ν) is attained at some γ, such optimal γ's are called optimal transport plans (for the cost c) between µ and ν. If there is an optimal γ which is induced by a transport map i.e. is of the form γ = (id, T ) # µ for some transport map T then T is obviously an optimal solution to Monge's problem. Another advantage of the linear relaxation is that it possesses a dual formulation that can be very useful. This dual formulation consists in maximising the linear form X ϕ dµ + Y ψ dν among all pairs (ϕ, ψ) ∈ C(X) × C(Y ) such that ϕ(x) + ψ(y) ≤ c(x, y), it is easy to see that this can be reformulated as a maximisation over ϕ only: Formula (A.3) is usually called Kantorovich duality formula and a maximiser ϕ in (A.3) is called a Kantorovich potential between µ and ν for the cost c. The existence of Kantorovich potentials under our assumptions is well-known (see [START_REF] Rachev | Mass Transportation problems[END_REF], [START_REF] Villani | Topics in optimal transportation[END_REF][START_REF] Villani | Optimal transport: old and new[END_REF]) and we observe that if ϕ is a Kantorovich potential then so is ϕ + C for every constant C.

(A.3) Wc(µ, ν) := sup
We have used in Section 3 the following result on the uniqueness of the Kantorovich potential and the differentiability of Wc(µ, ν) with respect to ν: Lemma A.1 Assume that X = Ω where Ω is some open bounded connected subset of R d with negligible boundary, that µ is equivalent to the Lebesgue measure on X (that is both measures have the same negligible sets) and that for every y ∈ Y , c(., y) is differentiable with ∇xc bounded on X × Y , let ν ∈ P(Y ) then there exists a unique (up to an additive constant) Kantorovich potential ϕ between µ and ν and for every ρ ∈ P(Y ) one has

lim ε→0 + Wc(µ, ν + ε(ρ -ν)) -Wc(µ, ν) ε = Y ϕ c d(ρ -ν).
Proof: The proof of the uniqueness of the Kantorovich potential ϕ between µ and ν up to an additive constant can be found for instance in [Carlier and Ekeland, 2007, Proposition 6.1]. As a normalisation we choose the potential ϕ such that ϕ(x0) = 0 where x0 is some given point of X.

To shorten notations, set νε = ν + ε(ρ -ν), thanks to Kantorovich duality formula (A.3) we have

(A.4) ε -1 [Wc(µ, νε) -Wc(µ, ν)] ≥ Y ϕ c d(ρ -ν)
and similarly if ϕε denotes the Kantorovich potential between µ and νε such that ϕε(x0) = 0, we have

ε -1 [Wc(µ, νε) -Wc(µ, ν)] ≤ Y ϕ c ε d(ρ -ν).
Now it is well-known that (ϕε)ε is bounded and uniformly equi-continuous uniformly with respect to ε hence, thanks to Ascoli's Theorem, up to a sub-sequence, it converges uniformly to some ϕ ∈ C(X) such that ϕ(x0) = 0 and it is easy to see that ϕ is a Kantorovich potential between µ and ν so that ϕ = ϕ and ϕ c ε converges to ϕ c . We then have lim sup

ε→0 + ε -1 [Wc(µ, νε) -Wc(µ, ν)] ≤ Y ϕ c d(ρ -ν).
The desired result thus follows from (A.4).

Q.E.D.

When X = Y and denoting by d the distance on Y , for p ∈ [1, +∞[, the p-Wasserstein distance between µ ∈ P(X) and ν ∈ P(X) is by definition

(A.5) Wp(µ, ν) := inf γ∈Π(µ,ν) X×Y d(x, y) p dγ(x, y) 1/p
The Wasserstein distances are indeed distances and they metrise the weak * topology of P(Y ). For p = 1, it is well-known that the Kantorovich duality formula can be rewritten as

W1(µ, ν) = sup X ϕ d(µ -ν) : ϕ 1-Lipschitz
so that for every Lipschitz continuous function ϕ on X, one has

X ϕ d(µ -ν) ≤ Lip(ϕ, X)W1(µ, ν),
an inequality we will use several times later on. As a simple illustration of the interest of the distance W1, let us equip Y m with the distance (x, y) = (x1, . . . , xm, y1, . . . ym) → dm(x, y) := m k=1 d(x k , y k ) for ν and θ in P(Y ), let γ be an optimal transport plan between µ and ν for W1 then since γ ⊗m := γ ⊗ • • • ⊗ γ has marginals ν ⊗m and θ ⊗m , we have

W1(ν ⊗m , θ ⊗m ) ≤ Y m ×Y m dm dγ ⊗m = mW1(ν, θ) .
Which shows in particular that if (νn)n weakly * converges to ν then (ν ⊗m n )n weakly * converges to

ν ⊗m i.e. Y m ϕ dν ⊗m n converges to Y m ϕ dν ⊗m as n → ∞ for every φ ∈ C(Y m ).
Of particular interest is also the quadratic case p = 2 in an euclidean setting for which a brief summary of the main results used in the paper is given in the next paragraphs.

The quadratic case and Monge-Ampère equation

We now restrict ourselves to the quadratic case, the solution of the quadratic optimal transport problem is due to Yann Brenier whose path-breaking paper [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF] totally renewed the field of optimal transport and was the starting point of an extremely active stream of research since the 90's.

Theorem A.2 (Brenier's theorem) Let µ ∈ P(R d ) be absolutely continuous with respect to the Lebesgue measure and compactly supported and ν ∈ P(R d ) be compactly supported, then the quadratic optimal transport problem

W 2 2 (µ, ν) := inf γ∈Π(µ,ν) R d ×R d |x -y| 2 dγ(x, y)
possesses a unique solution γ which is in fact a Monge solution γ = (id, T ) # µ. Moreover T = ∇u µ-a.e. for some convex function u and ∇u is the unique (up to µ-a.e. equivalence) gradient of a convex function transporting µ to ν; T = ∇u is called the Brenier map between µ and ν.

In fact the previous theorem holds under much more general assumptions (it is enough that µ and ν have finite second moments and that µ does not charge sets of Hausdorff dimension less than d -1, see [START_REF] Mccann | Existence and uniqueness of monotone measure-preserving maps[END_REF] or [START_REF] Villani | Topics in optimal transportation[END_REF]). Brenier's theorem roughly says that there is a unique optimal transport for the quadratic cost and that it is characterised by the fact that it is of the form ∇u with u convex, in other words, solving ∇u # µ = ν with u convex determines ∇u uniquely µ-a.e.. When we have additional regularity, i.e. when µ and ν have regular densities (still denoted µ and ν) and ∇u is a diffeomorphism between the support of µ and that of ν, thanks to the change of variables formula, we find that u solves the Monge-Ampère partial differential equation:

(A.6) µ = ν(∇u) det(D 2 u).
A deep regularity theory due to Luis Caffarelli: Caffarelli [1992a,b] implies that the Brenier map is a smooth diffeomorphism when in addition µ and ν are smooth, bounded away from 0 and have convex supports, in particular the Monge-Ampère equation is satisfied in this case which justifies the computations of Section 4.3.

where ϕ is a Kantorovich potential between µ and ν (it is unique up to an additive constant by Lemma A.1 and this constant plays no role since ρ and ν have the same mass). Minimising the left-hand side of (B.1), with respect to L p probabilities ρ, this yields that ν-a.e.

ϕ c + V[ν] = inf ρ∈D Y (ϕ c + V[ν]) dρ = M
where M := Essinf(ϕ c +V[ν]) denotes the essential infimum of (ϕ c +V[ν]) i.e. the largest constant that bounds (ϕ c +V[ν]) from below m0-a.e.. Since γ is an optimal transport plan we have c(x, y) = ϕ(x) + ϕ c (y) γ-a.e. whereas by definition c(x, z) ≥ ϕ(x) + ϕ c (z) for all (x, z) ∈ X × Y . We thus have

c(x, z) + V[ν](z) ≥ M + ϕ(x) for all x ∈ X and m0-a.e. z ∈ Y c(x, y) + V[ν](y) = M + ϕ(x)
for γ-a.e. (x, y).

This proves that γ is a Cournot-Nash equilibrium.

B.2. Proof of Corollary 3.3

Thanks to Theorem 3.2, it is enough to prove that (3.5) admits solutions and to recall that the set Πo(µ, ν) is nonempty. Let (νn)n be a minimising sequence of (3.5). Thanks to the growth condition (3.2), (νn)n is bounded in L p (m0). It thus admits a (not relabelled) sub-sequence that converges weakly in L p (m0) (and thus in particular weakly * in P(Y )) to some ν ∈ L p (m0). By the convexity of F (y, .), the first term in E is lower-semi continuous for the weak topology of L p (m0). By the continuity of φ the second term in E is continuous for the weak- * topology of P(Y ). Finally, the lower-semi continuity of Wc(µ, .) for the weak- * topology straightforwardly follows from the Kantorovich duality formula (A.3). We thus have inf

ν Jµ(ν) = lim inf n {Wc(µ, νn) + E [νn]} ≥ Wc(µ, ν) + E [ν] .
So that ν solves (3.5).

B.3. Proof of Proposition 3.5

Assume that γ is an equilibrium and let ν be its second marginal. Let then ϕ be a Kantorovich potential between µ and ν such that By convexity of E and (3.4), we obtain

E [ρ] -E [ν] ≥ Y V[ν] d(ρ -ν)
hence, finally using (B.2) and the fact that ρ is absolutely continuous with respect to m0, we get The existence of a minimiser is similar to the proof of Corollary 3.3 since the coercivity of F ensures that minimising sequences are uniformly integrable and its convexity guarantees sequential weak lower semi continuity of Jµ.

Let ν solve (3.5) and let us prove that it is bounded away from zero. Let us assume by contraction that m0({ν ≤ λ}) > 0 for every λ > 0. Let δ0 > 0, δ ∈ (0, δ0) (to be chosen later on). Let A := {x : δ0 ≤ ν(x) ≤ M } where M > 0 is large enough so that m0(A) > 0. For small ε > 0 such that ε < δm0(A)/2 then define νε := ν + ε(uA δ -uA)

where A δ := {ν ≤ δ} and for m0(B) > 0, uB denotes the (sort of uniform probability on B) uB := m0(B) -1 1lB. Since ε < δm0(A)/2 and ν > δ on A, νε is a probability measure. By optimality of ν, we then have Denoting by ϕε a Kantorovich potential between µ and νε, we first have

Wc(µ, νε) -Wc(µ, ν) ≤ ε Y ϕ c ε d(uA δ -uA) .
And since ϕ c ε has a modulus of continuity that is uniform with respect to ε (that of c) and ϕ c ε can be normalised so as to vanish at the same point, ϕ c ε is uniformly bounded independently of ε and δ. So that Wc(µ, ν) -Wc(µ, νε) ≤ C1ε for some constant C1. In a similar way, one finds a constant C2 such that for ε small enough and uniformly in δ one has since F is Lipschitz on [δ0/2, M ] the second term can be bounded from above by C3ε for a constant C3 again independent of δ and ε. Now let C := C1 + C2 + C3 thanks to Inada's condition there is some δ1 < δ0/2 such that f ≤ -C -1 on (0, 2δ1]. Choosing δ ≤ δ1 and ε small enough so that εm0(A δ ) -1 ≤ δ1, we then have

A δ F (ν + εm0(A δ ) -1 ) -F (ν) dm0 ≤ (-C -1)ε.
Putting everything together, the latter inequality gives the desired contradiction to (B.3).

The proof of the upper bound is similar: one assumes that ν / ∈ L ∞ (m0) and then considers a perturbation of the form νε := ν +ε(uC -uC M ) with CM := {ν ≥ M }, M large and C well chosen, the computations are the same as before and the contradiction comes from the Inada condition at +∞: F ′ (ν) = f (ν) → +∞ as ν → +∞.

B.5. Proof of Theorem 4.3

Uniqueness of a minimiser follows directly from the strict convexity along generalised geodesics with base µ of Jµ which follows from McCann's condition (4.2), the convexity of v and φ and the strict convexity along generalised geodesics with base µ of W 2 2 (µ, .). Let us now assume that γ is an equilibrium and that ν is its second marginal, then, for some constant M , we have: and is therefore continuous.

Let us now prove that ν solves (4.1), let ρ be another probability measure (which we can assume to have a positive and continuous density as well), and let t ∈ [0, 1] → νt denote the generalised geodesic with base µ joining ν and ρ, i.e. νt = Tt # µ := ((1 -t)T0 + tT1) # µ where T0 (resp. T1) denotes the Brenier map between µ and ν (resp. ρ). Since (T0, T0 + t(T1 -T0)) # µ has marginals ν and νt we have By the convexity of Jµ along generalised geodesics with base µ, setting g(t) = Jµ[νt] and using g(t) ≤ (1 -t)g(0) + tg(1) for all t ∈ (0, 1) we have:

X

  ϕ(x) dµ(x) + Y ϕ c (y) dν(y) where ϕ c denotes the c-transform of ϕ i.e. (2.3) ϕ c (y) := min x∈X {c(x, y) -ϕ(x)} .

  Theorem 3.2 (Minimisers are equilibria) Assume that V[ν] satisfies (3.4) with D = P(Y ) ∩ L p (m 0 ) for some p ∈ [1, +∞[ and that the assumptions of Lemma A.1 hold true. If ν solves (3.5) and γ ∈ Π o (µ, ν) then γ is a Cournot-Nash equilibrium.

Figure 1 .-

 1 Figure 1.-The distribution µ of the agents is dash line and the solution ν to (4.7) in the case f (x) = x 8 , φ(z) = 10 -4 |z| 2 and v = (x -10) 4 .

Figure 2 .-

 2 Figure 2.-The optimum in continuous line and the equilibrium in dash line on the left. The corresponding taxes on the right.

  Figure 3. The descent algorithm is very fast and the computed equilibrium is very stable with respect to the initial density for the gradient descent as shown in the left hand figure of Figure 4.

Figure 3 .-Figure 4 .-

 34 Figure 3.-Convergence and stabilisation toward the equilibrium in the case of a logarithmic congestion, cubic interaction, and a potential v(x) := (x -5) 3 with 1l [0,1] as initial guess. The inverse of the cumulative function on the left and the corresponding density on the right.

  for every Borel subset B of Y or equivalently by the change of variables formula (A.1) Y ϕ dT # µ = X ϕ(T (x)) dµ(x), ∀ϕ ∈ C(X).

  where ϕ c is the c-concave transform of ϕ i.e. ϕ c (y) := min x∈X {c(x, y) -ϕ(x)}, ∀y ∈ Y.

  (B.2) ϕ c + V[ν] ≥ 0 m0-a.e. ϕ c + V[ν] = 0 ν-a.e.Let ρ ∈ D, thanks to the Kantorovich duality formula (A.3), we first haveWc(µ, ρ) -Wc(µ, ν) ≥ Y ϕ c d(ρ -ν) .

  + V[ν]) d(ρ -ν) ≥ 0 which means that ν solves (3.5). A. BLANCHET & G. CARLIER B.4. Proof of Lemma 3.7

  (B.3) 0 ≤ Jµ[νε] -Jµ[ν] = Wc(µ, νε) -Wc(µ, ν) + E [νε] -E [ν].

F

  z) dνε(y) d(νε -ν)(z) ≤ C2ε.Now it remains to estimate the last term namelyY [F (νε) -F (ν)] dm0 = A δ F (ν + εm0(A δ ) -1 ) -F (ν) dm0 + A (ν -εm0(A) -1 ) -F (ν) dm0

(

  z) dν(z) + ϕ c ≥ M a.e. with an equality ν-a.e.Thanks to Inada's condition and the fact that the right hand side is continuous, this implies ν is bounded away from zero so that (B.4) actually is an equality (Lebesgue) almost everywhereV[ν] + ϕ c = M and thus ν satisfies ν(y) = f -1 M -v(y) -Y φ(y, z) dν(z) -ϕ c (y)

(

  dµ ≤ t diam(Y ) .

  νt] -Jµ[ν]).Let us write νt as νt := ν + tht, by the (usual) convexity of F and that of W 2 2 (µ, .), ν) + ϕ c ]ht dm0Let us now expand Y 2 φ(y, z) dνt(y) dνt(z) in powers of t as 1 2tY 2 φ(y, z) d [(ν + tht)(y)(ν + tht)(z) -ν(y)ν(z)] = Y 2 φ(y, z) dht(y) dν(z) + Rt z) d(νt-ν)(y) d(νt-ν)(z) ≤ 1 t W1(νt, ν)×Lip(ψt, Y ) ≤ diam(Y ) Lip(ψt, Y )where the last inequality follows from (B.5) and ψt is defined byψt(y) := Y φ(y, z) d(νt -ν)(z) = Y (φ(y, Tt(x)) -φ(y, T0(x))) dµ(x) .Since ∇φ is locally Lipschitz and Tt -T0 is uniformly bounded by t diam(Y ) we find that there is a constant C such that Lip(ψt, Y ) ≤ Ct so that Rt = O(t). Putting everything together and using (z) dν(z) + ϕ c (y)) dht(y) + Rt = Y (V[ν] + ϕ c ) dht + Rt ≥ M Y dht + Rt = Rt → 0 as t → 0 +

The role of the reference measure m0 is here to capture purely local congestion effects as in the examples below. In other words, we will require action distributions to be absolutely continuous with respect to m0. This departs from the common assumption that the cost is well defined for every action distribution and satisfies some strong continuity/semi-continuity with respect to the weak * topology of measures as in[START_REF] Mas-Colell | On a theorem of Schmeidler[END_REF] or[START_REF] Kahn | On Cournot-Nash equilibrium distributions for games with a nonmetrizable action space and upper semi-continuous payoffs[END_REF].

Note the analogy with the variational approach of[START_REF] Monderer | Potential games[END_REF] for potential games, i.e. games whose equilibria can be obtained by minimising some potential function.

Let us remark that in the case of a finite number of players, the role of symmetry for the potential approach to work was already pointed out in LeBreton and Weber[2011].

see[START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF],[START_REF] Villani | Topics in optimal transportation[END_REF][START_REF] Villani | Optimal transport: old and new[END_REF] for a rigorous justification.

Again the choice of a logarithmic congestion function is not so essential, power congestion functions can be considered as well.

Actually, in our simulations, we have relaxed the condition that the support of ν is [0, 1] i.e. the boundary conditions G(0) = 0, G(1) = 1.

Acknowledgements The authors acknowledge the support of the Agence Nationale de la Recherche through the projects ANR-09-JCJC-0096-01 EVaMEF and ANR-07-BLAN-0235 OTARIE. The authors wish to thank Jocelyn Donze, André Grimaud, Michel Le Breton, Jérôme Renault, François Salanié and the participants to the IAST LERNA -Eco/Biology and Chicago University Seminars for many interesting and fruitful discussions about the present work.

APPENDIX A: THE OPTIMAL TRANSPORT TOOLBOX

This appendix just gives some basic results from optimal transport theory that we have used in the paper, for a detailed exposition of this rich and rapidly developing subject, we refer the interested reader to the very accessible textbook [START_REF] Villani | Topics in optimal transportation[END_REF] or [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF], [START_REF] Villani | Optimal transport: old and new[END_REF] or, the more probability-oriented textbook [START_REF] Rachev | Mass Transportation problems[END_REF].

Convexity along generalised geodesics

The last ingredient from optimal transport theory that we have used (in Section 4) is the powerful notion of displacement convexity along generalised geodesics due to Ambrosio, Gigli-Savaré [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF] 9 . As in Section 4.2, we assume that X = Y = Ω where Ω is some open bounded convex subset of R d , that the cost is quadratic, that m0 is the Lebesgue measure on X and that µ is absolutely continuous with respect to m0 and has a positive density on Ω. In particular for every ν ∈ P(X), the Brenier's map between µ and ν is well-defined. Generalised geodesics with base µ for the Wasserstein distance W2 and the corresponding notion of convexity are defined as follows Definition A.3 (Convexity along generalised geodesics) Let ν ∈ P(X), ρ ∈ P(X), let T0 be the Brenier's map between µ and ν and let T1 be Brenier's map between µ and ρ, the generalised geodesic with base µ between ν and ρ is the curve of measures t ∈ [0, 1] → νt := ((1-t)T0+tT1) # µ. The functional J : P(X) → R ∪ {+∞} is called convex along generalised geodesics with base µ if for every pair of endpoints ν and ρ in P(X) and for every t ∈ [0, 1], one has

If, in addition, the previous inequality is strict for t ∈ (0, 1) and ρ = ν, J is called strictly convex along generalised geodesics with base µ.

In Section 4.2, we were interested in the strict convexity along generalised geodesics with base µ of the functional Jµ defined by ( 4

) is convex non-increasing on (0, +∞), it follows from [Ambrosio et al., 2005, Proposition 9.3.9]. Finally, in the functional Jµ defined by (4.1), we had the term W 2 2 (µ, .), to see that it is strictly convex along generalised geodesics with base µ, we can proceed exactly as we did for Wc(µ, .) in dimension one in Paragraph 4.1.

APPENDIX B: PROOFS OF THE RESULTS

B.1. Proof of Theorem 3.2

Let ν be a solution of (3.5), ρ ∈ D and ε ∈ (0, 1), we then have ε

Using the fact that V[ν] is the first variation of E and Lemma A.1, we thus get

9 Actually this notion of convexity is a slight variant of the notion of displacement convexity which first appeared in the seminal work of [START_REF] Mccann | A convexity principle for interacting cases[END_REF]. It is known that W 2 2 (µ, ν), as a function of ν is not displacement convex in the sense of McCann (see example 9.1.5. in [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF]) and this is the very reason why, following Ambrosio, Gigli and Savaré, we consider convexity along generalised geodesics with base µ rather than the initial notion of McCann. Let us however indicate that, in dimension one, both notions coincide.


Mathieu Giraud

Richard Groult

Florence Levé

Subject and counter-subject detection for analysis of the Well-Tempered Clavier fugues

Keywords: symbolic music analysis, contrapuntal music, fugue analysis, repeating patterns

come

Introduction

Contrapuntal music is a polyphonic music where each individual line bears interest in its own. Bach fugues are a particularly consistent model of contrapuntal music. The fugues of Bach's Well-Tempered Clavier are composed of two to five voices, appearing successively, each of these voices sharing the same initial melodic material: a subject and, in most cases, a counter-subject. These patterns, played completely during the exposition, are then repeated all along the piece, either in their initial form or more often altered or transposed, building a complex harmonic network.

To analyze symbolic scores with contrapuntal music, one can use generic tools detecting repeating patterns or themes, possibly with approximate occurrences. Similarity between a pattern and several parts of a piece may be computed by the Mongeau-Sankoff algorithm [START_REF] Mongeau | Comparaison of musical sequences[END_REF] and its extensions or by other methods for approximate string matching [START_REF] Clifford | Approximate string matching for music analysis[END_REF][START_REF] Crawford | String matching techniques for musical similarity and melodic recognition[END_REF], allowing a given number of restricted mismatches. Several studies focus on finding maximal repeating patterns, limiting the search to non-trivial repeating patterns, that is discarding patterns that are a sub-pattern of a larger one with the same frequency [START_REF] Hsu | Efficient repeating pattern finding in music databases[END_REF][START_REF] Karydis | Finding maximum-length repeating patterns in music databases[END_REF][START_REF] Liu | Efficient theme and non-trivial repeating pattern discovering in music databases[END_REF][START_REF] Yu | Fault tolerant non-trivial repeating pattern discovering for music data[END_REF]. Other studies try to find musically significant themes, with algorithms considering the number of occurrences [START_REF] Smith | Discovering themes by exact pattern matching[END_REF], but also the melodic contour or other features [START_REF] Meek | Automatic thematic extractor[END_REF]. Some MIR studies already focused on contrapuntal music. The study [START_REF] Weng | Automatic musical form analysis[END_REF] builds a tool to decide if a piece is a fugue or not, but no details are given on the algorithm. The bachelor thesis [START_REF] Browles | Creating a tool to analyse contrapuntal music[END_REF] contains a first approach to analyze fugues, including voice separation. For sequence analysis, it proposes several heuristics to help the selection of repeating patterns inside the algorithms of [START_REF] Hsu | Efficient repeating pattern finding in music databases[END_REF] which maximizes the number of occurrences. The website [START_REF] Hakenberg | The Pirate Fugues[END_REF] also produces an analysis of fugues, extracting sequences of some repeating patterns, but without precise formal analysis nor precise bounds.

One can take advantage of the apparently simple structure of a fugue: as the main theme -the subject -always begins at only one voice, this helps the analysis. But a good understanding of the fugue requires to find where the subject exactly ends. In this work, we start from a symbolic score which is already track-separated, and we propose an algorithm to sketch the plan of the fugue. The algorithm tries to retrieve the subjects and the counter-subjects, precisely determining the ends of such patterns. We tested several substitution functions to have a sensible and specific approximate matching. Our best results use a simple diatonic similarity between pitch intervals [START_REF] Cambouropoulos | A general pitch interval representation: Theory and applications[END_REF] combined with a strict length matching for all notes, except for the first and the last one.

The paper is organized as follows. Section 2 gives definitions and some background on fugues, Section 3 details the problem of the bounds of such patterns, Section 4 presents our algorithm, and Section 5 details the results on the 24 fugues of the first book of Bach's Well-Tempered Clavier. These results were evaluated against a reference musicological book [START_REF] Siglind Bruhn | s Well-Tempered Clavier[END_REF]. The algorithm predicts two thirds of the subjects with a musically relevant end, and finally retrieves 85% of the subject occurrences, with almost no false positives.

Preliminaries

A note x is described by a triplet (p, o, ), where p is the pitch, o the onset, and the length. The pitches can describe diatonic (based on note names) or semitone information. We consider ordered sequence of notes x 1 . . . x m , that is 1). The sequence is monophonic if there are never two notes sounding at the same onset, that is, for every i with 1 ≤ i < m, o i + i ≤ o i+1 . To be able to match transposed patterns, we consider relative pitches, also called intervals: the interval sequence is defined as ∆ x 2 . . . ∆ x m , where ∆ x i = ( ∆ p i , o i , i ) and

x 1 = (p 1 , o 1 , 1 ), . . . , x m = (p m , o m , m ), where 1 ≤ o 1 ≤ o 2 ≤ . . . ≤ o m (see Fig.
∆ p i = p i -p i-1 .
We now introduce some notions about fugue analysis (see for example [START_REF] Siglind Bruhn | s Well-Tempered Clavier[END_REF][START_REF] Prout | Analysis of J.S. Bach's forty-eight fugues (Das Wohltemperirte Clavier)[END_REF] for a complete musicological analysis). These concepts are illustrated by an example on Fugue #2, which has a very regular construction.

A fugue is given by a set of voices, where each voice is a monophonic sequence of notes. In Bach's Well-Tempered Clavier, the fugues have between 2 and 5 voices, and Fugue #2 is made of 3 voices.

The fugue is built on a theme called subject (S). The three first occurrences of the subject in Fugue #2 are detailed in Fig. 2: the subject is exposed at one voice (the alto), beginning by a C, until the second voice enters (the soprano, measure 3). The subject is then exposed at the second voice, but is now transposed to G. Meanwhile, the first voice continues with the first counter-subject (CS) that combines with the subject. Fig. 3 shows a sketch of the entire fugue. The fugue alternates between other instances of the subject together with counter-subjects (8 instances of S, 6 instances of CS, and 5 instances of the second counter-subject CS2) and development on these same patterns called episodes (E).

All these instances are not exact ones -the patterns can be transposed or altered in various ways. As an example, Fig. 4 shows the five complete occurrences of CS. For these occurrences, the patterns can be (diatonically) transposed, and the lengths are conserved except for the first and last note. 3 Where does the subject end?

A fundamental question concerns the precise length of the subject and of any other interesting pattern. The subject is heard alone at the beginning of the first voice, until the second voice enters. However, this end is generally not exactly at the start of the second voice. Formally, let us suppose that the first voice is x 1 , x 2 , ..., and the second one is y 1 , y 2 , ..., with x i = (p i , o i , i ) and y j = (p j , o j , j ). Let x z be the last note of the first voice heard before or at the start of the second voice, that is

z = max{i | o i ≤ o 1 }.
The end of the subject is roughly at x z . Table 2, at the end of the paper, lists the exact values of g such that the true subject ends at x z+g : Fig. 4. The 5 complete occurrences of the first counter-subject into Fugue #2 in C minor (BWV 847). (Note that this counter-subject actually has a latter occurrence, split between two voices.) In these occurrences, all notes -except the first and the last ones -have exactly the same length. The values in the occurrences indicate the intervals, in number of semitones, inside the counter-subject. Only occurrences #2 and #5 have exactly the same intervals. The occurrence #4 is almost identical to occurrence #1, except that it lacks the octave jump (+3 instead of +15). Between groups {#1, #4}, {#2, #5}, and {#3}, the intervals are not exactly the same. However, all these intervals (except the lack of the octave jump in #4) are equal when one considers only diatonic information (bottom small staff): clef, key and alterations are here deliberately omitted, as semitone information is not considered.

in the first book of Bach's Well-Tempered Clavier, we notice that g is always between -8 and +6, and, in the majority of cases, between -4 and +1 .

For example, in the Fugue #2, the subject has 20 notes, ending on alto note x 20 (the first sixteenth of the third measure, E , first circled note on Fig. 2), that is 2 notes before the start of the soprano voice (g = -2). This can be deduced from many observations in this third measure:

metrically, the phrase ends on a strong beat; harmonically, the five preceding notes "F G A G F" suggest a 9th dominant chord, which resolves on the E suggesting the C minor tonic; moreover, the subject ends with a succession of sixteenths with small intervals, whereas the following note x 21 (C) belongs to CS with the line of falling sixteenths.

Determining the precise end of the subject is thus an essential step in the analysis of the fugue: it will help to localize the counter-subject and build the structure with all occurrences of these patterns, but also to understand the rhythm, the harmony and the phraseology of the whole piece.

We could use generic algorithms to predict the subject end. For example, the "stream segment detection" described in [START_REF] Rafailidis | Detection of stream segments in symbolic musical data[END_REF] considers melody, pitch and rhythm informations. Many different features are also discussed in [START_REF] Meek | Automatic thematic extractor[END_REF] for theme extraction. However, in the following, we will show that a simple algorithm, only based on similarities, is able to detect precisely most of the subject ends.

Algorithm

Starting from track-separated data, we propose here to detect the subject as a repeating pattern finishing approximatively at the start of the second voice, under a substitution function considering a diatonic similarity for pitch intervals, and enforcing length equalities of all notes except the first one and the last one.

The similarity score between a pattern and the rest of the fugue piece can be computed via dynamic programming by the Mongeau-Sankoff equation [START_REF] Mongeau | Comparaison of musical sequences[END_REF]. The alignment can then be retrieved through backtracking in the dynamic programming table .  As almost all the content of a fugue is somewhat derived from a subject or some counter-subject, any part will match a part of the subject or of another base pattern within a given threshold. Here, we will use very conservative settingsonly substitution errors, and strict length requirements -to have as few false positives as possible, still keeping a high recognition rate.

Subject identification. To precisely find the end the subject, we thus want to test patterns finishing at notes x z+g , where g ∈ [g min , g max ] = [8, +6]. Each one of these candidates is matched against all the voices. In this process, we use a substitution cost function able to match the first and the last notes of the subject independently of their lengths.

Let S(a, b) be the best number of matched intervals when aligning the start x 1 . . . x a of a pattern (the subject) against a part of a given voice finishing at y b , and S f (a, b) the best number of matched intervals when aligning a complete pattern (the complete candidate subject) x 1 . . . x a against the same part. These tables S and S f may be computed by the following dynamic programming equation:

         S(1, b) = 0 ∀a ≥ 2, S(a, b) = S(a -1, b -1) + δ( ∆ x a , ∆ y b ) (match, substitution) ∀a ≥ 2, S f (a, b) = S(a -1, b -1) + δ f ( ∆ x a , ∆ y b ) (finishing)
The substitution functions δ and δ f are detailed on Fig. 5: δ checks pitch intervals and lengths, whereas δ f only considers pitch intervals. The length of the first notes (x 1 and y 1 ) is neither checked, as the algorithm actually compares ∆ x 2 . . . ∆ x a against ∆ y 2 . . . ∆ y b . Finally, notice that these equations only use substitution operations, but can be extended to consider other edit operations. As in [START_REF] Hsu | Efficient repeating pattern finding in music databases[END_REF], we only compute once each table (for a given voice), then we scan the table S f to find the occurrences: given a sequence x and a threshold τ , the candidate finishing at x z+g occurs in the sequence y if for some position i in the text, S f (z + g, i) ≥ τ . The best candidate is selected on the total number of matched intervals in all occurrences. We call x s its last note, so the subject is defined to be x 1 ...x s . The whole algorithm is in O(mn), where m = z + g max .

δ(( ∆ p, o, ), ( ∆ p , o , )) =    +1 if ∆ p ≈ ∆ p and = 0 if ∆ p ≈ ∆ p and = -∞ otherwise δ f (( ∆ p, o, ), ( ∆ p , o , )) = +1 if ∆ p ≈ ∆ p 0 otherwise
For example, on the Fugue #2, the algorithm correctly selects the note x 20 as the end of the subject (see Table 1). The score is the sum of the S f (z + g, i) values at least equal to τ : it is the total number of intervals exactly matched on all occurrences. Here this end corresponds to the "nontrivial maximal-length repeating pattern" for most occurrences, but it is not always the case.

Interval similarities and diatonic matching. The equation Fig. 5 needs to have a similarity relation ≈ on pitch intervals. Between a strict pitch equality and very relaxed "up/down" classes defining the contour of some melody [START_REF] Ghias | Query by humming: musical information retrieval in an audio database[END_REF], some intermediary interval classes may be defined as "step/leap intervals" [START_REF] Cambouropoulos | A pattern extraction algorithm for abstract melodic representations that allow partial overlapping of intervallic categories[END_REF] or "quantized partially overlapping intervals" (QPI) [START_REF] Lemström | Musical information retrieval using musical parameters[END_REF].

We propose here to use a similarity on diatonic pitches. Such a pitch representation is often mentioned [START_REF] Orpen | Measurement of similarity in music: A quantitative approach for non-parametric representations[END_REF][START_REF] Perttu | Combinatorial pattern matching in musical sequences[END_REF] and was studied in [START_REF] Cambouropoulos | A general pitch interval representation: Theory and applications[END_REF][START_REF] Cambouropoulos | Influence of musical similarity on melodic segmentation: Representations and algorithms[END_REF][START_REF] Hiraga | Structural recognition of music by pattern matching[END_REF]. A diatonic model is very relevant for tonal music: it is sensible enough to allow mode changes, while remaining specific enough -a scale will always match only a scale. For example, with diatonic similarity, all occurrences but one of the counter-subject on the Fig. 4 can be retrieved exactly, and the occurrence #4 with only one substitution.

Counter-subjects identification. The same method as for subject identification is used to retrieve the first counter-subject. The first occurrence of the countersubject starts right after the subject (at x s+1 ), and its length is approximatively equal to the length of the subject. We thus have a rough end of the countersubject at x w , where w = max{i | o i -o s+1 ≤ o s -o 1 }, and the same procedure refines the bound to find the good last note x cs , where again cs -w is in a given interval. To prevent detection of non-relevant patterns, the counter-subject is marked as not detected if the above procedure leads to more occurrences than the subject occurrences.

Results and discussion

We tested the algorithm of the previous section on the 24 fugues of the first book of Bach's Well-Tempered Clavier, starting from Humdrum files where the voices are separated, available for academic purposes at http://kern.humdrum.org/. The pitches were encoded according to two frameworks: MIDI encoding, and Base40 encoding [START_REF] Walter | A base-40 number-line representation of musical pitch notation[END_REF]. While the first one only counts semitones, the second one allows to discriminate enharmonic pitches, thus allowing a precise diatonic match as described in the previous section.

We ran the algorithm on the 24 fugues3 , and manually checked all results and occurrences. Results (with diatonic similarity) are summarized on Table 2. We fixed a minimum threshold of τ = 0.9z -3, where z is the number of notes defined in Section 3.

Subject lengths. We searched for end of subjects in the range [g min , g max ] = [-8, +6], that are the observed values. In 16 of the 24 fugues, the algorithm retrieves precisely the ends of the subjects. To our knowledge, this is the first algorithm able to correctly detect the ends of the subject: In [START_REF] Browles | Creating a tool to analyse contrapuntal music[END_REF], the subjects found are said to be "missing or including an extra 1 to 4 notes", and the ends of the subjects on [START_REF] Hakenberg | The Pirate Fugues[END_REF] are also very approximate.

Fugue #8 shows why the proposed algorithm does not always find the correct length of the subject. In this fugue, a subject of length 9 notes is found instead of 13 notes: there are several truncated occurrences of the subject, and the algorithm chooses the end that provides the best match throughout the piece (Fig. 6).

The algorithm already considers the last note in a special way (and the former notes can be handled through substitution errors in the pitch intervals). It is possible to adapt the matching to be even more relaxed towards the end of the pattern, but we did not see a global improvement in the detection of subject lengths.

False positives. There are very few false positives among the subjects found (specificity of 90%), even when the length of the subject is badly predicted. The false positives appear in only two fugues: Fig. 6. Some subject occurrences in Fugue #8 in D# minor. The occurrence #1 is the first one, and is similar to 16 occurrences, sometimes with diatonic transpositions. In the occurrence #2, the last but one note of the subject (circled E) has not the same length than in the other occurrences (and this is forbidden by our substitution function δ). In the occurrence #3, a supplementary note (circled G) is inserted before the end of the subject, again preventing the detection if the true length of the subject is considered. Moreover, the occurrences #3, #4 and #5 are truncated to the head of the subject, and lead to a false detection of subject length.

in Fugue #19, the 5 false positives correspond to 4 extended subjects [START_REF] Siglind Bruhn | s Well-Tempered Clavier[END_REF],

and one almost complete subject. in Fugue #5, the length of the subject is wrongly selected to the first 9 notes (8 first thirty-second notes and a final note), and this head of the subject matches the 11 true occurences, but also 24 false positives.

False negatives. The algorithm correctly retrieves about 85% of the subject occurrences. The false negatives are occurrences that are too much altered: insertions, deletions, or too many substitutions compared to the threshold.

Inverted and augmented subjects. In some fugues, the subject appears upside down (all its intervals are reversed) or augmented (all lengths are doubled). Once the subject is known, the same matching algorithm can thus be applied to the inversion or the augmentation of the subject. This method never produced a false positive, and was able to recover 72% (26/36) of the complete inverted and augmented subjects reported in [START_REF] Siglind Bruhn | s Well-Tempered Clavier[END_REF].

Counter-subjects. Counter-subjects were detected with the same algorithm within the range [g min , g max ] = [-2, +4]. In 40% of the fugues, the algorithm correctly detects the exact length of the CS or the absence of a CS.

In 9 fugues, the algorithm predicts the absence of CS. This was expected for Fugues #1, #8 (no CS), #15 (the CS occurs completely only once) #19 (late exposition of CS) and #20 (there is no real "characteristic and independent counter-subject" according to [START_REF] Siglind Bruhn | s Well-Tempered Clavier[END_REF]). As in the case of the subjects, there are false negatives due to the bad recogntion of altered patterns. Moreover, when the subject is badly detected, the detection of the counter-subject end fails in the majority of the cases.

The algorithm retrieves correctly about the half of the CS occurrences, with more than 80% specificity. Pitch interval similarities. We compared the diatonic matching against a simple exact matching on MIDI semitones, possibly adapting the error threshold. As expected, diatonic similarity has a better performance, because such a relaxed similarity is able to match approximate occurrences as the countersubjects shown on Fig. 4.

Starting from MIDI pitches, an idea could be thus to use pitch spelling methods as [START_REF] Meredith | Pitch spelling algorithms[END_REF]: such methods are almost perfect and provide the diatonic spelling of some pitches. However, we also tested a pseudo-diatonic matching on semitone information -considering as similar the intervals that differ from at most 1 semitone. The results (not shown) are very similar to those with true diatonic matching.

Other edit operations. Finally, we also tested other edit operations. The equations of Fig. 5 consider only substitutions, and can be simply extended to include the full Mongeau-Sankoff edit operations [START_REF] Mongeau | Comparaison of musical sequences[END_REF]. For instance, using insertions and allowing rhythm substitutions will, starting from the true subject, retrieve the occurrences #2 and #3 in Fig. 6. However, in the general case, insertions or deletions destroy the measure, leading to bad results on the predicted subject lengths.

More musical operations (fragmentation, consolidation), with fine-tuned costs, give a slight advantage in some of the 24 fugues, but this has not been reported here to keep the simplicity of the algorithm.

Conclusions

A complete fugue analysis tool should use any available information, including pattern repetition, harmonic analysis and phrasing considerations.

In this work, we focused only on pattern repetition. Our simple algorithm, based on the total number of matched intervals in all occurrences of patterns, allows to find precise ends of subjects and first counter-subjects in the majority of cases. This model considers a unique substitution operation with a diatonic similarity, enforcing the equality of lengths for all notes except the first and the last ones.

Extensions could include a study on the second counter-subject and on other inferred patterns. Combined with other techniques, this algorithm could lead to a more complete automatic fugue analysis pipeline.

Track-separated data. The current algorithm works on track-separated data. Starting from plain MIDI files, we could use voice separating algorithms. Although it would be a challenging problem to adapt our algorithm to directly treat standard polyphonic MIDI files, we first want to improve the current approach to complete our comprehension of any fugue.

Studying other fugues. Finally, it would be interesting to study the efficiency of our algorithm on other fugues than Bach's Well-Tempered Clavier, keeping in mind some pratical limitations (availability of track-separated files, ground truth). As far as the fugues keep the strict structure with a clear subject exposition, we are confident that our algorithm should give good results. As an example, the website http://www.lifl.fr/ ~giraud/fugues shows the output on the fugue in Mozart's Adagio and Fugue in C minor, K 546. We plan to further experiment it on other baroque or classical fugues, or on more recent corpus such as the Shostakovitch preludes and fugues (op. 87). The correct end for the CS is detected, but the actual CS begins after a small codetta. #8: The correct end for the CS is detected, but the actual CS begins after a small codetta.

Moreover, two incomplete inverted subject (also noted [START_REF] Siglind Bruhn | s Well-Tempered Clavier[END_REF]) are detected on measure 54. #9: The values for CS ( ) are not counted in the total, as the CS is presented in a segmented form in almost all measures of the fugue [?].

Table 2. Results of the proposed algorithm on the 24 fugues of the first book of Bach's Well-Tempered Clavier. We take as a truth the analysis of [START_REF] Siglind Bruhn | s Well-Tempered Clavier[END_REF], keeping here only the complete occurrence of each pattern. The columns "occ" lists the number of occurrences of Subjects and Counter-Subjects. The values s and cs indicate the index of the note ending the true subject and the counter-subject, whereas s and cs are the values predicted by the algorithm. See Section 3 for a definition of g. All false positives (FP) are counted in the remarks.

Fig. 1 .

 1 Fig. 1. A monophonic sequence of notes (start of Fugue #2, see Fig. 2), represented by (p, o, ) or ( ∆ p, o, ) triplets. In this example, onsets and lengths are counted in sixteenths, and pitches and intervals are counted in semitones through the MIDI standard.

2

 2 Fig. 1. A monophonic sequence of notes (start of Fugue #2, see Fig. 2), represented by (p, o, ) or ( ∆ p, o, ) triplets. In this example, onsets and lengths are counted in sixteenths, and pitches and intervals are counted in semitones through the MIDI standard.

  Fig. 1. A monophonic sequence of notes (start of Fugue #2, see Fig. 2), represented by (p, o, ) or ( ∆ p, o, ) triplets. In this example, onsets and lengths are counted in sixteenths, and pitches and intervals are counted in semitones through the MIDI standard.

Fig. 2 .

 2 Fig. 2. Start of Fugue #2 in C minor (BWV 847).

Fig. 3 .

 3 Fig. 3. Analysis of Fugue #2 in C minor (BWV 847). Top: diagram summarizing the analysis by S. Bruhn, used with permission [2], [3, p. 80]. Bottom: output of the proposed algorithm, retrieving all occurrences of S (and their degrees) and all but one occurrences of CS. The numbers indicate the pitch intervals exactly matching (in a diatonic way) those of the patterns (out of 19 for S). The two S17 occurrences correspond thus to approximate matches of the subject (tonal answers).

Fig. 5 .

 5 Fig. 5. Substitution operations between intervals. The actual comparison of length ( = ) also checks the equality of the rests that may be immediately before the compared notes. The relation ≈ is a similarity relation on pitch intervals.

Table 1 .

 1 Occurrences and scores when matching all candidate subjects in Fugue #2.

	z + g 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
	g	-8 -7 -6 -5 -4 -3 -2 -1 0 +1 +2 +3 +4 +5 +6
	occ. 8	8	8	8	8	8	8	3	3	3	3	2	2	2	2
	score 100 108 116 124 132 140 148 59 61 63 65 48 50 52 54

A part of the output of the algorithm is shown at the bottom of Fig.3, and the full output on all the 24 fugues is available at http://www.lifl.fr/ ~giraud/fugues.