匈牙利算法解决加权二分图问题

匈牙利方法是一种组合优化算法，它在多项式时间内解决了赋值问题，广泛应用于多目标跟踪的关联问题中。

图1：（a）二分图，（b）边权重矩阵，（c）边成本的替代表示形式

动机：分配问题

假设有 $n$n 辆卡车每辆可装载一种产品以及 $n$n 家商店，这些商店愿意以矩阵 $W$W 代表的不同价格购买 $n$n 种不同的产品。分配问题：若商店 $y_i$yi​ 提出以 $W_{ij}$Wij​ 美元从卡车 $x_i$xi​ 购买产品，我们如何指派每辆卡车 $x_i$xi​ 去商店 $y_j$yj​，以便在所有可能的任务中最大化合并利润？

1. 一般问题：给定 $X = {x_1, \dots , x_n}$X=x1​,…,xn​，$Y = {y_1, \dots , y_n}$Y=y1​,…,yn​，矩阵 $W$W 中 $W_{ij} = \mathrm{weight}(x_i, y_j )$Wij​=weight(xi​,yj​) 是将 $x_i$xi​ 分配给 $y_i$yi​ 的权重，找到将每个 $x_i$xi​ 分配给某个 $y_j$yj​ 的匹配，使得总权重最大化。

2. 朴素的算法为：遍历所有 $n!$n! 种可能的分配，选择得分最高的。

3. 假设：$\forall i, j \in {1, \dots , n} : Wij \geq 0$∀i,j∈1,…,n:Wij≥0。

4. 可以将问题视为完全加权的二分图 $G = (V, E)$G=(V,E)：
   $\begin{aligned} &V = X \cup Y \\ &E = {(x_i, y_j )}_{x_i\in X,y_j\in Y} \\ &\mathrm{weight}(x_i, y_j) = W_{ij} \end{aligned}$​V=X∪YE=(xi​,yj​)xi​∈X,yj​∈Y​weight(xi​,yj​)=Wij​​

5. 分配是完美匹配：问题简化为最大化权重，找到完美匹配。

匈牙利算法的基础（Kuhn-Munkres）

定义

1. 图 $G = (V, E)$G=(V,E) 的标记是函数 $l : V \rightarrow R$l:V→R，这样：
   $\forall(u, v) \in E : l(u) + l(v) \geq \mathrm{weight}((u, v))$∀(u,v)∈E:l(u)+l(v)≥weight((u,v))

2. 相等子图为子图 $G_l = (V, E_l) \subseteq G = (V, E)$Gl​=(V,El​)⊆G=(V,E)，固定标签函数 $l$l，则

$E_l = {(u, v) \in E : l(u) + l(v) = \mathrm{weight}((u, v))}$El​=(u,v)∈E:l(u)+l(v)=weight((u,v))

Kuhn-Munkres 定理

定理2.1：给定标记 $l$l，如果 $M$M 是 $G_l$Gl​ 上的完美匹配，则 $M$M 是 $G$G 的最大权重匹配。

