Exterior powers of finitely presented modules

[HechtiDerLachs]

This is a suggestions for how we could do exterior powers of modules. See the test file for how this can be used.

[fingolfin]

@jankoboehm it would be great if you could have a look at this and comment on it. We want to make sure that potential applications can actually use this interface; and if not, revise it. But for this of course we need to know what is needed / missing / ...

[HechtiDerLachs]

@fieker : Is this like you had intended?

[HechtiDerLachs]

@jankoboehm : Sorry, this was not intended to be here. It was a little play-around for me needed to take some steps towards implementing sheaf cohomology. But I guess this might be useful also for you and the graded modules in general. Would you like this to be put in a separate PR?

[HechtiDerLachs]

Since some people like performance comparisons:

julia> n = 12
12

julia> R, x = polynomial_ring(QQ, "x" => 1:n)
(Multivariate polynomial ring in 12 variables over QQ, QQMPolyRingElem[x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10], x[11], x[12]])

julia> F = FreeMod(R, n)
Free module of rank 12 over Multivariate polynomial ring in 12 variables over QQ

julia> v = sum(a*e for (a, e) in zip(x, gens(F)))
x[1]*e[1] + x[2]*e[2] + x[3]*e[3] + x[4]*e[4] + x[5]*e[5] + x[6]*e[6] + x[7]*e[7] + x[8]*e[8] + x[9]*e[9] + x[10]*e[10] + x[11]*e[11] + x[12]*e[12]

julia> K1 = koszul_complex(x)
K1_-1 <---- K1_0 <---- K1_1 <---- K1_2 <---- K1_3 <---- K1_4 <---- K1_5 <---- K1_6 <---- K1_7 <---- K1_8 <---- K1_9 <---- K1_10 <---- K1_11 <---- K1_12 <---- K1_13

julia> K2 = koszul_complex(v)
K2_-1 <---- K2_0 <---- K2_1 <---- K2_2 <---- K2_3 <---- K2_4 <---- K2_5 <---- K2_6 <---- K2_7 <---- K2_8 <---- K2_9 <---- K2_10 <---- K2_11 <---- K2_12 <---- K2_13

julia> @time koszul_complex(x) # the Singular implementation
  7.038626 seconds (23.58 M allocations: 1.750 GiB, 26.06% gc time)
C_-1 <---- C_0 <---- C_1 <---- C_2 <---- C_3 <---- C_4 <---- C_5 <---- C_6 <---- C_7 <---- C_8 <---- C_9 <---- C_10 <---- C_11 <---- C_12 <---- C_13

julia> @time koszul_complex(v) # the new OSCAR implementation
  1.954575 seconds (5.20 M allocations: 278.098 MiB, 31.56% gc time)
C_-1 <---- C_0 <---- C_1 <---- C_2 <---- C_3 <---- C_4 <---- C_5 <---- C_6 <---- C_7 <---- C_8 <---- C_9 <---- C_10 <---- C_11 <---- C_12 <---- C_13

[jankoboehm]

It looks good to me, shall we call the function mult_map rather pure and the inverse inv_pure as in the tensor product? See also in the tensor product code:

set_attribute!(F, :tensor_pure_function => pure, :tensor_generator_decompose_function => inv_pure)

Note that these functions are created, when creating the module, so they do not create a parent on the fly.

Perhaps add a show function similar to the one for the tensor product of modules, see:

set_attribute!(F, :show => Hecke.show_tensor_product, :tensor_product => G)

With regard to the show function of the elements, you can set the field F.S in a free module F, which should have the symbols for the generators.

We need not only the case of free modules but also subquos. That we can do in a seperate effort.

[lgoettgens]

I think it would make sense to put OrderedMultiIndex into src/Combinatorics/ as there surely are other places where this can be used. What do you think?

[HechtiDerLachs]

I think it would make sense to put OrderedMultiIndex into src/Combinatorics/ as there surely are other places where this can be used. What do you think?

Agreed and done.

@lgoettgens : There seems to be some problem with an existing function is_exterior_power from the Lie algebras. I don't know what's the problem and what should be done to resolve this. Could you maybe have a look?

[lgoettgens]

