Course on Optimal Transport

Mathématiques - Informatique
Bruno Lévy
ALICE Géométrie & Lumière
CENTRE INRIA Nancy Grand-Est
Symposium on Geometry Processing – Paris – 2018
Course on Numerical Optimal Transport – Bruno Lévy

Part. 1. Goals and Motivations
Part. 2. Introduction to Optimal Transport
Part. 3. Semi-Discrete Optimal Transport
Part. 4. Applications in Computational Physics
OVERVIEW

Part. 1 Optimal Transport
Goal #1: “Understanding”

What I can’t create
I don’t understand
Richard Feynman

Jos Stam,
Stable Fluids, 1999
The art of fluid sim.
Understand fluids
Explain
Program

I wrote code

Your mission statement:
1. Understand the stuff
2. Explain it in simple terms
Be a good teacher, to others and to yourself
Know what you know and what you don’t know
Try to know what you don’t know
3. Program it

Measuring distances between functions

Measuring distances between function

Interpolating functions

Gaspard Monge - 1784

Gaspard Monge – geometry and light
ANR TOMMI Workshop

Cédric Villani
Optimal Transport Old & New
Topics on Optimal Transport
Yann Brenier
The polar factorization theorem
(Brenier Transport)
Monge-Brenier-Villani, the french connection

[Castro, Merigot, Thibert 2014]
Optimal transport
geometry and light
[Caffarelli, Kochengin, and Oliker 1999]

Part. 1 Optimal Transport – Image Processing
Barycenters / mixing textures
Video-style transfer,
A.I., “data sciences”
[Nicolas Bonneel, Julien Rabin, Gabriel
Peyŕe, Hanspeter Pfister]
[Nicolas Bonneel, Kalyan Sunkavalli, Sylvain
Paris, Hanspeter Pfister]
[Marco Cuturi, Gabriel Peyré]

Optimal transport - geometry and light
[Chwartzburg, Testuz, Tagliasacchi, Pauly, SIGGRAPH 2014]

Part. 1. Motivations
Discretization of functionals involving the Monge-Ampère operator,
Benamou, Carlier, Mérigot, Oudet
arXiv:1408.4536
The variational formulation of the Fokker-Planck equation
Jordan, Kinderlehrer and Otto
SIAM J. on Mathematical Analysis
Geostrophic current

How to “morph” a shape into another one of same mass
while minimizing the “effort” ?
?
T

How to “morph” a shape into another one of same mass
while minimizing the “effort” ?
?
T
The “effort” of the best T defines a distance between the shapes

How to “morph” a shape into another one
while preserving mass and minimizing the effort ?
?
T

?
T
“minimum action principle”

?
T
“minimum action principle”“conservation law”

“minimum action principle subject to conservation law”
OT=
Yann Brenier:
“Each time the Laplace operator is used in a PDE,
it can be replaced with the Monge-Ampère operator”

OT=
Yann Brenier:
New ways of simulating physics with a computer

OT=
Yann Brenier:
Fast Fourier Transform

OT=
Yann Brenier:
Fast Fourier Transform Fast OT algo. ???

Optimal Transport
an elementary introduction
2

Part. 2 Optimal Transport – Monge’s problem
(X;μ) (Y;ν)
Two measures μ, ν such that ∫dμ(x) = ∫dν(x)
X Y

A map T is a transport map between μ and ν if
μ(T-1(B)) = ν(B) for any Borel subset B of Y
(X;μ) (Y;ν)
(Borel subset = subset that can be measured)

B
(X;μ) (Y;ν)

B
T-1(B)
(X;μ) (Y;ν)

Notation: if T is a transport map between μ and ν
then one writes ν = T#μ (ν is the pushforward of μ)
(X;μ) (Y;ν)

Monge’s problem:
Find a transport map T that minimizes C(T) = ∫X || x – T(x) ||2 dμ(x)
(X;μ) (Y;ν)

Monge’s problem:
• Difficult to study
• If μ has an atom (isolated Dirac),
it can only be mapped to another Dirac
(T needs to be a map)

Monge’s problem:
Transport from a measure concentrated on L1 onto another one concentrated on L2 and L3

Monge’s problem:
Transport from a measure concentrated on L1 onto another one concentrated on L2 and L3
The infimum is never realized by a map, need for a relaxation

Part. 2 Optimal Transport – Kantorovich
Monge’s problem:
Kantorovich’s problem (1942):
Find a measure γ defined on X x Y
such that ∫x in X dγ(x,y) = dν(y)
and ∫y in Y dγ(x,y) = dμ(x)
that minimizes ∫∫X x Y || x – y ||2 dγ(x,y)

Monge’s problem:
Kantorovich’s problem:
“γ(x,y)” :
How much sand goes from x to y

Monge’s problem:
Everything that is
transported from x sums to “μ(x)”

Monge’s problem:
Everything that is
transported to y sums to “ν(y)”

Transport plan – example 1/4 : translation of a segment

Transport plan – example 2/4 : spitting a segment

