Boost.Graphを使って巡回セールスマン問題を解いてみる

プログラム 数学

この記事はC++ Advent Calendar 2012の10日目の記事です。

はじめに

Boost.Graphでは様々なグラフの問題を扱うことができます。
例えば、グラフの経路探索を抽象化した、幅優先探索や深さ優先探索、最短経路問題を効率よく解くためのA*アルゴリズムや、平面グラフに対する彩色問題を解くためのアルゴリズムも提供されています。

さて、これらのコンテンツを眺めていると、面白そうなテーマとして、TSPやSmall Worldといった用語を見つけることができます。
この中から、今回はTSP(Traveling Salesman Problems, 巡回セールスマン問題)についてのBoost.Graph内での扱いを調べてみました。
それに付随して、その中で使われているBoost.Graphの機能とか考え方についてを簡単に述べます。

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

巡回セールスマン問題（じゅんかいセールスマンもんだい、Traveling Salesman Problem、TSP）は、都市の集合と各2都市間の移動コスト（たとえば距離）が与えられたとき、全ての都市をちょうど一度ずつ巡り出発地に戻る巡回路の総移動コストが最小のものを求める（セールスマンが所定の複数の都市を1回だけ巡回する場合の最短経路を求める）組合せ最適化問題である。 - Wikipedia, 巡回セールスマン問題

問題としては非常に有名な問題で、厳密解を得る場合はNP困難*1な問題であることが知られています。
これらの「困難な」問題に対するアプローチとしてはいくつかあります。
例えば、遺伝的アルゴリズムを用いて探索を行い、厳密解を求める事を放棄しつつもより良い解を求めていく方法や、解の質をある程度保証しつつも、より現実的な時間で何らかの経路を導きだす方法があります。
もちろん、厳密解を出す事も理論上不可能ではありませんが、都市数が増えた場合はNPの壁に引っかかるため、P v.s. NP問題がP = NPとして証明されない限り*2は、「どのような場合においても*3その厳密解を得る」ことは非常に難しいと考えられています。

metric_tsp_approx

この関数がBoost.Graph内でTSPを解くために提供されている関数です。
この関数は巡回セールスマン問題の解（候補）である経路を返します。
関数の説明としては、以下の内容が与えられています。

This is a traveling salesperson heuristic for generating a tour of vertices for a fully connected undirected graph with weighted edges that obey the triangle inequality. The current implementation is based off of the algorithm presented in Introduction to Algorithms on page 1029. This implementation generates solutions that are twice as long as the optimal tour in the worst case. The pseudo-code is listed below. http://www.boost.org/doc/libs/1_52_0/libs/graph/doc/metric_tsp_approx.html

上記の説明から、Boost.Graphを使って巡回セールスマン問題を解くために提供されている関数は以下の特性を持っている事が分かります。

  • 解を導けるグラフに以下の制約がある
    • 無向グラフ(undirected graph)であり、完全グラフ(fully connected graph)であること。
    • 辺の重み（距離）の性質として、三角不等式(triangle inequality)*4を満たす事。
  • この実装によって得られる解は厳密解ではない（かもしれない）が、どのような場合でも厳密解の２倍以下の値であること。（近似率２であること）

関数の名前にapprox(近似の)という名前が入っていることからも本当は分かるのですが、上記の制約を満たす問題・解で満足できるのであればBoost.Graphを使って巡回セールスマン問題の（近似）解を得る事ができます。

また、アルゴリズムの実装としてアルゴリズムイントロダクションのp.1029にある実装を元にしている、とマニュアルに記述があるのですが、記事の執筆段階では残念ながら私は所有していないため、実際に書籍にどのような記述があるかは分かりません。
ただ、マニュアルの同じページにアルゴリズムの概略が示してあります。

TSP-APPROX(G, w)
　select a vertex r in V[G]
　compute a minimum spanning tree T for G from root r
　　using PRIM-MST(G, r, w)
　let L be the list of vertices visited in a preorder traversal of T
　return (the Hamiltonian cycle H that visites the vertices in the order L)

http://www.boost.org/doc/libs/1_52_0/libs/graph/doc/metric_tsp_approx.html

簡単に上記の概略を読み取ると、以下のようになります。

TSP-APPROX は 対象となるグラフGと辺の重み（距離情報）wを受け取る
　頂点rをグラフGの頂点から1つ選択する
　rをrootとするGの最小全域木Tをプリム法*5によって求める
　Lを木Tを先降順*6で訪問した順の頂点リストとする
　上記で求めたLの中から、ハミルトン経路Hを返す

