快速迭代收缩阈值算法 FISTA 原理与实现

如此迭代产生序列 x_k 收敛于全局最优解。当 f 为 Lipschitz 光滑凸函数且 g 为凸函数时，ISTA 可保证 F(x_k) 以 O(1/k) 的速率逼近最优值。然而，ISTA 在实践中往往收敛较慢，需要大量迭代。

快速迭代收缩 - 阈值算法 (FISTA) 是 Beck 和 Teboulle (2009) 提出的对 ISTA 的加速改进。它借鉴了 Nesterov 在 1983 年提出的加速梯度思想，引入了动量项来加速收敛。具体来说，FISTA 并不直接对上一次的解 x_{k-1} 进行近端梯度，而是构造一个加速序列 y_k，使其包含了前两次迭代解的线性组合信息。FISTA 的迭代步骤如下：

1. 初始化： 取初始解 x_0（可取零向量），设 y_1 = x_0，并设动量参数 t_1 = 1。
2. 循环迭代： 对于 k = 1, 2, 3, ...:
   • (近端梯度步) 先对当前辅助变量 y_k 执行与 ISTA 相同的近端梯度更新：x_k = prox_{1/L g}(y_k - 1/L grad f(y_k))。对于 L1 情况，这一步仍是 x_k = S_{lambda/L}(y_k - 1/L grad f(y_k))。
   • (动量更新步) 计算新的动量系数：t_{k+1} = (1 + sqrt(1 + 4*t_k^2)) / 2。然后更新加速序列：y_{k+1} = x_k + (t_k - 1)/t_{k+1} * (x_k - x_{k-1})。这一步利用了当前解 x_k 和前一步解 x_{k-1} 的线性组合来得到新的搜索点 y_{k+1}。

可以看出，与 ISTA 每次仅利用 x_{k-1} 来迭代不同，FISTA 在 y_{k+1} 中加入了 (x_k - x_{k-1}) 项，也就是利用最近两次迭代的解的差值作为'动量'。这种设计看似'魔法般', 实则源自对迭代过程几何性质的精巧分析。通过恰当选择 t_k 序列（t_k 的递推形式来源于求解二次方程，保证后续分析中一个递推不等式的最优收敛速率），FISTA 能在理论上将收敛速度提升一个数量级。

FISTA 保持了与 ISTA 相近的每步计算复杂度，但其全局收敛速率提高到 O(1/k^2)。也就是说，达到同样精度所需的迭代次数大大减少。具体定理表明：对于 k >= 1，F(x_k) - F(x*) <= 2L ||x_0 - x*||^2 / (k+1)^2，其中 x* 为最优解，L 为 grad f 的 Lipschitz 常数。相较之下，ISTA 只有 O(1/k) 的渐进速率。因此在大量迭代后，FISTA 会显著优于 ISTA。

引入的辅助序列 y_k 可以看作在常规梯度下降基础上添加了一种'加速动量'策略，类似于物理中的动量累积，使算法在谷底附近不至于过于保守地缓慢逼近，而是带着惯性前进，从而快速逼近最优解。这与 Nesterov 加速梯度法在纯光滑优化中的作用类似。

下面用伪代码总结 FISTA 的核心步骤：

% 给定 A, b, 正则权重 lambda, Lipschitz 常数 L, 初始化 x0
x = x0;
y = x0;
t = 1;
for k = 1:MaxIter
    % 近端梯度步：梯度下降 + 软阈值
    grad = A'*(A*y - b); % 计算 grad f(y) = A^T(Ay - b)
    x_new = sign(y - grad/L) .* max(abs(y - grad/L) - lambda/L, 0);
    % 动量更新步
    t_new = (1 + sqrt(1 + 4*t^2)) / 2;
    y = x_new + ((t - 1) / t_new) * (x_new - x);
    % 更新迭代变量
    x = x_new;
    t = t_new;
end

上述代码实现了不带回溯的 FISTA 算法框架：首先计算基于当前 y 的梯度下降点，然后对其应用软阈值（对应 prox_{lambda/L} 操作），接着根据动量公式更新 y。经过多次迭代，x 将收敛到稀疏解。值得注意的是，该算法需要估计合理的 L（例如，可取 L = |A^T A| 的最大特征值）。若 L 不易获得，也可以在每步迭代中采用回溯线搜索调整步长，这也是 FISTA 论文中提到的变体，但基本思想不变。

