Intro GE Triangular Errors Special Case Software Summary Solving triangular systems Solve Ly = b : b(1) = b(1)/L(1,1); for k=2 to n tmp = b(k) - L(k,1:k-1)*b(1:k-1); b(k) = tmp/L(k,k); end k = 3 1 0 0 7 L = b = , − 0 . 3 1 0 6 . 1 0 . 5 − 25 1 155 Solve Ux = b : Similar (backward).
Intro GE Triangular Errors Special Case Software Summary Operation counts LU decomposition: for k=1 to n-1 A(k+1:n,k) = A(k+1:n,k)/A(k,k); A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n); end � n − 1 ( n − k ) + 2 ( n − k ) 2 dk ≈ 2 3 n 3 . 1
Intro GE Triangular Errors Special Case Software Summary Operation counts Solving triangular systems: b(1) = b(1)/L(1,1); for k=2 to n tmp = b(k) - L(k,1:k-1)*b(1:k-1); b(k) = tmp/L(k,k); end � n 2 ( k − 1 ) dk ≈ 2 n 2 . 2 2
Intro GE Triangular Errors Special Case Software Summary Operation counts Question How long does it take to solve a system of order 30 on a computer that can perform two billion floating-point addition or multiplication operations (flops) per second? Answer: ( 0 . 7 × 30 3 + 2 × 30 2 ) / ( 2 × 10 9 ) ≈ 10 − 5 seconds!
Intro GE Triangular Errors Special Case Software Summary What is pivoting? Change the (2,2)-entry of A slightly from 2 to 2.099 and b 2 in b accordingly so that the exact solution is unchanged. x 1 10 − 7 0 7 x 2 = − 3 2 . 099 6 3 . 901 x 3 5 − 1 5 6 Exact solution: ( 0 , − 1 , 1 ) T .
Intro GE Triangular Errors Special Case Software Summary What is pivoting? (cont.) Forward elimination ( β = 10, p = 4) x 1 10 − 7 0 7 x 2 = 0 − 0 . 001 6 6 . 001 x 3 0 2 . 5 5 2 . 5 pivot: 10, multipliers: 0.3, − 0 . 5 x 1 10 − 7 0 7 x 2 = 0 − 0 . 001 6 6 . 001 x 3 1 . 501 × 10 4 1 . 500 × 10 4 0 0 pivot: − 0 . 001; multiplier: 2500.
Intro GE Triangular Errors Special Case Software Summary What is pivoting? (cont.) Backward substitution x 3 = 1 . 500 × 10 4 / 1 . 501 × 10 4 = 9 . 993 × 10 − 1 x 2 = ( 6 . 001 − 6 x 3 ) / ( − 0 . 001 ) = 5 . 0 What went wrong? Step 1: multipliers m 21 = 0 . 3, m 31 = − 0 . 5 1 . 000 × 10 1 − 7 . 000 0 7 . 000 − 1 . 000 × 10 − 3 0 6 . 000 6 . 001 2 . 500 0 + 2 . 500 5 . 000 Exact, no rounding errors.
Intro GE Triangular Errors Special Case Software Summary What is pivoting? (cont.) Step 2: multiplier m 31 = 2 . 500 E + 3, exact. 1 . 000 × 10 1 − 7 . 000 0 7 . 000 − 1 . 000 × 10 − 3 6 . 001 0 6 . 000 1 . 500 × 10 4 1 . 501 × 10 4 0 0 The rounding error in 1 . 501 × 10 4 ( A ( 3 , 3 ) , the exact result is 15 , 005) or 1 . 500 × 10 4 ( b 3 , the exact result is 15 , 005) equals half of its ulp ( ( 0 . 001 × 10 4 ) / 2 = 5), which has the same size as the size of the solution.
Intro GE Triangular Errors Special Case Software Summary What is pivoting? (cont.) Back solve: computed exact x 3 1 . 001 1.0 x 2 6 . 001 − 6 × 1 . 001 6 . 001 − 6 × 1 . 0 = − 5 . 000 = − 1 . 0 − 0 . 001 − 0 . 001 Catastrophic cancellation.
Intro GE Triangular Errors Special Case Software Summary What is pivoting? (cont.) Error in the result can be as large as half of the ulp of the largest intermediate results. Solution: Avoid large intermediate results (entries in the lower-right submatrices). Avoid small pivots (causing large multipliers). How? Interchange equations (rows).
Intro GE Triangular Errors Special Case Software Summary Pivoting Forward elimination step 2: x 1 10 − 7 0 7 x 2 = 0 − 0 . 001 6 6 . 001 x 3 0 2 . 5 5 2 . 5 Matrix form: 1 0 0 A ← M 1 A , M 1 = 0 . 3 1 0 − 0 . 5 0 1
Intro GE Triangular Errors Special Case Software Summary Pivoting (cont.) Interchange equations (rows) (2) and (3) to avoid small pivot, x 1 10 − 7 0 7 x 2 = 0 2 . 5 5 2 . 5 x 3 0 − 0 . 001 6 6 . 001 Matrix form: 1 0 0 A ← P 2 M 1 A , P 2 = 0 0 1 0 1 0 pivot: 2.5, multiplier: 4 × 10 − 4
Intro GE Triangular Errors Special Case Software Summary Pivoting (cont.) The updated system: x 1 10 − 7 0 7 x 2 = 0 2 . 5 5 2 . 5 x 3 0 0 6 . 002 6 . 002 Matrix form: 1 0 0 A ← M 2 P 2 M 1 A , M 2 = 0 1 0 4 × 10 − 4 0 1
Intro GE Triangular Errors Special Case Software Summary Pivoting (cont.) The updated system: x 1 10 − 7 0 7 x 2 = 0 2 . 5 5 2 . 5 x 3 0 0 6 . 002 6 . 002 Backward substitution: x 3 = 6 . 002 / 6 . 002 = 1 . 000 x 2 ( 2 . 5 − 5 x 3 ) / 2 . 5 = − 1 . 000 = x 1 ( 7 − ( − 7 ) x 2 − 0 x 3 ) / 10 = 0 =
Intro GE Triangular Errors Special Case Software Summary LU decomposition M 2 P 2 M 1 A = U , ( M − 1 1 P − 1 2 M − 1 2 ) U = A 1 0 0 1 0 0 M − 1 P 2 = P − 1 = − 0 . 3 1 0 = 0 0 1 1 2 0 . 5 0 1 0 1 0 1 0 0 M − 1 = 0 1 0 2 − 4 × 10 − 4 0 1 10 − 7 0 U = 0 2 . 5 5 0 0 6 . 002
Intro GE Triangular Errors Special Case Software Summary LU decomposition (cont.) But 1 0 0 M − 1 1 P 2 M − 1 − 4 × 10 − 4 = − 0 . 3 1 2 0 . 5 1 0 is not lower triangular. However, 1 0 0 ( P 2 M − 1 1 P 2 ) M − 1 = 0 . 5 1 0 2 − 4 × 10 − 4 − 0 . 3 1 is lower triangular and elementary! Call it L . ( P 2 M − 1 1 P 2 M − 1 2 ) U = P 2 A is an LU decomposition of P 2 A , (row) permuted A .
Intro GE Triangular Errors Special Case Software Summary LU decomposition (cont.) Consider 1 0 0 M − 1 ˆ := P 2 M − 1 1 P 2 = 0 . 5 1 0 1 − 0 . 3 0 1 equivalent to interchanging m 21 and m 31 . In general, suppose that M 2 P 2 M 1 P 1 A = U , that is, A = P 1 M − 1 1 P 2 M − 1 2 U , then P 2 P 1 A = (( P 2 M − 1 1 P 2 ) M − 1 2 ) U An LU decomposition of permuted A .
Intro GE Triangular Errors Special Case Software Summary LU decomposition (cont.) In n -dimensional case, M n − 1 P n − 1 · · · M 2 P 2 M 1 P 1 A = U . Then A P 1 M − 1 1 P 2 M − 1 · · · P n − 2 M − 1 n − 2 P n − 1 M − 1 n − 1 U = 2 P 1 M − 1 1 P 2 M − 1 · · · P n − 3 M − 1 n − 3 P n − 2 P n − 1 = 2 ( P n − 1 M − 1 n − 2 P n − 1 ) M − 1 n − 1 U P 1 ... P n − 1 ( P n − 1 ... P 2 M − 1 1 P 2 ... P n − 1 ) · · · = ( P n − 1 M − 1 n − 2 P n − 1 ) M − 1 n − 1 U LU decomposition P n − 1 ... P 1 A = LU of permuted A . In programming, A is overwritten by L and U and L is stored in the lower part of A . When we interchange rows of A , we also interchange corresponding (entire) rows of L .
Intro GE Triangular Errors Special Case Software Summary Pivoting (cont.) Algorithm. Gaussian elimination with partial pivoting. for k=1 to n-1 find the pivot: max(|A(k:n,k)|); record the pivoting row index in p(k); interchange rows A(k,:) and A(p(k),:); A(k+1:n,k) = A(k+1:n,k)/A(k,k); A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n); end
Intro GE Triangular Errors Special Case Software Summary Example Original − 3 2 . 099 6 1 p = 10 − 7 0 2 5 − 1 5 1
Intro GE Triangular Errors Special Case Software Summary Example k = 1, pivot − 3 2 . 099 6 2 p = 10 − 7 0 2 5 − 1 5 1
Intro GE Triangular Errors Special Case Software Summary Example k = 1, permute P 1 10 − 7 0 2 p = − 3 2 . 099 6 2 5 − 1 5 − 1
Intro GE Triangular Errors Special Case Software Summary Example k = 1, Multipliers M 1 10 − 7 0 2 p = − 0 . 3 2 . 099 6 2 0 . 5 − 1 5 − 1
Intro GE Triangular Errors Special Case Software Summary Example k = 1, update 10 − 7 0 2 p = − 0 . 3 − 0 . 001 6 2 0 . 5 2 . 5 5 − 1
Intro GE Triangular Errors Special Case Software Summary Example k = 2, pivot 10 − 7 0 2 p = − 0 . 3 − 0 . 001 6 3 0 . 5 2 . 5 5 − 1
Intro GE Triangular Errors Special Case Software Summary Example k = 2, permute P 2 10 − 7 0 2 p = 0 . 5 2 . 5 5 3 − 0 . 3 − 0 . 001 6 1
Intro GE Triangular Errors Special Case Software Summary Example k = 2, multiplier M 2 10 − 7 0 2 p = 0 . 5 2 . 5 5 3 − 0 . 3 − 0 . 0004 6 1
Intro GE Triangular Errors Special Case Software Summary Example k = 2, update 10 − 7 0 2 p = 0 . 5 2 . 5 5 3 − 0 . 3 − 0 . 0004 6 . 002 1
Intro GE Triangular Errors Special Case Software Summary Example 1 0 0 10 − 7 0 LU = 0 . 5 1 0 0 2 . 5 5 − 0 . 3 − 0 . 0004 1 0 0 6 . 002 10 − 7 0 , = 5 − 1 5 − 3 2 . 099 6 permuted A − 3 2 . 099 6 2 P 2 P 1 p = 10 − 7 0 3 5 − 1 5 1 det ( A ) = p ( n ) U ( 1 , 1 ) ... U ( n , n ) = 150 . 05
Intro GE Triangular Errors Special Case Software Summary Basic Linear Algebra Subroutines (BLAS) Operations in GE with pivoting: imax: find the index of max ( | A ( k : n , k ) | ) swap: interchange A ( k , :) and A ( p ( k ) , :) scal: scalar-vector multiplication: A ( k , k ) − 1 A ( k + 1 : n , k ) axpy: α x + y , vector update: A ( k + 1 : n , j ) − A ( k , j ) A ( k + 1 : n , k ) , for j = k + 1 : n
Intro GE Triangular Errors Special Case Software Summary Problem setting Given the decomposition LU = PA , solve Ax = b or solve PAx = Pb Solve two triangular systems: stage 1, solve Ly = Pb (forward solve) stage 2, solve Ux = y (back solve) LUx = Pb Consider Ux = y , given U and y , the solution x overwrites y .
Intro GE Triangular Errors Special Case Software Summary Two versions Row version (C): y ( n ) = y ( n ) / U ( n , n ) ; for i = n − 1 : − 1 : 1; y ( i ) = y ( i ) − U ( i , i + 1 : n ) y ( i + 1 : n ) ; y ( i ) = y ( i ) / U ( i , i ) ; end for i
Intro GE Triangular Errors Special Case Software Summary Column version (FORTRAN) y ( n ) = y ( n ) / U ( n , n ) ; for i = n − 1 : − 1 : 1; y ( 1 : i ) = y ( 1 : i ) − y ( i + 1 ) U ( 1 : i , i + 1 ) ; y ( i ) = y ( i ) / U ( i , i ) ; end for i 10 − 7 0 7 0 2 . 5 5 2 . 5 0 0 6 . 002 6 . 002
Intro GE Triangular Errors Special Case Software Summary Column version (FORTRAN) y ( n ) = y ( n ) / U ( n , n ) ; for i = n − 1 : − 1 : 1; y ( 1 : i ) = y ( 1 : i ) − y ( i + 1 ) U ( 1 : i , i + 1 ) ; y ( i ) = y ( i ) / U ( i , i ) ; end for i 10 − 7 0 7 0 2 . 5 5 2 . 5 0 0 6 . 002 1 . 0
Intro GE Triangular Errors Special Case Software Summary Column version (FORTRAN) y ( n ) = y ( n ) / U ( n , n ) ; for i = n − 1 : − 1 : 1; y ( 1 : i ) = y ( 1 : i ) − y ( i + 1 ) U ( 1 : i , i + 1 ) ; y ( i ) = y ( i ) / U ( i , i ) ; end for i 10 − 7 0 7 0 2 . 5 5 − 2 . 5 0 0 6 . 002 1 . 0
Intro GE Triangular Errors Special Case Software Summary Column version (FORTRAN) y ( n ) = y ( n ) / U ( n , n ) ; for i = n − 1 : − 1 : 1; y ( 1 : i ) = y ( 1 : i ) − y ( i + 1 ) U ( 1 : i , i + 1 ) ; y ( i ) = y ( i ) / U ( i , i ) ; end for i 10 − 7 0 7 0 2 . 5 5 − 1 . 0 0 0 6 . 002 1 . 0
Intro GE Triangular Errors Special Case Software Summary Column version (FORTRAN) y ( n ) = y ( n ) / U ( n , n ) ; for i = n − 1 : − 1 : 1; y ( 1 : i ) = y ( 1 : i ) − y ( i + 1 ) U ( 1 : i , i + 1 ) ; y ( i ) = y ( i ) / U ( i , i ) ; end for i 10 − 7 0 0 0 2 . 5 5 − 1 . 0 0 0 6 . 002 1 . 0
Intro GE Triangular Errors Special Case Software Summary Column version (FORTRAN) y ( n ) = y ( n ) / U ( n , n ) ; for i = n − 1 : − 1 : 1; y ( 1 : i ) = y ( 1 : i ) − y ( i + 1 ) U ( 1 : i , i + 1 ) ; y ( i ) = y ( i ) / U ( i , i ) ; end for i 10 − 7 0 0 0 2 . 5 5 − 1 . 0 0 0 6 . 002 1 . 0
Intro GE Triangular Errors Special Case Software Summary Introduction Estimating the error in the computed solution. x Computed solution: ˆ Exact solution: x Relative forward error: � x − ˆ x � / � ˆ x � But we usually don’t know x . Check the residual r = b − A ˆ x ? A contrived example. β = 10, p = 3 � � x 1 � 1 . 15 � 2 . 15 � � 1 . 00 = x 2 1 . 41 1 . 22 2 . 63
Intro GE Triangular Errors Special Case Software Summary Introduction (cont.) Gaussian elimination with partial pivoting Interchange two rows; multiplier: 1 . 15 / 1 . 41 = 0 . 816; � � x 1 � 1 . 41 � 2 . 63 � � 1 . 22 = x 2 0 0 . 004 0 . 00 x 2 = 0, ˆ x 1 = 2 . 63 / 1 . 41 = 1 . 87. ˆ Residual: r 1 = 0 . 01, r 2 = 0. Exact solution: x 1 = x 2 = 1. Small residual does not imply small error in solution.
Intro GE Triangular Errors Special Case Software Summary Introduction (cont.) Check the pivots? If the pivots are small, A is nearly singular; however, A might be nearly singular but none of the pivots are small. How do we estimate the error in the computed solution? We know that the relative forward error: � x − ˆ x � ≤ cond · backward error . x � � ˆ
Intro GE Triangular Errors Special Case Software Summary Estimating backward error x is the exact solution of Suppose that the computed solution ˆ the perturbed system: ( A + ∆ A )ˆ x = b , then the relative backward error is � ∆ A � � A � . How do we get � ∆ A � ?
Intro GE Triangular Errors Special Case Software Summary Estimating backward error Since � b − A ˆ x � = � ∆ A ˆ x � x � x � ≤ � ∆ A � � A � , � A � � ˆ � A � � ˆ in practice, the relative residual � b − A ˆ x � � A � � ˆ x � can be used as an estimate for the backward error. In other words, the relative residual is an indication of the stability of the algorithm.
Intro GE Triangular Errors Special Case Software Summary Condition of a matrix The sensitivity of the solution x to the perturbations on A and b . Measure of nearness of singularity. Vector norms: � x � > 0, if x � = 0 � 0 � = 0 � cx � = | c | � x � , for all scalars c � x + y � ≤ � x � + � y � Examples. � x � 2 = i | x i | 2 � 1 / 2 � � � x � 1 = � i | x 1 | � x � ∞ = max i ( | x i | )
Intro GE Triangular Errors Special Case Software Summary Condition of a matrix (cont.) Think of a matrix as a linear transformation between two vector spaces. The range of possible change: � Ax � M � x � = � A � = max x � = 0 � Ax � m = min � x � x � = 0 When A is singular, m = 0. Examples of matrix norms. � A � 1 = max j i | a ij | � � � � A � ∞ = max i j | a ij | � � �
Intro GE Triangular Errors Special Case Software Summary Condition of a matrix Measurement: cond ( A ) = M / m If A is singular, m = 0, cond ( A ) = ∞ . If A is nonsingular, � x � � A − 1 y � m − 1 = max = � A − 1 � � Ax � = max � y � x � = 0 y � = 0 Condition for inverting A cond ( A ) = � A � � A − 1 �
Intro GE Triangular Errors Special Case Software Summary Condition of a matrix (cont.) Perturbing b : A ( x + ∆ x ) = b + ∆ b Since Ax = b and A (∆ x ) = ∆ b , � b � ≤ M � x � � ∆ b � ≥ m � ∆ x � and Thus � ∆ x � ≤ cond ( A ) � ∆ b � � x � � b � A relative error magnification factor. How do we get � A − 1 � ?
Intro GE Triangular Errors Special Case Software Summary Estimating � A − 1 � 1 Basic idea: Determine a vector e (all components 1 or − 1) so that the solution for A T Az = e is heuristically large. Given PA = LU ( A T = U T L T P and A T A = U T L T LU ), determine e so that the solution w for U T w = e is heuristically large; solve for y in L T y = w ; solve for z in Az = y ; normalize � z � 1 / � y � 1 ≈ � A − 1 � 1 . Cost: If the LU decomposition PA = LU ( O ( n 3 ) ) is available, it requires solving four triangular systems ( O ( n 2 ) ).
Intro GE Triangular Errors Special Case Software Summary Example revisited � 1 . 15 � 2 . 15 � � 1 . 00 A = b = , 1 . 41 1 . 22 2 . 63 Exact solution: x = [ 1 1 ] T x = [ 1 . 87 0 . 00 ] T , β = 10, t = 3. Computed solution: ˆ
Intro GE Triangular Errors Special Case Software Summary Example revisited The computed solution is the exact solution of � 215 / 187 � 1 . 00 A = b ˆ , and 263 / 187 1 . 22 The perturbation � 2 . 67 × 10 − 4 � 0 ∆ A = A − ˆ A ≈ 3 . 58 × 10 − 3 0 Relative change in A : � ∆ A � � A � = 1 . 5 × 10 − 3 := ρ u
Intro GE Triangular Errors Special Case Software Summary Example revisited Relative error in the computed solution: � x − ˆ x � ≈ 0 . 7088 x � � ˆ Almost 100 % error, that is, zero digit accuracy. The condition number cond ( A ) ≈ 8 . 26 × 10 2 Almost u − 1 = β t . � x − ˆ x � ≤ ρ cond ( A ) u x � � ˆ The relative change in A is magnified by cond ( A ) .
Intro GE Triangular Errors Special Case Software Summary Remarks The solution computed by GE with partial pivoting can be viewed as the exact solution of a slightly perturbed coefficient matrix. The relative perturbation is usually ρ u , where ρ is a constant of the same size as β . In other words, in practice, GE with partial pivoting is stable. The relative error in the computed solution (by GE with partial pivoting) is roughly of the size cond ( A ) u . In practice, the entries in the coefficient matrix A and the right-hand-side vector b contain measurement errors. In the computed solution, The measurement error is roughly magnified by cond ( A ) . For example, if the measurement accuracy is four decimal digits and the condition number is about 10 2 , then we expect the computed solution has two decimal digit accuracy.
Intro GE Triangular Errors Special Case Software Summary Symmetric positive definite systems Symmetric: A T = A , or a ij = a ji Positive definite: x T Ax > 0, for any x � = 0 Properties: nonsingular all principal submatrices are symmetric and positive definite The method of choice: Cholesky factorization A = LL T L : lower triangular
Intro GE Triangular Errors Special Case Software Summary Cholesky factorization Idea: An equation a 11 a 21 a 31 l 11 l 11 l 21 l 31 0 0 a 21 a 22 a 32 l 21 l 22 l 22 l 32 = 0 0 a 31 a 32 a 33 l 31 l 32 l 33 l 33 0 0
Intro GE Triangular Errors Special Case Software Summary Cholesky factorization a 11 = l 11 l 11 a 11 a 21 a 31 l 11 l 11 l 21 l 31 0 0 a 21 a 22 a 32 l 21 l 22 l 22 l 32 = 0 0 a 31 a 32 a 33 l 31 l 32 l 33 l 33 0 0
Intro GE Triangular Errors Special Case Software Summary Cholesky factorization a 21 = l 21 l 11 a 11 a 21 a 31 l 11 l 11 l 21 l 31 0 0 a 21 a 22 a 32 l 21 l 22 l 22 l 32 = 0 0 a 31 a 32 a 33 l 31 l 32 l 33 l 33 0 0
Intro GE Triangular Errors Special Case Software Summary Cholesky factorization a 22 = l 21 l 21 + l 22 l 22 a 11 a 21 a 31 l 11 l 11 l 21 l 31 0 0 a 21 a 22 a 32 l 21 l 22 l 22 l 32 = 0 0 a 31 a 32 a 33 l 31 l 32 l 33 l 33 0 0
Intro GE Triangular Errors Special Case Software Summary Cholesky factorization a 31 = l 31 l 11 a 11 a 21 a 31 l 11 l 11 l 21 l 31 0 0 a 21 a 22 a 32 l 21 l 22 l 22 l 32 = 0 0 a 31 a 32 a 33 l 31 l 32 l 33 l 33 0 0
Intro GE Triangular Errors Special Case Software Summary Cholesky factorization a 32 = l 31 l 21 + l 32 l 22 a 11 a 21 a 31 l 11 l 11 l 21 l 31 0 0 a 21 a 22 a 32 l 21 l 22 l 22 l 32 = 0 0 a 31 a 32 a 33 l 31 l 32 l 33 l 33 0 0
Intro GE Triangular Errors Special Case Software Summary Cholesky factorization a 33 = l 31 l 31 + l 32 l 32 + l 33 l 33 a 11 a 21 a 31 l 11 l 11 l 21 l 31 0 0 a 21 a 22 a 32 l 21 l 22 l 22 l 32 = 0 0 a 31 a 32 a 33 l 31 l 32 l 33 l 33 0 0
Intro GE Triangular Errors Special Case Software Summary Cholesky factorization In general, a 11 a 12 a 1 n · · · a 12 a 22 a 2 n · · · . . . . . . . . . . . . a 1 n a 2 n a nn · · · l 11 l 11 l 21 l 1 n 0 · · · l 21 l 22 l 22 l 2 n · · · = . . . . ... ... . . . . . . l n 1 l n 2 l nn l nn · · ·
Intro GE Triangular Errors Special Case Software Summary Cholesky factorization (cont.) Consider the first row of L : l 11 = √ a 11 . Now we consider the general case: the i th row. Suppose we have computed the rows 1 through i − 1 of L , then j − 1 l ij = ( a ji − l jk l ik ) / l jj j = 1 : i − 1 � k = 1 � i − 1 � l ii = � � a ii − l 2 � ik k = 1
Intro GE Triangular Errors Special Case Software Summary Algorithm. Cholesky factorization (scalar) This algorithm computes a lower triangular L from an n -by- n symmetric and positive definite A so that A = LL T . for i=1 to n for j=1 to i % compute L(i,1:i) s = A(j,i); for k=1 to j-1 s = s - L(j,k)*L(i,k); end if j<i L(i,j) = s/L(j,j); else % j=i L(i,i) = sqrt(s); end end end
Intro GE Triangular Errors Special Case Software Summary Example 25 10 10 A = 10 53 32 10 32 36 Initial 25 L = 10 53 10 32 36
Intro GE Triangular Errors Special Case Software Summary Example (cont.) for i=1 to n for j=1 to i L(i,j) = L(i,j) - L(j,1:j-1)’*L(i,1:j-1); if j < i L(i,j) = L(i,j)/L(j,j); else L(i,i) = sqrt(L(i,i)); end end end i = 1, j = 1 5 L = 10 53 10 32 36
Intro GE Triangular Errors Special Case Software Summary Example (cont.) for i=1 to n for j=1 to i L(i,j) = L(i,j) - L(j,1:j-1)’*L(i,1:j-1); if j < i L(i,j) = L(i,j)/L(j,j); else L(i,i) = sqrt(L(i,i)); end end end i = 2, j = 1 5 L = 2 53 10 32 36
Intro GE Triangular Errors Special Case Software Summary Example (cont.) for i=1 to n for j=1 to i L(i,j) = L(i,j) - L(j,1:j-1)’*L(i,1:j-1); if j < i L(i,j) = L(i,j)/L(j,j); else L(i,i) = sqrt(L(i,i)); end end end i = 2, j = 2 5 L = 2 7 10 32 36
Intro GE Triangular Errors Special Case Software Summary Example (cont.) for i=1 to n for j=1 to i L(i,j) = L(i,j) - L(j,1:j-1)’*L(i,1:j-1); if j < i L(i,j) = L(i,j)/L(j,j); else L(i,i) = sqrt(L(i,i)); end end end i = 3, j = 1 5 L = 2 7 2 32 36
Intro GE Triangular Errors Special Case Software Summary Example (cont.) for i=1 to n for j=1 to i L(i,j) = L(i,j) - L(j,1:j-1)’*L(i,1:j-1); if j < i L(i,j) = L(i,j)/L(j,j); else L(i,i) = sqrt(L(i,i)); end end end i = 3, j = 2 2 3 5 L = 2 7 4 5 2 4 36
Intro GE Triangular Errors Special Case Software Summary Example (cont.) for i=1 to n for j=1 to i L(i,j) = L(i,j) - L(j,1:j-1)’*L(i,1:j-1); if j < i L(i,j) = L(i,j)/L(j,j); else L(i,i) = sqrt(L(i,i)); end end end i = 3, j = 3 (final) 2 3 5 L = 2 7 4 5 2 4 4
Intro GE Triangular Errors Special Case Software Summary Cholesky factorization storage n ( n + 1 ) / 2 flops n 3 / 3 + O ( n 2 ) stable the cheapest way to check if a symmetric matrix is positive definite
Intro GE Triangular Errors Special Case Software Summary Software packages Direct methods for general linear systems NETLIB LAPACK: sgetrf, sgetrs, sgecon Direct methods for symmetric and positive definite systems NETLIB LAPACK: spotrf, spotrs Direct methods for symmetric and indefinite systems NETLIB LAPACK: ssytrf, ssytrs Direct methods for sparse systems NETLIB SuperLU, SPARSE MATLAB colmmd, symmmd, symrcm
Intro GE Triangular Errors Special Case Software Summary Summary Gaussian elimination with pivoting ( decomp , solve ): Working on one matrix, matrix update. Improve instability by avoiding small pivots (controlling the sizes of intermediate results) Error estimates: Condition number of a matrix. Special systems: Triangle, symmetric and positive definite, Cholesky factorization.
Recommend
More recommend
