Kabsch アルゴリズムによる最小二乗重ね合わせと RMSD 算出：数学的導出・歴史的背景・実利的応用

11942 words

60 minutes

Kabsch アルゴリズムによる最小二乗重ね合わせと RMSD 算出：数学的導出・歴史的背景・実利的応用

2026-03-16

Computational Chemistry

Structural Biology

Kabsch Algorithm

RMSD

SVD

Molecular Superposition

last_modified: 2026-03-16

生成AIによる自動生成記事に関する免責事項: 本記事は、Kabsch アルゴリズムおよびその関連理論に関する一次文献・教科書・実装を基に、生成AI（Claude, Anthropic）によって自動構成された技術解説記事です。内容の正確性を保持するよう努めていますが、厳密な理論的証明や数値データの解釈については、必ず原著論文を参照してください。本記事に含まれる数式・コード・参考文献の引用は、学術的な解説・教育目的のものであり、特定の研究成果の保証を意図するものではありません。

1. 序論：分子構造比較の問題設定と RMSD の中心的役割#

計算化学・構造生物学・創薬科学において、二つの分子構造が「どの程度一致しているか」を定量化することは、様々な下流タスク――タンパク質立体構造予測の精度評価、分子動力学（MD）シミュレーションにおける構造的安定性のモニタリング、量子化学計算における重複構造の除去、ならびにドラッグデザインにおける活性コンフォメーションの照合――の根幹をなす基本操作である。

この定量化の主要な指標として、原子間距離の二乗平均平方根（Root Mean Square Deviation; RMSD）が広く用いられている。RMSD は直観的な解釈が可能であり、単位が座標と同一（Å など）であるため、分子のスケールと直接対応する物理量として認識されやすい。しかしながら、二つの分子構造を比較するためには、それらを最適に重ね合わせた（superpose した）後のRMSDを算出しなければならない。なぜなら、剛体的な並進・回転によって RMSD は大きく変化するからである。

「最適な重ね合わせ」とは、並進と回転の組み合わせの中で RMSD を最小化するものを意味する。この問題は、最適剛体重ね合わせ問題（Optimal Rigid-Body Superposition Problem）と呼ばれ、行列解析の観点からは直交 Procrustes 問題（Orthogonal Procrustes Problem）の一例として定式化される。

本稿が主題とするのは、この最適化問題に対して Wolfgang Kabsch が 1976 年および 1978 年に提案した、特異値分解（Singular Value Decomposition; SVD）に基づく解法——一般に「Kabsch アルゴリズム」と呼称される手法——の理論的基礎・歴史的背景・数値実装上の諸問題、および実際の計算化学プログラムにおける具体的な利用形態についてである。特に、以下のコードに示される _kabsch_rmsd メソッドが実装する考え方を中心に据えながら、その数学的根拠を丁寧に展開する。

@staticmethod
def _kabsch_rmsd(pa: np.ndarray, pb: np.ndarray) -> float:
    """Minimum RMSD between pa and pb over all proper rotations.

    Chirality guard (d < 0 → return inf):
        When the optimal superposition requires an improper rotation
        (det < 0), the classical Kabsch algorithm still returns a finite
        RMSD after sign correction, potentially mis-identifying enantiomers
        for flat chiral molecules.
        Fix: return ``inf`` immediately when d < 0. The decision is based
        on *state* (whether a reflection is required), not on numerical
        magnitude, so it is independent of molecular shape. This has no
        effect on non-chiral molecules (H₂O, meso compounds) because their
        mirror images can be superimposed with a true rotation (d = +1) for
        at least one candidate.
    """
    if len(pa) == 0:
        return 0.0
    U, S, Vt = np.linalg.svd(pb.T @ pa)
    d = np.linalg.det(U) * np.linalg.det(Vt)

    if d < 0:
        return float("inf")

    E0 = np.sum(pa ** 2) + np.sum(pb ** 2)
    rmsd_sq = max(0.0, E0 - 2.0 * np.sum(S)) / len(pa)
    return float(np.sqrt(rmsd_sq))

2. 問題の形式的定義：最適剛体重ね合わせと RMSD#

2.1 RMSD の数学的定義#

$N$ 個の対応する原子対を持つ二つの点集合 $\mathbf{P} = \{\mathbf{p}_i\}_{i=1}^{N}$ および $\mathbf{Q} = \{\mathbf{q}_i\}_{i=1}^{N}$ （ $\mathbf{p}_i, \mathbf{q}_i \in \mathbb{R}^3$ ）が与えられたとき、RMSD は次のように定義される。

