Kuhn-Munkres算法

发表于 更新于 分类于 阅读次数： Valine：
本文字数： 5.4k 阅读时长 ≈ 5 分钟

Kuhn-Munkres算法

本文接着上一篇博客《匈牙利算法》，前篇讨论的是不带权的二分图求最大匹配，本篇讨论带权的二分图求最大权匹配——KM算法。实际上前者可以看作是后者的特殊情况，即权值只有0和1的情形。

几个基本事实

  • WLOG，通过人为添加权值为0的边，任一个二分图\(G\)都可以被改造成完全图而不影响我们所讨论的问题，因此以下的讨论全都视作完全二分带权图（只考虑非负权值）。

  • 如前篇所述，可以通过取反增广路的操作使得匹配边数+1

  • 最大权匹配一定是完备匹配（如果最大权匹配不完备的话，说明\(X\)和\(Y\)至少各有一个点未被匹配，加上这条边之后权和不可能更小）

可行顶标与相等子图

顶标：顶标是一个结点函数\(l:V\to \mathbb{R}\)

可行顶标：可行顶标是满足下面不等式的顶标（即边权不能超过两端点的顶标和） \[ l(x)+l(y)\ge w(x,y), \forall x \in X, y \in Y \] 相等子图：图\(G=(V,E)\)的（相对于顶标\(l\)的）相等子图\(G_l=(V,E_l)\)包含图\(G\)的所有顶点，但只包含边权等于两端点顶标和的边。即 \[ E_l={(x,y):l(x)+l(y)=w(x,y)} \]

Kuhn-Munkres定理：如果\(l\)是可行顶标，\(M\)是相等子图\(G_l\)里的完美匹配，那么\(M\)是图\(G\)的最大权匹配。

证明：用\(e=(e_x,e_y)\)来记图中的边，\(e_x\)和\(e_y\)分别是两端点。任取图\(G\)的一个完美匹配\(M'\)，容易的到下面的不等式 \[ w(M')=\sum_{e\in M'} w(e) \le \sum_{e\in M'} l(e_x)+l(e_y) = \sum_{v\in V} l(v) \] 这个不等式说明任何一个完美匹配权和的上界就是所有顶标和\(\sum_{v\in V} l(v)\)，现在\(M\)是相等子图\(G_l\)的完美匹配，从而也是原图\(G\)的完美匹配，所以有 \[ w(M)=\sum_{e\in M} w(e) = \sum_{v\in V} l(v) \] 最大权匹配一定是一个完备匹配，易知所有完备匹配的权和上界也是\(\sum_{v\in V}l(v)\)，所以\(M\)就是最大权匹配。

  • KM定理将求最大权匹配问题转化成了求完美匹配问题，这是组合优化问题里经典的技巧

  • 注意KM定理的证明里实际上证明了任意一个匹配\(M\)和任意一个可行顶标\(l\)都满足这个不等式

\[ w(M)\le \sum_{v\in V}l(v) \]

  • 和最大流最小割定理有哪些相似之处？

KM算法

算法框架

  1. 给定初始可行顶标\(l\)和初始匹配\(M\)
  2. 当\(M\)不是\(E_l\)的完美匹配时，重复以下操作
    1. 在相等子图\(E_l\)里寻找增广路，然后通过取反操作使匹配数+1
    2. 如果找不到增广路，将可行顶标\(l\)改进为\(l'\)使得\(E_l \subset E_{l'}\)，跳转至1
  • 如果\(G\)有完美匹配，因为每次迭代要么匹配数+1要么相等子图扩大，最终算法一定会停止，而且停止时\(M\)是\(E_l\)的完美匹配。根据KM定理，\(M\)就是\(G\)的最大权匹配。

初始可行顶标和初始匹配

最简单的可行顶标就是取\(X\)的顶标为邻接边的最大权值，\(Y\)的顶标取0 \[ \forall y \in Y, l(y) = 0 \\ \forall x \in X, l(x) = \max_{y\in Y}\{w(x,y)\} \] 这样显然满足可行顶标的要求 \[ \forall x \in X, y \in Y, w(x) \le l(x) + l(y) \] 最简单的初始匹配就是\(M=\phi\)

改进可行顶标

