LU・Cholesky・LDL 分解とピボット戦略

概要

連立方程式 $A\mathbf{x} = \mathbf{b}$ を解く最も基本的な経路が、ガウス消去法とその因数分解版である LU 分解である。 対称行列に対してはコストを半減できる Cholesky / LDL 分解が存在し、 係数行列の構造（正定値・不定値・疎・整数）に応じて適切な分解を選ぶことで、 計算量・メモリ・数値安定性を同時に最適化できる。

本記事では LU・Cholesky・LDL の数学的構造、ピボット選択戦略の意味、 Bareiss 法による整数行列式、ブロック LU の BLAS-3 化、Bunch-Kaufman 対角ピボット、 不完全 Cholesky 前処理子を順に述べ、sangi がどの場面でどの実装にディスパッチするかを解説する。

API の使い方については LinAlg API — LU 分解、 Cholesky 分解、 LDL 分解、 ガウス消去法、 ピボット戦略 を参照。

ガウス消去法と LU 分解の関係

ガウス消去法は、$A\mathbf{x} = \mathbf{b}$ の拡大係数行列 $[A \mid \mathbf{b}]$ に対して 前進消去で上三角化し、後退代入で解を得る古典的手法である。 この前進消去過程を行列の積として記述したものが LU 分解にほかならない。

第 $k$ 段で列 $k$ の対角下成分を消去する変換は、Frobenius 変換行列

$$M_k = I - \mathbf{m}_k \mathbf{e}_k^\top, \quad m_{ik} = a_{ik}^{(k)} / a_{kk}^{(k)} \;\; (i > k)$$

と書ける。$n-1$ 段の消去を合成すると

$$M_{n-1} \cdots M_1 A = U \quad\Longleftrightarrow\quad A = (M_1^{-1} \cdots M_{n-1}^{-1}) U = LU$$

となり、各 $M_k^{-1} = I + \mathbf{m}_k \mathbf{e}_k^\top$ の積として下三角行列 $L$ が組み上がる。 ピボット選択を加えると、$PA = LU$ という形になり、$P$ は行交換の置換行列である。

ガウス消去法と LU 分解は同じ計算を別の言語で記述したものであり、 $O(n^3)$ のフロップ数は変わらない。 LU 分解の利点は、係数行列が同じで右辺が異なる複数の連立方程式 $A\mathbf{x}_1 = \mathbf{b}_1, A\mathbf{x}_2 = \mathbf{b}_2, \ldots$ を解くときに、 分解結果を再利用して $O(n^2)$ で各右辺を処理できる点にある。 sangi の lu_decomposition はこの形で結果を返し、 gaussian_elimination は内部で同じ消去過程を実行して直接解を返す。

ピボット戦略

消去ステップ $k$ で対角要素 $a_{kk}^{(k)}$ が零あるいは小さい場合、$m_{ik}$ が発散・桁落ちし数値解が壊れる。 この回避策がピボット選択であり、sangi は PivotStrategy 列挙型で 3 通り提供する。

None（ピボットなし）

列内交換を行わず対角要素をそのまま使う。 対角優位行列や、構造的に $a_{kk} \neq 0$ が保証される場合のみ安全。 コストは最小だが、悪条件行列では完全に崩壊する。

Partial（部分ピボット）

第 $k$ 段で列 $k$ の対角以下から絶対値最大の要素を選び、その行を $k$ 行と交換する:

$$p = \arg\max_{i \geq k} \, |a_{ik}^{(k)}|$$

各段のコストは $O(n)$ で、トータル $O(n^2)$ の追加コストにとどまる。 実用上ほとんどの行列で十分な安定性が得られるため、sangi の既定値となっている。

Full（完全ピボット）

残り部分行列 $A^{(k)}_{[k:n,k:n]}$ 全体から絶対値最大の要素を選び、 対応する行と列の双方を交換する:

$$(p, q) = \arg\max_{i,j \geq k} \, |a_{ij}^{(k)}|$$

列交換が入るため逆置換 $Q$ を別途記録し、解は $\mathbf{x} = Q\,(LU)^{-1} P \mathbf{b}$ で得る。 各段のコストは $O((n-k)^2)$ で、トータル $O(n^3)$ となり Partial の 2 倍程度。 悪条件行列やランク欠落に近い行列で部分ピボットが失敗するときの保険として用いる。

