PCA算法

发表于 | 分类于 | 阅读次数:

综述

PCA算法，即主成分分析法（Principal Components Analysis, PCA），是常用的降维方法。总的来说，降维有以下优点，这些优点通常也是我们选择降维的原因。降维的优点如下：

  1. 避免维数灾难。当数据维度增加时，要维持相同程度的密采样所需要的样本数量会成指数级增长。比如，在[0, 1]区间内每隔0.01单位采样，仅需要100个采样点。而在三维空间中，如果将0.01小方格为采样间隔，则需要$10^6$个采样点。当维度更高时，采样点数量将极其庞大。
  2. 去除噪声。采样数据往往会带有噪声。以x、s、n分别表示采样数据、原始信号、噪声信号。通常可记其关系如下：$x = s + n$。由于PCA降维是保留方差最大的方向，而相对而言噪声信号n的方差要小于真实信号s的方差，因此关于噪声信号的相关信息会被滤除。
  3. 去除冗余。数据的各个维度可能并不完全独立，通过降维通过保留主要的维度从而去除冗余。

PCA算法原理及推导

从数据压缩的角度来看，我们希望压缩后的数据能够尽可能的保留更多的信息。以下为推导过程

样本数据集{$x_i$, i=1, 2, 3, …, m}，其中$x_i$为n维的列向量, 且已经去均均值化。${w_i, i=1, 2, 3, …, n}$表示特征空间的一组正交基, 其中$w_i$为n维列向量。令$W=[w_1, w_2, …, W_d]$表示由正交基中的d个基构成的矩阵，则W为$n*d$的矩阵。此时，

$$z_i=W^Tx_i \tag {1}$$
$$\hat{x_i}=Wz_i \tag {2}$$
最小化信息损失，则有

$$min \sum_{i=1}^m {\left|\right|Wz_i-x_i\left|\right|}^2 \tag {3}$$
将上式分解开来，则有

$$min -\sum_{i=1}^m z_i^Tz_i \tag {4}$$
此处分解忽视了只有$x_i$所引入的常量。由内积与矩阵的迹的对应关系$z_i^Tz_i=tr(z_iz_i^T)$，即可得

$$min -\sum_{i=1}^m tr(z_iz_i^T) \tag {5}$$

$$min -\sum_{i=1}^m tr(W^T(x_ix_i^T)W) \tag {6}$$
由于W由正交基构成，所以有

$$W^TW=I$$
又$\sum_i^mx_ix_i^T=XX^T$，即有如下优化

$$min -tr(W^TXX^TW) \tag {7}$$
$$W^TW=I \tag {8}$$
由拉格朗日乘子法可得

$$XX^TW=\lambda W \tag {9}$$

从方差最大化的角度，即我们希望降维后的样本尽可能的分开，可以得到同样的结果。从以上的推导过程中可以很容易的得到PCA的算法实现过程。

  1. 对数据去均值化
  2. 求协方差矩阵$A=\sum_i^mx_ix_i^T=XX^T$
  3. 对协方差矩阵A进行奇异值分解，选取特征值较大的d个特征值所对应的特征向量，即可得到W
  4. 由W可求出降维后的数据$Z=W^TX$

核化PCA算法

以上的PCA算法为线性的，有时我们需要引入非线性关系进行降维，这时，核技巧就派上用场了。下面我们来介绍核化PCA降维。

首先引入非线性关系$\phi$，由(9)式可得

$$XX^TW=(\sum_{i=1}^{m}\phi(x_i)\phi(x_i)^T)W=\lambda W \tag {10}$$

则有，

$$
W=\frac{1}{\lambda}(\sum_{i=1}^{m}\phi(x_i)\phi(x_i)^T)W
=\sum_{i=1}^{m}\phi(x_i)\frac{\phi(x_i)^TW}{\lambda}
\tag {11}
$$

记$\alpha_i=\frac{1}{\lambda}\phi(x_i)^TW$，则$\alpha_i$为d维行向量，则有

$$W=\sum_{i=1}^{m}\phi(x_i)\alpha_i \tag {12}$$

核技巧往往是已知核化后的内积，而非知道核本身，即

$$\kappa(x_i, x_j)=\phi(x_i)^T\phi(x_j)=\phi(x_j)^T\phi(x_i) \tag {13}$$

将式(12), (13)带入(10)式即可得到

$$
\sum_{i=1}^m\phi(x_i)(\sum_{j=1}^m\kappa(x_i, x_j)\alpha_j-\lambda\alpha_i)=0
\tag {14}
$$

令$K_{ij}=\kappa(x_i, x_j)$，$A=[\alpha_i;…;\alpha_m]$，则A为m*d的矩阵，有

$$KA=\lambda A \tag {15}$$
将K进行奇异值分解即可得到${\alpha_i, i=1,2,…,m}$，下面来整理一下维度，$x_i$为n维列向量，则X为$n*m$的矩阵，$\alpha_i$为d维行向量，则A为$m*d$的矩阵，K为$m*m$的矩阵，则KA仍然为$m*d$的矩阵。
令$z_i=W^T\phi(x_i)$, 则$z_i$为$d*1$的向量。核化PCA将$n*1$的向量降维为$d*1$的向量。则

$$
z_j=W^T\phi(x_j)=\sum_{i=1}^{m} \alpha_i^T \phi(x_i)^T x_j=\sum_{i=1}^{m} \alpha_i^T \kappa(x_i, x_j) \tag {16}
$$
将其向量化即可得到

$$Z=A^TK \tag {17}$$

以上即为KPCA算法的原理及推导，其算法流程如下：

  1. 对数据去均值化
  2. 求协方差矩阵$K$，$K_{ij}=\kappa(x_i, x_j)$
  3. 对协方差矩阵K进行奇异值分解，选取特征值较大的d个特征值所对应的特征向量，即可得到A
  4. 由W可求出降维后的数据$Z=A^TK$

PCA及KPCA实现

根据以上的算法过程，佐以numpy的函数库， 可以很轻易地实现PCA算法，其Python实现如下：

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
def PCA(X, dim):
    '''
    PCA算法实现
    :param X:X.shape = (m, n), n表示特征维度，m表示样本数量
    :param dim:表示转换后的维度
    :return:W(n*d), Z(m*d) W为转换矩阵，Z为转换后的数据
    '''
    # 均值化
    X = np.array(X) - np.mean(X, 0)
    A = np.matmul(X.T, X)
    u, s, v = np.linalg.svd(A)
    W = v[:, 0:dim]
    Z = np.matmul(X, W)
    return W, Z


def KPCA(X, dim, kernel, sigma):
    '''
    KPCA算法实现
    :param X: X.shape = (m, n), n表示特征维度，m表示样本数量
    :param dim: 表示转换后的特征维度
    :param kernel: 核方法
    :param sigma: 此处使用高斯核，因此直接显式表示出来了
    :return: A(m*d), Z(m*d) A为转换矩阵，Z为转换后的数据
    '''
    X = np.array(X) - np.mean(X, 0)
    K = [[kernel(xi, xj, sigma) for xi in X] for xj in X]
    u, s, v = np.linalg.svd(K)
    A = v[:, 0:dim]
    Z = np.matmul(K, A)
    return A, Z

def gaussian(xi, xj,sigma):
    return 1./sigma/np.sqrt(2\*np.pi)\*np.exp(-np.sum(np.square(np.subtract(xi, xj)))/sigma\*\*2)

参考文献

主要参考文献为周志华老师的《机器学习》西瓜书