原子范数前向后贪婪算法（FoBa-AN）实现

最新推荐文章于 2026-06-20 22:17:21 发布

原创最新推荐文章于 2026-06-20 22:17:21 发布 · 194 阅读

5 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

标签

#算法

原子范数前向后贪婪算法（FoBa-AN）实现，以线谱估计为例（原子为复指数向量）

代码 `foba_an_line_spectrum.m`

function [x_est, S, w] = foba_an_line_spectrum(A, b, lambda, Kmax, gamma, eps)
% FOBA_AN_LINE_SPECTRUM  原子范数前向后贪婪算法（线谱原子）
%
% 输入：
%   A      : m×n 测量矩阵（实或复）
%   b      : m×1 观测向量
%   lambda : 正则化参数
%   Kmax   : 最大原子数（可选，默认 20）
%   gamma  : 前向停止阈值（可选，默认 1e-3）
%   eps    : 后向删除容忍度（可选，默认 1e-6）
%
% 输出：
%   x_est  : n×1 恢复信号
%   S      : 所选原子的频率索引（归一化频率 0~1）
%   w      : 对应的系数向量
%
% 原子集：a(f) = [1, exp(1i*2π*f), ..., exp(1i*2π*f*(n-1))]' / sqrt(n)
%         其中 f ∈ [0,1) 连续

% --- 默认参数 ---
if nargin < 4, Kmax = 20; end
if nargin < 5, gamma = 1e-3; end
if nargin < 6, eps = 1e-6; end

[m, n] = size(A);

% --- 初始化 ---
S = [];            % 已选原子频率列表
w = [];            % 对应系数
x_est = zeros(n,1);
r = b;             % 残差

