行列計算の Python 実現ガイド

【概要】行列計算の基礎と実装について、NumPyライブラリを中心に解説する。行列の積、LU分解、逆行列、行列式、畳み込みなどの基本的な行列演算について説明し、Pythonによる実装例を示す。

【サイト内のPython関連主要ページ】

【外部リソース】

用語リスト

本ガイドで使用する用語を以下に定義する。

基礎概念

評価指標

PythonのNumPyライブラリを使用することで、効率的な行列計算を実行できる。

1. 行列の積

本セクションでは、行列の積の計算方法とNumPyによる実装を説明する。

行列の積C = A × Bにおいて、結果の各要素C[i][j]は、Aのi行目とBのj列目の対応する要素の積の総和として計算される。m×n行列Aとn×p行列Bの積は、m×p行列となる。

NumPyでは@演算子またはnp.dot()関数を使用する。内部では最適化されたBLAS（Basic Linear Algebra Subprograms：基本線形代数サブプログラム）が使用され、高速な演算が実現される。 

import numpy as np

A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
C = A @ B  # または C = np.dot(A, B)
print("行列A:\n", A)
print("行列B:\n", B)
print("A×B:\n", C)

実行結果： 

行列A:
 [[1 2]
 [3 4]]
行列B:
 [[5 6]
 [7 8]]
A×B:
 [[19 22]
 [43 50]]

上記の結果において、C[0][0] = 19は、1×5 + 2×7 = 19として計算される。同様に、C[1][1] = 50は、3×6 + 4×8 = 50として計算される。

2. LU分解と計算精度の評価

本セクションでは、LU分解の原理と実装、および計算精度の評価方法を説明する。

LU分解は行列Aを下三角行列（L）と上三角行列（U）の積に分解する数値計算手法である。実装では、数値安定性を高めるため、置換行列（P）も含めてPA = LUの形で分解される。

置換行列Pは行の入れ替えを表す行列である。LU分解において、対角成分に0または小さな値が現れると数値的に不安定になるため、行を入れ替えて対角成分に絶対値の大きな要素を配置する。この手法を部分ピボット選択と呼ぶ。

LU分解が連立方程式Ax = bを解くのに有用な理由は以下の通りである。

scipy.linalgモジュールのlu()関数を使用することで、LU分解を数値的に安定して効率的に計算できる。

以下のコードでは、LU分解の実行に加えて、誤差、残差、行列の階数を計算し、計算精度を評価する。 

import numpy as np
from scipy import linalg

# 3×3の行列を作成
A = np.array([[2, -1, 1],
              [4, 1, -1],
              [1, 1, 1]])

# LU分解
P, L, U = linalg.lu(A)
print("元の行列A:\n", A)
print("置換行列P:\n", P)
print("下三角行列L:\n", L)
print("上三角行列U:\n", U)

# LU分解の検証（誤差の測定）
print("\nP×L×U:\n", P @ (L @ U))
error_lu = np.max(np.abs(A - P @ (L @ U)))
print("元の行列との誤差:", error_lu)

# LU分解を使用した連立方程式の解法（残差の測定）
b = np.array([1, 2, 3])
x = linalg.lu_solve(linalg.lu_factor(A), b)
print("\n連立方程式Ax = bの解:", x)
residual = np.linalg.norm(A @ x - b)
print("残差 ||Ax - b||:", residual)

# 行列の階数
rank_A = np.linalg.matrix_rank(A)
print("行列Aの階数:", rank_A)

実行結果： 

元の行列A:
 [[ 2 -1  1]
 [ 4  1 -1]
 [ 1  1  1]]
置換行列P:
 [[0. 1. 0.]
 [1. 0. 0.]
 [0. 0. 1.]]
下三角行列L:
 [[1.   0.   0.  ]
 [0.5  1.   0.  ]
 [0.25 0.5  1.  ]]
上三角行列U:
 [[ 4.   1.  -1. ]
 [ 0.  -1.5  1.5]
 [ 0.   0.   0.5]]

P×L×U:
 [[ 2. -1.  1.]
 [ 4.  1. -1.]
 [ 1.  1.  1.]]
元の行列との誤差: 0.0

連立方程式Ax = bの解: [0.  0.5 2.5]
残差 ||Ax - b||: 0.0
行列Aの階数: 3

上記の結果において、置換行列Pは1行目と2行目を入れ替える操作を表している。これにより、Uの対角成分に絶対値の大きな要素（4, -1.5, 0.5）が配置され、数値的に安定した分解が実現される。

3. 逆行列と行列式

本セクションでは、逆行列と行列式の計算方法、およびその検証方法を説明する。

逆行列A⁻¹は、元の行列との積が単位行列となる行列である。線形変換の逆変換を表現する概念である。

