Hungarian Algorithm
Hungarian algorithm is one of the standard operations research algorithms. It’s also known as Munkres algorithm, or Kuhn-Munkres algorithm, to commemorate it’s Hungarian creators. It solves assignment problems in polynomial time, more can be read at Wikipedia.
NOTE: THIS IS NOT COMPLETE
The implementation here does have a fatal flaw (it lacks the marriage coupling algorithm implementation) and will not work in many cases. I just leave it here, in hope that this will make me fix it in the future. Or someone will contribute, what I somewhat doubt.
Do not try to use it in any real applications, as it now has a stupid heuristic plugged instead of the lacking algorithm and then brute force to check the correctness.
You’ve been warned…
Hungarian algorithm is one of the standard operations research algorithms. It’s also known as Munkres algorithm, or Kuhn-Munkres algorithm, to commemorate it’s Hungarian creators. It solves assignment problems in polynomial time, more can be read at Wikipedia.
A quick remark: for convenience reasons I use a simple class that encapsulates a two dimensional array and some helpers. Nothing fancy:
class Hungarian
{
int t[MAXN][MAXN],
n,
minr[MAXN],
minc[MAXN];
public:
void read (void);
void dump (void);
void subtract (void);
void cross_out (void);
void subtract2 (const Hungarian&);
vector<int> > solutions (const Hungarian&);
};
Main function represents what typical pseudocode shows. The dump() and read() simply input/output the array. Other functions shown below.
int main (void)
{
Hungarian o,
c,
a;
vector<int> > sol;
o.read();
c = o;
c.dump();
c.subtract();
c.dump();
do
{
a = c;
a.cross_out();
c.dump();
sol = o.solutions(c);
c.subtract2(a);
c.dump();
}
while(sol.size() == 0);
display(sol);
return 0;
}
First subtraction is the core idea. It takes advantage of the fact that when you subtract a certain value from all elements in a set, their relations are intact. Therefore throughout the algorithm all that you have to do is to subtract such values until you can make the assignment using only 0-cost pairs.
void Hungarian::subtract(void)
{
for(int i = 0; i < n; i++)
minr[i] = minc[i] = INF;
// Row-wise subtract
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
if(t[i][j] < minr[i])
minr[i] = t[i][j];
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
t[i][j] -= minr[i];
// Column-wise subtract
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
if(t[i][j] < minc[j])
minc[j] = t[i][j];
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
t[i][j] -= minc[j];
// Update the minimums
for(int i = 0; i < n; i++)
minr[i] = minc[i] = 0;
}
While crossing out rows and columns it’s important to keep track of double-crossed elements. They’ll be treated differently in the second subtraction.
void Hungarian::cross_out(void)
{
int row_zeroes[MAXN] = {},
col_zeroes[MAXN] = {};
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
if(t[i][j] == 0)
{
row_zeroes[i]++;
col_zeroes[j]++;
}
for(int done = 0; done < n; done++)
{
bool max_is_row;
int max = 0,
max_index;
for(int i = 0; i < n; i++)
{
if(max < row_zeroes[i])
{
max_is_row = true;
max = row_zeroes[i];
max_index = i;
}
if(max < col_zeroes[i])
{
max_is_row = false;
max = col_zeroes[i];
max_index = i;
}
#ifdef DEBUG
printf("erasing %s nr %d with %d zeroes \n", max_is_row ? "row" : "col", max_index, max);
#endif
if(max <1)
break;
if(max_is_row)
{
for(int i = 0; i < n; i++)
if(t[max_index][i] == -1)
t[max_index][i] = -2;
else
t[max_index][i] = -1;
row_zeroes[max_index] = -1;
}
else
{
for(int i = 0; i < n; i++)
if(t[i][max_index] == -1)
t[i][max_index] = -2;
else
t[i][max_index] = -1;
col_zeroes[max_index] = -1;
}
}
}
}
The second subtraction differs from the first one. It’s not done for individual rows/columns. Instead the global minimum of the whole array is subtracted from all the non-crossed elements. And it’s added to double-crossed ones.
void Hungarian::subtract2 (const Hungarian &aux)
{
int minimal = INF;
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
if(t[i][j] < minimal && aux.t[i][j] > 0)
minimal = t[i][j];
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
if(aux.t[i][j] > 0)
t[i][j] -= minimal;
else
if(aux.t[i][j] == -2)
t[i][j] += minimal;
}
To fit into polynomial time constraints it’s easiest to implement some marriage coupling algorithm. Here, I do test whether the result of Hungarian algorithm is the same as the result of brute-force search (in other words - if it’s correct). Therefore it’s non-polynomial, as the brute force is.
// Checks whether the outcome of Hungarian algorithm is valid, by comparing
// with brute force outcome. Returns all the possible solutions.
vector<int> > Hungarian::solutions (const Hungarian &aux)
{
vector<int> initial(n),
current;
vector<int> > result;
int cost,
best = INF;
bool should_be_best;
for(int i = 0; i < n; i++)
initial[i] = i;
current = initial;
do
{
cost = 0;
should_be_best = true;
for(int i = 0; i < n; i++)
{
cost += t[i][current[i]];
if(aux.t[i][current[i]] != 0)
should_be_best = false;
}
if(cost < best)
{
result.clear();
best = cost;
}
if(should_be_best)
{
assert(cost == best);
result.push_back(current);
}
next_permutation(current.begin(), current.end());
}
while(current != initial);
printf("best cost %d\n", best);
return result;
}
Finally there’s how I convert the one dimensional vector of assignments into two dimensional array, that’s normally the Hungarian algorithm output.
void display (vector<int> > sol)
{
for(int s = 0; s < sol.size(); s++)
{
printf("solution %d\n", s);
for(int i = 0; i < sol[s].size(); i++)
{
for(int j = 0; j < sol[s].size(); j++)
if(sol[s][i]==j)
printf("1 ");
else
printf("0 ");
printf("\n");
}
printf("\n");
}
}