PCA算法的数学原理以及Python实现

ＰＣＡ算法的数学原理，请查看这篇博客

博客中的笔记：

降维当然意味着信息的丢失，不过鉴于实际数据本身常常存在的相关性，我们可以想办法在降维的同时将信息的损失尽量降低。

根据相关性来讲信息的损失量降到最低

更正式的说，向量(x,y)实际上表示线性组合：

x(1,0)?+y(0,1)?

不难证明所有二维向量都可以表示为这样的线性组合。此处(1,0)和(0,1)叫做二维空间中的一组基。

一般的，如果我们有M个N维向量，想将其变换为由R个N维向量表示的新空间中，那么首先将R个基按行组成矩阵A，然后将向量按列组成矩阵B，那么两矩阵的乘积AB就是变换结果，其中AB的第m列为A中第m列变换后的结果。

协方差（Covariance）在概率论和统计学中用于衡量两个变量的总体误差。而方差是协方差的一种特殊情况，即当两个变量是相同的情况。

那么如何选择这个方向（或者说基）才能尽量保留最多的原始信息呢？一种直观的看法是：希望投影后的投影值尽可能分散。

方差

上面的问题（尽可能投影的分散）被形式化表述为：寻找一个一维基，使得所有数据变换为这个基上的坐标表示后，方差值最大。

协方差

从直观上说，让两个字段尽可能表示更多的原始信息，我们是不希望它们之间存在（线性）相关性的，因为相关性意味着两个字段不是完全独立，必然存在重复表示的信息。

数学上可以用两个字段的协方差表示其相关性，由于已经让每个字段均值为0，则：

 

Cov(a,b)=1m∑i=1maibi

当协方差为0时，表示两个字段完全独立。为了让协方差为0，我们选择第二个基时只能在与第一个基正交的方向上选择。因此最终选择的两个方向一定是正交的

至此，我们得到了降维问题的优化目标：将一组N维向量降为K维（K大于0，小于N），其目标是选择K个单位（模为1）正交基，使得原始数据变换到这组基上后，各字段两两间协方差为0，而字段的方差则尽可能大（在正交的约束下，取最大的K个方差

）。

协方差矩阵求法：

接下来讲解一下我对于PCA算法的理解：

 

协方差矩阵对角化

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

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

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

2）设特征向量λ

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

 

PCA算法

总结一下PCA的算法步骤：

设有m条n维数据。

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

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

3）求出协方差矩阵C=1mXX?

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

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

6）Y=PX

即为降维到k维后的数据

PCA本质上是将方差最大的方向作为主要特征，并且在各个正交方向上将数据“离相关”，也就是让它们在不同正交方向上没有相关性。

Python代码：

转载自https://www.cnblogs.com/wzjhoutai/p/6819611.html

'''
@author: Garvin
'''
from numpy import *
import matplotlib.pyplot as plt

def loadDataSet(fileName, delim='\t'):
    fr = open(fileName)
    dataArr=[]
    for line in fr.readlines():
        stringArr= line.strip().split(delim)
        colArr=[]
        for numString in stringArr:
            colArr.append(float(numString))
        dataArr.append(colArr)
    return mat(dataArr)

def pca(dataMat, topNfeat=9999999):
    meanVals = mean(dataMat, axis=0)#获取平均值。
    meanRemoved = dataMat - meanVals #remove mean
    covMat = cov(meanRemoved, rowvar=0)#协方差矩阵
    eigVals,eigVects = linalg.eig(mat(covMat))#求特征值，特征向量
    eigValInd = argsort(eigVals)     #特征值排序       #sort, sort goes smallest to largest
    print(eigValInd)
    eigValInd = eigValInd[:-(topNfeat+1):-1]  #cut off unwanted dimensions减掉不想要的特征值
    print(eigValInd)
    redEigVects = eigVects[:,eigValInd]    #识别最大到最小的特征向量   #reorganize eig vects largest to smallest
    lowDDataMat = meanRemoved * redEigVects #将数据转到新的维度中#transform data into new dimensions
    reconMat = (lowDDataMat * redEigVects.T) + meanVals
    return lowDDataMat, reconMat #转化回来的矩阵，以及最终的降维矩阵

##画出最合适的方向。
def plotBestFit(dataSet1,dataSet2):      
    dataArr1 = array(dataSet1)
    dataArr2 = array(dataSet2)
    n = shape(dataArr1)[0] 
    n1=shape(dataArr2)[0]
    xcord1 = []; ycord1 = []
    xcord2 = []; ycord2 = []
    xcord3=[];ycord3=[]
    j=0
    for i in range(n):
            xcord1.append(dataArr1[i,0]); ycord1.append(dataArr1[i,1])
            xcord2.append(dataArr2[i,0]); ycord2.append(dataArr2[i,1])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
    ax.scatter(xcord2, ycord2, s=30, c='green')
    plt.xlabel('X1'); plt.ylabel('X2');
    plt.show()    

if __name__=='__main__':
     mata=loadDataSet('testSet.txt')
     #print(mata)
     a,b= pca(mata, 2)
     plotBestFit(a,b)