机器学习笔记七之主成分分析法（PCA）、人脸识别应用

在机器学习的实际使用中，我们都希望有足够多的样本数据，并且有足够的特征来训练我们的模型，所以高维特征数据是经常会用到的，但是高维特征数据同样会带来一些问题：

机器学习算法收敛速度下降。
特征难于分辨，很难第一时间认识某个特征代表的意义。
会产生冗余特征，增加模型训练难度，比如说某一品牌型号汽车的特征数据，有从中国采集的，也有从国外采集的，那么就会产生公里/小时和英里/小时这种特征，但其实这两个特征代表的意义是一样的。
无法通过可视化对训练数据进行综合分析。

以上问题都是高维特征数据带来的普遍问题，所以将高维特征数据降为低维特征数据就很重要了。这篇笔记主要讲解机器学习中经常用到的降维算法PCA。

PCA是英文Principle Component Analysis的缩写，既主成分分析法。该算法能从冗余特征中提取主要成分，在不太损失模型质量的情况下，提升了模型训练速度。

理解PCA算法降维的原理

我们从二维降一维的场景来理解PCA降维的原理。上面的图示显示了一个二维的特征坐标，横坐标是特征1，纵座标是特征2。图中的五个点就表示了五条特征数据。我们先来想一下最简单粗暴的降维方式就是丢弃掉其中一个特征。

如上图中显示，将特征2抛弃，这里大家先注意一下这五个点落在特征1轴上的间距。

或者如上图所示抛弃特征1，大家再注意一下这五个点落在特征2轴上的间距。能很明显的发现，抛弃特征2，落在特征1轴上的五个点之间间距比较大，并且分布均匀。而抛弃特征1，落在特征2轴上的五个点之间间距大多都比较小，并且分布不均匀。

就这两种简单粗暴的降维方式而言，哪种更好一些呢？这里我们先来看看方差的概念，方差描述的是随机数据的离散程度，也就是离期望值（不严谨的说，期望值等同于均值）的距离。所以方差越大，数据的离散程度越高，约分散，离均值的距离越大。方差越小，数据的离散程度越小，约聚合，离均值的距离约小。那么大家可以想想作为机器学习算法训练的样本数据，每组特征应该尽可能的全，在该特征的合理范围内尽可能的广，这样才能更高的代表真实性，也就是每组特征数据的方差应该尽可能的大才好。所以就上面两种情况来看，抛弃特征2的降维方式更好一些。

但是简单粗暴的完全丢弃一个特征自然是不合理的，这极大的影响了训练出模型的正确性。所以，按照方差最大的理论，我们应该再找到一个特征向量，使样本数据落在这个特征向量后的方差最大，那么这个特征向量代表的特征就是我们需要的降维后的特征。这就是支撑PCA算法的理论之一。

如上图所示，降维后的特征方差明显大于抛弃特征1或抛弃特征2后的方差。

我们在使用PCA之前首先需要对样本数据进行特征去均值化，也就是将每个特征的值减去该维特征的平均值。去均值化的目的是去除均值对数据变化的影响，也就是避免第一主成分收到均值的影响。

在第二篇笔记中，我们提到过方差，它的公式为：
V a r ( X ) = ∑ i = 1 m ( X i − X ¯ ) 2 m " role="presentation">Var(X)=∑mi=1(Xi−¯¯¯¯¯X)2m V a r ( X ) = ∑ i = 1 m ( X i − X ¯ ) 2 m

那么当数据去均值化后，公式中的 X ¯ " role="presentation"> $\overline X$ 就成了0，所以去均值后的方差为：
V a r ( X ) = ∑ i = 1 m X i 2 m " role="presentation">Var(X)=∑mi=1X2im V a r ( X ) = ∑ i = 1 m X i 2 m

此时 X i " role="presentation"> $X_i$ 就是降维后的特征，我们记为 X p ( i ) " role="presentation"> $X_p^{(i)}$ ，那么降维后特征值的方差公式为：
V a r ( X p ) = ∑ i = 1 m ( X p ( i ) ) 2 m " role="presentation">Var(Xp)=∑mi=1(X(i)p)2m V a r ( X p ) = ∑ i = 1 m ( X p ( i ) ) 2 m

因为在高维特征下， X ( i ) " role="presentation"> $X^{(i)}$ 和 X p ( i ) " role="presentation"> $X_p^{(i)}$ 都是向量，所以求方差时候应该对他们取模：
V a r ( X p ) = ∑ i = 1 m | | X p ( i ) | | 2 m " role="presentation">Var(Xp)=∑mi=1||X(i)p||2m V a r ( X p ) = ∑ i = 1 m | | X p ( i ) | | 2 m

所以我们就是要求上面这个公式的最大值。那么首先如何求得这个 X p ( i ) " role="presentation"> $X_p^{(i)}$
呢？我们来具体看一下：

