help-octave
[Top][All Lists]

## Re: reduced row echelon form

 From: Paul Kienzle Subject: Re: reduced row echelon form Date: Mon, 16 Oct 2000 13:36:22 +0100 User-agent: Mutt/1.2.5i

```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
-----------------------------------------------------------------------

```

reply via email to