第二讲线性方程组直接方法Linear algebrain particular, the solution oflinear systems of equations —lies at the heart of mostcalculations in scientific computing.Dongarra&Eijkhout [32], 2000.考虑线性方程组(2.1)Ar=b,AeRnxn,bERn线性方程组的求解有着非常广泛的应用背景，科学计算中的很多问题最后都可能归结为求解一个或多个线性方程组.从纯数学角度来看，这个问题已经得到了完美的解决，因为它的解可以通过行列式直接表示出来，即Cramer法则.但在实际计算中，由于运算量增长速度太快，当n较大时，用Cramer法则解线性方程组是不可行的.另外，由于实际计算中的舍人误差，可能会导致一系列非常严重的问题一个从纯数学角度看似非常简单的问题，实际计算时可能会非常困难，有时甚至可能是一个无法解决的问题一般来说，求解线性方程组的数值方法可以分为两类：直接法与选代法.本讲介绍直接法，即Gauss消去法.直接法具有良好的稳定性和健壮性，而且在有限步内终止，因此在工程界很受欢迎.但由于运算量是O(n3),对于大规模问题，所需时间会很长(这里n表示未知量的个数).目前Gauss消去法是求解中小规模线性方程组或某些具有特殊结构的大规模稀疏线性方程组的首选方法凸Gauss消去法的思想可以追溯到公元一世纪左右的《九章算术》，Newton,Gauss,Lagrange等数学家都对该方法做出了贡献，相关历史可以参见[63]关于线性方程组直接法的相关参考文献G.H.Golub and C.F.Van Loan, Matrix Computations, 4th, 2013.[57]J.W.Demmel, Applied Numerical Linear Algebra, 1997. [30]L. N. Trefethen and D. Bau, II, Numerical Linear Algebra, 1997. [125]I. S.Duff, A. M.Erisman and J.K.Reid, Direct Methods for Sparse Matrices, 2nd, 2017.[38]T. A.Davis, Direct Methods for Sparse Linear Systems, SIAM, 2006. [29]后面两个文献主要是介绍大规模稀疏线性方程组的直接解法在本讲中，我们总是假定系数矩阵A是非奇异的，即线性方程组（2.1)的解存在且唯一，另外，为了讨论方便，我们只考虑实数情形，对于复系数线性方程组，其求解方法是类似的，41
第二讲 线性方程组直接方法 Linear algebra — in particular, the solution of linear systems of equations — lies at the heart of most calculations in scientific computing. — Dongarra & Eijkhout [32], 2000. 考虑线性方程组 Ax = b, A ∈ R n×n , b ∈ R n . (2.1) 线性方程组的求解有着非常广泛的应用背景, 科学计算中的很多问题最后都可能归结为求解一个 或多个线性方程组. 从纯数学角度来看, 这个问题已经得到了完美的解决, 因为它的解可以通过行 列式直接表示出来, 即 Cramer 法则. 但在实际计算中, 由于运算量增长速度太快, 当 n 较大时, 用 Cramer 法则解线性方程组是不可行的. 另外, 由于实际计算中的舍入误差, 可能会导致一系列非 常严重的问题. b 一个从纯数学角度看似非常简单的问题, 实际计算时可能会非常困难, 有时甚至可能是一 个无法解决的问题. 一般来说, 求解线性方程组的数值方法可以分为两类: 直接法与迭代法. 本讲介绍直接法, 即 Gauss 消去法. 直接法具有良好的稳定性和健壮性, 而且在有限步内终止, 因此在工程界很受欢 迎. 但由于运算量是 O(n 3 ), 对于大规模问题, 所需时间会很长 (这里 n 表示未知量的个数). 目前, Gauss 消去法是求解中小规模线性方程组或某些具有特殊结构的大规模稀疏线性方程组的首选方 法. b Gauss 消去法的思想可以追溯到公元一世纪左右的《九章算术》, Newton, Gauss, Lagrange 等 数学家都对该方法做出了贡献, 相关历史可以参见 [63]. 关于线性方程组直接法的相关参考文献 ▶ G. H. Golub and C. F. Van Loan, Matrix Computations, 4th, 2013. [57] ▶ J.W. Demmel, Applied Numerical Linear Algebra, 1997. [30] ▶ L. N. Trefethen and D. Bau, III, Numerical Linear Algebra, 1997. [125] ▶ I. S. Duff, A. M. Erisman and J. K. Reid, Direct Methods for Sparse Matrices, 2nd, 2017. [38] ▶ T. A. Davis, Direct Methods for Sparse Linear Systems, SIAM, 2006. [29] 后面两个文献主要是介绍大规模稀疏线性方程组的直接解法. 在本讲中, 我们总是假定系数矩阵 A 是非奇异的, 即线性方程组 (2.1) 的解存在且唯一. 另外, 为了讨论方便, 我们只考虑实数情形, 对于复系数线性方程组, 其求解方法是类似的. 41
:42第二讲线性方程组直接方法2.1LU分解与Gauss消去法2.1.1LU分解Gauss消去法本质上就是对系数矩阵A进行LU分解，即将A分解成两个矩阵的乘积A= LU,(2.2)其中L是单位下三角矩阵，U为非奇异上三角矩阵，然后再对两个三角方程进行求解：假定矩阵A存在LU分解（2.2)，则方程组(2.1)就转化为求解下面两个三角方程组[ Ly = b,Ur=y显然，上式中的两个三角方程组都非常容易求解将一个复杂问题分解成若干相对简单的问题，是矩阵分解的一个主要目的基于LU分解的Gauss消去法描述如下：算法2.1.Gauss消去法1：对A进行LU分解：A=LU，其中L为单位下三角矩阵，U为非奇异上三角矩阵；2:利用向前回代,求解Ly=b,即得y=L-1b;3:利用向后回代,求解Uc=y,即得a=U-1y=(LU)-1b=A-1b.我们知道，当系数矩阵A非奇异时，方程组（2.1）总是存在唯一解.但是，并不是每个非奇异矩阵都存在LU分解定理2.1(LU分解的存在性和唯一性）矩阵AERnxn存在LU分解（即存在单位下三角矩阵L和非奇异上三角矩阵U使得A=LU)的充要条件是A的所有顺序主子矩阵都非奇异.进一步若A存在LU分解，则分解是唯一的(板书)证明.必要性：设A11是A的k阶顺序主子矩阵，其中1<k≤n.将A=LU写成分块形式，即[U1[L110U12[LnUiA11A12LiiU12U22A21A22L22]0L21U11 L2iU12+L22U22L21可得A11=L11U11.由于L11和Ui1均非奇异，所以A11也非奇异充分性：用归纳法，当n=1时，结论显然成立假设结论对所有n－1阶矩阵都成立，即对任意n－1阶矩阵，如果其所有的顺序主子矩阵都非奇异，则存在LU分解.考虑n阶的矩阵A，写成分块形式AtA22A21
· 42 · 第二讲 线性方程组直接方法 2.1 LU 分解与 Gauss 消去法 2.1.1 LU 分解 Gauss 消去法本质上就是对系数矩阵 A 进行 LU 分解, 即将 A 分解成两个矩阵的乘积 A = LU, (2.2) 其中 L 是单位下三角矩阵, U 为非奇异上三角矩阵, 然后再对两个三角方程进行求解: 假定矩阵 A 存在 LU 分解 (2.2), 则方程组 (2.1) 就转化为求解下面两个三角方程组 Ly = b, Ux = y. 显然, 上式中的两个三角方程组都非常容易求解. b 将一个复杂问题分解成若干相对简单的问题, 是矩阵分解的一个主要目的. 基于 LU 分解的 Gauss 消去法描述如下: 算法 2.1. 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 非奇异时, 方程组 (2.1) 总是存在唯一解. 但是, 并不是每个非奇异矩 阵都存在 LU 分解. 定理 2.1 (LU 分解的存在性和唯一性) 矩阵 A ∈ R n×n 存在 LU 分解 (即存在单位下三角矩阵 L 和非奇异上三角矩阵 U 使得 A = LU) 的充要条件是 A 的所有顺序主子矩阵都非奇异. 进一步, 若 A 存在 LU 分解, 则分解是唯一的. (板书) 证明. 必要性: 设 A11 是 A 的 k 阶顺序主子矩阵, 其中 1 ≤ k ≤ n. 将 A = LU 写成分块形式, 即 " A11 A12 A21 A22# = " L11 0 L21 L22# "U11 U12 0 U22# = " L11U11 L11U12 L21U11 L21U12 + L22U22# . 可得 A11 = L11U11. 由于 L11 和 U11 均非奇异, 所以 A11 也非奇异. 充分性: 用归纳法. 当 n = 1 时, 结论显然成立. 假设结论对所有 n − 1 阶矩阵都成立, 即对任意 n − 1 阶矩阵, 如果其所有的顺序主子矩阵都 非奇异, 则存在 LU 分解. 考虑 n 阶的矩阵 A, 写成分块形式 A = " A11 A12 A21 A22#
2.1LU分解与Gauss消去法:43.其中 A11ER(n-1)x(n-1)是A的 n－1 阶顺序主子矩阵.由归纳假设可知,A11 存在LU分解,即存在单位下三角矩阵L11和非奇异上三角矩阵U11使得An= LnUl.令L21=A21U7l,U12=LilA12，U22=A22-L21U12则[LiiUn[Ui1U12LiiU12[A11A12L110U22].0[L21U11 U22+L21U12][A21 A22]L21[L11 0][U11U12为单位下三角矩阵，U会因此可得A的LU分解A=LU,其中L会为非L211U220奇异的上三角矩阵（U的非奇异性可由A的非奇异性可得）由归纳法可知，结论成立，下面证明唯一性.设A存在两个不同的LU分解：A=LU=LU其中L和L为单位下三角矩阵，U和U为非奇异上三角矩阵.则有L-1L=UU-1该等式左边为下三角矩阵，右边为上三角矩阵，所以只能是对角矩阵.又单位下三角矩阵的逆仍然是单位下三角矩阵，所以L-1i的对角线元素全是1,故L-i = I,口即L=L,U=U.唯一性得证记D为U的对角线部分，则A=LDü,其中ü=D-1U是单位上三角矩阵.因此我们就有下面的LDU分解推论2.2(LDU分解）设AERn×n的所有顺序主子矩阵都非奇异,则A存在LDU分解，即存在单位下三角矩阵L和单位上三角矩阵U，以及非奇异对角矩阵D,使得A=LDU,其中L,U,D都是唯一的.反之，若A存在LDU分解，则A的所有顺序主子矩阵都非奇异对角占优情形一般的非奇异矩阵不一定存在LU分解，但如果A是对角占优的，则存在LU分解.我们可以证明，严格对角占优矩阵在LU分解中保持严格对角占优性（见课后练习），因此通过数学归纳法可以证明严格对角占优矩阵一定存在LU分解.事实上，只要矩阵列对角占优且非奇异，则一定存在LU分解[57].定理2.3[57] 设AeRnxn非奇异且列对角占优,则A存在 LU分解且 L中的元素的绝对值都不超过1.（留作练习，数学归纳法）
2.1 LU 分解与 Gauss 消去法 · 43 · 其中 A11 ∈ R (n−1)×(n−1) 是 A 的 n − 1 阶顺序主子矩阵. 由归纳假设可知, A11 存在 LU 分解, 即 存在单位下三角矩阵 L11 和非奇异上三角矩阵 U11 使得 A11 = L11U11. 令 L21 = A21U −1 11 , U12 = L −1 11 A12, U22 = A22 − L21U12, 则 " L11 0 L21 1 # "U11 U12 0 U22# = " L11U11 L11U12 L21U11 U22 + L21U12# = " A11 A12 A21 A22# = A. 因此可得 A 的 LU 分解 A = LU, 其中 L ≜ " L11 0 L21 1 # 为单位下三角矩阵, U ≜ " U11 U12 0 U22# 为非 奇异的上三角矩阵 (U 的非奇异性可由 A 的非奇异性可得). 由归纳法可知, 结论成立. 下面证明 唯一性. 设 A 存在两个不同的 LU 分解: A = LU = L˜U, ˜ 其中 L 和 L˜ 为单位下三角矩阵, U 和 U˜ 为非奇异上三角矩阵. 则有 L −1L˜ = UU˜ −1 , 该等式左边为下三角矩阵, 右边为上三角矩阵, 所以只能是对角矩阵. 又单位下三角矩阵的逆仍然 是单位下三角矩阵, 所以 L −1L˜ 的对角线元素全是 1, 故 L −1L˜ = I, 即 L˜ = L, U˜ = U. 唯一性得证. □ 记 D 为 U 的对角线部分, 则 A = LDU˜, 其中 U˜ = D−1U 是单位上三角矩阵. 因此我们就有 下面的 LDU 分解. 推论 2.2 (LDU 分解) 设 A ∈ R n×n 的所有顺序主子矩阵都非奇异, 则 A 存在 LDU 分解, 即存在 单位下三角矩阵 L 和单位上三角矩阵 U, 以及非奇异对角矩阵 D, 使得 A = LDU, 其中 L, U, D 都是唯一的. 反之, 若 A 存在 LDU 分解, 则 A 的所有顺序主子矩阵都非奇异. 对角占优情形 一般的非奇异矩阵不一定存在 LU 分解, 但如果 A 是对角占优的, 则存在 LU 分解. 我们可以 证明, 严格对角占优矩阵在 LU 分解中保持严格对角占优性 (见课后练习), 因此通过数学归纳法可 以证明严格对角占优矩阵一定存在 LU 分解. 事实上, 只要矩阵列对角占优且非奇异, 则一定存在 LU 分解 [57]. 定理 2.3 [57] 设 A ∈ R n×n 非奇异且列对角占优, 则 A 存在 LU 分解且 L 中的元素的绝对值都 不超过 1. (留作练习, 数学归纳法)
-44.第二讲线性方程组直接方法需要注意的是，在实际计算中，由于舍入误差的影响，即使是严格对角占优矩阵，如果对角占优性不是很强的话，很有可能会失去对角占优性，从而导致LU分解失败2.1.2LU分解的实现矩阵的LU分解可以通过初等行变换来实现.给定矩阵a11a12+++aina21a22a2nERnxnA=.....anlan2...ann·第一步：假定a11≠0,构造矩阵..0[10021100.ail01..0其中li=Li =131.i=2,3,...,n.a11..：[Ini001...易知Li的逆为01000-12110L=1 =０10-131...:.-ln100***1用L-1左乘 A,并将所得到的矩阵记为 A(1),则a1la12...ainaga0..A(I) = L-"A...a.. am0即左乘L-1后，A的第一列中除第一个元素外其它都变为0.·第二步：类似地，我们可以将上面的操作作用在A(1)的子矩阵A(1)(2：n,2：n）上，将其第一列除第一个元素外都变为0.也就是说，假定≠0,构造矩阵[1000...0010...(1)a2013210其中L2 =li2 =0 ...："...[0n20..1用 L=1左乘 A(1),并将所得到的矩阵记为A(2),则
· 44 · 第二讲 线性方程组直接方法 需要注意的是, 在实际计算中, 由于舍入误差的影响, 即使是严格对角占优矩阵, 如果对角占 优性不是很强的话, 很有可能会失去对角占优性, 从而导致 LU 分解失败. 2.1.2 LU 分解的实现 矩阵的 LU 分解可以通过初等行变换来实现. 给定矩阵 A = a11 a12 · · · a1n a21 a22 · · · a2n . . . . . . an1 an2 · · · ann ∈ R n×n . • 第一步: 假定 a11 ̸= 0, 构造矩阵 L1 = 1 0 0 · · · 0 l21 1 0 · · · 0 l31 0 1 · · · 0 . . . . . . ln1 0 0 · · · 1 , 其中 li1 = ai1 a11 , i = 2, 3, . . . , n. 易知 L1 的逆为 L −1 1 = 1 0 0 · · · 0 −l21 1 0 · · · 0 −l31 0 1 · · · 0 . . . . . . −ln1 0 0 · · · 1 . 用 L −1 1 左乘 A, 并将所得到的矩阵记为 A(1) , 则 A (1) = L −1 1 A a11 a12 · · · a1n 0 a (1) 22 · · · a (1) 2n . . . . . . . . . 0 a (1) n2 · · · a (1) nn . 即左乘 L −1 1 后, A 的第一列中除第一个元素外其它都变为 0. • 第二步: 类似地, 我们可以将上面的操作作用在 A(1) 的子矩阵 A(1)(2 : n, 2 : n) 上, 将其第一 列除第一个元素外都变为 0. 也就是说, 假定 a (1) 22 ̸= 0, 构造矩阵 L2 = 1 0 0 · · · 0 0 1 0 · · · 0 0 l32 1 · · · 0 . . . . . . . . . 0 ln2 0 · · · 1 , 其中 li2 = a (1) i2 a (1) 22 , i = 3, 4, . . . , n. 用 L −1 2 左乘 A(1) , 并将所得到的矩阵记为 A(2) , 则
2.1LU分解与Gauss消去法-45.aina11a12a13...oa0...(2)A(2) = L,IA(1) = L,"L-1 A =00...a33........t.a]ae00...α(k-1) + 0 (k = 3,4, .., n -·依此类推，假定ak-1)，则我们可以构造一系列的矩阵L3,L4....Ln-1,使得a11a12a13ainaaad)0aaa..00L--....L"L-"A=:...:.alm-1)000+.为一个上三角矩阵.我们将这个上三角矩阵记为U，并记[1000...1l2100...1321..0L≤LL2...Ln-1=131(2.3)目..:IniIn2ln3...1则可得A= LU,这就是A的LU分解将上述过程写成算法，描述如下：算法2.2.LU分解1: Set L=I,U =0%将L设为单位矩阵，U设为零矩阵2: for k=1 to n-1 do3:for j = k to n do4:%更新U的第k行ukj =akjend for5:6:fori=k+1 to n do7:lk=aikakk%计算8:for j=k+ltondo9:%更新A(i,k+1:n)a=a-lukj10:end for
2.1 LU 分解与 Gauss 消去法 · 45 · A (2) = L −1 2 A (1) = L −1 2 L −1 1 A = a11 a12 a13 · · · a1n 0 a (1) 22 a (1) 23 · · · a (1) 2n 0 0 a (2) 33 · · · a (2) 3n . . . . . . . . . . . . 0 0 a (2) n3 · · · a (2) nn . • 依此类推, 假定 a (k−1) kk ̸= 0 (k = 3, 4, . . . , n − 1), 则我们可以构造一系列的矩阵 L3, L4, . . . , Ln−1, 使得 L −1 n−1 · · ·L −1 2 L −1 1 A = a11 a12 a13 · · · a1n 0 a (1) 22 a (1) 23 · · · a (1) 2n 0 0 a (2) 33 · · · a (2) 3n . . . . . . . . . . . . 0 0 0 · · · a (n−1) nn 为一个上三角矩阵. 我们将这个上三角矩阵记为 U, 并记 L ≜ L1L2 · · ·Ln−1 = 1 0 0 · · · 0 l21 1 0 · · · 0 l31 l32 1 · · · 0 . . . . . . . . . ln1 ln2 ln3 · · · 1 , (2.3) 则可得 A = LU, 这就是 A 的 LU 分解. 将上述过程写成算法, 描述如下: 算法 2.2. LU 分解 1: Set L = I, U = 0 % 将 L 设为单位矩阵, U 设为零矩阵 2: for k = 1 to n − 1 do 3: for j = k to n do 4: ukj = akj % 更新 U 的第 k 行 5: end for 6: for i = k + 1 to n do 7: lik = aik/akk % 计算 lik 8: for j = k + 1 to n do 9: aij = aij − likukj % 更新 A(i, k + 1 : n) 10: end for
: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 分解. 因此, 由数学归纳法可知, 结论成立. □