企业官网建设流程全解析

搞稀疏表示这个方向也有一阵子了。最初手里那个信号恢复Demo是用L1正则化做的，正经高斯随机字典下表现还行，但一旦把字典换成过完备DCT或者让列与列之间有相关性，误差就开始飘，支撑集也经常认错。后来我把整个方案切到贝叶斯视角，用稀疏贝叶斯学习（SBL）重新写了一遍，效果稳了不少，代价是前前后后花了一晚上把公式重新推了一遍。这篇实践漫谈就想把这段经历做个系统记录，围绕基于贝叶斯方法的稀疏表示学习展开，重点讲清楚为什么稀疏表示问题天然适合贝叶斯处理、SBL的分层先验和证据最大化到底在干什么，以及MATLAB R2018上落地时那些论文里不会写的数值细节和调参经验。想用MATLAB快速搭一个能跑的贝叶斯稀疏恢复框架的同学，应该能从里面直接抄到能用的代码和参数。

1. 折腾稀疏表示的那些日子：为什么最后停在贝叶斯这条路上

1.1 稀疏表示要解决的到底是什么

稀疏表示问题的场景其实很直白：我们有一个观测向量y，一个字典矩阵A，假设y可以近似表示成字典列的线性组合，即 y = A x + 噪声，其中系数向量x里大部分分量是零或者接近零。任务就是从少数观测里把x恢复出来。压缩感知、图像去噪、阵列信号处理、甚至某些反演问题，最后都会落到这个模型上。

这个问题的难点在于，当字典A的行数N远小于列数M时，方程是欠定的，直接求最小二乘解没有任何唯一性。传统做法是加一个稀疏惩罚，比如让x的L1范数尽量小，得到类似 min ||y - Ax||² + λ||x||₁ 的优化问题，也就是Basis Pursuit或LASSO这一类。我最初用这类方法时，最大的麻烦在于选那个正则化系数λ。λ小了不够稀疏，λ大了把真实信号也压没了，交叉验证又要多跑一堆实验。而字典原子相关性一高，L1解还容易出现支撑集抖动，明明该是第4列起作用，结果第5列甚至第3列也分到不少能量，导致恢复出来的信号形态不对。

当然也可以走贪心路线，比如OMP这类匹配追踪算法，一步步挑和残差最相关的列。这种思路简单直观，但每一步的选择都是“现在看起来最优”，缺乏全局视角，一旦某一步选错列，后面基本没有纠错机制。尤其在低信噪比场景，OMP很容易把噪声原子当成真实支撑集的一部分，后面再怎么迭代也拉不回来。

1.2 传统方法与贝叶斯方法的思路分岔

我后来意识到，上面这些方法的共同问题是：它们都在找一个“最优解”，但这个解本身没有任何不确定性信息。换句话说，它们只回答“x最可能是什么”，不回答“x可能是什么范围”。而现实中噪声是客观存在的，观测数据有限，我们对x的认知本来就该是一个分布，而不是一个点。

贝叶斯方法的思路完全不同。它把系数x、噪声精度、甚至稀疏程度都当成随机变量，先给它们施加合理的先验分布，然后通过观测数据y计算后验分布。在稀疏表示这个具体问题上，贝叶斯处理有几个天然优势：第一，稀疏性不是靠外部正则化系数强加进来的，而是通过分层先验从数据中自适应学出来的，省掉了交叉验证最头疼的环节；第二，得到了后验均值之外还有后验方差，可以直接给出每个系数的置信区间；第三，对噪声的处理是概率意义上的，低信噪比下不容易出现“硬判断”带来的连锁错误。

真正让我下决心切换路线的是一次对比实验：在字典列相关性达到0.6左右、SNR只有15dB的情况下，我那个L1方案的恢复MSE已经接近原始信号能量的20%，而换用SBL重跑同样数据，误差降到10%以内。差距不是一点点，而是肉眼可见的。从那时起我就把贝叶斯方法作为默认选项，传统方法只作为实时性要求极高时的备选。

2. 稀疏贝叶斯学习SBL：分层先验与证据最大化的完整推导

2.1 把稀疏性写进先验：分层高斯模型

想用贝叶斯方法解决稀疏表示，第一步是要把“大部分系数为零”这个直觉转化成数学上可计算的形式。最简单的想法是给每个系数x_i单独设一个先验分布，比如拉普拉斯分布，它确实在零点有尖峰，能促进稀疏，但后验推断会涉及L1范数的非共轭性，计算起来比较麻烦。

