线性代数是机器学习领域的基础，其中一个最重要的概念是奇异值分解（SVD），本文尽可能简洁的介绍SVD（奇异值分解）算法的基础理解，以及它在现实世界中的应用。
SVD是最广泛使用的无监督学习算法之一，它在许多推荐系统和降维系统中居于核心位置，这些系统是全球公司如谷歌、Netflix、Facebook、YouTube等的核心技术。
简单来说，SVD是将一个任意矩阵分解为三个矩阵。所以如果我们有一个矩阵A，那么它的SVD可以表示为：
A是矩阵，U是的正交矩阵，，的正交矩阵。
U也被称为左奇异向量，S为奇异值，V为右奇异向量。
带维度的奇异值分解：
用矩阵表示奇异值分解：
我们通常将具有较大特征值的向量排列在前，而较小特征值的向量则排在后面。
特征值与向量的对应关系：
与特征值分解相比，奇异值分解可以应用于非方阵。在SVD中，U和 V 对于任何矩阵都是可逆的，并且它们是正交归一的，这是我们所喜爱的特性。虽然这里不进行证明，但我们可以告诉你，奇异值比特征值在数值上更稳定。
为了更好地理解，我们通过一个例子演示SVD
假设我们有非方阵A：
我们计算矩阵与转置矩阵的乘积，有：
求解的特征值和特征向量：
求解的特征值和特征向量：
奇异值是正特征值的平方根，即5和3。因此非方阵A的SVD分解为：
SVD分解证明
为了得到矩阵A的SVD分解：
我们需要求解正交矩阵U和V
其中I是单位矩阵。
上面3个方程有3个未知数，因此是可以求解这3个未知数的。
矩阵A的转置是：
由于正交矩阵V和U满足：
我们计算：
最后一个方程等价于求矩阵的特征向量，我们只需将所有特征向量放入一个矩阵中，矩阵S则是包含特征值的对角矩阵。
即，有:
其中是的特征值。
计算矩阵的表达式：
其中矩阵V包含了的所有特征向量，矩阵S包含了的所有特征值的平方根。我们可以对重复相同的工程，并得到一个类似的方程。
现在我们知道了通过分别得到了矩阵U和V，其中U和V包含了各自的特征向量，S的特征值是矩阵特征值的平方根。因此3个未知量都求解了，即可得到矩阵A的SVD分解：
SVD的另一种表述
由于矩阵V是正交的，等于单位矩阵I。我们可以将SVD方程重新写为：
矩阵V：
其中：
矩阵U：
其中：
矩阵S：
我们取第1个特征向量：
得：
可以通用化上式：
因此A的SVD分解可以看做是的一系列外积：
SVD降维
若我们取前k个奇异值，根据上式得到新的矩阵。
矩阵是降维后的数值。若大家还有疑问，欢迎留言。
SVD应用
1.图像降维
我们可以将图像存储在一个矩阵中。每个图像由一组像素组成，这些像素是构成该图像的基本单元。每个像素表示图像中某个特定位置的颜色或光强度。在 PNG 格式的灰度图像中，每个像素的值在 0 和 1 之间，其中 0 对应黑色，1 对应白色。因此一个具有 像素的灰度图像可以存储在一个 矩阵或 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()
左图是原始图像，右图是取前30个奇异值重建的图像。
原始矩阵是 480×423 的，因此我们需要存储 480×423=203040 个值。经过 SVD 后，每个具有 480 个元素，每个 具有 423 个元素。为了使用前 30 个奇异值重建图像，我们只需要保留前 30 个 ，这意味着只需要存储 30×(1+480+423)=27120 个值。这大约是原始图像所需值数量的 13%。因此，通过使用 SVD，我们可以得到原始图像的良好近似，并节省大量内存。
下图计算了前 6 个奇异值对应的矩阵。每个矩阵 的秩为 1，并且具有与原始矩阵相同的行数和列数。
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()
请注意，与原始灰度图像不同，这些秩为 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()
现在我们绘制前 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()
从每个奇异值的图像可以看出，前两个奇异值构建的矩阵大致捕捉了四个圆形作为四个矩形，而最后四个奇异值构建的矩阵则添加了更多细节。
如下图，图像中有 6 根柱子，第一个奇异值对应的矩阵可以捕捉到原始图像中的柱子数量。
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()
在前一个示例中，我们将原始图像存储在一个矩阵中，然后使用 SVD 对其进行分解。在这里，我们采用另一种方法。我们有400 张图像并给每张图像分配一个从 1 到 400 的标签。现在我们使用独热编码来表示这些标签，使用一个包含 400 个元素的列向量。对于每个标签 k，除了第 k个元素，所有元素都是零。因此，标签 k将由以下向量表示：
现在我们将每张图像存储在一个列向量中。每张图像有 64 × 64 = 4096 个像素。因此，我们可以将每张图像展平并将像素值放入一个具有 4096 个元素的列向量中，如下图所示：
因此，每张标签为k的图像将存储在向量中，我们需要 400 个 向量来保存所有图像。现在我们定义一个变换矩阵，它将标签向量 转换为其对应的图像向量。这些 向量将作为矩阵 的列：
这个矩阵有 4096 行和 400 列。我们可以简单地使用 来找到每个标签对应的图像（可以是任何向量 ，而 将是对应的）。例如，对于这个数据集中的第三张图像，标签是 3， 的所有元素都是零，除了第三个元素是 1。现在记住分块矩阵的乘法，当我们用 乘以 时，除了第三列 之外， 的所有列都被乘以零，所以：
现在我们对矩阵进行SVD分解。
我们计算并显示前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()
上图就是我们要介绍的特征脸，每个特征脸捕捉了图像向量的一些信息。例如特征脸主要反映了眼睛信息，特征脸捕捉了鼻子信息。
上述代码的数组s具有400个元素，因此我们有400个非零的奇异值，矩阵的矩为400。因此我们需要前 400 个 向量来完全重建矩阵，使用这些基向量轻松重建其中一张图像。
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个奇异值重建图像）。
上图反映的信息与特征脸步骤的细节基本一致，如最高奇异值（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()
首先，我们加载图像并添加一些噪声。然后使用前 20、55 和 200 个奇异值重建图像。如上图所示，随着重建矩阵的秩增加，噪声的量也随之增加。因此如果我们使用较低的秩，比如 20，可以显著减少图像中的噪声。
结论
我真的觉得奇异值分解（SVD）被低估了。它是线性代数中一个非常重要的基础概念，而且它的应用非常酷！相信我，我们看到的只是 SVD 众多用途中的一小部分。有什么问题，欢迎讨论！
欢迎扫码关注：
