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
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
