引言

在互联网大数据场景下，我们经常需要面对高维数据，在对这些数据做分析和可视化的时候，我们通常会面对「高维」这个障碍。在数据挖掘和建模的过程中，高维数据也同样带来大的计算量，占据更多的资源，而且许多变量之间可能存在相关性，从而增加了分析与建模的复杂性。

我们希望找到一种方法，在对数据完成降维「压缩」的同时，尽量减少信息损失。由于各变量之间存在一定的相关关系，因此可以考虑将关系紧密的变量变成尽可能少的新变量，使这些新变量是两两不相关的，那么就可以用较少的综合指标分别代表存在于各个变量中的各类信息。机器学习中的降维算法就是这样的一类算法。

主成分分析（Principal Components Analysis，简称PCA）是最重要的数据降维方法之一。在数据压缩消除冗余和数据噪音消除等领域都有广泛的应用。本篇我们来展开讲解一下这个算法。

（本篇降维算法部分内容涉及到机器学习基础知识，没有先序知识储备的宝宝可以查看ShowMeAI的文章 图解机器学习 | 机器学习基础知识。

1.PCA与最大可分性

对于X=[x1x2...xn]X = \begin {bmatrix} x_1 \\ x_2 \\ ... \\ x_n \end{bmatrix}X=?????x1?x2?...xn????????，X∈RnX \in R^nX∈Rn。我们希望XXX从nnn维降到n′n^{'}n′维，同时希望信息损失最少。比如，从n=2n = 2n=2维降到n′=1n^{'} = 1n′=1。

左图为一个典型的例子，假如我们要对一系列人的样本进行数据降维（每个样本包含「身高」「体重」两个维度）。右图我们既可以降维到第一主成分轴，也可以降维到第二主成分轴。

哪个主成分轴更优呢？从直观感觉上，我们会认为「第一主成分轴」优于「第二主成分轴」，因为它比较大程度保留了数据之间的区分性（保留大部分信息）。

对PCA算法而言，我们希望找到小于原数据维度的若干个投影坐标方向，把数据投影在这些方向，获得压缩的信息表示。下面我们就一步一步来推导一下PCA算法原理。

2.基变换

先来复习一点点数学知识。我们知道要获得原始数据XXX新的表示空间YYY，最简单的方法是对原始数据进行线性变换（也叫做基变换）Y=PXY = PXY=PX。其中，XXX是原始样本，PPP是基向量，YYY是新表达。

数学表达为：

[p1p2?pr]r×n[x1x2?xm]n×m=[p1x1p1x2?p1xmp2x1p2x2?p2xm????prx1prx2?prxm]r×m\begin{bmatrix} p_1 \\ p_2 \\ \vdots \\ p_r \end{bmatrix}_{r \times n} \begin{bmatrix} x_1 & x_2 & \cdots & x_m \end{bmatrix}_{n \times m} = \begin{bmatrix} p_1 x_1 & p_1 x_2 & \cdots & p_1 x_m \\ p_2 x_1 & p_2 x_2 & \cdots & p_2 x_m \\ \vdots & \vdots & \ddots & \vdots \\ p_r x_1 & p_r x_2 & \cdots & p_r x_m\end{bmatrix}_{r\times m}???????p1?p2??pr?????????r×n?[x1??x2????xm??]n×m?=???????p1?x1?p2?x1??pr?x1??p1?x2?p2?x2??pr?x2???????p1?xm?p2?xm??pr?xm?????????r×m?

• 其中pip_ipi?是行向量，表示第iii个基；
• xjx_jxj?是一个列向量，表示第jjj个原始数据记录。

当r<nr < nr<n时，即「基的维度<数据维度」时，可达到降维的目的，即X∈Rn×m→Y∈Rr×mX \in R^{n \times m} \rightarrow Y \in R^{r \times m }X∈Rn×m→Y∈Rr×m。

以直角坐标系下的点(3,2)(3,2)(3,2)为例，要把点(3,2)(3,2)(3,2)变换为新基上的坐标，就是用(3,2)(3,2)(3,2)与第一个基做内积运算，作为第一个新的坐标分量，然后用(3,2)(3,2)(3,2)与第二个基做内积运算，作为第二个新坐标的分量。

上述变化，在线性代数里，我们可以用矩阵相乘的形式简洁的来表示：

[1212?1212][32]=[52?12]\begin{bmatrix}\frac{1}{\sqrt 2} & \frac{1}{\sqrt 2} \\ -\frac{1}{\sqrt 2} & \frac{1}{\sqrt 2} \end{bmatrix} \begin{bmatrix} 3 \\ 2\end{bmatrix} = \begin{bmatrix} \frac{5}{\sqrt 2} \\ - \frac{1}{\sqrt 2} \end{bmatrix}[2?1??2?1??2?1?2?1??][32?]=[2?5??2?1??]

