Improve edges/ridges for simple/simplicial polytopes

[kliem]

In #30040 we have improved the face iterator for simple/simplicial polytopes, but have not applied this to our heuristic, whether we should find edges/ridges in dual/non-dual mode.

This is noticeable for edges of the permutahedron (and theoretically also for other edges of simple polytopes) and does not make a difference otherwise.

Before:

sage: C = polytopes.hypercube(12).combinatorial_polyhedron()
sage: %time _ = C.edges()
CPU times: user 86.1 ms, sys: 9 µs, total: 86.1 ms
Wall time: 86 ms
sage: %time _ = C.ridges()
CPU times: user 617 µs, sys: 3 µs, total: 620 µs
Wall time: 622 µs
sage: C = polytopes.cyclic_polytope(6, 20).combinatorial_polyhedron()
sage: %time _ = C.edges()
CPU times: user 114 µs, sys: 0 ns, total: 114 µs
Wall time: 117 µs
sage: %time _ = C.ridges()
CPU times: user 1.09 ms, sys: 5 µs, total: 1.09 ms
Wall time: 1.1 ms
sage: C = polytopes.permutahedron(7).combinatorial_polyhedron()
sage: %time _ = C.edges()
CPU times: user 16.1 s, sys: 0 ns, total: 16.1 s
Wall time: 16.1 s
sage: %time _ = C.ridges()
CPU times: user 2.35 ms, sys: 0 ns, total: 2.35 ms
Wall time: 2.36 ms

After:

sage: C = polytopes.hypercube(12).combinatorial_polyhedron()
sage: %time _ = C.edges()
CPU times: user 39.5 ms, sys: 0 ns, total: 39.5 ms
Wall time: 39.4 ms
sage:  %time _ = C.ridges()
CPU times: user 725 µs, sys: 60 µs, total: 785 µs
Wall time: 791 µs
sage: C = polytopes.cyclic_polytope(6, 20).combinatorial_polyhedron()
sage: %time _ = C.edges()
CPU times: user 154 µs, sys: 12 µs, total: 166 µs
Wall time: 168 µs
sage: %time _ = C.ridges()
CPU times: user 1.16 ms, sys: 95 µs, total: 1.25 ms
Wall time: 1.26 ms
sage: C = polytopes.permutahedron(7).combinatorial_polyhedron()
sage: %time _ = C.edges()
CPU times: user 63.2 ms, sys: 0 ns, total: 63.2 ms
Wall time: 63.2 ms
sage: %time _ = C.ridges()
CPU times: user 4.56 ms, sys: 0 ns, total: 4.56 ms
Wall time: 4.57 ms

The heuristic isn't perfect, but if it fails, you can now select a better algorithm:

sage: P = polytopes.permutahedron(7) 
sage: P1 = P.stack(next(P.face_generator(1)))
sage: %time P1.vertex_adjacency_matrix()
CPU times: user 20.6 s, sys: 7.96 ms, total: 20.6 s
Wall time: 20.6 s
5041 x 5041 dense matrix over Integer Ring (use the '.str()' method to see the entries)
sage: P1 = P.stack(next(P.face_generator(1)))
sage: %time P1.vertex_adjacency_matrix(algorithm='primal')
CPU times: user 206 ms, sys: 16 ms, total: 222 ms
Wall time: 222 ms
5041 x 5041 dense matrix over Integer Ring (use the '.str()' method to see the entries)

Depends on #33644
Depends on #33646

CC: @tscrim @yuan-zhou

Component: geometry

Author: Jonathan Kliem

Branch/Commit: 3bccc27

Reviewer: Travis Scrimshaw, Yuan Zhou

Issue created by migration from https://trac.sagemath.org/ticket/33611

[kliem]

Branch: u/gh-kliem/improved_edges_ridges_heuristic

[kliem]

New commits:

a0b541b

apply #30040 to edges/ridges heuristic

[kliem]

Changed branch from u/gh-kliem/improved_edges_ridges_heuristic to u/gh-kliem/improved_edge_ridge_heuristic

[kliem]

Commit: a0b541b

[yuan-zhou]

comment:5

Maybe irrelevant to this ticket, but are there better criteria than n_Vrepresentation > (n_facets)² (or n_Vrepresentation² < n_facets) in which cases one knows that the primal (or the dual) algorithm is faster?

"Wild estimates" plus special treatment for simple/simplicial polytopes does not seem optimal to me.