以上の方法で近似率２となる学術的理由については、上記の本がないため分かりませんが、やっていること自体は単純なようです。

マニュアルを読み進めていくと、以下の記述が見つかります。

IN: const Graph& g
　An undirected graph. The type Graph must be a model of Vertex List Graph and Incidence Graph [2].

この事から、TSPを解くためのグラフに無向グラフかつ、何らかのモデルが必要となることが分かります。
このモデル、ならびにGraphのアルゴリズムに対して、どのように解法を追加するかを少し見ていきます。

Graph Concept

グラフのデータ構造には様々なものがあります。
例えば、pair, set > >の様に頂点集合と辺集合を直接持つことで表現する事もできますし、行列を使った表現方法もあります。
これらの表現方法は、解きたい問題によって有利不利があります。
例えば、「この頂点から出ている辺の本数」を求める場合なら、前者の例では辺集合を全て探す必要がありますが、後者であれば該当する辺の行を１行走査するだけですみます。(もちろん、工夫によってある程度は改善できますが、ここでは単純なケースのみを考えます)
ならば全ての機能を備えておけばいいように思いますが、そうした場合はメモリ容量を食いますし、グラフの変更時などに時間がかかりますので、結局のところトレードオフとなるのです。

そこで、Boost.Graphではデータの特徴に対してコンセプトが定義されています。
その中の１つ、VertexListGraphでは「グラフ中の全ての頂点を効率よく走査できること」を示す構造に対して設定されるコンセプトです。
アルゴリズムを実施する際に上記の性質が必要であるなら、これを要求できます。
事実、metrix_tsp_approxのインタフェースとして、VertexListGraphが要求されていることを確認できます。

template
void metric_tsp_approx_from_vertex(
　　const VertexListGraph& g,
　　typename graph_traits::vertex_descriptor start,
　　WeightMap weightmap,
　　VertexIndexMap indexmap,
　　TSPVertexVisitor vis)
　
http://www.boost.org/doc/libs/1_52_0/libs/graph/doc/metric_tsp_approx.html

VertexListGraphのモデル（具体的な実装）としては、Boost.Graphでおなじみのadjacency_listがあります。
しかし、metric_tsp_approxはVertexListGraphとIncidenceGraphの両方の特性を持つグラフを要求しています。
ここのところは自信がないのですが、続く注釈[2]と実装のインタフェースを鑑みるに、実際はVertexListGraphのグラフを渡してやればよいのですが、adjacency_listを渡した場合で特定の場合に処理効率が悪くなるため、別途IncidenceGraphの性質を要求しているのではないのかと思います。
なお、この両方の特性を満たすモデルとしては、adjacency_matrix(隣接行列)がありますので、今回の実装ではこれを利用します。

visitor

Boost.Graphにはvisitorという概念が存在します。
これは、STLでいうfunctor(関数オブジェクト)のようなものです。
私たちは、問題を解くためにグラフを作り、既存のアルゴリズムの中に注入しますが、グラフの問題を解いているいずれかのタイミングで何らかの処理をしたくなることがあります。
このように、問題の決められたタイミングでコールバックを行ってくれる関数オブジェクトを登録することができ、問題を拡張することができるようになっています。
そのため、functorは大体１つの呼び出し(operator())しか持ちませんが、visitorは複数の関数を持つ場合があります。

ここでは一例として、BFS Visitorを挙げます。
このvisitorは幅優先探索(Breadth First Search, BFS)を実行する関数breadth_first_searchのVisitorとして提供されています。
詳しくはマニュアルを見て欲しいのですが、「処理開始時」「アルゴリズムの実行中、まだ未探索のノードを見つけたとき」などなどの様々なタイミングで処理ができるように関数が定義されています。

TSP Tour Visitor

では、上記のmetric_tsp_approx関数にもVisitorは存在するのでしょうか？
それはあり、TSP Tour Visitorというものが定義されています。
が、このVisitorにはvisit_vertexという関数が１つだけしか定義されていません。
実装を見てみたところ、上記のアルゴリズムでハミルトン経路が求められた後、その経路を順に呼び出しているだけ見たいです。( 一番下の関数metric_tsp_approx_from_vertexの実装内を参照)

なので、上記Visitorのコンセプトを満たしているものとしてtsp_visitorが定義されていますが、専ら経路・距離の出力用に使用されています。
解を求めた後に順番に呼び出しているので、この場合のVisitorを活用して問題を拡張するのは難しいように思います。

