一类求解约束离散不适定问题的积极集随机迭代方法

摘要

许多科学和工程领域的应用问题都可以归结为线性离散不适定问题的求解。考虑大规模带盒子约束的线性离散不适定问题的求解，提出一类基于积极集策略的随机内外迭代方法。基于积极集策略的内外迭代法在外层迭代上更新积极集和对应的非积极集，并采用投影算子，将不在可行域中的数值解分量投影到可行域边界上，同时在内层迭代上采用Krylov子空间方法求解无约束子问题。提出一类积极集迭代法，在内层迭代上采用高性能随机算法，依照概率分布选取子问题系数矩阵的列进行更新，并利用Armijo下降准则对迭代步长进行选择，这样就可以保证目标函数值随着迭代步数的增加而单调下降。在图像复原问题的数值实验中，验证所构造算法的高效性。在偏差准则的收敛条件下，新的积极集内外迭代法所利用的计算量、迭代步数和CPU时间都比前人提出的算法更少。

关键词

线性离散不适定问题; 盒子约束; 随机迭代; 积极集方法; 内外迭代算法

线性离散不适定问题来源于许多科学计算和工程应用领域^［1］。例如第一类Fredholm积分方程数值离散后的问题；例如图像去模糊复原问题，已知观测到的模糊带噪声图像和模糊算子，需要尽可能复原得到清晰图像。也就是说，在数学上考虑一个不相容的线性方程组

$$\mathit{A}\mathit{x}=\mathit{b}$$

（1）

其中$\mathit{A}\in {\mathbf{R}}^{m\times n}$是一个大型稀疏的矩阵，表示离散的模糊算子，可以是方阵，也可以是长方形矩阵。右端项$\mathit{b}\in {\mathbf{R}}^{m\times \mathrm{1}}$表示观测得到的数据，即带有噪声的模糊图像。式（1）和通常的线性方程组不同之处在于，矩阵A是长方形矩阵，且往往行数m大于列数n，因此式（1）是一个超定的线性系统。从数值角度上分析，矩阵A来源于带积分项的模糊算子的离散，往往大型稀疏，且奇异值聚集在原点附近。这样的矩阵往往带有坏条件数，甚至接近奇异矩阵。观测的数据b往往带有观测误差。

$$\mathit{b}=\widehat{\mathit{b}}+\mathit{e}$$

（2）

式中：$\widehat{\mathit{b}}$为未知的没有噪声的向量，且$\widehat{\mathit{b}}=\mathit{A}\widehat{\mathit{x}}$，$\widehat{\mathit{x}}$为未知的清晰的真实图像；向量e表示某种误差或噪声，其中每个分量都独立服从相同的随机分布。因此，线性系统（1）是不相容的，只能求解其对应的最小二乘解。

令${\mathit{A}}^{†}$为矩阵A的摩尔-彭若斯广义逆，则不适定问题（1）的最小二乘解为

$$\mathit{x}={\mathit{A}}^{†}\mathit{b}=\widehat{\mathit{x}}+{\mathit{A}}^{†}\mathit{e}$$

然而，由于不适定问题中矩阵A病态，式（1）的解对右端项b中的误差e非常敏感，由噪声引起的分量${\mathit{A}}^{†}\mathit{e}$会将所需要的信息$\widehat{\mathit{x}}$掩盖，使得解${\mathit{A}}^{†}\mathit{b}$没有意义。为了解决此问题，需要正则化方法。常用的正则化方法分为2类：当噪声e的取值范围已知时，可以利用偏差原理（discrepancy principle）方法；当噪声e未知时，可以利用Tikhonov正则化方法、广义交叉校验准则和L曲线准则等确定正则化参数。本文采用Tikhonov正则化方法，即最小化带有惩罚项的目标函数

$$\mathrm{m}\mathrm{i}\mathrm{n}f(x)={‖\mathit{A}\mathit{x}-\mathit{b}‖}_{\mathrm{2}}^{\mathrm{2}}+{\mu }^{\mathrm{2}}{‖Lx‖}_{\mathrm{2}}^{\mathrm{2}}$$

（3）

