Linear Systems II CS3220 Summer 2008 Jonathan Kaldor Revisiting - - PowerPoint PPT Presentation

Linear Systems II

CS3220 Summer 2008 Jonathan Kaldor

Revisiting the LU Factorization

• Goal: Solve Ax=b
• Represents n linear equations, with n

unknowns

• Assume A is full rank (invertible,

nonsingular)

• Idea: Simplify without changing the solution
• What form is simpler to solve?

Triangular Matrices

Triangular Matrices

• Upper Triangular Matrix

a11 a12 .... a1n 0 a22 .... a2n . . . . . . . . . . an-1,n-1 an-1,n 0 0 .... ann x1 x2 . . xn-1 xn b1 b2 . . bn-1 bn =

Triangular Matrices

• Upper Triangular Matrix

a11 a12 .... a1n 0 a22 .... a2n . . . . . . . . . . an-1,n-1 an-1,n 0 0 .... ann x1 x2 . . xn-1 xn b1 b2 . . bn-1 bn = Solve for xn

Triangular Matrices

• Upper Triangular Matrix

a11 a12 .... a1n 0 a22 .... a2n . . . . . . . . . . an-1,n-1 an-1,n 0 0 .... ann x1 x2 . . xn-1 xn b1 b2 . . bn-1 bn = Solve for xn Substitute xn, solve for xn-1

Triangular Matrices

• Upper Triangular Matrix

a11 a12 .... a1n 0 a22 .... a2n . . . . . . . . . . an-1,n-1 an-1,n 0 0 .... ann x1 x2 . . xn-1 xn b1 b2 . . bn-1 bn = Solve for xn Substitute xn, solve for xn-1 Repeat for each variable

Triangular Matrices

• Upper Triangular Matrix
• “Backward Substitution”

a11 a12 .... a1n 0 a22 .... a2n . . . . . . . . . . an-1,n-1 an-1,n 0 0 .... ann x1 x2 . . xn-1 xn b1 b2 . . bn-1 bn = Solve for xn Substitute xn, solve for xn-1 Repeat for each variable

LU Factorization

• Convert Ax=b into MAx=Mb, where

MA is upper triangular

• M = Mn-1...M2M1
• M1 - matrix that eliminates first

variable from equations 2...n, M2 is...

• Property: If C and D are both upper (or

both lower) triangular, CD is also upper (or lower) triangular (so M is lower tri)

LU Factorization

• So after n-1 steps, have:

Mn-1...M2M1Ax = Mn-1...M2M1b

• Define U = Mn-1...M2M1A.
• Upper triangular by construction
• Multiply by L = M1-1M2-1...Mn-1-1
• Then LU = A

Gaussian Elimination

• What do the update matrices Mi look like?

1 0 ... 0 0 ... 0 0 1 ... 0 0 ... 0 ⋮ ⋮ ⋱ ⋮ ⋮ ⋱ ⋮ 0 0 ... 1 0 ... 0 0 0 ... -ai+1,i/ai,i 1 ... 0 ⋮ ⋮ ⋱ ⋮ ⋮ ⋱ ⋮ 0 0 ... -an,i/ai,i 0 ... 1 Ones along the diagonal Entries in the i-th column, below the diagonal, chosen to cancel i-th variable in equations i+1:n

Gaussian Elimination

• What do the inverse matrices Mi-1 look like?

Gaussian Elimination

• What do the inverse matrices Mi-1 look like?

1 0 ... 0 0 ... 0 0 1 ... 0 0 ... 0 ⋮ ⋮ ⋱ ⋮ ⋮ ⋱ ⋮ 0 0 ... 1 0 ... 0 0 0 ... ai+1,i/ai,i 1 ... 0 ⋮ ⋮ ⋱ ⋮ ⋮ ⋱ ⋮ 0 0 ... an,i/ai,i 0 ... 1 Note: entries below diagonal have flipped signs

The L in LU

• So each Mi-1 is lower triangular as well
• In fact, unit lower triangular (ones on

diagonal)

• How about the product M1-1M2-1...Mn-1-1?

Solving with LU factorization

• Have LUx = b
• Replace Ux with y
• Solve Ly = b (forward substitution)
• Then solve Ux = y (backward

substitution)

What does L look like?

What does L look like?

• We know: L is lower triangular, consists of

product of Mi-1 matrices

What does L look like?

• We know: L is lower triangular, consists of

product of Mi-1 matrices

• Do we need to multiply Mi-1 together to

form L?

What does L look like?

• We know: L is lower triangular, consists of

product of Mi-1 matrices

• Do we need to multiply Mi-1 together to

form L?

• No!

What does L look like?

• We know: L is lower triangular, consists of

product of Mi-1 matrices

• Do we need to multiply Mi-1 together to

form L?

• No!
• Do we need to scale/update rows like U?

What does L look like?

• We know: L is lower triangular, consists of

product of Mi-1 matrices

• Do we need to multiply Mi-1 together to

form L?

• No!
• Do we need to scale/update rows like U?
• No! (chalkboard)

