oscar-system / Oscar.jl

A comprehensive open source computer algebra system for computations in algebra, geometry, and number theory.
https://www.oscar-system.org
Other
344 stars 126 forks source link

questions concerning default rings #932

Closed micjoswig closed 2 years ago

micjoswig commented 2 years ago

https://github.com/oscar-system/Oscar.jl/blob/5fa3a5c7e538f11eca4972202ff88572cc18d128/src/ToricVarieties/NormalToricVarieties/attributes.jl#L73-L75

to be compared with:

https://github.com/oscar-system/Oscar.jl/blob/5fa3a5c7e538f11eca4972202ff88572cc18d128/src/Combinatorics/SimplicialComplexes.jl#L223-L227

I have several recommendations / questions concerning this code:

  1. I guess we should always allow for an optional ring argument.
  2. Do we want to take QQ or ZZ as a default coefficient ring here? I think ZZ is more natural (and more general, as it can be tensored into whatever).
  3. How do we want to pick the default variable names? Implicitly, as in the simplicial complex code, or explicitly by specifying a vector of variables?
micjoswig commented 2 years ago

We now have two implementations of Stanley-Reisner ideals and rings. Here is a comparsion of the running times (while this is only one example, the behavior is typical):

julia> P = polarize(Polyhedron(Polymake.polytope.rand_sphere(5,60)))
A polyhedron in ambient dimension 5

julia> V = NormalToricVariety(P)
A normal toric variety

julia> @time stanley_reisner_ideal(V);
342.250888 seconds (7.97 M allocations: 689.097 MiB, 0.03% gc time, 0.00% compilation time)

julia> function allrows(I::IncidenceMatrix)
           m,n = size(I)
           return collect( Vector{Int}(Polymake.row(I,k)) for k = 1:m )
       end
allrows (generic function with 1 method)

julia> I = P.pm_polytope.FACETS_THRU_VERTICES;

julia> K = SimplicialComplex(allrows(I))
Abstract simplicial complex of dimension 4 on 60 vertices

julia> @time stanley_reisner_ideal(K);
  0.024923 seconds (49.94 k allocations: 3.336 MiB)

Does the toric implementation go beyond simplicial fans? If not, you probably want to refactor through simplicial complexes.

The two implementations differ as far as the base rings are concerned; see above. Yet this should not matter much for the running times.

HereAround commented 2 years ago
  1. I guess we should always allow for an optional ring argument.

For clarity - is this optional argument supposed to be QQ or ZZ? Or rather the Cox ring?

If the Cox ring is meant, then providing its grading needs to be considered: We could either accept a non-graded ring and compute a grading (a.k.a. one of the many cokernel projection of the map from the lattice into the torus-invariant divisor group) or - and that happens to be fairly convenient/necessary for my physics applications - provide the graded ring.

  1. Do we want to take QQ or ZZ as a default coefficient ring here? I think ZZ is more natural (and more general, as it can be tensored into whatever).

I am happy with both.

For completeness, let me mention that I asked a similar question in Slack a while back (about October or November). At that time, I was trying to gauge if people would prefer a (formal - if it exists) ring of complex numbers or QQ. At that point, I received suggestions for QQ, which is therefore the current choice for toric varieties.

  1. How do we want to pick the default variable names? Implicitly, as in the simplicial complex code, or explicitly by specifying a vector of variables?

I would think that all of the following cases need to be considered:

  1. Default variable names (e.g. x1, x2, x3,...)
  2. All variable names specified by the user (e.g. u1, u2, w1, w2, w3, ...)
  3. For some special constructions, a more sophisticated choice is needed.

The first special case that comes to my mind are del Pezzo surfaces. It seem standard to label the blowup/exceptional coordinate with e (at least based on my experience/community). Of course, this name could be changed with yet another optional argument - the current choice in toric varieties however is e. Indeed, as part of the caching for toric varieties, I have implemented this so that we now have the following:

dP3 = del_pezzo(3) A normal, non-affine, smooth, projective, gorenstein, q-gorenstein, fano, 2-dimensional toric variety without torusfactor

