机器学习入门（八）：主成分分析（PCA）

PCA是机器学习中常用的方法，其主要作用是降维。因为做运算的时候会遇到维度特别大的情况，如果蛮力求解会导致维度灾难。而通过降维可以有效避免这些情况的产生，同时减少运算开销。

要了解PCA首先要了解特征值及特征向量

特征值与特征向量

定义：有一个n * n的矩阵，如果存在一个非零向量$x$使得$Ax=\lambda x$，则称标量$\lambda$为特征值（Eigenvalue），而x为特征向量（Eigenvector）。

光看定义其实很抽象，到底大家常说的特征值和特征向量的本质是什么？先看这么一个问题:

某个城镇，每年30%的已婚女性离婚，且20%的单婚女性结婚。假定共有8000名已婚和2000名未婚女性，并且总人口保持不变。我们研究结婚率和离婚率保持不变时将来长时间的期望问题。

首先，很简单，一年后的女性人口比例为：
$$w_1=A{w}_0=\left[\begin{array}{cc}0.7& 0.2\\ 0.3& 0.8\end{array}\right]\left[\begin{array}{c}8000\\ 2000\end{array}\right]=\left[\begin{array}{c}6000\\ 4000\end{array}\right]$$如果觉得矩阵看着不太方便的，直接列式计算也是一样。

同理，第二年$w_2=A{w}_1=A^2{w}_0$

可以用python进行验证以下几个结果
$w_{10}=\left[\begin{array}{c}4004\\ 5996\end{array}\right]$

$w_{20}=\left[\begin{array}{c}4000\\ 6000\end{array}\right]$

$w_{30}=\left[\begin{array}{c}4000\\ 6000\end{array}\right]$

发现过了某个点之后，人口会一直保持不变，（其实这就是个马尔可夫链）。实际上，从$w_{12}$之后就一直是${\left[\begin{array}{cc}4000& 6000\end{array}\right]}^T$,并且:
$$A{w}_{12}=\left[\begin{array}{cc}0.7& 0.2\\ 0.3& 0.8\end{array}\right]\left[\begin{array}{c}4000\\ 6000\end{array}\right]=\left[\begin{array}{c}4000\\ 6000\end{array}\right]$$

A乘上一个向量，等于这个向量本身，这是不是有点眼熟，正是$Ax=\lambda x$ 只不过这里的$\lambda =1$

这里还是没有说明特征值到底是什么，接着看。

这个收敛的过程，对任意的人口分布{10000-p， p}都是成立的（p是单身人口）,也就是说，无论初始向量是否相等，最后都会得到同样的稳态向量（这里不证明这个，可自行尝试其他初始值）。选择稳态向量的倍数$x_1=(2,3{)}^T$作为一个基向量，则:
$$A{x}_1=\left[\begin{array}{cc}0.7& 0.2\\ 0.3& 0.8\end{array}\right]\left[\begin{array}{c}2\\ 3\end{array}\right]=\left[\begin{array}{c}2\\ 3\end{array}\right]=x_1$$
这里的$x1$也是个稳态向量，但是不能把同一个向量的倍数当作第二个向量，所以暂时还只有一个稳态向量。

另外一个稳态向量$x_2=(-1,1{)}^T$同样满足条件使得:
$$A{x}_2=\left[\begin{array}{cc}0.7& 0.2\\ 0.3& 0.8\end{array}\right]\left[\begin{array}{c}-1\\ 1\end{array}\right]=\left[\begin{array}{c}\frac{-1}{2}\\ \\ \frac{1}{2}\end{array}\right]=\frac{1}{2}\left[\begin{array}{c}-1\\ 1\end{array}\right]=\frac{1}{2}{x}_2$$
于是$x_1$ 和$x_2$就是一对基向量。如果把整个过程看成一个线性变换的话，那么标量1和1/2就是线性变换的自然频率。如果还是不懂，可以把这个结婚离婚的事件看成一个振动过程，无论初始状态如何，最后都会停止，也就是收敛。

只要有振动，就会有特征值，即振动的自然频率。。

主成分分析

主成分分析，Principal Component Analysis，这里的主成分其实就是指最后求到的最大的K个特征值，而这个K是我们想要达到的维度。

其主要步骤如下：
现有数据集X={x1,x2,x3...xn}X=\{x_1,x_2,x_3... x_n\}X={x1​,x2​,x3​...xn​}，我们打算将数据降到K维度.

1）去掉平均值(中心化)得到新的X

2）计算协方差矩阵$\frac{1}{n-1}X{X}^T$或者$X{X}^T$

3）获得协方差矩阵$\frac{1}{n-1}X{X}^T$或者$X{X}^T$的特征值与特征向量

4）找到最大的K个特征值，将对应的特征向量按行排列成特征向量矩阵P

5）将数据转换到K个特征向量构建的维度中，即$Y=PX$

