☰
主成分分析MATLAB手写实现:协方差矩阵、特征分解与降维实战
2026/10/13 18:42:49 网站建设 项目流程

简介:一份面向数学建模与机器学习学习者的主成分分析(PCA)Matlab源代码资料包,内含带详细注释的PCA实现代码及多组例题数据,可完整演示从数据预处理、协方差矩阵计算、特征值分解到投影与逆变换的全流程。代码结构清晰,逐步注释关键计算环节,适合初学者对照入门,也便于在降维、聚类、回归等场景中直接复用。压缩包共9个文件,包含1个m源文件、4个xlsx表格数据、2个mat矩阵数据,以及dta和sav格式数据,覆盖统计软件常用数据接口,整体仅40KB,轻量易下载。数据内容涉及棉花产量论文作业、1999年全国31个省份城镇居民消费性支出、课堂奶粉数据等真实案例,可用于验证PCA在回归与聚类分析中的实际效果。已有7825人学习下载,是快速掌握PCA原理与Matlab实现的实用参考资料,适合课程设计、数模竞赛与科研入门。

1. 主成分分析为什么值得在 MATLAB 里自己写一遍

标题里的「主成分分析 matlab 源代码(带注释,带例题数据)」听起来像是一个现成的代码包,其实更像一条从原理走向落地的主线。在 MATLAB 里只调pca()一行函数,换来的只是一个黑匣子:特征向量符号为什么时正时负、数据没标准化会有什么后果、解释方差比应该怎么支撑决策,这些细节封装函数不会替你回答。这份笔记写一套带注释的源代码,配一份 10 个样本、5 个变量的例题数据,从手写协方差矩阵和特征分解开始,完整走一遍主成分分析流程,最后与内置函数互检、把常见坑交代清楚。适合两类人:刚接触 PCA 想一次把原理落到代码的初学者,以及会调用函数但想搞清边界与隐患的实战派。

2. PCA 原理与数据准备:从几何直觉到矩阵布局

2.1 主成分的几何直觉:旋转坐标轴,让投影方差最大

想象一堆点散布在一个二维平面上,沿着某个方向把这些点投影到一条直线上,投影点散得最开的方向,就是第一主成分。垂直于第一主成分再投影一次,剩下的点散得最开的方向就是第二主成分。三维以上依此类推。所以在直觉层面,PCA 就是在做一次正交旋转,把原始坐标系换成一个方差从大到小排列的新坐标系。每一根新坐标轴都是原始变量的线性组合,组合系数叫载荷。

