降维-PCA算法理论详解与实现

Posted by mudux on February 20, 2020

1. 介绍

 PCA(Principal Component Analysis)是一种常见的数据分析方式,常用于高维数据的降维,可用于提取数据的主要特征分量。
 PCA的数学推导可以从最大可分性和最近重构性两方面进行,前者的优化条件为划分后方差最大,后者的优化条件为点到划分超平面的距离最小,这里将从两个方面分别介绍。

2. 基于最大可分性

2.1 向量表示与基变换

 先复习下线性代数的基本知识

2.1.1 內积

 设有两个$n$维向量$A,B$,$A=(a_1,a_2,\cdots,a_n)^T,B=(b_1,b_2,\cdots,b_n)^T$,这两个向量的內积形式是这样:

 內积运算将两个向量映射为实数。接下来从几何角度来分析其物理意义,简单起见,假设$A,B$均为二维向量,则:

 其几何表示如下:

2.1

 可看出$A$与$B$的內积等于$A$到$B$的投影长度乘以$B$的模。
 如果$B$的模为1,即$|B|=1$,那么就相当于:

 也就是说,$A$与$B$的內积等于$A$向$B$所在直线投影的标量长度。

2.1.2 基

 在我们常说的坐标系中,向量$(3,2)$隐式引入了一个定义:以$x$轴和$y$轴正方向长度为1的向量为标准。向量$(3,2)$实际是说在$x$轴投影为3而$y$轴投影为2。
 所以对于向量$(3,2)$来说,如果我们想求它在$(1,0),(0,1)$这组基下的坐标,分别內积即可,內积完了还是$(3,2)$。
 所以可看出要准确描述向量,首先要确定一组基,然后给出在基所在的各个直线上的投影值就可以了。为了方便求坐标,希望这组基向量模长为1。因为向量的內积运算,当模长为1时,內积可以直接表示投影。然后还需要这组基是线性无关的。

2.1.3 基变换的矩阵表示

 先给出一个例子,对于向量$(3,2)$来说,在$(\frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}})$和$(-\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})$这组基下的是:

 左边矩阵的两行分别为两个正交基,乘以原向量,其结果刚好为新基的坐标。推广一下,如果有$m$个二维向量,只要将二维向量按列排成一个两行$m$列矩阵,然后用“基矩阵”乘以这个矩阵就可以得到所有这些向量在新基下的值。可以把它写成通用的表示形式:

 其中$p_i$是一个行向量,表示第$i$个基,$a_j$是一个列向量,表示第$j$个原始数据记录。

2.2 最大可分性

 上面讨论了选择不同的基可以对同样一组数据给出不同的表示,如果基的数量少于向量本身的维数,则可以达到降维的效果。降维肯定会导致信息的丢失。所以一个关键的问题是:如何选择基才是最优的。如果有一组$N$维向量,现在要将其降到$K$维($K$小于$N$),那么应该如何选择$K$个基才能最大程度的保留原有的信息?
 一种直观的看法是:希望投影后的投影值尽可能分散,因为如果重叠就会有样本消失。当然这个也可以从熵的角度进行理解,熵越大所含信息越多。

2.2.1 方差

 方差定义:

 为了方便处理,将变量的均值化为0,方差表示为:

 这时上面的问题被形式化表述为:寻找一组基,使得所有数据变换为这个基上的坐标表示后,方差值最大。

2.2.2 协方差

 用方差表示单个变量的分散程度,而协方差可以表示两个变量的相关性。注意这里的变量是一个序列,即一个特征。 协方差公式为:

 由预处理后均值为0,所以这里协方差可以表示为:

 为了让两个变量尽可能表示更多的原始信息,希望它们之间不存在线性相关性,因为相关性意味着两个变量不是完全独立,必然存在重复表示的信息,即协方差为0。为了让协方差为0,选择第二个基时只能在与第一个基正交的方向上进行选择。
 至此,得到了降维问题的优化目标:将一组$N$维向量降为$K$维,其目标是选择$K$个单位正交基,使得原始数据变换到这组基上后,各变量两两间协方差为0,而变量方差则尽可能大(在正交的约束下,取最大的$K$个方差)

