Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Aug 2016
Lecture topic: LU Decomposition

LU Decomposition: Brute Force

The process of LU decomposition attempts to turn a matrix into the product of a lower triangular matrix and an upper triangular matrix. The procedure is outlined here.

Brute Force Approach Leads to (n^2) Equations in (n^2)+n Unknowns

Here is a brute force attempt to decompose a matrix into upper and lower parts:

Clear[i,j]
n=4;
A=Array[a,{n,n}];
MatrixForm[A]

I hate this new command for defining a lower diagional matrix, but it works:

L =Normal[SparseArray[{i_,j_}/;i≥j->l[i,j],{n,n}]];
MatrixForm[L]

U =Normal[SparseArray[{i_,j_}/;i≤j->u[i,j],{n,n}]];
MatrixForm[U]

LU=L.U;
MatrixForm[%]

Now, just use brute force substitution to see what the relationship must be between elements of the matrices.

eqlist={};
varlist={};
Do[eqlist=Append[eqlist,LU[[i,j]]==A[[i,j]]],{i,1,n},{j,1,n}]
Do[varlist=Append[varlist,L[[j,i]]],{i,1,n},{j,i,n}]
Do[varlist=Append[varlist,U[[i,j]]],{i,1,n},{j,i,n}]

TableForm[eqlist]
Length[eqlist]
Length[varlist]

So we have (n^2) equations in (n^2)+n unknowns. If we impose different additional n conditions to get (n^2)+n equations in (n^2)+n unknowns, we will get different methods of determining an LU decomposition. Choosing wisely will lead to simple systems that are solvable by backsubstitution, which is somewhat surprising since the system of equations is nonlinear.

One's on Diagonal of Lower Triangular Matrix (Doolittle's Method)

Here is Doolittle's method to LU decompose a matrix.

n=4;
A={{1,4,5,-3},{4,5,6,-9},{8,7,1/3,9},{-1,1,-1,1}};
MatrixForm[A]

Here we have added the additional n conditions that the lower diagonal matrix has 1's on the diagonal.

L =Normal[SparseArray[{i_,j_}/;i≥j->l[i,j],{n,n}]];
Do[L[[i,i]]=1,{i,1,n}]
MatrixForm[L]

U = Normal[SparseArray[{i_,j_}/;i≤j->u[i,j],{n,n}]];
MatrixForm[U]

LU=L.U;
eqlist={};
varlist={};
Do[eqlist=Append[eqlist,LU[[i,j]]==A[[i,j]]],{i,1,n},{j,1,n}]
Do[varlist=Append[varlist,L[[j,i]]],{i,1,n},{j,i+1,n}]
Do[varlist=Append[varlist,U[[i,j]]],{i,1,n},{j,i,n}]

TableForm[eqlist]
Length[eqlist]
Length[varlist]

So we have (n^2) equations in (n^2) unknowns. These equations can be solved by back-substitution.

sol=Solve[eqlist,varlist]

Here we see that we have successfully factored the original matrix:

MatrixForm[L/.Flatten[sol]]
MatrixForm[U/.Flatten[sol]]
MatrixForm[L.U/.Flatten[sol]]

One's on Diagonal of Upper Triangular Matrix (Crout's Method)

Here is Crout's Method to LU decompose a Hilbert Matrix:

n=4;
A={{1,4,5,-3},{4,5,6,-9},{8,7,1/3,9},{-1,1,-1,1}};
MatrixForm[A]

Here we have added the additional n conditions that the upper diagonal matrix has 1's on the diagonal.

L =Normal[SparseArray[{i_,j_}/;i≥j->l[i,j],{n,n}]];
MatrixForm[L]

U =Normal[SparseArray[{i_,j_}/;i≤j->u[i,j],{n,n}]];
Do[U[[i,i]]=1,{i,1,n}]
MatrixForm[U]

LU=L.U;
eqlist={};
varlist={};
Do[eqlist=Append[eqlist,LU[[i,j]]==A[[i,j]]],{i,1,n},{j,1,n}]
Do[varlist=Append[varlist,L[[j,i]]],{i,1,n},{j,i,n}]
Do[varlist=Append[varlist,U[[i,j]]],{i,1,n},{j,i+1,n}]

TableForm[eqlist]
Length[eqlist]
Length[varlist]

So we have (n^2) equations in (n^2) unknowns, which again can be solveds with back-substitution.

sol=Solve[eqlist,varlist]

Here we see that we have successfully factored the original matrix:

MatrixForm[L/.Flatten[sol]]
MatrixForm[U/.Flatten[sol]]
MatrixForm[L.U/.Flatten[sol]]