如上图所示，蓝色向量是特征值原始维度的向量， w" role="presentation"> $w$ 黑色向量就是我们要寻求的新维度的向量，绿色的点就是蓝色点在新维度上的投影点，红色向量的长度就是投影点的特征值。下面我们先来看看初高中数学中学过的知识。

首先向量的点乘有代数定义，也有几何定义。在几何定义下，两个向量的点乘为这两个向量的长度相乘，再乘以这两个向量夹角的余弦：
a → ⋅ b → = | a → | | b → | cos ⁡ θ" role="presentation">→a⋅→b=|→a||→b|cosθ a → ⋅ b → = | a → | | b → | cos ⁡ θ

所以从上图来看：
X ( i ) ⋅ w = | | X ( i ) | | ⋅ | | w | | ⋅ cos ⁡ θ" role="presentation">X(i)⋅w=||X(i)||⋅||w||⋅cosθ X ( i ) ⋅ w = | | X ( i ) | | ⋅ | | w | | ⋅ cos ⁡ θ

因为我们需要的 w" role="presentation"> $w$ 只是方向，所以它可以是一个方向向量，既长度为1，所以上面的公式就变为：
X ( i ) ⋅ w = | | X ( i ) | | ⋅ cos ⁡ θ" role="presentation">X(i)⋅w=||X(i)||⋅cosθ X ( i ) ⋅ w = | | X ( i ) | | ⋅ cos ⁡ θ

然后由三角函数可知，夹角的余弦等于邻边除以斜边。上图中 θ" role="presentation"> $\theta$ 角的斜边就是 | | X ( i ) | | " role="presentation"> $||X^{(i)}||$ ，邻边就是 | | X p ( i ) | | " role="presentation"> $||X_p^{(i)}||$ ，所以：
| | X p ( i ) | | = | | X ( i ) | | ⋅ cos ⁡ θ" role="presentation">||X(i)p||=||X(i)||⋅cosθ | | X p ( i ) | | = | | X ( i ) | | ⋅ cos ⁡ θ

此时我们就知道了，我们要求得的红色向量的长度，既：
| | X p ( i ) | | = X ( i ) ⋅ w" role="presentation">||X(i)p||=X(i)⋅w | | X p ( i ) | | = X ( i ) ⋅ w

代入上面的方差公式为：
V a r ( X p ) = ∑ i = 1 m | | X ( i ) ⋅ w | | 2 m " role="presentation">Var(Xp)=∑mi=1||X(i)⋅w||2m V a r ( X p ) = ∑ i = 1 m | | X ( i ) ⋅ w | | 2 m

因为两个向量的点乘是一个标量，所以最终公式为：
V a r ( X p ) = ∑ i = 1 m ( X ( i ) ⋅ w ) 2 m " role="presentation">Var(Xp)=∑mi=1(X(i)⋅w)2m V a r ( X p ) = ∑ i = 1 m ( X ( i ) ⋅ w ) 2 m

那么我们的目标就是求 w" role="presentation"> $w$ 向量，使得上面的这个公式最大。上一篇笔记我们讲了用梯度下降法求函数极小值，那么这里我们就要用到梯度上升法求函数的极大值。

使用梯度上升法解决主成分分析问题

我们先将上面的公式展开（w是一个列向量）：
V a r ( X p ) = ∑ i = 1 m ( X ( i ) ⋅ w ) 2 m = 1 m ∑ i = 1 m ( X 1 ( i ) w 1 + X 2 ( i ) w 2 + … + X n ( i ) w n ) 2 " role="presentation">Var(Xp)=∑mi=1(X(i)⋅w)2m =1mm∑i=1(X(i)1w1+X(i)2w2+…+X(i)nwn)2 V a r ( X p ) = ∑ i = 1 m ( X ( i ) ⋅ w ) 2 m = 1 m ∑ i = 1 m ( X 1 ( i ) w 1 + X 2 ( i ) w 2 + … + X n ( i ) w n ) 2

既然是梯度上升，那么第一步当然是求梯度了，这和梯度下降是一样的，结合第五篇笔记中的梯度推导可得，上面公式的梯度为：
∇ f ( w ) = [ ∂ L ∂ w 1 ∂ L ∂ w 2 … ∂ L ∂ w n ] = 2 m [ ∑ i = 1 m 2 ( X ( i ) w ) X 1 ( i ) ∑ i = 1 m 2 ( X ( i ) w ) X 2 ( i ) … ∑ i = 1 m 2 ( X ( i ) w ) X n ( i ) ] " role="presentation">∇f(w)=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∂L∂w1∂L∂w2…∂L∂wn⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦=2m⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∑mi=12(X(i)w)X(i)1∑mi=12(X(i)w)X(i)2…∑mi=12(X(i)w)X(i)n⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ ∇ f ( w ) = [ ∂ L ∂ w 1 ∂ L ∂ w 2 … ∂ L ∂ w n ] = 2 m [ ∑ i = 1 m 2 ( X ( i ) w ) X 1 ( i ) ∑ i = 1 m 2 ( X ( i ) w ) X 2 ( i ) … ∑ i = 1 m 2 ( X ( i ) w ) X n ( i ) ]

