📍

巡回セールスマン問題のいくつかの解法を試してみた

2024/03/24に公開

Golang

巡回セールスマン問題を勉強したので、教科書に載っているいくつかの解法を実際に試してみました。
言語はGoで書いています。

教科書は次のものを使っています。

AI時代の離散数学

巡回セールスマン問題とは

ネットワーク[G;d]の全ての点において, すべての点を一度ずつ訪問して元に戻る巡回路の中で最短のものを求める

問題です。今回はとくに、平面上の巡回セールスマン問題と呼ばれる、平面上にn個の点が存在し、それらの二点の距離がユークリッド距離で定められているような場合を考えます。つまり、Gは完全グラフです。

アルゴリズム

巡回セールスマン問題は一般のグラフに対しては、NP困難です。(このあたりはよくわかってないです...)
ですので、近似的なアルゴリズムでよりよい解を求めることになります。
教科書に載っている近似アルゴリズムは次の通りです。

ランダム法(random method)
最近近傍法(nearest neighbor method)
挿入法(insertion method)
- 最近挿入法(nearest insertion method)
- 最遠挿入法(farthest insertion method)
- 最廉挿入法(cheapest insertion method)
(局所探索法(local search)による改善)

この記事では局所探索法による改善は行わないことにします(そこが一番大切な気もしますが)。

グラフの頂点数 $N$ について、 $N \leq 15$ 程度の場合は、実際に動的計画法を使って解くことができるので、その場合に最適解にどれくらい近い解が得られるかを調べてみようと思います。

問題例

競技プログラミングの鉄則の中のB23 - Traveling Salesman Problemを解くことを考えます。とくに、テストケースのin14.txtをためします。

input.txt

想定される出力は

output.txt

3441.451976970622

です。( $10^{-5}$ 以下の誤差は許容されます)

実装

簡単な説明とコード、計算量を記載していきます。より高速化できる実装法があれば教えてください。
基本的にmain関数の中身だけを記載します。入出力のために色々追加しているので詳しくはテンプレートをご確認ください。

AtCoder用テンプレート

package main

import (
	"bufio"
	"fmt"
	"os"
	"strconv"
	"strings"
)

func main() {
}

func sum[T int | float64](arr ...T) T {
	sum := arr[0]
	for i := 1; i < len(arr); i++ {
		sum += arr[i]
	}
	return sum
}

func pow(x, n int) int {
	res := 1
	for n > 0 {
		if n&1 == 1 {
			res = res * x
		}
		x = x * x
		n >>= 1
	}
	return res
}

func lcm(a, b int) int {
	return a / gcd(a, b) * b
}

func gcd(a, b int) int {
	if b == 0 {
		return a
	}
	return gcd(b, a%b)
}

func abs[T int | float64](x T) T {
	if x < 0 {
		return -x
	}
	return x
}

func min[T int | float64](arr ...T) T {
	min := arr[0]
	for _, v := range arr {
		if v < min {
			min = v
		}
	}
	return min
}

func max[T int | float64](arr ...T) T {
	max := arr[0]
	for _, v := range arr {
		if v > max {
			max = v
		}
	}
	return max
}

func ceilDiv(a, b int) int {
	if b < 0 {
		return ceilDiv(-a, -b)
	}
	if a < 0 {
		return a / b
	}
	return (a + b - 1) / b
}

func floorDiv(a, b int) int {
	if b < 0 {
		return floorDiv(-a, -b)
	}
	return (a - positiveMod(a, b)) / b
}

func positiveMod(a, b int) int {
	return (a%b + b) % b
}

var _stdInReader = bufio.NewReader(os.Stdin)

func readLine() string {
	line, _ := _stdInReader.ReadString('\n')
	return strings.TrimSpace(line)
}

func readString() string {
	return readLine()
}

func readInt() int {
	strs := readString()
	num, _ := strconv.Atoi(strs)
	return num
}

func readLong() int64 {
	strs := readString()
	num, _ := strconv.ParseInt(strs, 10, 64)
	return num
}

func readStrings() []string {
	line := readLine()
	return strings.Split(line, " ")
}

func readInts() []int {
	strs := readStrings()
	arr := make([]int, len(strs))
	for i := 0; i < len(strs); i++ {
		arr[i], _ = strconv.Atoi(strs[i])
	}
	return arr
}

func readLongs() []int64 {
	strs := readStrings()
	arr := make([]int64, len(strs))
	for i := 0; i < len(strs); i++ {
		arr[i], _ = strconv.ParseInt(strs[i], 10, 64)
	}
	return arr
}

