One question, four paradigms: is this polynomial Hurwitz?

Download: Jupyter notebook, Julia script, Maple version (same conclusions). To install PACE and run them, see Running the worksheets.

A characteristic polynomial with exact coefficients is given. Are all its roots in the open left half-plane? We answer with floating-point arithmetic, exact arithmetic, interval arithmetic and certified root isolation, and compare.

The example is built to be hard for floating point: five times the same pair of roots, close to the imaginary axis, $q(s) = \bigl(s^2 + a\,s + 1\bigr)^5, \qquad a = 2^{-10}.$ The answer is yes: every root has real part $-a/2 = -2^{-11} \approx -4.9\cdot 10^{-4}$. Every coefficient of $q$ is exactly representable in Float64, so no error comes from reading the data.

using LinearAlgebra, Nemo, PACE, PaceControl

a = 1 // 2^10
R, s = polynomial_ring(QQ, "s")
q = (s^2 + a * s + 1)^5
ce = Rational{BigInt}.(reverse([coeff(q, k) for k in 0:degree(q)]))   ## highest degree first
@assert Rational{BigInt}.(Float64.(ce)) == ce                          ## exact in Float64

1. Floating point

The usual way: the eigenvalues of the companion matrix, which is what roots does in Matlab.

function max_real_part(coefs)
    n = length(coefs) - 1
    C = zeros(n, n)
    C[1, :] = -Float64.(coefs[2:end]) ./ Float64(coefs[1])
    for i in 2:n
        C[i, i-1] = 1
    end
    maximum(real.(eigvals(C)))
end

max_real_part(ce)

6.580252729979086e-5

The largest real part is positive: the computed answer is unstable, and it is wrong. A root of multiplicity $m$ is perturbed by about $\varepsilon^{1/m}$, here $(2^{-53})^{1/5}\approx 6\cdot 10^{-4}$, more than the distance $2^{-11}$ to the axis.

Changing the algorithm does not help. The Hurwitz criterion asks that the leading principal minors of the Hurwitz matrix be positive. Computed on the same Float64 coefficients, by elimination without pivoting, some come out negative:

function leading_minors(H)
    H = copy(H)
    m = size(H, 1)
    D = similar(H, m)
    p = one(eltype(H))
    for k in 1:m
        p *= H[k, k]
        D[k] = p
        for i in k+1:m
            f = H[i, k] / H[k, k]
            for j in k:m
                H[i, j] -= f * H[k, j]
            end
        end
    end
    D
end

D = leading_minors(PACE.hurwitz_matrix(Float64.(ce)))
(sign.(D), all(>(0), D))

([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0], false)

The answer is unstable again, and wrong: the exact minors are all positive, the last ones of order $10^{-43}$ and $10^{-70}$. A Routh table in Float64 happens to answer correctly here, with pivots of order $10^{-27}$, but nothing tells which one to trust.

2. Exact arithmetic

The same criterion on rational numbers decides without any tolerance:

all(>(0), leading_minors(PACE.hurwitz_matrix(ce)))

true

hurwitz_stable of PaceControl (Liénard-Chipart criterion, exact) agrees:

hurwitz_stable(ce)

true

3. Interval arithmetic with adaptive precision

Each number is replaced by an interval that contains it, and each operation rounds outwards. A sign is decided when the interval does not contain $0$. We compute the leading principal minors of the Hurwitz matrix by elimination without pivoting.

function hurwitz_minors(coefs, prec)
    setprecision(BigFloat, prec) do
        H = PACE.hurwitz_matrix([BigInterval(BigFloat(x, RoundDown), BigFloat(x, RoundUp)) for x in coefs])
        m = size(H, 1)
        D = Vector{BigInterval}(undef, m)
        p = one(BigInterval)
        for k in 1:m
            p *= H[k, k]
            D[k] = p
            for i in k+1:m
                f = H[i, k] / H[k, k]
                for j in k:m
                    H[i, j] -= f * H[k, j]
                end
            end
        end
        D
    end
end

signs(D) = join(has_zero(d) ? "?" : (d > 0 ? "+" : "-") for d in D)

function decide(coefs)
    trace = String[]
    prec = 53
    while true
        D = hurwitz_minors(coefs, prec)
        push!(trace, lpad(prec, 4) * " bits: " * signs(D))
        all(d -> !has_zero(d), D) && return all(d -> d > 0, D), trace
        prec *= 2
    end
end

stable, trace = decide(ce)
trace

3-element Vector{String}:
 "  53 bits: ++++??????"
 " 106 bits: ++++++????"
 " 212 bits: ++++++++++"

At low precision some signs are undecided ("?"): intervals do not give a wrong answer, they say they do not know. Doubling the precision until every sign is decided gives the right answer, at a cost that adapts to the difficulty of the input.

4. Certified root isolation

