机器学习算法那些事 | 这是我见过最通俗易懂的SVD（奇异值分解）算法介绍

本文来源公众号“机器学习算法那些事”，仅用于学术分享，侵权删，干货满满。

原文链接：这是我见过最通俗易懂的SVD（奇异值分解）算法介绍

线性代数是机器学习领域的基础，其中一个最重要的概念是奇异值分解（SVD），本文尽可能简洁的介绍SVD（奇异值分解）算法的基础理解，以及它在现实世界中的应用。

SVD是最广泛使用的无监督学习算法之一，它在许多推荐系统和降维系统中居于核心位置，这些系统是全球公司如谷歌、Netflix、Facebook、YouTube等的核心技术。

简单来说，SVD是将一个任意矩阵分解为三个矩阵。所以如果我们有一个矩阵A，那么它的SVD可以表示为：

image-20240808230323995

与特征值分解相比，奇异值分解可以应用于非方阵。在SVD中，U和 V 对于任何矩阵都是可逆的，并且它们是正交归一的，这是我们所喜爱的特性。虽然这里不进行证明，但我们可以告诉你，奇异值比特征值在数值上更稳定。

为了更好地理解，我们通过一个例子演示SVD

奇异值是正特征值的平方根，即5和3。因此非方阵A的SVD分解为：

SVD分解证明

image-20240810203329639

SVD的另一种表述

SVD降维

SVD应用

1.图像降维

我们可以将图像存储在一个矩阵中。每个图像由一组像素组成，这些像素是构成该图像的基本单元。每个像素表示图像中某个特定位置的颜色或光强度。在 PNG 格式的灰度图像中，每个像素的值在 0 和 1 之间，其中 0 对应黑色，1 对应白色。因此一个具有mxn像素的灰度图像可以存储在一个 mxn 矩阵或 NumPy 数组中。在这里，我们使用 imread() 函数加载一张爱因斯坦的灰度图像，其分辨率为 480 × 423 像素，并将其存储为一个二维数组。然后，我们使用 SVD 来分解该矩阵，并使用前 30 个奇异值来重建图像。

# Reading the image
mat = plt.imread("Picture.png")

# SVD 
U, s, VT = LA.svd(mat)

Sigma = np.zeros((mat.shape[0], mat.shape[1]))
Sigma[:min(mat.shape[0], mat.shape[1]), :min(mat.shape[0], mat.shape[1])] = np.diag(s)

# Reconstruction of the matrix using the first 30 singular values
k = 30  # 取前30个特征值，重建图像
mat_approx = U[:, :k] @ Sigma[:k, :k] @ VT[:k, :]

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,8))
plt.subplots_adjust(wspace=0.3, hspace=0.2)

ax1.imshow(mat, cmap='gray')
ax1.set_title("Original image")

ax2.imshow(mat_approx, cmap='gray')
ax2.set_title("Reconstructed image using the \n first {} singular values".format(k))
plt.show()

image-20240810223105337

左图是原始图像，右图是取前30个奇异值重建的图像。

fig, axes = plt.subplots(2, 3, figsize=(10,8))
plt.subplots_adjust(wspace=0.3, hspace=0.2)

