文章来源于网络收集而来，版权归原创者所有，如有侵权请及时联系！

奇异值分解

发布于 2025-01-01 12:38:39 字数 9908 浏览 0 评论 0 收藏 0

SVD 介绍

“将矩阵分解为我们关心的更简单，更有意义的片段的便捷方式” - 大卫奥斯汀

“我忘了学习的最重要的线性代数概念” - Daniel Lemire

SVD 的应用：

主成分分析
数据压缩
伪逆
协同过滤
主题建模
背景消除
删除损坏的数据

U, s, V = np.linalg.svd(M, full_matrices=False)

这非常慢，因此你可能希望保存结果以便将来使用。

np.save("movie/U.npy", U)
np.save("movie/s.npy", s)
np.save("move/V.npy", V)

将来，你只需加载已保存的内容：

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

U, s, V 是什么样子？

U.shape, s.shape, V.shape

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

练习

检查它们是否是 M 的分解。

# Exercise: 

# True

他们正是。

np.allclose(M, reconstructed_matrix)

# True

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

`s` 的属性

np.diag(s[:6])

你看到 s 的顺序了吗？

s[0:2000:50]

'''
array([ 1341720.,    10528.,     6162.,     4235.,     3174.,     2548.,
           2138.,     1813.,     1558.,     1346.,     1163.,     1001.,
            841.,      666.,        0.,        0.,        0.,        0.,
              0.,        0.,        0.,        0.,        0.,        0.,
              0.,        0.,        0.,        0.,        0.,        0.,
              0.,        0.,        0.,        0.,        0.,        0.,
              0.,        0.,        0.,        0.])
'''

len(s)

# 6000

s[700]

# 3.2309523518534773e-10

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

U 是个巨大的矩阵，所以仅仅查看一小部分：

U[:5,:5]

'''
array([[-0.0083, -0.0009, -0.0007,  0.003 , -0.0002],
       [-0.0083, -0.0013, -0.0005,  0.0034, -0.0001],
       [-0.0084, -0.0012,  0.0002,  0.0045, -0.0003],
       [-0.0085, -0.0011,  0.0001,  0.0044, -0.    ],
       [-0.0086, -0.0013, -0.0002,  0.004 ,  0.0001]])
'''

寻找背景

U.shape, s.shape, V.shape

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

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');

高分辨率版本。

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

制作视频

我受到了 fast.ai 学生萨米尔·穆萨（Samir Moussa）的启发来制作人物视频。

from moviepy.video.io.bindings import mplfig_to_npimage

def make_video(matrix, dims, filename):
    mat_reshaped = np.reshape(matrix, (dims[0], dims[1], -1))

    fig, ax = plt.subplots()
    def make_frame(t):
        ax.clear()
        ax.imshow(mat_reshaped[...,int(t*fps)])
        return mplfig_to_npimage(fig)

    animation = mpe.VideoClip(make_frame, duration=int(10))
    animation.write_videofile('videos/' + filename + '.mp4', fps=fps)

make_video(M - low_rank, dims, "figures2")

'''
[MoviePy] >>>> Building video videos/figures2.mp4
[MoviePy] Writing video videos/figures2.mp4

100%|█████████▉| 600/601 [00:39<00:00, 15.22it/s]

[MoviePy] Done.
[MoviePy] >>>> Video ready: videos/figures2.mp4 
'''

mpe.VideoFileClip("videos/figures2.mp4").subclip(0,10).ipython_display(width=300)

# 100%|█████████▉| 600/601 [00:00<00:00, 858.48it/s]

SVD 分解不同尺寸矩阵的速度

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

一个视频中的两个背景

我们现在使用 BMC 2012 背景建模挑战数据集中的真实视频 008，以及上面使用的 003。

from moviepy.editor import concatenate_videoclips

video2 = mpe.VideoFileClip("movie/Video_008.avi")

concat_video = concatenate_videoclips([video2.subclip(0,20), video.subclip(0,10)])
concat_video.write_videofile("movie/concatenated_video.mp4")

'''
[MoviePy] >>>> Building video movie/concatenated_video.mp4
[MoviePy] Writing video movie/concatenated_video.mp4

100%|█████████▉| 300/301 [00:00<00:00, 481.76it/s]

[MoviePy] Done.
[MoviePy] >>>> Video ready: movie/concatenated_video.mp4 
'''

concat_video = mpe.VideoFileClip("movie/concatenated_video.mp4")

现在回到我们的背景消除问题

concat_video.ipython_display(width=300, maxduration=160)

# 100%|█████████▉| 300/301 [00:00<00:00, 533.88it/s]

scale = 0.5   # 调整比例来更改图像的分辨率
dims = (int(240 * scale), int(320 * scale))

N = create_data_matrix_from_video(concat_video, fps, scale)
# N = np.load("low_res_traffic_matrix.npy")
np.save("med_res_concat_video.npy", N)

N.shape

# (19200, 1800)

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

U_concat, s_concat, V_concat = np.linalg.svd(N, full_matrices=False)

这很慢，因此你可能希望保存结果以便将来使用。

np.save("movie/U_concat.npy", U_concat)
np.save("movie/s_concat.npy", s_concat)
np.save("movie/V_concat.npy", V_concat)

将来，你只需加载已保存的内容：

U_concat = np.load("movie/U_concat.npy")
s_concat = np.load("movie/s_concat.npy")
V_concat = np.load("movie/V_concat.npy")

low_rank = U_concat[:,:10] @ np.diag(s_concat[:10]) @ V_concat[:10,:]

U 的最小主成分：

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

# <matplotlib.image.AxesImage at 0x7f92bcf47da0>

plt.imshow(np.reshape(U_concat[:, 2], dims), cmap='gray')

# <matplotlib.image.AxesImage at 0x7f92bc691a90>

plt.imshow(np.reshape(U_concat[:, 3], dims), cmap='gray')

# <matplotlib.image.AxesImage at 0x7f92bc5aa240>

背景移除：

plt.imshow(np.reshape((N - low_rank)[:, -40], dims), cmap='gray')

# <matplotlib.image.AxesImage at 0x7f92bc540908>

plt.imshow(np.reshape((N - low_rank)[:, 240], dims), cmap='gray')

# <matplotlib.image.AxesImage at 0x7f92bc4d7f28>

数据压缩旁注

假设我们采用 700 个奇异值（记住，总共有 10000 个奇异值）。

s[0:1000:50]

'''
array([ 1341719.6552,    10527.5148,     6162.0638,     4234.9367,
           3174.0389,     2548.4273,     2138.1887,     1812.9873,
           1557.7163,     1345.805 ,     1163.2866,     1000.5186,
            841.4604,      665.7271,        0.    ,        0.    ,
              0.    ,        0.    ,        0.    ,        0.    ])
'''

k = 700
compressed_M = U[:,:k] @ np.diag(s[:k]) @ V[:k,:]

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

# <matplotlib.image.AxesImage at 0x7fefa0076ac8>

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

np.allclose(compressed_M, M)

# True

np.linalg.norm(M - compressed_M)

# 2.864899899979104e-09

U[:,:k].shape, s[:k].shape, V[:k,:].shape

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

节省的空间为对于 700 个奇异值的 U, s, V 中的数据比原始矩阵。

((19200 + 1 + 6000) * 700) / (19200 * 6000)

# 0.1531310763888889

我们只需要存储 15.3% 的数据，并且可以将精度保持在 1e-5 ！很棒！

很漂亮！！！但... SVD 的运行时复杂度为 O(min(m^2 n, m n^2)) 缺点：这真的很慢（同样，我们摒弃了很多计算）。

%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)

分享到QQ

分享到微博