Journal of Image and Signal Processing
Vol. 11  No. 03 ( 2022 ), Article ID: 53890 , 9 pages
10.12677/JISP.2022.113014

基于共轭梯度的脉冲噪声去噪算法

孙敏1,田茂英2

1枣庄学院,数学与统计学院,山东 枣庄

2山东煤炭卫生学校,山东 枣庄

收稿日期:2022年6月30日;录用日期:2022年7月11日;发布日期:2022年7月21日

摘要

共轭梯度法是数值计算领域的一类重要方法。本文基于Schmidt正交化方法,提出了两种改进的Hestenes-Stiefel共轭梯度法。无需线搜索条件,它们自动满足共轭梯度法的两个重要性质,即充分下降性和共轭条件。在适当的条件下,本文证明了所设计算法的全局收敛性。最后,将这两种算法应用到脉冲噪声去噪中,初步的数值试验表明了算法的有效性。

关键词

共轭梯度,收敛性,脉冲噪声去噪

Impulse Noise Denoising Algorithm Based on Conjugate Gradient

Min Sun1, Maoying Tian2

1School of Mathematics and Statistics, Zaozhuang University, Zaozhuang Shandong

2Shandong Coal Medical School, Zaozhuang Shandong

Received: Jun. 30th, 2022; accepted: Jul. 11th, 2022; published: Jul. 21st, 2022

ABSTRACT

Conjugate gradient method is an important method in numerical computation. Based on Schmidt orthogonalization, two improved Hestenes-Stiefel conjugate gradient methods are proposed in this paper. Without line search conditions, they automatically satisfy two important properties of the conjugate gradient method, namely, sufficient descent and conjugate conditions. Under appropriate conditions, the global convergence of the proposed algorithm is proved. Finally, the two algorithms are applied to impulse noise denoising, and preliminary numerical experiments show the effectiveness of the algorithm.

Keywords:Conjugate Gradient, Convergence, Impulse Noise Denoising

Copyright © 2022 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

考虑无约束最优化问题(记为UCP):

min x R n f ( x ) (1)

其中 f : R n R 是一个光滑函数,且有下界。UCP是最优化领域中最重要和最基本的问题之一,学者们对其进行了系统而深入的研究。有许多经典的迭代方法来求解UCP,包括最速下降法、共轭梯度法、信赖域法、牛顿法、拟牛顿法等 [1] [2]。尽管信赖域方法、牛顿方法和拟牛顿方法具有很好的收敛性,但它们不适合大规模的UCP,因为它们需要计算Hessian矩阵 2 f ( x k ) 或其近似值,从而导致大量的计算和存储。最速下降法和共轭梯度法在迭代方案中只包括梯度 f ( x k ) ,因此,它们适用于求解大规模的UCP。

共轭梯度法的基本思想可以追溯到Hestenes和Stiefel求解线性方程组的研究成果。后来,Fletcher和Reeves将其推广到解决非线性优化问题。共轭梯度法自出现以来,一直受到学者们的密切关注。从初始迭代 x 0 开始,共轭梯度法的迭代格式为

x k + 1 = x k + α k d k , k = 0 , 1 , (2)

其中 α k 是由某些线搜索计算的步长, d k 是生成的搜索方向,其定义下式

