推土机距离到Wassertein距离-程序员宅基地
推土机距离（Earth Mover’s Distance）
对于离散概率分布，推土机距离又称为 W a s s e r s t e i n \mathrm{Wasserstein} Wasserstein距离。如果将不同的概率分布看成不同的沙土堆，则推土机距离就是将一个沙土堆转换成另一个沙土堆所需的最小总工作量。假定有两个离散的概率分布 x ∼ P r x\sim P_r x∼Pr和 y ∼ P θ y\sim P_\theta y∼Pθ，其中每个概率分布都有 l l l种可能结果，如下图所示给出的两个概率分布特例
计算推土机距离是一个优化问题，从一个沙土堆到另一个沙土堆无数种方案进沙土传输迁移，所以需要找到一个最佳的传输方案 γ ( x , y ) \gamma(x,y) γ(x,y)。根据上图实例，则需要如下约束条件 { ∑ x γ ( x , y ) = P θ ( y ) ∑ y γ ( x , y ) = P r ( x ) \left\{\begin{aligned}\sum\limits_{x}\gamma(x,y)=P_\theta (y)\\\sum\limits_{y}\gamma(x,y)=P_r(x)\end{aligned}\right. ⎩⎪⎪⎪⎨⎪⎪⎪⎧x∑γ(x,y)=Pθ(y)y∑γ(x,y)=Pr(x) γ ( x , y ) ∈ Π ( P r , P θ ) \gamma(x,y)\in \Pi(P_r,P_\theta) γ(x,y)∈Π(Pr,Pθ)为联合概率分布，并且 Π ( P r , P θ ) \Pi(P_r,P_\theta) Π(Pr,Pθ)的边缘分布为 P r , P θ P_r,P_\theta Pr,Pθ。此时推土机距离为每一个 γ ( x , y ) \gamma(x,y) γ(x,y)乘以 x x x到 y y y的欧式距离之和的最小值，具体公式为 E M D ( P r , P θ ) = inf γ ∈ Π ∑ x , y ∥ x − y ∥ γ ( x , y ) = inf γ ∈ Π E ( x , y ) ∼ γ ∥ x − y ∥ \mathrm{EMD}(P_r,P_\theta)=\inf\limits_{\gamma\in \Pi}\sum\limits_{x,y}\|x-y\|\gamma(x,y)=\inf\limits_{\gamma\in \Pi}\mathbb{E}_{(x,y)\sim\gamma}\|x-y\| EMD(Pr,Pθ)=γ∈Πinfx,y∑∥x−y∥γ(x,y)=γ∈ΠinfE(x,y)∼γ∥x−y∥进一步令 Γ = γ ( x , y ) \Gamma=\gamma(x,y) Γ=γ(x,y)， D = ∥ x , y ∥ D=\|x,y\| D=∥x,y∥，其中 Γ , D ∈ R l × l \Gamma,D\in \mathbb{R}^{l\times l} Γ,D∈Rl×l，则上式可以化简为 E M D ( P r , P θ ) = inf γ ∈ Π ⟨ D , γ ⟩ F \mathrm{EMD}(P_r,P_\theta)=\inf\limits_{\gamma\in \Pi}\langle D,\gamma \rangle_{F} EMD(Pr,Pθ)=γ∈Πinf⟨D,γ⟩F其中 ⟨ , ⟩ \langle , \rangle ⟨,⟩表示斐波那契内积，即对应元素相乘再相加，矩阵 Γ \Gamma Γ和 D D D的热力示意图如下所示
线性规划求解
求解推土机距离中的最优传输方案可以利用线性规划标准型来求解。给定一个向量 x ∈ R n x\in \mathbb{R}^n x∈Rn，线性目标标准型的优化形式如下所示 min x z = c ⊤ x s . t . { A x = b x ≥ 0 \begin{array}{rl}\min\limits_{x}&z=c^{\top}x\\\mathrm{s.t.}&\left\{\begin{aligned}Ax&=b\\x&\ge 0\end{aligned}\right.\end{array} xmins.t.z=c⊤x{
Axx=b≥0其中 c ∈ R n c\in\mathbb{R}^n c∈Rn， A = ∈ R m × n A=\in \mathbb{R}^{m\times n} A=∈Rm×n， b ∈ R m b\in \mathbb{R}^m b∈Rm。根据以上实例将推土机距离转化为线性规划标准型，首先将矩阵 Γ \Gamma Γ和 D D D进行向量化，则有 x = v e c ( Γ ) ∈ R n c = v e c ( D ) ∈ R n \begin{aligned}x&=\mathrm{vec}(\Gamma)\in\mathbb{R}^{n}\\c&=\mathrm{vec}(D)\in\mathbb{R}^{n}\end{aligned} xc=vec(Γ)∈Rn=vec(D)∈Rn并且有 n = l 2 n=l^2 n=l2，将目标分布进行拼接则有 b = [ P r P θ ] b=\left[\begin{array}{c}P_r\\P_\theta\end{array}\right] b=[PrPθ]其中 m = 2 l m=2l m=2l，方程组 A x = b Ax=b Ax=b的具体形式如下所示 { γ ( x 1 , y 1 ) + γ ( x 1 , y 2 ) + γ ( x 1 , y 3 ) + γ ( x 1 , y 4 ) + γ ( x 1 , y 5 ) + γ ( x 1 , y 6 ) + γ ( x 1 , y 7 ) + γ ( x 1 , y 8 ) + γ ( x 1 , y 9 ) + γ ( x 1 , y 10 ) = P r ( x 1 ) γ ( x 2 , y 1 ) + γ ( x 2 , y 2 ) + γ ( x 2 , y 3 ) + γ ( x 2 , y 4 ) + γ ( x 2 , y 5 ) + γ ( x 2 , y 6 ) + γ ( x 2 , y 7 ) + γ ( x 2 , y 8 ) + γ ( x 2 , y 9 ) + γ ( x 2 , y 10 ) = P r ( x 2 ) γ ( x 3 , y 1 ) + γ ( x 3 , y 2 ) + γ ( x 3 , y 3 ) + γ ( x 3 , y 4 ) + γ ( x 3 , y 5 ) + γ ( x 3 , y 6 ) + γ ( x 3 , y 7 ) + γ ( x 3 , y 8 ) + γ ( x 3 , y 9 ) + γ ( x 3 , y 10 ) = P r ( x 3 ) γ ( x 4 , y 1 ) + γ ( x 4 , y 2 ) + γ ( x 4 , y 3 ) + γ ( x 4 , y 4 ) + γ ( x 4 , y 5 ) + γ ( x 4 , y 6 ) + γ ( x 4 , y 7 ) + γ ( x 4 , y 8 ) + γ ( x 4 , y 9 ) + γ ( x 4 , y 10 ) = P r ( x 4 ) γ ( x 5 , y 1 ) + γ ( x 5 , y 2 ) + γ ( x 5 , y 3 ) + γ ( x 5 , y 4 ) + γ ( x 5 , y 5 ) + γ ( x 5 , y 6 ) + γ ( x 5 , y 7 ) + γ ( x 5 , y 8 ) + γ ( x 5 , y 9 ) + γ ( x 5 , y 10 ) = P r ( x 5 ) γ ( x 6 , y 1 ) + γ ( x 6 , y 2 ) + γ ( x 6 , y 3 ) + γ ( x 6 , y 4 ) + γ ( x 6 , y 5 ) + γ ( x 6 , y 6 ) + γ ( x 6 , y 7 ) + γ ( x 6 , y 8 ) + γ ( x 6 , y 9 ) + γ ( x 6 , y 10 ) = P r ( x 6 ) γ ( x 7 , y 1 ) + γ ( x 7 , y 2 ) + γ ( x 7 , y 3 ) + γ ( x 7 , y 4 ) + γ ( x 7 , y 5 ) + γ ( x 7 , y 6 ) + γ ( x 7 , y 7 ) + γ ( x 7 , y 8 ) + γ ( x 7 , y 9 ) + γ ( x 7 , y 10 ) = P r ( x 7 ) γ ( x 8 , y 1 ) + γ ( x 8 , y 2 ) + γ ( x 8 , y 3 ) + γ ( x 8 , y 4 ) + γ ( x 8 , y 5 ) + γ ( x 8 , y 6 ) + γ ( x 8 , y 7 ) + γ ( x 8 , y 8 ) + γ ( x 8 , y 9 ) + γ ( x 8 , y 10 ) = P r ( x 8 ) γ ( x 9 , y 1 ) + γ ( x 9 , y 2 ) + γ ( x 9 , y 3 ) + γ ( x 9 , y 4 ) + γ ( x 9 , y 5 ) + γ ( x 9 , y 6 ) + γ ( x 9 , y 7 ) + γ ( x 9 , y 8 ) + γ ( x 9 , y 9 ) + γ ( x 9 , y 10 ) = P r ( x 9 ) γ ( x 10 , y 1 ) + γ ( x 10 , y 2 ) + γ ( x 10 , y 3 ) + γ ( x 10 , y 4 ) + γ ( x 10 , y 5 ) + γ ( x 10 , y 6 ) + γ ( x 10 , y 7 ) + γ ( x 10 , y 8 ) + ⋯ + γ ( x 10 , y 10 ) = P r ( x 10 ) γ ( x 1 , y 1 ) + γ ( x 2 , y 1 ) + γ ( x 3 , y 1 ) + γ ( x 4 , y 1 ) + γ ( x 5 , y 1 ) + γ ( x 6 , y 1 ) + γ ( x 7 , y 1 ) + γ ( x 8 , y 1 ) + γ ( x 9 , y 1 ) + γ ( x 10 , y 1 ) = P θ ( y 1 ) γ ( x 1 , y 2 ) + γ ( x 2 , y 2 ) + γ ( x 3 , y 2 ) + γ ( x 4 , y 2 ) + γ ( x 5 , y 2 ) + γ ( x 6 , y 2 ) + γ ( x 7 , y 2 ) + γ ( x 8 , y 2 ) + γ ( x 9 , y 2 ) + γ ( x 10 , y 2 ) = P θ ( y 2 ) γ ( x 1 , y 3 ) + γ ( x 2 , y 3 ) + γ ( x 3 , y 3 ) + γ ( x 4 , y 3 ) + γ ( x 5 , y 3 ) + γ ( x 6 , y 3 ) + γ ( x 7 , y 3 ) + γ ( x 8 , y 3 ) + γ ( x 9 , y 3 ) + γ ( x 10 , y 3 ) = P θ ( y 3 ) γ ( x 1 , y 4 ) + γ ( x 2 , y 4 ) + γ ( x 3 , y 4 ) + γ ( x 4 , y 4 ) + γ ( x 5 , y 4 ) + γ ( x 6 , y 4 ) + γ ( x 7 , y 4 ) + γ ( x 8 , y 4 ) + γ ( x 9 , y 4 ) + γ ( x 10 , y 4 ) = P θ ( y 4 ) γ ( x 1 , y 5 ) + γ ( x 2 , y 5 ) + γ ( x 3 , y 5 ) + γ ( x 4 , y 5 ) + γ ( x 5 , y 5 ) + γ ( x 6 , y 5 ) + γ ( x 7 , y 5 ) + γ ( x 8 , y 5 ) + γ ( x 9 , y 5 ) + γ ( x 10 , y 5 ) = P θ ( y 5 ) γ ( x 1 , y 6 ) + γ ( x 2 , y 6 ) + γ ( x 3 , y 6 ) + γ ( x 4 , y 6 ) + γ ( x 5 , y 6 ) + γ ( x 6 , y 6 ) + γ ( x 7 , y 6 ) + γ ( x 8 , y 6 ) + γ ( x 9 , y 6 ) + γ ( x 10 , y 6 ) = P θ ( y 6 ) γ ( x 1 , y 7 ) + γ ( x 2 , y 7 ) + γ ( x 3 , y 7 ) + γ ( x 4 , y 7 ) + γ ( x 5 , y 7 ) + γ ( x 6 , y 7 ) + γ ( x 7 , y 7 ) + γ ( x 8 , y 7 ) + γ ( x 9 , y 7 ) + γ ( x 10 , y 7 ) = P θ ( y 7 ) γ ( x 1 , y 8 ) + γ ( x 2 , y 8 ) + γ ( x 3 , y 8 ) + γ ( x 4 , y 8 ) + γ ( x 5 , y 8 ) + γ ( x 6 , y 8 ) + γ ( x 7 , y 8 ) + γ ( x 8 , y 8 ) + γ ( x 9 , y 8 ) + γ ( x 10 , y 8 ) = P θ ( y 8 ) γ ( x 1 , y 9 ) + γ ( x 2 , y 9 ) + γ ( x 3 , y 9 ) + γ ( x 4 , y 9 ) + γ ( x 5 , y 9 ) + γ ( x 6 , y 9 ) + γ ( x 7 , y 9 ) + γ ( x 8 , y 9 ) + γ ( x 9 , y 9 ) + γ ( x 10 , y 9 ) = P θ ( y 9 ) γ ( x 1 , y 10 ) + γ ( x 2 , y 10 ) + γ ( x 3 , y 10 ) + γ ( x 4 , y 10 ) + γ ( x 5 , y 10 ) + γ ( x 6 , y 10 ) + γ ( x 7 , y 10 ) + γ ( x 8 , y 10 ) + ⋯ + γ ( x 10 , y 10 ) = P θ ( y 10 ) \left\{\begin{aligned}\gamma(x_1,y_1)+\gamma(x_1,y_2)+\gamma(x_1,y_3)+\gamma(x_1,y_4)+\gamma(x_1,y_5)+\gamma(x_1,y_6)+\gamma(x_1,y_7)+\gamma(x_1,y_8)+\gamma(x_1,y_9)+\gamma(x_1,y_{10})&=P_r(x_1)\\\gamma(x_2,y_1)+\gamma(x_2,y_2)+\gamma(x_2,y_3)+\gamma(x_2,y_4)+\gamma(x_2,y_5)+\gamma(x_2,y_6)+\gamma(x_2,y_7)+\gamma(x_2,y_8)+\gamma(x_2,y_9)+\gamma(x_2,y_{10})&=P_r(x_2)\\ \gamma(x_3,y_1)+\gamma(x_3,y_2)+\gamma(x_3,y_3)+\gamma(x_3,y_4)+\gamma(x_3,y_5)+\gamma(x_3,y_6)+\gamma(x_3,y_7)+\gamma(x_3,y_8)+\gamma(x_3,y_9)+\gamma(x_3,y_{10})&=P_r(x_3)\\\gamma(x_4,y_1)+\gamma(x_4,y_2)+\gamma(x_4,y_3)+\gamma(x_4,y_4)+\gamma(x_4,y_5)+\gamma(x_4,y_6)+\gamma(x_4,y_7)+\gamma(x_4,y_8)+\gamma(x_4,y_9)+\gamma(x_4,y_{10})&=P_r(x_4)\\\gamma(x_5,y_1)+\gamma(x_5,y_2)+\gamma(x_5,y_3)+\gamma(x_5,y_4)+\gamma(x_5,y_5)+\gamma(x_5,y_6)+\gamma(x_5,y_7)+\gamma(x_5,y_8)+\gamma(x_5,y_9)+\gamma(x_5,y_{10})&=P_r(x_5)\\\gamma(x_6,y_1)+\gamma(x_6,y_2)+\gamma(x_6,y_3)+\gamma(x_6,y_4)+\gamma(x_6,y_5)+\gamma(x_6,y_6)+\gamma(x_6,y_7)+\gamma(x_6,y_8)+\gamma(x_6,y_9)+\gamma(x_6,y_{10})&=P_r(x_6)\\ \gamma(x_7,y_1)+\gamma(x_7,y_2)+\gamma(x_7,y_3)+\gamma(x_7,y_4)+\gamma(x_7,y_5)+\gamma(x_7,y_6)+\gamma(x_7,y_7)+\gamma(x_7,y_8)+\gamma(x_7,y_9)+\gamma(x_7,y_{10})&=P_r(x_7)\\\gamma(x_8,y_1)+\gamma(x_8,y_2)+\gamma(x_8,y_3)+\gamma(x_8,y_4)+\gamma(x_8,y_5)+\gamma(x_8,y_6)+\gamma(x_8,y_7)+\gamma(x_8,y_8)+\gamma(x_8,y_9)+\gamma(x_8,y_{10})&=P_r(x_8)\\\gamma(x_{9},y_1)+\gamma(x_{9},y_2)+\gamma(x_{9},y_3)+\gamma(x_{9},y_4)+\gamma(x_{9},y_5)+\gamma(x_{9},y_6)+\gamma(x_{9},y_7)+\gamma(x_{9},y_8)+\gamma(x_{9},y_9)+\gamma(x_{9},y_{10})&=P_r(x_{9})\\\gamma(x_{10},y_1)+\gamma(x_{10},y_2)+\gamma(x_{10},y_3)+\gamma(x_{10},y_4)+\gamma(x_{10},y_5)+\gamma(x_{10},y_6)+\gamma(x_{10},y_7)+\gamma(x_{10},y_8)+\cdots+\gamma(x_{10},y_{10})&=P_r(x_{10})\\\gamma(x_1,y_1)+\gamma(x_2,y_1)+\gamma(x_3,y_1)+\gamma(x_4,y_1)+\gamma(x_5,y_1)+\gamma(x_6,y_1)+\gamma(x_7,y_1)+\gamma(x_8,y_1)+\gamma(x_9,y_1)+\gamma(x_{10},y_{1})&=P_\theta(y_1)\\\gamma(x_1,y_2)+\gamma(x_2,y_2)+\gamma(x_3,y_2)+\gamma(x_4,y_2)+\gamma(x_5,y_2)+\gamma(x_6,y_2)+\gamma(x_7,y_2)+\gamma(x_8,y_2)+\gamma(x_9,y_2)+\gamma(x_{10},y_{2})&=P_\theta(y_2)\\\gamma(x_1,y_3)+\gamma(x_2,y_3)+\gamma(x_3,y_3)+\gamma(x_4,y_3)+\gamma(x_5,y_3)+\gamma(x_6,y_3)+\gamma(x_7,y_3)+\gamma(x_8,y_3)+\gamma(x_9,y_3)+\gamma(x_{10},y_{3})&=P_\theta(y_3)\\\gamma(x_1,y_4)+\gamma(x_2,y_4)+\gamma(x_3,y_4)+\gamma(x_4,y_4)+\gamma(x_5,y_4)+\gamma(x_6,y_4)+\gamma(x_7,y_4)+\gamma(x_8,y_4)+\gamma(x_9,y_4)+\gamma(x_{10},y_{4})&=P_\theta(y_4)\\\gamma(x_1,y_5)+\gamma(x_2,y_5)+\gamma(x_3,y_5)+\gamma(x_4,y_5)+\gamma(x_5,y_5)+\gamma(x_6,y_5)+\gamma(x_7,y_5)+\gamma(x_8,y_5)+\gamma(x_9,y_5)+\gamma(x_{10},y_{5})&=P_\theta(y_5)\\\gamma(x_1,y_6)+\gamma(x_2,y_6)+\gamma(x_3,y_6)+\gamma(x_4,y_6)+\gamma(x_5,y_6)+\gamma(x_6,y_6)+\gamma(x_7,y_6)+\gamma(x_8,y_6)+\gamma(x_9,y_6)+\gamma(x_{10},y_{6})&=P_\theta(y_6)\\\gamma(x_1,y_7)+\gamma(x_2,y_7)+\gamma(x_3,y_7)+\gamma(x_4,y_7)+\gamma(x_5,y_7)+\gamma(x_6,y_7)+\gamma(x_7,y_7)+\gamma(x_8,y_7)+\gamma(x_9,y_7)+\gamma(x_{10},y_{7})&=P_\theta(y_7)\\\gamma(x_1,y_8)+\gamma(x_2,y_8)+\gamma(x_3,y_8)+\gamma(x_4,y_8)+\gamma(x_5,y_8)+\gamma(x_6,y_8)+\gamma(x_7,y_8)+\gamma(x_8,y_8)+\gamma(x_9,y_8)+\gamma(x_{10},y_{8})&=P_\theta(y_8)\\\gamma(x_1,y_9)+\gamma(x_2,y_9)+\gamma(x_3,y_9)+\gamma(x_4,y_9)+\gamma(x_5,y_9)+\gamma(x_6,y_9)+\gamma(x_7,y_9)+\gamma(x_8,y_9)+\gamma(x_9,y_9)+\gamma(x_{10},y_{9})&=P_\theta(y_9)\\\gamma(x_1,y_{10})+\gamma(x_2,y_{10})+\gamma(x_3,y_{10})+\gamma(x_4,y_{10})+\gamma(x_5,y_{10})+\gamma(x_6,y_{10})+\gamma(x_7,y_{10})+\gamma(x_8,y_{10})+\cdots+\gamma(x_{10},y_{10})&=P_\theta(y_{10})\end{aligned}\right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧γ(x1,y1)+γ(x1,y2)+γ(x1,y3)+γ(x1,y4)+γ(x1,y5)+γ(x1,y6)+γ(x1,y7)+γ(x1,y8)+γ(x1,y9)+γ(x1,y10)γ(x2,y1)+γ(x2,y2)+γ(x2,y3)+γ(x2,y4)+γ(x2,y5)+γ(x2,y6)+γ(x2,y7)+γ(x2,y8)+γ(x2,y9)+γ(x2,y10)γ(x3,y1)+γ(x3,y2)+γ(x3,y3)+γ(x3,y4)+γ(x3,y5)+γ(x3,y6)+γ(x3,y7)+γ(x3,y8)+γ(x3,y9)+γ(x3,y10)γ(x4,y1)+γ(x4,y2)+γ(x4,y3)+γ(x4,y4)+γ(x4,y5)+γ(x4,y6)+γ(x4,y7)+γ(x4,y8)+γ(x4,y9)+γ(x4,y10)γ(x5,y1)+γ(x5,y2)+γ(x5,y3)+γ(x5,y4)+γ(x5,y5)+γ(x5,y6)+γ(x5,y7)+γ(x5,y8)+γ(x5,y9)+γ(x5,y10)γ(x6,y1)+γ(x6,y2)+γ(x6,y3)+γ(x6,y4)+γ(x6,y5)+γ(x6,y6)+γ(x6,y7)+γ(x6,y8)+γ(x6,y9)+γ(x6,y10)γ(x7,y1)+γ(x7,y2)+γ(x7,y3)+γ(x7,y4)+γ(x7,y5)+γ(x7,y6)+γ(x7,y7)+γ(x7,y8)+γ(x7,y9)+γ(x7,y10)γ(x8,y1)+γ(x8,y2)+γ(x8,y3)+γ(x8,y4)+γ(x8,y5)+γ(x8,y6)+γ(x8,y7)+γ(x8,y8)+γ(x8,y9)+γ(x8,y10)γ(x9,y1)+γ(x9,y2)+γ(x9,y3)+γ(x9,y4)+γ(x9,y5)+γ(x9,y6)+γ(x9,y7)+γ(x9,y8)+γ(x9,y9)+γ(x9,y10)γ(x10,y1)+γ(x10,y2)+γ(x10,y3)+γ(x10,y4)+γ(x10,y5)+γ(x10,y6)+γ(x10,y7)+γ(x10,y8)+⋯+γ(x10,y10)γ(x1,y1)+γ(x2,y1)+γ(x3,y1)+γ(x4,y1)+γ(x5,y1)+γ(x6,y1)+γ(x7,y1)+γ(x8,y1)+γ(x9,y1)+γ(x10,y1)γ(x1,y2)+γ(x2,y2)+γ(x3,y2)+γ(x4,y2)+γ(x5,y2)+γ(x6,y2)+γ(x7,y2)+γ(x8,y2)+γ(x9,y2)+γ(x10,y2)γ(x1,y3)+γ(x2,y3)+γ(x3,y3)+γ(x4,y3)+γ(x5,y3)+γ(x6,y3)+γ(x7,y3)+γ(x8,y3)+γ(x9,y3)+γ(x10,y3)γ(x1,y4)+γ(x2,y4)+γ(x3,y4)+γ(x4,y4)+γ(x5,y4)+γ(x6,y4)+γ(x7,y4)+γ(x8,y4)+γ(x9,y4)+γ(x10,y4)γ(x1,y5)+γ(x2,y5)+γ(x3,y5)+γ(x4,y5)+γ(x5,y5)+γ(x6,y5)+γ(x7,y5)+γ(x8,y5)+γ(x9,y5)+γ(x10,y5)γ(x1,y6)+γ(x2,y6)+γ(x3,y6)+γ(x4,y6)+γ(x5,y6)+γ(x6,y6)+γ(x7,y6)+γ(x8,y6)+γ(x9,y6)+γ(x10,y6)γ(x1,y7)+γ(x2,y7)+γ(x3,y7)+γ(x4,y7)+γ(x5,y7)+γ(x6,y7)+γ(x7,y7)+γ(x8,y7)+γ(x9,y7)+γ(x10,y7)γ(x1,y8)+γ(x2,y8)+γ(x3,y8)+γ(x4,y8)+γ(x5,y8)+γ(x6,y8)+γ(x7,y8)+γ(x8,y8)+γ(x9,y8)+γ(x10,y8)γ(x1,y9)+γ(x2,y9)+γ(x3,y9)+γ(x4,y9)+γ(x5,y9)+γ(x6,y9)+γ(x7,y9)+γ(x8,y9)+γ(x9,y9)+γ(x10,y9)γ(x1,y10)+γ(x2,y10)+γ(x3,y10)+γ(x4,y10)+γ(x5,y10)+γ(x6,y10)+γ(x7,y10)+γ(x8,y10)+⋯+γ(x10,y10)=Pr(x1)=Pr(x2)=Pr(x3)=Pr(x4)=Pr(x5)=Pr(x6)=Pr(x7)=Pr(x8)=Pr(x9)=Pr(x10)=Pθ(y1)=Pθ(y2)=Pθ(y3)=Pθ(y4)=Pθ(y5)=Pθ(y6)=Pθ(y7)=Pθ(y8)=Pθ(y9)=Pθ(y10)其中矩阵 A A A是一个大的 0 0 0和 1 1 1的二值稀疏矩阵为 A = [ 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 e 1 e 1 e 1 e 1 e 1 e 1 e 1 e 1 e 1 e 1 e 2 e 2 e 2 e 2 e 2 e 2 e 2 e 2 e 2 e 2 e 3 e 3 e 3 e 3 e 3 e 3 e 3 e 3 e 3 e 3 e 4 e 4 e 4 e 4 e 4 e 4 e 4 e 4 e 4 e 4 e 5 e 5 e 5 e 5 e 5 e 5 e 5 e 5 e 5 e 5 e 6 e 6 e 6 e 6 e 6 e 6 e 6 e 6 e 6 e 6 e 7 e 7 e 7 e 7 e 7 e 7 e 7 e 7 e 7 e 7 e 8 e 8 e 8 e 8 e 8 e 8 e 8 e 8 e 8 e 8 e 9 e 9 e 9 e 9 e 9 e 9 e 9 e 9 e 9 e 9 e 10 e 10 e 10 e 10 e 10 e 10 e 10 e 10 e 10 e 10 ] A=\left[\begin{array}{cccccccccc}\bf{1}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}\\\bf{0}&\bf{1}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}\\\bf{0}&\bf{0}&\bf{1}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}\\\bf{0}&\bf{0}&\bf{0}&\bf{1}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}\\\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{1}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}\\\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{1}&\bf{0}&\bf{0}&\bf{0}&\bf{0}\\\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{1}&\bf{0}&\bf{0}&\bf{0}\\\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{1}&\bf{0}&\bf{0}\\\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{1}&\bf{0}\\\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{0}&\bf{1}\\{\bf{e}}_1&{\bf{e}}_1&{\bf{e}}_1&{\bf{e}}_1&{\bf{e}}_1&{\bf{e}}_1&{\bf{e}}_1&{\bf{e}}_1&{\bf{e}}_1&{\bf{e}}_1\\{\bf{e}}_2&{\bf{e}}_2&{\bf{e}}_2&{\bf{e}}_2&{\bf{e}}_2&{\bf{e}}_2&{\bf{e}}_2&{\bf{e}}_2&{\bf{e}}_2&{\bf{e}}_2\\{\bf{e}}_3&{\bf{e}}_3&{\bf{e}}_3&{\bf{e}}_3&{\bf{e}}_3&{\bf{e}}_3&{\bf{e}}_3&{\bf{e}}_3&{\bf{e}}_3&{\bf{e}}_3\\{\bf{e}}_4&{\bf{e}}_4&{\bf{e}}_4&{\bf{e}}_4&{\bf{e}}_4&{\bf{e}}_4&{\bf{e}}_4&{\bf{e}}_4&{\bf{e}}_4&{\bf{e}}_4\\{\bf{e}}_5&{\bf{e}}_5&{\bf{e}}_5&{\bf{e}}_5&{\bf{e}}_5&{\bf{e}}_5&{\bf{e}}_5&{\bf{e}}_5&{\bf{e}}_5&{\bf{e}}_5\\{\bf{e}}_6&{\bf{e}}_6&{\bf{e}}_6&{\bf{e}}_6&{\bf{e}}_6&{\bf{e}}_6&{\bf{e}}_6&{\bf{e}}_6&{\bf{e}}_6&{\bf{e}}_6\\{\bf{e}}_7&{\bf{e}}_7&{\bf{e}}_7&{\bf{e}}_7&{\bf{e}}_7&{\bf{e}}_7&{\bf{e}}_7&{\bf{e}}_7&{\bf{e}}_7&{\bf{e}}_7\\{\bf{e}}_8&{\bf{e}}_8&{\bf{e}}_8&{\bf{e}}_8&{\bf{e}}_8&{\bf{e}}_8&{\bf{e}}_8&{\bf{e}}_8&{\bf{e}}_8&{\bf{e}}_8\\{\bf{e}}_9&{\bf{e}}_9&{\bf{e}}_9&{\bf{e}}_9&{\bf{e}}_9&{\bf{e}}_9&{\bf{e}}_9&{\bf{e}}_9&{\bf{e}}_9&{\bf{e}}_9\\{\bf{e}}_{10}&{\bf{e}}_{10}&{\bf{e}}_{10}&{\bf{e}}_{10}&{\bf{e}}_{10}&{\bf{e}}_{10}&{\bf{e}}_{10}&{\bf{e}}_{10}&{\bf{e}}_{10}&{\bf{e}}_{10}\end{array}\right] A=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡1000000000e1e2e3e4e5e6e7e8e9e100100000000e1e2e3e4e5e6e7e8e9e100010000000e1e2e3e4e5e6e7e8e9e100001000000e1e2e3e4e5e6e7e8e9e100000100000e1e2e3e4e5e6e7e8e9e100000010000e1e2e3e4e5e6e7e8e9e100000001000e1e2e3e4e5e6e7e8e9e100000000100e1e2e3e4e5e6e7e8e9e100000000010e1e2e3e4e5e6e7e8e9e100000000001e1e2e3e4e5e6e7e8e9e10⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤其中向量 1 \bf{1} 1， 0 \bf{0} 0和 e i {\bf{e}}_i ei具体取值如下所示 1 = [ 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 ] ∈ R 1 × 10 {\bf{1}}=[1,1,1,1,1,1,1,1,1,1]\in \mathbb{R}^{1\times 10} 1=[1,1,1,1,1,1,1,1,1,1]∈R1×10 0 = [ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ] ∈ R 1 × 10 {\bf{0}}=[0,0,0,0,0,0,0,0,0,0]\in \mathbb{R}^{1\times 10} 0=[0,0,0,0,0,0,0,0,0,0]∈R1×10 e i [ l ] = { 1 , l = i 0 , l ≠ i , e i ∈ R 1 × 10 , i ∈ { 1 , ⋯ , 10 } {\bf{e}}_i[l]=\left\{\begin{array}{ll}1,&l=i\\0,&l\ne i\end{array},\quad {\bf{e}}_i\in \mathbb{R}^{1\times 10}, \quad i\in\{1,\cdots,10\}\right. ei[l]={
1,0,l=il=i,ei∈R1×10,i∈{
1,⋯,10}通过求解线性规划问题得到的最优传输方案的示意图如下所示
对偶线性规划求解
当随机变量的离散状态数量与输入变量的维数呈指数关系时，线性规划标准型求解在这些情况下就很不实用，例如在深度学习图像领域中，用 G A N \mathrm{GAN} GAN生成图片其输入的维度可以很容易达到数千维，此时就不会用线性规划标准型进行求解矩阵 Γ \Gamma Γ。 G A N \mathrm{GAN} GAN的核心目的利用分布 P r P_r Pr生成一个与真实分布尽可能相同的分布 P θ P_\theta Pθ，所以并不需要关注联合概率分布 γ \gamma γ，而是关注生成分布 P r P_r Pr和真实分布 P θ P_\theta Pθ的 E M D \mathrm{EMD} EMD数值（推土机距离），然后将该距离看做成损失函数进而求解梯度 ∇ P θ E M D \nabla_{P_\theta}\mathrm{EMD} ∇PθEMD并对网络进行训练。任何线性规划的初始形式都有其对偶形式如下所示 p r i m a l f o r m : m i n i m i z e z = c ⊤ x s . t . { A x = b x ≥ 0 ∣ d u a l f o r m : m a x i m i z e z ^ = b ⊤ y s . t . A ⊤ . y ≤ c \left. {\bf{primal\text{ } form:\quad}}\begin{array}{rc}\mathrm{minimize}&z=c^{\top}x\\\mathrm{s.t.}&\left\{\begin{aligned}Ax&=b\\x &\ge 0\end{aligned}\right.\end{array}\right|{\bf{dual \text{ }form:}}\quad\begin{array}{rc}\mathrm{maximize}&\hat{z}=b^{\top}y\\\mathrm{s.t.}&A^{\top}.y \le c\\&\end{array} primal form:minimizes.t.z=c⊤x{ Axx=b≥0∣∣∣∣∣∣∣dual form:maximizes.t.z^=b⊤yA⊤.y≤c将求解初始问题的最小值转化为求解其对偶问题的最大值，目标 z ^ \hat{z} z^直接依赖于 b b b，其中 b b b包含概率分布 P r P_r Pr和 P θ P_\theta Pθ。初始线性规划 z z z的最小值是对偶线性规划 z ^ \hat{z} z^最大值的上界，具体则有 z = c ⊤ x ≥ y ⊤ A x = y ⊤ b = z ^ z=c^{\top}x\ge y^{\top}Ax=y^{\top}b=\hat{z} z=c⊤x≥y⊤Ax=y⊤b=z^以上不等式称为弱对偶定理。强对偶定理则是初始线性规划 z z z的最小值等于对偶线性规划 z ^ \hat{z} z^最大值即 z = z ^ z=\hat{z} z=z^。利用对偶形式计算推土机距离 E M D \mathrm{EMD} EMD，假定 y ∗ = [ f g ] y^{*}=\left[\begin{array}{c}{\bf{f}}\\{\bf{g}}\end{array}\right] y∗=[fg]其中 f , g ∈ R l × 1 {\bf{f,g}}\in\mathbb{R}^{l\times 1} f,g∈Rl×1，进一步根据上公式则有 E M D ( P r , P θ ) = f ⊤ P r + g ⊤ P θ \mathrm{EMD}(P_r,P_\theta)={\bf{f}}^{\top}P_r+{\bf{g}}^{\top}P_\theta EMD(Pr,Pθ)=f⊤Pr+g⊤Pθ其中向量 f {\bf{f}} f和 g {\bf{g}} g的值分别由函数 f ( ⋅ ) f(\cdot) f(⋅)和 g ( ⋅ ) g(\cdot) g(⋅)可得 { f i = f ( x i ) g i = g ( x i ) , i ∈ { 1 , ⋯ , n } \left\{\begin{aligned}{\bf{f}}_i&=f(x_i)\\{\bf{g}}_i&=g(x_i)\end{aligned}\right.,\quad i\in\{1,\cdots,n\} { figi=f(xi)=g(xi),i∈{ 1,⋯,n}对偶形式的约束为 A ⊤ y ≤ c A^{\top}y\le c A⊤y≤c，以上约束条件展开则有 f ( x i ) + g ( x j ) ≤ D i j = ∥ x i − x j ∥ , i , j ∈ { 1 , ⋯ , n } f(x_i)+g(x_j)\le D_{ij}=\|x_i-x_j\|,\quad i,j\in\{1,\cdots,n\} f(xi)+g(xj)≤Dij=∥xi−xj∥,i,j∈{ 1,⋯,n}此时对偶形式可以整理为如下形式 E M D 1 ( P r , P θ ) = max f , g { ∑ i = 1 n [ f ( x i ) P r ( x i ) + g ( x i ) P θ ( x i ) ] ∣ ∀ i , j , f ( x i ) + g ( x j ) ≤ D i j } \mathrm{EMD}^{1}(P_r,P_\theta)=\max\limits_{f,g}\left\{\left.\sum\limits_{i=1}^n[f(x_i)P_r(x_i)+g(x_i)P_{\theta}(x_i)]\right|\forall i,j, \quad f(x_i)+g(x_j)\le D_{ij} \right\} EMD1(Pr,Pθ)=f,gmax{ i=1∑n[f(xi)Pr(xi)+g(xi)Pθ(xi)]∣∣∣∣∣∀i,j,f(xi)+g(xj)≤Dij}又因为当 i = j i=j i=j时，则有 f ( x i ) + g ( x i ) ≤ D i i = 0 f(x_i)+g(x_i)\le D_{ii}=0 f(xi)+g(xi)≤Dii=0进而可知 ∑ i = 1 n [ f ( x i ) P r ( x i ) + g ( x i ) P θ ( x i ) ] ≤ ∑ i = 1 n [ f ( x i ) P r ( x i ) − f ( x i ) P θ ( x i ) ] \begin{aligned}\sum\limits_{i=1}^n[f(x_i)P_r(x_i)+g(x_i)P_{\theta}(x_i)]& \le \sum\limits_{i=1}^n[f(x_i)P_r(x_i)-f(x_i)P_{\theta}(x_i)]\end{aligned} i=1∑n[f(xi)Pr(xi)+g(xi)Pθ(xi)]≤i=1∑n[f(xi)Pr(xi)−f(xi)Pθ(xi)]所以当函数 g = − f g=-f g=−f时，则有 E M D 1 ( P r , P θ ) ≤ E M D 2 ( P r , P θ ) = max f , − f { ∑ i = 1 n [ f ( x i ) P r ( x i ) − f ( x i ) P θ ( x i ) ] ∣ ∀ i , j , f ( x i ) − f ( x j ) ≤ D i j } \mathrm{EMD}^1(P_r,P_{\theta})\le\mathrm{EMD}^{2}(P_r,P_\theta)=\max\limits_{f,-f}\left.\left\{\sum\limits_{i=1}^n[f(x_i)P_r(x_i)-f(x_i)P_{\theta}(x_i)]\right|\forall i,j,\quad f(x_i)-f(x_j)\le D_{ij}\right\} EMD1(Pr,Pθ)≤EMD2(Pr,Pθ)=f,−fmax{ i=1∑n[f(xi)Pr(xi)−f(xi)Pθ(xi)]∣∣∣∣∣∀i,j,f(xi)−f(xj)≤Dij}又因为 E M D 2 ( P r , P θ ) \mathrm{EMD}^2(P_r,P_\theta) EMD2(Pr,Pθ)中的函数取值范围 f , − f f,-f f,−f是 E M D 1 ( P r , P θ ) \mathrm{EMD}^1(P_r,P_\theta) EMD1(Pr,Pθ)中函数取值范围 f , g f,g f,g的一个特例，则有 E M D 1 ( P r , P θ ) ≥ E M D 2 ( P r , P θ ) \mathrm{EMD}^1(P_r,P_\theta)\ge \mathrm{EMD}^2(P_r,P_\theta) EMD1(Pr,Pθ)≥EMD2(Pr,Pθ)，综上所述则有 E M D ( P r , P θ ) = E M D 1 ( P r , P θ ) = E M D 2 ( P r , P θ ) \mathrm{EMD}(P_r,P_\theta)=\mathrm{EMD}^1(P_r,P_\theta)= \mathrm{EMD}^2(P_r,P_\theta) EMD(Pr,Pθ)=EMD1(Pr,Pθ)=EMD2(Pr,Pθ)又因为 ∣ f ( x i ) − f ( x j ) ≤ ∣ ∣ x i − x j ∣ ∣ |f(x_i)-f(x_j)\le ||x_i-x_j|| ∣f(xi)−f(xj)≤∣∣xi−xj∣∣所以 f f f是 1 1 1- L i p s c h i z \mathrm{Lipschiz} Lipschiz连续，所以可以将推土机距离 E M D ( P r , P θ ) \mathrm{EMD}(P_r,P_\theta) EMD(Pr,Pθ)对偶形式表示为 E M D ( P r , P θ ) = max ∥ f ∥ ≤ 1 ∑ i = 1 n [ f ( x i ) P r ( x i ) − f ( x i ) P θ ( x i ) ] = max ∥ f ∥ ≤ 1 E x ∼ P r ( x ) [ f ( x ) ] − E x ∼ P θ ( x ) [ f ( x ) ] \begin{aligned}\mathrm{EMD}(P_r,P_\theta)&=\max\limits_{\|f\|\le 1}\sum\limits_{i=1}^n[f(x_i)P_r(x_i)-f(x_i)P_{\theta}(x_i)]\\&=\max\limits_{\|f\|\le 1}\mathbb{E}_{x\sim P_{r}(x)}[f(x)]-\mathbb{E}_{x\sim P_{\theta}(x)}[f(x)]\end{aligned} EMD(Pr,Pθ)=∥f∥≤1maxi=1∑n[f(xi)Pr(xi)−f(xi)Pθ(xi)]=∥f∥≤1maxEx∼Pr(x)[f(x)]−Ex∼Pθ(x)[f(x)]
W a s s e r s t e i n \mathrm{Wasserstein} Wasserstein距离
当 W a s s e r s t e i n \mathrm{Wasserstein} Wasserstein距离取 1 1 1范数的时候则为推土机距离，以上考虑的概率分布是离散的情况。考虑随机变量是连续的情况，给定连续随机变量的分布 p r p_r pr和 p θ p_\theta pθ，并且它们是联合分布 π ( p r , p θ ) \pi(p_r,p_\theta) π(pr,pθ)的边缘分布，则 W a s s e r s t e i n \mathrm{Wasserstein} Wasserstein距离表示为 W ( p r , p θ ) = inf γ ∈ π ∫ x ∫ y ∥ x − y ∥ γ ( x , y ) d x d y = inf γ ∈ π E x , y ∼ γ [ ∥ x − y ∥ ] W(p_r,p_\theta)=\inf_{\gamma \in \pi}\int\limits_x\int\limits_{y}\|x-y\|\gamma(x,y)dxdy=\inf\limits_{\gamma \in \pi}\mathbb{E}_{x,y\sim\gamma}[\|x-y\|] W(pr,pθ)=γ∈πinfx∫y∫∥x−y∥γ(x,y)dxdy=γ∈πinfEx,y∼γ[∥x−y∥]通过引入一个额外的函数 f : x ⟼ k ∈ R f:x\longmapsto k\in \mathbb{R} f:x⟼k∈R，可以消除掉关于联合分布 γ \gamma γ所有的约束条件 W ( p r , p θ ) = inf γ ∈ π E x , y ∼ γ [ ∥ x − y ∥ ] = inf γ ∈ π E x , y ∼ γ [ ∥ x − y ∥ + sup f E s ∼ p r [ f ( s ) ] − E t ∼ p θ [ f ( t ) ] − ( f ( x ) − f ( y ) ) ] = inf γ sup f E x , y ∼ γ [ ∥ x − y ∥ + E s ∼ p r [ f ( s ) ] − E t ∼ p θ [ f ( t ) ] − ( f ( x ) − f ( y ) ) ] \begin{aligned}W(p_r,p_\theta)&=\inf\limits_{\gamma \in \pi}\mathbb{E}_{x,y\sim \gamma}[\|x-y\|]\\&=\inf\limits_{\gamma \in \pi}\mathbb{E}_{x,y\sim \gamma}[\|x-y\|+\sup\limits_{f} \mathbb{E}_{s \sim p_r}[f(s)]-\mathbb{E}_{t\sim p_\theta}[f(t)]-(f(x)-f(y))]\\&=\inf\limits_{\gamma}\sup\limits_{f}\mathbb{E}_{x,y\sim \gamma}[\|x-y\|+\mathbb{E}_{s\sim p_r}[f(s)]-\mathbb{E}_{t\sim p_{\theta}}[f(t)]-(f(x)-f(y))]\end{aligned} W(pr,pθ)=γ∈πinfEx,y∼γ[∥x−y∥]=γ∈πinfEx,y∼γ[∥x−y∥+fsupEs∼pr[f(s)]−Et∼pθ[f(t)]−(f(x)−f(y))]=γinffsupEx,y∼γ[∥x−y∥+Es∼pr[f(s)]−Et∼pθ[f(t)]−(f(x)−f(y))]其中有 sup f E s ∼ p r [ f ( s ) ] − E t ∼ p θ [ f ( t ) ] − ( f ( x ) − f ( y ) ) ] = { 0 , γ ∈ π + ∞ , o t h e r w i s e \sup\limits_{f} \mathbb{E}_{s \sim p_r}[f(s)]-\mathbb{E}_{t\sim p_\theta}[f(t)]-(f(x)-f(y))]=\left\{\begin{array}{rl}0,& \gamma \in \pi \\ +\infty,& \mathrm{otherwise}\end{array}\right. fsupEs∼pr[f(s)]−Et∼pθ[f(t)]−(f(x)−f(y))]={ 0,+∞,γ∈πotherwise以上问题是一个极大极小值双层优化问题，求解以上问题需要用到极小极大原理，在不改变解的前提下颠倒求解顺序，则有对偶形式如下所示 W ( p r , p θ ) = sup f inf γ E x , y ∼ γ [ ∥ x − y ∥ + E s ∼ p r [ f ( s ) ] − E t ∼ p θ [ f ( t ) ] − ( f ( x ) − f ( y ) ) ] = sup f E s ∼ p r [ f ( s ) ] − E t ∼ p θ [ f ( t ) ] + inf γ E x , y ∼ γ [ ∥ x − y ∥ − ( f ( x ) − f ( y ) ) ] \begin{aligned}W(p_r,p_\theta)=&\sup\limits_{f}\inf\limits_{\gamma} \mathbb{E}_{x,y\sim \gamma}[\|x-y\|+\mathbb{E}_{s\sim p_r}[f(s)]-\mathbb{E}_{t \sim p_{\theta}}[f(t)]-(f(x)-f(y))]\\=&\sup\limits_{f}\mathbb{E}_{s\sim p_r}[f(s)]-\mathbb{E}_{t\sim p_{\theta}}[f(t)]+\inf\limits_{\gamma}\mathbb{E}_{x,y \sim \gamma}[\|x-y\|-(f(x)-f(y))]\end{aligned} W(pr,pθ)==fsupγinfEx,y∼γ[∥x−y∥+Es∼pr[f(s)]−Et∼pθ[f(t)]−(f(x)−f(y))]fsupEs∼pr[f(s)]−Et∼pθ[f(t)]+γinfEx,y∼γ[∥x−y∥−(f(x)−f(y))]又因为 inf γ E x , y ∼ γ [ ∥ x − y ∥ − ( f ( x ) − f ( y ) ) ] = { 0 , ∥ f ∥ L ≤ 1 − ∞ , o t h e r w i s e \inf\limits_{\gamma}\mathbb{E}_{x,y \sim \gamma}[\|x-y\|-(f(x)-f(y))]=\left\{\begin{array}{rl}0,& \|f\|_L \le 1\\-\infty,& \mathrm{otherwise}\end{array}\right. γinfEx,y∼γ[∥x−y∥−(f(x)−f(y))]={ 0,−∞,∥f∥L≤1otherwise则可得最后对偶形式如下所示 W ( p r , p θ ) = sup ∥ f ∥ L ≤ 1 E s ∼ p r [ f ( s ) ] − E t ∼ p θ [ f ( t ) ] W(p_r,p_\theta)=\sup\limits_{\|f\|_L \le 1}\mathbb{E}_{s\sim p_r}[f(s)]-\mathbb{E}_{t\sim p_\theta}[f(t)] W(pr,pθ)=∥f∥L≤1supEs∼pr[f(s)]−Et∼pθ[f(t)]该函数 f f f适合用神经网络来逼近，而且这种方法的优点是，只需夹紧权值即可实现 L i p s c h i t z \mathrm{Lipschitz} Lipschitz连续性。
相关代码
本文中所涉及线性规划求解，以及相关的示意图程如如下所示
import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import linprog
from matplotlib import cm
import matplotlib.colors as colors
l = 10
P_r = np.array([12,7,4,1,19,14,9,6,3,2])
P_t = np.array([1,5,11,17,13,9,6,4,3,2])
P_r = P_r / np.sum(P_r)
P_t = P_t / np.sum(P_t)
plt.subplot(1,2,1)
plt.bar(range(l), P_r, 1, color='limegreen', edgecolor='black',linewidth=5, alpha=1)
plt.axis('off')
plt.ylim(0, 0.5)
plt.subplot(1,2,2)
plt.bar(range(l), P_t, 1, color='fuchsia', edgecolor='black',linewidth=5, alpha=1)
plt.axis('off')
plt.ylim(0, 0.5)
plt.tight_layout()
plt.show()
D = np.ndarray(shape=(l, l))
for i in range(l):
for j in range(l):
D[i,j] = abs(range(l)[i] - range(l)[j])
A_r = np.zeros((l, l, l))
A_t = np.zeros((l, l, l))
for i in range(l):
for j in range(l):
A_r[i, i, j] = 1
A_t[i, j, i] = 1
# print(A_r)
A = np.concatenate((A_r.reshape((l, l**2)), A_t.reshape((l, l**2))), axis=0)
print("A: \n", A.shape, "\n")
b = np.concatenate((P_r, P_t), axis=0)
c = D.reshape((l**2))
opt_res = linprog(c, A_eq = A, b_eq = b, bounds=[0, None])
emd = opt_res.fun
gamma = opt_res.x.reshape((l,l))
print('EMD:', emd, "\n")
plt.imshow(gamma, cmap=cm.rainbow, interpolation='nearest')
plt.axis('off')
plt.show()
plt.imshow(D, cmap=cm.rainbow, interpolation='nearest')
plt.axis('off')
plt.show()
opt_res = linprog(-b, A.T, c, bounds=(None,None))
emd = -opt_res.fun
f = opt_res.x[0:l]
g = opt_res.x[l:]
print('dual EMD:', emd)
print('f: \n', f)
plt.plot(range(l), f)
plt.show()
print('g: \n', g)
plt.plot(range(l), g)
plt.show()
plt.bar(range(l), np.multiply(P_r, f), 1, color='blue', alpha=1)
plt.ylim(-0.5, 0.5)
plt.axis('off')
plt.show()
plt.bar(range(l), np.multiply(P_t, g), 1, color='green', alpha=1)
plt.axis('off')
plt.ylim(-0.5, 0.5)
plt.show()
#check sum
emd = np.sum(np.multiply(P_r, f)) + np.sum(np.multiply(P_t, g))
print("emd: ", emd)
cNorm = colors.Normalize(vmin=0, vmax=l)
colorMap = cm.ScalarMappable(norm=cNorm, cmap=cm.terrain)
current_bottom = np.zeros(l)
r = range(l)
for i in r.__reversed__():
plt.bar(r, gamma[r, i], 1, color=colorMap.to_rgba(r), bottom=current_bottom)
current_bottom = current_bottom + gamma[r, i]
plt.axis('off')
plt.ylim(0, 0.5)
plt.show()
current_bottom = np.zeros(l)
for i in r:
plt.bar(r, gamma[i, r], 1, color=colorMap.to_rgba(i), bottom=current_bottom)
current_bottom = current_bottom + gamma[i, r]
plt.axis('off')
plt.ylim(0, 0.5)
plt.show()
智能推荐
freenas php5.56,FreeNAS中文网-程序员宅基地
文章浏览阅读296次。我们很高兴地宣布FreeNAS 11.3-BETA1全面上市。11.3系列代表了长达一年的开发和测试工作。与先前的测试版发布周期不同，此测试版已由iXsystems QA团队内部全面审核，用户应注意其稳定性和可用性方面的显着改进。在进行更新之前，请仔细阅读这些发行说明，以熟悉此更新所带来的许多新功能的潜在影响。请向https://jira.ixsystems.com/projects/NAS报告任..._freenas中文官网
Java设计模式（十）装饰器模式_装饰'模式java代码-程序员宅基地
文章浏览阅读868次。装饰器模式_装饰'模式java代码
Java学习手册：（数据结构与算法-字符串）如何输出字符串的所有组合？_输出a bbb ccccc的java语句源代码-程序员宅基地
文章浏览阅读533次。问题描述：假设字符串中所有的字符串都不重复，如何输出字符串的所有组合？例如，输入字符串“abc”，则输出a、b、c、ab、ac、bc、abc，共7种组合。方法一：递归法方法二：构造一个长度为n的01字符串(或二进制数)表示输出结果中是否包含某个字符。原题要求则变为输出“000”到“111”这2^n-1个组合对应的字符串。方法一代码如下：package com.haobi;/*..._输出a bbb ccccc的java语句源代码
oracle的in集合,oracle中in与not in集合中有空值问题-程序员宅基地
文章浏览阅读638次。IN和NOT IN本质上都是OR运算，因而计算逻辑OR时处理NULL的方式不同，产生的结果也不同。下面我们分析一下前面的三条语句SQL> select deptno2 from dept3 where deptno in (10,50,null);这里可以等价于where deptno=10 or deptno=50 or deptno=null，由于是or相连接，那么只要有一个条件为T..._oracle in 有空值
cdma特有效应_CDMA九问九答-程序员宅基地
文章浏览阅读198次。CDMA技术九问九答一、什么是CDMA技术CDMA直译为码分多址，是在数字通信技术的分支扩频通信的基础上发展起来的一种技术。所谓扩频，简单地说就是把频谱扩展。CDMA技术采用的是直接序列扩频方式，就是用具有噪声特性的载波以及比简单点到几点通信所需带宽宽得多的频带去传输相同的数据。同调频、调幅技术一样，直接序列扩频是一种调制技术，它采用一个码序列(高速)去调制原始数据信息(低速)，这样调制后的信息就..._cdma有阴影效应吗
场景文本检测与识别：最新进展及未来趋势_icdar 2003介绍-程序员宅基地
文章浏览阅读5k次，点赞10次，收藏42次。Scene text detection and recognition: recent advances and future trends文章目录摘要1、介绍2.场景文本检测与识别的最新进展2.1 场景文本检测相关工作2.2 场景文本识别相关工作2.3 端到端文本识别相关工作2.4 场景文本检测与识别的相关应用及产品3 基准数据集和评估协议3.1基准数据集3.2评估协议3.2.1文本检测算法的评估协议3.2.2 文本识别算法的评估协议4讨论5 结论摘要文本作为人类最具影响力的发明之一，在人类生_icdar 2003介绍
随便推点
网编SEO进阶:JavaScript对SEO的影响及解决之道-程序员宅基地
文章浏览阅读73次。JavaScript在SEO中是一个很头疼的问题,一方面我们在网页制作中需要使用JavaScript来实现绚丽的特效,而一方面JavaScript又会对搜索引擎的抓取分析造成不好的影响。...yahoo在官方文档中也有类似的强调:尽量使用搜索引擎能够识别的文本信息,避免过多的JavaScript、Cookie、框架、DHTML或Flash等复杂技术。...更多内容，请访问SEO助手博客转载于:h..._javasrciptseo代码
Python爬虫：Scrapy中间件Middleware和Pipeline_scrapy middleware pipeline-程序员宅基地
文章浏览阅读2.9k次，点赞2次，收藏3次。Scrapy提供了可自定义2种中间件，1个数据处理器名称作用用户设置数据收集器（Item-Pipeline)处理item覆盖下载中间件（Downloader-Middleware）处理request/response合并爬虫中间件（Spider-Middleware）处理item/response/request合并解释：用户设置：是指cu..._scrapy middleware pipeline
菜鸟小万的每天一道ccf题Day22（202104-01灰度直方图）_问题描述 一幅长宽分别为 个像素和 个像素的灰度图像可以表示为一个 大小的矩阵-程序员宅基地
文章浏览阅读347次。欢迎大家留言一起探讨更好的解法，有什么错误也欢迎指正呀。呜呜呜呜，dbq，连上21天课我要废了，时间只能让我完成第一题了。忙完这周必然开始第二题（flag）。标题已经加了题号了，方便查找。202104-01灰度直方图问题描述一幅长宽分别为 n 个像素和 m 个像素的灰度图像可以表示为一个 n×m 大小的矩阵 A。其中每个元素 Aij（0≤i<n、0≤j<m）是一个 [0,L) 范围内的整数，表示对应位置像素的灰度值。具体来说，一个 8 比特的灰度图像中每个像素的灰度范围是 [0,12_问题描述 一幅长宽分别为 个像素和 个像素的灰度图像可以表示为一个 大小的矩阵
校园网络的组建与设计——计算机网络大作业【计算机高级应用】-程序员宅基地
文章浏览阅读3.6w次，点赞102次，收藏700次。校园网络的组建与设计——计算机网络大作业【计算机高级应用】校园网络设计.pkt文件已在百度网盘中共享_计算机网络大作业
python读取数据校验数据_Python验证数据的抽样分布类型-程序员宅基地
文章浏览阅读155次。假如要对一份统计数据进行分析，一般其来源来自于社会调研/普查，所以数据不是总体而是一定程度的抽样。对于抽样数据的分析，就可以结合上篇统计量及其抽样分布的内容，判断数据符合哪种分布。使用已知分布特性，可以完成对总体的统计分析。本文使用python函数判断数据集是否符合特定抽样分布。数据来源本次试验使用kagglehttps://www.kaggle.com/datasets上的公开数据集，可以通过搜..._python读取数据校验数据
单词九连猜python编程_python实现猜单词小游戏-程序员宅基地
文章浏览阅读85次。Python初学者小游戏：猜单词，供大家参考，具体内容如下游戏逻辑：就像我们曾经英语学习机上的小游戏一样，电脑会从事先预置的词库中抽取单词，然后给出单词的字母数量，给定猜解次数，然后让玩家进行猜测，并给出每次猜测的正确字母与错误字母。涉及知识点：random.randint()，print()，input()（raw_input()）参考实现代码：#!/usr/bin/python# -*- co..._python 单词九连猜