式中：L为一个离散微分算子；$\mu$为正则化参数，$\mu >\mathrm{0}$。

在实际应用中，数据往往是在一定的取值范围内选取的。例如，图像的像素值是非负的；灰度图像的像素值则在［0，255］之间。因此，考虑一类带盒子约束的不适定问题^［2］，即

$$\mathrm{m}\mathrm{i}\mathrm{n}f(\mathit{x})$$ 使得 $$\mathit{l}\le \mathit{x}\le \mathit{u}$$

（4）

其中，$\mathit{l}\in {\mathbf{R}}^{n\times \mathrm{1}}$和$\mathit{u}\in {\mathbf{R}}^{n\times \mathrm{1}}$分别是下界向量和上界向量。

针对大规模带约束不适定问题（4），Bertero和Boccacci^［3］提出一类投影Landweber方法。Hanke等^［4］考虑变量无上界的特殊情形，利用显性变量替换将原问题转换成一类非线性方程组，然后采用拟牛顿法进行求解。Nagy和Strakos^［5］将最速下降法与投影技巧结合，在目标函数下降的同时，迭代解满足约束条件。Rojas和Steihaug^［6］从优化算法出发，提出惩罚函数法，使得数值解序列在可行区域内部收敛到最优点。Calvetti等^［2］将正则化方法和Krylow子空间迭代法相结合，提出一类内外迭代法。文献［6-7］在此基础上引入积极集策略，外层迭代快速更新积极集，而内层迭代则在子空间通过迭代求解一系列仅在非积极集上有约束的问题。数值实验验证此算法的有效性，并且数值效率高于前人提出的算法，所需的计算复杂度和CPU时间都更少。但此方法的缺点在于缺乏理论收敛性，尤其是无法保证目标函数值随着迭代步数逐渐单调下降。

针对内层迭代中需要求解的无约束最小二乘问题，近年来诸多学者提出高性能随机迭代方法。Strohmer和Vershynin^［8］提出随机Kaczmarz方法在求解相容系统上的收敛理论。Leventhal和Levis^［9］提出求解不相容系统，即最小二乘问题的随机坐标下降方法，并给出当系数矩阵列满秩时的收敛理论。随着新的有效概率准则的构造，更多高性能的随机迭代方法相继被提出^［9-11］，在理论和数值实验中均展示出比传统确定性方法更高效的计算性能。

本文提出一类基于积极集策略的内外迭代法。该算法在内层迭代上利用随机坐标下降方法求解无约束最小二乘问题，并利用Armijo下降准则对迭代步长进行选择，这样就可以保证目标函数值随着迭代步数的增加而单调下降；同时在外层迭代上更新积极集和对应的非积极集，并采用投影算子，将不在可行域中的数值解分量投影到可行域边界上。数值实验验证此算法的高效性。在偏差准则的收敛条件下，新的积极集内外迭代法所利用的计算量、迭代步数和CPU时间都比前人提出的算法更少。

1 基于积极集策略的内外迭代法

简单介绍求解问题（4）的2种内外迭代法。第1种方法是基于在内迭代中将数值解投影到盒子约束可行域$\Omega =\{x\in {\mathbf{R}}^n:l\le x\le u\}$中。第2种方法是基于积极集策略，通过在每次外迭代中更新积极集，从而在内迭代中只需要求解一个规模较小的子问题。

容易验证，问题（4）中的目标函数等价于

$$f(x)={‖\left[\begin{array}{c}\mathit{A}\\ \mu L\end{array}\right]x-\left[\begin{array}{c}\mathit{b}\\ \mathrm{0}\end{array}\right]‖}_{\mathrm{2}}^{\mathrm{2}}={‖\tilde{\mathit{A}}\mathit{x}-\tilde{\mathit{b}}‖}_{\mathrm{2}}^{\mathrm{2}}$$

如果不考虑对变量的盒子约束，极小化$f(x)$在数学上等价于求解一个正规方程组

$${\tilde{\mathit{A}}}^{\mathrm{T}}\tilde{\mathit{A}}\mathit{x}={\tilde{\mathit{A}}}^{\mathrm{T}}\tilde{\mathit{b}}$$

