Processing math: 76%

First time here? Check out the FAQ!

Ask Your Question
3

Quick (trivial) Groebner basis but too long lift of one

asked 3 years ago

updated 3 years ago

I have two (equivalent) sub-systems of polynomial equations of 10 variables and 12 equations, linked with one extra equations (so a total of 20 variable and 25 equations).

I can very quicky compute the Groebner basis of each sub-system. Then I can put them together, with the extra equation, and then I can quicky get the Groebner basis of the full system, which turns out to be trivial (total computation time: less than 3 minutes). See all the details in Appendix.

A direct computation of the Groebner basis of the full system is too long (I stopped after 1 day). The problem is that I need to compute the lift of 1 of the original full system (i.e. express 1 as a linear combination (with polynomial coefficients) of the original generators, by using the lift command list(R.one().lift(Id))), but it is also too long...

Now I expect that the application of a strategy involving the two sub-systems (as above) should exist to quickly get the lift of 1, but how?

I tried to compute this lift in two times, using the Groebner basis of the sub-systems, but it is also too long...


Appendix

First sub-system (10 variables and 12 equations)

sage: R.<u0,u1,u2,v0,v1,v2,v3,v4,v5,v6>=PolynomialRing(QQ,10)
sage: A1=[u0 + 7/5*u1 + 7/5*u2 - 4/125,
....:   5*v0 + 5*v1 + 7*v3 + 7*v5 + 1/5,
....:   25*v0^2 + 25*v1^2 + 35*v3^2 + 35*v5^2 - 4/5,
....:   5*v0^3 + 5*v1^3 + 7*v3^3 + 7*v5^3 - v0^2 + 1/125,
....:   5*v0*v1^2 + 5*v1*v2^2 + 7*v3*v4^2 + 7*v5*v6^2 + 1/125,
....:   5*u0*v1 - v1^2 + 7*u1*v3 + 7*u2*v5 + 1/125,
....:   5*v1 + 5*v2 + 7*v4 + 7*v6 + 1/5,
....:   25*v0*v1 + 25*v1*v2 + 35*v3*v4 + 35*v5*v6 + 1/5,
....:   5*v0^2*v1 + 5*v1^2*v2 + 7*v3^2*v4 + 7*v5^2*v6 - v1^2 + 1/125,
....:   25*v1^2 + 25*v2^2 + 35*v4^2 + 35*v6^2 - 4/5,
....:   5*v1^3 + 5*v2^3 + 7*v4^3 + 7*v6^3 - u0 + 1/125,
....:   5*u0*v2 - v2^2 + 7*u1*v4 + 7*u2*v6 + 1/125]
sage: Id=R.ideal(A1)
sage: %time G1=Id.groebner_basis()
CPU times: user 1.33 s, sys: 0 ns, total: 1.33 s
Wall time: 1.33 s
sage: C1=[g for g in G1]

Second (equivalent) sub-system:

sage: R.<y0,y1,y2,z0,z1,z2,z3,z4,z5,z6>=PolynomialRing(QQ,10)
sage: A2=[y0 + 7/5*y1 + 7/5*y2 - 4/125,
....:   5*z0 + 5*z1 + 7*z3 + 7*z5 + 1/5,
....:   25*z0^2 + 25*z1^2 + 35*z3^2 + 35*z5^2 - 4/5,
....:   5*z0^3 + 5*z1^3 + 7*z3^3 + 7*z5^3 - y0 + 1/125,
....:   5*y0*z0 - z0^2 + 7*y1*z3 + 7*y2*z5 + 1/125,
....:   5*z0*z1^2 + 5*z1*z2^2 + 7*z3*z4^2 + 7*z5*z6^2 - z1^2 + 1/125,
....:   5*z1 + 5*z2 + 7*z4 + 7*z6 + 1/5,
....:   25*z0*z1 + 25*z1*z2 + 35*z3*z4 + 35*z5*z6 + 1/5,
....:   5*z0^2*z1 + 5*z1^2*z2 + 7*z3^2*z4 + 7*z5^2*z6 + 1/125,
....:   5*y0*z1 - z1^2 + 7*y1*z4 + 7*y2*z6 + 1/125,
....:   25*z1^2 + 25*z2^2 + 35*z4^2 + 35*z6^2 - 4/5,
....:   5*z1^3 + 5*z2^3 + 7*z4^3 + 7*z6^3 - z2^2 + 1/125]
sage: Id=R.ideal(A2)
sage: %time G2=Id.groebner_basis()
CPU times: user 2.08 s, sys: 0 ns, total: 2.08 s
Wall time: 2.08 s
sage: C2=[g for g in G2]

