非负矩阵分解Matlab完整实现：从乘法更新到稀疏正则与调参避坑

简介：一套面向MATLAB使用者的非负矩阵分解完整代码包，适用于机器学习、图像处理、文本挖掘和生物信息学方向的研究者与学生。代码从非负矩阵分解的基本原理出发，涵盖随机初始化、乘法更新规则、交替最小二乘等主流迭代优化策略，同时加入稀疏NMF、半监督NMF等扩展方法，便于读者在理解算法细节的基础上，将已有数据快速代入实验。压缩包内共26个文件，以19个m源码文件为核心，另有3个mat示例数据、2个fig结果图形以及说明文档，整体约490KB，模块划分清晰，适合按需调用和二次开发。附带脑影像主题分解、文本主题发现等演示脚本，可观察不同场景下算法收敛与分解效果。已有2031人学习下载，对需要实现NMF、开展降维或特征提取实验的开发者具有直接参考价值。

1. 非负矩阵分解的Matlab代码不难找，但一套能直接跑通又能改的全套实现很少见

非负矩阵分解（NMF）的Matlab代码在各类论坛和工具包里都能找到片段，但要么只有迭代两行公式，要么封装得太重看不懂内部逻辑。这篇把一套完整的非负矩阵分解Matlab代码从头到尾拆开：从目标函数和初始化写起，给出可直接运行的主函数，再扩展出带稀疏约束的正则版本，最后用合成数据验证实现有没有写错。适合正在做数据挖掘实验、需要把NMF当基线的同学，也适合想搞清楚乘法更新每一步在算什么、遇到NaN和不收敛时知道去哪查的人。读完你能拿到一份属于自己的、可复现可改参数的NMF工具箱，而不是只能复制粘贴的一次性脚本。

2. NMF的目标函数与乘法更新规则：为什么迭代公式长这样

2.1 两种目标函数：Frobenius范数与KL散度怎么选

NMF把m×n的非负矩阵V分解成m×k的W和k×n的H，两个因子都非负，k远小于m和n。这里最关键的选择是“怎么衡量重构误差”。最常见的两个目标函数：

• Frobenius范数版本：min ||V-WH||_F^2，对应高斯噪声假设。计算简单，误差有明确的“能量”含义，适合连续测量数据。
• 广义KL散度版本：min Σ(V_ij log(V_ij/(WH)_ij) - V_ij + (WH)_ij)，对应泊松噪声假设，适合计数数据、文档词频、图像像素强度这类非负整数。

我一般默认写F范数版本，原因有三个：第一，它的梯度表达式最干净，乘法更新公式只有两项乘法一项除法；第二，重构误差可以直接用相对误差百分比报告，和业务方解释时不用引入信息散度概念；第三，绝大多数论文的baseline实验用的都是F范数，对比实验更方便。如果你的数据是稀疏计数矩阵，且每个元素都不大，比如文档-词频矩阵，可以换成KL版本，收敛行为和误差曲线的解释逻辑会不一样，排查时别拿F范数的经验硬套。

2.2 乘法更新公式的直觉：不需要设置学习率的关键

如果用梯度下降解带非负约束的优化问题，每次迭代要做“减去学习率乘以梯度”的操作，结果很容易出现负值，需要额外投影到非负象限，学习率也跟着麻烦。乘法更新的聪明之处在于，把梯度下降改写成逐元素乘法的形式：

H ← H ⊙ (W^T V) ⊙ (W^T W H + ε)^(-1)

W ← W ⊙ (V H^T) ⊙ (W H H^T + ε)^(-1)

其中⊙表示逐元素乘法，除法也是逐元素。这个公式看似天降，实际可以证明乘性更新保证目标函数单调不增。工程上理解更直接：初值非负时，乘法和除法都不会产生负数，非负约束自动满足，而且不需要手动调学习率。括号里的比值可以当作“自适应学习率矩阵”：

• 当(W^T W H)_ij大于(W^T V)_ij，说明H_ij当前估计偏大，分母大分子小，比值小于1，H_ij被缩小；
• 反过来比值大于1，H_ij被放大。

