From help-octave-request at bevo dot che dot wisc dot edu  Thu Nov 13 04:00:17 1997
Subject: Re: Reduced Echelon Form...
From: Dirk Laurie <dirk at calvyn dot puk dot ac dot za>
To: jwe at bevo dot che dot wisc dot edu (John W. Eaton)
Date: Thu, 13 Nov 1997 08:21:59 +0200 (SAT)

> On 12-Nov-1997, Ryan Kirkpatrick <rkirkpat at nag dot cs dot colorado dot edu> wrote:
> 
> | 	I have a quick question for all of you Octave people... Is there
> | a function / way to get the reduced echelon form of a matrix using Octave?

What do you mean by "the" reduced echelon form?

In general you get from a matrix to echelon form by a sequence of
  A(i,:) -= e(i,j)*A(j,:)
and
  A([i,j],:) = A([j,i],:)
operations.  In infinite precision you could in principle make this
unique by swapping only when necessary, and using the smallest available
index, but in floating point it is customary to swap so abs(A(i,i))
becomes as large as possible.  Then what about ties, etc.?

If all you need is _some_ echelon form, it can easily be done with existing
functions.

>> % Take as example a 5x4 matrix of rank 3
>> A=rand(5,3)*rand(3,4)
A =

  0.57946  0.47497  0.28647  0.69047
  1.14007  0.58709  0.55125  1.29656
  1.04450  0.68101  0.59835  1.35022
  0.50407  0.24968  0.39805  0.81140
  0.91741  0.42059  0.46076  1.06360
>> % Make it square by adding zeros and call lu
>> [L,U,P]=lu([A zeros(5,1)]);
>> U
U =

   1.14007   0.58709   0.55125   1.29656   0.00000
   0.00000   0.17657   0.00629   0.03147   0.00000
   0.00000   0.00000   0.15467   0.23990   0.00000
   0.00000   0.00000   0.00000  -0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000
>> % Diagnose which rows of U are zero
>> sum(abs(U'))
ans =

  3.57497  0.21433  0.39457  0.00000  0.00000
>> echelon=U(ans>eps,1:4)
echelon =

  1.14007  0.58709  0.55125  1.29656
  0.00000  0.17657  0.00629  0.03147
  0.00000  0.00000  0.15467  0.23990

You can see from this that the final echelon form depends on tricy
decisions whether a floating point number is zero.

Dirk