压缩感知稀疏贝叶斯：SBL到TMSBL原理与Python实现

简介：一套包含SBL、TSBL与TMSBL的压缩感知稀疏贝叶斯算法代码包，面向信号处理、压缩感知方向的研究者与工程师，解决低采样率下高维稀疏信号重构问题。三种算法各有侧重：SBL基于贝叶斯稀疏先验，TSBL引入时间相关性，TMSBL进一步支持多尺度分析，适合动态与非平稳信号。代码包共15个文件，以11个MATLAB脚本为主，含SBL/TSBL/TMSBL演示与实验脚本，另附2个PDF快速使用指南、1份ReadMe说明和1个docx文档，压缩包约479KB，结构直观，便于按需调用。作者已亲自测试，确保代码可直接运行。目前已有929人学习下载，适合需要复现算法、做对比实验或开展课程设计的读者。

1. 压缩感知稀疏贝叶斯：先想清楚它解决的是哪个环节

压缩感知里最常被问的一句话是：L1 范数都解得好好的，为什么还要碰稀疏贝叶斯？答案藏在两组场景里。一类是观测维度极低、信噪比又差，L1 类方法重建出来到处是伪峰；另一类是你要同时恢复多段共享稀疏结构的信号，比如多通道 EEG、阵列多快拍、宽带频谱感知，逐条用 L1 解不仅慢，还把通道间相关性丢光了。SBL（Sparse Bayesian Learning）在这两类场景下往往比 OMP、Lasso 更稳。它不是某种启发式阈值，而是把稀疏当作先验概率结构来学习，超参数由数据自己推断出来。TSBL 和 TMSBL 则是从单任务伸展出的多任务版本，前者共享稀疏支撑，后者额外建模时域相关性。这篇文章按 SBL → TSBL → TMSBL 的顺序，把数学上为什么有效、代码上怎么写、参数上怎么调一次讲透，适合正在做压缩感知重构、信道估计或阵列信号处理的人。

2. SBL 的核心机制：层级先验和自动相关性确定（ARD）

2.1 压缩感知重构问题与最大后验视角

压缩感知的重构目标是给定观测向量 y 和感知矩阵 A，求稀疏解 x，满足 y = Ax + n。A 是 M×N 的矩阵，M 通常远小于 N，n 是噪声。传统做法把问题写成 L1 正则化：

min_x ||y - Ax||_2^2 + λ||x||_1

从贝叶斯视角看，L1 等价于给 x 加了拉普拉斯先验。问题在于拉普拉斯先验只有一个尺度参数 λ，它对所有系数一视同仁，而且 λ 要手工调。稀疏贝叶斯换了一条路：不给 x 直接定一个稀疏分布，而是给 x 的方差设一个更上层的先验。x 还是零均值高斯分布，但每个系数的方差参数不同，由数据估计。这样差一点都不需要事先指定稀疏度，也不需要调 L1 惩罚系数。

这种分层建模在机器学习里叫层级先验（hierarchical prior），它的关键性质是：在边缘化掉方差参数之后，x 的边缘分布会变成重尾分布。重尾意味着大部分系数被压到接近零，少数系数可以很大，这正是稀疏信号该有的样子。SBL 的训练过程就是不断地根据当前 x 的后验估计去更新这些方差参数，让无关维度的方差收敛到零。

2.2 数学模型和概率图

SBL 的层级模型通常这样定义。第一层，x 服从零均值高斯分布：

p(x; γ) = ∏ N(x_i | 0, γ_i) i = 1..N

其中 γ_i 是第 i 个系数的先验方差，也叫超参数。第二层，γ 本身没有强约束，常被当作确定性参数用 EM 来估计。噪声也建模为高斯分布，方差为 σ²。于是后验分布 p(x | y; γ, σ²) 是高斯的，均值和协方差有闭式解：

μ = σ^{-2} Σ A^T y Σ = (σ^{-2} A^T A + Γ^{-1})^{-1}

其中 Γ = diag(γ_1, ..., γ_N)。这就是 SBL 的推断核心。给定当前的 γ 和 σ²，x 的后验均值 μ 就是重构结果，后验协方差 Σ 反映不确定性。更新 γ 时用 EM 或直接最大化边缘似然，常见迭代式是：

γ_i^(new) = μ_i^2 + Σ_ii

也就是后验二阶矩。Σ_ii 的存在让方差更新不会完全塌缩到零，数值上也更稳定。当某个 γ_i 持续变小时，对应列 A_i 对重构几乎没有贡献，这一点在 ARD（自动相关性确定）里被形式化描述。SBL 在压缩感知的语境下做的就是这个事：不是直接判断哪个系数稀疏，而是让证据决定哪些系数该从模型中移除。

