【4 - 降维算法PCA和SVD - 原理部分】菜菜sklearn机器学习_菜菜 线性回归大家族

课程地址：《菜菜的机器学习sklearn课堂》_哔哩哔哩_bilibili

• 第一期：sklearn入门 & 决策树在sklearn中的实现
• 第二期：随机森林在sklearn中的实现
• 第三期：sklearn中的数据预处理和特征工程
• 第四期：sklearn中的降维算法PCA和SVD
• 第五期：sklearn中的逻辑回归
• 第六期：sklearn中的聚类算法K-Means
• 第七期：sklearn中的支持向量机SVM（上）
• 第八期：sklearn中的支持向量机SVM（下）
• 第九期：sklearn中的线性回归大家族
• 第十期：sklearn中的朴素贝叶斯
• 第十一期：sklearn与XGBoost
• 第十二期：sklearn中的神经网络

---

目录

概述

（一）维度

数组 vs 特征矩阵

（二）sklearn中的降维算法（模块decomposition）

PCA与SVD

（一）概述

1. 降维的实现

2. PCA & SVD

3. PCA和特征选择的异同 

（二）PCA 

1. 重要参数 n_components（降维后需要的维度/保留的特征数量）

2. 鸢尾花数据集（高维数据）的可视化

（三）PCA中的SVD

1. 为什么PCA的类里会包含控制SVD分解器的参数？

2. 重要参数 svd_solver 与 random_state 

3. 重要属性 components_

4. 重要接口inverse_transform

---

前三周：

• 随机森林：通过随机抽取特征来建树，以避免高维计算
• sklearn中导入特征矩阵必须至少二维
• 特征工程：特征选择的目的是通过降维来降低算法的计算成本

本周： 

• PCA：主成分分析
• SVD：奇异值分解

---

概述

（一）维度

对于数组和Series来说，维度就是shape返回的结果，shape中返回了几个数字就是几维

• 一维（不分行列）：此时shape返回唯一的维度上的数据个数
• 二维（有行列之分，表）：shape返回几行几列，一张表最多二维
• 三维：当一个数组中存在2张3行4列的表时，shape返回(2,3,4)
• 四维：当数组中存在2组2张3行4列的表时，shape返回(2,2,3,4)
• 更高维：shape返回(更高维,行,列)

数组 vs 特征矩阵

• 高维矩阵 = 数组，数组中的每一张表，都可以是一个特征矩阵 
• 特征矩阵特指二维，只有行和列，其中行是样本，列是特征。针对每一张表，维度指的是特征的数量

对图像来说，维度就是图像中特征向量的数量（特征向量可以理解为坐标轴）

• 一个特征向量定义一条直线，是一维
• 两个相互垂直的特征向量定义一个平面，即一个直角坐标系，是二维
• 三个相互垂直的特征向量定义一个空间，即一个立体直角坐标系，是三维
• 三个以上的特征向量相互垂直，定义人眼无法看见也无法想象的高维空间

降维算法中的降维，指的是降低特征矩阵中特征的数量

• 让算法运算更快、效果更好
• 数据可视化：三维及以下的特征矩阵是可以被可视化的，便于理解数据分布

（二）sklearn中的降维算法（模块decomposition）

decomposition：本质是一个矩阵分解模块

• 可以用在降维、深度学习、聚类分析、数据预处理、低纬度特征学习、推荐系统、大数据分析等领域
• 奇异值分解SVD和主成分分析PCA都属于矩阵分解算法中的入门算法，都是通过分解特征矩阵来降维

---

PCA与SVD

在高维数据中，必然有一些特征是不带有有效信息的（比如噪音），或有一些特征带有的信息和其他一些特征是重复的（比如一些特征可能会线性相关）

衡量特征上所带的信息量，让降维过程中能够既减少特征的数量，又保留大部分有效信息

（一）概述

如果一个特征的方差很大，则说明这个特征上带有大量信息

PCA使用的信息量衡量指标 —— 样本方差/可解释性方差（方差越大，特征所带的信息量越多）

