LionKing数据科学专栏

购买普通会员高级会员可以解锁网站精华内容且享受VIP服务的优惠

想要查看更多数据科学相关的内容请关注我们的微信公众号知乎专栏

主成分分析(principal component analysis; PCA)

主成分分析的原理

主成分分析(principal component analysis),简称PCA,是对高维数据进行降维的一种方法。

主成分分析寻找$K$个两两正交的向量$v_1, \ldots, v_K$,使得数据投影到这$K$个向量组成的子空间后,投影的距离平方和最小。

下图为二维数据第一个主成分的例子:

对于数据$x$,我们将$x$替换为其在$v_1, \ldots, v_K$张成的子空间的投影$$\widetilde{x} = \sum_{i=1}^{K}\langle x, v_i\rangle v_i$$其中$$\langle u, v \rangle = u^Tv$$为内积(inner product)。

主成分分析只需要数据$X$本身,不需要目标变量$y$,因此是非监督学习(unsupervised learning)的一种。

主成分分析的应用

主成分分析有以下应用场景:

算法

假设数据为$X = (x^{(1)}, \ldots, x^{(n)})$

  1. 对样本进行中心化$x^{(i)} = x^{(i)} - \frac{1}{n}\sum_{j=1}^{n}x^{(j)}$
  2. 计算协方差矩阵$X^TX \in \mathbb{R}^{p \times p}$的特征值分解
  3. 取出最大的$K$个特征值对应的特征向量$v_1, \ldots, v_K$,将每个特征向量标准化后得到特征向量矩阵$W \in \mathbb{R}^{p \times K}$
  4. 将数据转化为$X' = XW \in \mathbb{R}^{n \times K}$

有时还需要对数据进行归一化(rescaling)再进行主成分分析,是否使用归一化取决于应用场景:如果需要保留数据的量级则不进行归一化,例如对基因数据;如果需要尽量让每个特征量级接近则需要进行归一化,例如特征中有年龄这样取值范围较大的,也有是否工作的二进制变量这样取值范围较小的。

数值上,现有的奇异值分解(singular value decomposition, SVD)算法效率很高,因此大多数主成分分析的实现的第一步是计算$X$的奇异值分解$X = UDV^T$,而无需计算出$XX^T$和其特征值分解。

主成分分析和奇异值分解的联系

主成分分析和奇异值分解联系紧密。事实上,奇异值分解是主成分分析数值计算中的核心步骤。

考虑奇异值分解$$X = UDV^T$$其中$U \in \mathbb{R}^{n \times n}, V \in \mathbb{R}^{p \times p}$正交(orthogonal),$D \in \mathbb{R}^{n \times p}$为对角阵。

则$$X^TX = VD^TU^TUDV^T = V(D^TD)V^T$$

注意$D^TD \in \mathbb{R}^{p \times p}$为对角阵,其元素为$D$的对角线元素的平方。

且$V(D^TD)V^T$为$X^TX$的特征值分解,$V$的列为$X^TX$的特征向量,$D^TD$的元素为特征值。

我们考虑$V$的前$K$列$V_{1, \ldots, K}$,则$\widetilde{X} = XV_{1, \ldots, K} \in \mathbb{R}^{n \times K}$即数据的主成分。

$K$的选择

选择主成分个数$K$没有统一的方法,需要结合实际业务需求和数据的表现。

对于大多数问题,选择$K$往往是满足业务需求和保留原始信息的一个权衡(tradeoff)。

特点 降维 压缩 降噪 可视化
$K$小 降维后的数据维数小,模型容易训练 压缩率高 噪声小 容易解释和观察
$K$大 保留尽可能多的信息

常用的选择$K$的方式是在了解业务需求(例如,需要把数据占用空间压缩到多少MB以内、维数需要到多少以内、希望观察多少个分量)的基础上,绘制数据的方差被前$K$个主成分解释的比例关于$K$的一个折线图。

对于$K$个分量,我们将数据从$X \in \mathbb{R}^{n \times p}$变换为了$\widetilde{X} \in \mathbb{R}^{n \times K}$。

考虑协方差矩阵(covariance matrix)$\Sigma = X^TX \in \mathbb{R}^{p \times p}$的特征值(eigenvalue),从大到小为$\lambda_1, \ldots, \lambda_p$。则$\widetilde{X}$相当于只使用了其中的前$K$个特征值$\lambda_1, \ldots, \lambda_K$。我们认为数据的总方差为$tr(\Sigma) = \lambda_1 + \ldots + \lambda_p$,而$K$个分量对应的方差为$tr(\widetilde{X}^T\widetilde{X}) = \lambda_1 + \ldots + \lambda_K$,因此解释方差比例为$$\frac{\lambda_1 + \ldots + \lambda_K}{\lambda_1 + \ldots + \lambda_p}$$

注意前$K$个分量解释的方差比例随着$K$的增加单调增加,当$K = p$时达到理论最大值1。

