Opened 4 years ago

Closed 19 months ago

#20145 closed defect (fixed)

Hilbert series bug

Reported by: stumpc5 Owned by:
Priority: major Milestone: sage-8.4
Component: commutative algebra Keywords: Hilbert series, polynomial ring
Cc: Merged in:
Authors: Simon King Reviewers: Travis Scrimshaw
Report Upstream: Reported upstream. Developers acknowledge bug. Work issues:
Branch: a1545cc (Commits) Commit: a1545cc6c843fcdb85cbca3c85ca8ce49de28616
Dependencies: #26243 Stopgaps:

Description (last modified by chapoton)

UPSTREAM: https://www.singular.uni-kl.de:8005/trac/ticket/750

This works correctly:

sage: n=4;m=10;P = PolynomialRing(QQ,n*m,"x"); x = P.gens(); M = Matrix(n,x)
sage: M
[ x0  x1  x2  x3  x4  x5  x6  x7  x8  x9]
[x10 x11 x12 x13 x14 x15 x16 x17 x18 x19]
[x20 x21 x22 x23 x24 x25 x26 x27 x28 x29]
[x30 x31 x32 x33 x34 x35 x36 x37 x38 x39]
sage: I = P.ideal(M.minors(2))
sage: I.hilbert_series().numerator()
-84*t^3 - 108*t^2 - 27*t - 1
sage: factor(I.hilbert_series().denominator())
(t - 1)^13

but if I increase the number of columns by one, I get

sage: n=4;m=11;P = PolynomialRing(QQ,n*m,"x"); x = P.gens(); M = Matrix(n,x)
sage: M
[ x0  x1  x2  x3  x4  x5  x6  x7  x8  x9 x10]
[x11 x12 x13 x14 x15 x16 x17 x18 x19 x20 x21]
[x22 x23 x24 x25 x26 x27 x28 x29 x30 x31 x32]
[x33 x34 x35 x36 x37 x38 x39 x40 x41 x42 x43]
sage: I = P.ideal(M.minors(2))
sage: I.hilbert_series().numerator()
120*t^33 - 3465*t^32 + 48180*t^31 - 429374*t^30 + 2753520*t^29 - 13522410*t^28 + 52832780*t^27 - 168384150*t^26 + 445188744*t^25 - 987193350*t^24 + 1847488500*t^23 + 1372406746*t^22 - 403422496*t^21 - 8403314*t^20 - 471656596*t^19 + 1806623746*t^18 + 752776200*t^17 + 752776200*t^16 - 1580830020*t^15 + 1673936550*t^14 - 1294246800*t^13 + 786893250*t^12 - 382391100*t^11 + 146679390*t^10 - 42299400*t^9 + 7837830*t^8 - 172260*t^7 - 468930*t^6 + 183744*t^5 - 39270*t^4 + 5060*t^3 - 330*t^2 + 1
sage: factor(I.hilbert_series().denominator())
(t - 1)^44

which is clearly not correct, as it can be checked with

sage: from sage.libs.singular.function import singular_function
sage: hilb = singular_function("hilb")
sage: std = singular_function("std")
sage: hilb(std(I))
// ** _ is no standard basis
...
//         1 t^0
//        30 t^1
//       135 t^2
//       120 t^3
// dimension (proj.)  = 13
// degree (proj.)   = 286

Change History (66)

comment:1 Changed 4 years ago by dimpase

Where is the error?

44 is the number of generators of P in n=4;m=11;P = PolynomialRing(QQ,n*m,"x");, and Sage simply returns (1-t)^P.ngens() as the denominator. Is this wrong mathematically?

comment:2 Changed 4 years ago by stumpc5

The correct answer, as given in the second computation hilb(std(I)), is

( 120*t^3 + 135*t^2 + 30*t + 1 ) / (t-1)^14

which is not what Sage outputs. (But the total degree -11 is still correct...)

comment:3 follow-ups: Changed 4 years ago by dimpase

well, for m=10 you only get the output above due to the fact that numerator has a factor (t-1)27 that gets cancelled.

And for m=11 the numerator is irreducible, so no cancellation happens. Does this mean that it's only due to this cancellation accident in the m=10 case you get the same output in Sage as you expect?

Numerator is computed by Singular and so it should not give any bug, I guess...

comment:4 in reply to: ↑ 3 ; follow-up: Changed 4 years ago by stumpc5

Replying to dimpase:

well, for m=10 you only get the output above due to the fact that numerator has a factor (t-1)27 that gets cancelled.

I am not sure I follow what you mean -- in the m=10 case, the (afaik) correct Hilbert series is

( 84*t^3 + 108*t^2 + 27*t + 1 ) / (1 - t)^13

which is what I get above. But in the m=11 case, the (afaik) correct Hilbert series is

( 120*t^3 + 135*t^2 + 30*t + 1 ) / (1 - t)^14

which is not what I get above. In particular, the numerator output of singular should be divisible by (1-t)^30, but - as you write - it is irreducible instead.

comment:5 in reply to: ↑ 3 Changed 4 years ago by stumpc5

Replying to dimpase:

Numerator is computed by Singular and so it should not give any bug, I guess...

