主成分分析（PCA）原理总结

1. PCA的思想

PCA顾名思义，就是找出数据里最主要的方面，用数据里最主要的方面来代替原始数据。具体的，假如我们的数据集是n维的，共有 $m$ 个数据 $(x^{(1)},x^{(2)},\dots,x^{(m)})$ 。我们希望将这 $m$ 个数据的维度从 $n$ 维降到 $n'$ 维，希望这 $m$ 个 $n'$ 维的数据集尽可能的代表原始数据集。我们知道数据从 $n$ 维降到 $n'$ 维肯定会有损失，但是我们希望损失尽可能的小。那么如何让这 $n'$ 维的数据尽可能表示原来的数据呢？

$\qquad$ 我们先看看最简单的情况，也就是 $n=2，n'=1$ ,也就是将数据从二维降维到一维。数据如下图。我们希望找到某一个维度方向，它可以代表这两个维度的数据。图中列了两个向量方向， $u_1$ 和 $u_2$ ，那么哪个向量可以更好的代表原始数据集呢？从直观上也可以看出， $u_1$ 比 $u_2$ 好。
主成分分析（PCA）原理总结
$\qquad$ 为什么 $u_1$ 比 $u_2$ 好?可以有两种解释，第一种解释是样本点到这个直线的距离足够近，第二种解释是样本点在这个直线上的投影能尽可能的分开。

$\qquad$ 假如我们把 $n'$ 从1维推广到任意维，则我们的希望降维的标准为：样本点到这个超平面的距离足够近,或者说样本点在这个超平面上的投影能尽可能的分开。

$\qquad$ 基于上面的两种标准，我们可以得到PCA的两种等价推导。

2. PCA的推导:基于最小投影距离

$\qquad$ 我们首先看第一种解释的推导，即样本点到这个超平面的距离足够近。

$\qquad$ 假设m个n维数据 $(x^{(1)},x^{(2)},\dots,x^{(m)})$ 都已经进行了中心化，即 $\sum\limits_{i=1}^mx^{(i)} = 0$ ，经过投影变换后得到的新坐标系为 $(w_1,w_2,\dots,w_n)$ ,其中 $w$ 是正交基，即 $||w||^2=1,w_i^Tw_j=0$

$\qquad$ 如果我们将数据从 $n$ 维降到 $n'$ 维，即丢弃新坐标系中的部分坐标，则新的坐标系为 $(w_1,w_2,\dots,w_{n'})$ ，样本点 $x^{(i)}$ 在 $n'$ 维坐标系中的投影为： $z^{(i)} = (z_1^{(i)},z_2^{(2)},\dots,z_{n'}^{(i)})^T$ ，是一个列向量。其中， $z_j^{(i)}=w_j^Tx_j^{(i)}$ 是 $x^{(i)}$ 在低维坐标系里第j维的坐标， $w_j$ 为单位正交列向量。

$\qquad$ 如果我们用 $z^{(i)}$ 来恢复原始数据 $x^{(i)}$ ,则得到的恢复数据 $\overline{x}^{(i)} = \sum\limits_{j=1}^{n'}z_j^{(i)}w_j = Wz^{(i)}$ ，其中， $W$ 为标准正交基组成的矩阵,维度为 $(n,n')$ 。

$\qquad$ 现在我们考虑整个样本集，我们希望所有的样本到这个超平面的距离足够近，即最小化下式：
$\sum_{i=1}^m ||\overline{x}^{(i)} - x^{(i)}||^2$
$\qquad$ 将这个式子进行整理，可以得到:
$\begin{aligned} \sum_{i=1}^m ||\overline{x}^{(i)} - x^{(i)}||^2 &= \sum_{i=1}^m ||Wz^{(i)} - x^{(i)}||^2 \\ &= \sum_{i=1}^m [(Wz^{(i)})^T(Wz^{(i)}) - 2(Wz^{(i)})^Tx^{(i)} + x^{(i)T}x^{(i)}] \\ &=\sum_{i=1}^m [z^{(i)T}W^TWz^{(i)} - 2z^{(i)T}W^Tx^{(i)} + x^{(i)T}x^{(i)}] \\ &=\sum_{i=1}^m [z^{(i)T}z^{(i)} - 2z^{(i)T}z^{(i)} + x^{(i)T}x^{(i)}] \\ &=- \sum_{i=1}^m z^{(i)T}z^{(i)} + \sum_{i=1}^m x^{(i)T}x^{(i)} \\ &= - \sum_{i=1}^m tr(z^{(i)}z^{(i)T}) + C \\ &=-tr(W^Tx^{(i)}x^{(i)T}W) + C \end{aligned}$