SBL采用了一种更精巧的做法：分层高斯模型。第一层，给定每个系数的精度alpha_i，x_i服从零均值的高斯分布，方差就是1/alpha_i。注意这里的alpha_i是逐分量独立的，每个系数都有自己的收缩强度。第二层，再给这些alpha_i施加一个Gamma分布先验。两层合在一起，得到的边际分布是一族具有重尾性质的分布，在零点附近有很高的概率密度，同时允许少数系数取较大的值。这就非常贴合稀疏信号的真实特性——绝大多数系数压向零，少量系数可以较活跃。

这个设计的精妙之处在于“共轭性”。高斯似然配高斯先验，后验仍然是高斯分布，所有推断都有解析表达式。alpha_i的存在相当于为每个系数按了一个自动调节的“收缩旋钮”：如果某个系数被数据强烈支持，它对应的alpha_i会变小，方差变大，允许系数自由取值；如果某个系数不被数据支持，alpha_i会变得很大，方差趋近于零，把系数硬压到零。所谓“稀疏性从数据中自动学出来”，说的就是这个机制。

2.2 后验推断与超参数更新公式

模型设定好之后，剩下就是标准的贝叶斯推断流程。先给定超参数alpha和噪声精度beta，系数x的后验分布可以解析写出：

• 后验协方差：Sigma = (beta * A' * A + diag(alpha))^(-1)
• 后验均值：mu = beta * Sigma * A' * y

这里的mu就是我们最终对系数向量的估计，Sigma则提供了不确定性度量。但问题在于，alpha和beta本身也是未知的，怎么定？SBL的做法是最大化证据函数，也就是观测数据y在超参数下的边际似然 p(y | alpha, beta)。把x积分掉之后，这个边际似然有显式形式，可以用EM算法或者直接对alpha求导得到固定点迭代公式。

标准的更新公式如下：

• 令 gamma_i = 1 - alpha_i * Sigma_ii
• 更新 alpha_i = gamma_i / (mu_i² + eps)
• 更新噪声精度 beta = (N - sum(gamma)) / (||y - A*mu||² + eps)

这些公式初看有点抽象，但物理含义很清楚。gamma_i介于0和1之间，可以理解为“第i个系数被数据支撑的程度”。当某个系数后验均值mu_i很小时，alpha_i就会变大，把系数推向零；当mu_i显著时，alpha_i保持在一个较小的水平，保留这个系数。整个过程里没有人为指定的稀疏度，支撑集的确定是迭代过程中自然涌现的。

2.3 迭代收敛到哪里：证据回路的实际行为

理论公式归理论，实际跑起来观察到的行为很有意思。我第一次用MATLAB把EM循环写出来，发现大概前20次迭代里所有alpha_i都变化得很平缓，就像在“观望”，之后突然进入加速阶段——非支撑集的alpha_i快速增长，支撑集的alpha_i缓慢下降。再过几十次迭代，alpha_i的分布就呈现出明显的两极化：要么极大，要么极小。这时候mu在支撑集上的取值基本稳定。

另一个让我印象深刻的点是，负证据函数并不总是严格下降，偶尔会有一个微小的反弹，然后又继续下降。这其实是因为EM更新用的是固定点迭代而非完全的坐标下降，数值上只要下降趋势总体符合预期就不必担心。我后来在收敛判断里增加了“连续若干次迭代内alpha的最大变化小于阈值”的条件，避免被这种小波动干扰。

还有一点需要提醒：不要指望一次EM迭代就收敛。按照我通常的实验规模，N=30、M=64、稀疏度3到5的情况下，50到300次迭代内基本能稳定，但保守起见上限设在1000次。后面会专门讨论收敛判据的工程细节。

3. MATLAB R2018实现要点：矩阵化迭代与数值稳定处理

3.1 核心迭代代码结构

MATLAB R2018在语法层面跑SBL没有任何障碍，关键是要把迭代写成矩阵形式，避免逐元素循环。下面是我在实践里打磨过的核心循环骨架，去掉了项目里的私有封装，保留了最本质的步骤：

% y: N x 1 观测, A: N x M 字典, 通常 N << M % alpha: M x 1 系数精度, prec: 噪声精度(scalar) [N, M] = size(A); y = y(:); alpha = ones(M, 1) * 1e-2; prec = 1 / (0.1 * var(y) + eps); maxIter = 1000; tol = 1e-6; for iter = 1:maxIter Lambda = spdiags(alpha, 0, M, M); AtA = A' * A; Sigma = (prec * AtA + Lambda) \ eye(M); mu = prec * (Sigma * (A' * y)); gamma = max(0, 1 - alpha .* diag(Sigma)); alpha_new = gamma ./ (mu.^2 + 1e-12); prec_new = (N - sum(gamma)) / (norm(y - A * mu)^2 + 1e-12); delta = max(abs(alpha_new - alpha) ./ (abs(alpha) + eps)); alpha = alpha_new; prec = prec_new; if delta < tol break; end end

