机器学习算法之PCA算法

作者: GiantPandaCV | 来源:发表于2019-11-15 09:54 被阅读0次

前言

在机器学习中降维是我们经常需要用到的算法，在降维的众多方法中PCA无疑是最经典的机器学习算法之一，最近准备撸一个人脸识别算法，也会频繁用到PCA，本文就带着大家一起来学习PCA算法。

前置内容

要学会PCA算法，首先需要了解矩阵分解算法。而矩阵分解算法又分为特征值分解和SVD(奇异值)分解，这两个算法的目的都是提取出一个矩阵最重要的特征。

特征值分解

特征值，特征向量

如果一个向量 $v$ 是矩阵A的特征向量，则一定可以表示为下面的形式：
$Av=\lambda v$
其中 $\lambda$ 是特征向量 $v$ 对应的特征值，一个矩阵的一组特征向量是一组正交向量。这里大概可以理解为矩阵A与向量v相乘实际上就是对向量v进行了一次线性变换(旋转或者拉伸)，而该变换的效果为常数 $\lambda$ 乘以向量v。当我们求特征值和特征向量的时候就是为了求矩阵A能使哪些向量(特征向量)只发生伸缩变换，而变换的成都可以用特征值 $\lambda$ 表示。

特征值与特征向量的几何意义

上面提到一个矩阵乘以一个向量就相当于一个线性变化可能难以理解，我们举个例子：
$M = \left[ \begin{matrix} 3 & 0\\ 0 & 1 \\ \end{matrix} \right]$
然后向量 $X=\left[ \begin{matrix} x\\ y \\ \end{matrix} \right]$
所以M乘以一个向量(x,y)就表示为：
$\left[ \begin{matrix} 3 & 0\\ 0 & 1\\ \end{matrix} \right]$ * $\left[ \begin{matrix} x\\ y \\ \end{matrix} \right]$ = $\left[ \begin{matrix} 3x\\ y \\ \end{matrix} \right]$

我们在平面上化出这个变换过程如下：

在这里插入图片描述同时，注意到上面的矩阵M是对称的，如果不是对称的，例如
那么M*X的变换就可以用下图来表示：

在这里插入图片描述其中蓝色箭头指的是一个最主要的变换方向。如果要描述好一个变换，我们就需要描述好这个变换最主要的变化方向。

特征值分解

对于矩阵A，有一组特征向量v，将这组向量进行正交化单位化，就能得到一组正交单位向量。特征值分解，就是将矩阵A分解为如下式：
$A=Q\sum Q^{-1}$
其中，Q是矩阵A的特征向量组成的矩阵， $\sum$ 则是一个对角阵，对角线上的元素就是特征值。我们来分析一下特征值分解的式子，分解得到的 $\sum$ 矩阵是一个对角矩阵，里面的特征值是由大到小排列的，这些特征值所对应的特征向量就是描述这个矩阵变换方向（从主要的变化到次要的变化排列）。
当矩阵是高维的情况下，那么这个矩阵就是高维空间下的一个线性变换，这个线性变换可能没法通过图片来表示，但是可以想象，这个变换也同样有很多的变化方向，我们通过特征值分解得到的前N个特征向量，就对应了这个矩阵最主要的N个变化方向。我们利用这前N个变化方向，就可以近似这个矩阵变换。也就是之前说的：提取这个矩阵最重要的特征。

特征值分解举例

这里我们用一个简单的方阵来说明特征值分解的步骤。我们的方阵A定义为：

