简介压缩感知稀疏贝叶斯算法实现包面向信号处理、压缩感知与贝叶斯推断方向的研究者或工程师聚焦SBL、TSBL、TMSBL从有限测量中恢复稀疏动态信号附带作者实测可用的完整Matlab代码。压缩包共15个文件主体为11个m脚本覆盖MSBL、MFOCUSS、TSBL、TMSBL等核心算法与多组演示脚本另有2份PDF快速上手指南、1份txt说明和1份docx配套文档总大小479KB轻量易用适合课程设计、科研复现与算法对比。已有929人学习下载。脚本按场景组织包含不同信噪比、多尺度时间序列、平稳向量与时变信号等演示用例可还原多种重构效果结合辅助函数便于对照论文验证、调整参数并迁移至图像处理、医学成像、通信信号检测等实际任务。使用者还可借助说明文档与快速上手指南快速理解核心流程并根据自身数据修改输入矩阵和稀疏约束。1. 压缩感知里的稀疏贝叶斯算法SBL 能解决什么问题适合谁做压缩感知恢复时我第一次把 L1、OMP 和稀疏贝叶斯算法SBL放在同一批低测量数数据上对比SBL 的表现让我直接放弃了继续调正则参数。它不需要预知稀疏度 K也不像 OMP 那样一次只挑一个原子而是把“哪些位置非零”当成参数估计出来。标题里的三个算法是一家人SBL 处理单个测量向量TSBL 把时间相关性写进先验TMSBL 同时利用多测量向量和时间结构。下面这篇文章按我复现并实测跑通的路子来写从原理、代码到参数和坑照做就能用。适合刚接触压缩感知恢复、或在 L1 和贝叶斯方法之间犹豫的从业者。2. SBL 的信号模型与 EM 迭代先验怎么选更新公式怎么落到代码稀疏贝叶斯学习的观测模型是 y Φx n其中 Φ 是 M×N 的测量矩阵M 远小于 Nx 本身只有 K 个非零元。贝叶斯路线不直接对 x 做 MAP 估计而是给每个 x_i 配一个独立的先验方差 γ_ix_i ~ N(0, γ_i)。γ_i 不是手工给的是靠观测数据用 Type-II 极大似然估出来的。估计完成后大部分 γ_i 被压到接近零对应位置的后验均值也趋近零这就是“稀疏”的来源。和 L1 类方法对比差别很实际L1 或近端类方法是在固定惩罚系数 λ 下做凸优化λ 要靠交叉验证挑SBL 的所有未知方差都交给数据不需要调稀疏惩罚强度。代价是每次迭代要做一个 N×N 的协方差求逆计算量比 L1 高一个量级。后面会讲到怎么用 Woodbury 恒等式把这个瓶颈消掉。2.1 两层先验高斯尺度混合怎么把“稀疏”写进贝叶斯框架把先验写成 x_i sqrt(γ_i) · z_iz_i ~ N(0,1)γ_i 是尺度因子。从宏观上看x_i 服从一个由 γ_i 调制的高斯分布边缘分布是重尾的这正是稀疏信号需要的特性大部分 x_i 接近于零少数 x_i 幅度大。实际实现时γ 初始化为全 1 也能收敛我试过初始化为随机正数最后支撑集位置基本一致稳定性比 OMP 好。OMP 对测量矩阵列之间的相关性很敏感两列高度相关时它可能选错原子SBL 的支撑判断依赖后验方差天然带有“这些列都不确定”的加权抖动小得多。噪声也要建模n ~ N(0, σ²I)。σ² 是噪声方差同样用数据估计。信噪比高时 σ² 的更新会往零跑如果不加下界求逆矩阵会变得病态建议在代码里把 σ² 截断到 1e-8 这类小值。2.2 EM 迭代γ 和 σ² 的更新公式与边界给定 γ 和 σ²x 的后验还是高斯分布两个关键量是Σ ( ΦᵀΦ / σ² Γ⁻¹ )⁻¹μ Σ Φᵀ y / σ²其中 Γ diag(γ)。EM 的 E 步就是用当前参数计算 Σ 和 μM 步更新超参数γ_i ← μ_i² Σ_iiσ² ← ( ||y - Φμ||² trace(ΣΦᵀΦ) ) / M注意 γ_i 的更新同时包含后验均值的平方和后验方差。这意味着即使某个位置的后验均值比较大如果它的方差也大也就是这个位置不确定下一轮也不会被激进放大这是它比简单阈值法稳的原因。参数上需要关注两处γ_i 建议加一个下界 1e-8防止某次迭代把 γ_i 更新到 0 后永远锁死σ² 更新公式里的 trace(ΣΦᵀΦ) 可以换写成 trace((ΦΣ)Φᵀ)把 N×N 矩阵乘法降成两个 M×N 乘 M×M 的运算。2.3 最小实现80 行的 numpy SBL 求解器下面这个实现是我在测试里一直用的版本去掉注释不到 40 行把原理完整走了一遍import numpy as np def sbl(y, Phi, max_iter200, tol1e-3, sigma2_initNone, gamma_initNone): SBL 求解 y Phi x n返回后验均值 mu超参数 gamma 和噪声方差 sigma2。 M, N Phi.shape # 超参数初始化经验值 gamma np.ones(N) if gamma_init is None else gamma_init sigma2 (np.var(y) / 10.0) if sigma2_init is None else sigma2_init PhiT Phi.T for _ in range(max_iter): # 用 1/gamma 直接构造对角阵避免对 diag(gamma) 求逆 Gamma_inv np.diag(1.0 / gamma) Sigma np.linalg.inv(PhiT Phi / sigma2 Gamma_inv) mu Sigma PhiT y / sigma2 # 更新 gamma_i mu_i^2 Sigma_ii gamma_new mu**2 np.diag(Sigma) # 更新噪声方差 resid y - Phi mu PhiSigma Phi Sigma sigma2_new (resid resid np.trace(PhiSigma PhiT)) / M # 用 gamma 的相对变化判断收敛 rel_change np.max(np.abs(gamma_new - gamma)) / (np.max(gamma) 1e-12) gamma, sigma2 gamma_new, sigma2_new if rel_change tol: break return mu, gamma, sigma2代码里有两个地方值得说明。第一初始化 σ² 用观测方差 var(y) 的十分之一这是一个经验值太大会让前几十轮迭代像在做最小二乘收敛慢太小会过早压制支撑集。第二相对变化阈值 to1 设置为 1e-3我测试中再往下调到 1e-5 对恢复结果影响很小只会拖慢速度所以工程上 1e-3 够用。提示这个版本每次迭代都做一次 N×N 的求逆N ≤ 512 时没问题N 到几千后要用 Woodbury 恒等式把求逆降到 M×M最后一章给代码。3. 从单快照到多快照TSBL 和 TMSBL 的时间相关性是怎么建模的上一章的 SBL 一次只吃一个 y。雷达、通信、脑电里更多时候拿回来的是 L 个快照它们共享同一个稀疏支撑集只是幅度不同。如果每个快照单独跑 SBL恢复出来的支撑集会抖动把 L 路拼成 Y ΦX N对 X 的行做联合稀疏约束恢复概率会明显改善。TSBL 和 TMSBL 都在这个框架里。先说清楚一个事情TSBL 在文献里的缩写没有唯一严格含义我按工程口径理解——当 L 个快照之间有明确时间先后关系且幅值随时间连续变化时把相邻快照的相关性写进先验。TMSBL 则在此基础上更进一步多测量向量和时间相关性同时建模相关矩阵 B 不从外部指定而是由算法自己学出来。很多场景下先跑通“共享支撑的多快照 SBL”再去升级到 TMSBL是最稳的路径。3.1 从单任务到多任务TSBL 的时域相关性怎么建模多测量向量模型写为 Y ΦX NX 是 N×L行向量表示某个原子在所有 L 个快照里的幅度。联合稀疏的含义是一行要么整行接近于零要么整行非零但非零行内部各快照的幅度可以不同。最朴素的多快照扩展是共享一套 γ每个快照单独算后验但 γ 的更新把 L 路的信息平均起来而不是每个快照各估各的。这相当于假设快照之间相互独立但支撑集完全一致。我一般把这个版本叫共享支撑 MMV-SBL它其实是 TSBL 的基础版代码量很小def mmv_sbl(Y, Phi, max_iter100, tol1e-3): 共享支撑多快照 SBL。 Y: M×LPhi: M×NX 的每一行共享同一个 γ。 M, N Phi.shape L Y.shape[1] gamma np.ones(N) sigma2 np.var(Y) / 10.0 PhiT Phi.T for _ in range(max_iter): Gamma_inv np.diag(1.0 / gamma) Sigma np.linalg.inv(PhiT Phi / sigma2 Gamma_inv) U Sigma PhiT Y / sigma2 # N×L后验均值 # 每个原子在所有快照上的平均能量 gamma_new np.mean(U**2, axis1) np.diag(Sigma) R Y - Phi U sigma2_new (np.sum(R**2) L * np.trace((Phi Sigma) PhiT)) / (M * L) rel_change np.max(np.abs(gamma_new - gamma)) / (np.max(gamma) 1e-12) gamma, sigma2 gamma_new, sigma2_new if rel_change tol: break return U, gamma, sigma2这段代码里U 的每一列对应一个快照的恢复结果γ 的更新公式把 L 列的后验均值平方取平均。注意 Σ_ii 对所有快照共用因为所有快照共享同一套测量矩阵和 γ。这意味着如果一个原子在单快照上不稳定多快照平均后它的 γ 会被压得更准这是多快照方案的核心收益。TSBL 在共享支撑的基础上加入时间相关时常见的实现是在 γ 更新时不做简单平均而是用相邻快照的幅值平滑项或者给行先验加上一个 L×L 的相关矩阵 B。工程上如果相邻快照的幅值变化是连续的这一项提升很明显如果快照本身是独立采样的加了反而过拟合。3.2 TMSBL把多测量向量与时间结构揉在一起学习跨任务相关矩阵TMSBL 的做法是把每一行的先验从 γ_j · I 升级为 γ_j · BB 是 L×L 的相关矩阵。γ_j 控制第 j 行整体是否活跃B 控制这个活跃行在 L 个快照之间的相关模式。当 B 退化为单位阵时TMSBL 就退化成上一节的共享支撑版本当 B 的非对角元被学到非零值时表示相邻快照存在时间相关。B 的更新在 EM 框架里是一个加权外积累加对每一个被判定为活跃的原子用它的恢复行向量 u_j 做外积并按 γ_j 加权。我在实现里不会每一轮都更新 B而是先固定 B I 跑前 50 轮或直到支撑集大致稳定然后每 5 轮再更新一次 B并在对角上加 1e-6 的扰动避免病态def update_B(U, gamma, eps1e-6): TMSBL 中相关矩阵 B 的工程近似更新。 U: N×L 后验均值gamma: 当前稀疏超参数。 L U.shape[1] B np.zeros((L, L)) for j in range(gamma.size): # 只让活跃原子参与避免零行把相关矩阵冲淡 if gamma[j] 1e-4 * np.max(gamma): uj U[j].reshape(-1, 1) B (uj uj.T) / gamma[j] B / gamma.size # 对角加载保证正定性 eig_min np.linalg.eigvalsh(B)[0] if eig_min eps: B eps * np.eye(L) return B这个版本是工程近似没有把后验方差的块修正项加进去精度上会比论文原版略差但在我测过的场景里支撑集恢复能力已经很接近完整实现了。如果你要做严格对比再去补 Σ_j 对应的块修正。3.3 选型判断什么场景用 SBL / TSBL / TMSBL场景推荐算法选择依据单测量向量静态稀疏K 未知SBL不需要稀疏度先验不用调 λ多测量向量各快照独立共享支撑MMV-SBL共享 γ快照间没有时间相关B 无意义多测量向量快照间幅值连续变化TMSBLB 矩阵自动学习时间相关单快照流式恢复支撑集随时间缓慢漂移TSBL相邻帧约束用前一帧结果初始化 γ加时域平滑我个人的判断标准是先算一下快照矩阵 Y 的行间自相关。如果相邻快照的相关系数低于 0.3就不要用 TMSBL硬加 B 只会引入偏差如果相关系数长期在 0.7 以上TMSBL 的收益非常显著。4. 完整测试流程合成稀疏信号、压缩观测与三算法对比写算法实现只是第一步真正判断能不能用要靠一份可复现的测试流程。我通常固定随机种子先生成已知支撑集和幅值的稀疏信号再用随机高斯测量矩阵做压缩观测最后用不同算法恢复看误差和支撑集重合度。4.1 合成数据与压缩观测保证实验可复现的生成规则生成测试数据的关键是让稀疏度 K、信号长度 N、观测数 M 三者关系明确。M/K 的比值决定问题的病态程度一般 M/K 在 8 到 15 之间才有稳定的恢复效果。我常用 N256、M80、K10这一组参数下 SBL 类算法能稳定恢复但又不是简单到看不出差异。rng np.random.default_rng(42) N, M, K 256, 80, 10 # 随机选支撑集幅值用标准正态 supp_true rng.choice(N, sizeK, replaceFalse) x_true np.zeros(N) x_true[supp_true] rng.standard_normal(K) # 高斯随机测量矩阵列归一化到能量 1 Phi rng.standard_normal((M, N)) Phi / np.linalg.norm(Phi, axis0) # 加噪噪声方差 0.01 对应较高 SNR y Phi x_true 0.01 * rng.standard_normal(M)测量矩阵列归一化这步很关键。如果不归一化Φ 各列能量差异会让 γ 初始化时的量级判断失效SBL 的前几十轮迭代会浪费在调整 γ 上。噪声方差 0.01 是我测试里比较舒服的水平能看出算法差异又不会让所有算法都失败。调低到 0.1 时SBL 的优势会被噪声盖住这也是注意点贝叶斯方法的噪声方差估计能力在中等信噪比下才发挥得出来。4.2 评估指标NMSE、支撑集重合度与成功率评估一个恢复结果不能只看一张图我常用三个指标。NMSE 归一化均方误差反映幅值整体拟合程度支撑集重合度反映“选对了哪些原子”成功率则用来统计多次实验的稳定性。def nmse(x_hat, x_true): return np.linalg.norm(x_hat - x_true)**2 / np.linalg.norm(x_true)**2 def support_overlap(x_hat, x_true, rtol1e-2): # 用最大幅度的百分之一作为非零判决阈值 thr rtol * np.max(np.abs(x_hat)) supp_hat np.where(np.abs(x_hat) thr)[0] supp_true np.where(np.abs(x_true) 0)[0] return len(set(supp_hat) set(supp_true)) / len(supp_true) # 单次实验判成功NMSE 0.1 且支撑集重合 1 def is_success(x_hat, x_true): return nmse(x_hat, x_true) 0.1 and support_overlap(x_hat, x_true) 1.0注意支撑判决阈值不要写死成一个绝对数值。不同应用的信号尺度差异很大用相对幅度的百分之一比绝对阈值通用。成功率判断也比较严格允许 NMSE 小于 0.1但支撑集必须完全重合这个标准能过滤掉“幅值拟合好但支撑集差一个原子”的侥幸恢复。4.3 三算法对比测试脚本结构与输出判读我在对比 SBL、TSBL、TMSBL 时会分别构造单快照和多快照两组数据。单快照用 sbl 函数多快照用 mmv_sbl 和加上了关联矩阵学习的版本# 多快照数据8 个快照共享同一个支撑集幅值随时间缓慢变化 L 8 X_true np.zeros((N, L)) X_true[supp_true, :] rng.standard_normal((K, L)) # 用一阶平滑让相邻快照相关 for n in range(N): X_true[n, 1:] 0.7 * X_true[n, :-1] 0.3 * X_true[n, 1:] Y Phi X_true 0.01 * rng.standard_normal((M, L)) U_mmv, gamma_mmv, _ mmv_sbl(Y, Phi)输出判读时我会先看 gamma 的分布理想情况下活跃位置的 γ 比非活跃位置大几个量级。如果 γ 全部集中在同一个量级说明迭代没有进入稀疏状态这时先检查 σ² 是否被压得过小或测量矩阵没归一化。多快照场景下TMSBL 学出来的 B 矩阵非对角元如果接近零就说明数据其实没有时间相关性这时 TMSBL 的结果不会比共享支撑版本好。5. 避坑与排查收敛缓慢、支撑集漂移、病态矩阵的定位与修复下面这几条是我在测试中真正翻过车并定位到原因的按现象、原因、解决写可以直接对照排查。5.1 γ 不稀疏恢复出一堆小幅度值现象迭代结束后 gamma 没有形成两极分化x_hat 几乎所有位置都有非零小值。原因噪声方差 σ² 初值给得太大EM 早期把 x 的后验分布拉得很宽γ 更新一直被噪声主导另一个常见原因是 Φ 列未归一化量级差了几个数量级的列会把 γ 更新带偏。解决把 σ² 初值设成 var(y) / 10对 Φ 做列归一化给 γ_i 加一个 1e-8 的下界。改完这两处gamma 通常会在 50 次迭代内开始分化。5.2 σ² 迭代到零矩阵求逆出现奇异或 NaN现象迭代到中途 np.linalg.inv 报 LinAlgError 或返回 NaN。原因测量矩阵行数 M 较小、信噪比高时σ² 的 EM 更新会持续下降一旦小于浮点数精度ΦᵀΦ/σ² 就溢出。解决在 σ² 更新后加一行sigma2 max(sigma2, 1e-8)。更稳的做法是把 σ² 和 γ 的更新解耦γ 每轮更新σ² 每 10 轮更新一次实测能缓解早期震荡。5.3 复数数据处理结果全是错的现象把 SBL 直接用到雷达 IQ 或通信基带数据上恢复结果完全不对但同一段代码在实数仿真里正常。原因所有转置都要换成共轭转置。复数场景下 Φᵀ 应该是 Φᴴ Φ.conj().T后验均值、残差计算都受影响。解决把代码中的.T全部替换为.conj().Tγ 更新里的 mu**2 改成mu.real**2 mu.imag**2。另外噪声方差估计里残差要用np.sum(np.abs(R)**2)不是R R这种实数的点积写法。5.4 TMSBL 的 B 矩阵病态和单快照版本比反而变差现象加入相关矩阵学习后多快照恢复效果不如直接用共享支撑的 MMV-SBL。原因B 在迭代早期支撑集还没稳定时就开始学习错误的活跃原子把相关矩阵污染了或者是快照数 L 小于稀疏度 KB 的自由度过高估计方差极大。解决前 50 轮固定 B I等支撑集稳定后再启动 B 更新每次更新 B 后做对角加载加上 1e-6 的单位阵扰动。如果 L 小于 K建议退回共享支撑版本B 矩阵学不出来的。5.5 单次实验看着不错换一组随机 Φ 就翻车现象固定一个随机种子测试SBL 恢复得很漂亮把 Φ 换成另一组随机实现成功率断崖式下跌。原因压缩感知恢复单次结果本来就带方差少数幸运的测量矩阵会把问题变得异常简单单次实验不能代表算法水平。解决至少跑 100 次蒙特卡洛统计成功率或者 NMSE 中位数。我在测试脚本里固定一组随机种子做回归另开一组不固定种子的跑统计两个结果同时看。6. 进阶用法Woodbury 加速、超参数初始化与蒙特卡洛验证当 N 从几百涨到几千核心瓶颈是每次迭代的 N×N 协方差求逆。用 Woodbury 恒等式把求逆降到 M×MSBL 就能在中等维度下保持实用# 用 Woodbury 恒等式避免 N×N 直接求逆 D np.diag(gamma) PhiD Phi D A sigma2 * np.eye(M) PhiD Phi.T # M×M 求逆 A_inv np.linalg.inv(A) Sigma D - D Phi.T A_inv PhiD mu (D Phi.T) A_inv y这个版本里每次迭代要重建 D但求逆对象始终是 M×MN4096、M256 时速度提升非常明显。注意减法抵消问题当 γ_i 小于 1e-12 时Σ 的对角元可能出现微小负值把它截断到零再继续迭代。超参数初始化我也会调。γ 全 1 是安全默认值如果想让早期收敛更快可以把 γ 初始化为np.mean(Phi.T y)**2的每个列投影能量。σ² 初值我固定用 var(y)/10信噪比极高时改成 var(y)/100 也行但要配合 5.2 的截断。最后是验证习惯。我测试压缩感知算法时一定会跑蒙特卡洛固定 N、M、K随机生成 100 组不同的 Φ 和 x画一条成功率随 M/K 变化的曲线。单次实验的成功没有意义成功率曲线才决定这个算法能不能用。这个习惯帮我筛掉过不少看似效果惊艳、换参数就崩的改进版本。写到这里我最想强调的其实是最后一条贝叶斯方法的好处是参数自适应但代价是每一组数据都要真正迭代收敛才算数。建议你先把第 2 章的 sbl 函数跑通再按第 4 章的测试流程做一次成功率统计然后再进入 TSBL 和 TMSBL 的改造。希望帮到你。本文还有配套的精品资源点击获取