2.3 为什么 SBL 容易收敛到稀疏解

仁科（Wipf）做过一个经典分析：SBL 的边缘似然函数在可辨识的情况下，局部最优解恰好在稀疏支撑集上，而 L1 正则化没有这种保证。更通俗的解释是，L1 惩罚是凸的但会让幅度整体收缩，稀疏度强的时候偏差明显；SBL 的非凸性虽然让优化路径更复杂，但其最优解的结构更干净。

实际信号处理里还有一个感知上的差异。OMP 这类贪婪算法一旦选错一列，后续很难修正；L1 内点法在 A 高相关时会出现支撑集抖动。SBL 由于在后验均值里保留了 Σ 的信息，相当于每个支撑候选都带着置信度，更新过程是连续平滑的，不会突然抛弃某个候选列。这在高相关感知矩阵下特别明显。

3. 手写一个能跑的 SBL：Python 实现、参数设置和收敛判断

3.1 最小可复现的 SBL 实现

下面给出一段可以直接跑的 Python 代码。它实现的是 EM 式的 SBL 迭代，没有调用任何专用稀疏贝叶斯库，依赖只有 NumPy。这里选择 EM 版本，一是容易对齐公式，二是可以清楚看到 γ 和 σ² 的更新位置。

import numpy as np def sbl_reconstruct(A, y, max_iter=200, tol=1e-8, noise_var=None): """ 压缩感知SBL重构 A: M x N 感知矩阵 y: M x 1 观测向量 max_iter: 最大迭代次数 tol: gamma变化相对容差 noise_var: 噪声方差，None表示自动估计 """ M, N = A.shape # 初始化：每个系数的先验方差设为1 gamma = np.ones((N, 1)) if noise_var is None: noise_var = np.var(y) * 0.1 # 保守初值 for it in range(max_iter): # 计算后验方差 Sigma gamma_inv = 1.0 / gamma Sigma = np.linalg.inv((A.T @ A) / noise_var + np.diag(gamma_inv[:, 0])) # 计算后验均值 mu mu = (Sigma @ A.T @ y) / noise_var # 更新 gamma：后验二阶矩 gamma_new = (mu ** 2 + np.diag(Sigma).reshape(-1, 1)) # 更新噪声方差（EM步骤） if noise_var is not None: residual = y - A @ mu noise_var = (np.linalg.norm(residual) ** 2 + np.trace(Sigma @ (A.T @ A))) / M # 判断收敛 delta = np.linalg.norm(gamma_new - gamma) / np.linalg.norm(gamma) gamma = gamma_new if delta < tol: break return mu, gamma, noise_var

这个实现里，后验均值和协方差完全按前面的闭式公式计算。矩阵求逆用了np.linalg.inv，当 N 上千时这一步会是瓶颈。gamma 初始化为全 1 是常见做法，实践中如果知道信号幅度量级，可以更接近真实方差，但初始化对所有维度一样，算法依然能自动把不相关维度压到零。噪声方差的更新里多了一项trace(Sigma @ (A.T @ A))，这是对残差的补偿，没这一项噪声方差会越估越小，导致过拟合。

3.2 参数怎么设才能不踩坑

SBL 的核心参数其实只有三个：最大迭代次数、收敛容差和噪声方差初值。迭代次数建议 100 到 500，因为一个维度从 1 收敛到接近 0 需要若干轮，支持集在大维度下是渐进激活的。收敛容差设 1e-8 到 1e-10 较为安全，如果设太大可能支撑集还没定型就停了。噪声方差初值影响路径但不太影响终值，一般取观测信号方差的 0.1 到 0.01 倍。

一个更微妙的点：感知矩阵 A 的列归一化。如果某些列能量远高于其他列，SBL 会偏向这些列，因为同样的 γ 变化在似然里产生的效果被放大了。我一般会在预处理里对 A 做列归一化，重构完成后把支撑系数按列缩放恢复回去。这个操作对 SBL 比对 L1 影响更明显，因为 γ 更新的分母对列灵敏度不同。表格列出几个关键控制项：

3.3 慢收敛与数值崩溃的处理

亲手跑 SBL 最容易遇到两个现象：gamma 在某几步剧烈震荡，或者一串维度持续不归零。震荡通常发生在噪声方差更新过快时，可以给噪声方差加指数加权，例如noise_var = 0.7 * new + 0.3 * old。持续不归零则要检查 A 是否存在近似线性相关的列，这会使得 Σ 对角项永远偏大，γ 的更新始终大于零。遇到这种情况可以对 γ 设一个下限，比如 1e-12，同时把最终支撑选择交给后处理的峰值检测。

