主成分分析PCA降维方法

主成分分析(principal component analysis)是无监督学习方法中一种常用于数据降维的方法。它通过正交线性变换把原始数据变换到新的坐标空间中，并保留少数几个方差最大的正交(不相关)坐标维度作为主成分，从而在较少的维度中保留原始数据绝大多数的信息。本文简述PCA的相关理论，并用numpy实现样本主成分分析方法。

降维即将庞大复杂的高维数据转换到低维空间进行表示，并且要求能够保留原始数据中的大部分信息。由于PCA将尽量保留最多的方差信息作为降维的目标，而方差表示随机变量的离散程度的衡量，因此其比较适用于发现数据中的基本结构，即数据中变量之间的关系，是数据分析的有力工具，常用于其他机器学习方法的前处理。

对于高维样本数据来说，可以通过计算高维随机变量的协方差矩阵对应的特征值和特征向量，找到方差最大的若干个主成分维度。具体地，对于总体主成分，有如下定理：

设 \(x\) 是 \(m\) 维随机变量，\(\Sigma\) 是 \(x\) 的协方差矩阵，\(\Sigma\) 的特征值分别是 \(\lambda_1 \ge \lambda_2 \ge \cdot \cdot \cdot \ge \lambda_m \ge 0\)，特征值对应的单位特征向量分别是 \(\alpha_1, \alpha_2, \cdot \cdot \cdot, \alpha_m\)，则 \(x\) 的第 \(k\) 主成分是

$$y_k = \alpha_k^T x = \alpha_{1k} x_1 + \alpha_{2k} x_2 + \cdot \cdot \cdot + \alpha_{mk} x_m, \ \ \ k=1,2,\cdot \cdot \cdot,m$$

\(x\) 的第 \(k\) 主成分的方差是

$$var(y_k) = \alpha_k^T \Sigma \alpha_k = \lambda_k, \ \ \ k=1,2,\cdot \cdot \cdot,m$$

即协方差矩阵 \(\Sigma\) 的第 \(k\) 个特征值。

由以上定理可知，为了求解主成分，仅需求解原始数据的协方差矩阵前几个最大的特征值对应的特征向量即可。这些特征向量即为从原始数据到主成分的线性转换矩阵。

接下来使用numpy实现PCA的协方差矩阵特征值分解算法。

读取需要降维的高维数组。注意在我的代码中，每一列代表一个高维样本数据，每一行代表一个数据维度。

import numpy as np

# 输入数据，每一列代表一个样本，每一行代表一个数据维度。此处4个维度x6个样本(m*n)
X = [[65, 70, 80, 79, 78, 66],
     [96, 87, 97, 86, 99, 77],
     [69, 63, 65, 66, 68, 64],
     [78, 82, 65, 55, 81, 60]]
X = np.array(X)

将样本数据进行规范化处理，即令样本中每个维度的数据均值为 \(0\)，方差为 \(1\)。具体地，对样本矩阵作如下变换：

$$x_{ij}^* = \frac{x_{ij} – \bar x_i}{\sqrt{s_{ii}}}, \ \ \ i=1,2,\cdot \cdot \cdot,m; \ \ \ j=1,2,\cdot \cdot \cdot,n$$

其中 \(x_{ij}^*\) 指规范化后的样本点，\(\bar x_i\)为第 \(i\) 个维度所有样本的均值，\(s_{ii}\) 为第 \(i\) 个维度所有样本的无偏方差（自由度为样本个数 \(n-1\)）。

在numpy代码中，需要指定axis=1，且在计算标准差时指定自由度ddof=1。

# 样本规范化（令每一个维度均值为0，方差为1）
xi_mean = np.mean(X, axis=1, keepdims=True)
xi_std = np.std(X, axis=1, ddof=1, keepdims=True)
x = (X - xi_mean) / xi_std

接下来计算规范化后样本的协方差矩阵 \(R\)，仍然注意方差计算时需要考虑自由度来进行无偏估计。

$$R = \frac{1}{n-1}XX^T$$

代码中我们既可以按照上述公式直接进行计算，也可以调用np.cov()方法，返回的结果一致。方法中rowvar参数指定是否用行表示维度，bias指定计算时采用无偏估计。

# 求协方差矩阵R(相关矩阵) (m*m)
cov_mat = np.cov(x, rowvar=True, bias=False)  # R = x @ x.T * (1/(n-1))

求解协方差矩阵的特征值和单位特征向量。这一步也很简单，numpy给我们提供了现成的方法，该方法返回两个值，分别为特征值和对应的单位特征向量矩阵（一列表示一个单位特征向量）。

# 求R的特征值和单位特征向量
eigVals, eigVects = np.linalg.eig(cov_mat)

