圖像處理——主成分分析

時間 2019-11-11

標籤圖像處理成分分析简体版

原文原文鏈接

1 引言html

　　1.1 維度災難python

　　　　分類爲例：如最近鄰分類方法（基本思想：以最近的格子投票分類）git

　　　　問題：當數據維度增大，分類空間爆炸增加。如圖1所示，github

　　　　　　　　　　　　　　　　　　　　　　圖1 維度增長示意圖shell

　　1.2 解決方法數組

　　　　緩解維度遭難的一個重用途徑是降維。降維是經過某種數學變換，將原始高維屬性空間轉換爲一個低維「子空間」，在這個子空間中樣本機器學習

　　密度大幅度提升，距離計算也變得更爲容易。函數

　　1.3 降維的可行性學習

　　　　數據樣本雖然是高維的，但與咱們關心的也許僅是某個低維分佈，於是能夠對數據進行有效的降維。測試

　　1.4 主成分分析

　　1.4.1 概念

　　　　主成分分析(Principal Component Analysys,簡稱PCV)，是將具備相關性的高維指標,經過線性變換，化爲少數相互獨立的低維綜合指標的

　　一種多元統計方法。又稱爲主份量分析。使用PCV降維時應該知足：

　　　　　　　　　　　　　　　（1）使降維後樣本的方差儘量大。

　　　　　　　　　　　　　　　（2）使降維後數據的均方偏差儘量小。

　　1.4.2 基本計算步驟

　　　　（1）計算給定樣本{x_n},n=1,2,...,N的均值mean_x和方差S；

　　　　（2）計算S的特徵向量與特徵值，X = UΛU^T；

　　　　（3）將特徵值從小到大排列，前M個特徵值λ1，λ1，...,λM所對應的特徵向量u₁,u₂,...,u_M構成投影矩陣。

2 矩陣特徵值及特徵向量求解

　　2.1 定義

　　　　A爲n階矩陣，若數λ和n維非0列向量x知足Ax=λx，那麼數λ稱爲A的特徵值，x稱爲A的對應於特徵值λ的特徵向量。

　　　　式Ax=λx也可寫成( A-λE)x=0，而且|λE-A|叫作A 的特徵多項式。

　　　　|λE-A| = 0，稱爲A的特徵方程（一個齊次線性方程組），求解特徵值的過程其實就是求解特徵方程的解。

　　2.2 特徵值分解

　　　　特徵值分解是將一個矩陣分解成下面的形式：

　　　　其中Q是這個矩陣A的特徵向量組成的矩陣。Σ = diag(λ₁, λ₂, ..., λ_n)是一個對角陣，每個對角線上的元素就是一個特徵值。

3 奇異值分解（SVD）

　　特徵值分解只適用於方陣。而在現實的世界中，咱們看到的大部分矩陣都不是方陣，咱們怎樣才能像描述特徵值同樣描述這樣通常矩陣的重要特徵呢？奇異值分解（Singular Value Decomposition, 簡稱SVD）可用來解決這個問題。SVD比特徵值分解的使用範圍更廣，是一個適用於任一矩陣分解的方法。

　　假設X是一個m * n的矩陣，那麼獲得的U是一個m * m的方陣（裏面的向量是正交的，U裏面的向量稱爲左奇異向量）,

　　　　　　　　　　　　　　　　　　　Σ是一個m * n的實數對角矩陣（對角線之外的元素都是0，對角線上的元素稱爲奇異值）,

　　　　　　　　　　　　　　　　　　　V^T(V的轉置)是一個n * n的矩陣，裏面的向量也是正交的，V裏面的向量稱爲右奇異向量，

　　從圖片來反映幾個相乘的矩陣的大小可得下面的圖片,

　　將一個矩陣X的轉置 X^T * X，將會獲得 X^TX 是一個方陣，咱們用這個方陣求特徵值能夠獲得：

　(A^TA)v_i = λ_iv_i

　這裏獲得的V，就是咱們上面的右奇異向量。這裏的σ_i 就是就是上面說的奇異值λ_i，u_i就是上面說的左奇異向量。常見的作法是將奇異值由大而小排列。如此Σ便能由M惟一肯定了。

　　奇異值σ跟特徵值相似，在矩陣Λ中也是從大到小排列，並且σ的減小特別的快，在不少狀況下，前10%甚至1%的奇異值的和就佔了所有的奇異值之和的99%以上了。也就是說，咱們也能夠用前r大的奇異值來近似描述矩陣，這裏定義一下部分奇異值分解：

