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

置顶 john_xia

于 2021-05-25 14:07:53 发布

阅读量5.8k

点赞数 22

分类专栏： matlab 文章标签： matlab 算法匹配

本文链接：https://blog.csdn.net/john_xia/article/details/117247980

版权

matlab 专栏收录该内容

15 篇文章 19 订阅

订阅专栏

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

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

或者，你是租车行的老板，你手上有N台车，各自配置不同，有豪华跑车，有越野型的，有商务SUV等等，每一种类型可能不止一辆。今天来了M位顾客，每位顾客的要求也明确，有的就要租越野车，有的要租SUV。你现在面临的问题就是怎么把顾客和车进行配对，从而使得收益最大化。
上述的两个问题都是二部图匹配的问题，具体说是从一个二部图中找出它的最大匹配。这种问题都可以通过上一篇文章中的匈牙利算法来得到解决。

但是现实生活中有些问题比上述问题更加复杂，比如租车问题中好几台豪华跑车都能满足同一个顾客的要求，但是出于对外观的好恶，顾客对这几台车的出价可能不同。而且每个顾客对每辆车的出价都不同。换句话说，二部图中连接左右节点的边加上了权重。那么，我们现在如何让自己收益最大化。显然匈牙利算法没法解决这个问题，因为它平等的看待每一条边。在此介绍一种新的算法，KM(Kuhn-Munkres Algorithm)算法来解决这种带权二部图中的最优匹配问题。

还是和上一篇文章一样从一个具体的例子来看算法在运行过程中的步骤和思想。
从租车问题说起，现在你有5辆车， $y_1$ 到 $y_5$ ，并且现在来了5位顾客 $x_1$ 到 $x_5$ 。每位顾客给每辆车的一天租金的报价列在右边的标中，如给2号顾客分配4号车，他一天只愿意出3单位的钱，但如果给他分配5号车，他一天愿意出8单位的钱。我们需要计算怎么分配车辆达到最大收益。
在这里插入图片描述
KM算法我也看了几天，最难理解的对我而言就是右节点的label的意义。因此这里先不从标准的KM算法讲起，而是从我简化过的"伪KM算法"作为引入来讲解。

伪KM算法

对于每位顾客，首先我们给他们贴上一个label(记为 $L(x_i)$ )，记录该顾客身上有多少钱，假设每个顾客只会带自己最高出价那么多的钱，因为再带多了也没用。因此 $L(x_i)=max(weight(x_i, y_j))$ ，其中j可以取任意有效值。理想情况下在把每个人的荷包都清空，我们就达到了最大收益，但是经常会遇见冲突，如1号顾客和2号顾客都给5号车出价最高，没法同时满足，因此就要进行妥协。实际上KM算法就是从最理想的状态不断进行妥协的过程。
在这里插入图片描述

第一步

和匈牙利的结构类似，还是对于左边的节点一个一个进行处理。首先对 $x_1$ ，直接让他匹配他现在的钱包能承受的最贵的车 $y_5$ ，即要求 $L(x_i)=weight(x_i,y_j)$ 才能匹配，j可以取任意有效值，此时 $y_5$ 没有被匹配，将 $x_1$ 与 $y_5$ 进行匹配，完成对 $x_1$ 的处理。
在这里插入图片描述

第二步

对 $x_2$ 进行匹配，让他匹配他现在的钱包能承受的最贵的车 $y_5$ ，但是现在 $y_5$ 已经被 $x_1$ 匹配了。现在需要让 $x_1$ 或 $x_2$ 其中一人进行妥协，但是 $x_1$ 表示他只盯着他出价最高的 $y_5$ ，同样 $x_2$ 也只盯着他出价最高的 $y_5$ ，两人都没有其它选择，杠上了。如何对他俩进行仲裁？
既然你俩对于 $y_5$ 的归属有争议，那么就进行“竞标”， $x_1$ 和 $x_2$ 都砸出 $d$ 单位金钱，同时把自己对于 $y_5$ 的出价调低 $d$ 单位（保证自己永远能买得起 $y_5$ ），看谁最先耐不住。同时不能让其它人占便宜，因此需要大家一起把对 $y_5$ 的出价调低 $d$ 单位。