下面对上面的公式再进行向量化，这里我再推导一遍，首先我们将 X ( i ) w" role="presentation"> $X^{(i)}w$ 看成是一个1行m列的行矩阵中的元素：
[ X ( 1 ) w X ( 2 ) w X ( 3 ) w … X ( m ) w ) ] " role="presentation">[X(1)wX(2)wX(3)w…X(m)w)] [ X ( 1 ) w X ( 2 ) w X ( 3 ) w … X ( m ) w ) ]

然后将它和一个m行n列的矩阵相乘：
[ X ( 1 ) w X ( 2 ) w X ( 3 ) w … X ( m ) w ) ] ⋅ [ X 1 ( 1 ) X 2 ( 1 ) X 3 ( 1 ) … X n ( 1 ) X 1 ( 2 ) X 2 ( 2 ) X 3 ( 2 ) … X n ( 2 ) X 1 ( 3 ) X 2 ( 3 ) X 3 ( 3 ) … X n ( 3 ) … X 1 ( m ) X 2 ( m ) X 3 ( m ) … X n ( m ) ] " role="presentation">[X(1)wX(2)wX(3)w…X(m)w)]⋅⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣X(1)1X(1)2X(1)3…X(1)nX(2)1X(2)2X(2)3…X(2)nX(3)1X(3)2X(3)3…X(3)n…X(m)1X(m)2X(m)3…X(m)n⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ [ X ( 1 ) w X ( 2 ) w X ( 3 ) w … X ( m ) w ) ] ⋅ [ X 1 ( 1 ) X 2 ( 1 ) X 3 ( 1 ) … X n ( 1 ) X 1 ( 2 ) X 2 ( 2 ) X 3 ( 2 ) … X n ( 2 ) X 1 ( 3 ) X 2 ( 3 ) X 3 ( 3 ) … X n ( 3 ) … X 1 ( m ) X 2 ( m ) X 3 ( m ) … X n ( m ) ]

因为X是一个m行n列矩阵，w是一个n行1列的矩阵，所以X乘w是一个m行1列的矩阵，上面我们将其转换为了1行m列的矩阵，所以上面的公式简写为 ( X w ) T X" role="presentation"> $(Xw)^TX$ ，相乘后的结果是一个1行n列的矩阵：
( X w ) T X = [ ∑ i = 1 m ( X ( i ) w ) X 1 ( i ) ∑ i = 1 m ( X ( i ) w ) X 2 ( i ) … ∑ i = 1 m ( X ( i ) w ) X n ( i ) ] " role="presentation">(Xw)TX=[∑mi=1(X(i)w)X(i)1∑mi=1(X(i)w)X(i)2…∑mi=1(X(i)w)X(i)n] ( X w ) T X = [ ∑ i = 1 m ( X ( i ) w ) X 1 ( i ) ∑ i = 1 m ( X ( i ) w ) X 2 ( i ) … ∑ i = 1 m ( X ( i ) w ) X n ( i ) ]

那我们的梯度是一个n行1列的矩阵，所以将上面的矩阵再做转置：
( ( X w ) T X ) T = X T ( X w )" role="presentation">((Xw)TX)T=XT(Xw) ( ( X w ) T X ) T = X T ( X w )

所以最终主成分分析的梯度向量化后为：
∇ f = 2 m X T ( X w )" role="presentation">∇f=2mXT(Xw) ∇ f = 2 m X T ( X w )

代码实现PCA梯度上升

首先我们构建样本数据，其中有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)

在上面的公式推导的过程中提到过，我们期望的向量 w" role="presentation"> $w$ 是一个单位向量，所以在代码实现计算的时候需要将传入的初始向量 w" role="presentation"> $w$ 和计算后的新 w" role="presentation"> $w$ 向量都转换为单位向量（向量的单位向量为该向量除以该向量的模）：

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中，初始向量 w" role="presentation"> $w$ 不能为0，因为方差公式里 w" role="presentation"> $w$ 在分子，所以如果为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()

这样我们就求出了样本数据的第一个降维到的向量，我们称为样本的第一主成分。

求数据的其他主成分

在上节中我们使用的样本数据是在二维空间的，求出的第一主成分其实可以看作是将坐标轴旋转后横轴或纵轴，我们降维后的数据其实是新的坐标轴上某个轴的分量，那么另外一个分量自然是降维到垂直于第一主成分向量的向量，既新坐标轴的另外一个轴。该向量是第一主成分向量的正交线。那么下面我们来看一下第一主成分的正交线如何求：

