回顧了下PCA的步驟，并用python實現(xiàn)。深刻的發(fā)現(xiàn)當年學的特征值、特征向量好強大。

Introduction to PCA

PCA是一種無監(jiān)督的學習方式，是一種很常用的降維方法。在數(shù)據(jù)信息損失最小的情況下，將數(shù)據(jù)的特征數(shù)量由n，通過映射到另一個空間的方式，變?yōu)閗(k<n)。

Sample Data

這里用一個2維的數(shù)據(jù)來說明PCA，選擇2維的數(shù)據(jù)是因為2維的比較容易畫圖。
這是數(shù)據(jù)：

data

import numpy as np
x=np.array([2.5,0.5,2.2,1.9,3.1,2.3,2,1,1.5,1.1])
y=np.array([2.4,0.7,2.9,2.2,3,2.7,1.6,1.1,1.6,0.9])
data=np.matrix([[scaled_x[i],scaled_y[i]] for i in range(len(scaled_x))])

Step 1: 求平均值以及做normalization

mean_x=np.mean(x)
mean_y=np.mean(y)
scaled_x=x-mean_x
scaled_y=y-mean_y
data=np.matrix([[scaled_x[i],scaled_y[i]] for i in range(len(scaled_x))])

畫個圖看看分布情況：

import matplotlib.pyplot as plt
plt.plot(scaled_x,scaled_y,'o')    

data distribution

Step 2: 求協(xié)方差矩陣(Covariance Matrix)

• 協(xié)方差矩陣
  方差的定義為：

variance

協(xié)方差的定義為：

covariance

假設(shè)n為數(shù)據(jù)的特征數(shù)，那么協(xié)方差矩陣M, 為一個nn的矩陣，其中Mij為第i和第j個特征的協(xié)方差，對角線是各個特征的方差。
在我們的數(shù)據(jù)中，n=2，所以協(xié)方差矩陣是22的，
通過numpy我們可以很方便的得到：

cov=np.cov(scaled_x,scaled_y)

得到cov的結(jié)果為：
array([[ 0.61655556, 0.61544444],
[ 0.61544444, 0.71655556]])

• 散度矩陣 Scatter Matrix
  另一種選擇是用scatter matrix，定義為：

sactter matrix

由于我們之前已經(jīng)做過normalization，因此對于我們來說，
這個矩陣就是 data*data的轉(zhuǎn)置矩陣。

np.dot(np.transpose(data),data)

得到結(jié)果：
matrix([[ 5.549, 5.539],
[ 5.539, 6.449]])

我們發(fā)現(xiàn)，其實協(xié)方差矩陣和散度矩陣關(guān)系密切，散度矩陣 就是協(xié)方差矩陣乘以（總數(shù)據(jù)量-1）。因此他們的特征根和特征向量是一樣的。這里值得注意的一點就是，散度矩陣是SVD奇異值分解的一步，因此PCA和SVD是有很大聯(lián)系的，他們的關(guān)系這里就不詳細談了，以后有機會再寫下。

Step 3: 求協(xié)方差矩陣的特征根和特征向量

用numpy計算特征根和特征向量很簡單，

eig_val, eig_vec = np.linalg.eig(cov)

但是他們代表的意義非常有意思，讓我們將特征向量加到我們原來的圖里：

plt.plot(scaled_x,scaled_y,'o',)
xmin ,xmax = scaled_x.min(), scaled_x.max()
ymin, ymax = scaled_y.min(), scaled_y.max()
dx = (xmax - xmin) * 0.2
dy = (ymax - ymin) * 0.2
plt.xlim(xmin - dx, xmax + dx)
plt.ylim(ymin - dy, ymax + dy)
plt.plot([eig_vec[:,0][0],0],[eig_vec[:,0][1],0],color='red')
plt.plot([eig_vec[:,1][0],0],[eig_vec[:,1][1],0],color='red')

data distribution with eig_vec

其中紅線就是特征向量。有幾點值得注意：

• 特征向量之間是正交的，PCA其實就是利用特征向量的這個特點，重新構(gòu)建新的空間體系

• 特征向量代表著數(shù)據(jù)的pattern(模式),比如一條代表著y隨著x的增大而增大的趨勢，而另外一條，則是代表數(shù)據(jù)也有該方面的變化。所以特征向量的命名是很科學的，他代表著矩陣的特征。
  如果我們將數(shù)據(jù)直接乘以特征向量矩陣，那么其實我們只是以特征向量為基底，重新構(gòu)建了空間，畫個圖，感受下吧：

  new_data=np.transpose(np.dot(eig_vec,np.transpose(data)))
  