$\mathrm{RMSD}(\mathbf{P}, \mathbf{Q}) = \sqrt{ \frac{1}{N} \sum_{i=1}^{N} \| \mathbf{p}_i - \mathbf{q}_i \|^2 }$

この定義は各原子対間の三次元ユークリッド距離の二乗の平均値の平方根であり、構造的類似性の直感的な尺度となる。ただし、上記の式は原子の対応関係が既知であり、かつ両構造が適切な座標系に配置されていることを前提とする。

2.2 最適重ね合わせ問題の定式化#

実際の構造比較においては、まず両構造の重心を合わせる並進操作を施した後、最適な回転行列 $\mathbf{R} \in SO(3)$ （行列式 $+1$ の直交行列の集合）を求める問題を解く。

並進の除去: 両構造の重心 $\bar{\mathbf{p}} = \frac{1}{N}\sum_{i=1}^{N}\mathbf{p}_i$ 、 $\bar{\mathbf{q}} = \frac{1}{N}\sum_{i=1}^{N}\mathbf{q}_i$ を減算して重心中心化を行う。

$\mathbf{p}'_i = \mathbf{p}_i - \bar{\mathbf{p}}, \quad \mathbf{q}'_i = \mathbf{q}_i - \bar{\mathbf{q}}$

最適回転の探索: 重心中心化後の点集合に対して、最小 RMSD を実現する回転行列を求める。

$\mathbf{R}^* = \arg\min_{\mathbf{R} \in SO(3)} \sum_{i=1}^{N} \| \mathbf{p}'_i - \mathbf{R} \mathbf{q}'_i \|^2$

この問題は直交 Procrustes 問題の特殊ケースであり、行列解析・数値最適化の文脈で古くから研究されてきた。

2.3 最適化目的関数の展開#

目的関数を展開すると、その本質的な構造が明らかになる。

$\sum_{i=1}^{N} \| \mathbf{p}'_i - \mathbf{R} \mathbf{q}'_i \|^2 = \sum_{i=1}^{N} \left( \|\mathbf{p}'_i\|^2 + \|\mathbf{q}'_i\|^2 - 2 \mathbf{p}'^{\top}_i \mathbf{R} \mathbf{q}'_i \right)$

回転行列 $\mathbf{R}$ は $\| \mathbf{R} \mathbf{q}'_i \| = \| \mathbf{q}'_i \|$ を保つ（等長変換）ため、上式の第一・第二項は $\mathbf{R}$ に依存しない。したがって、目的関数の最小化は、交差項 $\sum_{i=1}^{N} \mathbf{p}'^{\top}_i \mathbf{R} \mathbf{q}'_i$ の最大化と等価になる。

これをまとめると、次の等価な問題として定式化できる。

$\mathbf{R}^* = \arg\max_{\mathbf{R} \in SO(3)} \mathrm{tr} \left( \mathbf{R} \mathbf{H} \right)$

ただし、 $\mathbf{H}$ は共分散行列（クロス共分散行列）と呼ばれる $3 \times 3$ 行列で、次式で定義される。

$\mathbf{H} = \sum_{i=1}^{N} \mathbf{q}'^{\top}_i \mathbf{p}'_i = \mathbf{Q}'^{\top} \mathbf{P}'$

ここで $\mathbf{P}'$ および $\mathbf{Q}'$ はそれぞれ中心化された座標を行方向に並べた $N \times 3$ 行列である。このように、最適重ね合わせ問題の核心は、与えられた $3 \times 3$ 行列に対してトレースを最大化する $SO(3)$ 元を求める問題に帰着する。

3. 歴史的背景：Kabsch アルゴリズムの成立過程#

3.1 前史：Wahba 問題と直交 Procrustes 問題#

最適重ね合わせ問題の起源は、計算化学よりも古く、統計学・測量学の問題に端を発する。1965 年、Grace Wahba は人工衛星の姿勢決定（attitude determination）という工学的文脈において、雑音を含む方向ベクトルの観測値から最適な回転行列を推定する問題を定式化した。これは「Wahba 問題」として知られ、最適制御・ロボティクスの分野においても広く参照される古典的問題である。

一方、行列解析の分野では、直交 Procrustes 問題が 1966 年に Schönemann によって SVD を用いて解かれた。Procrustes 問題（プロクルステス問題）という名称は、ギリシャ神話に登場する宿の主人プロクルステス（Procrustês）に由来する――彼は旅人を鉄製の寝台に強制的に合わせたとされる。この名称は、一方のデータ集合を他方に最適に「合わせる」という操作の比喩として採用されており、現在は多変量統計・形状解析において標準的な用語となっている。

$\min_{\mathbf{R}: \mathbf{R}^{\top}\mathbf{R} = \mathbf{I}} \| \mathbf{A} - \mathbf{B}\mathbf{R} \|_F^2$

