机器学习——数据降维——主成分分析(PCA)和奇异值分解(SVD)

一、主成分分析(PCA)

主成分分析,Principal Component Analysis (PCA),是现代数据分析的标准工具,它可以把庞大复杂的高维数据集(变换之前的数据集各维度之间存在相关性),通过数学变换,转化成较低维度的数据集,并去除掉维度之间的相关性。

PCA的优势和劣势都在于,它是一个非参数的分析。没有需要调整的参数,它的答案是唯一的,独立于用户特性的。

PCA和LDA的区别可见《多元统计分析——数据降维——Fisher线性判别分析(LDA)》中第五节。

主成分分析的应用领域:

  • 构建综合指标:主成分分析主要用于构建综合指标来区分目标群体,例如构建顾客各种消费行为的综合指标来进行客户分级。

  • 数据降维:通常我们会用少于原始变量数的主成分来描述尽可能多的数据差异,特别是当原始变量维度很高时,可达到降维目的。

  • 数据可视化:当原始变量维度很高时,可采用第一、第二主成分散点图来直观表述数据特征,例如数据聚类信息等。

  • 变量压缩、重构:由“重要的”主成分重构原始变量(即原始数据X\rightarrow Z之后,选择主要方差的特征之后,通过线性变化回原维度数据X',以实现降噪的目的),可以去除原始数据中冗余的噪音,突出数据的特征,例如人脸识别。

1、PCA原理

主成分分析的原理非常简单,概括来说就是选择包含信息量大的维度,去除信息量少的“干扰”维度。

注意:这边所谓的“维度”不是原始数据的某个特征,而是原始数据各个特征的组合。

具体如下:

  • 数据从原来的坐标系X转换到新的坐标系\mu ^{T}X,新坐标系的选择是由数据本身决定的。第一个新坐标轴\mu _{1}选择的是原始数据中方差最大(\mu^{T}_{1}\Sigma \mu_{1})的方向,第二个新坐标轴\mu _{2}选择与第一个新坐标轴\mu _{1}正交(\mu^{T}_{1}\Sigma \mu_{2}=0且具有最大方差(\mu^{T}_{2}\Sigma \mu_{2}的方向,第三个新坐标轴\mu _{3}选择与第一个新坐标轴\mu _{1}、第二个新坐标轴\mu _{2}正交(\mu_{1}^{T}\Sigma \mu_{3}=0,\mu^{T}_{2}\Sigma \mu_{3}=0)且具有最大方差(\mu^{T}_{3}\Sigma \mu_{3}的方向........以此类推,第p个新坐标轴\mu _{p}选择与第一个新坐标轴\mu _{1}、第二个新坐标轴\mu _{2}、...、第p-1个新坐标轴\mu _{p-1}正交(\mu^{T}_{j}\Sigma \mu_{p}=0,j=1,2,...,p-1)且具有最小方差(\mu^{T}_{p}\Sigma \mu_{p}的方向,共建立与原始数据特征数目相等的新坐标轴。
  • 我们会发现,大部分方差都包含在最前面的几个新坐标轴中,因此我们可以忽略余下的坐标轴,从而实现降维。(方差大代表不同数据之间的差异大,即,包含的可区分信息量大)

注意,PCA降维后,原始数据被映射到了新坐标系,不是原始值了。

pca降维实现的步骤:

1)把数据集组织成一个 m\times n 的矩阵X

2)将X的每个变量(列)分别进行去中心化X_{j}-\mu _{j},j=1,....,p,得到新的数据集W(注意:有时需要做标准化,见下方第四节——基于标准化变量的主成分分析)。

3)求出协方差矩阵\Sigma =\frac{1}{n}WW^{T}

4)求出协方差矩阵的特征值\lambda _{j}和对应的特征向量e _{j}

5)将特征向量按对应特征值大小从上到下按行排列成矩阵,取前 k行组成矩阵P(注意:k如何取值请见下方第5节——主成分个数的选择)。

6) Y=PW即为降维到 k 维后的数据。

【知识点1】:

An阶方阵,如果标量\lambdan维非零列向量x使关系式Ax=\lambda x成立,则称\lambda是方阵A的特征值,非零列向量x称为A的对应于特征值\lambda的特征向量。

Ax=\lambda x可改写为(A-\lambda E)x=0

这是n个未知数n个方程的齐次线性方程组,它有非零解的充要条件是\left | A-\lambda E \right |=0,即