$A = \begin{pmatrix} -1&1&0\\ -4&3&0\\ 1&0&2 \end{pmatrix}$
首先，由方阵A的特征方程，求出特征值。
$|A-\lambda E|= \begin{vmatrix} -1-\lambda & 1&0 \\ -4 & 3-\lambda &0 \\ 1&0&2-\lambda \end{vmatrix}=(2-\lambda)\begin{vmatrix} -1-\lambda & 1 \\ -4 & 3-\lambda \end{vmatrix}=(2-\lambda)(\lambda-1)^2=0$
解方程得 $\lambda=2，1$ （重数为2）。
然后，把每个特征值 $\lambda$ 代入到线性方程组 $(A-\lambda E)x=0$ 里面，求出特征向量。
当 $\lambda=2$ 时，解线性方程组 $(A-2E)x=0$
$(A-2E)=\begin{pmatrix} -3&1&0\\ -4&1&0\\ 1&0&0 \end{pmatrix}->\begin{pmatrix} 1&0&0\\ 0&1&0\\ 0&0&0 \end{pmatrix}$ ，解得 $x_1=0,x_2=0$ 。特征向量为：
$p1=\begin{pmatrix} 0\\ 0\\ 1 \end{pmatrix}$

同理，当 $\lambda=1$ 时，解线性方程组 $(A-E)x=0$
$(A-E)=\begin{pmatrix} -2&1&0\\ -4&2&0\\ 1&0&1 \end{pmatrix}->\begin{pmatrix} 1&0&1\\ 0&1&2\\ 0&0&0 \end{pmatrix}$
$x_1+x_3=0,x_2+2x_3=0$ ，特征向量为： $p2=\begin{pmatrix} -1\\ -2\\ 1 \end{pmatrix}$
最后，方阵A的特征值分解为：

$A=Q\sum Q^{-1}=\begin{pmatrix} 0&-1&-1\\ 0&-2&-2\\ 1&1&1 \end{pmatrix}\begin{pmatrix} 2&0&0\\ 0&1&0\\ 0&0&1 \end{pmatrix}\begin{pmatrix} 0&-1&-1\\ 0&-2&-2\\ 1&1&1 \end{pmatrix}^{-1}$

奇异值分解

上面讲解的特征值分解在实际应用的时候有一个最致命的缺点，就是只能用于方阵，也即是n*n的矩阵，而我们实际应用中要分解的矩阵大多数都不是方阵，所以SVD出现了。

什么是奇异值分解

奇异值分解是一个能适用于任意矩阵的一种分解的方法，对于任意矩阵A总是存在一个奇异值分解：
$A=U\sum V^T$
假设A是一个m*n的矩阵，那么得到的U是一个m*m的方阵，U里面的正交向量被称为左奇异向量。 $\sum$ 是一个m*n的矩阵， $\sum$ 除了对角线其它元素都为0，对角线上的元素称为奇异值。 $V^T$ 是V的转置矩阵，是一个n*n的矩阵，它里面的正交向量被称为右奇异值向量。而且一般来讲，我们会将 $\sum$ 上的值按从大到小的顺序排列。
从上面的奇异值分解公式来看，我们是不知道如何拆分矩阵A的。但我们可以想办法把奇异值和特征值联系起来。首先，我们用矩阵A的转置乘以矩阵A，得到一个方阵，用这个方阵进行特征值分解，得到的特征值和特征向量满足下面的等式：
$(A^TA)v_i=\lambda_iv_i$
这里的 $v_i$ 就是我们要求的右奇异向量。
其次，我们将A和A的转置做矩阵的乘法，得到一个方阵，用这样的方阵进行特征分解，得到的特征和特征向量满足下面的等式：
$(AA^T) \mu _i = \lambda_i \mu_i$
这里的 $\mu_i$ 就是左奇异向量。

上面直接给出了 $\mu_i$ 是左奇异向量， $v_i$ 是右奇异向量，根据是什么呢？我们尝试以V矩阵举例来证明一下：
$A=U\sum V^T=>A^T=V\sum^TU^T=>A^TA=V\sum^TU^TU\sum V^T=V\sum^2V^T(*)$
上式证明种使用了 $U^TU=I,\sum^T\sum=\sum^2$ 。可以看出， $A^TA$ 的特征向量组成的矩阵就是我们SVD种的V矩阵，而 $AA^T$ 得特征矩阵就是我们SVD中的U矩阵。
所以，我们可以得到奇异值的两种求法：