数值崩溃则是 N 大于 2000 时的常客。直接求逆一个 N×N 矩阵在每次迭代里做一遍，复杂度 O(N³)。这里有两个思路：第一，如果感知矩阵是正交基或傅里叶部分采样，用矩阵恒等式把 Σ 的计算降到 M×M；第二，每次迭代后把 γ 小于某个相对阈值的维度直接剪掉，只在活跃集上求逆。这两种方式在工程里都是常见做法，前者适合结构化的 A，后者适合任意 A。

4. 从 SBL 扩展到 TSBL 和 TMSBL：多测量向量与时域结构建模

4.1 多任务压缩感知和 TSBL 的共享支撑假设

很多实际压缩感知场景不是解一个 y，而是解一组观测。比如一个传感器阵列在 L 个快拍下得到 y_1, ..., y_L，它们对应的稀疏信号 x_1, ..., x_L 在不同时刻幅度不同，但非零位置高度一致。逐列用 SBL 解会得到 L 个独立的支撑集，轻微噪声就会让支撑集抖动，重构结果无谓地出现虚假峰。TSBL（Temporal Sparse Bayesian Learning，或 T 取多任务 Multi-Task 语义）的核心是把这 L 个任务联合建模，共享同一个超参数 γ。

TSBL 的模型假设每个任务的先验方差都由同一个 γ 控制：

p(x_l | γ) = N(x_l | 0, Γ), l = 1..L

不同任务的观测方程分别是 y_l = A_l x_l + n_l。这里各任务感知矩阵可以不同，只要稀疏支撑一致。联合边缘似然写成所有任务后验的乘积，EM 更新时 γ 的迭代式为：

γ_i = (1/L) Σ_l ( μ_{l,i}^2 + Σ_{l,ii} )

也就是说，某个维度是否保留，要看它在所有任务里的平均能量。一个任务里噪声引起的单次伪峰，不会让该维度混进支撑集，因为其他任务会给它投票。这种多任务联合重构比逐列 SBL 在支撑检测上的优势，在小快拍、低信噪比下尤为突出。

4.2 TMSBL 在 TSBL 上加了什么

TSBL 假设了共享支撑，却没有建模同一维度在时间上的幅度相关性。实际信号里，相邻快拍的稀疏系数往往是平滑变化的，比如慢变信道或窄带信号。TMSBL（Temporal Multiple Sparse Bayesian Learning）把每个维度的 L 个快拍看作一个时间序列，给这一整条序列设定一个 L×L 的协方差矩阵 B：

p(x^i | γ_i, B) = N(x^i | 0, γ_i B), i = 1..N

这里的 x^i 是第 i 个系数在所有快拍上的取值向量。B 反映了该系数随时间的相关结构，比如指数相关 B_{pq} = ρ^{|p-q|}。与 TSBL 相比，TMSBL 多了一个 B 的学习或设定，γ 的更新不再是对每个任务后验方差直接平均，而是考虑时间协方差造成的耦合。B 的引入带来了更紧凑的表示：如果时间上高度相关，等效的自由参数更少，重构所需的观测总量可以进一步下降。

TMSBL 的 EM 更新里，γ 和 B 交替更新。B 的估计通常借助残差协方差和中间量构造，常见更新涉及对时间和任务两个维度的矩阵运算。一个直观印象是，TMSBL 可以理解为在共享支撑之上加一层时域去相关，它比 TSBL 更擅长利用信号本身的平滑性，因此适合观测快拍较长、时间相关性强的数据。

4.3 TSBL 和 TMSBL 的 Python 骨架

直接给一个能看懂的 TSBL 核心迭代骨架。它与单任务 SBL 的差别集中在 gamma 更新和多个任务后验统计量汇总上。

def tsbl_reconstruct(A_list, y_list, max_iter=200, tol=1e-8): L = len(y_list) N = A_list[0].shape[1] gamma = np.ones((N, 1)) # 每个任务独立维护后验统计量 mu_list = [] Sigma_list = [] for _ in range(max_iter): mu_list = [] Sigma_list = [] gamma_inv = 1.0 / gamma for l in range(L): A = A_list[l] y = y_list[l] Sigma = np.linalg.inv(A.T @ A + np.diag(gamma_inv[:, 0])) mu = Sigma @ A.T @ y mu_list.append(mu) Sigma_list.append(Sigma) gamma_new = np.zeros((N, 1)) for l in range(L): gamma_new += mu_list[l] ** 2 + np.diag(Sigma_list[l]).reshape(-1, 1) gamma_new /= L delta = np.linalg.norm(gamma_new - gamma) / np.linalg.norm(gamma) gamma = gamma_new if delta < tol: break return mu_list, gamma

