Siknhorn算法介绍

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

最优传输问题介绍

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

        是传输代价矩阵 

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

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

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

Sinkhorn正则化的形式

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

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

算法核心思想

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

        初始化：构造一个权重矩阵 K，其元素为： 

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

        迭代更新：

 

        其中 / 表示逐元素相除.

        重复迭代直到收敛。

算法步骤

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

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

        循环：

         

        检查收敛：判断  是否满足精度 τ。 

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

                误差=

                是当前矩阵的行和；

                P 是目标行和；

                是当前矩阵的列和；

                Q是目标列和；

                表示向量的范数（通常为 ℓ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}")