主成分分析(PCA)与矩阵的奇异值分解(SVD)

2020-03-02 本文已影响0人 Ice_spring

文章首发于本人CSDN博客，原文地址：ice_spring_的主页
文章内容概要：

矩阵奇异值分解SVD
主成分分析PCA及其应用
SVD与PCA之间的关系

1 矩阵奇异值分解SVD

1.1 矩阵奇异值分解的数学原理

在关于SVD(Singular Value Decomposition)的讲解中将涉及稍微多一点的数学推导。
定义：设 $A$ 是秩为 $r$ 的 $m×n$ 矩阵， $n$ 阶对称方阵 $A^TA$ 的特征值为 $\lambda_i$ $(i=1,2,...,n)$ ，且有

$\lambda_1\geq\lambda_2\geq...\geq\lambda_r>0，\lambda_{r+1}=\lambda_{r+2}=...=\lambda_{n}=0$

则称

$\sigma_i=\sqrt{\lambda_i}$

为矩阵 $A$ 的奇异值。
奇异值分解定理：设 $A$ 是秩为 $r$ 的 $m×n$ 矩阵，则存在 $m$ 阶正交矩阵 $U$ 和 $n$ 阶正交矩阵 $V$ ，使得

$A=U \Sigma V^T$

其中 $m×n$ 矩阵 $\Sigma$ 的分块形式为

$\Sigma= \left(\begin{array}{ccc} D & 0 \\0 & 0 \\ \end{array}\right)$

$D$ 为A的奇异值 $\sigma_i$ 构成的对角矩阵。
下面我们来证明这个定理，并从这个过程中理解 $U， \Sigma， V^T$ 代表的含义。
证明： $A^TA$ 是一个实对称矩阵，且有 $rank(A)=rank(A^TA)$ ，故存在 $n$ 阶正交矩阵 $V$ ，使得

$V^TA^TAV=diag(\lambda_1,\lambda_1,...,\lambda_r,0,...,0)=\Sigma^2$

其中

$V=(v_1,v_2,...,v_r,|v_{r+1},...,v_n)=(V_1,V_2)$

于是

$V_1^TA^TAV_1=D^2,V_2^TA^TAV_2=0$

因此有

$D^{-1}V_1^TA^TAV_1D^{-1}=E$

令 $U_1=AV_1D^{-1}$ ，则上式即为 $U^TU=E$ ，故 $U$ 的列向量是单位正交向量组，将 $U_1$ 扩充为 $U=(U_1|U_2)$ ，则有

$U\Sigma V^T=(U_1,U_2) \left(\begin{array}{ccc} D & 0 \\0 & 0 \\ \end{array}\right) \left(\begin{array}{ccc} V_1^T \\ V_2^T \\ \end{array}\right) =U_1DV_1^T=A$

证毕。
下面我们来探索一下奇异值分解中各个量的含义，首先：

$||Av_i||^2=(Av_i)^T(Av_i)=v_i^T(A^TAv_i)=v_i^T\lambda_iv_i=\lambda_i||v_i||^2$

于是

$||Av_i||=\sigma_i$

这表示 $A$ 对矩阵 $V$ 的分量作用后，得到的就是相应的奇异值。
另外， $\left\{Av_i|Av_i\neq0\right\}$ 是 $A$ 的列向量空间的一组正交基，证明如下：
首先 $(Av_i)^T(Av_j)=v_i^T(A^TAv_i)=0$ ，这证明了其正交性。另外，设

$A=(\alpha_1,\alpha_2,...,\alpha_n)$

向量

$y=Ax=x_1\alpha_1+x_2\alpha_2+...+x_n\alpha_n$

由于 $\left\{v_i\right\}$ 是 $n$ 维空间的一组标准正交基，故

$x=k_1v_1+k_2v_2+...+k_nv_n$

于是

$y=Ak_1v_1+Ak_2v_2+...+Ak_nv_n=k_1Av_1+k_2Av_2+...+k_nAv_n=\\k_1Av_1+k_2Av_2+...+k_rAv_r+0+...+0$

