%-*-mode: text-*-
\documentstyle[12pt]{article}
\setlength{\unitlength}{0.1in}

\newtheorem{theorem}{Theorem}
\newtheorem{lemma}[theorem]{Lemma}

\title{Hierarchical Reasoning about Inequalities}
\author{Elisha Sacks\footnotemark\\
        MIT Laboratory for Computer Science\\
        545 Technology Square, Room 370\\
        Cambridge, MA 02139, USA}

\begin{document}

\maketitle

\begin{abstract}
This paper describes a program called {\sc bounder} that proves
inequalities between functions over finite sets of constraints.
Previous inequality algorithms perform well on some subset of the
elementary functions, but poorly elsewhere.  To overcome this problem,
{\sc bounder} maintains a hierarchy of increasingly complex algorithms.
When one fails to resolve an inequality, it tries the next.  This
strategy resolves more inequalities than any single algorithm.  It also
performs well on hard problems without wasting time on easy ones.  The
current hierarchy consists of four algorithms: bounds propagation,
substitution, derivative inspection, and iterative approximation.
Propagation is an extension of interval arithmetic that takes linear
time, but ignores constraints between variables and multiple occurrences
of variables.  The remaining algorithms consider these factors, but
require exponential time.  Substitution is a new, provably correct,
algorithm for utilizing constraints between variables.  The final two
algorithms analyze constraints between variables.  Inspection examines
the signs of partial derivatives.  Iteration is based on several earlier
algorithms from interval arithmetic. 
\end{abstract}

\section{Introduction}

\footnotetext{This research was supported (in part) by National
Institutes of Health Grant No. R01 LM04493 from the National Library of
Medicine and National Institutes of Health Grant No. R24 RR01320 from
the Division of Research Resources.}

This paper describes a program called {\sc bounder} that proves
inequalities between functions over all points satisfying a finite set
of {\it constraints:} equalities and inequalities between functions.  
{\sc Bounder} manipulates {\it extended elementary functions:} 
polynomials and compositions of exponentials, logarithms, trigonometric
functions, inverse trigonometric functions, absolute values, maxima, and
minima.  It tests whether a set of constraints, $C$, implies an
inequality $a\leq b$ between the extended elementary functions $a$ and
$b$ by calculating upper and lower bounds for $a-b$ over all points
satisfying $C$.  It proves the inequality when the upper bound is
negative or zero, refutes it when the lower bound is positive, and fails
otherwise.  

Previous bounding algorithms perform well on some subset of the extended
elementary functions, but poorly elsewhere.  For this reason, 
{\sc bounder} maintains a hierarchy of increasingly complex bounding
algorithms.  When one fails to resolve an inequality, it tries the next.
Although complex algorithms derive tighter bounds than simple ones for
most functions, exceptions exist.  Hence, {\sc bounder}'s hierarchy of
algorithms derives tighter bounds than even its most powerful component.
It also performs well on hard problems without wasting time on easier
ones.

The purpose of {\sc bounder} is to resolve inequalities that arise in
realistic modeling problems efficiently, not to derive deep theoretical
results.  It is an engineering utility, rather than a theorem-prover for
pure mathematics.  For this reason, it only addresses universally
quantified inequalities, which make up the majority of practical
problems, while ignoring the complexities of arbitrary quantification.
{\sc Bounder} helps {\sc plr} \cite{sacks87c} explore the qualitative
behavior of dynamic systems, such as stability and periodicity.  For
example, suppose a linear system contains symbolic parameters.  Given
constraints on the parameters, one can use {\sc bounder} to reason
about the locations of the system's poles and zeroes.  {\sc Plr} also
enhances the performance of {\sc qmr} \cite{sacks85}, a program that
derives the {\it qualitative properties} of parameterized functions:
signs of the first and second derivatives, discontinuities,
singularities, and asymptotes.

{\sc Bounder} consists of an inequality prover, a context manager, and
four bounding algorithms: bounds propagation, substitution, derivative
inspection, and iterative approximation.  The prover uses the bounding
algorithms to resolve inequalities, as described above.  First, it
reduces the original inequality to an equivalent but simpler one by
canceling common terms and replacing monotonic functions with their
arguments.  For example, $x+1\leq y+1$ simplifies to $x\leq y$,
$-x\leq-y$ to $x\geq y$, and $e^x\leq e^y$ to $x\leq y$.  The prover
only cancels multiplicands whose signs it can determine by bounds
propagation. 

