y_uti のブログ

統計、機械学習、自然言語処理などに興味を持つエンジニアの技術ブログです

巡回セールスマン問題の近似解法を眺める

今年の 8 月頃に、paiza のプログラミングの問題に取り組みました。巡回セールスマン問題を解くものです*1
paiza.jp

巡回セールスマン問題の近似解法として、さまざまなバリエーションが知られているようです。そこで今回は、いくつかの簡単な方法を PHP で実装して、近似解を求める様子を確認してみます。アルゴリズムの実装にあたっては以下の資料を参考にしました。前者は Lecture Notes の "Networks 2", "TSP Heuristics" にさまざまなアルゴリズムの説明があります。後者からは講義資料 17 ページの Greedy (Kruskal の最小全域木構成アルゴリズム) を実装しました。

  1. Logistical and Transportation Planning Methods | Civil and Environmental Engineering | MIT OpenCourseWare
  2. 情報数理工学演習第二B (東京大学での 2009 年の講義のようです)
Nearest Neighbour 法

Nearest Neighbour 法は、未訪問の都市のなかで、現在の都市から一番近いものを選んで移動する方法です。いわゆる貪欲法として最も単純な解き方かもしれません。計算部分の実装は次のとおりです。$data は平面上の座標点の集合を保持する二次元配列で、これを訪問する順番に並べ替えて返します*2

<?php
function solve($data) {
    array_unshift($data, [0, 0]);
    $n = count($data);
    for ($i = 1; $i < $n; ++$i) {
        $nearest = find_nearest($data[$i - 1], $data, $i, $n);
        swap($data, $i, $nearest);
    }
    array_shift($data);
    return $data;
}

function find_nearest($curr, $data, $beg, $end) {
    $min_index = -1;
    $min_distance = INF;
    for ($i = $beg; $i < $end; ++$i) {
        if (($d = distance($curr, $data[$i])) < $min_distance) {
            $min_index = $i;
            $min_distance = $d;
        }
    }
    return $min_index;
}

function distance($a, $b) {
    return sqrt(($b[0] - $a[0]) ** 2 + ($b[1] - $a[1]) ** 2);
}

function swap(&$array, $i, $j) {
    [$array[$i], $array[$j]] = [$array[$j], $array[$i]];
}

ランダムに 100 点のデータを作成して、Nearest Neighbour 法で解を求める様子をアニメーションにしてみます。PHP のプログラムにデバッグコードを追加して途中の状態を出力させ、それを Python で読み込んで matplotlib で描画しています*3
f:id:y_uti:20171104101740g:plain

Insertion 法

Insertion 法は、未訪問の都市を一つ選んで、現在の訪問ルートのなかで距離の増加が最も少ない箇所に挿入することで訪問ルートを更新していく方法です。未訪問の都市をどのように選択するかによって、さらにいくつかのバリエーションがあります。最後の Cheapest Insertion は他の方法よりも計算量のオーダーが大きくなります。

アルゴリズム 選択方法
Random Insertion 特に基準を設けず適当に選択する
Nearest Insertion 現在のルートに含まれる都市に最も近いものを選択する
Farthest Insertion 現在のルートに含まれる都市から最も遠いものを選択する
Cheapest Insertion すべての候補のなかで距離の増加が最も少ないものを選択する

Random Insertion 法での計算方法は以下のとおりです。各ステップで未訪問の都市を無作為に選択するので、ここでは単純に $data の順番に選択しています。

<?php
function solve($data) {
    array_unshift($data, [0, 0]);
    $answer = [$data[0]];
    $n = count($data);
    for ($i = 1; $i < $n; ++$i) {
        $pos = find_insert_position($data[$i], $answer);
        array_splice($answer, $pos, 0, [$data[$i]]);
    }
    array_shift($answer);
    return $answer;
}