Transport plan – example 3/4 : splitting a Dirac into two Diracs

Transport plan – example 3/4 : splitting a Dirac into two Diracs
(No transport map)

Transport plan – example 4/4 : splitting a Dirac into two segments

Transport plan – example 4/4 : splitting a Dirac into two segments
(No transport map)

Part. 2 Optimal Transport – Duality

Duality is easier to understand with a discrete version
Then we’ll go back to the continuous setting.

(DMK):
Min <c, γ>
P1 γ = u
s.t. P2 γ = v
γ ≥ 0

(DMK):
Min <c, γ>
P1 γ = u
s.t. P2 γ = v
γ ≥ 0
γ =
γ11
γ12
…
γ1n
γ22
…
γ2n
…
…
γmn
Є IRmn

(DMK):
Min <c, γ>
P1 γ = u
s.t. P2 γ = v
γ ≥ 0
γ =
γ11
γ12
…
γ1n
γ22
…
γ2n
…
…
γmn
c =
c11
c12
…
c1n
c22
…
c2n
…
…
cmn
Є IRmnЄ IRmn

(DMK):
Min <c, γ>
P1 γ = u
s.t. P2 γ = v
γ ≥ 0
γ =
γ11
γ12
…
γ1n
γ22
…
γ2n
…
…
γmn
c =
c11
c12
…
c1n
c22
…
c2n
…
…
cmn
cij = || xi – yj ||2
Є IRmnЄ IRmn

(DMK):
Min <c, γ>
P1 γ = u
s.t. P2 γ = v
γ ≥ 0
γ =
γ11
γ12
…
γ1n
γ22
…
γ2n
…
…
γmn
c =
c11
c12
…
c1n
c22
…
c2n
…
…
cmn
cij = || xi – yj ||2
mn X m
Є IRmnЄ IRmn

(DMK):
Min <c, γ>
P1 γ = u
s.t. P2 γ = v
γ ≥ 0
γ =
γ11
γ12
…
γ1n
γ22
…
γ2n
…
…
γmn
c =
c11
c12
…
c1n
c22
…
c2n
…
…
cmn
cij = || xi – yj ||2
mn X n
mn X m
Є IRmnЄ IRmn

(DMK):
Min <c, γ>
P1 γ = u
s.t. P2 γ = v
γ ≥ 0
γ =
γ11
γ12
…
γ1n
γ22
…
γ2n
…
…
γmn
c =
c11
c12
…
c1n
c22
…
c2n
…
…
cmn
Є IRmnЄ IRmn
mn X n
mn X m

Consider L(φ,ψ) = <c, γ> - <φ, P1 γ - u> - <ψ, P2 γ - v>
(DMK):
Min <c, γ>
P1 γ = u
s.t. P2 γ = v
γ ≥ 0
< u, v > denotes the dot product between u and v

Remark: Sup[ L(φ,ψ) ] = < c, γ> if P1 γ = u and P2 γ = v
φ Є IRm
ψ Є IRn
(DMK):
Min <c, γ>
P1 γ = u
s.t. P2 γ = v
γ ≥ 0

φ Є IRm
ψ Є IRn
= +∞ otherwise
(DMK):
Min <c, γ>
P1 γ = u
s.t. P2 γ = v
γ ≥ 0

φ Є IRm
ψ Є IRn
= +∞ otherwise
Consider now: Inf [ Sup[ L(φ,ψ) ] ]
φ Є IRm
ψ Є IRn
γ ≥ 0
(DMK):
Min <c, γ>
P1 γ = u
s.t. P2 γ = v
γ ≥ 0

φ Є IRm
ψ Є IRn
= +∞ otherwise
Consider now: Inf [ Sup[ L(φ,ψ) ] ] = Inf [ < c, γ> ]
φ Є IRm
ψ Є IRn
γ ≥ 0 γ ≥ 0
P1 γ = u
P2 γ = v
(DMK):
Min <c, γ>
P1 γ = u
s.t. P2 γ = v
γ ≥ 0

φ Є IRm
ψ Є IRn
= +∞ otherwise
Consider now: Inf [ Sup[ L(φ,ψ) ] ] = Inf [ < c, γ> ] (DMK)
φ Є IRm
ψ Є IRn
γ ≥ 0 γ ≥ 0
P1 γ = u
P2 γ = v
(DMK):
Min <c, γ>
P1 γ = u
s.t. P2 γ = v
γ ≥ 0

φ Є IRm
ψ Є IRn
γ ≥ 0
Inf [ Sup[ <c, γ> - <φ, P1 γ - u> - <ψ, P2 γ - v> ] ]
(DMK):
Min <c, γ>
P1 γ = u
s.t. P2 γ = v
γ ≥ 0

φ Є IRm
ψ Є IRn
γ ≥ 0
φ Є IRm
ψ Є IRn
γ ≥ 0
Sup[ Inf[ <c, γ> - <φ, P1 γ - u> - <ψ, P2 γ - v> ] ]
Exchange Inf Sup
(DMK):
Min <c, γ>
P1 γ = u
s.t. P2 γ = v
γ ≥ 0

