PCA主成分分析

1、概述

    PCA（principle component analysis），也就是主成分分析，其主要目的或者作用是，将原来的数据变换到一个新的空间中，得到数据更好的表达方式，去除影响较小的部分，简化数据的维度。

    在看《机器学习实战》时候，一次看下，完全不知道他是怎么做的，网上搜索无数帖子，基本都在讨论或者阐述PCA实现的流程，先怎么样然后怎么样，很少看到讨论我们为什么要这么做，这么做的原因到底是什么。为什么要中心化？为什么要求相关系数矩阵？为什么相关系数矩阵的特征值大的包含的信息量越大（主成分），然后原始数据乘以该特征值对应的特征向量就可以得到原始数据到新的空间的变换？如此等等，都很少有人去解释，解释了也没有解释清楚，或者我等学渣看的云里雾里。

    最后发现一篇文章，可谓白话PCA：PCA的数学原理。需要纠正的是，文章中使用了协方差矩阵，而不是相关系数矩阵，应该使用相关系数矩阵，至于为什么，就要讨论标准化的作用了，参考 通俗理解协方差和相关系数

    这里对PCA数学原理一文和近来所看的一些博文进行整理，透析PCA的本质。

2、需要使用的数学工具

    （1）协方差和相关系数矩阵

    （2）矩阵的行列式计算

    （3）矩阵的特征值和特征向量

    （4）相似矩阵、矩阵的对角化

    逐一翻阅《线性代数》都有，呵呵，都忘了。

3、PCA原理

3.1 向量和基变换

    特征一般用向量来表示，一行表示一个采样，每一列是一个维度，每个维度上的数其实就是向量在该维上的投影大小。一般我们使用笛卡尔直角坐标系，但如果将坐标系旋转我们可以得到很多坐标系。一组特征在笛卡尔坐标系中，可能很不好表达，直观上不好理解数据。

    比如上图，是二维笛卡尔直角坐标系中的一些点，他们在横轴和纵轴上的投影分布差不多，因此需要要两个维度来表示。但如果我们将坐标系逆时针旋转45°，使得这些点的“主方向”落在横轴上呢，那在纵轴上的数据就显得简单汗多了，可能只是一些噪声，或许我们可以只用横轴就可以保留数据尽可能大的信息量。

    向量和速度一样，是相对的，可以用不同的坐标系作为参考来观察他们。而PCA的任务就是找到合适的坐标系来表示数据，这个合适，就是要减小每一个维度的相关性，因为两个维度相关了，数据就会冗余，比如上图中，x增加的时候y也会相应增加。如果每个维度之间的相关度减小了，那我们就希望数据会落在尽量少的维度中，比如上图，在牺牲一点信息量的条件下，我们可以只用一个维度来表示数据。而新的坐标满足的条件就是：数据在该维上的方差尽量大，如上图，我们直观上就会将主方向如下选择：

因为这个方向上数据分散，波动大，影响也大，关于这点，可以参考：学习成绩的分析，其中阐述了某一个科目成绩对总学习成绩的影响。

    将数据变换到新的坐标系的过程叫做基变换，向量的基可以是相关的，但这与我们这里的目的相悖，所以我们要寻找的基应该是正交基，且对基单位化（长度是1）

因此PCA可以变为如下问题：寻找新的正交坐标系表达数据，使得数据在新的维度上方差尽可能大，且每个维度之间互不相关。

3.2协方差和相关系数

    上面给出了通俗理解协方差和相关系数矩阵的链接，这里简单总结一下。

    协方差是描述两组数据直接相关程度的参数，它不考虑数据的量纲，而相关系数则考虑数据的量纲，可以证明，相关系数是[-1,1]的数，而相关系数矩阵，则是以矩阵形式表示两组数据的相关性:

2020-2-14注：此处有坑，np的Cov有参数bias，默认是False，计算协方差时候除以(N-1)而不是N，因此正常要把bias设置为True才可以得到样本的方差和协方差

    相关系数大于0，表示正相关，趋势相同，否则表示负相关。numpy的cov返回的是一个对称矩阵，斜对角是每一组数据的方差的无偏估计，(i,j)表示数据i和数据j的协方差。

    设有一组数据X，有a和b两个维度，共m组数据：

    

    用X乘以X.T，可以得到a和b两个维度的协方差矩阵。

    设已经对X进行了标准化，则上面求出来的就是相关系数矩阵。一般地，设我们有m个n维数据记录，将其按列排成n乘m的矩阵X（默认已经标准化），那么它的相关系数矩阵为：

    所以我们先对数据标准化，然后用numpy的cov可以得到相关系数矩阵。

