version 1.4.0.0 (3.27 KB) by
Yi Cao

An extremely fast implementation of the Hungarian algorithm on a native Matlab code.

This is an extremely fast implementation of the famous Hungarian algorithm (aslo known as Munkres' algorithm). It can solve a 1000 x 1000 problem in about 20 seconds in a Core Duo (T2500 @ 2.00GHz) XP laptop with Matlab 2008a, which is about 2.5 times faster than the mex code "assignmentoptimal" in FEX ID 6543, about 6 times faster than the author's first version in FEX ID 20328, and at least 30 times faster than other Matlab implementations in the FEX.

The code can also handle rectangular prolems and problems with forbiden allocations.

The new version (V2.3)is able to conduct a partial assignment if a full assignment is not feasible.

For more details of the Hungarian algorithm, visit http://csclab.murraystate.edu/bob.pilgrim/445/munkres.html

Yi Cao (2020). Hungarian Algorithm for Linear Assignment Problems (V2.3) (https://www.mathworks.com/matlabcentral/fileexchange/20652-hungarian-algorithm-for-linear-assignment-problems-v2-3), MATLAB Central File Exchange. Retrieved .

Created with
R2011a

Compatible with any release

**Inspired by:**
assignprob.zip, Functions for the rectangular assignment problem, Munkres Assignment Algorithm

**Inspired:**
Hungarian Algorithm for Linear Sum Assignment Problem, Minimum Cost Constrained Input-Output and Control Configuration Co-Design Problem, Eigenshuffle, LAPJV - Jonker-Volgenant Algorithm for Linear Assignment Problem V3.0, Hungarian based particle linking, simpletracker, Smooth Point-set Registration using Neighboring Constraints, TACTICS Toolbox

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

SergeWould be handy to have visual 2D example, eg:

n = 200; %number of points in set A

i = randperm(n,50); %random subset for set B

A = rand(n,2); %set A (can do ND instead of 2D)

B = A(i,:) + randn(numel(i),2)/50; %set B + noise

C = pdist2(A,B,'euclidean'); %cost matrix

[~,a,b] = find(munkres(C)); %associate A(a,:)<=>B(b,:)

%plot

clf, hold on, plot(A(:,1),A(:,2),'ob',B(:,1),B(:,2),'or'); %plot A(blue) B(red)

quiver(A(i,1),A(i,2),B(:,1)-A(i,1),B(:,2)-A(i,2),0,'g') %plot true associations (magenta)

quiver(A(a,1),A(a,2),B(b,1)-A(a,1),B(b,2)-A(a,2),0,'m') %plot Munkres associations over the top (green)

legend({'A' 'B' 'Truth' 'Munkres'}) %label

Raúl Ibáñez Couoh2Lindsay2I've noticed this code outputs junk costs and non-optimal assignments under certain circumstances, for example, this cost matrix yields the wrong assignment and a cost that is unrelated to the assignment it outputs (Cost=[9 5 8 9; 9 6 8 3; 8 1 7 1; 7 3 5 2]). Beware if you choose to use this.

Tianmin WangDanielJiancao HouProbably worth to mention that the algorithm wouldn't give correct assignment when costMat = [0 1 Inf; Inf Inf 0; Inf Inf 1].

fudongJohn MalikI get the the following error running MATLAB 2018b.

Error using mex

MEX cannot find library 'ut', specified with the -l option.

MEX searched for a file with one of the following names:

libut.lib

ut.lib

Verify the library name is correct. If the library is not

on the existing path, specify the path with the -L option.

chuckMuhammad UmerKOOKMIN UNIVallafi Omrani have matrix with multiple solution how can i get all the solution from your code please

echolove it .

David FrancoIt worked very well.

Guillaume BergerValentin MocanuValWarning, this is highly buggy when working with zero-weight edges. Example:

x = [1 1; 2 4; 2 5; 2 6; 7 7];

costMat = inf(max(x(:,1)), max(x(:,2)));

for itEdge = 1:size(x, 1)

costMat(x(itEdge, 1), x(itEdge, 2)) = 0;

end

[assignment,cost] = munkres(costMat);

matchedEdges = sum(assignment>0)

Here, we should see three matched edges instead of two. When changing the edge costs to 1 instead of 0, it works

PS: Your old version works fine.

ValWarning, this is highly buggy when working with zero-weight edges. Example:

x = [1 1; 2 4; 2 5; 2 6; 7 7];

costMat = inf(max(x(:,1)), max(x(:,2)));

for itEdge = 1:size(x, 1)

costMat(x(itEdge, 1), x(itEdge, 2)) = 0;

end

[assignment,cost] = munkres(costMat);

matchedEdges = sum(assignment>0)

Here, we should see three matched edges instead of two. When changing the edge costs to 1 instead of 0, it works.

Rui WuUttaran BhattacharyaRupert ThomasDavidHi people, who know the maximal number this file can deal with?

Roy SungPrakhar SinhaBiaobin Jiangarnaud ensbergbe careful it gives completely false results for 40*40 matrix !

MatteoHi people,

I a question about the hungarian algorithm. this algorithm is optimal algorithm for the assignment problem, and the time complexity is O(n^3), right?

But , if the input is the multidimensional matrix, it's possible to use the hungarian algorithm? how does it change the algorithm and the time complexity ?

thanks

KuiNote that only after I posted did I see the discussion about JIT and profiler results... I got my numbers by looking at the profiler and not with tic/toc or external timer. Perhaps the profiler is misleading?

KuiI have found a code optimization for the outerplus function using bsxfun. The following code is equivalent but about an order of magnitude faster:

function [minval,rIdx,cIdx]=outerplus(M,x,y)

xPlusYc = bsxfun(@plus, x, y);

M = M - xPlusYc;

minval = min(M(:));

[rIdx, cIdx] = find(M == minval);

For my 60x60 test cases, this reduces the runtime of the overall algorithm by about 30%.

Alexander FarleyMUHAMMAD UZAIRkindly help me by writing short and generalized code for assignment method of scheduling...

Shengjie GuoArun Kumar ChithanarI was wrong about the above comment. It completes, but takes a very long time. For instance, it took 287s for a 100 x 983 matrix. This matrix was from my specific problem. All the weights were real and there were no zero weights.

On the other hand, on the same system it completes a 400 x 400 random in 4 seconds.

Could the specific values in the input matrix cause such a drastic difference in performance?

I checked

Wasit Limprasert>> munkres([1 0])

ans =

1 0

Could you fix this problem?

SyedI am thankful to you for pointing this out as I was comparing your method with another algorithm and it was giving output columnwise.

regards

Yi CaoWell, I can see what you try to do is to increase the cost of selected assignment then to find next best assignment. However, you made a wrong change. The assignment results in dicated row 1 assigned with colume 3, but you miss understood as column 1 assigned with row 3. Wish this helps.

SyedThanks for a nice implementation

I used the code for following matrix

I have the following matrix

M=

67 98 65 95 79 82

50 40 91 66 57 35

61 74 85 112 39 79

41 63 72 97 39 56

56 55 83 91 59 50

66 34 98 70 69 54

Result=munkres(M);

gives me:

Result= 3 4 5 1 6 2

I change the Matrix to

Mp=

67 98 65 205 79 82

50 40 91 66 57 145

171 74 85 112 39 79

41 173 72 97 39 56

56 55 193 91 59 50

66 34 98 70 179 54

Result=munkres(Mp);

gives me:

Result= 3 4 5 1 6 2

cost is also same in both cases.

could u please tell me the reason.

Ricky WongHi, Yi. Thanks for your effort.

However, the function as below is missing

taht prevents me from using your code on Matlab...

bsxfun(@minus, dMat, minR)

Ricky

DQ LIUGreat work, Thanks a lot~~

BTW, I compared your implementation with Niclas Borlin's implementation (http://www.mathworks.com/matlabcentral/fileexchange/94-assignprob-zip). The results are not the same. Yours seems better. Given his is a classic implementation, do you know why? or did you did some improvements over the classic ones?

Turkay YILDIZYi CaoAnoop,

Thanks for pointing out the JV algorithm. You may wish to know that I have implemented a Matlab version of the JV algorithm in

http://www.mathworks.com/matlabcentral/fileexchange/26836

Certainly, it will not be as fast as the mex code, but it has no limitation on the cost matrix to be integer. If you test it, please let me know how it works with your application.

Yi

Anoop Korattikara BalanThanks... You might be right about the memory part. Although I was cleaning the workspace, I was using a MEX code which might have a memory leak (dont know if this is possible).

I ended up using the Jonker-Volgenant shortest augmenting path algorithm from 'http://boguslawobara.net/index.php?page=software.php'. This solves the same problem and appeared way faster.

Yi CaoHi Anoop Balan,

Thanks for commenting on the code. The algorithm is polynomial. This means the computation time will growth about n^p with n the size of the problem and p some constant depending on CPU speed, memory and software implementation. On my PC, Intel Core2 Quad CPU (Q9300) 2.5 GHz with 4 GB RAM, with XP and Matlab 2009b, I got the following results:

n = 1000, cpu time = 18 sec

n = 2000, cpu time = 150 sec

n = 4000, cpu time = 1216 sec

Another machin

n = 1000, cpu time = 128 sec

n = 2000, cpu time = 1024 sec

So roughly, p = 3 (2^3 = 8). I do not know why your machin has so large p (>5). When you solve a large size problem, make sure you start with a clean workspace to avoid any unnecessary overhead for swaping memory with hard disk.

Yi

Anoop Korattikara BalanThanks a lot for for sharing this code!

This works very well for sizes around 1000x1000 (around 35 secs). However, for a 1500x1500 problem it takes around 260 secs and for a 2000x2000 problem it takes around 1670 seconds.

My problem sets are of size 3000x3000 - 4000x4000. Would you know of any approximation algorithms that work well with such sizes?

JamesV. PoorOliver WoodfordExcellent

John D'Erricosplendid!

Yi CaoOh, yes. It was my mistake. A=-PROFIT gives the maximum of sum(PROFIT), but A=1./PROFIT results in the maximum of 1/sum(1./PROFIT), which is different from sum(PROFIT). Sorry for this.

ek deHi Yi,

Using COST = 1./PROFIT doesn't seem to give the same results as with COST = -PROFIT. Try out the code below. In some cases the solutions result with different profits.

%%

clc

clear

for i = 1:100

m = ceil(rand*20)+1;

n = ceil(rand*20)+1;

a = rand(m,n)+eps;

[assign1 cost1] = munkres(-a);

[assign2 cost2] = munkres(1./a);

if ~all(assign1 == assign2)

disp('Different assignments');

if ~all(sort(assign1) == sort(assign2))

disp('Assignment vectors do not agree on permutations');

disp([assign1;assign2]);

end

assign1=assign1(assign1~=0);m1 = length(assign1);

assign2=assign2(assign2~=0);m2 = length(assign2);

if m1 ~= m2

disp('Assignment Vecs not compatible');

continue;

end

cost1 = sum(sum((a.*accumarray([(1:m1)' assign1'],ones(m1,1),[m n]))));

cost2 = sum(sum((a.*accumarray([(1:m1)' assign2'],ones(m1,1),[m n]))));

disp(sprintf('Cost difference = %f',abs(cost1 - cost2)/cost1));

end

end

Yi CaoYes, the code works with negative cost. You can either use negative cost or use reciprocal, i.e. COST = 1./PROFIT if you do not have zero profit elements.

ek deHi Yi,

Will this work if I want to solve the maximum weight matching? That is, if we have a profit matrix rather than a cost, and we want to maximize the profit rather than minimize the cost. I assume that negating the cost matrix should work but was wondering if you could confirm this.

BTW, this is really fast!

Immanuel WeberHey Yi

that was a fast bugfix for your fast algorithm.

switched back to it.

Thank you!

Yi CaoThanks Immanuel. The bug has been fixed.

Immanuel WeberHi Yi, I used your algorithm and it did good work until I encountered a possible bug.

As long as the input matrix ist big enough it works fine, but when the matrix consists of only 2 entries the algorithm creates a wrong match. For example the matrix [1 0], the matching should be 1->2, a return of [2], but the algorithm returns 1->1 ([1]).

Caused by lack of time I switched for the moment to another algorithm, but it would be nice if you could post a bugfix.

In addition, I use it on R2006b, I can't test it on newer versions, so maybe this is related to the old one.

search searchcool, really fast. Thank you.

I test it using Matlab 2007a, will 2008a faster?

Yi CaoHi, Laszlo:

Thanks for comments.

The outerplus function uses JIT acceleration. If you use an old version Matlab, which does not support JIT, then you have to convert it to a vectorized version. Note, if you use profile view, this function will always be slower than vectorized. One way to do vectorization is as you described. Alternatively, you can use bsxfun, which should be faster than repmat.

Laszlo SragnerHi, good stuff, but it is much slower on my machine. It seems outerplus() is the culprit.

Is it possible that this implementation has some 64bit issues?

This helped me a bit:

function [minval,rIdx,cIdx]=outerplus3(M,x,y)

M=M-repmat(x,1,numel(y))-repmat(y,numel(x),1);

minval=min(M(:));

[rIdx,cIdx]=find(M==minval);