Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Exterior powers of finitely presented modules #2879

Merged
merged 36 commits into from
Oct 18, 2023

Conversation

HechtiDerLachs
Copy link
Collaborator

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

@fingolfin
Copy link
Member

@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 / ...

Comment on lines 134 to 144
@testset "multiplication map" begin
R, (x, y, u, v, w) = QQ[:x, :y, :u, :v, :w]

F = FreeMod(R, 5)
F3, mm = Oscar.exterior_power(F, 3)
v = (F[1], F[3], F[4])
u = (F[1], F[4], F[3])
@test mm(v) == -mm(u)
w = mm(v)
@test preimage(mm, w) == v
end
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@fieker : Is this like you had intended?

Comment on lines 72 to 83
@testset "all_monomials for graded free modules" begin
R, (x, y, u, v, w) = QQ[:x, :y, :u, :v, :w]
S, (x, y, u, v, w) = grade(R)

F = graded_free_module(S, [-1, 2])

for d in -4:4
amm = Oscar.all_monomials(F, d)
@test d < -1 || !isempty(amm)
@test all(x->degree(x) == grading_group(F)([d]), amm)
end
end
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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
Copy link
Collaborator Author

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
Copy link
Contributor

jankoboehm commented Oct 6, 2023

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
Copy link
Member

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
Copy link
Collaborator Author

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
Copy link
Member

@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.

return first(koszul_duals([v], cached=cached))
end

function induced_map_on_exterior_power(phi::FreeModuleHom{<:FreeMod, <:FreeMod, Nothing}, p::Int;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
Copy link
Collaborator Author

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.

Copy link
Member

@fingolfin fingolfin left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

# generators of M. But they can also be used for other purposes
# like enumerating subsets of p elements out of a set with n elements.
########################################################################
mutable struct OrderedMultiIndex{IntType<:IntegerUnion}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is this struct mutable?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

src/Combinatorics/OrderedMultiIndex.jl Outdated Show resolved Hide resolved
@assert bound(a) == bound(b) "multiindices must have the same bounds"

# in case of a double index return zero
any(x->(x in indices(b)), indices(a)) && return 0, a
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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[]

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

src/Combinatorics/OrderedMultiIndex.jl Outdated Show resolved Hide resolved

########################################################################
# A data type to facilitate iteration over all ordered multiindices
# of the form 0 < i₁ < i₂ < … < iₚ ≤ n for fixed 0 ≤ p ≤ n.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
Copy link
Collaborator Author

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
Copy link
Collaborator

wdecker commented Oct 12, 2023

@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.

R = parent(V[1])
@assert all(x->parent(x) == R, V)
n = length(V)
KM = [koszul_matrix(V, i) for i=n:-1:1]
KM = [_koszul_matrix_from_singular(V, i) for i=n:-1:1]
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

@HechtiDerLachs
Copy link
Collaborator Author

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
Copy link
Collaborator Author

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
Copy link
Collaborator Author

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
Copy link
Contributor

jankoboehm commented Oct 13, 2023

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
Copy link
Member

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
Copy link
Collaborator Author

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
Copy link
Collaborator Author

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

@wdecker wdecker merged commit f9d8294 into oscar-system:master Oct 18, 2023
20 of 28 checks passed
@HechtiDerLachs HechtiDerLachs deleted the exterior_powers branch October 19, 2023 09:57
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

5 participants