这个自适应特性让乘法更新在实际使用中相当稳，也是NMF相比其它非负优化算法更普及的原因。对不熟悉公式推导的读者，只需要记住：分子是“数据告诉我们的方向”，分母是“当前模型已有的方向”，两者相除就是对模型的修正。

2.3 先算公共中间量：每轮迭代的计算流程

一次迭代中有三个关键的矩阵乘法组，表达式里W^T W H、W H H^T都会反复出现。我在实现里总是先算中间量再逐元素更新，避免重复计算：

% 一次完整迭代的骨架，后续章节会封装成完整函数 WtV = W' * V; % k×n WtWH = W' * (W * H); % k×n，注意先算 W*H 再左乘 W' % 更新 H H = H .* WtV ./ (WtWH + epsVal); VHt = V * H'; % m×k WHHt = W * (H * H'); % m×k % 更新 W，注意这里 H 已经更新过 W = W .* VHt ./ (WHHt + epsVal);

逻辑说明：先更新H时使用当前W，再更新W时使用刚更新的H，这是标准的交替更新策略。两个中间量WtV和VHt主导了计算量，整体复杂度是O(mnk)。代码里的epsVal是防除零常数，取1e-12而不是1e-16，原因是double精度下分母接近1e-16时，除法结果会剧烈震荡甚至翻车。

参数说明：epsVal取太小会出问题，在第5章展开讲；如果你用的是单精度计算，建议直接取1e-6。这个骨架的效率关键在于避免在循环里重复构造W*H，有人会把W'*W先算成k×k矩阵再乘H，这样数学等价，但内存访问模式不如“先算WH再左乘W'”直观，数据量大时速度差异明显。

3. 完整可运行的Matlab实现：主函数、初始化、迭代与收敛控制

3.1 函数接口设计：输入输出参数怎么定

完整代码不能只写迭代公式，还要把初始化、收敛判定、耗时统计这些都放进一个可复用的函数里。我习惯用一个主函数加一个初始化子函数，主函数接口如下表：

输出设计为W、H、histErr、elapsed四个变量，其中histErr是每次迭代的相对重构误差序列，用来画收敛曲线判断终止阈值是否合理；elapsed返回耗时，方便对比不同初始化或不同秩的性能。函数内部先把opts里的字段逐个取出并赋予默认值，这样调用方只写nmf_fro(V, k)就能跑通，老手又可以通过opts结构体自由改参数。

3.2 NNDSVD初始化：让起点离好解更近

随机初始化简单但方差大，同一个数据集跑五次可能得到五种结果。NNDSVD的核心思想是用SVD的低秩逼近结果来构造非负初始值：取V的前k个奇异值分解，用奇异向量的绝对值构造W和H，这样初始点已经捕捉了大体结构，后续迭代更稳定。

function [W, H] = init_nndsvd(V, k, seed) % 基于SVD的非负初始化，简化的绝对投影版本 % 数据有效秩不足k时回退到随机初始化 rng(seed); [U, S, VT] = svds(V, k); % V ≈ U * S * VT' if size(S, 1) < k % 数据有效秩不足k时，回退到随机初始化 W = rand(size(V,1), k) + 1e-4; H = rand(k, size(V,2)) + 1e-4; return; end W = abs(U) * sqrt(S); % m×k，取U绝对值乘奇异值平方根 H = sqrt(S) * abs(VT'); % k×n，取VT转置的绝对值 % 处理全零列/全零行，防止后续除零 zeroColW = sum(W, 1) < 1e-12; zeroRowH = sum(H, 2) < 1e-12; if any(zeroColW) W(:, zeroColW) = rand(size(W,1), sum(zeroColW)) + 1e-4; end if any(zeroRowH) H(zeroRowH, :) = rand(sum(zeroRowH), size(H,2)) + 1e-4; end end

逻辑说明：svds返回的VT是n×k的右奇异向量矩阵，转置后得到k×n；乘上sqrt(S)让能量分布到两个因子上，避免W过大H过小或反过来。标准NNDSVD会保留正部和负部分开构造，我这里用abs投影的简化版本，优点是代码短、不会因为分量符号相消产生零列，缺点是理论保证比原版弱一点，但作为乘法更新的起点实测效果足够好。

