Efficiency Issues: 3/4 Efficiency Issues: 3/4 � The following can be proved with induction: The following can be proved with induction: 1 + + = + � n 1 2 ( 1 ) n n 2 ( ) 1 + + + + = + + 2 2 2 � n 2 3 2 1 2 3 6 2 3 n n n � Therefore, we have (verify yourself): − − − ( ( ) ) 1 1 1 n n n 1 ∑ ∑ ∑ ∑ ∑ ∑ − − + = − + − = − 2 2 3 3 ( ( )( )( 1 1 ) ) ( ( ) ) ( ( ) ) n i n i i i n i i n i i n n 3 = = − 1 1 1 i i i � This is an O ( n 3 ) algorithm ( i.e ., number of � Thi i 3 ) l O ( ith ( i b f multiplications proportional to n 3 ). 16
Efficiency Issues: 4/4 Efficiency Issues: 4/4 � In the backward substitution step, since we need to In the backward substitution step, since we need to compute the sum a i , i +1 x i +1 + a i , i +2 x i +2 +…+ a i , n x n , n - i multiplications are needed. multiplications are needed. � Since i runs from n -1 down to 1, the total number − n 1 1 ∑ ∑ − = = − of multiplications is of multiplications is ( ( ) ) ( ( 1 1 ) ) n n i i n n n n 2 = 1 i � In summary, the total number of multiplications used is dominated by the elimination step, which d i d i t d b th li i ti t hi h is proportional to n 3 /3, or O ( n 3 ) ! � Exercise: How many divisions are needed in the elimination and the backward substitution stages? 17
A Serious Problem A Serious Problem � There is a big � There is a big big problem When eliminating a big problem. When eliminating a k , i using a i , i , – a k , i / a i , i is computed. � What if a i , i ≈ 0? Division-by-zero or overflow � Wh t if 0? Di i i b fl in – a k , i / a i , i or truncation in (– a k , i / a i , i ) × a i , j + a k , j may occur because of (– a k , i / a i , i ) × a i,j >> a k,j ! Wh b f ( / ) ! Wh Wh ? Why? � Therefore, the use of Gaussian elimination is risky and some modifications would be needed. � However, if the following holds ( i.e ., diagonal , g ( , g dominant ), Gaussian elimination works fine. > ∑ n ∑ > | | | | | | | | a a a a (f (for every i ) i ) , , i i i j = ≠ 1 , 18 j j i
Partial Pivoting: 1/6 Partial Pivoting: 1/6 � To overcome the possible x /0 issue, one may use p , y an important technique: pivoting . � There are two types of pivoting, partial (or row ) yp p g, p ( ) and complete . For most cases, partial pivoting usually works efficiently and accurately. � Important Observation Important Observation: Swapping two equations does not change the solutions! � If a i , i ≈ 0, one may swap some equation k with equation i so that the new a i , i would not be zero. , � But, which k should be swapped? 19
Partial Pivoting: 2/6 Partial Pivoting: 2/6 � One of the best candidates is the a k , i such that | a k , i | i th is the maximum of | a i , i |, | a i +1, i |, …, | a n , i |. i f | | | | | | � Why? Why? If equation k with maximum | a k , i | is swapped to equation i , then, all | a j , i / a i , i | ≤ 1 for i ≤ j | ≤ 1 f i ≤ j d t ti i th ll | / ≤ n . � � a a a a , + 11 1 , 1 1 1 i i n � � � � 0 a a a + , , 1 , i i i i i n � 0 a a a + + + + + + + + 1 1 , 1 1 , 1 1 1 1 , i i i i i i i i i i n n find the maximum here � � � � followed by a swap � � 0 0 a a a a a a + , , 1 1 , n i n i n n 20
Partial Pivoting: 3/6 Partial Pivoting: 3/6 � The original Gaussian elimination is modified to i include pivoting: find the maximum | a k , i | of a i , i , l d i ti fi d th i | | f a i +1, i , …, a n , i , and do a row swap, including b i and b k . The remaining is the same! The remaining is the same! b DO i = 1, n-1 ! going through rows i+1 to n Max = i ! assume |a(i,i)| is the max DO k = i+1, n ! find the largest in Max IF (ABS(a(Max,i)) < ABS(a(k,i)) Max = k END DO DO j = i, n ! swap row Max and row i swap a(Max,j) and a(i,j) END DO O swap b(Max) and b(i) ! don’t forget to swap b(i) … do the elimination step … END DO END DO 21
Partial Pivoting: 4/6 Partial Pivoting: 4/6 � Blue dashed line: current column, Red circle: max A b A b 1 -2 3 1 3 3 -2 1 5 7 -2 1 -2 -1 -4 0 -1/3 -4/3 7/3 2/3 3 -2 1 5 7 0 -4/3 8/3 -2/3 2/3 1 -1 5 3 8 1 1 5 3 8 0 -1/3 14/3 4/3 17/3 0 1/3 14/3 4/3 17/3 row swap elimination elimination A b × (2/3)+ × (2/3)+ × (-1/3)+ × ( 1/3)+ × (-1/3)+ × ( 1/3)+ 3 -2 3 -2 1 5 7 1 5 7 -2 1 -2 -1 -4 1 -2 3 1 3 1 -1 5 3 8 22
Partial Pivoting: 5/6 Partial Pivoting: 5/6 � Blue dashed line: current column, Red circle: max A b A b 3 -2 1 5 7 3 -2 1 5 7 0 0 -1/3 -4/3 7/3 2/3 1/3 4/3 7/3 2/3 0 0 -4/3 8/3 -2/3 2/3 4/3 8/3 2/3 2/3 0 -4/3 8/3 -2/3 2/3 0 0 -2 5/2 1/2 0 -1/3 14/3 4/3 0 -1/3 14/3 4/3 17/3 17/3 0 0 0 0 4 4 3/2 11/2 3/2 11/2 row swap row swap elimination elimination A b 3 -2 3 -2 1 5 7 1 5 7 × (-1/4)+ × (-1/4)+ 0 -4/3 8/3 -2/3 2/3 0 -1/3 -4/3 7/3 2/3 0 -1/3 14/3 4/3 17/3 23
Partial Pivoting: 6/6 Partial Pivoting: 6/6 � Blue dashed line: current column, Red circle: max A b A b 3 -2 1 5 7 3 -2 1 5 7 0 -4/3 8/3 -2/3 2/3 0 -4/3 8/3 -2/3 2/3 0 0 -2 5/2 1/2 0 0 4 3/2 11/2 0 0 4 3/2 11/2 0 0 0 13/4 13/4 x 1 = x 2 = x 3 = x 4 = 1 elimination elimination row swap A b 3 3 -2 1 5 7 2 1 5 7 0 -4/3 8/3 -2/3 2/3 11/2 × (1/2)+ ( ) 0 0 0 0 4 4 3/2 11/2 3/2 0 0 -2 5/2 1/2 24
Complete Pivoting: 1/4 Complete Pivoting: 1/4 � What if a i , i , a i +1, i , …, a n , i are all very close to zero , , , when doing pivoting for column i ? � If this happens, partial pivoting has a problem, pp , p p g p , because no matter which row is swapped, there is a chance to have 0/0, overflow or , truncation/rounding. � In this case, complete pivoting may help. � In this case, complete pivoting may help. 25
Complete Pivoting: 2/4 Complete Pivoting: 2/4 � With complete pivoting, the block defined by a i , i and a n , n is searched for a maximum | a p , q |. d i h d f i | | � Then, a row swap (row i and row p ) and a column swap (column i and column q ) are required. ( l d l ) i d i � After swapping, do the usual elimination. � � a a a a , + 11 1 , 1 1 1 i i n � � � row i � 0 a a a + , , 1 , i i i i i n � 0 0 a a a a a a + + + + 1 1 , 1 1 , 1 1 1 1 , i i i i i n row p � � � � � � 0 0 a a a a a a + , , 1 , n i n i n n 26 column i column q
Complete Pivoting: 3/4 Complete Pivoting: 3/4 � While swapping rows does not change the While swapping rows does not change the solutions, swapping column i and column q changes the positions of x i and x q . changes the positions of x i and x q . � To overcome this problem, we have to keep track the swapping operations so that column swapping the swapping operations so that column swapping will not affect the solutions. � O � One may use an index array ! ! i d 27
Complete Pivoting: 4/4 Complete Pivoting: 4/4 � Index array idx() is initialized to 1, 2, …, n . y , , , � If columns i and q are swapped, the content of the i -th and q -th entries of idx() are also q swapped ( i.e ., swapping idx(i) and idx(q) ). � At the end, idx(k) = h means column k is , the solution of x h . For example, if idx(1)=4 , idx(2)=3 , idx(3)=1 and idx(4)=2 , this means columns 1, 2, 3 and 4 contain the solutions to x 4 , x 3 , x 1 and x 2 . � Sometimes, this index array is referred to as a permutation array . 28
Efficiency Concerns Efficiency Concerns � Elimination with pivoting does not increase the p g number of multiplications; however, it does use CPU time for comparisons and swapping. � Although compared with multiplications and divisions swapping may be insignificant, pivoting does add speed penalty to the efficiency of the methods. � One may use index arrays to avoid actually carrying out swapping. � Exercise: how can this be done? Refer to the index array method for complete swapping. 29
Is Pivoting Really Necessary? 1/3 Is Pivoting Really Necessary? 1/3 � Consider the following without pivoting. × (-1/0.0003)+ 0.0003 x + 3.000 y = 2.0001 1 0000 1.0000 x + 1.000 y = 1.0000 1 000 1 0000 elimination 0.0003 x + 3.000 y = 2.0001 exact arithmetic yields x = 1/3 and y = 2/3 -9999 y = -6666 y = 6666 9999 / − × 2 0001 . 30000 . ( 6666 9999 / ) x = 00003 00003 . 30
Is Pivoting Really Necessary? 2/3 Is Pivoting Really Necessary? 2/3 � Without pivoting: Without pivoting: Possible cancellation here, as 3 × (6666/9999) ≈ 2.0001! y = 6666 9999 / − × 2 0001 . 30000 . ( 6666 9999 / ) x = 00003 00003 . Precision y x 4 0.6667 0 4 0 6667 0 5 0.66667 0.3 6 0.666667 0.33 inaccurate inaccurate 7 0.6666667 0.333 31
Is Pivoting Really Necessary? 3/3 Is Pivoting Really Necessary? 3/3 � With pivoting: With pivoting: elimination 1.0000 x + 1.0000 y = 1.0000 1.0000 x + 1.0000 y = 1.0000 0 0003 x + 3 0000 y = 2 0001 0.0003 x + 3.0000 y = 2.0001 2 9997 y = 1 9998 2.9997 y = 1.9998 y = 19998 . y Backward substitution Backward substitution 2 9997 . = − 10000 . x y Precision y x 4 4 0.6667=2/3 0.3333 0 6667=2/3 0 3333 5 0.66667 0.33333 6 0.666667 0.333333 7 0.6666667 0.3333333 32
Pitfalls of Elimination: 1/2 Pitfalls of Elimination: 1/2 � The first one is, of course, a i i ≈ 0 when , , i , i computing – a k , i / a i , i . This can be overcome with partial or complete pivoting. � However, singular systems cannot be solved by elimination ( e.g ., two parallel lines do not have a solution). � Rounding can be a problem. Even with pivoting, rounding is still there and could propagate from earlier stages to the next. � In general, n ≤ 100 is OK. Otherwise, consider using other ( e.g ., indirect) methods. 33
Pitfalls of Elimination: 2/2 Pitfalls of Elimination: 2/2 � Ill-Conditioned systems are trouble makers. Ill Conditioned systems are trouble makers. � Ill-Conditioned systems are systems very sensitive to small changes to coefficients and sensitive to small changes to coefficients, and there could be many seemingly correct “solutions.” � Since these “solutions” seem to satisfy the � Si th “ l ti ” t ti f th systems, one may be misled to believe they are “ “correct” solutions. t” l ti 34
LU-Decompositions: 1/8 LU Decompositions: 1/8 � Idea: Idea: The basic idea of LU-decomposition is to Idea: Idea: The basic idea of LU decomposition is to “decompose” the given matrix A=[ a i , j ] of A • x=b into a product of a lower triangular and an into a product of a lower triangular and an upper triangular matrices ( i.e ., A = L • U). � The lower triangular matrix has all diagonal � The lower triangular matrix has all diagonal elements being 1’s ( i.e ., Doolittle form). lower triangular lower triangular upper triangular upper triangular ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ � � � 1 0 0 a a a u u u 1,1 1,2 1, n 1,1 1,2 1, n ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ � � � 1 0 0 a a a l u u ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ 2,1 2 1 2 2 2,2 2, 2 = 2 1 2,1 • 2 2 2,2 2, 2 n n n n ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ � � � � � � � � � � � � ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ � � � 1 0 0 a a a l l u ⎣ ⎣ ⎦ ⎦ ⎣ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎦ ,1 1 ,2 2 , ,1 1 ,2 2 , n n n n n n n n n n n n n n n n A = L • U 35
LU Decompositions: 2/8 LU-Decompositions: 2/8 � If A has been decomposed as A=L ⋅ U, then If A has been decomposed as A L U, then solving A ⋅ x = b becomes solving (L ⋅ U) ⋅ x = b. � (L ⋅ U) ⋅ x = b can be rewritten as L ⋅ (U ⋅ x) = b � (L ⋅ U) ⋅ x = b can be rewritten as L ⋅ (U ⋅ x) = b. � Let y = U ⋅ x. Then, L ⋅ (U ⋅ x) = b becomes tw o tw o systems of linear equations: L ⋅ y = b and U ⋅ x = y. t f li ti L b d U � Since L and b are known, one can solve for y. Once y becomes available, it is used in U ⋅ x = y to solve for x. � This is the key concept of LU-decomposition. 36
LU-Decompositions: 3/8 LU Decompositions: 3/8 � Would it make sense when one system A ⋅ x = b Would it make sense when one system A x b becomes two L ⋅ y = b and U ⋅ x = y? � It depends; however both systems are very easy � It depends; however, both systems are very easy to solve if A = L ⋅ U is available. � Can use backward substitution to solve U ⋅ x = y . � C b k d b tit ti t l U ⎡ � � ⎤ ⎡ ⎤ ⎡ ⎤ u u u u x y 1 1 , 1 2 , 1 , 1 , i n 1 1 ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ � � 0 u u u x y ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 2 2 , 2 , 2 , i n 2 2 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ � � � � � � � ⋅ ⎢ = ⎢ ⎢ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ � � 0 0 u u x y ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ , , i i i n i i ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ � � � � � � � ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ � � ⎢ ⎥ ⎣ 0 0 0 u ⎦ ⎣ ⎦ ⎣ ⎦ x y , 37 n n n n
LU-Decompositions: 4/8 LU Decompositions: 4/8 � The L ⋅ y = b system is also easy to solve. y y y � From row 1, we have y 1 = b 1 . Row 2 is l 2,1 y 1 + y 2 = b 2 . Row i is 2 + + + + = � l y l y l y y b − − ,1 1 ,2 2 , 1 1 i i i i i i i − 1 ( ) i − ∑ � Hence, = − + + + = � y b l y l y l y b l y − − ,1 1 ,2 2 , 1 1 , i i i i i i i i i k k = 1 k ⎡ � � ⎤ ⎡ ⎤ ⎡ ⎤ 1 0 0 0 y b 1 1 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ � � 1 0 0 l y b ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 2 1 , 2 2 ⎢ ⎢ ⎥ ⎥ ⋅ ⎢ ⎢ ⎥ ⎥ = ⎢ ⎢ ⎥ ⎥ � � � � � � � ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ � � 1 0 l l y b row i ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ , 1 , 2 i i i i ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ � � � � � � � � � � � � � � ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ � � 38 ⎢ ⎥ ⎣ 1 l l l ⎦ ⎣ ⎦ ⎣ ⎦ y b , 1 , 2 , n n n i n n
LU-Decompositions: 5/8 LU Decompositions: 5/8 � The following code is based on the formula The following code is based on the formula mentioned earlier − 1 ( ( ) ) i − ∑ ∑ = − + + + + + + = � y y b b l y l y l y l y l l y y b b l l y y − − ,1 1 1 1 ,2 2 2 2 , 1 1 1 1 , i i i i i i i i i i i i i i i i i k i k k k = 1 k � This is referred to as forward substitution since y 1 , y 2 , …, y i -1 are used to compute y i . d t t y 1 = b 1 y 1 b 1 DO i = 2, n y i = b i DO k = 1, i-1 This is an O ( n 2 ) method y i = y i – L i,k *y k Do it yourself. END DO END DO END DO 39
LU-Decompositions: 6/8 LU Decompositions: 6/8 � In summary, LU-decomposition methods have In summary, LU decomposition methods have the following procedure: � From A in A ⋅ x = b find a LU-decomposition � From A in A ⋅ x = b find a LU-decomposition A = L ⋅ U � Solve for y with forward substitution from � S l f ith f d b tit ti f L ⋅ y = b � Solve for x with backward substitution from U ⋅ x = y 40
LU-Decompositions: 7/8 LU Decompositions: 7/8 � Why Why LU Why LU-decomposition Why LU decomposition decomposition w hen decomposition w hen w hen hen Gaussian elimination is available? Gaussian elimination is available? � The reason is simple: saving time � The reason is simple: saving time. � Suppose we need to solve k systems of linear equations like this: A ⋅ x 1 = b 1 , A ⋅ x 2 = b 2 , A ⋅ x 3 = ti lik thi A b A b A b 3 , …, A ⋅ x k =b k . Note that they share the same A A and not all b i ’s are available at the same time. d t ll b ’ il bl t th ti � Without a LU-decomposition, Gaussian elimination would be repeated k times, one for each system. This is time consuming. 41
LU Decompositions: 8/8 LU-Decompositions: 8/8 � Since each Gaussian elimination requires O ( n 3 ) Since each Gaussian elimination requires O ( n ) multiplications, solving k systems requires O ( k × n 3 ) multiplications. O ( k × n ) multiplications. � A LU-decomposition decomposes A = L ⋅ U. � For each b i , applying a forward followed by a � F h b l i f d f ll d b backward substitution yields the solution x i . � Since each backward and forward substitution requires O ( n 2 ) multiplications, solving k systems requires O ( n 3 + k × n 2 ) multiplications. 3 2 � Therefore, LU-decomposition is faster! elimination is still needed 42
How to Decompose: 1/4 How to Decompose: 1/4 � LU-decomposition is easier than you thought! LU decomposition is easier than you thought! � Gaussian elimination generates an upper triangular matrix which is the U in A = L ⋅ U triangular matrix, which is the U in A = L ⋅ U. � More importantly, the elimination process also produces lower triangular matrix L, although we t i L lth h d l t i l never thought about this possibility. � The next few slides shows how to recover this lower triangular matrix L during the elimination process. 43
How to Decompose: 2/4 How to Decompose: 2/4 ⎡ ⎤ 1 ⎢ ⎥ � When handling g a ⎢ ⎢ ⎥ ⎥ 2,1 , 1 1 row i , - a k , i / a i , i is ⎢ ⎥ a 1,1 ⎢ ⎥ multiplied to row column i ⎢ ⎥ a a 3,1 3,2 1 ⎢ ⎢ ⎥ ⎥ i and the result is i and the result is a a ⎢ ⎥ 1,1 2,2 added to row k . ⎢ ⎥ � � � ⎢ ⎥ � The entry on row The entry on row a a ⎢ ⎢ ⎥ ⎥ ,1 1 ,2 2 i i i i 1 1 ⎢ ⎥ k and column i of a a ⎢ 1,1 2,2 ⎥ = L L, where k > i , is ⎢ ⎥ a a a + + + 1,1 1,2 1, i i i i 1 ⎢ ⎢ ⎥ ⎥ a k , j / a j , j . / a a a ⎢ ⎥ 1,1 2,2 , i i ⎢ ⎥ � Thus, L can be � � � � ⎢ ⎥ generated on-the- generated on-the- ⎢ ⎢ ⎥ ⎥ a a a a a a ,1 1 ,2 2 , k k k k k i k i 1 ⎢ ⎥ fly as a k , i / a i , i is a a a a row k ⎢ ⎥ 1,1 2,2 , i i ⎢ ⎥ multiplier! � � � � � ⎢ ⎢ ⎥ ⎥ a a a a ⎢ ⎥ , ,1 ,2 , n j � n n n i 1 ⎢ ⎥ 44 a a a a ⎣ ⎦ 1,1 2,2 , , i i j j
How to Decompose: 3/4 How to Decompose: 3/4 � During Gaussian elimination, the lower During Gaussian elimination, the lower triangular part is set to zero. � One may use this portion for the lower � One may use this portion for the lower triangular matrix L without its diagonal; however the diagonal is part of U rather than L however, the diagonal is part of U rather than L. DO i = 1, n-1 ! using row i, i.e., a(i,i) DO k = i+1, n ! we want to eliminate a(k,i) S = a(k,i)/a(i,i) DO j = i+1, n ! for each entry on row k a(k,j) = a(k,j) - S*a(i,j) ! update its value (k j) (k j) * (i j) d i l END DO a(k,i) = S ! save this “multiplier” END DO ! don’t forget to pdate b(k) END DO ! don’t forget to update b(k) END DO 45
How to Decompose: 4/4 How to Decompose: 4/4 � After the decomposition A = L ⋅ U, matrix A is After the decomposition A L U, matrix A is replaced by U in the upper triangular part and L in the lower triangular part without the diagonal . in the lower triangular part without the diagonal . matrix U ⎡ ⎤ � � u u u u u 1,1 1,2 1,3 1, 1, i n ⎢ ⎢ ⎥ ⎥ � � l u u u u ⎢ ⎥ 2,1 2,2 2,3 2, 2, i n ⎢ ⎥ � � l l u u u 3 1 3,1 3 2 3,2 3,3 3 3 3, 3 3, 3 i i n n ⎢ ⎢ ⎥ ⎥ = ⎢ � � � � � � A ⎥ ⎢ ⎢ ⎥ ⎥ � � l l l l l l u u u u ⎢ ⎥ ,1 ,2 ,3 , , i i i i i i n � � � ⎢ ⎥ ⎢ ⎢ ⎥ ⎥ � � l l l l l l l l u ⎣ ⎦ ,1 ,2 ,3 , , n n n n i n n 46 matrix L
Example: 1/4 Example: 1/4 � Blue: values of the lower diagonal portion ( i.e ., g p ( , matrix L without the diagonal) × (-3/4)+ × (-0/4)+ × (-3/4)+ 4 0 1 1 3 1 3 1 0 1 2 0 3 2 4 1 4 0 1 1 ¾ 1 9/4 ¼ 0 1 2 0 column 1 ¾ ¾ 2 13/4 1/4 2 13/4 1/4 47
Example: 2/4 Example: 2/4 � Blue: values of the lower diagonal portion ( i.e ., g p ( , matrix L without the diagonal) 4 0 1 1 × (-1/1)+ × (-2/1)+ ¾ 1 9/4 ¼ 0 1 2 0 ¾ 2 13/4 1/4 4 4 0 1 1 0 1 1 ¾ 1 9/4 ¼ column 2 0 0 1 1 -¼ ¼ ¼ -¼ ¾ 2 -5/4 -¼ 48
Example: 3/4 Example: 3/4 � Blue: values of the lower diagonal portion ( i.e ., g p ( , matrix L without the diagonal) 4 0 1 1 ¾ 1 9/4 ¼ × (-5)+ 0 1 -¼ -¼ ¾ 2 -5/4 -¼ 4 4 0 1 1 0 1 1 ¾ 1 9/4 ¼ column 3 0 0 1 1 -¼ ¼ ¼ -¼ ¾ 2 5 1 49
Example: 4/4 Example: 4/4 � Verification: � Verification: ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ 1 4 0 1 1 4 0 1 1 ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ 3 4 / 1 1 9 4 / 1 4 / 3 1 3 1 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ ⋅ ⎥ = ⎢ ⎢ ⎥ ⎢ ⎢ − − − − ⎥ ⎢ ⎢ ⎥ ⎥ 0 0 1 1 1 1 1 4 1 4 / / 1 4 1 4 / / 0 0 1 1 2 2 0 0 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ 3 4 / 2 5 1 ⎦ ⎣ 1 ⎦ ⎣ 3 2 4 1 ⎦ L • U = A 50
How to Solve: 1/2 How to Solve: 1/2 � The following is a forward substitution solving The following is a forward substitution solving for y from L and b ( i.e ., L ⋅ y = b): ⎡ ⎤ � � u u u u u y 1 = b 1 1,1 1,2 1,3 1, 1, i n ⎢ ⎥ � � l u u u u DO i 2, n DO i = 2, n ⎢ ⎢ ⎥ ⎥ 2 1 2,1 2 2 2,2 2 3 2,3 2, 2 2 2, i i n n ⎢ ⎥ � � y i = b i l l u u u 3,1 3,2 3,3 3, 3, ⎢ i n ⎥ DO k = 1, i-1 = ⎢ � � � � � � A ⎢ ⎥ ⎥ ⎢ ⎥ y i = y i – a i,k *y k � l l l u u ⎢ ⎥ ,1 ,2 ,3 , , i i i i i i n END DO � � � ⎢ ⎥ END DO END DO ⎢ ⎢ ⎥ ⎥ � � ⎣ l l l l u ⎦ ,1 ,2 ,3 , , n n n n i n n a i , k is actually l i , k 51
How to Solve: 2/2 How to Solve: 2/2 � The following is a backward substitution solving The following is a backward substitution solving for x from U ⋅ x = y. ⎡ ⎤ � � u u u u u DO i = n, 1, -1 1,1 1,2 1,3 1, 1, i n ⎢ ⎥ � � S = y i S = y l l u u u u u u u u ⎢ ⎢ ⎥ ⎥ 2 1 2,1 2,2 2 2 2,3 2 3 2, 2 2, 2 i n ⎢ ⎥ � � DO k = i+1,n l l u u u 3,1 3,2 3,3 3, 3, i n ⎢ ⎥ = ⎢ S = S - a i k *x k � � � � � � A ⎢ ⎥ ⎥ i,k k ⎢ ⎥ END DO � l l l u u ⎢ ⎥ i ,1 i ,2 i ,3 i i , i n , x i = S/a i,i � � � ⎢ ⎢ ⎥ ⎥ END DO ⎢ ⎢ ⎥ ⎥ � � l l l l u ⎣ ⎦ ,1 ,2 ,3 , , n n n n i n n a i,k here is u i,k ! 52
Why Is L Correct? 1/10 Why Is L Correct? 1/10 � For row 1, the elimination step multiplies every For row 1, the elimination step multiplies every element on row 1 by – a i ,1 / a 1,1 and adds the results to row i . This is a matrix multiplication! to row i . This is a matrix multiplication! ⎛ ⎞ a − × + ⎜ ,1 ⎟ i ⎡ � � ⎤ a a 1 0 0 0 0 ⎜ ⎟ 1,1 ,1 i ⎢ ⎥ a ⎝ ⎠ 1,1 a a ⎢ ⎢ ⎥ ⎥ − 2,1 2 1 � � 1 0 0 0 ⎢ ⎥ ⎡ ⎤ ⎡ ⎤ a � � � � a a a a a a a a a a 1,1 ⎢ ⎥ 1,1 1,2 1,3 1, 1, 1,1 1,2 1,3 1, 1, i n i n ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ � � � � a 0 a a a a a a a a a ⎢ ⎥ ⎢ ⎥ − 3,1 � � 2,1 2,2 2,3 2, 2, 2,2 2,3 2, 2, i n i n 0 1 0 0 ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ � � � � a 0 0 a a a a a a a a a ⎢ ⎥ 1,1 3,1 3,2 3,3 3, 3, 3,2 3,3 3, 3, ⎢ i n ⎥ ⎢ i n ⎥ ⎢ ⎥ = i � � � � � � � � � � � � � � � ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ � � 0 a a a a a a a a a a ⎢ ⎥ − ⎢ ⎢ ,1 ,2 ,3 , , ⎥ ⎥ ⎢ ⎢ ,2 ,3 , , ⎥ ⎥ ,1 i i i i i i n i i i i i n i � � 0 0 0 0 1 1 0 0 ⎢ ⎢ ⎥ ⎥ � � � � � � � � � � � � � � ⎢ ⎥ ⎢ ⎥ a ⎢ ⎥ 1,1 ⎢ ⎥ ⎢ ⎥ � � � � 0 ⎢ � � � ⎥ a a a a a a a a a ⎣ ⎦ ⎣ ⎦ ,1 ,2 ,3 , , ,2 ,3 , , n n n n i n n n n n i n n ⎢ ⎥ a ⎢ ⎢ ⎥ ⎥ − new values l n n ,1 1 � � 0 0 0 0 0 0 1 1 ⎢ ⎥ ⎣ a ⎦ row i 1,1 53
Why Is L Correct? 2/10 Why Is L Correct? 2/10 � Define matrix E 1 as follows. Then, E 1 ⋅ A sets all Define matrix E 1 as follows. Then, E 1 A sets all entries on column 1 below a 1,1 to zero. ⎡ ⎤ � � 1 0 0 0 0 ⎢ ⎥ − � � / 1 0 0 0 a a ⎢ ⎢ ⎥ ⎥ 2 1 2,1 1,1 11 ⎢ ⎥ − � � / 0 1 0 0 a a 3,1 1,1 ⎢ ⎥ = ⎢ � � � � � � E E ⎢ ⎥ ⎥ E i l E 1 is lower triangular! t i l ! 1 ⎢ ⎥ − � � / 0 0 1 0 a a ⎢ ⎢ ⎥ ⎥ ,1 1,1 i � � � ⎢ ⎥ ⎢ ⎥ − � � 1 / 0 0 0 1 a a ⎣ ⎣ ⎦ ⎦ ,1 , n n n n n n 54
Why Is L Correct? 3/10 Why Is L Correct? 3/10 � Define E 2 as follows. The same calculation shows Define E 2 as follows. The same calculation shows E 2 ⋅ (E 1 ⋅ A) sets all entries of E 1 ⋅ A below a 2,2 to 0. ⎡ ⎤ � � 1 0 0 0 0 ⎢ ⎥ � � � � 0 0 1 1 0 0 0 0 0 0 ⎢ ⎢ ⎥ ⎥ ⎢ ⎥ − � � 0 / 1 0 0 a a 3,2 2,2 ⎢ ⎢ ⎥ ⎥ = ⎢ � � � � � � � � E 2 is lower triangular! E i l t i l ! E ⎥ 2 ⎢ ⎥ − � � 0 / 0 1 0 a a ⎢ ⎢ ⎥ ⎥ ,2 2,2 i � � � � ⎢ ⎥ ⎢ ⎥ − � � 0 0 / / 0 0 0 0 1 1 a a a a ⎣ ⎣ ⎦ ⎦ 55 ,2 2 2,2 2 2 n
Why Is L Correct? 4/10 Why Is L Correct? 4/10 � For column i , lower triangular matrix E i is For column i , lower triangular matrix E i is defined as follows. It has all the “multipliers”. ⎡ ⎤ � � 1 0 0 0 0 ⎢ ⎥ � � � 0 0 0 0 ⎢ ⎢ ⎥ ⎥ ⎢ ⎥ � � 0 1 0 0 ⎢ ⎥ a + ⎢ ⎥ − � � i 1, i � � � � ⎢ ⎢ ⎥ ⎥ a ⎢ ⎥ , i i = E ⎢ ⎥ i a + − i 2, i � � 0 1 0 ⎢ ⎢ ⎥ ⎥ a ⎢ ⎥ , i i ⎢ ⎥ � � � � ⎢ ⎥ ⎢ ⎢ ⎥ ⎥ a − , � � n i 0 0 1 ⎢ ⎥ a ⎣ ⎦ 56 , i i
Why Is L Correct? 5/10 Why Is L Correct? 5/10 � Since E i 1 ⋅ E i 2 ⋅ … ⋅ E 2 ⋅ E 1 ⋅ A sets the lower diagonal Since E i -1 E i -2 … E 2 E 1 A sets the lower diagonal part of column 1 to column i -1 to 0, the new matrix E i ⋅ (E i 1 ⋅ E i 2 ⋅ … ⋅ E 2 ⋅ E 1 ⋅ A ) eliminates the matrix E i (E i -1 E i -2 … E 2 E 1 A ) eliminates the lower diagonal part of a i , i . ⎡ � � ⎤ 1 0 0 0 0 ⎢ ⎥ � � � 0 0 0 0 ⎢ ⎥ ⎡ ⎤ ⎡ ⎤ � � � � a a a a a a a a a a ⎢ ⎥ � � + + + + 1,1 1, 1, 1 1, 2 1, 1,1 1, 1, 1 1, 2 1, 0 1 0 0 i i i n i i i n ⎢ ⎥ ⎢ ⎥ ⎢ ⎢ ⎥ ⎥ � � � � � � � � � � � � � � � � � � � � � � 0 0 a a ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ 2, n ⎢ a ⎥ + − � i 1, i � � ⎢ ⎥ ⎢ ⎥ � � � � 0 0 a a a a a a a a ⎢ ⎥ + + + + , , 1 , 2 , , , 1 , 2 , ⎢ i i i i i i i n ⎥ = ⎢ i i i i i i i n ⎥ a ⎢ ⎥ i i , � � i 0 a a a a a a a ⎢ ⎥ ⎢ ⎥ + + + + + + + + + + + ⎢ ⎥ 1, 1, 1 1, 2 1, 1, 1 1, 2 1, i i i i i i i n i i i i i n a ⎢ ⎥ ⎢ ⎥ + − i 2, i � � 0 1 0 0 0 0 ⎢ ⎢ ⎥ ⎥ a a a a a a a + + + + + + + + + + + + + + + + + + + + + + ⎢ ⎢ 2, 2 2 2, 1 1 2 2, 2 2 2 2, ⎥ ⎥ ⎢ ⎢ 2, 2 1 1 2 2, 2 2 2, 2 ⎥ ⎥ i i i i i i i i i i i i i i n n i i i i i i i i i i n n a ⎢ ⎥ i i , � � � � � � � � � � � � ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ � � � � ⎢ ⎥ ⎢ ⎥ � � � � ⎢ ⎥ 0 0 0 a a a a a a a ⎣ ⎦ ⎣ ⎦ + + + + , , 1 , 2 , , 1 , 2 , n i n i n i n n n i n i n n ⎢ ⎥ a − E i -1 ⋅ E i -2 ⋅ … ⋅ E 2 ⋅ E 1 ⋅ A E ⋅ E ⋅ ⋅ E ⋅ E ⋅ A E ⋅ E ⋅ E ⋅ E i ⋅ E i -1 ⋅ E i -2 ⋅ … ⋅ E 2 ⋅ E 1 ⋅ A ⋅ E ⋅ E ⋅ A n i , � � 0 0 1 ⎢ ⎢ ⎥ ⎥ ⎣ ⎣ a a ⎦ ⎦ i i , E i 57
Why Is L Correct? 6/10 Why Is L Correct? 6/10 � Repeating this process, matrices E 1 , E 2 , …, E n -1 p g p 1 2 n 1 are constructed so that E n -1 ⋅ E n -2 ⋅ … ⋅ E 2 ⋅ E 1 ⋅ A is upper diagonal. � N t th t � Note that only E 1 , E 2 , …, E n -1 are needed, l E E E d d because the lower diagonal part has n -1 columns. � Therefore the elimination process is completely � Therefore, the elimination process is completely described by the product of matrices E i ’s. � Now, we have E n 1 ⋅ E n 2 ⋅ … ⋅ E 2 ⋅ E 1 ⋅ A = U, where U Now, we have E n -1 E n -2 … E 2 E 1 A U, where U is an upper triangular matrix, and A = (E n -1 ⋅ E n - 2 ⋅ … ⋅ E 2 ⋅ E 1 ) -1 ⋅ U, where T -1 means the inverse matrix of T ( i.e ., T ⋅ T -1 = I, the identify matrix). T T -1 t i f T ( i I th id tif t i ) � What is (E n -1 ⋅ E n -2 ⋅ … ⋅ E 2 ⋅ E 1 ) -1 ? In fact, this is matrix L we want! matrix L we want! 58
Why Is L Correct? 7/10 Why Is L Correct? 7/10 � In linear algebra, the inverse of A ⋅ B, (A ⋅ B) -1 , is In linear algebra, the inverse of A B, (A B) , is computed as (A ⋅ B) -1 = B -1 ⋅ A -1 . ⋅ E ⋅ E ) -1 = E -1 ⋅ E -1 ⋅ � (E n -1 ⋅ E n -2 ⋅ … ⋅ E 2 ⋅ E 1 ) = E 1 ⋅ E 2 ⋅ … ⋅ E n -2 ⋅ E n -1 . ⋅ E ⋅ ⋅ E -1 ⋅ E � (E -1 -1 is very easy to compute! See below. � E i column i column i ⎡ ⎤ ⎡ ⎤ 1 1 ⎢ ⎥ ⎢ ⎥ � � � � ⎢ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ 1 1 + ⎢ ⎥ ⎢ ⎥ − 1 1 m m ⎢ ⎥ ⎢ ⎥ + + ⋅ ⋅ = = 1 1 i i I I ⎢ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ � � � � n row k ⎢ ⎥ ⎢ ⎥ − ⎢ 1 ⎥ ⎢ 1 ⎥ m m k k ⎢ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ � � � � � � � � ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ ⎢ ⎥ ⎣ 1 ⎦ ⎣ 1 ⎦ m m 59 n n
Why Is L Correct? 8/10 Why Is L Correct? 8/10 -1 is obtained by removing the � Therefore, E i Therefore, E i is obtained by removing the negative signs from E i ! − 1 ⎡ ⎤ ⎡ ⎤ 1 1 ⎢ ⎥ ⎢ ⎥ � � ⎢ ⎥ ⎢ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ 1 1 1 1 ⎢ ⎥ ⎢ ⎥ a a ⎢ ⎥ ⎢ ⎥ + + − 1, 1, i i i i 1 1 ⎢ ⎥ ⎢ ⎥ a a ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ i i i i , i i i i , = ⎢ � � ⎥ ⎢ � � ⎥ ⎢ ⎥ ⎢ ⎥ a a ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ − , , k i k i 1 1 1 1 ⎢ ⎥ ⎢ ⎥ a a i i , i i , ⎢ ⎥ ⎢ ⎥ � � � � ⎢ ⎥ ⎢ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ a a ⎢ ⎥ ⎢ ⎥ − , , n i n i 1 1 60 ⎢ ⎥ ⎢ ⎥ a a ⎣ ⎦ ⎣ ⎦ , , i i i i
Why Is L Correct? 9/10 Why Is L Correct? 9/10 -1 is also easy -1 ⋅ E 2 -1 ⋅ … ⋅ E n 2 -1 ⋅ E n 1 � Additionally, E 1 Additionally, E 1 E 2 … E n -2 E n -1 is also easy to compute. -1 It � The following shows the results of E n -2 ⋅ E n -1 . It -1 ⋅ E � The following shows the results of E is basically “filling” the entries. ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ 1 1 1 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ � � � � � � ⎢ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = i 1 1 0 1 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 1 1 0 0 1 1 1 1 p p ⎢ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ − − 1 1 n n ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ 0 1 ⎦ ⎣ 0 1 ⎦ ⎣ 1 ⎦ p m p m n n n n 61
Why Is L Correct? 10/10 Why Is L Correct? 10/10 ⎡ ⎤ 1 -1 ⋅ E 2 -1 ⋅ … ⋅ E n -2 -1 ⋅ E n -1 � E 1 ⎢ ⎥ -1 a a 1 2 n -2 n -1 ⎢ ⎢ ⎥ ⎥ 2 1 2,1 1 is computed by ⎢ ⎥ a 1,1 ⎢ ⎥ ⎢ ⎥ a a removing the signs g g 3,1 3,2 1 ⎢ ⎢ ⎥ ⎥ a a ⎢ ⎥ 1,1 2,2 from all E i ’s and ⎢ ⎥ � � � ⎢ ⎥ collect them all a a ⎢ ⎢ ⎥ ⎥ ,1 ,2 i i 1 1 ⎢ ⎢ ⎥ ⎥ a a together into a ⎢ ⎥ − − − = 1,1 2,2 1 i 1 i�i 1 E E E − 1 2 1 n ⎢ ⎥ a a a + + + 1,1 1,2 1, lower triangular g i i i i 1 ⎢ ⎥ a a a a a a ⎢ ⎢ ⎥ ⎥ 11 1,1 2,2 2 2 i i i i , matrix . ⎢ ⎥ � � � � ⎢ ⎥ ⎢ a a a ⎥ � Verify this yourself! y y ,1 ,2 , k k k i 1 ⎢ ⎢ ⎥ ⎥ a a a ⎢ ⎥ 1,1 2,2 , i i ⎢ ⎥ Note that the a i , j ’s are not the � � � � � ⎢ ⎥ original of matrix A. They are a a a a ⎢ ⎢ ⎥ ⎥ ,1 ,2 , , n n n i n j � � 1 1 ⎢ ⎥ computed intermediate results a a a a ⎣ ⎦ 1,1 2,2 , , i i j j in the elimination process 62
Pivoting in LU Decomposition: 1/7 Pivoting in LU-Decomposition: 1/7 � LU-Decomposition may also use pivoting. LU Decomposition may also use pivoting. � Although partial pivoting in Gaussian elimination does not affect the order of the elimination does not affect the order of the solutions, it does in LU-decomposition because LU-decomposition only processes matrix A! LU-decomposition only processes matrix A! � Therefore, when using partial pivoting with LU- decomposition, an index array is needed to save d iti i d i d d t row swapping activities so that the same row swapping can also apply to matrix b later. i l l t t i b l t 63
Pivoting in LU-Decomposition: 2/7 Pivoting in LU Decomposition: 2/7 � To keep track partial pivoting , an index array To keep track partial pivoting , an index array idx() is needed, which will be used to swap elements of matrix b. elements of matrix b. � For example, if idx(1)=4 , idx(2)=1 , idx(3)=3 and idx(4)=2 this means after row idx(3)=3 and idx(4)=2 , this means after row swapping, row 1 is equation 4, row 2 is equation 1, row 3 is equation 3, and row 4 is equation 2. 1 row 3 is equation 3 and row 4 is equation 2 � Before using forward/backward substitution, we need to swap b 4 to the first position, b 1 , b 3 and b 2 d t b t th fi t iti d b b b to the 2 nd , 3 rd and 4 th , respectively. 64
Pivoting in LU-Decomposition: 3/7 Pivoting in LU Decomposition: 3/7 � Blue: values of the lower diagonal portion ( i.e ., g p ( , matrix L without the diagonal) index array array × (-3/4)+ × (-0/4)+ × (-3/4)+ 1 4 0 1 1 2 3 1 3 1 3 3 0 2 2 0 0 2 2 0 4 3 3 4 1 4 0 1 1 ¾ 1 9/4 ¼ 0 2 2 0 column 1 ¾ ¾ 3 13/4 ¼ 3 13/4 ¼ 65
Pivoting in LU-Decomposition: 4/7 Pivoting in LU Decomposition: 4/7 � Blue: values of the lower diagonal portion ( i.e ., g p ( , matrix L without the diagonal) index index array array array 1 4 0 1 1 1 4 0 1 1 2 ¾ 1 9/4 ¼ × (-2/3)+ × (-1/3)+ 4 ¾ 3 13/4 ¼ 3 0 2 2 0 3 0 2 2 0 4 ¾ 3 13/4 ¼ 2 ¾ 1 9/4 ¼ index array column 2 1 1 4 4 0 1 1 0 1 1 4 ¾ 3 13/4 ¼ 3 3 0 0 2/3 -1/6 2/3 1/6 1/6 -1/6 2 ¾ 1/3 7/6 1/6 66
Pivoting in LU-Decomposition: 5/7 Pivoting in LU Decomposition: 5/7 � Blue: values of the lower diagonal portion ( i.e ., g p ( , matrix L without the diagonal) index index array array 1 4 0 1 1 1 4 0 1 1 4 ¾ 3 13/4 ¼ 4 ¾ 3 13/4 ¼ × (1/7)+ 3 0 2/3 -1/6 -1/6 2 ¾ 1/3 7/6 1/6 2 ¾ 1/3 7/6 1/6 3 0 2/3 -1/6 -1/6 index array 1 1 4 4 0 0 1 1 1 1 4 ¾ 3 13/4 ¼ column 3 2 2 ¾ 1/3 7/6 1/6 ¾ 1/3 7/6 1/6 3 0 2/3 -1/7 -3/21 67
Pivoting in LU-Decomposition: 6/7 Pivoting in LU Decomposition: 6/7 � Verification: � Verification: ⎡ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎤ 1 4 0 1 1 index ⎢ ⎥ ⎢ ⎥ array 3 13 1 ⎡ ⎡ ⎤ ⎤ ⎢ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ 4 0 1 1 1 1 1 1 3 3 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 4 4 4 3 3 4 1 4 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⋅ = ⎢ 3 1 7 1 ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ 1 1 3 3 1 1 3 3 1 1 2 2 ⎢ ⎥ ⎢ 4 3 ⎥ ⎢ 6 6 ⎥ 3 ⎣ ⎦ 0 2 2 0 ⎢ ⎥ ⎢ ⎥ − − 2 1 3 ⎢ ⎢ 0 0 1 1 ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ ⎣ ⎦ ⎣ ⎦ 3 7 21 L • L U U = A A 68
Pivoting in LU-Decomposition: 7/7 Pivoting in LU Decomposition: 7/7 � LU-decomposition may use complete pivoting. LU decomposition may use complete pivoting. � A second index array is needed to keep track the swapping of variables swapping of variables. � After the decomposition A = L ⋅ U is found, the partial pivoting index array is used to swap ti l i ti i d i d t matrix b’s elements, then use the new b to solve f for the x. th � If complete pivoting is used, the second array is used to get the correct order of variables back. � See Complete Pivoting Complete Pivoting slides for the details. 69
Iterative Methods: 1/2 Iterative Methods: 1/2 � In addition to Gaussian elimination, LU- In addition to Gaussian elimination, LU decomposition, and other methods that guarantee to compute the solution in a fixed guarantee to compute the solution in a fixed number of steps, there are iterative methods . � Like we learned in solving non-linear equations � Like we learned in solving non-linear equations, iterative methods start with an initial guess x 0 of A • x = b and successively compute x 1 , x 2 , …, x k A • x = b and successively compute x x x until the error vector e k = A • x k -b is small enough. � In many cases ( e.g ., very large systems), iterative � I ( l t ) it ti methods may be more effective than the non- it iterative methods. ti th d 70
Iterative Methods: 2/2 Iterative Methods: 2/2 � There are many iterative methods, from the very There are many iterative methods, from the very classical ones ( e.g ., Jacobi and Gauss-Seidel) to very modern ones ( e.g ., conjugate gradient). very modern ones ( e.g ., conjugate gradient). � We only discuss the classical ones, Jacobi and Gauss-Seidel Gauss-Seidel. � Typical iterative methods have a general form similar to the following: i il t th f ll i + = = ⋅ ⋅ + + x x C x C x d d 1 k k 71
Jacobi Method: 1/4 Jacobi Method: 1/4 � Equation i in A • x=b has the following form: Equation i in A x b has the following form: + + + + + = � � a x a x a x a x b , 1 1 , 2 2 , , i i i i i i n n i � Solving for x i yields the following if a i , i ≠ 0: ≠ 0 � S l i f i ld th f ll i if ⎡ ⎤ [ [ ] ] ( ( ) ) n 1 1 ∑ ∑ = − + + + + + + + + + + = − � � ⎢ ⎢ ⎥ ⎥ x x b b a x a x a a x x a a x x a a x x b b a a x x − − + + , 1 1 1 1 , 1 1 1 1 , 1 1 1 1 , , i i i i i i i i i i i i i i i i i i i n i n i i i k i k k k ⎣ ⎦ a a = ≠ 1 , , , k k i i i i i � Jacobi method goes as follows: Jacobi method goes as follows: � Start with an initial guess x 0 = [ x 1 , x 2 , …, x n ] � Plug this x i into the above equations to compute x i +1 � Plug this x i into the above equations to compute x i +1 � Repeat this process until A • x k ≈ b for some k . � See next few slides for more details � See next few slides for more details. 72
Jacobi Method: 2/4 Jacobi Method: 2/4 � Let the system be: Let the system be: -5 x – y + 2 z = 1 x = -(1 + y - 2 z )/5 transformation 2 x + 6 y - 3 z = 2 y = (2 - 2 x + 3 z )/6 2 x + y + 7 z = 32 z = (32 - 2 x - y )/7 � Suppose the initial values are x = y = z =0 ( i.e ., x 0 = [0,0,0]) � Plug x 0 into the equations: x 1 =[-1/5,1/3,32/7]. 0 1 � Plug x 1 into the equations: ⎡ ⎤ ⎛ ⎛ ⎞ ⎞ 1 1 1 32 − + + − × × ⎜ ⎜ 5 1 2 2 ⎢ ⎢ ⎥ ⎥ ⎝ ⎠ 3 7 − ⎡ ⎤ ⎢ ⎥ 15619 . × − ⎛ ⎞ Converge in approx. 10 iterations ⎢ ⎥ 1 1 32 ⎢ ⎥ = − + × = ⎜ ⎟ 6 2 2 3 2 6857 . x 2 ⎢ ⎥ ⎢ ⎥ with x ≈ 0.998, y ≈ 1.997, z ≈ 3.991 ⎝ ⎠ , y , 5 7 ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ ⎣ ⎣ 4 5810 4 5810 . ⎦ ⎦ × − ⎛ ⎞ 1 1 1 ⎢ − − ⎥ ⎜ ⎟ 7 32 2 ⎝ ⎠ ⎢ ⎥ ⎣ ⎦ 73 5 3
Jacobi Method: 3/4 Jacobi Method: 3/4 � Ste Step 1: Initialization p 1: Initialization ⎡ ⎡ ⎤ ⎤ a a a n − − − 1 2 1 2 , 1 3 1 3 , 1 1 , � n 0 0 ⎢ ⎢ ⎥ ⎥ a a a ⎢ ⎥ 1 1 , 1 1 , 1 1 , ⎢ ⎥ a a a − − − 2 1 , 2 3 , 2 , � n ! Given system is A • x = b 0 ⎢ ⎥ a a a ⎢ ⎢ ⎥ ⎥ 2 2 2 2 , 2 2 2 2 , 2 2 2 2 , ! C(:,:) and d(:) are two = − C ⎢ a a a ⎥ − − 3 1 , 3 2 , 3 , � n 0 ⎢ ⎥ ! working arrays a a a ⎢ ⎥ 3 3 , 3 3 , 3 3 , � � � � ⎢ ⎥ ⎢ ⎢ ⎥ ⎥ a a a DO i = 1, n − − − , 1 , 2 , 3 n n n � 0 ⎢ ⎥ ⎣ a a a ⎦ DO j = 1, n , , , n n n n n n C(i,j) A(i,j)/A(i,i) C(i,j) = -A(i,j)/A(i,i) ⎡ ⎡ ⎤ ⎤ b b 1 ⎢ ⎥ END DO a ⎢ ⎥ 1 1 , C(i,i) = 0 ⎢ ⎥ b 2 ⎢ ⎢ ⎥ ⎥ d(i) = b(i)/A(i i) d(i) = b(i)/A(i,i) a a ⎢ ⎥ 2 2 , = END DO d b ⎢ ⎥ 3 ⎢ ⎥ Now we have x = C • x + d a 3 3 , ⎢ ⎢ ⎥ ⎥ � � ⎢ ⎢ ⎥ ⎥ b ⎢ ⎥ n row and column swaps may be needed 74 ⎢ ⎥ ⎣ a ⎦ , n n before starting iteration!
Jacobi Method: 4/4 Jacobi Method: 4/4 � Step Step 2: Step Step 2: 2: Iteration 2: Iteration Iteration Iteration DO x new = C*x + d x_new = C*x + d IF (ABS(A*x_new - b) < Tol) EXIT x = x_new END DO END DO ABS(A* ABS(A*x_new-b) < Tol means the absolute value b) < T l th b l t l of each entry of A*x_new – b is less than Tol . You may just use the equations for computation instead of the above matrix form. 75
Gauss-Seidel Method: 1/4 Gauss Seidel Method: 1/4 � Gauss-Seidel method is slightly different from g y Jacobi method. � With Jacobi method, all new x values are , computed and be used in the next iteration. � Gauss-Seidel method uses the same set of equations and the same initial values. � However , the new x 1 computed from equation 1 , p q 1 is used in equation 2 to compute x 2 , the new x 1 and x 2 are used in equation 3 to compute x 3 . � In general, the new x 1 , x 2 , …, x i -1 are used in equation i to compute the new x i . 76
Gauss-Seidel Method: 2/4 Gauss Seidel Method: 2/4 � Example: p -5 x – y + 2 z = 1 x = -(1 + y - 2 z )/5 transformation 2 x + 6 y - 3 z = 2 y = (2 - 2 x + 3 z )/6 2 2 x + y + 7 z = 32 7 32 z = (32 - 2 x - y )/7 � Initial value is x = y = z =0. � Plugging y and z into equation 1 yields the new � Plugging y and z into equation 1 yields the new x = -1/5. Now, we have x = -1/5, y = z =0. � Plugging x and z into equation 2 yields the new � Plugging x and z into equation 2 yields the new y = 2/5. Now, we have x =-1/5, y =2/5, z =0 � Plugging x and y into equation 3 yields the new � Plugging x and y into equation 3 yields the new z = 32/7. Now, we have x =-1/5, y =2/5, z =32/7. � This completes the first iteration! � This completes the first iteration! 77
Gauss-Seidel Method: 3/4 Gauss Seidel Method: 3/4 � Example: p -5 x – y + 2 z = 1 x = -(1 + y - 2 z )/5 transformation 2 x + 6 y - 3 z = 2 y = (2 - 2 x + 3 z )/6 2 2 x + y + 7 z = 32 7 32 z = (32 - 2 x - y )/7 � Iteration 2 starts with x =-1/5, y =2/5, z =32/7. � Plugging y and z into equation 1 yields the new x = � Plugging y and z into equation 1 yields the new x = 271/175 ≈ 1.5486, and x = 271/175, y =2/5, z =32/7. � Plugging x and z into equation 2 yields the new y = � Plugging x and z into equation 2 yields the new y = 368/175 ≈ 2.1029, and x =271/175, y =368/175, z =32/7 � Plugging x and y into equation 3 yields the new z = � Plugging x and y into equation 3 yields the new z = 134/35 ≈ 3.8286, and x =271/175, y =368/175, z =134/35. � This completes the second iteration! � This completes the second iteration! 78
Gauss-Seidel Method: 4/4 Gauss Seidel Method: 4/4 � Algorithm: � Algorithm: Algorithm: Algorithm: � The initialization part is the same as the Jacobi method since both methods use the same set of method, since both methods use the same set of equation transformations. DO DO DO i = 1, n ! Update x i EQN_i = 0 DO j = 1, n EQN_i = EQN_i + C(i,j)*x(j) END DO x(i) = EQN_i + d(i) END DO IF (ABS(A*x – b) < Tol) EXIT ( ( ) ) END DO 79
Convergence Convergence � Jacobi and Gauss-Seidel methods may not Jacobi and Gauss Seidel methods may not converge. � Even though they may converge they may be � Even though they may converge, they may be very slow. � If the system is diagonal dominant , both methods � If th t i di t b th th d l d i converge. > ∑ n ∑ > | | | | | | | | a a a a (for every i ) (for every i ) , , i i i j = ≠ 1 , j j i 80
Geometric Meaning: 1/4 Geometric Meaning: 1/4 � Suppose we have two lines, #1 and #2. Suppose we have two lines, #1 and #2. � Jacobi method uses y to find a new x with equation #1 and uses x to find a new y with equation #1, and uses x to find a new y with equation #2. #1 #1 x 2 x x 1 x 0 #2 81
Geometric Meaning: 2/4 Geometric Meaning: 2/4 � If the given system is If the given system is 3 x + y =3 x =1- y /3 diagonal dominant , the x +2 y =2 y =1- x /2 Jacobi method converges. Jacobi method converges. #1 � Let X 0 =[3,2]. � Then, X 1 =[1/3,-1/2]. The � Th X [1/3 1/2] Th X 0 right diagram shows how to find the x and y coordinates. fi d th d di t � X 2 =[7/6,5/6] and X 2 X 3 =[13/18,5/12]. X 3 � The solution is X * =[4/5,3/5] #2 #2 X 1 82
Geometric Meaning: 3/4 Geometric Meaning: 3/4 � Let us try the Gauss-Seidel method. Let us try the Gauss Seidel method. � This method uses y to compute a new x with equation 1 which replaces the old x and uses equation 1, which replaces the old x , and uses this new x to compute a new y . #1 #1 x 1 x 2 x 0 #2 #2 83
Geometric Meaning: 4/4 Geometric Meaning: 4/4 � Use Gauss-Seidel method Use Gauss Seidel method 3 x + y =3 x =1- y /3 with X 0 =[3,2]. x +2 y =2 y =1- x /2 � Since x =1-2/3=1/3 and y =1- � Since x =1-2/3=1/3 and y =1- #1 (1/3)/2=5/6, X 1 =[1/3,5/6]. � Si � Since x =1-(5/6)/3=13/18 and 1 (5/6)/3 13/18 d X 0 y =1-(13/18)/2=23/36, X 2 =[13/18,23/36] . X [13/18 23/36] � Since x =1-(23/36)/3=85/108 X 1 X 1 X 2 and y =1-(85/108)/2=131/216, X 3 =[85/108,131/216]. #2 #2 � Gauss-Seidel is faster ! 84
Iterative Refinement: 1/5 Iterative Refinement: 1/5 � Due to rounding error accumulation, elimination Due to rounding error accumulation, elimination methods may not be accurate enough. � This is especially true if matrix A in A ⋅ x = b is ill- � This is especially true if matrix A in A ⋅ x = b is ill- conditioned ( e.g ., nearly singular), and, as a result the computed solution may still be far result, the computed solution may still be far away from the actual solution. � It � Iterative refinement techniques can improve the i th ti fi t t h i accuracy of elimination methods. � To use iterative refinement methods, one has to preserve matrix A and computes the LU- decomposition A = L ⋅ U. 85
Iterative Refinement: 2/5 Iterative Refinement: 2/5 � Here is the algorithm. In each step, the Here is the algorithm. In each step, the “residual” is computed and used to solve for a “correction” to update the solution X. correction to update the solution X. Make copies of A to A’ and A” and b to b’ U Use A” to compute A” = L*U ” t t ” L*U ! A” is destroyed ! ” i d t d Apply an elimination method to solve A’*x=b’ Let the solution be x DO r = b – A * x ! compute residual IF (ABS(r) < Tol) EXIT Solve for ∆ from (L*U)* ∆ = r ! compute correction x = x + ∆ ! update x END DO END DO 86
Iterative Refinement: 3/5 Iterative Refinement: 3/5 � Consider solving the following system: Consider solving the following system: x + y = 2 2 + 3 2 x + 3 y = 5 5 � We have A, L, U and b matrices as follows: A = ⎡ ⎡ ⎥ = ⎡ ⎤ ⎤ ⎡ ⎤ ⎥ ⋅ ⎡ ⎤ ⎡ ⎤ ⎤ ⎡ ⎤ ⎡ ⎤ 1 1 1 1 1 2 = ⎢ ⎥ ⎢ ⎢ ⎢ ⎥ b ⎣ 2 3 ⎦ ⎣ 2 1 ⎦ ⎣ 1 ⎦ ⎣ ⎦ 5 L L U U � If a method provides a “solution” x = 0.9 and y = 1 3 the residual vector r is 1.3, the residual vector r is ⎥ = − = ⎡ ⎥ − ⎡ ⎤ ⎤ ⎥ ⋅ ⎡ ⎤ ⎡ ⎤ 2 1 1 09 . 02 . = − ⋅ r b A x ⎢ ⎢ ⎦ ⎥ ⎢ ⎢ ⎥ ⎣ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎥ − ⎣ ⎣ 5 5 ⎦ ⎣ ⎣ 2 2 3 3 ⎦ ⎦ ⎣ 13 13 . ⎦ ⎦ ⎣ ⎣ 07 07 . ⎦ ⎦ 87
Iterative Refinement: 4/5 Iterative Refinement: 4/5 � Now, we solve for ∆ from A ⋅∆ =r. Now, we solve for ∆ from A ∆ r. � Since A=L ⋅ U, we have L ⋅ (U ⋅∆ )=r. � Since L and r are known, we may solve for T in � Si L d k l f T i L ⋅ T = r (hence U ⋅∆ = T) as follows: ⎥ = − ⎡ ⎡ ⎤ ⎥ ⋅ ⎡ ⎤ ⎡ ⎤ ⎤ ⎡ ⎡ ⎤ ⎤ 1 1 02 02 . t t 1 ⎢ ⎢ ⎢ ⎥ − ⎣ 2 1 ⎦ ⎣ ⎦ ⎣ 07 . ⎦ t 2 L L � Forward substitution yields t 1 = -0.2 and t 2 = -0.3 − ⎡ ⎡ ⎤ ⎤ 0.2 0 2 = ⎢ ⎥ T − ⎣ ⎣ 0.3 ⎦ ⎦ 88
Iterative Refinement: 5/5 Iterative Refinement: 5/5 � We have U ⋅∆ = T as follows: We have U ∆ T as follows: ∆ ⎥ = − ⎡ ⎥ ⋅ ⎡ ⎤ ⎤ ⎡ ⎤ 1 1 02 . 1 ⎢ ⎢ ⎢ ⎥ ∆ ∆ − ⎣ ⎣ 1 1 ⎦ ⎦ ⎣ ⎣ ⎦ ⎦ ⎣ ⎣ 03 03 . ⎦ ⎦ 2 U � Backward substitution yields ∆ as follows: ⎡ ⎡ ⎤ ⎤ 0.1 ∆ = ⎢ ⎥ − ⎣ ⎣ ⎦ ⎦ 0.3 � The new x is ⎡ ⎡ ⎤ ⎤ ⎡ ⎡ ⎤ ⎤ ⎡ ⎤ ⎡ ⎤ 0.9 0 9 0 1 0.1 1 1 = + ∆ = + = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ x x − new old ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ 1.3 0.3 1 � Since A ⋅ x = b, we have the solution in 1 iteration! 89
Matrix Inversion: 1/9 Matrix Inversion: 1/9 � The inverse matrix of a n × n matrix A, A -1 , The inverse matrix of a n × n matrix A, A , satisfies A ⋅ A -1 = A -1 ⋅ A = I, where I is the n × n identity matrix. identity matrix. � Not all matrices are invertible. A matrix that does not an inverse is a singular matrix and has does not an inverse is a singular matrix, and has zero determinant. � Given � Given Given A how Given A , , how how do how do do w e do w e w e find w e find find its find its its inverse? its inverse? inverse? nverse? � If you know how to solve systems of linear equations, matrix inversion is just a little more complex. 90
Matrix Inversion: 2/9 Matrix Inversion: 2/9 � Let X be a matrix such that A ⋅ X = I. Let X be a matrix such that A X I. � Consider column i of X, X i , and column i of I, I i . � We have A ⋅ X i = I i as shown below. � W h A X I h b l � Since X i is unknown and I i is known, solving for X i gives column i of matrix X. ⎡ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎤ � � ⎡ ⎤ ⎡ ⎤ 0 0 a a a a a a x x 1,1 1, 1, 1, i n i ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ � � � � � � ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ ⋅ = a a a x 1 i -th position ,1 , , , i i i i n i i ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ � � � � � � ⎢ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ � � a a a x ⎣ ⎦ 0 ⎣ ⎦ ⎣ ⎦ 91 ,1 ,1 , , n n n n n i
Matrix Inversion: 3/9 Matrix Inversion: 3/9 � Therefore, we may use the same A and solve X 1 Therefore, we may use the same A and solve X 1 from I 1 , X 2 from I 2 , …, X n from I n . � It requires the execution of a linear system solver � It requires the execution of a linear system solver n times with the same A. � W � We may use LU-decomposition to save time! LU d iti t ti ! ⎡ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎤ � � ⎡ ⎤ ⎡ ⎤ 0 0 a a a a a a x x 1,1 1, 1, 1, i n i ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ � � � � � � ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ ⋅ = a a a x 1 i -th position ,1 , , , i i i i n i i ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ � � � � � � ⎢ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ � � a a a x ⎣ ⎦ 0 ⎣ ⎦ ⎣ ⎦ 92 ,1 ,1 , , n n n n n i
Matrix Inversion: 4/9 Matrix Inversion: 4/9 � First, do a LU-decomposition A = L ⋅ U. First, do a LU decomposition A L U. � For each i from 1 to n , let I i be the column vector with a 1 in the i -th position and 0 elsewhere with a 1 in the i -th position and 0 elsewhere. � We have (L ⋅ U) ⋅ X i = I i . Applying forward and backward substitutions yields X i , the i -th column b k d b tit ti i ld X th i th l of the inverse matrix X = A -1 . � Note that if partial (complete) pivoting is used in the construction of A = L ⋅ U, one must swap rows of I k ’s before solving, and swap rows (and columns) after the solution is obtained. 93
Matrix Inversion: 5/9 Matrix Inversion: 5/9 � Consider finding the inversion of the following Consider finding the inversion of the following matrix A: ⎡ ⎡ ⎤ ⎤ 4 0 1 ⎢ ⎥ = ⎢ 3 1 3 A ⎥ ⎢ ⎢ ⎥ ⎥ ⎣ ⎣ 0 0 1 1 2 2 ⎦ ⎦ � First, find A’s LU-decomposition as follows: ⎡ ⎤ ⎢ ⎥ ⎡ ⎤ 1 4 0 1 ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ 3 3 9 9 = ⋅ = ⋅ ⎢ ⎥ ⎢ 1 ⎥ 1 A L U 4 4 ⎢ ⎥ ⎢ ⎥ − 0 0 1 1 1 1 1 1 ⎣ ⎣ ⎦ ⎦ ⎢ ⎢ ⎥ ⎥ ⎢ ⎥ ⎣ ⎦ 4 94
Matrix Inversion: 6/9 Matrix Inversion: 6/9 � Now, find the first column of the inverse. Now, find the first column of the inverse. � The right-hand side is [1,0,0] T as follows: ⎡ ⎡ ⎤ ⎤ ⎢ ⎥ ⎡ ⎤ 1 ⎡ ⎤ 4 0 1 1 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 3 9 ⎥ ⋅ ⋅ ⎥ ⋅ ⋅ = = ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ 1 1 1 1 0 0 X X ⎢ ⎢ ⎥ ⎥ 1 4 4 ⎢ ⎥ ⎢ ⎥ ⎢ ⎦ ⎥ − ⎣ 0 0 1 1 1 ⎣ ⎦ ⎢ ⎥ ⎢ ⎢ ⎥ ⎥ ⎣ ⎣ ⎦ ⎦ 4 4 � Go through a forward and a backward substitution yields: substitution yields: ⎡ ⎡ ⎤ ⎤ 1 ⎢ ⎥ = ⎢ 6 X ⎢ ⎥ ⎥ 1 ⎢ ⎥ − ⎣ 3 ⎦ 95
Matrix Inversion: 7/9 Matrix Inversion: 7/9 � Then, find the second column of the inverse. Then, find the second column of the inverse. � The right-hand side is [0,1,0] T as follows: ⎡ ⎡ ⎤ ⎤ ⎢ ⎥ ⎡ ⎤ 1 ⎡ ⎤ 4 0 1 0 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 3 9 ⎥ ⋅ ⋅ ⎥ ⋅ ⋅ = = ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ 1 1 1 1 1 1 X X ⎢ ⎢ ⎥ ⎥ 2 4 4 ⎢ ⎥ ⎢ ⎥ ⎢ ⎦ ⎥ − ⎣ 0 0 1 1 1 ⎣ ⎦ ⎢ ⎥ ⎢ ⎢ ⎥ ⎥ ⎣ ⎣ ⎦ ⎦ 4 � Go through a forward and a backward substitution yields: substitution yields: − ⎡ ⎡ ⎤ ⎤ 1 ⎢ ⎥ = − 8 X ⎢ ⎢ ⎥ ⎥ 2 ⎢ ⎥ ⎣ 4 ⎦ 96
Matrix Inversion: 8/9 Matrix Inversion: 8/9 � Finally, find the third column of the inverse. Finally, find the third column of the inverse. � The right-hand side is [0,0,1] T as follows: ⎡ ⎡ ⎤ ⎤ ⎢ ⎥ ⎡ ⎤ 1 ⎡ ⎤ 4 0 1 0 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 3 9 ⎥ ⋅ ⋅ ⎥ ⋅ ⋅ = = ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ 1 1 1 1 0 0 X X ⎢ ⎢ ⎥ ⎥ 3 4 4 ⎢ ⎥ ⎢ ⎥ ⎢ ⎦ ⎥ − ⎣ 1 0 1 1 1 ⎣ ⎦ ⎢ ⎥ ⎢ ⎢ ⎥ ⎥ ⎣ ⎣ ⎦ ⎦ 4 � Go through a forward and a backward substitution yields: substitution yields: ⎡ ⎡ ⎤ ⎤ 1 ⎢ ⎥ = ⎢ 9 X ⎢ ⎥ ⎥ 3 ⎢ ⎥ − ⎣ 4 ⎦ 97
Matrix Inversion: 9/9 Matrix Inversion: 9/9 � Therefore, the inverse of A is [X 1 | X 2 | X 3 ] : Therefore, the inverse of A is [X 1 | X 2 | X 3 ] : − ⎡ ⎤ 1 1 1 ⎢ ⎢ ⎥ ⎥ A − = − 1 1 6 6 8 8 9 9 A ⎢ ⎥ ⎢ ⎥ − − ⎣ 3 4 4 ⎦ � Let us verity it: − ⎡ ⎡ ⎤ ⎤ ⎡ ⎡ ⎤ ⎤ ⎡ ⎡ ⎤ ⎤ 4 4 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ − 1 ⋅ = ⋅ − = 3 1 3 6 8 9 0 1 0 A A ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ − − ⎣ 0 1 2 ⎦ ⎣ 3 4 4 ⎦ ⎣ 0 0 1 ⎦ 98
Determinant: 1/3 Determinant: 1/3 � The determinant of a square matrix is also easy The determinant of a square matrix is also easy to compute. � If A = L ⋅ U then the determinant of A is the � If A = L ⋅ U, then the determinant of A is the product of the diagonal elements of U. � If the construction of A = L ⋅ U uses pivoting, the � If th t ti f A L U i ti th total number of row and and column swaps matters. If the total is odd, the product should be multiplied by -1. � This is because swapping two rows (or columns) changes the sign of the determinant. g g 99
Determinant: 2/3 Determinant: 2/3 � Here is an example with complete pivoting. Here is an example with complete pivoting. 1 × (-1/6)+ 6 0 2 4 1 1 0 6 row swap column swap 0 2 4 4 2 0 4 2 0 1 1 4 4 2 2 1 1 0 0 6 6 2 2 4 4 1 1 product is 74 6 0 1 6 0 1 6 0 1 11 11 0 0 4 4 × (-1/2)+ × ( 1/2)+ 11 row swap 6 0 2 4 0 4 6 37 0 0 11 11 0 0 4 4 0 0 2 2 4 4 12 12 6 (-1) 3 74=-74 3 swaps means (-1) 3 100
Recommend
More recommend
