武汉大学数学与统计学院 数值实验报告

1、武汉大学数学与统计学院数值分析实验报告实验名称关于正定矩阵cholesky分解的研究实验时间2006年 12月 16 日姓名雷锦江 班级05级数类2班学号200531000174成绩一、实验目的，内容 二、相关背景知识介绍 三、代码四、数值结果 五、计算结果的分析 六、计算中出现的问题，解决方法及体会一、 实验目的，内容通过对正定矩阵cholesky分解问题的深入研究，掌握对矩阵进行cholesky分解的不同方法。同时对C语言的编程方法进行练习，在此基础上，学会使用数学软件matlab对相关问题进行处理。二、 相关背景知识介绍定理(Cholesky分解)如果A是对称正定的，则存在惟一的一个对角

2、元全部大于零的下三角 阵，满足。证明:由对称矩阵的分解定理,存在单位下三角阵L和对角阵D=diag使得A=LD.由于大于零，则矩阵是对角元大于零的实下三角阵.它同时满足.惟一性由分解的惟一性可推得。分解被称为Cholesky分解,G被称为Cholesky三角阵.如果我们计算Cholesky分解,然后解三角形方程组Gy=b和x=y,则b=Gy=G(x)=Gx=Ax.在上述定理中,Cholesky分解的证明是构造性的.然而,可以通过利用方程来得到计算Cholesky三角阵的更有效的方法.例：矩阵是正定的.常规Cholesky分解对于n阶正定矩阵A，可导出按列的次序来计算其Cholesky分解因子G

3、的元素的计算公式。对第k=1,2,，n列，i=k+1,k+2,n 基于Gaxpy的Cholesky分解我们首先导出一个含有大量Gaxpy运算的Cholesky分解的实现方法.比较等式的第j列可得A(:,j)=.也就是说,G(j,j)(:,j)=A(:,j)-v.(*)如果G的第j-1列已知,则可计算出v.由(*)中各元素间的相等关系推出G(j:n,j)=v(j:n)/.这是数乘的Gaxpy运算。于是得到基于运算Gaxpy的Cholesky分解计算方法:For j =1:nv(j:n)=A(j:n,j)For k=1:j-1 v(j:n)=v(j:n)-G(j,k)G(j:n,k)EndG(j:

4、n,j)=v(j:n)/ End可以在计算过程中用G覆盖A的下三角部分.算法(Cholesky分解:基于Gaxpy运算)给出对称正定阵,本算法计算出一个下三角阵满足.对所有ij,G(I,j)覆盖A(i,j).For j=1:nIf j1A(j:n,j)=A(j:n,j)-A(j:n,1:j-1)EndA(j:n,j)=A(j:n,j)/End本算法需flops.基于外积的Cholesky分解另一种基于外积(秩为1)的Cholesky分解可通过对矩阵做如下划分得到: (*)这里,且因是A正定的,故.而是的主子阵,其中故它也是正定的.如果存在Cholesky分解= ,则由(4.2.6)有,其中.因

5、此可通过反复利用(*)来得到最终的Cholesky分解,其方向和kji形式的高斯消去法一样.算法(Cholesky分解:基于外积运算)给定一对称正定阵.本算法计算满足的下三角阵,对所有的ij,G(i,j)覆盖A(i,j)For k=1:nA(k,k)=A(k+1:n,k)=A(k+1:n,k)/A(k,k)For j=k+1:nA(j:n,j)=A(j:n,j)-A(j:n,k)A(j,k)End End本算法需flops.注意,其中的j循环计算外积的下三角部分:A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k).回想在Gaxpy运算和外积运算的比较,容易得知Ga

6、xpy运算中读写向量的次数要比外积运算中少一半.基于分块点积的Cholesky分解假定对称正定,将A=()和其Cholesky因子看作含有方块对角块的N*N分块阵.取出等式中关于第(I,j)块(ij)的等式有=定义 S=-则有,当i=j时,当ij时, .通过合适的排序,这些方程可用来求得所有的算法(Cholesky:基于分块点积运算)给定一个对称正定阵,本算法可求得一个下三角阵满足,A的下三角部分被G覆盖,A被看作是含方对角块的N*N分块阵.For j=1:NFor i=j:NS=-If i=j计算Cholesky分解S=Else从=S解出End用覆盖End End整个算法需flops,与前述

