Halko-Martinsson-Tropp Randomized Subspace Iteration and Low-Rank Matrix Decomposition

Stage A constructs an orthonormal range basis Q via Gaussian projection Y = A*Omega (Algorithm 4.1), adaptive a posteriori error estimation terminating when r consecutive projected Gaussian vectors satisfy ||(I - Q^{(j)}(Q^{(j)})^*)A*omega^{(i)}|| <= eps / (10*sqrt(2/pi)) with failure probability at most min{m, n}*10^{-r} (Algorithm 4.2, Lemma 4.1), QR-stabilized subspace iteration Y_0 = A*Omega = Q_0*R_0, Y_tilde_j = A^* Q_{j-1} = Q_tilde_j*R_tilde_j, Y_j = A*Q_tilde_j = Q_j*R_j (Algorithm 4.4), or Subsampled Random Fourier Transform Omega = sqrt(n/l)*D*F*R in O(m*n*log(l) + l^2*n) flops (Algorithm 4.5); Stage B then extracts Direct SVD B = Q^* A = U_tilde*Sigma*V^*, U = Q*U_tilde preserving ||A - U*Sigma*V^*|| <= eps (Algorithm 5.1) or PSD Nystrom Cholesky-SVD A approx (A*Q)*(Q^* A*Q)^{-1}*(A*Q)^* whose spectral error never exceeds ||A - Q*Q^* A|| (Algorithm 5.5).

Using the Loewner inverse perturbation I - (I + F^* F)^{-1} <= F^* F (Proposition 8.2) on the flattened range matrix Z = [I; F] with F = Sigma_2*Omega_2*Omega_1^dagger*Sigma_1^{-1}, Theorem 9.1 proves |||(I - P_Y)A|||^2 <= |||Sigma_2|||^2 + |||Sigma_2*Omega_2*Omega_1^dagger|||^2 in both spectral and Frobenius norms; combined with Proposition 8.6 (||P*M|| <= ||P*(M*M^*)^q*M||^{1/(2q+1)}), Proposition 10.2 (E||Omega_1^dagger||_F^2 = k/(p - 1), E||Omega_1^dagger|| <= e*sqrt(k + p)/p), and Gordon's inequality, this yields E||(I - P_Y)A||_F <= sqrt(1 + k/(p - 1))*||Sigma_2||_F (Theorem 10.5), E||(I - P_Y)A|| <= (1 + sqrt(k/(p - 1)))*sigma_{k+1} + (e*sqrt(k + p)/p)*||Sigma_2||_F (Theorem 10.6), and exponential power-iteration decay in (2q + 1) (Theorem 9.2, Corollary 10.10).

Applied to discretized single-layer Laplace contour operators and 140,000 x 140,000 lattice transfer matrices, nonlocal 9025 x 9025 patch-similarity graph Laplacians for image processing, and out-of-core 98,304 x 7,254 FERET eigenface matrices (5.4 GB uncompressed), the framework achieves machine-precision operator compression and 3-5 pass spectral convergence with multi-fold speedups over deterministic rank-revealing QR and LAPACK dgesdd.

Executing unnormalized power iteration Y = (A*A^*)^q A*Omega (Algorithm 4.3) in floating-point arithmetic with machine precision mu extinguishes all singular modes smaller than mu^{1/(2q+1)}*||A|| unless intermediate QR orthonormalizations are performed between matrix-vector products (Remark 4.3, Algorithm 4.4); additionally, single-pass streaming recovery B_approx*(Q^* Omega) approx Q^* Y (Algorithm 5.6) amplifies error by 1 / tau_min(Q^* Omega) when Q^* Omega is ill-conditioned (Remark 5.4), and SRFT matrices require O((k + log n)*log k) samples in worst-case coupon-collector decimation regimes with O(k^{-1}) failure probability (Remark 11.2, Theorem 11.2).