矩阵QR分解MATLAB自编程序.pdf

《矩阵QR分解MATLAB自编程序.pdf》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、3矩阵的QR分解3.1基于Schmidt正交化的QR分解3.1.1原理Schmidt正交化方法是矩阵的QR分解最常用的方法.主要依据下面的两个结论:结论1设?是n阶实非奇异矩阵,则存在正交矩阵?和实非奇异上三角矩阵?使?有??分解;且除去相差一个对角元素的绝对值（模）全等于1的对角矩阵因子外,分解是唯一的.结论2设A是m×n实矩阵,且其n个列向量线性无关,则A有分解A=QR,其中Q是m×n实矩阵,且满足QHTQ=E,R是n阶实非奇异上三角矩阵该分解除去相差一个对角元素的绝对值（模）全等于1的对角矩阵因子外是唯一的.用Schmidt正交化分解

2、方法对矩阵进行QR分解时，所论矩阵必须是列满秩矩阵。为方便与后续基于Householder变换的QR分解法对比，这里以方阵为例.3.1.2算法1写出矩阵的列向量;2列向量按照Schmidt正交化正交;3得出矩阵的?′,?′;4对?′的列向量单位化得到?，?′的每行乘?′每列的模得?.3.1.3流程图输入方阵A列向量Schmidt正交化整理得Q’R’’Q列向量单位化得Q’’R每行乘Q每列的模得R结束3.1.4程序function[X,Q,R]=QRDecomsch(A,b)%方阵的QR的Gram-Schmidt正交化分解法，并用于求解AX=b

3、方程组[m,n]=size(A);ifm~=ndisp('不满足QR分解要求');endQ=zeros(m,n);X=zeros(n,1);R=zeros(n);fork=1:nR(k,k)=norm(A(:,k));ifR(k,k)==0break;endQ(:,k)=A(:,k)/R(k,k);forj=k+1:nR(k,j)=Q(:,k)'*A(:,j);A(:,j)=A(:,j)-R(k,j)*Q(:,k);endendifnargin==2b=Q'*b;X(n)=b(n)/R(n,n);fori=n-1:-1:1X(i)=(b(i

4、)-sum(R(i,i+1:n).*X(i+1:n)'))/R(i,i);endelseX=[];end3.2基于Householder变换的QR分解3.2.1原理设A为任一n阶方阵，则必存在n阶酉矩阵Q和n阶上三角阵R，使得?=??nH设w∈C是一个单位向量，令HI2则称H是一个Householder矩阵或Householder变换。则对于任nax意的C存在Householder矩阵H，使得Hx=-au。其中3.2.2算法第一步，将矩阵A按列分块写成A=(α1，α2，…，αn）.如果α1≠0，则可得，存在n阶householde

5、r矩阵H1使得nH1α1=-a1e1，

6、a1

7、=

8、

9、α1

10、

11、，e1∈Ca*1于是有H1A=（H1α1，H1α2，…，H1αn）=0An1如果α1=0，则直接进行下一步，此时相当于取H1=In，而a1=0.第二步，将矩阵An-1按列分块写成An-1=（αi,α2，…,αn-1）。如果’’α1≠0，则可得，存在n-1阶householder矩阵H2使得H2α2=-a2e1，

12、a2n

13、=

14、

15、α2

16、

17、，e1∈C’’2’a2*于是有H2An-1=（Hα2，…，H2αn-1）=0An2T10此时，令H2=0'H

18、2则H2是n阶Householder矩阵，且使a**1HHA=210*a20An2如果α1=0，则直接进行下一步。第三步，对n-2阶矩阵继续进行类似的变换，如此下去，之多在第n-1步，我们可以找到Householder矩阵H1，H2，…，Hn-1使得a***10a**2Hn-1…H2H1A=00*000'nn令Q=Hn-1…H2H1，则Q是酉矩阵之积，从而必有酉矩阵并且A=QR3.2.3流程图3.2.4程序function[X,Q,R]=QRDecomhouse(A,b)%用H

19、ouseholder变换将方阵A分解为正交Q与上三角矩阵R的乘积，并用于求解AX=b方程组[n,n]=size(A);E=eye(n);X=zeros(n,1);R=zeros(n);P1=E;fork=1:n-1%构造w,使Pk=I-2ww's=-sign(A(k,k))*norm(A(k:n,k));R(k,k)=-s;ifk==1w=[A(1,1)+s,A(2:n,k)']';elsew=[zeros(1,k-1),A(k,k)+s,A(k+1:n,k)']';R(1:k-1,k)=A(1:k-1,k);endifnorm(w)~=0

20、w=w/norm(w);endP=E-2*w*w';A=P*A;P1=P*P1;R(1:n,n)=A(1:n,n);endQ=P1';ifnargin==2b=P1*b;X(n)