Schönemann の 1966 年の結果は、この問題の解が $\mathbf{B}^{\top}\mathbf{A}$ の SVD から直接得られることを示した。しかし、この時点では行列式の符号（ $\det(\mathbf{R}) = +1$ vs $-1$ ）の制約、すなわち $SO(3)$ （真の回転）と $O(3)$ （回転＋反転）の区別が必ずしも明示的に扱われていなかった。

3.2 Kabsch アルゴリズムの提案（1976–1978）#

化学結晶学者 Wolfgang Kabsch は、1976 年に Acta Crystallographica Section A に発表した論文「A solution for the best rotation to relate two sets of vectors」において、SVD を用いた最適回転の導出を結晶化学の文脈で独立に提案し、反転操作（improper rotation）の処理を明示的に組み込んだ。これが今日「Kabsch アルゴリズム」として広く参照される手法の原形である。

1978 年には同じく Acta Crystallographica Section A に続報「A discussion of the solution for the best rotation to relate two sets of vectors」を発表し、数値安定性の側面や特殊ケース（共線性など）の処理についての補足的な議論を行った。これにより、実装レベルでの利用可能性が大きく向上した。

Kabsch の貢献の中心は次の二点にある。

SVD による最適回転の構成的導出: 共分散行列 $\mathbf{H}$ の SVD $\mathbf{H} = \mathbf{U} \boldsymbol{\Sigma} \mathbf{V}^{\top}$ に基づく最適回転行列 $\mathbf{R} = \mathbf{V}\mathbf{U}^{\top}$ の明示的な表現の提示。
行列式条件による improper rotation の排除: $\det(\mathbf{V}\mathbf{U}^{\top}) = -1$ となる場合に、最小特異値に対応する符号を反転させることで真の回転（proper rotation）を保証する手続きの組み込み。

これらの寄与は、当時まだ計算機利用が一般化しつつあった結晶構造解析の実務において、数値的に安定かつ解釈の明確なアルゴリズムを提供したという点で大きな意義を持った。

3.3 並行する研究：McLachlan 法・Diamond 法#

Kabsch と同時代あるいはその後に、同種の問題に対する異なる解法が複数提案された。

McLachlan 法（1979, 1982）: Andrew D. McLachlan は、タンパク質構造比較の文脈において、共分散行列の固有値分解（eigendecomposition）に基づく最適回転行列の構成法を提案した。行列 $\mathbf{H}^{\top}\mathbf{H}$ の固有分解を利用するこの手法は、Kabsch 法と数値的に等価であるが、実装上の観点から一定の違いを持つ。特に McLachlan は、タンパク質の $C_\alpha$ 炭素のみを用いた骨格構造の比較に本手法を適用し、タンパク質フォールドの体系的分類における RMSD の有用性を示した。

Diamond 法（1988）: Robert Diamond は、四元数（quaternion）表現を用いて最適回転を記述する手法を Acta Crystallographica Section A に発表した。四元数による表現は、 $SO(3)$ の三重被覆（2対1射影）の観点から幾何学的に明快であり、計算機グラフィックス・ロボティクスの分野において好まれる定式化である。この方法では、 $4 \times 4$ 対称行列 $\mathbf{K}$ の最大固有値・固有ベクトルを計算することで最適回転が得られる。数値計算上の観点からは、小規模システムに対して特に有利であるとされるが、SVD ベースの Kabsch 法と数学的に等価である。

3.4 1990 年代以降の定着と標準化#

1980 年代後半から 1990 年代にかけて、タンパク質構造データバンク（Protein Data Bank; PDB）が整備され、分子動力学シミュレーションソフトウェア（AMBER, CHARMM, GROMACS 等）が広く普及するとともに、RMSD に基づく構造比較は計算生物学の標準的な操作として確立した。この過程で Kabsch アルゴリズムは、その明快な SVD による定式化と実装の容易さから、最も広く利用される手法の一つとなった。

現代では、Python の numpy.linalg.svd、あるいは scipy.spatial.transform.Rotation などのライブラリにより、数行のコードで Kabsch アルゴリズムを実装できる。後述する _kabsch_rmsd の実装は、この流れを受け継ぎ、NumPy の SVD 関数を直接利用した簡潔な実装例となっている。

4. 数学的基礎：特異値分解（SVD）の理論#

4.1 SVD の定義と存在定理#

任意の実行列 $\mathbf{A} \in \mathbb{R}^{m \times n}$ （ $m \geq n$ を仮定）に対して、特異値分解は次の形式で一意的（ただし特異値の順序と符号の扱いを除く）に存在する。

