参考リンクまとめ

Pythonコード


import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import qr, eig, solve, svd, LinAlgError

# -----------------------------
# LU分解：連立一次方程式の解法（Ax = b）
# -----------------------------
A_lu = np.array([[3, 1], [1, 2]])
b = np.array([9, 8])

try:
    x_lu = solve(A_lu, b)
except LinAlgError:
    x_lu = None

# -----------------------------
# QR分解：最小二乗法による直線フィッティング
# -----------------------------
x_data = np.linspace(0, 10, 10)
y_data = 2 * x_data + 1 + np.random.normal(0, 1, len(x_data))  # y = 2x + 1 + noise

# A = [x 1] の形でQR分解
A_qr = np.vstack([x_data, np.ones(len(x_data))]).T
Q, R = qr(A_qr)
coeffs_qr = solve(R, Q.T @ y_data)  # 最小二乗解
y_fit = A_qr @ coeffs_qr

# -----------------------------
# 固有値・固有ベクトル（スペクトル分解）
# -----------------------------
A_eig = np.array([[2, -1], [-1, 2]])
eigvals, eigvecs = eig(A_eig)

# -----------------------------
# 特異値分解（SVD）
# -----------------------------
A_svd = np.array([[3, 1], [1, 3]])
U, S, VT = svd(A_svd)  # S: 特異値の配列（Σの対角成分）

# -----------------------------
# 可視化
# -----------------------------
fig, axs = plt.subplots(1, 4, figsize=(22, 5))

# LU分解の結果
axs[0].set_title("LU Decomposition (Solve Ax = b)")
axs[0].axis('off')
axs[0].text(0.1, 0.5, f"A = {A_lu.tolist()}\nb = {b.tolist()}\nx = {x_lu.tolist() if x_lu is not None else '解なし'}",
            fontsize=12)

# QR分解の最小二乗法プロット
axs[1].set_title("QR Decomposition (Least Squares Fit)")
axs[1].scatter(x_data, y_data, label='Data')
axs[1].plot(x_data, y_fit, color='r', label='Fitted Line')
axs[1].legend()
axs[1].set_xlabel('x')
axs[1].set_ylabel('y')

# 固有値・固有ベクトルの視覚化
axs[2].set_title("Spectral Decomposition")
origin = [0], [0]
for i in range(len(eigvals)):
    axs[2].quiver(*origin, eigvecs[0, i], eigvecs[1, i], angles='xy', scale_units='xy', scale=1, label=f"λ={eigvals[i]:.2f}")
axs[2].set_xlim(-2, 2)
axs[2].set_ylim(-2, 2)
axs[2].axhline(0, color='black', linewidth=0.5)
axs[2].axvline(0, color='black', linewidth=0.5)
axs[2].grid()
axs[2].legend()
axs[2].set_aspect('equal')

# 特異値分解（SVD）の視覚化
axs[3].set_title("SVD (Singular Value Decomposition)")
axs[3].bar([f"σ{i+1}" for i in range(len(S))], S, color='purple')
axs[3].set_ylabel('Singular Values')
axs[3].set_ylim(0, max(S)*1.2)
axs[3].grid(axis='y')

plt.tight_layout()
plt.show()
結果