巡回セールスマン問題の答えを求めるぞ(心意気だけ)

これまでは物理計算について色々やってましたが、今回は巡回セールスマン問題をごり押しでなんとかできないかなと考えていたところ乱数を使ってできそうなことがあったのでプログラムしてみました。さながら、風車の巨人に立ち向かうドン・キホーテといったところですが。。。

そもそも巡回セールスマン問題って

ツイッターのフォロワーさんから聞いて知ったのですが、巡回セールスマン問題という数学(離散最適化)の有名な問題がありまして、その内容は以下です。

都市の数 n と，二つの都市 i と j の間を移動するのにかかる時間 dij が与えられる。このとき，全ての都市を巡ってもとの都市に戻ってくるような経路のうち，移動時間が最小なものを求めよ。
引用元 https://mathtrain.jp/tsp

　一筆書きの最短経路を求める問題なので一見非常に簡単に思えます。しかし、。総当たりで求めようとすると$(n-1)!/2$ 通りの検討が必要となるので、都市の数$n$が増えれば増えるほど、生きている間に計算が終わらなくなってしまう恐れがあるのです。(例えば$n=50000$とした時はおよそ$6.695 \times 10^{213231}/2$ 通りの検討が必要になるようです。数が大きすぎてよく分からん！)
　さて、どうやって効率的に最短経路を求めましょうかね。

巡回セールスマン問題を力技で屈服させるためのアイデア

　プログラムはどうとでもできるので、まずは最短経路をどうやって計算させるかを考えてみました。そこで思いついたのが、二つの都市$i$と$j$の間を移動するのにかかる時間$d_{ij}$を確率の重みづけに用いて、現在位置からなるべく時間のかからない次の位置を選択する割合が高くなるようにし試行を繰り返して、最も短かった経路を探し出すのはどうか、ということでした。言葉で言ってもわかりづらいので、以下に数式で説明します。

　都市の数を$n$として、$n$行$n$列の行列$\bf{X}$を作ります。$\bf{X}$の各成分$x_{ij}$は以下のようにします。

x_{ij} = \begin{cases} 0 \ (i = j) & \\ 1/d_{ij} \ (i \neq j) & \end{cases}

つまりは対角要素は$0$として、非対角要素には経過時間$d_{ij}$の逆数を格納します。従って、行列$\bf{X}$の各要素は移動時間が短ければ短いほど大きく、長ければ長いほど小さくなるような値となっています。また、後で分かりますが、0は移動できない場所を表すようになっています。この行列を重み付けとして使用するのが今回の肝です。

　さて、ではどのように重み付けに使用するかを、ある場所$i’$に到着したとして説明します。
　先ほどの行列$\bf{X}$の$i’$行を取り出し、それから$i’$行の既に通った経路の値を$0$とした後に、各要素の総和が$1$となるように調整したベクトルを$\bf{y}$とします。($\sum_0^n y_i =1 $)

　ここで、お得意の$0$~$1$の間の値のでるサイコロを振り、出た値$r$が$ \sum_0^{j'-1} y_j \le r < (\sum_0^{j'} y_j)$ の範囲であれば次の移動先を$j'$とします。$\bf{y}$の各要素の値は移動時間の逆数が元となっているので、移動時間が短い方になるべく移動しやすくなっている(重み付けができている)はずです！

　こうして移動を繰り返せば、そう多くない試行回数でいつかは最短経路を求められる気がする!?さて結果はどうでしょうか。

力技でコードを作成！

　例の如くpythonを使用して実装してみます。コードは以下のようになりました。

salesman_question


import numpy as np

# n行n列の行列xを与えると、移動位置のリストlist_k, 移動時間のリストと
# 合計list_t, np.sum(list_t)を返す。
def salesman(x):
    dx = x.copy()
    x_size = dx.shape[0]
    
    #開始場所を0とする
    list_k = np.array([0])
    #移動時間格納用のリスト
    list_t = np.array([])
    #あらかじめサイコロを振っておく
    list_p = np.random.rand(x_size-1)
    
    for p in list_p:
        #現在地の行を取り出す
        x_k = dx[int(list_k[list_k.shape[0]-1]), :]
        #既に通った経路の値を0とする
        x_k[list_k] = 0
        #各要素の総和を1にする
        x_choice = x_k/np.sum(x_k)
        #上記ベクトルから要素が0の場所を除く
        x_choice2 = x_choice[x_choice>0]
        #次の移動先の決定
        j = 0
        for num in x_choice2:
            if p >= j and p < j + num:
                index = np.where(x_choice == num)
                list_k = np.append(list_k,index)
                list_t = np.append(list_t, 1/x_k[index])
                break
            else:
                j += num
    #初期位置0に戻す
    list_k = np.append(list_k,list_k[0])
    list_t = np.append(list_t, 1/dx[list_k[list_k.shape[0]-2],list_k[0]])
    return list_k, list_t, np.sum(list_t)

