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
343 stars 126 forks source link

`grassmann_pluecker_ideal` saved GB wrong #4018

Closed YueRen closed 1 month ago

YueRen commented 2 months ago

I'm very happy that OSCAR has a constructor for Pluecker ideals, but the output is not very usable. For example, asking for its dim raises and error and asking for its degree returns something wrong:

julia> I = grassmann_pluecker_ideal(2,5)
Ideal generated by
  x[1]*x[8] - x[2]*x[6] + x[3]*x[5]
  x[1]*x[9] - x[2]*x[7] + x[4]*x[5]
  x[1]*x[10] - x[3]*x[7] + x[4]*x[6]
  x[2]*x[10] - x[3]*x[9] + x[4]*x[8]
  x[5]*x[10] - x[6]*x[9] + x[7]*x[8]

julia> dim(I)
ERROR: Not a Groebner basis
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] dimension(I::Singular.sideal{Singular.spoly{Singular.n_Q}})
   @ Singular ~/.julia/packages/Singular/A4riV/src/ideal/ideal.jl:105
 [3] __compute_dim__(I::MPolyIdeal{MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}})
   @ Oscar ~/.julia/dev/Oscar/src/Rings/mpoly-ideals.jl:1811
 [4] #379
   @ ~/.julia/packages/AbstractAlgebra/tABFQ/src/Attributes.jl:357 [inlined]
 [5] get!(default::Oscar.var"#379#380"{MPolyIdeal{MPolyDecRingElem{…}}}, h::Dict{Symbol, Any}, key::Symbol)
   @ Base ./dict.jl:479
 [6] get_attribute!
   @ ~/.julia/packages/AbstractAlgebra/tABFQ/src/Attributes.jl:230 [inlined]
 [7] dim(I::MPolyIdeal{MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}})
   @ Oscar ~/.julia/dev/Oscar/src/Rings/mpoly-ideals.jl:1806
 [8] top-level scope
   @ REPL[19]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> degree(I)
1
thofma commented 2 months ago

We started saving knowledge about a groebner basis in https://github.com/oscar-system/Oscar.jl/pull/3382. Maybe related to this?

YueRen commented 2 months ago

Good point, so either the fact that the generating set is a GB is wrong, or the GB does not match the ordering.

To me, the generators look pretty unassuming (the default "easy" generators that I would expect). I don't think they form a GB.

wdecker commented 2 months ago

Sorry, I do not understand the discussion, may be since I did not follow the history of this: We can easily check whether the stored generating set is a GB with respect to some ordering. If not, why not storing the correct GB?

YueRen commented 2 months ago

@wdecker The error is exactly because the stored generating set is no GB with the stored ordering. The construction of the generating set follows the usual formula for Pluecker relations, and there are two options:

  1. the resulting generating said is no GB with respect to any ordering (then https://github.com/oscar-system/Oscar.jl/pull/4019 should be merged)
  2. the ordering is wrong (then the right ordering should be provided)

I think (1) is the case, but only because if the standard construction would result in a GB in some ordering, I would have learned it by now.

wdecker commented 2 months ago

I did not have the time to check details here. But here are a guess or rather some questions to @antonydellavecchia:

  1. From the documentation, I guess that you make use of the straightening algorithm. True?
  2. If true, why do you use variables indexed by just one number as default: julia> I = grassmann_pluecker_ideal(2, 5) Ideal generated by x[1]x[8] - x[2]x[6] + x[3]x[5] x[1]x[9] - x[2]x[7] + x[4]x[5] x[1]x[10] - x[3]x[7] + x[4]x[6] x[2]x[10] - x[3]x[9] + x[4]x[8] x[5]x[10] - x[6]x[9] + x[7]x[8] This does not reflect the straightening algorithm. Here, you do better: julia> II = flag_pluecker_ideal([2], 5) Ideal generated by -x[[1, 4]]x[[2, 3]] + x[[2, 4]]x[[1, 3]] - x[[3, 4]]x[[1, 2]] -x[[1, 5]]x[[2, 3]] + x[[2, 5]]x[[1, 3]] - x[[3, 5]]x[[1, 2]] -x[[2, 5]]x[[3, 4]] + x[[3, 5]]x[[2, 4]] - x[[4, 5]]x[[2, 3]] -x[[1, 5]]x[[3, 4]] + x[[3, 5]]x[[1, 4]] - x[[4, 5]]x[[1, 3]] -x[[1, 5]]x[[2, 4]] + x[[2, 5]]x[[1, 4]] - x[[4, 5]]x[[1, 2]]
  3. The result of the straightening algorithm is indeed a degrevlex GB but only when you order the variables in a specific way. Quite likely, this is not the case in the polymake implementation. True?
benlorenz commented 2 months ago

This does look very much like a groebner basis:

julia> I = grassmann_pluecker_ideal(2,5);

julia> Icopy = ideal(gens(I));

julia> Icopy.gens.isGB
false

julia> is_groebner_basis(Icopy.gens)
true

julia> I.gens
Gröbner basis with elements
  1 -> x[3]*x[5] - x[2]*x[6] + x[1]*x[8]
  2 -> x[4]*x[5] - x[2]*x[7] + x[1]*x[9]
  3 -> x[4]*x[6] - x[3]*x[7] + x[1]*x[10]
  4 -> x[4]*x[8] - x[3]*x[9] + x[2]*x[10]
  5 -> x[7]*x[8] - x[6]*x[9] + x[5]*x[10]
with respect to the ordering
  degrevlex([x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10]])

