🏆
Kaggle Santa 2022 コンペ参加記

iizuka kota
2023/04/21に公開
Kaggle
最適化
tech
 はじめに2022年12月に開催された Kaggle コンテスト「 Santa 2022 - The Christmas Card Conundrum 」にて、私 (kota.iizuka) は 19位 で銀メダルを獲得し、 Kaggle Master へ昇格しました！ この記事では、コンテストの概要と解法について紹介します。

 コンテストの概要Santa コンテストは例年12月に Kaggle で開催される数理最適化コンテストで、今回は「与えられた画像の各ピクセルをなるべく無駄なく塗りつぶすような経路」を求める問題でした。

 入力と出力クリスマスコンペということもあり、入力画像はクリスマスにちなんで次の画像 $ I $ が与えられました。サイズは 257 × 257 で中心座標が (0, 0) 、 (x, y) の値域は [-128, 128] になっています。
出力は、サイズが N × 8 × 2 の配列 X で、次の（やや複雑な）条件を満たすことが必要です：
sum(X, axis=1) で得られる N × 2 の配列 R について、任意の整数 x, y\in [-128, 128] に対して R[i] == [x, y] になるような i が存在する
max(abs(X), axis=2) で得られる N × 8 の配列 A について、すべての i で A[i] == [64, 32, 16, 8, 4, 2, 1, 1] が成り立つ。
すべての i に対して、 max(abs(X[i+1] - X[i])) &lt;= 1 が成り立つ。
X[0] と X[-1] は [[64, 0], [-32, 0], [-16, 0], [-8, 0], [-4, 0], [-2, 0], [-1, 0], [-1, 0]] である。
条件1で現れる配列 R は、各画素を頂点とした巡回路を表しています。条件4と組み合わせると、始点と終点が (0, 0) であるような巡回路を表していることが分かります。また、2回以上同じ点を踏むようなパスは制限されていません。
条件2の配列 A は、各頂点を指すためのアームの配置に関する制約です。例として X[i] = [[64, 30], [20, -32], [-6, -16], [8, 3], [4, -2], [2, 0], [1, 1], [0, -1]] であるような配列 X を考えてみます。 X[i] を画像の上に描画すると次の太線のようになります。
一方で、この配列は A[i] == [max(64, 30), max(20, 32), max(6, 16), max(8, 3), max(4, 2), max(2, 0), max(1, 1), max(0, 1)] と計算できるので条件2を満たしています。この条件は、それぞれのアームの端点を、前のアームの端点から正方形の領域（図で細い線で描いた位置）に制約する、という条件であると解釈できます。
条件3は、隣り合うアームの配置が大きく変わるような構成を許さないような制約です。とくに、条件4から始点と終点が (0, 0) なので、パスの最初と最後の近くではかならず (0, 0) の周辺を通るパスになります。

 評価関数最小化したい評価関数は、アームの移動によるコスト cost_move と頂点の色変化によるコスト cost_color の和で表されます。各頂点のコストはそれぞれ下記のように定式化できます。
cost_move(X[i], X[i+1]) = sqrt(sum(abs(X[i+1] - X[i])))
アームを n 個動かすときに \sqrt{n} のコストがかかる。条件3から1回の移動あたりのコストはたかだか \sqrt{8} 。
cost_color(R[i], R[i+1]) = 3 * sum(color(R[i]) - color(R[i+1]))
RGB 値（値域は [0, 1] ）の差がそのままコストになる。1回の移動あたりのコストは最大で 3\times (1+1+1) = 9。

 手法私が序盤に検討していたのは局所探索に基づく手法（手法1）でしたが、検討する中で全体の最適化もできそうであると考えて方針転換しました（手法2）。

 手法1. 四分割して局所最適パスを求める方法まず、アームの位置によっては移動できない方向が存在することについて考えました。たとえば、初期位置 [[64, 0], [-32, 0], [-16, 0], [-8, 0], [-4, 0], [-2, 0], [-1, 0], [-1, 0]] から x\lt 0 の方向に移動しようとするとアームの長さの制約に反してしまうので、最初の移動は必ず y 方向に行われる必要があります。
