Faster Poyhedron.graph() #18860

nathanncohen · 2015-07-07T16:35:48Z

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

The 'definition' used to compute adjacent vertices can be found in [comment:10]

And of course the following works:

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

Depends on #18779

CC: @videlec @dimpase @vbraun @novoselt

Component: geometry

Author: Nathann Cohen

Branch/Commit: 5223cda

Reviewer: Dima Pasechnik

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

The text was updated successfully, but these errors were encountered:

nathanncohen · 2015-07-07T16:36:06Z

Branch: public/18860

sagetrac-git · 2015-07-07T16:36:47Z

Commit: 3b5eabb

sagetrac-git · 2015-07-07T16:36:47Z

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()`

dimpase · 2015-07-07T18:03:25Z

comment:4

+        # 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).

nathanncohen · 2015-07-07T19:29:55Z

comment:5

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

dimpase · 2015-07-08T09:10:05Z

comment:6

Replying to @nathanncohen:

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 ;-).

nathanncohen · 2015-07-08T11:41:25Z

comment:7

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

dimpase · 2015-07-08T14:24:27Z

comment:8

Replying to @nathanncohen:

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.)

nathanncohen · 2015-07-08T14:56:33Z

comment:9

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

nathanncohen · 2015-07-08T16:00:28Z

comment:10

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

For every inequality I, build the set set(I) of points for which it is an equality
For any pair of points i,j, list all inequalities I_1, ..., I_k that contain them both
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

sagetrac-git · 2015-07-08T16:00:55Z

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

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

sagetrac-git · 2015-07-08T16:00:55Z

Changed commit from 3b5eabb to d084b1f

nathanncohen · 2015-07-08T16:44:10Z

comment:12

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.

sagetrac-git · 2015-07-08T16:45:42Z

Changed commit from d084b1f to 74d8412

sagetrac-git · 2015-07-08T16:45:42Z

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

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

vbraun · 2015-07-08T19:47:46Z

comment:15

The adjacency matrix for unbound polyhedra is essentially the projectvization; See also the adjacency matrix docs for a more down-to-earth definition.

nathanncohen · 2015-07-08T20:00:29Z

comment:16

The adjacency matrix for unbound polyhedra is essentially the projectvization; See also the adjacency matrix docs for a more down-to-earth definition.

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

vbraun · 2015-07-08T20:55:30Z

comment:17

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.

nathanncohen · 2015-07-08T20:58:46Z

comment:18

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

dimpase · 2015-07-08T23:10:59Z

comment:19

Replying to @nathanncohen:

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

For every inequality I, build the set set(I) of points for which it is an equality

For any pair of points i,j, list all inequalities I_1, ..., I_k that contain them both

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

nathanncohen · 2015-07-09T07:21:39Z

comment:20

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

dimpase · 2015-07-09T10:37:59Z

comment:21

Replying to @nathanncohen:

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?

nathanncohen · 2015-07-09T10:44:08Z

comment:22

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

nathanncohen · 2015-07-09T10:45:53Z

comment:23

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

dimpase · 2015-07-09T11:08:29Z

comment:24

Replying to @nathanncohen:

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.

nathanncohen · 2015-07-09T11:11:29Z

comment:25

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

Nathann

dimpase · 2015-07-09T19:00:42Z

comment:26

Replying to @nathanncohen:

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.

nathanncohen · 2015-07-09T19:21:32Z

comment:27

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

dimpase · 2015-07-31T10:58:36Z

comment:28

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

nathanncohen · 2015-07-31T11:19:10Z

comment:29

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

dimpase · 2015-07-31T12:54:19Z

Reviewer: Dima Pasechnik

dimpase · 2015-07-31T12:54:19Z

comment:30

OK. Merges in 6.8, no problem...

nathanncohen · 2015-07-31T12:55:25Z

comment:31

Thanks!

vbraun · 2015-08-01T09:11:42Z

comment:32

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 item had failures:
   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...

sagetrac-git · 2015-08-01T09:27:47Z

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`

sagetrac-git · 2015-08-01T09:27:47Z

Changed commit from 74d8412 to 5223cda

dimpase · 2015-08-01T09:39:46Z

comment:35

Replying to @nathanncohen:

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

nathanncohen · 2015-08-01T09:44:51Z

comment:36

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

vbraun · 2015-08-01T09:47:40Z

comment:37

Vertices are definitely reordered in Polyhedron.

dimpase · 2015-08-01T10:14:05Z

comment:38

OK, looks good.

vbraun · 2015-08-01T23:07:52Z

Changed branch from public/18860 to 5223cda

nathanncohen mannequin added this to the sage-6.8 milestone Jul 7, 2015

nathanncohen mannequin added c: geometry labels Jul 7, 2015

This comment has been minimized.

Sign in to view

nathanncohen mannequin added the s: needs review label Jul 8, 2015

dimpase added s: positive review and removed s: needs review labels Jul 31, 2015

vbraun added s: needs work and removed s: positive review labels Aug 1, 2015

nathanncohen mannequin added s: needs review and removed s: needs work labels Aug 1, 2015

dimpase added s: positive review and removed s: needs review labels Aug 1, 2015

vbraun removed the s: positive review label Aug 1, 2015

vbraun closed this as completed in 4f379dc Aug 1, 2015

Faster Poyhedron.graph() #18860

Faster Poyhedron.graph() #18860

Comments

nathanncohen mannequin commented Jul 7, 2015

nathanncohen mannequin commented Jul 7, 2015

sagetrac-git mannequin commented Jul 7, 2015

sagetrac-git mannequin commented Jul 7, 2015

dimpase commented Jul 7, 2015

nathanncohen mannequin commented Jul 7, 2015

dimpase commented Jul 8, 2015

nathanncohen mannequin commented Jul 8, 2015

dimpase commented Jul 8, 2015

nathanncohen mannequin commented Jul 8, 2015

nathanncohen mannequin commented Jul 8, 2015

sagetrac-git mannequin commented Jul 8, 2015

sagetrac-git mannequin commented Jul 8, 2015

nathanncohen mannequin commented Jul 8, 2015

This comment has been minimized.

sagetrac-git mannequin commented Jul 8, 2015

sagetrac-git mannequin commented Jul 8, 2015

vbraun commented Jul 8, 2015

nathanncohen mannequin commented Jul 8, 2015

vbraun commented Jul 8, 2015

nathanncohen mannequin commented Jul 8, 2015

dimpase commented Jul 8, 2015

nathanncohen mannequin commented Jul 9, 2015

dimpase commented Jul 9, 2015

nathanncohen mannequin commented Jul 9, 2015

nathanncohen mannequin commented Jul 9, 2015

dimpase commented Jul 9, 2015

nathanncohen mannequin commented Jul 9, 2015

dimpase commented Jul 9, 2015

nathanncohen mannequin commented Jul 9, 2015

dimpase commented Jul 31, 2015

nathanncohen mannequin commented Jul 31, 2015

dimpase commented Jul 31, 2015

dimpase commented Jul 31, 2015

nathanncohen mannequin commented Jul 31, 2015

vbraun commented Aug 1, 2015

sagetrac-git mannequin commented Aug 1, 2015

sagetrac-git mannequin commented Aug 1, 2015

dimpase commented Aug 1, 2015

nathanncohen mannequin commented Aug 1, 2015

vbraun commented Aug 1, 2015

dimpase commented Aug 1, 2015

vbraun commented Aug 1, 2015