\begin{vmatrix} a_{11}-\lambda & a_{12} & ... & a_{1n} \\ a_{21} & a_{22}-\lambda & ... & a_{2n} \\ .& .& .&. \\ .& .& .&. \\ .& .& .&. \\ a_{n1}& a_{n2} & ... & a_{nn}-\lambda \end{vmatrix}=0

以上步骤的关键点在于:为何协方差的特征值对应的特征向量构成的矩阵能够实现降维的效果?以下我们通过结合案例数据和数据分布的图形进一步进行讲解。

案例分析:

案例1:有以下二维数据:

 房价(百万)面积(平方米)
a109
b23
c12
d76.5
e32.5

1)对数据进行去中心化——均值为0

去中心化之后的数据分布如下所示:O为坐标原点。

 【小贴士1】

观察点X_{i}有时用经过中心化的X_{i}-\overline{X}代替。这并不会对样本协方差矩阵或者相关系数矩阵有影响,所以不会改变最后得到的系数向量e_{j}

主成分的几何解释?

经过中心化的一系列主成分Z_{ij}=e_{j}^{T}(X_{i}-\overline{X}),j=1,...,p,i=1,...,n可看作将原来坐标系的原点平移到\overline{X},然后旋转坐标轴,使得坐标轴与数据最大方差方向一致,如下图所示。

 2)数据从原来的坐标系转换到新的坐标系。第一个新坐标轴选择的是原始数据中方差最大的方向(即下图中,数据在蓝色线的投影后的方差最大),第二个新坐标轴选择与第一个新坐标轴正交且具有最大方差的方向(即下图中,数据在红色线的投影后的方差最大)。

我们知道,给定一个单位向量\mu和空间内一点X_{i}X_{i}根据空间维度可以表示为[X_{i1},X_{i2},X_{i3},.....X_{ip}]p代表维度),X_{i}\mu上的投影距离原点的距离为\mu ^{T}X_{i},同样,如果X_{i}为数据集合XX可表示为矩阵\begin{bmatrix} X_{11} & X_{12} & ... & X_{1p} \\ X_{21} & X_{22} & ... & X_{2p} \\ .& .& .&. \\ .& .& .&. \\ .& .& .&. \\ X_{i1}& X_{i2} & ... & X_{ip} \end{bmatrix},也可以表示为列向量\begin{bmatrix} X_{1} \\ X_{2} \\ . \\ . \\ . \\ X_{i} \end{bmatrix})中的一点,因此为了最大化投影点的方差,我们需要选择一个单位向量\mu将如下式子最大化:

\frac{1}{n}\sum ^{n}_{i=1}(\mu ^{T}X_{i})^{2}=\frac{1}{n}\sum ^{n}_{i=1}(\mu ^{T}X_{i})(\mu ^{T}X_{i})^{T}=\frac{1}{n}\sum ^{n}_{i=1}(\mu ^{T}X_{i}X_{i}^{T}\mu)=\mu^{T}\frac{1}{n}\sum ^{n}_{i=1}( X_{i}X_{i}^{T})\mu,其中n代表样本点的个数。

注意:实际上样本方差和协方差的公式应该是n-1,此处为方便用n,对结果没有影响。

其中X=\begin{bmatrix} X_{1} \\ X_{2} \\ . \\ . \\ . \\ X_{i} \end{bmatrix}X^{T}=\begin{bmatrix} X_{1}^{T} & X_{2}^{T} & . & . & . & X_{i}^{T} \end{bmatrix}XX^{T}=\sum ^{n}_{i=1}X_{i}X_{i}^{T},代入到\mu^{T}\frac{1}{n}\sum ^{n}_{i=1}( X_{i}X_{i}^{T})\mu,得:

\mu^{T}\frac{1}{n}\sum ^{n}_{i=1}( X_{i}X_{i}^{T})\mu=\mu^{T}\frac{1}{n}( XX^{T})\mu

\Sigma =\frac{1}{n}( XX^{T}),则上式最终简化为\mu^{T}\Sigma \mu\Sigma即为原始变量X=\begin{bmatrix} X_{11} & X_{12} & ... & X_{1p} \\ X_{21} & X_{22} & ... & X_{2p} \\ .& .& .&. \\ .& .& .&. \\ .& .& .&. \\ X_{i1}& X_{i2} & ... & X_{ip} \end{bmatrix}的协方差矩阵。

【知识点2】:

转置矩阵性质:(AB)^{T}=B^{T}A^{T}

