1、1/41初等變分原理初等變分原理最速下降法最速下降法共軛梯度法共軛梯度法數(shù)值試驗(yàn)算例數(shù)值試驗(yàn)算例數(shù)值分析 1118:05從瞎從瞎子下山子下山到最優(yōu)化方法到最優(yōu)化方法18:05Science 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;設(shè)設(shè)A是是 n 階對(duì)稱正定陣階對(duì)稱正定陣( Ax, y ) = ( x, Ay )= (A y, x )= (

2、y, Ax );( Ax，x ) 0, 且且( Ax, x) = 0 x = 0 預(yù)備知識(shí)預(yù)備知識(shí)I設(shè)設(shè) , 記記 ( x , y) = xT ynRyx ,18:054/41預(yù)備知識(shí)預(yù)備知識(shí) II例如例如 f(x1,x2,x3)=x12x22x32梯度梯度: 12( )grad ( ),nTfffxxxf xf x1232212322123221232( )= 22fxfxfxx x xf xx x xx x x 18:05222112221Hessnnnffxx xfffxxx Hessian矩陣矩陣:222223123123222212313123222212312312244Hess4

3、24442x xx x xx x xfx x xx xx x xx x xx x xx x 5/41預(yù)備知識(shí)預(yù)備知識(shí) III2112,()1)( )()() ()()kikijnkkiiinkkiijjiff xxjxxxf xf xxxxxxxR 12( )()()(grad ()()Hess)kkTkkkkTf xff xfxxxxxxxxR 泰勒展式泰勒展式:18:052( )()()()()()kkkkkf xf xxxxxfRfxx 12221112222112212111111221222()()()()2(2)( )()()() ()()()() ()()kkkkkf xf xx

4、xf xfkkkkkkxxxxxfkkxxxkf xf xxxxxxxxxxxxxxxxxR 一階導(dǎo)數(shù)推廣一階導(dǎo)數(shù)推廣 二階導(dǎo)數(shù)推廣二階導(dǎo)數(shù)推廣6/41預(yù)備知識(shí)預(yù)備知識(shí) IV12,111( )(,)( , )2 =nnijijiii jif xAx xb xa x xb x gradfAxb222112221Hessnnnffxx xfAffxxx 18:057/41費(fèi)馬引理費(fèi)馬引理:000( ),( ), ()=0 xf xf xxfx 設(shè)設(shè)是是的的一一個(gè)個(gè)極極值值點(diǎn)點(diǎn) 且且在在處處導(dǎo)導(dǎo)數(shù)數(shù)存存在在 則則18:05注釋注釋: 費(fèi)馬引理的價(jià)值在于將極值問(wèn)題轉(zhuǎn)化為費(fèi)馬引理的價(jià)值在于將極值問(wèn)題轉(zhuǎn)化

