本文章是看了B站大佬关于主成分分析的讲解后做的笔记，如果看不懂，建议移步观看大佬的视频https://www.bilibili.com/video/BV1E5411E71z

前言

我们先来看个二维的数据D，它有两个维度x和y，降维的一个准则就是数据在新的维度上要尽可能的分散。观察到原始数据，无论是在x轴方向上，还是y轴方向，数据的分散程度都不是最大的。因此，我们需要找到一个新的坐标系，使得数据在新坐标系上的一个投影方向离散度最大。PCA的目标就是找到这个新的坐标系，离散度最大的投影方向就是PCA的第一主成分。

数学定义

令 D = [ x 1 x 2 . . . x n y 1 y 2 . . . y n ] D=\\begin{bmatrix}x_1 & x_2 & … & x_n \\\\ y_1 & y_2 & … & y_n \\\\ \\end{bmatrix} D=[x1​y1​​x2​y2​​......​xn​yn​​], 去中心化之后 D ′ = [ x 1 ′ x 2 ′ . . . x n ′ y 1 ′ y 2 ′ . . . y n ′ ] D\’=\\begin{bmatrix}x_1\’ & x_2\’ & … & x_n\’ \\\\ y_1\’ & y_2\’ & … & y_n\’ \\\\ \\end{bmatrix} D′=[x1′​y1′​​x2′​y2′​​......​xn′​yn′​​]，其中 x i ′ = x i − x ‾ , y i ′ = y i − y ‾ . x_i\’=x_i-\\overline x, y_i\’=y_i-\\overline y. xi′​=xi​−x,yi′​=yi​−y​.

PCA的第一步——去中心化，目的是为了把数据的中心移到原点，这样我们在寻找坐标系的时候会更方便一些，后续的证明过程也会得到简化。
下面我们来看几个数学定义：

协 方 差 ： c o v ( x , y ) = ∑ i = 1 n ( x i − x ‾ , y i − y ‾ ) n − 1 方 差 ： c o v ( x , x ) = ∑ i = 1 n ( x i − x ‾ ) 2 n − 1 协 方 差 矩 阵 ： C = [ c o v ( x , x ) c o v ( x , y ) c o v ( y , x ) c o v ( y , y ) ] \\begin{aligned} 协方差&：cov(x,y)=\\frac {\\sum_{i=1}^n(x_i-\\overline x,y_i-\\overline y)}{n-1} \\\\ 方差&：cov(x,x)=\\frac {\\sum_{i=1}^n(x_i-\\overline x)^2}{n-1}\\\\ 协方差矩阵&：C=\\begin{bmatrix}cov(x,x) & cov(x,y) \\\\ cov(y,x) & cov(y,y)\\\\ \\end{bmatrix}\\end{aligned} 协方差方差协方差矩阵​：cov(x,y)=n−1∑i=1n​(xi​−x,yi​−y​)​：cov(x,x)=n−1∑i=1n​(xi​−x)2​：C=[cov(x,x)cov(y,x)​cov(x,y)cov(y,y)​]​去中心化之后，均值为0，那么可以得到协方差矩阵的一个特殊形式为：

C = [ ∑ i = 1 n ( x i ) 2 n − 1 ∑ i = 1 n ( x i , y i ) n − 1 ∑ i = 1 n ( y i , x i ) n − 1 ∑ i = 1 n ( y i ) 2 n − 1 ] = 1 n − 1 [ ∑ i = 1 n ( x i ) 2 ∑ i = 1 n ( x i , y i ) ∑ i = 1 n ( y i , x i ) ∑ i = 1 n ( y i ) 2 ] = 1 n − 1 W W T C=\\begin{bmatrix}\\frac {\\sum_{i=1}^n(x_i)^2}{n-1} & \\frac {\\sum_{i=1}^n(x_i,y_i)}{n-1} \\\\ \\frac {\\sum_{i=1}^n(y_i,x_i)}{n-1} & \\frac {\\sum_{i=1}^n(y_i)^2}{n-1}\\\\ \\end{bmatrix} =\\frac {1}{n-1}\\begin{bmatrix}\\sum_{i=1}^n(x_i)^2 & \\sum_{i=1}^n(x_i,y_i) \\\\ \\sum_{i=1}^n(y_i,x_i) & \\sum_{i=1}^n(y_i)^2\\\\ \\end{bmatrix} = \\frac {1}{n-1}WW^T C=[n−1∑i=1n​(xi​)2​n−1∑i=1n​(yi​,xi​)​​n−1∑i=1n​(xi​,yi​)​n−1∑i=1n​(yi​)2​​]=n−11​[∑i=1n​(xi​)2∑i=1n​(yi​,xi​)​∑i=1n​(xi​,yi​)∑i=1n​(yi​)2​]=n−11​WWT

