【最优传输二十九】Wasserstein Barycenterand Its Application to Texture Mixing

motivation

本文提出了离散概率分布的平均作为Monge-Kantorovich最优传输空间重心的新定义。为了克服数值求解这类问题所涉及的时间复杂性，原始的Wasserstein度量被一维分布上的切片近似所取代。这使我们能够引入一种新的快速梯度下降算法来计算点云的Wasserstein质心。

这个概率重心的新概念可能会在计算机视觉中找到应用，在那里人们想要平均定义为分布的特征。我们展示了纹理合成和混合的应用，其中纹理的特征是对多尺度定向滤波器组的响应分布。这就提供了一种简单的方法来导航颜色纹理的凸域。

1.introduce

本文考虑了Monge-Kantorovich最优输运理论[1]在图像合成中的应用。自Rubner等人[2]的开创性工作以来，最优运输成本已被广泛用作比较直方图特征的度量，通常被称为“地球移动者的距离”。这里研究的另一个有趣的方面是运输映射本身。事实上，它允许各种图像修改，例如颜色转移[3]，纹理映射[4]，或视频的对比度均衡[5](如其他应用[6])。

这项工作的主要理论贡献是引入了Monge-Kantorovich空间中统计分布重心的新定义。我们还提出了一个近似的定义，更适合于依赖于一维投影的数值计算。然后引入牛顿梯度下降算法来计算相应的Wasserstein质心。本文的最后一个贡献是颜色纹理统计合成的一般框架，它包含了几个现有的纹理模型作为特殊情况。这个通用框架，连同Wasserstein质心计算，允许执行颜色纹理混合。

2. Wasserstein距离及其近似

本文考虑在 $\mathbb{R}^{d}$ 中离散密度分布，用点云 $X=\left \{ X_{i} \right \}_{i\in I}$ , $I=\left \{ 1,...,N \right \}$ 是点索引的列表。由于X的任何排列都对应于相同的分布，因此需要考虑度量:

其中 $\Sigma _{N}$ 为N个元素的所有排列的集合。本文提出的方法可以推广到加权点云和密度在 $\mathbb{R}^{d}$ 上连续定义。

2.1 Wasserstein Distance

定义两个相同大小 $\left | I \right |=N$ 的点云之间的二次Wasserstein距离 $W_{2}(X,Y)$ ，为:

本方法可以推广到任意严格凸距离，如对于p > 1, $\left \| X_{i}-Y_{\sigma (i)} \right \|^{p}$ 。可以证明 $W_{p}$ 在离散分布集合上定义了一个度量[X]。

线性规划公式。计算这个距离需要计算最小化 $W_{\sigma }(X,Y)$ in(2)的最优分配 $i \mapsto \sigma ^{\star }(i)$ 。可以将这个问题重新定义为一个线性规划问题:

其中 $P_{N}$ 为双随机矩阵的集合，即行和列之和为1的非负矩阵。问题(3)可以用标准线性规划算法和更专用的方法在 $O(N^{2.5}\log (N))$ 次操作中解决(参见[28])。

一维的情况。d = 1有一些特殊的结构，允许更快的解决方案。事实上，如果用 $\sigma _{X}$ 和 $\sigma _{Y}$ 表示点的排列顺序

最优排列 $\sigma ^{\star }$ 使公式(2)最小化

使得点 $X_{\sigma _{X}(i)}$ 分配给点 $X_{\sigma _{Y}(i)}$ 。因此，使用快速排序算法，可以在 $O(N\log (N))$ 次操作中计算出Wasserstein距离和最优分配。

2.2 Sliced Wasserstein Distance

然而，Wasserstein距离W2的计算对于我们所想到的图像处理应用程序来说计算要求太高，其中N可能相当大。此外，在涉及Wasserstein距离函数的点云优化问题中，W2太难处理。

由于这些原因，我们现在考虑一种分布之间的替代度量，它基于一维投影之间的运输成本。我们用 $SW_{2}$ 表示切片Wasserstein距离，它定义为投影点云之间的一维二次瓦瑟斯坦距离之和

其中 $\Omega =\left \{ \theta \in \mathbb{R}^{d}\setminus \left \| \theta \right \| =1\right \}$ 是单位球。切片Wasserstein距离允许我们使用一维分配的特殊情况，这种情况可以使用(5)以封闭形式轻松解决。

3 Barycenter in Wasserstein Space

3.1 Wasserstein Barycenter