The context manager organizes constraint sets in the format required by
the bounding algorithms.  The bounding algorithms derive upper and lower
bounds for a function over all points satisfying a constraint set.  The
context manager and bounding algorithms are described in the next two
sections.  The final two sections contain a review of literature and
conclusions.  I argue that current inequality provers are weak, brittle,
or inefficient because they process all inputs uniformly, whereas 
{\sc bounder} avoids these shortcomings with its hierarchical strategy.
While this paper discusses only inequality constraints and non-strict
inequalities, {\sc bounder} implements boolean combinations of
inequality constraints and strict inequalities analogously.

\section{The Context Manager}\label{s-cm}

The context manager derives an upper (lower) bound for a variable $x$
from an inequality $L\leq R$ by reformulating it as $x\leq U$ 
($x\geq U$) with $U$ free of $x$.  It derives upper and lower bounds for
$x$ from an equality $L=R$ by reformulating it as $x=U$.  Inequality
manipulation may depend on the signs of the expressions involved.  For
example, the constraint $ax\leq b$ can imply $x\leq b/a$ or $x\geq b/a$
depending on the sign of $a$.  In such cases, the context manager
attempts to derive the relevant signs from other members of the
constraint set using bounds propagation.  If it fails, it ignores the
constraint.  Constraints whose variables cannot be isolated, such as
$x\leq 2^x$, are ignored as well.  The number of variables in a
constraint is linear in its length and each variable requires linear
time to isolate.  Isolation may require deriving the signs of all the
subexpressions in the constraint.  Theorem~\ref{t-bp} implies that this
process takes linear time.  All told, processing each constraint
requires quadratic time in its length.  Subsequent complexity results
exclude this time.

Two pairs of functions form the interface between the context manager
and the bounding algorithms.  Given a variable $x$ and a set of constraints
$C$, the functions $\mbox{\sc var-lb}_C(x)$ and $\mbox{\sc var-ub}_C(x)$
return the maximum of $x$'s numeric lower bounds in $C$ and the minimum
of its numeric upper bounds.  The functions $\mbox{\sc lower}_C(x)$ and
$\mbox{\sc upper}_C(x)$ return the maximum over all lower bounds,
symbolic and numeric, and the minimum over all upper bounds.  Both 
{\sc var-lb} and {\sc lower} derive lower bounds for $x$, whereas both
{\sc var-ub} and {\sc upper} derive upper bounds.  However, {\sc lower}
and {\sc upper} produce tighter bounds then {\sc var-lb} and 
{\sc var-ub} because they take symbolic constraints into account. 
Examples of these functions appear in Table~\ref{t-constraints}.  All
four functions run in constant time once the contexts are constructed.

\begin{table}[htbp]
\caption{Bounds of \mbox{$\{a\geq 1,\,b\leq 0,\,b\geq -2,\,ab\geq -4,\,c=b\}$}}
        \label{t-constraints}
\vspace{.5ex}
$\begin{array}{l|llll}
  & \mbox{\sc var-lb} & \mbox{\sc var-ub} & \sc lower & \sc upper\\
\hline
a & 1 & \infty & 1 & -4/b\\
b & -2 & 0 & \max\{-2,-4/a,c\} & \min\{0,c\}\\
c & -\infty & \infty & b & b
\end{array}$
\end{table}

\section{Bounding Algorithms}
This section contains the details of the bounding algorithms.  Each
derives tighter bounds than its predecessor, but takes more time.  Each
invokes all of its predecessors for subtasks, except that derivative
inspection never calls bounds propagation.  The bounding algorithms
define the extended elementary functions on the extended real numbers in
the standard fashion, so that $1/\pm\infty=0$, $\,\log 0=-\infty$,
$\,2^\infty=\infty$, and so on.  Throughout this paper, ``number'' refers
to an extended real number.