1. 假设 $M_0$M0​ 是 $G$G 中的任意完美匹配。通过定义标记函数，由于 $M_0$M0​ 是完美的，
   $\mathrm{weight}(M') =\sum_{(u,v)\in M'} \mathrm{weight}((u, v)) \leq \sum_{(u,v)\in M'}l(u) + l(v) = \sum_{v \in V}l(v)$weight(M′)=(u,v)∈M′∑​weight((u,v))≤(u,v)∈M′∑​l(u)+l(v)=v∈V∑​l(v)
2. 这意味着：$\sum_{v\in V}l(v)$∑v∈V​l(v) 是 $G$G 的任何完美匹配 $M'$M′ 的上界。
3. 现在看看匹配 $M$M 的权重：
   $\mathrm{weight}(M) = \sum_{(u,v)\in M} \mathrm{weight}((u, v)) = \sum_{(u,v)\in M}l(u) + l(v)= \sum_{v \in V}l(v)$weight(M)=(u,v)∈M∑​weight((u,v))=(u,v)∈M∑​l(u)+l(v)=v∈V∑​l(v)
4. 通过1. 和3. 中等式可得 $G$G 的所有完美匹配 $M'$M′ 满足：
   $\mathrm{weight}(M) \geq \mathrm{weight}(M')$weight(M)≥weight(M′)

关键点：通过 Kuhn-Munkres 定理，找到最大权重分配的问题简化为找到正确的标记函数和相应的相等子图上的任何完美匹配。

增加匹配

给定标记 $l$l，$G_l= (V,E_l)$Gl​=(V,El​)，$G_l$Gl​ 上的一些匹配 $M$M，未匹配的顶点 $u\in V, u\notin M$u∈V,u∈/​M。

1. 如果路径在 $E_l-M$El​−M 和 $M$M 之间交替，并且路径的第一个和最后一个顶点在 $M$M 中未匹配，则该路径是 $M$M 在 $G_l$Gl​ 上的增广路径。记录从 $u$u 开始的“近似”增广路径。
2. 如果我们能找到一个不匹配的顶点 $v$v，那么我们创建从 $u$u 到 $v$v 的增广路径 $\alpha$α。
3. 通过将 $M$M 中的边替换为增广路径中 $E_l-M$El​−M 中的边来翻转匹配。
4. 由于起始和终止顶点未匹配，这增加了匹配的大小 $\implies | M'|> | M |$⟹∣M′∣>∣M∣

改进标签

1. $S\subseteq X$S⊆X，$T\subseteq Y$T⊆Y 且 $S$S、$T$T “几乎”代表当前匹配 $M$M 与外部其他边 $E_{l-M}$El−M​ 之间的增广交替路径。
2. Let Nl(S) be the neighbors to each node in S along El. Nl(S) = {v|∀u ∈ S : (u, v) ∈ El}
   设 $N_l(S)$Nl​(S) 为沿 $E_l$El​ 的 $S$S 中每个节点的邻居。
   $N_l(S) = \{ v\mid \forall u \in S : (u, v) \in E_l\}$Nl​(S)={v∣∀u∈S:(u,v)∈El​}
3. 如果 $N_l(S) = T$Nl​(S)=T 我们不能增加交替路径并增广，所以我们必须改进标签！
4. 计算：$\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))}
5. 改进 $l\rightarrow l'$l→l′：
   $\begin{aligned} l'(r) = \begin{cases} l(r) − \delta_l &\text{if } r \in S \\ l(r) + \delta_l &\text{if } r \in T \\ l(r) &\text{otherwise } \end{cases} \end{aligned}$l′(r)=⎩⎪⎨⎪⎧​l(r)−δl​l(r)+δl​l(r)​if r∈Sif r∈Totherwise ​​
6. 声明：$l'$l′ 是有效标签且 $E_l\subset E_{l'}$El​⊂El′​。
7. 通过检验元素 $u \in S, u \notin S, v \in T, v \notin T$u∈S,u∈/​S,v∈T,v∈/​T 所有可能的情况证明。

Kuhn-Munkres 算法

1. 从一些匹配 $M$M 开始，并且有效标记
   $l ::= \forall x \in X, y \in Y : l(y) = 0, l(x) = \max_{y'\in Y} (\text{weight}(x, y'))$l::=∀x∈X,y∈Y:l(y)=0,l(x)=y′∈Ymax​(weight(x,y′))
2. 执行以下操作直到 $M$M 完美匹配：
   (a) 寻找增广路径
   (b) 如果不存在增广路径，则改进 $l\rightarrow l'$l→l′ 并转到步骤 (a)。

复杂度

1. 每个步骤(a)或(b)每轮增加1个匹配的边，因此总轮数为 $O(|V | = 2n) = O(n)$O(∣V∣=2n)=O(n)。
2. 增加 $M$M：找到正确的顶点（如果存在）需要 $O(|V |)$O(∣V∣) ，翻转匹配为 $O(|V |)$O(∣V∣)。
3. 改进标签：找到 $\delta_l$δl​ 并更新标签的计算量为 $O(|V |)$O(∣V∣) 。如果没有找到增广路径，则改进标签可能发生 $O(|V |)$O(∣V∣) 次。所以在一轮中总共为 $O(|V |^2)$O(∣V∣2)。
4. $O(|V|)$O(∣V∣) 轮次且每次执行 $O(|V |^2)$O(∣V∣2)，总运行时间为 $O(|V|^3) = O(n^3)$O(∣V∣3)=O(n3)。

max_cost_assignment

C++SORT 使用 Dlib 中的 max_cost_assignment 函数完成关联匹配。

