Numerical
Optimal Transport
http://optimaltransport.github.io

Gradient Flows
Gabriel Peyré
www.numerical-tours.com
ÉCOLE NORMALES U P É R I E U R E

Overview• Gradient Flows• Lagrangian Discretization• Eulerian Discretization and EntropicRegularization• Unbalanced OT• Generalized Sinkhorn

Implicit Euler Stepping

Metric space (X , d), minimize F (x) on X .

p=1

Implicit Euler step:def.

xk+1 = argmin d(xk , x)2 + ⌧ F (x)x2X

12=+xk+2p

p=xk
F (x) = ||x||2 on (X = R2 , || · ||p )xk+1{x ; d(xk , x) ⇠ ⌧ }xk+3

Implicit vs. Explicit Stepping

Metric space (X , d), minimize F (x) on X .

Explicit

Implicit

xk+1 = argmin d(xk , x)2 + ⌧ hrF (xk ), xi

xk+1 = argmin d(xk , x)2 + ⌧ F (x)

x2X

x2X

p=1

p

F (x) = ||x||2 on (X = R2 , || · ||p )

1+=

2
=
p

p=+1

2
p=

p=1

Wasserstein Gradient Flows
Implicit Euler step:[Jordan, Kinderlehrer, Otto 1998]
def.
W2µt+1 = Prox⌧ f (µt ) = argmin W2 (µt , µ) + ⌧ f (µ)µ2M+ (X)

Formal limit ⌧ ! 0: @t µ = div (µr(f 0 (µ)))
R(heat di↵usion)f (µ) = log( dµ@t µ = µdx )dµR
f (µ) = wdµ(advection)@t µ = div(µrw)
R dµ m 1m (non-linear di↵usion)
1@t µ = µf (µ) = m 1 ( dx )dµrw

Evolution µtrw

Evolution µt

Euclidean L2 flow
@f= E 0 (f )@t
RE(f ) = ||rf ||2

Optimal transport flow
@f= div(f r(E 0 (f )))@t
RE(f ) =f (log(f ) 1)

min E(f )fp=1

Heat equation:
p=2

E(f ) =R@f=@tf

fp 1 pf p 1

Richard Jordan

David Kinderlehrer

p=4@fPorous medium:=@tfp
Felix OttoLinear vs Non-linear Diffusions

Overview• Gradient Flows• Lagrangian Discretization• Eulerian Discretization and EntropicRegularization• Unbalanced OT• Generalized Sinkhorn

Eulerian vs. Lagrangian Discretization

Lagrangian Discretization of Entropyµ̂n = n1Pn

i=1 xin ! +1µµ̂n
def. PĤ(µ̂n ) = i log(minj6=i ||xidef.

xj ||)

n ! +1

RH(µ) =log( dµdx (x))dµ(x)Lagrangian Discretization of Gradient Flowsdef.

min E(⇢) =⇢ZV (x)⇢(x)dx +⇢(x) log(⇢(x))dx⇢t + r(V ⇢t )V (x) = ||x||2d⇢t=
Wasserstein flow of E:
dtZ

time

Overview• Gradient Flows• Lagrangian Discretization• Eulerian Discretization and EntropicRegularization• Unbalanced OT• Generalized Sinkhorn

Gradient Flows: Crowd Motion
def.

µt+1 = argminµ W↵↵ (µt , µ) + ⌧ f (µ)

Congestion-inducing function:f (µ) = ◆[0,] (µ) + hw, µi[Maury, Roudne↵-Chupin, Santambrogio 2010]

Proposition: Prox 1" f (µ) = min(e "w µ, ) = ||µt=0 ||1 = 2||µt=0 ||1rw

 = 4||µt=0 ||1Crowd Motion on a Surface
X = triangulated mesh.

 = ||µt=0 ||1 = 6||µt=0 ||1Potential cos(w)Gradient Flows: Crowd Motion with Obstacles
X = sub-domain of R2 .

 = ||µt=0 ||1 = 2||µt=0 ||1  = 4||µt=0 ||1  = 6||µt=0 ||1 Potential cos(w)Crowd of Sheeps

https://www.youtube.com/watch?v=tDQw21ntR64
Tim Whittaker (New Zealand)(↵t+⌧ ,22t+⌧ ) = argmin W2 (↵t , ↵) + W2 ( t ,(↵, )6C↵t↵t+⌧...
ttt+⌧t+⌧) + ⌧ W22 (↵, )

↵t+2⌧

...

t+2⌧

tt + 2⌧

Overview• Gradient Flows• Lagrangian Discretization• Eulerian Discretization and EntropicRegularization• Unbalanced OT• Generalized Sinkhorn

