参考リンクまとめ
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()
結果