LU Decomposition

1.4.1 LU Decomposition In section 1.3.6, we have seen that calculating the inverse of a matrix A involves solving N linear problems: ~ [k] A a = e[k] (1.4.1-1) ~ [1] ~ [N] For the N column vectors a , … , a of A-1. Each of these N problems involve the same matrix A, so the Gaussian elimination procedure for each is the same. Often in the practice of numerical calculations for more advanced problems, we have to repeatedly solve a linear system with the same matrix A, but with different right hand side (RHS) vectors b[1], b[2],… Ax[k] = b[k] (1.4.1-2) LU decomposition is a technique that allows us to “remember” all of the row eliminations that we must perform to solve a linear problem with the matrix A. In the future, we therefore need only perform the N2 operations required for substitution to solve the system with a new RHS vector b[k]. Let us consider LU decomposition 1st for Gaussian elimination without partial pivoting. Consider the system Ax = b with the augmented matrix a 11 a  21 (A, b) = a 31   : a N1   a 22 ... a 2N b 2  a 32 ... a 3N b 3   : : :  a N2 ... a NN b N  a 12 ... a 1N b1 (1.4.1-3) We now perform the row operation a (2,1) ← a 2j − λ 21a 1j , 2j λ 21 = a 21 , b 2 ← b 2 − λ 21 b1 a 11 (1.4.1-4) to obtain the equivalent system a 11 a 12 ... a 1N b1   (2, 1) (2, 1) (2, 1)  0 a 22 ... a 2N b 2  (A(2,1), b(2, 1)) = a 31 a 32 ... a 3N b 3    : : :  :  a a ... a  NN b N   N1 N2 (1.4.1-5) We know from our choice of λ 21 that the (2, 1) position of A(2, 1) must contain a zero. Therefore, we are free to use the location in memory “wasted” in storing the known value of zero to store something else – for example the valyue of λ 21 . In this case, the contents stored in memory will actually be a 11 a 12 ... a 1N b1   (2, 1) (2, 1) (2, 1)  λ 21 a 22 ... a 2N b 2   (1.4.1-6) (A(2,1), b(2, 1))= a 31 a 32 ... a 3N b 3   : : :  :  a a ... a  NN b N  N1 N2  Next, we perform a row operation to “zero” the (3, 1) position. The operations required are a λ 31 = 31 , a 3j(3,1) ← a 3j − λ 31a 1j (j = 1, 2, ..., N), a 11 b 3 ← b 3 − λ 31 b1 (1.4.1-7) (3,1) Since we know that a 31 =0 by the coice of λ 31 , we can use this position in memory to store λ 31 . The contents of memory are then a 11  λ 21 (3,1) (3, 1) (A , b ) = λ 31   :  a  N1   a ... a b  (3, 1) (3, 1) a 32 ... a 3N b 3(3, 1)   : : :   a N2 ... a NN b N   a 12 ... a 1N b1 (2, 1) 22 (2, 1) 2N (2, 1) 2 (1.4.1-8) We then continue this process for the remainder of the Gaussian elimination algorithm. At the end, our augmented matrix and all of the λ ij factors are stored in memory as a 11 a 12 a 13 ... a 1N b1    (2, 1) (2, 1) (2, 1) (2, 1) λ 21 a 22 a 23 ... a 2N b 2  (N,N-1) (N, N-1) (3, 2) (3, 1) (3, 1)   (A ,b ) = λ 31 λ 32 a 33 ... a 3N b 3   : : :  :   (N, N -1) (N, N -1)  b NN  λ N1 λ N2 λ N3 ... a NN Because we now store all of the (1.4.1-9) λ ij factors, we have perfect memory of the row operations that are required to solve Ax[k] = b[k] with new RHS vectors b[k]. This saves us the work of performing Gaussian elimination each time. We 1st state how this knowledge of the λ ij ’s can be used to save the work of Gaussian elimination for repeated problems with the same matrix A. Then, we demonstrate that this approach is valid. From the final augmented matrix obtained by Gaussian elimination, let us extract the lower and upper triangular parts into separate matrices.  1  λ 1   21  U= L = λ 31 λ 32 1   :   : λ N1 λ N2 λ N3 ...1  a 11 a 12  1) a (2, 22       a 13 ... a 1N 1) (2, 1) a (2, 23 ... a 2N (3, 2) (3, 2) a 33 ... a 3N : : : N -1) a (N, NN         (1.4.1-10) Where the Gaussian elimination has been performed without pivoting. Then, LU Decomposition of A yields A = LU (1.4.1-11) This will be demonstrated shortly, but 1st let us see how we may use LU decomposition to avoid repeated Gaussian eliminations when solving Ax[k] = b[k]. We substitute A = LU in this problem to obtain Ax[k] = LU x[k] = b[k] (1.4.1-12) If we define Ux[k] = c[k], we can solve Ax[k] = b[k] by solving successively the 2 triangular problems Lc[k] = b[k] Ux[k] = c[k] (1.4.1-13) The 1st problem involves a lower triangular matrix L, λ 11 λ λ  21 22 λ 31 λ 32 λ 33  :  : λ N1 λ N2 λ N3 ...λ NN          c 1  c   2 c 3  =   ::  c N     b1  b   2 b 3    ::  b N    (1.4.1-14) This problem can be solved in N2 FLOP’s using forward substitution, L11c1 = b1 => c1 = b1/L11 L21c1 + L22c2 = b2 => c2 = (b2-L21c1)/L22 L31c1 + L32c2 + L33c3 = b3 => c3 = (b3 – L32c2 – L31c1)/L33 : : Note that L11=L22=L33 = 1, so there is no problem here with division by zero. The 2nd problem Ux[k] = c[k] involves an upper-triangular matrix and so may be solves in N2 FLOPs’ using backward substitution. Therefore, after we have performed LU decomposition on A in N3 FLOP’s once, we can solve successive problems Ax[k] = b[k] on only 2N2 < < N3 FLOP’s. LU decomposition (factorization) is a very useful and efficient technique that is used very often in practice. Now, let us show that the LU decomposition with L and U defined by (1.4.1-10) does in fact satisfy A = LU. To do so, let us define the matrix M as M = LU.  1  λ 1   21  M = λ 31 λ 32 1   :   : λ N1 λ N2 λ N3 ...1  a 11 a 12  1) a (2, 22       a 13 ... a (2, 1) 23 a a 1N ... a (3, 2) 33 (2, 1) 2N (3, 2) ... a 3N : : : N -1) a (N, NN         (1.4.1-16) By the rules of matrix multiplication, we have M11 = a11 M12 = a12 ….. M1N = a1N (1.4.1-17) So the 1st row of M equals the 1st row of A. For the 2nd row, M21 = λ 21 a11 (1.4.1-18) But since λ 21 = a21/a11, M21 = a21 (1.4.1-19) Next, for the (2, 2) position, using (1.4.1-4) 1) M22 = λ 21 a12 + a (2, = λ 21 a12 + [a22 - λ 21 a12] = a22 22 (1.4.1-20) And similarly for j = 3, …, N 1) M2j = λ 21 a1j + a (2, 2j = λ 21 a1j + [a2j - λ 21 a1j] = a2j (1.4.1-21) Therefore, the 2nd row of M equals the 2nd row of A. Continue with 3rd row, M31 = λ 31 a11 = (a31/a11)a11 = a31 1) M32 = λ 31 a12 + λ 32 a (2, 22 (1.4.1-22) (1.4.1-23) With (3, 1) 1) (3, 1) 1) λ 32 = a 32 / a (3, = a 32 / a (2, 22 22 (3, 1) (2, 1) (2, 1) and a 32 = a 32 - λ 31 a 12 (1.4.1-24) (1.4.1-25) we have (3, 1) 1) (2, 1) (3, 1) M32 = λ 31 a12 + [ a 32 / a (2, = λ 31 a12+ a 32 22 ] a 22 (2, 1) (2, 1) = λ 31 a12 + [ a 32 - λ 31 a 12 ] = λ 31 a12 +[a32 - λ 31 a12] M32 = a32 (1.4.1-26) And for j = 3, …, N 1) M3j = λ 31 a1j + λ 32 a (2, + a 3j(3, 2) 2j (1.4.1-27) Using 1) 1) a 3j(3, 2) = a 3j(3, 1) - λ 32 a (3, = a 3j(3, 1) - λ 32 a (2, 2j 2j a 3j(3, 1) = a 3j(2, 1) - λ 31 a 1j(2, 1) = a3j - λ 31 a1j 1) we have a 3j(3, 2) = a3j - λ 31 a1j - λ 32 a (2, 2j (1.4.1-28) (1.4.1-29) (1.4.1-30) therefore, 1) M3j = λ 31 a1j + λ 32 a (2, + a 3j(3, 2) 2j (1.4.1-31) 1) 1) = λ 31 a1j + λ 32 a (2, + [a3j - λ 31 a1j - λ 32 a (2, 2j 2j ] M3j = a3j The 3rd rows of M and A are equal. We can continue this process to demonstrate that M = LU satisfies M = A. 1) As an example, let us consider the system of a3j - λ 31 a1j - λ 32 a (2, 2j (1.2.1-15), page (1.2.1-6), 1 1 1  x 1  4 2 1 3  x  = 7    2    3 1 6   x 3  2 (1.4.1-32) We performed Gaussian elimination without pivoting on this system, and write the results as the L and U matrices 1   L = 2 1  3 2 1 1 1 1 U =  - 1 1  1 (1.4.1-33) Consult pages 1.2.1-6 to 1.2.1-11 to see the calculation performed during Gaussian elimination for this system. Then if M = LU, we have M11 = 1 M21 = (2)(1) = 2 M31 = (3)(1) = 3 M12 = (1)(1) = 1 M13 = (1)(1) = 1 M22 = (2)(1) + (1)(-1) = 1 M23 = (2)(1) + (1)(1) = 3 M32 = (3)(1) + (2)(-1) = 1 M33 = (3)(1) + (2)(1) +(1)(1) = 6 (1.4.1-34) Then, 1 1 1 M = LU = 2 1 3 . So we see that for this example A = LU. 3 1 6  (1.4.1-35) What form does an LU decomposition take when pivoting is performed? The book keeping is more complex, but the routine for performing LU decomposition will return a lower-triangular matrix L, an upper-triangular matrix U, and a permutation matrix P such that PA = LU (1.4.1-36) Then to solve Ax[k] = b[k], we pre-multiply by P, P Ax[k] = Pb[k] (1.4.1-37) Substitute LU = PA, Lux[k] = Pb[k] (1.4.1-38) We now solve the successive triangular problems Lc[k] = Pb[k] Ux[k] = c[k] (1.4.1-39) By permutation matrix, we mean only matrix that is obtained from the identity matrix by a sequence of row or column interchanges. 1 0 0 For example, the permutation matrix P = 0 0 1 0 1 0 1 0 0  v1   v1  Pv = 0 0 1  v 2  =  v 3  0 1 0  v 3   v 2  has the effect (1.4.1-40) (1.4.1-41) P, obtained from I by interchanging row #2 and 3, therefore acts on a vector v by interchanging the 2nd and 3rd components. Since P performs a permutation of v, it is called a permutation matrix. Note: As det(I) = 1, for any permutation matrix P, det(P) = ± 1 (1.4.1-42)