r是一個遠小於m、n的數，這樣矩陣的乘法看起來像是下面的樣子：

　　右邊的三個矩陣相乘的結果將會是一個接近於X的矩陣，在這兒，r越接近於n，則相乘的結果越接近於A。而這三個矩陣的面積之和（在存儲觀點來講，矩陣面積越小，存儲量就越小）要遠遠小於原始的矩陣A，咱們若是想要壓縮空間來表示原矩陣A，咱們存下這裏的三個矩陣：U、Λ、V就行了。

4 特徵分解與奇異值分解

　　4.1 低秩近似（Low-rank Approximation）

　　4.2 關聯

5 pca(matrix)函數代碼

 1 def pca(X):
 2     """    主成分分析(PCA)
 3         輸入: X, （每行存儲訓練數據（一張圖像）的扁平數組）
 4         返回: 投影矩陣(V,重要維度優先)、方差(S)和均值(mean_X)
 5     """
 6     
 7     # 得到維度：num_data,dim = 圖片數量 *  圖片大小       ------>圖片大小 = 長*寬 
 8     num_data,dim = X.shape
 9     
10     # 中心數據
11     mean_X = X.mean(axis=0)    # 每列求平均值 獲得 1 * dim
12     X = X - mean_X             # X.shape = num_data * dim
13     
14     if dim > num_data:
15         # PCA - 緊湊的使用技巧
16         M = dot(X,X.T) 　　　　　# 協方差矩陣
17         e,EV = linalg.eigh(M) 　# e,EV = M的特徵值，特徵向量
18         tmp = dot(X.T,EV).T 　　# 使用緊湊技巧？
19         V = tmp[::-1] 　　　　　　# 爲了獲得最後的特徵向量
20         S = sqrt(e)[::-1]　　　  # 反向，由於特徵值是遞增的
21         for i in range(V.shape[1]):   # V.shape[1] = ？？？
22             V[:,i] /= S
23     else:
24         # PCA - 使用奇異值分解(SVD)
25         U,S,V = linalg.svd(X)
26         V = V[:num_data] 　　　　　# 只有返回前num_data個纔有意義       可是我認爲此處不該該爲V[:num_data]而應該爲 V[:dim]  ------>由於num_data >= dim
27     
28     # 返回n the projection matrix, the variance and the mean
29     return V,S,mean_X

6 簡單實例說明PCV實現的過程

6.1 列數大於行數時

 1 # 初始化矩陣X  （列數大於行數時）
 2 X =  [[1 0 1 0 0 0]
 3         [0 1 0 0 0 0]
 4         [1 1 0 0 0 0]
 5         [1 0 0 1 1 0]
 6         [0 0 0 1 0 1]]
 7 #計算每列的平均值
 8 mean_X =X.mean(axis=0)= [0.6 0.4 0.2 0.4 0.2 0.2]
 9 # 計算中心數據