逆行列は行列式が0でない場合（正則行列の場合）にのみ存在する。実用的な計算では、数値安定性を考慮して逆行列の直接計算を避けることが推奨される。

行列式には以下の用途がある。

以下のコードでは、逆行列と行列式を計算し、A×A⁻¹ = Iが成り立つことを誤差の測定によって確認する。 

import numpy as np

A = np.array([[4, 7], [2, 6]])

# 行列式
det_A = np.linalg.det(A)
print("行列A:\n", A)
print("行列式:", det_A)

# 逆行列
A_inv = np.linalg.inv(A)
print("\n逆行列:\n", A_inv)

# 検証: A×A^(-1) = I（誤差の測定）
print("A×A^(-1):\n", A @ A_inv)
error_inv = np.max(np.abs(A @ A_inv - np.eye(2)))
print("逆行列の精度（誤差）:", error_inv)

実行結果： 

行列A:
 [[4 7]
 [2 6]]
行列式: 10.000000000000002

逆行列:
 [[ 0.6 -0.7]
 [-0.2  0.4]]
A×A^(-1):
 [[1.0000000e+00 0.0000000e+00]
 [1.1102230e-16 1.0000000e+00]]
逆行列の精度（誤差）: 1.1102230246251565e-16

上記の結果において、行列式が10.000000000000002と表示され、誤差として1.11×10⁻¹⁶程度の値が現れている。これは浮動小数点演算の特性に起因する。コンピュータは実数を有限桁の2進数で近似的に表現するため、計算過程で微小な丸め誤差が蓄積される。この程度の誤差（10⁻¹⁵〜10⁻¹⁶）は倍精度浮動小数点数の計算機イプシロンに相当し、数値計算として正常な結果である。

4. 畳み込み

本セクションでは、畳み込み演算の定義と実装、および出力サイズの計算方法を説明する。

畳み込みは、信号処理や画像処理における基本的な演算である。特徴抽出やフィルタリングに広く用いられる。

2次元畳み込みは、入力画像f(x, y)とカーネルg(x, y)に対して次の式で定義される。

$$(f * g)(x, y) = \sum_{i} \sum_{j} f(i, j) \cdot g(x - i, y - j)$$

この演算では、カーネルを入力画像上でスライドさせながら、対応する要素の積の総和を計算する。

scipy.signal.convolve2d()関数のmodeパラメータには以下の選択肢があり、用途に応じて選択する。

出力サイズは入力サイズとカーネルサイズから以下の式で計算される。入力サイズを(H, W)、カーネルサイズを(Kh, Kw)とすると、各modeの出力サイズは次の通りである。 

import numpy as np
from scipy import signal

# 2D配列（画像を想定）
image = np.array([[1, 2, 3, 4],
                  [5, 6, 7, 8],
                  [9, 10, 11, 12],
                  [13, 14, 15, 16]])

# カーネル（エッジ検出フィルタ）
kernel = np.array([[-1, -1, -1],
                  [-1,  8, -1],
                  [-1, -1, -1]])

print("入力画像:\n", image)
print("\nカーネル（エッジ検出フィルタ）:\n", kernel)

# 異なるモードでの畳み込み演算
modes = ['full', 'same', 'valid']
for mode in modes:
    conv_result = signal.convolve2d(image, kernel, mode=mode)
    print(f"\n畳み込み結果（mode='{mode}'）:\n", conv_result)
    print(f"出力サイズ: {conv_result.shape}")

実行結果： 

入力画像:
 [[ 1  2  3  4]
 [ 5  6  7  8]
 [ 9 10 11 12]
 [13 14 15 16]]

カーネル（エッジ検出フィルタ）:
 [[-1 -1 -1]
 [-1  8 -1]
 [-1 -1 -1]]

畳み込み結果（mode='full'）:
 [[ -1  -3  -5  -7  -7  -4]
 [ -6 -16 -24 -32 -26 -18]
 [-14 -24   0   0 -16 -26]
 [-22 -24   0   0  -8 -34]
 [-27 -40 -16  -8 -18 -38]
 [-13 -27 -29 -31 -33 -16]]
出力サイズ: (6, 6)

畳み込み結果（mode='same'）:
 [[-16 -24 -32 -26]
 [-24   0   0 -16]
 [-24   0   0  -8]
 [-40 -16  -8 -18]]
出力サイズ: (4, 4)

畳み込み結果（mode='valid'）:
 [[0 0]
 [0 0]]
出力サイズ: (2, 2)

上記の例では、入力サイズ(4, 4)とカーネルサイズ(3, 3)に対して、'full'モードでは(4+3-1, 4+3-1) = (6, 6)、'same'モードでは(4, 4)、'valid'モードでは(4-3+1, 4-3+1) = (2, 2)となる。

これらのコードを実行するには、次のコマンドで必要なライブラリをインストールする。 

pip install numpy scipy