上面这段假设噪声方差已知并做了单位化简化，真实场景要像单任务 SBL 一样加噪声方差估计。gamma_new的求法值得注意：它是所有任务二阶矩的平均，而不是.mean(axis=0)到别的统计量。这种平均给了所有任务对支撑的否决权，如果一个维度在某个任务里方差很大但在其他任务为零，平均值会把它拉低，这就是多任务鲁棒性的来源。

TMSBL 的代码骨架比 TSBL 多一个 B 矩阵。B 的更新公式依赖后验协方差块与均值外积的时域聚合，通常写成：

B = (1/N) * sum_i ( Mu[i] @ Mu[i].T + Sigma_i )

这里的 Mu[i] 是第 i 个系数在所有快拍上的后验均值拼接。实现时把 x 组织成 N×L 的矩阵，方便按行计算时域外积。B 初始化成单位阵时退化成 TSBL，迭代若干轮后 B 的非对角元素会体现时域相关强度。

4.4 三种算法怎么选

一句话：信号只有单次观测用 SBL；多次观测、支撑一致但幅度独立用 TSBL；多次观测、幅度在时间上平滑变化用 TMSBL。很多人误以为 TMSBL 在所有场景都优于 TSBL，实际如果时间相邻快拍之间的相关性弱，B 估计会带来额外方差，在多观测但独立场景下 TSBL 反而更稳。选型时先做一步简单的相关性分析：把逐任务 SBL 的结果按维度对齐，计算相邻快拍间稀疏系数的相关系数，超过 0.5 再上 TMSBL。另外，TMSBL 对 B 的维度 L×L 求逆发生在每次迭代里，L 很大的时候要注意矩阵求逆的开销，必要时把 B 限制为托普利兹结构。

5. 验证 SBL/TSBL/TMSBL 的对比实验与三个提效习惯

5.1 最小验证流程

拿到代码第一件事不是在真实数据上跑，而是先做合成实验。最常见的做法是生成高斯随机感知矩阵 A，随机生成 K 稀疏的 x，计算 y = Ax + n。以 M=64，N=256，K=8 为一组参数，信噪比 20dB。对单任务 SBL，计算重构误差norm(x_true - x_hat)/norm(x_true)，同时记录迭代轮数。对 TSBL 和 TMSBL，生成 L=8 个快拍，共享支撑，幅度在时间上做一阶自回归相关，然后对比三种算法的支撑恢复率与 MSE。支撑恢复率定义为正确找回的非零索引比例，比 MSE 更能反映稀疏重构质量。

5.2 三个工程化提效习惯

第一个习惯是强制 gamma 剪枝。SBL 类算法迭代过程中大量 γ_i 会缓慢逼近零，但不会真正等于零。每 20 轮做一次活跃集检测，把 γ_i 小于最大 γ 的 1e-4 倍对应列直接从 A 里移除，之后只在活跃集上更新。这个操作能让 N=5000 的问题迭代速度提升一个数量级，且几乎不影响重构精度。注意剪枝要设置回退机制，防止活跃集被过早固定。第二个习惯是监控边缘似然的增量而不是 gamma 变化。gamma 的变化尺度在不同维度不同，直接比较范数容易误判收敛。计算 log 边缘似然的增量，连续几次小于 1e-3 就停止，更稳定。第三个习惯是给感知矩阵做 QR 预分解。如果多次使用同一个 A，把A.T @ A和A.T @ y预先算好缓存，迭代里只做矩阵乘和求逆，省去反复计算 Gram 矩阵的开销。

5.3 调参失败的快速定位

重构结果稀疏度够了但幅度偏小，先查噪声方差估计，大概率是噪声方差估得过大，后验均值被压向零；支撑集频繁在相邻两个维度间跳变，查感知矩阵的列相关性，这是在高相关列下 SBL 支撑不稳定，给 γ 更新加动量或降低学习率；多任务联合重构比单任务还差，先确认各任务感知矩阵是否对齐，感知矩阵不一致时 TSBL 的共享支撑假设本身就不成立。SBL 类算法的收敛没有 L1 那么直白，把每一轮 γ 变化画成热力图，一眼就能看出维度是逐步归零还是抖动，这是排查参数最有效的工具。

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