\subsection{Bounds Propagation}\label{s-bp}
The bounds propagation algorithm ({\sc bp}) bounds a compound function
by bounding its components recursively and combining the results.  For
example, the upper bound of a sum is the sum of the upper bounds of its
addends.  The recursion terminates when it reaches numbers and
variables.  Numbers are their own bounds, while {\sc var-lb} and 
{\sc var-ub} bound variables.  Figures~\ref{f-ub} and~\ref{f-expt-ub} 
contain the upper bound algorithm, $\mbox{\sc ub}_C(e)$, for a function
$e$ over a set of constraints $C$.  The lower bound algorithm,
$\mbox{\sc lb}_C(e)$, is analogous.  One can represent $e$ as an
expression in its variables $x_1,\ldots,x_n$ or as a function $e(x)$ of
the vector $x=(x_1,\ldots,x_n)$.  From here on, these forms are used
interchangeably.  The {\sc trig-ub} algorithm, not shown here, uses
periodicity and monotonicity information to derive upper bounds for
trigonometric functions and their inverses.

\begin{figure}[htbp]
\begin{tabular}{ll}
\underline{$e$ is} & \underline{$\mbox{\sc ub}_C(e)$}\\
a number & $e$\\
a variable & $\mbox{\sc var-ub}_C(e)$\\
$a+b$ & $ub_a+ub_b$\\
$ab$ & $\max\left\{lb_a lb_b,lb_a ub_b,ub_a lb_b,ub_a ub_b\right\}$\\
$a^b$ & $\mbox{\sc expt-ub}_C(a,b)$\\
$\min\{a,b\}$ & $\min\left\{ub_a,ub_b\right\}$\\
$\max\{a,b\}$ & $\max\left\{ub_a,ub_b\right\}$\\
$\log a$ & $\log ub_a$\\
$|a|$ & $\max\left\{|lb_a|,|ub_a|\right\}$\\
trigonometric & $\mbox{\sc trig-ub}_C(e)$
\end{tabular}

\caption{The $\mbox{\sc ub}_C(e)$ algorithm; $lb_e$ and $ub_e$ abbreviate 
         $\mbox{\sc lb}_C(e)$ and $\mbox{\sc ub}_C(e)$.}\label{f-ub} 
\end{figure}

\begin{figure}[htbp]
\begin{tabular}{lll}
\multicolumn{2}{l}{\underline{Case}} & \underline{$\mbox{\sc expt-ub}_C(a,b)$}\\
\multicolumn{2}{l}{\rule{0em}{1.2em}$\mbox{\sc ub}_C(a)>0$} 
   & $\displaystyle e^{\mbox{\sc ub}_C(b\log a)}$\\
\multicolumn{2}{l}{$b=\frac{p}{q}$ with $p,q$ integers}\\
  & $p,q$ odd and positive & $\left[\mbox{\sc ub}_C(a)\right]^b$\\[.5ex]
  & $p,q$ odd and $\mbox{\it ub}_a<0$ & $\left[\mbox{\sc lb}_C(a)\right]^b$\\[.5ex]
  & $p$ even & $\displaystyle e^{\mbox{\sc ub}_C(b\log |a|)}$\\
  & else & $\infty$\\
\multicolumn{2}{l}{else} & $\infty$
\end{tabular}

\caption{The $\mbox{\sc expt-ub}_C(a,b)$ algorithm}\label{f-expt-ub}
\end{figure}

The correctness and complexity of {\sc bp} are summarized in the
theorem:
\begin{theorem}\label{t-bp}
For any extended elementary function $e(x)$ and set of constraints $C$,
bounds propagation derives numbers $lb_e$ and $ub_e$ satisfying
\begin{equation}
\forall x. \mbox{\it satisfies\/}(x,S)\Rightarrow 
           lb_e\leq e(x) \leq ub_e
\end{equation}
in time proportional to $e$'s length.

\end{theorem}
The proof is by induction on $e$'s length.  It appears in a longer
version of this paper \cite{sacks87}, as do all subsequent proofs.

Bounds propagation achieves linear time-complexity by ignoring
constraints among variables or multiple occurrences of a variable in
an expression.  It derives excessively loose bounds when these factors
prevent all the constituents of an expression from varying independently
over their ranges.  For instance, the constraint $a\leq b$ implies that
$a-b$ cannot be positive.  Yet given only this constraint, {\sc bp}
derives an upper bound of $\infty$ for $a-b$ by adding the upper bounds
of $a$ and $-b$, both $\infty$.  As another example, when no constraints
exist, the joint occurrence of $x$ in the constituents of $x^2+x$
implies a global minimum of $-1/4$.  Yet {\sc bp} deduces a lower bound
of $-\infty$ by adding the lower bounds of $x^2$ and $x$, 0 and
$-\infty$.  Subsequent bounding algorithms derive optimal bounds for
these examples.  Substitution analyzes constraints among variables and
the final two algorithms handle multiple occurrences of variables.  All
three obtain better results than {\sc bp}, but pay an exponential
time-complexity price.