2.2.3 协方差矩阵

 我们最终要达到的目标与变量内方差及变量间协方差有密切关系。因此希望能将两者统一表示,仔细观察发现,两者均可以表示为內积的形式,而內积又与矩阵相乘密切相关。
 假设只有$a$和$b$两个变量,将它们按行组成矩阵$X$:

然后

 矩阵对角线上是两个变量的方差,而其它元素是协方差。推广到一般情况:设有$m$个$n$维数据样本,将其排列成矩阵$X_{n\times m}$,设$C=\frac{1}{m}XX^T$,则$C$是一个对称矩阵,其对角线分别对应各个变量的方差,而第$i$行$j$列和第$j$行$i$列元素相同,表示第$i$和$j$两个变量的协方差

2.2.4 矩阵对角化

 根据优化条件,需要将除对角线外的其它元素化为0,并且在对角线上将元素按大小从上到下排列(变量方差尽可能大),这样就达到了优化目的。
 设原始数据矩阵$X$对应的协方差矩阵为$C$,而$P$是一组基按行组成的矩阵,设$Y=PX$,则$Y$为$X$对$P$做基变换后的数据。设$Y$的协方差矩阵为$D$,推导下$D$与$C$的关系:

 这样就很清楚了,要找的是能让原始协方差矩阵对角化的$P$。优化目标变成寻找一个矩阵$P$,满足$PCP^T$是一个对角矩阵,并且对角元素按从大到小依次排列,那么$P$的前$K$行就是要寻找的基,用$P$的前$K$行组成的矩阵满足乘以$X$就使得$X$从$N$维降到了$K$维并满足上述优化条件。
 至此,我们离PCA仅一步之遥,还需要完成对角化。
 由上文知道,协方差矩阵$C$是一个对称矩阵,在线性代数中实对称矩阵有一系列性质:

  • 实对称矩阵不同特征值对应的特征向量必然正交;
  • 实对称矩阵特征值$\lambda$的几何重数与代数重数相等,设代数重数为$r$,则必然存在$r$个线性无关的特征向量对应$\lambda$,因此可以将这$r$个特征向量单位正交化。

 由上面两条可知,一个$n$行$n$列的实对称矩阵一定可以找到$n$个单位正交特征向量,设这$n$个特征向量为$(e_1,e_2,\cdots,e_n)$,将其按列组成矩阵: $E=(e_1,e_2,\cdots,e_n)$。
 则协方差矩阵$C$有如下结论:

 其中$\Lambda$为对角矩阵,其对角元素为各特征向量对应的特征值。
 到这里,发现我们已经找到了需要的矩阵$P:P=E^T$。
 $P$是协方差矩阵的特征向量单位化后按行排列出的矩阵,其中每一行都是$C$的一个特征向量。如果$P$按照$\Lambda$中特征值的从大到小将特征向量从上到下排列,则用$P$的前$K$行组成矩阵乘以原始数据矩阵$X$,就得到了我们需要的降维后的数据矩阵$Y$。

3. 最近重构性

 假设$m$个$n$维数据$\left(x^{(1)},x^{(2)},\cdots,x^{(m)}\right)$都已经进行了中心化,即$\sum\limits_{i=1}^mx^{(i)}=0$。经过投影变换后得到的新坐标系为${w_1,w_2,\cdots,w_n}$,其中$w$是标准正交基,即
 如果我们将数据从$n$维降到$n^{‘}$维,即丢弃新坐标系中的部分坐标,则新的坐标系为${w_1,w_2,\cdots,w_{n^{‘}}}$,样本点$x^{i}$在$n^{‘}$维坐标系中的投影为:$z^{(i)}=\left(z_1^{(i)},z_2^{(i)},\cdots,z_{n^{‘}}^{(i)}\right)$。其中,$z_j^{(i)}=w_j^Tx^{(i)}$是$x^{(i)}$在低维坐标系里第$j$维的坐标。
 用$z^{(i)}$来恢复原始数据$x^{(i)}$,则得到的恢复数据${\overline x}^{(i)}=\sum\limits_{j=1}^{n^{‘}} z_j^{(i)}w_j=Wz^{(i)}$,其中,$W$为标准正交基组成的矩阵。
现在考虑整个数据集,希望所有样本到这个超平面的距离足够近,即最小化下式:

 将这个式子进行整理,可以得到:

注意这里$X$每个样本按列排列,注意$\sum\limits_{i=1}^m x^{(i)T}x^{(i)}$是一个常量。最小化上式等价于:

 利用拉格朗日函数可以得到:

 对$W$求导有$-XX^TW + \lambda W=0$,整理下即为:

 这样可清楚看出,$W$为$XX^T$的$n^{‘}$个特征向量组成的矩阵,而$\lambda$为$XX^T$的若干特征值组成的矩阵,特征值在注对角线上,其余位置为0。为了使距离最近,选特征值最大的$n^{‘}$个特征值对应的特征向量组成投影矩阵

4. 最近重构性另一个解释

 假设$m$个$n$维数据$\left(x^{(1)},x^{(2)},\cdots,x^{(m)}\right)$都已经进行了中心化,即$\sum\limits_{i=1}^mx^{(i)}=0$。经过投影变换后得到的新坐标系为${w_1,w_2,\cdots,w_n}$,其中$w$是标准正交基,即
 如果我们将数据从$n$维降到$n^{‘}$维,即丢弃新坐标系中的部分坐标,则新的坐标系为${w_1,w_2,\cdots,w_{n^{‘}}}$,样本点$x^{i}$在$n^{‘}$维坐标系中的投影为:$z^{(i)}=\left(z_1^{(i)},z_2^{(i)},\cdots,z_{n^{‘}}^{(i)}\right)$。其中,$z_j^{(i)}=w_j^Tx^{(i)}$是$x^{(i)}$在低维坐标系里第$j$维的坐标。
 对于任意一个样本$x^{(i)}$,在新的坐标系中的投影为$W^Tx^{(i)}$,在新坐标系中的投影方差为$x^{(i)T}WW^Tx^{(i)}$,要使所有的样本的投影方差和最大,也就是最大化$\sum\limits_{i=1}^m W^T x^{(i)} x^{(i)T}$的迹,即:

 和第三节基于最近投影的优化目标完全一样,只是符号问题。
 利用拉格朗日函数可以得到:

 对$W$求导有$XX^TW + \lambda W=0$,整理下即为:

5. PCA算法流程

 输入:$n$维样本集$D=\left(x^{(1)},x^{(2)},\cdots,x^{(m)}\right)$,要降维到$K$维;
 输出: 降维后的样本集$D^{‘}$

  1. 对所有样本进行中心化:$x^{(i)} = x^{(i)} - \frac{1}{m}\sum\limits_{i=1}^m x^{(j)}$;
  2. 计算样本的协方差矩阵$XX^T$;
  3. 对矩阵$XX^T$进行特征值分解;
  4. 取出最大的$K$个特征值对应的特征向量$(w_1,w_2,\cdots,w_K)$,将所有的特征向量标准化后,组成特征向量矩阵$W$;
  5. 输出样本集$D^{‘}=W^TD$。

 有时候,我们不指定降维后的$K$的值,而是换种方式,指定一个降维到的主成分比重阈值$t$。这个阈值$t$在(0,1]之间。假如我们的$n$个特征值为$\lambda_1 \ge \lambda_2 \ge \cdots \ge \lambda_n$,则$K$可以通过下式得到:

6. 核主成分分析PCA

 在上面的PCA算法中,我们假设存在一个线性的超平面,可以让我们对数据进行投影。但是有些时候数据不是线性的,不能直接进行PCA降维。这里就需要用到和支持向量机一样的核函数的思想,先把数据集从$n$维映射到线性可分的高维$N>n$,然后再从$N$维降维到一个低维度$n^{‘}$,这里的维度之间满足$n^{‘}<n<N$。
 使用了核函数的主成分分析一般称之为核主成分分析(Kernelized PCA,以下简称KPCA)。假设高维空间的数据是由$n$维空间的数据通过映射$\phi$产生。
 则对于$n$维空间的特征分解:

 映射为:

 通过在高维空间进行协方差矩阵的特征值分解,然后用和PCA一样的方法进行降维。一般来说,映射$\phi$不用显式的计算,而是在需要计算的视乎通过核函数完成。

6. PCA代码实现

6.1 python实现

 PCA.py文件

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
import numpy as np