$\mathbf{A} = \mathbf{U} \boldsymbol{\Sigma} \mathbf{V}^{\top}$

ここで、

$\mathbf{U} \in \mathbb{R}^{m \times m}$ は直交行列（ $\mathbf{U}^{\top}\mathbf{U} = \mathbf{I}_m$ ）
$\boldsymbol{\Sigma} \in \mathbb{R}^{m \times n}$ は対角行列で、対角成分 $\sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_n \geq 0$ を特異値と呼ぶ
$\mathbf{V} \in \mathbb{R}^{n \times n}$ は直交行列（ $\mathbf{V}^{\top}\mathbf{V} = \mathbf{I}_n$ ）

SVD の存在は、実対称行列のスペクトル定理（固有値分解の存在）から導くことができる。具体的には、 $\mathbf{A}^{\top}\mathbf{A}$ は半正定値実対称行列であるから、固有値分解 $\mathbf{A}^{\top}\mathbf{A} = \mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{\top}$ （ $\boldsymbol{\Lambda}$ は非負固有値の対角行列）が存在する。 $\sigma_i = \sqrt{\lambda_i}$ とおくと、 $\mathbf{U}$ の各列ベクトル（左特異ベクトル）が構成でき、SVD が得られる。

4.2 フロベニウスノルム最適化との関係#

直交 Procrustes 問題は、フロベニウスノルムを用いた行列最適化問題として書き直せる。

$\min_{\mathbf{R} \in O(n)} \| \mathbf{A} - \mathbf{B}\mathbf{R} \|_F^2$

フロベニウスノルムの展開から

$\| \mathbf{A} - \mathbf{B}\mathbf{R} \|_F^2 = \mathrm{tr}(\mathbf{A}^{\top}\mathbf{A}) + \mathrm{tr}(\mathbf{R}^{\top}\mathbf{B}^{\top}\mathbf{B}\mathbf{R}) - 2\mathrm{tr}(\mathbf{B}^{\top}\mathbf{A}\mathbf{R}^{\top})$

$\mathbf{R}$ が直交行列であるときに $\mathrm{tr}(\mathbf{R}^{\top}\mathbf{B}^{\top}\mathbf{B}\mathbf{R}) = \mathrm{tr}(\mathbf{B}^{\top}\mathbf{B})$ が成立するため（トレースの循環性： $\mathrm{tr}(\mathbf{ABC}) = \mathrm{tr}(\mathbf{BCA})$ ）、問題は次の最大化問題に等価になる。

$\max_{\mathbf{R} \in O(n)} \mathrm{tr}(\mathbf{B}^{\top}\mathbf{A}\mathbf{R}^{\top})$

$\mathbf{M} = \mathbf{B}^{\top}\mathbf{A}$ と記し、SVD $\mathbf{M} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{\top}$ を用いると、

$\mathrm{tr}(\mathbf{M}\mathbf{R}^{\top}) = \mathrm{tr}(\mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{\top}\mathbf{R}^{\top}) = \mathrm{tr}(\boldsymbol{\Sigma}\mathbf{V}^{\top}\mathbf{R}^{\top}\mathbf{U})$

$\mathbf{Z} = \mathbf{V}^{\top}\mathbf{R}^{\top}\mathbf{U}$ と置くと、 $\mathbf{R}$ が直交行列なら $\mathbf{Z}$ もまた直交行列である。したがって

$\mathrm{tr}(\boldsymbol{\Sigma}\mathbf{Z}) = \sum_i \sigma_i Z_{ii}$

直交行列の各対角成分は $|Z_{ii}| \leq 1$ を満たす（これは直交行列の行ノルム・列ノルムが 1 であることから従う）。特異値 $\sigma_i \geq 0$ を仮定すれば、 $\sum_i \sigma_i Z_{ii}$ は $Z_{ii} = 1$ のとき最大値 $\sum_i \sigma_i$ を取る。すなわち $\mathbf{Z} = \mathbf{I}$ （単位行列）のときが最適であり、このとき

$\mathbf{V}^{\top}\mathbf{R}^{\top}\mathbf{U} = \mathbf{I} \implies \mathbf{R} = \mathbf{V}\mathbf{U}^{\top}$

以上が、SVD に基づく最適回転行列の導出である。ただし、 $O(n)$ の制約（直交行列全体）の下では $\mathbf{R} = \mathbf{V}\mathbf{U}^{\top}$ が最適解であるが、 $SO(3)$ （真の回転のみ）の制約を課す場合には $\det(\mathbf{R}) = +1$ の条件が追加される。

4.3 行列式条件と improper rotation の処理#

