机器学习——主成分分析(PCA)

1 从KNN算法到PCA

1.1 KNN算法引入

KNN算法，也称为K近邻算法，是一种非常常用的监督学习方法，其核心思想十分简单，给定某个训测试样本，基于某种距离计算方法来计算与其最近的K个邻居，然后根据K个邻居的分类来决定测试样本的分类。 K近邻算法特点包括以下几个：

1. 无需进行事先的模型训练。训练的过程就是将训练样本保存起来。
2. 因为没有训练过程，所以预测的过程要比其他的模型时间更长，复杂度更高。

1.2 简单实现一个KNN

def KNN(D,X,labels,k):
    '''
    D是训练样本集合，包含的是各个训练样本
    X 是测试样本
    labels 是训练样本对应的标签的集合
    k:选定的最近的邻居的个数
    '''
    #首先计算的是测试样本对于各个训练样本的距离
    distance = {}
    for item in D:
        dis = 0.0
        for i in range(len(X)):
            dis += (X[i][0] - item[i][0])*(X[i][0] - item[i][0])
        distance[i]=dis
    sorted(distance,key=lambda x:distance[x],reverse=True)
    #取最近的k个邻居的下标
    label_index = []
    for i in range(k):
        temp = distance.keys()
        temp = list(temp)
        label_index.append(temp[i])
    #获取最近的k个邻居的分类
    label = []
    for i in label_index:
        label.append(labels[i])
    
    label = collections.Counter(label)
    label = label.keys()
    label = list(label)
    return label[0]

1.3 KNN存在的问题

在KNN算法中，我们发现，对于一个测试样本X而言，其包含N个属性值，也就是$X=x_1,x_2,...x_N$，同样对于训练样本集合D中每一个训练样本d而言，其也包含N个属性，那么我们在计算的时候，需要计算X中的每一个属性和d中的每一个属性的距离值，那么问题来了，这种距离计算有两个缺点：

1. 当N过大，也就是属性值过多的时候，计算的复杂度过高。
2. 对于X和d中的属性而言，某些属性对于当前的分类并没有作用，甚至于因为计算这些属性，还会引起分类的错误。

那么，当属性过多的时候，我们就需要去掉一些属性，那么我们应该如何减少这些属性呢？这就引出来来的我们下面的算法——主成分分析法(PCA)

2 主成分分析法(PCA)

2.1 PCA算法引入

首先，假设我们有一些样本点，每一个样本点具有两个属性，也就是两个特征，这些样本点在二维平面上的分布如图所示：

我们现在，想要样本点的特征缩减到一个，也就是将数据的维度从2降到1，最常见的思想，就是直接去掉某一个维度。如下图所示：


上面两张图中，我们分别去掉了特征2和特征1.我们比较一下上面的两种降维方案，似乎第一种更好一些，因为方案一中各个点的距离更大，更好区分。那么，有没有一种方案，使得降维之后的各个点能够保持的距离最大，最好区分呢？

我们看这样的一条直线，当我们把点映射到这条直线的时候，也就是如下图所示：

我们可以发现，这种映射相对于方案一，各个点的距离更大，更好区分，与原来的样本分布并没有太大的区别。

那么，如何找到使得样本间间距最大的轴？

在统计学中，方差代表的是数据之间得稀疏程度，方差越大，数据越稀疏，方差越小，数据越密集。那么，根据方差的性质，我们可以选择方差最大的情况来衡量降维之后各个样本的点的分布情况，也就是在降维之后，能够使得各个样本点之间方差最大的情况就是我们要找的对应的降维的轴。方差对应的公式为：
$$Var(X_i)=\frac{1}{p}\sum _{i=1}^p(X_i-X_{mean}{)}^2$$
其中$X^{mean}$表示样本的均值。在我们的例子中，我们想要确定一个轴在二维平面上的方向$(w_1,w_2)$，使我们样本在映射之后的方差最大，用公式表示为：
$$max(Var(X_i^{new})=\frac{1}{p}\sum _{i=1}^p||X_i^{new}-X_{mean}^{new}|{|}^2)$$
其中p表示样本的总个数。

2.2 PCA的数学原理

2.1 数据取均值

对于一个样本集合而言，为了便于数据的运算,我们让
$$X=X-X_{mean}$$
这种类似于归一的操作，其目的是将在向量空间中的向量移动位置，是再次求均值的时候，结果为0，便于计算。并且不会改变向量分布的特点。

2.2 推导过程

