《Power Voronoi图的数学原理》
Power Voronoi图的数学原理
Power Voronoi图(也称为加权Voronoi图或幂图)是标准Voronoi图的推广形式,在最优传输理论、计算几何和机器学习中有着重要应用。以下系统阐述其数学原理:
1. 基本定义与几何解释
1.1 标准Voronoi图回顾
给定欧氏空间 Rn\mathbb{R}^nRn 中的一组点 {y1,y2,…,yk}\{y_1, y_2, \dots, y_k\}{y1,y2,…,yk},标准Voronoi图将空间划分为 kkk 个区域:
Vj={m∈Rn∣∥m−yj∥≤∥m−yi∥,∀i≠j}V_j = \{m \in \mathbb{R}^n \mid \|m - y_j\| \leq \|m - y_i\|, \quad \forall i \neq j\}Vj={m∈Rn∣∥m−yj∥≤∥m−yi∥,∀i=j}
每个区域 VjV_jVj 包含所有离 yjy_jyj 比离其他点更近的点。
1.2 Power Voronoi图的定义
给定一组点 {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 Voronoi图定义为:
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}
几何解释:
- ∥m−yj∥2−rj2\|m - y_j\|^2 - r_j^2∥m−yj∥2−rj2 称为幂距离(power distance)
- 每个单元 VjV_jVj 由不等式定义,表示点 mmm 到 yjy_jyj 的"修正距离"最小
- 当 rj=0r_j = 0rj=0 时,退化为标准Voronoi图
1.3 等价形式
通过代数变换,可得:
∥m−yj∥2−rj2≤∥m−yi∥2−ri2\|m - y_j\|^2 - r_j^2 \leq \|m - y_i\|^2 - r_i^2∥m−yj∥2−rj2≤∥m−yi∥2−ri2
⇒mTyj−12(yjTyj+rj2)≤mTyi−12(yiTyi+ri2)\Rightarrow m^T y_j - \frac{1}{2}(y_j^T y_j + r_j^2) \leq m^T y_i - \frac{1}{2}(y_i^T y_i + r_i^2)⇒mTyj−21(yjTyj+rj2)≤mTyi−21(yiTyi+ri2)
令 hj=−12(∥yj∥2+rj2)h_j = -\frac{1}{2}(\|y_j\|^2 + r_j^2)hj=−21(∥yj∥2+rj2),则等价于:
mTyj+hj≤mTyi+hi,∀i≠jm^T y_j + h_j \leq m^T y_i + h_i, \quad \forall i \neq jmTyj+hj≤mTyi+hi,∀i=j
2. 与凸函数的联系
2.1 分段线性凸函数
定义分段线性凸函数:
θh(x)=max{⟨x,yj⟩+hj},j=1,…,k\theta_h(x) = \max\{\langle x, y_j \rangle + h_j\}, \quad j = 1, \dots, kθh(x)=max{⟨x,yj⟩+hj},j=1,…,k
关键性质:
- θh(x)\theta_h(x)θh(x) 是凸函数
- 其次微分 ∂θh(x)\partial \theta_h(x)∂θh(x) 是单值的,且 ∂θh(x)=yj\partial \theta_h(x) = y_j∂θh(x)=yj 当 x∈Vj(h)x \in V_j(h)x∈Vj(h)
- 梯度映射 ∇θh\nabla \theta_h∇θh 将每个Voronoi单元 VjV_jVj 映射到点 yjy_jyj
2.2 Alexandrov定理
论文中引用的Alexandrov定理指出:
假设 Ω\OmegaΩ 是 Rn\mathbb{R}^nRn 中具有非空内部的紧凸多面体,y1,…,yk⊂Rny_1, \dots, y_k \subset \mathbb{R}^ny1,…,yk⊂Rn 是 kkk 个不同点,ν1,…,νk>0\nu_1, \dots, \nu_k > 0ν1,…,νk>0 满足 ∑j=1kνj=vol(Ω)\sum_{j=1}^k \nu_j = \text{vol}(\Omega)∑j=1kνj=vol(Ω)。存在唯一的向量 h=(h1,…,hk)Th = (h_1, \dots, h_k)^Th=(h1,…,hk)T(在平移意义下),使得分段线性凸函数 θh(x)=max{⟨x,yj⟩+hj}\theta_h(x) = \max\{\langle x, y_j \rangle + h_j\}θh(x)=max{⟨x,yj⟩+hj} 满足:
vol(x∈Ω∣∇θh(x)=yj)=νj\text{vol}(x \in \Omega \mid \nabla \theta_h(x) = y_j) = \nu_jvol(x∈Ω∣∇θh(x)=yj)=νj
物理意义:
- 该定理确保了给定目标测度 ν\nuν,存在唯一的power Voronoi图实现测度保持映射
2.3 Brenier定理
Brenier定理建立了最优传输与凸函数的关系:
梯度映射 ∇θh\nabla \theta_h∇θh 提供了Monge问题的解,即最小化传输成本 ∫Ω∥x−θh(x)∥2\int_\Omega \|x - \theta_h(x)\|^2∫Ω∥x−θh(x)∥2
关键结论:
- 最优传输映射 TTT 可表示为凸函数的梯度:T=∇ϕT = \nabla \phiT=∇ϕ
- 在离散情况下,ϕ\phiϕ 是分段线性凸函数 θh\theta_hθh
- Power Voronoi图是该凸函数的次微分分解
3. 在最优传输中的应用
3.1 半离散最优传输
在变分Wasserstein聚类中,考虑半离散最优传输问题:
- 源分布 X(x,μ)X(x, \mu)X(x,μ):连续或离散经验分布
- 目标分布 Y(y,ν)Y(y, \nu)Y(y,ν):稀疏离散测度 {(yj,νj)}\{(y_j, \nu_j)\}{(yj,νj)}
Power Voronoi图提供了测度保持映射 π:X→Y\pi: X \to Yπ:X→Y,使得:
μ(Vj(h))=νj,j=1,…,k\mu(V_j(h)) = \nu_j, \quad j = 1, \dots, kμ(Vj(h))=νj,j=1,…,k
3.2 能量函数与优化
定义能量函数:
E(h)=∫Ωθh(x)μ(x)dx−∑j=1kνjhjE(h) = \int_\Omega \theta_h(x)\mu(x)dx - \sum_{j=1}^k \nu_j h_jE(h)=∫Ωθh(x)μ(x)dx−j=1∑kνjhj
梯度:
∇E(h)=(w1(h)−ν1,…,wk(h)−νk)T\nabla E(h) = (w_1(h) - \nu_1, \dots, w_k(h) - \nu_k)^T∇E(h)=(w1(h)−ν1,…,wk(h)−νk)T
其中 wj(h)=∑x∈Vjμ(x)w_j(h) = \sum_{x \in V_j} \mu(x)wj(h)=∑x∈Vjμ(x) 是第 jjj 个单元的总质量。
Hessian矩阵:
H=∂2E(h)∂hi∂hj={∑l∫filμ(x)dx∥yl−yi∥,i=j,∀l,s.t. fil≠∅−∫fijμ(x)dx∥yj−yi∥,i≠j,fij≠∅0,i≠j,fij=∅H = \frac{\partial^2 E(h)}{\partial h_i \partial h_j} =
\begin{cases}
\sum_l \frac{\int_{f_{il}} \mu(x)dx}{\|y_l - y_i\|}, & i = j, \forall l, \text{s.t. } f_{il} \neq \emptyset \\
-\frac{\int_{f_{ij}} \mu(x)dx}{\|y_j - y_i\|}, & i \neq j, f_{ij} \neq \emptyset \\
0, & i \neq j, f_{ij} = \emptyset
\end{cases}H=∂hi∂hj∂2E(h)=⎩⎨⎧∑l∥yl−yi∥∫filμ(x)dx,−∥yj−yi∥∫fijμ(x)dx,0,i=j,∀l,s.t. fil=∅i=j,fij=∅i=j,fij=∅
其中 fijf_{ij}fij 是相邻单元 ViV_iVi 和 VjV_jVj 的交集。
优化过程:
通过牛顿法迭代求解:
h(t+1)←h(t)−λH−1∇E(h)h^{(t+1)} \leftarrow h^{(t)} - \lambda H^{-1} \nabla E(h)h(t+1)←h(t)−λH−1∇E(h)
4. 几何性质与计算
4.1 凸分割性质
Power Voronoi图具有以下重要性质:
- 每个单元 VjV_jVj 是凸集
- 单元边界是超平面(在2D中是直线,3D中是平面)
- 整个划分是凸分割
4.2 边界计算
相邻单元 ViV_iVi 和 VjV_jVj 之间的边界由以下方程定义:
∥m−yi∥2−ri2=∥m−yj∥2−rj2\|m - y_i\|^2 - r_i^2 = \|m - y_j\|^2 - r_j^2∥m−yi∥2−ri2=∥m−yj∥2−rj2
⇒2mT(yj−yi)=∥yj∥2−∥yi∥2+ri2−rj2\Rightarrow 2m^T(y_j - y_i) = \|y_j\|^2 - \|y_i\|^2 + r_i^2 - r_j^2⇒2mT(yj−yi)=∥yj∥2−∥yi∥2+ri2−rj2
这是一个超平面方程,法向量为 yj−yiy_j - y_iyj−yi。
4.3 体积计算
单元 VjV_jVj 的体积(或质量)为:
wj(h)=∫Vjμ(x)dx=∑x∈Vjμ(x)(离散情况)w_j(h) = \int_{V_j} \mu(x)dx = \sum_{x \in V_j} \mu(x) \quad (\text{离散情况})wj(h)=∫Vjμ(x)dx=x∈Vj∑μ(x)(离散情况)
在优化过程中,需要计算相邻单元交集 fijf_{ij}fij 的质量:
∫fijμ(x)dx\int_{f_{ij}} \mu(x)dx∫fijμ(x)dx
5. 与k-means聚类的联系
5.1 Wasserstein均值问题
当目标分布 YYY 是稀疏的,Wasserstein距离最小化问题退化为:
infY∈P(M)W22(X,Y)=infY∈P(M),π∈P(M×M)∑yj=π(xi)μi∥xi−yj∥2\inf_{Y \in P(M)} W_2^2(X, Y) = \inf_{Y \in P(M), \pi \in P(M \times M)} \sum_{y_j = \pi(x_i)} \mu_i \|x_i - y_j\|^2Y∈P(M)infW22(X,Y)=Y∈P(M),π∈P(M×M)infyj=π(xi)∑μi∥xi−yj∥2
这正是加权k-means聚类问题。
5.2 变分Wasserstein聚类
在变分Wasserstein聚类中:
- Power Voronoi图定义了样本到簇的分配
- 通过调整 hhh,使每个单元的总质量 wj(h)w_j(h)wj(h) 等于目标测度 νj\nu_jνj
- 簇中心 yjy_jyj 通过加权平均更新:
yj=∑x∈Vjμixi∑x∈Vjμiy_j = \frac{\sum_{x \in V_j} \mu_i x_i}{\sum_{x \in V_j} \mu_i}yj=∑x∈Vjμi∑x∈Vjμixi
6. 算法实现
6.1 更新Power Voronoi图的算法
def Variational_OT(X, Y, mu, nu, epsilon):h = np.zeros(k) # 初始化参数while True:# 1. 更新Power Voronoi图V = compute_power_diagram(Y, h)# 2. 计算单元权重w = [sum(mu[x] for x in V[j]) for j in range(k)]# 3. 计算梯度和Hessiangrad = w - nuH = compute_hessian(Y, V, mu)# 4. 牛顿法更新hdelta_h = -np.linalg.solve(H, grad)h = h + delta_h# 5. 检查收敛if np.linalg.norm(grad) < epsilon:breakreturn V, h
6.2 变分Wasserstein聚类算法
def Variational_Wasserstein_Clustering(X, mu, k):# 初始化簇中心Y = initialize_centroids(X, k)nu = compute_target_measure(k) # 通常为均匀分布while not converged:# 1. 更新Power Voronoi图V, h = Variational_OT(X, Y, mu, nu, epsilon)# 2. 更新簇中心for j in range(k):Y[j] = sum(mu[x] * x for x in V[j]) / sum(mu[x] for x in V[j])return Y, V
7. 物理意义与应用
7.1 测度保持映射
Power Voronoi图实现了从源分布到目标分布的测度保持映射:
- 每个Voronoi单元 VjV_jVj 的总质量等于目标测度 νj\nu_jνj
- 没有质量分裂,符合Monge型最优传输
7.2 几何变形解释
考虑二维平面上的变形:
- 当 wj(h)>νjw_j(h) > \nu_jwj(h)>νj:单元 VjV_jVj 质量过大,需要缩小
- 减小 hjh_jhj,使边界向 yjy_jyj 移动
- 当 wj(h)<νjw_j(h) < \nu_jwj(h)<νj:单元 VjV_jVj 质量过小,需要扩大
- 增加 hjh_jhj,使边界远离 yjy_jyj
7.3 应用场景
- 域适应:将源域映射到目标域,同时保持分布特性
- 网格变形:根据曲率重新分布顶点,高曲率区域顶点更密集
- 表示学习:将高维数据嵌入低维空间,保留重要结构
8. 与标准Voronoi图的对比
特性 | 标准Voronoi图 | Power Voronoi图 |
---|---|---|
定义 | ∣m−yj∣≤∣m−yi∣|m - y_j| \leq |m - y_i|∣m−yj∣≤∣m−yi∣ | ∣m−yj∣2−rj2≤∣m−yi∣2−ri2|m - y_j|^2 - r_j^2 \leq |m - y_i|^2 - r_i^2∣m−yj∣2−rj2≤∣m−yi∣2−ri2 |
参数 | 仅中心点 yjy_jyj | 中心点 yjy_jyj 和权重 rj2r_j^2rj2 |
自由度 | 低 | 高(可通过 hjh_jhj 调整) |
测度保持 | 不能 | 可以(通过调整 hjh_jhj) |
应用场景 | 一般空间划分 | 最优传输、测度保持映射 |
9. 总结
Power Voronoi图是标准Voronoi图的推广,其核心数学原理包括:
- 定义基础:通过幂距离 ∥m−yj∥2−rj2\|m - y_j\|^2 - r_j^2∥m−yj∥2−rj2 定义单元边界
- 凸函数联系:与分段线性凸函数 θh(x)=max{⟨x,yj⟩+hj}\theta_h(x) = \max\{\langle x, y_j \rangle + h_j\}θh(x)=max{⟨x,yj⟩+hj} 直接相关
- 最优传输角色:实现Monge型测度保持映射,解决半离散最优传输问题
- 变分原理:通过优化能量函数 E(h)E(h)E(h) 调整划分以满足目标测度
- 几何性质:保持凸分割,边界为超平面
在变分Wasserstein聚类中,Power Voronoi图是连接最优传输与k-means聚类的关键桥梁,它不仅提供了样本到簇的分配,还确保了测度守恒,使算法能够同时优化聚类质量和Wasserstein距离。其数学优雅性和计算效率使其成为处理分布匹配和聚类问题的有力工具。
Kantorovich将Monge型最优传输问题松弛为寻找满足边缘约束的耦合集合中的最小化问题,这一松弛使得问题更具普适性和可解性。以下是详细解释:
1. 问题背景
Monge型最优传输要求找到一个测度保持映射 T:X→YT: X \to YT:X→Y,使得:
infT#μ=ν∫Xc(x,T(x))dμ(x),
\inf_{T_{\#}\mu = \nu} \int_X c(x, T(x)) d\mu(x),
T#μ=νinf∫Xc(x,T(x))dμ(x),
其中 T#μ=νT_{\#}\mu = \nuT#μ=ν 表示 TTT 是测度保持映射(即 μ(T−1(B))=ν(B),∀B⊂Y\mu(T^{-1}(B)) = \nu(B), \forall B \subset Yμ(T−1(B))=ν(B),∀B⊂Y)。然而,这类问题可能存在无解的情况(例如当 XXX 和 YYY 的维度不同或测度不满足正则性条件时)。
2. Kantorovich的松弛
Kantorovich提出将问题转化为寻找联合概率测度 π∈Π(μ,ν)\pi \in \Pi(\mu, \nu)π∈Π(μ,ν),其中:
- 耦合集合 Π(μ,ν)\Pi(\mu, \nu)Π(μ,ν) 定义为所有满足以下边缘约束的联合测度:
- 行约束:π(⋅,Y)=μ\pi(\cdot, Y) = \muπ(⋅,Y)=μ(对任意 A⊂XA \subset XA⊂X,∫Aπ(x,y)dy=μ(A)\int_A \pi(x, y) dy = \mu(A)∫Aπ(x,y)dy=μ(A))
- 列约束:π(X,⋅)=ν\pi(X, \cdot) = \nuπ(X,⋅)=ν(对任意 B⊂YB \subset YB⊂Y,∫Bπ(x,y)dx=ν(B)\int_B \pi(x, y) dx = \nu(B)∫Bπ(x,y)dx=ν(B))
目标函数变为:
infπ∈Π(μ,ν)∫X×Yc(x,y)dπ(x,y),
\inf_{\pi \in \Pi(\mu, \nu)} \int_{X \times Y} c(x, y) d\pi(x, y),
π∈Π(μ,ν)inf∫X×Yc(x,y)dπ(x,y),
其中 c(x,y)c(x, y)c(x,y) 是传输成本函数(通常取为距离的 ppp 次幂,如 c(x,y)=d(x,y)pc(x, y) = d(x, y)^pc(x,y)=d(x,y)p)。
3. 关键思想
- 质量分割:允许一个源点 xxx 的质量被分配到多个目标点 yyy,而非必须全部映射到一个点。
- 联合测度:π(x,y)\pi(x, y)π(x,y) 表示将质量从 xxx 分配到 yyy 的比例。
- 边缘约束:确保总质量守恒:
- 对所有 x∈Xx \in Xx∈X,∫Yπ(x,y)dy=μ(x)\int_Y \pi(x, y) dy = \mu(x)∫Yπ(x,y)dy=μ(x)(源测度的局部质量守恒)。
- 对所有 y∈Yy \in Yy∈Y,∫Xπ(x,y)dx=ν(y)\int_X \pi(x, y) dx = \nu(y)∫Xπ(x,y)dx=ν(y)(目标测度的局部质量守恒)。
4. 数学表达
在离散情况下,若 μ=∑i=1nμiδxi\mu = \sum_{i=1}^n \mu_i \delta_{x_i}μ=∑i=1nμiδxi,ν=∑j=1mνjδyj\nu = \sum_{j=1}^m \nu_j \delta_{y_j}ν=∑j=1mνjδyj,则:
- 耦合 π\piπ 是一个 n×mn \times mn×m 矩阵 {πij}\{\pi_{ij}\}{πij},满足:
∑j=1mπij=μi,∀i, \sum_{j=1}^m \pi_{ij} = \mu_i, \quad \forall i, j=1∑mπij=μi,∀i,
∑i=1nπij=νj,∀j. \sum_{i=1}^n \pi_{ij} = \nu_j, \quad \forall j. i=1∑nπij=νj,∀j. - 目标函数变为:
minπ∑i=1n∑j=1mc(xi,yj)πij, \min_{\pi} \sum_{i=1}^n \sum_{j=1}^m c(x_i, y_j) \pi_{ij}, πmini=1∑nj=1∑mc(xi,yj)πij,
这是一个线性规划问题。
5. 与Monge问题的关系
- 包含关系:Monge问题的解(若存在)对应于Kantorovich问题的一个特殊耦合(即 π(x,y)=δ[y=T(x)]dμ(x)\pi(x, y) = \delta[y = T(x)] d\mu(x)π(x,y)=δ[y=T(x)]dμ(x))。
- 存在性:Kantorovich问题始终存在解(由Kantorovich-Rubinstein定理保证),而Monge问题可能无解。
- 几何意义:Kantorovich解允许质量分裂,更适合处理高维或非光滑问题,但失去了Monge解的显式映射特性。
6. 在Wasserstein距离中的应用
Kantorovich松弛直接导出了Wasserstein距离的定义:
Wp(μ,ν)=(infπ∈Π(μ,ν)∫X×Yd(x,y)pdπ(x,y))1/p.
W_p(\mu, \nu) = \left( \inf_{\pi \in \Pi(\mu, \nu)} \int_{X \times Y} d(x, y)^p d\pi(x, y) \right)^{1/p}.
Wp(μ,ν)=(π∈Π(μ,ν)inf∫X×Yd(x,y)pdπ(x,y))1/p.
这是概率测度空间上的有效距离度量,具有良好的数学性质(如三角不等式)。
7. 实际意义
- 计算优势:Kantorovich问题可通过线性规划求解,而Monge问题通常需要非凸优化。
- 灵活性:适用于更广泛的测度类型(如连续测度、混合型测度)。
- 应用广泛性:在图像处理、机器学习(如GANs)、经济学等领域有广泛应用。
总结
Kantorovich的松弛通过引入耦合测度,将严格的测度保持映射扩展为允许质量分割的联合测度优化,从而解决了Monge问题的局限性。这一松弛不仅保证了问题的可行性,还提供了强大的数学工具(如线性规划)和理论保障(如Wasserstein距离的良定义),成为现代最优传输理论的核心。