new_data

藍色的三角形就是經(jīng)過坐標變換后得到的新點，其實他就是紅色原點投影到紅線、藍線形成的。

Step 4: 選擇主要成分

得到特征值和特征向量之后，我們可以根據(jù)特征值的大小，從大到小的選擇K個特征值對應(yīng)的特征向量。
這個用python的實現(xiàn)也很簡單：

eig_pairs = [(np.abs(eig_val[i]), eig_vec[:,i]) for i in range(len(eig_val))]
eig_pairs.sort(reverse=True)

從eig_pairs選取前k個特征向量就行。這里，我們只有兩個特征向量，選一個最大的。

feature=eig_pairs[0][1]

Step 5: 轉(zhuǎn)化得到降維的數(shù)據(jù)

主要將原來的數(shù)據(jù)乘以經(jīng)過篩選的特征向量組成的特征矩陣之后，就可以得到新的數(shù)據(jù)了。

new_data_reduced=np.transpose(np.dot(feature,np.transpose(data)))

output：

new_data

數(shù)據(jù)果然變成了一維的數(shù)據(jù)。
最后我們通過畫圖來理解下數(shù)據(jù)經(jīng)過PCA到底發(fā)生了什么。

plt.plot(scaled_x,scaled_y,'o',color='red')
plt.plot([eig_vec[:,0][0],0],[eig_vec[:,0][1],0],color='red')
plt.plot([eig_vec[:,1][0],0],[eig_vec[:,1][1],0],color='blue')
plt.plot(new_data[:,0],new_data[:,1],'^',color='blue')
plt.plot(new_data_reduced[:,0],[1.2]*10,'*',color='green')

final_result

綠色的五角星是PCA處理過后得到的一維數(shù)據(jù)，為了能跟以前的圖對比，將他們的高度定位1.2，其實就是紅色圓點投影到藍色線之后形成的點。這就是PCA,通過選擇特征根向量，形成新的坐標系，然后數(shù)據(jù)投影到這個新的坐標系，在盡可能少的丟失信息的基礎(chǔ)上實現(xiàn)降維。

總結(jié)

通過上述幾步的處理，我們簡單的實現(xiàn)了PCA第一個2維數(shù)據(jù)的處理，但是原理就是這樣，我們可以很輕易的就依此實現(xiàn)多維的。

PCA的python3實現(xiàn)

import numpy as np
def pca(X,k):#k is the components you want
  #mean of each feature
  n_samples, n_features = X.shape
  mean=np.array([np.mean(X[:,i]) for i in range(n_features)])
  #normalization
  norm_X=X-mean
  #scatter matrix
  scatter_matrix=np.dot(np.transpose(norm_X),norm_X)
  #Calculate the eigenvectors and eigenvalues
  eig_val, eig_vec = np.linalg.eig(scatter_matrix)
  eig_pairs = [(np.abs(eig_val[i]), eig_vec[:,i]) for i in range(n_features)]
  # sort eig_vec based on eig_val from highest to lowest
  eig_pairs.sort(reverse=True)
  # select the top k eig_vec
  feature=np.array([ele[1] for ele in eig_pairs[:k]])
  #get new data
  data=np.dot(norm_X,np.transpose(feature))
  return data

用sklearn的PCA與我們的pca做個比較：

from sklearn.decomposition import PCA
X = np.array([[-1, 1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]])
pca=PCA(n_components=1)
pca.fit(X)
pca.transform(X)

得到結(jié)果：

array([[-0.50917706],
   [-2.40151069],
   [-3.7751606 ],
   [ 1.20075534],
   [ 2.05572155],
   [ 3.42937146]])

用我們的pca試試

pca(X,1)

得到結(jié)果：

array([[-0.50917706],
   [-2.40151069],
   [-3.7751606 ],
   [ 1.20075534],
   [ 2.05572155],
   [ 3.42937146]])

完全一致，完美~
值得一提的是，sklearn中PCA的實現(xiàn)，用了部分SVD的結(jié)果，果然他們因緣匪淺。

To Be Continued:

1. The relationship between PCA and SVD
2. The implementation of PCA, and the comparison between the results of my PCA and sklearn's PCA (completed)