再稍微推广一下，假如我们有m个二维向量，只要将二维向量按列排成一个两行m列矩阵，然后用「基矩阵」乘以这个矩阵，就得到了所有这些向量在新基下的值。例如(1,1)、(2,2)、(3,3)，想变换到刚才那组基上，可以如下这样表示：

[1212?1212][123123]=[224262000]\begin{bmatrix}\frac{1}{\sqrt 2} & \frac{1}{\sqrt 2} \\ -\frac{1}{\sqrt 2} & \frac{1}{\sqrt 2} \end{bmatrix} \begin{bmatrix} 1 & 2 & 3 \\ 1 & 2 & 3\end{bmatrix} = \begin{bmatrix} 2\sqrt 2 & 4\sqrt2 & 6\sqrt2 \\ 0 & 0 & 0 \end{bmatrix}[2?1??2?1??2?1?2?1??][11?22?33?]=[22?0?42?0?62?0?]

3.方差

在本文的开始部分，我们提到了，降维的目的是希望压缩数据但信息损失最少，也就是说，我们希望投影后的数据尽可能分散开。在数学上，这种分散程度我们用「方差」来表达，方差越大，数据越分散。

• 定义方差VarVarVar：对于单一随机变量aaa，Var(a)=1m∑i=1m(ai?μ)2Var(a) = \frac{1}{m} \sum_{i = 1}^m (a_i - \mu)^2Var(a)=m1?∑i=1m?(ai??μ)2
• 对数据做去中心化（方便后面操作）：Var(a)=1m∑i=1mai2Var(a) = \frac{1}{m} \sum_{i = 1}^m a_i ^2Var(a)=m1?∑i=1m?ai2?

Var(a)Var(a)Var(a)表示aaa的取值与其数学期望之间的偏离程度。若Var(a)Var(a)Var(a)较小，意味着aaa的取值主要集中在期望μ\muμ也就是E(a)E(a)E(a))的附近；反之，若Var(a)Var(a)Var(a)较大，意味着aaa的取值比较分散。

我们来看一个具体的例子。假设我们5个样本数据，分别是x1=[11]x_1 = \begin{bmatrix} 1 \\ 1 \end{bmatrix}x1?=[11?]、x2=[13]x_2 = \begin{bmatrix} 1 \\ 3\end{bmatrix}x2?=[13?]、x3=[23]x_3 = \begin{bmatrix} 2 \\ 3\end{bmatrix}x3?=[23?]、x4=[44]x_4 = \begin{bmatrix} 4 \\ 4\end{bmatrix}x4?=[44?]、x5=[24]x_5 = \begin{bmatrix} 2 \\ 4 \end{bmatrix}x5?=[24?]，将它们表示成矩阵形式：

X=[1124213344]X = \begin{bmatrix} 1 & 1 & 2 & 4 & 2 \\ 1 & 3 & 3 & 4 & 4 \end {bmatrix}X=[11?13?23?44?24?]

为了后续处理方便，我们首先将每个字段内所有值都减去字段均值，其结果是将每个字段都变为均值为0。

我们看上面的数据，设第一个特征为aaa，第二个特征为bbb，则某个样本可以写作xi=[ab]x_i = \begin{bmatrix} a \\ b \end {bmatrix}xi?=[ab?] 且特征aaa的均值为2，特征bbb的均值为3。所以，变换后

X=[?1?1020?20011]X = \begin{bmatrix} -1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1 \end{bmatrix}X=[?1?2??10?00?21?01?]

Var(a)=65Var(a ) = \frac{\sqrt 6} {5}Var(a)=56??

Var(b)=65Var(b ) = \frac{\sqrt 6} {5}Var(b)=56??

4.协方差

协方差（Covariance）在概率和统计学中用于衡量两个变量的总体误差。比如对于二维随机变量xi=[ab]x_i = \begin{bmatrix} a \\ b \end{bmatrix}xi?=[ab?]，特征a、ba、ba、b除了自身的数学期望和方差，还需要讨论a、ba、ba、b之间互相关系的数学特征。

定义协方差CovCovCov：

Cov(a,b)=1m∑i=1maibiCov(a, b) = \frac{1}{m}\sum_{i = 1}^m a_i b_iCov(a,b)=m1?i=1∑m?ai?bi?