从上图可以看到， X ‘ ( i ) " role="presentation"> $X^{‘(i)}$ 就是第一主成分向量 w" role="presentation"> $w$ 的正交线，由向量的加减法可得：

X ‘ ( i ) = X ( i ) − X p ( i ) " role="presentation">X‘(i)=X(i)−X(i)p X ‘ ( i ) = X ( i ) − X p ( i )

因为上文推导过：

X p ( i ) = | | X p ( i ) | | w" role="presentation">X(i)p=||X(i)p||w X p ( i ) = | | X p ( i ) | | w

| | X p ( i ) | | = X ( i ) ⋅ w" role="presentation">||X(i)p||=X(i)⋅w | | X p ( i ) | | = X ( i ) ⋅ w

所以得：

X ‘ ( i ) = X ( i ) − ( X ( i ) ⋅ w ) w" role="presentation">X‘(i)=X(i)−(X(i)⋅w)w X ‘ ( i ) = X ( i ) − ( X ( i ) ⋅ w ) 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])]

高维数据向低维数据映射

我们再来回顾一下PCA降维的基本原理，首先要做的事情就是对样本数据寻找另外一个坐标系，这个坐标系中的每一个轴依次可以表达样本数据的重要程度，既主成分。我们取出前k个主成分，然后就可以将所有的样本数据映射到这k个轴上，既获得了一个低维的数据信息。

X = [ X 1 ( 1 ) X 2 ( 1 ) … X n ( 1 ) X 1 ( 2 ) X 2 ( 2 ) … X n ( 2 ) … X 1 ( m ) X 2 ( m ) … X n ( m ) ] " role="presentation">X=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣X(1)1X(1)2…X(1)nX(2)1X(2)2…X(2)n…X(m)1X(m)2…X(m)n⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ X = [ X 1 ( 1 ) X 2 ( 1 ) … X n ( 1 ) X 1 ( 2 ) X 2 ( 2 ) … X n ( 2 ) … X 1 ( m ) X 2 ( m ) … X n ( m ) ]

上面的 X" role="presentation"> $X$ 是样本数据，该样本数据有m行，n个特征，既是一个n维的样本数据。

W k = [ W 1 ( 1 ) W 2 ( 1 ) … W n ( 1 ) W 1 ( 2 ) W 2 ( 2 ) … W n ( 2 ) … W 1 ( k ) W 2 ( k ) … W n ( k ) ] " role="presentation">Wk=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣W(1)1W(1)2…W(1)nW(2)1W(2)2…W(2)n…W(k)1W(k)2…W(k)n⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ W k = [ W 1 ( 1 ) W 2 ( 1 ) … W n ( 1 ) W 1 ( 2 ) W 2 ( 2 ) … W n ( 2 ) … W 1 ( k ) W 2 ( k ) … W n ( k ) ]

假设上面的 W" role="presentation"> $W$ 是样本数据 X" role="presentation"> $X$ 的主成分向量矩阵，每一行代表一个主成分向量，一共有k个主成分向量，每个主成分向量上有n个值。

我们已经推导出了求映射后的向量的大小，也就是每一行样本数据映射到该主成分上的大小为：

| | X p ( i ) | | = X ( i ) w" role="presentation">||X(i)p||=X(i)w | | X p ( i ) | | = X ( i ) w

那如果将一行有n个特征的样本数据分别映射到k个主成分上，既得到有k个值的新向量，既降维后的，有k个特征的新样本数据。所以我们需要的就是 X" role="presentation"> $X$ 矩阵的第一行和 W" role="presentation"> $W$ 矩阵的每一行对应元素相乘然后再相加， X" role="presentation"> $X$ 矩阵的第二行和 W" role="presentation"> $W$ 矩阵的每一行对应元素相乘然后再相加，以此类推就可以求出降维后的，m行k列的新矩阵数据：

X ‘ = X W k T " role="presentation">X‘=XWTk X ‘ = X W k T

X ‘ " role="presentation"> $X^{‘}$ 就是降维后的数据，既然可以降维，那么我们也可从数学的角度将降维后的数据还原回去。 X ‘ " role="presentation"> $X^{‘}$ 是m行k列的矩阵， W k " role="presentation"> $W_k$ 是k行n列的矩阵，所以 X ‘ W k " role="presentation"> $X^{‘} W_k$ 就是还原后的ｍ行ｎ列的原矩阵。那为什么说是从数学角度来说呢，因为毕竟已经从高维降到了低维，那势必会有丢失的数据信息，所以还原回去的数据也不可能和原始数据一样的。

在PyCharm中封装PCA

我们在myML中新建一个类PCA：

import numpy as np