这表明，任何一个可由 $A$ 的列向量组线性表出的向量同样可以由 $\left\{Av_i|Av_i\neq0\right\}$ 线性表出。故 $\left\{Av_i|Av_i\neq0\right\}$ 是 $A$ 的列向量空间的一组正交基。
又由于 $||Av_i||=\sigma_i$ ，则 $\left\{\frac{Av_i}{\sigma_i}|Av_i\neq0\right\}$ 其实就是A的列空间的一组标准正交基，而 $U_1$ 正是由这组标准正交基组成的矩阵。

1.2 SVD分解实例

设有矩阵
$A= \left( \begin{array}{ccc} 1 & 2 \\ 3 & 4 \\ 5 & 6 \end{array} \right)$

试对其进行奇异值分解。
这里就不用手计算了(这个矩阵编的，懒得算了)，我们使用 $numpy$ 和 $scipy$ 来帮助计算，同样也用其中的矩阵乘法来验证分解的正确性。

import numpy as np
from scipy.linalg import svd
A = np.array([
    [1,2],
    [3,4],
    [5,6]
])
U, s, VT = svd(A)

注意中间的我用的是小写 $s$ ，因为 $scipy$ 计算 $SVD$ 分解结果的 $Sigma$ 矩阵返回的是非零奇异值列表，而不是真正的 $Sigma$ 矩阵。

SVD分解

验证分解的正确性即验证
$U.Sigma.VT=A$

#先将s转化为Sigma矩阵
Sigma = np.zeros(A.shape)
for i in range(len(s)):
    Sigma[i][i] = s[i]
#print(Sigma)
# 验证
print(U.dot(Sigma).dot(VT))

验证

可见分解是正确的。

2 主成分分析

主成分分析(Principal Component Analysis )是迄今为止机器学习领域最流行的降维算怯。它先是找出最接近数据特征的超平面，然后将数据投影其上，以达到降维目的。这样说还是有点抽象，下面一步一步来讲解PCA的数学原理，用数学知识推导一遍PCA。

2.1 向量与內积

首先来看一个中学就学过的向量內积运算：

內积

两个向量的內积定义为：

$a • b = |a||b|cosθ$

这样的內积有什么含义呢？下面将內积换一种写法：

$a • b =( |a|cosθ)|b|$

可以看到向量 $a$ 与 $b$ 的內积结果为 $a$ 在 $b$ 上投影的长度再乘以 $b$ 的模，如果让 $b$ 的模为 $1$ ，则：

$a • b =|a|cosθ$

即当 $b$ 模为 $1$ 时， $a$ 与 $b$ 的內积就是 $a$ 在 $b$ 上投影的长度。
推广至 $n$ 维空间，这样，当 $b$ 为某组单位正交基中的一个时， $a$ 与 $b$ 的內积就是 $a$ 在该单位正交基上投影的长度。也就是相应基方向上的坐标。

2.2 基、坐标及其变换的矩阵表示

下面继续讨论，为了方便，我们以二维平面中的点或向量为例。

坐标
在上图中，假设点(向量)的坐标为，这表示点(向量)在轴的投影距离原点为个单位长度，在轴的投影距离原点个单位长度，可见坐标其实是被我们赋予了含义的，显式表示这层含义，点可以表示为：

$2 \begin{pmatrix} 1 \\ 0 \end{pmatrix} +4\begin{pmatrix} 0 \\ 1 \end{pmatrix}$

二维平面中的所有坐标 $(x,y)$ 都可以表示为：

$x\begin{pmatrix} 1 \\ 0 \end{pmatrix}+y\begin{pmatrix} 0 \\ 1 \end{pmatrix}$

