机器学习数学基础：从奇异值分解 SVD 看 PCA 的主成分-技术圈

今天我们来看一个在数据分析和机器学习领域中常用的降维方法，即主成分分析（PCA）。它是探索性数据分析（EDA）和机器学习算法对数据的基本处理方法。

1引言

首先，我们来看一下机器学习中数据的表示形式。如下图所示，一般用一个矩阵表示数据（特别是表格型数据），

可以是一个目标向量，表示数据的标签。

为什么要对降维呢？我们可以找到多个理由，至少可以这么看

嫌长得太胖穿得太多，一下看不清。让它减减肥，好观察啊。

用行话说就是，特征太多，让人迷失，需要精简，才能在屏幕上展现数据的整体样貌。

嫌长得没啥特点，没有棱角，总之看了没留下啥印象。

用行话可以说是，特征不典型，特征之间线性相关，不利于对数据的后续处理和分析。

总之，嫌数据矩阵的列多并且不好。那怎么让它脱胎换骨，变得更好更美呢？！

2主成分分析

今天主要看的 PCA 方法就是试图在原始高维数据矩阵中找到特征的线性组合以构造更具代表性的数据特征表示形式，在降维的同时让数据更有精气神。

降维问题本身可以看作最优化问题，但本篇主要是从奇异值分解的角度来解读 PCA，因此对于降维问题的描述不作详细展开。简而言之，PCA 降维的目标为，

一方面为了减少数据的特征数，因此要挑选出最具代表性的一些特征。
保持特征的可分性，即在原来空间中有明显差异的数据在降维后也希望尽量保持差异。换句话说，保持数据间的线性结构，比如在原始空间中离得远的数据在降维后也希望离得远，近的自然也一样。

我们结合一个简单例子来阐述一下 PCA。如上图所示，给你一百个数据，每个数据只用两个特征（比如身高和体重），这里为了可视化，分别将它们看成和坐标。

PCA 的目的之一是从这些数据自身出发，找出新的相互正交的坐标轴（将它们看成新的特征），希望这些特征之间不线性相关。这些正交轴也称为数据特征空间中的主成分（PC），最后将数据点投影到这些 PC 上得到新的坐标表示。

我们假设只用一个特征来刻画上面这批数据，也就是说只需要找出一个坐标轴，即一个 PC。这些数据坐标往这个 PC 投影得到新的坐标，而且每个数据只有一个坐标了。那么你怎么选这个 PC 呢？

为了让新坐标尽量能够区分开不同的数据，我们希望沿着这个 PC 能捕获到数据中最大的差异。看下面的动图，我们假设 PC 经过数据集的中心点，然后这些数据点往不同的 PC 上投影，最终都投到一条直线上，那么哪条线将这些数据分开得最好呢？

相信大家能感觉出哪个方向好。但是问题来了，怎么用数学公式来量化不同方向的好差呢？

对了，可以用方差。再加上如果考虑多个 PC 之间的相关性的话，那就是协方差矩阵。

.协方差矩阵

在机器学习中，计算两个特征、（都看成列向量的形式）的协方差公式为

\operatorname{Cov}(\mathbf{x}, \mathbf{y})=\frac{\sum_{{i}=1}^{m}\left(\mathbf{x}_{{i}}-\overline{\mathbf{x}}\right)\left(\mathbf{y}_{{i}}-\overline{\mathbf{y}}\right)}{m-1},

式中表示个样本，上面除以是为了无偏估计。

而由两个特征、定义的协方差矩阵为，

\left[\begin{array}{ll} \operatorname{Cov}(\mathbf{x}, \mathbf{x}) & \operatorname{Cov}(\mathbf{x}, \mathbf{y}) \\ \operatorname{Cov}(\mathbf{y}, \mathbf{x}) & \operatorname{Cov}(\mathbf{y}, \mathbf{y}) \end{array}\right].

如果有多个特征，记为，，... ，，它们的协方差矩阵为，