\subsection{Substitution}
The substitution algorithm constructs bounds for an expression by
replacing some of its variables with their bounds in terms of the other
variables.  Substitution exploits all solvable constraints, whereas
bounds propagation limits itself to constraints between variables and
numbers.  In our previous example, substitution derives an upper bound
of 0 for $a-b$ from the constraint $a\leq b$ by bounding $a$ from above
with $b$, that is $a-b\leq b-b=0$.  Substitution is performed by the
algorithms $\mbox{\sc sup}_C(e,H)$ and $\mbox{\sc inf}_C(e,H)$, which
calculate upper and lower bounds on $e$ over the constraint set $C$ in
terms of the variable set $H$.  When $H$ is empty, the bounds reduce to
numbers.

Figures~\ref{f-sup} and~\ref{f-supp} contain the {\sc sup} function
and its auxiliary, {\sc supp}.  The auxiliary functions {\sc expt-sup}
and {\sc trig-sup} are derived from {\sc bp}'s exponential and
trigonometric bounding algorithms by replacing $\mbox{\sc ub}_C(a)$ with
$\mbox{\sc sup}_C(a,H)$, $\mbox{\sc lb}_C(a)$ with 
$\mbox{\sc inf}_C(a,H)$, and so on for $b$.  The expression $v(e)$
denotes the variables contained in $e$ and $f(b,a,H)$ abbreviates
$v(b)-v(a)\subseteq H$.  In the remainder of this section, we will focus
on {\sc sup}.  {\sc Inf} is analogous.

{\setlength{\arraycolsep}{1.5pt}
\begin{figure}[htb]
$\begin{array}{lll}
  & \underline{\mbox{$e$ is}} & \underline{\mbox{\sc sup}_C(e,H)}\\
1 & v(e)\subseteq H & e\\
2 & \mbox{a variable} & \mbox{\sc supp}_C(e,s(\mbox{\sc upper}_C(e), H\cup\{e\}))\\ 3 & a+b\\
  & \mbox{3.1 } f(b,a,H) & s(a,H)+s(b,H)\\
  & \mbox{3.2 else} & s(a+s(b,H\cup v(a)),H)\\
4 & ab\\
  & \mbox{4.1 } \mbox{\sc lb}_C(a)\geq 0\\
  &\hspace{1em}f(b,a,H) & \max\left\{s(a,H)s(b,H), i(a,H)s(b,H)\right\}\\
  &\hspace{1em}else & s(a\,s(b,H\cup v(a)),H)\\
  & \mbox{4.2 } \mbox{\sc ub}_C(a)\leq 0\\
  &\hspace{1em}f(b,a,H) & \max\left\{s(a,H)i(b,H), i(a,H)i(b,H)\right\}\\
  &\hspace{1em}else & s(a\,i(b,H\cup v(a)),H)\\
  & \mbox{4.3 else}
            & \!\!\begin{array}[t]{l}\max\left\{s(a,H)s(b,H), s(a,H)i(b,H),\right.\\
		                      \hfill\left.i(a,H)s(b,H), i(a,H)i(b,H)\right\}
                   \end{array}\\
5 & a^b & \mbox{\sc expt-sup}_C(a,b,H)\\
6 & \min\{a,b\} & \min\left\{s(a,H), s(b,H)\right\}\\
7 & \max\{a,b\} & \max\left\{s(a,H), s(b,H)\right\}\\
8 & \log a & \log s(a,H)\\
9 & |a| & \max\left\{\left|i(a,H)\right|,
                        \left|s(a,H)\right|\right\}\\
10 & \mbox{trig} & \mbox{\sc trig-sup}_C(e,H)
\end{array}$
\caption{The $\mbox{\sc sup}_C(e,H)$ algorithm.  The symbols $s$ and $i$ abbreviate 
$\mbox{\sc sup}_C$ and $\mbox{\sc inf}_C$ respectively.}\label{f-sup}
\end{figure}}

\begin{figure}[htb]

$\begin{array}{lll}
  & \underline{\mbox{case}} & \underline{\mbox{\sc supp}_C(x,B)}\\
