课堂上听老师讲PCA总觉得思路特别清楚标准化、算协方差、求特征值特征向量。可一旦自己动手算尤其是要从一堆矩阵运算里把结果一步步推出来经常对不上答案。后来我重新把PCA当成线性代数题来做才意识到所谓计算步骤其实是一条完整的矩阵链条从数据矩阵出发得到中心化矩阵再到协方差矩阵最后解特征方程、拼出载荷矩阵、求得主成分得分。这篇文章就把这条链路上的每一个环节都拆开讲配一个能从头算到尾的数值例子帮你把PCA的计算步骤真正吃透。适合正在复习线性代数、准备数据分析考试、或者刚接触多元统计打算搞懂原理的读者。1. 把PCA翻译成一件事求协方差矩阵的特征基PCA的计算步骤从数学上看非常统一。它的目标不是“降维”这个动作本身而是找到一组新的正交坐标轴让数据在这些轴上的投影方差依次达到最大。如果把它拆成线性代数语言每个主成分其实就是数据协方差矩阵的一个特征向量而主成分的方差就是对应的特征值。1.1 行样本、列变量的约定先约定记号。设原始数据矩阵为 X维度是 n×pn 代表样本数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变量2112223335444556变量1的均值是 3变量2的均值是 4于是[ X_c X - \begin{bmatrix} 1 \ 1 \ 1 \ 1 \ 1 \end{bmatrix} \begin{bmatrix} 3 4 \end{bmatrix} ]结果就是每列都减去对应均值样本变量1变量21-2-22-1-1301410522中心化后的数据每一列的均值都是 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-14)得到[ 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 原始数据与均值中心化原始数据样本变量1变量2A12B23C35D44E56均值向量[ \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-10120)第二列 (-2-11020)通过。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 的时候不要把步骤割裂成“算特征值”“算特征向量”两个任务而要把它们看作同一条链路上的产物——特征值决定方差大小特征向量决定投影方向得分矩阵才是最终用于分析的数据。我自己的习惯是每做完一步都用迹、列和、正交性这些线性代数性质验算一遍这些小检查比反复看教程管用得多。