0 引言

【研究意义】在理论研究和工程实践中,许多问题的解决需要根据观测数据来恢复或重构物理介质的结构信息,通常称这类问题为逆问题.逆问题分为线性逆问题和非线性逆问题,其中线性逆问题广泛出现在地球物理学、图像处理、地下水文学、全球海洋模拟、肿瘤检测、非破坏性检测等领域中^[1].【前人研究进展】对于线性逆问题的研究，已有不少数值解法，文献^[2]提出迭代Tikhonov正则化方法求解该问题，取得了较好的数值效果.文献^[3]提出了Levenberg-Marquardt方法求解该问题，且在数值试验中，用偏差原则作为停止准则，但该方法所需CPU时间较长.文献^[4]提出并讨论了一类NCPs逆问题.【本研究切入点】谱共轭梯度法是一类非常重要的优化方法，他是将谱梯度法和共轭梯度法的思想相结合的一种方法，具有迭代简单、易于编程以及存储要求低等优点，所以一般适用于求解大规模无约束优化问题.近年来，文献^[5-7]用谱共轭梯度法求解线性和非线性优化问题，得到很好的数值效果.因此，本文希望把谱共轭梯度法应用到线性逆问题中.【拟解决的关键问题】针对线性逆问题，提出一个新的谱共轭梯度法，在理论上证明新提共轭梯度法的下降性和全局收敛性，并把该方法应用到卷积方程反演和图像处理的反演中.

1 新的谱共轭梯度法

在许多科学研究和工程应用中，很多问题可以转化为下面形式的算子方程

$Hx=u.$

(1)

我们研究的问题是已知H和u反解x，这里H表示线性算子，u为实验测得的数据，因而不可避免地带有一定的噪声^[1]，即有

${{u}_{e}}=u+\eta ,$

(2)

这里，η表示与u维数相同的噪声。即问题(1)转化为求解

$Hx={{u}_{e}}.$

(3)

为了便于求解，将问题(3)等价转化为求解下面的最小二乘问题:

$min~f\left( x \right)=\|Hx-{{u}_{e}}{{\|}^{2}}.$

(4)

由于H是一个条件数很大的病态矩阵，所以问题(3)是一个不适定问题，因此没有唯一解.为了克服问题(3)的不适定性，本文在式(4)的基础上通过增加一个Tikhonov正则项‖x‖²，将不适定问题转化为适定问题，即求解下面的无约束优化问题：

$min~f\left( x \right)=\|Hx-u{{\|}^{2}}+\alpha \|x{{\|}^{2}},$

(5)

其中，α为正则化参数，而且该参数的选取对问题的求解发挥重要作用.

对于无约束优化问题:

$min~\{f\left( x \right)|x\in {{R}^{n}}\},$

(6)

其中f(x)为式(5)的目标函数.共轭梯度法是求解问题(5)的有效方法^[8-10]，其标准迭代形式为

${{x}_{k+1}}={{x}_{k}}+{{\alpha }_{k}}{{d}_{k}},$

(7)

其中α_k是Wolfe搜索产生的步长，d_k为搜索方向.d_k定义为