cox_ring(dP3) Multivariate Polynomial Ring in 6 variables x[1], e[3], x[2], e[1], ..., e[2] over Rational Field graded by x[1] -> [1 0 0 0] e[3] -> [0 1 0 0] x[2] -> [0 0 1 0] e[1] -> [0 0 0 1] x[3] -> [1 1 0 -1] e[2] -> [-1 0 1 1]

Note that the current variable names in toric varieties are always x[ 1 ]. This is somewhat cumbersome and I was about to prepare a PR to change this to x1 and alike.

A similar situation arises whenever we compute the blowup of a toric variety. Then again, I would find it most convenient to label the new homogeneous variable/exceptional coordinate with e (or e2 if e is already taken...).

Yet another case happens if we compute the Cartesian product of two toric varieties A, B. If the homogeneous coordinates of A are a1, a2, ... and those of B are b1, b2, ... then it would seem natural to have homogeneous coordinate a1, a2,..., b1, b2, ... for A*B (and one can then easily infer the Stanley-Reisner ideal of A*B from those of A and B).

I am not sure how to choose the variable names for A*B if say x1 is both a variable of A and B. Likewise, a choice/convention has to be made if only A has variable names but B has none. (If both have none, then we could fall back to xi).

Of course, there can (and most likely are) other special cases/constructions which require special care.

Does the toric implementation go beyond simplicial fans? If not, you probably want to refactor through simplicial complexes.

@micjoswig It does not go beyond simpicial. So I agree that this should be refactored. In fact - great that your implementation is so much faster! Are there plans for similar implementations of the irrelevant ideal?

fieker commented 2 years ago

On Sun, Jan 02, 2022 at 04:42:50AM -0800, Martin Bies wrote:

  1. I guess we should always allow for an optional ring argument.

For clarity - is this optional argument supposed to be QQ or ZZ? Or rather the Cox ring?

If the Cox ring is meant, then providing its grading needs to be considered: We could either accept a non-graded ring and compute a grading (a.k.a. one of the many cokernel projection of the map from the lattice into the torus-invariant divisor group) or - and that happens to be fairly convenient/necessary for my physics applications - provide the graded ring.

  1. Do we want to take QQ or ZZ as a default coefficient ring here? I think ZZ is more natural (and more general, as it can be tensored into whatever).

I am happy with both. What do you want to do with the ring? QQ[x] is much larger and mostly mathemtically and algorithmically easier, e.g. Groebner in ZZ is a complete pain.

For completeness, let me mention that I asked a similar question in Slack a while back (about October or November). At that time, I was trying to gauge if people would prefer a (formal - if it exists) ring of complex numbers or QQ. At that point, I received suggestions for QQ, which is therefore the current choice for toric varieties.

  1. How do we want to pick the default variable names? Implicitly, as in the simplicial complex code, or explicitly by specifying a vector of variables?

I would think that all of the following cases need to be considered:

  1. Default variable names (e.g. x1, x2, x3,...)
  2. All variable names specified by the user (e.g. u1, u2, w1, w2, w3, ...)
  3. For some special constructions, a more sophisticated choice is needed. I'd add a generic function set_names!(Rx, ...) to allow the user to retroactively change the printing names - on non-cached rings.

This is only relevant if you're expected to print lots of polynomials. For interactive sessions it should be possible to extend our "hack" from number fields: if the ring is assigned to a variable, say R, then we could print the generators, by default as R[1], R[2], ... This would allow for copy&paste.

The first special case that comes to my mind are del Pezzo surfaces. It seem standard to label the blowup/exceptional coordinate with e (at least based on my experience/community). Of course, this name could be changed with yet another optional argument - the current choice in toric varieties however is e. Indeed, as part of the caching for toric varieties, I have implemented this so that we now have the following:

dP3 = del_pezzo(3) A normal, non-affine, smooth, projective, gorenstein, q-gorenstein, fano, 2-dimensional toric variety without torusfactor

