2018/12/14

Dimensionality Reduction with PCA (主成份分析於縮減維度)

15 Dimensionality Reduction with PCA (主成份分析於縮減維度)

所謂的維度的魔咒 (Curse of dimensionality) 指資料集的維度增加會帶來災難,因當維度增加所需的資料集的實例也必需增加,這會使計算的複雜度快速成長,有時會變成無法處理 (infeasible)。因而 有些時候降低資料集的維度是一個很重要的課題。Principle component analysis (PCA, 主成份分析) 是一種常見用於縮減維度 的方法,要了解PCA 必需先了解一些基本的多變量統計的觀念。

⬛ 多變量統計
變量(variable)或稱變數是指在資料裡的一個特徵,而多變量是指要考慮很多的特徵,統計是指要求出這些變量彼此間存在的統計特性。

►平均值 (Mean)
平均值是種最基本的統計量,就是將所有數相加求平均, $$\bar{x}=\frac{1}{n}\sum_{i=1}^{n}x_i$$ 上式表示有$n$相加,所以加後再除$n$。 Numpy 的 mean()可用來求平均值,如下
import numpy as np
v = np.array([1,2,3,4,5,6])
#  v  = [1, 2, 3, 4, 5, 6]
print(v)
result = np.mean(v)
print(result)
[1, 2, 3, 4, 5, 6]
3.5
註解掉的第3行是Python list 形式,和一維的 np.array 通用,所以可取代第2行。若是輸入矩陣形式,如
v = np.array([[1,2],[3,4],[5,6]])
print(v)
result0 = v.mean(0)
result1 = v.mean(1)
print(result0)
print(result1)
[[1 2]
 [3 4]
 [5 6]]
[3. 4.]
[1.5 3.5 5.5]
第3行的 mean(0)是指對每一行所有元素求平均,結果的行數不變,列數只剩1,若求每一列所有元素平均就要用 mean(1),結果是列數不變,行數只剩1。記法是0,1和平時用以index 的用法相反。若都不給參數是求矩陣內的所有元素的平均。

► 變異數 (Variance)和標準偏差(Standard deviations)
在機率理論裡變異數是求$(x-\mu)^2$的期望值,$\mu$是平均值,這稱為母體變異數 (Population variance),但在統計分析無法知道「真正的」機率分佈,只有取樣的資料,因而定義變異數如下 $$\begin{equation} \sigma^2=\frac{1}{n-1}\sum_{i=1}^{n}(x_i-\bar{x})^2 \end{equation}$$ 其中分母的$(n-1)$是為了要校正偏差(bias)。變異數常常用Var() 來表示。在numpy裡有個var()函數可求變異數,其中有個參數是 ddof (Delta Degrees of Freedom),預設是0,但因上式中有$(n-1)$,所以要設ddof=1,如
v = np.array([1,2,3,4,5,6])
result = np.var(v, ddof=1)
print(result)
3.5
若沒設ddof=1則輸出會是2.916。和求平均值相似,也可用矩陣形式輸入,如要求每一列的變異數
v = np.array([[1,2,3,4,5,6],[1,2,3,10,11,12]])
result = np.var(v, ddof=1, axis=1)
print(result)
[ 3.5 25.1]
第一列所有元素的變異數和先前求 的一樣, 第二列25.1。需注意,不管求行或列,numpy的輸出都是以列的形式表示。另外Eq.(1) 中的 $\sigma$ (沒有平方)就定義標準偏差。Numpy 也有個std()函數可,用法和var()類似。

