第3章: 固有値・固有ベクトルと対角化

import numpy as np import matplotlib.pyplot as plt # 2×2行列 A = np.array([[4, 1], [2, 3]]) # 固有値・固有ベクトルの計算 eigenvalues, eigenvectors = np.linalg.eig(A) print("固有値・固有ベクトルの計算:") print(f"A =\n{A}\n") print(f"固有値: {eigenvalues}") print(f"固有ベクトル:\n{eigenvectors}\n") # 検証: Av = λv for i in range(len(eigenvalues)): lam = eigenvalues[i] v = eigenvectors[:, i] Av = A @ v lam_v = lam * v print(f"λ{i+1} = {lam:.4f}") print(f" Av = {Av}") print(f" λv = {lam_v}") print(f" Av = λv? {np.allclose(Av, lam_v)}\n")

# 固有ベクトルの可視化 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6)) # ランダムなベクトルと固有ベクトル random_vec = np.array([1, 1]) eigen_vec1 = eigenvectors[:, 0] eigen_vec2 = eigenvectors[:, 1] # 変換前のベクトル origin = [0, 0] ax1.quiver(*origin, *random_vec, angles='xy', scale_units='xy', scale=1, color='blue', width=0.01, label='一般ベクトル') ax1.quiver(*origin, *eigen_vec1, angles='xy', scale_units='xy', scale=1, color='red', width=0.01, label=f'固有ベクトル1 (λ={eigenvalues[0]:.2f})') ax1.quiver(*origin, *eigen_vec2, angles='xy', scale_units='xy', scale=1, color='green', width=0.01, label=f'固有ベクトル2 (λ={eigenvalues[1]:.2f})') ax1.set_xlim(-1, 5) ax1.set_ylim(-1, 5) ax1.set_xlabel('x') ax1.set_ylabel('y') ax1.set_title('変換前') ax1.legend() ax1.grid(True, alpha=0.3) ax1.axhline(y=0, color='k', linewidth=0.5) ax1.axvline(x=0, color='k', linewidth=0.5) # 変換後のベクトル random_transformed = A @ random_vec eigen_transformed1 = A @ eigen_vec1 eigen_transformed2 = A @ eigen_vec2 ax2.quiver(*origin, *random_transformed, angles='xy', scale_units='xy', scale=1, color='blue', width=0.01, label='一般ベクトル（方向変化）', alpha=0.7) ax2.quiver(*origin, *eigen_transformed1, angles='xy', scale_units='xy', scale=1, color='red', width=0.01, label='固有ベクトル1（方向不変）') ax2.quiver(*origin, *eigen_transformed2, angles='xy', scale_units='xy', scale=1, color='green', width=0.01, label='固有ベクトル2（方向不変）') # 元のベクトルを薄く表示 ax2.quiver(*origin, *random_vec, angles='xy', scale_units='xy', scale=1, color='blue', width=0.005, alpha=0.2) ax2.quiver(*origin, *eigen_vec1, angles='xy', scale_units='xy', scale=1, color='red', width=0.005, alpha=0.2) ax2.quiver(*origin, *eigen_vec2, angles='xy', scale_units='xy', scale=1, color='green', width=0.005, alpha=0.2) ax2.set_xlim(-1, 5) ax2.set_ylim(-1, 5) ax2.set_xlabel('x') ax2.set_ylabel('y') ax2.set_title('行列Aによる変換後') ax2.legend(fontsize=8) ax2.grid(True, alpha=0.3) ax2.axhline(y=0, color='k', linewidth=0.5) ax2.axvline(x=0, color='k', linewidth=0.5) plt.tight_layout() plt.show()

import sympy as sp # 記号変数 lam = sp.Symbol('lambda') # 行列Aを記号で定義 A_sym = sp.Matrix([[4, 1], [2, 3]]) # I - 単位行列 I = sp.eye(2) # 特性行列 A - λI char_matrix = A_sym - lam * I print("特性方程式の導出:") print(f"A - λI =") sp.pprint(char_matrix) # 特性方程式 det(A - λI) = 0 char_poly = char_matrix.det() print(f"\n特性多項式: {char_poly} = 0") # 固有値を求める eigenvals_sym = sp.solve(char_poly, lam) print(f"固有値: {eigenvals_sym}") # NumPyの結果と比較 print(f"\nNumPyの固有値: {eigenvalues}")

# 固有値を対角成分に持つ対角行列D D = np.diag(eigenvalues) # 固有ベクトルを列に持つ行列P P = eigenvectors # P^(-1) P_inv = np.linalg.inv(P) # 検証: P^(-1) A P = D result = P_inv @ A @ P print("対角化の検証:") print(f"P (固有ベクトル行列) =\n{P}\n") print(f"D (固有値対角行列) =\n{D}\n") print(f"P^(-1) A P =\n{result}\n") print(f"D に一致? {np.allclose(result, D)}") # 対角化の応用: A^n の高速計算 # A^n = P D^n P^(-1) n = 10 A_n_fast = P @ np.linalg.matrix_power(D, n) @ P_inv A_n_direct = np.linalg.matrix_power(A, n) print(f"\nA^{n} の計算:") print(f"対角化利用: {A_n_fast[0,0]:.4f}") print(f"直接計算: {A_n_direct[0,0]:.4f}") print(f"一致? {np.allclose(A_n_fast, A_n_direct)}")

