Issues with Computing Groebner Basis for High-Order Runge Kutta Systems with 8 Stages and Order 7

[jpcurbelo]

I'm having trouble computing the Groebner basis for high-order Runge Kutta systems. Specifically, I'm working with a problem that has 8 stages and order 7, and the system of equations is given in this file. However, when I try to compute the Groebner basis for this problem, it seems to require more than 1Tb of RAM, which is causing issues.

For lower-order problems, it seems to work just fine. However, for the higher-order problem, I'm not sure if there's some workaround in the inputs that I might be missing, or if this behavior is to be expected for this type of problem. If anyone has any suggestions or insights into this issue, I would greatly appreciate it. Thank you!

[sumiya11]

Hi @jpcurbelo , thanks for sharing this interesting system!
I'll have a look. At first glance, I think we don't compress monomials here, but we totally should

[sumiya11]

You mentioned that you also have smaller similar systems. If you expect those to have a "similar" groebner basis as this large one, feel free to share the small systems as well, it will perhaps help

[jpcurbelo]

Hi @sumiya11. Thanks for your quick response. In fact, the system for RK(s=8, p=7) doesn't have a solution - that's what I expected to obtain from this.

For example, RK(s=4, p=4) has a solution (multiple solutions, in fact) and I got this Grobner basis.

On the other hand, RK(s=6, p=6) doesn't have a solution, and the groebner(system) gives us {"groebner_basis":["1"]}. From what I've understood, this means no Groebner basis (no solution), right?

Please, let me know if I didn't make myself clear enough :)

Thanks again!

[jpcurbelo]

Hey @sumiya11. I was wondering if you have had the time to look at this issue (a significant amount of memory RAM is required for specific problems). Please, let me know if you think there is something I could (try to) do on my end to somehow improve the performance of the jobs I'm running. Many thanks for your efforts!

[sumiya11]

Hi @jpcurbelo. I have looked at the issue, the problem is that in F4 algorithm we don't limit the maximum number of critical pairs that can be handled at the same time, rather strangely though

I am at it ( #64 ), but haven't figured out the best fix yet. Will try again this weekend. Sorry for this taking so long

[sumiya11]

I think I messed up git, so the issue got closed automatically 😅 .
Reopening now

[sumiya11]

Hi @jpcurbelo , so there is good news and bad news:

• Good news is that groebner now has the maxpairs keyword argument, which can be used to limit the matrix size in the internal F4 implementation. It is available in the latest Groebner.jl version.
  For example, if system is RK(s=6, p=6) that you have kindly provided, then groebner(system, maxpairs=1000) should allocate a bit less than the baseline.
• Bad news is that the big RK(s=8, p=7) is still not tractable. A complication is that we don't know the optimal values for maxpairs (it's tricky to control the matrix size in F4), and for small maxpairs the computation expectedly becomes much-much slower.

On my machine, while computing RK(s=8, p=7) with maxpairs=500, the process allocates 100 GB of RAM before it gets killed.
I have experimented with other values for maxpairs, but seemingly to no avail 😞 .

The commands std and slimgb from Singular are not competitive on RK(s=6, p=6), so I haven't checked them on the bigger system.

The Groebner command in Maple reports Error, (in Groebner:-Basis) Segmentation Violation occurred in external routine after using up around 25 GB of RAM (I used this Maple script).

I will work on reducing the overall memory footprint of Groebner.jl this summer, but I think at the moment Groebner.jl unfortunately cannot tackle such systems

[sumiya11]

As a sidenote, if you do not care about the basis monomial ordering, you can explicitly set it to DegRevLex(). Usually, it is more efficient than the default one:

# system is RK(s=6, p=6)

@time groebner(system);
# 258.834883 seconds (661.72 M allocations: 154.943 GiB, 5.61% gc time)
# Max. RAM:   9100.656 MiB

@time groebner(system, ordering=DegRevLex());
# 28.632552 seconds (7.53 M allocations: 2.620 GiB, 2.63% gc time)
# Max. RAM:   1676.414 MiB

@time groebner(system, ordering=DegRevLex(), maxpairs=1000);
# 12.556571 seconds (6.62 M allocations: 1.900 GiB, 4.33% gc time)
# Max. RAM:   1574.133 MiB

[sumiya11]

I think there is a chance that incremental signature GB algorithms can deal with these problems (e.g., https://www.singular.uni-kl.de/Manual/4-0-3/sing_391.htm). I'll check if I can make it work via Singular.jl

[jpcurbelo]

Hi @sumiya11! Many thanks. Your efforts are highly appreciated. I will take a look at the updates and materials you mentioned. I'll keep my finger crossed so you can find time to work on reducing the memory footprint. Once again, thank you so much. Your package has already helped me with other problems as well.

[pogudingleb]

@jpcurbelo
I wonder if there is a simple way to describe how the system is formed?
The reason I am asking is that, from the shape of polynomials, it looks like the process involves some substitutions. In my experience, it is sometimes better not to perform all the substitution but keep extra variables and relations (for example, if you would like to replace c with a + b, you can just add c - a - b to the system).

[jpcurbelo]

Hi @pogudingleb. I don't think there is a simple way to describe how the systems are formed but you are right, there are key substitutions involved.
The systems are basically based on Butcher's Tableu (see the image) and, as we might see, the c terms are not in the systems based on the substitutions:
$$c_i = \sum_{i}^s a_{i,j}$$
That's a good hint. I'll try to generate the systems without the substitutions and will generate the systems without tl let you know how it went.