► Covariance (共變異數) and Covariance Matrix (共變異矩陣)
在機率理論裡共變異數定義成兩隨機變數減去個自的平均值後相乘再取期望值,表示成 $$E[(X-\mu _x)(Y-\mu _y)]$$ 共變異數是用來了解兩機變數之間的關聯性,如其值為一個正數表示當一個變數增加則另外那個變數也是增加,如果是0 則表示彼此沒有關聯性(uncorrelated)。在統計可由取樣到的資料求共變異數,定義成 $$cov(X,Y)=\frac{1}{n-1}\sum_{i=1}^{n}(X_i-\bar{x})(Y_i-\bar{y})$$ 上式是指兩個變數的共異數,但實際應用中變數會遠大於2,所以常以陣形式表示,就稱為共變異矩陣。如果有$n$個變數,$x_i,i=0\cdots,n-1$, 將共變異矩陣表成 $$\begin{bmatrix} var(x_0)&cov(x_0,x_1) &\cdots & cov(x_0,x_{n-1})\\ cov(x_1,x_0) & var(x_1) & \cdots &cov(x_1,x_{n-1} \\ \vdots & \vdots & \vdots & \vdots\\ cov(x_{n-1},x_0)& cov(x_{n-1},x_1) & \cdots & var(x_{n-1}) \end{bmatrix}$$ 可以看出在矩陣中,若index $i \neq j$ 則該元素就是$cov(x_i,x_j)$,若$i=j$則是$x_i$ 的變異數。Numpy 的 cov() function可計算共變異數,此函數預設已含偏差修正,所以不用ddof=1的參數。如下列程式
x =np.array([1,2,3,4,5,6])
y = np.array([1,2,3,10,11,12])
Sigma = np.cov(x,y)
print(Sigma)
[[ 3.5  8.9]
 [ 8.9 25.1]]
從結果看出變異數的部份和之前求出的一樣,而共變異 8.9,因共變異數矩陣是對稱矩陣,以index $i,j$ 和$j,i$的值相同。

⬛  Eigenvalues (特徵值) ,Eigenvectors (特徵向量)
在矩陣的運算中,對某一方矩陣若 $$\begin{equation} \textbf{Ax}=\lambda \textbf{x} \end{equation}$$ 則$\lambda$是$\textbf{A}$的一個特徵,$\textbf{x}$是$\textbf{A}$的一個特徵向量,相關細節請參考線性代數課本。NumPy 有一 eig() 函數可用來求矩陣這兩個量。如
from numpy.linalg import eig
A = np.array([
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9]])
values, vectors = eig(A)
print(values)
print(vectors)
[ 1.61168440e+01 -1.11684397e+00 -9.75918483e-16]
[[-0.23197069 -0.78583024  0.40824829]
 [-0.52532209 -0.08675134 -0.81649658]
 [-0.8186735   0.61232756  0.40824829]]
若將矩陣和第一個特徵向量相乘要等於第1個特徵值和第1個特徵向量相乘,
B = A.dot(vectors[:, 0])
print(B)
C = vectors[:, 0] * values[0]
print(C)
[ -3.73863537 -8.46653421 -13.19443305]
[ -3.73863537 -8.46653421 -13.19443305]
結果顯示如同(2)所示,這兩者是相等的。
由線性代數中我們知可將矩陣$\textbf{A}$表示成 $$\begin{equation} \textbf{A}=\textbf{P}^{-1}\textbf{D}\textbf{P}\end{equation} $$ 下列程式驗證(3)式
from numpy.linalg import inv
Q = vectors
R = inv(Q)
L = np.diag(values)
B = Q.dot(L).dot(R)
print(B)
[[1. 2. 3.]
 [4. 5. 6.]
 [7. 8. 9.]]
顯示正確求得。

⬛  Principle component analysis (PCA, 主成份分析)
有兩種方式可以處理 PCA,一種就是求資料集的共變異矩陣,再用其特徵向量對資料進行轉換,另一種是較複雜的奇異值分解,我們先看前項。
import numpy as np
xdata = np.array([[0.9, 1], [2.4, 2.6],[1.2, 1.7], [0.5, 0.7],
             [0.3, 0.7],[1.8, 1.4],[0.5, 0.6],[0.3, 0.6],
                  [2.5, 2.6],[1.3, 1.1]])
X0 = xdata[:,0]
X1 = xdata[:,1]
x_cov = np.cov(xdata.T) 
print(x_cov)
w, v = np.linalg.eig(np.array(x_cov))
print(w)
print(v)
a= v[:,0]
a=a.reshape(2,1)
yyy = np.matmul(xdata,a)
print(yyy)
 [[0.68677778 0.60666667]
 [0.60666667 0.59777778]]
