行列計算の Python 実現ガイド

用語リスト

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

基礎概念

• 行列の積（Matrix Multiplication）：2つの行列を乗算する演算
• LU分解（LU Decomposition）：行列を下三角行列Lと上三角行列Uの積に分解する手法
• 逆行列（Inverse Matrix）：元の行列との積が単位行列となる行列
• 行列式（Determinant）：正方行列に対して定義されるスカラー値
• 固有値・固有ベクトル（Eigenvalue / Eigenvector）：Av = λv を満たすスカラー λ とベクトル v
• 特異値分解（SVD：Singular Value Decomposition）：行列を A = UΣVᵀ の形に分解する手法
• 条件数（Condition Number）：入力の誤差が解にどれだけ増幅されるかを表す指標
• 畳み込み（Convolution）：信号処理や画像処理で用いられる演算
• 行列の階数（Rank）：線形独立（どのベクトルも他のベクトルの線形結合では表現できない関係）な行または列の最大数

評価指標

• 残差（Residual）：||Ax - b|| で表される値。連立方程式の解の精度を評価する指標
• 誤差（Error）：計算結果と理論値との差

本ガイドの演習を実行するには，次のコマンドで必要なライブラリをインストールする。

pip install numpy scipy

演習1．行列の積

行列の積 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：行列・ベクトル演算の標準ライブラリ）が使用される）。

手順

1. 2×2 の行列 A と B を作成する。
2. @演算子を用いて積 C = A @ B を計算する。
3. A，B，C を表示する。

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)

ヒント

• C[i][j] は，A の i 行目と B の j 列目の対応する要素の積の総和である。
• C[0][0] は A の0行目 [1, 2] と B の0列目 [5, 7] から計算する。

考察ポイント

• C[0][0] が 1×5 + 2×7 で得られる値と一致するか確認する。
• C[1][1] が 3×6 + 4×8 で得られる値と一致するか確認する。

演習2．LU分解と計算精度の評価

LU分解は，行列 A を下三角行列（L：対角成分より上側がすべて0の行列）と上三角行列（U：対角成分より下側がすべて0の行列）の積に分解する手法である。実装では，数値安定性を高めるため，置換行列（P：行の入れ替えを表す行列）も含めて A = PLU の形で分解する。分解の各段階で対角成分（ピボット）を割り算に用いるが，対角成分に0または絶対値の小さな値が現れると数値的に不安定になるため，行を入れ替えて，対象列の対角成分以下の要素のうち絶対値が最大の要素をピボットとして選択する。これを部分ピボット選択（partial pivoting）と呼ぶ。

scipy.linalg モジュールの lu()関数により，LU分解を計算できる。lu() は分解結果の確認に用い，連立方程式の求解には，分解と求解を分けて扱う lu_factor() と lu_solve() を用いる。

手順

1. 3×3 の行列 A を作成する。
2. linalg.lu() で LU分解し，P，L，U を表示する。
3. P @ L @ U を計算し，元の行列 A との誤差を測定する。
4. lu_factor() と lu_solve() で連立方程式 Ax = b を解き，残差を測定する。
5. matrix_rank() で行列の階数を求める。

import numpy as np
from scipy import linalg

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)

ヒント

• 置換行列 P は，どの行とどの行を入れ替えたかを表す。
• A の第1列で絶対値が最大の要素がどの行にあるかに注目する。
• 誤差は np.max(np.abs(...))，残差は np.linalg.norm(...) で測定する。

考察ポイント

• P がどの行を入れ替えたかを，A の第1列の値と対応づけて読み取る。
• P×L×U が元の行列 A に一致し，誤差が0に近い値（計算機イプシロン程度）になるか確認する。
• 求めた解 x を元の式に代入し，残差が0に近い値になるか確認する。
• 階数3が，この行列が正則（行列式が0でない）であることと整合するか確認する。

演習3．逆行列と行列式

逆行列 A⁻¹ は，元の行列との積が単位行列（対角成分が1で他がすべて0の正方行列）となる行列であり，線形変換の逆変換を表す。逆行列は行列式が0でない場合（正則行列の場合）にのみ存在する。連立方程式 Ax = b を解く場合は，逆行列を計算して A⁻¹b を求めるよりも，LU分解を用いる方法が計算コストおよび数値精度の点で優れている。本演習では，逆行列および行列式の計算手順を確認する目的で np.linalg.inv() を使用する。

