> On 12-Nov-1997, Ryan Kirkpatrick <

[hidden email]> 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