# Which data structure to use to solve equations?

Let's say I have two equations for a geometric object (a rectangle):

$\left\{ \begin{array}{l} x \ge 0 \\ y \ge 0 \\ A \ge 0 \\ P \ge 0 \\ A = x*y \\ P = 2*x + 2*y \end{array} \right.$

Now I would like to set the values of some variables, and compute the values of the other ones.

At first, I thought it could be done with simple tree transformations on the AST, but I realized that I was more or less building a full solver.

The easy case is when one side of the equality consists only of constant terms, as it consists of a simple recursive tree evaluation. But the interesting case is when there are several variable bound together in a subtree. For example, for $(A,P)=(30,11)$ the solution is $(x,y)=(5,6)$ or $(6,5)$

Can arbitrary complex instances of this problem solved just with tree transformations ? Or do I need something more elaborate ?

• To solve linear systems like this uses matrices (and perhaps linear programming). The case of quadratic (or in general polynomial) systems of equations is much harder. I don't think you get a simple answer. Use a computer algebra system? Ask over at math.stackexchange.com about techniques to solve this type of systems of equations? – vonbrand Apr 7 '13 at 19:44
• Are you looking to write something or do you just want to solve this sort of problem? If the latter, then try an SMT solver. Z3 is one example; it's free and open-source so you might learn something about attacking this sort of problem by examining its source code. – Kyle Jones Apr 7 '13 at 21:20

Which data structure to use to solve equations?

I assume that we are talking about polynomial equations here, because I want to base my answer on a Gröbner basis method. The example problem given in the question contains equations, inequalities, variables and parameters. We can replace inequalities $f(x_1,\ldots,x_n) \geq 0$ by equations with slack variables $f(x_1,\ldots,x_n) - e^2 = 0$. We can convert parameters $P$ into variables by adding equations enforcing their value $P-11 = 0$. We can convert equations $f(x_1,\ldots,x_n)=g(x_1,\ldots,x_n)$ to the form $f(x_1,\ldots,x_n)-g(x_1,\ldots,x_n) = 0$.

$\left\{ \begin{array}{r} A - 30 = 0\\ P - 11 = 0\\ x - e_x^2 = 0\\ y - e_y^2 = 0\\ A - e_A^2 = 0\\ P - e_P^2 = 0\\ x*y - A = 0 \\ 2*x + 2*y - P = 0 \end{array} \right.$

Each equation can be represented as a list of its monomials and their coefficients. A monomial $M$ is a product $M=x_1^{a_1}\cdots x_n^{a_n}$, where the $a_i$ are nonnegative integers. The vector $[a_1, \ldots ,a_n]$ is called the exponent vector of $M$. A monomial can be represented by its exponent vector. This statement and the preceding description assumes that we have already fixed an enumeration of the variables. There are certainly more efficient ways to store a single equation, but these are just optimizations.

Can arbitrary complex instances of this problem solved just with tree transformations ?

It depends a bit on what we consider to be a tree transformation, and whether we mean by solved just the symbolic precomputation, or the entire solution process. A Gröbner basis method is able to transform the intial system into a triangular form, with respect to any admissible monomial order. In our case, we will probably want to use an elimination order such that the parameters are smaller than all other monomials, and the slack variables are bigger than all other monomials. The example system will then be transformed to

$\left\{ \begin{array}{r} e_x^2 - x = 0 \\ e_y^2 - y = 0 \\ e_A^2 - A = 0 \\ e_P^2 - P = 0 \\ x + y - \frac{1}{2}P = 0 \\ y*y - \frac{1}{2}P*y + A = 0 \\ A - 30 = 0 \\ P - 11 = 0 \\ \end{array} \right.$

This would be the results of the preprocessing. We can now solve this system numerically by back substitution. So we set the parameters to $P=11$ and $A=30$, solve the quadratic equation for $y$, then get $x$ from the linear equation, and then check the positivity constraints for $P$, $A$, $y$ and $x$. There are two solutions for $y$, so we have to do part of this process twice.

The are only two operations used during the transformation. The first operation is to reduce a polynomial by other polynomials, that is $f$ is replaced by $f-\sum_{g\in G}q_gg$ where the $q_g$ are polynomials. (This should actually be a multivariate division with respect to a monomial order.) The other operation is to add the S-polynomial of two polynomials. The S-polynomial of two polynomials $f$ and $g$ with leading terms $l_f$ and $l_g$ (with respect to a monomial order) is given by $\frac{a}{l_f}f-\frac{a}{l_g}g$, where $a$ is the least common multiple of the leading terms $l_f$ and $l_g$. In a certain sense, one can really call these operations tree transformations, because we just transform polynomials into other polynomials.

One interesting lesson for myself here was that adding the inequalities and parameters didn't really affect the outcome of the Gröbner basis method. We had to choose an appropriate elimination order for the monomial order. However, there are well know fill-in reduction strategies, and any sufficiently complete strategy would have automatically chosen an appropriate order. (The equation system can be seen as a bipartite graph between the variables and the equations, where a variable is connected to an equation if it occurs in the equation. This bipartite graph already contains enough information to determine an appropriate elimination order such that any harm which could be caused by the inequalities and parameters is avoided.)

Or do I need something more elaborate ?

One major problem is the order in which to apply the transformations. This already starts with the determination of an appropriate monomial order, but doesn't end there. Some of the methods developed for the direct solution of sparse linear systems can also be applied here. In addition, the reduction to an univariate polynomial of high degree is also a bit dangerous, from a numerical robustness point of view.

Even so there might be no general (efficient, reliable, numerically robust, consistent) solution, problems like the one given in the question arise for many applications. The efficiency, robustness and consistency of the solution methods may have improved significantly over time, but even the imperfect, slow and inconsistent initial solution approaches were good enough for many "real world" use cases. I find the dissertation from 1999 about "reinventing" the kernel of the interactive geometry software Cinderella quite illuminating in this respect.

The methods used in the publication above are limited to geometrical constructions, but they remind me of homotopy continuation methods. This overview presentation also refers to some symbolic solution methods:

and the well known Newton-Raphson method. The wikipedia article on system of polynomial equations paints a similar picture.

Kyle Jones mentioned in a comment Z3 as a SMT solver that might be able to solve this kind of problem. Frankly, I don't have practical experience with any of these methods, and never even heard of Wu's method or SMT solvers before. But there exist solvers or libraries for nearly all the mentioned solution approaches, like PHCpack (written in Ada) for homotopy continuation methods and FGb for Gröbner basis computations.

I may later add another answer about how to specifically apply a Gröbner basis method to the example problem, but let me briefly show just the solution, to give a feeling how it works:

The system $\left\{ \begin{array}{l} A = x*y \\ P = 2*x + 2*y \end{array} \right.$ gets reduced to $\left\{ \begin{array}{r} x + y - \frac{1}{2}P = 0 \\ y*y - \frac{1}{2}P*y + A = 0 \\ \end{array} \right.$. This system is in triangular form, and can be solved by back substitution. One obvious problem is that you still have to solve polynomial equations in a single variable. And already the two different solutions of quadratic equations can quickly lead to an explosion of the runtime, when you have a triangular system of $n$ quadratic equations. (Each equation has two different solutions as a function of the previously chosen solutions, so you have to try $2^n$ solutions in the worst case.)