Constraint propagation for univariate quadratic constraints

Arnold Neumaier

Institut für Mathematik, Universität Wien Strudlhofgasse 4, A-1090 Wien, Austria email: Arnold.Neumaier@univie.ac.at WWW: http://www.mat.univie.ac.at/∼neum/

January 5, 2005

Abstract.

We present formulas for rigorous constraint propagation of quadratic equality or inequality constraints involving a single nonlinear variable. Since the analysis is very elementary, probably everything in here has been known for a long time. The present approach, based on directed rounding only, provides efficient alternatives to the procedures discussed by Hansen & Walster [1] (who only treat the solution of a quadratic equation with interval coefficients), which employ interval arithmetic.

In view of pending patent applications by these authors, who by these activities threaten to curb the freedom of research on interval methods, the following is explicitly stated:

Various modifications to the methods described will be readily apparent to those skilled in the art, and the general principles defined herein may be applied to such modifications without departing from the spirit and scope of the present methods. Thus, the present methods are not intended to be limited to the formulas shown, but are to be accorded the widest scope consistent with the principles and features disclosed herein.

Notation is as in my book Neumaier [2].

1 Bounds for quadratic expressions

To find a rigorous upper bound on

u = sup {ax^{2} + bx | x ∈ x},

we note that

u = max \{\underline{x}(a\underline{x} + b), \overline{x}(a\overline{x} + b)\},

except in case that ax^{2} + bx attains its global maximum in the interior of x. This is the case iff a < 0 and t = −b/(2a) is in the interior of x, in which case u = b^{2}/(−4a), attained at t.

If \underline{x} ≥ 0, we get a rigorous upper bound in finite precision arithmetic by computing with upward rounding as follows (xl = \underline{x}, xu = \overline{x}):

  roundup;
  u=max(xl*(a*xl+b),xu*(a*xu+b));
  s=b/2; t=s/(-a);
  if t>xl, r=(-2*a)*xu;
    if r>b, u=max(u,s*t); end;
  end;

With some extra analysis, it could be determined in most cases which of the three cases is the worst case; however, if the unconstrained maximum of the quadratic is very close to a bound (or to both bounds), two (or three) of the cases might apply due to uncertainty caused by rounding errors.

Finding a rigorous enclosure for the interval

c = sup {ax^{2} + bx | x ∈ x, a ∈ a, b ∈ b}

can be reduced to the above for \underline{x} ≥ 0, using

\overline{c} = sup \{\overline{a}x^{2} + \overline{b}x | x ∈ x\}, \underline{c} = − sup \{−\underline{a}x^{2} − \underline{b}x | x ∈ x\}.

The case \overline{x} ≤ 0 can be reduced to this by changing the sign of x, and the general case by splitting x at zero if necessary.

Essentially the same analysis holds for rigorous upper bounds on

u = sup { ∑_{i=1}^{n} a_{i}x^{i} | x ∈ x }

and for rigorous enclosures of

c = sup { ∑_{i=1}^{n} a_{i}x^{i} | x ∈ x, a ∈ a },

except that finding the interior extrema is more involved. It can be done with closed formulas for n ≤ 5 (though already n = 4 is quite cumbersome and not recommended), and in general (recommended for n > 3) using a root enclosure algorithm for the derivative, such as that in Neumaier [3].

2 Solving quadratic constraints

To find the set

X = {x ≥ 0 | ax^{2} + 2bx ≥ c},

we proceed as follows. If a = 0, the constraint is in fact linear, and we have