方差计算公式中为什么除数是n-1？

为了得到样本方差的无偏估计 

1. 降维的实现

例子：

• 图1：一组数据有两个特征x1、x2（分别作为两个特征向量），这组数据每个特征的均值都为2，方差都为1，数据的方差总和为2

目标：只用一个特征向量来描述这组数据，即将二维数据降为一维数据，并尽可能保留信息量（信息量用样本方差衡量），即让数据总方差接近2

• 图2：x2*上的数值都变成了0，故x2*不带有任何有效信息了（x2*的均值和方差都为0）；而x1*特征上的数据均值是，方差为2
• 图3：根据信息含量的排序，取信息含量最大的一个特征

通过旋转原有特征向量组成的坐标轴来找到新特征向量和新坐标平面，将三个样本点的信息压缩到了一条直线上，实现了二维变一维，并尽量保留原始数据的信息

在步骤3中，找出n个新特征向量，让数据能够被压缩到少数特征上并且总信息量不损失太多的技术——矩阵分解 

2. PCA & SVD

PCA和SVD是两种不同的降维算法，但都遵循上面的五步来实现降维，只是两种算法中矩阵分解的方法（步骤3）不同，且信息量的衡量指标不同

无论是PCA还是SVD都需要遍历所有的特征和样本来计算信息量指标，并且在矩阵分解过程中会产生比原来的特征矩阵更大的矩阵，还需要产生协方差矩阵去计算更多的信息，故降维算法的计算量很大，运行比较缓慢

（1）PCA：

• 使用方差作为信息量的衡量指标
• 特征值分解来找出空间V

（2）SVD：

• 使用奇异值来衡量特征上的信息量
• 奇异值分解来找出空间V

3. PCA和特征选择的异同 

（1）同：都是特征工程的一部分

特征工程有三种方式：特征提取、特征创造、特征选择

（2）异：以PCA为代表的降维算法是特征创造的一种

• 特征选择：从已存在的特征中选取携带信息最多的，选完之后的特征依然具有可解释性，依然知道这个特征在原数据的哪个位置，代表着原数据上的什么含义
• PCA：将已存在的特征进行压缩，降维后的特征不是原本特征矩阵中的任何一个特征，而且通过某些方式组合起来的新特征。在新的特征矩阵生成之前，无法知晓PCA都建立了怎样的新特征向量，新特征向量生成之后也不具有可读性（无法判断新特征矩阵的特征是从原数据中的什么特征组合而来，新特征虽然带有原始数据的信息，却已经不是原数据上代表的含义了）

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

（二）PCA 

 sklearn.decomposition.PCA — scikit-learn 1.2.0 documentation

1. 重要参数 n_components（降维后需要的维度/保留的特征数量）

降维流程中第二步里需要确认的k一般输入 [0, min(X.shape)] 中的整数，是一个超参数，值会影响模型表现

• 若留下的特征太多，就达不到降维的效果；
• 若留下的特征太少，新特征向量可能无法容纳原始数据集中的大部分信息

如果希望可视化一组数据来观察数据分布，往往将数据降到三维以下，很多时候是二维（即n_components = 2）

如何选择最好的n_components  —— 累计可解释方差贡献率曲线

累计可解释方差贡献率曲线：一条以降维后保留的特征个数为横坐标，降维后新特征矩阵捕捉到的可解释方差贡献率为纵坐标的曲线

当参数n_components中不填写任何值，则默认返回 min(X.shape) 个特征

一般来说，样本量都会大于特征数目，所以什么都不填就相当于：转换了新特征空间，但没有减少特征的个数

可以使用这种输入方式来画累计可解释方差贡献率曲线，以此选择最好的 n_components 整数取值（选曲线转折点的那个特征数量）

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
 
pca_line = PCA().fit(X)  # 不填参数，默认为特征个数
pca_line.explained_variance_ratio_   # array([0.92461872, 0.05306648, 0.01710261, 0.00521218])  
# .sum()为1
np.cumsum(pca_line.explained_variance_ratio_ )    # array([0.92461872, 0.97768521, 0.99478782, 1.        ])
# np.cumsum()累加
 
