主成分分析

主成分分析

概述

换一个角度去看世界，世界会变得大不一样。
对于 nn 维的数据而言，在原先的 nn 维坐标系中的基底下，在各个轴的方向上，可能看不出来有什么差别。但是如果换另一个坐标系去看，就会发现它们是沿着某一个或者某几个轴分布的。更换合适的坐标系以后我们会发现，沿着某些轴时，它们的分布范围很大；沿着另一些轴时，它们的分布范围很小。在某些极端情况下，我们会发现数据严格分布在某个平面或者超平面上，以下图的三维数据为例：

20150304195708775.png

当我们在原始的坐标系下看的时候，它们是没有规律的。然而当我们把图中 的PC1⃗\vec{PC_1}，PC2⃗\vec{PC_2} 以及另外一个垂直于二者但是没有画出来的向量放在一起组成新的坐标系的时候，数据的特征已经呼之欲出了：它们完全分布在由相互正交的 PC1⃗\vec{PC_1} 和 PC2⃗\vec{PC_2} 确定的平面上。在 PC1⃗\vec{PC_1} 方向，数据分布范围较大，我们换一种说法，叫做方差较大，对应地，在 PC2⃗\vec{PC_2} 方向，数据的方差较小。这两个方向上，我们都可以看出数据点之间的差别，只是大小有别。但是在垂直于两个向量的方向上，数据没有方差，也就是说，我们看不出它们的差别。 在这个维度上的数据完全没有作用。

**主成分分析（PCA）的作用，就是找出一个能让数据在某些轴方向上的差别尽可能大的坐标系来。**至于这个坐标系的维度，取决于新选取的坐标系上能够体现出明显的数据差异的轴的数量，也可以说是维度的数量。当我们选择了合适的维度，运气也足够好（数据集的维度里面正好有“没什么卵用”的维度），就可以看出，坐标系重选之后，例如对于一个新的五维坐标系，它的正交基底 e⃗=e1⃗+...+e5⃗\vec{e}=\vec{e_1}+...+\vec{e_5} ，可能在 e1⃗\vec{e_1} 方向上数据方差很大，在 e5⃗\vec{e_5} 上数据方差却很小，这个时候，我们就要考虑放弃 e5⃗\vec{e_5}，因此数据实际上变成了 4 维。

作用

主成分分析的主要目的是，尽量保留有用的数据，并降低数据的维数。

主成分分析的主要应用范围

图像压缩

图像也是数据集。对于一个 200*400 的灰度图像来说，你可以理解成它由 400 个 200 维的向量组成。它组成了一个 200*400 的矩阵。这个时候就可以对图像进行主成分分析，降维，实现压缩的效果。通常，主成分分析对图像进行处理可以在几乎还原全图的情况下达到 90% 以上的压缩率。

信号处理

一般认为信号的方差是很大的（参考正弦函数），噪声的方差是很小的。那么主成分分析就可以帮助我们提取信号，去掉噪声。

机器学习上的作用

主成分分析用在机器学习上，主要还是用来给数据降维，避免高维灾难。
高维灾难是指，当数据的维数变高的时候，数据会变得及其稀疏。
比如说我们要用 K-Means 聚类算法。高维的时候，利用欧氏距离来度量两点间的距离就不像二维情况那样合适了。我们需要降维。这个时候就可以用主成分分析来降低维数，降低聚类的难度。

推导

那么如何找到这个合适的坐标系呢？
首先我们找原点。事实上，将原点放在样本中心是坠吼的。所以首先要找到样本数据的中心点，然后进行坐标变换，把原点移到这个点去，实际应用中把每个数据点进行一次变换即可。
由前面的叙述可以得知，样本方差越大越好。那么我们现在要选出样本方差最大的方向。
因为数据已经中心化，那么所有方向投影的均值都为0。
在任意一个方向 u\boldsymbol{u} 上，第 ii 个数据 x(i)\boldsymbol{x}^{(i)} 的投影长度为 x(i)⋅u\boldsymbol{x}^{(i)}\cdot \boldsymbol{u}。
注意u\boldsymbol{u}是单位向量。
201104182110524172.png

故样本方差可以表示为：