% --- 主循环 ---
while length(S) < Kmax

    %% ── 1. 前向步：找与残差最相关的原子频率 ──
    % 计算相关度函数 g(f) = |<A'*r, a(f)>|
    % 首先用 FFT 粗略定位峰值（精度 1/n），再牛顿精化
    Ar = A' * r;                        % n×1
    % 对 Ar 补零做高密度 FFT 找初始峰值（此处用 n*8 点）
    Nfft = n * 32;                      % 过采样因子
    X = fft(Ar, Nfft);
    X = X(:);
    [~, idx] = max(abs(X));
    f0 = (idx-1)/Nfft;                  % 粗略频率

    % 牛顿法精化（最大化 |<Ar, a(f)>|^2）
    % 目标函数：J(f) = |a(f)' * Ar|^2
    % 导数与海森矩阵可用解析公式（略去细节，用有限差分）
    f = f0;
    for iter = 1:10
        a = exp(1i*2*pi*f*(0:n-1)') / sqrt(n);
        grad = 2*real( conj(a'*Ar) * (a'*Ar) * (-1i*2*pi*(0:n-1)' .* conj(a)) );
        hess = 2*real( conj(a'*Ar) * (a'*Ar) * (-4*pi^2*(0:n-1)'.^2 .* conj(a)) ) ...
               - 2*abs(a'*Ar)^2 * 4*pi^2 * mean((0:n-1).^2);
        df = -grad / (hess + eps);      % 牛顿步
        f = mod(f + df, 1);
        if abs(df) < 1e-12, break; end
    end

    % 检查是否值得加入
    a_new = exp(1i*2*pi*f*(0:n-1)') / sqrt(n);
    corr_val = abs(a_new' * Ar);
    if corr_val < gamma
        break;
    end

    % 加入新原子
    S = [S; f];

    %% ── 2. 凸拟合：在已选原子集上解 LASSO ──
    % 构建字典 D (n × |S|)
    D = exp(1i*2*pi*S(:)' .* (0:n-1)') / sqrt(n);  % n×|S|
    AD = A * D;                                      % m×|S|

    % 坐标下降解 LASSO: min_w 0.5*||AD*w - b||^2 + lambda*||w||_1
    w = lasso_cd(AD, b, lambda, w);

    % 剔除接近零的系数（数值稳定性）
    tol_zero = 1e-8;
    keep = abs(w) > tol_zero;
    S = S(keep);
    w = w(keep);

    %% ── 3. 后向步：尝试删除冗余原子 ──
    changed = true;
    while changed
        changed = false;
        best_idx = [];
        best_df = inf;
        for j = 1:length(S)
            S_try = S; S_try(j) = [];
            D_try = exp(1i*2*pi*S_try(:)' .* (0:n-1)') / sqrt(n);
            AD_try = A * D_try;
            w_try = lasso_cd(AD_try, b, lambda);
            pred_try = AD_try * w_try;
            f_try = 0.5 * norm(b - pred_try)^2 + lambda * sum(abs(w_try));

            % 当前目标值
            pred_cur = AD * w;
            f_cur = 0.5 * norm(b - pred_cur)^2 + lambda * sum(abs(w));

            df = f_try - f_cur;
            if df < best_df
                best_df = df;
                best_idx = j;
            end
        end
        if best_df < eps && ~isempty(best_idx)
            % 删除该原子
            S(best_idx) = [];
            D = exp(1i*2*pi*S(:)' .* (0:n-1)') / sqrt(n);
            AD = A * D;
            w = lasso_cd(AD, b, lambda);
            changed = true;
        end
    end

    %% ── 更新残差 ──
    x_est = D * w;
    r = b - A * x_est;

    % 显示进度
    fprintf('Iter %d: |S|=%d, obj=%.4e\n', length(S), length(S), ...
        0.5*norm(r)^2 + lambda*sum(abs(w)));
end

% --- 辅助函数：坐标下降解 LASSO ---
function w = lasso_cd(AD, b, lambda, w_init)
% 简单坐标下降，适用于复数
    [m, p] = size(AD);
    if nargin < 4 || isempty(w_init)
        w = zeros(p,1);
    else
        w = w_init;
    end
    r = b - AD * w;
    max_iter = 500;
    for it = 1:max_iter
        w_old = w;
        for j = 1:p
            aj = AD(:,j);
            rho = aj' * aj;               % 标量
            if abs(rho) < 1e-14, continue; end
            z = aj' * (r + aj * w(j));     % 去掉当前变量的贡献
            w(j) = soft_thresh(z / rho, lambda / rho);
            r = r - aj * (w(j) - w_old(j));
        end
        if norm(w - w_old) / (norm(w)+eps) < 1e-8
            break;
        end
    end
end

function s = soft_thresh(z, tau)
% 复数软阈值
    if abs(z) <= tau
        s = 0;
    else
        s = (abs(z) - tau) * sign(z);
    end
end

end

使用示例（测试脚本）

%% 测试 FoBa-AN 线谱恢复
clc; clear; close all;

% --- 生成真实信号 ---
n = 64;                     % 信号长度
true_freqs = [0.13, 0.37, 0.52];  % 三个频率
true_coeffs = [2, 1.5, 1];       % 幅度（复）
x_true = zeros(n,1);
for k = 1:length(true_freqs)
    a = exp(1i*2*pi*true_freqs(k)*(0:n-1)') / sqrt(n);
    x_true = x_true + true_coeffs(k) * a;
end

% --- 随机测量矩阵（欠采样）---
m = 48;                     % 测量数（小于 n，压缩感知）
A = randn(m,n) + 1i*randn(m,n);
A = A / sqrt(2*m);          % 归一化行范数
b = A * x_true + 0.02 * (randn(m,1) + 1i*randn(m,1)); % 加噪

% --- 调用 FoBa-AN ---
lambda = 0.1;               % 正则化参数
[x_est, S, w] = foba_an_line_spectrum(A, b, lambda, 10, 1e-3, 1e-6);

% --- 结果展示 ---
fprintf('\n真实频率: '); disp(true_freqs);
fprintf('恢复频率: '); disp(S(:)');

figure;
subplot(2,1,1);
stem(true_freqs, abs(true_coeffs), 'r', 'filled', 'DisplayName','真实');
hold on;
stem(S, abs(w), 'b', 'o', 'DisplayName','恢复');
xlabel('归一化频率'); ylabel('幅度');
legend; title('频率估计');
grid on;

subplot(2,1,2);
plot(1:n, real(x_true), 'r', 'LineWidth',1.5); hold on;
plot(1:n, real(x_est), 'b--', 'LineWidth',1.2);
xlabel('样本序号'); ylabel('实部');
legend('真实信号','恢复信号');
title('时域信号恢复');
grid on;