凌云时刻 · 技术
导读：这篇笔记主要讲解机器学习中经常用到的降维算法PCA。
作者 | 计缘
来源 | 凌云时刻（微信号：linuxpk）
前言
在机器学习的实际使用中，我们都希望有足够多的样本数据，并且有足够的特征来训练我们的模型，所以高维特征数据是经常会用到的，但是高维特征数据同样会带来一些问题：
机器学习算法收敛速度下降。
特征难于分辨，很难第一时间认识某个特征代表的意义。
会产生冗余特征，增加模型训练难度，比如说某一品牌型号汽车的特征数据，有从中国采集的，也有从国外采集的，那么就会产生公里/小时和英里/小时这种特征，但其实这两个特征代表的意义是一样的。
无法通过可视化对训练数据进行综合分析。
以上问题都是高维特征数据带来的普遍问题，所以将高维特征数据降为低维特征数据就很重要了。这篇笔记主要讲解机器学习中经常用到的降维算法PCA。
PCA是英文Principle Component Analysis的缩写，既主成分分析法。该算法能从冗余特征中提取主要成分，在不太损失模型质量的情况下，提升了模型训练速度。
理解PCA算法降维的原理
我们从二维降一维的场景来理解PCA降维的原理。上面的图示显示了一个二维的特征坐标，横坐标是特征1，纵座标是特征2。图中的五个点就表示了五条特征数据。我们先来想一下最简单粗暴的降维方式就是丢弃掉其中一个特征。
如上图中显示，将特征2抛弃，这里大家先注意一下这五个点落在特征1轴上的间距。
或者如上图所示抛弃特征1，大家再注意一下这五个点落在特征2轴上的间距。能很明显的发现，抛弃特征2，落在特征1轴上的五个点之间间距比较大，并且分布均匀。而抛弃特征1，落在特征2轴上的五个点之间间距大多都比较小，并且分布不均匀。
就这两种简单粗暴的降维方式而言，哪种更好一些呢？这里我们先来看看方差的概念，方差描述的是随机数据的离散程度，也就是离期望值（不严谨的说，期望值等同于均值）的距离。所以方差越大，数据的离散程度越高，约分散，离均值的距离越大。方差越小，数据的离散程度越小，约聚合，离均值的距离约小。那么大家可以想想作为机器学习算法训练的样本数据，每组特征应该尽可能的全，在该特征的合理范围内尽可能的广，这样才能更高的代表真实性，也就是每组特征数据的方差应该尽可能的大才好。所以就上面两种情况来看，抛弃特征2的降维方式更好一些。
但是简单粗暴的完全丢弃一个特征自然是不合理的，这极大的影响了训练出模型的正确性。所以，按照方差最大的理论，我们应该再找到一个特征向量，使样本数据落在这个特征向量后的方差最大，那么这个特征向量代表的特征就是我们需要的降维后的特征。这就是支撑PCA算法的理论之一。
如上图所示，降维后的特征方差明显大于抛弃特征1或抛弃特征2后的方差。
我们在使用PCA之前首先需要对样本数据进行特征去均值化，也就是将每个特征的值减去该维特征的平均值。去均值化的目的是去除均值对数据变化的影响，也就是避免第一主成分收到均值的影响。
在第二篇笔记中，我们提到过方差，它的公式为：
那么当数据去均值化后，公式中的 就成了0，所以去均值后的方差为：
此时 就是降维后的特征，我们记为 ，那么降维后特征值的方差公式为：
因为在高维特征下， 和 都是向量，所以求方差时候应该对他们取模：
所以我们就是要求上面这个公式的最大值。那么首先如何求得这个 呢？我们来具体看一下：
如上图所示，蓝色向量是特征值原始维度的向量， 黑色向量就是我们要寻求的新维度的向量，绿色的点就是蓝色点在新维度上的投影点，红色向量的长度就是投影点的特征值。下面我们先来看看初高中数学中学过的知识。
首先向量的点乘有代数定义，也有几何定义。在几何定义下，两个向量的点乘为这两个向量的长度相乘，再乘以这两个向量夹角的余弦：
所以从上图来看：
因为我们需要的 只是方向，所以它可以是一个方向向量，既长度为1，所以上面的公式就变为：
然后由三角函数可知，夹角的余弦等于邻边除以斜边。上图中 角的斜边就是 ，邻边就是 ，所以：
此时我们就知道了，我们要求得的红色向量的长度，即：
代入上面的方差公式为：
因为两个向量的点乘是一个标量，所以最终公式为：
那么我们的目标就是求 向量，使得上面的这个公式最大。上一篇笔记我们讲了用梯度下降法求函数极小值，那么这里我们就要用到梯度上升法求函数的极大值。
使用梯度上升法解决主成分分析问题
我们先将上面的公式展开（w是一个列向量）：
既然是梯度上升，那么第一步当然是求梯度了，这和梯度下降是一样的，结合第五篇笔记中的梯度推导可得，上面公式的梯度为：
下面对上面的公式再进行向量化，这里我再推导一遍，首先我们将 看成是一个1行m列的行矩阵中的元素：
然后将它和一个m行n列的矩阵相乘：
因为X是一个m行n列矩阵，w是一个n行1列的矩阵，所以X乘w是一个m行1列的矩阵，上面我们将其转换为了1行m列的矩阵，所以上面的公式简写为 ，相乘后的结果是一个1行n列的矩阵：
那我们的梯度是一个n行1列的矩阵，所以将上面的矩阵再做转置：
所以最终主成分分析的梯度向量化后为：
首先我们构建样本数据，其中有100条样本数据，每条样本数据中有2个特征：
import numpy as np
import matplotlib.pyplot as plt
# 构建样本数据
# 构建一个100行2列的矩阵
X = np.empty((100, 2))
# 第一个特征为0到100的随机分布
X[:, 0] = np.random.uniform(0., 100., size=100)
# 第二个特征和第一个特征有一定线性关系，并且增加了0到10的正态分布的噪音
X[:, 1] = X[:, 0] * 0.75 + 3. + np.random.normal(0, 10, size=100)
# 将特征绘制出来
plt.scatter(X[:, 0], X[:, 1])
plt.show()
接下来根据上文中讲到，下一步要对样本数据的每一个特征进行均值归0操作，也就是每一列的元素减去这一列所有元素的均值：
def demean(X):
# 对矩阵X在横轴方向上求均值，既求每一列的均值
return X - np.mean(X, axis=0)
# 均值归0化后的样本数据
X_demean = demean(X)
plt.scatter(X_demean[:, 0], X_demean[:, 1])
plt.show()
可以看到均值归0化后，样本数据的分布形态是没有变化的，但是坐标轴往右上移动了，既0点现在在样本数据的中间。下面来定义目标函数和梯度函数：
# 目标函数
def f(X, w):
return np.sum(X.dot(w)**2) / len(X)
# 梯度函数
def df(X, w):
return X.T.dot(X.dot(x)) / 2 * len(X)
在上面的公式推导的过程中提到过，我们期望的向量ww是一个单位向量，所以在代码实现计算的时候需要将传入的初始向量ww和计算后的新ww向量都转换为单位向量（向量的单位向量为该向量除以该向量的模）：
def direction(w):
# np.linalg.norm(w)为求向量的模
return w / np.linalg.norm(w)
接下来的梯度上升计算和梯度下降计算基本是大同小异的：
# 参数传入样本矩阵，初始向量，步长，查找循环次数，两次方差的最小差值
def gradient_ascent(X, initial_w, eta, n_iters=1e4, different=1e-8):
# 将向量转换为单位向量
w = direction(initial_w)
i_iters = 0
while i_iters < n_iters:
gradient = df(X, w)
last_w = w
w = w + eta * gradient
w = direction(w)
if(abs(f(X, w) - f(X, last_w)) < different):
break
i_iters += 1
return w
在使用我们定义的方法前，有两点需要注意的是，一点是在PCA中，初始向量ww不能为0，因为方差公式里ww在分子，所以如果为0了，那么方差始终为0，所以每次我们随机初始化一个不为0的向量。另外一点是在PCA中我们不对样本数据做归一化或标准化处理，因为一旦做了归一化处理，样本数据的方差就变成了1，自然求不到最大方差了。
# 初始化随机向量
initial_w = np.random.random(X.shape[1])
# 设置步长
eta = 0.01
# 梯度上升
w = gradient_ascent(X_demean, initial_w, eta)
# 绘制w向量
plt.scatter(X_demean[:, 0], X_demean[:, 1])
plt.plot([0, w[0]*30], [0, w[1]*30], color='r')
plt.show()
这样我们就求出了样本数据的第一个降维到的向量，我们称为样本的第一主成分。
在上节中我们使用的样本数据是在二维空间的，求出的第一主成分其实可以看作是将坐标轴旋转后横轴或纵轴，我们降维后的数据其实是新的坐标轴上某个轴的分量，那么另外一个分量自然是降维到垂直于第一主成分向量的向量，既新坐标轴的另外一个轴。该向量是第一主成分向量的正交线。那么下面我们来看一下第一主成分的正交线如何求：
从上图可以看到， 就是第一主成分向量w的正交线，由向量的加减法可得：
因为上文推导过：
所以得：
这就相当于原始样本数据减去投影到第一主成分上的数据，我们对去掉第一主成分数据的样本数据再求第一主成分数据，那么就相当于在求原始样本数据的第二主成分了，以此类推就可以求得样本数据的n个主成分。
下面我们来用代码实现，首先我们算出样本数据在新坐标轴上的另一个分量，根据上面推导出的公式可得：
X2 = X - X.dot(w).reshape(-1, 1) * w
# 绘制X2及第一主成分向量w
plt.scatter(X2[:, 0], X2[:, 1])
plt.plot([0, w[0]*30], [0, w[1]*30], color='r')
plt.show()
首先可以看到，当样本数据去掉第一主成分数据后，另一个分量的数据其实就是在正交于第一主成分向量的轴上，所以所有点都在一条直线上。
然后将之前的gradient_ascent方法改个名称，因为它就是求第一主成分的方法，所以改名为first_component，然后求出X2的第一主成分：
def first_component(X, initial_w, eta, n_iters=1e4, different=1e-8):
w = direction(initial_w)
i_iters = 0
while i_iters < n_iters:
gradient = df(X, w)
last_w = w
w = w + eta * gradient
w = direction(w)
if(abs(f(X, w) - f(X, last_w)) < different):
break
i_iters += 1
return w
w2 = first_component(X2, initial_w, eta)
由向量的正交定理知道，垂直的向量点乘结果为0，所以我们来验证一下w和w2之间的点乘结果：
w.dot(w2)
# 结果
3.2666630511712924e-10
可以看到结果基于趋近于0。
下面我们封装一个计算n个主成分的方法：
def first_n_component(n, X, eta=0.01, n_iters=1e4, different=1e-8):
# 拷贝原始样本数据
X_pca = X.copy()
# 对样本数据进行均值归一化
X_pca = demean(X_pca)
# 存储结果数组
res = []
# 希望计算几个主成分就循环几次
for i in range(n):
# 每次随机一个初始向量
initial_w = np.random.random(X_pca.shape[1])
# 通过获取主成分方法计算出主成分向量
w = first_component(X_pca, initial_w, eta)
res.append(w)
# 每次从原始样本数据中除去主成分数据
X_pca = X_pca - X_pca.dot(w).reshape(-1, 1) * w
return res
first_n_component(2, X)
# 结果
[array([ 0.77899988, 0.62702407]), array([-0.62702407, 0.77899988])]
END
往期精彩文章回顾
机器学习笔记（十二）：随机梯度下降
机器学习笔记（十一）：优化梯度公式机器学习笔记（十）：梯度下降
机器学习笔记（九）：多元线性回归
机器学习笔记（八）：线性回归算法的评测标准机器学习笔记（七）：线性回归
机器学习笔记（六）：数据归一化
机器学习笔记（五）：超参数
机器学习笔记（四）：kNN算法
机器学习笔记（三）：NumPy、Matplotlib、kNN算法长按扫描二维码关注凌云时刻
每日收获前沿技术与科技洞见