1 & x\not\in v(B) & B\\
2 & B=rx+A\\
  & r\in\Re,\;v\not\in v(A)\\
  & 2.1\;\; r\geq 1 & \infty\\
  & 2.2\;\; r<1 & \frac{A}{1-r}\\
3 & B=\min\{C,D\} & \min\left\{\mbox{\sc supp}_C(x,C),\mbox{\sc supp}_C(x,D)\right\}\\
4 & B=\max\{C,D\} & \max\left\{\mbox{\sc supp}_C(x,C),\mbox{\sc supp}_C(x,D)\right\}\\
5 & \mbox{else} & \mbox{\sc ub}_C(B)
\end{array}$

\caption{The $\mbox{\sc supp}_C(x,B)$ algorithm}\label{f-supp}
\end{figure}

In step~1, {\sc sup} calculates the upper bounds of numbers and of
variables included in $H$.  It analyzes a variable, $x$, not in $H$ by
constructing an intermediate bound
\begin{equation}
B=\mbox{\sc sup}_C(\mbox{\sc upper}_C(x),H\cup\{x\})
\end{equation}
for $x$ and calling {\sc supp} to derive a final bound.  If possible,
{\sc supp} derives an upper bound for $x$ in $H$ directly from the
inequality $x\leq B$.  Otherwise, it applies bounds propagation to $B$.
For instance, the inequality $x\leq 1-x$ yields a bound of $1/2$, but 
$x\leq x^2-1$ does not provide an upper bound, so {\sc supp} returns
$\mbox{\sc ub}_C(x^2-1)$.

{\sc Sup} exploits constraints among variables to improve its bounds
on sums and products.  If $b$ contains variables that $a$ lacks, but
which have bounds in $a$'s variables, {\sc sup} constructs an
intermediate upper bound, $U$, for $a+b$ or $ab$ by replacing $b$ with
these bounds.  A recursive application of {\sc sup} to $U$ produces a
final upper bound.  (Although not indicated explicitly in
Figure~\ref{f-sup}, these steps are symmetric in $a$ and $b$.)  If $a$
and $b$ have the same variables, {\sc sup} bounds $a+b$ and $ab$ by
recursively bounding $a$ and $b$ and applying bounds propagation to the
results.  For example, given the constraints $c\geq 1$, $d\geq 1$, and
$cd\leq 4$, {\sc sup} derives an intermediate bound of $3c/4$ for
$c-1/d$ by replacing $-1/d$ with $-c/4$, its upper bound in $c$.  This
bound is derived as follows:
\begin{equation}
\mbox{\sc sup}(-\frac{1}{d},\{c\})=
-\mbox{\sc inf}(\frac{1}{d},\{c\})=
\frac{-1}{\mbox{\sc sup}(d,\{c\})}=
-\frac{c}{4}
\end{equation}
{\sc Sup} uses the recursive call $\mbox{\sc sup}(c,\{\})=4$ to
derive a final bound of 3 for $c-1/d$.  The following theorem
establishes the correctness of substitution:
\begin{theorem}\label{t-sup}
For every extended elementary function $e$, variable set $H$, and
constraint set $C$, the expressions \mbox{$i=\mbox{\sc inf}_C(e,H)$} and 
$s=\mbox{\sc sup}_C(e,H)$ satisfy the conditions:
\begin{eqnarray}
\mbox{$i$ and $s$ are expressions in $H$}\\
\forall x. \mbox{\it satisfies\/}(x,S)\Rightarrow 
	         i(x)\leq e(x)\leq s(x)
\end{eqnarray}
\end{theorem}

Substitution utilizes constraints among variables to improve on the
bounds of {\sc bp}, but ignores constraints among multiple occurrences
of variables.  It performs identically to {\sc bp} on the example of
$x^2+x$, deriving a lower bound of $-\infty$.  Yet that bound is overly
pessimistic because no value of $x$ minimizes both addends
simultaneously. The last two bounding algorithms address this
shortcoming.

