ToricVarieties: Interplay with polyhedral geometry and performance

[HereAround]

I am currently using Oscar to triangulate a (facette of a) polytope, construct all corresponding toric varieties, compute their Stanley-Reisner ideals and finally intersect these ideals. This seems to pose both a performance and convenience test. I currently use the following code:

#(1) Set variable names of our choice
vars_dict = Dict()
vars_dict[matrix(ZZ,[-1 -1 -1])] = "x0"
vars_dict[matrix(ZZ,[2 -1 -1])] = "x1"
vars_dict[matrix(ZZ,[-1 2 -1])] = "x2"
vars_dict[matrix(ZZ,[-1 -1 5])] = "x3"
vars_dict[matrix(ZZ,[-1 -1 0])] = "x4"
vars_dict[matrix(ZZ,[-1 -1 1])] = "x5"
vars_dict[matrix(ZZ,[-1 -1 2])] = "x6"
vars_dict[matrix(ZZ,[-1 -1 3])] = "x7"
vars_dict[matrix(ZZ,[-1 -1 4])] = "x8"
vars_dict[matrix(ZZ,[-1 0 -1])] = "x9"
vars_dict[matrix(ZZ,[-1 0 0])] = "x10"
vars_dict[matrix(ZZ,[-1 0 1])] = "x11"
vars_dict[matrix(ZZ,[-1 0 2])] = "x12"
vars_dict[matrix(ZZ,[-1 0 3])] = "x13"
vars_dict[matrix(ZZ,[-1 1 -1])] = "x14"
vars_dict[matrix(ZZ,[-1 1 0])] = "x15"
vars_dict[matrix(ZZ,[-1 1 1])] = "x16"
vars_dict[matrix(ZZ,[0 -1 -1])] = "x17"
vars_dict[matrix(ZZ,[0 -1 0])] = "x18"
vars_dict[matrix(ZZ,[0 -1 1])] = "x19"
vars_dict[matrix(ZZ,[0 -1 2])] = "x20"
vars_dict[matrix(ZZ,[0 -1 3])] = "x21"
vars_dict[matrix(ZZ,[0 0 -1])] = "x22"
vars_dict[matrix(ZZ,[0 0 1])] = "x23"
vars_dict[matrix(ZZ,[0 1 -1])] = "x24"
vars_dict[matrix(ZZ,[1 -1 -1])] = "x25"
vars_dict[matrix(ZZ,[1 -1 0])] = "x26"
vars_dict[matrix(ZZ,[1 -1 1])] = "x27"
vars_dict[matrix(ZZ,[1 0 -1])] = "x28"

#(2) Form point set to be triangulated
f = convex_hull( [-1 -1 -1; 2 -1 -1; -1 2 -1])
zero = [0 for i in 1:ambient_dim(f)]
V = [Vector{fmpq}(v) for v in lattice_points(f) if !iszero(v)]
pts = vcat(matrix(QQ, transpose(zero)), matrix(QQ, transpose(hcat(V...))))

#(3) Triangulate
trias = star_triangulations(pts; full=true, regular=true)
trias_julia = [[[c[k]-1 for k in 2:length(c)] for c in t] for t in trias]

#(4) Compute rays of all toric varieties
points = lattice_points(f)
points = matrix(ZZ, transpose(hcat(points...)))
fan_rays = zeros(Int, nrows(points), ncols(points))
for i in 1:nrows(points)
    for j in 1:ncols(points)
        fan_rays[i,j] = ZZ(points[i,j])
    end
end

#(5) Construct ring which contains all Stanley-Reisner ideals
max_cones = IncidenceMatrix(vcat(trias_julia[1]))
variety = NormalToricVariety(PolyhedralFan(fan_rays, max_cones))
reference_rays = matrix(ZZ, rays(variety))
my_vars = [vars_dict[reference_rays[i,:]] for i in 1:nrows(reference_rays)]
set_coordinate_names(variety, my_vars)
R = cox_ring(variety)

