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

1、1/41初等变分原理初等变分原理最速下降法最速下降法共轭梯度法共轭梯度法数值试验算例数值试验算例数值分析 1112:09从瞎从瞎子下山子下山到最优化方法到最优化方法12:09Science of Better3/41( x, y ) = ( y, x );( tx, y ) = (x, ty )= t ( x, y);( x+ y, z ) = ( x, z ) + ( y, z ); ( x, x) 0, 且且( x, x) = 0 x = 0;(x,y)=|x|2|y|2cos;设设A是是 n 阶对称正定阵阶对称正定阵( Ax, y ) = ( x, Ay )= (A y, x )= (

2、y, Ax );( Ax，x ) 0, 且且( Ax, x) = 0 x = 0 预备知识预备知识I设设 , 记记 ( x , y) = xT ynRyx ,12:094/41预备知识预备知识 II例如例如 f(x1,x2,x3)=x12x22x32梯度梯度: 12( )grad ( ),nTfffxxxf xf x1232212322123221232( )= 22fxfxfxx x xf xx x xx x x 12:09222112221Hessnnnffxx xfffxxx Hessian矩阵矩阵:222223123123222212313123222212312312244Hess4

3、24442x xx x xx x xfx x xx xx x xx x xx x xx x 5/41预备知识预备知识 III2112,()1)( )()() ()()kikijnkkiiinkkiijjiff xxjxxxf xf xxxxxxxR 12( )()()(grad ()()Hess)kkTkkkkTf xff xfxxxxxxxxR 泰勒展式泰勒展式:12:092( )()()()()()kkkkkf xf xxxxxfRfxx 12221112222112212111111221222()()()()2(2)( )()()() ()()()() ()()kkkkkf xf xx

4、xf xfkkkkkkxxxxxfkkxxxkf xf xxxxxxxxxxxxxxxxxR 一阶导数推广一阶导数推广 二阶导数推广二阶导数推广6/41预备知识预备知识 IV12,111( )(,)( , )2 =nnijijiii jif xAx xb xa x xb x gradfAxb222112221Hessnnnffxx xfAffxxx 12:097/41费马引理费马引理:000( ),( ), ()=0 xf xf xxfx 设设是是的的一一个个极极值值点点 且且在在处处导导数数存存在在 则则12:09注释注释: 费马引理的价值在于将极值问题转化为费马引理的价值在于将极值问题转化

5、为非线性方程的求解问题。非线性方程的求解问题。8/41定理定理4.10(初等变分原理初等变分原理) 设设A =(aij )nn为实为实对称正对称正定矩阵定矩阵, 则则 x是二次函数是二次函数 ),(),(21)(xbxAxxf 的极小值点的极小值点 x 是线性方程组是线性方程组 Ax = b 的解的解。 证明证明: 设设 u 是是 Ax = b的解的解 Au = b ),(21),(),(21)()(uAuxbxAxufxf 0)(),(21 uxuxA对任意对任意 xR n , 只须证明只须证明 f (x) f (u) 0),(21)(uAuuf nRbx ,12:099/41 设设 u 使

6、使 f(x) 取极小值。取非零向量取极小值。取非零向量 xR n, 对任意对任意 tR , 有有1( )()( (),)( ,)2g tf utxA utxutxb utx),(2),()(2xAxtxbAutuf 当当 t=0 时时, g(0)= f(u)达到极小值达到极小值, 所以所以 g (0) =0 ,即即( Au b , x ) = 0 Au b = 0所以所以u 是方程组是方程组 Ax = b 的解。的解。12:09瞎子与计算机瞎子与计算机 瞎子瞎子: 能感觉到脚下的坡度能感觉到脚下的坡度(这是海拔函数这是海拔函数在当前点的梯度值在当前点的梯度值),但不知道山上其它点但不知道山上其

7、它点的任何情况的任何情况 计算机计算机: 计算目标函数在该点的信息计算目标函数在该点的信息(如函如函数值和梯度值数值和梯度值), 但不知道其它点的信息但不知道其它点的信息 12:0911/41最速下降法最速下降法(Gradient Descent)从初值点从初值点 x(0) 出发出发,以负梯度方向以负梯度方向 r 为搜索方向为搜索方向在在 x 处处,梯度方向是梯度方向是 f(x) 增长最快方向增长最快方向负梯度方向是负梯度方向是 f(x) 下降最快方向下降最快方向选择步长选择步长 t0, 使使 x(1) = x(0) + t0r 为为 f(x) 极小值点极小值点( )()(grad () ()

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

9、 求点求点: x(1) = x(0) + t0r0 使得使得2/ ),(),()()(0020)0()0(rArtrbAxtxftg 记记12:0913/41解对称正定方程组解对称正定方程组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

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

11、 -x1) )2+ +100( (x2- -x12) )212:09Barzilai-Borwein方法方法 方向方向(最速下降最速下降- -最好方向最好方向) (best rk) 步长步长(上一次精确搜索步长上一次精确搜索步长) (best tk-1) 最好的最好的rk +上一步最好的上一步最好的 tk-1 最好最好参考参考: : Two-Point Step Size Gradient Method12:09f(x1, ,x2)= =100 x12+ +x22最速下降法最速下降法12:09f(x1, ,x2)= =100 x12+ +x22BB方法方法12:0920/41 最速下降法思想简

