考虑如下鞍点优化问题
其中$ X\subseteq R^n, Y\subseteq R^m $是两个闭凸集合, $ Q(x, y) $是光滑的凸–凹函数.该模型在图像处理, 高维统计, 机器学习中有着广泛应用.注意到鞍点问题实际上可以转化为单调变分不等式问题, 那么求解单调变分不等式的数值方法[1-3]都可以用来求解问题(1.1).令$ v = (x, y)^T $, $ F(v) = (\nabla_x Q(x, y), -\nabla_y Q(x, y))^T $, $ K = X\times Y $, 那么寻找极小–极大问题(1.1)的鞍点$ (x^*, y^*) $就可以转换为如下形式的变分不等式问题
易知, $ v^* $是上述变分不等式的解当且仅当其满足不定点方程$ v^* = P_K(v^*-\alpha\nabla F(v^*)), $其中$ P $是投影算子, $ \alpha > 0 $是步长.鉴于约束集合$ K $的可分结构, 不动点方程可以具体写为
基于此, Bonettini和Ruggiero在文献[4]中提出了如下外梯度方案求解问题(1.1)
其中$ \alpha_k $是自适应步长, 由线搜索确定.该方法中, 原始变量和对偶变量共用一个步长, 且能自适应确定, 这是该方法的优点, 不必费大量时间调参.收敛性也能在较弱的条件下证明.该方法用在全变分图像恢复问题中, 数值表现较好.
此外, 问题(1.1)有着特殊的结构可以利用.因此, 许多高效的数值算法[5-10]提了出来.特别地, 即$ Q(x, y) = f(x)+\langle Bx, y\rangle-g(y). $当$ f(x) $是二次数据拟合项, $ g(y) $为示性函数, $ B $是梯度算子时, 全变分图像恢复模型可以看作是问题(1.1)的一个特殊形式.著名的原始–对偶混合梯度算法(PDHG)在文献[5]中提了出来.该方法虽然数值表现良好, 但收敛性还没被确定.随后, 许多学者开始关注此类型算法. He等学者在文献[6]中证明了如果目标函数中有一个函数是强凸的, 那么PDHG是收敛的.如果算法中的步长序列满足一定的条件, Bonettini等在文献[7]中亦证明了PDHG的收敛性. Chambolle和Pock在文献[8]中对PDHG算法进行了改进, 并得到了收敛性结论和次线性收敛速率.利用相关的不动点理论, Chen等在文献[9]中提出了一个原始–对偶不动点算法.当$ X = R^n $时, Zhang等在文献[10]中给出了一个固定步长的原始–对偶方法且原始变量的步长与对偶变量的步长不同, 算法描述如下:
其中$ \beta, \gamma > 0 $为步长.可以看出该方法每步都有闭示解, 子问题无需内迭代求解及矩阵求逆, 适合应用于计算机断层成像和并行磁共振成像等大规模问题上.
可以看到算法(1.2)和(1.3)区别在于步长的选取不同.算法(1.2)中, 每次迭代的步长相同由线搜索确定.在大规模问题中线搜索可能会增加一些额外的计算时间.算法(1.3)可以取不同的固定步长, 但与(1.2)相比, 它不能有效地处理变量全有约束的鞍点问题.本文希望结合两个算法的优点, 给出一个新算法.结合投影算子的性质和凸分析的相关知识, 给出算法的收敛性分析及收敛速率.最后, 将所设计的算法用于含有泊松噪音的图像恢复问题上, 给出相应的数值实验结果.
本节给出一个新的外梯度方法, 并分析算法的收敛性和收敛速率.结合算法(1.2)和(1.3), 一个简单的思路是用固定步长代替算法(1.2)中线搜索确定的步长, 以期待减少算法的计算时间, 且能处理约束情形下的问题(1.1).提出的方法形式简单, 有如下的迭代步骤
其中$ \beta, \gamma > 0 $满足一定的条件, 这将会在下面给出具体的说明.
接下来的部分, 具体分析新方法(2.1)产生的迭代序列收敛于问题(1.1)的鞍点.在证明之前, 列举需要用到的假设及投影算子的一些性质, 可在相关的文献[11, 12]中查到.
假设2.1 存在一个常数$ L $, 使得下面三个不等式成立
引理2.1 令$ \Omega $为非空闭凸集, $ w, z\in R^n $, $ u\in \Omega $.
(ⅰ)投影算子是非扩张的, 即$ \|P_\Omega(w)-P_\Omega(z)\|\leq\|w-z\|. $
(ⅱ)令$ G(z) $为凸可微函数, $ w = P_\Omega(z-\theta\nabla G(z)) $, 那么下式成立
(ⅲ) $ \|w-u\|^2\leq\|z-u\|^2-\|w-z\|^2+2\theta\langle\nabla G(w)-\nabla G(z), w-z\rangle+2(G(u)-G(w)). $
基于以上假设和引理, 下面给出算法(2.1)的收敛性证明.
定理2.1 如果$ \beta, \gamma $满足$ \frac1{2\beta}-L-\gamma L^2 > 0 $和$ \frac1{2\gamma}-L > 0 $, 那么新方法(2.1)生成的迭代序列$ \{(x_k, y_k)\} $收敛于优化问题(1.1)的鞍点$ (x^*, y^*) $.
证 对式(2.1)中的第三个式子应用引理2.1中的性质(ⅱ), 对任意的$ y\in Y $, 可得
同理, 令$ y = y_{k+1} $, 由(2.1)中的第一个式子得到
在式(2.3)中的右式加上和减去$ 2\gamma\langle \nabla_{y}Q(x_{k+1}, y_k), y_{k+\frac12}-y_{k+1}\rangle $, 整理后有
式(2.2)和式(2.4)相加后下式成立
又因$ Q $关于变量$ y $是凹函数, $ \langle\nabla_yQ(x_{k+1}, y_k), y-y_k\rangle\geq Q(x_{k+1}, y)-Q(x_{k+1}, y_k). $故式(2.5)的右边加上和减去$ 2\gamma\langle\nabla_yQ(x_{k+1}, y_k), y_k\rangle $得
现在考虑式(2.1)中的第二个式子, 利用引理2.1中的性质(ⅲ), 对任意的$ x\in X $, 有
(2.6)和(2.7)式相加可得
由$ Q $关于变量$ y $为凹函数, 可知
不等式(2.8)进一步转化为
利用Cauchy-Schwartz不等式(2.9)写为
由引理2.1中的性质(ⅰ)知
再由假设2.1和不等式(2.11)得到下式
令$ x = x^* $, $ y = y^* $, 且由鞍点性质$ Q(x^*, y_{k+\frac12})-Q(x_{k+1}, y^*)\leq0 $和$ \beta, \gamma $满足$ \frac1{2\beta}-L-\gamma L^2 > 0 $和$ \frac1{2\gamma}-L > 0 $, 可知
由文献[13]的结论可证序列$ \{(x_k, y_k)\} $收敛于优化问题(1.1)的鞍点$ (x^*, y^*) $.
接下来分析算法(2.1)的次线性收敛速率.为了确定算法的复杂性, 需要一些额外的记号.令$ N\geq1 $为正整数, $ x^N = \frac{1}{N+1}\sum\limits_{k = 0}^Nx_{k+1} $, $ y^N = \frac{1}{N+1}\sum\limits_{k = 0}^Ny_{k+\frac12} $.
定理2.2 如果$ \beta, \gamma $满足$ \frac1{2\beta}-L-\gamma L^2 > 0 $和$ \frac1{2\gamma}-L > 0 $, 那么
证 由式(2.12)可得到
对式(2.14)从$ k = 0, 1, \cdots, N $相加, 那么
再由函数$ Q(x, y) $的凸凹性和Jensen不等式可知
结合式(2.15)和(2.16), 易知要证的结论成立.
本节给出泊松噪音下图像去模糊的数值表现.令$ g\in R^m $是观测数据, 每个观测值$ g_i $由含有泊松随机变量的$ (Hx+b)_i $期望值确定, 其中$ x $是原始图像, $ H $可认为是采集系统造成的失真模糊矩阵, $ b\in R^m $是一个正的常数背景项.在泊松噪音的假设下, Kullback–Leibler函数经常用来作为数据拟合项.由于问题的病态性, 选取能保持图像边缘的全变分[14]作为正则项.令外, 也可以增加一些物理约束, 如像素的非负性.综上, 图像去模糊问题转化为如下的最优化问题
其中$ (\nabla x)_i $是$ x $在像素$ i $处的梯度.如果$ g_i = 0 $, 那么$ g_i\text{ln}g_i = 0 $.
模型(3.1)是非光滑凸优化问题, 直接求解有一定的难度.利用全变分函数的对偶公式, 可转化为光滑的鞍点优化问题
易知模型(3.2)是问题(1.1)的特殊形式, 且关于原始变量$ x $为凸的, 关于对偶变量$ y $是凹的.当$ Q(x, y) $由式(3.2)定义, 其梯度表达式为$ \nabla_xQ(x, y) = e_n-H^TZ(x)^{-1}g-\alpha\nabla\cdot y, \; \; \nabla_yQ(x, y) = \alpha\nabla x, $其中$ e_n = (1, 1, \cdots, 1)^T\in R^n $, $ \nabla\cdot $是散度算子; $ Z(x) $是对角矩阵, 对角元素由$ (Hx+b)_i $给定.当$ b\neq0 $时, 由文献[12]中引理4.6知$ \nabla_xQ(x, y) $李布希兹连续, 可知该问题满足假设2.1.故所提算法可以应用于问题(3.2).
记所提算法为算法1, 对比算法为文献[4]中的交替外梯度方法(AEM), 可在www.unife.it /prin/software下载, 以及算法(1.3), 按照文献[10]建议记为LPD.计算环境如下, Matlab版本为R2011a, 内存为2GB, 处理器为i3的个人电脑.在两个实验中, 图像分别为$ 128\times128 $的micro数据和$ 256\times256 $的circles数据, 详细可见文献[4]的说明.背景数据项均为$ b = 3 $.根据文献[12]知, 李布希兹常数为$ \frac{\|g\|_{\infty}\|H\|^2}{b^2} $.由此可计算出$ L $在两个问题中的值分别为0.14和0.1.据此, 测试了满足定理2.1条件下的不同参数, 选取了恢复效果比较好的参数, 具体如下:在第一个实验中, 算法1参数设置为$ \alpha = 0.09 $, $ \beta = 1.2 $, $ \gamma = 3.5 $; 第二个实验中, 参数设置为$ \alpha = 0.25 $, $ \beta = 1.3 $, $ \gamma = 1.6 $.其它两个算法按其文章里的建议设置.三个算法的停止准则为
为了评价算法恢复图像的质量, 采用相对误差(RE)作为指标, 定义为$ RE = \frac{\|x_k-x\|}{\|x\|}, $其中$ x_k $为算法恢复的图像, $ x $为原始图像.它能衡量算法恢复的图像与真实图像的接近程度.显然, 相对误差的值越小, 图像恢复的质量越好.
表 1给出了算法在两个含有泊松噪音数据上的去模糊表现.在micro问题中, 虽然算法1的迭代步数多于AEM, 但发现计算时间稍微少于AEM.这说明由线搜索确定的步长在一定程度上会消耗一点时间. LPD算法无论是迭代步数, 计算时间和相对误差均不占优势.在circles问题中, 无论迭代步数和计算时间都少于AEM. LPD方法虽然相对误差较小, 但计算步数较多, 计算时间较长.图 1和图 2显示了两个算法恢复的效果图, 可以看到算法基本上都能很好的恢复模糊的图像.这也可以从表 1的评价标准相对误差中得到, 因为相对误差的值相差不大.为更清楚地说明算法的表现, 在图 3中画出了目标函数值随计算时间的曲线图.从曲线下降趋势看, 算法1在迭代的前几步或十几步函数值下降快于AEM, 随后两个算法函数值基本保持不变.表明算法1可以更快的收敛于最优解.虽然算法LPD的目标函数值在micro问题中下降较快, 但它是在目标函数(3.2)在$ X = R^n $的情况下.由于无相关物理条件下的约束, 故其相对误差较大.这也说明了如果有先验知识的加入, 如非负约束, 则可能得到较好的恢复结果.最后, 和相关算法对比的数值结果表明了新方法是有效的.
对带有凸集约束的光滑鞍点优化问题, 提出了一个简单的原始–对偶梯度方法.算法每步具有显示解, 不需要内迭代求解子问题.新方法应用在全变分正则化泊松噪音图像恢复问题上, 结果表明在计算时间上具有一定的优势.值得注意的是该方法也可以拓展求解其它噪音类型的图像恢复问题.