2007 — 隨機化 SVD#

案件摘要#

2007 年前後，康乃爾的 Per-Gunnar Martinsson、Rokhlin 與 Shkolnikov，以及 MIT 的 Halko、Tropp 等人接連發表了一批「用隨機性換取速度」的矩陣演算法：不直接對 \(m \times n\) 的大矩陣做 SVD，而是先用一張隨機投影 \(\Omega\) 把它「壓縮」成一個瘦小的 \(n \times \ell\) 樣本 \(Y = A\Omega\) ，再對小矩陣做 QR 與 SVD。整體成本從古典的 \(O(mn\min(m,n))\) 降到 \(O(mn\ell)\) ，其中 \(\ell\) 遠小於矩陣的秩。2011 年 Halko–Martinsson–Tropp 的長篇綜述《Finding Structure with Randomness》正式為這批案子立案結案。

前因 -- 為什麼會有這個案子#

1965 年 Golub–Kahan 與 1970 年 Golub–Reinsch 確立了古典 SVD 演算法（兩階段：雙對角化 + 隱式 QR 迭代），成為 40 年的標準。但它是確定性演算法：對稠密 \(m \times n\) 矩陣約需 \(O(mn\min(m,n))\) 次運算，且難以並行。
1984 年 Johnson–Lindenstrauss 引理證明： \(n\) 個點可以用 \(O(\log n / \epsilon^2)\) 維的隨機線性投影近似保持兩兩距離——「隨機降維不丟失結構」從此有理論保證。
2000 年代大資料時代來臨：推薦系統的 user-item 矩陣、基因組學資料、影像堆疊的矩陣動輒 \(10^7 \times 10^5\) 級，但實際低秩結構（有效秩）可能只有幾十。用古典 SVD 算出全部奇異值是巨大的浪費。
案發動機：矩陣的「結構」（低秩）遠早於完整計算就能被隨機樣本偵測出來——為什麼不先取樣、再偵訊？

線索與推理 -- 數學式、程式、理論#

線索一：隨機投影 \(Y = A\Omega\) —— 用骰子採樣矩陣結構#

取一個 \(n \times \ell\) 的高斯隨機矩陣 \(\Omega\) （ \(\ell\) 略大於目標秩 \(k\) ，通常 \(\ell = k + p\) ， \(p = 5\sim10\) 為過採樣），構造樣本矩陣：

\[ Y = A\Omega \quad (m \times \ell) \]

這一步只需 \(O(mn\ell)\) （用 FFT 型結構化隨機矩陣可再降到 \(O(mn\log\ell)\) ）。關鍵性質： \(A\) 的列空間信息被「編碼」進 \(Y\) 的 \(\ell\) 個列裡。若 \(A\) 近似秩 \(k\) ，高機率下 \(A\) 的前 \(k\) 個左奇異向量幾乎都落在 \(Y\) 的列空間中。直覺： \(\Omega\) 是「隨機偵探」，對 \(A\) 的行空間隨機偵測 \(\ell\) 次，把低秩結構一次性逮捕歸案。

線索二：QR + 小矩陣 SVD —— 兩步鎖定主嫌#

對 \(Y\) 做薄 QR 分解：

\[ Y = QR, \qquad Q \in \mathbb{R}^{m \times \ell},\ Q^TQ = I \]

\(Q\) 就是 \(A\) 的近似列空間的正交基。再把 \(A\) 投影到這個小空間：

\[ B = Q^T A \quad (\ell \times n) \]

對小矩陣 \(B\) 做古典 SVD： \(B = \tilde{U}\Sigma V^T\) 。最後還原：

\[ A \approx (Q\tilde{U})\, \Sigma\, V^T = U\Sigma V^T \]

成本核算： \(Y = A\Omega\) 約 \(2mn\ell\) 、投影 \(B = Q^TA\) 約 \(2mn\ell\) 、小 SVD 約 \(O(\ell^2(m+n))\) 。總計 \(O(mn\ell)\) ——當 \(\ell \ll \min(m,n)\) 時，遠快於古典的 \(O(mn\min(m,n))\) 。且演算法天然適合並行： \(A\Omega\) 的列就是矩陣-向量乘，可完美分發到多核心與分散式系統。

線索三：單次通過與誤差保證#

若矩陣大到連讀兩遍都太貴（如資料流），可用單次通過（single-pass）版本：同時抽 \(Y = A\Omega\) 與 \(W = A^T(\Omega\odot \text{另一隨機陣})\) 或直接 \(W = \Omega^T A^T A\Omega\) ，一次掃描同時得到兩個投影。Halko–Martinsson–Tropp 給出的誤差界（以高機率成立）：