$\det(\mathbf{V}\mathbf{U}^{\top}) = \det(\mathbf{V})\det(\mathbf{U}^{\top}) = \det(\mathbf{V})\det(\mathbf{U})$ であり、直交行列の行列式は $\pm 1$ なので $\det(\mathbf{R}) = \pm 1$ となる。

$\det(\mathbf{R}) = +1$ の場合： $\mathbf{R}$ は真の回転（proper rotation）であり、 $SO(3)$ に属する。
$\det(\mathbf{R}) = -1$ の場合： $\mathbf{R}$ は反転を含む不正則変換（improper rotation）であり、 $O(3) \setminus SO(3)$ に属する。

分子構造の比較において、鏡像異性体（enantiomer）を同一視しないためには、 $SO(3)$ への制約が不可欠である。 $\det(\mathbf{R}) = -1$ となる場合の古典的な対処法（Kabsch 1978 の提案）は、最小特異値に対応する $\mathbf{V}$ の列ベクトルの符号を反転させて行列式を $+1$ に補正することである。すなわち

$\mathbf{R}^* = \mathbf{V} \begin{pmatrix} 1 & & \\ & 1 & \\ & & \det(\mathbf{V}\mathbf{U}^{\top}) \end{pmatrix} \mathbf{U}^{\top}$

この操作は、最小特異値 $\sigma_3$ を実効的に $-\sigma_3$ として扱うことに相当し、 $\mathrm{tr}(\boldsymbol{\Sigma}\mathbf{Z})$ の値はわずかに低下する（ $\sigma_3$ 分だけ RMSD が増加する）。これは improper rotation を許容した場合の RMSD よりも真の回転による RMSD が大きいことを意味し、鏡像関係にある構造は有限の RMSD 距離を持つという数学的保証となる。

5. Kabsch アルゴリズムの詳細な導出#

5.1 アルゴリズムの手順#

以下に、Kabsch アルゴリズムの完全な手順を示す。

入力: 中心化済み座標行列 $\mathbf{P}' \in \mathbb{R}^{N \times 3}$ （参照構造）および $\mathbf{Q}' \in \mathbb{R}^{N \times 3}$ （移動構造）

Step 1: 共分散行列の計算

