sklean实战04：降维算法PCA和SVD

1 PCA与SVD

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

Var代表一个特征的方差，n代表样本量，xi代表一个特征中的每个样本取值，xhat代表这一列样本的均值

1.1 sklearn.decomposition.PCA

class sklearn.decomposition.PCA (n_components=None
								, copy=True
								, whiten=False
								, svd_solver=’auto’
								, tol=0.0,
								iterated_power=’auto’
								, random_state=None)

在步骤3当中，我们用来找出n个新特征向量，让数据能够被压缩到少数特征上并且总信息量不损失太多的技术就是矩阵分解。PCA和SVD是两种不同的降维算法.
PCA使用方差作为信息量的衡量指标，并且特征值分解来找出空间V。降维完成之后，PCA找到的每个新特征向量就叫做“主成分”，而被丢弃的特征向量被认为信息量很少，这些信息很可能就是噪音。Q 和Q-1是辅助的矩阵，Σ是一个对角矩阵（即除了对角线上有值，其他位置都是0的矩阵），其对角线上的元素就是方差

而SVD使用奇异值分解来找出空间V，其中Σ也是一个对角矩阵，不过它对角线上的元素是奇异值，这也是SVD中用来衡量特征上的信息量的指标。U和V^{T}分别是左奇异矩阵和右奇异矩阵，也都是辅助矩阵

1.2 重要参数n_components

n_components是我们降维后需要的维度,一般输入[0, min(X.shape)]范围中的整数

1.2.1 案例：高维数据的可视化

import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
import pandas as pd
import numpy as np
iris = load_iris()
y = iris.target
X = iris.data
#建模
X_dr = PCA(2).fit_transform(X)
#画图
plt.figure()
color = ['red', 'black', 'orange']
for i in [0,1,2]:
    plt.scatter(X_dr[y==i,0]
             ,X_dr[y==i,1]
             ,alpha=0.7
             ,color = color[i]
             ,label=iris.target_names[i])
plt.legend()
plt.show()

#探索降维后的数据
#属性explained_variance_，查看降维后每个新特征向量上所带的信息量大小（可解释性方差的大小）
pca = PCA(2).fit(X)
pca.explained_variance_ #array([4.22824171, 0.24267075])

属性explained_variance_ratio，查看降维后每个新特征向量所占的信息量占原始数据总信息量的百分比
#又叫做可解释方差贡献率,#大部分信息都被有效地集中在了第一个特征上
pca.explained_variance_ratio_  #array([0.92461872, 0.05306648])
pca.explained_variance_ratio_.sum() 

选择最好的n_components：累积可解释方差贡献率曲线
当参数n_components中不填写任何值，则默认返回min(X.shape)个特征，一般来说，样本量都会大于特征数目，所以什么都不填就相当于转换了新特征空间，但没有减少特征的个数。一般来说，不会使用这种输入方式。但我们却可以使用这种输入方式来画出累计可解释方差贡献率曲线，以此选择最好的n_components的整数取值。累积可解释方差贡献率曲线是一条以降维后保留的特征个数为横坐标，降维后新特征矩阵捕捉到的可解释方差贡献率为纵坐标的曲线，能够帮助我们决定n_components最好的取值。

#选择最好的n_components：累积可解释方差贡献率曲线
pca_line = PCA().fit(X)
plt.figure()
plt.plot([1,2,3,4],np.cumsum(pca_line.explained_variance_ratio_))
plt.xticks([1,2,3,4])

从图中可以看出当特征选2个或3个时，信息保留的量大于98%，因此在这个数据集上n_components可以选2或3。当特征很多的时候，我们找曲线突然转折的那个点。

1.2.2 最大似然估计自选超参数

让PCA用最大似然估计(maximum likelihood estimation)自选超参数的方法，输入“mle”作为n_components的参数输入

pca_mle = PCA(n_components="mle")
pca_mle = pca_mle.fit(X)
X_mle = pca_mle.transform(X)
X_mle 
pca_mle.explained_variance_ratio_.sum() # 0.9947878161267246

1.2.3 按信息量占比选超参数

输入[0,1]之间的浮点数，并且让参数svd_solver ==‘full’，表示希望降维后的总解释性方差占比大于n_components指定的百分比，即是说，希望保留百分之多少的信息量。比如说，如果我们希望保留97%的信息量，就可以输入n_components = 0.97，PCA会自动选出能够让保留的信息量超过97%的特征数量。

