二阶变分图像恢复模型的重启动快速ADMM方法

宋田田，潘振宽，魏伟波，李青

青岛大学计算机科学技术学院,青岛 266071

0 引 言

经过近30年的发展，变分方法已经成为图像处理与分析的重要数学基础(Scherzer，2015；Paragios等，2006；Aubert和Kornprobst，2006；Chan和Shen，2005)，其中，图像恢复的变分模型具有基础地位。该类模型一般包含数据项和规则项，出于对图像特征保持的考虑，其规则项往往是非线性、非光滑，甚至是非凸的。这些特点不仅导致数值算法设计的困难，且其数值方法的计算效率往往较低，成为变分图像处理与分析模型工程应用的主要瓶颈之一(Kimmel和Tai，2019；Glowinski等，2016)。

以总变差(total variation，TV)(Rudin等，1992)为代表的一阶变分模型可以保持图像的边缘特征，是变分图像处理与分析问题的基础模型。为了克服基于原变量的时间步进方法(Rudin等，1992)和固定点迭代方法(Vogel和Oman，1996)计算效率低的问题，Chan等人(1999)提出了原—对偶变量方法(primal-dual method，PD)，Chambolle(2004)提出了经典的对偶方法(dual method，DM)，Goldstein和Osher(2009)提出了SB(split bregman)方法。Wu和Tai(2010)证明，SB算法与增广Lagrange方法(augmented Lagrangian method，ALM)或交替方向乘子法(alternating direction methods of multipliers，ADMM)(Glowinski和Le Tallec，1989)是等价的。Daubechies等人(2004)以交替迭代优化为框架提出了经典的软阈值迭代算法(iterative shirinkage-thresholding algorithm，ISTA)。上述方法均为基于一阶导数或梯度的一阶方法，其收敛率为O (1/k)，k为迭代次数。

早在20世纪80年代初，Nesterov(1983)针对光滑凸目标函数的优化问题，基于函数外插或惯性加速的思想，通过设计适当的惯性参数，将基于梯度的方法的算法收敛率提高至O (1/k2)，为变分图像处理模型的加速算法设计奠定了基础。Beck和Teboulle(2009a，b)针对组合目标函数结合Nesterov方法与ISTA方法提出了快速ISTA方法，即快速软阈值迭代算法(fast iterative shirinkage-thresholding algorithm，FISTA)。基于相似的思路，Beck和Teboulle(2014)提出了快速对偶方法(fast dual，FD)、Chambolle和Pock(2011)提出了快速原—对偶变量方法(fast PD，FPD)、Zhu和Chan(2008)提出了原—对偶变量混合梯度方法(primal-dual hybrid gradient，PDHG)，Goldstein等人(2014)提出了快速ADMM方法。该类加速方法设计的关键是最优惯性参数的设计，该参数依赖于目标函数的光滑度与凸性，但通常的变分图像处理模型往往局部强凸或完全非凸，导致这些参数估计困难或过分耗时，其惯性加速类算法会引起涟漪或振荡现象(O’Donoghue和Candès，2015)，达不到预期加速效果。改进的单调算法(Beck和Teboulle，2009b)、回溯算法(Goldstein等，2014)或重启动算法(O’Donoghue和Candès，2015)的研究成为近年惯性加速方法实用化的重要方向之一。单调类算法依据相邻两次迭代结果构造惯性加速，回溯方法是通过向后寻找能够使得能量值下降的最小步长来保证算法能量值下降，重启动算法则在能量非下降时重置初始参数，这些方案均能避免振荡现象，并保持算法收敛率仍为O(1/k2)，该研究正在成为凸优化计算领域的热点问题之一(Kim和Fessler，2018；Calatroni和Chambolle，2019；Buccini等，2020)。但相关研究和应用在计算机视觉领域还限于基于一阶导数的变分图像处理模型及经典的Lasso模型(陈少利，2017；李启朋，2018；李星 等，2019；潘树林 等，2019)。