さきほど説明した通りの内容ですね。うまくいくかはよく分からないですが、細かいことは置いていて実践してみましょう！実践あるのみ！

いざ実践!

　さて実践してみましょうか。いちいち経路を考えるのは面倒なので、今回はnumpy.random.rand()でランダムな行列$\bf{X}$を作成します。例えばこんな感じで$5 \times 5$の行列を作りました。

In:
n = 5
x = np.random.rand(n, n)
x[np.eye(x.shape[0], dtype=bool)] = 0
x

Out:
array([[0.        , 0.8748199 , 0.6166693 , 0.03229406, 0.84126589],
       [0.82060506, 0.        , 0.74087907, 0.45018579, 0.41579812],
       [0.49399052, 0.94061259, 0.        , 0.01392263, 0.87988552],
       [0.82709557, 0.96295332, 0.49125203, 0.        , 0.52172075],
       [0.96941202, 0.58100512, 0.60506291, 0.47758473, 0.        ]])

説明したとおり、対角要素は0でその他要素が移動時間の逆数と思ってください。

さて、100回ほど試行を繰り返してみて、最小となる結果を出力すると例えばこんな感じになりました。

In:
test_No = 100
list_K = np.zeros([test_No, x.shape[0]+1])
list_T = np.zeros([test_No, x.shape[0]])
list_sum_t = np.zeros(test_No)

for i in range(test_No):
    list_k,list_t, sum_t = salesman(x)
    list_K[i,:] = list_K[i,:] + list_k
    list_T[i,:] = list_T[i,:] + list_t
    list_sum_t[i] = sum_t

print(np.min(list_sum_t))
index = np.where(list_sum_t == np.min(list_sum_t))
print(index)
print(list_sum_t[index])
print(list_K[index,:])

Out:
6.932271276090873
(array([ 1, 10, 18, 23, 25, 30, 35, 38, 40, 52, 63, 64, 76, 86],
      dtype=int64),)
[6.93227128 6.93227128 6.93227128 6.93227128 6.93227128 6.93227128
 6.93227128 6.93227128 6.93227128 6.93227128 6.93227128 6.93227128
 6.93227128 6.93227128]
[[[0. 1. 2. 4. 3. 0.]
  [0. 1. 2. 4. 3. 0.]
  [0. 1. 2. 4. 3. 0.]
  [0. 1. 2. 4. 3. 0.]
  [0. 1. 2. 4. 3. 0.]
  [0. 1. 2. 4. 3. 0.]
  [0. 1. 2. 4. 3. 0.]
  [0. 1. 2. 4. 3. 0.]
  [0. 1. 2. 4. 3. 0.]
  [0. 1. 2. 4. 3. 0.]
  [0. 1. 2. 4. 3. 0.]
  [0. 1. 2. 4. 3. 0.]
  [0. 1. 2. 4. 3. 0.]
  [0. 1. 2. 4. 3. 0.]]]

今回の試行では、どうやら[0. 1. 2. 4. 3. 0.]の順番で行くのが最短経路となるようです。n=5くらいなら簡単に結果が出せたようですね。yattane!(冷静に考えたら試行回数100回は$5!/2 = 60$通りを超えていたので、20回の試行で計算しなおしたら同じ答えがでました。安心。)

そして結論へ……。

　そして、それから、調子に乗って$n$を増やしてみたところ、$n=10000$くらいまではそれなりに計算ができました。ただ、試行回数を稼ぎ解を正確にするには結局それなりの時間が必要そうなので今回はここまでとします。。。効率的にするのって難しいですね……。
　個人的には、私にしてはなかなか発想はよかったんじゃないかなと自負しています(自画自賛)。for文の中身をもう少し見直せばもうちょっと計算が早い気がするので、もうちょっとnumpyの勉強が必要だなと感じました。とても面白い課題でやりがいがありましたね。ありがとうございました。