课堂上听老师讲PCA，总觉得思路特别清楚：标准化、算协方差、求特征值特征向量。可一旦自己动手算，尤其是要从一堆矩阵运算里把结果一步步推出来，经常对不上答案。后来我重新把PCA当成线性代数题来做，才意识到，所谓计算步骤其实是一条完整的矩阵链条：从数据矩阵出发，得到中心化矩阵，再到协方差矩阵，最后解特征方程、拼出载荷矩阵、求得主成分得分。这篇文章就把这条链路上的每一个环节都拆开讲，配一个能从头算到尾的数值例子，帮你把PCA的计算步骤真正吃透。适合正在复习线性代数、准备数据分析考试、或者刚接触多元统计打算搞懂原理的读者。

1. 把PCA翻译成一件事：求协方差矩阵的特征基

PCA的计算步骤，从数学上看非常统一。它的目标不是“降维”这个动作本身，而是找到一组新的正交坐标轴，让数据在这些轴上的投影方差依次达到最大。如果把它拆成线性代数语言，每个主成分其实就是数据协方差矩阵的一个特征向量，而主成分的方差就是对应的特征值。

1.1 行样本、列变量的约定

先约定记号。设原始数据矩阵为 X，维度是 n×p，n 代表样本数，p 代表变量数，第 i 行是第 i 个样本，第 j 列是第 j 个变量。这个约定和大多数机器学习库一致：sklearn、R 的 prcomp、Python 里手写矩阵运算，都默认行是样本、列是特征。

有的教材习惯把数据写成 p×n，变量在行、样本在列，这会导致协方差矩阵的表示不同。你做 PCA 前第一步要做的，不是写公式，而是先明确自己的 X 到底是哪个方向。否则后面所有矩阵乘法方向都会混乱。

记列均值向量为：

[ \bar{x} = \frac{1}{n} X^T \mathbf{1}_n ]

(\mathbf{1}_n) 是 n 维全 1 列向量。那么把数据中心化，得到：

[ X_c = X - \mathbf{1}_n \bar{x}^T ]

这个记号看起来很抽象，实际就是每一列减去该列的均值。用矩阵视角看，均值中心化是在列方向施加一个投影矩阵：

[ H = I_n - \frac{1}{n} \mathbf{1}_n \mathbf{1}_n^T ]

[ X_c = HX ]

(H) 是一个对称幂等矩阵，也叫中心化投影矩阵，它把数据每列的均值“归零”。这个公式在推导协方差矩阵时很好用，因为 (X_c) 的任何线性变换仍保持零均值。

1.2 主成分方向就是最大方差方向

假设我们想找第一个主成分方向，用一个单位向量 (w) 表示，p 维。所有样本投影到 (w) 上的得分就是：

[ z = X_c w ]

这个投影向量的样本方差是：

[ \mathrm{Var}(z) = \frac{1}{n-1} w^T X_c^T X_c w = w^T S w ]

其中：

[ S = \frac{1}{n-1} X_c^T X_c ]

是 p×p 的样本协方差矩阵。所以问题变成：

[ \max_{w, ; w^T w = 1} w^T S w ]

用拉格朗日乘子法，求导可得：

[ S w = \lambda w ]

也就是说，使投影方差最大的方向 (w)，恰好是协方差矩阵 (S) 的特征向量，而最大方差值就是对应的特征值 (\lambda)。第二个主成分方向，需要与第一个方向正交，且同样最大化投影方差；由于 (S) 是对称矩阵，不同特征值的特征向量天然正交，所以按特征值从大到小取前 k 个特征向量即可。

这里的核心结论是：PCA 的计算过程，本质上就是求 (S) 的特征值分解。后面所有步骤，都是围绕“如何从 X 算出 S、如何分解 S”展开的。

2. 动手前的关键一步：中心化与标准化不能省

很多新手拿到数据就直接求协方差矩阵，发现计算结果和软件对不上，问题往往出在预处理上。中心化是必须的，标准化则要视情况而定。

2.1 中心化的矩阵实现

中心化就是让每个变量的均值变成 0。比如数据 X 是：

变量1的均值是 3，变量2的均值是 4，于是：

[ X_c = X - \begin{bmatrix} 1 \ 1 \ 1 \ 1 \ 1 \end{bmatrix} \begin{bmatrix} 3 & 4 \end{bmatrix} ]

结果就是每列都减去对应均值：

中心化后的数据，每一列的均值都是 0，这是一个非常实用的检查点。你任何时候手算完中心化，都应该先列求和验证：如果结果不是 0，说明减错了。