pca_f = PCA(n_components=0.97,svd_solver="full")
pca_f = pca_f.fit(X)
X_f = pca_f.transform(X)
pca_f.explained_variance_ratio_

实际中可以先用信息量占比尝试一下，推荐0.8，看看模型什么效果

1.3 PCA中的SVD

1.3.1 PCA中的SVD哪里来？

sklearn将降维流程拆成了两部分：一部分是计算特征空间V，由奇异值分解完成，另一部分是映射数据和求解新特征矩阵，由主成分分析完成，实现了用SVD的性质减少计算量，却让信息量的评估指标是方差，具体流程如下图：

在transform过程之后，fit中奇异值分解的结果除了V(k,n)以外，就会被舍弃，而V(k,n)会被保存在属性components_ 当中，可以调用查看。

PCA(2).fit(X).components_
PCA(2).fit(X).components_.shape

1.3.2 重要参数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，可以加快运算速度，适合特征矩阵很大的时候，但一般用于特征矩阵为稀疏矩阵的情况，此过程包含一定的随机性。
"randomized"：通过Halko等人的随机方法进行随机SVD。在"full"方法中，分解器会根据原始数据和输入的n_components值去计算和寻找符合需求的新特征向量，但是在"randomized"方法中，分解器会先生成多个随机向量，然后一一去检测这些随机向量中是否有任何一个符合我们的分解需求，如果符合，就保留这个随机向量，并基于这个随机向量来构建后续的向量空间。这个方法已经被Halko等人证明，比"full"模式下计算快很多，并且还能够保证模型运行效果。适合特征矩阵巨大，计算量庞大的情况。

而参数random_state在参数svd_solver的值为"arpack" or "randomized"的时候生效，可以控制这两种SVD模式中
的随机模式。通常我们就选用”auto

1.3.3 重要属性components_

在矩阵分解时，PCA是有目标的：在原有特征的基础上，找出能够让信息尽量聚集的新特征向量。在sklearn使用的PCA和SVD联合的降维方法中，这些新特征向量组成的新特征空间其实就是V(k,n)。当V(k,n)是数字时，我们无法判断V(k,n)和原有的特征究竟有着怎样千丝万缕的数学联系。但是，如果原特征矩阵是图像，V(k,n)这个空间矩阵也可以被可视化

faces = fetch_lfw_people(min_faces_per_person=60)
face.images.shape #(1277,62,47):1277是矩阵中图形的个数，62是每个图像的特征矩阵的行，47是每个图形的特征矩阵的列
face.data.shape # (1277,2914):1277是矩阵中图形的个数,2914=62*47,因为fit只接受二维数组
X = face.data
#看看图像什么样？将原特征矩阵进行可视化
fig,axes = plt.subplots(4,5
                       ,figsize=(8,5)
                       ,subplot_kw = {'xticks':[],'yticks':[]})# 不要x轴和y轴坐标
