机器学习手记[11]---PCA变换基础

本文链接：https://blog.youkuaiyun.com/a_31415926/article/details/41089283

1. 初衷

1 达到降维的目的

2 将原来的维度进行转换，转换到新的维度组上，维度之间相互正交，且方差大小依次递减。

2. 需求

数据样本数有40个，特征有3个，如何降维到2个维度，且两个维度相互正交线性无关.

目的：将一组特征向量，通过线性变换，转换为一组线性无关的变量。

3. 步骤

1. 生成【3*40】的实验数据

维度为3，样本数为40，我们生成两个类的数据，每个类的样本是20个，但在PCA变换把整个3*40均当作一个类进行处理。

def loadData():
    mu_vec1 = np.array([0,0,0])
    cov_mat1 = np.array([[1,0,0],[0,1,0],[0,0,1]]) 
    class1_sample = np.random.multivariate_normal(mu_vec1, cov_mat1, 20).T
    assert class1_sample.shape == (3,20), "The matrix has not the dimensions 3x20"
    
    mu_vec2 = np.array([1,1,1])
    cov_mat2 = np.array([[1,0,0],[0,1,0],[0,0,1]])
    class2_sample = np.random.multivariate_normal(mu_vec2, cov_mat2, 20).T
    assert class2_sample.shape == (3,20), "The matrix has not the dimensions 3x20"
            
    all_samples=np.concatenate((class1_sample,class2_sample),axis=1)
    
    assert  all_samples.shape==(3,60),"The matrix has not the dimensions 3x40"
    return all_samples

2. 求特征协方差矩阵【3*40】

计算样本集在三个维度的平均值。
计算样本集在三个维度上的协方差之和

def scatterMatCreate(dataArray):
    all_samples=dataArray
    meanFirst=np.mean(all_samples[0])
    meanSecond=np.mean(all_samples[1])
    meanThird=np.mean(all_samples[2])
    
    meanVector=np.array([meanFirst,meanSecond,meanThird])
    
    scatterMat=np.zeros((3,3))
    for i in range(all_samples.shape[1]):
        matPartA=(all_samples[:,i]-meanVector).reshape(3,1)
        matPartB=(all_samples[:,i]-meanVector).reshape(1,3)
        scatterMat+=matPartA * matPartB
        scatterMat/float(all_samples.shape[1])    #将各个样本的方差进行累加
    return scatterMat

3 根据得到的【3*3】的协方差矩阵，求解他的特征向量和特征值。

特征向量组中包含的向量有三个，其中对应特征值最大的特征向量是主特征向量，这个向量对应的方向反应了整个数据集具有最大方差的方向，特征值次大的特征向量对应的则反应了次大的方差方向，依次类推。

在特征向量组【3*3】里面，每个特征向量的维度是【3*1】，因此我们筛选的两个特征向量组的大小是【3*2】

这三个特征向量是相互正交的。

def eigenCompute(scatterMat):        
    eigenValArray,eigenVecArray=np.linalg.eig(scatterMat)
    sortedIndex=np.argsort(eigenValArray)
    
    topVecNum=2
    eigenFinalList=[]
    for i in range(topVecNum):
        eigenVec=eigenVecArray[:,sortedIndex[i]] #注意计算好eigen的大小shape   eigenVec 3*1 
        eigenVec=list(eigenVec)
        eigenFinalList.append(eigenVec)
    eigenMat=np.array(eigenFinalList).reshape(2,3)
    return eigenMat

4 根据特征值大小，选取前两个特征值对应的特征向量，进行PCA特征变换

利用特征变换向量组进行特征变换。 EigenMat【2*3】*Data【3*40】=newData【2*40】

<span style="font-family: Arial, Helvetica, sans-serif;"></span><pre name="code" class="python">def PCArun():
    dataArray=loadData() #源数据
    scatterMat=scatterMatCreate(dataArray) #协方差矩阵
    eigenMat=eigenCompute(scatterMat)  #计算特征值和特征向量
    newMat=eigenMat.dot(dataArray)  #进行PCA特征变换
    
    return newMat

特征变换后的效果

5 恢复源维度信息

PCA实现了三维到二维的降维，而如果我们想从二维返回到三维信息呢？

5.1 如果我们不降维的话

我们直接将特征向量组进行特征变换，那么我们得到的PCAData数据也是三维的，三个维度相互正交。

DataMat*EigenMat=PCAData

如果要返回到来的数据，我们求逆。

DataMat=PCAData*inverse(EigenMat)

但是由于特征矩阵比较特殊是正交矩阵, 它的逆矩阵，等于，它的转置矩阵。

因此DataMat=PCAData*transpose(EigenMat)

这样我们求的DataMat和原始数据是完全一样的，因为我们没有降维，没有信息损失。也就是说我们还是用三个维度去描述我们的数据，只不过这三个轴的位置不一样了，但是数据本身并没有变化。

5.2 如果我们降维的话

降维的话，其实特征矩阵也满足这个条件，也就是说，我们还可以利用矩阵的逆等于矩阵的转置这个条件，进行计算。

PCA变换：

DataMat*EigenMat=PCAData 这个时候PCAData是二维的，DataMat是三维的。

PCA逆变换

DataMat=PCAData*transpose(EigenMat)

这个时候DataMat相对于PCAData肯定是有信息损失的。

可视化如下：

从视角A我们看到是散乱均匀分布的，但是当我们旋转的时候，我们就会发现，PCA变换恢复得到矩阵其实是少了一个维度的信息的。

整体的源代码如下，

