同伦摄动稀疏正则化方法及其在.doc

同伦摄动稀疏正则化方法及其在
初始地形地貌重构中的应用
摘要：稀疏正则化方法在参数重构中起到了越来越重要的作用。与传统的正则化方法相比，稀疏正则化方法能较好地重构稀疏变量。由于稀疏正则化的不可微性，需要对已有的经典算法进行改进。本文构建同伦摄动稀疏正则化方法克服标准稀疏正则化的不可微性，并应用该方法应用到项目问题中，能够有效地重构山体初始表面。数值实验表明，所提出的方法是收敛和稳定的。
关键词：稀疏正则化；同伦摄动稀疏正则化方法；参数重构；重构山体初始表面

1. 引言
参数重构在许多应用中起到重要作用，如波动率和政策参数重构[1, 2]，在其他工程实践领域也有重要的应用，例如：图像去噪[3, 4]，地震信号重构[5, 6]以及心电信号重构 [7, 8]。一般情况下，参数重构问题是不适定的，也就是说，即使测量数据的噪声水平很小，也可能导致重构结果严重偏离真实参数[9]。正则化方法主要就是为了克服参数重构的不适定性，通过选取合适的正则化方法能够抑制观测数据中误差对于参数重构的不良影响，获得较为准确的重构值。
使用最为普遍的正则化方法是吉洪诺夫正规化方法，它的目标函数是由拟合项和惩罚项构成。吉洪诺夫正则化方法已被用于许多参数辨识问题中，很多学者对其数值计算方法进行了研究，如Landweber方法[10]，高斯牛顿法[11]，Newton-Kaczmarz方法[12]和多尺度平滑方法[13]，这些方法能够有效地重构光滑参数。
随着经济和金融理论的发展，波动率和经济政策参数的重构已经广泛应用于许多实际问题中。在实际应用中，很多需要重构的参数都是稀疏的，即参数的非零元素的数目非常有限，远远小于零元素的数目。即使重构参数稀疏程度不够，也可以利用小波和曲波变换使参数稀疏化。本文先对该方法在经济和金融领域中的应用进行研究，因为这两个领域中的参数通常可以分成已知部分和未知部分。已知部分通常与已有的经济和金融政策相对应，而未知部分通过参数重构，再将结果进行经济学和金融学解释，能够为政策制定部门提供行之有效的对策建议。通过这两类问题的研究表明该方法能够应用到参数分解成已知部分和未知部分的问题中。然后再将该方法应用到较为复杂的山体表面重构问题中。
传统的吉洪诺夫正则化方法对于重构稀疏参数效果很不好，而稀疏正则化方法却能较好地重构稀疏参数，但是稀疏正则化是不可微的，因此需要采用一些技巧来克服这一困难，典型的方法是Bregman迭代[14, 15-17]。本文构造同伦摄动稀疏正则化方法，达到提高算法精度和提高计算效率的目的。数值实验表明同伦摄动稀疏正则化方法是收敛和稳定的。
2. 稀疏正则化方法
在实际应用中，对于光滑参数重构，吉洪诺夫正则化方法具有良好的收敛性和稳定性。但是，在稀疏参数重构时，吉洪诺夫正则化方法的重构效果很差，很难满足工程实践的要求。因此，需要采用稀疏正则化方法进行参数重构。
2.1 吉洪诺夫正则化方法
参数重构问题可以归结为非线性算子方程(1)的形式：

(1)
其中， 分别代表非线性算子、需要重构的参数和观测数据。在实际问题中，往往还需要考虑观测数据和理想数据之间的噪音水平，即：

(2)
其中，分别代表真实的观测数据和噪音水平。
参数重构的难点在于不适定性，很小的噪音水平 也会使得重构结果严重背离真实的物理参数，从而造成结果的无意义。解决这一难点最重要的方法就是吉洪诺夫正则化方法，与之对应的吉洪诺夫正则化目标函数定义为：

(3)
其中， 是在范数意义下的数据拟合项, 是起到稳定作用的惩罚项, 参数为正则化参数, 该参数主要是起到平衡数据拟合项和惩罚项的作用。
吉洪诺夫正则化方法主要是求解下面的优化问题
.
(4)
最小值满足如下的欧拉方程：