设\(l\)是一个可行顶标，顶点\(u\in V\)，点集\(S\subseteq V\)，定义点\(u\)的邻接集\(N_l(u)\)和\(S\)的邻接集\(N_l(S)\) \[ N_l(u) = \{v:(u,v)\in E_l\} \\ N_l(S) = \bigcup_{u\in S} N_l(u) \]

引理：设\(S\subseteq X, T=N_l(S)\ne Y\)，令\(\alpha_l\)为

\[ \alpha_l = \min_{x\in S,y\notin T} \{l(x)+l(y)-w(x,y)\} \] 令顶标\(l'\)为 \[ l'(v) = \begin{cases} l(v)-\alpha_l, & v\in S \\ l(v)+\alpha_l, & v\in T \\ l(v), & otherwise \end{cases} \] 那么\(l'\)也是一个可行顶标，并且

  1. 若\((x,y)\in E_l,x\in S,y\in T\)，那么\((x,y)\in E_{l'}\)
  2. 若\((x,y)\in E_l,x\notin S,y\notin T\)，那么\((x,y)\in E_{l'}\)
  3. 满足\(x\in S, y\notin T\)的边\((x,y)\)可能属于新的相等子图\(E_{l'}\)
  • 这个引理说明：通过这种方式改进可行顶标，能够使原先在相等子图里的边仍然在相等子图里，有一端在\(S\)里，另一端不在\(T\)里的边有可能被纳入新的相等子图

伪代码

  1. 初始化可行顶标\(l\)和相等子图\(E_l\)的初始匹配\(M\) \[ \forall y \in Y, l(y) = 0 \\ \forall x \in X, l(x) = \max_{y\in Y}\{w(x,y)\} \\ M=\phi \]

  2. 若\(M\)已经是完美匹配，停止算法；否则，选取未匹配点\(u\in X\)，令\(S = \{u\}， T=\phi\)

  3. 若\(N_l(S) = T\)，更新可行顶标 \[ \alpha_l = \min_{x\in S,y\notin T} \{l(x)+l(y)-w(x,y)\} \\ l'(v) = \begin{cases} l(v)-\alpha_l, & v\in S \\ l(v)+\alpha_l, & v\in T \\ l(v), & otherwise \end{cases} \]

  4. 若\(N_l(S)\ne T\)，选取\(y\in N_l(S) - T\)

    1. 若\(y\)是未匹配点，那么从\(u\)到\(y\)就找到了一条增广路，取反使匹配数+1，跳转至2
    2. 若\(y\)是匹配点，设\(y\)目前与\(z\)匹配，令\(S=S\cup \{z\},T=T\cup \{y\}\)，跳转至3

算法正确性

由引理保证了算法在执行过程一致保持着顶标的可行性，算法通过改变可行顶标从而逐渐扩大相等子图。每次找到一个\(Y\)里的未匹配点，就使得匹配数+1。最终算法会在达到完美匹配的时候终止，根据KM定理，就得到了原图的最大权匹配。

时间复杂度分析

在具体实现中，通常为每个\(y\)结点定义一个松弛度存放在数组里，这样可以减少\(\alpha_l\)的计算量 \[ slack_y = \min_{x\in S} \{l(x) + l(y) - w(x,y)\} \] 首先将算法运行的过程分为\(|M|\)个阶段，使匹配数+1的若干次迭代记作一个阶段，每个阶段都经历了若干次步骤2、3、4，下面来具体分析每个步骤的时间复杂度。

步骤2：每个阶段都只会经历一次步骤2。初始化\(S\)和\(T\)显然是\(O(1)\)的，但此时还要初始化\(slack\)数组。因为\(S\)里只有1个点，只需要\(|S|\)次计算就可初始化完毕，所以步骤2总共是\(O(|V|)\)的。

步骤3：计算\(\alpha_l=\min_{y\in T} slack_y\)这一步是\(O(|T|)\)的，更新可行顶标\(l'\)是\(O(|T|+|S|)\)的，更新完可行顶标之后还要重新更新\(slack\)数组，根据定义不难推出这里只需要执行\(\forall y\notin T, slack_y = slack_y - \alpha_l\)，所以是\(O(|Y-T|)\)的。要注意的是步骤3在每个阶段可能经历多次，但最多只会经历\(O(|Y|)\)次（因为每经历一次步骤4如果需要回到步骤3的话，一定使\(T\)增加一个点）。上面的复杂度全都可以放缩成\(O(|V|)\)，所以步骤3总共是\(O(|V|^2)\)的。

步骤4：第1种情况在每个阶段只会经历一次，取反操作是\(O(|V|)\)的。第2种情况同步骤3的分析，至多经历\(O(|V|)\)次，更新集合\(S\)和\(T\)只需要\(O(1)\)的时间。但由于集合\(S\)多加了一个点，需要更新\(slack\)数组，\(\forall y\notin T, slack_y = \min \{slack_y, l(z)+l(y)-w(z,y)\}\)，这需要\(O(|Y|)\)的时间。步骤4总共是\(O(|V|^2)\)的。

最后，因为匹配数\(|M|\)不可能超过\(|V|/2\)，所以阶段数也是\(O(|V|)\)的，得出整个算法的时间复杂度是\(O(|V|^3)\)

cpp实现

由于时间紧迫，用松弛度数组优化后的\(O(|V|^3)\)的算法暂时先挖个坑以后再填。这里给出原始的\(O(|V|^4)\)实现。

用邻接矩阵存储二分图，两个数组存储顶标，match数组用来记录匹配（记录增广路），S和T两个数组来记录算法里的两个集合\(S\)和\(T\)

1
2
3
4
int Weight[maxm][maxn];
int Lx[maxm], Ly[maxn]; // 顶标
int match[maxn];    // 记录匹配
bool S[maxm], T[maxn];  // 算法中的两个集合S和T

步骤1的初始化可行顶标和初始化匹配写成一个函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
void Init()
{
    // 将X集合的顶标设为最大边权，Y集合的顶标设为0
    for (int i = 1; i <= m; i++)
    {
        Lx[i] = 0;
        for (int j = 1; j <= n; j++)
        {
            match[j] = 0;   // match记录的是Y集合里的点与谁匹配
            Ly[j] = 0;
            Lx[i] = max(Lx[i], Weight[i][j]);
        }
    }
}

步骤3的更新顶标就直接用没优化前的原始实现了

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
void update() // 更新顶标
{
    // 计算a
    int a = 1 << 30;
    for (int i = 1; i <= m; i++)
        if (S[i])
            for (int j = 1; j <= n; j++)
                if (!T[j])
                    a = min(a, Lx[i] + Ly[j] - Weight[i][j]);

    // 修改顶标
    for (int i = 1; i <= m; i++)
        if (S[i])
            Lx[i] -= a;
    for (int j = 1; j <= n; j++)
        if (T[j]) 
            Ly[j] += a;
}

步骤4可以用递归来实现，写成了DFS（这里和前篇的匈牙利算法实现非常相似，因为都是在寻找增广路）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
bool findPath(int i)
{
    S[i] = true;
    for (int j = 1; j <= n; j++)
    {
        if (Lx[i] + Ly[j] == Weight[i][j] && !T[j]) // 找出在相等子图里又还未被标记的边
        {
            T[j] = true;
            if (!match[j] || findPath(match[j])) // 未被匹配，或者已经匹配又找到增广路
            {
                match[j] = i;
                return true;
            }
        }
    }
    return false;
}

然后就是整个算法的合写

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
void KM()
{
    Init();

    for (int i = 1; i <= m; i++)
    {
        while (true)
        {
            for (int i = 1; i <= m; i++)
                S[i] = 0;
            for (int j = 1; j <= n; j++)
                T[j] = 0;
            if (!findPath(i))
                update();
            else
                break;
        }
    }
}

最后调用即可，完成后查看match数组就能查看匹配结果

参考文献

受益匪浅的一篇博客

另一篇受益匪浅的博客

又一篇受益匪浅的博客

形象描述算法运行过程

线性规划的观点值得一看

一篇介绍的非常详细的博客

本文撰写时的主要参考来源

一篇Murray State University的KM算法教程

最后放上Kuhn和Munkres的论文

[1] Munkres, James. "Algorithms for the assignment and transportation problems." Journal of the society for industrial and applied mathematics 5.1 (1957): 32-38.

[2] Kuhn, Harold W. "The Hungarian method for the assignment problem." Naval research logistics quarterly 2.1‐2 (1955): 83-97.