其他算法-SVD奇异值分解

奇异值分解即 Singular Value Decomposition，简称SVD，SVD广泛应用于数据分析的降维处理

SVD原理

特征值分解EVD

开始了解SVD前应当先了解特征值（本征值）分解，特征值分解 Eigen Value Decomposition，简称EVD；

首先回顾特征值和特征向量，对于 $Ax=\lambda x$，$A$是一个$n\times n$的矩阵，$x$是一个$n$维向量，则$\lambda$是矩阵$A$的一个特征值，而$x$被称为特征值$\lambda$对应的特征向量；

假设矩阵$A$是一个$m\times m$的实对称矩阵（即满足$A^T=A$），则矩阵$A$可以被分解成：

其中，$Q$为标准正交阵（即满足$Q{Q}^T=I$），$\sum$为对角矩阵，$Q$与$\sum$的维度均为$m\times m$，${\lambda }_i$为特征值，$Q$中的列向量$q_i$为特征值${\lambda }_i$对应的特征向量；

---

在线性代数中，单位矩阵可以用$I$或者$E$表示

---

对于分解结果，有以下关系：
$$A{q}_i={\lambda }_i{q}_i,\left\{q_i^T{q}_j=0|i\ne j\right\}$$

奇异值分解SVD

在实际问题中，很难遇到一个矩阵是实对称的情况，大多数的矩阵都是一般的形状为$m\times n$的矩阵，对于这样的矩阵，依然想把它分解成类似EVD的形式：
$$A=U\Sigma V^T$$
其中，$U,V$均为标准正交阵（$U{U}^T=I,V{V}^T=I$），$U$称为左奇异矩阵；$V$称为右奇异矩阵；$\sum$仅在主对角线上有值，这些值被称为奇异值，其他位置的元素均为0；

矩阵的维度分别为：
$$U\in {\mathbb{R}}^{m\times m},\Sigma \in {\mathbb{R}}^{m\times n},V\in {\mathbb{R}}^{n\times n}$$
一般，$\sum$的形式如下：

关于奇异值的计算
根据标准正交阵的性质，有以下化简过程：
$$A{A}^T=U\Sigma V^{T}V{\Sigma }^T{U}^T=U\Sigma {\Sigma }^T{U}^T$$
$$A^{T}A=V{\Sigma }^T{U}^{T}U\Sigma V^T=V{\Sigma }^T\Sigma V^T$$
注意，$\Sigma {\Sigma }^T\in {\mathbb{R}}^{m\times m}$与${\Sigma }^T\Sigma \in {\mathbb{R}}^{n\times n}$是不相等的，但它们在主对角线上的值是相等的，即：

由于$A{A}^T$和$A^{T}A$是对称阵，可以利用EVD对其进行分解，对$\Sigma {\Sigma }^T$分解得到特征矩阵$U$，对${\Sigma }^T\Sigma$分解得到特征矩阵$V$，对$\Sigma {\Sigma }^T$或者${\Sigma }^T\Sigma$开方即得到所有奇异值

---

演示：$A{A}^T$和$A^{T}A$是对称阵

---

SVD的性质

对于奇异值，这与EVD中的特征值类似，对矩阵$\Sigma$的奇异值按照大小排序，可以发现奇异值大小减少得很快，在很多情况下，前10%甚至1%的奇异值的和就占了全部的奇异值之和的99%以上的比例；也就是说，我们也可以用最大的 $k$ 个的奇异值和对应的左右奇异向量来近似描述矩阵：
$$A_{m\times n}=U_{m\times m}{\Sigma }_{m\times n}{V}_{n\times n}^T\approx U_{m\times k}{\Sigma }_{k\times k}{V}_{k\times n}^T$$
即原本很大的矩阵，现在可以用三个较小的矩阵近似：

基于这个性质，SVD常被用于数据降维，实现数据压缩和数据去噪

基于Numpy的SVD图像压缩

由于numpy已提供了SVD，因此，可以通过numpy直接使用SVD；首先导入图像，并查看信息：

import matplotlib.pyplot as plt
# 用于读入图像
import matplotlib.image as mpimg
import numpy as np

img_eg = mpimg.imread("./caffe.jpg")

print(img_eg.shape)
print(type(img_eg))
# (500, 333, 3)
# <class 'numpy.ndarray'>

将RGB图像的三通道数据转为矩阵，再使用SVD分解：

img_temp = img_eg.reshape(img_eg.shape[0],-1)
print(img_temp.shape) # (500, 999)

U,Sigma,VT = np.linalg.svd(img_temp)

# numpy的Sigma矩阵直接用奇异值组成向量表示
print(U.shape,Sigma.shape,VT.shape) # (500, 500) (500,) (999, 999)

---

numpy.diag用于提取对角线或构造对角线数组，补充numpy.diag的操作：

• 在 np.diag(array) 中，array是一个1维数组时，结果形成一个以一维数组为对角线元素的矩阵；
• array是一个二维矩阵时，结果输出矩阵的对角线元素

>>> a
array([1, 2, 3])
>>> b
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])
       
>>> np.diag(a)
array([[1, 0, 0],
       [0, 2, 0],
       [0, 0, 3]])
>>> np.diag(b)
array([1, 5, 9])

---

提取前部分的奇异值重构图像：

# 取前60个奇异值
sval_nums = 60
img_restruct1 = (U[:, 0:sval_nums]).dot(np.diag(Sigma[0:sval_nums])).dot(VT[0:sval_nums, :])
img_restruct1 = img_restruct1.reshape(img_eg.shape[0],img_eg.shape[1],img_eg.shape[2])

# 取前120个奇异值
sval_nums = 120
img_restruct2 = (U[:, 0:sval_nums]).dot(np.diag(Sigma[0:sval_nums])).dot(VT[0:sval_nums, :])
img_restruct2 = img_restruct2.reshape(img_eg.shape[0],img_eg.shape[1],img_eg.shape[2])

可视化结果：

fig, ax = plt.subplots(1, 3, figsize=(14, 22))
ax[0].imshow(img_eg)
ax[0].set(title="Original image")
ax[1].imshow(img_restruct1.astype(np.uint8))
ax[1].set(title="nums of sigma = 60")
ax[2].imshow(img_restruct2.astype(np.uint8))
ax[2].set(title="nums of sigma = 120")
plt.show()