10 X =X - mean_X =  [[ 0.4 -0.4  0.8 -0.4 -0.2 -0.2]
11  [-0.6  0.6 -0.2 -0.4 -0.2 -0.2]
12  [ 0.4  0.6 -0.2 -0.4 -0.2 -0.2]
13  [ 0.4 -0.4 -0.2  0.6  0.8 -0.2]
14  [-0.6 -0.4 -0.2  0.6 -0.2  0.8]]
15 # 計算協方差矩陣
16 M =dot(X, X.T)=[[ 1.20000000e+00 -4.00000000e-01  5.55111512e-17 -2.00000000e-01
17   -6.00000000e-01]
18  [-4.00000000e-01  1.00000000e+00  4.00000000e-01 -8.00000000e-01
19   -2.00000000e-01]
20  [ 5.55111512e-17  4.00000000e-01  8.00000000e-01 -4.00000000e-01
21   -8.00000000e-01]
22  [-2.00000000e-01 -8.00000000e-01 -4.00000000e-01  1.40000000e+00
23    0.00000000e+00]
24  [-6.00000000e-01 -2.00000000e-01 -8.00000000e-01  0.00000000e+00
25    1.60000000e+00]]
26 # 計算協方差的特徵值和特徵向量
27 e, EV = numpy.linalg.eigh(M) 
28 e =  [-6.04093951e-16  2.66097365e-01  1.17478886e+00  2.00000000e+00
29   2.55911378e+00]
30 EV =  [[ 4.47213595e-01 -1.38203855e-01  6.92423635e-01  5.00000000e-01
31   -2.26824169e-01]
32  [ 4.47213595e-01 -6.14411448e-01 -1.73966345e-01 -5.00000000e-01
33   -3.77139608e-01]
34  [ 4.47213595e-01  6.98980014e-01 -3.02870626e-01  4.38587603e-16
35   -4.68717745e-01]
36  [ 4.47213595e-01 -2.11286151e-01 -5.40988322e-01  5.00000000e-01
37    4.61183041e-01]
38  [ 4.47213595e-01  2.64921441e-01  3.25401658e-01 -5.00000000e-01
39    6.11498480e-01]]
40 #行之間反轉
41 tmp = dot(X.T, EV).T =[[-3.33066907e-16 -3.33066907e-16  4.16333634e-16 -3.33066907e-16
42   -3.88578059e-16  1.66533454e-16]
43  [ 3.49490007e-01  8.45685660e-02 -1.38203855e-01  5.36352895e-02
44   -2.11286151e-01  2.64921441e-01]
45  [-1.51435313e-01 -4.76836971e-01  6.92423635e-01 -2.15586664e-01
46   -5.40988322e-01  3.25401658e-01]
47  [ 1.00000000e+00 -5.00000000e-01  5.00000000e-01 -7.21644966e-16
48    5.00000000e-01 -5.00000000e-01]
49  [-2.34358872e-01 -8.45857353e-01 -2.26824169e-01  1.07268152e+00
50    4.61183041e-01  6.11498480e-01]]
51 #默認返回的是按特徵值升序，上下調轉一下
52 V =tmp[::-1] = [[-2.34358872e-01 -8.45857353e-01 -2.26824169e-01  1.07268152e+00
53    4.61183041e-01  6.11498480e-01]
54  [ 1.00000000e+00 -5.00000000e-01  5.00000000e-01 -7.21644966e-16
55    5.00000000e-01 -5.00000000e-01]
56  [-1.51435313e-01 -4.76836971e-01  6.92423635e-01 -2.15586664e-01
57   -5.40988322e-01  3.25401658e-01]
58  [ 3.49490007e-01  8.45685660e-02 -1.38203855e-01  5.36352895e-02
59   -2.11286151e-01  2.64921441e-01]
60  [-3.33066907e-16 -3.33066907e-16  4.16333634e-16 -3.33066907e-16
61   -3.88578059e-16  1.66533454e-16]]
62 #sigma = 特徵值開根號
63 #使特徵值降序排序
64 S =sqrt(e)[::-1]= [1.59972303 1.41421356 1.08387677 0.51584626        nan]
65 #特徵向量歸一化
66 for i in range(V.shape[1]):
67      V[:, i] /= S
68 -----------------------------------------------------------------------------     
69 V =  [[-1.46499655e-01 -5.28752375e-01 -1.41789650e-01  6.70542025e-01
70    2.88289305e-01  3.82252720e-01]
71  [ 7.07106781e-01 -3.53553391e-01  3.53553391e-01 -5.10280049e-16
72    3.53553391e-01 -3.53553391e-01]
73  [-1.39716356e-01 -4.39936516e-01  6.38839814e-01 -1.98903298e-01
74   -4.99123458e-01  3.00220160e-01]
75  [ 6.77508074e-01  1.63941415e-01 -2.67916753e-01  1.03975338e-01
76   -4.09591321e-01  5.13566659e-01]
77  [            nan             nan             nan             nan
78               nan             nan]]
79 S =  [1.59972303 1.41421356 1.08387677 0.51584626        nan]

6.2 列數不大於行數時

 1  1 #初始化數組X
 2  2 X =  [[ 1  3]
 3  3  [ 4  7]
 4  4  [ 3  2]
 5  5  [ 1 23]]
 6  6 num_data, dim = X.shape
 7  7 X.shape= (4, 2)
 8  8 U, S, V = linalg.svd(X)
 9  9 U.shape, S.shape, V.shape=(4, 4) (2,) (2, 2)
