Opened 6 years ago

Closed 6 years ago

# Faster Poyhedron.graph()

Reported by: Owned by: ncohen major sage-6.8 geometry vdelecroix, dimpase, vbraun, novoselt Nathann Cohen Dima Pasechnik N/A 5223cda 5223cda9f15b2a417b550039b978a07e126f9053 #18779

With this branch, `polytopes.Gosset_3_21().graph()` takes 300ms, against 16s before.

The 'definition' used to compute adjacent vertices can be found in 10

And of course the following works:

```sage: g=graphs.GossetGraph()
sage: g.is_isomorphic(polytopes.Gosset_3_21().graph())
True
```

### comment:1 Changed 6 years ago by ncohen

• Branch set to public/18860

### comment:2 Changed 6 years ago by git

• Commit set to 3b5eabb9735a867678c2c6a13dc7700c24d82732

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

 ​4a4ee0c `trac #18779: polytopes.gosset_3_21 and graphs.GossetGraph` ​b498a9e `trac #18779: Reviewer's (dead right) comments` ​6527f4b `capitalise Gosset` ​3b5eabb `trac #18860: Faster Poyhedron.graph()`

### comment:4 follow-up: ↓ 5 Changed 6 years ago by dimpase

```+        # Why 26 ?? I was expecting 'd-1'. Are there some useless inequalities
+        # in the list?
```

well, you have 26 facets on each edge - none are more useless than the others. You can do the following (note that 0th and 1st vertices are adjacent).

```sage: p=polytopes.Gosset_3_21()
sage: M=p.incidence_matrix()
sage: a=0; b=1;
sage: vv=[i for i in [0..702] if M[a,i]*M[b,i]>0]
sage: matrix([p.inequalities()[j].vector() for j in vv]).T.kernel()
Free module of degree 9 and rank 1 over Integer Ring
Echelon basis matrix:
[0 1 0 0 0 0 0 0 0]
```

or just

```sage: m=matrix([p.inequalities()[j].vector() for j in vv])
sage: m.rank()
8
```

which makes sense, as `m` is of size 27x9. (27, as there is not full-dimensional).

General code for adjacency of two vertices `a` and `b` of a polytope in `d`-dimensional space should check that `m` has corank `1`.

(There could be special cases when one of the vertices is the origin, or some inequalities are homogeneous).

### comment:5 in reply to: ↑ 4 ; follow-up: ↓ 6 Changed 6 years ago by ncohen

General code for adjacency of two vertices `a` and `b` of a polytope in `d`-dimensional space should check that `m` has corank `1`.

Would you know how to adapt this code for that? Hoping that it will remain faster than what we currently have?

Nathann

### comment:6 in reply to: ↑ 5 ; follow-up: ↓ 7 Changed 6 years ago by dimpase

General code for adjacency of two vertices `a` and `b` of a polytope in `d`-dimensional space should check that `m` has corank `1`.

Would you know how to adapt this code for that? Hoping that it will remain faster than what we currently have?

My understanding is that hacking `hasse_diagram.py` to stop earlier is what one actually needs to do. Except that it goes from bottom to top, which probably means that one gets adjacency of facets as the first step. But it should be possible to change it that goes the other way around. We can write a paper about it ;-).

### comment:7 in reply to: ↑ 6 ; follow-up: ↓ 8 Changed 6 years ago by ncohen

My understanding is that hacking `hasse_diagram.py` to stop earlier is what one actually needs to do.

I do not know how to do that. If you know how to fix this implementation, at least it will be done. Otherwise chances are that nothing will happen.

Nathann

### comment:8 in reply to: ↑ 7 Changed 6 years ago by dimpase

My understanding is that hacking `hasse_diagram.py` to stop earlier is what one actually needs to do.

I do not know how to do that. If you know how to fix this implementation, at least it will be done. Otherwise chances are that nothing will happen

the whole Polyhedron code does a potentially very slow thing: enumerating the inequalities at construction time, instead of doing it optionally or/and lazily.

Indeed, e.g. the problem at hand, deciding if two vertices are adjacent, does not need pre-computed inequalities. Namely, checking that v_1 and v_2 are adjacent can be done by solving the following LP:

