numpy实现鸢尾花数据集PCA降维

PCA降维过程

在前面的一篇博客中我已经从数学角度解释了PCA降维的原理，我们从中也可以得到PCA降维的过程
1）将原始数据做转置运算，每一行代表一个维度
2）每一行（代表一个属性字段）进行零均值化，即减去这一行的均值
3）得到原始数据的协方差矩阵
4）求出协方差矩阵的特征值及对应的特征向量的单位向量
5）将特征向量按对应特征值大小从上到下按行排列成矩阵，取前k行组成矩阵P
6）用上面得到矩阵P和标准化后数据相乘，即可得到降维到k维后的数据

数据集

数据集我们使用sklearn库中的iris数据集，数据集中的每个样本有4个特征参数

原始数据格式：
5.1, 3.5, 1.4, 0.2
4.9, 3.0, 1.4, 0.2
4.7, 3.2, 1.3, 0.2
4.6, 3.1, 1.5, 0.2
5.0, 3.6, 1.4, 0.2
5.4, 3.9, 1.7, 0.4
4.6, 3.4, 1.4, 0.3
5.0, 3.4, 1.5, 0.2
4.4, 2.9, 1.4, 0.2
4.9, 3.1, 1.5, 0.1

numpy实现PCA降维

加载数据集

数据我们使用sklearn库中的iris中的数据集，所以需要导入sklearn库

   data = datasets.load_iris()["data"]
   print(data)

我们看一下加载进来的数据

[[ 5.1  3.5  1.4  0.2]
 [ 4.9  3.   1.4  0.2]
 [ 4.7  3.2  1.3  0.2]
 [ 4.6  3.1  1.5  0.2]
 [ 5.   3.6  1.4  0.2]
 [ 5.4  3.9  1.7  0.4]
 ...
  [ 6.7  3.3  5.7  2.5]
 [ 6.7  3.   5.2  2.3]
 [ 6.3  2.5  5.   1.9]
 [ 6.5  3.   5.2  2. ]
 [ 6.2  3.4  5.4  2.3]
 [ 5.9  3.   5.1  1.8]]

数据标准化

原始数据每一列是同一个维度特征，在标准化时候我们需要的也是对维度进行数据标准化处理，所以需要按列取数据,另外因为后面我们需要计算协方差，所以对数据进行标准化方式是去均值

 #axis = 0,按列取值求均值
 mean_vector=np.mean(data,axis=0)
 print("均值向量为：%s\n标准化数据：%s"% (mean_vector,data - mean_vector))

标准化后的数据为：

均值向量为：[ 5.84333333  3.054       3.75866667  1.19866667] 
标准化数据：
[[ -7.43333333e-01   4.46000000e-01  -2.35866667e+00  -9.98666667e-01]
 [ -9.43333333e-01  -5.40000000e-02  -2.35866667e+00  -9.98666667e-01]
 [ -1.14333333e+00   1.46000000e-01  -2.45866667e+00  -9.98666667e-01]
 [ -1.24333333e+00   4.60000000e-02  -2.25866667e+00  -9.98666667e-01]
 [ -8.43333333e-01   5.46000000e-01  -2.35866667e+00  -9.98666667e-01]
 ......
  [  4.56666667e-01  -5.54000000e-01   1.24133333e+00   7.01333333e-01]
 [  6.56666667e-01  -5.40000000e-02   1.44133333e+00   8.01333333e-01]
 [  3.56666667e-01   3.46000000e-01   1.64133333e+00   1.10133333e+00]
 [  5.66666667e-02  -5.40000000e-02   1.34133333e+00   6.01333333e-01]]

协方差矩阵

协方差计算的是不同维度之间的协方差，不是不同样本的，首先我们要确定我们数据维度是按行还是列，我们的数据是一列是一个维度，并且要对标准化后的数据进行求协方差，这里我们使用numpy中的 cov() 函数

 # rowvar=0表示数据的每一列代表一个维度
 cov_mat =np.cov(standData,rowvar=0)
print("协方差矩阵:\n%s"%cov_mat)

原始数据中一共四个维度，所以我们得到的协方差矩阵是一个4*4的矩阵

协方差矩阵:
[[ 0.68569351 -0.03926846  1.27368233  0.5169038 ]
[-0.03926846  0.18800403 -0.32171275 -0.11798121]
[ 1.27368233 -0.32171275  3.11317942  1.29638747]
[ 0.5169038  -0.11798121  1.29638747  0.58241432]]

协方差矩阵的特征值和对应的特征向量的计算

在另一篇博客中也说了协方差矩阵是实对称方阵，所以是一定可以做矩阵分解的
在numpy的库函数中，linalg.eig()可以用来计算计算特征值和对应的特征向量，输入参数是一个方阵，得到两个值：w，v
w：特征值。每个特征值根据它的多重性重复。这个数组将是复杂类型，除非虚数部分为0。当传进的参数a是实数时，得到的特征值是实数
v：特征向量。每个特征值对应的特征向量，相同特征值可能对应的特征向量不同

fvalue,fvector = np.linalg.eig(covMat)
print("特征值为：%s\n特征向量为：\n%s" %(fvalue,fvector))

可以得到四个特征值，和四个特征向量的单位向量

特征值为：[ 4.22484077  0.24224357  0.07852391  0.02368303]
特征向量为：
[[ 0.36158968 -0.65653988 -0.58099728  0.31725455]
[-0.08226889 -0.72971237  0.59641809 -0.32409435]
[ 0.85657211  0.1757674   0.07252408 -0.47971899]
[ 0.35884393  0.07470647  0.54906091  0.75112056]]

向量矩阵

前面我们已经得到协方差矩阵的特征值和特征向量，接下来需要按照特征值的大小