这里 $\begin{pmatrix} 1 \\ 0 \end{pmatrix}$ 和 $\begin{pmatrix} 0 \\ 1 \end{pmatrix}$ 就是二维平面的一组基，为了表示二维平面的一个向量，首先要确定一组基，然后给出在基所在的各个直线上的投影值，也就是坐标，不过由于平时见到的几乎都是平面直角坐标系中，于是都默认 $\begin{pmatrix} 1 \\ 0 \end{pmatrix}$ 和 $\begin{pmatrix} 0 \\ 1 \end{pmatrix}$ 为基。
实际上，任意两个线性无关的向量都可以作为二维平面的一组基。不过为了研究问题的方便，一般我们希望基向量的模可以是 $1$ ，而且如果是一组基，它们之间最好是两两正交的，因为投影得到各个基方向的坐标的投影方式都是垂直投影，这样的基有非常良好的性质。从前面对內积的探讨中也可以看出，选用模为 $1$ 的基向量，我们的原向量坐标只需与新选取的基做內积运算就可以得到该向量在新基下的坐标，理解了这一点对理解下面推导PCA原理非常重要。
如我们选取标准正交基 $(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})^T$ 和 $(- \frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}})^T$ ，则 $(2,4)$ 在这组基下的坐标不难通过內积计算得到 $(\frac{6}{\sqrt{2}},\frac{2}{\sqrt{2}})$ ，图示如下：

坐标变换

为了表示方便，我们使用矩阵这一数学工具，将坐标变换表示为矩阵形式：

$(2, 4) \begin{pmatrix} \frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \end{pmatrix} = (\frac{6}{\sqrt{2}},\frac{2}{\sqrt{2}})$

而且还能同时计算多个坐标变换后的结果：

$\left( \begin{array}{ccc} 4 & 1 \\ 2 & 3 \\ 1 & 2 \\ \end{array} \right) \begin{pmatrix} \frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \end{pmatrix} =\begin{pmatrix} \frac{5}{\sqrt{2}} & -\frac{3}{\sqrt{2}} \\ \frac{5}{\sqrt{2}} & \frac{1}{\sqrt{2}} \\ \frac{3}{\sqrt{2}} & \frac{1}{\sqrt{2}} \end{pmatrix}$
其中等号左侧左边矩阵的行向量为直角坐标系下的坐标，基(列向量)拼接为变换矩阵，乘法的结果矩阵就是原直角坐标在新基下变换后的坐标。如果我们不取这组基中全部的基向量(这里二维情形我们只取一个)，这个乘法(內积)操作就相当于只在这个基上做了投影：

$\left( \begin{array}{ccc}4 & 1 \\ 2 & 3 \\ 1 & 2 \\ \end{array} \right) \begin{pmatrix} \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{pmatrix} =\begin{pmatrix} \frac{5}{\sqrt{2}} \\ \frac{5}{\sqrt{2}} \\ \frac{3}{\sqrt{2}} \end{pmatrix}$

如果把各个坐标看做是特征矩阵中的某个数据项，这便达到了降维的目的，上面的例子是把二维特征降维到了一维。(注意：与数学中维度考虑秩不同，我们这里说的特征维度指的是特征的个数。)
更一般地，设有特征矩阵 $X_{m×n}$ ，其中 $m$ 表示数据项的个数， $n$ 表示每个数据项拥有的特征个数，想将其降维为 $Y_{m×r}$ ，只需要乘以行数等于 $n$ 列数等于 $r$ 的变换矩阵 $P_{n×r}$ ：

$X_{m×n}•P_{n×r}=Y_{m×r}$

注意这里的 $r$ 可以小于 $n$ ， $r$ 的大小就决定了原特征矩阵变换后的特征维度。也就是说我们可以把原本高维的特征变换到低纬度空间去研究，当然前提是不损失特征包含的信息。
从上述分析中也可以看到，两个矩阵相乘的意义是将左边矩阵中的每一个行向量变换到右边矩阵中以列向量为一组基所表示的空间中去，当然实际上也可以看做对右边矩阵中列向量的变换，不过这里讨论PCA，为了迁就特征矩阵表示上的习惯使用了这一方式。

2.3 PCA降维问题的优化目标

通过上面的讨论我们知道可以选取新的基对数据重新进行表示，当基的数目小于原特征数目时便可以达到数据降维的目的，不过降维后的数据应尽可能保留原数据中的信息，这样的降维才有意义，现在的问题就是我们怎么选取一组基才能将数据降维并尽量保持原数据中的信息呢？要数学化这个问题并不那么简单，为了避免满篇数学公式的推导，我们用一个具体的例子，设有5个数据项2个特征的特征矩阵如下：

$\left( \begin{array}{ccc} 1 & 1 \\1 & 3 \\2 & 3\\4 & 4\\2 & 4 \end{array} \right)$