@lgoettgens : There seems to be some problem with an existing function is_exterior_power from the Lie algebras. I don't know what's the problem and what should be done to resolve this. Could you maybe have a look?

I think you create a function with that name in the global Oscar module, right? In this case, to make the two same-named functions the same, you would need to remove the two exports in experimental/LieAlgebras/src/LieAlgebras.jl and instead add is_exterior_power to the long import list at the top of the file.

[lgoettgens]

A similar function already exists for LieAlgebraModule.jl under the name hom_exterior_power. It would be great to have a unified interface for both types of modules.
This is not a request for you to change your proposal, I just wanted to mention it.
@fieker suggested to name this method simply hom and supply both exterior powers as arguments, as well as the home on the base modules

[HechtiDerLachs]

If anyone found the time to approve and merge this, that would be great! I'm afraid the green lights of the tests might be corrupted by future merge conflicts.

[fingolfin]

I started to look at this but have to run to other appointments, so just some quick comments, I hope I can look at the rest later.

[fingolfin]

Why is this struct mutable?

[HechtiDerLachs]

Good question. It probably shouldn't be, unless we want to allow in-place iteration. Do we? Is there a potential speedup from making it immutable?

[fingolfin]

This makes this function _mult not type stable. Maybe

Suggested change

	any(x->(x in indices(b)), indices(a)) && return 0, a
	any(x->(x in indices(b)), indices(a)) && return 0, indices(a)

or (IMHO nicer)

Suggested change

	any(x->(x in indices(b)), indices(a)) && return 0, a
	any(x->(x in indices(b)), indices(a)) && return 0, T[]

[HechtiDerLachs]

Thanks for the catch. But the second needs to reallocate memory, doesn't it?

[fingolfin]

So... this is essentially the parent of OrderedMultiIndex{T}, except that it does not specify T...? Or "hardcodes" T to be Int (?)

[HechtiDerLachs]

I was not so sure about this, to be honest. It's not actually a parent in the Oscar sense, but in spirit: yes. And it's true: It doesn't know about the T. But it's hardcoded only in the sense that iterating over this will always produce an OrderedMultiIndex with T = Int64.

At some point I thought of removing the T altogether, but then I didn't see any real problem in keeping it either. It's just a bit weird that OrderdMultiIndexSet doesn't know about it. But since it's not an actual parent: Do you see a real problem here?

[HechtiDerLachs]

As discussed with @wdecker: I deprecated some of the singular methods around koszul complexes etc. Speed tests revealed that singular's depth computation is significantly faster. The same might hold for koszul_homology. Looking at the code for koszul_complex again, I was wondering: Maybe importing from singular is only so slow because it allocates all the matrices for the maps?

Either way: All methods imported from singular are now still available as internal methods. In particular for the depth computation I still use it whenever the coefficient ring is a field. For the Koszul homology I'm not so sure what the best solution is: If we were using singular's version, then the user might be wondering why the output is not compatible with the output of koszul_complex. For the moment I'm using the exclusive Oscar version of koszul_homology.

[wdecker]

@ederc, @hannes14: Can you see why importing the functions mentioned by @HechtiDerLachs are slower on the OSCAR side as compared to the Singular side? See homological-algebra.jl.

[HechtiDerLachs]

@wdecker: I suppose it's because of this line here. The morphisms of the Koszul complex have huge matrices which are sparse. This forces Oscar to allocate the full matrix as a dense matrix. If we had the possibility to pull the morphisms from the Singular side as a list of module elements (images of the generators), we would probably be better off.

[jankoboehm]

Even if not urgent since anyway seemingly fast, we should use a ModuleGens to transfer sparse data.

[HechtiDerLachs]

I'm taking back my doubts about the speed of Oscar vs. Singular in koszul_homology. Look at this:

julia> n = 5
5

julia> R, _ = polynomial_ring(QQ, "x"=>1:n)
(Multivariate polynomial ring in 5 variables over QQ, QQMPolyRingElem[x[1], x[2], x[3], x[4], x[5]])

julia> FR = FreeMod(R, 1)
Free module of rank 1 over Multivariate polynomial ring in 5 variables over QQ

julia> IR = ideal(R, gens(R))
ideal(x[1], x[2], x[3], x[4], x[5])

