TropicalGeometry: new linear_space algorithm for triv val on QQ (osca…

…r-system#3759)

docs/oscar_references.bib

@@ -1808,6 +1808,20 @@ @Article{RSS03
   year          = {2003}
 }
 
+@Article{Rin13,
+  author        = {Rinc{\'o}n, Felipe},
+  title         = {Computing tropical linear spaces},
+  zbl           = {1319.14060},
+  journal       = {J. Symb. Comput.},
+  fjournal      = {Journal of Symbolic Computation},
+  volume        = {51},
+  pages         = {86--98},
+  year          = {2013},
+  doi           = {10.1016/j.jsc.2012.03.008},
+  language      = {English},
+  zbmath        = {6143064}
+}
+
 @Book{SS03,
   author        = {Semple, Charles and Steel, Mike},
   title         = {Phylogenetics},

docs/src/TropicalGeometry/linear_space.md

@@ -10,10 +10,12 @@ Objects of type `TropicalLinearSpace` need to be embedded, abstract tropical lin
 
 ## Constructors
 In addition to converting from `TropicalVariety`, objects of type `TropicalLinearSpace` can be constructed from:
-- Pluecker vectors over a tropical semiring,
-- Pluecker vectors over a field and a tropical semiring map,
-- matrices over a tropical semiring,
-- matrices over a field and a tropical semiring map.
+1. Pluecker vectors over a tropical semiring: uses a low-level implementation in `polymake`
+2. Pluecker vectors over a field and a tropical semiring map: computes coordinatewise valuation and uses constructor (1.)
+3. matrices over a tropical semiring: computes tropical minors and uses constructor (1.)
+4. matrices over a field and a tropical semiring map.
+    - if matrix over `QQ` and tropical semiring map is trivial, uses an implementation of Rincon's algorithm [Rin13](@cite) in `polymake`
+    - for general input, computes minors and uses constructor (2.)
 ```@docs
 tropical_linear_space
 ```

src/TropicalGeometry/linear_space.jl

@@ -11,6 +11,7 @@
     multiplicities::Vector{ZZRingElem}
 
     # tropical linear spaces need to be embedded
+    # Note: multiplicities needs to be passed explicitly to avoid complications with a third optional parameter (tropical semiring map)
     function TropicalLinearSpace{minOrMax,true}(Sigma::PolyhedralComplex, multiplicities::Vector{ZZRingElem}) where {minOrMax<:Union{typeof(min),typeof(max)}}
         return new{minOrMax,true}(Sigma,multiplicities)
     end
@@ -65,7 +66,7 @@ function add_missing_lineality_from_polymake(Sigma::PolyhedralComplex)
     VRmat = matrix(QQ,Vector{QQFieldElem}.(collect(VR)))
     L = vcat(matrix(QQ,lineality_space(Sigma)),
              matrix(QQ,ones(Int,1,ambient_dim(Sigma))))
-    return polyhedral_complex(IM,VRmat,rayIndices,L)
+    return polyhedral_complex(IM,VRmat,rayIndices,L; non_redundant=true)
 end
 
 
@@ -208,10 +209,25 @@ function tropical_linear_space(k::Int, n::Int, plueckerVector::Vector, nu::Union
 end
 
 
+# helper function that takes a matrix (tropical or algebraic)
+# and returns the non-zero minors and their indices
+function compute_pluecker_indices_and_vector(A::MatElem)
+    n = ncols(A)
+    k = nrows(A)
+    @req n>=k "matrix for Pluecker vector cannot have more rows than columns"
+    plueckerIndices = AbstractAlgebra.combinations(1:n,k)
+    plueckerVector = AbstractAlgebra.minors(A,k)
+    nonZeroIndices = findall(!iszero,plueckerVector)
+    plueckerIndices = plueckerIndices[nonZeroIndices]
+    plueckerVector = plueckerVector[nonZeroIndices]
+    return plueckerIndices, plueckerVector
+end
+
+
 @doc raw"""
     tropical_linear_space(A::MatElem{<:TropicalSemiringElem}; weighted_polyhedral_complex_only::Bool=false)
 
-Return a tropical linear space whose Pluecker vector are the tropical minors of `A`.  If `weighted_polyhedral_complex==true`, will not cache any extra information.
+Return a tropical linear space whose Pluecker vector are the tropical minors of `A`.  Assumes that `ncols(A)>=nrows(A)`.  If `weighted_polyhedral_complex==true`, will not cache any extra information.
 
 # Examples
 ```jldoctest
@@ -224,22 +240,21 @@ Min tropical linear space
 
 ```
 """
-function tropical_linear_space(A::Union{Matrix{T},MatElem{T}}; weighted_polyhedral_complex_only::Bool=false) where {T<:TropicalSemiringElem}
-    # convert to Oscar matrix (necessary for AbstractAlgebra.minors)
-    A isa Matrix && (A = matrix(parent(first(A)),A))
+function tropical_linear_space(A::MatElem{T}; weighted_polyhedral_complex_only::Bool=false) where {T<:TropicalSemiringElem}
+    @req ncols(A)>=nrows(A) "matrix for Pluecker vector cannot have more rows than columns"
 
-    n = max(nrows(A), ncols(A))
-    k = min(nrows(A), ncols(A))
-    plueckerIndices = AbstractAlgebra.combinations(1:n,k)
-    plueckerVector = AbstractAlgebra.minors(A, k)
+    plueckerIndices,plueckerVector = compute_pluecker_indices_and_vector(A)
     TropL = tropical_linear_space(plueckerIndices, plueckerVector,
                                   weighted_polyhedral_complex_only=weighted_polyhedral_complex_only)
     if !weighted_polyhedral_complex_only
         set_attribute!(TropL,:tropical_matrix,A)
     end
     return TropL
 end
-
+function tropical_linear_space(A::Matrix{T}; weighted_polyhedral_complex_only::Bool=false) where {T<:TropicalSemiringElem}
+    # convert to Oscar matrix (necessary for AbstractAlgebra.minors)
+    return tropical_linear_space(matrix(parent(first(A)),A), weighted_polyhedral_complex_only=weighted_polyhedral_complex_only)
+end
 
 @doc raw"""
     tropical_linear_space(A::MatElem,nu::TropicalSemiringMap; weighted_polyhedral_complex_only::Bool=false)
@@ -261,19 +276,13 @@ Min tropical linear space
 ```
 """
 function tropical_linear_space(A::MatElem, nu::Union{Nothing,TropicalSemiringMap}=nothing; weighted_polyhedral_complex_only::Bool=false)
-    # if nu unspecified, initialize as the trivial valuation + min convention
-    isnothing(nu) && (nu=tropical_semiring_map(base_ring(A)))
-
     # compute reduced row echelon form of A
     # and remove all zero rows so that matrix is of full rank
     _,A = rref(A)
     nonzeroRowIndices = findall(!iszero,[A[i,:] for i in 1:nrows(A)])
     A = A[nonzeroRowIndices,:]
 
-    n = max(nrows(A), ncols(A))
-    k = min(nrows(A), ncols(A))
-    plueckerIndices = AbstractAlgebra.combinations(1:n,k)
-    plueckerVector = AbstractAlgebra.minors(A, k)
+    plueckerIndices,plueckerVector = compute_pluecker_indices_and_vector(A)
     TropL = tropical_linear_space(plueckerIndices, plueckerVector, nu,
                                   weighted_polyhedral_complex_only=weighted_polyhedral_complex_only)
     if !weighted_polyhedral_complex_only
@@ -282,15 +291,52 @@ function tropical_linear_space(A::MatElem, nu::Union{Nothing,TropicalSemiringMap
     end
     return TropL
 end
+function tropical_linear_space(A::Matrix, nu::Union{Nothing,TropicalSemiringMap}=nothing; weighted_polyhedral_complex_only::Bool=false)
+    return tropical_linear_space(matrix(parent(first(A)),A), nu, weighted_polyhedral_complex_only=weighted_polyhedral_complex_only)
+end
+
+# Special implementation for trivial valuation over the rationals using Polymake.tropical.matroid_fan for better performance
+function tropical_linear_space(A::QQMatrix, nu::TropicalSemiringMap{QQField, Nothing, minOrMax}=tropical_semiring_map(QQ); weighted_polyhedral_complex_only::Bool=false) where {minOrMax<:Union{typeof(min),typeof(max)}}
+    # compute reduced row echelon form of A
+    # and remove all zero rows so that matrix is of full rank
+    _,A = rref(A)
+    nonzeroRowIndices = findall(!iszero,[A[i,:] for i in 1:nrows(A)])
+    A = A[nonzeroRowIndices,:]
+
+    # Compute Bergman fan in polymake and convert to polyhedral complex
+    Sigma = PolyhedralComplex{QQFieldElem}(Polymake.tropical.matroid_fan{convention(nu)}(A))
+
+    # Add the missing all ones vector to the lineality
+    Sigma = add_missing_lineality_from_polymake(Sigma)
+    multiplicities = ones(ZZRingElem, n_maximal_polyhedra(Sigma))
+    TropL = tropical_linear_space(Sigma,multiplicities,convention(nu))
+    if !weighted_polyhedral_complex_only
+        set_attribute!(TropL,:algebraic_matrix,A)
+        set_attribute!(TropL,:tropical_semiring_map,nu)
+    end
+    return TropL
+end
 
+# helper function that constructs a basis of the vanishing set of a linear ideal
+# returns a matrix whose rows are the basis vectors
+function basis_of_vanishing_of_linear_ideal(I::MPolyIdeal)
+    # check that the input ideal is homogeneous and of degree 1
+    @req all(isone, total_degree.(Iterators.flatten(monomials.(gens(I))))) "generators of ideal must be linear and homogeneous"
+
+    x = gens(base_ring(I))
+    G = gens(I)
+    macaulayMatrix = matrix([[coeff(g,xi) for xi in x] for g in G])
+    A = transpose(kernel(macaulayMatrix, side = :right))
+
+    return A
+end
 
 @doc raw"""
     tropical_linear_space(I::MPolyIdeal, nu::TropicalSemiringMap; weighted_polyhedral_complex_only::Bool=false)
 
 Return the tropicalization of the vanishing set of `I` with respect to the tropical semiring map `nu`.  Requires the generators of `I` to be linear and homogeneous.  If `weighted_polyhedral_complex==true`, will not cache any extra information.
 
 # Examples
-
 ```jldoctest
 julia> R,(x1,x2,x3,x4) = polynomial_ring(QQ,4);
 
@@ -308,16 +354,7 @@ Min tropical linear space
 ```
 """
 function tropical_linear_space(I::MPolyIdeal, nu::Union{Nothing,TropicalSemiringMap}=nothing; weighted_polyhedral_complex_only::Bool=false)
-    # initialize nu as the trivial valuation if not specified by user
-    isnothing(nu) && (nu=tropical_semiring_map(coefficient_ring(I)))
-
-    # check that the input ideal is homogeneous and of degree 1
-    @req all(isone, total_degree.(Iterators.flatten(monomials.(gens(I))))) "generators of ideal must be linear and homogeneous"
-
-    x = gens(base_ring(I))
-    G = gens(I)
-    macaulayMatrix = matrix([[coeff(g,xi) for xi in x] for g in G])
-    A = transpose(kernel(macaulayMatrix, side = :right))
+    A = basis_of_vanishing_of_linear_ideal(I)
     TropL = tropical_linear_space(A,nu,
                                   weighted_polyhedral_complex_only=weighted_polyhedral_complex_only)
     if !weighted_polyhedral_complex_only
@@ -334,15 +371,37 @@ end
 #  Properties
 #
 ###############################################################################
-
 @doc raw"""
     pluecker_indices(TropL::TropicalLinearSpace)
 
 Return the Pluecker indices used to construct `TropL`.  Raises an error, if it is not cached.
 """
 function pluecker_indices(TropL::TropicalLinearSpace)
-    @req has_attribute(TropL,:pluecker_indices) "no pluecker indices cached"
-    return get_attribute(TropL,:pluecker_indices)
+    if has_attribute(TropL,:pluecker_indices)
+        return get_attribute(TropL,:pluecker_indices)
+    end
+
+    # there are two ways to compute the pluecker indices:
+    # - from the algebraic matrix
+    # - from the tropical matrix
+    # (we generally assume that a tropical linear space has at most one of the two)
+
+    # check if tropical matrix is available and, if yes, compute the pluecker indices from it
+    if has_attribute(TropL,:tropical_matrix)
+        A = get_attribute(TropL,:tropical_matrix)
+        plueckerIndices,plueckerVector = compute_pluecker_indices_and_vector(A)
+        set_attribute!(TropL,:pluecker_indices,plueckerIndices)
+        set_attribute!(TropL,:tropical_pluecker_vector,plueckerVector)
+        return plueckerIndices
+    end
+
+    # otherwise, get the algebraic matrix and compute the pluecker indices from it
+    # (will rightfully error if algebraic matrix cannot be found or computed)
+    A = algebraic_matrix(TropL)
+    plueckerIndices,plueckerVector = compute_pluecker_indices_and_vector(A)
+    set_attribute!(TropL,:pluecker_indices,plueckerIndices)
+    set_attribute!(TropL,:algebraic_pluecker_vector,plueckerVector)
+    return plueckerIndices
 end
 
 
@@ -352,8 +411,16 @@ end
 Return the tropical Pluecker vector of `TropL`.  Raises an error, if it is not cached.
 """
 function tropical_pluecker_vector(TropL::TropicalLinearSpace)
-    @req has_attribute(TropL,:tropical_pluecker_vector) "no tropical pluecker vector cached"
-    return get_attribute(TropL,:tropical_pluecker_vector)
+    if has_attribute(TropL,:tropical_pluecker_vector)
+        return get_attribute(TropL,:tropical_pluecker_vector)
+    end
+
+    A = tropical_matrix(TropL)
+    plueckerIndices,plueckerVector = compute_pluecker_indices_and_vector(A)
+    set_attribute!(TropL,:pluecker_indices,plueckerIndices)
+    set_attribute!(TropL,:tropical_pluecker_vector,plueckerVector)
+
+    return plueckerVector
 end
 
 
@@ -363,8 +430,16 @@ end
 Return the Pluecker vector over a valued field used to construct `TropL`.  Raises an error, if it is not cached.
 """
 function algebraic_pluecker_vector(TropL::TropicalLinearSpace)
-    @req has_attribute(TropL,:algebraic_pluecker_vector) "no algebraic pluecker vector cached"
-    return get_attribute(TropL,:algebraic_pluecker_vector)
+    if has_attribute(TropL,:algebraic_pluecker_vector)
+        return get_attribute(TropL,:algebraic_pluecker_vector)
+    end
+
+    A = algebraic_matrix(TropL)
+    plueckerIndices,plueckerVector = compute_pluecker_indices_and_vector(A)
+    set_attribute!(TropL,:pluecker_indices,plueckerIndices)
+    set_attribute!(TropL,:algebraic_pluecker_vector,plueckerVector)
+
+    return plueckerVector
 end
 
 
@@ -396,8 +471,14 @@ end
 Return the matrix over a valued field used to construct `TropL`.  Raises an error, if it is not cached.
 """
 function algebraic_matrix(TropL::TropicalLinearSpace)
-    @req has_attribute(TropL,:algebraic_matrix) "no algebraic matrix cached"
-    return get_attribute(TropL,:algebraic_matrix)
+    if has_attribute(TropL,:algebraic_matrix)
+        return get_attribute(TropL,:algebraic_matrix)
+    end
+
+    I = algebraic_ideal(TropL)
+    A = basis_of_vanishing_of_linear_ideal(I)
+    set_attribute!(TropL,:algebraic_matrix,A)
+    return A
 end
 
 

test/TropicalGeometry/intersection.jl

@@ -17,5 +17,4 @@
         @test dim(TropLLL)<0
     end
 
-
 end

test/TropicalGeometry/linear_space.jl

@@ -0,0 +1,35 @@
+@testset "src/TropicalGeometry/linear_space.jl" begin
+
+    @testset "det(A::Generic.MatSpaceElem{<:TropicalSemiringElem})" begin
+        M = matrix(QQ,[0 -3 0 -3 0 3 9 9 4 -9 -10 -6;
+                       -8 -7 -10 -10 0 -4 0 -2 -8 0 1 0;
+                       -2 7 0 1 -5 0 -5 0 0 -10 -3 6;
+                       9 2 -6 5 0 -6 -5 -6 -2 -2 4 2;
+                       -1 0 5 7 9 1 -5 2 -6 0 0 0;
+                       0 -9 8 3 -10 -9 4 -2 8 -9 -6 2])
+
+        TropL = tropical_linear_space(M)
+        @test dim(TropL) == rank(M)
+        @test ambient_dim(TropL) == ncols(M)
+        @test lineality_dim(TropL) == 1
+    end
+
+    @testset "tropical linear space attributes" begin
+        R,(w,x,y,z) = QQ["w","x","y","z"]
+        I = ideal([w+x+y+z])
+        TropL = tropical_linear_space(I)
+        @test has_attribute(TropL, :tropical_pluecker_vector)
+
+        plueckerVector = algebraic_pluecker_vector(TropL)
+        delete!(TropL.__attrs, :algebraic_pluecker_vector)
+        @test issetequal(plueckerVector,algebraic_pluecker_vector(TropL))
+
+        TT = tropical_semiring()
+        M = matrix(TT,[1 2 3; 2 3 4])
+        TropL = tropical_linear_space(M)
+        plueckerVector = tropical_pluecker_vector(TropL)
+        delete!(TropL.__attrs, :tropical_pluecker_vector)
+        @test issetequal(plueckerVector,tropical_pluecker_vector(TropL))
+    end
+
+end

test/book/specialized/holt-ren-tropical-geometry/linearSpace1.jlcon

@@ -9,15 +9,15 @@ Min tropical linear space
 
 julia> vertices_and_rays(TropL)
 5-element SubObjectIterator{Union{PointVector{QQFieldElem}, RayVector{QQFieldElem}}}:
- [0, 0, 0, 0]
- [0, -1, -1, -1]
  [0, 1, 0, 0]
+ [0, -1, -1, -1]
  [0, 0, 1, 0]
  [0, 0, 0, 1]
+ [0, 0, 0, 0]
 
 julia> maximal_polyhedra(IncidenceMatrix,TropL)
 4×5 IncidenceMatrix
-[1, 2]
-[1, 3]
-[1, 4]
 [1, 5]
+[2, 5]
+[3, 5]
+[4, 5]