我们是要求方差最大化,如果\mu越大的话,总体方差也就越大,鉴于我们只想求\mu的方向,而不想依赖于\mu的长度,我们让\mu为一个标准化的向量(单位向量)。接下来, 我们可以通过拉格朗日乘子法构造下面方程:

max\,\, \, \, \mu^{T}\Sigma \mu

\\ s.t\, \, \, \, \mu ^{T} \mu=1

关于拉格朗日乘子法详细解析请见【如何理解拉格朗日乘子法? 

对其变形:F=\mu^{T}\Sigma \mu+\lambda (\mu^{T} \mu -1)

对其求\mu的偏导:\frac{\partial F}{\partial \mu }=2\Sigma \mu -2\lambda \mu=0\Rightarrow \Sigma \mu=\lambda \mu即当\Sigma \mu=\lambda \mu时,\mu^{T}\Sigma \mu取最大值。

根据上面的【知识点1】,我们可以得出:\lambda是方阵\Sigma的特征值,\mu称为\Sigma的对应于特征值\lambda的特征向量。

此时最大方差为\mu^{T}\Sigma \mu=\mu^{T}\lambda \mu=\lambda\lambda是一个常数,可以提取出来\lambda \mu^{T}\mu,又因为\mu ^{T} \mu=1,得\mu^{T}\lambda \mu=\lambda)。

即得出结论,PCA降维,即为求协方差矩阵的特征值和特征向量,\Sigma =\frac{1}{n}( XX^{T})即为原始数据的协方差矩阵。

【定理1】对于任意向量\mup\times p对称矩阵\Sigma,记(\lambda _{j},e_{j}),j=1,...,p\Sigma的特征对,其中\lambda _{1}\geq ...\geq \lambda _{p}\geq 0,则

极值极值点
\max_{\mu \neq 0}\frac{\mu^{T}\Sigma \mu}{\mu^{T}\mu}=\lambda _{1}\mu_{1} =e_{1}
\max_{\mu \neq 0,\mu_{j} \perp e_{1},..., e_{j-1}}\frac{\mu^{T}\Sigma \mu}{\mu^{T}\mu}=\lambda _{j},j=2,...,p-1\mu_{j} =e_{j}
\min_{\mu \neq 0}\frac{\mu^{T}\Sigma \mu}{\mu^{T}\mu}=\lambda _{p}\mu_{p} =e_{p}

注意:以上\mu_{j} \perp e_{1},..., e_{j-1}的限制条件其实就是\mu^{T}_{j}\Sigma \mu_{k}=0,k=1,2,...,j-1。我们以\mu^{T}_{2}\Sigma \mu_{1}=0为例,对求\mu _{2}而言,我们已经求得\mu _{1}=e_{1}(\lambda _{1},e_{1})\Sigma的一个特征对,\Sigma \mu_{1}=\lambda_{1}e_{1} ,则\mu^{T}_{2}\Sigma \mu_{1}=0\Rightarrow \mu^{T}_{2}\lambda _{1}e_{1}=0。因为\lambda _{1}是一个常数,可以提到最前,得:\lambda _{1}\mu^{T}_{2}e_{1}=0\Rightarrow \mu^{T}_{2}e_{1}=0,最后得到:\mu_{2}e_{1}是垂直的(\mu_{2} \perp e_{1},因为\mu _{2}^{T}\mu _{2}可以看做空间中两个相对于原点对称的点,\mu_{2}^{T} \perp e_{1}\mu_{2} \perp e_{1})。

 【定理2】记(\lambda _{1},e_{1}),...,(\lambda _{p},e_{p})为协方差矩阵\Sigma的特征值-特征向量,\lambda _{1}\geq ...\geq \lambda _{p}\geq 0并且特征向量e_{1},...,e_{p}是标准化的特征向量。

则变量X_{1},...,X_{p}的第j个主成分由下式给出:

Z_{j}=e_{j}^{T}X=e_{j1}X_{1}+e_{j2}X_{2}+...+e_{jp}X_{p},j=1,...,p

性质1:var(Z_{j})=e_{j}^{T}\Sigma e_{j}=\lambda _{j}

性质2:cov(Z_{j},Z_{k})=e_{j}^{T}\Sigma e_{k}=0

性质3:\sum _{j=1}^{p}var(Z_{j})=\sum _{j=1}^{p}var(Y_{j}),即变换后的数据还原了原始数据的信息(方差总和相同)