代数上,PCA 最终归结为对协方差矩阵做特征分解。设 X 是已经中心化的样本矩阵,每行一个样本、每列一个变量,那么协方差矩阵 C = (X'X)/(n-1)。C 的特征向量就是主成分方向,特征值是该方向上的方差。特征值越大,这个方向承载的信息越多。在 MATLAB 里做特征分解有eig()和svd()两条路线,数值稳定性不同,后面都会讲到。

这里最需要记住的一句话是:PCA 输出的是新特征,不是从原始变量里挑出几个。决定要用哪几个主成分,看的是解释方差;决定每个主成分由哪些原始变量主导,看的是载荷。这两样东西在第 4 章专门讲。

2.2 例题数据布局:一行样本、一列变量,先检查尺寸

数据矩阵的摆放方向是第一道容易被忽略的坑。MATLAB 的pca()期望输入是“行样本、列变量”,但很多数学推导习惯把样本写成列向量。我有一次把数据按“列是样本”传进去,输出的主成分数直接变成样本数,整个分析方向都偏了,排查半天才发现是转置问题。

我的建议是给自己定一条死的约定:每行是样本,每列是变量,所有下游步骤都围绕这个布局展开。例题数据我用 10 个产品批次、5 个质量指标来构造,变量名暂定为 x1 到 x5,分别代表硬度、工艺温度、压力、浓度和反应时间。数值是随意构造的,只是为了把代码跑通、把读数逻辑讲明白。

% 示例数据矩阵 X % 每行是一个样本(产品批次),每列是一个变量(质量指标) % 变量 x1~x5 分别代表:硬度、温度、压力、浓度、反应时间 X = [82 380 1.20 0.13 4.8 86 410 1.05 0.15 5.1 78 350 1.35 0.11 4.5 85 400 1.18 0.14 5.3 80 395 0.98 0.16 4.9 90 420 1.22 0.12 5.6 77 360 1.28 0.17 4.6 88 405 1.02 0.13 5.2 84 390 1.10 0.15 5.0 81 370 1.32 0.12 4.7];

数据录进去之后,第一件事是确认尺寸,不要靠数数。

[n, p] = size(X); % 期望 n=10, p=5 fprintf('样本数 n=%d,变量数 p=%d\n', n, p);

这行代码花不到一秒钟,却能挡住一半的翻车事故。接下来观察一下数值范围:温度是 350~420,压力是 0.98~1.35,浓度只有 0.11~0.17,量纲差距明显。如果直接拿原始矩阵跑 PCA,第一主成分必然被温度和硬度这种大数变量带跑偏,这就是下一小节要处理的问题。

2.3 标准化:PCA 的高压线,什么情况可以跳过

标准化的标准写法只有一行:

X_std = zscore(X); % 每列减去均值、除以标准差 zscore 会逐列计算,把每个变量转成均值为 0、标准差为 1 的尺度。

用zscore的理由很实在:协方差矩阵的数值由变量方差主导,方差大的变量在矩阵里“嗓门”更大。温度列标准差有四十几,浓度列标准差可能只有零点零几,不缩放的话,浓度对主成分的贡献几乎可以忽略。

什么情况可以不标准化?只有当所有变量本身是同一个物理量、同一个单位时,跳过才合理。比如 5 个温度传感器测同一条反应釜的温度曲线,量纲相同,保留原始方差是说得通的。除此之外,默认 zscore,不要犹豫。

还要区分一个常被混为一谈的概念:中心化和标准化。中心化只减均值,标准化是减均值后再除以标准差。协方差矩阵的计算本身要求数据已中心化,但如果不做方差缩放,高方差变量依然占主导。zscore 一步把两件事全做了,之后协方差矩阵的对角线就是各变量标准化后的方差,数值上全部为 1,PCA 才有公平的起点。

2.4 降维的本质:丢掉的是方差小的方向,不一定是噪声

降维到底丢掉了什么?很多人理所当然地认为后几个主成分就是噪声。这个理解在信噪比高的时候勉强成立,但 PCA 本身并不区分信号和噪声,它只按方差排序。如果原始数据本身就含有一部分可区分的结构,恰好落在方差小的方向上,强行丢掉它就会损失信息。

PCA 真正擅长压缩的是冗余。当原始变量高度相关时,协方差矩阵的有效秩比变量数低,少数主成分就能覆盖绝大部分方差,丢掉后面的维度相当于丢掉重复信息,风险很小。实际项目里,判断依据通常是两条:累计解释方差到达 85%~95%,或者碎石图上出现明显的“拐弯”。这两种做法我在第 4 章再展开,现在的关键结论是——主成分个数不能写死,要交给数据和业务阈值一起决定。

3. 手写 PCA 的 MATLAB 源代码:协方差矩阵、特征分解与内置互检

3.1 最小实现:eig() 版的九行核心代码

有了标准化的数据,PCA 主体只剩下协方差矩阵和特征分解两件事。

% 1. 协方差矩阵:注意除以 n-1 而不是 n,这是无偏估计 C = (X_std' * X_std) / (n - 1); % 2. 特征分解:V 的每一列是一个特征向量,D 的对角元素是特征值 [V, D] = eig(C); % 3. eig() 默认把特征值从小到大排,必须手动改成降序 [lambda, idx] = sort(diag(D), 'descend'); V = V(:, idx); % 特征向量按同样的顺序重排 lambda = lambda'; % 转成行向量,方便后续输出 % 4. 主成分得分:标准化后的数据在新坐标轴上的投影 Score = X_std * V; % 5. 解释方差百分比和累计贡献率 explained = lambda / sum(lambda) * 100; cum_explained = cumsum(explained);

几个关键点需要说明。协方差矩阵除以 n-1 是约定俗成的无偏估计,用cov(X_std)会更省事,但这里手写是为了让你看清矩阵长什么样。eig()返回的 V 每一列是一个特征向量,对应 D 对角线上相同位置的特征值;但编译器不会帮你排顺序,必须用sort(..., 'descend')手动把特征值从大到小排列,特征向量的列也要同步重排,否则后续主成分顺序全反。

主成分得分Score = X_std * V是 n×p 的矩阵,第一列就是所有样本在第一主成分上的投影坐标。这个积一定要写在纸面上理解:V 是 p×p 的载荷矩阵,X_std 是 n×p 的数据矩阵,相乘得到 n×p 的新坐标。很多新手在这里少转置一次,维度直接对不上。

3.2 内置函数对照:pca() 的输出怎么对齐

新版 MATLAB 推荐用pca(),旧版里常见的princomp()已经进入弃用流程。pca()的标准调用方式如下:

[coeff, score, latent, ~, explained] = pca(X_std); % coeff: p×p 载荷矩阵,每列是一个主成分 % score: n×p 主成分得分,和手写的 Score 对应 % latent: p×1 特征值向量 % explained: p×1 解释方差百分比

拿这个结果和手写版本互检,大多数情况下你会发现两个版本的解释方差几乎一致,但载荷矩阵的某些列可能出现整体正负号相反。这不是代码错误。特征分解得到的特征向量方向本身就不唯一,v 和 -v 都是特征向量,任何数值库都可以选其中任意一个返回。

所以互检的正确姿势是比较解释方差百分比,或者比较投影后点云的分布形状,而不是逐元素比较向量符号。如果非要让符号对齐,可以自己做一层归一化:把每列绝对值最大的元素强制设成正数。这个技巧在 5.1 小节给出具体代码。

3.3 更稳的 SVD 版本:什么时候该选它

eig(C)的实现是先算协方差矩阵再分解。当 p 非常小、n 非常大时,这是最快的路径;但当变量高度共线,协方差矩阵接近奇异时,eig()会有数值精度问题。这时候用 SVD 直接分解标准化后的数据矩阵,稳定性更好。

% SVD 版本:直接对 X_std 做奇异值分解 [U, S, V] = svd(X_std); % X_std = U * S * V',所以主成分得分 = U * S,载荷 = V score_svd = U * S; % n×p 矩阵 % 奇异值平方再除以 n-1 才是特征值 lambda_svd = (diag(S).^2) / (n - 1); explained_svd = lambda_svd / sum(lambda_svd) * 100;

这里最容易踩的坑是把奇异值当特征值直接用。S 的对角元素是奇异值,PCA 的特征值是奇异值的平方除以 n-1。如果直接用奇异值去算解释方差,比例会是错的。SVD 版本和 eig 版本的结果在大多数数据集上数值非常接近,只是符号可能再次不一致,这是正常的,不需要恐慌。样本量不大、没有严重共线问题的场景,用哪个都一样;有强迫症或者数据病态时选 SVD。

3.4 完整脚本:一份可以直接运行的 .m 文件

把上面所有片段合并成一个脚本,方便直接复制粘贴运行。数据读入、标准化、特征分解、结果打印、内置对照,五个部分集中在同一个文件里。

% ============================================= % 主成分分析完整脚本(带注释、带例题数据) % 数据布局:X 为 n 行 × p 列,每行一个样本 % ============================================= clear; clc; % ---- 1. 数据 ---- X = [82 380 1.20 0.13 4.8 86 410 1.05 0.15 5.1 78 350 1.35 0.11 4.5 85 400 1.18 0.14 5.3 80 395 0.98 0.16 4.9 90 420 1.22 0.12 5.6 77 360 1.28 0.17 4.6 88 405 1.02 0.13 5.2 84 390 1.10 0.15 5.0 81 370 1.32 0.12 4.7]; [n, p] = size(X); fprintf('样本数 n=%d,变量数 p=%d\n', n, p); % ---- 2. 标准化 ---- X_std = zscore(X); % ---- 3. 协方差矩阵与特征分解 ---- C = (X_std' * X_std) / (n - 1); [V, D] = eig(C); [lambda, idx] = sort(diag(D), 'descend'); V = V(:, idx); lambda = lambda'; % ---- 4. 主成分得分 ---- Score = X_std * V; % ---- 5. 解释方差 ---- explained = lambda / sum(lambda) * 100; cum_explained = cumsum(explained); % ---- 6. 打印结果 ---- fprintf('\n手写 PCA 结果:\n'); for i = 1:p fprintf('PC%d | 特征值 %.4f | 解释方差 %.2f%% | 累计 %.2f%%\n', ... i, lambda(i), explained(i), cum_explained(i)); end % ---- 7. 与内置 pca() 对照 ---- [~, ~, latent, ~, explained_pca] = pca(X_std); fprintf('\n内置 pca() 对照:\n'); for i = 1:p fprintf('PC%d | 特征值 %.4f | 解释方差 %.2f%%\n', ... i, latent(i), explained_pca(i)); end

脚本第 6 步的fprintf循环是调试时最常用的输出方式。特征值、解释方差、累计贡献率三列并排打印,一眼就能看出该取几个主成分。对照部分如果发现特征值完全一致但载荷符号不同,直接忽略符号差异即可。

4. 结果解读与可视化:解释方差、载荷矩阵与三张图

4.1 解释方差表:主成分个数不靠猜,靠阈值

跑完上面的脚本,控制台会打出一张类似下面的表:

主成分特征值解释方差 %累计贡献率 %
PC12.8757.457.4
PC21.3226.483.8
PC30.6112.296.0
PC40.122.498.4
PC50.081.6100.0

以上数据是示意,实际数值取决于你的 MATLAB 版本和数据处理细节,但结构是这样。看这张表时,真正要决策的是“取几个主成分”。最常见的做法是取累计贡献率达到 85%~95% 的个数,上表里 PC1 加 PC2 到 83.8%,加 PC3 到 96%,所以取前 3 个主成分就能覆盖绝大多数信息。

另一个办法是看特征值之间的落差。如果某个主成分的特征值突然比前一个小很多,说明新增方向提供的信息量骤降,这个位置就是天然的分界。两种方法可以结合用,先用累计贡献率定下限,再用特征值落差确认上限。

4.2 载荷矩阵:每个主成分的“配方表”

载荷矩阵 V 的每一列是一个主成分,每一行对应一个原始变量。某个变量在某主成分上的载荷绝对值越大,说明它对这个主成分的贡献越大。

% 打印第一主成分的载荷 fprintf('PC1 载荷向量:\n'); for i = 1:p fprintf('x%d = %.4f\n', i, V(i, 1)); end

假设输出是 x1 = 0.52、x2 = -0.41、x3 = 0.63、x4 = 0.18、x5 = -0.31,那么第一主成分主要由 x3 和 x1 主导,x4 几乎可以忽略。这只是一种可能的结果;实际数据不同,载荷分布也不同。

读载荷时要注意符号:正负只代表方向,不代表好坏。某个变量载荷为负,表示它与该主成分的运动方向相反,绝对值大小才是权重。想把第一主成分的“综合指标”含义讲给业务方听,就看绝对值排名前两到三个变量。

4.3 三张图:碎石图、载荷柱状图、得分散点图

三张图能从不同角度验证 PCA 结果,初学阶段建议全部画一遍。

% 图1:碎石图,横轴为主成分序号,纵轴为特征值 figure; plot(1:p, lambda, 'o-', 'LineWidth', 1.5); xlabel('主成分序号'); ylabel('特征值'); title('碎石图'); % 图2:第一主成分载荷柱状图 figure; bar(1:p, V(:,1)); xlabel('原始变量'); ylabel('载荷'); title('第一主成分各变量载荷'); set(gca, 'XTick', 1:p); set(gca, 'XTickLabel', {'x1','x2','x3','x4','x5'}); % 图3:样本在前两个主成分上的得分散点图 figure; scatter(Score(:,1), Score(:,2), 'filled'); xlabel('PC1'); ylabel('PC2'); title('样本主成分得分'); % 给每个点加上样本编号 for i = 1:n text(Score(i,1), Score(i,2), sprintf(' #%d', i)); end

碎石图的“阅读重点”是找拐点。前几个特征值陡降,后面进入平缓区,拐弯处的序号就是建议的主成分个数。载荷柱状图能直观看出哪些变量主导某个主成分。得分散点图把所有样本投影到前两个主成分平面上,是最常见的探索工具:点在图上聚集几团,就提示数据中可能存在几类结构。

如果在上面的得分图上看到某个点单独飞出,先别急着删数据。检查它的原始行取值,确认不是录入错误。PCA 对离群点很敏感,一个极端样本就可能把主成分方向整体拉过去,所以做 PCA 前多画几次箱线图是值得的。

5. PCA 避坑与常见问题排查:四个翻车现场的根因与解法

5.1 现象:两次跑 PCA,特征向量的符号对不上

这是最多人问的问题。昨天跑出来 coeff 第一列是正的,今天跑变成负的,代码没改,数据没换,结果却“变了”。原因在于特征分解的结果本来就不唯一。对于特征值 λ,如果 v 是特征向量,那么 -v 也是特征向量,数值库可以任意选择返回哪个方向。eig()、svd()、pca()对符号的选择逻辑不同,甚至同一函数在不同版本下都可能给出不同符号。

解决方法是做一次符号对齐。把每列中绝对值最大的元素调整为正值,这样不管底层库怎么选方向,最终结果都一致。

% 符号对齐:将每列绝对值最大的元素统一改为正号 for j = 1:p [~, max_idx] = max(abs(V(:, j))); if V(max_idx, j) < 0 V(:, j) = -V(:, j); end end

需要强调的是,这只影响表达方式,不影响解释方差和投影点云的相对位置。散点图整体翻转后,样本之间的远近关系不变,结论也不变。

5.2 现象:没做标准化,第一主成分被量纲劫持

症状很典型:主成分分析结果打印出来,PC1 解释方差超过 90%,载荷向量里某个变量的系数接近 1,其他变量几乎都接近 0。这说明第一主成分几乎等于那个变量本身,PCA 白做了。

原因就是量纲差异。温度列标准差几十,浓度列标准差零点零几,协方差矩阵里温度的对角线元素远大于浓度,特征分解自然优先沿着温度方向展开。PCA 变成了“找最大方差变量”,而不是“找变量组合”。

解决方法是默认对数据做zscore()。做完之后每个变量标准差都是 1,协方差矩阵对角线全是 1,变量之间的相关性才真正浮出水面。如果因为业务原因做不了标准化,至少要把每个变量的方差打印出来看一眼,确认没有某一个变量独占鳌头。

5.3 现象:样本数少于变量数,协方差矩阵奇异

当 n 小于 p 时,协方差矩阵的秩最多是 n-1,必然奇异,eig()可能计算出负的特征值或直接警告矩阵接近奇异。这个问题在基因表达谱、光谱分析、图像像素这类“高维小样本”数据里特别常见。

解决思路有三条:第一,在数学上,PCA 仍然可以做,但主成分个数最多只能有 n-1 个;第二,换 SVD 路径,直接分解 X_std,绕开构造协方差矩阵这一步,数值上更稳;第三,先做特征筛选,把明显无关的变量删掉,让 p 小于 n。第三条是最务实的工程做法,毕竟用 PCA 处理严重病态的数据,结果可信度也不高。

5.4 现象:得分图和原始数据对不上,坐标范围很别扭

有读者把Score的散点图画出来之后,发现坐标值跟原始数据的取值范围完全不一样,比如原始温度是三百多,得分图上却只有正负几,怀疑自己算错了。实际上这是预期行为。Score = X_std * V是用标准化后的数据去投影,原始数据的量纲和中心位置在标准化时已经被移除。

拿到得分图之前,心里要默默过一遍流程:原始数据 X 先减去均值、除以标准差,得到 X_std;然后 X_std 在载荷方向 V 上投影,得到主成分得分。得分是标准空间里的坐标,不属于原始量纲。如果想验证计算正确,可以用X_std(1,:) * V(:,1)和Score(1,1)对比,两个值应当完全相等。这一步手动验算能让人安心。

6. PCA 的进阶用法与验证习惯:从降维到异常检测

6.1 定主成分个数的两个可靠办法

累计解释方差阈值和碎石图拐点是最常用的两个判据。阈值法适合向业务方汇报,95% 比 90% 更保守;碎石图拐点适合探索阶段,快速决定主成分数量。两者不一致时,优先相信累计解释方差,因为拐点的位置受数据噪声影响较大。

6.2 把主成分得分当成特征喂给下游模型

降维后的Score矩阵可以直接作为分类器或回归模型的输入。这样做的好处是去掉冗余变量、降低过拟合风险、提升训练速度。但要注意,PCA 是在标准化的数据上计算的,所以下游模型的输入也是得分矩阵,而不是原始数据。保存模型时,标准化的均值、标准差和载荷矩阵都要一起保存,否则新数据无法完成同样的变换。

6.3 用重构误差做异常检测

PCA 还有一个容易被忽略的用法:把降维后的得分投影回原始空间,再和标准化数据比较,重构误差大的样本往往是离群点。

% 取前 k 个主成分做重构,k 根据累计贡献率选定 k = 3; X_recon = Score(:,1:k) * V(:,1:k)'; err = sum((X_std - X_recon).^2, 2); % 每个样本的重构误差

误差明显大于其他样本的那些点,就是 PCA 视角下的异常。这个方法不需要打标签,适合做数据质量筛查的第一步。我自己这几年用下来,只要是和“先跑通、再上生产”沾边的活,都保留同一个习惯:PCA 结果出来先打印解释方差和载荷矩阵,再画得分图,三样东西看完才决定下一步动作。这个顺序帮我避开过不少只看图说话的弯路。希望帮到你。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询