带权二部图匹配(KM算法)讲解及Matlab实现

在现实生活中经常会遇见匹配问题需要解决,我们手中有两种资源如何让他们之间的组合最优。

举个例子,你是联谊会的执行人,联谊上有一群男生和一群女生,每个人都对自己心目中的对象有明确要求,如身高籍贯等,你已知了每个人的个人信息和他们提出的要求,怎么才能凑出更多对符合对方期望的配对。

或者,你是租车行的老板,你手上有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算法学习笔记

  • 22
    点赞
  • 88
    收藏
    觉得还不错? 一键收藏
  • 16
    评论
很抱歉,刚才我回答的似乎是重复的内容。关于二分图最大匹配算法的应用及Matlab实现,我再给你详细讲解一下。 二分图最大匹配算法是解决二分图匹配问题的一种经典算法。在二分图中,每个顶点可以分为两个集合,一个是左侧顶点集合,一个是右侧顶点集合。匹配是指在这两个集合中每个顶点之间建立一条边,使得每个顶点只与一个顶点相邻。二分图最大匹配问题就是要求在这样的二分图中,找到最大的匹配。 二分图最大匹配算法的应用非常广泛,比如在稳定婚姻问题中,可以将男女分为两个集合,在这两个集合中找到最佳的匹配方案。在工人分配问题中,可以将工人和工作分为两个集合,找到最佳的工作分配方案。在学生选课问题中,可以将学生和课程分为两个集合,找到最佳的选课方案。 在Matlab中,可以使用bipartite_matching函数来实现二分图最大匹配算法。该函数的输入参数是一个二分图的邻接矩阵,输出参数是最大匹配的值和匹配矩阵。你可以通过调用该函数来解决二分图最大匹配问题。 另外,你也可以自己编写二分图最大匹配算法Matlab代码。其中,常用的算法包括匈牙利算法KM算法等。这些算法实现方法可以在网上找到相关的资料和代码。 希望我的回答能够帮助到你。如果你还有其他问题,可以继续问我。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 16
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值