(5)
其中， 是F-导数。在求解方程(5)的时候，采用的最普遍的方法是Landweber方法，该方法可以写成下面的表达式：
,
(6)
其中，表示迭代次数。方程(6)是著名的Landweber迭代格式，该数值格式的显著优点是稳定性特别好，但是，收敛速度很慢，不适合应用到大型实际问题中。另一个重要的方法就是高斯-牛顿方法，该方法收敛速度较快，但是，不如Landweber方法稳定。本文主要是在Landweber方法基础上研究同伦摄动稀疏正则化方法，整个过程可以推广到高斯-牛顿方法上。
2.2 稀疏正则化方法
为了能够有效地重构稀疏变量，将标准的吉洪诺夫正则化方法进行改进，使得吉洪诺夫正则化目标函数转换为如下形式:
,
(7)
其中，表示向量的非零元素的个数。近些年，由于稀疏正则化方法能够有效地重构稀疏变量，使得反问题领域学者越来越多地重视该方法在实际问题中的应用。该方法的难点在于，泛函(7)的惩罚项是不可微的，而且该问题还是NP不可解问题。为了解决NP不可解难点，采用下面的形式进行替代：
,
(8)
其中，范数表示. 在泛函(8)中, 利用范数代替了原有的范数，这样的改进使得计算效率得到了显著的提高。
尽管进行了上述的改进，但是(8)的惩罚项依然是不可微的，因此，还要在此基础上进一步改进，将(8)中的惩罚项改为下面的带有阻尼系数的范数：
,
(9)
其中，利用代替范数, 辅助参数是一个正实数。与吉洪诺夫正则化方法相仿，将Landweber方法应用到(9)中，得到下面的迭代格式：
.
(10)
3. 同伦摄动稀疏正则化方法
为了提高数值迭代格式(10)的收敛速度，将同伦摄动反演算法与稀疏正则化方法有机地融合到一起，达到提高重构效果和收敛速度的目的。
构造如下的同伦映射
,
,
(11)
其中，代表嵌入参数, 迭代算法的初始值。因此，有下面的关系式：

(12)
将写成关于嵌入参数的幂级数形式：

(13)
目标函数(9)的解可以表示为：
.
(14)
将方程(11)中的在点写成泰勒级数的形式：

(15)
可以得到下面的关系式：
.
(16)
和
.
(17)
根据(17)，对于含有噪音的观测数据，采用前两项构造稀疏变量:
,
(18)
其中，表示迭代次数。
采用第一项构造稀疏变量，则有下面的迭代格式：
,
(19)
迭代格式(19)正是前面介绍的Landweber方法，而迭代格式(18)可以看做是在其上面进行的修正，将更多的信息应用到迭代算法中，可以提高计算效率。
4 数值算例
4.1 重构隐含波动率
在金融领域中，基于布莱克一斯科尔斯期权定价模型(Black-Scholes (B-S) option pricing model)重构隐含波动率一直是重要的研究方向。利用有限差分方法求解该模型的正演问题，在重构隐含波动率的时候采用同伦摄动稀疏正则化方法，达到提高重构精度和计算效率的目的。
布莱克一斯科尔斯期权定价模型表明金融衍生品的价格变化满足标准几何布朗运动。该模型的边界条件随着金融衍生品的变化而变化。当边界条件给定的时候，金融衍生品的价格可以通过求解布莱克一斯科尔斯期权定价模型获得。当金融衍生品选定为期权时，本文定义正演问题为求解期权的价格，反演则是重构隐含波动率。根据期权的行权方式的不同，可以分为欧式期权和美式期权，本节以欧式期权为研究对象。
在时间上的欧式看涨期权, 设 是欧式看涨期权在区域上的价格, 满足下面的方程：

,
(20)

其中，: 股票价格, : 执行价格, : 利率, : 红利, : 到期日, : 时间, : 隐含波动率。.
反问题可以表述为: 从观测数据, 重构隐含波动率.
定义非线性向量值函数, 即：. 为了测试本文提出的方法和标准的吉洪诺夫正则化方法之间的差别，设置如下参数：到期日, 股票价格，利率, 执行价格. 为了测试算法的抑制噪音的能力，在观测数据上添加1%的高斯噪音。真实的隐含波动率如下：

(21)
将真实的波动率分解为两部分：
,
(22)
其中，分别表示常值隐含波动率和稀疏隐含波动率。设是已知的, 需要重构部分，该部分为：

(23)
目标函数(9) 可以写成：
,
(24)
采用同伦摄动稀疏正则化方法重构结果为：

