数值分析11(共轭梯度法)

1、1/34初等变分原理初等变分原理最速下降法最速下降法共轭梯度法共轭梯度法数值试验算例数值试验算例数值分析 11*2/34( x, x) 0, 当且仅当当且仅当x=0时等号成立时等号成立;( x, y ) = ( y, x );(kx+l y, z ) = k( x, z ) + l( y, z ); (x,y)=|x|2|y|2cos。预备知识预备知识设设 , 则实数则实数 ( x , y)= xT y= x1y1+ +xnynnRyx ,*3/34设设A是是 n 阶对称正定阵阶对称正定阵(Ax, x) 0, 当且仅当当且仅当x=0时等号成立时等号成立;( Ax, y ) = ( x, Ay

2、)= (A y, x )= ( y, Ax );(kAx+lA y, z ) = k( Ax, z ) + l(Ay, z ). 预备知识预备知识矩阵矩阵A正定正定,如果对于任意非零向量如果对于任意非零向量x满足满足 xTAx0.*4/34预备知识预备知识 例如例如 f(x1,x2,x3)=x12x22x32梯度梯度(多元函数的一阶导数信息多元函数的一阶导数信息): 12( )grad ( ),nTfffxxxf xf x 1232212322123221232( )= 22fxfxfxx x xf xx x xx x x *思考思考:多元函数的二阶导数信息？多元函数的二阶导数信息？5/34预

3、备知识预备知识 (21)2( )()()(| )kinf xkkkiiixf xf xxxOxx 22(grad ()( )()()(| )kkkkTf xf xxxxxfOx 泰勒展式泰勒展式(数值分析的基石数值分析的基石):*2( )()()() )kkkkf xf xxxOxfxx 12()(211222)( )()()()(|-| )kkf xfkkkxxkxf xf xxxxxOx x 6/34单变量函数极值点单变量函数极值点(费马引理费马引理):00( )()=0 xf xfx 设设 是是的的一一个个极极值值点点的的必必要要条条件件是是。*注释注释: 费马引理的价值在于将极值问题转