class PCA():
    def __init__(self, n_components=2):
        self.n_components = n_components

    def fit(self, X):
        if not isinstance(X, (np.ndarray, np.generic)):
            X = np.array(X)
        X = X.T # instance of input X is represented by each row, transpose it
        meanVals = np.mean(X, axis=1, keepdims=True)
        self.mean_ = meanVals
        meanRemoved = X - meanVals
        covMat = np.cov(meanRemoved, rowvar=True, bias=False)
        eigVals, eigVec = np.linalg.eig(covMat) # each column of eigVec is corresponding eigen vector
        eigValInd = np.argsort(eigVals)
        eigValInd = eigValInd[:-(self.n_components+1):-1]
        selectedEigVec = eigVec[:, eigValInd]
        self.eigVec_ = selectedEigVec # each column is eigen vector


    def tranform(self, X):
        self.check_is_fitted("mean_")
        X = X.T
        meanRemoved = X - self.mean_
        reducedX = np.dot(self.eigVec_.T, meanRemoved)
        # reconX = np.dot(self.eigVec_, reducedX)
        return reducedX.T   


    def fit_transform(self, X):
        self.fit(X)
        X = X.T
        meanRemoved = X - self.mean_
        reducedX = np.dot(self.eigVec_.T, meanRemoved)
        # reconX = np.dot(self.eigVec_, reducedX)
        return reducedX.T

    def inverse_transform(self, reducedX):
        self.check_is_fitted("mean_")
        reducedX = reducedX.T
        reconX = np.dot(self.eigVec_, reducedX)
        reconX = reconX + self.mean_
        return reconX.T
        

    def check_is_fitted(self, attributes, msg=None, all_or_any=all):
        if msg is None:
           msg = ("This %(name)s instance is not fitted yet. Call 'fit' with "
                   "appropriate arguments before using this method.")
        if not hasattr(self, 'fit'):
            raise TypeError("This is not an estimator instance.")
        if not isinstance(attributes, (list, tuple)):
            attributes = [attributes]
        if not all_or_any([hasattr(self, attr) for attr in attributes]):
            raise TypeError(msg % {'name': type(self).__name__})

测试一下

1
2
3
4
5
6
7
8
9
import numpy as np
from PCA import PCA

X = np.array([[-1, -1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]])
pca = PCA(n_components=2)
reducedX = pca.fit_transform(X)
print("my pca reducedX", reducedX)
reconX = pca.inverse_transform(reducedX)
print("my pca reconX", reconX)

 输出为:

6.1

6.2 sklearn库

1
2
3
4
5
6
7
8
9
import numpy as np
from sklearn.decomposition import PCA

X = np.array([[-1, -1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]])
pca = PCA(n_components=2)
reducedX = pca.fit_transform(X)
print("sklearn pca reducedX\n", reducedX)
reconX = pca.inverse_transform(reducedX)
print("sklearn pca reconX\n", reconX)

 输出为:

6.2

可以发现两个输出只有符号不一样,这和特征向量的符号有关,不影响投影,重构后结果是一样的。

7. PCA总结

 作为一个非监督学习的降维方法,它只需要特征值分解,就可以对数据进行压缩,去噪。因此在实际场景应用很广泛。为了克服PCA的一些缺点,出现了很多PCA的变种,比如为解决非线性降维的KPCA,还有解决内存限制的增量PCA方法Incremental PCA,以及解决稀疏数据降维的PCA方法Sparse PCA等。
 PCA算法的主要优点有:

  • 缓解维度灾难:PCA 算法通过舍去一部分信息之后能使得样本的采样密度增大(因为维数降低了),这是缓解维度灾难的重要手段;
  • 降噪:当数据受到噪声影响时,最小特征值对应的特征向量往往与噪声有关,将它们舍弃能在一定程度上起到降噪的效果; 
  • 特征独立:各主成分之间正交,可消除原始数据成分间的相互影响的因素;
  • 计算简单:主要运算是特征值分解,易于实现。

 PCA算法的主要缺点有:

  • 主成分各个特征维度的含义具有一定的模糊性,不如原始样本特征的解释性强。
  • 方差小的非主成分也可能含有对样本差异的重要信息,因降维丢弃可能对后续数据处理有影响。

8. Reference

1 降维——PCA-阿泽
2 主成分分析(PCA)原理总结-刘建平
3 机器学习第十章-周志华
4 机器学习实战第十三章