If we don't know a better criterion, does it make sense to leave the choice to the user, as an optional input?

[tscrim]

Reviewer: Travis Scrimshaw

[tscrim]

comment:6

I also think it would be good to answer that question and +1 for adding an option for the user to select. However, that can be left to another ticket.

For the change in this ticket. LGTM modulo the order of the checks. How expensive is the is_simple() compared to the other and how frequently do you expect it to get to that test? This is quite a micro-optimization, so if you don't want to be bothered with it, you can set a positive review.

[kliem]

comment:7

I think it is a good idea to expose the primal/dual or automatic options. I really don't know, when the algorithm is faster which way. Like I know for fixed polyhedra, but I have no clue in general. It is surprising how the primal algorithm even to generate edges is much better than the dual algorithm in many cases.

What do you think is better:

• a keyword algorithm, with options 'primal', 'dual' or None or
• a kewword dual with options True, False and None?

At the moment the Polyhedron_base method face_generator already exposes the keyword dual. But if this is not intuitive we should maybe switch to algorithm= instead.

[kliem]

comment:9

is_simple() is not cached (but could be). It has complexity of the incidence matrix, which is worse than the other checks, which are pretty much constant at this stuff needs to be computed anyway.

Maybe it makes sense to cache is_simple:

sage: C = polytopes.hypercube(12).combinatorial_polyhedron()
sage: %timeit C.is_simple()
24.5 µs ± 1.72 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
sage: %time C.f_vector()
CPU times: user 253 ms, sys: 11.7 ms, total: 265 ms
Wall time: 187 ms
(1, 4096, 24576, 67584, 112640, 126720, 101376, 59136, 25344, 7920, 1760, 264, 24, 1)

As for the complexity of the primal and dual algorithm (for edges):

• The dual algorithm takes about n_Vrepresentation ** 3 * n_facets.
• The primal algorithm in in the non-simple case takes sum(f_vector) * n_facets ** 2 * n_Vrepresentation.
• The primal algorithm in in the simple case takes sum(f_vector) * n_facets * n_Vrepresentation.

To estimate sum(f_vector) one can use the upper bound theorem. I came up with something in the complexity of dimension * binomial(min(n_facets, n_Vrepresentation), dimension // 2).

One could use those things to get a better estimate, but I don't know if this would be any good. Exposing the algorithm might be the better choice here.

[yuan-zhou]

Changed reviewer from Travis Scrimshaw to Travis Scrimshaw, Yuan Zhou

[yuan-zhou]

comment:10

Thank you for providing the complexity details! The formulas look very neat. I think it would be nice to use them in the code, especially when _f_vector was cached, because otherwise the users are unlikely to know such things.

I like the keyword dual with options True, False and None.

+1 for caching is_simple and is_simplicial.

+1 for "This is quite a micro-optimization, so if you don't want to be bothered with it, you can set a positive review."

[tscrim]

comment:11

comment:7: I would use algorithm as it is more flexible for the future.

comment:9: Thank you. I agree with Yuan that they should be cached. From what you're saying, I agree that they should be last in the checking order. Could there be a typo in your formulas? I don't understand what *** is. Perhaps this is my naïvity, but aren't facets in bijection with a minimal number of hyperplanes defining the polytope? Then I would expect the dual algorithm to almost always be the faster algorithm since squares grow much slow than binomials (after dividing by common factors). Unless the sum(f_vector) grows much slower than the upper bound theorem implies. (Sorry if these are more basic facts you need to explain to me; it's been a little while since I have done anything nontrivial in this area).

[kliem]

Dependencies: #33644

[kliem]

comment:13

Replying to @tscrim:

comment:7: I would use algorithm as it is more flexible for the future.

I also think this is a better name, but it involves deprecation. But well, that's life I guess.

comment:9: Thank you. I agree with Yuan that they should be cached. From what you're saying, I agree that they should be last in the checking order. Could there be a typo in your formulas?

I don't understand what *** is.

A typo, fixed now. Also the simple/simplical complexity had a typo.

Perhaps this is my naïvity, but aren't facets in bijection with a minimal number of hyperplanes defining the polytope? Then I would expect the dual algorithm to almost always be the faster algorithm since squares grow much slow than binomials (after dividing by common factors). Unless the sum(f_vector) grows much slower than the upper bound theorem implies. (Sorry if these are more basic facts you need to explain to me; it's been a little while since I have done anything nontrivial in this area).

