[Python] 什么是PCA降维技术以及scikit-learn中PCA类使用案例（图文教程，含详细代码）

什么是维度？

对于Numpy中数组来说，维度就是功能shape返回的结果，shape中返回了几个数字，就是几维。索引以外的数据，不分行列的叫一维（此时shape返回唯一的维度上的数据个数），有行列之分叫二维（shape返回行x列），也称为表。一张表最多二维，复数的表构成了更高的维度。当一个数组中存在2张3行4列的表时，shape返回的是(更高维，行，列)。当数组中存在2组2张3行4列的表时，数据就是4维，shape返回(2,2,3,4)。

数组中的每一张表，都可以是一个特征矩阵，这些结构永远只有一张表，所以一定有行列，其中行是样本，列是特征。针对每一张表，维度指的是样本的数量或特征的数量，一般无特别说明，指的都是特征的数量。除了索引之外，一个特征是一维，两个特征是二维，n个特征是n维。

对图像来说，维度就是图像中特征向量的数量。特征向量可以理解为是坐标轴，一个特征向量定义一条直线，是一维，两个相互垂直的特征向量定义一个平面，即一个直角坐标系，就是二维，三个相互垂直的特征向量定义一个空间，即一个立体直角坐标系，就是三维。三个以上的特征向量相互垂直，定义人眼无法看见，也无法想象的高维空间。

什么是PCA降维技术？

PCA（principal components analysis）即主成分分析技术，又称主分量分析，旨在利用降维的思想，把多特征指标转化为少数几个综合指标。

在统计学中，主成分分析PCA是一种简化数据集的技术。它是一个线性变换。这个变换把数据变换到一个新的坐标系统中，使得任何数据投影的第一大方差在第一个坐标(称为第一主成分)上，第二大方差在第二个坐标(第二主成分)上，依次类推。主成分分析经常用于减少数据集的维数，同时保持数据集的对方差贡献最大的特征。

降维算法中的”降维“，指的是降低特征矩阵中特征的数量，降维的目的是为了让算法运算更快，效果更好，还有另一种需求，就是数据可视化。从上面的图我们其实可以看得出，图像和特征矩阵的维度是可以相互对应的，即一个特征对应一个特征向量，对应一条坐标轴。所以，三维及以下的特征矩阵，是可以被可视化的，这可以帮助我们很快地理解数据的分布，而三维以上特征矩阵的则不能被可视化，数据的性质也就比较难理解。

什么是SVD降维技术？

SVD降维技术是一种基于奇异值分解（Singular Value Decomposition，SVD）的数据降维方法。SVD降维技术是一种矩阵分解技术，用于将一个矩阵分解为三个矩阵的乘积。在降维过程中，SVD将原始数据矩阵分解为三个矩阵，分别为左奇异向量矩阵、奇异值矩阵和右奇异向量矩阵。

SVD降维技术的关键思想是将原始数据矩阵投影到奇异值矩阵上，通过保留最重要的奇异值和对应的奇异向量，将数据的维度减少到较低的维度。这样可以去除数据中的噪声和冗余信息，保留了数据的主要特征。

SVD降维技术在特征提取和数据压缩等领域广泛应用。通过降低维度，可以减少计算复杂度和存储空间，提高算法的效率和性能。同时，降维还可以帮助可视化高维数据和发现数据的内在结构。

PCA降维究竟是怎样实现？

我们来看一组简单的二维数据的降维：

我们现在有一组简单的数据，有特征x1和x2，三个样本数据的坐标点分别为(1,1)，(2,2)，(3,3)。我们可以让x1和 x2分别作为两个特征向量，很轻松地用一个二维平面来描述这组数据。这组数据现在每个特征的均值都为2，方差则等于：

每个特征的数据一模一样，因此方差也都为1，数据的方差总和是2。现在我们的目标是：只用一个特征向量来描述这组数据，即将二维数据降为一维数据，并且尽可能地保留信息量，即让数据的总方差尽量靠近2。于是，我们将原本的直角坐标系逆时针旋转45°，形成了新的特征向量x1*和x2*组成的新平面，在这个新平面中，三个样本数据的坐标点可以表示为 $(\sqrt{2},0)$ ， $(2\sqrt{2}, 0)$ ， $(3\sqrt{2},0)$ 。可以注意到，x2*上的数值此时都变成了0，因此x2*明显不带有任何有效信息了（此时x2*的方差也为0了）。此时， x1*特征上的数据均值是 $2\sqrt{2}$ ，而方差则可表示成:

此时，我们根据信息含量的排序，取信息含量最大的一个特征，因为我们想要的是一维数据。所以我们可以将x2* 删除，同时也删除图中的x2*特征向量，剩下的x1*就代表了曾经需要两个特征来代表的三个样本点。通过旋转原有特征向量组成的坐标轴来找到新特征向量和新坐标平面，我们将三个样本点的信息压缩到了一条直线上，实现了二维变一维，并且尽量保留原始数据的信息。一个成功的降维，就实现了。