那么 $d$ 如何取值？
首先 $d$ 从0开始慢慢变大，变到多少才能停？就目前的情况来看， $x_1$ 和 $x_2$ 都有望获得 $y_5$ ，但同样谁也没有把握一定得到 $y_5$ ，因此 $x_1$ 和 $x_2$ 需要给自己留后路，即保证自己能租得起冲突的车除外的所有车。
用公式表达就是 $d=min(L(x_i)-weight(x_i, y_j))$ ， $x_i$ 为起冲突的顾客， $y_j$ 为不冲突的车辆。
就这种条件来看，对 $x_1$ 而言，最多 $d$ 能取3，否则他就租不起 $y_3$ 。对 $x_2$ 而言，最多 $d$ 能取2，否则他就租不起 $y_1$ 。因此总结下来 $d$ 取2，即从 $x_1$ 和 $x_2$ 的钱包中各押上2单位金钱，同时所有人对于 $y_5$ 的出价调低2单位。
在这里插入图片描述
这时候 $x_1$ 还能继续把 $d$ 加大。但 $x_2$ 发现自己已经受不了了，并且现在对他而言 $y_1$ 和 $y_5$ 的出价都一样，他不用再死盯着 $y_5$ ，可以转投 $y_1$ 了。
重新对 $x_2$ 进行匹配，此时 $y_1$ 没有被匹配，将 $x_2$ 与 $y_1$ 进行匹配，因此完成对 $x_2$ 的处理。
在这里插入图片描述

第三步

对 $x_3$ 进行匹配，让他匹配他现在的钱包能承受的最贵的车 $y_1$ ，但是现在 $y_1$ 已经被 $x_2$ 匹配了。与第二步刚开始不同，现在 $x_2$ 还有选择，即 $y_5$ ， $x_3$ 请求 $x_2$ 去选择 $y_5$ ，把 $y_1$ 腾出来。 $y_1$ 再回去跟 $x_1$ ‘“竞标”一次争夺 $y_5$ ，但同样这次 $x_1$ 表示他只盯着他出价最高的 $y_5$ 不放，导致 $x_2$ 让不出 $y_1$ ，三人陷入僵局。那么现在要开始对他们进行仲裁。
$x_1$ 、 $x_2$ 和 $x_3$ 都砸出 $d$ 单位金钱，然后大家一起把受争议的 $y_1$ 和 $y_5$ 的出价调低 $d$ 单位。
基于 $x_1$ 考虑，为保证之后还能选择 $y_3$ ， $d$ 最大取1；
基于 $x_2$ 考虑，为保证之后还能选择 $y_3$ ， $d$ 最大取1；
基于 $x_3$ 考虑，为保证之后还能选择 $y_2$ ， $d$ 最大取2；
综上 $d$ 取1。
调整后变为：
在这里插入图片描述
重新对 $x_3$ 进行匹配，此时步骤类似于匈牙利算法， $x_3$ 要 $x_2$ 挪窝， $x_2$ 发现自己现在可以选 $y_3$ 了，并且不冲突，于是挪过去让出 $y_1$ 给 $x_3$ 。

第四步

对 $x_4$ 进行匹配，让他匹配他现在的钱包能承受的最贵的车 $y_1$ ，但是现在 $y_1$ 已经被 $x_3$ 匹配了。 $x_3$ 只盯着对他的最高价 $y_1$ 不让步。仲裁。
$x_3$ 和 $x_4$ 都砸出 $d$ 单位金钱，然后大家一起把受争议的 $y_1$ 的出价调低 $d$ 单位。
基于 $x_3$ 考虑，为保证之后还能选择 $y_2$ ， $d$ 最大取1；
基于 $x_4$ 考虑，为保证之后还能选择 $y_2$ ， $d$ 最大取2；
综上 $d$ 取1。
调整后变为：
在这里插入图片描述
再次对 $x_4$ 进行匹配， $x_4$ 叫 $x_3$ 去选 $y_2$ ，问题解决。

第五步