为了后续处理方便，先将这组数据进行中心化，也就是让每个特征的均值为0，这只需特征的每个数据减去他们当前的均值，这里第一个特征均值为 $2$ ，第二个特征均值为 $3$ ，于是将特征矩阵处理为：

$\left( \begin{array}{ccc} -1 & -2 \\-1 & 0 \\0 & 0\\2 & 1\\0 & 1\\ \end{array} \right)$

这样处理的好处我们会在后面看到。现在将这组数据在二维直角坐标系中画出来：

特征可视化

现在要对这组数据进行特征降维，一个直观的想法是投影到 $x$ 轴，也就是只取数据 $x$ 轴的坐标(只保留一个特征)，不过此时我们发现有数据点重合，原本不同的数据无法再区分出来，这个降维方案是失败的。显然我们应该选择某个投影方向，使得数据点在该方向投影后(降维后)有最大的区分度，也就是投影后数值尽可能在该方向轴上分散。如将其投影到方向 $(1,1)$ (一三象限角平分线)上貌似就是个不错的选择：

投影

在数学上我们有对这种分散程度更精确的刻画，这种刻画就是方差。

2.3.1 方差与协方差

某个特征(特征矩阵的某一列)的方差计算公式为：

$var(x)=\frac{1}{m}\sum_{i=1}^{m}(x_i-\overline{x})^2$

其中 $m$ 为数据项个数， $\overline{x}$ 为他们的样本均值。
需要指出的是，在统计学上，样本方差计算是除以 $m-1$ 的，这是为了保证样本方差是总体方差的无偏估计，不过在机器学习中大样本情况下差别不大，所以本文讨论中都采取除以 $m$ 这种简单的写法。
实际运用中，都先把数据中心化，均值 $\overline{x}=0$ ，这样方差计算将简便的多：

$var(x)=\frac{1}{m}\sum_{i=1}^{m}x_i^2$

于是上述二维降一维的问题就转化为，寻找一个一维基，使得特征数据变换为这个基上的坐标表示后，方差值最大。
那如果是更高维的特征呢，现在考虑三维降二维问题。与之前相同，首先我们希望找到一个方向使得投影后方差最大，这样就完成了第一个方向的选择，接着我们选择第二个投影方向。如果我们还是单纯只选择方差最大的方向，很明显，这个方向与第一个方向应该是“几乎重合在一起”，显然这样的维度是没有用的，因此，应该有其他约束条件。从直观上说，让两个特征尽可能表示更多的原始信息，我们是不希望它们之间存在(线性)相关性的，因为相关性意味着两个特征不是完全独立，存在重叠表示的信息。数学上可以用两个字段的协方差表示其相关性：

$cov(a,b)=E(ab)-EaEb$

由于我们已经让特征 $a$ 和 $b$ 的均值为 $0$ ，于是：

$cov(a,b)=E(ab)=\frac{1}{m}\sum_{i=1}^{m}a_{i}b_{i}$

当协方差为0时，表示两个字段完全不相关。为了让协方差为0，我们选择的第二个基时只能在与第一个基正交的方向上选择。因此最终选择的两个方向一定是正交的。于是三维降二维的最优化问题便转化为寻找两个互相正交的基，使得原始三维特征变换到二维特征后，特征间协方差为 $0$ ，而特征内的方差值则尽可能大。
至此我们不难得到更一般的降维问题的最优化目标：将一组 $n$ 维特征降维为 $r$ 维，其目标是选择一组 $r$ 个单位正交基，使得原始数据变换到这组基上后，各特征间两两协方差为0，而特征内的方差则尽可能大(在正交的约束下，取最大的 $r$ 个方差)。

2.3.2 协方差矩阵及其对角化

在上面的讨论中我们看到，我们要达到的最优化目标与特征间协方差和特征内方差都有密切关系，如何将它们统一表示呢，那就是协方差矩阵，协方差矩阵的 $(i,j)$ 位置是特征 $i$ 与 $j$ 的协方差， $(i,i)$ 位置是特征 $i$ 自身的方差，由于已经事先让各个特征均值为0，无论方差还是协方差都是向量內积的形式，设原特征矩阵经过降维变换后的矩阵为 $Y$ ，此时协方差矩阵其实就是：

