線形連立方程式を解く

$A$ は，逆行列を持つ $n \times n$ の行列， $b$ は $n$ 次のベクトルとして，

Ax = b

を満たすベクトル $x \in \mathbb{R}^{n}$ を求める．

ガウスの消去法のアルゴリズム

ガウスの消去法は，連立方程式 $Ax=b$ を上三角行列 $U$ に関する連立方程式 $Ux = c$ に変換して解く方法である．

$A = (a_{ij})_{ij})$ ， $b = (b_1, \ldots, b_n)$ ， $x = (x_1, \ldots, x_n)$ として，連立方程式 $Ax = b$ を成分ごとに書くと，

\begin{align*} a_{11} x_1 + a_{12} x_2 + \cdots + a_{1n} x_n &= b_1 \\ a_{21} x_1 + a_{22} x_2 + \cdots + a_{2n} x_n &= b_2 \\ \hspace{1.5em} \vdots & \hspace{1.5em} \vdots\\ a_{n1} x_1 + a_{n2} x_2 + \cdots + a_{nn} x_n &= b_n \end{align*}

第1式を使って，第2式以降の $x_1$ を消去すると，

\begin{align*} a_{11} x_1 + a_{12} x_2 + \cdots + a_{1n} x_n &= b_1 \\ a_{22}^{(1)} x_2 + \cdots + a_{2n}^{(1)} x_n &= b_2^{(1)} \\ \hspace{1.5em} \vdots & \hspace{1.5em} \vdots\\ a_{n2}^{(1)} x_2 + \cdots + a_{nn}^{(1)} x_n &= b_n^{(1)} \end{align*}

となる．ただし， $a_{ij}^{(1)} = a_{ij} - (a_{i1}/a_{11})a_{1j}$ ， $b_i^{(1)} = b_i - (a_{i1}/a_{11})b_1$

あとは，これ↑を繰り返すだけである．

ガウスの消去法のアルゴリズムは以下のようになる． $A \in \mathbb{R}^{n\times n}$ に対して， $\mathcal{O}(n^3)$ のアルゴリズムになる．

Input: $A \in \mathbb{R}^{n\times n}$ , $b \in \mathbb{R}^n$

for $k = 1$ to $n-1$ do:

$d := 1/a_{kk}$
for $i = k + 1$ to $n$ do:
- $a_{ik} := a_{ik} \cdot d$
- $b_i := b_i - b_k \cdot a_{ik}$
- for $j=k+1$ to $n$ do:
  - $a_{ij} = a_{ij} - a_{ik} \cdot a_{kj}$

ガウスの消去法では，対角成分 $a_{ii}$ が $0$ だと計算が破綻する．
また， $0$ ではなくても，とても小さい数だと不可逆な(逆を辿っても $A$ に戻らない)変換になってしまう．

ガウスの消去法の反復の $k$ 回目では，行列 $A$ は，

A^{(k)} = \begin{pmatrix} a_{11} & \ast & \ast & & \cdots & \ast\\ 0 & a_{22}^{(2)} & \ast & & \cdots & \ast \\ \\ & & \ddots & & \\ \\ 0 & 0 & & a_{kk}^{(k)} & \cdots & a_{kn}^{(k)} \\ \\ \vdots & \vdots & & \vdots & & \vdots \\ \\ 0 & 0 & & a_{nk}^{(k)} & \cdots & a_{nn}^{(k)} \end{pmatrix}

のようになっている．ここで， $a_{kk}^{(k)}$ が $0$ や $0$ に近い値の時，ガウスの消去法は破綻する．
そこで， $k$ 行を他の行 $l \in \{ k, k+1, \ldots, n \}$ と入れ替える．

入れ替える行 $l$ は，

l = \argmax_{k \leq i \leq n} |a_{ik}^{(k)}|

と選ぶ．すなわち， $k$ 列め一番大きいものを探す．

$A^{(k)}$ の行の入れ替えを行ったら， $b^{(k)}$ の行の入れ替えも忘れずに行う．

以上の手続きは，行列 $A$ を下三角行列 $L$ と上三角行列 $U$ に分解することに対応する．

PA = LU

$P$ はピボットに対応する置換行列である．

$A$ が対称行列の時は，

A = LDL^\mathrm{T}

と下三角行列 $L$ と対角行列 $D$ でLU分解できる．
これをコレスキー分解という．

$A$ が正定値 ( $x^\mathrm{T} A x > 0$ ) だとコレスキー分解は必ず存在する．

ここで， $D^{1/2} = \mathrm{diag} (\sqrt{d_{11}}, \dots, \sqrt{d_{nn}})$ とすると，

A = LD^{1/2}(LD^{1/2})^\mathrm{T} \equiv \tilde{L} \tilde{L}^\mathrm{T}

と表せる．これを成分で表すと，

a_{ij} = \sum_{k=1}^{j} l_{ik}l_{jk} = \sum_{k=1}^{j-1} l_{ik}l_{jk} + l_{ij} l_{jj}

$l_{ij}$ について解けば，

\begin{align*} l_{11} &= \sqrt{a_{11}} \\ l_{jj} &= \left(a_{jj} - \sum_{k=1}^{j-1} l_{jk}^2 \right)^{1/2} \\ l_{ij} &= \frac{1}{l_{jj}} \left(a_{ij} - \sum_{k=1}^{j-1} l_{ik}l_{jk} \right) \end{align*}

完全には解けていないが， $j=1,2,\ldots, n$ と反復的に求めることができる．