# plt.plot(横坐标取值, 纵坐标)
plt.plot([1,2,3,4],np.cumsum(pca_line.explained_variance_ratio_))  
plt.xticks([1,2,3,4])   #这是为了限制坐标轴显示为整数
plt.xlabel("number of components after dimension reduction")
plt.ylabel("cumulative explained variance ratio")
plt.show()

（1）最大似然估计自选超参数

让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  #3列的数组
#可以发现，mle为我们自动选择了3个特征

pca_mle.explained_variance_ratio_.sum()   #0.9947878161267246
#得到了比设定2个特征时更高的信息含量，对于鸢尾花这个很小的数据集来说，3个特征对应这么高的信息含量，并不需要去纠结于只保留2个特征，毕竟三个特征也可以可视化

（2）按信息量占比选超参数

输入 [0,1] 之间的浮点数，并让参数 svd_solver = "full" ，表示希望降维后的总解释性方差占比＞ n_components 指定的百分比，即希望保留百分之多少的信息量 

svd_solver 是奇异值分解器

# 希望保留97%的信息量
pca_f = PCA(n_components=0.97,svd_solver="full")  # svd_solver="full"不能省略
pca_f = pca_f.fit(X)
X_f = pca_f.transform(X)
X_f   # PCA会自动选出能够让保留的信息量超过97%的特征数量

pca_f.explained_variance_ratio_   #array([0.92461872, 0.05306648])

2. 鸢尾花数据集（高维数据）的可视化

import matplotlib.pyplot as plt
from sklearn.datasets import load_iris   # 鸢尾花数据集
from sklearn.decomposition import PCA
 
iris = load_iris()
y = iris.target
X = iris.data   # 四个特征
 
# 作为数组，X是几维？   二维数组(150,4)，150个样本数量，4个特征
X.shape  #(150, 4)
 
# 作为数据表或特征矩阵，X是几维？   四维特征矩阵（4列特征）
import pandas as pd
pd.DataFrame(X).head()
 
# 数组的维度和矩阵的维度不同
 
# 建模，调用PCA
pca = PCA(n_components=2)           #实例化，降到二维
pca = pca.fit(X)                    #拟合模型，输入特征矩阵
X_dr = pca.transform(X)             #获取新矩阵
 
X_dr
 
#也可以fit_transform一步到位
#X_dr = PCA(2).fit_transform(X)

#要将三种鸢尾花的数据分布显示在二维平面坐标系中，对应的两个坐标（两个特征向量）应该是三种鸢尾花降维后的x1和x2，怎样才能取出三种鸢尾花下不同的x1和x2呢？
 
X_dr[y == 0, 0] #这里是布尔索引（返回为True的行），取出标签为0的鸢尾花对应的第一个特征的数据 

#要展示三种分类的分布，需要对三种鸢尾花分别绘图
#可以写成三行代码，也可以写成for循环
"""
plt.figure()  
plt.scatter(X_dr[y==0, 0], X_dr[y==0, 1], c="red", label=iris.target_names[0])
plt.scatter(X_dr[y==1, 0], X_dr[y==1, 1], c="black", label=iris.target_names[1])
plt.scatter(X_dr[y==2, 0], X_dr[y==2, 1], c="orange", label=iris.target_names[2])
plt.legend()
plt.title('PCA of IRIS dataset')
plt.show()
"""
 
colors = ['red', 'black', 'orange']
iris.target_names   # array(['setosa', 'versicolor', 'virginica'], dtype='<U10')
 
plt.figure()  # 画布
for i in [0, 1, 2]:
    # plt.scatter(第一个坐标轴上的值, 第二个坐标轴上的值, 颜色, label=这一组点的名称)画散点图
    plt.scatter(X_dr[y == i, 0]
                ,X_dr[y == i, 1]
                ,alpha=.7   # 点的透明度
                ,c=colors[i]
                ,label=iris.target_names[i]
               )
