Opened 4 years ago

Closed 3 years ago

#18860 closed enhancement (fixed)

Faster Poyhedron.graph()

Reported by: ncohen Owned by:
Priority: major Milestone: sage-6.8
Component: geometry Keywords:
Cc: vdelecroix, dimpase, vbraun, novoselt Merged in:
Authors: Nathann Cohen Reviewers: Dima Pasechnik
Report Upstream: N/A Work issues:
Branch: 5223cda (Commits) Commit: 5223cda9f15b2a417b550039b978a07e126f9053
Dependencies: #18779 Stopgaps:

Description (last modified by ncohen)

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

Change History (39)

comment:1 Changed 4 years ago by ncohen

  • Branch set to public/18860

comment:2 Changed 4 years ago by git

  • Commit set to 3b5eabb9735a867678c2c6a13dc7700c24d82732

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

4a4ee0ctrac #18779: polytopes.gosset_3_21 and graphs.GossetGraph
b498a9etrac #18779: Reviewer's (dead right) comments
6527f4bcapitalise Gosset
3b5eabbtrac #18860: Faster Poyhedron.graph()

comment:3 Changed 4 years ago by ncohen

  • Cc vbraun added

comment:4 follow-up: Changed 4 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: Changed 4 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: Changed 4 years ago by dimpase

  • Cc novoselt added

Replying to 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?

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: Changed 4 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 4 years ago by dimpase

Replying to 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

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 4 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 4 years ago by ncohen (previous) (diff)

comment:10 follow-up: Changed 4 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 4 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:

d084b1ftrac #18860: Faster Poyhedron.graph()

comment:12 follow-up: Changed 4 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 4 years ago by ncohen (previous) (diff)

comment:13 Changed 4 years ago by ncohen

  • Status changed from new to needs_review

comment:14 Changed 4 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:

74d8412trac #18860: Faster Poyhedron.graph()

comment:15 follow-up: Changed 4 years ago by vbraun

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

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

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

comment:17 follow-up: Changed 4 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 4 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: Changed 4 years ago by dimpase

Replying to 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.

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 4 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: Changed 4 years ago by dimpase

Replying to ncohen:

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 4 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: Changed 4 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: Changed 4 years ago by dimpase

Replying to 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.

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: Changed 4 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 4 years ago by dimpase

Replying to ncohen:

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 4 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: Changed 3 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 3 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 3 years ago by dimpase

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

OK. Merges in 6.8, no problem...

comment:31 Changed 3 years ago by ncohen

Thanks!

comment:32 Changed 3 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 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...

comment:33 Changed 3 years ago by git

  • Commit changed from 74d8412010748142918edf6178636bffdb14e60a to 5223cda9f15b2a417b550039b978a07e126f9053

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

5dfc776trac #18860: Merged with 6.9.beta0
5223cdatrac #18860: Remove vertex labels

comment:34 follow-up: Changed 3 years ago by ncohen

  • Status changed from needs_work to needs_review

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

Replying to ncohen:

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

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

comment:36 Changed 3 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 3 years ago by vbraun

Vertices are definitely reordered in Polyhedron.

comment:38 Changed 3 years ago by dimpase

  • Status changed from needs_review to positive_review

OK, looks good.

comment:39 Changed 3 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.