2008/11/18 Robert Young <

[hidden email]>:

> 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

_______________________________________________

Numpy-discussion mailing list

[hidden email]
http://projects.scipy.org/mailman/listinfo/numpy-discussion