主成分分析
- 一、概述
- 1.1 问题提出
- 1.2 降维的作用
- 二、主成分分析(PCA)主要思想
- 三、相关数学知识
- 四、PCA实现步骤
- 4.1 特征值分解矩阵
- 4.2 SVD分解协方差矩阵
- 五、python程序实现
- 5.1 利用数学公式实现
- 5.2 使用sklearn实现
一、概述
1.1 问题提出
在实际问题研究中,多变量问题是经常会遇到的。变量太多,无疑会增加分析问题的难度与复杂性,而且在许多实际问题中,多个变量之间是具有一定的相关关系的。因此,人们会很自然地想到,能否在相关分析的基础上, 用较少的新变量代替原来较多的旧变量,而且使这些较少的新变量尽可能多地保留原来变量所反映的信息?事实上,这种想法是可以实现的,主成分分析方法就是综合处理这种问题的一种强有力的工具。主成分分析是把原来多个变量划为少数几个综合指标的一种统计分析方法。从数学角度来看,这是一种降维处理技术。
1.2 降维的作用
降维是将高维度的数据(指标太多)保留下最重要的一些特征,去除噪声和不重要的特征,从而实现提升数据处理速度的目的。在实际的生产和应用中,降维在一定的信息损失范围内,可以为我们节省大量的时间和成本。降维也成为应用非常广泛的数据预处理方法。
降维具有如下一些优点:
- 使得数据集更易使用;
- 降低算法的计算开销;
- 去除噪声;
- 使得结果容易理解。
降维的算法有很多,比如奇异值分解(SVD)、主成分分析(PCA)、因子分析(FA)、独立成分分析(ICA)。
二、主成分分析(PCA)主要思想
PCA的主要思想是将n维特征映射到维上,这维是全新的正交特征也被称为主成分,是在原有维特征的基础上重新构造出来的维特征。
PCA的工作就是从原始的空间中顺序地找一组相互正交的坐标轴,新的坐标轴的选择与数据本身是密切相关的。其中,第一个新坐标轴选择是原始数据中方差最大
的方向,第二个新坐标轴选取是与第一个坐标轴正交的平面中使方差最大
的,第三个轴是与第1,2个轴正交的平面中方差最大
的。依次类推,可以得到个这样的坐标轴。通过这种方式获得的新的坐标轴,我们发现,大部分方差都包含在前面个坐标轴中,后面的坐标轴所含的方差几乎为0。于是,我们可以忽略余下的坐标轴,只保留前面个含有绝大部分方差的坐标轴。事实上,这相当于只保留包含绝大部分方差的维度特征,而忽略包含方差几乎为0的特征维度,实现对数据特征的降维处理。
思考:我们如何得到这些包含最大差异性的主成分方向呢?
答案:事实上,通过计算数据矩阵的协方差矩阵,然后得到协方差矩阵的特征值特征向量,选择特征值最大(即方差最大)的k个特征所对应的特征向量组成的矩阵。这样就可以将数据矩阵转换到新的空间当中,实现数据特征的降维。
由于得到协方差矩阵的特征值特征向量有两种方法:特征值分解协方差矩阵、奇异值分解协方差矩阵,所以PCA算法有两种实现方法:基于特征值分解协方差矩阵实现PCA算法
、基于SVD分解协方差矩阵实现PCA算法
。
三、相关数学知识
方差:
最先提到的一个概念,也是旋转坐标轴的依据。之所以使用方差作为旋转条件是因为:最大方差给出了数据的最重要的信息。
或
由方差可推出标准差:
协方差:
用来衡量两个变量的总体误差,方差是协方差的一种特殊情况,即当两个变量相同。可以通俗的理解为:两个变量在变化过程中的方向,度量各个维度偏离其均值的程度,定义为:
由协方差的定义可以推出两个性质:
协方差矩阵:
协方差只能处理二维问题(即两个特征X、Y),维数多了自然需要计算多个协方差矩阵,以三维数据为例,协方差矩阵表示为:
可见,协方差矩阵是一个对称的矩阵,且对角线是各维度上的方差。正是由于协方差矩阵为对称矩阵所以矩阵分解后特征值所对应的特征向量 一定无线性关系,且相互之间一定正交,即内积为零。
方差和协方差的除数是n-1,这是为了得到方差和协方差的无偏估计。
特征值与特征向量:
- 特征值与特征向量
如果一个向量 是矩阵的特征向量,将一定可以表示成下面的形式:
其中, 是特征向量对应的特征值,一个矩阵的一组特征向量是一组正交向量(实对称阵特征向量正交
)。 - 特征值分解矩阵
对于矩阵 ,有一组特征向量 ,将这组向量进行正交化单位化,就能得到一组正交单位向量。特 征值分解,就是将矩阵分解为如下式:
其中, 是矩阵 的特征向量组成的矩阵,
SVD分解矩阵:
奇异值分解是一个能适用于任意矩阵的一种分解的方法,对于任意矩阵总是存在一个奇异值分解:
假设 是一个 的矩阵,那么得到的 是一个 的方阵,U里面的正交向量被称为左奇异向 转置矩阵,是一个
- 求 的特征值和特征向量,用单位化的特征向量构成
- 求 的特征值和特征向量,用单位化的特征向量构成
- 将 或者 的特征值求平方根,然后构成 。
四、PCA实现步骤
4.1 特征值分解矩阵
输入:数据集 ,需要降到
- 去平均值(即去中心化),即每一位特征减去各自的平均值。
- 计算协方差矩阵 (
这里除或不除样本数量n或n-1,其实对求出的特征向量没有影响
) - 用特征值分解方法求协方差矩阵
- 对特征值从大到小排序,选择其中最大的 个。然后将其对应的 个特征向量分别作为行向量组 成特征向量矩阵
- 将数据转换到 个特征向量构建的新空间中,即
为什么数据要减去均值?
把数据中心化,减少过度拟合的可能性。
为什么选取前k个最大的特征值?
k维特征是将n维样本点转换为k维后,每一维上的样本方差都很大。特征值越大,说明矩阵在对应的特征向量上的方差越大,样本点越离散,越容易区分,信息量也就越多。因此,特征值最大的对应的特征向量方向上所包含的信息量就越多,如果某几个特征值很小,那么就说明在该方向的信息量非常少,我们就可以删除小特征值对应方向的数据,只保留大特征值方向对应的数据,这样做以后数据量减小,但有用的信息量都保留下来了。
在n维中选多少个k最合适?
研究表明,所选的主轴总长度占所有主轴长度之和的大约 即可如下图公式。
为样本数,
计算过程:
此处参考数学建模清风老师课件
假设有 个样本, 个指标, 则可构成大小为 的样本矩阵 :
- 我们首先对其进行去中心:
按列计算均值 计算得 去中心化数据 , 原始样本矩阵经过去中心化变为: - 计算处理样本的协方差矩阵
其中
(注意: 如果进行标准化处理,上面 12 两步可直接合并为一步:直接计算 矩阵的样本相关系数矩阵) - 计算 的特征值和特征向量
特征值: 是半正定矩阵, 且
特征向量: - 计算主成分贡献率以及累计贡献率
贡献率 累计贡献率 - 写出主成分
一般取累计贡献率超过 的特征值所对应的第一、第二、 第 个主成分。
第 个主成分: - 根据系数分析主成分代表的意义
对于某个主成分而言, 指标前面的系数越大, 代表该指标对于该主成分的影响越大。 - 利用主成分的结果进行后续的分析
(1) 主成分得分:可用于评价类模型吗? (千万别用!!!)
(2) 主成分可用于聚类分析(方便画图)
(3) 主成分可用于回归分析
4.2 SVD分解协方差矩阵
输入: 数据集 ,需要降到维。
- 去平均值,即每一位特征減去各自的平均值。
- 计算协方差矩阵。
- 通过SVD计算协方差矩阵的特征值与特征向量。
- 对特征值从大到小排序,选择其中最大的k个。然后将其对应的k个特征向量分别作为列向量組 成特征向量矩阵。
- 将数据转换到
在PCA降维中,我们需要找到样本协方差矩阵 的最大 个特征向量,然后用这最大的 个特 征向量组成的矩阵来做低维投影降维。可以看出,在这个过程中需要先求出协方差矩阵 , 当样本数多、样本特征数也多的时候,这个计算还是很大的。当我们用到SVD分解协方差矩阵的时候,SVD有两个好处:
- 有一些SVD的实现算法可以先不求出协方差矩阵 也能求出我们的右奇异矩阵 。也就是说,我们的PCA算法可以不用做特征分解而是通过SVD来完成,这个方法在样本量很大的时候很有效。实际上,
scikit-learn
的PCA算法的背后真正的实现就是用的SVD,而不是特征值分解。 - 注意到PCA仅仅使用了我们SVD的左奇异矩阵,没有使用到右奇异值矩阵,那么右奇异值矩阵有 什么用呢?
假设我们的样本是 的矩阵X,如果我们通过SVD找到了矩阵 最大的 个特征向量组成的 的矩阵 ,则我们可以做如下处理:
可以得到一个 的矩阵 , , 文个矩阵和我们原来 的矩阵 相比,列数从 減到了
五、python程序实现
5.1 利用数学公式实现
- 部分数据
32.502345269453031,31.70700584656992
53.426804033275019,68.77759598163891
61.530358025636438,62.562382297945803
47.475639634786098,71.546632233567777
59.813207869512318,87.230925133687393
55.142188413943821,78.211518270799232
52.211796692214001,79.64197304980874
39.299566694317065,59.171489321869508
48.10504169176825,75.331242297063056
52.550014442733818,71.300879886850353
45.419730144973755,55.165677145959123
54.351634881228918,82.478846757497919
44.164049496773352,62.008923245725825
58.16847071685779,75.392870425994957
56.727208057096611,81.43619215887864
48.955888566093719,60.723602440673965
44.687196231480904,82.892503731453715
60.297326851333466,97.379896862166078
45.618643772955828,48.847153317355072
38.816817537445637,56.877213186268506
66.189816606752601,83.878564664602763
65.41605174513407,118.59121730252249
47.48120860786787,57.251819462268969
41.57564261748702,51.391744079832307
51.84518690563943,75.380651665312357
59.370822011089523,74.765564032151374
57.31000343834809,95.455052922574737
63.615561251453308,95.229366017555307
- 第一种实现方式(去中心化)
import numpy as np
import matplotlib.pyplot as plt
data = np.genfromtxt("data.csv", delimiter=',')
x_data = data[:, 0]
y_data = data[:, 1]
plt.scatter(x_data, y_data)
# 数据中心化
def zeroMean(dataMat):
# 按列求平均,及几个特征的平均
meanval = np.mean(dataMat, axis=0)
newData = dataMat - meanval
return newData, meanval
newData, meanval = zeroMean(data)
print(newData)
# 协方差矩阵生成,n个变量,生成n*n矩阵
covMat = np.cov(newData, rowvar=False)
# 计算特征值、特征向量
eigvals, eigvects = np.linalg.eig(covMat)
print("特征值:\n",eigvals,"\n特征向量:\n",eigvects)
# 对特征值从大到小排序 np.argsort(a)从小到大排序 -a为从大到小排序
eigvalindice = np.argsort(-eigvals)
# 最大的1个特征值的下标
n_eigvalindice = eigvalindice[:1] # 返回[1]
# 最大的特征值对应的特征向量 n_eigVect[1]
n_eigVect = eigvects[:, n_eigvalindice]
print("最大特征值对应特征向量:\n",n_eigVect)
# 低维特征空间的数据
# np.dot()函数主要有两个功能,向量点积和矩阵乘法
lowDDataMat = np.dot(newData, n_eigVect)
print(lowDDataMat,lowDDataMat.shape)
# 利用低纬度数据来重构数据
reconMat = np.dot(lowDDataMat, n_eigVect.T) +meanval # (100,2)
x_data = np.array(reconMat)[:, 0]
y_data = np.array(reconMat)[:, 1]
plt.scatter(x_data, y_data, c='b')
plt.savefig("PCA降维.png")
plt.show()
结果显示:
- 第二种实现方式(标准化)
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn import preprocessing
# 导入数据集
data = np.genfromtxt("data_1.csv", delimiter=',')
data = data[1:, :]
print(data)
# 标准化
new_data = preprocessing.scale(data)
# 求相关系数矩阵 相关系数为行变量
covX = np.around(np.corrcoef(new_data.T), decimals=3)
# 求解特征值、特征变量
featValue, featVec = np.linalg.eig(covX)
print("特征值:\n",featValue,"\n特征向量:\n",featVec)
# 降序排列
# featValue = sorted(featValue, reverse=True)
# [::-1] 顺序相反操作
featValue = sorted(featValue)[::-1]
print("特征值降序排列:\n",featValue)
# 绘制散点图与折线图
plt.scatter(range(1, data.shape[1] + 1), featValue)
plt.plot(range(1, data.shape[1] + 1), featValue)
plt.title("Scree Plot")
plt.xlabel("Factors")
plt.ylabel("Eigenvalue")
plt.savefig("Factors")
plt.grid()
plt.show()
# 求特征值的贡献度
gx = featValue / np.sum(featValue)
print("特征值的贡献度:\n",gx)
# 累计贡献度
lg = np.cumsum(gx)
print("累计贡献度\n",lg)
plt.scatter(range(1, data.shape[1] + 1), lg)
plt.plot(range(1, data.shape[1] + 1), lg)
plt.axhline(y=0.85, c='r', xmin=0.1, xmax=0.9)
plt.text(1,0.86,"水平线")
plt.savefig("累计贡献度")
# 选出主成分对应下标
k = [i for i in range(len(lg)) if lg[i] < 0.85]
print("主成分对应下标\n",k)
# 选出主成分对应的特征向量
Vec = featVec.T[k].T
print("主成分对应的特征向量\n",Vec)
# 降维新数据
finalData = np.dot(new_data, Vec)
print("低维数据:\n",finalData)
plt.figure(figsize=(14, 14))
ax = sns.heatmap(Vec, annot=True, cmap="BuPu")
# 设置y轴字体大小
ax.yaxis.set_tick_params(labelsize=15)
plt.title("Factor Analysis", fontsize="xx-large")
# 设置y轴标签
plt.ylabel("Sepal Width", fontsize="xx-large")
# 保存图片
plt.savefig("factorAnalysis", dpi=500)
# 显示图片
plt.show()
特征值:
累计贡献度:
特征向量:
5.2 使用sklearn实现
一、参数说明 (Parameters)
sklearn.decomposition.PCA(n_components=None, copy=True, whiten=False)
- n_components: int, float, None or str
意义 :代表返回的主成分的个数,也就是你想把数据降到几维
n_components 代表返回前 2 个主成分
n_components <1代表满足最低的主成分方差男计贡献率
n_components ,指返回满足主成分方差男计贡献率达到 的主成分
n_components=None,返回所有主成分
components=‘mle’,将自动选取主成分个数 - copy: bool类型, False/True 默认是True
意义:在运行的过程中,是否将原数据复制。由于你在运行的过程中,是在降维,数据会变动。
这copy主要影响的是,调用显示降维后的数据的方法不同。
copy=True时,直接 fit_transform (X),就能够显示出降维后的数据。
copy=False时,需要 fit (X).transform (X) ,才能够显示出降维后的数据。
(fit_transform()方法后面会讲到!) - whiten: bool类型,False/True 默认是False
意义:白化。白化是一种重要的预处理过程,其目的就是降低输入数据的冗余性,使得经过白化处理的输入数据具有如下性质:(i)特征之间相关性较低;(ii)所有特征具有相同的方差。 - svd_solver: str类型,str {‘auto’, ‘full’, ‘arpack’, ‘randomized’}
意义:定奇异值分解 SVD 的方法。
svd_solver=auto:PCA 类自动选择下述三种算法权衡。
svd_solver= ‘full’:传统意义上的 SVD,使用了 scipy 库对应的实现。
svd_solver=‘arpack’:直接使用 scipy 库的 sparse SVD 实现,和 randomized 的适用场景类似。
svd_solver=‘randomized’:适用于数据量大,数据维度多同时主成分数目比例又较低的 PCA 降维。
二、属性 (Attributes)
- components_: 返回最大方差的主成分。
- explained_variance_: 它代表降维后的各主成分的方差值。方差值越大,则说明越是重要的主成分。
- explained_variance_ratio_: 它代表降维后的各主成分的方差值占总方差值的比例,这个比例越大,则越是重要的主成分。(主成分方差贡献率)
- singular_values_: 返回所被选主成分的奇异值。
实现降维的过程中,有两个方法,一种是用特征值分解,另一种用奇异值分解,前者限制比较多,需要矩阵是方阵,而后者可以是任意矩
阵,而且计算量比前者少,所以说一般实现 - mean_: 每个特征的经验平均值,由训练集估计。
- n_features_: 训练数据中的特征数。
- n_samples_: 训练数据中的样本数量。
- noise_variance_: 噪声协方差
三、方法 (Methods)
- fit(self, X,Y=None) 模型训练,由于PCA是无监督学习,所以Y=None,没有标签
- fit_transform(self, :将模型与X进行训练,并对X进行降维处理,返回的是降维后的数据。
- get_covariance(self) 获得协方差数据
- get_params(self,deep 返回模型的参数
- get_precision(self) 计算数据精度矩阵 ( 用生成模型)
- inverse_transform(self, X) 将降维后的数据转换成原始数据,但可能不会完全一样
- score(self, X, Y= None) 计算所有样本的log似然平均值
- transform(X) 将数据X转换成降维后的数据。当模型训练好后,对于新输入的数据,都可以用transform方法来降维。
# 用sklearn的PCA
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
import numpy as np
data = np.genfromtxt("data.csv",delimiter=",")
plt.scatter(data[:,0],data[:,1],c="b")
plt.show()
# pca = PCA(n_components=2,whiten=True) whiten=True标准化
pca=PCA(n_components=1)
pca.fit(data)
main_var = pca.explained_variance_ # 特征值
print(main_var)
pca_data = pca.transform(data)
print(pca_data)
sklearn降维部分结果显示:
[-44.02694787]
[ -1.49722533]
[ -3.35564513]
[ -1.73205523]
[ 17.84406034]
[ 7.68710859]
[ 7.6311404 ]
[-16.4703207 ]
[ 1.92574891]
[ 0.35289859]
[-17.26071108]
[ 11.13030667]
[-11.73358623]
[ 6.54975245]
[ 11.27989566]
[-10.70315359]
[ 7.11092921]
[ 27.10646295]
[-22.80011879]
sklearn中的PCA是通过svd_flip函数实现的,sklearn对奇异值分解结果进行了一个处理,因为,也就是和同时取反得到的结果是一样的,而这会导致通过PCA降维得到不一样的结果(虽然都是正确的)。
本文为学习过程记录,难免出现错误,欢迎指正。