在步骤3当中，我们用来找出n个新特征向量，让数据能够被压缩到少数特征上并且总信息量不损失太多的技术就是矩阵分解。PCA和SVD是两种不同的降维算法，但他们都遵从上面的过程来实现降维，只是两种算法中矩阵分解的方法不同，信息量的衡量指标不同罢了。PCA使用方差作为信息量的衡量指标，并且特征值分解来找出空间V。降维完成之后，PCA找到的每个新特征向量就叫做“主成分”，而被丢弃的特征向量被认为信息量很少，这些信息很可能就是噪音。

降维和特征选择都是特征工程技术，它们有什么不同？

特征工程中有三种方式：特征提取，特征创造和特征选择。

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

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

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

scikit-learn中PCA类介绍

sklearn中降维算法都被包括在模块decomposition中，

API Reference — scikit-learn 1.4.0 documentation

SVD（singular value decomposition）和主成分分析PCA（principal components analysis）都属于矩阵分解算法中的入门算法，都是通过分解特征矩阵来进行降维，它可以将高维数据映射到低维空间中，同时保留尽可能多的数据信息。

PCA类原型

PCA类原型如下：

Parameters
----------
n_components : int, float or 'mle', default=None
    Number of components to keep.
 
    if n_components is not set all components are kept::

        n_components == min(n_samples, n_features)

    If ``n_components == 'mle'`` and ``svd_solver == 'full'``, Minka's
    MLE is used to guess the dimension. Use of ``n_components == 'mle'``
    will interpret ``svd_solver == 'auto'`` as ``svd_solver == 'full'``.

    If ``0 < n_components < 1`` and ``svd_solver == 'full'``, select the
    number of components such that the amount of variance that needs to be
    explained is greater than the percentage specified by n_components.

    If ``svd_solver == 'arpack'``, the number of components must be
    strictly less than the minimum of n_features and n_samples.

    Hence, the None case results in::

        n_components == min(n_samples, n_features) - 1

copy : bool, default=True
    If False, data passed to fit are overwritten and running
    fit(X).transform(X) will not yield the expected results,
    use fit_transform(X) instead.

whiten : bool, default=False
    When True (False by default) the `components_` vectors are multiplied
    by the square root of n_samples and then divided by the singular values
    to ensure uncorrelated outputs with unit component-wise variances.

    Whitening will remove some information from the transformed signal
    (the relative variance scales of the components) but can sometime
    improve the predictive accuracy of the downstream estimators by
    making their data respect some hard-wired assumptions.

svd_solver : {'auto', 'full', 'arpack', 'randomized'}, default='auto'
    If auto :
        The solver is selected by a default policy based on `X.shape` and
        `n_components`: if the input data is larger than 500x500 and the
        number of components to extract is lower than 80% of the smallest
        dimension of the data, then the more efficient 'randomized'
        method is enabled. Otherwise the exact full SVD is computed and
        optionally truncated afterwards.
    If full :
        run exact full SVD calling the standard LAPACK solver via
        `scipy.linalg.svd` and select the components by postprocessing
    If arpack :
        run SVD truncated to n_components calling ARPACK solver via
        `scipy.sparse.linalg.svds`. It requires strictly
        0 < n_components < min(X.shape)
    If randomized :
        run randomized SVD by the method of Halko et al.

    .. versionadded:: 0.18.0

tol : float, default=0.0
    Tolerance for singular values computed by svd_solver == 'arpack'.
    Must be of range [0.0, infinity).

    .. versionadded:: 0.18.0

iterated_power : int or 'auto', default='auto'
    Number of iterations for the power method computed by
    svd_solver == 'randomized'.
    Must be of range [0, infinity).

    .. versionadded:: 0.18.0

n_oversamples : int, default=10
    This parameter is only relevant when `svd_solver="randomized"`.
    It corresponds to the additional number of random vectors to sample the
    range of `X` so as to ensure proper conditioning. See
    :func:`~sklearn.utils.extmath.randomized_svd` for more details.

    .. versionadded:: 1.1

power_iteration_normalizer : {'auto', 'QR', 'LU', 'none'}, default='auto'
    Power iteration normalizer for randomized SVD solver.
    Not used by ARPACK. See :func:`~sklearn.utils.extmath.randomized_svd`
    for more details.

    .. versionadded:: 1.1

random_state : int, RandomState instance or None, default=None
    Used when the 'arpack' or 'randomized' solvers are used. Pass an int
    for reproducible results across multiple function calls.
    See :term:`Glossary <random_state>`.

    .. versionadded:: 0.18.0

类参数说明

n_components：指定降维后的维度数。如果未指定，则保持所有特征。
copy：默认为True。如果为True，则将原始数据的副本进行降维，否则直接在原始数据上进行计算。
默认为False。如果为True，则对数据进行白化处理，保留所有特征的单位方差。
默认为'auto'。指定SVD求解器的类型。可选值为'auto'、'full'、'arpack'和'randomized'。
默认为0.0。指定迭代过程的收敛阈值。
默认为'auto'。指定幂迭代算法的迭代次数。可选值为'auto'和'exact'。
默认为None。指定随机数生成器的种子。