実装

これまでの説明をふまえて、実際にTSPアルゴリズムを使ってみました。
ソースコードは以下の通りです。
Githubにも問題とあわせてアップロードしてあるので、そちらもご参照ください。
# サンプルなのでusing namespaceの使い方が酷いですが、目をつぶってください。

#include <iostream>
#include <vector>
#include <fstream>
#include <set>

#include <boost/foreach.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/integer_traits.hpp>
#include <boost/graph/adjacency_matrix.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/simple_point.hpp>
#include <boost/graph/metric_tsp_approx.hpp>

using namespace boost;
using namespace std;

// simple_pointはboost/graph/simple_point.hppで定義されるx,yでアクセスできるペア構造体
typedef vector<simple_point<double> > PosVector;
typedef adjacency_matrix<undirectedS, no_property,
	property <edge_weight_t, double> > Graph;
typedef graph_traits<Graph>::vertex_descriptor Vertex;
typedef graph_traits<Graph>::edge_descriptor Edge;
typedef vector<Vertex> Container;
typedef property_map<Graph, edge_weight_t>::type WeightMap;
typedef property_map<Graph, vertex_index_t>::type VertexMap;

void usage(char* argv[]) {
	cout << "usage: " << argv[0] << " input.txt" << endl
		 << "input.txt has point on line. for example" << endl
		 << "10 30" << endl;
	return;
}

void print_graph(Graph& g, PosVector& pos) {

	WeightMap wmap(get(edge_weight, g));
	VertexMap vmap(get(vertex_index, g));

	cout << "All Vertices: ";
	BOOST_FOREACH(auto v, vertices(g)) {
		auto p = pos[vmap[v]];
		cout << "[" << vmap[v] <<  ":(" << p.x << "," << p.y << ")] ";
	}
	cout << endl << endl;

	cout << "All Edges: ";
	BOOST_FOREACH(auto e, edges(g)) {
		cout << "(" << vmap[target(e, g)] << "-" << vmap[source(e, g)] 
			<< ": " << wmap[e] << ") ";
	}
	cout << endl;

}

// ファイル入力 -> 座標変換
PosVector get_input(char* filename) {
	ifstream fin(filename);
	if (!fin) {
		cout << "ifstream error." << endl;
		exit(1);
	}

	PosVector vs;
	string line;
	// 取得方法はサンプルをほぼ流用
	while (getline(fin, line)) {
		simple_point<double> vertex;

		size_t idx(line.find(","));
		string xStr(line.substr(0, idx));
		string yStr(line.substr(idx + 1, line.size() - idx));

		vertex.x = lexical_cast<double>(xStr);
		vertex.y = lexical_cast<double>(yStr);
		vs.push_back(vertex); 
	}
	fin.close();
	return vs;
}

// 座標 -> 完全グラフ化（距離込)
Graph make_graph(PosVector vs) {
	// 隣接行列なので、頂点数のみでエッジなしのグラフが構成される
	Graph g(vs.size());
	WeightMap wmap(get(edge_weight, g));
	VertexMap vmap(get(vertex_index, g));

	// pair<頂点iter, 頂点iter>
	auto verts = vertices(g);
	for (auto src = verts.first; src != verts.second; src++) {
		for (auto dest = src; dest != verts.second; dest++) {
			// 自己ループのエッジはなし
			if (dest == src) continue;
			double dx = static_cast<double>(vs[vmap[*src]].x - vs[vmap[*dest]].x);
			double dy = static_cast<double>(vs[vmap[*src]].y - vs[vmap[*dest]].y);
			double weight = sqrt(pow(dx, 2.0) + pow(dy, 2.0));

			// 重み付き辺の作成
			Edge e;
			bool inserted;
			boost::tie(e, inserted) = add_edge(*src, *dest, g);
			wmap[e] = weight;
		}
	}
	return g;
}

void my_tsp_tour(Graph& g, PosVector& pos) {
	vector<Vertex> res; 
	VertexMap vmap(get(vertex_index, g));
	WeightMap wmap(get(edge_weight, g));

	// PATH (ONLY)
	// metric_tsp_approx_tour(g, back_inserter(res)); 

	// PATH with distance
	double dist = 0.0;
	res.clear();
	metric_tsp_approx(g, make_tsp_tour_len_visitor(
		g, back_inserter(res), dist, wmap));
	// result output
	cout << "PATH DISTANCE: " << dist << endl;
	cout << "PATH: ";
	BOOST_FOREACH(auto v, res) {
		cout << vmap[v] << " -> ";
	}
	cout << endl << endl;

}

