基于MATLAB的快速迭代收缩阈值算法（FISTA）实现

一、FISTA算法原理

FISTA通过引入动量项将ISTA的收敛速度从O(1/k)提升至O(1/k²)，其核心迭代公式为：

其中：

• $S_{λt}$为软阈值算子

• 为动态步长

二、MATLAB实现代码

1. 核心函数实现

function [x, out] = FISTA(y, A, lambda, maxIter, tol)
    % 参数说明：
    % y: 观测信号 (M×1)
    % A: 测量矩阵 (M×N)
    % lambda: 正则化参数
    % maxIter: 最大迭代次数
    % tol: 收敛阈值

    [M, N] = size(A);
    x = zeros(N,1);  % 初始化解
    y_aux = x;       % 辅助变量
    t = 1;           % 初始步长

    for k = 1:maxIter
        % 计算梯度
        grad = (1/M) * (A' * (A * y_aux - y)); 

        % 动量更新
        x_new = softThreshold(y_aux - t*grad, lambda*t); 

        % 动态步长更新
        t_new = (1 + sqrt(1 + 4*t^2)) / 2; 

        % 动量加速
        y_aux = x_new + ((t - 1)/t_new) * (x_new - x); 

        % 收敛判断
        if norm(x_new - x) < tol
            break;
        end
        x = x_new;
        t = t_new;
    end

    out.iter = k;
    out.fval = 0.5*norm(A*x - y)^2 + lambda*norm(x,1);
end

function x = softThreshold(x, kappa)
    % 软阈值操作
    x = sign(x) .* max(abs(x) - kappa, 0);
end

2. 参数设置与调用示例

%% 参数设置
lambda = 0.1;     % 正则化强度
maxIter = 1000;   % 最大迭代次数
tol = 1e-6;       % 收敛阈值

%% 数据生成
n = 100;          % 信号长度
k = 10;           % 稀疏度
A = randn(50,n);  % 测量矩阵
x_true = zeros(n,1);
x_true(randperm(n,k)) = randn(k,1);
y = A*x_true + 0.1*randn(50,1);  % 含噪声观测

%% 算法运行
tic;
[x_recon, out] = FISTA(y, A, lambda, maxIter, tol);
time_cost = toc;

%% 结果展示
figure;
subplot(2,1,1);
stem(x_true,'b','LineWidth',1.5); hold on;
stem(x_recon,'r--','LineWidth',1.5);
legend('真实信号','重构信号');
title(sprintf('FISTA重构结果 (迭代次数=%d, 耗时=%.2fs)', out.iter, time_cost));

subplot(2,1,2);
semilogy(1:out.iter, out.fval, 'r-o');
xlabel('迭代次数'); ylabel('目标函数值');
title('收敛曲线');

三、关键优化

1. 自适应步长选择（BB步长法）

function t = bb_step(A, x_prev, x_curr, y_prev, y_curr)
    d = x_curr - x_prev;
    g = y_curr - y_prev;
    t = (d'*d)/(d'*g);
end

% 在FISTA主循环中替换原有步长计算：
t = bb_step(A, y_aux, x_new, grad_prev, grad);

2. 并行计算加速

% 使用parfor加速梯度计算（适用于大规模矩阵）
parfor i = 1:M
    grad(i) = A(i,:) * (A * y_aux - y);
end
grad = grad / M;

3. 动态正则化参数

% 迭代衰减策略
lambda = lambda0 * (1 - exp(-iter/100));

四、实验对比（vs ISTA）

五、应用场景扩展

1. 图像去噪

% 图像去噪实现
denoised_img = FISTA(noisy_img(:), measurement_matrix, lambda, 1000, 1e-6);
denoised_img = reshape(denoised_img, [256,256]);

2. 动态MRI重建

% 动态MRI加速重建
for frame = 1:num_frames
    y = measure_frame(frame);
    x(:,:,frame) = FISTA(y, A, lambda, 500, 1e-5);
end

参考代码 MATLAB快速迭代收缩阈值算法 www.youwenfan.com/contentalh/59575.html

六、常见问题处理

1. 不收敛问题：

   • 检查梯度计算准确性
   • 调整步长范围（建议初始t=1，范围[1e-3,1e3]）
   • 增加最大迭代次数

2. 计算效率优化：

   • 使用GPU加速（需安装Parallel Computing Toolbox）

   gpu_A = gpuArray(A);
   gpu_y = gpuArray(y);
   x_gpu = FISTA(gpu_y, gpu_A, lambda, maxIter, tol);
   x = gather(x_gpu);
   

3. 稀疏性控制：

   • 调整λ值（λ越大稀疏性越强）
   • 结合OMP算法进行二次稀疏优化

七、参考文献

1. Beck A, Teboulle M. Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems[J]. SIAM Journal on Imaging Sciences, 2009.
2. Nesterov Y. A method for solving the convex programming problem with convergence rate O(1/k²)[C]. Dokl. Akad. Nauk SSSR, 1983.
3. 深度学习中的优化算法实践（李沐，2021）