1. 引言
近年来数学工具不断被引用到诸多领域。例如,通过地震波的测量来判断地下矿藏的位置或地球内部的结构;用X光分层扫描构像作医学诊断等等,都是在不能达到或直接接触到研究对象的情况下,运用特殊的手段来获得相关解的某些信息,进而转化为数学上的反问题来求解。而反问题最显著的特性就是不适定性和非线性性 [1] 。
考虑非线性不适定算子方程
(1)
其中,
,
为Hilbert空间。在实际问题中,不能得到y的精确值,只能得到误差水平为
的扰动数据
,即
(2)
针对方程(1)的求解问题,需要引入正则化方法获得满足要求的近似解。其中,正则化方法包括Tikhonov正则化和迭代正则化等。
1951年,Landweber提出了解决线性问题的Landweber迭代正则化方法 [2]
(3)
其中,k表示迭代步数,
为初始值,
为Fréchet导数,*表示算子的伴随。同年,Hanke等人 [3] 证明了Landweber迭代正则化方法是求解非线性不适定问题的一种稳定方法。
Landweber迭代正则化方法因其具有良好的稳定性特点,常常用于求解大规模不适定问题。然而,由于Landweber迭代正则化方法的迭代序列往往收敛速度较慢,限制了Landweber迭代正则化方法在实际问题中的广泛应用。因此,在求解不适定问题中,提高Landweber迭代正则化方法的收敛速度是非常有必要的。近年来国内外学者对Landweber迭代正则化方法的加速做了大量工作,2018年,Simon Hubmer和Ronny Ramlau [4] 利用Nesterov加速梯度法对Landweber迭代正则化方法进行加速,理论上给出了该方法在不适定问题上的收敛性分析,基于非线性对角算子和自卷积的反问题从数值上验证了该方法的有效性。2021年,Guangyu Gao等人 [5] 结合投影两点梯度法和Kaczmarz方法,提出了一种加速的Kaczmarz迭代方法用于求解多种类型的非线性不适定问题。2021年,Yuxin Xia等人 [6] 提出了基于Nesterov策略的加速homotopy扰动迭代法,基于椭圆参数识别问题验证了该方法的高效性和有效性。本文借鉴Nesterov加速策略思想,基于Landweber迭代正则化方法,引入自适应BB步长加速技巧,提出了加速的Landweber迭代正则化方法。
2. 基于自适应BB步长的修正Landweber迭代算法及收敛性分析
经典的Landweber迭代正则化方法可以看成极小化二次泛函
的步长为1的最速下降法。在优化理论中,步长的选取对于方法的收敛速度有着非常重要的影响。步长的选取一直是优化理论研究的热点和难点问题,很多学者对基于梯度信息构造步长做了大量工作。1988年,Borwein与Barzilai [7] 提出了Barzilai-Borwein (BB)梯度法,这种方法的主要思想是利用前一步迭代的信息来构建当前迭代步的步长。BB步长受到了许多学者的关注,袁亚湘院士等人 [8] [9] 将最速下降步长与BB步长相结合,给出了保持单调并且可与BB方法计算效果相媲美的梯度算法,并且对梯度法中步长因子的选取提出了许多新的改进 [10] [11] [12] [13] 。2006年,Bin Zhou等人 [13] 基于求解二维二次方程BB方法的超线性性质,提出了自适应BB方法,即在每次迭代中自适应的选择大步长和小步长,且通过数值实验表明该方法优于BB方法。2020年,王彦飞等人 [14] 在利用X射线计算机断层成像技术进行图像恢复问题中,引入自适应BB步长加速技巧提出快速梯度下降法求解X射线CT成像的大规模问题。
BB步长格式如下:
(4)
其中,
,
。
自适应BB(ABB)步长为:
(5)
其中,
且逼近于
,
。当ABB步长使用较小步长时可以诱导一个良好的下降步长,而较大的步长有利于产生足够的下降量 [8] 。
借鉴ABB步长的良好性质,将ABB步长引入到Landweber迭代正则化方法中,得到加速的Landweber迭代正则化方法
(6)
为了防止ABB步长
在收敛点附近过大,在算法中对ABB步长进行截断处理,通过引入截断参数
得到Landweber-ABB (LABB)迭代正则化方法
(7)
其中,
为 [15]
算法1所示为LABB迭代正则化方法的算法描述。
3. 收敛性分析
在非线性不适定问题(1)中,假设算子F在
的邻域
内Fréchet导数存在,一般情况下,仅局限于算子F的Fréchet导数的一致有界是不够的,同时要求算子F满足 [16]
(8)
由假设(8)易得
(9)
迭代中,从解的精度着眼,需要迭代次数充分大,但从稳定性考虑又要求迭代次数不能太大,因此迭代终止判据的给出是必然的 [16] 。根据广义偏差原理
(10)
其中,
为依赖于
的正数,满足
(11)
LABB迭代正则化方法的迭代值
满足上述迭代停止准则时,则停止迭代,迭代值
即为所求方程满足条件的近似值。
基于非线性条件(8)可得LABB迭代正则化方法的迭代误差的单调性。为了方便读者,在此给出简要证明,更多细节参考文献 [15] 。
定理3.1假设非线性条件(8)成立,
是非线性不适定问题(1)在邻域
内的解,方程右端项扰动数据
满足扰动误差
,
为满足(10)和(11)的最小迭代步,则有
(12)
此外,对
有
(13)
其中
(14)
证明:迭代停止条件
,根据(2),(8)及LABB迭代正则化方法格式,运用数学归纳法可以得到迭代值
,且
由偏差原理(10)和(14)可得
(15)
又有
因此(15)化简为
由(10)可知,
,有
因此可得迭代误差是递减的,即证。
为了证明收敛性,
必须满足连续且
(16)
下面将给出LABB迭代正则化方法的收敛性证明。
定理3.2假设条件(8)和(16)成立,非线性不适定问题(1)在邻域
内有解。令
(误差水平
),即迭代步
收敛于问题(1)的真实解。
表示问题(1)在邻域
内的最小二乘解,若对
,都有
成立,那么
。
证明:本定理证明见参考文献 [15] 。
在上诉定理中,仅对误差水平
时,LABB迭代正则化方法的收敛性进行了讨论。当存在误差水平时,如果LABB迭代正则化方法按照偏差原理停止迭代,那么可以得到如下定理,即能说明LABB迭代正则化方法是求解非线性不适定问题的正则化方法 [17] 。
定理3.3假设条件(8)和(16)成立,非线性不适定问题(1)在邻域
内有解。
为满足(10)和(11)迭代停止的最小迭代步,
表示问题(1)在邻域
内的最小二乘解,则迭代步
收敛于问题(1)的真实解。若对
,都有
成立,那么
。
证明:本定理证明见参考文献 [15] 。
4. 数值实验
为了从数值计算角度验证本文所提出的方法的有效性,考虑如下基于Dirichet边界条件的椭圆方程参数c的识别问题,且参数c满足
(17)
其中,
,
,
,
为有界域,u在边界
处Lipschitz有界,c为重构参数。在偏微分方程理论中,当参数c满足如下条件,(17)存在唯一解:
其中,
。定义非线性算子
为
于是椭圆参数识别问题转化为如下非线性方程的求解问题
(18)
参数识别问题(18)是一个非线性不适定问题,因此,在求解过程中需要引入正则化方法。
为了提高方法对解的不连续信息或边界的重构分辨率,在数值实验中加入TV正则化,即在第k步时极小化问题
其中,
表示重构问题的计算区域,
表示模型噪声方差,在本文实验中取0.001,
表示第k步LABB迭代正则化方法的迭代值。假设(17)的解为
,于是LABB迭代正则化方法的下一次迭代初值为
。在问题(17)的求解过程中,本文运用了主对偶混合梯度法(PDHG) [18] 。
为了定量分析计算效率和精度上的重构性能,从算法的计算时间(Time)、迭代次数(N)、模型相对误差
(
)以及数据残差(
)等角度对算法进行综合对比分析。由于经典Landweber
迭代正则化方法收敛速度过慢,因此,为了说明LABB迭代正则化方法的有效性,重点介绍该方法与修正的Landweber迭代正则化方法比较,修正的步长为 [19] :
其中,本文实验选取
,
,
,
。
4.1. 一维问题
在模拟中,考虑边界条件
,
和区间
的一维问题,并假定真参数为
其中,
,在u上添加随机高斯噪声,噪声数据
满足
。选取初始值,且偏差原理(10)中对所有的方法均选
。为了进行计算,将区间
分成128个等长的子区间,其中所有偏微分方程均用有限差分法近似求解 [6] 。为了防止因为误差所引起的ABB步长过大问题,当数据误差发生震荡时用下山法获取步长。根据ABB步长的特点,在实验中将另设一个停止条件,即满足条件
或偏差原理时停止迭代。
利用Matlab进行数值试验,取不同误差水平进行独立的实验,实验中总的计算(Time),最大迭代步数(N),模型残差(RES)和模型相对误差(ERR)如表1所示。从表1可以看出,相对于修正的Landweber迭代正则化方法,LABB迭代正则化方法在计算时间,最大迭代步数及收敛速度方面具有较强的竞争力。
Table 1. Numerical results of one-dimensional elliptic parameter identification
表1. 一维椭圆参数识别数值结果
图1展示了噪声水平为0.001的LABB迭代正则化方法的迭代结果,图1(a)为修正Landweber迭代正则化方法重构结果,图1(b)为LABB迭代正则化方法重构结果,表明LABB迭代正则化方法的迭代值可以很好的近似真实解;如图1(c)所示,随着迭代的进行,LABB迭代正则化方法较修正Landweber迭代正则化方法收敛快;为了看清前几步的收敛曲线,故将前50步的残差另作图(图1(d)),如图1(d)所示,LABB迭代正则化方法比修正Landweber迭代正则化方法收敛快。
4.2. 二维问题
在模拟中,区间为
的二维问题,真实值设定为
源项
,噪声数据
满足
。选取初始值
。为了便于比较,对于所有方法,偏差原理(10)中的参数
均设为1.1。将
分成64 × 64个等长小方块,用有限差分法近似求解所有涉及的偏微分方程 [6] 。同理,实验中为了防止因为误差所引起的ABB步长过大问题,当数据误差发生震荡时进行步长线搜索。根据ABB步长的特点,在实验中另设一个停止条件,即满足条件
或偏差原理时停止迭代。
利用Matlab进行数值试验,取不同误差水平进行独立的实验,实验中总的计算(Time),最大迭代步数(N),模型残差(RES)和模型相对误差(ERR)如表2所示。从表2可以看出,相对于修正的Landweber迭代正则化方法,LABB迭代正则化方法在计算时间,最大迭代步数及收敛速度方面具有较强的竞争力。
Table 2. Numerical results of two-dimensional elliptic parameter identification
表2. 二维椭圆参数识别数值结果
图2展示了噪声水平为0.001的LABB迭代正则化方法的迭代结果,图2(c)为LABB迭代正则化方法重构结果,表明LABB迭代正则化方法的迭代值可以很好的近似真实解;如图2(d)所示,随着迭代的进行,LABB迭代正则化方法较修正Landweber迭代正则化方法收敛快。
5. 总结
经典的Landweber迭代正则化方法具有计算稳定、算法实现简单等特点。但由于其收敛速度过慢影响了算法的广泛应用。本文基于Landweber迭代正则化方法,引入ABB方法,提出LABB迭代正则化方法。理论上证明了该方法的收敛性,并基于非线性不适定椭圆参数识别方程,对修正的Landweber迭代正则化方法以及LABB迭代正则化方法进行对比数值实验,结果表明本文提出的LABB迭代正则化方法是收敛的,和理论验证一致;在达到相同收敛精度下,LABB迭代正则化方法运算时间较少,表明了LABB迭代正则化方法在计算效率上的优势。