cox_ring(dP3) Multivariate Polynomial Ring in 6 variables x[1], e[3], x[2], e[1], ..., e[2] over Rational Field graded by x[1] -> [1 0 0 0] e[3] -> [0 1 0 0] x[2] -> [0 0 1 0] e[1] -> [0 0 0 1] x[3] -> [1 1 0 -1] e[2] -> [-1 0 1 1]

Note that the current variable names in toric varieties are always x[ 1 ]. This is somewhat cumbersome and I was about to prepare a PR to change this to x1 and alike.

Do whatever you like... however, as a user: alternating x and e means I have to explicitly do R, x1, e3, x2, e2, x3, e1 = ... instead of R, x, e = .... the latter also would allow PolynomialRing(QQ, :x=>1:3, :e=>3:-1:1) if you want.

In general x1 looks better than x[1] x[1] can be used in copy&paste easily, as PolynomialRing(QQ, :x=>1:10) will return a vector x, so printing x[1] will translate into the correct object.

x1, x2, ... can not be generated programmatically, so when dealing with an unknown number of variables, anything but a vector (or array) is infeasable.

A similar situation arises whenever we compute the blowup of a toric variety. Then again, I would find it most convenient to label the new homogeneous variable/exceptional coordinate with e (or e2 if e is already taken...).

Yet another case happens if we compute the Cartesian product of two toric varieties A, B. If the homogeneous coordinates of A are a1, a2, ... and those of B are b1, b2, ... then it would seem natural to have homogeneous coordinate a1, a2,..., b1, b2, ... for A*B (and one can then easily infer the Stanley-Reisner ideal of A*B from those of A and B).

I am not sure how to choose the variable names for A*B if say x1 is both a variable of A and B. Likewise, a choice/convention has to be made if only A has variable names but B has none. (If both have none, then we could fall back to xi). well, here you have a different problem: unless A and B are created with cached = false, they are likely to be the identical ring, so possibly AB = AA = A here?

Of course, there can (and most likely are) other special cases/constructions which require special care.

There is absolutely no chance at all to make sure that names can be chosen well and automatically. You have a couple of options

To me, the printing names are mostly irrelevant as most elements of most rings will never get printed anyway...

Does the toric implementation go beyond simplicial fans? If not, you probably want to refactor through simplicial complexes.

@micjoswig It does not go beyond simpicial. So I agree that this should be refactored. In fact - great that your implementation is so much faster! Are there plans for similar implementations of the irrelevant ideal?

-- Reply to this email directly or view it on GitHub: https://github.com/oscar-system/Oscar.jl/issues/932#issuecomment-1003709991 You are receiving this because you were assigned.

Message ID: @.***>

fieker commented 2 years ago

On Thu, Dec 30, 2021 at 05:22:37AM -0800, Michael Joswig wrote:

We now have two implementations of Stanley-Resiner ideals and rings. Here is a comparsion of the running times (while this is only one example, the behavior is typical):


julia> P = polarize(Polyhedron(Polymake.polytope.rand_sphere(5,60)))
A polyhedron in ambient dimension 5

julia> V = NormalToricVariety(P)
A normal toric variety

julia> @time stanley_reisner_ideal(V);
342.250888 seconds (7.97 M allocations: 689.097 MiB, 0.03% gc time, 0.00% compilation time)
I got:
julia> @time I1 = stanley_reisner_ideal(V);           
ERROR: InexactError: Int64(70323215040313106432)
Stacktrace:
[1] Type
@ ./gmp.jl:363 [inlined]
[2] Int64(int::Polymake.IntegerAllocated)
@ Polymake ~/.julia/dev/Polymake/src/integers.jl:4
...
it appears that s.o. is trying to create an Int64 matrix when the
entries can be (are) huge...

function map_from_character_to_principal_divisors(v::AbstractNormalToricVariety) mat = transpose(Matrix{Int}(Polymake.common.primitive(pm_object(v).RAYS))) matrix = AbstractAlgebra.matrix(ZZ, mat) return hom(character_lattice(v), torusinvariant_divisor_group(v), matrix) end

is the culprit. Is there a mathematical reason why this should be small ints?

julia> function allrows(I::IncidenceMatrix) m,n = size(I) return collect( Vector{Int}(Polymake.row(I,k)) for k = 1:m ) end allrows (generic function with 1 method)

