#10803 closed defect (fixed)
critical bug in real_roots
Reported by: | zimmerma | Owned by: | jason, jkantor |
---|---|---|---|
Priority: | critical | Milestone: | sage-5.4.1 |
Component: | numerical | Keywords: | real roots, sd40.5 |
Cc: | cwitty, was | Merged in: | sage-5.4.1.rc0 |
Authors: | Paul Zimmermann | Reviewers: | Jeroen Demeyer |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description (last modified by )
This was reported on http://groups.google.com/group/sage-support/msg/6a7c93ad16b892bf:
sage: from sage.rings.polynomial.real_roots import real_roots sage: x = polygen(QQ) sage: f = 2503841067*x^13 - 15465014877*x^12 + 37514382885*x^11 - 44333754994*x^10 + 24138665092*x^9 - 2059014842*x^8 - 3197810701*x^7 + 803983752*x^6 + 123767204*x^5 - 26596986*x^4 - 2327140*x^3 + 75923*x^2 + 7174*x + 102 sage: len(real_roots(f)), len(real_roots(f,strategy='warp')) (11, 13)
It seems different seeds give different numbers of roots:
sage: len(real_roots(f,seed=1)) 11 sage: len(real_roots(f,seed=3)) 13
Attachments (1)
Change History (29)
comment:1 follow-up: ↓ 2 Changed 11 years ago by
comment:2 in reply to: ↑ 1 Changed 11 years ago by
Replying to dsm:
Note that it's not just that real_roots hasn't fully resolved the roots, and so an interval contains more than one root with the right multiplicity being reported, which would also lower len(real_roots(f)). In the above case two roots are simply missing, as can be confirmed by looking at the reported roots themselves.
good point. I should have written:
sage: sum(x[1] for x in real_roots(f)) 11 sage: sum(x[1] for x in real_roots(f,seed=2)) 13
Paul
comment:3 Changed 11 years ago by
- Description modified (diff)
comment:4 Changed 11 years ago by
Let me point out that f.real_roots()
does seem to work (and uses a different algorithm).
comment:5 Changed 11 years ago by
Here's what might be a clue: somehow calling the roots method on the "ocean" class fixes the problem (I don't yet have any idea why):
from sage.rings.polynomial.real_roots import * x = polygen(QQ) f = 2503841067*x^13 - 15465014877*x^12 + 37514382885*x^11 - 44333754994*x^10 + 24138665092*x^9 - 2059014842*x^8 - 3197810701*x^7 + 803983752*x^6 + 123767204*x^5 - 26596986*x^4 - 2327140*x^3 + 75923*x^2 + 7174*x + 102
ctx = context(False, 1, 32) b, other = to_bernstein(f, -1, 7) oc1 = ocean(ctx, bernstein_polynomial_factory_ratlist(b), linear_map(-1, 7)) oc1.find_roots() len(oc1.roots())
outputs 11.
oc = ocean(ctx, bernstein_polynomial_factory_ratlist(b), linear_map(-1, 7)) _ = oc.roots() oc.find_roots() print len(oc.roots())
outputs 13.
comment:6 Changed 11 years ago by
Sorry, that was misleading. Its not calling roots(), its the re-use of the context object that seems to be fixing things.
comment:7 Changed 11 years ago by
Presumably the re-use of the context object means the seed is not the same, so I guess that's a red herring - my apologies.
comment:8 follow-up: ↓ 9 Changed 11 years ago by
FYI, here's the smallest case I know:
sage: f = 535944*x^12 - 3312356*x^11 + 8042902*x^10 - 9521317*x^9 + 5204589*x^8 - 461471*x^7 - 682572*x^6 + 174775*x^5 + 25797*x^4 - 5793*x^3 - 473*x^2 + 17*x + 1 sage: sum(r[1] for r in real_roots(f,seed=1)) 10 sage: len(f.real_roots()) 12
It seems pretty sensitive to the initial bounds: (-1,7) fails, as does (-1,7+1/2^k)
for k=13,14,15,17,21,.. etc.
Unfortunately debugging is incredibly slow as real_roots.pyx takes forever to compile for me. :-(
comment:9 in reply to: ↑ 8 Changed 11 years ago by
Replying to dsm:
FYI, here's the smallest case I know:
sage: f = 535944*x^12 - 3312356*x^11 + 8042902*x^10 - 9521317*x^9 + 5204589*x^8 - 461471*x^7 - 682572*x^6 + 174775*x^5 + 25797*x^4 - 5793*x^3 - 473*x^2 + 17*x + 1 sage: sum(r[1] for r in real_roots(f,seed=1)) 10 sage: len(f.real_roots()) 12
FYI,
Mathematica 7.0 for Sun Solaris x86 (64-bit) Copyright 1988-2009 Wolfram Research, Inc. In[1]:= f = 535944*x^12 - 3312356*x^11 + 8042902*x^10 - 9521317*x^9 + 5204589*x^8 - 461471*x^7 - 682572*x^6 + 174775*x^5 + 25797*x^4 - 5793*x^3 - 473*x^2 + 17*x + 1 2 3 4 5 6 Out[1]= 1 + 17 x - 473 x - 5793 x + 25797 x + 174775 x - 682572 x - 7 8 9 10 11 > 461471 x + 5204589 x - 9521317 x + 8042902 x - 3312356 x + 12 > 535944 x In[2]:= NSolve[f==0,x] Out[2]= {{x -> -0.279345}, {x -> -0.169895}, {x -> -0.0836443}, > {x -> -0.0391184}, {x -> 0.0525278}, {x -> 0.217266}, {x -> 0.608742}, > {x -> 0.743349}, {x -> 0.955059}, {x -> 0.956605}, {x -> 1.40057}, > {x -> 1.8183}} In[3]:= Length[%] Out[3]= 12
We should add a doctest for this, to check the bug does not come back once it gets fixed.
Dave
comment:10 Changed 11 years ago by
- Milestone changed from sage-4.7 to sage-4.7.1
comment:11 Changed 11 years ago by
- Milestone changed from sage-4.7.2 to sage-4.7.3
- Priority changed from blocker to critical
comment:12 follow-up: ↓ 13 Changed 11 years ago by
- Cc was added
I'm concerned by the fact that this ticket was delayed twice (first from 4.7 to 4.7.1, now from 4.7.2 to 4.7.3, btw why is there no trace of the change from 4.7.1 to 4.7.2) and moreover the priority is now changed from blocker to critical, without any justification.
Is Carl Witty (who wrote the real_roots code as far as I know) still out there?
Paul
comment:13 in reply to: ↑ 12 Changed 11 years ago by
Replying to zimmerma:
I'm concerned by the fact that this ticket was delayed twice (first from 4.7 to 4.7.1, now from 4.7.2 to 4.7.3, btw why is there no trace of the change from 4.7.1 to 4.7.2) and moreover the priority is now changed from blocker to critical, without any justification.
Well, it seems nobody really cares much to solve this problem. I have glanced at the code and it's also hard to understand what is going on.
Personally I don't think this should a "blocker". Simple mathematical bugs are almost never blockers. This real_roots()
code is also not the default real-root computing code. If you simply do pol.roots(ring=RR)
, this will use PARI and not this buggy code (there could still be bugs in PARI though).
Is Carl Witty (who wrote the real_roots code as far as I know) still out there?
AFAIK, he is not really into Sage development anymore.
comment:14 Changed 11 years ago by
Well, it seems nobody really cares much to solve this problem.
unless someone looks at it first, I will work on this ticket during the SageFlint? days in December.
Paul
comment:16 Changed 10 years ago by
- Milestone set to sage-5.1
Is it likely that anyone's going to look at this? If not, we should either put in a stopgap (now that we have them) or take a more drastic measure.
comment:17 Changed 10 years ago by
with the example in the description, over all seeds in [0,999], the number of roots is incorrect only for 1 and 304. Unfortunately 1 is the default seed!
Paul
comment:18 Changed 10 years ago by
- Keywords sd40.5 added
- Status changed from new to needs_review
I found out that reducing the Bernstein coefficients to integers fixes the two instances of this bug (the one in the description and the one from comment 8). Maybe it does not solve all instances, but it is better than nothing, and it should not hurt to work with integers instead of rationals.
I suggest reviewers try to find other instances, for example with:
from sage.rings.polynomial.real_roots import real_roots R.<x> = ZZ[] B=10^11 for i in range(10^6): f = R.random_element((i%13)+1,-B,B) print f sys.stdout.flush() n = len(real_roots(f,seed=-1)) for s in range(1000): if len(real_roots(f,seed=s))<>n: print f, n, s sys.stdout.flush() raise ValueError
Paul
comment:19 Changed 10 years ago by
Given how sensitive the problem is I think we're going to need to know why whatever we find fixes the problem before we're confident that we've fixed it. Here I'm not sure that it does. Digging up some of my old notes, even with the patch, we have problems we shouldn't:
sage: from sage.rings.polynomial.real_roots import real_roots sage: R.<x> = ZZ[] sage: f = 535944*x^12 - 3312356*x^11 + 8042902*x^10 - 9521317*x^9 + 5204589*x^8 - 461471*x^7 - 682572*x^6 + 174775*x^5 + 25797*x^4 - 5793*x^3 - 473*x^2 + 17*x + 1 sage: len(f.roots(RR)) 12 sage: len(real_roots(f, seed=1)) 12 sage: # but.. sage: f.roots(RR) [(-0.279345282327294, 1), (-0.169895157384578, 1), (-0.0836442589128367, 1), (-0.0391184461806829, 1), (0.0525277889691767, 1), (0.217266273158562, 1), (0.608742206578326, 1), (0.743349076348040, 1), (0.955058967243006, 1), (0.956604779629507, 1), (1.40057197815310, 1), (1.81829644637645, 1)] sage: [len(real_roots(f, seed=1, bounds=(-1, 7+1/2**k))) for k in [0..200]] [12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 10, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12]
And we see a clearly periodic pattern. For me real_roots.pyx takes a long time to compile so I got frustrated chasing this clue, but changing the bounds like this -- outside the range of all of the roots! -- should have no effect. If we could figure out why this happens I think we'd have it.
comment:20 Changed 10 years ago by
I've worked a little more on this. It seems the _count_variations
function is called with a wrong usign
(sign in x=1 if I understand correctly) in the following case:
enter _count_variations, self= <IBP: ((67412088, 41929949, 16504578, -2435470, -4640) + [0 .. 18683)) * 2^-10 over [3921/16384 .. 8017/32768]; usign 1; level 3; slope_err [-41131.637071203295 .. 41131.637071203295]>
and then returns 2 changes of sign whereas there is evidently only one in the coefficients.
Paul
comment:21 Changed 10 years ago by
I believe I've found the reason for the bug: the function down_degree
decreases the
degree of the (interval) polynomial, but keeps the value self.usign
. With the example in the description and seed=1
, down_degree
is called with:
enter down_degree, self= <IBP: ((34519150034, 25826190301, 17130808432, 9120071648, 2724478006, -830629355, 3694263) + [0 .. 21299)) * 2^-19 over [3921/16384 .. 8017/32768]; level 2; slope_err [-136.91992834614365 .. 136.91992834614365]> usign= 1
and the polynomial is reduced to:
coeffs= (67412088, 41929949, 16504578, -2435470, -4640) n= 5
where it is clear that usign=1
is wrong, since the sign of -4640
is -1.
Paul
Changed 10 years ago by
comment:22 Changed 10 years ago by
the new patch fixes the bug. Ready for review. Thanks Douglas for pushing me to really isolate the reason of the bug.
Paul
comment:23 Changed 10 years ago by
I just checked: the example from comment 19 passes too with the new patch.
Paul
comment:24 Changed 10 years ago by
Looks like nice work! I think we're going to have to add way more testing, though (I can do some of it as a reviewer patch this weekend.) I have a collection of failing cases and I want to make sure that more of them work; we should specifically doctest your diagnosis, not just the outcome; and I want to test that the bound stuff now works.
Now that we know the problem it should be much easier to construct more test cases, too.
comment:25 Changed 10 years ago by
- Milestone changed from sage-5.4 to sage-5.5
- Reviewers set to Jeroen Demeyer
- Status changed from needs_review to positive_review
Since the reviewer patch didn't happen, let's put positive_review here. You can still add more examples on a new ticket.
comment:26 Changed 10 years ago by
- Merged in set to sage-5.5.beta0
- Resolution set to fixed
- Status changed from positive_review to closed
comment:27 Changed 10 years ago by
- Description modified (diff)
- Merged in changed from sage-5.5.beta0 to sage-5.4.1.rc0
comment:28 Changed 10 years ago by
- Milestone changed from sage-5.5 to sage-5.4.1
Note that it's not just that real_roots hasn't fully resolved the roots, and so an interval contains more than one root with the right multiplicity being reported, which would also lower len(real_roots(f)). In the above case two roots are simply missing, as can be confirmed by looking at the reported roots themselves.