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

第二讲线性方程组直接方法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