基于二阶导数的变分模型可有效克服一阶导数模型引起的阶梯效应问题，其规则项的形式多样，且由于模型的非线性、非光滑性和非凸性，其计算效率低的问题更加突出。Bredies等人(2010)直接将Nesterov型加速方案应用于基于二阶导数的总广义变差(total generalized variation，TGV)模型计算；Liu等人(2012)设计了基于总海森(total Hessian，TH)变分模型的快速PDHG算法；Yashtini和Kang(2016)则用相同的方法设计了基于欧拉弹性能(Euler’s elastica，EE)的算法设计。Bredies等人(2010)的方法包含2个凸的非光滑规则项，Liu等人(2012)的方法包含1个凸的非光滑规则项，Yashtini和Kang(2016)的方法包含非凸、非光滑规则项，除快速算法本身可能引起的振荡问题，非凸性亦可引起能量上升，但前期研究均未考虑这些因素。重启动策略有望一并克服这些问题。

含二阶导数变分图像恢复模型的规则项主要包括总广义变差(TGV)(Bredies等，2010)、总拉普拉斯(total Laplacian，TL)(Chan等，2010)、总海森(total Hessian，TH)(Lysaker等，2003)、总平均曲率(total mean curvature，TMC)(Zhu等，2014)和欧拉弹性能(Euler’s elastica，EE)(Yashtini和Kang，2016；Tai等，2011；Kang等，2019)等形式，均为非强凸项，其中，TGV、TL和TH为凸非光滑项，TMC和EE为非凸非光滑项。

本文的目的是以快速ADMM方法为框架，探讨重启动快速算法在这些模型算法设计中应用的可能性。原始ADMM方法的收敛速率为O (1/k)，惯性加速方法如快速ADMM理论上可以将收敛速率提高至O (1/k2)。但该类加速方法过分依赖最优惯性参数的设计，一般的变分模型很难估计最优惯性参数，无法找到最优步长，并会产生能量振荡现象，降低收敛速率，无法达到预期效果。重启动算法可以避免振荡现象，从而提高计算效率。

本文针对二阶模型在加速后进行重启动，真正提高算法的收敛速率。重启动快速ADMM方法的时间复杂度的精确分析很困难，故本文选用两种典型的基于拉普拉斯以及基于曲率的二阶模型进行研究，以探讨重启动快速ADMM方法的计算效率问题。本文方法亦可自然拓展到基于类似模型的重启动快速算法。

1 TV模型的重启动快速ADMM算法

图像噪声去除是图像恢复的基本问题，图像噪声去除的TV模型(Rudin等，1992)可表达为能量泛函极值问题，具体为

(1)

(2)

式中，增广Lagrange函数为

在交替优化过程中的两个子优化问题分别为

(3)

(4)

Goldstein等人(2014)证明，采用该加速方案，其算法的收敛率可提高为O(1/k2)。为了避免涟漪现象，Goldstein等人(2014)设计了重启动方法，并提出判断重启动条件的简单形式的组合残差公式作为简化能量变化的依据。具体为

(5)

αk+1=1,wk+1←wk,βk+1←βk,ck+1←ck

(6)

算法1：TV模型的重启动快速ADMM算法

初始化：u0,w0,γ,β0,μ>0,

α0=1,c0>0,η∈(0,1)

fork=0，1，2，3，…，do

uk+1=argminE(u,wk)；

wk+1=argminE(uk+1,w)；

ifck+1<ηckthen

wk+1←wk+1+θk+1(wk+1-wk)；

βk+1←βk+1+θk+1(βk+1-βk)；

else;

αk+1=1；

wk+1←wk；

βk+1←βk；

end if;

直到E(u)收敛；

end for。

重启动快速ADMM方法的每次迭代属于以下3种类型之一：1)当算法中ck+1<ηck不等式未被满足时，发生“重启动”迭代。2)“未加速”迭代在“重启动”迭代之后立即发生。在这样的迭代中，αk=1，因此算法的加速步骤被禁用，使得迭代等同于原始ADMM。3)“加速”迭代是指没有任何“重启”或“未加速”的迭代，在这样的迭代中αk>1，算法中ck+1<ηck不等式被满足，调用算法的加速步骤。

