Gale-Shapleyアルゴリズム実装の高速化

　こんにちは、リブセンスでデータサイエンティストをしている北原です。前回に続き今回もGale-Shapleyアルゴリズムを扱います。前回の記事で紹介した実装のボトルネックを把握し、少し改良することで計算量を大幅に削減します。コードはJuliaです。なお、今回の記事は下記の前回記事を読んでいることを前提として書かれています。

analytics.livesense.co.jp

ボトルネック

　まず、ボトルネックを見つけましょう。以下は前回記事で紹介した実装です。よく見ると一見簡潔に書かれているものの計算時間がかかりそうなところが2箇所見つかります。

function gs(user1_prefs, user2_prefs, unmatched_user1_stack)
    matched_user2 = fill(-1, length(user2_prefs))
    while length(unmatched_user1_stack) > 0
        proposer = pop!(unmatched_user1_stack)
        if !isempty(user1_prefs[proposer])
            proposed = popfirst!(user1_prefs[proposer])
            # ボトルネック1（選好リストに含まれているかの判定）
            if proposer in user2_prefs[proposed]
                if matched_user2[proposed] < 0
                    matched_user2[proposed] = proposer
                else
                    # ボトルネック2(選好順位の取得)
                    current_user1_rank = findfirst(isequal(proposer), user2_prefs[proposed])
                    matched_user1_rank = findfirst(isequal(matched_user2[proposed]), user2_prefs[proposed])
                    if current_user1_rank < matched_user1_rank
                        push!(unmatched_user1_stack, matched_user2[proposed])
                        matched_user2[proposed] = proposer
                    else
                        push!(unmatched_user1_stack, proposer)
                    end
                end
            else
                push!(unmatched_user1_stack, proposer)
            end
        end
    end
    return matched_user2
end

　一つ目のボトルネックは、in演算子を利用して、プロポーズ中の人のIDがプロポーズされる側の選好リストに含まれているかをチェックするところです。マッチング可能な人数$N$が増えると選好リストの探索量も増えるので計算量は$O(N)$と考えられます。

　二つ目のボトルネックは、findfirst()を利用して現在プロポーズ中の人とマッチング済みの人の選好順位を比較するところです。ここも計算量は$O(N)$と考えられます。

　いずれのボトルネックも、プロポーズされる人の選好リストを探索するところです。プロポーズ中の人のIDを与えて、(1)選好リストに含まれているかを判定する、(2)選好順位を取得する、という2つのことを$O(1)$でできればボトルネックが解消できます。$O(1)$での探索というとハッシュが思い浮かびますが、キーが数値ですし、マッチング可能な人数が増えると衝突も増えそうです。

選好順位判定行列

　今回の実装では、ボトルネックの解消のために、プロポーズされる人からみた選好順位を判定する行列を利用しています。非常にシンプルな方法ですが計算時間におよぼす効果は大きいです。

　選好順位を判定する行列では、プロポーズされる人のIDを行に、プロポーズする人のIDを列に割り当て、プロポーズされる人からみたプロポーズする人の選好順位を要素としています。なお、行と列は逆でも問題ありません。このような行列を使うことで選好順位は$O(1)$で判定できます。さらに、プロポーズする人が選好リストに含まれていない場合は選好順位ではなく型の最大値を入れます。これにより、プロポーズ中の人がプロポーズされる人の選好リストに含まれているかは、要素が型の最大値より小さいかで$O(1)$で判定できます。

　簡単な例を示します。プロポーズされる人のIDが1から4まであり、選好の高い順にIDを並べたプロポーズされる人の選好リストが[[3, 2, 1], [4, 2], [2, 3, 1, 4]]だったとすると選好順位判定行列は次のようになります。ID2のプロポーズされる人の選好リストにID3が含まれていないことは、ID2の行のID3の列をチェックすることで判定できます。ID1のプロポーズされる人からみて、ID3のほうがID2より選好順位が高いことは、ID1の行のID2の列とID3の列の順位を比較することから判定できます。

	ID1	ID2	ID3	ID4
ID1	3	2	1	MAX
ID2	MAX	2	MAX	1
ID3	3	1	2	4

　選好順位判定行列のデメリットはメモリ使用量です。プロポーズする人とされる人のID種類数と同じ行数と列数が必要になるので、メモリ使用量は$O(N^2)$となります。そのため、メモリ使用量に制約がある場合は、あえて選好順位判定テーブルを使わないほうが現実的なこともあります。また、選好リストに含まれないIDが極めて多い場合は、前節で指摘したボトルネックの影響が小さくなるので選好順位判定テーブルを使わなくても計算量の急増は抑制されます。データの特徴によって適宜使い分けるのが適切です。

実装

　では、選好順位判定行列を使ってGale-Shapleyアルゴリズムを実装しましょう。サンプルデータの作成やコマンド実行は前回記事の内容とほぼ同じなので省略します。

　以下が選好順位判定行列を使った実装です。make_rank_matrix()が選好順位判定行列を作成する関数です。選好順に行列の対応する要素に順位を入れるだけですが、計算量は$O(N^2)$なので注意が必要です。 gs()の変更箇所は引数指定とボトルネック箇所2点のみです。

