本文章是看了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=[x1y1x2y2......xnyn], 去中心化之后 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)2n−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−11WWT
其中
W
W
W表示任意方向均值都为0的数据。
另外，定义缩放矩阵
S
=
[
a
0
0
b
]
S=\begin{bmatrix}a & 0 \\ 0 & b\\ \end{bmatrix}
S=[a00b]， 其中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
W
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−11WWT
白数据服从标准正态分布，故有
c
o
v
(
x
,
x
)
=
c
o
v
(
y
,
y
)
=
1
c
o
v
(
x
,
y
)
=
c
o
v
(
y
,
x
)
=
0
cov(x,x)=cov(y,y)=1 \\ cov(x,y)=cov(y,x)=0
cov(x,x)=cov(y,y)=1cov(x,y)=cov(y,x)=0
代入得
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)]=[1001]=I
故
C
W
=
1
n
−
1
W
W
T
=
I
C_W=\frac {1}{n-1}WW^T=I
CW=n−11WWT=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−11D′D′T=n−11(RSW)(RSW)T=n−11RSWWTSTRT=RS(n−11WWT)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
C_{D'}=RSS^TR^T=RSSR^{-1}
CD′=RSSTRT=RSSR−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=λ1v1CD′v2=λ2v2
写成矩阵的形式就是
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′[v1v2]=[v1v2][λ100λ2]
由于特征向量为单位向量，同时满足正交性，故可令旋转矩阵
R
=
[
v
1
v
2
]
R=\begin{bmatrix}v_1 & v_2\\ \end{bmatrix}
R=[v1v2]，同时
L
=
[
λ
1
0
0
λ
2
]
L= \begin{bmatrix}\lambda_1 & 0\\ 0 & \lambda_2\\ \end{bmatrix}
L=[λ100λ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
R
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，降维维度n
1、去中心化，得到D'
2、求协方差矩阵C'
3、求特征向量R和特征值L
4、按特征值从大到小排序，得到对应排序的特征向量R'
5、计算变换坐标系后的数据D_P=R'D'
输出：前n个主成分D_P[0:n]