这样，算法使用“重启动”方法增强函数值的单调性。这种类似的重启动规则已经用于针对无约束最小化的研究(Kim和Fessler，2018)。

2 TL模型的重启动快速ADMM算法

基于二阶导数图像恢复变分TL模型是为克服一阶TV模型的阶梯效应问题提出的(Lysaker等，2003)，其能量泛函极小值问题为

(7)

为了克服TL规则项变分引起的计算困难，引入了两组辅助变量、Lagrange乘子及惩罚参数设计其ADMM算法，以便将原复杂优化问题转化为交替优化的一系列简单子优化问题求解。然后，这些子问题可以通过快速傅里叶变换(fast Fourier transform，FFT)和解析形式的软阈值公式进行求解。

(8)

式中，增广Lagrange函数为

式(8)中的3个子优化问题分别为

(9)

使用标准变分方法求解式(9)第1行，可得关于u的Euler-Lagrange方程，具体为

(10)

(11)

同样，使用标准变分方法求解式(9)第2行，可得关于w的Euler-Lagrange方程，即

(12)

该方程组亦可采用FFT求解，即

(13)

对于式(9)第3行，其解v可用解析形式的软阈值公式表示，具体为

(14)

采用与Goldstein等人(2014)和Buccini等人(2020)相同的思路，可设计TL模型原变量与对偶变量组合残差公式，即

(15)

由此，得到应用于TL模型的重启动快速ADMM算法，描述如下。

算法2：TL模型的重启动快速ADMM算法

α0=1,c0>0,η∈(0,1);

fork=0，1，2，3，…，do

uk+1=argminE(u,wk,vk);

wk+1=argminE(uk+1,w,vk);

vk+1=argminE(uk+1,wk+1,vk);

ifck+1<ηckthen;

wk+1←wk+1+θk+1(wk+1-wk);

else

αk+1=1;

wk+1←wk;

end if；

直至E(u)收敛；

end for。

3 EE模型的重启动快速ADMM算法

EE模型结合了TV规则项和曲率项，使用欧拉的弹性项作为规则项，相应的图像噪声去除变分模型为

(16)

欧拉弹性能最早由Nitzberg等人(1993)应用于深度图像分割与虚幻轮廓恢复。Masnou和Morel(1998)将其推广到大破损图像修复，Zhu等人(2013)将其推广到图像分割变分模型。EE模型可以有效减少阶梯效应并产生更接近原图的图像，而且可以在去除噪声的同时保留图像的边缘信息。

(17)

|p|≤1，|w|-p·w≥0

(18)

为了进一步简化每个子优化问题的目标函数，再引入一个新的辅助向量m代替式(18)中的变量p。引入辅助向量m后，模型共有5个约束，即

(19)

在EE模型中，引入4个辅助变量w，p，v和m，4个Lagrange乘子β1，β2，β3和β4以及4个惩罚参数μ1，μ2，μ3和μ4，采用ADMM将式(17)转化为

(20)

式中，增广Lagrange函数为

应用交替优化技术，式(20)中的5个子优化问题分别为

(21)

使用标准变分方法求解式(21)第1行，得到关于u的Euler-Lagrange方程，即

(22)

(23)

对于式(21)第2行，其解w可用解析形式的广义软阈值公式表示，即

(24)

同样，使用标准变分方法求解式(21)第3行，可得关于p的Euler-Lagrange方程，即

(25)

该方程组亦可采用FFT求解，即

(26)

使用标准变分方法求解式(21)第4行，可得关于v的Euler-Lagrange方程，即

(27)

相应的解析解为

(28)

使用标准变分方法求解式(21)第5行，可得关于m的Euler-Lagrange方程，即

(29)

相应的解析解为

(30)

(31)

