Gröbner Bases for Polynomial Ideals and Applications

Authors: Christian Eder, Rafael Mohr and Mohab Safey El Din

How to use Gröbner bases in OSCAR explained on the computation of Hilbert series

julia> R, (X0, X1, X2, X3) = graded_polynomial_ring(QQ, ["X$i" for i in 0:3]);

julia> I = ideal(R, [X0*X2 - X1^2, X1*X3 - X2^2, X0*X3 - X1*X2]);

julia> A, _ = quo(R, I);

julia> hilbert_series(A)
(2*t^3 - 3*t^2 + 1, (-t + 1)^4)

julia> p = next_prime(rand(1:10000))
4783

julia> R, (X0, X1, X2, X3) = graded_polynomial_ring(GF(p), ["X$i" for i in 0:3]);

julia> I = ideal(R, [X0*X2 - X1^2, X1*X3 - X2^2, X0*X3 - X1*X2]);

julia> gb = groebner_basis_f4(forget_grading(I));

julia> lm_ideal = ideal(R, (leading_monomial).(gens(gb)));

julia> A, _ = quo(R, lm_ideal);

julia> hilbert_series(A)
(2*t^3 - 3*t^2 + 1, (-t + 1)^4)

Computing Gröbner bases with respect to different monomial orderings

Applying elimination

julia> R, (x, y, z) = polynomial_ring(QQ, ["x","y","z"]);

julia> o = degrevlex([x, y])*degrevlex([z]);

julia> I = ideal(R, [x^2 + y + z - 1, x + y^2 + z - 1, x + y + z^2 - 1]);

julia> groebner_basis(I, ordering = o)
Gröbner basis with elements
1 -> z^6 - 4*z^4 + 4*z^3 - z^2
2 -> 2*y*z^2 + z^4 - z^2
3 -> x + y + z^2 - 1
4 -> y^2 + x + z - 1
with respect to the ordering
degrevlex([x, y])*degrevlex([z])

julia> eliminate(I, [x,y]) # internally computes the above Gröbner basis
Ideal generated by
  z^6 - 4*z^4 + 4*z^3 - z^2

julia> groebner_basis(I)
Gröbner basis with elements
1 -> z^2 + x + y - 1
2 -> y^2 + x + z - 1
3 -> x^2 + y + z - 1
with respect to the ordering
degrevlex([x, y, z])

Using the FGLM algorithm for a change of ordering

julia> R, (x1, x2, x3, x4, x5, x6, x7, x8) = polynomial_ring(GF(65521), ["x1","x2","x3","x4","x5","x6", "x7","x8"]);

julia> I = ideal(R, [x1 + 2*x2 + 2*x3 + 2*x4 + 2*x5 + 2*x6 + 2*x7 + 2*x8 - 1, x1^2 - x1 + 2*x2^2 + 2*x3^2 + 2*x4^2 + 2*x5^2 + 2*x6^2 + 2*x7^2 + 2*x8^2, 2*x1*x2 + 2*x2*x3 - x2 + 2*x3*x4 + 2*x4*x5 + 2*x5*x6 + 2*x6*x7 + 2*x7*x8, 2*x1*x3 + x2^2 + 2*x2*x4 + 2*x3*x5 - x3 + 2*x4*x6 + 2*x5*x7 + 2*x6*x8, 2*x1*x4 + 2*x2*x3 + 2*x2*x5 + 2*x3*x6 + 2*x4*x7 - x4 + 2*x5*x8, 2*x1*x5 + 2*x2*x4 + 2*x2*x6 + x3^2 + 2*x3*x7 + 2*x4*x8 - x5, 2*x1*x6 + 2*x2*x5 + 2*x2*x7 + 2*x3*x4 + 2*x3*x8 - x6, 2*x1*x7 + 2*x2*x6 + 2*x2*x8 + 2*x3*x5 + x4^2 - x7]);

julia> @time groebner_basis_f4(I);
  0.010955 seconds (57.13 k allocations: 3.081 MiB)

julia> dim(I) == 0
true

julia> @time fglm(I, destination_ordering = lex(R));
  0.062046 seconds (130.52 k allocations: 3.313 MiB)

Compute real solutions

julia> R, (t1, t2) = polynomial_ring(QQ, ["t1", "t2"]);