The full system:

sage: R.<u0,u1,u2,v0,v1,v2,v3,v4,v5,v6,y0,y1,y2,z0,z1,z2,z3,z4,z5,z6>=PolynomialRing(QQ,20)
sage: C=C1+C2+[5*u0*z2 + 5*u1*z4 + 7*u2*z6 - u0 + 1/125]  # the last one is the extra one
sage: Id=R.ideal(C)
sage: %time G=Id.groebner_basis()
CPU times: user 2min 29s, sys: 16 ms, total: 2min 29s
Wall time: 2min 29s
sage: G
[1]
Preview: (hide)

Comments

It's not what you asked, and probably suboptimal, but: since the two individual ideals are zero-dimensional, and the respective varieties each contain 14 points over ˉQ, you could also loop through all 196 pairs of points and substitute them into the extra polynomial, and check that you never get zero.

rburing gravatar imagerburing ( 3 years ago )

A simpler variant of that idea (using linear algebra over Q instead of arithmetic over ˉQ) is to use that C1 + C2 is a Groebner basis of the ideal I generated by C1 + C2; hence there is a normal basis (basis of the quotient R/I as a Q-vector space) consisting of 196 monomials, and multiplication by f = 5*u0*z2 + 5*u1*z4 + 7*u2*z6 - u0 + 1/125 can be written as a 196x196 matrix over Q, and the values of f on V(I) are the eigenvalues of this matrix, so it suffices to show that the matrix does not have 0 as an eigenvalue, or that it has full rank, or that its determinant is nonzero (which is easy).

rburing gravatar imagerburing ( 3 years ago )
1

@rburing What you wrote in your first comment was done firstly, but that requires to use numerical solutions... The computation with Groebner bases only is cleaner (in my opinion) and still very quick in this case. Anyway, my problem is that I need the lift of 1, because the Groebner basis must be also trivial in char p>0, except for finitely many p, contained in the set of prime divisors of the denominators in the lift of 1, but I need to classify all p for which the Groebner basis is non-trivial: for p<1000000, there are exactly six one: 11, 29, 41, 467, 1423, 242057. With a lift of 1, I can easily classify all such p.

Sébastien Palcoux gravatar imageSébastien Palcoux ( 3 years ago )
1

Interesting problem. The numerator of the determinant of the aforementioned matrix has 785 digits and its factorization is 11 × 467 × 1423 × 242057 × 1477913 × 2673859 × 9528289 × 3869785541 × 715114922077295989 × 75154031037086343317 × 652748578372455171701 × (to be determined).

rburing gravatar imagerburing ( 3 years ago )

@rburing Very interesting!! If I am not mistaken, the primes that I am looking for are exactly the prime divisors of the determinant of this 196x196 matrix, right? If so, 29 and 41 should also be divisors. Did you miss them, or am I mistaken? Can you write the numerator of the determinant (with its 785 digits)?

Sébastien Palcoux gravatar imageSébastien Palcoux ( 3 years ago )

2 Answers

Sort by » oldest newest most voted
4

answered 3 years ago

rburing gravatar image

updated 3 years ago

In this rather special case we can find a lift of 1 quickly and efficiently using linear algebra.

