华东师范大学：《矩阵计算》课程教学资源（课堂讲义）第二讲 线性方程组直接方法

第二讲线性方程组直接方法一般来说，求解线性方程组的数值方法可以分为两类：直接法与选代法.本章介绍直接法，即Gauss消去法.直接法相对比较稳定，因此在工程计算中很受欢迎.但由于运算量是O（n3），当问题规模较大时，时间会很长（这里n表示未知量的个数）.目前，直接法主要用于小规模或中等规模线性方程组的数值求解2.1Gauss消去法和LU分解2.1.1LU分解考虑线性方程组Ar= b,(2.1)其中AERnxn非奇异，bER"为给定的右端项.Gauss消去法本质上就是对系数矩阵A进行LU分解，即将A分解成两个矩阵的乘积A = LU,(2.2)其中L是单位下三角矩阵，U为非奇异上三角矩阵.这个分解就称为LU分解假定矩阵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-1b3:利用向后回代，求解U=9,即得=U-1y=（LU)-1b=A-1b我们知道，当系数矩阵A非奇异时，方程组（2.1）总是存在唯一解但是，并不是每个非奇异矩阵都存在LU分解定理2.1（LU分解的存在性和唯一性）设AERnXn.则存在唯一的单位下三角矩阵L和非奇异上三角矩阵U，使得A=LU的充要条件是A的所有顺序主子矩阵A=A(1：k,1：k)都非奇异k=1,2,...,n.37

