E057. Венгерский Algoritmo решения задачи о назначениях

Источник: e-maxx.ru/algo, страница PDF 162.

Постановка задачи о назначениях

Tarea о назначениях ставится весьма естественно. Приведём несколько вариантов постановки (как легко видеть, все они эквивалентны друг другу):

● Есть

рабочих и

заданий. Для каждого рабочего известно, сколько денег он запросит за выполнение того или иного задания. Каждый рабочий может взять себе только одно задание. it is required распределить задания по рабочим так, чтобы минимизировать суммарные расходы.

● Дана матрица

размера

. it is required в каждой её строке выбрать по одному числу так, чтобы в любом столбце также было выбрано ровно по одному числу, и при этом сумма выбранных чисел была бы минимальной.

● Дана матрица

размера

. it is required find такую перестановку

длины

, что величина

— минимальна.

● Дан полный двудольный grafo с

vertexми; каждому ребру приписан некоторый вес. it is required find совершенное паросочетание минимального веса. Отметим, что все приведённые выше постановки "квадратны": в них обе размерности всегда совпадают (и равны ). На практике часто встречаются аналогичные "прямоугольные" постановки, когда

, и надо

выбрать

elementов. Впрочем, как легко заметить, от "прямоугольной" задачи всегда можно перейти к "квадратной", добавив строки/столбцы с нулевыми/бесконечными значениями соответственно. Также заметим, что по аналогии с поиском минимального решения также ставить задачу поиска максимального решения. Впрочем, эти две задачи эквивалентны друг другу: достаточно все веса

умножить на

.

Венгерский Algoritmo

Историческая справка

Algoritmo был разработан и опубликован Гарольдом Куном (Harold Kuhn) в 1955 г. Сам Кун дал Algoritmoу название "венгерский", потому что он был в значительной степени основан на более ранних работах двух венгерских математиков: Денеша Кёнига (Dénes Kőnig) и Эйгена Эгервари (Jenő Egerváry). В 1957 г. Джеймс Манкрес (James Munkres) показал, что этот Algoritmo работает за (строго) полиномиальное время

(т.е. за время порядка полинома от

, не зависящего от величины стоимостей). Поэтому в литературе данный Algoritmo известен не только как "венгерский", но и как "Algoritmo Куна-Манкреса" или "Algoritmo Манкреса". Впрочем, недавно (в 2006 г.) выяснилось, что точно такой же Algoritmo был изобретён за век до Куна немецким математиком Карлом Густавом Якоби (Carl Gustav Jacobi). Дело в том, что его работа "About the research of the order of a system of arbitrary ordinary differential equations", напечатанная посмертно в 1890 г., содержавшая помимо прочих результатов и полиномиальный Algoritmo решения задачи о назначениях, была написана на латыни, а её публикация прошла незамеченной среди математиков. Также стоит отметить, что первоначальный Algoritmo Куна имел асимптотику

, и лишь позже Джек

Эдмондс (Jack Edmonds) и Ричард Карп (Richard Karp) (и независимо от них Томидзава (Tomizawa))

показали, каким образом улучшить его до асимптотики

.

Построение Algoritmoа за

Сразу отметим во избежание неоднозначностей, что мы в основном рассматриваем здесь задачу о назначениях

в матричной постановке (т.е. дана матрица

, и надо выбрать из неё

ячеек, находящихся в разных cadenaх и

столбцах). Индексацию arregloов мы начинаем с единицы, т.е., наEjemplo, матрица

имеет

индексы

.

Также мы будем считать, что все числа в матрице

неотрицательны (если это не так, то всегда можно

перейти к неотрицательной матрице, прибавив ко всем числам некоторое number).

Назовём потенциалом два произвольных arregloа чисел

и

таких, что выполняется Enunciado:

(Как видно, числа

соответствуют cadenaм, а числа

— столбцам матрицы.)

Назовём значением

потенциала сумму его чисел:

С одной стороны, легко заметить, что стоимость искомого решения

не меньше значения любого потенциала:

(Prueba. Искомое Solución задачи представляет из себя

ячеек матрицы, и для каждой из них

выполняется Enunciado

. Поскольку все elementы находятся в разных cadenaх и столбцах,

то, суммируя эти неравенства по всем выбранным

, в левой части неравенства получаем

, а в правой —

,

что и требовалось доказать.)