```    variables: x_0..x_d,x_{d+1}
objective: max x_0
constraints:
-1 <= x_i <=1 for i in [1..d];
x_{d+1} >= 0;
sum_{j=1}^d v_{i,j}*x_j = x_{d+1},         for i=1,2;
sum_{j=1}^d v_{i,j}*x_j >= x_{d+1} + x_0,  for i>2;
```

if the optimum exists and is strictly positive, they form an edge, and otherwise they don't.

So we are trying to find an inequality that is valid for the convex hull of v_k's and turns into an equation on v_1 and v_2, but is a strict inequality on the rest of v_k's.

(The LP above will not work if there are more v_j's on the line joining v_1 and v_2 --- but this case can be taken care of trivially.)

### comment:9 Changed 6 years ago by ncohen

Oh. That's a nice formulation of the problem, thanks. I will not forget it.

Would you know a combinatorial way to compute it? It seems that we can already rely on cddlib (see #18861) in order to build the graph over RDF.

Nathann

Last edited 6 years ago by ncohen (previous) (diff)

### comment:10 follow-up: ↓ 19 Changed 6 years ago by ncohen

Here is a new 'combinatorial' version. It uses an 'idea', which I am not sure is a proper definition of the object, but well... You tell me `:-P`

1) For every inequality `I`, build the set `set(I)` of points for which it is an equality 2) For any pair of points `i,j`, list all inequalities `I_1, ..., I_k` that contain them both 3) Intersection `set(I_1), ..., set(I_k)`. If this intersection has cardinality 2, then ij is an edge.

Works for Gosset (at least) `:-P`

Nathann

### comment:11 Changed 6 years ago by git

• Commit changed from 3b5eabb9735a867678c2c6a13dc7700c24d82732 to d084b1f2d0a7fa39fd522079aba39494f1704526

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

 ​d084b1f `trac #18860: Faster Poyhedron.graph()`

### comment:12 follow-up: ↓ 21 Changed 6 years ago by ncohen

• Description modified (diff)

A bit cleaner: I moved the code to the `.graph()`, method for apparently the `_vertex_adjacency_matrix` is not what it claims to be. Its dimension may not even be the number of vertices (if you have rays n your polyhedron) [1].

With this, all tests pass in the `geometry` folder.

Nathann

[1] You can see what goes wrong by adding `return self.graph().adjacency_matrix()` in the first line of `_vertex_adjacency_matrix`.

Last edited 6 years ago by ncohen (previous) (diff)

### comment:13 Changed 6 years ago by ncohen

• Status changed from new to needs_review

### comment:14 Changed 6 years ago by git

• Commit changed from d084b1f2d0a7fa39fd522079aba39494f1704526 to 74d8412010748142918edf6178636bffdb14e60a

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

 ​74d8412 `trac #18860: Faster Poyhedron.graph()`

### comment:16 in reply to: ↑ 15 Changed 6 years ago by ncohen

I was not talking about the `.adjacency_matrix` method. I was talking about `._vertex_adjacency_matrix`, which should (from its name) have as many rows/cols as there are vertices.

Nathann

### comment:17 follow-up: ↓ 18 Changed 6 years ago by vbraun

Again its really the V-representation adjacency matrix; The non-compact case is just the projectivization. So in a sense its really the vertex adjacency matrix; I know its somewhat confusing but most people are interested in the compact case so being overly nitpicky about wording is going to be confusing too.

### comment:18 in reply to: ↑ 17 Changed 6 years ago by ncohen

Again its really the V-representation adjacency matrix; The non-compact case is just the projectivization. So in a sense its really the vertex adjacency matrix; I know its somewhat confusing but most people are interested in the compact case so being overly nitpicky about wording is going to be confusing too.

So was this line incorrect `return Graph(self.vertex_adjacency_matrix(), loops=False)` ?

I expect the `.vertex_graph()` to give me a graph whose vertices are the vertices of the polyhedron. And the incidence matrix of that graph is *not* the output of the `_vertex_adjacency_matrix` whose rows are not the vertices only. Do we agree?

Nathann

### comment:19 in reply to: ↑ 10 ; follow-up: ↓ 20 Changed 6 years ago by dimpase

Here is a new 'combinatorial' version. It uses an 'idea', which I am not sure is a proper definition of the object, but well... You tell me `:-P`

