前言
本文主要对PCA主成分分析中的数学原理进行介绍,将不涉及或很少涉及代码实现或应用,阅读前请确保已了解基本的机器学习相关内容和数学基础。
PCA的原理涉及大量线性代数知识,若没有了解或系统学习过线性代数,则不建议直接阅读此文章。
使用PCA算法固然简单,但若想合理高效地使用,了解其原理则是必不可少的一环,当前很少有博客对基本的PCA各种原理做很详细地介绍与解读,本文旨在一步步由深入浅、全方位地介绍其相关数学原理。
文章概述
PCA 主成分分析属于矩阵分解算法中的入门算法,通过分解特征矩阵来实现降维。虽是入门算法但并不代表着它的数学原理容易。
本文主要内容:
- 前置知识:样本方差、协方差、协方差矩阵、散度矩阵的简单介绍
- 特征值分解 (EVD) 和奇异值分解 (SVD) 的原理和流程
- PCA 的基本实现原理:分别基于 EVD 和 SVD 的 PCA 实现方法
- PCA的 应用以及对一些应用或说明的补充
前置知识
样本方差(Variance):样本方差反映一组数据变异程度或分散程度大小的指标;计算公式:
Var=1N−1∑i=1N(xi−xˉ)2Var = frac{1}{N-1}sum^N_{i=1}(x_i-bar{x})^2
- NN:某特征下的样本量
- xix_i:某特征下的第 ii 个样本值
- xˉbar{x}:某特征下样本的均值
为什么分母为 N−1N-1:这是在计算样本方差时对方差计算公式的一个修复,目的是得到总体方差的一个无偏估计,具体可参考下面文章中的”无偏估计与样本方差小节”:
协方差(Covariance):用来刻画两个随机变量 XX,YY之间的相关性;计算公式:
Cov(X,Y)=1N−1∑i=1N(xi−xˉ)(yi−yˉ)Cov(X,Y) = frac{1}{N-1}sum^N_{i=1}(x_i-bar{x})(y_i-bar{y})
- NN:特征 X,YX, Y 下的样本量
- xi,yix_i, y_i:特征 X,YX,Y 下的第 ii 个样本值
- xˉ,yˉbar{x},bar{y}:特征 X,YX, Y 下样本的均值
样本方差就是协方差的一种特殊形式,当两个变量相同时,协方差与样本方差相同。
如果两个变量的变化趋势一致,也就是说如果其中一个大于自身的期望值,另外一个也大于自身的期望值,那么两个变量之间的协方差就是正值,反之则为负值。
- 协方差为正或为负描述的是两组数据中的整体变化趋势相同或相反;协方差为正,则我们说两变量正相关,反之则为负相关。
协方差的绝对值越大,一定程度上两变量相关性越强
- 从公式上可以看出,若当 xix_i 大于 xˉbar{x} 时,yiy_i 总是大于或小于 yˉbar{y},则结果的绝对值是最大的(全为正数或全为负数然后累加),这种情况也就是变量 XX 与变量 YY 之间的变化趋势总是相同的,我们就认为这两个变量之间有着很强的相关性。
即协方差的正负表示两组数据的总体变化趋势,绝对值大小表示变化趋势符合程度。
协方差矩阵(Covariance Matrix):一个大小为 N∗NN*N 的对称矩阵,且是半正定矩阵,由某组数据中的每对变量之间的协方差组成,正对角线上元素为各变量的方差。
例:某数据中含有三个变量(特征),分别为 X,Y,ZX,Y,Z,则该数据的协方差矩阵为:
Σ=[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)]Sigma = begin{bmatrix} 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{bmatrix}
不难看出,Cov(X,X)Cov(X,X) 其实就是变量 XX 的方差,即正对角线上的三个元素分别为这三个变量本身的方差;由于 Cov(X,Y)=Cov(Y,X)Cov(X,Y) = Cov(Y,X),因此该矩阵为对称矩阵;共含有3个变量,因此该协方差矩阵的大小为 3∗33*3。
散度矩阵:散度矩阵和协方差矩阵类似,有时计算散度矩阵,有时计算协方差矩阵,二者意义差不多,散度矩阵公式:
S=(n−1)[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)]S = (n-1) begin{bmatrix} 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{bmatrix}
从公式也可以看出,散度矩阵其实就是协方差矩阵乘以 n−1n-1,不难看出两矩阵的特征值是相同的。
特征值分解与奇异值分解
EVD:特征值分解
对于矩阵 AA,有一组特征向量 VV,将这组向量进行正交化、单位化,就能得到一组正交单位向量。特征值分解,就是将矩阵 AA 分解为如下式:
A=QΛQ−1A=QLambda Q^{-1}
- AA:nn 阶方阵
- QQ:由矩阵 AA 对应的特征向量组成的矩阵
- ΛLambda:对角阵,正对角线上各元素由矩阵 AA 对应的特征向量组成,其中按特征值大小进行依次排列
线性代数中相似矩阵的定义:设 A,BA,B 是 nn 阶方阵,若存在可逆矩阵 PP,使得 B=P−1APB=P^{-1}AP,则称矩阵 AA 与 BB 相似,P−1APP^{-1}AP 是对 AA 做相似变换。
参考相似矩阵的定义,左乘 PP 右乘 P−1P^{-1} 和左乘 P−1P^{-1} 右乘 PP 只是参考矩阵 AA 和 BB 的角度不同,例如在 B=P−1APB= P^{-1}AP 中两边分别左乘一个矩阵 PP 右乘一个矩阵 P−1P^{-1} 后可以得到 A=PBP−1A=PBP^{-1}。
在特征值分解中,不难看出其本质就是矩阵 AA 与对角阵 ΛLambda 相似,也就是 矩阵AA 可相似对角化,ΛLambda 为矩阵 AA 的相似标准型,A=QΛQ−1A=QLambda Q^{-1},其中对角阵中各元素为矩阵 AA 对应的特征值,特征值对应的特征向量组成矩阵 QQ,即特征值分解其实就是矩阵 AA 的相似对角化。
特征值分解基本流程:已知 nn 阶方阵 AA,有以下式子:
Ax=λx→(A−λE)X=0→∣A−λE∣=0Ax=lambda x rightarrow (A-lambda E)X=0 rightarrow |A-lambda E|=0
以此求出方阵 AA 的特征值,然后将各个特征值代入到原式 ∣A−λE∣|A-lambda E| 中求出各特征向量,以特征值和特征向量得到 A=QΛQ−1A=QLambda Q^{-1} 中矩阵 QQ 和对角阵 ΛLambda。
一个高维矩阵本质上就是高维空间中的一个线性变换,该矩阵的各特征值反映该矩阵各个方向变换的程度,我们取从大到小特征值中的前 mm 个特征值,也就是提取了该矩阵变化程度最大的 mm 个方向,其它方向进行舍弃,也就提取了该矩阵中最重要的 mm 个特征,以此来近似原矩阵的线性变换。
特征值分解的局限性:特征值分解公式中的矩阵 AA 为 nn 阶方阵,但是在绝大多数情况下数据集中的特征和样本数都是不同的,特征值分解无法应用于这种情况。
SVD:奇异值分解
相比特征值分解,奇异值分解可以适用于任意形状的矩阵,公式如下:
A=UΣVTA=USigma V^T
- AA:形状为 m∗nm*n 的矩阵
- UU:形状为 m∗mm*m 的半正交矩阵,其中的各向量被称为左奇异向量;该矩阵也被称为左奇异矩阵
- ΣSigma:形状为 m∗nm*n 的矩阵,该矩阵除了正对角线以外其它位置的元素均为0,正对角线上的元素被称为奇异值;该矩阵也被称为对角奇异矩阵
- VV:形状为 n∗nn*n 的半正交矩阵,其中的各向量被称为右奇异向量;该矩阵也被称为右奇异矩阵
奇异值的数量:通过对角奇异矩阵的形状为 m∗nm*n 可以发现奇异值的数量为 min(m,n)min(m,n)
基于特征值分解的局限性,我们引入奇异值分解,奇异值分解的思路与其大同小异,假设有形状为 m∗nm*n 的矩阵 AA,则有形状为 n∗nn*n 的方阵 ATAA^TA 和 形状为 m∗mm*m 的方阵 AATAA^T
为什么要使用矩阵 ATAA^TA 和矩阵 AATAA^T?
- AATAA^T 为矩阵 AA 中各行向量的内积;而 ATAA^TA 为矩阵 AA 中各列向量的内积,两矩阵都为方阵,对角线上各元素为各向量与其本身的内积,非对角线元素则为向量与其它向量的内积。
- 以矩阵 AATAA^T 为例,矩阵 AATAA^T 中对角线为各行向量与其本身的内积,为各向量的模;各非对角线元素,也就是矩阵 AATAA^T 中的各个行向量则为各行向量与其它行向量的内积,表示每个向量与其它向量的相似度,也就是该矩阵中的各个行向量反映矩阵 AA 中该行向量本身和该行向量与矩阵 AA 中的其它行向量之间的关系。
- 同理,矩阵 ATAA^TA 中的各列向量则表示矩阵 AA 中各列向量本身信息和各列向量之间的关系。
- 内积在一定程度上可以表示两向量之间的相似度;两向量之间的余弦相似度为两向量之间的内积除以向量 aa 和向量 bb 的模,因此当将两向量都归一化在同一范围中(统一量纲,一般在[0,1]之间)之后,向量的内积也就能反应出余弦相似度了。
- 矩阵 ATAA^TA 和矩阵 AATAA^T 分别保留了原矩阵 AA 中各列向量和各行向量的信息,且为方阵。
基于与特征值分解中相同的求特征值的公式,有:
(ATA)vi=λvi→(ATA−λE)vi=0→∣ATA−λE∣=0(AAT)ui=λui→(AAT−λE)ui=0→∣AAT−λE∣=0(A^TA)v_i = lambda v_i rightarrow (A^TA – lambda E)v_i = 0 rightarrow |A^TA – lambda E| = 0 \
(AA^T)u_i = lambda u_i rightarrow (AA^T – lambda E)u_i = 0 rightarrow |AA^T – lambda E| = 0 \
分别求出方阵 ATAA^TA 和 方阵 AATAA^T 的各特征值和特征向量后,由于:
ATA=(UΣVT)T(UΣVT)=VΣTUTUΣVT=VΣ2VT=VΣ2V−1AAT=(UΣVT)(UΣVT)T=UΣVTVΣTUT=UΣ2UT=UΣ2U−1A^TA = (USigma V^T)^T(USigma V^T) = VSigma^T U^T U Sigma V^T = VSigma^2 V^T = VSigma^2 V^{-1} \
AA^T = (USigma V^T)(USigma V^T)^T = U Sigma V^T VSigma^T U^T = USigma^2 U^T = USigma^2 U^{-1}
- 性质一:ΣSigma 为对称矩阵,对于对称矩阵有 ΣTΣ=Σ2Sigma^T Sigma = Sigma^2
- 性质二:V,UV,U 为半正交矩阵,对于正交矩阵有 VTV=E(单位阵)V^T V = E(单位阵),VT=V−1V^T=V^{-1}
正交矩阵的定义:若 nn 阶矩阵满足 AAT=ATA=EAA^T=A^TA=E,则称 AA 为 nn 阶正交矩阵。
明显正交矩阵一定为方阵。但在此定义下正交矩阵的行向量和列向量都相互正交,若仅满足 AAT=EAA^T=E 或 ATA=EA^TA=E,则这时候仅满足各行向量或各向量相互正交,这时候并不一定要求矩阵一定为方阵,因此这里说矩阵 V,UV,U (不一定为方阵) 为半正交矩阵。
从中我们可以看出方阵 ATAA^TA 的各个特征向量组成矩阵 UU,而方阵 AATAA^T 的各个特征向量组成矩阵 VV,矩阵 U,VU,V 中的各特征向量又被分别称为左奇异向量和右奇异向量。
矩阵 ATAA^TA 的各特征值表示矩阵 AA 中各列向量的重要程度;矩阵 AATAA^T 的各特征值表示矩阵 AA 中各行向量的重要程度。
奇异值求解的思路有两种;第一种方法为:
A=UΣVT→AV=UΣVTV→AV=UΣ→Avi=σiui→σi=AviuiA=USigma V^T rightarrow AV = USigma V^T V rightarrow AV = USigma rightarrow Av_i = sigma_i u_i rightarrow sigma_i = frac{Av_i}{u_i}
- 这里的除法为对应位置元素相除
AV=UΣ→Avi=σiuiAV=USigma rightarrow Av_i = sigma_i u_i,展开这个矩阵方程可得:
AV=UΣ→A[v1⋯vn]=[u1⋯um][σ1σ2⋯]→Avi=σiuiAV=USigma rightarrow Abegin{bmatrix} v_1 cdots v_n end{bmatrix} = begin{bmatrix}u_1 \ cdots \ u_mend{bmatrix} begin{bmatrix} sigma_1 & & \ & sigma_2 & \ & & cdots end{bmatrix} rightarrow Av_i=sigma_iu_i
- viv_i:第 ii 个右奇异向量
- uiu_i:第 ii 个左奇异向量
由于我们上面方阵 ATAA^TA 和 AATAA^T 的展开式中的对角奇异矩阵为 Σ2Sigma^2,因此易得第二种方法:
σi=λisigma_i = sqrt{lambda_i}
- 这里的 λilambda_i 为第 ii 个特征值
- 方阵 AATAA^T 和 ATAA^TA 的各特征值相同,但特征值对应的特征向量不同,即左奇异向量和右奇异向量不同
左奇异矩阵 UU 的特征值的平方根和右奇异矩阵 VV 的特征值的平方根在大小上是相等的,并且它们都等于矩阵 AA 的奇异值
在矩阵 ΣSigma 中各奇异值一般也是从大到小排列的。
通过SVD的过程我们可以看出,奇异值分解相当于将特征值分解中限定的方阵拓展到了任意形状的矩阵,但这仅仅是SVD相对EVD的其中一个优点,我们这里介绍的奇异值分解本质上还是基于特征值分解的,但奇异值分解现在也有了很多相对特征值分解更高效的算法,使得计算效率得到大大提高。对于这些算法在这里就不进行详细说明了(文章内容过多),感兴趣可以自行查找资源。
PCA基本原理
注意:以下在介绍过程中若不声明则默认特征矩阵的各行向量为特征,列向量为不同样本数据
特征降维的目的:减少特征的数量的同时,又保留大部分有效信息。
- 将带有重复信息的特征进行合并,删除无效特征等,创造出能够代表原特征矩阵大部分信息的特征更少的新特征矩阵
PCA的本质就是在特征空间中通过某种线性变换将某些特征进行合并,合并后的新特征保留原来特征的大部分信息,即拥有最大的方差,被合并的特征在经过这种变换后方差很小,因此不再能提供有效信息,我们将其舍弃(被丢弃的特征向量被认为信息量很少, 这些信息很可能就是噪音),从而实现特征降维;此外,保留合并的新特征被称为 “主成分”。
如图所示,图中通过旋转坐标轴将特征 x1x_1 与 x2x_2 合并为 x1∗x_1^*,通过观察可以看出原始特征 x1x_1 和 x2x_2 的方差都为 23frac{2}{3},变换处理后的特征 x1∗x_1^* 和 x2∗x_2^* 的方差分别为 43frac{4}{3} 和 00,此时 x2∗x_2^* 的方差过小,我们不认为该特征能有效地提供信息,因此舍去,保留 x1∗x_1^*。
此案例数据比较完美,也就是合并后 x2∗x_2^* 的方差为0,x1∗x_1^* 的方差等于原始特征 x1x_1 和 x2x_2 的方差之和,这里PCA认为保留了原来100%的信息,具体方法下方会详细介绍。
PCA主成分分析分别基于特征值分解和奇异值分解,因此一般有两种实现方法。
PCA的两种实现方法
基于特征值分解:假设数据集 X=[x1x2⋯xn]X=begin{bmatrix}x_1 & x_2 & cdots & x_n end{bmatrix},共 nn 个特征,需要降至 kk 维,一般流程如下:
- 去中心化处理,即每个特征下的各个样本值减去各自均值
- 计算协方差矩阵,即 1nXXTfrac{1}{n}XX^T
- 基于特征值分解方法求协方差矩阵 1nXXTfrac{1}{n}XX^T 的特征值和特征向量(这里除不除以 1nfrac{1}{n} 对结果没有影响)
- 对特征值从大到小排序,选择其中最大的 kk 个特征值
- 将选中的 kk 个特征向量分别作为行向量组成新的特征向量矩阵 PP
- 将原数据映射新空间中,即 Y=PXY=PX,其中 YY 为降维后的特征矩阵
X(m,n)→XX(m,m)T→P(k,m)→Y(m,k)=P(k,n)X(m,n)X_{(m,n)} rightarrow XX^T_{(m,m)} rightarrow P_{(k,m)} rightarrow Y_{(m,k)}=P_{(k,n)}X_{(m,n)}
基于奇异值分解:假设数据集 X=[x1x2⋯xn]X=begin{bmatrix}x_1 & x_2 & cdots & x_n end{bmatrix},共 nn 个特征,需要降至 kk 维,一般流程如下:
- 去中心化处理,即每个特征下的各个样本值减去各自均值
- 计算协方差矩阵,即 1nXXTfrac{1}{n}XX^T
- 基于奇异值分解方法求协方差矩阵 1nXXTfrac{1}{n}XX^T 的特征值和特征向量(这里除不除以 1nfrac{1}{n} 对结果没有影响)
- 对左奇异矩阵从大到小排序,选择其中最大的 kk 个特征值
- 将选中的 kk 个特征向量分别作为列向量组成新的特征向量矩阵 PP
- 将原数据映射新空间中,即 Y=PXY=PX,其中 YY 为降维后的特征矩阵
对矩阵 XX 经过去中心化处理然后乘以其本身的转置后再乘以 1nfrac{1}{n} 就得到了原特征矩阵 XX 对应的协方差矩阵;其中 XXTXX^T 中 XX 中的各行表示各特征,若各列表示各特征(我们平时遇见的大多情况),则为 XTXX^TX,因为我们要取每个特征的信息。
使用协方差矩阵 XXTXX^T 的原因:
可以发现,无论是基于EVD还是基于SVD,计算的都是协方差矩阵 XXTXX^T(忽略这里的 1nfrac{1}{n}),对于EVD,EVD只能处理方阵,而我们的数据集的大小绝大多数情况都不是方阵,为了基于EVD对特征矩阵进行降维,参考我们在介绍SVD时的 AATAA^T,同理这里的 XXTXX^T 中的各个行向量表示原特征矩阵 XX 中各行向量的信息(该向量本身的信息和该向量与其它行向量的相似度),即矩阵 XXTXX^T 相比原特征矩阵 XX 保留了各行向量(特征)的重要信息且将原矩阵变换为了方阵,我们可以通过这种方式使得基于EVD来处理非方阵的特征矩阵;此外,若为SVD则直接使用 XXTXX^T 进行后续操作即可,无需再次进行变换。
此外,协方差矩阵 XXTXX^T 为实对称矩阵,实对称矩阵的相似对角化后对应于不同特征值的各特征向量之间相互正交,也就是原矩阵中各特征向量表示在高维空间中各方向的变换。
这里的介绍为原因之一,下方会进行补充说明
基于EVD的PCA为什么可以使用协方差矩阵 XXTXX^T 来计算:参考SVD的原理,在SVD中 XXTXX^T 对应的各特征向量组成的矩阵也就是左奇异矩阵,我们是根据左奇异矩阵来对原特征矩阵进行变换的,因此使用 XXTXX^T 也能达到变换的效果,这也是在基于SVD的PCA中我们不用再重复计算 XXTXX^T 的原因。
基于SVD相比基于EVD在结果上的不同:我们在讲解SVD的时候,求了左奇异矩阵 UU 和右奇异矩阵 VV,正如我们上面所讲的,左奇异矩阵中的各特征向量和对应特征值表示的是原特征矩阵中各行向量(特征数据)的重要程度,同理,右奇异矩阵中的各特征向量和对应特征值表示的是原特征矩阵中各列向量(样本数据)的重要程度,因此相比基于EVD,基于SVD的方法中左奇异矩阵用于压缩行(特征),而左奇异矩阵用于压缩列,即使用SVD分解协方差的PCA可以实现两个方向的降维处理。
此外,由于现在的编程语言在大型矩阵的运算上都不是很擅长,需要占用大量的计算资源,相比基于EVD的PCA,而一些SVD算法实现可以不求协方差矩阵 XXTXX^T,直接求出左奇异矩阵,大大提高了计算效率,降低了时间成本,因此当前主流的PCA算法更多地基于SVD来进行实现。
最大方差理论
我们基于最大方差理论来进行介绍。
在信号处理中认为信号具有较大的方差,噪声有较小的方差,信噪比就是信号与噪声的方差比,越大越好。例如上图中样本在 u1u_1(红线)方向上的投影方差较大,在 u2u_2(绿线)方向上的投影方差较小,那么可认为 u2u_2 上的投影是由噪声引起的;因此认为特征的方差越大,对建模越有用,最好的 kk 维特征是将 nn 维样本数据降至 kk 维后,每个特征的样本方差都很大。
假如说我要要从图中的 u1u_1 和 u2u_2 选一条来做投影,那按照最大方差理论肯定是选 u1u_1 更好。
PCA使用的信息量衡量指标:样本方差,又称可解释性方差;一定程度上,某特征的方差越大,说明样本区分度越高,提供的信息越多;在PCA中,认为方差越大则特征所能表达的信息就越多。
对于EVD和SVD的补充:相似对角化后得到的对角阵 ΣSigma 正对角线上的元素由矩阵的各特征值组成,我们使用的原矩阵不是原特征矩阵,而是原特征矩阵对应的协方差矩阵,协方差矩阵的各向量为某特征与其它各个特征的协方差,所求的该矩阵的每个特征值则表示原矩阵在每个主成分方向,也就是每个特征向量上的方差,我们基于特征值也就是基于方差来选择最大的k个主成分方向。
此外,基于方差作为衡量标准也是我们在进行PCA处理时选择原特征矩阵对应的协方差矩阵进行分解的重要原因,我们基于协方差矩阵得到由各特征向量组成的高维投影矩阵,其中各个主成分方向的重要性(特征值)也正是由方差来决定的,基于PCA我们可以用这种方式发现这种潜在的性质,在选择前k个特征向量组成的投影矩阵进行变换时,被舍弃的方向则会投影到其它方向上,这也就是我们在最开始介绍的将含有相关性的特征进行融合的方式,虽看似简单粗暴,但事实证明这确实很有效。
重建原始数据:例如将原始的 mm 维数据降维至 kk 维数据,我们降维的公式为 Y=PXY=PX,等式两边同时左乘矩阵 PTP^T,得到逆转公式:
Xapprox(m,n)=P(m,k)TY(k,n)X_{approx(m,n)}=P^T_{(m,k)}Y_{(k,n)}
- XapproxX_{approx}:逆转后得到的矩阵
- PP:右奇异矩阵中由前 kk 个特征值对应的特征向量作为列向量组成的矩阵
- YY:降维后矩阵
重建原始数据得到的 XapproxX_{approx} 和原来的矩阵 XX 并不会完全相同,因为部分数据在降维过程中已经丢失了,对于丢失的那部分数据,在回复后该部分数据仍然符合整体分布,PP 矩阵矩阵仅包含了主要的规律特征信息,丢失了随机扰动(噪声),因此,很容易可以联想到该功能可以用于数据降噪,在接下来会展开说明。
平均平方映射误差(Average Squared Projection Error) :原始数据 xix_i 和映射值 xapprox(i)x_{approx(i)} 之间的差。即 xix_i 和其在低维表面上的映射点之间的距离的平方;描述两向量之间的相关性。
ASPE=1N∑i=1N∣∣xi−xapprox(i)∣∣2ASPE = frac{1}{N}sum^{N}_{i=1}||x_i – x_{approx(i)}||^2
数据的总方差 (Total Variation) :这些样本向量 xix_i 与本身的内积(方差)的均值;描述数据本身的整体变异程度(平均来看,我的训练样本距离零向量多远)。
Var=1N∑i=1N∣∣xi∣∣2Var = frac{1}{N}sum^N_{i=1}||x_i||^2
PCA所做的事是尽可能最小化平均平方映射误差,对于降至目标维度 kk,超参数 kk 的选择可以依据下方的经验法则:
平均平方映射误差数据的总方差=1N∑i=1N∣∣xi−xapprox(i)∣∣21N∑i=1N∣∣xi∣∣2≤0.01frac{平均平方映射误差}{数据的总方差}=frac{frac{1}{N}sum^{N}_{i=1}||x_i – x_{approx(i)}||^2}{frac{1}{N}sum^N_{i=1}||x_i||^2} le 0.01
数据的总方差越大越好,平均平方映射误差越小越好,其比值越小越好,我们一般用以下式子:
∑i=1kSii∑i=1NSii≥0.99frac{sum_{i=1}^k S_{ii}}{sum_{i=1}^N S_{ii}} ge 0.99
即协方差矩阵中选择的 kk 个特征值对应的特征向量的方差之和与原来全部特征向量的总方差的比值,若结果大于0.99,我们则认为降维后的矩阵保留了原来矩阵99%以上的信息。
解释方差(explained variance)与解释方差比(explained variance ratio):对于使用PCA降维后的数据,计算各个主成分的方差,每个主成分对应的方差被称为解释方差,将每个主成分对应的方差分别除以原来数据的总方差得到各个主成分对应的解释方差比,将降维后的各个主成分对应的解释方差比求和得到降维后数据的总解释方差比,该值反映保留了多少原数据信息。
- 例如降维前数据各特征总方差为100,降维后各主成分(共3个主成分)对应方差分别为80,10,5,则这里的方差就是解释方差,对应的比例 0.8, 0.1, 0.05 就是各主成分的解释方差比,总解释方差比为0.95,是原数据总方差的95%,则我们认为降维后数据保留了原数据95%的信息。
Whiten:数据的白化处理
数据的白化处理是数据预处理中常见的方法之一,一些PCA算法在使用PCA处理数据之前可以选择地先对数据进行白化处理,若不想深入了解,则可跳过本段。
对数据进行白化处理的原因:在机器学习的过程中,输入数据常为各种测量结果,相邻采样点之间的数据具有很强的相关性。比如一幅图像中,非边缘部分的相邻像素与中心像素之间差异一般很小,如果全部输入网络中进行训练,就会有很大的数据冗余,使得整个网络的训练效率不高。如果对输入数据进行白化处理,就能降低数据之间的相关度,不同数据所蕴含的信息之间的重复性就会降低,网络的训练效率就会提高。
“白化”的由来(通俗的定义):”白化”一词最早来源于信号处理领域,跟其中最常见的一种噪声——白噪声有很大的联系;在信号处理理论中,白噪声指的是一种在不同频率下都有相同功率的随机信号,即其功率谱密度为常数,功率与频率无关;”白噪声” 这个名字则来源于白光,由于白光包含了光谱中所有的颜色,且其功率谱密度也呈平坦状;类似地,”白噪声”这种噪声中包含了所有频率的噪声,因此被冠上了”白”的称号。其他不满足该功率谱密度特性的噪声也被称为有色噪声,类似于非白光的其他颜色光。
“白噪声” 的数学定义:若一个随机向量 ww 为一个白色随机向量,则其平均值函数为 00 且自相关函数为一个单位矩阵的倍数。即该信号的平均值为0,且各个分量之间互不相关。
数据的白化处理:了解以上定义后,不难理解所谓数据白化处理,不过就是将原数据进行零均值化(均值为0)和白化变换(各维度间不相关)
零均值化(中心化):某变量下所有样本值减去该变量的均值
白化变换:白化变换有多种算法,在这里我们仅介绍PCA白化这种变换原理
PCA白化的实现与基于EVD的PCA降维类似,即在求得特征矩阵对应的协方差矩阵的特征值与特征向量以后,用特征向量矩阵的转置左乘原始数据矩阵以实现对数据的旋转变换,再对变换后数据矩阵每一维除以对应方差(即特征值)即可。
使得特征矩阵对应的协方差矩阵变换为单位矩阵(正对角线元素为1,非对角线元素为0),
要求白化后的方差为1的原因:例如在信号处理中,一个随机信号 u 通过一个放大率为100的线性放大器得到另一个随机变量 w ,从方差的角度看其方差变为原来的10000倍,看起来是另一个与 u 完全不同的信号,但是实际上只是数值按比例放大,比如采用了不同的计量单位来量化同一个物理信号,两个看似不同的信号却描述了同一个物理现象。为了避免这种情况,最好在白化变换中将信号的整体方差也归一化为1,这样就可以避免类似的问题。
PCA算法的相关使用
对于PCA算法,除了使用其进行特征降维处理,也常用于高维数据的可视化和数据降噪处理,在这里,本文针对PCA举几个常见案例进行讲解。
高维数据可视化
我们常常无法理解四维及四维以上空间,但我们可以通过降维算法将高维数据降至三维以下然后对其进行可视化处理,通过使用PCA可以将高维数据 “压缩” 到低维,在保留原特征中大部分信息的情况下便于进行我们可以理解的可视化展示。
相关代码说明(可忽略并折叠):
from sklearn.decomposition import PCA
from sklearn.datasets import load_iris
import matplotlib as mpl
import matplotlib.pyplot as plt
plt.style.use('seaborn')
mpl.rcParams.update(rc)
iris = load_iris()
pca = PCA(n_components=2)
data_pca = pca.fit_transform(iris.data)
plt.figure(figsize=(10,6))
for i in set(iris.target):
plt.scatter(data_pca[iris.target==i, 0], data_pca[iris.target==i, 1], s=40, edgecolor='white', linewidth=1, alpha=0.8, label=iris.target_names[i])
plt.title('Iris Dataset Visonlization')
plt.legend(frameon=True)
plt.show()
如图所示,以特征数量为4,标签分类数量为3,样本数量为150的鸢尾花数据集为例,我们通过使用PCA将其降至2维,保留了原来约 97.77% 的信息,将二维数据可视化为散点图,其中同一颜色的样本表示所属同一类别。
PCA降维后的信息保存量
以LFW数据集为例,该数据集中包含各人脸图像数据,对于图像数据,我们也能通过可视化右奇异矩阵来理解右奇异矩阵到底在做什么。
实验设计:原数据中每张图像形状为 (62, 47),在统计机器学习中常将其展平处理,每个像素作为一个特征,则特征维度为 62 * 47 = 2914,使用PCA算法将其降至150维,由于我们对列进行降维,则有 Y(m,150)=X(m,2914)V(2914,150)TY_{(m,150)}=X_{(m,2914)}V_{(2914,150)}^T,两边同时右乘矩阵 V(150,2914)V_{(150, 2914)} 可得到逆转公式以重建数据,可视化展示部分右奇异矩阵 V(150,2914)V_{(150, 2914)},最后将经过PCA处理的低维数据经过逆转公式重新投影到原来的高维并进行可视化,观察结果。
相关代码说明(可忽略并折叠):
from sklearn.datasets import fetch_lfw_people
import matplotlib as mpl
import matplotlib.pyplot as plt
plt.style.use('seaborn')
rc = {
'font.family': 'Microsoft YaHei',
'mathtext.fontset': 'stix',
}
mpl.rcParams.update(rc)
# 原数据
lfw = fetch_lfw_people(min_faces_per_person=60) # (1348, 2914)
# 降维后得到右奇异矩阵
pca = PCA(n_components=150)
lfw_pca = pca.fit_transform(lfw.data) # (1348, 150)
V = pca.components_ # (150, 2914)
image_shape = lfw.images.shape # (1348, 62, 47)
# 将降维后的图像数据重新投影到高维
lfw_inverse = pca.inverse_transform(lfw_pca)
# 每行10张图,第一行为原图像,第二行为右奇异矩阵的可视化结果
fig, axes = plt.subplots(3,10,figsize=(8,3), subplot_kw={'xticks': [], 'yticks': []})
for i, ax in enumerate([*axes.flat]):
if i < 10:
ax.imshow(lfw.images[i, :, :], cmap='gray')
elif i < 20:
ax.imshow(V[i-10, :].reshape(62, 47))
else:
ax.imshow(lfw_inverse[i, :].reshape(62, 47), cmap='gray')
text_lst = ['原图像数据, shape: (62, 47)', '右奇异矩阵, shape: (62, 47)', '将低维数据重新投影到高维']
for i, text in enumerate(text_lst):
fig.text(0.95, 0.75 - i * 0.27,
text,
style = 'italic',
fontsize = 11,
color = "black")
plt.show()
上图中第一行为原数据集中的前十张图像,第二行为进行PCA处理后得到的对应的右奇异矩阵,第三行为使用逆转公式重新将降维后的图像投影到原来的高维得到的图像数据。
右奇异矩阵的可视化:右奇异矩阵的形状为 (150,2914)(150, 2914),也就是共150行,2914列,其中2914列表示的是特征数据,而这150行则分别表示150个主成分方向,降维后的数据为k维,也就是150个主成分,每个主成分都对应一个主成分方向,每个主成分方向反映其对应的主成分在 “关注” 图像中的哪部分,若为一些结构化数据,我们则不能直接对该矩阵进行可视化,因为我们知道它表示的是主成分方向,但却不知道这些矩阵中的数字有什么含义,我们这里以图像数据为例的目的就是这种情况下右奇异矩阵中的主成分方向可以被很容易地可视化出来,我们也就很方便就观察到这些主成分方向到底 “长什么样子”。
重建数据结果分析:对于我们重建的结果(第三行数据),观察可知虽然图像形状与原数据相同,但却变得更模糊了,这就是我们上方所介绍的,由于部分数据已经在降维过程中丢失了,我们使用这种方法恢复数据只能将未丢失的数据恢复,丢失的数据则符合与未丢失的数据相同的分布,但并不可能和原来的数据一模一样,但同时也达到了一定的降噪效果。
数据降噪
实验设计:以手写数字数据集为例,手动为数据集添加噪声,再使用PCA进行降噪处理,查看效果。
![image.png](https://p3-juejin.byteimg.com/tos-cn-i-k3u1fbpfcp/2d9c7df940174b92b61b8b159fa51e17~tplv-k3u1fbpfcp-jj-mark:0:0:0:0:q75.image#?w=1099&h=194&s=36053&e=png&b=e9e9e9)
from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
import numpy as np
digits = load_digits()
# 用正态分布将原来的像素值作为均值,设置方差,抽取的样本值作为新的像素值,以此来给图片添加噪声
X_noisy = np.random.normal(digits.data,2)
# 降维和降噪处理
pca = PCA(0.5).fit(X_noisy)
X_pca = pca.transform(noisy)
X_inverse = pca.inverse_transform(X_pca)
fig, axes = plt.subplots(2,10,figsize=(10,2)
,subplot_kw = {"xticks":[],"yticks":[]}
)
for i, ax in enumerate(axes.flat):
if i < 10:
ax.imshow(X_noisy[i].reshape(8,8),cmap="binary")
else:
ax.imshow(X_inverse[i-10].reshape(8, 8), cmap='binary')
ax.imshow
fig.text(0.92, 0.7, '降噪处理前', fontsize=10, color='black')
fig.text(0.92, 0.28, '降噪处理前', fontsize=10, color='black')
观察结果图像可知,降噪处理后的数据相对处理前的数据效果很显著,明显看到噪音数据变少了,以我们在介绍最大方差理论时的图像为例,方差最大的方向为信号,方差小的方向为噪音,则PCA其实就是在一定程度上将方差小的方向的数据去除掉,然后再加上与去掉的这部分数据数量相同的、与方差最大方向数据同一分布的数据,也就达到了降噪的效果。
补充说明
① 运算复杂:现在的编程语言在大型矩阵的运算上都不是很擅长,因此PCA往往需要占用很大的算力资源,计算较为缓慢,但由于SVD算法不仅仅局限于我们今天介绍的基于EVD的算法,它也衍生出很多高效算法,使得基于SVD的PCA在计算上更加高效,但总体来说仍然需要占用较大的算力资源。
② 主成分不具有可解释性:PCA是将原始特征进行压缩,降维后得到的主成分与原始特征不同,是通过某种变换组合起来的新特征。通常来说,在新的特征矩阵生成之前,我们不知道PCA建立了哪些新特征向量,新特征矩阵生成之后也不具有可读性。因此PCA一般不适用于探索特征和标签之间的关系的模型(如线性回归),因为无法解释的新特征和标签之间的关系不具有意义。
③ SVD中的左/右奇异矩阵:例如对于右奇异矩阵以及PCA的逆转公式,在处理图像时,我们在计算得到右奇异矩阵后,选择前 kk 个特征值对应的特征向量组成新矩阵 V(k,n)V_{(k,n)},当处理纯数值时,我们无法观察该矩阵的含义,但是若我们处理的是图像,由前 kk 个特征值对应的特征向量组成的该矩阵表示的就是从原图像中提取的压缩后的主要信息,也就是各主成分方向,我们则可以通过可视化该空间矩阵,并将可视化结果与原图像进行对比以查看PCA主要从原图像中提取哪些信息。
⑤ 噪音过滤:PCA的逆转公式并不能真正完全将降维后的数据映射成和原来一模一样的特征矩阵,因为原来的信息在合并时就已经丢失了,降维并不是完全可逆的,可逆转是因为采用某种策略或算法基于降维后的数据对原来数据的预测;此外,降维的目的之一就是希望抛弃掉对模型带来负面影响的特征,带有效信息的特征的方差应该是远大于噪音的,所以相比噪音,有效的特征所带的信息应该不会在PCA过程中被大量抛弃,PCA的逆转能够在不恢复原始数据的情况下,将降维后的数据返回到原本的高维空间,因此,利用这个性质我们可以实现数据噪音的过滤。
Reference
- 主成分分析(PCA)原理详解 – 知乎 (zhihu.com)
- mp.weixin.qq.com/s/Dv51K8JET…
- 菜菜的sklearn机器学习
- 为什么特征值是“沿对应的特征向量的数据的方差”?_百度知道 (baidu.com)
- 数据预处理之白化(Whitening transformation)_数据白化_执笔论英雄的博客-CSDN博客
- 白化变换:PCA白化、ZCA白化 – 知乎 (zhihu.com)