φ Є IRm
ψ Є IRn
γ ≥ 0
φ Є IRm
ψ Є IRn
γ ≥ 0
φ Є IRm
ψ Є IRn
γ ≥ 0
Sup[ Inf[ < γ,c-P1
t φ – P2
t ψ > + <φ,u> + <ψ, v> ] ]
Exchange Inf Sup
Expand/Reorder/Collect
(DMK):
Min <c, γ>
P1 γ = u
s.t. P2 γ = v
γ ≥ 0

φ Є IRm
ψ Є IRn
γ ≥ 0
φ Є IRm
ψ Є IRn
γ ≥ 0
φ Є IRm
ψ Є IRn
γ ≥ 0
Sup[ Inf[ < γ,c-P1
t φ – P2
t ψ > + <φ,u> + <ψ, v> ] ]
Exchange Inf Sup
Expand/Reorder/Collect
Interpret
(DMK):
Min <c, γ>
P1 γ = u
s.t. P2 γ = v
γ ≥ 0

φ Є IRm
ψ Є IRn
γ ≥ 0
Sup[ Inf[ < γ,c-P1
t φ – P2
t ψ > + <φ,u> + <ψ, v> ] ]
Interpret
Sup[ <φ,u> + <ψ, v> ]
φ Є IRm
ψ Є IRn
P1
t φ + P2
t ψ ≤ c
(DMK):
Min <c, γ>
P1 γ = u
s.t. P2 γ = v
γ ≥ 0
(DDMK)

φ Є IRm
ψ Є IRn
γ ≥ 0
Sup[ Inf[ < γ,c-P1
t φ – P2
t ψ > + <φ,u> + <ψ, v> ] ]
Interpret
Sup[ <φ,u> + <ψ, v> ]
φ Є IRm
ψ Є IRn
P1
t φ + P2
t ψ ≤ c
φi + ψj ≤ cij (i,j)
A
(DMK):
Min <c, γ>
P1 γ = u
s.t. P2 γ = v
γ ≥ 0
(DDMK)

Part. 2 Optimal Transport – Kantorovich dual
such that ∫x in X dγ(x,y) = dμ(x)
and ∫y in Y dγ(x,y) = dν(x)
Dual formulation of Kantorovich’s problem (Continuous):
Find two functions φ in L1(μ) and ψ in L1(ν)
Such that for all x,y, φ(x) + ψ(y) ≤ ½|| x – y ||2
that maximize ∫X φdμ + ∫Y ψdν

Dual formulation of Kantorovich’s problem:
Your point of view:
Try to minimize transport cost

Point of view of a “transport company”:
Try to maximize transport price
Your point of view:

that maximize ∫X φ(x)dμ + ∫Y ψ(y)dν
What they charge for loading at x
Your point of view:

What they charge for loading at x What they charge for unloading at y
Your point of view:

What they charge for loading at x What they charge for unloading at y
Price (loading + unloading) cannot
be greater than transport cost
(else you do the job yourself)
Your point of view:

Part. 2 Optimal Transport – c-conjugate functions

If we got two functions φ and ψ that satisfy the constraint
Then it is possible to obtain a better solution by replacing ψ with φc defined by:
For all y, φc(y) = inf x in X ½|| x – y ||2 - φ(y)

If we got two functions φ and ψ that satisfy the constraint
Then it is possible to obtain a better solution by replacing ψ with φc defined by:
For all y, φc(y) = inf x in X ½|| x – y ||2 - φ(y)
• φc is called the c-conjugate function of φ
• If there is a function φ such that ψ = φc then ψ is said to be c-concave
• If ψ is c-concave, then ψ cc = ψ

Find a c-concave function ψ
that maximizes ∫X ψ(x)dμ + ∫Y ψc(y)dν

ψ is called a “Kantorovich potential”

Part. 2 Optimal Transport – c-subdifferential
What about our initial problem ?

What about our initial problem ? (i.e., this is T() that we want to find …)

Proof: see OTON, chap. 10.

Heuristic argument (at the beginning of the same chapter):

w

-w