julia> I = P.pm_polytope.FACETS_THRU_VERTICES;

julia> K = SimplicialComplex(allrows(I)) Abstract simplicial complex of dimension 4 on 60 vertices

julia> @time stanley_reisner_ideal(K); 0.024923 seconds (49.94 k allocations: 3.336 MiB)


Does the toric implementation go beyond simplicial fans?  If not, you probably want to refactor through simplicial complexes. 

The two implementations differ as far as the base rings are concerned; see above. Yet this should not matter much for the running times.

-- 
Reply to this email directly or view it on GitHub:
https://github.com/oscar-system/Oscar.jl/issues/932#issuecomment-1003027093
You are receiving this because you were assigned.

Message ID: ***@***.***>
fieker commented 2 years ago

On Thu, Dec 30, 2021 at 05:22:37AM -0800, Michael Joswig wrote:

We now have two implementations of Stanley-Resiner ideals and rings. Here is a comparsion of the running times (while this is only one example, the behavior is typical): But, while I have no idea what they are supposed to be doing, they are computing different objects.

julia> P = polarize(Polyhedron(Polymake.polytope.rand_sphere(5,60)))
A polyhedron in ambient dimension 5

julia> V = NormalToricVariety(P)
A normal toric variety

julia> @time stanley_reisner_ideal(V);
This is an ideal in a Graded Ring, the grading has insanely large
entries, the coeff. ring is QQ
342.250888 seconds (7.97 M allocations: 689.097 MiB, 0.03% gc time, 0.00% compilation time)

julia> function allrows(I::IncidenceMatrix)
           m,n = size(I)
           return collect( Vector{Int}(Polymake.row(I,k)) for k = 1:m )
       end
allrows (generic function with 1 method)

julia> I = P.pm_polytope.FACETS_THRU_VERTICES;

julia> K = SimplicialComplex(allrows(I))
Abstract simplicial complex of dimension 4 on 60 vertices

julia> @time stanley_reisner_ideal(K);
This is an ungraded ideal over ZZ
Is the time in the grading?

  0.024923 seconds (49.94 k allocations: 3.336 MiB)

Does the toric implementation go beyond simplicial fans? If not, you probably want to refactor through simplicial complexes.

The two implementations differ as far as the base rings are concerned; see above. Yet this should not matter much for the running times.

-- Reply to this email directly or view it on GitHub: https://github.com/oscar-system/Oscar.jl/issues/932#issuecomment-1003027093 You are receiving this because you were assigned.

Message ID: @.***>

fieker commented 2 years ago

Looking at the implementation.... my guess is that the "new" implementation is "just" using Polymake in a suboptimal way. The (a) bottleneck is "primitive_collections" of the "fan". The fast version is using "minimal_nonfaces"

micjoswig commented 2 years ago

The fast version is using "minimal_nonfaces"

Yip, that implementation was tuned quite a bit in the past.

micjoswig commented 2 years ago
  1. I guess we should always allow for an optional ring argument.

For clarity - is this optional argument supposed to be QQ or ZZ? Or rather the Cox ring?

Yes, I meant the Cox ring (which is called "total coordinate ring" in CLS, by the way). Please add that name to the doc, too, with a ref.

If the Cox ring is meant, then providing its grading needs to be considered: We could either accept a non-graded ring and compute a grading (a.k.a. one of the many cokernel projection of the map from the lattice into the torus-invariant divisor group) or - and that happens to be fairly convenient/necessary for my physics applications - provide the graded ring.

Understood. I agree that this should be discussed.

  1. Do we want to take QQ or ZZ as a default coefficient ring here? I think ZZ is more natural (and more general, as it can be tensored into whatever).

I am happy with both.

For completeness, let me mention that I asked a similar question in Slack a while back (about October or November). At that time, I was trying to gauge if people would prefer a (formal - if it exists) ring of complex numbers or QQ. At that point, I received suggestions for QQ, which is therefore the current choice for toric varieties.

