ADMM,ISTA,FISTA算法步骤详解,MATLAB代码,求解LASSO优化问题
原创文章!转载需注明来源:©️ Sylvan Ding’s Blog ❤️
实验目的
- 了解 ADMM, ISTA, FISTA 算法的基本原理、收敛性和复杂度;
- 使用上述三种算法,解决 LASSO 问题;
- 分析三种算法的表现情况。
实验环境
MATLAB R2021a
文章目录
- ADMM,ISTA,FISTA算法步骤详解,MATLAB代码,求解LASSO优化问题
- 实验目的
- 实验环境
- 实验内容
- LASSO Problem
- ADMM
- ISTA
- FISTA
- 算法描述
- ADMM
- [F]ISTA
- 实验步骤
- ADMM [6]
- [F]ISTA [7]
- 绘图
- 结果分析
- 参考文献
实验内容
ADMM 和 [F]ISTA 的收敛性证明在 [1] 和 [2] 中给出,ADMM 和 ISTA 的全局收敛速率(global convergence rate)为 O(1/k)O(1/k)O(1/k) [Eckstein and Bertsekas, 1990; Deng and Yin, 2012; He and Yuan, 2012],然而 FISTA 的全局收敛速率为 O(1/k2)O(1/k^2)O(1/k2) 。但是,[2] 得出了在某些情况下,ISTA 却比 FISTA 更快的结论,并给出了具体原因。
LASSO Problem
许多机器学习和数据拟合问题都可以看成是最小二乘问题,添加正则项从而避免过拟合。在许多应用中,使用 l1l1l1-范数 作为正则项可以带来不错的泛化效果,所以我们求解含有 l1l1l1-正则项 的线性最小二乘问题,也称为 LASSO 问题 [2],其一般形式为:(P)
minx∈Rn12∥Ax−b∥2+λ∥x∥1.\min \limits _{x\in \mathbb{R} ^n} \frac{1}{2} \Vert Ax-b \Vert ^2 + \lambda \Vert x \Vert _1.x∈Rnmin21∥Ax−b∥2+λ∥x∥1.
其中,A∈Rm×nA\in \mathbb{R} ^{m\times n}A∈Rm×n 是一个行满秩矩阵,n>mn>mn>m ;bbb 是一个已知向量; λ>0\lambda >0λ>0 是一标量;l1l1l1-正则项 ∥x∥1\Vert x \Vert _1∥x∥1 会产生一稀疏解,在降低代价的同时,避免发生过拟合 [Tibshirani, 1996]。
在实验一和实验三中,分别使用了不同的正则项实现了信号的降噪,体现了 l1l1l1-正则项 的另一个优势,就是相较于 l2l2l2-正则项 来说,l1l1l1-正则项 对边界值(outliers)不敏感,这一特点有利于图像降噪中对锐利边缘(sharp edges)的处理。
在实验三中,曾将(P)视为盒状约束(box-constrained)的二次优化问题,使用梯度投影法进行求解。本文要求使用 ADMM, ISTA, FISTA 这三种算法,求解 LASSO 问题。
ADMM
ADMM 算法(Alternating Direction Method of Multipliers)[Gabay, Mercier, Glowinski, Marrocco, 1976] 将原问题的变量 xxx 分解为两个变量 xxx 和 zzz ,在最小二乘的损失函数里使用 xxx ,在 l1l1l1-范数 的正则项里使用 zzz 。这样,外加一个等式约束,就可构造增广拉格朗日函数,进而分别求解两个变量的最优值,使得计算较为容易。
问题(P)等价于:
minx12∥Ax−b∥2+λ∥z∥1s.t.x−z=0\begin{array}{l} \min \limits _x & \frac{1}{2} \Vert Ax-b \Vert ^2 + \lambda \Vert z \Vert _1 \\ \mathrm{s.t.} & x-z=0 \end{array}xmins.t.21∥Ax−b∥2+λ∥z∥1x−z=0
构造增广拉格朗日函数(augmented Lagrangian) LρL_\rhoLρ :
Lρ(x,z,μ)=f(x)+g(z)+μT(x−z)+ρ2∥x−z∥2.L_\rho (x,z,\mu) = f(x)+g(z)+\mu ^T(x-z) + \frac{\rho }{2} \Vert x-z \Vert ^2.Lρ(x,z,μ)=f(x)+g(z)+μT(x−z)+2ρ∥x−z∥2.
其中,f(x)=12∥Ax−b∥2f(x)=\frac{1}{2} \Vert Ax-b \Vert ^2f(x)=21∥Ax−b∥2 ,g(z)=λ∥z∥1g(z)=\lambda \Vert z \Vert _1g(z)=λ∥z∥1 ,ρ\rhoρ 是惩罚参数。
不难发现,最优化 xxx 时,是一个二次优化问题,而最优化 zzz 时,是一个软阈(soft-thresholding)问题,那么,ADMM 的迭代形式为:
xk+1:=(ATA+ρI)−1(ATb+ρzk−μk)x.minimizationzk+1:=Sλ/ρ(xk+1+μk/ρ)z.minimizationμk+1:=μk+ρ(xk+1−zk+1)dualupdate\begin{array}{l} x^{k+1} &:= (A^TA+\rho I)^{-1}(A^Tb+\rho z^k-\mu ^k) & \mathrm{ x.minimization } \\ z^{k+1} &:= S_{\lambda / \rho} (x^{k+1} + \mu ^{k}/ \rho) & \mathrm{ z.minimization }\\ \mu ^{k+1} &:= \mu ^k + \rho (x^{k+1}-z^{k+1}) & \mathrm{ dual \quad update } \end{array}xk+1zk+1μk+1:=(ATA+ρI)−1(ATb+ρzk−μk):=Sλ/ρ(xk+1+μk/ρ):=μk+ρ(xk+1−zk+1)x.minimizationz.minimizationdualupdate
其中,SτS_\tauSτ 是软阈函数(又称为Shrinkage),Sτ(g):=sign(g)⋅(∣g∣−τ)+S_\tau (g):=sign(g)\cdot (|g|-\tau)_+Sτ(g):=sign(g)⋅(∣g∣−τ)+ [3],也即 [5]
[Sτ(g)]j={gj−τg>τ0−τ≤g≤τgj+τg<−τ,j=1,2,…,n\left [S_\tau (g) \right ]_j= \left\{\begin{matrix} g_j-\tau & g> \tau \\ 0 & -\tau \le g \le \tau \\ g_j+\tau & g<-\tau \end{matrix}\right. , \quad j=1,2,\dots ,n[Sτ(g)]j=⎩⎨⎧gj−τ0gj+τg>τ−τ≤g≤τg<−τ,j=1,2,…,n
实际上,ADMM 算法往往在一系列迭代后,能获得较为精确的解,但是所需要的迭代次数是大量的 [5]。
ρ\rhoρ 的选取会在极大程度上影响 ADMM 的收敛性。若 ρ\rhoρ 太大,则对于优化 f+gf+gf+g 的能力会下降;反之,则会削弱约束条件 x=zx=zx=z. Boyd et al. (2010) 给出了选 ρ\rhoρ 的策略,效果不错,但不能保证一定收敛。
ISTA
现在流行的一种解决问题(P)的方法是 ISTA(iterative shrinkage-thresholding algorithms)[Daubechies et al, 2004],该方法在每次迭代时会计算矩阵和向量的乘积,随后是收缩步骤(shrinkage/soft-threshold)。
对于连续可微函数 f:Rn→Rf:\mathbb{R} ^n \to \mathbb{R}f:Rn→R 的无约束优化问题 min{f(x):x∈Rn}\min \{ f(x): x\in \mathbb{R} ^n \}min{f(x):x∈Rn} ,解决此类问题的一种简单的方法是使用梯度算法,通过 xk=xk−1−tk∇f(xk−1)x_k=x_{k-1}-t_k\nabla f(x_{k-1})xk=xk−1−tk∇f(xk−1) 生成序列 {xk}\{ x_k \}{xk}. 这种梯度迭代法可以被视为一种线性函数 fff 在点 xk−1x_{k-1}xk−1 处的近端正则化(proximal regularization)[Martinet, Bernard ,1970],即
xk=argminx{f(xk−1)+⟨x−xk−1,∇f(xk−1)⟩+12tk∥x−xk−1∥2}.x_k=\arg \min \limits _x \left \{ f(x_{k-1} ) + \left \langle x-x_{k-1},\nabla f(x_{k-1}) \right \rangle + \frac{1}{2t_k} \Vert x-x_{k-1} \Vert ^2 \right \}.xk=argxmin{f(xk−1)+⟨x−xk−1,∇f(xk−1)⟩+2tk1∥x−xk−1∥2}.
其中,⟨x,y⟩=xTy\left \langle x,y \right \rangle =x^Ty⟨x,y⟩=xTy 表示两向量的内积。将这种梯度思想运用到非光滑的含 l1l1l1-正则项 的优化问题(P)上,那么可以得到迭代策略:
xk=argminx{f(xk−1)+⟨x−xk−1,∇f(xk−1)⟩+12tk∥x−xk−1∥2+λ∥x∥1}.x_k=\arg \min \limits _x \left \{ f(x_{k-1} ) + \left \langle x-x_{k-1},\nabla f(x_{k-1}) \right \rangle + \frac{1}{2t_k} \Vert x-x_{k-1} \Vert ^2 + \lambda \Vert x \Vert _1 \right \}.xk=argxmin{f(xk−1)+⟨x−xk−1,∇f(xk−1)⟩+2tk1∥x−xk−1∥2+λ∥x∥1}.
合并同类项,忽略常数项后,我们能得到如下形式:
xk=argminx{12tk∥x−(xk−1−tk∇f(xk−1))∥2+λ∥x∥1}.x_k = \arg \min \limits _x \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 \}.xk=argxmin{2tk1∥x−(xk−1−tk∇f(xk−1))∥2+λ∥x∥1}.
这是一种特殊情况,因为 l1l1l1-范数 是可以分离的,所以计算 xkx_kxk 简化为解决一维的优化问题,即
xk=τλtk(xk−1−tk∇f(xk−1)).x_{k}=\tau _{\lambda t_k} (x_{k-1} -t_k\nabla f(x_{k-1})).xk=τλtk(xk−1−tk∇f(xk−1)).
针对(P)问题时,ISTA 的迭代形式为:(PistaP_{ista}Pista)
xk+1=τλtk+1(xk−tk+1AT(Axk−b)).x_{k+1}=\tau _{\lambda t_{k+1}} (x_k -t_{k+1}A^T(Ax_k-b)).xk+1=τλtk+1(xk−tk+1AT(Axk−b)).
其中,tk+1t_{k+1}tk+1 为步长,步长 tk∈(0,1/∥ATA∥)t_k\in (0,1/\Vert A^TA \Vert )tk∈(0,1/∥ATA∥) 确保了 xkx_kxk 可以收敛到 x∗x^*x∗ 。τα\tau _{\alpha}τα 为收缩算子(shrinkage operator),和 shrinkage function 无异。[1] 中证明了下降和收敛性。
这一算法的思想可以追溯到 proximal forwardbackward iterative scheme。最近,一种由 proximal forward-backward algorithms 产生 {xk}\{ x_k \}{xk} 序列的算法也作出了贡献 [1]。
对于(PistaP_{ista}Pista),步长 ttt 的确定有两种方法,分别是定步长法和回溯法(backtracking)。在使用定步长法时,我们需要设定 tˉ=1/L(f)\bar{t}=1/L(f)tˉ=1/L(f) ,L(f)L(f)L(f) 是 ∇f\nabla f∇f 的 Lipschitz 常数。对于(P),L(f)=2λmax(ATA)L(f)=2\lambda _{\max} (A^TA)L(f)=2λmax(ATA) ,这就导致对于大规模问题(AAA 维数过大),求解 L(f)L(f)L(f) 往往是困难的,而且,很多时候,L(f)L(f)L(f) 不可求。所以,我们使用回溯法求解合适的步长。(R)
虽然 ISTA 是一种非常简单的方法,但是它的收敛速度较慢。最近的研究表明,对于一些特定的 AAA ,ISTA 生成的 {xk}\{ x_k \}{xk} 序列有着同样慢速的渐进收敛速率,这很糟糕 [Bredies and D. Lorenz, 2008]。
FISTA
那么,现在需要找到一种形式上和 ISTA 一样简单,在速度上却优于 ISTA 的算法。Beck A, Teboulle M 在 2009 年提出 FISTA [1],FISTA 和 ISTA 主要的不同就在于,收缩算子 τ\tauτ 没有用在 xk−1x_{k-1}xk−1 上,而是用在 yky_kyk 上,yky_kyk 以线性方式结合了前两次的迭代结果 xk−1x_{k-1}xk−1 和 xk−2x_{k-2}xk−2 ,即
xk=τλ/Lk(yk−1LkAT(Ayk−b))x_{k}=\tau _{\lambda /L_{k}} (y_k -\frac{1}{L_k}A^T(Ay_k-b))xk=τλ/Lk(yk−Lk1AT(Ayk−b))
yk+1=xk+(tk−1tk+1)(xk−xk−1).y_{k+1}=x_k+\left ( \frac{t_k-1}{t_{k+1}} \right )(x_k - x_{k-1}).yk+1=xk+(tk+1tk−1)(xk−xk−1).
此外,tkt_ktk 的迭代策略为:
tk+1=1+1+4tk22.t_{k+1}=\frac{1+\sqrt{1+4t^2_k}}{2}.tk+1=21+1+4tk2.
注意,这里 tk+1t_{k+1}tk+1 并非步长,而是 yk+1y_{k+1}yk+1 线性结合前两次迭代结果的系数。此时的步长为 LLL ,LLL 恰好和之前的步长互为倒数。和上述原因(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
- 初始化 x0,z0,μ0=0x_0,z_0,\mu _0 = \mathrm{0}x0,z0,μ0=0 [6]
- 更新变量(k≥0k \ge 0k≥0):
xk+1:=(ATA+ρI)−1(ATb+ρzk−μk)x.minimizationzk+1:=Sλ/ρ(xk+1+μk/ρ)z.minimizationμk+1:=μk+ρ(xk+1−zk+1)dualupdate\begin{array}{l} x^{k+1} &:= (A^TA+\rho I)^{-1}(A^Tb+\rho z^k-\mu ^k) & \mathrm{ x.minimization } \\ z^{k+1} &:= S_{\lambda / \rho} (x^{k+1} + \mu ^{k}/ \rho) & \mathrm{ z.minimization }\\ \mu ^{k+1} &:= \mu ^k + \rho (x^{k+1}-z^{k+1}) & \mathrm{ dual \quad update } \end{array}xk+1zk+1μk+1:=(ATA+ρI)−1(ATb+ρzk−μk):=Sλ/ρ(xk+1+μk/ρ):=μk+ρ(xk+1−zk+1)x.minimizationz.minimizationdualupdate
[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
-
取 L0>0,η>1,x0∈RnL_0>0 ,\eta>1,x_0 \in \mathbb{R} ^nL0>0,η>1,x0∈Rn, 设置 y1=x0,t1=1y_1=x_0,t_1=1y1=x0,t1=1 ;
-
第 k≥1k \ge 1k≥1 步,寻找最小非负整数 iki_kik , 使得对于 Lˉ=ηikLk−1\bar{L}=\eta ^{i_k}L_{k-1}Lˉ=ηikLk−1, 有
F(pLˉ(yk))≤QLˉ(pLˉ(yk),yk).F(p_{\bar{L}}(y_k))\le Q_{\bar{L}} (p_{\bar{L}}(y_k), y_k).F(pLˉ(yk))≤QLˉ(pLˉ(yk),yk). [1]其中,
F(x)≡f(x)+g(x),F(x) \equiv f(x)+g(x),F(x)≡f(x)+g(x),
∇f(y)=AT(Ay−b),\nabla f(y) = A^T(Ay-b) ,∇f(y)=AT(Ay−b),
pLˉ(yk)=τλ/Lˉ(yk−1LˉAT(Ayk−b)),p_{\bar{L}} (y_k) = \tau _{\lambda /\bar{L}} (y_k -\frac{1}{\bar{L}}A^T(Ay_k-b)) ,pLˉ(yk)=τλ/Lˉ(yk−Lˉ1AT(Ayk−b)),
QLˉ(x,y):=f(y)+⟨x−y,∇f(y)⟩+Lˉ2∥x−y∥2+g(x).Q_{\bar{L}}(x,y):=f(y)+\left \langle x-y,\nabla f(y) \right \rangle + \frac{\bar{L}}{2} \Vert x-y \Vert ^2 + g(x) .QLˉ(x,y):=f(y)+⟨x−y,∇f(y)⟩+2Lˉ∥x−y∥2+g(x).
-
设置 Lk=ηikLk−1L_k=\eta ^{i_k} L_{k-1}Lk=ηikLk−1 ,更新如下变量:[2]
xk=pLk(yk),x_{k}= p_{L_k} (y_k),xk=pLk(yk),
tk+1=1+1+4tk22,t_{k+1}=\frac{1+\sqrt{1+4t^2_k}}{2},tk+1=21+1+4tk2,
δk=0forISTA,ortk−1tk+1forFISTA,\delta _k=0 \ \mathrm{for \ ISTA},\ or\ \frac{t_k-1}{t_{k+1}}\ \mathrm{for \ FISTA},δk=0 for ISTA, or tk+1tk−1 for FISTA,
yk+1=xk+δk(xk−xk−1).y_{k+1} = x_{k} + \delta _k (x_k-x_{k-1}) .yk+1=xk+δk(xk−xk−1).
注:Remark 3.2 [1] 给出了 LkL_kLk 的取值范围:
L0≤Lk≤ηL(f).L_0\le L_k \le \eta L(f) .L0≤Lk≤ηL(f).
实验步骤
设定 m=1500,n=5000,p=0.02m=1500,n=5000,p=0.02m=1500,n=5000,p=0.02 ,ppp 是真值 uuu 的稀疏度(sparsity density),A∈Rm×nA\in \mathbb{R} ^{m\times n}A∈Rm×n. b = A*u + sqrt(0.001)*randn(m,1);
在 bbb 上添加了随机噪声,初值 x0=0x_0=0x0=0 ,λ=0.15\lambda=0.15λ=0.15 . 误差的估计使用 E(x^)=∥x^−u∥1E(\hat{x}) = \Vert \hat{x}-u \Vert _1E(x^)=∥x^−u∥1 ,x^\hat{x}x^ 是预测值。
clc;clearrandn('seed', 0);
rand('seed',0);m = 1500; % number of examples
n = 5000; % number of features
p = 100/n; % sparsity densityu = 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]
设置 ρ=0.7\rho = 0.7ρ=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_sx_s = [x_s, x];% minimize x,z,ux = (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% updatez_old = z_new;u_old = u_new;
end
[F]ISTA [7]
设置 L0=1.05,η=1.01L_0=1.05,\eta=1.01L0=1.05,η=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) + g(x);x_s = [];
x_old = x0;
y_old = x0;
L_new = L0;
t_old = 1;
%% MAIN LOOP
for ii = 1:iter% find i_kj = 1;while trueL_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;breakelsej = j + 1;endendx_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_sx_s = [x_s, x_new];% check stop criteria% e = norm(x_new-x_old,1)/numel(x_new);% if e < eps% break% end% updatex_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 中,设置 ρ=0.7\rho = 0.7ρ=0.7 ,ρ\rhoρ 较小 ,虽然获得了不错的收敛速率,但可能在一定程度上削弱了约束条件 x=zx=zx=z. 此外,ADMM 在第二次迭代时(所有方法都存在这个问题,暂时不清楚原因),E(xˉ)E(\bar{x})E(xˉ) 波动幅度较大。此外,FISTA 的收敛速率明显优于 ISTA,在 20 次迭代左右时,基本达到的最优值。
在实验时,ADMM 的求解速度明显慢于 ISTA 和 FISTA,时间主要消耗在计算 xk+1x_{k+1}xk+1 时的求逆操作上。
原创文章!转载需注明来源:©️ Sylvan Ding’s Blog ❤️
参考文献
- 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.
- 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.
- Selesnick I. A derivation of the soft-thresholding function[J]. Polytechnic Institute of New York University, 2009.
- S. Boyd. Alternating Direction Method of Multipliers. Available online at https://web.stanford.edu/class/ee364b/lectures/admm_slides.pdf.
- Ryan Tibshirani. Alternating Direction Method of Multipliers. Available online at https://stat.cmu.edu/~ryantibs/convexopt/lectures/admm.pdf.
- 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.
- Tiep Vu. fista_backtracking.m. Available online at https://github.com/tiepvupsu/FISTA.