Munkres 分配算法

匈牙利方法（或 Kuhn 算法）是由4个基本步骤组成的迭代过程。该方法使用“最小行集”覆盖“操纵”成本矩阵的零点，当所需的“最小行集”等于给定成本矩阵的维数时，过程终止。

Munkres 算法是一种最优分配算法，可以看作是匈牙利算法的一个变种。它的不同之处在于处理过程中：

• 包含所有零的“最小行集”
• 独立零点的“最大集”

该算法的一个重要特征是它的发明人在介绍该算法的论文中还导出了一个方程来计算“完全解决任何 $N\times N$ 指派问题所需的最大运算数”。该公式由以下公式给出：
$$N_{\max }=(n/6)(11n^2+12n+31)$$
这比完全枚举方法所需的操作（$!n$）小的多。

Munkres 算法

算法流程如下图所示：

以下6步算法是原始 Munkres 的分配算法（有时称为匈牙利算法）的改进形式。该算法描述了如何通过对 $0$ 加星和加撇以及覆盖和揭开行和列来手动操作二维矩阵。这是因为，在出版时（1957年），很少有人可以使用计算机，并且算法是手工操作的。

• 步骤0：创建一个称为成本矩阵的 $n\times m$ 矩阵，其中每个元素代表将 $n$ 个工人之一分配给 $m$ 个工作之一的成本。旋转矩阵，使列至少与行一样多，并令 $k=\min (n,m)$。
• 步骤1：对于矩阵的每一行，找到最小的元素，并从该行的每个元素中减去它。转到步骤2。
• 步骤2：在结果矩阵中找到零（$z$）。如果其所在行和列中没有加星号的零，则将 $z$ 标星。对矩阵中的每个元素重复此操作。转到步骤3。
• 步骤3：覆盖包含加星零的列。
  • 如果总共覆盖了 $k$ 列，则加星号的零表示完整的唯一分配集。在这种情况下，请转到“完成”；
  • 否则，请转到步骤4。
• 步骤4：
  • 找到一个未覆盖的零并将其加撇。
  • 如果在包含该加撇零的行中没有加星号的零，请转到步骤5。
  • 否则，请覆盖该行并揭开包含加星号的零的列。
  • 以这种方式继续，直到没有剩余的零为止。
  • 保存最小的未覆盖值，然后转到步骤6。
• 步骤5：按以下步骤构造一系列交替的加撇和加星号的零。令 $z_0$ 表示在步骤4中发现的未覆盖的加撇零。令 $z_1$ 表示 $z_0$ 所在列中的加星的零（如果有）。令 $z_2$ 表示 $z_1$ 所在行中的加撇零（总会有一个）。继续该序列直到加撇零所在列中没有加星的零时终止。取消序列中零的星号，并对加撇零加星，擦除所有撇号并揭开矩阵中的每一行。返回步骤3。
• 步骤6：将在步骤4中找到的值加到每个覆盖行的每个元素，并从每个未覆盖列的每个元素中减去它。返回第4步，而不更改任何星标、撇号或覆盖线。
• 完成：分配对由成本矩阵中加星号的零的位置指示。如果 $C(i,j)$ 是加星零，则将与行 $i$ 关联的元素分配给与列 $j$ 关联的元素。

在第4步中，可能的情况是存在一个未覆盖且加撇的零，如果在其行中没有加星号的零，则程序进入步骤5；如果根本没有未覆盖的零，程序转到步骤6。 步骤5为增广路径算法，步骤6修改成本矩阵。

Munkres’ Assignment Algorithm Modified for Rectangular Matrices 中提供的 C# 程序质量上乘。但考虑到 mcximing/hungarian-algorithm-cpp 应用较为广泛还是进行一下分析。

HungarianAlgorithm

问题来了，既然定义了类为何不将所用变量定义为成员？

public:
	HungarianAlgorithm();
	~HungarianAlgorithm();
	double Solve(vector <vector<double> >& DistMatrix, vector<int>& Assignment);

private:
	void assignmentoptimal(int *assignment, double *cost, double *distMatrix, int nOfRows, int nOfColumns);
	void buildassignmentvector(int *assignment, bool *starMatrix, int nOfRows, int nOfColumns);
	void computeassignmentcost(int *assignment, double *cost, double *distMatrix, int nOfRows);
	void step2a(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim);
	void step2b(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim);
	void step3(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim);
	void step4(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim, int row, int col);
	void step5(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim);

HungarianAlgorithm::Solve