虽然解释方差比例越高越好,也意味着$K$越大,数据中的信息尽可能被保留,但是在真实数据中,往往在$K$很小的时候就已经可以解释很大比例的方差。例如$K = 5$时,可解释的方差比例为90%;$K = 10$时,可解释的方差比例为99%;而$K = 100$时,可解释的方差比例为100%。此时,结合业务需求选择$5 \leqslant K \leqslant 10$就已经能够在不损失太多数据信息的前提下获得$K$较小的好处。

Python实现

import numpy as np
from sklearn import decomposition
from matplotlib import pyplot as plt

# 生成数据
alpha = np.pi / 10
thetas = np.arange(0, np.pi * 2, np.pi / 30)
xs = 10 * np.cos(thetas)
ys = 2 * np.sin(thetas)
xs, ys = np.cos(alpha) * xs - np.sin(alpha) * ys, np.cos(alpha) * ys + np.sin(alpha) * xs
X = np.c_[(xs, ys)]

# 创建计算主成分分析的对象
pca_computer = decomposition.PCA(n_components=1)
# 进行主成分的计算并返回主成分系数
X_transformed = pca_computer.fit_transform(X)
# 解释方差比例
explained_var_ratio = pca_computer.explained_variance_ratio_[-1]
print('解释方差的比例为%.3f' % (explained_var_ratio, ))
# 主成分
principal_components = pca_computer.components_
# 计算数据投影到第一个主成分的坐标
X_projected = X_transformed.dot(principal_components)

plt.close()
plt.plot(X[:, 0], X[:, 1], '.', label='original data')
plt.plot(X_projected[:, 0], X_projected[:, 1], '.', label='projected data')
plt.xlim((-10, 10))
plt.ylim((-10, 10))
plt.legend()
plt.savefig('pca-example.png')
        
输出:

解释方差的比例为0.962
	

图中展示了原始的二维数据,和进行在第一个主成分上的投影。

我们再看一个例子

import numpy as np
from sklearn import decomposition
from matplotlib import pyplot as plt
from scipy.stats import ortho_group

np.random.seed(10)
# 数据量为n
n = 100
# 维度为p
p = 15
# 有效维度为K_true
K_true = 5
# 生成随机正交矩阵
U = ortho_group.rvs(dim=n)
V = ortho_group.rvs(dim=p)
# 生成对角阵
lambdas = np.random.uniform(-2, 0, p)
lambdas[:K_true] = np.random.uniform(2, 4, K_true)
lambdas = sorted(lambdas, reverse=True)
lambdas = np.exp(lambdas)
D = np.zeros((n, p))
D[:p, ] = np.diag(lambdas)
# 生成数据
X = U.dot(D).dot(V.T)

# 创建主成分分析的对象
pca_computer = decomposition.PCA(n_components=p)
# 进行主成分的计算
pca_computer.fit(X)
# 解释方差比例
evr_numerical = np.cumsum(pca_computer.explained_variance_ratio_)
print('根据数值结果,前K个分量对方差的解释比例为%s' % (evr_numerical, ))
# 理论值
var_tot = np.sum(lambdas ** 2)
var_components = np.cumsum(lambdas ** 2)
evr_theoretical = var_components / var_tot
print('根据理论,前K个分量对方差的解释比例为%s' % (evr_theoretical, ))

plt.close()
plt.plot(np.arange(1, p + 1), evr_numerical * 100.0, '-')
plt.xlabel('K')
plt.ylabel('Percentage of variance explained (%)')
plt.savefig('pca-example2.png')
        
输出:

根据数值结果,前K个分量对方差的解释比例为[0.40891091 0.75343889 0.93729569 0.97861421 0.99938924 0.99955085
 0.99968549 0.9997985  0.99986364 0.99990803 0.99994931 0.99996767
  0.99998132 0.99999445 1.        ]
根据理论,前K个分量对方差的解释比例为[0.40970147 0.75386118 0.93705987 0.97861547 0.99938809 0.99954917
 0.99968457 0.9997971  0.9998632  0.9999077  0.99994893 0.99996741
  0.99998115 0.99999429 1.        ]
	

可以看到对前$K$个分量解释的方差比例的理论计算结果和实际数值结果相吻合。

图中展示了前$K$个分量解释的方差比例关于$K$的折线图,可以看到当$K$超过5时,所能解释的方差比例增加较小,这意味着数据的大部分信息蕴藏在前5个主成分分量之中。

常见面试问题

Q:给定数据$X$,数值上如何求主成分分析?

Q:主成分分析和奇异值分解(singular value decomposition)有什么关系?

Q:如何选择$K$?

更多机器学习相关问题见本网站论坛机器学习理论版面机器学习实践版面

更多面试问题见面试真题汇总

想要查看更多数据科学相关的内容请关注我们的微信公众号知乎专栏