当Cov(a,b)=0Cov(a, b) = 0Cov(a,b)=0时，变量a、ba、ba、b完全独立，这也是我们希望达到的优化目标。方差是协方差的一种特殊情况，即当两个变量是相同的情况Cov(a,a)=Var(a)Cov(a, a) = Var(a)Cov(a,a)=Var(a)。

5.协方差矩阵

对于二维随机变量xi=[ab]x_i = \begin{bmatrix} a \\ b \end {bmatrix}xi?=[ab?]，定义协方差矩阵C=[Var(a)Cov(a,b)Cov(b,a)Var(b)]C = \begin{bmatrix} Var(a) & Cov(a, b) \\ Cov(b, a) &Var(b)\end{bmatrix}C=[Var(a)Cov(b,a)?Cov(a,b)Var(b)?]。

对于nnn维随机变量 xi=[x1x2?xn]x_{i}=\left[\begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{array}\right]xi?=???????x1?x2??xn?????????

C=[Var(x1)Cov(x1,x2)?Cov(x1,xn)Cov(x2,x1)Var(x2)?Cov(x1,xn)????Cov(xn,x1)Cov(xn,x2)?Var(xn)]C = \begin{bmatrix} Var(x_1) & Cov(x_1, x_2) &\cdots & Cov(x_1, x_n)\\ Cov(x_2, x_1)& Var(x_2) & \cdots & Cov(x_1, x_n)\\ \vdots & \vdots & \ddots & \vdots \\ Cov(x_n, x_1) & Cov(x_n, x_2) & \cdots &Var(x_n) \end{bmatrix}C=???????Var(x1?)Cov(x2?,x1?)?Cov(xn?,x1?)?Cov(x1?,x2?)Var(x2?)?Cov(xn?,x2?)??????Cov(x1?,xn?)Cov(x1?,xn?)?Var(xn?)????????

我们可以看到，协方差矩阵是nnn行nnn列的对称矩阵，主对角线上是方差，而协对角线上是协方差。

我们再来用一个示例对应讲解一下。还是同样的5个样本数据

• x1=[11]x_1 = \begin{bmatrix} 1 \\ 1 \end{bmatrix}x1?=[11?]
• x2=[13]x_2 = \begin{bmatrix} 1 \\ 3\end{bmatrix}x2?=[13?]
• x3=[23]x_3 = \begin{bmatrix} 2 \\ 3\end{bmatrix}x3?=[23?]
• x4=[44]x_4 = \begin{bmatrix} 4 \\ 4\end{bmatrix}x4?=[44?]
• x5=[24]x_5 = \begin{bmatrix} 2 \\ 4 \end{bmatrix}x5?=[24?]

去中心化后表示成矩阵

X=[?1?1020?20011]X = \begin{bmatrix} -1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1 \end{bmatrix}X=[?1?2??10?00?21?01?]

那如果有mmm个样本的话，X=[a1a2?amb1b2?bm]X =\begin{bmatrix} a_1 & a_2 & \cdots &a_m \\ b_1 & b_2 & \cdots & b_m\end{bmatrix}X=[a1?b1??a2?b2?????am?bm??]。对XXX做一些变换，用XXX乘以XXX的转置，并乘上系数1/m1/m1/m：

1mXXT=1m[a1a2?amb1b2?bm][a1b1a2b2??ambm]=[1m∑i=1mai21m∑i=1maibi1m∑i=1maibi1m∑i=1mbi2]\frac{1}{m}XX^T = \frac{1}{m}\begin{bmatrix} a_1 & a_2 & \cdots &a_m \\ b_1 & b_2 & \cdots & b_m\end{bmatrix}\begin{bmatrix} a_1 & b_1 \\ a_2 & b_2 \\ \vdots & \vdots \\ a_m &b_m \end{bmatrix} = \begin{bmatrix} \frac{1}{m} \sum_{i = 1}^m a_i ^2 & \frac{1}{m}\sum_{i = 1}^m a_i b_i \\ \frac{1}{m}\sum_{i = 1}^m a_i b_i& \frac{1}{m} \sum_{i = 1}^m b_i^2 \end{bmatrix}m1?XXT=m1?[a1?b1??a2?b2?????am?bm??]???????a1?a2??am??b1?b2??bm?????????=[m1?∑i=1m?ai2?m1?∑i=1m?ai?bi??m1?∑i=1m?ai?bi?m1?∑i=1m?bi2??]

这正是协方差矩阵！我们归纳得到：设我们有mmm个nnn维数据记录，将其按列排成nnn乘mmm的矩阵XXX，设C=1mXXTC = \frac{1}{m}XX^TC=m1?XXT，则CCC是一个对称矩阵，其对角线分别个各个特征的方差，而第iii行jjj列和jjj行iii列元素相同，表示iii和jjj两个特征之间的协方差。