# 対称行列 A_sym = np.array([[3, 1], [1, 3]]) # 固有値・固有ベクトル eigenvals_sym, eigenvecs_sym = np.linalg.eigh(A_sym) # 対称行列専用 print("対称行列の固有値分解:") print(f"A (対称) =\n{A_sym}\n") print(f"固有値: {eigenvals_sym}") print(f"固有ベクトル:\n{eigenvecs_sym}\n") # 固有ベクトルの直交性を確認 v1 = eigenvecs_sym[:, 0] v2 = eigenvecs_sym[:, 1] inner_product = np.dot(v1, v2) print(f"固有ベクトルの内積: {inner_product:.10f}") print(f"直交? {np.abs(inner_product) < 1e-10}") # 直交行列性の確認: Q^T Q = I Q = eigenvecs_sym QTQ = Q.T @ Q print(f"\nQ^T Q =\n{QTQ}") print(f"単位行列? {np.allclose(QTQ, np.eye(2))}")

from sklearn.datasets import make_blobs # サンプルデータ生成 np.random.seed(42) X, _ = make_blobs(n_samples=100, n_features=2, centers=1, cluster_std=2.0) # データを中心化 X_centered = X - np.mean(X, axis=0) # 共分散行列 cov_matrix = np.cov(X_centered.T) # 固有値分解 eigenvals_pca, eigenvecs_pca = np.linalg.eigh(cov_matrix) # 固有値の大きい順にソート idx = eigenvals_pca.argsort()[::-1] eigenvals_pca = eigenvals_pca[idx] eigenvecs_pca = eigenvecs_pca[:, idx] print("PCA（主成分分析）:") print(f"共分散行列:\n{cov_matrix}\n") print(f"固有値: {eigenvals_pca}") print(f"寄与率: {eigenvals_pca / eigenvals_pca.sum()}") # 可視化 plt.figure(figsize=(10, 8)) plt.scatter(X_centered[:, 0], X_centered[:, 1], alpha=0.5) # 主成分方向を描画 origin = [0, 0] scale = np.sqrt(eigenvals_pca) for i in range(2): vec = eigenvecs_pca[:, i] * scale[i] * 3 plt.quiver(*origin, *vec, angles='xy', scale_units='xy', scale=1, color=['red', 'blue'][i], width=0.01, label=f'PC{i+1} (λ={eigenvals_pca[i]:.2f})') plt.xlabel('Feature 1') plt.ylabel('Feature 2') plt.title('主成分分析（PCA）') plt.legend() plt.grid(True, alpha=0.3) plt.axis('equal') plt.show()

# 第1主成分のみに射影（2D → 1D） PC1 = eigenvecs_pca[:, 0].reshape(-1, 1) # 射影 X_pca = X_centered @ PC1 # 元の空間に戻す（近似再構成） X_reconstructed = X_pca @ PC1.T # 再構成誤差 reconstruction_error = np.mean(np.linalg.norm(X_centered - X_reconstructed, axis=1)**2) print(f"\n次元削減:") print(f"元の次元: {X_centered.shape[1]}D") print(f"削減後: 1D") print(f"再構成誤差: {reconstruction_error:.4f}") print(f"第1主成分の寄与率: {eigenvals_pca[0]/eigenvals_pca.sum()*100:.1f}%") # 可視化 plt.figure(figsize=(10, 8)) plt.scatter(X_centered[:, 0], X_centered[:, 1], alpha=0.3, label='元データ') plt.scatter(X_reconstructed[:, 0], X_reconstructed[:, 1], alpha=0.5, label='再構成データ', s=10) plt.plot([0, PC1[0,0]*scale[0]*3], [0, PC1[1,0]*scale[0]*3], 'r-', linewidth=3, label='PC1') plt.xlabel('Feature 1') plt.ylabel('Feature 2') plt.title('第1主成分への射影と再構成') plt.legend() plt.grid(True, alpha=0.3) plt.axis('equal') plt.show()

# 2自由度連成振動系 # m1 = m2 = 1 kg, k1 = k2 = k3 = 1 N/m # 運動方程式: M x'' + K x = 0 # 固有値問題: det(K - ω² M) = 0 M = np.array([[1, 0], [0, 1]]) # 質量行列 K = np.array([[2, -1], [-1, 2]]) # 剛性行列 # 一般化固有値問題を解く eigenvals_vibration, eigenvecs_vibration = np.linalg.eigh(K, M) # 固有角振動数 ω = sqrt(λ) omega = np.sqrt(eigenvals_vibration) print("振動モード解析:") print(f"質量行列 M:\n{M}\n") print(f"剛性行列 K:\n{K}\n") print(f"固有値 λ: {eigenvals_vibration}") print(f"固有角振動数 ω: {omega} rad/s") print(f"固有振動数 f: {omega/(2*np.pi)} Hz\n") print("振動モード（固有ベクトル）:") for i in range(2): print(f"モード{i+1} (f={omega[i]/(2*np.pi):.3f} Hz):") print(f" 質点1: {eigenvecs_vibration[0,i]:.4f}") print(f" 質点2: {eigenvecs_vibration[1,i]:.4f}") if eigenvecs_vibration[0,i] * eigenvecs_vibration[1,i] > 0: print(f" → 同位相振動\n") else: print(f" → 逆位相振動\n")