3.3 PCA的推导

    回顾上面PCA演变出来的问题：寻找新的正交坐标系表达数据，使得数据在新的维度上方差尽可能大，且每个维度之间互不相关。结合相关系数矩阵的知识，我们可以把问题重新定义如下：寻找新的正交坐标系表达数据，使得数据的相关系数矩阵是对角矩阵。

    为什么是对角矩阵呢？因为我们知道，A的相关系数矩阵的对角线上式每一维数据的方差的无偏估计，而其他位置的数据表示维度与维度之间的相关性，我们让这些位置为0，那么就保证了每一个维度之间不存在相关性，既然如此，那么维度上的方差就能达到最大。我们只要选择前N个大的方差（前提是前N个大的方差占据了主导优势，剩下的那些方差都相对很小）对应的维度来表示数据，就能够达到简化数据维度的目的。注意这里延伸出来的问题，我们最终是要寻找对应的维度，而不是方差，也就是寻找用来表达数据的新的坐标系（一组基）

    这里注意：我们使用的是相关系数矩阵，不是协方差矩阵，但是描述和推导的时候还是说协方差矩阵，因为我们已经对数据进行了标准化。下面参考《PCA的数学原理》对PCA进行推导。

 

    设原始数据矩阵X对应的协方差矩阵为C，而P是一组基按行组成的矩阵

    设Y=PX，则Y为X对P做基变换后的数据。设Y的协方差矩阵为D，我们推导一下D与C的关系：

    

现在事情很明白了！我们要找的P不是别的，而是能让原始协方差矩阵对角化的P。换句话说，优化目标变成了寻找一个矩阵P，满足 是一个对角矩阵，并且对角元素按从大到小依次排列，那么P的前K行就是要寻找的基，用P的前K行组成的矩阵乘以X就使得X从N维降到了K维并满足上述优化条件。

 

由上文知道，协方差矩阵C是一个是对称矩阵，在线性代数上，实对称矩阵有一系列非常好的性质：

1）实对称矩阵不同特征值对应的特征向量必然正交。

2）设特征向量λ重数为r，则必然存在r个线性无关的特征向量对应于λ，因此可以将这r个特征向量单位正交化。

由上面两条可知，一个n行n列的实对称矩阵一定可以找到n个单位正交特征向量，设这n个特征向量为e1,e2,⋯,en，我们将其按列组成矩阵：

 

则对协方差矩阵C有如下结论：

其中Λ为对角矩阵，其对角元素为各特征向量对应的特征值（可能有重复）。

以上结论不再给出严格的数学证明，对证明感兴趣的朋友可以参考线性代数书籍关于“实对称矩阵对角化”的内容。到这里，我们发现我们已经找到了需要的矩阵P:

 

 

 

P是协方差矩阵的特征向量单位化后按行排列出的矩阵，其中每一行都是C的一个特征向量。如果设P按照Λ中特征值的从大到小，将特征向量从上到下排列，则用P的前K行组成的矩阵乘以原始数据矩阵X，就得到了我们需要的降维后的数据矩阵Y。

3.4 PCA手算示例

    为了便于查阅，这里直接摘抄《PCA的数据原理》中的示例：

PCA算法   

 总结一下PCA的算法步骤：

设有m条n维数据。

1）将原始数据按列组成n行m列矩阵X

2）将X的每一行（代表一个属性字段）进行零均值化，即减去这一行的均值

3）求出协方差矩阵

4）求出协方差矩阵的特征值及对应的特征向量

5）将特征向量按对应特征值大小从上到下按行排列成矩阵，取前k行组成矩阵P

6）Y=PX即为降维到k维后的数据

4、PCA的Python实现

  实验数据来自《机器学习实战》第十三章，原始数据在3.1节已经给出。

 