参数说明：加上的1e-4是为了让零列变成全正列，防止第一步迭代就除零；rng(seed)放在函数开头而不是调用处，是为了让初始化结果可以复现。如果你的矩阵是稀疏的，svds内部会自动利用稀疏结构，比直接svd快很多。

3.3 主迭代循环：乘法更新与提前终止

主函数代码是这个方案的核心，完整贴出来可以直接保存成nmf_fro.m使用。

function [W, H, histErr, elapsed] = nmf_fro(V, k, opts) % 非负矩阵分解，Frobenius范数乘法更新 if nargin < 2, error('至少需要V和k'); end if nargin < 3, opts = struct(); end % 参数默认值 maxIter = 200; tol = 1e-4; init = 'nndsvd'; seed = 0; verbose = false; keepStep = true; if isfield(opts,'maxIter'), maxIter = opts.maxIter; end if isfield(opts,'tol'), tol = opts.tol; end if isfield(opts,'init'), init = opts.init; end if isfield(opts,'seed'), seed = opts.seed; end if isfield(opts,'verbose'), verbose = opts.verbose; end if isfield(opts,'keepStep'), keepStep = opts.keepStep; end if any(V(:) < 0), error('V必须非负'); end epsVal = 1e-12; % 初始化 if strcmp(init, 'nndsvd') [W, H] = init_nndsvd(V, k, seed); else rng(seed); W = rand(size(V,1), k); H = rand(k, size(V,2)); end normV = norm(V, 'fro'); prevErr = inf; histErr = zeros(maxIter, 1); tStart = tic; for iter = 1:maxIter % 更新H H = H .* (W' * V) ./ (W' * (W * H) + epsVal); % 更新W，用更新后的H W = W .* (V * H') ./ (W * (H * H') + epsVal); err = norm(V - W*H, 'fro') / normV; relChange = abs(err - prevErr) / max(err, 1e-12); prevErr = err; if keepStep, histErr(iter) = err; end if verbose fprintf('iter %4d, relErr %.6f, relChg %.2e\n', iter, err, relChange); end if relChange < tol, break; end end if keepStep histErr = histErr(1:iter); % 截断到实际迭代次数 else histErr = err; end elapsed = toc(tStart); end

逻辑说明：主循环严格按“先H后W”的顺序，H用当前W更新，W用新H更新，这是标准交替更新。relErr用相对重构误差归一化，除以normV后不同数据集之间可以横向比较。relChange衡量的是本次误差相对上次的变化幅度，当它小于tol时认为收敛，提前退出，避免白跑满maxIter。

参数说明：tol=1e-4适合大多数实验场景，如果你要的是高精度分解，设到1e-6，但迭代次数会显著增加；verbose打开后能看到迭代进程，我建议第一次跑数据时先开verbose，确认误差曲线形态正常再关掉。histErr如果keepStep=false，返回最终误差；如果true，返回的是截断后的收敛曲线。

3.4 调用示例：跑一个带噪声的合成数据

用一个简单脚本演示怎么调用、怎么画误差曲线，验证代码能跑通：

% 生成测试数据，真实秩为3 rng(42); Wt = rand(50, 3); Ht = rand(3, 200); V = Wt * Ht + 0.01 * randn(50, 200); V(V < 0) = 0; opts = struct(); opts.maxIter = 500; opts.tol = 1e-6; opts.init = 'nndsvd'; opts.seed = 42; opts.verbose = true; [W, H, histErr, elapsed] = nmf_fro(V, 3, opts); figure; semilogy(histErr, 'LineWidth', 1.5); xlabel('迭代次数'); ylabel('相对重构误差'); title('NMF收敛曲线'); grid on; fprintf('耗时 %.2f秒，最终误差 %.6f\n', elapsed, histErr(end));

逻辑说明：生成数据时故意用随机矩阵相乘作为已知真值，再用噪声污染，这样后面可以验证恢复效果。画误差曲线用semilogy，因为乘法更新前期误差下降快、后期慢，线性坐标会看不清尾部的缓慢下降段。