证明:我们知道变化之后的方差之和为:\sum _{j=1}^{p}var(Z_{j})=\lambda _{1}+\lambda _{2}+...+\lambda _{p},原始数据的方差之和\sum _{j=1}^{p}var(Y_{j})为协方差矩阵\Sigma的对角线元素之和,即为\Sigma的迹tr(\Sigma );我们知道\Sigma的迹tr(\Sigma )即为协方差矩阵\Sigma所有特征值之和)

关于矩阵迹和特征值得关系可见《线性代数——韦达定理、矩阵行列式、矩阵的迹、矩阵特征值及关系

在这边,特征向量我们可以理解为转换的方向,特征值,我们可以理解为平移的速度。

 【知识点3】:

教科书上求n阶方阵的特征值和特征向量的步骤如下:

  • 求出n阶方阵A的特征多项式\left | A-\lambda E \right |

  • 求出特征方程\left | A-\lambda E \right |=0的全部根,\lambda _{1}\lambda _{2},......,\lambda _{n},即为A的特征值。

  • 把每个特征值\lambda _{i}代入线性方程组( A-\lambda E)x=0,求出基础解系,就是A对应于\lambda _{i}的特征向量,基础解析的线性组合(零向量外)就是A对应于\lambda _{i}的全部特征向量。

我们可以对上面的例子求下解。

  • 对原始数据集去中心化得到:\begin{pmatrix} 5.4 & 4.4\\ -2.6 & -1.6\\ -3.6 &-2.6 \\ 2.4&1.9 \\ -1.6& -2.1 \end{pmatrix},我们以X表示。

  • 求出协方差矩阵\Sigma =\frac{1}{n}X^{T}X=\begin{pmatrix} 5.4 &-2.6 & -3.6 & 2.4& -1.6\\4.4& -1.6& -2.6&1.9& -2.1 \end{pmatrix}\begin{pmatrix} 5.4 & 4.4\\ -2.6 & -1.6\\ -3.6 &-2.6 \\ 2.4&1.9 \\ -1.6& -2.1 \end{pmatrix}=\begin{pmatrix} 11.44 & 9.04\\9.04 & 7.34 \end{pmatrix}

  • 求解特征多项式\left | A-\lambda E \right |=\begin{vmatrix} 11.44-\lambda & 9.04\\9.04 & 7.34-\lambda \end{vmatrix}=0,得\lambda _{1}=18.65953\lambda _{2}=0.12047

  • \lambda _{1}=18.65953( A-\lambda E)x=\begin{pmatrix} -7.21953 & 9.04\\9.04 & -11.31953 \end{pmatrix}\begin{pmatrix} x_{1}\\x_{2} \end{pmatrix}=0,当\lambda _{2}=0.12047( A-\lambda E)x=\begin{pmatrix} 11.31953& 9.04\\9.04 &7.21953 \end{pmatrix}\begin{pmatrix} x_{1}\\x_{2} \end{pmatrix}=0,再结合求得的向量为单位向量\left \| \mu \right \|=1\Rightarrow \mu \mu ^{T}-1=0,求得\lambda _{1}的特征向量为\begin{pmatrix} 0.78139452\\0.62403734 \end{pmatrix}\lambda _{2}的特征向量为\begin{pmatrix} 0.62403734\\ -0.78139452\end{pmatrix}

  • 按照特征值大小排列,得主成分矩阵为\begin{pmatrix} 0.78139452&0.62403734\\0.62403734& -0.78139452\end{pmatrix}

2、PCA算法的python实现

import numpy as np
dataMat=pd.DataFrame([[10,9],
[2,3],
[1,2],
[7,6.5],
[3,2.5]])

meanVals=np.mean(dataMat,axis=0) #求dataMat各列均值
meanRemoved=dataMat-meanVals #减去原始数据中的均值,避免协方差计算中出现乘以0的情况
covMat=np.cov(meanRemoved,rowvar=0) #rowvar=0-->以列代表一个变量,计算各列之间的协方差
eigVals,eigVects=np.linalg.eig(np.mat(covMat)) #协方差矩阵的特征值和特征向量

eigVals
eigVects

输出:

array([23.32440667,  0.15059333])

matrix([[ 0.78139452, -0.62403734],
        [ 0.62403734,  0.78139452]])

去中心化的数据集乘以主成分矩阵,即为降维后的数据,注意:此时eigVects中一列是代表着一个特征向量。 

meanRemoved*eigVects

输出:

matrix([[ 6.96529471,  0.06833426],
        [-3.0300855 ,  0.37226585],
        [-4.43551736,  0.21490867],
        [ 3.0610178 , -0.01304002],
        [-2.56070965, -0.64246875]])

 3、利用sklearn中的PCA进行降维

import numpy as np
from sklearn.decomposition import PCA

#导入数据
x=np.array([[10,9],
[2,3],
[1,2],
[7,6.5],
[3,2.5]])
x

x_t=x-x.mean()  #去中心化
x_t

#降维
model_pca=PCA()  #建立PCA模型对象
model_pca.fit(x1)   #将数据输入模型
x_=model_pca.transform(x1)   #将数据集进行转换映射,x_为进行降维过后的数据集,形状和原数据集一样
components=model_pca.components_  #获得转换后的所有主成分
components_var=model_pca.explained_variance_    #获得各主成分的方差
components_var_ratio=model_pca.explained_variance_ratio_   #获得各主成分的方差占比
x_
components
components_var
components_var_ratio

输出:

 【参数解析】:

  • x:原始数据集
  • x_t:原始数据集去中心化后的数据集
  • x_:降维过后的数据集
  • components:转换过后的所有主成分,算法的重点在于求出此主成分(矩阵);原数据集去中心化之后,乘以该主成分(矩阵),即为降维过后的数据集。有如下关系

       \begin{pmatrix} 5.4 & 4.4\\ -2.6&-1.6 \\ -3.6 & -2.6\\ 2.4& 1.9\\ -1.6 & -2.1 \end{pmatrix}\begin{pmatrix} 0.78139452 &0.62403734 \\ 0.62403734& -0.78139452 \end{pmatrix}=\begin{pmatrix} 6.96529471 & -0.06833425\\ -3.030085496&-0.372265852 \\ -4.435517356 & -0.214908672\\ 3.061017794& 0.013040027\\ -2.560709646 & 0.642468748 \end{pmatrix}

  • components_var:转换过后的数据集x_各维度(注:形状和原始数据集一样,但是维度和原始数据集的完全不一样)的数据方差,即x_.var()
  • components_var_ratio:components_var中各维度的方差,占全部维度方差之和的比例。

4、基于标准化变量的主成分分析

如果变量X=(X_{1},...,X_{p})的波动性(由于度量单位不同等原因)差距过大,直接用协方差矩阵\Sigma生成的主成分会由方差大的变量主导。

这个时候,我们可对每个变量分别进行标准化:W_{j}=(X_{j}-\mu _{j})/\sqrt{\sigma _{jj}}j=1,...,p\sqrt{\sigma _{jj}}为单个变量的标准差。

则标准化之后的数据W=(W_{1},...,W_{p})=D_{s}^{-1}(X-\mu ),其中D_{s}=diag(\sqrt{\sigma _{11}},...,\sqrt{\sigma _{pp}})X为原始数据,\mu为原始数据的均值向量。

所以对标准化后的数据进行主成分分析,即求W的协方差矩阵cov(W),进而求cov(W)的特征值和特征向量。

cov(W)=cov(D_{s}^{-1}(X-\mu ))=D_{s}^{-1}(X-\mu )(X-\mu )^{T}(D_{s}^{-1})^{T},其中(X-\mu )(X-\mu )^{T}即为原数据的协方差矩阵\SigmaD_{s}为对角阵,则(D_{s}^{-1})^{T}=D_{s}^{-1},代入得:cov(W)=D_{s}^{-1}\Sigma D_{s}^{-1}=P(即为X的相关系数矩阵) 。

总结:由于cov(W)=D_{s}^{-1}\Sigma D_{s}^{-1}=P(即为X的相关系数矩阵),则对W_{1},...,W_{p}进行主成分分析,就相当于基于原变量X_{1},...,X_{p}的相关系数矩阵进行主成分分析。

【定理3】记(\widetilde{\lambda _{1}},\widetilde{e_{1}}),...,(\widetilde{\lambda _{p}},\widetilde{e_{p}})X的相关系数矩阵P的特征对,\widetilde{\lambda _{1}}\geq ...\geq \widetilde{\lambda _{p}}\geq 0并且特征向量\widetilde{e_{1}},...,\widetilde{e_{p}}是标准化的特征向量。

则对W_{1},...,W_{p}进行主成分分析,相当于基于原变量X_{1},...,X_{p}的相关系数矩阵P进行主成分分析,第j个主成分由下式给出:

V_{j}=\widetilde{e}_{j}^{T}W=\widetilde{e}_{j1}W_{1}+\widetilde{e}_{j2}W_{2}+...+\widetilde{e}_{jp}W_{p},j=1,...,p