第一种：
$A=U\sum V^T=>AV=U\sum V^TV=>AV=U\sum=>Av_i=\delta_iu_i=>\delta_i=\frac{A_vi}{u_i}$
第二种：
通过上面*式的证明，我们还可以看出，特征值矩阵等于奇异值矩阵的平方，也就是说特征值和奇异值满足如下关系：

$\delta_i=\sqrt{\lambda_i}$
其中， $\sigma_i$ 就是奇异值，奇异值 $\sigma_i$ 和特征值相似，在矩阵 $\sigma$ 种也是从大到小排列。
上面了我们将了奇异值分解的原理，那么我们就可以直接按照上面的公式来计算了吗？
考虑一下直接使用上面的公式来计算，我们可以计算一下算法的复杂度，我们将m*n的矩阵分解为 $m\times m$ 的矩阵U, $m\times n$ 的矩阵 $\sum$ 和 $n*n$ 的矩阵 $V^T$ ，这个矩阵的复杂度大概是O(N^3)级别的，所以在求奇异值的时候需要耗费非常多的时间，整个算法的开销就过大了。
幸运的是，当我们将分解矩阵 $\sum$ 种的奇异值按照从大到小的顺序排列之后，奇异值 $\delta_i$ 从大到小的顺序减小的特别快。在很多情况下，前10%甚至1%的奇异值的和就占了全部的奇异值之和的99%以上。也就是说，剩下的90%甚至99%的奇异值几乎没有什么作用。因此，我们可以用前面r个大的奇异值来近似描述矩阵，于是奇异值分解公式可以写成如下：

$A_{m*n}\approx U_{m*r}\sum_{r*r}V_{r*n}^T$

其中r是一个远远小于m和n的数，右边的三个矩阵的相乘结果将会得到一个接近A的矩阵。如果r的取值远远小于n，从计算机内存的角度来说，右边三个矩阵的存储内存要远远小于矩阵A的。所以在奇异值分解中r的取值很重要，就是在计算精度和时间空间之间做选择。

奇异值分解举例

定义矩阵A为：
$A=\begin{pmatrix} 0&1\\ 1&1\\ 1&0 \end{pmatrix}$
首先，我们先求出 $A^TA$ 和 $AA^T$ ：
$A^TA=\begin{pmatrix} 0&1&1\\ 1&1&0\\ \end{pmatrix}*\begin{pmatrix} 0&1\\ 1&1\\ 1&0 \end{pmatrix}=\begin{pmatrix} 2&1\\ 1&2\\ \end{pmatrix}$

$AA^T=\begin{pmatrix} 0&1\\ 1&1\\ 1&0 \end{pmatrix}*\begin{pmatrix} 0&1&1\\ 1&1&0\\ \end{pmatrix}=\begin{pmatrix} 1&1&0\\ 1&2&1\\ 0&1&1 \end{pmatrix}$
然后，求出 $A^TA$ 和 $AA^T$ 的特征值和特征向量：
$A^TA$ 的特征向量和特征值:
$\lambda_1=3, v_1=\begin{pmatrix} \frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}\\ \end{pmatrix}$
$\lambda_2=1, v_2=\begin{pmatrix} \frac{-1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}\\ \end{pmatrix}$
$AA^T$ 的特征向量和特征值:

