numpy实战，PCA降维（特征值分解方法，不讲奇异值分解SVD）

PCA介绍

在实践中，获取的数据维度都比较高。因为很多时候会把离散变量使用读入编码弄成多维空间，这样数据可以很稀疏，也会包含一些噪声。此时可以运用PCA降维，使特征之间更加独立，也能去除噪声减小计算量。
PCA（Principle Component Analysis）即主成分分析，不仅可对高维数据降维，更重要的是经过降维，去除了噪声，从而发现数据中的一些固有的模式。
PCA把原先的N个特征数用数目更少的M个特征来代替，新特征是旧特征的线性组合，这些组合最大化了样本方差，尽可能使新的M个特征互不相关。

PCA降维的主要目的：

1. 降低了特征的个数，同时也减少了计算量；
2. 降低了特征之间的相关性，使得特征之间更加独立；
3. 减少噪声数据的影响，使模型更稳定；
4. 方便对数据进行可视化

在工程领域，PCA的实现一般有两种，一种是特征值分解，一种是奇异值分解（SVD）。比如，在scikit-learn中就是使用SVD来实现的。为更好地使用numpy 实现，本文以特征值分解为例，其通常包含以下步骤：

1. 对数据进行零均值化处理；
2. 计算均值化后的数据的协方差矩阵；
3. 计算协方差矩阵的特征值和特征向量；
4. 找出所需的K个特征值及其对应的特征向量；
5. 将数据映射到K个特征向量的空间中，实现降维。

下面就按这几个步骤一步步使用numpy来编码，目的是对sklearn（scikit-learn）中自带的Iris数据集进行降维，Iris数据集是4维的数据，使用PCA将其降到2维，以方便进行数据的可视化。

1.数据均值化

加载sklearn的数据集，便可使用np.mean方法计算每个特征值的均值。
零均值化就是求每个特征值的平均值，然后用该特征的所有数都减去这个均值。也就是说，这里的零均值化是对每一个特征（每一列）而言的。因此在np.mean方法中必须加一个axis=0的参数。零均值化后，每个特征的均值变成0。

# 求矩阵每一列的均值
meanVals = np.mean(dataMat, axis=0)

2.协方差矩阵

协方差（Covariance），意为协同方差，指两个变量之间的关系。如果结果为正值，则说明两者是正相关的，结果为负值则是负相关。如果为0就是统计上说的相互独立。
协方差矩阵是对多维特征而言的，计算两两特征之间的协方差，从而形成一个矩阵。如n个特征的数组，其协方差矩阵为：

从上面公式中可以看出，矩阵是一个对称的方阵，因为cov里的两个变量换一下位置也相同。
要计算数据的协方差，在numpy中可以使用cov函数。参数rowvar很重要，若rowvat=True（默认值），说明每行为一个特征。若rowvar=False，说明每列代表一个特征，因此rowvar需要设置为False。
在零均值化的数据上计算协方差矩阵，代码为cov_mat = np.cov(meaned_data, rowvar=False)sd

3.特征值和向量

定义：A为n阶方阵，若数λ和n维非0列向量x满足Ax=λx，那么数λ称为A的特征值，x称为A的对应于特征值λ的特征向量。式Ax=λx也可写成( A-λE)x=0，并且|λE-A|叫做A 的特征多项式。当特征多项式等于0的时候，称为A的特征方程，特征方程是一个齐次线性方程组，求解特征值的过程其实就是求解特征方程的解。方阵A的所有特征值，叫做A的谱。

numpy中有一个字库linag，专门用于矩阵线性代数运算，计算特征值和特征向量就可以使用其中的eig方法，具体代码如下：

#求特征值和特征向量
eig_vals, eig_vects = np.linalg.eig(np.mat(cov_mat))

#对特征值从大到小排序
sorted_index = np.argsort(-eig_vals)

#取最大的2个特征
topn_index = sorted_index[:2]

#最大的n个特征值对应的特征向量
topn_vects = eig_vects[:, topn_index]

