Lay Chapter 2,
LU Decomposition

LU Factorization with 2x2 Elementary Matrices

Use premultiplication by elementary matrices to produce the LU factorization.

In[1]:=

Clear[A, L, U] ; 

A = {{6, 9}, {4, 5}} ;

%//MatrixForm

Out[3]//MatrixForm=

( {{6, 9}, {4, 5}} )

In[4]:=

e1 = ( {{1, 0}, {-2/3, 1}} ) ; 

U = e1 . A ;

%//MatrixForm

Out[6]//MatrixForm=

( {{6, 9}, {0, -1}} )

In[7]:=

L = ( {{1, 0}, {2/3, 1}} ) ;

Check.

In[8]:=

L . UA

Out[8]=

True

LU Factorization with 3x3 Elementary Matrices

Use elementary matrices to produce the LU factorization.

In[9]:=

Clear[A, L, U] ; 

A = {{2, -4, 4, -2}, {6, -9, 7, -3}, {-1, -4, 8, 0}} ;

%//MatrixForm

Out[11]//MatrixForm=

( {{2, -4, 4, -2}, {6, -9, 7, -3}, {-1, -4, 8, 0}} )

In[12]:=

e1 = ({{1, 0, 0}, {-3, 1, 0}, {1/2, 0, 1}}) ;

e1 . A//MatrixForm

Out[13]//MatrixForm=

( {{2, -4, 4, -2}, {0, 3, -5, 3}, {0, -6, 10, -1}} )

In[14]:=

e2 = ({{1, 0, 0}, {0, 1, 0}, {0, 2, 1}}) ;

e2 . e1 . A//MatrixForm

Out[15]//MatrixForm=

( {{2, -4, 4, -2}, {0, 3, -5, 3}, {0, 0, 0, 5}} )

In[16]:=

U = e2 . e1 . A ;

L = ({{1, 0, 0}, {3, 1, 0}, {-1/2, -2, 1}}) ;

Check.

In[18]:=

L . UA

Out[18]=

True

Solving Ax=b with LU Factorization

Construct the matrices L and U.

In[19]:=

Clear[A, L, U, b, y] ; 

A = {{3, -7, -2}, {-3, 5, 1}, {6, -4, 0}} ;

%//MatrixForm<br />

b = {{-7}, {5}, {2}} ;

Out[21]//MatrixForm=

( {{3, -7, -2}, {-3, 5, 1}, {6, -4, 0}} )

In[23]:=

e1 = ({{1, 0, 0}, {1, 1, 0}, {-2, 0, 1}}) ;

e1 . A ;

%//MatrixForm

Out[25]//MatrixForm=

( {{3, -7, -2}, {0, -2, -1}, {0, 10, 4}} )

In[26]:=

e2 = ({{1, 0, 0}, {0, 1, 0}, {0, 5, 1}}) ;

U = e2 . e1 . A ;

%//MatrixForm

Out[28]//MatrixForm=

( {{3, -7, -2}, {0, -2, -1}, {0, 0, -1}} )

In[29]:=

L = ({{1, 0, 0}, {-1, 1, 0}, {2, -5, 1}}) ;

Check.

In[30]:=

AL . U

Out[30]=

True

Solve A.x==b by LU Factorization

In[31]:=

(* Axb = > LUxb = > Lyb and Uxy *)

In[32]:=

y = {y1, y2, y3} ; 

