企业官网建设流程全解析

简介：面向需要掌握PCA原理与MATLAB实现的开发者，这份源代码覆盖了主成分分析的完整流程，包括数据标准化、协方差矩阵求解、特征值分解、按方差占比选取主成分，以及将数据投影到新坐标系并重建。资源共3个文件，含2个m脚本和1个xls示例数据文件，m脚本按步骤实现算法，xls用于验证降维结果，压缩包仅8KB，轻量易用。目前已有1137人学习浏览，实用性得到初步验证。通过学习这些代码，可理解PCA降维的完整逻辑，掌握zscore、cov、eig或svd等关键函数的实际用法，并可直接套用到特征提取、数据可视化、降噪等日常数据处理任务中，适合科研与工程场景快速参考。

1. 为什么要自己在 MATLAB 里写主成分分析源代码

主成分分析（PCA）是入门机器学习时几乎绕不开的降维手段。很多工程师点开 MATLAB 直接输入pca，一秒钟就能得到降维结果，但等到要解释载荷矩阵为什么是这个符号、主成分个数到底按什么标准取、换一批数据结果怎么就不稳定时，黑盒子就变成黑锅。自己维护一份主成分分析源代码的价值在于，把“数据中心化、协方差、特征值、得分”这条链路完整摆在面前，遇到异常时可以逐行排查，还能根据需要把函数复用到测试集投影、重建误差计算、批量交叉验证里。这篇内容适合正在用 MATLAB 做主成分分析、但又不想止步于调包的人。代码会从数学推导开始写，最后落成一个可直接调用的函数，附带参数边界和常见坑。

2. 主成分分析的数学骨架：从协方差矩阵到 SVD

2.1 标准化：让所有变量在源代码里站在同一条起跑线

PCA 的核心目标是寻找方差最大的方向。但这里的“方差”对量纲极其敏感。假设数据里一个变量是“长度”，单位是米；另一个变量是“时间”，单位是秒。长度的数值范围可能是几百万，时间可能只有零点几。如果不做任何处理，协方差矩阵里长度带来的方差会压过时间，第一主成分几乎由长度决定，结果完全没有解释性。

因此，自写主成分分析源代码时，我一般把标准化作为默认前置步骤。标准化的含义有两个层次：一是中心化，让每个变量的均值为 0；二是缩放，让每个变量的方差为 1。也可以只做中心化。区别在于，只去均值时 PCA 寻找的是“围绕原点的最大方差方向”，同时除以标准差时每个变量对距离的贡献权重相同。

下面是最基础的标准化代码：

X = meas; % 假设 meas 是 n×p 原始数据矩阵 mu = mean(X, 1); % 沿第 1 维求均值，返回 1×p 向量 Xc = X - mu; % 去均值中心化 sigma = std(Xc, 0, 1); % 计算每个变量的样本标准差 sigma(sigma < eps) = 1; % 防止零方差变量导致除零 Xs = Xc ./ sigma; % 除以标准差，完成标准化

这里mean(X, 1)的第二个参数是关键。MATLAB 默认对列求均值，但写清楚维度能避免数据行列方向搞反。std(Xc, 0, 1)中的0表示除以 n-1，对应样本标准差。对 PCA 来说，用总体标准差还是样本标准差只会让特征值整体缩放，不影响特征向量方向，因此不必过度纠结。

2.2 用协方差矩阵的 eig 分解得到主成分得分

标准化完成后，下一步是计算协方差矩阵，并求其特征值和特征向量。数学上，协方差矩阵是

C = Xs' * Xs / (n-1)

在 MATLAB 中直接使用cov(Xs)更方便，它会自动执行中心化和除以 n-1。特征值表示对应特征向量方向上的方差，特征向量就是载荷向量，也就是主成分的方向。将数据投影到这些特征向量上，就得到主成分得分。

