ECCV-2018《Variational Wasserstein Clustering》

核心思想

该论文提出了一个基于最优传输(optimal transportation) 理论的新型聚类方法，称为变分Wasserstein聚类(Variational Wasserstein Clustering, VWC)。其核心思想有三点：

1. 建立最优传输与k-means聚类的联系：作者指出k-means聚类问题本质上等价于求解一个特殊的Wasserstein重心问题(Wasserstein barycenter problem)，当目标是一个单变量测度时，这被称为Wasserstein均值问题(Wasserstein means problem)。

2. 采用Monge-Brenier最优传输视角：与主流的Kantorovich最优传输方法不同，本文采用Monge-Brenier方法，将最优传输视为测度保持映射(measure-preserving mapping)，即一个样本不能被分割到多个位置，这与k-means聚类的特性更加吻合。

3. 利用power diagrams作为传输计划：通过变分原理(variational principle)求解最优传输问题，使用power Voronoi图作为传输计划，将任意域聚集成固定数量的簇，同时在移动簇中心点的过程中保持最小聚类能量。

这种方法的优势在于：(1)它是局部微分同胚；(2)不需要预先计算成对距离；(3)避免在乘积空间中搜索，从而大幅减少参数数量。

目标函数

论文的目标函数基于2-Wasserstein距离，旨在找到一组稀疏的簇中心$Y(y,\nu )$，使其与目标分布$X(x,\mu )$之间的Wasserstein距离最小化：

$$\underset{Y\in P(M)}{\inf }{W}_2^2(X,Y)=\underset{Y\in P(M),\pi \in P(M\times M)}{\inf }\sum _{y_j=\pi (x_i)}{\mu }_i\Vert x_i-y_j{\Vert }^2$$

其中$P(M)$表示度量空间$M$上的所有Borel概率测度集合。

在固定测度$\nu$的情况下，该问题等价于：

$$f(h,y)=\sum _{j=1}^k\sum _{x_i\in V_j(h)}{\mu }_i\Vert x_i-y_j{\Vert }^2$$

这里：

• $h=(h_1,\dots ,h_k{)}^T$是power diagram的参数向量
• $y=(y_1,\dots ,y_k)$是簇中心
• $V_j(h)$是由power diagram定义的第$j$个Voronoi单元
• ${\mu }_i$是样本$x_i$的测度

power Voronoi图的定义为：
Vj≜{m∈M∣∥m−yj∥2−rj2≤∥m−yi∥2−ri2},∀j≠iV_j \triangleq \{m \in M \mid \|m - y_j\|^2 - r_j^2 \leq \|m - y_i\|^2 - r_i^2\}, \forall j \neq iVj​≜{m∈M∣∥m−yj​∥2−rj2​≤∥m−yi​∥2−ri2​},∀j=i

其中$r_j$与单元的总质量相关。

目标函数的优化过程

论文提出了迭代测度保持映射(Iterative Measure-Preserving Mapping) 算法来优化目标函数，该算法交替执行两个步骤：

1. 固定簇中心$y$，更新power diagram（即更新$h$）

通过求解变分最优传输问题，最小化能量函数：
$$E(h)={\int }_{\Omega }{\theta }_h(x)\mu (x)dx-\sum _{j=1}^k{\nu }_j{h}_j$$

其中${\theta }_h(x)=\max \{⟨x,y_j⟩+h_j\}$是一个分段线性凸函数。

使用牛顿法求解：

