一、PCA（Principal Component Analysis）介紹
PCA是數據處理中的一個常用方法，用於數據降維，特征提取等，實際上是將在原有的特征空間中分布的數據映射到新的特征空間（或者說，將原有到正交坐標系進行旋轉，使得在旋轉后的坐標系下，在某幾根坐標軸上數據分布的方差比較大。在這些特征上，不同數據區分度比較高，最高的我們稱為第一主成份，按照大小順序依次第二主成份，第三主成份）。對於高維數據，我們只需要保留k個主要的投影結果，即模最大也就是對應奇異值最大的投影方向，因此投影矩陣也只需要保留前k個主成分對應的特征向量，故 V 也只需要n×k維即可。對數據降維。經過PCA分解以后的分量都是正交的，因此是不相關的，但不一定獨立。經過重新映射后的數據不再具有最初的物理意義（都打亂了）。
思考：我們如何得到這些包含最大差異性的主成分方向呢？
答案：計算數據矩陣去均值后的協方差矩陣，得到協方差矩陣的特征值、特征向量，選擇特征值最大(即方差最大)的k個特征所對應的特征向量組成的矩陣（映射矩陣），將映射矩陣轉置后乘以去均值的數據矩陣得到降維的數據。
二、數學推導
得到協方差矩陣的特征值和特征向量有兩種方法：1特征值分解協方差矩陣 。2奇異值分解協方差矩陣。
所以，PCA算法有兩種實現方法：基於特征值分解協方差矩陣實現PCA算法、基於SVD分解協方差矩陣實現PCA算法。(矩陣論有關知識)
步驟：
輸入：數據集 
，需要降到k維。
1) 去平均值(即去中心化)，即每一位特征減去各自的平均值。
2) 計算協方差矩陣
（注：這里除或不除樣本數量n或n-1,其實對求出的特征向量沒有影響）。
3) 用特征值分解/SVD分解方法求協方差矩陣的特征值與特征向量。
4) 對特征值從大到小排序，選擇其中最大的k個。對應的k個特征向量分別作為行向量組成特征向量矩陣P。
5) 將數據轉換到k個特征向量構建的新空間中，即去中心化的數據*協方差矩陣的特征向量。
注釋：
1 協方差矩陣
我們可以將矩陣看作作用於所有數據並朝向某個方向的力。同時我們還知道了變量間的相關性可以由方差和協方差表達，並且我們希望保留最大方差以實現最優的降維。因此我們希望能將方差和協方差統一表示，並且兩者均可以表示為內積的形式，而內積又與矩陣乘法密切相關。因此我們可以采用矩陣乘法的形式表示。若輸入矩陣 X 有兩個特征 a 和 b，且共有 m 個樣本，那么有：
如果我們用 X 左乘 X 的轉置，那么就可以得出協方差矩陣
這個矩陣對角線上的兩個元素分別是兩特征的方差，而其它元素是 a 和 b 的協方差。兩者被統一到了一個矩陣的，因此我們可以利用協方差矩陣描述數據點之間的方差和協方差，即經驗性地描述我們觀察到的數據。
2 PCA降維之前為什么要先標准化？
PCA(主成分分析)所對應的數學理論SVD(矩陣的奇異值分解本身是完全不需要對矩陣中的元素做標准化或者去中心化的。
但是對於機器學習，我們通常會對矩陣（也就是數據）的每一列先進行標准化。
PCA通常是用於高維數據的降維，它可以將原來高維的數據投影到某個低維的空間上並使得其方差盡量大。如果數據其中某一特征（矩陣的某一列）的數值特別大，那么它在整個誤差計算的比重上就很大，那么可以想象在投影到低維空間之后，為了使低秩分解逼近原數據，整個投影會去努力逼近最大的那一個特征，而忽略數值比較小的特征。因為在建模前我們並不知道每個特征的重要性，這很可能導致了大量的信息缺失。為了“公平”起見，防止過分捕捉某些數值大的特征，我們會對每個特征先進行標准化處理，使得它們的大小都在相同的范圍內，然后再進行PCA。
此外，從計算的角度講，PCA前對數據標准化還有另外一個好處。因為PCA通常是數值近似分解，而非求特征值、奇異值得到解析解，所以當我們使用梯度下降等算法進行PCA的時候，我們最好先要對數據進行標准化，這是有利於梯度下降法的收斂。
三、PCA實現
1.python
（1）PCA的Python實現：
##Python實現PCA 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 #輸入數據 X = np.array([[-1, 1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]]) print(pca(X,1))
結果：
（2）sklearn的PCA實現：
1.sklearn.decomposition.PCA參數介紹
1）n_components：指定PCA降維后的特征維數。n_components是一個（0，1]之間的數。將參數設置為"mle", 此時PCA類會用MLE算法根據特征的方差分布情況自己去選擇一定數量的主成分特征來降維。默認值為：n_components=min(樣本數，特征數)。
2）whiten ：判斷是否進行白化。所謂白化，就是對降維后的數據的每個特征進行歸一化，讓方差都為1.對於PCA降維本身來說，一般不需要白化。如果你PCA降維后有后續的數據處理動作，可以考慮白化。默認值是False，即不進行白化。
3）svd_solver：指定奇異值分解SVD的方法。一般的PCA庫都是基於SVD實現的。有4個可以選擇的值：{‘auto’, ‘full’, ‘arpack’, ‘randomized’}。randomized一般適用於數據量大，數據維度多同時主成分數目比例又較低的PCA降維，它使用了一些加快SVD的隨機算法。 full則是傳統意義上的SVD，使用了scipy庫對應的實現。arpack和randomized的適用場景類似，區別是randomized使用的是scikit-learn自己的SVD實現，而arpack直接使用了scipy庫的sparse SVD實現。默認是auto，即PCA類會自己去在前面講到的三種算法里面去權衡，選擇一個合適的SVD算法來降維。一般來說，使用默認值就夠了。
explained_variance_：它代表降維后的各主成分的方差值。方差值越大，則說明越是重要的主成分。
explained_variance_ratio_：它代表降維后的各主成分的方差值占總方差值的比例，這個比例越大，則越是重要的主成分。
另：
fit(): Method calculates the parameters μ and σ and saves them as internal objects.
解釋：求得訓練集X的均值，方差，最大值，最小值,這些訓練集X固有的屬性。
transform(): Method using these calculated parameters apply the transformation to a particular dataset.
解釋：在fit的基礎上，進行標准化，降維，歸一化等操作（看具體用的是哪個工具，如PCA，StandardScaler等）。
fit_transform(): joins the fit() and transform() method for transformation of dataset.
解釋：fit_transform是fit和transform的組合，既包括了訓練又包含了轉換。
transform()和fit_transform()二者的功能都是對數據進行某種統一處理（比如標准化~N(0,1)，將數據縮放(映射)到某個固定區間，歸一化，正則化等）
fit_transform(trainData)對部分數據先擬合fit，找到該part的整體指標，如均值、方差、最大值最小值等等（根據具體轉換的目的），然后對該trainData進行轉換transform，從而實現數據的標准化、歸一化等等。
2.程序
##用sklearn的PCA import numpy as np 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) newX = pca.fit_transform(X) invX = pca.inverse_transform(X) # 將降維后的數據轉換成原始數據 print(newX) print(invX) print(pca.explained_variance_ratio_)
結果：
對比結果：
sklearn中的PCA是通過svd_flip函數實現的，sklearn對奇異值分解結果進行了一個處理，因為ui*σi*vi=(-ui)*σi*(-vi)，也就是u和v同時取反得到的結果是一樣的，而這會導致通過PCA降維得到不一樣的結果（雖然都是正確的）。具體了解可以分析一下sklearn中關於PCA的源碼。
出現問題：路徑種含有中文
解決方法：windows10系統修改c盤user文件夾下的計算機名稱
2.matlab
（1）PCA的matlab實現：
注意：直接調用PCA函數時不包括對數據標准化。
load('data.mat');
X=traces(1:1000,:);%1000X6000
[m,n]=size(X); %計算矩陣的行m和列n
%-------------第一步：標准化矩陣-----------------%
mv=mean(X); %計算各變量的均值
st=std(X); %計算各變量的標准差
X=(X-repmat(mv,m,1))./repmat(st,m,1); %標准化矩陣X
%z=zscore(X)
%-------------第二步：計算相關系數矩陣-----------------%
% R1=X'*X/(m-1); %方法一：協方差矩陣計算公式
% R2=cov(X); %方法二：協方差矩陣計算函數
R=corrcoef(X); %方法三：相關系數矩陣函數
%-------------第三步：計算特征向量和特征值-----------------%
[V,D]=eig(R); %計算矩陣R的特征向量矩陣V和特征值矩陣D,特征值由小到大
V=(rot90(V))'; %將特征向量矩陣V從大到小排序
D=rot90(rot90(D)); %將特征值矩陣由大到小排序
E=diag(D); %將特征值矩陣轉換為特征值向量
%-------------第四步：計算貢獻率和累計貢獻率-----------------%
ratio=0; %累計貢獻率
for k=1:n
r=E(k)/sum(E); %第k主成份貢獻率
ratio=ratio+r; %累計貢獻率
if(ratio>=0.9) %取累計貢獻率大於等於90%的主成分
break;
end
end
%-------------第五步：計算得分-----------------%
S=X*V;
load('atk_data.mat');
X=atk_traces;%能量跡條數 eps一個很小的正數:2.2204e-16
meanValue=ones(size(X,1),1)*mean(X);
X=X-meanValue;%每個維度減去該維度的均值
C=X'*X/(size(X,1)-1);%計算協方差矩陣
%計算特征向量，特征值
[V,D]=eig(C);
%將特征向量按降序排序
[dummy,order]=sort(diag(-D));
V=V(:,order);%將特征向量按照特征值大小進行降序排列
d=diag(D);%將特征值取出，構成一個列向量
newd=d(order);%將特征值構成的列向量按降序排列
T=V(:,1:600);%取前cols個特征向量，構成變換矩陣T
newX=X*T;%用變換矩陣T對X進行降維
（2）直接調用PCA的matlab實現：
load('data.mat');
X=traces(1:1000,:);
[coeff, score, latent, tsquared, explained, mu]=pca(X);
a=cumsum(latent)./sum(latent); % 計算特征的累計貢獻率
% explained和latent均可用來計算降維后取多少維度能夠達到自己需要的精度，且效果等價。
% explained=100*latent./sum(latent);
idx=find(a>0.9); % 將特征的累計貢獻率不小於0.9的維數作為PCA降維后特征的個數
k=idx(1);
Feature=score(:,1:k); % 取轉換后的矩陣score的前k列為PCA降維后特征
（3）fastPCA的matlab實現：
function [pcaA, V] = fastPCA(A,k)
% 快速PCA
% 輸入：A --- 樣本矩陣，每行為一個樣本
% k --- 降維至 k 維
% 輸出：pcaA --- 降維后的 k 維樣本特征向量組成的矩陣，每行一個樣本，列數 k 為降維后的樣本特征維數
% V --- 主成分向量
[r, ~] = size(A);
% 樣本均值
meanVec = mean(A);
% 計算協方差矩陣的轉置 covMatT
Z = (A-repmat(meanVec, r, 1));
covMatT = Z * Z';
% 計算 covMatT 的前 k 個本征值和本征向量
[V, ~] = eigs(covMatT, k);
% 得到協方差矩陣 (covMatT)' 的本征向量
V = Z' * V;
% 本征向量歸一化為單位本征向量
for i=1:k
V(:,i)=V(:,i)/norm(V(:,i));
end
% 線性變換（投影）降維至 k 維
pcaA = Z * V;
% 保存變換矩陣 V 和變換原點 meanVec
四、遇到的問題
1 如果出現NaN(Not-a-Number的簡寫，中文譯為“非數”，表示計算結果為不定),可能是因為分母為0，可以在數據后面加上 eps，不影響結果。
X=traces+ eps;%能量跡條數 eps一個很小的正數:2.2204e-16
2 未解決：調用函數是會莫名奇妙的小於想要的維數，用自定義才能求出來。 導致原因應該是協方差的計算方法不一樣。
% R1=X'*X/(size(x,1)-1); %方法一：協方差矩陣計算公式 % R2=cov(X); %方法二：協方差矩陣計算函數 % R=corrcoef(X); %方法三：相關系數矩陣函數
3.eig排序到底是從大到小還是從小到大 ？？？？
[V,D]=eig(R); %計算矩陣R的特征向量矩陣V和特征值矩陣D
文獻來源
1.理論：主成分分析（PCA）原理詳解_Microstrong0305的博客-CSDN博客_pca
2.代碼:12-.ipynb · ni1o1/pygeo-tutorial - Gitee.com