1. 简介

主成分分析是指将数据中相关性很高的属性/变量转化成彼此相互独立或不相关的新属性/变量，利用较少的新属性/变量（主成分）去解释原来数据中的大部分属性/变量的一种降维方法。

2. 原理

2.1 基础思想

以学校课程为例，用 ${x_1,x_2,\cdots,x_p}$ 表示 p 门课程的成绩变量， ${c_1,c_2,\cdots,c_p}$ 表示各门课程的权重，则成绩的加权和为

$$\begin{array}{c} s = c_1x_1 + c_2x_2 + \cdots + c_px_p \end{array}$$

现在希望能够选择合适的权重来更好的区分学生的成绩 ${s_1,s_2,\cdots,s_n}$（ $n$ 为学生人数且 ${n \gt p}$）。一般来说，如果这些值很分散，则表明区分的好。故需要寻找这样的加权，使得 ${s_1,s_2,\cdots,s_n}$ 能够尽可能的分散。

而数据的分散程度可以利用方差的大小来体现，设 ${X_1,X_2,\cdots,X_p}$ 表示以 ${x_1,x_2,\cdots,x_p}$ 为样本观测值的随机变量，规定（否则权值可以无穷大而失去意义） ${c_1^2 + c_2^2 + \cdots + c_p^2 = 1}$ 。在此约束下，若能找到 ${c_1,c_2,\cdots,c_p}$ 使得 ${Var(c_1X_1 + c_2X_2 + \cdots + c_pX_p)}$ 的值达到最大，则表明此时学生的成绩区分达到最好。

可见，以上问题本质为一个求最优解的优化问题：

$$\begin{array}{c} c_1^2 + c_2^2 + \cdots + c_p^2 = 1 \\ \max \ Var(c_1X_1 + c_2X_2 + \cdots + c_pX_p) \end{array}$$

求出的最优解 ${c_1,c_2,\cdots,c_p}$ 是一个单位向量，代表一个主成分方向，对应的主成分为 ${Z = c_1X_1 + c_2X_2 + \cdots + c_pX_p}$ 。
一个主成分不足以代表原来的 $p$ 个变量，故需要寻找第二、三、··· 个主成分。且后面的主成分不应该再包含前面已经求出来的主成分的信息，统计上理解就是主成分间两两之间的协方差为零，几何上理解就是主成分间两两正交。
设 ${Z_i}$ 表示第 $i$ 个主成分，${i = 1,2,\cdots,p}$，则

$$\begin{array}{c} Z_1 = c_{11}X_1 + c_{12}X_2 + \cdots + c_{1p}X_p \\ Z_2 = c_{21}X_1 + c_{22}X_2 + \cdots + c_{2p}X_p \\ \vdots \\ Z_p = c_{p1}X_1 + c_{p2}X_2 + \cdots + c_{pp}X_p \end{array}$$

其中，对于每一个 $i$，均有 ${c_{i1}^2 + c_{i2}^2 + \cdots + c_{ip}^2 = 1}$。且 ${(c_{11},c_{12},\cdots,c_{1p})}$ 使得 ${Var(Z_1)}$ 的值达到最大；${(c_{21},c_{22},\cdots,c_{2p})}$ 不仅垂直于 ${(c_{11},c_{12},\cdots,c_{1p})}$，而且使得 ${Var(Z_2)}$ 的值达到最大；${(c_{31},c_{32},\cdots,c_{3p})}$ 不仅垂直于 ${(c_{11},c_{12},\cdots,c_{1p})}$和${(c_{21},c_{22},\cdots,c_{2p})}$，而且使 ${Var(Z_3)}$ 的值达到最大；以此类推得到全部 $p$ 个主成分。

【注意事项】

• 主成分分析的结果受量纲的影响，故实际应用中先将各变量的数据标准化，然后使用协方差矩阵或相关系数矩阵进行分析。
• 使用相关系数矩阵求主成分时，Kaiser 主张将特征值小于 1 的主成分予以放弃（这也是 SPSS 软件的默认值）
• 在实际研究中，由于主成分的目的是为了降维，减少变量个数，故一般选取少量的主成分（不超过 5 或 6 个），只要这些主成分的累计贡献率（定义见下文）能够达到 70%~80% 就行了。

设 $p$ 个变量 ${X_1,X_2,\cdots,X_p}$ 在第 $i$ 次试验中的取值为 ${x_{i1},x_{i2},\cdots,x_{ip} (i = 1,2,\cdots,n)}$，则原数据写成矩阵形式为

