简介：这份非负矩阵分解MATLAB代码包较为完整地涵盖了基础NMF及常见扩展与优化实现，适合机器学习、图像处理、文本挖掘和生物信息学方向的研究者、学生及工程人员使用。包内围绕V=WH的分解目标，提供了乘法更新、交替最小二乘等主流优化算法实现，并配有脑成像与文本主题分解的演示脚本和可视化图表，方便读者对照理解从初始化、迭代更新到收敛判断的完整流程。压缩包内共26个文件，核心为19个m源码文件，另有3个mat示例数据、2个fig结果图、1个txt说明文档及1个asv自动保存文件，rar包整体约490KB，结构紧凑、便于按需筛选。目前已有2031人学习下载，说明其在实际场景中具有一定参考价值。读者可借助这些代码直接运行示例、替换自有数据，也可在此基础上理解稀疏NMF、半监督NMF等变体的实现思路，适合用于算法研究、课程作业或项目原型验证。
1. 非负矩阵分解是什么：一个让结果可解释的矩阵拆解工具
作为工程师，我在图像处理、语音分离和文本主题建模里都遇到过同一类问题：手里有一堆非负数据，想拆成「几个组成部分 + 每组权重」，又希望拆出来的结果每项都能解释。非负矩阵分解（NMF）就是为这个需求设计的——它把一个大矩阵拆成两个小矩阵的乘积，约束条件是所有元素不小于零，代价是分解结果天然带「加性组合」的意义，不像 PCA 或 SVD 那样会出现正负抵消的难解释分量。这篇文章把一份能直接跑的完整 NMF MATLAB 代码拆开讲清楚：从目标函数、更新规则到秩选择、初始化、稀疏约束与稳定性验证，照着调参就能用在自己的数据上。适合刚接触矩阵分解的研究生，也适合需要把 NMF 落地到实际系统的工程师。
2. NMF 的目标函数与更新规则：为什么乘法更新是首选
2.1 两种目标函数：欧氏距离与 KL 散度怎么选
NMF 要解决的问题用一句话概括：给定非负矩阵 V（m×n），找非负矩阵 W（m×r）和非负矩阵 H（r×n），让 W×H 尽可能逼近 V。这里 r 远小于 m 和 n，所以本质是在做降维，只不过降维后的基向量和系数都被迫保持非负。
逼近程度怎么衡量，直接决定了更新规则的形式。最常用的是欧氏距离，也就是 Frobenius 范数的平方：min ||V - W×H||²。这个目标函数对 W 和 H 分别求梯度，梯度里会出现 V 和 W×H 的差，更新规则自然从「梯度下降」推出来。它的优点是数学上干净、收敛曲线好画，适合图像重构、数据压缩这类数值型任务。
另一种是 KL 散度：min Σ(Vij × log(Vij / (WH)ij) - Vij + (WH)ij)。它假设数据服从某种计数分布，对量级差异更敏感。比如文本词频矩阵、基因表达计数矩阵这类稀疏整数数据，用 KL 散度通常比欧氏距离效果好。选哪个，我一般看数据本身：如果元素值跨数量级大、稀疏度高，先试 KL 散度；如果是稠密连续值，直接用欧氏距离，代码和调参都简单。
2.2 乘法更新规则的直觉：除法找尺度，乘法保非负
有了目标函数，最优解可以通过梯度下降迭代求解，但普通梯度下降有两个毛病：步长难调，小了收敛慢、大了发散；更新后可能出现负值，还得每次强截到零，破坏了梯度方向的合理性。
经典的乘法更新规则避开了这两个问题。先看欧氏距离版本对 H 的更新：
H ← H .× (W'×V) ./ (W'×W×H + eps)
拆开看这个式子的结构。分子 W'×V 可以理解为「数据在基空间上的投影」，分母 W'×W×H 是「当前重构结果在基空间上的投影」。两者相除得到一个逐元素的缩放因子，H 的每一项都由它自己乘上这个因子得到。除法算出方向和步长，乘法保证非负——因为因子永远是正的。W 的更新完全对称：
W ← W .× (V×H') ./ (W×H×H' + eps)
这个规则最省心的地方是无需调学习率，每一步都在缩小目标函数，而且天然满足非负约束。实际实现时，分子分母里都要加一个极小的 eps（比如 1e-9），防止零除以零产生 NaN。这个 eps 不能加太大，否则收敛变慢，分解精度下降。
2.3 为什么在 MATLAB 里要自己写而不是直接用内置函数
MATLAB 自带 nnmf 函数，的确能跑，但它是个黑匣子：初始化策略固定、更新规则固定、迭代终止条件固定，一旦结果不理想，你能调的东西非常有限。我最早贪图省事直接用 nnmf，后来遇到一个语音分离的任务，需要基向量满足某种稀疏形态，内置函数给不了这种灵活性，只能自己写。
自己写还有一个实际好处：你能在每一轮迭代里拿到中间量。目标函数值、W 和 H 的变化量、基向量的稀疏度，这些信息对排查「为什么不收敛」「为什么结果不稳定」至关重要。内置函数只给你最终结果，中间过程全在底层。做工程排查问题时，看不到中间状态等于瞎猜。所以我的习惯是：先用自己写的版本跑通、摸清数据特性，再回头决定要不要用内置函数做性能优化。
3. 手写一份可运行的 NMF MATLAB 代码：从最小版到完整版
3.1 最小可用版：十几行跑通分解
先给一个能直接运行的乘法更新实现，输入 V、秩 r、最大迭代次数，输出 W 和 H。这个版本不做任何花哨处理，但足以验证算法通路。
function [W, H] = nmf_min(V, r, maxiter) % 最小版 NMF：欧氏距离 + 乘法更新 % V: m×n 非负矩阵 % r: 分解秩 % maxiter: 最大迭代次数 % 输出 W: m×r, H: r×n，均为非负矩阵 [m, n] = size(V); % 初始化：均匀随机数，加一个小量避免零值 W = rand(m, r) + 0.01; H = rand(r, n) + 0.01; eps_val = 1e-9; for iter = 1:maxiter % 更新 H：分子是 W' * V，分母是 W' * W * H H = H .* (W' * V) ./ (W' * W * H + eps_val); % 更新 W：分子是 V * H'，分母是 W * H * H' W = W .* (V * H') ./ (W * H * H' + eps_val); % 每轮对 W 做列归一化，防止尺度漂移 col_norm = sqrt(sum(W.^2, 1)) + eps_val; W = W ./ col_norm; H = H .* col_norm'; end end逻辑说明：循环里先更新 H 再更新 W，顺序固定，这是乘法更新规则的标准做法。更新公式直接对应 2.2 节推导出的形式，没有任何额外项。W 归一化这一步容易被忽略，但非常重要——如果不做，W 的列范数会在迭代中不断变化，H 也会跟着缩放，虽然乘积不变，但数值范围会漂移，最终可能导致精度损失甚至溢出。
参数说明：初始化用 rand(m,r)+0.01，加 0.01 是为了避免零矩阵死锁，这个细节在避坑章节还会展开。eps_val 是防除零的修正项，选 1e-9 比较稳妥。maxiter 一般从 100 起步，数据量大或精度要求高时加到 500，但不要盲目加到几千，后期迭代收益极低。
3.2 完整版：收敛判据与目标函数记录
最小版能跑通，但工程上不够用。实际项目中，你需要知道「什么时候该停」而不仅仅是「跑满迭代次数」。下面是完整版，增加了目标函数计算、相对变化判据、中间过程记录。
function [W, H, info] = nmf_full(V, r, opts) % 完整版 NMF：带收敛判据和历史记录 % opts 支持字段: % maxiter 最大迭代次数，默认 500 % tol 相对变化容忍度，默认 1e-4 % verbose 是否打印每轮信息，默认 false % 默认参数填充 if ~isfield(opts, 'maxiter'), opts.maxiter = 500; end if ~isfield(opts, 'tol'), opts.tol = 1e-4; end if ~isfield(opts, 'verbose'), opts.verbose = false; end [m, n] = size(V); % 推荐用基于 SVD 的初始化，见 4.2 节 [W, H] = nmf_init_nndsvd(V, r); eps_val = 1e-9; % 记录目标函数值和相对变化量 obj_hist = zeros(opts.maxiter, 1); rel_hist = zeros(opts.maxiter, 1); for iter = 1:opts.maxiter % 计算当前目标函数值 R = V - W * H; obj_hist(iter) = sum(R(:).^2); % 保存旧值用于相对变化判断 H_old = H; W_old = W; % 乘法更新 H = H .* (W' * V) ./ (W' * W * H + eps_val); W = W .* (V * H') ./ (W * H * H' + eps_val); % 归一化 W 并调整 H 保持乘积不变 col_norm = sqrt(sum(W.^2, 1)) + eps_val; W = W ./ col_norm; H = H .* col_norm'; % 计算相对变化量 delta = norm(W - W_old, 'fro') / (norm(W_old, 'fro') + eps_val) + ... norm(H - H_old, 'fro') / (norm(H_old, 'fro') + eps_val); rel_hist(iter) = delta; if opts.verbose fprintf('iter %4d, obj %.6e, rel %.6e\n', iter, obj_hist(iter), delta); end % 收敛判断：连续 3 轮相对变化小于 tol 才停 if iter > 3 && max(rel_hist(iter-2:iter)) < opts.tol break; end end % 截断历史记录到实际迭代轮数 info.obj_hist = obj_hist(1:iter); info.rel_hist = rel_hist(1:iter); info.iter = iter; end逻辑说明：目标函数值 obj_hist 在更新前计算，记录的是「当前 W、H 对应的重构误差」，这样你能画出随迭代下降的曲线，判断是否收敛。相对变化量 delta 综合了 W 和 H 的变化幅度，比只看单个矩阵更可靠。收敛判据用了「连续 3 轮都小于 tol」，避免某一轮偶然波动导致提前停止。
参数说明：tol 取 1e-4 是一个通用折中，追求精度可以压到 1e-6，但要接受迭代次数翻倍。verbose 在正式跑批时关掉，调试时打开，每轮打印能直观看到目标函数的下降趋势。这里相对变化用了两个项相加，实际操作中也可以只取 H 的变化量——如果 W 已经归一化稳定，H 的变化更能反映收敛状态。
3.3 用合成数据验证分解结果
代码写完，第一件事是用已知结构的合成数据验证，而不是直接上真实数据。我习惯这样构造：随机生成一个 W_true 和一个 H_true，稀疏且结构清晰，乘出 V，再加少量噪声，然后用 NMF 分解 V，看能否恢复出 W_true 和 H_true 的结构。
% 构造验证数据 rng(42); % 固定随机种子，保证可复现 m = 100; n = 80; r = 5; W_true = rand(m, r) .* (rand(m, r) > 0.6); % 60% 稀疏 H_true = rand(r, n) .* (rand(r, n) > 0.7); % 70% 稀疏 V = W_true * H_true + 0.01 * randn(m, n); % 加 1% 噪声 V = max(V, 0); % 保证非负 % 跑分解 opts.maxiter = 500; opts.tol = 1e-5; opts.verbose = true; [W, H, info] = nmf_full(V, r, opts); % 画目标函数收敛曲线 figure; semilogy(info.obj_hist); xlabel('迭代轮数'); ylabel('目标函数值'); title('重构误差随迭代下降'); grid on; % 检查恢复质量：W_true 与 W 的列能否一一对应 % 归一化后计算余弦相似度 Wn = W ./ sqrt(sum(W.^2, 1)); Wt = W_true ./ sqrt(sum(W_true.^2, 1)); sim = abs(Wn' * Wt); [max_sim, match] = max(sim, [], 2); fprintf('W 恢复质量（余弦相似度）：%.3f\n', mean(max_sim));逻辑说明：构造数据时用点乘做稀疏化，让 W_true 和 H_true 有清晰的零元素模式，这样验证才有意义——如果随机稠密矩阵也能被完美恢复，反而说明不了问题。加噪声后用 max(V,0) 保证非负，否则分解前就得先处理负值。恢复质量通过余弦相似度衡量，因为 NMF 的基向量可以任意换序，所以要取相似度矩阵每行的最大值，看平均恢复水平。
参数说明：rng(42) 的本质作用在验证场景里是把随机过程变成可复现的实验。稀疏比例 0.6 和 0.7 是经验值，太稀疏会让 V 变得几乎不可分，太稠密则验证不出算法的区分能力。噪声幅度 0.01 对应信噪比约 40dB，属于「干净但非理想」的场景，适合做第一轮验证。如果这个实验都恢复不好，先不要怀疑数据，回去检查代码。
4. NMF 的三个关键参数与选型经验：秩、初始化与稀疏约束
4.1 秩 r 怎么定：重构误差拐点法
NMF 的秩 r 是用户给的，不像 PCA 可以按特征值贡献率自动定。r 太小时，W×H 表达不了 V 的结构，重构误差大；r 太大时，分解会开始拟合噪声，而且基向量之间高度相关，分解结果不稳定。这是个典型的 bias-variance 权衡，实际定秩我用拐点法。
做法很简单：对同一个 V，跑一组不同的 r，比如 2 到 30，记录每个 r 收敛后的重构误差，然后画一条「r vs 误差」的曲线。曲线会有一个明显的肘部：r 小于某个值误差降得快，过了这个值误差下降变缓。这个拐点就是合理的秩。
% 用拐点法定秩 r_list = 2:30; err_list = zeros(length(r_list), 1); for i = 1:length(r_list) opts.maxiter = 300; opts.tol = 1e-4; [W, H, ~] = nmf_full(V, r_list(i), opts); R = V - W * H; err_list(i) = sum(R(:).^2); end % 画拐点曲线 figure; plot(r_list, err_list, 'o-'); xlabel('秩 r'); ylabel('重构误差'); title('选择拐点作为秩'); grid on; % 数值方法找拐点：计算连续差分，取最大变化点 diff_err = abs(diff(err_list)); [~, idx] = max(diff_err); fprintf('建议秩: %d\n', r_list(idx));逻辑说明：这里用 Frobenius 范数平方作为重构误差指标，和 NMF 的目标函数一致。数值上找拐点用的是「最大一阶差分」，即误差下降最快的位置，这个位置通常就是合理秩的上限。需要注意，拐点法给的是一个区间而不是一个精确定值，工程上会在拐点附近再试 2 到 3 个值，结合下游任务效果确定最终秩。
参数说明：r_list 的取值范围要覆盖数据的真实维度。图像数据我一般从 5 试到 50；文本主题建模从 10 试到 100；语音频谱数据从 10 试到 60。每档只跑 300 轮，不需要完全收敛，因为拐点位置对不同迭代轮数不敏感，这是个能省不少算力的经验。
4.2 初始化方法对比：随机、NNDSVD 与 k-means
初始化对 NMF 的影响经常被低估。乘法更新规则本质是局部最优求解器，随机初始化落在哪个局部最优完全看运气，结果可能差出一大截。我在一个图像分解项目里做过对比：同样的 V 和同样的 r，10 次随机初始化，重构误差相差超过 20%，这个波动在工程上是不可接受的。
随机初始化最省事，但只适合快速验证代码通路。NNDSVD 是更可靠的默认选择：先对 V 做奇异值分解，用 SVD 左奇异向量和右奇异向量的非负部分构造 W 和 H 的初始值。这样初始点就有全局结构信息，收敛快，结果稳定性好。实现并不复杂：
function [W, H] = nmf_init_nndsvd(V, r) % 基于 SVD 的非负初始化 % 取 SVD 的左右奇异向量的正部分构造初值 [U, S, Vt] = svds(V, r); S = diag(S); m = size(V, 1); n = size(V, 2); W = zeros(m, r); H = zeros(r, n); for j = 1:r % 取左奇异向量的正部 u_pos = max(U(:, j), 0); u_neg = max(-U(:, j), 0); % 取右奇异向量的正部 v_pos = max(Vt(j, :)', 0); v_neg = max(-Vt(j, :)', 0); % 用能量较大的那个方向作为初始基向量 if norm(u_pos)^2 >= norm(u_neg)^2 W(:, j) = sqrt(S(j)) * u_pos; else W(:, j) = sqrt(S(j)) * u_neg; end if norm(v_pos)^2 >= norm(v_neg)^2 H(j, :) = sqrt(S(j)) * v_pos'; else H(j, :) = sqrt(S(j)) * v_neg'; end end % 加小量避免零行（SVD 稀疏时可能出现零值） W = W + 1e-4; H = H + 1e-4; end逻辑说明：svds(V, r) 返回 top-r 奇异值和对应的左右奇异向量。每个奇异向量都有正部分和负部分，NNDSVD 的思路是保留能量更大的那一部分作为初始基向量，乘以对应奇异值的平方根。这保证了初始化有正确的尺度量级，同时满足非负约束。后面加 1e-4 是防御性处理，防止极稀疏情况下出现整列为零。
参数说明：这段代码是 NNDSVD 的简化实现，原始方法还包含对零元素列的再处理，实际中加小量已经够用。k-means 初始化是另一个选择：先对 V 的列做 k-means 聚类，聚类中心作为 W 的初始值，再投影求 H。这个方法在 V 的列天然分簇时效果好，但聚类本身要调 k 值和迭代次数，不如 NNDSVD 省心。我的默认选择是 NNDSVD，只在 NNDSVD 效果不理想时才考虑 k-means 或多次随机初始化取最优。
4.3 稀疏约束：什么时候加、加到什么程度
很多落地场景里，直接跑标准 NMF 得到的基向量是稠密的，每个分量都有一堆非零值，解释起来很费劲。比如做文本主题提取，一个主题下 5000 个词几乎全都有权重，这还能叫主题吗？这时候需要加稀疏约束，强制 W 或 H 的大部分元素往零靠。
加稀疏约束最直接的方式是在目标函数里加 L1 正则项。以收紧 H 为例，新的目标函数是：min ||V - W×H||² + λ×sum(H(:))。λ 是稀疏强度，越大 H 越稀疏，但重构误差也会越大。对应的乘法更新规则可以近似写成：
% 在完整版代码的 H 更新处替换成下面这行 H = H .* (W' * V) ./ (W' * W * H + lambda + eps_val);逻辑说明：分母上加一个正的 λ，等于把除法因子整体压低，H 的每一项都会比不加正则时更小，多次迭代后小值趋近于零。这是个工程上常用的近似处理——严格来说，L1 正则对应带次梯度的求解，但乘法更新里的这个近似在实践中效果不错，收敛也稳定。注意只加在分母上不会完全归零，需要配合阈值截断，比如每 10 轮把小于 1e-6 的 H 元素直接置零。
参数说明：λ 的量级取决于 V 的数值范围。我一般先跑一个不加正则的版本，记录目标函数值，然后从 λ = 0.01 × (初始目标函数 / sum(H(:))) 起步，逐步增大，每次观察 H 的非零比例和重构误差。一个经验参考：「重构误差上升不超过不加正则时的 10%，H 非零比例降到 40% 以下」通常是个好的平衡点。只管往大加 λ 会让分解变成纯粹的聚类，失去 NMF 的表示能力。
5. NMF 落地避坑：5 个常见问题与排查思路
5.1 迭代到一半 W 或 H 变成 NaN
现象：代码跑了几十轮，目标函数值直接变成 Inf 或 NaN，verbose 打印显示 obj 最后一次正常值后一下跳变。排查后发现 W 或 H 矩阵里有 NaN。
原因：最常见的是分子分母同时为零，在除零处产生 NaN。即使加了 eps，如果 V 是稀疏矩阵，某列全是零，V 里该列对应的 H 更新会分子为零、分母也为零，eps 只防了分母为零，但分子为零除以 eps 还是零，按理不该出 NaN——真正的坑是 eps 加得太小，而 W'WH 中某些元素是极小的正数，在浮点精度下被当成零，分子也是极小值，两者相除溢出。另一个原因是 V 里含负值或 NaN，随机初始化时 rand 不会产生 NaN，所以问题几乎都在输入数据。
解决：在分解前做一次数据清洗。具体做法是 V = max(V, 0) 强制截断负值，V(isnan(V)) = 0 清零 NaN，V(isinf(V)) = 0 清零 Inf。对稀疏列，可以在初始化后的第一次迭代前，给 V 整个加上一个极小偏移量，比如 V = V + 1e-4，保证没有全零列。这不是改数据语义，纯粹是数值稳定性处理。做完这两步，99% 的 NaN 问题都能解决。
5.2 分解结果每次跑都不一样
现象：同样的 V、同样的 r，连续跑三次，三次的 W 列顺序不同，数值也不同，重构误差有时差得还挺大。最典型的是在多次随机初始化时发现最优值和最差值的误差能差 20% 以上。
原因：乘法更新规则是非凸问题的局部寻优算法，初始化落在哪个吸引域就停在哪个局部最优。随机初始化每次点的位置不同，结果自然不同。这个现象不是代码 bug，是算法本质。
解决：分两层处理。第一层，默认初始化从随机换成 NNDSVD，这一步能把结果波动缩小一大半。第二层，如果 NNDSVD 还不够稳定，就做多次初始化取最优：比如跑 5 次不同的随机初始化，选目标函数值最小的结果。这个策略简单粗暴但非常有效。我见过一项评测显示多次启动取最优能把重构误差的波动从 20% 压到 3% 以内。注意每次跑之前固定 rng 种子，保证实验可复现，否则连「取最优」这个操作本身都不可复现。
5.3 目标函数不降反升
现象：verbose 打印的目标函数值在迭代到某个阶段后开始上升，或者呈锯齿状波动，完全没有单调下降趋势。还有人遇到这样的情况：前 50 轮正常下降，50 到 100 轮突然涨回去，然后又降下来。
原因：目标函数上升的最常见原因是 W 的列归一化和 H 的尺度补偿没有同步。如果在某一轮 W 归一化后没有正确调整 H，W 的列变小了，H 还没来得及放大，乘积整体变小，重构误差自然变大。另一个原因是学习步长过大——但在乘法更新规则里没有显式步长，这个现象通常来自归一化的时序错配。
解决：检查归一化代码，确保 W 归一化和 H 缩放补偿是原子操作，中间不能有目标函数计算。我遇到的情况里，有一版代码把目标函数计算放在了 W 归一化之后、H 补偿之前，导致那一轮算出的误差虚高。修正方式是先完成 H 更新，再完成 W 归一化和 H 补偿，最后算目标函数。如果排查代码没问题，那就是 tol 设得太紧，模型在收敛点附近震荡，适当放宽 tol 到 1e-4 或 1e-3 就能消掉锯齿。
5.4 稀疏数据下迭代极慢
现象：V 是一个稀疏矩阵，非零元素占比不到 5%，跑 300 轮耗时惊人，而且收敛曲线下降缓慢，像没吃饱饭一样一格格往下挪。
原因：NMF 的乘法更新是对稠密矩阵运算，每一步 W'WH 都是稠密矩阵乘。稀疏 V 在计算分子 W'*V 时可以通过稀疏乘法加速，但分母 W'WH 里 W'*W 是稠密的，乘以 H 也是稠密的，计算量跟 V 的稀疏度无关。换句话说，V 再稀疏，每一步迭代的代价都按满矩阵算，自然慢。
解决：有两招。第一招，利用结构分解稀疏乘法：分子 W'*V 改用 spfun 或手动索引只算非零位置，分母保持原样——这能省掉分子的计算量，通常占到总耗时三到四成。第二招，降低迭代精度需求，稀疏数据的目标函数下降本来就平缓，tol 用 1e-3 而不是 1e-4，能少跑几十轮。如果数据量级上万维，考虑先做降维预处理再用 NMF，比如先用随机投影把维度压到 2000 以内，再跑分解。这一步在文本数据上效果尤其明显。
5.5 基向量高度相似，分解失去意义
现象：分解出的 W 各列之间余弦相似度很高，比如两列的相似度超过 0.9，看起来像是同一个模式被拆成了两份。H 里对应的行也相似，整个分解结果无法解释。
原因：秩 r 定得太大了。模型有太多自由度去拟合 V，但 V 本身的可分性没那么高，多余的自由度就会分解出「同一个基向量的微扰版本」。这和 PCA 中主成分末尾出现噪声分量是同一逻辑，只是 NMF 没有特征值排序，你一眼看不出来哪些分量是多余的。
解决：先按 4.1 节的拐点法重新定秩，通常把 r 降到拐点附近就能消除相似基向量。但也要注意另一种情况：基空间本身就高度相关。比如图像数据里不同物体有相似纹理时，NMF 会合理地给出相近基向量。这时先用原始 V 看 W 的相似度矩阵，再对一个按经验降秩后的 V 重新分解，对比两组 W 的相似度是否有显著变化，就能判断是秩的问题还是数据本身的问题。如果降秩后相似度依然很高，那就是数据可分性差，考虑改用投影梯度法或加入正交约束的变体。
6. 分解质量的验证技巧：用重构误差、稀疏度与稳定性说话
6.1 三条曲线判断收敛
代码跑完别只看最终结果，先把三条曲线画出来：目标函数值曲线、相对变化量曲线、W 列归一化后的稀疏度曲线。第一条看整体下降趋势是否正常；第二条看迭代在什么位置真正稳定；第三条看基向量是否从一开始就保持了合理的稀疏形态。
如果目标函数曲线尾部还带着明显下降斜率就停了，说明 tol 设得太松；如果稀疏度曲线在迭代后期还在持续变化，说明模型还在重组基向量，离稳定还远。这三个信号都正常后，再谈结果可不可用。可视化这一步花的五分钟，能省掉后面分析垃圾结果的数小时。
6.2 稳定性检验：多次分解的基向量相似度
一次分解结果好不代表算法可靠，尤其在上生产环境之前，要做稳定性检验。做法是固定参数，换多个随机种子分别初始化，跑多次分解，然后计算这几次 W 两两之间的最大余弦相似度，取平均值。这个均值越接近 1 越稳定。
我一般跑 10 次，算出一个「平均结构相似度」。大于 0.9 说明分解结构稳定，结果可信；0.6 到 0.9 说明结果有波动，需要检查初始化策略和秩的选择；小于 0.6 说明分解基本不可复现，这样的模型无法上线。这个指标在我做过的项目里比任何单一重构误差都更能反映分解质量。
6.3 一个值得保留的检查习惯
最后提一个我长期保留的检查习惯：每次跑完 NMF，把 W 的每列单独画出来或打印成热图，肉眼扫一遍再决定是否进入下游。这个习惯救过我很多次——数值指标全部正常时，基向量可能仍然是噪声，只有人眼对这个领域足够敏感，能一眼看出「这一列不合理」。用这个习惯配合前面的曲线和稳定性检验，能覆盖绝大多数 NMF 部署质量问题，把玄学变成可控的工程过程。希望帮到你。
本文还有配套的精品资源，点击获取