#argsort函数返回的是数组值排序的索引值
#将特征值按从大到小顺序排列的索引值
fvaluesort=np.argsort(-fvalue)
print(fvaluesort)
#特征值[ 4.22484077  0.24224357  0.07852391  0.02368303]
#output :[0 1 2 3]

argsort函数返回的是排序的索引值，默认是从小到大
eg：
x = np.array([3, 1, 2])
np.argsort(x)
out： array([1, 2, 0])
其他用法可以参考np.argsort的用法
我们假设现在要将数据降到二维，那么就取最大的两个特征值对应的特征向量

#先获取特征值较大的两个索引以便于下一步得到对应的特征向量
fValueTopN = fValueSort[:2 ]
#得到较大的特征值对应的特征向量
newdata = fVector[:,fValueTopN]
print(newdata)
#output：
[[ 0.36158968 -0.65653988]
 [-0.08226889 -0.72971237]
 [ 0.85657211  0.1757674 ]
 [ 0.35884393  0.07470647]]

我们和上面得到的特征向量和特征值做个比较，确定拿到的是较大的两个特征值对应的特征向量

降维后数据

拿上面得到的前2个较大特征值对应的特征向量和标准化后数据做内积运算，便可得到降维后的数据

newdata = np.dot(data,vectorMatrix)
print(newdata)
#output
降维后的数据：
[[-2.68412563 -0.31939725]
[-2.71414169  0.17700123]
[-2.88899057  0.14494943]
[-2.74534286  0.31829898]
[-2.72871654 -0.32675451]
.......
[ 1.52716661  0.37531698]
[ 1.76434572 -0.07885885]
[ 1.90094161 -0.11662796]
[ 1.39018886  0.28266094]]

至此，我们的降维过程已经完成，我们再试着恢复数据，看看和原始数据有没有很大的差别

print("最终重构结果为:\n{}".format(np.mat(newdata) * fvectormat.T + mean_vector))
#output
复原的数据：
[[5.08303897 3.51741393 1.40321372 0.21353169]
[4.7462619  3.15749994 1.46356177 0.24024592]
[4.70411871 3.1956816  1.30821697 0.17518015]
[4.6422117  3.05696697 1.46132981 0.23973218]
[5.07175511 3.52655486 1.36373845 0.19699991]
[5.50581049 3.79140823 1.67552816 0.32616959] 
......
  [6.95201347 3.04358556 5.90548444 2.09665999]
[6.91756285 3.07544671 5.77722508 2.04293748]
[6.66904015 3.02994114 5.39094874 1.88173174]
[6.14880195 2.65421139 5.13134845 1.77482994]
[6.53272206 2.96578609 5.25579114 1.825527  ]
[6.60688475 2.98181821 5.3662607  1.87161698]
[6.16013695 2.73344296 4.99793961 1.71875852]] 
原始数据;
[[ 5.1  3.5  1.4  0.2]
[ 4.9  3.   1.4  0.2]
[ 4.7  3.2  1.3  0.2]
[ 4.6  3.1  1.5  0.2]
[ 5.   3.6  1.4  0.2]
[ 5.4  3.9  1.7  0.4]
...
 [ 6.7  3.3  5.7  2.5]
[ 6.7  3.   5.2  2.3]
[ 6.3  2.5  5.   1.9]
[ 6.5  3.   5.2  2. ]
[ 6.2  3.4  5.4  2.3]
[ 5.9  3.   5.1  1.8]] 

和原始数据做个比较就会发现和原始数据差别很小，因为做降维相对来说会损失信息量，所以很难做到完整还原原始数据

完整代码

from sklearn import datasets
import numpy as np

class PCAtest():
   def __init__(self,k):
       #降到k维
       self.k = k
       
   # 加载鸢尾花数据集中的特征作为PCA的原始数据集 
   def loadIris(self):
       data = datasets.load_iris()["data"]
       return data
   #数据标准化（去均值）
   def stand_data(self,data):
       #axis = 0,按列取值求均值
       mean_vector = np.mean(data,axis=0)
       return mean_vector,data -mean_vector

   # 计算协方差矩阵
   def getCovMat(self,standData):
       # rowvar=0表示数据的每一列代表一个维度
       return np.cov(standData,rowvar=0)

   # 计算协方差矩阵的特征值和特征向量
   def getFValueAndFVector(self,covMat):
       fValue,fVector = np.linalg.eig(covMat)
       return fValue,fVector

   # 得到特征向量矩阵
   def getVectorMatrix(self,fValue,fVector):
       #从大到小排序，并返回排序后的原索引值
       fValueSort = np.argsort(-fValue)
       #print(fValueSort)
       fValueTopN = fValueSort[:self.k]
       #print(fValueTopN)
       return fVector[:,fValueTopN]

   # 得到降维后的数据
   def getResult(self,data,vectorMat):
       return np.dot(data,vectorMat)


if __name__=="__main__":
   pca = PCAtest(2)
   data = pca.loadIris()
   print("原始数据：\n%s"%data)
   (mean_vector ,standdata)= pca.stand_data(data)
   print("均值向量为：%s \n标准化数据：\n%s" % (mean_vector, standdata))
   cov_mat = pca.getCovMat(standdata)
   print("协方差矩阵:\n%s"%cov_mat)
   fvalue,fvector = pca.getFValueAndFVector(cov_mat)
   print("特征值为：%s\n特征向量为：\n%s" %(fvalue,fvector))
   fvectormat = pca.getVectorMatrix(fvalue,fvector)
   print("最终需要的特征向量：\n%s"%fvector)
   newdata = pca.getResult(standdata,fvectormat)
   print("降维后的数据：\n%s"%newdata)
   print("最终重构结果为:\n{}".format(np.mat(newdata) * fvectormat .T + mean_vector))