class PCA:

	# 初始化PCA
	def __init__(self, n_components):
		assert n_components >= 1, "至少要有一个主成分"
		self.n_components = n_components
		self.component_ = None

	# 训练主成分矩阵
	def fit(self, X, eta=0.01, n_iters=1e4):
		assert self.n_components <= X.shape[1], "主成分数要小于等于样本数据的特征数"

		# 均值归一化
		def demean(X):
			return X - np.mean(X, axis=0)

		# 目标函数
		def f(w, X):
			return np.sum((X.dot(w)**2)) / len(X)
		# 梯度
		def df(w, X):
			return (X.T.dot(X.dot(w)) * 2) / len(X)

		# 求单位向量
		def direction(w):
			return w / np.linalg.norm(w)

		def first_component(X, initial_w, eta=0.01, n_iters=1e4, different=1e-8):

			# 转换初始向量为单位向量，既只表明方向
			w = direction(initial_w)
			cur_iters = 0

			while cur_iters < n_iters:
				# 求出梯度
				gradient = df(w, X)
				# 记录上一个方向向量
				last_w = w
				# 通过梯度上升求下一个方向向量
				w = w + eta * gradient
				# 将新求出的方向向量单位向量化
				w = direction(w)

				if(abs(f(w, X) - f(last_w, X)) < different):
					break

				cur_iters += 1

			return w

		# 对样本数据的特征数据均值归一化
		X_pca = demean(X)
		# 构建一个空的主成分矩阵，大小和样本数据保持一致
		self.component_ = np.empty(shape=(self.n_components, X.shape[1]))
		for i in range(self.n_components):
			# 随机生成一个初始向量
			initial_w = np.random.random(X_pca.shape[1])
			# 求第一主成分
			w = first_component(X_pca, initial_w, eta, n_iters)
			# 存储主成分
			self.component_[i, :] = w

			X_pca = X_pca - X_pca.dot(w).reshape(-1, 1) * w

		return self

	# 根据主成分矩阵降维样本数据
	def transform(self, X):
		assert X.shape[1] == self.component_.shape[1], "样本数据的列数，既特征数要等于主成分矩阵的列数"

		return X.dot(self.component_.T)

	# 根据主成分矩阵还原样本数据
	def inverse_transform(self, X_pca):
		assert X_pca.shape[1] == self.component_.shape[0], "降维后的样本数据特征数要等于主成分矩阵的行数"

		return X_pca.dot(self.component_)

	def __repr__(self):
		return "PCA(n_components=%d)" % self.n_components

在Jupyter Notebook中使用封装的PCA

首先构建样本数据：

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()

然后导入我们封装好的PCA类，训练主成分并根据主成分对样本数据降维：

# 导入我们封装的PCA类
from myML.PCA import PCA
# 初始化PCA类，给定只训练一个主成分
pca = PCA(n_components=1)
# 训练主成分矩阵
pca.fit(X)
# 查看主成分
pca.component_
# 结果
array([[ 0.7756218 ,  0.63119792]])

# 根据主成分对样本数据进行降维
X_reduction = pca.transform(X)
# 降维后的数据为一个100行1列的矩阵
X_reduction.shape

看到我们非常简单地就把一个二维特征的样本数据根据主成分映射为了一维特征的样本数据。同时我们还可以将其恢复二维特征数据：

# 恢复样本数据维度
X_restore = pca.inverse_transform(X_reduction)
X_restore.shape
# 结果
(100, 2)

# 将原始样本数据和从低维恢复后的样本数据绘制出来进行对比
plt.scatter(X[:, 0], X[:, 1])
plt.scatter(X_restore[:, 0], X_restore[:, 1], color='r')
plt.show()

在前面提到过，从高维降到低维就已经有部分信息丢失了，所以再恢复回去后势必不会和原始数据一样。从上图中可以看到，恢复后的二维特征数据其实是在一条直线上，而这条直线其实就是原始样本数据的主成分。

Scikit Learn中的PCA

这一节我们来看看Scikit Learn中封装的PCA如何使用，

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)

# 导入Scikit Learn中的PCA
from sklearn.decomposition import PCA
pca = PCA(n_components=1)
pca.fit(X)
X_reduction = pca.transform(X)
X_reduction.shape
# 结果
(100, 1)

X_restore = pca.inverse_transform(X_reduction)
X_restore.shape
# 结果
(100, 2)

plt.scatter(X[:, 0], X[:, 1])
plt.scatter(X_restore[:, 0], X_restore[:, 1], color='r')
plt.show()

可以看到，我们封装PCA类时，使用标准的机器学习算法的模式，所以在使用Scikit Learn提供的PCA时，几乎是一样的。

使用真实的数据

这一节我们使用真实的数据来体会一下PCA的威力。我们使用Scikit Learn中提供的手写数字数据：

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets

digits = datasets.load_digits()
X = digits.data
y = digits.target
X.shape
# 结果
(1797, 64)

可以看到，Scikit Learn提供的手写数据是一个64维特征的样本数据，一共有1797行，也就是一个1797行，64列的矩阵。

我们先用KNN分类算法来计算这个样本数据：