plt.legend()  #图例
plt.title('PCA of IRIS dataset')
plt.show()

探索降维后的数据：

#属性explained_variance_，查看降维后每个新特征向量上所带的信息量大小（可解释性方差的大小）
pca.explained_variance_   #查看方差是否从大到小排列，第一个最大，依次减小   array([4.22824171, 0.24267075])
 
#属性explained_variance_ratio，查看降维后每个新特征向量所占的信息量占原始数据总信息量的百分比，又叫做可解释方差贡献率
pca.explained_variance_ratio_   #array([0.92461872, 0.05306648])
#大部分信息都被有效地集中在了第一个特征上
 
# .sum()查看降维后的特征向量所带的所有信息在原始总信息量上的占比
pca.explained_variance_ratio_.sum()   #0.9776852063187949

（三）PCA中的SVD

1. 为什么PCA的类里会包含控制SVD分解器的参数？

Answer：把SVD当做PCA的一种求解方法，即在矩阵分解时不使用PCA本身的特征值分解，而是用SVD的奇异值分解来减少计算量

PCA和SVD涉及了大量的矩阵计算，但SVD可以跳过数学神秘的宇宙，不计算协方差矩阵，直接找出一个新特征向量组成的n维空间，而这个n维空间就是奇异值分解后的右矩阵

若X为(m,n)特征矩阵， 就是(n,n)，取其前k行（切片），即将V转换为(k,n)，原特征矩阵X与  相乘，即可得到降维后的特征矩阵（m×k）

奇异值分解可以不计算协方差矩阵等结构复杂、计算冗长的矩阵，就直接求出新特征空间和降维后的特征矩阵，故SVD在矩阵分解中的过程比PCA简单快速。但SVD的信息量衡量指标（奇异值）比PCA的（方差）复杂

因此，sklearn将降维流程拆成了两部分，实现了用SVD的性质减少计算量，却让信息量的评估指标是方差：

• 一部分是计算新的特征空间V，由奇异值分解SVD完成
• 另一部分是映射数据和求解新特征矩阵，由主成分分析PCA完成

PCA(2).fit(X).components_  # 2指 n_components=2
# array([[ 0.36138659, -0.08452251,  0.85667061,  0.3582892 ],
#        [ 0.65658877,  0.73016143, -0.17337266, -0.07548102]])
 
# X.shape()是(m,n) = (150,4)
PCA(2).fit(X).components_.shape  #(2, 4)  V(k,n)即切片后/降维后的新特征空间

2. 重要参数 svd_solver 与 random_state 

参数 svd_solver 是在降维过程中用来控制矩阵分解的一些细节的参数，有四种模式：

• "auto"：默认。数据量小时为full（精确完整的SVD将被计算，截断将会在矩阵被分解完成后有选择地发生），数据量大时为randomized
• "full"：生成精确完整的SVD，适合数据量比较适中、计算时间充足的情况
• "arpack"：分解时就将特征数量降到n_components中输入的数值k，可以加快运算速度，适合特征矩阵很大时，但一般用于特征矩阵为稀疏矩阵（每一列都是由0和1组成的，且大部分为0）的情况，此过程包含一定的随机性
• "randomized"：随机SVD，比"full"模式下计算快很多，且还能够保证模型运行效果，适合特征矩阵巨大、计算量庞大的情况

通常选用"auto"，算不出来的话换"randomized"

参数 random_state 在参数 svd_solver 的值为 "arpack" 或 "randomized" 时生效，可以控制这两种SVD模式中的随机模式

3. 重要属性 components_

V(k,n) 是新特征空间，是要将原始数据进行映射的那些新特征向量组成的矩阵，用它来计算新特征矩阵，用属性components_调用查看

在矩阵分解时，PCA是有目标的：在原有特征的基础上，找出能够让信息尽量聚集的新特征向量。在sklearn使用的PCA和SVD联合的降维方法中，新特征向量组成的新特征空间就是 V(k,n)

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

例：人脸识别数据集中属性components_的运用

