https://smooth-pudding.hatenablog.com/entry/2024/06/05/235137

つい先日、技術書典16のオフライン会場に立ち寄りました。
techbookfest.org

前回の記事でも宣伝しましたが、同イベントで私も書籍を出しているので、ぜひご購入ください。
techbookfest.org

さて、そのイベント会場で、以下の書籍『炭素と水素による構造式』に出会いました（もちろんこちらも発売中です）。「アルゴリズムを工夫して、頑張って炭化水素を列挙するぞ！！」という、明らかに私に刺さるテーマの本だったので、気づいたら買っていました。
techbookfest.org

いろいろ工夫が満載で面白かったので、自分でもやってみました。というのが今回の記事の内容です。

例によって Rust で書いているので、ごちゃごちゃ御託はいいからコードを見せろ！！という人は以下のリンクから読んでみてください。この記事はこのコードの解説にあたります。
github.com

※本文中に引用したコードは最新版と異なる場合があります。

「炭化水素を列挙する」とは

炭化水素は、炭素と水素からなる分子です。メタン、エタン、ベンゼンなどがこれにあたります。化学的な性質は置いておいて、とにかく考えうる限りの炭化水素を列挙するぞ！というのが目標です。

もう少しルールを丁寧に書くとこんな感じです。

炭素は4本の手、水素は1本の手を持つ。
炭素同士や炭素と水素は手を繋いで結合する。
すべての手はいずれかの手とつながっている。
自分の手同士でつなぐことはない。
一つの構造ではすべての原子がつながっている。
立体的な構造は無視し、単につながりの関係だけに着目する。つまり立体異性体や鏡像異性体などは同一視する。

もちろん炭素や水素の数を増やしていけばいくらでも構造が考えられるので、ここでは「与えられた炭素の数の構造を全列挙する」という問題を考えることにします。

『炭素と水素による構造式』によれば、まずは単結合のみをもつ構造（=分子内のどの2原子も高々1本の手でしかつながっていないもの）を全列挙することがポイントです。そこでこの記事では、このような分子を全列挙する方法について考えます。

グラフ理論と隣接行列

よく考えてみると、炭化水素の構造を列挙することは、炭素のつながりを列挙することに他ならないことが分かります。なぜなら、炭素のつながりを考えたあとに、余った手にそれぞれ水素を結合させたものが炭化水素だからです。

炭素がいくつかあって、それらのつながり方を列挙したい。そのような状況にピッタリ当てはまる、グラフ理論という数学の理論があります。グラフ理論では、その名の通りグラフを扱うのですが、このグラフはいわゆる棒グラフや円グラフとは異なります。

グラフとは、いくつかの頂点（ノード）とそれらを結ぶ辺（エッジ）があるものが対象となります。下の図はグラフの一種です。このグラフはノードが4個、エッジが4本ありますね。

グラフのうちノードが炭素、エッジが炭素同士の結合と考えれば、これはひとつの炭化水素の構造を与えることになります。つまり、炭化水素の全列挙をグラフ理論の言葉で表現すれば、以下の条件を満たすグラフを列挙する問題となります。

ノードの数が与えられた個数である。
各ノードにつながるエッジの本数は高々4本である。
グラフ全体は分かれておらず、一つながりになっている（連結）。

さて、各ノードに 0, 1, 2, 3 の番号をつけてみます。

するとこのグラフは (0, 1), (1, 3), (0, 3), (2, 3) のつながりによって特徴づけられることが分かります。これを行列の形で次のように表現することができます。

$\begin{pmatrix}0 & 1 & 0 & 1 \\ 1 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 \\ 1 & 1 & 1 & 0 \end{pmatrix}$
これは各行・列に 0, 1, 2, 3 と番号を付けたうえで、「 $i$ 行 $j$ 列の成分が1」=「ノード $i$ とノード $j$ がつながっている」として表現したものです。このように、グラフの結合状態を表す行列のことを、隣接行列と呼びます。すなわち、グラフを列挙する問題は、隣接行列を列挙する問題になりました。

