您好,欢迎访问三七文档
当前位置:首页 > 商业/管理/HR > 商业计划书 > 矩阵QR分解MATLAB自编程序
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正交化分解方法对矩阵进行QR分解时,所论矩阵必须是列满秩矩阵。为方便与后续基于Householder变换的QR分解法对比,这里以方阵为例.3.1.2算法1写出矩阵的列向量;2列向量按照Schmidt正交化正交;3得出矩阵的𝑸′,𝑹′;4对𝑸′的列向量单位化得到𝑸,𝑹′的每行乘𝑸′每列的模得𝑹.3.1.3流程图输入方阵A列向量Schmidt正交化整理得Q’R’Q’列向量单位化得QR’每行乘Q’每列的模得R结束3.1.4程序function[X,Q,R]=QRDecomsch(A,b)%方阵的QR的Gram-Schmidt正交化分解法,并用于求解AX=b方程组[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)-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,使得𝑨=𝑸𝑹设w∈Cn是一个单位向量,令2HHI则称H是一个Householder矩阵或Householder变换。则对于任意的Cn存在Householder矩阵H,使得Hx=-au。其中3.2.2算法第一步,将矩阵A按列分块写成A=(α1,α2,…,αn).如果α1≠0,则可得,存在n阶householder矩阵H1使得H1α1=-a1e1,|a1|=||α1||,e1∈Cn于是有H1A=(H1α1,H1α2,…,H1αn)=11*0naA如果α1=0,则直接进行下一步,此时相当于取H1=In,而a1=0.第二步,将矩阵An-1按列分块写成An-1=(αi,α2,…,αn-1)。如果α1≠0,则可得,存在n-1阶householder矩阵H’2使得H’2α2=-a2e1,|a2|=||α2||,e1∈Cn于是有H’2An-1=(H’2α2,…,H’2αn-1)=22*0naA此时,令H2=2100'TH则H2是n阶Householder矩阵,且使H2H1A=122**0*0naaA如果α1=0,则直接进行下一步。第三步,对n-2阶矩阵继续进行类似的变换,如此下去,之多在第n-1步,我们可以找到Householder矩阵H1,H2,…,Hn-1使得axHn-1…H2H1A=12***0**00*000'nnaa令Q=Hn-1…H2H1,则Q是酉矩阵之积,从而必有酉矩阵并且A=QR3.2.3流程图3.2.4程序function[X,Q,R]=QRDecomhouse(A,b)%用Householder变换将方阵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)~=0w=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)=b(n)/R(n,n);fori=n-1:-1:1X(i)=(b(i)-sum(R(i,i+1:n).*X(i+1:n)'))/R(i,i);endelseX=[];end3.3运行验证与分析A=12101−1211−1−1112−31,b=[1011]T,求A得QR分解,并解方程组A𝒙=𝒃.3.3.1MATLAB自带QR分解函数计算在MATLAB中调用:[Q0,R0]=qr(A),运行结果如下:clearallA=[1,2,1,-1;1,0,2,1;1,-1,1,2;-1,1,-3,1];b=[1011]';X0=inv(A)*b[Q0,R0]=qr(A)Q0=-0.50000.8165-0.0577-0.2828-0.5000-0.00000.17320.8485-0.5000-0.4082-0.7506-0.14140.50000.4082-0.63510.4243R0=-2.00000-3.5000-0.500002.4495-0.8165-1.2247001.4434-1.90530001.2728X0=2.0000-0.0000-1.000003.3.2Schmidt正交化计算结果clearallA=[1,2,1,-1;1,0,2,1;1,-1,1,2;-1,1,-3,1];b=[1011]';[X1,Q1,R1]=QRDecomsch(A,b)X1=2.00000.0000-1.00000.0000Q1=0.50000.8165-0.0577-0.28280.500000.17320.84850.5000-0.4082-0.7506-0.1414-0.50000.4082-0.63510.4243R1=2.000003.50000.500002.4495-0.8165-1.2247001.4434-1.90530001.27283.3.3Householder变换运行结果clearallA=[1,2,1,-1;1,0,2,1;1,-1,1,2;-1,1,-3,1];b=[1011]';[X2,Q2,R2]=QRDecomhouse(A,b)X2=2.00000.0000-1.00000.0000Q2=0.50000.81650.0577-0.28280.50000.0000-0.17320.84850.5000-0.40820.7506-0.1414-0.50000.40820.63510.4243R2=2.000003.50000.500002.4495-0.8165-1.224700-1.44341.90530001.27283.3.4数据分析由1、2、3可知,Gram-Schmidt正交化法与Householder变换法所得Q、R以及方程组解与MATLAB自带函数结果一致,可见程序运行结果正确。
本文标题:矩阵QR分解MATLAB自编程序
链接地址:https://www.777doc.com/doc-1897048 .html