平衡二分探索木を使った2-optの効率化について

こんにちは。 平衡二分探索木で2-optを高速化する方法について研究したので、本記事ではそれについて書こうと思います。

2-optとは

2-optは巡回セールスマン問題(TSP)の近似解を求めるヒューリスティックアルゴリズムです。
下図のように、サイクルに対してランダムに二つのパスを選び、パスを組み替えた方がサイクルの長さが短くなる場合、パスを組み替えます。具体的には、A→B + C→D > A→C + B→D の時、辺ABと辺CDを辺ACと辺BDに組み替えます。
f:id:aspirator:20210325193100j:plain しかし、辺を張り直すだけではサイクルにならないので、B~CまたはD~Aの間の順序を逆にする必要があります。
2-opt法では順路を配列などで管理することでこれを実現しています。
f:id:aspirator:20210325193138j:plain そしてこれが`測定で用いた2-optの実装です。

#include<iostream>
#include<algorithm>
#include<vector>
#include<iterator>
#include<random>
#include<chrono>
using namespace std;
using time_point = chrono::system_clock::time_point;


time_point get_time() {
    return chrono::system_clock::now();
}
int time_diff(chrono::system_clock::time_point start) {
    return static_cast<int>(chrono::duration_cast<chrono::milliseconds>(chrono::system_clock::now() - start).count());
}

std::random_device seed;
std::mt19937 random_number(seed());

struct point {
    int x, y;
};
vector<point>points;
double dist(point a, point b) { //点a, b間の距離
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
int loop_count = 0;
void two_opt(vector<int>& route, int limit_ms) {
    time_point start_time = get_time();
    int a, b, size = route.size();
    while (time_diff(start_time) < limit_ms) {
        loop_count++;
        a = random_number() % size;
        b = random_number() % size;
        if (a == b)continue;
        if (b < a) swap(a, b);
        if (dist(points[route[a]], points[route[(a + 1) % size]]) + dist(points[route[b]], points[route[(b + 1) % size]]) >
            dist(points[route[a]], points[route[b]]) + dist(points[route[(a + 1) % size]], points[route[(b + 1) % size]])) {
            reverse(route.begin() + a + 1, route.begin() + b + 1);
        }
    }
}

int main() {
    const int N = 10000;
    points.resize(N);
    std::vector<int>route(N);
    for (int i = 0; i < N; i++) {
        points[i] = { random_number() % 10000,random_number() % 10000 };
        route[i] = i;
    }
    two_opt(route, 2000);
    cout << loop_count << endl;
}

しかし、順序を反転させる操作にはO(N)かかってしまいます。
この記事では平衡二分探索木を用いることで反転操作をO(\log N)にします。
今回は平衡二分探索木としてtreapを使いました。
長いので別に載せておきます。
treapのライブラリ
そして測定で用いたコードがこちらです。

#include <iostream>
#include <vector>
#include <utility>
#include <cmath>
#include <random>
#include <chrono>
#include "reversibleArray.cpp"
using namespace std;
using time_point = chrono::system_clock::time_point;


time_point get_time() {
    return chrono::system_clock::now();
}
int time_diff(chrono::system_clock::time_point start) {
    return static_cast<int>(chrono::duration_cast<chrono::milliseconds>(chrono::system_clock::now() - start).count());
}

std::random_device seed;
std::mt19937 random_number(seed());

struct point {
    int x, y;
};
vector<point>points;
double dist(point a, point b) { //点a, b間の距離
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
int loop_count = 0;
void original_two_opt(vector<int>& route, int limit_ms) {
    time_point start_time = get_time();
    int l, r, size = route.size();
    reversible_array opt(route);
    while (time_diff(start_time) < limit_ms) { 
        loop_count++;
        l = random_number() % size;
        r = random_number() % size;
        if (l == r)continue;
        if (r < l) swap(l, r);
        int c = opt[l], d = opt[(l + 1) % size],
            e = opt[r], f = opt[(r + 1) % size];
        if (dist(points[c], points[d]) + dist(points[e], points[f]) >
            dist(points[c], points[e]) + dist(points[d], points[f])) {/*組み替えた方が短くなる場合反転*/
            opt.reverse(l, r);
        }
    }
    for (int i = 0; i < size; i++) { /*平衡二分探索木で出た結果をrouteに戻す*/
        route[i] = opt[i];
    }
}

int main() {
    const int N=10000;
    points.resize(N);
    std::vector<int>route(N);
    for (int i = 0; i < N; i++) {
        points[i] = { random_number() % 10000,random_number() % 10000 };
        route[i] = i;
    }
    original_two_opt(route, 2000);
    cout << loop_count << endl;
}

計測

計測方法

計測はwhileループを2000ms回してループした回数を数えます。
それを5回繰り返して平均値をとります。

計測結果

頂点数とループ回数

頂点数 10^{2} 10^{2.5} 10^{3} 10^{3.5} 10^{4} 10^{4.5} 10^{5} 10^{5.5} 10^{6}
2-opt 1.07\cdot 10^{6} 1.20\cdot 10^{6} 1.17\cdot 10^{6} 1.13\cdot 10^{6} 9.56\cdot 10^{5} 4.67\cdot 10^{5} 9.68\cdot 10^{4} 3.04\cdot 10^{4} 9.84\cdot 10^{3}
オリジナル 1.69\cdot 10^{5} 1.37\cdot 10^{5} 1.12\cdot 10^{5} 9.29\cdot 10^{4} 7.61\cdot 10^{4} 6.41\cdot 10^{4} 5.80\cdot 10^{4} 5.14\cdot 10^{4} 4.48\cdot 10^{4}

f:id:aspirator:20210830120559j:plain

考察

通常の2-optでは頂点数が10^{2}~10^{4}の範囲でループ回数が1.1\cdot 10^{6}に固まっています。
これはある程度解が改善するとそれ以上改善されにくくなり、ボトルネックであるrevesrseパートが行われにくくなるからだと考えられます。
また、頂点数が10^{5}未満の時、通常の2-optにかなわなかったのも、同様の理由から取得O(\log N)の平衡二分探索木より取得O(1)vectorの方が高速になったからだと考えられます。
結果の表やグラフから、頂点数が多い場合やタイムリミットが短い場合は平衡二分探索木を用いた2-optを使った方が効率が良く、頂点数が少ない場合やタイムリミットが長い場合は普通の2-optを使う方が効率が良いことが分かりました。

途中まで平衡二分探索木を用いて、ある程度解が良くなったところで通常の配列を用いた2-optを使うと、さらに効率が良くなるかもしれません。

感想

自分からこんな研究をするなんて自分でもびっくりしました。
ラソンコンテストへの採用は出来なさそうだったので、またいつか細かいところの定数倍高速化をしたり、続きの研究をしてリベンジしたいと思います。

夏休みの自由研究の課題としてこれを提出しました。