Question

我们今天的目标：

加载和格式化数据

让我们使用 BMC 2012 背景模型挑战数据集中的真实视频 003。

导入所需的库：

import imageio
imageio.plugins.ffmpeg.download()

import moviepy.editor as mpe
import numpy as np
import scipy

%matplotlib inline
import matplotlib.pyplot as plt

scale = 0.50   # 调整比例来更改图像的分辨率
dims = (int(240 * scale), int(320 * scale))
fps = 60      # 每秒的帧

M = np.load("movie/med_res_surveillance_matrix_60fps.npy")

print(dims, M.shape)

# (120, 160) (19200, 6000)

plt.imshow(np.reshape(M[:,140], dims), cmap='gray');

plt.figure(figsize=(12, 12))
plt.imshow(M, cmap='gray')

# <matplotlib.image.AxesImage at 0x7f601f315fd0>

将来，你可以加载已保存的内容：

U = np.load("U.npy")
s = np.load("s.npy")
V = np.load("V.npy")

U, S, V 是什么样呢？

U.shape, s.shape, V.shape

# ((19200, 6000), (6000,), (6000, 6000))

检查它们是 M 的分解。

reconstructed_matrix = U @ np.diag(s) @ V

np.allclose(M, reconstructed_matrix)

# True

是的。

移除背景

low_rank = np.expand_dims(U[:,0], 1) * s[0] * np.expand_dims(V[0,:], 0)

plt.figure(figsize=(12, 12))
plt.imshow(low_rank, cmap='gray')

# <matplotlib.image.AxesImage at 0x7f1cc3e2c9e8>

plt.imshow(np.reshape(low_rank[:,0], dims), cmap='gray');

我们如何获取里面的人？

plt.imshow(np.reshape(M[:,0] - low_rank[:,0], dims), cmap='gray');

SVD 对不同大小的矩阵的速度

s 是对角矩阵的对角线。

np.set_printoptions(suppress=True, precision=4)

import timeit
import pandas as pd

m_array = np.array([100, int(1e3), int(1e4)])
n_array = np.array([100, int(1e3), int(1e4)])

index = pd.MultiIndex.from_product([m_array, n_array], names=['# rows', '# cols'])

pd.options.display.float_format = '{:,.3f}'.format
df = pd.DataFrame(index=m_array, columns=n_array)

# %%prun
for m in m_array:
    for n in n_array:      
        A = np.random.uniform(-40,40,[m,n])  
        t = timeit.timeit('np.linalg.svd(A, full_matrices=False)', number=3, globals=globals())
        df.set_value(m, n, t)

df/3

	100	1000	10000
100	0.006	0.009	0.043
1000	0.004	0.259	0.992
10000	0.019	0.984	218.726

很好！！！但是...

缺点：这真的很慢（同样，我们摒弃了很多计算）。

%time u, s, v = np.linalg.svd(M, full_matrices=False)

'''
CPU times: user 5min 38s, sys: 1.53 s, total: 5min 40s
Wall time: 57.1 s
'''

M.shape

# (19200, 6000)

随机化 SVD 的最简单版本

想法：让我们使用更小的矩阵！

我们还没有找到更好的通用 SVD 方法，我们只会使用我们在较小矩阵上使用的方法，该矩阵与原始矩阵的范围大致相同。

def simple_randomized_svd(M, k=10):
    m, n = M.shape
    transpose = False
    if m < n:
        transpose = True
        M = M.T

    rand_matrix = np.random.normal(size=(M.shape[1], k))  # short side by k
    Q, _ = np.linalg.qr(M @ rand_matrix, mode='reduced')  # long side by k
    smaller_matrix = Q.T @ M                              # k by short side
    U_hat, s, V = np.linalg.svd(smaller_matrix, full_matrices=False)
    U = Q @ U_hat

    if transpose:
        return V.T, s.T, U.T
    else:
        return U, s, V

%time u, s, v = simple_randomized_svd(M, 10)

'''
CPU times: user 3.06 s, sys: 268 ms, total: 3.33 s
Wall time: 789 ms
'''

