Computational Optimal Transport 笔记——第三章（2）

常用定义

树的等价命题

3.Algorithmic Foundations

3.5 A Heuristic Description of the Network Simplex

network simplex依赖于两个简单的原则：每个可行的原始解可以对应一个互补对$(f,g)$。如果这个互补对是可行的，则我们达到最优值。如果不是，可以考虑更新$P$，使得$P$仍然是可行的，同时它的互补对$(f,g)$也被更新，它更接近可行区域。

3.5.1 Obtaining a Dual Pair Complementary to P

Simplex 首先将任意极值解$P$关联一对互补对偶变量(f, g)。通过找两个向量$f$和$g$，对于$S(P)$中的任意$(i,j^{\prime })$，$f_i+g_j=C_{i,j}$。并且保证$(f,g)$是可行的。
令$s$是$S(P)$的势。因为$P$是极值点，因此$s\le n+m-1$。由于$G(P)$没有圈，$G(P)$是一个树或者一个森林，如图3.3所示。

对于 $(f,g)$ 的目标是与$P$互补，即应该满足在$n+m$个变量上的$s$个线性等式约束

其中$S(P)$的元素被列举为$(i_1,j_1^{\prime }),...,(i_s,j_s^{\prime })$。
因为$s\le n+m-1<n+m$，线性系统（3.6）总是不确定的。即在$G(P)$中对偶变量是不确定的。
考虑在$G(P)$中的一棵树。假设这个树在源点中有$k$个节点$i_1,...,i_k$，在目标节点处有$l$个节点$j_1^{\prime },...,j_l^{\prime }$，令$r\stackrel{def.}{=}k+l$，则有$r-1$条边，对应于$f$中的$k$个变量、$g$中的$l$个变量，$r-1$个线性等式。为了消除不确定性，我们可以在树中选择任意的根节点，并且将对应的对偶变量赋值为$0$。
此时一棵树中有$r-1$个变量&$r-1$个约束，可以求解出树中所有的对偶变量。
从根节点开始我们可以使用广度优先或深度优先搜索遍历树，以获得一系列简单的变量赋值，这些赋值决定了树中所有其他双变量的值，例如图 3.4。这个过程对于图$P$中所有的树都可重复进行，以获得与$P$互补的对偶变量对$(f,g)$。

3.5.2 Network Simplex Update

如果通过3.5.1的算法得到的对偶对$(f,g)$是可行的，也就说对于所有的$i,j$，我们都有$f_i+g_j\le C_{i,j}$，根据Proposition 3.3我们得到的结果是最优的。如果存在$i,j$有$f_i+g_j>C_{i,j}$，则网络单纯形算法发挥作用。我们首先初始化图G为$G(P)$，对应于可行解$P$并且增加了不满足条件的边$(i,j^{\prime })$。会有接下来两种情况：

