文章目录
说明:下文中,粗斜体字母均表示矩阵(如 A \boldsymbol A A);为不引起歧义,列向量也均加箭头表示(如 a ⃗ \vec a a)
概述
PCA的目的
假设现在有这样一个情景:现在要统计并可视化分析男大学生体测成绩,如果只参考立定跳远和1000m成绩两项指标,我们可以以立定跳远成绩作为 x x x轴,1000m成绩作为 y y y轴做出散点图,每个点代表一个学生;若统计三项指标,我们也可以在三维空间中做出散点图;但如果要统计四项乃至更多的指标,我们就无法再以此方法进行数据的可视化。
而主成分分析(Principal Component Analysis,PCA)的方法,可以将具有多个观测变量的高维数据集降维,使人们可以从事物之间错综复杂的关系中找出一些主要的方面,从而能更加有效地利用大量统计数据进行定量分析,并可以更好地进行可视化、回归等后续处理。
PCA的几何意义
先将问题简化为二维情形。有 N N N个样品,具有两个观测变量 X 1 , X 2 X_1,\ X_2 X1, X2,做出散点图(如下图中的蓝色点),这样,在由 X 1 , X 2 X_1,\ X_2 X1, X2组成的坐标空间中, N N N个样品的分布情况如带状。现在问:如果现在要将两个观测变量缩减为一个,应该如何选取?
可以直观地看出,这 N N N个样品无论沿 X 1 X_1 X1轴还是沿 X 2 X_2 X2轴方向,均有较大的离散性(其离散程度可以分别用观测变量 X 1 X_1 X1的方差和 X 2 X_2 X2的方差定量表示),也就是说,只考虑 X 1 , X 2 X_1,\ X_2 X1, X2的其中一个,原始数据均会有较大损失。
现在考虑以下线性组合,变换坐标空间:
{
T
1
=
X
1
cos
θ
+
X
2
sin
θ
T
2
=
−
X
1
sin
θ
+
X
2
cos
θ
\begin{cases} T_1 = X_1 \cos \theta + X_2 \sin \theta \\ T_2 = -X_1 \sin \theta + X_2 \cos \theta \end{cases}
{T1=X1cosθ+X2sinθT2=−X1sinθ+X2cosθ
即
[ T 1 T 2 ] = [ X 1 X 2 ] [ cos θ − sin θ sin θ cos θ ] = [ X 1 X 2 ] W (1) \begin{bmatrix} T_1 & T_2 \end{bmatrix}=\begin{bmatrix} X_1 & X_2 \end{bmatrix} \begin{bmatrix} \cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{bmatrix} = \begin{bmatrix} X_1 & X_2 \end{bmatrix} \boldsymbol W \tag{1} [T1T2]=[X1X2][cosθsinθ−sinθcosθ]=[X1X2]W(1)
式中 W \boldsymbol W W为旋转变换矩阵,有如下性质:
- W \boldsymbol W W的第 i i i列组成的列向量(也就是后面将要提到的特征向量),代表的就是新坐标空间的基底 T i T_i Ti在原坐标空间中的坐标,且其为单位向量;
- W \boldsymbol W W为正交阵,即满足: W T W = I \boldsymbol W^{\rm T} \boldsymbol W=\boldsymbol I WTW=I。
经过旋转, N N N个数据点在 T 1 T_1 T1轴上的离散程度最大,因而变量 T 1 T_1 T1代表了原始数据的绝大部分信息,这样,即使不考虑变量 T 2 T_2 T2也不会损失太多数据信息。这个 T 1 T_1 T1即为第一主成分(Principal Component 1,PC1),如图中箭头所示。若将所有数据点投影到 T 1 T_1 T1轴上(图中橙色点),就得到了降维后的数据。若有多个主成分,则:
- 这些主成分之间相互独立,即没有重叠的信息,亦即这些特征向量之间正交, c o v ( T i , T j ) = 0 , i ≠ j {\rm cov} \left( T_i, T_j \right) = 0,\quad i \ne j cov(Ti,Tj)=0,i=j;
- 主成分的方差依次递减, v a r ( T 1 ) ≥ v a r ( T 2 ) ≥ ⋯ {\rm var}\left( T_1 \right) \ge {\rm var}\left( T_2 \right) \ge \cdots var(T1)≥var(T2)≥⋯
也就是说,PCA并不会对原有数据做任何的改变,而只是将“观看”原有数据的视角转换了,即,在原有数据空间中的数据的相对位置,与在主成分空间(Principal Component Space)中的相对位置是完全相同的,相当于只是更换了原有数据的基底。
原理与步骤简述
算法一:特征分解(Eigen Decomposition)
假设有一 n × m n\times m n×m维的数据矩阵 A = [ a ⃗ 1 T a ⃗ 2 T ⋮ a ⃗ n T ] \boldsymbol A = \begin{bmatrix} \vec a_1^{\rm T} \\ \vec a_2^{\rm T} \\ \vdots \\ \vec a_n^{\rm T} \end{bmatrix} A=⎣⎢⎢⎢⎡a1Ta2T⋮anT⎦⎥⎥⎥⎤,其中 n n n为样本量, m m m为观测变量的数量。PCA的步骤如下:
-
先对 A \boldsymbol A A进行中心化(整体平移数据,使数据中心在 ( 0 , 0 ) (0,0) (0,0)):
- 对 A \boldsymbol A A求列上的平均值: a ⃗ T ‾ = 1 n ∑ i = 1 n a ⃗ i T \overline {\vec a^{\rm T}} = \dfrac 1 n \sum _{i=1}^n \vec a_i ^{\rm T} aT=n1∑i=1naiT(结果为一行向量);
- 记 A ˉ = [ 1 1 ⋮ 1 ] a ⃗ T ‾ = [ a ⃗ T ‾ a ⃗ T ‾ ⋮ a ⃗ T ‾ ] \boldsymbol {\bar A} = \begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{bmatrix}\overline {\vec a^{\rm T}} = \begin{bmatrix} \overline {\vec a^{\rm T}} \\ \overline {\vec a^{\rm T}} \\ \vdots \\ \overline {\vec a^{\rm T}} \end{bmatrix} Aˉ=⎣⎢⎢⎢⎡11⋮1⎦⎥⎥⎥⎤aT=⎣⎢⎢⎢⎡aTaT⋮aT⎦⎥⎥⎥⎤;
- 中心化后的数据矩阵 X = A − A ˉ \boldsymbol X = \boldsymbol A - \boldsymbol {\bar A} X=A−Aˉ( n × m n\times m n×m维)。
-
计算 X \boldsymbol X X的协方差矩阵 C \boldsymbol C C:
C = X T X (2) \boldsymbol C = \boldsymbol X^{\rm T}\boldsymbol X \tag{2} C=XTX(2) -
对 X \boldsymbol X X做特征分解,即求解特征方程
∣ C − λ I ∣ = 0 (3) \left| \boldsymbol C - \lambda \boldsymbol I\right|=0 \tag{3} ∣C−λI∣=0(3)
可得到 m m m个特征值(Eigenvalues) λ i \lambda_i λi。再解方程
( C − λ i I ) w ⃗ i = 0 (4) \left( \boldsymbol C - \lambda _i \boldsymbol I\right)\vec w_i =0 \tag{4} (C−λiI)wi=0(4)
其中 w ⃗ i = [ w 1 i w 2 i ⋮ w m i ] , i = 1 , 2 , ⋯ , m \vec w_i= \begin{bmatrix} w_{1i} \\ w_{2i} \\ \vdots \\ w_{mi}\end{bmatrix},\quad i=1,2,\cdots ,m wi=⎣⎢⎢⎢⎡w1iw2i⋮wmi⎦⎥⎥⎥⎤,i=1,2,⋯,m,得到 m m m个特征向量(Eigenvectors) w ⃗ i \vec w_i wi,将它们组成矩阵 W \boldsymbol W W。可以验证, ∑ j = 1 m w j i = 1 \sum_{j=1}^m w_{ji}=1 ∑j=1mwji=1。 -
将特征值降序排列,其对应的特征向量也排列到对应位置(调换 W \boldsymbol W W的列)。我们这样做的原因是, v a r ( T i ) = λ i {\rm var} \left( T_i\right)=\lambda_i var(Ti)=λi。
进行特征还原:
T = X W (5) \boldsymbol T = \boldsymbol X \boldsymbol W \tag{5} T=XW(5)
其中:- T \boldsymbol T T, n × m n\times m n×m维,称为主成分得分(principal component scores),即为新坐标空间中的数据点
- W \boldsymbol W W, m × m m\times m m×m维,为特征向量组成的矩阵(称为loadings)
-
我们可以只取 W \boldsymbol W W的前 r r r列,即将 m × m m\times m m×m维矩阵缩减为 m × r m\times r m×r维矩阵,记作 W r \boldsymbol W_r Wr,则有:
T r = X W r (6) \boldsymbol T_r = \boldsymbol X \boldsymbol W_r \tag{6} Tr=XWr(6)
T r \boldsymbol T_r Tr同样为 T \boldsymbol T T从 n × m n\times m n×m维缩减为 n × r n\times r n×r维的结果,相当于将原有的 m m m个观测变量缩减为最主要的 r r r个,即达到了我们的目的——降维。
算法二:奇异值分解(Singular Value Decomposition,SVD)
对矩阵
X
\boldsymbol X
X进行奇异值分解:
X
=
U
Σ
V
∗
=
U
Σ
V
T
\boldsymbol X = \boldsymbol U \boldsymbol \Sigma \boldsymbol V^* = \boldsymbol U \boldsymbol \Sigma \boldsymbol V^{\rm T}
X=UΣV∗=UΣVT
其中:
- U \boldsymbol U U( n × n n\times n n×n 维), V \boldsymbol V V( m × m m\times m m×m 维)分别称为左奇异向量和右奇异向量;
- V ∗ \boldsymbol V^* V∗表示矩阵 V \boldsymbol V V的共轭转置矩阵,因为 X \boldsymbol X X为实数矩阵,所以可写为 V ∗ = V T \boldsymbol V^*=\boldsymbol V^{\rm T} V∗=VT;
- Σ \boldsymbol \Sigma Σ为一矩形对角矩阵( n × m n\times m n×m维),其对角线元素称为奇异值(singular value)。
在PCA的问题中恒有: W = V \boldsymbol W = \boldsymbol V W=V。与 W \boldsymbol W W相似, V \boldsymbol V V具有以下性质:
- Σ \boldsymbol \Sigma Σ的对角线元素也满足 σ 1 > σ 2 > ⋯ \sigma_1> \sigma_2 > \cdots σ1>σ2>⋯,即,也是降序排列的;
- U \boldsymbol U U和 V \boldsymbol V V的第 i i i列,对应于 Σ \boldsymbol \Sigma Σ的第 i i i个元素(第 i i i大元素) σ i \sigma_i σi;
- U \boldsymbol U U和 V \boldsymbol V V满足: U ∗ U = I \boldsymbol U^* \boldsymbol U = \boldsymbol I U∗U=I, V ∗ V = I \boldsymbol V^* \boldsymbol V = \boldsymbol I V∗V=I。
故,公式 ( 5 ) (5) (5)可写为:
T = X W = X V = U Σ V ∗ V = U Σ \begin{aligned} \boldsymbol T & = \boldsymbol X \boldsymbol W\\ & = \boldsymbol X \boldsymbol V\\ & = \boldsymbol U \boldsymbol \Sigma \boldsymbol V^*\boldsymbol V \\ & = \boldsymbol U \boldsymbol \Sigma \\ \end{aligned} T=XW=XV=UΣV∗V=UΣ
即:
T = U Σ (7) \boldsymbol T = \boldsymbol U \boldsymbol \Sigma \tag{7} T=UΣ(7)
我们同样可以对 Σ \boldsymbol \Sigma Σ只取第一个 r × r r \times r r×r的块,记作 Σ r \boldsymbol \Sigma _ r Σr,相应地 U \boldsymbol U U也只取前 r r r列,记作 U r \boldsymbol U_r Ur,则有
T r = U r Σ r (8) \boldsymbol T_r = \boldsymbol U_r \boldsymbol \Sigma _r \tag{8} Tr=UrΣr(8)
r r r的选取标准
计算方差的累积贡献率:
f
(
k
)
=
∑
i
=
1
i
λ
k
∑
i
=
1
m
λ
i
,
k
=
1
,
2
,
⋯
(9)
f(k)=\dfrac{\sum _{i=1}^i \lambda_k}{\sum_{i=1}^m \lambda_i},\quad k = 1,2,\cdots \tag{9}
f(k)=∑i=1mλi∑i=1iλk,k=1,2,⋯(9)
作出其图像。因为
λ
1
>
λ
2
>
⋯
\lambda_1> \lambda_2 > \cdots
λ1>λ2>⋯,故
f
(
k
)
f(k)
f(k)为一单调递增的函数,且其递增速度随着
k
k
k增加逐渐降低。
一般地,我们可以取使得 f ( r ) ≥ f(r) \ge f(r)≥某一阈值(如 95 % 95\% 95%)的最小的 r r r,这样最多只会损失掉5%的方差。
对于SVD法,将公式 ( 9 ) (9) (9)中的 λ \lambda λ换为 σ \sigma σ即可。
两种算法的比较
在采用特征分解法时,我们无法避免计算 X T X \boldsymbol X^{\rm T}\boldsymbol X XTX,而在观测变量数 m m m非常大时,这一算法的劣势将被无限放大(协方差矩阵为 m × m m\times m m×m维)。
而采用SVD算法,则只需要计算
T
r
=
U
r
Σ
r
\boldsymbol T_r = \boldsymbol U_r \boldsymbol \Sigma _r
Tr=UrΣr,而
Σ
r
\boldsymbol \Sigma _r
Σr为对角阵,显然这一算法的计算量要小很多(这一点类似于DFT与FFT之间的比较)。默认情况下,MATLAB中的pca
函数也会使用SVD算法。
MATLAB的实现方法
我们先载入MATLAB自带的数据集fisheriris
(该数据集统计了三种鸢尾花的花萼长、花萼宽、花瓣长、花瓣宽),然后进行中心化处理,并计算协方差矩阵:
load fisheriris;
X = meas; % n = 150, m = 4
% 中心化
meanX = ones(size(X,1), 1) * mean(X);
centredX = X - meanX;
C = cov(centredX); % 直接调用cov直接计算协方差矩阵即可
特征分解法:利用eig
函数
[W, Lambda] = eig(C); % W是特征向量组成的矩阵(4×4),Lambda是特征值组成的对角矩阵
ev = (diag(Lambda))'; % 提取特征值
ev = ev(:, end:-1:1); % eig计算出的特征值是升序的,这里手动倒序(W同理)
W = W(:, end:-1:1);
sum(W.*W, 1) % 可以验证每个特征向量各元素的平方和均为1
Wr = W(:, 1:2); % 提取前两个主成分的特征向量
Tr = centredX * Wr; % 新坐标空间的数据点
% 作图
figure;
stairs(cumsum(ev)/sum(ev), 'LineWidth',1.5);
axis([1 4 0 1]);
xlabel('$ k $', 'Interpreter', 'latex');
ylabel('$ f(k)=\frac{\sum _{i=1}^i \lambda_k}{\sum_{i=1}^m \lambda_i} $',...
'Interpreter', 'latex');
hold on;
plot([1 4], [0.95 0.95], '--'); % 从图中可以看出,取r = 2即可
figure;
scatter(Tr(:,1), Tr(:,2), 130, categorical(species), '.');
colormap(winter);
xlabel('Principal Component 1');
ylabel('Principal Component 2');
SVD法:利用svd
函数
[U, Sigma, V] = svd(X); % 可以检验,W和V完全相同(向量的正负号不影响)
Vr = V(:, 1:2); % 提取前两个主成分的特征向量
Tr = X * Vr; % 新坐标空间的数据点
% 画图部分同上,略
利用pca
函数
pca
的常用调用格式如下:
[loadings, scores] = pca(X, 'NumComponents', r);
其中:
loadings
为 W r \boldsymbol W_r Wr矩阵( m × r m\times r m×r维),即主成分系数;scores
为 T r \boldsymbol T_r Tr矩阵( n × r n\times r n×r维);eigenvalues
为所有特征值组成的列向量。
且默认情况下,pca
会自动将数据X
中心化。
在本例中,我们可以略去中心化的步骤,直接调用该函数:
[Wr, Tr, ev] = pca(X, 'NumComponents',2);
% 画图部分略
应用
聚类分析
如上面鸢尾花的例子中,降维后的数据仍可以清晰地分为三类。这样,当我们拿到一种鸢尾花,计算相应的 T 1 T_1 T1和 T 2 T_2 T2,将结果画在散点图中,我们就可以判断出其属于哪一种鸢尾花。
例如,我们在电商平台浏览并购买商品时,平台就会收集你的年龄、性别、购买商品平均价格、购买频率、最初浏览商品直到最终购买之间的时间间隔等等大量、多维度的信息,然后进行降维,将你归于“大学生”“白领”“一家三口”等类别,然后定向为你推送促销商品的通知。
图像压缩
假设有一张 n × m n \times m n×m的图片 X \boldsymbol X X,根据 ( 6 ) (6) (6)式,我们可以利用矩阵 W r \boldsymbol W_r Wr将其降至 n × r n\times r n×r维。如果我们只传输 T r \boldsymbol T_r Tr和 W r \boldsymbol W_r Wr,就可以反推还原出 X \boldsymbol X X,而压缩比可达 n m r ( n + m ) ≤ n m 2 r \dfrac {nm}{r(n+m)} \le \dfrac {\sqrt{nm}} {2r} r(n+m)nm≤2rnm
人脸检测与匹配
假设有 n n n个人脸训练样本,每个样本共 m m m个像素,每个样本由其像素灰度值展开组成一个行向量,按列组成矩阵 X \boldsymbol X X。
同样先求其平均向量(称为“平均脸”),中心化后求协方差矩阵,并进行特征分解。任何一幅人脸图像都可以变换到主成分空间,得到“特征脸”(Eigenfaces)。
这样,若将待识别的人脸做同样的变换,遍历已有的特征脸中,寻找最为接近的特征脸,即完成了匹配。