主成分分析(PCA)原理详解
1.PCA的概念
PCA(Principal Component Analysis),即主成分分析方法,是一种使用最广泛的数据降维算法。PCA的主要思想是将n维特征映射到k维上,这k维是全新的正交特征也被称为主成分,是在原有n维特征的基础上重新构造出来的k维特征。PCA的工作就是从原始的空间中顺序地找一组相互正交的坐标轴,新的坐标轴的选择与数据本身是密切相关的。其中,第一个新坐标轴选择是原始数据中方差最大的方向,第二个新坐标轴选取是与第一个坐标轴正交的平面中使得方差最大的,第三个轴是与第1,2个轴正交的平面中方差最大的。依次类推,可以得到n个这样的坐标轴。通过这种方式获得的新的坐标轴,我们发现,大部分方差都包含在前面k个坐标轴中,后面的坐标轴所含的方差几乎为0。于是,我们可以忽略余下的坐标轴,只保留前面k个含有绝大部分方差的坐标轴。事实上,这相当于只保留包含绝大部分方差的维度特征,而忽略包含方差几乎为0的特征维度,实现对数据特征的降维处理。
通过计算数据矩阵的协方差矩阵,然后得到协方差矩阵的特征值特征向量,选择特征值最大(即方差最大)的k个特征所对应的特征向量组成的矩阵。这样就可以将数据矩阵转换到新的空间当中,实现数据特征的降维。
得到协方差矩阵的特征值特征向量有两种方法:特征值分解协方差矩阵、奇异值分解协方差矩阵,所以PCA算法有两种实现方法:
- 基于特征值分解协方差矩阵实现PCA算法
- 基于SVD分解协方差矩阵实现PCA算法
2.协方差矩阵和散度矩阵
-
样本均值:
x ~ = 1 n ∑ i = 1 n a i \tilde{x}=\frac {1} {n} \sum_{i=1}^{n} {a_i} x~=n1i=1∑nai -
样本方差:
s 2 = 1 n − 1 ∑ i = 1 1 ( x i − x ~ ) 2 s^2=\frac{1} {n-1} \sum_{i=1}^{1}(x_i- \tilde{x})^2 s2=n−11i=1∑1(xi−x~)2 -
样本x和y的协方差:
c o v ( x , y ) = E [ ( x − E ( x ) ) ( y − E ( y ) ) ] cov(x,y)=E[(x-E(x))(y-E(y))] cov(x,y)=E[(x−E(x))(y−E(y))]
= 1 n − 1 ∑ i = 1 n ( x i − x ~ ) ( y i − y ~ ) = \frac {1} {n-1} \sum_{i=1}^{n} (x_i- \tilde{x})(y_i- \tilde y) =n−11i=1∑n(xi−x~)(yi−y~)- c o v ( x , y ) > 0 , x , y cov(x,y)>0,x,y cov(x,y)>0,x,y是正相关
- c o v ( x , y ) = 0 , x , y cov(x,y)=0,x,y cov(x,y)=0,x,y相互独立
- c o v ( x , y ) < 0 , x , y cov(x,y)<0,x,y cov(x,y)<0,x,y是负相关
-
当样本是 n n n维矩阵时,它们的协方差实际上是协方差矩阵(对称方阵)
-
对于三维矩阵 ( x , y , z ) (x,y,z) (x,y,z),它们的协方差是:
c o v ( x , y , z ) = ( c o v ( x , x ) c o v ( x , y ) c o v ( x , z ) c o v ( y , x ) c o v ( y , y ) c o v ( y , z ) c o v ( z , x ) c o v ( z , y ) c o v ( z , z ) ) cov(x,y,z)=\begin{pmatrix} cov(x,x) & cov(x,y) & cov(x,z) \\ cov(y,x) & cov(y,y) & cov(y,z) \\ cov(z,x) & cov(z,y) & cov(z,z) \\\end{pmatrix} cov(x,y,z)=⎝⎛cov(x,x)cov(y,x)cov(z,x)cov(x,y)cov(y,y)cov(z,y)cov(x,z)cov(y,z)cov(z,z)⎠⎞ -
散度矩阵
S = ∑ k − 1 n ( x k − x ~ ) ( x k − x ~ ) T S= \sum_{k-1}^{n}(x_k - \tilde x) (x_k - \tilde x)^T S=k−1∑n(xk−x~)(xk−x~)T
对于数据x的散度矩阵
X
X
T
XX^T
XXT
散度矩阵=协方差矩阵
×
\times
×(总数据量 -1)
因此散度矩阵和协方差矩阵的特征值相同,特征向量相同
其中散度矩阵也是SVD奇异分解的第一步
3.特征值分解矩阵
3.1特征值与特征向量
A
ν
=
λ
ν
A \nu=\lambda \nu
Aν=λν
其中:
ν
\nu
ν表示特征向量,
λ
\lambda
λ表示特征值
矩阵的一组特征向量是一组正交向量
3.2特征值分解矩阵
对于矩阵A,有一组特征向量
ν
\nu
ν,将这组特征向量进行正交化单位化,就能得到一组正交单位向量,特征值分解就是将矩阵A分解为如下格式:
A
=
Q
Σ
Q
A=Q \Sigma Q
A=QΣQ-1
其中,
Q
Q
Q是矩阵A的特征向量组成的矩阵,
Σ
\Sigma
Σ是一个对角阵,对角线上的元素就是特征值。
特征值分解的例子
分解如下矩阵:
A
=
(
−
1
1
0
−
4
3
0
1
0
2
)
A=\begin{pmatrix} -1 & 1 & 0 \\ -4 & 3 & 0 \\ 1 & 0 & 2 \\\end{pmatrix}
A=⎝⎛−1−41130002⎠⎞
首先有方阵的特征方程求出特征值.
∣
λ
E
−
A
∣
=
∣
−
1
−
λ
1
0
−
4
3
−
λ
0
1
0
2
−
λ
∣
|\lambda E-A|=\begin{vmatrix} -1-\lambda & 1 &0 \\ -4 & 3-\lambda & 0 \\ 1 & 0 & 2- \lambda \\ \end{vmatrix}
∣λE−A∣=∣∣∣∣∣∣−1−λ−4113−λ0002−λ∣∣∣∣∣∣
=
(
2
−
λ
)
(
λ
−
1
)
2
=
0
=(2-\lambda)(\lambda-1)^2=0
=(2−λ)(λ−1)2=0
特征值
λ
=
2
,
1
\lambda=2,1
λ=2,1(重根为2)
然后把每个特征值代入线性方程组
(
A
−
λ
E
)
x
=
0
(A-\lambda E)x=0
(A−λE)x=0,求出特征向量
当
λ
=
2
\lambda=2
λ=2时,解线性方程组
(
A
−
2
E
)
x
=
0
(A-2E)x=0
(A−2E)x=0
A
−
2
E
=
(
−
3
1
0
−
4
1
0
1
0
0
)
>
>
>
(
1
0
0
0
1
0
0
0
0
)
A-2E=\begin{pmatrix} -3 & 1 & 0 \\ -4 & 1 & 0 \\ 1 & 0 & 0 \\\end{pmatrix} >>> \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \\\end{pmatrix}
A−2E=⎝⎛−3−41110000⎠⎞>>>⎝⎛100010000⎠⎞
解得:
x
1
=
0
,
x
2
=
0
x_1=0,x_2=0
x1=0,x2=0
特征向量
p
1
=
(
0
0
0
)
p_1= \begin{pmatrix} 0 \\ 0 \\ 0 \\\end{pmatrix}
p1=⎝⎛000⎠⎞
当
λ
=
1
\lambda=1
λ=1时,解线性方程组:
(
A
−
2
E
)
x
=
0
(A-2E)x=0
(A−2E)x=0
(
A
−
2
E
)
=
(
−
2
1
0
−
4
2
0
1
0
1
)
>
>
>
(
1
0
1
0
1
2
0
0
0
)
(A-2E)=\begin{pmatrix} -2 & 1 & 0 \\ -4 & 2 & 0 \\ 1 & 0 & 1 \\\end{pmatrix} >>>\begin{pmatrix} 1 & 0 & 1 \\ 0 & 1 & 2 \\ 0 & 0 & 0 \\\end{pmatrix}
(A−2E)=⎝⎛−2−41120001⎠⎞>>>⎝⎛100010120⎠⎞
所以有:
{
x
1
+
x
3
=
0
x
2
+
2
x
3
=
0
\begin{cases} x_1+x_3=0 \\ x_2+2x_3=0 \end{cases}
{x1+x3=0x2+2x3=0
特征向量为
p
2
=
(
−
1
−
2
1
)
p_2=\begin{pmatrix} -1\\ -2 \\ 1 \\\end{pmatrix}
p2=⎝⎛−1−21⎠⎞
最后方阵A的特征值分解为:
A
=
Q
Σ
Q
−
1
=
(
1
0
1
0
1
2
0
0
0
)
(
2
0
0
0
1
0
0
0
1
)
(
0
−
1
−
1
0
−
2
−
2
1
1
1
)
A=Q \Sigma Q^{-1}=\begin{pmatrix} 1 & 0 & 1 \\ 0 & 1 & 2 \\ 0 & 0 & 0 \\\end{pmatrix} \begin{pmatrix} 2 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\\end{pmatrix} \begin{pmatrix} 0 & -1 & -1 \\ 0 & -2 & -2 \\ 1 & 1 & 1 \\\end{pmatrix}
A=QΣQ−1=⎝⎛100010120⎠⎞⎝⎛200010001⎠⎞⎝⎛001−1−21−1−21⎠⎞
4.SVD分解矩阵原理
特征值分解只适用于方阵,即
n
×
n
n \times n
n×n的矩阵,而在实际应用中,大部分都不是方阵
A
=
U
Σ
V
T
A=U \Sigma V^{T}
A=UΣVT
其中:
A:
m
×
n
m \times n
m×n矩阵
U
:
m
×
m
U:m \times m
U:m×m方阵,里面的正交向量被称为左奇异向量
Σ
:
m
×
n
\Sigma:m \times n
Σ:m×n矩阵,除对角线外,其他元素都是0,对角线上的元素称为奇异值
V
T
:
V
V^T:V
VT:V的转置矩阵,是
n
×
n
n \times n
n×n矩阵,里面的正交向量称为右奇异向量
一般来说,将
Σ
\Sigma
Σ上的值按从小到大排列
S V D SVD SVD分解矩阵A的步骤
- 求 A A T AA^T AAT的特征值和特征向量,用单位化的特征向量构成 U U U
- 求 A T A A^TA ATA的特征值和特征向量,用单位化的特征向量构成 V V V
- 将 A A T AA^T AAT或 A T A A^TA ATA的特征值求平方根,构成 Σ \Sigma Σ
SVD的例子
A
=
(
0
1
1
1
1
0
)
A=\begin{pmatrix} 0 & 1\\ 1 & 1 \\ 1 & 0 \\\end{pmatrix}
A=⎝⎛011110⎠⎞
首先计算
A
A
T
=
(
0
1
1
1
1
0
)
(
0
1
1
1
1
0
)
=
(
2
1
1
2
)
AA^T=\begin{pmatrix} 0 & 1 & 1\\ 1 & 1 & 0 \end{pmatrix} \begin{pmatrix} 0 & 1\\ 1 & 1 \\ 1 & 0 \\\end{pmatrix} =\begin{pmatrix} 2 & 1\\ 1 & 2 \end{pmatrix}
AAT=(011110)⎝⎛011110⎠⎞=(2112)
A
T
A
=
(
0
1
1
1
1
0
)
(
0
1
1
1
1
0
)
=
(
1
1
0
1
2
1
0
1
1
)
A^TA=\begin{pmatrix} 0 & 1\\ 1 & 1 \\ 1 & 0 \\\end{pmatrix} \begin{pmatrix} 0 & 1 & 1\\ 1 & 1 & 0 \end{pmatrix} =\begin{pmatrix} 1 & 1 & 0 \\ 1 & 2 & 1 \\ 0 & 1 & 1 \\\end{pmatrix}
ATA=⎝⎛011110⎠⎞(011110)=⎝⎛110121011⎠⎞
然后计算
A
A
T
AA^T
AAT或
A
T
A
A^TA
ATA的特征值和特征向量:
A
A
T
AA^T
AAT的特征值和特征向量:
λ
1
=
3
,
V
1
=
(
1
2
1
2
)
\lambda_1=3,V_1=\begin{pmatrix} \frac{1}{\sqrt2} \\ \frac{1}{\sqrt2} \end{pmatrix}
λ1=3,V1=(2121)
λ
2
=
1
,
V
2
=
(
−
1
2
1
2
)
\lambda_2=1,V_2=\begin{pmatrix} -\frac{1}{\sqrt2} \\ \frac{1}{\sqrt2} \end{pmatrix}
λ2=1,V2=(−2121)
A T A A^TA ATA的特征值和特征向量:
λ
1
=
3
,
U
1
=
(
1
6
2
6
1
6
)
\lambda_1=3,U_1=\begin{pmatrix} \frac{1}{\sqrt6} \\ \frac{2}{\sqrt6} \\\frac{1}{\sqrt6} \end{pmatrix}
λ1=3,U1=⎝⎜⎛616261⎠⎟⎞
λ
2
=
1
,
U
2
=
(
1
2
0
−
1
2
)
\lambda_2=1,U_2=\begin{pmatrix} \frac{1}{\sqrt2} \\ 0 \\ -\frac{1}{\sqrt2} \end{pmatrix}
λ2=1,U2=⎝⎛210−21⎠⎞
λ
2
=
1
,
U
2
=
(
1
3
−
1
3
1
3
)
\lambda_2=1,U_2=\begin{pmatrix} \frac{1}{\sqrt3} \\ -\frac{1}{\sqrt3} \\ \frac{1}{\sqrt3} \end{pmatrix}
λ2=1,U2=⎝⎜⎛31−3131⎠⎟⎞
然后,利用
A
V
i
=
σ
i
u
i
,
i
−
1
,
2
AV_i=\sigma_i u_i,i-1,2
AVi=σiui,i−1,2,求奇异值
(
0
1
1
1
1
0
)
(
1
2
1
2
)
=
σ
1
(
1
6
2
6
1
6
)
\begin{pmatrix} 0 & 1\\ 1 & 1 \\ 1 & 0 \\\end{pmatrix} \begin{pmatrix} \frac{1}{\sqrt2} \\ \frac{1}{\sqrt2} \end{pmatrix}=\sigma_1 \begin{pmatrix} \frac{1}{\sqrt6} \\ \frac{2}{\sqrt6} \\\frac{1}{\sqrt6} \end{pmatrix}
⎝⎛011110⎠⎞(2121)=σ1⎝⎜⎛616261⎠⎟⎞
解得:
σ
1
=
3
\sigma_1=\sqrt3
σ1=3
(
0
1
1
1
1
0
)
(
−
1
2
1
2
)
=
σ
1
(
1
3
−
1
3
1
3
)
\begin{pmatrix} 0 & 1\\ 1 & 1 \\ 1 & 0 \\\end{pmatrix} \begin{pmatrix} -\frac{1}{\sqrt2} \\ \frac{1}{\sqrt2} \end{pmatrix}=\sigma_1 \begin{pmatrix} \frac{1}{\sqrt3} \\ -\frac{1}{\sqrt3} \\ \frac{1}{\sqrt3} \end{pmatrix}
⎝⎛011110⎠⎞(−2121)=σ1⎝⎜⎛31−3131⎠⎟⎞
解得:
σ
2
=
1
\sigma_2=1
σ2=1
【注】这一步也可以用
σ
i
=
λ
i
\sigma_i=\sqrt{\lambda_i}
σi=λi直接计算求出奇异值
3
\sqrt 3
3和
1
1
1
最后得到A的奇异值分解为:
A
=
U
Σ
V
T
=
(
1
6
1
2
1
3
2
6
0
−
1
3
1
6
1
2
1
3
)
(
3
0
0
1
0
0
)
(
1
2
1
2
−
1
2
1
2
)
A=U \Sigma V^{T}=\begin{pmatrix} \frac{1}{\sqrt6} &\frac{1}{\sqrt2} & \frac{1}{\sqrt3} \\ \frac{2}{\sqrt6} & 0 & -\frac{1}{\sqrt3} \\ \frac{1}{\sqrt6} & \frac{1}{\sqrt2} & \frac{1}{\sqrt3} \end{pmatrix}\begin{pmatrix} \sqrt3 & 0\\ 0 & 1 \\ 0 & 0 \end{pmatrix}\begin{pmatrix} \frac{1}{\sqrt2} & \frac{1}{\sqrt2} \\ -\frac{1}{\sqrt2} & \frac{1}{\sqrt2} \end{pmatrix}
A=UΣVT=⎝⎜⎛6162612102131−3131⎠⎟⎞⎝⎛300010⎠⎞(21−212121)
5.PCA算法实现的两种方法
(1) 基于特征值分解协方差矩阵实现PCA算法
输入:数据集 x = { x 1 , x 2 , x 3 , ⋯ , x n } x=\{x_1,x_2,x_3,\cdots,x_n\} x={x1,x2,x3,⋯,xn},需要降到 k k k维
- 去平均值(即去中心化),每一维特征减去各自的平均值
- 计算协方差矩阵
1
n
X
X
T
\frac{1}{n}XX^T
n1XXT
【注】这里除 n n n或者 n − 1 n-1 n−1,对计算结果没有影响 - 用特征值分解法求协方差矩阵 1 n X X T \frac{1}{n}XX^T n1XXT的特征值和特征向量
- 对特征值从小到大排序,选择其中最大的 k k k个,然后将其对应的 k k k个特征向量分别作为行向量组成特征向量矩阵 P P P
- 将数据转换到 k k k个特征向量构建的新空间中,即 Y = P X Y=PX Y=PX
【例】
若有
X
X
X矩阵如下:请使用PCA方法将
X
X
X降为1行数据
X
=
(
−
1
−
1
0
2
0
−
2
0
0
1
1
)
X=\begin{pmatrix} -1 & -1 & 0 & 2 & 0\\ -2 & 0 & 0 & 1 & 1 \end{pmatrix}
X=(−1−2−10002101)
(1)因为
X
X
X矩阵每行已经是零均值,因此不需要再去求均值
(2)求协方差矩阵
C
=
1
5
(
−
1
−
1
0
2
0
−
2
0
0
1
1
)
(
−
1
−
1
0
2
0
−
2
0
0
1
1
)
T
=
(
6
5
4
5
4
5
6
5
)
C=\frac{1}{5}\begin{pmatrix}-1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1\end{pmatrix} \begin{pmatrix}-1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1\end{pmatrix}^T=\begin{pmatrix} \frac{6}{5} & \frac{4}{5} \\ \frac{4}{5} & \frac{6}{5} \end{pmatrix}
C=51(−1−2−10002101)(−1−2−10002101)T=(56545456)
(3)求协方差矩阵的特征值和特征向量
求解后的特征值为:
λ
1
=
2
,
λ
2
=
2
5
\lambda_1=2,\lambda_2=\frac{2}{5}
λ1=2,λ2=52
对应的特征向量为一个通解,
C
1
,
C
2
C_1,C_2
C1,C2可取任意实数,标准化后的特征向量为:
(
1
2
1
2
)
,
(
−
1
2
1
2
)
\begin{pmatrix} \frac{1}{\sqrt2} \\ \frac{1}{\sqrt2} \end{pmatrix},\begin{pmatrix} -\frac{1}{\sqrt2} \\ \frac{1}{\sqrt2} \end{pmatrix}
(2121),(−2121)
(4)矩阵
P
P
P为
P
=
(
1
2
1
2
−
1
2
1
2
)
P=\begin{pmatrix} \frac{1}{\sqrt2} & \frac{1}{\sqrt2} \\ -\frac{1}{\sqrt2} & \frac{1}{\sqrt2} \end{pmatrix}
P=(21−212121)
(5)最后,我们用
P
P
P的第一行乘以数据矩阵
X
X
X,就得到了降维后的表示:
Y
=
(
1
2
1
2
)
(
−
1
−
1
0
2
0
−
2
0
0
1
1
)
=
(
−
3
2
−
1
2
0
3
2
1
2
)
Y=\begin{pmatrix} \frac{1}{\sqrt2} \frac{1}{\sqrt2} \end{pmatrix} \begin{pmatrix} -1 & -1 & 0 & 2 & 0\\ -2 & 0 & 0 & 1 & 1 \end{pmatrix}=\begin{pmatrix} -\frac{3}{\sqrt2} & -\frac{1}{\sqrt2} & 0 & \frac{3}{\sqrt2} & \frac{1}{\sqrt2} \end{pmatrix}
Y=(2121)(−1−2−10002101)=(−23−2102321)
注意:通过特征值分解协方差矩阵,那么我们只能得到一个方向的PCA降维,这个方向就是对数据矩阵X从行(或列)方向上的压缩存储
(2)基于SVD分解协方差矩阵实现PCA算法
输入:数据集 x = { x 1 , x 2 , x 3 , ⋯ , x n } x=\{x_1,x_2,x_3,\cdots,x_n\} x={x1,x2,x3,⋯,xn},需要降到 k k k维
- 去平均值,每一维特征减去各自的平均值
- 计算协方差矩阵
- 通过SVD计算协方差矩阵的特征值和特征向量
- 对特征值从大到小排序,选取其中最大的K个,然后将其对应的K个特征向量分别作为列向量组成特征向量矩阵
- 将数据转换到K个特征向量构建的新空间中
【注】
(1)当样本量很大的时候,往往采用SVD实现PCA算法,因为SVD可以不计算特征值,大量的样本的特征值的计算量往往很大
(2)右奇异矩阵的用处:
假设样本是 m × n m×n m×n的矩阵X,通过SVD找到矩阵 X T X X^TX XTX最大的K个特征值组成的矩阵 V T V^T VT,则我们做如下处理:
X m × n ′ = X m × n V n × k T > > > X_{m×n}^{'}=X_{m×n}V_{n×k}^{T} >>> Xm×n′=Xm×nVn×kT>>>列数从n减到k即进行了列压缩。
左奇异矩阵用于对行数的压缩,右奇异矩阵用于对列数的压缩
6.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))
上面代码实现了对数据X进行特征的降维。结果如下:
参考博文:https://blog.csdn.net/program_developer/article/details/80632779