在现实生活中经常会遇见匹配问题需要解决,我们手中有两种资源如何让他们之间的组合最优。
举个例子,你是联谊会的执行人,联谊上有一群男生和一群女生,每个人都对自己心目中的对象有明确要求,如身高籍贯等,你已知了每个人的个人信息和他们提出的要求,怎么才能凑出更多对符合对方期望的配对。
或者,你是租车行的老板,你手上有N台车,各自配置不同,有豪华跑车,有越野型的,有商务SUV等等,每一种类型可能不止一辆。今天来了M位顾客,每位顾客的要求也明确,有的就要租越野车,有的要租SUV。你现在面临的问题就是怎么把顾客和车进行配对,从而使得收益最大化。
上述的两个问题都是二部图匹配的问题,具体说是从一个二部图中找出它的最大匹配。这种问题都可以通过上一篇文章中的匈牙利算法来得到解决。
但是现实生活中有些问题比上述问题更加复杂,比如租车问题中好几台豪华跑车都能满足同一个顾客的要求,但是出于对外观的好恶,顾客对这几台车的出价可能不同。而且每个顾客对每辆车的出价都不同。换句话说,二部图中连接左右节点的边加上了权重。那么,我们现在如何让自己收益最大化。显然匈牙利算法没法解决这个问题,因为它平等的看待每一条边。在此介绍一种新的算法,KM(Kuhn-Munkres Algorithm)算法来解决这种带权二部图中的最优匹配问题。
还是和上一篇文章一样从一个具体的例子来看算法在运行过程中的步骤和思想。
从租车问题说起,现在你有5辆车,
y
1
y_1
y1到
y
5
y_5
y5,并且现在来了5位顾客
x
1
x_1
x1到
x
5
x_5
x5。每位顾客给每辆车的一天租金的报价列在右边的标中,如给2号顾客分配4号车,他一天只愿意出3单位的钱,但如果给他分配5号车,他一天愿意出8单位的钱。我们需要计算怎么分配车辆达到最大收益。
KM算法我也看了几天,最难理解的对我而言就是右节点的label的意义。因此这里先不从标准的KM算法讲起,而是从我简化过的"伪KM算法"作为引入来讲解。
伪KM算法
对于每位顾客,首先我们给他们贴上一个label(记为
L
(
x
i
)
L(x_i)
L(xi)),记录该顾客身上有多少钱,假设每个顾客只会带自己最高出价那么多的钱,因为再带多了也没用。因此
L
(
x
i
)
=
m
a
x
(
w
e
i
g
h
t
(
x
i
,
y
j
)
)
L(x_i)=max(weight(x_i, y_j))
L(xi)=max(weight(xi,yj)),其中j可以取任意有效值。理想情况下在把每个人的荷包都清空,我们就达到了最大收益,但是经常会遇见冲突,如1号顾客和2号顾客都给5号车出价最高,没法同时满足,因此就要进行妥协。实际上KM算法就是从最理想的状态不断进行妥协的过程。
第一步
和匈牙利的结构类似,还是对于左边的节点一个一个进行处理。首先对
x
1
x_1
x1,直接让他匹配他现在的钱包能承受的最贵的车
y
5
y_5
y5,即要求
L
(
x
i
)
=
w
e
i
g
h
t
(
x
i
,
y
j
)
L(x_i)=weight(x_i,y_j)
L(xi)=weight(xi,yj)才能匹配,j可以取任意有效值,此时
y
5
y_5
y5没有被匹配,将
x
1
x_1
x1与
y
5
y_5
y5进行匹配,完成对
x
1
x_1
x1的处理。
第二步
对
x
2
x_2
x2进行匹配,让他匹配他现在的钱包能承受的最贵的车
y
5
y_5
y5,但是现在
y
5
y_5
y5已经被
x
1
x_1
x1匹配了。现在需要让
x
1
x_1
x1或
x
2
x_2
x2其中一人进行妥协,但是
x
1
x_1
x1表示他只盯着他出价最高的
y
5
y_5
y5,同样
x
2
x_2
x2也只盯着他出价最高的
y
5
y_5
y5,两人都没有其它选择,杠上了。如何对他俩进行仲裁?
既然你俩对于
y
5
y_5
y5的归属有争议,那么就进行“竞标”,
x
1
x_1
x1和
x
2
x_2
x2都砸出
d
d
d单位金钱,同时把自己对于
y
5
y_5
y5的出价调低
d
d
d单位(保证自己永远能买得起
y
5
y_5
y5),看谁最先耐不住。同时不能让其它人占便宜,因此需要大家一起把对
y
5
y_5
y5的出价调低
d
d
d单位。
那么
d
d
d如何取值?
首先
d
d
d从0开始慢慢变大,变到多少才能停?就目前的情况来看,
x
1
x_1
x1和
x
2
x_2
x2都有望获得
y
5
y_5
y5,但同样谁也没有把握一定得到
y
5
y_5
y5,因此
x
1
x_1
x1和
x
2
x_2
x2需要给自己留后路,即保证自己能租得起冲突的车除外的所有车。
用公式表达就是
d
=
m
i
n
(
L
(
x
i
)
−
w
e
i
g
h
t
(
x
i
,
y
j
)
)
d=min(L(x_i)-weight(x_i, y_j))
d=min(L(xi)−weight(xi,yj)),
x
i
x_i
xi为起冲突的顾客,
y
j
y_j
yj为不冲突的车辆。
就这种条件来看,对
x
1
x_1
x1而言,最多
d
d
d能取3,否则他就租不起
y
3
y_3
y3。对
x
2
x_2
x2而言,最多
d
d
d能取2,否则他就租不起
y
1
y_1
y1。因此总结下来
d
d
d取2,即从
x
1
x_1
x1和
x
2
x_2
x2的钱包中各押上2单位金钱,同时所有人对于
y
5
y_5
y5的出价调低2单位。
这时候
x
1
x_1
x1还能继续把
d
d
d加大。但
x
2
x_2
x2发现自己已经受不了了,并且现在对他而言
y
1
y_1
y1和
y
5
y_5
y5的出价都一样,他不用再死盯着
y
5
y_5
y5,可以转投
y
1
y_1
y1了。
重新对
x
2
x_2
x2进行匹配,此时
y
1
y_1
y1没有被匹配,将
x
2
x_2
x2与
y
1
y_1
y1进行匹配,因此完成对
x
2
x_2
x2的处理。
第三步
对
x
3
x_3
x3进行匹配,让他匹配他现在的钱包能承受的最贵的车
y
1
y_1
y1,但是现在
y
1
y_1
y1已经被
x
2
x_2
x2匹配了。与第二步刚开始不同,现在
x
2
x_2
x2还有选择,即
y
5
y_5
y5,
x
3
x_3
x3请求
x
2
x_2
x2去选择
y
5
y_5
y5,把
y
1
y_1
y1腾出来。
y
1
y_1
y1再回去跟
x
1
x_1
x1‘“竞标”一次争夺
y
5
y_5
y5,但同样这次
x
1
x_1
x1表示他只盯着他出价最高的
y
5
y_5
y5不放,导致
x
2
x_2
x2让不出
y
1
y_1
y1,三人陷入僵局。那么现在要开始对他们进行仲裁。
x
1
x_1
x1、
x
2
x_2
x2和
x
3
x_3
x3都砸出
d
d
d单位金钱,然后大家一起把受争议的
y
1
y_1
y1和
y
5
y_5
y5的出价调低
d
d
d单位。
基于
x
1
x_1
x1考虑,为保证之后还能选择
y
3
y_3
y3,
d
d
d最大取1;
基于
x
2
x_2
x2考虑,为保证之后还能选择
y
3
y_3
y3,
d
d
d最大取1;
基于
x
3
x_3
x3考虑,为保证之后还能选择
y
2
y_2
y2,
d
d
d最大取2;
综上
d
d
d取1。
调整后变为:
重新对
x
3
x_3
x3进行匹配,此时步骤类似于匈牙利算法,
x
3
x_3
x3要
x
2
x_2
x2挪窝,
x
2
x_2
x2发现自己现在可以选
y
3
y_3
y3了,并且不冲突,于是挪过去让出
y
1
y_1
y1给
x
3
x_3
x3。
第四步
对
x
4
x_4
x4进行匹配,让他匹配他现在的钱包能承受的最贵的车
y
1
y_1
y1,但是现在
y
1
y_1
y1已经被
x
3
x_3
x3匹配了。
x
3
x_3
x3只盯着对他的最高价
y
1
y_1
y1不让步。仲裁。
x
3
x_3
x3和
x
4
x_4
x4都砸出
d
d
d单位金钱,然后大家一起把受争议的
y
1
y_1
y1的出价调低
d
d
d单位。
基于
x
3
x_3
x3考虑,为保证之后还能选择
y
2
y_2
y2,
d
d
d最大取1;
基于
x
4
x_4
x4考虑,为保证之后还能选择
y
2
y_2
y2,
d
d
d最大取2;
综上
d
d
d取1。
调整后变为:
再次对
x
4
x_4
x4进行匹配,
x
4
x_4
x4叫
x
3
x_3
x3去选
y
2
y_2
y2,问题解决。
第五步
对
x
5
x_5
x5进行匹配,让他匹配他现在的钱包能承受的最贵的车
y
1
y_1
y1,但是现在
y
1
y_1
y1已经被
x
4
x_4
x4匹配了。
x
4
x_4
x4只盯着对他的最高价
y
1
y_1
y1不让步。仲裁。
x
4
x_4
x4和
x
5
x_5
x5都砸出
d
d
d单位金钱,然后大家一起把受争议的
y
1
y_1
y1的出价调低
d
d
d单位。
基于
x
4
x_4
x4考虑,为保证之后还能选择
y
2
y_2
y2,
d
d
d最大取1;
基于
x
5
x_5
x5考虑,为保证之后还能选择
y
3
y_3
y3,
d
d
d最大取1;
综上
d
d
d取1。
调整后变为:
再次对
x
5
x_5
x5进行匹配,类似和匈牙利算法一样,
x
5
x_5
x5首先叫
x
4
x_4
x4去选
y
2
y_2
y2,
x
4
x_4
x4叫
x
3
x_3
x3挪窝腾出
y
2
y_2
y2遭拒,信息返回到
x
5
x_5
x5。
x
5
x_5
x5再去选
y
3
y_3
y3,叫
x
2
x_2
x2挪窝去选
y
5
y_5
y5,
x
2
x_2
x2叫
x
1
x_1
x1挪窝腾出
y
5
y_5
y5,但
x
5
x_5
x5占了
y
3
y_3
y3又没法让。所有人都陷入僵局,争执
y
1
y_1
y1、
y
2
y_2
y2、
y
3
y_3
y3和
y
5
y_5
y5的所有权。
所有人都砸出
d
d
d单位金钱,然后大家一起把受争议的
y
1
y_1
y1、
y
2
y_2
y2、
y
3
y_3
y3和
y
5
y_5
y5的出价调低
d
d
d单位。
基于
x
1
x_1
x1考虑,为保证之后还能选择
y
4
y_4
y4,
d
d
d最大取2;
基于
x
2
x_2
x2考虑,为保证之后还能选择
y
4
y_4
y4,
d
d
d最大取2;
基于
x
3
x_3
x3考虑,为保证之后还能选择
y
4
y_4
y4,
d
d
d最大取1;
基于
x
4
x_4
x4考虑,为保证之后还能选择
y
4
y_4
y4,
d
d
d最大取1;
基于
x
5
x_5
x5考虑,为保证之后还能选择
y
4
y_4
y4,
d
d
d最大取1;
综上
d
d
d取1。
调整后变为:
重新对
x
5
x_5
x5进行匹配,此时步骤同匈牙利算法,得到:
最终结果
所有点都处理完成,算法结束,最终结果:
至此基本上理清了“伪KM算法”的思想,主要是每个人都只能选择现在对各自权重最高的点,以及如果出现冲突后
d
d
d的计算方法,其它的基本上与匈牙利算法基本一致。
KM算法
那么再推广到真正的KM算法上去。之前为了便于理解,直接对右边的权重表进行操作,十分繁复,并且可以看到,每次操作对于同一列都是同时减去同一个值,因此,其实我们可以把对权重表格的操作转换为右侧
y
y
y节点的属性,这就是右节点的label。在KM算法中,只有满足
L
(
x
i
)
+
L
(
y
j
)
=
w
e
i
g
h
t
(
x
i
,
y
j
)
L(x_i)+L(y_j)=weight(x_i,y_j)
L(xi)+L(yj)=weight(xi,yj)时
i
i
i顾客才能选
j
j
j号车。
并且
d
d
d的计算公式变为
d
=
m
i
n
(
L
(
x
i
)
+
L
(
y
j
)
−
w
e
i
g
h
t
(
x
i
,
y
j
)
)
d=min(L(x_i)+L(y_j)-weight(x_i,y_j))
d=min(L(xi)+L(yj)−weight(xi,yj))
左节点的label初始化方法还是一样,右节点的label都初始化为0。每次对左节点的label减去d,对右节点的label加上d。
因此具体算法步骤不再赘述,可以自己去想一想。这里放出每一步过后的结果:
最后附上可运行的KM算法代码。
以及说明一个很重要的问题:KM算法要求左右节点一样多,在我们要解决的问题中左右节点个数不同怎么办,8台车遇见9位顾客。
这时可以通过添加虚拟节点补足差数,把虚拟节点的所有连接权重设成一个特别小的值使它不会对其它正常点的匹配产生影响并在最后去除它就行。
% 带权二部图匹配(KM算法)
clc
clear
close all
global N
global adj_matrix
global label_left
global label_right
global match_right
global visit_left
global visit_right
% 左右各有N个点
% KM算法要求左右两边的节点数相等,可以通过添加虚拟节点的方法实现
N = 5;
adj_matrix = [3 4 6 4 9;
6 4 5 3 8;
7 5 3 4 2;
6 3 2 2 5;
8 4 5 4 7];
% N = 9;
% adj_matrix = round(rand(N)*100);
% adj_matrix(adj_matrix<70) = 0;
% 初始化顶标
label_left = max(adj_matrix, [], 2);
label_right = zeros(N, 1);
% 初始化匹配结果
match_right = ones(N, 1) * nan;
% 初始化辅助变量
visit_left = ones(N, 1) * false;
visit_right = ones(N, 1) * false;
res = KM();
% KM主函数
function res = KM()
global N
global adj_matrix
global label_left
global label_right
global match_right
global visit_left
global visit_right
graph_num = 1;
display_graph(graph_num, '原始二部图');
% 对左边的点依次进行处理
for i = 1: N
while 1
% 重置辅助变量
visit_left = ones(N, 1) * false;
visit_right = ones(N, 1) * false;
% 能找到可行匹配
if find_path(i)
break;
end
% 不能找到可行匹配,修改顶标
% (1)将所有在增广路中的X方点的label全部减去一个常数d
% (2)将所有在增广路中的Y方点的label全部加上一个常数d
d = Inf;
for j = 1: N
if visit_left(j)
for k = 1: N
if ~visit_right(k)
% 左边的点中已经访问过的点,即已经匹配过的点可能需要重新匹配以得到更大的总权值,
% 所以修改顶标,往子图中添加一条边,重新寻找增广路看能不能增广
% 取与左边的点相邻的未匹配边中跟当前存在子图中的以该点为端点的边相差最小的两条边
% 这样才能保持总权值最大
d = min(d, label_left(j) + label_right(k) - adj_matrix(j, k));
end
end
end
end
for k = 1: N
if visit_left(k)
label_left(k) = label_left(k) - d;
end
if visit_right(k)
label_right(k) = label_right(k) + d;
end
end
end
graph_num = graph_num + 1;
display_graph(graph_num, ['第' num2str(i) '步']);
end
graph_num = graph_num + 1;
display_graph(graph_num, '最终结果');
res = 0;
for j = 1: N
if match_right(j) >=0 && match_right(j) < N
res = res + adj_matrix(match_right(j), j);
end
end
end
% 寻找增广路,深度优先
function result = find_path(i)
global N
global adj_matrix
global label_left
global label_right
global match_right
global visit_left
global visit_right
visit_left(i) = true;
for j = 1: length(adj_matrix(i, :))
match_weight = adj_matrix(i, j);
if visit_right(j)
% 已被匹配(解决递归中的冲突)
continue;
end
gap = label_left(i) + label_right(j) - match_weight;
% 当 gap == 0 时 x_i 和 y_j 之间存在一条边,且该边是当前 x_i 可以匹配的权值最大的边
if gap == 0
% 找到可行匹配
visit_right(j) = true;
% j未被匹配,或虽然j已被匹配,但是j的已匹配对象有其他可选备胎
% 此处同匈牙利算法
if isnan(match_right(j)) || find_path(match_right(j))
match_right(j) = i;
result = true;
return;
end
end
end
result = false;
return;
end
function display_graph(graph_num, title_name)
global N
global adj_matrix
global label_left
global label_right
global match_right
global visit_left
global visit_right
figure(graph_num);
set(gcf, 'Position', [100, 100, 1000, 500]);
subplot(1,2,1);
cla
set( gca, 'XTick', [], 'YTick', [] );
set( gca, 'TickLength', [0 0]);
box on
hold on
xlim([-1, 2]);
ylim([-N, 1]);
temp = N - 1;
% 绘制匹配边
for j = 1: length(match_right)
if ~isnan(match_right(j))
plot([0, 1], [(match_right(j) - 1) * -temp/(N-1), (j - 1) * -temp/(N-1)], 'r', 'LineWidth', 2);
end
end
% 绘制点
scatter(zeros(1, N), 0:-temp/(N-1):-temp, 20, [217/255 83/255 25/255], 'filled');
scatter(ones(1, N), 0:-temp/(N-1):-temp, 20, [0 114/255 189/255], 'filled');
for j = 1: N
text(-0.2, (j - 1) * -temp/(N-1), ['x_' num2str(j)]);
end
for j = 1: N
text(1.1, (j - 1) * -temp/(N-1), ['y_' num2str(j)]);
end
for j = 1: N
text(-0.5, (j - 1) * -temp/(N-1), num2str(label_left(j)), 'Color', [217/255 83/255 25/255], 'FontSize', 15);
end
for j = 1: N
text(1.4, (j - 1) * -temp/(N-1), num2str(label_right(j)), 'Color', [0 114/255 189/255], 'FontSize', 15);
end
title(title_name);
% 边权重
subplot(1,2,2);
cla
set( gca, 'XTick', [], 'YTick', [] );
set( gca, 'TickLength', [0 0]);
box on
hold on
xlim([-1, N]);
ylim([-N, 1]);
for j = 0 : N-1
plot([j, j], [-N, 1], 'k') ;
end
for j = -N+1 : 0
plot([-1, N], [j, j], 'k') ;
end
for j = 1: N
display_text(j, 0, ['x_' num2str(j)], 0);
end
for j = 1: N
display_text(0, j, ['y_' num2str(j)], 0);
end
for i = 1: N
for j = 1: N
if match_right(j) == i
display_text(i, j, num2str(adj_matrix(i,j)), 1);
else
display_text(i, j, num2str(adj_matrix(i,j)), 0);
end
end
end
% saveas(gcf, [num2str(graph_num) '.png']);
end
% 在x行y列显示t文本
function display_text(x, y, t, bold)
if bold
text(y - 0.5, -x + 0.5, t, 'FontSize', 15, 'FontWeight', 'bold', 'Color', 'r');
else
text(y - 0.5, -x + 0.5, t, 'FontSize', 15, 'FontWeight', 'normal', 'Color', 'k');
end
end
参考资料:
知乎用户洪九(李戈)的专栏文章:策略算法工程师之路-图优化算法(一)(二分图&最小费用最大流);
夜空骑士的CSDN博文:带你入门多目标跟踪(三)匈牙利算法&KM算法;
丹追兵的KM算法学习笔记;