\documentstyle[11pt]{article} \setlength{\oddsidemargin}{0in} \setlength{\evensidemargin}{0in} \setlength{\textwidth}{6.5in} \setlength{\topmargin}{0in} \setlength{\textheight}{8.5in} \newcommand{\bigsqcap}{ \mbox{\Large $ \sqcap$}} \hbadness=10000 \newtheorem{definition_int}{Definition} \newenvironment{definition}{\begin{definition_int}\begin{rm}}% {\end{rm}\end{definition_int}} \newtheorem{specification_int}{Specification} \newenvironment{specification}{\begin{specification_int}\begin{rm}}% {\end{rm}\end{specification_int}} \newtheorem{implementation_int}{Implementation} \newenvironment{implementation}{\begin{implementation_int}\begin{rm}}% {\end{rm}\end{implementation_int}} \newtheorem{proposition_int}{Proposition} \newenvironment{proposition}{\begin{proposition_int}\begin{rm}}% {\end{rm}\end{proposition_int}} \newtheorem{example_int}{Example} \newenvironment{example}{\begin{example_int}\begin{rm}}% {\end{rm}\end{example_int}} \newtheorem{lemma_int}{Lemma} \newenvironment{lemma}{\begin{lemma_int}\begin{rm}}% {\end{rm}\end{lemma_int}} \newtheorem{remark_int}{Remark} \newenvironment{remark}{\begin{remark_int}\begin{rm}}% {\end{rm}\end{remark_int}} \newenvironment{proof}{\noindent {\bf Proof }}% {$\Box$ \vspace{.4cm}} \newtheorem{theorem_int}{Theorem} \newenvironment{theorem}{\begin{theorem_int}\begin{rm}}% {\end{rm}\end{theorem_int}} \newcommand{\separe}[1]{\noindent \begin{minipage}{1cm} \setlength{\unitlength}{1cm} \begin{picture}(.1,.1)(0,0) \put(.4,0){\line(1,0){15}} \end{picture} \end{minipage} } \newcommand{\R}{\Re} \newcommand{\Rn}{\Re^{n}} \newcommand{\Ri}{\Re^{\infty}} \newcommand{\Rp}{\Re^{+}} \newcommand{\F}{{\cal F}} \newcommand{\I}{{\cal I}} \newcommand{\In}{{\cal I}^{n}} \newcommand{\ap}[1]{\overline{{#1}}} \newcommand{\floor}[1]{\lfloor {#1} \rfloor} \newcommand{\ceil}[1]{\lceil {#1} \rceil} \newcommand{\LEFT}{{\it left}} \newcommand{\RIGHT}{{\it right}} \newcommand{\CENTER}{{\it center}} \newcommand{\fto}{\rightarrow} \newcommand{\lbra}{{\rm [\hspace{-.02in}[}} \newcommand{\rbra}{{\rm ]\hspace{-.02in}]}} \newcommand{\sem}[1]{\lbra {#1} \rbra} \newcommand{\Snseq}{D_1,\ldots,D _n} \newcommand{\xnseq}{x_1,\ldots,x_n} \newcommand{\rnseq}{r_1,\ldots,r_n} \newcommand{\Xnseq}{x_1,\ldots,x_n} \newcommand{\Inseq}{I_1,\ldots,I_n} \newcommand{\IInseq}{I_1,\ldots,I_{i-1},\ap{r_i},I_{i+1},\ldots,I_n} \newcommand{\Sonseq}{D^0_1,\ldots,D^0_n} \newcommand{\Ionseq}{I^0_1,\ldots,I^0_n} \newcommand{\Pmseq}{P_1,\ldots,P_m} \newcommand{\IE}{{\it IE}} \newcommand{\NE}{{\it NE}} \newcommand{\EF}{{\it EF}} \newcommand{\EE}{{\it EE}} \newcommand{\TE}{{\it TE}} \newcommand{\DE}{{\it DE}} \newcommand{\newton}{{\tt Newton}} \newcommand{\Bool}{{\it Bool}} \newcommand{\bo}[1]{{\it box}\{{#1}\}} \newcommand{\Constraint}{{\it Constraint}} \newcommand{\Function}{{\it Function}} \newcommand{\CONSTRAINT}{{\it CONSTRAINT}} \newcommand{\FUNCTION}{{\it FUNCTION}} \newcommand{\UDIV}{{\it UDIV}} \newcommand{\CS}{{\cal S}} \newcommand{\CON}{{\it cond}} %____________________________________________________________________ \begin{document} \Large \noindent {\bf Solving Polynomial Systems Using a Branch and Prune Approach} \normalsize \\ \\ \parbox{2.1in}{ \begin{center} {\bf P. Van Hentenryck} \\ Brown University \\ Box 1910 \\ Providence, RI 02912 \\ {\tt pvh@cs.brown.edu} \\ \end{center} } \ \ \parbox{2.1in}{ \begin{center} {\bf D. McAllester} \\ MIT AI Lab \\ Technology Square, 545 \\ Cambridge, USA \\ {\tt dam@ai.mit.edu} \\ \end{center} } \ \ \parbox{2.3in}{ \begin{center} {\bf D. Kapur} \\ SUNY at Albany \\ Dep. of Computer Science \\ Albany, NY-12222 \\ {\tt kapur@cs.albany.edu} \end{center} } \section*{Abstract} This paper presents {\tt Newton}, a branch \& prune algorithm to find all isolated solutions of a system of polynomial constraints. {\tt Newton} can be characterized as a global search method which uses intervals for numerical correctness and for pruning the search space early. The pruning in {\tt Newton} consists in enforcing at each node of the search tree a unique local consistency condition, called box-consistency, which approximates the notion of arc-consistency well-known in artificial intelligence. Box-consistency is parametrized by an interval extension of the constraint and can be instantiated to produce Hansen-Segupta narrowing operator (used in interval Newton methods) as well as new operators which are more effective when the computation is far from a solution. {\tt Newton} has been evaluated on a variety of benchmarks from kinematics, chemistry, combustion, economics, and mechanics. On these benchmarks, it outperforms the interval Newton methods we are aware of and compares well with state-of-the-art continuation methods. Limitations of {\tt Newton} (e.g. a sensibility to the size of the initial intervals on some problems) are also discussed. Of particular interest is the mathematical and programming simplicity of the method. \section{Introduction} Many applications in science and engineering (e.g. chemistry, robotics, economics, mechanics) require finding all isolated solutions to a system of polynomial constraints over real numbers. This problem is difficult due to its inherent computational complexity (i.e., it is NP-hard) and due to the numerical issues involved to guarantee correctness (i.e., finding all solutions) and to ensure termination. Several interesting methods have been proposed in the past for this task, including two fundamentally different methods: interval Newton methods (e.g. \cite{Moore66,Hansen81,Hansen83,Krawczyk69}) and continuation methods (e.g. \cite{Morgan87,Cools94}). Continuation methods have been shown to be effective for problems for which the total degree is not too high, since the number of paths explored depends on the estimation of the number of solutions. Interval Newton methods are generally robust but tend to be slow. The purpose of this paper is to propose and to study a novel algorithm called {\tt Newton}. From a user standpoint, {\tt Newton} receives as input a system of polynomial constraints over, say, variables $x_1,\ldots,x_n$ and a box, i.e., an interval tuple $\langle I_1,\ldots,I_n \rangle$ specifying the initial range of these variables; it returns a set of boxes of specified accuracy containing all solutions. Operationally, {\tt Newton} is a branch \& prune algorithm which was inspired by the traditional branch and bound approach used to solve combinatorial optimization problems. {\tt Newton} uses intervals to address the two fundamental problems listed above. Numerical reliability is obtained by evaluating functions over intervals using outward rounding (as in interval Newton methods). The complexity issue is addressed by using constraints to reduce the intervals early in the search. The pruning in {\tt Newton} is achieved by enforcing a unique local consistency condition, called box-consistency, at each node of the search tree. Box-consistency is an approximation of arc-consistency, a notion well-known in artificial intelligence \cite{Mackworth77,Montanari74} and used to solve discrete combinatorial problems in several systems (e.g., \cite{ANNALSOR,book}). Box-consistency is parametrized by an interval extension operator for the constraint and can be instantiated to produce various narrowing operators. In particular, box-consistency on the Taylor extension of the constraint produces a generalization of Hansen-Segupta operator \cite{Hansen81} well-known in interval Newton methods. In addition, box-consistency on the natural extension produces narrowing operators which are more effective when the algorithm is not near a solution. {\tt Newton} has the following properties: \begin{itemize} \item {\bf Correctness:} {\tt Newton} finds all isolated solutions to the system in the following sense: if $\langle v_1,\ldots,v_n \rangle$ is a solution, then {\tt Newton} returns at least one box $\langle I_1',\ldots,I_n' \rangle$ such that $v_i \in I_i \;\; (1 \leq i \leq n)$. In addition, {\tt Newton} may guarantee the existence of a unique solution in some or all the boxes in the result. \item {\bf Termination:} {\tt Newton} always terminates in finite time. \item {\bf Effectiveness:} {\tt Newton} has been evaluated on a variety of benchmarks from kinematics, chemistry, combustion, economics, and mechanics. It outperforms the interval Newton methods we are aware of and compares well with state-of-the-art continuation methods on many problems. Interestingly, {\tt Newton} solves the Broyden banded function problem \cite{Hansen81} and Mor\'{e}-Cosnard discretization of a nonlinear integral equation \cite{More79} for several hundred variables. \item {\bf Simplicity and Uniformity:} {\tt Newton} is based on simple mathematical results and is easy to use and to implement. It is also based on a single concept: box-consistency. \end{itemize} The rest of this paper is organized as follows. Section \ref{overview} gives an overview of the approach. Section \ref{background} contains the preliminaries. Section \ref{pruning} presents an abstract version of the branch \& prune algorithm. Section \ref{implementation} discusses the implementation of the algorithm. Section \ref{experimental} describes the experimental results. Section \ref{related} discusses related work and the development of the ideas presented here. Section \ref{conclusion} concludes the paper. \section{Overview of The Approach} \label{overview} As mentioned, {\tt Newton} is a global search algorithm which solves a problem by dividing it into subproblems which are solved recursively. In addition, {\tt Newton} is a branch \& prune algorithm which means that it is best viewed as an iteration of two steps \begin{enumerate} \item pruning the search space; \item making a nondeterministic choice to generate two subproblems \end{enumerate} until one or all solutions to a given problem are found. The pruning step is responsible to make sure that some local consistency conditions are satisfied. It consists of reducing the intervals associated with the variables so that every constraint appears to be locally consistent. The local consistency condition of {\tt Newton} is called box-consistency, an approximation of arc-consistency, a notion well-known in artificial intelligence \cite{Mackworth77,Montanari74} and used in many systems (e.g. \cite{book,CCFD,SCREAMER}) to solve discrete combinatorial search problems. Informally speaking, a constraint is arc-consistent if for each value in the range of a variable there exist values in the ranges of the other variables such that the constraint is satisfied. {\tt Newton} approximates arc-consistency which cannot be computed on real numbers in general. The pruning step either fails, showing the absence of solution in the intervals, or succeeds in enforcing the local consistency condition. Sometimes, local consistency also implies global consistency as in the case of the Broyden banded function, i.e. \[ f_i(x_1,\ldots,x_n) = x_i (2 + 5 x_i^2) + 1 - \sum_{j \in J_i} x_j (1 + x_j) \;\;\; (1 \leq i \leq n) \] where $J_i = \{ j \ \mid \ j \neq i \ \& \ {\it max}(1,i-5) \leq j \leq {\it min}(n,i+1)\}$. The pruning step of {\tt Newton} solves this problem in essentially linear time for initial intervals of the form $[-10^8,10^8]$ and always proves the existence of a solution in the final box. However, in general, local consistency does not imply global consistency either because they are multiple solutions or simply because the local consistency condition is too weak. Consider the intersection of a circle and of a parabola: \[ \left\{ \begin{array}{c} x_1^2 + x_2^2 = 1 \\ x_1^2 + x_2 = 0 \end{array} \right. \] with initial intervals in $[-10^8,10^8]$. The pruning step returns the intervals \begin{quote} $[-1.0000000000012430057,+1.0000000000012430057] $ \\ $[-0.0000000000000000000,+1.0000000000012430057] $ \end{quote} with 8 digits of accuracy. Branching on $x_1$ directly leads to the two solutions \begin{quote} $[-0.7861513777574234974,-0.7861513777574231642] $ \\ $[+0.6180339887498946804,+0.6180339887498950136] $ \end{quote} and \begin{quote} $[+0.7861513777574231642,+0.7861513777574233864] $ \\ $[+0.6180339887498946804,+0.6180339887498950136]$. \end{quote} Note that, in this case, {\tt Newton} makes the smallest number of choices to isolate the solutions. To conclude this motivating section, let us illustrate {\tt Newton} on a larger example which describes the inverse kinematics of an elbow manipulator \cite{Hong95}: \[ \left\{ \begin{array}{l} s_2 c_5 s_6 - s_3 c_5 s_6 - s_4 c_5 s_6 + c_2 c_6 + c_3 c_6 + c_4 c_6 = 0.4077 \\ c_1 c_2 s_5 + c_1 c_3 s_5 + c_1 c_4 s_5 + s_1 c_5 = 1.9115\\ s_2 s_5 + s_3 s_5 + s_4 s_5 = 1.9791\\ c_1 c_2 + c_1 c_3 + c_1 c_4 + c_1 c_2 + c_1 c_3 + c_1 c_2 = 4.0616\\ s_1 c_2 + s_1 c_3 + s_1 c_4 + s_1 c_2 + s_1 c_3 + s_1 c_2 = 1.7172\\ s_2 + s_3 + s_4 + s_2 + s_3 + s_2 = 3.9701\\ s_i^2 + c_i^2 = 1 \;\;\;\;\;\; (1 \leq i \leq 6). \end{array} \right. \] and assumes that the initial intervals are in $[-10^8,10^8]$ again. The pruning step returns the intervals \begin{quote} $[-1.0000000000000000000,+1.0000000000000000000]$ \\ $[-1.0000000000000000000,+1.0000000000000000000]$ \\ $[+0.3233666666666665800,+1.0000000000000000000]$ \\ $[-1.0000000000000000000,+1.0000000000000000000]$ \\ $[-0.0149500000000000189,+1.0000000000000000000]$ \\ $[-1.0000000000000000000,+1.0000000000000000000]$ \\ $[-0.0209000000000001407,+1.0000000000000000000]$ \\ $[-1.0000000000000000000,+1.0000000000000000000]$ \\ $[+0.6596999999999998420,+1.0000000000000000000]$ \\ $[-0.7515290480087772896,+0.7515290480087772896]$ \\ $[-1.0000000000000000000,+1.0000000000000000000]$ \\ $[-1.0000000000000000000,+1.0000000000000000000]$ \end{quote} showing already some interesting pruning. After exactly 12 branchings and in less than a second, {\tt Newton} produces the first box with a proof of existence of a solution in the box. \section{Preliminaries} \label{background} In this section, we review some basic concepts needed for this paper, including interval arithmetic and the representation of constraints. More information on interval arithmetic can be found in many places (e.g., \cite{Alefeld83,Hansen81,Hansen83,Moore66,Moore79}). Our definitions are slightly non-standard. \subsection{Interval Arithmetic} We consider $\Ri = \R \cup \{-\infty,\infty\}$ the set of real numbers extended with the two infinity symbols and the natural extension of the relation $<$ to this set. We also consider a finite subset $\F$ of $\Ri$ containing $-\infty,\infty,0$. In practice, $\F$ corresponds to the floating-point numbers used in the implementation. \begin{definition}[Interval] An interval $[a,b]$ with $a,b \in \F$ is the set of real numbers $$\{ r \in \R \ \mid \ a \leq r \leq b \}.$$ The set of intervals is denoded by $\I$ and is ordered by set inclusion. \end{definition} \begin{definition}[Approximation] Let $S$ be a subset of $\R$. The approximation of $S$, denoted by $\ap{S}$ or ${\it box\{S\}}$, is the smallest interval $I$ such that $S \subseteq I$. We often write $r$ instead of $\ap{\{r\}}$. We also write $I_1 \uplus I_2$ to denote ${\it box}\{I_1 \cup I_2\}$. \end{definition} \noindent In the following, we denote real numbers by the letters $r,v$, $\F$-numbers by the letters $a,b,l,m,u$, intervals by the letter $I$, real functions by the letters $f,g$ and interval functions by the letters $F,G$, all possibly subscripted. We use $a^+$ (resp. $a^-$) to denote the smallest (resp. largest) $\F$-number strictly greater (resp. smaller) than the $\F$-number $a$. To capture outward rounding, we use $\ceil{r}$ (resp. $\floor{r}$) to return the smallest (resp. largest) $\F$-number greater (resp. smaller) or equal to the real number $r$. We also use $\vec{I}$ to denote a box $\langle I_1, \ldots, I_n \rangle$ and $\vec{r}$ to denote a tuple $\langle r_1, \ldots, r_n \rangle$. Finally, we use the following notations. \begin{tabbing} 12345\= \kill \> $ \LEFT([l,u]) = l$ \\ \> $ \RIGHT([l,u]) = u$ \\ \> $\CENTER([l,u]) = \floor{(l+u)/2}$ \end{tabbing} \noindent The fundamental concept of interval arithmetic is the notion of interval extension. \begin{definition}[Interval Extension] $F: \In \fto \I$ is an interval extension of $f: \Rn \fto \R$ iff \[ \forall I_1 \ldots \forall I_n: r_1 \in I_1, \ldots, r_n \in I_n \Rightarrow f(r_1,\ldots,r_n) \in F(I_1,\ldots,I_n). \] An interval relation $C: \In \fto {\it Bool}$ is an interval extension of a relation $c: \Rn \fto {\it Bool}$ iff \[ \forall I_1 \ldots \forall I_n: r_1 \in I_1, \ldots, r_n \in I_n \Rightarrow [c(r_1,\ldots,r_n) \Rightarrow C(I_1,\ldots,I_n)]. \] \end{definition} \begin{example} The interval function $\oplus$ defined as \[ [a_1,b_1] \oplus [a_2,b_2] = [\floor{a_1 + a_2}, \ceil{b_1 + b_2}] \] is an interval extension of addition of real numbers. The interval relation $\approx$ defined as \[ I_1 \approx I_2 \; \Leftrightarrow \; (I_1 \cap I_2 \neq \emptyset) \] is an interval extension of the equality relation on real numbers. \end{example} \noindent In the following, we assume fixed interval extensions for the basic real operators $+,-,\times$ and exponentiation (for instance, the interval extension of $+$ is defined by $\oplus$) and the basic real relations $=,\geq$. In addition, we overload the real symbols and use them for their interval extensions. Finally, we denote relations by the letter $c$ possibly subscripted, interval relations by the letter $C$ possibly subscripted. Note that constraints and relations are used as synonyms in this paper. \subsection{Unions of Intervals} \label{union} It is sometimes useful to work with unions of intervals instead of with single intervals to approximate sets of real numbers, although unions of intervals are less effective in general. A typical use of unions of intervals in interval arithmetic is to take the intersection of an interval $I$ with the result of an operation of the form $I_1 - I_2/I_3$. When $0 \notin I_3$, it is easy to compute the resulting interval by using an interval extension of division. Moreover, computing $I_1 - I_2/I_3$ with unions of intervals makes it possible to obtain interesting results even when $0 \in I_3$. The computation returns either a single interval $I_r$ or two intervals $I_n$ and $I_p$ and the final result is either $I_r \cap I$ or $(I_n \cap I) \uplus (I_p \cap I)$. Note that the result of an interval division $[a,b] / [c,d]$ using unions of intervals where $c \leq 0 \leq d$ and $c < d$ is specified as follows: \begin{tabbing} 12345\=1234567890123456789012345667890\=\kill \> $[\floor{b/c},\infty]$ \> if $b \leq 0$ and $d = 0$ \\ \> $[-\infty,\ceil{b/d}] \cup [\floor{b/c},\infty]$ \> if $b \leq 0$ and $c < 0 < d$ \\ \> $[-\infty,\ceil{b/d}]$ \> if $b \leq 0$ and $c = 0$ \\ \> $[-\infty,\infty]$ \> if $a < 0 < b$ \\ \> $[-\infty,\ceil{a/c}]$ \> if $a \geq 0$ and $d = 0$ \\ \> $[-\infty,\ceil{a/c}] \cup [\floor{b/c},\infty]$ \> if $a \geq 0$ and $c < 0 < d$ \\ \> $[\floor{a/d},\infty]$ \> if $a \geq 0$ and $c = 0$. \end{tabbing} while $[a,b]/[0,0] = [-\infty,\infty]$. \subsection{Constraint Representations} It is well-known that different computer representations of a real function produce different results when evaluated with floating-point numbers on a computer. As a consequence, the way constraints are written may have an impact on the behaviour on the algorithm. For this reason, a constraint or a function in this paper is considered to be an expression written in terms of the real functional and relational symbols. This can be formalized by assuming that all constraints used in this paper are written in a language whose abstract syntax is specified by the following grammar: \\ \\ % % \begin{tabular}{lll} $c$&$ \in$&$ \Constraint$\\ $f$&$ \in$&$ \Function$\\ $x_i$&$ \in$&$ {\it RealVariables} = \{\xnseq\}$\\ $q$&$ \in$&$ {\cal Q}$\\ $n$&$ \in$&$ {\cal N}$\\ \end{tabular} \\ \\ \noindent\begin{tabular}{lll} $c$&$ ::=$&$ f = 0 \ | \ f \geq 0 $ \\ $f$&$ ::=$&$ q \ | \ x_i \ | \ f + f \ | \ f - f \ | \ f \times f \ | \ f^n \ | \ (f) $ \end{tabular} \\ \noindent It is easy to extend the language to accommodate more functions (e.g. {\it sin, cos, $\ldots$}). Note that we assume for simplicity that all constraints are written using a finite (but arbitrary large) set of variables $\{\xnseq\}$. Let $e$ be a tuple $(\rnseq)$ and $e_{\mid i}$ be $r_i$ $(1 \leq i \leq n)$. The semantics of the language is given by the following semantic functions whose signatures are \begin{tabbing} 12345\=\kill \> ${\cal S}_c: \Constraint \fto \Rn \fto {\it Bool}$ \\ \> ${\cal S}_f: \Function \fto \Rn \fto \R$ \end{tabbing} and whose semantic equations are specified as follows: \\ \\ \begin{tabular}{lll} ${\cal S}_c \sem{f = 0} $& = &$\lambda e. \ {\cal S}_f \sem{f} e \ = \ 0$ \\ ${\cal S}_c \sem{f \geq 0} $& = &$\lambda e. \ {\cal S}_f \sem{f} e \ \geq \ 0$ \\ ${\cal S}_f \sem{q} $& = &$\lambda e. \ q$\\ ${\cal S}_f \sem{x_i}$& = &$\lambda e. \ e_{\mid i}$\\ ${\cal S}_f \sem{f_1 + f_2} $& = &$\lambda e. \ {\cal S}_f \sem{f_1} e + {\cal S}_f \sem{f_2} e $\\ ${\cal S}_f \sem{f_1 - f_2} $& = &$\lambda e. \ {\cal S}_f \sem{f_1} e - {\cal S}_f \sem{f_2} e $\\ ${\cal S}_f \sem{f_1 \times f_2} $& = &$\lambda e. \ {\cal S}_f \sem{f_1} e \times {\cal S}_f \sem{f_2} e $\\ ${\cal S}_f \sem{f^n} $& = &$\lambda e. \ ({\cal S}_f \sem{f} e)^n$\\ ${\cal S}_f \sem{(f)} $& = &$\lambda e. \ {\cal S}_f \sem{f} e$ \end{tabular} \\ \noindent In the following, we often abuse notation and use $f$ (resp. $c$) to denote ${\cal S}_f \sem{f}$ (resp. ${\cal S}_c \sem{c}$ ) and vice-versa. However, the meaning should be clear from the context. We also assume that interval functions and interval relations are defined using a similar abstract syntax and semantics except that real variables (i.e. $x_i$) are replaced by interval variables (i.e. $X_i$) and rational numbers by intervals. We use $\CONSTRAINT$ and $\FUNCTION$ to denote the set of interval constraints and functions. For simplicity of exposion, we restrict attention to equations. It is straightforward to generalize our results to inequalities (see Section \ref{tricks}). \section{The Branch \& Prune Algorithm} \label{pruning} This section describes the branch \& prune algorithm {\tt Newton}. Section \ref{box} defines box-consistency, the key concept behing our algorithm. Section \ref{extensions} shows how box-consistency can be instantiated to produce various pruning operators achieving various tradeoffs between accuracy and efficiency. Section \ref{cond} defines a conditioning operator used in {\tt Newton} to improve the effectiveness of box-consistency. Section \ref{newton-consistency} specifies the pruning in {\tt Newton}. Section \ref{algo} describes the algorithm. Recall that we assume that all constraints are defined over variables $\xnseq$. \subsection{Box Consistency} \label{box} Box-consistency \cite{ILPS-94-CO} is an approximation of arc-consistency, a notion well-known in artificial intelligence \cite{Mackworth77} which states a simple local condition on a constraint $c$ and the set of possible values for each of its variables, say $D_1,\ldots,D_n$. Informally speaking, a constraint $c$ is arc-consistency if none of the $D_i$ can be reduced by using projections of $c$. \begin{definition}[Projection Constraint] A projection constraint $\langle c, i \rangle$ is the association of a constraint $c$ and of an index $i$ ($1 \leq i \leq n$). Projection constraints are denoted by the letter $P$, possibly subscripted. \end{definition} \begin{definition}[Arc-Consistency] A projection constraint $\langle c, i\rangle$ is arc-consistent wrt $\langle D_1, \ldots,D_n \rangle$ iff $ D_i = D_i \cap \{ r_i \ \mid \ \exists r_1 \in D_1, \ldots, \exists r_{i-1} \in D_{i-1}, \ldots, \exists r_{i+1} \in D_{i+1}, \ldots, \exists r_{n} \in D_n: c(r_1,\ldots,r_n) \}. $ A constraint $c$ is arc-consistent wrt $\langle D_1,\ldots,D_n \rangle$ if each of its projections is arc-consistent wrt $\langle D_1,\ldots,D_n \rangle$. A system of constraints ${\cal S}$ is arc-consistent wrt $\langle D_1,\ldots,D_n \rangle$ if each constraint in ${\cal S}$ is arc-consistent wrt $\langle D_1,\ldots,D_n \rangle$. \end{definition} \noindent Arc-consistency cannot be computed in general when working with real numbers and polynomial constraints and simple approximations to capture machine precision are very expensive to compute. For instance, a simple approximation of arc-consistency consists in working with intervals and approximating the set computed by arc-consistency to return an interval, i.e. \[ I_i = \bo{I_i \cap \{ \ r_i \ \mid \ \exists r_1 \in I_1, \ldots, \exists r_{i-1} \in I_{i-1}, \ldots, \exists r_{i+1} \in I_{i+1}, \ldots, \exists r_{n} \in I_n: c(r_1,\ldots,r_n) \ \}}. \] This condition, used in systems like \cite{BNR,CLP-BNR}, is easily enforced on simple constraints such as \[ x_1 = x_2 + x_3, \;\; x_1 = x_2 - x_3, \;\; x_1 = x_2 \times x_3 \] but it is also computationally very expensive for complex constraints with multiple occurrences of the same variables. Moreover, decomposing complex constraints into simple constraints entails a substantial loss in pruning, making this approach unpractical on many applications. The notion of box-consistency introduced in \cite{ILPS-94-CO} is a coarser approximation of arc-consistency which provides a much better trade-off between efficiency and pruning. It consists in replacing the existential quantification in the above condition by the evaluation of an interval extension of the constraint on the intervals of the existential variables. Since different interval extensions will produce different results, it is necessary and useful to parametrize the definition of box-consistency with an operator which maps a constraint into one of its extensions. \begin{definition}[Extension Operator] An extension operator is a function of signature \[ \Constraint \fto \CONSTRAINT \] \end{definition} \begin{definition}[Box-Consistency] An projection constraint $\langle c, i \rangle$ is box-consistent wrt $\langle \Inseq \rangle$ and an extension operator \IE{} iff \[ \IE(c)(I_1,\ldots,I_{i-1},[l,l^+],I_{i+1},\ldots,I_n) \; \wedge \; \IE(c)(I_1,\ldots,I_{i-1},[u^-,u],I_{i+1},\ldots,I_n). \] where $l = \LEFT(I_i)$ and $u = \RIGHT(I_i)$. A constraint is box-consistent wrt $\vec{I}$ and \IE{} if each of its projections is box-consistent wrt $\vec{I}$ and \IE{}. A system of constraints is box-consistent wrt $\vec{I}$ and \IE{} iff each constraint in the system is box-consistent wrt $\vec{I}$ and \IE{}. \end{definition} \subsection{Interval Extensions for Box Consistency} \label{extensions} Box-consistency strongly depends on the extension operator $\IE$ and different choices of operators can produce very different (often incomparable) tradeoffs between pruning and computational complexity. In this section, we consider the three extension operators used in \newton{}: natural interval extension, expanded interval extension, and Taylor interval extension. \subsubsection{Natural Interval Extension} The easiest extension operator produces the natural interval extension of a constraint. Informally speaking, it consists in replacing each real number by an interval, each real variable by an interval variable, each real operation by its fixed interval extension. Formally, it is captured by the following definition. \begin{definition}[Natural Interval Extension] The natural interval extension of a constraint $c$, denoted by $\NE(c)$ and of signature $\NE: \Constraint \fto \CONSTRAINT$ is defined inductively as follows: \\ \\ \begin{tabular}{lll} $\NE \sem{f = 0} $& = &$\NE \sem{f} = \ap{0}$ \\ $\NE \sem{q} $& = &$\ap{q}$\\ $\NE \sem{x_i} $& = &$X_i$\\ $\NE \sem{f_1 + f_2} $& = &$\NE \sem{f_1} + \NE \sem{f_2} $\\ $\NE \sem{f_1 - f_2} $& = &$\NE \sem{f_1} - \NE \sem{f_2} $\\ $\NE \sem{f_1 \times f_2} $& = &$\NE \sem{f_1} \times \NE \sem{f_2} $\\ $\NE \sem{f^n} $& = &$(\NE \sem{f} )^n$\\ $\NE \sem{(f)} $& = &$\NE \sem{f} $ \end{tabular} \\ \end{definition} The advantage of this extension is that it preserves the way constraints are written and hence users of the system can choose constraint representations particularly appropriate for the problem at hand. A very nice application where this extension is fundamental is the Mor\'{e}-Cosnard discretization of a nonlinear integral equation (See Section \ref{more-sec}). \subsubsection{Expanded Interval Extension} The second interval extension used by \newton{} does not preserve the way constraints are written but uses an expanded form of the constraints. The key advantage of this extension is that it allows the algorithm to enforce box-consistency by applying Newton interval method on univariate real functions. The real functions are derived from univariate interval constraints obtained by replacing all but one variable by their intervals. As a consequence, applying box-consistency will be particularly efficient, although the pruning may be weaker than for the natural extension due to the dependency problem of interval arithmetics.\footnote{Note that it is not always necessary to go through the expanded form to obtain the above property but \newton{} adopts it for simplicity.} The main use of this extension is to speed up box-consistency on the natural extension but it may happen that it gives more precision as well. \begin{definition}[Expanded Form] A constraint $c$ is in expanded form iff it is of the form \[ p_1 + \ldots + p_n \; = \; 0 \] where each $p_i$ is of the form $q x_1^{m_1} \ldots x_n^{m_n}$ with $q \in {\cal Q}$ and $m_i \in {\cal N}$ and any pair of products $(p_i,p_j)$ ($i \neq j$) is such that $p_i$ and $p_j$ have different sequences of exponents $m_1,\ldots,m_n$. Given a constraint $c$, we denote its expanded version by $\EF(c)$.\footnote{The expanded version can easily be turned into a canonical representation for constraints.} \end{definition} \begin{definition}[Expanded Interval Extension] The expanded interval extension of a constraint $c$, denoted by $\EE(c)$ and of signature $\EE: \Constraint \fto \CONSTRAINT$, is \[ \NE(\EF(c)). \] \end{definition} \subsubsection{Taylor Interval Extension Generator} The last interval extension we introduce is based on the taylor expansion around a point. Instead of presenting a single extension, we present a generator of extensions, since it will be important to generate an extension for a box under consideration. The generator has signature \[ \In \fto \Constraint \fto \CONSTRAINT \] and it assumes that the constraint which it is applied to is of the form $f \ = \ 0$ where $f$ denotes a function which has continuous partial derivatives. Given these assumptions, the key idea behind the operator is to apply a Taylor expansion of the function around the center of the box and to bound the rest of the series using the box. \begin{definition}[Taylor Interval Extension Generator] Let $c$ be a constraint $f = 0$, $f$ be a function with continuous partial derivatives, $\vec{I}$ be a box $(\Inseq)$, and $c_i$ be the center of $I_i$. The Taylor interval extension generator $\TE(\vec{I})(c)(X_1,\ldots,X_n)$ is defined \[ \NE(f)(\ap{c_1},\ldots,\ap{c_n}) + \sum_{i=1}^{n} \DE(f,i)(I_1,\ldots,I_n) \; (X_i - \ap{c_i}) \;\; = \;\; \ap{0}. \] where $\DE(f,i)$ returns an interval extension of $\frac{\partial f}{\partial x_i}$. \end{definition} \noindent In the current version of our system, the operator $\DE$ of signature $(\Function \times {\cal N}) \fto \FUNCTION$ is defined as $\NE(D(f,i))$ where $D: (\Function \times {\cal N}) \fto \FUNCTION$ is specified as follows: \\ \\ \begin{tabular}{lll} $D \sem{q,i} $& = &$0$\\ $D \sem{x_i,i} $& = &$1$\\ $D \sem{x_j,i} $& = &$0$ $\;\;$ $(i \neq j)$\\ $D \sem{f_1 + f_2,i} $& = &$D \sem{f_1,i} + D \sem{f_2,i} $\\ $D \sem{f_1 - f_2,i} $& = &$D \sem{f_1,i} - D \sem{f_2,i} $\\ $D \sem{f_1 \times f_2,i} $& = &$f_1 \times D \sem{f_2,i} + f_2 \times D \sem{f_1,i} $\\ $D \sem{f^n,i} $& = &$n \times f^{n-1} \times D \sem{f,i}$\\ $D \sem{(f),i} $& = &$D \sem{f,i} $\\ \end{tabular} \\ \subsection{Conditioning} \label{cond} It is interesting to note that box-consistency on the Taylor interval extension is closely related to Hansen-Segupta's operator \cite{Hansen81}, which is an improvement over Krawczyk's operator \cite{Krawczyk69}.\footnote{Note that the operators of Krawczyk and Hansen-Segupta can be used to prove the existence of a solution. We will come back to this point in the implementation section.} Hansen and Smith \cite{Hansen67} also argued that these operators are more effective for a system $\{ f_1 = 0, \ldots, f_n = 0 \}$ wrt a box $\langle \Inseq \rangle$ when the interval matrix \[ M_{ij} = \DE(f_i,j)(I_1,\ldots,I_n) \; (1 \leq i, j \leq n) \] is diagonally dominant, i.e., \[ mig(M_{i,j}) \geq \sum_{j=1, j\neq i}^{n} mag(M_{i,j}) \] where \[ {\it mig} ([l,u]) = {\it min}(\mid l \mid, \mid u \mid) \mbox{ and } {\it mag} ([l,u]) = {\it max}(\mid l \mid, \mid u \mid). \] They also suggest a conditioning which consists in multiplying the linear relaxation by a real matrix which is the inverse of the matrix obtained by taking the center of $M[i,j]$. The resulting system is generally solved through Gauss-Seidel iterations, giving Hansen-Segupta's operator. {\tt Newton} exploits this idea to improve the effectiveness of box-consistency on the Taylor interval extension. The conditioning of {\tt Newton} is abstracted by the following definition. \begin{definition}[Conditioning] Let $S = \{ f_1 = 0, \ldots, f_n = 0\}$. A conditioning of $S$ is a system $S' = \{ f_1' = 0, \ldots, f_n' = 0\}$ where \[ f_i' = \sum_{k=1}^{n} a_{ik} f_k \] where $a_{ik} \in {\cal Q}$. \end{definition} \noindent In its present implementation, {\tt Newton} uses a conditioning {\it cond}($\{ f_1 = 0, \ldots, f_n = 0\}$,$\vec{I}$) which returns a system $\{ f_1' = 0, \ldots, f_n' = 0\}$ such that \[ f_i' = \sum_{k=1}^{n} A_{ik} f_k \] where \[ \begin{array}{lll} M_{ij} & = & \DE(f_i,j)(I_1,\ldots,I_n) \; (1 \leq i, j \leq n) \\ B_{ij} & = & \CENTER(M_{ij}) \\ A_{ij} & = & \left\{ \begin{array}{ll} B^{-1}_{ij} & \mbox{ if $B$ is not singular } \\ {\it if} \; i = j \; {\it then} \; 1 \; {\it else} \; 0 & \mbox{ otherwise.} \end{array} \right. \end{array} \] Note that the computation of the inverse of $B$ is obtained by standard floating-point algorithms. \subsection{Pruning in {\tt Newton}} \label{newton-consistency} We are now in position to specify the pruning of {\tt Newton}. \begin{definition}[Newton Consistency] A system of constraints ${\cal S}$ is Newton-consistent wrt $\vec{I}$ if $\CS$ is box-consistent wrt $\vec{I}$ and $\NE$, $\CS$ is box-consistent wrt $\vec{I}$ and $\EE$, and $\CON(\CS,\vec{I})$ is box-consistent wrt $\vec{I}$ and $\TE(\vec{I})$. \end{definition} \noindent The key idea behind {\tt Newton} is to apply Newton-consistency at each node of the search tree, i.e., {\tt Newton} reduces the current of values for the variables in such a way that the constraint system is Newton-consistent wrt the reduced intervals and no solution is removed. This is captured by the following definition. \begin{definition}[Maximal Newton-Consistency] $\vec{I}$ is maximally Newton-consistent wrt a system of constraints ${\cal S}$ and an initial box $\vec{I_0}$ iff $\vec{I}$ is the largest box included in $\vec{I_0}$ such that ${\cal S}$ is Newton-consistent wrt $\vec{I}$. In the following, we denote the box $\vec{I}$ by {\tt PRUNE}($\CS,\vec{I_0}$) \end{definition} \noindent The following soundness result follows from the definition of interval extensions and the soundness of arc-consistency. \begin{proposition}[Soundness of Newton-Consistency] Let $\vec{I}$ be maximally Newton-consistent wrt a system of constraints ${\cal S}$ and an initial box $\vec{I_0}$ and let $\vec{r} \in \vec{I_0}$. If $\vec{r}$ is a solution of ${\cal S}$, then $\vec{r} \in \vec{I}$. \end{proposition} \subsection{The Branch and Prune Algorithm {\tt Newton}} \label{algo} We are now ready to present the branch and prune algorithm which is depicted in Figure \ref{mainalgo}. The algorithm makes use of a a function {\tt BRANCH} is used to split a box along one dimension to produce two subboxes. Note a box is empty as soon as one of its components is empty. The key idea behind the algorithm is extremely simple. It applies operation {\tt PRUNE} on the initial box. If the resulting box is empty, there is no solution. If it is small enough, the box is returned as result. Otherwise, the algorithm splits the box into two parts and applies the algorithm recursively. Our experimental results assume a simple round-robin heuristics to choose the variable to be split in operation {\tt BRANCH}. \begin{figure} \begin{tt} \begin{tabbing} 123\=123\=123\=123\=123\=123\=\kill {\bf procedure} BranchAndPrune(${\cal S}$: Set of Constraint; $\vec{I_0}: \In$): Set of $\In$; \\ {\bf begin} \\ \> $\vec{I}$ := PRUNE(${\cal S}$,$\vec{I_0}$); \\ \> {\bf if} $\neg$ IsEmpty($\vec{I}$) {\bf then} \\ \>\> {\bf if} IsSmallEnough($\vec{I}$) {\bf then} \\ \>\>\> {\bf return} $\{\vec{I}\}$ \\ \>\> {\bf else} \\ \>\>\> $\langle \vec{I_1}, \vec{I_2} \rangle$ := BRANCH($\vec{I}$);\\ \>\>\> {\bf return} BranchAndPrune(${\cal S}$,$\vec{I_1}$) $\cup$ BranchAndPrune(${\cal S}$,$\vec{I_2})$ \\ \>\> {\bf endif} \\ \> {\bf else} \\ \>\> {\bf return} $\emptyset$ \\ \> {\bf endif} \\ {\bf end} \end{tabbing} \end{tt} \caption{The Branch and Prune Algorithm} \label{mainalgo} \end{figure} \section{Implementation} \label{implementation} The purpose of this section is to show how to implement Newton-consistency for the various interval extensions. The key idea behind the implementation is to define narrowing operators for each projection. Given these operators, it is easy to enforce maximal Newton-consistency by applying them until no reduction occurs. We first describe a simple fixpoint algorithm using narrowing operators. We then show how box-consistency on the interval extensions can be implemented. Implementation issues are also discussed. \subsection{The Fixpoint Algorithm for Newton-Consistency} \begin{figure} \begin{tt} \begin{tabbing} 12345\=12345\=12345\=\kill {\bf procedure} PRUNE({\bf in} ${\cal S}$: Set of Constraint; {\bf inout} $\vec{I}: \In$) \\ {\bf begin} \\ \> {\bf repeat} \\ \>\> $\vec{I_p}$ := $\vec{I}$; \\ \>\> BOX-PRUNE($\CS$,$\CON(\CS,\vec{I})$,$\vec{I}$); \\ \> {\bf until} $\vec{I} = \vec{I_p}$; \\ {\bf end} \\ \\ {\bf procedure} BOX-PRUNE(${\cal S},{\cal S'}$: Set of Constraint; {\bf inout} $\vec{I}: \In$) \\ {\bf begin} \\ \> {\bf repeat} \\ \>\> $\vec{I_p}$ := $\vec{I}$; \\ \>\> $\vec{I}$ := $\bigcap \{ \mbox{{\tt BOX-NARROW}}((c,i),\vec{I},\EE) \ \mid \ c\in {\cal S} \; \& \; 1 \leq i \leq n \}$; \\ \>\> $\vec{I}$ := $\bigcap \{ \mbox{{\tt BOX-NARROW}}((c,i),\vec{I},\NE) \ \mid \ c\in {\cal S} \; \& \; 1 \leq i \leq n \}$; \\ \>\> $\vec{I}$ := $\bigcap \{ \mbox{{\tt BOX-NARROW}}((c,i),\vec{I},\TE(\vec{I})) \ \mid \ c\in {\cal S'} \; \& \; 1 \leq i \leq n \}$; \\ \> {\bf until} $\vec{I} = \vec{I_p}$; \\ {\bf end} \end{tabbing} \end{tt} \caption{A Fixpoint Algorithm for Operation {\tt PRUNE}} \label{fixpoint} \end{figure} \noindent Newton-consistency can be enforced easily by using narrowing operators deduced from the definition of box-consistency. These operators are used to reduce the interval of a variable using a projection constraint. \begin{definition}[Box-Narrowing] Let $\langle c,i \rangle$ be a projection constraint, $\langle \Inseq \rangle$ be a box, and $\IE$ be a interval extension operator. The narrowing operator {\tt BOX-NARROW} is defined as \begin{tt} \begin{tabbing} BOX-NARROW($\langle c, i\rangle$,$\langle \Inseq \rangle$,\IE) = $\langle I_1,\ldots,I_{i-1},I,I_{i+1},\ldots,I_n \rangle$ \end{tabbing} \end{tt} where $I$ is defined as the largest set included in $I_i$ such that $\langle c,i\rangle$ is box-consistent with respect to $\langle I_1, \ldots, I_{i-1}, I , I_{i+1}, \ldots, I_{n} \rangle$ and $\IE$. \end{definition} Maximal Newton-consistency is essentially obtained by applying the narrowing operators of each projection until no further reduction occurs. This fixpoint algorithm is depicted in Figure \ref{fixpoint}. Termination of the algorithm is guaranteed since the set $\F$ is finite and thus the intervals can only be reduced finitely often. \subsection{Box-Consistency on the Natural Interval Extension} \begin{figure} \begin{tt} \begin{tabbing} 123\=123\=123\=123\=123\=123\=\kill {\bf function} {\tt LNAR}($F,F': \I \rightarrow \I$, $I: \I$): $\I$; \\ {\bf begin} \\ \> $I$ := $N^*(F,F',I)$; \\ \> {\bf if} 0 $\in$ $F([\LEFT(I),\LEFT(I)^+])$ {\bf then} \\ \>\> {\bf return} $I$ \\ \> {\bf else if} 0 $\notin$ $F(I)$ {\bf then} \\ \>\> {\bf return} $\emptyset$ \\ \> {\bf else if} \\ \>\> $\langle I_1, I_2 \rangle$ := SPLIT($I$); \\ \>\> {\bf if} ${\tt LNAR}(F,F',I_1) \neq \emptyset$ {\bf then} \\ \>\>\> {\bf return} ${\tt LNAR}(F,F',I_1) \cup I_2$ \\ \>\> {\bf else} \\ \>\>\> {\bf return} ${\tt LNAR}(F,F',I_2)$ \\ \> {\bf endif} \\ {\bf end;} \end{tabbing} \end{tt} \caption{The function {\tt LNAR}} \label{LNAR} \end{figure} Box-consistency can often be reduced to two subproblems which, informally speaking, consist in shrinking the left (resp. the right) of an interval $I'$ to the leftmost (resp. rightmost) zero of a univariate interval function $F$ in $I'$. Moreover, for polynomial systems, we will always have at our disposal a function $F'$ such that, for all $I \subseteq I'$, \[ 0 \in F(I) \; \Rightarrow \; 0 \in N(F,F',I) \] where \[ N(F,F',I) = I \ \cap \ \ap{\CENTER(I)} - \frac{F(I)}{F'(I)}. \] Note that the division and the top-level substraction are evaluated with unions of intervals and the intersection returns a single interval (see Section \ref{union}), The shrinkings can now be obtained easily using a variation of the interval Newton method. Figure \ref{LNAR} depicts a simple procedure {\tt LNAR} to shrink to the leftmost zero (procedure {\tt RNAR} is defined similarly). It make uses of $N^*(F,F',I) = \bigcap_{i=0}^{\infty} I_i$ where \[ \begin{array}{lll} I_0 & = & I \\ I_{i+1} & = & N(F,F',I_i) \; \; (0 \leq i) \\ \end{array} \] and of a function {\tt SPLIT} which, given an interval $I$, returns two intervals $I_1$ and $I_2$ such that $I = I_1 \cup I_2$ and $\LEFT(I_2) \geq \RIGHT(I_1)$. We are now in position to define the narrowing operator for the natural interval extension. \begin{definition}[Narrowing Operator for the Natural Interval Extension] \label{nar-nat} The narrowing operator for the natural interval extension is defined as follows: \begin{tt} \begin{tabbing} 12345\=123\=\kill BOX-NARROW($\langle f = 0, i \rangle,\langle \Inseq \rangle, \NE$) = $ \langle I_1,\ldots,I_{i-1},I,I_{i+1},\ldots,I_n \rangle$ \\ {\bf where} \\ \> $I$ \> = {\tt LNAR}$(F,F',{\tt RNAR}(F,F',I_i))$ \\ \> $F$ \> = $\lambda X. \ \NE(f)(I_1,\ldots,I_{i-1},X,I_{i+1},\ldots,I_{n})$ \\ \> $F'$ \> = $\lambda X. \ \DE(f,i)(I_1,\ldots,I_{i-1},X,I_{i+1},\ldots,I_{n})$ \end{tabbing} \end{tt} \end{definition} It is important to note that box-consistency on the natural extension (and on the expanded extension as well) can be applied even if the function is not differentiable. It suffices to omit the application of operator $N^*$ in the functions {\tt LNAR} and {\tt RNAR}. \subsection{Box-Consistency on the Expanded Interval Extension} \noindent We now turn to the narrowing operator for the expanded interval extension. Box-consistency on the expanded interval extension can be enforced by using the interval function \[ F = \lambda X. \ \EE(f)(I_1,\ldots,I_{i-1},X,I_{i+1},\ldots,I_{n}) \] and by applying similar techniques as in the previous section. However, since the function is in expanded form, it is possible to do better by using an idea from \cite{Hong95}. The key insight is to sandwitch $F$ exactly between two real functions $f_l$ and $f_u$ defined as follows: \[ \begin{array}{lll} f_l & = & \lambda x. {\it left}(F(\ap{x})) \\ f_u & = & \lambda x. {\it right}(F(\ap{x})). \end{array} \] Box-consistency can now be enforced by searching the leftmost and rightmost zeros of these two functions using, say, the interval Newton method. The key advantage of the expanded interval extension is that it is easy to define the function $f_l$ and $f_u$ constructively. Let $F$ be of the form \[ F = \lambda X. \ I_1 X^{n_1} + \ldots + I_p X^{n_p}. \] Function $f_l$ is defined as \[ f_l = \lambda x. \ {\it low}(I_1,x,n_1) + \ldots + {\it low}(I_p,x,n_p) \] where \\ \[ {\it low}(I,x,n) = \left\{ \begin{array}{ll} {\it left}(I) \ x^n & \mbox{ if } x \geq 0 \vee j \mbox{ is even } \\ {\it right}(I) \ x^n & \mbox{ otherwise }\\ \end{array} \right. \] Function $f_u$ is defined as \[ f_u = \lambda x. \ {\it high}(I_1,x,n_1) + \ldots + {\it high}(I_p,x,n_p) \] where \\ \[ {\it high}(I,x,n) = \left\{ \begin{array}{ll} {\it right}(I) \ x^n & \mbox{ if } x \geq 0 \vee j \mbox{ is even } \\ {\it left}(I) \ x^n & \mbox{ otherwise }\\ \end{array} \right. \] It is easy to see that the two definitions of $f_l$ and $f_u$ are equivalent. The narrowing operator can now be obtained by using Newton interval method on the above two functions. Some care must be applied obviously since $f_l$ and $f_u$ are not differentiable at 0. The method is more efficient than applying the interval Newton method on the interval function, since intervals have been replaced by numbers increasing the precision of the Newton operator $N^*$. \begin{definition}[Narrowing Operator for the Expanded Interval Extension] Let $\langle f = 0, i \rangle$ be a projection constraint. Let $F$ be $\lambda X. \ \EE(f)(I_1,\ldots,I_{i-1},X,I_{i+1},\ldots,I_{n})$, $f_l$ be $\lambda x. {\it left}(F(\ap{x}))$, $f_u$ be $\lambda x. {\it right}(F(\ap{x}))$. The narrowing operator for the expanded interval extension is defined as follows: \begin{tt} \begin{tabbing} 123\=123\=\kill BOX-NARROW($\langle (f = 0, \EE), i \rangle,\langle \Inseq \rangle$) = $ \langle I_1,\ldots,I_{i-1},[l_i,u_i],I_{i+1},\ldots,I_n \rangle$ \\ {\bf where} \\ \> $l_i$ \> = $\left\{ \begin{array}{ll} l & \mbox{ if } 0 \in F([l,l^+]) \\ {\it left}({\tt LNAR}(F_l,F_l',[l,0]) \cup {\tt LNAR}(F_l,F_l',[0,u]) & \mbox{ if } F([l,l^+]) \geq \ap{0} \\ {\it left}({\tt LNAR}(F_u,F_u',[l,0]) \cup {\tt LNAR}(F_u,F_u',[0,u]) & \mbox{ otherwise } \\ \end{array} \right.$ \\ \> $u_i$ \> = $\left\{ \begin{array}{ll} u & \mbox{ if } 0 \in F([u^-,u]) \\ {\it right}({\tt RNAR}(F_l,F_l',[l,0]) \cup {\tt RNAR}(F_l,F_l',[0,u]) & \mbox{ if } F([l,l^+]) \geq \ap{0} \\ {\it right}({\tt RNAR}(F_u,F_u',[l,0]) \cup {\tt RNAR}(F_u,F_u',[0,u]) & \mbox{ otherwise } \\ \end{array} \right.$ \\ \> $F_l$ \> = $\NE(f_l)$ \\ \> $F_u$ \> = $\NE(f_u)$ \\ \> $F_l'$ \> = $\DE(f_l)$ \\ \> $F_u'$ \> = $\DE(f_u)$ \\ \> $l$ \> = ${\it left}(I_i)$ \\ \> $u$ \> = ${\it right}(I_i)$ \end{tabbing} \end{tt} \end{definition} \subsection{Taylor Interval Extension} We conclude by presenting the narrowing operator for the Taylor interval extension. Box-consistency on the Taylor interval extension can be enforced by using the interval function \[ F = \lambda X. \ \TE(f,\langle \Inseq \rangle)(I_1,\ldots,I_{i-1},X,I_{i+1},\ldots,I_{n}) \] and by applying a simple implementation of {\tt LNAR} and {\tt RNAR} where the Newton operator $N^*$ is omitted. However, it is possible to do better by noticing that the constraint is of the form \begin{small} \[ \NE(f)(\ap{c_1},\ldots,\ap{c_n}) + \sum_{j=1}^{i-1} \DE(f,j)(\vec{I}) (I_j - \ap{c_j}) + \DE(f,j)(\vec{I}) (X_i - \ap{c_i}) + \sum_{j=i+1}^{n} \DE(f,j)(\vec{I}) (I_j - \ap{c_j}) = \ap{0} \] \end{small} where $c_i = \CENTER(I_i)$ and contains a single variable which can be isolating to compute box-consistency directly.\footnote{Strictly speaking, this technique can increase the precision of box-consistency by one floating point number since we remove an arithmetic operation. It could be formalized by associating interval extensions with projections. Doing this would obscure the paper without adding any real insight.} \begin{definition}[Narrowing Operator for the Taylor Interval Extension] The narrowing operator for the Taylor interval extension is defined as follows: \[ \mbox{{\tt BOX-NARROW}}(\langle f = 0, i \rangle,\langle \Inseq \rangle,\TE) = \langle I_1,\ldots,I_{i-1},I,I_{i+1},\ldots,I_n \rangle \] where \[ I = I_i \cap \ap{c_i} - \frac{1}{\DE(f,i)(I_1,\ldots,I_n)} [\sum_{j=1, j \neq i}^{n} \DE(f,j)(I_1,\ldots,I_n) \; (I_j - \ap{c_i}) - \NE(f)(\ap{c_1},\ldots,\ap{c_n})]. \] and $c_i = \CENTER(I_i)$. \end{definition} \noindent Note that the above expression is evaluated using unions of intervals (see Section \ref{union}). \subsection{Implementation Issues} \label{tricks} We now review some implementation issues which arise in programming the algorithm. \paragraph{Priorities} In the fixpoint algorithm, it is important for efficiency reasons to use a priority queue to ensure that projections over the expanded interval extension be selected first and projections over the Taylor interval extension be selected last. This is motivated by the fact that box-consistency over the expanded interval extension can be enforced efficiently and produces significant pruning when far from a solution, while box-consistency over the Taylor interval extension requires the computation of the Jacobian of the systems of constraints and a preconditioning of the system and is mainly effective when near a solution. {\tt Newton} also does not enforce box-consistency on the natural extension when the constraints are given in expanded form and does not enforce box-consistency on the expanded version whenever it is believed that the expansion would lose too much precision (e.g. an expression raised to some power). \paragraph{Precision} In practice, it is often sufficient to return intervals whose widths\footnote{The width of $[l,u]$ is $u-l$.} are below at certain level instead of returning intervals of the form $[l,l^+]$. It is easy to modify the {\tt BRANCH} operation to split only intervals whose widths is above the required accuracy. Our system allows users to specify the accuracy. \paragraph{Improvement Factor} Box-consistency can sometimes take much time to remove small parts of the intervals. In these cases, it is probably more cost-effective to branch. Once again, it is easy to modify the algorithm to avoid this problem by making sure that the narrowing operators do not update the intervals unless some significant reduction has taken place. Since the notion of significant reduction may be problem-dependent, our system lets users specify the improvement factor necessary to update an interval in a projection. \paragraph{Automatic Differentiation} As mentioned, our algorithm takes a very simple approach to obtain partial derivatives, i.e., no effort is spent in factoring common expressions to reduce the dependency problem of interval arithmetic. The main reason comes from the fact that we are using automatic differentiation \cite{Rall81} to evaluate the derivatives together with the functions. This choice may be reconsidered in a further version of the system. \paragraph{Existence Proof of a Solution} As we mentioned earlier, box-consistency on a system of equations can be viewed as a generalization of Hansen-Segupta operator used in interval Newton methods. This operator can be used to prove the existence of solutions. Our algorithm is organized to exploit this property so that it is able to report this information to users. See \cite{Krawczyk86} for more discussion on this topic. \paragraph{Inequalities} It is simple to generalize the above algorithms for inequalities. In general, it suffices to test if the inequality is satisfied at the end points of the interval. If it is not, then the problem reduces once again to finding the leftmost and/or rightmost zeros. \section{Experimental Results} \label{experimental} This section reports experimental results of {\tt Newton} on a variety of standard benchmarks. The benchmarks were taken from papers on numerical analysis \cite{More79}, interval analysis \cite{Hansen81,Hong95,Moore77a}, and continuation methods \cite{Cools94,Morgan87,Morgan-kin,Morgan90}. We also compare {\tt Newton} with a traditional interval method using Hansen-Segupta's operator, range testing, and branching. This method uses the same implementation technology as {\tt Newton} and is denoted by {\tt HRB} in the following.\footnote{Some interval methods such as \cite{Hansen83} are more sophisticated than {\tt HRB} but the sophistication aims at speeding up the computation near a solution. Our main contribution is completely orthogonal and aims at speeding up the computation when far from a solution and hence comparing it to {\tt HRB} is meaningful.} Finally, we compare {\tt Newton} with a state-of-the-art continuation method \cite{Cools94}, denoted by {\tt CONT} in the following. Note that all results given in this section were obtained by running {\tt Newton} on a {\tt Sun Sparc 10} workstation to obtain all solutions. In addition, the final intervals must have widths smaller than $10^{-8}$ and {\tt Newton} always uses an improvement factor of 10\%. The results are summarized in Table \ref{summary}. For each benchmark, we give the number of variables ($n$), the total degree of the system ($d$), the initial range for the variables, and the results of each method in seconds. Note that the times for the continuation method are on a {\tt DEC 5000/200}. A space in a column means that the result is not available for the method. A question mark means that the method does not terminate in a reasonable time ($>$ 1 hour). The rest of the section describes each benchmark and the results in much more detail. \begin{table} \small \centering \begin{tabular}{||l|c|c|c||r|r|r||} \hline \hline Benchmarks & $v$ & $d$ & range & {\tt Newton} & {\tt HRB} & {\tt CONT} \\ \hline Broyden & 10 & $3^{10}$ & [-1,1] & 1.65 & 18.23 & \\ Broyden & 20 & $3^{20}$ & [-1,1] & 4.25 & ? & \\ Broyden & 320 & $3^{320}$ & [-1,1] & 113.71 & ? & \\ Broyden & 320 & $3^{320}$ & $[-10^8,10^8]$ & 143.40 & ? & \\ Mor\'{e}-Cosnard & 20 & $3^{20}$ & $[-4,5]$ & 24.49 & 968.25 & \\ Mor\'{e}-Cosnard & 40 & $3^{40}$ & $[-4,5]$ & 192.81 & ? & \\ Mor\'{e}-Cosnard & 80 & $3^{80}$ & $[-4,5]$ & 1752.64 & ? & \\ Mor\'{e}-Cosnard & 80 & $3^{80}$ & $[-10^8,0]$ & 1735.09 & ? & \\ i1 & 10 & $3^{10}$ & [-2,2] & 0.06 & 14.28 & \\ i2 & 20 & $3^{20}$ & [-1,2] & 0.30 & 1821.23 & \\ i3 & 20 & $3^{20}$ & [-2,2] & 0.31 & 5640.80 & \\ i4 & 10 & $6^{10}$ & [-1,1] & 73.94 & 445.28 & \\ i5 & 10 & $11^{10}$ & [-1,1] & 0.08 & 33.58 & \\ kin1 & 12 & 4608 & $[-10^8,10^8]$ & 14.24 & 1630.08 & \\ kin2 & 8 & 256 & $[-10^8,10^8]$ & 353.06 & 4730.34 & 35.61 \\ eco & 4 & 18 & $[-10^8,10^8]$ & 0.60 & & 1.13 \\ eco & 5 & 54 & $[-10^8,10^8]$ & 3.35 & & 5.87 \\ eco & 6 & 162 & $[-10^8,10^8]$ & 22.53 & & 50.18 \\ eco & 7 & 486 & $[-10^8,10^8]$ & 127.65 & & 991.45 \\ eco & 8 & 1458 & $[-10^8,10^8]$ & 915.24 & & \\ eco & 9 & 4374 & $[-10^8,10^8]$ & 8600.28 & & \\ combustion & 10 & 96 & $[-10^8,10^8]$ & 9.94 & & 57.40 \\ chemistry & 5 & 108 & $[0,10^8]$ & 6.32 & & 56.55 \\ neuro & 6 & 1024 & $[-10,10]$ & 0.91 & & 5.02 \\ neuro & 6 & 1024 & $[-1000,1000]$ & 172.71 & & 5.02 \\ \hline \end{tabular} \caption{Summary of the Experimental Results} \label{summary} \end{table} \subsection{Broyden Banded Functions} This is a traditional benchmark of interval techniques and was used for instance in \cite{Hansen83}. It consists in finding the zeros of the functions \[ {\textstyle f_i(x_1,\ldots,x_n) = x_i (2 + 5 x_i^2) + 1 - \sum_{j \in J_i} x_j (1 + x_j) \;\;\; (1 \leq i \leq n)} \] where $J_i = \{ j \ \mid \ j \neq i \ \& \ {\it max}(1,i-5) \leq j \leq {\it min}(n,i+1)\}$. One of the interesting features of this benchmark is that it is easy to scale up to an arbitrary dimension and hence provides a good basis to compare various methods. Table \ref{broyden-newton} reports the results of our algorithm for various sizes assuming initial intervals $[-1,1]$. For each size, we report the CPU times in seconds, the growth of the CPU time, the number of branch operations {\it branching}, the number of narrowings on the various extensions {\it na-ne, na-ee, na-te}, the total number of narrowings {\it na-tot}, the number of function evaluations (including evaluation of derivatives) for each of the extensions {\it fe-ne, fe-ee, fe-te} and the total number of function evaluations {\it fe-tot}. We also indicate %the growth of the function evaluations {\it fe-grow} and the number of preconditionings by {\it pr-con} and whether the algorithm can prove the existence of the solutions in the resulting intervals by {\it proof}. \begin{table} \small \centering \begin{tabular}{||l||r|r|r|r|r|r|r||} \hline \hline & 5 & 10 & 20 & 40 & 80 & 160 & 320 \\ \hline {\it CPU time} & 0.20 & 1.65 & 4.25 & 9.79 & 22.13 & 48.30 & 113.71 \\ {\it growth} & & 8.25 & 2.57 & 2.30 & 2.26 & 2.18 & 2.35 \\ {\it branching} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline {\it na-ne} & 57 & 260 & 661 & 1607 & 4351 & 8096 & 17126 \\ {\it na-ee} & 1226 & 8334 & 21236 & 48540 & 102797& 206926 & 414798 \\ {\it na-te} & 35 & 110 & 260 & 560 & 1200 & 2400 & 4480 \\ {\it na-tot} & 1318 & 8704 & 22157 & 50707 & 108348& 217422 & 436404 \\ \hline {\it fe-ne} & 81 & 1828 & 2943 & 5103 & 11993 & 20431 & 42125 \\ {\it fe-ee} & 3462 & 21518& 53722 & 121984& 257398& 517640 & 1036210 \\ {\it fe-te} & 95 & 320 & 920 & 2720 & 8800 & 30400 & 111360 \\ {\it fe-tot} & 3638 & 23666& 57585 & 129807& 278191& 568471 & 1189695 \\ \hline {\it pr-con} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline {\it proof} & yes & yes & yes & yes & yes & yes & yes \\ \hline \end{tabular} \caption{{\tt Newton} on the Broyden Banded functions with initial intervals $[-1,1]$} \label{broyden-newton} \end{table} \begin{table} \small \centering \begin{tabular}{||l||r|r|r|r|r|r|r||} \hline \hline & 5 & 10 & 20 & 40 & 80 & 160 & 320 \\ \hline {\it CPU time} & 0.31 & 2.15 & 5.09 & 12.49 & 27.69 & 61.60 & 143.40 \\ {\it growth} & & 6.93 & 2.36 & 2.45 & 2.21 & 2.22 & 2.32 \\ {\it branching} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline {\it pr-con} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline {\it proof} & yes & yes & yes & yes & yes & yes & yes \\ \hline \end{tabular} \caption{{\tt Newton} on the Broyden Banded functions with initial intervals $[-10^8,10^8]$} \label{broyden-newton-W} \end{table} The results indicate that {\tt Newton} solves the problem using only constraint propagation: no branching is needed. In addition, the growth of the computation times is very low and indicates that {\tt Newton} is essentially linear and can thus solve very large instances of this problem. Finally, {\tt Newton} proves the existence of a solution in the final intervals. To our knowledge, no other algorithm has all these functionalities. Table \ref{broyden-newton-W} shows the same results when the initial intervals are $[-10^8,10^8]$. They indicate that the CPU time increases only slightly in this problem when the initial intervals become substantially larger. It is interesting to note that substantial pruning is obtained by box-consistency on the natural and expanded extensions alone. For $n=10$, maximal box-consistency on these two extensions produces the intervals \begin{quote} $[-0.4283028737061274627,-0.4283028534683728794]$ \\ $[-0.4765964317901201786,-0.4765964169224605195]$ \\ $[-0.5196524683730758821,-0.5196524589206473754]$ \\ $[-0.5580993358758108425,-0.5580993137885511545]$ \\ $[-0.5925061654931400579,-0.5925061481657747375]$ \\ $[-0.6245036923913307448,-0.6245036720076052594]$ \\ $[-0.6232394806883442274,-0.6232394621928379896]$ \\ $[-0.6213938520278742273,-0.6213938315652728361]$ \\ $[-0.6204536054436834425,-0.6204535878744913413]$ \\ $[-0.5864692773020701023,-0.5864692641387999616]$ \end{quote} which have widths lower than $10^{-6}$. Note that the Hansen-Segupta's operator alone does not produce any pruning initially and returns the initial intervals whether they be of the form $[-10^8,+10^8]$ or $[-1,1]$. This indicates that box-consistency on the natural and expanded interval extensions are particularly effective when far from a solution while box-consistency on the Taylor extension (and the Hansen-Segupta's operator) is effective when near a solution. It is also interesting to stress the importance of box-consistency on the natural extension in this example to reduce the growth factor. Without it, the algorithm takes about 48 and 440 seconds instead of 27 and 61 for {\tt Newton} for $n = 80$ and $n = 160$, since the expanded interval extension loses precision due to the dependency problem. Finally, it is interesting to compare {\tt Newton} with traditional interval Newton methods. {\tt HRB} takes 0.34 seconds on $n=5$ with 18 branchings, about 18 seconds for $n=10$ with about 300 branchings, and does not return after more than an hour on $n=20$. \subsection{Discretization of a Nonlinear Integral Equation} \label{more-sec} This example comes from \cite{More79} and is also a standard benchmark for nonlinear equation solving . It consists in finding the root of the functions ${\textstyle f_k(x_1,\ldots,x_m)}$ $( 1 \leq k \leq m)$ defined as \[ x_k + \frac{1}{2(m+1)} [(1 - t_k) \sum_{j=1}^{k} t_j (x_j + t_j + 1)^3 + t_k \sum_{j=k+1}^{m} (1 -t_j)(x_j + t_j + 1)^3 ] \] where $t_j = jh$ and $h = 1/(m+1)$. These functions come from the discretization of a nonlinear integral equation, giving a constraint system denser than the sparse constraint system for the Broyden banded functions. The variables $x_i$ were given initial domains $[-4,5]$ as in \cite{Rokne88} and the computation results are given in Table \ref{more-newton}. \begin{table} \small \centering \begin{tabular}{||l||r|r|r|r|r||} \hline \hline & 5 & 10 & 20 & 40 & 80 \\ \hline {\it CPU time} & 0.75 & 4.07 & 24.49 & 192.81 & 1752.64 \\ {\it growth} & & 5.42 & 6.02 & 7.87 & 9.08 \\ {\it branching} & 0 & 0 & 0 & 0 & 0 \\ \hline {\it na-ne} & 3663 & 12616& 46555 & 213949 & 1236532 \\ {\it na-te} & 104 & 255 & 505 & 1005 & 2807 \\ {\it na-tot} & 3767 & 12871& 47060 & 214954 & 1239339 \\ \hline {\it fe-ne} & 8775 & 31837& 107977& 466595 & 2586907 \\ {\it fe-te} & 884 & 3111 & 11211 & 42411 & 166415 \\ {\it fe-tot} & 9659 & 34948& 119188& 509006 & 2753322 \\ {\it fe-grow} & & 3.62 & 3.41 & 4.27 & 5.40 \\ \hline {\it pr-con} & 1 & 1 & 1 & 1 & 1 \\ \hline {\it proof} & yes & no & no & no & no \\ \hline \end{tabular} \caption{{\tt Newton} on the Mor\'{e}-Cosnard nonlinear integral Equation with initial intervals in $[-4,5]$} \label{more-newton} \end{table} \begin{table} \small \centering \begin{tabular}{||l||r|r|r|r|r||} \hline \hline & 5 & 10 & 20 & 40 & 80 \\ \hline {\it CPU time} & 0.70 & 3.82 & 20.81 & 189.94 & 1735.09 \\ {\it growth} & & 5.45 & 5.44 & 9.12 & 9.13 \\ {\it branching} & 0 & 0 & 0 & 0 & 0 \\ \hline {\it pr-con} & 0 & 0 & 0 & 0 & 0 \\ \hline {\it proof} & yes & no & no & no & no \\ \hline \hline \end{tabular} \caption{{\tt Newton} on the Mor\'{e}-Cosnard nonlinear integral Equation with initial intervals in $[-10^8,0]$} \label{more-newton-W} \end{table} \begin{table} \small \centering \begin{tabular}{||l||r|r|r|r|r||} \hline \hline & 5 & 10 & 20 & 40 & 80 \\ \hline {\it CPU time} & 0.66 & 7.76 & 968.25 & ? & ? \\ {\it branching} & 5 & 24 & 508 & ? & ? \\ \hline {\it fe-tot} & 3709 & 20194 & 1285764 & ? & ? \\ \hline {\it pr-con} & 7 & 32 & 667 & ? & ? \\ \hline {\it proof} & yes & no & no & ? & ? \\ \hline \hline \end{tabular} \caption{The {\tt HRB} algorithm on the Mor\'{e}-Cosnard nonlinear integral Equation with initial intervals in $[-4,5]$} \label{more-hansen} \end{table} Once again, it is interesting to note that {\tt Newton} is completely deterministic on this problem, i.e. it does not do any branching. {\tt Newton} is probably cubic in the number of variables for this problem. It is important to point out the critical role of box-consistency on the natural extension to solve this problem efficiently. {\tt Newton} without the natural extension would not be deterministic and would slow down exponentially, since box-consistency on the expanded extension loses too much precision due to the dependency problem (multiple occurrences of the same variable) and box-consistency on the Taylor interval extension is not helpful initially. Once again, we observe that box-consistency over the natural extension is helpful when far from a solution while box-consistency on the Taylor extension is useful to terminate the search quickly. Table \ref{more-newton-W} gives the result for the initial intervals of size $[-10^8,0]$, which shows that the algorithm continues to perform well in this case. Finally, Table \ref{more-hansen} gives the results for the {\tt HRB} algorithm on this problem. Once again, {\tt Newton} outperforms the {\tt HRB} method substantially. \subsection{Interval Arithmetic Benchmarks} This section considers standard benchmarks from interval arithmetic papers \cite{Moore79,Hong95}. Benchmark {\tt i1} is the following set of equations \begin{small} \[ \left\{ \begin{array}{l} 0 = x_{1} - 0.25428722 - 0.18324757 \ x_{4} \ x_{3} \ x_{9} \\ 0 = x_{2} - 0.37842197 - 0.16275449 \ x_{1} \ x_{10} \ x_{6} \\ 0 = x_{3} - 0.27162577 - 0.16955071 \ x_{1} \ x_{2} \ x_{10} \\ 0 = x_{4} - 0.19807914 - 0.15585316 \ x_{7} \ x_{1} \ x_{6} \\ 0 = x_{5} - 0.44166728 - 0.19950920 \ x_{7} \ x_{6} \ x_{3} \\ 0 = x_{6} - 0.14654113 - 0.18922793 \ x_{8} \ x_{5} \ x_{10} \\ 0 = x_{7} - 0.42937161 - 0.21180486 \ x_{2} \ x_{5} \ x_{8} \\ 0 = x_{8} - 0.07056438 - 0.17081208 \ x_{1} \ x_{7} \ x_{6} \\ 0 = x_{9} - 0.34504906 - 0.19612740 \ x_{10} \ x_{6} \ x_{8} \\ 0 = x_{10} - 0.42651102 - 0.21466544 \ x_{4} \ x_{8} \ x_{1} \end{array} \right. \] \end{small} with initial intervals $[-2,2]$. Benchmark {\tt i2} is the set of equations \begin{small} \[ \left\{ \begin{array}{l} 0 = x_{1} - 0.24863995 - 0.19594124 \ x_{7} \ x_{10} \ x_{16} \\ 0 = x_{2} - 0.87528587 - 0.05612619 \ x_{18} \ x_{8} \ x_{11} \\ 0 = x_{3} - 0.23939835 - 0.20177810 \ x_{10} \ x_{7} \ x_{11} \\ 0 = x_{4} - 0.47620128 - 0.16497518 \ x_{12} \ x_{15} \ x_{1} \\ 0 = x_{5} - 0.24711044 - 0.20198178 \ x_{8} \ x_{9} \ x_{16} \\ 0 = x_{6} - 0.33565227 - 0.15724045 \ x_{16} \ x_{18} \ x_{11} \\ 0 = x_{7} - 0.13128974 - 0.12384342 \ x_{12} \ x_{13} \ x_{15} \\ 0 = x_{8} - 0.45937304 - 0.18180253 \ x_{19} \ x_{15} \ x_{18} \\ 0 = x_{9} - 0.46896600 - 0.21241045 \ x_{13} \ x_{2} \ x_{17} \\ 0 = x_{10} - 0.57596835 - 0.16522613 \ x_{12} \ x_{9} \ x_{13} \\ 0 = x_{11} - 0.56896263 - 0.17221383 \ x_{16} \ x_{17} \ x_{8} \\ 0 = x_{12} - 0.70561396 - 0.23556251 \ x_{14} \ x_{11} \ x_{4} \\ 0 = x_{13} - 0.59642512 - 0.24475135 \ x_{7} \ x_{16} \ x_{20} \\ 0 = x_{14} - 0.46588640 - 0.21790395 \ x_{13} \ x_{3} \ x_{10} \\ 0 = x_{15} - 0.10607114 - 0.20920602 \ x_{1} \ x_{9} \ x_{10} \\ 0 = x_{16} - 0.26516898 - 0.21037773 \ x_{4} \ x_{19} \ x_{9} \\ 0 = x_{17} - 0.20436664 - 0.19838792 \ x_{20} \ x_{10} \ x_{13} \\ 0 = x_{18} - 0.56003141 - 0.18114505 \ x_{6} \ x_{13} \ x_{8} \\ 0 = x_{19} - 0.92894617 - 0.04417537 \ x_{7} \ x_{13} \ x_{16} \\ 0 = x_{20} - 0.57001682 - 0.17949149 \ x_{1} \ x_{3} \ x_{11} \end{array} \right. \] \end{small} \noindent with initial intervals $[-1,2]$. Benchmark {\tt i3} has the same set of equations as {\tt i2} but has initial intervals $[-2,2]$. Benchmark {\tt i4} has the set of equations \begin{small} \[ \left\{ \begin{array}{l} 0 = x_1^2 - 0.25428722 - 0.18324757 \ x_4^2 \ x_3^2 \ x_9^2 \\ 0 = x_2^2 - 0.37842197 - 0.16275449 \ x_1^2 \ x_{10}^2 \ x_6^2 \\ 0 = x_3^2 - 0.27162577 - 0.16955071 \ x_1^2 \ x_2^2 \ x_{10}^2 \\ 0 = x_4^2 - 0.19807914 - 0.15585316 \ x_7^2 \ x_1^2 \ x_6^2 \\ 0 = x_5^2 - 0.44166728 - 0.19950920 \ x_7^2 \ x_6^2 \ x_3^2 \\ 0 = x_6^2 - 0.14654113 - 0.18922793 \ x_8^2 \ x_5^2 \ x_{10}^2 \\ 0 = x_7^2 - 0.42937161 - 0.21180486 \ x_2^2 \ x_5^2 \ x_8^2 \\ 0 = x_8^2 - 0.07056438 - 0.17081208 \ x_1^2 \ x_7^2 \ x_6^2 \\ 0 = x_9^2 - 0.34504906 - 0.19612740 \ x_{10}^2 \ x_6^2 \ x_8^2 \\ 0 = x_{10}^2 - 0.42651102 - 0.21466544 \ x_4^2 \ x_8^2 \ x_1^2 \end{array} \right. \] \end{small} \noindent and initial intervals $[-1,1]$. The number of solutions must be a multiple of 1024. Benchmark {\tt i5} has the following set of equations \begin{small} \[ \left\{ \begin{array}{l} 0 = x_1 - 0.25428722 - 0.18324757 \ x_4^3 \ x_3^3 \ x_9^3 \ + \ x_3^4 \ x_9^7 \\ 0 = x_2 - 0.37842197 - 0.16275449 \ x_1^3 \ x_{10}^3 \ x_6^3 \ + \ x_{10}^4 \ x_6^7 \\ 0 = x_3 - 0.27162577 - 0.16955071 \ x_1^3 \ x_2^3 \ x_{10}^3 \ + \ x_2^4 \ x_{10}^7 \\ 0 = x_4 - 0.19807914 - 0.15585316 \ x_7^3 \ x_1^3 \ x_6^3 \ + \ x_1^4 \ x_6^7 \\ 0 = x_5 - 0.44166728 - 0.19950920 \ x_7^3 \ x_6^3 \ x_3^3 \ + \ x_6^4 \ x_3^7 \\ 0 = x_6 - 0.14654113 - 0.18922793 \ x_8^3 \ x_5^3 \ x_{10}^3 \ + \ x_5^4 \ x_{10}^7 \\ 0 = x_7 - 0.42937161 - 0.21180486 \ x_2^3 \ x_5^3 \ x_8^3 \ + \ x_5^4 \ x_8^7 \\ 0 = x_8 - 0.07056438 - 0.17081208 \ x_1^3 \ x_7^3 \ x_6^3 \ + \ x_7^4 \ x_6^7 \\ 0 = x_9 - 0.34504906 - 0.19612740 \ x_{10}^3 \ x_6^3 \ x_8^3 \ + \ x_6^4 \ x_8^7 \\ 0 = x_{10} - 0.42651102 - 0.21466544 \ x_4^3 \ x_8^3 \ x_1^3 \ + \ x_8^4 \ x_1^7 \end{array} \right. \] \end{small} \noindent and initial intervals $[-1,1]$. \begin{table} \small \centering \begin{tabular}{||l||r|r|r|r|r||} \hline \hline &{\tt i1} &{\tt i2} &{\tt i3}&{\tt i4}&{\tt i5} \\ \hline {\it CPU time} & 0.06 & 0.30 & 0.31 & 73.94 & 0.08 \\ {\it branching} & 0 & 0 & 0 & 1023 & 0 \\ \hline {\it na-ee} & 625 & 2107 & 1949 & 290776 & 381 \\ {\it na-te} & 0 & 80 & 80 & 37930 & 21 \\ {\it na-tot} & 625 & 2187 & 2029 & 328706 & 401 \\ \hline {\it fe-ee} & 1760 & 5698 & 5318 & 752220 & 992 \\ {\it fe-te} & 0 & 560 & 560 & 269760 & 140 \\ {\it fe-tot} & 1760 & 6258 & 5878 & 1021980 & 1132 \\ \hline {\it pr-con} & 0 & 1 & 1 & 1939 & 1 \\ \hline {\it proof} & yes & yes & yes & yes & yes \\ \hline \end{tabular} \caption{{\tt Newton} on the Traditional Interval Arithmetic Benchmarks} \label{trad-newton} \end{table} {\tt Newton} solves all the problems with one solution without branching and solves the problem having 1024 solutions with 1023 branchings. Note also that box-consistency on the expanded extension solves benchmark {\tt i1} alone. The results once again confirm our observation on when the various extensions are useful. Closely related results were observed in \cite{Hong95} on these benchmarks (see the related work section for a more detailed comparison) but our algorithm is in general about 4 times faster (assuming similar machines) and does not do any branching on {\tt i5}. Table \ref{trad-hansen} also describes the results for the traditional interval arithmetic method. The importance of box-consistency on the expanded extension can easily be seen from these results. Note also that {\tt Newton} (and interval methods) can prove the existence of a solution in the final intervals for all these problems. \begin{table} \small \centering \begin{tabular}{||l||r|r|r|r|r||} \hline \hline &{\tt i1} &{\tt i2} &{\tt i3}&{\tt i4}&{\tt i5} \\ \hline {\it CPU time} & 14.28 & 1821.23 & 5640.80& 445.28 & 33.58 \\ {\it branching} & 498 & 9031 & 36933 & 11263 & 1173 \\ \hline {\it fe-tot} & 77380 & 6441640 & 19979025 &2554066 & 154948 \\ \hline {\it pr-con} & 586 & 12817 & 42193 & 14335 & 1211 \\ \hline {\it proof} & yes & yes & yes & yes & yes \\ \hline \end{tabular} \caption{{\tt HRB} on the Traditional Interval Arithmetic Benchmarks} \label{trad-hansen} \end{table} \subsection{Kinematics Applications} We now describe the performance of {\tt Newton} on two kinematics examples. Application {\tt kin1} comes from robotics and describes the inverse kinematics of an elbow manipulator \cite{Hong95}. It consists of a sparse system with 12 variables and the set of equations is as follows: \[ \left\{ \begin{array}{l} s_2 c_5 s_6 - s_3 c_5 s_6 - s_4 c_5 s_6 + c_2 c_6 + c_3 c_6 + c_4 c_6 = 0.4077 \\ c_1 c_2 s_5 + c_1 c_3 s_5 + c_1 c_4 s_5 + s_1 c_5 = 1.9115\\ s_2 s_5 + s_3 s_5 + s_4 s_5 = 1.9791\\ c_1 c_2 + c_1 c_3 + c_1 c_4 + c_1 c_2 + c_1 c_3 + c_1 c_2 = 4.0616\\ s_1 c_2 + s_1 c_3 + s_1 c_4 + s_1 c_2 + s_1 c_3 + s_1 c_2 = 1.7172\\ s_2 + s_3 + s_4 + s_2 + s_3 + s_2 = 3.9701\\ s_i^2 + c_i^2 = 1 \;\;\;\;\;\; (1 \leq i \leq 6). \end{array} \right. \] The second benchmark, denoted by {\tt kin2}, is from \cite{Morgan-kin} and describes the inverse position problem for a six-revolute-joint problem in mechanics. The equations which describe a denser constraint system are as follows: \[ \left\{ \begin{array}{l} x_i^2 + x_{i+1}^2 - 1 = 0 \;\;\;\;\; (1 \leq i \leq 4) \\ a_{1i} x_1 x_3 + a_{2i} x_1 x_4 + a_{3i} x_2 x_3 + a_{4i} x_2 x_4 + a_{5i} x_5 x_7 + a_{6i} x_5 x_8 + a_{7i} x_6 x_7 + a_{8i} x_6 x_8 \\ a_{9i} x_1 + a_{10i} x_2 + a_{11i} x_3 + a_{12i} x_4 + a_{13i} x_5 a_{14i} x_6 + a_{15i} x_7 + a_{16i} x_8 + a_{17i} = 0 \;\; (1 \leq i \leq 4) \end{array} \right. \] where the coefficients $a_{ki}$ are given in table \ref{kin-coef}. In both examples, the initial intervals were given as $[-10^8,10^8]$. \begin{table} \small \centering \begin{tabular}{||r|r|r|r||} \hline \hline - 0.249150680 & + 0.125016350 & - 0.635550070 & + 1.48947730 \\ + 1.609135400 & - 0.686607360 & - 0.115719920 & + 0.23062341 \\ + 0.279423430 & - 0.119228120 & - 0.666404480 & + 1.32810730 \\ + 1.434801600 & - 0.719940470 & + 0.110362110 & - 0.25864503 \\ + 0.000000000 & - 0.432419270 & + 0.290702030 & + 1.16517200 \\ + 0.400263840 & + 0.000000000 & + 1.258776700 & - 0.26908494 \\ - 0.800527680 & + 0.000000000 & - 0.629388360 & + 0.53816987 \\ + 0.000000000 & - 0.864838550 & + 0.581404060 & + 0.58258598 \\ + 0.074052388 & - 0.037157270 & + 0.195946620 & - 0.20816985 \\ - 0.083050031 & + 0.035436896 & - 1.228034200 & + 2.68683200 \\ - 0.386159610 & + 0.085383482 & + 0.000000000 & - 0.69910317 \\ - 0.755266030 & + 0.000000000 & - 0.079034221 & + 0.35744413 \\ + 0.504201680 & - 0.039251967 & + 0.026387877 & + 1.24991170 \\ - 1.091628700 & + 0.000000000 & - 0.057131430 & + 1.46773600 \\ + 0.000000000 & - 0.432419270 & - 1.162808100 & + 1.16517200 \\ + 0.049207290 & + 0.000000000 & + 1.258776700 & + 1.07633970 \\ + 0.049207290 & + 0.013873010 & + 2.162575000 & - 0.69686809 \\ \hline \hline \end{tabular} \caption{Coefficients for the Inverse Kinematics Example.} \label{kin-coef} \end{table} The results of {\tt Newton} on these two benchmarks are given in Table \ref{kin-newton}. {\tt Newton} is fast on the first benchmark and does not branch much to obtain all solutions. The algorithm in \cite{Hong95} branches more (the reported figure is 257 branches but it is not really comparable due to the nature of the algorithm) and is about 16 times slower on comparable machines. We are not aware of the results of continuation methods on this problem. {\tt Newton} is slower on the second application and takes about 6 minutes. The continuation method described in \cite{Cools94} requires about 30 seconds on a {\tt DEC 5000/200}. This method exploits the fact that the Newton polytopes for the last 4 equations are the same. Note that {\tt HRB} requires about 1630 and 4730 seconds on these examples. Note also that {\tt Newton} can prove the existence of solutions in the final intervals for these problems. \begin{table} \small \centering \begin{tabular}{||l||r|r|r|r|r||} \hline \hline &{\tt kin1} &{\tt kin2} \\ \hline {\it CPU time} & 14.240 & 353.06 \\ {\it branching} & 89 & 5693 \\ \hline {\it na-ee} & 17090 & 784687 \\ {\it na-te} & 10176 & 123032 \\ {\it na-tot} & 27266 & 907719 \\ \hline {\it fe-ee} & 45656 & 1714779 \\ {\it fe-te} & 62080 & 854384 \\ {\it fe-tot} & 107736 & 2569163 \\ \hline {\it pr-con} & 163 & 9505 \\ \hline {\it proof} & yes & yes \\ \hline \end{tabular} \caption{{\tt Newton} on the Kinematics Benchmarks} \label{kin-newton} \end{table} \subsection{An Economics Modelling Application} \begin{table} \small \centering \begin{tabular}{||l||r|r|r|r|r|r||} \hline \hline & 4 & 5 & 6 & 7 & 8 & 9 \\ \hline {\it CPU time} & 0.21 & 1.22 & 8.20 & 46.59 & 352.80 & 3311.42 \\ {\it growth} & & 5.80 & 6.72 & 5.68 & 7.57 & 9.38 \\ {\it branching} & 24 & 119 & 517 & 2231 & 12248 & 82579 \\ \hline {\it na-ee} & 834 & 4701 & 42481 & 214430& 1622417 & 14031838 \\ {\it na-te} & 480 & 1976 & 6325 & 29238 & 157402 & 1219960 \\ {\it na-tot} & 1314 & 6677 & 48806 & 243668& 1779819 & 15251798 \\ \hline {\it fe-ee} & 2220 & 116628& 99994 & 489894& 3626847 & 30782836 \\ {\it fe-te} & 1293 & 6304 & 29825 & 160284& 1062789 & 9118448 \\ {\it fe-tot} & 3513 & 17932 & 129819 & 650178& 4689636 & 39901284 \\ \hline {\it pr-con} & 37 & 147 & 687 & 2828 & 15265 & 104352 \\ \hline {\it proof} & yes & yes & yes & yes & yes & yes \\ \hline \end{tabular} \caption{{\tt Newton} on the economics modelling problem with initial intervals in $[-100,100]$} \label{eco-newton} \end{table} \begin{table} \small \centering \begin{tabular}{||l||r|r|r|r|r|r||} \hline \hline & 4 & 5 & 6 & 7 & 8 & 9 \\ \hline {\it CPU time} & 0.60 & 3.35 & 22.53 & 127.65 & 915.24 & 8600.28 \\ {\it growth} & & 5.58 & 6.72 & 5.66 & 7.16 & 9.39 \\ {\it branching} & 102 & 500 & 1778 & 7527 & 38638 & 244263 \\ \hline {\it na-ee} & 2689 & 15227 & 122662 & 606805& 4413150 & 36325819 \\ {\it na-te} & 978 & 3296 & 14055 & 64632 & 366436 & 2647016 \\ {\it na-tot} & 3667 & 18523 & 136717 & 671437 & 4779586 & 38972835 \\ \hline {\it fe-ee} & 7140 & 38160 & 291206 & 1400216& 9913850 & 80264897 \\ {\it fe-te} & 3288 & 15024 & 80110 & 429396 & 2825368 & 22672208 \\ {\it fe-tot} & 10428 & 53184 & 371716 &1829612 & 12739218& 102937105 \\ \hline {\it pr-con} & 148 & 527 & 2080 & 8337 & 42704 & 271534 \\ \hline {\it proof} & yes & yes & yes & yes & yes & yes \\ \hline \end{tabular} \caption{{\tt Newton} on the economics modelling problem with initial intervals in $[-10^8,10^8]$} \label{eco-newton-W} \end{table} \noindent The following example is taken from \cite{Morgan87}. It is a difficult economic modelling problem that can be scaled up to arbitrary dimensions. For a given dimension $n$, the problem can be stated as the system \[ \left\{ \begin{array}{l} (x_k + \sum_{i = 1}^{n-k-1} x_i x_{i+k}) x_n - c_k = 0 \;\;\; (1 \leq k \leq n-1) \\ \sum_{l=1}^{n-1} x_l + 1 = 0 \end{array} \right. \] and the constants can be chosen at random. Table \ref{eco-newton} reports the results for various values of $n$ with an initial interval of $[-100,100]$. It is interesting to compare those results with the continuation methods presented in \cite{Cools94}. \cite{Cools94} reports times (on a DEC-5000/200) of about 1 second for $n = 4$, 6 seconds for $n=5$, 50 seconds for $n=6$ and 990 seconds for $n=7$. {\tt Newton} is substantially faster on this problem than this continuation method, since it takes about 47 seconds for $n =7$. More importantly, the growth factor seems much lower in {\tt Newton}. The continuation method has growths of about 8 and 20 when going from 5 to 6 and 6 to 7, while {\tt Newton} has growths of about 6.72 and 5.68. Table \ref{eco-newton-W} gives the same results for initial intervals in $[-10^8,10^8]$. It is interesting to note that the computation times increase by less than a factor 3 and that the growth is consistent with the results on the smaller intervals. Note also the existence proof of solutions for these problems. \subsection{Combustion Application} \noindent This problem is also from Morgan's book \cite{Morgan87} and represents a combustion problem for a temparature of $3000^{\circ}$. The problem is described by the following sparse systems of equations \[ \left\{ \begin{array}{l} x_2 + 2 \ x_6 + x_9 + 2\ x_{10} = 10^{-5} \\ x_3 + x_8 = 3 \ 10^{-5} \\ x_1 + x_3 + 2 \ x_5 + 2 \ x_8 + x_9 + x_{10} = 5 \ 10^{-5} \\ x_4 + 2 \ x_7 = 10^{-5} \\ 0.5140437 \ 10^{-7} \ x_5 = x_1^2 \\ 0.1006932 \ 10^{-6} \ x_6 = 2 \ x_2^2 \\ 0.7816278 \ 10^{-15} \ x_7 = x_4^2 \\ 0.1496236 \ 10^{-6} \ x_8 = x_1\ x_3 \\ 0.6194411 \ 10^{-7} \ x_9 = x_1\ x_2 \\ 0.2089296 \ 10^{-14} \ x_{10} = x_1\ x_2^2 \\ \end{array} \right. \] which is typical of chemical equilibrium systems. Table \ref{combustion-newton} describes the results of {\tt Newton} on for the initial intervals $[-1,1]$ and $[-10^8,10^8]$. {\tt Newton} behaves well on this example, since the continuation method of \cite{Cools94} takes about 57 seconds. Note once again that a substantial increase in the size of the initial intervals only induces a slowdown of about 2.5 for {\tt Newton}. Note also the proof of existence of the solutions. \begin{table} \small \centering \begin{tabular}{||l||r|r||} \hline \hline & [-1,1] & [-$10^8$,$10^8$] \\ \hline {\it CPU time} & 4.06 & 9.94 \\ {\it branching} & 183 & 523 \\ \hline {\it na-ee} & 10676 & 28473 \\ {\it na-te} & 5248 & 7952 \\ {\it na-tot} &15924 & 36425 \\ \hline {\it fe-ee} & 28928 & 77508 \\ {\it fe-te} & 22464 & 49632 \\ {\it fe-tot} & 51392 & 127140 \\ \hline {\it pr-con} & 187 & 527 \\ \hline {\it proof} & yes & yes \\ \hline \end{tabular} \caption{{\tt Newton} on the Combustion Problem.} \label{combustion-newton} \end{table} \subsection{Chemical Equilibrium Application} \begin{table} \small \centering \begin{tabular}{||r|r||r|r|r||r|r|r||r|r||} \hline \hline {\it CPU time} & {\it branching} & {\it na-ee} & {\it na-te} & {\it na-tot} & {\it fe-ee} & {\it fe-te} & {\it fe-tot} & {\it pr-con} & {\it proof} \\ \hline 6.32 & 256 & 13725 & 3400 & 17125 & 34811 & 17425 & 52236 & 425 & yes \\ \hline \hline \end{tabular} \caption{{\tt Newton} on the Chemestry Problem with Initial Intervals $[0,10^8]$} \label{chemestry-newton} \end{table} \noindent This problem originates from \cite{Morgan90} and describes a chemical equilibrium system. The set of equations is as follows: \[ \left\{ \begin{array}{l} R = 10 \\ R_5 = 0.193 \\ R_6 = 0.002597 / \sqrt{40} \\ R_7 = 0.003448 / \sqrt{40} \\ R_8 = 0.00001799 / 40 \\ R_9 = 0.0002155 / \sqrt{40} \\ R_{10} = 0.00003846 / 40 \\ x_1 \ x_2 + x_1 - 3 \ x_5 = 0 \\ 2 \ x_1 \ x_2 + x_1 + x_2 \ x_3^2 + R_8 \ x_2 - R \ x_5 + 2 \ R_{10} \ x_2^2 + R_7 \ x_2 \ x_3 + R_9 \ x_2 \ x_4 = 0 \\ 2 \ x_2 \ x_3^2 + 2 \ R_5 \ x_3^2 - 8 \ x_5 + R_6 \ x_3 + R_7 \ x_2 \ x_3 = 0 \\ R_9 \ x_2 \ x_4 + 2 \ x_4^2 - 4 \ R \ x_5 = 0 \\ x_1 \ x_2 + x_1 + R_{10} \ x_2^2 + x_2 \ x_3^2 + R_8 \ x_2 + R_5 \ x_3^2 + x_4^2 - 1 + R_6 \ x_3 + R_7 \ x_2 \ x_3 + R_9 \ x_2 \ x_4 = 0 \end{array} \right. \] and all $x_i$ must be positive. The results are depicted in Table \ref{chemestry-newton} for an initial interval $[0,10^8]$. They indicate that {\tt Newton} is particularly effective on this problem, since it takes about 6 seconds and proves the existence of a solution in the final intervals. Note that the continuation method of \cite{Cools94} takes about 56 seconds on this problem. \subsection{A Neurophysiology Application} We conclude the experimental section by showing an example illustrating the limitations of {\tt Newton}. The application is from neurophysiology \cite{Cools94} and consists of the following system of equations: \[ \left\{ \begin{array}{l} x_1^2 + x_3^2 = 1 \\ x_2^2 + x_4^2 = 1 \\ x_5 \ x_3^3 + x_6 \ x_4^3 = c_1 \\ x_5 \ x_1^3 + x_6 \ x_2^3 = c_2 \\ x_5 \ x_1 \ x_3^2 + x_6 \ x_4^2 \ x_2 = c_3 \\ x_5 \ x_1^2 \ x_3 + x_6 \ x_2^2 \ x_4 = c_4 \end{array} \right. \] No initial intervals for the variables were given and the constants $c_i$ can be chosen at random. The continuation method of \cite{Cools94} solves this problem in about 6 seconds. The results of {\tt Newton} are depicted in Table \ref{neuro-newton} for various initial intervals. {\tt Newton} is fast when the initial intervals are small (i.e., $[-10,10]$). Unfortunately, the running time of the algorithm increases linearly with the size of the initial intervals, showing a limitation of the method on this example. \begin{table} \small \centering \begin{tabular}{||l||r|r|r|r||} \hline \hline & $[-10,10]$ & $[-10^2,10^2]$ & $[-10^3,10^3]$ & $[-10^4,10^4]$ \\ \hline {\it CPU time} & 0.91 & 11.69 & 172.71 & 2007.51 \\ {\it growth} & & 12.84 & 14.77 & 11.62 \\ {\it branching} & 52 & 663 & 9632 & 115377 \\ \hline {\it na-ee} & 4290 & 57224 & 645951 & 6541038 \\ {\it na-te} & 810 & 6708 & 96888 & 1173456 \\ {\it na-tot} & 5100 & 63932 & 742839 & 7714494 \\ \hline {\it fe-ee} & 11012 & 144843 & 1620538 & 16647442 \\ {\it fe-te} & 4104 & 48804 & 769056 & 9376632 \\ {\it fe-tot} & 15116 & 193647 & 2389594 & 26024074 \\ \hline {\it pr-con} & 69 & 983 & 15980 & 195270 \\ \hline {\it proof} & yes & yes & yes & yes \\ \hline \end{tabular} \caption{{\tt Newton} on the Neurophysiology Problem.} \label{neuro-newton} \end{table} \section{Related Work and Discussion} \label{related} The research described in this paper originated in an attempt to improve the efficiency of constraint logic programming languages based on intervals such as {\tt BNR-Prolog} \cite{BNR} and {\tt CLP(BNR)} \cite{CLP-BNR}. These Constraint Logic Programming languages use constraint solving as basic operation and they were based on the simple generalization of arc-consistency described previously, i.e. \[ I_i = \bo{I_i \cap \{ \ r_i \ \mid \ \exists r_1 \in I_1, \ldots, \exists r_{i-1} \in I_{i-1}, \ldots, \exists r_{i+1} \in I_{i+1}, \ldots, \exists r_{n} \in I_n: c(r_1,\ldots,r_n) \ \}}. \] This approximation was enforced on simple constraints such as \[ x_1 = x_2 + x_3, \;\; x_1 = x_2 - x_3, \;\; x_1 = x_2 \times x_3 \] and complex constraints were decomposed in terms of these simple constraints. As mentioned previously, this approach is not very effective and our main goal was to design new approximations of arc-consistency that could make use of existing interval Newton methods. The main problem was the difficulty in characterizing the pruning of the Newton operator $N^*$ in a declarative way (in order to introduce it nicely in the above programming languages) and box-consistency emerged as an attempt to generalize the operator to make sure that the bounds of the interval were locally consistent. Subsequent research made us realize that box-consistency is independent of the Newton operator and can be enforced even if the functions are not continuous or differentiable. In addition, the value of applying box-consistency on several extensions became clear. On the one hand, box-consistency on the Taylor extension generalizes interval methods based on Gauss-Seidel iterations and enables us to capture nicely Hansen-Segupta's operator. On the other hand, box-consistency on the natural and expanded extensions is really orthogonal to the pruning obtained from the Taylor expansion, producing a particularly effective algorithm. It is also worth pointing out that {\tt Newton} spends most time in the natural and expanded extensions. However, for many applications, the use of the Taylor interval extension is critical to terminate the search quickly and to avoid generating many small intervals around the solutions. As a general observation, box-consistency on the natural and expanded extensions seem effective when far from a solution while box-consistency on the Taylor expansion seems effective when near a solution. It is worth mentioning that that the interval community has spent much effort to design additional techniques to speed up further the computation when near a solution but have not considered techniques to improve pruning when far from a solution. It is interesting to note that the idea of using approximations of arc-consistency was also used independently by Hong and Stahl \cite{Hong95}, who were also exposed to research on Constraint Logic Programming. Their use of projections is however quite different from ours. The key idea is to work with a set of boxes and to use projections to split a box into several subboxes by isolating all zeros of a projection. This gives an algorithm of a very different nature which cannot easily be characterized as a branch \& prune algorithm since constraints are used to branch. Our approach seems to be more effective in practice, since their use of projections may generate many subboxes that may all need to be pruned away later on, implying much redundant work. Our approach postpones the branching until no pruning takes place and generates only subboxes when they are strictly necessary to progress. It is also very interesting to report that on all benchmarks that we tested the projection never isolates more than two zeros, showing that box-consistency may be the right trade-off here. Finally, note that their approach seems to use implicitly an expanded extension\footnote{The idea of sandwitching the interval function in between two real functions is described there.} but they do not make use of the natural extension which is very important for some applications. The research described here also provides a uniform framework to integrate these techniques in Constraint Logic Programming, to understand the importance of the various pruning operators and their relationships and to suggest further research directions. For instance, higher notions of consistency such as path-consistency \cite{Montanari74} may be worth investigating for some applications. \section{Conclusion} \label{conclusion} In this paper, we presented a branch \& prune algorithm to find all isolated solutions to a system of polynomial constraints over the reals. The algorithm is based on a single concept, box-consistency, which is an approximation of arc-consistency, a notion well-known in artificial intelligence. Box-consistency can be instantiated by to produce Hansen-Segupta operator as well as other narrowing operators which are more effective when the computation is far from a solution. The algorithm and its mathematical foundations are simple. Moreover, the algorithm is shown to behave well on a variety of benchmarks from kinematics, mechanics, chemistry, combustion, and economics. It outperforms the interval Newton methods we know of and compares well with continuation methods on their benchmarks. In addition, problems such as the Broyden banded function and the Mor\'{e}-Cosnard discretization of a nonlinear integral equation can be solved for several hundred variables. Limitations of the method (e.g. a sensitivity to the size of the initial intervals on some problems) have also been identified. \section*{Acknowledgments} We would like to thank F. Benhamou, A. Colmerauer, and B. Le Charlier for many interesting discussions on this topic. This research was partly supported by the Office of Naval Research under grant N00014-91-J-4052 ARPA order 8225, the National Science Foundation under grant numbers CCR-9357704, a NSF National Young Investigator Award. \bibliographystyle{plain} \bibliography{pascal} \newpage \appendix \end{document} Note that working with unions of intervals is called extended interval arithmetic in the interval arithmetic literature. \begin{figure} \begin{tt} \begin{tabbing} 123\=123\=123\=123\=123\=123\=\kill {\bf function} {\it ShrinkL}( $F$: $\I \rightarrow \I$; $F'$: $\I \rightarrow \I$; $I: \I$)$\fto \I$; \\ {\bf begin} \\ \> success := false; \\ \> rightI := {\it right(I)}; \\ \> {\it {\tt PUSH}(I)}; \\ \> {\bf while} $\neg$ EMPTY\_STACK $\wedge$ $\neg$ success {\bf do} \\ \>\> {\it I} := POP\_STACK; \\ \>\> $I$ := $N^*(F,F',I)$; \\ \>\> {\bf if} $0 \in F(I)$ {\bf then} \\ \>\>\> {\bf if} {\it LeftBounded(F,I)} {\bf then} \\ \>\>\>\> success := true \\ \>\>\> {\bf else} \\ \>\>\>\> PUSH([{\it center(I)},{\it right(I)}]); \\ \>\>\>\> PUSH([{\it left(I)},{\it center(I)}]);\\ \>\>\> {\bf endif} \\ \>\> {\bf endif} \\ \> {\bf endwhile} \\ \> {\bf if} $\neg$ success {\bf then} {\bf return} $\emptyset$ {\bf else return} [{\it left(I)},rightI] \\ {\bf end} \end{tabbing} \end{tt} \caption{Implementation of Function {\it ShrinkL} using Newton Iteration} \label{shr} \end{figure} In the following, we use an operation $\UDIV(I_c,I_a,I_d)$ which computes $I_c - I_a/I_d$ using unions of intervals. In addition, we take the convention that $\UDIV(I_c,I_a,I_d)$ always returns a union $I_n \cup I_p$, where $I_n$ is empty or contains only nonpositive numbers and $I_p$ is empty or contains only nonnegative numbers. \begin{definition}[Newton Reduction] The Newton reduction operator $N: ((\I \fto \I) \times (\I \fto \I) \times \I) \fto (\I \times \I)$ \[ N(F,F',I) = \UDIV(\ap{c},F(\ap{c}),F'(I)). \] where $c = \CENTER(I)$. The iterated Newton reduction operator $N^*: (\I \fto \I \times \I \fto \I \times \I) \fto \I$ is defined as \begin{tt} \begin{tabbing} 12345\= \kill $N^*(F,F',I) = I_n$ $(n \geq 1)$ where \\ \> $I_0 = I$ \\ \> $I_n \cup I_p = \UDIV(F,F',I_i)$ \\ \> $I_{i+1} = (I_i \cap I_n) \cup (I_i \cap I_p)$ \\ \> $I_n = I_{n-1}$. \end{tabbing} \end{tt} \end{definition} The following property is a condition to guarantee the correctness of the Newton reduction operators. \begin{definition}[Newton Compatibility] $F': \I \fto \I$ is are Newton-compatible with $F: \I \fto \I$ and $I$ iff \[ \forall I' \subseteq I \ \forall r \in I': \ 0 \in F(\ap{r}) \Rightarrow r \in N(F,F',I'). \] \end{definition} \noindent Assuming a function $F'$ that is Newton-compatible with $F$ and $I$, the implementation of $ShrinkL$ can be refined as shown in Figure \ref{shr}. The only change is the addition of the iterated Newton reduction operator to reduce the interval which is removed from the stack. \begin{quote} $[-0.3894150838568418705,-0.3894150838568416483]$ \\ $[-0.9210623716474195533,-0.9210623716474192201]$ \\ $[+0.5646048742322073987,+0.5646048742322289372]$ \\ $[-0.8253613366236840499,-0.8253613366236680626]$ \\ $[+0.7174034286244790337,+0.7174034286245180026]$ \\ $[-0.6966579652870119599,-0.6966579652869752115]$ \\ $[+0.8414785200543410947,+0.8414785200543567490]$ \\ $[-0.5402905702371235686,-0.5402905702370986995]$ \\ $[+0.9320048415873238845,+0.9320048415873262160]$ \\ $[-0.3624458238934301080,-0.3624458238934238907]$ \\ $[-0.8559151179237649565,-0.8559151123762789303]$ \\ $[-0.5171163509365193090,-0.5171163417544921703]$ \end{quote} When $n=5$, {\tt Newton} returns the intervals \begin{quote} $[-0.4283028646427008979,-0.4283028646427007313]$ \\ $[-0.4765965315010954329,-0.4765965315010952663]$ \\ $[-0.5196377221007547620,-0.5196377221007544289]$ \\ $[-0.5588619565270254164,-0.5588619565270250832]$ \\ $[-0.5588619565270254164,-0.5588619565270250832]$ \end{quote}