This raises an interesting design question. I guess, if you are interested in the cohomology of toric varieties, then QQ is one natural choice. Yet, the ideal always has integer coefficients (for arbitrary simplicial complexes). Typical questions which occur there involve, e.g., Cohen-Macaulayness over some field, and here QQ makes as much sense as any finite field of prime order.

My current implementation for simplicial complexes insists on integer coefficients, but I would prefer that the user can also choose the coefficient ring (together with the variables).

  1. How do we want to pick the default variable names? Implicitly, as in the simplicial complex code, or explicitly by specifying a vector of variables?

I would think that all of the following cases need to be considered:

1. Default variable names (e.g. `x1, x2, x3,...`)

Yip.

2. All variable names specified by the user (e.g. `u1, u2, w1, w2, w3, ...`)

Yip.

3. For some special constructions, a more sophisticated choice is needed.

The first special case that comes to my mind are del Pezzo surfaces. It seem standard to label the blowup/exceptional coordinate with e (at least based on my experience/community). Of course, this name could be changed with yet another optional argument - the current choice in toric varieties however is e. Indeed, as part of the caching for toric varieties, I have implemented this so that we now have the following:

dP3 = del_pezzo(3) A normal, non-affine, smooth, projective, gorenstein, q-gorenstein, fano, 2-dimensional toric variety without torusfactor

cox_ring(dP3) Multivariate Polynomial Ring in 6 variables x[1], e[3], x[2], e[1], ..., e[2] over Rational Field graded by x[1] -> [1 0 0 0] e[3] -> [0 1 0 0] x[2] -> [0 0 1 0] e[1] -> [0 0 0 1] x[3] -> [1 1 0 -1] e[2] -> [-1 0 1 1]

Note that the current variable names in toric varieties are always x[ 1 ]. This is somewhat cumbersome and I was about to prepare a PR to change this to x1 and alike.

A similar situation arises whenever we compute the blowup of a toric variety. Then again, I would find it most convenient to label the new homogeneous variable/exceptional coordinate with e (or e2 if e is already taken...).

Sure this makes sense. But isn't this a special case of the previous? Where the user specifies some way of naming the variables?

Yet another case happens if we compute the Cartesian product of two toric varieties A, B. If the homogeneous coordinates of A are a1, a2, ... and those of B are b1, b2, ... then it would seem natural to have homogeneous coordinate a1, a2,..., b1, b2, ... for A*B (and one can then easily infer the Stanley-Reisner ideal of A*B from those of A and B).

I am not sure how to choose the variable names for A*B if say x1 is both a variable of A and B. Likewise, a choice/convention has to be made if only A has variable names but B has none. (If both have none, then we could fall back to xi).

Of course, there can (and most likely are) other special cases/constructions which require special care.

Does the toric implementation go beyond simplicial fans? If not, you probably want to refactor through simplicial complexes.

@micjoswig It does not go beyond simpicial. So I agree that this should be refactored. In fact - great that your implementation is so much faster! Are there plans for similar implementations of the irrelevant ideal?

We do not have this yet. But this is much simpler to compute than the Stanley-Reisner ideal. It would make sense to have some technology in OSCAR for dealing with monomials ideals in general.

micjoswig commented 2 years ago

What do you want to do with the ring? QQ[x] is much larger and mostly mathemtically and algorithmically easier, e.g. Groebner in ZZ is a complete pain.

As I sketched above: integer coefficients are the natural choice (for simplicial complexes; one could argue QQ is better for toric varieties), but generally the user should have an option to choose other coefficients.

micjoswig commented 2 years ago

I feel this is a bit too complicated to discuss here. It would be good to have a short session for coding up some sketches together. After the next OSCAR Friday meeting? Some quick sketches, aided by experts on rings in OSCAR, should be enough. Then everybody can fill in some part of the blanks.

HereAround commented 2 years ago

I agree that this might be too involved to discuss in large detail here. A first draft towards resolving the points raised (or at least some of them) can be found here: https://github.com/oscar-system/Oscar.jl/pull/939

In the current form, this PR is bad/not satisfying (more details in the PR) and should NOT be merged. It should however help the discussion. Hopefully, this can eventually be "the" solution.