[1.25057433 0.03398123]
[[ 0.73251454 -0.68075138]
 [ 0.68075138  0.73251454]]
[[1.34001447]
 [3.5279885 ]
 [2.0362948 ]
 [0.84278324]
 [0.69628033]
 [2.27157811]
 [0.7747081 ]
 [0.62820519]
 [3.60123995]
 [1.70109543]]
程式中先輸入一筆$10\times 2$的資料,之後求 其共變異數的特徵值及特徵向量,特徵值為[1.25057433 0.03398123],因第二個值很小可捨棄,再用主特徵值對應的特徵向量進行轉換,轉換之後變成 $10\times 1$的資料,所以是將資料縮減維度的一種方法。

在scikit-learn 中有一個PCA class,其中 fit_transform() method 可將輸入資料降維,如
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.datasets import load_iris
data = load_iris()
y = data.target
X = data.data
pca = PCA(n_components=2)
reduced_X = pca.fit_transform(X)
red_x, red_y = [], []
blue_x, blue_y = [], []
green_x, green_y = [], []
for i in range(len(reduced_X)):
    if y[i] == 0:
        red_x.append(reduced_X[i][0])
        red_y.append(reduced_X[i][1])
    elif y[i] == 1:
        blue_x.append(reduced_X[i][0])
        blue_y.append(reduced_X[i][1])
    else:
        green_x.append(reduced_X[i][0])
        green_y.append(reduced_X[i][1])
plt.scatter(red_x, red_y, c='r', marker='x')
plt.scatter(blue_x, blue_y, c='b', marker='D')
plt.scatter(green_x, green_y, c='g', marker='.')
plt.show()
降維讓資料可視
程式中第3行載入scikit-learn 內建的 iris (鳶尾花)資料集,此集含有三種不同亞種的鳶尾花,因原本的特徵維度是4,將其降為2維,每一亞種用不同點色畫在圖上。此資料集進一步資訊如下
import numpy as np
from sklearn.datasets import load_iris
iris = load_iris()
iris.target[[10, 25, 50]]
list(iris.target_names)
# print(iris.DESCR)
print(iris.target)
print(iris.target_names)
X_iris, y_iris = iris.data, iris.target
print (X_iris.shape, y_iris.shape)
print (X_iris[49], y_iris[49])
 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 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 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 2 2 2 2 2 2 2 2 2 2 2
 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 2 2]
['setosa' 'versicolor' 'virginica']
(150, 4) (150,)
[5.  3.3 1.4 0.2] 0
共有150個實例,每一實例有4個 features和1個target。
三種鳶尾花

⬛  Singular Value Decomposition (奇異值分解 )
前述可將一方矩陣分解成三個矩陣相乘,奇異值分解也是用來將矩陣(不限定是方陣)分解成三個矩陣相乘, $$\begin{equation} \textbf{A}=\textbf{U}\boldsymbol{\Sigma} \textbf{V}^T \end{equation}$$ 其中$\textbf{A}$ 是一個實數 $n\times m$矩陣,$\textbf{U}$是$n\times n$的方陣,$\boldsymbol{\Sigma}$是$n\times m$的對角矩陣,$\textbf{V}$是 $m\times m$。$\boldsymbol{\Sigma}$裡的值是矩陣$\textbf{A}$的奇異值,$\textbf{U}$裡的行向量稱為$\textbf{A}$的左奇異向量,而$\textbf{V}$的行是右奇異向量,相關理論及如何分解的原請參考較進階的線代課本。

在Numpy中的svc() 函數可用來求此分解,傳回和(4)中右邊的三個矩陣有關的量。
from scipy.linalg import svd
A = np.array([
    [1, 2],
    [3, 4],
    [5, 6]])
U, s, V = svd(A)
print(U)
print(s)
print(V)
[[-0.2298477   0.88346102  0.40824829]
 [-0.52474482  0.24078249 -0.81649658]
 [-0.81964194 -0.40189603  0.40824829]]