$\lambda_1=3, u_1=\begin{pmatrix} \frac{1}{\sqrt{6}}\\ \frac{2}{\sqrt{6}}\\ \frac{1}{\sqrt{6}} \end{pmatrix}$
$\lambda_2=1, u_2=\begin{pmatrix} \frac{1}{\sqrt{2}}\\ 0\\ \frac{-1}{\sqrt{2}}\\ \end{pmatrix}$
$\lambda_3=0, u_3=\begin{pmatrix} \frac{1}{\sqrt{3}}\\ \frac{-1}{\sqrt{3}}\\ \frac{1}{\sqrt{3}}\\ \end{pmatrix}$
然后，我们再利用 $Av_i=\sigma_i\mu_i,i=1,2$ 求奇异值：
$\begin{pmatrix} 0&1\\ 1&1\\ 1&0\\ \end{pmatrix}*\begin{pmatrix} \frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}\\ \end{pmatrix}=\sigma_1\begin{pmatrix} \frac{1}{\sqrt{6}}\\ \frac{2}{\sqrt{6}}\\ \frac{1}{\sqrt{6}} \end{pmatrix}=>\sigma_1=\sqrt{3}$
$\begin{pmatrix} 0&1\\ 1&1\\ 1&0\\ \end{pmatrix}*\begin{pmatrix} \frac{-1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}\\ \end{pmatrix}=\sigma_2\begin{pmatrix} \frac{1}{\sqrt{2}}\\ 0\\ \frac{-1}{\sqrt{2}} \end{pmatrix}=>\sigma_2=1$

这一步也可以直接用 $\sigma_i=\sqrt{\lambda_i}$ 求出奇异值 $\sqrt{3}$ 和1。
最后，我们得到 $A$ 的奇异值分解为:
$A=U\sum V^T=\begin{pmatrix} \frac{1}{\sqrt{6}}&\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{3}}\\ \frac{2}{\sqrt{6}}&0&\frac{-1}{\sqrt{3}}\\ \frac{1}{\sqrt{6}}&\frac{-1}{\sqrt{2}}&\frac{-1}{\sqrt{3}}\\ \end{pmatrix}*\begin{pmatrix} \sqrt{3}&0\\ 0&1\\ 0&0 \end{pmatrix}*\begin{pmatrix} \frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}\\ \frac{-1}{\sqrt{2}}&\frac{1}{\sqrt{2}}\\ \end{pmatrix}$

容易发现，经过SVD分解后，可以用前r个非零奇异值对应的奇异向量表示矩阵A的主要特征，这样就把矩阵A进行了降维。

协方差和散度矩阵

样本均值：
$\hat{x}=\frac{1}{n}\sum_{i=1}^Nx_i$
样本方差：
$S^2=\frac{1}{n-1}\sum_{i=1}^n(x_i-\hat{x})^2$
样本X和样本Y的协方差：
$Cov(X,Y)=E[(X-E(X))(Y-E(Y))]=\frac{1}{n-1}\sum_{i=1}^n(x_i-\hat{x})(y_i-\hat{y})$
由上面的公式，我们可以得到以下结论：

(1) 方差的计算公式是针对一维特征，即针对同一特征不同样本的取值来进行计算得到；而协方差则必须要求至少满足二维特征；方差是协方差的特殊情况。

(2) 方差和协方差的除数是n-1,这是为了得到方差和协方差的无偏估计。
协方差为正时，说明X和Y是正相关关系；协方差为负时，说明X和Y是负相关关系；协方差为0时，说明X和Y是相互独立。Cov(X,X)就是X的方差。当样本是n维数据时，它们的协方差实际上是协方差矩阵(对称方阵)。例如，对于3维数据(x,y,z)，计算它的协方差就是：

在这里插入图片描述散度矩阵定义为：

在这里插入图片描述对于数据的散度矩阵为。协方差矩阵和散度矩阵关系密切，散度矩阵就是协方差矩阵乘以。因此它们的特征值和特征向量是一样的。这里值得注意的是散度矩阵就是上面讲的SVD分解的第一步，因此PCA和SVD也是有相关性的。

PCA算法

PCA即(Principal Component Analysis)主成分分析算法，是机器学习种应用得最广泛的数据降维算法。PCA的思想是将原始n维的数据映射到k维上(k<n)，这k维是全新的正交特征，也叫主成分。PCA的工作就是在原始的数据空间种顺序的找一组相互正交的坐标轴，新的坐标轴和数据本身是密切相关的。其中第一个坐标轴是原始数据中方差最大的方向，第二个坐标轴是和第一个坐标轴相交的坐标轴种最大的，以此内推，k个坐标轴是完全正交的。研究发现，大部分方差都包含在k个坐标轴种，后面的坐标轴所含的方差几乎维0。所以可忽略不计，以实现对数据的降维处理。

