On Thu, Sep 8, 2011 at 7:09 PM, Dima Pasechnik <dimp...@gmail.com> wrote: > > > On Thursday, 8 September 2011 23:55:00 UTC+8, William wrote: >> >> On Thu, Sep 8, 2011 at 8:28 AM, Rado <rki...@gmail.com> wrote: >> > sage: L = [[2,4,16],[1/2,1,2],[1,1,1]] >> > sage: A = matrix(RDF, L) >> > sage: xvalues = [sqrt(2), 2, 4] >> > sage: L2 = [[i^2 for i in xvalues], [log(i, 2) for i in xvalues], [1 for >> > i >> > in xvalues]] >> > sage: A2 = matrix(RDF, L2) >> > sage: print (A - A2).norm() >> > 1.11022302463e-16 >> > sage: b = vector(RDF, [5, 20, 106]) >> > sage: print A.solve_right(b) >> > (153.25, -37.875, -9.375) >> > sage: print A.inverse() * b >> > (153.25, -37.875, -9.375) >> > sage: print A2.solve_right(b) >> > (128.0, -32.0, -9.375) >> > sage: print A2.inverse() * b >> > (153.25, -37.875, -9.375) >> > The third output is quite off and I don't think it can be attributed to >> > numerical error. Any idea what is happening? >> > Rado >> >> Numerical error. >> >> For double precision matrices, at least in sage-4.7.2, solve_right >> (and solve_left) still just calls the completely generic algorithm >> (from school), which is designed for exact rings... in particular, it >> creates the augmented matrix C: >> >> sage: C = A2.augment(b); C >> [ 2.0 4.0 16.0 5.0] >> [ 0.5 1.0 2.0 20.0] >> [ 1.0 1.0 1.0 106.0] >> >> Then, it calls echelon on that, which just runs a very simple >> "10-line" standard Gauss elimination, which is designed assuming that >> the base ring is exact. This is "def >> _echelon_in_place_classical(self):" in matrix/matrix2.pyx. > > even over exact rings, this is by far not the right thing to do, as it can > lead to explosions in the size. > one should do Bareiss, if possible. > http://en.wikipedia.org/wiki/Bareiss_algorithm > http://www.singular.uni-kl.de/Manual/latest/sing_186.htm
There are far more sophisticated algorithms than even Bareiss in Sage (in the matrix directory). For example, there is a complete general implementation of echelon form via a form of Strassen recursive decomposition, which really does do the general case of nxm matrices -- this was a painful bit of work by David Harvey and Robert Bradshaw, and is available over any base ring, and reduces the problem to matrix multiplication. There's also code I wrote on top of IML (integer matrix library) for computing echelon forms very efficiently of certain shape matrices over QQ, which uses p-adic system solving (e.g., Dickson's method) under the hood. And there's also all of linbox, which implements a range of sophisticated algorithms over certain base rings (mostly finite fields and QQ and ZZ). -- William > > Just in case. > Dima > >> >> sage: C.echelon_form() >> [ 1.0 0.0 0.0 128.0] >> [ 0.0 1.0 0.0 -32.0] >> [ 0.0 0.0 1.0 -9.375] >> >> And there's your mystery "answer" in the last column. >> >> Remember that pretty much everything is nonzero and nonequal in RDF, >> which confuses Gauss elimination. I put a print statement in the >> function mentioned above, and get this output during the >> echelonization. >> >> sage: C.echelon_form() >> Pivoting on self[0,0] = 2.0 >> Now self = >> [ 1.0 2.0 8.0 >> 2.5] >> [ 0.0 -2.22044604925e-16 -2.0 >> 18.75] >> [ 0.0 -1.0 -7.0 >> 103.5] >> Pivoting on self[1,1] = -2.22044604925e-16 >> Now self = >> [ 1.0 0.0 -1.80143985095e+16 >> 1.68884986026e+17] >> [ 0.0 1.0 9.00719925474e+15 >> -8.44424930132e+16] >> [ 0.0 0.0 9.00719925474e+15 >> -8.44424930132e+16] >> Pivoting on self[2,2] = 9.00719925474e+15 >> Now self = >> [ 1.0 0.0 0.0 128.0] >> [ 0.0 1.0 0.0 -32.0] >> [ 0.0 0.0 1.0 -9.375] >> >> Note that one source of confusion is missing digits in printing. >> E.g., if you paste the numbers back in, >> then you get the right answer: >> >> sage: C = matrix(RDF,3,4,[2.0, 4.0, 16.0, 5.0, 0.5, 1.0, 2.0, 20.0, >> 1.0, 1.0, 1.0, 106.0]) >> sage: C.echelon_form() > > -- > To post to this group, send an email to sage-devel@googlegroups.com > To unsubscribe from this group, send an email to > sage-devel+unsubscr...@googlegroups.com > For more options, visit this group at > http://groups.google.com/group/sage-devel > URL: http://www.sagemath.org > -- William Stein Professor of Mathematics University of Washington http://wstein.org -- To post to this group, send an email to sage-devel@googlegroups.com To unsubscribe from this group, send an email to sage-devel+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-devel URL: http://www.sagemath.org
