《矩阵计算》课程教学资源（讲义）第二讲 线性方程组直接方法

:46第二讲线性方程组直接方法end for11:12:end for评价算法的一个主要指标是执行时间，但这依赖于计算机硬件和编程技巧等，因此直接给出算法执行时间是不太现实的.所以我们通常是统计算法中算术运算（加减乘除）的次数在矩阵计算中，大多仅仅涉及加减乘除和开方运算，一般情况下，加减运算次数与乘法运算次数具有相同的量级，而除法运算和开方运算次数具有更低的量级凸为了尽可能地减少运算量，在实际计算中，数，向量和矩阵做乘法运算时的先后执行次序为：先计算数与向量的乘法，然后计算矩阵与向量的乘法，最后才计算矩阵与矩阵的乘法，比如计算αAB&，其中α是数，A,B是矩阵，是向量，如果按照从左往右计算的话，则运算量为O(n3)，但是如果先计算αa,然后计算B(αa)，最后再计算A(B(αa)的话，运算量则为O(n2),相差一个量级.矩阵L和U的存储当A的第i列（严格下三角部分）被用于计算L的第列后，在后面的计算中不再被使用而A的第行（上三角部分）更新后就是U的第i行.因此，为了节省存储空间，我们可以在计算过程中将L的第i列存放在A的第i列（严格下三角部分，L的对角线全部为1,不需要存储)，将U的第i行存放在A的第i行（上三角部分），这样就不需要另外分配空间存储L和U.计算结束后，A的上三角部分为U，其严格下三角部分为L的绝对下三角部分.此时算法可以描述为：算法2.3. LU分解(用 A存储L和U)1: for k = 1 to n - 1 do2:fori=k+ltondo3:aik=aik/akk4:for j= k + 1 to n do5:ai=aij-aakj6:end for7:end for8:end forLU分解的运算量由算法2.2可知，LU分解的运算量为1-m-)=(2n3 -3n2+n)=gn3+0(n2).·乘法次数：Tp=k-1k=1 i=k+1 j=k+1

· 46 · 第二讲 线性方程组直接方法 11: end for 12: end for b 评价算法的一个主要指标是执行时间, 但这依赖于计算机硬件和编程技巧等, 因此直接给 出算法执行时间是不太现实的. 所以我们通常是统计算法中算术运算 (加减乘除) 的次数. 在矩阵计算中, 大多仅仅涉及加减乘除和开方运算. 一般情况下, 加减运算次数与乘法运算 次数具有相同的量级, 而除法运算和开方运算次数具有更低的量级. b 为了尽可能地减少运算量, 在实际计算中, 数, 向量和矩阵做乘法运算时的先后执行次序为: 先计算数与向量的乘法, 然后计算矩阵与向量的乘法, 最后才计算矩阵与矩阵的乘法. 比如 计算 αABx, 其中 α 是数, A, B 是矩阵, x 是向量, 如果按照从左往右计算的话, 则运算量 为 O(n 3 ), 但是如果先计算 αx, 然后计算 B(αx), 最后再计算 A(B(αx)) 的话, 运算量则为 O(n 2 ), 相差一个量级. 矩阵 L 和 U 的存储 当 A 的第 i 列 (严格下三角部分) 被用于计算 L 的第 i 列后, 在后面的计算中不再被使用. 而 A 的第 i 行 (上三角部分) 更新后就是 U 的第 i 行. 因此, 为了节省存储空间, 我们可以在计算过程 中将 L 的第 i 列存放在 A 的第 i 列 (严格下三角部分, L 的对角线全部为 1, 不需要存储), 将 U 的 第 i 行存放在 A 的第 i 行 (上三角部分), 这样就不需要另外分配空间存储 L 和 U. 计算结束后, A 的上三角部分为 U, 其严格下三角部分为 L 的绝对下三角部分. 此时算法可以描述为： 算法 2.3. LU 分解 (用 A 存储 L 和 U) 1: for k = 1 to n − 1 do 2: for i = k + 1 to n do 3: aik = aik/akk 4: for j = k + 1 to n do 5: aij = aij − aikakj 6: end for 7: end for 8: end for LU 分解的运算量 由算法 2.2 可知, LU 分解的运算量为 • 乘法次数: Tp = nX−1 k=1 Xn i=k+1 Xn j=k+1 1 = nX−1 k=1 (n − k) 2 = 1 6 (2n 3 − 3n 2 + n) = 1 3 n 3 + O(n 2 )

