JP / EN

Gaussian Elimination & LU Decomposition

連立一次方程式の数値解法アルゴリズム

1. ガウスの消去法 (Gaussian Elimination)

連立一次方程式 $\boldsymbol{A}\boldsymbol{x} = \boldsymbol{b}$ を解くための最も基本的な直接法です。行列を上三角化する「前進消去」と、解を求める「後退代入」で構成されます。

① 前進消去 (Forward Elimination)

第 $k$ ステップにおいて、第 $k$ 行を用いてそれより下の行（$i > k$）の要素を消去します。

倍率 $m_{ik}$ の計算：

\[ m_{ik} = \frac{a_{ik}^{(k)}}{a_{kk}^{(k)}} \]

要素の更新：

\[ a_{ij}^{(k+1)} = a_{ij}^{(k)} - m_{ik} a_{kj}^{(k)} \quad (j = k+1, \dots, n) \] \[ b_{i}^{(k+1)} = b_{i}^{(k)} - m_{ik} b_{k}^{(k)} \]

② 後退代入 (Backward Substitution)

変形された上三角行列から、下の変数から順に値を決定します。

\[ x_n = \frac{b_n}{a_{nn}} \] \[ x_i = \frac{1}{a_{ii}} \left( b_i - \sum_{j=i+1}^{n} a_{ij} x_j \right) \quad (i = n-1, \dots, 1) \]

2. ピボット演算 (Pivoting)

数値計算の精度と安定性を保つための必須の処理です。

部分選択法 (Partial Pivoting):
前進消去の各ステップ $k$ で、第 $k$ 列の $|a_{kk}|$ から $|a_{nk}|$ の中で絶対値が最大の要素を持つ行を探し、$k$ 行目と入れ替えます。これにより、小さな値での除算による桁落ち（丸め誤差）を防ぎます。

3. LU分解 (LU Decomposition)

ガウスの消去法の計算過程を整理し、行列 $\boldsymbol{A}$ を下三角行列 $\boldsymbol{L}$ と上三角行列 $\boldsymbol{U}$ の積に分解する手法です（$\boldsymbol{A} = \boldsymbol{LU}$）。

Doolittle法による漸化式：

\[ u_{ij} = a_{ij} - \sum_{k=1}^{i-1} l_{ik} u_{kj} \quad (i \le j) \] \[ l_{ij} = \frac{1}{u_{jj}} \left( a_{ij} - \sum_{k=1}^{j-1} l_{ik} u_{kj} \right) \quad (i > j) \]

一度LU分解を行えば、右辺ベクトル $\boldsymbol{b}$ が変わっても代入操作のみで高速に解けるため、ニュートン法などの反復解法に最適です。

4. スパース行列処理の数学的背景

回路行列などの大規模系では、非ゼロ要素の配置（構造）を数学的に最適化することが計算効率の鍵となります。

① グラフ理論による構造解析

行列 $\boldsymbol{A}$ のスパース構造は、隣接行列を持つ有向グラフ $G(\boldsymbol{A}) = (V, E)$ として表現できます。消去過程におけるピボット選択は、グラフにおける「ノードの除去」と、それに伴う「エッジの追加（フィルイン）」に対応します。

定理：消去グラフ (Elimination Graph)

ノード $v_k$ を消去する際、そのノードに隣接していた全てのノード間にエッジが存在しない場合、それらを結ぶ新しいエッジ（フィルイン）が生成される。このグラフの完全グラフ化を避けることが、スパース性維持の数学的本質である。

② Markowitz Criterion の数理

第 $k$ ステップにおいて、要素 $a_{ij}$ をピボットに選んだ際のフィルイン発生数の上界は、以下の **Markowitz Product** によって与えられます。

\[ M_{ij} = (r_i^{(k)} - 1)(c_j^{(k)} - 1) \]

ここで、$r_i^{(k)}$ は第 $i$ 行の非ゼロ要素数、$c_j^{(k)}$ は第 $j$ 列の非ゼロ要素数です。消去ごとに $M_{ij}$ を最小化する選択（局所最適化）を行うことで、全体としての $O(N^3)$ の計算量を劇的に削減します。

③ シンボリック分解 (Symbolic Factorization)

数値計算（定数計算）を開始する前に、行列の構造（インデックス）のみを用いて、どの位置にフィルインが発生するかをあらかじめ決定します。これにより、以下のメリットが得られます：

  • 静的メモリ割り当て：
    実際の数値計算中に動的なメモリ確保を排除し、キャッシュ効率を最大化する。
  • 演算のスケジューリング：
    依存関係のない演算を特定し、並列計算への最適化を可能にする。
Numerical Stability vs Sparsity:
数学的には Markowitz Product が最小の要素を選ぶべきですが、数値的安定性（ピボットの絶対値が小さすぎないこと）も考慮する必要があります。実装上は、安定性を担保する閾値を設けつつ、その範囲内で最小の Markowitz 積を持つ要素を選択する「Threshold Pivoting」が一般的です。

5. 具体的な計算例 (2変数)

以下の系を、ピボット選択を含むガウスの消去法で解きます。

\[ \begin{cases} 1x_1 + 2x_2 = 5 \\ 3x_1 + 4x_2 = 11 \end{cases} \]

Step 1: ピボット選択
$|3| > |1|$ のため、1行目と2行目を入れ替えます。

\[ \begin{bmatrix} 3 & 4 & | & 11 \\ 1 & 2 & | & 5 \end{bmatrix} \]

Step 2: 前進消去
1行目を $1/3$ 倍して2行目から引きます。

\[ \begin{bmatrix} 3 & 4 & | & 11 \\ 0 & 0.667 & | & 1.333 \end{bmatrix} \]

Step 3: 後退代入
$0.667x_2 = 1.333 \Rightarrow x_2 = 2$
$3x_1 + 4(2) = 11 \Rightarrow 3x_1 = 3 \Rightarrow x_1 = 1$

結果： $x_1 = 1, x_2 = 2$