手順

1. 2×2 の行列 A を作成する。
2. np.linalg.det() で行列式を計算する。
3. np.linalg.inv() で逆行列を計算する。
4. A @ A_inv を計算し，単位行列との誤差を測定する。

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)

ヒント

• 2×2 行列 [[a, b], [c, d]] の行列式は ad - bc である。
• np.eye(2) は 2×2 の単位行列を表す。
• 表示される行列式や誤差に，浮動小数点演算（実数を有限桁の2進数で近似的に表現して行う演算）に由来する微小な値が現れる。

考察ポイント

• 行列式が0でないこと（正則であること）を確認する。
• A×A⁻¹ が単位行列に近づき，誤差が計算機イプシロン（倍精度では約 2.22×10⁻¹⁶）程度になるか確認する。
• 行列式や誤差に現れる微小な値（10⁻¹⁵〜10⁻¹⁶程度）が，丸め誤差の蓄積によるものであり，数値計算として正常な範囲であることを読み取る。

演習4．連立一次方程式の求解

連立一次方程式 Ax = b を解く場合，逆行列を計算して A⁻¹b を求めるのではなく，np.linalg.solve()関数を用いる。solve() は内部で LU分解に基づく直接法を用いており，逆行列を経由するよりも計算コストおよび数値精度の点で優れている。演習2の lu_factor()／lu_solve() は，同じ係数行列 A に対して異なる b を繰り返し解く場合に用いる。

手順

1. 2×2 の行列 A とベクトル b を作成する。
2. np.linalg.solve() で連立方程式 Ax = b を解く。
3. 残差を測定する。

import numpy as np

A = np.array([[3, 2], [1, 2]])
b = np.array([2, 0])

# 連立方程式 Ax = b を解く
x = np.linalg.solve(A, b)
print("行列A:\n", A)
print("ベクトルb:", b)
print("解 x:", x)

# 検証（残差の測定）
residual = np.linalg.norm(A @ x - b)
print("残差 ||Ax - b||:", residual)

ヒント

• 残差は np.linalg.norm(...) で測定する。
• 得られた解 x を元の式 Ax に代入すると b に一致する。

考察ポイント

• 解 x を元の式に代入した結果が b と一致するか確認する。
• 残差が0に近い値になるか確認する。

演習5．固有値と固有ベクトル

正方行列 A に対して，Av = λv を満たすスカラー λ を固有値，0でないベクトル v を固有ベクトルと呼ぶ。固有値・固有ベクトルは，主成分分析，振動解析，グラフ解析などに応用される。一般の正方行列には np.linalg.eig()関数を用いる。行列が対称行列（A = Aᵀ）またはエルミート行列の場合は，np.linalg.eigh()関数を用いると，実数の固有値が昇順で得られ，数値的にも安定する。固有ベクトルの符号や全体の定数倍には自由度があり，環境により符号が異なる場合がある。

手順

1. 2×2 の対称行列 A を作成する。
2. np.linalg.eigh() で固有値と固有ベクトルを求める。
3. A @ eigenvectors と eigenvectors * eigenvalues の誤差を測定する。

import numpy as np

# 対称行列
A = np.array([[2, 1], [1, 2]])

# 固有値・固有ベクトル（対称行列なので eigh を使用）
eigenvalues, eigenvectors = np.linalg.eigh(A)
print("行列A:\n", A)
print("固有値:", eigenvalues)
print("固有ベクトル（各列が1つの固有ベクトル）:\n", eigenvectors)

# 検証: A v = λ v（誤差の測定）
error_eig = np.max(np.abs(A @ eigenvectors - eigenvectors * eigenvalues))
print("固有値分解の誤差:", error_eig)

ヒント

• 固有ベクトルは行列 eigenvectors の各列に対応する（1列目が1つ目の固有値に対応する固有ベクトル）。
• eigenvectors * eigenvalues は，各列の固有ベクトルに対応する固有値を掛けた結果（列ごとのスケーリング）を表す。
• eigh() の固有値は昇順で返される。

考察ポイント

• 固有値が昇順で返されているか確認する。
• A @ eigenvectors と eigenvectors * eigenvalues の誤差が0に近い値になり，Av = λv が成り立つことを確認する。

演習6．特異値分解（SVD）

