I've opened https://trac.sagemath.org/ticket/34724 to deal with this
On Sat, Nov 5, 2022 at 11:08 AM Dima Pasechnik <dimp...@gmail.com> wrote:
>
> Well, applying a naive exact linear algebra routine to inexact data,
> and that's what Sage is doing here, is prone to errors.
> sage: A=m.augment(identity_matrix(RR,2))
> sage: A
> [-6.12323399573677e-17 -1.72508242466029 1.00000000000000
> 0.000000000000000]
> [ 0.579682446302195 6.12323399573677e-17 0.000000000000000
> 1.00000000000000]
> sage: A.echelonize(algorithm="classical");A # OOOPS! - that's the default here
> [ 1.00000000000000 0.000000000000000 4.00000000000000
> 1.72508242466029]
> [ 0.000000000000000 1.00000000000000 -0.579682446302195
> -6.12323399573676e-17]
> sage: A=m.augment(identity_matrix(RR,2))
> sage: A.echelonize(algorithm='scaled_partial_pivoting');A # that's how
> it should be
> [ 1.00000000000000 0.000000000000000 6.12323399573677e-17
> 1.72508242466029]
> [ 0.000000000000000 1.00000000000000 -0.579682446302195
> -6.12323399573677e-17]
>
> On Sat, Nov 5, 2022 at 10:21 AM Emmanuel Charpentier
> <emanuel.charpent...@gmail.com> wrote:
> >
> > something is definitely unhinged here : On 9.8.beta3 running on Debian
> > testing on core i7 + 16 GB RAM, after running :
> >
> > a = RR(-4967757600021511 / 2**106)
> > b = RR(-7769080564883485 / 2**52)
> > c = RR( 5221315298319565 / 2**53)
> > m = matrix([[a, b], [c, -a]])
> > M = matrix([[var("p%d%d"%(u, v), latex_name="p_{%s,%d}"%(u, v))
> > for v in range(2)]
> > for u in range(2)])
> > S = dict(zip(M.list(), [a, b, c, -a]))
> > MN = M.apply_map(lambda u:u.subs(S))
> >
> > one gets :
> >
> > sage: m.parent()
> > Full MatrixSpace of 2 by 2 dense matrices over Real Field with 53 bits of
> > precision
> > sage: m*~m
> > [ 1.00000000000000 -1.23259516440783e-32]
> > [ 2.31872978520878 1.00000000000000]
> > sage: (~m)*m
> > [ 1.00000000000000 -6.90032969864117]
> > [6.16297582203915e-33 1.00000000000000]
> > sage: MN.parent()
> > Full MatrixSpace of 2 by 2 dense matrices over Symbolic Ring
> > sage: MN*~MN
> > [ 1.00000000000000 -1.23259516440783e-32]
> > [ 2.31872978520878 1.00000000000000]
> > sage: (~MN)*MN
> > [ 1.00000000000000 -6.90032969864117]
> > [6.16297582203915e-33 1.00000000000000]
> >
> > all being wrong, wrong, wrong…
> >
> > However :
> >
> > sage: (M*~M).apply_map(lambda u:u.subs(S))
> > [ 1.00000000000000 0]
> > [-3.54953126192945e-17 1.00000000000000]
> > sage: ((~M)*M).apply_map(lambda u:u.subs(S))
> > [ 1.00000000000000 1.05630833481279e-16]
> > [ 0 1.00000000000000]
> >
> > both being acceptable.
> >
> > One also notes that the form of :
> >
> > sage: ~M
> > [1/p00 - p01*p10/(p00^2*(p01*p10/p00 - p11))
> > p01/(p00*(p01*p10/p00 - p11))]
> > [ p10/(p00*(p01*p10/p00 - p11))
> > -1/(p01*p10/p00 - p11)]
> > sage: (~M).apply_map(simplify)
> > [1/p00 - p01*p10/(p00^2*(p01*p10/p00 - p11))
> > p01/(p00*(p01*p10/p00 - p11))]
> > [ p10/(p00*(p01*p10/p00 - p11))
> > -1/(p01*p10/p00 - p11)]
> >
> > is somewhat unexpected ; one expects :
> >
> > sage: (~M).apply_map(lambda u:u.simplify_full())
> > [-p11/(p01*p10 - p00*p11) p01/(p01*p10 - p00*p11)]
> > [ p10/(p01*p10 - p00*p11) -p00/(p01*p10 - p00*p11)]
> >
> > which is also the form returned by maxima :
> >
> > sage: maxima_calculus.invert(M).sage()
> > [-p11/(p01*p10 - p00*p11) p01/(p01*p10 - p00*p11)]
> > [ p10/(p01*p10 - p00*p11) -p00/(p01*p10 - p00*p11)]
> >
> > giac :
> >
> > sage: giac.inverse(giac(M)).sage()
> > [[-p11/(p01*p10 - p00*p11), p01/(p01*p10 - p00*p11)],
> > [p10/(p01*p10 - p00*p11), -p00/(p01*p10 - p00*p11)]]
> >
> > fricas :
> >
> > sage: fricas.inverse(M._fricas_()).sage()
> > [-p11/(p01*p10 - p00*p11) p01/(p01*p10 - p00*p11)]
> > [ p10/(p01*p10 - p00*p11) -p00/(p01*p10 - p00*p11)]
> >
> > mathematica :
> >
> > sage: mathematica.Inverse(M).sage()
> > [[-p11/(p01*p10 - p00*p11), p01/(p01*p10 - p00*p11)],
> > [p10/(p01*p10 - p00*p11), -p00/(p01*p10 - p00*p11)]]
> >
> > and (somewhat un-backconvertible) :
> >
> > sage: sympy.sympify(M)^-1._sage_()
> > Matrix([
> > [ p11/(p00*p11 - p01*p10), -p01/(p00*p11 - p01*p10)],
> > [-p10/(p00*p11 - p01*p10), p00/(p00*p11 - p01*p10)]])
> >
> > This is, IMNSHO, a critical bug. Could you open a tichet for this, and mark
> > it as such ?
> >
> > Le samedi 5 novembre 2022 à 07:59:27 UTC+1, Håkan Granath a écrit :
> >>
> >> Hi, there seems to be a problem with inverses of matrices with elements in
> >> RR. It only occurs very sporadically for me, but here is an example:
> >>
> >> a = RR(-4967757600021511 / 2**106)
> >> b = RR(-7769080564883485 / 2**52)
> >> c = RR( 5221315298319565 / 2**53)
> >>
> >> m = matrix([[a, b], [c, -a]])
> >>
> >> print(m)
> >> print()
> >> print(~m)
> >>
> >> On my machines it produces the output
> >>
> >> [-6.12323399573677e-17 -1.72508242466029]
> >> [ 0.579682446302195 6.12323399573677e-17]
> >>
> >> [ 4.00000000000000 1.72508242466029]
> >> [ -0.579682446302195 -6.12323399573676e-17]
> >>
> >> Clearly the element 4 is wrong (the correct inverse is -m). Is this a
> >> known bug?
> >>
> >> Some system information:
> >>
> >> SageMath version 9.7, using Python 3.10.5
> >> OS: Ubuntu 20.04.5 LTS
> >> CPU: Intel(R) Core(TM) i7-7700 CPU @ 3.60GHz
> >>
> >> Best regards,
> >>
> >> Håkan Granath
> >
> > --
> > You received this message because you are subscribed to the Google Groups
> > "sage-devel" group.
> > To unsubscribe from this group and stop receiving emails from it, send an
> > email to sage-devel+unsubscr...@googlegroups.com.
> > To view this discussion on the web visit
> > https://groups.google.com/d/msgid/sage-devel/bec2d18e-8034-440c-8076-9a5e9ec93d2cn%40googlegroups.com.
--
You received this message because you are subscribed to the Google Groups
"sage-devel" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to sage-devel+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/sage-devel/CAAWYfq0d71J6eLmCratONnrDYV2VbvbCC2c0_5Ng2nAgaA5V0w%40mail.gmail.com.