12、单最速下降法思想简单,但是收敛速度慢。本但是收敛速度慢。本质上是因为负梯度方向函数下降快是局部性质。质上是因为负梯度方向函数下降快是局部性质。全局思想全局思想: :局部思想局部思想: : N 维空间的任意向量可以由维空间的任意向量可以由N个线性无关的个线性无关的向量线性表示。向量线性表示。12:0921/4112:09 共轭梯度法共轭梯度法的关键是构造一组两两共的关键是构造一组两两共轭的方向轭的方向( (即一组线性无关向量即一组线性无关向量) )。巧妙的。巧妙的是是, 共轭方向可以由上次搜索方向和当前的共轭方向可以由上次搜索方向和当前的梯度方法之组合来产生。梯度方法之组合来产生。pk+1 :=

13、 rk+1 + beta*pk The Best of the 20th Century: Editors Name Top 10 Algorithms, SIAM News 现代迭代方法现代迭代方法: Krylov子空间方法子空间方法22/41 共轭共轭 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 ,

14、pn-1 线性无关线性无关两个向量两个向量 p1, p2关于关于 A共轭共轭:12:0923/41 共轭共轭 更加整体地考虑搜索方向的选择更加整体地考虑搜索方向的选择, 选择一组选择一组关于关于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 系系数数计计算算12:0924/41 共轭向量组构造共轭向量组构造 思路思

15、路: 借鉴借鉴Gram-Schmidt正交化过程正交化过程()kkrbAxA 将将残残差差向向量量组组改改造造成成 共共轭轭向向量量组组T(,)(,)ikkiTi ki kiiiiiikkkAp Arrppp rpAppAppr 12:09121113231122112112323(,)(,)(,)(,)(,)(,)(,)(,),ikiiu vu uu vuvu uu viukuuikukuuuvvvuuuuvu 25/41第一步第一步:取初值取初值 x(0)R(n) , 0, 计算计算 r0 = b Ax(0), 若若 | r0| 结束结束; 否则否则 p0r0, k 1, 转第二步转第二步

16、;原始共轭梯度算法原始共轭梯度算法第二步第二步: 计算计算 = (pk-1,b ) / (Apk-1, pk-1) x(k) = x(k-1) + pk-1 ; 第三步第三步: 如果如果 k = n, 则结束则结束; 否则计算否则计算 rk = b Ax(k), 转第四步转第四步;1k 1k 12:09第四步第四步: 如果如果 |rk| , 则则结束结束;否则计算否则计算 = -(Apj , rk ) / (Apj , pj ) , ( j = 0, k-1) pk = rk + ( p0+ + pk-1 ) k k+ 1,转第二步转第二步。j 0 1k 26/41定理定理4.12 A 是是

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

