Matrix 行列式と逆行列
概要
$N \times N$ 正方行列 $A$ に対する行列式 $\det A$ と逆行列 $A^{-1}$ は、線形代数の中で最も「数学的な定義」と「計算手順」が乖離する対象である。 数学的には余因子展開と adjugate 行列で定義されるが、その通り実装すると $O(N!)$ になる。 実用では LU 分解を経由するのが標準で、$O(N^3)$ で済む。
本記事では sangi の det と inverse が
サイズ・要素型・利用文脈に応じてどのアルゴリズムを選ぶかを解説する。
分解そのものの詳細は Matrix API リファレンス および
個別の分解解説記事 (LU / QR / Cholesky / SVD) を参照。
余因子展開 (小行列限定)
余因子展開は行列式の定義そのものである:
$$\det A = \sum_{j=0}^{N-1} (-1)^{i+j} a_{ij} \det M_{ij}$$
($M_{ij}$ は $i$ 行 $j$ 列を除いた $(N-1) \times (N-1)$ 部分行列)
素直に再帰すると計算量は $O(N!)$ で爆発する。 sangi は $N \leq 4$ に限って閉形式・展開済みコード経路を提供する:
| $N$ | 計算量 | 備考 |
|---|---|---|
| 1 | $O(1)$ | $a_{00}$ そのまま |
| 2 | $O(1)$ | $a_{00} a_{11} - a_{01} a_{10}$ |
| 3 | $O(1)$ | Sarrus 規則 (6 項) |
| 4 | $O(1)$ | $2 \times 2$ ブロックの余因子による高速展開 |
| $\geq 5$ | LU 経路に切替 | $O(N^3)$ |
$N \in \{2, 3, 4\}$ は SE(3) 変換 ($4 \times 4$ ホモジニアス) や $3 \times 3$ 回転で頻出するため、 専用展開で完全アンロールしている。 コンパイラの定数畳み込みが効き、レジスタ内で完結する。
関連記事: 余因子展開 / Cramer の公式 / 行列式の導出
LU 分解による det
$N \geq 5$ では LU 分解 (部分ピボッティング付き) を経由する。 $PA = LU$ に分解できれば:
$$\det A = \det(P^{-1}) \cdot \det L \cdot \det U = (-1)^s \cdot \prod_{i=0}^{N-1} u_{ii}$$
($s$ はピボット行交換の回数。$\det L = 1$ は対角が 1 の単位下三角行列だから。)
計算量は LU 分解そのものが $O(N^3)$、対角成分の積取りが $O(N)$。 単純な余因子展開の $O(N!)$ から劇的に改善される。
符号 (sign) の追跡
部分ピボッティングでは各ステップで行交換が発生し得る。
各交換で $\det$ の符号が反転するため、sangi は交換回数を int カウンタで保持し、
最後に対角積へ $(-1)^s$ を掛ける。
オーバーフロー耐性
$\prod u_{ii}$ は $N$ が大きいと容易にオーバーフローまたはアンダーフローする。 対数領域で計算する手法は次の形になる:
$$\log |\det A| = \sum_i \log |u_{ii}|, \quad \mathrm{sign}(\det A) = (-1)^s \cdot \prod_i \mathrm{sign}(u_{ii})$$
統計学 (尤度関数)、機械学習 (Gaussian の対数尤度) では対数領域での累積が標準である。
専用関数 logdet の public API としての追加は今後の課題である。
関連記事: 行基本変形による行列式 / LU 分解 / LU 分解 (数値解析)
Bareiss 法 (整数領域)
要素型が Int (sangi 多倍長整数) や Rational の場合、
LU 分解の浮動小数点版は不適切である。
ピボット成分による除算が分数になり、要素サイズが指数的に膨らみ得る (中間係数膨張、intermediate expression swell)。
Bareiss 法 (1968) はこの問題を解決するアルゴリズムで、 $A^{(k)}$ の更新式に過去のピボットによる正確除算を組み込み、要素を常に整数に保つ:
$$a^{(k+1)}_{ij} = \frac{a^{(k)}_{kk} a^{(k)}_{ij} - a^{(k)}_{ik} a^{(k)}_{kj}}{a^{(k-1)}_{k-1,k-1}}$$
($a^{(-1)}_{-1,-1} := 1$ と約束)
分子は明らかに整数で、分母も「Bareiss の定理」によって割り切れることが保証される。 中間係数の絶対値は Hadamard の不等式により $\leq N^{N/2} \cdot M^N$ ($M = \max |a_{ij}|$) で抑えられる。
最終的に $\det A = a^{(N-1)}_{N-1,N-1}$ がそのまま得られる。 浮動小数点演算は一切使わず、丸め誤差ゼロ。
計算量は $O(N^3)$ の多倍長除算で、要素サイズを $L$ ビットとすると $O(N^3 L M(L) / L) = O(N^3 M(L))$ ($M(L)$ は $L$ ビット乗算のコスト) になる。 浮動小数点 LU より定数倍は大きいが、整数 / 有理数領域での正確性は計算機代数で必須である。
sangi の det は要素型を見て、浮動小数点系なら LU、整数 / 有理数系なら Bareiss を自動選択する。
閉形式の逆行列
$N = 2, 3$ では adjugate (余因子転置) 行列による閉形式が高速かつ数値的に十分。
2×2
$$A^{-1} = \frac{1}{\det A} \begin{pmatrix} a_{11} & -a_{01} \\ -a_{10} & a_{00} \end{pmatrix}$$
3×3
各成分は $2 \times 2$ 余因子で表される:
$$(A^{-1})_{ij} = \frac{(-1)^{i+j}}{\det A} \det M_{ji}$$
($M_{ji}$ は $j$ 行 $i$ 列を除いた小行列で、添字が転置していることに注意)
sangi は $N = 2, 3$ で完全アンロールされた閉形式を用いる。 $4 \times 4$ も閉形式は存在するが、項数が膨れ ($24$ 項以上)、LU 経路と同等以下の速度になるため、 $N \geq 4$ では LU 経路に切り替える設計である。
Gauss-Jordan による逆行列
Gauss-Jordan 消去は $[A \mid I]$ の拡大行列を行基本変形で $[I \mid A^{-1}]$ に持ち込む手順である:
- $A$ の左半分を行基本変形で行階段標準形に
- 同じ操作を右半分の $I$ に適用
- 左半分が $I$ になったとき、右半分が $A^{-1}$
計算量は $O(N^3)$。LU 分解より大まかには 2 倍遅い (前進消去と後退代入の両方を行うため)。
現行 sangi の inverse は次節の LU 経路をデフォルトとし、Gauss-Jordan 専用の inverse_gauss_jordan は未提供である (今後追加予定)。
関連記事: Gauss-Jordan 消去 / ピボット選択
LU 経由の逆行列
$A^{-1}$ は方程式 $AX = I$ の解である。 $X$ を列ごと $\mathbf{x}_j$ に分解すると、各列は $A \mathbf{x}_j = \mathbf{e}_j$ の解となる。
$PA = LU$ を先に求めておけば、$N$ 本の三角系を解くだけで全列が得られる:
$$L \mathbf{y}_j = P \mathbf{e}_j, \quad U \mathbf{x}_j = \mathbf{y}_j$$
計算量の内訳
- LU 分解: $O(N^3 / 3)$ 程度
- $N$ 本の三角系: 各 $O(N^2)$、合計 $O(N^3)$
- 合計: $O(N^3)$
右辺が単位行列という構造を利用すれば、前進代入で $L \mathbf{y}_j = P \mathbf{e}_j$ を解くとき
最初の数行をスキップできるため、定数倍がさらに改善される。
sangi の inverse はこの最適化を含めて実装している。
特異性と条件数
$A$ が特異 (非可逆) であれば $A^{-1}$ は存在しない。 LU 分解の途中でピボットがゼロになった瞬間、$\det A = 0$ かつ $A$ は特異と判定できる。
実数値の場合、ピボットが厳密にゼロになる確率はゼロだが、 非常に小さい値で「ほぼ特異」な行列が現れる。 この状況を数値的に評価するのが条件数である:
$$\kappa(A) = \|A\| \cdot \|A^{-1}\|$$
$\kappa(A) \gg 1$ なら $A^{-1}$ の計算結果はわずかな入力誤差で大きく変動し、信頼性が低い。 IEEE 754 倍精度 ($\epsilon \approx 10^{-16}$) で $\kappa(A) > 10^{14}$ ならば、 $A^{-1}$ の結果は数値的に意味のある精度を持たない。
条件数の正確計算には SVD が必要で $O(N^3)$ かかるが、
LU 経由の安価な推定 (LAPACK の dgecon 相当) もある。
sangi では条件数の詳細を別記事で扱う。
なぜ逆行列を作らず solve を使うか
線形方程式 $A \mathbf{x} = \mathbf{b}$ を解くために、 数学の教科書では「両辺に $A^{-1}$ を掛けて $\mathbf{x} = A^{-1} \mathbf{b}$」と書く。 しかし数値計算では逆行列を陽に作るのは ほぼ常に悪手である。 理由は以下のとおり。
コスト
- $A^{-1}$ を作って掛ける: $O(N^3)$ + $O(N^2) \approx O(N^3)$
- LU で直接 solve: $O(N^3 / 3) + O(N^2) \approx O(N^3 / 3)$
逆行列構築は 3 倍程度のコストを支払った上で、得られた $A^{-1}$ は捨てられることが多い。
数値精度
$A^{-1}$ の各要素は LU の三角系を解くより誤差が大きくなる。 $A^{-1} \mathbf{b}$ の積は浮動小数点演算の縮約で更に誤差が乗る。 LU 分解 + 三角系 solve なら、前進・後退代入の誤差解析が古典的に整理されており、 後退安定性 (backward stability) が保証される。
例外
同じ $A$ に対して数千本の右辺ベクトルを順に解く場合、
$A^{-1}$ を一度だけ作って各 $\mathbf{b}_i$ に作用させる方が安いケースはある。
ただしこの場合も LU.solve(b) を繰り返す方が精度・速度ともに優れ、
$A^{-1}$ を作る正当な理由は外部ライブラリインタフェースで明示的に要求された場合に限る。
sangi の inverse 関数は提供されているが、API 解説では
solve(A, b) または A.lu().solve(b) を優先するよう推奨している。
関連記事: ガウス消去法
設計判断
小サイズで閉形式、中以上で LU
$N \leq 3$ は閉形式 (Sarrus、adjugate) が高速で、$N = 4$ は閉形式と LU が拮抗するが 実装の単純さから sangi は $N \leq 4$ までを閉形式、$N \geq 5$ を LU としている。 SE(3) ($4 \times 4$) を完全に閉形式で扱える設計上のメリットは大きい。
整数領域では Bareiss
計算機代数の用途では、浮動小数点 LU は精度保証が崩れる。 Bareiss 法は中間係数膨張を制御する古典的解で、sangi は要素型を見て自動切替する。 Rational 要素の場合も Bareiss を適用し、分数簡約のコストを最小化できる。
solve を推奨する API 設計
inverse は教科書互換のために提供するが、ドキュメントでは solve を優先するよう繰り返し書く。
多くのユーザは数学の表記からそのままコードに翻訳するため、
「$A^{-1} b$ と書きたければ solve(A, b)」というメッセージを
ドキュメント・サンプルコード両方で徹底する必要がある。
参考文献
- Bareiss, E. H. (1968). "Sylvester's Identity and Multistep Integer-Preserving Gaussian Elimination". Mathematics of Computation, 22 (103), 565–578.
- Golub, G. H. and Van Loan, C. F. Matrix Computations. 4th ed. Johns Hopkins University Press, 2013.
- Higham, N. J. Accuracy and Stability of Numerical Algorithms. 2nd ed. SIAM, 2002. (逆行列 vs solve の数値安定性議論)
- Trefethen, L. N. and Bau, D. Numerical Linear Algebra. SIAM, 1997. (条件数と数値線形代数の基礎)