\left[\begin{array}{ccccc} \operatorname{Cov}(\mathbf{X}_1,\mathbf{X}_1) & \operatorname{Cov}(\mathbf{X}_1,\mathbf{X}_2) & \ldots & \operatorname{Cov}(\mathbf{X}_1, \mathbf{X}_{n}) \\ \operatorname{Cov}(\mathbf{X}_2,\mathbf{X}_1) & \operatorname{Cov}(\mathbf{X}_2,\mathbf{X}_2) & \ldots & \operatorname{Cov}(\mathbf{X}_2, \mathbf{X}_{n}) \\ \vdots & \vdots & \ddots & \vdots \\ \operatorname{Cov}(\mathbf{X}_{n}, \mathbf{X}_1) & \operatorname{Cov}(\mathbf{X}_{n}, \mathbf{X}_2) & \ldots & \operatorname{Cov}(\mathbf{X}_{n}, \mathbf{X}_{n}) \end{array}\right].

协方差矩阵衡量的是数据的不同坐标分量两两之间一起变化的程度。如果它是一个对角矩阵，说明特征之间没有线性相关性，这也是 PCA 所追求的目标之一。

简化形式

如果将这些特征先分别零中心化，即每列元素减去该列的平均值，再按列组装成矩阵，那么此时的协方差矩阵可以简写为，

\mathbf{C}=\frac{\mathbf{X}^{\top} \mathbf{X}}{m-1}

接下来我们用两种方法来得到 PC。

.特征分解

由于协方差矩阵是一个半正定对称矩阵，因此我们可以对它作特征分解，

\begin{array}{lll} \mathbf{C}\!\!\!&=\dfrac{\mathbf{X}^{\top} \mathbf{X}}{m-1}\\[1em]&=\mathbf{Q \Lambda} \mathbf{Q}^{-1} \end{array}

可以看到，协方差矩阵被分解为由特征向量按列组成的矩阵以及由特征值组成的对角矩阵。而这些特征向量就是我们要找的 PC。

这里需要注意的是，这些特征向量和特征值不是随便排列的，例如，特征值在对角线上是按从大到小排列的。

然后，我们可以将数据投影到这些 PC 上去，即。此时，我们再计算一下新特征间的协方差矩阵，

\begin{array}{lll} \mathbf{C}_{\mathbf{Y}}\!\!\!&=\dfrac{\mathbf{Y}^{\top} \mathbf{Y}}{m-1}\\[1em]&=\dfrac{\mathbf{Q}^{\top}\mathbf{X}^{\top} \mathbf{XQ}}{m-1}\\[1em]&=\mathbf{Q}^{\top}\mathbf{Q \Lambda} \mathbf{Q}^{-1}\mathbf{Q}\\[1em]&=\mathbf{\Lambda} \end{array}

新特征的协方差矩阵是对角矩阵，而且，对角线上的方差是从大到小排列的，这重要吗？是的，不要忘了，我们还要降维，因此将数据点投影到前个 PC 上，从而得到数据的新的坐标，

\mathbf{X\;\!Q}_k=\mathbf{X}\!\left[\mathbf{q}_1\;,\mathbf{q}_2\;\!,\;\!\cdots\;\!,\;\!\mathbf{q}_k\right]

而且，这些保留的前个 PC 上的方差是从大到小排列的。回到上面那个平面数据点的例子，如果最后只保留一个 PC 的话，那就是导致方差最大的那个方向了。

.奇异值分解

我们也可以用奇异值分解来计算 PC，但不是分解协方差矩阵，而是分解特征矩阵。

先对作 SVD 分解，得

\mathbf{X}=\mathbf{U} \mathbf{\Sigma V}^{\top}

然后代入协方差矩阵，

\begin{array}{lll} \mathbf{C}\!\!\!&=\dfrac{\mathbf{X}^{\top} \mathbf{X}}{m-1}\\[1em]&=\dfrac{\left(\mathbf{V} \mathbf{\Sigma U}^{\top}\right)\left(\mathbf{U \Sigma V}^{\top}\right)}{m-1}\\[1em]&=\dfrac{\mathbf{V \Sigma}^{2} \mathbf{V}^{\top}}{m-1}\\[1em]&=\dfrac{\mathbf{V \Sigma}^{2} \mathbf{V}^{-1}}{m-1} \end{array}

现在的 PC 是矩阵的列。将数据点投影到前个 PC 上，从而得到数据的新的坐标，

\mathbf{X\;\!V}_k=\mathbf{X}\!\left[\mathbf{v}_1\;,\mathbf{v}_2\;\!,\;\!\cdots\;\!,\;\!\mathbf{v}_k\right].

