First, let’s motivate the study of Gröbner bases through an example of solving a system of linear equations. Consider the following set of linear equations

which can be written in the matrix form

By performing Gaussian Elimination, we arrive at the reduced row echelon form

You can check this yourself by running the following Python code using SymPy:

import sympy as sp
def gaussian_elimination(M):
    return M.rref()[0]

gaussian_elimination(sp.Matrix([[2, 3, -1, 0],
                                [1, 1, 0, 1],
                                [1, 0, 1, 3]]))

The results should look like this:

Matrix([[1, 0, 1, 3], [0, 1, -1, -2], [0, 0, 0, 0]])

From the above, we can see that the solution to the system is

which is equivalent to a line in \(\mathbb{R}^{3}\).

Now, can we do something similar to a system of polynomial equations? In other words, we want to develop a method to determine the solution to

The method of Gröbner bases allows us to do so.

I will start with a concrete example of using Gröbner Bases to solve a system of polynomial equations. Words and phrases in bold are new concepts which I will define later.

Consider the following system of equations:

We first define the ideal

where the angle brackets represent a set of polynomials following

Right now treat this ideal roughly as the set of polynomials of which the solutions to \(x\), \(y\) and \(z\) satisfies the original polynomial system. The Gröbner bases for this ideal \(I\) with respect to lex order is the following:

Note that \(g(4)\) contains only \(z\), and

which becomes zero when \(z = 0, 1, -1 + \sqrt{2}, -1 - \sqrt{2}\). Substituting them back to \(g_{2} = 0\) and \(g_{3} = 0\), we can determine \(y = 1, 0, -1 + \sqrt{2}, -1 - \sqrt{2}\). Then using \(g_{1}\), we can get the corresponding values for \(x\), which leads to the following five possible solutions to the original system:

Notice how the existence of \(g(4)\), a polynomial with only one variable \(z\), enables us to solve for all variables through back-substitution. To understand how this process works, we need to answer the following two questions:

1. What is the relationship of the ideal with respect to the original system of equations?
2. What are Gröbner bases?
3. How do Gröbner bases help us solve polynomial systems?

Now we will make an attempt to answer these three questions.

Unfortunately we need to go over a few definitions so that we have a common ground for discussing Gröbner bases. Concepts like polynomials should be familiar, but we need some more definitions in order to discuss them in a precise manner. We will see relevant SymPy code as well to play around with some of the concepts.

We now try to answer Question 1 we proposed at the end of Section A Taste of Gröbner Bases. In particular, we define what exactly are the solutions to a system of polynomials, and how ideals relate to them.

Gröbner bases, introduced by B. Buchberger, are the tools we need. In A Theoretical Basis For the Reduction Of Polynomials To Canonical Forms (1976)⁵, Buchberger says the following regarding Gröbner bases:

Our algorithm proceeds by constructing a new basis for a given ideal from which the answer to the computability and decidability problems may be easily read off.

We now discuss the Gröbner bases in brief (answering Question 2), and learn how they connect to solving systems of polynomial equations.

The definition of Gröbner bases is the following: given a monomial order for a polynomial ring \(k[x_{1}, \ldots, x_{n}]\), a finite non-zero subset \(G\) of the ideal \(I \subseteq k[x_{1}, \ldots, x_{n}]\) is a Gröbner basis if

where \(\text{LT}(I)\) is the set of leading terms of nonzero elements in \(I\) and defined as

To check whether a set is a Gröbner basis, one way is to see whether the definition above holds. Consider the ideal \(J = \langle x+z, y-z \rangle\) with bases \(x+z\) and \(y-z\). We claim that \(x+z\) and \(y-z\) constitute a Gröbner basis using lex order. Then \(\langle \text{LT}(x+z), \text{LT}(y-z) \rangle = \langle x, y \rangle\).