julia> IR = IR*IR
ideal(x[1]^2, x[1]*x[2], x[1]*x[3], x[1]*x[4], x[1]*x[5], x[2]^2, x[2]*x[3], x[2]*x[4], x[2]*x[5], x[3]^2, x[3]*x[4], x[3]*x[5], x[4]^2, x[4]*x[5], x[5]^2)

julia> @time koszul_homology(gens(IR), FR, 5);
 14.186331 seconds (32.92 M allocations: 1.295 GiB, 6.16% gc time)

julia> @time Oscar._koszul_homology_from_singular(gens(IR), FR, 5);
197.457997 seconds (332.20 M allocations: 10.411 GiB, 0.16% gc time)

[HechtiDerLachs]

With the depth, however, I do not get as good results as with singular. I temporarily rewinded the depth method to use the Oscar version and I got this:

julia> n = 4
4

julia> R, _ = polynomial_ring(QQ, "x"=>1:n)
(Multivariate polynomial ring in 4 variables over QQ, QQMPolyRingElem[x[1], x[2], x[3], x[4]])

julia> FR = FreeMod(R, 1)
Free module of rank 1 over Multivariate polynomial ring in 4 variables over QQ

julia> IR = ideal(R, gens(R))
ideal(x[1], x[2], x[3], x[4])

julia> IR = IR*IR
ideal(x[1]^2, x[1]*x[2], x[1]*x[3], x[1]*x[4], x[2]^2, x[2]*x[3], x[2]*x[4], x[3]^2, x[3]*x[4], x[4]^2)

julia> @time depth(IR, FR) # Oscar version
  0.409409 seconds (1.72 M allocations: 69.002 MiB)
4

julia> @time Oscar._depth_from_singular(IR, FR)
  0.188377 seconds (419 allocations: 21.906 KiB)
4

julia> @time depth(IR, FR)
  0.374969 seconds (1.72 M allocations: 68.984 MiB)
4

julia> @time Oscar._depth_from_singular(IR, FR)
  0.184218 seconds (417 allocations: 21.875 KiB)
4

[HechtiDerLachs]

I can not find the reason why Singular's depth computation is so much faster. The only thing I could imagine is that the modulo command used here is simply faster than the homology method we're using in Oscar. Other than that, I have no idea. @jankoboehm : Do you think this might be the reason?

[jankoboehm]

At the moment the homology computes kernel (which might be expensive), and image (which is free of cost) creates a subquo (which is free of cost), and then does membership to determine whether the generators are (hence the module is) zero, which triggers a Gröbner basis computation of the image. One could go via a presentation to determine whether the module is zero (which would rely on modulo). The other thing is, that we do a lot of bookkeeping in Oscar to have things free of an ordering of a basis, which might also cost a bit and is not needed to find that something is zero.

I think as a first step, we should restructure the is_zero for subquos a bit relying rather on division of one ModuleGens by another, than checking is_zero for individual generators, to avoid copying things around. I will have a look. If that does nor work, we can try to use the presentation command. If (and only if) that is slower than in Singular, we can map it to modulo.

[lgoettgens]

I think that this PR tries to achieve many different things, in particular, it does a lot more than written in the PR title. What do you think about splitting it into one that just implements exterior powers for different types of modules, and a second one with everything about homology happening here that builds on the first one?
This would make reviewing this PR easier, we can probably merge the first one a lot faster, and in the case of future bisections to find bugs there is a finer step size in changes.

[HechtiDerLachs]

I think that this PR tries to achieve many different things, in particular, it does a lot more than written in the PR title. What do you think about splitting it into one that just implements exterior powers for different types of modules, and a second one with everything about homology happening here that builds on the first one? This would make reviewing this PR easier, we can probably merge the first one a lot faster, and in the case of future bisections to find bugs there is a finer step size in changes.

No, I strongly object to this proposal.

We have discussed splitting this at earlier points and decided against it. Splitting it now will produce several hours of more work for the person performing the split with relatively little benefit. Reviewing will probably become a little more structured, but the actual amount of work will not decrease. Moreover, the review has already taken place incrementally here and I don't see the point in doing these discussions again on different PRs.

If anything, I can change the title again.

[HechtiDerLachs]

I spoke with @fingolfin and he has no objections to merging this anymore. If needed, we can clean up things later.