凌云时刻 · 技术
导读： 这篇笔记将继续讲解机器学习中经常用到的降维算法PCA及PCA降噪的作用。
作者 | 计缘
来源 | 凌云时刻（微信号：linuxpk）
高维数据向低维数据映射
我们再来回顾一下PCA降维的基本原理，首先要做的事情就是对样本数据寻找另外一个坐标系，这个坐标系中的每一个轴依次可以表达样本数据的重要程度，既主成分。我们取出前k个主成分，然后就可以将所有的样本数据映射到这k个轴上，既获得了一个低维的数据信息。
上面的
是样本数据，该样本数据有m行，n个特征，既是一个n维的样本数据。
假设上面的
是样本数据
的主成分向量矩阵，每一行代表一个主成分向量，一共有k个主成分向量，每个主成分向量上有n个值。
我们已经推导出了求映射后的向量的大小，也就是每一行样本数据映射到该主成分上的大小为：
那如果将一行有n个特征的样本数据分别映射到k个主成分上，既得到有k个值的新向量，既降维后的，有k个特征的新样本数据。所以我们需要的就是
矩阵的第一行和
矩阵的每一行对应元素相乘然后再相加，
矩阵的第二行和
矩阵的每一行对应元素相乘然后再相加，以此类推就可以求出降维后的，m行k列的新矩阵数据：
就是降维后的数据，既然可以降维，那么我们也可从数学的角度将降维后的数据还原回去。
是m行k列的矩阵，
是k行n列的矩阵，所以
就是还原后的ｍ行ｎ列的原矩阵。那为什么说是从数学角度来说呢，因为毕竟已经从高维降到了低维，那势必会有丢失的数据信息，所以还原回去的数据也不可能和原始数据一样的。
我们在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
首先构建样本数据：
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()
其实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的威力所在。
Scikit Learn中的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降噪的作用。
END
往期精彩文章回顾
机器学习笔记（十三）：主成分分析法（PCA）
机器学习笔记（十二）：随机梯度下降
机器学习笔记（十）：梯度下降
机器学习笔记（九）：多元线性回归
机器学习笔记（七）：线性回归
机器学习笔记（六）：数据归一化
机器学习笔记（五）：超参数
机器学习笔记（四）：kNN算法
长按扫描二维码关注凌云时刻
每日收获前沿技术与科技洞见
在PyCharm中封装PCA
2
0
已为社区贡献12条内容
所有评论(0)