Re: Reduced Echelon Form...

[Dirk Laurie]

> On 12-Nov-1997, Ryan Kirkpatrick <address@hidden> 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