类的主要方法

除了初始化参数外，PCA类中还提供了以下方法：

fit(X[, y])：计算PCA模型。
fit_transform(X[, y])：计算PCA模型并进行降维。
transform(X)：将数据进行降维。
inverse_transform(X)：将降维后的数据恢复到原始维度。
get_covariance()：返回数据的协方差矩阵。
get_params([deep])：获取模型的参数。
get_precision()：返回数据的精度矩阵。
score(X[, y])：计算样本的似然值或平均对数似然值。
score_samples(X)：计算样本的似然值或平均对数似然值。
set_params(**params)：设置模型的参数。

重要参数 - n_components

n_components是我们降维后需要的维度，即降维后需要保留的特征数量，降维流程中第二步里需要确认的k值，一般输入[0, min(X.shape)]范围中的整数。一说到K，大家可能都会想到，类似于KNN中的K和随机森林中的 n_estimators，这是一个需要我们人为去确认的超参数，并且我们设定的数字会影响到模型的表现。如果留下的特征太多，就达不到降维的效果，如果留下的特征太少，那新特征向量可能无法容纳原始数据集中的大部分信息，因此，n_components既不能太大也不能太小。那怎么办呢？可以先从我们的降维目标说起：如果我们希望可视化一组数据来观察数据分布，我们往往将数据降到三维以下，很多时候是二维，即n_components的取值为2。

参考使用案例中的 “高维数据的可视化”、“累积可解释方差贡献率曲线”、“最大似然估计自选超参数”和“按信息量占比选超参数”。

重要参数 - svd_solver

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

scikit-learn中PCA函数使用案例

高维数据的可视化

# 导入包模块
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
import pandas as pd

# 导入鸢尾花数据集
iris_ds = load_iris()
y = iris_ds.target
X = iris_ds.data
#作为数据表或特征矩阵，X是几维？
print(X.shape)
pd.DataFrame(X).head()

#调用PCA
pca = PCA(n_components=2) #实例化, 将4维的列特征降维2维
pca = pca.fit(X) #拟合模型
X_dr = pca.transform(X) #获取新矩阵
#也可以fit_transform一步到位
#X_dr = PCA(2).fit_transform(X)
print(X_dr.shape)
pd.DataFrame(X_dr).head()

# 特征可视化 
# plt.figure()
# plt.scatter(X_dr[y==0, 0], X_dr[y==0, 1], c="red", label=iris_ds.target_names[0])
# plt.scatter(X_dr[y==1, 0], X_dr[y==1, 1], c="black", label=iris_ds.target_names[1])
# plt.scatter(X_dr[y==2, 0], X_dr[y==2, 1], c="orange", label=iris_ds.target_names[2])
# plt.legend()
# plt.title('PCA of IRIS dataset')
# plt.show()

# 或者
colors = ['red', 'black', 'orange']
iris_ds.target_names
plt.figure()
for i in [0, 1, 2]:  # 0, 1, 2为鸢尾花的分类id
    plt.scatter(X_dr[y == i, 0]  # x轴为降维后的第一列的特征
               ,X_dr[y == i, 1]  # y轴为降维后的第二列的特征
               ,alpha=.7
               ,c=colors[i]
               ,label=iris_ds.target_names[i]
               )
plt.legend()
plt.title('PCA of IRIS dataset')
plt.show()

鸢尾花的分布被展现在我们眼前了，从上图可以看出明显是一个分簇的分布，并且每个簇之间的分布相对比较明显，也许 versicolor和virginia这两种花之间会有一些分类错误，但setosa肯定不会被分错。这样的数据很容易分类，可以预见，KNN，随机森林，神经网络，朴素贝叶斯，Adaboost这些分类器在鸢尾花数据集上，未调整的时候都可以有 95%上下的准确率。

# 探索降维后的数据

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

累积可解释方差贡献率曲线

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

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

iris_ds = load_iris()
y = iris_ds.target
X = iris_ds.data

pca_line = PCA().fit(X)
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")
plt.show()

最大似然估计自选超参数

除了输入整数，n_components还有哪些选择呢？矩阵分解的理论发展在业界独树一帜，数学大神Minka, T.P.在麻省理工学院媒体实验室做研究时找出了让PCA用最大似然估计(maximum likelihood estimation)自选超参数的方法，输入“mle”作为n_components的参数输入，就可以调用这种方法。

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

iris_ds = load_iris()
y = iris_ds.target
X = iris_ds.data
pca_mle = PCA(n_components="mle")
pca_mle = pca_mle.fit(X)
X_mle = pca_mle.transform(X)
pd.DataFrame(X_mle).head()
#可以发现，mle为我们自动选择了3个特征

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

按信息量占比选超参数

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

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

iris_ds = load_iris()
y = iris_ds.target
X = iris_ds.data
pca_f = PCA(n_components=0.97,svd_solver="full")
pca_f = pca_f.fit(X)
X_f = pca_f.transform(X)
pd.DataFrame(X_f).head()