func readDoubles() []float64 {
	strs := readStrings()
	arr := make([]float64, len(strs))
	for i := 0; i < len(strs); i++ {
		arr[i], _ = strconv.ParseFloat(strs[i], 64)
	}
	return arr
}

func scanStringVariables(args ...*string) {
	n := len(args)
	scanned := 0
	for n > scanned {
		inputSlice := readStrings()
		m := len(inputSlice)
		if m == 0 || m > n-scanned {
			fmt.Print("Something wrong in scanStringVariables !!!")
			return
		}
		for i := 0; i < m; i++ {
			*args[i+scanned] = inputSlice[i]
		}
		scanned += m
	}
}

func scanIntVariables(args ...*int) {
	n := len(args)
	scanned := 0
	for n > scanned {
		inputSlice := readInts()
		m := len(inputSlice)
		if m == 0 || m > n-scanned {
			fmt.Printf("m %d n %d scanned %d. ", m, n, scanned)
			fmt.Print("Something wrong in scanIntVariables !!!")
			return
		}
		for i := 0; i < m; i++ {
			*args[i+scanned] = inputSlice[i]
		}
		scanned += m
	}
}

func scanLongVariables(args ...*int64) {
	n := len(args)
	scanned := 0
	for n > scanned {
		inputSlice := readLongs()
		m := len(inputSlice)
		if m == 0 || m > n-scanned {
			fmt.Print("Something wrong in scanIntVariables !!!")
			return
		}
		for i := 0; i < m; i++ {
			*args[i+scanned] = inputSlice[i]
		}
		scanned += m
	}
}

func write(arg ...interface{})                 { fmt.Print(arg...) }
func writeLine(arg ...interface{})             { fmt.Println(arg...) }
func writeFormat(f string, arg ...interface{}) { fmt.Printf(f, arg...) }

DPによる解法(最適解)

const[n][i]: すでに訪れた都市の集合が $n$ であって、最後にいる都市が $i$ であるときの、合計コストの最小値。

と定義することで、DPで解くことができます。 $n$ は、集合を二進数とみることで数値として扱えます。

DPのコード

func main() {
	n := readInt()
	XY := make([][2]int, n)
	for i := 0; i < n; i++ {
		scanIntVariables(&XY[i][0], &XY[i][1])
	}
	D := make([][]float64, n)
	for i := 0; i < n; i++ {
		D[i] = make([]float64, n)
		for j := 0; j < n; j++ {
			D[i][j] = math.Sqrt(float64((XY[i][0]-XY[j][0])*(XY[i][0]-XY[j][0]) + (XY[i][1]-XY[j][1])*(XY[i][1]-XY[j][1])))
		}
	}
	inf := float64(1 << 60)
	cost := make([][]float64, 1<<n)
	for i := 0; i < 1<<n; i++ {
		cost[i] = make([]float64, n)
		for j := 0; j < n; j++ {
			cost[i][j] = inf
		}
	}
	cost[0][0] = 0
	for b := 0; b < 1<<n; b++ {
		for i := 0; i < n; i++ {
			for j := 0; j < n; j++ {
				if i == j { // 同じ頂点には移動しない
					continue
				}
				if (b>>j)&1 == 1 { // すでに j に訪問済み
					continue
				}
				s := b | 1<<j // j を訪問済みにする
				cost[s][j] = min(cost[s][j], cost[b][i]+D[i][j])
			}
		}
	}
	writeLine(cost[(1<<n)-1][0])
}

計算量は、 $\text{O}(2^NN^2)$ です。

ランダム法

ランダムに頂点の順序を決め、その順序通りに巡回したときの距離を求めるという操作を有限時間繰り返し、その中で最小のものを近似解として出力します。

ここでは、制限時間は1980msとします。

ランダム法のコード

func main() {
	n := readInt()
	XY := make([][2]int, n)
	for i := 0; i < n; i++ {
		scanIntVariables(&XY[i][0], &XY[i][1])
	}
	D := make([][]float64, n)
	for i := 0; i < n; i++ {
		D[i] = make([]float64, n)
		for j := 0; j < n; j++ {
			D[i][j] = math.Sqrt(float64((XY[i][0]-XY[j][0])*(XY[i][0]-XY[j][0]) + (XY[i][1]-XY[j][1])*(XY[i][1]-XY[j][1])))
		}
	}
	arr := make([]int, n)
	for i := 0; i < n; i++ {
		arr[i] = i
	}
	s := time.Now()
	ans := math.Inf(1)
	for {
		rand.Shuffle(n, func(i, j int) {
			arr[i], arr[j] = arr[j], arr[i]
		})
		tmp := 0.0
		for i := 0; i < n; i++ {
			tmp += D[arr[i]][arr[(i+1)%n]]
		}
		ans = min(ans, tmp)
		if time.Since(s).Milliseconds() > 1980 {
			break
		}
	}
	writeLine(ans)
}