julia> Icopy.gens
Gröbner basis with elements
  1 -> x[3]*x[5] - x[2]*x[6] + x[1]*x[8]
  2 -> x[4]*x[5] - x[2]*x[7] + x[1]*x[9]
  3 -> x[4]*x[6] - x[3]*x[7] + x[1]*x[10]
  4 -> x[4]*x[8] - x[3]*x[9] + x[2]*x[10]
  5 -> x[7]*x[8] - x[6]*x[9] + x[5]*x[10]
with respect to the ordering
  degrevlex([x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10]])

julia> dim(Icopy)
7

julia> degree(Icopy)
5

The same works for other dimensions as well.

The reference in the polymake code is:

Sturmfels: Algorithms in invariant theory, Springer, 2nd ed., 2008

There is also code that generates variable names but that function is not straightforward to use from Oscar, for the $2,5$ case it looks like this:

p_(0,1) p_(0,2) p_(0,3) p_(0,4) p_(1,2) p_(1,3) p_(1,4) p_(2,3) p_(2,4) p_(3,4)

But these are 0-based and not 1-based as one would expect in Oscar.

YueRen commented 2 months ago

This is weird, as dim claims that it isn't. Here is an example summing up the issue:

R,p = polynomial_ring(QQ,"p"=>sort(AbstractAlgebra.combinations(5,2)))
I = grassmann_pluecker_ideal(R,2,5)
dim(I) # error due to missing GB
issetequal(groebner_basis(ideal(gens(I))),gens(I)) # confirms that gens are GB
wdecker commented 2 months ago

p(0,1) p(0,2) p(0,3) p(0,4) p(1,2) p(1,3) p(1,4) p(2,3) p(2,4) p(3,4)

Yes, precisely, the indices must be lex ordered. How about using p12 etc in Oscar

YueRen commented 2 months ago

There is another problem with the output of grassamann_pluecker_ideal, it doesn't live in the polynomial ring it is supposed to:

julia> R,x = polynomial_ring(QQ,10);

julia> I = grassmann_pluecker_ideal(R,2,5); # should live in R

julia> phi = hom(R,R,gens(R)); # identitiy on R

julia> phi(I) # Error
ERROR: ideal not defined over the domain of the map
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] (::Oscar.MPolyAnyMap{QQMPolyRing, QQMPolyRing, Nothing, QQMPolyRingElem})(I::MPolyIdeal{MPolyDecRingElem{…}})
   @ Oscar ~/.julia/dev/Oscar/experimental/Schemes/src/IdealSheaves.jl:1094
 [3] top-level scope
   @ REPL[24]:1
Some type information was truncated. Use `show(err)` to see complete types.
benlorenz commented 2 months ago

There is another problem with the output of grassamann_pluecker_ideal, it doesn't live in the polynomial ring it is supposed to:

help?> grassmann_pluecker_ideal
...
 If the input ring is not graded, return the ideal in the ring with the standard grading.

See #3765 and #3632

wdecker commented 2 months ago

@benlorenz Thx for the explanation. We will check and correct the problem.

joschmitt commented 2 months ago

I tried this with Singular.jl master, so including oscar-system/Singular.jl#820, and this is not completely fixed.

julia> I = grassmann_pluecker_ideal(2,5)
Ideal generated by
  x[1]*x[8] - x[2]*x[6] + x[3]*x[5]
  x[1]*x[9] - x[2]*x[7] + x[4]*x[5]
  x[1]*x[10] - x[3]*x[7] + x[4]*x[6]
  x[2]*x[10] - x[3]*x[9] + x[4]*x[8]
  x[5]*x[10] - x[6]*x[9] + x[7]*x[8]