5、為非線性方程的求解問(wèn)題。非線性方程的求解問(wèn)題。8/41定理定理4.10(初等變分原理初等變分原理) 設(shè)設(shè)A =(aij )nn為實(shí)為實(shí)對(duì)稱正對(duì)稱正定矩陣定矩陣, 則則 x是二次函數(shù)是二次函數(shù) ),(),(21)(xbxAxxf 的極小值點(diǎn)的極小值點(diǎn) x 是線性方程組是線性方程組 Ax = b 的解的解。 證明證明: 設(shè)設(shè) u 是是 Ax = b的解的解 Au = b ),(21),(),(21)()(uAuxbxAxufxf 0)(),(21 uxuxA對(duì)任意對(duì)任意 xR n , 只須證明只須證明 f (x) f (u) 0),(21)(uAuuf nRbx ,18:059/41 設(shè)設(shè) u 使

6、使 f(x) 取極小值。取非零向量取極小值。取非零向量 xR n, 對(duì)任意對(duì)任意 tR , 有有1( )()( (),)( ,)2g tf utxA utxutxb utx),(2),()(2xAxtxbAutuf 當(dāng)當(dāng) t=0 時(shí)時(shí), g(0)= f(u)達(dá)到極小值達(dá)到極小值, 所以所以 g (0) =0 ,即即( Au b , x ) = 0 Au b = 0所以所以u(píng) 是方程組是方程組 Ax = b 的解。的解。18:05瞎子與計(jì)算機(jī)瞎子與計(jì)算機(jī) 瞎子瞎子: 能感覺(jué)到腳下的坡度能感覺(jué)到腳下的坡度(這是海拔函數(shù)這是海拔函數(shù)在當(dāng)前點(diǎn)的梯度值在當(dāng)前點(diǎn)的梯度值),但不知道山上其它點(diǎn)但不知道山上其

7、它點(diǎn)的任何情況的任何情況 計(jì)算機(jī)計(jì)算機(jī): 計(jì)算目標(biāo)函數(shù)在該點(diǎn)的信息計(jì)算目標(biāo)函數(shù)在該點(diǎn)的信息(如函如函數(shù)值和梯度值數(shù)值和梯度值), 但不知道其它點(diǎn)的信息但不知道其它點(diǎn)的信息 18:0511/41最速下降法最速下降法(Gradient Descent)從初值點(diǎn)從初值點(diǎn) x(0) 出發(fā)出發(fā),以負(fù)梯度方向以負(fù)梯度方向 r 為搜索方向?yàn)樗阉鞣较蛟谠?x 處處,梯度方向是梯度方向是 f(x) 增長(zhǎng)最快方向增長(zhǎng)最快方向負(fù)梯度方向是負(fù)梯度方向是 f(x) 下降最快方向下降最快方向選擇步長(zhǎng)選擇步長(zhǎng) t0, 使使 x(1) = x(0) + t0r 為為 f(x) 極小值點(diǎn)極小值點(diǎn)( )()(grad () ()

8、()(grad (),kkTkkTkkkkkkf xf xf xxxf xtf xddxxt 其其中中為為單單位位向向量量為為步步長(zhǎng)長(zhǎng)。最速下降方向最速下降方向: r = f = b Ax ( , ) |cos,a baba ba bab其其中中表表示示向向量量 和和 的的夾夾角角。18:0512/41求解得求解得 t0 = ( r0 , r0) / (Ar0 , r0)0),(),()(000)0( rArtrbAxtg為選取最佳步長(zhǎng)為選取最佳步長(zhǎng) t0 ，令令取初值點(diǎn)取初值點(diǎn) x(0), 取負(fù)梯度方向取負(fù)梯度方向 r0 = b A x(0)(min)(0)0(00)0(rtxfrtxfRt

9、 求點(diǎn)求點(diǎn): x(1) = x(0) + t0r0 使得使得2/ ),(),()()(0020)0()0(rArtrbAxtxftg 記記18:0513/41解對(duì)稱正定方程組解對(duì)稱正定方程組Ax = b 的的最速下降算法最速下降算法:第一步第一步: 取初值取初值 x(0)R(n) , 0,計(jì)算計(jì)算 r0 = b Ax(0), k 0;第二步第二步: 計(jì)算計(jì)算 tk = (rk ,rk ) / (Ark , rk) x(k+1) = x(k) + tk rk ; rk+1 = b Ax(k+1) ;第三步第三步: k k+ 1, 如果如果 |rk| ,轉(zhuǎn)第二步轉(zhuǎn)第二步; 否則輸出否則輸出 x(k

10、), 結(jié)束。結(jié)束。18:05注釋注釋: 最速下降算法思想簡(jiǎn)單且容易實(shí)現(xiàn)最速下降算法思想簡(jiǎn)單且容易實(shí)現(xiàn),是是求解無(wú)約束優(yōu)化問(wèn)題的經(jīng)典方法。求解無(wú)約束優(yōu)化問(wèn)題的經(jīng)典方法。最好最好 + + 最好最好 = = 最好最好 ? 方向方向( (最速下降最速下降) ) (best rk) 步長(zhǎng)步長(zhǎng)( (精確搜索精確搜索) ) (best tk) x(k+1) = x(k) + tkrk 是否最好是否最好 ?18:05Rosenbrock 方程方程f( (x1,x2)=()=(1- -x1) )2+ +100( (x2- -x12) )218:05Rosenbrock 方程方程f( (x1,x2)=()=(1-

11、 -x1) )2+ +100( (x2- -x12) )218:05Barzilai-Borwein方法方法 方向方向(最速下降最速下降- -最好方向最好方向) (best rk) 步長(zhǎng)步長(zhǎng)(上一次精確搜索步長(zhǎng)上一次精確搜索步長(zhǎng)) (best tk-1) 最好的最好的rk +上一步最好的上一步最好的 tk-1 最好最好參考參考: : Two-Point Step Size Gradient Method18:05f(x1, ,x2)= =100 x12+ +x22最速下降法最速下降法18:05f(x1, ,x2)= =100 x12+ +x22BB方法方法18:0520/41 最速下降法思想簡(jiǎn)

12、單最速下降法思想簡(jiǎn)單,但是收斂速度慢。本但是收斂速度慢。本質(zhì)上是因?yàn)樨?fù)梯度方向函數(shù)下降快是局部性質(zhì)。質(zhì)上是因?yàn)樨?fù)梯度方向函數(shù)下降快是局部性質(zhì)。全局思想全局思想: :局部思想局部思想: : N 維空間的任意向量可以由維空間的任意向量可以由N個(gè)線性無(wú)關(guān)的個(gè)線性無(wú)關(guān)的向量線性表示。向量線性表示。18:0521/4118:05 共軛梯度法共軛梯度法的關(guān)鍵是構(gòu)造一組兩兩共的關(guān)鍵是構(gòu)造一組兩兩共軛的方向軛的方向( (即一組線性無(wú)關(guān)向量即一組線性無(wú)關(guān)向量) )。巧妙的。巧妙的是是, 共軛方向可以由上次搜索方向和當(dāng)前的共軛方向可以由上次搜索方向和當(dāng)前的梯度方法之組合來(lái)產(chǎn)生。梯度方法之組合來(lái)產(chǎn)生。pk+1 :=