この制約をはじめとして、アームの動きが可能であるかどうかは頂点だけでなくアーム全体に依存しているので、いったんこれをルールベースで回避する方法を考えました。左上（x\geq 0, y\geq 0）のみ考えることにすると、最初のアームを x=64, 2番目以降を y = [32, 16, 8, 4, 2, 1, 1] が成り立つようにとり、 y 方向の移動は1番目のアーム、 x 方向の移動は 2番目以降のアームを使って、任意の座標から8近傍に移動できるようなアームの構成が得られます。（文章だけだと分かりにくいですが こちらの記事 のアニメーションが分かりやすいです）
この構成で、単純に画面全体を走査するだけでも 81147 点となり、公式のベースライン （166305 点）と比べて半分程度までスコアを更新できます。各辺のコストを色で表した図は次のようになっています。右上の雪だるまの白と黒の境界や、右上のトナカイと背景の境界などのコストが大きいので、距離のコストよりRGB値によるコストが全体的に目立ちます。
次に、4分割したそれぞれに対してコストを最小化するパスを探します。今回は、問題に合ったコスト行列を作成し、 LKH-3 という巡回セールスマン問題のソルバに入力しました。（LKH-3 はコストを整数で表現する必要があるので、画素値が正確に計算できるように全体を 255 倍しました。辺のコストが \\sqrt{2} など 1/255 の整数倍にはなっていないためここで多少の誤差が発生しますが今回は無視しました。また移動不可能な辺間には大きいコストを与えてその辺が採用されないようにしています）
この方法で、 kaggle notebook 上で30分程度かかりますが 74257 点という銀メダル相当の解が得られます。実際に動作する notebook は こちら から確認できます。先ほどの単純な解に比べてRGB値が異なるパスをなるべく通過しないようになっています。

 手法2. 最適パスを求めたあとにアームを構成する方法手法1では x=0 または y=0 を通るようなパスが人工的に制限されていますが、これをどの程度まで緩めることができるかについて考察しました。その材料として、まずは y=0 の周辺で x の符号が変わるような移動を制限するようなパスを LKH-3 で作成し、そのパスに対してうまくアームをあてはめられるかを検討しました。
6万頂点以上あるグラフについて LKH-3 で辺のコストを具体的に入力すると数十GBにもなってしまうので、 LKH-3 のソースコードを書き換えて辺のコストを実行時に計算できるように修正します（RGB値のコストを3倍するのは配列の定義時に行っています）。さらに、 -15 \lt y \lt 15 で y 軸を横切るようなパスを制約する項を追加します。コードとしては下記のようになります。
with open("LKH-3.0.8/SRC/Distance_SPECIAL.c", "w") as f:
    print('#include "LKH.h"', file=f)
    print("const int image[3*257*257] = {", file=f)
    for i, c in enumerate("rgb"):
        print(f"// color {c}", file=f)
        for j in range(257):
            print("".join([f"{int(255*3*image[j, k, i]):3}," for k in range(257)]), file=f)
    print("""};

const int pos_cost[] = {
    0,
    255,
    (int)(255 * sqrt(2)),
    (int)(255 * sqrt(3)),
    (int)(255 * sqrt(4)),
    (int)(255 * sqrt(5)),
    (int)(255 * sqrt(6)),
    (int)(255 * sqrt(7)),
    (int)(255 * sqrt(8)),
};

int Distance_SPECIAL(Node * Na, Node * Nb)
{
    const int ax = (int)(Na->X);
    const int ay = (int)(Na->Y);
    const int bx = (int)(Nb->X);
    const int by = (int)(Nb->Y);
    int dx = ax - bx;
    int dy = ay - by;
    if (dx < 0) dx = -dx;
    if (dy < 0) dy = -dy;
    // アームの差が長すぎるパスは作れない
    if (dx + dy > 8) return 10000;
    // 自由度が低いパスは作れない
    if ((((ax >= 0) && (bx < 0)) || ((ax < 0) && (bx >= 0))) && (((-15 < ay) && (ay < 15)) || ((-15 < by) && (by < 15)))) {
        return 10000;
    }

    // color_cost
    const int ia = 257*(ax + 128) + (ay + 128);
    const int ib = 257*(bx + 128) + (by + 128);
    int dr = image[257*257*0 + ia] - image[257*257*0 + ib];
    int dg = image[257*257*1 + ia] - image[257*257*1 + ib];
    int db = image[257*257*2 + ia] - image[257*257*2 + ib];
    if (dr < 0) dr = -dr;
    if (dg < 0) dg = -dg;
    if (db < 0) db = -db;
    const int color_cost = dr + dg + db;

    return pos_cost[dx + dy] + color_cost;
}""", file=f)
このコスト関数を使用してパスを計算すると、次の図のような結果が得られます（上の図と異なり、パス全体における位置を色で表示しています）。原点付近で x\lt 0 へ向かうパスが制限されていることが図から確認できます。このパスに合うアームの構成が発見できたとするとスコアは 74082.4 と 16 位に相当します。
あとはこのパスに合ったアームの構成を探すだけ、ということでパスを探索するアルゴリズムを考えました。アームの配置の評価関数としては、単純に各方向に進める数を計算した arm_capacity 関数のほかに、（たとえばスタート地点のような）斜辺が少なく同じ向きの辺が多いと動きにくいことから、 diagness 関数と same_arg_penalty 関数が考えられます。
def arm_capacity(arm, delta=1):
    # 左(x-)、右(x+)、上(y+)、下(y-)にあと何マス進めるか（sum, min が大きいほうが良い）
    dl, dr, dt, db = 0, 0, 0, 0
    for r, (x, y) in zip(TORUS_RADIUS, arm):
        if abs(x) == r:
            dt += min(delta, r - y)
            db += min(delta, r + y)
        if abs(y) == r:
            dl += min(delta, r - x)
            dr += min(delta, r + x)
    return dl, dr, dt, db

