Угорський алгоритм для розв'язання задачі про призначення

Постановка задачі про призначення

Існує кілька стандартних формулювань задачі про призначення (усі вони по суті еквівалентні). Ось деякі з них:

• Є $n$ робіт і $n$ працівників. Кожен працівник вказує суму грошей, яку очікує за конкретну роботу. Кожного працівника можна призначити лише на одну роботу. Мета — розподілити роботи між працівниками так, щоб мінімізувати сумарну вартість.

• Дано матрицю $A$ розміру $n \times n$. Потрібно вибрати по одному числу з кожного рядка так, щоб з кожного стовпця було вибрано рівно одне число, а сума вибраних чисел була мінімальною.

• Дано матрицю $A$ розміру $n \times n$. Потрібно знайти перестановку $p$ довжини $n$ таку, що величина $\sum A[i]\left[p[i]\right]$ мінімальна.

• Розглянемо повний двочастковий граф з $n$ вершинами в кожній частині, де кожному ребру призначено вагу. Мета — знайти досконале парування з мінімальною сумарною вагою.

Важливо зауважити, що всі наведені вище сценарії є "квадратними" задачами, тобто обидві розмірності завжди дорівнюють $n$. На практиці часто трапляються подібні "прямокутні" формулювання, де $n$ не дорівнює $m$, і потрібно вибрати $\min(n,m)$ елементів. Однак можна помітити, що "прямокутну" задачу завжди можна звести до "квадратної", додавши рядки або стовпці зі значеннями нуль чи нескінченність відповідно.

Також зауважимо, що за аналогією з пошуком мінімального розв'язку можна поставити й задачу пошуку максимального розв'язку. Однак ці дві задачі еквівалентні одна одній: достатньо помножити всі ваги на $-1$.

Коли підходить цей алгоритм?

• Це задача про призначення: дано квадратну матрицю вартостей, і потрібно вибрати по одному елементу з кожного рядка та стовпця з мінімальною сумою (досконале парування мінімальної ваги у дводольному графі)?
• Ребра мають ваги (вартості)? (якщо ваг немає і потрібне лише максимальне за потужністю парування → алгоритм Куна)
• Потрібен компактний прямий алгоритм за $O(n^3)$? (ту саму задачу можна звести до потоку мінімальної вартості, але це зазвичай повільніше й громіздкіше)

Угорський алгоритм

Історична довідка

Алгоритм був розроблений і опублікований Гарольдом Куном (Harold Kuhn) у 1955 році. Сам Кун дав йому назву "угорський", оскільки той ґрунтувався на попередніх роботах угорських математиків Денеша Кеніга (Dénes Kőnig) та Єне Егерварі (Jenő Egerváry).
У 1957 році Джеймс Манкрес (James Munkres) показав, що цей алгоритм працює за (строго) поліноміальний час, незалежно від вартостей.
Тому в літературі цей алгоритм відомий не лише як "угорський", а й як "алгоритм Куна–Манкреса" чи "алгоритм Манкреса".
Однак нещодавно, у 2006 році, було виявлено, що той самий алгоритм винайшов за століття до Куна німецький математик Карл Густав Якобі (Carl Gustav Jacobi). Його праця Про дослідження порядку системи довільних звичайних диференціальних рівнянь, опублікована посмертно у 1890 році, серед інших результатів містила поліноміальний алгоритм розв'язання задачі про призначення. На жаль, оскільки публікація була латиною, вона залишилася непоміченою серед математиків.

Варто також зазначити, що оригінальний алгоритм Куна мав асимптотичну складність $\mathcal{O}(n^4)$, і лише пізніше Джек Едмондс (Jack Edmonds) і Річард Карп (Richard Karp) (а незалежно від них — Томідзава) показали, як покращити її до асимптотичної складності $\mathcal{O}(n^3)$.

Алгоритм за $\mathcal{O}(n^4)$

Щоб уникнути неоднозначності, відразу зауважимо, що ми переважно маємо справу із задачею про призначення в матричному формулюванні (тобто дано матрицю $A$, і потрібно вибрати з неї $n$ комірок, що лежать у різних рядках і стовпцях). Масиви ми індексуємо з $1$, тобто, наприклад, матриця $A$ має індекси $A[1 \dots n][1 \dots n]$.

