Numerical Methods & Scientific Computing: lecture notes
Linear Systems
LU factorization
Many software libraries contain no provision for Gauss Elimination — WHY NOT?
if you solve Ax = b1 (⇡ 1 n3 operations) and, some time later, 133
solve Ax = b2 (⇡ 3 n operations) most of these operations have been done before! (all but n2 operations)
=) try to store information so you don’t waste this e↵ort.
In practice, we store the multipliers lij in a lower triangular matrix
l43 1
lij i>j L=:1i=j 0i
usually consider only symmetric positive definite matrices
arise naturally in some applications e.g. Least Squares fitting same as having all positive eigenvalues
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
Special matrices
Choleski factorization
Theorem
A is symmetric positive definite is equivalent to A has a symmetric triangular factorization
A = RT R
where R is a nonsingular upper (right) triangular matrix.
Called the Choleski factorization of A after Choleski: †1918 RT plays the same role as L; R plays the same role as U. The Choleski factorization is unique.
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
Special matrices
Proof.
To show A = RT R =) A is positive definite: xTAx = xTRTRx = (Rx)TRx = ||Rx||2 > 0 8x 6= 0 To show A = RTR =) A is symmetric:
AT =(RTR)T =RT(RT)T =RTR=A
To show converse: by induction.
This factorization takes ⇡ 1 n3 ⇥ / ÷ +np operations 6
Turns out that Choleski factorization is not vulnerable to subtractive cancellation, so don’t need to pivot.
Moral: If you know A is symmetric positive definite, can save 50% of work by Choleski factorization.
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
Special matrices
Banded matrices
A matrix is banded if it has nonzeroes only near the main diagonal: Aij =0for |i j|>k
bandwidth of matrix = 2k + 1
i.e. a regular pattern of zeroes, a generalization of diagonal matrices.
Example
k = 1 tridiagonal matrices
k = 2 pentadiagonal matrices
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
Special matrices
Fast factorization
If no pivoting is required, the LU factors inherit this banded structure ! can be factored very fast
Example
the LU factors of a tridiagonal matrix are bidiagonal
hence a tridiagonal matrix can be factored without pivoting in ⇡ 2n ⇥ /÷ operations
VERY FAST!
Use any special structure of the matrix to save time/storage.
If pivoting is required, it destroys the banded structure, so the factors are not banded!
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
Special matrices
What does MATLAB’s backslash do?
To solve the system Ax = b, in MATLAB just type x=A\b;
to remember this, think of it as premultiplying by the inverse: x = A 1b
This backslash command (amongst other things):
solves by substitution if A is triangular
attempts Cholesky factorization if A is symmetric
does GEPP by LU factorization if A is a general dense matrix
So, easiest way to solve with LU if you want to keep the LU factors is [L,U,P]=lu(A); x=U\(L\(P*b))
similarly for Cholesky:
R=chol(A); x=R\(R’\b)
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
Special matrices
End of Week 6