$D=\frac{1}{m}Y^TY$

原特征对应的协方差矩阵：
$C=\frac{1}{m}X^TX$

那么 $D$ 与 $C$ 有什么关系呢？下面我们就来计算一下，由于 $Y$ 是特征降维后的矩阵，故 $D$ 一定是对角矩阵(协方差为 $0$ )，于是：

$D=\frac{1}{m}Y^TY = \frac{1}{m}(XP)^T(XP) = \frac{1}{m}P^TX^TXP=P^TCP$

这下事情清楚了，我们要找的一组正交基组成的变换矩阵 $P$ 正是可以使得原协方差矩阵 $C$ 对角化所对应的正交矩阵。而原协方差矩阵 $C=X^TX$ 是一个实对称矩阵，实对称矩阵拥有非常良好的性质：

实对称矩阵一定可以对角化；
实对称矩阵不同特征值对应的特征向量互相正交；
实对称矩阵的k重特征值对应k个线性无关的特征向量；

我们可以通过实对称矩阵对角化理论求出我们需要的变换矩阵。
原特征矩阵 $C$ 对角化后的矩阵 $D$ 的对角线元素正是 $C$ 的特征值(设这些特征值按从大到小的顺序排列)，由上面我们知道这里的特征值就是变换后各特征的方差，其含义就是将原数据投影到该特征值对应的特征向量方向后，新数据的重要程度，方差越大则数据越重要。
一般我们只取 $P$ 的前 $r$ 列 $P_{n×r}$ ，代表变换后前 $r$ 个最重要的特征所在的方向，而丢弃剩下的列

$X_{m×n}•P_{n×r}=Y_{m×r}$

这便达到了降维的效果。
以上就是PCA的原理。怎么样，够清晰吧。

2.4 PCA降维在手写数字数据集应用实例

接下来，用一个机器学习实例来看一下PCA的强大之处。
我们使用手写数字数据集，这个数据集共有 $42000$ 个样本，每个样本都拥有784个特征(每个样本就是一张 $28×28$ 的手写数字图像)。

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

data = pd.read_csv("digit recognizor.csv")
X = data.iloc[:,1:]
y = data.iloc[:,0]
X.shape # (42000,784)

可视化一部分样本(40个)

def plot(mnist):
    fig, axes = plt.subplots(4,10,figsize=(10,4),subplot_kw={'xticks':[],'yticks':[]},
                            gridspec_kw=dict(hspace=0.1,wspace=0.1))
    for i,ax in enumerate(axes.flat):
        ax.imshow(mnist[i].reshape(28,28),cmap='binary')
    plt.show()
X_mnixt_example = X[:40]
X_mnixt_example = np.array(X_mnixt_example, dtype='int32')# 整理成numpy格式数据
plot(X_mnixt_example)# 调用函数画图

mnist

接着用 $KNN$ 算法训练一遍原数据集，使用交叉验证评价其表现：

# Time using 1.5h
from sklearn.neighbors import KNeighborsClassifier as KNN
cross_val_score(KNN(5), X, y, cv=5, n_jobs=-1).mean() # 98.3%

$KNN$ 虽然简单但却也是非常强大的机器学习算法，不过由于其每次计算往往都要遍历全部样本，这样当数据特别多，数据维度特别大时将非常消耗时间。对于 $KNN$ 这样的分类器，我们希望在提升模型表现和训练模型时间消耗上达到一个均衡。这也正是降维的意义之一。
可以看到 $KNN$ 算法在 $5$ 折交叉验证下在数据集上达到平均 $98.3\%$ 的准确率。但是时间消耗为1.5小时，这是因为数具量实在太大，当然也和 $KNN$ 这种算法有关。同样的数据在随机森林上只需 $10$ 分钟左右。不过这里为了展示PCA的强大之处，使用了 $KNN$ 。
接下来用 $PCA$ 对数据进行降维，为了能有更好的降维效果，首先绘制维度的学习曲线，大致确定一个不至于使模型表现太坏的维度。由于数据有一定的稳定性，在 $RF$ 表现好的数据在 $KNN$ 也极大可能会表现好，为了效率先用 $RF$ 分类器作训练来绘制维度的学习曲线：

