数值代数 夏银华 中国科学技术大学
扩充矩阵 给定矩阵 A 和向量 b a 11 a 12 a 13 b 1 A = a 21 a 22 a 23 b = b 2 a 31 a 32 a 33 b 3 定义扩充矩阵 ( A b ) = a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33 b 1 b 2 b 3
初等变换矩阵 放缩 (scaling): 第 i 个方程 (E i ) 乘上任意非零常数 λ, 记作 (λe i ) (E i ) 行加法 (scaled addition): 第 i 个方程 (E i ) 加上第 j 个方程 (E j ) 的 λ 倍, 记作 (E i + λe j ) (E i ) 行交换 (reordering): 第 i 个方程 (E i ) 与第 j 个方程 (E j ) 交换位置, 记作 (E j ) (E i )
高斯消元 目的是将扩充矩阵 高斯消元与 ( A b ) 利用初等变换, 成为上三角矩阵 ã 11 ã 12 ã 13 0 ã 22 ã 23 0 0 ã 33 b 1 b 2 b 3
高斯消元 目的是将扩充矩阵 高斯消元与 ( A b ) 利用初等变换, 成为上三角矩阵 ã 11 ã 12 ã 13 0 ã 22 ã 23 0 0 ã 33 b 1 b 2 b 3 再利用向后替换算法求解 x 3 = b 3 /ã 33 x 2 = ( b 2 ã 23 x 3 )/ã 22 x 1 = ( b 1 ã 12 x 2 ã 13 x 3 )/ã 11
通过行向量之间的初等变换, 将 A C m n 转化成上三角阵 U
通过行向量之间的初等变换, 将 A C m n 转化成上三角阵 U 下三角阵 L i 将第 i 列对角元素以下变成 0: 其中 L = L 1 1 L 1 2 L 1 m 1 L m 1 L 2 L }{{} 1 A = U A = LU L 1
通过行向量之间的初等变换, 将 A C m n 转化成上三角阵 U 下三角阵 L i 将第 i 列对角元素以下变成 0: 其中 L = L 1 1 L 1 2 L 1 m 1 三角三角化 L m 1 L 2 L }{{} 1 A = U A = LU L 1
rows Each L i introduces zeros below diagonal of column i: 高斯消元与 L m 1 L 2 L }{{} 1 A = U = A = LU where L = L 1 1 L 1 2 L 1 m 1 L 1 L 1 0 0 0 L 2 0 0 L 3 0 A L 1 A L 2 L 1 A L 3 L 2 L 1 A Triangular triangularization 第 k 步, 2 a k = ( a 1k a kk a k+1,k a mk ) T L k a k = ( a 1k a kk 0 0 ) T
令 l jk = a jk /a kk, 1 L k = 1 l k+1,k 1 l mk 1 L k = I l k e k l k = ( 0,, 0, l k+1,k, l m,k ) T
L = L 1 1 L 1 2 L 1 m 1 1 l 21 1 L = l 31 l 32 1........ l m1 l m2 l m,m 1 1 L 1 k = I + l k e k
Factorize A C m m into A = LU : 对 A C m m 进行 A = LU: 高斯消元与 Gaussian Elimination without Pivoting Algorithm: Gaussian Elimination (no pivoting) U = A, L = I for k = 1 to m 1 for j = k + 1 to m l jk = u jk /u kk u j,k:m = u j,k:m l jk u k,k:m 工作量 m 1 2(m k)(m k) 2m 3 /3 The inner loop can be written using matrix operations instead of for-loop k=1 Operation count m k=1 2(m k)(m k) 2 m k=1 k2 2m 3 /3 5
Pivoting 高斯消元与 At step k, we used matrix element k, k as pivot and introduced zeros in 第 entry k 步时 k of,x remaining kk 作为主元 rows (pivot) 对第 k 列进行消元 x kk x kk 0 0 0 But any other element i k in column k can be used as pivot: x ik 0 0 x ik 0 6
Pivoting 高斯消元与 At At step step k, k, we we used used matrix matrix element element k, k as as pivot pivot and and introduced introduced zeros zeros in in 第 entry entry k 步时 k of of,x remaining remaining kk 作为主元 rows rows (pivot) 对第 k 列进行消元 x kk kk kk x kk 0 0 0 第 But But k 步时 any any other other,x element element ik, (i i k) k in in 也可以作为主元 column column k can can be be used used as as pivot: pivot: ik x ik 0 0 x ik ik 0 6
Pivoting Also, any other column 同样, 第 k 步时,x j ij, (i, k can be used: j k) 也可以作为主元 x ij 0 0 x ij 0 Choosing different pivots means we can avoid zero or very small pivots Instead of using pivots at different entries, change rows or columns and use the standard triangular algorithm (pivoting) A computer code might account for the pivoting indirectly instead of actually moving the data
Pivoting Also, any other column 同样, 第 k 步时,x j ij, (i, k can be used: j k) 也可以作为主元 x ij 0 0 x ij 0 选择不同的主元目的是为了避免 0 或者很小的主元 Choosing different pivots means we can avoid zero or very small pivots Instead of using pivots at different entries, change rows or columns and use the standard triangular algorithm (pivoting) A computer code might account for the pivoting indirectly instead of actually moving the data
Pivoting Also, any other column 同样, 第 k 步时,x j ij, (i, k can be used: j k) 也可以作为主元 x ij 0 0 x ij 0 选择不同的主元目的是为了避免 0 或者很小的主元 Choosing different pivots means we can avoid zero or very small pivots 选主元的方法可以通过对矩阵进行交换行 列位置达到 (pivoting) Instead of using pivots at different entries, change rows or columns and use the standard triangular algorithm (pivoting) A computer code might account for the pivoting indirectly instead of actually moving the data
高斯消元与 完全选主元法 (complete pivoting): 搜索最 好 的主元, 但代价太大
高斯消元与 完全选主元法 (complete pivoting): 搜索最 好 的主元, 但 Partial Pivoting 代价太大 部分选主元法 Searching among (partial all validpivoting): pivots is expensive 只考虑第 (complete k 列的元素作主 pivoting) 元, 只对行交换顺序 Consider pivots in column k only and interchange rows (partial pivoting) P 1 L 1 x ik x ik x ik 0 0 0 Pivot selection Row interchange Elimination In terms of matrices: L m 1 P m 1 L 2 P 2 L 1 P 1 A = U
高斯消元与 完全选主元法 (complete pivoting): 搜索最 好 的主元, 但 Partial Pivoting 代价太大 部分选主元法 Searching among (partial all validpivoting): pivots is expensive 只考虑第 (complete k 列的元素作主 pivoting) 元, 只对行交换顺序 Consider pivots in column k only and interchange rows (partial pivoting) P 1 L 1 x ik x ik x ik 0 0 0 Pivot selection 写成矩阵形式 : In terms of matrices: Row interchange Elimination L m 1 L m 1 P m 1 P L 2 P 2 2 L 1 P 1 A = U
重新写成 : (L m 1 L 2L 1)(P m 1 P 2 P 1 )A = U 其中 L k = P m 1 P k+1 L k P 1 k+1 P 1 m 1 此时 : PA = LU
Gaussian Elimination with Partial Pivoting Factorize A C m m into P A = LU : 部分选主元法 Algorithm: Gaussian Elimination (partial pivoting) U = A, L = I, P = I for k = 1 to m 1 Select i k to maximize u ik u k,k:m u i,k:m (interchange two rows) l k,1:k 1 l i,1:k 1 p k,: p i,: for j = k + 1 to m l jk = u jk /u kk u j,k:m = u j,k:m l jk u k,k:m 10
完全选主元法, 主元可能来自不同的列, 需要列置换矩阵 Q k 定义 L m 1 P m 1 L 2 P 2 L 1 P 1 AQ 1 Q 2 Q m 1 = U (L m 1 L 2L 1)(P m 1 P 2 P 1 )AQ 1 Q 2 Q m 1 = U L = (L m 1 L 2L 1) 1 P = P m 1 P 2 P 1 Q = Q 1 Q 2 Q m 1 完全选主元法可以写成 PAQ = LU
没有选主元的 : 稳定性 LŨ = A + δa, δa L U = O(ɛ M) 注意 : 是对 LŨ 的误差, 而不是对 L 或 Ũ 单个的误差
没有选主元的 : 稳定性 LŨ = A + δa, δa L U = O(ɛ M) 注意 : 是对 LŨ 的误差, 而不是对 L 或 Ũ 单个的误差 注意 : L U 在分母上, 而不是 A
没有选主元的 : 稳定性 LŨ = A + δa, δa L U = O(ɛ M) 注意 : 是对 LŨ 的误差, 而不是对 L 或 Ũ 单个的误差 注意 : L U 在分母上, 而不是 A L 和 U 可能任意大
没有选主元的 : 稳定性 LŨ = A + δa, δa L U = O(ɛ M) 注意 : 是对 LŨ 的误差, 而不是对 L 或 Ũ 单个的误差 注意 : L U 在分母上, 而不是 A L 和 U 可能任意大 因此, 算法不稳定
高斯消元与 选主元的 :L 中元素的绝对值都 1, 所以 L = O(1)
选主元的 :L 中元素的绝对值都 1, 所以 L = O(1) U 的大小, 引入增长因子 因此 U = O(ρ A ) ρ = max u i,j max a ij
选主元的 :L 中元素的绝对值都 1, 所以 L = O(1) U 的大小, 引入增长因子 因此 U = O(ρ A ) 对 PA = LU ρ = max u i,j max a ij LŨ = PA + δa, δa A = O(ρɛ M)
选主元的 :L 中元素的绝对值都 1, 所以 L = O(1) U 的大小, 引入增长因子 因此 U = O(ρ A ) 对 PA = LU ρ = max u i,j max a ij LŨ = PA + δa, δa A = O(ρɛ M) 如果 ρ = O(1), 算法向后稳定
增长因子 ρ 2 m 1, 对固定 m m 矩阵是 O(1)
增长因子 ρ 2 m 1, 对固定 m m 矩阵是 O(1) 因此, 按照定义算法向后稳定, 但是结果可能没用
增长因子 ρ 2 m 1, 对固定 m m 矩阵是 O(1) 因此, 按照定义算法向后稳定, 但是结果可能没用 实际应用中, 出于某种原因增长因子总是很小
HPD 矩阵 对称矩阵 : A = A T,A R m m Hermitian 矩阵 : A = A,A C m m
HPD 矩阵 对称矩阵 : A = A T,A R m m Hermitian 矩阵 : A = A,A C m m 对称正定矩阵 (SPD): 对称矩阵满足 x T Ax > 0, x 0 Hermitian 正定矩阵 (HPD):Hermitian 矩阵满足 x Ax > 0, x 0
HPD 矩阵 对称矩阵 : A = A T,A R m m Hermitian 矩阵 : A = A,A C m m 对称正定矩阵 (SPD): 对称矩阵满足 x T Ax > 0, x 0 Hermitian 正定矩阵 (HPD):Hermitian 矩阵满足 x Ax > 0, x 0 正定矩阵具有实的特征根和正交的特征向量
Cholesky 分解 Cholesky Factorization Eliminate below pivot and to the right of pivot: 同时消去主元以下列和右边行 : A = a 11 w = α 0 α w /α w K w/α I 0 K ww /a 11 = α 0 1 0 α w /α = R1A 1 R 1 w/α I 0 K ww /a 11 0 I 其中 where α α = aa 11 11 KK www /a /a 11 11 是 is principal R1 AR 1 submatrix 1 的主子式 of PD matrix R its upper-left entry is positive 1 AR1 1, therefore
Cholesky 分解 递归作用得到 : A = (R 1 R 2 R m)(r m R 2 R 1 ) = R R, r jj > 0 每个正定矩阵都有唯一的 Cholesky 分解 Cholesky 分解算法不会崩溃 α = a 11 唯一确定, 因此 w /α 给定
The Cholesky Factorization Algorithm Cholesky 分解 Factorize hermitian positive definite A C m m into A = R R: Cholesky 分解 Algorithm: Cholesky Factorization R = A for k = 1 to m k=1 j=k+1 for j = k + 1 to m R j,j:m = R j,j:m R k,j:m R kj /R kk R k,k:m = R k,k:m / R kk 工作量 Operation : m m count 2(m j + 1) O( m3 3 ) k=1 j=k+1 m m m k m 2(m j) 2 j k=1 j=1 k=1 k 2 m3 3
Cholesky 分解 稳定性 :Cholesky 分解得到的 R 满足 R R = A + δa, δa A = O(ɛ M) 即算法向后稳定
Cholesky 分解 稳定性 :Cholesky 分解得到的 R 满足 R R = A + δa, δa A = O(ɛ M) 即算法向后稳定 但是向前误差 R R / R = O(κ(A)ɛ M )
Cholesky 分解 稳定性 :Cholesky 分解得到的 R 满足 R R = A + δa, δa A = O(ɛ M) 即算法向后稳定 但是向前误差 R R / R = O(κ(A)ɛ M ) 求解 Ax = b (A 正定 ): 利用 Cholesky 分解 + 向前替换 + 向后替换算法 工作量 O( m3 3 ) 且向后稳定 (A + A) x = b, A A = O(ɛ M)