$$\begin{array}{c} X = \left( \begin{matrix} x_{11} & x_{12} & \cdots & x_{1p} \\ x_{21} & x_{22} & \cdots & x_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ x_{n1} & x_{n2} & \cdots & x_{np} \end{matrix} \right) \end{array}$$

为了保证主成分分析不受量纲的影响，将 $X$ 进行标准化变换

$$\begin{array}{c} \tilde{x_{ij} } = (x_{ij}-\bar{x_{j} })/s_j \end{array}$$

其中，$\bar{x_j}$、$s_j$ 分别为 $X$ 的第 $j$ 列的均值和标准差。
变换后的数据矩阵为

$$\begin{array}{c} \tilde{X} = \left( \begin{matrix} \tilde{x_{11} } & \tilde{x_{12} } & \cdots & \tilde{x_{1p} } \\ \tilde{x_{21} } & \tilde{x_{22} } & \cdots & \tilde{x_{2p} } \\ \vdots & \vdots & \ddots & \vdots \\ \tilde{x_{n1} } & \tilde{x_{n2} } & \cdots & \tilde{x_{np} } \end{matrix} \right) \end{array}$$

2.2 确定主成分变量

对于自变量的任意一个线性组合 ：

$$\begin{array}{c} z = c_1\tilde{x_1} + c_2\tilde{x_2} + \cdots + c_p\tilde{x_p}, \sum_{j=1}^{p}c_j^2 = 1 \end{array}$$

其在第 $i$ 次试验中的取值为

$$\begin{array}{c} z_i = c_1\tilde{x_{i1} } + c_2\tilde{x_{i2} } + \cdots + c_p\tilde{x_{ip} }\ (i = 1,2,\cdots,n) \end{array}$$

由于 $\tilde{X}$ 已经标准化了，因此

$$\begin{array}{c} \bar{z}={1 \over n}\sum_{i=1}^{n}z_i = {1 \over n}\sum_{i=1}^{n}\sum_{j=1}^{p}c_j\tilde{x_{ij} } = {1 \over n}\sum_{j=1}^{p}c_j\sum_{i=1}^{n}\tilde{x_{ij} } = 0 \end{array}$$

记 ${w = (c_1,c_2,\cdots,c_p)^T}$，则

$$\begin{array}{c} M^* = Var(z) = {1 \over n}\sum_{i=1}^{n}(z_i - \bar{z})^2 = {1 \over n}\sum_{i=1}^{n}z_i^2 = {1 \over n}(\tilde{X}w)^T(\tilde{X}w) \end{array}$$

对于新变量 ${z}$ 来说，如果在 $n$ 次试验下它的取值变化不大，即 ${z}$ 的 ${M^*}$ 较小，则这个新变量 ${z}$ 可以去掉。反之，如果 ${z}$ 的 ${M^*}$ 较大，则说明该新变量 ${z}$ 的作用比较明显。因此，我们希望所选择的 ${c_i\ (i=1,2,\cdots,p)}$ 能使 ${M^*}$ 达到最大，即使得新变量 $z$ 有较大作用。

如果 ${X^TX}$ 的特征值 ${\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_p}$，它们对应的标准化正交特征向量为 ${\eta_1,\eta_2,\cdots,\eta_p}$，则 ${M^* = (Xw)^T(Xw)/n}$ 的最大值在 ${w = \eta_1}$ 时达到，且最大值为 ${\lambda_1 / n}$ 。此时新变量 ${z}$ 为

$$\begin{array}{c} z = \tilde{x}^T\eta_1 \end{array}$$

其中 ${\tilde{x} = (\tilde{x_1},\tilde{x_2},\cdots,\tilde{x_p})}$。

常记 ${Z_1 = \tilde{x}^T\eta_1}$ ，并称 ${Z_1}$ 为变量的第一主成分。
一般的，如果已经确定了 $k$ 个主成分：

$$\begin{array}{c} Z_i = \tilde{x}^T\eta_i\ \ (i = 1,2,\cdots,k) \end{array}$$

则第 ${k+1}$ 个主成分 ${Z_{k+1} = \tilde{x}^Tw}$ 可由以下两个条件决定：

• ${w^T\eta_i = 0,\ i = 1,2,\cdots,k,\ w^Tw = 1}$；
• 在以上条件下，使得 ${M^*}$ 达到最大。

由二次型的条件极值可知，第 ${k+1}$ 个主成分就是 ${Z_{k+1} = \tilde{x}^T\eta_{k+1}\ \ (i = 1,2,\cdots,p)}$ 。如此总共可以取到 $p$ 个主成分 ${Z_i = \tilde{x}^T\eta_i \ (i = 1,2,\cdots,p)}$ 。