(Subscript[l, ii])=(Subscript[u, ii]) (Choleski's Method)

Here is Choleski's method to LU decompose a matrix

n=4;
A={{1,4,5,-3},{4,5,6,-9},{8,7,1/3,9},{-1,1,-1,1}};
MatrixForm[A]

Here we have added the additional n conditions that the lower diagonal elements equal the upper diagonal.

L =Normal[SparseArray[{i_,j_}/;i≥j->l[i,j],{n,n}]];
MatrixForm[L]
U = Normal[SparseArray[{i_,j_}/;i≤j->u[i,j],{n,n}]];
Do[U[[i,i]]=L[[i,i]],{i,1,n}]
MatrixForm[U]

LU=L.U;
eqlist={};
varlist={};
Do[eqlist=Append[eqlist,LU[[i,j]]==A[[i,j]]],{i,1,n},{j,1,n}]
Do[varlist=Append[varlist,L[[j,i]]],{i,1,n},{j,i,n}]
Do[varlist=Append[varlist,U[[i,j]]],{i,1,n},{j,i+1,n}]

TableForm[eqlist]
varlist
Length[eqlist]
Length[varlist]

So we have (n^2) equations in (n^2) unknowns, again solvable though back-substitution.

sol=Solve[eqlist,varlist];

Here we see that we have successfully factored the original matrix:

MatrixForm[L/.Flatten[sol]]
MatrixForm[U/.Flatten[sol]]
MatrixForm[L.U/.Flatten[sol]]

Notice this decomposition resulted in a compex valued matrix, since the system of equations to solve was nonlinear.

LU Decomposition: Doolittle’s Method with Matrices

In General

Let’s take a closer look at the LU decomposition process for a general 3×3 matrix.

Clear[M,m,L,l,U,u]
n=3;
A=Array[a,{n,n}];
MatrixForm[A]

Here is the first M matrix.

M1=Normal[SparseArray[{i_,j_}/;i≥j->m1[i,j],{n,n}]];
(*assign ones to diagonal*)
Do[M1[[i,i]]=1,{i,1,n}]
(*assign nonzero values to lower column in row 1*)
k=1;
Do[M1[[j,k]]=-A[[j,k]]/A[[k,k]],{j,k+1,n}]
(*assign zero to lower rest of columns*)
Do[
If[kk≠k,Do[M1[[j,kk]]=0,{j,kk+1,n}]]
,{kk,1,n}]
MatrixForm[M1]

Notice that it has the proper inverse.

Inverse[M1]//MatrixForm

Now the problem is that to form the second matrix M2, we need to have completed the first step of the Guassian elimination on A. But this is exactly what the matrix M1 does!

A2=M1.A;
MatrixForm[A2]

M2=Normal[SparseArray[{i_,j_}/;i≥j->m2[i,j],{n,n}]];
(*assign ones to diagonal*)
Do[M2[[i,i]]=1,{i,1,n}]
(*assign nonzero values to lower column in row 1*)
k=2;
Do[M2[[j,k]]=-A2[[j,k]]/A2[[k,k]],{j,k+1,n}]
(*assign zero to lower rest of columns*)
Do[
If[kk≠k,Do[M2[[j,kk]]=0,{j,kk+1,n}]]
,{kk,1,n}]
MatrixForm[M2]

A3=M2.A2;
MatrixForm[A3]//Simplify

Notice this is an upper triangular matrix!

Specific Example

Let’s take a closer look at the LU decomposition process for a specific matrix.

Clear[M,m,L,l,U,u]
n=4;
A={{1,4,5,-3},{4,5,6,-9},{8,7,1/3,9},{-1,1,-1,1}};
MatrixForm[A]

Here is the first matrix.

M1=Normal[SparseArray[{i_,j_}/;i≥j->m1[i,j],{n,n}]];
(*assign ones to diagonal*)
Do[M1[[i,i]]=1,{i,1,n}]
(*assign nonzero values to lower column in row 1*)
k=1;
Do[M1[[j,k]]=-A[[j,k]]/A[[k,k]],{j,k+1,n}]
(*assign zero to lower rest of columns*)
Do[
If[kk≠k,Do[M1[[j,kk]]=0,{j,kk+1,n}]]
,{kk,1,n}]
MatrixForm[M1]

Notice A2 form:

A2=M1.A;
MatrixForm[A2]

Compute the second matrix M2

M2=Normal[SparseArray[{i_,j_}/;i≥j->m2[i,j],{n,n}]];
(*assign ones to diagonal*)
Do[M2[[i,i]]=1,{i,1,n}]
(*assign nonzero values to lower column in row 1*)
k=2;
Do[M2[[j,k]]=-A2[[j,k]]/A2[[k,k]],{j,k+1,n}]
(*assign zero to lower rest of columns*)
Do[
If[kk≠k,Do[M2[[j,kk]]=0,{j,kk+1,n}]]
,{kk,1,n}]
MatrixForm[M2]

Notice form of A3

A3=M2.A2;
MatrixForm[A3]//Simplify

Compute the third matrix M3

M3=Normal[SparseArray[{i_,j_}/;i≥j->m3[i,j],{n,n}]];
(*assign ones to diagonal*)
Do[M3[[i,i]]=1,{i,1,n}]
(*assign nonzero values to lower column in row 1*)
k=3;
Do[M3[[j,k]]=-A3[[j,k]]/A3[[k,k]],{j,k+1,n}]
(*assign zero to lower rest of columns*)
Do[
If[kk≠k,Do[M3[[j,kk]]=0,{j,kk+1,n}]]
,{kk,1,n}]
MatrixForm[M3]

A4=M3.A3;
MatrixForm[A4]//Simplify

Notice this is an upper triangular matrix! And is exactly the same one we found when we did Doolittle’s method by brute force above (reproduced below):

The lower triangular matrix we get from the inverse matrices (which just don’t have the minus sign, so are easy to track). Since I wasn’t keeping track of the inverses above, I will use Mathematica to compute the inverses.

Inverse[M3].Inverse[M2].Inverse[M1]//MatrixForm

Notice that if we track the inverses, we build up the matrix L row by row as we compute the matrices M

Inverse[M1]//MatrixForm
Inverse[M2].Inverse[M1]//MatrixForm
Inverse[M3].Inverse[M2].Inverse[M1]//MatrixForm

LU Decomposition (Doolittle's) Module to do the LU Decomposition

Here is a module we can create that implements this process. Notice it is fairly simple to code!

MatLUDoo[A0_]:=
Module[ {A=A0,p,k,i,m,U,L,n=Length[A0]}, 
(*
Input matrix A0;
Output two matrices, the LU decomposition
*)
(*initialize the LU matrices*)
U=A0; 
L=IdentityMatrix[n];
(*end iniitialization*)
(*do the Guassian elimination, and at each step make the changes to compute L and U*) 
For[ p=1, p≤n-1, p=p+1, 
For[ i=p+1, i≤n, i=i+1, 
(*matrix m is created at each step, update L and U*)
(*notice the lack of minus sign, this this is the inverse matrix which is used to create L*)
m=((A[[i,p]])/(A[[p,p]]));
L[[i,p]]=m; 
A[[i]] = A[[i]] - m A[[p]];
U = A; ]
 ];
(*All steps completed, L U have their final values*)
Return[{L,U}];]

L = First[MatLUDoo[A]]//MatrixForm
U= Last[MatLUDoo[A]]//MatrixForm

LU Module for Solving Linear System

Here is a module that takes input as the matrix A and vector b, and returns the solution to the associated system of linear equations.

SolveSystemLUDoo[A0_,b0_]:=
Module[ {A=A0,b=b0,p,k,i,m,U,L,n=Length[A0],x,y,f}, 
U=A0; 
L=IdentityMatrix[n]; 
For[ p=1, p≤n-1, p=p+1, 
For[ i=p+1, i≤n, i=i+1, 
m=((A[[i,p]])/(A[[p,p]]));
L[[i,p]]=m; 
A[[i]] = A[[i]] - m A[[p]];
U = A; ]
 ];

y=Array[f,n];
Do[y[[i]]=((b[[i]]-(∑_(j=1)^(i-1))y[[j]]L[[i,j]])/(L[[i,i]])),{i,1,n}];
x=Array[f,n];
Do[x[[i]]=((y[[i]]-(∑_(j=i+1)^n)x[[j]]U[[i,j]])/(U[[i,i]])),{i,n,1,-1}];
Return[x];]

Let's use it to solve the linear system of equations. First, solve L y = b for y via backsubstitution:

SolveSystemLUDoo[A,b]

Check our result:

Table[Sum[A[[j,i]]SolveSystemLUDoo[A,b][[i]],{i,1,n}]==b[[j]],{j,1,n}]

This module can be used to replace Mathematica's built-in Solve command if the equations are linear.

x={x1,x2,x3,x4};
eqs=Table[Sum[A[[j,i]]x[[i]],{i,1,n}]==b[[j]],{j,1,n}]
Solve[eqs,x]

which is really just using the Mathematica's built-in LinearSolve command:

LinearSolve[A,b]