第二讲 线性方程组直接方法 一般来说, 求解线性方程组的数值方法可以分为两类: 直接法与迭代法. 本章介绍直接法, 即 Gauss 消去法. 直接法相对比较稳定, 因此在工程计算中很受欢迎. 但由于运算量是 O(n 3 ), 当问题规模较大 时, 时间会很长 (这里 n 表示未知量的个数). 目前, 直接法主要用于小规模或中等规模线性方程组的数 值求解. 2.1 Gauss 消去法和 LU 分解 2.1.1 LU 分解 考虑线性方程组 Ax = b, (2.1) 其中 A ∈ R n×n 非奇异, b ∈ R n 为给定的右端项. Gauss 消去法本质上就是对系数矩阵 A 进行 LU 分解, 即将 A 分解成两个矩阵的乘积 A = LU, (2.2) 其中 L 是单位下三角矩阵, U 为非奇异上三角矩阵. 这个分解就称为 LU 分解. 假定矩阵 A 存在 LU 分解 (2.2), 则方程组 (2.1) 就转化为求解下面两个三角方程组 { Ly = b, Ux = y. 显然, 这两个方程组都非常容易求解. 基于 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. 则存在唯一的单位下三角矩阵 L 和非奇异上 三角矩阵 U, 使得 A = LU 的充要条件是 A 的所有顺序主子矩阵 Ak = A(1 : k, 1 : k) 都非奇异, k = 1, 2, . . . , n. 37

:38.第二讲线性方程组直接方法证明.必要性：设A11是A的k阶顺序主子矩阵，将A=LU写成分块形式LnU12A11A12LnUiiL11U11U12L21U11L21U12+L22U22A21A22L210U22可得A11=LiiUi1.由于L11和Uii均非奇异，所以A11也非奇异充分性：用归纳法当n=1时，结论显然成立假设结论对n－1阶矩阵都成立，即对任意n－1阶矩阵，如果其所有的顺序主子矩阵都非奇异，则存在LU分解考虑n阶的矩阵A,写成分块形式21其中A11ER(n-1)x(n-1)是A的n-1阶顺序主子矩阵.由归纳假设可知,A11存在LU分解,即存在单位下三角矩阵L11和非奇异上三角矩阵U11使得A1 = LiuUi.令L21=A21U-l,U12=LiA12,U22=A22-L21U12则[A11.A12[LuUn[L1L11U120U11ALU1A21A22L2iUi1U22+L21U12L210U22易知U非奇异，所以A存在LU分解，下面证明唯一性设A存在两个不同的LU分解A=LU=LU其中L和工为单位下三角矩阵，U和ü为非奇异上三角矩阵.则有L-1L=UU-1该等式左边为下三角矩阵，右边为上三角矩阵，所以只能是对角矩阵，由于单位下三角矩阵的逆仍然是单位下三角矩阵，所以L-1L的对角线元素全是1，故L-1i= I,即i=L,U=U.口由归纳法可知，结论成立

· 38 · 第二讲 线性方程组直接方法 证明. 必要性: 设 A11 是 A 的 k 阶顺序主子矩阵, 将 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] , 其中 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, 则 A = [ A11 A12 A21 A22] = [ L11U11 L11U12 L21U11 U22 + L21U12] = [ L11 0 L21 1 ] [U11 U12 0 U22] ≜ LU. 易知 U 非奇异, 所以 A 存在 LU 分解. 下面证明唯一性. 设 A 存在两个不同的 LU 分解: A = LU = L˜U, ˜ 其中 L 和 L˜ 为单位下三角矩阵, U 和 U˜ 为非奇异上三角矩阵. 则有 L −1L˜ = UU˜ −1 , 该等式左边为下三角矩阵, 右边为上三角矩阵, 所以只能是对角矩阵. 由于单位下三角矩阵的逆仍然是 单位下三角矩阵, 所以 L −1L˜ 的对角线元素全是 1, 故 L −1L˜ = I, 即 L˜ = L, U˜ = U. 由归纳法可知, 结论成立. □

2.1Gauss消去法和LU分解·39 .2.1.2LU分解的实现给定一个矩阵a11α12aina21a22a2nERnxnA=:..[anl...an2ann我们可以通过矩阵初等变换来构造A的LU分解，·第一步：假定a11≠0,构造矩阵[100...012110..0Qil01.0其中l31Li =lil:.i=2,3,...,n.a11.[ni00.1易知L的逆为0001121100Lil =l3101...0:..-ln1001...用L-1左乘A,并将所得到的矩阵记为A(1),则a11α12ain(1)(1)0022A(1) = L-"A.....a0即左乘L-1后，A的第一列中除第一个元素外其它都变为0.·第二步：类似地，我们可以将上面的操作作用在A(1）的子矩阵A(1)(2：n.2：n）上，将其第一列除第一个元素外都变为0.也就是说，假定a≠0.构造矩阵100700100al)..13210C其中L2 =1.2 ...：..01n20...1

2.1 Gauss 消去法和 LU 分解 · 39 · 2.1.2 LU 分解的实现 给定一个矩阵 A =        a11 a12 · · · a1n a21 a22 · · · a2n . . . . . . an1 an2 · · · ann        ∈ R n×n . 我们可以通过矩阵初等变换来构造 A 的 LU 分解. • 第一步: 假定 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

:40 .第二讲线性方程组直接方法用L-1左乘A(1)并将所得到的矩阵记为A(2)，则[a11a12aina13.e号ageen000A(2) = L-1 A(1) = L-1L-1 A =...:aea2]Lo0....依此类推，假定ak-1)0(k=3,4,n-1),则我们可以构造一系列的矩阵L3,L4....，Ln-1,使得[a11α12α13ain.ago0...00L---. L?"L-'A -..........alm-1)]L 000...为一个上三角矩阵.我们将这个上三角矩阵记为U，并记n000..10l21...013113210...L = L,L2-.- Ln-1 -(2.3)......"..[n1In3...In21则可得A= LU,这就是A的LU分解将上面的过程写成算法，描述如下：算法2.2.LU分解1:SetL=I,U=0%将L设为单位矩阵，U设为零矩阵2:fork=1ton-1do3:fori=k+1tondo4:%计算L的第K列li=ai/akkend for5:6:for j=k ton do7:%计算U的第人行ukj = akjend for8:for i=k+1 to n do9:for j=k+1 to ndo10:11:%更新A(k+1:n,k+1:n)aij = aij -likukj12:end for

· 40 · 第二讲 线性方程组直接方法 用 L −1 2 左乘 A(1) , 并将所得到的矩阵记为 A(2) , 则 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 i = k + 1 to n do 4: lik = aik/akk % 计算 L 的第 k 列 5: end for 6: for j = k to n do 7: ukj = akj % 计算 U 的第 k 行 8: end for 9: for i = k + 1 to n do 10: for j = k + 1 to n do 11: aij = aij − likukj % 更新 A(k + 1 : n, k + 1 : n) 12: end for

2.1Gauss消去法和LU分解41.13:end for14:end forGauss消去法的运算量由算法2.2可知，LU分解的运算量（含加减乘除）为n-17121+Z(n -i+2(n-i) =n3+0(n2)13j=i+1k=i+1i=1i=i+1n3 +0(n2)由于回代过程的运算量为O(n2)，所以Gauss消去法的总运算量为G十评价算法的一个主要指标是执行时间，但这依赖于计算机硬件和编程技巧等，因此直接给出算法执行时间是不太现实的.所以我们通常是统计算法中算术运算（加减乘除）的次数在数值算法中，大多仅仅涉及加减乘除和开方运算，一般地，加减运算次数与乘法运算次数具有相同的量级而除法运算和开方运算次数具有更低的量级十为了尽可能地减少运算量，在实际计算中，数，向量和矩阵做乘法运算时的先后执行次序为：先计算数与向量的乘法，然后计算矩阵与向量的乘法，最后才计算矩阵与矩阵的乘法矩阵L和U的存储当A的第i列被用于计算L的第i列后，在后面的计算中不再被使用.同样地，A的第i行被用于计算U的第行后，在后面的计算中也不再被使用，因此，为了节省存储空间，我们可以在计算过程中将L的第i列存放在A的第i列，将U的第i行存放在A的第i行，这样就不需要另外分配空间存储L和U.计算结束后，A的上三角部分为U，其绝对下三角部分为L的绝对下三角部分（L的对角线全部为1，不需要存储).此时算法可以描述为：算法2.3.LU分解（用A存储L和U)l:fork=1ton-1do2:fori=k+1tondo3:aik = aik/akk4:for j=k+1 to n do5:aij=aj-aiakj6:endforend for7:8:endfor十根据指标的循环次序，算法2.3也称为KIJ型LU分解.在实际计算中，我们一般不建议使用这个算法.因为对于指标k的每次循环，都需要更新A的第k十1至第n行.这种反复读取数据的做法会使得计算效率大大降低

2.1 Gauss 消去法和 LU 分解 · 41 · 13: end for 14: end for Gauss 消去法的运算量 由算法 2.2 可知, LU 分解的运算量 (含加减乘除) 为 n∑−1 i=1   ∑n j=i+1 1 + ∑n j=i+1 ∑n k=i+1 2   = n∑−1 i=1 ( n − i + 2(n − i) 2 ) = 2 3 n 3 + O(n 2 ). 由于回代过程的运算量为 O(n 2 ), 所以 Gauss 消去法的总运算量为 2 3 n 3 + O(n 2 ). † 评价算法的一个主要指标是执行时间, 但这依赖于计算机硬件和编程技巧等, 因此直接给出算法 执行时间是不太现实的. 所以我们通常是统计算法中算术运算 (加减乘除) 的次数. 在数值算法 中, 大多仅仅涉及加减乘除和开方运算. 一般地, 加减运算次数与乘法运算次数具有相同的量级, 而除法运算和开方运算次数具有更低的量级. † 为了尽可能地减少运算量, 在实际计算中, 数, 向量和矩阵做乘法运算时的先后执行次序为: 先计 算数与向量的乘法, 然后计算矩阵与向量的乘法, 最后才计算矩阵与矩阵的乘法. 矩阵 L 和 U 的存储 当 A 的第 i 列被用于计算 L 的第 i 列后, 在后面的计算中不再被使用. 同样地, A 的第 i 行被用于 计算 U 的第 i 行后, 在后面的计算中也不再被使用. 因此, 为了节省存储空间, 我们可以在计算过程中将 L 的第 i 列存放在 A 的第 i 列, 将 U 的第 i 行存放在 A 的第 i 行, 这样就不需要另外分配空间存储 L 和 U. 计算结束后, A 的上三角部分为 U, 其绝对下三角部分为 L 的绝对下三角部分 (L 的对角线全部为 1, 不需要存储). 此时算法可以描述为： 算法 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 † 根据指标的循环次序, 算法 2.3 也称为 KIJ 型 LU 分解. 在实际计算中, 我们一般不建议使用这个 算法. 因为对于指标 k 的每次循环, 都需要更新 A 的第 k + 1 至第 n 行. 这种反复读取数据的做 法会使得计算效率大大降低