指派问题算法

指派问题

指派问题数学模型

指派问题的一般形式如下:
KaTeX parse error: No such environment: align* at position 8: \begin{̲a̲l̲i̲g̲n̲*̲}̲ \min\ &z=\sum_…
特别地,当 m = n m=n m=n时,简化为如下形式:
KaTeX parse error: No such environment: align* at position 8: \begin{̲a̲l̲i̲g̲n̲*̲}̲ \min\ &z=\sum_…

  • 从数学模型上看,指派问题是一个特殊的整数线性规划问题
  • 求最大化的指派问题可以转化为求最小化的指派问题求解

二分图的最大匹配

二分图的最大匹配问题和指派问题的关系

二分图的最大匹配问题等价于如下的指派问题:
KaTeX parse error: No such environment: align* at position 8: \begin{̲a̲l̲i̲g̲n̲*̲}̲ \max\ &z=\sum_…
特别地,当 m = n m=n m=n时,简化为如下形式:
KaTeX parse error: No such environment: align* at position 8: \begin{̲a̲l̲i̲g̲n̲*̲}̲ \max\ &z=\sum_…

其中的 C C C矩阵是0-1矩阵,表示二分图的邻接矩阵,例如
C = [ 0 1 1 0 0 1 0 0 0 ] C = \begin{bmatrix} 0&1&1\\ 0&0&1\\ 0&0&0 \end{bmatrix} C=000100110
不过需要注意的是最后解得的 x i j = 1 x_{ij}=1 xij=1时,我们需要回到邻接矩阵中检查对应的 c i j c_{ij} cij是否为1,为1代表最大匹配包含 i i i j j j的匹配,为0则可以忽略

求解二分图的最大匹配问题

二分图匹配有如下性质:

  • 如果一个匹配是完备的(complete),即不可以再通过transfer(也称增广路径,见下图)的方式扩大匹配了,则它就是最大匹配(具体证明见Kuhn, H. W. . “The Hungarian method for the assignment problem.” Naval Research Logistics 52.1-2(2010):7–21.

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-LzLeZDsQ-1618162893824)(imgs/bigraph.png)]

所以从任意一个匹配出发,如果能通过transfer的方式扩大匹配,则继续扩大,直到不能继续扩大时,便得到一个最大匹配

带权二分图的最大匹配

带权二分图的最大匹配问题和指派问题的关系

二分图的最大匹配问题等价于如下的指派问题:
KaTeX parse error: No such environment: align* at position 8: \begin{̲a̲l̲i̲g̲n̲*̲}̲ \max\ &z=\sum_…
特别地,当 m = n m=n m=n时,简化为如下形式:
KaTeX parse error: No such environment: align* at position 8: \begin{̲a̲l̲i̲g̲n̲*̲}̲ \max\ &z=\sum_…

其中的 C C C矩阵是非负矩阵,表示二分图的带权邻接矩阵,例如
C = [ 0 1 2 0 0 3 0 0 0 ] C = \begin{bmatrix} 0&1&2\\ 0&0&3\\ 0&0&0 \end{bmatrix} C=000100230
不过需要注意的是最后解得的 x i j = 1 x_{ij}=1 xij=1时,我们需要回到邻接矩阵中检查对应的 c i j c_{ij} cij是否大于0,大于0代表最大匹配包含 i i i j j j的匹配,为0则可以忽略

求解带权二分图的最大匹配问题

注:本节如下内容均需要假定 m = n m=n m=n,即左右节点数目相同

可行顶标:给每个节点分配一个权值 ,对于所有边满足 w ( u , v ) ≤ l ( u ) + r ( v ) w(u,v)\le l(u)+r(v) w(u,v)l(u)+r(v), l ( ⋅ ) , r ( ⋅ ) l(\cdot),r(\cdot) l(),r()分别表示对二分图的左边节点和右边节点分配权值的函数

相等子图:在一组可行顶标下原图的生成子图,包含所有点但只包含满足 w ( u , v ) = l ( u ) + r ( v ) w(u,v)=l(u)+r(v) w(u,v)=l(u)+r(v)的边

  • 如果存在邻接矩阵为0的情况,此时看作连接了一条权值为0的边,最后求解出来完整的最大匹配后再去掉权值为0的匹配即可

