1、三、四章內(nèi)容提要 典型例題分析,數(shù)值分析典型例題 II, ,11:38,化難為易,化繁為簡(jiǎn),化繁為簡(jiǎn),初等行變換不改變方程組的解,1.交換矩陣的第i行與第j行的位置,2.以非零數(shù)k乘以矩陣的第i行的每個(gè)元素,3.把矩陣的第i行的每個(gè)元素的k倍加到第j行的對(duì)應(yīng)元素上去,A(n 1) = Fn-1Fn-2F1 A,其中Fk 為 Frobenius矩陣。,A=F1-1F2-1 Fn-1-1 A(n 1),直接方法: 高斯消元法,L U,高斯消元法本質(zhì)是矩陣的分解。,矩陣分解(Top 10 Algorithms),(1)特征值分解: A=CDC, C,D=eig(A),(3) LU分解: PA=LU,

2、 L,U,P=lu(A),(4)Cholesky分解: A=LLT, R=chol(A),11:38,(2)奇異值分解: A=USV , U,S,V = svd(A),(5) 非負(fù)矩陣分解 Learning the Parts of Objects by Non-negative Matrix Factorization, Nature, 1999,Demo1,I=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

3、*Snew*V,) n2=20; Snew=diag(s(1:n2);zeros(size(s,1)-n2,1); figure, imshow(U*Snew*V,),11:38,迭代法思想:,11:38,Iterate:To say or doagainoragain and again,迭代背后的思想是一種與傳統(tǒng)思維模式截然不同的方式,傳統(tǒng)思維方式往往希望一遍做好,一次成功;但是迭代開(kāi)發(fā)意味著反復(fù)地做,不斷地根據(jù)反饋進(jìn)行調(diào)整。,11:38,迭代格式構(gòu)造,收斂條件(局部vs全局),中止準(zhǔn)則,加速(松弛思想),Aitken加速方法 超松弛加速方法,11:38,共軛梯度法的關(guān)鍵是構(gòu)造一組兩兩共軛

4、的方向(第k步迭代生成共軛方向張成k維子空間)。巧妙的是共軛方向可以由上次搜索方向和當(dāng)前的梯度方向組合產(chǎn)生。,現(xiàn)代迭代方法 (Top 10 Algorithms),最速下降法思想簡(jiǎn)單,但收斂速度慢。本質(zhì)上因?yàn)樨?fù)梯度方向函數(shù)下降快是局部性質(zhì)。,復(fù)雜性:,高斯消元法共用乘法和除法次數(shù)為n3/3+ n2-n/3,常用記號(hào)O表示是多少階的,則O (n3/3)。,注釋2: 復(fù)雜性對(duì)估計(jì)求解大型方程組所需的時(shí)間有用。例如在一臺(tái)特定的計(jì)算機(jī)上求解n=500個(gè)方程的方程組所需的時(shí)間我們可以通過(guò)求解一個(gè)n=50個(gè)方程的方程組得到一個(gè)很好的猜測(cè),即對(duì)用掉的時(shí)間按比例放大1000倍。,注釋1:按O的記法,把n的不同

5、次冪相加的結(jié)果僅保留了最高次冪,因?yàn)樽罡叽蝺鐩Q定了當(dāng)n趨近無(wú)窮時(shí)的極限形態(tài)。換而言之,對(duì)于大的n ,低階項(xiàng)對(duì)算法的執(zhí)行時(shí)間的估計(jì)沒(méi)有太大影響, 僅需要近似估計(jì)執(zhí)行時(shí)間時(shí)可以忽略不計(jì)。,常用的范數(shù):,Hilbert矩陣條件數(shù): for i=1:10 c(i)=cond(hilb(i),2);%vander(1:i) end,plot(1:10,c),范數(shù)(全局),問(wèn)題的好與壞,算法的快與慢,| X(k+1) X* | |B| | X(k) X* |,范數(shù)的威力和魅力:,ekT= 0 0 1 0 0 ,= I mkekT,Fk1 = I + mk ekT,( k = 1, 2, , n 1),F1

6、-1F2-1 Fn-1-1 =, ,Fk1 Fj1 = I + mk ekT+mj ejT kj,F1-1F2-1 Fn-1-1 = I + m1 e1T+mn-1 en-1T,例1.,ekTmj =0, kj,例2.設(shè)A為對(duì)稱矩陣。高斯消元法一步后,A約化為,證明 B 也是對(duì)稱矩陣。,約化主元不為零的判斷,定理3.1 約化主元 的充分必要條件是矩陣A各階順序主子式 不等于零。即,例4.嚴(yán)格對(duì)角占優(yōu)矩陣的約化主元ak,k(k-1) 0 (k=1,n) 。,例3.嚴(yán)格主對(duì)角占優(yōu)矩陣一定是非奇異的。,嚴(yán)格對(duì)角占優(yōu)矩陣:高斯消元法中約化主元不等于零,Jacobi方法, GS方法和SOR方法收斂。,思