用前面提到的投影矩阵 (H) 同样能得到 (X_c)，两种方式等价。但在推导 PCA 的一些性质时，例如证明主成分得分均值为 0，用 (H) 更简洁。

2.2 标准化到底在做什么

标准化是在中心化之后，每一列再除以该列的标准差。设每列标准差为 (s_j)，定义对角矩阵：

[ D = \mathrm{diag}(s_1, s_2, \dots, s_p) ]

标准化后的数据矩阵为：

[ X_{std} = X_c D^{-1} ]

这其实是对数据列方向做了一次伸缩变换。为什么要做这一步？因为协方差矩阵对变量的量纲极其敏感。假设第一个变量单位是米，取值在个位数级别，第二个变量单位是毫米，取值在几千的级别，那么第二个变量的方差会远远大于第一个变量，PCA 会把全部注意力放在它上面，第一个变量的信息几乎被忽略。

这往往不是你想要的。如果你希望每个变量在 PCA 中具备同等地位，就需要标准化，之后计算的其实是相关矩阵：

[ R = \frac{1}{n-1} X_{std}^T X_{std} ]

相关矩阵的对角线全是 1，特征值之和等于变量数 p。协方差矩阵的特征值之和则等于原始变量的方差之和，保留的是原始尺度信息。

2.3 什么时候用协方差矩阵，什么时候用相关矩阵

我的经验是：如果各变量的测量单位相同、量级接近，比如几个物理尺寸读数，直接用协方差矩阵；如果变量单位不同或量级差异大，比如身高、体重、收入、年龄混在一起，必须先标准化，用相关矩阵。

严格说，标准化的本质是人为地给每个变量相同的权重。这个选择没有绝对对错，取决于你对变量的态度。有些人会说“PCA 之前必须标准化”，这句话不准确；真实情况是“当你需要消除量纲影响时，标准化是必要步骤”。教材里的例子通常会用协方差矩阵，因为你跟着原数据计算特征值，能直观看到方差占比；实战中多为相关矩阵。

3. 特征值到手的方法：行列式、齐次方程组与单位化

中心化和标准化完成，数据矩阵变成 (X_c)（或 (X_{std})），接下来就是 PCA 计算的核心：对协方差矩阵做特征值分解。

3.1 第一步：写出协方差矩阵

对中心化矩阵 (X_c)，协方差矩阵是：

[ S = \frac{1}{n-1} X_c^T X_c ]

它是一个 (p \times p) 的对称方阵。以第 2 小节的例子为例：

[ X_c = \begin{bmatrix} -2 & -2 \ -1 & -1 \ 0 & 1 \ 1 & 0 \ 2 & 2 \end{bmatrix} ]

先算：

[ X_c^T X_c = \begin{bmatrix} 10 & 9 \ 9 & 10 \end{bmatrix} ]

再除以 (n-1=4)，得到：

[ S = \begin{bmatrix} 2.5 & 2.25 \ 2.25 & 2.5 \end{bmatrix} ]

这个矩阵的对称性是一个天然检查点：右上角和左下角必须相等。如果不等，几乎可以断定前面的数据矩阵构造出了问题。

3.2 第二步：解特征方程

要求特征值，解：

[ |S - \lambda I| = 0 ]

对于上面的 2×2 矩阵：

[ \begin{vmatrix} 2.5 - \lambda & 2.25 \ 2.25 & 2.5 - \lambda \end{vmatrix} = (2.5 - \lambda)^2 - 2.25^2 = 0 ]

展开：

[ \lambda^2 - 5\lambda + 1.1875 = 0 ]

解得：

[ \lambda_1 = 4.75, \quad \lambda_2 = 0.25 ]

对 2×2 的情况，用求根公式即可。对更高维矩阵，手算行列式会很痛苦，教材一般也只要求 3×3 以下。真实计算中不会展开行列式，而是采用数值迭代算法，比如幂迭代法或 QR 算法，但求“特征多项式→求根”的思路是一致的。

这里有一个非常有用的验证：特征值之和等于矩阵迹，即对角线元素之和。本例中 (4.75 + 0.25 = 5)，而 (\mathrm{trace}(S) = 2.5 + 2.5 = 5)，符合。这也是 PCA 中“累计解释方差”计算的基础：每个特征值除以所有特征值之和，就是该主成分解释的方差比例。

3.3 第三步：回代求特征向量

对每个特征值 (\lambda)，求满足下式的单位向量：

[ (S - \lambda I) v = 0 ]

这一步要解齐次线性方程组。以 (\lambda_1 = 4.75) 为例：

