[Numpy-discussion] Reduced row echelon form

[Anne Archibald]

2008/11/18 Robert Young <rob at roryoung.co.uk>:

> Is there a method in NumPy that reduces a matrix to it's reduced row echelon
> form? I'm brand new to both NumPy and linear algebra, and I'm not quite sure
> where to look.

Unfortunately, reduced row-echelon form doesn't really work when using
approximate values: any construction of the reduced row-echelon form
forces you many times to ask "is this number exactly zero or just very
small?"; if it's zero you do one thing, but if it's very small you do
something totally different - usually divide by it. With
floating-point numbers, every calculation is approximate, and such a
method will blow up completely. If you really need reduced row echelon
form, you have to start with exact numbers and use a package that does
exact computations (I think SymPy might be a place to start).

In practice, numerical linear algebra is rather different from linear
algebra as presented in math classes. In particular, problems that you
might solve on the chalkboard with row reduction or the like are
instead solved by matrix factorizations of special forms. For example
LU factorization writes a matrix as a product of a lower-triangular
matrix and an upper-triangular matrix. This allows, for example, very
easy calculation of determinants. (It also allows fast solution of
linear equations, just like reduced row echelon form.) But LU
factorization is much more resistant to the problems involved in
working with approximate numbers.

If you have a problem that is classically solved with something like
reduced row echelon form, you first need to think about how to make it
make sense in an approximate setting. For example, the rank of a
matrix: if two rows are exactly equal, the matrix is singular. But if
the rows are even slightly different, the matrix is non-singular.
There's just no way to make this work precisely using approximate
numbers. (Sometimes you can rephrase the problem in a way that does
work; singular value decomposition lets you deal with ranks in a
sensible fashion, by giving a reasonable criterion for when you want
to consider a linear combination equal to zero.) If, however, your
question makes sense with approximate numbers (solving Ax=b usually
does, for example) but your algorithm for getting the answer doesn't
work, look for a numerical matrix decomposition that will work. The
singular value decomposition is the swiss army knife for this sort of
problem, but others can work better in some cases.

Numerical Recipes is the traditional book to recommend in this sort of
case. Their algorithms may not be the best, and don't use their code,
but their descriptions do instill a sensible caution.

Good luck,
Anne