We now need to show leading terms of all nonzero elements of \(\langle x+z, y-z \rangle\) lie within \(\langle x, y \rangle\), which is equivalent to showing they are all divisible by either \(x\) or \(y\). We can prove this by contradiction. Assume we have an \(f = A(x+z) + B(y-z) \in J\), and that \(f\) is nonzero and \(\text{LT}(f)\) is divisible by neither \(x\) nor \(y\). Hence \(f\) is a polynomial in \(z\) only. Since \(f \in J\), \(f\) is zero on all points in \(\mathbf{V}(x+z, y-z)\). Note that \((-t, t, t)\) is in \(\mathbf{V}(x+z, y-z)\), and \(f\), a polynomial in \(z\) alone, has to be the zero polynomial, which is a contradiction. Hence \(x+z\) and \(y-z\) form a Gröbner basis for \(J\).

The “proper” and more mechanical way to check whether a set of polynomials is a Gröbner basis (which also forms the basis of the algorithm to construct Gröbner basis) is to use the so-called Buchberger’s Criterion:

There are two concepts in this definition that we haven’t discussed yet. The first one is the division of a polynomial by \(G\). The actual algorithm³ matters less in this case; we only need to know that there exists an algorithm that can rewritten a polynomial \(f\) in \(k[x_{1}, \ldots, x_{n}]\) as

where \((g_{1}, \ldots, g_{s})\) is an ordered set of polynomials, \(q_i \in k[x_{1}, \ldots, x_{n}]\), and either \(r = 0\) or \(r\) is a polynomial which has no monomial that is divisible by any of \(\text{LT}(g_{1}), \ldots, \text{LT}(g_{s})\).

\(S(f, g)\) refers to S-polynomials, which is defined as

where \(x^{\gamma}\) is the least common multiple of \(\text{LM}(f)\) and \(\text{LM}(g)\) (a monomial with power of each variable equal to the largest power of the corresponding variable in \(\text{LM}(f)\) and \(\text{LM}(g)\)). We can define s_poly(f, g) using SymPy to calculate S-polynomials:

import sympy as sp

def s_poly(f, g, *gens):
    """Calculate the S-polynomial for f and g.
    Note that this uses the default lex order.
    """
    lcm = sp.lcm(sp.LM(f, *gens), sp.LM(g, *gens))
    s = sp.simplify(lcm * (f / sp.LT(f, *gens) - g / sp.LT(g, *gens)))
    return s

The definition of S-polynomials may seem arbitrary. One intuition behind it is that S-polynomials cancel out the leading terms. Take \(\{ xy + 2x - z, x^{2} + 2y - z\}\) as an example. In this case, \(x^{\gamma} = x^{2} y\), and we have

Notice how the leading terms of \(xy + 2x - z\) and \(x^{2} + 2y - z\) got canceled out. The precise way to describe such cancellation is that \(\text{multideg}(S(f, g)) \leq \gamma\) (where \(\gamma\) follows the definition in Theorem 4.)

Now we can test out the Buchberger’s Criterion. Again use \(\{ xy + 2x - z, x^{2} + 2y - z\}\) as an example. We have \(S( xy + 2x - z, x^{2} + 2y - z ) = 2x^{2} - xz - 2y^{2} + yz\). We can then check the remainder of \(2x^{2} - xz - 2y^{2} + yz\) with division by \(\{ xy + 2x - z, x^{2} + 2y - z \}\), which is \(-x*z - 2*y**2 + y*z - 4*y + 2*z\), hence the two polynomials do not form a Gröbner basis.

Buchberger’s Criterion leads nicely to the algorithm that we can use to actually construct Gröbner bases. Given a list of polynomials \(f_{1}, \ldots, f_{s}\) and let \(I = \langle f_{1}, \ldots, f_{s} \rangle\), the Buchberger’s Algorithm constructs a Gröbner basis \(G = (g_{1}, \ldots, g_{t})\) for \(I\). Below is a crude implementation in Python.

import copy
import itertools

