一、PCA算法的原理
PCA（principle component analysis），即主成分分析法，是一個非監督的機器學習算法，是一種用於探索高維數據結構的技術，主要用於對數據的降維，通過降維可以發現更便於人理解的特征，加快對樣本有價值信息的處理速度，此外還可以應用於可視化（降到二維）和去噪。
1、PCA與LDA算法的基本思想
數據從原來的坐標系轉換到新的坐標系，新坐標系的選擇是由數據本身決定的。第一個新坐標軸選擇的是原始數據中方差最大的方向，第二個新坐標軸選擇和第一個坐標軸正交且具有最大方差的方向。該過程一直重復，重復次數為原始數據中特征的數目。我們會發現，大部分方差都包含在最前面的幾個新坐標軸中。因此，我們可以忽略余下的坐標軸，即對數據進行降維處理。
2、數學推導過程
PCA本質上是將方差最大的方向作為主要特征，並且在各個正交方向上將數據“離相關”，也就是讓它們在不同正交方向上沒有相關性。
求解思路：用方差來定義樣本的間距，方差越大表示樣本分布越稀疏，方差越小表示樣本分布越密集。
方差的公式如下：
在求解最大方差前，為了方便計算，可以先對樣本進行demean（去均值）處理，即減去每個特征的均值，這種處理方式不會改變樣本的相對分布（效果就像坐標軸進行了移動）。去均值后，樣本x每個特征維度上的均值都是0，方差的公式轉換下圖的公式：
在這里，
代表已經經過映射后的某樣本。
對於只有2個維度的樣本，現在的目標就是：求一個軸的方向w=（w1,w2），使得映射到w方向后，方差最大。
目標函數表示如下：
為求解此問題，需要使用梯度上升算法，梯度的求解公式如下：
3、PCA算法流程:
（1）去平均值，即每一位特征減去各自的平均值；
（2）計算協方差矩陣；
（3）計算協方差矩陣的特征值與特征向量；
（4）對特征值從大到小排序；
（5）保留最大的個特征向量；
（6）將數據轉換到個特征向量構建的新空間中。
4、PCA算法實現一般流程：
（1）對數據進行歸一化處理；
（2）計算歸一化后的數據集的協方差矩陣；
（3）計算協方差矩陣的特征值和特征向量；
（4）保留最重要的k個特征（通常k要小於n）；
（5）找出k個特征值相應的特征向量
（6）將m * n的數據集乘以k個n維的特征向量的特征向量（n * k）,得到最后降維的數據。
5、PCA降維准則：
（1） 最近重構性：樣本集中所有點，重構后的點距離原來的點的誤差之和最小。
（2） 最大可分性：樣本在低維空間的投影盡可能分開。
6、PCA算法優點：
（1）使得數據集更易使用；
（2）降低算法的計算開銷；
（3）去除噪聲；
（4）使得結果容易理解；
（5）完全無參數限制。
7、PCA算法缺點：
（1）如果用戶對觀測對象有一定的先驗知識，掌握了數據的一些特征，卻無法通過參數化等方法對處理過程進行干預，可能會得不到預期的效果，效率也不高；
（2） 特征值分解有一些局限性，比如變換的矩陣必須是方陣；
（3） 在非高斯分布情況下，PCA方法得出的主元可能並不是最優的。
8、PCA算法應用：
（1）高維數據集的探索與可視化。
（2）數據壓縮。
（3）數據預處理。
（4）圖象、語音、通信的分析處理。
（5）降維(最主要)，去除數據冗余與噪聲。
二、代碼實現
1.自己實現的PCA算法（不使用sklearn）
import numpy as np import matplotlib.pyplot as plt X=np.empty((100,2)) X[:,0]=np.random.uniform(0,100,size=100) X[:,1]=0.75*X[:,0]+3+np.random.normal(0,10,size=100) plt.scatter(X[:,0],X[:,1]) def demean(X): return X-np.mean(X,axis=0) X_demean=demean(X) plt.figure(2) plt.scatter(X_demean[:,0],X_demean[:,1]) #print(np.mean(X[:,0])) #print(np.mean(X_deman[:,0])) #print(np.mean(X_deman[:,1])) def f(w,X): return np.sum((X.dot(w)**2))/len(X) def df_math(w,X): return X.T.dot(X.dot(w))*2/len(X) def direction(w): return w / np.linalg.norm(w) def gradient_ascent(df, X, initial_w, eta, n_iters = 1e4, epsilon=1e-8): w = direction(initial_w) cur_iter = 0 while cur_iter < n_iters: gradient = df(w, X) last_w = w w = w + eta * gradient w = direction(w) # 注意1：每次求一個單位方向 if(abs(f(w, X) - f(last_w, X)) < epsilon): break cur_iter += 1 return w initial_w = np.random.random(X.shape[1]) # 注意2：不能用0向量開始 eta = 0.001 w = gradient_ascent(df_math, X_demean, initial_w, eta) plt.figure(3) plt.scatter(X_demean[:,0], X_demean[:,1]) #plt.plot([0, w[0]*30], [0, w[1]*30], color='r') plt.plot([0, w[0]*50], [0 , w[1]*50], color='r')
輸出結果：
2、PCA分類
import numpy as np import pandas as pd import matplotlib.pyplot as plt #計算均值,要求輸入數據為numpy的矩陣格式，行表示樣本數，列表示特征 def meanX(dataX): return np.mean(dataX,axis=0)#axis=0表示依照列來求均值。假設輸入list,則axis=1 """ 參數： - XMat：傳入的是一個numpy的矩陣格式，行表示樣本數，列表示特征 - k：表示取前k個特征值相應的特征向量 返回值： - finalData：參數一指的是返回的低維矩陣，相應於輸入參數二 - reconData：參數二相應的是移動坐標軸后的矩陣 """ def pca(XMat, k): average = meanX(XMat) m, n = np.shape(XMat) data_adjust = [] avgs = np.tile(average, (m, 1)) data_adjust = XMat - avgs covX = np.cov(data_adjust.T) #計算協方差矩陣 featValue, featVec= np.linalg.eig(covX) #求解協方差矩陣的特征值和特征向量 index = np.argsort(-featValue) #依照featValue進行從大到小排序 finalData = [] if k > n: print("k must lower than feature number") return else: #注意特征向量時列向量。而numpy的二維矩陣(數組)a[m][n]中，a[1]表示第1行值 selectVec = np.matrix(featVec.T[index[:k]]) #所以這里須要進行轉置 finalData = data_adjust * selectVec.T reconData = (finalData * selectVec) + average return finalData, reconData
#輸入文件的每行數據都以\t隔開 def loaddata(datafile): return np.array(pd.read_csv(datafile,sep=" ",header=-1)).astype(np.float)
def plotBestFit(data1, data2): dataArr1 = np.array(data1) dataArr2 = np.array(data2) m = np.shape(dataArr1)[0] axis_x1 = [] axis_y1 = [] axis_x2 = [] axis_y2 = [] for i in range(m): axis_x1.append(dataArr1[i,0]) axis_y1.append(dataArr1[i,1]) axis_x2.append(dataArr2[i,0]) axis_y2.append(dataArr2[i,1]) fig = plt.figure() ax = fig.add_subplot(111) ax.scatter(axis_x1, axis_y1, s=50, c='red', marker='s') ax.scatter(axis_x2, axis_y2, s=50, c='blue') plt.xlabel('x1'); plt.ylabel('x2'); plt.savefig("outfile.png") plt.show()
#依據數據集data.txt def main(): datafile = "data.txt" XMat = loaddata(datafile) k = 2 return pca(XMat, k)
if __name__ == "__main__": finalData, reconMat = main() plotBestFit(finalData, reconMat)
運行結果：
