序
我觉得我现在还可以
安利一下 rsx's blog
正文
本文权当作学习二分图匹配的一些笔记,个人风格极为严重,请勿深究语文
匈牙利
本文不介绍算法的步骤,因为可以 google(逃,写的都是一些自己的总结
匈牙利的算法主要功能是可以做二分图的最大匹配。它是一种交替的找增广路的算法,其本质就是逆流修改(Singercoder 说的)
其核心的就是 re 这个数组不会回溯 (当然每次要清)
为什么呢,其实我觉得和当前当前弧优化有点联系。你想想,他要是回溯的话,都不是多项式算法了(笑)。但是其算法正确性也是有的。re 表示的是这个结点是否寻找过增广路。
理论上最坏时间复杂度会到 (O(VE)) 还行吧。一般卡不了。
算法伪代码
main()
{
last...
for 1 to n
fill re with 0
/* 每一次都要清空 re */
if find(i) ans++;
next...
}
find(u)
{
for any edge of u
v = e.v;
/* 与暴力算法的区别,和当前弧异曲同工吧 */
if re[v]
continue
/* 打上标记 */
re[v] = 1
if p[v] == NULL or find(p[v])
/* 逻辑要清晰 */
p[v] = u
return true
return false
}
写的不是很好,但是还算清楚吧。毕竟代码实现的能力也很重要的(误)。
然后通过这个操作,就可以求出二分图的最大匹配了。ヾ(≧∇≦*)ゝ
当然了,谈到二分图匹配自然会有 Singercoder 冒出来,你匈牙利是 (O(nm)) ,我最大流跑是 (O(sqrt n m)) (没有严格证明,听dalao说的),客观的说,其实各有千秋吧。
- 运行时间:几乎相等,因为没有毒瘤可以把匈牙利卡到 (O(nm)),而且 ISAP 的常数肯定比匈牙利大呀.
- 代码实现难度:匈牙利完爆各种最大流
- 写完的成就感:最大流完爆匈牙利
- 运行空间:会卡这个?
真正的代码
bool find(int u) {
for (int i = h[u]; i != -1; i = edge[i].lac)
{
int to = edge[i].to;
if (re[to]) continue;
re[to] = 1;
if (!p[to] || find(p[to])) { p[to] = u; return 1; }
}
return 0;
}
然后我们就可以用它来搞各种与二分图有关的东西。
下文介绍几种别的扩展
- 最少点覆盖
- 最少边覆盖
- 最大独立子集
- 最少不相交路径覆盖
最小点覆盖
反证一下,他应该就是最大匹配数,因为如果有边并未覆盖。那么他的左右端点应该都没被选,所以就不是最大匹配了,矛盾
所以:最小点覆盖 = 最大匹配数
模板:poj1325(有坑)
最少边覆盖
这么思考,贪心的想,先把所有匹配的边选了,然后对于剩下的点随便选一个边
于是:最少边覆盖 = 点数 - 最大匹配数
模板:poj3020
最大独立子集
理解为一个删点和与其相连的所有边,那么剩下的点就是独立子集。删去最小点覆盖即可,
于是最大独立子集 = 点数 - 最小点覆盖
模板:poj1466(据说数据有坑但是没有掉进去...但是还是很坑
最少不相交路径覆盖
ps:我实现这个一般都是网络流..因为 luogu 的模板要输出路径...
这么想,我们最开始有 n 个不相交路径(n 个点),然后要进行连边,连尽可能多的边,就可以减少路径数...
然后,把一个点分为前点,后点做二分图匹配即可。
最少路径数 = 点数 - 最大匹配数
模板:poj1422
关于算法选择的闲话
我们可以选择匈牙利来跑,因为 (O(nm)) 的时间复杂度的确可以接受,如果觉得不保险其实跑一个 isap 也未尝不可,所以我还是更喜欢最大流。
KM算法
KM算法是可以解决完备匹配下的最大权匹配,个人认为是对匈牙利进行了一番操作)
完备匹配
完备匹配就是要求匹配数为 n
书面的说:所谓的完备匹配就是在二部图中,x 点集中的所有点都有对应的匹配且 y 点集中所有的点都有对应的匹配,则称该匹配为完备匹配。
浅析算法正确性
该算法是通过给每个顶点一个标号(叫做顶标)来把求最大权匹配的问题转化为求完备匹配的问题的。设顶点 (xi) 的顶标为 (lx[i]),顶点 (y_j) 的顶标为 (ly[j]) ,边权为 (mp[i][j]) 。在算法执行过程中的任一时刻,对于任一条Edge (<i,j>) 有 (lx[i] + ly[j] ge mp[i][j]) 始终成立。
给出一条定理:若由二分图中所有满足(lx[i] + ly[j] = mp[i][j]) 的边 (<i,j>) 构成的子图(称做相等子图)有完备匹配,那么这个完备匹配就是二分图的最大权匹配。
这是显然的,因为对于每个匹配他其他的选择的边权一定满足 (mp[i][j] le lx[i] + ly[j]) 要保证其有最大值,只能选相等子图的
其最佳答案即为 (sum lx[i] +ly[i]) 所以
初始时为了使 (lx[i] + ly[j] ge mp[i][j]) 恒成立,令 (lx[i]) 为所有与顶点 (x_i) 关联的边的最大权,(ly[j] = 0)。如果当前的相等子图没有完备匹配,就按下面的方法修改顶标以使扩大相等子图,直到相等子图具有完备匹配为止。
我们求当前相等子图的完备匹配失败了,是因为对于某个 (x) 顶点,我们找不到一条从它出发的交错路。这时我们获得了一棵交错树,它的叶子结点全部是 (x) 顶点。我们把交错树中 (x) 顶点的顶标全都减小某个值 (d), (y) 顶点的顶标全都增加同一个值 (d) ,那么我们会发现:
1)两端都在交错树中的边 (<i,j>) (lx + ly)(以下简称xy)的值没有变化。也就是说,它原来属于相等子图,仍属于相等子图。
2)两端都不在交错树中的边 (<i, j>) ,(xy) 都没有变化。也就是说,它原来属于(或不属于)相等子图,仍属于(或不属于)相等子图。
3)(x) 端不在交错树中,(y) 端在交错树中的边 (<i,j>) ,它的 (xy) 的值有所增大。它原来不属于相等子图,仍不属于相等子图。
4)(x) 端在交错树中,(y) 端不在交错树中的边 (<i,j>) ,它的 (xy) 的值有所减小。它原来不属于相等子图,可能进入了相等子图,因而使相等子图得到了扩大。
5)到最后,(x) 端每个点至少有一条线连着,(y) 端每个点有一条线连着,说明最后补充完的相等子图一定有完备匹配。所以,最大权子图子图是一步一步变大的。
而显然的为了满足性质有: (d = min(lx[i]+ly[j]-mp[i][j])) 其中 $ x_i$ 在交替树里,(y_i) 不在。
为啥?我们要扩大相等子图。
对于 (d) 的求解可以在求解交替树时完成
为啥我上面的都划掉了,ssinger 提出了一种反例,观察第三个,我们可以知道,还会有的边原来在子图里,但是他的x端减小了!所以就会出子图,但对算法的正确性没有什么影响,因为
代码
bool find(int u)
{
visx[u] = 1;
for(int i = 1; i <= n; ++i)
{
if(visy[i]) continue;
if(lx[u] + ly[i] == mp[u][i])
{
visy[i] = 1;
if(!p[i] || find(p[i]))
{
p[i] = u;
return 1;
}
}
else d = min(d, lx[u] + ly[i] - mp[u][i]);
/* 在这里计算d */
}
return 0;
}
/* 当然了,如果可以匹配成功,那就不用d了 */
main:
for(int k = 1; k <= n; ++k)
{
while(1)
{
memset(visx, 0, sizeof visx);
memset(visy, 0, sizeof visy);
d = 1 << 30;
if(find(k)) break;
for(int i = 1; i <= n; ++i)
{
if(visx[i]) lx[i] -= d;
if(visy[i]) ly[i] += d;
}
}
}
简单口胡下时间复杂度:
循环 (n) 次 每次找增广修改 (n) 次顶标 用 (n) 次算 (d) 所以为 (O(n^3)) 但是怎么会跑到...
之前口胡错了。应该是 (n^4) 因为不是 (n) 次算 (d)
所以似乎之前写的不是多好用了有一个 (d) 是全程量,但是这并不是很好,因为计算的不是很准确。在下文会介绍一个使用 (slac[ ]) 数组来更精准的来算出delta(lamda)。
话说,为啥不用费用流来跑,因为 ek 和 zkw 都挺慢的,KM 和 zkw 复杂度在上界上都为 (O(n^4)),但是 KM 一般会更快
在随机数据下,使用 slac 数组可以在随机数据下平均复杂度约为 (O(n^3))
介绍一下,定义 (slac[v]) 表示对于你 (find()) 函数中把 (d) 改为 (slac[v]) 主要原因是要精准,那么我们的最后的 (d) 是啥,应该是 (!visy[v]) 时里的最小值。就是说你要填几个边才行..就是要找不在交替树的最小 slac 从 x 加边
代码:(hdu3435)
#define _CRT_SECURE_NO_WARNINS
#include <bits/stdc++.h>
using namespace std;
template <typename T>
inline T read()
{
T x = 0;
char ch = getchar();
bool f = 0;
while(ch < '0' || ch > '9')
{
f = (ch == '-');
ch = getchar();
}
while(ch <= '9' && ch >= '0')
{
x = (x << 1) + (x << 3) + (ch - '0');
ch = getchar();
}
return f? -x : x;
}
template <typename T>
void put(T x)
{
if(x < 0)
{
x = -x;
putchar('-');
}
if(x < 10) {
putchar(x + 48);
return;
}
put(x / 10);
putchar(x % 10 + 48);
return ;
}
const int Maxn = 1001, inf = 0x3f3f3f3f;
int cnt, h[Maxn], lx[Maxn], ly[Maxn], slac[Maxn], p[Maxn], t, n, m;
bool visx[Maxn], visy[Maxn];
struct Edge
{
int to, lac, wg;
void insert(int x, int y, int z) { to = y; lac = h[x]; h[x] = cnt++; wg = z; }
}edge[Maxn << 6];
void add_edge(int w, int v, int u)
{
edge[cnt].insert(u, v, w);
edge[cnt].insert(v, u, w);
lx[u] = max(lx[u], w);
lx[v] = max(lx[v], w);
}
bool find(int u)
{
visx[u] = 1;
for(int i = h[u]; i != -1; i = edge[i].lac)
{
int to = edge[i].to;
if(visy[to]) continue;
// 注意这里的位置
if(lx[u] + ly[to] == edge[i].wg)
{
visy[to] = 1;
if(!p[to] || find(p[to]))
{
p[to] = u;
return 1;
}
}
else slac[to] = min(slac[to], lx[u] + ly[to] - edge[i].wg);
}
return 0;
}
int KM()
{
memset(p, 0, sizeof p);
for(int k = 1; k <= n; ++k)
{
memset(slac, 63, sizeof slac);
while(1)
{
memset(visx, 0, sizeof visx);
memset(visy, 0, sizeof visy);
int minz = 0x3f3f3f3f;
if(find(k)) break;
for(int i = 1; i <= n; ++i) if(!visy[i] && minz > slac[i]) minz = slac[i];
if(minz == 0x3f3f3f3f) return -1;
for(int i = 1; i <= n; ++i)
{
if(visx[i]) lx[i] -= minz;
if(visy[i]) ly[i] += minz;
else slac[i] -= minz;
}
}
}
int ans = 0;
for(int i = 1; i <= n; ++i) ans += lx[i] + ly[i];
return -ans;
}
int main()
{
freopen("in.txt", "r", stdin);
int cnt1 = 0;
t = read <int> ();
while(t--)
{
n = read <int> ();
m = read <int> ();
memset(h, -1, sizeof h); cnt = 0;
for(int i = 1; i <= n; ++i) lx[i] = -inf, ly[i] = 0;
for(int i = 1; i <= m; ++i) add_edge(-read <int> (), read <int> () + 1, read <int> () + 1);
int ans = KM();
if(ans == -1) printf("Case %d: NO
", ++cnt1);
else printf("Case %d: %d
", ++cnt1, ans);
}
return 0;
}
关于判无解,我本人提出了一种想法,就是先跑匈牙利(笑),singercoder 提出了一种想法,当 d = inf 时,但是由于他那个神奇的常数问题(所以导致了死循环)....神奇吧....其实这个想法是对的-.-
所以真的是很神奇,KM 的坑实在是太多了。
当然了,由于 (O(n^4)) 的复杂度无法接受,有的人研究出了一种 bfs 的操作
这个码是 uoj#80 的
#define _CRT_SECURE_NO_WARNINS
#include <bits/stdc++.h>
using namespace std;
template <typename T>
inline T read()
{
T x = 0;
char ch = getchar();
bool f = 0;
while(ch < '0' || ch > '9')
{
f = (ch == '-');
ch = getchar();
}
while(ch <= '9' && ch >= '0')
{
x = (x << 1) + (x << 3) + (ch - '0');
ch = getchar();
}
return f? -x : x;
}
template <typename T>
void put(T x)
{
if(x < 0)
{
x = -x;
putchar('-');
}
if(x < 10) {
putchar(x + 48);
return;
}
put(x / 10);
putchar(x % 10 + 48);
return ;
}
const int Maxn = 404;
int mp[Maxn][Maxn], n, bl, br, m, lx[Maxn], ly[Maxn], link[Maxn], p[Maxn], slac[Maxn], pre[Maxn], mat[Maxn];
bool visx[Maxn], visy[Maxn];
// 把增广路上的取反
void aug(int x)
{
if(!x) return;
link[x] = pre[x];
aug(mat[pre[x]]);
mat[pre[x]] = x;
return;
}
// 通过 bfs 来寻找增广路
void bfs(int u)
{
int delta;
memset(slac, 63, sizeof slac);
memset(visx, 0, sizeof visx);
memset(visy, 0, sizeof visy);
queue <int> q;
q.push(u);
while(1)
{
// 通过广搜来寻找
while(!q.empty())
{
int fr =q.front();
q.pop();
visx[fr] = 1;
for(int i = 1; i <= n; ++i)
{
// 这里就跟dfs的差不多,但是需要用pre来记录路径
// 在这里 visy 的位置十分关键
if(visy[i]) continue;
delta = lx[fr] + ly[i] - mp[fr][i];
if(!delta)
{
visy[i] = 1;
pre[i] = fr;
if(!link[i]) { aug(i); return ; }
q.push(link[i]);
}
else if(slac[i] > delta) slac[i] = delta, pre[i] = fr;
}
}
delta = 0x3f3f3f3f;
for(int i = 1; i <= n; ++i) if(!visy[i]) delta = min(delta, slac[i]);
if(delta == 0x3f3f3f3f) return;
for(int i = 1; i <= n; ++i)
{
if(visx[i]) lx[i] -= delta;
// 左边的点顶标减去 delta
if(visy[i]) ly[i] += delta;
// 遍历的点加上 delta
else slac[i] -= delta;
// 记住没有打标记的点就 slac -= delta
}
// 把新增的点添到队列中
for(int i = 1; i <= n; ++i)
{
if(!visy[i] && !slac[i])
{
visy[i] = 1;
if(!link[i]) { aug(i); return ; }
// 如果可以就返回
q.push(link[i]);
}
}
}
}
void KM()
{
for(int i = 1; i <= n; ++i)
{
lx[i] = 0;
for(int j = 1; j <= n; ++j) lx[i] = max(lx[i], mp[i][j]);
}
for(int i = 1; i <= n; ++i)
bfs(i);
long long ans = 0;
for(int i = 1; i <= br; ++i)
{
if(!mp[link[i]][i]) continue;
ans += mp[link[i]][i];
p[link[i]] = i;
}
put <long long> (ans); putchar('
');
for(int i = 1; i <= bl; ++i) put <int> (p[i]), putchar(' ');
return ;
}
int main()
{
#ifdef _DEBUG
freopen("in.txt", "r", stdin);
#endif
bl = read <int> ();
br = read <int> ();
n = max(bl, br);
m = read <int> ();
while(m--) mp[read <int> ()][read <int> ()] = read <int> ();
KM();
return 0;
}
这个时间复杂度是 (O(n^3))
嵬
我没有那么好的天份 ,我只是多下点力气
关于时间复杂度的证明,好像 algorithm 里给了一个主定理之类的东西吧。去学习下,要肝日校了,可能会比较晚写