I have a collection of polynomials with non-negative integer coefficients that I know to be real-rooted with no repeated roots by examining Sturm chains for each. I would like to show that they form an interlacing family (hence their convex combinations are real-rooted). Coming from Maple, I would use the realroot command (based on Descartes’ rule of signs) to find rationally bounded isolating intervals for the real roots of the polynomials. Is there a package that implements something equivalent in Julia? For a ballpark idea of scale, I’d like to be able to work with 10^5 polynomials with degree 80, and coefficients on the order of 2^480, but something working out-of-the-box on smaller instances would still be helpful.
You will likely find IntervalRootFinding to be helpful here. It will give rigorous bounds. (The issue with roots from Polynomials is the high degree will likely yield numeric issues). Here is a sample usage, though it sounds like you might use the Polynomial(coeffs) constructor:
julia> using Polynomials, IntervalRootFinding, IntervalArithmetic
julia> x = [2,3,4,5,7,8,9];
julia> p = fromroots(Polynomial, x)
Polynomial(-60480 + 100536*x - 68078*x^2 + 24433*x^3 - 5036*x^4 + 598*x^5 - 38*x^6 + x^7)
julia> IntervalRootFinding.roots(x->p(x), -20..20)
7-element Array{Root{Interval{Float64}},1}:
Root([2.99999, 3.00001], :unique)
Root([1.99999, 2.00001], :unique)
Root([7.99999, 8.00001], :unique)
Root([4.99999, 5.00001], :unique)
Root([6.99999, 7.00001], :unique)
Root([8.99999, 9.00001], :unique)
Root([3.99999, 4.00001], :unique)
The :unique shows there is a root and only 1 in the specified interval.
The choice of 20..20 can be widened without issue, and could be computed using some bounds on real roots based on coefficients.
To show the nesting, you might do something like this:
xs = IntervalRootFinding.roots(x->p(x), -20..20)
ys = IntervalRootFinding.roots(x->Polynomials.derivative(p)(x), -20..20)
xs = sort([i.interval for i in xs])# grab intervals only
ys = sort([i.interval for i in ys]);
all(xs[i] < ys[i] < xs[i+1] for i in 1:length(xs)-1) # true
3 Likes
IntervalRootFinding.roots doesn’t seem to be able to take advantage of the fact that the coefficients are integers. I’m finding it chokes on relatively low degree polynomials.
using IntervalRootFinding, IntervalArithmetic, Polynomials
bp = [
942438915208811912419937422298363203125
164182217245953398816894035758761902846875
4584900574568933770264468813466870772155175
48995332393110515074735075708247882042540865
266674183150777010544241114017741621823207005
852443280934837985352088128423887894438557515
1738546146302892245736990844587487000484756535
2381558158813900978436742173305983349418813145
2262003889258248241081177038309445610985409335
1516025051068561122302699855213604175083575145
720810987764796866354279087485114863858085005
241341213302325116160821160849326697595681275
55563538585119205063994483187179167406616375
8363912237256094118085070946083688310200625
740493466743082745510080711751444519503125
29215606371473169285018060091249259296875]
intervalroots(a) = IntervalRootFinding.roots(x->a(x), -big(4)..big(0))
poly = Polynomials.Polynomial(bp)
@btime intervalroots(poly)
this takes minutes at 100% CPU usage, though relative low memory. In contrast, in Maple, I can run
cl := [ 942438915208811912419937422298363203125,
164182217245953398816894035758761902846875,
4584900574568933770264468813466870772155175,
48995332393110515074735075708247882042540865,
266674183150777010544241114017741621823207005,
852443280934837985352088128423887894438557515,
1738546146302892245736990844587487000484756535,
2381558158813900978436742173305983349418813145,
2262003889258248241081177038309445610985409335,
1516025051068561122302699855213604175083575145,
720810987764796866354279087485114863858085005,
241341213302325116160821160849326697595681275,
55563538585119205063994483187179167406616375,
8363912237256094118085070946083688310200625,
740493466743082745510080711751444519503125,
29215606371473169285018060091249259296875];
poly := add(cl[i]*x^(i-1), i=1..nops(cl));
realroot(poly,1/2^200);
to get interval that isolate to intervals of width 2^(-200) in on the order of 1/50th of a second.
I have a naive, hacked together, Sturm chain approach in Julia that gets me much wider intervals in a couple of seconds, but I can’t see how I’m going to make up that many orders of magnitude of difference.
I’m way out of my depth here, but have you looked at Real-root isolation - Wikipedia?
1 Like
That looks like about what I need. The Recent Improvements section there, sure makes it look like there should be some interval arithmetic trick close by. I guess tomorrow’s project is to try and limit the number of intermediate values for all those change of variables clean when working with BigInts.
I don’t see how the fact that the coefficients are integers will help? But it’s very interesting that Maple is so fast.
You might also be interested in HomotopyContinuation.jl
1 Like
It seems that you can get some initial information quickly using Float64:
julia> bp2 = Interval{Float64}.(bp);
julia> poly2 = Polynomials.Polynomial(bp2);
julia> roots(x -> poly2(x), -Inf..Inf, Bisection, 1.0)
22-element Vector{Root{Interval{Float64}}}:
Root([-6.95784, -6.40872], :unknown)
Root([-9.129, -8.57124], :unknown)
[etc]
Then hopefully refining with BigFloat will be faster.
For fast polynomial manipulations in Julia I believe you can use https://github.com/Nemocas/AbstractAlgebra.jl
If you have an approximate root, you can verify it, e.g.:
julia> using Roots
julia> rts = find_zeros(x -> poly(x), -10, 10); # approximate: Roots.jl
julia> big_roots = roots(x -> poly(x), big(rts[3] ± 1e-5), Newton, 1e-50) # rigorous: IntervalRootFinding.jl
1-element Vector{Root{Interval{BigFloat}}}:
Root([-3.39911, -3.3991]₂₅₆, :unique)
julia> diam(interval(big_roots[1]))
7.879417721966192043308957487195098830393286784342275257087148835580431618280549e-70
This finds an interval of with ~10^(-69) in which a unique root is proved to exist.
You can also use Hecke/Nemo. Here are the timings on my very dated laptop:
julia> using Hecke
julia> Qx, x = PolynomialRing(QQ)
(Univariate Polynomial Ring in x over Rational Field, x)
julia> Qx, x = PolynomialRing(QQ);
julia> f = Qx(bp);
julia> Hecke._roots(f, 200) # the balls will have radious < 2^(-200)
15-element Array{acb,1}:
[-3.920616711153724253790204590529063281656171296624620351746200485310655120855543613013541271680907275607223979079 +/- 4.84e-112] + i*0
[-3.706111517064614641012473797293410340554769655675971950139715152216337896330719648922654909467811143497620428030 +/- 7.80e-112] + i*0
[-3.399101345956385184809590396663823245233910492969600172371668237761238821534096232299586114103416940774067560924 +/- 8.79e-112] + i*0
[-3.035412954945010421040655247164632911783923123071964813560272462311419306648072410827429032399236488318520657029 +/- 9.05e-112] + i*0
[-2.642241820837020072256321446404090852721123523348319017917035488525807038562259125380960186537327785108752839951 +/- 5.03e-112] + i*0
[-2.240386668347612606849921371192140625120311772030273906303162991172574142894086307045961801604703579352289007669 +/- 3.04e-112] + i*0
[-1.8459392275279762188118911328074682418106652190523122728768478452769447419809494076465189839407551707688825849549 +/- 9.17e-113] + i*0
[-1.4714364265485984347457103358546916040139602018877260092947033620076406186360124337515593788219523840023676498250 +/- 3.09e-113] + i*0
[-1.12667548938228648371659194465060333123610781848166932944102272185644329811383422404304887752307268162794786601481 +/- 6.41e-114] + i*0
[-0.819305999801997298040859319753994185201032318699276601025148883053951656453384749232848418410546346810549845454155 +/- 6.49e-115] + i*0
[-0.5552656322189970256982534594202473203239575869892228528269814900273624281956924496928014910474082527863381439566599 +/- 8.75e-116] + i*0
[-0.33910202510401947795464384494964958658273285483978651160525210708434910971533965735125598235860614612011137885597364 +/- 4.76e-117] + i*0
[-0.174209253846368310831930193848706110784111374105168275102590356950016921390554301061214369088829872621188752931599096 +/- 4.12e-118] + i*0
[-0.06299838108334910672357500092170790024222765409891867028284948240638078820689793485367480324463392976269126659634877413 +/- 9.09e-120] + i*0
[-0.007015398195279156379264993894770551409664575420748371105657755069724334773977198396431121866041220489066572460776894302 +/- 2.04e-121] + i*0
julia> @btime Hecke._roots(f, 200);
4.501 ms (3677 allocations: 204.58 KiB)
This is internally using http://arblib.org/.
P.S. To convert to BigFloat you can do BigFloat(real(x)), but you must handle the global precision yourself.
2 Likes
I’m not really clear myself, I just know from browsing the Maple documentation, that they rely on integer (or rational) coefficients, presumably because it simplifies something. Page 2 of https://arxiv.org/pdf/1605.00410.pdf seems to suggest that one of the problems is being able to guarantee that it’s possible to determine the sign of a polynomial near its roots. It also seems (page 8, paragraph 2), that there might be some implications to a heuristic that switches between MPFI and exact arithmetic depending on which is expected to be cheaper for a given subproblem.
Another guarantee from integer coefficients is that it’s possible to pre-filter and ensure no rational roots, so we don’t have to worry about corner cases where a root falls exactly at an end-point of an interval.
2 Likes
Thanks. That’s definitely fast enough to fit into my pipeline.
I wanted to see if my initial suggestion to use IntervalRootFinding could be sped up paying more attention to the polynomial evaluation time and the underlying numbers used, so tried this:
using IntervalRootFinding, IntervalArithmetic, Polynomials
bp = [
942438915208811912419937422298363203125
164182217245953398816894035758761902846875
4584900574568933770264468813466870772155175
48995332393110515074735075708247882042540865
266674183150777010544241114017741621823207005
852443280934837985352088128423887894438557515
1738546146302892245736990844587487000484756535
2381558158813900978436742173305983349418813145
2262003889258248241081177038309445610985409335
1516025051068561122302699855213604175083575145
720810987764796866354279087485114863858085005
241341213302325116160821160849326697595681275
55563538585119205063994483187179167406616375
8363912237256094118085070946083688310200625
740493466743082745510080711751444519503125
29215606371473169285018060091249259296875]
bps = [convert(ArbFloat, string(b)) for b in bp]
p = ImmutablePolynomial(bps)
# error with show
# Error showing value of type Tuple{ArbFloat{116},ArbFloat{116}}:
# ERROR: MethodError: no method matching Array{String,1}(::Int64)
a,b = -convert(ArbFloat, "4"), convert(ArbFloat, "0");
IntervalRootFinding.roots(f,a..b) # error about a..b
#ERROR: MethodError: no method matching ArbFloat{116}(::String)
IntervalRootFinding.roots(f,Interval(a,b))
# ERROR: MethodError: *(::ArbFloat{116}, ::Interval{ArbFloat{116}}) is ambiguous. Candidates:
So a bunch of errors, the last one non-circumventable. I could file issues with ArbFloats for the first two, but the method ambiguity doesn’t have a natural place to resolve.