Opened 11 years ago

Closed 10 years ago

Last modified 10 years ago

#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:

Status badges

Description (last modified by jdemeyer)

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)

trac_10803.patch (1.6 KB) - added by zimmerma 10 years ago.

Download all attachments as: .zip

Change History (29)

comment:1 follow-up: Changed 11 years ago by 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.

comment:2 in reply to: ↑ 1 Changed 11 years ago by zimmerma

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 jdemeyer

  • Description modified (diff)

comment:4 Changed 11 years ago by jdemeyer

Let me point out that f.real_roots() does seem to work (and uses a different algorithm).

comment:5 Changed 11 years ago by mhampton

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 mhampton

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 mhampton

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: Changed 11 years ago by 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

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 drkirkby

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 jdemeyer

  • Milestone changed from sage-4.7 to sage-4.7.1

comment:11 Changed 11 years ago by jdemeyer

  • Milestone changed from sage-4.7.2 to sage-4.7.3
  • Priority changed from blocker to critical

comment:12 follow-up: Changed 11 years ago by zimmerma

  • 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 jdemeyer

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 zimmerma

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:15 Changed 11 years ago by jdemeyer

  • Milestone sage-4.7.3 deleted

Milestone sage-4.7.3 deleted

comment:16 Changed 10 years ago by dsm

  • 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 zimmerma

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 zimmerma

  • Authors set to Paul Zimmermann
  • 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 dsm

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 zimmerma

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 zimmerma

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 zimmerma

comment:22 Changed 10 years ago by zimmerma

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 zimmerma

I just checked: the example from comment 19 passes too with the new patch.

Paul

comment:24 Changed 10 years ago by dsm

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 jdemeyer

  • 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 jdemeyer

  • 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 jdemeyer

  • Description modified (diff)
  • Merged in changed from sage-5.5.beta0 to sage-5.4.1.rc0

comment:28 Changed 10 years ago by jdemeyer

  • Milestone changed from sage-5.5 to sage-5.4.1
Note: See TracTickets for help on using tickets.