# How to use Mathematica's GroebnerBasis to automate solving system of polynomials

Currently, Only way I know to use GroebnerBasis to solve system of polynomials require lots of manual steps. I'd like to know how to automate this process and if Mathematica has commands to help with this.

For an example, given these 4 polynomials in x,y,z,r

vars = {x,y,z,r}
eq1  = 3 x^2+2 y z-2 x r==0;
eq2  = 2 x z - 2 y r==0;
eq3  = 2 x y-2 z-2 z r==0;
eq4  = x^2+y^2+z^2-1==0;
opt={MonomialOrder->Lexicographic};


One can solve these using Solve

Solve[{eq1,eq2,eq3,eq4},vars]


And get these solutions

 {{x->-1,y->0,z->0,r->-(3/2)},
{x->-(2/3),y->-(1/3),z->-(2/3),r->-(4/3)},
etc....
}


To obtain the same solutions using GroebnerBasis, I start by finding basis

 g=GroebnerBasis[{eq1,eq2,eq3,eq4},vars,opt]


which gives

36 r-225 r^2-493 r^3-116 r^4+212 r^5+96 r^6
-4 z+25 r z+53 r^2 z+24 r^3 z
-189 r-4338 r^2-2076 r^3+1928 r^4+960 r^5+1105 z^2
5 r y+5 r^2 y+2 z+r z-r^2 z
4347 r+6291 r^2+12 r^3-2796 r^4-864 r^5+2210 y z+2210 r y z
-9945+10782 r+55171 r^2+25232 r^3-22556 r^4-13344 r^5+9945 y^2
13747 r+23859 r^2+4788 r^3-10604 r^4-5280 r^5+3315 x+3315 y z


There are 7 basis polynomials. Here comes the manual steps.

I look to see which polynomial above has only one variable in it. In this case, one is lucky. There is one. The first one depends only on r. So use that to solve for r

 Solve[g[[1]]==0,r]


Now I look to see which other polynomial in the basis has r and only one other variable to solve for.

This will be the second one in the list. So use it to solve for z using the first r value found above. (this process is very much like Gaussian elimination, the back substitution phase)

 Solve[(g[[2]]/.r->-3/2)==0,z]


Now we know r and z. So now will look for a basis polynomial which has r and z in it and only one more unknown to solve for. This will be the 4th one. Hence

 Solve[(g[[4]]/.{r->-3/2,z->0})==0,y]


Now we know r,z,y. Then look for polynomial which has these and one more unknown, which is x. This will be the last basis. Hence

 Solve[(g[[7]]/.{r->-3/2,z->0,y->0})==0,x]


So the one of the solution from above is

 {x->-1,y->0,z->0,r->-(3/2)}


This is one of the solutions found by Solve. To find the others, I repeat the above process, but now starting with the second r solution found in the first step of the process.

But this manual process is time consuming.

Does Mathematica have build in functions to automate this process once the GroebnerBasis are found?

• In what way will this be different than what you get from Solve? I.e., when you ask "Does Mathematica have build in functions to automate this process (...)", isn't the answer "Yes, the function is called Solve."? – Marius Ladegård Meyer Jan 26 '18 at 14:13
• @MariusLadegårdMeyer That's a bit different, because Solve is a general equation solver, it doesn't only deal with polynomial, so it'll have to first "check" the type of the system before solving it, which is just a waste of time if we've already know the type of the system. The time wasting can be considerable large… – xzczd Jan 27 '18 at 11:37
• @MariusLadegårdMeyer …Though I'm not familiar with polynomial equation, I've played quite a bit with linear equations, consider the following example about Solve vs LinearSolve: n = 1000; m = RandomReal[1, {n, n}]; b = RandomReal[1, n]; test2 = LinearSolve[m, b]; // AbsoluteTiming var = x /@ Range@n; test = var /. First@Solve[m.var == b // Thread, var]; // AbsoluteTiming Notice if one just has the equations at hand, m and b can be extracted as follows: {b2, m2} = {-1, 1} CoefficientArrays[m.var == b // Thread, var]; – xzczd Jan 27 '18 at 11:38