FISTA 的详尽收敛性证明较为复杂，涉及构造恰当的辅助序列和能量函数。其核心在于证明 y_k 的选取使得 F(x_k) - F(x*) 符合所需的 O(1/k^2) 上界。这里不赘述完整推导，但强调两点：(1) t_k 的递推形式 (1 + sqrt(1 + 4*t_k^2))/2 是精心设计的，使得动量项权重逐渐减小，避免过冲；(2) FISTA 的形式也可从将 Nesterov 加速应用于 ISTA 的等价形式推导得到，一些后续研究表明 FISTA 可以被视为对近端梯度法在最坏情况下加速的'最优'方法。

3. Matlab 实现

假设我们要解决一个简单的压缩感知重建问题：Ax ≈ b，其中 A 是测量矩阵，b 是观测到的信号，我们希望通过最小化 F(x) = 1/2 ||Ax - b||_2^2 + lambda ||x||_1 来求得稀疏解 x。FISTA 可高效地求解此问题。

Matlab 函数构造：

function x = FISTA_L1(A, b, lambda, L, maxIter)
    % 初始化
    x = zeros(size(A,2), 1); % 初始解 x0
    y = x;
    t = 1;
    for k = 1:maxIter
        % 计算梯度并执行梯度下降一步
        grad = A' * (A * y - b);
        z = y - (1/L) * grad;
        % 执行 L1 近端（软阈值）操作
        x_new = sign(z) .* max(abs(z) - lambda/L, 0);
        % 更新动量系数和辅助变量
        t_new = (1 + sqrt(1 + 4 * t^2)) / 2;
        y = x_new + ((t - 1) / t_new) * (x_new - x);
        % 准备下次迭代
        x = x_new;
        t = t_new;
    end
end

上述函数 FISTA_L1 实现了 FISTA 的主要迭代流程。其中：

• 初始化部分，将初始解设为全零向量，并初始化辅助变量 y = x 及动量参数 t = 1。
• 每次迭代首先计算当前 y 下的梯度 grad，对应 grad f(y) = A^T(Ay - b)。
• 然后令 z = y - (1/L)*grad，这一步是普通梯度下降的结果。
• 接着对 z 执行软阈值：x_new = sign(z).*max(abs(z)-lambda/L, 0)，得到当前迭代的稀疏估计 x_k。
• 计算新的动量系数 t_new 并更新辅助变量：y = x_new + ((t-1)/t_new)*(x_new - x)，对应 y_{k+1} = x_k + (t_k-1)/t_{k+1}*(x_k - x_{k-1})。
• 循环更新直到达到最大迭代次数（或其它收敛判据，如迭代前后 |x| 相对变化很小）。

使用该函数时，需要提供合理的 L（一般可取为 A^TA 的最大特征值近似）。这里我们可以测试一个小例子：生成随机高斯矩阵 A 以及稀疏的真实信号 x_true，计算 b = Ax_true（可加噪声），然后调用 FISTA_L1(A,b,lambda,L,maxIter) 重建，并与真值比较。

%% FISTA_L1_demo.m
% 使用 FISTA + L1 进行稀疏信号重建。
clear; clc; close all;
%% 准备数据
rng('default'); % 固定随机种子，方便复现
m = 80; % 测量数
n = 100; % 信号长度
A = randn(m, n); % 随机生成测量矩阵 A
% 生成一个稀疏的 x_true
k0 = 10; % 希望的非零项数
permIdx = randperm(n); % 随机打乱下标
supp = permIdx(1:k0);
x_true = zeros(n,1);
x_true(supp)= randn(k0,1); % 在这 k0 个位置上放一些随机值
% 生成观测 b
b = A * x_true;
%% 设置 FISTA 参数
lambda = 0.1; % L1 惩罚系数
[U, S, ~] = svd(A, 'econ');
L = max(diag(S))^2; % 步长倒数，近似取 A 的最大奇异值平方
maxIter = 100; % 最大迭代次数
%% 调用 FISTA_L1 求解
x_est = FISTA_L1(A, b, lambda, L, maxIter);
%% 比较结果
recovery_error = norm(x_est - x_true) / norm(x_true);
fprintf('重建相对误差 = %.4e\n', recovery_error);
%% 作图比较 x_true 与 x_est
figure('Name', 'FISTA-L1 Reconstruction', 'NumberTitle', 'off');
plot(1:n, x_true, '-o', 'LineWidth', 1.5);
hold on;
plot(1:n, x_est, '--', 'LineWidth', 1.5);
grid on;
legend('x_true', 'x_est', 'Location', 'best');
xlabel('Index (n)');
ylabel('Signal Value');
title('FISTA + L1 Sparse Reconstruction');
% 设置坐标轴字体
set(gca, 'FontName', 'Times New Roman', 'FontSize', 13);

用到的函数如下