[ \begin{bmatrix} 2.5 - 4.75 & 2.25 \ 2.25 & 2.5 - 4.75 \end{bmatrix}

\begin{bmatrix} -2.25 & 2.25 \ 2.25 & -2.25 \end{bmatrix} ]

解：

[ -2.25 v_1 + 2.25 v_2 = 0 ]

得到 (v_1 = v_2)。也就是说，特征向量方向是 ((1, 1)) 的倍数。单位化后：

[ v_1 = \begin{bmatrix} 0.7071 \ 0.7071 \end{bmatrix} ]

对 (\lambda_2 = 0.25)：

[ \begin{bmatrix} 2.25 & 2.25 \ 2.25 & 2.25 \end{bmatrix} ]

解：

[ 2.25 v_1 + 2.25 v_2 = 0 ]

得到 (v_1 = -v_2)。单位化后选择：

[ v_2 = \begin{bmatrix} 0.7071 \ -0.7071 \end{bmatrix} ]

两个特征向量点积为 0，正交性成立。它们是 PCA 的载荷向量，也就是新的坐标轴方向。注意特征向量是可以同时乘以 -1 的，这不改变方向，所以符号正负没有本质影响。

3.4 第四步：按特征值降序排列，计算得分矩阵

特征值表示对应主成分的重要性，所以按 (\lambda_1 > \lambda_2) 排列。载荷矩阵 (V) 的每一列是一个特征向量：

[ V = \begin{bmatrix} 0.7071 & 0.7071 \ 0.7071 & -0.7071 \end{bmatrix} ]

主成分得分就是中心化数据在载荷矩阵上的投影：

[ Z = X_c V ]

Z 的第 j 列就是第 j 个主成分的得分。实际应用中，如果只保留前 k 个主成分，就取 V 的前 k 列。这个矩阵乘法是整个 PCA 计算链条的收尾，也是我下一节要完整演示的部分。

4. 一遍 5 个样本的完整计算：从中心化到得分矩阵

前面零散讲了每一步的原理，这节用一个完整例子把所有步骤串起来。数据还是第 2 节那 5 个样本、2 个变量。我特意选 5 个样本，是为了演示“除以 n-1”时不会出现尴尬的分数，同时样本量小，每个数字都能检验。

4.1 原始数据与均值中心化

原始数据：

均值向量：

[ \bar{x}_1 = 3, \quad \bar{x}_2 = 4 ]

中心化矩阵：

[ X_c = \begin{bmatrix} -2 & -2 \ -1 & -1 \ 0 & 1 \ 1 & 0 \ 2 & 2 \end{bmatrix} ]

验证列和：第一列 (-2-1+0+1+2=0)，第二列 (-2-1+1+0+2=0)，通过。

4.2 协方差矩阵的手算过程

计算 (X_c^T X_c) 时，逐个元素来：

[ S_{11} = \frac{(-2)^2 + (-1)^2 + 0^2 + 1^2 + 2^2}{4} = \frac{10}{4} = 2.5 ]

[ S_{12} = \frac{(-2)(-2) + (-1)(-1) + 0 \times 1 + 1 \times 0 + 2 \times 2}{4} = \frac{9}{4} = 2.25 ]

[ S_{21} = S_{12} ]

[ S_{22} = \frac{(-2)^2 + (-1)^2 + 1^2 + 0^2 + 2^2}{4} = \frac{10}{4} = 2.5 ]

所以：

[ S = \begin{bmatrix} 2.5 & 2.25 \ 2.25 & 2.5 \end{bmatrix} ]

两个变量的方差相同，都是 2.5，协方差是 2.25，说明呈正相关关系。协方差矩阵的迹为 5，总方差是 5。

4.3 特征值与解释方差占比

特征方程已经解过，简化如下：

[ (2.5 - \lambda)^2 - 2.25^2 = 0 ]

[ \lambda_1 = 4.75, \quad \lambda_2 = 0.25 ]

每个主成分解释的方差比例为：

[ \frac{4.75}{5} = 95%, \quad \frac{0.25}{5} = 5% ]

这意味着，只用第一个主成分就能保留原始数据 95% 的总方差。二维数据降到一维，损失只有 5%，是非常理想的情况。

4.4 特征向量与主成分载荷

前面已求出：

[ v_1 = \begin{bmatrix} 0.7071 \ 0.7071 \end{bmatrix}, \quad v_2 = \begin{bmatrix} 0.7071 \ -0.7071 \end{bmatrix} ]

载荷矩阵按特征值降序排列：

[ V = \begin{bmatrix} 0.7071 & 0.7071 \ 0.7071 & -0.7071 \end{bmatrix} ]