1. G仍然是一个森林（ 当增加的$(i,j^{\prime })$存在于两棵子树之间）。在3.5.1中提出的方法仍然可以被使用在图G上，来获得新的互补对偶向量$(f,g)$。可以注意到这个增加的边（约束）移除了在$n+m$个对偶变量之间的部分不确定性，没有对原始变量$P$进行任何改变。这个更新被称作"degenerate"（退化），因为虽然$P_{i,j}$仍然为0，但是$(i,j^{\prime })$进入了图G中。$G(P)$包含在$G$中。
2. G有一个圈。我们需要去移除$G$中的边，确保$G$ 仍然是一个森林，同时也更新$P$使得$P$是可行的，并且$G(P)$包含在$G$中。通过如下方法实现：增加$P_{i,j}$（不满足条件的边 $(i,j^{\prime })$对应的权重）的值并且改进在圈中的$P$中的其他元素，与Proposition 3.4中证明的过程类似的。具体的过程如下。我们记图$G$中的圈为$(i_1,j_1^{\prime }),(j_1^{\prime },i_2),(i_2,j_2^{\prime }),...,(i_l,j_l^{\prime }),(j_l^{\prime },i_{l+1})$，其中$i_1=i_{l+1}=i$，$j_1=j$，来强调圈是从增加的边 $(i,j)$ 处开始的.
   接下来，增加所有正向边 { i , j } \{i,j\} {i,j}，减少所有反向边 { j k ′ , i k + 1 } \{ j'_k, i_{k+1} \} {jk′​,ik+1​}，获得一个更新的结果$\tilde{P}$。下式中$\theta ={\min }_k{P}_{i_{k+1},j_k}$，令$k^{\ast }=\arg {\min }_k{P}_{i_{k+1},j_k}$。相当于在图$G$中移除边$(i_{k^{\ast }+1},j_{k^{\ast }})$，使用3.5.1的方法计算新的对偶变量。
   
   

3.5.3 Inprovement of the Primal solution

虽然这未必是我们最初的动机，但可以看出，以上的操作只会提高P的成本。如果加入的边没有形成一个圈，也就是第一种情况，原始的解并未改变，当形成了圈，也就是第二种情况，$P$被更新为$\tilde{P}$，并且如下有如下公式

我们现在使用在上节迭代中形成的对偶向量$(f,g)$。他们满足$f_{i_k}+g_{j_k}=C_{i_k,j_k}$和$f_{i_{k+1}}+g_{i_k}=C_{i_{k+1},j_k}$对于所有G中的初始边都成立，并且有

因为$C_{i,j}<f_i+g_j$，因此这一项是负的。因此如果$\theta >0$，我们有

如果$\theta =0$，图$G$被改变，但是$P$并未改变。
$networksimplexalgorithm$ 可以被总结如下：初始化算法，使用NW corner rule 得到极值解$P$，用$G(P)$初始化图$G$。计算与$P$互补的对偶变量$(f,g)$（使用3.5.1描述的树结构求解）。接下来分两种情况讨论：

1. 寻找不满足约束$C-f\oplus g\ge 0$的索引对。如果没有则$P$是最优解，停止算法。如果有，进行第二步。
2. 将满不足约束的一对边$(i,j^{\prime })$加入到$G$中。如果$G$仍然没有圈，则按照$3.5.1$更新$(f,g)$。如果有圈，在圈内操作，保证$(i,j^{\prime })$被标记为正，移除有最小流值得边，如图3.5一般更新$P$和$G$，建立相应的互补对$(f,g)$。返回第一步。
   被证明，此算法是多项式时间复杂度算法。

3.6 Dual Ascent Methods

不同于 network simplex，dual ascent methods 在每次迭代中保持对偶可行解，并且通过在$f$和$g$上加稀疏向量，使得对偶问题的目标函数值逐渐增加。
Definition 3.2 对于$S\subset [n]$，$S^{\prime }\subset [m{]}^{\prime }\stackrel{def.}{=}{1}^{\prime },...,m^{\prime }$，$1_s\in R^n$索引在S中的元素为1其余元素为0，对于$1_{S^{\prime }}\in R^m$有相同的定义。
在下文中，$(f,g)$是在$R(C)$中可行对偶对，即对于所有的$(i,j^{\prime })\in [n]\times [m{]}^{\prime },f_i+g_j\le C_{i,j}$。我们称$(i,j^{\prime })$是平衡对，如果$f_i+g_j=C_{i,j}$，其他的称为不活跃的。接下来我们说明一个可行对偶对$(f,g)$可以被以$S$和$S^{\prime }$的集合为索引的稀疏向量扰动，仍然是可行的。
Proposition 3.5 如果对于所有的$i\in S$，$(i,j^{\prime })$是平衡的（意味着$j^{\prime }\in S^{\prime }$）
,$对于足够小的$\varepsilon>0$是可行的对偶变量，
都有 $(\tilde{f},\tilde{g})\stackrel{def.}{=}(f,g)+\epsilon (1_S,-1_{S^{\prime }})$.
在3.5.1中 network simplex的迭代背后的动力被获得。为了取得这个目标，$P$逐步被更新以至于它的互补对偶对达到可行解。一个对称的方法，从可行对偶变量开始来获得可行初始问题的解$P$，推动了dual ascent methods。接下来的proposition是dual ascent methods 的主要推动，它保证了$(f,g)$的可行上升方向的存在性。
Proposition 3.6要不然$(f,g)$是问题3.4的最优解

否则存在$S\subset [n],S^{\prime }\subset [m{]}^{\prime }$使得$(\tilde{f},\tilde{g})\stackrel{def.}{=}(f,g)+\epsilon (1_S,-1_{S^{\prime }})$ 是可行的 （对于任意小的$\epsilon$）并且有更好的目标函数值。

dual ascent method 通过由集合$S,S^{\prime }$产生的任意向量更新可行解$(f,g)$。当集合$(S,S^{\prime })$由Proposition 3.6的证明过程构造。步长$\epsilon$在Proposition 3.5被定义。我们推导了一种方法称为 primal-dual method。对于匹配问题，此方法退化为 Hungarian algorithm。
两种算法的不同点：

• Simplex-type methods总是保证目前的解是可行集合的极值点，寻找满足$R(C)$的对偶解。然而dual ascent 没有做这个假设，迭代可以自由生成位于可行集内部的解。
• dual network simplex通过更新$(f,g)$来产生一个初始问题的解$P$满足边界约束但是在收敛时是非负的。dual ascent 建立$P$，它是非负的，但是不需要满足约束条件。

3.7 Auction Algorithm（拍卖算法）

Complementary slackness（互补松弛性）可以注意到在optimal assignment problem，对于 optimal transport problem提出的原始对偶条件是更加容易公式化，因为任何极值解$P$都是排列矩阵$P_{\sigma }$，对于给定的$\sigma$。给定最优的原始的$P_{{\sigma }^{\ast }}$和对偶解$f^{\ast },g^{\ast }$，我们有：

由于在3.2中阐明的 C-transforms 原则， 一个人可以选择 $f^{\ast }=g^{\bar{C}}$。因此，我们有

相反，如果存在满足满足如下条件的向量$g$和排列$\sigma$，则他们都是最优的，也就是说$\sigma$是最优分配，$g^{\bar{C}}$和$g$是最优的对偶对。

Partial assignment and $\epsilon -$ complementary slackness auction algorithm 的目标是迭代地更新三元组 $S,\xi ,g$，其中$S$是$[n]$的子集，$\xi$是一个部分分配向量，也就是从$S$到$[n]$的单射，$g$是一个对偶向量。对偶向量$g$收敛于满足近似互补松弛条件（3.8），$S$增长逐步包含$[n]$，$\xi$收敛于一个排列。在每次迭代以后算法保持一下三个性质：

算法更新 给定一个点$j$，auction algorithm 不仅使用了在一般 C-transform 中的最优值，而且使用了 second best 来定义接下来在$g$上的更新，对于索引$i\notin S$，$S$和$\xi$。

1. 更新$g$：将$g$中$j_i^1$的元素移动$\epsilon$、最小和第二小的调整的cost （${C_{i,j}-g_j}_j$）的加和
   
2. 更新$S$和$\xi$：如果存在一个索引$i^{\prime }\in S$有${\xi }_{i^{\prime }}=j_i^1$，移除它
   
   设${\xi }_i=j_i^2$并且在$S$中加入$i$。（感觉原始文献中存在问题，自己修改了一下公式）
   
   算法性质 算法从空的分配点集$S=\varphi$、空的部分分配向量$\xi$，和$g=0_n$开始，当$S=[n]$时停止，循环执行上述两个步骤，直到终止。从公式（3.9）中可以明显地看出特性（b）和（c）成立。$\epsilon -$互补松弛在第一次迭代（$S=\varphi$）是成立。通过下述的proposition可以知道特性（a）成立。

Proposition 3.7 auction algorithm在每步迭代中保存了$\epsilon -$ 互补松弛。

Proposition 3.8 auction algorithm最多执行$N=n\Vert C{\Vert }_{\infty }/\epsilon$步。

remark 3.3这个算法终止前进行了$N^3\Vert C{\Vert }_{\infty }/\epsilon$此简单的操作。当使用 $\epsilon -$scaling method时，复杂度会减少$N^3\Vert C{\Vert }_{\infty }$（在每次迭代中减少$\epsilon$的数值）。

Proposition 3.9 auction algorithm 找到了一个分配，它的费用是$n\epsilon$次优的。

综上所述 auction algorithm可以被视为一个使用 C-transform机制的方法。接下来我们将探讨另一种基于正则化的方法，即所谓的Sinkhorn算法，它与auction算法也有相似之处。最后请注意，在欧几里得空间的低维规则网格上，可以将这些经典线性求解器与多尺度策略相结合，以获得显著的加速效果