PCA的原理与实现

PCA原理

引入

PCA主要用于数据降维。

假设将二维数据降到一维，如下图所示，则希望变换后的散度越大越好，如C和D因为太靠近就会重合，导致丢失信息。那就是使方差最大化。我们还希望不同特征之间不相关（如果两个特征线性相关，说明其中一个可以用另一个来表达，那信息就重复了），而要线性无关，即协方差为0。

所以我们需要让每一个特征内方差最大，不同特征间协方差为0。

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

目的是让X从m×n降维到m×r。

将该矩阵变换到新的坐标空间就是右乘新坐标空间的基向量组成的矩阵：

求方差

每一列的方差为：

这是没有对X进行均值化的情况，如果对X均值化：

X均值化后，Y=X·P也均值化了：

现在；Y的每一列的方差为：

写成矩阵形式：

求协方差

Y中任意两列的协方差为：

协方差公式中还包含均值，这里Y已经0均值化，所以就没有了。

上式写成矩阵形式：

满足两个条件

现在，我们要让Y的每一列列内的方差最大，不同列间的协方差为0。

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

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

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

因为\(\frac{1}{m-1}X^T\cdot X\)是实对称矩阵，所以一定可以相似对角化，且特征值\(\lambda_i≥0,i\in(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矩阵。

说明下，虽然样本方差的分母是应该为m-1，但是上面所以式子分母都可以用m，分母采用m是因为这样算出来的样本方差Var(X)为一致估计量，不会太影响计算结果并且可以减小运算负担。

算法实现

算法流程

① 对原数据集零均值化。代码是：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('\t')
        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)

原文地址：https://www.cnblogs.com/pengweiblog/p/12572257.html