U_rand, s_rand, V_rand = simple_randomized_svd(M, 10)

low_rank = np.expand_dims(U_rand[:,0], 1) * s_rand[0] * np.expand_dims(V_rand[0,:], 0)

plt.imshow(np.reshape(low_rank[:,0], dims), cmap='gray');

我们如何获取里面的人？

这个方法在做什么

rand_matrix = np.random.normal(size=(M.shape[1], 10))

rand_matrix.shape

# (6000, 10)

plt.imshow(np.reshape(rand_matrix[:4900,0], (70,70)), cmap='gray');

temp = M @ rand_matrix; temp.shape

# (19200, 10)

plt.imshow(np.reshape(temp[:,0], dims), cmap='gray');

plt.imshow(np.reshape(temp[:,1], dims), cmap='gray');

Q, _ = np.linalg.qr(M @ rand_matrix, mode='reduced'); Q.shape

# (19200, 10)

np.dot(Q[:,0], Q[:,1])

# -3.8163916471489756e-17

plt.imshow(np.reshape(Q[:,0], dims), cmap='gray');

plt.imshow(np.reshape(Q[:,1], dims), cmap='gray');

smaller_matrix = Q.T @ M; smaller_matrix.shape

# (10, 6000)

U_hat, s, V = np.linalg.svd(smaller_matrix, full_matrices=False)

U = Q @ U_hat

plt.imshow(np.reshape(U[:,0], dims), cmap='gray');

reconstructed_small_M = U @ np.diag(s) @ V

以及人。

plt.imshow(np.reshape(M[:,0] - reconstructed_small_M[:,0], dims), cmap='gray');

时间比较

from sklearn import decomposition
import fbpca

完整的 SVD：

%time u, s, v = np.linalg.svd(M, full_matrices=False)

'''
CPU times: user 5min 38s, sys: 1.53 s, total: 5min 40s
Wall time: 57.1 s
'''

我们的（过度简化）的 randomized_svd ：

%time u, s, v = simple_randomized_svd(M, 10)

'''
CPU times: user 2.37 s, sys: 160 ms, total: 2.53 s
Wall time: 641 ms
'''

Sklearn：

%time u, s, v = decomposition.randomized_svd(M, 10)

'''
CPU times: user 19.2 s, sys: 1.44 s, total: 20.7 s
Wall time: 3.67 s
'''

来自 Facebook fbpca 库的随机 SVD：

%time u, s, v = fbpca.pca(M, 10)

'''
CPU times: user 7.28 s, sys: 424 ms, total: 7.7 s
Wall time: 1.37 s
'''

我会选择 fbpca，因为它比 sklearn 更快，比我们简单的实现更健壮，更准确。

以下是 Facebook Research 的一些结果：

k 变化下的时间和准确度

import timeit
import pandas as pd

U_rand, s_rand, V_rand = fbpca.pca(M, 700, raw=True)
reconstructed = U_rand @ np.diag(s_rand) @ V_rand

np.linalg.norm(M - reconstructed)

# 1.1065914828881536e-07

plt.imshow(np.reshape(reconstructed[:,140], dims), cmap='gray');

pd.options.display.float_format = '{:,.2f}'.format
k_values = np.arange(100,1000,100)
df_rand = pd.DataFrame(index=["time", "error"], columns=k_values)

# df_rand = pd.read_pickle("svd_df")

for k in k_values:
    U_rand, s_rand, V_rand = fbpca.pca(M, k, raw=True)
    reconstructed = U_rand @ np.diag(s_rand) @ V_rand
    df_rand.set_value("error", k, np.linalg.norm(M - reconstructed))
    t = timeit.timeit('fbpca.pca(M, k)', number=3, globals=globals())
    df_rand.set_value("time", k, t/3)

df_rand.to_pickle("df_rand")

df_rand

	100	200	300	400	500	600	700	800	900	1000
time	2.07	2.57	3.45	6.44	7.99	9.02	10.24	11.70	13.30	10.87
error	58,997.27	37,539.54	26,569.89	18,769.37	12,559.34	6,936.17	0.00	0.00	0.00	0.00