(25)
采用标准吉洪诺夫正则化方法重构结果为：

(26)
从上面的重构结果可以看出，同伦摄动稀疏正则化方法能够提高标准吉洪诺夫正则化方法对于重构稀疏隐含波动率的精度。
4.2 托达罗模型的政策参数重构
托达罗模型是著名的描述人口流动的经济学模型，在劳动经济学领域具有重要的地位。托达罗模型将人口流动的数量和不同地区的收入差距联系到一起：

(27)
其中，分别表示人口流动的数量和不同地区之间的收入差距。函数是单调递增的，即. 随着经济的发展，标准的托达罗模型(27)已经不能准确地描述人口流动和收入差距之间的关系，还需要考虑政策参数的作用，将方程(27)改为下面的形式：

(28)
其中，表示政府部门的政策参数。
设有个农村地区和个城市地区，第个农村区域和第个城市区域之间的收入差距表示为(), 第个和 第个城市区域之间收入差距表示为(). 人口转移到第个城市区域的数量用表示。 对于第城市区域, 政策参数可以分解为两部分, 其中分别表示政策参数对于农村和城市的作用效果。为了简化讨论，假设从城市区域转移到农村区域的人口数量为零。
在实际应用中，函数有多种表达方式，由于本文主要是测试同伦摄动稀疏正则化方法的有效性，因此，函数选取为线性函数的形式。托达罗模型可以写成下面的形式：

(29)
其中，、和 已知，政策参数向量是未知的。
令
,
(30)
和
,
.
(31)
方程 (30) 可以写成下面的形式：
.
(32)
值得注意的是方程(32)是欠定的，因此，需要正则化方法求解该方程。将政策参数进行分解 , 其中分别表示以前的政策参数(已知)和稀疏政策参数(未知)。目标函数改写为：
.
(33)
在数值算例中， 设, 的每个元素都是. 为了测试算法的抑制噪音的能力，在观测数据上添加1%的高斯噪音。重构5个稀疏政策参数，每一个 有且只有一个非零元素，大小为0.5. 同伦摄动稀疏正则化方法能够精确地将非零元素的位置判断出来，但是标准的吉洪诺夫正则化方法就不具备这样的能力。表1给出了两种不同正则化方法重构的均方差。
表1：两种不同方法的均方差
Table 1: Mean Square Errors
非零元素的位置
同伦摄动稀疏正则化方法
吉洪诺夫正则化方法
2
5.25%
10.07%
4
5.19%
12.96%
6
4.69%
15.78%
8
4.92%
14.06%
10
4.71%
9.93%
4.3重构初始地形地貌
地形地貌演变模型将引起地球表面变化的内因和外因有机的结合到一起，
影响地形地貌演变的内因是地球内部板块速度场, 而外因则非常多. 主要有以
下三个因素：(I) 地球表面扩散过程(diffusion of hillslope topography): 这里所说的扩散过程是一个包含很多复杂因素的统称，例如：风化作用(weathering), 坡面冲刷(slope wash), 地面水流(overland flow), 土壤滑动(soil creep)和基岩滑坡引起的质量损失(mass wasting by bedrock-involved landsliding); (II) 河流冲刷对于基床的切割(bedrock incision)；(III) 冰川的融化影响(melting effect).建立地形地貌演变数学模型是属于地形地貌学的范畴，描述上面因素的数学公式一般都是高度非线性的，将这些具体模型的模型进行统一考虑，可以得到如下模型：
(34)
模型正演定义为： 给 定 地 球 表 面 的 速 度 场和 初 始 地 形 地 貌 分 布 函数求解地形地貌随时间变化关系，写成数学表达式为：
(36)
模型反演的主要目的是利用现在的地形地貌，重构初始地形地貌，可以表示为：
(37)
值得注意的是我们这里假设山体的随时间运动的速度场是已知的，将研究的重点放到了地形地貌运动模型上，并没有放到地球内部热场分布模型上，该部分的研究已经在上一年度完场。
方程 (37) 可以写成下面的线性形式：
.
(38)
值得注意的是方程(38)是不适定的，因此，需要正则化方法求解该方程。将初始地形地貌进行分解 , 其中分别表示通过其他探测手段已知的地形地貌成分(尺度较大)和细节地形地貌成分(尺度较小)。目标函数改写为：
.
(39)
在数值算例中，，在小波基下表示为稀疏系数，重构这些稀疏系数。 为了测试算法的抑制噪音的能力，在观测数据上添加1%的高斯噪音。同伦摄动稀疏正则化方法能够精确地将非零元素的位置判断出来，但是标准的吉洪诺夫正则化方法就不具备这样的能力。表2给出了两种不同正则化方法重构的均方差。从表2中还可以看出新方法更加稳定。
表2：两种不同方法的均方差
Table2: Mean Square Errors
实验
次数
同伦摄动稀疏正则化方法
吉洪诺夫正则化方法
1
2.21%
8.07%
2
2.16%
8.16%
3
2.23%
7.98%
4
2.39%
9.01%
5
2.31%
8.58%

