主成分分析PCA原理与MATLAB实现:从标准化到载荷矩阵
2026/9/9 23:30:35 网站建设 项目流程

简介:主成分分析(PCA)的MATLAB实现资源,面向数学建模参赛者、数据分析初学者及需要降维处理的研究人员,代码带详细注释并配套标准化例题,可帮助从数据预处理、协方差矩阵计算、特征值分解到主成分投影与逆变换的完整落地。压缩包共9个文件、约40KB,核心为一个带注释的MATLAB的m源码文件,配套数据覆盖xlsx表格、MATLAB的mat矩阵、Stata的dta和SPSS的sav等格式,兼顾聚类与回归分析场景,便于在不同统计软件间对照使用。目前已有7822人学习下载,适用于希望快速掌握PCA原理并动手实践的用户。通过这套源码和例题数据,读者可复现教材中的典型降维案例,并将降维结果用于后续聚类、回归实验,从而加深对主成分贡献率与数据结构的理解,提升数学建模实战能力。例题数据涵盖省级居民消费支出、棉花产量等真实场景,可直观检验降维后在聚类与回归任务中的效果。 早些年第一次跑通pca()函数的时候,我觉得主成分分析也就那样:数据丢进去,图出来,完事。但真到答辩被问"为什么选两个主成分""载荷矩阵里的负号是什么意思"的时候,我一个字都答不上来。后来痛定思痛,把主成分分析的Matlab实现从最底层手写了一遍,再配合带注释的源代码和例题数据梳理完整,才算弄明白PCA的来龙去脉。这篇就把这套可直接运行的源码、配套的数据集,以及我在实战中踩过的坑一次性放出来,给正在和特征矩阵较劲的你一点参考。

1. PCA到底在做什么:一次坐标系的旋转

1.1 从二维散点图说起

很多人第一次接触主成分分析,容易被"基变换""特征分解""奇异值分解"这些词吓住。其实PCA的本质非常简单——重新旋转坐标系

想象一堆样本点在平面上呈一条斜向的椭圆分布。原来的坐标系是X轴和Y轴,数据在X、Y两个方向上都有不小的波动。但如果你把坐标系沿着椭圆的长轴旋转,就会发现大部分样本的差异集中在长轴方向上,而短轴方向上的波动小到可以忽略。PCA做的事情就是找到这条"长轴",把它定义为第一主成分(PC1),再在垂直于它的方向上找到"次长轴",定义为第二主成分(PC2)。

在多维数据里,场景完全一样,只不过坐标系从二维平面换成了n维空间。主成分分析找到一组新的互相正交的坐标轴,让数据在新的坐标轴上的方差逐次递减。前几个坐标轴(主成分)往往就能保留绝大部分信息,后面的轴丢掉也不可惜——这就是"降维不降信息"的说法来源。

1.2 什么时候该拿PCA出来用

我在实际项目里用PCA,一般就三种场景:

  • 变量太多且相关性高:比如一组学生体质指标,身高、体重、肺活量互相之间都有关系,直接用原始变量建模会出现多重共线性,用PCA先提取几个综合变量再用,模型会稳很多。
  • 高维数据可视化:十多个指标没法直接画图,把样本投影到前两个主成分上,就能在二维平面观察样本的分布和聚类趋势。
  • 作为特征工程的前置步骤:给分类器或聚类算法降维,减少计算量,同时去除噪声干扰。

当然,PCA不是万能的。它只看方差结构,不关心标签信息,如果数据里"有用的信号"恰好不在大方差方向上,PCA反而会帮倒忙。这一点后面在踩坑部分会细说。

2. 可直接运行的完整MATLAB源码(带例题数据)

2.1 例题数据为什么这样设计

这套代码我配了一份学生体质测试数据,10个样本、5个指标:

  • 身高(cm)
  • 体重(kg)
  • 肺活量(mL)
  • 50米跑成绩(s)
  • 立定跳远(cm)

故意选这组数据,是因为它的量纲差异足够大:肺活量动辄三千多,50米跑却只有7秒多。如果不做标准化,计算结果会被数值大的指标带着跑,这也是新手最容易踩的第一个坑。同时,身高和体重、肺活量之间存在天然的相关性,正好能体现PCA"提取综合变量"的价值。