{
grad ψ(x) with ψ(x) := (½ x2- ψ(x))

Part. 2 Optimal Transport – convexity
{
ψ is convex

Part. 2 Optimal Transport – no collision
If T(.) exists, then
T(x) = x – grad ψ(x) = grad (½ x2- ψ(x) )
{grad ψ (x)
ψ is convex
Two transported particles cannot “collide”

Part. 2 Optimal Transport – Monge-Ampere
What about our initial problem ? If T(.) exists, then one can show that:
{
for all borel set A, ∫A dμ = ∫T(A) (|JT|) dν (change of variable)
Jacobian of T (1st order derivatives)

What about our initial problem ? If T(.) exists, then one can show that:
{
for all borel set A, ∫A dμ = ∫T(A) (|JT|) dν = ∫T(A) (H ψ ) dν
Det. of the Hessian of ψ (2nd order derivatives)

What about our initial problem ?
{
When μ and ν have a density u and v, (H ψ(x)). v(grad ψ(x)) = u(x)
Monge-Ampère
equation
for all borel set A, ∫A dμ = ∫T(A) (|JT|) dν = ∫T(A) (H ψ ) dν
{grad ψ(x) with ψ(x) := (½ x2- ψ(x))

Part. 2 Optimal Transport – summary

After several rewrites and under some conditions….
(Kantorovich formulation, dual, c-convex functions)

Monge-Ampère equation
(When μ and ν have a density u and v resp.)
Solve (H ψ(x)). v(grad ψ(x)) = u(x)

Brenier, Mc Cann, Trudinger: The optimal transport map is then given by:
T(x) = grad ψ(x)
Solve (H ψ(x)). v(grad ψ(x)) = u(x) Monge-Ampère equation
(When μ and ν have a density u and v resp.)

Semi-Discrete Optimal Transport
3

Part. 3 Optimal Transport – how to program ?
Continuous
(X;μ) (Y;ν)

Continuous
Semi-discrete
(X;μ) (Y;ν)

Continuous
Semi-discrete
Discrete
(X;μ) (Y;ν)

Part. 3 Optimal Transport – semi-discrete
∫X ψc (x)dμ + ∫Y ψ(y)dν
Sup
ψ Є ψc
(DMK)
(X;μ) (Y;ν)

(X;μ) (Y;ν)
∑j ψ(yj) vj
Sup
ψ Є ψc
(DMK)

Sup
ψ Є ψc
(DMK)
∑j ψ(yj) vj

Sup
ψ Є ψc
(DMK)
∫X inf yj ЄY [ || x – yj ||2 - ψ(yj) ] dμ
∑j ψ(yj) vj

Sup
ψ Є ψc
(DMK)
∑j ∫Lagψ(yj) || x – yj ||2 - ψ(yj) dμ
∫X inf yj ЄY [ || x – yj ||2 - ψ(yj) ] dμ
∑j ψ(yj) vj

G(ψ) =∑j ∫Lag ψ(yj) || x – yj ||2 - ψ(yj) dμ + ∑j ψ(yj) vj
Sup
ψ Є ψc
(DMK)
Where: Lag ψ(yj) = { x | || x – yj ||2 – ψ(yj) < || x – yj ||2 - ψ(yj’) } for all j’ ≠ j

Sup
ψ Є ψc
(DMK)
Laguerre diagram of the yj’s
(with the L2 cost || x – y ||2 used here, Power diagram)

Sup
ψ Є ψc
(DMK)
Weight of yj in the power diagram

Sup
ψ Є ψc
(DMK)
Weight of yj in the power diagram
ψ is determined by the
weight vector [ψ(y1) ψ(y2) … ψ(ym)]

Voronoi diagram: Vor(xi) = { x | d2(x,xi) < d2(x,xj) }
Part. 3 Power Diagrams

Power diagram: Pow(xi) = { x | d2(x,xi) – ψi < d2(x,xj) – ψj }
Voronoi diagram: Vor(xi) = { x | d2(x,xi) < d2(x,xj) }
Part. 3 Power Diagrams

Theorem: (direct consequence of MK duality
alternative proof in [Aurenhammer, Hoffmann, Aronov 98] ):
Given a measure μ with density, a set of points (yj), a set of positive coefficients vj
such that ∑ vj = ∫ dμ(x), it is possible to find the weights W = [ψ(y1) ψ(y2) … ψ(ym)]
such that the map TS
W is the unique optimal transport map
between μ and ν = ∑ vj δ(yj)

Theorem: (direct consequence of MK duality
alternative proof in [Aurenhammer, Hoffmann, Aronov 98] ):
Given a measure μ with density, a set of points (yj), a set of positive coefficients vj
such that ∑ vj = ∫ dμ(x), it is possible to find the weights W = [ψ(y1) ψ(y2) … ψ(ym)]
such that the map TS
W is the unique optimal transport map
between μ and ν = ∑ vj δ(yj)
Proof: G(ψ) =∑j ∫Lag ψ(yj) || x – yj ||2 - ψ(yj) dμ + ∑j ψ(yj) vj
Is a concave function of the weight vector [ψ(y1) ψ(y2) … ψ(ym)]

Part. 3 Optimal Transport – the AHA paper
Idea of the proof
Consider the function fT(W) = ∫(|| x– T(x) ||2 – ψ(T(X)))dμ(x)
The (unknown) weights W = [ψ(y1) ψ(y2) … ψ(ym)]

Idea of the proof
T : an arbitrary but fixed assignment.

Idea of the proof
W
fT(W)
fT(W)

Idea of the proof
fT(W) is linear in W
fT(W)
W
fT(W)
Fixed T

Idea of the proof
f TW
(W) : defined by Laguerre diagram
fT(W)
fixed W
W

Idea of the proof
f TW
(W) = minT fT(W)
fT(W)
fixed W
W

Idea of the proof
f: W → f TW
(W) is concave !!
(because its graph is the lower
enveloppe of linear functions)
fT(W)
W
fT(W)
f TW
(W) = minT fT(W)

Idea of the proof
f: W → f TW
(W) is concave
fT(W)
W
fT(W)
f TW
(W) = minT fT(W)
Consider g(W) = f TW
(W) + ∑vj ψj

fT(W)
W
fT(W)
f TW
(W) = minT fT(W)
Consider g(W) = f TW
(W) + ∑vj ψj
vj -∫Lag
ψ
(yj) dμ(x) and g is concave.∂ g / ∂ ψ j =
Idea of the proof
f: W → f TW
(W) is concave

Part. 3 Optimal Transport – the algorithm
Semi-discrete OT Summary:
Sup
ψ Є ψc
(DMK) G(ψ) =

G(ψ) = g(W) =∑j ∫Lag
ψ
(yj) || x – yj ||2 - ψ(yj) dμ + ∑j ψ(yj) vj is concave
Sup
ψ Є ψc
(DMK) G(ψ) =

ψ
vj -∫Lag(yj)dμ(x) (= 0 at the maximum)∂ G / ∂ ψ j =
Sup
ψ Є ψc
(DMK) G(ψ) =

ψ
vj -∫Lag(yj)dμ(x) (= 0 at the maximum)∂ G / ∂ ψ j =
Sup
ψ Є ψc
(DMK) G(ψ) =
Desired mass at yj Mass transported to yj

Part. 3 Optimal Transport – the Hessian
vj -∫Lag(yj)dμ(x)∂ G / ∂ ψ j =

∂
2
G / ∂ ψ i ψ j = - ∂ / ∂ ψ j ∫Lag(yj)dμ(x)

∂
2
yi
yj

∂
2
yi
yj
ψ j ψ j + δψ j

∂
2
yi
yj
ψ j ψ j + δψ j
Reynold’s thm:
∂ / ∂ ψ j ∫Lag(yj)dμ(x) = ∫∂ Lag(yj) v.n dμ(x)

yi
yj
Reynold’s thm:
fij(x) = 0

yi
yj
Reynold’s thm:
fij(x) = 0
c(x,yi) - c(x,yj) + ψ j - ψ i = 0

Part. 3 the Hessian
yi
yj
Reynold’s thm:
fij(x) = 0
c(x,yi) - c(x,yj) + ψ j - ψ i = 0
δx
δψ j
dfij = gradx(c(x,yi) - c(x,yj))dx + dψ j
n

Part. 3 the Hessian
yi
yj
Reynold’s thm:
fij(x) = 0
c(x,yi) - c(x,yj) + ψ j - ψ i = 0
δx
δψ j
n
δx = δh n = δh gradx fij(x) / || gradx fij(x)||

Part. 3 the Hessian
yi
yj
Reynold’s thm:
fij(x) = 0
c(x,yi) - c(x,yj) + ψ j - ψ i = 0
δx
δψ j
n
δh

Part. 3 the Hessian
yi
yj
Reynold’s thm:
fij(x) = 0
c(x,yi) - c(x,yj) + ψ j - ψ i = 0
δx
δψ j
n
δh
∂ h / ∂ ψ j = -1/ || gradx c(x,yi) – gradx c(x,yj) ||

Part. 3 the Hessian
yi
yj
Reynold’s thm:
fij(x) = 0
c(x,yi) - c(x,yj) + ψ j - ψ i = 0
δx
δψ j
n
δh
∂ h / ∂ ψ j = -1/ || gradx c(x,yi) – gradx c(x,yj) ||
∂ / ∂ ψ j ∫Lag(yj)dμ(x) = ∫Lag(yi) ∩Lag(yj) -1/|| gradx c(x,yi) – gradx c(x,yj) ||dμ(x)

Part. 3 the Hessian
∂
2
/ ∂ ψ i ∂ ψ j F = ∫Lag(yi) ∩ Lag(yj) -1/|| gradx c(x,yi) – gradx c(x,yj) ||dμ(x)
∂
2
/ ∂ ψ i
2 F = - ∑ ∂
2
/ ∂ ψ i ∂ ψ j

Part. 3 the Hessian
∂
2
∂
2
/ ∂ ψ i
2 F = - ∑ ∂
2
/ ∂ ψ i ∂ ψ j
c(x,y) = || x – y ||2
∂
2
/ ∂ ψ i ∂ ψ j F = ∫Lag(yi) ∩ Lag(yj) 1 / || xj – xi || dμ(x)

Part. 3 the Hessian
∂
2
∂
2
/ ∂ ψ i
2 F = - ∑ ∂
2
/ ∂ ψ i ∂ ψ j
c(x,y) = || x – y ||2
∂
2
/ ∂ ψ i ∂ ψ j F = ∫Lag(yi) ∩ Lag(yj) 1 / || xj – xi || dμ(x)
IP1 FEM Laplacian (not a big surprise)

Let’s program it !
Hierarchical algorithm [Mérigot]
Geometry, 3D [L], [L, Schwindt]
Damped Newton algorithm, [Kitagawa, Mérigot, Thibert]

Optimal Transport
applications in computational physics
4

Euler
Lagrange
The Least Action Principle
Hamilton,
Legendre,
Maupertuis
Axiom 1: There exists a function L(x,x,t) that describes the state
of a physical system
Short summary of the 1st chapter of Landau,Lifshitz Course of Theoretical Physics

Euler
Lagrange
Hamilton,
Legendre,
Maupertuis
position

Euler
Lagrange
Hamilton,
Legendre,
Maupertuis
position
speed

Euler
Lagrange
Hamilton,
Legendre,
Maupertuis
position
speed
time

Euler
Lagrange
∫t1
t2
L(x,x,t) dt
Hamilton,
Legendre,
Maupertuis
Axiom 2: The movement (time
evolution) of the physical system
minimizes the following integral

∫t1
t2
L(x,x,t) dt

∫t1
t2
L(x,x,t) dt
Theorem 1: (Lagrange equation):
∂L
∂x
d
dt
∂L
∂x
=
t=0
t=1

∫t1
t2
L(x,x,t) dt
Axiom 1: There exists L
Axiom 2: The movement minimizes
∂L
∂x
d
dt
∂L
∂x
=
Axiom 3:
Invariance w.r.t. change of
Gallileo frame + hom. + isotrop. :
x’
t’ =
x+vt
t

∫t1
t2
L(x,x,t) dt
∂L
∂x
d
dt
∂L
∂x
=
Axiom 3:
x’
t’ =
x+vt
t
Theorem 2:
∂L
∂x
x - L = cte

∫t1
t2
L(x,x,t) dt
∂L
∂x
d
dt
∂L
∂x
=
Axiom 3:
x’
t’ =
x+vt
t
Theorem 2:
∂L
∂x
x - L = cte
Homogeneity of time →
Preservation of energy

∫t1
t2
L(x,x,t) dt
∂L
∂x
d
dt
∂L
∂x
=
Axiom 3:
x’
t’ =
x+vt
t
Theorem 2:
∂L
∂x
x - L = cte
Homogeneity of space →
Preservation of momentum

∫t1
t2
L(x,x,t) dt
∂L
∂x
d
dt
∂L
∂x
=
Axiom 3:
x’
t’ =
x+vt
t
Theorem 2:
∂L
∂x
x - L = cte
Isotropy of space →
Preservation of angular momentum

∫t1
t2
L(x,x,t) dt
∂L
∂x
d
dt
∂L
∂x
=
Axiom 3:
x’
t’ =
x+vt
t
Theorem 2:
∂L
∂x
x - L = cte
Preserved quantities
“Integrals of Motion”
Noether’s theorem

Axiom 2: The movement minimizes ∫ L
∂L
∂x
d
dt
∂L
∂x
=
Axiom 3:
x’
t’ =
x+vt
t
Theorem 2:
∂L
∂x
x - L = cte
Free particle:
Theorem 3: v = cte (Newton’s law I)

∂L
∂x
d
dt
∂L
∂x
=
Axiom 3:
x’
t’ =
x+vt
t
Theorem 2:
∂L
∂x
x - L = cte
Free particle:
Expression of the Lagrangian:
L = ½ m v2

∂L
∂x
d
dt
∂L
∂x
=
Axiom 3:
x’
t’ =
x+vt
t
Theorem 2:
∂L
∂x
x - L = cte
Free particle:
L = ½ m v2
Expression of the Energy:
E = ½ m v2

∂L
∂x
d
dt
∂L
∂x
=
Axiom 3:
x’
t’ =
x+vt
t
Particle in a field:
L = ½ m v2 – U(x)
Free particle:
L = ½ m v2
E = ½ m v2

∂L
∂x
d
dt
∂L
∂x
=
Axiom 3:
x’
t’ =
x+vt
t
L = ½ m v2 – U(x)
E = ½ m v2 + U(x)
Free particle:
L = ½ m v2
E = ½ m v2

∂L
∂x
d
dt
∂L
∂x
=
Axiom 3:
x’
t’ =
x+vt
t
L = ½ m v2 – U(x)
E = ½ m v2 + U(x)
Theorem 4:
mx = - U (Newton’s law II)
Free particle:
L = ½ m v2
E = ½ m v2

∂L
∂x
d
dt
∂L
∂x
=
Axiom 3:
Invariance w.r.t. Lorentz change of
frame
x’
t’ =
(x-vt) x γ
(t – vx/c2) x γ
γ = 1 / √ ( 1 – v2 / c2)
(relativistic setting – just for fun…)

(relativistic setting – just for fun…)
∂L
∂x
d
dt
∂L
∂x
=
Axiom 3:
Invariance w.r.t. Lorentz change of
frame
x’
t’ =
(x-vt) x γ
(t – vx/c2) x γ
γ = 1 / √ ( 1 – v2 / c2)
Theorem 5:
E = ½ γ m v2 + mc2

(quantum physics setting – just for fun…)

?
T
Fluids – Benamou Brenier
ρ1 ρ2

?
T
Minimize
A(ρ,v) =
∫t1
t2
(t2-t1)
∫Ω
ρ(x,t) ||v(t,x)||2
dxdt
s.t. ρ(t1,.) = ρ1 ; ρ(t2,.) = ρ2 ; d ρ
dt
= - div(ρv)
ρ1 ρ2

?
T
∫t1
t2
(t2-t1)
∫Ω
ρ(x,t) ||v(t,x)||2
dxdt
s.t. ρ(t1,.) = ρ1 ; ρ(t2,.) = ρ2 ; d ρ
dt
= - div(ρv)
ρ1 ρ2
Minimize C(T) =
∫Ω
|| x – T(x) ||2 dx
s.t. T is measure-preserving
ρ1(x)
Minimize
A(ρ,v) =

Le schéma [Mérigot-Gallouet]
Applications en graphisme: [De Goes et.al] (power particles)
Part. 4 Optimal Transport – Fluids

Let’s code it !

Part. 4 To infinity and beyond…

René Descartes - 1663
Vortices in “ether” ?
Part. 4 To infinity and beyond

The Data: #1: the Cosmic Microwave Background
COBE 1992Part. 4 To infinity and beyond…

WMAP 2003
2006
2008
2010
The Data: #1 the Cosmic Microwave Background

The Data: #2 redshift acquisition surveys

The millenium simulation project, Max Planck Institute fur Astrophysik
pc/h : parsec (= 3.2 années lumières)Part. 4 To infinity and beyond…

The universal swimming pool

Coop. with MOKAPLAN & Institut d’Astrophysique de Paris & Observatoire de Paris
Time = Now
Early Universe Reconstruction

. . . . .

Time = BigBang
(- 13.7 billion Y)

“Time-warped” map of the
universe

Cosmic Microware Background:
“Fossil light” emitted 380 000 Y after BigBang
and measured now
universe

Cosmic Microware Background:
“Fossil light” emitted 380 000 Y after BigBang
and measured now
Do they match ?
universe

Conclusions
Open Questions
References
Online resources

Conclusions – Open questions
* Connections with physics, Legendre transform and entropy ?
[Cuturi & Peyré] – regularized discrete optimal transport – why does it work ?
Hint 1: Minimum action principle subject to conservation laws
Hint 2: Entropy = dual of temperature ; Legendre = Fourier[(+,*) → (Max,+)]…
* More continuous numerical algorithms ?
[Benamou & Brenier] fluid dynamics point of view – very elegant, but 4D problem !!
FEM-type adaptive discretization of the subdifferential (graph of T) ?
* Can we characterize OT in other semi-discrete settings ?
measures supported on unions of spheres
piecewise linear densities
* Connections with computational geometry ?
Singularity set [Figalli] = set of points where T is discontinuous
Looks like a “mutual power diagram”, anisotropic Voronoi diagrams

Conclusions - References
A Multiscale Approach to Optimal Transport,
Quentin Mérigot, Computer Graphics Forum, 2011
Variational Principles for Minkowski Type Problems, Discrete Optimal Transport,
and Discrete Monge-Ampere Equations
Xianfeng Gu, Feng Luo, Jian Sun, S.-T. Yau, ArXiv 2013
Minkowski-type theorems and least-squares clustering
AHA! (Aurenhammer, Hoffmann, and Aronov), SIAM J. on math. ana. 1998
Topics on Optimal Transportation, 2003
Optimal Transport Old and New, 2008
Cédric Villani

Conclusions - References
Polar factorization and monotone rearrangement of vector-valued functions
Yann Brenier, Comm. On Pure and Applied Mathematics, June 1991
A computational fluid mechanics solution of the Monge-Kantorovich mass transfer
problem, J.-D. Benamou, Y. Brenier, Numer. Math. 84 (2000), pp. 375-393
Pogorelov, Alexandrov – Gradient maps, Minkovsky problem (older than AHA
paper, some overlap, in slightly different context, formalism used by Gu & Yau)
Rockafeller – Convex optimization – Theorem to switch inf(sup()) – sup(inf())
with convex functions (used to justify Kantorovich duality)
Filippo Santambrogio – Optimal Transport for Applied Mathematician,
Calculus of Variations, PDEs and Modeling – Jan 15, 2015
Gabriel Peyré, Marco Cuturi, Computational Optimal Transport, 2018

Online resources
All the sourcecode/documentation available from:
http://alice.loria.fr/software/geogram
Demo: www.loria.fr/~levy/GLUP/vorpaview
* L., A numerical algorithm
for semi-discrete L2 OT in 3D,
ESAIM Math. Modeling
and Analysis, 2015
* L. and E. Schwindt,
Notions of OT and how to
implement them on a computer,
Computer and Graphics, 2018.
J. Lévy – 1936-2018 - To you Dad, I miss you so much.

Bonus Slides
The Isoperimetric Inequality

The isoperimetric inequality
For a given volume,
ball is the shape that minimizes border area

∫| grad f | ≥ n Vol(B2
n)1/n (∫f n/(n-1))(n-1)/n
L1 Sobolev inegality: Given f: IRn → IR sufficiently regular
Explanation in [Dario Cordero Erauquin] course notes

n)1/n (∫f n/(n-1))(n-1)/n
Consider a compact set Ω such that Vol(Ω) = Vol(B2
3)
and f = the indicatrix function of Ω

n)1/n (∫f n/(n-1))(n-1)/n
3)
Vol(∂Ω) ≥ n Vol(B2
3)1/3 Vol(B2
3)2/3

n)1/n (∫f n/(n-1))(n-1)/n
3)
Vol(∂Ω) ≥ n Vol(B2
3)1/3 Vol(B2
3)2/3
Vol(∂Ω) ≥ 4 π = Vol(∂B2
3)

n)1/n (∫f n/(n-1))(n-1)/n
L1 Sobolev inegality: a proof with OT [Gromov]

n)1/n (∫f n/(n-1))(n-1)/n
We suppose w.l.o.g. that ∫f n/(n-1) = 1

