[Python知识库] Python学习（3）---主成分分析（PCA）的基本原理及其Python实现

开发: C++知识库 Java知识库 JavaScript Python PHP知识库人工智能区块链大数据移动开发嵌入式开发工具数据结构与算法开发测试游戏开发网络协议系统运维
教程: HTML教程 CSS教程 JavaScript教程 Go语言教程 JQuery教程 VUE教程 VUE3教程 Bootstrap教程 SQL数据库教程 C语言教程 C++教程 Java教程 Python教程 Python3教程 C#教程
数码: 电脑笔记本显卡显示器固态硬盘硬盘耳机手机 iphone vivo oppo 小米华为单反装机图拉丁

-> Python知识库 -> Python学习（3）---主成分分析（PCA）的基本原理及其Python实现 -> 正文阅读

[Python知识库]Python学习（3）---主成分分析（PCA）的基本原理及其Python实现

? ? ?一、降维的基本概念

????????对于实际分析过程中的高维数据，在进行具体的数据分析和特征建模之前，需要进行数据降维处理。降维是指通过某种方法从原始数据的N个特征中选取K个（K<N）进行数据表示，在减少数据信息丢失的前提下实现原始数据的压缩表示，其主要目的包括以下几点：

? ? ? ? （1）提高数据的使用率；

? ? ? ? （2）减低算法的计算开销（软件和硬件）；

? ? ? ? （3）去除干扰噪声；

? ? ? ? （4）提高分析结果的解释性。

????????? 现有的数据降维方法很多，按照个人的了解（存在认识不足），根据输入变量X和目标变量的关联关系，可分为：输入变量数据降维（主成分分析PCA、独立主成分分析ICA、因子分析FA等）和关联数据降维（偏最小二乘PLS、Lasso、逐步线性回归、回归数（回归数也可以做数据降维）等），在上述方法中应用最广、传播度最大的属于PCA，本次仅对其基本原理和分析思路进行介绍，最后通过python编程，感受一下其巨大的魅力，最后分析过程的局限性。

? ? ? ? 二、PCA

? ? ? ? 在周志华老师的《机器学习》书中写到：方差越大，数据所含表示信息量越多，意味着数据集不同样本、不同维度之间的变化能够完全表征样本特征的可能性越高，反之亦然。这句话是什么意思呢？就是某一维变量的方差越大时，其所能包含的信息量越大，表明其在整个数据中的‘话语权’越大，通俗的说就是见多识广，知识范围越广，越能恰当、清楚地表述一件事物各方面的特征。PCA的主要思想就是基于上述原理，即通过选取K个能够最大量度覆盖原始数据的特征（N个），实现数据降维，进而提高模型解释性及进行可视化分析等工作。

? ? ? ? 在分析PCA基本原理之前，需要理解以下几点概念：

? ? ? ? 1. 向量内积的几何意义

? ? ? ? 假设A $(x_1, y_1)$ 、B $(x_2, y_2)$ 为二维平面上的两个非零向量，其内积可表示为：

$AB = \left | A \right |\left | B \right |cos\theta =x_1 * x_2 + y_1 * y_2$

坐标轴表示为：

?其中， $\left | A \right |cos\theta$ 实际上表示向量A在向量B方向上的投影。当B为单位向量时，A和B的内积表示A向量在B向量上的投影，且满足上式。

? ? ? ? 2. 基的变换

