Венгерский алгоритм решения задачи о назначениях

Венгерский алгоритм (англ. Hungarian algorithm) — алгоритм, решающий задачу о назначениях за полиномиальное время. Оригинальная версия была придумана и разработана Х. Куном в 1955 году и имела асимптотику [math] O(n^4) [/math], но позже Эдмонс и Карп (а также, независимо от них, Томидзава) показали, что можно улучшить ее до [math] O(n^3) [/math].

Задача:

Пусть дан взвешенный полный двудольный граф c целыми весами ребер [math] K_{n, n} [/math], нужно найти в нем полное паросочетание минимального веса. Вес паросочетания определяется как сумма весов его ребер. Далее будем обозначать левую и правую доли графа за [math] X [/math] и [math] Y [/math] соответственно, вес ребра [math] xy [/math] — как [math] c(xy) [/math].

Вспомогательные леммы

Лемма:

Если веса всех ребер графа, инцидентных какой-либо вершине, изменить (увеличить или уменьшить) на одно и то же число, то в новом графе оптимальное паросочетание будет состоять из тех же ребер, что и в старом.

Доказательство:

[math]\triangleright[/math]

Полное паросочетание для каждой вершины содержит ровно одно ребро, инцидентное этой вершине. Указанная операция изменит на одно и то же число вес любого паросочетания. Значит, ребро, которое принадлежало оптимальному паросочетанию в старом графе, в новом графе тоже будет ему принадлежать.

[math]\triangleleft[/math]

Далее будем рассматривать только графы с неотрицательной весовой функцией, так как, согласно этой лемме, задачу о назначениях на остальных графах можно свести к задаче о назначениях на них.

Лемма:

Выделим в множествах [math]X[/math] и [math]Y[/math] подмножества [math]X', Y'[/math]. Пусть [math]d = \min \{c(xy) \mid x \in X \setminus X', y \in Y'\}[/math]. Прибавим [math] d [/math] ко всем весам ребер, инцидентных вершинам из [math]X'[/math]. Затем отнимем [math] d [/math] от всех весов ребер, инцидентных вершинам из [math]Y'[/math] (далее для краткости эта операция обозначается как [math] X' \uparrow\downarrow Y' [/math]). Тогда: Веса всех ребер графа останутся неотрицательными. Веса ребер вида [math]xy[/math], где [math]x \in X', y \in Y'[/math] или [math]x \in X \backslash X', y \in Y \backslash Y'[/math], не изменятся.

Доказательство:

[math]\triangleright[/math]

Рассмотрим матрицу весов графа. Не умаляя общности, можно сказать, что множества [math] X' [/math] и [math] Y' [/math] состоят из первых элементов множеств [math] X [/math] и [math] Y [/math] соответственно (мы упорядочиваем множества по номерам вершин). Тогда вся матрица делится на 4 блока: [math] Y' [/math] [math] Y \backslash Y' [/math] [math] X' [/math] [math] A + d - d [/math] [math] C + d [/math] [math] X \backslash X' [/math] [math] B - d [/math] [math] D [/math] Веса группы [math] A [/math] будут сначала увеличены, а потом уменьшены на [math] d [/math], поэтому они не изменятся, веса группы [math] D [/math] вообще изменяться не будут. Все веса группы [math] B [/math] будут уменьшены на [math] d [/math], но [math] d [/math] — минимум среди этих весов, поэтому они останутся неотрицательными.

[math]\triangleleft[/math]

Лемма:

Если веса всех ребер графа неотрицательны и некоторое полное паросочетание состоит из ребер нулевого веса, то оно является оптимальным.

Доказательство:

[math]\triangleright[/math]

Действительно, паросочетание с какими-то другими весами ребер имеет больший вес и оптимальным не является.

[math]\triangleleft[/math]

Общий метод

Доказанные ранее утверждения позволяют придумать схему алгоритма, решающего задачу о назначениях: нужно найти полное паросочетание из ребер нулевого веса в графе, полученном из исходного преобразованиями, описанными в первых двух леммах.

Алгоритм, решающий задачу, работает с графом, как с матрицей весов.

1. Вычитаем из каждой строки значение ее минимального элемента. Теперь в каждой строке есть хотя бы один нулевой элемент.
2. Вычитаем из каждого столбца значение его минимального элемента. Теперь в каждом столбце есть хотя бы один нулевой элемент.
3. Ищем в текущем графе полное паросочетание из ребер нулевого веса:
4. • Если оно найдено, то желаемый результат достигнут, алгоритм закончен.
   • В противном случае, покроем нули матрицы весов минимальным количеством строк и столбцов (это не что иное, как нахождение минимального вершинного покрытия в двудольном графе). Пусть [math] X_c [/math] и [math] Y_c [/math] — множества вершин минимального вершинного покрытия из левой и правой долей (то есть, строк и столбцов) соответственно, тогда применим преобразование [math] X_c \uparrow\downarrow (Y \setminus Y_c) [/math]. Для этого преобразования [math] d [/math] будет минимумом по всем ребрам между [math] X \setminus X_c [/math] и [math] Y \setminus Y_c [/math], то есть, ребер нулевого веса здесь нет, поэтому, после его выполнения в матрице весов появится новый нуль. После этого перейдем к шагу 1.