[9.52551809 0.51430058]
[[-0.61962948 -0.78489445]
 [-0.78489445  0.61962948]]
因原矩陣是$3\times 2$,所以$\textbf{U}$ 是$3\times 3$,而$\boldsymbol{\Sigma}$ 應該是要$3\times 2$,但是只傳回其中非 0的兩個元素,而傳回的第三個量是$\textbf{V}^T$,是$2\times 2$, 若要原始的$\textbf{V}$需再轉置!我們可利用所傳回的矩陣再行組成原來的矩陣,如
Sigma = np.zeros((A.shape[0], A.shape[1]))
Sigma[:A.shape[1], :A.shape[1]] = np.diag(s)
y1  = np.matmul(Sigma,V)
y2  = np.matmul(U,y1)
y2
array([[1., 2.],
       [3., 4.],
       [5., 6.]])
顯示的結果是正確的。

► SVD 縮減維度
SVD 的一個很重要的應用就是用來縮減維度,首先給一個$3\times 10$的矩陣,
from scipy.linalg import svd
A = np.array([
[1,2,3,4,5,6,7,8,9,10],
[11,12,13,14,15,16,17,18,19,20],
[21,22,23,24,25,26,27,28,29,30]])
之後執行svd(),如
U, s, V = svd(A)
Sigma = np.zeros((A.shape[0], A.shape[1]))
Sigma[:A.shape[0], :A.shape[0]] = np.diag(s)
y1  = np.matmul(Sigma,V)
y2  = np.matmul(U,y1)
print(y2)
s
[[ 1.  2.  3.  4.  5.  6.  7.  8.  9. 10.]
 [11. 12. 13. 14. 15. 16. 17. 18. 19. 20.]
 [21. 22. 23. 24. 25. 26. 27. 28. 29. 30.]]
array([9.69657342e+01, 7.25578339e+00, 5.82172067e-15])
確認可轉回 ,但第三個奇異值相當小,幾可忽略,所以取前兩個來近似,
n_elements = 2
Sigma = Sigma[:, :n_elements]
V = V[:n_elements, :]
print(Sigma)
print(V)
B = np.matmul(U,np.matmul(Sigma,V))
print(B)
[[96.96573419  0.        ]
 [ 0.          7.25578339]
 [ 0.          0.        ]]
[[-0.24139304 -0.25728686 -0.27318068 -0.2890745  -0.30496832 -0.32086214
  -0.33675595 -0.35264977 -0.36854359 -0.38443741]
 [-0.53589546 -0.42695236 -0.31800926 -0.20906617 -0.10012307  0.00882003
   0.11776313  0.22670623  0.33564933  0.44459242]]
[[ 1.  2.  3.  4.  5.  6.  7.  8.  9. 10.]
 [11. 12. 13. 14. 15. 16. 17. 18. 19. 20.]
 [21. 22. 23. 24. 25. 26. 27. 28. 29. 30.]]
可見維後再轉回的矩陣是完全一樣。

► scikit-learn  TruncatedSVD()
scikit-learn 提供一個class 可執行上述降維的操作,
from sklearn.decomposition import TruncatedSVD
A = np.array([
[1,2,3,4,5,6,7,8,9,10],
[11,12,13,14,15,16,17,18,19,20],
[21,22,23,24,25,26,27,28,29,30]])
svd = TruncatedSVD(n_components=2)
svd.fit(A)
result = svd.transform(A)
print(result)
AA = svd.inverse_transform(result)
print(AA)
[[18.52157747  6.47697214]
 [49.81310011  1.91182038]
 [81.10462276 -2.65333138]]
[[ 1.  2.  3.  4.  5.  6.  7.  8.  9. 10.]
 [11. 12. 13. 14. 15. 16. 17. 18. 19. 20.]
 [21. 22. 23. 24. 25. 26. 27. 28. 29. 30.]]
程式將給的矩陣轉成$3\times 2$,之後再轉回原來的$3\times 10$,結果也是一樣。

參考文獻
Gavin Hackeling, Mastering Machine Learning with scikit-learn, 2nd, Packt Publishing, 2017

沒有留言:

張貼留言