# 首先将样本数据拆分为训练数据集和测试数据集
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)

X_train.shape
# 结果
(1347, 64)

# 然后使用Scikit Learn的KNN算法训练分类模型，
from sklearn.neighbors import KNeighborsClassifier
knn_clf = KNeighborsClassifier()
# 在训练的时候计一下时
%time knn_clf.fit(X_train, y_train)
# 结果
CPU times: user 4.61 ms, sys: 3.25 ms, total: 7.86 ms
Wall time: 38.1 ms
KNeighborsClassifier(algorithm='auto', leaf_size=30, metric='minkowski', metric_params=None, n_jobs=1, n_neighbors=5, p=2, weights='uniform')

# 查看分类准确率
knn_clf.score(X_test, y_test)
# 0.98666666666666669

从上面的代码可以看出，使用KNN算法对样本数据进行训练时通过网格搜索的邻近点为5个，使用了明可夫斯基距离，但是p是2，所以其实还是欧拉距离，并且没有使用距离权重。训练后的分类准确率为98.7%，在我的电脑上耗时38.1毫秒。

下面我们先简单粗暴的将这个64维特征的样本数据降至2维特征数据，然后再用KNN算法训练一下看看情况：

# 使用Scikit Learn提供的PCA
from sklearn.decomposition import PCA
# 设置主成分为2
pca = PCA(n_components=2)
pca.fit(X_train)
# 将训练数据和测试数据分别降至2维
X_train_reduction = pca.transform(X_train)
X_test_reduction = pca.transform(X_test)
X_train_reduction.shape
# 结果
(1347, 2)

# 对降维后的数据进行KNN分类训练
knn_clf = KNeighborsClassifier()
%time knn_clf.fit(X_train_reduction, y_train)
# 结果
CPU times: user 1.75 ms, sys: 1.02 ms, total: 2.77 ms
Wall time: 1.77 ms

# 查看分类准确率
knn_clf.score(X_test_reduction, y_test)
# 结果
0.60666666666666669

从上面的代码和结果可以看到，首先使用KNN算法训练的耗时从64维时的38.1毫秒降至了1.77毫秒，所以这验证了PCA降维的其中的减少计算时间的作用。但是当我们查看分类准确率的时候发现非常低，所以说明我们降维度的降的太低，丢失了太多的数据信息。那么PCA中的超参数n_components应该如何取呢？其实Scikit Learn的PCA提供了一个方法就是可以计算出每个主成分代表的方差比率：

pca.explained_variance_ratio_
# 结果
array([ 0.14566817,  0.13735469])

比如通过explained_variance_ratio_我们可以知道通过PCA分析出的手写数据的前两个主成分的方差比率为14.6%和13.7%，加起来既标识降维后的数据只能保留了原始样本数据38.3%的数据信息，所以自然分类准确率很差了。那么如果我们使用PCA将64维数据计算出64个主成分，然后看看每个主成分的方差比率是如何变化的：

pca = PCA(n_components=X_train.shape[1])
pca.fit(X_train)
pca.explained_variance_ratio_
# 结果
array([  1.45668166e-01,   1.37354688e-01,   1.17777287e-01,
		 8.49968861e-02,   5.86018996e-02,   5.11542945e-02,
		 4.26605279e-02,   3.60119663e-02,   3.41105814e-02,
		 3.05407804e-02,   2.42337671e-02,   2.28700570e-02,
		 1.80304649e-02,   1.79346003e-02,   1.45798298e-02,
		 1.42044841e-02,   1.29961033e-02,   1.26617002e-02,
		 1.01728635e-02,   9.09314698e-03,   8.85220461e-03,
		 7.73828332e-03,   7.60516219e-03,   7.11864860e-03,
		 6.85977267e-03,   5.76411920e-03,   5.71688020e-03,
		 5.08255707e-03,   4.89020776e-03,   4.34888085e-03,
		 3.72917505e-03,   3.57755036e-03,   3.26989470e-03,
		 3.14917937e-03,   3.09269839e-03,   2.87619649e-03,
		 2.50362666e-03,   2.25417403e-03,   2.20030857e-03,
		 1.98028746e-03,   1.88195578e-03,   1.52769283e-03,
		 1.42823692e-03,   1.38003340e-03,   1.17572392e-03,
		 1.07377463e-03,   9.55152460e-04,   9.00017642e-04,
		 5.79162563e-04,   3.82793717e-04,   2.38328586e-04,
		 8.40132221e-05,   5.60545588e-05,   5.48538930e-05,
		 1.08077650e-05,   4.01354717e-06,   1.23186515e-06,
		 1.05783059e-06,   6.06659094e-07,   5.86686040e-07,
		 7.44075955e-34,   7.44075955e-34,   7.44075955e-34,
		 7.15189459e-34])

