# Gröbner bases are in ‘qspray’

**Saturn Elephant**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

There is something new in the **qspray** package:
*Gröbner bases*. Gröbner bases have many applications, and we
will see one of them now. My package **qspray** deals with
multivariate polynomials with rational coefficients.

Consider the following system of three polynomial equations:

\[\begin{align} x^2 + y + z & = 1 \\ x + y^2 + z & = 1 \\ x + y + z^2 & = 1 \end{align}\]

The goal is to find all solutions of this system, that is to say, the simultaneous roots of the three polynomials \(x^2 + y + z – 1\), \(x + y^2 + z – 1\), \(x + y + z^2 – 1\). Not easy, is it?

*Fact:* a Gröbner basis of (the ideal generated by) a list of
polynomials is a list of polynomials with the same simultaneous roots.

library(qspray) x <- qlone(1); y <- qlone(2); z <- qlone(3) f1 <- x^2 + y + z - 1 f2 <- x + y^2 + z - 1 f3 <- x + y + z^2 - 1 gb <- groebner(list(f1, f2, f3)) # Gröbner basis of f1, f2, f3 ( gbstrings <- lapply(gb, prettyQspray, vars = c("x", "y", "z")) ) ## [[1]] ## [1] "y + x + z^2 - 1" ## ## [[2]] ## [1] "((-9)*z^3)/2 - 2*z^5 - z^6 + 6*z^4 - y + z + z^7/2 + y^2" ## ## [[3]] ## [1] "y*z^2 - z^2/2 + z^4/2" ## ## [[4]] ## [1] "z^6 - z^2 + 4*z^3 - 4*z^4"

Do you see something nice? The fourth polynomial,
\(z^6 - z^2 + 4z^3 - 4z^4\), is
univariate: it depends only on the variable
\(z\). We could numerically solve it
with the R function `polyroot`

, but we want exact roots. It
is easy to see that \(0\) and
\(1\) are two roots of this polynomial.
So we can divide this polynomial by
\(z(z-1)\), and it will remain a factor
of degree \(4\). But I know, you are
lazy. Me too. That’s why I prefer to use the
**Ryacas** package to perform this task:

library(Ryacas) yac_str("Factor(z^6 - z^2 + 4*z^3 - 4*z^4)") ## [1] "(2*z+z^2-1)*(z-1)^2*z^2"

So it remains to find the roots of
\(2z + z^2 - 1\). We use
**Ryacas** again, as we’re lazy:

yac_str("PSolve(2*z + z^2 - 1, z)") ## [1] "{(Sqrt(8)-2)/2,(-(Sqrt(8)+2))/2}"

There are two roots: \(\sqrt{2} - 1\) and \(-\sqrt{2} - 1\).

But let’s consider the two friendly roots first: \(0\) and \(1\). If you plug \(z=0\) and \(z=1\) in the three other polynomials of the Gröbner basis, you easily find three simultaneous roots: \((1, 0, 0)\), \((0, 1, 0)\) and \((0, 0, 1)\).

You should do this exercise, because now we will plug \(z=\sqrt{2}-1\) and \(z=-\sqrt{2}-1\), and this will be less easy. In fact I will only treat \(z=\sqrt{2}-1\), and leave the other root as another exercise for you.

Let’s go. We plug \(z=\sqrt{2}-1\) in
the third polynomial. Clearly, we will obtain an univariate polynomial
in \(y\), and we will solve it. Or
rather, **Ryacas** will solve it:

yac_str(sprintf("p3 := %s", gbstrings[[3L]])) ## [1] "y*z^2-z^2/2+z^4/2" yac_str("PSolve(Eliminate(z, Sqrt(2)-1, p3), y)") ## [1] "-((Sqrt(2)-1)^4/2-(Sqrt(2)-1)^2/2)/(Sqrt(2)-1)^2"

Yacas is not always nice. It does not simplify this expression. So let’s look at its numerical approximation:

yac_str("N(PSolve(Eliminate(z, Sqrt(2)-1, p3), y))") ## [1] "0.4142135623"

We recognize \(\sqrt{2}-1\). Now we plug \(z=\sqrt{2}-1\) into the second polynomial and we solve the resulting polynomial of degree \(2\) in \(y\):

yac_str(sprintf("p2 := %s", gbstrings[[2L]])) ## [1] "((-9)*z^3)/2-2*z^5-z^6+6*z^4-y+z+z^7/2+y^2" yac_str("N(PSolve(Eliminate(z, Sqrt(2)-1, p2), y))") ## [1] "{0.5857864376,0.4142135623}"

There are two roots: \(y=\sqrt{2}-1\) and the other one has no importance, since we didn’t find it before.

So finally we just have to plug \(y=\sqrt{2}-1\) and \(z=\sqrt{2}-1\) into the first polynomial and we solve the resulting univariate polynomial in \(x\):

yac_str(sprintf("p1 := %s", gbstrings[[1L]])) ## [1] "y+x+z^2-1" yac_str( "N(PSolve(Eliminate(y, Sqrt(2)-1, Eliminate(z, Sqrt(2)-1, p1)), x))" ) ## [1] "0.4142135623"

Again? Well, we find \(x = \sqrt{2}-1\). Finally we found one simultaneous root: \[(x, y, z) = (\sqrt{2}-1, \sqrt{2}-1, \sqrt{2}-1).\] Now I leave you with the promised exercise: treat the case \(z=-\sqrt{2}-1\). I give you the solution, you will find: \[(x, y, z) = (-\sqrt{2}-1, -\sqrt{2}-1, -\sqrt{2}-1).\]

Isn’t it impressive? I started to learn about Gröbner bases only a few
days ago, and I still have to learn. I hope some applications will be
implemented in **qspray** in the future, and/or exposed on
this blog.

The new version of **qspray** is not on CRAN yet. If you
are in a hurry, you can install the development version hosted
on Github.

**leave a comment**for the author, please follow the link and comment on their blog:

**Saturn Elephant**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.