Using Gröbner bases for solving polynomial equations

Question

In my attempts to understand just how computer algebra systems "do things", I tried to dig around a bit on Gröbner bases, which are described almost everywhere as "a generalization of the Euclidean algorithm and Gaussian elimination". I've tried to look for examples of Gröbner bases in action, but have been unable to find any (that can be easily understood).

I could go ahead and just ask for an explanation with examples from people, but I'll go one step further. General plane conics can be represented by the Cartesian equation

$$ax^2+2bxy+cy^2+dx+fy+g=0$$

One common problem in dealing with conics is finding out if two conics intersect, and if so, where the intersection point(s) are. Usually one can do all this by eliminating variables accordingly.

How would, say, Buchberger's method, proceed on determining if two given conics intersect, and then find where they intersect?

Answer 1

In order to find the solutions of the system of conics you mention, it suffices to give a procedure to find the projection of the simultaneous vanishing set of the two conics onto some axis, for instance, the $y$-axis: the coordinates of the projection onto the $y$-axis are the $y$-coordinates of the intersection points. Knowing them allows you to substitute back into the equations the values of $y$ and solve a system of equations in a single variable $x$, thus making the problem simpler. In terms of ideals, suppose that we can find a non-zero element $r$ in the ideal generated by the two conics, that depends only on a single variable, say $y$. This means that every solution of the system has $y$-coordinate satisfying the polynomial $r$. Thus we have severely limited the choices for the $y$-coordinates of the intersection (namely, they all have to satisfy the polynomial $r$) and, provided we can actually solve the polynomial $r$ in one variable, we can then substitute the various values of $y$ we found back into the initial equations and solve for $x$, again using our algorithm for solving polynomials in a single variable. So far, so good, I hope! The question is how to produce the element $r$. This is where Gröbner bases come in.

In order to do the computation you ask, one would have to choose an appropriate monomial order. I will gloss over this, and simply "do the obvious". It is an exercise for you to figure out which order to use so that the computation I am going to make is actually a Gröbner bases computation. Scattered throughout the computation there will be also a couple of special cases that I will not deal with: again, you can treat those as exercises for you!

