import matplotlib.pyplot as plt
import matplotlib.image as mpimg #if you want to convert your own image
import numpy as np

np.random.seed(1234) # for reproducing results

gray = np.loadtxt("titan_matrix.txt")
plt.title("Original Gray Image")
plt.imshow(gray,cmap = plt.get_cmap('gray'))
plt.show()

#an idea of what the data looks like
print(gray[:10,:10]) 

[[64.222 64.222 65.222 65.222 66.222 67.222 67.934 67.934 64.934 69.646]
 [64.222 65.222 65.222 66.222 67.222 68.222 67.934 68.934 67.934 70.646]
 [65.222 66.222 66.222 67.222 68.222 69.222 68.934 69.934 70.934 71.646]
 [66.222 66.222 67.222 68.222 69.222 69.222 69.934 69.934 71.934 71.646]
 [66.222 66.222 67.222 68.222 69.222 69.222 69.934 69.934 70.934 69.646]
 [65.222 66.222 66.222 67.222 68.222 69.222 68.934 69.934 67.934 68.646]
 [64.222 65.222 65.222 66.222 67.222 68.222 67.934 68.934 65.934 68.646]
 [64.222 64.222 65.222 65.222 66.222 67.222 67.934 67.934 64.934 69.646]
 [66.222 66.222 66.222 67.222 67.222 68.222 67.934 67.934 65.934 65.646]
 [66.222 66.222 66.222 67.222 67.222 68.222 67.934 67.934 67.934 67.646]]

Correct SVD

%%time
u, s, vh = np.linalg.svd(gray,full_matrices=False)
newGray = (u*s)@vh

CPU times: user 739 ms, sys: 52.5 ms, total: 791 ms
Wall time: 483 ms

plt.plot(s)
plt.yscale('log')
plt.title("Magnitude of Singular Values")
plt.show()

k=10
u, s, vh = np.linalg.svd(gray,full_matrices=False)
s[-(len(s)-k):]=0

newGray = (u*s)@vh
plt.imshow(newGray, cmap = plt.get_cmap('gray'))
plt.title("True k=%s"%(k))
plt.show()

Random SVD implementations

#straight forward implementation

def easy_randomSVD(M,k):
    m,n = M.shape
    
    if m < n:
        M=M.T
    Ω   = np.random.normal(0,1,(M.shape[1],k))
    Q,R = np.linalg.qr(M @ Ω)
    B   = Q.T @ M
    u,S,Vt = np.linalg.svd(B,full_matrices=False)
    
    U = Q @ u
    if m<n:
        return Vt.T,S.T,U.T
    else:
        return U,S,Vt

# includes projection, but does it in a bad way
def full_unstable_randomSVD(M,k,q):
    m,n = M.shape
    
    if m < n:
        M=M.T
    Ω = np.random.normal(0,1,(M.shape[1],k))
    Y = M @ Ω
    MM = (M@M.T) #cache
    for j in range(q): 
        Y = MM @ Y
    Q,R = np.linalg.qr(Y)
    B   = Q.T @ M
    u,S,Vt = np.linalg.svd(B,full_matrices=False)
    
    U = Q @ u
    if m<n:
        return Vt.T,S.T,U.T
    else:
        return U,S,Vt   

#randomSVD with a more stable solution (but slow QR)
def full_randomSVD(M,k,q):
    m,n = M.shape
    
    if m < n:
        M=M.T
    Ω   = np.random.normal(0,1,(M.shape[1],k))
    Q,R = np.linalg.qr(M @ Ω)
    #Randomized Subspace Iteration
    for j in range(q):
        Q,R = np.linalg.qr(M.T @ Q )
        Q,R = np.linalg.qr(M @ Q  )
    B   = Q.T @ M
    u,S,Vt = np.linalg.svd(B,full_matrices=False)
    
    U = Q @ u
    if m<n:
        return Vt.T,S.T,U.T
    else:
        return U,S,Vt

%%time
k=10

U,S,Vt = easy_randomSVD(gray,k)
newGray = (U*S)@Vt

CPU times: user 13.7 ms, sys: 2.56 ms, total: 16.3 ms
Wall time: 9.92 ms

print("Approx error is:",np.linalg.norm(gray -newGray))
plt.imshow(newGray, cmap = plt.get_cmap('gray'))
plt.title("k=%s"%(k))
plt.show()

Approx error is: 19015.6076409738