d k = { g k if k = 0 , g k + β k d k 1 , if k 1 , (3)

其中, g k = f ( x k ) β k 是共轭参数,这是不同共轭梯度方法之间的根本区别。比较著名的共轭参数包括:Hestenes-Stiefiel (HS),Polak-Ribiere-Polyak (PRP),Fletcher-Reeves (FR),与Dai-Yuan (DY)。具体定义如下 [3]:

β k HS = g k ( g k g k 1 ) d k 1 ( g k g k 1 ) , β k PRP = g k ( g k g k 1 ) g k 1 2 ,

β k FR = g k 2 g k 1 2 , β k DY = g k 2 d k 1 ( g k g k 1 ) .

如果UCP中的 f ( x ) 是一个强二次函数,并且(2)中的步长 α k 是通过精确的线搜索生成的,则上述四个 β k 公式是等价的,否则可以得到四个著名的共轭梯度方法,即HS方法、PRP方法、FR方法和DY方法。这四个方法具有不同的收敛结果和数值性能。FR方法和DY方法分别在强Wolfe线搜索和Wolfe线搜索下具有全局收敛性。对于HS方法和PRP方法,它们通常被认为是实际应用中最有效的两种共轭梯度法,但即使步长 α k 是精确步长,它们也可能发散。

大量数值试验表明,高效的共轭梯度法产生的搜索方向满足两个性质:充分下降条件和共轭条件,也就是

g k d k c g k 2 , c > 0 (4)

y k 1 d k = 0 (5)

其中 y k 1 = g k g k 1 。近几年学者们设计了多种满足(4)或(5)的共轭梯度法。比如,Hager与Zhang [4] 将(3)中的参数 β k 取为

β k HZ = 1 d k 1 , y k 1 y k 1 2 d k 1 y k 1 2 d k 1 , y k 1 , g k .

如果 d k 1 , y k 1 0 ,则有

g k d k 7 8 g k 2 .

文献 [5] [6] 提出了两种修正的PRP共轭梯度法,其方向取为

d k ZTPRP = g k + β k PRP d k 1 θ k y k 1 ,

d k CTPRP = ( 1 + β k PRP g k d k 1 g k 2 ) g k + β k PRP d k 1 ,

其中 θ k = g k d k 1 g k 1 2 。容易验证上面两类方向都满足 g k d k = g k 2 ,并且容易看出 d k CTPRP 的设计动机来自于Schmidt正交化方法。基于这个设计思想, [7] 提出了一类满足夹角性质的记忆梯度法,其产生的方向同时满足(4)与夹角条件 cos g k , d k 1 / 3

本文继续对这一方向进行研究,提出了两种同时满足充分下降性和共轭条件的改进HS共轭梯度法。本文的其余部分组织如下:第2节提出了HS方法的两种变体,并证明了生成序列的一些重要性质。第3节证明了它们的全局收敛性。在第4节中,我们给出了一些图像去噪的试验结果。第5节给出一些未来的研究方向。

2. 两类改进的HS方法

在这一节中,我们设计两种改进的HS方法的迭代格式。在无需线搜索的条件下,它们产生的方向满足充分下降条件和共轭条件。

正如在第1节中所分析,由下式产生的方向 d k 满足充分下降条件:

d k = ( 1 + β k g k d k 1 g k 2 ) g k + β k d k 1 (6)

因此,为了确保 d k 满足共轭条件,我们只需将其代入(5)式,并求解得到的线性方程。也就是

y k 1 ( ( 1 + β k g k d k 1 g k 2 ) g k β k d k 1 ) = 0

由此得

β k HS1 = g k 2 g k y k 1 g k 2 y k 1 d k 1 HS1 ( g k d k 1 HS1 ) ( g k y k 1 ) (7)

其中 d k 1 替换成了 d k 1 HS1 。然后,将(7)代入(6),得到以下迭代方向:

d k HS1 = ( 1 + β k HS1 g k d k 1 g k 2 ) g k + β k HS1 d k 1 HS1 . (8)

引理2.1由(8)定义的方向 d k HS1 满足充分下降条件和共轭条件。此外,如果步长 α k 1 是通过精确的线搜索计算的,并且 d k 1 HS1 = d k 1 HS ,那么 d k HS1 = d k HS ,其中 d k HS 是HS方法的方向,即

d k HS = g k + β k HS d k 1 HS .

证明 从 d k HS1 的公式可得

g k d k HS1 = g k 2 , y k 1 d k HS1 = 0. (9)

因此,方向 d k HS1 满足充分下降条件和共轭条件。

由于 α k 1 是通过精确的行搜索计算的,因此我们有 g k d k 1 HS1 = 0 。这与(7)、(8)表明如果 d k 1 HS1 = d k 1 HS ,则 β k HS1 = β k HS d k HS1 = d k HS 。证毕。

参数 β k HS1 由以下两个步骤生成:首先通过Schmidt正交化使 d k 满足充分下降条件,然后通过求解相应的线性方程调整得到的 d k 满足共轭条件。显然,上述两个步骤可以交换顺序。也就是说,可以先通过施密特正交化使 d k 满足共轭条件,然后调整得到的 d k 以满足充分下降条件。具体地,为了使(3)定义的 d k 满足共轭条件,通过Schmidt正交化得

d k = g k + β k d k 1 g k y k 1 + β k d k 1 y k 1 y k 1 2 y k 1 .

然后将上述等式代入等式 g k d k = g k 2 ,得

β k HS2 = ( g k y k 1 ) 2 ( g k y k 1 ) ( y k 1 d k 1 HS2 ) y k 1 2 g k d k 1 HS2 (10)

其中 d k 1 替换成了 d k 1 HS2 。于是,得到了一个新的方向

d k HS2 = g k + β k HS2 d k 1 HS2 g k y k 1 + β k HS2 ( d k 1 HS2 ) y k 1 y k 1 2 y k 1 . (11)

引理2.2由(11)定义的方向 d k HS2 满足充分下降条件和共轭条件。此外,如果步长 α k 1 是通过精确的线搜索计算的,并且 d k 1 HS2 = d k 1 HS ,那么 d k HS2 = d k HS

证明 从 d k HS2 的公式可得

g k d k HS2 = g k 2 , y k 1 d k HS2 = 0. (12)

因此,方向 d k HS2 满足充分下降条件和共轭条件。

由于 α k 1 是通过精确的行搜索计算的,因此我们有 g k d k 1 HS2 = 0 。这与(9)、(10)表明如果 d k 1 HS2 = d k 1 HS ,则 β k HS2 = β k HS d k HS2 = d k HS 。证毕。

在下文中,对问题(1)做出以下限制。假设:

(H1) 水平集 Ω = { x R n | f ( x ) f ( x 0 ) } 有界;

(H2) 在 Ω 的某个邻域中,梯度 g ( x ) 是Lipschitz连续的,即存在一个常数 L > 0 ,使得

g ( x ) g ( y ) L x y , x , y N ;

(H3) 函数 f ( x ) Ω 上一致凸,即存在一个常数 μ > 0 ,使得

( g ( x ) g ( y ) ) ( x y ) μ x y 2 , x , y Ω . (13)

假设(H1)与(H2)表明存在正常数 γ 与B使得

g ( x ) γ , x Ω , (14)

x y B , x , y Ω . (15)

在共轭梯度法中,步长 α k 通常由非精确线性搜索确定,比如 Armijo线搜索、Wolfe线搜索、Goldstein线搜索等。本文采用Armijo线搜搜,即步长 α k 是集合 { 1 , ρ , ρ 2 , } 中满足下面不等式的最大的 α

f ( x k + α d k ) f ( x k ) σ α 2 d k 2 , (16)

其中 ρ ( 0 , 1 ) σ > 0

由(13)得

d k 1 y k 1 = d k 1 ( g k g k 1 ) α k 1 d k 1 2 0. (17)

如果 y k 1 = 0 ,则由充分下降性及(17)得 g k 1 = 0 ,这说明 x k 1 是问题(1)的一个稳定点。

基于上面的分析,下面给出求解问题(1)的两项HS共轭梯度法。

TTMHSCG1方法(两项修正HS共轭方法1)

步0. 选取正常数 ν ρ σ 满足 0 < ρ < σ < 1 。选取初始迭代点 x 0 R n ,计算 g 0 。令 d 0 HS1 = g 0 k : = 0

步1. 如果 g k = 0 ,停止否则转步2。

步2. 由Armijo线搜索(16)确定步长 α k

步3. 令 x k + 1 = x k + α k d k ,计算 g k + 1

步4. 如果 g k + 1 y k g k + 1 d k HS1 的符号相同且

| g k + 1 2 y k d k HS1 ( g k + 1 d k HS1 ) ( g k + 1 y k ) | < ν d k HS1

则令 d k + 1 HS1 = g k + 1 ;否则 d k + 1 HS1 由(8)式确定。

步5. 令 k = k + 1 ,返回步1。

类似地,基于 d k HS2 可得另一类HS共轭梯度法。

TTMHSCG2方法(两项修正HS共轭方法2)

步0~步3. 与TTMHSCG1方法一样。

步4. 如果 g k + 1 y k g k + 1 d k HS2 的符号相同且

| ( g k + 1 y k ) ( y k d k HS2 ) y k 2 g k + 1 d k HS2 | < ν d k HS2

则令 d k + 1 HS2 = g k + 1 ;否则 d k + 1 HS2 由(11)式确定。

步5. 令 k = k + 1 ,返回步1。

3. 全局收敛性

本文讨论TTMHSCG1与TTMHSCG2的全局收敛性。下面的引理给出了Armijo线搜索所生成的步长序列 { α k } 存在一个正下界。

引理3.1 如果函数 f ( ) 满足假设(H1)与(H2),步长 α k 由Armijo线搜索确定,则

α k min { 1 , ρ g k 2 ( L + σ ) d k 2 } . (18)

证明 如果 α k 1 ,这说明 α k = α k / ρ 不满足(16),即

f ( x k + α k d k ) > f ( x k ) σ ( α k ) 2 d k 2 . (19)

由中值定理及(H2)得,存在常数 θ k ( 0 , 1 ) 使得

f ( x k + α k d k ) f ( x k ) = α k g ( x k + θ k α k d k ) d k = α k g k d k + ( g ( x k + θ k α k d k ) g k ) d k α k g k 2 + ( α k ) 2 L d k 2 ,

由此式及(19)容易得到(18)成立。证毕。

定理3.1 假设(H1)~(H3)成立,令 { x k } 是算法TTMHSCG1与TTMHSCG2生成的点列,则

lim k g k = 0.

证明 我们只证明了TTMHSCG1的全局收敛性。TTMHSCG2的证明与其类似。事实上,由(H2)可得 y k 1 L x k x k 1 ;由(H3)可得 y k 1 d k 1 HS1 μ α k 1 d k 1 HS1 2 > 0

以下证据分为三种情况。第一种情况,当 g k y k 1 g k d k 1 HS1 有不同的符号时,有

| β k HS1 | = g k 2 | g k y k 1 | g k 2 y k 1 d k 1 HS1 + | ( g k d k 1 HS1 ) ( g k y k 1 ) | | g k y k 1 | y k 1 d k 1 HS1 L g k x k x k 1 μ α k 1 d k 1 HS1 2 = L g k μ d k 1 HS1 ,

由此可得

d k HS1 g k + 2 | β k HS1 | d k 1 HS1 ( 1 + 2 L μ ) g k . (20)

第二种情况,当 | g k 2 y k 1 d k 1 HS1 ( g k d k 1 HS1 ) ( g k y k 1 ) | ν d k 1 ,并且 g k y k 1 g k d k 1 HS1 符号相同时,有

| β k HS1 | = g k 2 | g k y k 1 | | g k 2 y k 1 d k 1 HS1 ( g k d k 1 HS1 ) ( g k y k 1 ) | g k 2 | g k y k 1 | ν d k 1 HS1 L D γ 2 g k d k 1 HS1 ,

其中不等式是根据(14)、(15)得到的。于是

d k HS1 g k + 2 | β k HS1 | d k 1 HS1 ( 1 + 2 L D γ 2 ) g k . (21)

第三种情况,当 | g k 2 y k 1 d k 1 HS1 ( g k d k 1 HS1 ) ( g k y k 1 ) | < ν d k 1 ,并且 g k y k 1 g k d k 1 HS1 符号相同时,由TTMHSCG1的第4步得

d k HS1 = g k . (22)

因此,从(20)~(22),可得

d k τ g k τ γ , (23)

其中 τ = 1 + 2 L max { 1 / μ , D γ 2 }

另一方面,从Armijo线搜索(17),我们有

lim k α k d k = 0.

由这个式子与充分下降条件可得

lim k α k g k lim k α k d k = 0. (24)

下面利用反证法证明结论。假设 lim k g k = 0 不成立,于是存在常数 ϵ > 0 与无穷点列K使得

g k ϵ , k K .

由此式及(18) (23) (24)得

lim k , k K α k = 0 ,

以及

α k min { 1 , ρ ϵ 1 2 ( L + σ ) τ 2 γ 2 } , k K .

显然上面两个式子是矛盾的,因此假设不成立。证毕。

4. 数值实验

本节利用TTMHSCG1与TTMHSCG2求解 [8] 中提出的图像处理的两阶段法。令 X = [ x i , j ] M × N 表示一个有 M × N 个像素的图像,对于每个 ( i , j ) A : = { 1 , 2 , , M } × { 1 , 2 , , N } ,令 V i , j 表示 ( i , j ) 的相邻点,即

V i , j = { ( i , j 1 ) , ( i , j + 1 ) , ( i 1 , j ) , ( i + 1 , j ) } .

同时令 y i , j 表示点 ( i , j ) 的观测像素。下面简要回顾一下两阶段法。在第一阶段,两阶段法利用自适应中值滤波器(AMF)探测出椒盐噪声的位置,并令 N A 表示椒盐噪声的坐标集合;第二阶段,两阶段法通过求解下面的最优化问题填补椒盐噪声点处的像素值。

min u G α ( u ) = ( i , j ) N { ( m , n ) V i , j \ N φ α ( u i , j y m , n ) + 1 2 ( m , n ) V i , j N φ α ( u i , j y m , n ) } ,

其中 α 是正则化参数, φ α 是保边函数, u = [ u i , j ] ( i , j ) N 是一个列向量(按字典序排列)。

图1图2给出利用TTMHSCG1与TTMHSCG2还原Lena (256 × 256)的结果,其中图1的PSNR由13.8402增加到29.0134,图2的PSNR由13.8571增加到28.4543。

Figure 1. The denosing results of TTMHSCG1 (left: original, middle: nosed, right: recovered)

图1. TTMHSCG1去噪效果(左:原始图,中:含噪图,右:去噪图)

Figure 2. The denosing results of TTMHSCG2 (left: original, middle: nosed, right: recovered)

图2. TTMHSCG2去噪效果(左:原始图,中:含噪图,右:去噪图)

图1图2的还原效果可以看出,TTMHSCG1与TTMHSCG2成功地将脉冲噪声去除,因此其在脉冲噪声去噪方面表现良好。

5. 结论

本文设计了两种修正的HS共轭梯度法,它们同时满足充分下降性与共轭条件,同时这两个性质是不依赖于线搜索的。同时本文证明了两种修正HS共轭梯度法的全局收敛性。最后初步的数值实验表明所设计方法在脉冲噪声去噪方面是有效的。

致谢

感谢审稿人对本文提出的宝贵意见。

基金项目

枣庄学院博士科研启动基金、枣庄学院教学改革重点项目(XJG20019)。

文章引用

孙 敏,田茂英. 基于共轭梯度的脉冲噪声去噪算法
Impulse Noise Denoising Algorithm Based on Conjugate Gradient[J]. 图像与信号处理, 2022, 11(03): 125-133. https://doi.org/10.12677/JISP.2022.113014

参考文献

  1. 1. 马昌凤. 最优化方法及其Matlab程序设计[M]. 北京: 科学出版社, 2010: 87-108.

  2. 2. 王宜举, 修乃华. 非线性最优化理论与方法[M]. 北京: 科学出版社, 2011: 112-134.

  3. 3. 戴彧虹, 袁亚湘. 非线性共轭梯度法[M]. 上海: 上海科学技术出版社, 2000: 1-32.

  4. 4. Hager, W.W. and Zhang, H.C. (2005) A New Conjugate Gradient Method with Guaranteed Descent and an Efficient Line Search. SIAM Journal on Optimization, 16, 170-192.
    https://doi.org/10.1137/030601880

  5. 5. Zhang, L., Zhou, W.J. and Li, D.H. (2006) A Descent Modified Polak-Ribière-Polyak Conjugate Gradient Method and Its Global Convergence. IMA Journal of Numerical Analysis, 26, 629-640.
    https://doi.org/10.1093/imanum/drl016

  6. 6. Cheng, W.Y. (2007) A Two-Term PRP-Based Descent Method. Numerical Functional Analysis and Optimization, 28, 1217-1230.
    https://doi.org/10.1080/01630560701749524

  7. 7. Sun, M. and Bai, Q.G. (2011) A New Descent Memory Gradient Method and Its Global Convergence. Journal of Systems Science and Complexity, 24, Article No. 784.
    https://doi.org/10.1007/s11424-011-8150-0

  8. 8. Cai, J.F., Chan, R.H. and Fiore, C.D. (2007) Minimization of a Detail-Preserving Regularization Functional for Impulse Noise Removal. Journal of Mathematical Imaging and Vision, 29, 79-91.
    https://doi.org/10.1007/s10851-007-0027-4

期刊菜单