主成分分析：用于降维的主成分分析（PCA）

主成分分析（PCA）降维的实现

> from mlxtend.feature_extraction import PrincipalComponentAnalysis

# 主成分分析（PCA）
主成分分析（PCA）是一种统计技术，用于将数据集降维，同时保留数据中最重要的特征。

## PCA的步骤
1. 标准化数据。
2. 计算协方差矩阵。
3. 计算特征值和特征向量。
4. 选择主成分。
5. 转换数据集。

通过以上步骤，我们可以将高维数据投影到低维空间中，同时尽可能保留数据的变化信息。

## 示例
下面的代码示例演示了如何使用`mlxtend`库进行主成分分析。

```python
from mlxtend.feature_extraction import PrincipalComponentAnalysis

# 假设我们有一个数据集 X
pca = PrincipalComponentAnalysis(n_components=2)
X_pca = pca.fit_transform(X)

# 现在 X_pca 是压缩后的数据集



## 概述


现代时代数据的庞大规模不仅是计算机硬件的挑战，也是许多机器学习算法性能的主要瓶颈。PCA分析的主要目标是识别数据中的模式；PCA旨在检测变量之间的相关性。如果变量之间存在强相关性，那么尝试降低维度是有意义的。简而言之，这就是PCA的核心内容：在高维数据中寻找最大方差的方向，并将其投影到较小维度的子空间，同时尽可能保留大部分信息。


### PCA和降维

通常，期望的目标是通过将$d$维数据集投影到一个$(k)$维子空间（其中$k\;<\;d$）来减少维度，以提高计算效率，同时保留大部分信息。一个重要的问题是“什么是能够‘很好’表示数据的$k$的大小？”

稍后，我们将计算数据集的特征向量（主成分）并将它们收集在一个投影矩阵中。每个特征向量都与一个特征值相关联，该特征值可以被解释为对应特征向量的“长度”或“大小”。如果一些特征值的大小显著大于其他特征值，那么通过丢弃“信息较少”的特征对将数据集通过PCA降到较小维数的子空间就是合理的。


### PCA 方法总结

- 标准化数据。
- 从协方差矩阵或相关矩阵中获取特征向量和特征值，或者执行奇异值分解。
- 将特征值按降序排列，并选择与 $k$ 个最大特征值对应的 $k$ 个特征向量，其中 $k$ 是新特征子空间的维度数（$k \le d$）。
- 从选定的 $k$ 个特征向量构造投影矩阵 $\mathbf{W}$。
- 通过 $\mathbf{W}$ 转换原始数据集 $\mathbf{X}$，以获得 $k$ 维特征子空间 $\mathbf{Y}$。


### 参考文献

- Pearson, Karl. "LIII. [关于与空间中点系统的最接近拟合线和平面。](https://www.tandfonline.com/doi/abs/10.1080/14786440109462720?journalCode=tphm17)" 伦敦、爱丁堡和都柏林哲学杂志与科学期刊 2.11 (1901): 559-572.


## 示例 1 - 在鸢尾花数据集上进行PCA分析



```python
from mlxtend.data import iris_data
from mlxtend.preprocessing import standardize
from mlxtend.feature_extraction import PrincipalComponentAnalysis

X, y = iris_data()
X = standardize(X)

pca = PrincipalComponentAnalysis(n_components=2)
pca.fit(X)
X_pca = pca.transform(X)

import matplotlib.pyplot as plt

with plt.style.context('seaborn-whitegrid'):
    plt.figure(figsize=(6, 4))
    for lab, col in zip((0, 1, 2),
                        ('blue', 'red', 'green')):
        plt.scatter(X_pca[y==lab, 0],
                    X_pca[y==lab, 1],
                    label=lab,
                    c=col)
    plt.xlabel('Principal Component 1')
    plt.ylabel('Principal Component 2')
    plt.legend(loc='lower center')
    plt.tight_layout()
    plt.show()

png

示例 2 - 绘制解释方差比率

from mlxtend.data import iris_data
from mlxtend.preprocessing import standardize

X, y = iris_data()
X = standardize(X)

pca = PrincipalComponentAnalysis(n_components=None)
pca.fit(X)
X_pca = pca.transform(X)

pca.e_vals_

array([2.91081808, 0.92122093, 0.14735328, 0.02060771])

pca.e_vals_normalized_

array([0.72770452, 0.23030523, 0.03683832, 0.00515193])

import numpy as np

tot = sum(pca.e_vals_)
var_exp = [(i / tot)*100 for i in sorted(pca.e_vals_, reverse=True)]
cum_var_exp = np.cumsum(pca.e_vals_normalized_*100)

with plt.style.context('seaborn-whitegrid'):
    fig, ax = plt.subplots(figsize=(6, 4))
    plt.bar(range(4), var_exp, alpha=0.5, align='center',
            label='individual explained variance')
    plt.step(range(4), cum_var_exp, where='mid',
             label='cumulative explained variance')
    plt.ylabel('Explained variance ratio')
    plt.xlabel('Principal components')
    plt.xticks(range(4))
    ax.set_xticklabels(np.arange(1, X.shape[1] + 1))
    plt.legend(loc='best')
    plt.tight_layout()

