https://smooth-pudding.hatenablog.com/entry/2024/06/18/001629

この記事は以下の記事の続編です。もし読んでいない方は先に前回をお読みください。
smooth-pudding.hatenablog.com

前回「やり切ったぜ～」的なことを言ってたんですが、引き続きいじっていたら結構速くなりました。約3倍です。
前回↓

今回↓

本当にいろいろな工夫をしたので、すべてを書くわけにはいきませんが、今回はその一部を紹介したいと思います。

特徴量の工夫

前回の記事では、以下の組み合わせを使いました。

Morgan 特徴量
固有多項式の $N-3$ 次の係数 ( $N$ は行列のサイズ)

これの精度を上げるため、少し工夫します。

まず、特徴量とはどういうものか、数学的に整理しておきます。 $e_{i}$ $(i = 0, 1, \dots, N - 1)$ を標準基底とします。すなわち

$e_{0}=\begin{pmatrix}1 \\ 0 \\ 0 \\ \vdots \\ 0\end{pmatrix}, e_{1} = \begin{pmatrix}0 \\ 1 \\ 0 \\ \vdots \\ 0\end{pmatrix}, e_{2} = \begin{pmatrix}0 \\ 0 \\ 1 \\ \vdots \\ 0\end{pmatrix}, \dots, e_{N-1} = \begin{pmatrix}0 \\ 0 \\ 0 \\ \vdots \\ 1\end{pmatrix}$ とします。対称群 $S_{N}$ （すなわち $\{ 0, 1, \dots, N - 1 \}$ の並べ替え方全体の集合）の元 $\sigma$ に対し, $(e_{0}, e_{1}, \dots, e_{N-1}) \mapsto (e_{\sigma(0)}, e_{\sigma(1)}, \dots, e_{\sigma(N-1)})$ という写像を定義することができます。気持ちとしては、これがノードの番号の付け替えに対応します。

隣接行列全体の集合を $M$ とします。隣接行列に対する写像 $c \colon M \to \mathbb{Z}^{N}$ が不変量であるとは、次の意味で基底の入れ替えに対して不変であることをいいます。すなわち $\sigma \in S_{N}$ に対して基底変換の行列 $P_{\sigma}$ を

$P_{\sigma} = \begin{pmatrix}e_{\sigma(0)} & e_{\sigma(1)} & \dots & e_{\sigma({N-1})}\end{pmatrix}$ と定義するとき、任意の $A \in M$ に対して $c_{\sigma(i)}(A) = c_{i}(P_{\sigma}^{T} A P_{\sigma}) \quad (\text{$c_{i}$ は $c$ の第 $i$ 成分})$ が成り立つときをいいます。直観的にいえば、ノードの番号の付け方を入れ替えたとしても、そのノードに対応する値は変わらない、という定義です。前回の記事で導入した特徴量は、この不変量という言葉で言い換えられます。

前回紹介した Morgan 特徴量が不変量であることを確認しておきます。 $s \geqq 1$ を固定します。 $A \in M$ に対して

$c(A) = A^{s} \begin{pmatrix}1 \\ 1 \\ \vdots \\ 1\end{pmatrix}$ と定義すると $c_{i}(A) = e_{i}^{T} A^{s} \displaystyle\sum\limits_{j = 0}^{N - 1} e_{j}$ と書き直せます。このとき $\begin{align}c_{\sigma(i)}(A) &= e_{\sigma(i)}^{T} A^{s} \sum\limits_{j = 0}^{N - 1} e_{j}\\ &= e_{\sigma(i)}^{T} A^{s} \sum\limits_{j = 0}^{N - 1} e_{\sigma(j)}\\ &=e_{i}^{T} P_{\sigma}^{T} A^{s} P_{\sigma} \sum\limits_{j = 0}^{N - 1} e_{j} \\ &= c_{i}(P_{\sigma}^{T} A^{s} P_{\sigma})\end{align}$ が成り立ちます。この議論は $s$ に依らないので、無事証明できました。

今回のコードでは、以下の2つの不変量を利用しました。

$c_{i}(A) = e_{i}^{T} A^{s} e_{i}$
$c_{i}(A) = \displaystyle\sum_{j = 0}^{N - 1} \left( e_{i}^{T} A^{s} e_{j} \right)^{2}$

それぞれ不変量であることを示しておきます*1。

