Transcript LU Factorization of A

Numerical Computation
Lecture 6: Linear Systems – part II
United International College
Review Presentation
• 10 minutes for review presentation.
Review
• During our Last Two Classes we covered:
– Linear Systems: Gaussian Elimination
– Section 3.1 in Pav
– Sections 2.1, 2.2, 2.3 in Moler
Today
• We will cover:
– Linear Systems: LU Factorization (or
Decomposition)
– Section 3.3 in Pav
– Sections 2.4, 2.5 in Moler
– Linear Algebra Theorems
LU Factorization
• Today we will show that the Gaussian
Elimination method can be used to factor (or
decompose) a square, invertible, matrix A into
two parts, A = LU, where:
– L is a lower triangular matrix with 1’s on the
diagonal, and
– U is an upper triangular matrix
LU Factorization
• A= LU
• We call this a LU Factorization of A.
LU Factorization Example
• Consider the augmented matrix for Ax=b:
• Gaussian Elimination (with no row swaps)
yields:
LU Factorization Example
• Let’s Consider the first step in the process –
that of creating 0’s in first column:
• This can be achieved by multiplying the
original matrix [A|b] by M1
LU Factorization Example
• Thus, M1 [A|b] =
• In the next step of Gaussian Elimination, we
pivot on a22 to get:
• We can view this as
multiplying M1 [A|b] by
LU Factorization Example
• Thus, M2 M1 [A|b] =
• In the last step of Gaussian Elimination, we
pivot on a33 to get:
• We can view this as
multiplying M2 M1 [A|b] by
LU Factorization Example
• Thus, M3 M2 M1 [A|b] =
• So, the original problem of solving Ax=b can
now be thought of as solving
M3 M2 M1 Ax = M3 M2 M1 b
LU Factorization Example
• Note that M3 M2 M1 A is upper triangular:
• Thus,
• We hope that L is lower triangular
LU Factorization Example
• Recall: In Exercise 3.6 you proved that the
inverse to the matrix
• Is the matrix
LU Factorization Example
• Thus, L = M1-1 M2-1 M3-1
LU Factorization Example
• So, A can be factored as A=LU
Solving Ax=b using LU Factorization
• To solve Ax=b we can do the following:
– Factor A = LU
– Solve the two equations:
• Lz = b
• Ux = z
– This is the same as L(Ux) = b or Ax=b.
LU Factorization Example
• Our example was Ax=b:
• Factor:
LU Factorization Example
• Solve Lz=b using forward substitution:
• We get:
LU Factorization Example
• Solve Ux=z using back substitution:
• We get the solution vector for x:
(Note typo in Pav, page 36
bottom of page “z” should be “x”)
LU Factorization Theorem
• Theorem 3.3. If A is an nxn matrix, and
Gaussian Elimination does not encounter a
zero pivot (no row swaps), then the algorithm
described in the example above generates a
LU factorization of A, where L is a lower
triangular matrix (with 1’s on the diagonal),
and U is an upper triangular matrix.
• Proof: (omitted)
LU Factorization Matlab Function
(1 of 2)
function [ L U ] = lu_gauss(A)
% This function computes the LU factorization of A
% Assumes A is not singular and that Gauss Elimination requires no row swaps
[n,m] = size(A); % n = #rows, m = # columns
if n ~= m; error('A is not a square matrix');
end
for k = 1:n-1 % for each row (except last)
if A(k,k) == 0, error('Null diagonal element');
end
for i = k+1:n % for row i below row k
m = A(i,k)/A(k,k); % m = scalar for row i
A(i,k) = m; % Put scalar for row i in (i,k) position in A.
% We do this to store L values. Since the lower
% triangular part of A gets zeroed out in GE, we
% can use it as a storage place for values of L.
LU Factorization Matlab Function
(1 of 2)
for j = k+1:n % Subtract m*row k from row i -> row i
% We only need to do this for columns k+1 to n
% since the values below A(k,k) will be zero.
A(i,j) = A(i,j) - m*A(k,j);
end
end
end % At this point, A should be a matrix where the upper triangular
% part of A is the matrix U and the rest of A below the diagonal
% is L (but missing the 1's on the diagonal).
L = eye(n) + tril(A, -1); % eye(n) is a matrix with 1's on diagonal
U = triu(A);
end
LU Factorization Matlab
>> A=[2 1 1 3; 4 4 0 7; 6 5 4 17; 2 -1 0 7]
A=
2 1 1 3
4 4 0 7
6 5 4 17
2 -1 0 7
>> [l u] = lu_gauss(A);
LU Factorization Matlab
>> >> l
l=
1 0 0 0
2 1 0 0
3 1 1 0
1 -1 -1 1
>> u
u=
2 1 1 3
0 2 -2 1
0 0 3 7
0 0 0 12
LU Factorization Matlab Solve
Function
• Class Discussion: What must we do to modify
the function lu_gauss so that we can compute
the solution to Ax=b using the LU
factorization?
• New function: lu_solve on handouts. Discuss
and consider Matlab output ->
LU Factorization Matlab
>>> [x z] = lu_solve(A, b);
>> z
z=
7
-3
13
18
>> x
x=
1.5417
-1.4167
0.8333
1.5000
Linear Algebra Review
• Definition: Two matrices A and B are equal if
they have the same number of rows and
columns and if aij = bij .
• Definition: If A and B are nxm matrices, the
sum A + B is the nxm matrix with entries aij +
bij .
• Definition: If A is nxm and c is a real number,
the scalar multiplication cA is the nxm matrix
with entries c aij .
Linear Algebra Review
• Note: In Matlab we multiply two matrices by using
‘*’
• >> A=[1 2; 3 4];
• >> B=[1 1; 2 2];
• >> A*B
Linear Algebra Review
• Note: In Matlab we multiply two matrices by using
the ‘*’ operator:
>> A=[1 2; 3 4];
>> B=[1 1; 2 2];
>> A*B
ans =
5 5
11 11
Linear Algebra Review
Linear Algebra Review
• Note: In Matlab, we can create an identity
matrix by typing eye(n)
>> eye(4)
ans =
1 0
0 1
0 0
0 0
0
0
1
0
0
0
0
1
Linear Algebra Review
Linear Algebra Review
Linear Algebra Review
• Note: In Matlab, we use the ' symbol to find
the transpose:
>> A=[1 2; 3 4]
A=
1 2
3 4
>> A'
ans =
1 3
2 4