function x = FISTA_L1(A, b, lambda, L, maxIter)
    % FISTA_L1 使用 FISTA 迭代来求解
    % min_x 0.5*||A*x - b||^2 + lambda * ||x||_1
    %
    % 输入：
    % A, b 线性系统和观测
    % lambda L1 惩罚系数
    % L 步长因子，通常可取近似的最大特征值 A'*A
    % maxIter 最大迭代次数
    %
    % 输出：
    % x 稀疏重建后的解向量
    % 初始化
    x = zeros(size(A,2),1); % 初始解 x0
    y = x; % 初始化辅助变量 y
    t = 1; % 动量参数 t0
    for k = 1:maxIter
        % 计算梯度
        grad = A'*(A*y - b);
        grad = A' * (A*y - b);
        % 常规梯度下降一步：
        z = y - (1/L)*grad;
        z = y - (1/L)*grad;
        % 执行软阈值(shrinkage)操作，得到 x_new
        x_new = sign(z) .* max(abs(z) - lambda/L, 0);
        % 更新动量参数 t_new
        t_new = (1 + sqrt(1 + 4*t^2)) / 2;
        % 更新辅助变量 y
        y = x_new + ((t - 1)/t_new) * (x_new - x);
        % 准备下一次迭代
        x = x_new;
        t = t_new;
    end
end

4. FISTA 在图像处理与压缩感知中的应用

FISTA 自提出以来已在图像处理领域的诸多逆问题中显示出强大能力，包括去去噪、欠采样重建等。其优势在于算法简单、易于实现，却又兼具接近于最优的一阶收敛速度，使其非常适合求解大规模问题。下面结合具体应用谈几点：

4.1. 基于小波稀疏先验的图像去噪

我们将把去噪问题写成 min_x 1/2 ||x - y||_2^2 + lambda ||W x||_1 其中

• y 是加噪图像
• W 是二维小波变换算子
• lambda 是稀疏惩罚参数

对于这个问题，

• f(x) = 1/2 ||x-y||^2，梯度 grad f(x) = x-y，Lipschitz 常数 L = 1。
• g(x) = lambda ||W x||_1，其近端算子就是对小波系数做软阈值。

实现上述去噪问题的主函数为

function x = fista_wavelet_denoise(y, lambda, waveletName, levels, maxIter)
    %FISTA_WAVELET_DENOISE 用 FISTA+小波稀疏做图像去噪
    % 输入：
    % y - 加噪图像（二维灰度，double）
    % lambda - 惩罚系数
    % waveletName - 小波类型，如 'db4'
    % levels - 分解层数
    % maxIter - 最大迭代次数
    % 输出：
    % x - 去噪后图像
    % 参数检查
    if nargin<5, maxIter = 200; end
    % 初始化
    xk = y;
    yk = xk;
    tk = 1;
    L = 1; % f 的 Lipschitz 常数
    % 预分解一次（仅获取尺寸）
    [C0, S] = wavedec2(yk, levels, waveletName);
    for k = 1:maxIter
        % —— 1. 梯度步 ——
        grad = yk - y; % grad f(yk)
        zk = yk - (1/L)*grad; % 梯度下降
        % —— 2. 近端步 —— 小波域软阈值 ——
        C = wavedec2(zk, levels, waveletName);
        Cthr = soft_thresh(C, lambda/L);
        x_next = waverec2(Cthr, S, waveletName);
        % —— 3. FISTA 加速步 ——
        t_next = (1 + sqrt(1 + 4*tk^2))/2;
        y_next = x_next + ((tk-1)/t_next)*(x_next - xk);
        % 更新
        xk = x_next;
        yk = y_next;
        tk = t_next;
    end
    x = xk;
end

上面用到的软阈值函数为：

function z = soft_thresh(v, thresh)
    %SOFT_THRESH 对向量 v 执行元素级软阈值
    z = sign(v) .* max(abs(v) - thresh, 0);
end

然后我们就可以测试这个去噪程序了，如下

%% demo_fista_denoise.m
clc; clear; close all;
% 1. 读取并归一化
I = im2double(imread('cameraman.tif'));
% 2. 加高斯噪声
sigma = 0.05;
yn = I + sigma*randn(size(I));
% 3. 参数设置
lambda = 0.1; % 惩罚强度，可根据噪声调节
waveletName = 'db4'; % Daubechies4 小波
levels = 3; % 分解层数
maxIter = 100; % 迭代次数
% 4. 运行 FISTA 去噪
xd = fista_wavelet_denoise(yn, lambda, waveletName, levels, maxIter);
% 5. 显示对比
figure('Position',[100 100 800 300]);
subplot(1,3,1);
imshow(I);
title('原始图像');
subplot(1,3,2);
imshow(yn);
title(['加噪图像 (σ=',num2str(sigma),')']);
subplot(1,3,3);
imshow(xd);
title('FISTA+小波 去噪结果');