给定一个点云族 $Y=\left \{ Y^{j} \right \}_{j\in J}$ ，我们感兴趣的是计算一个加权平均点云 $X^{\star }$ ，即通过类比欧几里得设置来定义为问题(B)的最小值。

$\left \{ \rho _{j}\geqslant 0 \right \}_{j\in J}$ 是一个权值的集合，它被约束满足 $\sum _{j}\rho _{j}=1$ 。

除了正态分布的特殊情况[30]外，这个问题没有已知的封闭形式解(7)。与我们的工作不同，Agueh和Carlier对这个问题进行了数学分析[31]。它们证明了连续分布的解和对偶公式的存在性。

一维的情况。在一维情况下，Wasserstein重心可以用 $O(N\log (N))$ 次运算来计算，使用排列 $\sigma _{Y^{j}}$ ，将值 $Y^{j}\subset \mathbb{R}$ 的集合排序，如(4)所示。然后重心读取

切片Wasserstein重心。如果考虑任意密度(不一定是固定数目N个点的云，甚至不一定是离散的云)上的问题(B)， Agueh和Carlier在[31]中证明，可以通过求解线性问题来找到质心。可以证明，如果N个离散点支持Y的密度，则质心最多支持 $N^{\left | J \right |}$ 个点。这种方法的困难在于计算时间是 $N^{\left | J \right |}$ 的多项式，这对于成像应用来说是禁止的，并且存在重心解的基数的组合爆炸。

此外，问题(B)对点云(式(7))的限制可以归结为一个多任务问题，不幸的是，这个问题是NP-hard。

为了同时解决这两个问题，我们定义一个近似于原质心的“切片Wasserstein质心”

注意到，SEY的最小化很容易，并且与度量集 $Y=\left \{ Y^{j} \right \}_{j\in J}$ 具有相同维度的点云的期望属性。在实践中，如第3.2节所述，通过执行该能量的梯度下降，可以快速计算出近似质心，从而得到局部最小值。

3.2 Gradient Descent Algorithm

通过最小化(8)找到重心对应于非凸泛函的最小化。该能量用有限的ω方向集进行离散，其中|ω| > d。在数值实验中，我们使用从 $\mathbb{R}^{d}$ 的单位球上的均匀分布中得到的随机方向。

用梯度下降算法可以找到该泛函的一个平稳点。

牛顿梯度下降法。该算法从某个初始点云 $X^{\left ( 0 \right )}\subset \mathbb{R}^{d}$ 开始，该点云可以选择为任意云 $Y^{j}$ , k = 0。在每次迭代k时，点云 $X^{(k)}$ 以以下牛顿下降步长更新:

其中 $\triangledown SE_{Y_{\theta }}(X^{(k)})$ 为 $SE_{Y_{\theta }}$ 在点 $X^{(k)}$ 处的梯度，H∈Rd×d为Hessian矩阵，η > 0为固定步长。 $H^{-1}$ 是Hessian矩阵的逆，它是可逆的，因为我们选择了ω s.t。|ω| > d。

梯度计算。对于每个θ∈ω，计算梯度 $\triangledown SE_{Y_{\theta }}(X^{(k)})$ 要求对于每个j∈J，计算最优的1-D分配 $\sigma_{\theta }^{j}$

通过对 $X_{\theta }^{(k)}$ 和 $Y_{\theta }^{j}$ 的值进行排序，需要 $O(N\log (N))$ 次操作来详细计算(5)。

梯度的每个元素i可以表示为

Hessian矩阵

观察到Hessian矩阵与点X无关，因此可以预先计算。

3.3 计算分布上的投影

在许多应用中，人们不仅对计算Wasserstein距离 $W_{2 }(X,Y)$ 感兴趣，而且对最优排列 $\sigma ^{\star }$ 使 $W_{\sigma }(X,Y)$ in(2)最小化。这种最优排列允许计算正交投影

其中，与Y表示相同统计分布的所有点云的集合[Y]定义于(1)。

切片投影。还可以观察到，在初始化 $X^{(0)}:=Y^{1}$ ，权值为{ρ1 = 0， ρ2 = 1}的两个分布(|J| = 2)的特殊情况下，我们的算法类似于[29]中提出的计算点云 $Y^{1}$ 在 $Y^{2}$ 上的近似投影(即近似最优分配)的原始算法。唯一不同的是，如[29]所述，我们使用了随机梯度下降(即在每次迭代k时使用随机的方向集 $\omega _{k}\in \Omega$ )。需要指出的是，当考虑到我们的质心能量时，这种梯度下降策略并不稳定。