PCA算法实现

基于特征值分解协方差矩阵实现PCA算法

输入数据集 $X={x_1,x_2,x_3,...,x_n}$ ，需要降维到k维。
1）去均值，即将每一维特征减掉各自的平均值。
2）计算协方差矩阵 $\frac{1}{n}XX^T$ ，注：里除或不除样本数量n或n-1,其实对求出的特征向量没有影响。
3）用特征值分解方法求协方差矩阵 $\frac{1}{n}XX^T$ 的特征值与特征向量。
4）对特征值从大到小排序，选择其中最大的k个。然后将其对应的k个特征向量分别作为行向量组成特征向量矩阵P。
5）将数据转换到k个特征向量构建的新空间中，即Y=PX。

举个例子：
我们要将之前提到的矩阵 $A=\begin{pmatrix} -1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1 \end{pmatrix}$ 降到1维。
因为这个矩阵的每行已经是零均值，这里我们直接求协方差矩阵：
$C=\frac{1}{5}\begin{pmatrix} -1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1 \end{pmatrix}\begin{pmatrix} -1 & -2 \\ -1 & 0 \\ 0 & 0 \\ 2 & 1 \\ 0 & 1 \end{pmatrix}=\begin{pmatrix} \frac{6}{5} & \frac{4}{5} \\ \frac{4}{5} & \frac{6}{5} \end{pmatrix}$
然后求其特征值和特征向量，得到：
$\lambda_1=2,\lambda_2=2/5$
其对应的特征向量分别是：
$c_1\begin{pmatrix} 1 \\ 1 \end{pmatrix},c_2\begin{pmatrix} -1 \\ 1 \end{pmatrix}$
其中对应的特征向量分别是一个通解, $c1$ 和 $c2$ 可以取任意实数。那么标准化后的特征向量为：
$\begin{pmatrix} 1/\sqrt{2} \\ 1/\sqrt{2} \end{pmatrix},\begin{pmatrix} -1/\sqrt{2} \\ 1/\sqrt{2} \end{pmatrix}$
因此我们的矩阵P是：
$P=\begin{pmatrix} 1/\sqrt{2} & 1/\sqrt{2} \\ -1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}$
可以验证协方差矩阵C的对角化：
$PCP^\mathsf{T}=\begin{pmatrix} 1/\sqrt{2} & 1/\sqrt{2} \\ -1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}\begin{pmatrix} 6/5 & 4/5 \\ 4/5 & 6/5 \end{pmatrix}\begin{pmatrix} 1/\sqrt{2} & -1/\sqrt{2} \\ 1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}=\begin{pmatrix} 2 & 0 \\ 0 & 2/5 \end{pmatrix}$
最后我们用P的第一行乘以数据矩阵，就得到了降维后的表示：
$Y=\begin{pmatrix} 1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}\begin{pmatrix} -1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1 \end{pmatrix}=\begin{pmatrix} -3/\sqrt{2} & -1/\sqrt{2} & 0 & 3/\sqrt{2} & -1/\sqrt{2} \end{pmatrix}$
数据矩阵X降维投影结果为：

在这里插入图片描述

基于SVD实现PCA算法

输入数据集 $X={x_1,x_2,x_3,...,x_n}$ ，需要降维到k维。

去平均值，即每一位特征减去各自的平均值。
计算协方差矩阵。
通过SVD计算协方差矩阵的特征值与特征向量。
对特征值从大到小排序，选择其中最大的k个。然后将其对应的k个特征向量分别作为列向量组成特征向量矩阵。
将数据转换到k个特征向量构建的新空间中。
和利用特征值分解实现PCA算法类似，我们需要找到样本协方差矩阵 $XX^T$ 的最大k个特征向量，然后用这最大的k个特征向量组成的矩阵来做低维投影降维。可以看出，在这个过程中需要先求出协方差矩阵 $XX^T$ 当样本数多、样本特征数也多的时候，这个计算还是很大的。当我们用到SVD分解协方差矩阵的时候，SVD有两个好处：