2.3 选择主成分变量

确定了 $p$ 个主成分后，将标准化后的数据 ${\tilde{x_1},\tilde{x_2},\cdots,\tilde{x_p} }$ 转换成主成分 ${z_1,z_2,\cdots,z_p}$，令

$$\begin{array}{c} Z = \tilde{X}(\eta_1,\eta_2,\cdots,\eta_p) = \left( \begin{matrix} z_{11} & z_{12} & \cdots & z_{1p} \\ z_{21} & z_{22} & \cdots & z_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ z_{n1} & z_{n2} & \cdots & z_{np} \end{matrix} \right) \end{array}$$

记 ${Q = (\eta_1,\eta_2,\cdots,\eta_p)_{p \times p} }$，$Q$ 为标准正交阵，且 ${Z = XQ}$，则

$$\begin{array}{c} Z^TZ = Q^TX^TXQ = Q^T(X^TX)Q = \left( \begin{matrix} \lambda_1 & & 0 \\ & \ddots & \\ 0 & & \lambda_p \end{matrix} \right) \end{array}$$

由上式可知，${X^TX}$ 的特征值 ${\lambda_i}$ 度量了第 $i$ 个主成分 ${z_i}$ 在第 $n$ 试验中取值变化的大小。如果 ${\lambda_i \approx 0}$，则说明该主成分 ${z_i}$ 在 $n$ 次实验中的取值变化很小，可以选择剔除该主成分。

【筛选方法】

• 首先将 ${X^TX}$ 的特征值按由大到小的次序排序。
• 然后剔除这些特征值 ${\lambda_{r+1},\lambda_{r+2},\cdots,\lambda_{p} }$对应的主成分，这些特征值满足

$$\begin{array}{c} \sum_{i=r+1}^{p}\lambda_i \lt 15\%\sum_{i=1}^{p}\lambda_i \end{array}$$

即筛选留下的主成分对应的特征值之和所占比重（定义为累计贡献率）要超过 ${85\%}$。

• 单纯考虑累计贡献率有时是不够的，还需考虑选择的主成分对原始变量的贡献值，定义贡献值为相关系数的平方和。如果选择的主成分为 ${z_1,z_2,\cdots,z_p}$ ，则它们对原始变量 ${\tilde{x_i} }$ 的贡献值为

$$\begin{array}{c} \rho_i = \sum_{j=1}^{r}r^2(z_j,x_i) \end{array}$$

这里 ${r(z_j,x_i)}$ 表示 ${z_j}$ 与 ${x_i}$ 的相关系数。

3. 步骤

假设进行主成分分析的指标变量有 $m$ 个：${x_1,x_2,\cdots,x_m}$；共有 $n$ 个评价对象，第 $i$ 个评价对象的第 $j$ 个指标的取值为 ${a_{ij} }$，评价矩阵为

$$\begin{array}{c} A = (a_{ij})_{n \times m} \end{array}$$

3.1 对原始数据进行标准化

将各指标值转换成标准指标 ${\tilde{a_{ij} }}$ ：

$$\begin{array}{c} \tilde{a_{ij} } = { {a_{ij} - \mu_j} \over s_j}\ (i=1,2,\cdots,n;\ j=1,2,\cdots,m) \end{array}$$

其中

$$\begin{array}{c} \mu_j = {1 \over n}\sum_{i=1}^na_{ij} \\ s_j = {1 \over {n-1} }\sum_{i=1}^n(a_{ij}-\mu_j)^2 \\ j = 1,2,\cdots,m \end{array}$$

即 ${\mu_j}$，${s_j}$ 为第 $j$ 个指标的标准均值和样本标准差。对应地，称

$$\begin{array}{c} \tilde{x_i}={ {x_i-\mu_i} \over s_i} \ \ (i = 1,2,\cdots,m) \end{array}$$

为标准化指标变量。标准化后的评价矩阵为

$$\begin{array}{c} \tilde{A} = (\tilde{a_{ij} })_{n \times m} \end{array}$$

3.2 计算相关系数矩阵 $R$

相关系数矩阵 ${R=(r_{ij})_{m \times m} }$，其中

$$\begin{array}{c} r_{ij} = { {\sum_{k=1}^n\tilde{a_{ki} } \cdot \tilde{a_{kj} } } \over {n-1} }\ \ (i,j=1,2,\cdots,m) \end{array}$$