Анализ времени работы

Поиск максимального паросочетания или минимального вершинного покрытия в двудольном графе совершается за [math] O(n^3) [/math] операций. При каждом повторении шагов 1-4 в матрице весов появляется новый нуль. Этот нуль соответствует некоторому новому ребру между вершинами из множеств [math] X \setminus X_c [/math] и [math] Y \setminus Y_c [/math]. Всего в графе n^2 ребер, значит, всего будет совершено не более [math] O(n^2) [/math] итераций внешнего цикла. Поэтому, верхняя оценка времени работы данного метода — [math] O(n^5) [/math]. Более точная оценка довольно сложна и зависит от порядка чисел в матрице весов графа.

Алгоритм за [math] O(n^3) [/math]

Общая идея

Будем добавлять в рассмотрение строки матрицы одну за одной, а не рассматривать их все сразу.

Описание алгоритма

• Добавляем в рассмотрение очередную строку матрицы [math] \mathtt{a} [/math].
• Пока нет увеличивающей цепи, начинающейся в этой строке, пересчитываем потенциал.
• Как только появляется увеличивающая цепь, чередуем паросочетание вдоль неё (включая тем самым последнюю строку в паросочетание), и переходим к началу (к рассмотрению следующей строки).

Реализация

[math] \mathtt{a[1 \dots n][1 \dots m]} [/math] — прямоугольная входная матрица, где [math] \mathtt{n \leqslant m} [/math]. Матрица хранится в 1-индексации.

[math] \mathtt{u[0 \dots n], v[0 \dots n]} [/math] — массивы потенциалов.

[math] \mathtt{p[0 \dots m]} [/math] — массив паросочетания. Для каждого стобца [math] \mathtt{i = 0 \dots m} [/math] он хранит номер соответствующей выбранной строки [math] \mathtt{p[i]} [/math] (или 0, если ничего не выбрано). Полагаем, что [math] \mathtt{p[0]} [/math] равно номеру рассматриваемой строки.

[math] \mathtt{minv[1 \dots m]} [/math] — массив, хранящий для каждого столбца j вспомогательные минимумы, необходимые для быстрого пересчета потенциала.

[math] \mathtt{minv[j] = \min_{i \in Z_1}\{a[i][j] - u[i] - v[j]\}} [/math]

[math] \mathtt{way[1 \dots m]} [/math] — массив, содержащий информацию о том, где эти минимумы достигаются, чтобы мы могли впоследствии восстановить увеличивающую цепочку.

function hungarianAlgorithm(a):
  for i = 1 to n
    p[0] = i
    j0 = 0
    заполняем массивы minv — [math] \infty [/math], used — false
    while true
      used[j0] = true
      i0 = p[j0]
      [math] \delta [/math] = [math] \infty [/math]
      for j = 1 to m
        if used[j] == 0
          cur = a[i0][j] - u[i0] - v[j]
          if cur < minv[j]
            minv[j] = cur
            way[j] = j0
          if minv[j] < [math] \delta [/math]
            [math] \delta [/math] = minv[j]
            j1 = j
      for j = 0 to m
        if used[j]
          u[p[j]] += [math] \delta [/math]
          v[j] -= [math] \delta [/math]
        else
          minv[j] -= [math] \delta [/math]
      j0 = j1
      if p[j0] != 0
        break
    while true
      j1 = way[j0]
      p[j0] = p[j1]
      j0 = j1
      if j0 [math] \neq [/math] 0
        break

Время работы

Оценим время работы алгоритма. Во внешнем цикле мы добавляем в рассмотрение строки матрицы одну за другой. Каждая строка обрабатывается за время [math] O(n^2) [/math], поскольку при этом могло происходить лишь [math] O(n) [/math] пересчётов потенциала (каждый — за время [math] O(n) [/math]), для чего за время [math] O(n^2) [/math] поддерживается массив [math] \mathtt{minv[]} [/math]; алгоритм Куна суммарно отработает за время [math] O(n^2) [/math] (поскольку он представлен в форме [math] O(n) [/math] итераций, на каждой из которых посещается новый столбец).

Итоговая асимптотика составляет [math] O(n^3) [/math].

См. также

• Алгоритм Куна для поиска максимального паросочетания
• Связь максимального паросочетания и минимального вершинного покрытия в двудольных графах

Источники информации

• Асанов М., Баранский В., Расин В. — Дискретная математика: Графы, матроиды, алгоритмы — 2010, 368 стр.
• Венгерский алготитм в Википедии
• Визуализатор алгоритма
• Реализация венгерского алгоритма на C++
• Венгерский алгоритм решения задачи о назначениях