この記事はなに？

私は現在4x4x4ルービックキューブ(ルービックリベンジ)を解くロボットを作り、世界記録を目指しています。ロボットを作る上で一番のハードルであったのが4x4x4キューブを現実的な時間に、なるべく少ない手数で解くアルゴリズムを実装することです。
調べてみると4x4x4については文献や先人が少ないことがわかると思います。そんな貴重な資料の一つになれば嬉しいと思い、この記事を書いています。
GitHubはこちら
この記事の内容は完全ではありません。一部効率の悪いところを含んでいるかもしれません。今後改良したら随時追記していきます。

↓競技用4x4x4ルービックキューブとプログラム制作のために番号を振られた競技用4x4x4ルービックキューブ

全貌

この記事集は全部で3つの記事から構成されています。

概要
アルゴリズム
実装(本記事)

この記事で話すこと

この記事では前回に説明したアルゴリズムを実際に実装するときのやり方、およびその際のTIPSをお話しします。実装について手取り足取り説明するのは冗長になってしまうので避けます。こちらの記事に2x2x2ルービックキューブについての解説を詳しく書きました(実際のコードも載っています)のでこちらも参照いただければと思います。

インデックス化

細かくフェーズを分けて探索することで、探索すべき領域が非常に小さくなりました。そのおかげで考えうるそのフェーズの状態に順番に番号をつけられるようになります。これによっていちいちパズルの動きをシミュレートするなどというコストの大きい作業をする必要がなくなり、ただただ配列参照をするだけでパズルを動かすのと等価な作業が行えます。
ただし例えばEPとCPを一緒に考えると$10^6-10^7$個前後の要素数の配列に収まらず、メモリをとても食べてしまいます。そこで例えばEPだけ、CPだけで番号付けをし、2つの数字でパズルの状態を表したりするフェーズもあります。

番号の付け方はなんでも良いのですが、わかりやすいようにEPとCPであればパズルの順列番号(階乗進数)、EOとCOであればそれぞれ2進数、3進数、センターであれば重複あり階乗進数で表すと良いと思います。階乗進数についてはこちらがわかりやすいです。

IDA*探索

実際の探索ではIDA探索(Iterative Deepening A)と呼ばれる探索を使うことになると思います。これまで幾度となく私の記事の中で登場している文句ですが、こちらの記事から引用すると、

IDAを一言で表すと、「最大深さを制限した深さ優先探索(DFS)を、深さを増やしながら繰り返す」です。
IDAのからくりは、一度探索したノードを忘れることにあります。ノードを忘れてしまえば、メモリを解放できますよね。

IDA*では深さ$N-1$までで深さ優先探索を行って解が見つからなかったら、最大深さを$N$に増やしてまた一から探索をやり直します。こうすることで、

返される結果は必ず最短手数
メモリ使用量がとても少ない

という恩恵があります。

ただし、深さの浅いノード(状態)については何度も同じ探索を繰り返してしまうため、計算量は若干増大します。しかしパズルの状態は深さに対して指数関数として増大するため、計算量の増大はそこまで大きくありません。

IDA*のミソは、現在の状態からフェーズ完成までの「距離」を推測することにあります。この「距離」は今回の場合は解けるまでの手数です。具体的には、パズルの状態をインデックス$i$で表すとすれば、

$depth \geq f(i) = g(i) + h(i)$

となるように探索を進めます。詳しく説明しましょう。

$depth$は探索する際の深さ(=手数)の上限です。
$f(i)$はインデックス$i$の状態を経過する解が何手かかるのかの推定値です。
$g(i)$は初期状態からインデックス$i$の状態に至るまでにかかった手数です。
$h(i)$はインデックス$i$の状態からフェーズ完成状態までにかかる手数の推定値です。

この式を満たしながら$depth$を0から1ずつ増やしていけば、手数の短い解が見つかります。
なお、インデックス$i$の状態からフェーズ完成にかかる実際の手数を$h^\star(i)$としたとき、常に$h(i)\leq h^\star(i)$を満たす場合に最適解が見つかります。つまり(そのフェーズ内では)最小手数の解が見つかります。この場合をadmissibleであると言います。

同じ状態を重複して探索しないようにする

同じ状態を何回も探索したら非効率ですよね。でも状態をどこかに保存しておくと膨大なメモリを食ってしまいます。そこで、ルービックキューブを解くにあたっては手順を最適化したら同じ手順になる場合を回避することでこれを実装します。
もちろんこれだけでは「違う手順を回したが同じ状態になる」場合を回避できません。しかし、これだけでもかなりの確率で同じ状態を探索せずに済むとのことです(3x3についてですがこちらに書いてあります。)。

4x4x4ルービックキューブの場合には、以下の2つの場合に探索をしないようにすると良いです。

同じ面を連続して回す場合
同じ軸を回す場合、内層->外層(または外層->内層: どっちにするかは最初に決めておく)の順番で回すとき

完成までの手数の推定値を返す関数