? ? ? ? 基变换就是坐标轴变换（PCA的主要操作），通过新的正交轴去代替原始坐标轴（XOY坐标轴）。此处新的坐标轴的基向量的模必须为1且两两正交，主要原因是由两两正交的单位向量组合的正交矩阵在将原始坐标轴正交变换为新坐标轴时，不会改变向量基本特征属性（模、夹角、距离）。此外，正交矩阵的基本特性（ $A^{'}A = I$ ）便于矩阵归约处理。

? ? ? ? 对于向量 $A (x_a, y_a)$ ，假设新的坐标轴（注：新的坐标轴的基向量表示仍然是基于旧坐标轴而言的）可分别表示为：x轴 $\vec{i}: (x_i, y_i)$ 和y轴 $\vec{j}: (x_j, y_j)$ 。那么，按照第1节表述的内积的基本应用，向量A在新坐标轴上可表示为：

$(x_a, y_a)\begin{pmatrix} x_i &x_j \\ y_i&y_j \end{pmatrix}=(x_a * x_i +y_a * y_i, x_a * x_j +y_a * y_j)=(x_{new}, y_{new})$

其中， $\begin{pmatrix} x_i &x_j \\ y_i&y_j \end{pmatrix}$ 为以列向量排列的矩阵，即为变换矩阵； $(x_{new}, y_{new})$ 为在新坐标轴下的坐标表示。式中右乘的原因是因为变换矩阵由列向量组成，每一个列表示一个新的坐标轴，按照内积的计算方法，因此选择右乘。

? ? ? ? 举个例子，对于XOY坐标系中的向量（3， 2）, 假设选定新的坐标轴其基向量分别为 $(1/\sqrt{2}, 1/\sqrt{2})$ 和 $(-1/\sqrt{2}, 1/\sqrt{2})$ ，按照上节表述的内积的意义，可得向量在新的的坐标轴的表示为：

? $(x_{new}, y_{new}) =(3, 2 )\begin{pmatrix} 1/\sqrt{2} &-1/\sqrt{2} \\ 1/\sqrt{2}& 1/\sqrt{2} \end{pmatrix}=(5/\sqrt{2}, -1/\sqrt{2})$

具体表示为：

? ? ? ? 3. 协方差矩阵

? ? ? ? 前面讲到：方差越大，信息量越大，对于下图所示数据，其XOY坐标表示的方向均不是方差最大的方向。结合坐标变换，我们可以重新构建一个基于上述理念的坐标系（红线），让样本信息尽量集中在少数几个坐标轴中，剩余坐标轴几乎不含信息，即可去掉此坐标轴实现特征降维。

? ? ? ? 现在的关键问题是如何衡量确定坐标轴所能表示的信息量。首先考虑的是方差，对于转换后的坐标，假设数据在某一坐标轴上的方差越大，则可认为数据在该坐标轴上所含的信息量越大，这样我们就可以根据方差的大小依次选取所需要的的特征（或坐标轴）进行降维。但是（最怕这个），我们在做数据分析的时候不得不考虑一个问题，就是数据的相关性问题。

? ? ? ? 假设选择方差最大的轴作为第一维（第一变换方向），当变量之间关联性较低或者不相关时，可考虑依次考虑方差递减对应轴作为候选变换方向；但是当候选第一维和第二维之间存在较强关联性时，即两者表示的信息有很大部分相同的（此处的相同不是说值的大小，而是所表征的样本特征变化情况，例如具有同时变化趋势），此时选择第二维并不能最大化信息表示，反而会降低样本信息的有效表示。因此，需要考虑其他信息考量指标。

? ? ? ? 协方差矩阵表征了两两变量之间是否同时偏离均值，其中对角线变量表示每个特征的方差，非对角线变量表示不同特征之间的相关性，绝对值越大，表明变量之间相关性越强。协方差矩阵结合了方差越大，信息量越大的特点，同时表征了不同维度之间的相关性，因此可用于变量信息量度量。

? ? ? ? 假设经过归一化处理后的输入数据为 $X\in R^{m\times n}$ （归一化用于降低不同维度取值范围差异的干扰），则其协方差矩阵为：

$C_X =\frac{1}{m}X^TX =\begin{bmatrix} Cov(x_1, x_1) & Cov(x_1, x_2) & \cdots &Cov(x_1, x_m) \\ Cov(x_2, x_1) & Cov(x_2, x_2) & \cdots & Cov(x_2, x_m)\\ \cdots & \cdots & \cdots& \cdots\\ Cov(x_m, x_1)& Cov(x_m, x_2) & \cdots & Cov(x_m, x_m) \end{bmatrix}$

?假设有变换矩阵 $P\in R^{n\times k}(k<n)$ ，根据矩阵投影，所得 $Y=XP$ 可实现n维到k维的降维处理。其中Y矩阵的协方差为：

$C_Y=\frac{1}{m}(XP)^T(XP)=\frac{1}{m}P^TX^TXP=P^T(\frac{1}{m}X^TX)P=P^TC_XP$

?根据矩阵迹以及协方差矩阵的定义，最大化保留数据信息等价于最大化协方差矩阵的迹，即优化目标可表示为：

$J=\left\{\begin{matrix} max(tr(P^TC_XP))\\ P^TP=I\end{matrix}\right.$

式中，? $tr(C_Y)=tr(P^TC_XP)=Cov(y_1,y_1)+Cov(y_2,y_2)+\cdots +Cov(y_m,y_m)$ 表示变换后坐标协方差矩阵的迹， $P^TP=I$ 是变换矩阵为正交矩阵的约束条件。根据拉格朗日方法，有：

$J(P)=tr(P^TC_XP)+\lambda (I-P^TP)$

因为前面我们采用列向量表示样本，此处求导过程采用分母布局，有：

$\frac{\partial J}{\partial P}=\frac{\partial tr(P^TC_XP)}{\partial P}+\frac{\partial (-\lambda P^TP)}{\partial P}=\frac{\partial tr(PP^TC_X)}{\partial P}- \lambda P$

?上式用到以下性质： $\frac{\partial X^TX}{\partial X}=X$ ， $tr(AB)=tr(BA)$ 。

进一步，有：

$\frac{\partial J}{\partial P}=\frac{\partial tr(PP^TC_X)}{\partial P}- \lambda P=(P^TC_X)^T- \lambda P=C_X^TP- \lambda P=C_XP- \lambda P$

?上式用到的性质有 $\frac{\partial tr(AB)}{\partial A}=B^T$ ，令偏导表达式为0，则有：

$\frac{\partial J}{\partial P}\Rightarrow C_XP-\lambda P=0 \Rightarrow C_XP=\lambda P$

根据特征值分析定义可知，最终随求结果满足特征向量和特征值的关系表达式，即由 $C_X$ 矩阵的特征向量组成的基矩阵，就是所求的变换矩阵P。根据特征变换的基本意义（特征值 $\lambda$ 表示矩阵 $P$ 能够将矩阵 $C_X$ 进行保留特征变换的倍数值），数据降维过程中信息量保存能力最大的基向量一定是输入矩阵X协方差矩阵的最大特征值对应的特征向量，其保持信息量的大小（大小不太准确，应该是比例）就是其对应特征值的绝对值。因此，我们可以对输入矩阵X的协方差矩阵进行特征值求解，根据特征值绝对值大小依次选择对应变换向量，进而实现数据降维处理。

????????至此，我们介绍了PCA的基本处理过程和分析原理，简单表述如下：

? ? ? ? step1：输入矩阵归一化处理；

? ? ? ? step2：计算样本协方差矩阵；

? ? ? ? step3：求解协方差矩阵指定数目最大特征值对应特征向量；

? ? ? ? step4：确定变换矩阵，求解降维数据。

? ? ? ? 三、Python编程实现PCA分析

? ? ? ? 基于上述步骤，我们依次进行子函数编程（个人喜好matlab，很多表述用于python不太准确，敬请谅解），有：

#/usr/nom/env python
# _*_coding:utf-8_*_
# @Time      :2021/9/3 10:04
# @Author    :A bigfish
# @FileName  :maindemo13.py
# @Software  :PyCharm

import matplotlib.pyplot as plt
import numpy as np
from pylab import *

# 首先导入数据，此部分为从存储列表或单元中读取分析数据
def loadDataSet(filename, delim='\t'):    #此处的'\t'表示不同变量间的分隔符，t表示tab键键入的空格
    fr = open(filename)
    stringArr = [line.strip().split(delim) for line in fr.readlines()]
    dataArr = [list(map(float, line)) for line in stringArr]
    return np.mat(dataArr)


# 定义pca分析函数
def pca(dataset, topNfeat = 99999):        #topNfeat最大特征值数目，通常不用设置，因为后续要进行可视化分析                
    meanVals = np.mean(dataset, axis=0)    #求均值
    meanRemoved = dataset - meanVals       #预处理
    covMat = np.cov(meanRemoved, rowvar=0) #求解输入数据协方差矩阵
    eigVals, eigVects = np.linalg.eig(np.mat(covMat))    #求解特征值，特征向量
    eigVaInd = np.argsort(eigVals)         #对特征值进行排序处理，默认为升序
    eigVaInd = eigVaInd[-1:-(topNfeat):-1] #根据指定数目进行逆序处理
    redEigVects = eigVects[:,eigVaInd]     #选取对应特征向量
    lowDataMat = meanRemoved * redEigVects #数据降维X*P
    recontMat = (lowDataMat * redEigVects.T) + meanVals #c处理进行了数据重构，非必须选项
    return lowDataMat, recontMat, eigVals  #返回数据

# 定义特值值绘制函数
def plotEig(dataset, numFeat=20):            
    mpl.rcParams['font.sans-serif'] = ['Times NewRoman']
    sumData = np.zeros((1, numFeat))
    dataset = dataset / sum(dataset)
    for i in range(numFeat):
        sumData[0, i] = sum(dataset[0:i])

    X = np.linspace(1, numFeat, numFeat)
    fig = plt.figure()
    ax = fig.add_subplot(211)
    ax.plot(X, (sumData*100).T, 'r-+')
    mpl.rcParams['font.sans-serif'] = ['SimHei']
    plt.ylabel('累计方差百分比')

    ax2 = fig.add_subplot(212)
    ax2.plot(X.T, (dataset[0:numFeat].T)*100, 'b-*')
    plt.xlabel('主成分数')
    plt.ylabel('方差百分比')
    plt.show()

# 定义原始数据及第一主成分绘制函数
def plotData(OrigData, recData):
    import matplotlib.pyplot as plt
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(OrigData[:,0].flatten().A[0], OrigData[:, 1].flatten().A[0], c='blue',marker='^', s=90)
    ax.scatter(recData[:, 0].flatten().A[0], recData[:, 1].flatten().A[0], c='red', marker='o',s=90)
    plt.show()

在main窗扣，有：


if __name__ == '__main__':
    import maindemo13
    dataMat = maindemo13.loadDataSet('secom.txt', ' ') 
    newData = maindemo13.replaceNan(dataMat)    #此部分为缺失值Nan处理函数，未在上面进行分析
    lowData, recData, eigVals = maindemo13.pca(dataMat, 50)    #选取前50个特征
    maindemo13.plotData(dataMat, recData)    #绘制原始数据及第一主成分方向
    maindemo13.plotEig(eigVals, 50)    #绘制特征值累计贡献率

所得原始数据（蓝色）和第一主成分分析方向（红色）如下：