# -*- coding: utf-8 -*-
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d import proj3d
from matplotlib.patches import FancyArrowPatch

def loadData():
    mu_vec1 = np.array([0,0,0])
    cov_mat1 = np.array([[1,0,0],[0,1,0],[0,0,1]])
    class1_sample = np.random.multivariate_normal(mu_vec1, cov_mat1, 20).T
    assert class1_sample.shape == (3,20), "The matrix has not the dimensions 3x20"
    
    mu_vec2 = np.array([1,1,1])
    cov_mat2 = np.array([[1,0,0],[0,1,0],[0,0,1]])
    class2_sample = np.random.multivariate_normal(mu_vec2, cov_mat2, 20).T
    assert class2_sample.shape == (3,20), "The matrix has not the dimensions 3x20"
            
    all_samples=np.concatenate((class1_sample,class2_sample),axis=1)
    
    assert  all_samples.shape==(3,40),"The matrix has not the dimensions 3x40"
    return all_samples
    
#协方差矩阵
def scatterMatCreate(dataArray):    
    all_samples=dataArray
    meanFirst=np.mean(all_samples[0])
    meanSecond=np.mean(all_samples[1])
    meanThird=np.mean(all_samples[2])
    
    meanVector=np.array([meanFirst,meanSecond,meanThird])
    
    scatterMat=np.zeros((3,3))
    for i in range(all_samples.shape[1]):
        matPartA=(all_samples[:,i]-meanVector).reshape(3,1)
        matPartB=(all_samples[:,i]-meanVector).reshape(1,3)
        scatterMat+=matPartA * matPartB
        scatterMat/float(all_samples.shape[1])    #将各个样本的方差进行累加
    return scatterMat    
    
#特征向量提取        
def eigenCompute(scatterMat):        
    eigenValArray,eigenVecArray=np.linalg.eig(scatterMat)
    sortedIndex=np.argsort(-eigenValArray) #降序排列
    
    topVecNum=2
    eigenFinalList=[]
    for i in range(topVecNum):
        eigenVec=eigenVecArray[:,sortedIndex[i]]  
        eigenVec=list(eigenVec)
        eigenFinalList.append(eigenVec)
    eigenMat=np.array(eigenFinalList).reshape(2,3)
    return eigenMat

#显示三维数据
def imgShow3DData(dataArray):
    dim=dataArray.shape[0]
    fig=plt.figure(figsize=(10,10))  #plt-->figure-->subplot-->plot
    ax=fig.add_subplot(111,projection='3d')
    ax.plot(dataArray[0,:],dataArray[1,:],dataArray[2,:],'o',markersize=8,color='green',alpha=0.5, label='Original Points')
    
def imgShowPCAData(dataArray):
    dim=dataArray.shape[0]
    fig=plt.figure(figsize=(10,10))  #plt-->figure-->subplot-->plot
    ax=fig.add_subplot(111)
    ax.plot(dataArray[0,:],dataArray[1,:],'o',markersize=8,color='green',alpha=0.5, label='PCAed Points')


#三维矢量箭头的显示方法
class Arrow3D(FancyArrowPatch):
    def __init__(self,xs,ys,zs,*args,**kwargs):
        FancyArrowPatch.__init__(self,(0,0),(0,0),*args,**kwargs)
        self._verts3d=xs,ys,zs
    
    def draw(self,renderer):
        xs3d,ys3d,zs3d=self._verts3d
        xs,ys,zs=proj3d.proj_transform(xs3d,ys3d,zs3d,renderer.M)
        self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
        FancyArrowPatch.draw(self,renderer)

#显示原始数据和两个特征向量方向  
def imgShowEigenVec(dataArray,eigenMat):
    fig=plt.figure(figsize=(10,10))  #plt-->figure-->subplot-->plot
    ax=fig.add_subplot(111,projection='3d')
    ax.plot(dataArray[0,:],dataArray[1,:],dataArray[2,:],'o',markersize=8,color='green',alpha=0.5, label='Original Points')
    
    meanFirst=np.mean(dataArray[0])
    meanSecond=np.mean(dataArray[1])
    meanThird=np.mean(dataArray[2])
    meanVector=np.array([meanFirst,meanSecond,meanThird])
    for v in eigenMat:
        a=Arrow3D([meanVector[0],v[0]],[meanVector[1],v[1]],[meanVector[2],v[2]],
                mutation_scale=20, lw=3,arrowstyle="-|>",color="r")
        ax.add_artist(a)

#PCA变换半恢复                
def restoredDataArray(PCADataArray, eigenMat):
    restoredDataArray=eigenMat.T.dot(PCADataArray)
    imgShow3DData(restoredDataArray)


#主函数 
def PCArun():
    dataArray=loadData()
    scatterMat=scatterMatCreate(dataArray)      #协方差矩阵
    eigenMat=eigenCompute(scatterMat)           #特征向量提取
    PCADataArray=eigenMat.dot(dataArray)        #PCA变换
      
    imgShow3DData(dataArray)                    #显示原始三维数据
    imgShowEigenVec(dataArray,eigenMat)         #显示原始数据和两个特征向量方向
    imgShowPCAData(PCADataArray)                #显示PCA变换后的2维数据
    restoredDataArray(PCADataArray, eigenMat)   #显示PCA变换半恢复后的数据

改编自:

可视化和数据源参考 http://sebastianraschka.com/Articles/2014_pca_step_by_step.html#introduction

PCA变换后的半恢复参考
http://www.cs.otago.ac.nz/cosc453/student_tutorials/principal_components.pdf