通过上面的描述，我们已经能够确定PCA算法的目标以及如何衡量。下面我们将X的属性从2维扩展到N维。根据上面的描述，我们可以知道，PCA的本质就是选择m个方向向量，将N维的数据映射到m维度的空间。我们用数学化的形式来描述一下就是:
$$X_{new}=AX+B$$
其中X是原始的向量(N ,1)，A表示用于降维的的矩阵(N,m)，B是一个偏置项(m,1)。通过A矩阵的映射，将p维的X向量映射为m维的$X_{new}$向量。下面，我们具体来分析映射矩阵A，将A矩阵展开之后有如下的表示：
$$A=\begin{array}{ccccccc}{a}_{11}& & a_{12}& & ...& & a_{1N}\\ a_{21}& & a_{22}& & ...& & a_{2N}\\ ......& & .......& & ........& & ......\\ a_{m1}& & a_{m2}& & ...& & a_{mN}\end{array}$$
我们将每一行用$A_i$来表示有
$$A=\begin{array}{c}{A}_1\\ \\ A_2\\ \\ ...\\ A_m\end{array}$$

根据矩阵乘法有$A_{i}X=X_{new}^i$，也就是说，每一行$A_i$代表的是将X映射到新空间某一个维度的方向。整体的数学化表示就是：
$$Y=\begin{array}{c}{Y}_1\\ \\ Y_2\\ \\ ...\\ Y_m\end{array}=\begin{array}{c}{A}_1(X-X_{mean})\\ \\ A_2(X-X_{mean})\\ \\ ...\\ A_m(X-X_{mean})\end{array}$$

其中$Y_i$表示新样本集合的第i个样本点。$Y_{mean}$表示各个样本X在映射为Y之后的，新产生的各个样本点的均值。下面，我们将新样本点i进行展开有：
$$Y_i=\begin{array}{c}{A}_1(X_i-X_{mean})\\ \\ A_2(X_i-X_{mean})\\ \\ ...\\ A_m(X_i-X_{mean})\end{array}=\begin{array}{c}{Y}_{i1}\\ \\ Y_{i2}\\ \\ ...\\ Y_{im}\end{array}其中i\in [1,p]$$
其中$X_i$表示原始样本集合的第i个样本点。

根据上面的描述，我们的目标是要保证在映射之后，各个样本的点新坐标的方差最大。那么也就是：
$$max\sum _{i=1}^p||Y_i-Y_{mean}|{|}^2$$

$Y_{mean}$表示各个样本X在映射为Y之后的均值。我们从新空间的逐个方向来进行展开，首先是$A_1$方向，在该方向中，我们要求方差最大，也就是：
$$max\sum _{i=1}^p(Y_{i1}-Y_{mean1}{)}^2$$
其中$Y_{mean1}$表示的是新的样本点集合中的所有样本在第一个方向上的均值。数学化表达就是：
$$Y_{mean1}=\frac{1}{p}\sum _{i=1}^p{Y}_{i1}=\frac{1}{p}\sum _{i=1}^p{A}_1(X_i-X_{mean})=0$$
进而推导出：
$$max\sum _{i=1}^p(Y_{i1}{)}^2=\sum _{i=1}^p||A_1(X_i-X_{mean})|{|}^2$$
根据矩阵向量乘法有：
$$\begin{gathered} \sum _{i=1}^p||A_1(X_i-X_{mean})|{|}^2=\sum _{i=1}^p[A_1(X_i-X_{mean})][A_1(X_i-X_{mean}){]}^T \\ =\sum _{i=1}^p{A}_1(X_i-X_{mean})(X_i-X_{mean}{)}^T{A}_1^T \end{gathered}$$
再进一步，我们们能够发现，对于
$$\epsilon =(X_i-X_{mean})(X_i-X_{mean}{)}^T$$
其表示的就是原样本集合的协方差矩阵。这最终原始可以化简为：

$$max\sum _{i=1}^p(Y_{i1}{)}^2=\sum _{i=1}^p{A}_1\epsilon A_1^T$$

在之前我们提到过，任意的一个$A_i$表示的都是新空间中的一个方向，则为了便于计算，我们将这个表示方向的参数进行归一化。也就是我们约定：
$$A_1{A}_1^T=||A_1|{|}^2=1$$

所以，我们可以总结一下对于方向$A_1$的目标函数以及其限制条件：
$$\begin{gathered} 目标函数：max\sum _{i=1}^p(Y_{i1}{)}^2=\sum _{i=1}^p{A}_1\epsilon A_1^T \\ 限制条件：A_1{A}_1^T=||A_1|{|}^2=1 \end{gathered}$$

下面我们补充一个数学知识：

1. 拉格朗日乘子法：设给定函数z=ƒ(x)和附加条件φ(x)=0，为寻找z=ƒ(x)在附加条件下的极值点，先做拉格朗日函数 ，其中λ为参数。
   $F(x)=f(x)+\lambda \phi (x)$
   令F(x,λ)对x和y和λ的一阶偏导数等于零，即
   F’(x)=ƒ’x(x)+λφ’x(x)=0
   F’(λ)=φ(x)=0
   由上述方程组解出x,y及λ，如此求得的(x,y)，就是函数z=ƒ(x,y)在附加条件φ(x,y)=0下的可能极值点。
   若这样的点只有一个，由实际问题可直接确定此即所求的点。（内容源于百度百科）