function find_insert_position($node, $data) {
    $n = count($data);
    $min_cost = distance($data[$n - 1], $node);
    $pos = $n;
    for ($i = 1; $i < $n; ++$i) {
        $c = distance($data[$i - 1], $node) + distance($node, $data[$i]) - distance($data[$i - 1], $data[$i]);
        if ($c < $min_cost) {
            $min_cost = $c;
            $pos = $i;
        }
    }
    return $pos;
}

先ほどと同様に、解を求める様子をアニメーションにしてみます。Nearest Neighbour 法と異なり、Insertion 法ではルートの一部が更新されていく形で計算が進む様子がわかります。
f:id:y_uti:20171104134800g:plain

Nearest Insertion 法での計算方法は以下のとおりです。select 関数で都市を選択します。solve 関数の省略部分や find_insert_position 関数は、Random Insertion 法の実装と同じです。

<?php
function solve($data) {
    ...
    for ($i = 1; $i < $n; ++$i) {
        $next = select($data, $i);
        $pos = find_insert_position($data[$next], $answer);
        array_splice($answer, $pos, 0, [$data[$next]]);
        swap($data, $i, $next);
    }
    ...
}

function select($data, $k) {
    $result = -1;
    $min = INF;
    $n = count($data);
    for ($i = $k; $i < $n; ++$i) {
        for ($j = 0; $j < $k; ++$j) {
            if (($d = distance($data[$i], $data[$j])) < $min) {
                $result = $i;
                $min = $d;
            }
        }
    }
    return $result;
}

Farthest Insertion 法は、Nearest Insertion 法の select 関数を置き換えるだけです。

<?php
function select($data, $k) {
    $result = -1;
    $max = 0;
    $n = count($data);
    for ($i = $k; $i < $n; ++$i) {
        $min = INF;
        for ($j = 0; $j < $k; ++$j) {
            if (($d = distance($data[$i], $data[$j])) < $min) {
                $min = $d;
            }
        }
        if ($min > $max) {
            $result = $i;
            $max = $min;
        }
    }
    return $result;
}

Cheapest Insertion 法では、find_cheapest 関数で都市と挿入位置を同時に選びます。挿入位置の計算は find_insert_position 関数を利用しますが、挿入位置だけではなく距離の増分も返すように実装を変更します。

<?php
function solve($data) {
    ...
    for ($i = 1; $i < $n; ++$i) {
        [$next, $pos] = find_cheapest($data, $i, $answer);
        array_splice($answer, $pos, 0, [$data[$next]]);
        swap($data, $i, $next);
    }
    ...
}

function find_cheapest($data, $k, $answer) {
    $next = -1;
    $pos = -1;
    $min = INF;
    $n = count($data);
    for ($i = $k; $i < $n; ++$i) {
        [$p, $c] = find_insert_position($data[$i], $answer);
        if ($c < $min) {
            $next = $i;
            $pos = $p;
            $min = $c;
        }
    }
    return [$next, $pos];
}

function find_insert_position($node, $data) {
    ...
    return [$pos, $min_cost];
}

それぞれの方法で解を求める様子は下図のとおりです。左から順に、Nearest, Farthest, Cheapest です。
f:id:y_uti:20171104111710g:plain:w200 f:id:y_uti:20171104111738g:plain:w200 f:id:y_uti:20171104112238g:plain:w200

Greedy 法 (Kruskal の最小全域木構成アルゴリズム)

参考にした講義資料で説明されている Greedy 法は、都市間の距離が小さな順にルートの断片を組み立てていく方法です。これは最小全域木の特殊な場合 (ノードの次数が 2 を超えない) と考えられるため、このような名前で説明されているのだと思います。PHP で次のように実装してみました。make_distances 関数で都市のペアを距離の昇順に列挙し、connect 関数によってルートに追加していきます。最後に build_answer 関数で、(0, 0) から順に都市間の接続情報を辿って最終的なルートを求めます。