${{d}_{k}}=\left\{ \begin{matrix} -{{g}_{k}},k=1, \\ -{{g}_{k}}+{{\beta }_{k}}{{d}_{k-1}},k\ge 2, \\ \end{matrix} \right.$

(8)

其中g_k=Δf(x_k)，β_k为方向调控参数.不同的β_k对应不同的共轭梯度法，其中著名的β_k 计算公式^[11-12]有

$\beta _{k}^{~CD}=-\frac{g_{k}^{T}{{g}_{k}}}{d_{k-1}^{T}{{g}_{k-1}}},\beta _{k}^{LS}=-\frac{g_{k}^{T}({{g}_{k}}-{{g}_{k-1}})}{g_{k-1}^{T}{{d}_{k-1}}}.$

(9)

上述由参数β_k构成的方法分别称为CD方法和LS方法.

Birgin等^[13]结合谱梯度的思想，提出了下面形式的谱共轭梯度法：

${{d}_{k}}=\left\{ \begin{matrix} -{{g}_{k}},k=1, \\ -{{\theta }_{k}}{{g}_{k}}+{{\beta }_{k}}{{d}_{k-1}},k>1, \\ \end{matrix} \right.$

(10)

并取

${{\beta }_{k}}=\frac{{{\left( {{\theta }_{k}}{{y}_{k-1}}-{{s}_{k-1}} \right)}^{T}}{{g}_{k}}}{d_{k-1}^{T}{{y}_{k-1}},}$

(11)

${{\theta }_{k}}=\frac{s_{k-1}^{T}{{s}_{k-1}}}{s_{k-1}^{T}{{y}_{k-1}}},$

(12)

其中y_k－1=g_k－g_k－1，s_k－1=x_k－x_k－1.基于迭代格(10)，本文给出新的β_k和θ_k(文献^[14])分别为

${{\beta }_{k}}=\frac{\|{{g}_{k}}{{\|}^{2}}+1}{d_{k-1}^{T}({{g}_{k}}-{{g}_{k-1}})},$

(13)

$\frac{\left( 1+\frac{1}{{{\left\| {{g}_{k}} \right\|}^{2}}} \right)g_{k}^{T}{{d}_{k-1}}-\left( 1+\frac{2}{{{\left\| {{g}_{k}} \right\|}^{2}}} \right)g_{k-1}^{T}{{d}_{k-1}}}{d_{k-1}^{T}\left( {{g}_{k}}-{{g}_{k-1}} \right)}.$

(14)

称式(10),(13),(14)为新的谱共轭梯度法(SCG).以下结合Wolfe线搜索建立新的谱共轭梯度算法，并分析算法的收敛性质.

算法1

步骤1 给定初值x₁∈Rⁿ,0<δ<1/2<σ<1,d₁:=－g₁,ε≥0.如果‖g₁‖≤ε,则停止迭代.

步骤2 由Wolfe线搜索准则计算步长α_k，即α_k满足：

$f\left( {{x}_{k}}+{{\alpha }_{k}}{{d}_{k}} \right)-{{f}_{k}}\le \delta {{\alpha }_{k}}g_{k}^{T}{{d}_{k}},$

(15)

$g{{\left( {{x}_{k}}+{{\alpha }_{k}}{{d}_{k}} \right)}^{T}}{{d}_{k}}\ge \sigma g_{k}^{T}{{d}_{k}}.$

(16)

步骤3 计算x_k+1=x_k+α_kd_k.若‖g_k+1‖≤ε，则停止迭代.

步骤4 由式(13)计算β_k+1，由式(14)计算θ_k+1，由d_k+1=－θ_k+1g_k+1+β_k+1d_k计算d_k+1.

步骤5 k:=k+1,转步骤2.

以下均假设‖g_k‖≠0，否则算法因已找到稳定点而停止.

引理1 设{g_k}和{d_k}为算法1生成的序列，则

$g_{k}^{T}{{d}_{k}}<0,\forall k\ge 1.$

(17)

证明 用数学归纳法证明.当k=1时，g₁^Td₁=g₁^T(－g₁)=－‖g₁‖²<0.现假设对于k－1的情形有g^T_k－1d_k－1<0成立，由式(16)可得

$d_{k-1}^{T}({{g}_{k}}-{{g}_{k-1}})\ge \left( \sigma -1 \right)g_{k-1}^{T}{{d}_{k-1}}>0.$

(18)

则对于k的情形就有

$\begin{align} & g_{k}^{T}{{d}_{k}}=g_{k}^{T}\left( -{{\theta }_{k}}{{g}_{k}}+{{\beta }_{k}}{{d}_{k-1}} \right)= \\ & \frac{-(1+\|{{g}_{k}}{{\|}^{2}})g_{k}^{T}{{d}_{k-1}}+(\|{{g}_{k}}{{\|}^{2}}+2)g_{k-1}^{T}{{d}_{k-1}}}{d_{k-1}^{T}({{g}_{k}}-{{g}_{k-1}})}+ \\ & \frac{(\|{{g}_{k}}{{\|}^{2}}+1)g_{k}^{T}{{d}_{k-1}}}{d_{k-1}^{T}({{g}_{k}}-{{g}_{k-1}})}\text{ }=\frac{(\|{{g}_{k}}{{\|}^{2}}+2)g_{k-1}^{T}{{d}_{k-1}}}{d_{k-1}^{T}({{g}_{k}}-{{g}_{k-1}})}<0. \\ \end{align}$

(19)

从而由数学归纳法知引理1成立.

引理1 说明算法1产生的搜索方向是下降方向.

引理2 如果步长因子α_k满足Wolfe线搜索条件,则有

$0<{{\beta }_{k}}<\frac{g_{k}^{T}{{d}_{k}}}{g_{k-1}^{T}{{d}_{k-1}}}.$

(20)

证明 由式(13)可得

${{\beta }_{k}}=\frac{\|{{g}_{k}}{{\|}^{2}}+1}{d_{k-1}^{T}({{g}_{k}}-{{g}_{k-1}})}>0.$

(21)

对式(19)两边同时除以g_k－1^Td_k－1,可得

$\frac{g_{k}^{T}{{d}_{k}}}{g_{k-1}^{T}{{d}_{k-1}}}=\frac{\|{{g}_{k}}{{\|}^{2}}+2}{d_{k-1}^{T}({{g}_{k}}-{{g}_{k-1}})}>$

(22)

$\begin{align} & \frac{\|{{g}_{k}}{{\|}^{2}}+1}{d_{k-1}^{T}\left( {{g}_{k}}-{{g}_{k-1}} \right)}={{\beta }_{k}}< \\ & \frac{g_{k}^{T}{{d}_{k}}}{g_{k-1}^{T}{{d}_{k-1}}}. \\ \end{align}$

(23)

所以，引理2得证.

2 算法的全局收敛性

为了证明算法1的全局收敛性，作如下假设^[14-16]：

(H₁)目标函数f(x)在其水平集Ω={x∈Rⁿ|f(x)≤f(x₁)}上有界.

(H₂)设f(x)的梯度g(x)在Ω上Lipschitz连续,即存在L>0，使

$\|g\left( y \right)-g\left( x \right)\|\le L\|y-x\|,\forall x,y\in \Omega .$

(24)

引理3 假设(H₁)，(H₂)成立，{g_k}和{d_k}为算法 1 生成的序列，则

$\sum\limits_{k\ge 1}{\frac{{{(g_{k}^{T}{{d}_{k}})}^{2}}}{\|{{d}_{k}}{{\|}^{2}}}}<+\infty .$

(25)

证明 由式(7)和式(16)及假设(H₂)，可得

$L{{\alpha }_{k}}\|{{d}_{k}}{{\|}^{2}}\ge d_{k}^{T}({{g}_{k+1}}-{{g}_{k}})\ge \left( \sigma -1 \right)g_{k}^{T}{{d}_{k}},$

(26)

从而得到

${{\alpha }_{k}}\ge \frac{\left( \sigma -1 \right)}{L}\frac{g_{k}^{T}{{d}_{k}}}{\|{{d}_{k}}{{\|}^{2}}}.$

(27)

由Wolfe搜索准则和d_k的下降性及假设(H₁)知，{f_k}为单调递减的收敛数列，结合式(27)可得

${{f}_{k}}-{{f}_{k+1}}\ge -\delta {{\alpha }_{k}}g_{k}^{T}{{d}_{k}}\ge \frac{\delta \left( 1-\sigma \right)}{L}\frac{{{(g_{k}^{T}{{d}_{k}})}^{2}}}{\|{{d}_{k}}{{\|}^{2}}}.$

(28)

对式(28)左右两端分别从k=1,2,…求和，再利用{f_k}的收敛性，可得式(25)成立,则引理3得证.

定理1 假设(H₁)，(H₂)成立，{g_k}为算法1生成的序列，则

$\underset{k\to \infty }{\mathop{lim}}\,~inf~\|{{g}_{k}}\|=0.$

(29)

证明 若定理1不成立，则必存在常数r>0，使得对任意的k≥1，有

$\|{{g}_{k}}\|\ge r.$

(30)

式(10)移项，得d_k+θ_kg_k=β_kd_k－1，等式两端取模平方，移项并利用式(20)，可得

$\begin{align} & \|{{d}_{k}}{{\|}^{2}}={{\left( {{\beta }_{k}} \right)}^{2}}\|{{d}_{k-1}}{{\|}^{2}}-2{{\theta }_{k}}g_{k}^{T}{{d}_{k}}- \\ & \theta _{k}^{2}\|{{g}_{k}}{{\|}^{2}}~<\frac{{{(g_{k}^{T}{{d}_{k}})}^{2}}}{{{(g_{k-1}^{T}{{d}_{k-1}})}^{2}}}\|{{d}_{k-1}}{{\|}^{2}}-2{{\theta }_{k}}g_{k}^{T}{{d}_{k}}- \\ & \theta _{k}^{2}\|{{g}_{k}}{{\|}^{2}}. \\ \end{align}$

(31)

式(31)两端同时除以(g_k^Td_k)²，可得

$\begin{align} & \frac{\|{{d}_{k}}{{\|}^{2}}}{{{\left( g_{k}^{T}{{d}_{k}} \right)}^{2}}}<\frac{\|{{d}_{k-1}}{{\|}^{2}}}{{{\left( g_{k-1}^{T}{{d}_{k-1}} \right)}^{2}}}-\frac{2{{\theta }_{k}}}{g_{k}^{T}{{d}_{k}}}-\frac{\theta _{k}^{2}\|{{g}_{k}}{{\|}^{2}}}{{{\left( g_{k}^{T}{{d}_{k}} \right)}^{2}}}~= \\ & \frac{\|{{d}_{k-1}}{{\|}^{2}}}{{{\left( g_{k-1}^{T}{{d}_{k-1}} \right)}^{2}}}+\frac{1}{\|{{g}_{k}}{{\|}^{2}}}-{{\left( \frac{1}{\|{{g}_{k}}\|}+\frac{{{\theta }_{k}}\|{{g}_{k}}\|}{\left( g_{k}^{T}{{d}_{k}} \right)} \right)}^{2}}< \\ & \frac{\|{{d}_{k-1}}{{\|}^{2}}}{{{(g_{k-1}^{T}{{d}_{k-1}})}^{2}}}+\frac{1}{\|{{g}_{k}}{{\|}^{2}}}. \\ \end{align}$

(32)

由式(32)的递推，并利用d₁=－g₁，可得

$\begin{align} & \frac{\|{{d}_{k}}{{\|}^{2}}}{{{\left( g_{k}^{T}{{d}_{k}} \right)}^{2}}}<\frac{\|{{d}_{1}}{{\|}^{2}}}{\|{{g}_{1}}{{d}_{1}}{{\|}^{2}}}+ \\ & \sum\limits_{i=2}^{k}{\frac{1}{\|{{g}_{i}}{{\|}^{2}}}}=\sum\limits_{i=1}^{k}{\frac{1}{\|{{g}_{i}}{{\|}^{2}}}}. \\ \end{align}$

(33)

结合式(30)，可得

$\frac{{{(g_{k}^{T}{{d}_{k}})}^{2}}}{\|{{d}_{k}}{{\|}^{2}}}>\frac{{{r}^{2}}}{k}.$

(34)

对式(34)两端分别求和可得

$\sum\limits_{k\ge 1}{\frac{{{(g_{k}^{T}{{d}_{k}})}^{2}}}{\|{{d}_{k}}{{\|}^{2}}}}>{{r}^{2}}\sum\limits_{k\ge 1}{\frac{1}{k}}=+\infty ,$

(35)

显然，式(35)与式(25)矛盾，所以定理1成立.

3 数值实验

算法测试的环境为Matlab 2013a，Windows 7操作系统，Intel(R)Core(TM) i3-2370M CPU @6.40 GHz.用Matlab的rand函数在真实解上添加一定水平的噪声.参数δ=0.025，σ=0.9，α=0.014.为了评价不同方法的求解效果，引入真实解和数值解的均方根误差(RMSE)的指标^[17-18]来衡量：

$RMSE=\frac{\|x-{{x}_{k}}{{\|}_{2}}}{\|x{{\|}_{2}}},$

(36)

其中，x和x_k分别是真实解和数值解，显然RMSE越小，数值解对应的算法越好.停止准则是

$\|{{g}_{k}}\|\le {{10}^{-5}}.$

(37)

例1 卷积方程的反演

考察在位势理论中有重要应用的卷积方程^[1]：

$\begin{align} & -\frac{1}{2\pi }\int_{0}^{2\pi }{(1-ln(4si{{n}^{2}}\left( t-t\prime \right)/2))x\left( t\prime \right)d\left( t\prime \right)} \\ & =u\left( t \right),0\le t\le 2\pi , \\ \end{align}$

(38)

其中u(t)与位势方程边界条件有关,弱奇异性核H(t,t′)=1－ln(4sin²(t－t′)/2).现将式(38)记为算子形式：

$\int_{0}^{2\pi }{H\left( t,t\prime \right)x\left( t\prime \right)dt\prime =u\left( t \right),0\le t\le 2\pi }.$

(39)

用配置法离散式(39)得到离散的算子方程H(t)x(t)=u(t).

假设观测数据是在真实数据上添加一定噪声水平的高斯白噪声^[19]，则有

${{u}_{e}}=u+e*rand(size\left( u \right)),$

(40)

其中，e表示噪声水平，rand(·)表示维数与u相同的高斯白噪声，即有

$H\left( t \right)x\left( t \right)={{u}_{e}}\left( t \right).$

(41)

问题(41)等价于求解下列无约束优化问题:

$min~f\left( x\left( t \right) \right)=\|H\left( t \right)x\left( t \right)-{{u}_{e}}\left( t \right){{\|}^{2}}.$

(42)

为了克服问题(41)的不适定性，引入Tikhonov正则项，即该问题转化为求解下述问题:

$\begin{align} & min~f\left( x\left( t \right) \right)=\|H\left( t \right)x\left( t \right)-{{u}_{e}}\left( t \right){{\|}^{2}} \\ & +\gamma \|x\left( t \right){{\|}^{2}}, \\ \end{align}$

(43)

这里，α表示正则化参数.

选择核函数H(t)和真实解x_true(x):

${{x}_{true}}=exp~(3sin~t),$

(44)

$H\left( t \right)=\frac{1}{0.05\sqrt{\pi }}exp~(-\frac{{{\left( t-1/200 \right)}^{2}}}{{{0.05}^{2}}}).$

(45)

在实验中，离散点设置为256个单位，即矩阵H(t)∈R^256×256，真实解x_ture(w,t)∈R²⁵⁶，观测数据u_e(t)的噪声水平设置为e=0.01，现用观测数据u_e(t)来反演x(t).

为了更公平的测试算法的速度，每个实验重复做10次，结果的平均值见表 1.从表 1可以看出本文提出的谱共轭梯度法(SCG)的REMS值相对较小,且求解过程所需的迭代次数和CPU时间也是最小的，说明SCG方法所需更少的计算代价.

从图 1~4可以看出SCG方法具有一定的优势.从图 5的效果图也可以得出SCG方法是有效的.

例2 图像处理的反演问题

考察简单的一维情形图像识别问题,如方程：

$\int_{a}^{b}{h\left( t-t\prime \right)x\left( t\prime \right)=u\left( t \right)},$

(46)

其中x(t′)代表光源强度，u(t)代表观测的图像，通常是一个模糊核的连续图像；h(t－t′)表示成像的过程，通常是一个点扩展函数(PSF).问题给定光源函数x(t′)和点扩展函数h(t－t′)来求观测图像u(t).与此相关的反问题：给定点扩展函数h(t－t′)和观测图像u(t)来确定光源强度x(t′).

由有限差分法离散式(39)得到离散的算子方程

$H\left( t \right)x\left( t \right)=u\left( t \right),$

(47)

假设观测数据为在真实数据上添加imnoise函数产生的高斯噪声，则有

${{u}_{e}}=u+e*imnoise(size\left( u \right)),$

(48)

其中，e表示噪声水平，imnoise(·)表示维数与u相同的高斯噪声，即有

$H\left( t \right)x\left( t \right)={{u}_{e}}\left( t \right).$

(49)

问题(49)等价于求解下列无约束优化问题:

$min~f\left( x\left( t \right) \right)=\|H\left( t \right)x\left( t \right)-{{u}_{e}}\left( t \right){{\|}^{2}}.$

(50)

为了克服问题(49)的不适定性，引入Tikhonov正则项，即转该问题化为求解下述问题:

$\begin{align} & min~f\left( x\left( t \right) \right)=\|H\left( t \right)x\left( t \right)-{{u}_{e}}\left( t \right){{\|}^{2}} \\ & +\gamma \|x\left( t \right){{\|}^{2}}, \\ \end{align}$

(51)

这里，γ表示正则化参数.

选择核函数H(t)和真实解x_true(x):

${{x}_{true}}=si{{n}^{4}}(2\pi t),$

(52)

$H\left( t \right)=\frac{1}{1+100{{\left( t-0.5 \right)}^{2}}~}.$

(53)

在实验中，离散点设置为512个单位，即矩阵H(t)∈R^512×512，真实解x_ture(w,t)∈R⁵¹²，观测数据u_e(t)的噪声水平设置为e=0.02，现用观测数据u_e(t)来反演x(t).

为了更公平的测试算法的速度，每个实验重复做10次(表 2).从表 2可以看出,SCG方法的RMSE值是0.012，在4个算法中相对较小，而且求解过程所需的迭代次数和CPU时间都相对较小，因此SCG方法所需求解成本相对较低.

从图 6~9可以看出SCG方法的数值解比较理想的逼近真实解.从图 10的效果图也可以得出,SCG方法比其他3种算法更有效.

4 结论

本文提出一个求解线性逆问题的新谱共轭梯度法，证明了该算法的全局收敛性,并利用该方法对卷积方程的反演问题和图像处理的反演问题进行数值实验，发现该方法比Landweber，TSVD和TV方法更有效.