png

示例 3 - 通过奇异值分解进行主成分分析 (PCA)

虽然协方差或相关矩阵的特征分解可能更直观，但大多数PCA实现采用奇异值分解（SVD）以提高计算效率。使用SVD的另一个优点是结果往往具有更好的数值稳定性，因为我们可以直接分解输入矩阵，而不需要额外的协方差矩阵步骤。

from mlxtend.data import iris_data
from mlxtend.preprocessing import standardize
from mlxtend.feature_extraction import PrincipalComponentAnalysis

X, y = iris_data()
X = standardize(X)

pca = PrincipalComponentAnalysis(n_components=2,
                                 solver='svd')
pca.fit(X)
X_pca = pca.transform(X)

import matplotlib.pyplot as plt

with plt.style.context('seaborn-whitegrid'):
    plt.figure(figsize=(6, 4))
    for lab, col in zip((0, 1, 2),
                        ('blue', 'red', 'green')):
        plt.scatter(X_pca[y==lab, 0],
                    X_pca[y==lab, 1],
                    label=lab,
                    c=col)
    plt.xlabel('Principal Component 1')
    plt.ylabel('Principal Component 2')
    plt.legend(loc='lower center')
    plt.tight_layout()
    plt.show()

png

如果我们将这个PCA投影与示例1中的先前图进行比较，我们会注意到它们是彼此的镜像。这并不是因为这两个实现中的任何一个出现错误，而是由于这个差异的原因是，取决于特征求解器，特征向量可以具有正或负符号。

例如，如果$v$是矩阵$\Sigma$的一个特征向量，我们有

$$\Sigma v = \lambda v,$$

其中$\lambda$是我们的特征值。

那么$-v$也是一个具有相同特征值的特征向量，因为

$$\Sigma(-v) = -\Sigma v = -\lambda v = \lambda(-v)。$$

示例 4 - 因子负载量

在调用 fit 方法后，因子载荷可以通过 loadings_ 属性获得。简单来说，载荷是特征向量的未标准化值。换句话说，我们可以将载荷解释为输入特征与主成分（或特征向量）之间的协方差（如果我们标准化了输入特征，则为相关性），这些主成分已经缩放到单位长度。

通过缩放载荷，它们在幅度上变得可比较，我们可以评估在一个成分中有多少方差归因于输入特征（因为成分仅仅是输入特征的加权线性组合）。

from mlxtend.data import iris_data
from mlxtend.preprocessing import standardize
from mlxtend.feature_extraction import PrincipalComponentAnalysis
import matplotlib.pyplot as plt

X, y = iris_data()
X = standardize(X)

pca = PrincipalComponentAnalysis(n_components=2,
                                 solver='eigen')
pca.fit(X);

xlabels = ['sepal length', 'sepal width', 'petal length', 'petal width']

fig, ax = plt.subplots(1, 2, figsize=(8, 3))

ax[0].bar(range(4), pca.loadings_[:, 0], align='center')
ax[1].bar(range(4), pca.loadings_[:, 1], align='center')

ax[0].set_ylabel('Factor loading onto PC1')
ax[1].set_ylabel('Factor loading onto PC2')

ax[0].set_xticks(range(4))
ax[1].set_xticks(range(4))
ax[0].set_xticklabels(xlabels, rotation=45)
ax[1].set_xticklabels(xlabels, rotation=45)
plt.ylim([-1, 1])
plt.tight_layout()

png

例如，我们可以说，第一个主成分中的大部分方差归因于花瓣特征（尽管萼片长度在主成分1上的载荷也不逊色）。相比之下，主成分2捕捉到的剩余方差主要是由于萼片宽度。请注意，从示例2中我们知道主成分1解释了大部分方差，并且基于载荷图中的信息，我们可以说花瓣特征和萼片长度的组合可能解释了数据中的大部分分布。

示例 5 - 特征提取管道

from sklearn.pipeline import make_pipeline
from sklearn.neighbors import KNeighborsClassifier
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from mlxtend.data import wine_data

X, y = wine_data()
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=123, test_size=0.3, stratify=y)

pipe_pca = make_pipeline(StandardScaler(),
                         PrincipalComponentAnalysis(n_components=3),
                         KNeighborsClassifier(n_neighbors=5))

pipe_pca.fit(X_train, y_train)


print('Transf. training accyracy: %.2f%%' % (pipe_pca.score(X_train, y_train)*100))
print('Transf. test accyracy: %.2f%%' % (pipe_pca.score(X_test, y_test)*100))

Transf. training accyracy: 96.77%
Transf. test accyracy: 96.30%

示例 6 - 白化处理

某些算法要求对数据进行白化。这意味着特征具有单位方差，且非对角线的值均为零（即，特征之间不相关）。主成分分析（PCA）已经确保特征是不相关的，因此，我们只需要对变换后的数据进行简单的缩放以实现白化。

例如，对于给定的变换特征 $X'_i$，我们将其除以对应特征值的平方根 $\lambda_i$：