对 $x_5$ 进行匹配，让他匹配他现在的钱包能承受的最贵的车 $y_1$ ，但是现在 $y_1$ 已经被 $x_4$ 匹配了。 $x_4$ 只盯着对他的最高价 $y_1$ 不让步。仲裁。
$x_4$ 和 $x_5$ 都砸出 $d$ 单位金钱，然后大家一起把受争议的 $y_1$ 的出价调低 $d$ 单位。
基于 $x_4$ 考虑，为保证之后还能选择 $y_2$ ， $d$ 最大取1；
基于 $x_5$ 考虑，为保证之后还能选择 $y_3$ ， $d$ 最大取1；
综上 $d$ 取1。
调整后变为：
在这里插入图片描述
再次对 $x_5$ 进行匹配，类似和匈牙利算法一样， $x_5$ 首先叫 $x_4$ 去选 $y_2$ ， $x_4$ 叫 $x_3$ 挪窝腾出 $y_2$ 遭拒，信息返回到 $x_5$ 。 $x_5$ 再去选 $y_3$ ，叫 $x_2$ 挪窝去选 $y_5$ ， $x_2$ 叫 $x_1$ 挪窝腾出 $y_5$ ，但 $x_5$ 占了 $y_3$ 又没法让。所有人都陷入僵局，争执 $y_1$ 、 $y_2$ 、 $y_3$ 和 $y_5$ 的所有权。
所有人都砸出 $d$ 单位金钱，然后大家一起把受争议的 $y_1$ 、 $y_2$ 、 $y_3$ 和 $y_5$ 的出价调低 $d$ 单位。
基于 $x_1$ 考虑，为保证之后还能选择 $y_4$ ， $d$ 最大取2；
基于 $x_2$ 考虑，为保证之后还能选择 $y_4$ ， $d$ 最大取2；
基于 $x_3$ 考虑，为保证之后还能选择 $y_4$ ， $d$ 最大取1；
基于 $x_4$ 考虑，为保证之后还能选择 $y_4$ ， $d$ 最大取1；
基于 $x_5$ 考虑，为保证之后还能选择 $y_4$ ， $d$ 最大取1；
综上 $d$ 取1。
调整后变为：
在这里插入图片描述
重新对 $x_5$ 进行匹配，此时步骤同匈牙利算法，得到：

最终结果

所有点都处理完成，算法结束，最终结果：
在这里插入图片描述
至此基本上理清了“伪KM算法”的思想，主要是每个人都只能选择现在对各自权重最高的点，以及如果出现冲突后 $d$ 的计算方法，其它的基本上与匈牙利算法基本一致。

KM算法

那么再推广到真正的KM算法上去。之前为了便于理解，直接对右边的权重表进行操作，十分繁复，并且可以看到，每次操作对于同一列都是同时减去同一个值，因此，其实我们可以把对权重表格的操作转换为右侧 $y$ 节点的属性，这就是右节点的label。在KM算法中，只有满足 $L(x_i)+L(y_j)=weight(x_i,y_j)$ 时 $i$ 顾客才能选 $j$ 号车。
并且 $d$ 的计算公式变为
$d=min(L(x_i)+L(y_j)-weight(x_i,y_j))$
左节点的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算法学习笔记；

john_xia

关注

22
点赞
踩
88

收藏

觉得还不错? 一键收藏
16
评论
带权二部图匹配（KM算法）讲解及Matlab实现

在现实生活中经常会遇见匹配问题需要解决，我们手中有两种资源如何让他们之间的组合最优。举个例子，你是联谊会的执行人，联谊上有一群男生和一群女生，每个人都对自己心目中的对象有明确要求，如身高籍贯等，你已知了每个人的个人信息和他们提出的要求，怎么才能凑出更多对符合对方期望的配对。或者，你是租车行的老板，你手上有N台车，各自配置不同，有豪华跑车，有越野型的，有商务SUV等等，每一种类型可能不止一辆。今天来了M位顾客，每位顾客的要求也明确，有的就要租越野车，有的要租SUV。你现在面临的问题就是怎么把顾客和车进行配
复制链接

扫一扫