Unbalanced Transport⇣ ⌘RRdef.
d⇠2(⇠, µ) 2 M+ (X) , KL(⇠|µ) = X log dµ dµ + X (dµd⇠)def.

W Fc (µ, ⌫) = min hc, ⇡i + KL(P1] ⇡|µ) + KL(P2] ⇡|⌫)⇡[Liereo, Mielke, Savaré 2015]⇡Ifc(x,y)=log(cos(min(d(x,y),))Proposition:2
1/2
then W Fc is a distance on M+ (X).

[Liereo, Mielke, Savaré 2015] [Chizat, Schmitzer, Peyré, Vialard 2015]! “Dynamic” Benamou-Brenier formulation.

[Liereo, Mielke, Savaré 2015] [Kondratyev, Monsaingeon, Vorotnikov, 2015][Chizat, Schmitzer, P, Vialard 2015]

Balanced OT

Unbalanced OTOverview• Gradient Flows• Lagrangian Discretization• Eulerian Discretization and EntropicRegularization• Unbalanced OT• Generalized SinkhornGeneralized Entropic Regularization
Primal: min hdp , ⇡i + f1 (P1] ⇡) + f2 (P2] ⇡) + "KL(⇡|⇡0 )⇡Dual:

maxu,vf1⇤ (u)f2⇤ (v)⇡(x, y) = a(x)K(x, y)b(y)Block coordinates maxu
relaxation:
maxv

u""hedef.

, Ke(a, b) = (ef1⇤ (u)f2⇤ (v)"he"heu
"v
"v
"u"

v
",e, Ke⇤i,K ev")iu"i(Iu )
(Iv )Proposition: the solutions of (Iu ) and (Iv ) read:KL⇤
Prox(Ka)
ProxKL(Kb)
f2 /"
f1 /"
b=a=K ⇤aKbdef.
KLProxf1 /" (µ) = argmin ⌫ f1 (⌫) + "KL(⌫|µ)! Only matrix-vector multiplications. ! Highly parallelizable.
! On regular grids: only convolutions! Linear time iterations.

Generalized Sinkhorn

Generalized Sinkhorn

regularized OT problem (4.2) is a special
case ofa structured ofconvex
optimization
EntropicRegularization
OptimalTranspor
lem of the form82

ÿ

Entropic Regularization of Optimal Transport

shown in [Peyré,
et al.,+2015,
et Tal.,anmin 2015,
Ci,jFrognerPi,j ≠ ÁH(P)F (P1Chizat1n ).2017, Karlsson(4.49m ) + G(PPi,j 2015, Frogner
Asshownin[Peyré,al., hence
2015, Chizat
et al., 2017,Karlssonand an2016], Sinkhorn iterations(4.15)etcan
be extended
to thisproblem,Ringh, 2016], Sinkhorn iterations (4.15) can hence be extended to this problem, anded, defining F = ÿ{a} andG = ÿ{b} , where the indicatorof a closed convex
TKL function
readProxKL(Kv)
Prox(Ku)
F KL
KLG T
C isProx(Kv)Proxu)uΩandv
Ω,(4.51FG (KITu Ω Kv
and v Ω,(4.51)
0 if x œ C, KTKu u
KvÿCfor(x)=KL(4.50the proximal
operatorforthedivergenceiswherethe proximal
operator
theKL +Œdivergenceisotherwise,N NKLÕKLÕ ÕÕ’ uRœR,Prox(u)=argminKL(u|u)+F(u).
(4.52) (4.52’u œ,Prox(u)=argminKL(u|u)+F(u).
+F+FNN
retrieves the hard marginal constraints
defining
U(a, b). The proof of ProposiuuÕÕœRœR++4.3 somecarries
to this F,
moreproblem(4.49),
so thatsolutionof (4.49Forfunctions
G itgeneralis possible
to provethe linear
ratetheof uniqueconvergence
for iterme functions F, G it is possible to prove the linear rate of convergence for iter
has theformand(4.12).
ations(4.51),
these schemes can be generalized to arbitrary measures, see [Chizat(4.51),
andtheseschemes can be generalized to arbitrary measures, see [Chizaet al., 2017] for more details.
KL2017]for moredetails.
Iterations(4.51)
are thus interesting in the cases where ProxKL
F and ProxG can beKLrations(4.51)
are formthusorinterestingin thecasesProxandcomputedin closedvery efficiently.
Thisis in whereparticular
theKLcasefor Prox
separable
FG can bqfunctionsof the form
form For(u)= efficiently.
in thiscaseparticular the case for separabl
ted in closed
veryThisisini Fi (ui ) since
12q KLKL
ons of the form F (u) = Prox
in Fthis= Prox(uicase)i ) sincei FiF(u(u)
i