可以看到上面64个方差比率是从大到小排序的，而且后面的方差率越来越小，所以从这个数据我们其实已经可以计算出一个合适的主成分个数，使其方差比率之和达到一个极大值。我们将维数和方差率绘制出来看看：

# 横轴为维数，纵轴为每个维度对应的之前所有维度的方差率之和
plt.plot([i for i in range(X_train.shape[1])], [np.sum(pca.explained_variance_ratio_[:i+1]) for i in range(X_train.shape[1])])
plt.show()

从图中可以看到，当维度数在30左右的时候，方差率上升已经很平缓了，所以从这个图中都可以目测出，我们将64维特征的样本数据降维至30维左右是比较合适的。

其实Scikit Learn的PCA提供了一个参数，就是我们期望达到的总方差率为多少，然后会帮我们自动计算出主成分个数：

# 传入一个0到1的小数，既总方差率
pca = PCA(0.95)
pca.fit(X_train)
# 查看主成分数
pca.n_components_
# 结果
28

可以看到，我们期望的总方差率为95%时的主成分数为28。然后我们再使用KNN来训练一下降为28维特征的样本数据，看看准确率和时间为多少：

# 对原始训练数据和原始测试数据进行降维
X_train_reduction = pca.transform(X_train)
X_test_reduction = pca.transform(X_test)

X_train_reduction.shape
# 结果
(1347, 28)

# 对28维特征的数据进行KNN训练
knn_clf = KNeighborsClassifier()
%time knn_clf.fit(X_train_reduction, y_train)
# 结果
CPU times: user 2.25 ms, sys: 1.54 ms, total: 3.79 ms
Wall time: 2.44 ms

# 查看分类准确率
knn_clf.score(X_test_reduction, y_test)
# 结果
0.97999999999999998

从上面代码的结果可以看到，在使用KNN训练28维特征的数据时耗时也只有2.44毫秒，但是分类准确率达到了98%。比64维特征的数据耗时减少了15倍，但是准确率只减少了0.6%。这个性价比是非常之高的，这就是PCA的威力所在。

使用PCA对数据进行降噪

这张图是之前小节中生成的，其中蓝色的点是我们构建的原始样本数据，红色的点是经过PCA降维后，又通过PCA还原维度的样本数据。对这个图我们可以这样理解，原始样本数据的分布都在这条红色点组成的直线上下，而导致这些蓝色点没有落在红色直线上的原因就是因为数据有噪音，所以通过PCA降维其实是去除了数据的噪音。但这些噪音也是也是数据信息，所以通常我们说使用PCA对数据进行降维后会丢失一些数据信息。

下面我们通过一个实际的例子来看一下PCA的降噪过程。我们依然使用手写识别的例子，我们手写识别的样本数据中加一些噪音，然后看PCA如何去除这些噪音：

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets

digits = datasets.load_digits()
X = digits.data
y = digits.target

# 构建噪音数据，将原始样本数据增加上方差为4的噪音数据
nosiy_digits = X + np.random.normal(0, 4, size=X.shape)

# 绘制手写图片
# 取手写数字为0的前十条样本数据
example_digits = nosiy_digits[y == 0, :][:10]

# 累加循环取手写数字为1到9的前十条样本数据
for num in range(1, 10):
	num_digits = nosiy_digits[y == num, :][:10]
	example_digits = np.vstack([example_digits, num_digits])

# 目前example_digits的结构为，10条手写数字0，10条手写数字1，...， 10条手写数字9。一共100行。 
example_digits.shape
# 结果
(100, 64)

# 定义一个显示手写数字图片的方法
def plot_digits(data):
	# 构建尺寸为10 X 10的画图区域，既每行10个图片，一共10行
	plt.figure(figsize=(10, 10))
	# 遍历加了噪音的手写数字数据
	for index, imageData in enumerate(data):
		# 图片分布为10行10列，正好100个图片，第三个参数是每张图片的位置
		plt.subplot(10, 10, index + 1)
		# 将有64个元素的数组转换为8 X 8的矩阵，既8 X 8个像素的图片
		image = np.reshape(imageData, (8, 8))
		# 图片的ColorMap用灰度显示，为了能更好的观察
		plt.imshow(image, cmap = plt.cm.gray)
	
	plt.show()
	
# 显示加了噪音的手写数字
plot_digits(example_digits)

从图中可以看出，手写数字的识别度非常差。下面我们使用PCA对example_digits进行降噪处理：

from sklearn.decomposition import PCA

# 只保留50%的主成分，因为我们之前添加的噪音数据方差比较大，也就是认为噪音比较多，既有50%的噪音
pca = PCA(0.5)
pca.fit(nosiy_digits)
pca.n_components_
# 结果
12

example_components = pca.transform(example_digits)
example_components.shape
# 结果
(100, 12)

当我们只保留50%主成分的时候，特征维度从64维降到了12维。然后我们再将其还原为64维，既过滤掉了噪音：

