主成分分析（Principal Components Analysis, PCA）

有效解决问题的前提是选择一个看待问题的好视角，使问题展现其最清晰简单的一面。所谓主成分分析，就是选择一个可以将各待分析样本显著区分开的视角。

一、目标

设有$m$个$n$维向量 ${\mathbf{x}}_1，{\mathbf{x}}_2，\cdots {\mathbf{x}}_m$，代表$m$个样本，其中每个样本都由$n$个指标来描述，即
$$\mathbf{X}=\left[\begin{array}{cccc}{\mathbf{x}}_1& {\mathbf{x}}_2& \cdots & {\mathbf{x}}_m\end{array}\right]=\left[\begin{array}{cccc}{x}_{11}& x_{21}& \cdots & x_{m1}\\ x_{12}& x_{22}& \cdots & x_{m2}\\ ⋮& ⋮& ⋮& ⋮\\ x_{1n}& x_{2n}& \cdots & x_{mn}\end{array}\right]$$而一个指标存在的意义就在于有效的把各样本区别开，如果对于某项指标，所有样本的取值都几乎相同，即通过该指标无法把各样本区分开，那说明该指标对于分析这组样本没什么意义。因此用于标识样本的一组指标最好在任一指标维度上都能尽量把各样本区分开，这就需要找到一个观察各样本的理想视角，在这个视角下，各样本的关键差异被体现的淋漓尽致。如下图所示，三维空间中的每个点代表一个用三个指标描述的样本，但其实只需要一个新指标就可以把所有这些样本有效区分开。

用数学语言表述就是在$n$维空间中找到使各样本向量在其上投影方差尽可能大的一系列基向量。

二、途径

（一）将各指标中心化

就是对每个样本的每个指标维度都减去该指标的各样本平均值。从而使得所有样本向量的和是一个0向量，即
$$\sum _{i=1}^m{\mathbf{x}}_i=0$$目的还是为了只抓关键。

（二）投影到标准正交基

$n$维空间内的任一个$k$维子空间（$k<n$）的一组标准正交基可表示为：
$$\mathbf{e}=\left[\begin{array}{c}{\mathbf{e}}_1\\ {\mathbf{e}}_2\\ ⋮\\ {\mathbf{e}}_k\end{array}\right]=\left[\begin{array}{cccc}{e}_{11}& e_{12}& \cdots & e_{1n}\\ e_{21}& e_{22}& \cdots & e_{2n}\\ ⋮& ⋮& ⋮& ⋮\\ e_{k1}& e_{k2}& \cdots & e_{kn}\end{array}\right]$$其中每个基向量都是一个$n$为向量，但向量组的秩是$k$。任一样本向量投影到这组正交基各维度的坐标向量为：
$${\mathbf{y}}_i=\mathbf{e}{\mathbf{x}}_i$$即将原来用$n$维向量$\mathbf{x}$表示的各样本改用新的$k$维向量$\mathbf{y}$表示。从而实现换一个视角来看待各样本。

（三）计算各基向量上的投影方差