A Vectorized U update

• How can we quickly apply Mi to partially

computed U matrix?

Pivoting

• Gaussian (and Gauss-Jordan) elimination

both have failure cases

• Consider the following matrix:

1 2 1

• 1 -2 0

0 1 -1

Pivoting

• Gaussian (and Gauss-Jordan) elimination

both have failure cases

• Consider the following matrix:

1 2 1

• 1 -2 0

0 1 -1 1 2 1 0 0 1 0 1 -1 Scale first equation by -(-1/1), add to second equation

Pivoting

• Gaussian (and Gauss-Jordan) elimination

both have failure cases

• Consider the following matrix:

1 2 1 0 0 1 0 1 -1 Scale first equation by -(-1/1), add to second equation

Pivoting

• Gaussian (and Gauss-Jordan) elimination

both have failure cases

• Consider the following matrix:

1 2 1 0 0 1 0 1 -1

Pivoting

• Gaussian (and Gauss-Jordan) elimination

both have failure cases

• Consider the following matrix:

1 2 1 0 0 1 0 1 -1 Scale second equation by -(2/0), add to first equation

Pivoting

• Gaussian (and Gauss-Jordan) elimination

both have failure cases

• Consider the following matrix:

1 2 1 0 0 1 0 1 -1 Scale second equation by -(2/0), add to first equation

?

Pivoting

• Gaussian (and Gauss-Jordan) elimination

both have failure cases

• Consider the following matrix:

1 2 1 0 0 1 0 1 -1 Scale second equation by -(2/0), add to first equation

?

Division by zero!

Pivoting

1 2 1 0 0 1 0 1 -1

Pivoting

• Does that mean that the system is singular?

1 2 1 0 0 1 0 1 -1

Pivoting

• Does that mean that the system is singular?
• Does the order we specify equations

matter? 1 2 1 0 0 1 0 1 -1

Pivoting

• Observe: if we swap the second and third

equations, Gaussian elimination works as expected

• We can reorder equations as-needed in
• rder to get “good” pivots (pivot: Ajj)

1 2 1 0 0 1 0 1 -1 1 2 1 0 1 -1 0 0 -1

What makes a good pivot?

• Again, back to our example. Suppose

A(2,2) = 1e-15

• Instead of dividing by zero, dividing by a very

small number

• Almost as bad
• Observation: Want to swap equations so

that diagonal element is as large as possible in absolute value

Partial Pivoting

• At i-th step
• Look at elements in i-th column below the

diagonal, find largest in absolute value (say, j-th)

• Swap i-th and j-th rows. Largest entry

now in A(i,i)

• Advantage: all multipliers in L ≤ 1 (recall,

multipliers are A(j,i)/A(i,i) )

Permutation Matrices

• Want to represent these row swaps using

matrix notation

• Suppose we want to swap rows i and j in A
• Take identity matrix
• Swap rows i and j, call this matrix P
• Then PA is the same as A with rows i and

j interchanged

Permutation Matrices

• Note: P-1 = PT
• PA : swaps rows
• AP : swaps columns
• Like Mi matrices, never form them

explicitly

• Oftentimes, don’t even need to explicitly

swap rows

Pivoting in Matrix Formulation

• Introduce Pi matrices during computation

Mn-1Pn-1...M2P2M1P1Ax = Mn-1Pn-1...M2P2M1P1b

• Job of each Pi : move row with largest

element in i-th column to the i-th row

Pivoting in Matrix Formulation

• U = Mn-1Pn-1...M2P2M1P1A
• L is still formed from Mi, but multipliers

below diagonal get permutated

• P = Pn-1...P2P1
• End up with LU = PA

An Example

1 2 2 4 4 2 4 6 4

LU with Partial Pivoting

• Factor A into PA = LU
• To solve Ax = b :
• Note: A = PTLU
• So we have PTLUx = b
• Multiply by P to get LUx = Pb

LU with Partial Pivoting

• Now have LUx = Pb
• Solve Ly = Pb for the vector y using

forward substitution

• Solve Ux = y for the vector x using

backward substitution

LU with Partial Pivoting

• What happens if all the elements in a

column below the diagonal are zero?

Solving in MATLAB

• MATLAB provides an operator to solve

linear systems

• A \ b solves equation Ax = b
• Can also specify matrix B
• Be careful: also / operator (solves AB-1)
• Solves more than just n x n systems

(more to come!)

Solving in MATLAB

• To compute LU factorization of matrix A,

use function lu(A)

• Syntax: [L, U] = lu(A)
• L is not lower triangular...

‘psychologically lower triangular’

• Or: [L,U,P] = lu(A)
• P: permutation matrix

Wait! ‘Partial’ Pivoting?!

• We reordered equations to get a good

pivot element... could we also relabel our unknowns to get a better element?

• Short answer: Yes
• Longer answer: Yes, but typically partial

pivoting is good enough

LU Factorization: A Summary

• Generalized solution method for n x n

systems of linear equations in n unknowns

• Requires pivoting in order to avoid certain

failure cases