ySoln = Solve[L . y (b//Flatten), y]//Flatten

Out[33]=

{y1 -7, y2 -2, y36}

In[34]:=

x = {x1, x2, x3} ; 

xSoln = Solve[U . xy/.ySoln, x]//Flatten

Out[35]=

{x13, x24, x3 -6}

In[36]:=

ans = x/.xSoln

Out[36]=

{3, 4, -6}

Solve A.x==b by row reduction.

In[37]:=

<<LinearAlgebra`MatrixManipulation`

In[38]:=

b = {{-7}, {5}, {2}} ; 

aug = AppendRows[A, b] ;

%//MatrixForm

Out[40]//MatrixForm=

( {{3, -7, -2, -7}, {-3, 5, 1, 5}, {6, -4, 0, 2}} )

In[41]:=

rr = RowReduce[aug] ;

%//MatrixForm

Out[42]//MatrixForm=

( {{1, 0, 0, 3}, {0, 1, 0, 4}, {0, 0, 1, -6}} )

In[43]:=

ans2 = rr[[All, 4]]

Out[43]=

{3, 4, -6}

In[44]:=

ansans2

Out[44]=

True

makeMatrix

Let's write a procedure which creates matrices.

In[45]:=

makeMatrix[n_, m_] := Table[Random[Integer, {-9, 9}], {n}, {m}]

Generate a few matrices.
New matrices will be generated each time these cells are evaluated.

In[46]:=

Clear[a, b, c] ; 

a = makeMatrix[2, 2] ;

%//MatrixForm

Out[48]//MatrixForm=

( {{6, -9}, {-3, -2}} )

In[49]:=

b = makeMatrix[3, 4] ;

%//MatrixForm

Out[50]//MatrixForm=

( {{5, -3, -8, 1}, {-5, -1, -8, -2}, {6, -3, -8, 0}} )

In[51]:=

c = makeMatrix[4, 1] ;

%//MatrixForm

Out[52]//MatrixForm=

( {{-5}, {-5}, {2}, {-6}} )

Solving a 3x3 System of Equations with the LU Decomposition

We wish to solve the linear system represented by the following augmented matrix.

In[53]:=

Clear[a, m, b, lud] ; 

a = makeMatrix[3, 4] ;

%//MatrixForm

Out[55]//MatrixForm=

( {{-1, 6, -4, 9}, {-2, 8, 2, 9}, {6, -4, 1, 8}} )

Here is an efficient way to do that -- using the LU decomposition of a square matrix m.
First, extract the relevant components, the matrix m and the vector b.

In[56]:=

<<LinearAlgebra`MatrixManipulation`

In[57]:=

m = TakeColumns[a, 3] ;

%//MatrixForm

Out[58]//MatrixForm=

( {{-1, 6, -4}, {-2, 8, 2}, {6, -4, 1}} )

In[59]:=

b = TakeColumns[a, -1] ;

%//MatrixForm

Out[60]//MatrixForm=

( {{9}, {9}, {8}} )

Now use the LU decomposition of the square matrix m to solve the original system.

In[61]:=

lud = LUDecomposition[m] ;

Map[MatrixForm, lud]

Out[62]=

{( {{-1, 6, -4}, {2, -4, 10}, {-6, -8, 57}} ), ( {{1}, {2}, {3}} ), 1}

In[63]:=

x = LUBackSubstitution[lud, b] ;

%//MatrixForm

Out[64]//MatrixForm=

( {{293/114}, {413/228}, {-10/57}} )

Check.

In[65]:=

m . xb

Out[65]=

True

Mathematica's LUDecomposition

LUDecomposition returns a triple containing a matrix, a vector, and a condition number.

The matrix returned by LUDecomposition contains both L and U in it.
The nonzero entries of U are on and above its diagonal, and
the entries below its diagonal are the corresponding entries of L.

The procedures Upper and Lower, which are defined below, have been modified from those given in the Mathematica Help Browser.

In[66]:=

Upper[LU_ ? MatrixQ] := LU Table[If[i≤j, 1, 0], {i, Length[LU]}, {j, Length[LU]}]

Lower[LU_ ? MatrixQ] := LU - Upper[LU] + IdentityMatrix[Length[LU]] 

LU = lud[[1]] ; 

{L = Lower[LU], U = Upper[LU]} ;

Map[MatrixForm, %]

Out[70]=

{( {{1, 0, 0}, {2, 1, 0}, {-6, -8, 1}} ), ( {{-1, 6, -4}, {0, -4, 10}, {0, 0, 57}} )}

The vector returned by LUDecomposition indicates the permutation of the rows of m that will bring it into conformity with L.U.

In[71]:=

P = lud[[2]] ; 

Part[m, P] - Lower[LU] . Upper[LU] ;

%//MatrixForm

Out[73]//MatrixForm=

( {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}} )


Created by Mathematica  (February 10, 2005) Valid XHTML 1.1!