Polynomials with uncertain coefficients

Download: Jupyter notebook, Julia script, Maple script (PACE called from Maple through UDX, see Using PACE.jl in Maple). To install PACE and run them, see Running the worksheets.

The coefficients of a polynomial often come from measurements: they are known up to an error, that is, they lie in intervals. Such a polynomial is a family of polynomials, and its "roots" are the set of the roots of all its members. We show how to enclose this set, with a proof, by combining an exact computation and interval arithmetic:

  1. solve exactly an instance of the family, with small rational coefficients;
  2. prove with the Kantorovich test that each polynomial of the family has one root, and only one, in a ball around the exact root;
  3. shrink this ball by interval Newton to an enclosure of the roots of the family.

The example is $x^3 - 2x - 5$, the polynomial on which Newton presented his method, with an uncertain constant term $-5 \pm e$.

using PACE, AbstractAlgebra
using PACE.DynamicPolynomials: @polyvar

1. An exact instance

The instance with $e = 0$ has rational coefficients: its real roots are isolated and refined exactly (see the worksheet Real roots of a univariate polynomial).

R, x = polynomial_ring(QQ, "x")
coeffs(p) = Rational{BigInt}.(collect(coefficients(p)))

function real_roots(p, prec)
    iso = isolation(coeffs(p))
    PACE.refine!(iso, prec)
    isolate(iso)
end

f0 = x^3 - 2x - 5
real_roots(f0, 64)
1-element Vector{LACE.Root}:
 Root(2.0945514815423265 .. 2.0945514815423265)

One real root. Its approximation x0 is a floating-point number, at 64 bits:

x0 = BigFloat(PACE.lo(real_roots(f0, 64)[1]); precision = 128)
2.094551481542326591481353393123243655793

2. The family

A polynomial with interval coefficients is written with DynamicPolynomials, the coefficients being intervals of MPFI (BigInterval). With $e = 1/1000$:

@polyvar y
family(e) = y^3 - 2y + (BigInterval(-5; precision = 128) +
                        BigInterval(e; precision = 128) * BigInterval(-1, 1))
F = family(1//1000)

\[ ([-5.001000000000000000000000000000000000000003030571373213709981481898885026102138, -4.998999999999999999999999999999999999999996969428626786290018518101114973897862]) + ([-2.0, -2.0])y + y^{3} \]

3. The Kantorovich test

kantorovich_test computes, in interval arithmetic, two numbers at x0: a radius $\rho$ and a number $\mu$. When $\mu < 1$, Newton's method converges from x0 to a root, unique in the ball of radius $\rho$ around x0. The computation being done on the intervals of the coefficients, the bounds hold for every polynomial of the family.

k = kantorovich_test([F], [x0], 128, 64)
2-element Vector{BigInterval}:
 [0.0, 0.000179188384956238985318]
 [0.0, 0.000201775836515190224103]

Here $\mu \approx 2 \cdot 10^{-4} < 1$: each polynomial of the family has exactly one root in $[x_0 - \rho, x_0 + \rho]$.

rho = right(k[1])
X0 = BigInterval(x0 - rho, x0 + rho)
[2.094372293157370352496035292088152370004738855868708924390375614166259765625, 2.094730669927282830466671494158334941582122468162197037599980831146240234375]

4. Interval Newton

interval_newton starts from this box. Each step intersects the Newton image with the current box, so that the roots in the box are never lost:

status, sol, _ = interval_newton([F], [X0], 128, 128, 64, 30; verbose = Cint(0))
status, sol[1]
(2, [2.09446187831001952805, 2.09464108477463365507])

The status 2 says that the box no longer shrinks: with uncertain coefficients the roots fill an interval of positive width, and the box cannot be smaller than it. As the starting box was proved to hold the roots, the result is an enclosure of them.

Here the exact set of roots is known: the root increases when the constant term decreases, so the roots of the family go from the root of $x^3 - 2x - 5 + e$ to that of $x^3 - 2x - 5 - e$. Both are isolated exactly:

e = 1//1000
a = real_roots(f0 + e, 80)[1]
b = real_roots(f0 - e, 80)[1]
Float64(PACE.lo(a)), Float64(PACE.hi(b))
(2.0944618828303576, 2.0946410712160968)

and the enclosure holds them, with little overestimation:

left(sol[1]) <= PACE.lo(a) && PACE.hi(b) <= right(sol[1])
true

5. Larger uncertainties

The same three steps for several widths:

function enclose(F, x0)
    k = kantorovich_test([F], [x0], 128, 64)
    mu, rho = right(k[2]), right(k[1])
    mu < 1 || return (Float64(mu), "not certified")
    _, sol, _ = interval_newton([F], [BigInterval(x0 - rho, x0 + rho)], 128, 128, 64, 30;
                                verbose = Cint(0))
    return (Float64(mu), sol[1])
end

[(e, enclose(family(e), x0)...) for e in (1//10^6, 1//10^3, 1//10, 1//2, 1)]
5-element Vector{Tuple{Real, Float64, BigInterval}}:
 (1//1000000, 2.0175859337701702e-7, [2.09455139194812507502, 2.09455157113652810788])
 (1//1000, 0.00020177583651519022, [2.09446187831001952805, 2.09464108477463365507])
 (1//10, 0.02034846159740268, [2.0855000129189569468, 2.1036029501656962361])
 (1//2, 0.1051943877018348, [2.04726516743844730653, 2.14183779564620587659])
 (1, 0.2190189746907231, [1.9938035002471511545, 2.19529946283750202851])

The certificate holds up to $e = 1$: all the polynomials of the family still have a single real root, and the enclosure $[1.994, 2.196]$ holds the exact set of these roots, $[2, 2.18]$. The overestimation grows with the width of the coefficients, as always with interval arithmetic.

6. When the uncertainty is too large

\[y^2 - 2y + c\]

with $c = 1 - 10^{-4}$ has two real roots, $1 \pm 1/100$. With an uncertainty $e$ on $c$, the family contains $(y-1)^2$ as soon as $e \ge 10^{-4}$: two roots merge, and beyond, some members have no real root.

c = 1 - 1//10^4
close_family(e) = y^2 - 2y + (BigInterval(c; precision = 128) +
                              BigInterval(e; precision = 128) * BigInterval(-1, 1))
starts = [BigFloat(99//100; precision = 128), BigFloat(101//100; precision = 128)]
[(e, [enclose(close_family(e), s) for s in starts]...) for e in (1//10^6, 5//10^5, 1//10^4, 2//10^4)]
4-element Vector{Tuple{Rational{Int64}, Tuple{Float64, Any}, Tuple{Float64, Any}}}:
 (1//1000000, (0.01, [0.989949747468305832657, 0.99005025253169416736]), (0.01, [1.00994974746830583259, 1.0100502525316941674]))
 (1//20000, (0.5, [0.984999999999999999959, 0.995000000000000000004]), (0.5, [1.005, 1.0150000000000000001]))
 (1//10000, (1.0, "not certified"), (1.0, "not certified"))
 (1//5000, (2.0, "not certified"), (2.0, "not certified"))

Up to $e = 5 \cdot 10^{-5}$ both roots are enclosed, in disjoint intervals: every polynomial of the family has two simple real roots, one in each. At $e = 10^{-4}$ the test fails ($\mu = 1$), as it must: no ball can hold one root of each member, since one of them has a double root. The failure of the certificate is the information: the uncertainty on the data is too large to tell the two roots apart.