对于上第$j$个新基向量，各样本在其上投影的方差为：
$${\lambda }_j=\frac{1}{m}\sum _{i=1}^m({\mathbf{e}}_j{\mathbf{x}}_i-{\mathbf{e}}_j\overline{\mathbf{x}}{)}^2=\frac{1}{m}\sum _{i=1}^m({\mathbf{e}}_j{\mathbf{x}}_i{)}^2=\frac{1}{m}\sum _{i=1}^m{\mathbf{e}}_j{\mathbf{x}}_i{\mathbf{x}}_i^T{\mathbf{e}}_j^T={\mathbf{e}}_j\left(\frac{1}{m}\sum _{i=1}^m{\mathbf{x}}_i{\mathbf{x}}_i^T\right){\mathbf{e}}_j^T$$其中，
$$\begin{array}{rl}\frac{1}{m}\sum _{i=1}^m{\mathbf{x}}_i{\mathbf{x}}_i^T& =\frac{1}{m}\sum _{i=1}^m\left[\begin{array}{cccc}{x}_{i1}{x}_{i1}& x_{i1}{x}_{i2}& \cdots & x_{i1}{x}_{in}\\ x_{i2}{x}_{i1}& x_{i2}{x}_{i2}& \cdots & x_{i2}{x}_{in}\\ ⋮& ⋮& \ddots & ⋮\\ x_{in}{x}_{i1}& x_{in}{x}_{i2}& \cdots & x_{in}{x}_{in}\end{array}\right]\\ & =\frac{1}{m}\mathbf{X}{\mathbf{X}}^T\\ & =\left[\begin{array}{cccc}\frac{\sum _{i=1}^m(x_{i1}-0)(x_{i1}-0)}{m}& \frac{\sum _{i=1}^m(x_{i1}-0)(x_{i2}-0)}{m}& \cdots & \frac{\sum _{i=1}^m(x_{i1}-0)(x_{in}-0)}{m}\\ \frac{\sum _{i=1}^m(x_{i2}-0)(x_{i1}-0)}{m}& \frac{\sum _{i=1}^m(x_{i2}-0)(x_{i2}-0)}{m}& \cdots & \frac{\sum _{i=1}^m(x_{i2}-0)(x_{in}-0)}{m}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{\sum _{i=1}^m(x_{in}-0)(x_{i1}-0)}{m}& \frac{\sum _{i=1}^m(x_{in}-0)(x_{i2}-0)}{m}& \cdots & \frac{\sum _{i=1}^m(x_{in}-0)(x_{in}-0)}{m}\end{array}\right]\\ & =\left[\begin{array}{cccc}\frac{\sum _{i=1}^m(x_{i1}-{\overline{x}}_1)(x_{i1}-{\overline{x}}_1)}{m}& \frac{\sum _{i=1}^m(x_{i1}-{\overline{x}}_1)(x_{i2}-{\overline{x}}_2)}{m}& \cdots & \frac{\sum _{i=1}^m(x_{i1}-{\overline{x}}_1)(x_{in}-{\overline{x}}_n)}{m}\\ \frac{\sum _{i=1}^m(x_{i2}-{\overline{x}}_2)(x_{i1}-{\overline{x}}_1)}{m}& \frac{\sum _{i=1}^m(x_{i2}-{\overline{x}}_2)(x_{i2}-{\overline{x}}_2)}{m}& \cdots & \frac{\sum _{i=1}^m(x_{i2}-{\overline{x}}_2)(x_{in}-{\overline{x}}_n)}{m}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{\sum _{i=1}^m(x_{in}-{\overline{x}}_n)(x_{i1}-{\overline{x}}_1)}{m}& \frac{\sum _{i=1}^m(x_{in}-{\overline{x}}_n)(x_{i2}-{\overline{x}}_2)}{m}& \cdots & \frac{\sum _{i=1}^m(x_{in}-{\overline{x}}_n)(x_{in}-{\overline{x}}_n)}{m}\end{array}\right]\end{array}$$可见，上式正好是原始样本各指标之间的协方差矩阵，可记为$C$，则第$j$个新基向量上的样本投影方差为：
$${\lambda }_j={\mathbf{e}}_{j}C{\mathbf{e}}_j^T$$由于${\mathbf{e}}_j$是单位向量，所以有：
$$\begin{array}{rl}{\mathbf{e}}_j{\mathbf{e}}_j^T{\lambda }_j& ={\mathbf{e}}_{j}C{\mathbf{e}}_j^T\\ {\mathbf{e}}_j^T{\lambda }_j& =C{\mathbf{e}}_j^T\\ {\lambda }_j{\mathbf{e}}_j^T& =C{\mathbf{e}}_j^T\end{array}$$说明新基各维度上的方差正好是原样本指标中心化后的协方差矩阵特征值。关于矩阵特征值及特征向量的概念和解法，可参考《简述矩阵运算》相关内容。

（四）选出投影方差较大的几个基向量

至此，选出投影方差较大的几个新基向量，就是选出样本协方差矩阵几个较大特征值对应的特征向量，将原始样本向量投影到由这些特征向量组成的新基中即实现了将数据由$n$维到$k$维的降维变换。

三、示例

（一）构造数据

# -*- coding: utf-8 -*-
"""
Created on Wed Apr 27 17:13:05 2022

@author: wquan
"""

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import numpy.matlib as matlib