2.1LU分解与Gauss消去法:47.n-1TE111132Nn3 +0(n2)7·加法次数：Ta=2n3=1i=k+1j=k+1121-11112gn(n-1) =+ 0(n)·除法次数：Ta=-nSn-2k=1i=k+121122-3n3+ O(n2).因此总运算量（含乘法、除法与加法）为-n6n=33数据更新顺序与计算效率根据指标的循环次序，算法2.3也称为KII型LU分解，在实际计算中，我们一般不建议使用这个算法.因为对于指标k的每次循环，都需要更新A的第k十1至第n行.这种反复读取数据的做法会使得计算效率大大降低如果数据是按行存储的，如C/C++，为了提高计算效率，我们可以采用下面的IKI型LU分解算法.图IKI型LU分解算法可以用下图来描述算法2.4.LU分解(IKJ型)([113, p. 291])1: fori = 2 to n do2:for k=l to i-l doAccessed but notmodified3:aik=aik/akkNfor j = k + 1 to n doiccessedand5:aij=aj-aiakimodified6:endforNotaccessedend for7:8: end for思考：如果数据是按列存储的，如FORTRAN或MATLAB，则怎样设计算法比较好？2.1.3Gauss消去法假定矩阵A存在LU分解A=LU，则方程组Ac=b就转化为求解下面两个三角方程组Ly=b, Ua=y.这两个方程组都非常容易求解，于是Gauss消去法描述如下：算法2.5.Gauss消去法1:将A进行LU分解：A=LU,其中L为单位下三角矩阵，U为非奇异上三角矩阵；2:利用向前回代，求解Ly=b,即得y=L-1b;3:利用向后回代，求解U&=y,即得=U-1y=(LU)-1b=A-1b.得到A的LU分解后，我们最后需要用回代法求解两个三角方程组，计算过程描述如下

2.1 LU 分解与 Gauss 消去法 · 47 · • 加法次数: Ta = nX−1 k=1 Xn i=k+1 Xn j=k+1 1 = 1 3 n 3 − 1 2 n 2 + 1 6 n = 1 3 n 3 + O(n 2 ). • 除法次数: Td = nX−1 k=1 Xn i=k+1 1 = 1 2 n(n − 1) = 1 2 n 2 − 1 2 n = 1 2 n 2 + O(n). 因此总运算量 (含乘法、除法与加法) 为 2 3 n 3 − 1 2 n 2 − 1 6 n = 2 3 n 3 + O(n 2 ). 数据更新顺序与计算效率 根据指标的循环次序, 算法 2.3 也称为 KIJ 型 LU 分解. 在实际计算中, 我们一般不建议使用 这个算法. 因为对于指标 k 的每次循环, 都需要更新 A 的第 k + 1 至第 n 行. 这种反复读取数据的 做法会使得计算效率大大降低. 如果数据是按行存储的, 如 C/C++, 为了提高计算效率, 我们可以采用下面的 IKJ 型 LU 分解 算法. 算法 2.4. LU 分解 (IKJ 型) 1: for i = 2 to n do 2: for k = 1 to i − 1 do 3: aik = aik/akk 4: for j = k + 1 to n do 5: aij = aij − aikakj 6: end for 7: end for 8: end for  IKJ 型 LU 分解算法可以用下图来描述. ([113, p. 291]) 思考：如果数据是按列存储的, 如 FORTRAN 或 MATLAB, 则怎样设计算法比较好? 2.1.3 Gauss 消去法 假定矩阵 A 存在 LU 分解 A = LU, 则方程组 Ax = b 就转化为求解下面两个三角方程组 Ly = b, Ux = y. 这两个方程组都非常容易求解, 于是 Gauss 消去法描述如下: 算法 2.5. Gauss 消去法 1: 将 A 进行 LU 分解: A = LU, 其中 L 为单位下三角矩阵, U 为非奇异上三角矩阵; 2: 利用向前回代, 求解 Ly = b, 即得 y = L −1 b; 3: 利用向后回代, 求解 Ux = y, 即得 x = U −1y = (LU) −1 b = A−1 b. 得到 A 的 LU 分解后, 我们最后需要用回代法求解两个三角方程组, 计算过程描述如下

-48第二讲线性方程组直接方法算法2.6.回代求解Ly=b和Ua=y1:y1=b1/l11%向前回代求解Ly=b2: fori=2:n do3: for j=l:i-1do4:bi=bi-lijyjend for5:6:yi= bi/lu7:end for%向后回代求解Ur=y8: In = yn/unn9: fori= n-1:-1:1do10:fori=n:-l:i+1do11:=yi-uirjend for12:13:Ti=yi/ui14:end for如果数据是按列存储的，则采用列存储方式效率会高一些。下面是以U=y为例，描述按列存储时的回代求解过程算法2.7.向后回代求解UT=u（列存储方式）1: for k = n:-1: 1 do2:k=yk/ukk3:for i=k-1:-1:1do4:i=yi-rkuik5:end for6:end for计算复杂度这两个算法的运算量：乘法和加法均n(n-1)次，除法2n次.加上LU分解的运算量，Gauss消去法的总运算量为n3+O(n2).可以证明，以上两个算法都是向后稳定的（componentwisebackwardstable)[70]