6.协方差矩阵对角化

再回到我们的场景和目标：

• 现在我们有mmm个样本数据，每个样本有nnn个特征，那么设这些原始数据为XXX，XXX为nnn行mmm列的矩阵。
• 想要找到一个基PPP，使Yr×m=Pr×nXn×mY_{r \times m} = P_{r \times n}X_{n \times m}Yr×m?=Pr×n?Xn×m?，其中$r，达到降维的目的。

设XXX的协方差矩阵为CCC，YYY的协方差矩阵为DDD，且Y=PXY = PXY=PX。

• 我们的目的变为：对原始数据XXX做PCA后，得到的YYY的协方差矩阵DDD的各个方向方差最大，协方差为0。

那么CCC与DDD是什么关系呢？

D=1mYYT=1m(PX)(PX)T=1mPXXTPT=1mP(XXT)PT=PCPT=P[1m∑i=1mai21m∑i=1maibi1m∑i=1maibi1m∑i=1mbi2]PT\begin{aligned} D & =\frac{1}{m} Y Y^{T} \\ & =\frac{1}{m}(P X)(P X)^{T} \\ & =\frac{1}{m} P X X^{T} P^{T} \\ & =\frac{1}{m} P\left(X X^{T}\right) P^{T} \\ & =P C P^{T} \\ & =P\left[\begin{array}{cc} \frac{1}{m} \sum_{i=1}^{m} a_{i}^{2} & \frac{1}{m} \sum_{i=1}^{m} a_{i} b_{i} \\ \frac{1}{m} \sum_{i=1}^{m} a_{i} b_{i} & \frac{1}{m} \sum_{i=1}^{m} b_{i}^{2} \end{array}\right] P^{T} \end{aligned}D?=m1?YYT=m1?(PX)(PX)T=m1?PXXTPT=m1?P(XXT)PT=PCPT=P[m1?∑i=1m?ai2?m1?∑i=1m?ai?bi??m1?∑i=1m?ai?bi?m1?∑i=1m?bi2??]PT?

我们发现，要找的PPP不是别的，而是能让原始协方差矩阵对角化的PPP。

换句话说，优化目标变成了寻找一个矩阵PPP，满足PCPPCP^PCPT是一个对角矩阵，并且对角元素按从大到小依次排列，那么P的前K行就是要寻找的基，用P的前K行组成的矩阵乘以X就使得X从N维降到了K维并满足上述优化条件。

最终我们聚焦在协方差矩阵对角化这个问题上。

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

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

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

由上面两条可知，一个nnn行nnn列的实对称矩阵一定可以找到nnn个单位正交特征向量，设这nnn个特征向量为e1,e2,?,ene_1,e_2,?,e_ne1?,e2?,?,en?，我们将其按列组成矩阵：

E=[e1e2? en]E = \begin{bmatrix} e_1 & e_2 & \cdots \ e_n\end{bmatrix}E=[e1??e2??? en??]

则对协方差矩阵CCC有如下结论：

ETCE=Λ=[λ1λ2?λn]E^T C E = \Lambda = \begin{bmatrix} \lambda_1 \\ & \lambda_2 \\ &&\ddots \\ &&&\lambda_n\end {bmatrix}ETCE=Λ=??????λ1??λ2????λn????????

其中Λ\LambdaΛ为对角矩阵，其对角元素为各特征向量对应的特征值（可能有重复）。 结合上面的公式：

D=PCPTD = PCP^TD=PCPT

其中，DDD为对角矩阵，我们可以得到：

P=ETP = E^TP=ET

PPP是协方差矩阵CCC的特征向量单位化后按行排列出的矩阵，其中每一行都是CCC的一个特征向量。如果设PPP按照Λ\LambdaΛ中特征值的从大到小，将特征向量从上到下排列，则用PPP的前KKK行组成的矩阵乘以原始数据矩阵XXX，就得到了我们需要的降维后的数据矩阵YYY。

7.PCA算法

总结一下PCA的算法步骤：

设有mmm条nnn维数据。

1）将原始数据按列组成nnn行mmm列矩阵XXX

2）将XXX的每一行（代表一个特征）进行零均值化，即减去这一行的均值

3）求出协方差矩阵C=1mXXC=\frac{1}{m}XX^C=m1?XXT

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

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

6）Y=PXY=PXY=PX即为降维到kkk维后的数据

8.PCA代码实践

我们这里直接使用python机器学习工具库scikit-learn来给大家演示PCA算法应用（相关知识速查可以查看ShowMeAI文章AI建模工具速查|Scikit-learn使用指南），sklearn工具库中与PCA相关的类都在sklearn.decomposition包里，最常用的PCA类就是sklearn.decomposition.PCA。