# gs2.jl
user1_fname = ARGS[1]
user2_fname = ARGS[2]

function read_pref_file(fname)
    user_prefs = (Vector{Int32})[]
    id_vec = (Int32)[]
    f = open(fname, "r")
    for l_ in eachline(f)
        l = parse(Int32, l_)
        if l != -1
            push!(id_vec, l)
        else
            push!(user_prefs, id_vec)
            id_vec = (Int32)[]
        end
    end
    close(f)
    return user_prefs
end

@time user1_prefs = read_pref_file(user1_fname)
@time user2_prefs = read_pref_file(user2_fname)

# マッチ相手のいないuser1のスタック
@time unmatched_user1_stack = collect(1:length(user1_prefs))

# 低優先ユーザーの高優先ーザーへの選好順位判定行列
# 各行： 低優先ユーザー
# 各列： 高優先ユーザー
REJECT_RANK = typemax(Int32)
function make_rank_matrix(user2_prefs, REJECT_RANK)
    user2_pref_rank = fill(REJECT_RANK, length(user2_prefs), length(user1_prefs))
    rank = 1
    for i in eachindex(user2_prefs)
        for j in user2_prefs[i]
            user2_pref_rank[i, j] = rank
            rank += 1
        end
        rank = 1
    end
    return user2_pref_rank
end
@time user2_pref_rank = make_rank_matrix(user2_prefs, REJECT_RANK)

function gs(user1_prefs, user2_pref_rank, unmatched_user1_stack, REJECT_RANK)
    # user2のマッチ相手を格納する配列（未マッチは-1）
    matched_user2 = fill(-1, size(user2_pref_rank)[1])

    while length(unmatched_user1_stack) > 0
        proposer = pop!(unmatched_user1_stack)
        if !isempty(user1_prefs[proposer])
            proposed = popfirst!(user1_prefs[proposer])
            if user2_pref_rank[proposed, proposer] < REJECT_RANK
                # proposerがuser2の選好リストに存在する場合
                if matched_user2[proposed] < 0
                    # user2が未マッチの場合
                    matched_user2[proposed] = proposer
                else
                    # user2がマッチ済みの場合
                    if user2_pref_rank[proposed, proposer] < user2_pref_rank[proposed, matched_user2[proposed]]
                        # proposerの方が優先時順位が高い場合
                        push!(unmatched_user1_stack, matched_user2[proposed])
                        matched_user2[proposed] = proposer
                    else
                        # proposerの方が優先時順位が低い場合
                        push!(unmatched_user1_stack, proposer)
                    end
                end
            else
                # proposerがuser2の選好リストに存在しない場合
                push!(unmatched_user1_stack, proposer)
            end
        end
    end
    return matched_user2
end
@time matched_user2 = gs(user1_prefs, user2_pref_rank, unmatched_user1_stack, REJECT_RANK)
# println(matched_user2)

実行時間の計測

　筆者の環境で@timeを使ってgs()の実行時間を計測し、どの程度改善されたか調べました。選好順位判定行列の作成が新たなボトルネックになる可能性もあるのでmake_rank_matrix()の実行時間も計測します。また、選好順位判定行列を利用する前との比較用も行います。データ量の表記は前回の記事と同じく、プロポーズする側x人される側y人のとき、nxmyと表記しています。

　まず、選好順位判定行列を利用する前との比較結果を以下に示します。make_rank_matrix()の実行時間は無視できないくらい大きいものの、それ以上にgs()の実行時間が短縮されています。make_rank_matrix()の実行時間を含めても改善後のほうが実行時間が短くなっており、データ量が増大するにつれて差が大きくなっていることが確認できます。

データ量	改善前のgs()実行時間([sec])	改善後のgs()実行時間([sec])	make_rank_matrix()実行時間([sec])
n2500m2000	6	0.117	0.07
n3750m3000	20.5	0.27	0.117
n5000m4000	47.6	0.5	0.21

　次に、より大きデータで実行時間を調べた結果を以下に示します。データ量が少ない場合は$O(N^2)$で実行時間が増えますが、データ量が多くなると$O(N^3)$を超える実行時間の増加が見られます。アルゴリズムだけでなく、メモリ使用量増加などの影響もありそうです。合計5万人前後の規模のデータ量になると、入力データのファイルサイズは一つあたりおよそ2.5〜4GBになるので、実行時間よりメモリ使用量が問題になってきます。合計5〜6万人以上のデータを扱う場合は、別のアプローチを検討するのが現実的と思われます。

データ量	改善後のgs()実行時間([sec])	make_rank_matrix()実行時間([sec])
n10000m8000	2.5	0.87
n15000m12000	7.17	3.06
n20000m16000	13.2	5.12
n25000m20000	31.2	13.81
n30000m24000	70.6	25

まとめ

　今回は前回紹介したGale-ShapleyアルゴリズムのJulia実装を改善し高速化する方法を紹介しました。プロポーズされる人からみた選好順位を判定する行列を導入し、わずかな修正をするだけで、計算時間が数十倍以上短縮されました。ボトルネックを発見し適切な対処を行うことで、高度な技術を使わなくとも大幅な性能改善が図れることがあります。