def buchberger(F, *gens):
    """Buchberger's Algorithm
    Note that this is slightly different from the pesudocode
    provided in Cox et. al. 2015.
    """
    G = copy.deepcopy(F)
    pqs = set(itertools.combinations(G, 2))
    while pqs:
        p, q = pqs.pop()
        s = s_poly(p, q, *gens)
        _, h = sp.reduced(s, G, *gens)
        if h != 0:
            for g in G:
                pqs.add((g, h))
            G.append(h)
    return G

The s_poly function is used to generate S-polynomials. Intuitively, we are extending the original basis repeatedly with the nonzero remainders. For a complete proof of this algorithm, please refer to Chapter 2 of Ideals, Varieties, and Algorithms ³.

We can check this implementation against Sympy:

f1 = x**3 - 2*x*y
f2 = x**2*y - 2*y**2 +x
G = [f1, f2]

S = s_poly(f1, f2)
print(f"S Poly: = {S}")

a, r = sp.reduced(s, G)
print(f"Remainder: {r}")

# test Buchberger
G_basis = buchberger(G)
print(f"Groebner basis: {G_basis}")

G_basis_sympy = sp.groebner(G, x, y, z, order='lex')
print(f"Groebner basis (SymPy): {G_basis_sympy}")

Running this we will see something like this:

S Poly: = x**3*y*(-(x**2*y + x - 2*y**2)/(x**2*y) + (x**3 - 2*x*y)/x**3)
Remainder: -x**2
Groebner's basis:
[x**3 - 2*x*y,
 x**2*y + x - 2*y**2,
 -x**2,
 x - 2*y**2,
 -4*y**3]
Groebner's basis (SymPy):
GroebnerBasis(
[x - 2*y**2,
 y**3], x, y, z, domain='ZZ', order='lex')

It seems like the basis returned by SymPy is a subset of our basis (with a factor of \(-4\) in front of \(y^3\)). It turns out that the groebner(G) function in SymPy is implemented so that it returns the reduced Gröbner basis, whereas our implementation only returns one Gröbner basis. For a nonzero polynomial ideal, the reduced Gröbner basis is unique, but there might be infinitely many Gröbner bases. We define reduced Gröbner bases as the following:

A common way to obtain a reduced Gröbner basis is to (1) obtain a minimal basis by ensuring leading monomials of the elements do not divide each other and (2) replace each element in the basis with the remainder of its reduction by other elements in the basis, and then divide each element by the coefficient of its leading term. One may also define the reduced Gröbner basis without the condition on its leading coefficient, and the uniqueness of the basis is therefore up to a multiplicative factor. Here’s the final function that we can use to compute a reduced Gröbner basis:

def groebner(F, *gens):
    """Calculate a reduced Groebner basis for F.
    Use the default lex order.
    """
    F_polys, opt = sp.parallel_poly_from_expr(F, *gens)
    domain = sp.EX
    ring = sp.polys.rings.PolyRing(gens, domain=domain)

    # buchberger
    G = buchberger(F_polys, *gens)

    # minimal
    temp = copy.deepcopy(G)
    G_minimal = []
    while temp:
        f0 = temp.pop()
        if not any(sp.polys.monomials.monomial_divides(f.LM(), f0.LM()) for f in temp + G_minimal):
            G_minimal.append(f0)

    # reduce
    G_reduced = []
    for i, g in enumerate(G_minimal):
        _, remainder = sp.reduced(g, G_reduced[:i] + G_minimal[i+1:])
        if remainder != 0:
            G_reduced.append(remainder)

    # sort
    # courtesy of SymPy buchberger implementation
    polys, opt = sp.parallel_poly_from_expr(G_reduced, *gens)
    polys = [ring.from_dict(poly.rep.to_dict()) for poly in polys if poly]
    G_reduced = sorted(polys, key=lambda f: f.LM, reverse=True)
    return sp.parallel_poly_from_expr([x.monic().as_expr() for x in G_reduced], *gens)[0]

