机器学习（四）——PCA

在机器学习（三）——梯度下降中记录了一种寻找函数局部最小值的方法：梯度下降法，在本篇中将记录一种非监督的机器学习算法——主成分分析PCA。

1.引言

有时我们收集到的数据有很多特征，可能有成百上千个特征，那么使用这些数据进行训练就会非常耗时并且会占用大量的内存。此时就可以对数据进行降维。

1.1 数据降维

从低维度来理解数据降维，可以把数据降维理解为投影：
（1）比如有三个特征的样本数据，现希望将其降维到二维，就可以理解成这些三维空间中的点投影到一个二维平面上，形成新的二维平面上的二维样本数据；
借用吴恩达老师PPT上的图片：

（2）再比如有两个特征的样本数据，现希望将其降维到一维，就可以理解成这些二维平面上的点投影到一条轴上，形成新的一维轴向上的一维样本数据。
借用吴恩达老师PPT上的图片：

1.2 数据降维方法

大致了解了数据降维，咱来了解一下数据降维的方法。
最常用的数据降维方法就是PCA（主成分分析），PCA是一种非监督的机器学习方法。通过数据降维，不仅可以将数据压缩，也有利于数据的可视化。
顾名思义，主成分分析就是将n维的数据降到k维，即提取出前k个最主要的特征。下面就介绍一下如何利用主成分分析进行数据降维。

2.PCA算法原理

由于吴恩达老师的课程中没有详细介绍PCA算法的具体数学原理，我刚开始一直搞不懂为什么要计算协方差矩阵、为什么要计算协方差矩阵的特征向量，为了更深入的理解PCA算法，我查阅了一些关于PCA算法数学原理的资料，如果对其数学原理感兴趣的可以参考这篇文章：链接
以下内容是我依据文章内容的个人的理解和总结，如果不想了解PCA的数学原理，想直接从PCA的实现开始的话，可以直接从3.PCA算法步骤开始看。

2.1 基

假设二维向量$\vec{n}$ = (1, 2),我们如何确定这个向量在二维平面上的表示呢？我们一般都会将向量表示到二维笛卡尔直角坐标系中，我们会认为$\vec{n}$ = (1, 2)是一个从(0,0)指向(1,2)的有向线段，其中1就是向量在x轴上的投影值，2就是向量在y轴上的投影值。其实这是我们默认了以(1,0)和(0,1)为基来表示向量。

什么是基？二维向量(x,y)可以表示为这样的线性组合：$x(1,0{)}^T+y(0,1{)}^T$，把(1,0)和(0,1)称作二维平面上的一组基。因此$\vec{n}$ = (1, 2)其实就是在基(1,0)(0,1)上的投影值组合。

因此，在描述向量时，我们首先要确定一组基，在基所在的各个直线方向上给出投影值，就可以确定一个向量。只不过根据我们的习惯，通常会略过选择基这一步骤，通常会默认以(1,0)和(0,1)为基，因为这样比较方便直观，这样描述的向量就和直角坐标系中的坐标值一一对应。

2.2 基变换

在二维平面中，任意两个线性无关的二维向量都可以作为一组基（线性无关可以理解为这俩向量所在的直线不重合）。一般我们会希望基的模为1，因为这样的话，通过向量点乘基就可以获得向量在新的基上的新坐标了。

例如：选择的基为($\frac{1}{\sqrt{2}}$, $\frac{1}{\sqrt{2}}$)和(-$\frac{1}{\sqrt{2}}$, $\frac{1}{\sqrt{2}}$)，原始向量为(3,2)，那么通过向量点乘基就可以获得新的坐标：
$$\left(\begin{array}{cc}\frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}\\ -\frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}\end{array}\right)\left(\genfrac{}{}{0ex}{}{3}{2}\right)=\left(\genfrac{}{}{0ex}{}{\frac{5}{\sqrt{2}}}{-\frac{1}{\sqrt{2}}}\right)$$