$\qquad$ 其中使用到推导公式：

$\qquad$ (1） $\overline{x}^{(i)} = Wz^{(i)}$ , $z^{(i)} = W^Tx^{(i)}$ , (2) 平方和展开公式, (3）矩阵转置公式 $(AB)^T = B^TA^T$ , (4) $W^TW = I$ ,(5) 由矩阵的迹 $tr(AB) = tr(BA)$ ， $tr(z^{(i)} z^{(i)T}) = tr(z^{(i)T} z^{(i)})=z^{(i)T} z^{(i)}$

$\qquad$ 注意到 $\sum\limits_{i=1}^m x^{(i)}x^{(i)T}$ 是数据集的协方差矩阵，W的每一个向量 $w_j$ 是标准正交基。而 $\sum\limits_{i=1}^m x^{(i)Tx^{(i)}}$ 是一个常量C,最小化上式等价于：
$\underbrace{arg\ min}_W -tr(W^TXX^TW),\ \ \ s.t. \ W^TW =I$

$\qquad$ 这个最小化不难，直接观察也可以发现最小值对应的W由协方差矩阵 $XX^T$ 最大的 $n'$ 个特征值对应的特征向量组成。当然用数学推导也很容易。利用拉格朗日函数可以得到。
$J(W) = -tr(W^TXX^TW + \lambda(W^TW-I))$
$\qquad$ 对W求导有 $-XX^TW + \lambda W = 0$ 整理下即为：
$XX^TW = \lambda W$

$\qquad$ 这样可以更清楚的看出， $W$ 为 $XX^T$ 的 $n'$ 个特征向量组成的矩阵，而 $\lambda$ 为 $XX^T$ 的若干特征值组成的矩阵，特征值在主对角线上，其余位置为0。当我们将数据集从n维降到 $n'$ 维时，需要找到最大的 $n'$ 个特征值对应的特征向量。这 $n'$ 个特征向量组成的矩阵W即为我们需要的矩阵。对于原始数据集，我们只需要用 $z^{(i)} = W^Tx^{(i)}$ ,就可以把原始数据集降维到最小投影距离的 $n'$ 维数据集。

$\qquad$ 如果你熟悉谱聚类的优化过程，就会发现和PCA的非常类似，只不过谱聚类是求前k个最小的特征值对应的特征向量，而PCA是求前k个最大的特征值对应的特征向量。

3. PCA的推导:基于最大投影方差

$\qquad$ 现在我们再来看看基于最大投影方差的推导。

$\qquad$ 对于任意一个样本 $x^{(i)}$ ，在新的坐标系中的投影为 $W^Tx^{(i)}$ ，在新坐标系中的投影方差为 $(W^Tx^{(i)})^T(W^Tx^{(i)})$ ,要使所有的样本的投影方差和最大，也就是最大化 $\sum\limits_{i=1}^m (W^Tx^{(i)})^T(W^Tx^{(i)})$ ，由协方差矩阵的迹可知:
$\sum\limits_{i=1}^m (W^Tx^{(i)})^T(W^Tx^{(i)}) = tr((W^TX)^T(W^TX)) =tr(W^TXX^TW)$
即：
$\underbrace{arg\ max}_W \ \ tr(W^TXX^TW) \ \ \ s.t. \ \ \ W^TW = I$
$\qquad$ 观察第二节的基于最小投影距离的优化目标，可以发现完全一样，只是一个是加负号的最小化，一个是最大化。

$\qquad$ 利用拉格朗日函数可以得到:
$J(W) = tr(W^TXX^TW +\lambda(W^TW-I))$
$\qquad$ 对 $W$ 求导有 $XX^TW+\lambda W = 0$ ，整理下即为：
$XX^TW = -\lambda W$
$\qquad$ 和上面一样可以看出， $W$ 为 $XX^T$ 的 $n'$ 个特征向量组成的矩阵，而 $-\lambda$ 为 $XX^T$ 的若干特征值组成的矩阵，特征值在主对角线上，其余位置为0。当我们将数据集从 $n$ 维降到 $n'$ 维时，需要找到最大的 $n'$ 个特征值对应的特征向量。这 $n'$ 个特征向量组成的矩阵 $W$ 即为我们需要的矩阵。对于原始数据集，我们只需要用 $z^{(i)}=W^Tx^{(i)}$ ,就可以把原始数据集降维到最小投影距离的 $n'$ 维数据集。

4. PCA算法流程

$\qquad$ 从上面两节我们可以看出，求样本 $x^{(i)}$ 的 $n'$ 维的主成分其实就是求样本集的协方差矩阵 $XX^T$ 的前 $n'$ 个特征值对应特征向量矩阵 $W$ ，然后对于每个样本 $x^{(i)}$ ,做如下变换 $z^{(i)}=W^Tx^{(i)}$ ，即达到降维的PCA目的。

$\qquad$ 下面我们看看具体的算法流程。

$\qquad$ 输入： $n$ 维样本集 $D = (x^{(1)},x^{(2)},\dots,x^{(n)})$ ,要降维到的维数 $n'$ .

$\qquad$ 输出：降维后的样本集 $D'$

$\qquad$ 1) 对所有的样本进行中心化： $x^{(i)} = x^{(i)} -\dfrac{1}{m} \sum\limits_{j=1}^mx^{(j)}$

$\qquad$ 2) 计算样本的协方差矩阵 $XX^T$

