摘要

PCA的实现一般有两种，一种是用特征值分解去实现的，一种是用奇异值分解去实现的。在上篇文章中便是基于特征值分解的一种解释。特征值和奇异值在大部分人的印象中，往往是停留在纯粹的数学计算中。而且线性代数或者矩阵论里面，也很少讲任何跟特征值与奇异值有关的应用背景。

奇异值分解是一个有着很明显的物理意义的一种方法，它可以将一个比较复杂的矩阵用更小更简单的几个子矩阵的相乘来表示，这些小矩阵描述的是矩阵的重要的特性。就像是描述一个人一样，给别人描述说这个人长得浓眉大眼，方脸，络腮胡，而且带个黑框的眼镜，这样寥寥的几个特征，就让别人脑海里面就有一个较为清楚的认识，实际上，人脸上的特征是有着无数种的，之所以能这么描述，是因为人天生就有着非常好的抽取重要特征的能力，让机器学会抽取重要的特征，SVD是一个重要的方法。

奇异值分解（SVD）是一种矩阵因子分解方法，是线性代数的概念，但在统计学习中被广泛使用，成为其重要工具，其中主成分分析、潜在语义分析都用到奇异值分解。

任意一个$m\times n$矩阵，都可以表示为三个矩阵的乘积（因子分解）形式，分别是m阶正交矩阵、由降序排列的非负的对角线元素组成的$m\times n$矩形对角矩阵、n阶正交矩阵，称为该矩阵的奇异值分解。矩阵的奇异值分解一定存在，但不唯一。奇异值分解可以看作是矩阵数据压缩的一种方法，即用因子分解的方式近似地表示原始矩阵，这种近似是在平方损失意义下的最优近似。

1 奇异值分解的定义与定理

1.1 奇异值分解的定义

矩阵的奇异值分解是指将一个非零的$m\times n$实矩阵$A$，$A\in R^{m\times n}$，表示为以下三个实矩阵乘积形式的运算，即进行矩阵的因子分解：$$A=U\Sigma V^T$$，其中U是m阶正交矩阵，V是n阶正交矩阵，$\Sigma$是由降序排列的非负的对角线元素组成的$m\times n$矩形对角矩阵，满足
$$U{U}^T=I$$
$$V{V}^T=I$$
$$\Sigma =diag({\sigma }_1,{\sigma }_2,\cdots ,{\sigma }_p)$$
$${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_p\ge 0$$
$$p=min(m,n)$$
$A=U\Sigma V^T$称为矩阵A的奇异值分解，${\sigma }_i$称为矩阵A的奇异值，U的列向量称为A的左奇异向量，V的列向量称为右奇异向量。

1.2 奇异值分解的基本定理

若$A$为一$m\times n$实矩阵，$A\in R^{m\times n}$，则$A$的奇异值分解存在$$A=U\Sigma V^T$$，其中U是m阶正交矩阵，V是n阶正交矩阵，$\Sigma$是$m\times n$矩形对角矩阵，其对角线元素非负，且按降序排列。

任意给定一个实矩阵，奇异值分解矩阵是否一定存在呢？由奇异值分解的基本定理可知，答案是肯定的。

1.3 奇异值分解的几何解释

从线性变换的角度理解奇异值分解，$m\times n$矩阵A表示从n维空间$R^n$到m维空间$R^m$的一个线性变换，$$T:x\to Ax$$，$x\in R^n,Ax\in R^m$，x和Ax分别是各自空间的向量。

线性变换可以分解为三个简单的变换：一个坐标系的旋转或反射变换、一个坐标轴的缩放变换、另一个坐标系的旋转或反射变换。奇异值定理保证这种分解一定存在。

2 紧奇异值分解和截断奇异值分解

奇异值分解定义给出的奇异值分解$$A=U\Sigma V^T$$，又称为矩阵的完全奇异值分解。实际常用的是奇异值分解的紧凑形式和截断形式。紧奇异值分解是与原始矩阵等秩的奇异值分解，截断奇异值分解是比原始矩阵低秩的奇异值分解。

2.1 紧奇异值分解

设有$m\times n$实矩阵$A$，其秩为$rank(A)=r,r\le min(m,n)$，则称$U_r{\Sigma }_r{V}_r^T$为A的紧奇异值分解，即$$A=U_r{\Sigma }_r{V}_r^T$$
，其中$U_r$是$m\times n$矩阵，$V_r$是$n\times r$矩阵，${\Sigma }_r$是$r$阶对角矩阵；矩阵$U_r$由完全奇异值分解中U的前r列、矩阵$V_r$由V的前r列、矩阵${\Sigma }_r$由$\Sigma$的前r个对角线元素得到。紧奇异值分解的对角矩阵${\Sigma }_r$的秩与原始矩阵A的秩相等。

2.2 截断奇异值分解

在矩阵的奇异值分解中，只取最大的k个奇异值（$k<r,r$为矩阵的秩）对应的部分，就得到矩阵的截断奇异值分解。实际应用中提到矩阵的奇异值分解时，通常指截断奇异值分解。

设A是$m\times n$实矩阵，其秩为$rank(A)=r,0<k<r$，则称$U_k{\Sigma }_k{V}_k^T$为矩阵A的截断奇异值分解$$A\approx U_k{\Sigma }_k{V}_k^T$$。其中$U_k$是$m\times k$矩阵，$V_k$是$n\times k$矩阵，${\Sigma }_k$是$k$阶对角矩阵；矩阵$U_k$由完全奇异值分解中U的前k列，矩阵$V_k$由V的前k列，矩阵${\Sigma }_k$由$\Sigma$的前k个对角线元素得到。对角矩阵${\Sigma }_k$的秩比原始矩阵A的秩低。

3 奇异值分解的计算

矩阵A的奇异值分解可以通过求对称矩阵$A^{T}A$的特征值和特征向量得到。$A^{T}A$的特征向量构成正交矩阵V的列；$A^{T}A$的特征值${\lambda }_j$的平方根为奇异值${\sigma }_i$，即：$${\sigma }_j={\sqrt{\lambda }}_j,j=1,2,\cdots ,n$$对其由大到小排列作为对角线元素，构成对角矩阵$\Sigma$；求正奇异值对应的左奇异向量，再求扩充的$A^T$的标准正交基，构成正交矩阵U的列。从而得到A的奇异值分解$$A=U\Sigma V^T$$

4 特征值分解和奇异值分解

• 只有方阵才能进行特征值分解。

4.1 特征值分解（EIG）

4.1.1 特征值分解的定义

如果说一个向量$v$是方阵$A$的特征向量，即可以表示为下面形式$$Av=\lambda v$$，此时$\lambda$称为特征向量$v$对应的特征值，一个矩阵的一组特征向量是一组正交向量，

特征值分解是将一个矩阵分解为下面形式：$$A=Q\Sigma 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 np

a = 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\Sigma V^T$$

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 np

a = 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 np

a = 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]]
验证完全奇异值分解 True

