#15712: Fix determinant for RR
----------------------------------+------------------------
Reporter: ppurka | Owner:
Type: defect | Status: new
Priority: major | Milestone: sage-6.1
Component: linear algebra | Resolution:
Keywords: | Merged in:
Authors: | Reviewers:
Report Upstream: N/A | Work issues:
Branch: | Commit:
Dependencies: | Stopgaps:
----------------------------------+------------------------
Description changed by ppurka:
Old description:
> Refer to [https://groups.google.com/d/topic/sage-
> devel/Oeg6dpPkfxw/discussion this sage-devel post]
>
> {{{
> A = matrix([[ 1.0, -1.50614628068, 2.26847661882,
> -3.41665762226, 5.14598617013, -7.7506079306, 11.6735493077,
> -17.5820728722],
> [ 1.0, -0.936702701875, 0.877411951699,
> -0.821874145813, 0.769851732984, -0.721122198329, 0.675477111557,
> -0.632721235449],
> [ 1.0, -0.443181140009, 0.19640952286,
> -0.0870449962496, 0.03857670067, -0.0170964661807,
> 0.00757683137208, -0.00335790876514],
> [ 1.0, 0.352786603689, 0.124458387743,
> 0.0439072519123, 0.0154898902795, 0.00546462578321, 0.00192784677049,
> 0.000680118514595],
> [ 1.0, 0.647213396311, 0.418885180364,
> 0.271108100248, 0.175464794329, 0.11356316547, 0.07349960202,
> 0.0475699270508],
> [ 1.0, 1.44318114001, 2.08277180288,
> 3.00581698486, 4.33793838286, 6.26043086067, 9.03493574645,
> 13.0390488705],
> [ 1.0, 1.93670270187, 3.75081735545,
> 7.26421810653, 14.0686308339, 27.2467553477, 52.7688646993,
> 102.197602838],
> [ 1.0, 2.50614628068, 6.28076918019,
> 15.7405263208, 39.4480614948, 98.8626125954, 247.764168855,
> 620.933250262]])
>
> B = A.change_ring(RDF)
> print "det(A) = {}, det(B) = {}".format(A.determinant(), B.determinant())
>
> det(A) = -4.19430400000000e6, det(B) = 16801.7979988
> }}}
>
> According to Peter Bruin, a possible fix is to use pari:
>
> ''Well, it should also be fixed for a RealField of higher precision. An
> easy solution for that is to use PARI, which uses a numerically more
> stable algorithm (Gaussian elimination, choosing pivots of maximal
> absolute value; I don't know about proven error bounds). Example:''
> {{{
> sage: A._pari_().matdet()
> 16801.7979988279 # same as when doing the computation over QQ
> }}}
> ''Sage's determinant() already uses PARI over Z/nZ for n less than the
> machine word size; it would be trivial to adapt it to work also over the
> reals.''
New description:
Refer to [https://groups.google.com/d/topic/sage-
devel/Oeg6dpPkfxw/discussion this sage-devel post]
{{{
A = matrix([[ 1.0, -1.50614628068, 2.26847661882,
-3.41665762226, 5.14598617013, -7.7506079306, 11.6735493077,
-17.5820728722],
[ 1.0, -0.936702701875, 0.877411951699,
-0.821874145813, 0.769851732984, -0.721122198329, 0.675477111557,
-0.632721235449],
[ 1.0, -0.443181140009, 0.19640952286,
-0.0870449962496, 0.03857670067, -0.0170964661807, 0.00757683137208,
-0.00335790876514],
[ 1.0, 0.352786603689, 0.124458387743,
0.0439072519123, 0.0154898902795, 0.00546462578321, 0.00192784677049,
0.000680118514595],
[ 1.0, 0.647213396311, 0.418885180364,
0.271108100248, 0.175464794329, 0.11356316547, 0.07349960202,
0.0475699270508],
[ 1.0, 1.44318114001, 2.08277180288,
3.00581698486, 4.33793838286, 6.26043086067, 9.03493574645,
13.0390488705],
[ 1.0, 1.93670270187, 3.75081735545,
7.26421810653, 14.0686308339, 27.2467553477, 52.7688646993,
102.197602838],
[ 1.0, 2.50614628068, 6.28076918019,
15.7405263208, 39.4480614948, 98.8626125954, 247.764168855,
620.933250262]])
B = A.change_ring(RDF)
print "det(A) = {}, det(B) = {}".format(A.determinant(), B.determinant())
det(A) = -4.19430400000000e6, det(B) = 16801.7979988
}}}
According to Peter Bruin, a possible fix is to use pari:
''Well, it should also be fixed for a RealField of higher precision. An
easy solution for that is to use PARI, which uses a numerically more
stable algorithm (Gaussian elimination, choosing pivots of maximal
absolute value; I don't know about proven error bounds). Example:''
{{{
sage: A._pari_().matdet()
16801.7979988279 # same as when doing the computation over QQ
}}}
''Sage's determinant() already uses PARI over Z/nZ for n less than the
machine word size; it would be trivial to adapt it to work also over the
reals.''
According to Dima, the fix should be:
''Is the default choice of the algorithm the right one?
One can see that ''
{{{
sage: A.determinant(algorithm="hessenberg")
16801.7979988558
}}}
''is quite good... ''
''
Sage computes det() by computing charpoly(0), too... The division-free
algorithm is obviously meant for more general setting of rings, not
fields. I don't know why it was made default here, perhaps just an
oversight.''
--
--
Ticket URL: <http://trac.sagemath.org/ticket/15712#comment:1>
Sage <http://www.sagemath.org>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at http://groups.google.com/group/sage-trac.
For more options, visit https://groups.google.com/groups/opt_out.
