グラフ+bitDP:巡回セールスマン問題

今回は巡回セールスマン問題(Traveling Salesman Problem: TSP)のコードを書きます。グラフ上のすべてのノードをめぐる最小距離を求める問題です。どちらかというとグラフの探索アルゴリズムというよりかは、DP 的側面の強いアルゴリズムという印象。

素朴に書こうとすると、ノード数を  n としたとき、

  1. グラフ上のすべての2点間の最短距離を求める(素朴に書くと  O(n^{2}) の計算量)

  2. 巡回する順序を全パターン総当たりする( O(n!) の計算量)

総和としては、 O(n^{2} + n!) = O(n!) の計算量となります。

今回ご紹介するスクリプトを使用すると、計算量は  O(n^{2} \times 2^{n}) になります。まあまあましになりますね。 n自然数としたとき、 n \geq 8 O(n^{2} \times 2^{n}) の方が小さくなります。

こちらが出発点は0番ノード固定、かつ出発点に戻ってくる必要のない場合の巡回セールスマン問題を解くスクリプトです。出発点が0番ノード以外であったり、出発点に戻ってくる必要がある場合はスクリプトを改変すれば対応できるので、臨機応変にお使いください。

import math

class TravelingSalesman:

    def __init__(self, costs):
        self.costs = costs
        self.n = len(costs)
        self.s = 2 ** self.n

    def solve(self):
        dp = [[math.inf for _ in range(self.n) ] for _ in range(self.s)]

        dp[1][0] = 0

        for s in range(1, self.s):
            for u in range(self.n):
                u_bit = 2 ** u 
                if s & u_bit == u_bit:
                    for v in range(self.n):
                        s2 = s | 2 ** v
                        dp[s2][v] = min(dp[s2][v], dp[s][u] + self.costs[u][v])

        return min(dp[self.s - 1])


使い方は以下のようになるかと思います。

n = int(input())
m = int(input())

costs = [[math.inf for _ in range(n)] for _ in range(n)]
for _ in range(m):
    ai, bi, ci = map(int, input().split())
    costs[ai][bi] = ci
    costs[bi][ai] = ci
tsp = TravelingSalesman(costs)
print(tsp.solve())


以下スクリプトの解説です。

import math

TravelingSalesman クラスの solve メソッド内で、DP テーブルを math.inf で初期化したかったがための import です。解く DP によっては初期値が math.inf でないこともあるので、適宜書き換えてください。

class TravelingSalesman:

    def __init__(self, costs):
        self.costs = costs
        self.n = len(costs)
        self.s = 2 ** self.n

コンストラクタです。使い方にもある通り、引数には重み付きグラフの行列表現を渡してください。渡されたグラフの配列長(= グラフのノード数)と、それを2の指数に持ってきたものを格納しておきます。役割は後続のステップにて。

    def solve(self):
        dp = [[math.inf for _ in range(self.n) ] for _ in range(self.s)]

        dp[1][0] = 0

        for s in range(1, self.s):
            for u in range(self.n):
                u_bit = 2 ** u 
                if s & u_bit == u_bit:
                    for v in range(self.n):
                        s2 = s | 2 ** v
                        dp[s2][v] = min(dp[s2][v], dp[s][u] + self.costs[u][v])

        return min(dp[self.s - 1])

巡回セールスマン問題を解くメソッドです。細かく分解してみていきましょう。

        dp = [[math.inf for _ in range(self.n) ] for _ in range(self.s)]

DP テーブルを初期化します。
テーブルサイズは、グラフのノード数を  n として  n \times 2^{n} です。
1次元目( 2^{n} の方)は、今までに到達したノードを表現します。ノードが  n 個ある場合、そのうちどれに到達したかは boolean を  n 個保持すればよいので、 2^{n} で表現可能です。
2次元目( n の方)は、今どのノードにいるかを表現します。ノードが  n 個なので、 n で表現可能です。

        dp[1][0] = 0

DP テーブルに開始点の情報を入力します。
0番ノードから開始するため、 dp \lbrack (000 \ldots 001)_{2} \rbrack \lbrack 0 \rbrackに累計コスト  0 を設定します。

        for s in range(1, self.s):
            for u in range(self.n):
                u_bit = 2 ** u 
                if s & u_bit == u_bit:
                    for v in range(self.n):
                        s2 = s | 2 ** v
                        dp[s2][v] = min(dp[s2][v], dp[s][u] + self.costs[u][v])

DP の遷移です。配る DP を行います。テーブルを辿っていく順番は、
 dp \lbrack (000 \ldots 001)_{2} \rbrack \lbrack 0 \rbrack,
 dp \lbrack (000 \ldots 001)_{2} \rbrack \lbrack 1 \rbrack,
 dp \lbrack (000 \ldots 001)_{2} \rbrack \lbrack 2 \rbrack,
 \vdots
 dp \lbrack (000 \ldots 001)_{2} \rbrack \lbrack n \rbrack,
 dp \lbrack (000 \ldots 010)_{2} \rbrack \lbrack 0 \rbrack,
 dp \lbrack (000 \ldots 010)_{2} \rbrack \lbrack 1 \rbrack,
 \vdots
 dp \lbrack (000 \ldots 010)_{2} \rbrack \lbrack n \rbrack,
 \vdots
 dp \lbrack (111 \ldots 111)_{2} \rbrack \lbrack n \rbrack
です。

                u_bit = 2 ** u 
                if s & u_bit == u_bit:

 s u の組み合わせによってはありえない組み合わせがあることもあるので、そこを省略します。例えば  dp \lbrack (000 \ldots 010)_{2} \rbrack \lbrack 0 \rbrack は、1次元目が「1番のノードにのみ到達済み」と言っているのに2次元目が「今0番ノードにいる」と言っており不整合が起きています。そのような DP テーブル内の値からの遷移は考えないものとします。

                    for v in range(self.n):
                        s2 = s | 2 ** v
                        dp[s2][v] = min(dp[s2][v], dp[s][u] + self.costs[u][v])

 v が次に移動する先のノード番号、 s2 v 到達後の到達済みノードの集合のバイナリ表現です。もし、到達済みノード  s の状態で  u から  v に移動した結果の総コストが、今  dp \lbrack s2 \rbrack \lbrack v \rbrack に格納されている値より小さかった場合、  dp \lbrack s2 \rbrack \lbrack v \rbrack を更新します。

        return min(dp[self.s - 1])

self.s - 1 の値は  (111 \ldots 111)_{2} です。そして、 dp \lbrack (111 \ldots 111)_{2} \rbrack \lbrack i \rbrack に格納されている値は、「全ノード到達済みで、現在  i 番目のノードにいるときの最小総コスト」です。
そのため、この式で返される値は、全ノードを辿る場合の最小総コストとなります。