13、 rk+1 + beta*pk The Best of the 20th Century: Editors Name Top 10 Algorithms, SIAM News 現(xiàn)代迭代方法現(xiàn)代迭代方法: Krylov子空間方法子空間方法22/41 共軛共軛 A是是 n 階對(duì)稱正定矩陣階對(duì)稱正定矩陣,非零向量非零向量 p1, p2Rn(Ap1, p2)=0n 個(gè)向量個(gè)向量 p0, p1 , pn-1 共軛共軛:(Api , pj )=0(ij; i, j = 0,1,n-1 )非零向量非零向量 p0, p1 , pn-1 Rnp0, p1 , pn-1 關(guān)于關(guān)于 A 共軛共軛 p0, p1 ,

14、pn-1 線性無(wú)關(guān)線性無(wú)關(guān)兩個(gè)向量?jī)蓚€(gè)向量 p1, p2關(guān)于關(guān)于 A共軛共軛:18:0523/41 共軛共軛 更加整體地考慮搜索方向的選擇更加整體地考慮搜索方向的選擇, 選擇一組選擇一組關(guān)于關(guān)于A共軛的向量共軛的向量:n 個(gè)向量個(gè)向量 p0, p1 , pn-1 共軛共軛:(Api , pj )=0(ij; i, j = 0,1,n-1 )1100*nniiiiiixpb AxAp我我們們有有且且 = =10TT*TT, nkkkikikkkipp bp Axp App Ap 對(duì)對(duì)任任意意的的TTkkkkp bp Ap 系系數(shù)數(shù)計(jì)計(jì)算算18:0524/41 共軛向量組構(gòu)造共軛向量組構(gòu)造 思路思

15、路: 借鑒借鑒Gram-Schmidt正交化過(guò)程正交化過(guò)程()kkrbAxA 將將殘殘差差向向量量組組改改造造成成 共共軛軛向向量量組組T(,)(,)ikkiTi ki kiiiiiikkkAp Arrppp rpAppAppr 18:05121113231122112112323(,)(,)(,)(,)(,)(,)(,)(,),ikiiu vu uu vuvu uu viukuuikukuuuvvvuuuuvu 25/41第一步第一步:取初值取初值 x(0)R(n) , 0, 計(jì)算計(jì)算 r0 = b Ax(0), 若若 | r0| 結(jié)束結(jié)束; 否則否則 p0r0, k 1, 轉(zhuǎn)第二步轉(zhuǎn)第二步