\subsection{Derivative Inspection}
Derivative inspection calculates bounds for a function over a constraint
set $C$ from the signs of its partial derivatives.  Let us define the
{\it range} of $x_i$ in $C$ as the interval  
\begin{equation}
X_i=\left[\mbox{\sc inf}_C(x_i,\{\}), \mbox{\sc sup}_C(x_i,\{\})\right]
\end{equation}
and the {\it range} of $x=(x_1,\ldots,x_n)$ in $C$ as the Cartesian
product $X=X_1\times\cdots\times X_n$ of its components' ranges.
Derivative inspection splits the range of a function $f(x)$ into
subregions by dividing the range of each $x_i$ into maximal intervals on
which $\partial f/\partial x_i$ is non-negative, non-positive, or of
unknown sign.  The maximum upper bound over all subregions bounds $f$
from above on $X$.  This bound is valid over all points satisfying $C$
by  Theorem~\ref{t-sup}.  Each region can be collapsed to the upper
(lower) bound of $x_i$ in every dimension $i$ where
$\partial f/\partial x_i$ is non-negative (non-positive) without
altering $f$'s upper bounds.  An analogous procedure derives lower
bounds.  

Derivative inspection takes time proportional to the number of regions
into which $f$'s domain splits.  For this reason, it only applies to
functions whose partial derivatives all have finitely many zeroes in
$X$.  When the signs of all partial derivatives are known, derivative
inspection yields optimal bounds directly, since all regions reduce to
points.  For example, it derives an optimal lower bound of $-1/4$ for
$x^2+x$ because the derivative of $x^2+x$ is non-positive on
$[-\infty,-1/2]$ and non-negative on $[-1/2,\infty]$.  Otherwise, one
must use a second bounding algorithm to calculate bounds on the
non-trivial subregions.  This two-step approach generally yields tighter
bounds than applying the second algorithm directly on $f$'s entire
domain, since the subregions are smaller and often reduce to points
along some dimensions. 

\subsection{Iterative Approximation} 
Iterative approximation, like derivative inspection, reduces the errors
in bounds propagation and substitution caused by multiple occurrences of
variables.  Instead of bounding a function over its entire range
directly, it subdivides the regions under consideration and combines the
results.  Intuitively, {\sc bp}'s choice of multiple worst case values
for a variable causes less damage on smaller regions because all these
values are less far apart.  Figure~\ref{f-iademo} illustrates this idea
for the function $x^2-x$ on the interval $[0,1]$.  Part~(a) demonstrates
that {\sc bp} derives an overly pessimistic lower bound on $[0,1]$
because it minimizes both $-x$ and $x^2$ independently.  Part~(b) shows
that this factor is less significant on smaller intervals: the maximum
of the two lower bounds, $-3/4$, is a tighter bound for $x^2-x$ on
$[0,1]$ than that of part~(a).  One can obtain arbitrarily tight bounds
by constructing sufficiently fine partitions. 

\begin{figure}[htbp]
\begin{picture}(34,15)(8,0)

\thicklines

\put(20,0){(a)}
\put(46,0){(b)}
\put(20,5){\makebox(0,0){$-1$}}
\put(40,5){\makebox(0,0){$-\frac{1}{2}$}}
\put(52,5){\makebox(0,0){$-\frac{3}{4}$}}

\put(10,10){\line(1,0){20}}
\put(10,9){\makebox(0,0)[t]{0}}
\put(10,10){\circle*{.5}}
\put(10,11){\makebox(0,0)[b]{$\bf n$}}
\put(30,9){\makebox(0,0)[t]{1}}
\put(30,10){\circle*{.5}}
\put(30,11){\makebox(0,0)[b]{$\bf m$}}

\put(35,10){\line(1,0){10}}
\put(35,9){\makebox(0,0)[t]{0}}
\put(35,10){\circle*{.5}}
\put(35,11){\makebox(0,0)[b]{$\bf n$}}
\put(45,9){\makebox(0,0)[t]{$\frac{1}{2}$}}
\put(45,10){\circle*{.5}}
\put(45,11){\makebox(0,0)[b]{$\bf m$}}

\put(47,10){\line(1,0){10}}
\put(47,9){\makebox(0,0)[t]{$\frac{1}{2}$}}
\put(47,10){\circle*{.5}}
\put(47,11){\makebox(0,0)[b]{$\bf n$}}
\put(57,9){\makebox(0,0)[t]{1}}
\put(57,10){\circle*{.5}}
\put(57,11){\makebox(0,0)[b]{$\bf m$}}

\end{picture}
\caption{Illustration of iterative approximation on $[0,1]$.  The
symbols $m$ and $n$ mark the values of $x$ that minimize $-x$ and $x^2$.
The numbers below are $\mbox{\sc lb}(x^2-x)$.}\label{f-iademo}
\end{figure}