Let’s test it out against our previous example.

G_basis_reduced = groebner(G, x, y, z)
print(f"Reduced Groebner basis: {G_basis_reduced}")

You should see outputs similar to this:

Reduced Groebner basis: [Poly(x - 2*y**2, x, y, z, domain='ZZ'),
                         Poly(y**3, x, y, z, domain='ZZ')]

Now we have a reduced Gröbner basis! Notice how there is an element in the basis that has only \(y\). It turns out this property of Gröbner bases will pave the way for us to solve systems of polynomial equations.

We are now ready to answer Question 3. We have learned about Gröbner basis, and understands that a reduced Gröbner basis uniquely determines an ideal. We have established that solutions of a polynomial system correspond to variety of the ideal generated by the system (Proposition 3). In this section, we discuss the elimination property of Gröbner basis, and use it to solve systems of polynomial equations.

Let’s go back to the original system at the beginning,

We can find its reduced Gröbner basis using the tools we have developed.

F = [x**2 + y + z - 1, x + y**2 + z - 1, x + y + z**2 - 1]
G = groebner(F, x, y, z)
print(f"Reduced Groebner basis: {G}")

You should see outputs like the ones listed below:

Reduced Groebner basis:
[Poly(x + y + z**2 - 1, x, y, z, domain='QQ'),
 Poly(y**2 - y - z**2 + z, x, y, z, domain='QQ'),
 Poly(y*z**2 + 1/2*z**4 - 1/2*z**2, x, y, z, domain='QQ'),
 Poly(z**6 - 4*z**4 + 4*z**3 - z**2, x, y, z, domain='QQ')]

The last element of the reduced Gröbner basis is \(z^{6} - z^{4} + 4z^{3} - 2z^{2} + 4z\), which is in \(z\) only. In A Taste of Gröbner Bases, we then proceed to solve for \(z\) and back substitute the value of \(z\) to solve for \(x\) and \(y\).

Roughly, we can divide this process into three steps:

1. Find the reduced Gröbner basis, locate the element that has only one variable and solve for that variable.
2. Substitute that variable into the rest of the basis.
3. Repeat the previous two steps with the new basis.

It turns out that this elimination of variables under lex order is a property of Gröbner bases. This result is the so-called Elimination Theorem:

\(I_l\) is defined as \(I \cap k\left[x_{l+1}, \ldots, x_n\right]\). The use of \(\cap k\left[x_{l+1}, \ldots, x_n\right]\) (intersection) may seem a bit confusing, but it is essentially a formal way to say that we have eliminated \(x_{1}, \ldots, x_{l}\). This theorem essentially allows us to recursively eliminates variables one by one.

The final piece that we need is the Extension Theorem, which gives us a way to understand when can we successfully extend a partial solution to a full solution.

The part about \(\left(a_2, \ldots, a_n\right) \notin \mathbf{V}\left(c_1, \ldots, c_s\right)\) essentially states that a partial solution can be extended if it does not cause the leading coefficients (\(c_{i}\)) go to zero. Proofs for Theorem 6 and Theorem 7 can be found in Chapter 3 of Cox et. al³.

We are now ready to implement a solver for system of polynomials. First, we define two helper functions (is_univariate and subs_root) that determines whether a polynomial is in one variable only, and substitute in a solution for a variable.

def is_univariate(f):
    """Returns True if 'f' is univariate in its last variable.
    Based on SymPy solve_generic
    SymPy License: https://github.com/sympy/sympy/blob/master/LICENSE
    """
    for monom in f.monoms():
        if any(monom[:-1]):
            return False
    return True

def subs_root(f, gen, zero):
    """ Substitute in a solution for a generator
    Based on SymPy solve_generic
    SymPy License: https://github.com/sympy/sympy/blob/master/LICENSE
    """
    p = f.as_expr({gen: zero})

    if f.degree(gen) >= 2:
        p = p.expand(deep=False)

    return p