A couple of simple/simplicial examples:

sage: C = polytopes.hypercube(10).combinatorial_polyhedron()
sage: C.f_vector()
(1, 1024, 5120, 11520, 15360, 13440, 8064, 3360, 960, 180, 20, 1)
sage: %time sum(1 for _ in C.face_iter(1, dual=True))
CPU times: user 175 ms, sys: 0 ns, total: 175 ms
Wall time: 175 ms
5120
sage: %time sum(1 for _ in C.face_iter(1, dual=False))
CPU times: user 3.71 ms, sys: 0 ns, total: 3.71 ms
Wall time: 3.71 ms
5120
sage: C.n_vertices() ** 3 * C.n_facets()  # complexity dual
21474836480
sage: sum(C.f_vector()) * C.n_facets() * C.n_vertices()  # complexity primal
1209344000

sage: C = polytopes.hypercube(12).combinatorial_polyhedron()
sage: C.f_vector()
(1, 4096, 24576, 67584, 112640, 126720, 101376, 59136, 25344, 7920, 1760, 264, 24, 1)
sage: %time sum(1 for _ in C.face_iter(1, dual=True))
CPU times: user 9.89 s, sys: 154 µs, total: 9.89 s
Wall time: 9.89 s
24576
sage: %time sum(1 for _ in C.face_iter(1, dual=False))
CPU times: user 36.7 ms, sys: 0 ns, total: 36.7 ms
Wall time: 36.3 ms
24576
sage: C.n_vertices() ** 3 * C.n_facets()  # complexity dual
1649267441664
sage: sum(C.f_vector()) * C.n_facets() * C.n_vertices()  # complexity primal
52242874368

sage: C = polytopes.cyclic_polytope(10, 20).combinatorial_polyhedron()
sage: C = C.polar()
sage: C.f_vector()
(1, 4004, 20020, 42315, 49140, 34470, 15504, 4845, 1140, 190, 20, 1)
sage: %time sum(1 for _ in C.face_iter(1, dual=True))
CPU times: user 9.55 s, sys: 18 µs, total: 9.55 s
Wall time: 9.55 s
20020
sage: %time sum(1 for _ in C.face_iter(1, dual=False))
CPU times: user 19.9 ms, sys: 0 ns, total: 19.9 ms
Wall time: 19.8 ms
20020
sage: C.n_vertices() ** 3 * C.n_facets()  # complexity dual
1283843841280
sage: sum(C.f_vector()) * C.n_facets() * C.n_vertices()  # complexity primal
13745732000

A non-simple and non-simplicial case:

sage: P = polytopes.permutahedron(7) 
sage: P1 = P.stack(next(P.face_generator(1)))
sage: C = CombinatorialPolyhedron(P1)
sage: %time C.f_vector()
CPU times: user 184 ms, sys: 0 ns, total: 184 ms
Wall time: 28.9 ms
(1, 5041, 16251, 19761, 11144, 2860, 267, 1)
sage: %time sum(1 for _ in C.face_iter(1, dual=True))
CPU times: user 20.1 s, sys: 0 ns, total: 20.1 s
Wall time: 20.1 s
16251
sage: %time sum(1 for _ in C.face_iter(1, dual=False))
CPU times: user 28.1 ms, sys: 0 ns, total: 28.1 ms
Wall time: 27.9 ms
16251
sage: C.n_vertices() ** 3 * C.n_facets()  # complexity dual
34202775806907
sage: sum(C.f_vector()) * C.n_facets() ** 2 * C.n_vertices()  # complexity primal
19882385613774

Note that the primal algorithm appears to be much better than our estimate. The reason is that the farther we go into the tree, the better the actual runtime will differ from the worst case scenario.

I also find the runtimes quite counter-intuitive. It seems that obtaining the f-vector of a simple/simplicial polytope with this iterator is much better than using the h-vector method.

[tscrim]

comment:14

On the internal function, having it be an int is much faster than a string for checking cases. You can also change that without any deprecation.