根据，我们拉格朗日乘子法以及我们的目标函数，我们可以再次定义我们的目标函数：
$$E(A_1)=A_1\epsilon A_1^T-{\lambda }_1(A_1{A}_1^T-1)$$
为了求极值，我们可以对于$A_1$进行求导了：
$$\frac{\partial E(A_1)}{\partial A_1}=\frac{\partial (A_1\epsilon A_1^T)}{\partial A_1}-\frac{\partial ({\lambda }_1(A_1{A}_1^T-1))}{\partial A_1}$$

这个具有的求导过程我就不详细展开了，需要记住一点的就是协方差矩阵ε是对称矩阵。解出式子为：
$$\frac{\partial E(A_1)}{\partial A_1}=(2\epsilon A_1^T{)}^T-2{\lambda }_1{A}_1$$
当导数为0的时候有：
$$\epsilon A_1^T={\lambda }_1{A}_1^T$$

下面，我们再来补充一个数学知识：

特征值与特征向量：

设 A 是n阶方阵，如果存在数m和非零n维列向量 x，使得 Ax=mx 成立，则称 m 是矩阵A的一个特征值，而x称为特征向量。(来源：百度百科)

好了，根据数学定义，我们可以知道$A_1^T和{\lambda }_1$是ε的特征向量和特征值。在回到我们的最开始的目标函数：
$$\sum _{i=1}^p{A}_1\epsilon A_1^T$$
于是就有了：
$$\sum _{i=1}^p{A}_1\epsilon A_1^T=A_1(\epsilon A_1^T)=A_1{\lambda }_1{A}_1^T={\lambda }_1(A_1{A}_1^T)=\lambda$$

于是，我们就确定了${\lambda }_1$是ε最大的特征值，$A_1$是ε的最大特征向量。且有限制条件为$A_1{A}_1^T=1$

然后，我们来看第二个方向$A_2$，对于方向$A_2$，根据之前的过程，我们可以知道其目标函数和限制条件：
$$\begin{gathered} 目标函数：max\sum _{i=1}^p{A}_2\epsilon A_2^T \\ 限制条件(1)：A_2{A}_2^T=||A_2|{|}^2=1 \\ 限制条件(2)：A_1{A}_2^T=A_2{A}_1^T=0 \end{gathered}$$

这里，我解释以下限制条件2，对于原始的样本点$X_i=[x_1,x_2,....x_n]$而言，其每一个属性值是相对于一个坐标轴而言的，就类似于我们之前给出图示：

对于每一个特征点，其坐标是根据特征1和特征2的坐标轴而言的。那么对于我们映射之后的样本点$Y_i$，其每一个特征值也是根据坐标轴来的，对于一个坐标系而言，我们要求改坐标系中的每一个坐标轴之间是垂直的，那么也就是不同坐标轴对应的方向向量的内积为0，这也就是限制条件2产生的原因。

对于$A_2$的目标函数，我们也采用拉格朗日乘子法对其进行变换：
$$E(A_2)=E(A_1)=A_2\epsilon A_2^T-{\lambda }_2(A_2{A}_2^T-1)-\beta A_1{A}_2^T$$
对$A_2$求导之后有：
$$\frac{\partial E(A_2)}{\partial A_2}=0则有：(2\epsilon A_2^T{)}^T-2{\lambda }_2{A}_2-\beta A_1=0$$

根据ε是一个对称矩阵，我们继续推导有：
$$(2\epsilon A_2^T{)}^T-2{\lambda }_2{A}_2-\beta A_1=0$$
$$(2A_2{\epsilon }^T-2{\lambda }_2{A}_2-\beta A_1)A_1^T=0$$
$$2A_2(\epsilon A_1^T)-2{\lambda }_2(A_2{A}_1^T)-\beta A_1{A}_1^T=0$$
$$2A_2({\lambda }_1{A}_1^T)-\beta =0$$
$$2{\lambda }_1{A}_2{A}_1^T=\beta$$
则有β=0。则对于原式子有：
$$(2\epsilon A_2^T{)}^T-2{\lambda }_2{A}_2-\beta A_1=(2\epsilon A_2^T{)}^T-2{\lambda }_2{A}_2=0$$
则有：
$$\epsilon A_2^T={\lambda }_2{A}_2^T$$
则求解出${\lambda }_2,A_2^T$是ε的第二大特征值和对应的特征向量。

进而推导出
${\lambda }_3,{\lambda }_4,.....{\lambda }_m，A_3^T,A_4^T,....A_m^T$分别是其第3,4，…，m大的特征值和对应特征向量。

2.3 PCA算法流程

1. 求协方差矩阵$\epsilon ={\sum }_{i=1}^P(X_i-X_{mean})(X_i-X_{mean}{)}^T$
2. 求ε的特征值，并且按照大小进行排序${\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4,.....{\lambda }_m$和对应的特征向量$A_1^t,A_2^T,A_3^T,A_4^T,....A_m^T$
3. 归一化所有的$A_i$，使得$A_i{A}_i^T=1$
4. 将$A_i$组合成对应的A矩阵
5. 映射X到Y。