Attempts at fixing regression for the OFPModules

[HechtiDerLachs]

I do this PR to investigate #6088 and the aftermath of #6052. It seems that we need to tune the behavior of presentation and _presentation_graded further. The changes which I disable here, were introduced in #5990 in order for some surfaces code of @wdecker to work. However, in another project with Simon Felten this now led to a mysterious regression, which I do not fully understand.

To reproduce our example, do the following:

P, (t,) = polynomial_ring(QQ, [:t])
S, (x, y, z, u, w) = graded_polynomial_ring(P, [:x, :y, :z, :u, :w])
F0 = graded_free_module(S, [2 for _ in 1:10])
e = gens(F0)
rels1 = FreeModElem{MPolyDecRingElem{QQMPolyRingElem, AbstractAlgebra.Generic.MPoly{QQMPolyRingElem}}}[u*e[1] - w*e[2] + 2*z*e[6] + t*z*e[7], u*e[4] - w*e[5] + 2*y*e[6] + t*y*e[7], y*e[1] + w*e[3] - z*e[4], y*e[2] + u*e[3] - z*e[5], z*e[3] - w*e[5] + y*e[6], u*e[1] + z*e[6] + t*y*e[8], z*e[1] + w*e[6] + t*y*e[9], z*e[2] - u*e[6] - t*u*e[7] + t*y*e[10], 2*x*e[6] + t*x*e[7] + u*e[9] - w*e[10], x*e[1] + w*e[8] - z*e[9], x*e[2] + u*e[8] - z*e[10], x*e[3] + z*e[7] - y*e[8], x*e[4] + w*e[7] - y*e[9], x*e[5] + u*e[7] - y*e[10], x*e[6] - z*e[8] + u*e[9], -y*z*e[3] + y*u*e[4] + y^2*e[6] + 1//3*x*z*e[7] + (1//3*x*y - t*z^2)*e[8] - t*w^2*e[9] - t*u^2*e[10], t*z^2*e[3] + t*w^2*e[4] + t*u^2*e[5] + 2//3*y*z*e[7] + (2*x*z - 1//3*y^2)*e[8] - x*u*e[9] - x*w*e[10], -t*w^2*e[1] - t*u^2*e[2] + t*y^2*e[3] + 1//3*z^2*e[7] + (t*x^2 - 2//3*y*z)*e[8], t*z^2*e[1] + t*y^2*e[4] - 2*t*u^2*e[6] + (1//3*z*w - t^2*u^2)*e[7] + 1//3*y*w*e[8] + (t*x^2 - y*z)*e[9], t*z^2*e[2] + t*y^2*e[5] + 2*t*w^2*e[6] + (1//3*z*u + t^2*w^2)*e[7] + 1//3*y*u*e[8] + (t*x^2 - y*z)*e[10], (-1//3*t*x*y*z - 1//3*z^3 + 1//3*z*w*u)*e[7] + (-1//3*t*x*y^2 - 1//3*y*z^2 + 1//3*y*w*u)*e[8]]
M = quo_object(F0, sub_object(F0, rels1));

G0 = graded_free_module(S, [0])
e = gens(G0)
rels2 = FreeModElem{MPolyDecRingElem{QQMPolyRingElem, AbstractAlgebra.Generic.MPoly{QQMPolyRingElem}}}[(-t*x^3 + x*y*z - t*y^3 - t*z^3 - t*w^3 - t*u^3)*e[1], (-t*x*y - z^2 + w*u)*e[1]]
N = quo_object(G0, sub_object(G0, rels2));

@time hom(M, N)

On the current master I did not see this code finish. With the changes in this PR it takes only about 40s.

I tried to investigate this a bit and suspected that the denominators in the coefficients might be the cause for the difference in runtimes. But eliminating those denominators manually did not solve the problem. It seems that the seemingly redundant production of the presentation of an already presented module is leading to a significantly better choice of a generating set. And by significant I mean the difference between running for 40 seconds or running forever. That is odd to me and I wanted to put this out for other people to have a look and share their insights. @jankoboehm ? @wdecker ?

It would be nice if we could complement this with @wdecker 's example and try to find the compromise which hopefully resolves both.

[wdecker]

It would be nice if we could complement this with @wdecker 's example and try to find the compromise which hopefully resolves both.

I had no particular examples in mind. I just thought that recomputing the presentation in cases where the presentation can be directly read off from how the module is given is redundant (costly) and should be avoided.

Why this is bad for the subsequent computations in one special example should be carefully examined. The problem seems to be this computation

kDelta = kernel(delta)

in this functon

function hom(M::OFPModule, N::OFPModule; algorithm::Symbol=:maps)

from

Hom_and ext.jl.

[HechtiDerLachs]

Thank you for taking a look, @wdecker . Indeed, I spotted the same place where it hangs, but did not find a good reason to explain the different behaviors.

[HechtiDerLachs]

Quick update: The problem seems to not arise in positive characteristic. The following code works for me both on master and on this PR:

p = next_prime(32024)
kk = GF(p)
Pp, to_Pp = change_base_ring(c->kk(numerator(c))*inv(kk(denominator(c))), P)
Sp, to_Sp = change_base_ring(to_Pp, S)
Mp, to_Mp = change_base_ring(to_Sp, M)
Np, to_Np = change_base_ring(to_Sp, N)
@time hom(Mp, Np)

In this case it runs even faster on master. So it might actually be a problem about coefficient swell over QQ. And, hence, a matter of choices of representatives. As I already mentioned, I played around a bit to make smart, but cheap choices without recomputing the presentation, but no luck so far. Any experience/advice/expertise here would be greatly appreciated!

[JohnAAbbott]

I confirm that on master @HechtiDerLachs example was still running after about 600s CPU (vmem was 3-4 Gbyte). After commenting out the code from #5990 the same example finished in under 40s. Trying to investigate...

[JohnAAbbott]

Minor update I modified the input by multiplying some of the rels1 elements by 3 to clear denominators. The code (incl. the change from #5990) produced an answer in about 1700s (roughly 30 mins). The coefficients in the intermediate value kDelta (see file Hom_and_ext.jl near line 85) were noticeably larger than without the change from #5990; the reported time was about 200s.

[HechtiDerLachs]

Thanks for the update! I tried a similar thing to multiplying by 3, but did probably not wait the full half hour. Thanks for the report. Do I understand correctly that 200s were achieved on this PR in conjunction with multiplication by 3? Then this multiplication is still making things worse compared to simply applying Singular's implementation for a new presentation.

[JohnAAbbott]

Another brief update: I printed out the values of kDelta and did some minimal "analysis" of them.
With the revision in PR #5990 the printed value was about 940Mbytes; the largest coefficient has about 28000 digits.
Without PR #5990 the value was about 640Mbytes, and the largest coefficient has about 2600 digits (ten times smaller). Not all polynomials had sch large coefficients.

[JohnAAbbott]

A further update: I ran the modified input over a finite field (Native.GF(32003)); below are the salient results.
With PR 5990 printed form of kDelta occupies about 3.5Mbytes; reported time 23s (incl 44% compilation).
Without PR 5990 printed form of kDelta occupies about 12Mbytes; reported time 43s (incl 24% compilation).

Here the revision from PR #5990 leads to faster computation and a smaller result; obviously, over a finite field there is no problem of "coefficient growth".

[JohnAAbbott]

Here is a small example to show potential problems with computing syzygy modules (which is related to kernel computation):

julia> P, (x,y,z,t) = QQ["x","y","z","t"];
julia> LL = [7*x^2*y - 2*x^2 + 2*x*y^2 + 8*x*y*z + y + 5*z*t + 5*z + 8*t^3, x^2*y - x^2 + 5*x*y^2 - x*y*t - x*z*t + 2*x - 2*y*z^2 - z*t + 6*t^2, 4*x^3 + 8*x^2 - 7*x*y*z - x*t + 9*y^3 - 5*y^2*z + 5*y^2 + 2*y*z*t + 3*t^2 + 4*t];
julia> syzygy_generators(reverse(LL))
julia> syzygy_generators(LL)

The sizes of the results are wildly different -- I didn't measure computation time.
The example was found via a random search.

ADDENDUM a more extreme is the following:

L = [-x^2*y + 2*x*y^2 - 2*x*y*z + 2*x*y + x*z^2 - 3*y^3 + 9*y^2*z - 8*y^2*t + 4*y*z*t + y*z - 8*t^3 - 5*t, -7*x^3 - 7*x^2*y - 8*x^2 - 8*x*y*z + 9*x*z - 4*x*t - 3*y^2*t - 8*y*t + 9*y - 5*t, 7*x^2*y + 4*x^2 + 2*x*t - 3*y^2*t - 9*y*z^2 - 2*y*z*t - 8*y*z + y*t^2 + 9*z + 4*t^2 + t];

(somewhat surprisingly) The computation times did not differ greatly despite a significant difference in the sizes of the result.

[jankoboehm]

In the shortcut introduced by #5990 a check=false seems to be missing. But that seems to make no difference with the issue.

The problem is that the given presentation is used and not recomputed as before #5990, hence simplifications are not done leading to longer computation times in commands like kernel when using the presentation. So this is dependent on the example, in some cases the recomputation of the presentation is a problem (which is why #5990 was introduced by @wdecker to make computations finish), in other cases the simplification is the problem. Since the current state of the presentation code is returning a correct result (and is faster) it seems to be reasonable to actually keep it, and rather request simplification/recomputation explicitely where needed? I would still fix the check thing.

[wdecker]

check=false seems

In which line?

[jankoboehm]

In aug = hom(F, SQ, gens(SQ)). However this would only skip grading checks here, which can be expensive, but not in the current example. The other thing is to set generators_map_to_generators. Again, this can be relevant, but not in the current example.

[fingolfin]

suggestion:
add a way to take f:A->B and Hom(B,C) to produce Hom(A,C).
Similar for g:C->D should be able to get Hom(B,D) (or overall: Hom(A,C))

You need that anyway internally, this is just about putting it into nice helpers.

Some suggestions for syntax: f*Hom(B,C)*g or compose(f,Hom(B,C)) or ...

(I'll try to clean up this comment later)

[jankoboehm]

One other aspect: Simplify or not depends actually to a large extent on what you want to do with it afterwards. For example in a Hom computation it is reasonable to simplify. If simplification is expensive then the Hom would most likely be even more. So this cannot only depend on a heuristic on the module.