PCA原理与实现

Posted 2020-11-26 pengweiblog

参考：
如何理解主元分析（PCA）？
PCA的数学原理
转自：python实现PCA

★ PCA思路：

我们的初始矩阵为X，它是m×n维的矩阵，其中：m是该数据集有m条记录，n是每条记录中有n个特征，X的基本格式为：

我们看出矩阵X每一行就是一条记录，每一列就是特征，我们想要对它降维，也就是想让它从m×n维降到m×r维(n>r)，那到底删除哪几列损失信息才最少？或者通过乘以一个满秩矩阵映射到其他坐标系让它特征变少？其实无论你是直接在原数据上删除某几列还是想通过坐标变换映射到其他基坐标系，都可以通过让X右乘一个满秩矩阵的形式来实现。

设变换矩阵为P，P为n×r的满秩矩阵，原始数据集X经过P变换，成为Y(m×r维)，即：Y=X?P，我们知道右乘一个满秩矩阵来进行坐标变换，其实就是数据向这些新的基向量上投影，而这个满秩矩阵P就是新的基向量组成的矩阵：

P = [ P1, P2,...Pr ]??? 其中Pi为一个列向量，也就是一个基向量

原始数据集X我们是以单位矩阵E(n×n维)的每一列为基向量的，而Y我们是以{ P1,P2,...Pr }为基的。其实任何一个m×r维的满秩矩阵，都可以使原始矩阵X降维成m×r维，但是我们对降维后的数据有两个要求。

(1.) 第一个要求：

我们要求这变换后的矩阵Y，让它损失的信息最少，因为Y比原始矩阵X维数低，肯定在变换的过程中损失了一些信息，现在我们要求损失的信息最少，怎么衡量损失信息最少呢？我们想到：让变换后的数据散度最大。

上图中是二维数据降维到一维：左图是选x轴为基向量降到一维，可以看到四个点投影到x轴，因为重合还剩两个点了，大量信息丢失了；左图还画出了以y轴为基向量，四个点往y轴上投影，因为重合所以还剩3个点了，也是丢失了一部分信息。而右图，我们往该基向量上投影，发现四个点的散度比左图的大(虽然不是最大)，点没有重合，保留了4个点，因此，我们希望降维后的数据能够散度最大，才能尽量的减少信息损失。

在数学中，我们描述分散程度用方差来表示，在变换后的数据集矩阵Y(m×r维)中，一个字段的方差可以看做是每个元素与字段均值的差的平方和的均值，把数据集Y按列表示：Y=[ b1, b2,...br]? 其中(b_i)(1≤i≤r)是一个列向量。

再求出Y的每一列的均值为：(u1, u2, ...un)，故根据公式，每一列的方差为：

我们看到上述公式不是很简洁，因为还有个均值ui这一项，这是因为我们没对原数据集X的每一列进行0均值化处理，假如对原数据X进行了0均值化处理了，那么经过P变换(Y=X?P)形成的Y矩阵的每一列也会0均值化(这个不信的话大家可以举个例子看看)，即Y的每一列的均值就为0了，上述公式也就没有ui这一项了，因此我们需要先对原数据集X进行0均值化，0均值化大家都熟悉，就是每个数据减去它所在列的均值：

而Y为：

X零均值化后，经过P变换(Y=X?P)后的Y也是零均值化了，现在Y矩阵的每一列的方差公式就可以去掉均值(因为均值为0)，Y每一列的方差为：

Y均值化后，上述公式可写为矩阵形式：

(2.)? 第二个要求：

我们要求变换后的Y矩阵每一列(一列就是一个字段)之间都是不相关的，在数学中，描述变量之间的相关性可用协方差来描述，只要让协方差为0，就说明这两个变量(字段)线性无关，其实这个要求也好理解，假如降维后的Y矩阵，有某两列线性相关，那么其中一列肯定可以用令一列表示出来，那么给人的感觉就是这两列表达的信息重复了，既然是重复的列完全可以删除，所以才要求任意两列都是线性无关的，即任意两列协方差为0，比如Y中的任意两个字段(列)bp、bq，这两个字段的协方差为：

? 注意：因为Y已经被我们零均值化了，因此协方差公式中的均值是0，我们就不必要写出了，所以上述公式才会变得如此简洁，这就是零均值化的好处。

上式可以写为矩阵形式：

如果p，q你看不习惯，你可以换成i，j，只是表示符号不同而已：

现在我们来综合一下对降维后的矩阵Y的两个要求：一是、让Y每一列的自身方差最大；二是、让Y任意两列的协方差为0。

我们正好想到，Y的协方差矩阵，不正好就是这两个要求吗?因为Y被我们零均值化，所以，它的协方差矩阵，我们可以写为矩阵形式：