这里最重要的两个操作是Sigma的求解和mu的计算。我建议不要直接用inv(prec * AtA + Lambda)，而是用反斜杠运算符配合单位矩阵，让MATLAB底层选择合适的高斯消元策略，对小规模矩阵两者差距不算大，但代码健壮性好很多。

3.2 矩阵求逆的稳定替代方案

当M变大，比如到了几百甚至上千，直接构建M×M矩阵再求逆就变得很笨重，数值上还容易出现警告。这时一个经典手法是换用Woodbury恒等式，把求逆从M维空间降回N维空间：

Sigma = Lambda^(-1) - Lambda^(-1) * A' * (A * Lambda^(-1) * A' + (1/prec) * I)^(-1) * A * Lambda^(-1)

由于N远小于M，中间那个N×N矩阵的逆要便宜得多。mu也可以利用这个结构直接算：mu = prec * Sigma * A' * y，相当于先算A'*y，再通过逆矩阵作用一次，避免了显式构造完整Sigma。在我的实验里，M=256、N=32时，Woodbury版本比直接求逆快差不多一个数量级，而且数值稳定性肉眼可见地好。

另外，如果担心矩阵条件数过大，可以考虑对A的列做归一化处理。这一步放在预处理阶段，效果比在迭代里折腾数值更直接。归一化之后，alpha_i的尺度也会统一很多，后续设阈值时不用每个分量单独适配。

3.3 收敛判据设计的细节

工程实现里收敛判据比公式本身更影响使用体验。我的经验是不要只看单一指标。alpha的相对变化是一个不错的指标，但有些场景下alpha已经稳定而mu还在缓慢漂移，反过来也存在。更稳妥的办法是同时监控两项：alpha的最大相对变化和mu的最大相对变化，只要任意一项连续几次低于阈值，就认为收敛。

阈值也不要设得太激进。R2018的双精度下，1e-6到1e-8之间是合理区间。再小的话，迭代后期会陷入纯粹的数值震荡，白白浪费时间。另一个不太优雅但非常有效的办法是固定迭代上限，比如500或1000，保证最坏情况下程序也能退出。我通常会在循环外记录实际退出时的迭代次数，方便批量实验时后续分析。

如果要做多个独立实验的批量仿真，可以在循环外加一层for或者parfor。R2018的Parallel Computing Toolbox支持parfor，但需要确认并行池已开启。实测下来，并行开销在小规模矩阵上有时反而比串行慢，所以只有M和N都比较大时才建议开启。

4. 仿真实测：稀疏度、字典相关性和噪声水平下的真实表现

4.1 实验设定与评价指标

理论说得再好，最终要拿数据说话。我构建了一组标准实验，参数如下：观测维度N=30，字典原子数M=64，真值x的稀疏度K分别取3、5、8，非零位置随机，幅度随机取正负。字典A除了标准高斯随机矩阵之外，我额外构造了一组相关性增强的字典：先随机生成高斯矩阵，再对每列施加一个平滑核，使得相邻列的相关性升高到0.5到0.8之间。噪声按信噪比SNR从5dB到30dB分档添加。

评价指标我用三个：第一个是恢复信号与真值的均方误差MSE，归一化到信号能量；第二个是支撑集识别准确率，定义为零位置判零、非零位置判非零的逐元素正确率；第三个是成功恢复率，即支撑集与真支撑集完全一致的比例，这个指标非常苛刻，但最能反映实际工程里“找对地方”的能力。