1m∑i=1m(x(i)u)2=1muTx(i)x(i)Tu=uT(1m∑i=1mx(i)x(i)T)u(1)\begin{aligned} \frac{1}{m}\sum_{i=1}^{m}(\boldsymbol{x}^{(i)}\boldsymbol{u})^2 &=\frac{1}{m}\boldsymbol{u}^T\boldsymbol{x}^{(i)}{\boldsymbol{x}^{(i)}}^T\boldsymbol{u} \\ &=\boldsymbol{u}^T \biggl(\frac{1}{m} \sum_{i=1}^{m}\boldsymbol{x}^{(i)}{\boldsymbol{x}^{(i)}}^T \biggr) \boldsymbol{u} \end{aligned} \tag{1}

用λ\lambda 表示样本方差，则问题转化为求 max⁡ λ=1m∑i=1m(x(i)u)2\max\ \lambda = \frac{1}{m}\sum_{i=1}^{m}(\boldsymbol{x}^{(i)}\boldsymbol{u})^2 的最优化问题。

令Σ\Sigma等于括号里的部分,因为 ∀i,E(x(i))=E(x(i)T)=0\forall i, E(\boldsymbol{x}^{(i)})=E({\boldsymbol{x}^{(i)}}^T)=0,所以有

Σ=1m∑i=1mx(i)x(i)T=1m∑i=1m[x(i)−E(x(i))][x(i)T−E(x(i)T)](2)\begin{aligned} \Sigma &= \frac{1}{m} \sum_{i=1}^{m}\boldsymbol{x}^{(i)}{\boldsymbol{x}^{(i)}}^T \\ &=\frac{1}{m} \sum_{i=1}^{m} \biggl[\boldsymbol{x}^{(i)}-E(\boldsymbol{x}^{(i)})\biggr] \biggl[{\boldsymbol{x}^{(i)}}^T-E({\boldsymbol{x}^{(i)}}^T)\biggr]\\ \end{aligned} \tag{2}

显而易见， Σ\Sigma 就是样本的协方差矩阵。（一般协方差矩阵前面是1m−1\frac{1}{m-1}，但是这不重要，具体请看下去。）
现在公式(1)可以化为 λ=uTΣu\lambda=\boldsymbol{u}^T \Sigma \boldsymbol{u}。
因为uu是单位向量，所以uTu=1\boldsymbol{u}^T \boldsymbol{u} = 1。两边左乘 u\boldsymbol{u}，得 Σu=λu\Sigma \boldsymbol{u} = \lambda \boldsymbol{u}。
因此我们得到约束条件：λ\lambda 就是 Σ\Sigma（样本协方差矩阵）的特征值，u\boldsymbol{u} 是特征向量。
这大大缩小了 λ\lambda 的解空间。
更重要的是，**协方差矩阵是半正定实对称矩阵，这意味着它所有的特征向量一定是彼此正交的。**这满足了坐标系的基底所需要的“彼此正交”的条件。
因此，我们可以得出结论：协方差矩阵 Σ\Sigma 最大的特征值对应的特征向量，就是样本方差最大的方向，即新坐标系里对数据影响最大的方向；第二大特征值对应的特征向量，是次主要的方向，以此类推。

现在我们就知道为什么前面说1m\frac{1}{m}和1m−1\frac{1}{m-1}的差别不重要了。因为更换这个数字也只是对坐标系进行线性变换而已，不影响方向。所以我们直接拿协方差矩阵进行运算也是可以的。

至此，主成分分析的方法已经确定下来了。我们将它归纳在下一节。

算法

  1. 将数据中心化；
  2. 求出样本矩阵 XX 的协方差矩阵 Σ\Sigma；
  3. 求 Σ\Sigma 的特征值矩阵，按照从大到小排序，选出前 kk 个特征值 (k<nk<n) 对应的特征向量，得到向量矩阵U=[u1u2⋯uk] U= \begin{bmatrix} \boldsymbol{u}_1 & \boldsymbol{u}_2 & \cdots & \boldsymbol{u}_k \end{bmatrix}，其中 ui\boldsymbol{u}_i 是列向量。
  4. 新的样本 Y=UTXY = U^TX。（UTU^T 是 k×nk\times n 矩阵， XX 是 n×mn\times m 矩阵，所以 YY 是 k×mk\times m 矩阵，与 XX 相比，由 nn 维降到了 kk 维。）
  5. 通过选取最大的 kk 个特征值对应的特征向量，我们就将方差较小的特征（比如噪声）抛弃掉了。

数据还原

原理

新样本 YY 所在的坐标系是新选取的坐标系。在一些情况下（比如图像压缩），我们还要将这个数据还原到原来的空间。