First, I am going to choose a generic basis, so that the first conic has the form $$ x^2 + \alpha x + \beta , $$ where $\alpha$ and $\beta$ are polynomials in $y$ only (I am trying to simplify the notation as much as possible; the only "assumption that I have made is that the coefficient of $x^2$ is non-zero, which can be arranged unless the "conic" is defined by a polynomial of degree at most one).

In this basis, the second equation can be assumed to have no $x^2$ term (indeed, eliminating the $x^2$ term using the first equation would be the first step in any reasonable Gröbner basis computation in which the term $x^2$ is the highest term in sight). Thus I am going to write the second conic as $$ \gamma x + \delta $$ where, as before, $\gamma$ and $\delta$ are polynomials in $y$ only. Again, just to fix on a definite case, I am going to assume that the polynomial $\gamma$ is non-zero. (If $\gamma$ were zero, then we would have found a polynomial in the ideal generated by the two conics which only depends on $y$: this was our goal at the start anyway! Obviously, you should worry about the case $\gamma=\delta=0$, but I won't.)

To compute a Gröbner basis, you would now compute S-polynomials: let's do it here as well. I am going to try to eliminate the $x^2$ term from the first equation, by using the second equation. This is easy: multiply the first equation by $\gamma$ and the second one by $x$ and take the difference: we are left with the equation $$ (\alpha \gamma - \delta) x + \beta \gamma . $$ Now we are going to eliminate the $x$ term from this last equation using again the second equation: multiply the second equation by $(\alpha \gamma - \delta)$, the last equation by $\gamma$ and subtract to obtain $$ \delta^2 - \alpha \gamma \delta + \beta \gamma^2 . $$ We found an expression independent of $x$!! We are done... provided this expression is not identically zero. You can figure out what this would mean and what happens in this case. Note also that the final expression is what we would have obtained if, at the very beginning, we had "solved" $x=-\frac{\delta}{\gamma}$ using the second equation, substituted in the first equation and cleared the denominators.

I hope that this "hybrid" computation explains what is going on: Gröbner bases and Buchberger's algorithm are a systematic way of "solving" systems of equations. You do not have to do any thinking, once you set up the problem. But you need to set up the problem so that it computes what you want. In this case, you could have used several shortcuts to get to the answer, without following all the steps. In more complicated situations, Buchberger's algorithm might be the best way of keeping track of all the steps to be taken.

Let me also comment that, except in the case in which you have a computer doing the computations for you, it is highly unlikely that Gröbner bases will help you with a specific question, unless you could have also found simple tricks to solve it right away.

Answer 2

You have to read through the book [David A. Cox, John B. Little and Don O'Shea, Ideals, Varieties, Algorithms].

In any case, if you start with a system of polynomial equations and compute a Groebner basis for the ideal they generate, you get a "maximally triangular" system of equations which is equivalent to the original one---that is why Groebner bases generalize Gaussian elimination.

Let me do a simple example using Macaulay 2. Consider the ring $\mathbb Q[x,y,z]$:

i1 : R = QQ[x, y, z, MonomialOrder => Lex]

o1 = R

o1 : PolynomialRing

the Fermat quintic surface $x^5+y^5+z^5=1$, whose ideal is

i2 : surface = ideal (x^5 + y^5 + z^5 - 1)

             5    5    5
o2 = ideal(x  + y  + z  - 1)

o2 : Ideal of R

and the curve $x^2=y^2$, $y^3+1=z^3$:

i3 : curve = ideal (x^2 - y^2, z^3 - y^3 - 1)

             2    2     3    3
o3 = ideal (x  - y , - y  + z  - 1)

o3 : Ideal of R

The ideal of the intersection of the surface and the curve is the sum of the ideals of the intersecands:

i4 : intersection = curve + surface

              2    2     3    3       5    5    5
o4 = ideal (x  - y , - y  + z  - 1, x  + y  + z  - 1)

o4 : Ideal of R

The number of points on the intersection, counting multiplicities, is the degree of the ideal:

i5 : degree intersection

o5 = 30

if we now compute the degree of the radical of the ideal, we get a lower number, so not all the points are simple:

i6 : degree radical intersection

o6 = 25

Now look at the lexicographic Groebner basis of the ideal of the intersection:

i7 : ideal gens gb intersection

              17      16      15      14      13      12       11       10      9      8       7       6      5      4      3  
o7 = ideal (9z   + 27z   + 54z   + 50z   + 15z   - 63z   - 104z   - 108z   - 35z  + 35z  + 108z  + 104z  + 63z  - 15z  - 50z  -
     ----------------------------------------------------------------------------------------------------------------------------
        2                 5            16      15      14      13      12       11      10      9       8      7       6       5
     54z  - 27z - 9, 20y*z  - 20y + 27z   + 18z   + 45z   - 75z   - 35z   - 164z   + 79z   - 10z  + 230z  - 10z  + 119z  - 164z 
     ----------------------------------------------------------------------------------------------------------------------------
          4       3      2              3    3                                       16        15        14        13        12  
     - 35z  - 155z  + 45z  + 18z + 67, y  - z  + 1, 50x*z - 50x + 50y*z - 50y + 1242z   + 2718z   + 4770z   + 2220z   - 1235z   -
     ----------------------------------------------------------------------------------------------------------------------------
          11        10        9        8        7        6        5        4        3        2                 2    2
     7979z   - 7481z   - 6010z  + 1810z  + 5240z  + 9394z  + 5346z  + 1240z  - 4030z  - 4005z  - 2657z - 583, x  - y )

o7 : Ideal of R

The first generator depends only on $z$. The second one and the third, only on $z$ and $y$, and the fourth depends (linearly) on $x$ too. This is a "triangular" system of equations, from which you can compute the 30 points of intersection, assuming you get to solve the first equation for $z$, which is of degree 17...

Answer 3

The example that you propose is a bit too general to serve as an introductory pedagogical example (compare e.g. the famous Kahan SIGSAM problem 9 of determining conditions for an ellipse to lie inside the unit circle). Instead, begin with one of the earliest historical examples: Gauss's proof that every symmetric polynomial can be written uniquely as a polynomial in the elementary symmetric polynomials. This is one of the first and the simplest applications of rewriting reduction via a lexicographic order. For a nice presentation see Chapter 7 of Cox, Little, O'Shea: Ideals, Varieties and Algorithms. They also present generalizations to the ring of invariants of a finite matrix group $\rm G \subset GL(n,k)$. You can find online copies of said book via various ebook databases.

Answer 4

This is a simplistic version of the Gröbner basis computation and complements Mariano's answer. The idea behind the Buchberger algorithm is to build a basis so that you are to be able to do "unique division" between nonlinear polynomials. This is done by computing (what is called the) S-polynomial between any 2 given polynomials until there are no more produced.

Decide some ordering (grlex, lex, revlex)
Initialize I = <f1, f2> 
where we use the polynomials you gave above.
f_1=ax^2+2bxy+cy^2+dx+fy+g
f_2=ax^2+2bxy+cy^2+dx+fy+g
Set n=2

Repeat the following in a loop
    Iteratively choose any two polynomials ( f_i, f_j ) from set I 
    n = n + 1
    Compute f[n] = S-polynomial( f_i, f_j )
Until all S-polynomials return 0

At the end of the computation, you will have a set I = {f1,f2, ..., fk} such that one of the polynomials will be in only one term (say 'x'). Solve that equation to get the one coordinate of the intersection, and then work upwards through the rest of the polynomials.

The general question is related to Elimination theory, and AFAIK there are three basic methods ( I can't speak of their complexity )

1. Gröbner basis, as seen above
2. Use Resultants (example below)
3. Wu's method of characteristic sets, which I can't say much about.

If you wanted to use resultants to find the intersection of say

E={f1(x, y,z), f2(x, y, z), f3(x, y, z)},

you would first compute

E1 = {g1 = Resultant_in_x (f1, f2), g2 = Resultant_in_x(f2, f3)} 

and then

E2 = { h1 = Resultant_in_y(g1, g2} }

The last polynomial h1 would be in one variable and amenable to solution.

Answer 5

At the OP's request here is a short description of how to use Groebner bases to generate thermodynamical identities. The basic setup for the thermodynamics of a Gibbsian substance is a pair of equations $T=f(p,V)$ and $S=g(p,V)$ which express the temperature and entropy as functions of $p$ and $V$ (the example that everybody knows is the ideal gas which, omitting physical constants, has $f(p,V) = pV$, $g(p,V)= \frac 1 {\gamma-1} (\ln p + \gamma \ln V)$). One can then express all thermodynamical quantities (e.g., the heat capacities--$C_p=\frac{fg_p}{f_p}$ and $C_V=\frac{fg_V}{f_V}$) in terms of the functions $f$, $g$ and their partials. This firstly provides a unified and systematic method to verify such identities but also one to generate them. The basic reason for this is that one can express a very large number of thermodynamic quantities (it runs into many tens of thousands) as simple algebraic expressions of a very few terms so that there are bound to be multiple relations between them (the law of small numbers). This is a situation which Groebner bases are tailor-made to handle. The details can be found in my article "A Systematic Approach to Thermodynamical Identities" (arXiv 1108.4760) which also contains a Mathematica programme which allows one to compute all thermodynamical quantities such as the the heat capacities at the click of a button (doing this by hand can be rather tedious for more elaborate models such as the van der Waals gas or the Feynman gas---it is usually only done in the secondary literature for very simple quantities and models---typically the ideal gas).

Answer 6

You might find Lecture 3 of James Carlson's CIMAT Lectures of value as well as other material on this page:

https://web.archive.org/web/20160909151728/http://www.math.utah.edu/~carlson/cimat/

(original link http://www.math.utah.edu/~carlson/cimat/ is broken.)

Answer 7

Here is a Mathematica example taken from Golden Fields: A case for the heptagon by Peter Steinbach

By applying Ptolemy's theorem to the diagonals ($a$ and $b$) of a regular heptagon one arrives at the equations:

$$ \begin{align} a^2-(1+b) &= 0, \\ a b-(a+b) &= 0, \\ b^2-(1+a b)&= 0 \end{align} $$

Applying GroebnerBasis gives:

GroebnerBasis[{a^2-(1+b),a b-(a+b),b^2-(1+a b)},{b,a}]

{1-2a-a^2+a^3,1-a^2+b}

That is,

$$ \begin{align} 1-2a-a^2+a^3 &= 0, \\ 1-a^2+b &= 0 \end{align} $$

The first of which can be solved to find a:

{1.80194,-1.24698,0.445042}

Whence the second equation can be used to find $b$.

Answer 8

Learning computation in M2, you find step-by-step examples in the first two books while the Cox book is more mathematically oriented. I provided below an example about maximal tringulization in terms of GR basis that generalises the Gaussian elimination as explained in the awesome answers such as Mariano. SandeepJ method outlining contains GR basis, Resultants and Wu's method. The Bill's notice "generalizations to the ring of invariants of a finite matrix group $G\subset GL(n,k)$" may motivate to investigate this further: How to analyse a sparse adjacency matrix?

Books

• Computations in algebraic geometry with Macaulay 2

• Macaulay 2 Tutorials

• Ideals, Varieties and Algorithms by Cox et all (more mathematical, no focus on M2)

Example with GR basis in M2

 i89 : R=QQ[x,y,z,MonomialOrder=>Lex]; I=ideal(x^2+y^2+z^2-1, x^2+z^2-y, x-z); transpose gens gb I


 o90 : Ideal of R

 o91 = {-4} | 4z4+2z2-1 |
       {-2} | y-2z2     |
       {-1} | x-z       |

               3       1
 o91 : Matrix R  <--- R

 i92 : degree radical I == degree I

 o92 = true

The GR basis method simplifies the original problem of polynomials to more easier polynomials. This has more focus on the mathematical structure. In doing this example, I found this procedure Solving equations useful and other answers above such as checking multiplicities, radical versus non-radical roots line.