数值分析7选主元消

《数值分析》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