#(6) Compute all ideals
sr_ideals_julia = []
for n in 1:length(trias_julia)
    max_cones = IncidenceMatrix(vcat(trias_julia[n]))
    variety = NormalToricVariety(PolyhedralFan(fan_rays, max_cones))
    new_rays = matrix(ZZ, rays(variety))
    convert = Dict()
    for k in 1:nrows(new_rays)
        for l in 1:nrows(reference_rays)
            if new_rays[k,:] == reference_rays[l,:]
                convert[l] = k
            end
        end
    end
    mnf = Oscar._minimal_nonfaces(variety)
    nr = Polymake.nrows(mnf)
    nc = Polymake.ncols(mnf)
    my_gens = [R([1], [[Int(mnf[i, convert[j]]) for j in 1:nc]]) for i in 1:nr]
    push!(sr_ideals_julia, ideal(my_gens))
end

#(7) Finall intersect all ideals
I1 = intersect(sr_ideals_julia...)

Points regarding convenience:

1. I want to enforce variable names to be used for the indeterminates. As a consequence of https://github.com/oscar-system/Oscar.jl/issues/1208, this requires the use of the dictionaries vars_dict and convert. If the rays can be fixed (which I would very much appreciate), the code in step 6 could be simplified a lot. Namely, the for-loop would more or less just execute stanley_reisner_ideal(R, variety).
2. In step 4, I would ike to say fan_rays = lattice_points(f). Sadly, this is not accepted as input for PolyhedralFan. Maybe this can be changed?
3. In steps 2 and 3, I once add the origin by hand just so I can eventually remove it(s index) again in the list of triangulations obtained in step 3. Maybe the addition of the origin could be done internally? I am also wondering if it might make sense to remove the origin from the list of star triangulations before returning it. Maybe this requires to introduce a data_type StarTriangulation with an attribute star_center which could be set to the origin before return? (Aside: This might also be a good extension to allow for star centers other than the origin. But I do not have application of that in mind just yet.)

Regarding performance:

1. The above example code seems to run in about 1 to 2 minutes. That is ok.
2. If we replace f = convex_hull( [-1 -1 -1; 2 -1 -1; -1 2 -1]) by f = convex_hull( [-1 -1 -1; 2 -1 -1; -1 -1 5]), then this code still completes, but it seems to take between one and two hours. The bottle-neck are not the triangulations, but rather a combination of the computation of all Stanley-Reisner ideals and possibly their intersection. It feels like this could be improved.

Suggestions and hints very much appreciated.

[Edit: I will add more accurate run times shortly.]

[benlorenz]

Some simplifications for the matrix conversions:

#(2) Form point set to be triangulated
f = convex_hull( [-1 -1 -1; 2 -1 -1; -1 2 -1])
points = lattice_points(f)
V = filter(!iszero, points)
# using reduce instead of the splash operator (V...) is supposedly faster
pts = vcat(zero_matrix(ZZ, 1, ambient_dim(f)), matrix(ZZ, transpose(reduce(hcat,V))))

#(3) Triangulate
trias = star_triangulations(pts; full=true, regular=true)
trias_julia = [[[c[k]-1 for k in 2:length(c)] for c in t] for t in trias]

#(4) Compute rays of all toric varieties
fan_rays = matrix(ZZ, points)

#(5) Construct ring which contains all Stanley-Reisner ideals
max_cones = IncidenceMatrix(vcat(trias_julia[1]))
variety = NormalToricVariety(PolyhedralFan(fan_rays, max_cones))

lattice_points(P) returns an object encoding points while PolyhedralFan expects rays. So the best way to do this is to convert the points to a matrix, i.e. PolyhedralFan(matrix(ZZ,lattice_points(P)), some_cones) should work fine. (Internally the lattice points are also encoded as rays but with an additional homogenization coordinate, to avoid any confusion whether to keep or drop the homogenization there is no direct conversion)

PS: I improved the code display in your post to make it easier to read.

[HereAround]

Thank you @benlorenz (for both the comment and the editing)!

Any chance for a convenience method so that one can directly put the points into PolyhedralFan?

Here are a few timings

1. Oscar, f = convex_hull( [-1 -1 -1; 2 -1 -1; -1 2 -1]) and without the improvements by @benlorenz: BenchmarkTools finds 7.805 s (2357225 allocations: 40.11 MiB)
2. Oscar, f = convex_hull( [-1 -1 -1; 2 -1 -1; -1 2 -1]) and with the improvements by @benlorenz: BenchmarkTools finds 7.534 s (2361068 allocations: 39.87 MiB)
3. Oscar, f = convex_hull( [-1 -1 -1; 2 -1 -1; -1 -1 5]) and with the above improvement by @benlorenz: BenchmarkTools finds 2648.960 s (967571326 allocations: 21.53 GiB)

