On Nov 3, 2009, at 8:13 PM, Jason Grout wrote: > 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.
I'd also like to point out that we don't just want to fall back and do everything over the rationals (even though any finite decimal expansion is rational) as things get much slower due to coefficient explosion. For example sage: m = random_matrix(ZZ, 10); m [ -1 -1 -1 -3 -23 -1 -1 0 1 0] [ -1 -1 -1 -1 3 -2 3 2 -1 0] [ 0 13 -1 0 -25 -4 -2 -1 1 -1] [ 1 7 -1 -1 0 -1 3 0 1 2] [ 2 4 -2 1 3 0 -3 1 -2 0] [ -7 -1 6 -3 0 1 -1 -60 0 -1] [ 0 2 0 -3 1 -1 -1 46 -1 -46] [ -19 -380 1 -6 2 1 3 1 -2 -1] [ 0 -1 0 -2 -2 0 -2 -5 0 0] [ 1 8 -2 1 2 0 2 0 0 8] Looks tame enough, right? sage: m^-1 [ 640562611/3385102976 -63055861/846275744 -549008479/3385102976 1157497801/1692551488 780319975/1692551488 -140389705/3385102976 -208906825/1692551488 -2793317/423137872 -1842638119/3385102976 -1535072739/1692551488] [ -5600347/1692551488 3587565/423137872 4483047/1692551488 -25271681/846275744 -16141679/846275744 2475841/1692551488 3754913/846275744 -559627/211568936 36871727/1692551488 28063787/846275744] [ 523145421/3385102976 160100405/846275744 -550981409/3385102976 306090295/1692551488 164983577/1692551488 -267861687/3385102976 -279353655/1692551488 -9830251/423137872 -1509145049/3385102976 -1738898909/1692551488] [ -360747883/3385102976 -118724611/846275744 405509879/3385102976 -128430449/1692551488 29384129/1692551488 116143857/3385102976 44713841/1692551488 3489709/423137872 -824715777/3385102976 323560411/1692551488] [ -140903663/1692551488 324025/423137872 46587723/1692551488 -36980957/846275744 -50688435/846275744 -1387523/1692551488 13708285/846275744 360201/211568936 240384563/1692551488 91072991/846275744] [ 447725143/1692551488 -73220577/423137872 -377118163/1692551488 -11011595/846275744 81470715/846275744 62756171/1692551488 15510539/846275744 -964865/211568936 -478364923/1692551488 71808441/846275744] [ 649005011/3385102976 51552491/846275744 -463894527/3385102976 487440489/1692551488 278630215/1692551488 59739863/3385102976 -67119209/1692551488 -1656309/423137872 -1667806919/3385102976 -533883395/1692551488] [ -169891/1692551488 12587301/423137872 -7854769/1692551488 -51955289/846275744 -38099159/846275744 -35116903/1692551488 -1753223/846275744 -398835/211568936 56466071/1692551488 22643/846275744] [ -185218317/423137872 -41492209/105784468 131301281/423137872 -25994167/211568936 -111115473/211568936 8319751/423137872 30984231/211568936 1213689/52892234 335120761/423137872 193991029/211568936] [ 15210787/3385102976 42522395/846275744 -36094351/3385102976 -104937127/1692551488 -71997097/1692551488 -83127577/3385102976 -46898457/1692551488 -1191381/423137872 179148969/3385102976 -39909971/1692551488] This is the typical behavior. Note also that 0.3 and 3/10 have different behaviors in Sage, sometimes you want one, sometimes you want the other. For example sage: 0.3^1000 # fast, approximate to a fixed number of digits 1.32207081948076e-523 sage: (3/10)^1000 # exact and memory intensive 132207081948080663689045525975214436596542203275214816766492036822682859 734670489954077831385060806196390977769687258235595095458210061891186534 272525795367402762022519832080387801477422896484127439040011758861804112 894781562309443806156617305408667449050617812548034440554705439703889581 746536825491613622083026856377858229022841639830788789691855640408489893 760937324217184635993869551676501894058810906042608967143886410281435038 5648747165832010614366132173102768902855220001/1000000000000000000000000 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 Note that 0.3 can't even be represented exactly as a (binary) floating point number. - Robert --~--~---------~--~----~------------~-------~--~----~ To post to this group, send email to sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support-unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-support URL: http://www.sagemath.org -~----------~----~----~----~------~----~------~--~---
