PythonでbitDPを使い巡回セールスマン問題を解く #Python

PythonでbitDPを使って巡回セールスマン問題を解く記事を探したところ、あまり見当たらなかったため備忘録がてら投稿します。また、重み付き有向グラフについては知っていることを前提としています。
注意点として、自分にとってわかりやすいかを基準にしているため、厳密でない、わかりにくい場合があります...

問題リンク

AOJ Traveling Salesman Problem
もしAtCoderにも同様の問題がありましたら教えていただけると幸いです...

問題の要点

最初に問題の要点だけ書いておきます。

まず、タイムリミットは1 sec, メモリの上限は131072 KBです。2 secではありません。

重み付きの有向グラフが与えられる。このとき、以下の二つを満たす経路のうち、通った辺の重みの総和が最小のものを求めよ（最小の重みを出力する）。ただし、そのような経路が存在しない場合は-1を出力せよ。
・ある頂点から出発し、最後には戻ってくる経路である
・すべての頂点をちょうど1回通る経路である（すべての頂点を通る必要があるが、2回以上同じ頂点を通ってはいけない）

入力形式は

V E
s₀ t₀ d₀
s₁ t₁ d₁
:
s_E-1 t_E-1 d_E-1

となっており、Vは頂点の数、Eは辺の数、s_i, t_i, d_iは頂点s_iから頂点t_iに向かって重みがd_iの辺が張られていることを表しています。

制約は

2 ≤ V ≤ 15
0 ≤ E ≤ V(V-1)
0 ≤ s_i, t_i ≤ V-1 (s_i ≠ t_i)
0 ≤ d_i ≤ 1000
与えられるグラフに多重辺はない

となっています。Eの制約は明記されていませんが、与えられるグラフは多重辺を持たないと書いてあるので、V(V-1)（無向グラフの場合はV(V-1)/2だが、有向グラフなので、s_iとt_iを入れ替えた辺が存在する可能性もあるため/2する必要がない）であると言えます。（実際、E > V(V-1)の場合exitする処理を追加したコードでもACが取れました。また、E = 0の場合も同様にexitするコードを提出した際WAになったため、この制約でおそらく間違いはないでしょう）

bitDPって何？

bitDPは、集合に対するDPと言われることがあります。少し具体的にいうと、N個の要素の順列の中で最も適切な順列を効率的に求めたい時に使うことが多いようです。N個の要素の順列を愚直に処理すると、順列はN!通りあるため、計算量はO(N×N!)となり、この問題はタイムリミットが1 secなので、PyPyではNが10のときが限界なってしまいます。しかし、例えばこの問題ではbitDPを用いると計算量をO(n^2×2^n)まで落とすことができ、Nが19まではギリギリ通すことができます。今回はN ≤ 15なので、余裕を持って通すことができます。
今回の問題では実際にどんなDPになるかというと、

DP[S][v] := 集合の全体のうちの部分集合Sの順列の中で、vが末項となるもののうち、最も最適なもののコスト（順列の評価方法は、その順列の通りの順番で頂点を通った時の重みの総和の低さで、最適な順列とは条件を満たすすべての順列の中で最も重みの総和が低いもの）

となります。このとき、Sとはbitで表された集合です。集合をbitで表すのでbitDPと呼ばれているわけですね。ですので、Sの上限は2^V-1となります。また、末項で分ける理由として、例えば集合{1, 2, 5}の順列のうち末項が2であるものは{1, 5, 2}, {5, 1, 2}の2つですが、どちらとも状況としては「頂点1と5はすでに通っており、現在は頂点2にいる」と言えるため、このうち最もコストが低いものだけを考えることで効率化できるためです。

注意点

さて、プログラムに落とし込む前に1つ注意点があります。それはどの頂点から出発したかの判定の部分です。これは単純で、常に頂点0から出発したことにすればよいです。
例えば、上の図のような順番で頂点を通ったとき、出発地点が4であっても2であっても0であっても（勿論他の頂点であったとしても）、出発地点まで戻るまでのコストは変わりません。なので、出発地点を0で固定したとしても問題ありません。

実際のプログラム

まず実際のプログラムにコメントを追加したものを載せておきます。今回は配るDPを使用しました。また、①〜⑤はそれぞれプログラムの下に少し詳しく説明を書いています。

V, E = map(int,input().split())
G = [[float('inf')]*V for i in range(V)] # 存在しないパスはinfになるように、最初にすべてinfにしておく
for i in range(E):
    s,t,d = map(int,input().split())
    G[s][t] = d # s,tは0以上V-1以下なので、デクリメントの必要はない
dp = [[float('inf')]*V for i in range(2**V)] # dpの長さは2^V必要
dp[0][0] = 0 # 0から出発するのでdp[0][0]を0にしておく


for S in range(2**V): # Sは集合をbitで表している
    for v in range(V): # vは配られる側の要素を表している
        for u in range(V): # uは配る側の要素を表している
            if not (S >> u) & 1 and S != 0: # ①
                continue
            if (S >> v) & 1 == 0: # ②
                if dp[S][u] + G[u][v] < dp[S | (1 << v)][v]:
                    dp[S | (1 << v)][v] = dp[S][u] + G[u][v] # ③

# 全ての要素が含まれていて、末項が0のものの最小コストを出力する
# infだった場合は到達できないということなので-1を出力する
if dp[2**V-1][0] != float('inf'):
    print(dp[2**V-1][0])
else:
    print(-1)

このプログラムでは、vは配られる側の要素を表し、uは配る側の要素を表しています。今回、要素とは頂点のことであるため、v,uには0〜V-1が入ります。

①このプログラムではuが配る側の要素を表しているため、そもそもSのu番目のbitが立っていない、つまりu番目の頂点には到達していない場合、u番目の頂点から配ることはできません。しかし、今回は配るDPを使用しているため、sが0の時は例外的に処理をする必要があります。そうしないと、dp[0][0]以外全てinfのままになってしまいます。ですので、この部分は日本語で、「u番目の頂点に訪れている、あるいはどこの頂点にも訪れていない場合は続きの処理をし、そうでない場合は今回のループではこれ以上の処理をせず、次のループに移る」と表すことができます。

②今回、vは配られる側であり、Sのv番目のbitが立っているということはすでにv番目の頂点には訪れているということであるため、出発地点を除いた頂点を2回以上訪れてはいけないという条件から、処理をする必要はありません。ですので、この部分ではSのv番目のbitが立っていない場合のみ続きの処理をするようにしています。

③これは、dpを更新しているだけです。dp[S][v]にはbitで表された集合Sの順列のうち、末項がvであるものの最小コストが代入されているため、そのような順列のなかで、これまで記録されていた最小コストよりも少ないコストの順列を見つけることができた場合そのコストを記録しています。先にも説明しましたが、順列を保存する必要はなく、集合Sと末項vのみがわかればいいため、このようなdpリストになっています。

最後に

今回、初めての競プロ記事ということもありかなりゴチャゴチャした記事になってしまっています（極力スッキリはさせたつもりですが...）。ですが、もしどなたかの参考になったなら幸いです。

参考文献

今回の記事を書くにあたって、以下の二つの記事を参考にさせていただきました。ありがとうございます。

ビットDP(bit DP)の考え方 ~集合に対する動的計画法~

ビット演算 (bit 演算) の使い方を総特集！〜マスクビットから bit DP まで〜