數(shù)值分析7選主元消

《數(shù)值分析》6列主元消元法直接三角分解法特殊矩陣的分解矩陣分解的應(yīng)用例1例1經(jīng)過第一輪Gauss消元增廣矩陣為

如果

ε非常小,則1/ε非常大,矩陣不接近奇異(行列式為-1-ε),也沒有成千上萬次運(yùn)算造成的誤差累積。20:08部分選主元的高斯消元法:絕對(duì)值小的主元可能產(chǎn)生麻煩,選取消元后的低階矩陣內(nèi)絕對(duì)值最大的元素作為主元,高斯消元法具有更好的數(shù)值穩(wěn)定性。例2交換第一行和第三行20:08第二行減去(-1/2)*第一行第三行減去

1/2*第一行

交換第二行和第三行第三行減去(-1/2)*第二行

部分選主元的高斯消元法是否可以表示為L(zhǎng)U分解的形式？LU選主元高斯消元法矩陣形式:PA=LUL?UMatlab:[L,U,P]=lu(A)計(jì)算機(jī)的早期開創(chuàng)者馮.諾依曼和圖靈等,主要關(guān)心的是用高斯消元法求解大規(guī)模線性方程組時(shí)舍入誤差累積是否會(huì)使結(jié)果誤差很大,這一點(diǎn)上,開始的研究結(jié)果是非常悲觀的,但是不久后計(jì)算實(shí)踐驚人地表明了此方法的穩(wěn)定性和準(zhǔn)確性。威爾金森提出了向后誤差分析方法,系統(tǒng)地研究了矩陣計(jì)算的誤差,對(duì)(選主元)高斯消元法的這種良好特性做出了解釋和證明。由于威爾金森在發(fā)展數(shù)值計(jì)算技術(shù)和方法上的杰出貢獻(xiàn),他在1970年被授予圖靈獎(jiǎng),他在圖靈獎(jiǎng)?lì)C獎(jiǎng)大會(huì)上作了題為SomeCommentsfromanumericalanalyst的演講。

Ref:美國(guó)SIAM主席、英國(guó)皇家學(xué)會(huì)院士提的兩個(gè)數(shù)值代數(shù)問題NickTrefethen,TheSmartMoneyIsonNumericalAnalysts,SIAMNews(翻譯)20:08lugui\%BackSlash(DirectMethod)具體參考:C.Moler,Matlab數(shù)值計(jì)算20:08矩陣滿足什么條件存在LU分解？矩陣滿足什么條件存在PA=LU分解？參考：

或NecessaryAndSufficientConditionsForExistenceoftheLUFactorizationofanArbitraryMatrix2=u11,4=u12,4=u13,2=u143

=l21u11,2=

l31u11,

4=

l41u11l21=3/2,l31=2/2=1,

l41=

4/2=2a22

=l21u12+u22,a23=

l21u13+u23,

a24=l21u14+u24u22=a22-l21u12=-3,u23=a23-l21u13=6,u24=a24-l21u14=3a32

=l31u12+l32u22,a42=l41u12+l42u22

l32=(a32-l31u12)/u22=0,l42=(a42-l41u12)/u22=2①更新順序:先行后列②列除行不除③舊元素減去其所在行和列前(k-1)個(gè)元素的對(duì)應(yīng)乘積然后求和20:08矩陣的Doolittle分解

A=LU,L為單位下三角矩陣,U為上三角矩陣.矩陣的Doolittle分解a11=u11,···,a1n=

u1na21

=l21u11,···,an1=

ln1u11l21u12+u22=a22,···,l21u1n+u2n=a2n

u22=a22-l21u12,···,

u2n=a2n-l21u1n

l31u12+l32u22=a32,···,ln1u12+ln2u22=an2

l32=(a32-l31u12)/u22,···,ln2=(an2-ln1u12)/u22矩陣L和矩陣U的元素計(jì)算公式計(jì)算順序123456對(duì)于akj

