Opened 8 years ago
Last modified 5 years ago
#15712 new defect
Fix determinant for RR
Reported by: | ppurka | Owned by: | |
---|---|---|---|
Priority: | major | Milestone: | sage-6.4 |
Component: | linear algebra | Keywords: | |
Cc: | Merged in: | ||
Authors: | Reviewers: | ||
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description (last modified by )
Refer to 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.
Change History (5)
comment:1 Changed 8 years ago by
- Description modified (diff)
comment:2 Changed 8 years ago by
- Milestone changed from sage-6.1 to sage-6.2
comment:3 Changed 8 years ago by
- Milestone changed from sage-6.2 to sage-6.3
comment:4 Changed 8 years ago by
- Milestone changed from sage-6.3 to sage-6.4
now (sage 7.6 beta 4)
also returns the faulty -4.19430400000000e6