下面将特征值按照从大到小的顺序进行排序，同时需要获取排序后对应的特征向量。采用np.argsort方法返回升序排序好的特征值索引，用np.flipud方法反转顺序使其成为降序，最后用它作为切片获取排序好的特征值和单位特征向量。

# 将特征值由大到小排序，求排序后特征值对应的单位特征向量
eigVal_indice = np.flipud(np.argsort(eigVals))
eigVals = eigVals[eigVal_indice]
eigVects = eigVects[:, eigVal_indice]

下面我们引入方差贡献率的概念，将其作为具体选择主成分数量 \(k\) 的依据。

第 \(k\) 主成分 \(y_k\) 的方差贡献率定义为 \(y_k\) 的方差与所有方差之和的比，记作 \(\eta_k\)

$$\eta_k = \frac{\lambda_k}{\sum_{i=1}^m \lambda_i}$$

\(k\) 个主成分 \(y_1,y_2,\cdot \cdot \cdot,y_k\) 的累计方差贡献率定义为 \(k\) 个方差之和与所有方差之和的比

$$\sum_{i=1}^k \eta_i = \frac{\sum_{i=1}^k \lambda_i}{\sum_{i=1}^m \lambda_i}$$

通常取 \(k\) 使得累计方差贡献率达到规定的百分比以上，如75%以上。累计方差贡献率反映了主成分保留信息的比例。

代码实现中，首先设定累计方差贡献率的预定值，累加计算每个主成分的贡献率直到满足预定值的条件，最终k表示需要的主成分个数，var_contri_ratio为主成分的方差贡献率列表。

# 设定累计方差贡献率阈值
var_perc = 0.75

# 求累计方差贡献率达到预定值的主成分个数k
eigVals_sum = np.sum(eigVals)
tmpSum = 0
k = 0
var_contri_ratio = []
for i in eigVals:
    var_contri_ratio.append(i/eigVals_sum)
    tmpSum += i
    k += 1
    if tmpSum >= eigVals_sum * var_perc:
        break
# 前k个主成分的累计方差贡献率
var_contri_ratio.append([tmpSum / eigVals_sum])
print("各主成分方差贡献率和累计方差贡献率:", var_contri_ratio)

各主成分方差贡献率和累计方差贡献率: [0.5040585481009499, 0.3111015984253003, [0.8151601465262501]]

此处得到 \(k=2\) 接下来便可以轻易地筛选出前 \(k\) 个单位特征向量，即各主成分对应的线性转换矩阵，每个主成分 \(y_i = \alpha_i^T x\)，\(\alpha_i\) 为第 \(i\) 个单位特征向量。

# 获取前k个特征向量
k_eigVects = eigVects[:, :k]

最终将线性转换矩阵应用到原始规范化样本数据中，得到最终结果：样本主成分 \(y = \alpha^T x\)，\(\alpha\) 为 \(k\) 个单位特征向量组成的特征向量矩阵。

# 求出k个样本主成分
y = k_eigVects.T @ x
print(y)

[[-1.30072829  0.56250854 -0.31497481  0.73513588 -1.76225232  2.080311  ]
 [ 1.27239194  0.92472781 -1.23198346 -1.33561276 -0.24930643  0.61978289]]

可见原始数据被成功降维成了仅包含两维的低维数据。

我们也可以通过降维后的数据重构出对应的高维复杂数据，但由于主成分数量小于原始数据维度，重构的数据包含信息损失，与原始数据并不一致。

# 反转换重构原始矩阵
print((k_eigVects @ y) * xi_std + xi_mean)

[[ 67.41286447  67.07746901  80.3539056   79.64299561  76.58562861  66.9271367 ]
 [ 95.97399406  85.82632544  93.85818027  88.0096006  100.74144814  77.5904515 ]
 [ 67.61526013  65.18997936  66.13196194  64.77305031  68.07729051  63.21245775]
 [ 84.92243683  72.91922638  64.20280608  58.00950079  77.93928044  63.00674948]]

至此所有工作全部完成，下面附上sklearn中封装的PCA方法的调用代码。sklearn默认以一行代表一个样本，为了使输出结果一致，在调用时应用了转置矩阵。需要注意的是，其仅将原始矩阵均值规范化为 \(0\)，未对方差进行规范化处理，这里我直接传入了已严格规范化好的矩阵，因此结果没有影响。

# 使用sklearn验证（其默认仅归一化均值为0，未将方差设为1）
from sklearn.decomposition import PCA
pca = PCA(n_components=k)  # 指定保留的主成分数量
sk_y = pca.fit_transform(x.T)
print(sk_y.T)  # 输出主成分
print(pca.explained_variance_ratio_)  # 输出主成分方差贡献率
print(pca.inverse_transform(sk_y).T * xi_std + xi_mean)  # 输出重构矩阵

输出的主成分、方差贡献率和重构矩阵与之前的结果完全一致。