由DistMatrix构造distMatrixIn，输入矩阵为 MATLAB 的列优先存储。

	unsigned int nRows = DistMatrix.size();
	unsigned int nCols = DistMatrix[0].size();

	double *distMatrixIn = new double[nRows * nCols];
	int *assignment = new int[nRows];
	double cost = 0.0;

	// Fill in the distMatrixIn. Mind the index is "i + nRows * j".
	// Here the cost matrix of size MxN is defined as a double precision array of N*M elements. 
	// In the solving functions matrices are seen to be saved MATLAB-internally in row-order.
	// (i.e. the matrix [1 2; 3 4] will be stored as a vector [1 3 2 4], NOT [1 2 3 4]).
	for (unsigned int i = 0; i < nRows; i++)
		for (unsigned int j = 0; j < nCols; j++)
			distMatrixIn[i + nRows * j] = DistMatrix[i][j];

HungarianAlgorithm::assignmentoptimal 中为主要实现。
求出的结果给Assignment，为何不直接传入？

	// call solving function
	assignmentoptimal(assignment, &cost, distMatrixIn, nRows, nCols);

	Assignment.clear();
	for (unsigned int r = 0; r < nRows; r++)
		Assignment.push_back(assignment[r]);

	delete[] distMatrixIn;
	delete[] assignment;
	return cost;

HungarianAlgorithm::assignmentoptimal

assignment为每个轨迹的匹配结果。

	double *distMatrix, *distMatrixTemp, *distMatrixEnd, *columnEnd, value, minValue;
	bool *coveredColumns, *coveredRows, *starMatrix, *newStarMatrix, *primeMatrix;
	int nOfElements, minDim, row, col;

	/* initialization */
	*cost = 0;
	for (row = 0; row<nOfRows; row++)
		assignment[row] = -1;

distMatrixEnd为矩阵末尾。
由distMatrixIn构造distMatrix，防止数据被覆盖。
丧心病狂，newStarMatrix一路传入 HungarianAlgorithm::step4 才使用。

	/* generate working copy of distance Matrix */
	/* check if all matrix elements are positive */
	nOfElements = nOfRows * nOfColumns;
	distMatrix = (double *)malloc(nOfElements * sizeof(double));
	distMatrixEnd = distMatrix + nOfElements;

	for (row = 0; row<nOfElements; row++)
	{
		value = distMatrixIn[row];
		if (value < 0)
			cerr << "All matrix elements have to be non-negative." << endl;
		distMatrix[row] = value;
	}


	/* memory allocation */
	coveredColumns = (bool *)calloc(nOfColumns, sizeof(bool));
	coveredRows = (bool *)calloc(nOfRows, sizeof(bool));
	starMatrix = (bool *)calloc(nOfElements, sizeof(bool));
	primeMatrix = (bool *)calloc(nOfElements, sizeof(bool));
	newStarMatrix = (bool *)calloc(nOfElements, sizeof(bool)); /* used in step4 */

从数量少的方向开始，寻找距离最小值。按列查找会面临访问不连续的问题。为什么不用二维数组访问？

直接修改distMatrix，每行（列）减去其中的最小值。如不存在冲突（一对多），将coveredColumns相应位置置为true，相当于匹配 $M$。

distMatrixTemp名字略显怪异，仅是一个访问元素的指针。
starMatrix记录当前的匹配。
能使用 std::vector 管理内存吗？
不对称，行更多时设置coveredRows[row]为false。

	/* preliminary steps */
	if (nOfRows <= nOfColumns)
	{
		minDim = nOfRows;

		for (row = 0; row<nOfRows; row++)
		{
			/* find the smallest element in the row */
			distMatrixTemp = distMatrix + row;
			minValue = *distMatrixTemp;
			distMatrixTemp += nOfRows;
			while (distMatrixTemp < distMatrixEnd)
			{
				value = *distMatrixTemp;
				if (value < minValue)
					minValue = value;
				distMatrixTemp += nOfRows;
			}

			/* subtract the smallest element from each element of the row */
			distMatrixTemp = distMatrix + row;
			while (distMatrixTemp < distMatrixEnd)
			{
				*distMatrixTemp -= minValue;
				distMatrixTemp += nOfRows;
			}
		}

		/* Steps 1 and 2a */
		for (row = 0; row<nOfRows; row++)
			for (col = 0; col<nOfColumns; col++)
				if (fabs(distMatrix[row + nOfRows*col]) < DBL_EPSILON)
					if (!coveredColumns[col])
					{
						starMatrix[row + nOfRows*col] = true;
						coveredColumns[col] = true;
						break;
					}
	}
	else /* if(nOfRows > nOfColumns) */
	{
		minDim = nOfColumns;

		for (col = 0; col<nOfColumns; col++)
		{
			/* find the smallest element in the column */
			distMatrixTemp = distMatrix + nOfRows*col;
			columnEnd = distMatrixTemp + nOfRows;

			minValue = *distMatrixTemp++;
			while (distMatrixTemp < columnEnd)
			{
				value = *distMatrixTemp++;
				if (value < minValue)
					minValue = value;
			}

			/* subtract the smallest element from each element of the column */
			distMatrixTemp = distMatrix + nOfRows*col;
			while (distMatrixTemp < columnEnd)
				*distMatrixTemp++ -= minValue;
		}

		/* Steps 1 and 2a */
		for (col = 0; col<nOfColumns; col++)
			for (row = 0; row<nOfRows; row++)
				if (fabs(distMatrix[row + nOfRows*col]) < DBL_EPSILON)
					if (!coveredRows[row])
					{
						starMatrix[row + nOfRows*col] = true;
						coveredColumns[col] = true;
						coveredRows[row] = true;
						break;
					}
		for (row = 0; row<nOfRows; row++)
			coveredRows[row] = false;

	}