输入参数cost_是只读的，所以复制一份。
检查矩阵元素是否为整型。

        const_temp_matrix<EXP> cost(cost_);
        typedef typename EXP::type type;
        // This algorithm only works if the elements of the cost matrix can be reliably 
        // compared using operator==. However, comparing for equality with floating point
        // numbers is not a stable operation. So you need to use an integer cost matrix.
        COMPILE_TIME_ASSERT(std::numeric_limits<type>::is_integer);
        DLIB_ASSERT(cost.nr() == cost.nc(),
            "\t std::vector<long> max_cost_assignment(cost)"
            << "\n\t cost.nr(): " << cost.nr()
            << "\n\t cost.nc(): " << cost.nc()
            );

参考链接已打不开，算法复杂度为 $O(n^3)$O(n3)。

        using namespace dlib::impl;
        /*
            I based the implementation of this algorithm on the description of the
            Hungarian algorithm on the following websites:
                http://www.math.uwo.ca/~mdawes/courses/344/kuhn-munkres.pdf
                http://www.topcoder.com/tc?module=Static&d1=tutorials&d2=hungarianAlgorithm
            Note that this is the fast O(n^3) version of the algorithm.
        */

        if (cost.size() == 0)
            return std::vector<long>();

lx即 $l(x)$l(x)，ly即 $l(y)$l(y)。xy和yx表示 $G_l$Gl​ 的完美匹配，即匹配 $M$M。$S \subseteq X, T\subseteq Y$S⊆X,T⊆Y，$S, T$S,T 表示匹配 $M$M 和其他边 $E_l-M$El​−M 之间的当前“近似”增广交替路径。

        std::vector<type> lx, ly;
        std::vector<long> xy;
        std::vector<long> yx;
        std::vector<char> S, T;
        std::vector<type> slack;
        std::vector<long> slackx;
        std::vector<long> aug_path;

初始时 $M = \emptyset$M=∅。

        // Initially, nothing is matched. 
        xy.assign(cost.nc(), -1);
        yx.assign(cost.nc(), -1);
        /*
            We maintain the following invariant:
                Vertex x is matched to vertex xy[x] and
                vertex y is matched to vertex yx[y].
                A value of -1 means a vertex isn't matched to anything.  Moreover,
                x corresponds to rows of the cost matrix and y corresponds to the
                columns of the cost matrix.  So we are matching X to Y.
        */