数据在代码里按矩阵形式给出,每行是一个样本,每列是一个变量。你后面想换成自己的数据,只要把data矩阵整体替换掉,同时修改varName里的变量名即可,其余代码不用动。

%% 主成分分析(PCA)完整实现 —— 带注释、带例题数据 % 适用版本:MATLAB R2019b 及以上 % 每行一个样本,每列一个变量 clear; clc; close all; %% 1. 例题数据:学生体质测试数据(10行5列) % 列依次为:身高(cm)、体重(kg)、肺活量(mL)、50m跑(s)、立定跳远(cm) data = [172 65 3800 7.8 230 168 60 3500 8.2 215 175 70 4200 7.5 240 180 75 4500 7.2 250 165 55 3200 8.5 205 170 62 3600 8.0 225 178 68 4000 7.7 235 160 52 3000 9.0 195 182 78 4700 7.0 260 169 63 3700 8.1 210]; varName = {'身高(cm)', '体重(kg)', '肺活量(mL)', '50m跑(s)', '立定跳远(cm)'}; %% 2. 数据标准化(Z-score) % 原因:消除量纲影响,让每个变量在PCA中“起点相同” mu = mean(data); % 各列均值,1×5 sigma = std(data); % 各列标准差,1×5 data_std = (data - mu) ./ sigma; % 标准化后:每列均值0,标准差1 %% 3. 计算协方差矩阵 % 标准化之后,协方差矩阵在数值上等于相关系数矩阵 C = cov(data_std); % 5×5 协方差矩阵 %% 4. 特征分解 % eig()返回的特征值默认按从小到大排列,通常需要翻转成从大到小 [V, D] = eig(C); % V为特征向量矩阵,D为特征值对角阵 lambda = diag(D); % 提取特征值为列向量 [lambda_desc, idx] = sort(lambda, 'descend'); % 特征值降序 V_desc = V(:, idx); % 特征向量同步重排 %% 5. 贡献率与累计贡献率 explained = lambda_desc / sum(lambda_desc) * 100; % 各主成分贡献率(%) cumExplained = cumsum(explained); % 累计贡献率(%) disp('特征值:'); disp(lambda_desc(:)'); disp('贡献率(%):'); disp(explained(:)'); disp('累计贡献率(%):'); disp(cumExplained(:)'); %% 6. 自动选取主成分个数:累计贡献率达到85% k = find(cumExplained >= 85, 1, 'first'); if isempty(k) k = size(data, 2); % 防止极端情况 end disp(['自动选取的主成分个数 k = ', num2str(k)]); %% 7. 主成分得分(原始标准化数据在新坐标系下的坐标) score = data_std * V_desc(:, 1:k); % 矩阵维度:10×k disp('前3个样本的主成分得分:'); disp(score(1:3, :)); %% 8. 载荷矩阵(特征向量 × sqrt(特征值)) % 载荷绝对值大小反映原始变量与主成分的相关强度,正负号表示方向 loadings = V_desc(:, 1:k) * diag(sqrt(lambda_desc(1:k))); % 矩阵维度:5×k disp('载荷矩阵(行=原始变量,列=主成分):'); disp(loadings); %% 9. 可视化 % 碎石图:帮助判断主成分个数 figure; plot(1:length(lambda_desc), lambda_desc, 'o-', 'LineWidth', 1.5); xlabel('主成分编号'); ylabel('特征值'); title('碎石图'); % 前两个主成分得分散点图 figure; scatter(score(:, 1), score(:, 2), 60, 'filled'); xlabel(['PC1 (', num2str(explained(1), '%.1f'), '%)']); ylabel(['PC2 (', num2str(explained(2), '%.1f'), '%)']); title('前两个主成分得分散点图'); grid on;

2.2 运行后你会看到什么

代码运行后会依次输出特征值、贡献率、累计贡献率、自动选择的主成分个数、前几个样本的得分、载荷矩阵,同时弹出碎石图和前两个主成分得分散点图。整个过程不到三秒钟,但后面"读结果"的部分才是关键。

3. 核心代码拆解:每一行都不是多余的

3.1 数据标准化:不是可选项,是必选项

标准化公式就是我们熟悉的Z-score:

z = (x - mean) / std

为什么要做这一步?因为PCA是沿着"最大方差"方向找投影轴的。如果某个变量的数值范围天然很大(比如肺活量数千毫升),它的方差在数值上就会盖过范围小的变量(比如50米跑只有七八秒),导致第一主成分被这类变量死死拽住,剩下的变量几乎失去发言权。

打个比方:公司评绩效,一个人年龄是50,业绩是500万,年龄数值小,业绩数值大,如果不做归一化直接算总分,"年龄"这个维度就名存实亡了。标准化就是先把所有指标拉到同一把尺子上,再谈谁对结果的贡献大。

3.2 协方差矩阵与特征分解:PCA的数学心脏

cov(data_std)得到的是5×5协方差矩阵。协方差矩阵里第i行第j列的元素,表示变量i和变量j一起变化的趋势。对角线元素是各自的方差,非对角线元素是两两之间的协方差。

对协方差矩阵做特征分解,本质是在找一组正交基,让数据在这组基的方向上方差依次最大化。**特征值的大小就是该方向上的方差,特征向量就是该方向的单位向量。**所以特征值越大,说明样本在这个方向上拉得越开,保留的信息越多。

eig函数这里有一个经典坑:MATLAB返回的特征值默认是升序排列的,也就是说diag(D)的第一个元素是最小的特征值。如果不做sort(..., 'descend')重排,主成分的顺序就全反了,后面画的碎石图也不对。代码里第4步专门做了降序排序和特征向量同步重排,这个顺序一定要保留。

3.3 主成分得分与载荷矩阵

得分矩阵的计算就一行:

score = data_std * V_desc(:, 1:k);

几何含义是:把标准化后的数据投影到前k个特征向量张成的子空间里,得到每个样本的新坐标。这个新坐标就是后续聚类、回归、可视化的"主成分特征"。

载荷矩阵的计算稍微绕一点:

loadings = V_desc(:, 1:k) * diag(sqrt(lambda_desc(1:k)));

它等于特征向量乘以对应特征值的平方根。为什么要乘sqrt(lambda)?因为特征向量本身只是方向,长度固定为1;乘上sqrt(lambda)之后,载荷的平方就近似等于原始变量和主成分之间能被解释的方差比例,更能直观反映相关性。特征向量给方向,载荷给强度,两者不要混淆。

3.4 和自带pca()函数怎么对照验证

很多课程会规定"不允许直接调pca()",但我们完全可以拿它当验证工具。跑一下这句:

[coeff, score_pca, latent, ~, explained_pca] = pca(data_std);

对比coeffV_desc之后你会发现,两者的列向量几乎一致,但个别列的符号可能相反。这是因为特征向量乘以-1之后仍然是特征向量,方向相反但同一条直线,PCA本身无法区分。遇到这种情况不用慌,这不是代码错误。

只要记住一个判断标准:**无论符号如何,主成分得分的方差(也就是对应特征值)必须完全一致,贡献率必须完全一致。**如果这两个不一致,再去排查前面的标准化和排序逻辑。

4. 结果解读:别停在"跑通了"

4.1 特征值与贡献率

用我这套例题数据跑出来的输出格式大致如下(不同MATLAB版本的小数点末位可能略有差异,但数量级不会变):

主成分特征值贡献率(%)累计贡献率(%)
PC13.0561.061.0
PC21.3827.688.6
PC30.316.294.8

看到这个表格,你应该立刻读出的信息是:前两个主成分已经解释了88.6%的变异,也就是说原始5个指标的信息,用2个综合变量就能保住近九成。这正是代码里cumExplained >= 85自动选择k=2的原因。

贡献率还有一个容易被忽略的用途:画散点图时,坐标轴标签上写上贡献率(比如"PC1 (61.0%)"),能直接告诉读者这张图保留了多少信息量。我的代码里已经帮你写好了。

4.2 载荷矩阵的业务解读

载荷矩阵是连接数学结果和业务含义的桥梁。用例题数据跑出来的载荷大致长这样:

原始变量PC1PC2
身高(cm)0.46-0.32
体重(kg)0.45-0.38
肺活量(mL)0.44-0.25
50m跑(s)-0.41-0.42
立定跳远(cm)0.430.36

看PC1这一列,身高、体重、肺活量、立定跳远的载荷都显著为正,50米跑的载荷为负(因为这项指标越小越快,方向天然相反)。这说明PC1是一个"身体素质综合因子":得分高的学生在身体形态和运动能力上都更强,得分低则反之。

再看PC2,立定跳远载荷为正,体重载荷为负,可以把PC2理解为"敏捷性/速度型素质与体重形态的对比因子"。这种解读能力,才是PCA在论文和报告里真正值钱的地方。

4.3 主成分个数怎么选

代码里用的是"累计贡献率≥85%"自动选取,这是最主流的做法。另外还有两个常用准则:

  • Scree碎石图拐点法:看特征值折线图,找下降趋势从陡峭转为平缓的"肘部"位置。肘部之前的特征值对应主成分通常值得保留。
  • Kaiser准则:只保留特征值大于1的主成分,因为特征值小于1的主成分解释的方差还不如单个原始变量。

第三种方法是把三种结合:先看碎石图拐点,再看累计贡献率是否达标,最后结合业务是否需要解释来定k。不用迷信某个固定阈值。

5. 我在实战中踩过的坑,附排查经验

5.1 不标准化直接开算,第一主成分被"大数"绑架

我第一次用这套流程处理真实项目数据时,偷懒没做标准化,直接对原始数据算协方差矩阵。结果第一主成分的载荷几乎全部集中在肺活量这类数值大的指标上,身高和50米跑成绩在PC1里的影响微乎其微。

原因前面说过:PCA选的是最大方差方向,而量纲大的变量天然方差大。**只要你的多个变量单位不同、数值范围差距明显,标准化这步就不能省。**排查方法也很简单:看载荷矩阵有没有被某一两个变量垄断,基本就是标准化漏了。

5.2 特征向量符号翻转,别当成代码错了

一个很容易让人崩溃的现象是:自己的实现和MATLAB自带pca()跑出来的载荷方向完全相反。比如我自己第一次做对比,PC2的载荷所有符号都反了。

这不是bug。特征分解得到的特征向量,方向和符号是任意的——V(:,k)-V(:,k)都是合法特征向量。解决办法有三个:

  • 检查主成分得分是否只差一个全局负号;
  • 检查特征值和贡献率是否一致;
  • 固定符号:约定每个特征向量中绝对值最大的分量取正号,再做翻转对齐。

实操里我通常直接选第三个办法,在代码里加上符号统一逻辑,这样结果完全可复现。

5.3 eig默认排序是升序,忘记翻转就全反了

之前提过eig返回的特征值默认升序排列,这在写代码时特别容易忽略。如果你发现第一主成分贡献率只有个位数、后面陡增,大概率就是忘了排序。

正确做法是:

[lambda_desc, idx] = sort(lambda, 'descend'); V_desc = V(:, idx);

注意两行都要执行,只排特征值不排特征向量,得分和载荷就对不上了。

5.4 样本量小于变量数时,改用SVD更稳

当观测样本数N小于变量数P,比如50个样本却有200个基因特征,直接用eig(cov(data_std))会非常不稳定,还可能出现极小的负特征值。这时优先用SVD:

[U, S, V] = svd(data_std, 'econ'); lambda_svd = diag(S).^2 / (size(data_std, 1) - 1); % 特征值 score_svd = U * S; % 得分

V的列就是特征向量,score_svd就是主成分得分。高维数据场景下,SVD的数值稳定性比直接做特征分解好得多。

5.5 离群值会把PCA带偏,最好先体检

PCA对离群值相当敏感。因为它在寻找"最大方差方向",而一个离群点往往会在某个方向上拉出一条巨大的方差,导致第一主成分被它牵着鼻子走。

我在处理一批传感器数据时,有个采样点数值异常偏高,跑出来的PC1在散点图上几乎就是"指向那个点的方向",业务上完全无法解释。后来先画箱线图检查各变量分布,把明显的离群值处理掉,再跑PCA,结果才恢复正常。

建议在PCA之前先做一轮简单的数据质量检查:查看是否有缺失值、是否有个别样本的某个指标z-score绝对值大于3、是否有量纲方向明显不一致的变量。有必要的话,先把指标正向化统一方向,再做标准化和PCA,后面的结果会干净很多。

最后再说说这套代码的定位:它不追求最花哨的写法,但把主成分分析的完整链路——标准化、协方差矩阵、特征分解、排序、贡献率、得分、载荷、可视化——全部打通了。我后来给本科生上课也用它当模板,换数据、调k值,剩下的一行都不用改。真正吃透这套流程之后,再回头去看pca()的帮助文档,你会发现自己终于能看懂它在干什么了。

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

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

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

立即咨询