有一些SVD的实现算法可以先不求出协方差矩阵 $XX^T$ 也能求出我们的右奇异矩阵V。也就是说，我们的PCA算法可以不用做特征分解而是通过SVD来完成，这个方法在样本量很大的时候很有效。实际上，scikit-learn的PCA算法的背后真正的实现就是用的SVD，而不是特征值分解。
注意到PCA仅仅使用了我们SVD的左奇异矩阵，没有使用到右奇异值矩阵，那么右奇异值矩阵有什么用呢？
假设我们的样本是 $m\times n$ 的矩阵X，如果我们通过SVD找到了矩阵 $X^TX$ 最大的k个特征向量组成的k*n的矩阵 $V^T$ ，则我们可以做如下处理：
$X'_{m*k}=X_{m*n}V_{n*k}^T$
可以得到一个 $m\times k$ 的矩阵X',这个矩阵和我们原来 $m\times n$ 的矩阵X相比，列数从n减到了k，可见对列数进行了压缩。也就是说，左奇异矩阵可以用于对行数的压缩；右奇异矩阵可以用于对列(即特征维度)的压缩。这就是我们用SVD分解协方差矩阵实现PCA可以得到两个方向的PCA降维(即行和列两个方向)。

代码实现

import numpy as np
# 零均值化，即中心化，是数据的预处理方法
def zero_centered(data):
    matrix_mean = np.mean(data, axis=0)
    return data - matrix_mean

def pca_eig(data, n):
    new_data = zero_centered(data)
    conv_mat = np.dot(new_data.T, new_data) #也可以用np.cov()方法
    eig_values, eig_vectors = np.linalg.eig(np.mat(conv_mat))
    # 求特征值和特征向量，特征向量是列向量
    value_indices = np.argsort(eig_values) #将特征值从小到大排序
    n_vectors = eig_vectors[:, value_indices[-1: -(n+1): -1]]
    # 最大的n个特征值对应的特征向量
    return new_data * n_vectors #返回低维特征空间的数据

def pca_svd(data, n):
    new_data = zero_centered(data)
    cov_mat = np.dot(new_data.T, new_data)
    U, s, V = np.linalg.svd(cov_mat) #将协方差矩阵奇异值分解
    pc = np.dot(new_data, U) #返回矩阵的第一个列向量即是降维后的结果
    return pc[:, 0]

def unit_test():
    data = np.array([[2.5, 2.4], [0.5, 0.7], [2.2, 2.9], [1.9, 2.2], [3.1, 3.0], 
                     [2.3, 2.7], [2, 1.6], [1, 1.1], [1.5, 1.6], [1.1, 0.9]])
    result_eig = pca_eig(data, 1)
    # 使用常规的特征值分解法，将二维数据降到一维
    print(result_eig)
    result_svd = pca_svd(data, 1)
    # 使用奇异值分解法将协方差矩阵分解，得到降维结果
    print(result_svd)
unit_test()

得到的降维结果如下：

在这里插入图片描述

可以看到两种方法的降维结果是一致的。

参考文章

https://mp.weixin.qq.com/s/Dv51K8JETakIKe5dPBAPVg
https://blog.csdn.net/program_developer/article/details/80632779
http://blog.codinglabs.org/articles/pca-tutorial.html

欢迎关注我的微信公众号GiantPadaCV，期待和你一起交流机器学习，深度学习，图像算法，优化技术，比赛及日常生活等。

图片.png

机器学习算法之PCA算法

前言

前置内容

特征值分解

特征值，特征向量

特征值与特征向量的几何意义

特征值分解

特征值分解举例

奇异值分解

什么是奇异值分解

奇异值分解举例

协方差和散度矩阵

PCA算法

PCA算法实现

基于特征值分解协方差矩阵实现PCA算法

基于SVD实现PCA算法

代码实现

参考文章

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读