令$\widehat{\mathit{x}}$为此方程组的解，则通常$\widehat{\mathit{x}}\notin \Omega$。定义一个正交投影算子$P_{\Omega }$ ，则$x=P_{\Omega }(\widehat{x})$为

$$x_j=\left\{\begin{array}{c}{l}_j,\mathrm{当}{\widehat{x}}_j<l_j\\ u_j,\mathrm{当}{\widehat{x}}_j>u_j\\ {\widehat{x}}_j,\mathrm{当}{l}_j\le {\widehat{x}}_j\le u_j\end{array}\right.$$

其中$j$ =1，2，3…，$n$，这样可以保证$x\in \Omega$。通过反复迭代，使得数值解对应的目标函数值逐渐下降，从而得到所需要的解。

投影内外迭代算法^［3］求解带盒子约束的问题（4）如下：

算法1  ^［3］     投影内外迭代法。①选取一个迭代向量${\mathit{x}}^{\mathrm{0}}$ 。②计算$r^{\mathrm{0}}=\tilde{b}-\tilde{\mathit{A}}{\mathit{x}}^{\mathrm{0}}$ 。③开始迭代k=0，1，2，…，直到收敛。④通过迭代法求得如下问题的数值解$w^k$：

$$\underset{}{\underset{w\in R^{\breve{n}}}{\mathrm{m}\mathrm{i}\mathrm{n}}{‖\tilde{\mathit{A}}w-r^k‖}_{\mathrm{2}}}$$

（5）

⑤令${\widehat{\mathit{x}}}_{}^{k+\mathrm{1}}={\mathit{x}}_{}^k+w^k$并投影到可行区域中${\mathit{x}}^{k+\mathrm{1}}=P_{\Omega }({\widehat{\mathit{x}}}^{k+\mathrm{1}})$。⑥计算$r^{k+\mathrm{1}}=\tilde{\mathit{b}}-\tilde{\mathit{A}}{\mathit{x}}^{k+\mathrm{1}}$。⑦结束迭代。其中，最小二乘问题（5）的求解为内层迭代，k的更新为外层迭代，因此算法1为内外迭代法。对于式（5）的求解，可以利用CGLS算法^［1］，即共轭梯度法求解其正规方程组：${\tilde{\mathit{A}}}^{\mathbf{T}}\tilde{\mathit{A}}w=\tilde{\mathit{A}}{r}^k$。如果选取内迭代的初始向量为零向量，则得到的数值解$w^k$在Krylov子空间$K({\tilde{\mathit{A}}}^{\mathrm{T}}\tilde{\mathit{A}},{\tilde{\mathit{A}}}^{\mathrm{T}}\mathit{b})$=span｛${\tilde{\mathit{A}}}^{\mathrm{T}}\mathit{b},{\tilde{\mathit{A}}}^{\mathrm{T}}\tilde{\mathit{A}}{\tilde{\mathit{A}}}^{\mathrm{T}}\mathit{b},({\tilde{\mathit{A}}}^{\mathrm{T}}{\tilde{\mathit{A}})}^{\mathrm{2}}{\tilde{\mathit{A}}}^{\mathrm{T}}\mathit{b}$，…｝中。

Morigi等^［7］在投影内外迭代算法2的基础上引入积极集策略。令拉格朗日算子${\mathit{\lambda }}^k=-{\tilde{\mathit{A}}}^{\mathrm{T}}{\mathit{r}}^k$，其中${\mathit{x}}^k$ 是第k步的数值解，${\mathit{r}}^k=\tilde{\mathit{b}}-\tilde{\mathit{A}}{\mathit{x}}^k$是${\mathit{x}}^k$对应的残量。根据问题（4）的Karush-Kuhn-Tucker条件可知，${\mathit{x}}^{*}$是问题（4）的最优解当且仅当其满足

$$x_j^{*}\ge l_j,{\lambda }_j\ge \mathrm{0},{\lambda }_j(x_j^{*}-l_j)=\mathrm{0}$$

（6）

$$x_j^{*}\le u_j,{\lambda }_j\le \mathrm{0},{\lambda }_j(u_j-x_j^{*})=\mathrm{0}$$

（7）

对于$x^k$，定义积极集为$B(x^k)=\{j:x_j^k=l_j,{\lambda }_j^k\ge \mathrm{0}\}\bigcup \{j:x_j^k=u_j,{\lambda }_j^k\le \mathrm{0}\}$，其对应的补集自由变量为$F(x^k)=\{\mathrm{1,2},...,n\}\backslash B(x^k)$。基于积极集算法的思想是：如果已知$x^{*}$的积极集，则约束最优化问题（4）可以转化成无约束优化问题$\mathrm{m}\mathrm{i}\mathrm{n}{‖{\tilde{\mathit{A}}}_F{\mathit{x}}_F-\stackrel{-}{\mathit{b}}‖}_{\mathrm{2}}^{\mathrm{2}}$ 使得 $l_F\le x_F\le u_F$，其中$\stackrel{\_}{\mathit{b}}=\tilde{\mathit{b}}-{\tilde{\mathit{A}}}_B{\mathit{x}}_B$，${\tilde{\mathit{A}}}_F$是矩阵$\tilde{\mathit{A}}$中列指标属于集合$F$的子矩阵。由于$x^{*}$未知，所以先选取初始迭代向量不断更新积极集，直到收敛为止。

基于积极集的投影内外迭代算法^［5］求解带盒子约束的问题（4）如下。

算法2  ^［5］     基于积极集的投影内外迭代法。①选取一个迭代初始向量${\mathit{x}}^{\mathrm{0}}$。②计算${\mathit{r}}^{\mathrm{0}}=\tilde{\mathit{b}}-\tilde{\mathit{A}}{\mathit{x}}^{\mathrm{0}}$。③开始迭代k=0，1，2，…，直到收敛。④计算拉格朗日算子${\lambda }^k=-{\tilde{\mathit{A}}}^{\mathrm{T}}{\mathit{r}}^k$。⑤更新积极集$B(x^k)$和自由变量集$F(x^k)$。⑥通过迭代法求得如下问题的数值解$w^k$：$\underset{}{\underset{w\in {\mathbf{R}}^{\breve{n}}}{\mathrm{m}\mathrm{i}\mathrm{n}}{‖{\tilde{A}}_{F}w-r^k‖}_{\mathrm{2}}}$。⑦令${\widehat{x}}_F^{k+\mathrm{1}}=x_F^k+w^k$，${\widehat{x}}_B^{k+\mathrm{1}}=x_B^k$，并投影到可行区域中，即$x^{k+\mathrm{1}}=P_{\Omega }({\widehat{x}}^{k+\mathrm{1}})$。⑧计算${\mathit{r}}^{k+\mathrm{1}}=\tilde{\mathit{b}}-\tilde{\mathit{A}}{\mathit{x}}^{k+\mathrm{1}}$。⑨结束迭代。

在算法2中，外层迭代快速更新积极集，而内层迭代则在子空间通过迭代求解一系列仅在非积极集上有约束的问题。和算法1类似，内层迭代可以采用CGLS算法快速求解。数值实验^［5］表明，基于积极集的算法收敛速度更快，且所需的迭代步数和CPU时间都更少。

算法2的缺点在于理论上无法保证目标函数值随着迭代单调下降。因此，在第2节中构造新的基于积极集的算法。

2 结合Armijo策略的新积极集算法

将Armijo充分下降准则和积极集思想相结合，同时，在内层无约束最小二乘问题的求解上采用随机坐标下降策略。

算法3  　基于积极集和Armijo条件的随机投影内外迭代法。①选取一个迭代初始向量${\mathit{x}}^{\mathrm{0}}$。②计算${\mathit{r}}^{\mathrm{0}}=\tilde{\mathit{b}}-\tilde{\mathit{A}}{\mathit{x}}^{\mathrm{0}}$。③开始迭代$k=\mathrm{0,1},\mathrm{2},...,$ 直到收敛。④计算拉格朗日算子${\lambda }^k=-{\tilde{\mathit{A}}}^{\mathrm{T}}{\mathit{r}}^k$。⑤更新积极集$B(x^k)$和自由变量$F(x^k)$。⑥通过随机坐标下降迭代法求得如下问题的数值解$w^k$：$\underset{}{\underset{w\in {\mathbf{R}}^{\breve{n}}}{\mathrm{m}\mathrm{i}\mathrm{n}}{‖{\tilde{A}}_{F}w-r^k‖}_{\mathrm{2}}}$。⑦令${\widehat{x}}_F^{k+\mathrm{1}}=x_F^k+{\beta }^m{w}^k,{\widehat{x}}_B^{k+\mathrm{1}}=x_B^k$，其中m是满足如下Armijo充分下降准则的最小非负整数：

$${‖\tilde{\mathit{A}}{\mathit{x}}^{k+\mathrm{1}}-\tilde{\mathit{b}}‖}_{\mathrm{2}}^{\mathrm{2}}\le {‖\tilde{\mathit{A}}{\mathit{x}}^k-\tilde{\mathit{b}}‖}_{\mathrm{2}}^{\mathrm{2}}+\mathrm{2}\alpha ({\mathit{\lambda }}^k{)}^{\mathrm{T}}({\mathit{x}}^{k+\mathrm{1}}-{\mathit{x}}^k)$$

（8）

并投影到可行区间中

$$x^{k+\mathrm{1}}=P_{\Omega }({\widehat{x}}^{k+\mathrm{1}})$$

⑧计算$r^{k+\mathrm{1}}=\tilde{b}-\tilde{A}{x}^{k+\mathrm{1}}$。⑨结束迭代。在算法2的步骤⑦中，参数$\mathrm{0}\le \alpha <\mathrm{1},\mathrm{0}<\beta <\mathrm{1}$。

3 数值实验

用图像去模糊复原的数值实验来验证基于积极集和Armijo条件的投影内外迭代法求解带盒子约束的不适定问题的可行性和有效性。

在数值实验中，采用的收敛准则为偏差原理。另外，选取二阶微分算子为正则化光滑子

$$L=\left[\begin{array}{cccccc}\mathrm{1}& -\mathrm{2}& \mathrm{1}& & & \\ & \mathrm{1}& -\mathrm{2}& \mathrm{1}& & \\ & & \ddots & \ddots & \ddots & \\ & & & \mathrm{1}& -\mathrm{2}& \mathrm{1}\end{array}\right]$$

在图1中分别画出图像复原算例“satellite”的原始图像、模糊带噪声图像、算法2复原图像、算法3复原图像。图2、图3分别为相对误差和相对残量随着迭代步数变化的曲线。从图2可以看出，2类积极集内外迭代算法随着迭代步数的增加，相对误差都先下降后增加。最优外迭代步数在5步附近。从图3可以看出，相比于算法2，算法3的相对残量下降更快，振荡更小，这说明随机坐标下降方法结合Armijo下降准则可以使计算更快更平稳地收敛，所需要的计算量更少。

图1 图像复原算例“satellite”

Fig. 1 Restored image of ‘statellite’

图2 图像复原算例“satellite”相对误差随迭代步数变化

Fig. 2 Relative error of restored image of ‘satellite’ versus the number of iterative steps

图3 图像复原算例“satellite”相对残量随迭代步数变化

Fig. 3 Relative residue of restored image ‘satellite’ versus the number of iteration steps

4 结论和展望

考虑带盒子约束的不适定问题，构造一类基于积极集策略的内外迭代法。此算法的创新之处在于，在内层迭代上利用随机坐标下降方法，同时采用Armijo下降准则对迭代步长进行选择，这样就可以保证目标函数值随着迭代步数的增加而单调下降。同时在外层迭代上更新积极集和对应的非积极集，并采用投影算子，将不在可行域中的数值解分量投影到可行域边界上。数值实验验证此算法的高效性。在偏差准则的收敛条件下，新的积极集内外迭代法所用的计算量、迭代步数和CPU时间都比前人提出的算法更少。在后续工作中，将考虑稀疏约束条件下的最小二乘问题，并构造类似的积极集随机迭代算法进行求解。

作者贡献声明

郑 宁：负责算法构造、数值实验和分析、论文撰写等。

殷俊锋：负责算法设计和数值案例分析。