Iterative approximation generalizes interval subdivision to multivariate
functions and increases its efficiency, using ideas from Moore
\cite{moore} and Asaithambi et al.\ \cite{asaithambi82}.  As an
additional optimization, it bounds functions over the regions generated
by derivative inspection, rather than over their entire domains.  Let
$f(x_1,\ldots,x_n)$ be continuously differentiable on a region $X$ and
let $w_i$ denote the width of the interval $X_i$.  For every positive
$\epsilon$, iterative approximation derives an upper bound for $f$ on
$X$ that exceeds the least upper bound by at most $\epsilon$ within 
\begin{equation}
\left(\frac{L}{\epsilon}\right)^n \prod_{i=1}^n w_i 
\end{equation}
iterations, where the constant $L$ depends on $f$ and $X$.

\section{Related Work}
In this section, I discuss, in order of increasing generality, existing
programs that derive bounds and prove inequalities.  As one would
expect, the broader the domain of functions and constraints, the slower
the program.  The first class of systems bounds linear functions subject
to linear constraints.  Vald\'{e}s-P\'{e}rez \cite{Perez86} analyzes
sets of {\it simple linear inequalities} of the form  $x-y\geq n$ with
$x$ and $y$ variables and $n$ a number.  He uses graph search to test
their consistency in $cv$ time for $c$ constraints and $v$ variables.
Malik and Binford \cite{malik83} and Bledsoe \cite{bledsoe75} check sets
of general linear constraints for consistency and calculate bounds on
linear functions over consistent sets of constraints.  Both methods
require exponential time.%
\footnote{The simplex algorithm often performs better in practice.
Also, a polynomial alternative exists.}
The former uses the Simplex algorithm, whereas the latter introduces
preliminary versions of {\sc bounder}'s substitution algorithms.
Bledsoe defines {\sc sup, supp, inf,} and {\sc inff} for linear
functions and constraints and proves the linear version of
Theorem~\ref{t-sup}.  In fact, these algorithms produce {\it exact}
bounds, as Shostak \cite{shostak} proves.

The next class of systems bounds nonlinear functions, but allows only
range constraints.  All resemble {\sc bounder}'s bounds propagation and
all stem from Moore's \cite{moore} {\it interval arithmetic}.  Moore
introduces the rules for bounding elementary functions on finite domains
by combining the bounds of their constituents.  His algorithm takes
linear time in the length of its input.  Bundy \cite{bundy84} implements
an {\it interval package} that resembles {\sc bp} closely.  It
generalizes the combination rules of interval arithmetic to any function
that has a finite number of extrema.  If the user specifies the sign of
a function's derivative over its domain, Bundy's program can perform
interval arithmetic on it.  Unlike {\sc bounder}'s derivative inspection
algorithm, it cannot derive this information for itself.  Many other
implementations of interval arithmetic exist, some in hardware.  

Moore also proposes a simple form of iterative approximation, which
Skelboe \cite{skelboe74}, Asaithambi et al.\ \cite{asaithambi82}, and
Ratschek and Rokne \cite[ch.\ 4]{ratschek84} improve.  {\sc Bounder}'s
iterative approximation algorithm draws on all these sources.

Simmons \cite{simmons86} handles functions and constraints containing
numbers, variables, and the four arithmetic operators.  He augments
interval arithmetic with simple algebraic simplification and inequality
information.  For example, suppose $x$ lies in the interval $[-1,1]$.
Simmons simplifies $x-x$ to 0, whereas interval arithmetic produces the
range $[-2,2]$.  He also deduces that $x\leq z$ from the constraints
$x\leq y$ and $y\leq z$ by finding a path from $x$ to $z$ in the graph
of known inequalities.  The algorithm is linear in the total number of
constraints.  Although more powerful than {\sc bounder}'s bounds
propagation, Simmons's program is weaker than substitution.  For
example, it cannot deduce that $x^2\geq y^2$ from the constraints
$x\geq y$ and $y\geq 0$.

