Ch5 解线性方程组的直接方法解线性方程组的直接方法求解求解 高斯消元法：高斯消元法：思思路路首先将首先将A化为上三角阵化为上三角阵 /* upper-triangular matrix */，再回代求解，再回代求解 /* backward substitution */。=记记Step 1：设设 ，计算因子，计算因子将增广矩阵将增广矩阵/* augmented matrix */ 第第 i 行行 mi1 第第1 1行行，得到得到其中其中Step k：设设 ，计算因子，计算因子且计算且计算共进行共进行 ？ 步步n 1What if ?No unique solution exists.What if ?Then we must find the smallest integer k i with , and interchange the k-th row with the i-th row. What if we cant find such k ?No unique solution exists.定理定理 若若A的所有的所有顺序主子式顺序主子式 /* determinant of leading principal submatrices */ 均不为均不为0，则高斯消元无需换行即可，则高斯消元无需换行即可进行到底，得到唯一解。进行到底，得到唯一解。注注注注：事实上，只要事实上，只要 A 非奇异，即非奇异，即 A 1 存在，则可通过逐次存在，则可通过逐次消元及行交换，将方程组化为三角形方程组，求出唯消元及行交换，将方程组化为三角形方程组，求出唯一解。一解。1 Gaussian Elimination The Method 选主元消去法选主元消去法例：例：单精度解方程组单精度解方程组/* 精确解为精确解为 和和 */8个个8个个用用Gaussian Elimination计算：计算：8个个小主元小主元 /* Small pivot element */ 可能导致计可能导致计算失败。算失败。 全主元消去法全主元消去法 /* Complete Pivoting */每一步选绝对值最大的元素为主元素，保证每一步选绝对值最大的元素为主元素，保证 。Step k: 选取选取 If ik k then 交换第交换第 k 行与第行与第 ik 行行; If jk k then 交换第交换第 k 列与第列与第 jk 列列; 消元消元注注注注：列交换改变了列交换改变了 xi 的顺序，须记录的顺序，须记录交换次序交换次序，解完后再，解完后再换回来。换回来。 列主元消去法列主元消去法 /* Partial Pivoting, or maximal column pivoting */省去换列的步骤，每次仅选一列中最大的元。省去换列的步骤，每次仅选一列中最大的元。例：例： 注注注注：列主元法没有全主元法稳定。列主元法没有全主元法稳定。例：例：注意：这两个方程组在注意：这两个方程组在数学上数学上严格等价严格等价。 标度化列主元消去法标度化列主元消去法 /* Scaled Partial Pivoting */对每一行计算对每一行计算 。为省时间，。为省时间，si 只在初始时计只在初始时计算一次。以后每一步考虑子列算一次。以后每一步考虑子列 中中 最大的最大的 aik 为主元。为主元。注注注注：稳定性介于列主元法和全主元法之间。稳定性介于列主元法和全主元法之间。1 Gaussian Elimination Pivoting Strategies2 三角分解法三角分解法 /* Matrix Factorization */ 高斯消元法的矩阵形式高斯消元法的矩阵形式 /* Matrix Form of G.E. */：Step 1:记记 L1 =，则，则Step n 1:其中其中 Lk =2 Matrix Factorization Matrix Form of G.E.记为记为L单位下三角阵单位下三角阵/* unitary lower-triangular matrix */记记 U =A 的的 LU 分解分解/* LU factorization */Hey hasnt GE given me enough headache? Why do I have to know its matrix form?!When you have to solve the system for different with a fixed A.Could you be more specific, please?Factorize A first, then for every you only have to solve two simple triangular systems and.2 Matrix Factorization Matrix Form of G.E.定理定理 若若A的所有顺序主子式的所有顺序主子式 /* determinant of leading principal submatrices */ 均不为均不为0，则，则 A 的的 LU 分解唯一（其分解唯一（其中中 L 为为单位单位下三角阵）。下三角阵）。证明：证明：由由1中定理可知，中定理可知，LU 分解存在。下面证明唯一性。分解存在。下面证明唯一性。若不唯一，则可设若不唯一，则可设 A = L1U1 = L2U2 ，推出，推出Upper-triangularLower-triangularWith diagonal entries 1注注注注： L 为一般下三角阵而为一般下三角阵而 U 为为单位单位上三角阵的分解称为上三角阵的分解称为Crout 分解分解。 实际上只要考虑实际上只要考虑 A* 的的 LU 分解，即分解，即 ，则，则 即是即是 A 的的 Crout 分解。分解。2 Matrix Factorization Doolittle 道立特分解法道立特分解法 /* Doolittle Factorization */： LU 分解的紧凑格式分解的紧凑格式 /* compact form */反复计算反复计算,很浪费哦很浪费哦 通过比较法直接导出通过比较法直接导出L 和和 U 的计算公式。的计算公式。思思路路2 Matrix Factorization Choleski 平方根法平方根法 /* Choleskis Method */： 对称对称 /* symmetric */ 正定正定 /* positive definite */ 矩阵的分解法矩阵的分解法一个矩阵一个矩阵 A = ( aij )n n 称为称为对称阵对称阵，如果，如果 aij = aji 。一个矩阵一个矩阵 A 称为称为正定阵正定阵，如果，如果 对任意非对任意非零向量零向量 都成立。都成立。回顾：回顾：对称正定阵的几个重要性质对称正定阵的几个重要性质 A 1 亦对称正定，且亦对称正定，且 aii 0若不然，则若不然，则存在非零解，即存在非零解，即存在非零解。存在非零解。 对任意对任意 , 存在存在 , 使得使得 ,即即 。 其中其中第第 i 位位 A 的顺序主子阵的顺序主子阵 /* leading principal submatrices */ Ak 亦对亦对 称正定称正定对称性显然。对任意对称性显然。对任意 有有 , 其中其中 。 A 的特征值的特征值 /* eigen value */ i 0 设对应特征值设对应特征值 的非零特征向量的非零特征向量为为 ，则，则 。 A 的全部顺序主子式的全部顺序主子式 det ( Ak ) 0因为因为2 Matrix Factorization Choleski将将对称对称 正定阵正定阵 A 做做 LU 分解分解U =uij=u11uij / uii111u22unn记为记为 A 对称对称即即记记 D1/2 =Why is uii 0?Since det(Ak) 0则则 仍是下三角阵仍是下三角阵定理定理 设矩阵设矩阵A对称正定，则存在非奇异下三角阵对称正定，则存在非奇异下三角阵 使得使得 。若限定。若限定 L 对角元为正，则分解唯一。对角元为正，则分解唯一。注注注注： 对于对称正定阵对于对称正定阵 A ，从，从 可知对任意可知对任意k i 有有 。即即 L 的元素不会增大，误差可控，的元素不会增大，误差可控，不需选主元。不需选主元。2 Matrix Factorization Choleski Algorithm: Choleskis MethodTo factor the symmetric positive definite n n matrix A into LLT, where L is lower triangular.Input: the dimension n; entries aij for 1 i, j n of A.Output: the entries lij for 1 j i and 1 i n of L. Step 1 Set ;Step 2 For j = 2, , n, set ;Step 3 For i = 2, , n 1, do steps 4 and 5Step 4 Set ; Step 5 For j = i+1, , n, set ; Step 6 Set ;Step 7 Output ( lij for j = 1, , i and i = 1, , n ); STOP. 因为因为A 对称，所以只需存半个对称，所以只需存半个 A，即，即其中其中运算量为运算量为 O(n3/6)， 比普通比普通LU分解少一半，但有分解少一半，但有 n 次开方。用次开方。用A = LDLT 分解，可省开方时间分解，可省开方时间(p.50-51)。HW: p.54 #2, #5, #62 Matrix Factorization Tridiagonal System 追赶法解追赶法解三对角三对角方程组方程组 /* Crout Reduction for Tridiagonal Linear System */Step 1: 对对 A 作作Crout 分解分解直接比较等式两边直接比较等式两边的元素，可得到计的元素，可得到计算公式算公式（p.52）。Step 2: 追追即解即解 ：Step 3: 赶赶即解即解 ：与与G.E.类似，一旦类似，一旦 i = 0 则算法中断，故并非任何则算法中断，故并非任何三对角阵都可以用此方法分解。三对角阵都可以用此方法分解。2 Matrix Factorization Tridiagonal System定理定理 若若 A 为为对角占优对角占优 /* diagonally dominant */ 的三对角的三对角阵，且满足阵，且满足 ，则，则追赶法可解以追赶法可解以 A 为系数矩阵的方程组。为系数矩阵的方程组。Hey, what does diagonally dominant mean? It means that the diagonal entries of the matrix are very LARGE.Well, how large is LARGE? They satisfy the following inequality:注注注注： 如果如果 A 是是严格严格对角占优阵，则不要求三对角线上的对角占优阵，则不要求三对角线上的所有元素非零。所有元素非零。 根据不等式根据不等式 可知：分解过程中，矩阵元素不会过分增大，算法可知：分解过程中，矩阵元素不会过分增大，算法保证稳定。保证稳定。 运算量为运算量为 O(6n)。课堂作业课堂作业1.矩阵范数矩阵范数误差分析误差分析矩阵的条件数矩阵的条件数精确解精确解 (2, 0)精确解精确解 (1, 1)下面两个定理描述了当系矩阵或是常数项下面两个定理描述了当系矩阵或是常数项有微小变化对解的影响有微小变化对解的影响矩阵的正交三角化及应用矩阵的正交三角化及应用初等反射阵的几何意义初等反射阵的几何意义wvyxv平面转转矩阵或吉文斯（平面转转矩阵或吉文斯（Givens）变换）变换 第第i列列第第j列列第第i行行第第j行行