7、其他形式的Cholesky算法相同.假定A被适当分块,则本算法中含有大量的矩阵乘法运算.例如,如果n=rN,且每个都是R*R的,则3级运算量约占1-(1/.).由于没有给出基于gaxpy运算的Cholesky算法得出。将基于gaxpy的Cholesky分解执行r步后，我们已知中的和。再接下来我们不对A，而是对已明显形成具有可利用的对称结构的约化后的矩阵进行另外的r步基于gaxpy的Cholesky计算。按此方法进行下去，我们得到基于分块的Cholesky计算。按此方法进行下去，我们得到基于分块的Cholesky分解算法，去第k不是对n-(k-1)r阶的矩阵进行r次基于gaxpy的Cholesk

8、y分解，接着是阶为n-kr的三级算法。如果，则三级运算量约占1-3/（2N）。Cholesky分解的稳定性在精确运算下,我们知对称正定阵存在Cholesky分解.反之,如果一个Cholesky分解过程能够顺利进行完且得到严格大于零的平方根,那么A是正定的.因此,判断是否正定,我们只需用上述任一方法来试着计算其Cholesky分解.有舍入误差是情形更为有意思. Cholesky算法的数值稳定性大致可从不等式:导出.该不等式说明Cholesky三角阵因子有很好的界.由也可推出相同的结论.Wilkinson(1968)在其经典的论文中对Cholesky分解的舍入误差作了深入研究.利用该论文的结果可证

9、明,如果是通过上述任一Cholesky分解求得的Ax=b的计算解,则满足扰动方程组(A+E)=b,其中=,是由决定的小常数. Wilkinson进一步指出,如果,其中是另一个小常数,则分解能够执行到底,而不会出现对负数开平方根.三、代码常规Cholesky分解#define N 3#includemath.hmain()float aNN,gNN,*pa,temp1,temp2; int i,j,k; pa=a0; printf(Please input the maxim:n); for(i=0;iN;i+) pa=pa+i; for(j=i;jN;j+)scanf(%f,pa+); for

10、(i=1;iN;i+) for(j=0;ji;j+) aij=aji; for(k=0;kN;k+) temp1=0; for(j=0;jk;j+) temp1=temp1+gkj*gkj; gkk=sqrt(akk-temp1); for(i=k+1;iN;i+) temp2=0; for(j=0;jk;j+) temp2=temp2+gij*gkj; gik=(aik-temp2)/gkk; for(i=0;iN-1;i+) for(j=i+1;jN;j+) gij=0; for(i=0;iN;i+) printf(n); for(j=0;jN;j+)printf(%6.2f ,gij);

11、 getch(); 基于Gaxpy的Cholesky分解#define N 5#includemath.hmain()float aNN,*pa,temp; int i,j,k; pa=a0; printf(Please input the maxim:n); for(i=0;iN;i+) pa=pa+i; for(j=i;jN;j+)scanf(%f,pa+); for(i=1;iN;i+) for(j=0;ji;j+) aij=aji; for(j=0;jN;j+) temp=0; for(k=0;kj;k+)temp=temp+ajk*ajk; ajj=sqrt(ajj-temp); f

12、or(i=j+1;iN;i+)temp=0; for(k=0;kj;k+) temp=temp+ajk*aik; aij=(aij-temp)/ajj; for(i=0;iN-1;i+) for(j=i+1;jN;j+) aij=0; for(i=0;iN;i+) printf(n); for(j=0;jN;j+)printf(%6.2f ,aij); getch(); 基于外积的Cholesky分解#define N 4#includemath.hmain()float aNN,*pa; int i,j,p,q; pa=a0; printf(Please input the maxim:n)

13、; for(i=0;iN;i+) pa=pa+i; for(j=i;jN;j+)scanf(%f,pa+); for(i=1;iN;i+) for(j=0;ji;j+) aij=aji; for(j=0;jN;j+) for(p=j+1;pN;p+)for(q=j+1;qN;q+) apq=apq-apj*aqj/ajj; ajj=sqrt(ajj); for(i=j+1;iN;i+)aij=aij/ajj; for(i=0;iN-1;i+) for(j=i+1;jN;j+) aij=0; for(i=0;iN;i+) printf(n); for(j=0;jN;j+)printf(%6.2f

14、 ,aij); getch(); 基于分块点积的Cholesky分解#define N 5#includemath.hmain()float aNN,*pa,temp; int i,j,k,p,q; pa=a0; printf(Please input the maxim:n); for(i=0;iN;i+) pa=pa+i; for(j=i;jN;j+)scanf(%f,pa+); for(i=1;iN;i+) for(j=0;ji;j+) aij=aji; p=N/3; for(q=0;qp;q+) for(j=3*q;j3*q+3;j+)temp=0; for(k=0;kj;k+) te

