Systems 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 worksheet Polynomials with uncertain coefficients enclosed the roots of a family of univariate polynomials. We do the same for systems: the coefficients lie in intervals, and we look for boxes that hold, for every system of the family, exactly one solution. The three steps are the same:

  1. solve exactly an instance of the family, with rational coefficients (worksheet Real solutions of a zero-dimensional system);
  2. at each real solution, prove with the Kantorovich test that every system of the family has one solution, and only one, in a ball around it;
  3. shrink these balls by interval Newton.

A last section shows what this approach cannot see.

using PACE, AbstractAlgebra
using PACE.DynamicPolynomials: @polyvar

1. The exact instance

The system of the previous worksheet, the cubic $y^2 = x^3 - x$ and the circle $x^2 + y^2 = 3$:

R, (x, y) = polynomial_ring(QQ, ["x", "y"])
sys0 = [y^2 - x^3 + x, x^2 + y^2 - 3]
sols = real_solutions(sys0; prec = 64)
2-element Vector{Vector{BigInterval}}:
 [[1.359304085971776420730660392788800573278, 1.359304085971776420730660392812266330097], [-1.07344883523176522996096257952156660135, -1.073448835231765229960962579519779626699]]
 [[1.359304085971776420730660392788800573278, 1.359304085971776420730660392812266330097], [1.073448835231765229960962579519779626699, 1.07344883523176522996096257952156660135]]

2. The family

The coefficient $1$ of $x$ in the first equation and the radius $3$ in the second are now known up to $\pm e$. A system with interval coefficients is written with DynamicPolynomials:

@polyvar X Y
uncertain(q, e) = BigInterval(q; precision = 128) + BigInterval(e; precision = 128) * BigInterval(-1, 1)
family(e) = [Y^2 - X^3 + uncertain(1, e) * X, X^2 + Y^2 - uncertain(3, e)]
F = family(1//100)
2-element Vector{DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, BigInterval}}:
 ([0.9899999999999999999999999999999999999999926531603073607030751953966423609645136, 1.010000000000000000000000000000000000000007346839692639296924804603357639035486])X + Y² + -X³
 ([-3.010000000000000000000000000000000000000007346839692639296924804603357639035486, -2.989999999999999999999999999999999999999992653160307360703075195396642360964514]) + Y² + X²

3. Kantorovich, then Newton

At the center x0 of each exact solution, the Kantorovich test gives a radius $\rho$ and a number $\mu$; when $\mu < 1$, every system of the family has a unique solution in the ball of radius $\rho$ around x0. Interval Newton started from this ball shrinks it:

function enclose(F, sol; flags = 0)
    x0 = [BigFloat(mid(c); precision = 128) for c in sol]
    k = kantorovich_test(F, x0, 128, 64; flags = flags)
    mu, rho = right(k[2]), right(k[1])
    mu < 1 || return (mu = Float64(mu), box = nothing)
    _, box, _ = interval_newton(F, [BigInterval(c - rho, c + rho) for c in x0], 128, 128, 64, 30;
                                verbose = Cint(0), flags = flags)
    return (mu = Float64(mu), box = box)
end

enclosures = [enclose(F, s) for s in sols]
2-element Vector{@NamedTuple{mu::Float64, box::Vector{BigInterval}}}:
 (mu = 0.08120587234951192, box = [[1.35368646867140421028, 1.36492170327214863115], [-1.07883218480638047392, -1.06806548565714998602]])
 (mu = 0.08120587234951192, box = [[1.35368646867140421028, 1.36492170327214863115], [1.06806548565714998602, 1.07883218480638047392]])

Both solutions are certified: every system of the family has a solution in each box.

4. A check

The four systems at the corners of the family, $a = 1 \pm e$ and $c = 3 \pm e$, are solved exactly; each of their real solutions lies in one of the boxes:

e = 1//100
inbox(s, B) = all(is_inside(s[j], B[j]) for j in eachindex(B))
[(a, c, [any(inbox(s, E.box) for E in enclosures)
         for s in real_solutions([y^2 - x^3 + a * x, x^2 + y^2 - c]; prec = 64)])
 for a in (1 - e, 1 + e) for c in (3 - e, 3 + e)]