第一个主成分是变量1和变量2的等权平均方向，第二个主成分是两者之差的方向。之所以第一个方向是“平均方向”，是因为两个变量高度正相关，数据的整体差异主要体现在它们一起变大或变小；第二个方向体现的是变量之间的相对差异。

4.5 主成分得分矩阵

计算 (Z = X_c V)：

[ X_c V = \begin{bmatrix} -2 & -2 \ -1 & -1 \ 0 & 1 \ 1 & 0 \ 2 & 2 \end{bmatrix} \begin{bmatrix} 0.7071 & 0.7071 \ 0.7071 & -0.7071 \end{bmatrix} ]

逐个样本计算：

样本 A：

[ z_{A1} = -2 \times 0.7071 + (-2) \times 0.7071 = -2.8284 ]

[ z_{A2} = -2 \times 0.7071 + (-2) \times (-0.7071) = 0 ]

样本 B：

[ z_{B1} = -1 \times 0.7071 + (-1) \times 0.7071 = -1.4142 ]

[ z_{B2} = -1 \times 0.7071 + (-1) \times (-0.7071) = 0 ]

样本 C：

[ z_{C1} = 0 \times 0.7071 + 1 \times 0.7071 = 0.7071 ]

[ z_{C2} = 0 \times 0.7071 + 1 \times (-0.7071) = -0.7071 ]

样本 D：

[ z_{D1} = 1 \times 0.7071 + 0 \times 0.7071 = 0.7071 ]

[ z_{D2} = 1 \times 0.7071 + 0 \times (-0.7071) = 0.7071 ]

样本 E：

[ z_{E1} = 2 \times 0.7071 + 2 \times 0.7071 = 2.8284 ]

[ z_{E2} = 2 \times 0.7071 + 2 \times (-0.7071) = 0 ]

最终得分矩阵：

[ Z = \begin{bmatrix} -2.8284 & 0 \ -1.4142 & 0 \ 0.7071 & -0.7071 \ 0.7071 & 0.7071 \ 2.8284 & 0 \end{bmatrix} ]

4.6 用方差做反向验证

计算 Z 第一列的均值：

[ \frac{-2.8284 - 1.4142 + 0.7071 + 0.7071 + 2.8284}{5} = 0 ]

样本方差：

[ \frac{(-2.8284)^2 + (-1.4142)^2 + 0.7071^2 + 0.7071^2 + 2.8284^2}{4} = \frac{19}{4} = 4.75 ]

正好等于 (\lambda_1)。Z 第二列的样本方差：

[ \frac{0 + 0 + 0.7071^2 + 0.7071^2 + 0}{4} = 0.25 ]

正好等于 (\lambda_2)。到这一步，整个计算闭环了：主成分得分的方差等于对应特征值，这正是 PCA 要求的最大方差性质。

再加上正交性检查，两列得分的内积为 0，说明两个主成分完全不相关：

[ (-2.8284)(0) + (-1.4142)(0) + (0.7071)(-0.7071) + (0.7071)(0.7071) + (2.8284)(0) = 0 ]

5. 当计算结果和软件不一致时，先查这三处

手算与软件结果不一致，是第一遍练习里最常见的事。我见过不少同学在 PCA 作业里算到凌晨，最后发现不是思路错，而是下面三个细节之一出了问题。

5.1 特征向量符号正负不同

特征向量乘以 -1 后仍是特征向量。软件包之间、软件不同版本之间，特征向量的符号都可能不同。例如 (v_1 = (0.7071, 0.7071)) 和 ((-0.7071, -0.7071)) 表示同一个方向。你手算取正号，软件取负号，载荷矩阵的每一列可能整体变号，但主成分得分只是跟着整体变号，数据的投影位置完全不变。

这不算错误，但很多人会误以为自己算错了。判断方法很简单：看你的载荷向量和软件输出的载荷向量是否成倍数关系，且比例为 ±1。如果比例是 1 或 -1，就说明特征向量本质相同。

5.2 协方差矩阵的除数是 n 还是 n-1

样本协方差矩阵有两种写法：(\frac{1}{n-1}X_c^T X_c) 或 (\frac{1}{n}X_c^T X_c)。不同教材、不同软件采用不同约定。用 4 除还是 5 除，会改变特征值的具体数值，但特征向量的方向不变。

给一个具体对应关系：如果用 (\frac{1}{n})，本例特征值会变成 3.8 和 0.2；用 (\frac{1}{n-1}) 则是 4.75 和 0.25。主成分得分满足的关系仍然是“方差等于特征值”，因为两个定义是等比缩放。如果你拿手算结果和某软件输出对比，先确认软件的方差选项是总体方差还是样本方差。

