随机奇异值分解（Randomized SVD, rSVD）
假设 X ∈ R n × m X \in \mathbb{R}^{n\times m} X∈Rn×m 为瘦高矩阵， n > m n > m n>m
rSVD 算法流程如下：
# Define randomized SVD function
def rSVD(X,r,q=0,p=0):
# Step 1: Sample column space of X with P matrix
ny = X.shape[1]
P = np.random.randn(ny,r+p)
Z = X @ P
for k in range(q):
Z = X @ (X.T @ Z)
Q, R = np.linalg.qr(Z,mode='reduced')
# Step 2: Compute SVD on projected Y = Q.T @ X
Y = Q.T @ X
UY, S, VT = np.linalg.svd(Y,full_matrices=0)
U = Q @ UY
return U, S, VT
代码详解
案例
原始图片，大小： (3207, 2260)
A = imread(os.path.join('..','DATA','jupiter.jpg'))
X = np.mean(A,axis=2) # Convert RGB -> grayscale
U, S, VT = np.linalg.svd(X,full_matrices=0) # Deterministic SVD, [time spent: 8.07s]
r = 400 # Target rank
q = 1 # Power iterations
p = 5 # Oversampling parameter
rU, rS, rVT = rSVD(X,r,q,p) # [time spent: 0.46s]
SVD 耗时 8.07 s
RSVD 耗时 0.46s，很快啊
## Reconstruction
XSVD = U[:,:(r+1)] @ np.diag(S[:(r+1)]) @ VT[:(r+1),:] # SVD approximation
errSVD = np.linalg.norm(X-XSVD,ord=2) / np.linalg.norm(X,ord=2)
XrSVD = rU[:,:(r+1)] @ np.diag(rS[:(r+1)]) @ rVT[:(r+1),:] # SVD approximation
errSVD = np.linalg.norm(X-XrSVD,ord=2) / np.linalg.norm(X,ord=2)
## Plot
fig, axs = plt.subplots(1,3,figsize=(20,10))
plt.set_cmap('gray')
axs[0].imshow(X)
axs[0].axis('off')
axs[0].set_title('original image')
axs[1].imshow(XSVD)
axs[1].axis('off')
axs[1].set_title('SVD')
axs[2].imshow(XrSVD)
axs[2].axis('off')
axs[2].set_title('rSVD')
plt.show()
参考文献
- http://databookuw.com/
- 《Data-Driven Science and Engineering —— Machine Learning, Dynamical Systems, and Control》