为了保证结论可靠，每个参数组合我都跑了至少200次蒙特卡洛仿真，记录均值与方差。批处理就是用前面说的parfor实现的，R2018跑起来很顺畅，前提是先把A和y在parfor外定义好，否则worker之间反复传输大矩阵反而拖慢速度。

4.2 与OMP/L1方法的对比观察

结果有几个让我印象深刻的点。第一，在高斯随机字典、SNR=20dB、K=5的场景下，SBL、OMP、L1方法（我用几轮坐标下降实现的同规模L1求解）都能恢复到不错的状态，SBL的优势并不夸张，只比L1低零点几个dB。但一旦切到相关性增强的字典，差距就拉开了：OMP的成功恢复率跌到70%左右，经常把相邻列错认成支撑集；L1方法略好但恢复幅度明显有收缩偏差；SBL的成功率还在95%附近，恢复出来的幅度也更接近真值。

第二个有趣的现象是收敛路径。SBL的误差曲线在前期有一段相当长的“潜伏期”，大概几十次迭代内改善非常缓慢，然后突然进入快速下降阶段，最终收敛到一个比L1更低的水平。这个行为有点像模拟退火，前面是在探索正确的超参数区域，后面才集中优化。理解了这一点，就不会看到前期曲线平缓就误以为算法卡住而提前终止。

第三个关于计算成本的观察也很有价值。OMP在这种规模下毫秒级就结束，L1也就几十毫秒，而SBL通常需要一两秒。但考虑到很多场景的离线处理并不在意这一两秒，而准确率提升是实打实的，这个取舍我认为非常值得。如果对实时性要求极高，可以把SBL当作“离线训练”阶段：先在小规模样本上调好先验分布参数，再用OMP做在线快速恢复。

4.3 支撑集识别的失败模式

失败模式比成功数据更值得研究。我专门把所有仿真里恢复失败的样本挑出来看，发现典型的坏情况有这么几类。

一是低信噪比下出现“幽灵支撑”，也就是把噪声原子的alpha也压得比较小，系数mu非零，导致支撑集多出来几个假位置。这种情况在SNR低于10dB时尤其明显，本质上是因为数据信息量不足，后验分布本身就模糊，任何算法也没办法百分之百恢复真实支撑集。SBL的优势在于它能给出后验方差，此时方差通常很大，等于在主动提醒“这部分不可靠”。

二是字典相关性极强时出现支撑集“成团”现象。两个相邻原子高度相似，SBL会把两个都保留下来，各自分到一部分权重，而不是果断扔掉一个。我最初觉得这是个缺陷，后来发现这其实是概率推断的诚实反映——数据确实无法分辨这两个原子，保留两个比强行二选一在统计上更合理。实际使用中，如果必须得到“干净”的支撑集，可以再加一道后处理，比如对成团的原子做一次原子合并，或者用相关系数聚类后再选择代表原子。

三是初始化设置不当导致的“过压缩”。当alpha初始值取得过大，相当于先验强烈认为所有系数都是零，迭代就很难把真正有信息的原子拉回来，最后所有系数都被压成零。这个现象在做高稀疏度实验时特别容易出现，后面专门讲怎么规避。

5. 调参与避坑：从初始化到阈值选择的实战心得

5.1 超参数初始化的影响

对付“过压缩”最直接的手段是合理的alpha初始化。我试过多种取值：统一取0.01是万金油，在大多数场景下能正常工作，但在K=8的中等稀疏度下偶尔会收敛到局部最优。后来我改用一组更聪明的策略：在没有先验信息时，把alpha初始化为1 / (var(y) / M) 的数量级，这相当于先假设信号能量均匀分布在各原子间，再由迭代自己去浓缩到真正的支撑集上。实测下来，这种初始化在K=8时的成功恢复率提升了大约5个百分点。