$\mathrm{weight}(M') =\sum_{(u,v)\in M'} \mathrm{weight}((u, v)) \leq \sum_{(u,v)\in M'}l(u) + l(v) = \sum_{v \in V}l(v)$weight(M′)=(u,v)∈M′∑​weight((u,v))≤(u,v)∈M′∑​l(u)+l(v)=v∈V∑​l(v)

        // Create an initial feasible labeling.  Moreover, in the following
        // code we will always have: 
        //     for all valid x and y:  lx[x] + ly[y] >= cost(x,y)
        lx.resize(cost.nc());
        ly.assign(cost.nc(),0);
        for (long x = 0; x < cost.nr(); ++x)
            lx[x] = max(rowm(cost,x));

逐顶点搜索。队列q用于存储 BFS 的变量。
每次搜索前重置 $S, T$S,T，以及slack、slackx。
slack存储顶点可行的标签。S和T记录是否选中顶点。

        // Now grow the match set by picking edges from the equality subgraph until
        // we have a complete matching.
        for (long match_size = 0; match_size < cost.nc(); ++match_size)
        {
            std::deque<long> q;

            // Empty out the S and T sets
            S.assign(cost.nc(), false);
            T.assign(cost.nc(), false);

            // clear out old slack values
            slack.assign(cost.nc(), std::numeric_limits<type>::max());
            slackx.resize(cost.nc());
            /*
                slack and slackx are maintained such that we always
                have the following (once they get initialized by compute_slack() below):
                    - for all y:
                        - let x == slackx[y]
                        - slack[y] == lx[x] + ly[y] - cost(x,y)
            */

            aug_path.assign(cost.nc(), -1);

遍历列，尝试找到一个未匹配的顶点x，调用 compute_slack 函数。这里是不是应该使用cost.nr()？

            for (long x = 0; x < cost.nc(); ++x)
            {
                // If x is not matched to anything
                if (xy[x] == -1)
                {
                    q.push_back(x);
                    S[x] = true;

                    compute_slack(x, slack, slackx, cost, lx, ly);
                    break;
                }
            }

从q中取出一个x，尝试找到一个未匹配的y。

            long x_start = 0;
            long y_start = 0;

            // Find an augmenting path.  
            bool found_augmenting_path = false;
            while (!found_augmenting_path)
            {
                while (q.size() > 0 && !found_augmenting_path)
                {
                    const long x = q.front();
                    q.pop_front();
                    for (long y = 0; y < cost.nc(); ++y)
                    {
                        if (cost(x,y) == lx[x] + ly[y] && !T[y])
                        {
                            // if vertex y isn't matched with anything
                            if (yx[y] == -1) 
                            {
                                y_start = y;
                                x_start = x;
                                found_augmenting_path = true;
                                break;
                            }

对于已匹配的y，沿其路继续搜寻。

                            T[y] = true;
                            q.push_back(yx[y]);

                            aug_path[yx[y]] = x;
                            S[yx[y]] = true;
                            compute_slack(yx[y], slack, slackx, cost, lx, ly);
                        }
                    }
                }

                if (found_augmenting_path)
                    break;

未找到增广路径，改进 $l\rightarrow l'$l→l′：
$\begin{aligned} l'(r) = \begin{cases} l(r) − \delta_l &\text{if } r \in S \\ l(r) + \delta_l &\text{if } r \in T \\ l(r) &\text{otherwise } \end{cases} \end{aligned}$l′(r)=⎩⎪⎨⎪⎧​l(r)−δl​l(r)+δl​l(r)​if r∈Sif r∈Totherwise ​​

                // Since we didn't find an augmenting path we need to improve the 
                // feasible labeling stored in lx and ly.  We also need to keep the
                // slack updated accordingly.
                type delta = std::numeric_limits<type>::max();
                for (unsigned long i = 0; i < T.size(); ++i)
                {
                    if (!T[i])
                        delta = std::min(delta, slack[i]);
                }
                for (unsigned long i = 0; i < T.size(); ++i)
                {
                    if (S[i])
                        lx[i] -= delta;

                    if (T[i])
                        ly[i] += delta;
                    else
                        slack[i] -= delta;
                }

改进标签后，如果松弛量为0的边未匹配，则找到增广路径；否则，继续沿其搜索。

                q.clear();
                for (long y = 0; y < cost.nc(); ++y)
                {
                    if (!T[y] && slack[y] == 0)
                    {
                        // if vertex y isn't matched with anything
                        if (yx[y] == -1)
                        {
                            x_start = slackx[y];
                            y_start = y;
                            found_augmenting_path = true;
                            break;
                        }
                        else
                        {
                            T[y] = true;
                            if (!S[yx[y]])
                            {
                                q.push_back(yx[y]);

                                aug_path[yx[y]] = slackx[y];
                                S[yx[y]] = true;
                                compute_slack(yx[y], slack, slackx, cost, lx, ly);
                            }
                        }
                    }
                }
            } // end while (!found_augmenting_path)

沿着增广路径反转边。

            // Flip the edges along the augmenting path.  This means we will add one more
            // item to our matching.
            for (long cx = x_start, cy = y_start, ty; 
                 cx != -1; 
                 cx = aug_path[cx], cy = ty)
            {
                ty = xy[cx];
                yx[cy] = cx;
                xy[cx] = cy;
            }

        }


        return xy;

compute_slack

compute_slack 函数对slack进行更新，计算顶点x的松弛量并记录对应的y。
$\delta_l= \min_{u\in S,v \notin T }{l(u) +l(v)−weight((u,v))}$δl​=u∈S,v∈/​Tmin​l(u)+l(v)−weight((u,v))

            for (long y = 0; y < cost.nc(); ++y)
            {
                if (lx[x] + ly[y] - cost(x,y) < slack[y])
                {
                    slack[y] = lx[x] + ly[y] - cost(x,y);
                    slackx[y] = x;
                }
            }

参考资料：

• 二分图的最大匹配、完美匹配和匈牙利算法
• The Optimal Assignment Problem
• assignment-problem-and-hungarian-algorithm
• The Hungarian Algorithm for Weighted Bipartite Graphs
• 匈牙利算法详解（含时间复杂度）
• The Hungarian algorithm: Kuhn-Munkres theorem