数据集 fetch_lfw_people 介绍： 

sklearn.datasets.fetch_lfw_people — scikit-learn 1.2.0 documentation

机器学习笔记（十五）：人脸识别_云布道师的博客-CSDN博客

这个数据集很大，有200+M，直接用fetch命令下载却失败了，手动下载数据

链接：https://pan.baidu.com/s/11ebeCTH7E24XAgYVL7y_-A

提取码：3gut

然后在使用 fetch_lfw_people() 时增加 data_home 参数，指定存放数据集的目录，sklearn会去指定的目录下解压 lfw-funneled 压缩包

（1）实例化数据集，探索数据 

from sklearn.datasets import fetch_lfw_people   #4个人的1000多张人脸图片组成的一组人脸数据
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
import numpy as np
 
# 实例化数据集，探索数据
faces = fetch_lfw_people(data_home = 'D:\Anaconda\Lib\site-packages\sklearn\datasets\data\lfw-funneled',min_faces_per_person=60)  #实例化   min_faces_per_person=60：每个人取出60张脸图
faces   #一个字典形式的数据

faces.data.shape  #（820,2914）   行是样本，列是样本相关的所有特征：2914 = 62 * 47
 
faces.images.shape  #（820,62,47）  820是矩阵中图像的个数，62是每个图像的特征矩阵的行，47是每个图像的特征矩阵的列
# 若是三维，先写第三维，再是行列，即(第三维,行,列)
 
# faces.data是特征矩阵，但faces.images才是用来画图的矩阵
X = faces.data

（2）将原特征矩阵进行可视化（将像素转化为图像） 

#数据本身是图像，和数据本身只是数字，使用的可视化方法不同
 
#创建画布和子图对象
'''
plt.figure()无法画多个并列的图，除非给figure的尺寸
plt.subplots()专门用来画子图和建立子图画布的
'''
fig, axes = plt.subplots(4,5  # 4行5列个图，即20张子图
                        ,figsize=(8,4)  # figsize指的是图的尺寸和比例
                        ,subplot_kw = {"xticks":[],"yticks":[]}  #不要显示坐标轴
                        )

'''
axes.shape是(4,5)，是二维结构，可以有两种循环方式：
一种是使用索引，循环一次同时生成一列上的四个图，循环五次即可；
另一种是把数据拉成一维，循环一次只生成一个图，需要循环20次
在这里，究竟使用哪一种循环方式，是要看我们要画的图的信息储存在一个怎样的结构里
我们使用 子图对象.imshow 来将图像填充到空白画布上
而imshow要求的数据格式必须是一个(m,n)格式的矩阵，即每个数据都是一张单独的图
因此我们需要遍历的是faces.images，其结构是(1277, 62, 47)
要从一个数据集中取出20个图，明显是一次性的循环切片[i,:,:]来的便利（i指第i张图，后面两个数是每张图的特征矩阵的行和列）
因此我们要把axes的结构拉成一维来循环
'''
axes[0][0].imshow(faces.images[0,:,:])   # 第一张图片
 
axes.flat  #降低一个维度，从二维到一维。是一个惰性对象
 
[*axes.flat] # 2维
# 惰性对象用[*对象]就可以看到内容
 
len([*axes.flat]) # 20

enumerate(axes.flat)  # 惰性对象
 
[*enumerate(axes.flat)]  # 给之前的20个对象分别加了索引，组成元组(索引，对象)，放在列表里

#填充图像
for i, ax in enumerate(axes.flat):  # i是索引（用来控制循环20次），ax是画图对象
    ax.imshow(faces.images[i,:,:] 
              ,cmap="gray" #选择色彩的模式
            )
 
# cmap参数取值选择各种颜色：https://matplotlib.org/tutorials/colors/colormaps.html

（3）建模降维，提取新特征空间矩阵 

#原本有 62*47=2914维，我们现在来降到150维
pca = PCA(150).fit(X)  #这里X = faces.data，不是faces.images,因为sklearn只接受2维数组降，不接受高维数组降
 
