Ask Your Question
3

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

asked 2021-05-20 12:28:20 +0200

updated 2021-05-20 12:30:07 +0200

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]
edit retag flag offensive close merge delete

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 $\bar{\mathbb{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 ( 2021-05-20 16:00:32 +0200 )edit

A simpler variant of that idea (using linear algebra over $\mathbb{Q}$ instead of arithmetic over $\bar{\mathbb{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 $\mathbb{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 $\mathbb{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 ( 2021-05-20 20:27:03 +0200 )edit
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 ( 2021-05-20 20:53:56 +0200 )edit
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 ( 2021-05-21 12:37:33 +0200 )edit

@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 ( 2021-05-21 21:47:23 +0200 )edit

2 Answers

Sort by » oldest newest most voted
4

answered 2021-05-24 01:25:28 +0200

rburing gravatar image

updated 2021-05-24 10:58:49 +0200

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 $\mathbb{C}$, $I \subset R$ a zero-dimensional radical ideal and $f \in R$. To show $ \langle 1 \rangle = I + \langle f \rangle$ it suffices (by Hilbert's Nullstellensatz) to show $\varnothing = V(I + \langle f \rangle) = V(I) \cap 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 $\mathbb{C}$-vector space.

  • Consider the $\mathbb{C}$-linear "multiplication by $f$" map $M_f: g \mapsto f \cdot g$ on $R/I$. Let $g$ be an eigenvector of $M_f$ with eigenvalue $\lambda$, so $f \cdot g = \lambda \cdot g$ or $(f-\lambda) g = 0$ holds in $R/I$, i.e. $(f-\lambda)g \in I$. The fact that $g \neq 0$ in $R/I$ means $g \not\in I$, hence there is a point $p \in V(I)$ with $g(p) \neq 0$. Now from $(f-\lambda)g \in I$ it follows that $\lambda = f(p)$ is the value of $f$ at the point $p$.

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

So the values of $f$ on $V(I)$ are exactly the eigenvalues of $M_f$, 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 $\mathbb{C}$-vector space (called the normal basis). The matrix of $M_f$ 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 $\varnothing = V(I + \langle f \rangle)$ and $\langle 1 \rangle = I + \langle f \rangle$.

Since $M_f$ is invertible, we can find a polynomial representative of $f^{-1} + 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 = f^{-1} \cdot f$ in $R/I$, i.e. $1 = f^{-1} \cdot f + g$ for some $g \in I$, and we can express $g$ in terms of generators by applying the division algorithm to $1 - f^{-1} \cdot f$ (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
edit flag offensive delete link more

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(M_f)$.

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 ( 2021-05-24 04:59:30 +0200 )edit

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 ( 2021-05-24 10:47:25 +0200 )edit
1

@rburing This is a very nice answer.

mwageringel gravatar imagemwageringel ( 2021-05-24 11:06:41 +0200 )edit

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

Sébastien Palcoux gravatar imageSébastien Palcoux ( 2021-07-02 15:20:30 +0200 )edit

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

Sébastien Palcoux gravatar imageSébastien Palcoux ( 2021-07-09 11:30:05 +0200 )edit
3

answered 2021-05-22 23:52:48 +0200

mwageringel gravatar image

updated 2021-05-23 12:54:56 +0200

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
edit flag offensive delete link more

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 ( 2021-05-23 16:50:10 +0200 )edit

Very interesting and useful answer!

Sébastien Palcoux gravatar imageSébastien Palcoux ( 2021-05-24 04:25:06 +0200 )edit

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 ( 2021-05-24 10:54:42 +0200 )edit

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 ( 2021-05-24 11:19:02 +0200 )edit

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

Sébastien Palcoux gravatar imageSébastien Palcoux ( 2021-07-02 11:24:53 +0200 )edit

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: 2021-05-20 12:28:20 +0200

Seen: 420 times

Last updated: May 24