# coding: utf-8 # # Randomized SVD # In[2]: import numpy as np import numpy.linalg as la import matplotlib.pyplot as pt # ## Obtaining a low-rank matrix # Let's begin by making a low-rank matrix: # In[3]: n = 100 A0 = np.random.randn(n, n) U0, sigma0, VT0 = la.svd(A0) print(la.norm((U0*sigma0).dot(VT0) - A0)) pt.plot(sigma0) # A randomly drawn Gaussian matrix: Emphatically not low-rank. Let's swap out the singular values with cleanly exponentially decaying ones: # In[4]: sigma = np.exp(-np.arange(n)) pt.semilogy(sigma) A = (U0 * sigma).dot(VT0) # ## Find the approximate range # Let's fix parameters first. What accuracy should we obtain for the values of \$k\$ and \$p\$ below? (where `p` represents the 'extra' dimensions) # In[5]: k = 10 p = 5 # Draw a random Gaussian matrix `Omega` and obtain orthogonal columns in a matrix `Q` spanning the range: # In[6]: Omega = np.random.randn(n, k+p) Y = A.dot(Omega) Q, _ = la.qr(Y) # As an alternative to the above, use a few iterations of the power method on `Omega`: # In[7]: Omega = np.random.randn(n, k+p) Y = A.dot(Omega) Y = A.dot(A.T.dot(Y)) Q, _ = la.qr(Y) # * Observations about associativity? # ## Convert to factorization form # Reconstruct \$C\$ in the factorization-form LRA \$A\approx BC\$: # # (Recall \$A\approx QQ^TA\$, \$B=Q\$) # In[8]: C = Q.T.dot(A) print(C.shape) # Sanity-check that form: # In[9]: Aapprox1 = Q.dot(C) la.norm(A - Aapprox1, 2) # ## Compute the approximate SVD # Compute the SVD of \$C=U_C \Sigma_C V_C^T\$: # # (Make sure to pass `full_matrices=False` to the SVD.) # In[10]: UC, sigmaC, VTC = la.svd(C, full_matrices=False) # Reconstruct the SVD of \$A\$: \$A\approx QU_C \Sigma_C V_C^T\$ # In[11]: UAapprox = Q.dot(UC) sigmaAapprox = sigmaC VTAapprox = VTC # Compare the 2-norm of the reconstructed \$A\$ with the original \$A\$: # In[12]: Aapprox = (UAapprox*sigmaAapprox).dot(VTAapprox) la.norm(A - Aapprox, 2) # ## Estimate the error # # Compute an a-posteriori estimate of approximation error in the spectral norm: # In[18]: omega = np.random.randn(n) Aomega = A.dot(omega) err = Aomega - Q.dot(Q.T.dot(Aomega)) la.norm(err, 2) / la.norm(omega, 2) # * Is this the right direction for the error estimator? # * Is the estimator supposed to be conservative? Can it be?