隣接行列を列挙する際は以下に注意します。

対称行列*1
各成分は 0 または 1（単結合のものを考えるため）*2
対角成分は 0（なんでもよいが、便宜上 0 としておく）
各行（および各列）の和は4以下
ノードの番号の対応を入れ替えて行列が一致するものは同じものとみなす
連結なグラフを表現している

最後の二点については、後ほどどのように判定するのか述べます。

隣接行列のビット表現

0 と 1 しか成分がないと聞いて思い浮かぶのは、そう、ビット表現ですね。例えば先ほどの行列であれば、4bit の整数4つで表現できます。

冒頭のリポジトリでは matrix.rs の中で定義している SymmetricBitMatrix<N> がその行列です。

#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
pub struct SymmetricBitMatrix<const N: usize> {
    rows: [u16; N],
}

内部データは16ビット整数 N 個（=行列のサイズ）になっていて、各整数の下 N 桁をデータとして用いる設計にしています。例えば先ほどの隣接行列は以下のような配列に対応します。

$\begin{bmatrix} 1010_{(2)} \\ 1001_{(2)} \\ 1000_{(2)} \\ 0111_{(2)} \end{bmatrix}$
ここで左右が入れ替わっているのは、実装上その方が分かりやすかったからです。 $i$ 行 $j$ 列が、第 $i$ 成分の下から $j$ bit 目に対応する、とすると分かりやすいですね。

次節で実際に列挙の戦略について述べますが、これらもビット演算によって表現することになります。それぞれどういう手続きに対応するかの詳細は省略するので、実際にコードを読んでぜひ確認してみてください。

列挙の戦略

以下の3つのステップに分かれます。

明らかに不要なものを除いて、とにかく隣接行列を列挙する
特徴量を用いてざっくり分類する
同じ特徴量を持つもの同士を比較して、異なる構造のみを残す

順に説明していきます。

明らかに不要なものを除いて、とにかく隣接行列を列挙する

基本は全列挙です。ただ本当に $N \times N$ の行列をすべて列挙すると $2^{N^{2}}$ 個列挙することになり、とてつもない時間がかかることになります*3。そのため、いかにして明らかに不要なものをスキップするかが肝心です。

まず、隣接行列が対称行列で、対角成分が 0 であることを使います。これにより、右上の三角の部分が決まれば、左下の三角の部分は決まることになります。これで $2^{N(N - 1)/2}$ 個まで減りました。

次に、番号を割り当てる際に、ノードに集まっているエッジの本数（次数と言います）が多いもの順（つまり小さい番号はエッジ数が多い）になるようにします。仮に次数が1, 2, 3, 4 のノードがそれぞれ $N/4$ ずつあるとした場合、何も考えない場合の番号の付け方は $N!$ 通り、次数を考慮した場合は $((N/4)!)^{4}$ 通りになるので、ざっくり $((N/4)!)^{4} / N!$ 倍になります。実際に $2^{N(N - 1)/2}$ にこの数を掛けた数字を計算してみると、 $N = 8$ のときで $10^{5}$ 程度となります。 $2^{N (N - 1) / 2}$ が $10^{8}$ ぐらいなので、かなり絞ることができていることが分かります。

search.rsの MatrixSearcher<N> ではさらに、「次数が同じノード同士の場合、それらの二進数表現（つまり SymmetricBitMatrix<N> の成分）が大きい方が番号が小さい」という順序を満たすように列挙しています。具体的には MatrixSearcher<N> のメソッドの is_fine_up_to でこの判定を実装しています。これによりさらに劇的に列挙する数を減らすことに成功しています。

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()
}

実は、どのようなグラフであってもこのような番号付けが可能かどうかは、まだ証明できていません。もし誰か証明できたという方がいればこっそり教えてください*4。下記はその問題を述べた投稿です。