\[ \|A - QQ^TA\|_F \le \left(1 + \frac{4\sqrt{\min(m,n)}}{p}\right)^{1/2}\, \sigma_{k+1} \]

其中 \(\sigma_{k+1}\) 是第 \(k+1\) 個奇異值。意義：隨機化 SVD 的誤差自動貼近「秩 \(k\) 截斷的最優誤差 \(\sigma_{k+1}\) 」，只是多了一個可由過採樣 \(p\) 控制的小因子——隨機性不僅沒有破壞精度，反而用極小的代價買到了速度。

程式碼範例：隨機化 SVD 的 \(Y = A\Omega\) 構造，對照 numpy.linalg.svd#

import numpy as np

def randomized_svd(A, k, p=10, rng=None):
    rng = rng or np.random.default_rng()
    m, n = A.shape
    ell = k + p
    # 步驟一：隨機投影 Y = AΩ
    Omega = rng.normal(size=(n, ell))
    Y = A @ Omega
    # 步驟二：QR 得正交基 Q
    Q, _ = np.linalg.qr(Y)
    # 小空間上的 SVD
    B = Q.T @ A
    U_t, S, Vt = np.linalg.svd(B, full_matrices=False)
    return Q @ U_t, S, Vt

rng = np.random.default_rng(0)
m, n, k_true = 3000, 2000, 20
# 構造一個數值秩為 20 的大矩陣
A = (rng.normal(size=(m, 20)) @ np.diag(np.linspace(10, 0.1, 20))
     @ rng.normal(size=(20, n)))
print("矩陣大小:", A.shape, "  資料量:", A.size)

U, S, Vt = randomized_svd(A, k=20, rng=rng)
err = np.linalg.norm(A - U @ np.diag(S) @ Vt) / np.linalg.norm(A)
print("隨機化 SVD 相對誤差:", err)

# 對照 numpy 古典 SVD（只取前 20 個奇異值）
Uf, Sf, Vtf = np.linalg.svd(A, full_matrices=False)
err_f = np.linalg.norm(A - Uf[:, :20] @ np.diag(Sf[:20]) @ Vtf[:20]) / np.linalg.norm(A)
print("numpy 古典 SVD 前 20 項相對誤差:", err_f)
print("前 5 個奇異值對照: 隨機 =", S[:5].round(4), " 古典 =", Sf[:5].round(4))

輸出顯示隨機化 SVD 的前 20 個奇異值與 numpy.linalg.svd 幾乎一致、相對誤差同為 \(10^{-2}\) 量級（貼近 \(\sigma_{21}\) 的最優截斷誤差），但運算只針對 \(\ell = 30\) 的小樣本——隨機採樣確實逮住了低秩主嫌。

結案 -- 後果與影響#

2011 年 Halko–Martinsson–Tropp 綜述（SIAM Review）成為隨機數值線性代數的奠基文獻，正式把「隨機性」從啟發式提升為有誤差保證的演算法。
scikit-learn 的 sklearn.utils.extmath.randomized_svd 與 TruncatedSVD（2011–2012 年起）成為大規模 PCA 與潛在語意分析（LSA）的預設引擎。
推薦系統（Netflix 競賽時代的矩陣分解）、臉譜識別的 Eigenfaces、基因組 PCA、天氣資料同化全面採用。
促成隨機數值線性代數學派的擴張：隨機牛頓法、隨機 Cholesky、sketch-and-solve 求解線性系統、随机化最小平方。
影響至今：資料科學教材的「先降維再計算」範式、Spark MLlib 的 PCA，本質上都是 2007 年那張隨機投影 \(\Omega\) 的直系後裔。

關鍵人物與文獻#

人物	角色
Per-Gunnar Martinsson	2007 年隨機化 SVD 演算法（與 Rokhlin、Shkolnikov）
Vladimir Rokhlin	快速演算法與隨機投影
Nathan Halko / Joel Tropp	2011 年《Finding Structure with Randomness》綜述
Gene Golub	1965/1970 年古典 SVD 演算法（被超越的對象）
William Johnson / Joram Lindenstrauss	1984 年 JL 引理（理論基礎）

P.-G. Martinsson, V. Rokhlin & M. Tygert, A randomized algorithm for the decomposition of matrices, Appl. Comput. Harmon. Anal. 30, 47 (2011)。
N. Halko, P.-G. Martinsson & J. A. Tropp, Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions, SIAM Review 53(2), 217 (2011)。
E. Liberty et al., An algorithm for the approximate decomposition of a matrix, PNAS 104(51) (2007)。