性质:\sum _{j=1}^{p}var(V_{j})=\sum _{j=1}^{p}var(W_{j})=p(注意:相关系数矩阵P的对角线都为1)

 【小贴士2】

  • 基于相关系数矩阵P的主成分分析不受度量单位影响。

  • P\Sigma的特征值、特征向量是不同的,因此两组主成分的数值和方向都是不同的。一组主成分并不是另一组主成分的简单变换。

  • 如果变量数值差异很大,推荐进行标准化来获得一个较平衡的数据集。

 5、主成分个数的选择

在主成分分析中,为了有效地归纳数据、达到降维目的,应该使用多少个主成分?

5.1、百分比截点(Percentage cutoff)

使用足够多的主成分来反映一定百分比(比如80%)的总方差\frac{\sum ^{k}_{j=1}\lambda _{j}}{\sum ^{p}_{j=1}\lambda _{j}}=80\%

5.2、平均比截点(Average cutoff)

使用特征值大于平均特征值\frac{\sum ^{p}_{j=1}\lambda _{j}}{p}的主成分;对使用相关矩阵做的主成分分析来说,这个平均值为1。

平均截点准则广为使用,是很多软件包的默认准则。

当数据能够在相对较小的维度下被很好地归纳时,特征值可以明显地地分为“大特征值”和“小特征值”。

5.3、碎石图 (Scree graph):

画出\lambda _{j}关于j的图像,从图上寻找能区分“大特征值”和“小特征值”的分界点。

6、案例分析:主成分的解读

例:下表记录了52位学生6门功课的考试分数的部分数据,Y_{1}Y_{6}依次表示数学、物理、化学、语文、历史、英语成绩,部分数据如下:

6.1、计算相关系数矩阵

corr=pd.DataFrame(np.corrcoef(data,rowvar=0))
corr

输出:

6.2、数据标准化

#数据标准化
data_=(data-data.mean(0))/data.std()
data_

输出:

6.3、sklearn实现PCA降维

from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

from sklearn.decomposition import PCA
model_pca=PCA()  #建立PCA模型对象
model_pca.fit(data_)   #将数据输入模型
x_=model_pca.transform(data_)   #将数据集进行转换映射,x_为进行降维过后的数据集,形状和原数据集一样
components=model_pca.components_  #获得转换后的所有主成分
components_var=model_pca.explained_variance_    #获得各主成分的方差
components_var_ratio=model_pca.explained_variance_ratio_   #获得各主成分的方差占比

components
components_var
components_var_ratio

输出:

array([[-0.41205198, -0.38117795, -0.33213465,  0.46118457,  0.4205876 ,
         0.43013716],
       [-0.37597731, -0.356706  , -0.5626165 , -0.27852313, -0.41478358,
        -0.40650217],
       [ 0.21582978, -0.80555264,  0.46743533, -0.04426879, -0.25039004,
         0.14612244],
       [-0.78801362,  0.11755209,  0.58763655, -0.02783261, -0.03376008,
        -0.13410793],
       [ 0.0205822 , -0.21203603,  0.0333622 , -0.59904487,  0.73843439,
        -0.22218   ],
       [ 0.14450829, -0.14061074,  0.09068468,  0.59003773,  0.20479353,
        -0.74902427]])

array([3.70990438, 1.26248122, 0.44083653, 0.27050177, 0.16951589,
       0.14676021])

array([0.6183174 , 0.21041354, 0.07347275, 0.04508363, 0.02825265,
       0.02446003])

其中,components代表的是所有的主成分,第一行代表的是第一个主成分e_{1},第二行代表的是第一个主成分e_{2},以此类推。由components_var_ratio我们知道前两个主成分对应的方差之和为0.8287309(超过80%)

6.3、画碎石图,选择主成分个数

import matplotlib.pyplot as plt
plt.plot(components_var,marker='o')

输出:

我们也可以通过观测碎石图,得出结论:取前两个主成分进行进一步的研究。

6.4、解读

根据以上两个主成分,我们是可以写出主成分的构造方式,如下:

Z_{1}=-0.412Y_{1}-0.381Y_{2}-0.332Y_{3}+0.461Y_{4}+0.421Y_{5}+0.430Y_{6}

Z_{2}=-0.376Y_{1}-0.357Y_{2}-0.563Y_{3}-0.279Y_{4}-0.415Y_{5}-0.407Y_{6}