For comparison Sage, which can do the same thing of course.

1. Sage, f = convex_hull( [-1 -1 -1; 2 -1 -1; -1 2 -1]) and computation of triangulations: 1.79s
2. Sage, f = convex_hull( [-1 -1 -1; 2 -1 -1; -1 2 -1]) and reading triangulations from file: 1.78s
3. Sage, f = convex_hull( [-1 -1 -1; 2 -1 -1; -1 -1 5]) and computation of triangulations: 1050s
4. Sage, f = convex_hull( [-1 -1 -1; 2 -1 -1; -1 -1 5]) and reading triangulations from file: 875s

So for the simple computation, Sage is about 4 times faster. For the more involved one, Sage is about 3 times faster. While my code above is far from fully optimized, this still feels a lot.

I will try to investigate more and find out what the bottleneck is. If anybody has hint/suggestions or comments, that would be appreciated.

[benlorenz]

I had another look, this time more on the performance side: most of the time seems to be spent in the redundancy elimination of the fan in the loop. This is exactly what Alexej suggested in the other ticket with a non_redundant flag for the constructor, I have used such a construction (via a polymake object) in the code below. Using the same construction inside the loop with the remaining code as in the original post makes it already run a lot faster. But in this case the fan construction inside the loop should not even be necessary as the triangulation should be enough to compute the minimal nonfaces and the ideals.

#(2) Form point set to be triangulated
f = convex_hull( [-1 -1 -1; 2 -1 -1; -1 2 -1]);
#f = convex_hull( [-1 -1 -1; 2 -1 -1; -1 -1 5]);
points = lattice_points(f);
V = filter(!iszero, points);
# using reduce instead of the splash operator (V...) is supposedly faster
pts = vcat(zero_matrix(ZZ, 1, ambient_dim(f)), matrix(ZZ, transpose(reduce(hcat,V))));

#(3) Triangulate
@time trias = star_triangulations(pts; full=true, regular=true);
trias_julia = [[[c[k]-1 for k in 2:length(c)] for c in t] for t in trias];

#(4) Compute rays of all toric varieties
fan_rays = matrix(ZZ, points);

#(5) Construct ring which contains all Stanley-Reisner ideals
max_cones = IncidenceMatrix(vcat(trias_julia[1]));
# these two lines do the `non_redundant` construction mentioned in the other ticket
pmfan = Polymake.fan.PolyhedralFan(RAYS=fan_rays, MAXIMAL_CONES=max_cones);
variety = NormalToricVariety(PolyhedralFan{fmpq}(pmfan))
my_vars = [vars_dict[fan_rays[i,:]] for i in 1:nrows(fan_rays)];
set_coordinate_names(variety, my_vars);
R = cox_ring(variety);

#(6) Compute all ideals
sr_ideals_julia = []
@time for n in 1:length(trias_julia)
    id = stanley_reisner_ideal(R, SimplicialComplex(trias_julia[n]));
    push!(sr_ideals_julia,id)
end;

#(7) Finall intersect all ideals
@time I1 = reduce(intersect,sr_ideals_julia)

In my code the non-redundant fan construction is needed to make sure the indices of the triangulation match the indices of the variables which is constructed from the variety, this could probably be avoided by computing one triangulation at the beginning, then constructing the variety and the ring, and later computing the triangulations from the rays of the fan.

Some timings from my laptop:

1. Original code, redundant input, small example: 5.323264 seconds (2.33 M allocations: 37.987 MiB)

2. Original code but with non_redundant fan construction in the loop, small example: 0.251202 seconds (489.74 k allocations: 22.697 MiB)

3. Original code but with non_redundant fan construction in the loop, large example: 167.594815 seconds (310.76 M allocations: 15.486 GiB, 2.91% gc time, 0.04% compilation time)

4. My code, small example: 0.190431 seconds (361.13 k allocations: 13.841 MiB)

5. My code, large example: 144.285265 seconds (290.37 M allocations: 12.196 GiB, 3.20% gc time)

In the large example the time is distributed roughly as follows: 33 seconds for the triangulations via topcom, 64 seconds for the minimal-nonfaces computations in the loop, and 47 seconds for the final intersection.

[HereAround]

Thank you all. I think we can close this.