16、;原始共軛梯度算法原始共軛梯度算法第二步第二步: 計(jì)算計(jì)算 = (pk-1,b ) / (Apk-1, pk-1) x(k) = x(k-1) + pk-1 ; 第三步第三步: 如果如果 k = n, 則結(jié)束則結(jié)束; 否則計(jì)算否則計(jì)算 rk = b Ax(k), 轉(zhuǎn)第四步轉(zhuǎn)第四步;1k 1k 18:05第四步第四步: 如果如果 |rk| , 則則結(jié)束結(jié)束;否則計(jì)算否則計(jì)算 = -(Apj , rk ) / (Apj , pj ) , ( j = 0, k-1) pk = rk + ( p0+ + pk-1 ) k k+ 1,轉(zhuǎn)第二步轉(zhuǎn)第二步。j 0 1k 26/41定理定理4.12 A 是是

17、n 階對(duì)稱正定矩陣階對(duì)稱正定矩陣, p0, p2 , pn-1 是關(guān)于是關(guān)于 A 共軛的向量組共軛的向量組, 任取任取 x(0)Rn , 計(jì)算計(jì)算 = (b, pk-1) / (Apk-1, pk-1)x(k) = x(k 1) + pk-1 ( k = 1,2, n )則有則有 Ax(n) = b。 共軛梯度方法理論上有限步終止共軛梯度方法理論上有限步終止, 故僅僅故僅僅被視為是直接法被視為是直接法, 所以在其后的所以在其后的20年沒(méi)有受到年沒(méi)有受到重視。重視。1972年共軛梯度法被首次作為迭代法研年共軛梯度法被首次作為迭代法研究究,很有新意。第很有新意。第k步迭代迭代步迭代迭代, p0,

18、p2 , pk-1張張成成k維子空間維子空間,故此類方法稱為子空間方法。故此類方法稱為子空間方法。1k 1k 18:0527/41重要性質(zhì)重要性質(zhì):1 ( ) (,)=0 ()ijAp pij 18:05(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+118:0529/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;

20、q=A* 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;end18:0530/41注注1.111)(),kkkkkkkkrbAxrrAprbAx 11TTTTTTTTTT()()kkkkkkkkkkkkkkkkkkkkkkkprAxp rrprp App Appr rp AAp bApppp 110110( )()( )( )kkkkkkiiixxpx

21、xp 1111111111111TTTTTTT()()()kkkkkkkkkkkkkkkkkkkkkkkkp Arprrrrrrrrr rpAprrrAp 18:051011kkkiiikkkpprrp 31/41注注3. 共軛梯度法誤差分析所用范數(shù)為共軛梯度法誤差分析所用范數(shù)為2/1*)()(|xxAxxxxkTkAk AknnAkxxxx|)(2|*011* 注注4. .設(shè)設(shè) x* 是方程組是方程組Ax = b的解的解, ,A的特征值為的特征值為： 1 2 n 0共軛梯度法迭代向量共軛梯度法迭代向量xk 誤差估計(jì)結(jié)果為誤差估計(jì)結(jié)果為注注2. 共軛梯度法適用于求解對(duì)稱正定矩陣方程組。共軛梯度

22、法適用于求解對(duì)稱正定矩陣方程組。18:0532/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,

23、 1,0nju ui, 0 = 0, ui, n+1 = 01,11,140iji jijiji juuuuu ( i, ,j = 1, , ,n )算例算例234/41 A = gallery(poisson, n);35/41Axb 總結(jié)總結(jié): :18:05直接法直接法現(xiàn)代迭現(xiàn)代迭代法代法古典迭古典迭代法代法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

24、 U高斯消元法本質(zhì)是矩陣的分解。 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)18:05(4)QR分解分解: A=QR, Q,R=qr(A)(5)奇異值分解奇異值分解: A=USVH, U,S,V = svd(A) (6) 非負(fù)矩陣分解非負(fù)矩陣分解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);