我们可以将其理解为两个综合指标,那么我们如何来解读呢?

由下图的分解,我们可以将Z_{1}指标理解为文理科的偏科情况,Z_{2}可以理解为全部科目的平均分。

 有了这两个指标的意义,进一步探索一些典型学生的样本主成分取值/得分?

x_.round(2)[[5,6,44],:]
data.iloc[[5,6,44],:]

输出: 

例如:5,6,44号同学,利用上面的结论:第一列代表的是偏科情况(数值越大文科越好,越小理科越好),说明这三位同学是理科成绩要好于文科(注意:“各科主成分得分”是经PCA降维之后的数据x_),对应我们查看其原始成绩数据:发现确实其理科成绩要优于文科成绩。

x_.round(2)[[29,48],:]
data.iloc[[29,48],:]

 输出:

29,48号同学,情况正好相反:文科成绩要好于理科。

print('-----各科主成分得分-----')
x_.round(2)[[25,32],:]
print('-----对应原始成绩-----')
data.iloc[[25,32],:]

输出:

25,32号同学,利用上面的结论:第二列代表的是均衡情况(数值越大总成绩越小,越小总成绩越高),说明这二位同学总成绩是比较高的(各科成绩相对均衡)。

print('-----各科主成分得分-----')
x_.round(2)[[7],:]
print('-----对应原始成绩-----')
data.iloc[[7],:]

输出:

7号同学,情况正好相反:总成绩比较低(不理想)。

二、奇异值分解(SVD)

1、SVD定义

SVD可以看成是对PCA主特征向量的一种解法,在上述PCA介绍过程中,为了求数据X的主特征方向,我们通过求协方差矩阵\Sigma =\frac{1}{n}X^{T}X的特征向量来表示样本数据X的主特征向量,但其实我们可以通过对X进行奇异值分解得到主特征方向,下面我们首先比较一下特征值分解和奇异值分解,然后分析一下特征值和奇异值的关系以及SVD在PCA中的应用。

SVD也是对矩阵进行分解,但是和特征分解不同,SVD并不要求要分解的矩阵为方阵。假设我们的矩阵A是一个m行×n列的矩阵(注意:只有行数m>列数n的矩阵才能进行svd分解),那么我们定义矩阵A的SVD为:

A=U\Sigma V^{T}

其中 U 是一个 m\times m 的矩阵, \Sigma 是一个 m\times n 的矩阵,除了主对角线上的元素以外全为0,主对角线上的每个元素都称为奇异值, V 是一个 n\times n 的矩阵。 U和 V都是酉矩阵,即满足U^{T}U=EV^{T}V=E。下图可以很形象的看出上面SVD的定义:

 2、如何求出SVD分解后的U,Σ,V这三个矩阵

  • V矩阵:将A的转置和A做矩阵乘法,得到n×n的一个方阵A^{T}A,得到的特征值和特征向量满足下式:(A^{T}A)v_{i}=\lambda _{i}v_{i}。这样可以得到矩阵A^{T}A 的n个特征值\lambda _{i}和对应的n个特征向量v_{i}了。将A^{T}A的所有特征向量张成一个n×n的矩阵V,就是我们SVD公式里面的V矩阵了。一般我们将V中的每个特征向量叫做A的右奇异向量。

  • U矩阵矩阵:将A和A的转置做矩阵乘法,得到m×m的一个方阵AA^{T},得到的特征值和特征向量满足下式:(AA^{T})u_{i}=\lambda _{i}u_{i}。这样可以得到矩阵AA^{T} 的m个特征值\lambda _{i}和对应的m个特征向量u_{i}了。将AA^{T}的所有特征向量张成一个m×m的矩阵U,就是我们SVD公式里面的U矩阵了。一般我们将U中的每个特征向量叫做A的左奇异向量。

  • \Sigma矩阵:由于Σ除了对角线上是奇异值其他位置都是0,那我们只需要求出每个奇异值σ就可以了。A=U\Sigma V^{T}\Rightarrow AV=U\Sigma V^{T}V=U\Sigma \Rightarrow Av_{i}=\sigma _{i}u_{i}\Rightarrow \sigma _{i}=Av_{i}/u_{i},根据\sigma _{i}=Av_{i}/u_{i}算出每一个奇异值\sigma _{i},组成的\Sigma矩阵等于\begin{pmatrix}\sigma_{1} & 0 & ... & 0 \\ 0 & \sigma_{2} & ... & 0 \\ .& .& .&. \\ .& .& .&. \\ .& .& .&. \\ 0& 0 & ... & \sigma_{n}\\ .& . & .& .\\ 0& 0 & ... & 0 \end{pmatrix}\Sigma矩阵为m行*n行的矩阵。

