奇異值分解（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\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,?,{\sigma }_p)$$
$${\sigma }_1\ge {\sigma }_2\ge ?\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,?,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 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\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 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]] 驗證完全奇異值分解 True

5 奇異值分解的應用

圖像壓縮（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=?\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}?$$
當然，最簡單的方法就是給這個矩陣直接發過去，這是一個5x5的矩陣，你至少需要發送25個數字。

但是我們可以把這個矩陣分解為兩個矩陣的乘積，這樣只需要發送10個數字。

$$A=?\begin{array}{c}1\\ 1\\ 1\\ 1\\ 1\end{array}?\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 = 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）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。