Gauss-Jordanの消去法

Gauss-Jordanの消去法#

ガウスの消去法 （Gauss elimination）、 ガウス＝ジョルダンの消去法 （Gauss-Jordan elimination）、あるいは 掃き出し法 と呼ばれる方法。

\(n\) 個の未知数 \(x_j(j=1, \cdots, n)\) と \(n\) 個の方程式からなる連立 1 次方程式

\[ \sum_{j=1}^n a_{i j} x_j=b_i \quad(i=1, \cdots, n) \]

について考える。

連立1次方程式を

\[\begin{split} % 太字のalias \newcommand{\b}[1]{\boldsymbol{#1}} % 演算子の定義 \DeclareMathOperator{\im}{ \text{Im} } \DeclareMathOperator{\rank}{ \text{rank} } \DeclareMathOperator{\span}{ \text{Span} } \DeclareMathOperator{\Ker}{ \text{Ker} } % \begin{align} & a_{11}^{(1)} x_1+a_{12}^{(1)} x_2+\cdots+a_{1 n}^{(1)} x_n=b_1^{(1)}\\ & a_{21}^{(1)} x_1+a_{22}^{(1)} x_2+\cdots+a_{2 n}^{(1)} x_n=b_2^{(1)}\\ &\quad \vdots\\ & a_{n 1}^{(1)} x_1+a_{n 2}^{(1)} x_2+\cdots+a_{n n}^{(1)} x_n=b_n^{(1)} \\ \end{align} \end{split}\]

の形に書く。

Gaussの消去法では、前進消去と後退代入という操作によって方程式を解く。

前進消去#

(1) 1番目の未知数\(x_1\)に着目し、2~\(n\)番目の方程式から消去する。

\(i=2,\cdots,n\)について、\(i\)番目の式から\(1\)番目の式の\(m_{i1} = a_{i1}^{(1)} / a_{11}^{(1)}\)倍を引けばよい（\(a_{11}^{(1)}\neq 0\)と仮定しておく）

この操作のあとの連立1次方程式は、

\[\begin{split} \begin{align} a_{i j}^{(2)}=a_{i j}^{(1)}-m_{i 1} a_{1 j}^{(1)} & \quad(i=2, \cdots, n ; j=2, \cdots, n) \\ b_i^{(2)}=b_i^{(1)}-m_{i 1} b_1^{(1)} & \quad (i=2, \cdots, n) \end{align} \end{split}\]

とおくと、

\[\begin{split} \begin{aligned} a_{11}^{(1)} x_1+a_{12}^{(1)} x_2+\cdots+a_{1 n}^{(1)} x_n=b_1^{(1)}&\\ a_{22}^{(2)} x_2+\cdots+a_{2 n}^{(2)} x_n=b_2^{(2)}&\\ \vdots \quad &\\ a_{n 2}^{(2)} x_2+\cdots+a_{n n}^{(2)} x_n=b_n^{(2)}&\\ \end{aligned} \end{split}\]

となる。

この消去操作において着目した行列要素の位置 \((1,1)\)を枢軸（pivot）、その要素\(a_{11}^{(1)}\)を 枢軸要素 、 \(m_{i1}\)を乗数と呼ぶ。

(2) 2番目の未知数\(x_2\)に着目し、3~\(n\)番目の方程式から消去する。

以下、このような操作を\(n-1\)回まで繰り返す。そして以下を得る

\[\begin{split} \begin{aligned} a_{11}^{(1)} x_1+a_{12}^{(1)} x_2+a_{13}^{(1)} x_3+\cdots+a_{1 n}^{(1)} x_n =b_1^{(1)}&\\ a_{22}^{(2)} x_2+a_{23}^{(2)} x_3+\cdots+a_{2 n}^{(2)} x_n=b_2^{(2)}&\\ a_{33}^{(3)} x_3+\cdots+a_{3 n}^{(3)} x_n=b_3^{(3)}&\\ \ddots \quad \quad \vdots \quad \quad \vdots \quad&\\ a_{n n}^{(n)} x_n=b_n^{(n)} &\\ \end{aligned} \end{split}\]

ただし、\(a_{11}^{(1)} \neq 0, \quad a_{22}^{(2)} \neq 0, \quad \cdots, \quad a_{n n}^{(n)} \neq 0\)と仮定している。

後退代入#

\(x_n\)については、

\[ a_{n n}^{(n)} x_n=b_n^{(n)} \]

なので、

\[ x_n = \frac{ b_n^{(n)} }{ a_{n n}^{(n)} } \]

として解が求まる。

\(x_{n-1}\)については

\[ a_{n-1, n-1}^{(n-1)} x_{n-1} + a_{n n}^{(n-1)} x_n = b_{n-1}^{(n-1)} \]

となっているので

\[ x_{n-1} = \frac{1}{a_{n-1, n-1}^{(n-1)}} (b_{n-1}^{(n-1)} - a_{n n}^{(n-1)} x_n ) \]

として解が求まる。

これを一般化すると、解は

\[ x_i=\frac{1}{a_{i i}^{(i)}}\left(b_i^{(i)}-\sum_{j=i+1}^n a_{i j}^{(i)} x_j\right) \quad(i=n, n-1, \cdots, 1) \]

によって求まる。これを 後退代入 と呼ぶ。

# データの作成
import numpy as np

A = np.array([
    [1,2,3],
    [4,5,6],
    [7,8,0],
])
x_true = np.array([1, 2, 3])
b = A @ x_true

def forward_elimination(A, b):
    # 前進消去
    n = len(b)
    for k in range(n):
        for i in range(k + 1, n):
            m = A[i,k] / A[k,k]
            # numpyのベクトル演算に頼っているが、本来はここもA[i,]の各j要素にわたってfor j in range(n)が必要 → 計算量はO(n^3)
            A[i,k:] = A[i,k:] - m * A[k,k:]
            b[i] = b[i] - m * b[k]
            # print(f"{i=}, {j=}, {m=}")
            # print(A)
    return A, b

A, b = forward_elimination(A, b)
print(f"{A=}")
print(f"{b=}")

A=array([[ 1,  2,  3],
       [ 0, -3, -6],
       [ 0,  0, -9]])
b=array([ 14, -24, -27])

def backward_substitution(A, b):
    # 後退代入
    n = len(b)
    x = np.zeros(shape=(n,))
    for i in range(n-1, -1, -1):
        x[i] = (1 / A[i,i]) * (b[i] - A[i,i+1:] @ x[i+1:])
    return x

x = backward_substitution(A, b)
print(f"{x=}")

x=array([1., 2., 3.])

消去法が使える行列#

以下の条件のいずれかを満たす行列\(A\)は、\(a_{11}^{(1)} \neq 0, \quad a_{22}^{(2)} \neq 0, \quad \cdots, \quad a_{n n}^{(n)} \neq 0\)を満たし、消去法が適用可能

行方向に一般化狭義優対角.
列方向に一般化狭義優対角.
(正則な) \(\mathrm{M}\) 行列.
対称部分が正定値: \(\boldsymbol{x}^{\top} A \boldsymbol{x}=\frac{1}{2} \boldsymbol{x}^{\top}\left(A+A^{\top}\right) \boldsymbol{x}>0 \quad(\forall \boldsymbol{x} \neq \mathbf{0})\).

一般化狭義優対角行列

各行において, 対角要素の絶対値が非対角要素の絶対値の和に比べて大きい行列, すなわち,

\[ \left|a_{i i}\right| \geq \sum_{j \neq i}\left|a_{i j}\right| , \quad (1 \leq i \leq n) \]

を満たす行列 \(A\) を(行方向の) 狭義優対角行列 と呼ぶ.

また、正の実数\(d_1, \cdots, d_n\) に対して

\[ \left|a_{i i}\right| d_i>\sum_{j \neq i}\left|a_{i j}\right| d_j \quad(i=1, \cdots, n) \]

が成り立つとき， \(A\) を(行方向の) 一般化狭義優対角行列 と呼ぶ

M行列

正則な実行列\(A=(a_{ij})\)について、非対角要素がすべて非正\(a_{ij} \leq 0 (i \neq j)\)であり、逆行列\(A^{-1}\)が非負行列（各要素が非負の実数である行列）であるとき、\(A\)を M行列 と呼ぶ。

枢軸選択#

枢軸要素\(a_{kk}^{(k)}\)が0になると次の段に進めなくなるし、0に近い値だと誤差の拡大が起きて精度の良い解が得られない。

第\(k\)段の消去の前に、\(k\)番目以降の方程式を入れ替えて、枢軸要素の絶対値\(|a_{kk}^{(k)}|\)がそれより下の位置にある要素の絶対値\(|a_{ik}^{(k)}| (k\leq i \leq n)\)の中で最大になるようにすることができる。この操作を 枢軸選択 （pivoting）あるいは 部分枢軸選択 （partial pivoting）と呼ぶ。

import numpy as np

A = np.array([
    [-0.001, 6],
    [3, 5],
])
x_true = np.array([-1, 1])
b = A @ x_true
print(f"{b=}")

b=array([6.001, 2.   ])

A, b = forward_elimination(A, b)
print(f"{A=}")
print(f"{b=}")

x = backward_substitution(A, b)
print(f"{x=}")

A=array([[-1.0000e-03,  6.0000e+00],
       [ 0.0000e+00,  1.8005e+04]])
b=array([6.0010e+00, 1.8005e+04])
x=array([-1.,  1.])

スケーリング#

各方程式にゼロでない定数を掛けることをスケーリングと呼ぶ。

数学的には解は不変のはずだが、枢軸選択の結果はスケーリングに依存してしまう問題がある。

参考文献#

杉原正顯, & 室田一雄. (2009). 線形計算の数理. 岩波書店.