我们的要求是让Y每一列自身方差最大，任意两列协方差为0，正好不就是让Y的协方差矩阵Cov(Y)的对角线元素最大，让非对角元素为0，就是让(Y^Tcdot Y)是对角阵，且让这个矩阵的对角元素最大。也就是我们期望Cov(Y)是这个形式：

我们回想一下，P是m×r的满秩矩阵，我们的目的是找到这么一个矩阵P，使得(Y^Tcdot Y)满足上述要求，而上述公式根据大学学的线性代数，这个问题非常类似实矩阵对角化问题。

因为(frac{1}{m-1}X^Tcdot X)是实对称矩阵，所以一定可以相似对角化，且特征值(lambda_i≥0,iin(0,r))，即存在满秩矩阵Q(n×n维)使得：

而这个正交矩阵Q是由特征向量组成的矩阵[ Q1, Q2,...Qn ]，且Qi是与的顺序对应的，我们假设对角阵中的从大到小排列起来了，即(lambda_1≥lambda_2≥...≥lambda_n)，那么对应的特征向量Qi也要根据(lambda_i)的顺序排列好，我们观察目标矩阵Cov(Y)和对角化矩阵(Lambda)，发现Cov(Y)矩阵是(Lambda)的一部分，我们需要Cov(Y)的对角元素最大，也就是要取(Lambda)中的(lambda_i)最大的几个，而P也是根据选取的(lambda_i)而从Q中选取。比如我们要降维到m×r维，只需取(Lambda)中前r个特征值，然后从Q中找到这r个特征值对应的特征向量[ Q1,Q2,...Qr ]，这些特征向量组成的矩阵就是我们要找的P矩阵。

算法流程：

① 对原数据集零均值化。代码是：meanRemoved = dataMat - mean(dataMat,axis=0)

② 求出均值化X的协方差矩阵：公式是：，代码是：covMat = cov(meanRemoved,rowvar=0)

③ 求这个协方差矩阵的特征值，特征向量，代码是：eigVals, eigVects = linalg.eig(mat(covMat))

④ 把这些特征值按从大到小排列，返回特征值的下标，代码是：eigValInd = argsort(-eigVals)

⑤ 选出前topNfeat个特征值，返回这些选中的特征值的下标，并根据下标从特征向量矩阵eigVects中取出这些选中的特征向量组成矩阵P，这就是我们要找的变换矩阵P，代码是：redEigVects = eigVects[:,eigValInd[:topNfeat] ]

⑥ 返回降维后的数据，公式是：Y=X?P，代码是：lowDDataMat = meanRemoved * redEigVects

⑦ 原数据映射到新的空间中。公式是：，代码是：reconMat = (lowDDataMat * redEigVects.T) + meanValues

import numpy as np
import matplotlib.pyplot as plt
 
def loadDataSet(fileName):
    dataSetList = []
    fr = open(fileName)
    for row in fr.readlines():
        cur_line = row.strip().split('	')
        proce_line = list(map(float,cur_line))
        dataSetList.append(proce_line)
    dataSetList = np.array(dataSetList)
    
    return dataSetList


def pca(dataMat, topNfeat = 999999):
    meanValues = np.mean(dataMat,axis=0) # 竖着求平均值，数据格式是m×n
    meanRemoved = dataMat - meanValues  # 0均值化  m×n维
    covMat = np.cov(meanRemoved,rowvar=0)  # 每一列作为一个独立变量求协方差  n×n维
    eigVals, eigVects = np.linalg.eig(np.mat(covMat)) # 求特征值和特征向量  eigVects是n×n维
    eigValInd = np.argsort(-eigVals)  # 特征值由大到小排序，eigValInd十个arrary数组 1×n维
    eigValInd = eigValInd[:topNfeat]  # 选取前topNfeat个特征值的序号  1×r维
    redEigVects = eigVects[:,eigValInd] # 把符合条件的几列特征筛选出来组成P  n×r维
    lowDDataMat = meanRemoved * redEigVects  # 矩阵点乘筛选的特征向量矩阵  m×r维 公式Y=X*P
    reconMat = (lowDDataMat * redEigVects.T) + meanValues  # 转换新空间的数据  m×n维

    return lowDDataMat, reconMat



def drawPoints(dataset1,dataset2):  # 画图，dataset1是没降维的数据，dataset2是数据映射到新空间的数据
    fig = plt.figure()
    ax1 = fig.add_subplot(211)
    ax2 = fig.add_subplot(212)
    ax1.scatter(dataset1[:,0],dataset1[:,1],marker='s',s=40,color='red')
    dataset2 = np.array(dataset2)
    ax2.scatter(dataset2[:,0],dataset2[:,1],s=40,color='blue')
    plt.savefig('dataset1.png')
    plt.show()


if __name__ == '__main__':
    data = loadDataSet('C:/Users/Peng Wei/Desktop/pca_dataset2.txt')
    proccess_data, reconMat = pca(data,1)
    drawPoints(data,reconMat)