$$X'_{\text{whitened}} = \frac{X'_i}{\sqrt{\lambda_i}}.$$

通过在初始化时设置 whitening=True，可以实现 PrincipalComponentAnalysis 的白化。让我们通过一个示例来演示这一点。

from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from mlxtend.data import wine_data

X, y = wine_data()
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=123, test_size=0.3, stratify=y)

常规主成分分析（PCA）

sc = StandardScaler()

pca1 = PrincipalComponentAnalysis(n_components=2)

X_train_scaled = sc.fit_transform(X_train)
X_train_transf = pca1.fit(X_train_scaled).transform(X_train_scaled)


with plt.style.context('seaborn-whitegrid'):
    plt.figure(figsize=(6, 4))
    for lab, col in zip((0, 1, 2),
                        ('blue', 'red', 'green')):
        plt.scatter(X_train_transf[y_train==lab, 0],
                    X_train_transf[y_train==lab, 1],
                    label=lab,
                    c=col)
    plt.xlabel('Principal Component 1')
    plt.ylabel('Principal Component 2')
    plt.legend(loc='lower center')
    plt.tight_layout()
    plt.show()

png

np.set_printoptions(precision=1, suppress=True)

print('Covariance matrix:\n')
np.cov(X_train_transf.T)

Covariance matrix:






array([[4.9, 0. ],
       [0. , 2.5]])

如我们所见，经过转换后特征之间没有相关性，但方差不是单位方差。

PCA与白化处理

sc = StandardScaler()

pca1 = PrincipalComponentAnalysis(n_components=2, whitening=True)

X_train_scaled = sc.fit_transform(X_train)
X_train_transf = pca1.fit(X_train_scaled).transform(X_train_scaled)


with plt.style.context('seaborn-whitegrid'):
    plt.figure(figsize=(6, 4))
    for lab, col in zip((0, 1, 2),
                        ('blue', 'red', 'green')):
        plt.scatter(X_train_transf[y_train==lab, 0],
                    X_train_transf[y_train==lab, 1],
                    label=lab,
                    c=col)
    plt.xlabel('Principal Component 1')
    plt.ylabel('Principal Component 2')
    plt.legend(loc='lower center')
    plt.tight_layout()
    plt.show()

png

np.set_printoptions(precision=1, suppress=True)

print('Covariance matrix:\n')
np.cov(X_train_transf.T)

Covariance matrix:






array([[1., 0.],
       [0., 1.]])

正如我们上面所看到的，白化实现了所有特征现在具有单位方差。即，变换特征的协方差矩阵变为单位矩阵。

API

PrincipalComponentAnalysis(n_components=None, solver='svd', whitening=False)

Principal Component Analysis Class

Parameters

n_components : int (default: None)

The number of principal components for transformation. Keeps the original dimensions of the dataset if None.
solver : str (default: 'svd')

Method for performing the matrix decomposition. {'eigen', 'svd'}
whitening : bool (default: False)

Performs whitening such that the covariance matrix of the transformed data will be the identity matrix.

Attributes

w_ : array-like, shape=[n_features, n_components]

Projection matrix
e_vals_ : array-like, shape=[n_features]

Eigenvalues in sorted order.
e_vecs_ : array-like, shape=[n_features]

Eigenvectors in sorted order.
e_vals_normalized_ : array-like, shape=[n_features]

Normalized eigen values such that they sum up to 1. This is equal to what's often referred to as "explained variance ratios."
loadings_ : array_like, shape=[n_features, n_features]

The factor loadings of the original variables onto the principal components. The columns are the principal components, and the rows are the features loadings. For instance, the first column contains the loadings onto the first principal component. Note that the signs may be flipped depending on whether you use the 'eigen' or 'svd' solver; this does not affect the interpretation of the loadings though.

Examples

For usage examples, please see https://rasbt.github.io/mlxtend/user_guide/feature_extraction/PrincipalComponentAnalysis/

Methods

fit(X, y=None)

Learn model from training data.

Parameters

X : {array-like, sparse matrix}, shape = [n_samples, n_features]

Training vectors, where n_samples is the number of samples and n_features is the number of features.

Returns

self : object

get_params(deep=True)

Get parameters for this estimator.

Parameters

deep : boolean, optional

If True, will return the parameters for this estimator and contained subobjects that are estimators.

Returns

params : mapping of string to any

Parameter names mapped to their values.'

adapted from https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/base.py Author: Gael Varoquaux gael.varoquaux@normalesup.org License: BSD 3 clause

set_params(params)

Set the parameters of this estimator. The method works on simple estimators as well as on nested objects (such as pipelines). The latter have parameters of the form <component>__<parameter> so that it's possible to update each component of a nested object.

Returns

self

adapted from https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/base.py Author: Gael Varoquaux gael.varoquaux@normalesup.org License: BSD 3 clause

transform(X)

Apply the linear transformation on X.

Parameters

X : {array-like, sparse matrix}, shape = [n_samples, n_features]

Training vectors, where n_samples is the number of samples and n_features is the number of features.

Returns

X_projected : np.ndarray, shape = [n_samples, n_components]

Projected training vectors.