sumiya11 / Groebner.jl

Groebner bases in (almost) pure Julia
https://sumiya11.github.io/Groebner.jl/
GNU General Public License v2.0
66 stars 13 forks source link

Calculate transformation matrix #130

Open martinkjlarsson opened 8 months ago

martinkjlarsson commented 8 months ago

Can we add an option for calculating the transformation matrix/change of basis matrix? That is, the matrix that maps the generators of the ideal to the Gröbner basis (see for example Oscar.jl or Macaualay2).

I don't know how much work this is, but it would be helpful for what I'm trying to do.

sumiya11 commented 8 months ago

Hi @martinkjlarsson ,

This is perhaps a bit tricky implementation-wise, but should be possible to add. I will look into it.

Two questions for the moment:

martinkjlarsson commented 8 months ago

Hi @sumiya11, the reason I asked was that I'm currently using Oscar.jl, but it's quite a large package and also doesn't work on Windows (except WSL). Another option for me would be to work directly against Singular.jl.

To answer your questions:

sumiya11 commented 8 months ago

@martinkjlarsson thanks for the clarification.

I have implemented the first version of this. For now, you can do add Groebner#master in Pkg mode to try it.

An example:

using Groebner, AbstractAlgebra

R, (x, y, z) = polynomial_ring(GF(2^30 + 3), ["x", "y", "z"], internal_ordering=:degrevlex)

f = [x * y * z - 1, x * y + x * z + y * z, x + y + z]
g, m = Groebner.groebner_with_change_matrix(f)

@show m
3×3 Matrix{AbstractAlgebra.Generic.MPoly{AbstractAlgebra.GFElem{Int64}}}:
 0  0             1
 0  1073741826    y + z
 1  1073741826*z  z^2

where the following holds:

@assert isgroebner(g)
@assert m * f == g

groebner_with_change_matrix works over GF and QQ.

At the moment, it supports degree reverse lex. ordering (you can either create polynomial rings with internal_ordering=:degrevlex as in the example, or use Groebner.groebner_with_change_matrix(f, ordering=DegRevLex()))

sumiya11 commented 8 months ago

groebner_with_change_matrix has noticeable performance overhead over groebner on some examples. I have not compared it against Singular or Macaulay2 yet (perhaps will do..).

If it becomes a bottleneck in your code, please let me know, I think I know how to speed it up

martinkjlarsson commented 8 months ago

Sweet! You were super fast with this. Thanks a lot!

sumiya11 commented 8 months ago

Out of curiosity, do you use these matrices in some application ?

martinkjlarsson commented 8 months ago

I'm experimenting with developing a Julia package for solving polynomial equation systems. It uses Gröbner basis methods to reduce a system to an eigendecomposition problem. In any case, there is a need to express some polynomials in the original generators rather than the Gröbner basis. The transformation matrix deals with this mapping.

sumiya11 commented 8 months ago

Thank you for the explanation! Do you mean using something akin to the Stickelberger’s theorem ?

I am also interested in solving polynomial systems. If you would need help with Groebner.jl or if you would like to discuss your experiments, feel free to let me know.

BTW, in case it can be useful, there is SemialgebraicSets.jl, which uses multiplication matrices with the Schur decomposition for solving systems.

martinkjlarsson commented 8 months ago

The method is used for producing a solver for a family of polynomial systems. For every instance of this family, one could in principle produce multiplication matrices and find the roots, but this involves Gröbner basis calculations which are costly. However, by performing some offline calculations first, we can replace the Gröbner basis calculations with a linear system - the so-called elimination template. The resulting solver will just contain a linear system and an eigendecomposition.

If you're curious here is a recent paper, but I'd recommend this one first. More fleshed out reading can be found in this thesis from p. 29.

It's too early for me to share results or indeed promise a package, but I'll let you know if I need help with Groebner.jl. Thanks!