score = []
for i in range(1,101,10):
    X_dr = PCA(i).fit_transform(X)
    once = cross_val_score(RFC(n_estimators=10,random_state=0) ,X_dr,y,cv=5).mean()
    score.append(once)
plt.figure(figsize=[20,5])
plt.plot(range(1,101,10),score)
plt.show()

学习曲线

可以发现维度在 $10$ 之后对准确率的提升已经很小了。我们选取降维到 $21(784->21)$ 。这里需要解释的是，在 $RF$ 上准确率仅在 $91\%$ 上下，是因为我们没有对 $RF$ 进行调参，如果增加其中决策树的数量是可以看到准确率在 $93\%$ 以上的，另外，每种模型往往也都有自己更擅长的数据集， $RF$ 也许对这种数据就不太擅长。接下来再在 $KNN$ 上看一下降维数据的表现，由于数据已降维到 $21$ ，相比之前运算量已大幅下降，故使用 $KNN$ 也是可以的：

# Time using 5min
X_dr = PCA(21).fit_transform(X)
X_dr.shape #(42000, 21)
cross_val_score(KNN(5), X_dr, y, cv=5, n_jobs=-1).mean()

准确率

可以看到虽然准确率不如使用全部数据，但也达到了 $97\%$ ，原数据有 $784$ 个维度，而我们只使用了降维的 $21$ 维度数据就有接近原数据的效果！这就是 $PCA$ 的强大之处。

3 PCA与SVD之间的关系

我们经常听到别人说 $PCA$ 的求解用到了 $SVD$ 分解，这是什么意思呢？首先我们来看二者之间的关系。
在 $PCA$ 的求解中
$X_{m×n}•P_{n×r}=Y_{m×r}$
变换矩阵 $P_{n×n}$ 是原特征矩阵的协方差矩阵 $X^TX$ 的正交特征向量按列拼成的矩阵。而在 $SVD$ 分解中
$A=U \Sigma V^T$
矩阵 $V$ 是 $A^TA$ 的正交特征向量按列组成的矩阵。
这下我们知道了，若 $A=X$ ，那么 $SVD$ 分解中的正交矩阵 $V$ 正是PCA中的正交矩阵 $P$ 。这就是二者之间的数学关系。
而 $SVD$ 分解已经有更快的数值解法了，不需要直接计算 $A^TA$ 的特征值和特征向量(对于大矩阵如果这样计算将会有非常非常大的计算量)， $SVD$ 可以不计算协方差矩阵等复杂的计算过程，就直接求出新特征空间和降维后的特征矩阵。所以 $PCA$ 的求解通常都是用 $SVD$ 分解得到 $V(P)$ 矩阵。
例如在 $sklearn$ 中，其封装的 $PCA$ 就是基于 $SVD$ 求解的，而矩阵 $U$ 和 $Σ$ 虽然会被计算出来，但完全不会被用到。

from sklearn.decomposition import PCA
pca_f = PCA(n_components=0.97,svd_solver="full")
pca_f = pca_f.fit(X)
X_f = pca_f.transform(X)
pca_f.explained_variance_ratio_ #0.978

如上面一段代码，在 $n\_components$ 位置输入 $[0,1]$ 之间的浮点数，并且让参数 $svd\_solver ='full'$ (该参数就是指定 $SVD$ 的计算方式的)，表示希望降维后的总解释性方差占比大于 $n\_components$ 指定的百分比，即是说，希望保留百分之多少的信息量。 $PCA$ 按照 $SVD$ 分解过程会自动选出能够让保留的信息量超过我们指定百分比的特征数量。

sklearnPCA

注：本文在解释PCA原理部分的思路参考了知乎文章PCA的数学原理，由于已不知道文章最早的作者是谁，特此说明，侵删。

参考文献：
.
张贤达. 矩阵分析与应用(第二版). 北京：清华大学出版社，2013.
.
PCA的数学原理. 知乎文章. 郑申海：https://zhuanlan.zhihu.com/p/21580949，2016.
.
Aurélien Géron. Hands-On Machine Learning with Scikit-Learn and TensorFlow. Dimensionality Reduction

转载请注明出处，谢谢合作！