Thank you for the timings. I don't really have any intuition for what should be faster. So I won't be much help and will simply defer to your (and Yuan's) judgement.

[kliem]

Changed dependencies from #33644 to #33644, #33646

[kliem]

comment:16

Replying to @tscrim:

On the internal function, having it be an int is much faster than a string for checking cases. You can also change that without any deprecation.

Thank you for the timings. I don't really have any intuition for what should be faster. So I won't be much help and will simply defer to your (and Yuan's) judgement.

In #33646 I propose to use the keyword algorithm and face_generator over face_iter, but without deprecation.

The problem with exponential size is the following: Given a polytope as finite intersection of closed halfspaces. The number of vertices can be expected to be exponential with respect to the number of facets.

However, if you have a non-simple and non-simplicial polytope things might be a lot different. Those polytopes are less-generic but certainly more interesting.

[kliem]

New commits:

d64cf60	cache is_simple and is_simplicial for CombinatorialPolyhedron
01f30c5	allow algorithm for Polyhedron_base.face_generator
88f4bcd	prefer face_generator over face_iter and algorithm= over dual=
ed15041	deprecate keyword dual
89f9278	Merge branch 'public/33646' of git://trac.sagemath.org/sage into public/33611-new
dd05cc2	expose algorithm
860d335	use better heuristic to select dual/non-dual algorithm for edges/ridges
eb9e554	fix two doctests

[kliem]

Changed commit from a0b541b to eb9e554

[kliem]

Changed branch from u/gh-kliem/improved_edge_ridge_heuristic to public/33611

[yuan-zhou]

comment:18

In docstring of facet_graph, incidences --> indices

If ``names`` is ``False``, the ``vertices`` of the graph  will
         be the incidences of the facets in the Hrepresentation.

[yuan-zhou]

comment:19

In the docstring of def edges(self, names=True, algorithm=None):,
``names`` -- boolean (default: ``True``); if ``True``, ... --> if False ?

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

d27aad6

typo

[sagetrac-git]

Changed commit from eb9e554 to d27aad6

[yuan-zhou]

comment:21

Maybe change all face iterator in the docstrings to face generator as well?

[sagetrac-git]

Changed commit from d27aad6 to f986979

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

f986979

use "face generator" in the documentation

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

3bccc27

fix documentation error

[sagetrac-git]

Changed commit from f986979 to 3bccc27

[yuan-zhou]

comment:24

Regarding

     sage: C.edges()[:2]
-    ((A vertex at (-1, -1, -1, 1), A vertex at (-1, -1, -1, -1)),
-     (A vertex at (-1, 1, -1, -1), A vertex at (-1, -1, -1, -1)))
+    ((A vertex at (1, -1, -1, -1), A vertex at (-1, -1, -1, -1)),
+     (A vertex at (-1, -1, -1, 1), A vertex at (-1, -1, -1, -1)))
     sage: C.edges(names=False)[:2]
-    ((14, 15), (10, 15))
+    ((6, 15), (14, 15))

Are the vertices and facets of a CombinatorialPolyhedron always stored in a fixed order?
I vaguely recall that there was an ordering problem for the representations of a Polyhedron with different backends. I'm not sure if this is relevant here.

---

New commits:

f986979	use "face generator" in the documentation
3bccc27	fix documentation error

[kliem]

comment:25

Replying to @yuan-zhou:

Regarding

     sage: C.edges()[:2]
-    ((A vertex at (-1, -1, -1, 1), A vertex at (-1, -1, -1, -1)),
-     (A vertex at (-1, 1, -1, -1), A vertex at (-1, -1, -1, -1)))
+    ((A vertex at (1, -1, -1, -1), A vertex at (-1, -1, -1, -1)),
+     (A vertex at (-1, -1, -1, 1), A vertex at (-1, -1, -1, -1)))
     sage: C.edges(names=False)[:2]
-    ((14, 15), (10, 15))
+    ((6, 15), (14, 15))

Are the vertices and facets of a CombinatorialPolyhedron always stored in a fixed order?
I vaguely recall that there was an ordering problem for the representations of a Polyhedron with different backends. I'm not sure if this is relevant here.

---

New commits:

f986979	use "face generator" in the documentation
3bccc27	fix documentation error

The order of Vrepresentation and Hrepresentation is taken from the input and is then fixed.
What is happening here is that the default algorithm for the edges changed and then this changed the order edges.

[yuan-zhou]

comment:26

Thank you very much! All looks good to me.

@tscrim: Would you mind taking a look at the Patchbot failing? (I don't understand it.)

[kliem]

comment:27

At the moment there is only one failure and looks like a mathematica time out, so unrelated.

[tscrim]

comment:29

Indeed, it is unrelated. The latest patchbot run also came back morally green as well.

[kliem]

comment:30

Thank you.

[vbraun]

Changed branch from public/33611 to 3bccc27