С другой стороны, оказывается, что всегда существует Solución и потенциал, на которых это неравенство обращается в равенство. Венгерский Algoritmo, описанный ниже, будет конструктивным Pruebaм этого факта. Пока же лишь обратим внимание на то, что если какое-либо Solución имеет стоимость, равную по величине какому-либо потенциалу, то это Solución — оптимально.

Зафиксируем некоторый потенциал. Назовём edge

жёстким, если выполняется: Вспомним об альтернативной постановке задачи о назначениях, с помощью двудольного grafoа. Обозначим через двудольный grafo, составленный только из жёстких рёбер. Фактически, венгерский Algoritmo поддерживает для текущего потенциала максимальное по количеству рёбер паросочетание

grafoа

: и

как только это паросочетание станет содержать

рёбер, рёбра этого паросочетания и будут являться

искомым оптимальным Soluciónм. Перейдём непосредственно к описанию Algoritmoа.

● В начале Algoritmoа потенциал полагается равным нулю

, и паросочетание

полагается пустым.

● Далее, на каждом шаге Algoritmoа мы пытаемся, не меняя потенциала, увеличить мощность текущего паросочетания

на единицу (напоминаем, паросочетание ищется в grafoе жёстких рёбер ). Для этого фактически используется обычный Algoritmo Куна поиска максимального паросочетания в двудольных grafoах. Напомним здесь этот Algoritmo.

Все рёбра паросочетания

ориентируются по направлению от второй доли к первой, все остальные рёбра grafoа ориентируются в противоположную сторону. Напомним (из терминологии поиска паросочетаний), что vertex называется насыщенной, если ей смежно edge из текущего паросочетания. vertex, которой не смежно ни одно edge из текущего паросочетания, называется ненасыщенной. Путь нечётной длины, в котором первое edge не принадлежит паросочетанию, а для всех последующих рёбер происходит чередование (принадлежит/не принадлежит) — называется увеличивающим путём. Из всех ненасыщенных вершин первой доли запускается обход в глубину/в ширину. Если в результате обхода удалось достигнуть ненасыщенной вершины второй доли, то это означает, что мы нашли увеличивающий путь из первой доли во вторую. Если прочередовать рёбра вдоль этого пути (т.е. первое edge включить в паросочетание, второе исключить, третье включить, и т.д.), то тем самым мы увеличим мощность паросочетания на единицу. Если же увеличивающего пути не было, то это означает, что текущее паросочетание

— максимально в grafoе

, поэтому в таком случае переходим к следующему пункту.

● Если на текущем шаге не удалось увеличить мощность текущего паросочетания, то производится некий

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

Обозначим через

множество вершин первой доли, которые были посещены обходом Algoritmoа Куна при

попытке поиска увеличивающей цепи; через

— множество посещённых вершин второй доли. Обозначим через

множество ненасыщенных вершин первой доли, через

— множество ненасыщенных вершин второй доли.

Посчитаем величину

: Эта величина строго положительна.

(Prueba. Предположим, что

. Тогда существует жёсткое edge

, причём

и

.

Из этого следует, что edge

должно было быть ориентированным от второй доле к первой, т.е. это жёсткое

edge

должно Entradaить в паросочетание

. Однако это невозможно, т.к. мы не могли попасть в

насыщенную вершину

, кроме как пройдя по ребру из

в

. Пришли к противоречию, значит,

.)

Теперь пересчитаем потенциал таким образом: для всех вершин

сделаем

, а для

всех вершин

— сделаем

. Получившийся потенциал по-прежнему останется

корректным потенциалом. (Prueba. Для этого надо показать, что по-прежнему для всех

и

выполняется: .

Для случаев, когда

или

— это так, поскольку для них сумма

и

не изменилась. Когда

неравенство только усилилилось. Наконец, для случая

— хотя левая часть неравенства и увеличивается, но всё равно неравенство сохраняется, поскольку величина , согласно её определению — это как раз максимальное увеличение, не приводящее к нарушению неравенства.)

Кроме того, старое паросочетание

из жёстких рёбер можно будет оставить, т.е. все рёбра паросочетания останутся жёсткими.