参数说明：这里的maxIter设到500是给足余量；seed=42保证每次跑的初始化都一样。V(V<0)=0把噪声产生的负值截断掉，保证输入满足非负条件。跑完观察两条关键信息：最终误差量级是否接近噪声水平（0.01的噪声对应相对误差约0.05以内），以及是否有迭代卡住的情况。

4. 从基础版到正则版：给NMF加稀疏约束的代码改动方案

4.1 什么场景需要加稀疏约束：当分解结果“糊成一团”时

基础NMF在数据本身缺乏明显结构时，分解出的W和H常常每个分量都混着多个主题，没有清晰语义。比如做文本主题分解，不加正则时每个主题都包含大量共有词；做音频源分离，基向量会变成混合音色。此时给目标函数加一项稀疏惩罚，让H的每一列只倾向于用少量基向量表示，分解才具有可解释性。稀疏约束的本质是对比“用很多小系数”和“用少数大系数”，后者在L1惩罚下更便宜。

4.2 正则应该加在H还是W：先想清楚你要什么

L1惩罚加在H上，得到的是“稀疏表示”——每个样本只用少数基向量组合，这是字典学习和源分离的常见设置。L1惩罚加在W上，会让基向量本身稀疏，比如某些基只在个别特征上有非零值，这在特征筛选场景有意义，但更容易因W过稀疏导致重构误差飙升。如果要限制W的尺度，更稳妥的是加L2惩罚，防止W元素无限增大。我的默认组合是：对H加L1稀疏，对W加L2平滑，兼顾可解释性和数值稳定性。

4.3 正则化后的更新公式与完整代码改动

目标函数变为||V-WH||_F^2 + αΣH_ij + β||W||_F^2。对H的L1项求梯度得到全1矩阵，所以更新H时分母额外加上α；W的L2项梯度是2W，所以更新W时分母加上βW：

function [W, H, histErr] = nmf_fro_reg(V, k, alpha, beta, opts) % 带稀疏约束的NMF % 目标: ||V-WH||_F^2 + alpha*sum(H(:)) + beta*||W||_F^2 if nargin < 5, opts = struct(); end if nargin < 4, beta = 0; end if nargin < 3, alpha = 0; end % 沿用第3章的默认参数设置逻辑，此处省略重复代码 optsDefault.maxIter = 200; optsDefault.tol = 1e-4; optsDefault.init = 'nndsvd'; optsDefault.seed = 0; if ~isfield(opts,'maxIter'), opts.maxIter = optsDefault.maxIter; end if ~isfield(opts,'tol'), opts.tol = optsDefault.tol; end if ~isfield(opts,'init'), opts.init = optsDefault.init; end if ~isfield(opts,'seed'), opts.seed = optsDefault.seed; end if strcmp(opts.init, 'nndsvd') [W, H] = init_nndsvd(V, k, opts.seed); else rng(opts.seed); W = rand(size(V,1), k); H = rand(k, size(V,2)); end normV = norm(V, 'fro'); prevErr = inf; histErr = zeros(opts.maxIter, 1); epsVal = 1e-12; for iter = 1:opts.maxIter % 更新H：分母多一项 alpha，对应L1稀疏惩罚 H = H .* (W' * V) ./ (W' * (W * H) + alpha + epsVal); % 更新W：分母多一项 beta*W，对应L2权重衰减 W = W .* (V * H') ./ (W * (H * H') + beta * W + epsVal); err = norm(V - W*H, 'fro') / normV; relChange = abs(err - prevErr) / max(err, 1e-12); prevErr = err; histErr(iter) = err; if relChange < opts.tol, break; end end histErr = histErr(1:iter); end

逻辑说明：这个函数和基础版唯一的差别就是更新公式分母加了alpha和beta*W，数学上对应把正则项梯度并入原目标函数梯度。alpha和beta是非负超参数，为0时完全退化为基础版，方便你对照测试正则项是否生效。

