Improvement for incidence matrix of polyhedra

[gh-kliem]

This ticket collects some improvements for obtaining incidence matrices:

• The method self._is_zero is called many times and we can save plenty of time by obtaining it beforehand.
• We can obtain the slack matrix by multiplying the Vrepresentation matrix with the homogenous Vrepresentation. From this we can obtain the incidence matrix (strangely this last step takes up most of the time).

Unfortunately, matrix multiplication for algebraic fields isn't very efficient. So we need to keep the old method in this case.

Before this ticket:

sage: P = polytopes.permutahedron(7)
sage: %time _ = P.incidence_matrix()
CPU times: user 751 ms, sys: 7.95 ms, total: 759 ms
Wall time: 758 ms
sage: P = polytopes.associahedron(['A', 9])
sage: %time _ = P.incidence_matrix()
CPU times: user 1.94 s, sys: 4.03 ms, total: 1.94 s
Wall time: 1.94 s
sage: P = polytopes.hypercube(14)
sage: %time _ = P.incidence_matrix()
CPU times: user 670 ms, sys: 0 ns, total: 670 ms
Wall time: 669 ms
sage: K.<sqrt2> = QuadraticField(2)
sage: P = polytopes.permutahedron(7).dilation(sqrt2)
sage: %time _ = P.incidence_matrix()
CPU times: user 3.82 s, sys: 0 ns, total: 3.82 s
Wall time: 3.82 s

With this ticket:

sage: P = polytopes.permutahedron(7)
sage: %time _ = P.incidence_matrix()
CPU times: user 64.4 ms, sys: 7 µs, total: 64.4 ms
Wall time: 63.6 ms
sage: P = polytopes.associahedron(['A', 9])
sage: %time _ = P.incidence_matrix()
CPU times: user 325 ms, sys: 12.1 ms, total: 337 ms
Wall time: 335 ms
sage: P = polytopes.hypercube(14)
sage: %time _ = P.incidence_matrix()
CPU times: user 139 ms, sys: 12 ms, total: 151 ms
Wall time: 149 ms
sage: K.<sqrt2> = QuadraticField(2)
sage: P = polytopes.permutahedron(7).dilation(sqrt2)
sage: %time _ = P.incidence_matrix()
CPU times: user 3.69 s, sys: 3.69 ms, total: 3.69 s
Wall time: 3.69 s

[gh-kliem]

New commits:

​42b5c42	determine is_zero beforehand
​02a8a2f	small fix
​6dbe5d2	use matrix multiplication to obtain incidence matrix

[gh-kliem]

This should be rebased on #29838 and another ticket that speeds up getting the incidence matrix from the slack matrix.

[gh-kliem]

New commits:

​43e6dfb	implement _zero_matrix
​05eba1e	fix bug
​7b70830	rename
​47a19c6	implement slack matrix
​54755d7	Merge branch 'public/29838' of git://trac.sagemath.org/sage into public/29837-reb
​a4576d6	set up incidence/adjacency matrix with GF(2)
​30c4715	Merge branch 'public/29840' of git://trac.sagemath.org/sage into public/29837-reb
​fcd62d2	determine is_zero beforehand
​13b8583	small fix
​89aeab3	obtain incidence matrix from slack matrix

[mkoeppe]

I don't think changing the interface so that GF(2) matrices are returned is a good idea. How about leaving it as ZZ by default, and using a keyword argument for the base ring so it can be overridden.

[gh-kliem]

Replying to mkoeppe:

I don't think changing the interface so that GF(2) matrices are returned is a good idea. How about leaving it as ZZ by default, and using a keyword argument for the base ring so it can be overridden.

I don't think the change to GF(2) was needed to make things better. So I can rewrite some stuff and remove #29839, that would also do it, I think. sage.matrix.matrix_integer_dense is probably just as fast or faster than the other one.

I really only thought that this is the more natural thing for a 0/1 matrix. If the more natural thing is base ring integer or GF(2) really performs awful, I can just revert that decision.

[mkoeppe]

See my comment in #29840

[gh-kliem]

Changing the base ring of incidence matrices was rejected in #29340.

[dimpase]

I wish this was implemented using numpy matrices (knows as arrays), which provide faster linear algebra and have built-in methods to manipulate zero patterns of their matrices.

[gh-kliem]

I don't know how fast conversion back and forth works. It is true, that numpy is a bit faster, but that doesn't help me if I cannot convert things.

E.g.

sage: P = polytopes.associahedron(['A', 9])
sage: %time a = np.array(P.Vrepresentation())
CPU times: user 383 ms, sys: 0 ns, total: 383 ms
Wall time: 382 ms

This is already longer than it takes to obtain the entire incidence matrix. Of course this is most likely not an ideal way to convert, but I don't know any better way at the moment.

And conversion back is also slow.

Replying to dimpase:

I wish this was implemented using numpy matrices (knows as arrays), which provide faster linear algebra and have built-in methods to manipulate zero patterns of their matrices.

[gh-kliem]

New commits:

​62f1ff9	use integer matrix for zero_pattern_matrix
​88f6025	Merge branch 'public/29838' of git://trac.sagemath.org/sage into public/29837-reb2
​7cee020	determine is_zero beforehand
​3d28d09	small fix

[git]

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

​7ef6a11

use slack matrix to obtain incidence matrix quicker