Brooks \cite[sec.\ 3]{brooks81} extends Bundy's {\sc sup} and {\sc inf}
to nonlinear functions and argues informally that Theorem~\ref{t-sup}
hold for his algorithms.  This argument must be faulty because his
version of $\mbox{\sc sup}_H(e,\{\})$ recurses infinitely when $e$
equals $x+1/x$ or $x+x^2$, for instance.  Brooks's program only exploits
constraints among the variables of sums $rx+B$ and of products $x^nB$
with $r$ real, $x$ a variable of known sign, $B$ an expression free of
$x$, and $n$ an integer.  In other cases, it adds or multiplies the
bounds of constituents, as in steps 3.1, 4.1.1, 4.2.1, and 4.3 of 
{\sc bounder}'s {\sc sup} (Figure~\ref{f-sup}).  These overly
restrictive conditions rule out legitimate substitutions that steps 3.2,
4.1.2, and 4.2.2 permit.  For example, {\sc bounder} can deduce that 
$1/x-1/y\geq 0$ from the constraints $y\geq x$ and $x\geq 1$, but
Brooks's algorithm cannot.  On some functions and non-empty sets $H$,
his algorithm makes recursive calls with $H$ empty.  This produces
needlessly loose bounds and sometimes causes an infinite recursion.

Bundy and Welham \cite{bundy79} derive upper bounds for a variable $x$
from an inequality $L\leq R$ by reformulating it as $x\leq U$ with $U$
free of $x$.  If $U$ contains a single variable, they try to find its
global maximum, $M$, by inspecting the sign of its second derivative at
the zeroes of its first derivative.  When successful, they bound $x$
from above with $M$.  Lower bounds and strict inequalities are treated
analogously.  They use a modified version of the {\sc press} equation
solver \cite{bundy81} to isolate $x$.  As discussed in
section~\ref{s-cm}, inequality manipulation depends on the signs of the
expressions involved.  When this information is required, they use
Bundy's interval package to try to derive it.  The complexity of this
algorithm is unclear, since {\sc press} can apply its simplification
rules repeatedly, possibly producing large intermediate expressions.
{\sc Bounder} contains both steps of Bundy and Welham's bounding
algorithm: its context manager derives bounds on variables from
constraints, while its derivative inspection algorithm generalizes
theirs to multivariate functions.  {\sc Press} may be able to exploit
some constraints that {\sc bounder} ignores because it contains a
stronger equation solver than does {\sc bounder}.

The final class of systems consists of theorem provers for predicate
calculus that treat inequalities specially.  These systems focus on
general theorem proving, rather than problem-solving.  They handle
more logical connectives than {\sc bounder}, including disjunction and 
existential quantification, but fewer functions, typically just
addition.  Bledsoe and Hines \cite{bledsoe80} derive a restricted form
of resolution that contains a theory of dense linear orders without
endpoints.  Bledsoe et al.\ \cite{bledsoe83} prove this form of
resolution complete.  Finally, Bledsoe et al.\
\cite{Bledsoe-Bruell-Shostak} extend a natural deduction system with
rules for inequalities.  Although none of these authors discuss
complexity, all their algorithms must be at least exponential.

\section{Conclusions}
Current inequality reasoners are weak, brittle, or inefficient because
they process all inputs uniformly.  Interval arithmetic systems, such as
Bundy's and Simmons's, run quickly, but generate exceedingly pessimistic
bounds when dependencies exist among the components of functions.  These
dependencies are caused by constraints among variables or multiple
occurrences of a variable, as discussed in Section~\ref{s-bp}.  The
upper bound of $a-b$ given $a\leq b$ demonstrates the first type, while
the lower bound of $x^2+x$ given no constraints demonstrates the second.
Each of the remaining systems is brittle because it takes only one type
of dependency into account.  Iterative approximation, suggested by
Moore, and derivative inspection, performed in the univariate case by
Bundy and Welham, address the second type of dependency, but ignore the
first.  Conversely, substitution, used (in a limited form) by Brooks and
Simmons, exploits constraints among variables, while ignoring multiple
occurrences of variables.  All these systems are inefficient because
they apply a complex algorithm to every input without trying a simple
one first.

{\sc Bounder} overcomes the limitations of current inequality reasoners
with its hierarchical strategy.  It uses substitution to analyze
dependencies among variables and derivative inspection and iterative
approximation to analyze multiple occurrences of variables.  Together,
these techniques cover far more cases than any single-algorithm system.
Yet unlike those systems, {\sc bounder} does not waste time applying
overly powerful methods to simple problems.  It tries bounds
propagation, which has linear time-complexity, before resorting to its
other methods.  An inequality reasoner like {\sc bounder} should be an
important component of future general-purpose symbolic algebra packages.

\bibliography{}
\bibliographystyle{plain}
\end{document}
