Abstract: The existence of a quantifier elimination algorithm for real arithmetic is one of the foundational results that enables formal reasoning and verification of CPS. Most of the well-known algorithms for quantifier elimination are extremely complicated, too inefficient to be used on even the simplest of formulae, or both. This makes studying quantifier elimination algorithms difficult. This project aims to rectify this problem by providing a writeup and implementation of the Cohen-Hörmander algorithm along with visualizations to aid understanding.
The modeling of cyber-physical systems (CPS), and the subsequent formal verification of the model, are made possible by a multitude of results. Foremost among these is the development of a logic (such as differential dynamic logic) that can express desirable properties of a CPS as logical formulae, along with a set of inference rules that can be used to construct proofs of these formulae. Manually constructing these computationally-verifiable proofs from the axioms of all the formal logics involved would be far too painful even for simple CPS. It thus becomes important to seek methods of automating the proof construction.
While it is impossible for many useful logics to fully automate the proof construction process [3], it is certainly possible to automate portions
of it. Rather surprisingly, a 1931 result by Tarski [7], along with related more recent developments [5], show that the entire proof construction
process can be automated once the desired goal has been reduced to proving a formula of real arithmetic. This result is absolutely foundational for CPS verification, both
from a theoretical and a practical perspective. The existence of an automatic verification procedure for formulae of real arithmetic allows one to abstract away the formal reasoning
process behind the real arithmetic proof goals, and simply give inference rules such as the following, which holds whenever
Practically speaking, the real arithmetic proof goals that result from attempts to prove properties of CPS are often prohibitively complex for manual methods.
Given the significance of this result, and the rather mysterious nature of the real arithmetic proof rule, the question "what is really going on here?" likely crosses many students' minds. A little research would reveal a number of algorithms for automatically deciding the truth of a sentence of real arithmetic (called quantifier elimination (QE) algorithms), but many of the choices have significant disadvantages for someone studing real QE for the first time:
- Tarski's original algorithm was a very important theoretical breakthrough, but complicated and so inefficient that it isn't useful for anything besides theoretical purposes [4]. Given the complexity of this algorithm, understanding it would be difficult, and given the inefficiency, interaction with an implementation (which would be very useful for understanding) would not be feasible.
- Cylindrical-Algebraic Decomposition (CAD) is the state of the art when it comes to practical QE, so it doesn't suffer from the inefficiency problem of Tarski's algorithm [2]. However, it is incredibly complicated: so much so that it took experts in the field 30 years to produce a working implementation (citation needed). As such, it is likely not suitable for a student in an introductory CPS class.
-
Virtual substitution is efficient [9], and simple enough to be part of CMU's introductory 15-424 Logical Foundations of Cyber-Physical Systems course. The only shortcoming of this algorithm is that it isn't complete, in the sense that there are theoretical limitations (given by the well-known Abel-Ruffini theorem) that prevent it from deciding the truth of arbitrary sentences of real arithmetic. Understanding this algorithm is thus not equivalent to understanding what's going on behind the scenes of the
$\R$ proof rule.
However, there is a (not too well-known) alternative that offers a reasonable balance: the Cohen-Hörmander algorithm [1]. It is simple enough to be described in this paper, complete in the sense that it can (in principle) decide the truth of any sentence of real arithmetic, and efficient enough to admit implementations that one can actually interact with. This work thus aims to introduce an audience of students taking introductory logic courses to real quantifier-elimination by providing a writeup, a number of visuals, and an implementation of the Cohen-Hörmander algorithm.
[10] develops a tool that can be used to visualize the data structures computed by the cylindrical algebraic decomposition algorithm. This present work is analogous to that one in the sense that we provide an implementation that allows users to view and understand the main data structure (the sign matrix) computed by the Cohen-Hörmander algorithm.
[4] presents the Cohen-Hörmander algorithm and provides an implementation in OCaml. Our presentation of the algorithm is based on this one. While the implementation provided by this book is likely more readable than ours, we improve accessibility and usability of the implementation by embedding it in a website and enabling verbose output that allows the reader to trace the operation of the algorithm.
Our presentation of the algorithm also differs from both of the above in that it includes animated visualizations intended to complement text-based explanations.
In this section, we very briefly review some of the background material that is necessary to properly define the problem solved by the Cohen-Hörmander algorithm.
The first-order theory of real closed fields is a formal language for stating properties of the real numbers. The language is built up recursively as follows:
Terms: Terms are the construct that the language uses to refer to real numbers, or to combine existing numbers into other ones. They are built up via the following inference rules
Formulae: Formulae are the construct that the language uses to express assertions about real numbers. The basic, or atomic, formulae are constructed as follows:
Formulae can be joined together using boolean connectives:
Finally, variables occuring in terms can be given meaning by way of quantifiers.
Real arithmetic is what we get when we give these syntactic constructs their natural meaning over
the real numbers. That is, the symbols
With these constructs, we can formally express many properties of the real numbers. Some examples include:
- Every number is positive, negative, or zero:
$\forall x, (x > 0 \vee x < 0 \vee x = 0)$ . - Every number has an additive inverse:
$\forall x, \exists y, (x + y = 0)$ . - Every nonzero number has a multiplicative inverse:
$\forall x, (x = 0 \vee \exists y, (x \cdot y = 1))$ .
However, not everything that we intuitively think of as a property of the real numbers can actually be accurately expressed in this language. A typical example is the supremum property: the assertion that every nonempty set of real numbers which is bounded above has a least upper bound has no equivalent in this language [8]. As we shall shortly see, the expressiveness (or lack thereof) of this language is key to the operation of the Cohen-Hörmander algorithm.
Now we can properly define what the Cohen-Hörmander algorithm actually does. It is a quantifier-elimination
algorithm: it takes as input a formula
-
$\exists y, (x < y \wedge y \leq 0)$ might be reduced to$x < 0$ , because regardless of which value in$\R$ we choose to assign to$x$ , the two formulaes are either both true or both false. -
$\forall x, (x > 0 \vee x < 0 \vee x = 0)$ might be reduced to$\top$ , the formula which is always true. Indeed, since the original formula has no free variables, the quantifier-eliminated formula also cannot have free variables, and thus must be equivalent to either$\top$ (true) or$\bot$ (false).
In this section, we list a few definitions and theorems of basic analysis/algebra that are useful in understanding the Cohen-Hörmander algorithm.
Definition (Sign):The sign of a real number
-
$+$ , or positive, when$x > 0$ -
$0$ when$x = 0$ -
$-$ , or negative, when$x < 0$
Theorem (Intermediate value): If
Theorem (Mean value): If
Definition (Polynomials with rational coefficients):
Theorem (Polynomial division): If
$a = bq + r$ -
$\mathrm{deg}(r) < \mathrm{deg}(b)$ .
Our approach will be to first treat the simpler univariate case: formulas of the form
A priori, deciding the truth of a formula of the form
The first step in unlocking a decision procedure for real arithmetic is to look closely at what formulae in this language can express. All formulae are ultimately built up from atomic formulae, and atomic formulae are built out of terms, so we'll start there.
Recall the inductive construction of terms: we started with rational constants
and variables (in our present case, only
Atomic formula were constructed as
This realization is key to transforming our infinite loop over
A quantifier-free formula
More formally, let
i.e, the singleton sets at the roots and the intervals between them. The entry of the sign matrix at row
Here's an example of a sign matrix for the set of polynomials
And here's an animation that illustrates the meaning of the sign matrix.
Finally, here's an animation that shows how this sign matrix can be used to compute the truth value
of the formula
One important thing to note is that while the sign matrix relies crucially on the ordering of the
roots of the polynomials involved, it doesn't actually contain any numerical information about the roots
themselves. In our toy example, it's easy to see that
Unfortunately, building the sign matrix for an arbitrary set of polynomials isn't as simple as telling
the computer to "draw a graph," as we did in the animation above. However, to compute the sign matrix
for the set of polynomials
$p_2, \dots, p_n$ $p_1'$ - The remainder
$r_1$ that results upon dividing$p_1$ by$p_1'$ - The remainders
$r_2, \dots, r_n$ that result upon dividing$p_1$ by$p_i$ , for$2 \leq i \leq n$
We recursively compute the sign matrix for this set, and use it to construct the sign matrix for
Including
The reason for including
So, if
The remainders are included to help us deduce the sign of
for every
But why bother with division for this? Indeed, if we added the polynomials
-
$p_1'$ has smaller degree than$p_1$ - Because
$r_1$ is the remainder upon dividing$p_1$ by$p_1'$ ,$r_1$ has smaller degree than$p_1'$ (which already has smaller degree than$p_1$ ). - Because
$r_i$ is the remainder upon dividing$p_1$ by$p_i$ ,$r_i$ has smaller degree than$p_i$ .
Now, if we choose
- Because
$r_i$ is the remainder upon dividing$p_1$ by$p_i$ ,$r_i$ has smaller degree than$\mathbf{p_1}$ .
Then, when we combine all the observations, we see that the input to the recursive call is constructed by removing
a polynomial
- Either we degree the maximum degree of the input set,
- Or we decrease the number of polynomials in the input having the maximum degree by
$1$
Any recursion having this property will terminate; indeed, the property implies that the set of pairs (# of polys with max degree, max degree)
form a strictly decreasing sequence in
The base case of the recursion can be taken to be any set of constant polynomials. Constant polynomials
never change their sign, so the sign matrix in this case has just one column:
Almost all the pieces are in place. We've described how to construct the input to the recursive call,
we've provided a base case, and we've argued that the recursion terminates. We now explain how to
construct the solution given the output of the recursive call. Most of the work was already done in
explaining the intuition behind choosing the input to the recursive call. The signs of
We start by going through the roots
The column DETERMINING SIGN of p_1 AT ROOTS of p_1', p_2, ..., p_n in the implementation.
From this point onwards, we don't need any of the information from the remainder polynomials, so we
implement a step which removes them from the sign matrix. In doing so, we may end up with "root" columns
which are no longer the root of any polynomial in the matrix (e.g. REMOVING REMAINDER INFORMATION AND FIXING MATRIX
in the provided implementation.
Since all roots where we didn't deduce a sign for
- Is there a root of
$p_1$ in the interval? Since$p_1'$ is part of the recursively computed sign matrix, it has no roots in any of the intervals. Our reasoning from before then shows that$p_1$ can have at most one root in each interval, so we only have to see if the number of roots is$0$ or$1$ . - What is the sign of
$p_1$ on the interval? Or, if there is a root of$p_1$ in the interval that splits the interval into pieces, what is the sign of each polynomial on the new pieces?
There are four different types of intervals: the bi-infinite interval
These signs can be computed from the sign of the derivative
- If
$p_1'$ is positive as$x \to \infty$ (i.e. if in the interval of the sign matrix which touches$\infty$ , the sign of$p_1'$ is positive), then$p_1$ must eventually increase forever. Since we're dealing with polynomials which cannot have asymptotic growth, increasing forever implies increasing without bound. This means that$p_1$ must be positive as$x \to \infty$ , so we say that the sign of$p_1$ at the pseudo-endpoint$\infty$ is positive. Likewise, if$p_1'$ is negative as$x \to \infty$ , the sign of$p_1$ at$\infty$ is negative. - If
$p_1'$ is positive as$x \to -\infty$ (i.e. if in the interval of the sign matrix which touches$-\infty$ , the sign of$p_1'$ is positive), then if we go far enough to the left (past all the roots of$p_1$ ),$p_1$ must eventually be decreasing as we move further left towards$-\infty$ (the sign of the derivative gets flipped here because$x$ is moving to the left, in the opposite of the canonical direction). So we say that the sign of$p_1$ at the pseudo-endpoint$-\infty$ is negative. Likewise, if$p_1'$ is negative as$x \to \infty$ , we say that the sign of$p_1$ at$-\infty$ is positive. - The last case, where
$p_1'$ is zero on an interval, implies that$p_1'$ is the zero polynomial. This means that$p_1$ , a polynomial of maximal degree, is constant. But then we wouldn't have entered this recursive case anyway - sets of constant polynomials are handled by the base case. So this case is impossible here.
Now that we have assigned pseudo-endpoints to the unbounded intervals, along with the sign of
Given an interval with endpoints
- If
$p_1(a)$ is positive while$p_1(b)$ is negative, the intermediate value theorem guarantees the existence of a root of$p_1$ in$(a, b)$ . Note that$(a, b)$ is an interval in a sign matrix containing$p_1'$ , so$p_1'$ has no roots in$(a, b)$ . By our earlier observation, this means that$p_1$ has at most one root in$(a, b)$ . Thus,$p_1$ has exactly one root in$(a, b)$ , which we call$c$ . This root splits the interval$(a, b)$ into three pieces:$(a, c)$ ,$c$ , and$(c, b)$ . The signs of all polynomials on these new pieces are determined as follows:- All the "old" polynomials
$p_1', p_2, \dots, p_n$ have invariant sign on the entire interval$(a, b)$ , so they certainly have invariant sign on the subsets of this interval$(a, c)$ ,$c$ , and$(c, b)$ . Thus, the sign of any of these polynomials on the pieces$(a, c)$ ,$c$ , and$(c, b)$ is the same as their sign on$(a, b)$ , and we can just copy the sign over. -
$p_1$ itself is certainly$0$ at$c$ , since$c$ is a root of$p_1$ . Moreover, we know that$c$ is the only root of$p_1$ in$(a, b)$ . This implies that there are no roots of$p_1$ in$(a, c)$ or in$(c, b)$ . Since$p_1$ is positive at$a$ , we have that$p_1$ must be positive over all of$(a, c)$ ; else there would be a root of$p_1$ in$(a, c)$ . Likewise, since$p_1$ is negative at$b$ ,$p_1$ must be negative over all of$(c, b)$ .
- All the "old" polynomials
- The case where
$p_1(a)$ is negative while$p_1(b)$ is positive is exactly dual to this. Again, there must exist a root$c$ of$p_1$ in$(a, b)$ , splitting the interval into three pieces$(a, c)$ ,$c$ , and$(c, b)$ . The signs of$p_1', p_2, \dots, p_n$ on the new pieces are copied from their signs on$(a, b)$ , and the signs of$p_1$ on$(a, c)$ ,$c$ , and$(c, b)$ are negative, zero, and positive respectively. - If
$p_1(a)$ is positive while$p_1(b)$ is either positive or zero, we look at the sign of$p_1'$ on the interval. As discussed before,$p_1$ is not constant, so$p_1'$ is either positive or negative on the interval, implying that$p_1$ is either strictly increasing or strictly decreasing on the interval. In the increasing case, since$x \in (a, b) \implies x > a$ , we have$p_1(x) > p_1(a) > 0$ for all$x \in (a, b)$ , and so$p_1$ has no roots in$(a, b)$ and is positive on$(a, b)$ . In the decreasing case, since$x \in (a, b) \implies x < b$ , we have$p_1(x) > p_1(b) \geq 0$ for all$x \in (a, b)$ , and so again$p_1$ has no roots in$(a, b)$ and is positive on$(a, b)$ . - The cases where
$p_1(a)$ is negative and$p_1(b)$ is either negative or zero is identical to the above. - The cases where we instead fix
$p_1(b)$ to be positive (resp. negative) and allow$p_1(a)$ to be either positive or zero (resp. negative or zero) use the same reasoning as well. - The last remaining case is of
$p_1(a) = p_1(b) = 0$ . This is impossible, as the mean value theorem would guarantee the existence of$c \in (a, b)$ with$p_1'(c) = (p_1(b) - p_1(a)) / (b - a) = 0$ , i.e. a root of$p_1'$ in$(a, b)$ , which we know cannot exist.
Combining the above reasoning with our pre-processing, we've described now how to compute the signs of FILTERING AND MERGING RESULT,
and is very similar to the previous REMOVING REMAINDER INFORMATION AND FIXING MATRIX step), and we output the result.
We provide an implementation of sign matrix computation via this algorithm here. The implementation is also embedded below. Enter a comma separated list of polynomials with rational coefficients. Enough output will be produced to trace the steps described above. Here are a couple of example inputs:
-
$p_1(x) = x^2 -1, p_2(x) = x^2 + 2x$ can be entered asx^2 - 1, x^2 + 2x. -
$p_1(x) = 4x^2 - 4$ ,$p_2(x) = (x + 1)^3$ , and$p_3(x) = -5x + 5$ can be entered as4x^2 - 4, x^3 + 3x^2 + 3x + 1, -5x + 5. Expansion must be done manually as the parser is not intelligent. - If you try your own inputs, please keep them relatively small in length and degree. The implementation is not optimized in the least, and will likely not run in time for large inputs.
Some notes on how to read the output:
- In the presentation above, the columns of the sign matrix were the roots and intervals, while the rows were the polynomials.
In the implementation it was more convenient to print the roots and intervals as the rows and have an association list mapping
each polynomial in each row to its sign. For example, here are the "mathematical" and "program" notations for the sign matrix
of
${p_1, p_2, p_3}$ , where$p_1(x) = 4x^2 - 4$ ,$p_2(x) = (x + 1)^3$ , and$p_3(x) = -5x + 5$ . \[ \begin{array}{cccccc} & (-\infty, x_1) & x_1 & (x_1, x_2) & x_2 & (x_2, \infty) \\\ p_1 & + & 0 & - & 0 & + \\\ p_2 & - & 0 & + & + & + \\\ p_3 & + & + & + & 0 & - \end{array} \]
neginf: (4x^2-4, +), (x^3+3x^2+3x+1, -), (-5x+5, +),
root: (4x^2-4, 0), (x^3+3x^2+3x+1, 0), (-5x+5, +),
interval: (4x^2-4, -), (x^3+3x^2+3x+1, +), (-5x+5, +),
root: (4x^2-4, 0), (x^3+3x^2+3x+1, +), (-5x+5, 0),
posinf: (4x^2-4, +), (x^3+3x^2+3x+1, +), (-5x+5, -),
- The row type
infdenotes the bi-infinite interval$(-\infty, \infty)$ . - In the
DETERMINING SIGN ON INTERVALSstep, each time we encounter an interval, we print the "context" - this is the information about the signs on the surrounding roots that is necessary to process the interval. - For reasons unknown, the creator of the polynomial/rational number library we use decided to output rational
numbers in decimal format. For example, the rational number
$1/6$ is outputted as0.1(6). Writing rational numbers in fraction form for the input seems to work though.
Comma-separated list of polynomials
Print recursively (may generate a lot of output)
Run
This script combines the sign matrix computation with the decision procedure we described earlier to evaluate the truth of univariate formulae. The input format is not very polished; here are a couple of examples and notes on the input format
- The formula
$\exists x, (x^2 + 1 \leq 0)$ is inputted as follows:
{
quantifier: "exists",
formula: {
poly: new Polynomial("x^2+1"),
signs: ["-", "0"]
}
}
- The formula
$\forall x, (\neg(x^2 + 1 \leq 0))$ is inputted as follows:
{
quantifier: "forall",
formula: {
conn: "not",
sf: {
poly: new Polynomial("x^2+1"),
signs: ["-", "0"]
}
}
}
- The formula
$\forall x, [(p_1(x) \geq 0 \wedge p_2(x) \geq 0) \vee (p_3(x) > 0)]$ , where$p_1(x) = 4x^2 - 4$ ,$p_2(x) = (x + 1)^3$ , and$p_3(x) = -5x + 5$ , is inputted as follows:
{
quantifier: "forall",
formula: {
conn: "or",
sf: [
{
conn: "and",
sf: [
{
poly: new Polynomial("4x^2-4"),
signs: ["+", "0"]
},
{
poly: new Polynomial("x^3+3x^2+3x+1"),
signs: ["+", "0"]
}
]
},
{
poly: new Polynomial("-5x+5"),
signs: ["+"]
}
]
}
}
- When inputting polynomials, it is important to use the variable
xand to leave no white space in the string. - Again, please stick to small examples (in terms of number of unique polynomials and their degree) in order for the computation to complete quickly.
Run
This webpage embeds all the deliverables, and serves as both the final project and the term paper. We summarize all the deliverables briefly here:
- The written explanation given in the above sections.
- The animations that complement the explanations, embedded in the appropriate locations. The code used to generate the animations and high-quality renders of the animations are also available; see here
- The implementations of the sign matrix calculation and univariate decision procedure are available here.
Thanks to Prof. Platzer for making me aware of the Cohen-Hörmander algorithm and providing resources to learn more about it.
Thanks to Grant Sanderson (3blue1brown) and the other manim developers for the framework, without which creating the animations seen above would have been impossible.
Thanks to Robert Eisele (GitHub user infusion) for the Fraction.js rational number library and the Polynomial.js univariate polynomial library.