共轭梯度法是介于最速下降法与牛顿法之间的一个方法。
它仅需要利用一阶导数信息,但克服了最速下降法收敛慢的缺点,又避免了牛顿法需要存储和计算Hesse 矩阵并求逆的缺点。
共轭梯度法不仅是解决大型线性方程组最有用的发发之一,也是解大型非线性最优化问题最有效的算法之一。
共轭梯度法最早是由计算数学家Hestenes 和几何学家Stiefel 在20世纪50年代初为求解线性方程组Ax b =n x R ∈而各自独立提出的。
他们合作的文章被公认为共轭梯度法的奠基之作。
该文详细讨论了求解线性方程组的共轭梯度法的性质以及它和其他方法的关系。
在A 为对称正定阵时,上述线性方程组等价于最优化问题1min 2n T Tx R x Ax b x ∈-。
由此,Hestenes 和Stiefel 的方法也可视为求二次函数。
提到最优化问题,这里首先介绍最速下降法。
考虑线性方程组Ax b =的求解问题,其中A 是给定的n 阶对称正定矩阵,b 是给定的n 维向量。
为此我们定义二次泛函()2T T x x Ax b x Φ=-对最速下降法做一简单分析就会发现,负梯度方向虽从局部来看是最佳的下山方向,但从整体来看并非最佳。
这就促使人们去寻求更好的下山方向,当然,我们自然希望每步确定新的下山方向付出的代价不要太大。
共轭梯度法就是根据这一意思设计的,其具体计算过程如下:给定初始向量0x ,第一步仍选负梯度方向为下山方向,即0p = 0r ,于是有0000TTr p Ap α=,1000x x p α=+,11r b Ax =-对以后各步,例如,第k+1步(1k ≥),下山方向就不再取k r ,而是在过点k x 由向量k r 和1k p -所张成的二维平面21{:,}k k k x x r p R πξηξη-==++∈内找出使函数Φ下降最快的方向作为新的下山方向k p 。
考虑Φ在2π上的限制:1111(,)()()()2()T T k k k k k k k k k k k k x r p x r p A x r p b x r p ψξηξηξηξηξη----=Φ++=++++-++直接计算可得:12()T T T k k k k k k r Ar r Ap r r ψξηξ-∂=+-∂1112()T Tk k k k r Ap p Ap ψξηη---∂=+∂ 其中最后一式用到了10T k k r p -=,这可由k r 的定义直接验证。
令ψξ∂∂=ψη∂∂=0,即知Φ在2π内有唯一的极小点~001k k k x x r p ξη-=++,其中0ξ和0η满足方程组001T T T k k k k k r Ar r Ap r r ξη-+=(1)01011T Tk k k k r Ap p Ap ξη---+=0 (2)上式蕴含着0k r ≠必有0ξ≠0,因此我们可取~101()k k k k p x x r p ηξξ-=-=+作为新的下山方向。
显然,这是在平面2π内可得到的最佳下山方向,令010k ηβξ-=,则由(2)式得到1111T k k k T k k r Ap p Ap β----=-。
这样确定的k p 满足1Tkk p Ap -=0,即所谓的k p 与1k p -是相互共轭的。
k p 确定以后k α的确定仍用公式k α=T k kT k k r p p Ap ,然后计算1k k k k x x p α+=+。
总结上述讨论,可得如下计算公式:k α=T k kT k kr p p Ap 1k k k k x x p α+=+ 11k k r b Ax ++=-1T k kk T k kr Ap p Ap β+=- 11k k k p r p β++=+在实际计算中,常将上述公式进一步简化,从而得到一个形式上更为简单而且对称的计算公式。
首先来简化1k r +的计算公式:11k k r b Ax ++=-=()k k k b A x p α-+=k k k r Ap α-(3)因为k Ap 在计算k α时已经求出,所以,计算1k r +时可以不必将1k x +代入方程去计算,而是由(3)得到。
再来简化k α和k β的计算公式。
我们需要用到下面的关系式:1110T T T k k k k k k r r r p r p +-+===,k=1,2, (4)这些关系式的证明包含在方程组(1)(2)的证明中。
由(3)(4)导出1111111()T T T k k k k k k k kkr Ap r r r r r αα+++++=-=-,111()TTTk k k k k k k kkp Ap p r r p r αα+=-=。
由此可得:T k k k T k k r r p Ap α=,11T k k k T k kr r r r β++=最优化方法的收敛性是算法研究领域的基本问题。
一种算法是否具备应用价值,取决于是否具备局部或全局收敛性,以及收敛速度的快慢。
共轭梯度法不同的方向修正公式,以及采用何种线搜索策略来确定步长,都取决于它的收敛性质。
下面重点简述两个共轭梯度法收敛性的成果。
1.FR 方法:早期对FR 方法的收敛性研究是建立在精确线搜索基础上的。
Powell 在精确线搜索下,得到FR 方法的一个很不利的性质,即:如果FR 方法在某一步产生一个很小的步长,则相继的许多步长也可能非常小。
同时,Powell 还给出了FR 方法最简单的全局有效性分析。
这些分析解释了为什么FR 方法在数值表现上并不十分满意。
尽管FR 方法可能收敛速度很慢,Zoutendijk 证明了采取精确线搜索的FR 方法对一般非凸函数总是收敛的。
在实际计算中,人们通常采用非精确线搜索,而不是使用精确线搜索。
最早的非精确线搜索的全局收敛结果是由A1.Baali 在1985年给出的,他证明了使用参数12σ<的强Wolfe 线搜索的FR 方法一定满足充分下降条件,而且是全局收敛的。
有人后将A1.Baali 的结果推广到了12σ=。
通过分析使用推广的Wolfe 线搜索的FR 方法,通过考虑相邻两个迭代点列,发现只要FR 方法的每个搜索方向下降,那么任意相邻两个迭代点列中至少有一个使得充分下降条件成立,从而较简单的证明了12σ≤时采用强Wolfe 线搜索的算法的全局收敛性。
对于12σ>,有人举出反例表明,FR 方法可能产生一个上升方向而导致失败,并且提出了一种广义线搜索策略,证明了FR 方法在广义Wolfe 线搜索的全局收敛性。
2.PRP 方法:PRP 方法是目前认为数值表现最好的共轭梯度算法之一。
当算法产生一个小步长时,由PRP 方法定义的搜索方向k d 自动靠近负梯度方向,从而较为有效地避免FR 方法可能连续产生小步长的缺陷。
基于此,Powell 证明了当步长1k k k s x x +=-趋于0时PRP 方法的全局收敛性。
进一步可以得到,采取精确线搜索的PRP 方法对一致凸函数的全局收敛性。
但对一般非凸函数,Powell 举出了一个三维的反例表明,即使按Curry 原则选取步长因子,PRP 方法可能在6个点附近进行循环,而其中任意一点都不是目标函数的稳定点。
在二维时,Powe11证明了采取Curry 原则的PRP 方法对一般非凸函数的收敛性。
如果使用非精确线搜索,戴彧虹举例表明,即使对于一致凸函数,采用强wolfe 线搜索,且参数充分接近于0,PRP 方法都有可能产生一个上升方向。
如果再要求每一个搜索方向下降,那么PRP 方法对凸函数是收敛的。
对一般非凸函数,Powell 建议限制PRP 方法中的参数PRP k β为非负max{,0}PRP k k ββ=这样做的目的是避免当║k d ║很大时,相邻两个搜索方向会产生相反的趋势。
Gilbert 和Nocedal 接受了Powell 的上述建议,并在适当的线搜索条件下,得到了上述修正PRP 方法对一般非凸函数的全局收敛性结果。
然而,Gilbert 和Nocedal 也举出反例表明,即使对一致凸的目标函数, PRP 也可能为负。
于是,Grippo 和Lucidi 设计了一种Armijo 型的线搜索,并证明了原始PRP 方法在该线搜索下,对一般非凸函数的全局收敛性。
他们的结果是应用当步长k s 很小时,PRP 方法的方向接近最速下降方向的性质。
戴彧虹和袁亚湘得出PRP 方法一个新的性质,即证明了取常数步长因子的PRP 方法在每次迭代都产生一个下降方向,而且全局收敛。
但是,他的这种常数步长因子的选取依赖于Lipschitz 常数L ,而L 在实际计算中很难预先估计,因此并不容易实现。
共轭梯度法在实际生活中的应用很广泛,对数学科学的研究、生产生 活有着重要意义,这里我们以一例进行说明分析。
例:用FR 共轭梯度法解极小化问题121222122123)(min x x x x x x f --+=这里我们将用到共轭梯度法相关的数学模型:共轭梯度法是一个典型的共轭方向法,它的每一个共轭方向是相互共轭的而这些搜索方向k d 仅仅是负梯度方向k g -与上一次迭代的搜索方向1-k d 的组合。
因此,存储量少计算方便记 11--+-=k k k k d g d β左乘G d Tk 1-并使得01=-k T k Gd d 1111----=k T k k Tk k Gd d Gdg β (Hestenes-Stiefel )我们可以将数值公式1111----=k T k k T k k Gd d Gd g β改写为()()()()eves Fletcher g g g g Wolfe Crowder g g d g g g k Tk kTkk k T k k kTkk Re ,111111-=---=------β注意到对于正定二次函数k k k r b Gx g ∆=-= 其中k r 是方程组b Gx k =的残量,以及kTk kT k k T k k T k k k k k k Gd d r r Gd d d g Gd r r =-==-+αα,1 下面给出关于正定二次函数极小化的共轭梯度法。
算法步骤:步1(初始步)给出0,0>εx ;计算b Gx r -=00令0:,00=-=k r d 步2 如果ε≤k r 停止 步3 计算kk k k k Tk kT k k d x x Gd d r r αα+==+1kk k k k Tk k T k k kk k k d r d r r r r Gd r r ββα+-==+=+++++11111步4 令,1+=k k 转步2通过Matlab 软件进行计算,函数从初始点T x )4,2(0-=开始迭代,计算两次,得出最优解T x )1,1(2=。