由于 UU 是正交矩阵，所以 UT=U−1U^T = U^{-1}，即UUT=UTU=EUU^T=U^TU=E（此处的作用是，利用正交矩阵性质，用转置代替求逆，从而节省计算量），所以有

UY=UUTX=EX=X⇒X^=UYUY=UU^TX=EX=X \Rightarrow\\ \hat{X} = UY

因此数据还原的过程就是：将 YY 左乘 UU。

用于图片压缩时的压缩比计算

原图像的数字数量：m×nm\times n；
压缩后的数字数量：k×mk \times m；
需要一并传输的向量矩阵UU的数字数量：k×nk\times n；
压缩比为 k(n+m)mn\frac{k(n+m)}{mn}。

实际效果

此处给出了误差随着 kk 的增加而变化的趋势，以及当 kk 分别为8, 32, 64, 128, 256, 512时，算法对于一个 512 X 512 的图片的恢复效果。
误差-维数曲线.png

效果对比.png

代码实现

主成分分析的算法主要还是靠矩阵运算，所以利用 numpy 库可以很简单地把这个程序写出来。
下面的代码实现了主成分分析类的封装。
踩坑记录： python 的 cov 函数把一列当做一个向量， 特征值函数eigh() 的返回值亦以一列为一个特征向量。如果错把每一行当成一个特征向量，还原图片数据时可以观察到规则的数条白色噪声条纹。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
import numpy as np


class PCA(object):
    """
    PCA 训练对象
    """
    def __init__(self, k):
        """
        :param k: 要压缩到的维数
        """
        self.k = k
        self.mean = 0

    def fit(self, x):
        """
        :param x: 数据集 Numpy array, 每一行看做一条数据
        :return y: TransformedPCA对象，包含 过后的数据，每一行看做一条数据（一个向量）,数据默认为 float64
        """
        if self.k > x.shape[1]:
            raise ValueError('要压缩的维数大于原数据维数({k}>{x})'.format(k=self.k, x=x.shape[1]))

        k = self.k
        # 中心化
        self.mean = np.array([np.mean(i) for i in x.T])
        mx = np.float64(x) - self.mean
        # 计算协方差矩阵及其特征向量
        cov_x = np.cov(mx.T)
        D, V = np.linalg.eigh(cov_x)

        # 注意此处的特征向量是按列排序的
        u = V[:, -k:]  # 特征值的排序默认从小到大，故特征向量取后 k 个
        y = np.dot(u.T, mx.T).T
        return PCA.__TransformedPCA(self.mean, u, y)

    def __repr__(self):
        return '<PCA Fit Object {0}>'.format(id(self))

    # 子类：转换后的对象
    class __TransformedPCA (object):
        """
        numpy array 转换为PCA以后的对象

        可调用的属性：
            y: PCA运算之后的 numpy array, 每一行被看做一条数据（一个向量）
            u: 前 k 个特征向量，每一行被看做一个向量
        """
        def __init__(self, mean, u, y):
            self.u = u
            self.__mean = mean
            self.y = y

        def inverse_pca(self):
            """
            将PCA过后的数据还原到原来的线性空间
            :return: 原线性空间下的 numpy array，数据默认为 float64
            """
            # n = self.u.shape[1]
            # k = self.u.shape[0]
            # u = np.row_stack((self.u, np.zeros((n - k, n))))
            # y = np.column_stack((self.y, np.zeros((n, n - k))))

            y1 = np.dot(self.u, self.y.T).T
            # 坐标还原
            y1 += self.__mean
            return y1

        def __repr__(self):
            return '<Transformed PCA Object {0}>'.format(id(self))

应用示例:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
from pca import PCA
import numpy as np
from PIL import Image

k = 256
# 灰度图像
origin_lena = np.float64(Image.open('lena_grey.png'))
# k 决定压缩至多少维
pca = PCA(k)
# pca.fit()函数返回 PCA 对象，其 y 属性是变换后的结果，inverse_pca()方法返回还原后的numpy array
solve = pca.fit(origin_lena).inverse_pca()
# 防止溢出
solve[solve > 255] = 0
solve[solve < 0] = 0
Image.fromarray(np.uint8(solve)).save('test.png')

参考资料及扩展阅读

#旧站迁移文章
主成分分析
http://blog.yotubird.club/posts/2017/7276e58.html
作者
nqr
发布于
2017年10月27日
许可协议