1) For every inequality `I`, build the set `set(I)` of points for which it is an equality 2) For any pair of points `i,j`, list all inequalities `I_1, ..., I_k` that contain them both 3) Intersection `set(I_1), ..., set(I_k)`. If this intersection has cardinality 2, then ij is an edge.

yes, that's correct for the polytopes - provided there is no redundancy. E.g. if a polytope is not full-dimensional, then adding an equation satisfied by all the vertices to an inequality is a "copy" of this inequality, vertex-vise.

Works for Gosset (at least) `:-P`

Nathann

### comment:20 in reply to: ↑ 19 Changed 6 years ago by ncohen

yes, that's correct for the polytopes - provided there is no redundancy.

Excellent news! Then we have a very nice speedup ahead `:-PPPP`

As for redundancy, it is true that this algorithm will not behave very well if there are several vertices at the same place which are actually the same. Redundancy of inequalities, however, is not a problem.

E.g. if a polytope is not full-dimensional, then adding an equation satisfied by all the vertices to an inequality is a "copy" of this inequality, vertex-vise.

Yeah but in terms of sets it produces the very same one. So it's cool.

Nathann

### comment:21 in reply to: ↑ 12 ; follow-ups: ↓ 22 ↓ 23 Changed 6 years ago by dimpase

A bit cleaner: I moved the code to the `.graph()`, method for apparently the `_vertex_adjacency_matrix` is not what it claims to be. Its dimension may not even be the number of vertices (if you have rays n your polyhedron) [1].

With this, all tests pass in the `geometry` folder.

Aren't we changing the behaviour of some functions for the non-compact case?

### comment:22 in reply to: ↑ 21 Changed 6 years ago by ncohen

Aren't we changing the behaviour of some functions for the non-compact case?

This patch only touches the function `.vertex_graph`. Its behaviour changes indeed for the non-compact case, but unless you can defend that the previous behavior was not a bug (there were more vertices in the graph than vertices in the polyhedron) I would say that this is not a problem.

Nathann

### comment:23 in reply to: ↑ 21 ; follow-up: ↓ 24 Changed 6 years ago by ncohen

Aren't we changing the behaviour of some functions for the non-compact case?

also note that not a single doctest changed. But to be honest, I still believe that the `._vertex_adjacency_matrix` should be either fixed or renamed. It is not what it claims to be.

Nathann

### comment:24 in reply to: ↑ 23 ; follow-up: ↓ 25 Changed 6 years ago by dimpase

Aren't we changing the behaviour of some functions for the non-compact case?

also note that not a single doctest changed. But to be honest, I still believe that the `._vertex_adjacency_matrix` should be either fixed or renamed. It is not what it claims to be.

presently we have

```adjacency_matrix = vertex_adjacency_matrix
```

in `geometry/polhedron/base.py`. I'd think we should rename `._vertex_adjacency_matrix` to something like `._Vrep_adjacency_matrix`, and, dually `._facet_adjacency_matrix` to `._Hrep_adjacency_matrix`.

And then set `adjacency_matrix = _Vrep_adjacency_matrix`, and provide a meaningful implementation of `vertex_adjacency_matrix`, as you propose here.

### comment:25 in reply to: ↑ 24 ; follow-up: ↓ 26 Changed 6 years ago by ncohen

Seems reasonable. Though I would very much like to see the definition of a `V_rep adjacency matrix`.

Nathann

### comment:26 in reply to: ↑ 25 Changed 6 years ago by dimpase

Seems reasonable. Though I would very much like to see the definition of a `V_rep adjacency matrix`.

well, check out `sage.geometry.polyhedron.constructor?`. In short, any polyhedron is the Minkowski sum of a polytope P (spanned by vertices), a pointed polyhedral cone (spanned by rays), and a subspace S (spanned by lines).

• If S is empty or of dimension bigger than 1, it can be ignored.
• Otherwise S is a line, and it is adjacent to all the vertices, and to no ray; and, counter-intuitively, no vertices are adjacent to each other.
• Rays are not adjacent to each other, and might be adjacent to some vertices.

How these general adjacencies are used, I don't know.

### comment:27 Changed 6 years ago by ncohen

Oh. Then it seems that it will be possible to rewrite that method using the algorithm implemented here, then. And everything will be faster.

Nathann

