Principal Component Analysis 主成分分析
1.PCA背景
Figure 1.
Principal Component Analysis(PCA)是一种数据降维方式,可以将n维数据降低到指定的k维。此时延伸出一个问题,为什么数据要进行降维操作。如 Figure 1 在我们常用的二维坐标系中,我们可以很容易的观察出二维数据的分布规律,从而判断出数据的特性。假设在此数据的基础上,加上z轴,每个数据点在z平面上的抖动非常小,则数据点在z轴的变化对于数据整体的特性影响很小。我们足以在x、y轴上观察到数据的主要规律,那么数据在z轴上的维度可以选择忽略。数据降维有以下优点:
- 降低计算成本:高维数据集通常需要更多计算资源来处理和分析,通过降低数据维度,可以减少计算成本和存储需求
- 可视化:高维数据难以直观理解和可视化,降维可以将数据投影到低维空间中,使得数据可以更容易地可视化和理解
- 特征提取:降维可以帮助识别和提取数据中最重要的特征,从而简化数据分析和模型训练过程
- 去除冗余信息:高维数据中可能存在大量冗余信息,降维可以去除这些冗余信息,提高数据的表征效率
- 减少相关性:在数据中存在高度相关的特征时,降维可以减少特征之间的相关性,提高模型的稳定性和可解释性
由此可以看出,数据降维是有很大意义的。PCA作为一种简单而有效的降维方法,具有线性变换、最大化方差、去相关性、数据可视化和保留数据结构等优势
2.PCA降维过程
Figure 2.
以二维数据为例,首先将原始数据中心化(如Figure 2.),是数据每个维度的均值为0,使其分布在原点周围。中心化可以消除数据的平移性影响,使得数据的中心位于原点附近,助于更好地理解数据的分布和结构。对数据
X
∈
R
n
×
2
\mathbf{X}\in\mathbb{R}^{n\times 2}
X∈Rn×2 的维度(x、y)之间计算协方差得到协方差矩阵
C
\mathbf{C}
C:
C
=
(
c
o
v
(
x
,
x
)
c
o
v
(
x
,
y
)
c
o
v
(
y
,
x
)
c
o
v
(
y
,
y
)
)
(1)
\mathbf{C}=\left(\begin{array}{cc}cov(x,x)&cov(x,y)\\cov(y,x)&cov(y,y)\end{array}\right) \tag{1}
C=(cov(x,x)cov(y,x)cov(x,y)cov(y,y))(1)
计算得到特征值(
λ
1
,
λ
2
\lambda_1,\lambda_2
λ1,λ2)与对应的特征向量(
ξ
1
,
ξ
2
\xi_1,\xi_2
ξ1,ξ2)。将特征值从大到小排列,取前k个特征值与相应的特征向量计算主要成分,以上图为例,k取1,我们得到了u1轴对应的
λ
,
ξ
\lambda,\xi
λ,ξ 。降维后的数据
X
′
∈
R
n
×
1
\mathbf{X}^{'}\in\mathbb{R}^{n\times 1}
X′∈Rn×1:
X
′
=
X
×
ξ
T
(2)
\mathbf{X}^{'}=\mathbf{X}\times \xi^{\mathbf{T}} \tag{2}
X′=X×ξT(2)
为了在二维空间中表示出来,我们将
X
′
\mathbf{X}^{'}
X′ 可视化为中心化数据在
ξ
\xi
ξ 上的投影点,如Figure 3.,向量投影过程:设
α
=
(
α
x
,
α
y
)
,
β
=
(
β
x
,
β
y
)
\alpha=(\alpha_x,\alpha_y),\beta=(\beta_x,\beta_y)
α=(αx,αy),β=(βx,βy),则
α
\alpha
α 在
β
\beta
β 上的投影表示为:
α
p
r
o
=
(
α
⋅
β
∣
∣
β
∣
∣
2
)
β
(3)
\alpha_{pro}=(\frac{\alpha\cdot\beta}{||\beta||^2})\beta \tag{3}
αpro=(∣∣β∣∣2α⋅β)β(3)
Figure 3.
已知我们选取的特征值是最大的特征值, C \mathbf{C} C 有两个特征值,如果取小的特征值对应的特征向量(u2轴)效果会什么样呢,如Figure 4.
Figure 4.
综上,我们发现在较大特征值的特征向量上的投影上的数据点要比较小的特征值对应的特征向量所得到的投影数据方差更大,更能表现出数据的特征规律,故应选取较大特征值对应的特征向量
3.PCA数学推理
给定一组数据
{
z
⃗
1
,
z
⃗
2
,
.
.
.
,
z
⃗
n
}
\{\vec{z}_1,\vec{z}_2,...,\vec{z}_n\}
{z1,z2,...,zn} ,数据中心位于
μ
⃗
=
1
n
∑
i
=
1
n
z
⃗
i
\vec{\mu}=\frac1n\sum_{i=1}^n\vec{z}_i
μ=n1∑i=1nzi ,中心化处理后数据为
{
x
⃗
1
,
x
⃗
2
,
.
.
.
x
⃗
n
}
=
{
z
⃗
1
−
μ
⃗
,
z
⃗
2
−
μ
⃗
,
.
.
.
,
z
⃗
n
−
μ
⃗
}
\{\vec{x}_1,\vec{x}_2,...\vec{x}_n\}=\{\vec{z}_1-\vec{\mu},\vec{z}_2-\vec{\mu},...,\vec{z}_n-\vec{\mu}\}
{x1,x2,...xn}={z1−μ,z2−μ,...,zn−μ} 。 中心化后的数据在第一主轴u1方向上分布散的最开,也就是说在u1方向上的投影的绝对值之和最大(方差最大),将x与u1做内积作为计算投影的方法,由于只需要求u1的方向,所以设u1也是单位向量,则转化为最大化
1
n
∑
i
=
1
n
∣
x
⃗
i
⋅
u
⃗
1
∣
→
1
n
∑
i
=
1
n
(
x
⃗
i
⋅
u
⃗
1
)
2
\frac1n{\sum_{i=1}^n}|\vec{x}_i\cdot\vec{u}_1|\to \frac1n{\sum_{i=1}^n(\vec{x}_i\cdot\vec{u}_1)}^2
n1∑i=1n∣xi⋅u1∣→n1∑i=1n(xi⋅u1)2 ,向量的乘法转化为内积
x
⃗
i
⋅
u
⃗
1
=
x
i
T
u
1
\vec{x}_i\cdot\vec{u}_1=x_i^Tu_1
xi⋅u1=xiTu1 ,故目标函数表示为
1
n
∑
i
=
1
n
(
x
i
T
u
1
)
2
\frac1n{\sum_{i=1}^n(x_i^Tu_1)^2}
n1∑i=1n(xiTu1)2 ,由于列向量转置以后是行向量,行向量乘以列向量得到一个数,一个数的转置还是其本身,所以又可以将目标函数化为:
1
n
∑
i
=
1
n
(
x
i
T
u
1
)
T
(
x
i
T
u
1
)
=
1
n
∑
i
=
1
n
u
1
T
x
i
x
i
T
u
1
(4)
\frac1n\sum_{i=1}^n(x_i^Tu_1)^T(x_i^Tu_1)=\frac1n\sum_{i=1}^nu_1^Tx_ix_i^Tu_1 \tag{4}
n1i=1∑n(xiTu1)T(xiTu1)=n1i=1∑nu1TxixiTu1(4)
令
X
=
{
x
1
,
x
2
,
.
.
.
,
x
n
}
X=\{x_1,x_2,...,x_n\}
X={x1,x2,...,xn} ,得到
X
X
T
=
∑
i
=
1
n
x
i
x
i
T
XX^T=\sum_{i=1}^nx_ix_i^T
XXT=∑i=1nxixiT ,目标函数转化为:
1
n
u
1
T
X
X
T
u
1
(5)
\frac1nu_1^TXX^Tu_1 \tag{5}
n1u1TXXTu1(5)
其中
u
1
T
X
X
T
u
1
u_1^TXX^Tu_1
u1TXXTu1 是二次型,设
X
X
T
XX^T
XXT 的某个特征值为
λ
\lambda
λ ,对应的特征向量为
ξ
\xi
ξ ,则有:
(
X
X
T
)
T
=
X
X
T
X
X
T
ξ
=
λ
ξ
(
X
X
T
ξ
)
T
ξ
=
(
λ
ξ
)
T
ξ
ξ
T
X
X
T
ξ
=
λ
ξ
T
ξ
ξ
T
X
X
T
ξ
=
(
X
T
ξ
)
T
(
X
T
ξ
)
=
∥
X
T
ξ
∥
2
=
λ
ξ
T
ξ
=
λ
∥
ξ
∥
2
≥
0
(6)
(XX^T)^T=XX^T \\ XX^T\xi=\lambda\xi \\ \left(XX^{T}\xi\right)^{T}\xi=\left(\lambda\xi\right)^{T}\xi \\ \xi^TXX^T\xi=\lambda\xi^T\xi \\ \xi^{T}XX^{T}\xi=\left(X^{T}\xi\right)^{T}\left(X^{T}\xi\right)=\left\|X^{T}\xi\right\|^{2}=\lambda\xi^{T}\xi=\lambda\|\xi\|^{2}\geq0 \tag{6}
(XXT)T=XXTXXTξ=λξ(XXTξ)Tξ=(λξ)TξξTXXTξ=λξTξξTXXTξ=(XTξ)T(XTξ)=
XTξ
2=λξTξ=λ∥ξ∥2≥0(6)
故
λ
≥
0
,
X
X
T
\lambda\geq0,XX^T
λ≥0,XXT 为半正定的对称矩阵,即
ξ
T
X
X
T
ξ
\xi^{T}XX^{T}\xi
ξTXXTξ 为半正定阵的二次型,目标函数存在最大值。由二范数
∥
x
∥
2
2
=
<
x
,
x
>
=
x
T
x
\left\|x\right\|_2^2=<x,x>=x^Tx
∥x∥22=<x,x>=xTx 我们求解目标函数最大值:
u
1
T
X
X
T
u
1
=
(
X
T
u
1
)
T
(
X
T
u
1
)
=
<
X
T
u
1
,
X
T
u
1
>
=
∥
X
T
u
1
∥
2
2
=
(
∥
X
T
u
1
∥
2
∥
u
1
∥
2
)
2
(7)
u_1^TXX^Tu_1=(X^Tu_1)^T(X^Tu_1)=<X^Tu_1,X^Tu_1>=\left\|X^Tu_1\right\|_2^2=(\frac{\left\|X^Tu_1\right\|_2}{\left\|u_1\right\|_2})^2 \tag{7}
u1TXXTu1=(XTu1)T(XTu1)=<XTu1,XTu1>=
XTu1
22=(∥u1∥2
XTu1
2)2(7)
把二次型化成一个范数的形式,由于u1取单位向量,最大化目标函数的基本问题也就转化为:对一个矩阵,它对一个向量做变换,变换前后的向量的模长伸缩尺度如何才能最大?我们有矩阵代数中的定理知,向量经矩阵映射前后的向量长度之比的最大值就是这个矩阵的最大奇异值,即:
∥
X
T
u
1
∥
2
∥
u
1
∥
2
≤
σ
1
(
X
T
)
(8)
\frac{\left\|X^Tu_1\right\|_2}{\left\|u_1\right\|_2}\leq\sigma_1(X^T) \tag{8}
∥u1∥2
XTu1
2≤σ1(XT)(8)
其中,
σ
1
(
⋅
)
\sigma_1(\cdot)
σ1(⋅) 表示矩阵的最大奇异值,等价于
X
X
T
,
X
T
X
XX^T,X^TX
XXT,XTX 的最大特征值开方,且
X
X
T
,
X
T
X
XX^T,X^TX
XXT,XTX 的各自的特征值都是大于等于0的,各自的特征向量之间都是正交的。当目标函数取最大值时,则符合以下:
u
1
T
X
X
T
u
1
=
(
∥
X
T
u
1
∥
2
∥
u
1
∥
2
)
2
=
σ
1
(
X
T
)
2
=
σ
1
(
X
X
T
)
(9)
u_1^TXX^Tu_1=(\frac{\left\|X^Tu_1\right\|_2}{\left\|u_1\right\|_2})^2 =\sigma_1(X^T)^2=\sigma_1(XX^T) \tag{9}
u1TXXTu1=(∥u1∥2
XTu1
2)2=σ1(XT)2=σ1(XXT)(9)
故,此时
u
1
u_1
u1 应取
σ
1
(
X
X
T
)
\sigma_1(XX^T)
σ1(XXT) 对应的的特征向量
综上,证明了数据降维的主成分是由协方差矩阵的由大到小的特征值对应的特征向量所得到的
4.PLT绘图代码
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
matplotlib.use('Agg')
# 示例二维数据
np.random.seed(0)
X = np.random.multivariate_normal([2, 3], [[1, 0.9], [0.9, 1]], 100)
# 计算PCA
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X)
# 获取主成分向量
pc1 = pca.components_[0]
pc2 = pca.components_[1]
pc=pc2
# 计算数据的均值
mean = np.mean(X, axis=0)
# 定义沿主成分方向的直线
line_points = np.array([mean + pc1 * t for t in np.linspace(-3, 3, 100)])
# 将原始数据投影到主成分上
projections = np.dot(X - mean, pc) # 计算投影长度
# 从主成分再投影回二维空间用于可视化
projected_points = np.outer(projections, pc)
# 绘制原始数据、主成分线和投影点
plt.scatter(X[:, 0]-mean[0], X[:, 1]-mean[1], label='centered data')
plt.plot([pc[0]*4*(pc[1]/pc[0]),2], [pc[1]*4*(pc[1]/pc[0]),2*(pc[1]/pc[0])], color='red', label='eigenvector')
plt.scatter(projected_points[:, 0], projected_points[:, 1], color='green', label='pro-dot')
for i in range(X.shape[0]):
plt.plot([X[i, 0]-mean[0], projected_points[i, 0]], [X[i, 1]-mean[1], projected_points[i, 1]], 'k--', linewidth=0.5)
plt.legend()
plt.xlabel('X')
plt.ylabel('Y')
plt.title('PCA')
plt.savefig('./pca_plot.png')
plt.show()