YoJlのポケット

勉強したことや学んだことを記録しています。

実装を用いた行列の演算と基本変形(階数・逆行列・連立一次方程式)

想定する読者層

Python (numpy)の実装環境がある.
・行列の演算と基本変形をだいたい理解している.

実装環境

Google Colab
 使い方はGoogle Colabの使い方を参照ください.

1. 行列の演算

1.1. 行列の和と差

 手計算ですと, 以下のような定義が成り立ちます.

f:id:YoJl:20201215224205p:plain:w300

f:id:YoJl:20201215224402p:plain:w300

 上のように, 行列の和は要素ごとの足し算となっています.
 また, 同様に行列の差でも要素ごとの引き算となっています.

【行列の和と差の実装】
 まず, 2×2の正方行列を2つ生成します.

# 2×2の正方行列を2つ生成する
import numpy as np
A = np.arange(1,5).reshape((2,2))
B = np.arange(11,15).reshape((2,2))
print(A)
print(B)

出力
[[1 2]
 [3 4]]
[[11 12]
 [13 14]]

# AとBの足し算を行う
C = A + B
print(C)

出力
[[12 14]
 [16 18]]

# AとBの引き算を行う
C = A - B
print(C)

出力
[[-10 -10]
 [-10 -10]]

1.2. 行列のスカラー

 手計算ですと, 以下のような定義が成り立ちます.

f:id:YoJl:20201215224553p:plain:w300

 上のように, 要素一つ一つをスカラー(k)倍します.

【行列のスカラー倍の実装】
 まず, 2×2の正方行列を生成します.

# 2×2の正方行列を生成する
import numpy as np
A = np.arange(1,5).reshape((2,2))
print(A)

出力
[[1 2]
 [3 4]]

# Aをk倍する
k = 2
D = k * A
print(D)

出力
[[2 4]
 [6 8]]

1.3. 行列の積

 手計算ですと, 以下のような定義が成り立ちます.

f:id:YoJl:20201219224549p:plain:w550

 上のように, Aの行ベクトルとBの列ベクトルを用いて計算します.

【行列の積の実装】
 まず, 2×2の正方行列を2つ生成します.

# 2×2の正方行列を2つ生成する
import numpy as np
A = np.arange(1,5).reshape((2,2))
B = np.arange(11,15).reshape((2,2))
print(A)
print(B)

出力
[[1 2]
 [3 4]]
[[11 12]
 [13 14]]

# AとBの積算
E = A @ B
print(E)

出力
[[37 40]
 [85 92]]

ちなみに, 要素ごとの積算をしたい場合は, A * Bとします.

2. 行列の基本変形

 行列の基本変形を行うことで, 行列の階数, 逆行列などを求めることが可能です.
 行列の基本変形の操作の種類は以下の通りです.

① 行(列)を入れ替える
② 行(列)を定数倍する. ただし, 定数は0以外に限る.
③ ある行(列)に他の行(列)の定数倍を加える.

また,
行に対して行う基本変形を特に行の基本変形,
列に対して行う基本変形を特に列の基本変形といいます.
基本変形は, 線形代数の基礎となります.

さらに, 行の基本変形はどのような問題でも利用可能ですが, 列の基本変形は使える問題が限られます.
一般的には, 列の基本変形は行列の階数を求める問題のみで利用可能と覚えると良いと思います.
※厳密には, 他にも列の基本変形を用いることができる問題もありますが, 行の基本変形で代用可能です.

2.1. 行列の階数

 行列の基本変形を用いて, 標準形を作ります.
 このとき, 現れた1の個数が階数です.
階数のイメージとしては, 次元数のことです.

【行列の階数の実装】

# 2×2の正方行列を生成する
import numpy as np
# 0から9の3×3の行列を生成
A = np.random.randint(0, 10, (3, 3))
print(A)

出力
[[4 2 5]
 [3 1 0]]
 [6 7 1]]

# 階数を求める
print(np.linalg.matrix_rank(A))

出力
3

2.2. 逆行列

n次の正方行列Aについて,

AX = E, XA = E

このとき, XA逆行列といいます.
また, 逆行列をもつ行列を正則行列といいます.
逆行列の求め方は,

f:id:YoJl:20201223232719p:plain:w280

とおき, 行の基本変形を行うことで,

f:id:YoJl:20201223232909p:plain:w220

という形にします. したがって、

f:id:YoJl:20201223233051p:plain:w200

【行列の階数の実装】

# ある一次独立の行列を用意
import numpy as np
A = np.array([[1, 2, 3], [2, 3, 1], [2, 4, 5]])
print(A)

出力
[[1 2 3]
 [2 3 1]]
 [2 4 5]]