采用与Goldstein等人(2014)和Buccini等人(2020)相同的策略，设计EE模型的原变量和对偶变量组合残差公式，即

(32)

由此，得到应用于EE模型的重启动快速ADMM算法，描述如下。

算法3：EE模型的重启动快速ADMM算法

(μ1,μ2,μ3,μ4)>0,α0=1,c0>0,η∈(0,1)

fork=0，1，2，3，…，do

uk+1=argminE(u,wk,pk,vk,mk);

wk+1=argminE(uk+1,w,pk,vk,mk);

pk+1=argminE(uk+1,wk+1,p,vk,mk);

vk+1=argminE(uk+1,wk+1,pk+1,v,mk);

mk+1=argminE(uk+1,wk+1,pk+1,vk+1,m);

ifck+1<ηckthen

wk+1←wk+1+θk+1(wk+1-wk);

else

αk+1=1;

wk+1←wk;

end if；

直到E(u)收敛；

end for。

4 数值实验

针对TL模型和 EE模型在噪声去除和边缘保持等方面的性能已经有了大量研究。本文研究的重点是在相关惩罚参数确定的情况下，在保持原有模型噪声去除性态的基础上，比较原始ADMM、快速ADMM 和重启动快速ADMM方法的计算效率。

由于TL模型和EE模型的解决方案通常不是唯一的，无法比较均方根差(root mean square error，RMSE)。实验时，使用相对能量误差|Ek-Ek-1|/Ek<0.001作为停止准则，Ek表示当前第k步迭代的原始能量，Ek-1表示前一步(第k-1步)迭代的原始能量。当相对能量误差无法达到停止标准时，算法将在迭代500次后自动停止。

实验使用3.19 GHz CPU，在Windows 10(64位)操作系统下用MATLAB R2016a版本进行。

4.1 Total Laplacian 模型

TL模型如式(7)所示。f表示含噪声的图像时，最小化式(7)的目的是找到与含噪声的图像f相似的去噪后的图像u，同时保证图像光滑、边缘保持，且不会出现边缘的阶梯效应。参数γ控制数据项和规则项之间的权衡。模型参数γ主要影响图像的质量，为了更清晰地展示惩罚参数μ1和μ2与算法性能之间的关系，固定γ=1。

快速算法的目标是在保证图像质量不降低的基础上，提高算法的收敛速率从而加快计算效率。以Texture图像(图1第3行)为例，图2展示了原始ADMM、快速ADMM 和重启动快速ADMM这3种算法应用到TL模型后输出的图像，3种算法的参数取值均为μ1=5,μ2=0.001,γ=1，峰值信噪比(peak signal-to-noise ratio，PSNR)分别为26.518、26.542和26.969。

图1 测试图像Fig.1 Test images ((a)original images;(b)Gaussian noisy images;(c)salt and pepper noisy images)

图2 去噪效果图Fig.2 Images denoising effect ((a)ADMM;(b)fast ADMM;(c)restart fast ADMM)

惩罚参数的数值理论上可以无穷大，但由于引入了Lagrange方程，惩罚可以通过Lagrange乘子β调节，所以惩罚参数不需要过大。惩罚参数的选择在理论上没有答案，目前都是通过实验调节选取。本实验选取了大量惩罚参数，在不同参数组合的情况下分析3种算法的计算效率。不同的参数对于快速ADMM方法的计算效率影响较大，当μ1=5时，快速ADMM方法均加速无效，对于重启动方法则基本无影响。

表1 TL模型消除高斯白噪声的运行结果Table 1 Running results of TL model to eliminate Gaussian noise

表2 TL模型消除椒盐噪声的运行结果Table 2 Running results of TL model to eliminate salt and pepper noise

为了可以更清晰地了解在计算过程中3种算法的能量值和相对能量误差的变化趋势，以TL模型消除高斯白噪声为例，对3种算法在TL模型中迭代40次的能量值进行比较，结果如图3所示。图4表示使用相对能量误差作为停止准则时，3种算法在TL模型中的收敛曲线。参数取值均为μ1=5、μ2=0.001、γ=1。

