Juliaで実装するハンガリアン法

Table of Content

概要
アルゴリズム
1. Juliaでの実装
2. ステップ3について
関連
参考

概要

N個のものをM個にそれぞれペアで関連づけるマッチングで，コストが最小になるように割り当てる
$O(N^3)$

アルゴリズム

行の最小値をその行から引く
- すべての行について行う
列の最小値をその列から引く
- すべての列について行う
0を各行, 各列から1つずつ選ぶことができるか判定 (行数個選択)
- 選べた場合それが解
0を含む行または列を取り除く
- 行数+列数が最小になるように選択する
- 最小数個だけ取り除くには，最大個数の0を含む行か列を選択していく。
残りの要素の最小値$m$を各要素から引く
最小値$m$を取り除いた行, 列が交わる要素に加える
3に戻る

Juliaでの実装

"
`M`: cost matrix
"""
function hungarian(M)
    # 各行の最小値で引く
    M = mapslices(row->row .- minimum(row), M, dims=2)
    # 各列の最小値で引く
    M = mapslices(col->col.- minimum(col), M, dims=1)
    # 各行各列に0が含まれるかcheck
    while true
        pairs = getpairs(M)
        # ペアが見つかれば終了
        !isnothing(pairs) && return pairs
        # ペアがないので行と列を消す
        zerocoords = findall(x->iszero(x), M)
        rows, cols = Int[], Int[]
        while 0 < length(zerocoords)
            target = pop!(zerocoords)
            nzerocol = count(coord->coord[2] == target[2], zerocoords) # 列の0
            nzerorow = count(coord->coord[1] == target[1], zerocoords) # 行の0
            if nzerocol < nzerorow
                zerocoords = filter(coord->coord[1] != target[1], zerocoords) # 同一行を取り除く
                push!(rows, target[1])
            else
                zerocoords = filter(coord->coord[2] != target[2], zerocoords) # 同一列を取り除く
                push!(cols, target[2])
            end
        end
        # 残りの領域
        indices = filter(coord->coord[1] ∉ rows && coord[2] ∉ cols, CartesianIndices(M))
        crosspoints = CartesianIndex.((Base.product(rows, cols)))
        x_min = minimum(M[indices])
        M[indices] .-= x_min
        M[crosspoints] .+= x_min
    end
end

function getpairs(M)
    selects = CartesianIndex[]
    # 0要素index取得
    indices = findall(x->iszero(x), M)
    nrow = size(M, 1)
    for i in 1:nrow
        zerocoords = filter(coord->coord[1] == i, indices)
        # zerocoordsがない場合はpair計算不可能
        if isempty(zerocoords)
            return nothing
        end
        # 行について，0がある列の0の個数を計算
        # 最小の0の個数の要素を選択する
        zeros = map(zerocoords) do target
            count(coord->coord[2] == target[2], indices)
        end
        # 選択する0要素
        selected = zerocoords[argmin(zeros)]
        push!(selects, selected)
        # 同一列を取り除く
            indices = filter(coord->coord[2]!=selected[2], indices)
    end
    return selects
end

# 例
testmat = [5 4 7 6; 6 7 3 2; 8 11 2 5; 9 8 6 8]
x = hungarian(testmat)

ライブラリは存在する
FugroRoames/Munkres.jl: Munkres algorithm for the optimal assignment problem (github.com)