跟上面特征分解的协方差矩阵比较可得，

\begin{array}{lll} \mathbf{\Lambda}\!\!\!&=\dfrac{\mathbf{\Sigma}^{2}}{m-1} \end{array}

但这里好像没有用到它们啊。其实不然，一般来说，大佬是在幕后操作的，这里处在 C 位的奇异值或者特征值们也一样。

要知道，我们最后得到的 PC 都是按大佬们的大小来排座位的。

.换个角度

除了从方差角度考虑外，也可以从数据点到 PC 轴的垂直距离来考虑问题。

你可以想象下图中的黑线是实心杆，而每条红线表示弹簧。弹簧的能量与它的长度平方成正比（物理学上称为胡克定律），因此，如果所有弹簧长度平方之和达到最小时，杆子将处于平衡状态。下图模拟了达到平衡状态的过程，

这里其实想说的是，PCA 方法可以从多个角度去理解，那么对应的数学问题也会有所不同，但这个不是本篇的重点。

3降维实战

import numpy as np
import seaborn as sns

.小型数据集

iris = sns.load_dataset('iris')

先看一眼数据，显示前五行。

iris.head()

	sepal_length	sepal_width	petal_length	petal_width	species
0	5.1	3.5	1.4	0.2	setosa
1	4.9	3.0	1.4	0.2	setosa
2	4.7	3.2	1.3	0.2	setosa
3	4.6	3.1	1.5	0.2	setosa
4	5.0	3.6	1.4	0.2	setosa

Iris 鸢尾花数据集是一个经典数据集。

数据集内包含 3 类共 150 条记录，每类各 50 个数据，每条记录都有 4 项特征：花萼长度、花萼宽度、花瓣长度、花瓣宽度，可以通过这 4 个特征预测鸢尾花卉的品种。上面五列对应的特征说明，

sepal_length 花萼长度
sepal_width 花萼宽度
petal_length 花瓣长度
petal_width 花瓣宽度
species 品种（共 setosa, versicolour, virginica 三类）

.特征矩阵

我们把最后一列去掉，只取前四列，可以拿来预测最后一列。因此，特征矩阵定义如下，

X = iris.values[:,:-1]
X = X.astype(np.float64)

# 查看前五行数据
X[0:5]

array([[5.1, 3.5, 1.4, 0.2],
       [4.9, 3. , 1.4, 0.2],
       [4.7, 3.2, 1.3, 0.2],
       [4.6, 3.1, 1.5, 0.2],
       [5. , 3.6, 1.4, 0.2]])

# 1、对特征矩阵零中心化
X -= np.mean(X, axis=0)

看一下协方差矩阵

cov_XTX = X.T@X
cov_XTX

array([[102.16833333,  -6.32266667, 189.873     ,  76.92433333],
       [ -6.32266667,  28.30693333, -49.1188    , -18.12426667],
       [189.873     , -49.1188    , 464.3254    , 193.0458    ],
       [ 76.92433333, -18.12426667, 193.0458    ,  86.56993333]])

.降维

# 2、对特征矩阵进行 SVD 分解
U, Sigma, VT = np.linalg.svd(X, full_matrices=False)

U.shape, Sigma.shape, VT.shape

((150, 4), (4,), (4, 4))

由于上面我们设置参数 full_matrices=False，所以得到的矩阵大小为而不是。

我们来构造新的特征，但要注意，上面函数 np.linalg.svd 返回的矩阵 VT 是指。所以，新的特征矩阵对应的代码为，

Y = X@VT.T

再看一下这个新特征矩阵的协方差矩阵

cov_YTY = Y.T@Y
cov_YTY

array([[ 6.30008014e+02, -8.12751455e-14, -1.46101800e-14,
        -6.74844752e-14],
       [-8.12751455e-14,  3.61579414e+01, -8.02807470e-15,
         9.98660131e-15],
       [-1.46101800e-14, -8.02807470e-15,  1.16532155e+01,
        -8.45168460e-16],
       [-6.74844752e-14,  9.98660131e-15, -8.45168460e-16,
         3.55142885e+00]])

非对角线上的协方差都是极小了吧，当零看，因此是个对角矩阵。