matlab基础贝叶斯变换的压缩感知
一、核心原理：从优化到概率推断
1.1 贝叶斯基本公式
传统压缩感知模型：y = Φx + e，其中 y 为观测向量 (M×1)，Φ 为测量矩阵 (M×N, M<<N)，x 为稀疏信号 (N×1)，e 为噪声。
贝叶斯方法将信号 x 和噪声方差 σ² 都视为随机变量，通过贝叶斯定理计算后验分布：
p(x, σ² | y) ∝ p(y | x, σ²) × p(x) × p(σ²)
其中：
- 似然函数
p(y | x, σ²)：描述观测模型，通常为高斯分布N(Φx, σ²I) - 先验分布
p(x)：编码信号的稀疏性知识（关键所在） - 超参数先验
p(σ²)：对噪声方差的先验假设
1.2 稀疏信号先验的两种主要建模方式
| 先验类型 | 数学形式 | 特点 | 适用场景 |
|---|---|---|---|
| 拉普拉斯先验 | `p(x_i) = (λ/2)exp(-λ | x_i | )` |
| 高斯尺度混合 | p(x_i) = ∫ N(0, γ_i)p(γ_i)dγ_i |
可学习个体稀疏性，更灵活 | 结构化稀疏信号 |
二、算法实现：稀疏贝叶斯学习（SBL）
SBL使用高斯尺度混合先验，为每个系数 x_i 分配独立方差 γ_i，并通过证据最大化自动学习这些参数。
2.1 基础SBL算法实现
function [x_est, gamma_est, sigma2_est] = SBL_Basic(y, Phi, max_iter, tolerance)
% 稀疏贝叶斯学习（SBL）基础实现
% 输入：
% y - 观测向量 (M×1)
% Phi - 测量矩阵 (M×N)
% max_iter - 最大迭代次数
% tolerance - 收敛阈值
% 输出：
% x_est - 信号估计 (N×1)
% gamma_est - 方差估计 (N×1)
% sigma2_est - 噪声方差估计
[M, N] = size(Phi);
% 初始化参数
gamma = ones(N, 1); % 信号方差的初始值
sigma2 = var(y) * 0.1; % 噪声方差的初始估计
x_est = zeros(N, 1); % 信号估计
% 预计算，提高效率
PhiT = Phi';
PhiT_y = PhiT * y;
% 迭代优化
for iter = 1:max_iter
% ========== E步：计算后验均值和方差 ==========
% 计算后验协方差矩阵 Σ = (σ^-2 Φ^T Φ + Γ^-1)^-1
% 使用高效计算避免直接求N×N矩阵的逆
Gamma_inv = diag(1 ./ gamma);
Sigma_inv = (PhiT * Phi) / sigma2 + Gamma_inv;
% 使用Cholesky分解稳定求逆
[R, p] = chol(Sigma_inv);
if p > 0
% 如果非正定，添加正则化项
Sigma_inv = Sigma_inv + eye(N) * 1e-8;
R = chol(Sigma_inv);
end
Sigma = inv_chol(R); % 使用Cholesky分解求逆
% 计算后验均值 μ = σ^-2 Σ Φ^T y
mu = Sigma * PhiT_y / sigma2;
% ========== M步：更新超参数 ==========
% 更新信号方差 γ_i
for i = 1:N
gamma(i) = Sigma(i, i) + mu(i)^2;
end
% 更新噪声方差 σ^2
y_pred = Phi * mu;
error = y - y_pred;
sigma2_new = (error' * error) / M + sigma2 * trace(Phi * Sigma * PhiT) / M;
% ========== 收敛检查 ==========
gamma_change = norm(gamma - gamma_old) / norm(gamma_old);
sigma_change = abs(sigma2_new - sigma2) / sigma2;
sigma2 = sigma2_new;
x_est = mu; % 当前信号估计
% 打印迭代信息
if mod(iter, 10) == 0
sparsity = sum(gamma > 1e-3 * max(gamma)) / N;
fprintf('迭代 %d: σ²=%.4e, 稀疏度=%.2f%%, 变化量=[γ:%.2e, σ²:%.2e]\n', ...
iter, sigma2, sparsity*100, gamma_change, sigma_change);
end
% 检查收敛
if max(gamma_change, sigma_change) < tolerance
fprintf('算法在第%d次迭代收敛\n', iter);
break;
end
end
gamma_est = gamma;
sigma2_est = sigma2;
end
% 使用Cholesky分解高效求逆
function A_inv = inv_chol(R)
% R是上三角Cholesky因子，满足 R' * R = A
N = size(R, 1);
% 解 R' * R * X = I
X = eye(N);
A_inv = R \ (R' \ X);
end
2.2 完整的仿真测试框架
%% 贝叶斯压缩感知完整仿真示例
clear; close all; clc;
% ========== 1. 生成稀疏测试信号 ==========
N = 256; % 信号长度
K = 12; % 稀疏度（非零元个数）
x_true = zeros(N, 1);
% 随机选择非零位置
nonzero_indices = randperm(N, K);
% 生成非零值（高斯随机或正负脉冲）
x_true(nonzero_indices) = randn(K, 1) .* (1 + 0.5*rand(K, 1));
% ========== 2. 压缩观测 ==========
M = 64; % 观测数（压缩比 M/N = 25%）
% 随机高斯测量矩阵
Phi = randn(M, N) / sqrt(M); % 归一化
% 生成无噪声观测
y_clean = Phi * x_true;
% 添加高斯噪声
SNR_dB = 20; % 信噪比
signal_power = norm(y_clean)^2 / M;
noise_power = signal_power / (10^(SNR_dB/10));
noise = sqrt(noise_power) * randn(M, 1);
y = y_clean + noise;
% ========== 3. 对比不同重建方法 ==========
fprintf('=== 压缩感知重建对比测试 ===\n');
fprintf('信号维度: N=%d, 观测数: M=%d, 稀疏度: K=%d, SNR: %d dB\n', N, M, K, SNR_dB);
% 方法1：传统L1最小化（BPDN）
fprintf('\n1. 运行L1最小化 (BPDN)...\n');
lambda = 0.1 * norm(Phi' * y, 'inf'); % 正则化参数
x_l1 = l1_bpdn(y, Phi, lambda);
% 方法2：正交匹配追踪（OMP）
fprintf('2. 运行OMP算法...\n');
x_omp = omp(y, Phi, K);
% 方法3：稀疏贝叶斯学习（SBL）
fprintf('3. 运行稀疏贝叶斯学习 (SBL)...\n');
tic;
[x_sbl, gamma_sbl, sigma2_sbl] = SBL_Basic(y, Phi, 200, 1e-6);
time_sbl = toc;
% ========== 4. 性能评估 ==========
% 计算重建误差
mse_l1 = 10*log10(mean((x_l1 - x_true).^2) / mean(x_true.^2));
mse_omp = 10*log10(mean((x_omp - x_true).^2) / mean(x_true.^2));
mse_sbl = 10*log10(mean((x_sbl - x_true).^2) / mean(x_true.^2));
% 计算支持集恢复精度
threshold = 0.1 * max(abs(x_true));
true_support = find(abs(x_true) > threshold);
est_support_l1 = find(abs(x_l1) > threshold);
est_support_omp = find(abs(x_omp) > threshold);
est_support_sbl = find(abs(x_sbl) > threshold);
% 计算精确恢复的指标
precision_l1 = length(intersect(true_support, est_support_l1)) / length(est_support_l1);
recall_l1 = length(intersect(true_support, est_support_l1)) / length(true_support);
precision_sbl = length(intersect(true_support, est_support_sbl)) / length(est_support_sbl);
recall_sbl = length(intersect(true_support, est_support_sbl)) / length(true_support);
% ========== 5. 可视化结果 ==========
figure('Position', [100, 100, 1400, 900]);
% 子图1：原始稀疏信号
subplot(3, 3, 1);
stem(x_true, 'b.', 'MarkerSize', 10); grid on;
xlim([0, N]); title(sprintf('原始信号 (N=%d, K=%d)', N, K));
xlabel('索引'); ylabel('幅值');
% 子图2：L1最小化重建
subplot(3, 3, 2);
stem(x_l1, 'r.', 'MarkerSize', 10); grid on;
xlim([0, N]); title(sprintf('L1最小化重建\n相对MSE: %.2f dB', mse_l1));
xlabel('索引'); ylabel('幅值');
% 子图3：OMP重建
subplot(3, 3, 3);
stem(x_omp, 'g.', 'MarkerSize', 10); grid on;
xlim([0, N]); title(sprintf('OMP重建\n相对MSE: %.2f dB', mse_omp));
xlabel('索引'); ylabel('幅值');
% 子图4：SBL重建
subplot(3, 3, 4);
stem(x_sbl, 'm.', 'MarkerSize', 10); grid on;
xlim([0, N]);
title(sprintf('SBL重建 (%.2f秒)\n相对MSE: %.2f dB', time_sbl, mse_sbl));
xlabel('索引'); ylabel('幅值');
% 子图5：方差估计（不确定性量化）
subplot(3, 3, 5);
bar(gamma_sbl, 'FaceAlpha', 0.6); grid on;
xlim([0, N]); title('SBL估计的方差 γ_i');
xlabel('索引'); ylabel('方差估计');
% 子图6：噪声方差学习曲线（可扩展部分）
subplot(3, 3, 6);
% 实际中需要记录每次迭代的sigma2值
plot(1:length(sigma2_history), 10*log10(sigma2_history), 'LineWidth', 2);
grid on; title('噪声方差 σ² 学习过程');
xlabel('迭代次数'); ylabel('σ² (dB)');
% 子图7-9：残差与不确定性分析
subplot(3, 3, 7);
% 计算后验标准差
posterior_std = sqrt(diag(Sigma));
errorbar(1:20, x_sbl(1:20), 2*posterior_std(1:20), 'o', 'LineWidth', 1.5);
hold on; plot(x_true(1:20), 'r*', 'MarkerSize', 8);
grid on; legend('估计值±2σ', '真实值', 'Location', 'best');
title('估计不确定性（前20个系数）'); xlabel('索引'); ylabel('幅值');
subplot(3, 3, 8);
residual = y - Phi * x_sbl;
histogram(residual, 30, 'Normalization', 'pdf'); grid on;
hold on;
% 绘制拟合的高斯分布
x_fit = linspace(min(residual), max(residual), 100);
pdf_fit = normpdf(x_fit, 0, sqrt(sigma2_sbl));
plot(x_fit, pdf_fit, 'r-', 'LineWidth', 2);
title(sprintf('残差分布 (估计σ²=%.3e)', sigma2_sbl));
xlabel('残差'); ylabel('概率密度');
subplot(3, 3, 9);
% 支持集恢复性能对比
methods = {'L1最小化', 'SBL'};
precision = [precision_l1, precision_sbl];
recall = [recall_l1, recall_sbl];
bar_data = [precision; recall]';
h = bar(1:2, bar_data);
set(h(1), 'FaceColor', [0.2, 0.6, 0.8]);
set(h(2), 'FaceColor', [0.8, 0.4, 0.2]);
set(gca, 'XTickLabel', methods);
ylabel('百分比'); title('支持集恢复性能');
legend('精确率', '召回率', 'Location', 'best');
grid on;
fprintf('\n=== 性能总结 ===\n');
fprintf('方法\t\t相对MSE(dB)\t精确率\t\t召回率\t\t计算时间\n');
fprintf('L1最小化\t%.2f\t\t%.2f\t\t%.2f\t\tN/A\n', mse_l1, precision_l1, recall_l1);
fprintf('OMP\t\t%.2f\t\tN/A\t\tN/A\t\tN/A\n', mse_omp);
fprintf('SBL\t\t%.2f\t\t%.2f\t\t%.2f\t\t%.2f秒\n', mse_sbl, precision_sbl, recall_sbl, time_sbl);
2.3 实用辅助函数
% L1基追踪去噪（BPDN）实现
function x = l1_bpdn(y, A, lambda, max_iter)
% 使用迭代软阈值算法（ISTA）
if nargin < 4, max_iter = 1000; end
[M, N] = size(A);
x = zeros(N, 1);
L = norm(A' * A, 2); % Lipschitz常数
for iter = 1:max_iter
grad = A' * (A * x - y);
x_temp = x - (1/L) * grad;
% 软阈值
x_new = sign(x_temp) .* max(abs(x_temp) - lambda/L, 0);
if norm(x_new - x) / norm(x) < 1e-6
break;
end
x = x_new;
end
end
% OMP算法实现
function x = omp(y, A, K)
[M, N] = size(A);
x = zeros(N, 1);
residual = y;
support = [];
for k = 1:K
% 找到与残差最相关的列
corr = abs(A' * residual);
[~, idx] = max(corr);
support = union(support, idx);
% 在选定的支持集上最小二乘
As = A(:, support);
x_temp = As \ y;
% 更新残差
residual = y - As * x_temp;
if norm(residual) < 1e-6
break;
end
end
x(support) = As \ y;
end
三、贝叶斯压缩感知的优势与扩展
3.1 与传统方法的对比
| 特性 | 传统L1最小化 | 稀疏贝叶斯学习（SBL） |
|---|---|---|
| 参数选择 | 需手动调λ | 自动学习所有超参数 |
| 不确定性量化 | 无 | 提供后验方差估计 |
| 计算复杂度 | O(N³) | O(N³)但可优化 |
| 抗噪能力 | 中等 | 强（自动估计噪声） |
| 灵活性 | 固定稀疏性 | 可学习结构化稀疏性 |
3.2 高级扩展：快速SBL算法
function [x_est, gamma_est] = Fast_SBL(y, Phi, max_iter)
% 使用Tipping的快速SBL算法
[M, N] = size(Phi);
% 初始化
gamma = ones(N, 1);
keep_list = (1:N)';
for iter = 1:max_iter
% 当前活动集
Phi_active = Phi(:, keep_list);
gamma_active = gamma(keep_list);
% 计算后验统计量
Sigma_inv = Phi_active' * Phi_active / sigma2 + diag(1./gamma_active);
mu = (Sigma_inv \ Phi_active') * y / sigma2;
% 更新gamma
gamma_active_new = abs(mu).^2 + diag(inv(Sigma_inv));
gamma(keep_list) = gamma_active_new;
% 剪枝小gamma对应的基
prune_threshold = 1e3 * max(gamma_active_new);
keep_indices = find(gamma_active_new < prune_threshold);
keep_list = keep_list(keep_indices);
if isempty(keep_list)
break;
end
end
x_est = zeros(N, 1);
x_est(keep_list) = mu;
gamma_est = gamma;
end
3.3 针对二维图像的应用调整
function [image_est, gamma_map] = SBL_2D(measurements, sensing_matrix, image_size)
% 二维图像的贝叶斯压缩感知
% sensing_matrix可以是随机投影、部分傅里叶等
% 将图像向量化
[M, N] = size(sensing_matrix);
% 使用小波稀疏基替代标准基
Psi = construct_wavelet_dict(image_size, 'db4', 3);
A = sensing_matrix * Psi; % 组合矩阵
% 运行SBL
[coeffs_est, gamma_est] = SBL_Basic(measurements, A, 150, 1e-6);
% 重建图像
image_est = reshape(Psi * coeffs_est, image_size);
gamma_map = reshape(gamma_est, image_size);
end
参考代码 基础贝叶斯变换的压缩感知 www.youwenfan.com/contentcno/96285.html
四、实际应用建议
4.1 关键参数设置指南
-
初始化策略：
gamma：初始化为ones(N,1)或从abs(Φ'y)估计sigma2：初始化为var(y)/10，避免过大
-
收敛准则：
- 相对参数变化<
1e-6 - 最大迭代100-200次（复杂场景可增加）
- 相对参数变化<
-
稀疏先验选择：
- 拉普拉斯先验：快速简单
- 高斯尺度混合：精度更高
- 伯努利-高斯混合：适合确切稀疏性
4.2 性能提升技巧
% 技巧1：预热策略
function result = warm_start_SBL(y, Phi, lambda_list)
% 使用多个lambda值预热
results = cell(length(lambda_list), 1);
for i = 1:length(lambda_list)
x_temp = l1_bpdn(y, Phi, lambda_list(i));
% 用L1结果初始化SBL
gamma_init = abs(x_temp).^2 + 1e-3;
results{i} = SBL_with_init(y, Phi, gamma_init);
end
% 选择最佳结果
[~, best_idx] = min([results{:}.mse]);
result = results{best_idx};
end
% 技巧2：分块处理大规模问题
function result = block_SBL(y, Phi, block_size)
% 分块处理大规模信号
N = size(Phi, 2);
num_blocks = ceil(N / block_size);
x_est = zeros(N, 1);
gamma_est = zeros(N, 1);
for b = 1:num_blocks
idx = (1:block_size) + (b-1)*block_size;
idx(idx > N) = [];
Phi_b = Phi(:, idx);
[x_b, gamma_b] = SBL_Basic(y, Phi_b, 100, 1e-5);
x_est(idx) = x_b;
gamma_est(idx) = gamma_b;
end
end
4.3 常见问题诊断
| 问题现象 | 可能原因 | 解决方案 |
|---|---|---|
| 收敛慢 | 初始值不合适 | 使用L1结果初始化 |
| 过拟合 | 噪声方差低估 | 增加sigma2先验的尺度 |
| 欠拟合 | 稀疏性过强 | 减小gamma的初始值 |
| 数值不稳定 | 矩阵病态 | 添加1e-8*eye(N)正则项 |
贝叶斯压缩感知通过概率框架将参数学习、不确定性量化和自适应稀疏性融为一体，特别适合信噪比较低、信号结构复杂的实际情况。虽然计算量通常大于传统优化方法，但其提供的后验分布和自动参数调节能力在许多应用中具有不可替代的价值。
对于您特定的应用场景，我可以进一步提供：
- 非高斯噪声模型下的贝叶斯压缩感知
- 结构化稀疏先验（如分组稀疏、树状稀疏）
- 在线/增量式贝叶斯压缩感知实现
浙公网安备 33010602011771号