Michael Orlitzky wrote: > Jason Grout wrote: > >> The problem here is that no one has really implemented a numerically >> stable version of echelon_form for RR. I thought we called scipy for >> rank calculations over RDF, but apparently we don't even do that! Scipy >> answers correctly: >> >> sage: import scipy >> sage: scipy.rank(m.numpy()) >> 2 >> >> >> Very little attention has been paid to numerically stable linear algebra >> for non-exact rings. We do get some numerically stable things for >> matrices over RDF and CDF because those are based on numpy and scipy >> (which call blas and lapack routines). However, apparently we don't >> call numpy/scipy to calculate the rank, but instead rely on our unstable >> echelon form computation! > > This was hinted at by the (m == n) example in my original message, but > if I can type it, it isn't an irrational number. Since any (decimal) > number entered manually is guaranteed to be rational, it can be > represented as a fraction, and those are handled correctly. > > This is probably already implemented to some extent, or else (m == n) > would not have returned True. > > >> If someone wanted to take make a good library of numerically stable >> linear algebra routines based on mpfr, I think it would be absolutely >> awesome. > > Irrationals have to be entered symbolically, and SAGE can already > manipulate symbols well-enough. I would like to say, "just compute the > rref symbolically," but that idea breaks down somewhat given two > irrationals that are linearly dependent. > > I should review my algebra before saying something dumb, but I think > that arbitrary precision would be useful in that case. >
I don't think it's an issue of irrational versus rational. It's numerical precision and inexact floating point numbers. This matrix is terribly ill-conditioned. It is right on the border line between being invertible or not, numerically speaking: sage: m.change_ring(RDF).SVD()[1] [ 0.772642968023 0.0 0.0] [ 0.0 0.450580563234 0.0] [ 0.0 0.0 3.13289758759e-17] As you probably know, the ratio between the smallest and largest eigenvalues being so high gives us an indication that this matrix is really a messy one numerically, and deserves the strictest of attention. If we are allowed to do the computations with exact arithmetic (e.g., using fractions instead of decimals), then we can compute its rank and echelon form exactly. You see the same problems in other math software too, with different matrices. If pressed, I could come up with an a very nice-looking matrix example from a linear algebra textbook problem that matlab totally gave the wrong answer to because it was a very ill-conditioned matrix. Thanks, Jason -- Jason Grout --~--~---------~--~----~------------~-------~--~----~ To post to this group, send email to [email protected] To unsubscribe from this group, send email to [email protected] For more options, visit this group at http://groups.google.com/group/sage-support URL: http://www.sagemath.org -~----------~----~----~----~------~----~------~--~---