Let R be a polynomial ring over C, IR a zero-dimensional radical ideal and fR. To show 1=I+f it suffices (by Hilbert's Nullstellensatz) to show =V(I+f)=V(I)V(f), i.e. that f does not vanish on V(I). Since I is radical, the ideal of polynomials in R vanishing on V(I) is exactly I (by the Nullstellensatz). By definition of zero-dimensional, the quotient R/I is a finite-dimensional C-vector space.

  • Consider the C-linear "multiplication by f" map Mf:gfg on R/I. Let g be an eigenvector of Mf with eigenvalue λ, so fg=λg or (fλ)g=0 holds in R/I, i.e. (fλ)gI. The fact that g0 in R/I means gI, hence there is a point pV(I) with g(p)0. Now from (fλ)gI it follows that λ=f(p) is the value of f at the point p.

  • Conversely, if P is the minimum polynomial of Mf then P(Mf)=0 implies P(f)=0 in R/I (by evaluating at 1R/I), hence P(f)I and it follows that 0=P(f)(p)=P(f(p)) for all pV(I), i.e. all values f(p) for pV(I) are zeros of P, and hence they are eigenvalues of Mf.

So the values of f on V(I) are exactly the eigenvalues of Mf, and it suffices to show that 0 is not one of them.

Given a Groebner basis G of I, the multivariate polynomial division algorithm implies that the set of all monomials not divisible by any leading terms of G is a basis of R/I as a C-vector space (called the normal basis). The matrix of Mf with respect to this basis is easily computed:

def multiplication_matrix(f, I):
    G = I.groebner_basis()
    f = f.reduce(G)
    NB = I.normal_basis()
    M = Matrix(f.base_ring(), len(NB))
    for (i, m1) in enumerate(NB):
        g = (f*m1).reduce(G)
        M.set_column(i, [g.monomial_coefficient(m2) for m2 in NB])
    return M

We apply it to the problem at hand:

sage: f = 5*u0*z2 + 5*u1*z4 + 7*u2*z6 - u0 + 1/125
sage: I = R.ideal(C1 +  C2)
sage: M_f = multiplication_matrix(f, I)
sage: M_f.det() != 0
True

Hence f does not attain the value 0 on V(I), so =V(I+f) and 1=I+f.

Since Mf is invertible, we can find a polynomial representative of f1+I:

sage: NB = list(I.normal_basis())
sage: len(NB), NB.index(1)
(196, 195)
sage: finv_coeffs = M_f \ vector([0]*195 + [1])
sage: finv = sum(c*m for (c,m) in zip(finv_coeffs, NB))
sage: (f*finv).reduce(I)
1

Then we have 1=f1f in R/I, i.e. 1=f1f+g for some gI, and we can express g in terms of generators by applying the division algorithm to 1f1f (or finding a lift of it to I):

sage: lift_coeffs = [finv] + list((1 - finv*f).lift(I))
sage: sum(c*m for c,m in zip(lift_coeffs, [f] + I.gens()))
1
Preview: (hide)
link

Comments

Very good!

Then about the existence of solutions in char. p>0, in this case:

  • first we need to consider independently the finite set S of prime divisors of the (reduced) denominators in the Groebner basis of I,

  • next, out of S, there is a solution in char. p if p divides the (reduced) numerator n of det.

This "if" is not "iff" because both 41 \not \in S and does not divides n. In fact, the ideal I is not zero-dimensional in char. 41 but positive-dimensional (in fact 2-dim). I guess that 41 is the only prime in this alternative case.

Sébastien Palcoux gravatar imageSébastien Palcoux ( 3 years ago )

What your answer give is the lift of 1 with respect to C_1 + C_2 + extra. What is asked is the lift of 1 with respect to A_1 + A_2 + extra (with a bit of luck, it could be easier...).

Sébastien Palcoux gravatar imageSébastien Palcoux ( 3 years ago )
1

@rburing This is a very nice answer.

mwageringel gravatar imagemwageringel ( 3 years ago )

@rburing Can your method by adapted to lift C_1 with respect to A_1?

Sébastien Palcoux gravatar imageSébastien Palcoux ( 3 years ago )

@rburing see this new post motivating the problem in my previous comment.

Sébastien Palcoux gravatar imageSébastien Palcoux ( 3 years ago )
3

answered 3 years ago

mwageringel gravatar image

updated 3 years ago

Note that computing the lift of 1 with respect to your C modulo a prime is fast – it takes only a few seconds. Doing this for many primes can allow to lift the result to ℚ. In this case, the result F (see attachment uploaded here (31MB), for now, as the file seems too large for an attachment) is huge, though, with coefficients as large as 4826 and 4818 digits for numerator and denominator.

sage: F, C = load("lift_one_result_F_C.sobj")  # saved with Sage 9.3
sage: sum(f*c for f,c in zip(F, C))  # 1 minute
1

I am not sure if this solution would be unique. Maybe there is a smaller one.

To obtain the result with respect to A1 and A2, you could try to lift C1 and C2 individually, but this might make the result even larger.


Update

When computing the lift of 1 mod p, most of the polynomials in the result have the same monomials for their support, if p varies. This makes it very likely that the correct representative over ℚ has those same monomials as support and that the polynomials we computed mod p are the mod-p versions of the desired polynomials over ℚ.

This can fail for some primes which we need to ignore. In the function likely_support, we use a simple ad-hoc heuristic to find the exponents that probably occur in the ℚ-result: if a monomial occurs in 50% of the mod-p polynomials (i.e. modulo 50% of the primes or more), we expect that it should occur in the ℚ-result as well. The 50% threshold does not matter much. We could have picked 10% or 90% instead, I suppose.

We only keep the primes p for which all the polynomials in the lift of 1 mod p have the expected supports. Out of the 2800 primes we started with, 2358 remain.

Next, we use the Chinese remainder theorem CRT_list to lift the coefficients to ℤ modulo the product of the primes and then rational_reconstruction to obtain rational coefficients, both implemented in lift_poly_crt. For rational_reconstruction to work, the product of the primes must have about twice as many digits as the numerators and denominators in the result over ℚ (which I do not know how to determine in advance). In our case, we obtain a modulus with 10538 digits.

I have used cached_function.precompute in order to parallelize some computations because it is a bit less awkward to use than the @parallel decorator, but also because it allows to increase the number of primes until rational reconstruction works, without redoing all the previous slow computations.

The slowest part of the computation is lifting 1 modulo the 2800 primes, so this benefits a lot from parallelization. Each prime takes about 10 seconds, so this would take about 7.5 hours on a single thread.

R1.<u0,u1,u2,v0,v1,v2,v3,v4,v5,v6>=PolynomialRing(QQ,10)
A1=[u0 + 7/5*u1 + 7/5*u2 - 4/125,
  5*v0 + 5*v1 + 7*v3 + 7*v5 + 1/5,
  25*v0^2 + 25*v1^2 + 35*v3^2 + 35*v5^2 - 4/5,
  5*v0^3 + 5*v1^3 + 7*v3^3 + 7*v5^3 - v0^2 + 1/125,
  5*v0*v1^2 + 5*v1*v2^2 + 7*v3*v4^2 + 7*v5*v6^2 + 1/125,
  5*u0*v1 - v1^2 + 7*u1*v3 + 7*u2*v5 + 1/125,
  5*v1 + 5*v2 + 7*v4 + 7*v6 + 1/5,
  25*v0*v1 + 25*v1*v2 + 35*v3*v4 + 35*v5*v6 + 1/5,
  5*v0^2*v1 + 5*v1^2*v2 + 7*v3^2*v4 + 7*v5^2*v6 - v1^2 + 1/125,
  25*v1^2 + 25*v2^2 + 35*v4^2 + 35*v6^2 - 4/5,
  5*v1^3 + 5*v2^3 + 7*v4^3 + 7*v6^3 - u0 + 1/125,
  5*u0*v2 - v2^2 + 7*u1*v4 + 7*u2*v6 + 1/125]
I1=R1.ideal(A1)
C1=list(I1.groebner_basis())

R2.<y0,y1,y2,z0,z1,z2,z3,z4,z5,z6>=PolynomialRing(QQ,10)
A2=[y0 + 7/5*y1 + 7/5*y2 - 4/125,
  5*z0 + 5*z1 + 7*z3 + 7*z5 + 1/5,
  25*z0^2 + 25*z1^2 + 35*z3^2 + 35*z5^2 - 4/5,
  5*z0^3 + 5*z1^3 + 7*z3^3 + 7*z5^3 - y0 + 1/125,
  5*y0*z0 - z0^2 + 7*y1*z3 + 7*y2*z5 + 1/125,
  5*z0*z1^2 + 5*z1*z2^2 + 7*z3*z4^2 + 7*z5*z6^2 - z1^2 + 1/125,
  5*z1 + 5*z2 + 7*z4 + 7*z6 + 1/5,
  25*z0*z1 + 25*z1*z2 + 35*z3*z4 + 35*z5*z6 + 1/5,
  5*z0^2*z1 + 5*z1^2*z2 + 7*z3^2*z4 + 7*z5^2*z6 + 1/125,
  5*y0*z1 - z1^2 + 7*y1*z4 + 7*y2*z6 + 1/125,
  25*z1^2 + 25*z2^2 + 35*z4^2 + 35*z6^2 - 4/5,
  5*z1^3 + 5*z2^3 + 7*z4^3 + 7*z6^3 - z2^2 + 1/125]
I2=R2.ideal(A2)
C2=list(I2.groebner_basis())

R.<u0,u1,u2,v0,v1,v2,v3,v4,v5,v6,y0,y1,y2,z0,z1,z2,z3,z4,z5,z6>=PolynomialRing(QQ,20)
C=C1+C2+[5*u0*z2 + 5*u1*z4 + 7*u2*z6 - u0 + 1/125]  # the last one is the extra one
I=R.ideal(C)

########################################################################

@cached_function
def lift_one_modp(Ip):
    Fp = Ip.ring()(1).lift(Ip.ideal())
    print("completed p = %s" % Fp.ring().characteristic())
    return Fp

@cached_function
def lift_poly_crt(fps):
    ps = [fp.parent().characteristic() for fp in fps]
    modulus = prod(ps)
    Ds = [fp.dict() for fp in fps]
    supp = set([e for Di in Ds for e in Di])
    return R({e: CRT_list([Di[e].lift() if e in Di else 0 for Di in Ds],
                          list(ps)).rational_reconstruction(modulus) for e in supp})

def likely_support(fps):
    supps = [set(fp.dict().keys()) for fp in fps]
    supp_union = set.union(*supps)
    supp = set()
    half = len(supps) * 1/2  # exact amount does not matter
    for e in supp_union:
        k = 0
        for s in supps:
            if e in s:
                k += 1
            if k > half:  # if exponent vector is in more than half the supports, we will probably want to keep it
                supp.add(e)
                break
    return supp


ps = sorted(RecursivelyEnumeratedSet([next_prime(16000)], lambda p: [next_prime(p)])[:2800])  # many primes
Rps = [R.change_ring(GF(p)) for p in ps]
Ips = [I.change_ring(Rp).gens() for Rp in Rps]

CPUs = 16  # choose based on hardware
lift_one_modp.precompute(Ips, num_processes=CPUs)
assert all(type(r) != str for r in lift_one_modp.cache.values())  # otherwise, there was an exception
Fps = [lift_one_modp(Ip) for Ip in Ips]

supps = [likely_support(fps) for fps in zip(*Fps)]

# keep indices of primes for which all polynomials have the expected support
idxs = [j for j, Fp in enumerate(Fps) if all(set(fp.dict()) == supp for fp, supp in zip(Fp, supps))]

# restrict the results to those primes
ps1 = tuple(ps[j] for j in idxs)
Fps1 = [Fps[j] for j in idxs]
fpss = [tuple(fps) for fps in zip(*Fps1)]
modulus1 = prod(ps1)
floor(log(modulus1, 10))  # 10538 (must be large enough for rational reconstruction: this needs about twice as many digits as the numerators or denominators in the result have)

lift_poly_crt.precompute([(fps,) for fps in fpss], num_processes=CPUs)
assert all(type(r) != str for r in lift_poly_crt.cache.values())  # otherwise, there was an exception (if so, remove it from cache by lift_poly_crt.cache.clear())
F = [lift_poly_crt(fps) for fps in fpss]  # the final result
Preview: (hide)
link

Comments

2

Nice work! The common denominator is divisible by the 785 digit number I found in the comments. The coefficients of the lift of 1 are not uniquely determined because you can always add syzygies (lifts of 0), and conversely if you have two different lifts of 1 then their difference is a syzygy.

rburing gravatar imagerburing ( 3 years ago )

Very interesting and useful answer!

Sébastien Palcoux gravatar imageSébastien Palcoux ( 3 years ago )

With a bit of luck, the lift of 1 with respect to A_1+A_2+ extra could be easier...

Sébastien Palcoux gravatar imageSébastien Palcoux ( 3 years ago )

The trouble is that lifting 1 with respect to A_1 and A_2 is too slow in general, even modulo p it does not finish within a couple of minutes.

I have also tried to lift the elements of C_1,C_2 with respect to A_1,A_2 modulo p – this finishes quickly, so doing this for many p should allow lifting to ℚ again. However, if you compose this result with the F from my answer, you get very large polynomials with about 500k terms. (Using lift_coeffs from @rburing's answer instead of F would be a bit better, of course.)

mwageringel gravatar imagemwageringel ( 3 years ago )

@mwageringel How can I use your code to lift C_1 with respect to A_1?

Sébastien Palcoux gravatar imageSébastien Palcoux ( 3 years ago )

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 3 years ago

Seen: 809 times

Last updated: May 24 '21