噪声精度prec的初始化同样关键。最差的做法是随意给一个值。更好的做法是用拟合残差做粗估计：先算一个最小二乘解 y_pseudo = pinv(A) * y，残差r = y - A * y_pseudo，然后令 prec_init = 1 / (r' * r / N + eps)。这种方法在SNR中等以上时非常稳，基本不需要再手动调整。如果连残差都因为矩阵病态算不准，就退回到 prec_init = 1 / (0.1 * var(y) + eps)。

我还做过一个低配版的“多起点”策略：分别用alpha_init = 1e-4、1e-2、1三组初始化跑同一个实验，最后按证据函数大小选最优结果。这个方法在批量离线实验里非常有用，虽然计算量翻了三倍，但对局部最优问题的缓解是立竿见影的。对于一个强调稳定性的系统，这三倍开销我完全能接受。

5.2 停止条件和阈值选择

停止条件我在3.3节已经说过要双指标监控，这里补充一个具体可用的组合：alpha相对变化小于1e-6，或者mu的最大变化小于1e-8，两者只要连续满足三次就退出。配合maxIter=1000上限，几乎不会出现死循环。

支撑集的提取阈值是另一个容易踩坑的地方。我最初试着用绝对阈值判断alpha_i是否“大”，比如alpha_i > 100就算作零。但alpha的尺度严重依赖数据规模和字典归一化方式，绝对阈值根本不通用。后来改成相对准则：计算alpha数组的最大值maxAlpha，某个原子如果满足 alpha_i > 1e4 * min(alpha_j)，就认为它不在支撑集内。这个准则几乎不需要针对不同数据单独调参。

更实用的一招是“SBL辨识支撑 + LS校正幅度”。SBL给出的后验均值mu自带一定的收缩效果，因为分层先验对大系数也有少许压缩。支撑集确定后，把支撑列抽出来组成A_s，再直接用普通最小二乘 x_refined = A_s \ y 重新估计一次幅度，可以让最终恢复精度再提升一截。我几乎所有仿真最后都用这个两步走流程，效果普遍比直接用mu好。

5.3 我在实践中踩过的几个坑

这几个坑都是真实花过时间才解决的问题，写出来给后来人省点力气。

第一个坑是最开始直接用了inv(prec * AtA + Lambda)，几十维矩阵就偶尔弹出“矩阵接近奇异或缩放错误”的警告。换成反斜杠后警告消失，但当时还没意识到性能问题，直到M扩到256才反应过来，改用Woodbury。现在我的默认实现是Woodbury路径，代码也不复杂，建议直接参考这个做法，一步到位。

第二个坑是忘了对字典列做归一化。有一段时间我用同一套调参逻辑在不同预处理的字典上跑，结果支撑集识别率忽高忽低，后来才定位到是列模长不一致导致alpha的阈值判断失效。归一化之后，所有原子在同一个能量尺度下参与竞争，问题迎刃而解。

第三个坑和低噪声极限有关。当SNR特别高，比如接近无噪声时，残差项norm(y - A * mu)会变得非常小，prec_new的分母接近零，更新出来的噪声精度数值极不稳定，反过来又影响Sigma的计算。我处理的办法是给分母加一个eps级的小常数，同时增加一个分支：如果当前相对残差已经低于预设下限，比如1e-10，直接认为收敛退出，不做这种无意义的超参数更新。

第四个坑发生在批量仿真里。我习惯把多组实验的结果矩阵预先分配好，结果发现某个循环里忘记在每次实验重置随机数种子，导致所有实验结果高度相关。R2018的随机数全局流是自动推进的，但如果误用了相同的种子，复现出来的“不同”实验其实是同一组随机样本的重复，统计结论当然失真。这个错误很隐蔽，因为我是在分析结果分布时才发现的。提醒大家批量实验前一定要明确是否固定种子，并且固定后要主动混合不同种子的结果。

最后一个经验是关于“SBL是不是越快越好”的心态。我见过不少朋友上来就追求把迭代次数压到最少，恨不得20次内出结果。但实际上SBL的不确定性学习天然需要一定迭代轮数来传导信息，过于激进地截断往往让alpha还没分化就被误判为已收敛。与其牺牲准确性，不如老老实实保留几百次迭代上限，把时间花在更值得优化的数据编码和特征工程上。

现在回头看，用贝叶斯方法重写整个稀疏表示学习流程，最大的收获不是那几个百分点的精度提升，而是开始习惯用分布而不是点估计去思考问题。alpha的演化过程本身就是一副“数据如何逐步确认哪些原子可信”的生动图景。在后续的项目里，我甚至直接把SBL当时的二值化支持信息作为特征做了可解释性分析，效果出奇地好。MATLAB R2018在这个流程里表现稳定，矩阵化迭代和Woodbury配合起来非常顺手，这套代码我就留在自己的工具箱里当基线了。