def load_dataset():
    f = open("D:\\leran\PCA\\testSet.txt")
    d = []
    for l in f.readlines():
        ls = l.strip().split("\t")
        d.append([float(ls[0]), float(ls[1])])
    return np.matrix(d)

def pca1():
    data = load_dataset()
    mean_val = np.mean(data, axis=0)
    print("mean_val:", mean_val)
    mean_dmat = data - mean_val  # 数据中心化
    std_dev = np.sqrt(np.var(mean_dmat, axis=0)).A1
    print("std dev:", std_dev)
    mean_dmat /= std_dev  # 数据标准化
    # cov_m = np.cov(mean_dmat, rowvar=0)  # 协方差矩阵
    cov_m = mean_dmat.T * mean_dmat / mean_dmat.shape[0]  # 计算协方差矩阵（协方差系数）
    print("cov_mat:", cov_m)
    eig_value, eig_vector = np.linalg.eig(cov_m)  # 特征值和特征向量
    print("eig_value:", eig_value)
    print("eig_vector:", eig_vector)

    top_n_feat = 1  # 留下方差（特征值）最大的前几个特征向量

    eigval_ind = np.argsort(eig_value)  # 按照特征值大小排序，输出值在列表中的下标位置
    print("***", eigval_ind)
    eigval_ind = eigval_ind[-1:-(1 + top_n_feat):-1]  # 从最后一个向前取top_n_feat个，也就是最大的top_n_feat个
    print("pc eig:", eigval_ind)
    red_eig_vectors = eig_vector[:, eigval_ind]  # 去除对应的特征向量
    print("pc vector:", red_eig_vectors)
    lowd_date_mat = mean_dmat * red_eig_vectors  # 将原始数据变换到新的坐标系中，实际应用时候要这么用得到新的特征向量
    print(lowd_date_mat.shape)
    # 将新坐标系中的数据反变换会原来的坐标系，由于部分信息已经丢失，这里的recon_mat并不完整
    recon_mat = (lowd_date_mat * red_eig_vectors.T)

    # recon_mat[:, 0] *= std_dev[0]
    # recon_mat[:, 1] *= std_dev[1]
    # recon_mat += mean_val
    print(recon_mat.shape)

    k1 = eig_vector[1, 0] / eig_vector[0, 0]
    k2 = eig_vector[1, 1] / eig_vector[0, 1]
    # print(k1*k2)
    # k1 *= std_dev[1] / std_dev[0]
    # k2 *= std_dev[1] / std_dev[0]
    # print("vectors k:", k1, k2, k1 * k2)
    x1 = np.linspace(-5, 5, 10)
    # plt.figure("pca1")
    plt.plot(x1, x1 * k1, "r")  # + mean_val[0, 0] + mean_val[0, 1]
    plt.plot(x1, x1 * k2, "b")  # + mean_val[0, 0] + mean_val[0, 1]
    # plt.scatter(eig_vector[:, 0].flatten(), eig_vector[:, 1].flatten(), marker="o", color="r", alpha=1)
    plt.scatter(recon_mat[:, 0].A1, recon_mat[:, 1].A1, marker="*", color="r", alpha=1)
    plt.scatter(mean_dmat[:, 0].A1, mean_dmat[:, 1].A1, marker="*", color="k", alpha=0.5)
    plt.scatter(data[:, 0].A1, data[:, 1].A1, marker="*", color="k", alpha=0.5)
    plt.show()

代码各部分的注释已经都写的很清楚，这里不再赘述。

最后的运行结果：

mean_val: [[ 9.06393644  9.09600218]]
std dev: [ 1.0251496  1.486423 ]
cov_mat: [[ 1.          0.73730306]
 [ 0.73730306  1.        ]]
eig_value: [ 1.73730306  0.26269694]
eig_vector: [[ 0.70710678 -0.70710678]
 [ 0.70710678  0.70710678]]
*** [1 0]
pc eig: [0]
pc vector: [[ 0.70710678]
 [ 0.70710678]]
(1000, 1)
(1000, 2)

5、小结

    PCA的代码比较简单，几行代码就完成了，但其后面蕴含的知识比较宽广，所以网上很多人都只是说一下PCA的公式在代码中怎么用，说其然不说其所以然。