降维 – PCA 解释

Question

我不认为我对 PCA 有很好的了解，有人可以帮我解决下面的困惑吗？

以鸢尾花数据集为例，我有4个协变量，x1:sepal length； x2：萼片宽度； x3：花瓣长度； x4：花瓣宽度。公式如下所示，a1,a2,a3,a4 是协变量的权重。并且 PCA 将尝试使用不同的线性变换来最大化方差。同时也遵循 a1^2 + a2^2 + a3^2 + a4^2=1 的规则。我有兴趣知道 a1,a2,a3,a4 的值。

a1*x1 + a2*x2 + a3*x3 + a4*x4

我在 python 上有以下我认为正确的代码？

# load libraries
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
import seaborn as sns
import pandas as pd
import numpy as np

iris = load_iris()
X = iris.data
df = pd.DataFrame(X,columns=iris.feature_names)

pca = decomposition.PCA(n_components = 4)
digits_pca_4 = pca.fit(X)
digits_pca_4.explained_variance_ratio_

结果是

array([0.92461872, 0.05306648, 0.01710261, 0.00521218])

我的问题是：

我是否正确假设 a1=sqrt(0.92), a2=sqrt(0.05), a3=sqrt(0.02), a4=sqrt(0.005)？

第二个问题：

如果我选择 a1=a2=a3=a4=0.5 的线性组合，与 PCA 的方差相比，它的方差是多少（我假设它小于 PCA 结果，因为 PCA 最大化了方差？）？如何在 python 中获得 a1=a2=a3=a4=0.5 时的方差？下面的代码是来自 PCA 的差异吗？

pca.explained_variance_.sum()

非常感谢！

Answer 1

直接回答您的问题：不，您最初的解释不正确。

说明

PCA 完成的实际投影是矩阵乘法Y = (X - u) W，其中u 是X 的平均值（u = X.mean(axis=0)），W 是PCA 找到的投影矩阵：n x p 正交矩阵，其中n 是原始数据维度，p 是所需的输出维度。您给出的表达式 (a1*x1 + a2*x2 + a3*x3 + a4*x4) 并不意味着所有值都是标量。充其量，它可能意味着计算单个组件，使用下面W 的一列j 作为a_k：Y[i, j] == sum(W[k, j] * (X[i, k] - u[k]) for k in range(n))。

在任何情况下，您都可以使用vars(pca) 检查pca = PCA.fit(...) 结果的所有变量。特别是，上述投影矩阵可以找到为W = pca.components_.T。可以验证以下陈述：

# projection
>>> u = pca.mean_
... W = pca.components_.T
... Y = (X - u).dot(W)
... np.allclose(Y, pca.transform(X))
True

>>> np.allclose(X.mean(axis=0), u)
True

# orthonormality
>>> np.allclose(W.T.dot(W), np.eye(W.shape[1]))
True

# explained variance is the sample variation (not population variance)
# of the projection (i.e. the variance along the proj axes)
>>> np.allclose(Y.var(axis=0, ddof=1), pca. explained_variance_)
True

图形演示

理解 PCA 的最简单方法是它纯粹是 n-D 中的 旋转（均值去除后），同时仅保留第一个 p 维。旋转使得您的数据的最大方差方向与投影中的自然轴对齐。

这是一个小演示代码，可帮助您直观地了解正在发生的事情。另请阅读Wikipedia page on PCA。

def pca_plot(V, W, idx, ax):
    # plot only first 2 dimensions of W along with axes W
    colors = ['k', 'r', 'b', 'g', 'c', 'm', 'y']
    u = V.mean(axis=0)  # n
    axes_lengths = 1.5*(V - u).dot(W).std(axis=0)
    axes = W * axes_lengths  # n x p
    axes = axes[:2].T  # p x 2
    ax.set_aspect('equal')
    ax.scatter(V[:, 0], V[:, 1], alpha=.2)
    ax.scatter(V[idx, 0], V[idx, 1], color='r')
    hlen = np.max(np.linalg.norm((V - u)[:, :2], axis=1)) / 25
    for k in range(axes.shape[0]):
        ax.arrow(*u[:2], *axes[k], head_width=hlen/2, head_length=hlen, fc=colors[k], ec=colors[k])

def pca_demo(X, p):
    n = X.shape[1]  # input dimension
    pca = PCA(n_components=p).fit(X)
    u = pca.mean_
    v = pca.explained_variance_
    W = pca.components_.T
    Y = pca.transform(X)
    assert np.allclose((X - u).dot(W), Y)
    
    # plot first 2D of both input space and output space
    # for visual identification: select a point that's as far as possible
    # in the direction of the diagonal of the axes cube, after normalization
    # Z: variance-1 projection
    Z = (X - u).dot(W/np.sqrt(v))
    idx = np.argmax(Z.sum(axis=1) / np.sqrt(np.linalg.norm(Z, axis=1)))

    fig, ax = plt.subplots(ncols=2, figsize=(12, 6))

    # input space
    pca_plot(X, W, idx, ax[0])
    ax[0].set_title('input data (first 2D)')

    # output space
    pca_plot(Y, np.eye(p), idx, ax[1])
    ax[1].set_title('projection (first 2D)')
    
    return Y, W, u, pca

示例

虹膜数据

# to better understand the shape of W, we project onto
# a space of dimension p=3
X = load_iris().data
Y, W, u, pca = pca_demo(X, 3)

请注意，投影实际上只是(X - u) W：

>>> np.allclose((X - u).dot(W), Y)
True

合成椭球数据

A = np.array([
    [20, 10, 7],
    [-1, 3, 7],
    [5, 1, 2],
])
X = np.random.normal(size=(1000, A.shape[0])).dot(A)
Y, W, u, pca = pca_demo(X, 3)