# 逆数を求める
print(np.linalg.inv(A))

出力
[[11. 2. -7.]
 [-8. -1. 5.]
 [ 2. 0. -1.]]

2.3. 行列の連立一次方程式
f:id:YoJl:20201221210505p:plain:w200

について,

f:id:YoJl:20201221210647p:plain:w450

とすると, A係数行列, b定数ベクトル, (Ab)係数拡大行列
といいます.

この係数拡大行列において,

f:id:YoJl:20201221210945p:plain:w500

となるように行の基本変形を行います. E単位行列です.
また, パラメータの個数が0のとき, 一次独立といい, そうでないときは, 一次従属といいます.
パラメータの個数は以下の方法で求めることができます.

パラメータの個数 = 未知数(x,yなど)の個数ーrank(階数)

【行列の連立一次方程式の実装~解が一つの値に定まる~】

import numpy as np
# ある係数拡大行列を用意
A = np.array([[2, 4], [3, 7]])
b = np.array([[16], [27]])
A_b = np.hstack([A, b])
print(A_b)

出力
[[2 4 16]
 [3 7 27]]

from numpy.linalg import solve
print(solve(A, b))

出力
[[2.]
 [3.]]

となり, x=2, y=3が得られます.
このようなx, yが一つに定まるような連立一次方程式は, numpyのsloveモジュールを用いることで解を得られます.
また、逆行列を用いて

Ax=b ⇔ A-1Ax=A-1b ⇔ x=A-1b

と求めることも可能です.

いずれにしても, 係数拡大行列が一次独立で逆行列が存在していてるもの(または, パラメータ数が0)に限られています.

【行列の連立一次方程式の実装~解が一つの値に定まらない~】
 次に, 解が一つに定まらない場合(一次従属, 逆行列が存在しない, パラメータ数>0)を考えます.
 手計算ですとパラメータ(tなど)を用いて以下のように示します.

f:id:YoJl:20201221211531p:plain:w300

実は, このような場合の便利なPythonモジュールは存在しません.
なぜならば, パラメータに置き換える文字によって答えが変わるためです.
ですので, 自分で関数を作って実装する必要があります.
しかし, 計算手段やその順番が関わってくるため, 機械にやらせるとなるととても長いコードを書く必要があります.
よって, 今回は全自動にせず, 手助け的な立ち位置で, 実装していきたと思います.

import numpy as np
# 一次従属の係数拡大行列を用意
A = np.array([[1, 2, 3], [2, 5, 10], [5, 12, 23]])
b = np.array([[4], [9], [22]])
A_b = np.hstack([A, b])
print(A_b)

出力
[[ 1 2 3 4]
 [ 2 5 10 9]
 [ 5 12 23 22]]

では, 基本変形を行うための関数を宣言します.

# 行の定数倍
def kakeru(row,k):
  return row*k

# ある行を他の行の定数倍加える
def tasu_row(row1, row2, k):
  return row1 + k*row2

以下, 係数拡大行列(Ab)の行ベクトルを順に, ①, ②, ③とします.

# ②-2×①
A_b[1,:] = tasu_row(A_b[1,:],A_b[0,:],-2)
print(A_b)

出力 [[ 1 2 3 4]
 [ 0 1 4 1]
 [ 5 12 23 22]]

# ③-5×①
A_b[2,:] = tasu_row(A_b[2,:],A_b[0,:],-5)
print(A_b)

出力
[[ 1 2 3 4]
 [ 0 1 4 1]
 [0 2 8 2]]

# ①-2×②と③×(1/2)
A_b[0,:] = tasu_row(A_b[0,:],A_b[1,:],-2)
A_b[2,:] = kakeru(A_b[2,:],(1/2))
print(A_b)

出力
[[ 1 0 -5 2]
 [ 0 1 4 1]
 [0 0 0 0]]

単位行列で, 零行列の形になりましたので, 基本変形は以上で終了です.
あとは, z = tとして, x =の形にすると,

f:id:YoJl:20201221211531p:plain:w300

となります.

以上, 実装を用いた行列の演算と基本変形(階数・逆行列・連立一次方程式)でした!

p.s. いつか、パラメータのある連立一次方程式を求められるようなスクリプトを書きたいとは思いますが...うーむ、手計算だと簡単な割に、スクリプトにするとなると、難しいような問題なので、そこまでやりたい!と思えないのが本音なのですよねぇ…うーむ。"(-""-)"

【参考文献】

目次の作り方
Google Colabの使い方
NumPyで階数(ランク)を求めるlinalg.matrix_rank関数の使い方
Python NumPy SciPy サンプルコード: 線形連立方程式, LU 分解
NumPyで逆行列を求めるlinalg.invの使い方