In [51]:
# !pip install japanize-matplotlib
In [52]:
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
In [53]:
# 行列 A を定義（今回の例は 2×3）
A = np.array([[1., 0., 1.],
[-1., 1., 0.]])
print("A =\n", A)
A = [[ 1. 0. 1.] [-1. 1. 0.]]
In [54]:
U, s, Vt = np.linalg.svd(A, full_matrices=False)
V = Vt.T
print("U=\n", U)
print("s=\n", s)
print("Vt=\n", Vt)
# ✅ 再構成（形が合う正しいやり方）
A_reconstructed = U @ np.diag(s) @ Vt
print("\nA_reconstructed=\n", A_reconstructed)
print("recon_error =", np.linalg.norm(A - A_reconstructed))
U= [[-0.70710678 0.70710678] [ 0.70710678 0.70710678]] s= [1.73205081 1. ] Vt= [[-8.16496581e-01 4.08248290e-01 -4.08248290e-01] [-1.19879883e-16 7.07106781e-01 7.07106781e-01]] A_reconstructed= [[ 1.00000000e+00 -1.21168839e-16 1.00000000e+00] [-1.00000000e+00 1.00000000e+00 4.47411937e-16]] recon_error = 6.882618060940576e-16
In [55]:
# 2) 固有値分解からSVDを求める流れ
print("\n--- 固有値分解から作るSVD ---")
# AA^T を作る
M = A @ A.T
lam, U2 = np.linalg.eigh(M) # 固有値分解（昇順）
s2 = np.sqrt(np.clip(lam, 0, None))
# 大きい順に並べ替え
idx = np.argsort(s2)[::-1]
s2, U2 = s2[idx], U2[:, idx]
# Σ の再構成
Sigma2 = np.zeros_like(A, dtype=float)
np.fill_diagonal(Sigma2, s2)
# V の計算：v_i = (A^T u_i) / σ_i
V2 = np.zeros((A.shape[1], A.shape[1]))
for i in range(len(s2)):
if s2[i] > 1e-12: # 特異値が0でない場合
vi = (A.T @ U2[:, i]) / s2[i]
V2[:, i] = vi / np.linalg.norm(vi)
print("U =\n", U2)
print("Sigma =\n", Sigma2)
print("V =\n", V2)
--- 固有値分解から作るSVD --- U = [[-0.70710678 -0.70710678] [ 0.70710678 -0.70710678]] Sigma = [[1.73205081 0. 0. ] [0. 1. 0. ]] V = [[-0.81649658 0. 0. ] [ 0.40824829 -0.70710678 0. ] [-0.40824829 -0.70710678 0. ]]
In [56]:
# 再構成誤差
A_reconstructed2 = U2 @ Sigma2 @ V2.T
print("\n再構成 A =\n", A_reconstructed2)
print("再構成誤差 =", np.linalg.norm(A - A_reconstructed2))
再構成 A = [[ 1.00000000e+00 -1.01465364e-17 1.00000000e+00] [-1.00000000e+00 1.00000000e+00 -1.01465364e-17]] 再構成誤差 = 1.434936932798653e-17
In [57]:
# 特異値（期待は √3 ≈ 1.732, と 1.0）
print("\n特異値 =", s)
特異値 = [1.73205081 1. ]
In [58]:
# --- プロット ---
fig, ax = plt.subplots(1, 3, figsize=(18,6))
# U のベクトル
origin = np.zeros(2)
colors = ["red","blue"]
for i in range(U.shape[1]):
vec = U[:, i]
ax[0].quiver(*origin, vec[0], vec[1], angles='xy', scale_units='xy',
scale=1, color=colors[i])
ax[0].text(vec[0]*1.1, vec[1]*1.1, f"u{i+1}", fontsize=12, color=colors[i])
ax[0].axhline(0, color="gray", linewidth=0.5)
ax[0].axvline(0, color="gray", linewidth=0.5)
ax[0].set_xlim(-1,1); ax[0].set_ylim(-1,1)
ax[0].set_aspect("equal")
ax[0].set_title("左特異ベクトル U")
# V のベクトル（第1,2列だけ）
colors = ["red","blue","green"]
for i in range(V.shape[1]):
vec = V[:, i][:2] # 2次元に射影
ax[1].quiver(*origin, vec[0], vec[1], angles='xy', scale_units='xy',
scale=1, color=colors[i])
ax[1].text(vec[0]*1.1, vec[1]*1.1, f"v{i+1}", fontsize=12, color=colors[i])
ax[1].axhline(0, color="gray", linewidth=0.5)
ax[1].axvline(0, color="gray", linewidth=0.5)
ax[1].set_xlim(-1,1); ax[1].set_ylim(-1,1)
ax[1].set_aspect("equal")
ax[1].set_title("右特異ベクトル V（PC1-PC2平面）")
# Σ の効果を楕円で描く
theta = np.linspace(0, 2*np.pi, 200)
circle = np.stack([np.cos(theta), np.sin(theta)], axis=0) # 単位円
ellipse = np.diag(s[:2]) @ circle # Σで伸縮
ax[2].plot(circle[0], circle[1], 'k--', label="単位円")
ax[2].plot(ellipse[0], ellipse[1], 'r-', label="Σによる伸縮(楕円)")
ax[2].set_aspect("equal")
ax[2].axhline(0, color="gray", linewidth=0.5)
ax[2].axvline(0, color="gray", linewidth=0.5)
ax[2].legend()
ax[2].set_title("Σ による伸縮")
plt.show()
SVD（回転 → 伸縮 → 回転）をアニメで可視化（2×2行列）¶
In [59]:
# SVD（回転 → 伸縮 → 回転）をアニメで可視化（2×2行列）
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
# ===== 1) 例として2×2行列を用意（好きな行列に変えてOK） =====
# 元の 2x3 行列 A の“雰囲気”に近い 2x2 を例示
M = np.array([[1.0, 1.0],
[-1.0, 1.0]])
# ===== 2) SVD: M = U Σ R^T =====
U, s, Rt = np.linalg.svd(M, full_matrices=False)
R = Rt.T
# 直交行列の向きを揃えて det=+1 に（反射が入ると回転補間が難しいため）
def fix_to_rotation(Q, Sigma):
if np.linalg.det(Q) < 0:
Q[:, 0] *= -1
Sigma[0] *= -1
return Q, Sigma
U, s = fix_to_rotation(U, s.copy())
R, s = fix_to_rotation(R, s.copy())
Sigma = np.diag(s)
# ===== 3) 2×2直交行列を“角度”で補間する関数 =====
def angle_from_rot2(Q):
# 2x2回転行列 Q = [[c,-s],[s,c]] の角度 θ を返す
return np.arctan2(Q[1,0], Q[0,0])
def rot2(theta):
c, s = np.cos(theta), np.sin(theta)
return np.array([[c, -s],
[s, c]])
theta_R = angle_from_rot2(R)
theta_U = angle_from_rot2(U)
# ===== 4) 伸縮（Σ）の補間：t=0 → I、t=1 → Σ =====
def sigma_interp(t, svals):
# 対角を 1 -> s_i に線形補間
return np.diag(1 + t*(svals - 1))
# ===== 5) 表示用の図形（単位円＆正方格子）を準備 =====
# 単位円
t = np.linspace(0, 2*np.pi, 200)
circle = np.stack([np.cos(t), np.sin(t)], axis=0) # shape=(2,N)
# 正方格子
grid_x = np.linspace(-1, 1, 11)
grid_y = np.linspace(-1, 1, 11)
H = []
for gx in grid_x:
H.append(np.stack([np.full_like(grid_y, gx), grid_y], axis=0))
for gy in grid_y:
H.append(np.stack([grid_x, np.full_like(grid_x, gy)], axis=0))
# ===== 6) アニメ設定 =====
fig, ax = plt.subplots(figsize=(5,5))
ax.set_aspect('equal', adjustable='box')
ax.set_xlim(-3, 3); ax.set_ylim(-3, 3)
ax.set_title("SVD Animation: rotate (V^T) → stretch (Σ) → rotate (U)")
circle_line, = ax.plot([], [], lw=2)
grid_lines = [ax.plot([], [], lw=0.5, alpha=0.6)[0] for _ in H]
txt_phase = ax.text(0.02, 0.98, "", transform=ax.transAxes, va='top')
# 1→3 の3フェーズ（各 0〜1 の補間）
def current_transform(t_all):
# t_all ∈ [0,1)
phase = int(3 * t_all) # 0:V^T, 1:Sigma, 2:U
local = 3 * t_all - phase # フェーズ内の 0→1
if phase == 0:
# V^T ~ R^T を 0→θ_R で補間
T = rot2(local * theta_R).T # R^T の補間
label = "Phase 1: rotate by V^T"
elif phase == 1:
# Σ を I→Σ に補間（回転は固定で R^T 済）
T = sigma_interp(local, s) # 伸縮のみ
label = "Phase 2: stretch by Σ"
else:
# U を 0→θ_U で補間（最後に掛ける）
T = rot2(local * theta_U)
label = "Phase 3: rotate by U"
return phase, local, T, label
def full_map(t_all, X):
# X: shape=(2,N)
# 全体の変換： (Phase1: R^T) → (Phase2: Σ) → (Phase3: U)
# それぞれを部分的に適用していく
phase = int(3 * t_all)
local = 3 * t_all - phase
# 1) R^T を全適用
Y = R.T @ X
if phase == 0:
# まだ途中
Y = rot2(local * theta_R).T @ X
return Y
# 2) Σ を適用
Y = Sigma @ Y
if phase == 1:
# 途中補間
Y = (sigma_interp(local, s) @ (R.T @ X))
return Y
# 3) U を適用
Y = U @ Y
if phase == 2:
Y = (rot2(local * theta_U) @ (Sigma @ (R.T @ X)))
return Y
return Y
def init():
circle_line.set_data([], [])
for gl in grid_lines:
gl.set_data([], [])
txt_phase.set_text("")
return [circle_line, *grid_lines, txt_phase]
def update(frame):
t_all = frame / 180 # 0→1
phase, local, Tlocal, label = current_transform(t_all)
# 図形を変換
C = full_map(t_all, circle)
circle_line.set_data(C[0], C[1])
for i, seg in enumerate(H):
S = full_map(t_all, seg)
grid_lines[i].set_data(S[0], S[1])
txt_phase.set_text(label)
return [circle_line, *grid_lines, txt_phase]
anim = FuncAnimation(fig, update, frames=180, init_func=init, interval=30, blit=True)
HTML(anim.to_jshtml()) # Colabで再生
Output hidden; open in https://colab.research.google.com to view.
SVD（回転 → 伸縮 → 回転）をアニメで可視化（2×3行列）¶
In [60]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
# ======================
# 1) 2×3 の元の行列
# ======================
A = np.array([[1., 0., 1.],
[-1., 1., 0.]])
# SVD（rank=2 を想定）
U, s, Vt = np.linalg.svd(A, full_matrices=False)
V = Vt.T # V: (3,2), U: (2,2), s: (2,)
Sigma = np.diag(s)
# 右特異ベクトルが張る 2D 平面の基底（3D中の 2 本の直交ベクトル）
# 任意の z∈R^2 は x = V z として 3D に埋め込まれる
# ======================
# 2) 可視化用の入力図形（2Dの z 座標）を準備
# ======================
# 単位円と格子（z平面上）
t = np.linspace(0, 2*np.pi, 200)
circle_z = np.stack([np.cos(t), np.sin(t)], axis=0) # shape=(2,N)
grid_x = np.linspace(-1.2, 1.2, 13)
grid_y = np.linspace(-1.2, 1.2, 13)
H = []
for gx in grid_x:
H.append(np.stack([np.full_like(grid_y, gx), grid_y], axis=0))
for gy in grid_y:
H.append(np.stack([grid_x, np.full_like(grid_x, gy)], axis=0))
# ======================
# 3) 補間用ユーティリティ
# ======================
def rot2(theta):
c, s = np.cos(theta), np.sin(theta)
return np.array([[c, -s],[s, c]])
# U は 2x2 直交なので角度に落とせる
theta_U = np.arctan2(U[1,0], U[0,0])
def U_interp(alpha):
# alpha∈[0,1] で 0→I, 1→U の回転補間
return rot2(alpha * theta_U)
def Sigma_interp(alpha):
# alpha∈[0,1] で I→Sigma（対角線形補間）
s_lin = 1 + alpha * (s - 1) # 1→s_i
return np.diag(s_lin)
# ======================
# 4) 描画セットアップ（左右2枚のパネル）
# 左：入力 V-平面の z（Phase 1〜2）
# 右：出力平面の y = Uα Σα z（Phase 3 で回転）
# ======================
fig, ax = plt.subplots(1, 2, figsize=(12,6))
for a in ax:
a.set_aspect('equal')
a.axhline(0, linewidth=0.5)
a.axvline(0, linewidth=0.5)
a.set_xlim(-2.5, 2.5)
a.set_ylim(-2.5, 2.5)
ax[0].set_title("入力側：V-平面の座標 z")
ax[1].set_title("出力側：y = U Σ z")
# 左：z平面の図形
circle_in_line, = ax[0].plot([], [], lw=2)
grid_in_lines = [ax[0].plot([], [], lw=0.5, alpha=0.7)[0] for _ in H]
# 右：出力側の図形
circle_out_line, = ax[1].plot([], [], lw=2)
grid_out_lines = [ax[1].plot([], [], lw=0.5, alpha=0.7)[0] for _ in H]
txt0 = ax[0].text(0.02, 0.98, "", transform=ax[0].transAxes, va='top')
txt1 = ax[1].text(0.02, 0.98, "", transform=ax[1].transAxes, va='top')
# ======================
# 5) アニメの設計
# 0→1 の間を 3 フェーズに分割
# P1: V^T（zの座標系合わせ）… zはそのまま（見かけ変化なし、テキストだけ）
# P2: Σ で z を伸縮 → z' = Σα z
# P3: U で y = Uα Σ z へ回転
# ======================
def phase_and_alpha(t_all):
# t_all∈[0,1]
phase = int(3 * t_all) # 0,1,2 のいずれか
alpha = 3 * t_all - phase # 各フェーズ内 0→1
if phase > 2:
phase = 2; alpha = 1.0
return phase, alpha
def update(frame):
t_all = frame / 180 # 0..1
phase, alpha = phase_and_alpha(t_all)
# --- 入力側（左） ---
# P1: 単位円＆格子をそのまま表示（座標合わせ中）
# P2: z を Σα で伸縮（楕円＆ひずんだ格子）
if phase == 0:
circ_in = circle_z
txt0.set_text("Phase 1: 観る向きを V に合わせる（zは変わらない）")
else:
Sigma_a = Sigma_interp(alpha if phase==1 else 1.0)
circ_in = Sigma_a @ circle_z
txt0.set_text("Phase 2: Σ で伸縮（z → Σ z）" if phase==1 else "Phase 3: Σ 完了")
circle_in_line.set_data(circ_in[0], circ_in[1])
for i, seg in enumerate(H):
seg_in = circ_in if seg.shape[1] == circle_z.shape[1] and np.allclose(seg, circle_z) else None
seg_in = Sigma_a @ seg if phase>=1 else seg
grid_in_lines[i].set_data(seg_in[0], seg_in[1])
# --- 出力側（右） ---
# P1: まだ何も適用しないので空・あるいは薄く表示
# P2: y = Σα z を仮表示（Uは未適用）
# P3: y = Uα Σ z
if phase == 0:
circle_out_line.set_data([], [])
for gl in grid_out_lines:
gl.set_data([], [])
txt1.set_text("Phase 1: 出力はまだ表示しない")
elif phase == 1:
Sigma_a = Sigma_interp(alpha)
circ_out = Sigma_a @ circle_z
circle_out_line.set_data(circ_out[0], circ_out[1])
for i, seg in enumerate(H):
seg_out = Sigma_a @ seg
grid_out_lines[i].set_data(seg_out[0], seg_out[1])
txt1.set_text("Phase 2: 出力側に Σ の効果だけ表示")
else:
U_a = U_interp(alpha)
circ_out = U_a @ (Sigma @ circle_z)
circle_out_line.set_data(circ_out[0], circ_out[1])
for i, seg in enumerate(H):
seg_out = U_a @ (Sigma @ seg)
grid_out_lines[i].set_data(seg_out[0], seg_out[1])
txt1.set_text("Phase 3: U で回転（最終的に y = U Σ z）")
return [circle_in_line, *grid_in_lines, circle_out_line, *grid_out_lines, txt0, txt1]
ani = FuncAnimation(fig, update, frames=181, interval=30, blit=False)
HTML(ani.to_jshtml())
Output hidden; open in https://colab.research.google.com to view.
In [60]: