QR分解的数值效果报告.doc

基于Householder、givens、gram-schmit、modified gram-schmit方法的QR分解数值效果分析报告班级：0811101学号：081110215姓名：桑树浩目 录1. 摘要32. 关键词33. 引言34. 理论基础45. 数值实验86. 总结与分析197.参考文献与附注.20QR分解数值效果分析报告基于数值试验的gram-schmit方法，modified gram-schmit方法，householder变换，givens变换计算矩阵QR分解。一摘要：由于工程需要引入了QR分解。可以利用多次正交变换（householder变换、givens变换、gram-schmit方法，modified gram-schmit方法）来实现QR分解，各个方法各有优劣和其意义，在此就分析一下其各自效果。经过本次数值实验知道givens变换所用时间较长，但是精确度还是比较高的，householder变换所用时间最短，gram-schmit方法，modified gram-schmit方法所用时间较短，精度最高。可以根据不同需要选用不同算法。二关键词：householder变换 givens变换 gram-schmit方法 modified gram-schmit方法 QR分解 数值效果比较三引言：很多情况下在工程中抽象出来的数学方程组是超定的，没有精确解，这样就需要找一个在某种意义下最接近精确解的解。设A是m*n的实矩阵，|Ax-b|2=|Q(Ax-b)|2,这样min|Ay-b|2就等价与min|QAx-Qb|2,为方便求解，需要QA是上三角矩阵，这样引入QR分解就比较求解方便。可以利用多次正交变换（householder变换、givens变换、gram-schmit方法，modified gram-schmit方法）来实现QR分解1. 豪斯霍尔德变换（Householder transformation）又称初等反射（Elementary reflection），最初由A.C Aitken在1932年提出1。Alston Scott Householder在1958年指出了这一变换在数值线性代数上的意义。这一变换将一个向量变换为由一个超平面反射的镜像，是一种线性变换。其变换矩阵被称作豪斯霍尔德矩阵，在一般内积空间中的类比被称作豪斯霍尔德算子。超平面的法向量被称作豪斯霍尔德向量。2. Givens变换：欲把一个向量中许多分量化为零，可以用Householder变换，例如前面所讲到的把一个向量中若干相邻分量化为零。如果只将其中一个分量化为零，则就采用Givens变换，它有如下形式：其中.易证是一个正交阵。Givens变换又叫平面旋转变换。3.Gram-Schmidt正交化方法：在线性代数中，如果内积空间上的一组向量能够组成一个子空间，那么这一组向量就称为这个子空间的一个基。GramSchmidt正交化提供了一种方法，能够通过这一子空间上的一个基得出子空间的一个正交基，并可进一步求出对应的标准正交基。Gram-Schmidt正交化的基本想法，是利用投影原理在已有正交基的基础上构造一个新的正交基。4.modified Gram-schmidt正交化方法：在Gram-schmidt正交化过程中产生了r(i,j)=q(i)*v(j)+ q(i)*r(k,j)*q(k)k=1,2,3,j-1后面这一项在计算过程中由于机器精度有限会产生计算误差，故可以优化为r(I,j)=q(i)*v1(j)v1(j)=v(j)+ r(k,j)*q(k)k=1,2,3,j-1四理论依据：1.householder变换：任意的一个向量x，经过一个正交变换H=I-2*w*w后总可以变为一个与之范数相等的另一个向量Hx。如上图中所示，记v=Hx-xw=v/|v|2H=I-2*w*w上述H既为所要求的householder变换。具体操作时将需要变化的矩阵的每一列当做一个向量，第一列变为除了第一个元素都为0的向量第i列变为除了前i个元素都为0的向量.为了保证在每次变化时不改变已经变好的0元素，第i次只变化每列第i个元素到第n个元素。每个变化矩阵的形式是I 0;0 H。算法如下：q=diag(ones(n,1); 初始化q向量for i=1:n-1 x=a(i:n,i); y=; y(1)=|x|2; for j=2:n-i+1 y(j)=0; end w=(y-x)/|y-x|2; b=diag(ones(n-i+1,1)-2*w*w; c=zeros(n); e=diag(ones(i-1,1); c(1:i-1,1:i-1)=e; c(i:n,i:n)=b; h=c; a=h*a; q=h*q;endq=q-1;r=a; %a即为变化之后的矩阵，即R2.givens变换：任意的一个向量x被其作用后，y=G*x,有yi=c*xi+s*xkyk=-s*xi+c*xkyj=xj, j=i,k。 将一个矩阵QR分解时，依次将该矩阵的第一列的第2个、第3个.第n个元素变为0，第二列第3、4.n个元素变为0第n-1列元素的第n个元素变为0.具体算法如下： q=diag(ones(n,1);%初始化Qfor i=1:n for j=i+1:n if a(j,i)=0; g=diag(ones(n,1); c=a(i,i)/sqrt(a(i,i)2+a(j,i)2); s=a(j,i)/sqrt(a(i,i)2+a(j,i)2); g(i,i)=c; g(j,j)=c; g(j,i)=-s; g(i,j)=s; a=g*a; q=g*q; end endendr=a;q=q-1;3Gram-Schmidt正交化基本想法，是利用投影原理在已有正交基的基础上构造一个新的正交基。设。是上的维子空间，其标准正交基为，且不在上。由投影原理知，与其在上的投影之差是正交于子空间的，亦即正交于的正交基。因此只要将单位化，即那么就是在上扩展的子空间的标准正交基。根据上述分析，对于向量组张成的空间 ()，只要从其中一个向量（不妨设为）所张成的一维子空间开始（注意到就是的正交基），重复上述扩展构造正交基的过程，就能够得到 的一组正交基。这就是Gram-Schmidt正交化。具体操作是：q(j)即是QR分解后的Q的列向量，r(i,j)为R的相应元素。r(1,1)*q(1)=v（1）=a(1)r(2,2)*q(2)=v(2)=a(2)-r(1,2)*q(1)r(j,j)*q(j)=v(j)=a(j)-r(1,j)*q(1)-r(j-1,j)*q(j-1)r(n,n)*q(n)=v(n)=a(n)-r(1,n)*q(1)-r(n-1,n)*q(n-1)取r(j,j)=|v(j)|2具体计算时，由于q(i)正交于q(j),i!=j,第j个等式两边用q(i)做内积，得到r(i,j)=q(i)*a(j)= q(i)*v(j)+ q(i)*r(k,j)*q(k)k=1,2,3,j-1。相应算法如下：for j=1:nv(j)=a(j)for i=1:j-1r(i,j)=q(i)a(j)v(j)=v(j)-r(i,j)*q(j)endr(j,j)=|v(j)|2q(j)=v(j)/r(j,j)end4. modified Gram-schmidt正交化方法：在Gram-schmidt正交化过程中产生了r(i,j)=q(i)*v(j)+ q(i)*r(k,j)*q(k)k=1,2,3,j-1后面这一项在计算过程中由于机器精度有限会产生计算误差，故可以优化为r(I,j)=q(i)*v1(j)v1(j)=v(j)+ r(k,j)*q(k)k=1,2,3,j-1具体算法如下：for j=1:nv(j)=a(j)endfor i=1:nr(i,i)=|v(i)|2q(i)=v(i)/r(i,i)for k=i+1:nr(i,k)=q(i)*v(k)v(k)=v(k)-r(i,k)*q(k)end forend for五数值实验结果：为了保证在一定阶数下，所做的数值试验具有代表性，我们随机产生9个n阶矩阵，随着n的不断变化，得到如下图（*）一系列图：6阶矩阵7阶矩阵8阶矩阵9阶矩阵10阶矩阵15阶矩阵16阶矩阵17阶矩阵20阶矩阵30阶矩阵40阶矩阵50阶矩阵60阶矩阵70阶矩阵80阶矩阵90阶矩阵100阶矩阵 图（*）从左到右依次是householder、gramschmidt、modified-gramschmidt、givens变换QR分解。第一个是9次实验相对误差图，第二个是相对时间图。下面是几种古典算法（householder、gramschmidt、modified-gramschmidt、givens）中用时最少的householder变换在矩阵阶数从100增加到1000时所用绝对时间：矩阵阶数1002003004005006007008009001000时间（s）0.04540.34971.25344.700310.709220.458937.432167.953999.9676176.8824 图（一）图（二）由上图可见所用时间随矩阵阶数增长速度略小于指数增长，下用三次样条差值拟合： 图（三）选取图（*）中对应的givens变换相对于householder变换所用相对时间的平均值绘制如下图，得到随着矩阵阶数的增加，givens变换所用时间是householder变换所用时间倍数的关系图： 图（四）六总结与分析：由上面的数值试验结果和绘图可知：1.gramschmidt、modified-gramschmidt正交化方法相对于householder方法的QR分解结果的误差比较小，是householder变换的QR分解误差的0.5倍以下，并且随着计算矩阵阶数的增加，其相对于householder变换的QR分解误差越来越小，矩阵阶数大于30阶是相对误差甚至下降到0.2倍以下。而且在矩阵阶数一定的条件小对不同矩阵相对稳定。所用时间是householder变换的QR分解所用时间的13倍。modified-gramschmidt正交化方法比gramschmidt正交化方法稍微误差小一点，时间也更少一点，但是差距不大。2.Givens变换相对于householder方法的QR分解结果的相对时间在矩阵阶数一定的条件下对于不同矩阵来说变化比较大，比较不稳定。随着矩阵阶数的增加，其相对于householder变换的QR分解的误差逐渐趋于稳定，该稳定值大约是0.8倍。矩阵阶数为10时，在半数以上的情况下误差比householder变化的QR分解误差小，但是在一些情况下比householder变化的QR分解误差大，但是差距不大，其误差是householder变换的QR分解的误差的0.71.5倍之间。所用时间一般情况下比householder变化要长，而且随着矩阵阶数的增加，其所用时间比较快速的增加（接近于线性增加，图（四）所示如），到100阶时已经飙升至5060倍。3.在矩阵阶数比较小（小于17阶）的情况下，givens变换所用时间比gramschmidt、modified-gramschmidt正交化方法短。当矩阵阶数较大时所用时间最长。4.随着矩阵阶数的增加（大于100阶后），所用时间最少的householder变换的QR分解所用时间迅速增长，在1001000阶内其增长速度略小于指数增长（见图（二）。当矩阵阶数是1000阶时householder变换的QR分解所用时间已将近3分钟。可见这几种QR分解方法在矩阵阶数在几百阶内时所用时间还可以接受，矩阵阶数很大时，如果只考虑时间这一因素，上述几种算法都不太理想。4.由以上分析可知在不同工程要求中可以采用不同的QR分解方法，如果工程比较注重于速度，后者先迅速估计一下计算解，可以采用householder变换的QR分解；当矩阵计算规模较小（小于17阶，情况较少）同时要求一定的计算精度时可以考虑用givens变换；当计算规模较大（但是矩阵阶数在几百阶之内时），同时要求高精度时可以采用gramschmidt、modified-gramschmidt正交化方法实现QR分解。七参考文献： 1数值线性代数（徐树方） 2 数值代数（张凯院）附：调用各种方法的函数：for j=10:10:100hold offfor i=1:9a=chooseA(j);q4,r4,l4=givensQR(a);q2,r2,l2=gramschmitQR(a);q1,r1,l1=householderQR(a);q3,r3,l3=MDgramschmitQR(a);v(1)=norm(q1*r1-a,fro);v(2)=norm(q2*r2-a,fro);v(3)=norm(q3*r3-a,fro);v(4)=norm(q4*r4-a,fro);v=v/v(1);subplot(2,1,1)hold onplot(v)plot(v,ob)text(1,v(1),householderQR)text(2,v(2),gramschmitQR)text(4,v(4),givensQR)text(3,v(3),MDgramschmitQR)m=v(1),v(1),v(1),v(1);plot(m,.)ylabel(各种算法相对于householder变化产生的相对误差)l=l1/l1,l2/l1,l3/l1,l4/l1;subplot(2,1,2)hold onplot(l)plot(l,ob)text(1,l(1),householderQR)text(2,l(2),gramschmitQR)text(4,l(4),givensQR)text(3,l(3),MDgramschmitQR)n=l(1),l(1),l(1),l(1);plot(n,.)ylabel(各种算法相对于householder变化所用相对时间)endpauseend