$\qquad$ 3) 对矩阵 $XX^T$ 进行特征值分解

$\qquad$ 4）取出最大的n’个特征值对应的特征向量 $(w_1,w_2,\dots,w_{n'})$ ，将所有的特征向量标准化后，组成特征向量矩阵 $W$ 。

$\qquad$ 5）对样本集中的每一个样本 $x^{(i)}$ ,转化为新的样本 $x^{(i)} = W^Tx^{(i)}$ 。

$\qquad$ 6) 得到输出样本集 $D = (z^{(1)},z^{(2)},\dots,z^{(m)})$

$\qquad$ 有时候，我们不指定降维后的 $n'$ 的值，而是换种方式，指定一个降维到的主成分比重阈值 $t$ 。这个阈值t在 $(0,1]$ 之间。假如我们的 $n$ 个特征值为 $\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_{n'}$ ，则 $n'$ 可以通过下式得到:
$\dfrac{\sum\limits_{i=1}^{n'} \lambda_i}{\sum\limits_{i=1}^{n} \lambda_i} \ge t$

5. PCA实例

$\qquad$ 下面举一个简单的例子，说明PCA的过程。
$\qquad$ 假设我们的数据集有10个二维数据(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)，需要用PCA降到1维特征。

$\qquad$ 首先我们对样本中心化，这里样本的均值为(1.81, 1.91),所有的样本减去这个均值后，即中心化后的数据集为(0.69, 0.49), (-1.31, -1.21), (0.39, 0.99), (0.09, 0.29), (1.29, 1.09), (0.49, 0.79), (0.19, -0.31), (-0.81, -0.81), (-0.31, -0.31), (-0.71, -1.01)。

import numpy as np

data = [(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)]

data = np.array([[t[0], t[1]] for t in data])
average = np.mean(data, axis=0)
data = data - average
print(data)
"""
[[ 0.69  0.49]
 [-1.31 -1.21]
 [ 0.39  0.99]
 [ 0.09  0.29]
 [ 1.29  1.09]
 [ 0.49  0.79]
 [ 0.19 -0.31]
 [-0.81 -0.81]
 [-0.31 -0.31]
 [-0.71 -1.01]]
"""

$\qquad$ 现在我们开始求样本的协方差矩阵，由于我们是二维的，则协方差矩阵为：
$XX^T =\left(\begin{matrix} cov(x_1,x_1) & cov(x_1,x_2) \\ cov(x_2,x_1) & cov(x_2,x_2) \end{matrix}\right)$
$\qquad$ 对于我们的数据，求出协方差矩阵为：
$XX^T =\left(\begin{matrix} 0.61655556 & 0.61544444 \\ 0.61544444 & 0.7165555 \end{matrix}\right)$

cov = np.cov(data[:, 0], data[:, 1])
print(cov)
"""
[[0.61655556 0.61544444]
 [0.61544444 0.71655556]]
"""

$\qquad$ 求出特征值为 $(0.490833989,1.28402771)$ ，对应的特征向量分别为： $(-0.73517866,0.6778734)^T$ 和 $(-0.6778734,-0.73517866)^T$

eigenvalue, featurevector = np.linalg.eig(cov)
print(eigenvalue)
"""
[0.0490834  1.28402771]
"""
print(featurevector) # 该特征向量为列向量
"""
[[-0.73517866 -0.6778734 ]
 [ 0.6778734  -0.73517866]]
"""

$\qquad$ 由于最大的 $k=1$ 个特征值为1.28402771，对于的 $k=1$ 个特征向量为 $(-0.6778734,-0.73517866)^T$ ，则我们的 $W=(-0.6778734,-0.73517866)^T$

$\qquad$ 我们对所有的数据集进行投影 $z^{(i)}=W^Tx^{(i)}$ ，得到PCA降维后的10个一维数据集为：(-0.827970186， 1.77758033， -0.992197494， -0.274210416， -1.67580142， -0.912949103， 0.0991094375， 1.14457216, 0.438046137， 1.22382056)

z = np.matmul(data, featurevector[:, 1].T)
print(z)

6. 核主成分分析KPCA介绍

$\qquad$ 在上面的PCA算法中，我们假设存在一个线性的超平面，可以让我们对数据进行投影。但是有些时候，数据不是线性的，不能直接进行PCA降维。这里就需要用到和支持向量机一样的核函数的思想，先把数据集从 $n$ 维映射到线性可分的高维 $N>n$ ,然后再从 $N$ 维降维到一个低维度 $n'$ , 这里的维度之间满足 $n'<n<N$ 。

$\qquad$ 使用了核函数的主成分分析一般称之为核主成分分析(Kernelized PCA, 以下简称KPCA。假设高维空间的数据是由n维空间的数据通过映射 $\phi$ 产生。

$\qquad$ 则对于n维空间的特征分解：
$\sum_{i=1}^mx^{(i)}x^{(i)T}W = \lambda W$

$\qquad$ 映射为：
$\sum_{i=1}^m\phi(x^{(i)})\phi(x^{(i)})^TW = \lambda W$

$\qquad$ 通过在高维空间进行协方差矩阵的特征值分解，然后用和PCA一样的方法进行降维。一般来说，映射 $\phi$ 不用显式的计算，而是在需要计算的时候通过核函数完成。由于KPCA需要核函数的运算，因此它的计算量要比PCA大很多。

7. PCA算法总结

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

$\qquad$ PCA算法的主要优点有：

$\qquad$ 1）仅仅需要以方差衡量信息量，不受数据集以外的因素影响。
　
$\qquad$ 2）各主成分之间正交，可消除原始数据成分间的相互影响的因素。

$\qquad$ 3）计算方法简单，主要运算是特征值分解，易于实现。

$\qquad$ PCA算法的主要缺点有：

$\qquad$ 1）主成分各个特征维度的含义具有一定的模糊性，不如原始样本特征的解释性强。

$\qquad$ 2）方差小的非主成分也可能含有对样本差异的重要信息，因降维丢弃可能对后续数据处理有影响。