x_dr = pca.transform(X)  # 降维后的数组
x_dr.shape  #(820,150)
 
V = pca.components_  #V(k,n)用来映射的新特征向量空间，决定了新特征有什么含义、叫什么名字、是什么方向
# V * 原有特征矩阵X = 降维后的矩阵
V.shape  # V（k，n）   (150, 2914)

（4）将新特征空间矩阵可视化

V[0].shape  #(2914,)
V[0].reshape(62,47).shape  #(62, 47)
 
fig, axes = plt.subplots(4,5,figsize=(8,4),subplot_kw = {"xticks":[],"yticks":[]})
 
for i, ax in enumerate(axes.flat):
    ax.imshow(V[i,:].reshape(62,47),cmap="gray")

比起降维前的数据，新特征空间可视化后的人脸非常模糊，这是因为原始数据还没有被映射到特征空间中。但是可以看出，整体比较亮的图片获取的信息比较多，整体比较暗的图片却只能看见黑漆漆一块；在比较亮的图片中，眼睛、鼻子、嘴巴都相对清晰，脸的轮廓、头发等比较模糊。这说明，新特征空间里的特征向量大部分是 “五官” 和 “亮度” 相关的向量，所以新特征向量上的信息肯定大部分是由原数据中和 “五官”、“亮度” 相关的特征中提取出来的

通过可视化新特征空间V解释了一部分降维后的特征：画出来的图表示这些特征是和 “五官”、“亮度” 相关的，说明PCA能够将原始数据集中重要的数据进行聚集

人脸识别的核心是瞳孔/虹膜识别

4. 重要接口inverse_transform

特征工程中学习了接口 inverse_transform，可以将归一化、标准化、做过哑变量的特征矩阵都还原回原始数据中的特征矩阵

原特征矩阵X(m,n) * 新特征空间矩阵V(k,n) 的转置 = 新特征矩阵X_dr(m,k)   

是否可让新特征矩阵X_dr还原为X？

Answer：inverse_transform没有将降维逆转，只是将数据重新映射到了原数据所在的特征空间中（降维不是完全可逆的）

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

（1）导入数据并探索

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(data_home = 'D:\Anaconda\Lib\site-packages\sklearn\datasets\data\lfw-funneled', min_faces_per_person=60)
faces.images.shape   # (820, 62, 47)
faces.data.shape   # (820, 2914)
X = faces.data

（2）建模降维，获取降维后的特征矩阵X_dr

pca = PCA(150)  #实例化
X_dr = pca.fit_transform(X)  #拟合+提取结果
X_dr.shape   # (820, 150)

（3）将降维后矩阵用 inverse_tansform 返回原空间 

X_inverse = pca.inverse_transform(X_dr)
# 期待X_inverse和原数据有相同的结果，若相同则说明inverse_transform实现了降维过程的逆转    
X_inverse.shape   #(820, 2914)

（4）将特征矩阵 X 和 X_inverse 可视化

fig, ax = plt.subplots(2,10,figsize=(10,2.5)
                      ,subplot_kw={"xticks":[],"yticks":[]}
                     )
 
#需要对子图对象进行遍历的循环，来将图像填入子图中
#那在这里，我们使用怎样的循环？
#现在我们的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")   # 第二行
# .imshow()只接受二维数组，且二维必须要是图像数据上的行特征和列特征

inverse_transform并没有实现数据的完全逆转，这是因为，在降维的时候部分信息已经被舍弃了，X_dr中往往不会包含原数据100%的信息，所以在逆转时即便维度升高，原数据中已经被舍弃的信息也不可能再回来了

inverse_transform的功能，是基于X_dr中的数据进行升维，将数据重新映射到原数据所在的特征空间中，而非恢复所有原有数据。但同时，降维后的数据确实保留了原数据的大部分信息，所以图像看起来才会和原数据高度相似，只是稍稍模糊 —— 提取和储存的信息可以很少，但还是可以有很高的匹配精度（火车站等地方人脸识别机器识别速度很快的原因）