$\mathbf{H} = (\mathbf{Q}')^{\top} \mathbf{P}' \in \mathbb{R}^{3 \times 3}$

上記の実装では pb.T @ pa に相当する。pa が $\mathbf{P}'$ （参照構造）、pb が $\mathbf{Q}'$ （移動構造）に対応する。

Step 2: $\mathbf{H}$ の SVD の計算

$\mathbf{H} = \mathbf{U} \boldsymbol{\Sigma} \mathbf{V}^{\top}$

np.linalg.svd(pb.T @ pa) は $\mathbf{U}$ （左特異ベクトル行列）、 $\mathbf{S}$ （特異値ベクトル）、 $\mathbf{V}^{\top}$ （右特異ベクトル行列の転置）を返す。

Step 3: 行列式の符号確認

$d = \det(\mathbf{U}) \cdot \det(\mathbf{V}^{\top})$

d = np.linalg.det(U) * np.linalg.det(Vt)

$d = +1$ であれば最適な重ね合わせは真の回転であり、通常の手順を続ける。 $d = -1$ であれば、improper rotation が必要であり、実装に応じて対処する。

Step 4（ $d = +1$ の場合）: RMSD の計算

$\mathrm{RMSD}^2 = \frac{1}{N} \max\left(0,\ E_0 - 2\sum_{i=1}^{3}\sigma_i\right)$

ここで $E_0 = \|\mathbf{P}'\|_F^2 + \|\mathbf{Q}'\|_F^2 = \sum_i \|\mathbf{p}'_i\|^2 + \sum_i \|\mathbf{q}'_i\|^2$ はノルムの二乗和である。

5.2 重心中心化の必要性#

Kabsch アルゴリズムを適用する前に、両構造の重心中心化が必要である。並進は RMSD の最小化問題において独立に解け、その最適解が重心を一致させることである。

$\bar{\mathbf{p}} = \frac{1}{N}\sum_{i=1}^N \mathbf{p}_i$ 、 $\bar{\mathbf{q}} = \frac{1}{N}\sum_{i=1}^N \mathbf{q}_i$ として、

$\sum_{i=1}^N \| \mathbf{p}_i - \mathbf{R}\mathbf{q}_i - \mathbf{t} \|^2$

の $\mathbf{t}$ についての最小化は $\mathbf{t}^* = \bar{\mathbf{p}} - \mathbf{R}\bar{\mathbf{q}}$ を与える。これを代入すると、残差は中心化座標のみに依存する形になり、回転の最適化が重心中心化後の座標のみを用いて独立に行えることが示される。

コード中の _kabsch_rmsd は、呼び出し前に座標が既に中心化済みであることを前提としており、実際に compute_rmsd メソッドにおいて

ca = coords_a - coords_a.mean(axis=0)
cb = coords_b - coords_b.mean(axis=0)

として重心中心化が行われている。

6. RMSD 計算式の解析的導出#

6.1 二乗 RMSD のトレース表現#

重心中心化座標 $\mathbf{P}', \mathbf{Q}'$ と回転行列 $\mathbf{R}$ を用いた二乗 RMSD を展開する。

$N \cdot \mathrm{RMSD}^2 = \| \mathbf{P}' - \mathbf{Q}' \mathbf{R}^{\top} \|_F^2$

$= \mathrm{tr}\left[(\mathbf{P}' - \mathbf{Q}'\mathbf{R}^{\top})^{\top}(\mathbf{P}' - \mathbf{Q}'\mathbf{R}^{\top})\right]$

$= \mathrm{tr}\left[\mathbf{P}'^{\top}\mathbf{P}' - \mathbf{R}\mathbf{Q}'^{\top}\mathbf{P}' - \mathbf{P}'^{\top}\mathbf{Q}'\mathbf{R}^{\top} + \mathbf{R}\mathbf{Q}'^{\top}\mathbf{Q}'\mathbf{R}^{\top}\right]$

トレースの線形性と直交行列の性質 $\mathrm{tr}(\mathbf{R}\mathbf{A}\mathbf{R}^{\top}) = \mathrm{tr}(\mathbf{A})$ を用いると、

$N \cdot \mathrm{RMSD}^2 = \mathrm{tr}(\mathbf{P}'^{\top}\mathbf{P}') + \mathrm{tr}(\mathbf{Q}'^{\top}\mathbf{Q}') - 2\mathrm{tr}(\mathbf{R}\mathbf{Q}'^{\top}\mathbf{P}')$

$= \|\mathbf{P}'\|_F^2 + \|\mathbf{Q}'\|_F^2 - 2\mathrm{tr}(\mathbf{R}\mathbf{H})$

ここで $\mathbf{H} = \mathbf{Q}'^{\top}\mathbf{P}'$ である。

6.2 特異値による RMSD の表現#

前節の結果から、二乗 RMSD の最小化は $\mathrm{tr}(\mathbf{R}\mathbf{H})$ の最大化に等価である。 $\mathbf{H} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{\top}$ の SVD を用い、最適回転 $\mathbf{R}^* = \mathbf{V}\mathbf{U}^{\top}$ （ $d = +1$ の場合）を代入すると、

$\mathrm{tr}(\mathbf{R}^* \mathbf{H}) = \mathrm{tr}(\mathbf{V}\mathbf{U}^{\top}\mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{\top}) = \mathrm{tr}(\mathbf{V}\boldsymbol{\Sigma}\mathbf{V}^{\top}) = \mathrm{tr}(\boldsymbol{\Sigma}) = \sum_{k=1}^{3} \sigma_k$

したがって、最小二乗 RMSD は次の閉形式で表される。

$\mathrm{RMSD}_{\min}^2 = \frac{1}{N}\left(\|\mathbf{P}'\|_F^2 + \|\mathbf{Q}'\|_F^2 - 2\sum_{k=1}^{3}\sigma_k\right)$

$\mathrm{RMSD}_{\min} = \sqrt{\frac{1}{N}\max\left(0,\ E_0 - 2\sum_{k=1}^{3}\sigma_k\right)}$

この式は、_kabsch_rmsd の実装に直接対応している。

E0 = np.sum(pa ** 2) + np.sum(pb ** 2)
rmsd_sq = max(0.0, E0 - 2.0 * np.sum(S)) / len(pa)
return float(np.sqrt(rmsd_sq))

np.sum(pa ** 2) は $\|\mathbf{P}'\|_F^2$ （中心化済み参照構造の全原子の二乗距離の和）
np.sum(pb ** 2) は $\|\mathbf{Q}'\|_F^2$ （中心化済み移動構造の同様の量）
np.sum(S) は特異値の和 $\sum_{k=1}^{3}\sigma_k$
max(0.0, ...) は数値誤差による負値の防止

6.3 `max(0, ...)` の必要性：数値安定性の観点#

理論的には $E_0 - 2\sum_k \sigma_k \geq 0$ が成立する。これは、 $\mathrm{RMSD}^2$ が非負の量であることから自明であるが、有限精度浮動小数点演算においては丸め誤差により極めて小さな負値が生じることがある。特に、両構造が完全一致する場合（ $\mathrm{RMSD} = 0$ ）に最も顕著であり、 $\sqrt{(-\epsilon)}$ （ $\epsilon \approx 10^{-15}$ ）のような NaN 生成を防ぐために max(0.0, ...) の保護が実装上必要となる。

6.4 `d * S[2]` を使わない理由#

実装コメントに以下の記述がある。

# Using `d * S[2]` (where d is a float ≈ 1.0 ± 1e-15) is technically
# incorrect and introduces unnecessary floating-point bias; np.sum(S)
# is exact.

一部の古い実装では、improper rotation 補正として最小特異値の項に符号 $d$ を掛けた式

$\mathrm{RMSD}^2 \sim \frac{1}{N}\left(E_0 - 2(\sigma_1 + \sigma_2 + d \cdot \sigma_3)\right)$

を用いる。しかし本実装では、 $d < 0$ の場合に inf を返す（後述）ため、この補正が不要となっている。 $d = +1$ が保証された状態では np.sum(S) を直接使用するのが正確であり、 $d$ は厳密には整数値 $\pm 1$ ではなく np.linalg.det が返す浮動小数点値であるため、余分な数値誤差を回避する観点からも np.sum(S) の使用が適切である。

7. キラリティ処理：improper rotation の検出と inf 返却#

7.1 キラリティの問題提起#

鏡像異性体（enantiomers）は、互いに重ね合わせることのできない鏡像関係にある分子である。Kabsch アルゴリズムの古典的な実装（Kabsch 1978）では、 $\det(\mathbf{R}) = -1$ となる場合に符号補正を行って有限の RMSD を返す。しかし、この有限の RMSD は「真の回転による最善の重ね合わせ」のコストではなく、「最も RMSD を最小化する improper rotation」のコストに補正操作を加えたものであり、鏡像関係の二構造に対して意図せず小さな RMSD を返す可能性がある。

特に、ほぼ平面的な分子（平面キラル分子）に対してこの問題は顕在化しやすい。平面分子では $z$ 軸方向の原子変位が小さく、鏡像操作（ $xy$ 平面を鏡面とする反転）が真の回転と構造的に区別しにくい状況が生じるためである。

7.2 本実装における対処： $d < 0$ ならば `inf` 返却#

本実装の _kabsch_rmsd は、 $d = \det(\mathbf{U}) \cdot \det(\mathbf{V}^{\top}) < 0$ のときに即座に float("inf") を返す。

if d < 0:
    return float("inf")

この設計判断の根拠をコメントは次のように述べている。

最適な重ね合わせが improper rotation（det < -1、鏡映）を要する場合、古典的な Kabsch アルゴリズムは符号補正後に有限の RMSD を返すため、平坦なキラル分子の鏡像異性体を同一と誤判定する可能性がある。修正として、 $d < 0$ のとき即座に inf を返す。この判断は（数値的な大きさではなく）「反転が必要か否か」という状態に基づくため、分子の形状に依存しない。

この設計の重要な利点は次の二点にある。

1. 形状非依存性: 従来の補正操作では、最小特異値 $\sigma_3$ の大きさによって RMSD の補正量が変化するため、分子形状（平面性・球対称性）によって挙動が異なる。inf 返却は「真の回転では重ね合わせられない」という状態を数値的大きさに依存せず明示的に宣言する。

2. 非キラル分子への無影響性: $\mathrm{H}_2\mathrm{O}$ やメソ化合物などのアキラルな分子に対しては、少なくとも一つの候補回転において $d = +1$ が成立するため、inf が返されることはない。キラル処理は必要な場合にのみ作用する。

7.3 エナンチオマー判定の文脈#

より上位の _try_candidates メソッドでは、複数の候補回転行列 $\mathbf{R}$ を試行する。

def _try_candidates(self, candidates, ca, cb, groups_a, groups_b):
    min_rmsd = float("inf")
    for R in candidates:
        cb_rot = cb @ R.T
        perm = self._optimal_mapping(ca, cb_rot, groups_a, groups_b)
        if perm is None:
            continue
        rmsd = self._kabsch_rmsd(ca, cb_rot[perm])
        if rmsd < min_rmsd:
            min_rmsd = rmsd
            if min_rmsd < self.rmsd_threshold:
                return min_rmsd
    return min_rmsd

候補回転行列のリストには真の回転（ $\det = +1$ ）のみが含まれており、_kabsch_rmsd 内では「候補回転適用後の構造 cb_rot と参照構造 ca の間にさらに improper rotation が必要かどうか」が検査される。候補回転が適切に選ばれている場合、エナンチオマーに対するすべての候補回転は $d = -1$ を返し、結果として min_rmsd = inf が維持される。これにより、エナンチオマーはいかなる条件下でも「同一構造」と判定されない。

8. 実装上の考察：数値計算の側面#

8.1 NumPy による実装#

以下に、_kabsch_rmsd の全コードを再掲する。

@staticmethod
def _kabsch_rmsd(pa: np.ndarray, pb: np.ndarray) -> float:
    """Minimum RMSD between pa and pb over all proper rotations."""
    if len(pa) == 0:
        return 0.0
    U, S, Vt = np.linalg.svd(pb.T @ pa)
    d = np.linalg.det(U) * np.linalg.det(Vt)

    if d < 0:
        return float("inf")

    E0 = np.sum(pa ** 2) + np.sum(pb ** 2)
    rmsd_sq = max(0.0, E0 - 2.0 * np.sum(S)) / len(pa)
    return float(np.sqrt(rmsd_sq))

pb.T @ pa: 共分散行列 $\mathbf{H} = \mathbf{Q}'^{\top}\mathbf{P}'$ の計算（ $3 \times N$ 行列と $N \times 3$ 行列の積、 $O(3 \times N \times 3) = O(N)$ ）
np.linalg.svd: $3 \times 3$ 行列に対する SVD（サイズ固定のため定数時間操作 $O(1)$ ）
np.linalg.det(U) * np.linalg.det(Vt): 行列式の積（各 $O(3^3) = O(1)$ ）
np.sum(pa ** 2) + np.sum(pb ** 2): フロベニウスノルムの二乗和（ $O(N)$ ）

全体の計算複雑度は $O(N)$ であり、原子数 $N$ に対して線形スケールする。

8.2 `np.linalg.svd` の内部実装と数値安定性#

numpy.linalg.svd は LAPACK の dgesdd（分割統治アルゴリズムに基づく SVD ルーチン）または dgesvd（ゴルブ-ラインシュ二重対角化アルゴリズム）を内部的に呼び出す。 $3 \times 3$ という小行列に対しては数値安定性の問題はほぼ生じないが、一般の $m \times n$ 行列に対して SVD の数値的安定性は QR 分解と並んで最良のクラスに属する（条件数の増幅が最小限に抑えられる）。

行列式の計算においては、np.linalg.det が直接用いられているが、 $3 \times 3$ 行列に対しては解析式で計算される。浮動小数点誤差の観点から $d$ は $+1$ または $-1$ に近い値を取るが、 $d = 0.9999...$ が $d < 0$ と誤判定されることはない。

8.3 `full_matrices=False` による計算省略#

本実装では np.linalg.svd(pb.T @ pa) をデフォルト引数で呼び出している。pb.T @ pa は $3 \times 3$ 行列であるため、完全 SVD（full_matrices=True）と圧縮 SVD（full_matrices=False）は等価であり、full_matrices=False による計算省略は不要である。 $N \gg 3$ となる場合（大規模タンパク質の比較など）には、 $N \times N$ の補完ベクトルの計算を回避するために full_matrices=False の指定が有効であるが、共分散行列は常に $3 \times 3$ であるため本実装では影響しない。

8.4 原子マッピングとの組み合わせ#

_kabsch_rmsd は、呼び出し元の _try_candidates において、ハンガリアンアルゴリズムによる最適原子マッピング（_optimal_mapping）の後に適用される。

for R in candidates:
    cb_rot = cb @ R.T
    perm = self._optimal_mapping(ca, cb_rot, groups_a, groups_b)
    if perm is None:
        continue
    rmsd = self._kabsch_rmsd(ca, cb_rot[perm])

このため _kabsch_rmsd に渡される pa（= ca）と pb（= cb_rot[perm]）は原子順序の対応が既に確立されており、Kabsch アルゴリズムは純粋に「固定された原子対応に対する最小 RMSD」を計算する役割を担う。

全体のアルゴリズムは「粗い回転探索 → ハンガリアン法による原子対応 → Kabsch 法による微細な最小化」の三段階となっており、PCA 整合後に候補回転を試行するという設計は、組み合わせ爆発を回避しつつ大域的最適解に近い結果を得る実用的なアプローチである。

9. 関連手法との比較#

9.1 四元数法（Diamond 1988, Horn 1987）#

四元数 $q = (q_0, q_1, q_2, q_3)$ （ $|q|=1$ ）による $SO(3)$ の表現を用いる手法では、 $4 \times 4$ 実対称行列 $\mathbf{K}$ の最大固有値・固有ベクトルを求めることで最適回転が得られる。 $\mathbf{K}$ の行列要素は共分散行列 $\mathbf{H}$ の成分から構成される。