快速迭代收缩-阈值算法（FISTA）

1. 数学与优化基础

在许多信号处理与机器学习问题中，我们希望获得稀疏解，即解向量中非零元素很少。直接以 L0 “范数” （非零元素个数）作为稀疏性度量会导致非凸优化难题。一个常用的替代是采用 L1 范数（||x||₁是各元素绝对值之和）作为正则项，它是 L0的凸松弛，可以鼓励稀疏解。与 L2 正则化不同，L1 正则化对小系数给予恒定的惩罚减少，这使得较小的系数更容易被压缩到0，从而产生稀疏性。简言之，L1 范数由于在零点处的尖锐拐角（梯度不连续）会诱导解的许多分量为零。

假设你有 5 个箱子，里面分别放了 $x_1,x_2,x_3,x_4,x_5$ 个苹果。L0“范数”就是这 5 个箱子里有苹果的箱子的数量。

• 若 $x_2=0,x_4=0$，其它 3 个箱子里苹果不为 0，那么 L0 “范数”就是 3。
• 想要让 L0 变小，就要“尽量让更多箱子变成空箱子”，从而达到稀疏(空箱子越多，非零元素越少)。
• L0“范数”本身是不连续、不光滑的(要么 0，要么 1)，在数学优化中很难直接求解(优化过程不易“沿着梯度”慢慢变化，一下跳零一下跳非零)。这会导致“最小化 L0”时比较棘手。

L1 范数是一个凸函数，有很好处理的数学性质(可以使用梯度或次梯度等优化算法)。虽然它并不会像 L0 一样“明确地”只算非零数量，但在最小化 L1 范数时，有很多解会自动出现“某些元素被压到正好为 0”——这就是所谓的“稀疏解”。依然用“箱子”来打比方。

• L1 范数就是把每个箱子里的苹果数量(绝对值)加起来之和，努力让这个总和尽量小。
• 当某个箱子里的苹果数量比较小，L1 优化就会倾向于“干脆把这个箱子的苹果数量压到 0”，以进一步降低总和。
• 于是就会出现类似 L0 的“有些箱子空了”的效果(稀疏解)。

L2 范数是各元素平方和的开方，即 $\Vert x{\Vert }_2=\sqrt{x_1^2+x_2^2+\cdots +x_n^2}$。它可以让解趋向于“总体较小”，但不会把任何一维直接压到完全 0。若还以箱子比方，L2 更像是在“平均地”减小每个箱子的苹果数量，而不会直接把其中某些箱子砍到空箱。

典型例子包括LASSO（最小绝对收缩选择算子）和**基础追踪（Basis Pursuit）**等，它们通过在最小二乘误差目标后加上 L1 范数惩罚来实现稀疏信号的估计。这种 L1 正则化问题属于凸优化，通常可转化为线性规划或二阶锥规划来求解。例如，有噪声情况下的基础追踪去噪（BPDN）和 LASSO 都考虑求解形如
$$\underset{x}{\min }\frac{1}{2}\Vert Ax-b{\Vert }_2^2+\lambda \Vert x{\Vert }_1,$$
的凸问题，其中 λ 为权衡稀疏度的正则化参数。

像上面这样的目标函数 $F(x)=f(x)+g(x)$ 中，$f(x)=\frac{1}{2}\Vert Ax-b{\Vert }_2^2$ 是光滑凸函数，具有Lipschitz连续的梯度，而 $g(x)=\lambda \Vert x{\Vert }_1$ 是连续凸函数但在0点不可微。

凸函数的一个重要性质是局部极小即全局极小，这为设计算法提供了可靠性保障。

梯度下降法是求解光滑凸优化的基本一阶方法：反复沿负梯度方向迭代。对于二次损失这样的$f(x)$，梯度为$\nabla f(x)=A^T(Ax-b)$，直接的梯度迭代为 $x\leftarrow x-\eta \nabla f(x)$。然而，当存在非光滑项$g(x)$时（例如 L1 正则化），我们不能简单对其求梯度，而需借助次梯度或近端算法。