(続き)このとき f:V→{1,2,...,|V|} であって, 以下をすべて満たすものが存在することを示せ.
・fは全単射
・v,w∊V に対し, |v|<|w| ならば f(v)>f(w)
・v,w∊V に対し, |v|=|w| かつ f(v)>f(w) ならば S_v≦S_w. ただし
S_v=Σ2^f(v') (v'はvに隣接するすべてのノードを走る)
とする(S_w も同様).
— すむーずぷりんちゃん🍮@技術書典16オンライン出店予定 (@mat_der_D) May 29, 2024

またその隣接行列が連結なグラフを表現していることを確認する必要があります。これは隣接行列を $A$ として

$\displaystyle\sum_{n=0}^{N - 1} A^{n} \begin{pmatrix}1 \\ 0 \\ \vdots \\ 0 \end{pmatrix}$
のすべての成分が 0 でないことを確かめればよいことがわかります。なぜなら、このベクトルのある成分が 0 でないことは、0番のノードから出発して $(N - 1)$ 本以下のエッジをうまく辿ればそのノードにたどり着けることを意味していますが、これはグラフが連結であることと同値だからです。この判定は SymmetricBitMatrix<N> の is_connected というメソッドで実装しています。

const REPUNIT: u16 = {
    let mut repunit = 0;
    let mut i = 0;
    while i < N {
        repunit |= 1 << i;
        i += 1;
    }
    repunit
};

pub fn is_connected(&self) -> bool {
    let mut connected = 1;
    for _ in 0..N {
        let tmp = connected;
        for i in 0..N {
            if tmp & 1 << i != 0 {
                connected |= self.rows[i];
            }
        }
        if connected == Self::REPUNIT {
            return true;
        }
    }
    false
}

特徴量による分類

隣接行列の候補を列挙する方法が分かったので、あとは比較して同じものを省く・・・としてもよいのですが、これはかなり計算時間がかかります。二つの隣接行列が与えられたとき、なんの前情報もない場合は、単純にすべての番号付けを入れ替えて一致するかどうか確かめる他ありません。これは $O(N!)$ の計算時間がかかるので、できるだけ避けたいです。

そこで、グラフの形状から定まる特徴量を計算する方法を考えます。Morgan のアルゴリズムとよばれる手法を用いると、各ノードに整数値を割り当てることができます。二つのグラフが与えられたとき、このアルゴリズムで割り当てられた数値が異なれば、最初からその二つのノードは対応関係にないことが分かる、という代物です。私のプログラムでは、オリジナルのものからは少し変えた、以下の手法を採用しています。

まず $s = 1, 2, \dots N$ に対して

$\begin{pmatrix}a_{0, s} \\ a_{1, s} \\ \vdots \\ a_{N - 1, s} \end{pmatrix} = A^{s} \begin{pmatrix}1 \\ 1 \\ \vdots \\ 1\end{pmatrix}$
を計算して
$\begin{pmatrix}a_{0,1} & a_{0,2} & \cdots & a_{0, N} \\ a_{1, 1} & a_{1, 2} & \cdots & a_{1, N} \\ \vdots & \vdots & \ddots & \vdots \\ a_{N-1, 1} & a_{N - 1, 2} & \cdots & a_{N-1, N} \end{pmatrix}$
という行列を計算します。この各 $i$ 行の行ベクトルがノード $i$ の特徴量となります。整数値にするため、行ベクトルの配列のハッシュ値をもって特徴量とします。

実装は SymmetricBitMatrix<N> の calc_morgan_hashes で与えています。

fn calc_morgan_hashes(&self) -> [u64; N] {
    let mut hash_array = [[0; N]; N];
    let mut hash = [1u32; N];
    for s in 0..N {
        let mut new_hash = [0; N];
        for (i, new_hash_i) in new_hash.iter_mut().enumerate() {
            for (j, hash_j) in hash.iter().enumerate() {
                if self.rows[i] & 1 << j != 0 {
                    *new_hash_i += hash_j;
                }
                hash_array[i][s] = *new_hash_i;
            }
        }
        hash = new_hash;
    }

    let mut ret = [0; N];
    for (ret_elem, arr) in ret.iter_mut().zip(hash_array.iter()) {
        let mut hasher = FxHasher::default();
        arr.hash(&mut hasher);
        *ret_elem = hasher.finish();
    }
    ret
}