从图3和图4可以看出，重启动快速ADMM算法的能量值和相对能量误差是单调下降的，且能量值达到稳定状态最快；快速ADMM算法会产生振荡。这是由于TL模型的非光滑性，使光滑参数和Lipschitz参数很难估计，导致振荡的周期很难估计。重启动快速ADMM算法在计算过程中根据组合残差的大小，自适应地调整步长从而消除振荡现象，提高了计算效率。

图3 3种算法在TL模型中迭代40次的能量值比较Fig.3 Comparison of energy values of the three algorithms iterations 40 times in TL model((a)Lena;(b)Castle;(c)Texture;(d)Barbara;(e)Peppers)

图4 3种算法在TL模型中的收敛曲线Fig.4 Convergence curves of the three algorithms in TL model((a)Lena;(b)Castle;(c)Texture;(d)Barbara;(e)Peppers)

4.2 Euler’s elastica 模型

表3 EE模型消除高斯白噪声的运行结果Table 3 Running results of EE model to eliminate Gaussian noise

表4 EE模型消除椒盐噪声的运行结果Table 4 Running results of EE model to eliminate salt and pepper noise

同样，为了可以更清楚地了解在计算过程中3种算法能量值和相对能量误差的变化趋势，以EE模型消除高斯白噪声为例，对3种算法在EE模型中迭代40次的能量值进行比较，结果如图5所示。图6表示使用相对能量误差作为停止准则时，3种算法在EE模型中的收敛曲线。参数取值均为a=0.2、b=2、μ1=0.06、μ2=2、μ3=2 000、μ4=400。

从图5和图6可以看出:1)通过函数的能量值曲线，重启动快速ADMM的函数值相较于其他两种算法能够更快地接近稳定。2)快速 ADMM算法在图5的5幅图像中都出现了函数值不单调下降的情况，这种缺陷极大降低了算法效率。这是因为快速ADMM无法计算最优步长，加速步长过大，导致错过了极小值从而出现了函数值不单调下降情况，因此产生振荡。3)重启动快速ADMM算法的能量值和相对能量误差是单调下降的，且能量值达到稳定状态最快。

图5 3种算法在EE模型中迭代40次的能量值比较Fig.5 Comparison of energy values of the three algorithms iterations 40 times in EE model((a)Lena;(b)Castle;(c)Texture;(d)Barbara;(e)Peppers)

图6 3种算法在EE模型中的收敛曲线Fig.6 Convergence curves of the three algorithms in EE model((a)Lena;(b)Castle;(c)Texture;(d)Barbara;(e)Peppers)

5 结 论

本文以图像噪声去除为背景，针对含二阶导数规则项的非线性、非光滑的TL模型及非线性、非光滑、非凸的EE模型，以ADMM方法为基础，引入Neserov惯性加速及重启动策略设计了相应的快速算法及重启动快速算法。对不同的变量进行加速，重启动快速算法的计算效率较原始的ADMM方法和快速ADMM方法均有较大提高，为更明显地显示实验效果，本文选取加速效率较为明显的变量进行加速实验。重启动算法的计算效率对所选惩罚参数具有鲁棒性，其计算效率的提高还得益于优化子问题的快速FFT求解及解析形式的软阈值公式的使用。这些有益的探索可为其他形式高阶变分图像处理与分析模型的快速求解提供借鉴。

但是，含高阶导数的非线性、非光滑和非凸变分模型的ADMM方法及其重启动快速算法设计还缺乏足够的理论支撑。目前的理论研究还局限于由两个函数组合的目标函数的优化问题，且确定性的理论成果集中于光滑、强凸及含一个线性约束的优化问题。对于计算机视觉中的含高阶导数的非光滑、非凸变分模型的快速算法研究，本文的工作仅限于尝试性算法设计及数值验证，后续复杂的算法理论分析对算法的推广将具有重要价值。