【问题】

 为什么 A^{T}A的特征向量组成的就是我们SVD中的V矩阵,而AA^{T}的特征向量组成的就是我们SVD中的U矩阵?

以V矩阵的证明为例。

A=U\Sigma V^{T}\Rightarrow A^{T}=V\Sigma U^{T}\Rightarrow A^{T}A=V\Sigma U^{T}U\Sigma V^{T}=V\Sigma ^{2}V^{T}

上式证明使用了 U^{T}U=EE^{T}=E 。可以看出 A^{T}A的特征向量组成的的确就是我们SVD中的V矩阵。类似的方法可以得到  的特征向量组成的就是我们SVD中的U矩阵。

进一步我们还可以看出我们的特征值矩阵等于奇异值矩阵的平方,也就是说特征值和奇异值满足如下关系:

\sigma _{I}=\sqrt{\lambda _{i}}

这样也就是说,我们可以不用 \sigma _{i}=Av_{i}/u_{i}来计算奇异值,也可以通过求出A^{T}A的特征值取平方根来求奇异值。

得到V,取V的前k个列向量即可,注意:np.linalg.svd()算法得到的是右奇异矩阵V的转置矩阵。

左奇异矩阵可以用于行数的压缩。

右奇异矩阵可以用于列数即特征维度的压缩,也就是我们的PCA降维。

3、SVD算法的python实现

import numpy as np
dataMat=pd.DataFrame([[10,9],
[2,3],
[1,2],
[7,6.5],
[3,2.5]])

meanVals=np.mean(dataMat,axis=0) #求dataMat各列均值
meanRemoved=dataMat-meanVals #减去原始数据中的均值,避免协方差计算中出现乘以0的情况
U,sigma,VT = np.linalg.svd(meanRemoved)  #svd分解
U
sigma
VT  #矩阵V的转置矩阵

输出:

array([[-0.72111445, -0.08804519,  0.3845965 , -0.3104834 ,  0.47741762],
       [ 0.31370366, -0.47964546, -0.36128726,  0.08356321,  0.73076601],
       [ 0.45920751, -0.2768988 ,  0.8388389 ,  0.09029741,  0.02551835],
       [-0.31690607,  0.0168014 ,  0.07527715,  0.94221156,  0.07654401],
       [ 0.26510935,  0.82778805,  0.11063315,  0.02647648,  0.48118945]])

array([9.65906966, 0.77612712])

array([[-0.78139452, -0.62403734],
       [ 0.62403734, -0.78139452]])

 A=U\Sigma V^{T}\begin{pmatrix} 5.4 & 4.4\\ -2.6&-1.6 \\ -3.6 & -2.6\\ 2.4& 1.9\\ -1.6 & -2.1 \end{pmatrix}=\begin{pmatrix} -0.72111445 & -0.08804519& 0.3845965& -0.3104834&0.47741762\\0.31370366 & -0.47964546& -0.36128726& 0.08356321 &0.73076601 \\ 0.45920751 & -0.2768988&0.8388389 & 0.09029741 &0.02551835\\ -0.31690607 & 0.0168014&0.07527715 & 0.94221156 &0.07654401\\ 0.26510935 & 0.82778805&0.11063315 & 0.02647648 & 0.48118945\end{pmatrix}\begin{pmatrix} 9.65906966 & 0\\ 0&0.77612712 \\0 & 0\\ 0& 0\\ 0 & 0 \end{pmatrix}\begin{pmatrix} -0.78139452 & -0.62403734\\ 0.62403734&-0.78139452 \end{pmatrix}

三、主成分分析的特点

  • 有多少个变量就会有多少个正交的主成分;
  • 主成分的变异(方差)之和等于原始变量的所有变异(方差);
  • 前几大主成分的变异(方差)累计可以解释原多元数据中的多数变异(方差);
  • 原始自变量间有较高线性相关的情况下可以考虑使用主成分分析,如果原始变量不相,则不需要做主成分分析;
  • 主成分分析适用于正态分布的变量,商业分析中右偏型数据可以通过取对数、Tukey正态性打分等方式可转换为接近正态分布。

 

 

 

 

 

 

 

  • 19
    点赞
  • 72
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

xia ge tou lia

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值