julia> dim(I)
7

julia> degree(I)
1

The dimension is fine, but the degree should be:

julia> degree(ideal(gens(I)))
5

If I compute the degree without computing the dimension first, I actually get an error:

julia> I = grassmann_pluecker_ideal(2,5)
Ideal generated by
  x[1]*x[8] - x[2]*x[6] + x[3]*x[5]
  x[1]*x[9] - x[2]*x[7] + x[4]*x[5]
  x[1]*x[10] - x[3]*x[7] + x[4]*x[6]
  x[2]*x[10] - x[3]*x[9] + x[4]*x[8]
  x[5]*x[10] - x[6]*x[9] + x[7]*x[8]

julia> degree(I)
ERROR: UndefRefError: access to undefined reference
Stacktrace:
  [1] getproperty
    @ ./Base.jl:37 [inlined]
  [2] getproperty
    @ ~/.julia/dev/Oscar/src/Rings/mpoly.jl:220 [inlined]
  [3] HilbertData(I::MPolyIdeal{MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}})
    @ Oscar ~/.julia/dev/Oscar/src/Rings/mpoly-graded.jl:1506
  [4] degree(A::MPolyQuoRing{MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}})
    @ Oscar ~/.julia/dev/Oscar/src/Rings/mpoly-affine-algebras.jl:461
  [5] __compute_degree__(I::MPolyIdeal{MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}})
    @ Oscar ~/.julia/dev/Oscar/src/Rings/mpoly-graded.jl:2439
  [6] #246
    @ ~/.julia/packages/AbstractAlgebra/bB5QU/src/Attributes.jl:357 [inlined]
  [7] get!(default::Oscar.var"#246#247"{MPolyIdeal{MPolyDecRingElem{…}}}, h::Dict{Symbol, Any}, key::Symbol)
    @ Base ./dict.jl:479
  [8] get_attribute!
    @ ~/.julia/packages/AbstractAlgebra/bB5QU/src/Attributes.jl:230 [inlined]
  [9] degree(I::MPolyIdeal{MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}})
    @ Oscar ~/.julia/dev/Oscar/src/Rings/mpoly-graded.jl:2432
 [10] top-level scope
    @ REPL[33]:1
Some type information was truncated. Use `show(err)` to see complete types.

This also happens without oscar-system/Singular.jl#820.

jankoboehm commented 2 months ago

We should change groebner_assure to the following (adding the missing else)

function groebner_assure(I::MPolyIdeal, complete_reduction::Bool = false, need_global::Bool = false)
  if !isempty(I.gb)
    for G in values(I.gb)
      need_global || return G
      is_global(G.ord) || continue
      complete_reduction || return G
      if !G.isReduced
        I.gb[G.ord] = _compute_standard_basis(G, G.ord, true)
        return I.gb[G.ord]
      else
        return G
      end
    end
  end
  ord = default_ordering(base_ring(I))
  (need_global <= is_global(ord)) || error("Monomial ordering must be global.")
  I.gb[ord] = groebner_assure(I, ord, complete_reduction)
  return I.gb[ord]
end
joschmitt commented 2 months ago

I'm sorry, but I don't understand. With this change

diff --git a/src/Rings/groebner.jl b/src/Rings/groebner.jl
index 8496956053..58007d35d5 100644
--- a/src/Rings/groebner.jl
+++ b/src/Rings/groebner.jl
@@ -41,6 +41,8 @@ function groebner_assure(I::MPolyIdeal, complete_reduction::Bool = false, need_g
       if !G.isReduced
         I.gb[G.ord] = _compute_standard_basis(G, G.ord, true)
         return I.gb[G.ord]
+      else
+        return G
       end
     end
   end

I still get the very same error(s):

julia> I = grassmann_pluecker_ideal(2, 5)
Ideal generated by
  x[1]*x[8] - x[2]*x[6] + x[3]*x[5]
  x[1]*x[9] - x[2]*x[7] + x[4]*x[5]
  x[1]*x[10] - x[3]*x[7] + x[4]*x[6]
  x[2]*x[10] - x[3]*x[9] + x[4]*x[8]
  x[5]*x[10] - x[6]*x[9] + x[7]*x[8]

julia> degree(I)
ERROR: UndefRefError: access to undefined reference
[...]

julia> dim(I)
ERROR: Not a Groebner basis
[...]

julia> degree(I)
1