Rook ピボットと閾値ピボット

Rook ピボットは部分ピボットと完全ピボットの中間で、 「同じ行内で絶対値最大」かつ「同じ列内で絶対値最大」となる要素（チェスのルークの動きから命名）を選ぶ。 平均的に $O(n^2)$ コストで完全ピボットに近い安定性が得られる。 閾値ピボット (threshold pivoting) は、現在のピボット候補が 最大候補の $\tau$ 倍以上（$0 < \tau \le 1$）あれば交換せずに採用する戦略で、 疎行列分解で fill-in を抑える目的で使われる。 sangi は密 LU では Partial / Full、疎 LU では閾値ピボットを内部で用いる。

Rational 型の特別扱い

$T = $ Rational の場合、絶対値最大ではなく高さ最小のピボットを選ぶ。 高さとは $h(p/q) = \max(|p|, |q|)$ であり、これが小さいピボットを選ぶことで 消去後の有理数演算で分子・分母が指数的に膨張するのを抑える。 浮動小数点とは異なり Rational は丸め誤差がないため、安定性ではなく演算コストの最適化が目的である。

Wilkinson の成長因子と数値安定性

ガウス消去で計算過程に現れる行列要素の最大値 $\|A^{(k)}\|_\infty$ と、 入力の最大値 $\|A\|_\infty$ の比を成長因子 (growth factor) と呼ぶ:

$$\rho_n(A) = \frac{\max_{i,j,k} \, |a_{ij}^{(k)}|}{\max_{i,j} \, |a_{ij}|}$$

後退誤差解析の枠組み (Wilkinson 1961) では、計算解 $\hat{\mathbf{x}}$ は摂動系 $(A + \Delta A) \hat{\mathbf{x}} = \mathbf{b}$ の厳密解であり、 $\|\Delta A\|_\infty$ は $\rho_n(A) \cdot u \cdot \|A\|_\infty$ のオーダーで抑えられる（$u$ はユニット丸め）。 従って成長因子が小さいほど後退誤差が小さい。

• 部分ピボット: 最悪 $\rho_n \le 2^{n-1}$ だが、実用上は $\rho_n \le 10$ 程度に収まる。Wilkinson 行列という人為的構成のみが指数増加を示す。
• 完全ピボット: $\rho_n \le n^{1/2} (2 \cdot 3^{1/2} \cdots n^{1/(n-1)})^{1/2}$ で、未だ未解決ながら多項式オーダーと予想されている。
• 対称正定値 (Cholesky): $\rho_n = 1$。成長は起こらず、ピボット選択そのものが不要。

悪条件行列の判定は条件数 $\kappa(A) = \|A\| \cdot \|A^{-1}\|$ を用いる。 $\kappa(A)$ が $1/u$ 程度に達するとピボット選択をしても解の有効桁が失われ、 SVD ベースの正則化または反復精化が必要になる。 sangi の svd_solve はこの場面で rcond によるダンピングを提供する。

Bareiss 法（整数係数の行列式）

整数係数の行列式を普通の浮動小数点 LU で計算すると、 除算で丸め誤差が混入する。 Rational で計算すれば厳密だが、各除算で分母子が膨張し $O(2^n)$ 規模の有理数が出現する。 Bareiss 法は除算なしに整数演算だけで行列式を得る古典的手法で、 Sylvester の恒等式に基づく。

第 $k$ 段の消去で更新後の要素を次のように定義する:

$$M_{ij}^{(k)} = \frac{M_{ij}^{(k-1)} M_{kk}^{(k-1)} - M_{ik}^{(k-1)} M_{kj}^{(k-1)}}{M_{k-1,k-1}^{(k-2)}}$$

右辺の分子は行列のミナーであり、Sylvester の決定要因恒等式によって $M_{k-1,k-1}^{(k-2)}$ で常に厳密に割り切れる。 従ってこの割り算は整数除算であり、Float も Rational も介在しない。 最終的に $M_{n-1,n-1}^{(n-2)}$ が $\det(A)$ そのものになる。