Ми також вважатимемо, що всі числа в матриці A невід'ємні (якщо це не так, матрицю завжди можна зробити невід'ємною, додавши до всіх чисел деяку константу).

Назвемо потенціалом два довільні масиви чисел $u[1 \ldots n]$ та $v[1 \ldots n]$ такі, що виконується наступна умова:

$$u[i]+v[j]\leq A[i][j],\quad i=1\dots n,\ j=1\dots n$$

(Як бачимо, $u[i]$ відповідає $i$-му рядку, а $v[j]$ — $j$-му стовпцю матриці).

Назвемо значенням $f$ потенціалу суму його елементів:

$$f=\sum_{i=1}^{n} u[i] + \sum_{j=1}^{n} v[j].$$

З одного боку, легко бачити, що вартість шуканого розв'язку $sol$ не менша за значення будь-якого потенціалу.

інформація

Лема. $sol\geq f.$

Доведення

Шуканий розв'язок задачі складається з $n$ комірок матриці $A$, тому для кожної з них $u[i]+v[j]\leq A[i][j]$. Оскільки всі елементи в $sol$ лежать у різних рядках і стовпцях, підсумовуючи ці нерівності по всіх вибраних $A[i][j]$, отримуємо $f$ у лівій частині нерівності та $sol$ — у правій.

З іншого боку, виявляється, що завжди існують розв'язок і потенціал, які перетворюють цю нерівність на рівність. Описаний нижче угорський алгоритм буде конструктивним доведенням цього факту. Поки що звернімо лише увагу на те, що якщо вартість якогось розв'язку дорівнює якомусь потенціалу, то цей розв'язок є оптимальним.

Зафіксуємо деякий потенціал. Назвемо ребро $(i,j)$ жорстким, якщо $u[i]+v[j]=A[i][j].$

Пригадаємо альтернативне формулювання задачі про призначення з використанням двочасткового графа. Позначимо через $H$ двочастковий граф, складений лише з жорстких ребер. Для поточного потенціалу угорський алгоритм підтримуватиме парування $M$ з максимальною кількістю ребер у графі $H$. Щойно $M$ міститиме $n$ ребер, розв'язком задачі буде саме $M$ (адже це буде розв'язок, вартість якого збігається зі значенням потенціалу).

Перейдімо безпосередньо до опису алгоритму.

Крок 1. На початку потенціал вважається нульовим ($u[i]=v[i]=0$ для всіх $i$), а парування $M$ вважається порожнім.

Крок 2. Далі на кожному кроці алгоритму ми намагаємося, не змінюючи потенціалу, збільшити потужність поточного парування $M$ на одиницю (нагадаємо, що парування шукається в графі жорстких ребер $H$). Для цього використовується звичайний алгоритм Куна для пошуку максимального парування у двочасткових графах. Нагадаємо цей алгоритм тут. Усі ребра парування $M$ орієнтуються в напрямку від правої частини до лівої, а всі інші ребра графа $H$ — у протилежному напрямку.

Пригадаємо (з термінології пошуку парувань), що вершина називається насиченою, якщо до неї прилягає ребро поточного парування. Вершина, до якої не прилягає жодне ребро поточного парування, називається ненасиченою. Шлях непарної довжини, у якому перше ребро не належить паруванню, а для всіх наступних ребер належність до парування чергується (належить/не належить), називається доповнювальним шляхом. З усіх ненасичених вершин лівої частини запускається пошук у глибину або пошук у ширину. Якщо внаслідок пошуку вдалося дістатися до ненасиченої вершини правої частини, то ми знайшли доповнювальний шлях з лівої частини до правої. Якщо ми включимо до парування непарні ребра шляху й вилучимо парні (тобто включимо до парування перше ребро, виключимо друге, включимо третє і т. д.), то збільшимо потужність парування на одиницю.

Якщо доповнювального шляху не було, то поточне парування $M$ є максимальним у графі $H$.

