Broadcasting for in-place operations on polynomials (blog post)

Specific Domains Optimization (Mathematical)
1

Some of you may be interested in this post that I just published:

TL;DR: I used the new broadcasting feature (thanks @mbauman + others!) to get a significant performance optimization in a setting without axes and indices.

5 Likes
2

Defining an alternative iterating style for terms is interesting. In MultivatiatePolynomials one needs to call terms on a polynomial to get an iterator over the terms.

We also have an implementation of Buchberger algo in SemialgebraicSets. It would be interesting to see a benchmark comparison.

1 Like
3

Thanks for your pointer @blegat ! It was very helpful for me to see your code as it reminded me to sort my input vector before starting the algorithm. That gives a good speedup in some of my test cases. I also like how we both take advantage of Julia’s unicode support to name the function gröbner :slight_smile:

I tried to do a side-by-side comparison, but that’s quite hard as Semialgebraicsets (and TypedPolynomials, on which it is based) make floating point approximations whereas PolynomialRings works over BigInt (or theoretically over any field; I’ve used Galois fields too).

That may theoretically give incorrect results, and unfortunately that also happens in practice for even moderately sized examples. Below, I included a code listing where Semialgebraicsets concludes that the cyclic-6 ideal is (1) whereas PolynomialRings returns a vector of 45 polynomials (as does Singular, which I guess we can use as a referee in this matter).

If this topic interests you, would you be willing to give PolynomialRings a try for your use case? It comes with a few other nice features:

  • works over any base ring
  • can also return the transformation that produces the Groebner basis from the inputs (“lift” as Singular calls it – this is a rather slow operation unfortunately)
  • it is much, much faster for the basic operations. Example:
julia> @btime (a+b+c+d+e+f)^10;                   # exponentiation in PolynomialRings
  25.418 ms (574591 allocations: 10.84 MiB)

julia> @btime (u+v+w+x+y+z)^10;                   # exponentiation in TypedPolynomials (probably power by squaring?)
  598.072 ms (171926 allocations: 4.22 GiB)

julia> @btime prod(a+b+c+d+e+f for _=1:10);       # repeated product in PolynomialRings
  2.712 ms (43144 allocations: 1.28 MiB)

I’m super willing to help make it work for you, either by writing APIs that you need or by accepting pull requests. Here’s the introductory docs.

Thanks again for your feedback!

Code listing for cyclic-6:

julia> import TypedPolynomials; import PolynomialRings; import SemialgebraicSets;

julia> cyclic6(a,b,c,d,e,f) = [
           a*b*c*d*e*f + -1,
           a + b + c + d + e + f,
           a*b + b*c + c*d + d*e + a*f + e*f,
           a*b*c + b*c*d + c*d*e + a*b*f + a*e*f + d*e*f,
           a*b*c*d + b*c*d*e + a*b*c*f + a*b*e*f + a*d*e*f + c*d*e*f,
           a*b*c*d*e + a*b*c*d*f + a*b*c*e*f + a*b*d*e*f + a*c*d*e*f + b*c*d*e*f,
       ]
cyclic6 (generic function with 1 method)

julia> vars_sas = TypedPolynomials.@polyvar u v w x y z
(u, v, w, x, y, z)

julia> PolynomialRings.@ring! ℤ[a,b,c,d,e,f]
ℤ[a,b,c,d,e,f]

julia> using BenchmarkTools

julia> @btime G_pr = PolynomialRings.gröbner_basis(cyclic6(a,b,c,d,e,f));
  1.316 s (18214531 allocations: 494.00 MiB)

julia> @btime G_sas = SemialgebraicSets.gröbnerbasis(cyclic6(u,v,w,x,y,z));
  7.002 s (1258677 allocations: 819.86 MiB)

julia> length(G_pr)
45

julia> length(G_sas)
1

julia> G_pr[1:2]
2-element Array{ℤ[a,b,c,d,e,f],1}:
 a + b + c + d + e + f                                       
 -b^2 + -b*d + c*d + -b*e + d*e + -2*b*f + -c*f + -d*f + -f^2

julia> G_sas[1]
 1.0
4

SemialgebraicSets is not based on TypedPolynomials, it can be used with any polynomial implementation that implements the MultivariatePolynomials API, for instance you can also use it with DynamicPolynomials. Moreover, you can use both TypedPolynomials and DynamicPolynomials with BigInt too :slight_smile:

I’d love to ! I try to make my code independent from a specific polynomial implementation by using MultivariatePolynomials and let the user choose with polynomial implementation he wants to use with the algorithms I implement so if PolynomialRings implements the MultivariatePolynomials API, that would automatically make it usable :slight_smile: You can find a list of packages built on top of the MultivariatePolynomials API here

5

Cool! I’ll make it so and post here with an update as soon as that is done. Thanks in advance!

2 Likes