critical bug in real_roots

[zimmerma]

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

[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.

[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

[jdemeyer]

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

[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.

[mhampton]

Sorry, that was misleading. Its not calling roots(), its the re-use of the context object that seems to be fixing things.

[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.

[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. :-(

[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

[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.

Is Carl Witty (who wrote the real_roots code as far as I know) still out there?

Paul

[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.

[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

[jdemeyer]

Milestone sage-4.7.3 deleted

[dsm]

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.

[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

[zimmerma]

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

[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.

[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

[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

[zimmerma]

the new patch fixes the bug. Ready for review. Thanks Douglas for pushing me to really isolate the reason of the bug.

Paul

[zimmerma]

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

Paul

[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.

[jdemeyer]

Since the reviewer patch didn't happen, let's put positive_review here. You can still add more examples on a new ticket.