Крок 3. Якщо на поточному кроці не вдається збільшити потужність поточного парування, то виконується перерахунок потенціалу так, щоб на наступних кроках було більше можливостей для збільшення парування.

Позначимо через $Z_1$ множину вершин лівої частини, які були відвідані під час останнього обходу алгоритму Куна, а через $Z_2$ — множину відвіданих вершин правої частини.

Обчислимо величину $\Delta$:

$$\Delta = \min_{i\in Z_1,\ j\notin Z_2} A[i][j]-u[i]-v[j].$$

інформація

Лема. $\Delta > 0.$

Доведення

Припустимо, що $\Delta=0$. Тоді існує жорстке ребро $(i,j)$ з $i\in Z_1$ та $j\notin Z_2$. Звідси випливає, що ребро $(i,j)$ має бути орієнтоване від правої частини до лівої, тобто $(i,j)$ має входити до парування $M$. Однак це неможливо, бо ми не могли дістатися до насиченої вершини $i$ інакше, ніж пройшовши ребром від j до i. Отже, $\Delta > 0$.

Тепер перерахуємо потенціал таким чином:

• для всіх вершин $i\in Z_1$ робимо $u[i] \gets u[i]+\Delta$,

• для всіх вершин $j\in Z_2$ робимо $v[j] \gets v[j]-\Delta$.

інформація

Лема. Отриманий потенціал залишається коректним потенціалом.

Доведення

Покажемо, що після перерахунку $u[i]+v[j]\leq A[i][j]$ для всіх $i,j$. Для всіх елементів $A$ з $i\in Z_1$ та $j\in Z_2$ сума $u[i]+v[j]$ не змінюється, тому нерівність залишається істинною. Для всіх елементів з $i\notin Z_1$ та $j\in Z_2$ сума $u[i]+v[j]$ зменшується на $\Delta$, тому нерівність так само залишається істинною. Для інших елементів, у яких $i\in Z_1$ та $j\notin Z_2$, сума зростає, але нерівність усе одно зберігається, оскільки величина $\Delta$ за означенням є максимальним приростом, який не порушує нерівність.

інформація

Лема. Старе парування $M$ жорстких ребер залишається коректним, тобто всі ребра парування лишаться жорсткими.

Доведення

Щоб якесь жорстке ребро $(i,j)$ перестало бути жорстким унаслідок зміни потенціалу, потрібно, щоб рівність $u[i] + v[j] = A[i][j]$ перетворилася на нерівність $u[i] + v[j] < A[i][j]$. Однак це може статися лише тоді, коли $i \notin Z_1$ та $j \in Z_2$. Але $i \notin Z_1$ означає, що ребро $(i,j)$ не могло бути ребром парування.

інформація

Лема. Після кожного перерахунку потенціалу кількість вершин, досяжних обходом, тобто $|Z_1|+|Z_2|$, строго зростає.

Доведення

По-перше, зауважимо, що будь-яка вершина, що була досяжною до перерахунку, лишається досяжною. Справді, якщо деяка вершина досяжна, то існує деякий шлях до неї з досяжних вершин, що починається з ненасиченої вершини лівої частини; оскільки для ребер виду $(i,j),\ i\in Z_1,\ j\in Z_2$ сума $u[i]+v[j]$ не змінюється, весь цей шлях збережеться після зміни потенціалу. По-друге, покажемо, що після перерахунку стане досяжною щонайменше одна нова вершина. Це випливає з означення $\Delta$: ребро $(i,j)$, до якого відноситься $\Delta$, стане жорстким, тож вершина $j$ стане досяжною з вершини $i$.

Через останню лему до знаходження доповнювального шляху та збільшення потужності парування $M$ може статися не більше $n$ перерахунків потенціалу. Таким чином, рано чи пізно буде знайдено потенціал, що відповідає досконалому паруванню $M^*$, і $M^*$ буде відповіддю до задачі. Якщо говорити про складність алгоритму, то вона становить $\mathcal{O}(n^4)$: загалом має бути не більше $n$ збільшень парування, перед кожним з яких — не більше $n$ перерахунків потенціалу, кожен з яких виконується за час $\mathcal{O}(n^2)$.