C = cov(Xs); % p×p 协方差矩阵 [V, D] = eig(C); % V 列为特征向量，D 为对角矩阵 [lambda, idx] = sort(diag(D), 'descend'); V = V(:, idx); % 按特征值从大到小排列载荷向量 score = Xs * V; % n×p 的主成分得分矩阵

这段代码的坑有两个。第一，eig返回的特征值和特征向量不保证按大小排列，必须手动排序。第二，特征向量的符号可以整体乘以 -1 而不改变方向，所以两次运行或与内置函数比较时，载荷列可能完全相反。这不影响降维结果，但如果要做严格的符号比较，需要统一符号规则，后面会专门说。

2.3 用 SVD 替代 eig：数值稳定性更高

在自写主成分分析源代码时，我更推荐使用奇异值分解而不是特征值分解。SVD 直接对数据矩阵做分解，不必先计算协方差矩阵。协方差矩阵的计算相当于对数据做了一次平方，当变量维度很高时，平方操作会放大浮点舍入误差，丢失小方差方向的信息。特别是在“变量数多于样本数”的 p >> n 场景下，协方差矩阵可能严重病态，eig的结果并不可靠。

SVD 的数学形式是

Xs = U * S * V'

其中 U 的列是左奇异向量，S 的对角元是奇异值，V 的列就是载荷向量，与协方差矩阵的特征向量一致。主成分得分可以直接用 U * S 得到。对应关系为：协方差矩阵的特征值等于奇异值的平方除以 n-1。

[U, S, V] = svd(Xs, 'econ'); % econ 是经济分解，节省内存 score_svd = U * S; lambda_svd = diag(S).^2 / (size(Xs, 1) - 1);

在 MATLAB 中，'econ'返回的矩阵尺寸尽量小。当 n ≥ p 时，U 是 n×p，S 和 V 是 p×p；当 n < p 时，U 是 n×n，S 是 n×n，V 是 p×n。无论哪种情况，score_svd列数都是 min(n, p)，正好对应最多能取到的主成分数量。

下面这张表总结了eig与svd在 PCA 实现中的差异：

理解了这些关系后，就可以开始写一个真正能复用的主成分分析源代码函数。

3. 用 MATLAB 写一个可复用的主成分分析源代码函数

3.1 函数输入输出与扩展参数设计

写源代码前要先定接口。我习惯把输入输出设计成和 MATLAB 内置pca接近但不完全一样，因为这样在迁移时直觉成本最低。

输入参数有三个：原始数据X，保留的主成分个数num_comp，以及一组可选的center/scale开关。输出参数包括主成分得分score、特征值latent、载荷矩阵coeff、训练集均值mu和标准差scalevec。后面两个参数是为了让新样本用同一个标准投到已训练的主成分空间里。

3.2 完整源代码与关键行解释

下面是完整的函数代码，可以直接复制保存为pca_manual.m：

function [score, latent, coeff, mu, scalevec] = pca_manual(X, num_comp, options) % 自写主成分分析源代码 % 输入： % X - n×p 数值矩阵，行为样本，列为变量 % num_comp - 保留的主成分个数，默认 min(n, p) % center - 是否去均值，默认 true % scale - 是否除以标准差，默认 false % 输出： % score - n×num_comp 主成分得分 % latent - num_comp×1 特征值，表示各主成分方差 % coeff - p×num_comp 载荷矩阵 % mu - 1×p 训练集均值 % scalevec - 1×p 训练集标准差，未缩放时全为 1 arguments X {mustBeNumeric} num_comp {mustBePositive, mustBeNumeric} = min(size(X, 1), size(X, 2)) options.center (1,1) logical = true options.scale (1,1) logical = false end [n, p] = size(X); Xc = X; mu = zeros(1, p); scalevec = ones(1, p); if options.center mu = mean(X, 1); Xc = Xc - mu; end if options.scale scalevec = std(Xc, 0, 1); scalevec(scalevec < eps) = 1; % 避免零方差噪声变量 Xc = Xc ./ scalevec; end [U, S, V] = svd(Xc, 'econ'); latent = diag(S).^2 / (n - 1); coeff = V; score = U * S; num_comp = min(num_comp, size(V, 2)); score = score(:, 1:num_comp); coeff = coeff(:, 1:num_comp); latent = latent(1:num_comp); end

