行列演算を純粋な Python で実装する｜Python

線形代数の計算には 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]]

行列の乗算

行列の積 $C = A B$ は、 $C_{ij} = \sum_{k} A_{ik} B_{k j}$ で定義される。 $A$ が $m \times n$ 、 $B$ が $n \times p$ のとき、 $C$ は $m \times p$ になる。

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]]

この三重ループは $O (mn p)$ の計算量を持つ。大きな行列では非効率だが、原理を理解するには十分だ。

転置

行と列を入れ替える。

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 列の行列として扱える。列ベクトル $x$ と行列 $A$ の積 $A x$ を計算する専用関数を作ると便利だ。

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

一般の $n \times n$ 行列の行列式は、余因子展開や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 で実装することで、行列演算のアルゴリズムを深く理解できる。次のステップとして、連立方程式の解法（ガウス消去法）を学ぶとよい。