ハンガリアン法(Hungarian Algorithm)
(flow/hungarian.cpp)

• View this file on GitHub
• Last update: 2026-03-08 22:25:54+09:00

説明

二部割当問題の最適解を $O(N^2 M)$ で求める。 N <= M の長方形行列に対応し、各行を異なる列へ 1 つずつ割り当てる。

できること

• hungarian<T, Minimize>(cost) コスト行列 cost に対する最適値、各行の割当先列、行側双対変数、列側双対変数を返す
• Minimize=true 最小コスト割当を求める
• Minimize=false 最大コスト割当を求める

返り値は tuple<T, vector<int>, vector<T>, vector<T>> で、順に次を表す。

• 最適値
• match[i]: 行 i を割り当てた列番号
• row[i]: 行側双対変数
• col[j]: 列側双対変数

Minimize=true では row[i] + col[j] <= cost[i][j] を満たし、最適値は双対和と一致する。

使い方

行数が列数以下のコスト行列を渡す。

vector<vector<long long>> cost = {
    {3, 5, 4},
    {6, 1, 9},
};
auto [ans, match, row, col] = hungarian<long long>(cost);

最大化したいなら hungarian<long long, false>(cost) を使う。

実装上の補足

空行列には 0 と空ベクトルを返す。 各行の長さが揃っていること、N <= M を assert する。

Verified with

• test/yosupo_assignment_hungarian.test.cpp

Code

template<class T, bool Minimize = true>
tuple<T, vector<int>, vector<T>, vector<T>> hungarian(const vector<vector<T>> &cost) {
    int n = cost.size();
    if (n == 0) return {T(0), {}, {}, {}};
    int m = cost[0].size();
    assert(n <= m);
    for (int i = 0; i < n; ++i) assert((int)cost[i].size() == m);

    vector<vector<T>> a(n + 1, vector<T>(m + 1));
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < m; ++j) {
            a[i + 1][j + 1] = Minimize ? cost[i][j] : -cost[i][j];
        }
    }

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

    for (int i = 1; i <= n; ++i) {
        p[0] = i;
        fill(minv.begin(), minv.end(), numeric_limits<T>::max());
        fill(used.begin(), used.end(), 0);
        int j0 = 0;
        while (p[j0] != 0) {
            used[j0] = 1;
            int i0 = p[j0], j1 = 0;
            T delta = numeric_limits<T>::max();
            for (int j = 1; j <= m; ++j) {
                if (used[j]) continue;
                T 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 (j0 != 0) {
            int j1 = way[j0];
            p[j0] = p[j1];
            j0 = j1;
        }
    }

    vector<int> match(n, -1);
    for (int j = 1; j <= m; ++j) {
        if (p[j] != 0) match[p[j] - 1] = j - 1;
    }

    vector<T> row(n), col(m);
    for (int i = 0; i < n; ++i) row[i] = u[i + 1];
    for (int j = 0; j < m; ++j) col[j] = v[j + 1];
    T ans = -v[0];
    if (!Minimize) {
        ans = -ans;
        for (int i = 0; i < n; ++i) row[i] = -row[i];
        for (int j = 0; j < m; ++j) col[j] = -col[j];
    }
    return {ans, match, row, col};
}

/**
 * @brief ハンガリアン法(Hungarian Algorithm)
 */

#line 1 "flow/hungarian.cpp"
template<class T, bool Minimize = true>
tuple<T, vector<int>, vector<T>, vector<T>> hungarian(const vector<vector<T>> &cost) {
    int n = cost.size();
    if (n == 0) return {T(0), {}, {}, {}};
    int m = cost[0].size();
    assert(n <= m);
    for (int i = 0; i < n; ++i) assert((int)cost[i].size() == m);

    vector<vector<T>> a(n + 1, vector<T>(m + 1));
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < m; ++j) {
            a[i + 1][j + 1] = Minimize ? cost[i][j] : -cost[i][j];
        }
    }

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

    for (int i = 1; i <= n; ++i) {
        p[0] = i;
        fill(minv.begin(), minv.end(), numeric_limits<T>::max());
        fill(used.begin(), used.end(), 0);
        int j0 = 0;
        while (p[j0] != 0) {
            used[j0] = 1;
            int i0 = p[j0], j1 = 0;
            T delta = numeric_limits<T>::max();
            for (int j = 1; j <= m; ++j) {
                if (used[j]) continue;
                T 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 (j0 != 0) {
            int j1 = way[j0];
            p[j0] = p[j1];
            j0 = j1;
        }
    }

    vector<int> match(n, -1);
    for (int j = 1; j <= m; ++j) {
        if (p[j] != 0) match[p[j] - 1] = j - 1;
    }

    vector<T> row(n), col(m);
    for (int i = 0; i < n; ++i) row[i] = u[i + 1];
    for (int j = 0; j < m; ++j) col[j] = v[j + 1];
    T ans = -v[0];
    if (!Minimize) {
        ans = -ans;
        for (int i = 0; i < n; ++i) row[i] = -row[i];
        for (int j = 0; j < m; ++j) col[j] = -col[j];
    }
    return {ans, match, row, col};
}

/**
 * @brief ハンガリアン法(Hungarian Algorithm)
 */

Back to top page