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