18、p2 , pk-1张张成成k维子空间维子空间,故此类方法称为子空间方法。故此类方法称为子空间方法。1k 1k 12:0927/41重要性质重要性质:1 ( ) (,)=0 ()ijAp pij 12:09(2) ( , )=0 ()ijr rij 01(3) ( ,)=0 (,)kjr pjk 28/41共轭梯度算法共轭梯度算法(The Conjugate Gradient Algorithm) 1. Start:x0:=0,r0:= b, p0:= r0, k:=0.2. Iterate: Until convergence do,(a) alpha := (rk,rk)/(Apk,pk)(

19、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+112:0929/41程序片段程序片段: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*

20、 p; alpha=rho/(q*p) ; x=x+alpha*p ; res = res - alpha *q; rho1= res*res; beta = rho1 / rho ; p = res + beta* p; iter=iter+1; relres=norm(res)/res0;end12:0930/41注注1.111)(),kkkkkkkkrbAxrrAprbAx 11TTTTTTTTTT()()kkkkkkkkkkkkkkkkkkkkkkkprAxp rrprp App Appr rp AAp bApppp 110110( )()( )( )kkkkkkiiixxpxxp 1

21、111111111111TTTTTTT()()()kkkkkkkkkkkkkkkkkkkkkkkkp Arprrrrrrrrr rpAprrrAp 12:091011kkkiiikkkpprrp 31/41注注3. 共轭梯度法误差分析所用范数为共轭梯度法误差分析所用范数为2/1*)()(|xxAxxxxkTkAk AknnAkxxxx|)(2|*011* 注注4. .设设 x* 是方程组是方程组Ax = b的解的解, ,A的特征值为的特征值为： 1 2 n 0共轭梯度法迭代向量共轭梯度法迭代向量xk 误差估计结果为误差估计结果为注注2. 共轭梯度法适用于求解对称正定矩阵方程组。共轭梯度法适用于

22、求解对称正定矩阵方程组。12:0932/41511223223435251622311311311131113113xxxxxx 算例算例133/410,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,0

23、nju ui, 0 = 0, ui, n+1 = 01,11,140iji jijiji juuuuu ( i, ,j = 1, , ,n )算例算例234/41 A = gallery(poisson, n);35/41Axb 总结总结: :12:09直接法直接法现代迭现代迭代法代法古典迭古典迭代法代法36/41A(n 1) = Fn-1Fn-2F1 AFk = I mkekT ( k = 1, 2, , n 1) 称称Fk 为为 Frobenius矩阵矩阵, Fk1 = I + mk ekTA=F1-1F2-1 Fn-1-1 A(n 1)直接方法直接方法: 高斯消元法高斯消元法 L U高斯

24、消元法本质是矩阵的分解。 1111,121nnnmmmL )1()1(2)1(2211211nnnnnaaaaaaU37/41矩阵分解矩阵分解(1)特征值分解特征值分解: A=CDCT, C,D=eig(A)(2) LU分解分解: PA=LU, L,U,P=lu(A)(3)Cholesky分解分解: A=RTR, R=chol(A)12:09(4)QR分解分解: A=QR, Q,R=qr(A)(5)奇异值分解奇异值分解: A=USVH, U,S,V = svd(A) (6) 非负矩阵分解非负矩阵分解Non-negative Matrix Factorization A=WH38/41Demo1I=imread(monalisa.pgm);U,S,V=svd(double(I);s=diag(S);n1=5;Snew=diag(s(1:n1);zeros(size(s,1)-n1,1);figure,imshow(U*Snew*V,)n2=20;Snew=diag(s(1:n2);