📂行列代数

ガウス消去法を使った逆行列の求め方アルゴリズム

アルゴリズム

Input

可逆行列 $M \in \mathbb{R}^{n \times n}$ が与えられているとする。

---

Step 1. 初期化

$M$ と同じサイズの単位行列 $X$ を作る。

---

Step 2. エシュロン形

ガウスの消去法を通じて $M$ をエシュロン形にする。本来なら上三角行列の形で十分だが、逆行列計算アルゴリズムではただ対角行列まで計算する。ただし、この時に $M$ に施された行操作を $X$ にも全く同様に施す。

---

Output

$MX = XM = E_{n}$ を満たす、つまり $M$ の逆行列である $X$ を得る。

例

詳細な説明を聞くよりも、一つの例を見る方が良い。

Input

行列 $M = \begin{bmatrix} 1 & 0 & 1 & 1\\ 2 & 0 & 1 & 0\\ -2& 3 & 4 & 0\\ -5& 5 & 6 & 0\end{bmatrix}$ が与えられているとする。

---

Step 1. 初期化

$$\begin{bmatrix} 1 & 0 & 1 & 1 \\ 2 & 0 & 1 & 0 \\ -2& 3 & 4 & 0 \\ -5& 5 & 6 & 0 \end{bmatrix} \cdots \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}$$

与えられた行列と同じサイズの単位行列を作る。上では左が $M$、右が後に逆行列となる $X$ だ。

---

Step 2. エシュロン形

厳密に言えば、このポストで紹介する方法でエシュロン形を作るわけではない。計算の便宜のために先に対角行列を作り、各行をピボットで割ることにする。

$$\begin{align*} & \begin{bmatrix} 1 & 0 & 1 & 1 \\ 2 & 0 & 1 & 0 \\ -2& 3 & 4 & 0 \\ -5& 5 & 6 & 0 \end{bmatrix} & \cdots & \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \\ \to & \begin{bmatrix} 1 & 0 & 1 & 1 \\ 0 & 0 & -1 & -2 \\ 0 & 3 & 6 & 2 \\ 0 & 5 & 11 & 5 \end{bmatrix} & \cdots & \begin{bmatrix} 1 & 0 & 0 & 0 \\ -2 & 1 & 0 & 0 \\ 2 & 0 & 1 & 0 \\ 5 & 0 & 0 & 1 \end{bmatrix} \end{align*}$$

$M$ に行操作を施すとき、$X$ にも全く同じように施すというのは、$M$ の一行目に $-2, 2, 5$ を掛けて二行目、三行目、四行目に加える操作を $X$ にも同様に行うということだ。

$$\begin{align*} & \begin{bmatrix} 1 & 0 & 1 & 1 \\ 0 & 0 & -1 & -2 \\ 0 & 3 & 6 & 2 \\ 0 & 5 & 11 & 5 \end{bmatrix} & \cdots & \begin{bmatrix} 1 & 0 & 0 & 0 \\ -2 & 1 & 0 & 0 \\ 2 & 0 & 1 & 0 \\ 5 & 0 & 0 & 1 \end{bmatrix} \\ \to & \begin{bmatrix} 1 & 0 & 1 & 1 \\ 0 & 3 & 6 & 2 \\ 0 & 0 & -1 & -2 \\ 0 & 5 & 11 & 5 \end{bmatrix} & \cdots & \begin{bmatrix} 1 & 0 & 0 & 0 \\ 2 & 0 & 1 & 0 \\ -2 & 1 & 0 & 0 \\ 5 & 0 & 0 & 1 \end{bmatrix} \end{align*}$$

ピボット^pivotが $0$ の場合、つまり $M_{ii} = 0$ であれば、 $j$ 番目の行と $i$ 番目の行を交換する。$M_{ji} \ne 0$ を満たす $j$ 行の存在は $M$ が可逆行列であるという仮定によって保証される。当然、この行交換操作も $X$ にまったく同様に施す必要がある。

$$\begin{align*} & \begin{bmatrix} 1 & 0 & 1 & 1 \\ 0 & 3 & 6 & 2 \\ 0 & 0 & -1 & -2 \\ 0 & 5 & 11 & 5 \end{bmatrix} & \cdots & \begin{bmatrix} 1 & 0 & 0 & 0 \\ 2 & 0 & 1 & 0 \\ -2 & 1 & 0 & 0 \\ 5 & 0 & 0 & 1 \end{bmatrix} \\ \to & \begin{bmatrix} 1 & 0 & 1 & 1 \\ 0 & 3 & 6 & 2 \\ 0 & 0 & -1 & -2 \\ 0 & 0 & 1 & 5/3 \end{bmatrix} & \cdots & \begin{bmatrix} 1 & 0 & 0 & 0 \\ 2 & 0 & 1 & 0 \\ -2 & 1 & 0 & 0 \\ 5/3 & 0 & -5/3 & 1 \end{bmatrix} \\ \to & \begin{bmatrix} 1 & 0 & 0 & 1 \\ 0 & 3 & 0 & -10 \\ 0 & 0 & -1 & -2 \\ 0 & 0 & 0 & -1/3 \end{bmatrix} & \cdots & \begin{bmatrix} -1 & 1 & 0 & 0 \\ -10 & 6 & 1 & 0 \\ -2 & 1 & 0 & 0 \\ -1/3 & 1 & -5/3 & 1 \end{bmatrix} \\ \to & \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 3 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & -1/3 \end{bmatrix} & \cdots & \begin{bmatrix} 0 & -2 & 5 & -3 \\ 0 & -24 & 51 & -30 \\ 0 & -5 & 10 & -6 \\ -1/3 & 1 & -5/3 & 1 \end{bmatrix} \end{align*}$$

