主成分分析法(Principal Component Analysis,PCA)是一种用于把高维数据降成低维,使分析变得更加简便的分析方法。比如我们的一个样本可以由维随机变量来刻画,运用主成分分析法,我们可以把这些分量用更少的、这个分量的线性组合来表示。本文多为学习后的个人理解,如有错误还请指出。
我们把降维后的变量称为主成分(Principal Component),设其为(我们并不取这全部的个变量,否则降维就没有意义了)。称为第主成分。每个主成分都是原来个变量的线性组合,即 或者写成更简便的线性代数形式:设则这个关系可被表示为 为了达到降维的目的,我们需要保证
是线性无关的,这要求是正交阵。 如果存在线性相关的关系( 不全为的使得),我们得到的结果中就存在着冗余信息(某个主成分可以由其它主成分表示),这种情况应该被排除。
选出个主成分中能显著代表原本变量的个,来实现对数据的降维。 这一条提出了一个问题:我们按照什么标准来衡量主成分的好坏关系呢?统计学认为,一组数据越分散,它的方差越大,它所包含的信息就越多。(知乎的这个问题讨论了这一点)因此,PCA选出这个主成分中方差最大的个作为新的变量。
设我们的个样本对应的数据矩阵为,每列代表一个样本;我们先将其标准化,使每行的均值为,并消除量纲的影响,便于进一步处理:
其中和分别为第行的均值和样本标准差,记处理后的矩阵为对该矩阵做线性变换的结果为
我们用表示维随机变量的协方差矩阵,那么有
由上一部分的条件,应当有线性不相关,即.而,即的方差(注意跟上文的的记号意义不同,上面的指的是协方差矩阵),因此就是一个对角矩阵,即
由协方差的定义,
由于经我们处理,对任意的都有,故
那么上面的还可表示为
又因为,这里表示维随机变量的协方差矩阵。
由于是实对称矩阵(这一点由协方差矩阵的定义即得:),由实对称矩阵的性质,一定可以相似对角化,即存在正交矩阵使得
其中为的个特征值。
这时候我们取出上面得到的式子
由基本思想部分的前提,应当是正交矩阵,于是我们得到
由线性代数定理(若阶方阵与对角矩阵相似,则对角线上的元素即为的个特征值),即更巧妙的是,我们可以求得变换矩阵,即将相似对角化所需的正交阵之转置。得到了这个变换矩阵,我们就能得到个变换后的主成分了。具体地说,若
则第个主成分前面的系数即的第行,的第列,正好是的第个特征值对应的特征向量(指的是相似对角化矩阵中标准的特征向量)。
经过上面一大串推导,我们得到如下结论:
若是标准化处理过的数据矩阵,那么的个特征值即为线性变换后的个随机变量(即我们提到的主成分)的方差;线性变换所需矩阵,其中为将相似对角化所需的正交阵(由线性代数知识,这也是的个特征向量组成的矩阵)。
这个结论使我们能够求出线性变换所需要的矩阵;此外我们可以根据特征值将个主成分排序,求得方差最大的个主成分。更具体地,求排序后的个主成分的算法如下:
举个例子:(数据来自https://zhuanlan.zhihu.com/p/454447043)
我们要将如下的数据中5个变量(设能力,品格,担保,资本,环境为)降维:
首先我们写出它对应的原始数据矩阵:
然后将其标准化:
求出的协方差矩阵:
求出的特征值(从大到小排列)和特征向量:
我们怎么确定最终取出几个主成分呢?一般认为当取出的个主成分方差贡献比例之和达到时就可以较好地代替原来的个变量了。因此我们还需要计算每个特征值所对应的方差贡献比例。由于特征值的方差占比即按照上述方法,计算方差占比如下:
特征值的方差贡献率,累计方差贡献率 特征值的方差贡献率,累计方差贡献率(已到达) 特征值的方差贡献率,累计方差贡献率 特征值的方差贡献率,累计方差贡献率 特征值的方差贡献率,累计方差贡献率
可以看到前个特征值对应的主成分即达到了的方差贡献率,因此我们可以把原来的个变量“浓缩”表示为下面的个主成分(系数即为特征值对应的特征向量的各分量):
这样,我们就成功实现了降维,以后我们分析数据的时候就可以用和来代替原来的五个指标了。
下面我们用python来实现一下上述的过程。
import numpy as np
from numpy import linalg
class PCA:
'''
dataset 形如array([样本1,样本2,...,样本m]),每个样本是一个n维的ndarray
'''
def __init__(self, dataset):
# 这里的参数跟上文是反着来的(每行是一个样本),需要转置一下
self.dataset = np.matrix(dataset, dtype='float64').T
'''
求主成分;
threshold可选参数表示方差累计达到threshold后就不再取后面的特征向量.
'''
def principal_comps(self, threshold = 0.85):
# 返回满足要求的特征向量
ret = []
data = []
# 标准化
for (index, line) in enumerate(self.dataset):
self.dataset[index] -= np.mean(line)
# np.std(line, ddof = 1)即样本标准差(分母为n - 1)
self.dataset[index] /= np.std(line, ddof = 1)
# 求协方差矩阵
Cov = np.cov(self.dataset)
# 求特征值和特征向量
eigs, vectors = linalg.eig(Cov)
# 第i个特征向量是第i列,为了便于观察将其转置一下
for i in range(len(eigs)):
data.append((eigs[i], vectors[:, i].T))
# 按照特征值从大到小排序
data.sort(key = lambda x: x[0], reverse = True)
sum = 0
for comp in data:
sum += comp[0] / np.sum(eigs)
ret.append(
tuple(map(
# 保留5位小数
lambda x: np.round(x, 5),
# 特征向量、方差贡献率、累计方差贡献率
(comp[1], comp[0] / np.sum(eigs), sum)
))
)
print('特征值:', comp[0], '特征向量:', ret[-1][0], '方差贡献率:', ret[-1][1], '累计方差贡献率:', ret[-1][2])
if sum > threshold:
return ret
return ret
测试一下刚才的例子:
p = PCA(
[[66, 64, 65, 65, 65],
[65, 63, 63, 65, 64],
[57, 58, 63, 59, 66],
[67, 69, 65, 68, 64],
[61, 61, 62, 62, 63],
[64, 65, 63, 63, 63],
[64, 63, 63, 63, 64],
[63, 63, 63, 63, 63],
[65, 64, 65, 66, 64],
[67, 69, 69, 68, 67],
[62, 63, 65, 64, 64],
[68, 67, 65, 67, 65],
[65, 65, 66, 65, 64],
[62, 63, 64, 62, 66],
[64, 66, 66, 65, 67]]
)
lst = p.principal_comps()
print(lst)
输出结果:
# 这部分是运行时输出的
特征值: 3.4531784074578318 特征向量: [0.48198 0.51227 0.45384 0.51336 0.18914] 方差贡献率: 0.69064 累计方差贡献率: 0.69064
特征值: 1.2230892804516 特征向量: [ 0.33297 0.13247 -0.39212 0.20476 -0.82213] 方差贡献率: 0.24462 累计方差贡献率: 0.93525
# 这部分是返回的结果,为了美观稍微调整了一下格式
[
(array([0.48198, 0.51227, 0.45384, 0.51336, 0.18914]), 0.69064, 0.69064),
(array([ 0.33297, 0.13247, -0.39212, 0.20476, -0.82213]), 0.24462, 0.93525)
]
我这里设置的返回结果是一个三元组,如果有需要也可以自己调整一下返回的结果格式。此外,通过调整threshold可选参数可以设置累计方差贡献率到达多少时停止取主成分向量(默认为0.85)。
这种经典的算法sklearn库也已经帮我们实现好了。对于上面的例子,其等价的使用sklearn库的代码如下:
import numpy as np
from sklearn.decomposition import PCA
X = np.array(
[[66, 64, 65, 65, 65],
[65, 63, 63, 65, 64],
[57, 58, 63, 59, 66],
[67, 69, 65, 68, 64],
[61, 61, 62, 62, 63],
[64, 65, 63, 63, 63],
[64, 63, 63, 63, 64],
[63, 63, 63, 63, 63],
[65, 64, 65, 66, 64],
[67, 69, 69, 68, 67],
[62, 63, 65, 64, 64],
[68, 67, 65, 67, 65],
[65, 65, 66, 65, 64],
[62, 63, 64, 62, 66],
[64, 66, 66, 65, 67]]
)
# n_components 指明了降到几维
pca = PCA(n_components = 2)
# 利用数据训练模型(即上述得出特征向量的过程)
pca.fit(X)
# 得出原始数据的降维后的结果;也可以以新的数据作为参数,得到降维结果。
print(pca.transform(X))
# 打印各主成分的方差占比
print(pca.explained_variance_ratio_)
运行结果:
[[-1.51394918 -0.21382815]
[ 0.25137676 -1.8134245 ]
[10.61577071 2.68155382]
[-6.48520841 -1.16575919]
[ 5.53026102 -1.52083322]
[ 0.70154125 -1.8544697 ]
[ 1.82460091 -1.29624147]
[ 2.44281085 -1.60484093]
[-1.40146605 -0.59189041]
[-7.76925956 3.34817657]
[ 1.8850487 0.61749314]
[-5.41819247 -0.9163256 ]
[-1.764172 0.155228 ]
[ 3.06230672 1.51679123]
[-1.96146925 2.65837042]]
# 下面是方差贡献率
[0.82399563 0.11748567]
我们发现这个方差贡献率(也就是特征值占比)跟我们手写的很不一样。(手写的是[0.69064, 0.24462])。这里可能会出现分歧的地方就是是否对原始数据除以样本标准差。当我把手写代码部分的
self.dataset[index] /= np.std(line, ddof = 1)
这一行注释掉后,发现运行结果与sklearn库的基本一致了:
特征值: 22.075235372070864 特征向量: [0.56177 0.58975 0.27868 0.50573 0.05644] 方差贡献率: 0.824 累计方差贡献率: 0.824
特征值: 3.1474971323292125 特征向量: [ 0.37826 -0.06431 -0.61334 0.06946 -0.68686] 方差贡献率: 0.11749 累计方差贡献率: 0.94148
因此可以得出sklearn库在训练时似乎没有消除量纲,即没有对数据除以其样本标准差。当然,这仅仅是个人理解,不过与sklearn库结果基本吻合大致上印证了这个猜想。在一篇文章里我找到了关于量纲的讨论:若各个变量的单位一致,则各个属性是可以比较的,此时可以直接求协方差;但当各个变量单位不同时(如身高/体重),这时不同变量之间没有可比性,这时就应该消除量纲的影响(即除以样本标准差)。
1. 2.https://www.zhihu.com/question/274997106/answer/1055696026 3.https://zhuanlan.zhihu.com/p/454447043 4.https://www.jianshu.com/p/c21c0e2c403a