Reflector antenna problem Boris Thibert LJK Universit e de - - PowerPoint PPT Presentation
Reflector antenna problem Boris Thibert LJK Universit e de Grenoble Joint work with Quentin M erigot and Pedro Machado Modelisation with optimal transport October 3-4, 2013 1 Far-Field Reflector Antenna Problem Punctual light at origin
1
Modelisation with optimal transport October 3-4, 2013
LJK Universit´ e de Grenoble Joint work with Quentin M´ erigot and Pedro Machado
2
S2
∞
Punctual light at origin o, µ measure on S2
∞
S2
ν
2
S2
∞
Punctual light at origin o, µ measure on S2
∞
S2
(S∞, ν) under reflection by Snell’s law.
µ ν
2
S2
∞
Punctual light at origin o, µ measure on S2
∞
S2
(S∞, ν) under reflection by Snell’s law.
◮ R is parameterized over S2
ν
2
S2
∞
Punctual light at origin o, µ measure on S2
∞
S2
(S∞, ν) under reflection by Snell’s law.
◮ R is parameterized over S2
TR : x ∈ S2
0 → y = x − 2x|nn
µ ν
2
S2
∞
Punctual light at origin o, µ measure on S2
∞
S2
(S∞, ν) under reflection by Snell’s law.
◮ R is parameterized over S2
TR : x ∈ S2
0 → y = x − 2x|nn
Brenier formulation T♯µ = ν µ ν
2
S2
∞
Punctual light at origin o, µ measure on S2
∞
S2
(S∞, ν) under reflection by Snell’s law.
◮ R is parameterized over S2
TR : x ∈ S2
0 → y = x − 2x|nn
Brenier formulation i.e. for every borelian B T♯µ = ν µ(T −1(B)) = ν(B) µ ν
2
S2
∞
Punctual light at origin o, µ measure on S2
∞
S2
(S∞, ν) under reflection by Snell’s law.
◮ R is parameterized over S2
TR : x ∈ S2
0 → y = x − 2x|nn
Brenier formulation i.e. for every borelian B Monge-Ampere equation g(T(x)) det(DT(x)) = f(x) T♯µ = ν µ(T −1(B)) = ν(B) If µ(x) = f(x)dx and ν(y) = g(y)dy ◮ highly non linear µ ν
2
S2
∞
Punctual light at origin o, µ measure on S2
∞
S2
(S∞, ν) under reflection by Snell’s law.
◮ R is parameterized over S2
TR : x ∈ S2
0 → y = x − 2x|nn
Brenier formulation i.e. for every borelian B Monge-Ampere equation g(T(x)) det(DT(x)) = f(x)
Caffarelli & Oliker 94
T♯µ = ν µ(T −1(B)) = ν(B) If µ(x) = f(x)dx and ν(y) = g(y)dy ◮ highly non linear ◮ Existence
Wang 96, Guan & Wang 98
◮ Regularity, uniqueness µ ν
3
S2
∞
Punctual light at origin o, µ measure on S2
3
S2
∞
Punctual light at origin o, µ measure on S2
∞
S2
3
S2
∞
Punctual light at origin o, µ measure on S2
∞
S2
y1 R : paraboloid of direction y1 and focal O
3
S2
∞
y1 y2 y3 Punctual light at origin o, µ measure on S2
i νiδyi on S2 ∞
3
y1 y2 y3 Punctual light at origin o, µ measure on S2
direction yi and focal distance κi R( κ) = ∂
i=1Pi(κi)
P2 µ P1 Prescribed far-field: ν =
i νiδyi on S2 ∞
3
y1 y2 y3 Punctual light at origin o, µ measure on S2
direction yi and focal distance κi R( κ) = ∂
i=1Pi(κi)
κ) = πS2
κ) ∩ ∂Pi(κi))
R( κ) ∩ ∂P3(κ3) PI3( κ)
P2 µ = directions that are reflected towards yi. Prescribed far-field: ν =
i νiδyi on S2 ∞
3
y1 y2 y3 Punctual light at origin o, µ measure on S2
direction yi and focal distance κi R( κ) = ∂
i=1Pi(κi)
κ) = πS2
κ) ∩ ∂Pi(κi)) Problem (FF): Find κ1, . . . , κN such that for every i, µ(PIi( κ)) = νi.
R( κ) ∩ ∂P3(κ3) PI3( κ) amount of light reflected in direction yi.
P2 µ = directions that are reflected towards yi. Prescribed far-field: ν =
i νiδyi on S2 ∞
4
Lemma: With c(x, y) = − log(1 − x|y), and ψi := log(κi), PIi( κ) = {x ∈ S2
0, c(x, yi) + ψi ≤ c(x, yj) + ψj
∀j}. P1
PI3( κ)
P2
Caffarelli-Oliker ’94
4
Proof: ∂Pi(κi) is parameterized in radial coordinates by Lemma: With c(x, y) = − log(1 − x|y), and ψi := log(κi), PIi( κ) = {x ∈ S2
0, c(x, yi) + ψi ≤ c(x, yj) + ψj
∀j}. P1
PI3( κ)
P2 ρi : x ∈ S2
κi 1−x|yi
Caffarelli-Oliker ’94
4
Proof: ∂Pi(κi) is parameterized in radial coordinates by ⇐ ⇒ i ∈ arg minj
κj 1−x|yj
x ∈ PIi( κ) Lemma: With c(x, y) = − log(1 − x|y), and ψi := log(κi), PIi( κ) = {x ∈ S2
0, c(x, yi) + ψi ≤ c(x, yj) + ψj
∀j}. P1
PI3( κ)
P2 ρi : x ∈ S2
κi 1−x|yi
Caffarelli-Oliker ’94
4
Proof: ∂Pi(κi) is parameterized in radial coordinates by ⇐ ⇒ i ∈ arg minj
κj 1−x|yj
x ∈ PIi( κ) ⇐ ⇒ i ∈ arg minj log(κj) − log(1 − x|yj) Lemma: With c(x, y) = − log(1 − x|y), and ψi := log(κi), PIi( κ) = {x ∈ S2
0, c(x, yi) + ψi ≤ c(x, yj) + ψj
∀j}. P1
PI3( κ)
P2 ρi : x ∈ S2
κi 1−x|yi
Caffarelli-Oliker ’94
4
Proof: ∂Pi(κi) is parameterized in radial coordinates by ⇐ ⇒ i ∈ arg minj
κj 1−x|yj
x ∈ PIi( κ) ⇐ ⇒ i ∈ arg minj log(κj) − log(1 − x|yj) ⇐ ⇒ i ∈ arg minj ψj + c(x, yj) Lemma: With c(x, y) = − log(1 − x|y), and ψi := log(κi), PIi( κ) = {x ∈ S2
0, c(x, yi) + ψi ≤ c(x, yj) + ψj
∀j}. P1
PI3( κ)
P2 ρi : x ∈ S2
κi 1−x|yi
Caffarelli-Oliker ’94
4
Proof: ∂Pi(κi) is parameterized in radial coordinates by ⇐ ⇒ i ∈ arg minj
κj 1−x|yj
x ∈ PIi( κ) ⇐ ⇒ i ∈ arg minj log(κj) − log(1 − x|yj) ⇐ ⇒ i ∈ arg minj ψj + c(x, yj) Lemma: With c(x, y) = − log(1 − x|y), and ψi := log(κi), PIi( κ) = {x ∈ S2
0, c(x, yi) + ψi ≤ c(x, yj) + ψj
∀j}. P1
PI3( κ)
P2 ρi : x ∈ S2
κi 1−x|yi
Wang ’04 Caffarelli-Oliker ’94
◮ An optimal transport problem
5
µ = probability measure on X ν = prob. measure on finite Y y with density, X = manifold =
y∈Y νyδy
5
µ = probability measure on X ν = prob. measure on finite Y T −1(y) y with density, X = manifold =
y∈Y νyδy
Transport map: T : X → Y s.t. ∀y ∈ Y, µ(T −1({y})) = νy in short: T#µ = ν.
5
µ = probability measure on X ν = prob. measure on finite Y T −1(y) y with density, X = manifold =
y∈Y νyδy
Transport map: T : X → Y s.t. ∀y ∈ Y, µ(T −1({y})) = νy in short: T#µ = ν. Cc(T) =
Cost function: c : X × Y → R
5
µ = probability measure on X ν = prob. measure on finite Y T −1(y) y with density, X = manifold =
y∈Y νyδy
Transport map: T : X → Y s.t. ∀y ∈ Y, µ(T −1({y})) = νy in short: T#µ = ν. Cc(T) =
Cost function: c : X × Y → R =
y
5
µ = probability measure on X ν = prob. measure on finite Y Monge problem: Tc(µ, ν) := min{Cc(T); T#µ = ν} T −1(y) y with density, X = manifold =
y∈Y νyδy
Transport map: T : X → Y s.t. ∀y ∈ Y, µ(T −1({y})) = νy in short: T#µ = ν. Cc(T) =
Cost function: c : X × Y → R =
y
Aurenhammer, Hoffman, Aronov ’98 Merigot ’2010
6
y We assume (Twist), i.e. c ∈ C∞ and ∀x ∈ X the map y ∈ Y → ∇xc(x, y) is injective.
Y finite set, ψ : Y → R
6
y T ψ
c (x) = arg miny∈Y c(x, y) + ψ(y)
We assume (Twist), i.e. c ∈ C∞ and ∀x ∈ X the map y ∈ Y → ∇xc(x, y) is injective.
Y finite set, ψ : Y → R
6
y Vorψ
c (y) = {x ∈ Rd; T ψ c (x) = y}
T ψ
c (x) = arg miny∈Y c(x, y) + ψ(y)
We assume (Twist), i.e. c ∈ C∞ and ∀x ∈ X the map y ∈ Y → ∇xc(x, y) is injective. = generalized weighted Voronoi cell
Y finite set, ψ : Y → R
6
y Vorψ
c (y) = {x ∈ Rd; T ψ c (x) = y}
T ψ
c (x) = arg miny∈Y c(x, y) + ψ(y)
We assume (Twist), i.e. c ∈ C∞ and ∀x ∈ X the map y ∈ Y → ∇xc(x, y) is injective. NB: Under (Twist), (Vorψ
c (y))y∈Y partitions X and T ψ c well-defined a.e.
= generalized weighted Voronoi cell
Y finite set, ψ : Y → R
6
y Vorψ
c (y) = {x ∈ Rd; T ψ c (x) = y}
T ψ
c (x) = arg miny∈Y c(x, y) + ψ(y)
Lemma: Given a measure µ with density and ψ : Y → R, the map T ψ
c is a c-optimal transport between µ and T ψ c#µ.
We assume (Twist), i.e. c ∈ C∞ and ∀x ∈ X the map y ∈ Y → ∇xc(x, y) is injective. NB: Under (Twist), (Vorψ
c (y))y∈Y partitions X and T ψ c well-defined a.e.
= generalized weighted Voronoi cell
Y finite set, ψ : Y → R
6
y Vorψ
c (y) = {x ∈ Rd; T ψ c (x) = y}
T ψ
c (x) = arg miny∈Y c(x, y) + ψ(y)
Lemma: Given a measure µ with density and ψ : Y → R, the map T ψ
c is a c-optimal transport between µ and T ψ c#µ.
We assume (Twist), i.e. c ∈ C∞ and ∀x ∈ X the map y ∈ Y → ∇xc(x, y) is injective. NB: Under (Twist), (Vorψ
c (y))y∈Y partitions X and T ψ c well-defined a.e.
◮ Note: T ψ
c#µ = y∈Y µ(Vorψ c (y))δy.
= generalized weighted Voronoi cell
Y finite set, ψ : Y → R
6
y Vorψ
c (y) = {x ∈ Rd; T ψ c (x) = y}
T ψ
c (x) = arg miny∈Y c(x, y) + ψ(y)
Lemma: Given a measure µ with density and ψ : Y → R, the map T ψ
c is a c-optimal transport between µ and T ψ c#µ.
We assume (Twist), i.e. c ∈ C∞ and ∀x ∈ X the map y ∈ Y → ∇xc(x, y) is injective. NB: Under (Twist), (Vorψ
c (y))y∈Y partitions X and T ψ c well-defined a.e.
◮ Converse ? ◮ Note: T ψ
c#µ = y∈Y µ(Vorψ c (y))δy.
= generalized weighted Voronoi cell
Y finite set, ψ : Y → R
7
Lemma: With c(x, y) = − log(1 − x|y), and ψi := log(κi), PIi( κ) = {x ∈ S2
0, c(x, yi) + ψi ≤ c(x, yj) + ψj
∀j}. P1
PI3( κ)
P2 y1 y2 y3
7
Lemma: With c(x, y) = − log(1 − x|y), and ψi := log(κi), PIi( κ) = {x ∈ S2
0, c(x, yi) + ψi ≤ c(x, yj) + ψj
∀j}. P1
PI3( κ)
P2 Optimal transport formulation y1 y2 y3
◮ PIi( κ) = Vorψ
c (yi).
◮ T ψ
c (x) = arg miny∈Y c(x, y) + ψ(y)
7
Lemma: With c(x, y) = − log(1 − x|y), and ψi := log(κi), PIi( κ) = {x ∈ S2
0, c(x, yi) + ψi ≤ c(x, yj) + ψj
∀j}. P1
PI3( κ)
P2 Optimal transport formulation The map T ψ
c is a c-optimal transport between µ and T ψ c#µ.
y1 y2 y3
◮ PIi( κ) = Vorψ
c (yi).
◮ T ψ
c (x) = arg miny∈Y c(x, y) + ψ(y)
7
Lemma: With c(x, y) = − log(1 − x|y), and ψi := log(κi), PIi( κ) = {x ∈ S2
0, c(x, yi) + ψi ≤ c(x, yj) + ψj
∀j}. P1
PI3( κ)
P2 Optimal transport formulation The map T ψ
c is a c-optimal transport between µ and T ψ c#µ.
y1 y2 y3
◮ PIi( κ) = Vorψ
c (yi).
◮ T ψ
c (x) = arg miny∈Y c(x, y) + ψ(y)
Problem (FF): Find ψ1, . . . , ψN such that T ψ
c#µ = ν.
8
Cafarelli-Kochengin-Oliker’99: coordinate-wise ascent, with minimum increment
8
∀y ∈ Y \ {y0}, µ(Vorψ
c (p)) ≤ νy + δ
Initialization: Fix y0 ∈ Y , let δ = ε/N and compute ψ s.t. While ∃y = y0 such that µ(Vorψ
c (y)) ≤ νy − δ, do:
decrease ψ(y) s.t. µ(Vorψ
c (y)) ∈ [νy, νy + δ],
Result: ψ s.t. for all y, |µ(Vorψ
c (y)) − νy| ≤ ε.
Cafarelli-Kochengin-Oliker’99: coordinate-wise ascent, with minimum increment
8
∀y ∈ Y \ {y0}, µ(Vorψ
c (p)) ≤ νy + δ
Initialization: Fix y0 ∈ Y , let δ = ε/N and compute ψ s.t. While ∃y = y0 such that µ(Vorψ
c (y)) ≤ νy − δ, do:
decrease ψ(y) s.t. µ(Vorψ
c (y)) ∈ [νy, νy + δ],
Result: ψ s.t. for all y, |µ(Vorψ
c (y)) − νy| ≤ ε.
Cafarelli-Kochengin-Oliker’99: coordinate-wise ascent, with minimum increment
8
∀y ∈ Y \ {y0}, µ(Vorψ
c (p)) ≤ νy + δ
Initialization: Fix y0 ∈ Y , let δ = ε/N and compute ψ s.t. While ∃y = y0 such that µ(Vorψ
c (y)) ≤ νy − δ, do:
decrease ψ(y) s.t. µ(Vorψ
c (y)) ∈ [νy, νy + δ],
Result: ψ s.t. for all y, |µ(Vorψ
c (y)) − νy| ≤ ε.
◮ Complexity of SP: N 2/ε steps Cafarelli-Kochengin-Oliker’99: coordinate-wise ascent, with minimum increment
8
∀y ∈ Y \ {y0}, µ(Vorψ
c (p)) ≤ νy + δ
Initialization: Fix y0 ∈ Y , let δ = ε/N and compute ψ s.t. While ∃y = y0 such that µ(Vorψ
c (y)) ≤ νy − δ, do:
decrease ψ(y) s.t. µ(Vorψ
c (y)) ∈ [νy, νy + δ],
Result: ψ s.t. for all y, |µ(Vorψ
c (y)) − νy| ≤ ε.
◮ Complexity of SP: N 2/ε steps Cafarelli-Kochengin-Oliker’99: coordinate-wise ascent, with minimum increment ◮ Generalization of Oliker–Prussner in R2 with c(x, y) = x − y2
8
∀y ∈ Y \ {y0}, µ(Vorψ
c (p)) ≤ νy + δ
Initialization: Fix y0 ∈ Y , let δ = ε/N and compute ψ s.t. While ∃y = y0 such that µ(Vorψ
c (y)) ≤ νy − δ, do:
decrease ψ(y) s.t. µ(Vorψ
c (y)) ∈ [νy, νy + δ],
Result: ψ s.t. for all y, |µ(Vorψ
c (y)) − νy| ≤ ε.
◮ Complexity of SP: N 2/ε steps ◮ Generalization: MTW+ costs
Kitagawa ’12
Cafarelli-Kochengin-Oliker’99: coordinate-wise ascent, with minimum increment ◮ Generalization of Oliker–Prussner in R2 with c(x, y) = x − y2
9
Aurenhammer, Hoffman, Aronov ’98
Theorem: κ solves (FF) iff ψ = log( κ) maximizes Φ(ψ) :=
i
c (yi)[c(x, yi) + ψi] d µ(x) −
i ψiνi
with c(x, y) = − log(1 − x|y).
9
Aurenhammer, Hoffman, Aronov ’98
◮ A consequence of Kantorovich duality. Theorem: κ solves (FF) iff ψ = log( κ) maximizes Φ(ψ) :=
i
c (yi)[c(x, yi) + ψi] d µ(x) −
i ψiνi
with c(x, y) = − log(1 − x|y).
10
10
◮ ∂+Φ(λ) = {v ∈ Rd, Φ(µ) ≤ Φ(λ) + µ − λ|v ∀µ ∈ Rd}.
10
◮ In this case : ∂+Φ(λ) = {∇Φ(λ)} a.e. ◮ ∂+Φ(λ) = {v ∈ Rd, Φ(µ) ≤ Φ(λ) + µ − λ|v ∀µ ∈ Rd}. ◮ Φ concave ⇔ ∀λ ∈ Rd ∂+Φ(λ) = ∅. ◮ λ maximum of Φ ⇔ 0 ∈ ∂+Φ(λ)
10
Φ(ψ) :=
i
c (yi)[c(x, yi) + ψi] d µ(x) −
i ψiνi
10
Φ(ψ) :=
i
c (yi)[c(x, yi) + ψi] d µ(x) −
i ψiνi
=
i ψiνi
10
Φ(ψ) :=
i
c (yi)[c(x, yi) + ψi] d µ(x) −
i ψiνi
=
i ψiνi
For all ϕ ∈ Rd min1≤i≤N[c(x, yi) + ϕi] ≤ [c(x, yTψ(x)) + ϕTψ(x)]
10
Φ(ψ) :=
i
c (yi)[c(x, yi) + ψi] d µ(x) −
i ψiνi
=
i ψiνi
For all ϕ ∈ Rd min1≤i≤N[c(x, yi) + ϕi] ≤ [c(x, yTψ(x)) + ϕTψ(x)] Tψ(x) = i ⇔ x ∈ Vorψ
c (yi)
10
Φ(ψ) :=
i
c (yi)[c(x, yi) + ψi] d µ(x) −
i ψiνi
=
i ψiνi
For all ϕ ∈ Rd min1≤i≤N[c(x, yi) + ϕi] ≤ [c(x, yTψ(x)) + ϕTψ(x)] ≤ [c(x, yTψ(x)) + ψTψ(x)] + ϕTψ(x) − ψTψ(x) Tψ(x) = i ⇔ x ∈ Vorψ
c (yi)
10
Φ(ψ) :=
i
c (yi)[c(x, yi) + ψi] d µ(x) −
i ψiνi
=
i ψiνi
For all ϕ ∈ Rd min1≤i≤N[c(x, yi) + ϕi] ≤ [c(x, yTψ(x)) + ϕTψ(x)] ≤ [c(x, yTψ(x)) + ψTψ(x)] + ϕTψ(x) − ψTψ(x) Φ(ψ) +
i ψiνi
Φ(ϕ) +
i ϕiνi
Tψ(x) = i ⇔ x ∈ Vorψ
c (yi)
10
Φ(ψ) :=
i
c (yi)[c(x, yi) + ψi] d µ(x) −
i ψiνi
=
i ψiνi
Φ(ψ) − Φ(ϕ) ≤
i(ϕi − ψi)νi
Tψ(x) = i ⇔ x ∈ Vorψ
c (yi)
10
Φ(ψ) :=
i
c (yi)[c(x, yi) + ψi] d µ(x) −
i ψiνi
=
i ψiνi
Φ(ψ) − Φ(ϕ) ≤
i(ϕi − ψi)νi
≤
c (yi)
d µ(x) − νi
Tψ(x) = i ⇔ x ∈ Vorψ
c (yi)
10
Φ(ψ) :=
i
c (yi)[c(x, yi) + ψi] d µ(x) −
i ψiνi
=
i ψiνi
Φ(ψ) − Φ(ϕ) ≤
i(ϕi − ψi)νi
≤
c (yi)
d µ(x) − νi
Tψ(x) = i ⇔ x ∈ Vorψ
c (yi)
= DΦ(ψ)|ϕ − ψ with DΦ(ψ) =
c (yi)) − νi
10
Φ(ψ) :=
i
c (yi)[c(x, yi) + ψi] d µ(x) −
i ψiνi
=
i ψiνi
Φ(ψ) − Φ(ϕ) ≤
i(ϕi − ψi)νi
≤
c (yi)
d µ(x) − νi
Tψ(x) = i ⇔ x ∈ Vorψ
c (yi)
= DΦ(ψ)|ϕ − ψ with DΦ(ψ) =
c (yi)) − νi
◮ DΦ(ψ) ∈ ∂+Φ(λ) ⇒ Φ concave. ◮ ψ maximum of Φ ⇔ µ(Vorψ
c (yi)) = νi ∀i
11
12
Computation of descent direction / time step
LBFGS: low-storage version of the BFGS quasi-Newton scheme
◮ Quasi-Newton scheme:
12
c (p) d µ(x)
c (y) c(x, y) d µ(x)
Computation of descent direction / time step
LBFGS: low-storage version of the BFGS quasi-Newton scheme
Main difficulty: computation of Vorψ
c (y)
◮ Evaluation of Φ and ∇Φ: ◮ Quasi-Newton scheme:
12
c (p) d µ(x)
c (y) c(x, y) d µ(x)
Computation of descent direction / time step
LBFGS: low-storage version of the BFGS quasi-Newton scheme
Main difficulty: computation of Vorψ
c (y)
◮ Evaluation of Φ and ∇Φ: ◮ Quasi-Newton scheme:
13
Definition: Given P = {pi}1≤i≤N ⊆ Rd and (ωi)1≤i≤N ∈ RN Powω
P (pi) := {x ∈ Rd; i = arg minj x − pj2 + ωj}
13
Definition: Given P = {pi}1≤i≤N ⊆ Rd and (ωi)1≤i≤N ∈ RN Powω
P (pi) := {x ∈ Rd; i = arg minj x − pj2 + ωj}
◮ Efficient computation of (Powω
P (pi))i using CGAL (d = 2, 3)
13
Lemma: With ψ = log( κ), pi := − yj
2κj and ωi := − yj 2κj 2 − 1 κj ,
Definition: Given P = {pi}1≤i≤N ⊆ Rd and (ωi)1≤i≤N ∈ RN Powω
P (pi) := {x ∈ Rd; i = arg minj x − pj2 + ωj}
Vorψ
c (yi) = Powω P (pi) ∩ S2
◮ Efficient computation of (Powω
P (pi))i using CGAL (d = 2, 3)
13
Proof: x ∈ Vorψ
c (yi) ⊆ S2
ψ = log( κ), pi := − yj
2κj and ωi := − yj 2κj 2 − 1 κj ,
Definition: Given P = {pi}1≤i≤N ⊆ Rd and (ωi)1≤i≤N ∈ RN Powω
P (pi) := {x ∈ Rd; i = arg minj x − pj2 + ωj}
Vorψ
c (yi) = Powω P (pi) ∩ S2
⇐ ⇒ i ∈ arg minj
κj 1−x|yj
◮ Efficient computation of (Powω
P (pi))i using CGAL (d = 2, 3)
13
Proof: x ∈ Vorψ
c (yi) ⊆ S2
ψ = log( κ), pi := − yj
2κj and ωi := − yj 2κj 2 − 1 κj ,
⇐ ⇒ i ∈ arg minjx| yj
κj − 1 κj
Definition: Given P = {pi}1≤i≤N ⊆ Rd and (ωi)1≤i≤N ∈ RN Powω
P (pi) := {x ∈ Rd; i = arg minj x − pj2 + ωj}
Vorψ
c (yi) = Powω P (pi) ∩ S2
⇐ ⇒ i ∈ arg minj
κj 1−x|yj
◮ Efficient computation of (Powω
P (pi))i using CGAL (d = 2, 3)
13
Proof: x ∈ Vorψ
c (yi) ⊆ S2
ψ = log( κ), pi := − yj
2κj and ωi := − yj 2κj 2 − 1 κj ,
⇐ ⇒ i ∈ arg minjx| yj
κj − 1 κj
⇐ ⇒ i ∈ arg minj x +
yj 2κj 2 − yj 2κj 2 − 1 κj
−pj ωj
Definition: Given P = {pi}1≤i≤N ⊆ Rd and (ωi)1≤i≤N ∈ RN Powω
P (pi) := {x ∈ Rd; i = arg minj x − pj2 + ωj}
Vorψ
c (yi) = Powω P (pi) ∩ S2
⇐ ⇒ i ∈ arg minj
κj 1−x|yj
◮ Efficient computation of (Powω
P (pi))i using CGAL (d = 2, 3)
13
Proof: x ∈ Vorψ
c (yi) ⊆ S2
ψ = log( κ), pi := − yj
2κj and ωi := − yj 2κj 2 − 1 κj ,
⇐ ⇒ i ∈ arg minjx| yj
κj − 1 κj
⇐ ⇒ i ∈ arg minj x +
yj 2κj 2 − yj 2κj 2 − 1 κj
−pj ωj
⇐ ⇒ x ∈ Powω
P (pi) ∩ S2
Definition: Given P = {pi}1≤i≤N ⊆ Rd and (ωi)1≤i≤N ∈ RN Powω
P (pi) := {x ∈ Rd; i = arg minj x − pj2 + ωj}
Vorψ
c (yi) = Powω P (pi) ∩ S2
⇐ ⇒ i ∈ arg minj
κj 1−x|yj
◮ Efficient computation of (Powω
P (pi))i using CGAL (d = 2, 3)
14
◮ in general, the cells Ci := Powω
P (pi) ∩ S2 can
be disconnected, have holes, etc.
14
◮ in general, the cells Ci := Powω
P (pi) ∩ S2 can
be disconnected, have holes, etc. ◮ boundary representation: a family of oriented cycles composed of circular arcs per cell.
14
◮ in general, the cells Ci := Powω
P (pi) ∩ S2 can
be disconnected, have holes, etc. ◮ boundary representation: a family of oriented cycles composed of circular arcs per cell. ◮ lower complexity bound: Ω(N log N).
14
Algorithm: for each cell Ci = Powω
P (pi) ∩ S2
◮ in general, the cells Ci := Powω
P (pi) ∩ S2 can
be disconnected, have holes, etc. ◮ boundary representation: a family of oriented cycles composed of circular arcs per cell. ◮ lower complexity bound: Ω(N log N).
14
Algorithm: for each cell Ci = Powω
P (pi) ∩ S2
every edge of Ci and S2. Set vertices V := { }. ◮ in general, the cells Ci := Powω
P (pi) ∩ S2 can
be disconnected, have holes, etc. ◮ boundary representation: a family of oriented cycles composed of circular arcs per cell. ◮ lower complexity bound: Ω(N log N).
14
Algorithm: for each cell Ci = Powω
P (pi) ∩ S2
every edge of Ci and S2. Set vertices V := { }.
and construct oriented edges E between vertices. ◮ in general, the cells Ci := Powω
P (pi) ∩ S2 can
be disconnected, have holes, etc. ◮ boundary representation: a family of oriented cycles composed of circular arcs per cell. ◮ lower complexity bound: Ω(N log N).
14
Algorithm: for each cell Ci = Powω
P (pi) ∩ S2
every edge of Ci and S2. Set vertices V := { }.
and construct oriented edges E between vertices.
◮ in general, the cells Ci := Powω
P (pi) ∩ S2 can
be disconnected, have holes, etc. ◮ boundary representation: a family of oriented cycles composed of circular arcs per cell. ◮ lower complexity bound: Ω(N log N).
14
Algorithm: for each cell Ci = Powω
P (pi) ∩ S2
every edge of Ci and S2. Set vertices V := { }.
and construct oriented edges E between vertices.
◮ in general, the cells Ci := Powω
P (pi) ∩ S2 can
be disconnected, have holes, etc. ◮ boundary representation: a family of oriented cycles composed of circular arcs per cell. ◮ lower complexity bound: Ω(N log N).
14
Algorithm: for each cell Ci = Powω
P (pi) ∩ S2
every edge of Ci and S2. Set vertices V := { }.
and construct oriented edges E between vertices.
◮ in general, the cells Ci := Powω
P (pi) ∩ S2 can
be disconnected, have holes, etc. ◮ boundary representation: a family of oriented cycles composed of circular arcs per cell. ◮ lower complexity bound: Ω(N log N). Complexity: O(N log N + C) where C = complexity of the Power diagram.
15
ν = N
i=1 νiδxi obtained by discretizing a picture of G. Monge.
µ = uniform measure on half-sphere S2
+
N = 1000 drawing of (Vorψ
c (yi)) (on S2 +) for ψ = 0
15
ν = N
i=1 νiδxi obtained by discretizing a picture of G. Monge.
µ = uniform measure on half-sphere S2
+
N = 1000 drawing of (Vorψ
c (yi)) (on S2 +) for ψsol
15
ν = N
i=1 νiδxi obtained by discretizing a picture of G. Monge.
µ = uniform measure on half-sphere S2
+
N = 1000 rendering of the image reflected at infinity (using LuxRender)
16
ν = N
i=1 νiδxi obtained by discretizing a picture of G. Monge.
µ = uniform measure on half-sphere S2
+
N = 15000 drawing of (Vorψ
c (yi)) (on S2 +) for ψsol
16
ν = N
i=1 νiδxi obtained by discretizing a picture of G. Monge.
µ = uniform measure on half-sphere S2
+
N = 15000 solution to the far-field reflector problem: R(κsol)
16
ν = N
i=1 νiδxi obtained by discretizing a picture of G. Monge.
µ = uniform measure on half-sphere S2
+
N = 15000 rendering of the image reflected at infinity (using LuxRender)
17
18
Theorem: For N paraboloids, the complexity of the diagram (PIi( κ))1≤i≤N is O(N).
18
Theorem: For N paraboloids, the complexity of the diagram (PIi( κ))1≤i≤N is O(N). Complexity: E + F + V , where E = # edges V = # vertices F = total # of connected components
18
Theorem: For N paraboloids, the complexity of the diagram (PIi( κ))1≤i≤N is O(N). Proof:
◮ F ≤ N
18
Theorem: For N paraboloids, the complexity of the diagram (PIi( κ))1≤i≤N is O(N). Proof:
R( κ) ∩ ∂P3(κ3) PI3( κ)
P1 P3
◮ F ≤ N
Pi Pj
Lemma: The projection of ∂Pi∩ Pj onto the plane {y⊥
i } is a disc.
{y3}⊥
18
Theorem: For N paraboloids, the complexity of the diagram (PIi( κ))1≤i≤N is O(N). Proof:
R( κ) ∩ ∂P3(κ3) PI3( κ)
P1 P3
◮ F ≤ N
Pi Pj
Lemma: The projection of ∂Pi∩ Pj onto the plane {y⊥
i } is a disc.
{y3}⊥
= ⇒ the projection of R( κ) ∩ ∂Pi on {yi}⊥ is convex
18
Theorem: For N paraboloids, the complexity of the diagram (PIi( κ))1≤i≤N is O(N). Proof:
R( κ) ∩ ∂P3(κ3) PI3( κ)
P1 P3
◮ F ≤ N
Pi Pj
Lemma: The projection of ∂Pi∩ Pj onto the plane {y⊥
i } is a disc.
{y3}⊥
= ⇒ the projection of R( κ) ∩ ∂Pi on {yi}⊥ is convex = ⇒ PIi( κ) is connected.
18
Theorem: For N paraboloids, the complexity of the diagram (PIi( κ))1≤i≤N is O(N). Proof:
◮ F ≤ N
◮ Every vertex has 3 edges ⇒ 3V ≤ 2E.
18
Theorem: For N paraboloids, the complexity of the diagram (PIi( κ))1≤i≤N is O(N). Proof:
◮ F ≤ N
◮ Every vertex has 3 edges ⇒ 3V ≤ 2E. ◮ Euler’s formula V − E + F = 2 implies V ≤ 2F − 4 and E ≤ 3F − 6.
19
Proposition: Computing (PIi( κ))i requires at least Ω(N log N) operations.
19
Proposition: Computing (PIi( κ))i requires at least Ω(N log N) operations. Proof: reduction to a sorting problem
19
Proposition: Computing (PIi( κ))i requires at least Ω(N log N) operations. Proof: reduction to a sorting problem
◮ Let t1, . . . , tN ∈ R ti
19
Proposition: Computing (PIi( κ))i requires at least Ω(N log N) operations. Proof: reduction to a sorting problem
◮ Let t1, . . . , tN ∈ R ◮ yi = ϕ(ti) ∈ S2 and κi = cste. ti yi ϕ
19
Proposition: Computing (PIi( κ))i requires at least Ω(N log N) operations. Proof: reduction to a sorting problem
◮ Let t1, . . . , tN ∈ R ◮ yi = ϕ(ti) ∈ S2 and κi = cste. ◮ PIi( κ) = Powω
P (−yi) ∩ S2
with pi = −yi and ωi = cste. ti yi ϕ pi
19
Proposition: Computing (PIi( κ))i requires at least Ω(N log N) operations. Proof: reduction to a sorting problem
◮ Let t1, . . . , tN ∈ R ◮ yi = ϕ(ti) ∈ S2 and κi = cste. ◮ PIi( κ) = Powω
P (−yi) ∩ S2
with pi = −yi and ωi = cste. ti yi ϕ pi
19
Proposition: Computing (PIi( κ))i requires at least Ω(N log N) operations. Proof: reduction to a sorting problem
◮ Let t1, . . . , tN ∈ R ◮ yi = ϕ(ti) ∈ S2 and κi = cste. ◮ PIi( κ) = Powω
P (−yi) ∩ S2
with pi = −yi and ωi = cste. ti yi ϕ pi
◮ There exists a cycle of size N in the dual
◮ There exists a cycle of size N in the dual
19
Proposition: Computing (PIi( κ))i requires at least Ω(N log N) operations. Proof: reduction to a sorting problem
◮ Let t1, . . . , tN ∈ R ◮ yi = ϕ(ti) ∈ S2 and κi = cste. ◮ PIi( κ) = Powω
P (−yi) ∩ S2
with pi = −yi and ωi = cste. ti yi ϕ pi
◮ There exists a cycle of size N in the dual
◮ There exists a cycle of size N in the dual
◮ Finding the cycle ⇔ sorting t1, . . . tN
20
21
y1 y2 y3 Punctual light at origin o, µ measure on S2
i νiδyi on S2 ∞
21
y1 y2 y3 Punctual light at origin o, µ measure on S2
i νiδyi on S2 ∞
P1
P2 µ Pi(κi) = convex hull of paraboloid with focal o, direction yi and focal distance κi
21
y1 y2 y3 Punctual light at origin o, µ measure on S2
i νiδyi on S2 ∞
P1
P2 µ Pi(κi) = convex hull of paraboloid with focal o, direction yi and focal distance κi R( κ) = ∂
i=1Pi(κi)
κ) ∩ ∂P3(κ3)
21
y1 y2 y3 Punctual light at origin o, µ measure on S2
i νiδyi on S2 ∞
Pi(κi) = convex hull of paraboloid with focal o, direction yi and focal distance κi R( κ) = ∂
i=1Pi(κi)
κ) ∩ ∂P3(κ3)
21
y1 y2 y3 Punctual light at origin o, µ measure on S2
i νiδyi on S2 ∞
Pi(κi) = convex hull of paraboloid with focal o, direction yi and focal distance κi R( κ) = ∂
i=1Pi(κi)
κ) = πS2
κ) ∩ ∂Pi(κi))
PU3( κ) R( κ) ∩ ∂P3(κ3)
21
y1 y2 y3 Punctual light at origin o, µ measure on S2
i νiδyi on S2 ∞
Pi(κi) = convex hull of paraboloid with focal o, direction yi and focal distance κi R( κ) = ∂
i=1Pi(κi)
κ) = πS2
κ) ∩ ∂Pi(κi)) Far-field reflector antenna problem:
PU3( κ) R( κ) ∩ ∂P3(κ3) Problem (FF’): Find κ1, . . . , κN such that for every i, µ(PUi( κ)) = νi.
22
with pi :=
yj 2κi and ωi := − yi 2κi 2 + 1 κi .
Lemma: PUi( κ) = Powω
P (pi) ∩ S2
y1 y2 y3
PU3( κ) R( κ) ∩ ∂P3(κ3)
22
with pi :=
yj 2κi and ωi := − yi 2κi 2 + 1 κi .
Lemma: PUi( κ) = Powω
P (pi) ∩ S2
y1 y2 y3
PU3( κ) R( κ) ∩ ∂P3(κ3)
Complexity:
22
with pi :=
yj 2κi and ωi := − yi 2κi 2 + 1 κi .
Lemma: PUi( κ) = Powω
P (pi) ∩ S2
y1 y2 y3
PU3( κ) R( κ) ∩ ∂P3(κ3)
Pi Pj
◮ Projection of ∂Pi∩Pj onto y⊥
i is a disc.
Complexity:
22
with pi :=
yj 2κi and ωi := − yi 2κi 2 + 1 κi .
Lemma: PUi( κ) = Powω
P (pi) ∩ S2
y1 y2 y3
PU3( κ) R( κ) ∩ ∂P3(κ3)
Pi Pj
◮ Projection of ∂Pi∩Pj onto y⊥
i is a disc.
◮ However, the projection of R( κ) ∩ ∂Pi is not always connected. Complexity:
22
with pi :=
yj 2κi and ωi := − yi 2κi 2 + 1 κi .
Lemma: PUi( κ) = Powω
P (pi) ∩ S2
y1 y2 y3
PU3( κ) R( κ) ∩ ∂P3(κ3)
Pi Pj
◮ Projection of ∂Pi∩Pj onto y⊥
i is a disc.
◮ However, the projection of R( κ) ∩ ∂Pi is not always connected. Complexity: ⇒ unknown complexity
23
y1 y2 y3 Punctual light at origin o, µ measure on S2
i νiδyi on R3
23
y1 y2 y3 Punctual light at origin o, µ measure on S2
i νiδyi on R3
and yi, and eccentricity ei
23
y1 y2 y3 Punctual light at origin o, µ measure on S2
i νiδyi on R3
µ
and yi, and eccentricity ei
23
y1 y2 y3 Punctual light at origin o, µ measure on S2
i νiδyi on R3
µ
and yi, and eccentricity ei R( e) = ∂
i=1Ei(ei)
23
y1 y2 y3 Punctual light at origin o, µ measure on S2
i νiδyi on R3
and yi, and eccentricity ei R( e) = ∂
i=1Ei(ei)
e) = πS2
e) ∩ ∂Ei(κi)) EI1( e)
23
y1 y2 y3 Punctual light at origin o, µ measure on S2
i νiδyi on R3
and yi, and eccentricity ei R( e) = ∂
i=1Ei(ei)
e) = πS2
e) ∩ ∂Ei(κi)) Near-field reflector antenna problem: Problem (NF): Find e1, . . . , eN such that for every i, µ(EIi( e)) = νi .
amount of light reflected to the point yi.
EI1( e)
Oliker ’04
24
ρi : x ∈ S2
di 1−eix|yi/yi where di = yi(1−e2
i )
2ei
. ◮ ∂Ei(ei) is parameterized in radial coordinates by
24
ρi : x ∈ S2
di 1−eix|yi/yi where di = yi(1−e2
i )
2ei
. ◮ ∂Ei(ei) is parameterized in radial coordinates by
e) ⇐ ⇒ i ∈ arg minj
dj 1−ejx|yj/yj
24
⇐ ⇒ i ∈ arg minj log(dj) − log(1 − ejx|yj/yj) ρi : x ∈ S2
di 1−eix|yi/yi where di = yi(1−e2
i )
2ei
. ◮ ∂Ei(ei) is parameterized in radial coordinates by
e) ⇐ ⇒ i ∈ arg minj
dj 1−ejx|yj/yj
24
⇐ ⇒ i ∈ arg minj log(dj) − log(1 − ejx|yj/yj) ρi : x ∈ S2
di 1−eix|yi/yi where di = yi(1−e2
i )
2ei
. non-linearity ◮ ∂Ei(ei) is parameterized in radial coordinates by
e) ⇐ ⇒ i ∈ arg minj
dj 1−ejx|yj/yj
24
⇐ ⇒ i ∈ arg minj log(dj) − log(1 − ejx|yj/yj) ρi : x ∈ S2
di 1−eix|yi/yi where di = yi(1−e2
i )
2ei
. non-linearity ◮ ∂Ei(ei) is parameterized in radial coordinates by
◮ x ∈ EIi( e) ⇐ ⇒ i ∈ arg minj
dj 1−ejx|yj/yj
25
y1 y2
e)
25
y1 y2
e) with pi := −
eiyi 2diyi and ωi := − e2
i
2d2
i − 1
di .
Lemma: EIi( e) = Powω
P (pi) ∩ S2
25
y1 y2
e) with pi := −
eiyi 2diyi and ωi := − e2
i
2d2
i − 1
di .
Lemma: EIi( e) = Powω
P (pi) ∩ S2
⇒ yi =
−4 (ωi+pi2)2−4pi2 pi and ei = 2pi −ωi−pi2 .
◮ Inversion of pi := −
eiyi 2diyi and ωi := −
2di
2 − 1
di :
25
y1 y2
e) with pi := −
eiyi 2diyi and ωi := − e2
i
2d2
i − 1
di .
Lemma: EIi( e) = Powω
P (pi) ∩ S2
⇒ yi =
−4 (ωi+pi2)2−4pi2 pi and ei = 2pi −ωi−pi2 .
◮ Inversion of pi := −
eiyi 2diyi and ωi := −
2di
2 − 1
di :
◮ ei ∈ (0, 1) ⇐ ⇒ ωi < −pi2 and ωi < 1 − (pi + 1)2. Can always be satisfied by adding a negative constant to ωi
25
y1 y2
e) with pi := −
eiyi 2diyi and ωi := − e2
i
2d2
i − 1
di .
Lemma: EIi( e) = Powω
P (pi) ∩ S2
⇒ yi =
−4 (ωi+pi2)2−4pi2 pi and ei = 2pi −ωi−pi2 .
◮ Inversion of pi := −
eiyi 2diyi and ωi := −
2di
2 − 1
di :
◮ ei ∈ (0, 1) ⇐ ⇒ ωi < −pi2 and ωi < 1 − (pi + 1)2. Can always be satisfied by adding a negative constant to ωi ◮ We can associate a family of ellipsoid to any (pi, ωi)i
26
Proposition: The complexity of (Powω
P (pi) ∩ S2) is Ω(N 2).
26
Proposition: The complexity of (Powω
P (pi) ∩ S2) is Ω(N 2).
Proof: construct (pi)1≤i≤N s.t. the diagram Vor(pi) ∩ S2 has N 2 edges.
26
radius 2 − ε length η Proposition: The complexity of (Powω
P (pi) ∩ S2) is Ω(N 2).
Proof: construct (pi)1≤i≤N s.t. the diagram Vor(pi) ∩ S2 has N 2 edges.
26
radius 2 − ε length η Proposition: The complexity of (Powω
P (pi) ∩ S2) is Ω(N 2).
Proof: construct (pi)1≤i≤N s.t. the diagram Vor(pi) ∩ S2 has N 2 edges. ◮ N points: k = ⌊N/2⌋ p1, . . . , pk : pk+1, . . . , pN :
26
radius 2 − ε length η Proposition: The complexity of (Powω
P (pi) ∩ S2) is Ω(N 2).
Proof: construct (pi)1≤i≤N s.t. the diagram Vor(pi) ∩ S2 has N 2 edges. ◮ N points: k = ⌊N/2⌋ ◮ ∀i, j, ∃ two distinct points in Vor(pi) ∩ Vor(pj) ∩ S2 pi pj p1, . . . , pk : pk+1, . . . , pN :
26
Proposition: The complexity of (Powω
P (pi) ∩ S2) is Ω(N 2).
Proof: construct (pi)1≤i≤N s.t. the diagram Vor(pi) ∩ S2 has N 2 edges. ◮ N points: k = ⌊N/2⌋ ◮ ∀i, j, ∃ two distinct points in Vor(pi) ∩ Vor(pj) ∩ S2 pi pj p1, . . . , pk : pk+1, . . . , pN :
26
Proposition: The complexity of (Powω
P (pi) ∩ S2) is Ω(N 2).
Proof: construct (pi)1≤i≤N s.t. the diagram Vor(pi) ∩ S2 has N 2 edges. ◮ N points: k = ⌊N/2⌋ ◮ ∀i, j, ∃ two distinct points in Vor(pi) ∩ Vor(pj) ∩ S2 p1, . . . , pk : pk+1, . . . , pN :
26
Proposition: The complexity of (Powω
P (pi) ∩ S2) is Ω(N 2).
Proof: construct (pi)1≤i≤N s.t. the diagram Vor(pi) ∩ S2 has N 2 edges. ◮ N points: k = ⌊N/2⌋ ◮ ∀i, j, ∃ two distinct points in Vor(pi) ∩ Vor(pj) ∩ S2 p1, . . . , pk : pk+1, . . . , pN : ⇒ O(N 2) edges
26
Proposition: The complexity of (Powω
P (pi) ∩ S2) is Ω(N 2).
Proof: construct (pi)1≤i≤N s.t. the diagram Vor(pi) ∩ S2 has N 2 edges. ◮ N points: k = ⌊N/2⌋ ◮ ∀i, j, ∃ two distinct points in Vor(pi) ∩ Vor(pj) ∩ S2 p1, . . . , pk : pk+1, . . . , pN : ⇒ O(N 2) edges Similar results for the complexity of the ellipsoid union.
27
A simple quasi-Newton scheme can be used to solve rather large (15k points) geometric instances of optimal transport.
27
A simple quasi-Newton scheme can be used to solve rather large (15k points) geometric instances of optimal transport. Power diagrams can be used to compute efficiently the c-Voronoi cells
27
A simple quasi-Newton scheme can be used to solve rather large (15k points) geometric instances of optimal transport.
◮ quantitative stability results ? ◮ Near field reflector problem Power diagrams can be used to compute efficiently the c-Voronoi cells ◮ complexity of paraboloid union ? higher dimension ?
27
A simple quasi-Newton scheme can be used to solve rather large (15k points) geometric instances of optimal transport.
◮ quantitative stability results ? ◮ Near field reflector problem
Power diagrams can be used to compute efficiently the c-Voronoi cells ◮ complexity of paraboloid union ? higher dimension ?
28
y1 y2 y3 Punctual light at origin o, µ measure on S2
i νiδyi on R3
µ
and yi, and eccentricity ei
28
y1 y2 y3 Punctual light at origin o, µ measure on S2
i νiδyi on R3
µ
and yi, and eccentricity ei R( e) = ∂
i=1Ei(ei)
28
y1 y2 y3 Punctual light at origin o, µ measure on S2
i νiδyi on R3
and yi, and eccentricity ei R( e) = ∂
i=1Ei(ei)
e) = πS2
e) ∩ ∂Ei(κi)) EU1( e)
28
y1 y2 y3 Punctual light at origin o, µ measure on S2
i νiδyi on R3
and yi, and eccentricity ei R( e) = ∂
i=1Ei(ei)
e) = πS2
e) ∩ ∂Ei(κi)) EU1( e) with pi :=
eiyi 2diyi and ωi := − e2
i
4d2
i + 1
di .
Lemma: EUi( e) = Powω
P (pi) ∩ S2
28
y1 y2 y3 Punctual light at origin o, µ measure on S2
i νiδyi on R3
and yi, and eccentricity ei R( e) = ∂
i=1Ei(ei)
e) = πS2
e) ∩ ∂Ei(κi)) EU1( e) with pi :=
eiyi 2diyi and ωi := − e2
i
4d2
i + 1
di .
Lemma: EUi( e) = Powω
P (pi) ∩ S2
Proposition: The complexity of (Powω
P (pi) ∩ S2) is Ω(N 2).
29
Tc(µ, ν) = min{
y
T −1(y) Kantorovich problem = linear relaxation of Monge’s problem:
29
Tc(µ, ν) = min{
y
= min{
y
y πy = µ}
T −1(y) Kantorovich problem = linear relaxation of Monge’s problem: Under the Twist condition. Kp(µ, ν)
29
Kantorovich duality: Kp(µ, ν) = max{
y∈Y ψ(y)νy; φ(x) − ψ(y) ≤ c(x, y)}
= maxψ
Y ψ(y)νy
Tc(µ, ν) = min{
y
= min{
y
y πy = µ}
T −1(y) Kantorovich problem = linear relaxation of Monge’s problem: Under the Twist condition. = maxψ
c (y)[c(x, y) + ψ(y)] d µ(x) −
Y ψ(y)νy
Kp(µ, ν)
29
Kantorovich duality: Kp(µ, ν) = max{
y∈Y ψ(y)νy; φ(x) − ψ(y) ≤ c(x, y)}
= maxψ
Y ψ(y)νy
:= Φ(ψ) Tc(µ, ν) = min{
y
= min{
y
y πy = µ}
T −1(y) Kantorovich problem = linear relaxation of Monge’s problem: Under the Twist condition. Unconstrained concave maximization problem. = maxψ
c (y)[c(x, y) + ψ(y)] d µ(x) −
Y ψ(y)νy
Kp(µ, ν)
30
erigot 2010
Assuming µ has density > 0 on a bounded smooth connected domain,
30
erigot 2010
νn :=
y∈Yn νn y δy
ν :=
y∈Y νyδy
weak
Assuming µ has density > 0 on a bounded smooth connected domain,
30
erigot 2010
ψn : Yn → R solution of c-OT νn :=
y∈Yn νn y δy
ν :=
y∈Y νyδy
weak
Assuming µ has density > 0 on a bounded smooth connected domain,
30
erigot 2010
ψn : Yn → R solution of c-OT νn :=
y∈Yn νn y δy
ν :=
y∈Y νyδy
weak
Then: for (yn) ∈ (Yn) converging to y ∈ Y , limn ψn(yn) = ψ∞(y) Assuming µ has density > 0 on a bounded smooth connected domain,
30
erigot 2010
ψn : Yn → R solution of c-OT νn :=
y∈Yn νn y δy
ν :=
y∈Y νyδy
weak
Then: for (yn) ∈ (Yn) converging to y ∈ Y , limn ψn(yn) = ψ∞(y) Assuming µ has density > 0 on a bounded smooth connected domain,
30
erigot 2010
ψn : Yn → R solution of c-OT νn :=
y∈Yn νn y δy
ν :=
y∈Y νyδy
weak
Initialization: Set ν0 = ν, and construct using k-means Then: for (yn) ∈ (Yn) converging to y ∈ Y , limn ψn(yn) = ψ∞(y) Assuming µ has density > 0 on a bounded smooth connected domain,
νi
30
erigot 2010
ψn : Yn → R solution of c-OT νn :=
y∈Yn νn y δy
ν :=
y∈Y νyδy
weak
Initialization: Set ν0 = ν, and construct using k-means Then: for (yn) ∈ (Yn) converging to y ∈ Y , limn ψn(yn) = ψ∞(y) Assuming µ has density > 0 on a bounded smooth connected domain, νi+1 approximation of νi with |Pi+1| ≤ |Pi|
c
νi+1 νi
30
erigot 2010
ψn : Yn → R solution of c-OT νn :=
y∈Yn νn y δy
ν :=
y∈Y νyδy
weak
Initialization: Set ν0 = ν, and construct using k-means Then: for (yn) ∈ (Yn) converging to y ∈ Y , limn ψn(yn) = ψ∞(y) Assuming µ has density > 0 on a bounded smooth connected domain, νi+1 approximation of νi with |Pi+1| ≤ |Pi|
c
Algorithm to solve OT between µ and ν:
νi+1 νi
30
erigot 2010
ψn : Yn → R solution of c-OT νn :=
y∈Yn νn y δy
ν :=
y∈Y νyδy
weak
Initialization: Set ν0 = ν, and construct using k-means Then: for (yn) ∈ (Yn) converging to y ∈ Y , limn ψn(yn) = ψ∞(y) Assuming µ has density > 0 on a bounded smooth connected domain, νi+1 approximation of νi with |Pi+1| ≤ |Pi|
c
Algorithm to solve OT between µ and ν:
νi+1 νi
ψℓ ← 0
30
erigot 2010
ψn : Yn → R solution of c-OT νn :=
y∈Yn νn y δy
ν :=
y∈Y νyδy
weak
Initialization: Set ν0 = ν, and construct using k-means Then: for (yn) ∈ (Yn) converging to y ∈ Y , limn ψn(yn) = ψ∞(y) Assuming µ has density > 0 on a bounded smooth connected domain, νi+1 approximation of νi with |Pi+1| ≤ |Pi|
c
νi
Algorithm to solve OT between µ and ν: ψi ← minimizer of Φi using LBFGS, For i = ℓ . . . 0:
νi+1 νi
starting from ψi+1 ◦ projPi+1. ψℓ ← 0