(Prueba. Чтобы некоторое жёсткое edge

перестало быть жёстким в результате изменения

потенциала, надо, чтобы равенство

превратилось в неравенство

. Однако левая часть могла уменьшиться только в одном случае: когда .

Но раз

, то это означает, что edge

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

,

строго увеличивается. (При этом нельзя утверждать, что увеличивается количество жёстких рёбер.) (Prueba. Во-первых, любая vertex, которая была достижимой, достижимой и останется. В самом деле, если некоторая vertex достижима, то до неё есть некоторый путь из достижимых вершин, начинающийся

в ненасыщенной вершине первой доли; а поскольку для рёбер вида

сумма

не меняется, то весь этот путь сохранится и после изменения потенциала, что и требовалось доказать. Во- вторых, покажем, что в результате пересчёта потенциала появилась хотя бы одна новая достижимая vertex. Но

это почти очевидно, если вернуться к определению

: в самом деле, если минимум был достигнут на ребре

, то

это edge теперь станет жёстким, а, значит, vertex

станет достижимой благодаря этому ребру и вершине

.)

Таким образом, всего может происходить не более

пересчётов потенциала, прежде чем обнаружится

увеличивающая цепочка и мощность паросочетания

будет увеличена. Таким образом, рано или поздно будет найден потенциал, которому соответствует совершенное паросочетание , являющееся ответом на задачу.

Если говорить об асимптотике Algoritmoа, то она составляет

, поскольку всего должно произойти

увеличений паросочетания, перед каждым из которых происходит не более

пересчётов потенциала, каждый

из которых выполняется за время

.

Реализацию за

мы здесь приводить не будем, поскольку она всё равно получится не короче, чем описанная

ниже Implementación за

.

Построение Algoritmoа за

(

)

Научимся теперь реализовывать тот же Algoritmo за асимптотику

(для прямоугольных задач

—

). Ключевая идея: теперь мы будем добавлять в рассмотрение строки матрицы одну за одной, а не рассматривать их все сразу. Таким образом, описанный выше Algoritmo примет вид:

● Добавляем в рассмотрение очередную строку матрицы

.

● Пока нет увеличивающей цепи, начинающейся в этой строке, пересчитываем потенциал.

● Как только появляется увеличивающая цепь, чередуем паросочетание вдоль неё (включая тем самым последнюю строку

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

время

(для прямоугольных задач — за

). Для этого мы вспомним два факта, доказанных нами выше:

● При изменении потенциала вершины, которые были достижимы обходом Куна, достижимыми и останутся.

● Всего могло произойти лишь

пересчётов потенциала, прежде чем будет найдена увеличивающая цепь. Отсюда вытекают ключевые идеи, позволяющие достичь требуемой асимптотики:

● Для проверки наличия увеличивающей цепочки нет необходимости запускать обход Куна заново после каждого

пересчёта потенциала. Вместо этого можно оформить обход Куна в итеративном виде: после каждого пересчёта потенциала мы просматриваем добавившиеся жёсткие рёбра и, если их левые концы были достижимыми, помечаем их правые концы также как достижимые и продолжаем обход из них.

● Развивая эту идею дальше, можно прийти к такому представлению Algoritmoа: это цикл, на каждом шаге которого

сначала пересчитывается потенциал, затем находится столбец, ставший достижимым (а таковой всегда найдётся, поскольку после пересчёта потенциала всегда появляются новые достижимые вершины), и если этот столбец был ненасыщен, то найдена увеличивающая цепь, а если столбец был насыщен — то соответствующая ему в паросочетании cadena также становится достижимой. Теперь Algoritmo принимает вид: цикл добавления столбцов, на каждом из которых сначала пересчитывается потенциал, а затем какой-то новый столбец помечается как достижимый.

● Чтобы быстро пересчитывать потенциал (быстрее, чем наивный вариант за

), надо

поддерживать вспомогательные минимумы по каждому из столбцов

:

Как легко видеть, искомая величина

выражается через них следующим образом:

Таким образом, нахождение

теперь можно произвести за

.

Поддерживать этот arreglo

необходимо при появлении новых посещённых строк. Это, очевидно, можно

сделать за

на одну добавляемую строку (что в сумме даст

). Также обновлять arreglo

надо

при пересчёте потенциала, что также делается за время

на один пересчёт потенциала (поскольку

меняется только для недостигнутых пока столбцов: а именно, уменьшается на ). Таким образом, Algoritmo принимает такой вид: во внешнем цикле мы добавляем в рассмотрение строки матрицы одну

за другой. Каждая cadena обрабатывается за время

, поскольку при этом могло происходить лишь

пересчётов потенциала (каждый — за время

), для чего за время

поддерживается arreglo

; Algoritmo Куна суммарно отработает за время

(поскольку он представлен в форме

итераций,

на каждой из которых посещается новый столбец).

Итоговая Asymptotic complexity составляет

— или, если Tarea прямоугольна,

.

Implementación венгерского Algoritmoа за

(

)

Приведённая Implementación фактически была разработана Андреем Лопатиным несколько лет назад. Её отличает удивительная лаконичность: весь Algoritmo помещается в 30 строк кода. Данная Implementación ищет Solución для прямоугольной Entradaной матрицы

, где

.

Матрица хранится в

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

arregloы

и

хранят потенциал. Изначально он нулевой, что верно для матрицы, состоящей из нуля строк. (Отметим, что для данной реализации не важно, имеются или нет в матрице

отрицательные числа.)

arreglo

содержит паросочетание: для каждого столбца

он хранит номер

соответствующей выбранной строки

(или

, если пока ничего не выбрано). При этом

для удобства

реализации полагается равным номеру текущей рассматриваемой строки.

arreglo

содержит для каждого столбца

вспомогательные минимумы, необходимые для

быстрого пересчёта потенциала:

arreglo

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

для каждого столбца

надо хранить номер строки, а также завести ещё один arreglo: для каждой строки запомнить номер столбца, из которого мы в неё пришли. Однако вместо этого можно заметить, что Algoritmo Куна всегда попадает в строки, проходя по ребру паросочетания из столбцов, поэтому номера строк для восстановления цепочки всегда можно взять

из паросочетания (т.е. из arregloа

). Таким образом,

для каждого столбца

содержит

номер предшествующего столбца (или

, если такого нет). Сам Algoritmo представляет из себя внешний цикл по cadenaм матрицы, внутри которого

происходит добавление в рассмотрение

-ой строки матрицы. Внутренняя часть представляет собой цикл "do-while (p[j0] !

= 0)", который работает, пока не будет найден свободный столбец

. Каждая итерация цикла помечает

посещённым новый столбец с номером

(посчитанным на прошлой итерации; а изначально равным нулю — т.

е. стартуем мы с фиктивного столбца), а также новую строку

— смежную ему в паросочетании (т.е.

;

а изначально при

берётся

-ая cadena). Из-за появления новой посещённой строки

нужно

соответствующим образом пересчитать arreglo

, заодно мы находим минимум в нём — величину

, и

в каком столбце

этот минимум был достигнут (заметим, что при такой реализации

могло оказаться

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

того). После этого производится пересчёт потенциала

, соответствующее изменение arregloа

. По окончании цикла "do-while" мы нашли увеличивающую цепочку, оканчивающуюся в столбце

, "раскрутить"

которую можно, пользуясь arregloом предков

.

Константа

— это "бесконечность", т.е. некоторое number, заведомо большее всех возможных чисел во

Entradaной матрице

.

vector<int> u (n+1), v (m+1), p (m+1), way (m+1);

for (int i=1; i<=n; ++i) {

p[0] = i;

int j0 = 0;

vector<int> minv (m+1, INF);

vector<char> used (m+1, false);

do {

used[j0] = true;

int i0 = p[j0],  delta = INF,  j1;

for (int j=1; j<=m; ++j)

if (!used[j]) {

int cur = a[i0][j]-u[i0]-v[j];

if (cur < minv[j])

minv[j] = cur, way[j] = j0;

if (minv[j] < delta)

delta = minv[j], j1 = j;

}

for (int j=0; j<=m; ++j)

if (used[j])

u[p[j]] += delta, v[j] -= delta;

else

minv[j] -= delta;

j0 = j1;

} while (p[j0] != 0);

do {

int j1 = way[j0];

p[j0] = p[j1];

j0 = j1;

} while (j0); } Восстановление ответа в более привычной форме, т.е. нахождение для каждой строки

номера выбранного

в ней столбца

, делается следующим образом:

vector<int> ans (n+1);

for (int j=1; j<=m; ++j)

ans[p[j]] = j;

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

содержит в себе сумму всех

величин

, т.е. суммарное изменение потенциала. Хотя при каждом изменении потенциала изменяться могли

сразу несколько величин

и

, суммарное изменение величины потенциала в точности равно

,

поскольку пока нет увеличивающей цеп

...

using System;
using System.Collections.Generic;
using System.Linq;

public static class AlgorithmDraft
{
    // Auto-generated C# draft from the original e-maxx C/C++ listing. Review before production use.
    List<int> u (n+1), v (m+1), p (m+1), way (m+1);
    for (int i=1; i<=n; ++i) {
            p[0] = i;
            int j0 = 0;
            List<int> minv (m+1, INF);
            List<char> used (m+1, false);
            do {
                    used[j0] = true;
                    int i0 = p[j0],  delta = INF,  j1;
                    for (int j=1; j<=m; ++j)
                            if (!used[j]) {
                                    int cur = a[i0][j]-u[i0]-v[j];
                                    if (cur < minv[j])
                                            minv[j] = cur,  way[j] = j0;
                                    if (minv[j] < delta)
                                            delta = minv[j],  j1 = j;
                            }
                    for (int j=0; j<=m; ++j)
                            if (used[j])
                                    u[p[j]] += delta,  v[j] -= delta;
                            else
                                    minv[j] -= delta;
                    j0 = j1;
            } while (p[j0] != 0);
            do {
                    int j1 = way[j0];
                    p[j0] = p[j1];
                    j0 = j1;
            } while (j0);
    }
    List<int> ans (n+1);
    for (int j=1; j<=m; ++j)
            ans[p[j]] = j;
    int cost = -v[0];
}

vector<int> u (n+1), v (m+1), p (m+1), way (m+1);
for (int i=1; i<=n; ++i) {
        p[0] = i;
        int j0 = 0;
        vector<int> minv (m+1, INF);
        vector<char> used (m+1, false);
        do {
                used[j0] = true;
                int i0 = p[j0],  delta = INF,  j1;
                for (int j=1; j<=m; ++j)
                        if (!used[j]) {
                                int cur = a[i0][j]-u[i0]-v[j];
                                if (cur < minv[j])
                                        minv[j] = cur,  way[j] = j0;
                                if (minv[j] < delta)
                                        delta = minv[j],  j1 = j;
                        }
                for (int j=0; j<=m; ++j)
                        if (used[j])
                                u[p[j]] += delta,  v[j] -= delta;
                        else
                                minv[j] -= delta;
                j0 = j1;
        } while (p[j0] != 0);
        do {
                int j1 = way[j0];
                p[j0] = p[j1];
                j0 = j1;
        } while (j0);
}
vector<int> ans (n+1);
for (int j=1; j<=m; ++j)
        ans[p[j]] = j;
int cost = -v[0];

import java.util.*;
import java.math.*;

public class AlgorithmDraft {
    // Auto-generated Java draft from the original e-maxx C/C++ listing. Review before production use.
    ArrayList<Integer> u (n+1), v (m+1), p (m+1), way (m+1);
    for (int i=1; i<=n; ++i) {
            p[0] = i;
            int j0 = 0;
            ArrayList<Integer> minv (m+1, INF);
            ArrayList<Character> used (m+1, false);
            do {
                    used[j0] = true;
                    int i0 = p[j0],  delta = INF,  j1;
                    for (int j=1; j<=m; ++j)
                            if (!used[j]) {
                                    int cur = a[i0][j]-u[i0]-v[j];
                                    if (cur < minv[j])
                                            minv[j] = cur,  way[j] = j0;
                                    if (minv[j] < delta)
                                            delta = minv[j],  j1 = j;
                            }
                    for (int j=0; j<=m; ++j)
                            if (used[j])
                                    u[p[j]] += delta,  v[j] -= delta;
                            else
                                    minv[j] -= delta;
                    j0 = j1;
            } while (p[j0] != 0);
            do {
                    int j1 = way[j0];
                    p[j0] = p[j1];
                    j0 = j1;
            } while (j0);
    }
    ArrayList<Integer> ans (n+1);
    for (int j=1; j<=m; ++j)
            ans[p[j]] = j;
    int cost = -v[0];
}