$c_{i}(A) = e_{i}^{T} A^{s} e_{i}$
次の計算によります。 $\begin{align}c_{\sigma(i)}(A) &= e_{\sigma(i)}^{T} A^{s} e_{\sigma(i)} \\ &= e_{i}^{T} P_{\sigma}^{T} A^{s} P_{\sigma} e_{i} \\ &= c_{i}(P_{\sigma}^{T} A P_{\sigma})\end{align}$
$c_{i}(A) = \displaystyle\sum_{j = 0}^{N - 1} \left( e_{i}^{T} A^{s} e_{j}\right)^{2}$
次の計算によります。 $\begin{align}c_{\sigma(i)}(A) &= \sum_{j = 0}^{N - 1} e_{\sigma(i)}^{T} A^{s} e_{j} e_{j}^{T} A^{s} e_{\sigma(i)} \\ &= e_{\sigma(i)}^{T} A^{s} \left( \sum_{j = 0}^{N - 1} e_{j} e_{j}^{T} \right) A^{s} e_{\sigma(i)} \\ &= e_{\sigma(i)}^{T} A^{s} \left( \sum_{j = 0}^{N - 1} e_{\sigma(j)} e_{\sigma(j)}^{T} \right) A^{s} e_{\sigma(i)}\\ &= \sum_{j = 0}^{N - 1} e_{\sigma(i)}^{T} A^{s} e_{\sigma(j)} e_{\sigma(j)}^{T} A^{s} e_{\sigma(i)} \\ &= \sum_{j = 0}^{N - 1} e_{i}^{T} P_{\sigma}^{T} A^{s} P_{\sigma} e_{j} e_{j}^{T} P_{\sigma}^{T} A^{s} P_{\sigma} e_{i} \\ &= c_{i}(P_{\sigma}^{T} A P_{\sigma}) \end{align}$ ただし3つ目の等式では $\sum_{j} e_{j} e_{j}^{T}$ と $\sum_{j} e_{\sigma(j)} e_{\sigma(j)}^{T}$ がいずれも単位行列であることを用いました。

今回のコードでは、この二つの不変量を $s = 1, 2, 3$ に対して計算し、3要素の列によって特徴量としました。前回は $s = 1, 2, \dots, 10$ で計算していたためこれだけでも省エネですが、さらに行列を分解する能力も高まりました。一度に比較しなければならない行列の数が抑制されたことにより、メモリ効率が良くなり、結果として計算時間が飛躍的に縮まりました。

枝切り判定の工夫

flamegraph を使って性能解析をしていると、どうやら全列挙の処理の時間のうち、枝切り判定の関数でかなりの時間を使っていることが分かりました。そこで枝切り判定部分にもメスを入れました。

枝切り判定の考え方については特に前回から変更していません。今回変更したのは、枝切りの評価を行う回数の部分です。よくコードを観察してみると、同じ行列に対して何度も判定を行っている箇所がありました。具体的には...

0 を入れた直後も判定していた
→ 全く同じ行列への評価なので省略
多くの 0 や 1 を埋めていった後の方になっても、最初の方の数字の正当性をチェックしていた
→ 最初の部分の確認はスキップ
状態の変わっていない行に対する count_bits を何度も計算していて重い
→ 可変借用の引数を使って共有
比較が重い
→二つの整数の比較を or で結ぶ代わりに、一つの整数にまとめて比較する形に変更

これらの工夫を行って、処理がそれなりに圧縮されました。

以下に変更前後の枝切り判定関数を引用しておきます。

// 変更前
fn is_fine_up_to(
    matrix: &SymmetricBitMatrix<N>, 
    idx_row: usize,
) -> bool {
    let mut prev_row = u16::MAX;
    let mut prev_nbits = 4;
    for (i, &row) in matrix.rows().iter().enumerate() {
        let nbits = row.count_ones();
        if nbits > prev_nbits || (nbits == prev_nbits && row > prev_row) {
            return false;
        }
        if i < idx_row {
            prev_row = row;
            prev_nbits = nbits;
        }
    }
    idx_row < N - 1 || matrix.is_connected()
}

// 変更後
fn is_fine_up_to(
    matrix: &SymmetricBitMatrix<N>,
    row: usize,
    next_row: usize,
    conf: &mut u32, // 最終確定行の情報を保持する. 前半はノードの次数, 後半は隣接ノードの情報
) -> bool {
    for (i, &row_bits) in matrix.rows().iter().enumerate().skip(row) {
        let this = (row_bits.count_ones() << 16) | (row_bits as u32);
        if this > *conf {
            return false;
        }
        if i < next_row {
            *conf = this;
        }
    }
    next_row < N - 1 || matrix.is_connected()
}