参数说明：alpha的单位取决于H的数值范围，H在1左右时alpha从0.01起步试；beta通常比alpha小一个量级，因为W的L2项作用在整矩阵上，太大会让W被压扁、重构误差爆炸。正则版收敛速度通常比基础版慢，maxIter可以放宽到500。

4.4 alpha怎么调：先扫对数网格再看两项指标

调alpha不能只看重构误差，因为加了正则后重构误差必然上升。正确做法是同时观察重构误差和稀疏度两个指标，找“误差上升可接受、稀疏度明显提高”的区域：

alphaList = [0, 1e-3, 1e-2, 1e-1, 1, 1e1]; opts.maxIter = 300; opts.init = 'nndsvd'; opts.seed = 1; for i = 1:length(alphaList) [W, H, ~] = nmf_fro_reg(V, 5, alphaList(i), 0, opts); reconErr(i) = norm(V - W*H, 'fro') / norm(V, 'fro'); thr = 1e-4 * max(H(:)); sparseRate(i) = sum(H(:) < thr) / numel(H); fprintf('alpha=%8.3f reconErr=%.4f sparsity=%.2f%%\n', ... alphaList(i), reconErr(i), 100*sparseRate(i)); end

逻辑说明：稀疏度统计H中所有小于阈值（H全局最大值的万分之一）的元素占比。alpha=0时的稀疏度是基准，一般很低；随着alpha增大稀疏度上升，重构误差也上升。选alpha的准则是“稀疏度曲线从陡峭变平缓的拐点附近”，或者直接选业务能接受的误差上限对应的alpha。

参数说明：阈值1e-4*max(H(:))是相对阈值，避免不同量级数据导致误判。这个扫描过程一定要在固定seed下进行，否则不同alpha之间混入初始化随机性，曲线会抖得没法看。

5. NMF避坑与排查实录：初始化、除零、收敛与内存的五个坑

5.1 迭代到一半出现NaN：eps太小与初始值里有零

现象：跑上几十次迭代后矩阵里突然冒出NaN，之后所有更新全变NaN，误差曲线断掉。原因有两个：一是epsVal取太小，比如取到1e-16，当分母里存在小于机器精度的数时除出来是无穷或NaN；二是初始化阶段有全零列，第一步更新就让分母出现零。解决：把epsVal固定为1e-12，单精度数据则取1e-6；初始化函数里做全零列检测，给零列注入随机正数。另外一个隐蔽来源是V里有NaN，输入检查务必把isnan(V)一并拦截。

5.2 重构误差震荡不下降：数据尺度与目标函数错配

现象：误差曲线不是平滑下降，而是小范围抖动甚至上升，relChange始终达不到设定阈值，直到跑满maxIter被截断。原因：乘法更新对F范数在理论上保证单调不增，但前提是目标函数定义与数据分布匹配。如果V的行之间量级差很大，比如某几行是上万、其余是零点几，F范数会把优化重心压在大尺度行上，小尺度行拟合得很差，甚至数值上出现更新震荡。解决：对V做预处理，每行除以该行的最大值，或者做min-max归一化，分解完再把尺度乘回去。连续性数据用F范数，计数数据改KL版本，不要在泊松数据上硬套F范数然后怪算法不收敛。

5.3 换一次随机种子结果差很多：未固定seed与初始化太差

现象：同一份数据、同一个k，换一个随机种子跑出来的W和H差异很大，重构误差也有明显差别。原因：NMF目标函数非凸，随机初始化会落到不同局部解。NNDSVD初始化能缓解但不能完全消除。解决：固定seed做实验，这是最基本的；对关键结论要跑5-10个seed，报告中位数和方差；对比实验时所有方法用同一组seed，保证可复现。如果某个seed下误差明显偏高，可以先看是不是初始化零列没处理干净。

5.4 加正则后误差暴涨且W失去意义：正则对象和系数选错

现象：alpha设到1甚至更大，重构误差直接从0.05涨到0.5，W列归一化后看每个基向量都是相同的稀疏模式，毫无判别力。原因：L1加的力度过大时，H被压成零矩阵，W只剩背景；更隐蔽的是正则加错对象，比如对W加L1导致基向量全变稀疏碎片。解决：把alpha从0.001开始做对数扫描，找到误差拐点；想约束W的尺度就用L2而不是L1；正则版要和基础版在同一个seed下对比，排除初始化干扰。