和akj其中j,i=k,…,n。例3求矩陣的Doolittle分解程序片段:Matlabcode:Doolittle分解functionA=Doolittle(A)[n,n]=size(A);fork=1:nA(k,k)=A(k,k)-A(k,1:k-1)*A(1:k-1,k);if(A(k,k)==0)fprintf(’Error:A(%d,%d)=0!\n’,i,i);return;end

A(k,k+1:n)=A(k,k+1:n)-A(k,1:k-1)*A(1:k-1,k+1:n);A(k+1:n,k)=(A(k+1:n,k)-A(k+1:n,1:k-1)*A(1:k-1,k))/A(k,k);end矩陣Crout分解計(jì)算順序123456例4求矩陣的Crout分解①更新順序:先行后列②列除行不除③舊元素減去其所在行和列前k-1個(gè)元素的對(duì)應(yīng)乘積然后求和Doolittle分解：①更新順序:先列后行②行除列不除③舊元素減去其所在行和列前k-1個(gè)元素的對(duì)應(yīng)乘積然后求和Crout分解：三對(duì)角矩陣(tridiagonalmatrix)

現(xiàn)實(shí)應(yīng)用中矩陣通常具有特殊的結(jié)構(gòu)。如何針對(duì)矩陣的特殊結(jié)構(gòu),定制快速的矩陣分解算法:稀疏矩陣(sparsematrix)帶狀矩陣(bandedmatrix)參考:1.特殊矩陣,陳景良2.矩陣分析與應(yīng)用,張賢達(dá)對(duì)稱矩陣(symmetrymatrix)托普利茲矩陣(Toeplitzmatrix)三對(duì)角矩陣分解三對(duì)角矩陣的非零系數(shù)在主對(duì)角線上和兩條次對(duì)角線上。三對(duì)角矩陣和更一般的帶狀矩陣來源于偏微分方程的差分方法等。

如何產(chǎn)生對(duì)角矩陣？T=diag([111])+diag([22],1)+diag([33],-1);函數(shù)diag(V,K)生成對(duì)角矩陣,其中向量V的元素出現(xiàn)在矩陣第K條對(duì)角線。特別的,K=0為主對(duì)角矩陣,

K>0為上對(duì)角線,

K<0為下對(duì)角線。LU分解三角矩陣的LU分解:=三對(duì)角矩陣單位下三角陣上三角陣三對(duì)角矩陣的LU分解(k=2,3,···,n)對(duì)稱矩陣:AT=A正定矩陣:AT=A且A的各階順序主子式大于零(1)對(duì)稱矩陣LDLT

分解20:08大一統(tǒng)的框架：只要會(huì)Doolittle分解Crout分解(轉(zhuǎn)置)三對(duì)角矩陣分解對(duì)稱矩陣LDLT

分解對(duì)稱正定矩陣Cholesky分解(2)對(duì)稱正定矩陣的BBT分解(Cholesky分解)例5.

計(jì)算如下矩陣的Cholesky分解解:Matlab:R=chol(A)Demochol([13;32])矩陣的Cholesky分解

A=LLT,L為下三角矩陣。思考:嘗試推導(dǎo)Cholesky分解的公式,并編程實(shí)現(xiàn)MatrixFactorizationJungleMATLAB中矩陣分解(1)特征值分解:AV=VD

%[V,D]=eig(A)(2)LU分解:PA=LU%[L,U,P]=lu(A)(3)Cholesky分解:A=RTR%R=chol(A)QR分解:

A=QR%[Q,R]=qr(A)(5)

奇異值分解:A=USVH%[U,S,V]=SVD(A)(6)非負(fù)矩陣分解:

A=WH20:08矩陣分解的應(yīng)用1: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*Snew*V',[])n2=20;Snew=diag([s(1:n2);zeros(size(s,1)-n2,1)]);figure,imshow