def diagness(arm, coeff=[1, 1, 1, 1, 1, 1, 1, 1]):
    # 斜辺がたくさんあるほうが良い（値が大きいほうが良い）
    ret = 0
    for r, (x, y), c in zip(TORUS_RADIUS, arm, coeff):
        ret += c * min(abs(x)/r, abs(y)/r)
    return ret

def same_arg_penalty(arm):
    # 同じ向きの辺が多いと詰まりやすいので避ける
    theta = np.array([np.arctan2(y, x) for x, y in arm])
    return -abs(np.average(np.exp(1j * theta)))
現在のアームの配置から次の座標を表現できるアームを全列挙した中で、これらの重みづけ和が最も大きいものを選択することを繰り返していきます。次のアームが存在しないような配置に当たってしまったら、少し（数百頂点程度）戻ってから重みを変えて再度探索を開始します。
このアルゴリズムを検証するために、まずは重み変更については手動でパスを見ながらやることにして実行確認をしたのですが、その段階でほとんどの場合動作することが確認できました（その意味ではアームに関する制約はそこまで強くはないことがわかります）。ただし、端から 2 列目（ x, y=\pm 127 ）を横切るようなパスについては何回繰り返しても正しいパスが見つかりませんでした。よく考えると、このパスを実現するようなアームの配置は（この方法に限らず、理論的に）存在しないことが分かります。
その分の修正を加えてパスの探索から再実行したものが次の図になります。制約が増えた分スコアは 74087.0 と下がってしまいましたが、それでも四分割する解法に比べて 150 点程度更新できています。

 上位解法上位4チーム（1位・2位・3位・4位）は同じスコアでおそらく最適解に到達しているのですが、4チームとも主に遺伝的アルゴリズム GA-EAX を利用した探索を行っており、私が使用した LKH-3 のような局所探索に基づく方法はあまり使われていませんでした。数万点を超えるような規模の巡回セールスマン問題で（今回の原点周辺の制約のように）複数の頂点に依存するような制約がある場合は、遺伝的アルゴリズムを使用したほうが良い傾向があるようです。私は解法が公開されるまでそういったアルゴリズムごとの傾向を知らず、以前から数理最適化コンペでよく使われている LKH-3 を使わないような方針についてあまり考えていませんでしたが、金メダルあるいは賞金圏内に入るためには、そういった先入観にとらわれずに適切なアルゴリズムを調査する技術を高めることが課題であると感じました。

 まとめ巡回セールスマン問題は一般的には NP 完全で難しい問題ですが様々なヒューリスティックが知られており、今回のように特定の制約が決まっているときに数万頂点で最適解を得られることがある非常に興味深い分野だと思いました。引き続き kaggle の様々な問題に取り組んで問題解決能力を高めていきたいです。

 参考リンクLKH-3 http://webhotel4.ruc.dk/~keld/research/LKH-3/
今回使用した巡回セールスマン問題のソルバです。局所探索に基づく Lin-Kernighan アルゴリズムが使われています。