Hungarian Algorithm for Linear Assignment Problems (V2.1)

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

Description

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 30 seconds in a Core Duo (T2500 @ 2.00GHz) XP laptop with Matlab 2008a, which is about 1.5 times faster than the mex code "assignmentoptimal" in FEX ID 6543, about 4 times faster than the author's first version in FEX ID 20328, and at least 20 times faster than other Matlab implementations in the FEX.

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

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

Acknowledgements

The author wishes to acknowledge the following in the creation of this submission:
assignprob.zip
This submission has inspired the following:
Eigenshuffle

• MATLAB release: MATLAB 7.6 (R2008a)

[Laszlo Sragner]

Hi, 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);

[Yi Cao]

Hi, 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.

[search search]

cool, really fast. Thank you.
I test it using Matlab 2007a, will 2008a faster?

[Immanuel Weber]

Hi 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.

[Yi Cao]

Thanks Immanuel. The bug has been fixed.

[Immanuel Weber]

Hey Yi
that was a fast bugfix for your fast algorithm.
switched back to it.
Thank you!

[ek de]

Hi 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!

[Yi Cao]

Yes, 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 de]

Hi 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 Cao]

Oh, 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.

[John D'Errico]

splendid!

[Oliver Woodford]

Excellent

[V. Poor]

 

[James]