二分图匹配---Munkres算法（匈牙利算法）
在任务指派问题（如n项工作由n个人承担，每个人完成不同工作所花时间不同，那如何分配使得花费的时间最少）以及一些多目标检测任务中的数据关联部分（如一个目标有多个特征点，有多个目标时检测到的特征点属于哪一个目标的问题）常常会看到Munkres算法，这里从原理及实现上简单介绍一下Munkres算法。
一、Munkres算法介绍
Munkres算法也称为匈牙利算法，匈牙利算法起源于1955年 Kuhn 提出的指派问题解法，kuhn引用匈牙利数学家康尼格关于矩阵中0元素的定理:系数矩阵中独立0元素的最多个数等于能覆盖所有0元素的最少直线数，该解法称为匈牙利方法（Hungarian Method），1956年MERRILL M.FLOOD 给出了匈牙利方法的算法实现步骤，1957年James Munkres针对该方法做了改进，后来大家习惯叫匈牙利算法或Kuhn-Munkres(KM)算法或者Munkres算法，后续一些改进都沿用了这个名称，其主要用来解决分配问题。
如下图a,b,c,d四个工人完成p,q,r,s四项工作，矩阵C(i, j)表示每个工人完成对应一项工作所需的时间，经过求解后工人a对应工作q，b对应s，c对应r，d对应p，最少是23个单位的时间。
二、Munkres算法条件
- 系数矩阵C为n阶方阵
- 系数矩阵每个元素$C_{ij} \geq 0 $
- 目标是求解最小值min
三、Munkres算法
接下来以实际例子来展开Munkres算法基本步骤。公司接到EJGR四项任务，如何给甲乙丙丁四名员工分配可以使得总花费时间最少？
构建成数学模型如下：
3.1 矩阵变换
- 步骤1：系数矩阵的每行减去改行的最小值
- 步骤2：系数矩阵的每列减去该列的最小值
这个利用了指派问题的一个性质：系数矩阵C的一行（列）各元素减去该行（列）最小元素后得到的新系数矩阵B，B的最优解和C的最优解相同
那上述例子系数矩阵就变为：
3.2 试指派
经过第一步变换后，系数矩阵中出现了0元素，如果能找到n个独立的0元素，那系数矩阵将这些独立零元素所在位置记为1，其余设为0，即可得到问题的最优解。试指派可以分为以下几步：
- 步骤1：从只有一个0元素的行开始，给这个0元素加圈\(\bigcirc\) （表示该行所代表的人只有一种任务可指派），然后划去该零元素所在列的其他0元素(表明该任务已指派)，标记为\(\Phi\)
- 步骤2：对只有一个0元素的列加圈 \(\bigcirc\) ，然后划去该行其他0元素，记为\(\Phi\)
- 步骤3：反复进行步骤1和步骤2，直到所有的0元素都被加圈或者划去为止
- 步骤4：如果经过步骤3之后还有未被加圈的0元素，且同一行至少有两个0元素（表明这个人可以指派多项任务），此时可以比较这些0元素对应列0元素较少的0元素加圈，然后划掉此列的0元素（表明选择性多的礼让选择性少的）进行试探，反复进行直到所有0元素都被加圈或者划掉
- 步骤5：如果加圈的0元素个数m刚好等于系数矩阵的阶数n，则问题得到求解。如果m < n，则继续下一大步骤。
上面的例子经过试指派前面三个步骤得到：
此时加圈的0元素个数刚好也为4，所以
也就是说甲分配任务R，乙分配任务J，丙分配任务E，丁分配任务G，共花费28小时。
现在公司接到新的任务ABCDE,且新来了一个员工，其系数矩阵如下：
经过步骤3.1和3.2后得到如下情况：
这时加圈的0元素个数为4 < 5(矩阵的阶数)，此时就是试分配的步骤5，需要按3.3 作最少的直线覆盖所有零元素
3.3 作最少的直线覆盖所有零元素
这一步利用了前面提到的0元素的定理，其大体可以分为以下几步：
- 步骤1： 对没有加圈\(\bigcirc\)零元素的行打上\(\surd\)
- 步骤2： 对打\(\surd\)的行中含有标记为\(\Phi\)的列打上\(\surd\)
- 步骤3： 再对打\(\surd\)的列中含有加圈\(\bigcirc\)零元素的行打上\(\surd\)
- 步骤4： 重复步骤2和步骤3，直到不能新增\(\surd\)的行和列为止
- 步骤5： 对没有打\(\surd\)的行画一横线，对打\(\surd\)的列画上一竖线，这就得到了覆盖所有零元素的最少直线数L
如果L<n，则说明需要对系数矩阵进一步进行变换，转到3.4 对系数矩阵变换以增加0元素
如果L==n 而m < n，则回到3.2的步骤4，另外选择一个零元素继续试探。
新的示例中系数矩阵经过3.3后变为：
此时直线数L=4 < 5(矩阵阶数)，所以需要执行3.4
3.4 系数矩阵变换增加0元素
这一步的目的是增加0元素个数，分为以下几个步骤：
- 步骤1： 在没有被直线覆盖的0元素中找出最小元素值min
- 步骤2： 对打上\(\surd\)的行的各元素都减去min
- 步骤3： 对打上\(\surd\)的列的各元素都加上min，这样可以保证原来的0元素不变。
- 步骤4： 对变换后的系数矩阵执行3.2和3.3，直到找到最优解
在上述例子中，未被直线覆盖的元素中最小值为2，所以经过变换后系数矩阵变为：
再经过3.2得到
此时可以找到5个加圈的0元素，所以将其变为解矩阵：
四、Munkres实现
python实现可以参考：https://software.clapper.org/munkres/
C++实现可以参考： https://brc2.com/the-algorithm-workshop/
注意：在C++实现与前面提到的基本步骤有些差别（这样实现更快），参考链接里面步骤5涉及到了增广路径，这里简单讲一下其含义：
增广路径
其作用是从工人出发找到一条到任务零元素路径（表述可能不太准确），举个例子应该就明白了：
如下图所示，中间*号就是前面提到的加圈\(\bigcirc\)，加一撇是前面提到的\(\Phi\)，此时q任务分配给了a，p任务分配给了b，但此时c只能分配p，但p已经分配给了b。这时就出现了"礼让"的情况，就是选择多的让给选择少的（此时系数都是0，所以让了之后不会改变最终解），因为b还可以执行q任务，b将任务p让给了c，那就变成了p>c, 同理a把任务让给b，最后就变为了r>a, q>b, p>c。 这中间就出现了一条增广路径c>p>b>q>r。
C++ Code
参照后自己用C++实现了下，这里给出具体的代码放在本文最后附录。
五、衍生问题
前面提到了采用Munkres算法需要满足3个条件，但实际情况下可能并不满足，这里参考：https://blog.csdn.net/Wonz5130/article/details/80678410，给出了几种衍生问题。
5.1 如果不是求最小值，是求最大值Max如何处理？
找出系数矩阵中的最大值max，然后令系数矩阵变为max-系数矩阵各元素值，得到新系数矩阵，再按照正常匈牙利法即可求到。
5.2 系数矩阵不是方阵如何处理？
- 工作 < 人数： 增加虚拟工作，虚拟工作对应的系数全设为0
- 人数 < 工作： 增加虚拟工人，虚拟工人对应的系数全设为0
5.3 一个人可以做几件事
如果一个人可以做k件事，那就把这个人的系数复制k份，相当于增加了虚拟人数，如果不是方阵，再构建虚拟工作
5.4 某个人不能做某项工作
将问题转换为求最小值后，如果第i个人不能做第j项工作，则将系数\(C_{ij}\)设置为一个比系数矩阵中最大值大一些的数。
六、其他扩展
二分图最小分配问题还可以采用分枝定界法
- 分支: 将一个问题不断细分为若干子问题, 之后逐个讨论子问题;
- 定界: 分支很多的情况下, 需要讨论的情况也随之增多, 这里就需要定界, 决定在什么时候不在进行分支;
定界的作用是剪掉没有讨论意义的分支, 只讨论有意义的分支
具体内容可参考：
https://www.geeksforgeeks.org/job-assignment-problem-using-branch-and-bound/
http://www.cs.umsl.edu/~sanjiv/classes/cs5130/lectures/bb.pdf.
七、参考链接
https://brc2.com/the-algorithm-workshop/
https://software.clapper.org/munkres/
https://blog.csdn.net/Wonz5130/article/details/80678410
https://zhuanlan.zhihu.com/p/62981901
https://www.geeksforgeeks.org/hungarian-algorithm-assignment-problem-set-1-introduction/
https://www.geeksforgeeks.org/hungarian-algorithm-for-assignment-problem-set-2-implementation/.
钱颂迪 《运筹学》第四版第6章
八、附录
C++代码
class Munkres
{
public:
Munkres(int nRow, int nCol);
~Munkres();
bool Build(const std::vector<std::vector<int> >& costs);
float Compute(const std::vector<std::vector<int> >& costs, std::vector<int>& pairs);
private:
void StepOne(int& step);
void StepTwo(int& step);
void StepThree(int& step);
void StepFour(int& step);
void StepFive(int& step);
void StepSix(int& step);
void StepSeven(std::vector<int>& pairs);
void FindAZero(int& row, int& col);
bool StarInRow(int row);
void FindStarInRow(int row, int& col);
void FindStarInCol(int col, int& row);
void FindPrimeInRow(int row, int& col);
void AugmentPath();
void ClearCovers();
void ErasePrimes();
void FindSmallest(float& minVal);
void ShowCostMatrix();
void ShowMaskMatrix();
private:
int nRow_;
int nCol_;
std::vector<std::vector<int>> costs_; // nRow_ * nCol_
std::vector<std::vector<int>> augPaths_; // nRow_ * 2
std::vector<char> rowCover_; // 1 * nRow_
std::vector<char> colCover_; // 1 * nCol_
std::vector<std::vector<char>> marks_; // nRow_ * nCol_
int pathRow0;
int pathCol0;
int pathCount;
};
/*
* Ref: https://brc2.com/the-algorithm-workshop/
*
* The following 6-step algorithm is a modified form of the original Munkres’ Assignment Algorithm (sometimes referred to as the Hungarian Algorithm). This algorithm describes to the manual manipulation of a two-dimensional matrix by starring and priming zeros and by covering and uncovering rows and columns. This is because, at the time of publication (1957), few people had access to a computer and the algorithm was exercised by hand.
Step 0: Create an nxm matrix called the cost matrix in which each element represents the cost of assigning one of n workers to one of m jobs. Rotate the matrix so that there are at least as many columns as rows and let k=min(n,m).
Step 1: For each row of the matrix, find the smallest element and subtract it from every element in its row. Go to Step 2.
Step 2: Find a zero (Z) in the resulting matrix. If there is no starred zero in its row or column, star Z. Repeat for each element in the matrix. Go to Step 3.
Step 3: Cover each column containing a starred zero. If K columns are covered, the starred zeros describe a complete set of unique assignments. In this case, Go to DONE, otherwise, Go to Step 4.
Step 4: Find a noncovered zero and prime it. If there is no starred zero in the row containing this primed zero, Go to Step 5. Otherwise, cover this row and uncover the column containing the starred zero. Continue in this manner until there are no uncovered zeros left. Save the smallest uncovered value and Go to Step 6.
Step 5: Construct a series of alternating primed and starred zeros as follows. Let Z0 represent the uncovered primed zero found in Step 4. Let Z1 denote the starred zero in the column of Z0 (if any). Let Z2 denote the primed zero in the row of Z1 (there will always be one). Continue until the series terminates at a primed zero that has no starred zero in its column. Unstar each starred zero of the series, star each primed zero of the series, erase all primes and uncover every line in the matrix. Return to Step 3.
Step 6: Add the value found in Step 4 to every element of each covered row, and subtract it from every element of each uncovered column. Return to Step 4 without altering any stars, primes, or covered lines.
DONE: Assignment pairs are indicated by the positions of the starred zeros in the cost matrix. If C(i,j) is a starred zero, then the element associated with row i is assigned to the element associated with column j.
Some of these descriptions require careful interpretation. In Step 4, for example, the possible situations are, that there is a noncovered zero which get primed and if there is no starred zero in its row the program goes onto Step 5. The other possible way out of Step 4 is that there are no noncovered zeros at all, in which case the program goes to Step 6.
*/
Munkres::Munkres(int nRow, int nCol): nRow_(nRow), nCol_(nCol), pathRow0(0), pathCol0(0), pathCount(0)
{
costs_.resize(nRow_);
for (int r = 0; r < nRow_; ++r) {
costs_[r].resize(nCol_, 0);
}
marks_.resize(nRow_);
for (int r = 0; r < nRow_; ++r) {
marks_[r].resize(nCol_); // 1 is star zero, 2 is prime zero
}
rowCover_.resize(nRow_); // record the row whether be convered
colCover_.resize(nCol_); // record the col whether be convered
augPaths_.resize(nRow_ * 2); // max save 2 in one row
for (int r = 0; r < nRow_ * 2; ++r) {
augPaths_[r].resize(2); // record the augment path
}
}
Munkres::~Munkres()
{
costs_.clear();
augPaths_.clear();
marks_.clear();
rowCover_.clear();
colCover_.clear();
}
bool Munkres::Build(const std::vector<std::vector<int> >& costs)
{
costs_ = costs;
return true;
}
float Munkres::Compute(const std::vector<std::vector<int> >& costs, std::vector<int>& pairs)
{
costs_ = costs;
bool done = false;
int step = 1;
while (!done) {
ShowCostMatrix();
ShowMaskMatrix();
std::cout << "step: " << step << std::endl;
switch (step) {
case 1:
StepOne(step);
break;
case 2:
StepTwo(step);
break;
case 3:
StepThree(step);
break;
case 4:
StepFour(step);
break;
case 5:
StepFive(step);
break;
case 6:
StepSix(step);
break;
case 7:
StepSeven(pairs);
done = true;
std::cout << "=========Done!===========" << std::endl;
break;
}
}
float minCost = 0.0;
for (int r = 0; r < pairs.size()-1; r += 2) {
minCost += costs[pairs[r]][pairs[r + 1]];
}
return minCost;
}
void Munkres::StepOne(int& step)
{
/*
For each row of the matrix, find the smallest element and subtract it from every element in its row.
Go to Step 2. We can define a local variable called minval that is used to hold the smallest value in a row.
Notice that there are two loops over the index j appearing inside an outer loop over the index i.
The first inner loop over the index j searches a row for the minval.
Once minval has been found this value is subtracted from each element of that row in the second inner loop over j.
The value of step is set to 2 just before stepone ends.
*/
int minInRow;
for (int r = 0; r < nRow_; ++r) {
minInRow = costs_[r][0];
for (int c = 1; c < nCol_; ++c) {
if (costs_[r][c] < minInRow) {
minInRow = costs_[r][c];
}
}
for (int c = 0; c < nCol_; ++c) {
costs_[r][c] -= minInRow;
}
}
step = 2;
}
void Munkres::StepTwo(int& step)
{
/*
Find a zero (Z) in the resulting matrix. If there is no starred zero in its row or column, star Z.
Repeat for each element in the matrix. Go to Step 3. In this step, we introduce the mask matrix M,
which in the same dimensions as the cost matrix and is used to star and prime zeros of the cost matrix.
If M(i,j)=1 then C(i,j) is a starred zero, If M(i,j)=2 then C(i,j) is a primed zero.
We also define two vectors R_cov and C_cov that are used to “cover” the rows and columns of the cost matrix C.
In the nested loop (over indices i and j) we check to see if C(i,j) is a zero value and if its column or row is not already covered.
If not then we star this zero (i.e. set M(i,j)=1) and cover its row and column (i.e. set R_cov(i)=1 and C_cov(j)=1).
Before we go on to Step 3, we uncover all rows and columns so that we can use the cover vectors to help us count the number of starred zeros.
*/
for (int r = 0; r < nRow_; ++r) {
for (int c = 0; c < nCol_; ++c) {
if (costs_[r][c] == 0 && rowCover_[r] == 0 && colCover_[c] == 0) {
marks_[r][c] = 1;
rowCover_[r] = 1;
colCover_[c] = 1;
}
}
}
ClearCovers();
step = 3;
}
void Munkres::StepThree(int& step)
{
/*
Cover each column containing a starred zero. If K columns are covered, the starred zeros describe a complete set of unique assignments.
In this case, Go to DONE, otherwise, Go to Step 4. Once we have searched the entire cost matrix, we count the number of independent zeros found.
If we have found (and starred) K independent zeros then we are done. If not we procede to Step 4.
*/
for (int r = 0; r < nRow_; ++r) {
for (int c = 0; c < nCol_; ++c) {
if (marks_[r][c] == 1) {
colCover_[c] = 1;
}
}
}
int colCount = 0;
for (int c = 0; c < nCol_; ++c) {
if (colCover_[c] == 1) {
colCount += 1;
}
}
if (colCount >= nCol_ || colCount >= nRow_) {
step = 7;
}
else {
step = 4;
}
}
void Munkres::StepFour(int& step)
{
/*
Find a noncovered zero and prime it. If there is no starred zero in the row containing this primed zero, Go to Step 5.
Otherwise, cover this row and uncover the column containing the starred zero. Continue in this manner until there are no uncovered zeros left.
Save the smallest uncovered value and Go to Step 6.
*/
int row = -1;
int col = -1;
bool done = false;
while (!done) {
FindAZero(row, col);
if (row == -1) {
done = true;
step = 6;
}
else {
marks_[row][col] = 2;
if (StarInRow(row)) {
FindStarInRow(row, col);
rowCover_[row] = 1;
colCover_[col] = 0;
}
else {
done = true;
step = 5;
pathRow0 = row;
pathCol0 = col;
}
}
}
}
void Munkres::StepFive(int& step)
{
/*
Construct a series of alternating primed and starred zeros as follows. Let Z0 represent the uncovered primed zero found in Step 4.
Let Z1 denote the starred zero in the column of Z0 (if any). Let Z2 denote the primed zero in the row of Z1 (there will always be one).
Continue until the series terminates at a primed zero that has no starred zero in its column.
Unstar each starred zero of the series, star each primed zero of the series, erase all primes and uncover every line in the matrix. Return to Step 3.
You may notice that Step 5 seems vaguely familiar.
It is a verbal description of the augmenting path algorithm (for solving the maximal matching problem) which we discussed in Lecture 3.
We decompose the operations of this step into a main procedure and five relatively simple subprograms.
*/
int r = -1;
int c = -1;
pathCount = 1;
augPaths_[pathCount - 1][0] = pathRow0;
augPaths_[pathCount - 1][1] = pathCol0;
bool done = false;
while (!done) {
FindStarInCol(augPaths_[pathCount - 1][1], r);
if (r > -1) {
pathCount += 1;
augPaths_[pathCount - 1][0] = r;
augPaths_[pathCount - 1][1] = augPaths_[pathCount - 2][1];
}
else {
done = true;
}
if (!done) {
FindPrimeInRow(augPaths_[pathCount - 1][0], c);
pathCount += 1;
augPaths_[pathCount - 1][0] = augPaths_[pathCount - 2][0];
augPaths_[pathCount - 1][1] = c;
}
}
AugmentPath();
ClearCovers();
ErasePrimes();
step = 3;
}
void Munkres::StepSix(int& step)
{
/*
Add the value found in Step 4 to every element of each covered row, and subtract it from every element of each uncovered column.
Return to Step 4 without altering any stars, primes, or covered lines. Notice that this step uses the smallest uncovered value in the cost matrix to modify the matrix.
Even though this step refers to the value being found in Step 4 it is more convenient to wait until you reach Step 6 before searching for this value.
It may seem that since the values in the cost matrix are being altered, we would lose sight of the original problem.
However, we are only changing certain values that have already been tested and found not to be elements of the minimal assignment.
Also we are only changing the values by an amount equal to the smallest value in the cost matrix,
so we will not jump over the optimal (i.e. minimal assignment) with this change.
*/
float minVal = std::numeric_limits<float>::max();
FindSmallest(minVal);
for (int r = 0; r < nRow_; ++r) {
for (int c = 0; c < nCol_; ++c) {
if (rowCover_[r] == 1) {
costs_[r][c] += minVal;
}
if (colCover_[c] == 0) {
costs_[r][c] -= minVal;
}
}
}
step = 4;
}
void Munkres::StepSeven(std::vector<int>& pairs)
{
for (int r = 0; r < nRow_; ++r) {
pairs.emplace_back(r);
int c = -1;
FindStarInRow(r, c);
pairs.emplace_back(c);
}
}
void Munkres::FindAZero(int& row, int& col)
{
int r = 0;
int c;
bool done = false;
row = -1;
col = -1;
while (!done) {
c = 0;
while (true) {
if (costs_[r][c] == 0 && rowCover_[r] == 0 && colCover_[c] == 0) {
row = r;
col = c;
done = true;
}
c += 1;
if (c >= nCol_ || done) {
break;
}
}
r += 1;
if (r >= nRow_) {
done = true;
}
}
}
bool Munkres::StarInRow(int row)
{
bool tmp = false;
for (int c = 0; c < nCol_; ++c) {
if (marks_[row][c] == 1) {
tmp = true;
break;
}
}
return tmp;
}
void Munkres::FindStarInRow(int row, int& col)
{
col = -1;
for (int c = 0; c < nCol_; ++c) {
if (marks_[row][c] == 1) {
col = c;
break;
}
}
}
void Munkres::FindStarInCol(int col, int& row)
{
row = -1;
for (int r = 0; r < nRow_; ++r) {
if (marks_[r][col] == 1) {
row = r;
break;
}
}
}
void Munkres::FindPrimeInRow(int row, int& col)
{
for (int c = 0; c < nCol_; ++c) {
if (marks_[row][c] == 2) {
col = c;
break;
}
}
}
void Munkres::AugmentPath()
{
for (int p = 0; p < pathCount; ++p) {
if (marks_[augPaths_[p][0]][augPaths_[p][1]] == 1) {
marks_[augPaths_[p][0]][augPaths_[p][1]] = 0;
}
else {
marks_[augPaths_[p][0]][augPaths_[p][1]] = 1;
}
}
}
void Munkres::ClearCovers()
{
for (int r = 0; r < nRow_; ++r) {
rowCover_[r] = 0;
}
for (int c = 0; c < nCol_; ++c) {
colCover_[c] = 0;
}
}
void Munkres::ErasePrimes()
{
for (int r = 0; r < nRow_; ++r) {
for (int c = 0; c < nCol_; ++c) {
if (marks_[r][c] == 2) {
marks_[r][c] = 0;
}
}
}
}
void Munkres::FindSmallest(float& minVal)
{
for (int r = 0; r < nRow_; ++r) {
for (int c = 0; c < nCol_; ++c) {
if (rowCover_[r] == 0 && colCover_[c] == 0) {
if (minVal > costs_[r][c]) {
minVal = costs_[r][c];
}
}
}
}
}
void Munkres::ShowCostMatrix()
{
std::cout << "-----Costs-----" << std::endl;
for (int r = 0; r < nRow_; ++r) {
for (int c = 0; c < nCol_; ++c) {
std::cout << costs_[r][c] << " ";
}
std::cout << std::endl;
}
}
void Munkres::ShowMaskMatrix()
{
std::cout << "*****Marks*****" << std::endl;
for (int r = 0; r < nRow_; ++r) {
for (int c = 0; c < nCol_; ++c) {
std::cout << int(marks_[r][c]) << " ";
}
std::cout << std::endl;
}
std::cout << "===========================" << std::endl;
}
int main()
{
int nRow = 3;
int nCol = 3;
std::vector<std::vector<int> > costs { {1, 2, 3}, {2, 4, 6}, {3, 6, 9} };
Munkres munkres = Munkres(nRow, nCol);
//munkres.Build(costs);
std::vector<int> pairs;
float cost = munkres.Compute(costs, pairs);
std::cout << "cost: " << cost << std::endl;
for (int i = 0; i < pairs.size() - 1; i += 2) {
std::cout << pairs[i] << "==>" << pairs[i + 1] << std::endl;
}
return 0;
}