其中$W$表示任意方向均值都为0的数据。
另外，定义缩放矩阵 S = [ a 0 0 b ] S=\\begin{bmatrix}a & 0 \\\\ 0 & b\\\\ \\end{bmatrix} S=[a0​0b​]， 其中a表示数据在x轴方向拉伸a倍，b表示数据在y轴方向拉伸b倍。

旋转矩阵 R = [ s i n ( θ ) − s i n ( θ ) s i n ( θ ) c o s ( θ ) ] R=\\begin{bmatrix}sin(\\theta) &-sin(\\theta) \\\\ sin(\\theta) & cos(\\theta)\\\\ \\end{bmatrix} R=[sin(θ)sin(θ)​−sin(θ)cos(θ)​]，代表数据以原点为中心旋转了$\theta$度。

我们假设白数据为矩阵$W$，其中白数据的特点是在任意方向上的分布都服从标准正态分布，同时认为中心位于原点的任意数据都可以由白数据经过线性变换（缩放、旋转）得到。那么去中心化之后的数据就可以用矩阵乘法来表示： D ′ = R S W D\’=RSW D′=RSW

理论推导

由之前的定义可知白数据的协方差矩阵：

C W = [ c o v ( x , x ) c o v ( x , y ) c o v ( y , x ) c o v ( y , y ) ] = 1 n − 1 W W T C_W=\\begin{bmatrix}cov(x,x) & cov(x,y) \\\\ cov(y,x) & cov(y,y)\\\\ \\end{bmatrix}=\\frac {1}{n-1}WW^T CW​=[cov(x,x)cov(y,x)​cov(x,y)cov(y,y)​]=n−11​WWT

白数据服从标准正态分布，故有
$$\begin{gathered} cov(x,x)=cov(y,y)=1 \\ cov(x,y)=cov(y,x)=0 \end{gathered}$$
代入得
C W = [ c o v ( x , x ) c o v ( x , y ) c o v ( y , x ) c o v ( y , y ) ] = [ 1 0 0 1 ] = I C_W=\\begin{bmatrix}cov(x,x) & cov(x,y) \\\\ cov(y,x) & cov(y,y)\\\\ \\end{bmatrix} =\\begin{bmatrix}1 & 0 \\\\ 0 & 1\\\\ \\end{bmatrix}=I CW​=[cov(x,x)cov(y,x)​cov(x,y)cov(y,y)​]=[10​01​]=I

