机器学习理论《统计学习方法》学习笔记：奇异值分解（SVD）
奇异值分解（SVD）
- 摘要
- 1 奇异值分解的定义与定理
- 1.1 奇异值分解的定义
- 1.2 奇异值分解的基本定理
- 1.3 奇异值分解的几何解释
- 2 紧奇异值分解和截断奇异值分解
- 2.1 紧奇异值分解
- 2.2 截断奇异值分解
- 3 奇异值分解的计算
- 4 特征值分解和奇异值分解
- 4.1 特征值分解（EIG）
- 4.1.1 特征值分解的定义
- 4.1.2 使用Python实现特征值分解
- 4.2 奇异值分解（SVD）
- 4.2.1 完全奇异值分解
- 4.2.2 紧奇异值分解
- 4.2.3 在python中实现奇异值分解
- 情况1：完全奇异值分解
- 情况2：紧奇异值分解
- 5 奇异值分解的应用
- 5.1 图像压缩
- 5.1.1 奇异值分解压缩的原理
- 5.1.2 Python实现图像压缩
- 5.2 图像去噪
- 5.2.1 什么是噪声
- 5.2.2 椒盐噪声
- 5.2.3 高斯噪声
- 5.2.4 奇异值分解去噪
摘要
PCA的实现一般有两种，一种是用特征值分解去实现的，一种是用奇异值分解去实现的。在上篇文章中便是基于特征值分解的一种解释。特征值和奇异值在大部分人的印象中，往往是停留在纯粹的数学计算中。而且线性代数或者矩阵论里面，也很少讲任何跟特征值与奇异值有关的应用背景。
奇异值分解是一个有着很明显的物理意义的一种方法，它可以将一个比较复杂的矩阵用更小更简单的几个子矩阵的相乘来表示，这些小矩阵描述的是矩阵的重要的特性。就像是描述一个人一样，给别人描述说这个人长得浓眉大眼，方脸，络腮胡，而且带个黑框的眼镜，这样寥寥的几个特征，就让别人脑海里面就有一个较为清楚的认识，实际上，人脸上的特征是有着无数种的，之所以能这么描述，是因为人天生就有着非常好的抽取重要特征的能力，让机器学会抽取重要的特征，SVD是一个重要的方法。
奇异值分解（SVD）是一种矩阵因子分解方法，是线性代数的概念，但在统计学习中被广泛使用，成为其重要工具，其中主成分分析、潜在语义分析都用到奇异值分解。
任意一个m×nm\times nm×n矩阵，都可以表示为三个矩阵的乘积（因子分解）形式，分别是m阶正交矩阵、由降序排列的非负的对角线元素组成的m×nm\times nm×n矩形对角矩阵、n阶正交矩阵，称为该矩阵的奇异值分解。矩阵的奇异值分解一定存在，但不唯一。奇异值分解可以看作是矩阵数据压缩的一种方法，即用因子分解的方式近似地表示原始矩阵，这种近似是在平方损失意义下的最优近似。
1 奇异值分解的定义与定理
1.1 奇异值分解的定义
矩阵的奇异值分解是指将一个非零的m×nm\times nm×n实矩阵AAA，A∈Rm×nA\in R^{m\times n}A∈Rm×n，表示为以下三个实矩阵乘积形式的运算，即进行矩阵的因子分解：A=UΣVTA=U\Sigma V^{T}A=UΣVT，其中U是m阶正交矩阵，V是n阶正交矩阵，Σ\SigmaΣ是由降序排列的非负的对角线元素组成的m×nm\times nm×n矩形对角矩阵，满足
UUT=IUU^{T}=IUUT=I
VVT=IVV^{T}=IVVT=I
Σ=diag(σ1,σ2,⋯,σp)\Sigma=diag(\sigma_1,\sigma_2,\cdots,\sigma_p)Σ=diag(σ1,σ2,⋯,σp)
σ1≥σ2≥⋯≥σp≥0\sigma_1\ge\sigma_2\ge\cdots\ge\sigma_p\ge0σ1≥σ2≥⋯≥σp≥0
p=min(m,n)p=min(m,n)p=min(m,n)
A=UΣVTA=U\Sigma V^{T}A=UΣVT称为矩阵A的奇异值分解，σi\sigma_iσi称为矩阵A的奇异值，U的列向量称为A的左奇异向量，V的列向量称为右奇异向量。
1.2 奇异值分解的基本定理
若AAA为一m×nm\times nm×n实矩阵，A∈Rm×nA\in R^{m\times n}A∈Rm×n，则AAA的奇异值分解存在A=UΣVTA=U\Sigma V^TA=UΣVT，其中U是m阶正交矩阵，V是n阶正交矩阵，Σ\SigmaΣ是m×nm\times nm×n矩形对角矩阵，其对角线元素非负，且按降序排列。
任意给定一个实矩阵，奇异值分解矩阵是否一定存在呢？由奇异值分解的基本定理可知，答案是肯定的。
1.3 奇异值分解的几何解释
从线性变换的角度理解奇异值分解，m×nm\times nm×n矩阵A表示从n维空间RnR^nRn到m维空间RmR^mRm的一个线性变换，T:x→AxT:x\rightarrow AxT:x→Ax，x∈Rn,Ax∈Rmx\in R^n,Ax\in R^mx∈Rn,Ax∈Rm，x和Ax分别是各自空间的向量。
线性变换可以分解为三个简单的变换：一个坐标系的旋转或反射变换、一个坐标轴的缩放变换、另一个坐标系的旋转或反射变换。奇异值定理保证这种分解一定存在。
2 紧奇异值分解和截断奇异值分解
奇异值分解定义给出的奇异值分解A=UΣVTA=U\Sigma V^TA=UΣVT，又称为矩阵的完全奇异值分解。实际常用的是奇异值分解的紧凑形式和截断形式。紧奇异值分解是与原始矩阵等秩的奇异值分解，截断奇异值分解是比原始矩阵低秩的奇异值分解。
2.1 紧奇异值分解
设有m×nm\times nm×n实矩阵AAA，其秩为rank(A)=r,r≤min(m,n)rank(A)=r,r\le min(m,n)rank(A)=r,r≤min(m,n)，则称UrΣrVrTU_r\Sigma_r V_r^TUrΣrVrT为A的紧奇异值分解，即A=UrΣrVrTA=U_r\Sigma_r V_r^TA=UrΣrVrT
，其中UrU_rUr是m×nm\times nm×n矩阵，VrV_rVr是n×rn\times rn×r矩阵，Σr\Sigma_rΣr是rrr阶对角矩阵；矩阵UrU_rUr由完全奇异值分解中U的前r列、矩阵VrV_rVr由V的前r列、矩阵Σr\Sigma_rΣr由Σ\SigmaΣ的前r个对角线元素得到。紧奇异值分解的对角矩阵Σr\Sigma_rΣr的秩与原始矩阵A的秩相等。
2.2 截断奇异值分解
在矩阵的奇异值分解中，只取最大的k个奇异值（k<r,rk<r,rk<r,r为矩阵的秩）对应的部分，就得到矩阵的截断奇异值分解。实际应用中提到矩阵的奇异值分解时，通常指截断奇异值分解。
设A是m×nm\times nm×n实矩阵，其秩为rank(A)=r,0<k<rrank(A)=r,0<k<rrank(A)=r,0<k<r，则称UkΣkVkTU_k\Sigma_kV_k^TUkΣkVkT为矩阵A的截断奇异值分解A≈UkΣkVkTA\approx U_k \Sigma_k V_k^TA≈UkΣkVkT。其中UkU_kUk是m×km\times km×k矩阵，VkV_kVk是n×kn\times kn×k矩阵，Σk\Sigma_kΣk是kkk阶对角矩阵；矩阵UkU_kUk由完全奇异值分解中U的前k列，矩阵VkV_kVk由V的前k列，矩阵Σk\Sigma_kΣk由Σ\SigmaΣ的前k个对角线元素得到。对角矩阵Σk\Sigma_kΣk的秩比原始矩阵A的秩低。
3 奇异值分解的计算
矩阵A的奇异值分解可以通过求对称矩阵ATAA^TAATA的特征值和特征向量得到。ATAA^TAATA的特征向量构成正交矩阵V的列；ATAA^TAATA的特征值λj\lambda_jλj的平方根为奇异值σi\sigma_iσi，即：σj=λj,j=1,2,⋯,n\sigma_j=\sqrt\lambda_j,j=1,2,\cdots,nσj=λj,j=1,2,⋯,n对其由大到小排列作为对角线元素，构成对角矩阵Σ\SigmaΣ；求正奇异值对应的左奇异向量，再求扩充的ATA^TAT的标准正交基，构成正交矩阵U的列。从而得到A的奇异值分解A=UΣVTA=U\Sigma V^TA=UΣVT
4 特征值分解和奇异值分解
- 只有方阵才能进行特征值分解。
4.1 特征值分解（EIG）
4.1.1 特征值分解的定义
如果说一个向量vvv是方阵AAA的特征向量，即可以表示为下面形式Av=λvAv=\lambda vAv=λv，此时λ\lambdaλ称为特征向量vvv对应的特征值，一个矩阵的一组特征向量是一组正交向量，
特征值分解是将一个矩阵分解为下面形式：A=QΣQ−1A=Q\Sigma Q^{-1}A=QΣQ−1，其中Q是该矩阵A的特征向量组成的矩阵，Σ\SigmaΣ是一个对角阵，每个对角线上的元素就是一个特征值。首先，要明确的是，一个矩阵其实就是一个线性变换，因为一个矩阵乘以一个向量后得到的向量，其实就相当于将这个向量进行了线性变换。
4.1.2 使用Python实现特征值分解
Numpy中的linalg已经实现了ELG，可以直接调用，具体为：
e_vals, e_vecs=np.linalg.eig(a)
输入参数：a为需要分解的参数
返回：
- e_vals: 由特征值构成的向量
- e_vecs: 由特征向量构成的矩阵
注意矩阵求逆可以使用np.linalg.inv(a)
import numpy as npa = np.random.randn(4, 4)
e_vals, e_vecs = np.linalg.eig(a)
print('矩阵分解得到的形状：\n', e_vals.shape, e_vecs.shape)
print('特征值：\n', e_vals)
print('特征向量矩阵：\n', e_vecs)
# smat = np.zeros((4, 4))
smat = np.diag(e_vals)
print('特征值对角矩阵：\n', smat)
# 验证特征值分解
# 对比两个矩阵的各个元素，一致则返回true
result = np.allclose(a, np.dot(e_vecs, np.dot(smat, np.linalg.inv(e_vecs))))
print('验证特征值分解：\n', result)输出结果
矩阵分解得到的形状：(4,) (4, 4)
特征值：[-1.15983308+0.j 1.47606642+1.78459968j 1.47606642-1.78459968j 0.52887458+0.j ]
特征向量矩阵：[[-0.88012147+0.j 0.13851919+0.00751309j 0.13851919-0.00751309j 0.27291217+0.j ][-0.0313536 +0.j -0.13663243-0.01138874j -0.13663243+0.01138874j -0.81808112+0.j ][-0.24771796+0.j -0.02890239+0.57829333j -0.02890239-0.57829333j 0.34163481+0.j ][-0.40378083+0.j -0.79164344+0.j -0.79164344-0.j -0.37356109+0.j ]]
特征值对角矩阵：[[-1.15983308+0.j 0. +0.j 0. +0.j 0. +0.j ][ 0. +0.j 1.47606642+1.78459968j 0. +0.j 0. +0.j ][ 0. +0.j 0. +0.j 1.47606642-1.78459968j 0. +0.j ][ 0. +0.j 0. +0.j 0. +0.j 0.52887458+0.j ]]
验证特征值分解：True
4.2 奇异值分解（SVD）
若有一个矩阵A，对其进行奇异值分解可以得到三个矩阵：A=UΣVTA=U\Sigma V^{T}A=UΣVT
4.2.1 完全奇异值分解
- 当进行完全奇异值分解时，三个矩阵的大小为下图所示
矩阵Sigma（上图U和V之间的矩阵）除了对角元素不为0，其他元素都为0，并且对角元素是从大到小排列的，这些对角元素就是奇异值。
4.2.2 紧奇异值分解
- Sigma中有n个奇异值，由于排在后面的大多数接近于0，所以仅保留比较大的r个奇异值，此时称为紧奇异值分解。
实际应用中，我们仅需要保留三个比较小的矩阵就可以表示A，不仅节省存储量，更减少了计算量
高清矢量图下载地址：https://download.csdn.net/download/qq_40507857/13124280
4.2.3 在python中实现奇异值分解
在Numpy中已经实现了SVD，可以直接调用，具体为：
U,S,Vh=np.linalg.svd(a,full_matrices=True,compute_uv=True)
输入参数：
a：要分解的矩阵，维数大于等于2
full_matrices：bool值，默认为True，此时为完全奇异值分解；若为False，此时为紧奇异值分解。
compute_uv：bool值，表示除了S之外，是否计算U和Vh，默认为True，即结果返回三个矩阵
返回值：完全奇异值分解的三个矩阵
情况1：完全奇异值分解
import numpy as npa = np.random.randn(9, 6)
u, s, vh = np.linalg.svd(a, full_matrices=True, compute_uv=True)
print('完全奇异值分解得到的形状：')
print('U:', u.shape, 'S:', s.shape, 'Vh:', vh.shape)
print('奇异值：\n', s)
smat = np.zeros((9, 6))
smat[:6, :6] = np.diag(s)
print('奇异矩阵：\n', smat)
# 验证奇异值分解
result = np.allclose(a, np.dot(u, np.dot(smat, vh)))
print('验证完全奇异值分解', result)输出结果
完全奇异值分解得到的形状：
U: (9, 9) S: (6,) Vh: (6, 6)
奇异值：[5.22112777 4.0473851 3.1832999 2.44399385 2.31411792 1.76042697]
奇异矩阵：[[5.22112777 0. 0. 0. 0. 0. ][0. 4.0473851 0. 0. 0. 0. ][0. 0. 3.1832999 0. 0. 0. ][0. 0. 0. 2.44399385 0. 0. ][0. 0. 0. 0. 2.31411792 0. ][0. 0. 0. 0. 0. 1.76042697][0. 0. 0. 0. 0. 0. ][0. 0. 0. 0. 0. 0. ][0. 0. 0. 0. 0. 0. ]]
验证完全奇异值分解 True情况2：紧奇异值分解
import numpy as npa = np.random.randn(9, 6)
u, s, vh = np.linalg.svd(a, full_matrices=False, compute_uv=True)
print('紧奇异值分解得到的形状：')
print('U:', u.shape, 'S:', s.shape, 'Vh:', vh.shape)
print('奇异值：\n', s)
smat = np.zeros((6, 6))
smat=np.diag(s)
print('奇异矩阵：\n', smat)
# 验证奇异值分解
result = np.allclose(a, np.dot(u, np.dot(smat, vh)))
print('验证完全奇异值分解', result)输出结果
紧奇异值分解得到的形状：
U: (9, 6) S: (6,) Vh: (6, 6)
奇异值：[5.49788835 4.6860516 3.80953519 3.4174992 2.45542127 0.80275215]
奇异矩阵：[[5.49788835 0. 0. 0. 0. 0. ][0. 4.6860516 0. 0. 0. 0. ][0. 0. 3.80953519 0. 0. 0. ][0. 0. 0. 3.4174992 0. 0. ][0. 0. 0. 0. 2.45542127 0. ][0. 0. 0. 0. 0. 0.80275215]]
验证完全奇异值分解 True5 奇异值分解的应用
图像压缩（image compression）：较少的奇异值就可以表达出图像中大部分信息，舍弃掉一部分奇异值来实现压缩。
图像降噪（image denoise）：噪声一般存在于图像高频部分，也表现在奇异值小的部分，故可以借助SVD实现去噪。
音频滤波（filtering）：Andrew Ng的机器学习课程上有个svd将混杂声音分离的例子，其实和噪声滤波类似。
求任意矩阵的伪逆（pseudo-inverse）：由于奇异矩阵或非方阵矩阵不可求逆，在特殊情况下需要广义求逆时可用svd方法。
模式识别（pattern recognition）：特征为矩阵，数据量较大时，可以用svd提取主要的成分。
潜在语义索引（Latent Semantic Indexing）：NLP中，文本分类的关键是计算相关性，这里关联矩阵A=USV’，分解的三个矩阵有很清楚的物理含义，可以同时得到每类文章和每类关键词的相关性。
5.1 图像压缩
5.1.1 奇异值分解压缩的原理
先看一个简单的例子，如果你想要在网络上给别人发送一段数据，数据的内容为
A=(1111111111111111111111111)A= \begin{pmatrix} 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ \end{pmatrix} A=⎝⎜⎜⎜⎜⎛1111111111111111111111111⎠⎟⎟⎟⎟⎞
当然，最简单的方法就是给这个矩阵直接发过去，这是一个5x5的矩阵，你至少需要发送25个数字。
但是我们可以把这个矩阵分解为两个矩阵的乘积，这样只需要发送10个数字。
A=(11111)(11111)A= \begin{pmatrix} 1\\ 1\\ 1\\ 1\\ 1\\ \end{pmatrix} \begin{pmatrix} 1&1&1&1&1\\ \end{pmatrix} A=⎝⎜⎜⎜⎜⎛11111⎠⎟⎟⎟⎟⎞(11111)
图像也可以被视为矩阵，图像的每一个点都是由RGB值定义的，所以每个图像可以被表示为三个巨型矩阵（分别是R，G，B矩阵）。
但图像所生成的矩阵显然不会像上面的例子那样简单的就被分解了。想要分解任意矩阵，这就需要用到SVD了。
SVD分解可以被认为是EVD（Eigen Value Decomposition 特征值分解）的延伸。特征值分解将一个矩阵分解为两组正交的特征向量和一个特征值对角线矩阵。
而特征值矩阵又是从大到小排列的，特征值大小的下降速度很快，我们可以通过丢弃一些特征值来压缩数据。对于压缩图像来说，只要人眼不可察觉便可以认为是成功的压缩。
简单来说，就是通过把一块大的数据分解为很多项，通过给数据的每个项的重要程度排序，挑选出一部分最重要的保留，丢弃一部分最不重要的，来实现数据压缩。
5.1.2 Python实现图像压缩
在python中使用SVD算法很容易，直接使用库函数即可。这里主要使用numpy库用来进行矩阵计算，matplotlib用来显示图像以及PIL库用来读取本地测试图片。
首先需要把测试图片导入进来，转换为numpy的矩阵。
img=Image.open(filename)
a=np.array(img)
这里图片转矩阵后的格式实际上是一个图片长乘宽的矩阵，这个矩阵的每一个项都包含3个数字，分别是R,G,B的值
SVD分解只需要一句话即可
u,sigma,v=np.linalg.svd(a[:,:,0])
这里的SVD分解返回三个执行后返回三个矩阵，分别是u,sigma和v
实现重建函数
# 重建图像函数
def rebuild_img(u, sigma, v, p):"""p为使用特征值的比例。通过改变p来比较特征值比例对图像的影像:param u::param sigma::param v::param p::return:"""# 首先计算出m和n，即图片矩阵的长和宽，然后创建一个零矩阵a作为组装场地。m = len(u)n = len(v)a = np.zeros((m, n))# count是所有特征值加起来的总和，用于后面计算比例使用count = (int)(sum(sigma))curSum = 0k = 0while curSum <= count * p:# uk和vk就是从参数u和v中取出，改变形式后形成的与当前特征值对应得一组特征向量。uk = u[:, k].reshape(m, 1)vk = v[k].reshape(1, n)# 不断地从参数中取出uk、vk和sigma，运算后叠加到a上去，直到满足一定的比例。a += sigma[k] * np.dot(uk, vk)curSum += sigma[k]k += 1a[a < 0] = 0a[a > 255] = 255return np.rint(a).astype(dtype=np.int32)
重建函数接受4个参数，u,sigma,v即重建矩阵所需的内容，p则为使用特征值的比例，我们将通过改变比例p来看使用特征值比例对画面的影响。
算法的步骤如下描述：
首先计算出m和n，即图片矩阵的长和宽，然后创建一个零矩阵a作为组装场地。
count是所有特征值加起来的总和，用于后面计算比例使用
uk和vk就是从参数u和v中取出，改变形式后形成的与当前特征值对应得一组特征向量。
然后不断地从参数中取出uk、vk和sigma，运算后叠加到a上去，直到满足一定的比例。
最后把所有矩阵内的项取整数退出即可。
有了分解与重建，现在可以设计数据压缩试验了。
这里我们控制特征值的使用比例，从0.1到1，每次步进0.1，然后分解重建，看看图像的显示情况。
for i in np.arange(0.1,1,0.1):u,sigma,v=np.linalg.svd(a[:,:,0])R=rebuild_img(u,sigma,v,i)u,sigma,v=np.linalg.svd(a[:,:,1])G=rebuild_img(u,sigma,v,i)u,sigma,v=np.linalg.svd(a[:,:,2])B=rebuild_img(u,sigma,v,i)I=np.stack((R,G,B),2)plt.subplot(330+i*10)plt.title(i)plt.imshow(I)plt.show()
【完整程序】
import numpy as np
from PIL import Image
import matplotlib.pyplot as plt# 重建图像函数
def rebuild_img(u, sigma, v, p):"""p为使用特征值的比例。通过改变p来比较特征值比例对图像的影像:param u::param sigma::param v::param p::return:"""# 首先计算出m和n，即图片矩阵的长和宽，然后创建一个零矩阵a作为组装场地。m = len(u)n = len(v)a = np.zeros((m, n))# count是所有特征值加起来的总和，用于后面计算比例使用count = (int)(sum(sigma))curSum = 0k = 0while curSum <= count * p:# uk和vk就是从参数u和v中取出，改变形式后形成的与当前特征值对应得一组特征向量。uk = u[:, k].reshape(m, 1)vk = v[k].reshape(1, n)# 不断地从参数中取出uk、vk和sigma，运算后叠加到a上去，直到满足一定的比例。a += sigma[k] * np.dot(uk, vk)curSum += sigma[k]k += 1a[a < 0] = 0a[a > 255] = 255return np.rint(a).astype(dtype=np.int32)def main():filepath = './dataset/images/gaoyuanyuan.jpg'# 首先需要把测试图片导入进来，转换为numpy的矩阵。img = Image.open(filepath)a = np.array(img)# 实现SVD分解for i in np.arange(0.1, 1.0, 0.1):u, sigma, v = np.linalg.svd(a[:, :, 0])R = rebuild_img(u, sigma, v, i)u, sigma, v = np.linalg.svd(a[:, :, 1])G = rebuild_img(u, sigma, v, i)u, sigma, v = np.linalg.svd(a[:, :, 2])B = rebuild_img(u, sigma, v, i)I = np.stack((R, G, B), 2)plt.subplot(330 + i * 10)title = int(i * 10) / 10plt.title(title)plt.imshow(I)plt.show()if __name__ == '__main__':main()【原图】来自于互联网，仅用于技术交流与分享
【不同比例的压缩图片】
可以看到，当sigma比例在0.5及以下时，能够明显察觉到图片被压缩的痕迹，但当sigma比例超过0.6时，细节的还原就比较好了，当0.7,0.8,0.9时，肉眼几乎无法发现压缩痕迹，证明了SVD作为图像压缩算法，在细节丢失方面是可以控制得比较好的。在保持细节的前提下，可以将数据压缩10%-30%左右。
下面程序实现单通道图像的压缩
import numpy as np
import matplotlib.pyplot as plt
from PIL import Imagedef svd_decompose(img, s_num):u, s, vt = np.linalg.svd(img)h, w = img.shape[:2]s_new = np.diag(s[:s_num], 0) # 用s_num个奇异值生成新对角矩阵u_new = np.zeros((h, s_num), float)vt_new = np.zeros((s_num, w), float)u_new[:, :] = u[:, :s_num]vt_new[:, :] = vt[:s_num, :]svd_img = u_new.dot(s_new).dot(vt_new)return svd_imgdef main():img = Image.open('./dataset/images/gaoyuanyuan.jpg') # (256,256)img = img.convert('L') # 转黑白图像img = np.array(img)print(img.shape)svd_decompose(img, 1)svd_1 = svd_decompose(img, 1)svd_5 = svd_decompose(img, 5)svd_10 = svd_decompose(img, 10)svd_20 = svd_decompose(img, 20)svd_50 = svd_decompose(img, 50)svd_100 = svd_decompose(img, 100)plt.figure(1)plt.subplot(331)plt.imshow(img, cmap='gray')plt.title('original')plt.xticks([])plt.yticks([])plt.subplot(332)plt.imshow(svd_1, cmap='gray')plt.title('1 Singular Value')plt.xticks([])plt.yticks([])plt.subplot(333)plt.imshow(svd_5, cmap='gray')plt.title('5 Singular Values')plt.xticks([])plt.yticks([])plt.subplot(335)plt.imshow(svd_10, cmap='gray')plt.title('10 Singular Values')plt.xticks([])plt.yticks([])plt.subplot(336)plt.imshow(svd_20, cmap='gray')plt.title('20 Singular Values')plt.xticks([])plt.yticks([])plt.subplot(338)plt.imshow(svd_50, cmap='gray')plt.title('50 Singular Values')plt.xticks([])plt.yticks([])plt.subplot(339)plt.imshow(svd_100, cmap='gray')plt.title('100 Singular Values')plt.xticks([])plt.yticks([])plt.show()if __name__ == '__main__':main()输出结果：同样，结果可见前50个特征就基本涵盖了原图所有信息。
5.2 图像去噪
5.2.1 什么是噪声
在噪声的概念中，通常采用信噪比（Signal-Noise Rate, SNR）衡量图像噪声。通俗的讲就是信号占多少，噪声占多少，SNR越小，噪声占比越大。
在信号系统中，计量单位为dB，为10lg(PS/PN), PS和PN分别代表信号和噪声的有效功率。在这里，采用信号像素点的占比充当SNR，以衡量所添加噪声的多少。
常见噪声有 椒盐噪声 和 高斯噪声 ，椒盐噪声可以理解为斑点，随机出现在图像中的黑点或白点；高斯噪声可以理解为拍摄图片时由于光照等原因造成的噪声。
5.2.2 椒盐噪声
椒盐噪声也称为脉冲噪声，是图像中经常见到的一种噪声，它是一种随机出现的白点（盐噪声）或者黑点（椒噪声），可能是亮的区域有黑色像素或是在暗的区域有白色像素，或是两者皆有。
成因：可能是影像讯号受到突如其来的强烈干扰而产生、类比数位转换器或位元传输错误等。例如失效的感应器导致像素值为最小值，饱和的感应器导致像素值为最大值。
图像模拟添加椒盐噪声原理：通过随机获取像素点并设置为高亮度点和低灰度点来实现的，简单说就是随机的将图像某些像素值改为0或255。
def sp_noise(image, prob):"""给图像加椒盐噪声:param image:图像:param prob:噪声比例:return:"""output = np.zeros(image.shape, np.uint8)thres = 1 - probfor i in range(image.shape[0]):for j in range(image.shape[1]):rdn = random.random()if rdn < prob:output[i][j] = 0elif rdn > thres:output[i][j] = 255else:output[i][j] = image[i][j]return output
5.2.3 高斯噪声
高斯噪声是指高斯密度函数服从高斯分布的一类噪声。特别的，如果一个噪声，它的幅度分布服从高斯分布，而它的功率谱密度有事均匀分布的，则称这个噪声为高斯白噪声。高斯白噪声二阶矩不相关，一阶矩为常数，是指先后信号在时间上的相关性。
高斯噪声包括热噪声和三里噪声。高斯噪声由它的事变平均值和两瞬时的协方差函数来确定，若噪声是平稳的，则平均值与时间无关，而协方差函数则变成仅和所考虑的两瞬时之差有关的相关函数，在意义上它等同于功率谱密度。高斯早生可以用大量独立的脉冲产生，从而在任何有限时间间隔内，这些脉冲中的每一个买充值与所有脉冲值得总和相比都可忽略不计。
高斯噪声是指它的概率密度函数服从高斯分布（即正态分布）的一类噪声。
def gauss_noise(image, mean, var=0.001):"""给图像添加高斯噪声:param image: 图像:param mean: 均值:param var: 方差:return:"""image = np.array(image / 255, dtype=np.float)noise = np.random.normal(mean, var ** 0.5, image.shape)out = image + noiseif out.min() < 0:low_clip = -1.else:low_clip = 0.out = np.clip(out, low_clip, 1.0)out = np.uint8(out * 255)return out
5.2.4 奇异值分解去噪
奇异值（Singular Value）往往对应着矩阵中的隐含的重要信息，且重要性与奇异值大小呈正相关。
一般来说，较少的奇异值就可以表达一个矩阵很重要的信息，所以我们可以舍掉一部分奇异值来实现压缩。
在图像处理中，奇异值小的部分往往代表着噪声，因此可以借助SVD算法来实现去噪。
def svd_denoise(img, radio=0.1):u, sigma, vt = np.linalg.svd(img)h, w = img.shape[:2]h_new = int(h * radio) # 取前10%的奇异值重构图像sigma_new = np.diag(sigma[:h_new], 0) # 用奇异值生成对角阵u_new = np.zeros((h, h_new), np.float)u_new[:, :] = u[:, :h_new]vt_new = np.zeros((h_new, w), np.float)vt_new[:, :] = vt[:h_new, :]return u_new.dot(sigma_new).dot(vt_new)
【完整程序】
import cv2
import numpy as np
import random
import matplotlib.pyplot as pltdef sp_noise(image, prob):"""给图像加椒盐噪声:param image:图像:param prob:噪声比例:return:"""output = np.zeros(image.shape, np.uint8)thres = 1 - probfor i in range(image.shape[0]):for j in range(image.shape[1]):rdn = random.random()if rdn < prob:output[i][j] = 0elif rdn > thres:output[i][j] = 255else:output[i][j] = image[i][j]return outputdef gauss_noise(image, mean, var=0.001):"""给图像添加高斯噪声:param image: 图像:param mean: 均值:param var: 方差:return:"""image = np.array(image / 255, dtype=np.float)noise = np.random.normal(mean, var ** 0.5, image.shape)out = image + noiseif out.min() < 0:low_clip = -1.else:low_clip = 0.out = np.clip(out, low_clip, 1.0)out = np.uint8(out * 255)return outdef svd_denoise(img, radio=0.1):u, sigma, vt = np.linalg.svd(img)h, w = img.shape[:2]h_new = int(h * radio) # 取前10%的奇异值重构图像sigma_new = np.diag(sigma[:h_new], 0) # 用奇异值生成对角阵u_new = np.zeros((h, h_new), np.float)u_new[:, :] = u[:, :h_new]vt_new = np.zeros((h_new, w), np.float)vt_new[:, :] = vt[:h_new, :]return u_new.dot(sigma_new).dot(vt_new)def main():img = cv2.imread('./dataset/images/gaoyuanyuan.jpg', cv2.IMREAD_COLOR)img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)print(img.shape)# 加噪声out1 = sp_noise(img, 0.05)out2 = gauss_noise(img, 0, 0.001)# 去噪声out1_denoise = svd_denoise(out1)out2_denoise = svd_denoise(out2)# 显示图像titles = [['Original', 'Add Salt and Pepper noise', 'Denoise'],['Original', 'Add Gaussian noise', 'Denoise']]images = [[img, out1, out1_denoise],[img, out2, out2_denoise]]plt.figure()plt.subplot(2, 3, 1)plt.imshow(images[0][0], 'gray')plt.title(titles[0][0])plt.xticks([])plt.yticks([])plt.subplot(2, 3, 2)plt.imshow(images[0][1], 'gray')plt.title(titles[0][1])plt.xticks([])plt.yticks([])plt.subplot(2, 3, 3)plt.imshow(images[0][2], 'gray')plt.title(titles[0][2])plt.xticks([])plt.yticks([])plt.subplot(2, 3, 4)plt.imshow(images[1][0], 'gray')plt.title(titles[1][0])plt.xticks([])plt.yticks([])plt.subplot(2, 3, 5)plt.imshow(images[1][1], 'gray')plt.title(titles[1][1])plt.xticks([])plt.yticks([])plt.subplot(2, 3, 6)plt.imshow(images[1][2], 'gray')plt.title(titles[1][2])plt.xticks([])plt.yticks([])plt.show()if __name__ == '__main__':main()
参考文献
- 李航《统计学习方法》第二版
- https://blog.csdn.net/C_chuxin/article/details/84898942
- https://zhuanlan.zhihu.com/p/110020411
- https://blog.csdn.net/index20001/article/details/73501632
- https://www.cnblogs.com/LeftNotEasy/archive/2011/01/19/svd-and-applications.html
- https://blog.csdn.net/Khan__Liu/article/details/54581075
- https://blog.csdn.net/qq_38395705/article/details/106311905
机器学习理论《统计学习方法》学习笔记：奇异值分解（SVD）相关推荐
- 统计学习方法 学习笔记（1）统计学习方法及监督学习理论
统计学习方法及监督学习理论 1.1.统计学习 1.1.1.统计学习的特点 1.1.2.统计学习的对象 1.1.3.统计学习的目的 1.1.4.统计学习的方法 1.1.5.统计学习的研究 1.1.6.统 ...
- 逻辑斯蒂回归_逻辑斯蒂回归详细解析 | 统计学习方法学习笔记 | 数据分析 | 机器学习...
本文包括: 重要概念 逻辑斯蒂回归和线性回归 二项逻辑斯谛回归模型 逻辑斯蒂回顾与几率 模型参数估计 多项逻辑斯谛回归 其它有关数据分析,机器学习的文章及社群 1.重要概念: 在正式介绍逻辑斯蒂回归模 ...
- 统计学习方法 学习笔记（十）：决策树
这一个学习笔记将要了解决策树,在研一上机器学习这门课的时候,老师在讲到这一节的时候,举了一个例子我现在还能记得:你们坐在这里上课,就像这个决策树一样,在你人生中的每一个重要结点,你都做出了选择,经过多 ...
- 统计学习方法 学习笔记（五）：支持向量机（下）
通过支持向量机(上)和支持向量机(中)的介绍,对支持向量机应该有点感性的认识啦!在这个学习笔记中,来继续探寻带核函数的支持向量机(解决如下图所示的问题) 对解线性分类问题,线性分类支持向量机是一种非常 ...
- 统计学习方法读书笔记（六）-逻辑斯蒂回归与最大熵模型（迭代尺度法（IIS））
全部笔记的汇总贴:统计学习方法读书笔记汇总贴 逻辑斯谛回归 (logistic regression )是统计学习中的经典分类方法.最大熵是概率模型学习的一个准则,将其推广到分类问题得到最大熵模型(m ...
- 统计学习方法读书笔记（九）-EM算法及其推广
全部笔记的汇总贴:统计学习方法读书笔记汇总贴 EM算法用于含有隐变量(hidden variable)的概率模型参数的极大似然估计,或极大后验概率估计.EM算法的每次迭代由两步组成:E步,求期望(ex ...
- 基于神经网络的机器阅读理解综述学习笔记
基于神经网络的机器阅读理解综述学习笔记 一.机器阅读理解的任务定义 1.问题描述 机器阅读理解任务可以形式化成一个有监督的学习问题:给出三元组形式的训练数据(C,Q,A),其中,C 表示段落,Q 表示 ...
- 统计学习方法-李航-第一章：统计学习方法概论-笔记1
文章目录 0 机器学习分类 0.1 监督学习 0.2 无监督学习 0.3 半监督学习 0.4 强化学习 1 统计学习方法概论 1.1 监督学习的步骤 1.2 统计学习三要素 1.3 模型评估 1.4 ...
- 李航《统计学习方法》笔记
虽然书名是统计学习,但是却是机器学习领域中和重要的一本参考书.当前的机器学习中机器指计算机,但是所运用的方法和知识是基于数据(对象)的统计和概率知识,建立一个模型,从而对未来的数据进行预测和分析(目的 ...
- 机器学习：《统计学习方法》笔记（一）—— 隐马尔可夫模型
参考:<统计学习方法>--李航:隐马尔可夫模型--码农场 摘要 介绍隐马尔可夫模型的基本概念.概率计算.学习方法.预测方法等内容. 正文 1. 基本概念 隐马尔可夫模型是关于时序的模型,描 ...
最新文章
- python自学教程读书导图-自学Python第一天：起点读书自动领取经验值(附思路讲解）...
- python怎么识别拼音-python获取一组汉字拼音首字母的方法
- PyTorch 1.10正式版上线了！附相关资源
- Apple 远程推送APNS 服务
- fanuc系统ug后处理_UG新版后置post configurator后处理配置器之备刀(预选刀)换刀不输出T问题处理方法...
- Nginx location执行顺序和匹配规则
- PLSQL 日期格式修改
- iOS开发UI篇—UIScrollView控件实现图片缩放功能
- [转]rails常用验证方法
- 在Linux中修复U盘
- （8）USB协议 —— 高速模式握手过程
- 智能配电房综合监控系统的探讨
- 计算机网络的定义及答案,计算机网络习题库
- 奇迹之剑萌新晋升大神辅助攻略 奇迹之剑游戏脚本挂机工具介绍
- c语言表达ch是大写英文字母,如何用C语言输出26个英文字母和其ascii码的对照表...
- cs224w（图机器学习）2021冬季课程学习笔记14 Reasoning over Knowledge Graphs
- 【学习OpenCV】使用OpenCV播放AVI视频
- 货币汇率换算器隐私协议
- 晒弟弟考取的教资证写的朋友圈文案
- CentOS7.3下Zabbix3.5之微信报警配置