IDA*探索では$h(i)$が重要だと話しました。ここでは$h(i)$の算出手法について(私自身研究中ですが)軽くお話しします。

まず、各インデックスについて完成状態から何手で揃うのかを事前計算してテーブルにしておきます。私の場合はこれをcsvにしました。
実際にこの前計算した値を使う場合、往々にしてインデックスは複数あります($n$個あるとします)から、事前計算したテーブルから取り出せる値は$n$個あります。この$n$個の「距離」を使って、最終的にその盤面が完成からどれくらい遠いのかを表す距離を1つ出力しなくてはいけません。

やってみたことをひたすら羅列します。$n$個の距離を$L=(l_0, l_1, \dots l_{n-1})$とします。また、$\mathrm{sd}(L)$は$L$の標準偏差を表すとします。

$h(i)=\max(L)$
$h(i)=\sum L$
$h(i)=\max(L)+a\ \mathrm{sd}(L)$
$h(i)=\sqrt{\sum_j l_j^2}$
$t=\mathrm{pow}\left(a,-\left(\frac{\max(L)}{b\ \mathrm{sd}(L)}-c\right)^d\right)$として$h(i)=(1-t)\max(L)+t\sqrt{\sum_j l_j^2}$また、$\mathrm{sd}(L)=0$のときは$h(i)=\max(L)$

1つ目はadmissibleなことが保証されますが$h(i)$が$h^\star(i)$を大きく下回ることが多いため、計算量が大きくなりがちです。
2つ目はマンハッタン距離を使っていますが、admissibleを破りすぎてしまいます。簡単な例では同じ回転Rを回せば揃う$n$個の評価項目$(l_0, l_1, \dots l_{n-1})$について、$h^\star(i)=1$であるにも関わらず$h(i)=n$と返してしまいます。admissibleを破ると解の手数が増えるだけでなく、探索する空間が大きくなってしまうので計算量の増大にも繋がりがちです。
3つ目は、$L$のばらつきが大きいと$h^\star(i)$は大きくなる傾向があるという仮説のもと考えた式です。ですがあまりうまくいきませんでした。
4つ目はユークリッド距離です。2つ目よりはマシですがadmissibleを破ってしまう可能性が残っています。
5つ目が現時点で使っている方法です。定数$a, b, c, d$を調整して、$0<t<1$を作ります。この$t$を使って$L$の最大値とユークリッド距離を内分します。$t$は、$\max(L)$に比べて$\mathrm{sd}(L)$が大きいときにユークリッド距離に寄った値を返すという方針で計算しています。この関数$h(i)$を勝手にNyanyan's Functionと言っています。

疑似コード

私がCythonで書いたプログラムのPythonチックな擬似コードです。

def phase_search(phase, indexes, depth, h_i): # phase: フェーズ, indexes: パズルの状態インデックス, depth: 残り回せる手数, h_i: indexesの状態でのh(indexes)
    global path
    if depth == 0: # 残り手数0の場合には解にたどり着いたかを返す
        return h_i == 0
    twist = successor[phase][0] # twistは回す手
    while twist <= successor[phase][-1] : # successorは回転の候補
        if skip(twist, path): # 前の手で回したのと同じ面を回すなどはしない
            twist = skip_axis[phase][twist] # 同じ面を回さなくなるまでtwistを進める
            continue
        next_indexes = move(indexes, phase, twist) # 配列参照によってパズルを動かす。indexesの定義はフェーズによって異なるのでフェーズを引数にとる
        next_h_i = h(indexes, phase) # h(i)を返す。h(i)もフェーズごとに定義が異なるのでフェーズを引数に取る。
        if next_h_i > depth or next_h_i > h_i + 1: # 明らかに遠回りをしようとしている場合はスキップ
            twist = skip_axis[phase][twist]
        path.append(twist) # ここまで回してきた手順にtwistを追加する
        if phase_search(phase, next_indexes, depth - 1, next_h_i): # 再帰で次の手を探索する
            return True # 解が見つかった
        path.pop() # ここまで回してきた手順から今回した手順を取り除く
        twist = next_twist(twist, phase) # 次の手
    return False # 解が見つからなかった

def solver(puzzle): # puzzle: パズルのすべての状態を表したクラスのオブジェクトなど
    global path
    solution = [] # 解
    for phase in range(6): # フェーズをforで回す
        indexes = initialize_indexes(puzzle, phase) # パズルの状態をインデックスに変換
        h_i = h(indexes, phase)
        for depth in range(60): # depthを回す。なお60は適当
            path = [] # pathを空にする
            if phase_search(phase, indexes, depth, h_i): # IDA*探索を行う
                for twist in path: # フェーズが終わった状態までパズルをシミュレート
                    puzzle = puzzle.move(twist)
                solution.extend(path) # solutionに今のフェーズの解を付け加える
                break # 解が見つかったのでbreakして次のフェーズへ
    return solution # 解を返す

まとめ

ここまで3回の記事で私が4x4x4ルービックキューブを解くプログラムを書く上で学んだことをまとめて書いてきました。皆様の参考になれば幸いです。