X = \{ \begin{matrix}∅ & if c > 0, b ≤ 0, \\ [0.5c/b, ∞] & if c > 0, b > 0, \\ [0, 0.5c/b] & if c ≤ 0, b < 0, \\ [0, ∞] & if c ≤ 0, b ≥ 0, \end{matrix}

which can be nested such that only two compares are needed in any particular case. For a rigorous enclosure in finite precision arithmetic, rounding must be downwards in the second case, and upwards in the third case.

If a ≠ 0, the behavior is governed by the zeros of the quadratic equation ax^{2} + 2bx − c = 0, given by

t_{1} = \frac{−b − \sqrt{Δ}}{a} = \frac{c}{b − \sqrt{Δ}}, t_{2} = \frac{−b + \sqrt{Δ}}{a} = \frac{c}{b + \sqrt{Δ}},

where Δ := b^{2} + ac. If Δ ≥ 0, the zeros are real and the nonnegative zeros determine

X = \{ \begin{matrix}[0, ∞] \ ]t_{1}, t_{2}[ & if a > 0, \\ [0, ∞] ∩ [t_{2}, t_{1}] & if a < 0. \end{matrix}

Depending on the signs of a, b and c we find

X = \{ \begin{matrix}∅ & if a < 0, b ≤ 0, c > 0, \\ [0, −(c/z)] & if a < 0, b ≤ 0, c ≤ 0, \\ [0, z/(−a)] & if a < 0, b ≥ 0, c ≤ 0, \\ [−((−c)/z), z/(−a)] & if a < 0, b ≥ 0, c > 0, \\ [0, −(c/z)] ∪ [z/a, ∞] & if a > 0, b ≤ 0, c ≤ 0, \\ [z/a, ∞] & if a > 0, b ≤ 0, c > 0, \\ [−((−c)/z), ∞] & if a > 0, b ≥ 0, c > 0, \\ [0, ∞] & if a > 0, b ≥ 0, c ≤ 0, \end{matrix}

where

z = |b| + \sqrt{Δ}.

These formulas are numerically stable, and can be nested such that only three compares are needed in any particular case. (There are avoidable overflow problems for huge |b|, which can be cured by using for huge |b| instead of \sqrt{b^{2} + ac} the formula |b|\sqrt{1 + ac/b^{2}}.)

Rigorous results in the presence of rounding errors are obtained if lower bounds are rounded downwards, and upper bounds are rounded upwards. With the bracketing as given, this happens if in cases 2,5 and 6 all computations (including those of Δ = \sqrt{b^{2} + ac} and z = |b| + \sqrt{Δ}) are done with rounding downwards, and in the other cases with rounding upwards. (However, this does not hold for the version guarded against overflow, where further care is needed for the directed rounding of \sqrt{Δ} = |b|\sqrt{1 + ac/b^{2}}.)

If (the exact) Δ is negative, there is no real solution, and X is empty if c > 0 and [0, ∞] otherwise. The case when the sign of Δ cannot be determined due to rounding errors needs special consideration. In the first and last case, the conclusion holds independent of the sign of Δ, so that the latter need only be computed for cases 2–7. In the cases 2, 3, 6, and 7 we have ac ≥ 0, so that Δ ≥ 0 automatically. This leaves cases 4 and 5. Now it is easily checked that with the recommended rounding and, in place of cases 4 and 5,

X = \{ \begin{matrix}∅ & if a < 0, b ≥ 0, c > 0, Δ < 0, \\ [−((−c)/z), z/(−a)] & if a < 0, b ≥ 0, c > 0, Δ ≥ 0, \\ [0, −(c/z)] ∪ [z/a, ∞] & if a > 0, b ≤ 0, c ≤ 0, Δ ≥ 0, \\ [0, ∞] & if a > 0, b ≤ 0, c ≤ 0, Δ < 0, \end{matrix}

a rigorous enclosure is computed in all cases.

Finding the set

X′ = {x ≥ 0 | ax^{2} + 2bx ∈ c for some a ∈ a, b ∈ b}

can be reduced to the previous task since

X′ = \{x ≥ 0 | \underline{a}x^{2} + 2\underline{b}x ≤ \overline{c}\} ∩ \{x ≥ 0 | \underline{a}x^{2} + 2\underline{b}x ≤ \overline{c}\}.

The sets

X′′ = {x ∈ x_{0} | ax^{2} + 2bx ≥ c}

and

X′′′ = {x ∈ x_{0} | ax^{2} + 2bx ∈ c for some a ∈ a, b ∈ b}

can be obtained by intersecting the result of the above tasks with x_{0} if \underline{x}_{0} ≥ 0, by negating x, x_{0}, and b if \overline{x}_{0} ≤ 0, and by splitting x_{0} at zero if 0 is in the interior of x_{0}. By modifying the code appropriately, one can also avoid computing roots which can be seen to lie outside x_{0}.

With minor changes, these formulas also apply for strict inequalities and interior enclosures. Also, it is clear that polynomial inequalities and inclusions of interval polynomials can be solved by a straightforward adaptation of the above arguments.

References

[1] E. R. Hansen and G. W. Walster, Sharp bounds on interval polynomial roots, Reliable Computing 8 (2002), 115–122.

[2] A. Neumaier, Interval Methods for Systems of Equations, Cambridge Univ. Press, Cambridge 1990.

[3] A. Neumaier, Enclosing clusters of zeros of polynomials, J. Comput. Appl. Math. 156 (2003), 389–401.