Re: reduced row echelon form

[Paul Kienzle]

I did some speed comparisons between your function and mine, and in the
process, discovered some problems:

1) does not handle empty matrices

        octave> rref([])
        error: matrix dimension mismatch solution of linear equations

2) does not handle square matrices.  Perhaps the final test should be 
   if length(j)<m rather than if length(k)<m ?

        octave> rref(eye(3))
        ans =

          1  0  0
          0  1  0
          0  0  0

3) k is wrong for [R,k] = rref(...). Should return j instead of k?

        octave> [R,k] = rref2([1 1 2; 0 0 1])
        R =

          1  1  0
          0  0  1

        k = 3


4) sometimes gives singular matrix warning and returns the wrong matrix

        octave> R = rref([1 1 2; 0 eps 1])
        warning: matrix singular to machine precision, rcond = 1.11022e-16
        R =

          0.50000  0.50000  1.00000
          0.50000  0.50000  1.00000

5) Fails when find(U(i,:)) is empty

        octave> [R,k] = rref2([1 1 2; 1 1 2; 1 1 4])
        error: single index only valid for row or column vector

Note that your function is up to 5x faster for square A, but can be
much, much slower for non-square A since it builds a square matrix
in the process of solving.  I don't know how rref is usually used, so
I don't know which approach is better.

Here is my function again, somewhat improved:

##rref   Reduced row echelon form
##      R = rref (A, tol) returns the reduced row echelon form of a.
##      tol defaults to eps * max (size (A)) * norm (A, inf)
##
##      [R, k] = rref (...) returns the vector of "bound variables",
##      which are those columns on which elimination has been performed.

## This code is in the public domain

## Adapted-by: Paul Kienzle

function [A, k] = rref (A, tolerance)

  ## Supress empty list warnings
  eleo = empty_list_elements_ok;
  unwind_protect
    empty_list_elements_ok = 1;
    
    [rows,cols] = size (A);

    if (nargin < 2)
      tolerance = eps * max (rows, cols) * norm (A, inf);
    endif

    used = zeros(1,cols);
    r = 1;
    for c=1:cols
      ## Find the pivot row
      [m, pivot] = max (abs (A (r:rows, c)));
      pivot = r + pivot - 1;
      
      if (m <= tolerance)
        ## Skip column c, making sure the approximately zero terms are
        ## actually zero.
        A (r:rows, c) = zeros (rows-r+1, 1);
      else
        used (1, c) = 1;
        ## Swap row p and row r
        A ([pivot, r], c:cols) = A ([r, pivot], c:cols);
        ## Divide by leading term
        A (r, c:cols) = A (r, c:cols) / A (r, c);
        ## Eliminate column c
        ridx = [1:r-1, r+1:rows];
        A (ridx, c:cols) = A (ridx, c:cols) - A (ridx, c) * A (r, c:cols);
        ## Check if done
        if (r++ == rows) break; endif
      endif
    endfor
    k = find(used);

  unwind_protect_cleanup
    ## Restore state
    empty_list_elements_ok = eleo;
  end_unwind_protect
endfunction


On Mon, Oct 16, 2000 at 08:55:36AM +0200, Dirk Laurie wrote:
> ben partan skryf:
> > Is there a built-in command for get a reduced row echelon form of a
> > matrix? Or is there an easy way to get one?
> 
> function [R,k]=rref(A,tol)
> % reduced row echelon form
> [m,n]=size(A); p=max(m,n);
> if nargin<2, tol=p*eps*norm(A,inf); end
>   % pad matrix with zeros to make it square
> if mŠþ=n, A(p,p)=0; end
> [L,U,P]=lu(A); r=sum(max(abs(U'))>tol);
> for i=1:r, k=find(U(i,:)); j(i)=k(1); end
> R=U(1:r,j)ŠÜU(1:r,:);
>   % remove added columns if any, add correct number of zero rows
> if n<p, R(:,n+1:p)=[]; end
> if length(k)<m, R(m,n)=0; end
> endfunction
> 
> Consider the above code as alpha tested only.
> 
> Dirk
> 
> 
> 
> -----------------------------------------------------------------------
> Octave is freely available under the terms of the GNU GPL.
> 
> Octave's home on the web:  http://www.che.wisc.edu/octave/octave.html
> How to fund new projects:  http://www.che.wisc.edu/octave/funding.html
> Subscription information:  http://www.che.wisc.edu/octave/archive.html
> -----------------------------------------------------------------------
> 
> 



-----------------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.che.wisc.edu/octave/octave.html
How to fund new projects:  http://www.che.wisc.edu/octave/funding.html
Subscription information:  http://www.che.wisc.edu/octave/archive.html
-----------------------------------------------------------------------