武汉大学数学与统计学院_数值分析实验报告.doc

武汉大学数学与统计学院 数值分析实验报告 实验名称 关于正定矩阵cholesky分解的研究 实验时间 2006年 12月 16 日 姓名 雷锦江 班级 05级数类2班 学号 200531000174 成绩 一、实验目的，内容 二、相关背景知识介绍 三、代码 四、数值结果 五、计算结果的分析 六、计算中出现的问题，解决方法及体会 实验目的，内容 通过对正定矩阵cholesky分解问题的深入研究，掌握对矩阵进行cholesky分解的不同方法。同时对C语言的编程方法进行练习，在此基础上，学会使用数学软件matlab对相关问题进行处理。 相关背景知识介绍 定理(Cholesky分解)如果A是对称正定的，则存在惟一的一个对角元全部大于零的下三角 阵，满足。 证明:由对称矩阵的分解定理,存在单位下三角阵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的元素的计算公式。对第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:n v(j:n)=A(j:n,j) For k=1:j-1 v(j:n)=v(j:n)-G(j,k)G(j:n,k) End G(j:n,j)=v(j:n)/ End 可以在计算过程中用G覆盖A的下三角部分. 算法(Cholesky分解:基于Gaxpy运算)给出对称正定阵,本算法计算出一个下三角阵满足.对所有i≥j,G(I,j)覆盖A(i,j). For j=1:n If j1 A(j:n,j)=A(j:n,j)-A(j:n,1:j-1) End A(j:n,j)=A(j:n,j)/ End 本算法需flops. 基于外积的Cholesky分解 另一种基于外积(秩为1)的Cholesky分解可通过对矩阵做如下划分得到: (**) 这里,且因是A正定的,故.而是的主子阵,其中 故它也是正定的.如果存在Cholesky分解= ,则由(4.2.6)有,其中 . 因此可通过反复利用(**)来得到最终的Cholesky分解,其方向和kji形式的高斯消去法一样. 算法(Cholesky分解:基于外积运算)给定一对称正定阵.本算法计算满足的下三角阵,对所有的ij,G(i,j)覆盖A(i,j) For k=1:n A(k,k)= A(k+1:n,k)=A(k+1:n,k)/A(k,k) For j=k+1:n A(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运算和外积运算的比较,容易得知Gaxpy运算中读写向量的次数要比外积运算中少一半. 基于分块点积的Cholesky分解 假定对称正定,将A=()和其Cholesky因子看作含有方块对角块的N*N分块阵.取出等式中关于第(I,j)块(i≥j)的等式有 = 定义 S=- 则有,当i=j时,当ij时, .通过合适的排序,这些方程可用来求得所有的 算法(Cholesky:基于分块点积运算)给定一个对称正定阵,本算法可求得一个下三角阵满足,A的下三角部分被G覆盖,A被看作是含方对角块的N*N分块阵. For j=1:N For i=j:N S=- If i=j 计算Cholesky分解S= Else 从=S解出 End 用覆盖 End End 整个算法需flops,与前述其他形式的Cholesky算法相同.假定A被适当分块,