当前位置：首页 > news >正文

Siknhorn算法介绍

news 2025/7/4 8:37:22

SiknHorn算法是一个快速求解离散优化问题的经典算法，特别适用于计算离散分布之间的**最优传输（Optimal Transport）**距离；

最优传输问题介绍

计算两个概率分布 P 和 Q 之间的传输成本，通常表示为：

$min_{\pi \in U(P,Q)} \sum_{i,j} \pi_{ij} C_{ij}$

$C_{ij}$ 是传输代价矩阵

π 是联合分布（运输计划），满足边缘分布等于 P和 Q；

U(P,Q) 是所有满足边缘分布的有效运输计划的集合;

直接求解此问题的复杂度较高，为 $O\{n^3\}$ 。Sinkhorn算法通过在目标函数中引入正则化项（如Kullback-Leibler散度）将问题转化为更易解的形式.

Sinkhorn正则化的形式

引入熵正则化后，问题变为：

$min_{\pi \in U(P,Q)} \sum_{i,j} \pi_{ij} C_{ij} + \epsilon \sum_{ij} \pi_{ij}log(\pi_{ij})$

其中 ϵ>0 是正则化参数，用来控制正则化项的权重。此时的优化目标是凸的，可以通过迭代方法快速求解。

算法核心思想

Sinkhorn算法利用行列缩放的思想（行列缩放的思想-CSDN博客），将优化问题转化为矩阵的归一化迭代：

初始化：构造一个权重矩阵 K，其元素为： $K_{ij} = e^{-\frac{C_ij}{\epsilon }}$

标量因子： 定义标量因子 u,v 来调整 K的行列和，使其分别等于分布 P和 Q：

迭代更新：

$u \leftarrow \frac{P}{Kv}$ $v \leftarrow \frac{Q}{K^Tu}$

其中 / 表示逐元素相除.

重复迭代直到收敛。

算法步骤

输入：代价矩阵 C，分布 P,Q, 正则化参数 ϵ，收敛阈值 τ

初始化：设置 u=1（全为1的向量），计算 K。

循环：

$u \leftarrow \frac{P}{Kv}$ $v \leftarrow \frac{Q}{K^Tu}$

检查收敛：判断 $\pi = diag(u) Kdiag(v)$ 是否满足精度 τ。

精度 τ是一个用于判断算法是否收敛的阈值。它控制的是最终结果与目标分布之间的误差大小：

误差= $||\pi_r - P|| + ||\pi_c - Q||$

$\pi_r$ 是当前矩阵的行和；

P 是目标行和；

$\pi_c$ 是当前矩阵的列和；

Q是目标列和；

$||\cdot ||$ 表示向量的范数（通常为 ℓ1 或 ℓ2 范数）。

输出：最终的传输计划 π 和传输成本。

import numpy as npdef sinkhorn_algorithm(C, r, c, epsilon=1e-3, max_iter=1000, tol=1e-6):"""Sinkhorn算法计算最优传输问题的近似解。参数：C (numpy.ndarray): 传输代价矩阵 (n, m)。r (numpy.ndarray): 源分布 (n,)。c (numpy.ndarray): 目标分布 (m,)。epsilon (float): 正则化参数，默认为 1e-3。max_iter (int): 最大迭代次数。tol (float): 收敛阈值，默认为 1e-6。返回：pi (numpy.ndarray): 近似的最优传输计划矩阵。transport_cost (float): 最优传输距离。"""# 确保分布为 numpy 数组并且是列向量形式r = np.array(r, dtype=np.float64)c = np.array(c, dtype=np.float64)# 初始化 K 矩阵，K[i, j] = exp(-C[i, j] / epsilon)K = np.exp(-C / epsilon)# 初始化缩放因子 u 和 vu = np.ones_like(r)v = np.ones_like(c)# 迭代更新 u 和 vfor iteration in range(max_iter):u_prev = u.copy()  # 保存上一轮的 u 以判断收敛u = r / (K @ v)  # 更新行缩放因子v = c / (K.T @ u)  # 更新列缩放因子# 判断是否收敛if np.allclose(u, u_prev, atol=tol):break# 计算最终的传输计划矩阵 pipi = np.diag(u) @ K @ np.diag(v)# 计算最优传输成本transport_cost = np.sum(pi * C)return pi, transport_cost# 示例用法
if __name__ == "__main__":# 定义代价矩阵 (3x3)C = np.array([[4, 8, 6],[3, 7, 5],[2, 4, 6]])# 定义源分布和目标分布r = np.array([0.5, 0.3, 0.2])  # 源分布c = np.array([0.4, 0.4, 0.2])  # 目标分布# 调用 Sinkhorn 算法pi, cost = sinkhorn_algorithm(C, r, c, epsilon=1e-2, max_iter=500, tol=1e-6)# 输出结果print("传输计划矩阵 pi:")print(pi)print(f"最优传输距离: {cost}")

查看全文

http://www.lryc.cn/news/496412.html