import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as pt
n = 100
A0 = np.random.randn(n, n)
U0, sigma0, VT0 = la.svd(A0)
print(la.norm((U0*sigma0).dot(VT0) - A0))
sigma = np.exp(-np.arange(n))
A = (U0 * sigma).dot(VT0)
pt.semilogy(sigma)
import scipy.linalg.interpolative as sli
Compute a fixed-rank factorization:
(There's also an adaptive, fixed-precision mode.)
k = 20
idx, proj = sli.interp_decomp(A, k)
What does numpy.argsort
do?
idx
sort_idx = np.argsort(idx)
idx[sort_idx]
Reconstruct the matrix:
B = A[:,idx[:k]]
P = np.hstack([np.eye(k), proj])[:,np.argsort(idx)]
Aapprox = np.dot(B, P)
la.norm(A - Aapprox, 2)
What's the structure of $P$?
(ignoring the column permuation)
pt.imshow(np.hstack([np.eye(k), proj]))