1）参数介绍

sklearn中的PCA类使用简单，基本无需调参，一般只需要指定需要降维到的维度，或者降维后的主成分的方差和占原始维度所有特征方差和的比例阈值就可以了。

下面是sklearn.decomposition.PCA的主要参数介绍：

• n_components：PCA降维后的特征维度数目。
• whiten：是否进行白化。所谓白化，就是对降维后的数据的每个特征进行归一化，让方差都为1，默认值是False，即不进行白化。
• svd_solver：奇异值分解SVD的方法，有4个可以选择的值：{‘auto’,‘full’,‘arpack’,‘randomized’}。

除上述输入参数，还有两个PCA类的成员属性也很重要：

• ① explained_variance_，它代表降维后的各主成分的方差值。
• ② explained_variance_ratio_，它代表降维后的各主成分的方差值占总方差值的比例。

2）代码实例

# 构建数据样本并可视化

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
from sklearn.datasets import make_blobs

# X为样本特征，Y为样本簇类别， 共1000个样本，每个样本3个特征，共4个簇
X, y = make_blobs(n_samples=10000, n_features=3, centers=[[3,3, 3], [0,0,0], [1,1,1], [2,2,2]], cluster_std=[0.2, 0.1, 0.2, 0.2], 
                  random_state =9)
fig = plt.figure()
ax = Axes3D(fig, rect=[0, 0, 1, 1], elev=30, azim=20)
plt.scatter(X[:, 0], X[:, 1], X[:, 2],marker='x')

复制代码

先不降维，只对数据进行投影，看看投影后的三个维度的方差分布，代码如下：

from sklearn.decomposition import PCA
pca = PCA(n_components=3)
pca.fit(X)
print(pca.explained_variance_ratio_)
print(pca.explained_variance_)
复制代码

输出如下：

[0.98318212 0.00850037 0.00831751]
[3.78521638 0.03272613 0.03202212]
复制代码

可以看出投影后三个特征维度的方差比例大约为98.3%：0.8%：0.8%。投影后第一个特征占了绝大多数的主成分比例。现在我们来进行降维，从三维降到2维，代码如下：

pca = PCA(n_components=2)
pca.fit(X)
print(pca.explained_variance_ratio_)
print(pca.explained_variance_)
复制代码

输出如下：

[0.98318212 0.00850037]
[3.78521638 0.03272613]
复制代码

这个结果其实可以预料，因为上面三个投影后的特征维度的方差分别为：[ 3.78521638 0.03272613]，投影到二维后选择的肯定是前两个特征，而抛弃第三个特征。为了有个直观的认识，我们看看此时转化后的数据分布，代码如下：

X_new = pca.transform(X)
plt.scatter(X_new[:, 0], X_new[:, 1],marker='x')
plt.show()
复制代码

从上图可以看出，降维后的数据依然清楚可见之前三维图中的4个簇。现在我们不直接指定降维的维度，而指定降维后的主成分方差和比例，来试验一下。

pca = PCA(n_components=0.9)
pca.fit(X)
print(pca.explained_variance_ratio_)
print(pca.explained_variance_)
print(pca.n_components_)
复制代码

我们指定了主成分至少占90%，输出如下：

[0.98318212]
[3.78521638]
1
复制代码

可见只有第一个投影特征被保留。这也很好理解，我们的第一个主成分占投影特征的方差比例高达98%。只选择这一个特征维度便可以满足90%的阈值。我们现在选择阈值99%看看，代码如下：

pca = PCA(n_components=0.99)
pca.fit(X)
print(pca.explained_variance_ratio_)
print(pca.explained_variance_)
print(pca.n_components_)
复制代码

此时的输出如下：

[0.98318212 0.00850037]
[3.78521638 0.03272613]
2
复制代码

这个结果也很好理解，因为我们第一个主成分占了98.3%的方差比例，第二个主成分占了0.8%的方差比例，两者一起可以满足我们的阈值。最后我们看看让MLE算法自己选择降维维度的效果，代码如下：

pca = PCA(n_components= 'mle',svd_solver='full')
pca.fit(X)
print(pca.explained_variance_ratio_)
print(pca.explained_variance_)
print(pca.n_components_)
复制代码

输出结果如下：

[0.98318212]
[3.78521638]
1
复制代码

可见由于我们的数据的第一个投影特征的方差占比高达98.3%，MLE算法只保留了我们的第一个特征。

更多无监督学习的算法模型总结可以查看ShowMeAI的文章 AI知识技能速查 | 机器学习-无监督学习。