4.2 压缩感知图像重建

在压缩感知(Compressed Sensing, CS) 框架中，我们以远低于奈奎斯特采样率的线性测量来重建图像。经典 CS 理论指出，当图像在某种域（如小波）是稀疏的时，可以通过求解 L1 正则化的凸优化精确重建图像。FISTA 正是求解这类问题的理想算法之一。举例来说，在 MRI 成像中采集不完备 k 空间数据时，常建立 min_x 1/2 ||Ax-b||^2 + lambda ||W(x)||_1（W(x) 表示小波变换系数）模型来重建图像；实践表明利用 FISTA 求解该模型能够在保持高重建质量的同时大幅减少迭代时间。即使在深度学习蓬勃发展的今天，研究者仍关注传统 CS 重建的潜力，通过合理优化，经典 L1-范数重建在某些情况下仍然有竞争力。这也证明了 FISTA 等算法在成像领域的持久价值。

我们以最简单的'稀疏先验 + 小波变换'为例，解决 min_x 1/2 ||Phi x - y||_2^2 + lambda ||W x||_1, 其中

• x ∈ R^N 是待重建的灰度图像（展开成向量），
• Phi ∈ R^{M×N} 是测量矩阵（M << N），
• y ∈ R^M 是测量向量，
• W 和 W^{-1} 分别是二维小波变换算子和其逆算子。

实现上述稀疏重建的主函数为

function x_rec = fista_cs_wavelet(Phi, y, lambda, waveletName, levels, imgSize, opts)
    %FISTA_CS_WAVELET 用 FISTA+小波稀疏做压缩传感重建
    % 输入：
    % Phi - M×N 测量矩阵
    % y - M×1 测量向量
    % lambda - 稀疏惩罚系数
    % waveletName - 小波类型，如 'db4'
    % levels - 小波分解层数
    % imgSize - [h, w] 原图高宽
    % opts - 可选结构体
    % .maxIter (默认 200)
    % .tol (默认 1e-4)
    % 输出：
    % x_rec - 重建后图像（h×w）
    if nargin<7, opts = struct(); end
    if ~isfield(opts,'maxIter'), opts.maxIter = 200; end
    if ~isfield(opts,'tol'), opts.tol = 1e-4; end
    % 初始化
    N = prod(imgSize);
    xk = zeros(N,1);
    yk = xk;
    tk = 1;
    L = norm(Phi,2)^2; % Lipschitz 常数
    % 预计算小波分解结构
    [~, S] = wavedec2(reshape(yk,imgSize), levels, waveletName);
    for k = 1:opts.maxIter
        % 1) 梯度下降
        grad = Phi'*(Phi*yk - y);
        z = yk - (1/L)*grad;
        % 2) 小波软阈值
        Cz = wavedec2(reshape(z,imgSize), levels, waveletName);
        Cthr = soft_thresh(Cz, lambda/L);
        x_next = waverec2(Cthr, S, waveletName);
        x_next = x_next(:);
        % 3) FISTA 加速
        t_next = (1 + sqrt(1 + 4*tk^2)) / 2;
        y_next = x_next + ((tk - 1)/t_next)*(x_next - xk);
        % 收敛判断
        if norm(x_next - xk) < opts.tol
            xk = x_next;
            break;
        end
        % 更新
        xk = x_next;
        yk = y_next;
        tk = t_next;
    end
    % 输出重建图像
    x_rec = reshape(xk, imgSize);
end

然后我们就可以用于演示一个采样率为 50% 的图像重建，

% demo_fista_cs.m
clc; clear; close all;
% 1. 原始图像及展平
I = im2double(imread('cameraman.tif'));
I = imresize(I, 0.25);
imgSize = size(I);
x_true = I(:);
% 2. 构造测量矩阵 Phi（随机高斯或 0/1 伯努利）
M = round(0.5 * numel(x_true)); % 50% 采样率
Phi = randn(M, numel(x_true));
Phi = bsxfun(@rdivide, Phi, sqrt(sum(Phi.^2,2))); % 每行归一化
% 3. 生成测量 y
noise_sigma = 1e-3;
y = Phi * x_true + noise_sigma*randn(M,1);
% 4. 重建参数
lambda = 0.01;
waveletName = 'db4';
levels = 3;
opts.maxIter= 300;
opts.tol = 1e-5;
% 5. 调用 FISTA-CS 重建
tic;
x_rec = fista_cs_wavelet(Phi, y, lambda, waveletName, levels, imgSize, opts);
toc;
% 6. 显示结果
figure('Position',[100 100 900 300]);
subplot(1,2,1);
imshow(I);
title('原始图像');
subplot(1,2,2);
imshow(x_rec,[]);
title('FISTA-CS 重建');