7.2 SVD 計算與低秩近似#

import numpy as np

A = np.array([[3., 1.],
              [2., 4.],
              [1., 1.]])              # 3 x 2

U, s, Vt = np.linalg.svd(A)

print(U.shape)    # (3, 3)——方陣！
print(s.shape)    # (2,)——一維向量，不是矩陣
print(Vt.shape)   # (2, 2)

U, s, Vt = np.linalg.svd(A, full_matrices=False)
print(U.shape)    # (3, 2)——只留會用到的欄

# 寫法一：拼出完整 Sigma（full 版）
Sigma = np.zeros(A.shape)
Sigma[:len(s), :len(s)] = np.diag(s)
A_rec = U @ Sigma @ Vt

# 寫法二：精簡版直接乘，最省事
U, s, Vt = np.linalg.svd(A, full_matrices=False)
A_rec = (U * s) @ Vt        # 廣播：把 s 乘到 U 的每一欄
print(np.allclose(A_rec, A))   # True

s = np.linalg.svd(A, compute_uv=False)

import numpy as np
import matplotlib.pyplot as plt

# 模擬「真實資料」：低秩結構 + 少量雜訊
rng = np.random.rand(200, 50)           # 這裡先示範用隨機矩陣
s = np.linalg.svd(rng, compute_uv=False)

plt.plot(s, "o-")
plt.xlabel("索引 $i$"); plt.ylabel(r"奇異值 $\sigma_i$")
plt.title("奇異值衰減曲線")
plt.show()

import numpy as np

rng = np.random.default_rng(0)
A = rng.random((30, 20))

U, s, Vt = np.linalg.svd(A, full_matrices=False)

k = 5
Ak = (U[:, :k] * s[:k]) @ Vt[:k, :]

# 譜範數誤差應等於 s[k]（第 k+1 大的奇異值）
print(np.linalg.norm(A - Ak, 2), s[k])
# Frobenius 誤差應等於剩餘奇異值的平方和開根號
print(np.linalg.norm(A - Ak, "fro"), np.sqrt((s[k:]**2).sum()))

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(1)
# 有明顯低秩結構的矩陣：兩層結構 + 雜訊
base = np.outer(rng.random(60), rng.random(40)) \
     + np.outer(rng.random(60), rng.random(40))
A = base + 0.05 * rng.random((60, 40))

U, s, Vt = np.linalg.svd(A, full_matrices=False)
energy = s**2 / (s**2).sum()                       # 每個方向的能量占比
cum_energy = np.cumsum(energy)                     # 前 k 個方向累積能量

k_all = np.arange(1, len(s) + 1)
err2 = np.sqrt((s**2).sum() - np.cumsum(s**2))     # Frobenius 誤差

fig, axes = plt.subplots(1, 2, figsize=(11, 4))
axes[0].plot(k_all, cum_energy, "o-")
axes[0].set_xlabel("$k$"); axes[0].set_ylabel("累積能量比例")
axes[0].set_title("前 $k$ 個奇異值保留多少能量")
axes[1].plot(k_all, err2, "o-", color="tomato")
axes[1].set_xlabel("$k$"); axes[1].set_ylabel(r"$\|A - A_k\|_F$")
axes[1].set_title("重建誤差隨 $k$ 衰減")
plt.tight_layout(); plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 延續上一格的 A, U, s, Vt
ks = [1, 2, 3, 5, 10, 40]

fig, axes = plt.subplots(2, 3, figsize=(12, 6))
for ax, k in zip(axes.ravel(), ks):
    Ak = (U[:, :k] * s[:k]) @ Vt[:k, :]
    im = ax.imshow(Ak, cmap="viridis")
    rel = np.linalg.norm(A - Ak, "fro") / np.linalg.norm(A, "fro")
    ax.set_title(f"$k={k}$，相對誤差 $={rel:.3f}$")
    ax.axis("off")
fig.suptitle("截斷 SVD：不同 $k$ 的低秩近似")
plt.tight_layout(); plt.show()