• 梯度：$\nabla E(h)=(w_1(h)-{\nu }_1,\dots ,w_k(h)-{\nu }_k{)}^T$
• Hessian矩阵：
  $$H=\frac{{\partial }^{2}E(h)}{\partial h_i\partial h_j}=\{\begin{array}{ll}{\sum }_l\frac{{\int }_{f_{il}}\mu (x)dx}{\Vert y_l-y_i\Vert },& i=j,\forall l,\text{s.t. }{f}_{il}\ne \varnothing \\ -\frac{{\int }_{f_{ij}}\mu (x)dx}{\Vert y_j-y_i\Vert },& i\ne j,f_{ij}\ne \varnothing \\ 0,& i\ne j,f_{ij}=\varnothing \end{array}$$

其中$w_j(h)={\sum }_{x\in V_j(h)}\mu (x)$是第$j$个单元的总质量，$f_{ij}$是相邻单元的交集。

更新规则：$h^{(t+1)}\leftarrow h^{(t)}-\lambda H^{-1}\nabla E(h)$

2. 固定power diagram（即固定$h$），更新簇中心$y$

更新规则为加权平均：
$$y_j^{(t+1)}=\frac{{\sum }_{x\in V_j}{\mu }_i{x}_i^{(t)}}{{\sum }_{x\in V_j}{\mu }_i}$$

论文证明了该算法具有以下性质：

• 单调收敛性：每次迭代都减少目标函数值
• 有限步收敛：在有限次迭代后收敛
• 唯一解：产生唯一的局部解

主要贡献点

1. 理论贡献：建立了最优传输与k-means聚类之间的理论联系，从Monge-Brenier角度重新解释了k-means聚类问题。

2. 算法创新：提出了变分Wasserstein聚类算法，利用power Voronoi图同时优化Wasserstein距离和聚类质量，实现了测度保持映射。

3. 计算效率：相比基于Kantorovich最优传输的方法，避免了在乘积空间中搜索，大大减少了参数数量和计算复杂度。

4. 多领域应用：成功将该方法应用于域适应、重新网格化和学习表示三个不同领域，展示了其广泛适用性。

5. 理论保证：证明了算法的收敛性和解的唯一性，为方法的可靠性提供了理论支持。

实验结果

论文在三个不同任务上验证了VWC方法的有效性：

1. 合成数据上的域适应

• 实验设置：源域包含两个高斯分布（各30个样本），目标域包含两个不同均值和方差的高斯分布（各1500个样本）
• 结果：
  • VWC在RBF核上的分类准确率达到99.31%，略高于D2(99.25%)和JDOT(99.23%)
  • 传统k-means++在没有先验知识（如线性偏移）的情况下表现极差（准确率仅50%）
  • VWC不需要预先知道两个域之间的关系，即可有效进行知识迁移

2. 三角网格变形

• 实验设置：将人脸表面的三角网格重新分布，使顶点向高曲率区域移动
• 结果：
  • 鼻子尖端等高曲率区域顶点更密集，额头等平坦区域顶点更稀疏
  • 通过将曲面映射到单位圆盘，应用VWC，再映射回原始曲面，实现了基于曲率的自适应网格划分
  • 该方法可用于计算机图形学中的网格优化

3. 脑图像表示学习

• 实验设置：在100个MRI脑图像上进行实验（50个阿尔茨海默病AD，50个正常对照NC）
• 结果：
  • 使用VWC学习的低维表示在SVM分类中表现显著优于PCA
  • 随着中心点数量增加，VWC的分类准确率稳定提高，而PCA表现波动
  • 验证了VWC在医学图像分析中的潜力，特别是对AD相关脑萎缩的表征能力

算法实现过程

VWC的完整算法流程如下：

1. 预处理阶段

• 初始化：从已知分布中采样得到初始测度$\nu$
• 域统一：如果源域$M$和目标域$N$不同，使用调和映射(Harmonic mapping)将它们映射到凸规范空间（通常是欧氏空间$R^n$或单位圆盘$D^n$）
  • 例如：将3D脑表面映射到单位球

2. 迭代测度保持映射

(a) 更新Voronoi划分（变分最优传输）

• 计算当前$h$对应的power diagram $V$
• 计算每个单元的权重w(h)={∑m∈Vjμ(m)}w(h) = \{\sum_{m \in V_j} \mu(m)\}w(h)={∑m∈Vj​​μ(m)}
• 计算能量函数的梯度$\nabla E(h)$和Hessian矩阵$H$
• 更新$h\leftarrow h-\lambda H^{-1}\nabla E(h)$
• 重复直到$\Vert \nabla E(h)\Vert <ϵ$

(b) 更新簇中心

• 对每个簇$j$，计算新的中心点：
  $$y_j=\frac{{\sum }_{x\in V_j}{\mu }_i{x}_i}{{\sum }_{x\in V_j}{\mu }_i}$$
• 这是基于测度${\mu }_i$的加权平均

© 收敛判断

• 重复(a)和(b)步骤，直到簇中心$y$的变化小于阈值

3. 输出结果

• 返回测度保持映射$\pi :X\to Y$，表示为$(y,V)$
• 其中$y$是最终的簇中心，$V$是对应的power Voronoi图

实现细节

• 使用Voro++库计算Voronoi图
• 对于高维数据，可能面临计算和内存挑战
• 可通过梯度下降替代完整Voronoi图计算，只需相邻单元交集来计算Hessian
• 每个样本的分配可通过最近邻搜索确定

总结

VWC方法将最优传输理论与k-means聚类巧妙结合，提供了一种全新的聚类视角。其核心优势在于能够同时优化聚类质量和Wasserstein距离，实现测度保持的映射。该方法在域适应、网格优化和医学图像分析等任务中展示了出色的性能，特别是对于需要保持分布特性的应用场景。相比传统聚类方法和基于Kantorovich最优传输的方法，VWC在计算效率和理论保证方面都有明显优势，为聚类分析和分布匹配问题提供了新的解决方案。

以下是对该公式的详细推导与解释：

---

1. 幂图（Power Diagram）的定义

给定一组中心点 {y1,y2,…,yk}⊂Rn\{y_1, y_2, \dots, y_k\} \subset \mathbb{R}^n{y1​,y2​,…,yk​}⊂Rn 和权重 {r12,r22,…,rk2}\{r_1^2, r_2^2, \dots, r_k^2\}{r12​,r22​,…,rk2​}，幂图（Power Diagram）将空间 ${\mathbb{R}}^n$ 划分为 $k$ 个区域 $V_1,V_2,\dots ,V_k$，其中：
$$V_j=\left\{m\in M|\Vert m-y_j{\Vert }^2-r_j^2\le \Vert m-y_i{\Vert }^2-r_i^2,\forall i\ne j\right\}.$$
每个区域 $V_j$ 包含所有点 $m$，使得对任意其他中心 $y_i$，加权距离 $\Vert m-y_j{\Vert }^2-r_j^2$ 不大于 $\Vert m-y_i{\Vert }^2-r_i^2$。

---

2. 不等式的推导

从幂图的定义出发，我们推导出更简洁的线性不等式形式：

步骤1：展开平方项

对任意 $i\ne j$，比较 $\Vert m-y_j{\Vert }^2-r_j^2$ 和 $\Vert m-y_i{\Vert }^2-r_i^2$：
$$\Vert m-y_j{\Vert }^2-r_j^2\le \Vert m-y_i{\Vert }^2-r_i^2.$$
展开平方项：
$$(m-y_j{)}^T(m-y_j)-r_j^2\le (m-y_i{)}^T(m-y_i)-r_i^2.$$

步骤2：化简表达式

展开后得到：
$$m^{T}m-2m^T{y}_j+y_j^T{y}_j-r_j^2\le m^{T}m-2m^T{y}_i+y_i^T{y}_i-r_i^2.$$
两边同时减去 $m^{T}m$，并整理：
$$-2m^T{y}_j+(y_j^T{y}_j-r_j^2)\le -2m^T{y}_i+(y_i^T{y}_i-r_i^2).$$

步骤3：移项得到线性不等式

将含 $m$ 的项移到左边，常数项移到右边：
$$-2m^T{y}_j+2m^T{y}_i\le (y_i^T{y}_i-r_i^2)-(y_j^T{y}_j-r_j^2).$$
提取公因子 $2$：
$$2m^T(y_i-y_j)\le (y_i^T{y}_i-r_i^2)-(y_j^T{y}_j-r_j^2).$$
两边除以 $2$，得到最终形式：
$$m^T{y}_j-\frac{1}{2}(y_j^T{y}_j+r_j^2)\le m^T{y}_i-\frac{1}{2}(y_i^T{y}_i+r_i^2).$$

---

3. 几何意义

上述不等式可以看作是超平面分割的条件：

• 左端：$m^T{y}_j-\frac{1}{2}(y_j^T{y}_j+r_j^2)$ 是关于 $m$ 的线性函数。
• 右端：$m^T{y}_i-\frac{1}{2}(y_i^T{y}_i+r_i^2)$ 同样是关于 $m$ 的线性函数。

因此，每个不等式对应一个超平面：
$$m^T(y_j-y_i)\le \frac{1}{2}\left(y_i^T{y}_i-y_j^T{y}_j+r_j^2-r_i^2\right).$$
这表明，幂图的每个区域 $V_j$ 是由多个超平面分割出的凸多面体，因此整个幂图是凸分割。

---

4. 与分段线性凸函数的关系

根据文献[27]（Aurenhammer, 1987），幂图与分段线性凸函数存在一一对应关系：

• 每个幂图区域 $V_j$ 对应函数 $u_h(x)$ 的一个线性片（affine piece）。
• 函数 $u_h(x)$ 在区域 $V_j$ 内的表达式为：
  $$u_h(x)=x^T{y}_j-\frac{1}{2}(y_j^T{y}_j+r_j^2).$$
• 因此，幂图是该分段线性凸函数的次微分分解（subdifferential decomposition）。

---

5. 应用场景

在论文《Variational Wasserstein Clustering》中，幂图被用于：

1. 最优传输映射：通过调整权重 $r_j^2$，将源分布 $X$ 映射到目标分布 $Y$。
2. 聚类优化：每个幂图区域 $V_j$ 对应一个簇，簇中心 $y_j$ 和权重 $r_j^2$ 通过迭代优化确定，以最小化Wasserstein距离。

---

总结

该不等式揭示了幂图的线性超平面结构，并通过分段线性凸函数建立了理论联系。这一结果为变分Wasserstein聚类算法提供了关键工具，使得复杂的最优传输问题可以通过凸优化求解。

---

1. 核心概念解析

(a) 概率测度空间 $\mathcal{P}(M)$

• 定义：$\mathcal{P}(M)$ 是定义在度量空间 $M$ 上的所有Borel概率测度的集合。
• 示例：
  • 离散测度：如 $X={\sum }_{i=1}^n{\mu }_i{\delta }_{x_i}$（样本点 $x_i$ 的加权组合）。
  • 连续测度：如高斯分布 $N(\mu ,{\sigma }^2)$。

(b) 测度保持映射 $T$

• 定义：若映射 $T:X\to Y$ 满足：
  $$\mu (T^{-1}(B))=\nu (B),\forall B\subset Y,$$
  则称 $T$ 是测度保持映射（measure-preserving mapping）。
• 几何意义：将 $X$ 的质量完美地“搬运”到 $Y$，不增不减。

© 耦合（Coupling）

• 定义：耦合 $\pi$ 是 $X$ 和 $Y$ 在乘积空间 $M\times M$ 上的联合概率测度，其边缘分布为：
  $$\mu =\pi (\cdot ,M),\nu =\pi (M,\cdot ).$$
• 物理意义：$\pi (x,y)$ 表示将质量从 $x$ 运输到 $y$ 的比例。

---

2. 最优传输问题

(a) 运输成本函数 $c(x,y)$

• 定义：通常取为测地距离的 $p$ 次幂：
  $$c(x,y)=d(x,y{)}^p.$$
• 常见选择：
  • $p=1$：对应Earth Mover’s Distance (EMD)。
  • $p=2$：对应最常见的Wasserstein距离。

(b) Wasserstein距离的定义

$$W_p(\mu ,\nu )={\left(\underset{\pi \in \Pi (\mu ,\nu )}{\inf }{\int }_{M\times M}c(x,y)d\pi (x,y)\right)}^{1/p},$$
其中 $\Pi (\mu ,\nu )$ 是所有满足边缘约束的耦合集合。

• 物理意义：$W_p$ 是将分布 $\mu$ 转换为 $\nu$ 的最小平均运输成本。
• 关键性质：
  • 满足距离公理（非负性、同一性、对称性、三角不等式）。
  • 对分布的形状和位置敏感，适合衡量高维数据差异。

---

3. Monge型与Kantorovich型最优传输

(a) Monge型传输

• 限制条件：每个质量点 $x$ 只能被映射到一个 $y$，即：
  $$d\pi (x,y)=d\mu (x)\delta [y=T(x)].$$
• 目标函数：
  $${\pi }_{T_{\text{opt}}}=T_{\text{opt}}=\underset{T}{\mathrm{argmin}}{\int }_{M}c(x,T(x))d\mu (x).$$
• 优点：直接给出显式的映射 $T$，便于几何解释。

(b) Kantorovich型传输

• 允许质量分裂：一个质量点 $x$ 可以被分配到多个 $y$，通过耦合 $\pi (x,y)$ 描述。
• 目标函数：
  $$W_p(\mu ,\nu )={\left(\underset{\pi \in \Pi (\mu ,\nu )}{\inf }{\int }_{M\times M}c(x,y)d\pi (x,y)\right)}^{1/p}.$$

---

4. 式(2)的推导

当研究Monge型传输时，由于每个 $x$ 只能映射到一个 $y$，耦合 $\pi$ 可表示为：
$$\pi (x,y)=\delta [y=T(x)]d\mu (x).$$
代入Wasserstein距离的定义：
$$W_p^p(\mu ,\nu )={\int }_{M}c(x,T(x))d\mu (x).$$
因此，寻找最优传输等价于：
$$T_{\text{opt}}=\underset{T}{\mathrm{argmin}}{\int }_{M}c(x,T(x))d\mu (x),$$
即式(2)所示。

---

5. 应用场景与意义

在《Variational Wasserstein Clustering》中，作者利用Monge型传输的特性：

1. 保持测度：确保聚类过程中质量守恒。
2. 显式映射：通过求解 $T_{\text{opt}}$ 直接获得簇中心与样本的对应关系。
3. 变分原理：将最优传输转化为能量最小化问题，便于数值求解。

---

总结

该段文字的核心是：

• 定义了最优传输问题，通过耦合 $\pi$ 描述质量转移。
• 区分了Monge型与Kantorovich型传输，前者更适合聚类任务。
• 推导了Wasserstein距离的表达式，并指出在Monge框架下可简化为式(2)。

这一理论基础为后续的变分Wasserstein聚类算法提供了支撑，使得聚类过程既保持测度又具有几何可解释性。

Wasserstein度量的数学原理

Wasserstein度量（也称为Earth Mover’s Distance，EMD）是基于最优传输理论的概率分布间距离度量。下面我将系统阐述其数学原理，从基础定义到高级理论。

1. 最优传输问题基础

1.1 Monge问题（1781年）

Monge提出的原始问题：给定两个概率测度 $\mu$ 和 $\nu$，寻找一个测度保持映射 $T:X\to Y$，使得传输成本最小化：

$$\underset{T_{\#}\mu =\nu }{\inf }{\int }_{X}c(x,T(x))d\mu (x)$$

其中：

• $T_{\#}\mu =\nu$ 表示 $T$ 是测度保持映射（即 $\mu (T^{-1}(B))=\nu (B),\forall B\subset Y$）
• $c(x,y)$ 是传输成本函数（通常为 $c(x,y)=d(x,y{)}^p$）
• $d(x,y)$ 是基础空间上的距离度量

关键限制：每个质量点 $x$ 只能被映射到一个 $y$（不能分割）

1.2 Kantorovich松弛（1941年）

Kantorovich将问题松弛为：

$$\underset{\pi \in \Pi (\mu ,\nu )}{\inf }{\int }_{X\times Y}c(x,y)d\pi (x,y)$$

其中 $\Pi (\mu ,\nu )$ 是所有满足边缘约束的耦合集合：

• $\pi (\cdot ,Y)=\mu$
• $\pi (X,\cdot )=\nu$

关键区别：允许质量分割，即一个 $x$ 可以被分配到多个 $y$

2. Wasserstein距离的严格定义

2.1 p-Wasserstein距离

对于 $p\ge 1$，p-Wasserstein距离定义为：

$$W_p(\mu ,\nu )={\left(\underset{\pi \in \Pi (\mu ,\nu )}{\inf }{\int }_{M\times M}d(x,y{)}^{p}d\pi (x,y)\right)}^{1/p}$$

其中：

• $M$ 是基础度量空间
• $d(x,y)$ 是 $M$ 上的度量
• $\Pi (\mu ,\nu )$ 是所有满足边缘约束的耦合集合

2.2 Wasserstein距离的性质

Wasserstein距离满足度量空间的所有公理：

1. 非负性：$W_p(\mu ,\nu )\ge 0$
2. 同一性：$W_p(\mu ,\nu )=0$ 当且仅当 $\mu =\nu$
3. 对称性：$W_p(\mu ,\nu )=W_p(\nu ,\mu )$
4. 三角不等式：$W_p(\mu ,\zeta )\le W_p(\mu ,\nu )+W_p(\nu ,\zeta )$

此外，它还具有以下重要特性：

• 对分布形状敏感：不仅考虑分布的均值，还考虑分布的整体形状
• 几何意义明确：可以理解为将一种分布"重塑"为另一种分布所需的最小"工作量"

3. Monge-Brenier最优传输理论

3.1 Brenier定理（1987年突破）

在欧氏空间 ${\mathbb{R}}^n$ 上，若 $\mu$ 是绝对连续测度，则存在唯一的最优传输映射 $T$，且 $T$ 可表示为凸函数的梯度：

$$T=\nabla \varphi$$

其中 $\varphi$ 是 ${\mathbb{R}}^n$ 上的凸函数。

证明关键：最优传输映射 $T$ 使传输成本 $\int \Vert x-T(x){\Vert }^{2}d\mu (x)$ 最小化，且 $T$ 是某个凸函数的梯度。

3.2 变分原理与Power Diagram

根据论文中的Alexandrov定理和Brenier定理：

• 给定凸多面体 $\Omega \subset {\mathbb{R}}^n$ 和 $k$ 个不同点 $y_1,\dots ,y_k\subset {\mathbb{R}}^n$，存在唯一的向量 $h=(h_1,\dots ,h_k{)}^T$（在平移意义下），使得分段线性凸函数：
  $${\theta }_h(x)=\max \{⟨x,y_j⟩+h_j\},j=1,\dots ,k$$
  满足：
  $$\text{vol}(x\in \Omega \mid \nabla {\theta }_h(x)=y_j)={\nu }_j$$

• 梯度映射 $\nabla {\theta }_h$ 提供了Monge问题的解，即最小化传输成本 ${\int }_{\Omega }\Vert x-{\theta }_h(x){\Vert }^2$

• Power Voronoi图：由凸函数 ${\theta }_h$ 诱导的凸分割，定义为：
  Vj={m∈M∣∥m−yj∥2−rj2≤∥m−yi∥2−ri2,∀i≠j}V_j = \{m \in M \mid \|m - y_j\|^2 - r_j^2 \leq \|m - y_i\|^2 - r_i^2, \quad \forall i \neq j\}Vj​={m∈M∣∥m−yj​∥2−rj2​≤∥m−yi​∥2−ri2​,∀i=j}
  其中 $r_j^2=-\Vert y_j{\Vert }^2-2h_j$

4. Wasserstein距离的计算方法

4.1 Kantorovich方法

• 将问题转化为线性规划
• 需要预先计算成对距离
• 在离散情况下，解耦合 $\pi$ 是 $n\times m$ 矩阵（$n$ 和 $m$ 分别是两个分布的样本数）
• 计算复杂度为 $O(n^3\log n)$

4.2 Monge-Brenier方法（变分方法）

• 通过变分原理求解：最小化能量函数
  $$E(h)={\int }_{\Omega }{\theta }_h(x)\mu (x)dx-\sum _{j=1}^k{\nu }_j{h}_j$$

• 梯度：$\nabla E(h)=(w_1(h)-{\nu }_1,\dots ,w_k(h)-{\nu }_k{)}^T$
  其中 $w_j(h)={\sum }_{x\in V_j(h)}\mu (x)$ 是第 $j$ 个单元的总质量

• Hessian矩阵：
  $$H=\frac{{\partial }^{2}E(h)}{\partial h_i\partial h_j}=\{\begin{array}{ll}{\sum }_l\frac{{\int }_{f_{il}}\mu (x)dx}{\Vert y_l-y_i\Vert },& i=j,\forall l,\text{s.t. }{f}_{il}\ne \varnothing \\ -\frac{{\int }_{f_{ij}}\mu (x)dx}{\Vert y_j-y_i\Vert },& i\ne j,f_{ij}\ne \varnothing \\ 0,& i\ne j,f_{ij}=\varnothing \end{array}$$
  其中 $f_{ij}$ 是相邻单元的交集

• 通过牛顿法迭代求解：$h^{(t+1)}\leftarrow h^{(t)}-\lambda H^{-1}\nabla E(h)$

5. Wasserstein距离与k-means聚类的联系

5.1 Wasserstein均值问题

当目标分布 $Y$ 是稀疏的（即 $Y$ 有有限支持），Wasserstein距离最小化问题退化为：

$$\underset{Y\in P(M)}{\inf }{W}_2^2(X,Y)=\underset{Y\in P(M),\pi \in P(M\times M)}{\inf }\sum _{y_j=\pi (x_i)}{\mu }_i\Vert x_i-y_j{\Vert }^2$$

这正是k-means聚类问题：

• $Y$ 是簇中心集合
• $\pi$ 确定了样本到簇的分配
• 目标是最小化加权平方距离和

5.2 变分Wasserstein聚类

如论文所述，变分Wasserstein聚类通过以下方式同时优化：

1. 聚类质量：${\sum }_{x_i\in V_j}{\mu }_i\Vert x_i-y_j{\Vert }^2$
2. Wasserstein距离：$W_2^2(X,Y)$

通过迭代更新：

• Power Voronoi图：固定簇中心 $y$，更新划分 $V$
• 簇中心：固定划分 $V$，更新 $y_j=\frac{{\sum }_{x\in V_j}{\mu }_i{x}_i}{{\sum }_{x\in V_j}{\mu }_i}$

6. Wasserstein度量的几何解释

6.1 测地距离视角

在概率测度空间 $({\mathcal{P}}_2(M),W_2)$ 上：

• $W_2$ 是测地距离
• 概率测度间的最短路径由McCann插值给出：
  $${\mu }_t=((1-t)I+tT{)}_{\#}{\mu }_0,t\in [0,1]$$
  其中 $T$ 是从 ${\mu }_0$ 到 ${\mu }_1$ 的最优传输映射

6.2 Barycenter（重心）

多个分布 ${\mu }_1,\dots ,{\mu }_n$ 的Wasserstein重心定义为：
$$\bar{\mu }=\arg \underset{\mu \in {\mathcal{P}}_2(M)}{\min }\sum _{i=1}^n{\lambda }_i{W}_2^2(\mu ,{\mu }_i)$$

当 $n=1$ 时，这就是Wasserstein均值问题，与k-means聚类等价。

7. 不同p值的Wasserstein距离

7.1 1-Wasserstein距离 (EMD)

$$W_1(\mu ,\nu )=\underset{\pi \in \Pi (\mu ,\nu )}{\inf }{\int }_{M\times M}d(x,y)d\pi (x,y)$$

• 优点：计算相对简单，有双对偶形式
• 应用：图像和形状比较

7.2 2-Wasserstein距离

$$W_2(\mu ,\nu )={\left(\underset{\pi \in \Pi (\mu ,\nu )}{\inf }{\int }_{M\times M}d(x,y{)}^{2}d\pi (x,y)\right)}^{1/2}$$

• 优点：
  • 具有良好的几何性质
  • 与欧氏距离有紧密联系
  • 可通过变分方法高效计算
• 应用：图像处理、几何处理、聚类

8. Wasserstein度量的优势与局限

8.1 优势

1. 对分布敏感：考虑分布的整体形状，而非仅统计矩
2. 几何意义明确：具有清晰的物理和几何解释
3. 连续性：当分布连续变化时，Wasserstein距离也连续变化
4. 适用于不同支撑集：即使两个分布的支撑集不相交，也能定义有意义的距离

8.2 局限

1. 计算复杂度高：尤其对于高维数据
2. 维度灾难：在高维空间中收敛速度较慢
3. 对离群点敏感：极端值可能显著影响距离值

9. 总结

Wasserstein度量基于最优传输理论，为概率分布间提供了一种具有几何意义的距离度量。其核心数学原理包括：

1. 最优传输问题：Monge问题和Kantorovich松弛
2. Brenier定理：将最优传输映射表示为凸函数的梯度
3. 变分原理：通过能量最小化求解最优传输
4. Power Voronoi图：作为传输计划的几何表示
5. 与k-means聚类的联系：Wasserstein均值问题等价于加权k-means聚类

Wasserstein度量不仅在理论上具有优雅的数学结构，在实际应用中也展现出强大能力，如论文中所示的域适应、网格变形和表示学习等任务。特别是2-Wasserstein距离，因其良好的几何性质和与聚类问题的自然联系，在机器学习和计算机视觉领域得到广泛应用。

---

如何更新划分 V：变分Wasserstein聚类中的关键步骤

在变分Wasserstein聚类(Variational Wasserstein Clustering, VWC)中，更新划分 $V$ 是算法的核心环节。下面我将详细解释这一过程的数学原理和实现细节。

1. 划分 $V$ 的数学定义

划分 $V$ 是由power Voronoi diagram定义的，其数学表达为：

Vj={m∈M∣∥m−yj∥2−rj2≤∥m−yi∥2−ri2},∀j≠iV_j = \{m \in M \mid \|m - y_j\|^2 - r_j^2 \leq \|m - y_i\|^2 - r_i^2\}, \quad \forall j \neq iVj​={m∈M∣∥m−yj​∥2−rj2​≤∥m−yi​∥2−ri2​},∀j=i

通过简单代数变换，可得等价形式：

$$m^T{y}_j-\frac{1}{2}(y_j^T{y}_j+r_j^2)\le m^T{y}_i-\frac{1}{2}(y_i^T{y}_i+r_i^2)$$

其中 $r_j^2$ 与单元总质量相关，且 $h_j=-\frac{1}{2}(\Vert y_j{\Vert }^2+r_j^2)$。

2. 更新划分 $V$ 的算法流程

更新划分 $V$ 的完整过程由Algorithm 1 (Variational-OT)实现，核心步骤如下：

(a) 初始化

• 设置初始参数 $h^{(0)}=0$

(b) 迭代更新 $h$

重复以下步骤直到收敛（$\Vert \nabla E(h)\Vert <ϵ$）：

1. 更新power diagram $V$：

   • 使用当前的 $(y,h)$ 计算power Voronoi图
   • 每个单元 $V_j(h)$ 由不等式 $m^T{y}_j+h_j\ge m^T{y}_i+h_i$ 定义

2. 计算单元权重：
   $$w_j(h)=\sum _{m\in V_j}\mu (m)$$

   • $w_j(h)$ 表示第 $j$ 个Voronoi单元的总质量

3. 计算梯度和Hessian：

   • 梯度：$\nabla E(h)=(w_1(h)-{\nu }_1,\dots ,w_k(h)-{\nu }_k{)}^T$
   • Hessian矩阵：
     $$H=\frac{{\partial }^{2}E(h)}{\partial h_i\partial h_j}=\{\begin{array}{ll}{\sum }_l\frac{{\int }_{f_{il}}\mu (x)dx}{\Vert y_l-y_i\Vert },& i=j,\forall l,\text{s.t. }{f}_{il}\ne \varnothing \\ -\frac{{\int }_{f_{ij}}\mu (x)dx}{\Vert y_j-y_i\Vert },& i\ne j,f_{ij}\ne \varnothing \\ 0,& i\ne j,f_{ij}=\varnothing \end{array}$$
     其中 $f_{ij}$ 是相邻单元 $V_i$ 和 $V_j$ 的交集

4. 牛顿法更新 $h$：
   $$h^{(t+1)}\leftarrow h^{(t)}-\lambda H^{-1}\nabla E(h)$$

   • $\lambda$ 是步长参数（通常设为1）
   • $H^{-1}$ 是Hessian矩阵的逆

© 返回结果

• 最终的划分 $V$ 和参数 $h$

3. 数学原理详解

(a) 能量函数 $E(h)$

更新划分 $V$ 的核心是优化能量函数：
$$E(h)={\int }_{\Omega }{\theta }_h(x)\mu (x)dx-\sum _{j=1}^k{\nu }_j{h}_j$$
其中 ${\theta }_h(x)=\max \{⟨x,y_j⟩+h_j\}$ 是分段线性凸函数。

物理意义：

• ${\int }_{\Omega }{\theta }_h(x)\mu (x)dx$：表示传输成本
• ${\sum }_{j=1}^k{\nu }_j{h}_j$：表示约束项
• 最小化 $E(h)$ 等价于找到最优传输映射

(b) 梯度的几何解释

$$\nabla E(h)=(w_1(h)-{\nu }_1,\dots ,w_k(h)-{\nu }_k{)}^T$$

• 当 $w_j(h)>{\nu }_j$：第 $j$ 个单元质量过大，需要缩小
• 当 $w_j(h)<{\nu }_j$：第 $j$ 个单元质量过小，需要扩大

通过调整 $h_j$，可以控制单元 $V_j$ 的大小：

• 增加 $h_j$：扩大 $V_j$
• 减少 $h_j$：缩小 $V_j$

© Hessian的几何解释

Hessian矩阵 $H$ 描述了单元边界的变化率：

• 对角线元素：${\sum }_l\frac{{\int }_{f_{il}}\mu (x)dx}{\Vert y_l-y_i\Vert }$

  • 表示单元 $V_i$ 的"刚度"，值越大越难变形
  • 与相邻边界的长度和质量成正比

• 非对角线元素：$-\frac{{\int }_{f_{ij}}\mu (x)dx}{\Vert y_j-y_i\Vert }$

  • 表示单元 $V_i$ 和 $V_j$ 之间的"耦合强度"
  • 值越大表示边界越容易移动

4. 实现细节

(a) Voronoi图计算

• 使用Voro++库计算power Voronoi图
• 对于每个样本点 $x_i$，确定其所属的Voronoi单元 $V_j$：
  j=arg⁡max⁡k{xiTyk+hk}j = \arg\max_k \{x_i^T y_k + h_k\}j=argkmax​{xiT​yk​+hk​}

(b) 边界积分计算

• 计算相邻单元交集 $f_{ij}$ 的质量：
  $${\int }_{f_{ij}}\mu (x)dx=\sum _{x\in f_{ij}}\mu (x)$$
• 在离散情况下，$f_{ij}$ 由共享边界的样本点组成

© 收敛条件

• 当 $\Vert \nabla E(h)\Vert <ϵ$ 时停止迭代
• 通常 $ϵ=10^{-6}$

5. 与聚类过程的整合

在完整的VWC算法中，更新划分 $V$ 是迭代测度保持映射的一部分：

1. 更新划分 $V$：

   • 通过Variational-OT算法，固定簇中心 $y$，更新划分 $V$
   • 目标：使每个单元的总质量 $w_j(h)$ 等于目标测度 ${\nu }_j$

2. 更新簇中心 $y$：

   • 固定划分 $V$，更新簇中心：
     $$y_j^{(t+1)}=\frac{{\sum }_{x\in V_j}{\mu }_i{x}_i}{{\sum }_{x\in V_j}{\mu }_i}$$
   • 这是基于测度 ${\mu }_i$ 的加权平均

这两个步骤交替进行，直到收敛，确保：

• 划分 $V$ 满足测度保持条件
• 簇中心 $y$ 最小化聚类能量

6. 几何直观解释

考虑二维空间中的一个简单例子：

• 有3个簇中心 $y_1,y_2,y_3$
• 初始 $h=(0,0,0{)}^T$ 产生标准Voronoi图
• 如果单元 $V_1$ 质量过大（$w_1>{\nu }_1$）
  • 减小 $h_1$，使 $V_1$ 缩小
  • 边界向 $y_1$ 移动
• 如果单元 $V_2$ 质量过小（$w_2<{\nu }_2$）
  • 增加 $h_2$，使 $V_2$ 扩大
  • 边界远离 $y_2$

这个过程持续进行，直到所有单元的质量匹配目标测度 $\nu$。

7. 总结

更新划分 $V$ 的核心是通过调整参数 $h$ 优化power Voronoi图，使其满足测度保持条件。这一过程：

1. 基于变分原理：最小化能量函数 $E(h)$
2. 使用牛顿法：高效求解凸优化问题
3. 保持几何结构：通过power diagram维持凸分割
4. 确保测度守恒：使每个单元的总质量等于目标测度

这种方法的优势在于：

• 避免了在乘积空间中搜索（相比Kantorovich方法）
• 保证了局部微分同胚性质
• 同时优化Wasserstein距离和聚类质量

在实际应用中，这一过程高效稳定，为域适应、网格变形和表示学习等任务提供了强大的理论和算法支持。

---

如何理解 Wasserstein 重心（Barycenter）

1. Wasserstein 重心的定义

多个概率测度 ${\mu }_1,\dots ,{\mu }_n$ 的 Wasserstein 重心（Barycenter）定义为：
$$\bar{\mu }=\arg \underset{\mu \in {\mathcal{P}}_2(M)}{\min }\sum _{i=1}^n{\lambda }_i{W}_2^2(\mu ,{\mu }_i),$$
其中：

• $W_2(\cdot ,\cdot )$ 是 2-Wasserstein 距离
• ${\lambda }_i\ge 0$ 是权重，满足 ${\sum }_{i=1}^n{\lambda }_i=1$
• ${\mathcal{P}}_2(M)$ 是所有具有有限二阶矩的概率测度集合

物理意义：Wasserstein 重心是使总加权传输成本最小化的“平均”分布，类似于几何中的质心概念，但扩展到概率测度空间。

---

2. 当 $n=1$ 时的特殊情形

当 $n=1$ 时，目标函数退化为：
$$\bar{\mu }=\arg \underset{\mu \in {\mathcal{P}}_2(M)}{\min }{\lambda }_1{W}_2^2(\mu ,{\mu }_1).$$
由于 ${\lambda }_1=1$（归一化条件），此问题等价于：
$$\bar{\mu }=\arg \underset{\mu \in {\mathcal{P}}_2(M)}{\min }{W}_2^2(\mu ,{\mu }_1).$$

关键结论

此时，Wasserstein 重心问题等价于Wasserstein 均值问题，其解是 ${\mu }_1$ 本身（因为最小距离在 $\mu ={\mu }_1$ 时取得）。然而，这一结论需要结合具体应用场景理解：

• 离散测度情况：若 ${\mu }_1$ 是经验分布（如 k-means 中的样本分布），则 Wasserstein 均值问题转化为寻找一个离散测度 $\bar{\mu }$（通常为单点 Dirac 测度），使其到 ${\mu }_1$ 的 Wasserstein 距离最小。
• k-means 聚类的联系：在 k-means 中，簇中心是样本的几何均值，而 Wasserstein 均值问题在离散情况下也寻求类似的结果。因此，当 $n=1$ 时，两者在数学形式上等价。

---

3. 与 k-means 聚类的等价性

3.1 k-means 目标函数

k-means 的目标是将样本划分为 $k$ 个簇，并最小化总平方误差：
$$\underset{C_1,\dots ,C_k}{\min }\sum _{j=1}^k\sum _{x\in C_j}\Vert x-m_j{\Vert }^2,$$
其中 $m_j$ 是第 $j$ 个簇的均值。

3.2 Wasserstein 均值问题

当 $n=1$ 且 ${\mu }_1$ 是经验分布（如样本分布）时，Wasserstein 均值问题可表示为：
$$\bar{\mu }=\arg \underset{\mu \in {\mathcal{P}}_2(M)}{\min }{W}_2^2(\mu ,{\mu }_1).$$
若 $\mu$ 是单点 Dirac 测度 ${\delta }_m$，则：
$$W_2^2({\delta }_m,{\mu }_1)={\int }_M\Vert x-m{\Vert }^{2}d{\mu }_1(x),$$
这正是 k-means 中簇中心 $m$ 的目标函数！

因此，当 $n=1$ 时，Wasserstein 均值问题等价于 k-means 聚类的单簇情况。

---

4. 数学推导

4.1 Wasserstein 距离的表达式

对于离散测度 ${\mu }_1=\frac{1}{N}{\sum }_{i=1}^N{\delta }_{x_i}$ 和单点测度 $\mu ={\delta }_m$，有：
$$W_2^2(\mu ,{\mu }_1)=\frac{1}{N}\sum _{i=1}^N\Vert x_i-m{\Vert }^2.$$
这与 k-means 的目标函数完全一致。

4.2 最优解

k-means 的最优解是样本均值：
$$m^{\ast }=\frac{1}{N}\sum _{i=1}^N{x}_i.$$
同样，Wasserstein 均值问题的最优解也是这个均值，因为：
$$\frac{\partial }{\partial m}\left(\frac{1}{N}\sum _{i=1}^N\Vert x_i-m{\Vert }^2\right)=-\frac{2}{N}\sum _{i=1}^N(m-x_i)=0\Rightarrow m^{\ast }=\frac{1}{N}\sum _{i=1}^N{x}_i.$$

---

5. 应用场景与意义

5.1 离散数据聚类

在 k-means 中，样本被分配到最近的簇中心，而 Wasserstein 均值问题通过最优传输理论给出了这一过程的理论依据。

5.2 概率分布的平均

Wasserstein 重心不仅适用于离散数据，还可处理连续分布。例如，在图像生成任务中，Wasserstein GAN（W-GAN）利用 Wasserstein 距离优化生成器，其目标之一是逼近真实数据分布的重心。

5.3 变分 Wasserstein 聚类

在论文《Variational Wasserstein Clustering》中，作者通过 power diagram 和变分原理求解 Wasserstein 重心，实现了同时优化聚类质量和测度保持映射。当 $n=1$ 时，这一方法退化为标准 k-means 聚类。

---

6. 总结

• Wasserstein 重心是概率测度空间中的“平均”概念，通过最小化总传输成本定义。
• 当 $n=1$ 时，Wasserstein 均值问题等价于 k-means 聚类的单簇情况，两者在数学形式和目标上完全一致。
• 这一等价性为将最优传输理论应用于聚类问题提供了理论基础，特别是在处理概率分布数据时（如图像、文本等）。