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])
%%time
u, s, vh = np.linalg.svd(gray,full_matrices=False)
newGray = (u*s)@vh
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()
#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
print("Approx error is:",np.linalg.norm(gray -newGray))
plt.imshow(newGray, cmap = plt.get_cmap('gray'))
plt.title("k=%s"%(k))
plt.show()
%%time
k=10
q=1
U,S,Vt = full_randomSVD(gray,k, q)
newGray = (U*S)@Vt
plt.imshow(newGray, cmap = plt.get_cmap('gray'))
plt.title("k=%s q=%s"%(k,q))
plt.show()
%%time
k=10
q=10
U,S,Vt = full_randomSVD(gray,k, q)
newGray = (U*S)@Vt
plt.imshow(newGray, cmap = plt.get_cmap('gray'))
plt.title("k=%s q=%s"%(k,q))
plt.show()
Now to see that the unstable version is in fact unstable
k=10
q=10
U,S,Vt = full_unstable_randomSVD(gray,k, q)
newGray = (U*S)@Vt
plt.imshow(newGray, cmap = plt.get_cmap('gray'))
plt.title("k=%s q=%s"%(k,q))
plt.show()
#===============================================
%load_ext watermark
%watermark