简介这是一份稀疏贝叶斯学习SBL算法的 MATLAB 代码实现包面向从事信号处理、模式识别等高维数据建模的研究者与工程开发者。包内实现了 TSBL、TMSBL、MSBL、MFOCUSS 等多种稀疏贝叶斯算法并附带多个 demo 演示脚本便于对照论文理解算法细节与参数作用。资源共 14 个文件以 11 个 .m 源码脚本为主配合 2 个 PDF 算法使用说明和 1 个 ReadMe 文本压缩包整体仅 462KB结构紧凑清晰适合快速阅读与二次开发。已有 1509 人学习。通过运行 demo_identicalVector、demo_time_varying、demo_fig6_SNR10 等示例可直观观察不同算法在静态、时变及不同信噪比场景下的重构效果进而迁移到自身任务中完成稀疏系数求解与特征选择。1. 稀疏贝叶斯学习代码复现为什么说它比 Lasso 更值得自己动手写做特征选择或稀疏信号重构的人多半在 Lasso 的网格搜索里吃过亏λ 调小了解里全是非零项λ 调大了真信号也被一起压没了。同样是做稀疏化稀疏贝叶斯学习落地成代码时走的不是这条反复试参的路线它用层级先验把每个权重的方差变成可学习参数在迭代中自动把不相关的项“关掉”。这篇笔记把 SBL 从公式推导翻译成可运行的 Python 代码并给出模拟数据验证和五条实战避坑记录适合被 L1 调参折磨、想复现相关向量机RVM或给压缩感知、DOA 估计写基线算法的工程师。读完你至少能独立写出一份 100 行以内的核心实现并知道怎么判断它到底有没有跑对。2. 把 SBL 公式翻译成代码先看懂稀疏贝叶斯在优化什么2.1 层级先验与自动相关性确定SBL 的稀疏性不是靠惩罚项线性的观测模型写成 $t \Phi w \varepsilon$其中 $t$ 是 $N$ 维观测向量$\Phi$ 是 $N \times M$ 的字典矩阵$w$ 是待求的 $M$ 维权重噪声 $\varepsilon \sim \mathcal{N}(0, \sigma^2 I)$。Lasso 的做法是在目标函数里加 $\lambda |w|_1$用凸优化的方式把一部分权重压到零。SBL 换了一条路对每个权重单独设置一个零均值高斯先验$$p(w_i \mid \alpha_i) \mathcal{N}(0, \alpha_i^{-1}),$$这里的 $\alpha_i$ 是权重精度的倒数它自己又服从 Gamma 先验。把 $\alpha_i$ 边缘化之后每个 $w_i$ 的边际先验会变成 Student-t 分布。这个分布的特点是尖峰厚尾中间非常集中尾部却比高斯厚。落在尖峰附近的权重被推向零而真正有信号的权重因为尾部够厚不会被过度压缩。这套机制在文献里叫自动相关性确定ARD也是相关向量机的核心。与 Lasso 的 $\ell_1$ 凸代理相比SBL 的求解目标是非凸的但在字典列高度相关、真实稀疏度未知的场景下它通常能给出更接近 $\ell_0$ 行为的解这也是它常被拿来当压缩感知和阵列信号处理基线的直接原因。代码实现时要注意一个容易被忽略的点$\alpha_i$ 是每个维度独立估计的。这个“独立”是 ARD 稀疏性的来源也是量纲问题的来源。特征列没标准化时数值范围大的列天然更容易被选中后面第 4 章的避坑清单里第一条就是这个。2.2 从后验到迭代EM 推导与迭代权重更新有了先验框架接下来就是标准的贝叶斯推断。记噪声精度 $\beta \sigma^{-2}$权重先验精度矩阵 $A \mathrm{diag}(\alpha_1, \dots, \alpha_M)$。给定观测 $t$权重的后验分布仍是高斯$$\Sigma (\beta \Phi^T \Phi A)^{-1}, \quad \mu \beta \Sigma \Phi^T t,$$其中 $\Sigma$ 是后验协方差$\mu$ 是后验均值。SBL 的训练目标不是直接最大化后验而是最大化边缘似然 $\ln p(t \mid \alpha, \beta)$也就是把权重 $w$ 积分掉之后的证据。这个目标没法闭式求解常见做法是 EM 迭代把 $w$ 当作隐变量交替更新后验参数和超参数 $\alpha, \beta$。一轮 EM 更新可以写成四个式子$$\Sigma (\beta \Phi^T \Phi A)^{-1}, \quad \mu \beta \Sigma \Phi^T t,$$$$\gamma_i 1 - \alpha_i \Sigma_{ii}, \quad \alpha_i^{\mathrm{new}} \frac{\gamma_i}{\mu_i^2},$$$$\beta^{\mathrm{new}} \frac{N - \sum_i \gamma_i}{|t - \Phi \mu|_2^2}.$$这里的 $\gamma_i$ 有明确的直观含义它表示第 $i$ 个基函数被数据“解释”的程度取值在 0 到 1 之间。$\sum_i \gamma_i$ 可以当作有效参数个数。当某个 $\mu_i$ 接近零时$\alpha_i$ 会趋向无穷大对应的 $A_{ii}$ 变得极大后续迭代中该维度的 $\Sigma_{ii}$ 和 $\mu_i$ 进一步被压到零这就是 SBL 能把不相关特征“关掉”的微观机制。实际代码里需要给 $\mu_i^2$ 加一个极小值下界防止除零。收敛判据建议监控对数边缘似然的变化量而不是只看权重。边缘似然有解析式$$\ln p(t \mid \alpha, \beta) -\frac{1}{2} \left[ N \ln 2\pi \ln |C| t^T C^{-1} t \right], \quad C \beta^{-1} I \Phi A^{-1} \Phi^T.$$实现的时候用np.linalg.slogdet算对数行列式避免直接算 $C$ 的行列式导致数值溢出。这个式子不求绝对准确但用来判断迭代是否有进展已经足够。2.3 为什么这里不用 MCMCEM 和坐标上升才是代码主角第一次接触 SBL 的人很容易问既然有完整的贝叶斯后验为什么不直接用 MCMC 采样真实项目里很少这么做。SBL 要解的是 Type-II 最大似然也就是把权重积分掉之后只对超参数做点估计这个目标函数本身是确定的不需要采样。MCMC 在高维稀疏问题里收敛诊断复杂跑一次要几万步而 EM 每一步都是解析更新利用的是高斯-伽马共轭结构几十轮就能到平稳点。另一种常见路径是用坐标上升固定 $\beta$ 逐维更新 $\alpha$在字典维度很高时比整体 EM 更稳。梯度下降和 Adam 也能优化边缘似然但 SBL 目标非凸无结构梯度法容易停在差的局部点这一点第 4 章会单独讲。初学者直接走 EM 路线是最省事的代码量小行为也最容易对齐论文里的公式。3. 用 Python 把稀疏贝叶斯学习跑起来最小实现与调参笔记3.1 从零实现 SBL 核心迭代30 行搞定 EM 更新下面是一个最小可跑的稀疏贝叶斯线性回归实现代码只有几十行核心就是上一章的四组更新公式。先不要加花哨的加速和剪枝把基线跑通再说。import numpy as np class SparseBayesianLR: 稀疏贝叶斯学习SBL最小实现线性模型 EM 更新 观测模型: t Phi w noise def __init__(self, max_iter500, tol1e-8, alpha_init1.0, noise_initNone): self.max_iter max_iter self.tol tol self.alpha_init alpha_init self.noise_init noise_init def _log_marginal_likelihood(self, Phi, t, alpha, beta): 计算 ln p(t|alpha,beta)用于收敛判断 N Phi.shape[0] C (1.0 / beta) * np.eye(N) Phi np.linalg.inv(np.diag(alpha)) Phi.T sign, logdet np.linalg.slogdet(C) return -0.5 * (N * np.log(2.0 * np.pi) logdet t np.linalg.solve(C, t)) def fit(self, Phi, t): N, M Phi.shape alpha np.full(M, self.alpha_init, dtypefloat) if self.noise_init is None: # 噪声方差初值取观测方差的 10%是一个不容易翻车的起点 beta 1.0 / (np.var(t) * 0.1) else: beta 1.0 / self.noise_init L_prev -np.inf for it in range(self.max_iter): # E 步计算后验均值与协方差 A np.diag(alpha) Sigma np.linalg.inv(beta * Phi.T Phi A) mu beta * Sigma Phi.T t # M 步更新权重精度和噪声精度 gamma 1.0 - alpha * np.diag(Sigma) alpha gamma / np.maximum(mu ** 2, 1e-12) beta (N - gamma.sum()) / np.sum((t - Phi mu) ** 2) L self._log_marginal_likelihood(Phi, t, alpha, beta) if np.abs(L - L_prev) self.tol: break L_prev L self.mu_ mu self.Sigma_ Sigma self.alpha_ alpha self.beta_ beta return self def predict(self, Phi): 用后验均值做预测 return Phi self.mu_这段代码里有几个参数值得说明。alpha_init1.0表示所有权重先验精度从 1 起步对应先验方差为 1这个值在输入特征做了标准化之后基本不用改tol1e-8是边缘似然变化的绝对阈值实际数据噪声大时没必要这么严调到1e-4或1e-5能省掉不少空转noise_init默认取观测方差的 10%这是经验上比较稳的起点。M 步里np.maximum(mu ** 2, 1e-12)是为了防止某个权重被压到极小后除以零这个下界太小会导致 $\alpha$ 爆炸太大会把稀疏性吃掉1e-12 对普通浮点精度足够安全。每次迭代都要重新计算Phi.T Phi这个矩阵乘法在 $N \times M$ 较大时是主要开销后面优化可以缓存结果但第一次实现不建议加缓存先把逻辑跑对。3.2 用模拟数据验证稀疏恢复效果示例代码与判据代码写完第一件事不是上真实数据而是构造一个稀疏系数已知的模拟问题验证算法能不能找回真实支撑集。这样可以确定问题是出在实现还是出在数据。示例代码讲解的验证路径如下随机生成 120 样本、500 维字典真实稀疏系数只有 10 个非零位置加少量噪声后跑 SBL。rng np.random.default_rng(42) N, M 120, 500 # 样本数远小于字典列数属于欠定稀疏问题 # 随机字典 列归一化避免量纲干扰 Phi rng.standard_normal((N, M)) Phi / np.linalg.norm(Phi, axis0, keepdimsTrue) # 构造稀疏真实权重10 个非零位置 w_true np.zeros(M) support rng.choice(M, size10, replaceFalse) w_true[support] rng.standard_normal(size10) * 3.0 # 观测加高斯噪声信噪比约 30dB t Phi w_true rng.standard_normal(N) * 0.1 # 训练 SBL model SparseBayesianLR(max_iter500, tol1e-8) model.fit(Phi, t) # 提取估计支撑集取 |mu| 超过最大值的 1% 的位置 threshold 0.01 * np.max(np.abs(model.mu_)) est_support np.where(np.abs(model.mu_) threshold)[0] true_set set(support.tolist()) est_set set(est_support.tolist()) overlap len(true_set est_set) / len(true_set) print(支撑集重合率:, overlap) print(恢复误差:, np.linalg.norm(model.mu_ - w_true) / np.linalg.norm(w_true))这个验证里有两个细节容易翻车。第一个是字典列归一化如果不做某些列二范数天然大SBL 会偏好它们支撑集重合率会明显下降第二个是支撑集提取阈值用绝对值固定阈值会被量纲绑架相对阈值更稳通常取最大估计值的 1% 到 5%。判据方面支撑集重合率达到 90% 以上、恢复误差在 0.1 以内说明核心实现没有大问题。真实场景里如果重合率长期低于 70%不要急着调超参数先回去检查字典的条件数和噪声初值。3.3 三个必调参数收敛容差、迭代上限、噪声初值SBL 的实际调参负担比 Lasso 小很多但有几个参数值得认真对待。我整理成一张参数表方便在真实项目里直接对照。参数典型取值影响调整建议alpha_init1.0太大起点稀疏性过强太小起点接近岭回归特征标准化后统一用 1.0多起点验证时取 0.1、1.0、10.0tol1e-6 ~ 1e-8太松会提前停太紧会白白迭代数据噪声大时放宽到 1e-4干净模拟数据用 1e-8noise_initvar(t) * 0.1最关键给太大会把信号当噪声给太小会不稀疏先按 0.1 倍方差起步非零项过多就往下调到 0.01 倍tol是边缘似然的绝对变化阈值不是权重变化。真实项目里常见误区是只盯mu的变化结果边缘似然还没平稳就停了。noise_init是这三个里影响最大的一个它直接决定 EM 迭代的起点在边缘似然曲面上的位置。把噪声初值设成观测方差的 0.1 倍是通俗做法但如果观测数据信噪比本身很低这个起点会让算法把噪声当成真实信号结构表现为支撑集过大这时先把噪声初值调小到 0.01 倍方差再跑。这三个参数在代码里实现非常简单但调起来需要配合边缘似然曲线一起看单看某个指标很容易被误导。4. 稀疏贝叶斯学习代码复现的五个常见坑现象、原因与解决4.1 特征没做标准化α 先验把大数值列全“选中”了现象数据里某些特征取值范围在 0 到 1另一些在 0 到 10000训练完成后非零权重几乎全落在大数值列上。原因所有 $\alpha_i$ 初值相同先验方差一致但大数值列在 $\Phi^T\Phi$ 里的贡献远大于小数值列。EM 更新时这些列的 $\mu_i$ 天然偏大$\gamma_i$ 也偏大算法等于给它们开了一条快速通道。解决建模前对 $\Phi$ 做 z-score 标准化或列归一化求解完成后再把估计权重按列范数还原。SBL 对特征尺度比 Lasso 更敏感因为稀疏性完全靠独立的 $\alpha_i$ 驱动量纲不一致直接扭曲 ARD 的公平性。4.2 噪声方差初值给太小稀疏解变成稠密解现象迭代收敛后非零项有上百个看起来像 Lasso 在 $\lambda$ 很小时的解完全没有稀疏性。原因$\beta$ 初值太大相当于噪声方差初值太小似然项在目标里占比过高后验均值几乎贴着最小二乘解走$\alpha_i$ 没有机会把弱权重压到零。解决把noise_init调回np.var(t) * 0.1起步如果确实需要更小的初值每次只降一个量级比如先 0.1 倍方差再 0.01 倍。同时观察边缘似然曲线如果迭代初期就出现大幅上升后停滞多半是噪声初值太乐观。4.3 高相关字典下迭代震荡边缘似然曲线来回跳现象对数边缘似然不单调每轮迭代上下波动$\mu$ 也跟着震荡收敛判据长时间不触发。原因字典列高度相关时$\Phi^T\Phi$ 的条件数很大$\Sigma (\beta\Phi^T\Phi A)^{-1}$ 的病态程度加剧EM 更新在 $\alpha$ 的坐标方向上出现振荡。解决先做列归一化并用np.linalg.cond(Phi.T Phi)检查条件数给 $\Sigma$ 的对角项加微小 jitter例如 1e-6能显著改善数值稳定性如果震荡仍然存在把 EM 改成固定 $\beta$、只对 $\alpha$ 做坐标上升的混合策略收敛行为更平稳。4.4 只看权重变化就停边缘似然还没到位现象$\mu$ 的变化已经小于 1e-6提前退出迭代但在验证集上预测误差仍然偏大。原因权重变化小不代表超参数已经收敛。EM 后期 $\mu$ 微调时$\alpha$ 和 $\beta$ 可能还在缓慢漂移尤其当 $\gamma$ 接近 0 或 1 时$\alpha$ 的更新步长会被放大。解决把主收敛条件改成边缘似然变化量$\mu$ 的变化只当辅助参考。边缘似然需要自己实现用第 2 章的公式即可数据量大时可以在每 10 轮算一次边缘似然降低日志行列式的计算次数。4.5 EM 收敛到局部极值支撑集只有一半现象模拟数据上支撑集重合率不到 50%但 MSE 不算离谱说明算法停在了边缘似然的某个局部极值。原因$\ln p(t \mid \alpha, \beta)$ 是非凸函数EM 的起点决定了走向。单一初值很容易停在一个“半稀疏”的鞍点附近。解决做多起点验证alpha_init分别取 0.1、1.0、10.0 跑三遍选边缘似然最大的结果还可以在迭代初期固定 $\beta$只让 $\alpha$ 先跑 20 轮给它一个更合理的搜索方向后再放开 $\beta$ 一起更新。后者在真实项目里往往比多起点更有效因为它直接避免了“噪声和信号互相竞争”的早期混乱。5. 用支持集恢复率给 SBL 代码做“出厂测试”一个比 MSE 更严格的验证技巧代码能跑通之后最后一步值得做的验证是扫描信噪比画支撑集恢复率曲线。MSE 低了不代表稀疏结构找对了在很多工程场景里你真正关心的是“哪几个传感器通道、哪几个特征在起作用”而不是权重数值有多接近。支持集恢复率直接回答这个问题也会暴露出 MSE 看不出的噪声初值问题和字典相关度问题。做法很简单保持真实支撑集不变把噪声标准差从 0.01 扫到 1.0在每个信噪比下重新生成观测、重新训练 SBL统计支撑集重合率。我一般用这样的小循环snr_levels [0.01, 0.05, 0.1, 0.2, 0.5, 1.0] recovery_rates [] for noise_std in snr_levels: rates [] for seed in range(10): rng np.random.default_rng(seed) t Phi w_true rng.standard_normal(N) * noise_std model SparseBayesianLR(max_iter500, tol1e-6) model.fit(Phi, t) est np.where(np.abs(model.mu_) 0.01 * np.max(np.abs(model.mu_)))[0] rates.append(len(set(est) set(support)) / len(support)) recovery_rates.append(np.mean(rates)) print(fnoise_std{noise_std:.2f}, recovery_rate{np.mean(rates):.2f})理想的 SBL 实现会在噪声标准差小于 0.1 时保持 90% 以上的恢复率随着噪声增大缓慢下降。如果曲线在低噪声段就掉到 70% 以下基本可以断定实现或初值设置有问题而不是 SBL 方法本身不行。这个测试也适合拿来对比不同字典归一化方式我自己的项目里列归一化加噪声初值 0.1 倍方差是组合里最稳的一组。最后留一个个人习惯每次跑真实数据之前我都先花十分钟跑一遍这个支持集恢复率测试。以前有次在结构健康监测数据上MSE 看起来很低实际恢复的支撑集和物理机理完全对不上折腾了一周才发现是噪声初值设得过于乐观。那之后这条流程就没断过。希望这个验证习惯能帮你少走一段弯路把时间花在真正该调整的地方。本文还有配套的精品资源点击获取