julia> psi = [rand(R, 6:6, 3:5, -10:10) for _ in 1:3]
3-element Vector{QQMPolyRingElem}:
 1//2*t1^5*t2^3 - 2*t1^3*t2^4 + 2//9*t1^3*t2^3
 -3//5*t1^5*t2^5 - 8//7*t1^3*t2^4 - 17//4*t1^3*t2^3
 -23//12*t1^5*t2^5 - 3//7*t1^4*t2^3 - 13//6*t1^3*t2^5 - 9*t1^3*t2^4

julia> D = sum([(f - rand(Int))^2 for f in psi]);

julia> J = jacobian_ideal(D);

julia> S, (t1, t2, s) = polynomial_ring(QQ, ["t1", "t2", "s"]);

julia> phi = hom(R, S, [t1, t2]);

julia> f = sum([rand(Int) * p for p in minors(jacobian_matrix(psi), 2)]);

julia> JJ = phi(J) + ideal(S, s*phi(f) - 1);

julia> sols, param = real_solutions(JJ);

julia> length(sols)
8

julia> degree(param.elim)
62

Computing non-proper points

julia> R, (x, y, z, t) = polynomial_ring(QQ, ["x", "y", "z", "t"]);

julia> graph_id = ideal(R, [x*y, x*z, y*z, x*(x-1), x-t]);

julia> w = [1 1 1 0]
1×4 Matrix{Int64}:
 1  1  1  0

julia> H = homogenizer(R, w, "w"); clo = H(graph_id)
Ideal generated by
  -x + t*w
  t^2 - t
  z*t
  y*t
  x*t - x
  y*z
  x*z
  x*y
  x^2 - x*w

julia> S = base_ring(clo)
Multivariate polynomial ring in 5 variables over QQ graded by
  x -> [1]
  y -> [1]
  z -> [1]
  t -> [0]
  w -> [1]

julia> x, y, z, t, w = gens(S);

julia> irrel_ideal = ideal(S, [x,y,z,w]);

julia> eliminate(saturation(clo + ideal(S,x), irrel_ideal), [x,y,z,w])
Ideal generated by
  t

julia> eliminate(saturation(clo + ideal(S,y), irrel_ideal), [x,y,z,w])
Ideal generated by
  t^2 - t

julia> eliminate(saturation(clo + ideal(S,z), irrel_ideal), [x,y,z,w])
Ideal generated by
  t^2 - t

julia> eliminate(saturation(clo + ideal(S,w), irrel_ideal), [x,y,z,w])
Ideal generated by
  t

Computing conormal spaces and Whitney stratifications

julia> R, (x, y, z, w) = graded_polynomial_ring(QQ, ["x","y","z","w"]);

julia> f = y^2*w^2 + z^3*w - x^2*z^2;

julia> J = jacobian_matrix(f)
[          -2*x*z^2]
[           2*y*w^2]
[-2*x^2*z + 3*z^2*w]
[     2*y^2*w + z^3]

julia> I_sing = ideal(R, [f, minors(J, 1)...])
Ideal generated by
  -x^2*z^2 + y^2*w^2 + z^3*w
  -2*x*z^2
  2*y*w^2
  -2*x^2*z + 3*z^2*w
  2*y^2*w + z^3

julia> minimal_primes(I_sing)
2-element Vector{MPolyIdeal{MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}}}:
 Ideal (w, z)
 Ideal (z, y)

julia> (comp -> comp[2]).(primary_decomposition(I_sing))
6-element Vector{MPolyIdeal{MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}}}:
 Ideal (w, z)
 Ideal (z, y)
 Ideal (w, z, y)
 Ideal (w, z, x)
 Ideal (z, y, x)
 Ideal (x, y, z, w)

julia> R, (x0,x1,x2,x3,a0,a1,a2,a3) = polynomial_ring(QQ, ["x0","x1","x2","x3","a0","a1","a2","a3"]);

julia> f = x0^6 + x0^4*x1*x2 + x2^3*x3^3;

julia> J = jacobian_matrix(f)[1:4, 1:1]
[6*x0^5 + 4*x0^3*x1*x2]
[              x0^4*x2]
[x0^4*x1 + 3*x2^2*x3^3]
[          3*x2^3*x3^2]

julia> JJ = hcat(J, R[a0; a1; a2; a3])
[6*x0^5 + 4*x0^3*x1*x2   a0]
[              x0^4*x2   a1]
[x0^4*x1 + 3*x2^2*x3^3   a2]
[          3*x2^3*x3^2   a3]

julia> id = ideal(R, minors(JJ, 2));

julia> I_sing = ideal(R, minors(J, 1));

julia> conorm = saturation(ideal(R, f) + id, I_sing);