由此可以进行推广：m个n维向量，将其变换为R个n维向量表示的新空间中，其转换步骤为：
（1）将R个基按行堆叠组成矩阵$A$ ($A$为Rn矩阵）；
（2）将原向量按列堆叠组成矩阵$B$ （$B$为nm矩阵）;
（3）$AB$就是变换结果，即m个n维向量在新空间$R_n$

我们还可以根据这个转换过程来理解矩阵相乘的几何表示。例如：$AB$就相当于将$B$中的每一列列向量变换到$A$中每一行行向量为基的新空间中。

2.3 如何降维和问题转换

了解了基和基变换之后，我们来考虑如何完成降维的问题。答案显而易见啦，当基的数量少于向量本身的维数时，就可以实现降维。比如基的个数为$k$，向量的维数为$n$，其中$k<n$，通过基变换就可以实现将向量从$n$维降到$k$维。

但是新的问题出现了：如何选择这$k$个基，从而使原始信息能最大程度的保留？

我们来看这几张图片：其中(a)表示原始数据，$x_1$和$x_2$表示特征，说明数据是有两个的特征的二维数据

（b）（c）都是将二维数据降到一维，但是我们会觉得（b）的降维效果更好，因为降维后的数据较为分散，数据间的可区分度较高（反观（c）其中有两个降维后的数据离得较近，可区分度较低）。再来看看（b）和（d），我们又会觉得（d）的效果更好，因为降维后的数据不仅比较分散，而且和原始数据的分布差别较小。

考虑了上述问题，首先我们来解决降维后数据的可区分度的问题。如何来描述数据的分散程度，是不是想到了方差，此处附上百度对方差的解释：方差是在概率论和统计方差衡量随机变量或一组数据时离散程度的度量。 哇那么问题一下子清晰明了了，我们就将原始问题转化为：寻找一组基，使得数据转换到这个基上的新坐标之后，方差最大 假设降维后的一个数据特征为$X$（一个数据特征对应原数据矩阵中的一列），则其方差表示为：
$$Var(X)=\frac{1}{m}\sum _{i=1}^m(x_i-\overline{X}{)}^2$$
为了方便处理，我们先处理数据使其均值为0，那么方差就变成了：
$$Var(X)=\frac{1}{m}\sum _{i=1}^m(x_i{)}^2$$
均值归0后的数据我们就可以形象的表示为上图中的（e）。

2.4 协方差

在$2.3$中我们考虑了二维数据降维到一维，我们只需要找到一个方向，使得投影后的数据方差最大。将问题推广到高维，我们在选择了第一个轴向之后，如何选择后续的轴向？如果我们继续选择方差最大的轴向，那么这就和第一个轴向几乎重合了，这是没有意义的。

我们当然希望降维后的数据能够最大程度的保留原始信息，这就要求在选择后续的轴向时，我们希望这些方向是无相关性的，因为有相关性就表明它们之间不是完全独立的，就会有重复表示的信息。此时引入协方差，在数学上用协方差来表示两个随机变量的相关性。此处引用百度上对协方差的解释：协方差为0的两个随机变量称为是不相关的。

协方差的计算公式为：

将其引入到我们的例子中，假设我们的两个数据为$X^{(1)}$和$X^{(2)}$，我们已将其均值归0，这时就显出了这样处理的好处啦，$X^{(1)}$和$X^{(2)}$的协方差就表示为：
$$Cov(X^{(1)},X^{(2)})=\frac{1}{m}\sum _{i=1}^m{x}_i^{(1)}{x}_i^{(2)}$$
这里插一句：那个上标表示的是第几个样本数据，下标表示样本数据的第几个特征值。
为了使两数据相互独立，则需使$Cov(X^{(1)},X^{(2)})=0$ ,那么第二个方向就要在第一个方向正交的方向中选择，也就是说第二个方向一定与第一个方向正交。

此时我们又可以将问题转化为：选择k个单位正交基（正交基的模均为1），使得原始数据变换到这组基后，各个特征间的协方差为0，而且特征的方差尽可能大。

通过上面的分析，我们引入了方差和协方差，我们现在就想能不能统一表示特征内方差和特征间协方差？接下来介绍协方差矩阵。

2.5 协方差矩阵

我们先来举个例子，假设现有m个二维数据组成的矩阵$X$:
$$X=\left(\begin{array}{cccc}{x}_1^{(1)}& x_1^{(2)}& ...& x_1^{(m)}\\ x_2^{(1)}& x_2^{(2)}& ...& x_2^{(m)}\end{array}\right)$$
我们做这样的运算，看看会发生什么：
$$C=\frac{1}{m}X{X}^T=\left(\begin{array}{cc}\frac{1}{m}{\sum }_{i=1}^m(x_1^{(i)}{)}^2& \frac{1}{m}{\sum }_{i=1}^m{x}_1^{(i)}{x}_2^{(i)}\\ \frac{1}{m}{\sum }_{i=1}^m{x}_1^{(i)}{x}_2^{(i)}& \frac{1}{m}{\sum }_{i=1}^m(x_2^{(i)}{)}^2\end{array}\right)$$
观察矩阵$C$我们发现，$C$的对角线上的元素对应的是第一个特征和第二个特征的方差，其他的元素对应了特征间的协方差。这不就正好符合咱们刚刚的想法嘛：统一表示特征内方差和特征间协方差。

其实矩阵$C$就是矩阵$X$的协方差矩阵，可参考百度对协方差的解释。协方差的每个元素对应$X$的分量$X_i$和$X_j$的协方差，当i=j时该元素就对应了分量的方差。

引入协方差矩阵之后，我们的目标是：
（1）协方差矩阵对角化
即除了对角线元素之外，其余元素为0。这就满足了我们希望特征间协方差为0的目的。
（2）对角线元素排序
将对角线上元素按行从大到小排序，这样前k行就是我们要选择的k个单位正交基

基于这两个目标，我们做进一步的推导。

2.6 协方差矩阵的特征向量

先做如下假设：
$C$：原始数据矩阵对应的协方差矩阵；
$P$：一组基按行排列组成的矩阵；
$Y=PX$：$Y$是$X$对$P$做基变换后的新数据；
$D$：$Y$的协方差矩阵。
来推导$C$与$D$之间的关系：

$$D=\frac{1}{m}Y{Y}^T=\frac{1}{m}(PX)(PX{)}^T=\frac{1}{m}PX{X}^T{P}^T=P(\frac{1}{m}X{X}^T)P^T=PC{P}^T$$
由此可知，我们要寻找一个矩阵$P$使得$PC{P}^T$是一个对角矩阵，而且该对角矩阵的元素按行从大到小依次排列，则$P$的前k行就是我们要找的k个基，用这k个基组成的新矩阵$P_k$点乘$X$，就使得$X$从n维降到了k维。

现在问题就变成了如何寻找矩阵$P$？

我们知道$C$为实对称矩阵，根据实对称矩阵的性质：
（1）实对称矩阵的不同特征值对应的特征向量是正交的。
（2）n阶实对称矩阵A必可相似对角化，且相似对角阵上的元素即为矩阵本身特征值。
可知：一个n*n的实对称矩阵，一定可以找到n个单位正交特征向量$\overrightarrow{e_1},\overrightarrow{e_2},...,\overrightarrow{e_n}$,使得$E=\left(\begin{array}{c}\overrightarrow{e_1},\overrightarrow{e_2},...,\overrightarrow{e_n}\end{array}\right)$，对$C$有
$$E^{T}CE=\left(\begin{array}{cccc}{\lambda }_1& & & \\ & {\lambda }_2& & \\ & & ...& \\ & & & {\lambda }_n\end{array}\right)$$
其中$E^{T}CE$为对角矩阵，对角元素为各个特征向量对应的特征值。

由$E^{T}CE=PC{P}^T$可得：$P=E^T$。将$P$按照矩阵中特征值从大到小排列，则前k行组成的新矩阵$P_k$就是目标矩阵。则$Y=P_{k}X$，$Y$就是$X$从n维降到k维的新数据。

3.PCA算法步骤

根据上面的分析，我们可以总结出PCA算法的步骤：
（1）将原始数据按列排列组成n*m的矩阵$X$；
（2）对$X$中的每一行数据进行均值归0处理，即对每个特征进行均值归0处理（一行对应一个特征）；
（3）求协方差矩阵$C=\frac{1}{m}{\sum }_{i=1}^n(x^{(i)})(x^{(i)}{)}^T$
（4）求$C$的特征值及对应的特征向量；
（5）根据特征值的大小从上到下按行排列获得矩阵$P$，取前k行构成矩阵$P_k$;
（6）$Y=P_{k}X$，$Y$就是$X$从n维降到k维的新数据。

4.PCA算法的实现

说明：代码编写及测试均在notebook中进行。

4.1 生成随机数据

随机生成100个二维数据：

X = np.empty((100, 2))
X[:, 0] = np.random.uniform(0., 100., size=100)
X[:, 1] = 0.5 * X[:, 0] + 5. + np.random.normal(0, 5, size=100)

plt.xlabel('x1')
plt.ylabel('x2')
plt.plot(X[:, 0], X[:, 1],'o')

print(X)

数据大致如下：

将数据表示到二维平面上：（可以理解为第一维数据x1就是第一个特征，第二维数据x2就是第二个特征）

4.2 均值归零处理

def demean(X):
    return X - np.mean(X, axis=0)

X_demean = demean(X)
print(X_demean[:,0].mean())
print(X_demean[:, 1].mean())

可以看到处理后特征值的均值为0：
这里要说一下，如果特征在不同的数量级上，我们还需要将其除以标准差：

X_std = np.std(X)
X_std = X_demean / X_std

查看处理后的数据在二维平面上的分布：

4.3 计算协方差矩阵

X = X_std
cov = np.cov(X[:, 0], X[:, 1])
print(cov)

4.4 计算协方差矩阵的特征值和特征向量

eig_val, eig_vec = np.linalg.eig(cov)

查看特征向量在二维平面上的表示：

plt.plot(X[:, 0], X[:, 1], 'o')

plt.plot([eig_vec[:,0][0],0],[eig_vec[:,0][1],0],color='red')
plt.plot([eig_vec[:,1][0],0],[eig_vec[:,1][1],0],color='red')

可以看到两个特征向量几乎是正交的：

验证一下是否正交：

eig_vec[:,0][0]*eig_vec[:,1][0] + eig_vec[:,0][1]*eig_vec[:,1][1]

结果为0，说明两个特征向量正交：

4.5 计算矩阵$P$和$P_k$

P = [(np.abs(eig_val[i]), eig_vec[:,i]) for i in range(len(eig_val))]
P.sort(reverse=True)
print(P)

获得的$P$是这样的，特征值从上到下按大小依次排列：

获取$P$的前k行构成矩阵$P_k$，由于我们这里处理的是二维数据，所以这里选择k为1，表示将二维数据转化为一维数据：

# 将array类型转换为matrix类型
P_k = np.matrix(P[0][1])
P_k.shape

$P_k$就表示为由一个二维基组成的矩阵。

4.6 计算降维后的数据

我们用符号X_reduction表示降维后的数据：

X = np.matrix(X)
X_reduction = (P_k.dot(X.T)).T
print(X_reduction.shape)

可以看到$X$降到一维了：

将数据表示到二维平面上：

plt.plot(X[:, 0], X[:, 1], 'o')

plt.plot([eig_vec[:,0][0],0],[eig_vec[:,0][1],0],color='red')
plt.plot([eig_vec[:,1][0],0],[eig_vec[:,1][1],0],color='green')
plt.plot(X_reduction[:,0], [1]*100, '+')

其中，蓝色的点是原始数据，红色和绿色的线分别代表两个特征向量，而且红色的线是我们选出来的投影方向，橙色的点（就是X_reduction对应的点）是蓝色点投影到红色方向上后的点，为了方便在图中查看，此处将橙色点的纵坐标都设置为1。需要注意的是，虽然这里将橙色点X_reduction表示在了二维平面中，但是这些橙色点其实是蓝色点在以$P_k$（就是红色方向）为基的新的空间中的新坐标。

通过以上实验步骤，我们可以看到PCA对二维数据的降维过程，接着我们尝试封装一下PCA算法。

4.7 PCA算法封装

4.7.1 封装

呜呼！经过上面的测试，我们现在来试着封装一下PCA：

import numpy as np
def my_PCA(X, k):
    # X的第一位代表样本个数，第二维代表样本的特征数量
    sample_num, feature_num = X.shape
    # 均值归零处理
    mean = np.array([np.mean(X[:, i]) for i in range(feature_num)])
    X_demean = X - mean
    X = X_demean
    # 求协方差矩阵
    # 这里注意：由于X的第一维表示样本个数，第二维表示样本特征数量，所以这里是X的转置点乘X
    cov = X.T.dot(X)
    # 求特征值和特征矩阵
    eig_val, eig_vec = np.linalg.eig(cov)
    # 求矩阵P
    P =[(np.abs(eig_val[i]), eig_vec[:,i]) for i in range(feature_num)]
    P.sort(reverse=True)    
    # 矩阵P_k，P_k是一个k*n的矩阵，n为样本特征数量
    P_k = np.array([p[1] for p in P[:k]])
    # 降维后的数据，X_reduction为m*k的矩阵，其中m为样本个数
    X_reduction = (P_k.dot(X.T)).T
    
    return X_reduction

4.7.2 测试

对封装的PCA算法做个测试：
（1）随机生成样本数据；
（2）使用sklearn中的PCA方法对样本数据进行降维；
（3）使用我们封装的PCA算法进行降维测试，并与（2）中的测试结果做比较。

from sklearn.decomposition import PCA
import numpy as np
X = np.empty((10, 2))
X[:, 0] = np.random.uniform(0., 10., size=10)
X[:, 1] = 0.5 * X[:, 0] + 5. + np.random.normal(0, 5, size=10)
pca=PCA(n_components=1)
pca.fit(X)
pca.transform(X)

在sklearn的PCA测试结果为：

使用我们封装的的PCA进行测试：

my_PCA(X, 1)

可以看到数据一致只是投影方向不同，可见我们封装的PCA算法是可行的。

接着再来测试一下三维数据降到二维的效果：
sklearn中的PCA：

from sklearn.decomposition import PCA
import numpy as np
X = np.empty((10, 3))
X[:, 0] = np.random.uniform(0., 10., size=10)
X[:, 1] = 0.5 * X[:, 0] + 5. + np.random.normal(0, 5, size=10)
X[:, 2] = 0.75 * X[:, 1] + 3. + np.random.normal(0, 10, size=10)
pca=PCA(n_components=2)
pca.fit(X)
pca.transform(X)

测试结果：

测试我们的PCA：

my_PCA(X, 2)

5.总结

至此，PCA算法的数学原理，PCA算法的步骤以及实现就基本总结完了。

（1）首先从空间转平面、平面转直线，我们理解了数据降维；
（2）之后从数据降维的角度来理解PCA，我们知道PCA算法可以用来数据降维，并了解了其数学原理，从向量的表示和转换来理解PCA的降维过程是如何实现的；
（3）了解了PCA算法的数学原理之后，就明白了为什么要计算协方差矩阵、为什么要计算协方差矩阵的特征向量，由此就可以总结出PCA算法的步骤；
（4）根据PCA算法的步骤，我们尝试用代码实现PCA算法的降维过程，封装了PCA算法并对封装的算法进行了测试，测试结果表明我们封装的算法是可行的。

至此本篇也接近尾声了，有几点要说明一下：
（1）因为本人也是初学者，文章内容可能会有错误和不足，非常欢迎大家指正和补充！
（2）PCA算法也不单单只用来进行数据降维，它还有很多别的用途，但在本文中就没有具体讲述了，希望了解PCA算法的其他用途可查阅相关资料。

好啦，希望这篇文章能够帮助大家理解PCA，也非常期待大家的指正和补充，期待与大家交流哟！

机器学习（五）——逻辑回归