シミュレーテッドアニーリングで巡回セールスマン問題

世の中にはシミュレーテッドアニーリングというメタヒューリスティックがります。焼きなまし法ともいわれます。例えば、ガラスについて説明します。ガラスは個体ではありません。結晶を作らないからです。高温の液体のガラスを徐々に冷やしていくと内部エネルギーが最低の状態になります。内部エネルギーが最低なのでガラスの分子が互に動けなくなります。これが私たちがよく見ているガラスです。このような現象をコンピューターで再現します。それがシミュレーテッドアニーリングです。量子アニーリングなんかもこの現象をもとにして実現されています。

#SAについてもう少し詳しく

普通の勾配法なんかは周囲を見渡して値が下がる方向へ最小値を探索しに行きます。しかし、それだと局所最小解に陥る場合があります。ですので、SAでは値が悪くなる場合でもある確率で更新を許します。これによって局所最小解に陥ることを防ぐことができます。具体的にどのような確率を使うかというと、

```math
P = e^{\frac{\Delta\beta}{T}}
```

というPを使います。Δβは値の変化量、Tは温度です。これは、最初に温度が大きい状態では高い確率で遷移を許すが、Tが徐々に小さくなるにつれて徐々に遷移し辛くなるということです。

#巡回セールスマン問題についてちょっとだけ

巡回セールスマン問題とはn個の地点を最短ですべての点を通る経路は何かということを考える問題です。全探索するとO(n!)かかってしまいますが、未だ世界中の私など遠く及ばない天才たちの知力の粋を集めてもなかなか全探索より大幅に効率の良い方法を見つけることができない状態です。

#ここからコードを書いていくよ

##まずいろいろ初期化

巡回セールスマン問題を解く上でグラフの表現が必要になってきます。今回は各辺の端点の番号を入れたスライスgraphとその点が座標のどこに位置しているかを示すplaneという変数を定義してあります。Goの乱数のシード値の初期値は決まっているので座標生成まではシード値を指定せずに同じものを生成してやります。

```go
package main

import (
	"time"
	"math/rand"
	"fmt"
	"math"
)

/*
シミュレーテッドアニーリングを実装
巡回セールスマン問題を解いてみる
*/
//辺を表す。辺の両端の点を書いとく{{0,1},{1,2},{2,0}}みたいに
var graph [][]int
//2次元平面上の点の位置を記録
var plane [][]int
//頂点数を表す
var vNum = 50
//温度
var t = 1000000
//初期化
func init(){
	graph = make([][]int,vNum)
	plane = make([][]int,vNum)
	for i:=0;i<vNum;i++{
		//辺を順番に入れる
		graph[i] = []int{i,(i + 1) % vNum}
		//頂点の座標を記録
		//31^2-1まで許してしまうとオーバーフローが怖そうなので
		//1000までにしときます。（三乗しても壊れない数）
		plane[i] = []int{rand.Intn(1001),rand.Intn(1001)}
	}
	rand.Seed(time.Now().UnixNano())
}
```

##コスト関数の定義
この組み合わせがどの程度あっているのかの指標が必要なのでコスト関数を定義します。今回は辺の長さとします。二つの関数に分かれてるのはコスト関数もうちょっと大きくなるかなって思ってたからです。

```go
//二つの頂点を渡すと間のユークリッド距離を出す関数
func measure(a int,b int) float64 {
	x := plane[a][0] - plane[b][0]
	y := plane[a][1] - plane[b][1]
	return math.Sqrt(float64(x*x) + float64(y*y))
}

//コスト関数
//今んとこ二つのユークリッド距離です。
//あ・え・て拡張性を残してます。
func cost(a int,b int) float64 {
	return measure(a,b)
}
```

##遷移する確率
上の方で書いた確率の式を関数にします。順当に小さくなる場合は1を出力します。

```go
//遷移する確率を返す
func prob(before float64,after float64) float64 {
	if before >= after{
		return 1.0
	}
	return math.Exp((before - after)/float64(t))
}
```

##関数をまとめる

上で書いた関数たちをまとめます。まず、ランダムに重複を許さずに二つの辺を選び、頂点を交換しない場合とした場合を考えます。遷移確率に従って頂点を交換するか交換しないかを決めます。これをやります。同じ頂点で辺を結ぶ（つまり、頂点１から頂点１へ辺をつなぐ）ということはないようにはじいてやるということも忘れずに