As noted in comments, Solve does this. An efficient "by hand" approach is to solve for the first equation, plug the solutions into the next, solve, join corresponding solutions. This can be done in two steps. First the example set-up.

vars = {x, y, z, r};
eq1 = 3 x^2 + 2 y z - 2 x r == 0;
eq2 = 2 x z - 2 y r == 0;
eq3 = 2 x y - 2 z - 2 z r == 0;
eq4 = x^2 + y^2 + z^2 - 1 == 0;
opt = {MonomialOrder -> Lexicographic};
g = GroebnerBasis[{eq1, eq2, eq3, eq4}, vars, opt]

(* Out[108]= {36 r - 225 r^2 - 493 r^3 - 116 r^4 + 212 r^5 +
96 r^6, -4 z + 25 r z + 53 r^2 z + 24 r^3 z, -189 r - 4338 r^2 -
2076 r^3 + 1928 r^4 + 960 r^5 + 1105 z^2,
5 r y + 5 r^2 y + 2 z + r z - r^2 z,
4347 r + 6291 r^2 + 12 r^3 - 2796 r^4 - 864 r^5 + 2210 y z +
2210 r y z, -9945 + 10782 r + 55171 r^2 + 25232 r^3 - 22556 r^4 -
13344 r^5 + 9945 y^2,
13747 r + 23859 r^2 + 4788 r^3 - 10604 r^4 - 5280 r^5 + 3315 x +
3315 y z} *)


Here is a helper function. Given a partial solution as a list of rules, and a polynomial, plug values into the equation and solve it. I am skipping details about whether/when one should expect the new polynomial to be univariate: a sufficient condition is that the system have a finite solution set (this is a consequence of the Gianni/Kalkbrenner theorems from I believe 1987).

solvenext[lastsol_, poly_] :=
Map[Join[lastsol, #] &, Solve[(poly /. lastsol) == 0]]


Now we deploy it, iterating from first polynomial in the Groebner basis to last. Fold is handy for this (okay, I first did it with a Do loop, then decided not to act my age and went with a functional programming construct instead).

Fold[
Apply[Join, Table[solvenext[prevsol, #2], {prevsol, #1}]] &, {{}}, g]

(* Out[149]= {{r -> -(3/2), z -> 0, y -> 0, x -> -1}, {r -> -(4/3),
z -> -(2/3), y -> -(1/3), x -> -(2/3)}, {r -> -(4/3), z -> 2/3,
y -> 1/3, x -> -(2/3)}, {r -> -1, z -> -1, y -> 0,
x -> 0}, {r -> -1, z -> -1, y -> 0, x -> 0}, {r -> -1, z -> 1,
y -> 0, x -> 0}, {r -> -1, z -> 1, y -> 0, x -> 0}, {r -> 0, z -> 0,
y -> -1, x -> 0}, {r -> 0, z -> 0, y -> 1, x -> 0}, {r -> 1/8,
z -> -(Sqrt[(11/2)]/8), y -> (3 Sqrt[11/2])/8,
x -> -(3/8)}, {r -> 1/8, z -> Sqrt[11/2]/8,
y -> -((3 Sqrt[11/2])/8), x -> -(3/8)}, {r -> 3/2, z -> 0, y -> 0,
x -> 1}} *)

• Had to edit...juxtaposing a pair of digits, I had put Gianni and Kalkbrenner back in the 19th century. I hope they enjoyed their brief stay therein. – Daniel Lichtblau Jan 27 '18 at 17:02
• Is Roots a better choice compared to Solve in this case? – xzczd Jan 28 '18 at 6:12
• @xzczd Roots avoids a bit of overhead but the replacement rule form of the Solve result is more convenient for the substitutions. – Daniel Lichtblau Jan 28 '18 at 15:31