5 结论
本文构造了同伦摄动稀疏正则化方法，提高了标准吉洪诺夫正则化方法重构稀疏变量的精度，同时也提高了计算效率，并将该方法应用到重构初始地形地貌问题中。对于实际情况而言，绝大部分参数变量都可以分解为已知部分和稀疏部分，因此只需要重构对应的稀疏部分即可，说明该方法具有较好的应用前景。同时，数值算例表明本文构造的算法具有较高的精度和收敛速度。
参考文献
[1] Egger H, Engl H W. Tikhonov regularization applied to the inverse problem of option pricing: convergence analysis and rate[J]. Inverse Problems, 2005, 21: 1027–1045.
[2] Trong D D, Thanh D, Lan N N, Uyen P H. Calibration of the purely T-dependent Black-Scholes implied volatility[J]. Applicable Analysis, 2014, 93: 859-874.
[3] Wu C D. A new denoising method based on multistage median filter and nonsubsampled contourlet transform[J]. Journal of Computational Information Systems, 2012, 8(12): 4881-4888.
[4] Zhao Q, Ye B, Cao J. Image denoising based on improved non-local means and nonsubsampled contourlet transform wiener filtering[J]. Journal of Computational Information Systems, 2010, 6 (2): 601-609.
[5] Baig A M, M. Campillo M, Brenguier F. Denoising seismic noise cross correlations[J]. Journal of Geophysical Research, 2009, 114: doi: 10.1029/2008JB006085.
[6] Hennenfent G, Herrmann F J. Simply denoise: wavefield reconstruction via jittered undersampling[J]. Geophysics, 2008, 73 (3): V19-V28.
[7] Istepanian R H, Petrocian A A. Optimal zonal wavelet-based ECG data compression for a mobile telecardiology system[J]. IEEE Transactions on Information Technology in Biomedicine, 2000, 4 (3): 200-211.
[8] Sameni R, Shamsollahi M B, Jutten C, Clifford G.D. A nonlinear Bayesian filtering framework for ECG denoising[J]. IEEE Transactions on Biomedical Engineering, 2007, 54 (12): 2172-2185.
[9] Engl H W, Hanke M, Neubauer A. Regularization of Inverse Problems[M]. Kluwer: Dordrecht, 1996.
[10] Hanke M, Neubauer A, Scherzer O. A convergence analysis of the Landweber iteration for nonlinear ill-posed problems[J]. Numer. Math., 1995, 72: 21-37.
[11] Blaschke B, Neubauer A, Scherzer O. On convergence rates for the iteratively regularized Gauss-Newton method[J]. IMA J. Numer. Anal., 1997, 17: 421-436.
[12] Burger M, Kaltenbacher B. Regularizing Newton-Kaczmarz methods for nonlinear ill-posed problems[J]. SIAM J. Numer. Anal., 2006, 44: 153–182.
[13] Bredies K, Lorenz D, Maass P. Mathematical concepts of multiscale smoothing[J]. Appl. Comput. Harmon. Anal, 2005, 19: 141–161.
[14] Brune C, Sawatzky A, Burger M. Primal and dual Bregman methods with application to optical nanoscopy[J]. Int. J. Comput. Vis, 2011, 92:211-229.
[15] Burger M, Resmerita E, He L. Error estimation for Bregman iterations and inverse scale space methods in image restoration[J]. Computing, 2007, 81(2-3):109–135.
[16] Cai J F, Osher S, Shen Z. Convergence of the linearized Bregman iteration for l1-norm minimization[J]. Math. Comput. , 2009, 78: 2127–2136.
[17] Cai J F, Osher S, Shen Z. Linearized Bregman iterations for compressed sensing[J]. Math. Comput., 2009, 78: 1515–1536.