def randomLine(data, density, sigma=1):
    '''
    沿任意一组多维数据给出各数据点周围各随机离散点的各维度坐标并绘图。
    
    Parameters
    ----------
    data : list of np.arange
        一个由n个序列组成的列表，代表一组n维样本数据，其中各序列依次代表该样本各维度上坐标。
    density : int
        密度参数，代表沿线每个点上周围的随机离散点数量。
    sigma : float
        标准差，反映离散点的中心偏离度。
        
    Returns
    -------
    Rdata : np.arange
        一个n行数组，其中第i行代表各随机离散的点在第i个维度上的各坐标。默认是1。
    
    Examples
    --------
    import numpy as np
    Value = np.linspace(0, 8 * np.pi, 100)
    z = np.linspace(0, 4, 100)
    x = z * np.sin(Value)
    y = z * np.cos(Value)
    rData = randomLine([x,y,z], 10, 0.2)
    '''
    # 验证入参数据各维度坐标数量须一致
    i = 1
    while i < len(data):
        if len(data[i-1]) != len(data[i]):
            raise Exception("the lenth of each arange in the list must equal")
        else:
            i += 1
    Mdata = np.matrix(np.vstack(data)) # 将入参数据矩阵化
    n, m = np.shape(Mdata)             # 获取入参数据的维度和数量
    MRdata = matlib.zeros((n, 1))
    for i in range(m):                 # 对于入参数据中的每个点
        # 为每个入参点生成一团中心在原点的随机点，即density个n维坐标向量
        nr = sigma * matlib.randn(n, density)
        MRdata = np.c_[MRdata, (Mdata[:,i]*np.matlib.ones((1,density)) + nr)]
    # 将最终随机坐标矩阵删除第一列后转成list格式
    rData = np.delete(MRdata,0,1).tolist()

    if n == 2:                                # 如果入参数据是二维
        # 准备画布
        plt.figure()                          # 创建一个新绘图窗口
        plt.axis('scaled')                    # 横纵坐标尺度一致（即使在窗口缩放时）
        x_min = min(rData[0])-1               # X轴下限
        x_max = max(rData[0])+1               # X轴上限
        y_min = min(rData[1])-1               # Y轴下限
        y_max = max(rData[1])+1               # Y轴上限
        all_min = min(x_min, y_min)           # 全局下限
        all_max = max(x_max, y_max)           # 全局上限
        plt.xlim([all_min, all_max])          # 设置x轴坐标系范围
        plt.ylim([all_min, all_max])          # 设置y轴坐标系范围
        # 绘图
        plt.plot(data[0],data[1])             # 绘制沿线
        plt.plot(rData[0], rData[1], 'ro')    # 绘制沿线离散点
    
    if n == 3:                                # 如果入参数据是三维
        # 准备画布
        fig = plt.figure()                    # 获得一个图像对象
        ax = Axes3D(fig)                      # 获得图像对象中的三维坐标系对象
        plt.axis('scaled')                    # 限定坐标系尺度比例固定
        x_min = min(rData[0])-1               # X轴下限
        x_max = max(rData[0])+1               # X轴上限
        y_min = min(rData[1])-1               # Y轴下限
        y_max = max(rData[1])+1               # Y轴上限
        z_min = min(rData[2])-1               # Z轴下限
        z_max = max(rData[2])+1               # Z轴上限
        all_min = min(x_min, y_min, z_min)    # 全局下限
        all_max = max(x_max, y_max, z_max)    # 全局上限
        ax.set_xlim(all_min, all_max)         # 设置x轴坐标系范围
        ax.set_ylim(all_min, all_max)         # 设置y轴坐标系范围
        ax.set_zlim(all_min, all_max)         # 设置z轴坐标系范围
        # 绘图
        ax.plot(data[0], data[1], data[2])           # 绘制沿线
        ax.plot(rData[0], rData[1], rData[2], 'ro')  # 绘制沿线离散点
    
    return rData

# 准备数据
x = np.arange(-10, 10, 1)
y = 2*x+3
z = 0*x+2                                # 在三维空间中生成一条平行于XY平面的直线
rData = randomLine([x,y,z], 10)          # 生成200个在上述直线周围的随机点
mData = np.matrix(rData)                 # 矩阵化

上述代码即可得到第一章图中数据。

（二）降维变换

将构成出的试验样本数据rData按下述代码做操作：

import numpy as np
import numpy.matlib as matlib
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 各维数据减去该维均值使数据中心化
mData = np.matrix(rData)                 # 矩阵化
n, m = np.shape(mData)
cData = mData - mData*matlib.ones((m,1))/m

# 求各维的协方差矩阵
covarianceMatrix = (cData*cData.T)/m

# 求协方差矩阵的特征值和特征向量
coLambdas, coVectors = np.linalg.eig(covarianceMatrix)

# 取出对应于最大特征值的特征向量作为新基向量
E1 = coVectors[:,coLambdas.tolist().index(max(coLambdas))]

# 将原数据投影到新基向量上实现由3维数据降为1维数据
newData = E1.T * mData

# 绘图
newX = newData.tolist()[0]               # 将新生成的1维数据作为X轴坐标
newY = (0 * newData).tolist()[0]         # 将Y轴坐标都置为0
newZ = (0 * newData).tolist()[0]         # 将Z轴坐标都置为0
fig = plt.figure()                       # 获得一个图像对象
ax = Axes3D(fig)                         # 获得图像对象中的三维坐标系对象
plt.axis('scaled')                       # 限定坐标系尺度比例固定
all_min = min(newX)-1                    # 全局下限
all_max = max(newX)+1                    # 全局上限
ax.set_xlim(all_min, all_max)            # 设置x轴坐标系范围
ax.set_ylim(all_min, all_max)            # 设置y轴坐标系范围
ax.set_zlim(all_min, all_max)            # 设置z轴坐标系范围
ax.plot(newX, newY, newZ, 'ro')          # 绘制沿线离散点

即得到降维后的效果：

各样本都只有一个指标维度了。