df = pd.DataFrame(index=["error"], columns=k_values)

for k in k_values:
    reconstructed = U[:,:k] @ np.diag(s[:k]) @ V[:k,:]
    df.set_value("error", k, np.linalg.norm(M - reconstructed))

df.to_pickle("df")

fig, ax1 = plt.subplots()
ax1.plot(df.columns, df_rand.loc["time"].values, 'b-', label="randomized SVD time")
ax1.plot(df.columns, np.tile([57], 9), 'g-', label="SVD time")
ax1.set_xlabel('k: # of singular values')
# Make the y-axis label, ticks and tick labels match the line color.
ax1.set_ylabel('time', color='b')
ax1.tick_params('y', colors='b')
ax1.legend(loc = 0)

ax2 = ax1.twinx()
ax2.plot(df.columns, df_rand.loc["error"].values, 'r--', label="randomized SVD error")
ax2.plot(df.columns, df.loc["error"].values, 'm--', label="SVD error")
ax2.set_ylabel('error', color='r')
ax2.tick_params('y', colors='r')
ax2.legend(loc=1)

#fig.tight_layout()
plt.show()

数学细节

随机 SVD 背后的处理

下面是一个计算截断 SVD 的过程，在“ 带有随机性的搜索结构：用于构造近似矩阵分解的概率算法 ”中描述，并在 此博客文章 中总结：

1.计算 A 的近似范围。也就是说，我们希望 Q 具有 r 个正交列，使得：

2.构造 ，它比较小 r×n 。

3.通过标准方法计算 B 的 SVD（因为 B 小于 A 所以更快），。

1. 由于：

如果我们设置 U = QS ，那么我们有一个低秩的近似值 。

那么我们如何找到 Q （步骤 1）？

为了估计 A 的范围，我们可以只取一堆随机向量 ，来求解 形成的子空间。 我们可以用 作为列来形成矩阵 W 。 现在，我们采用 AW = QR 的 QR 分解，然后 Q 的列形成 AW 的标准正交基，这是 A 的范围。

由于乘积矩阵 AW 的行比列更多，因此列大致正交。 这是一个简单的概率 - 有很多行，列很少，列不可能是线性相关的。

为什么

我们试图找到矩阵 Q ，使得 。 我们对 M 的范围很感兴趣，我们称之为 MX 。 Q 有正交列，因此 （但 不是 I ，因为 Q 是矩形的）。

于是...

如果 X 是单位，我们就做成了（但是 X 会太大，我们不会得到我们想要的加速）。 在我们的问题中， X 只是一个小的随机矩阵。 Johnson-Lindenstrauss 引理为其工作原理提供了一些理由。

QR 分解

我们稍后将深入了解 QR 分解。 现在，你只需要知道 A = QR ，其中 Q 由正交列组成， R 是上三角形。 Trefethen 说 QR 分解是数值线性代数中最重要的思想！ 我们一定会将回顾它。

我们该如何选择 r ？ 假设我们的矩阵有 100 列，我们想要 U 和 V 中的 5 列。为了安全起见，我们应该将矩阵投影到正交基上，其中的行数和列数多于 5（让我们使用 15）。 最后，我们取 U 和 V 的前 5 列。

因此，即使我们的投影只是近似的，通过使它比我们需要的更大，我们可以弥补精度的损失（因为我们随后只采用了一个子集）。

这与随机均有何不同

test = M @ np.random.normal(size=(M.shape[1], 2)); test.shape

# (4800, 2)

随机均值：

plt.imshow(np.reshape(test[:,0], dims), cmap='gray');

均值图像：

plt.imshow(np.reshape(M.mean(axis=1), dims), cmap='gray')

# <matplotlib.image.AxesImage at 0x7f83f4093fd0>

ut, st, vt = np.linalg.svd(test, full_matrices=False)

plt.imshow(np.reshape(smaller_matrix[0,:], dims), cmap='gray');

plt.imshow(np.reshape(smaller_matrix[1,:], dims), cmap='gray');

plt.imshow(np.reshape(M[:,140], dims), cmap='gray');