5.5 数据一大就卡死或内存溢出：中间矩阵膨胀与稀疏丢失

现象：V只有5000×10000就报内存不足，或者迭代一次要几十秒。原因：m×n的V不是瓶颈，真正吃内存的是W'V、VH'这类中间结果，虽然尺寸只有k×n和m×k，但如果k取得较大或者误用了full()把稀疏矩阵变稠密，内存会成倍翻。再一个常见误用是在循环里反复计算W*H，每次迭代多一次O(mnk)乘法。解决：V用稀疏矩阵存储，乘法更新里保持稀疏；中间量按第2.3节的方式先算公共项；循环里避免任何多余矩阵乘。另外一个取巧办法：如果m远大于n，可以分块处理样本维，每次迭代只取一批行做局部更新，这在流式场景里也适用。

6. 合成数据验证与批量选秩：NMF代码上生产前必做的两件事

6.1 用已知的W和H反推V，验证你的实现没有bug

NMF的分解结果有置换和缩放自由度，直接拿恢复的W和H跟真值逐元素对比必然对不上。验证的正确思路是构造已知的W_true、H_true，相乘得到V，跑完分解后检查乘积WH是否近似等于V，以及能否通过列匹配找到与真值对应的基向量。

rng(1); m = 60; n = 100; k = 3; W_true = rand(m, k); H_true = rand(k, n); V = W_true * H_true; V = V + 0.001 * randn(size(V)); V(V < 0) = 0; opts.seed = 1; opts.init = 'nndsvd'; opts.maxIter = 500; opts.tol = 1e-8; [W, H, ~] = nmf_fro(V, k, opts); % 第一阶段：重构误差必须接近噪声水平 reconErr = norm(V - W*H, 'fro') / norm(V, 'fro'); fprintf('重构误差: %.6f\n', reconErr); % 第二阶段：列匹配，判断恢复出的基方向对不对 Wn = W ./ sqrt(sum(W.^2, 1)); Wt_n = W_true ./ sqrt(sum(W_true.^2, 1)); for j = 1:k sim = abs(Wt_n(:, j)' * Wn); [best, idx] = max(sim); fprintf('真实第%d列 匹配到 W 第%d列，|内积|=%.4f\n', j, idx, best); end

逻辑说明：如果实现正确，reconErr应该在0.001这个噪声水平附近，比如0.02以内；列匹配的|内积|应当接近1，说明恢复的基方向与真值一致。如果reconErr很小但列匹配失败，多半是k取得不对或者数据本身的基方向互相纠缠；如果reconErr本身偏大，需要回查迭代是否收敛、初始化是否出错。

6.2 批量扫rank找拐点，以及一个关于seed和记录的小习惯

选k是NMF落地时最常纠结的问题。一个可靠的办法是扫不同rank，记录每个rank下的最小重构误差，误差曲线的肘部通常就是合理的秩。因为NMF对同一rank也有多个局部解，每个rank重复跑几次取最小值，避免被初始化带偏。

rmax = 10; repNum = 5; bestErr = zeros(rmax, 1); for r = 1:rmax e = zeros(repNum, 1); for rep = 1:repNum opts.seed = rep; [W, H, ~] = nmf_fro(V, r, opts); e(rep) = norm(V - W*H, 'fro') / norm(V, 'fro'); end bestErr(r) = min(e); end figure; plot(1:rmax, bestErr, 'o-'); xlabel('rank k'); ylabel('最小相对重构误差'); grid on;

我自己的习惯是：每一轮实验前先跑合成数据验证，把“实现正确”这件事钉死，之后再调参、改正则、换数据，翻车时先查数据不查算法。固定seed并记录每次实验的seed、rank、alpha这些参数，能帮你省掉大量“上次那个结果是怎么跑出来”的后悔时间。NMF这个方向值得投入，它不像深度学习那样需要大规模算力，一套完整代码几百行就能覆盖从原理到实验验证的整个闭环。希望帮到你。

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