5 奇异值分解的应用

1. 图像压缩（image compression）：较少的奇异值就可以表达出图像中大部分信息，舍弃掉一部分奇异值来实现压缩。

2. 图像降噪（image denoise）：噪声一般存在于图像高频部分，也表现在奇异值小的部分，故可以借助SVD实现去噪。

3. 音频滤波（filtering）：Andrew Ng的机器学习课程上有个svd将混杂声音分离的例子，其实和噪声滤波类似。

4. 求任意矩阵的伪逆（pseudo-inverse）：由于奇异矩阵或非方阵矩阵不可求逆，在特殊情况下需要广义求逆时可用svd方法。

5. 模式识别（pattern recognition）：特征为矩阵，数据量较大时，可以用svd提取主要的成分。

6. 潜在语义索引（Latent Semantic Indexing）：NLP中，文本分类的关键是计算相关性，这里关联矩阵A=USV’，分解的三个矩阵有很清楚的物理含义，可以同时得到每类文章和每类关键词的相关性。

5.1 图像压缩

5.1.1 奇异值分解压缩的原理

先看一个简单的例子，如果你想要在网络上给别人发送一段数据，数据的内容为
$$A=\left(\begin{array}{ccccc}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{array}\right)$$
当然，最简单的方法就是给这个矩阵直接发过去，这是一个5x5的矩阵，你至少需要发送25个数字。

但是我们可以把这个矩阵分解为两个矩阵的乘积，这样只需要发送10个数字。

$$A=\left(\begin{array}{c}1\\ 1\\ 1\\ 1\\ 1\end{array}\right)\left(\begin{array}{ccccc}1& 1& 1& 1& 1\end{array}\right)$$
图像也可以被视为矩阵，图像的每一个点都是由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 = 0
    k = 0

    while 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 += 1
    a[a < 0] = 0
    a[a > 255] = 255
    return 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 = 0
    k = 0

    while 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 += 1
    a[a < 0] = 0
    a[a > 255] = 255
    return 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) / 10
        plt.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 Image


def 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_img


def 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 - prob

    for i in range(image.shape[0]):
        for j in range(image.shape[1]):
            rdn = random.random()
            if rdn < prob:
                output[i][j] = 0
            elif rdn > thres:
                output[i][j] = 255
            else:
                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 + noise
    if 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 plt


def sp_noise(image, prob):
    """
    给图像加椒盐噪声
    :param image:图像
    :param prob:噪声比例
    :return:
    """
    output = np.zeros(image.shape, np.uint8)
    thres = 1 - prob

    for i in range(image.shape[0]):
        for j in range(image.shape[1]):
            rdn = random.random()
            if rdn < prob:
                output[i][j] = 0
            elif rdn > thres:
                output[i][j] = 255
            else:
                output[i][j] = image[i][j]
    return output


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 + noise
    if 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


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)


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

参考文献

1. 李航《统计学习方法》第二版
2. https://blog.csdn.net/C_chuxin/article/details/84898942
3. https://zhuanlan.zhihu.com/p/110020411
4. https://blog.csdn.net/index20001/article/details/73501632
5. https://www.cnblogs.com/LeftNotEasy/archive/2011/01/19/svd-and-applications.html
6. https://blog.csdn.net/Khan__Liu/article/details/54581075
7. https://blog.csdn.net/qq_38395705/article/details/106311905