代码逻辑并不复杂，但有几点值得说明。

第一个是arguments块。它从 MATLAB R2019b 开始支持，作用是在入口处做类型检查和默认值赋值。num_comp的默认值使用了size(X, 1)和size(X, 2)中的较小值，这对应 SVD 最多能分解出的主成分数量。如果你还在维护 R2018a 等旧版本，需要把arguments块改成nargin判断，这里不展开写。

第二个是options.center和options.scale是逻辑标量。在调用时直接写'center', true或"center", true都可以。MATLAB 会把名称-值对映射到options结构体里。

第三个是scalevec(scalevec < eps) = 1。当某个变量在所有样本上的取值完全相同，标准差就是 0，除以 0 会产生 NaN，后面的 SVD 全部失效。把标准差替换成 1 等价于不做缩放，同时保留这条变量。这个保护在真实脏数据里非常常见。

3.3 与内置 pca 函数对齐和符号修正

写完自写函数后，第一件事是拿它和内置pca对比，确认结果在数值上没有大偏差。对比代码很简单：

[score_pca, ~, ~] = pca(X, 'Center', true, 'Scale', true, 'NumComponents', 4); [score_manual, latent_manual, coeff_manual] = pca_manual(X, 4, 'center', true, 'scale', true);

内置pca的输出和自写函数的输出大概率会有符号差异。因为特征向量 $v$ 和 $-v$ 对应同一个特征值，两者都是合法解。如果强迫它们逐列完全一致，可以对每一列检查首元素符号：

for k = 1:size(coeff_manual, 2) if coeff_manual(1, k) * coeff_pca(1, k) < 0 coeff_manual(:, k) = -coeff_manual(:, k); score_manual(:, k) = -score_manual(:, k); end end

符号对齐的目的是方便调试，实际建模时不需要做。更值得对比的是latent和主成分得分之间的重建关系。按行对比score_manual和score_pca，数值的绝对值应该一致，哪怕符号相反。

4. 用自写源代码做主成分个数选择和可视化

4.1 在 fisheriris 数据上跑通全流程

为了验证函数可用性，直接在 MATLAB 自带的鸢尾花数据上跑一遍。这个数据集是 150×4 的矩阵，四个变量分别是花萼长度、花萼宽度、花瓣长度、花瓣宽度，类别三组。它足够简单，方便观察每个中间量的含义。

load fisheriris; X = meas; % 150×4，行是样本，列是变量 [score, latent, coeff, mu, scalevec] = pca_manual(X, 4, 'center', true, 'scale', true); explained = latent ./ sum(latent) * 100; disp(explained);

score的第一列是第一主成分得分，latent的第一个值是第一主成分的方差，explained表示每个主成分解释的方差比例。对这个标准化后的数据，前两个主成分合计大约解释了 95% 以上的方差，这说明把 4 维压到 2 维是合理的。

4.2 用累计方差贡献率和碎石图确定保留主成分个数

保留几个主成分并没有固定答案。常见做法是画累计方差贡献率图，观察曲线在哪个位置开始变平。下面这段代码会生成一张碎石图：

figure; plot(1:numel(explained), cumsum(explained), '-o', 'LineWidth', 1.5); xlabel('主成分编号'); ylabel('累计方差解释百分比'); grid on;

从工程经验看，累计贡献率达到 80% 以上就可以接受。如果目标是可视化，取前两到三个主成分。如果目标是做后续回归或分类，还要结合交叉验证结果决定，不能只看方差。