計算量は $O(n^3)$、ただし中間整数が指数的に大きくなる可能性があるため 多倍長整数 (sangi の Int) で扱う必要がある。 Hadamard の不等式により $|\det(A)| \le \prod_i \|\mathbf{a}_i\|_2$ から 中間結果のビット幅は $O(n \log(n M))$（$M$ は入力の最大絶対値）で抑えられる。

sangi の Matrix::determinant() メンバ関数は、 numeric_traits<T>::is_integer が真の型に対しては自動的に bareiss_determinant を呼び出す。 浮動小数点型では LU 経路（lu_determinant、対角積×ピボット符号）が選ばれる。

ブロック LU と BLAS-3 最適化

素朴な LU 実装はループ内で axpy 相当の Level-2 BLAS 演算を多用するため、 メモリ帯域がボトルネックになる。 ブロック LU は行列を $b \times b$ のブロックに分割し、 各段の更新を行列乗算 (Level-3 BLAS, gemm) に帰着させてキャッシュ再利用を最大化する。

$A$ を次のように分割する:

$$A = \begin{pmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{pmatrix}$$

$A_{11}$ ($b \times b$) を素朴な LU で分解して $L_{11} U_{11}$ を得たあと:

$$L_{21} = A_{21} U_{11}^{-1}, \quad U_{12} = L_{11}^{-1} A_{12}, \quad A_{22}^{(1)} = A_{22} - L_{21} U_{12}$$

を計算する。$L_{21}, U_{12}$ は三角ソルバ (trsm, Level-3)、 $A_{22}^{(1)}$ の更新は $\mathrm{gemm}$ で、それぞれフロップに対するメモリアクセス比が $O(1/b)$ に下がる。残り $(n-b) \times (n-b)$ のシュア補元 $A_{22}^{(1)}$ に再帰的に同じ手続きを適用する。

ブロックサイズ $b$ は L1/L2 キャッシュサイズと演算強度から決まり、 典型的には $b = 32 \sim 128$。 sangi は MKL/OpenBLAS が利用可能なときは dgetrf 経由でブロック LU を呼び、 フォールバック時は内部の Strassen 風 panel-update 実装を用いる。

Cholesky 分解

$A$ が対称かつ正定値（任意の $\mathbf{x} \neq 0$ で $\mathbf{x}^\top A \mathbf{x} > 0$）のとき、 $A = LL^\top$ なる下三角行列 $L$ が一意に存在する。これが Cholesky 分解である。

対角要素と非対角要素の漸化式は次のとおり:

$$L_{kk} = \sqrt{A_{kk} - \sum_{j=1}^{k-1} L_{kj}^2}, \quad L_{ik} = \frac{1}{L_{kk}} \left(A_{ik} - \sum_{j=1}^{k-1} L_{ij} L_{kj}\right) \quad (i > k)$$

各段で対称性により上半分を計算せずに済み、フロップ数は LU の半分、すなわち $n^3/3$ となる。 正定値性により $L_{kk}^2 = A_{kk} - \sum L_{kj}^2 > 0$ が保証されピボット選択は不要、 さらに成長因子 $\rho_n = 1$ で数値的にも極めて安定である。

用途は最小二乗法の正規方程式 $A^\top A \mathbf{x} = A^\top \mathbf{b}$、 カルマンフィルタの共分散行列、ガウス過程の対数尤度、 二次計画法の KKT 系の対角ブロックなど、正定値性が構造的に保証される問題に広い。

ブロック化版（ブロック Cholesky）も同様に存在し、対角ブロックは小さな Cholesky、 オフ対角は trsm、シュア補元更新は syrk（対称ランク-$k$ 更新）で実装する。

LDL 分解と Bunch-Kaufman ピボット

対称行列が正定値とは限らない場合、Cholesky の平方根は破綻する。 この一般化が LDL 分解 $A = L D L^\top$ で、$L$ は単位下三角、$D$ は対角行列である:

$$D_{kk} = A_{kk} - \sum_{j=1}^{k-1} D_{jj} L_{kj}^2, \quad L_{ik} = \frac{1}{D_{kk}} \left(A_{ik} - \sum_{j=1}^{k-1} D_{jj} L_{ij} L_{kj}\right)$$

平方根が消えるため $\sqrt{\,}$ の数値コストと、負の引数による失敗の双方を回避できる。 計算量は Cholesky と同じ $n^3/3$。

Bunch-Kaufman 対角ピボット

純粋な LDL は対角要素が小さいまたは零の場合に破綻する。 Bunch-Kaufman ピボットは、対角に $1\times1$ か $2\times2$ のブロックを許す $A = L D L^\top$ 分解で、$D$ は対角ブロック行列となる:

$$D = \mathrm{blockdiag}(D_1, D_2, \ldots), \quad D_k \in \mathbb{R}^{1\times 1} \cup \mathbb{R}^{2\times 2}$$

各段で「対称交換のみで対角の $1\times1$ ピボットが安全か」を判定し、 危険な場合は隣接する 2 つの対角要素を $2\times2$ ブロックとしてまとめる。 この戦略により対称構造を保ったまま不定値・特異に近い対称行列を扱える。

用途は最適化問題の鞍点系（KKT 行列）、不定値ヘッシアン、 対称な内部点法の反復で生じるブロック対角サドル系などである。

不完全 Cholesky 前処理子

大規模疎行列に対する反復ソルバー（共役勾配 CG、MINRES 等）の収束を加速するために、 近似的な分解を前処理子 $M \approx A$ として用いる。 Cholesky の対称構造を保ちつつ fill-in を制限したものが不完全 Cholesky 分解 (IC) である。

レベル $0$ の不完全 Cholesky $\mathrm{IC}(0)$ は、$A$ の非零パターン $\mathcal{S} = \{(i,j) : A_{ij} \neq 0\}$ を保ち、 消去で発生する fill-in を捨てる:

$$\tilde{L}_{ik} = \begin{cases} (A_{ik} - \sum_{j<k} \tilde{L}_{ij} \tilde{L}_{kj}) / \tilde{L}_{kk} & (i,k) \in \mathcal{S} \\ 0 & \text{otherwise} \end{cases}$$

$\tilde{M} = \tilde{L} \tilde{L}^\top$ は厳密な $A$ ではないが、 前処理付き CG (PCG) の反復で $\tilde{M}^{-1} \mathbf{r}$ を計算すれば 条件数 $\kappa(\tilde{M}^{-1} A)$ が大幅に改善し収束が速くなる。 $\mathrm{IC}(p)$ は $p$ 段までの fill-in を許す変種、 modified IC は対角を行和で補正する変種で、それぞれ問題に応じて選ぶ。

不完全 Cholesky は対称正定値の場合のみ「分解失敗 (negative pivot) なし」が保証されるわけではなく、 対角に微小なシフト $\alpha I$ を足して破綻を回避する Manteuffel シフトが実用的に用いられる。

不完全 LU (ILU) は対称でない場合の対応物で、同じ概念を $L, U$ ペアに適用する。 sangi の疎行列反復ソルバーは 反復法と疎行列のページで扱う。

sangi の自動 dispatch

sangi LinAlg の分解関数は、入力行列の型と構造に応じて内部実装を選択する。 利用者は通常 lu_solve、cholesky_solve、ldl_solve 等の高水準 API のみを呼べばよい。

サイズに応じて、sangi は内部で BLAS-3 ブロック更新 (MKL/OpenBLAS フォールバック) に切り替える。 $n < 64$ では素朴 LU、$n \ge 64$ ではブロック LU、 $n \ge 1024$ では並列パネル分解を用いる。 この閾値はベンチマークに基づいており、ユーザーが意識する必要はない。

参考文献

• Golub, G. H. and Van Loan, C. F. (2013). Matrix Computations, 4th ed., Johns Hopkins University Press. Chapters 3–4.
• Wilkinson, J. H. (1961). "Error analysis of direct methods of matrix inversion". J. ACM, 8(3), 281–330.
• Bareiss, E. H. (1968). "Sylvester's identity and multistep integer-preserving Gaussian elimination". Math. Comp., 22(103), 565–578.
• Bunch, J. R. and Kaufman, L. (1977). "Some stable methods for calculating inertia and solving symmetric linear systems". Math. Comp., 31(137), 163–179.
• Higham, N. J. (2002). Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM.