PCA降维算法（内含PCA可视化迷你案例+PCA人脸识别降维案例+PCA逆转降噪案例）_pca降维实例

1、sklearn中的降维算法

PCA使用的信息量衡量指标，就是样本方差，又称可解释性方 差，方差越大，特征所带的信息量越多。

（1）PCA

我们希望能够找出一种办法来帮助我们衡量特征上所带的信息量，让我们在降维的过程中，能够即减少特征的数量，又保留大部分有效信息——将那些带有重复信息的特征合并，并删除那些带无效信息的特征等等——逐渐创造出能够代表原特征矩阵大部分信息的，特征更少的，新特征矩阵。因此，在降维中，PCA使用的信息量衡量指标，就是样本方差，又称可解释性方差，方差越大，特征所带的信息量越多。


有空可以补充一下这个知识点



将为过程的主要步骤如下：

（2）SVD

（3）思考

PCA和特征选择技术都是特征工程的一部分，它们有什么不同？

特征工程中有三种方式：特征提取，特征创造和特征选择。仔细观察上面的降维例子和上周我们讲解过的特征 选择，你发现有什么不同了吗?

特征选择是从已存在的特征中选取携带信息最多的，选完之后的特征依然具有可解释性，我们依然知道这个特 征在原数据的哪个位置，代表着原数据上的什么含义。

而PCA，是将已存在的特征进行压缩，降维完毕后的特征不是原本的特征矩阵中的任何一个特征，而是通过某 些方式组合起来的新特征。通常来说，在新的特征矩阵生成之前，我们无法知晓PCA都建立了怎样的新特征向 量，新特征矩阵生成之后也不具有可读性，我们无法判断新特征矩阵的特征是从原数据中的什么特征组合而 来，新特征虽然带有原始数据的信息，却已经不是原数据上代表着的含义了。以PCA为代表的降维算法因此是 特征创造（feature creation，或feature construction）的一种。

可以想见，PCA一般不适用于探索特征和标签之间的关系的模型（如线性回归），因为无法解释的新特征和标 签之间的关系不具有意义。在线性回归模型中，我们使用特征选择。

2、重要参数n_components

n_components是我们降维后需要的维度，即降维后需要保留的特征数量，降维流程中第二步里需要确认的k值， 一般输入[0, min(X.shape)]范围中的整数。K是一个需要我们人为去确认的超参数，并且我们设定的数字会影响到模型的表现。就达不到降维的效果，如果留下的特征太少，那新特征向量可能无法容纳原始数据集中的大部分信息，因此，n_components既不能太大也不能太小。那怎么办呢？

可以先从我们的降维目标说起：如果我们希望可视化一组数据来观察数据分布，我们往往将数据降到三维以下，很 多时候是二维，即n_components的取值为2。

迷你实例（鸢尾花降维可视化）

from sklearn.decomposition import PCA
from matplotlib import pyplot as plt
from sklearn.datasets import load_iris
import pandas as pd
import numpy as np

iris = load_iris()
#print(iris)

x = iris.data
y = iris.target
#print(x)
#print(y)
#print(pd.DataFrame(x))
#可以了解到鸢尾花数据集一共有4个特征


pca = PCA(n_components=2) #实例化
pca = pca.fit(x)          #弥合模型
result = pca.transform(x) #获取新矩阵
#print(result)   #新矩阵只有两个特征值

x_dr = PCA(n_components=2).fit_transform(x)
#print(x_dr)    #这两种方法结果是一样的

#可视化，将三种鸢尾花的数据分布显示在二维平面坐标系中，对应的两个字符欧标应该是三种鸢尾花降维之后的坐标位置
#选取每一个种类的特征值
#print(x_dr[y == 0, 0]) #选取标签为1的第一列值，第二列也是一个意思的
plt.figure()
plt.scatter(x_dr[y == 0, 0], x_dr[y == 0, 1], color='red', label=iris.target_names[0])
plt.scatter(x_dr[y == 1, 0], x_dr[y == 1, 1], color='black', label=iris.target_names[1])
plt.scatter(x_dr[y == 2, 0], x_dr[y == 2, 1], color='orange', label=iris.target_names[2])
plt.title('pca of IRIS dataset')
plt.legend()
plt.show()

#这两种方法效果是一样的，一个不写循环，一个写循环
colors = ['red', 'black', 'yellow']
for i in range(3):
    plt.scatter(x_dr[y == i, 0], x_dr[y == i, 1], alpha=0.7, color=colors[i],label=iris.target_names[i])

plt.legend()
plt.title('PCA of IRIS dataset')
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

#探索降维后的数据
print(pca.explained_variance_ ) 
#查看降维后每个特征向量上所带的信息大小（可解释方差的大小）
#[4.22824171 0.24267075]  #就是每个新特征的方差值

print(pca.explained_variance_ratio_)  #查看降维后每个新特征向量所占的信息量占原始数据总信息量的百分比
#[0.92461872 0.05306648]

print(pca.explained_variance_ratio_.sum())
#0.9776852063187949
1
2
3
4
5
6
7
8
9
10

pca_line = PCA().fit(x)  #当不填写n_compents的时候，会计算每个特征向量的贡献值
plt.plot([1, 2, 3, 4], np.cumsum(pca_line.explained_variance_ratio_))  #np.cumsum将每个特征信息比例进行累加
plt.xticks([1, 2, 3, 4])  #这是为了限制坐标轴显示为整数
plt.xlabel('number of components after dimension reduction')
plt.ylabel('cumulative explained variance ratio')
plt.show()
1
2
3
4
5
6

最大似然估计自选超参数

pca = PCA(n_components='mle')  
#让pca进行选择降维到几维
pca_mle = pca.fit(x)
x_mle = pca_mle.transform(x)
1
2
3
4

按信息量占比选超参数

pca_f = PCA(n_components=0.97, svd_solver='full')
pca_f.fit(x)
x_f = pca_f.transform(x)
#print(x_f)

print(pca_f.explained_variance_ratio_)
#[0.92461872 0.05306648]
1
2
3
4
5
6
7

3、PCA中的SVD

重要参数svd_solver 与 random_state

参数svd_solver是在降维过程中，用来控制矩阵分解的一些细节的参数。有四种模式可选：“auto”, “full”, “arpack”,“randomized”，默认”auto"。

“auto”：基于X.shape和n_components的默认策略来选择分解器：如果输入数据的尺寸大于500x500且要提 取的特征数小于数据最小维度min(X.shape)的80％，就启用效率更高的”randomized“方法。否则，精确完整 的SVD将被计算，截断将会在矩阵被分解完成后有选择地发生。

“full”：从scipy.linalg.svd中调用标准的LAPACK分解器来生成精确完整的SVD，适合数据量比较适中，计算时 间充足的情况，生成的精确完整的SVD的结构为：

“arpack”：从scipy.sparse.linalg.svds调用ARPACK分解器来运行截断奇异值分解(SVD truncated)，分解时就将特征数量降到n_components中输入的数值k，可以加快运算速度，适合特征矩阵很大的时候，但一般用于 特征矩阵为稀疏矩阵的情况，此过程包含一定的随机性。截断后的SVD分解出的结构为：

“randomized”，通过Halko等人的随机方法进行随机SVD。在"full"方法中，分解器会根据原始数据和输入的 n_components值去计算和寻找符合需求的新特征向量，但是在"randomized"方法中，分解器会先生成多个随机向量，然后一一去检测这些随机向量中是否有任何一个符合我们的分解需求，如果符合，就保留这个随 机向量，并基于这个随机向量来构建后续的向量空间。这个方法已经被Halko等人证明，比"full"模式下计算快 很多，并且还能够保证模型运行效果。适合特征矩阵巨大，计算量庞大的情况。

重要属性components_
通常来说，在新的特征矩阵生成之前，我们无法知晓PCA都建立了怎样的新特征向量，新 特征矩阵生成之后也不具有可读性，我们无法判断新特征矩阵的特征是从原数据中的什么特征组合而来，新特征虽 然带有原始数据的信息，却已经不是原数据上代表着的含义了。

但是其实，在矩阵分解时，PCA是有目标的：在原有特征的基础上，找出能够让信息尽量聚集的新特征向量。

如果原特征矩阵是图像，V(k,n)这 个空间矩阵也可以被可视化的话，我们就可以通过两张图来比较，就可以看出新特征空间究竟从原始数据里提取了什么重要的信息。

人脸识别案例

from sklearn.decomposition import PCA
from sklearn.datasets import fetch_lfw_people
from matplotlib import pyplot as plt
import numpy as np
import pandas as pd

faces = fetch_lfw_people(min_faces_per_person=60)  #实例化
#print(faces.images.shape)  #(1348, 62, 47)
#print(faces.data.shape)  #(1348, 2914)  2914 = 62 * 47
x = faces.data #设置特征矩阵

#画图，将原特征矩阵进行可视化
fig, axes = plt.subplots(4, 5, figsize=(8, 4), subplot_kw={'xticks': [], 'yticks': []}) #不要显示坐标轴
#print(fig)
#print(axes) ##不难发现，axes中的一个对象对应fig中的一个空格
#print(axes.shape)  #(4, 5) 四行五列图片

#二维结构，可以有两种循环方式，一种是使用索引，循环一次同时生成一列上的三个图
#另一种是把数据拉成一维，循环一次只生成一个图
#在这里，究竟使用哪一种循环方式，是要看我们要画的图的信息，储存在一个怎样的结构里
#我们使用 子图对象.imshow 来将图像填充到空白画布上
#而imshow要求的数据格式必须是一个(m,n)格式的矩阵，即每个数据都是一张单独的图
#因此我们需要遍历的是faces.images，其结构是(1277, 62, 47)
#要从一个数据集中取出24个图，明显是一次性的循环切片[i,:,:]来得便利
#因此我们要把axes的结构拉成一维来循环

#print(enumerate(axes.flat))
#填充图像
for i, ax in enumerate(axes.flat):
    ax.imshow(faces.images[i, :, :], cmap='gray')#选择色彩模式
#plt.show()
#print(fig)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

#建模降维，提取新特征空间矩阵
pca = PCA(150).fit(x) #将data原来的2914维降到150维
v = pca.components_
#print(v.shape)  #(150, 2914) （新特征维度，原特征维度）

#print(v[0, :].reshape(62, 47))
'''
[[-0.00579708 -0.00595363 -0.00615767 ... -0.00813068 -0.00789127
  -0.00791125]
 [-0.00586497 -0.00604798 -0.00647989 ... -0.00849601 -0.00826257
  -0.00840366]
 [-0.00622997 -0.00658633 -0.00714348 ... -0.00897475 -0.0085009
  -0.0089165 ]
 ...
 [-0.00270353 -0.00276129 -0.00370097 ... -0.01046781 -0.00896865
  -0.00820496]
 [-0.00298345 -0.003292   -0.00394306 ... -0.01013859 -0.00875235
  -0.00805416]
 [-0.00329591 -0.00362476 -0.00457148 ... -0.01000113 -0.00901094
  -0.00813919]]
'''

#print(v[0, :]) #取降维后的像素，新的特征将轴进行了变换，我们想看看变换后的图像是什么样的
'''
[-0.00579718 -0.00595365 -0.00615768 ... -0.01000111 -0.00901092
 -0.00813918]
'''

#将新特征空间矩阵可视化
fig, axes = plt.subplots(3, 8, figsize=(8, 4), subplot_kw={'xticks': [], 'yticks': []})

for i, ax in enumerate(axes.flat):
    ax.imshow(v[i, :].reshape(62, 47), cmap='gray')
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

提取的主要特征有光线和眼睛

4、重要接口inverse_transform

from sklearn.datasets import fetch_lfw_people
from sklearn.decomposition import PCA
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt

faces = fetch_lfw_people(min_faces_per_person=60)
x = faces.data
#print(x)
#print(x.shape) #(1348, 2914)

pca = PCA(150)
x_dr = pca.fit_transform(x)

#print(x_dr)
#print(x_dr.shape)

#print(x_dr[0])
x_inverse = pca.inverse_transform(x_dr)
#print(x_inverse)

fig, ax = plt.subplots(2, 10, figsize=(10, 2.5), subplot_kw={"xticks": [], "yticks": []})
#plt.show()
#print(ax)

for i in range(10):
    ax[0, i].imshow(faces.images[i, :, :], cmap='binary_r')
    ax[1, i].imshow(x_inverse[i].reshape(62, 47), cmap='binary_r')

plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

逆转换之后的图像并没有原来的图像清洗，是因为在经过PVC降维之后去掉了一些噪音，只提取了重要的特征，所以在逆转的过程中并不能逆转已经去掉的信息。

案例：用PCA做噪音过滤

from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt

digits = load_digits()
#print(digits)

#print(digits.data.shape)#(1797, 64)
#print(digits.images.shape) #(1797, 8, 8)
#print(set(digits.target.tolist())) #{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}
#是一个手写识别分类数据，标签是从0-9


