3×3行列の固有値分解¶
In [14]:
# !pip install japanize-matplotlib
In [15]:
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
In [16]:
from mpl_toolkits.mplot3d import Axes3D
固有値
In [3]:
# 行列
A = np.array([[ 0., 0., 1.],
[ 0., 0., -1.],
[ 1., -1., 1.]])
In [4]:
# 固有値・固有ベクトル（対称行列なので np.linalg.eigh を使用）
w, V = np.linalg.eigh(A)
In [5]:
print("固有値 w =", w)
固有値 w = [-1.00000000e+00 1.08691268e-16 2.00000000e+00]
固有ベクトル
In [7]:
# 正規化済みの固有ベクトル
order = [1, 2, 0] # w[1]=0, w[2]=2, w[0]=-1
lam = w[order]
P_raw = V[:, order]
In [8]:
print("固有値（並べ替え後） lam =", lam)
固有値（並べ替え後） lam = [ 1.08691268e-16 2.00000000e+00 -1.00000000e+00]
直交行列・対角行列¶
In [9]:
# 列ごとに正規化（eigh の列は既に単位長のはずだが、念のため）
P = P_raw / np.linalg.norm(P_raw, axis=0, keepdims=True)
D = np.diag(lam)
In [10]:
print("\n直交行列 P =\n", np.round(P, 6))
print("\n対角行列 D =\n", D)
直交行列 P = [[-0.707107 0.408248 -0.57735 ] [-0.707107 -0.408248 0.57735 ] [-0. 0.816497 0.57735 ]] 対角行列 D = [[ 1.08691268e-16 0.00000000e+00 0.00000000e+00] [ 0.00000000e+00 2.00000000e+00 0.00000000e+00] [ 0.00000000e+00 0.00000000e+00 -1.00000000e+00]]
In [12]:
# 検算：P^T P = I, A = P D P^T
I_check = P.T @ P
recon = P @ D @ P.T
In [13]:
print("\nP^T P =\n", np.round(I_check, 6))
print("\nP D P^T =\n", np.round(recon, 6))
print("\nA =\n", A)
P^T P = [[ 1. 0. -0.] [ 0. 1. -0.] [-0. -0. 1.]] P D P^T = [[ 0. -0. 1.] [-0. 0. -1.] [ 1. -1. 1.]] A = [[ 0. 0. 1.] [ 0. 0. -1.] [ 1. -1. 1.]]
グラフ化¶
In [23]:
# 単位球メッシュ
u = np.linspace(0, 2*np.pi, 40)
v = np.linspace(0, np.pi, 20)
xs = np.outer(np.cos(u), np.sin(v))
ys = np.outer(np.sin(u), np.sin(v))
zs = np.outer(np.ones_like(u), np.cos(v))
S = np.stack([xs, ys, zs], axis=0).reshape(3, -1) # (3, N)
# A を作用
E = (A @ S).reshape(3, xs.shape[0], xs.shape[1])
xe, ye, ze = E[0], E[1], E[2]
# 固有ベクトル（列）
u0, u2, um1 = P[:, 0], P[:, 1], P[:, 2] # λ=0,2,-1 の順
lam0, lam2, lamm1 = lam[0], lam[1], lam[2]
fig = plt.figure(figsize=(12, 6))
# 左：単位球と楕円体
ax1 = fig.add_subplot(1, 2, 1, projection='3d')
ax1.plot_wireframe(xs, ys, zs, color='0.6', linewidth=0.5, rstride=2, cstride=2, label='単位球')
ax1.plot_wireframe(xe, ye, ze, color='C3', linewidth=0.8, rstride=2, cstride=2)
# 固有ベクトル（主軸）
ax1.quiver(0, 0, 0, u0[0], u0[1], u0[2], color='C0', linewidth=2)
ax1.quiver(0, 0, 0, u2[0], u2[1], u2[2], color='C1', linewidth=2)
ax1.quiver(0, 0, 0, um1[0], um1[1], um1[2], color='C2', linewidth=2)
ax1.set_title("単位球 → Aで楕円体に（主軸＝固有ベクトル）")
ax1.set_box_aspect([1,1,1])
# 右：固有ベクトルと λv の伸び
ax2 = fig.add_subplot(1, 2, 2, projection='3d')
def draw3(ax, v, label, color='C0', scale=1.0, ls='-'):
ax.quiver(0, 0, 0, *(scale*v), color=color, linewidth=2)
ax.text(*(1.05*scale*v), label, color=color)
draw3(ax2, u0, f"u0 (λ={lam0:.0f})", color='C0')
draw3(ax2, lam0*u0, "λ u0", color='C0')
draw3(ax2, u2, f"u2 (λ={lam2:.0f})", color='C1')
draw3(ax2, lam2*u2, "λ u2", color='C1')
draw3(ax2, um1, f"u-1 (λ={lamm1:.0f})", color='C2')
draw3(ax2, lamm1*um1, "λ u-1", color='C2')
# パネルに P, D の数値
panel = "P =\n" + np.array2string(np.round(P, 3)) + "\n\nD =\n" + np.array2string(D)
ax2.text2D(0.98, 0.98, panel, transform=ax2.transAxes,
ha='right', va='top', bbox=dict(facecolor='white', alpha=0.85, edgecolor='0.6'))
ax2.set_title("固有ベクトルと 伸び（λv）")
ax2.set_box_aspect([1,1,1])
# 表示範囲を少し広めに
all_pts = np.column_stack([u0, u2, um1, lam0*u0, lam2*u2, lamm1*um1, E.reshape(3,-1)])
m = 1.15*np.max(np.abs(all_pts))
for ax in (ax1, ax2):
ax.set_xlim(-m, m); ax.set_ylim(-m, m); ax.set_zlim(-m, m)
ax.set_xlabel('X'); ax.set_ylabel('Y'); ax.set_zlabel('Z')
plt.tight_layout()
plt.show()
In [ ]: