数学概念表达

前言

这里记录一下学习过程中遇到的一些数学记号及其对应性质，在下次遇到后可以直接使用不至于陌生。对于成体系的数学记号，则单独归纳为一篇博客。

KL散度

参考自：KL散度

英文名：Kullback-Leibler divergence，是一种用于衡量两个概率分布之间差异的指标，其定义为：

$$D_\mathrm{KL}[P(x)||Q(x)]=\sum\limits_{x\in\mathcal{X}}P(x)\ln\Big[\frac{P(x)}{Q(x)}\Big]$$

KL散度描述的是使用概率分布$Q(x)$来描述分布$P(x)$的编码损失，其性质是非负的、仿射不变的且非对称的。

期望记法

在文献Auto-Encoding Variational Bayes和Denoising Diffusion Probabilistic Models等中经常见到的一种表达式，例如：$\mathbb{E}_q[\bullet]$和$\mathbb{E}_{q_\phi(\mathbf{z|x})}[\bullet]$，似乎是将概率密度当成了随机变量在求均值，乍一看难以理解，实际上，如果设随机变量$x$的事件全集为$\Omega$，则有：

$$\mathbb{E}_x[f(x)]=\sum_{x\in\Omega}xf(x)$$

式中$f$表示某一映射，将随机变量映射为确定的量以用于求解均值，常见的映射如概率密度函数$p(x)$、$p(x|y)$等，也经常在KL散度的表示式中出现对数概率密度$\log p(x)$等。如果在该式中出现了映射函数的乘积，例如$\mathbb{E}_x[f(x)g(x)]$，可进一步简写为：

$$\sum_{x\in\Omega}f(x)g(x)=\sum_{f(x)\in f(\Omega)}f(x)g(x)\triangleq\mathbb{E}_{f(x)}[g(x)]\triangleq\mathbb{E}_{x\sim f(x)}[g(x)]$$

其中后者在GAN网络的论文中曾经出现过，也是一种常见的期望记号。从而：

$$\begin{aligned} D_\mathrm{KL}[q(x)||p(x)]&=\mathbb{E}_{q(x)}[\log q(x)-\log p(x)]\\ &=\mathbb{E}_{x\sim q(x)}[\log q(x)-\log p(x)] \end{aligned}$$

条件熵

条件熵记号，其定义式为：

$$\forall X(e)\in\Omega,H(X|y)=-\sum\limits_{i=1}^np(x_i|y)\ln p(x_i|y)$$

该定义来自于信息论（自然对数或者常用对数都不影响最终结果），当$X$是预测分类的结果，$y$是对应的正确分类时，上式成为分类问题常用的交叉熵。

Fischer信息

参考自维基百科：Fischer信息

Fischer信息反映了参数估计的似然函数与实际的概率密度函数之间的差异，当二者完全相同时，Fischer信息为0，当二者差异越大时，Fischer信息越大。Fischer信息定义为似然函数得分score的方差，即为：

$$I(\theta)=\mathbb{E}\Big[\Big(\frac{\partial}{\partial\theta}\ln L(x;\theta)\Big)^2\Big]\triangleq\Big\langle\Big(\frac{\partial}{\partial\theta}\ln L(x;\theta)\Big)^2\Big\rangle$$

特别地，如果对数似然函数是二阶可微的，上式可以得到进一步化简为：

$$I(\theta)=-\mathbb{E}\Big[\frac{\partial^2}{\partial\theta^2}\ln L(x;\theta)\Big]\triangleq-\Big\langle\frac{\partial^2}{\partial\theta^2}\ln L(x;\theta)\Big\rangle$$

此时，认为似然函数是当前概率分布密度函数的近似，则有：

$$\int\limits_{-\infty}^{+\infty}\frac{\partial^2 L(x;\theta)}{\partial \theta^2}\mathrm{d}x=0$$

因此上式成立，请注意负号的存在。

资格迹

参考自知乎和medium