Ми не наводитимемо тут реалізацію алгоритму за $\mathcal{O}(n^4)$, оскільки вона виявиться не коротшою за реалізацію алгоритму за $\mathcal{O}(n^3)$, описану нижче.

Алгоритм за $\mathcal{O}(n^3)$

Тепер навчімося реалізовувати той самий алгоритм за $\mathcal{O}(n^3)$ (для прямокутних задач $n \times m$ — за $\mathcal{O}(n^2m)$).

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

1. Розглядаємо наступний рядок матриці $A$.

2. Поки немає доповнювального шляху, що починається в цьому рядку, перераховуємо потенціал.

3. Щойно знайдено доповнювальний шлях, розповсюджуємо парування вздовж нього (тим самим включаючи останнє ребро до парування) і починаємо з кроку 1 заново (щоб розглянути наступний рядок).

Щоб досягти потрібної складності, необхідно реалізувати кроки 2–3, які виконуються для кожного рядка матриці, за час $\mathcal{O}(n^2)$ (для прямокутних задач — за $\mathcal{O}(nm)$).

Для цього пригадаємо два доведені вище факти:

• Зі зміною потенціалу вершини, що були досяжні обходом Куна, лишаться досяжними.

• Загалом до знаходження доповнювального шляху могло статися лише $\mathcal{O}(n)$ перерахунків потенціалу.

Звідси випливають такі ключові ідеї, що дозволяють досягти потрібної складності:

• Щоб перевірити наявність доповнювального шляху, немає потреби запускати обхід Куна заново після кожного перерахунку потенціалу. Натомість можна виконувати обхід Куна в ітеративній формі: після кожного перерахунку потенціалу дивимося на додані жорсткі ребра і, якщо їхні ліві кінці були досяжними, позначаємо їхні праві кінці теж досяжними та продовжуємо обхід з них.