for i in range(0, 6):
    mat_i = s[i] * U[:,i].reshape(-1,1) @ VT[i,:].reshape(1,-1)
    axes[i // 3, i % 3].imshow(mat_i)
    axes[i // 3, i % 3].set_title("$\sigma_{0}\mathbf{{u_{0}}}\mathbf{{v_{0}}}^T$".format(i+1), fontsize=16)
    
plt.show()

image-20240810223817059

请注意，与原始灰度图像不同，这些秩为 1 的矩阵的元素值可以大于 1 或小于 0，它们不应被解释为灰度图像。因此，我没有使用 cmap='gray' 也没有将它们显示为灰度图像。在绘制这些矩阵时，我们不关心像素的绝对值，而是关注它们之间的相对值。

为了理解每个矩阵中如何存储图像信息，我们可以研究一个更简单的图像。下面的代码读取了一个包含五个简单形状的二进制图像（一个矩形和 4 个圆形），并取前2个、前4个和前6个奇异值重建图像。

mat = plt.imread("shapes.png")

# SVD 
U, s, VT = LA.svd(mat)

Sigma = np.zeros((mat.shape[0], mat.shape[1]))
Sigma[:min(mat.shape[0], mat.shape[1]), :min(mat.shape[0], mat.shape[1])] = np.diag(s)

fig, axes = plt.subplots(2, 2, figsize=(10,8))
plt.subplots_adjust(wspace=0.3, hspace=0.2)

axes[0, 0].imshow(mat, cmap='gray')
axes[0, 0].set_title("Original image")

for i in range(1, 4):
    k = i * 2
    # Reconstruction of the matrix using the first k singular values
    mat_approx = U[:, :k] @ Sigma[:k, :k] @ VT[:k, :]

    axes[i // 2, i % 2].imshow(mat_approx, cmap='gray')
    axes[i // 2, i % 2].set_title("Reconstructed image using the \n first {} singular values".format(k))

plt.show()

image-20240810224420549

现在我们绘制前 6 个奇异值对应的矩阵：

fig, axes = plt.subplots(2, 3, figsize=(10,6))
plt.subplots_adjust(wspace=0.3, hspace=0.05)

for i in range(0, 6):
    mat_i = s[i] * U[:,i].reshape(-1,1) @ VT[i,:].reshape(1,-1)
    #mat_i[mat_i < 1e-8] = 0
    axes[i // 3, i % 3].imshow(mat_i)
    axes[i // 3, i % 3].set_title("$\sigma_{0}\mathbf{{u_{0}}}\mathbf{{v_{0}}}^T$".format(i+1), fontsize=16)
    
plt.show()

image-20240810224526788

从每个奇异值的图像可以看出，前两个奇异值构建的矩阵大致捕捉了四个圆形作为四个矩形，而最后四个奇异值构建的矩阵则添加了更多细节。

如下图，图像中有 6 根柱子，第一个奇异值对应的矩阵可以捕捉到原始图像中的柱子数量。

image-20240810224851967

2.特征脸

在这个示例中，我们将使用 Scikit-learn 库中的 Olivetti faces 数据集。该数据集包含 400 张图像，图像展示了 40 位不同受试者的面部特征。对于一些受试者，图像是在不同时间拍摄的，具有不同的光照、面部表情和面部细节。这些图像为灰度图像，每张图像具有 64×64 像素。每个像素的强度值在区间 [0, 1] 上。首先，我们加载数据集：

data = fetch_olivetti_faces()
imgs = data.images
print(imgs.shape)

我们调用它来读取数据，并将图像存储在 imgs 数组中。这是一个形状为 (400, 64, 64) 的数组，包含 400 张 64×64 的灰度图像。我们可以在这里展示其中的一些作为示例：

fig, axes = plt.subplots(1, 5, figsize=(14, 8))
plt.subplots_adjust(wspace=0.4)

for i in range(0, 5):
    axes[i].imshow(imgs[i*40], cmap='gray')

plt.show()

image-20240810225330410

在前一个示例中，我们将原始图像存储在一个矩阵中，然后使用 SVD 对其进行分解。在这里，我们采用另一种方法。我们有400 张图像并给每张图像分配一个从 1 到 400 的标签。现在我们使用独热编码来表示这些标签，使用一个包含 400 个元素的列向量。对于每个标签 k，除了第 k个元素，所有元素都是零。因此，标签 k将由以下向量表示：

现在我们将每张图像存储在一个列向量中。每张图像有 64 × 64 = 4096 个像素。因此，我们可以将每张图像展平并将像素值放入一个具有 4096 个元素的列向量f中，如下图所示：

mage-20240810230051322

我们计算并显示前8个奇异值对应的U向量：

U, s, VT = LA.svd(M)

fig, axes = plt.subplots(2, 4, figsize=(10,6))
plt.subplots_adjust(wspace=0.3, hspace=0.1)

for i in range(0, 8):
    axes[i // 4, i % 4].imshow(U[:, i].reshape((64,64)))
    axes[i // 4, i % 4].set_title("$\mathbf{{u_{0}}}$".format(i+1), fontsize=16)
    
plt.show()

image-20240810231222673

x= np.zeros((400, 1))
x[160, 0] = 1
Sigma = np.zeros((M.shape[0], M.shape[1]))
Sigma[:min(M.shape[0], M.shape[1]), :min(M.shape[0], M.shape[1])] = np.diag(s)

fig, axes = plt.subplots(2, 4, figsize=(14, 8))
fig.suptitle("Reconstructed image using the first k singular values", fontsize=16)
plt.subplots_adjust(wspace=0.3, hspace=0.1)

axes[0, 0].imshow(imgs[160], cmap='gray')
axes[0, 0].set_title("Original image")

k_list = [1, 6, 10, 15, 20, 35, 80]
for i in range(1, 8):
    # Reconstruction of the matrix using the first k singular values
    k = k_list[i-1] 
    mat_approx = U[:, :k] @ Sigma[:k, :k] @ VT[:k, :] @ x

    axes[i // 4, i % 4].imshow(mat_approx.reshape((64,64)), cmap='gray')
    axes[i // 4, i % 4].set_title("k = {}".format(k))

plt.show()

这里我们使用不同个数的奇异值去重建第160号图像（k表示使用前k个奇异值重建图像）。

image-20240810232014008

上图反映的信息与特征脸步骤的细节基本一致，如最高奇异值（k=1）重建的图像反映了眼睛信息，同样地，第一个特征脸反映的也是眼睛信息。

3.降低噪声

SVD可以用于降低图像中的噪声。

mat = plt.imread("text.png")

# Adding noise
noise = np.random.rand(mat.shape[0], mat.shape[1])
mat[noise > 0.95] = 0

# SVD 
U, s, VT = LA.svd(mat)

Sigma = np.zeros((mat.shape[0], mat.shape[1]))
Sigma[:min(mat.shape[0], mat.shape[1]), :min(mat.shape[0], mat.shape[1])] = np.diag(s)

fig, axes = plt.subplots(2, 2, figsize=(10,8))
plt.subplots_adjust(wspace=0.2, hspace=0.1)

axes[0, 0].imshow(mat, cmap='gray')
axes[0, 0].set_title("Original image", y=1.08)

k_list = [20, 55, 200]
for i in range(1, 4):
    k = k_list[i-1]
    mat_rank_k = U[:, :k] @ Sigma[:k, :k] @ VT[:k, :]
    axes[i // 2, i % 2].imshow(mat_rank_k, cmap='gray')
    axes[i // 2, i % 2].set_title("Reconstructed image using the \n first {} singular values".format(k), y=1.08)

plt.show()

image-20240811135427217

首先，我们加载图像并添加一些噪声。然后使用前 20、55 和 200 个奇异值重建图像。如上图所示，随着重建矩阵的秩增加，噪声的量也随之增加。因此如果我们使用较低的秩，比如 20，可以显著减少图像中的噪声。

结论

我真的觉得奇异值分解（SVD）被低估了。它是线性代数中一个非常重要的基础概念，而且它的应用非常酷！相信我，我们看到的只是 SVD 众多用途中的一小部分。有什么问题，欢迎讨论！

THE END !

文章结束，感谢阅读。您的点赞，收藏，评论是我继续更新的动力。大家有推荐的公众号可以评论区留言，共同学习，一起进步。