計算量は、一回あたり $\text{O}(N)$ です。実際に計測してみると、1980msの間に約5300000回のループが行われていました。

5回試してみた結果は以下の通りです。

3886.654694326012
3857.1945450443095
3855.2302193326327
3798.6987806854304
3825.7162365425474

挿入法

こちらも貪欲法の一種で、部分巡回路 $(v_{j_1}, v_{j_2}, \dots , v_{j_l})$ に、点 $v_k$ を追加していくというものです。 $v_k$ を挿入する場所は、

\Delta_k(i) = d(v_{j_i}, v_k) + d(v_k, v_{j_{i+1}}) - d(v_{j_i}, v_{j_{i+1}})

を最小にするような $v_{j_i}$ の後ろです。また、挿入する $v_k$ の選び方によって

最近挿入法
最遠挿入法
最廉挿入法

があります。ここでは、最近挿入法を試してみようと思います。最近挿入法は、部分巡回路 $(v_{j_1}, v_{j_2}, \dots , v_{j_l})$ に対して、

D_\pi(k) = \underset{1 \leq i \leq l}{\min}{d(v_k, v_{j_i})}

を最小にする $v_k$ を用いるものです。

最近挿入法のコード

type Edge struct {
	from int
	to   int
	cost float64
}

type Node struct {
	value int
	prev  *Node
	next  *Node
}

type UniqueDoublyLinkedList struct {
	head     *Node
	tail     *Node
	accessor map[int]*Node
	len      int
}

func newLinkedList() *UniqueDoublyLinkedList {
	return &UniqueDoublyLinkedList{
		head:     nil,
		tail:     nil,
		accessor: make(map[int]*Node),
	}
}

func (ll *UniqueDoublyLinkedList) add(value int) {
	node := &Node{value: value}
	if ll.head == nil {
		ll.head = node
		ll.tail = node
	} else {
		ll.tail.next = node
		node.prev = ll.tail
		ll.tail = node
	}
	ll.accessor[value] = node
	ll.len++
}

func (ll *UniqueDoublyLinkedList) insertAfter(preValue int, value int) {
	preNode := ll.accessor[preValue]
	node := &Node{value: value}
	node.prev = preNode
	node.next = preNode.next
	if preNode == ll.tail {
		ll.tail = node
	} else {
		preNode.next.prev = node
	}
	preNode.next = node
	ll.len++
	ll.accessor[value] = node
}

func (ll *UniqueDoublyLinkedList) has(value int) bool {
	_, ok := ll.accessor[value]
	return ok
}

func main() {
	n := readInt()
	XY := make([][2]int, n)
	for i := 0; i < n; i++ {
		scanIntVariables(&XY[i][0], &XY[i][1])
	}
	D := make([][]float64, n)
	for i := 0; i < n; i++ {
		D[i] = make([]float64, n)
		for j := 0; j < n; j++ {
			D[i][j] = math.Sqrt(float64((XY[i][0]-XY[j][0])*(XY[i][0]-XY[j][0]) + (XY[i][1]-XY[j][1])*(XY[i][1]-XY[j][1])))
		}
	}
	G := make([][]Edge, n)
	for i := 0; i < n; i++ {
		for j := 0; j < n; j++ {
			G[i] = append(G[i], Edge{i, j, D[i][j]})
		}
	}

	res := 0.0
	nonVisited := make(map[int]bool)
	for i := 0; i < n; i++ {
		nonVisited[i] = true
	}
	pi := newLinkedList() // 部分巡回路
	s := 0
	v := s
	pi.add(v)
	delete(nonVisited, v)
	que := priorityqueue.New[Edge](func(a, b Edge) int {
		if a.cost < b.cost {
			return -1
		} else if a.cost > b.cost {
			return 1
		}
		return 0
	}) // 挿入する辺の候補
	for _, e := range G[v] {
		que.Push(e)
	}
	for !que.Empty() {
		e := que.Pop()
		if pi.has(e.to) {
			continue
		}
		if pi.len == 1 {
			res += e.cost
			pi.add(e.to)
		} else {
			u := pi.head
			minCost := math.Inf(1)
			minV := -1
			for u.next != nil {
				cost := D[u.value][e.to] + D[e.to][u.next.value] - D[u.value][u.next.value]
				if cost < minCost {
					minCost = cost
					minV = u.value
				}
				u = u.next
			}
			res += minCost
			pi.insertAfter(minV, e.to)
		}
		delete(nonVisited, e.to)
		v = e.to
		for u := range nonVisited {
			que.Push(G[v][u])
		}
	}
	res += D[pi.tail.value][s]

	writeLine(res)
}