• Розвиваючи цю ідею далі, можемо подати алгоритм так: на кожному кроці циклу перераховується потенціал. Після цього визначається стовпець, що став досяжним (а такий завжди існуватиме, оскільки після кожного перерахунку потенціалу з'являються нові досяжні вершини). Якщо стовпець ненасичений, то знайдено доповнювальний ланцюг. І навпаки, якщо стовпець насичений, то досяжним стає й рядок парування.

• Щоб швидко перераховувати потенціал (швидше за наївну версію за $\mathcal{O}(n^2)$), потрібно підтримувати допоміжні мінімуми для кожного зі стовпців:

  
  $minv[j]=\min_{i\in Z_1} A[i][j]-u[i]-v[j].$
  

  Легко бачити, що шукана величина $\Delta$ виражається через них так:

  
  $\Delta=\min_{j\notin Z_2} minv[j].$
  

  Таким чином, знаходження $\Delta$ тепер можна виконати за $\mathcal{O}(n)$.

  Необхідно оновлювати масив $minv$, коли з'являються нові відвідані рядки. Це можна зробити за $\mathcal{O}(n)$ для доданого рядка (що сумарно по всіх рядках дає $\mathcal{O}(n^2)$). Також необхідно оновлювати масив $minv$ при перерахунку потенціалу, що теж робиться за час $\mathcal{O}(n)$ ($minv$ змінюється лише для стовпців, які ще не були досягнуті: а саме, зменшується на $\Delta$).

Таким чином, алгоритм набуває наступного вигляду: у зовнішньому циклі ми розглядаємо рядки матриці по одному. Кожен рядок обробляється за час $\mathcal{O}(n^2)$, оскільки могло статися лише $\mathcal{O}(n)$ перерахунків потенціалу (кожен за час $\mathcal{O}(n)$), а масив $minv$ підтримується за час $\mathcal{O}(n^2)$; алгоритм Куна працюватиме за час $\mathcal{O}(n^2)$ (оскільки він поданий у формі $\mathcal{O}(n)$ ітерацій, кожна з яких відвідує новий стовпець).

Отримана складність становить $\mathcal{O}(n^3)$ або, якщо задача прямокутна, $\mathcal{O}(n^2m)$.

Реалізація угорського алгоритму

Наведену нижче реалізацію розробив Андрій Лопатін кілька років тому. Вона вирізняється дивовижною лаконічністю: весь алгоритм складається з 30 рядків коду.

Реалізація знаходить розв'язок для прямокутної матриці $A[1\dots n][1\dots m]$, де $n\leq m$. Матриця 1-індексована для зручності та стислості коду: ця реалізація вводить фіктивний нульовий рядок і нульовий стовпець, що дозволяє нам записувати багато циклів у загальному вигляді, без додаткових перевірок.

Масиви $u[0 \ldots n]$ та $v[0 \ldots m]$ зберігають потенціал. Спочатку вони встановлюються в нуль, що узгоджується з матрицею з нульових рядків (Зауважимо, що для цієї реалізації неважливо, чи містить матриця $A$ від'ємні числа).

Масив $p[0 \ldots m]$ містить парування: для кожного стовпця $j = 1 \ldots m$ він зберігає номер $p[j]$ вибраного рядка (або $0$, якщо ще нічого не вибрано). Для зручності реалізації $p[0]$ вважається рівним номеру поточного рядка.

Масив $minv[1 \ldots m]$ містить для кожного стовпця $j$ допоміжні мінімуми, потрібні для швидкого перерахунку потенціалу, як описано вище.

Масив $way[1 \ldots m]$ містить інформацію про те, де ці мінімуми досягаються, щоб ми могли пізніше відновити доповнювальний шлях. Зауважимо, що для відновлення шляху достатньо зберігати лише значення стовпців, оскільки номери рядків можна взяти з парування (тобто з масиву $p$). Таким чином, $way[j]$ для кожного стовпця $j$ містить номер попереднього стовпця в шляху (або $0$, якщо такого немає).

Сам алгоритм є зовнішнім циклом по рядках матриці, всередині якого розглядається $i$-й рядок матриці. Перший цикл do-while виконується, доки не буде знайдено вільний стовпець $j0$. Кожна ітерація циклу позначає відвіданим новий стовпець з номером $j0$ (обчислений на попередній ітерації; а спочатку рівний нулю — тобто ми починаємо з фіктивного стовпця), а також новий рядок $i0$ — суміжний з ним у паруванні (тобто $p[j0]$; а спочатку, коли $j0=0$, береться $i$-й рядок). Через появу нового відвіданого рядка $i0$ потрібно відповідно перерахувати масив $minv$ та $\Delta$. Якщо $\Delta$ оновлюється, то стовпець $j1$ стає мінімумом, який було досягнуто (зауважимо, що за такої реалізації $\Delta$ може виявитися рівним нулю, що означає, що потенціал не можна змінити на поточному кроці: вже є новий досяжний стовпець). Після цього перераховуються потенціал та масив $minv$. Наприкінці циклу "do-while" ми знайшли доповнювальний шлях, що закінчується в стовпці $j0$, який можна "розгорнути" за допомогою масиву предків $way$.

Константа INF — це "нескінченність", тобто деяке число, яке очевидно більше за всі можливі числа у вхідній матриці $A$.

C++

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<bool> 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);
}

Python

# Матриця A 1-індексована: A[1..n][1..m], n <= m.
# Масиви u, v, p, way, minv, used також 1-індексовані (індекс 0 — фіктивний).
u = [0] * (n + 1)
v = [0] * (m + 1)
p = [0] * (m + 1)
way = [0] * (m + 1)
for i in range(1, n + 1):
    p[0] = i
    j0 = 0
    minv = [INF] * (m + 1)
    used = [False] * (m + 1)
    while True:                       # do-while: шукаємо доповнювальний шлях
        used[j0] = True
        i0 = p[j0]
        delta = INF
        j1 = -1
        for j in range(1, m + 1):
            if not used[j]:
                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 j in range(0, m + 1):     # перерахунок потенціалу та minv
            if used[j]:
                u[p[j]] += delta
                v[j] -= delta
            else:
                minv[j] -= delta
        j0 = j1
        if p[j0] == 0:
            break
    while True:                       # розгортаємо доповнювальний шлях через way
        j1 = way[j0]
        p[j0] = p[j1]
        j0 = j1
        if j0 == 0:
            break