#为了进行可视化，需要将数据降到2维，所以 取前面两个最大的特征值与其对应的特征向量。在取对应向量的时候，必须注意是取前两个列的向量。

question:：python array行向量，列向量搞不清楚

4.数据映射将维

取出最大的两个特征向量后，直接将均值化后的数据与这两个特征向量进行矩阵相乘，即将数据投影到这个低维空间。若要还原进行逆向操作即可，与特征向量的转置矩阵相乘，再加上前面计算出来的每列的均值即可。
代码如下所示：

#投影到低维空间
pca_data = meaned_data * topn_vects

#重构数据，逆向还原
recon_data = (pca_data * topn_vects.T) + mean_val

自此，一个完整的PCA流程就完成了。将原始的4维数据降低到了2维数据，就可以对整个数据进行可视化了。
降维后的2维特征并不是从原来的4维中选择了主要的2维出来，而是映射到了新的二维空间中，所以有一定的信息丢失，当然逆向还原也有误差，与深度学习中的自编码器非常相似。

摘一段人家完整步骤的代码

import numpy as np
import matplotlib.pyplot as plt

"""
函数说明：解析文本数据

Parameters:
filename - 文件名
delim - 每一行不同特征数据之间的分隔方式，默认是tab键‘\t’

Returns:
j将float型数据值列表转化为矩阵返回

"""
def loadDataSet(filename, delim='\t'):
    fr = open(filename)
    stringArr = [line.strip().split(delim) for line in fr.readlines()]
    datArr = [list(map(float, line)) for line in stringArr]
    return np.mat(datArr)


"""
函数说明：PCA特征维度压缩函数

Parameters:
dataMat - 数据集数据
topNfeat - 需要保留的特征维度，即要压缩成的维度数，默认4096

Returns:
lowDDataMat - 压缩后的数据矩阵
reconMat - 压缩后的数据矩阵反构出原始数据矩阵

"""
def pca(dataMat, topNfeat=4096):
    # 求矩阵每一列的均值
    meanVals = np.mean(dataMat, axis=0)
    # 数据矩阵每一列特征减去该列特征均值
    meanRemoved = dataMat - meanVals
    # 计算协方差矩阵，处以n-1是为了得到协方差的无偏估计
    # cov(x, 0) = cov(x)除数是n-1(n为样本个数)
    # cov(x, 1)除数是n
    covMat = np.cov(meanRemoved, rowvar=0)
    # 计算协方差矩阵的特征值及对应的特征向量
    # 均保存在相应的矩阵中
    eigVals, eigVects = np.linalg.eig(np.mat(covMat))
    # sort():对特征值矩阵排序(由小到大)
    # argsort():对特征矩阵进行由小到大排序，返回对应排序后的索引
    eigValInd = np.argsort(eigVals)
    # 从排序后的矩阵最后一个开始自下而上选取最大的N个特征值，返回其对应的索引
    eigValInd = eigValInd[: -(topNfeat+1): -1]
    # 将特征值最大的N个特征值对应索引的特征向量提取出来，组成压缩矩阵
    redEigVects = eigVects[:, eigValInd]
    # 将去除均值后的矩阵*压缩矩阵，转换到新的空间，使维度降低为N
    lowDDataMat = meanRemoved * redEigVects
    # 利用降维后的矩阵反构出原数据矩阵(用作测试，可跟未压缩的原矩阵比对)
    # 此处用转置和逆的结果一样redEigVects.I
    reconMat = (lowDDataMat * redEigVects.T) + meanVals
    print(reconMat)
    # 返回压缩后的数据矩阵及该矩阵反构出原始数据矩阵
    return lowDDataMat, reconMat


if __name__ == '__main__':
    dataMat = loadDataSet('数据.txt')
    lowDmat, reconMat = pca(dataMat, 1)
    print(np.shape(lowDmat))
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(dataMat[:, 0].flatten().A[0], dataMat[:, 1].flatten().A[0], marker='^', s=90)
    ax.scatter(reconMat[:, 0].flatten().A[0], reconMat[:, 1].flatten().A[0], marker='o', s=90, c='red')
    plt.show()