10 10 V = [[ 0.10462808  0.99451142]
11 11  [ 0.99451142 -0.10462808]]
12 12 V = V[:dim]                                        #V = V[:dim]  or V =V[:num_data]    ???   注意：此時num_data > dim，因此我認爲應該是V[:dim]
13 13 V[dim-1] =V[ 1 ]= [ 0.99451142 -0.10462808]
14 14 #計算結果
15 15 U = [[ 0.12635702  0.149642   -0.4064109  -0.89245244]
16 16  [ 0.30196809  0.71358512 -0.49810453  0.38923441]
17 17  [ 0.09422707  0.60994996  0.75324035 -0.22740114]
18 18  [ 0.94019702 -0.31042648  0.13910799  0.0177182 ]]
19 19 S =  [24.43997403  4.54836997]
20 20 V =  [[ 0.10462808  0.99451142]
21 21  [ 0.99451142 -0.10462808]]

7 環境配置

7.1 須要準備的安裝包

　　下載python(x,y)2.7.x: http://www.softpedia.com/get/Programming/Other-Programming-Files/Python-x-y.shtml

　　下載PVC包：https://github.com/willard-yuan/pcv-book-code

7.2 安裝Python(x, y)

　　在Windows下，推薦安裝Python(x,y)2.7.x。Python(x,y)2.7.x是一個庫安裝包，除了Python自身外，還包含了第三方庫，下面是安裝Python(x, y)時的界面：

　　　　選擇Full模式，自動安裝的包見下圖，

從上面第二幅圖能夠看出，pythonxy不只包含了SciPy、NumPy、PyLab、OpenCV、MatplotLib,還包含了機器學習庫scikits-learn。爲避免出現運行實例時出現的依賴問題，譯者建議將上面的庫所有選上，也就是選擇「full」(譯者也是用的所有安裝的方式進行後面的實驗的)。安裝完成後，爲驗證安裝是否正確，能夠在Python shell裏確認一下OpenCV是否已安裝來進行驗證，在Python Shell裏輸入下面命令：

1 from cv2 import __version__
2 __version__

自動安裝的pyopencv有問題。須要你把opencv卸載後從新安裝就能夠了。從新安裝完成後，再次輸入上面命令，就能夠看到OpenCV的版本信息，則說明python(x,y)已安裝正確。

　　7.2 安裝PCV庫

　　7.2.1 安裝twine

　　　　cd 到對應python的scripts目錄中，執行如下命令

1 pip install twine　　
或者
2 python3 -m pip install twine

　　7.2.2 PVC的安裝

　　　　PCV庫是一個第三方庫。安裝PVC以前必須安裝twine，不然安裝會出錯。假設你已從上面網址上下載了中譯版源碼，從Windows cmd終端　　

　　進入PCV所在目錄：

1 cd PCV
2 python setup.py install

　　　　運行上面命令，便可完成PCV庫的安裝。爲了驗證PCV庫是否安裝成功，在運行上面命令後，能夠打開Python自帶的Shell，在Shell輸入：

1 import PCV

　　　　若是未報錯，則代表你已成功安裝了該PCV庫。

8 圖像處理案案例

8.1 圖像處理測試代碼

 1 # -*- coding: utf-8 -*-
 2 import pickle
 3 from PIL import Image
 4 from numpy import *
 5 from pylab import *
 6 from PCV.tools import imtools, pca
 7 
 8 # 獲取圖像列表和他們的尺寸
 9 imlist = imtools.get_imlist('../data/fontimages/a_thumbs')  # fontimages.zip is part of the book data set
10 im = array(Image.open(imlist[0]))  # open one image to get the size
11 m, n = im.shape[:2]  # get the size of the images
12 imnbr = len(imlist)  # get the number of images
13 print "The number of images is %d" % imnbr
14 
15 # Create matrix to store all flattened images
16 immatrix = array([array(Image.open(imname)).flatten() for imname in imlist], 'f')
17 
18 # PCA降維
19 V, S, immean = pca.pca(immatrix)
20 
21 # 保存均值和主成分
22 f = open('../data/fontimages/font_pca_modes.pkl', 'wb')
23 pickle.dump(immean, f)
24 pickle.dump(V, f)
25 f.close()
26 
27 # Show the images (mean and 7 first modes)
28 # This gives figure 1-8 (p15) in the book.
29 figure()
30 gray()
31 subplot(2, 4, 1)
32 axis('off')
33 imshow(immean.reshape(m, n))
34 for i in range(7):
35     subplot(2, 4, i+2)
36     imshow(V[i].reshape(m, n))
37     axis('off')
38 show()