example_components_inverse = pca.inverse_transform(example_components)

example_components_inverse.shape
# 结果
(100, 64)
# 显示降噪后的手写数字图片
plot_digits(example_components_inverse)

可以看到，此时图片中的手写数字的识别度有明显的提升。这就是PCA降噪的作用。

人脸识别

PCA有一个典型的实际应用，就是人脸识别。我们这一节就来简单看看PCA在人脸识别中的应用。

首先我们还是先从PCA的原理来说，PCA就将高维数据降至相对的低维数据，但是这些低维的数据却能反应了原始高维数据的绝大多数主要特征。那么由PCA训练出的这些主成分其实就代表了原始数据的主要特征。那么如果原始高维数据是一张张不同的人脸数据时，那么由PCA训练出的主成分其实就是这一张张人脸的主要特征，称为特征脸。好比，爷爷，爸爸，儿子三个人长的很像，其实是这三张人脸有共同的特征脸。下面我们就使用Scikit Learn中提供的人脸数据来看看特征脸的应用。

import numpy as np
import matplotlib.pyplot as plt

# 导入人脸数据集
from sklearn.datasets import fetch_lfw_people

# 导入人脸数据集(这里需要等一阵，因为需要下载数据)
faces = fetch_lfw_people()
faces.data.shape
# 结果
(13233, 2914)

faces.images.shape
# 结果
(13233, 62, 47)

如果网络不好的朋友，可以去这里手动下载All images aligned with funneling数据集，然后在使用fetch_lfw_people()时需要增加data_home这个参数，指定你存放数据集的目录，比如fetch_lfw_people(data_home='/XXX/XXX')然后执行，Scikit Learn就会去你指定的目录下解压lfw_funneled压缩包，抓取数据。从上面代码结果可以看到，这个人脸的数据集一共有13233张人脸图片，每张人脸图片有2914个特征，其实也就是由2914个不同值的像素组成，通过faces.images.shape我们知道这13233张图片的大小是62乘以47像素，相乘的值正好就是2914。

因为这个数据集里人脸的分布并不均匀，所以我们在用之前先将其进行随机处理：

# 求一个随机索引
random_indexs = np.random.permutation(len(faces.data))
X = faces.data[random_indexs]

# 取前36张脸进行展示
example_faces = X[:36, :]

# 定义一个显示人脸图片的方法
def plot_faces(data):
	# 构建尺寸为12 X 12的画图区域，既每行12个图片，一共12行
	plt.figure(figsize=(12, 12))
	# 遍历加了噪音的手写数字数据
	for index, imageData in enumerate(data):
		# 图片分布为6行6列，正好36张图片，第三个参数是每张图片的位置
		plt.subplot(6, 6, index + 1)
		# 将有2914个元素的数组转换为62 X 47的矩阵，既62 X 47个像素的图片
		image = np.reshape(imageData, (62, 47))
		plt.imshow(image, cmap = 'bone')
	
	plt.show()

然后我们通过PCA求出这个13233条样本数据，2914维特征数据集的所有主成分。所有主成分要等于小于样本数据的特征数量。下面来求这个样本数据集的最大维度：

from sklearn.decomposition import PCA
# 求出样本数据的最大主成分数
pca = PCA(svd_solver='randomized')
pca.fit(X)

pca.components_.shape
# 结果
(2914, 2914)

# 查看前36个主成分的方差比率
pca.explained_variance_ratio_[:36]

# 结果
array([ 0.23333138,  0.10693729,  0.08233336,  0.07026902,  0.03252576,
		0.03012665,  0.0210468 ,  0.01823955,  0.01764975,  0.0158582 ,
		0.01480933,  0.01337274,  0.01252589,  0.01161669,  0.01068777,
		0.00960774,  0.00873549,  0.00775879,  0.00771184,  0.00693537,
		0.00674195,  0.00626633,  0.00590671,  0.00531572,  0.00512598,
		0.00500887,  0.0048442 ,  0.00466389,  0.00435357,  0.0041365 ,
		0.00390986,  0.00387675,  0.00366964,  0.0035524 ,  0.00342629,
		0.00328901], dtype=float32)
		
# 绘制前36个主成分，既最能描述人脸特征的成分
plot_faces(pca.components_[:36, :])

从代码结果可以看出，主成分的方差比率由大到小排序，第一个主成分能表示23.33%的人脸特征，从图上也能看到，第一个图显示了一个人脸的整个脸庞，第二个图显示了人脸的五官位置等。图中的这些脸就是特征脸，可以清晰的看到特征脸依次显示人脸特征的变化过程。

申明：本文为慕课网liuyubobobo老师《Python3入门机器学习经典算法与应用》课程的学习笔记，未经允许不得转载。

今天看啥 - 高品质阅读平台
本文地址：http://www.jintiankansha.me/t/S8l5JqbAJ2