特異値分解は，行列 A を A = UΣVᵀ の形に分解する手法である。Σ の対角成分を特異値と呼び，非負の値が降順で並ぶ。特異値分解は，行列の階数の判定，低ランク近似，擬似逆行列の計算などに用いられる。演習2で用いた np.linalg.matrix_rank() も，内部では特異値を用いて階数を求めている。np.linalg.svd()関数の返り値は U，特異値の配列 s，Vᵀ の3つである。U と Vᵀ の符号には自由度があり，環境により符号が異なる場合があるが，特異値 s は一意に定まる。

手順

1. 2×2 の行列 A を作成する。
2. np.linalg.svd() で特異値分解する。
3. 特異値の配列 s を対角行列 Σ に変換し，U @ Σ @ Vᵀ で再構成して誤差を測定する。

import numpy as np

A = np.array([[2, 0], [0, 3]])

# 特異値分解
U, s, Vt = np.linalg.svd(A)
print("行列A:\n", A)
print("特異値:", s)

# 検証: A = U Σ Vᵀ（誤差の測定）
Sigma = np.diag(s)
A_reconstructed = U @ Sigma @ Vt
print("再構成した行列 U Σ Vᵀ:\n", A_reconstructed)
error_svd = np.max(np.abs(A - A_reconstructed))
print("特異値分解の誤差:", error_svd)

ヒント

• np.diag(s) は，特異値の配列 s を対角成分とする対角行列 Σ を作る。
• 特異値は降順で返される。

考察ポイント

• 特異値が降順で返されているか確認する。
• U Σ Vᵀ が元の行列 A に一致し，誤差が0に近い値になるか確認する。

演習7．条件数

条件数は，連立方程式 Ax = b において，入力（A や b）の小さな誤差が解 x にどれだけ増幅されて伝わるかを表す指標である。条件数が1に近いほど数値的に安定（良条件）であり，値が大きいほど誤差が増幅されやすく不安定（悪条件）である。np.linalg.cond()関数の既定では2ノルムが用いられ，このとき条件数は「最大特異値 ÷ 最小特異値」に等しい。

手順

1. 演習6と同じ 2×2 の行列 A を作成する。
2. np.linalg.cond() で条件数を計算する。

import numpy as np

A = np.array([[2, 0], [0, 3]])

# 条件数（既定は2ノルム）
cond_A = np.linalg.cond(A)
print("行列A:\n", A)
print("条件数:", cond_A)

ヒント

• 条件数は最大特異値を最小特異値で割った値である。
• 演習6で求めた特異値（3 と 2）を用いて手計算と照合できる。

考察ポイント

• 条件数が，演習6で求めた最大特異値と最小特異値の比（3 ÷ 2）に一致するか確認する。
• 得られた値が1に近いことから，この行列が良条件であると読み取れるか確認する。

演習8．畳み込み

畳み込みは，信号処理や画像処理で用いられる演算であり，特徴抽出やフィルタリングに適用される。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() は上記の定義に従いカーネルを180度反転して計算する（要素を反転せずに計算する「相関」とは異なる）。本演習のカーネルは点対称のため，反転の有無で結果は変わらない。

modeパラメータには以下の選択肢がある。

• 'full'：カーネルが入力に1要素でも重なる範囲をすべて計算する（出力は入力よりも大きくなる）
• 'same'：出力が入力と同じサイズになるように中央部分を取り出す
• 'valid'：カーネルが入力に完全に収まる位置のみを計算する

入力サイズを (H, W)，カーネルサイズを (Kh, Kw) とすると，各 mode の出力サイズは次の通りである。

• 'full'：(H + Kh - 1, W + Kw - 1)
• 'same'：(H, W)
• 'valid'：(H - Kh + 1, W - Kw + 1)

手順

1. 4×4 の入力画像と 3×3 のカーネルを作成する。
2. 'full'，'same'，'valid' の各モードで convolve2d() を実行する。
3. 各モードの畳み込み結果と出力サイズを表示する。

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}")

ヒント

• 入力サイズは (4, 4)，カーネルサイズは (3, 3) である。
• 各モードの出力サイズは，上記の出力サイズの式に値を代入して求められる。

考察ポイント

• 各モードの出力サイズが，出力サイズの式（'full' は (4+3-1, 4+3-1)，'same' は (4, 4)，'valid' は (4-3+1, 4-3+1)）と一致するか確認する。
• モードによって出力サイズと取り出される範囲がどう変わるかを読み取る。