lu-decomposition.mw

This file demonstrates the A=LU decomposition both with and without partial pivoting.
The determinant of the A matrix can be found from the U matrix.
 

Any section/theorem/algorithm references are for "Matrix Computations" by Golub and Loan. 

 

> with( LinearAlgebra ):
 

The input matrix which we are computing the LU decomposition and the determinant for. 

> n := 3;
A := RandomMatrix(n,n);
 

We use the native Maple implementation to verify that our implementation is correct.  
   Note that the determinant is found by multiplying the diagonal values of the U matrix. (Theorem 3.2.1)
   Note also that row switching by the LU decomposition is not yet taken into account with the determinant sign.
   The evalf statement forces the pivoting to use the largest magnitude value in each column (thanks to Robert Israel, comp.soft-sys.math.maple) 

> (mP, mL, mU) := LUDecomposition( evalf(A), output=['P','L','U'] );
mDet := Determinant( A );
mDetUpp := mul( i, i=Diagonal( mU ) );
 

LU Decomposition without Partial Pivoting.
   Section 3.2.4 "Upper Triangularizing"
   Section 3.2.7 "Where is L?" 

> g_vec := Vector(n):
L := IdentityMatrix(n,n,compact=false):
U := Matrix( A ):
for k from 1 to n do

   # Without partial pivoting, we cannot handle case where leading element on diagonal is zero
   if U[k,k] = 0 then break end if:


   # Find multipliers needed to make leading value of remaining rows zero (use to subtract in next step)
   g_vec[k+1..n] := U[k+1..n,k] / U[k,k]:


  # Ignoring rows already computed, subtract an outter product from rest of matrix
   U[k+1..n,1..n] := U[k+1..n,1..n] - g_vec[k+1..n] . U[k,1..n]:

   # Each column of L is made from the multipliers for the associated row of U
   L[k+1..n,k] := g_vec[k+1..n];

end do:

(L,U);
det := mul( i, i=Diagonal( U ) );
 

Outer Product Gaussian Elimination without Partial Pivoting
   Section 3.2.6 "Some Practical Details"
   Note that both the L and U matrices are packed in the overwritten input matrix 

> oA := Matrix( A ):

for k from 1 to (n-1) do

   # Without partial pivoting, we cannot handle case where leading element on diagonal is zero
   if oA[k,k] = 0 then break end if:

   # Stash first L column multipliers
   oA[k+1..n,k] := oA[k+1..n,k] / oA[k,k];

   # Apply multipliers to submatrix
   oA[k+1..n,k+1..n] := oA[k+1..n,k+1..n] - oA[k+1..n,k] . oA[k,k+1..n];
   
end do:

oA;
 

Gauss Elimination with Partial Pivoting
   Section 3.4.3 "Partial Pivoting Details" 

> gA := Matrix( A ):
p := Vector( n-1 ):

for k from 1 to (n-1) do
   
   # Find partial pivot (note that we pulled vector norm out of loop as invariant)
   cNorm := VectorNorm( gA[k..n,k], infinity );
   for mu from k to n do
       if abs(gA[mu,k]) = cNorm then break end if:
   end do;

   # Swap row with largest pivot to top of current submatrix
   tempRow := gA[mu,k..n];
   gA[mu,k..n] := gA[k,k..n];
   gA[k,k..n] := tempRow;

   # Update permutation vector with swap we just made
   p[k] := mu;

   # Use outer product gaussian elimination on remaining submatrix
   if gA[k,k] <> 0 then

       # Stash first L column multipliers
       gA[k+1..n,k] := gA[k+1..n,k] / gA[k,k];

       # Apply multipliers to submatrix
       gA[k+1..n,k+1..n] := gA[k+1..n,k+1..n] - gA[k+1..n,k] . gA[k,k+1..n];

   end if;

end do:

gA;
p;