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)

について考える。

連立1次方程式を

Undefined control sequence: \DeclareMathOperator at position 57: …{#1}}
% 演算子の定義
\̲D̲e̲c̲l̲a̲r̲e̲M̲a̲t̲h̲O̲p̲e̲r̲a̲t̲o̲r̲{\im}{ \text{Im…

% 太字の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}

の形に書く。

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{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}

(3)

とおくと、

\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}

(4)

となる。

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

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

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

\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}

(5)

ただし、 $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)}

(6)

なので、

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

(7)

として解が求まる。

$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)}

(8)

となっているので

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

(9)

として解が求まる。

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

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)

(10)

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

# データの作成
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.0000e-03,  6.0000e+00],
       [ 0.0000e+00,  1.8005e+04]])
b=array([6.0010e+00, 1.8005e+04])

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})$ .

枢軸選択¶

枢軸要素 $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). 線形計算の数理. 岩波書店.