The question can also be answered by locating the roots. Exact arithmetic first removes the multiplicities: the square-free part $q/\gcd(q,q')$ has the same roots, each simple.

qs = squarefree_part(q)

Its roots are then isolated in certified complex balls (PWPoly, through PACE). With integer coefficients (not on Windows, where PWPoly is not available):

den = lcm([denominator(coeff(qs, k)) for k in 0:degree(qs)])
ci = [BigInt(numerator(coeff(qs, k) * den)) for k in 0:degree(qs)]
rts, = isolate_complex_roots_sqrfree_poly(ci; log2eps = -80, acb = true)

(Nemo.ComplexFieldElem[[-0.00048828125000000000000 +/- 1e-28] + [0.99999988079070334379 +/- 5.46e-22]*im, [-0.00048828125000000000000 +/- 1e-28] + [-0.99999988079070334379 +/- 5.46e-22]*im], Bool[0, 0], 0)

all(r -> Nemo.is_negative(real(r)), rts)

true

Each ball contains exactly one root and lies in the open left half-plane: the polynomial is Hurwitz, and we also know where its roots are, and with which multiplicity.

5. When the data are not floating-point numbers

Take now $\delta = 2^{-20}$ and $p(s) = \bigl((s+\delta)^2+1\bigr)^4,$ Hurwitz again (real parts $-\delta$). Its coefficients need up to 161 bits:

p = ((s + 1 // 2^20)^2 + 1)^4
cp = Rational{BigInt}.(reverse([coeff(p, k) for k in 0:degree(p)]))
[ndigits(denominator(x), base = 2) for x in cp]

9-element Vector{Int64}:
   1
  18
  39
  58
  80
  98
 119
 138
 161

Reading them as Float64 changes the polynomial, and the exact test shows that the rounded polynomial is not Hurwitz:

cr = Rational{BigInt}.(Float64.(cp))
(hurwitz_stable(cp), hurwitz_stable(cr))

(true, false)

Here no floating-point algorithm, however good, can give the right answer: the question is lost before any computation. This is why the inputs of the exact and certified tools are rational numbers, and why a decimal coefficient such as $0.1$ should be read as $1/10$.

6. When the data are not known exactly

Coefficients coming from measurements are only known to lie in intervals. The question becomes: is every polynomial of the family Hurwitz? Take $(s+1)^6$, whose coefficients are $1, 6, 15, 20, 15, 6, 1$, and let each of them vary independently by $\pm 12.5\%$ ($\delta = 1/8$):

c6 = Rational{BigInt}[1, 6, 15, 20, 15, 6, 1]
δ = 1 // 8
lo, hi = c6 .* (1 - δ), c6 .* (1 + δ)

(Rational{BigInt}[7//8, 21//4, 105//8, 35//2, 105//8, 21//4, 7//8], Rational{BigInt}[9//8, 27//4, 135//8, 45//2, 135//8, 27//4, 9//8])

A first idea is to compute the Hurwitz minors with these interval coefficients, by the elimination of section 1. If all of them are positive intervals, every polynomial of the family is Hurwitz. Here some signs stay undecided, and raising the precision does not help: the width comes from the data, and interval arithmetic overestimates it.

interval_signs(prec) = setprecision(BigFloat, prec) do
    signs(leading_minors(PACE.hurwitz_matrix(
        [BigInterval(BigFloat(l, RoundDown), BigFloat(h, RoundUp)) for (l, h) in zip(lo, hi)])))
end
(interval_signs(53), interval_signs(1000))

("+++???", "+++???")

Kharitonov's theorem (1978): when the coefficients vary independently and the interval of the leading coefficient does not contain $0$, every polynomial of the family is Hurwitz if and only if four particular polynomials are. Their coefficients are bounds of the intervals: from the constant term upwards, lower ($l$) or upper ($u$) bounds in the patterns $l\,l\,u\,u$, $l\,u\,u\,l$, $u\,l\,l\,u$ and $u\,u\,l\,l$, repeated every four powers. With rational bounds, they are tested exactly.

function kharitonov(lo, hi)
    n = length(lo)
    pats = ((0, 0, 1, 1), (0, 1, 1, 0), (1, 0, 0, 1), (1, 1, 0, 0))
    [[pat[(n - i) % 4 + 1] == 0 ? lo[i] : hi[i] for i in 1:n] for pat in pats]
end

[hurwitz_stable(k) for k in kharitonov(lo, hi)]

4-element Vector{Bool}:
 1
 1
 1
 1

The four polynomials are Hurwitz: every polynomial of the family is, which the interval minors could not decide. PACE.hurwitz_stable on interval coefficients does the same:

PACE.hurwitz_stable([BigInterval(BigFloat(l), BigFloat(h)) for (l, h) in zip(lo, hi)])

true

With $\pm 25\%$ ($\delta = 1/4$), two of the four polynomials are not Hurwitz. They belong to the family, so the family is not robustly stable, and this is again a proof:

[hurwitz_stable(k) for k in kharitonov(c6 .* (1 - 1 // 4), c6 .* (1 + 1 // 4))]

4-element Vector{Bool}:
 1
 0
 1
 0

Summary

And when the data are rounded before the computation, even an exact method answers another question. When they are only known to lie in intervals, the minors computed in interval arithmetic may not decide, and Kharitonov's theorem reduces the question to four exact tests.