4-element Vector{Tuple{Rational{Int64}, Rational{Int64}, Vector{Bool}}}:
 (99//100, 299//100, [1, 1])
 (99//100, 301//100, [1, 1])
 (101//100, 299//100, [1, 1])
 (101//100, 301//100, [1, 1])

As $e$ grows, the boxes widen, and the certificate fails when $\mu$ reaches 1:

[(e, [E.mu for E in (enclose(family(e), s) for s in sols)]) for e in (1//10^6, 1//100, 1//10, 1//5)]
4-element Vector{Tuple{Rational{Int64}, Vector{Float64}}}:
 (1//1000000, [8.030002049513736e-6, 8.030002049513736e-6])
 (1//100, [0.08120587234951192, 0.08120587234951192])
 (1//10, [0.8971094607363902, 0.8971094607363902])
 (1//5, [1.9989361865884612, 1.9989361865884612])

5. More variables: Katsura-6

The Katsura-6 system has 32 real solutions (worksheet Real solutions of a zero-dimensional system). We give every coefficient a relative uncertainty $e$: the coefficient $c$ becomes $[c - e|c|, c + e|c|]$.

function katsura(n)
    K, u = polynomial_ring(QQ, ["u$i" for i in 0:n])
    U(i) = abs(i) <= n ? u[abs(i)+1] : zero(K)
    eqs = [sum(U(l) * U(m - l) for l in -n:n) - U(m) for m in 0:n-1]
    push!(eqs, sum(U(l) for l in -n:n) - 1)
    return eqs
end
k6 = katsura(6)
k6_sols = real_solutions(k6; prec = 64)
length(k6_sols)
32

The conversion of a polynomial of AbstractAlgebra into a polynomial with interval coefficients:

@polyvar v[1:7]
relative(p, e) = sum(uncertain(c, e * abs(c)) * prod(v[i]^m[i] for i in eachindex(m))
                     for (c, m) in zip(coefficients(p), exponent_vectors(p)))
relative(k6[end], 1//10^6)

\[ ([-1.000001000000000000000000000000000000000000000878176932010790960543050245092791, -0.999998999999999999999999999999999999999999999121823067989209039456949754907209]) + ([1.999997999999999999999999999999999999999999998243646135978418078913899509814418, 2.000002000000000000000000000000000000000000001756353864021581921086100490185582])v_{7} + ([1.999997999999999999999999999999999999999999998243646135978418078913899509814418, 2.000002000000000000000000000000000000000000001756353864021581921086100490185582])v_{6} + ([1.999997999999999999999999999999999999999999998243646135978418078913899509814418, 2.000002000000000000000000000000000000000000001756353864021581921086100490185582])v_{5} + ([1.999997999999999999999999999999999999999999998243646135978418078913899509814418, 2.000002000000000000000000000000000000000000001756353864021581921086100490185582])v_{4} + ([1.999997999999999999999999999999999999999999998243646135978418078913899509814418, 2.000002000000000000000000000000000000000000001756353864021581921086100490185582])v_{3} + ([1.999997999999999999999999999999999999999999998243646135978418078913899509814418, 2.000002000000000000000000000000000000000000001756353864021581921086100490185582])v_{2} + ([0.999998999999999999999999999999999999999999999121823067989209039456949754907209, 1.000001000000000000000000000000000000000000000878176932010790960543050245092791])v_{1} \]

Two boxes that overlap could hold the same solution of a system of the family: to count distinct solutions, the boxes must also be pairwise disjoint.

disjoint(A, B) = any(right(A[j]) < left(B[j]) || right(B[j]) < left(A[j]) for j in eachindex(A))

function certify(e; flags = 0)
    F = [relative(p, e) for p in k6]
    boxes = [E.box for E in (enclose(F, s; flags = flags) for s in k6_sols) if E.box !== nothing]
    overlaps = count(!disjoint(boxes[i], boxes[j]) for i in eachindex(boxes) for j in i+1:length(boxes))
    width = isempty(boxes) ? 0.0 : maximum(Float64(diam(c)) for B in boxes for c in B)
    return (e = e, certified = length(boxes), overlaps = overlaps, max_width = width)
end

[certify(e) for e in (1//10^12, 1//10^6, 1//10^5, 1//10^4, 1//10^3)]
5-element Vector{@NamedTuple{e::Rational{Int64}, certified::Int64, overlaps::Int64, max_width::Float64}}:
 (e = 1//1000000000000, certified = 32, overlaps = 0, max_width = 6.649556736366685e-9)
 (e = 1//1000000, certified = 32, overlaps = 0, max_width = 0.006965535808070613)
 (e = 1//100000, certified = 31, overlaps = 0, max_width = 0.08262257330400717)
 (e = 1//10000, certified = 19, overlaps = 0, max_width = 0.49477991162240065)
 (e = 1//1000, certified = 1, overlaps = 0, max_width = 0.008088522401773086)

With $e = 10^{-12}$ the 32 solutions are certified, but the widest box is already about $10^{-8}$ wide: some solutions amplify the uncertainty of the data by $10^4$. The number of certified solutions then drops: at $e = 10^{-3}$ a single one remains. A solution that is not certified may still exist for every system of the family; the test does not decide it.

6. Preconditioning options

kantorovich_test and interval_newton take a keyword flags, a sum of:

  • PACE.KANTO_HESSIAN: a bound independent of a diagonal scaling of the equations;
  • PACE.KANTO_PRECOND: the test is run on the system multiplied by $C$, an approximate inverse of the Jacobian at x0 (the inverse of its midpoint, in floating point at the working precision); the bound no longer depends on how the equations are written;
  • PACE.KANTO_SCHULZ, with PACE.KANTO_PRECOND: when $C$ is poor, it is improved by Schulz iterations $C \leftarrow C(2I - J C)$.

flags = 0 is the classical test, which returns $[2\eta, \mu]$. Otherwise the test returns $[2\eta, \mu, \mu_1, \mu_2, \delta]$: $\mu$ the best of the bounds computed, $\mu_1$ the one of KANTO_HESSIAN, $\mu_2$ the one of KANTO_PRECOND (Inf when not computed) and $\delta = \|I - CJ\|$, the defect of the preconditioner ($-1$ without it). In interval_newton, flags sets the Kantorovich test run at the starting point, whose constants are its third result; the Newton steps themselves do not depend on it.

Scaling the equations

Multiplying an equation by a constant changes neither the solutions nor the Newton iterates, but it changes the classical Kantorovich bound. With the second equation of section 2 multiplied by $s$:

scaled(s) = [F[1], s * F[2]]
[(s, enclose(scaled(s), sols[1]).mu) for s in (1, 1000, 10^6)]
3-element Vector{Tuple{Int64, Float64}}:
 (1, 0.08120587234951192)
 (1000, 7.413687480789473)
 (1000000, 7401.304198868133)

At $s = 1000$ the solution is no longer certified, although nothing has changed. With the preconditioned test:

[(s, enclose(scaled(s), sols[1]; flags = PACE.KANTO_PRECOND).mu) for s in (1, 1000, 10^6)]
3-element Vector{Tuple{Int64, Float64}}:
 (1, 0.018900049235540332)
 (1000, 0.018900049235540332)
 (1000000, 0.018900049235540332)

The bound no longer depends on $s$, and is smaller than the classical one even for $s = 1$.

An ill-conditioned system, from a point

Two curves that are almost the same near $(1, 1)$: their Jacobian there has a condition number of about $10^{18}$.

t = QQ(1, big(10)^18)
near = [x^2 + x + y - 3, x^2 + x + (1 + t) * y - 3 - t]
2-element Vector{AbstractAlgebra.Generic.MPoly{Rational{BigInt}}}:
 x^2 + x + y - 3
 x^2 + x + 1000000000000000001//1000000000000000000*y - 3000000000000000001//1000000000000000000

interval_newton is started from a point close to the solution, at a working precision of 64 bits, with each option:

x64 = [BigFloat(1; precision = 64) + BigFloat(2)^-30, BigFloat(1; precision = 64) - BigFloat(2)^-31]
options = ["classical" => 0, "HESSIAN" => PACE.KANTO_HESSIAN, "PRECOND" => PACE.KANTO_PRECOND,
           "PRECOND + SCHULZ" => PACE.KANTO_PRECOND | PACE.KANTO_SCHULZ]
function from_point(sys, x0, prec, flags)
    status, _, k = interval_newton(sys, x0, prec, prec, 100, 64, 30; verbose = Cint(0), flags = flags)
    return (status = status, mu = Float64(right(k[2])),
            delta = length(k) == 5 ? Float64(right(k[5])) : nothing)
end
[(name, from_point(near, x64, 64, f)...) for (name, f) in options]
4-element Vector{Tuple{String, Int64, Float64, Any}}:
 ("classical", 0, 1.178990514505501e28, nothing)
 ("HESSIAN", 0, 5.894952572527505e27, -1.0)
 ("PRECOND", 0, 0.7040000014752149, 0.375)
 ("PRECOND + SCHULZ", 0, 0.6400000013411045, 0.375)

The classical bound is about $10^{28}$; with the preconditioner it is below 1: the solution is certified, unique in the ball of radius $2\eta$ around the starting point. Here $\delta \approx 0.4 > 0.1$, so the Schulz iterations are run, and lower $\mu$ a little.

The Newton iteration itself fails: status 0, a singular interval Jacobian. Its steps solve the interval linear system by Gaussian elimination, without preconditioning, and at 64 bits the interval Jacobian on the box contains singular matrices. It needs more precision to conclude on its own:

[(prec, from_point(near, [BigFloat(c; precision = prec) for c in x64], prec,
                   PACE.KANTO_PRECOND | PACE.KANTO_SCHULZ).status) for prec in (64, 80, 128, 200)]
4-element Vector{Tuple{Int64, Int64}}:
 (64, 0)
 (80, 0)
 (128, 2)
 (200, 1)

At 200 bits Newton encloses the solution (status 1). At 64 bits the certificate is the ball of the preconditioned Kantorovich test.

Katsura-6

With uncertain coefficients, the starting balls of section 5 come from the Kantorovich test, and the preconditioned test certifies many more solutions when the uncertainty grows:

[(e, certify(e).certified, certify(e; flags = PACE.KANTO_PRECOND).certified)
 for e in (1//10^5, 1//10^4, 1//10^3)]
3-element Vector{Tuple{Rational{Int64}, Int64, Int64}}:
 (1//100000, 31, 32)
 (1//10000, 19, 31)
 (1//1000, 1, 14)

At $e = 10^{-3}$, 14 solutions are certified instead of one, in boxes that are still pairwise disjoint:

certify(1//10^3; flags = PACE.KANTO_PRECOND)
(e = 1//1000, certified = 14, overlaps = 0, max_width = 1.9181061907287857)

Some of them are wide (up to about 2): for these solutions an uncertainty of $10^{-3}$ on the data leaves a large uncertainty on the result, but each box still holds exactly one solution of every system of the family.

7. What this approach cannot see

The method follows the real solutions of the exact instance. It proves nothing about the solutions the instance does not have. The line $y = c$ and the unit circle, with $c = 1 + 1/1000$, do not meet:

c0 = 1 + 1//1000
real_solutions([x^2 + y^2 - 1, y - c0]; prec = 64)
Vector{BigInterval}[]

Nothing to enclose. Yet with an uncertainty $e = 1/100$ on $c$, the systems with $c \le 1$ have real solutions: at the corner $c = c_0 - e$,

real_solutions([x^2 + y^2 - 1, y - (c0 - 1//100)]; prec = 64)
2-element Vector{Vector{BigInterval}}:
 [[-0.1338618691039386316546826372555798740267, -0.1338618691039386316546826372555082097454], [0.9909999999999999999999999999989788563097, 0.9910000000000000000000000000010242575049]]
 [[0.1338618691039386316546826372555082097454, 0.1338618691039386316546826372555798740267], [0.9909999999999999999999999999989788563097, 0.9910000000000000000000000000010242575049]]

Two complex solutions of the instance became real. To describe all the systems of a family, the uncertain coefficients must be treated as parameters, and the parameter space split into regions where the number of real solutions is constant: this is the subject of the worksheet on systems depending on parameters.