数值分析知识内容 (7).pdf

3.3 矩阵的三角分解 高斯消去法有很多变形，有的是高斯消去法的改进、改写，有的是用于某一类特殊性质矩阵的高斯消去法的简化.本节介绍矩阵的三角分解(Matrix Factorization).3.3.1 直接三角分解法 一、高斯消去法的矩阵解释 由定理2可知，设nnRA，如果),2,1(0)(nkakkk，则可以对A实施高斯消元运算化为上三角形，即存在初等下三角阵121,nLLL，使UALLLn121，U为上三角矩阵.高斯消元过程矩阵描述如下.第一步消元：相当于对)1(A左乘矩阵1L，即)2()1(1AAL，其中)1(11)1(11)2()2(2)2(2)2(22)1(1)1(12)1(11)2(131211,1001011aalaaaaaaaAlllLiinnnnnn.第二步消元：对应于)3()2(2AAL.一般地，1,2,1,)1()(nkAALkkk，其中 nkkiaalllLkkkkikiknkkkk,2,1,10111)()(,1.整个消元过程为 UAALLLLnnn)(1221nnnnuuuuuu22211211，从而 ULULLLLULLLLAnnnn1112121111221)(，其中L是单位下三角矩阵，即 1,2,1,3,2,1111)()(21323121ijniaallllllLjjjjijijnn.【注】消元过程等价于 A 分解为 LU 的过程，回代过程是解上三角方程组的过程.二、矩阵的三角分解（1）杜利特尔分解（Doolittle Factorization）：若将 A 分解成 L U，即 A=L U，其中 L 为单位下三角矩阵（Lower Triangular Matrix）,U 为非奇异上三角矩阵（Upper Triangular Matrix）,则称之为对 A 的 Doolittle 分解.当 A 的顺序主子式都不为零时，消元运算可进行，从而 A 存在唯一的 Doolittle 分解.即，若存在两种分解 A=L1U1和 A=L2U2，则必有 L1=L2，U1=U2.证明 因为 L1U1=L2U2，而且 L1，L2都是单位下三角矩阵，U1,U2都是可逆上三角矩阵,所以有 112112UULL，因此)(112112单位矩阵IUULL，即 L1=L2，U1=U2.（2）克劳特分解（Crout Factorization）：若 L 是非奇异下三角矩阵，U 是单位上三角矩阵时，A 存在唯一的三角分解，A=LU，称其为 A 的 Crout 分解（对应于用列变换实施消元）.三、直接分解（LU分解）算法 对于 Doolittle 分解，可通过直接用 A 元素计算矩阵 A 的三角分解矩阵 L 和 U.这种直接计算 A 的三角分解的方法有实用上的好处.下面利用矩阵乘法规则来确定三角矩阵 L 和U.nnnnnnuuuuuulllllA22211211213231211111.第一步：利用 A 的第一行、第一列元素确定 U 的第一行、L 的第一列元素.由矩阵乘法，),2,1()0,0,()0,0,0,1(1211njuuuuajTijjjj，),3,2()0,0,()0,1,(111111211niululllaiTiiiii，得到),2,1(11njaujj,),3,2(/1111niualii.（3.7）设已经计算出 U 的第 1 至 r-1 行元素，L 的第 1 至 r-1 列元素，现在要计算 U 的第 r行元素及 L 的第 r 列元素.第 r 步：利用 A 的第 r 行、第 r 列剩下的元素确定 U 的第 r 行、L 的第 r 列元素),3,2(nr.由矩阵乘法，有 Tjjjjrrrrrjuuullla)0,0,()0,0,1,(21121),1,(11nrrjuulrjkjrkrk，得 U 的第 r 行元素为 11),3,2,1,(rkkjrkrjrjnrnrrjulau.（3.8）由Trrrriiiiiruuullla)0,0,()0,0,1,(21121),2,1(11nrriululrrirkrrkik，得),1,1,3,2(/)(11nrinruulalrrkrrkikirir.（3.9）公式（3.7）、（3.8）、（3.9）就是 LU 分解的一般计算公式，其结果与高斯消去法所得结果完全一样，但它却避免了中间过程的计算，所以称为A的直接分解公式.为便于记忆，将 LU 的元素写在一起，形成紧凑格式.由于 L 和 U 中的元素计算好后，A 中对应元素就不用了，因此计算好 L，U 的元素后就存放在 A 的相应位置.nnnnnnnnnnnnnnnnullluulluuuluuuuaaaaaaaaaaaaaaaaA321333323122322211131211321333323122322211131211，最后在存放 A 的矩阵中得到 L，U 的元素.在 LU 分解计算公式中，需要计算形如iiba的式子，可采用“双精度累加”，以提高精度.四、方程组的三角分解算法（LU 分解算法解方程组）对于方程组bAx，设由 Doolittle 分解得到 A=LU，则方程组bAx 的求解转化为两个三角形方程组bLy 和yUx 的求解.解方程组的 LU 分解算法步骤如下.（1）分解 L 和 U 的过程：对于nr,2,1，11),1,(rkkjrkrjrjnrrjulau，),1(/)(11nriuulalrrkrrkikirir；（2）求解单位下三角方程组bLy，得 1111),3,2(,ikkikiiniylbyby；（3）求解上三角方程组yUx，得)1,2,2,1(,/)(,/1nniuxuyxuyxiinikkikiinnnn.LU 分解算法大约需要3/3n次乘除法，与高斯消去法的计算量基本相同.例 5 用 LU 分解法求解方程组 713542774322322xxx.解 对系数矩阵 A 进行 LU 分解，1,2,3,2,23121131211lluuu.由jjjulau12122，有1,32322uu.2/)(2212313232uulal，6)(233213313333ululau.因此 613322121121LUA.解方程组bLy，得627,521,3213121yyyyyy.解方程组yUx，得22/)323(,23/)5(,1321323xxxxxx.解方程组的LU分解法MATLAB程序：LU_Factorization.mfunction x,L,U,index=LU_Factorization(A,b)%LU分解法解方程组，其中，%A为要分解的矩阵；%b为方程组的右端常数项；%x为方程组的解；%L为单位下三角阵；%U为上三角阵；%index为指标变量，index=0表示计算失败，index=1表示计算成功.n,m=size(A);nb=length(b);%当方程组行与列的维数不相等时，停止计算，并输出出错信息.if n=m error(The rows and columns of matrix A must be equal!);return;end%当方程组与右端项的维数不匹配时，停止计算，并输出出错信息.if m=nb error(The columns of A must be equal the length of b!);return;end%开始计算，先赋初值 L=eye(n);U=zeros(n);index=1;x=zeros(n,1);y=zeros(n,1);%矩阵三角分解过程 for k=1:n for j=k:n z=0;for q=1:k-1 z=z+L(k,q)*U(q,j);end U(k,j)=A(k,j)-z;end if abs(U(k,k)1e-10 index=0;return;end for i=k+1:n z=0;for q=1:k-1 z=z+L(i,q)*U(q,k);end L(i,k)=(A(i,k)-z)/U(k,k);end end%求解两个三角形方程组的过程 y(1)=b(1);for k=2:n z=0;for j=1:k-1 z=z+L(k,j)*y(j);end y(k)=b(k)-z;end x(n)=y(n)/U(n,n);for k=n-1:-1:1 z=0;for j=k+1:n z=z+U(k,j)*x(j);end x(k)=(y(k)-z)/U(k,k);end 调用函数 LU_Factorization.m 解例 5.输入 A=2 2 3;4 7 7;-2 4 5;b=3;1;-7;x,L,U,index=LU_Factorization(A,b)得到方程组的解及相应的 LU 分解矩阵：x=2 L=1 0 0 U=2 2 3 index=1 计算成功.-2 2 1 0 0 3 1 1 -1 2 1 0 0 6 3.3.2 列主元三角分解法 当用LU分解法解方程组时，由第r),2,1(nr步分解计算公式可知 11),1,(rkkjrkrjrjnrrjulau，),1(/)(11nriuulalrrkrrkikirir.当0rru时，计算将中断；或者当rru很小时，可能引起舍入误差的累积、扩大.但如果 A 非奇异，我们可通过交换 A 的行实现矩阵 PA 的 LU 分解.因此，可采用与列主元消去法类似方法，将直接三角分解法修改为列主元三角分解法（与列主元消去法在理论上是等价的），它通过交换 A 的行实现三角分解 PA=LU，其中 P 为初等置换阵.设第 1 步至 r-1 步分解计算己完成，则有 nnnrrnnrnrrrrnrrrnaallaaluuuluuA1,11,11,12221111.第 r),2,1(nr步计算：为了避免用绝对值很小的数作除数，引进中间量：11),(rkkrikirinriulaS，则有：),1(,nriSSlSuriirrrr.（1）选列主元，即确定inirirSSir max,使.（2）交换两行：当rir时，交换 A 的第 r 行与第ri行元素（相当于先交换原始矩阵A 第 r 行与第ri行元素后，再进行分解计算得到的结果，且1irl）.（3）进行第 r 步分解计算.经过带行交换的列主元 LU 分解后，矩阵 PA=LU，则解方程组bAx 即为求解方程组PbPAx，转化为解两个三角方程组yUxPbLy,.例 6 利用列主元 LU 分解法解方程组.11163,1852,174321321321xxxxxxxxx 解 对系数矩阵作列主元 LU 分解得到 PA=LU，其中 12/13/213/11L，13/1021163U，010001100P.解PbLy，得Ty)0,3/2,1(.解yUx，得Tx)0,3/1,3/1(.LU 分解的 MATLAB 函数 MATLAB 提供了列主元三角分解函数 LU()，用此函数求解例 5.输入 A=2 2 3;4 7 7;-2 4 5;b=3;1;-7;L,U=lu(a)y=Lb;x=Uy 得到方程组的解及相应的列主元 LU 分解矩阵：L=0.5000 -0.2000 1.0000 U=4.0000 7.0000 7.0000 x=2.0000 1.0000 0 0 0 7.5000 8.5000 -2.0000 -0.5000 1.0000 0 0 0 1.2000 1.0000 或者输入 A=2 2 3;4 7 7;-2 4 5;b=3;1;-7;L,U,P=lu(A)y=L(P*b);x=Uy 得到列主元 LU 分解矩阵、所用的排列矩阵及方程组的解：L=1.0000 0 0 -0.5000 1.0000 0 0.5000 -0.2000 1.0000 U=4.0000 7.0000 7.0000 0 7.5000 8.5000 0 0 1.2000 P=0 1 0 0 0 1 1 0 0 x=2.0000 -2.0000 1.0000 3.3.3 平方根法 在工程技术问题中，例如用有限元方法解结构力学中问题时，常常需要求解具有对称正定矩阵的方程组，对于这种具有特殊性质系数的矩阵，利用矩阵的三角分解法求解就得到解对称正定矩阵方程组的平方根法.平方根法是解对称正定矩阵方程组的有效方法，目前在计算机上被广泛应用.一、平方根法（乔莱斯基方法）设n阶方程组bAx，A是对称正定矩阵(Positive Definite Matrix)，则A有三角分解 nnnnnnuuuuuulllllLUA22211211213231211111.再将U分解为 nnnnuuuuuuU22211211022211111122211111UDuuuuuuuuunnnn,则0UDLA.（1）对称正定矩阵有唯一的分解TLDLA 这是由于0UDLA，TTTDLUA0，且对称阵AAT，则有 TTDLULDU00 再利用三角分解的唯一性，得LUT0.因此，对称正定矩阵有唯一的分解TLDLA.（2）D是正定对角阵（即niuii,2,1,0）由于A对称正定的充要条件是TBAB对称正定，其中B是n阶可逆方阵.取TLB，就推知D是正定对角阵.因此D的对角元素),2,1(0niuii，记2121DDD，其中 nnuuuD221121,则TTLDLDLDLDA)()(21212121.（3）乔莱斯基（Cholesky）分解 将21LD记为L，则TLLA 称为 Cholesky 分解.利用 Cholesky 直接分解公式，推导出的解方程组方法，称为 Cholesky 方法或平方根法.（4）解方程组的平方根法（Cholesky 方法）由 Cholesky 分解，有 TnnnnnnnnLLllllllllllllA2221211121222111.（3.10）利用矩阵乘法，逐步确定L的第i行元素),2,1(ijlij.由11jkjjijjkikijlllla（当kj 时，0jkl），有分解公式：对于ni,2,1 112/1211)(,)1,2,1(/)(ikikiiiijkjjjkikijijlalijlllal （3.11）将对称正定矩阵A作 Cholesky 分解后，则解方程组bAx 就转化为解两个三角方程组yxLbLyT,.【注】（1）平方根算法是一个数值稳定的方法.这是因为，由分解公式（11）可得，ikikiinila12),2,1(，于是，),2,1;,2,1(,max12ikniaaliiniiiik.这说明在不选主元的平方根法中，矩阵L元素有界，或者说在分解过程中产生的L元素ikl(即消元乘数)的数量级不会增长，且),2,1(0nilii.（2）平方根分解算法的计算量约为 6/3n次乘除法，大约为LU分解计算量的一半；存贮量也大约为LU分解法的一半.例 7 用 Cholesky 方法解方程组 25.7645.375.2175.225.41114x.解 对系数矩阵作 Cholesky 分解得到 15.125.05.0215.15.025.025.375.2175.225.41114TLL.解bLy，得Ty)1,5.3,2(.解yxLT，得Tx)1,1,1(.Cholesky 分解的 MATLAB 函数 MATLAB 提供了 Cholesky 分解函数 chol()，用此函数求解例 7.输入 A=4-1 1;-1 4.25 2.75;1 2.75 3.5;b=4;6;7.25;L=chol(A),y=Lb;x=Ly 得到 Cholesky 分解矩阵（上三角矩阵）及方程组的解：L=2.0000 -0.5000 0.5000 x=1 0 2.0000 1.5000 1 0 0 1.0000 1 对称正定矩阵的Cholesky分解算法 Andre-Louis Cholesky(1875-1918)是一位法国军官，他在第一次世界大战之前，参加了在克里特岛和北非的大地测量和调查.为了求解大地测量中出现的最小二乘数据拟合问题的法方程，他提出了后来以他的名字命名的方法.在他死后的1924年，一位叫Benoit的军官以 Cholesky的名义在Bulletin Geodesique（测地学快报）上发表了这项工作.二、改进的平方根法 由公式（3.11）可知，用平方根法解对称正定方程组时，计算L的元素iil需要用到开方运算.为避开开方运算，我们将平方根法改进得到改进的平方根法，该算法既适合A对称正定，也适合A对称且顺序主子式全不为零的情况.（1）TLDL分解算法 因为对称正定矩阵有唯一的分解TLDLA，即 11111111232131212121323121nnnnnllllldddlllllA.由矩阵乘法 nknkjkjijjkkikjkjjjijjkkikjkkikkjTiKijdlldlldlldlldlLLDa111111)()(，所以求L的元素和D的元素公式为：111ad，1,2,1,3,2/)(11ijnidldlaljjkjkkikijij，),3,2(112nidladikkikiii.（2）简化的TLDL分解算法 为避免重复计算，令jijijdlt，则得到简化的TLDL分解算法：111ad 对于ni,3,2，1111,1,2,1,/,ikikikiiijijijjkjkikijijltadijdtlltat （3.12）（3）改进的平方根法（方程组的TLDL分解算法）将对称正定矩阵A作TLDL分解后，则解方程组bAx 就转化为解两个三角方程组bLy 和yDxLT1，从而得到求解方程组的算法公式.先解bLy，即),3,2(,1111niylbybyikkikii.再解yDxLT1，即)1,2,1(/,/1nnixldyxdyxnikkkiiiinnn.例 8 解方程组 25.15.065.375.2175.225.41114x.解 由TLDL分解算法，得 175.025.0125.01L，144D.解bLy，得Ty)11,6(.解yDxLT1，得Tx)1,1,2(.【注】（1）改进的平方根分解算法的计算量约为 6/3n次乘除法，但没有开方运算.（2）平方根法或改进的平方根法是目前计算机上解对称正定线性方程组的一个有效方法，比用消去法优越.其计算量和存贮量都比用消去法大约节省一半左右，且数值稳定、不需要选主元，能求得较高精度的计算解.解方程组的改进平方根法MATLAB程序：LDL_Factorization.mfunction x,L,D,index=LDL_Factorization(A,b)%改进平方根法解方程组，其中，%A为要分解的矩阵；%b为方程组的右端常数项；%x为方程组的解；%L为下三角阵；%D为对角阵；%index为指标变量，index=0表示计算失败，index=1表示计算成功.n,m=size(A);nb=length(b);%当方程组行与列的维数不相等时，停止计算，并输出出错信息.if n=m error(The rows and columns of matrix A must be equal!);return;end%当方程组与右端项的维数不匹配时，停止计算，并输出出错信息.if m=nb error(The columns of A must be equal the length of b!);return;end%开始计算，先赋初值L=eye(n);D=zeros(n);d=zeros(1,n);T=zeros(n);index=1;x=zeros(n,1);y=zeros(n,1);%矩阵LDL的分解过程 d(1)=A(1,1);for i=2:n for j=1:i-1 T(i,j)=A(i,j);for k=1:j-1 T(i,j)=T(i,j)-T(i,k)*L(j,k);end L(i,j)=T(i,j)/d(j);end d(i)=A(i,i);for k=1:i-1 d(i)=d(i)-T(i,k)*L(i,k);end if abs(d(i)1e-10 index=0;return;end end D=diag(d);%求解两个三角形方程组的过程 y(1)=b(1);for i=2:n y(i)=b(i);for k=1:i-1 y(i)=y(i)-L(i,k)*y(k);end end x(n)=y(n)/d(n);for i=n-1:-1:1 x(i)=y(i)/d(i);for k=i+1:n x(i)=x(i)-L(k,i)*x(k);end end 调用函数 LDL_Factorization.m 解例 8.输入 A=4-1 1;-1 4.25 2.75;1 2.75 3.5;b=6-0.5 1.25;x,L,D,index=LDL_Factorization(A,b)得到方程组的解及相应的 LDLT分解矩阵：x=2 L=1.0000 0 0 D=4 0 0 index=1 计算成功.1 -0.2500 1.0000 0 0 4 0 -1 0.2500 0.7500 1.0000 0 0 1 3.4.4 追赶法 在许多实际问题中，如，常微分方程两点边值问题、三次样条插值方法等，往往遇到线性方程组fAx 的求解，其中 TnnnnnnfffbacbacbacbA),(,111122211.（3.13）称具有公式（3.13）形式的系数矩阵A为三对角阵,称相应的线性方程组为三对角方程组（Tridiagonal Linear Systems).具有这种形式的方程组在实际问题中是经常遇到的，而且往往是对角占优（Diagonally Dominant）的.A满足条件：011 cb，)1,3,2(0,0nicacabiiiii，0nnab.这类方程组fAx 的解存在唯一（A非奇异），可以直接利用高斯消去法或直接分解法，而其解答可以用极其简单的递推公式表示出来，即下面介绍的追赶法.追赶法通常是数值稳定的.一、分解算法 对A作 LU 分解（Doolitle 分解），可以发现 L、U 具有非常简单的形式.nnnluululmmmULA12211321111.由矩阵乘积，得 nnnnnnnnnnnlumlmuulumlmulbacbacbacbA1112212121111122211.比较等式两端，得到.,3,2 ,3,2,/,1,2,1,1111niumblnilamblnicuiiiiiiiii （3.14）二、解方程组的追赶法 因为上述分解LUA，则方程组fAx 的求解转化为解两个简单的三角方程组fLy 和yUx，从而得到求解方程组的算法公式.先解fLy，即),3,2(,111niymfyfyiiii.（3.15）再解yUx，即)1,2,1(/)(,/1nnilxuyxlyxiiiiinnn.（3.16）这种把三对角方程组的解用递推公式（3.14）、（3.15）、（3.16）表示出来的方法形象化地叫做追赶法，其中（3.14）、（3.15）是关于下标i由小到大的递推公式称为追的过程，而（16）却是下标i由大到小的递推公式称为赶的过程，一追一赶构成了求解fAx 的追赶法.【注】追赶法的优点：存贮单元少，计算量小，且舍入误差不增长，算法数值稳定.追赶法有条件0iica，如果有某0ia或0ic，则三对角方程组fAx 可化为两个低阶方程组，例如当03a时，43214321443322211ffffxxxxbacbcbacbAx.于是，解fAx 化为解两个方程组 43434433ffxxbacb和3221212211xcffxxbacb.例 9 用追赶法解三对角方程组 2314114114x.解 系数矩阵分解得到 7333.3175.314,12667.0125.01UL.解fLy，得Ty)8667.2,25.3,1(.解yUx，得Tx)7679.0,0714.1,5179.0(.调用函数 LU_Factorization.m 解例 9.输入 A=4-1 0;-1 4-1;0-1 4;b=1;3;2;x,L,U,index=LU_Factorization(A,b)得到方程组的解及相应的 LU 分解矩阵：x=0.5179 L=1.0000 0 0 U=4.0000 -1.0000 0 1.0714 -0.2500 1.0000 0 0 3.7500 -1.0000 0.7679 0 -0.2667 1.0000 0 0 3.7333