英文原名：Eligibility Traces，是强化学习中的一种训练方法（另外两种常见的分别是蒙特卡洛方法Monte Carlo methods和时间差分法Time Differetiating），类似于小批量梯度下降和随机梯度下降这种关系，资格迹方法是蒙特卡洛方法和时间差分法的结合。在有值函数和状态-动作值函数的情况下，资格迹方法可以用来更新值函数和状态-动作值函数的参数，如下式：

$$\theta \leftarrow \theta + \alpha \times e(s, a) \times \nabla_\theta Q(s, a)$$

式中，$\alpha$是学习率，$\theta$是值函数或状态-动作值函数的参数，$e(s, a)$是资格迹这里是与迭代步长呈现指数衰减的函数，$Q(s, a)$是值函数或状态-动作对所对应的值函数。

2D射影变换估计

参考自《计算视觉中的多视图几何》Multiple View Geometry in Computer Vision Second Edition。经常遇到一种奇特的代价函数，被称为是“reprojection error”，即重投影误差，今天在课上才发现居然此书中就有，因此学习之。

该问题可以归纳为，空间中两个点列的对应关系由通式：$\mathbf{x_i'=Hx_i}$给出，其中左右分别代表某个点的齐次坐标，变换矩阵是升一维的齐次变换阵。问题即在已知两个点列以及对应关系的情况下，求解这个齐次变换阵。

数学模型

以二维为例，齐次式含有隐式的齐次因子，且不同的点之间这个齐次因子有可能是不一样的（因为同一个齐次坐标可以有无数个表示方式，很有可能该齐次因子是不同的），所以需要在同一个点内部构造相等变换，消去齐次因子，这就需要使用向量外积，有：

$$\mathbf{x_i'}\times\mathbf{Hx_i}=0$$

这里，可以使用外积的矩阵法表示，即有：

$$\mathbf{x_i'x_i^\mathrm{T}H^\mathrm{T}}-\mathbf{Hx_ix_i'^\mathrm{T}}=0$$

如果以$\mathbf{h_1^\mathrm{T},h_2^\mathrm{T},h_3^\mathrm{T}}$作为$\mathbf{H}$矩阵的行向量，以$x_i=(x_i,y_i,z_i)^\mathrm{T}$作为向量的坐标，则上式可以写作：

$$\begin{cases} y_i'\mathbf{h_3^\mathrm{T}x}_i-z_i'\mathbf{h_2^\mathrm{T}x}_i=0\\\\ z_i'\mathbf{h_1^\mathrm{T}x}_i-x_i'\mathbf{h_3^\mathrm{T}x}_i=0\\\\ x_i'\mathbf{h_2^\mathrm{T}x}_i-y_i'\mathbf{h_1^\mathrm{T}x}_i=0 \end{cases}$$

该方程中实际只有两个是独立的，第三个方程可以由前两个方程推出，因此可以只使用前两个方程并改写为更容易记忆的形式：

$$\begin{bmatrix} 0&-z_i\mathbf{x}_i^\mathrm{T}&y_i\mathbf{x}_i^\mathrm{T}\\ z_i\mathbf{x}_i^\mathrm{T}&0&-x_i\mathbf{x}_i^\mathrm{T}\\ \end{bmatrix}\begin{bmatrix} \mathbf{h_1}\\ \mathbf{h_2}\\ \mathbf{h_3} \end{bmatrix}=0\triangleq\mathbf{A}_i\mathcal{V}(\mathbf{H})=0$$

式中，$\mathbf{A}_i$是一个$2\times9$的矩阵，$\mathcal{V}(\mathbf{H})$表示将矩阵$\mathbf{H}$按行展开为一个列向量。求解该方程，一般可以使用奇异值分解求最小特征值的方法，但是为了满足更高精度的要求，需要选用更加合适的向量范数以构造代价函数。另一种思路则是将齐次坐标非齐次化，即某维度归一，但是这样不一定能够得到稳定的解（比如变换出现在无穷远点），因此不推荐使用。

代价函数

1. 代数距离，对于两个向量$\mathbf{x}_1,\mathbf{x}_2$，可以定义其代数距离为：

$$d^2(\mathbf{x}_1,\mathbf{x}_2)=a_1^2+a_2^2\text{ with }\mathbf{a=x_\mathrm{1}\times x}_2$$

即代数距离只计算有效方程贡献的差值，不计算线性相关项。代数距离计算简单且快捷，但是计算结果没有具体的物理意义，因此多用于初次迭代。
2. 几何距离，对于两幅点相互对应的图像，在已知变换矩阵的情况下，可以计算两幅图像中对应点的几何距离，该几何距离就是欧氏距离。

转移误差

在动态视觉或者多视图领域，往往需要将一幅图向另一幅图进行对应，对应的过程引入了转换矩阵，而对应的结果需要进行评估，这就需要引入转移误差的概念，将转移误差优化至最小值，即为算法得到了最佳的转换矩阵估计值$\mathbf{\~H}$。

若设$\mathbf{x}_i$是第一幅图的点，$\mathbf{x}_i'$是第二幅图的点，转移矩阵为$\mathbf{H}$。

1. 单图转移误差
   仅从第一张图向第二张图转换，得到的转换误差可用全部点的几何距离描述，即有：

$$C(\mathbf{H})=\sum\limits_{i}d^2(\mathbf{Hx}_i,\mathbf{x}_i')=\|\mathbf{x}_i'-\mathbf{Hx}_i\|_2^2$$

2. 对称转移误差
   认为转移是两个方向的，即有：

$$C(\mathbf{H})=\sum\limits_{i}d^2(\mathbf{Hx}_i,\mathbf{x}_i')+\sum\limits_{i}d^2(\mathbf{H}^{-1}\mathbf{x}_i',\mathbf{x}_i)=\|\mathbf{x}_i'-\mathbf{Hx}_i\|_2^2+\|\mathbf{x}_i-\mathbf{H}^{-1}\mathbf{x}_i'\|_2^2$$

对称转移误差得到的结果仍然不能保证两幅图中的点是重合的，对于需要估计“校正”结果的情况，该方法不适用，为此，引入了重投影误差。

3. 重投影误差
   核心思想是：在两幅图上找到一组能够转换重合的点（称为完全匹配点，记为点$\mathbf{\~x}'=\mathbf{H\~x}$），计算各自图上点到该附加点的几何距离作为校正误差，称为重投影误差，即有：

$$C(\mathbf{H})=\sum\limits_{i}\Big[d^2(\mathbf{x}_i,\mathbf{\~x}_i)+d^2(\mathbf{x}_i',\mathbf{\~x}_i')\Big]$$

求解过程

求解基于几何距离的重投影误差是非常困难的，通常需要迭代法，另一种可行的方法是使用Sampson误差近似，将代价函数使用Taylor一阶展开式：

$$C(\mathbf{H,X+\delta X})\approx C(\mathbf{H,X})+\frac{C(\mathbf{H,X})}{\partial\mathbf{X}}\delta\mathbf{X}=\mathbf{E+J\delta X}$$

使代价函数为0，即有：

$$\mathbf{E+J\delta X}=0$$

使用拉格朗日乘子法并求解，可以得到Sampson距离，即有：

$$||\mathbf{\delta X}||_s^2=\mathbf{E^\mathrm{T}(JJ^\mathrm{T})^{-1}E}$$

用该距离代替几何距离以计算重投影误差，可以得到Sampson误差近似。

分类结果评价指标

参考自知乎

在机器学习的分类领域经常见到这个准则，即假阳性、真阳性判据，可用于更加全面地计算平均准确率mAP。需要对这个准则作一个总结。

四个基本量

1. TP：真阳性，即预测为正样本且实际为正样本的数量。
2. FP：假阳性，即预测为正样本但实际为负样本的数量。
3. TN：真阴性，即预测为负样本且实际为负样本的数量。
4. FN：假阴性，即预测为负样本但实际为正样本的数量。

三个导出量

1. 精确率：$\displaystyle P=\frac{TP}{TP+FP}$，即预测为正样本的样本中实际为正样本的比例。
2. 召回率：$\displaystyle R=\frac{TP}{TP+FN}$，即实际为正样本的样本中预测为正样本的比例。
3. 准确率：$\displaystyle A=\frac{TP+TN}{TP+TN+FP+FN}$，即预测正确的样本占总样本的比例。又简称为ACC。

两个评价方式

1. F1-score，即$\displaystyle F_1=\frac{(a^2+1)PR}{a^2(P+R)}$，是精确率和召回率的综合评估量，$a$值可经显著性水平查表得取得，可以用于评价分类器的整体性能。
2. ROC曲线，即以假阳性率为横轴，真阳性率为纵轴，绘制的曲线，可以用于评价分类器的整体性能。
   
   曲线越靠近左上（此图中），曲线越远离中间对角性，则分类性能越好。

几个有关图像相似性的评价指标

最原始的指标有逐像素的RME以及峰值信噪比PSNR，但是这两个指标都反映了整体上的概念而对细节上把握程度不够，不能够反映人眼对图像获取到的相似度感受，因此引入了以下几种指标，用于效仿从人眼观察图像相似的感受。

SSIM

参考自知乎
SSIM分为三个部分组成，分别是图像的亮度luminance、对比度contrast以及结构structure。这三个部分分别用$\mu_x,\mu_y$表示，$\sigma_x^2,\sigma_y^2$表示，$\sigma_{xy}$表示，即：

$$\begin{aligned} l(x,y)=\frac{\mu_x\mu_y+C_1}{\mu_x^2+\mu_y^2+C_1}\\ c(x,y)=\frac{\sigma_x\sigma_y+C_2}{\sigma_x^2+\sigma_y^2+C_2}\\ s(x,y)=\frac{\sigma_{xy}+C_3}{\sigma_x\sigma_y+C_3} \end{aligned}$$

因此，SSIM的定义为：

$$\mathrm{SSIM}(x,y)=l(x,y)^{\alpha}c(x,y)^{\beta}s(x,y)^{\gamma}$$

实际中常取$C_3=C_2/2$，因此上式可以进一步化简，此处不再赘述。

MS-SSIM

本指标又简称为MSSIM，是SSIM的改进版，特别之处在于SSIM使用了多次下采样，每次下采样图像尺寸便缩小一半，在下采样之前计算两幅图像的对比度相似度和结构相似度，经过$M-1$次下采样得到缩小为$2^{M-1}$大小的图像，用该图像计算亮度相似度，并将三个相似度按幂权重归一化乘积起来，即有：

FSIM

FSIM指标与以上两种基于“结构相似性度量”Structure Similarity Measurement的指标不同，本指标“首次”引入了特征的概念，使用手工设计的两个正交滤波器对分别对图像作一次处理，然后计算两个处理后图像的相位一致性值PC：

$$PC(\mathbf{x})=\frac{\sum_jE_{\theta(j)}(\mathbf{x})}{\epsilon+\sum_n\sum_jA_{n,\theta(j)}(\mathbf{x})}$$

其中的$A$是在各个像素位置$\mathbf{x}$处的，两个滤波器各自滤波结果的独立和（即平方和再开根号），$E$则是将同一个尺度Scale下全部滤波结果的先作相关和再作独立和（即先求和，然后再平方再开根号）。

图像的梯度幅度Gradient Magnitude计算，作者在原始文献中提供了三种正交算子，分别是Sobel、Prewitt和Scharr算子，使用其中任意一种算子对$x,y$两个方向上计算其偏导数，然后逐像素叠加即可得到在每个像素点处的梯度幅度，称为$G=\sqrt{G_x^2+G^2_y}$。

计算FSIM指标需要通过PC和GM，分为以下三步进行：

$$\begin{aligned} S_{PC}(\mathbf{x})&=\frac{2PC_1(\mathbf{x})\cdot PC_2(\mathbf{x})+T_1}{PC_1^2(\mathbf{x})+ PC_2^2(\mathbf{x})+T_1}\\ S_{G}(\mathbf{x})&=\frac{2G_1(\mathbf{x})\cdot G_2(\mathbf{x})+T_2}{G_1^2(\mathbf{x})+ G_2^2(\mathbf{x})+T_2}\\ S_L(\mathbf{x})&=S_{PC}^\alpha(\mathbf{x})S_G^\beta(\mathbf{x}) \end{aligned}$$

最终计算FSIM的公式即遍历整个位矢$\mathbf{x}$的取值空间（整个图像域）上的平均，即：

$$\mathrm{FSIM}=\frac{\sum_{x\in\Omega}S_L(\mathbf{x})\cdot PC_m(\mathbf{x})}{\sum_{x\in\Omega}PC_m{\mathbf{x}}}$$

实际上对图像相似度评价的指标还有好多种，像HDR-VDP等等，最新的文献The Unreasonable Effectiveness of Deep Features as a Perceptual Metric提出了一种基于特征的新思路，即不用单独引入外部评价方法，将参考图像与图像分别通过相同的网络（往往就是待训练的生成式网络），计算其每一层网络的输出特征图各自对应点的加权欧氏距离和，即：

$$d(\tilde{x},x)=\sum_l\frac{1}{H_lW_l}\sum_{h,w}||w_l\odot (\tilde{y}^l_{hw}-y^l_{hw})||_2^2$$

这里，$l$是代表第$l$层网络输出的结果，$w_l$是各一层网络不同的权重，这个权重可以层与层之间不同也可以相同，甚至可以同一层内也不同，当权重恒等于1时，该距离与计算图像的余弦距离相同。

二分图匹配算法

解决二分图中最优权匹配问题的算法，最著名的莫过于匈牙利算法，在最早提出解决方案的论文^[1]中，二分图匹配问题被称为是人事安排问题，即描述为：现有甲、乙、丙和丁四个人，以及A、B、C、D四件不同的工作，其中根据每个人是否有资格参与对应的工作，列出下表：

	A	B	C	D
甲	是	是	是	否
乙	否	否	是	是
丙	否	否	否	是
丁	否	否	否	是

如果用资格矩阵表示上述表格，可得：

$$\mathbf{Q}=\begin{bmatrix} 1 & 1 & 1 & 0 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 \end{bmatrix}$$

当前条件是，每件工作最多只能让一个人去执行，问题是：如何安排人事，使得被分配出去的工作数量达到最大？

用矩阵的话说就是：在保证只能从某一行或者某一列挑选的前提下，能够从资格矩阵中挑出1的最大数量？

一个简单的解法

按照论文的观点，很容易从没有被分配的工作中加入分配的人员，我们因此首先对两个“多面手”随便分配一个工作，如对甲分配C，对乙分配D，记分配矩阵修改为：

$$\mathbf{A}=\begin{bmatrix} 1 & 1 & 1* & 0 \\ 0 & 0 & 1 & 1* \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 \end{bmatrix}$$

其中"*"号表示的是该项工作被分配。显然该分配方案并不是最优分配，因为只需要将甲从工作C调整至A，将乙从工作D调整至C，就会多出一个空余的工作D，由丙和丁两人任意分配。

能够将工作分配至当前手上尚未有工作的人员的分配方案不是最优方案

上述结论是显然的，其否命题则是：无论如何调整当前人事安排，都不能将任意尚未有工作的人员分配一个工作，则当前分配方案取得最优。

二分图问题

更广义上讲，现给定一个$\mathbb{R}^{n\times n}$的正整数矩阵，要求找到一个行（或者列）序号的排列方式$\mathbf{P}$，使得从每一行（列）中选取一个数$r_{i,j_i}$出来，这些数的和最大，矩阵式的表示更加直观，来自于维基百科：

$$\argmax_{\mathbf{P}}\sum_{i=0}^n r_{i,j_i}=\argmin_{\mathbf{P}}\mathrm{tr}({\mathbf{PC}})$$

其中，$\mathbf{P}$是初等变换阵permutation matrix。参考自@Jinzhong Xu的博客和CSDN博客后，才算是正式弄懂了该算法。用匈牙利算法解决二分图匹配问题的关键是找到矩阵中含0行列的最小划线覆盖，下例将演示完全的匹配过程。

示例

现有成本矩阵如下，其中行代表人员，列代表工作，目标是在保证工作都被完成的前提下，取得最小的成本。

$$\begin{bmatrix} 16 & 12 & 38 & 24 \\ 42 & 23 & 43 & 63 \\ 33 & 54 & 15 & 25 \\ 42 & 25 & 27 & 37 \end{bmatrix}$$

用矩阵表述该问题是，如何找到一种分配方案（对应矩阵的初等变换），使得变换后矩阵的迹最小？

1. 找到各行矩阵的最小值，并用对应行减去（这一步称为归约）：

$$\mathrm{row}\min = \begin{bmatrix} 12\\23\\15\\25 \end{bmatrix}\implies \begin{bmatrix} 4 & 0 & 26 & 12 \\ 19 & 0 & 20 & 40 \\ 18 & 39 & 0 & 10 \\ 17 & 0 & 2 & 12 \end{bmatrix}\implies\mathrm{col}\min = \begin{bmatrix} 4&0&0&10 \end{bmatrix}\implies \begin{bmatrix} 0 & 0 & 26 & 2 \\ 15 & 0 & 20 & 30 \\ 14 & 39 & 0 & 0 \\ 13 & 0 & 2 & 2 \end{bmatrix}$$

2. 用最少行、列直线覆盖当前全部的零，统计当前的划线覆盖数为3，小于当前成本矩阵的最小维度4，因此需要继续归约。

$$\begin{bmatrix} - & + & - & - \\ 15 & | & 20 & 30 \\ - & + & - & - \\ 13 & | & 2 & 2 \end{bmatrix}$$

3. 对于没有被划线覆盖的元素，找到其中最小的值（这里是2）。用未被覆盖的元素减去该最小值，被划线的部分保持不变：

$$\begin{bmatrix} - & + & - & - \\ 13 & | & 18 & 28 \\ - & + & - & - \\ 11 & | & 0 & 0 \end{bmatrix}$$

4. 再划线并查看划线数是否等于最小维度，此时划线结果为4，已经等于最小维度，划线步骤停止：

$$\begin{bmatrix} - & + & - & - \\ 13 & | & 18 & 28 \\ - & + & - & - \\ - & + & - & - \end{bmatrix}$$

5. 从每行零数量最少的矩阵开始迭代，将当前行的0标记为$T_0$，然后将$T_0$所在的列中的其他零标记为$S_0$，如下图示：

$$\begin{bmatrix} 0 & 0 & 26 & 2 \\ 13 & 0 & 18 & 28 \\ 14 & 39 & 0 & 0 \\ 11 & 0 & 0 & 0 \end{bmatrix}\tag{*}\implies\begin{bmatrix} 0 & \blue{S_0} & 26 & 2 \\ 13 & \red{T_0} & 18 & 28 \\ 14 & 39 & 0 & 0 \\ 11 & \blue{S_0} & 0 & 0 \end{bmatrix}$$

6. 将$T_0,S_0$不再视为零，重复以上步骤直到全部的零都被标记，统计其中的$\red{T_0}$零个数为4，已经等于当前最小维度数了，此步骤结束后将停止迭代：

$$\begin{bmatrix} \red{T_0} & \blue{S_0} & 26 & 2 \\ 13 & \red{T_0} & 18 & 28 \\ 14 & 39 & 0 & 0 \\ 11 & \blue{S_0} & 0 & 0 \end{bmatrix}\implies\begin{bmatrix} \red{T_0} & \blue{S_0} & 26 & 2 \\ 13 & \red{T_0} & 18 & 28 \\ 14 & 39 & \red{T_0} & \blue{S_0}\\ 11 & \blue{S_0} & \blue{S_0} & 0 \end{bmatrix}\implies\begin{bmatrix} \red{T_0} & \blue{S_0} & 26 & 2 \\ 13 & \red{T_0} & 18 & 28 \\ 14 & 39 & \red{T_0} & \blue{S_0} \\ 11 & \blue{S_0} & \blue{S_0} & \red{T_0} \end{bmatrix}$$

7. 输出在$\red{T_0}$位置处的原始成本，即为最优安排：

$$\begin{bmatrix} \red{16} & 12 & 38 & 24 \\ 42 & \red{23} & 43 & 63 \\ 33 & 54 & \red{15} & 25 \\ 42 & 25 & 27 & \red{37} \end{bmatrix}$$

当前问题的最优成本为$16+23+15+37=91$，即人员1做工作1，人员2做工作2，人员3做工作3，人员4做工作4。可以验证，该方法确实是最优的。

当划线数等于最小维度数时，就一定能够保证有解存在，否则，需要反复执行步骤3

算法实现

作为非常经典的算法，匈牙利算法早已被C/C++、Python等语言实现过，这里以Wiki百科上提供的源代码为例：

template <class T> vector<T> hungarian(const vector<vector<T>> &C) {
    const int J = (int)size(C), W = (int)size(C[0]);
    assert(J <= W);
    // job[w] = job assigned to w-th worker, or -1 if no job assigned
    // note: a W-th worker was added for convenience
    vector<int> job(W + 1, -1);
    vector<T> h(W);  // Johnson potentials
    vector<T> answers;
    T ans_cur = 0;
    const T inf = numeric_limits<T>::max();
    // assign j_cur-th job using Dijkstra with potentials
    for (int j_cur = 0; j_cur < J; ++j_cur) {
        int w_cur = W;  // unvisited worker with minimum distance
        job[w_cur] = j_cur;
        vector<T> dist(W + 1, inf);  // Johnson-reduced distances
        dist[W] = 0;
        vector<bool> vis(W + 1);     // whether visited yet
        vector<int> prv(W + 1, -1);  // previous worker on shortest path
        while (job[w_cur] != -1) {   // Dijkstra step: pop min worker from heap
            T min_dist = inf;
            vis[w_cur] = true;
            int w_next = -1;  // next unvisited worker with minimum distance
            // consider extending shortest path by w_cur -> job[w_cur] -> w
            for (int w = 0; w < W; ++w) {
                if (!vis[w]) {
                    // sum of reduced edge weights w_cur -> job[w_cur] -> w
                    T edge = C[job[w_cur]][w] - h[w];
                    if (w_cur != W) {
                        edge -= C[job[w_cur]][w_cur] - h[w_cur];
                        assert(edge >= 0);  // consequence of Johnson potentials
                    }
                    if (ckmin(dist[w], dist[w_cur] + edge)) prv[w] = w_cur;
                    if (ckmin(min_dist, dist[w])) w_next = w;
                }
            }
            w_cur = w_next;
        }
        for (int w = 0; w < W; ++w) {  // update potentials
            ckmin(dist[w], dist[w_cur]);
            h[w] += dist[w];
        }
        ans_cur += h[w_cur];
        for (int w; w_cur != W; w_cur = w) job[w_cur] = job[w = prv[w_cur]];
        answers.push_back(ans_cur);
    }
    return answers;
}

对于Python，可以直接调用以下API：

scipy.optimize.linear_sum_assignment

亦有更精细的如对稀疏矩阵特别优化后的算法，仍然在scipy库中可以直接使用。

参考文献

[1] Kuhn H. W. .The hungarian method for the assignment Problem[J/OL].Nav. Res. Logist. Q.,1955,2(1-2):83-97