def plot_digits(data):
    fig, axes = plt.subplots(4, 10, figsize=(10, 4), subplot_kw={'xticks': [], 'yticks': []})
    for i, ax in enumerate(axes.flat):
        ax.imshow(data[i].reshape(8, 8), cmap="binary")
    plt.show()

plot_digits(digits.data)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

给手写数字加噪音，取出一部分服从正太分布的、方差为2的数据

rng = np.random.RandomState(42)
#给原始数据加上噪音

digits_noise = rng.normal(digits.data, 2)
#从原始数据集中取出来一个服从正太分布，方差为2的新数据集
#plot_digits(digits_noise)
1
2
3
4
5
6

加了噪音之后非常模糊

利用pca提取使得特征解释率在0.5上的特征即可，经过主要特征的特区进行逆转，逆转之后的图片只体现主要的信息，会稍微清晰一些

pca = PCA(0.5, svd_solver='full').fit(digits_noise)
x_dr = pca.transform(digits_noise)
print(x_dr.shape)  #(1797, 6)

x_inverse = pca.inverse_transform(x_dr)
#plot_digits(x_inverse)
1
2
3
4
5
6

5、重要接口，参数和属性总结

6、案例：PCA对手写数字数据集的降维

from sklearn.decomposition import PCA
from sklearn.ensemble import RandomForestClassifier as RFC
from sklearn.neighbors import KNeighborsClassifier as KNN
from sklearn.model_selection import cross_val_score
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import time

data = pd.read_csv('digit recognizor.csv')
#print(data.shape)
#print(data.head())
#print(data.info())

x = data.iloc[:, 1:] #(42000, 785)
y = data.iloc[:, 0]

pca_line = PCA().fit(x)
plt.figure(figsize=(20, 5))
plt.plot(np.cumsum(pca_line.explained_variance_ratio_)) #将维度解释率累加
plt.xlabel('number of components after dimension reduction')
plt.ylabel('cumulative explained variance ratio')
#plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

先通过解释率累加图来寻找一个合理的维度，因为785维实属太高维了！

我们一般拐点处的维度，大概看起来就是在1-101之间

scores = []
#继续缩小最佳维度的范围
for i in range(1, 101, 10):
    x_dr = PCA(i).fit_transform(x)
    once = cross_val_score(RFC(n_estimators=10, random_state=0), x_dr, y, cv=5).mean()
    scores.append(once)
plt.figure(figsize=(10, 8))
plt.plot(range(1, 101, 10), scores)
plt.show()
1
2
3
4
5
6
7
8
9

再次缩小范围

start_time = time.time()
scores = []
#继续缩小最佳维度的范围
for i in range(15, 30):
    x_dr = PCA(i).fit_transform(x)
    once = cross_val_score(RFC(n_estimators=10, random_state=0), x_dr, y, cv=5).mean()
    scores.append(once)
plt.figure(figsize=(10, 8))
plt.plot(range(15, 30), scores)
end_time = time.time()
print(end_time-start_time)
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12

于是，取维度26作为我们降维之后的降维之后的新特征维度，并输入随机森林中

x_dr = PCA(26).fit_transform(x)
score = cross_val_score(RFC(n_estimators=100, random_state=0), x_dr, y, cv=5).mean()
print(score)  #0.9466190476190477
1
2
3

此时，我们想到了KNN比随机森林更适合跑这个数据集，所以我们换一个算法再次验证

x_dr = PCA(26).fit_transform(x)
score_knn = cross_val_score(KNN(), x_dr, y, cv=5).mean()
print(score_knn)  #0.9710238095238095
1
2
3

因为这时的KNN（）里面是默认值，所以我们可以寻找一个更好的参数

x_dr = PCA(26).fit_transform(x)
score_1 = []
for i in range(10):
    once = cross_val_score(KNN(i+1), x_dr, y, cv=5).mean()
    score_1.append(once)
plt.figure(figsize=(10, 8))
plt.plot(range(10), score_1)
plt.show()
1
2
3
4
5
6
7
8

此时，KNN里面的参数取2和4都试试

x_dr = PCA(26).fit_transform(x)
score_knn = cross_val_score(KNN(3), x_dr, y, cv=5).mean()
print(score_knn)  #0.9708571428571429

x_dr = PCA(26).fit_transform(x)
score_r = cross_val_score(KNN(5), x_dr, y, cv=5).mean()
print(score_knn)  #0.9708571428571429
1
2
3
4
5
6
7

发现结果都不如默认的参数，那就还是使用默认参数。