これにより、各ノードに整数値を割り当てることができました。もし二つの隣接行列が同じグラフ（の番号違い）を表現しているのであれば、この特徴量の組み合わせも一致していなければなりません。つまり、特徴量の組が異なれば、比較するまでもなく別の隣接行列と判定できるわけです。

それだけではありません。仮に特徴量の組が同じ二つの隣接行列があったとします。このとき、同じグラフを表すかの判定で並べ替えなければならない数は $N!$ よりずっと少なく済む可能性があります。なぜなら、同じ特徴量を持つノード同士以外はそもそも対応するはずがないからです。

例えば、何も情報がない $N = 4$ の隣接行列同士の比較であれば $4! = 24$ 通りの入れ替えが必要ですが、特徴量が $(4, 4, 8, 8)$ と算出された隣接行列同士の比較であれば、4 同士のノード、8 同士のノードで比較するだけなので、 $2! \times 2! = 4$ 通りだけ入れ替えればよいのです。この機構は SymmetricBitMatrix<N> の partially_canonicalized および search.rs にある RowOrderStore<N> の実装によって実現しています（引用略）。

さらにダメ押しで、固有多項式の $(N - 3)$ 次の係数（から定まるハッシュ値）も計算しています（SymmetricBitMatrix<N> の calc_eig3_hash）。この量は固有値の対称式で表されるので、同じグラフを表す隣接行列同士では一致していなければなりません。

fn calc_eig3_hash(&self) -> u64 {
    let mut eig3 = 0;
    for (i, row_i) in self.rows.iter().enumerate() {
        for j in 0..i {
            if row_i & (1 << j) == 0 {
                continue;
            }
            for row_k in self.rows.iter().take(j) {
                eig3 += (row_k >> i) & (row_k >> j) & 1;
            }
        }
    }

    let mut state = FxHasher::default();
    eig3.hash(&mut state);
    state.finish()
}

以上を組み合わせて、最初から比較せずとも異なるグラフだと分かるものは仕訳けておくようにしています。matrix.rs の MatrixHash<N> は、これらの特徴量を保持するための構造体で、これをキーとした HashMap を利用しています。

#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
pub struct MatrixHash<const N: usize> {
    morgan_hashes: [u64; N],
    canonical_morgan_hashes: [u64; N],
    eig3_hash: u64,
}

同じ特徴量を持つもの同士を比較して、異なる構造のみを残す

ここまでなんとかして比較する量を減らそうと工夫をしてきました。しかし最後はやはり頑張って比較するより他ありません。

先ほど述べた特徴量を用いた仕訳けおよび入れ替え量の制限により、大幅に計算量は落とされていますが、それでもなお計算量を落とせないケースがあります。それは、どのノードの次数も同じ場合です。これについては頑張って $O(N!)$ の計算をして比較するしかありません。そして残念ながら、この部分が計算時間のボトルネックになっています。

ひとつ良いお知らせは、特徴量ごとに仕訳けたおかげで、並列に計算できることです。私の環境では、128 個ぐらいに分割するといい感じに計算時間が短くなりました。

列挙の結果

さて、プログラムを実行すると以下のようになります。

===== # of hydrocarbons with only single-bonds =====
#C ==> # of structures
 2 ==>        1
 3 ==>        2
 4 ==>        6
 5 ==>       21
 6 ==>       78
 7 ==>      353
 8 ==>     1929
 9 ==>    12207
10 ==>    89402

これらの数値は『炭素と水素による構造式』のプログラムの数値と一致しました。やったー！

なお Measure-Command コマンドを使って測った計算時間は以下の通り7.7秒程度となりました。『炭素と水素による構造式』によればこの計算に72時間程度かかったようなので、劇的な進歩なのではないでしょうか。