Gauss-Jordanの消去法

Gauss-Jordanの消去法#

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

$n$ 個の未知数 $x_{j} (j = 1, \dots, n)$ と $n$ 個の方程式からなる連立 1 次方程式

\sum_{j = 1}^{n} a_{i j} x_{j} = b_{i} (i = 1, \dots, n)

について考える。

連立1次方程式を

\begin{array}{r} \begin{aligned} a_{11}^{(1)} x_{1} + a_{12}^{(1)} x_{2} + \dots + a_{1 n}^{(1)} x_{n} = b_{1}^{(1)} \\ a_{21}^{(1)} x_{1} + a_{22}^{(1)} x_{2} + \dots + a_{2 n}^{(1)} x_{n} = b_{2}^{(1)} \\ ⋮ \\ a_{n 1}^{(1)} x_{1} + a_{n 2}^{(1)} x_{2} + \dots + a_{n n}^{(1)} x_{n} = b_{n}^{(1)} \end{aligned} \end{array}

の形に書く。

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

前進消去#

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

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

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

\begin{array}{r} \begin{aligned} a_{i j}^{(2)} = a_{i j}^{(1)} - m_{i 1} a_{1 j}^{(1)} & (i = 2, \dots, n; j = 2, \dots, n) \\ b_{i}^{(2)} = b_{i}^{(1)} - m_{i 1} b_{1}^{(1)} & (i = 2, \dots, n) \end{aligned} \end{array}

とおくと、

\begin{array}{r} \begin{aligned} a_{11}^{(1)} x_{1} + a_{12}^{(1)} x_{2} + \dots + a_{1 n}^{(1)} x_{n} = b_{1}^{(1)} \\ a_{22}^{(2)} x_{2} + \dots + a_{2 n}^{(2)} x_{n} = b_{2}^{(2)} \\ ⋮ \\ a_{n 2}^{(2)} x_{2} + \dots + a_{n n}^{(2)} x_{n} = b_{n}^{(2)} \end{aligned} \end{array}

となる。

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

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

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

\begin{array}{r} \begin{aligned} a_{11}^{(1)} x_{1} + a_{12}^{(1)} x_{2} + a_{13}^{(1)} x_{3} + \dots + a_{1 n}^{(1)} x_{n} = b_{1}^{(1)} \\ a_{22}^{(2)} x_{2} + a_{23}^{(2)} x_{3} + \dots + a_{2 n}^{(2)} x_{n} = b_{2}^{(2)} \\ a_{33}^{(3)} x_{3} + \dots + a_{3 n}^{(3)} x_{n} = b_{3}^{(3)} \\ ⋱ ⋮ ⋮ \\ a_{n n}^{(n)} x_{n} = b_{n}^{(n)} \end{aligned} \end{array}

ただし、 $a_{11}^{(1)} \neq 0, a_{22}^{(2)} \neq 0, \dots, 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)}} (b_{i}^{(i)} - \sum_{j = i + 1}^{n} a_{i j}^{(i)} x_{j}) (i = n, n - 1, \dots, 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, a_{22}^{(2)} \neq 0, \dots, a_{n n}^{(n)} \neq 0$ を満たし、消去法が適用可能

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

一般化狭義優対角行列

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

| a_{i i} | \geq \sum_{j \neq i} | a_{i j} |, (1 \leq i \leq n)

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

また、正の実数 $d_{1}, \dots, d_{n}$ に対して

| a_{i i} | d_{i} > \sum_{j \neq i} | a_{i j} | d_{j} (i = 1, \dots, n)

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

M行列

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

枢軸選択#

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

第 $k$ 段の消去の前に、 $k$ 番目以降の方程式を入れ替えて、枢軸要素の絶対値 $| a_{k k}^{(k)} |$ がそれより下の位置にある要素の絶対値 $| a_{i k}^{(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). 線形計算の数理. 岩波書店.