入れ替えパターン列挙の工夫

行列同士が同じ構造を表すかどうかの判定をするためには、最終的にノードの入れ替えを行う必要があります。どう入れ替えるべきかについては特徴量によって制限されているので、特徴量のパターンに沿った入れ替え方を考える必要があります。この入れ替え方を事前にすべて列挙しておいて、すべてのスレッドで使いまわすというのは、前回にすでに行っていた工夫でした。

ところで、この部分の生成には itertools という外部クレートにある、順列を求める関数を使っていました。
docs.rs

これ自体はよくできているクレートなので、普通の使い方をする分には問題ないものだと思います。今回は無茶苦茶に高速化しようとしているので、一応実装の中身を確認してみました。すると、もしかして一般的な状況に対応するために若干効率を犠牲にしているんじゃないか？という臭いがしました。

そこで、連続する整数に限って、すべての並べ替えを列挙するイテレータを自分で実装しなおすことにしました。といっても、元々どう実装すべきか分からなかったので、ChatGPT に相談しました。

いきなり正解を出してきました。これを読んで理解してから、いざ実装しようとしてみると、なんと関数名を書いたぐらいの時点で GitHub Copilot がほとんどコード生成してしまいした。恐るべし、AI-aided 開発。

これにより、最後の壁だった2秒の壁を越えることができました。

せっかくなので、半分（というか9割）AIが書いてくれた順列全列挙イテレータのコードを貼っておきます。

#[derive(Debug, Clone)]
struct PermutationGenerator {
    current: Vec<usize>,
    is_first: bool,
}

impl PermutationGenerator {
    fn from_start_count(start: usize, count: usize) -> Self {
        let current = (start..(start + count)).collect();
        Self {
            current,
            is_first: true,
        }
    }
}

impl Iterator for PermutationGenerator {
    type Item = Vec<usize>;

    fn next(&mut self) -> Option<Self::Item> {
        if self.is_first {
            self.is_first = false;
            return Some(self.current.clone());
        }

        let mut i = self.current.len() - 1;
        while i > 0 && self.current[i - 1] >= self.current[i] {
            i -= 1;
        }
        if i == 0 {
            return None;
        }

        let mut j = self.current.len() - 1;
        while self.current[j] <= self.current[i - 1] {
            j -= 1;
        }

        self.current.swap(i - 1, j);
        self.current[i..].reverse();
        Some(self.current.clone())
    }
}

さらなる高速化？

正直、ほぼ思いつく限りの高速化はやってしまったので、自分の力でこれ以上高速化できる余地があるのか自信がありません。唯一、今手を付けていない（そして付けるつもりもない）ところが、ハッシュ関数の工夫の部分です。現時点でも FxHasher というかなり高速（だがそれなりに高性能）なハッシュ関数を用いているので、これを越えるのはそう簡単ではないと思いますが、行列同士の比較の部分などで HashSet は大量に使っているので、うまくいけば効果はあるのではないでしょうか。私はやらないので物好きなアナタにお任せします。

今回の記事の計算で用いたコードは以下にあります。が今後変更するかもしれません。
github.com

なおこの記事を執筆している時点での最新のコミットはこちらです。
github.com

今度こそ終わりです。ではまた。

追記：
終わりだと言ったな。あれは嘘だ。

以下を追加で行いました:

行列の列挙を並列化
行入れ替えパターンの列挙を並列化
スレッド数のチューニング

そしたら1秒切って0.8秒台になりました。ほげぇ。

↑を行ったコードは以下に置きました。多分もう編集しないです。多分。知らんけど。
github.com

今度こそホントに終わりです。多分。知らんけど。

*1:実は任意の関数 $\psi, \phi$ および整数係数の多項式 $p$ に対して

$\displaystyle c_{i}(A) = \psi(e_{i}^{T} p(A) e_{i}) + \sum_{j \neq i} \phi( e_{i}^{T} p(A) e_{j} )$ が不変量であることが証明できます。対称群の元が互換の積に分解できることを利用し、各互換に対して不変であることを示せば従います。この記事で紹介している不変量は、いずれもこれの特別な形です。