举个例子
现在有:
$$X=\left[\begin{array}{ccccc}-1& 1& 2& -1& -1\\ -1& -1& 2& 2& -2\end{array}\right]$$
由于每一行均值都为0，不需要去平均值，则：
$$A=X{X}^T=\left[\begin{array}{ccccc}-1& 1& 2& -1& -1\\ -1& -1& 2& 2& -2\end{array}\right]\left[\begin{array}{cc}-1& -1\\ 1& -1\\ 2& 2\\ -1& 2\\ -1& -2\end{array}\right]=\left[\begin{array}{cc}8& 4\\ 4& 14\end{array}\right]$$

获得方阵之后再计算特征值与特征向量，如果给定的X已经是个方阵（n x n）那么可以直接计算特征值，不用计算协方差矩阵。

这里不会按照上面例子中的方式求，直接通过公式$Ax=\lambda x$则有
$$(A-\lambda I)x=0$$特征方程为：
$$\left|\begin{array}{cc}8-\lambda & 4\\ 4& 14-\lambda \end{array}\right|=0或{\lambda }^2-22\lambda +96=0$$
可得${\lambda }_1=16$, ${\lambda }_2=6$，无论计算中是否有$\frac{1}{n-1}$，只会使得特征值按倍数增减，但不会影响特征向量的值以及关于K的选取。

然后继续代回，现在的问题是他们对应的特征向量。

等价于求 $(A-16I)x=0$ 和 $(A-6I)x=0$ 中的x。

这里代入一下可以很快算出两个特征向量分别是
$$x_1=(1,2{)}^T,x_2=(2,-1{)}^T$$

因为任意$x_1,x_2$的倍数都可以是特征向量，所以标准化之后得：
$$x_1=(\frac{1}{\sqrt{5}},\frac{2}{\sqrt{5}}{)}^T,x_2=(\frac{2}{\sqrt{5}},-\frac{1}{\sqrt{5}}{)}^T$$

按列排列，得到特征矩阵：

$$P=\left[\begin{array}{cc}\frac{1}{\sqrt{5}}& \frac{2}{\sqrt{5}}\\ \frac{2}{\sqrt{5}}& -\frac{1}{\sqrt{5}}\end{array}\right]$$
其实到这里，都还是特征值与其向量求解的问题，接下来才是主成分分析。

假如我们想要将数据降到1维，即K=1，从大到小选取第K个特征。

于是，选取${\lambda }_1$对应的特征向量进行计算
$$y=px=(\frac{1}{\sqrt{5}},\frac{2}{\sqrt{5}})\left[\begin{array}{ccccc}-1& 1& 2& -1& -1\\ -1& -1& 2& 2& -2\end{array}\right]$$
$$=(-\frac{3}{\sqrt{5}},-\frac{1}{\sqrt{5}},\frac{6}{\sqrt{5}},\frac{3}{\sqrt{5}},-\frac{5}{\sqrt{5}})$$

这个操作就相当于把X中的坐标点$(-1,-1),(1,-1),(2,2),(-1,2),(-1,-2)$投影到了其中一个基向量上，在这个基向量（也可以叫主元）上的长度分别是$$-\frac{3}{\sqrt{5}},-\frac{1}{\sqrt{5}},\frac{6}{\sqrt{5}},\frac{3}{\sqrt{5}},-\frac{5}{\sqrt{5}}$$从而实现降维的目的(将二维投影到一维)。

代码验证一下：

import numpy as np
x=np.array([[-1,1,2,-1,-1],[-1,-1,2,2,-2]])
print(np.cov(x))

x_mean=x-np.mean(x)#减平均值
cov=(x_mean@x_mean.T)/(len(x[0])-1) #AAˆT/(n-1)
print(cov)

发现两种方式求出来相同，而我在上面用的$A{A}^T$正好是这个结果的4倍，也就是（n-1）。继续看特征值与特征向量：

a,b=np.linalg.eig(np.cov(x))
print(a) #eigen value 特征值
print(b) #eigen vector 特征向量

这个特征值1.5和4其实就是对应的6和16，正好（n-1）倍的关系，而同样地，它的特征向量也是把特征值小的放在前面，所以和结果有所出入，并且正负号也有点不同，因为特征向量并不是唯一的。

如果同样把特征值小的放前面，并且都取负，则:
$$P=\left[\begin{array}{cc}-\frac{2}{\sqrt{5}}& -\frac{1}{\sqrt{5}}\\ \frac{1}{\sqrt{5}}& -\frac{2}{\sqrt{5}}\end{array}\right]$$
这就正好对应代码结果，可以演算一下$-\frac{2}{\sqrt{5}}$就等于-0.89442719。

验证一下PCA:

from sklearn.decomposition import PCA
pca=PCA(n_components=1)
pca.fit(x.T)
print(pca.transform(x.T))

n_components设为1及以上时，表示要达到的维度，如果设为(0,1]区间的数，表示主成分要达到的比例。

结果如下：

这就对应着降维（投影）之后的在基向量上的长度：
$$-\frac{3}{\sqrt{5}},-\frac{1}{\sqrt{5}},\frac{6}{\sqrt{5}},\frac{3}{\sqrt{5}},-\frac{5}{\sqrt{5}}$$

至此，PCA的一些基本介绍大致完成，关于其他的一些推导证明问题暂不讨论，实际上sklearn中PCA的实现是基于SVD（奇异值分解）的，后面会继续探讨。