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

如果

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

1/2*第一行

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

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

Ref:美國SIAM主席、英國皇家學會院士提的兩個數值代數問題NickTrefethen,TheSmartMoneyIsonNumericalAnalysts,SIAMNews(翻譯)20:08lugui\%BackSlash(DirectMethod)具體參考:C.Moler,Matlab數值計算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)個元素的對應乘積然后求和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的元素計算公式計算順序123456對于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分解計算順序123456例4求矩陣的Crout分解①更新順序:先行后列②列除行不除③舊元素減去其所在行和列前k-1個元素的對應乘積然后求和Doolittle分解：①更新順序:先列后行②行除列不除③舊元素減去其所在行和列前k-1個元素的對應乘積然后求和Crout分解：三對角矩陣(tridiagonalmatrix)

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

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

K>0為上對角線,

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

分解20:08大一統的框架：只要會Doolittle分解Crout分解(轉置)三對角矩陣分解對稱矩陣LDLT

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

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

A=LLT,L為下三角矩陣。思考:嘗試推導Cholesky分解的公式,并編程實現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)非負矩陣分解:

A=WH20:08矩陣分解的應用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