15、mp=temp+ajk*ajk; ajj=sqrt(ajj-temp); for(i=j+1;iN;i+) temp=0; for(k=0;kj;k+)temp=temp+ajk*aik; aij=(aij-temp)/ajj; if(N-p*3)=2) for(j=N-2;jN;j+) temp=0; for(k=0;kj;k+)temp=temp+ajk*ajk; ajj=sqrt(ajj-temp); for(i=j+1;iN;i+)temp=0; for(k=0;kj;k+) temp=temp+ajk*aik; aij=(aij-temp)/ajj; if(N-p*3)=1) tem

16、p=0; for(k=0;kN-1;k+) temp=temp+aN-1k*aN-1k; aN-1N-1=sqrt(aN-1N-1-temp); for(i=0;iN-1;i+) for(j=i+1;jN;j+) aij=0; for(i=0;iN;i+) printf(n); for(j=0;jN;j+)printf(%6.2f ,aij); getch();四数值结果常规Cholesky分解Please input the maxim:4 -2 4 10 -2 8 2.00 0.00 0.00 -1.00 3.00 0.00 2.00 0.00 2.00基于Gaxpy的Cholesky分解

17、Please input the maxim:1 1 1 1 1 2 2 2 2 3 3 3 4 4 1.00 0.00 0.00 0.00 0.00 1.00 1.00 0.00 0.00 0.00 1.00 1.00 1.00 0.00 0.001.00 1.00 1.00 1.00 0.001.00 1.00 1.00 1.00 1.00基于外积的Cholesky分解Please input the maxim:1 1 1 1 2 2 2 3 3 4 1.00 0.00 0.00 0.00 1.00 1.00 0.00 0.00 1.00 1.00 1.00 0.00 1.00 1.00

18、 1.00 1.00 基于分块点积的Cholesky分解Please input the maxim9 6 3 12 3 29 12 13 12 41 30 47 37 40 64300 0.00 0.00 0.00 0.00 2.00 5.00 0.00 0.00 0.00 1.00 2.00 6.00 0.00 0.004.00 1.00 4.00 2.00 0.001.00 2.00 7.00 3.00 1.00五计算结果分析由于在编程中才有的是浮点数，所以由于系统的限制，该程序的精确度最高能达到6-7位有效数字，而且由于水平有限，该程序只是针对正定矩阵进行Cholesky分解的。对于标

19、准的正定矩阵，且在有效数字位数之内，该程序是有效的。在程序的输出中，为了使输出更为明了清晰，我选择了保留两位小数，再不影响精度的情况下，这种选取是可取的。六计算中出现的问题，解决方法及体会这篇论文从开始策划到最终完成，共花了我们一个多星期的时间，从最开始的担心无法完成课题的研究，到后来的查资料，编程，到最后的完成论文，我们体会到了其间的种种苦乐，当最终论文完成时，我们感到就如同完成了一件花了很多时间和精力完成的工艺品，也许有很多地方略显粗糙，但这毕竟是我们努力的结果，我们对于这样的结果是满意的，已经没有遗憾了。在进行Cholesky分解的研究过程中，我认为最关键的是对于算法的研究和程序的编写，

20、由于相关资料略显简略，我花了一定的时间将算法弄透，接下来就进入了程序的编写过程。在编程中，我确实遇到了一些困难，其中多次出现了参数设定的错误，而且该错误不能被低阶的矩阵所反映，所以我曾一度以为程序编写完成，而实际上程序中还存在着问题，当我对程序作最后测试的时候，才发现了问题，并花了相当长的时间完成了程序的矫正。从中我也得到了颇多的收获。同时我想指出的是，在考虑基于分块矩阵的Cholesky分解时，我曾一度想让计算机能够根据矩阵的阶数自动将其分解，其具体方法是设定一个素数库，然后在素数库中按从小到大的顺序，依次除以矩阵阶数，直到出现能整除的素数为止，其数学基础是应用了代数基本原则。但后来我又考虑

21、到如果矩阵的阶数是一个很大的素数，那么这种方法将转变为基于Gaxpy的Cholesky分解，而且由于素数库的限制，该算法不能将某些矩阵分解。于是我最终采取了前面提到的关于基于分块矩阵的Cholesky分解。最后，要感谢向老师对我们这个课题研究的关心，他在百忙之中抽空对我们的研究提出了宝贵的意见，在这里向他表示诚挚的感谢！七参考书目：1.MATRIX COMPUTATIONS（Third Edition） Gene H. Golub & Charles F. Van Loan2.数值计算方法 郑慧娆 陈绍林 莫忠息 黄象鼎（与邢周华同学合作完成）教 师 评 语指导教师： 年 月 日武汉大学数学与统计学院数值分析原始书记记录实验名称：关于正定矩阵cholesky分解的研究 实验时间