ADMM,ISTA,FISTA算法步骤详解，MATLAB代码，求解LASSO优化问题

ADMM,ISTA,FISTA算法步骤详解，MATLAB代码，求解LASSO优化问题

原创文章！转载需注明来源：©️ Sylvan Ding’s Blog ❤️

实验目的

• 了解 ADMM, ISTA, FISTA 算法的基本原理、收敛性和复杂度；
• 使用上述三种算法，解决 LASSO 问题；
• 分析三种算法的表现情况。

实验环境

MATLAB R2021a

实验内容

ADMM 和 [F]ISTA 的收敛性证明在 [1] 和 [2] 中给出，ADMM 和 ISTA 的全局收敛速率（global convergence rate）为 $O(1/k)$ [Eckstein and Bertsekas, 1990; Deng and Yin, 2012; He and Yuan, 2012]，然而 FISTA 的全局收敛速率为 $O(1/k^2)$​ 。但是，[2] 得出了在某些情况下，ISTA 却比 FISTA 更快的结论，并给出了具体原因。

LASSO Problem

许多机器学习和数据拟合问题都可以看成是最小二乘问题，添加正则项从而避免过拟合。在许多应用中，使用 $l1$-范数 作为正则项可以带来不错的泛化效果，所以我们求解含有 $l1$-正则项 的线性最小二乘问题，也称为 LASSO 问题 [2]，其一般形式为：（P）

$$\underset{x\in {\mathbb{R}}^n}{\min }\frac{1}{2}\Vert Ax-b{\Vert }^2+\lambda \Vert x{\Vert }_1.$$

其中，$A\in {\mathbb{R}}^{m\times n}$ 是一个行满秩矩阵，$n>m$ ；$b$ 是一个已知向量； $\lambda >0$ 是一标量；$l1$-正则项 $\Vert x{\Vert }_1$ 会产生一稀疏解，在降低代价的同时，避免发生过拟合 [Tibshirani, 1996]。

在实验一和实验三中，分别使用了不同的正则项实现了信号的降噪，体现了 $l1$-正则项 的另一个优势，就是相较于 $l2$-正则项 来说，$l1$-正则项 对边界值（outliers）不敏感，这一特点有利于图像降噪中对锐利边缘（sharp edges）的处理。

在实验三中，曾将（P）视为盒状约束（box-constrained）的二次优化问题，使用梯度投影法进行求解。本文要求使用 ADMM, ISTA, FISTA 这三种算法，求解 LASSO 问题。

ADMM

ADMM 算法（Alternating Direction Method of Multipliers）[Gabay, Mercier, Glowinski, Marrocco, 1976] 将原问题的变量 $x$ 分解为两个变量 $x$ 和 $z$ ，在最小二乘的损失函数里使用 $x$ ，在 $l1$-范数 的正则项里使用 $z$ 。这样，外加一个等式约束，就可构造增广拉格朗日函数，进而分别求解两个变量的最优值，使得计算较为容易。

问题（P）等价于：

$$\begin{array}{lc}\underset{x}{\min }& \frac{1}{2}\Vert Ax-b{\Vert }^2+\lambda \Vert z{\Vert }_1\\ s.t.& x-z=0\end{array}$$

构造增广拉格朗日函数（augmented Lagrangian） $L_{\rho }$ ：

$$L_{\rho }(x,z,\mu )=f(x)+g(z)+{\mu }^T(x-z)+\frac{\rho }{2}\Vert x-z{\Vert }^2.$$

其中，$f(x)=\frac{1}{2}\Vert Ax-b{\Vert }^2$ ，$g(z)=\lambda \Vert z{\Vert }_1$ ，$\rho$ 是惩罚参数。

不难发现，最优化 $x$ 时，是一个二次优化问题，而最优化 $z$ 时，是一个软阈（soft-thresholding）问题，那么，ADMM 的迭代形式为：

$$\begin{array}{lcc}{x}^{k+1}& :=(A^{T}A+\rho I{)}^{-1}(A^{T}b+\rho z^k-{\mu }^k)& x.minimization\\ z^{k+1}& :=S_{\lambda /\rho }(x^{k+1}+{\mu }^k/\rho )& z.minimization\\ {\mu }^{k+1}& :={\mu }^k+\rho (x^{k+1}-z^{k+1})& \mathrm{dual}\mathrm{update}\end{array}$$

其中，$S_{\tau }$ 是软阈函数（又称为Shrinkage），$S_{\tau }(g):=sign(g)\cdot (|g|-\tau {)}_{+}$ [3]，也即 [5]