故 C W = 1 n − 1 W W T = I C_W=\\frac {1}{n-1}WW^T=I CW​=n−11​WWT=I
然后我们现在来化简去中心化之后的数据的协方差矩阵：
C D ′ = 1 n − 1 D ′ D ′ T = 1 n − 1 ( R S W ) ( R S W ) T = 1 n − 1 R S W W T S T R T = R S ( 1 n − 1 W W T ) S T R T = R S S T R T \\begin{aligned} C_{D\’} &=\\frac {1}{n-1}D\'{D\’}^T \\\\ &=\\frac {1}{n-1}(RSW)(RSW)^T \\\\ &=\\frac {1}{n-1}RSWW^TS^TR^T \\\\ &=RS(\\frac {1}{n-1}WW^T)S^TR^T \\\\ &=RSS^TR^T \\end{aligned} CD′​​=n−11​D′D′T=n−11​(RSW)(RSW)T=n−11​RSWWTSTRT=RS(n−11​WWT)STRT=RSSTRT​
这样数据的协方差公式就可以用缩放矩阵和旋转矩阵来表示了，而这个旋转矩阵其实就是我们要找的新坐标系。
根据缩放矩阵和旋转矩阵的特性：
S T = S R T = R − 1 \\begin{aligned} S^T &=S \\\\ R^T &=R^{-1} \\end{aligned} STRT​=S=R−1​
有 C D ′ = R S S T R T = R S S R − 1 C_{D\’}=RSS^TR^T=RSSR^{-1} CD′​=RSSTRT=RSSR−1
令 L = S S T = S S = S 2 L=SS^T=SS=S^2 L=SST=SS=S2
得 C D ′ = R S S T R T = R S S R − 1 = R L R − 1 C_{D\’}=RSS^TR^T=RSSR^{-1}=RLR^{-1} CD′​=RSSTRT=RSSR−1=RLR−1
矩阵的特征向量定义为
C D ′ v = λ v C_{D\’}v=\\lambda v CD′​v=λv
以二维矩阵为例
C D ′ v 1 = λ 1 v 1 C D ′ v 2 = λ 2 v 2 C_{D\’}v_1=\\lambda_1 v_1 \\\\ C_{D\’}v_2=\\lambda_2 v_2 CD′​v1​=λ1​v1​CD′​v2​=λ2​v2​
写成矩阵的形式就是
C D ′ [ v 1 v 2 ] = [ v 1 v 2 ] [ λ 1 0 0 λ 2 ] C_{D\’} \\begin{bmatrix}v_1 & v_2\\\\ \\end{bmatrix}=\\begin{bmatrix}v_1 & v_2\\\\ \\end{bmatrix} \\begin{bmatrix}\\lambda_1& 0 \\\\ 0 & \\lambda_2\\\\ \\end{bmatrix} CD′​[v1​​v2​​]=[v1​​v2​​][λ1​0​0λ2​​]
由于特征向量为单位向量，同时满足正交性，故可令旋转矩阵 R = [ v 1 v 2 ] R=\\begin{bmatrix}v_1 & v_2\\\\ \\end{bmatrix} R=[v1​​v2​​]，同时 L = [ λ 1 0 0 λ 2 ] L= \\begin{bmatrix}\\lambda_1 & 0\\\\ 0 & \\lambda_2\\\\ \\end{bmatrix} L=[λ1​0​0λ2​​]
则有 C D ′ R = R L C_{D\’}R=RL CD′​R=RL，推出 C D ′ = R L R − 1 C_{D\’}=RLR^{-1} CD′​=RLR−1
故可得旋转矩阵$R$就是协方差矩阵的特征向量。

python代码实现

我们根据原理来实现一下PCA的过程：

def my_pca(X, n_components):    X = X - np.mean(X, axis=0)  # 减去每个方向的均值    n = X.shape[0]    X_cov = 1/(n-1)*X.T@X    L, R = np.linalg.eig(X_cov)  # 求协方差矩阵的特征值L和特征向量R    index = np.argsort(L)[::-1]    R_rank = R[:, index]  # 根据特征值大小对特征向量R进行排序    new_X = X@R_rank  # 矩阵相乘得到数据在新坐标系下的变换    return new_X[:, :n_components]  # 返回变换后的数据

然后随机生成100行4维的数据测试一下：

data = np.random.random([100, 4])X_1 = my_pca(data, n_components=2)  # 数据降到2维print(X_1[:5, :])

输出降维后的前5条数据：

[[-0.2606472  -0.17418523] [-0.42146202  0.02135193] [-0.32852065 -0.31256129] [ 0.00865476  0.0472365 ] [-0.57729244  0.33415261]]

调用sklearn包的PCA方法：

X_2 = PCA(n_components=2).fit_transform(data)print(X_2[:5, :])

可以看到结果是一样的，说明我们的计算没有错

[[-0.2606472  -0.17418523] [-0.42146202  0.02135193] [-0.32852065 -0.31256129] [ 0.00865476  0.0472365 ] [-0.57729244  0.33415261]]

总结

最后总结一下PCA的步骤：

输入：原始数据D，降维维度n1、去中心化，得到D\'  2、求协方差矩阵C\'3、求特征向量R和特征值L4、按特征值从大到小排序，得到对应排序的特征向量R\'5、计算变换坐标系后的数据D_P=R\'D\'输出：前n个主成分D_P[0:n]