Computational Optimal Transport 笔记——第三章(1)
符号说明
I n \mathcal{I}_n In :规模为 n n n 的单位阵
常用定义
3.Algorithmic Foundations
本章介绍了对于组合优化和线性规划常见的算法工具,用于解决优化运输的离散形式。回顾第二章的内容,离散问题的原始形式(2.11):
对偶形式(2.20):
线性规划算法的目标是求解约束条件和目标函数对于优化变量是线性的最优化问题。
3.1 The Kantorovich Linear Programs
原始的优化问题为
下面我们将这个问题变换成为线性规划的标准型。
令
I
n
\mathcal{I}_n
In 代表规模为
n
n
n 的单位阵,
⊗
\otimes
⊗ 为 Kronecker’s product。
(
n
+
m
)
×
n
m
(n+m) \times nm
(n+m)×nm的矩阵
被用来编码行和与列和限制。 简化矩阵
P
∈
R
n
×
m
P \in R^{n \times m}
P∈Rn×m为一个向量
p
∈
R
n
m
p \in R^{nm}
p∈Rnm,其中
p
p
p 的第
i
+
n
(
j
−
1
)
i+n(j-1)
i+n(j−1)个元素等于
P
i
,
j
P_{i,j}
Pi,j(将
P
P
P矩阵拉长以后的列向量)。可以证得
原始问题可以写为
其中
c
c
c为
C
C
C矩阵按照列拉伸以后的
n
m
nm
nm 维向量。
Remark 3.1 上面所使用的
n
+
m
n+m
n+m 个约束是冗余的,也就是说矩阵A的向量不是线性独立的。矩阵的前
n
n
n行相加和后
m
m
m行相加得到相同的向量。
因此移除
A
A
A的一行也会得到相同的线性系统。
对于公式(3.2)的对偶问题为
注意到,这个规划与公式(2.4)是等价的。
Remark 3.2 我们提供了对偶性结果的一个简单推导,可看为 Remark 2.21 的公式化说明。强烈的对偶性,即原始(3.2)和对偶(3.3)问题的最优性确实重合的事实,需要更长的证明。为了简化说明,我们记
现在考虑一个最优运输问题的松弛原始问题,其中约束
A
p
=
q
Ap=q
Ap=q用 cost
h
T
(
A
p
−
q
)
h^{T} (Ap-q)
hT(Ap−q)代替,
h
∈
R
n
+
m
h\in R^{n+m}
h∈Rn+m为任意的参数向量。这个松弛问题可以被记为(依赖于 cost vector
h
h
h)
此问题中求解的是
H
(
h
)
H(h)
H(h) 最小而不是
z
ˉ
=
L
C
(
a
,
b
)
\bar{z}=L_{C}(a,b)
zˉ=LC(a,b)。令
p
∗
p^{*}
p∗为原始问题(3.1)的解,我们有
上述方法定义了一个问题,可以用于计算原始问题(3.1)的上界,对于任何的 cost vector h。 那个函数被称作
L
L
L 的拉格朗日对偶函数。对偶理论的目标是对于任何 cost vector
h
h
h,通过最大化
H
H
H,计算最优下界
z
‾
\underline{z}
z
第二项是关于
p
p
p的最小化,如果
c
T
−
A
T
h
c^{T}-A^{T}h
cT−ATh的某坐标是负的,第二项会趋于
−
∞
-\infty
−∞。当尝试最大化
H
(
h
)
H(h)
H(h)的下界时,限制
h
h
h满足
A
T
h
≤
c
A^{T}h \leq c
ATh≤c才有意义,在这种情况下最好的下界为
因此,我们证明了一个弱对偶结果
3.2 C-Transforms
本节,我们展示了对偶优化运输问题(3.3)的重要性质,可以被用于 semidiscrete optimal transport problem中。
本节在原始公式(2.20)上建立的,它根据行和列的加和约束分割了对偶变量。
推导
考虑对偶可行对
(
f
,
g
)
(f,g)
(f,g)。如果我们固定
f
f
f的值,我们可以发现对于
g
g
g没有比
f
f
f的 C-transform vector 更好的向量解,用
f
C
∈
R
m
f^{C} \in R^{m}
fC∈Rm来表示
可以证明
(
f
,
f
c
)
∈
R
(
C
)
(f,f^{c}) \in R(C)
(f,fc)∈R(C),并且
f
C
f^{C}
fC是满足约束条件的最大的可能的向量(满足
f
⊕
g
≤
C
f \oplus g\leq C
f⊕g≤C,最大的
g
g
g)。我们因此会有
这一结果使我们可以将对偶问题重新表述为用单变量
f
f
f表示的分段仿射凹最大化问题
将上述结果放到一边,如果我们现在固定
g
g
g的数值,考虑
g
g
g的
C
ˉ
\bar{C}
Cˉ-transform,也就是向量
g
C
ˉ
∈
R
n
g^{\bar{C}} \in R^{n}
gCˉ∈Rn
有
从给定
f
f
f开始,可以交替
C
C
C和
C
ˉ
\bar{C}
Cˉ变换多次来改进
f
f
f。我们会有如下不等式
但是在每次迭代的时候不会严格递增,因为交替
C
C
C和
C
ˉ
\bar{C}
Cˉ变换会很快达到plateau.
Proposition 3.1 下面的恒等式,其中向量之间的不等号应该在元素层面被理解:
3.3互补松弛
原始问题(3.2)和对偶问题(3.3),(2.20)可以被独立求解获得 原始的优化解 P ∗ P^{*} P∗和 对偶问题的解 ( f ∗ , g ∗ ) (f^{*},g^{*}) (f∗,g∗)。下面的 proposition 说明了他们之间的关系。
Proposition 3.2 令
P
∗
P^{*}
P∗和
(
f
∗
,
g
∗
)
(f^{*},g^{*})
(f∗,g∗)分别是原始(2.11)和对偶(2.20)问题的最优解。则对于任何一对
(
i
,
j
)
∈
[
n
]
×
[
m
]
(i,j) \in [n] \times [m]
(i,j)∈[n]×[m],满足
P
i
,
j
∗
(
C
i
,
j
−
f
i
∗
−
g
i
,
j
∗
)
=
0
P^{*}_{i,j}(C_{i,j}-f^{*}_{i}-g^{*}_{i,j})=0
Pi,j∗(Ci,j−fi∗−gi,j∗)=0
换句话说,如果
P
i
,
j
>
0
P_{i,j}>0
Pi,j>0,则
f
i
∗
+
g
j
∗
=
C
i
,
j
f^{*}_{i}+g^{*}_{j}=C_{i,j}
fi∗+gj∗=Ci,j,如果
f
i
∗
+
g
j
∗
<
C
i
,
j
f^{*}_{i}+g^{*}_{j}<C_{i,j}
fi∗+gj∗<Ci,j,则
P
i
,
j
∗
=
0
P^{*}_{i,j}=0
Pi,j∗=0。
定义 3.1(松弛) 一个矩阵 P ∈ R n × m P \in R^{n \times m} P∈Rn×m 和一对向量 ( f , g ) (f,g) (f,g)关于 C C C是松弛的也就是说对于任意的索引对 ( i , j ) (i,j) (i,j),如果 P i , j > 0 P_{i,j}>0 Pi,j>0,则有 C i , j = f i + g j C_{i,j}=f_i+g_j Ci,j=fi+gj。
Proposition 3.3 如果 P P P和 ( f , g ) (f,g) (f,g)对于原始问题和对偶问题是互补的可行解,则 P P P和 ( f , g ) (f,g) (f,g)是原始和对偶问题的最优解。(Proposition 3.2的逆命题也是成立的)
3.4 运输多面体的顶点(极点)
本节是从图的角度分析运输优化问题
回忆一个凸集的顶点或极点是点
x
x
x,满足如果集合中存在
y
y
y和
z
z
z有
x
=
(
y
+
z
)
/
2
x=(y+z)/2
x=(y+z)/2则必有
x
=
y
=
z
x=y=z
x=y=z。具有非空有界可行集的线性规划在可行集的一个顶点(或极值点)处达到最小。因为原始的优化运输问题的可行解
U
(
a
,
b
)
U(a,b)
U(a,b)是有界的,可以限制
P
P
P的搜索在多面体
U
(
a
,
b
)
U(a,b)
U(a,b)的极点。在
U
(
a
,
b
)
U(a,b)
U(a,b)的极点处的矩阵
P
P
P有一个有趣的结构,这种结构需要使用二部图的形式来描述运输问题。
3.4.1 U ( a , b ) U(a,b) U(a,b)的所有顶点的支撑集的树形结构
令
V
=
(
1
,
2
,
.
.
.
,
n
)
V=(1,2,...,n)
V=(1,2,...,n)和
V
′
=
(
1
′
,
2
′
,
.
.
.
,
m
′
)
V'=(1',2',...,m')
V′=(1′,2′,...,m′)是两个点集。考虑他们的并集
V
∩
V
′
V \cap V'
V∩V′,有
n
+
m
n+m
n+m个点,并且他们之间
n
m
nm
nm条有向边集合
ε
\varepsilon
ε
将每条边
(
i
,
j
′
)
(i,j')
(i,j′)与 cost value
C
i
j
C_{ij}
Cij相关联。在
V
V
V和
V
′
V'
V′间的完全二部图
G
\mathcal{G}
G为
一个优化计划是在图上的一个流,需要满足 source(从节点
i
i
i流出量为
a
i
a_i
ai) 和 sink(流向节点
j
′
j'
j′的量为
b
j
b_j
bj) 约束。例如Figure 3.1.
在
U
(
a
,
b
)
U(a,b)
U(a,b)上的一个极点有接下来的性质。
Proposition 3.4(Extremal solutions) 令
P
P
P是多面体
U
(
a
,
b
)
U(a,b)
U(a,b)的一个极点。令
S
(
P
)
⊂
v
a
r
e
p
s
i
l
o
n
S(P) \subset varepsilon
S(P)⊂varepsilon 是边集的子集
{
(
i
,
j
′
)
,
i
∈
[
n
]
,
j
∈
[
m
]
,
并
且
P
i
,
j
>
0
}
\{(i,j'), i \in [n], j \in [m],并且 P_{i,j}>0\}
{(i,j′),i∈[n],j∈[m],并且Pi,j>0}。则图
没有圈。特别地,
P
P
P非零元的个数不超过
n
+
m
−
1
n+m-1
n+m−1个。
3.4.2 The North-West Corner Rule
东北角规则是一个启发式规则,可以在不大于
n
+
m
n+m
n+m次操作时,产生多面体
U
(
a
,
b
)
U(a,b)
U(a,b)的一个顶点。这个规则在初始化算法求解原始优化运输问题时发挥了重大作用。
Start
- 对 P 1 , 1 P_{1,1} P1,1赋最高可能值,通过令 P 1 , 1 = m i n ( a 1 , b 1 ) P_{1,1}=min(a_1, b_1) P1,1=min(a1,b1)。
Iteration
- 选择 P i , j P_{i,j} Pi,j使得第 i i i行或者第 j j j列的约束饱和,或者两者都饱和。
- 索引根据如下规则更新:第一种情况下增加 i i i,第二种情况下增加 j j j,第三种情况下同时增加 i i i和 j j j。
Stop
- 当 P n , m P_{n,m} Pn,m被赋值。
算法具体的工作流程如下:
- i i i和 j j j被赋值为1, r ← a 1 r \gets a_1 r←a1, c ← b 1 c \gets b_1 c←b1。
- 当 i ≤ n i \leq n i≤n 和 j ≤ m j \leq m j≤m,设 t ← m i n ( r , c ) t \gets min(r,c) t←min(r,c), P i , j ← t P_{i,j} \gets t Pi,j←t, r ← r − t r \gets r-t r←r−t, c → c − t c \to c-t c→c−t;
- 如果 r = 0 r=0 r=0则增加 i i i,更新 r ← a i r \gets a_i r←ai 判断 i ≤ n i \leq n i≤n;如果 c = 0 c=0 c=0,增加 j j j,更新 c ← b j c \gets b_j c←bj 判断 j ≤ n j \leq n j≤n;
- 重复此过程
通过这个启发式算法,我们可以得到
N
W
(
a
,
b
)
NW(a,b)
NW(a,b)。
可以注意到通过以任意的顺序重排
a
,
b
a, b
a,b(
σ
,
σ
′
\sigma, \sigma'
σ,σ′),计算对应的 corner table, 然后通过转化行列顺序(
σ
−
1
,
σ
′
−
1
\sigma^{-1}, \sigma'^{-1}
σ−1,σ′−1)得到
U
(
a
,
b
)
U(a,b)
U(a,b)中的矩阵,因此通过North-West Corner Rule我们可以获得大量的解。
令
N
(
a
,
b
)
\mathcal{N}(a,b)
N(a,b)代表通过这种方法获得的所有NW corner solution
构造出的所有 NW corner solution 至多有
n
+
m
−
1
n+m-1
n+m−1个非0元素。
N
(
a
,
b
)
\mathcal{N}(a,b)
N(a,b)是
U
(
a
,
b
)
U(a,b)
U(a,b)的极点的子集。