<?php
function solve($data) {
    array_unshift($data, [0, 0]);
    $n = count($data);

    // $links[$i] : $data[$i] と接続されているノード番号
    // $data[0] を端点とするためダミーノードと接続済みとして処理する
    $links = array_fill(0, $n, []);
    $links[0][0] = PHP_INT_MAX;

    // $edges[$i] :
    //   $data[$i] の次数が 0 (未接続) の場合は $i
    //   $data[$i] の次数が 1 (端点) の場合は反対側の端点のノード番号
    //   $data[$i] の次数が 2 (中間点) の場合は -1
    $edges = range(0, $n - 1);
    $edges[0] = PHP_INT_MAX;

    // $distances : ノード間の距離を保持する二次元配列
    //   $distances[*][0] 一方のノード番号
    //   $distances[*][1] もう一方のノード番号
    //   $distances[*][2] ノード間の距離
    $distances = make_distances($data);
    usort($distances, function ($a, $b) { return $a[2] <=> $b[2]; });
    foreach ($distances as [$i, $j, $distance]) {
        connect($data, $i, $j, $edges, $links);
    }

    return build_answer($data, $links);
}

function make_distances($data) {
    $distances = [];
    $n = count($data);
    for ($i = 0; $i < $n; ++$i) {
        for ($j = $i + 1; $j < $n; ++$j) {
            $distances[] = [$i, $j, distance($data[$i], $data[$j])];
        }
    }
    return $distances;
}

function connect($data, $i, $j, &$edges, &$links) {
    // 既存の経路の中間点から枝分かれしてはいけない
    if ($edges[$i] == -1 || $edges[$j] == -1) return;
    // 経路の両端を繋いでループにしてはいけない
    if ($edges[$i] == $j || $edges[$j] == $i) return;

    // $i と $j が繋がるので "$i の反対端 ($ei)" の反対端は "$j の反対端 ($ej)" になる ($j も同様)
    $ei = $edges[$i];
    $ej = $edges[$j];
    $edges[$ei] = $ej;
    $edges[$ej] = $ei;
    // この接続により $i が中間点になる場合を考慮する ($j も同様)
    // $i が未接続だった場合は $i == $ei なので手前の処理で正しく更新済み
    if ($ei != $i) $edges[$i] = -1;
    if ($ej != $j) $edges[$j] = -1;

    $links[$i][] = $j;
    $links[$j][] = $i;
}

function build_answer($data, $links) {
    $answer = [];
    $prev = PHP_INT_MAX;
    $curr = 0;
    while (count($links[$curr]) == 2) {
        if (($next = $links[$curr][0]) == $prev) $next = $links[$curr][1];
        $prev = $curr;
        $curr = $next;
        $answer[] = $data[$curr];
    }
    return $answer;
}

この方法で解を求める様子をアニメーションで見てみます。Nearest Neighbour 法とも Insertion 法とも異なり、断片的な経路が少しずつ接続されて全体の経路が求められる様子がわかります。
f:id:y_uti:20171104124102g:plain

その他の方法

今回紹介した方法のほかにも、Insertion 法で全都市の凸包から始める方法、最小全域木を作ってから往復になる箇所を短絡していく方法などがあるようです。また、これらの方法で経路を作成した後に、経路を局所的に入れ替えて距離を短縮していく方法も知られています。記事の最初に挙げた参考資料では、それらの方法についても詳しく説明されています。

なお、きっかけとなった paiza の問題では、Random Insertion の後に 2-opt (今回の記事では未紹介) で改善する方法で登録し、最短距離*4よりも 4 % ほど悪い結果でした。

*1:ただし、出発点が固定されている点と、最後に出発点に戻らなくてよい点が、通常の巡回セールスマン問題とは異なっています。

*2:paiza の問題設定では、出発点が (0, 0) に固定されていて、(0, 0) 自体は入力データに含まれないことになっています。そのため、solve 関数の最初と最後で [0, 0] を追加、削除して問題に合わせています。

*3:もちろん、このような面倒なことをしなくても最初から Python で書けばよいです。

*4:正確には優勝者の距離です。同一の距離で優勝している方が 20 名以上いるので、おそらく厳密解に到達できているのだと思います。