n)1/n (∫f n/(n-1))(n-1)/n
There exists an optimal transport T = gradΨ between
f n/(n-1)(x)dx and 1B2
n/Vol(B2
n)dx

n)1/n (∫f n/(n-1))(n-1)/n
n/Vol(B2
n)dx
Monge-Ampère equation: Vol(B2
n) fn/(n-1)(x) = det Hess Ψ

n)1/n (∫f n/(n-1))(n-1)/n
n/Vol(B2
n)dx
Arithmetico-geometric ineq: det (H) 1/n ≤ trace(H)/n if H positive

n)1/n (∫f n/(n-1))(n-1)/n
n/Vol(B2
n)dx
det (Hess Ψ) 1/n ≤ trace(Hess Ψ)/n

n)1/n (∫f n/(n-1))(n-1)/n
n/Vol(B2
n)dx
det (Hess Ψ) 1/n ≤ trace(Hess Ψ)/n
det (Hess Ψ) 1/n ≤ ∆ Ψ / n

n)1/n (∫f n/(n-1))(n-1)/n
det (Hess Ψ) 1/n ≤ (∆ Ψ)/n
Monge-Ampère equation:
Vol(B2

n)1/n (∫f n/(n-1))(n-1)/n
Vol(B2
n) = Vol(B2
n) ∫f n/(n-1)
= ∫f Vol(B2
n) f 1/(n-1)
Vol(B2

n)1/n (∫f n/(n-1))(n-1)/n
Vol(B2
n) = Vol(B2
n) ∫f n/(n-1)
= ∫f Vol(B2
n) f 1/(n-1)
≤ 1/n ∫ f ∆ Ψ
Vol(B2

n)1/n (∫f n/(n-1))(n-1)/n
Vol(B2
n) = Vol(B2
n) ∫f n/(n-1)
= ∫f Vol(B2
n) f 1/(n-1)
≤ 1/n ∫ f ∆ Ψ
∫f ∆ Ψ = - ∫ grad f . grad Ψ
Vol(B2

n)1/n (∫f n/(n-1))(n-1)/n
Vol(B2
n) = Vol(B2
n) ∫f n/(n-1)
= ∫f Vol(B2
n) f 1/(n-1)
≤ 1/n ∫ f ∆ Ψ
∫f ∆ Ψ = - ∫ grad f . grad Ψ ≤ ∫ | grad f | (T = grad Ψ Є B2
n )
Vol(B2

n)1/n (∫f n/(n-1))(n-1)/n
Vol(B2
n) = Vol(B2
n) ∫f n/(n-1)
= ∫f Vol(B2
n) f 1/(n-1)
≤ 1/n ∫ f ∆ Ψ
∫f ∆ Ψ = - ∫ grad f . grad Ψ ≤ ∫ | grad f | (T = grad Ψ Є B2
n )
∫ | grad f | ≥ n Vol(B2
n)1/n
■
Vol(B2

Bonus Slides
Plotting the potential & optics

The [AHA] paper summary:
• The optimal weights minimize a convex function
• The gradient and Hessian of this convex function is easy to compute
Note: the weight w(s) correspond to the Kantorovich potential ψ(x)
(solves a “discrete Monge-Ampere” equation)
The algorithm:
Summary:
The algorithm computes the weights wi such that the power cells associated with
the Diracs correspond to the preimages of the Diracs.
μ ν
Plotting the potential, “optics”

hi

Numerical Experiment: A disk becomes two disks

Course on Optimal Transport

Recommended

Recommended

More Related Content

What's hot

What's hot (20)

Similar to Course on Optimal Transport

Similar to Course on Optimal Transport (20)

More from Bruno Levy

More from Bruno Levy (12)

Recently uploaded

Recently uploaded (20)

Course on Optimal Transport