QR分解数值效果分析
报告
班级：0811101 学号：081110215 姓名：桑树浩
1
目 录
§1. 摘要 ...................................................................................................... 3
§2. 关键词 .................................................................................................. 3
§3. 引言 ...................................................................................................... 3
§4. 理论基础.............................................................................................. 4
§5. 数值实验.............................................................................................. 8
§6. 总结与分析 ....................................................................................... 19
§7.参考文献与附注…………………………………………………….20
2
QR分解数值效果分析报告
基于数值试验的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||Q’Ax-Q’b||2,为方便求解，需要Q’A是上三角矩阵，这样引入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变换，它有如下形式：
3
其中旋转变换。
3.Gram-Schmidt正交化方法：在线性代数中，如果内积空间上的一组向量能够组成一个子空间，那么这一组向量就称为这个子空间的一个基。Gram－Schmidt正交化提供了一种方法，能够通过这一子空间上的一个基得出子空间的一个正交基，并可进一步求出对应的标准正交基。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。
.易证
是一个正交阵。Givens变换又叫平面
如上图中所示，记v=Hx-x w=v/||v||2 H=I-2*w*w’
4
上述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;
end q=q^-1;
r=a; %a即为变化之后的矩阵，即R 2.givens变换：
任意的一个向量x被其作用后，y=G*x,有 yi=c*xi+s*xk yk=-s*xi+c*xk
yj=xj, j~=i,k。
将一个矩阵QR分解时，依次将该矩阵的第一列的第2个、第3个…….第n
5
个元素变为0，第二列第3、4…..n个元素变为0……第n-1列元素的第n个元素变为0. 具体算法如下：
q=diag(ones(n,1));%初始化Q for 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 end end r=a; q=q^-1;
3．Gram-Schmidt正交化基本想法，是利用投影原理在已有正交基的基础上
构造一个新的正交基。设且不在
。
是
上的维子空间，其标准正交基为
上的投影
，之差
上。由投影原理知，与其在
6
是正交于子空间即
的，亦即
正交于
的正交基
。因此只要将
单位化，
那么
就是
的标准正交基。
根据上述分析，对于向量组要从其中一个向量（不妨设为到到
就是
张成的空间
）所张成的一维子空间
(
)，只
开始（注意
在上扩展的子空间
的正交基），重复上述扩展构造正交基的过程，就能够得
的一组正交基。这就是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:n v(j)=a(j) for i=1:j-1
r(i,j)=q(i)’a(j) v(j)=v(j)-r(i,j)*q(j) end
r(j,j)=||v(j)||2 q(j)=v(j)/r(j,j) end
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]后面这一项在
7
各种算法相对于householder变化产生的相对误差各种算法相对于householder变化所用相对时间计算过程中由于机器精度有限会产生计算误差，故可以优化为
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:n v(j)=a(j) end
for i=1:n
r(i,i)=||v(i)||2 q(i)=v(i)/r(i,i) for k=i+1:n
r(i,k)=q(i)’*v(k) v(k)=v(k)-r(i,k)*q(k) end for end for
五．数值实验结果：
为了保证在一定阶数下，所做的数值试验具有代表性，我们随机产生9个n阶矩阵，随着n的不断变化，得到如下图（*）一系列图： 6阶矩阵
1.5givensQRgivensQRgivensQRgivensQRgivensQRgivensQRgivensQRgramschmitQRgramschmitQRgramschmitQRgramschmitQRgramschmitQR11.522.5MDgramschmitQRMDgramschmitQRMDgramschmitQRMDgramschmitQRMDgramschmitQR33.541householderQR0.50432gramschmitQR1householderQR0MDgramschmitQRgivensQRgivensQR4gramschmitQRMDgramschmitQR11.522.533.57阶矩阵
8
于householder变化产生的相于householder变化产生的相对误差对间对间相时相时法对法相对于householder变化法对算相算差算相种用种误种用各种算法相对于householder变化各所各对各所1.5givensQRgivensQRgivensQRgivensQRgivensQRgivensQRgramschmitQRgramschmitQRgramschmitQRgramschmitQRgramschmitQR11.522.5MDgramschmitQRMDgramschmitQRMDgramschmitQR33.541householderQR0.502gramschmitQRgramschmitQRgramschmitQR1.5MDgramschmitQRMDgramschmitQRgivensQRgivensQRgivensQRgivensQR1householderQR0.511.522.533.548阶矩阵
1.21householderQR0.80.60.40.211.5gramschmitQRgramschmitQRgramschmitQRgramschmitQRgramschmitQR22.5MDgramschmitQRMDgramschmitQRMDgramschmitQRMDgramschmitQRMDgramschmitQR33.5
givensQRgivensQRgivensQRgivensQRgivensQRgivensQRgivensQR421.81.61.41.21householderQR11.522.533.5gramschmitQRgramschmitQRMDgramschmitQRMDgramschmitQRMDgramschmitQRgivensQRgivensQR49阶矩阵
9
差误于householder变化产生的相对于householder变化产生的相对误差对间对间相时相时法对法对算相算法相对于householder变化算相种用种种用各种算法相对于householder变化各所各各所1.21householderQR0.80.60.40.211.5gramschmitQRgramschmitQRgramschmitQRgramschmitQRgramschmitQR22.5MDgramschmitQRMDgramschmitQRMDgramschmitQRMDgramschmitQRMDgramschmitQRMDgramschmitQR33.5givensQRgivensQRgivensQRgivensQRgivensQR421.81.61.41.21householderQR11.522.533.5gramschmitQRgramschmitQRgramschmitQRMDgramschmitQRMDgramschmitQRMDgramschmitQRMDgramschmitQRgivensQRgivensQRgivensQRgivensQRgivensQR410阶矩阵
1.5
1householderQRgivensQRgivensQRgivensQRgivensQRgivensQRgivensQRgivensQRgivensQRgramschmitQRgramschmitQRgramschmitQRgramschmitQRMDgramschmitQRMDgramschmitQRMDgramschmitQRMDgramschmitQR33.540.5011.522.532.521.51householderQR11.5gramschmitQRgramschmitQRgramschmitQRMDgramschmitQRMDgramschmitQR33.5givensQRgivensQRgivensQRgivensQR4gramschmitQRMDgramschmitQR22.515阶矩阵
10
于householder变化产生的相于householder变化产生的相对误差对间对间相时相时法对法相对于householder变化法对算相算差算相种用种误种用各种算法相对于householder变化各所各对各所1.5givensQRgivensQRgivensQRgivensQRgivensQRgivensQRgivensQRgramschmitQRgramschmitQRgramschmitQRgramschmitQRgramschmitQR11.522.5MDgramschmitQRMDgramschmitQRMDgramschmitQR33.541householderQR0.501.81.61.41.21householderQR11.5gramschmitQRgramschmitQRMDgramschmitQRMDgramschmitQRMDgramschmitQRgivensQRgivensQR22.533.5416阶矩阵
1.5
givensQR1householderQRgivensQRgivensQRgivensQRgivensQRgivensQRgramschmitQRgramschmitQR011.522.5MDgramschmitQRMDgramschmitQR33.540.52.52gramschmitQRgramschmitQRgramschmitQRgramschmitQRgramschmitQRMDgramschmitQRMDgramschmitQRMDgramschmitQRgivensQRgivensQRgivensQR1.51householderQR11.522.533.5417阶矩阵
11
于householder变化产生的相于householder变化产生的相对误差对间对间相时相时法对法相对于householder变化法对算相算差算相种用种误种用各种算法相对于householder变化各所各对各所1householderQR0.80.60.40.2011.5gramschmitQRgramschmitQRgramschmitQR22.5MDgramschmitQRMDgramschmitQRMDgramschmitQR33.5givensQRgivensQRgivensQRgivensQR432.521.51householderQR11.5gramschmitQRgramschmitQRMDgramschmitQRMDgramschmitQRgivensQRgivensQRgivensQR22.533.5420阶矩阵
1householderQR0.80.60.40.2011.5gramschmitQRgramschmitQRgramschmitQR22.5MDgramschmitQRMDgramschmitQRMDgramschmitQRMDgramschmitQR33.54
givensQRgivensQRgivensQRgivensQR3.532.521.51householderQR11.522.533.54gramschmitQRgramschmitQRgramschmitQRgramschmitQRgramschmitQRMDgramschmitQRMDgramschmitQRMDgramschmitQRMDgramschmitQRMDgramschmitQRgivensQRgivensQRgivensQRgivensQRgivensQRgivensQRgivensQR30阶矩阵
12
于householder变化产生的相对误差对间相时法对算相种用各所各种算法相对于householder变化产生的相对误差各种算法相对于householder变化各种算法相对于householder变化所用相对时间1householderQRgivensQRgivensQRgivensQRgivensQR0.80.60.40.2011.5gramschmitQR22.5MDgramschmitQRMDgramschmitQR33.548642householderQR011.522.533.54gramschmitQRgramschmitQRMDgramschmitQRMDgramschmitQRgivensQRgivensQRgivensQRgivensQR40阶矩阵
1householderQR
0.80.60.40.2011.5gramschmitQR22.5MDgramschmitQR33.5givensQRgivensQRgivensQRgivensQR415givensQRgivensQRgivensQR105gramschmitQRgramschmitQR0householderQR11.522.533.54MDgramschmitQR50阶矩阵
13
于householder变化产生的相于householder变化产生的相对误差对间对间相时相时法对法对算相算相种用种用各所各所法相对于householder变化算差种误各种算法相对于householder变化各对1householderQR0.80.60.40.2011.5gramschmitQRgramschmitQR22.5MDgramschmitQR33.54givensQRgivensQRgivensQRgivensQRgivensQR20151050householderQR11.5gramschmitQRgramschmitQR22.5givensQRgivensQRgivensQRgivensQRgivensQRMDgramschmitQRMDgramschmitQRMDgramschmitQRMDgramschmitQR33.5460阶矩阵
1householderQR0.80.60.40.2011.5gramschmitQR22.5MDgramschmitQRMDgramschmitQR33.54
givensQRgivensQRgivensQRgivensQRgivensQR2520151050householderQR11.5gramschmitQRgramschmitQRgramschmitQR22.5MDgramschmitQRMDgramschmitQRMDgramschmitQR33.5givensQRgivensQRgivensQRgivensQRgivensQRgivensQRgivensQR470阶矩阵
14
于householder变化产生的相于householder变化产生的相对误差对间对间相时相时法对法对算相算相种用种用各所各所法相对于householder变化算差种误各种算法相对于householder变化各对1householderQR0.80.60.40.2011.5gramschmitQR22.5MDgramschmitQR33.54givensQRgivensQRgivensQR3020givensQRgivensQRgivensQRgivensQRgivensQRgivensQR100householderQR11.5gramschmitQRgramschmitQRgramschmitQR22.5MDgramschmitQRMDgramschmitQR33.5480阶矩阵
1householderQR0.80.60.40.2011.5gramschmitQR22.5MDgramschmitQR33.54
givensQRgivensQRgivensQR403020100householderQR11.5gramschmitQRgramschmitQR22.5MDgramschmitQRMDgramschmitQR33.5givensQRgivensQRgivensQRgivensQRgivensQR490阶矩阵
15
于householder变化产生的相于householder变化产生的相对误差对间对间相时相时法对法对算相算相种用种用各所各所法相对于householder变化算差种误各种算法相对于householder变化各对1householderQR0.80.60.40.2011.5gramschmitQR22.5MDgramschmitQR33.54givensQRgivensQRgivensQR50403020100householderQR11.5gramschmitQRgramschmitQR22.5MDgramschmitQRMDgramschmitQR33.54givensQRgivensQRgivensQRgivensQRgivensQR100阶矩阵
1householderQR0.80.60.40.2011.5gramschmitQR22.5MDgramschmitQR33.54
givensQRgivensQRgivensQR60givensQRgivensQRgivensQRgivensQRgivensQRgivensQRgivensQRgivensQRgivensQR40200householderQR11.5gramschmitQR22.5MDgramschmitQR33.54 图（*）
从左到右依次是householder、gramschmidt、modified-gramschmidt、givens变换QR分解。第一个是9次实验相对误差图，第二个是相对时间图。
16
下面是几种古典算法（householder、gramschmidt、modified-gramschmidt、givens）中用时最少的householder变换在矩阵阶数从100增加到1000时所用绝对时间： 矩阵阶 100 数 200 300 400 500 600 700 800 900 1000 时间0.0454 （s） 0.3497 1.2534 4.7003 10.7092 20.4589 37.4321 67.9539 99.9676 176.8824 180160140120所用时间100806040200100200300400500600矩阵阶数7008009001000
图（一）
17
6543所用时间取对数210-1-2-3-4100200300400500600矩阵阶数7008009001000
图（二）
由上图可见所用时间随矩阵阶数增长速度略小于指数增长，下用三次样条差值拟合：
180160140120所用时间100806040200100200300400500600矩阵阶数7008009001000
18
图（三）
选取图（*）中对应的givens变换相对于householder变换所用相对时间的平均值绘制如下图，得到随着矩阵阶数的增加，givens变换所用时间是householder变换所用时间倍数的关系图：
40givens变换是householder变换所用时间的倍数353025201510500102030405060矩阵阶数708090100
图（四） 六．总结与分析：
由上面的数值试验结果和绘图可知：
1.gramschmidt、modified-gramschmidt正交化方法相对于householder方法的QR分解结果的误差比较小，是householder变换的QR分解误差的0.5倍以下，并且随着计算矩阵阶数的增加，其相对于householder变换的QR分解误差越来越小，矩阵阶数大于30阶是相对误差甚至下降到0.2倍以下。而且在矩阵阶数一定的条件小对不同矩阵相对稳定。所用时间是householder变换的QR分解所用时间的1~3倍。modified-gramschmidt正交化方法比gramschmidt正交化方法稍微误差小一点，时间也更少一点，但是差距不大。
2.Givens变换相对于householder方法的QR分解结果的相对时间在矩阵阶数一定的条件下对于不同矩阵来说变化比较大，比较不稳定。随着矩阵阶数的增加，其相对于householder变换的QR分解的误差逐渐趋于稳定，该稳定值大约是0.8倍。矩阵阶数为10时，在半数以上的情况下误差比householder变化的QR分解误差小，但是在一些情况下比householder变化的QR分解误差大，但是差距不大，其误差是householder变换的QR分解的误差的0.7~1.5倍之间。所用时间一般情况下比householder变化要长，而且随着矩阵阶数的增加，其所用时间比较快速的增加（接近于线性增加，图（四）所示如），到100阶时已经飙升至50~60倍。
19
3.在矩阵阶数比较小（小于17阶）的情况下，givens变换所用时间比gramschmidt、modified-gramschmidt正交化方法短。当矩阵阶数较大时所用时间最长。 4.随着矩阵阶数的增加（大于100阶后），所用时间最少的householder变换的QR分解所用时间迅速增长，在100~1000阶内其增长速度略小于指数增长（见图（二））。当矩阵阶数是1000阶时householder变换的QR分解所用时间已将近3分钟。可见这几种QR分解方法在矩阵阶数在几百阶内时所用时间还可以接受，矩阵阶数很大时，如果只考虑时间这一因素，上述几种算法都不太理想。
4.由以上分析可知在不同工程要求中可以采用不同的QR分解方法，如果工程比较注重于速度，后者先迅速估计一下计算解，可以采用householder变换的QR分解；当矩阵计算规模较小（小于17阶，情况较少）同时要求一定的计算精度时可以考虑用givens变换；当计算规模较大（但是矩阵阶数在几百阶之内时），同时要求高精度时可以采用gramschmidt、modified-gramschmidt正交化方法实现QR分解。
七．参考文献：
[1]数值线性代数（徐树方）
[2] 数值代数（张凯院） 附：调用各种方法的函数：
for j=10:10:100 hold off for i=1:9
a=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 on plot(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];
20
subplot(2,1,2) hold on plot(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变化所用相对时间') end pause end
21
因篇幅问题不能全部显示，请点此查看更多更全内容