julia> Y = radical(ideal(R,f) + I_sing)
Ideal generated by
  x0
  x2*x3

julia> id_fib_Y = saturation(conorm + Y, ideal(R, [x0,x1,x2,x3]));

julia> prim_comps = (comp -> comp[2]).(primary_decomposition(id_fib_Y));

julia> [eliminate(comp, [a0,a1,a2,a3]) for comp in prim_comps]
6-element Vector{MPolyIdeal{QQMPolyRingElem}}:
 Ideal (x2, x0)
 Ideal (x3, x2, x0)
 Ideal (x3, x0)
 Ideal (x3, x2, x0)
 Ideal (x3, x1, x0)
 Ideal (x2, x0, 4*x1^3 + 27*x3^3)

An example from combinatorics

julia> F = GF(rand_bits_prime(ZZ, 20))
Prime field of characteristic 646637

julia> R, (z0, z1, z2, t, u, x) = polynomial_ring(F, ["z0","z1","z2","t", "u","x"]);

julia> P = (u^3-t*(u^7+1))*x-u^3+t*z0+t*u*z1+inv(F(2))*t*u^2*z2;

julia> sys = [P,derivative(P,u),derivative(P,x)];

julia> I = saturation(ideal(R, sys), ideal(R, u));

julia> I = eliminate(ideal(R, sys), [x]);

julia> sys[3]
646636*t*u^7 + 646636*t + u^3

julia> S, (z0, z1, z2, t, u, u1, u2, x) = polynomial_ring(F, ["z0","z1","z2","t","u","u1","u2","x"]);

julia> phi = hom(R,S,[z0,z1,z2,t,u,x]);

julia> I = phi(I);

julia> I += ideal(S, [f(z0,z1,z2,t,u1,u1,u2,x) for f in gens(I)]);

julia> I += ideal(S, [f(z0,z1,z2,t,u2,u1,u2,x) for f in gens(I)]);

julia> I = saturation(I,ideal(S,(u-u1)*(u1-u2)*(u-u2)));

julia> res = eliminate(I, [z1,z2,u,u1,u2,x]);

julia> factor(first(gens(res)))
1 * (z0^35*t^35 + 646636*z0^31*t^28 + z0^30*t^28 + 646636*z0^29*t^28 + 5*z0^28*t^28 + 646636*z0^25*t^21 + z0^24*t^21 + 3*z0^23*t^21 + 646633*z0^22*t^21 + 10*z0^21*t^21 + z0^19*t^14 + 646636*z0^18*t^14 + 5*z0^17*t^14 + 3*z0^16*t^14 + 646631*z0^15*t^14 + 10*z0^14*t^14 + z0^13*t^7 + 646636*z0^12*t^7 + 3*z0^10*t^7 + z0^9*t^7 + 646633*z0^8*t^7 + 5*z0^7*t^7 + 646636*z0 + 1)

Checking for local complete intersections

julia> R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"]);

julia> I = ideal(R, [x*y, x*z, y*z]);

julia> conorm = subquotient(matrix(gens(I)), matrix(gens(I^2)))
Subquotient of Submodule with 3 generators
1 -> x*y*e[1]
2 -> x*z*e[1]
3 -> y*z*e[1]
by Submodule with 6 generators
1 -> x^2*y^2*e[1]
2 -> x^2*y*z*e[1]
3 -> x*y^2*z*e[1]
4 -> x^2*z^2*e[1]
5 -> x*y*z^2*e[1]
6 -> y^2*z^2*e[1]

julia> fitting_ideal(conorm, 1)
Ideal generated by
  y*z
  x*z
  x*y

julia> fitting_ideal(conorm, 2)
Ideal generated by
  z
  y
  x

julia> equidimensional_decomposition_weak(I)
1-element Vector{MPolyIdeal{QQMPolyRingElem}}:
 Ideal (x*y, x*z, y*z)

julia> F = fitting_ideal(ideal_as_module(I), 2);

julia> codim(F) >= codim(I) + 1
true

julia> I1 = ideal(R, [x-1, y-1]);

julia> I2 = ideal(R, [x, y, z]);

julia> I = intersect(I1, I2);

julia> conorm = subquotient(matrix(gens(I)), matrix(gens(I^2)));

julia> Pk(conorm, k) = saturation(I, fitting_ideal(conorm, k)) + fitting_ideal(conorm, k-1);

julia> Pk(conorm, 2) == ideal(R, [x-1,y-1])
true

julia> Pk(conorm, 3) == ideal(R, [z,y,x])
true