#填充图像
for i, ax in enumerate(axes.flat): #enumerate(axes.flat)等价于enumerate(axes.ravel())
    ax.imshow(faces.images[i,:,:]  #而imshow要求的数据格式必须是一个(m,n)格式的矩阵，即每个数据都是一张单独的图,因此我们需要遍历的faces.images，其结构是(1277, 62, 47),要从一个数据集中取出24个图，明显是一次性的循环切片[i,:,:]来得遍历
             ,cmap="gray" #选择色彩的模式)
#建模降维，提取新特征空间矩阵
#原本有2900维，我们现在来降到150维
pca = PCA(150).fit(X) 
V = pca.components_    
V.shape # (150,2914)
#将新特征空间矩阵可视化
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")

1.3.4 重要接口inverse_transform

在sklearn中，我们通过让原特征矩阵X右乘新特征空间矩阵V(k,n)来生成新特征矩阵X_dr，那理论上来说，让新特征矩阵X_dr右乘V(k,n)的逆矩阵 ，就可以将新特征矩阵X_dr还原为X。那sklearn是否这样做了呢？让我们来看看下面的案例。

1.3.4.1 案例：用人脸识别看PCA降维后的信息保存量

from sklearn.datasets import fetch_lfw_people
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
import numpy as np
#导入数据，探索数据
faces = fetch_lfw_people(min_faces_per_person=60)
faces.images.shape
#怎样理解这个数据的维度？
faces.data.shape
#换成特征矩阵之后，这个矩阵是什么样？
X = faces.data
#建模，获取降维后的特征矩阵X_dr
pca = PCA(150)
X_dr = pca.fit_transform(X)
X_dr.shape #（1277，150）
#将降维后的矩阵用inverse_transform返回原空间
X_inverse = pca.inverse_transform(X_dr)
X_inverse.shape #（1277，2914）
#将特征矩阵X和X_inverse可视化
fig, ax = plt.subplots(2,10,figsize=(10,2.5),subplot_kw={"xticks":[],"yticks":[]})
#和2.3.3节中的案例一样，我们需要对子图对象进行遍历的循环，来将图像填入子图中
#那在这里，我们使用怎样的循环？
#现在我们的ax中是2行10列，第一行是原数据，第二行是inverse_transform后返回的数据
#所以我们需要同时循环两份数据，即一次循环画一列上的两张图，而不是把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") #第二行画inverse后的数据

可以明显看出，这两组数据可视化后，由降维后再通过inverse_transform转换回原维度的数据画出的图像和原数据画的图像大致相似，但原数据的图像明显更加清晰。这说明，inverse_transform并没有实现数据的完全逆转。
这是因为，在降维的时候，部分信息已经被舍弃了，X_dr中往往不会包含原数据100%的信息，所以在逆转的时候，即便维度升高，原数据中已经被舍弃的信息也不可能再回来了。所以，降维不是完全可逆的。
Inverse_transform的功能，是基于X_dr中的数据进行升维，将数据重新映射到原数据所在的特征空间中，而并非恢复所有原有的数据。但同时，我们也可以看出，降维到300以后的数据，的确保留了原数据的大部分信息，所以
图像看起来，才会和原数据高度相似，只是稍稍模糊罢了。

1.3.4.2 用PCA做噪音过滤

降维的目的之一就是希望抛弃掉对模型带来负面影响的特征，而我们相信，带有效信息的特征的方差应该是远大于噪音的，所以相比噪音，有效的特征所带的信息应该不会在PCA过程中被大量抛弃。inverse_transform能够在不恢复原始数据的情况下，将降维后的数据返回到原本的高维空间，即是说能够实现”保证维度，但去掉方差很小特征所带的信息“。利用inverse_transform的这个性质，我们能够实现噪音过滤。

digits = load_digits()
digits.data.shape
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")
plot_digits(digits.data)
#为数据加上噪音
rng = np.random.RandomState(42) 
#在指定的数据集中，随机抽取服从正态分布的数据
#两个参数，分别是指定的数据集，和抽取出来的正太分布的方差
noisy = rng.normal(digits.data,2)
plot_digits(noisy)
pca = PCA(0.5).fit(noisy) #数据特征保留50%
X_dr = pca.transform(noisy)
X_dr.shape
#逆转降维结果，实现降噪
without_noise = pca.inverse_transform(X_dr)
plot_digits(without_noise)

2 PCA对手写数字数据集的降维

data = pd.read_csv("digit recognizor.csv") 
X = data.iloc[:,1:]
y = data.iloc[:,0] 
X.shape
#画累计方差贡献率曲线，找到最佳降维后维度的范围
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()
#继续缩小范围
score = []
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()
    score.append(once)
plt.figure(figsize=[20,5])
plt.plot(range(1,101,10),score)
plt.show()
#细化学习曲线，找出降维后的最佳维度
score = []
for i in range(10,35):
    x_dr = PCA(i).fit_transform(X)
    once = cross_val_score(RFC(n_estimators=10,random_state=0),x_dr,y,cv=5).mean()
    score.append(once)
plt.figure(figsize=[20,5])
plt.plot(range(10,35),score)
plt.show()
#导入找出的最佳维度进行降维，查看模型效果
X_dr = PCA(26).fit_transform(X)
cross_val_score(RFC(n_estimators=100,random_state=0),X_dr,y,cv=5).mean()
# 突发奇想，特征数量已经不足原来的3%，换模型怎么样？
cross_val_score(KNN(),X_dr,y,cv=5).mean()
#细化学习曲线，找出降维后的最佳维度
score = []
for i in range(10):
    once = cross_val_score(KNN(i+1),x_dr,y,cv=5).mean()
    score.append(once)
plt.figure(figsize=[20,5])
plt.plot(range(10),score)
plt.show()