5.3 软件底层用的是协方差矩阵还是相关矩阵

这是一个非常隐蔽的差异。有些软件 PCA 默认会把数据标准化，计算的是相关矩阵的特征分解；有些默认只做中心化，计算协方差矩阵的特征分解。同一份原始数据，两个选项得到的特征值、特征向量完全不同。

比如我上面这个例子，如果先标准化，每个变量的标准差都是 1，相关矩阵就是：

[ R = \begin{bmatrix} 1 & 0.9 \ 0.9 & 1 \end{bmatrix} ]

特征值变成 1.9 和 0.1，解释方差比例仍是 95% 和 5%，方向向量为 ((0.7071, 0.7071)) 和 ((0.7071, -0.7071))，与协方差矩阵分解结果相同。这是本例的特殊巧合，因为两个变量方差相等。一旦变量方差不等，协方差矩阵和相关矩阵的分解结果会有明显差异，所以看软件文档时，务必搞清楚默认选项。

5.4 载荷矩阵 V 的顺序或符号

还有一些软件会按特征值升序排列，把最小方差的成分放在第一列，如果你只取前 k 列，可能取到的是方差最小的方向。主成分分析的惯例是特征值降序排列，但程序内部不一定这么存，输出的顺序也可能和你手算习惯不一致。对比结果时，先看每个主成分解释的方差百分比，而不是硬对列号。

R 的princomp、Python 的PCA、手工用numpy.linalg.eig直接算，三者输出顺序都存在差异。最稳妥的做法是：先检查特征值序列，再看特征向量。

6. 从手算到数值方法：变量很多时 PCA 仍然只是线性代数

2×2 的例子适合理解原理，但真实场景中 p 可能有几十、几百甚至上千。这时候手写行列式展开是不现实的，必须理解软件在背后做的仍然是线性代数运算，只是换了更稳定的算法。

6.1 从协方差矩阵到特征值分解的数值算法

求矩阵特征值的方法有很多，最常见的是幂迭代法和 QR 算法。幂迭代法的思想很直观：从一个随机向量出发，不断左乘协方差矩阵并归一化，向量会逐渐收敛到最大特征值对应的特征向量。要想得到其他特征向量，可以从矩阵中减去已求出主成分的贡献，或者改用更完整的 QR 算法。

QR 算法的核心是反复做矩阵的 QR 分解，不断迭代，最终让矩阵逼近对角形式，对角线上的元素就是特征值。许多数值计算库的特征值函数底层采用的就是这类算法，在计算中不需要显式求解特征多项式，因为高次多项式求根在数值上非常不稳定。

6.2 SVD与特征值分解的关系

处理 PCA 常见的前沿操作，还有一条路线是奇异值分解。对中心化矩阵 (X_c) 做 SVD：

[ X_c = U \Sigma V^T ]

其中 (U) 是 n×p 左奇异向量矩阵，(\Sigma) 是奇异值对角矩阵，(V) 是 p×p 右奇异向量矩阵。奇异值与协方差矩阵特征值之间有关系：如果协方差矩阵是 (\frac{1}{n-1}X_c^T X_c)，那么它的特征值等于 (\frac{\sigma_j^2}{n-1})，其中 (\sigma_j) 是第 j 个奇异值。

主成分得分有两种等价算法：

[ Z = X_c V ]

或：

[ Z = U \Sigma ]

两个公式得到的数值在浮点误差范围内完全一致。SVD 的好处是数值稳定性更好，而且避免了先构造协方差矩阵 (S)。当变量数 p 很大时，(p \times p) 的协方差矩阵本身就可能很庞大；直接从 (X_c) 做 SVD，省去了协方差矩阵的构造步骤，因此很多现代 PCA 实现都走 SVD 路线。

从线性代数角度理解全部内容，其实是一幅很完整的画面：PCA 的全部信息都藏在矩阵的特征结构里。数据矩阵通过中心化、标准化变成形状规则的矩阵，通过矩阵乘法构成对称矩阵，通过特征值分解提取方向和重要程度，最后通过与载荷矩阵相乘得到投影之后的得分矩阵。无论算法如何变化，数学骨架始终不变。如果用人话总结这段经验：手算 PCA 的时候，不要把步骤割裂成“算特征值”“算特征向量”两个任务，而要把它们看作同一条链路上的产物——特征值决定方差大小，特征向量决定投影方向，得分矩阵才是最终用于分析的数据。我自己的习惯是，每做完一步都用迹、列和、正交性这些线性代数性质验算一遍，这些小检查比反复看教程管用得多。