[git]

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

​b583745	make methods more generic
​535dc85	note for developers
​0988730	Merge branch 'public/29839' of git://trac.sagemath.org/sage into public/29837-reb2
​759d1f2	account for change in #29839

[gh-kliem]

Last 10 new commits:

​9462376	Adding optimized get_is_zero_unsafe() to rational sparse matrices.
​94febde	Custom get_is_zero_unsafe() for modn sparse matrices.
​bf0352f	See if we can check trivial zeros first for symbolic matrices.
​2a08abe	small mistakes
​333c4ba	wrong doctest
​82c6b4a	docstring
​211f4f1	Merge branch 'public/29839' of git://trac.sagemath.org/sage into public/29837-reb3
​4c5433e	determine is_zero beforehand
​27a1b90	use slack matrix to obtain incidence matrix quicker
​764a99b	account for change in #29839

[gh-kliem]

Needs review.

The the ticket itself is only 3 commits on top of #29838 and #29839, which are both closed now.

[tscrim]

You might be able to get a little more speed out by doing

-return x == 0
+return (not x)

[tscrim]

However, whether you make that change or not, you can set a positive review on my behalf.

[gh-kliem]

About three percent for P = polytopes.permutahedron(7).dilation(sqrt2).

Thank you.

[git]

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

​f70e1c5

improvement suggested by reviewer

[fbissey]

I am getting this with this ticket in sage-on-gentoo

sage -t --long --warn-long 101.0 /usr/lib/python3.7/site-packages/sage/geometry/polyhedron/backend_cdd.py
**********************************************************************
File "/usr/lib/python3.7/site-packages/sage/geometry/polyhedron/backend_cdd.py", line 508, in sage.geometry.polyhedron.backend_cdd.Polyhedron_RDF_cdd.__init__
Failed example:
    TestSuite(p).run()
Expected nothing
Got:
    Failure in _test_an_affine_basis:
    Traceback (most recent call last):
      File "/usr/lib/python3.7/site-packages/sage/misc/sage_unittest.py", line 296, in run
        test_method(tester=tester)
      File "/usr/lib/python3.7/site-packages/sage/geometry/polyhedron/base.py", line 2122, in _test_an_affine_basis
        b = self.an_affine_basis()
      File "/usr/lib/python3.7/site-packages/sage/geometry/polyhedron/base.py", line 2085, in an_affine_basis
        chain = self.a_maximal_chain()[1:]  # we exclude the empty face
      File "/usr/lib/python3.7/site-packages/sage/geometry/polyhedron/base.py", line 3098, in a_maximal_chain
        comb_chain = self.combinatorial_polyhedron().a_maximal_chain()
      File "sage/misc/cachefunc.pyx", line 2310, in sage.misc.cachefunc.CachedMethodCallerNoArgs.__call__ (/dev/shm/portage/sci-mathematics/sage-9999/work/sage-9999/src-python3_7/build/cythonized/sage/misc/cachefunc.c:12875)
        self.cache = f(self._instance)
      File "/usr/lib/python3.7/site-packages/sage/geometry/polyhedron/base.py", line 3362, in combinatorial_polyhedron
        return CombinatorialPolyhedron(self)
      File "sage/geometry/polyhedron/combinatorial_polyhedron/base.pyx", line 367, in sage.geometry.polyhedron.combinatorial_polyhedron.base.CombinatorialPolyhedron.__init__ (/dev/shm/portage/sci-mathematics/sage-9999/work/sage-9999/src-python3_7/build/cythonized/sage/geometry/polyhedron/combinatorial_polyhedron/base.c:6786)
        data = data.incidence_matrix()
      File "sage/misc/cachefunc.pyx", line 2310, in sage.misc.cachefunc.CachedMethodCallerNoArgs.__call__ (/dev/shm/portage/sci-mathematics/sage-9999/work/sage-9999/src-python3_7/build/cythonized/sage/misc/cachefunc.c:12875)
        self.cache = f(self._instance)
      File "/usr/lib/python3.7/site-packages/sage/geometry/polyhedron/base.py", line 2749, in incidence_matrix
        is_zero == self._is_zero
    UnboundLocalError: local variable 'is_zero' referenced before assignment
    ------------------------------------------------------------
    The following tests failed: _test_an_affine_basis
**********************************************************************

seems circular.

[tscrim]

-is_zero == self._is_zero
+is_zero = self._is_zero

[gh-kliem]

The problem is that usually the cdd backend figures out the incidence matrix itself for inexact polyhedra. So you don't notice, because incidence_matrix is never run (except for the empty polyhedron, hence the error). I should add a doctest for this.

[gh-kliem]

I rebased as well.

---

New commits:

​169077a	Merge branch 'public/29837-reb3' of git://trac.sagemath.org/sage into public/29837-reb4
​26590f1	fix mistake
​da9ec61	add doctest

[tscrim]

Trivial tweak, but the formatting is better:

         Test that this method works for inexact base ring
-        (`cdd` sets the cache already)::
+        (``cdd`` sets the cache already)::

[git]

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

​955f12b

sphinx syntax

[gh-kliem]

Done.

Replying to tscrim:

Trivial tweak, but the formatting is better:

         Test that this method works for inexact base ring
-        (`cdd` sets the cache already)::
+        (``cdd`` sets the cache already)::

[tscrim]

Thanks.