int main(int argc, char* argv[]) {
	// input check
	if(argc < 2) {
		usage(argv);
		return 0;
	}

	auto in = get_input(argv[1]);
	Graph g = make_graph(in);
	print_graph(g, in);
	my_tsp_tour(g, in);

	return 0;
}

少々長くなっていますが、このアルゴリズムを用いる場合のポイントはmainに書かれている通りで、

さえ守ればOKです。

関数別に見ていきます

  • get_inputでファイルの入力内容を(x,y)のペアのvectorに変換しています。
  • make_graphでグラフを作成しています。
    • 今回のグラフは隣接行列なので、グラフのインスタンス化のためにはノード数を指定することで、そのノード数を持つ辺０のグラフができます。
    • vertices(g)が返すのはstd::pair<頂点.begin(), 頂点.end()>です。今回は自身をのぞく他の頂点と無向辺を結びます。
    • 三角不等式を守るために入力頂点を座標と見なして距離を算出します(これも頂点のプロパティマップとする方が良かったと思いますが、今回はしていないです）。
    • 辺のプロパティマップに２頂点間の距離を辺の重みとして設定します。
  • my_tsp_tourでmetric_tsp_approxを実施しています。
    • len_visitorを使用することで、経路とその経路の距離がまとめて取れます
    • 今回はグラフが作成できることを前提としているので、そのエラーチェックなどは行っていません。

さて、では実際に入力を行ってみます。

Run

１つ目は簡易なテストデータを用意しました。
下記の文字イメージが示すような格子イメージの頂点です。

座標: (0, 0), (2, 0), (4, 0), (0, 2), (2, 2), (4, 2), (0, 4), (2, 4), (4, 4)
ファイル内容:
0,0
2,0
4,0
0,2
2,2
4,2
0,4
2,4
4,4

All Vertices: [0:(0,0)] [1:(2,0)] [2:(4,0)] [3:(0,2)] [4:(2,2)] [5:(4,2)] [6:(0,4)] [7:(2,4)] [8:(4,4)]

All Edges: (0-1: 2) (0-2: 4) (1-2: 2) (0-3: 2) (1-3: 2.82843) (2-3: 4.47214) (0-4: 2.82843) (1-4: 2) (2-4: 2.82843) (3-4: 2) (0-5: 4.47214) (1-5: 2.82843) (2-5: 2) (3-5: 4) (4-5: 2) (0-6: 4) (1-6: 4.47214) (2-6: 5.65685) (3-6: 2) (4-6: 2.82843) (5-6: 4.47214) (0-7: 4.47214) (1-7: 4) (2-7: 4.47214) (3-7: 2.82843) (4-7: 2) (5-7: 2.82843) (6-7: 2) (0-8: 5.65685) (1-8: 4.47214) (2-8: 4) (3-8: 4.47214) (4-8: 2.82843) (5-8: 2) (6-8: 4) (7-8: 2)
PATH DISTANCE: 21.3006
PATH: 0 -> 1 -> 2 -> 5 -> 8 -> 4 -> 3 -> 6 -> 7 -> 0 ->

ということで、一応ハミルトン閉路の出力は行われています。
が、実際の最短パス（例）は [0 -> 1 -> 2 -> 5 -> 8 -> 4 -> 7 -> 6 -> 3 -> 0 (DIST: 16 + 2√2 < 20)] であるため、厳密解は得られていないことが分かります。
まあ、これは解の求め方の性質上だとは思うのですが、、、


TSPLIB

TSPは広く研究されており、作成したアルゴリズムの性能比較のための問題セットが用意されています。
http://elib.zib.de/pub/Packages/mp-testdata/tsp/tsplib/tsplib.html （など）
このページのトップにある画像は面白い問題で、最短経路を導けると「TSP」の文字が浮かび上がってくるようになっています。
この問題はtsp225と呼ばれる問題です。
実際にこの問題を解いてみました。ただし、問題入力の簡単化のため、端数は全て切り捨てて計算しているため、厳密に同じ問題を解いているわけではありません。
結果は次の通りでした。

PATH DISTANCE: 5244.07
PATH: 0 -> 2 -> 197 -> 3 -> 4 -> 5 -> 6 -> 7 -> 8 -> 42 -> 41 -> 43 -> 45 -> 193 -> 194 -> 217 -> 44 -> 47 -> 192 -> 195 -> 191 -> 190 -> 198 -> 132 -> 189 -> 224 -> 46 -> 1 -> 58 -> 57 -> 206 -> 48 -> 50 -> 49 -> 56 -> 55 -> 54 -> 51 -> 52 -> 53 -> 188 -> 26 -> 187 -> 116 -> 115 -> 222 -> 114 -> 113 -> 112 -> 111 -> 63 -> 62 -> 61 -> 60 -> 59 -> 64 -> 65 -> 66 -> 67 -> 68 -> 69 -> 70 -> 71 -> 72 -> 73 -> 74 -> 75 -> 215 -> 218 -> 216 -> 76 -> 29 -> 28 -> 27 -> 203 -> 25 -> 24 -> 207 -> 23 -> 22 -> 12 -> 11 -> 10 -> 9 -> 13 -> 14 -> 15 -> 16 -> 19 -> 20 -> 202 -> 18 -> 17 -> 21 -> 35 -> 36 -> 37 -> 38 -> 39 -> 40 -> 34 -> 32 -> 33 -> 201 -> 205 -> 30 -> 31 -> 77 -> 78 -> 96 -> 95 -> 94 -> 79 -> 80 -> 208 -> 93 -> 92 -> 91 -> 90 -> 89 -> 86 -> 209 -> 83 -> 82 -> 81 -> 84 -> 85 -> 130 -> 210 -> 129 -> 221 -> 128 -> 127 -> 126 -> 125 -> 124 -> 131 -> 88 -> 87 -> 220 -> 97 -> 98 -> 99 -> 100 -> 101 -> 102 -> 103 -> 219 -> 104 -> 105 -> 106 -> 107 -> 108 -> 110 -> 109 -> 117 -> 118 -> 184 -> 119 -> 174 -> 120 -> 121 -> 122 -> 123 -> 170 -> 169 -> 168 -> 167 -> 211 -> 149 -> 148 -> 147 -> 146 -> 145 -> 144 -> 152 -> 153 -> 154 -> 155 -> 156 -> 143 -> 142 -> 200 -> 151 -> 150 -> 213 -> 166 -> 165 -> 164 -> 163 -> 212 -> 157 -> 162 -> 136 -> 135 -> 182 -> 134 -> 137 -> 138 -> 139 -> 140 -> 141 -> 161 -> 160 -> 159 -> 158 -> 214 -> 133 -> 171 -> 172 -> 180 -> 173 -> 179 -> 178 -> 175 -> 176 -> 177 -> 181 -> 183 -> 185 -> 186 -> 223 -> 204 -> 196 -> 199 -> 0 ->

この問題の最適解は距離3919であるため、誤差を踏まえても最適解の２倍以下の距離の経路を解として出力している事が分かります。

まとめ

Boost.Graphに用意されているTSPの仕組みを使って、限定された形ではありますが巡回セールスマン問題の近似解を求めることがかなり簡単にできることが分かりました。
Graphは楽しいので、機会があれば別の機能も使ってみたいと思います。
また、間違いなどありましたらご指摘いただければ幸いです。

*1:NP完全問題よりも難しい問題。NP完全問題の例としては充足可能性問題(SAT)などがある。NP完全問題は解がYesかNoで答えられる決定問題が属するクラスであり、ある問題が決定問題より難しく、特殊例や部分問題がNP完全な問題である事を示せば、その問題はNP困難足りうる。巡回セールスマン問題に大しては、NP完全な問題であるハミルトン閉路問題を部分問題として含み、さらにそのコストが最小になるパターンを求めなくてはならないのでNP困難である。

*2:証明こそまだだが、P ≠ NPであるという主張の方が強い

*3:勘違いされやすいですが、ノード数が増え過ぎるような「任意の」問題に対して難しいのであり、状況を限定する場合は現実的な時間で厳密解を得られる可能性があります

*4:ユークリッド空間上での三角形の辺の長さに成り立つ関係。三角形の辺a,b,cに対して、長さがa

*5:http://ja.wikipedia.org/wiki/%E3%83%97%E3%83%AA%E3%83%A0%E6%B3%95

*6:ノードの子供より先に自分自身を評価して、次に子供を評価する方式。深さ優先・幅優先については記載がないため、ここの詳細はアルゴリズムイントロダクションにて確認したい