Then we can implement the solver that recursively eliminates variable and extends solution.

def solve_poly_system_recursive(F, gens, entry=False):
    """ Recursive helper function
    Based on SymPy solve_generic
    SymPy License: https://github.com/sympy/sympy/blob/master/LICENSE
    """
    basis = groebner(F, *gens)

    if len(basis) == 1 and basis[0].is_ground:
        if not entry:
            return []
        else:
            return None

    if len(basis) < len(gens):
        raise ValueError("System not zero-dimensional.")

    univar = [x for x in basis if is_univariate(x)]

    if len(univar) == 1:
        f = univar.pop()
    else:
        raise ValueError("System not zero-dimensional.")

    # find solutions for the univariate element
    gens = f.gens
    gen = gens[-1]

    zeros = list(sp.roots(f.ltrim(gen)).keys())

    if not zeros:
        # no solution
        return []

    # single variable
    if len(basis) == 1:
        return [(zero,) for zero in zeros]

    # recursively solve for the rest of the basis
    solutions = []
    for zero in zeros:
        new_system = []
        new_gens = gens[:-1]

        # back substitution of solution
        for b in basis[:-1]:
            eq = subs_root(b, gen, zero)

            if eq is not sp.core.S.Zero:
                new_system.append(eq)

        new_system = sp.parallel_poly_from_expr(new_system, *new_gens)[0]
        for solution in solve_poly_system_recursive(new_system, new_gens):
            solutions.append(solution + (zero,))

    if solutions and len(solutions[0]) != len(gens):
        raise ValueError("System not zero-dimensional.")

    return solutions


def solve_poly_system(F, *gens):
    """ Solve a system of polynomials with Groebner basis
    Based on SymPy solve_generic
    SymPy License: https://github.com/sympy/sympy/blob/master/LICENSE
    """
    result = solve_poly_system_recursive(F, gens, entry=True)
    return sorted(result , key=sp.default_sort_key)

The few ValueError("System not zero-dimensional.") represent various edge cases where the system has potentially infinite number of solutions. Please refer to Theorem 6, Chapter 5 of Cox et. al.³ for more details.

Let’s check our implementation against the official SymPy one:

F = [x**2 + y + z - 1, x + y**2 + z - 1, x + y + z**2 - 1]
F = list(map(lambda x : sp.Poly(x), F))

solution_sympy = sp.solve_poly_system(F, x, y, z)
print(f"Poly Sols (SymPy): {solution_sympy}")

solution = solve_poly_system(F, x, y, z)
print(f"Poly Sols: {solution}")

You should see outputs like these:

Poly Sols (SymPy): [(0, 0, 1),
(0, 1, 0), (1, 0, 0),
(-1 + sqrt(2), -1 + sqrt(2), -1 + sqrt(2)),
(-sqrt(2) - 1, -sqrt(2) - 1, -sqrt(2) - 1)]
Poly Sols: [(0, 0, 1),
(0, 1, 0), (1, 0, 0),
(-1 + sqrt(2), -1 + sqrt(2), -1 + sqrt(2)),
(-sqrt(2) - 1, -sqrt(2) - 1, -sqrt(2) - 1)]

Looks like we have the correct solutions!

In this post, we have focused on using Gröbner bases to solve polynomial systems. I hope you have learned something out of this long article.

Here are some additional resources if you want to dive deeper:

• You can use techniques related to Gröbner bases to prove geometric theorems automatically. See Chapter 6 of Cox et. al³.
• For a fast Gröbner basis library, you can use FGb.
• You can use Groebner.jl for calculating Gröbner bases in Julia.
• You can read about using reinforcement learning to select pairs of polynomials to compute S-polynomials in Buchberger’s algorithm here.
• You can read about strategies for selecting bases to speed up minimal solvers for computer vision here.