$M$ が対角行列になるまで繰り返す。元々は上三角行列の形にしてから行のピボットで割ってエシュロン形にする。その後、他の行のピボット以降の成分を消去するが、コード作成が面倒なため順序を変更している。この順序変更はガウスの消去法の理論的背景から問題ない。

$$\begin{align*} & \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 3 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & -1/3 \end{bmatrix} & \cdots & \begin{bmatrix} 0 & -2 & 5 & -3 \\ 0 & -24 & 51 & -30 \\ 0 & -5 & 10 & -6 \\ -1/3 & 1 & -5/3 & 1 \end{bmatrix} \\ \to & \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} & \cdots & \begin{bmatrix} 0 & -2 & 5 & -3 \\ 0 & -8 & 17 & -10 \\ 0 & 5 & -10 & 6 \\ 1 & -3 & 5 & -3 \end{bmatrix} \end{align*}$$

対角行列を完成させたら、各行をそれぞれのピボットで割る。これで、$M$ は単位行列の形になるべきだ。

---

Output

$$X = \begin{bmatrix} 0 & -2 & 5 & -3 \\ 0 & -8 & 17 & -10 \\ 0 & 5 & -10 & 6 \\ 1 & -3 & 5 & -3 \end{bmatrix}$$

$X$ は $M$ の逆行列だ。今までの過程をまとめると、次のようになる。

$$\begin{align*} & \begin{bmatrix} 1 & 0 & 1 & 1 \\ 2 & 0 & 1 & 0 \\ -2& 3 & 4 & 0 \\ -5& 5 & 6 & 0 \end{bmatrix} & \cdots & \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \\ \to & \begin{bmatrix} 1 & 0 & 1 & 1 \\ 0 & 0 & -1 & -2 \\ 0 & 3 & 6 & 2 \\ 0 & 5 & 11 & 5 \end{bmatrix} & \cdots & \begin{bmatrix} 1 & 0 & 0 & 0 \\ -2 & 1 & 0 & 0 \\ 2 & 0 & 1 & 0 \\ 5 & 0 & 0 & 1 \end{bmatrix} \\ \to & \begin{bmatrix} 1 & 0 & 1 & 1 \\ 0 & 3 & 6 & 2 \\ 0 & 0 & -1 & -2 \\ 0 & 5 & 11 & 5 \end{bmatrix} & \cdots & \begin{bmatrix} 1 & 0 & 0 & 0 \\ 2 & 0 & 1 & 0 \\ -2 & 1 & 0 & 0 \\ 5 & 0 & 0 & 1 \end{bmatrix} \\ \to & \begin{bmatrix} 1 & 0 & 1 & 1 \\ 0 & 3 & 6 & 2 \\ 0 & 0 & -1 & -2 \\ 0 & 0 & 1 & 5/3 \end{bmatrix} & \cdots & \begin{bmatrix} 1 & 0 & 0 & 0 \\ 2 & 0 & 1 & 0 \\ -2 & 1 & 0 & 0 \\ 5/3 & 0 & -5/3 & 1 \end{bmatrix} \\ \to & \begin{bmatrix} 1 & 0 & 0 & 1 \\ 0 & 3 & 0 & -10 \\ 0 & 0 & -1 & -2 \\ 0 & 0 & 0 & -1/3 \end{bmatrix} & \cdots & \begin{bmatrix} -1 & 1 & 0 & 0 \\ -10 & 6 & 1 & 0 \\ -2 & 1 & 0 & 0 \\ -1/3 & 1 & -5/3 & 1 \end{bmatrix} \\ \to & \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 3 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & -1/3 \end{bmatrix} & \cdots & \begin{bmatrix} 0 & -2 & 5 & -3 \\ 0 & -24 & 51 & -30 \\ 0 & -5 & 10 & -6 \\ -1/3 & 1 & -5/3 & 1 \end{bmatrix} \\ \to & \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} & \cdots & \begin{bmatrix} 0 & -2 & 5 & -3 \\ 0 & -8 & 17 & -10 \\ 0 & 5 & -10 & 6 \\ 1 & -3 & 5 & -3 \end{bmatrix} \end{align*}$$

コード

以下は、ガウスの消去法を使って逆行列を求めるジュリアのコードだ。

using LinearAlgebra

M = [
      1 0 1 1
      2 0 1 0
     -2 3 4 0
     -5 5 6 0
     ]
inv_M = inv(M)

n = size(M)[1]
M = Float64.(M)
X = Matrix(1.0I,n,n)
for j = 1:n
    if iszero(M[j,j])
        for i = j:n
            if !iszero(M[i,i])
                M[[j,i],:] = M[[i,j],:]
                X[[j,i],:] = X[[i,j],:]
                break
            end
        end
    end

    pivot = M[j,j]
    for i = 1:n
        if i == j continue end
        r = (- M[i,j] / pivot)
        M[i,:] += r*M[j,:]
        X[i,:] += r*X[j,:]
    end
end

for i = 1:n
    pivot = M[i,i]
    M[i,:] /= pivot
    X[i,:] /= pivot
end
M
X
inv_M