HungarianAlgorithm::step2b 统计覆盖列的数量。

HungarianAlgorithm::computeassignmentcost 计算分配方案的成本。

	/* move to step 2b */
	step2b(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);

	/* compute cost and remove invalid assignments */
	computeassignmentcost(assignment, cost, distMatrixIn, nOfRows);

	/* free allocated memory */
	free(distMatrix);
	free(coveredColumns);
	free(coveredRows);
	free(starMatrix);
	free(primeMatrix);
	free(newStarMatrix);

	return;

HungarianAlgorithm::step2b

统计覆盖列的数量。判断是否完成。

HungarianAlgorithm::buildassignmentvector 由星号获得结果。
HungarianAlgorithm::step3

	int col, nOfCoveredColumns;

	/* count covered columns */
	nOfCoveredColumns = 0;
	for (col = 0; col<nOfColumns; col++)
		if (coveredColumns[col])
			nOfCoveredColumns++;

	if (nOfCoveredColumns == minDim)
	{
		/* algorithm finished */
		buildassignmentvector(assignment, starMatrix, nOfRows, nOfColumns);
	}
	else
	{
		/* move to step 3 */
		step3(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
	}

HungarianAlgorithm::buildassignmentvector

如果starMatrix相应位置做了标记，则将匹配结果记录到assignment数组。

	int row, col;

	for (row = 0; row<nOfRows; row++)
		for (col = 0; col<nOfColumns; col++)
			if (starMatrix[row + nOfRows*col])
			{
#ifdef ONE_INDEXING
				assignment[row] = col + 1; /* MATLAB-Indexing */
#else
				assignment[row] = col;
#endif
				break;
			}

HungarianAlgorithm::step3

缩进结构太深了。
找到一个未覆盖的零并将其加撇。如果在包含该加撇零的行中没有加星号的零，调用 HungarianAlgorithm::step4。
否则，覆盖该行并揭开包含加星号的零的列。以这种方式继续，直到没有剩余的零为止。
保存最小的未覆盖值，然后调用 HungarianAlgorithm::step5。

	bool zerosFound;
	int row, col, starCol;

	zerosFound = true;
	while (zerosFound)
	{
		zerosFound = false;
		for (col = 0; col<nOfColumns; col++)
			if (!coveredColumns[col])
				for (row = 0; row<nOfRows; row++)
					if ((!coveredRows[row]) && (fabs(distMatrix[row + nOfRows*col]) < DBL_EPSILON))
					{
						/* prime zero */
						primeMatrix[row + nOfRows*col] = true;

						/* find starred zero in current row */
						for (starCol = 0; starCol<nOfColumns; starCol++)
							if (starMatrix[row + nOfRows*starCol])
								break;

						if (starCol == nOfColumns) /* no starred zero found */
						{
							/* move to step 4 */
							step4(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim, row, col);
							return;
						}
						else
						{
							coveredRows[row] = true;
							coveredColumns[starCol] = false;
							zerosFound = true;
							break;
						}
					}
	}

	/* move to step 5 */
	step5(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);

HungarianAlgorithm::step4

指定了行列索引。
newStarMatrix先整体拷贝starMatrix，略显低效。

	int n, starRow, starCol, primeRow, primeCol;
	int nOfElements = nOfRows*nOfColumns;

	/* generate temporary copy of starMatrix */
	for (n = 0; n<nOfElements; n++)
		newStarMatrix[n] = starMatrix[n];

	/* star current zero */
	newStarMatrix[row + nOfRows*col] = true;

	/* find starred zero in current column */
	starCol = col;
	for (starRow = 0; starRow<nOfRows; starRow++)
		if (starMatrix[starRow + nOfRows*starCol])
			break;

	while (starRow<nOfRows)
	{
		/* unstar the starred zero */
		newStarMatrix[starRow + nOfRows*starCol] = false;

		/* find primed zero in current row */
		primeRow = starRow;
		for (primeCol = 0; primeCol<nOfColumns; primeCol++)
			if (primeMatrix[primeRow + nOfRows*primeCol])
				break;

		/* star the primed zero */
		newStarMatrix[primeRow + nOfRows*primeCol] = true;

		/* find starred zero in current column */
		starCol = primeCol;
		for (starRow = 0; starRow<nOfRows; starRow++)
			if (starMatrix[starRow + nOfRows*starCol])
				break;
	}

	/* use temporary copy as new starMatrix */
	/* delete all primes, uncover all rows */
	for (n = 0; n<nOfElements; n++)
	{
		primeMatrix[n] = false;
		starMatrix[n] = newStarMatrix[n];
	}
	for (n = 0; n<nOfRows; n++)
		coveredRows[n] = false;

	/* move to step 2a */
	step2a(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);

HungarianAlgorithm::step5

计算： δ l = min ⁡ u ∈ S , v ∉ T { l ( u ) + l ( v ) − w e i g h t ( ( u , v ) ) } \delta_l = \min_{u\in S,v\notin T}\{l(u) + l(v) − \mathrm{weight}((u, v))\} δl​=minu∈S,v∈/T​{l(u)+l(v)−weight((u,v))}

找到最小的未发现元素h。

	double h, value;
	int row, col;

	/* find smallest uncovered element h */
	h = DBL_MAX;
	for (row = 0; row<nOfRows; row++)
		if (!coveredRows[row])
			for (col = 0; col<nOfColumns; col++)
				if (!coveredColumns[col])
				{
					value = distMatrix[row + nOfRows*col];
					if (value < h)
						h = value;
				}

改进标签 $l\to l^{\prime }$：
$$\begin{array}{r}{l}^{\prime }(r)=\{\begin{array}{ll}l(r)-{\delta }_l& \text{if }r\in S\\ l(r)+{\delta }_l& \text{if }r\in T\\ l(r)& \text{otherwise }\end{array}\end{array}$$
将h加到distMatrix每个覆盖的行，未覆盖的列则减去它。

	/* add h to each covered row */
	for (row = 0; row<nOfRows; row++)
		if (coveredRows[row])
			for (col = 0; col<nOfColumns; col++)
				distMatrix[row + nOfRows*col] += h;

	/* subtract h from each uncovered column */
	for (col = 0; col<nOfColumns; col++)
		if (!coveredColumns[col])
			for (row = 0; row<nOfRows; row++)
				distMatrix[row + nOfRows*col] -= h;

返回 HungarianAlgorithm::step3

	/* move to step 3 */
	step3(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);

HungarianAlgorithm::step2a

遍历每一列，如果元素中有加星的零则覆盖该列。

	bool *starMatrixTemp, *columnEnd;
	int col;

	/* cover every column containing a starred zero */
	for (col = 0; col<nOfColumns; col++)
	{
		starMatrixTemp = starMatrix + nOfRows*col;
		columnEnd = starMatrixTemp + nOfRows;
		while (starMatrixTemp < columnEnd){
			if (*starMatrixTemp++)
			{
				coveredColumns[col] = true;
				break;
			}
		}
	}

HungarianAlgorithm::step2b

	/* move to step 3 */
	step2b(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);

参考资料：

• 机器学习之线性回归法
• Functions for the rectangular assignment problem
• cpp/frams/model/similarity/hungarian/hungarian.cpp@869
• HungarianAlg.h
• Smorodov/Multitarget-tracker
• Fast algorithms for finding optimal matchings on weighted bipartite graphs
• assignment-problem-and-hungarian-algorithm
• santhosh-kumar/MultiviewAppearanceTracker
• Fast Track developer guide
• Assignment of Vehicles to Routes by Latex Foam Rubber Products Limited-Kumasi
• Segmentation of hands & Tracking bats and cells
• pmtk3/toolbox/Algorithms/optimization/assignmentoptimal.m
• assignmentoptimal.cpp
• Fast Track 4.2.0
• The Hungarian algorithm
• Munkres Assignment Algorithm
• jcheong0428/hungarian3.py
• Munkres’ Assignment Algorithm Modified for Rectangular Matrices
• The Assignment Problem: Searching for An Optimal and Efficient Solution
• Hungarian method cost assignment. - Dyalog APL
• Concurrent Implement at ion Of Munkres Algorithm