这里补充一个常见误用：很多人直接用latent > 1作为筛选条件，也就是 Kaiser 准则。这个准则只对标准化后的数据有参考意义，因为标准化后每个原始变量方差为 1，特征值大于 1 意味着该主成分解释的方差大于单个变量。但它只是一个经验门槛，不适用于所有场景。

4.3 把主成分得分和载荷画在同一张图里

降维后最常用的可视化是二维散点图。按类别着色可以看到数据分布是否可分。下面的代码使用gscatter，它需要 Statistics and Machine Learning Toolbox。如果没装，可以改用plot配合颜色向量替代。

figure; gscatter(score(:,1), score(:,2), species); xlabel('第一主成分'); ylabel('第二主成分'); grid on;

光有得分图还不够，最好把原始变量的载荷向量画上去。这样能直观看到每个原始变量对主成分的贡献方向。载荷向量从原点出发，在得分图上的坐标就是coeff(i,1)和coeff(i,2)。注意载荷的数值范围通常是 [-1,1]，为了和得分重叠显示，需要乘一个缩放系数。

hold on; scale_factor = 4; % 根据得分范围调整 for i = 1:size(coeff, 1) quiver(0, 0, coeff(i, 1) * scale_factor, coeff(i, 2) * scale_factor, 'k'); text(coeff(i, 1) * scale_factor, coeff(i, 2) * scale_factor, ... sprintf('var%d', i), 'FontSize', 10); end hold off;

载荷向量越长，说明该变量对当前主成分的影响越大；两个变量夹角接近 90 度，说明它们在这个主成分平面上近似正交，相关性弱。这种可视化能直接回答“哪些变量驱动了样本分离”的问题。

5. 把主成分分析源代码推进到重建、缺失值与批量场景

5.1 用主成分重建原始数据并计算误差

降维不是终点，很多场景还需要从主成分得分重构回原始数据。重建公式是原始标准化数据近似等于主成分得分乘载荷矩阵的转置，再乘以标准差并加回均值。

k = 2; X_hat = score(:, 1:k) * coeff(:, 1:k)' .* scalevec + mu; recon_error = mean((X(:) - X_hat(:)).^2);

这里scalevec就是第 3 章函数返回的训练集标准差。如果当时没有做标准化，scalevec全为 1，乘不乘都一样。重建误差可以两个方面用：一是作为自写源代码正确性的验证，保留全部主成分时误差应该接近 0；二是用来选择主成分个数，绘制“重建误差 vs 主成分个数”曲线，误差拐点往往对应合理压缩维度。

5.2 有 NaN 的矩阵要提前处理而不是硬算

SVD 遇到 NaN 会直接失败。自写主成分分析源代码没有内置缺失值处理能力，所以进入函数前必须清洗数据。最简单的做法是删除全为缺测的行，或者使用fillmissing插值：

X = fillmissing(X, 'linear');

但插值要谨慎。如果某个变量的缺失比例超过三成，插值结果会严重扭曲方差结构，这时不如直接删除该变量。真实项目中，建议先把缺失值情况做成一份报告，再决定用均值填充、线性插值还是行删除。不要指望 PCA 本身能神奇地应对不完整的矩阵。

5.3 训练集和测试集分离的投影写法

在交叉验证或上线部署时，一个常见错误是对训练集和测试集分别做标准化。正确做法是用训练集的mu和scalevec处理测试集，再乘训练好的载荷矩阵。下面这段代码给出了标准写法：

function test_score = pca_transform(newX, coeff, mu, scalevec) Xc = (newX - mu) ./ scalevec; test_score = Xc * coeff; end

这样写能保证测试样本和训练样本落到同一个主成分空间中。如果你从内置pca迁移到这个自写函数，记得把内置pca的第四输出mu和第二输出coeff对应传进来。调试时如果发现测试集投影结果与训练集分布明显脱节，第一步检查就是是否误用了测试集自己的均值去中心化。

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