易知式中 ${r_{ii} = 1}$，${r_{ij} = r_{ji} }$，${r_{ij} }$ 是第 $i$ 个指标与第 $j$ 个指标的相关系数。

3.3 计算特征值和特征向量

计算相关系数矩阵 $R$ 的特征值 ${\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_m \geq 0}$，及对应的特征向量 ${\mu_1,\mu_2,\cdots,\mu_m}$，其中 ${\mu_j = (\mu_{1j},\mu_{2j},\cdots,\mu_{mj})^T}$，由特征向量组成 $m$ 个新的指标变量

$$\begin{array}{c} y_1 = \mu_{11}\tilde{x_1} + \mu_{21}\tilde{x_2} + \cdots + \mu_{m1}\tilde{x_m} \\ y_2 = \mu_{12}\tilde{x_1} + \mu_{22}\tilde{x_2} + \cdots + \mu_{m2}\tilde{x_m} \\ \vdots \\ y_m = \mu_{1m}\tilde{x_1} + \mu_{2m}\tilde{x_2} + \cdots + \mu_{mm}\tilde{x_m} \end{array}$$

其中 ${y_i}$ 表示第 $i$ 主成分。即

$$\begin{array}{c} Y = \tilde{X}^T\mu \end{array}$$

其中，${\tilde{X} = (\tilde(x_1),\tilde{x_2},\cdots,\tilde{x_m})^T}$，${Y = (y_1,y_2,\cdots,y_m)^T}$，${\mu = (\mu_1,\mu_2,\cdots,\mu_m)}$ 。
标准化后的评价矩阵 ${\tilde{A} }$ 转化为新评价矩阵 ${\tilde{B} }$：

$$\begin{array}{c} \tilde{B} = \tilde{A}\mu \end{array}$$

3.4 选择 ${p \ (p \leq m)}$ 个主成分

计算特征值 ${\lambda_j \ (j = 1,2,\cdots,m)}$ 的信息贡献率和累计贡献率。称

$$\begin{array}{c} b_j = {\lambda_j \over {\sum_{k=1}^{m}\lambda_k} } \ \ (j = 1,2,\cdots,m) \end{array}$$

为主成分 ${y_j}$ 的信息贡献率；称

$$\begin{array}{c} \alpha_p = { {\sum_{k=1}^{p}\lambda_k} \over {\sum_{k=1}^m\lambda_k} } \end{array}$$

为主成分 ${y_1,y_2,\cdots,y_p}$ 的累计贡献率，当 ${\alpha_p}$ 接近1（${\alpha_p = 0.85 / 0.90 / 0.95}$）时，选择前 $p$ 个指标变量 ${y_1,y_2,\cdots,y_p}$ 作为 $p$ 个主成分，代替原来的 $m$ 个指标变量。选择该 $p$ 个主成分后，新评价矩阵 $\tilde{C}$ 为：

$$\begin{array}{c} \tilde{C} = \tilde{B}_{(n \times p)} \end{array}$$

3.5 计算综合评价值

定义综合得分为

$$\begin{array}{c} z = \sum_{j=1}^pb_jy_j \end{array}$$

其中 ${b_j}$ 为第 $j$ 个主成分的信息贡献；则 $n$ 个评价对象的综合评价值向量为

$$\begin{array}{c} Z = \tilde{C}b \end{array}$$

其中，${b = (b_1,b_2,\cdots,b_p)^T}$ 。
然后便可以依据该评价值向量对这 $n$ 个对象进行综合评价。

4. 代码实现

4.1 MatLab 代码实现如下

• 利用 pcacov 函数：

% 1. 数据标准化
data=zscore(data);
% 2. 计算相关系数矩阵
R=corrcoef(data);
% 3. PCA 分析
% x：特征向量矩阵；y：特征值向量；z：主成分贡献率向量（总和为 100 ）
[x,y,z]=pcacov(R);
% 4. 选择 5 个主成分
p = 5;
% 5. 计算新评价矩阵
C = data*x;
C = C(:,1:p);
% 6. 计算综合评价值
Z = C*z(1:p)/100;

• 利用 pca 函数：

% 1. 数据标准化
data=zscore(data);
% 2. PCA 分析
% x：特征向量矩阵；C：新评价矩阵；y：特征值向量
[x,C,y]=pca(data);
% 3. 选择 5 个主成分
p = 5;
% 4. 计算主成分贡献率（总和为 1 ）
z = y/sum(y);
% 5. 计算综合评价值
Z = C*z(1:p);