### comment:28 follow-up: ↓ 29 Changed 6 years ago by dimpase

Are we still in `needs review` here? I struggle to remember what was the sticking point.

### comment:29 in reply to: ↑ 28 Changed 6 years ago by ncohen

Are we still in `needs review` here? I struggle to remember what was the sticking point.

I do not think that there was any. We had discussions on the name of other functions, but this new implementation only touches the "vertex graph" and clearly improves the running time.

Nathann

### comment:30 Changed 6 years ago by dimpase

• Reviewers set to Dima Pasechnik
• Status changed from needs_review to positive_review

OK. Merges in 6.8, no problem...

Thanks!

### comment:32 Changed 6 years ago by vbraun

• Status changed from positive_review to needs_work
```sage -t --long src/sage/graphs/generators/random.py
**********************************************************************
File "src/sage/graphs/generators/random.py", line 806, in sage.graphs.generators.random.RandomTriangulation
Failed example:
for i in range(10):
g = graphs.RandomTriangulation(30,embed=True)
assert g.is_planar() and g.size() == 3*g.order()-6
Exception raised:
Traceback (most recent call last):
File "/mnt/highperf/buildbot/slave/sage_git/build/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 496, in _run
self.compile_and_execute(example, compiler, test.globs)
File "/mnt/highperf/buildbot/slave/sage_git/build/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 858, in compile_and_execute
exec(compiled, globs)
File "<doctest sage.graphs.generators.random.RandomTriangulation[3]>", line 3, in <module>
assert g.is_planar() and g.size() == Integer(3)*g.order()-Integer(6)
File "/mnt/highperf/buildbot/slave/sage_git/build/local/lib/python2.7/site-packages/sage/graphs/generic_graph.py", line 4035, in is_planar
planar = is_planar(G,kuratowski=kuratowski,set_pos=set_pos,set_embedding=set_embedding)
File "sage/graphs/planarity.pyx", line 99, in sage.graphs.planarity.is_planar (/mnt/highperf/buildbot/slave/sage_git/build/src/build/cythonized/sage/graphs/planarity.c:1800)
g.relabel(ffrom)
File "/mnt/highperf/buildbot/slave/sage_git/build/local/lib/python2.7/site-packages/sage/graphs/generic_graph.py", line 18006, in relabel
new_attr[perm[v]] = value
KeyError: 0
**********************************************************************
1 of   5 in sage.graphs.generators.random.RandomTriangulation
[71 tests, 1 failure, 5.56 s]
```

Would have been easy to catch if you had checked the patchbot output...

### comment:33 Changed 6 years ago by git

• Commit changed from 74d8412010748142918edf6178636bffdb14e60a to 5223cda9f15b2a417b550039b978a07e126f9053

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

 ​5dfc776 `trac #18860: Merged with 6.9.beta0` ​5223cda `trac #18860: Remove vertex labels`

### comment:34 follow-up: ↓ 35 Changed 6 years ago by ncohen

• Status changed from needs_work to needs_review

### comment:35 in reply to: ↑ 34 Changed 6 years ago by dimpase

Could you explain what happens in the last commit? (E.g.. why did it work earlier?)

Last edited 6 years ago by dimpase (previous) (diff)

### comment:36 Changed 6 years ago by ncohen

Oh, no problem: there is a difference between the previous output or `.vertex_graph` and the new version: in the former, the vertices of the graph were labelled on integers. Now, the vertices of that graph are the 'vertex' objects contained in the polyhedron.

I adapted the code, for it expected the keys of the position dictionary (the vertices) to be integers, instead of the new points.

Actually, it convinced me that the previous version of this function was incorrect in theory. It was assumed that the vertices of the polyhedron appeared in the very same order as how they were first given to the polyhedron constructor. I don't know if it is true in practice, but for sure there is no reason to assume that a priori. So you can see this as a bugfix `;-)`

Nathann

### comment:37 Changed 6 years ago by vbraun

Vertices are definitely reordered in `Polyhedron`.

### comment:38 Changed 6 years ago by dimpase

• Status changed from needs_review to positive_review

OK, looks good.

### comment:39 Changed 6 years ago by vbraun

• Branch changed from public/18860 to 5223cda9f15b2a417b550039b978a07e126f9053
• Resolution set to fixed
• Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.