TypeScript

// Матриця A 1-індексована: A[1..n][1..m], n <= m.
// Масиви u, v, p, way, minv, used теж 1-індексовані (індекс 0 — фіктивний).
const u = new Array<number>(n + 1).fill(0);
const v = new Array<number>(m + 1).fill(0);
const p = new Array<number>(m + 1).fill(0);
const way = new Array<number>(m + 1).fill(0);
for (let i = 1; i <= n; ++i) {
  p[0] = i;
  let j0 = 0;
  const minv = new Array<number>(m + 1).fill(INF);
  const used = new Array<boolean>(m + 1).fill(false);
  do {
    used[j0] = true;
    const i0 = p[j0];
    let delta = INF;
    let j1 = -1;
    for (let j = 1; j <= m; ++j) {
      if (!used[j]) {
        const 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 (let 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 {
    const j1 = way[j0];
    p[j0] = p[j1];
    j0 = j1;
  } while (j0);
}

Go

// Матриця A 1-індексована: A[1..n][1..m], n <= m.
// Зрізи u, v, p, way, minv, used теж 1-індексовані (індекс 0 — фіктивний).
u := make([]int, n+1)
v := make([]int, m+1)
p := make([]int, m+1)
way := make([]int, m+1)
for i := 1; i <= n; i++ {
    p[0] = i
    j0 := 0
    minv := make([]int, m+1)
    used := make([]bool, m+1)
    for j := range minv {
        minv[j] = INF
    }
    for { // do-while: шукаємо доповнювальний шлях
        used[j0] = true
        i0 := p[j0]
        delta := INF
        j1 := -1
        for j := 1; j <= m; j++ {
            if !used[j] {
                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 j := 0; j <= m; j++ { // перерахунок потенціалу та minv
            if used[j] {
                u[p[j]] += delta
                v[j] -= delta
            } else {
                minv[j] -= delta
            }
        }
        j0 = j1
        if p[j0] == 0 {
            break
        }
    }
    for { // розгортаємо доповнювальний шлях через way
        j1 := way[j0]
        p[j0] = p[j1]
        j0 = j1
        if j0 == 0 {
            break
        }
    }
}

Щоб відновити відповідь у більш звичному вигляді, тобто знайти для кожного рядка $i = 1 \ldots n$ номер $ans[i]$ вибраного в ньому стовпця, можна зробити так:

C++

vector<int> ans (n+1);
for (int j=1; j<=m; ++j)
    ans[p[j]] = j;

Python

ans = [0] * (n + 1)
for j in range(1, m + 1):
    ans[p[j]] = j

TypeScript

const ans = new Array<number>(n + 1).fill(0);
for (let j = 1; j <= m; ++j) ans[p[j]] = j;

Go

ans := make([]int, n+1)
for j := 1; j <= m; j++ {
    ans[p[j]] = j
}

Вартість парування можна просто взяти як потенціал нульового стовпця (взятий з протилежним знаком). Справді, як видно з коду, $-v[0]$ містить суму всіх значень $\Delta$, тобто сумарну зміну потенціалу. Хоча відразу могло змінитися кілька значень $u[i]$ та $v[j]$, сумарна зміна потенціалу точно дорівнює $\Delta$, оскільки доки існує доповнювальний шлях, кількість досяжних рядків точно на один більша за кількість досяжних стовпців (лише поточний рядок $i$ не має "пари" у вигляді відвіданого стовпця):

C++

int cost = -v[0];

Python

cost = -v[0]

TypeScript

const cost = -v[0];

Go

cost := -v[0]

Зв'язок з алгоритмом послідовних найкоротших шляхів

Угорський алгоритм можна розглядати як алгоритм послідовних найкоротших шляхів, адаптований до задачі про призначення. Не вдаючись у деталі, дамо інтуїтивне розуміння зв'язку між ними.

Алгоритм послідовних шляхів використовує модифіковану версію алгоритму Джонсона як техніку перезважування. Вона ділиться на чотири кроки:

• Використовуємо алгоритм Беллмана–Форда, починаючи зі стоку $s$, і для кожного вузла знаходимо мінімальну вагу $h(v)$ шляху від $s$ до $v$.

Для кожного кроку основного алгоритму:

• Перезважуємо ребра вихідного графа так: $w(u,v) \gets w(u,v)+h(u)-h(v)$.
• Використовуємо алгоритм Дейкстри для пошуку підграфа найкоротших шляхів вихідної мережі.
• Оновлюємо потенціали для наступної ітерації.

З огляду на цей опис можна помітити сильну аналогію між $h(v)$ та потенціалами: можна перевірити, що вони рівні з точністю до сталого зсуву. Крім того, можна показати, що після перезважування множина всіх ребер нульової ваги є підграфом найкоротших шляхів, у якому основний алгоритм намагається збільшити потік. Це саме відбувається і в угорському алгоритмі: ми створюємо підграф із жорстких ребер (тих, для яких величина $A[i][j]-u[i]-v[j]$ дорівнює нулю) і намагаємося збільшити розмір парування.

На кроці 4 оновлюються всі $h(v)$: щоразу, коли ми змінюємо потокову мережу, ми маємо гарантувати, що відстані від джерела правильні (інакше на наступній ітерації алгоритм Дейкстри може дати збій). Це нагадує оновлення, яке виконується над потенціалами, але в цьому випадку вони збільшуються неоднаково.

Щоб глибше зрозуміти потенціали, зверніться до цієї статті.

Приклади задач

Ось кілька прикладів, пов'язаних із задачею про призначення, від зовсім тривіальних до менш очевидних:

• Дано двочастковий граф, у якому потрібно знайти максимальне парування з мінімальною вагою (тобто, насамперед максимізується розмір парування, а вже потім мінімізується його вартість).
  Щоб розв'язати її, ми просто будуємо задачу про призначення, ставлячи число "нескінченність" на місце відсутніх ребер. Після цього розв'язуємо задачу угорським алгоритмом і вилучаємо з відповіді ребра нескінченної ваги (вони могли потрапити у відповідь, якщо задача не має розв'язку у вигляді досконалого парування).

• Дано двочастковий граф, у якому потрібно знайти максимальне парування з максимальною вагою.
  Розв'язок знову очевидний: усі ваги треба помножити на мінус одиницю.

• Задача виявлення рухомих об'єктів на зображеннях: було зроблено два знімки, унаслідок чого отримано дві множини координат. Потрібно співвіднести об'єкти на першому й другому зображеннях, тобто визначити для кожної точки другого зображення, якій точці першого зображення вона відповідала. При цьому потрібно мінімізувати суму відстаней між зіставленими точками (тобто ми шукаємо розв'язок, у якому об'єкти сумарно подолали найкоротший шлях).
  Для розв'язання ми просто будуємо й розв'язуємо задачу про призначення, де вагами ребер є евклідові відстані між точками.

• Задача виявлення рухомих об'єктів локаторами: є два локатори, які не можуть визначити положення об'єкта в просторі, а лише його напрямок. Обидва локатори (розташовані в різних точках) отримали інформацію у вигляді $n$ таких напрямків. Потрібно визначити положення об'єктів, тобто визначити очікувані положення об'єктів та відповідні їм пари напрямків так, щоб сума відстаней від об'єктів до променів напрямків була мінімальною.
  Розв'язок: знову ж таки, ми просто будуємо й розв'язуємо задачу про призначення, де вершинами лівої частини є $n$ напрямків від першого локатора, вершинами правої частини — $n$ напрямків від другого локатора, а вагами ребер — відстані між відповідними променями.

• Покриття орієнтованого ациклічного графа шляхами: дано орієнтований ациклічний граф, у якому потрібно знайти найменшу кількість шляхів (за рівності — з найменшою сумарною вагою) так, щоб кожна вершина графа лежала рівно в одному шляху.
  Розв'язок полягає в тому, щоб побудувати з даного графа відповідний двочастковий граф і знайти в ньому максимальне парування мінімальної ваги. Докладніше див. окрему статтю.

• Розфарбовування дерева. Дано дерево, у якому кожна вершина, окрім листків, має рівно $k-1$ дітей. Потрібно вибрати для кожної вершини один із $k$ доступних кольорів так, щоб жодні дві суміжні вершини не мали однакового кольору. Крім того, для кожної вершини та кожного кольору відома вартість фарбування цієї вершини цим кольором, і потрібно мінімізувати сумарну вартість.
  Щоб розв'язати цю задачу, скористаємося динамічним програмуванням. А саме, навчімося обчислювати величину $d[v][c]$, де $v$ — номер вершини, $c$ — номер кольору, а сама величина $d[v][c]$ — це мінімальна вартість, потрібна для розфарбовування всіх вершин у піддереві з коренем у $v$, причому саму вершину $v$ — кольором $c$. Щоб обчислити таку величину $d[v][c]$, необхідно розподілити решту $k-1$ кольорів між дітьми вершини $v$, а для цього необхідно побудувати й розв'язати задачу про призначення (у якій вершинами лівої частини є кольори, вершинами правої частини — діти, а вагами ребер — відповідні значення $d$).
  Таким чином, кожне значення $d[v][c]$ обчислюється за допомогою розв'язку задачі про призначення, що зрештою дає асимптотику $\mathcal{O}(nk^4)$.

• Якщо в задачі про призначення ваги знаходяться не на ребрах, а на вершинах, причому лише на вершинах однієї частини, то немає потреби застосовувати угорський алгоритм: достатньо відсортувати вершини за вагою й запустити звичайний алгоритм Куна.

• Розглянемо наступний окремий випадок. Нехай кожній вершині лівої частини призначено деяке число $\alpha[i]$, а кожній вершині правої частини — $\beta[j]$. Нехай вага будь-якого ребра $(i,j)$ дорівнює $\alpha[i]\cdot \beta[j]$ (числа $\alpha[i]$ та $\beta[j]$ відомі). Розв'яжіть задачу про призначення.
  Щоб розв'язати її без угорського алгоритму, спочатку розглянемо випадок, коли обидві частини мають по дві вершини. У цьому випадку, як легко бачити, краще з'єднувати вершини у зворотному порядку: з'єднати вершину з меншим $\alpha[i]$ з вершиною з більшим $\beta[j]$. Це правило легко узагальнюється на довільну кількість вершин: потрібно відсортувати вершини першої частини в порядку зростання значень $\alpha[i]$, другої частини — в порядку спадання значень $\beta[j]$, і з'єднати вершини попарно в цьому порядку. Таким чином, ми отримуємо розв'язок зі складністю $\mathcal{O}(n\log n)$.

• Задача про потенціали. Дано матрицю $A[1 \ldots n][1 \ldots m]$, потрібно знайти два масиви $u[1 \ldots n]$ та $v[1 \ldots m]$ такі, що для будь-яких $i$ та $j$ виконується $u[i] + v[j] \leq a[i][j]$, а сума елементів масивів $u$ та $v$ є максимальною.
  Знаючи угорський алгоритм, розв'язати цю задачу буде неважко: угорський алгоритм якраз і знаходить такий потенціал $u, v$, що задовольняє умову задачі. З іншого боку, без знання угорського алгоритму розв'язати таку задачу видається майже неможливим.

  Зауваження

  Цю задачу також називають двоїстою задачею до задачі про призначення: мінімізація сумарної вартості призначення еквівалентна максимізації суми потенціалів.

Література

• Ravindra Ahuja, Thomas Magnanti, James Orlin. Network Flows [1993]

• Harold Kuhn. The Hungarian Method for the Assignment Problem [1955]

• James Munkres. Algorithms for Assignment and Transportation Problems [1957]

Задачі для практики

• UVA - Crime Wave - The Sequel

• UVA - Warehouse

• SGU - Beloved Sons

• UVA - The Great Wall Game

• UVA - Jogging Trails

Відеоматеріали