At least luisfe on sage-devel (https://groups.google.com/d/msg/sage-devel/tyTDd5dEDGw/8FxH9hUzBwAJ) seems to confirm this to be a bug in Singular...

comment:6 in reply to: ↑ 4 ; follow-up: Changed 4 years ago by dimpase

Replying to stumpc5:

Replying to dimpase:

well, for m=10 you only get the output above due to the fact that numerator has a factor (t-1)27 that gets cancelled.

I am not sure I follow what you mean

You can read the code of I.hilbert_series to see that it takes I.hilbert_numerator(), which is computed by Singular (again, read the code), and divides it by (t-1)P.ngens().

-- in the m=10 case, the (afaik) correct Hilbert series is

( 84*t^3 + 108*t^2 + 27*t + 1 ) / (1 - t)^13

which is what I get above. But in the m=11 case, the (afaik) correct Hilbert series is

( 120*t^3 + 135*t^2 + 30*t + 1 ) / (1 - t)^14

which is not what I get above. In particular, the numerator output of singular should be divisible by (1-t)^30, but - as you write - it is irreducible instead.

So it is a bug in Singular, indeed.

comment:7 in reply to: ↑ 6 ; follow-up: Changed 4 years ago by stumpc5

Replying to dimpase:

Replying to stumpc5:

Replying to dimpase:

You can read the code of I.hilbert_series to see that it takes I.hilbert_numerator(), which is computed by Singular (again, read the code), and divides it by (t-1)P.ngens().

I actually did read the code, but

  • I only see "Singular is doing the computation" in sage.libs.singular.function_factory.ff.hilb,
  • the output is then going to Sage in a way I don't know, and finally
  • there is a polynomial division.

So a priori, the bug can be in either of these parts.

Anyway, it seems to be a bug in Singular, so I suggest to leave this ticket open until that is fixed, and then add a doctest for checking...

I would also appreciate if someone knowing how this Hilbert numerator is computed in Singular would send a bug report there (or shows me how to reproduce the bug in Singular).

comment:8 in reply to: ↑ 7 ; follow-up: Changed 4 years ago by dimpase

  • Report Upstream changed from N/A to Not yet reported upstream; Will do shortly.

Replying to stumpc5:

Anyway, it seems to be a bug in Singular, so I suggest to leave this ticket open until that is fixed, and then add a doctest for checking...

I would also appreciate if someone knowing how this Hilbert numerator is computed in Singular would send a bug report there (or shows me how to reproduce the bug in Singular).

I think you know that in the output of hilb(std(I)) (the part you omitted and indicated by ...) is the numerator that should be divisible by the appropriate power of t-1, but it is not.

As to writing up the corresponding Singular commands to show this... I can do this for you, although I didn't touch it for 5+ years...

comment:9 in reply to: ↑ 8 ; follow-up: Changed 4 years ago by stumpc5

Replying to dimpase:

Replying to stumpc5: I think you know that in the output of hilb(std(I)) (the part you omitted and indicated by ...) is the numerator that should be divisible by the appropriate power of t-1, but it is not.

Oh, I didn't know, thanks. If that is the not correct Singular output, how is the stuff below

//         1 t^0
//        30 t^1
//       135 t^2
//       120 t^3
// dimension (proj.)  = 13

then again correct (saying that it gives the numerator divided by the maximal power of (t-1)?

comment:10 in reply to: ↑ 9 ; follow-up: Changed 4 years ago by dimpase

Replying to stumpc5:

Replying to dimpase:

Replying to stumpc5: I think you know that in the output of hilb(std(I)) (the part you omitted and indicated by ...) is the numerator that should be divisible by the appropriate power of t-1, but it is not.

Oh, I didn't know, thanks. If that is the not correct Singular output, how is the stuff below

//         1 t^0
//        30 t^1
//       135 t^2
//       120 t^3
// dimension (proj.)  = 13

then again correct (saying that it gives the numerator divided by the maximal power of (t-1)?

It's all a bit Greek to me now, and https://www.singular.uni-kl.de/Manual/latest/sing_310.htm#SEC349 does not help much.

Perhaps there is more than one way to compute this info, and Singular does something leading to the correct answer, even though the data it supplies to Sage is not correct.

You seem to know a reference for the correct answers to this question for all m, no? This would be handy in the bug report upstream.

comment:11 in reply to: ↑ 10 Changed 4 years ago by stumpc5

Replying to dimpase:

You seem to know a reference for the correct answers to this question for all m, no? This would be handy in the bug report upstream.

This Hilbert series is known (though maybe implicit as some h vector of a simplicial complex), but I don't have a reference at hand. I will try to get one.

The reason it seemed suspicious is that I conjecture a term order for which its initial ideal is the Stanley-Reisner ideal of a simplicial sphere I was studying for completely different reasons. And this example was spit out to contradict that conjecture. But, as often, it was the code that's buggy and and not the conjectural description. That's anyway not yet in a form to be made public.

comment:12 Changed 4 years ago by dimpase

One would try to compute the same with Macaulay2. Did you try this?

comment:13 Changed 4 years ago by stumpc5

Yes, I did -- here it is the output again (with a few long outputs ...'ed):

i31 : R = QQ[x11,x12,x13,x14,x15,x16,x17,x18,x19,x110,x111,x21,x22,x23,x24,x25,x26,x27,x28,x29,x210,x211,x31,x32,x33,x34,x35,x36,x37,x38,x39,x310,x311,x41,x42,x43,x44,x45,x46,x47,x48,x49,x410,x411]

o31 = R

o31 : PolynomialRing

i32 : M = matrix{{x11,x12,x13,x14,x15,x16,x17,x18,x19,x110,x111},{x21,x22,x23,x24,x25,x26,x27,x28,x29,x210,x211},{x31,x32,x33,x34,x35,x36,x37,x38,x39,x310,x311},{x41,x42,x43,x44,x45,x46,x47,x48,x49,x410,x411}}

o32 = | x11 x12 x13 x14 x15 x16 x17 x18 x19 x110 x111 |
      | x21 x22 x23 x24 x25 x26 x27 x28 x29 x210 x211 |
      | x31 x32 x33 x34 x35 x36 x37 x38 x39 x310 x311 |
      | x41 x42 x43 x44 x45 x46 x47 x48 x49 x410 x411 |

              4       11
o32 : Matrix R  <--- R

i33 : I = minors(2,M) 

o33 = ...
o33 : Ideal of R

i34 : s = hilbertSeries (R/I)

o34 = ...
o34 : Expression of class Divide

i35 : denominator s

             44
o35 = (1 - T)

o35 : Expression of class Product

i36 : n = numerator s

o36 = ...
o36 : ZZ[T]

i37 : factor n

               30               2       3
o37 = (- 1 + T)  (1 + 30T + 135T  + 120T )

o37 : Expression of class Product

comment:14 Changed 4 years ago by dimpase

by the way, here is what they mean by the 1st and by the 2nd Hilbert series. https://www.singular.uni-kl.de/DEMO_HTML/Examples/Hilbert/index.html

Indeed, this alone seems enough for a bug report. (save the Singular code :-))

comment:15 Changed 4 years ago by stumpc5

Okay, I am sending a bug report to Singular now...

comment:16 Changed 4 years ago by dimpase

any progress on this?

comment:17 Changed 4 years ago by stumpc5

You find the bug report from Feb at https://www.singular.uni-kl.de:8005/trac/ticket/750, but no one seems to have looked at it so far...

comment:18 Changed 4 years ago by dimpase

  • Report Upstream changed from Not yet reported upstream; Will do shortly. to Reported upstream. No feedback yet.

comment:19 Changed 2 years ago by stumpc5

As rechecked at http://www.singular.uni-kl.de:8002/trac/ticket/750#comment:1, this appears to be an interface bug, and no bug directly in singular.

comment:20 Changed 2 years ago by chapoton

Looks like an overflow error

sage: QQ['z'](hilb(I, 2))*(1-z)**30-QQ['z'](hilb(I, 1))
// ** _ is no standard basis
// ** _ is no standard basis
-4294967296*z^22 + 4294967296*z^21 - 4294967296*z^20 + 4294967296*z^19 - 4294967296*z^18
sage: factor(4294967296)
2^32

comment:21 Changed 2 years ago by stumpc5

On the singular ticket they observe that what is shown here as the 2nd Hilbert series is actually their 1st Hilbert series.

comment:22 Changed 2 years ago by chapoton

But the long singular answer should factorise with a factor (1-t)^30, no ? And it does not!

Last edited 2 years ago by chapoton (previous) (diff)

comment:23 Changed 2 years ago by chapoton

So (I double checked) this is a bug in singular indeed. Someone should reopen ​http://www.singular.uni-kl.de:8002/trac/ticket/750

comment:24 Changed 2 years ago by chapoton

More precisely, the difference between the good result and the singular result is

(2^32) * t^18 * (t^4 - t^3 + t^2 - t + 1)

in singular 4.1.0 (which is the version used by sage and available online)

comment:25 Changed 2 years ago by stumpc5

Okay, I can follow your computation. Is this problem now just a coincidence, or how should it be related to the bug I reported?

The bug I had was that the Sage Hilbert series is the singular second Hilbert series. This is given correctly by the computation 4x10, while the (as you say wrongly computed) first Hilbert series was given for the 4x11 example.

comment:26 Changed 2 years ago by chapoton

Sage only ever calls the "first Hilbert series". If the result of singular for hilb(I, 1) was correct, the rational function in sage would simplify and its numerator would be the same result as the "second Hilbert series" of singular. *The problem is a bug in Singular*.

Now we could "fix" this bug by calling instead hilb(I, 2) and changing the sage code, but I would not call that a good solution.

comment:27 Changed 2 years ago by stumpc5

oh, I see, sorry for my ignorance.

comment:28 Changed 2 years ago by chapoton

  • Report Upstream changed from Reported upstream. No feedback yet. to Reported upstream. Developers acknowledge bug.

comment:29 Changed 2 years ago by chapoton

  • Description modified (diff)

comment:30 follow-up: Changed 19 months ago by SimonKing

Note that in the current version 4.1.1p2.p0 of Singular, there is not a wrong result due to integer overflow, but an honest "int overflow error" is raised.

In my computations of group cohomology, I suffered from these overflow errors as well, which is why I opened #26243 and provided an implementation of Hilbert series computations that goes beyond Singular's limitations.

My implementation is not much slower than Singular in the examples in which Singular works. However, in #26243, I do not propose to use my implementation as default way to compute Hilbert series (in particular, I didn't touch the .hilbert_numerator() method of multipolynomial ideals). However, if you consider to use it in the .hilbert_*() methods by default in order to fix *this* ticket: Go ahead!

EDIT: The failing example from here became doctest for the code from #26243.

Last edited 19 months ago by SimonKing (previous) (diff)

comment:31 follow-up: Changed 19 months ago by SimonKing

The code in #25243 has improved, it now beats Singular also speedwise.

So, I suggest to use #25243 here.

comment:32 in reply to: ↑ 31 Changed 19 months ago by tscrim

Replying to SimonKing:

The code in #25243 has improved, it now beats Singular also speedwise.

So, I suggest to use #25243 here.

Typo: #26243 :P

comment:33 Changed 19 months ago by SimonKing

  • Dependencies set to #26243

comment:34 in reply to: ↑ 30 Changed 19 months ago by john_perry

Replying to SimonKing:

However, in #26243, I do not propose to use my implementation as default way to compute Hilbert series (in particular, I didn't touch the .hilbert_numerator() method of multipolynomial ideals). However, if you consider to use it in the .hilbert_*() methods by default in order to fix *this* ticket: Go ahead!

+1 to this. I may work on it this coming week if no one else volunteers, but be advised that various obligations mean I'm slow.

comment:35 follow-up: Changed 19 months ago by SimonKing

If I see that correctly, the following functions need to be taken care of (all in sage.rings.polynomial.multi_polynomial_ideal):

  • hilbert_series: Here we can use what I called hilbert_poincare_series in sage.rings.hilbert (name according to both the book of Greuel and Pfister, and to the article of Bigatti).
  • hilbert_numerator: Here we can use what I called first_hilbert_series in sage.rings.hilbert (name according to the book of Greuel and Pfister).
  • hilbert_polynomial: This needs to be re-implemented.

If I understand correctly, one wouldn't provide degree weights for hilbert_polynomial, and I think it is easy enough to follow Definition 5.1.4 in the Greuel-Pfister book.

comment:36 in reply to: ↑ 35 ; follow-up: Changed 19 months ago by john_perry

Replying to SimonKing:

If I see that correctly, the following functions need to be taken care of (all in sage.rings.polynomial.multi_polynomial_ideal):

  • hilbert_series: Here we can use what I called hilbert_poincare_series in sage.rings.hilbert (name according to both the book of Greuel and Pfister, and to the article of Bigatti).
  • hilbert_numerator: Here we can use what I called first_hilbert_series in sage.rings.hilbert (name according to the book of Greuel and Pfister).

Thanks. I'll aim for that.

  • hilbert_polynomial: This needs to be re-implemented.

I have done that before, so I should be able to do it for this. :-)

If I understand correctly, one wouldn't provide degree weights for hilbert_polynomial...

Yes, the Hilbert series is non-polynomial for non-standard weights, and I understand that finding formulas for this is an area of research that no one is actually involved in, unless someone extended the work of Massimo Caboara's Master's student about 4-5 years ago.

comment:37 in reply to: ↑ 36 Changed 19 months ago by SimonKing

Replying to john_perry:

Replying to SimonKing:

If I see that correctly, the following functions need to be taken care of (all in sage.rings.polynomial.multi_polynomial_ideal):

  • hilbert_series: Here we can use what I called hilbert_poincare_series in sage.rings.hilbert (name according to both the book of Greuel and Pfister, and to the article of Bigatti).
  • hilbert_numerator: Here we can use what I called first_hilbert_series in sage.rings.hilbert (name according to the book of Greuel and Pfister).

Thanks. I'll aim for that.

  • hilbert_polynomial: This needs to be re-implemented.

I have done that before, so I should be able to do it for this. :-)

I'm almost done and will probably submit it later today. But of course I'd be glad about a review.

Last edited 19 months ago by SimonKing (previous) (diff)

comment:38 Changed 19 months ago by SimonKing

  • Branch set to u/SimonKing/Hilbert_series_bug

comment:39 Changed 19 months ago by SimonKing

  • Authors set to Simon King
  • Commit set to dcc3f1709927a1ea622e42b38543e42ed3fb0f16
  • Status changed from new to needs_review

With the new commit, the stuff from #26243 is used for all Hilbert series computations. The tests in sage/rings/polynomial pass, so, I guess it is "needs review" (I trust that the patchbot will run full tests).

For the record, here is a timing that comes out of some of my cohomology computations, in a size that Singular can manage:

sage: I
Ideal (a_1_0^2, a_1_0*b_1_1, b_2_2*a_1_0, b_2_3*a_1_0, b_2_1*b_1_1, a_1_0*b_3_4, a_1_0*b_3_5, a_1_0*b_3_6, b_2_1*b_2_2, b_1_1*b_3_5, b_4_10*a_1_0, b_2_1*b_3_4, b_2_1*b_3_6, b_2_2*b_3_5, b_4_10*b_1_1 + b_2_3*b_3_4 + b_2_3^2*b_1_1 + b_2_2*b_3_6, a_1_0*b_5_17, b_2_1*b_4_10 + b_2_1*b_2_3^2, b_3_4^2 + b_1_1^3*b_3_6 + b_2_2*b_1_1*b_3_4 + b_2_2*b_2_3*b_1_1^2 + b_2_2^3, b_3_4*b_3_5, b_3_4*b_3_6 + b_1_1*b_5_17 + b_2_3*b_1_1*b_3_6 + b_2_3*b_1_1*b_3_4 + b_2_2^2*b_2_3, b_3_5^2 + b_2_1*b_2_3^2, b_3_5*b_3_6, b_3_6^2 + b_2_3*b_1_1*b_3_6 + b_2_2*b_2_3^2 + c_4_11*b_1_1^2, b_2_1*b_5_17, b_4_10*b_3_4 + b_2_3*b_1_1^2*b_3_6 + b_2_3^2*b_3_4 + b_2_2*b_5_17 + b_2_2*b_2_3*b_3_6 + b_2_2*b_2_3^2*b_1_1, b_4_10*b_3_5 + b_2_3^2*b_3_5, b_4_10*b_3_6 + b_2_3*b_5_17 + b_2_3^2*b_3_4 + b_2_2*b_2_3*b_3_6 + b_2_2*c_4_11*b_1_1, b_4_10^2 + b_2_3^2*b_1_1*b_3_6 + b_2_3^4 + b_2_2*b_2_3*b_4_10 + b_2_2^2*c_4_11, b_3_4*b_5_17 + b_2_3*b_1_1*b_5_17 + b_2_3^2*b_1_1*b_3_6 + b_2_3^2*b_1_1*b_3_4 + b_2_2*b_1_1*b_5_17 + b_2_2^2*b_4_10 + c_4_11*b_1_1^4, b_3_5*b_5_17, b_3_6*b_5_17 + b_2_3^2*b_1_1*b_3_6 + b_2_2*b_2_3*b_4_10 + c_4_11*b_1_1*b_3_4 + b_2_3*c_4_11*b_1_1^2, b_4_10*b_5_17 + b_2_3^3*b_3_6 + b_2_3^3*b_3_4 + b_2_2*b_2_3*b_5_17 + b_2_2*b_2_3^2*b_3_6 + b_2_3*c_4_11*b_1_1^3 + b_2_2*c_4_11*b_3_4 + b_2_2*b_2_3*c_4_11*b_1_1, b_5_17^2 + b_2_3^3*b_1_1*b_3_6 + b_2_2*b_2_3*b_1_1*b_5_17 + b_2_2*b_2_3^2*b_1_1*b_3_6 + b_2_2*b_2_3^4 + b_2_2^2*b_2_3*b_4_10 + c_4_11*b_1_1^3*b_3_6 + b_2_3*c_4_11*b_1_1^4 + b_2_3^2*c_4_11*b_1_1^2 + b_2_2*c_4_11*b_1_1*b_3_4 + b_2_2*b_2_3*c_4_11*b_1_1^2 + b_2_2^3*c_4_11) of Multivariate Polynomial Ring in b_2_1, b_2_2, b_2_3, b_4_10, c_4_11, a_1_0, b_1_1, b_3_4, b_3_5, b_3_6, b_5_17 over Finite Field of size 2
sage: G = I.groebner_basis()    # caching the Gröbner bases

Now, with the new commit:

sage: %timeit I.hilbert_polynomial()
The slowest run took 9.05 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 686 µs per loop
sage: I.hilbert_polynomial()
2/3*t^3 + 4*t^2 + 16/3*t + 1

Without the new commit (i.e., based on using Singular):

sage: %timeit I.hilbert_polynomial()
The slowest run took 8.71 times longer than the fastest. This could mean that an intermediate result is being cached.
100 loops, best of 3: 2.33 ms per loop
sage: I.hilbert_polynomial()
2/3*t^3 + 4*t^2 + 16/3*t + 1

Last 10 new commits:

7698bb6Some code optimizations in Hilbert series computation
8ae070bUse flint polynomial boilerplate for Hilbert series computation
d997abeFix a corner case in Hilbert series computation
5af3439Fix the erroneous fix of some corner cases
138f85bTurn some functions of polynomial.hilbert to methods of ETuple
e900222Hilbert series: Remove unnecessary bound checks, simplify code branches
f0e5c85Add documentation to new ETuple methods
2d29b9bUsing slices, custom median, and formatting/P#P8 tweaks to Hilbert code.
07ff2e1Reverting use of get_exp for __getitem__.
dcc3f17Do not use Singular for Hilbert series computations

comment:40 follow-up: Changed 19 months ago by dimpase

Impressive speedup! Where does it come from? A better algorithm?

comment:41 in reply to: ↑ 40 Changed 19 months ago by SimonKing

Replying to dimpase:

Impressive speedup! Where does it come from? A better algorithm?

I am not totally sure. We have to distinguish: hilbert_numerator and hilbert_series on the one hand, and hilbert_polynomial on the other hand.

The computation of "first Hilbert series" (also known as "Hilbert numerator") and "Hilbert-Poincaré series" (also known as "Hilbert series") relies on a new implementation introduced at #26243. The original motivation was to allow for really big computations (>70 variables, >2000 ideal generators) that would fail in Singular because of its 32-bit integer limitation. But in the end, it was also about speed: The new implementation beats Singular in medium size examples.

The algorithm for the Hilbert-Poincaré series isn't exactly new, but differs from what is done in Singular as follows. Let I be a monomial ideal in a multivariate polynomial Ring R and let m be a monomial of degree d. Then, I think all implementations for the computation of the Hilbert-Poincaré series HP(R/I) rely on the equation HP(R/I) = HP(R/(I+m))+t^d*HP(R/(I:(m))). The differences lie in the choice of m. Some systems choose m as one of the generators of I and read the above equation backwards: HP(R/(I+m)) = HP(R/I) - t^d*HP(R/(I:(m))). Singular chooses m to be a variable that appears in at least one generator of degree >1 of I. I think CoCoA uses the gcd of two randomly chosen generators. My implementation chooses m to be the variable that appears most frequently in the generators of I, to the power of the median of all exponents of that variable in the generators of I.

I was told by John Perry that this choice is known in the literature, but I don't have a reference. Anyway, I think part of the speedup is due to good data structures and careful usage of Cython.

However, all of the above doesn't apply here. With the example of the above timing for hilbert_polynomial, the timings for hilbert_numerator and hilbert_series have hardly changed. With Singular:

sage: %timeit I.hilbert_numerator()
The slowest run took 4.14 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 270 µs per loop
sage: %timeit I.hilbert_series()
The slowest run took 5.27 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 375 µs per loop

With the new implementation:

sage: %timeit I.hilbert_numerator()
1000 loops, best of 3: 280 µs per loop
sage: %timeit I.hilbert_series()
The slowest run took 5.73 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 367 µs per loop

So, in this example, the change really is in the computation of hilbert_polynomial. Since my code is just straight forward, I really don't know why the old code was so much slower.

comment:42 Changed 19 months ago by git

  • Commit changed from dcc3f1709927a1ea622e42b38543e42ed3fb0f16 to 631824276e86028e86931364c97b46f88e325eff

Branch pushed to git repo; I updated commit sha1. New commits:

6318242Fix a sign in Hilbert polynomial computation

comment:43 Changed 19 months ago by SimonKing

By playing with a bigger example (where I is the ideal defined in the attachment to #26243), I found that I got the sign wrong, meaning that depending on the example I got the correct polynomial or its negative. The new commit fixes it.

With the new commit and the example from #26243, I get consistent results (i.e., the Hilbert polynomial seems to return the coefficients of the Hilbert-Poincaré series for all but finitely many degrees):

sage: P = PowerSeriesRing(QQ,'t',default_prec=50)
sage: hp = I.hilbert_series()
sage: hs = P(hp.numerator())/P(hp.denominator())
sage: hs
1 + 77*t + 436*t^2 + 1396*t^3 + 3446*t^4 + 7194*t^5 + 13376*t^6 + 22856*t^7 + 36626*t^8 + 55806*t^9 + 81644*t^10 + 115516*t^11 + 158926*t^12 + 213506*t^13 + 281016*t^14 + 363344*t^15 + 462506*t^16 + 580646*t^17 + 720036*t^18 + 883076*t^19 + 1072294*t^20 + 1290346*t^21 + 1540016*t^22 + 1824216*t^23 + 2145986*t^24 + 2508494*t^25 + 2915036*t^26 + 3369036*t^27 + 3874046*t^28 + 4433746*t^29 + 5051944*t^30 + 5732576*t^31 + 6479706*t^32 + 7297526*t^33 + 8190356*t^34 + 9162644*t^35 + 10218966*t^36 + 11364026*t^37 + 12602656*t^38 + 13939816*t^39 + 15380594*t^40 + 16930206*t^41 + 18593996*t^42 + 20377436*t^43 + 22286126*t^44 + 24325794*t^45 + 26502296*t^46 + 28821616*t^47 + 31289866*t^48 + 33913286*t^49 + O(t^50)
sage: p = I.hilbert_polynomial()
sage: p(t=1)  # the only exception
86
sage: p(t=2)
436
sage: p(t=3)
1396
sage: p(t=4)
3446
sage: p(t=5)
7194
sage: p(t=20)
1072294
sage: p(t=40)
15380594

comment:44 Changed 19 months ago by git

  • Commit changed from 631824276e86028e86931364c97b46f88e325eff to ae4a2fc0526e9e9a88bf010fd9aceacdefed8b2b

Branch pushed to git repo; I updated commit sha1. New commits:

ae4a2fcAdd consistency check for Hilbert polynomial and fix a corner case

comment:45 Changed 19 months ago by SimonKing

It was needed to special-case zero-dimensional ideals. Problem was: In that case, the pole order of the Hilbert-Poincaré series at t=1 is zero, which means one was dividing by the factorial of -1, which gave an error. In fact, if the Hilbert-Poincaré series is a polynomial, then all but finitely many slices of R/I are of dimension zero, which means that the Hilbert polynomial is zero.

While I was at it: In the first commit on this ticket, I thought I needed to amend the sign, but it turned out that I don't need. Just to be on the safe side, I am now asserting that the leading coefficient of the Hilbert polynomial is positive (it has to be, because if t is sufficiently large then the Hilbert polynomial is a natural number).

comment:46 follow-up: Changed 19 months ago by tscrim

Some micro-optimizations: if len(denom): -> if denom:, you could use for n, coeff in enumerate(second_hilbert), and you do not need to pass a list to sum and prod (which is slower because it has to create and populate the list).

Some stylistic things: PEP8 spaces are my friend and can be your friend too :P, and error messages in Python are causes (no initial capital letter and no period/full-stop).

I think it would be good to still have a way to access the Singular implementation through an algorithm keyword (of course, defaulting to 'sage'). (Is this there, and I missed it?) That way you can implement more robust tests by comparing the results and perhaps there are some cases where Singular does win that someone may need frequently.

comment:47 in reply to: ↑ 46 ; follow-up: Changed 19 months ago by SimonKing

Replying to tscrim:

Some micro-optimizations: if len(denom): -> if denom:,

OK. I was not sure if a factorisation object behaves the same as a list.

I think it would be good to still have a way to access the Singular implementation through an algorithm keyword (of course, defaulting to 'sage'). (Is this there, and I missed it?) That way you can implement more robust tests by comparing the results and perhaps there are some cases where Singular does win that someone may need frequently.

That's a very good idea! So, I should merge in the old implementation.

comment:48 in reply to: ↑ 47 Changed 19 months ago by tscrim

Replying to SimonKing:

Replying to tscrim:

Some micro-optimizations: if len(denom): -> if denom:,

OK. I was not sure if a factorisation object behaves the same as a list.

This is true in basically every case I've come across:

sage: f = factor(123521321)
sage: f
7 * 11 * 997 * 1609
sage: %timeit bool(len(f))
The slowest run took 35.73 times longer than the fastest. This could mean that an intermediate result is being cached.
1000000 loops, best of 3: 334 ns per loop
sage: %timeit bool(f)
The slowest run took 34.95 times longer than the fastest. This could mean that an intermediate result is being cached.
1000000 loops, best of 3: 314 ns per loop

comment:49 Changed 19 months ago by git

  • Commit changed from ae4a2fc0526e9e9a88bf010fd9aceacdefed8b2b to 1505cceba6f52b43b00d92b09fa087a8eaff07fb

Branch pushed to git repo; I updated commit sha1. New commits:

1505cceKeep Singular as optional backend for Hilbert series computations; some PEP8

comment:50 Changed 19 months ago by SimonKing

I changed the (old!) error message to lower case and removed the full-stop. And I hope the style is now enough "spacey"...

Moreover, I added an optional parameter algorithm.

However, one thing is very very strange now. A doctest that previously worked is now still giving the correct result, but also showing an odd comment; this does NOT happen when running the test interactively:

File "src/sage/rings/polynomial/multi_polynomial_ideal.py", line 2506, in sage.rings.polynomial.multi_polynomial_ideal.?.hilbert_series
Failed example:
    I.hilbert_series()
Expected:
    (t^4 + t^3 + t^2 + t + 1)/(t^2 - 2*t + 1)
Got:
       skipping text from `parameter`
    (t^4 + t^3 + t^2 + t + 1)/(t^2 - 2*t + 1)

WTH??


New commits:

1505cceKeep Singular as optional backend for Hilbert series computations; some PEP8

comment:51 Changed 19 months ago by SimonKing

Apparently the comment is an artefact from a previous libsingular error:

sage: n = 4; m = 11; P = PolynomialRing(QQ, n * m, "x"); x = P.gens(); M = Matrix(n, x)
sage: Minors = P.ideal(M.minors(2))
sage: hp = Minors.hilbert_polynomial(); hp
1/21772800*t^13 + 61/21772800*t^12 + 1661/21772800*t^11 + 26681/21772800*t^10 + 93841/7257600*t^9 + 685421/7257600*t^8 + 1524809/3110400*t^7 + 39780323/21772800*t^6 + 6638071/1360800*t^5 + 12509761/1360800*t^4 + 2689031/226800*t^3 + 1494509/151200*t^2 + 12001/2520*t + 1
sage: Minors.hilbert_polynomial(algorithm = 'singular')
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
RuntimeError: Error raised calling singular function
Exception RuntimeError: RuntimeError('Error raised calling singular function',) in 'sage.libs.singular.function.LibraryCallHandler.handle_call' ignored
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-2-bfe31fe9a36b> in <module>()
----> 1 Minors.hilbert_polynomial(algorithm = 'singular')

/home/king/Sage/git/sage/local/lib/python2.7/site-packages/sage/rings/polynomial/multi_polynomial_ideal.pyc in __call__(self, *args, **kwds)
    295         if not R.base_ring().is_field():
    296             raise ValueError("Coefficient ring must be a field for function '%s'."%(self.f.__name__))
--> 297         return self.f(self._instance, *args, **kwds)
    298 
    299 require_field = RequireField

/home/king/Sage/git/sage/local/lib/python2.7/site-packages/sage/rings/polynomial/multi_polynomial_ideal.pyc in hilbert_polynomial(self, algorithm)
   2466             hilbPoly = sage.libs.singular.function_factory.ff.poly__lib.hilbPoly
   2467 
-> 2468             hp = hilbPoly(self)
   2469             t = ZZ['t'].gen()
   2470             fp = ZZ(len(hp)-1).factorial()

/home/king/Sage/git/sage/local/lib/python2.7/site-packages/sage/libs/singular/function.pyx in sage.libs.singular.function.SingularFunction.__call__ (build/cythonized/sage/libs/singular/function.cpp:15394)()
   1328         if not (isinstance(ring, MPolynomialRing_libsingular) or isinstance(ring, NCPolynomialRing_plural)):
   1329             raise TypeError("Cannot call Singular function '%s' with ring parameter of type '%s'"%(self._name,type(ring)))
-> 1330         return call_function(self, args, ring, interruptible, attributes)
   1331 
   1332     def _instancedoc_(self):

/home/king/Sage/git/sage/local/lib/python2.7/site-packages/sage/libs/singular/function.pyx in sage.libs.singular.function.call_function (build/cythonized/sage/libs/singular/function.cpp:17499)()
   1526     if errorreported:
   1527         errorreported = 0
-> 1528         raise RuntimeError("error in Singular function call %r:\n%s" %
   1529             (self._name, "\n".join(error_messages)))
   1530 

RuntimeError: error in Singular function call 'hilbPoly':
int overflow in hilb 1
error occurred in or before poly.lib::hilbPoly line 58: `   intvec v=hilb(I,2);`
expected intvec-expression. type 'help intvec;'
sage: P.<x,y,z> = PolynomialRing(QQ)
sage: I = Ideal([x^3*y^2 + 3*x^2*y^2*z + y^3*z^2 + z^5])
sage: I.hilbert_series()
   skipping text from `parameter`
(t^4 + t^3 + t^2 + t + 1)/(t^2 - 2*t + 1)
sage: P.<x,y,z> = PolynomialRing(QQ)
sage: I = Ideal([x^3*y^2 + 3*x^2*y^2*z + y^3*z^2 + z^5])
sage: I.hilbert_series()
(t^4 + t^3 + t^2 + t + 1)/(t^2 - 2*t + 1)

What shall I do to fix it?

comment:52 Changed 19 months ago by SimonKing

The easiest solution would be to remove the test that demonstrates the shortcomings of Singular (maybe that would also be better for my karma). Do you agree it can be removed?

comment:53 follow-up: Changed 19 months ago by tscrim

Rather than remove it, I would mark it as either # known bug (20145) or # not tested (with a quick explanation why it is not tested).

comment:54 in reply to: ↑ 53 ; follow-up: Changed 19 months ago by SimonKing

Replying to tscrim:

Rather than remove it, I would mark it as either # known bug (20145) or # not tested (with a quick explanation why it is not tested).

"not tested" might be an option. If I understand correctly, "known bug" would still mean that it is tested, right? That would be bad, because of the side-effect of the error.

In any case, I will open a new ticket so that the side-effect can be fixed. And I would prefer to remove the test for now and introduce it at the new ticket (showing that the side-effect is fixed).

comment:55 in reply to: ↑ 54 Changed 19 months ago by tscrim

Replying to SimonKing:

Replying to tscrim:

Rather than remove it, I would mark it as either # known bug (20145) or # not tested (with a quick explanation why it is not tested).

"not tested" might be an option. If I understand correctly, "known bug" would still mean that it is tested, right? That would be bad, because of the side-effect of the error.

I don't remember off-hand. However, it is simple enough to test.

In any case, I will open a new ticket so that the side-effect can be fixed. And I would prefer to remove the test for now and introduce it at the new ticket (showing that the side-effect is fixed).

I would prefer # not tested (unless that is what you meant by remove). At least that gives an example to the user where it doesn't work and the side effect is the only reason it would not be tested.

comment:56 Changed 19 months ago by SimonKing

See #26300. One possibility would be to wait and see if #26300 can easily be fixed (say, within one day). Then, #26300 can be a dependency for this ticket (when I created #26300, I did it the other way around).

comment:57 Changed 19 months ago by git

  • Commit changed from 1505cceba6f52b43b00d92b09fa087a8eaff07fb to abe9c18f3441e29ed8d45258642dda7784092a94

Branch pushed to git repo; I updated commit sha1. New commits:

abe9c18Remove a test because of a side-effect

comment:58 Changed 19 months ago by SimonKing

I hope that the last commit complies with your suggestion: The test is still present, but marked # not tested, and it is commented that this is because of a side-effect that is investigated at #26300.

comment:59 Changed 19 months ago by tscrim

  • Branch changed from u/SimonKing/Hilbert_series_bug to u/tscrim/hilbert_series_bug-20145
  • Commit changed from abe9c18f3441e29ed8d45258642dda7784092a94 to ea5a851616b8c84fda8032d9f9c1219e50d03a99
  • Milestone changed from sage-7.1 to sage-8.4
  • Reviewers set to Travis Scrimshaw

Yes, that is good. I made a few other tweaks to match PEP8 and our coding/doc conventions in Sage. If my changes are good, then positive review.


New commits:

aac2eedMerge branch 'u/SimonKing/Hilbert_series_bug' of git://trac.sagemath.org/sage into u/tscrim/hilbert_series_bug-20145
ea5a851A few reviewer tweaks.

comment:60 Changed 19 months ago by SimonKing

  • Status changed from needs_review to positive_review

I trust that the tests still pass (didn't check myself). Your changes look good to me. However, although I did read PEP8, I don't understand why you added blank spaces around some operators and removed blank spaces around others. And where can I read about doc formatting conventions of Sage (in particular when it is about optional parameters)?

comment:61 follow-up: Changed 19 months ago by tscrim

PEP8 allows you to not write spaces for higher priority operators. So like in the exponentiation, that is higher than the addition, and thus, I don't need the spaces around the **. Of course, that is my own stylistic preference showing. :) For this:

-    def hilbert_series(self, grading = None, algorithm = 'sage'):
+    def hilbert_series(self, grading=None, algorithm='sage'):
-                return self.hilbert_numerator(algorithm = 'singular') / (1 - t)**n
+                return self.hilbert_numerator(algorithm='singular') / (1 - t)**n

PEP8 says arguments should not have space around the =.

Some of these things are not directly specified in the Sage doc conventions, but done as a general accepted-ish convention. At least, an optional to me means it does not have a default value (in some appropriate sense of the term for a user).

comment:62 in reply to: ↑ 61 ; follow-up: Changed 19 months ago by SimonKing

Replying to tscrim:

PEP8 allows you to not write spaces for higher priority operators. So like in the exponentiation, that is higher than the addition, and thus, I don't need the spaces around the **. Of course, that is my own stylistic preference showing. :)

My intention (although I quite often do typos) is to put spaces around...

  • assignments of variables (but not of optional arguments): bla = foo(bar)
  • comparisons: if bla <= foo:
  • + and -
  • after (but not before) a comma or semicolon

and to not put spaces around

  • *, / and **, because they have higher priority than + and -: (2*x - 1)**4)
  • assignments of optional arguments: foo(bla, bar=4)

Actually I dislike (4 * x - 1) ** 3, I find it difficult to read.

Would the above comply to PEP8?

Last edited 19 months ago by SimonKing (previous) (diff)

comment:63 Changed 19 months ago by git

  • Commit changed from ea5a851616b8c84fda8032d9f9c1219e50d03a99 to a1545cc6c843fcdb85cbca3c85ca8ce49de28616
  • Status changed from positive_review to needs_review

Branch pushed to git repo; I updated commit sha1 and set ticket back to needs_review. New commits:

a8d02ddFixing doctest output in palp_database.py due to signs.
7b2ed50Merge branch 'u/tscrim/hilbert_functions-26243' of git://trac.sagemath.org/sage into u/tscrim/hilbert_functions-26243
5637e07Fixing INPUT:: to INPUT: in hilbert.pyx.
a1545ccMerge branch 'u/tscrim/hilbert_functions-26243' of git://trac.sagemath.org/sage into u/tscrim/hilbert_series_bug-20145

comment:64 in reply to: ↑ 62 Changed 19 months ago by tscrim

Replying to SimonKing:

Replying to tscrim:

PEP8 allows you to not write spaces for higher priority operators. So like in the exponentiation, that is higher than the addition, and thus, I don't need the spaces around the **. Of course, that is my own stylistic preference showing. :)

My intention (although I quite often do typos) is to put spaces around...

  • assignments of variables (but not of optional arguments): bla = foo(bar)
  • comparisons: if bla <= foo:
  • + and -
  • after (but not before) a comma or semicolon

and to not put spaces around

  • *, / and **, because they have higher priority than + and -: (2*x - 1)**4)
  • assignments of optional arguments: foo(bla, bar=4)

Actually I dislike (4 * x - 1) ** 3, I find it difficult to read.

Would the above comply to PEP8?

By my understanding, yes. Of course, I don't always follow it as PEP8 are more guidelines, but I do try to stay as close to it as I can.

comment:65 Changed 19 months ago by tscrim

  • Status changed from needs_review to positive_review

Last commit was a trivial fix of a doctest and merging in a fix of a bad block from #26243. I'm allowing myself to reset to positive.

comment:66 Changed 19 months ago by vbraun

  • Branch changed from u/tscrim/hilbert_series_bug-20145 to a1545cc6c843fcdb85cbca3c85ca8ce49de28616
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.