· 48 · 第二讲 线性方程组直接方法 算法 2.6. 回代求解 Ly = b 和 Ux = y 1: y1 = b1/l11 % 向前回代求解 Ly = b 2: for i = 2 : n do 3: for j = 1 : i − 1 do 4: bi = bi − lijyj 5: end for 6: yi = bi/lii 7: end for 8: xn = yn/unn % 向后回代求解 Ux = y 9: for i = n − 1 : −1 : 1 do 10: for j = n : −1 : i + 1 do 11: yi = yi − uijxj 12: end for 13: xi = yi/uii 14: end for 如果数据是按列存储的, 则采用列存储方式效率会高一些. 下面是以 Ux = y 为例, 描述按列 存储时的回代求解过程. 算法 2.7. 向后回代求解 Ux = y (列存储方式) 1: for k = n : −1 : 1 do 2: xk = yk/ukk 3: for i = k − 1 : −1 : 1 do 4: yi = yi − xkuik 5: end for 6: end for 计算复杂度 这两个算法的运算量: 乘法和加法均 n(n − 1) 次, 除法 2n 次. 加上 LU 分解的运算量, Gauss 消 去法的总运算量为 2 3 n 3 + O(n 2 ). b 可以证明, 以上两个算法都是向后稳定的 (componentwise backward stable) [70]

2.1LU分解与Gauss消去法:49.2.1.4选主元LU分解在LU分解过程中，我们称ak-1）为主元。如果akt-1）=0,则算法就无法进行下去。即使at-1）不为零，但如果ax-1"]的值很小，由于舍人误差的原因，也可能会给计算结果带来很大的dkk误差此时我们就需要通过选主元来解决这个问题，[0.0261.361.5例2.1用LU分解求解线性方程组A=b,其中A：要求在运算过3.43-8.525.8(板书)程中保留3位有效数字解。根据LU分解待定系数法，设u11u12A= LU012100直接比较等式两边可得u11 = a11 = 0.02, u12 = a12 = 61.3l21 = a21/u11 = 3.43/0.02 ~ 172,u22=a22-l21u12~-8.5-1.05×104~-1.05×104解方程组Ly=b可得y1= 61.5,y2= b2-l211~-1.06×104解方程组Ur=y可得2= y2/u22~1.01,1=(y1 u122)/u11~-0.413/u11 ~-20.7口凸事实上，精确解为1=10.0和22=1.00.我们发现1的数值解误差非常大.出现这个问题的原因就是1a111太小，用它做主元时会放大舍人误差选主元是通过利用置换矩阵（也称排列矩阵）对行进行交换来实现的，首先介绍置换矩阵的些基本性质引理2.4设PeRnxn为置换矩阵，AERnxn为任意矩阵，则(1)PA相当于将A的行进行置换；AP相当于将A的列进行置换；(2)P-1 =PT,即 P 是正交矩阵;(3) det(P) = ±1;（4）置换矩阵的乘积仍然是置换矩阵定理2.5（选主元LU分解的存在性）设AERnxn非奇异，则存在置换矩阵PL，PR，以及单位下三角矩阵L和非奇异上三角矩阵U，使得PLAPR=LU

2.1 LU 分解与 Gauss 消去法 · 49 · 2.1.4 选主元 LU 分解 在 LU 分解过程中, 我们称 a (k−1) kk 为主元. 如果 a (k−1) kk = 0, 则算法就无法进行下去. 即使 a (k−1) kk 不为零, 但如果 |a (k−1) kk | 的值很小, 由于舍入误差的原因, 也可能会给计算结果带来很大的 误差. 此时我们就需要通过选主元来解决这个问题. 例 2.1 用 LU 分解求解线性方程组 Ax = b, 其中 A = " 0.02 61.3 3.43 −8.5 # , b = " 61.5 25.8 # , 要求在运算过 程中保留 3 位有效数字. (板书) 解. 根据 LU 分解待定系数法, 设 A = LU = " 1 0 l21 1 # "u11 u12 0 u22# . 直接比较等式两边可得 u11 = a11 = 0.02, u12 = a12 = 61.3, l21 = a21/u11 = 3.43/0.02 ≈ 172, u22 = a22 − l21u12 ≈ −8.5 − 1.05 × 104 ≈ −1.05 × 104 , 解方程组 Ly = b 可得 y1 = 61.5, y2 = b2 − l21y1 ≈ −1.06 × 104 . 解方程组 Ux = y 可得 x2 = y2/u22 ≈ 1.01, x1 = (y1 − u12x2)/u11 ≈ −0.413/u11 ≈ −20.7 □ b 事实上, 精确解为 x1 = 10.0 和 x2 = 1.00. 我们发现 x1 的数值解误差非常大. 出现这个问 题的原因就是 |a11| 太小, 用它做主元时会放大舍入误差. 选主元是通过利用置换矩阵 (也称排列矩阵) 对行进行交换来实现的. 首先介绍置换矩阵的 一些基本性质. 引理 2.4 设 P ∈ R n×n 为置换矩阵, A ∈ R n×n 为任意矩阵, 则 (1) P A 相当于将 A 的行进行置换; AP 相当于将 A 的列进行置换; (2) P −1 = P T, 即 P 是正交矩阵; (3) det(P) = ±1; (4) 置换矩阵的乘积仍然是置换矩阵. 定理 2.5 (选主元 LU 分解的存在性) 设 A ∈ R n×n 非奇异, 则存在置换矩阵 PL, PR, 以及单位下 三角矩阵 L 和非奇异上三角矩阵 U, 使得 PLAPR = LU

