随机奇异值分解（Randomized SVD, rSVD）

假设 $X\in {\mathbb{R}}^{n\times 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》