定理:对于某组可行顶标,如果其相等子图存在完美匹配,那么,该匹配就是原二分图的最大权完美匹配

证明:考虑原二分图任意一组完美匹配 M M M,其边权和为
∑ u , v ∈ M w ( u , v ) ≤ ∑ u , v ∈ M l ( u ) + r ( v ) ≤ ∑ u l ( u ) + ∑ v r ( v ) \sum_{u,v\in M}w(u,v)\le\sum_{u,v\in M}l(u)+r(v)\le\sum_{u}l(u)+\sum_vr(v) u,vMw(u,v)u,vMl(u)+r(v)ul(u)+vr(v)
任意一组可行顶标的相等子图的完美匹配 M ′ M' M的边权和为
∑ u , v ∈ M ′ w ( u , v ) = ∑ u , v ∈ M ′ l ( u ) + r ( v ) = ∑ u l ( u ) + ∑ v r ( v ) \sum_{u,v\in M'}w(u,v)=\sum_{u,v\in M'}l(u)+r(v)=\sum_{u}l(u)+\sum_vr(v) u,vMw(u,v)=u,vMl(u)+r(v)=ul(u)+vr(v)
于是可得
∑ u , v ∈ M w ( u , v ) ≤ ∑ u , v ∈ M ′ w ( u , v ) \sum_{u,v\in M}w(u,v)\le\sum_{u,v\in M'}w(u,v) u,vMw(u,v)u,vMw(u,v)
故该匹配 M ′ M' M就是原二分图的最大权完美匹配

算法(一般称之为匈牙利算法):

  1. 初始化 l ( u ) = max ⁡ j w ( i , j ) , r ( v ) = 0 l(u)=\max_{j}w(i,j),r(v)=0 l(u)=maxjw(i,j),r(v)=0

  2. 在相等子图中,从任意一个匹配出发,如果能通过transfer的方式扩大匹配,则继续扩大,直到不能继续扩大为止,如果此时已经达到完美匹配,该匹配就是原二分图的最大权完美匹配,算法停止,否则转3

  3. S S S T T T表示此时二分图中相等子图的无法通过transfer的方式扩大匹配的匹配(称之为交错树)的所有左边节点和右边节点, S ′ S' S T ′ T' T表示不在交错树的所有左边节点和右边节点,令
    α = min ⁡ { l ( x ) + r ( y ) − w ( x , y ) ∣ x ∈ S , y ∈ T ′ } \alpha = \min\{l(x)+r(y)-w(x,y)|x\in S,y\in T'\} α=min{l(x)+r(y)w(x,y)xS,yT}
    S S S中的节点的顶标均减 α \alpha α, T T T中节点顶标均加 α \alpha α, 此时 S S S中节点和 T ′ T' T中节点在原图中的边可能进入相等子图,从而此时交错树可能可以通过transfer的方式继续扩大匹配,如果不行,继续执行3,最后一定能够到达相等子图的完美匹配,该匹配就是原二分图的最大权完美匹配,算法停止

    注:此处 α \alpha α的选取本质上是使得 S S S中的点连接到 T ′ T' T中使得权值和下降最小的选取方式

求解指派问题的Munkres算法

指派问题的性质

  • m = n m=n m=n时,对 C C C矩阵的每行或每列同时加上或减去 t t t时,最优解不变,最优值会加上或减去 t t t(视前面执行的是加还是减而定)

  • m = n m=n m=n时,假定 C C C非负矩阵,我们对每行同时减去该行最小值,那么每行均会出现0元素,如果有 n n n个零元素它们当中任意两个均不同行不同列(此时称它们为独立零元素),那么这些0所对应的 n n n个位置就是原问题的一组最优解

    于是求解原指派问题转化为寻找 n n n个独立零元素

Munkres算法

Step 0: Create an nxm matrix called the cost matrix in which each element represents the cost of assigning one of n workers to one of m jobs. Rotate the matrix so that there are at least as many columns as rows and let k=min(n,m).

Step 1: For each row of the matrix, find the smallest element and subtract it from every element in its row. Go to Step 2.

Step 2: Find a zero (Z) in the resulting matrix. If there is no starred zero in its row or column, star Z. Repeat for each element in the matrix. Go to Step 3.

Step 3: Cover each column containing a starred zero. If K columns are covered, the starred zeros describe a complete set of unique assignments. In this case, Go to DONE, otherwise, Go to Step 4.

Step 4: Find a noncovered zero and prime it. If there is no starred zero in the row containing this primed zero, Go to Step 5. Otherwise, cover this row and uncover the column containing the starred zero. Continue in this manner until there are no uncovered zeros left. Save the smallest uncovered value and Go to Step 6.

Step 5: Construct a series of alternating primed and starred zeros as follows. Let Z0 represent the uncovered primed zero found in Step 4. Let Z1 denote the starred zero in the column of Z0 (if any). Let Z2 denote the primed zero in the row of Z1 (there will always be one). Continue until the series terminates at a primed zero that has no starred zero in its column. Unstar each starred zero of the series, star each primed zero of the series, erase all primes and uncover every line in the matrix. Return to Step 3.

Step 6: Add the value found in Step 4 to every element of each covered row, and subtract it from every element of each uncovered column. Return to Step 4 without altering any stars, primes, or covered lines.

DONE: Assignment pairs are indicated by the positions of the starred zeros in the cost matrix. If C(i,j) is a starred zero, then the element associated with row i is assigned to the element associated with column j.

Munkres算法具体流程参考资料:

  • https://brc2.com/the-algorithm-workshop/
  • https://users.cs.duke.edu/~brd/Teaching/Bio/asmb/current/Handouts/munkres.html

Munkres算法的几点说明

  1. Munkres算法本质上就是就是不断通过对 C C C矩阵的每行或每列同时加上或减去某个数,产生更多的0元素,然后一直在其中不断标记独立0元素的位置,当独立0元素的个数达到行数或列数的最小值时,这些0所对应的个位置就是原问题的一组最优解
  2. Munkres算法可以适用于 m × n m\times n m×n矩阵,不仅仅适用于方阵(更多细节见An extension of the Munkres algorithm for the assignment problem to rectangular matrices. Communications of the ACM, 142302-806, 1971.
    • 对于 m × n m\times n m×n矩阵的维数较小者,比如行,那么做行的加减一定还是同解的,对于列的话,需要注意的是Munkres算法只对没有被覆盖的列进行了减去最小值的变换,而减去最小值一定出来新的0,也就是说必有一列会加入最优解,而对于会匹配的列做加减一定还是同解的,而其他不会匹配的列加减也不影响结果
  3. Munkres算法在实现时可以修改使得可以适用于不允许某对 ( i , j ) (i,j) (i,j)匹配,可以将 C o s t ( i , j ) Cost(i,j) Cost(i,j)设置为无穷大,然后找到所有有限值的 C o s t ( i , j ) Cost(i,j) Cost(i,j),把 C o s t ( i , j ) Cost(i,j) Cost(i,j)为无穷大的值变为比这些有限值的 C o s t ( i , j ) Cost(i,j) Cost(i,j)的最大值还大的有限值,然后调用Munkres算法求解,对解出来的 ( i , j ) (i,j) (i,j)配对进行检查,如果它们在原始 C o s t ( i , j ) Cost(i,j) Cost(i,j)上是无穷大,则丢弃这对 ( i , j ) (i,j) (i,j)配对即可

值得一提的是Munkres算法和前面说的匈牙利算法本质上是等价的,匈牙利算法是由Kuhn整理发表的,所以Munkres算法有时也被称为Kuhn-Munkres算法,简称K-M算法,而Kuhn整理发表匈牙利算法的思想是来自K¨onig,而K¨onig又参考了一位匈牙利数学家的文章,并且这位匈牙利数学家的结论更具普遍性,所以Kuhn因此以匈牙利算法对其命名

In a footnote, K¨onig refers to a paper of E. Egerv´ary (in Hungarian), which seemed to contain the treatment of a more general case.

Munkres算法的Matlab代码如下:

function [assignment, cost] = assignmentoptimal(distMatrix)
%ASSIGNMENTOPTIMAL    Compute optimal assignment by Munkres algorithm
%		ASSIGNMENTOPTIMAL(DISTMATRIX) computes the optimal assignment (minimum
%		overall costs) for the given rectangular distance or cost matrix, for
%		example the assignment of tracks (in rows) to observations (in
%		columns). The result is a column vector containing the assigned column
%		number in each row (or 0 if no assignment could be done).
%
%		[ASSIGNMENT, COST] = ASSIGNMENTOPTIMAL(DISTMATRIX) returns the
%		assignment vector and the overall cost.
%
%		The distance matrix may contain infinite values (forbidden
%		assignments). Internally, the infinite values are set to a very large
%		finite number, so that the Munkres algorithm itself works on
%		finite-number matrices. Before returning the assignment, all
%		assignments with infinite distance are deleted (i.e. set to zero).
%
%		A description of Munkres algorithm (also called Hungarian algorithm)
%		can easily be found on the web.
%
%		<a href="assignment.html">assignment.html</a>  <a href="http://www.mathworks.com/matlabcentral/fileexchange/6543">File Exchange</a>  <a href="https://www.paypal.com/cgi-bin/webscr?cmd=_s-xclick&hosted_button_id=EVW2A4G2HBVAU">Donate via PayPal</a>
%
%		Markus Buehren
%		Last modified 05.07.2011

% Copyright (c) 2004, Markus Buehren
% All rights reserved.
% 
% Redistribution and use in source and binary forms, with or without 
% modification, are permitted provided that the following conditions are 
% met:
% 
%     * Redistributions of source code must retain the above copyright 
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright 
%       notice, this list of conditions and the following disclaimer in 
%       the documentation and/or other materials provided with the distribution
%       
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
% POSSIBILITY OF SUCH DAMAGE.

% save original distMatrix for cost computation
originalDistMatrix = distMatrix;

% check for negative elements
if any(distMatrix(:) < 0)
	error('All matrix elements have to be non-negative.');
end

% get matrix dimensions
[nOfRows, nOfColumns] = size(distMatrix);

% check for infinite values
finiteIndex   = isfinite(distMatrix);
infiniteIndex = find(~finiteIndex);
if ~isempty(infiniteIndex)
	% set infinite values to large finite value
	maxFiniteValue = max(max(distMatrix(finiteIndex)));
	if maxFiniteValue > 0
		infValue = abs(10 * maxFiniteValue * nOfRows * nOfColumns);
	else
		infValue = 10;
	end
	if isempty(infValue)
		% all elements are infinite
		assignment = zeros(nOfRows, 1);
		cost       = 0;
		return
	end	
	distMatrix(infiniteIndex) = infValue;
end

% memory allocation
coveredColumns = zeros(1,       nOfColumns);
coveredRows    = zeros(nOfRows, 1);
starMatrix     = zeros(nOfRows, nOfColumns);
primeMatrix    = zeros(nOfRows, nOfColumns);

% preliminary steps
if nOfRows <= nOfColumns
	minDim = nOfRows;
	
	% find the smallest element of each row
	minVector = min(distMatrix, [], 2);
	
	% subtract the smallest element of each row from the row
	distMatrix = distMatrix - repmat(minVector, 1, nOfColumns);
	
	% Steps 1 and 2
	for row = 1:nOfRows
		for col = find(distMatrix(row,:)==0)
			if ~coveredColumns(col)%~any(starMatrix(:,col))
				starMatrix(row, col) = 1;
				coveredColumns(col)  = 1;
				break
			end
		end
	end
	
else % nOfRows > nOfColumns
	minDim = nOfColumns;
	
	% find the smallest element of each column
	minVector = min(distMatrix);
	
	% subtract the smallest element of each column from the column
	distMatrix = distMatrix - repmat(minVector, nOfRows, 1);
	
	% Steps 1 and 2
	for col = 1:nOfColumns
		for row = find(distMatrix(:,col)==0)'
			if ~coveredRows(row)
				starMatrix(row, col) = 1;
				coveredColumns(col)  = 1;
				coveredRows(row)     = 1;
				break
			end
		end
	end
	coveredRows(:) = 0; % was used auxiliary above
end

if sum(coveredColumns) == minDim
	% algorithm finished
	assignment = buildassignmentvector__(starMatrix);
else
	% move to step 3
	[assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = step3__(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim); %#ok
end

% compute cost and remove invalid assignments
[assignment, cost] = computeassignmentcost__(assignment, originalDistMatrix, nOfRows);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function assignment = buildassignmentvector__(starMatrix)

[maxValue, assignment] = max(starMatrix, [], 2);
assignment(maxValue == 0) = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [assignment, cost] = computeassignmentcost__(assignment, distMatrix, nOfRows)

rowIndex   = find(assignment);
costVector = distMatrix(rowIndex + nOfRows * (assignment(rowIndex)-1));
finiteIndex = isfinite(costVector);
cost = sum(costVector(finiteIndex));
assignment(rowIndex(~finiteIndex)) = 0;

% Step 2: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = step2__(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim)

% cover every column containing a starred zero
maxValue = max(starMatrix);
coveredColumns(maxValue == 1) = 1;

if sum(coveredColumns) == minDim
	% algorithm finished
	assignment = buildassignmentvector__(starMatrix);
else
	% move to step 3
	[assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = step3__(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim);
end

% Step 3: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = step3__(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim)

zerosFound = 1;
while zerosFound
	
	zerosFound = 0;		
	for col = find(~coveredColumns)
		for row = find(~coveredRows')
			if distMatrix(row,col) == 0
				
				primeMatrix(row, col) = 1;
				starCol = find(starMatrix(row,:));
				if isempty(starCol)
					% move to step 4
					[assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = step4__(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, row, col, minDim);
					return
				else
					coveredRows(row)        = 1;
					coveredColumns(starCol) = 0;
					zerosFound              = 1;
					break % go on in next column
				end
			end
		end
	end
end

% move to step 5
[assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = step5__(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim);

% Step 4: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = step4__(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, row, col, minDim)

newStarMatrix          = starMatrix;
newStarMatrix(row,col) = 1;

starCol = col;
starRow = find(starMatrix(:, starCol));

while ~isempty(starRow)

	% unstar the starred zero
	newStarMatrix(starRow, starCol) = 0;
	
	% find primed zero in row
	primeRow = starRow;
	primeCol = find(primeMatrix(primeRow, :));
	
	% star the primed zero
	newStarMatrix(primeRow, primeCol) = 1;
	
	% find starred zero in column
	starCol = primeCol;
	starRow = find(starMatrix(:, starCol));
	
end
starMatrix = newStarMatrix;

primeMatrix(:) = 0;
coveredRows(:) = 0;

% move to step 2
[assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = step2__(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim);


% Step 5: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = step5__(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim)

% find smallest uncovered element
uncoveredRowsIndex    = find(~coveredRows');
uncoveredColumnsIndex = find(~coveredColumns);
%[s, index1] = min(distMatrix(uncoveredRowsIndex,uncoveredColumnsIndex));
[s, index1] = min(distMatrix(uncoveredRowsIndex,uncoveredColumnsIndex),[],1);
[s, index2] = min(s); %#ok
h = distMatrix(uncoveredRowsIndex(index1(index2)), uncoveredColumnsIndex(index2));

% add h to each covered row
index = find(coveredRows);
distMatrix(index, :) = distMatrix(index, :) + h;

% subtract h from each uncovered column
distMatrix(:, uncoveredColumnsIndex) = distMatrix(:, uncoveredColumnsIndex) - h;

% move to step 3
[assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = step3__(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim);
  • 0
    点赞
  • 5
    收藏
    觉得还不错? 一键收藏
  • 1
    评论

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

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

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值