Special Matrices

• Oftentimes, our linear system has certain

properties

• Useful to exploit these properties for

speed or space savings

Symmetric Matrices

• A matrix A is symmetric if AT = A
• Another way of putting it: Ai,j = Aj,i
• Only need to store half of the matrix
• (other half is the same)

LU Factorization of Symmetric Matrix

• (Example)

4 2 2 2 10 4 2 4 18

LU Factorization of Symmetric Matrix

• A was symmetric: ideally, there would be

some similarity between L and U

• Not immediately apparent from LU

matrices

LU Factorization of Symmetric Matrix

• What is the effect of multiplying a diagonal

matrix D times U?

• Scales the i-th row of U by Di,i
• Factor out diagonal matrix D from U so

that U is now unit upper triangular

LU Factorization of Symmetric Matrix

4 2 2 2 10 4 2 4 18

LDLT Factorization

• Rewriting U = DU’, where U’ is unit upper

triangular, shows that U’ = LT

• At least in our example
• True in general?

LDLT Factorization

• Let A be symmetric, so A = AT
• Substitute A = LU
• LU = (LU)T ⇒ LU = UTLT

⇒ U = L-1UTLT ⇒ U(LT)-1 = L-1UT

LDLT Factorization

• U(LT)-1 = L-1UT
• Note: LT is upper triangular, as is (LT)-1
• L-1 is lower triangular, as is UT
• Left hand side: upper triangular. Right hand

side: lower triangular

• Implies both sides are diagonal matrices

(why?)

LDLT Factorization

• U(LT)-1 = L-1UT
• So both sides are diagonal matrices. Call

this matrix D. Then D = L-1UT

• U(LT)-1 = L-1UT ⇒ U(LT)-1 = D

⇒ U = DLT

• So A = LU = LDLT

Symmetric Matrices

• Advantages: reduces amount of space

needed to store factorization

Positive Definite Matrices

• A matrix A is positive definite if, for all

vectors x ≠ 0, xTAx > 0

• Seems like a strange definition
• Very important (and fairly common) in

practice

Positive Definite Matrices

• Engineering applications: need to find

stresses on object with external forces

• Use finite element method
• Requires solving matrix equation, very

commonly positive definite

Positive Definite Matrices

• Simplest positive definite matrix: for any

nonsingular matrix A, ATA is positive definite

• Why?

Cholesky Factorization

• If A is symmetric and positive definite, then

all diagonal entries in D are positive

• Can take square root of D: square root
• f each diagonal element
• Leads to

LD1/2D1/2LT

(LD1/2)(LD1/2)T

(L’)(L’)T

Cholesky Factorization

• A = L’L’T
• Only works for positive definite matrices
• Note: L’ lower triangular, but not unit
• How do we know if a matrix is positive

definite?

• Run Cholesky algorithm. If it fails... not

positive definite (or close enough that it may as well not be)

Cholesky Factorization

• Can be derived directly from A = LLT
• Take 2 x 2 case:

Red: we don’t know what its value is yet A11 A12 A21 A22 L11 0 L21 L22 L11 L21 0 L22 =

Cholesky Factorization

A11 A12 A21 A22 L11 0 L21 L22 L11 L21 0 L22 = A11 = L11*L11 by rules of matrix multiplication So: L11 = sqrt(A11)

Cholesky Factorization

A11 A12 A21 A22 L11 0 L21 L22 L11 L21 0 L22 = A21 = L11 * L21, again by matrix mult. So L21 = A21 / L11 L11 = sqrt(A11)

Cholesky Factorization

A11 A12 A21 A22 L11 0 L21 L22 L11 L21 0 L22 = A22 = L21 * L21 + L22*L22 So L22 = sqrt(A22 - L21*L21) L11 = sqrt(A11) L21 = A21/L11

Cholesky Factorization

A11 A12 A21 A22 L11 0 L21 L22 L11 L21 0 L22 = L11 = sqrt(A11) L21 = A21/L11 L22 = sqrt(A22 - L21*L21)

Cholesky Example

4 2 2 2 10 4 2 4 18

Cholesky Factorization

• Find L by columns - find first column, then

second, ...

• First find diagonal entry Li,i using Ai,i and

Li,1:i-1 , which we already know

• Then find entries below diagonal, using Li,i

and entries in previous columns, both of which we already know

Cholesky Factorization

• In MATLAB: chol(A)
• Other benefits of Cholesky
• pivoting not required
• half storage space
• half of the running time of LU

Summary of Linear System Solvers

• Gauss-Jordan Elimination
• Reduce matrix to diagonal form
• Can compute inverse matrix
• Gaussian Elimination
• Reduce matrix to upper triangular form
• Faster than G-J : About 2/3n3 FLOPS

Summary of Linear System Solvers

• LU Factorization: Use Gaussian Elimination

to form two triangular systems

• Can factor without having b available
• Cholesky Factorization: Find A = LLT for

symmetric positive definite matrices

• About twice as fast as LU : 1/3n3 FLOPS
• Does not require pivoting