4、化为费马引理的价值在于将极值问题转化为非线性方程的求解问题。非线性方程的求解问题。7/3412,111( )(,)( , ) = 2nnijijiii jif xAx xb xa x xb x1111grad0niiinniinia xbfAxba xb *多变量函数极值点多变量函数极值点:00( ) grad ()0 xf xf x 设设 是是的的一一个个极极值值点点的的必必要要条条件件是是。8/34定理定理4.10(初等变分原理初等变分原理) 设设A =(aij )nn为实对称正为实对称正定矩阵定矩阵, 则则 x是二次函数是二次函数 ),(),(21)(xbxAxxf 的极小值点的极小值点

5、 x 是线性方程组是线性方程组 Ax = b 的解的解。 证明证明: 设设 u 是是 Ax = b的解的解 Au = b ),(21),(),(21)()(uAuxbxAxufxf 0)(),(21 uxuxA对任意对任意 xR n , 只须证明只须证明 f (x) f (u) 0),(21)(uAuuf nRbx ,*9/34 设设 u 是是 f(x) 极小值点。取非零向量极小值点。取非零向量 xR n, 对任意对任意 tR , 有有1( )()( (),)( ,)2g tf utxA utxutxb utx),(2),()(2xAxtxbAutuf 当当 t=0 时时, g(0)= f(u

6、)达到极小值达到极小值, 所以所以 g (0) =0,即即( Au b , x ) = 0 Au b = 0所以所以u 是方程组是方程组 Ax = b 的解。的解。*10/34最速下降方向最速下降方向从初值点从初值点 x(0) 出发出发,以负梯度方向以负梯度方向 r 为搜索方向为搜索方向在在 x 处处,梯度方向是梯度方向是 f(x) 增长最快方向增长最快方向负梯度方向是负梯度方向是 f(x) 下降最快方向下降最快方向选择步长选择步长 t0, 使使 x(1) = x(0) + t0r 为为 f(x) 极小值点极小值点( )()(grad () ()()(grad (),kkTkkTkkkkkkf

7、 xf xf xxxf xtf xddxxt 其其中中为为单单位位向向量量为为步步长长。最速下降方向最速下降方向: r = f = b Ax 22( , ) | | cos,a baba ba bab 其其中中表表示示向向量量 和和 的的夹夹角角。*11/34求解得求解得 t0 = ( r0 , r0) / (Ar0 , r0)0),(),()(000)0( rArtrbAxtg为选取最佳步长为选取最佳步长 t0,令令取初值点取初值点 x(0), 取负梯度方向取负梯度方向 r0 = b A x(0)(min)(0)0(00)0(rtxfrtxfRt 求点求点: x(1) = x(0) + t0

8、r0 使得使得2/ ),(),()()(0020)0()0(rArtrbAxtxftg 记记*精确搜索步长精确搜索步长12/34解对称正定方程组解对称正定方程组Ax = b 的的最速下降算法最速下降算法:第一步第一步: 取初值取初值 x(0)R(n) , 0,计算计算 r0 = b Ax(0), k 0;第二步第二步: 计算计算 tk = (rk ,rk ) / (Ark , rk) x(k+1) = x(k) + tk rk , rk+1 = b Ax(k+1)第三步第三步: k k+ 1, 如果如果 |rk| ,转第二步转第二步; 否则输出否则输出 x(k), 结束。结束。*注释注释: 最

9、速下降算法思想简单且容易实现最速下降算法思想简单且容易实现,是是求解无约束优化问题的经典方法。求解无约束优化问题的经典方法。最好最好 + + 最好最好 = = 最好最好 ? 方向方向( (最速下降最速下降) ) (best rk) 步长步长( (精确搜索精确搜索) ) (best tk) x(k+1) = x(k) + tkrk 是否最好是否最好 ?*Rosenbrock 方程方程f( (x1,x2)=()=(1- -x1) )2+ +100( (x2- -x12) )2*Rosenbrock 方程方程f( (x1,x2)=()=(1- -x1) )2+ +100( (x2- -x12) )2

10、*Barzilai-Borwein方法方法 方向方向(最速下降最速下降- -最好方向最好方向) (best rk) 步长步长(上一次精确搜索步长上一次精确搜索步长) (best tk-1) 最好的最好的rk +上一步最好的上一步最好的 tk-1 最好最好参考文献参考文献: : Two-Point Step Size Gradient Method*f(x1, ,x2)= =100 x12+ +x22最速下降法最速下降法*f(x1, ,x2)= =100 x12+ +x22BB方法方法*19/34 最速下降法思想简单最速下降法思想简单,但是收敛速度慢。本但是收敛速度慢。本质上是因为负梯度方向函数

11、下降快是局部性质。质上是因为负梯度方向函数下降快是局部性质。全局思想全局思想: :局部思想局部思想: : 在在n维空间中维空间中,任意任意n个线性无关的向量构成个线性无关的向量构成n维空间的基。换句话说维空间的基。换句话说, ,n维空间中任意向量均维空间中任意向量均可以由这组基线性表示。可以由这组基线性表示。*20/34* 共轭梯度法共轭梯度法的关键是构造一组两两共的关键是构造一组两两共轭的方向轭的方向( (即张成即张成n维空间的基维空间的基) )。巧妙的是。巧妙的是共轭方向可以由上次搜索方向和当前的梯共轭方向可以由上次搜索方向和当前的梯度方向线性组合产生。度方向线性组合产生。pk+1 :=

12、rk+1 + beta*pk The Best of the 20th Century: Editors Name Top 10 Algorithms, SIAM News 现代迭代方法现代迭代方法: 子空间方法子空间方法21/34 共轭共轭 A是是 n 阶对称正定矩阵阶对称正定矩阵,非零向量非零向量 p1, p2Rn(Ap1, p2)=0n 个向量个向量 p0, p1 , pn-1 共轭共轭:(Api , pj )=0(ij; i, j = 0,1,n-1 )非零向量非零向量 p0, p1 , pn-1 Rnp0, p1 , pn-1 关于关于 A 共轭共轭 p0, p1 , pn-1 线性

13、无关线性无关两个向量两个向量 p1, p2关于关于 A共轭共轭:*22/34 系数计算系数计算 更加整体地考虑搜索方向的选择更加整体地考虑搜索方向的选择, 选择一组关选择一组关于于A共轭的向量共轭的向量:n 个向量个向量 p0, p1 , pn-1 共轭共轭:(Api , pj )=0(ij; i, j = 0,1,n-1 )1100*nniiiiiixpb AxAp 我我们们有有且且 = =10TT*TT, nkkkikikkkipp bp Axp App Ap 对对任任意意的的TTkkkkp bp Ap : : 计计算算系系数数*23/34 共轭向量组构造共轭向量组构造 思路思路: Gra

14、m-Schmidt正交化过程正交化过程()kkrbAxA 一一系系列列最最速速下下降降方方向向改改造造成成 共共轭轭向向量量组组(,)(,)ikikkiikiAp rpApppr *121113231122123312121(,)(,)(,)(,)(,)(,)(,)(,)ikiiu vu uu vuvuu viuuuukuki kuuuuuuuuvvvv 24/34第一步第一步:取初值取初值 x(0) =0, 0, 计算计算 r0 = b Ax(0), 若若 | r0| 结束结束; 否则否则 p0r0, k 1, 转第二步转第二步;原始共轭梯度算法原始共轭梯度算法第二步第二步: 计算计算 =

15、(pk-1,b ) / (Apk-1, pk-1) x(k) = x(k-1) + pk-1 ; (张成张成k维子空间维子空间) 第三步第三步: 如果如果 k = n, 则结束则结束; 否则计算否则计算 rk = b Ax(k), 转第四步转第四步;1k 1k *第四步第四步: 如果如果 |rk| , 则则结束结束;否则计算否则计算 = (Apj , rk ) / (Apj , pj ) , ( j = 0, k-1) pk = rk + ( p0+ + pk-1 ) k k+1,转第二步转第二步。j 0 1k 25/34定理定理4.12 A 是是 n 阶对称正定矩阵阶对称正定矩阵, p0,

16、p1 , pn-1 是关于是关于 A 共轭的向量组共轭的向量组, 任取任取 x(0)Rn , 计算计算 = (pk-1, b) / (Apk-1, pk-1)x(k) = x(k 1) + pk-1 ( k = 1,2, n )则有则有 Ax(n) = b。 共轭梯度方法理论上有限步终止共轭梯度方法理论上有限步终止, 故仅被故仅被视为直接法视为直接法, 所以在其后的所以在其后的20年没有受到重视。年没有受到重视。1972年共轭梯度法被首次作为迭代法研究年共轭梯度法被首次作为迭代法研究,很有很有新意。第新意。第k步迭代步迭代, p0, p1 , pk-1张成张成k维子空间维子空间,故此类方法称为

17、子空间方法。故此类方法称为子空间方法。1k 1k *26/34共轭梯度算法共轭梯度算法(The Conjugate Gradient Algorithm) 1. Start:x0,r0:= b-Ax0, p0:= r0, k:=0.2. Iterate: Until convergence do,(a) alpha := (rk,rk)/(Apk,pk)(b) xk+1 := xk + alpha* pk(c) rk+1 := rk alpha*Apk(d) beta := (rk+1,rk+1)/(rk,rk)(e) pk+1 := rk+1 + beta*pk k:=k+1*27/3422

18、6253uv 例例1 用共轭梯度法求解用共轭梯度法求解0000603( ),xrp 00006633502122663325TTTTr rp Ap 1010750021570603( )( )xxp 127510002124762263253 rrAp 112772277110012160496633 TTTTr rr r 28/34226253uv 例例1 用共轭梯度法求解用共轭梯度法求解18017491611004921207496123 prp 12127724247711180180114949120120494971102225 TTTTr rp Ap 211018074971110

19、512074941( )( ) xxp 18012749721111024120749220250 rrAp 29/34程序片段程序片段:Matlab Code : CG方法function x,relres,iter=conjgrad(A,b,x,nmax,tol)res0=norm(b-A*x); iter=0;relres=1;res=b; p=res; rho1=res*res; while relrestol & iternmax rho=rho1; q=A* p; alpha=rho/(q*p) ; x=x+alpha*p ; res = res - alpha *q; r

20、ho1= res*res; beta = rho1 / rho ; p = res + beta* p; iter=iter+1; relres=norm(res)/res0;end*30/34511223223435251622311311311131113113xxxxxx 算例算例131/340,0,1(0, )( ,0)( ,1)(1, )0 xxyyuux yuyu xu xuy Possion方程:令令 h = 1/(n+1) , xj= jh, yj = jh ( i , , j = 0, ,1, , , , n+1 )记记 ui,j= u(xi , yj ), ( i , j = 0,1, , n+1 )线性方程组线性方程组041, 11, 1 jijiijjijiuuuuu( i, ,j = 1, , ,n )u0, j = 0, 1,0nju ui, 0 = 0, ui, n+1 = 01,11,140iji jijiji juuuuu ( i, ,j = 1, , ,n )算例算例232/34 A = gallery(poisson, n);33/34重要性质重要性质:1 ( ) (,)=0 ()ijAp pij *(2) ( , )=0 ()ijr rij 01(3) ( ,)=0 (,)kjr pjk34/34证明