近端梯度（Proximal Gradient）方法通过引入近端算子来处理非光滑项，每步对光滑部分做梯度更新后，对非光滑部分应用一个近端映射(即求解一个带L2范数惩罚的最优问题)。对于 L1 项，其近端映射就是**软阈值（soft-thresholding）**操作，其公式为：
$${\mathrm{prox}}_{\lambda \Vert \cdot {\Vert }_1}(v_i)=\{\begin{array}{ll}0,& |v_i|\le \lambda ,\\ v_i-\lambda \mathrm{sign}(v_i),& |v_i|>\lambda ,\end{array}$$
即将向量 $v$ 中每个分量按上述规则收缩。软阈值操作也通常写为 $S_{\lambda }(v)=\mathrm{sign}(v)\cdot \max (|v|-\lambda ,0)$，把绝对值低于阈值$\lambda$的系数设为0，其余系数减小$\lambda$。软阈值是 L1 非光滑项的近端，因为它正是最小化 $\frac{1}{2}\Vert x-v{\Vert }_2^2+\lambda \Vert x{\Vert }_1$ 所得的解。这个操作在稀疏信号重建中可视为一种“去噪”步骤：减小信号中的小幅成分，同时保留/突出大幅成分，从而促进稀疏性。

2. FISTA 算法的原理、推导与机制

在引入 FISTA 前，先看其基础算法——迭代收缩阈值算法 (ISTA)。ISTA 是一种近端梯度下降方法，每次迭代包括一个梯度下降步骤和一个近端（软阈值）步骤。对目标 $F(x)=f(x)+g(x)$，给定梯度Lipschitz常数为 $L$，ISTA 的每步更新可表示为：

• 梯度步： $y_k=x_{k-1}-\frac{1}{L}\nabla f(x_{k-1})$，即从上一步解 $x_{k-1}$ 沿负梯度方向前进一步长度 $1/L$。
• 近端步： $x_k={\mathrm{prox}}_{\frac{1}{L}g}(y_k)$，即对 $y_k$ 应用非光滑项 $g$ 的近端算子。如果 $g(x)=\lambda \Vert x{\Vert }_1$，这一步就是对 $y_k$ 进行软阈值 $S_{\lambda /L}(y_k)$。

如此迭代产生序列 $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,\dots$:
   • (近端梯度步) 先对当前辅助变量 $y_k$ 执行与 ISTA 相同的近端梯度更新： $x_k={\mathrm{prox}}_{\frac{1}{L}g}(y_k-\frac{1}{L}\nabla f(y_k)).$ 对于 L1 情况，这一步仍是 $x_k=S_{\lambda /L}(y_k-\frac{1}{L}\nabla f(y_k))$。
   • (动量更新步) 计算新的动量系数： $t_{k+1}=\frac{1+\sqrt{1+4t_k^2}}{2}$。然后更新加速序列： $y_{k+1}=x_k+\frac{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\ge 1$， $F(x_k)-F(x^{\ast })\le \frac{2L\Vert x_0-x^{\ast }{\Vert }^2}{(k+1{)}^2},$ 其中 $x^{\ast }$为最优解，$L$为$\nabla 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);                  % 计算∇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$ 的梯度下降点，然后对其应用软阈值（对应 ${\mathrm{prox}}_{\lambda /L}$ 操作），接着根据动量公式更新 $y$。经过多次迭代，$x$ 将收敛到稀疏解。值得注意的是，该算法需要估计合理的 $L$（例如，可取 $L=|A^{T}A|$ 的最大特征值）。若 $L$ 不易获得，也可以在每步迭代中采用回溯线搜索调整步长，这也是 FISTA 论文中提到的变体，但基本思想不变。

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

3. Matlab 实现

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

Matlab 函数构造：

function x = FISTA_L1(A, b, lambda, L, maxIter)% 初始化x = zeros(size(A,2), 1);    % 初始解 x0y = 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，对应 $\nabla 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+\frac{t_k-1}{t_{k+1}}(x_k-x_{k-1})$。
• 循环更新直到达到最大迭代次数（或其它收敛判据，如迭代前后$|x|$相对变化很小）。

使用该函数时，需要提供合理的 $L$（一般可取为 $A^{T}A$ 的最大特征值近似）。这里我们可以测试一个小例子：生成随机高斯矩阵 $A$ 以及稀疏的真实信号 $x_{\text{true}}$，计算 $b=A{x}_{\text{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);  % 初始解 x0y = x;                   % 初始化辅助变量 yt = 1;                   % 动量参数 t0for k = 1:maxIter% 计算梯度 grad = A'*(A*y - b)grad = A' * (A*y - b);% 常规梯度下降一步： z = y - (1/L)*gradz = y - (1/L)*grad;% 执行软阈值(shrinkage)操作，得到 x_newx_new = sign(z) .* max(abs(z) - lambda/L, 0);% 更新动量参数 t_newt_new = (1 + sqrt(1 + 4*t^2)) / 2;% 更新辅助变量 yy = x_new + ((t - 1)/t_new) * (x_new - x);% 准备下一次迭代x = x_new;t = t_new;end
end

结果如下：

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

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

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

我们将把去噪问题写成
$$\underset{x}{\min }\frac{1}{2}\Vert x-y{\Vert }_2^2+\lambda \Vert Wx{\Vert }_1$$
其中

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

对于这个问题，

• $f(x)=\frac{1}{2}\Vert x-y{\Vert }^2$，梯度 $\nabla f(x)=x-y$，Lipschitz 常数 $L=1$。
• $g(x)=\lambda \Vert Wx{\Vert }_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;          % ∇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;endx = 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\frac{1}{2}\Vert Ax-b{\Vert }^2+\lambda \Vert W(x){\Vert }_1$（$W(x)$表示小波变换系数）模型来重建图像；实践表明利用FISTA求解该模型能够在保持高重建质量的同时大幅减少迭代时间。即使在深度学习蓬勃发展的今天，研究者仍关注传统CS重建的潜力，通过合理优化，经典L1-范数重建在某些情况下仍然有竞争力。这也证明了FISTA等算法在成像领域的持久价值。

我们以最简单的“稀疏先验+小波变换”为例，解决
$$\underset{x}{\min }\frac{1}{2}\Vert \Phi x-y{\Vert }_2^2+\lambda \Vert Wx{\Vert }_1,$$
其中

• $x\in {\mathbb{R}}^N$ 是待重建的灰度图像（展开成向量），
• $\Phi \in {\mathbb{R}}^{M\times N}$ 是测量矩阵（$M\ll N$），
• $y\in {\mathbb{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(); endif ~isfield(opts,'maxIter'), opts.maxIter = 200; endif ~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.tolxk = 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. 构造测量矩阵 Φ（随机高斯或 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 重建');