```go
//2辺をランダムに選んで遷移確率に従って辺の組み合わせを変える
func change(){
	a := rand.Intn(vNum)
	b := 0
	for {
		b = rand.Intn(vNum)
		if b != a{break}
	}
	current := cost(graph[a][0],graph[a][1]) + cost(graph[b][0],graph[b][1])
	next := cost(graph[a][0],graph[b][1]) + cost(graph[a][1],graph[b][0])
	if prob(current,next) >= math.Abs(rand.Float64()) && graph[a][0] != graph[b][1] && graph[a][1] != graph[b][0]{
		graph[a][0],graph[b][0] = graph[b][0],graph[a][0]
	}
}
```

##結果

```go
func main()  {
	//最初のコストを見てみる
	sum := 0.0
	for i:=0;i<vNum;i++{
		sum += cost(graph[i][0],graph[i][1])
	}
	fmt.Println("最初のコスト",sum)
	//温度を下げながらループ
	for ;t > 1;t--{
		change()
	}
	sum = 0.0
	for i:=0;i<vNum;i++{
		sum += cost(graph[i][0],graph[i][1])
	}
	//終了後のコスト
	fmt.Println("終了後のコスト",sum)
	//終了後の組み合わせを見てみる
	fmt.Println("組み合わせ",graph)
}
```

最初のコスト 28546.149515622925
終了後のコスト 13487.002687185388
組み合わせ [[42 1] [33 2] [4 3] [44 4] [11 5] [40 6] [23 7] [31 8] [26 9] [49 10] [48 11] [30 12] [43 13] [7 14] [39 15] [20 16] [18 17] [22 18] [45 19] [3 20] [19 21] [36 22] [27 23] [8 24] [0 25] [16
26] [1 27] [24 28] [9 29] [32 30] [28 31] [15 32] [12 33] [10 34] [46 35] [37 36] [14 37] [5 38] [25 39] [6 40] [38 41] [47 42] [41 43] [13 44] [34 45] [17 46] [29 47] [2 48] [35 49] [21 0]]

コストを減らすことができました。大域最適解かどうかはわかりませんがとりあえずは成功したということで完

#さいごに
このシミュレーテッドアニーリングはボルツマンマシンというものに使われてたりします。ニューラルネットワークの欠点は勾配法を使っているので局所最適解に陥ることがあることです。ですので、ニューラルネットワークの中に今回取り上げたアルゴリズムを導入してあげることで局所最適解へ陥るのを抑制してやろうというものです。詳しい理論はまだよくわかってないので適当なこと言ってるかもしれません。申し訳ナス。


#追記
温度をめちゃめちゃゆっくり下げるようにしてみました。

```go
func main()  {
	sum := 0.0
	for i:=0;i<vNum;i++{
		sum += cost(graph[i][0],graph[i][1])
	}

	fmt.Println("最初のコスト",sum)
	n := 0
	for ;t > 1;{
		change()
		//カウンタnが10000で割り切れるときに温度を下げる
		if n % 10000 == 0{
			t --
			n = 0
		}
		n ++
	}

	sum = 0.0
	for i:=0;i<vNum;i++{
		sum += cost(graph[i][0],graph[i][1])
	}


	//終了後のコスト
	fmt.Println("終了後のコスト",sum)
	//終了後の組み合わせを見てみる
	fmt.Println("組み合わせ",graph)
}
```

温度の変化がゆっくりになったので初期の温度を100度にしてやりました。

結果
最初のコスト 28546.149515622925
終了後のコスト 4400.156529118279
組み合わせ [[23 1] [48 2] [39 3] [46 4] [38 5] [44 6] [14 7] [31 8] [24 9] [20 10] [17 11] [33 12] [15 13] [7 14] [13 15] [21 16] [30 17] [22
18] [37 19] [34 20] [16 21] [18 22] [47 23] [9 24] [45 25] [29 26] [28 27] [27 28] [26 29] [11 30] [8 31] [43 32] [40 33] [10 34] [49 35] [42
36] [0 37] [5 38] [3 39] [12 40] [32 41] [36 42] [41 43] [6 44] [25 45] [4 46] [1 47] [2 48] [35 49] [19 0]]

最初のやつより大幅にコストが減少しました。いろんな組み合わせでやってもこれ以上下がらなかったのでもしかしたらこれが最適解なのかもしれません。2000×2000の座標で終了コストが4400なので結構限界ギリギリまでコストを削減してる気がします。