線形代数の計算には NumPy を使うのが一般的だが、行列演算の原理を理解するために純粋な Python で実装してみよう。リストのリストで行列を表現し、基本的な演算を自前で書く。
行列の表現
行列は「リストのリスト」で表す。各要素のリストが行に対応する。
# 2×3 行列 A = [ [1, 2, 3], [4, 5, 6] ] # 要素へのアクセス(0 始まり) print(A[0][1]) # 2(1 行目 2 列目) print(A[1][2]) # 6(2 行目 3 列目)
行数と列数を取得する関数を用意しておくと便利だ。
def shape(A): """行列の形状 (行数, 列数) を返す""" return len(A), len(A[0]) A = [[1, 2, 3], [4, 5, 6]] print(shape(A)) # (2, 3)
ゼロ行列と単位行列
よく使う特殊な行列を生成する関数を作る。
def zeros(m, n): """m×n のゼロ行列を生成""" return [[0 for _ in range(n)] for _ in range(m)] def identity(n): """n×n の単位行列を生成""" return [[1 if i == j else 0 for j in range(n)] for i in range(n)] print(zeros(2, 3)) # [[0, 0, 0], [0, 0, 0]] print(identity(3)) # [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
行列の加算
同じサイズの行列同士を要素ごとに足す。
def add(A, B): """行列の加算 A + B""" m, n = shape(A) return [[A[i][j] + B[i][j] for j in range(n)] for i in range(m)] A = [[1, 2], [3, 4]] B = [[5, 6], [7, 8]] print(add(A, B)) # [[6, 8], [10, 12]]
スカラー倍
行列の全要素にスカラーを掛ける。
def scalar_mult(c, A): """スカラー倍 c * A""" m, n = shape(A) return [[c * A[i][j] for j in range(n)] for i in range(m)] A = [[1, 2], [3, 4]] print(scalar_mult(3, A)) # [[3, 6], [9, 12]]
行列の乗算
行列の積 は、 で定義される。 が 、 が のとき、 は になる。
def mult(A, B): """行列の乗算 A @ B""" m, n = shape(A) n2, p = shape(B) if n != n2: raise ValueError("行列のサイズが不一致") C = zeros(m, p) for i in range(m): for j in range(p): for k in range(n): C[i][j] += A[i][k] * B[k][j] return C A = [[1, 2], [3, 4]] B = [[5, 6], [7, 8]] print(mult(A, B)) # [[19, 22], [43, 50]]
この三重ループは の計算量を持つ。大きな行列では非効率だが、原理を理解するには十分だ。
転置
行と列を入れ替える。
def transpose(A): """転置行列 A^T""" m, n = shape(A) return [[A[i][j] for i in range(m)] for j in range(n)] A = [[1, 2, 3], [4, 5, 6]] print(transpose(A)) # [[1, 4], [2, 5], [3, 6]]
ベクトルと行列の積
ベクトルは 1 列の行列として扱える。列ベクトル と行列 の積 を計算する専用関数を作ると便利だ。
def matvec(A, x): """行列とベクトルの積 A @ x""" m, n = shape(A) if len(x) != n: raise ValueError("サイズ不一致") return [sum(A[i][j] * x[j] for j in range(n)) for i in range(m)] A = [[1, 2], [3, 4], [5, 6]] x = [1, 2] print(matvec(A, x)) # [5, 11, 17]
内積
ベクトルの内積(ドット積)も基本的な演算だ。
def dot(u, v): """ベクトルの内積 u · v""" return sum(u[i] * v[i] for i in range(len(u))) u = [1, 2, 3] v = [4, 5, 6] print(dot(u, v)) # 32(= 1*4 + 2*5 + 3*6)
行列式(2×2 と 3×3)
行列式は正方行列に対して定義されるスカラー値だ。2×2 と 3×3 の場合を実装する。
def det2(A): """2×2 行列の行列式""" return A[0][0] * A[1][1] - A[0][1] * A[1][0] def det3(A): """3×3 行列の行列式(サラスの公式)""" return (A[0][0] * A[1][1] * A[2][2] + A[0][1] * A[1][2] * A[2][0] + A[0][2] * A[1][0] * A[2][1] - A[0][2] * A[1][1] * A[2][0] - A[0][1] * A[1][0] * A[2][2] - A[0][0] * A[1][2] * A[2][1]) A = [[1, 2], [3, 4]] print(det2(A)) # -2 B = [[1, 2, 3], [4, 5, 6], [7, 8, 10]] print(det3(B)) # -3
一般の 行列の行列式は、余因子展開やLU分解を使って計算する。
逆行列(2×2)
2×2 行列の逆行列は公式で計算できる。
def inv2(A): """2×2 行列の逆行列""" d = det2(A) if d == 0: raise ValueError("逆行列が存在しない(特異行列)") return [[A[1][1] / d, -A[0][1] / d], [-A[1][0] / d, A[0][0] / d]] A = [[4, 7], [2, 6]] A_inv = inv2(A) print(A_inv) # [[0.6, -0.7], [-0.2, 0.4]] # 検算: A @ A^(-1) = I print(mult(A, A_inv)) # [[1.0, 0.0], [0.0, 1.0]]
純粋な Python
原理の理解に最適。小さな行列向き。大きな行列では遅い。
NumPy
高度に最適化されている。実務では必ず NumPy を使う。
純粋な Python で実装することで、行列演算のアルゴリズムを深く理解できる。次のステップとして、連立方程式の解法(ガウス消去法)を学ぶとよい。