.50.第二讲线性方程组直接方法（板书）证明.用归纳法当n=1时,取PL=PR=L=1,U=A即可假设结论对所有n－1阶矩阵都成立，考虑n阶非奇异矩阵AERnxn，易知A至少存在一个非零元，取置换矩阵PI和P2使得[a1 A12]PiAP2 :[A2 A22]其中 a11 ≠ 0, A22 ER(n-1)x(n-1),令U12=A12，L21=A2i/a11,U22=A22-L2iU12u11=a11,则u11≠0,且有A121 U12a11u11PIAP2.L210U22A21A22两边取行列式可得0 + det(PiAP2) = deta11 - det(U22).0Looo所以det(U22)≠0,即 U22 E R(n-1)x(n-1)非奇异.由归纳假设可知,存在置换矩阵PLEIR(n-1)x(n-1)和PRER(n-1)x(n-1),使得P,U22PR = L22U22,其中L22为单位下三角矩阵，U22为非奇异上三角矩阵.令10Pi,PR=P2PL=0 PL0PR则有101PiAPPLAPRPLD012/11:PD010U12u1PLOPRPLL210PTL22U2PT(1010U12U2PToPR[PLL2]PL] [o PTL22]010U12PRu11 LU,L22U22[PLL2]0其中L为单位下三角矩阵，U为非奇异上三角矩阵，即A存在LU分解因此，由数学归纳法可知，结论成立

· 50 · 第二讲 线性方程组直接方法 (板书) 证明. 用归纳法. 当 n = 1 时, 取 PL = PR = L = 1, U = A 即可. 假设结论对所有 n − 1 阶矩阵都成立. 考虑 n 阶非奇异矩阵 A ∈ R n×n , 易知 A 至少存在一个非零元, 取置换矩阵 P1 和 P2 使得 P1AP2 = " a11 A12 A21 A22# , 其中 a11 ̸= 0, A22 ∈ R (n−1)×(n−1) . 令 u11 = a11, U12 = A12, L21 = A21/a11, U22 = A22 − L21U12. 则 u11 ̸= 0, 且有 " 1 0 L21 I # "u11 U12 0 U22# = " a11 A12 A21 A22# = P1AP2. 两边取行列式可得 0 ̸= det(P1AP2) = det " 1 0 L21 I #! · det "u11 U12 0 U22#! = a11 · det(U22). 所以 det(U22) ̸= 0, 即U22 ∈ R (n−1)×(n−1) 非奇异. 由归纳假设可知, 存在置换矩阵 P˜ L ∈ R (n−1)×(n−1) 和 P˜R ∈ R (n−1)×(n−1) , 使得 P˜ LU22P˜R = L˜ 22U˜ 22, 其中 L˜ 22 为单位下三角矩阵, U˜ 22 为非奇异上三角矩阵. 令 PL = " 1 0 0 P˜ L # P1, PR = P2 " 1 0 0 P˜R # , 则有 PLAPR = " 1 0 0 P˜ L # P1AP2 " 1 0 0 P˜R # = " 1 0 0 P˜ L # " 1 0 L21 I # "u11 U12 0 U22# "1 0 0 P˜R # = " 1 0 P˜ LL21 P˜ L # "u11 U12 0 P˜T L L˜ 22U˜ 22P˜T R # "1 0 0 P˜R # = " 1 0 P˜ LL21 P˜ L # "1 0 0 P˜T 1 L˜ 22# "u11 U12 0 U˜ 22P˜T R # "1 0 0 P˜R # = " 1 0 P˜ LL21 L˜ 22# "u11 U12P˜R 0 U˜ 22 # ≜ LU, 其中 L 为单位下三角矩阵, U 为非奇异上三角矩阵, 即 A 存在 LU 分解. 因此, 由数学归纳法可知, 结论成立. □