$${\left[S_{\tau }(g)\right]}_j=\{\begin{array}{cc}{g}_j-\tau & g>\tau \\ 0& -\tau \le g\le \tau \\ g_j+\tau & g<-\tau \end{array},j=1,2,\dots ,n$$

实际上，ADMM 算法往往在一系列迭代后，能获得较为精确的解，但是所需要的迭代次数是大量的 [5]。

$\rho$ 的选取会在极大程度上影响 ADMM 的收敛性。若 $\rho$ 太大，则对于优化 $f+g$ 的能力会下降；反之，则会削弱约束条件 $x=z$. Boyd et al. (2010) 给出了选 $\rho$ 的策略，效果不错，但不能保证一定收敛。

ISTA

现在流行的一种解决问题（P）的方法是 ISTA（iterative shrinkage-thresholding algorithms）[Daubechies et al, 2004]，该方法在每次迭代时会计算矩阵和向量的乘积，随后是收缩步骤（shrinkage/soft-threshold）。

对于连续可微函数 $f:{\mathbb{R}}^n\to \mathbb{R}$ 的无约束优化问题 min ⁡ { f ( x ) : x ∈ R n } \min \{ f(x): x\in \mathbb{R} ^n \} min{f(x):x∈Rn} ，解决此类问题的一种简单的方法是使用梯度算法，通过 $x_k=x_{k-1}-t_k\nabla f(x_{k-1})$ 生成序列 { x k } \{ x_k \} {xk​}. 这种梯度迭代法可以被视为一种线性函数 $f$ 在点 $x_{k-1}$ 处的近端正则化（proximal regularization）[Martinet, Bernard ,1970]，即

$$x_k=\arg \underset{x}{\min }\left\{f(x_{k-1})+⟨x-x_{k-1},\nabla f(x_{k-1})⟩+\frac{1}{2t_k}\Vert x-x_{k-1}{\Vert }^2\right\}.$$

其中，$⟨x,y⟩=x^{T}y$ 表示两向量的内积。将这种梯度思想运用到非光滑的含 $l1$-正则项 的优化问题（P）上，那么可以得到迭代策略：

$$x_k=\arg \underset{x}{\min }\left\{f(x_{k-1})+⟨x-x_{k-1},\nabla f(x_{k-1})⟩+\frac{1}{2t_k}\Vert x-x_{k-1}{\Vert }^2+\lambda \Vert x{\Vert }_1\right\}.$$

合并同类项，忽略常数项后，我们能得到如下形式：

$$x_k=\arg \underset{x}{\min }\left\{\frac{1}{2t_k}\Vert x-(x_{k-1}-t_k\nabla f(x_{k-1})){\Vert }^2+\lambda \Vert x{\Vert }_1\right\}.$$

这是一种特殊情况，因为 $l1$-范数 是可以分离的，所以计算 $x_k$ 简化为解决一维的优化问题，即

$$x_k={\tau }_{\lambda t_k}(x_{k-1}-t_k\nabla f(x_{k-1})).$$

针对（P）问题时，ISTA 的迭代形式为：（$P_{ista}$）

$$x_{k+1}={\tau }_{\lambda t_{k+1}}(x_k-t_{k+1}{A}^T(A{x}_k-b)).$$

其中，$t_{k+1}$ 为步长，步长 $t_k\in (0,1/\Vert A^{T}A\Vert )$ 确保了 $x_k$ 可以收敛到 $x^{\ast }$ 。${\tau }_{\alpha }$ 为收缩算子（shrinkage operator），和 shrinkage function 无异。[1] 中证明了下降和收敛性。

这一算法的思想可以追溯到 proximal forwardbackward iterative scheme。最近，一种由 proximal forward-backward algorithms 产生 { x k } \{ x_k \} {xk​} 序列的算法也作出了贡献 [1]。

对于（$P_{ista}$），步长 $t$ 的确定有两种方法，分别是定步长法和回溯法（backtracking）。在使用定步长法时，我们需要设定 $\bar{t}=1/L(f)$ ，$L(f)$ 是 $\nabla f$ 的 Lipschitz 常数。对于（P），$L(f)=2{\lambda }_{\max }(A^{T}A)$ ，这就导致对于大规模问题（$A$ 维数过大），求解 $L(f)$ 往往是困难的，而且，很多时候，$L(f)$ 不可求。所以，我们使用回溯法求解合适的步长。（R）

虽然 ISTA 是一种非常简单的方法，但是它的收敛速度较慢。最近的研究表明，对于一些特定的 $A$ ，ISTA 生成的 { x k } \{ x_k \} {xk​} 序列有着同样慢速的渐进收敛速率，这很糟糕 [Bredies and D. Lorenz, 2008]。

FISTA

那么，现在需要找到一种形式上和 ISTA 一样简单，在速度上却优于 ISTA 的算法。Beck A, Teboulle M 在 2009 年提出 FISTA [1]，FISTA 和 ISTA 主要的不同就在于，收缩算子 $\tau$ 没有用在 $x_{k-1}$ 上，而是用在 $y_k$ 上，$y_k$ 以线性方式结合了前两次的迭代结果 $x_{k-1}$ 和 $x_{k-2}$ ，即

$$x_k={\tau }_{\lambda /L_k}(y_k-\frac{1}{L_k}{A}^T(A{y}_k-b))$$

$$y_{k+1}=x_k+\left(\frac{t_k-1}{t_{k+1}}\right)(x_k-x_{k-1}).$$

此外，$t_k$ 的迭代策略为：

$$t_{k+1}=\frac{1+\sqrt{1+4t_k^2}}{2}.$$

注意，这里 $t_{k+1}$ 并非步长，而是 $y_{k+1}$ 线性结合前两次迭代结果的系数。此时的步长为 $L$ ，$L$ 恰好和之前的步长互为倒数。和上述原因（R）相同，我们也使用回溯法求解合适的步长。

算法描述

ADMM

ADMM Algorithm for LASSO problem based on [4-5]

% Solves the following problem via ADMM:
% minimize 1/2*|| Ax - b ||_2^2 + \lambda || x ||_1

% INPUT
%=======================================
% A
% b
% rho ....... augmented Lagrangian parameter
% lambda .... coefficient of l1-norm
% iter ...... iteration number
% OUTPUT
%=======================================
% x_s ....... sequences {xk} generated by ADMM

1. 初始化 $x_0,z_0,{\mu }_0=0$ [6]
2. 更新变量（$k\ge 0$）：

$$\begin{array}{lcc}{x}^{k+1}& :=(A^{T}A+\rho I{)}^{-1}(A^{T}b+\rho z^k-{\mu }^k)& x.minimization\\ z^{k+1}& :=S_{\lambda /\rho }(x^{k+1}+{\mu }^k/\rho )& z.minimization\\ {\mu }^{k+1}& :={\mu }^k+\rho (x^{k+1}-z^{k+1})& \mathrm{dual}\mathrm{update}\end{array}$$

[F]ISTA

[F]ISTA with backtracking for LASSO problem based on [1]

% Solves the following problem via [F]ISTA:
% minimize 1/2*|| Ax - b ||_2^2 + \lambda || x ||_1

% INPUT
%=======================================
% A
% b
% x0......... initial point
% L0 ........ initial choice of stepsize
% eta ....... the constant in which the stepsize is multiplied
% lambda .... coefficient of l1-norm
% iter ...... iteration number
% opt ....... 0 for ISTA or 1 for FISTA
% eps ....... stop criterion
% OUTPUT
%=======================================
% x_s ....... sequences {xk} generated by [F]ISTA

1. 取 $L_0>0,\eta >1,x_0\in {\mathbb{R}}^n$, 设置 $y_1=x_0,t_1=1$ ;

2. 第 $k\ge 1$ 步，寻找最小非负整数 $i_k$ , 使得对于 $\bar{L}={\eta }^{i_k}{L}_{k-1}$, 有
   $F(p_{\bar{L}}(y_k))\le Q_{\bar{L}}(p_{\bar{L}}(y_k),y_k).$ [1]

   其中，

   $$F(x)\equiv f(x)+g(x),$$

   $$\nabla f(y)=A^T(Ay-b),$$

   $$p_{\bar{L}}(y_k)={\tau }_{\lambda /\bar{L}}(y_k-\frac{1}{\bar{L}}{A}^T(A{y}_k-b)),$$

   $$Q_{\bar{L}}(x,y):=f(y)+⟨x-y,\nabla f(y)⟩+\frac{\bar{L}}{2}\Vert x-y{\Vert }^2+g(x).$$

3. 设置 $L_k={\eta }^{i_k}{L}_{k-1}$ ，更新如下变量：[2]

   $x_k=p_{L_k}(y_k),$

   $$t_{k+1}=\frac{1+\sqrt{1+4t_k^2}}{2},$$

   $${\delta }_k=0\mathrm{for}\mathrm{ISTA},or\frac{t_k-1}{t_{k+1}}\mathrm{for}\mathrm{FISTA},$$

   $$y_{k+1}=x_k+{\delta }_k(x_k-x_{k-1}).$$

注：Remark 3.2 [1] 给出了 $L_k$ 的取值范围：

$$L_0\le L_k\le \eta L(f).$$

实验步骤

设定 $m=1500,n=5000,p=0.02$ ，$p$ 是真值 $u$ 的稀疏度（sparsity density），$A\in {\mathbb{R}}^{m\times n}$. b = A*u + sqrt(0.001)*randn(m,1); 在 $b$ 上添加了随机噪声，初值 $x_0=0$ ，$\lambda =0.15$ . 误差的估计使用 $E(\hat{x})=\Vert \hat{x}-u{\Vert }_1$ ，$\hat{x}$ 是预测值。

clc;clear

randn('seed', 0);
rand('seed',0);

m = 1500;       % number of examples
n = 5000;       % number of features
p = 100/n;      % sparsity density

u = sprandn(n,1,p);
A = randn(m,n);
A = A*spdiags(1./sqrt(sum(A.^2))',0,n,n); % normalize columns
b = A*u + sqrt(0.001)*randn(m,1);

iter = 70;
lambda = 0.15;
x0 = zeros(5000,1);

E = @(x) sum(abs(x-u)); % compute error

ADMM [6]

设置 $\rho =0.7$ .

%% ADMM DEMO
rho = 0.7;
x_admm = admm_lasso(A,b,rho,lambda,iter);
conv_admm = E(x_admm);

function x_s = admm_lasso(A,b,rho,lambda,iter)
% Solves the following problem via ADMM:
% minimize 1/2*|| Ax - b ||_2^2 + \lambda || x ||_1

% INPUT
%=======================================
% A
% b
% rho ....... augmented Lagrangian parameter
% lambda .... coefficient of l1-norm
% iter ...... iteration number
% OUTPUT
%=======================================
% x_s ....... sequences {xk} generated by ADMM

%% shrinkage operator
S = @(tau, g) max(0, g - tau) + min(0, g + tau);

x_s = [];
[~,n] = size(A);
I = eye(n);
x = zeros(n,1);
z_old = zeros(n,1);
u_old = zeros(n,1);
%% MAIN LOOP
for ii = 1:iter
    % record x_s
    x_s = [x_s, x];
    % minimize x,z,u
    x = (A'*A+rho*I) \ (A'*b+rho*z_old-u_old);
    z_new = S(lambda/rho, x+u_old/rho);
    u_new = u_old + rho*(x-z_new);
    % check stop criteria
    % e = norm(x_new-x_old,1)/numel(x_new);
    % if e < eps
        % break
    % end
    % update
    z_old = z_new;
    u_old = u_new;
end

[F]ISTA [7]

设置 $L_0=1.05,\eta =1.01$ .

%% [F]ISTA DEMO
L0 = 1.05;
eta = 1.01;
x_ista = fista_backtracking_lasso(A,b,x0,L0,eta,lambda,iter,0,0);
conv_ista = E(x_ista);
x_fista = fista_backtracking_lasso(A,b,x0,L0,eta,lambda,iter,1,0);
conv_fista = E(x_fista);

function x_s = fista_backtracking_lasso(A,b,x0,L0,eta,lambda,iter,opt, eps)
% Solves the following problem via [F]ISTA:
% minimize 1/2*|| Ax - b ||_2^2 + \lambda || x ||_1

% INPUT
%=======================================
% A
% b
% x0......... initial point
% L0 ........ initial choice of stepsize
% eta ....... the constant in which the stepsize is multiplied
% lambda .... coefficient of l1-norm
% iter ...... iteration number
% opt ....... 0 for ISTA or 1 for FISTA
% eps ....... stop criterion
% OUTPUT
%=======================================
% x_s ....... sequences {xk} generated by [F]ISTA

%% f = 1/2*|| Ax - b ||_2^2
f = @(x) 0.5 * norm(A*x-b)^2;

%% g = \lambda || x ||_1
g = @(x) lambda * norm(x,1);

%% the gradient of f
grad = @(x) A'*(A*x-b);

%% computer F
F = @(x) 0.5*(norm(A*x-b))^2 + lambda*norm(x,1);

%% shrinkage operator
S = @(tau, g) max(0, g - tau) + min(0, g + tau);

%% projection
P = @(L, y) S(lambda/L, y - (1/L)*grad(y));

%% computer Q
Q = @(L, x, y) f(y) + (x-y)'*grad(y) + 0.5*L*norm(x-y)^2 + g(x);

x_s = [];
x_old = x0;
y_old = x0;
L_new = L0;
t_old = 1;
%% MAIN LOOP
for ii = 1:iter
    % find i_k
    j = 1;
    while true
        L_bar = eta^j * L_new;
        if F(P(L_bar, y_old)) <= Q(L_bar, P(L_bar, y_old), y_old)
            L_new = L_new * eta^j;
            break
        else
            j = j + 1;
        end
    end
    x_new = P(L_new, y_old);
    t_new = 0.5 * (1+sqrt(1+4*t_old^2));
    del = opt * (t_old-1)/(t_new);
    y_new = x_new + del*(x_new-x_old);
    % record x_s
    x_s = [x_s, x_new];
    % check stop criteria
    % e = norm(x_new-x_old,1)/numel(x_new);
    % if e < eps
        % break
    % end
    % update
    x_old = x_new;
    t_old = t_new;
    y_old = y_new;
end

注：这里将 ISTA 和 FISTA 算法合并成一个，可以使用 opt 参数进行选择。在函数中，加了停止条件 check stop criteria，为方便测试，这里不使用，故注释掉。

绘图

%% draw
plot(conv_admm,'LineWidth',1,...
    'DisplayName','ADMM','Color','black')
hold on
plot(conv_ista,'--','LineWidth',1,...
    'DisplayName','ISTA','Color','blue')
plot(conv_fista,'-.','LineWidth',1,...
    'DisplayName','FISTA','Color','red')
hold off
ylim([0 270])
xlabel('iteration k')
ylabel('$E(\bar{x})$','Interpreter','latex')
legend('Location','northeast')
title('Convergence behavior in terms of error of iterates of ADMM, ISTA, FISTA.',...
    '$m=1500,n=5000,p=0.02,x_0=0,\lambda=0.15,\rho = 0.7,L_0=1.05,\eta=1.01$',...
    'Interpreter','latex')

结果分析

在此实验中，ADMM 算法的收敛速率最快，FISTA 的收敛速度和 ADMM 相当，ISTA 最慢。在 ADMM 中，设置 $\rho =0.7$ ，$\rho$ 较小 ，虽然获得了不错的收敛速率，但可能在一定程度上削弱了约束条件 $x=z$. 此外，ADMM 在第二次迭代时（所有方法都存在这个问题，暂时不清楚原因），$E(\bar{x})$ 波动幅度较大。此外，FISTA 的收敛速率明显优于 ISTA，在 20 次迭代左右时，基本达到的最优值。

在实验时，ADMM 的求解速度明显慢于 ISTA 和 FISTA，时间主要消耗在计算 $x_{k+1}$ 时的求逆操作上。

2023.05.25更新说明

感谢@ericlhh和@秋意思故人对本文提出的修改建议，即评论区置顶评论，现已将FISTA代码里面计算Q的部分改为了计算二范数|x-y|的平方。具体修改如下：

%% computer Q
- Q = @(L, x, y) f(y) + (x-y)'*grad(y) + 0.5*L*norm(x-y) + g(x);
+ Q = @(L, x, y) f(y) + (x-y)'*grad(y) + 0.5*L*norm(x-y)^2 + g(x);

原创文章！转载需注明来源：©️ Sylvan Ding’s Blog ❤️

参考文献

1. Beck A, Teboulle M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems[J]. SIAM journal on imaging sciences, 2009, 2(1): 183-202.
2. Tao S, Boley D, Zhang S. Convergence of common proximal methods for l1-regularized least squares[C]//Twenty-Fourth International Joint Conference on Artificial Intelligence. 2015.
3. Selesnick I. A derivation of the soft-thresholding function[J]. Polytechnic Institute of New York University, 2009.
4. S. Boyd. Alternating Direction Method of Multipliers. Available online at https://web.stanford.edu/class/ee364b/lectures/admm_slides.pdf.
5. Ryan Tibshirani. Alternating Direction Method of Multipliers. Available online at https://stat.cmu.edu/~ryantibs/convexopt/lectures/admm.pdf.
6. You are allowed to develop your solvers by adopting the ADMM codes available from the following Stanford University websites: https://web.stanford.edu/~boyd/papers/admm/lasso/lasso.html.
7. Tiep Vu. fista_backtracking.m. Available online at https://github.com/tiepvupsu/FISTA.