かなり雑な実装になりました。挿入を行うために連結リストを用いました。すべての頂点の値はユニークなので、挿入を定数時間で行うためにmapで直接ノードにアクセスできるようにしています。また、挿入する $v_k$ を毎回探すと時間がかかるため、優先度付きキューを用いています。

計算量は大きく見積もっても $\text{O}(N^3)$ (始点をランダムにひとつ選ぶ場合は $\text{O}(N^2)$ ) くらいだと思います。

今回の問題例に対する出力は

3441.451976970623

で、最適解を得ることができました。

おまけ( $N=150$ )

$N=15$ 程度だと実際にDPで解けるのであまり面白くないです。せっかくなので現実的な時間では解くことができない場合でどれくらい差が出るのか試してみたいと思います。

入力はランダムに作成した次を用います。

input

input.txt

150
3892 -5182
24 5203
-6981 1752
6011 1659
2037 4086
-4522 3455
825 -4313
3677 175
-1063 -6756
-9948 -4369
-1610 7499
2510 6130
-7403 -2122
-7770 -8959
1134 2583
9670 -635
6735 2212
-8641 8333
-8847 3569
8085 7334
-2093 2507
-8248 -6616
395 -3256
9188 8271
2388 6100
8940 -8988
-5426 -5545
8879 911
-4357 -864
2394 -4992
-1784 1219
3875 4033
6409 3843
2567 6673
-1936 8450
-9515 844
5417 1166
6281 3503
-9788 -4932
7533 -172
-9603 232
5399 3788
1495 6147
-3731 7457
-4988 -500
9199 4551
-8171 -5191
920 166
-5810 2833
4383 7259
7067 -3149
6753 7877
-730 -6274
-3100 -9303
3238 1466
9657 -4907
5177 -5193
8081 -1567
-9577 -3469
8393 -3330
-2914 -9827
-7304 1287
8859 9484
-8028 2896
1962 7024
-7150 -6095
-5960 5027
-4018 -493
-2573 -57
-3145 6131
2326 -4624
219 -9358
6334 -3410
9808 -629
490 3251
-4036 -50
5974 -6034
9533 8221
3568 9945
8690 2291
9444 3286
8356 -872
5213 9676
-6692 9072
3121 5334
4698 -6635
5765 9585
-6990 -614
-2988 5271
-4967 -12
-2604 8409
82 -5961
-1824 2646
-4322 -9678
-2198 -5055
-7389 3061
-4864 5168
-6223 -7519
9712 -9537
-448 468
-3291 5
3527 6021
-5999 80
-1004 5706
5632 5521
5931 -7466
1355 -1063
-429 1219
6093 4189
4374 -7568
8447 -8160
9871 -9711
4913 9974
-9245 -9697
-1781 -4598
7112 6853
-9424 3651
4842 -1225
-7435 6643
3630 3979
9636 2467
6412 -3221
9994 5081
1299 4130
-9370 5638
-9913 9787
7515 -4810
-2434 9842
2397 3112
-3490 5646
-1172 712
-7687 -8665
-7553 6628
2061 -8665
4324 2821
-5327 4737
3603 9989
-4455 -3441
8462 9621
7083 5086
-2537 -2589
-8554 7211
-9789 8131
5786 -7554
1517 -1680
5349 -8529
9014 5967
9313 -5376
-5860 -2507
-3568 9774

結果は以下の通りです。

アルゴリズム	実行時間	出力
ランダム法	1989ms	1334675.015424154
最近近傍法	11ms	222276.1101097709
最近挿入法	668ms	242387.82579042445

まとめ

最近挿入法の実装がかなり無駄が多い気がします。もっとシンプルにしたり、計算量の改善ができる方法があれば教えてください。

$N=150$ のときに最近近傍法が一番良い結果を得ることは意外でした。ある程度散らばりのあるデータの場合は、最後のボトルネックが小さくなるのかもしれません。

教科書を読むだけではイメージし辛いアルゴリズムも、実際に実装することで流れが理解できたので、これからも続けていこうと思います。

最後まで読んでいただきありがとうございました。見苦しいコードが多くてすみません。

巡回セールスマン問題とは

アルゴリズム

問題例

実装

DPによる解法(最適解)

ランダム法

最近近傍法

挿入法

おまけ(N=150)

まとめ

Discussion

おまけ( $N=150$ )