1、2 三角分解法三角分解法 /*Matrix Factorization*/高斯消元法的矩阵形式高斯消元法的矩阵形式 /*Matrix Form of G.E.*/：Step 1:)0(/111111 aaamii记记 L1=1.11121nmm ，则，则)1()1(1bAL)1(1)1(1)1(11.baan)2(A)2(bStep n 1:)()2(2)1(1)()2(2)2(22)1(1)1(12)1(11121.nnnnnnnnnbbbaaaaaabALLL其中其中 Lk=1.11,1knkkmm 2 Matrix Factorization Matrix Form of G.E.1kL
2、1.11,1knkkmm 111211.nLLL111jim,记为记为L单位下三角阵单位下三角阵/*unitary lower-triangular matrix*/记记 U=)()2(2)2(22)1(1)1(12)1(11.nnnnnaaaaaaLUA A 的的 LU 分解分解/*LU factorization*/Hey hasnt GE given me enough headache?Why do I have to know its matrix form?!When you have to solve the system for different with a fixed A
3、.bCould you be more specific,please?Factorize A first,then for every you only have to solve two simple triangular systems and.bbyL yxU 2 Matrix Factorization Matrix Form of G.E.证明证明：由由1 1中定理可知，中定理可知，LU LU 分解存在。事实上分解存在。事实上,kkkkALALLLA)(111211.)()()1(1)()()1(11knnkknnknkkkkkaaaaaaA其中其中 形如形如kA knkkkIHM
4、L0)(1111121kkkkmmmM其中其中 0001,1,1,21,21,11,1knnkkkkkkkmmmmmmH 若若A的所有顺序主子式的所有顺序主子式/*determinant of leading principal submatrices*/均不为均不为0，则，则 A 的的 LU 分解唯一（其分解唯一（其中中 L 为为单位单位下三角阵）。下三角阵）。定理定理记为分块形式记为分块形式 )()()()(22211211222112110kkkkknkkAAAAIHMAAAA)(1111kkAMA 显然显然)()2(22)1(11)(11)(1111)det()det()det()de
5、t(kkkkkkaaaAAMA 从而从而因此因此 即即A有有LU分解的充要条件是分解的充要条件是A的所有顺序主子式非奇异。的所有顺序主子式非奇异。,1,2,1,0)(nkakkk2 Matrix Factorization Matrix Form of G.E.其中其中 为为 的的k k阶顺序主子阵，阶顺序主子阵，为为 的的k k阶顺序主子阵阶顺序主子阵。)(11kA11AAkA2 Matrix Factorization Matrix Form of G.E.下面证明唯一性。下面证明唯一性。若不唯一，则可设若不唯一，则可设 A=L1U1=L2U2，推出，推出 121UU211122211LL
6、UULL Upper-triangularLower-triangularWith diagonal entries 1I L 为一般下三角阵而为一般下三角阵而 U 为为单位单位上三角阵的分解称为上三角阵的分解称为Crout 分解分解。实际上只要考虑实际上只要考虑 A*的的 LU 分解，即分解，即 ，则，则 即是即是 A 的的 Crout 分解。分解。ULA*LUA 2 Matrix Factorization Doolittle 道立特分解法道立特分解法 /*Doolittle Factorization*/：LU 分解的紧凑格式分解的紧凑格式 /*compact form*/反复计算反复计
7、算,很浪费哦很浪费哦 通过比较法直接导出通过比较法直接导出L 和和 U 的计算公式。的计算公式。思思路路 nnnnnnnnuuullaaaa.1.11.1111211111 ),min(1jikjkkiul jia2 Matrix Factorization Doolittle ),min(1jikjkkijiula固定固定 i:对对 j=i,i+1,n 有有ijkjikikijuula 11lii=1kjikikijijulau 11a将将 i，j 对换，对对换，对 j=i,i+1,n 有有iijikiikjkjiulula 11iiikkijkjijiuulal/)(11 b Algori
8、thm:Doolittle FactorizationStep 1:u1j=a1j;lj1=aj1/u11;(j=1,n)Step 2:compute and for i=2,n 1;Step 3:ab 11nkknnknnnnulau一般采用一般采用列主元列主元法增强稳定性。但注意法增强稳定性。但注意 也必须做相应的也必须做相应的行交换。行交换。b2 Matrix Factorization Cholesky 平方根法平方根法 /*Choleskys Method*/：对称对称/*symmetric*/正定正定/*positive definite*/矩阵的分解法矩阵的分解法回顾：回顾：对称
9、正定阵的几个重要性质对称正定阵的几个重要性质 A 1 亦对称正定，且亦对称正定，且 aii 0若不然，则若不然，则0 xA存在非零解，即存在非零解，即0 xAxT存在非零解。存在非零解。1111)()(,AAIAAIAATTT对任意对任意 ,存在存在 ,使得使得 ,即即 。0 x0 yxyA xAy1 011 yAyyAAAyxAxTTT 其中其中0 xAxaTiiTx)0.1.0(第第 i 位位 A 的顺序主子阵的顺序主子阵/*leading principal submatrices*/Ak 亦对亦对 称正定称正定对称性显然。对任意对称性显然。对任意 有有 ,其中其中 。kkRx 00 x
10、AxxAxTkkTknkRxx 00 A 的特征值的特征值/*eigen value*/i 0 设对应特征值设对应特征值 的非零特征向量的非零特征向量为为 ，则，则 。20 xxxxAxTT x A 的全部顺序主子式都有的全部顺序主子式都有 det(Ak)0因为因为 niiA1)det(定义定义一个矩阵一个矩阵 A=(aij)n n 称为称为对称阵对称阵，如果，如果 aij=aji。定义定义一个矩阵一个矩阵 A 称为称为正定阵正定阵，如果，如果 对任意非对任意非零向量零向量 都成立都成立。0 xAxTx2 Matrix Factorization Choleski将将对称对称 正定阵正定阵 A
11、 做做 LU 分解分解U=uij=u11uij/uii111u22unn记为记为UD A 对称对称TUL 即即TLDLA 记记 D1/2=11u22unnuWhy is uii 0?Since det(Ak)02/1LDL 则则 仍是下三角阵仍是下三角阵TLLA nnRL 定理定理 设矩阵设矩阵A对称正定，则存在非奇异下三角阵对称正定，则存在非奇异下三角阵 使得使得 。若限定。若限定 L 对角元为正，则分解唯一。对角元为正，则分解唯一。TLLA 对于对称正定阵对于对称正定阵 A，从，从 可知对任意可知对任意k i 有有 。即即 L 的元素不会增大，误差可控，不的元素不会增大，误差可控，不需选主
12、元。需选主元。ikikiila12iiikal|2 Matrix Factorization Cholesky Algorithm:Choleskys MethodTo factor the symmetric positive definite n n matrix A into LLT,where L is lower triangular.Input:the dimension n;entries aij for 1 i,j n of A.Output:the entries lij for 1 j i and 1 i n of L.Step 1 Set ;Step 2 For j=2,
13、n,set ;Step 3 For i=2,n 1,do steps 4 and 5Step 4 Set ;Step 5 For j=i+1,n,set ;Step 6 Set ;Step 7 Output(lij for j=1,i and i=1,n);STOP.1111al 1111/laljj 112ikikiiiilal iiikikjkjijilllal 11 112nknknnnnlal因为因为A 对称，所以只需存半个对称，所以只需存半个 A，即，即其中其中 nnnaaaaannA.,.,2/)1(1222111 jiiAaij 2/)1(运算量为运算量为 O(n3/6)，比普通
14、比普通LU分解少一半，但有分解少一半，但有 n 次开方。用次开方。用A=LDLT 分解，可省开方时间分解，可省开方时间(p.50-51)。2 Matrix Factorization Tridiagonal System 追赶法解追赶法解三对角三对角方程组方程组 /*Crout Reduction for Tridiagonal Linear System*/nnnnnnnfffxxxbacbacbacb212111122211Step 1:对对 A 作作LU分解，比如（分解，比如（Crout 分解、分解、Doolittle分解）分解）111121nnnA 直接比较等式两边直接比较等式两边的元
15、素，可得到计的元素，可得到计算公式算公式（p.182）。Step 2:追追即解即解 ：fyL,111 fy ),.,2()(1niyrfyiiiii Step 3:赶赶即解即解 ：yxU)1,.,1(,1 nixyxyxiiiinn 与与G.E.类似，一旦类似，一旦 i=0 则算法中断，故并非任何则算法中断，故并非任何三对角阵都可以用此方法分解。三对角阵都可以用此方法分解。2 Matrix Factorization Tridiagonal System定理定理 若若 A 为为对角占优对角占优/*diagonally dominant*/的三对角的三对角阵，且满足阵，且满足 ，则追赶，则追赶法
16、可解以法可解以 A 为系数矩阵的方程组。为系数矩阵的方程组。0,0,0|,0|11 iinncaabcb如果如果 A 是是严格严格对角占优阵，则不要求三对角线上的对角占优阵，则不要求三对角线上的所有元素非零。所有元素非零。根据不等式根据不等式 可知：分解过程中，矩阵元素不会过分增大，算法可知：分解过程中，矩阵元素不会过分增大，算法保证稳定。保证稳定。运算量为运算量为 O(6n)。|,1|1iiiiiiiiabbab 2 Matrix Factorization Tridiagonal SystemLab 06.Crout Reduction for Tridiagonal Linear Sys
17、temsApply Crout Reduction to solve a given nn tridiagonal linear systemInputThere are several sets of inputs.For each set:The 1st line contains an integer 100 n 0 which is the size of a matrix.n=1 signals the end of file.The 2nd line contains n 1 real numbers .The 3rd line contains n real numbers .T
18、he 4th line contains n 1 real numbers .The 5th line contains n real numbers .The numbers are separated by spaces.naaa.32 nnnnnnnfffxxxbacbacbacb212111122211nbbb.21121.ncccnfff.212 Matrix Factorization Tridiagonal SystemOutputEach entry of the solution is to be printed as in the C fprintf:fprintf(out
19、file,%16.8en,x);If the method fails to give a solution,print the message “TheCroutmethodfailed.n”./*here represents a space*/The outputs of two test cases must be seperated by a blank line.Sample InputSample Input (represents a space)represents a space)511114444411111002002002001003110330.518871Sample Output Sample Output(represents a space)represents a space)4.61538462e+0018.46153846e+0019.23076923e+0018.46153846e+0014.61538462e+001 The Crout method failed.