7、考: 若A是嚴(yán)格對(duì)角占優(yōu)矩陣, 當(dāng)0 w =1時(shí), SOR方法收斂。,例5.Ax=b,其中A對(duì)稱正定,問(wèn)解此方程的Jacobi迭代法是否一定收斂？,對(duì)稱正定矩陣:直接法高斯消元法中約化主元大于零,迭代法GS方法,SOR方法,最速下降方法和共軛梯度法收斂。,例6.,例7.,例8.,病因是條件數(shù)大,病癥是什么呢?,例9.矩陣的Doolittle分解,A = LU, L為單位下三角矩陣,U為上三角矩陣。,a11= u11, , a1n= u1n,a21 = m21u11, , an1 = mn1u11,11:38,對(duì)A的元素aij ,當(dāng) jk 和 ik+1時(shí),m21u12+ u22=a22, , m

8、21u1n+ u2n=a2n,u22=a22 - m21u12, , u2n=a2n- m21u1n,m31u12+ m32u22=a32, , mn1u12+ mn2u22=an2,m32=(a32- m31u12)/u22, , mn2=(an2- mn1u12)/u22,11:38,矩陣L和矩陣U的元素計(jì)算公式,計(jì)算次序,11:38,更新順序: 先行后列 列除行不除 舊元素減去所在行和列前k-1分量點(diǎn)乘的加和,Doolittle分解：,更新順序: 先列后行 行除列不除 舊元素減去所在行和列前k-1分量點(diǎn)乘的加和,Crout分解：,求矩陣的Doolittle分解,11:38,三對(duì)角矩陣分解

9、,11:38,步驟I:,三對(duì)角矩陣分解 A = L U,( k = 2, 3, , n ),11:38,下三角方程組 LY = f,步驟II:,上三角方程組 UX = Y,步驟III:,11:38,求解三對(duì)角方程組的追趕法,11:38,對(duì)稱正定矩陣的 Cholesky分解,(1) 對(duì)于非零向量, xTAx總是正的; (2) A的順序主子式全大于零; (3) A的特征值全為正實(shí)數(shù)。 help chol,對(duì)稱正定矩陣:,1) 序列收斂 2) 迭代矩陣譜半徑小于1 3) 迭代矩陣特征值全小于1 4) 迭代矩陣范數(shù)小于1 5) 反證法,迭代法收斂證明思路:,例10.設(shè)A是一個(gè)可逆矩陣, 矩陣序列滿足

10、Xk+1=Xk(2I A Xk )，（k =0,1,2,） 則當(dāng) 時(shí),證明：由Xk+1=Xk(2I A Xk )，得 I AXk+1 = I A Xk(2I A Xk )= (I A Xk )2 于是 I AXk =(I A Xk -1)2 =(I A Xk -2)22 = ,例11.,思路:,(1) A1 = B ( I + R + R2 + )； (2)任意給定n階矩陣X0,由迭代格式 Xk+1 = Xk R + B ( k = 0,1,2, ) 產(chǎn)生的矩陣序列 Xk 收斂到矩陣A-1； (3)對(duì)矩陣序列 Xk , 有誤差估計(jì)式,例12.設(shè)A是n階可逆矩陣,有A的一個(gè)近似逆B,令 R=I

11、AB。如果 | R | q 1,試證明,收斂的w取值范圍。,例13. 方程組Ax=b, 其中A是對(duì)稱正定陣,討論迭代格式,思路:,例14.,思路:,續(xù)例14.,思路:,例15.,思路:,定理4.1 對(duì)任意的f和任意的初始向量X(0)迭代法 X(k+1) =B X(k) +f 收斂的充分必要條件是,例16.,思路:,例17.,思路: -1/2a1/2收斂, -1/2a1時(shí)系數(shù)矩陣正定。,例18.,例18.,例18.,參考:Writting Fast Matlab Codes,1.中小規(guī)模線性方程組 x = Ab; % Solves A*x = b,2.(超)大規(guī)模線性方程組(Preconditi

12、oned cg),作業(yè),題目1: 從理論角度(復(fù)雜度和收斂性)比較各 種方法的優(yōu)劣。,11:38,題目2: 從數(shù)值角度比較各種方法的性能(公平),題目3: 科研或生活中遇到的線性方程組？,提示:,參考代碼: Matlab代碼Iterative Solver,11:38,1. 如何生成方程組 A = gallery(poisson, n); b=ones(size(A,1),1); x0=zeros(size(A,1),1);,2. 如何生成方程組 矩陣市場(chǎng): /MatrixMarket/ UF Sparse Matrix Collection: http:

13、//research/sparse/matrices/,提示:,11:38,3. 更具挑戰(zhàn)和趣味的例子 圖像編輯 參考: Matlab代碼Possion 圖像復(fù)原 參考: Deblurring Images: Matrices, Spectra, and Filtering Google PageRank 參考: Numerical Computing with MATLAB 數(shù)據(jù)挖掘和模式識(shí)別 參考: Matrix Methods in Data Mining and Pattern Recognition 高光譜圖像解混 參考: Sparse Unmixing of Hyperspect