faster matching polynomial

[735e67df-bbe3-427b-86a6-1ad4def7a38c]

Added matching_generating_poly, modified matching_polynomial, added an option to permanental_minor_polynomial, which for some graphs is faster than the previous algorithm.

matching_generating_poly computes the matching generating polynomial on (possibly weighted) graphs; this is used in permanental_minor_polynomial.

The algorithm (see [ButPer]) uses polynomials on nilpotent elements, implemented in the class Hobj in matchpoly.pyx. The idea is that the matching generating polynomial can be computed from the product of terms (1 + x w_{i j} \eta_i \eta_j), for all edges (i, j) where w_{i j} is the weight of the edge (i, j) and \eta_i^2 = \eta_j^2=0. The term 1 corresponds to the absence of the dimer (i,j), the other term to its presence. A configuration with two adjacent dimers does not contribute due to nilpotency. The matching generating polynomial is given by the sum of the coefficients of the non-vanishing products of \eta's.

While the result does not depend on the ordering of the edges in the above product, the speed of the algorithm depends on the ordering; in fact doing the above product from left to right, one can set \eta_i=1 as soon as \eta_i does not appear in the yet unused factors on the right, reducing in this way the number of terms. A greedy argument to order edges is contained in active_nodes.py

The Godsil algorithm is faster for small graphs, which take little time with either algorithms; it becomes progressively slower as the number of vertices increases; e.g. on my computer for graphs.KnightGraph([4, n]), the Godsil algorithm for n=5 is 25% faster, for n=6 5x slower, for n=9 4000x slower.

The new algorithm can be much faster than the rook algorithm for matching polynomials of bipartite graphs

sage: time BipartiteGraph(graphs.LadderGraph(30)).matching_polynomial()[0]
Wall time: 65.7 ms
1346269

It is also much faster than the previous algorithm (now called bipartite) for computing the sum of the permanental minors, in the case of randomized band matrices

sage: from sage.matrix.matrix_misc import permanental_minor_polynomial
sage: n, w = 20, 3
sage: m = matrix([[i*j + 1 if abs(i-j) <= w else 0 for i in range(n)] for j in range(n)])
sage: a = list(m); shuffle(a); b = zip(*a); shuffle(b); m1 = matrix(b)
sage: time p1 = permanental_minor_polynomial(m1, algorithm='matching')
Wall time: 174 ms

CC: @videlec @nathanncohen

Component: graph theory

Branch/Commit: u/pernici/ticket/17921 @ c87c039

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

[735e67df-bbe3-427b-86a6-1ad4def7a38c]

Branch: u/pernici/ticket/17921

[735e67df-bbe3-427b-86a6-1ad4def7a38c]

Commit: fe2fa06

[735e67df-bbe3-427b-86a6-1ad4def7a38c]

Description changed:

--- 
+++ 
@@ -1,12 +1,24 @@
 ADD DESCRIPTION

-Added ``matching_generating_poly``, modified ``matching_polynomial``
+Added ``matching_generating_poly``, modified ``matching_polynomial``, added an option
+to ``permanental_minor_polynomial``, which for some graphs is faster than the previous
+algorithm.

 ``matching_generating_poly`` computes the matching generating polynomial
-on (possibly weighted) graphs.
+on (possibly weighted) graphs; this is used in ``permanental_minor_polynomial``.

 The algorithm (see [ButPer]) uses polynomials on nilpotent elements,
 implemented in the class ``Hobj`` in matchpoly.pyx,
 and a greedy argument to order edges, contained in active_nodes.py

-This algorithm is faster than the Godsil algorithm on large graphs.
+This algorithm is much faster than the Godsil algorithm on large graphs.
+
+``
+sage: from sage.graphs.matchpoly import matching_generating_poly
+sage: g = graphs.BuckyBall()
+sage: time sum(matching_generating_poly(g).coefficients())
+CPU times: user 184 ms, sys: 0 ns, total: 184 ms
+Wall time: 184 ms
+1417036634543488
+``
+

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

2961204

small change in ``BipartiteGraph.matching_polynomial``

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from fe2fa06 to 2961204

[videlec]

comment:6

I clean the description of the ticket to make it readable. Could you try to make it understandable?

[videlec]

Description changed:

--- 
+++ 
@@ -1,24 +1,21 @@
-ADD DESCRIPTION
-
-Added ``matching_generating_poly``, modified ``matching_polynomial``, added an option
-to ``permanental_minor_polynomial``, which for some graphs is faster than the previous
+Added `matching_generating_poly`, modified `matching_polynomial`, added an option
+to `permanental_minor_polynomial`, which for some graphs is faster than the previous
 algorithm.

-``matching_generating_poly`` computes the matching generating polynomial
-on (possibly weighted) graphs; this is used in ``permanental_minor_polynomial``.
+`matching_generating_poly` computes the matching generating polynomial
+on (possibly weighted) graphs; this is used in `permanental_minor_polynomial`.

 The algorithm (see [ButPer]) uses polynomials on nilpotent elements,
-implemented in the class ``Hobj`` in matchpoly.pyx,
+implemented in the class `Hobj` in matchpoly.pyx,
 and a greedy argument to order edges, contained in active_nodes.py

 This algorithm is much faster than the Godsil algorithm on large graphs.

-``
+```
 sage: from sage.graphs.matchpoly import matching_generating_poly
 sage: g = graphs.BuckyBall()
 sage: time sum(matching_generating_poly(g).coefficients())
 CPU times: user 184 ms, sys: 0 ns, total: 184 ms
 Wall time: 184 ms
 1417036634543488
-``
-
+```

[735e67df-bbe3-427b-86a6-1ad4def7a38c]

Description changed:

--- 
+++ 
@@ -6,8 +6,23 @@
 on (possibly weighted) graphs; this is used in `permanental_minor_polynomial`.

 The algorithm (see [ButPer]) uses polynomials on nilpotent elements,
-implemented in the class `Hobj` in matchpoly.pyx,
-and a greedy argument to order edges, contained in active_nodes.py
+implemented in the class `Hobj` in matchpoly.pyx.
+The idea is that the matching generating polynomial can be computed from
+the product of terms
+`(1 + x w_{i j} \eta_i \eta_j)`, for all edges `(i, j)`
+where `w_{i j}` is the weight of the edge `(i, j)` and `\eta_i^2 = \eta_j^2=0`.
+The term `1` corresponds to the absence of the dimer `(i,j)`,
+the other term to its presence. A configuration with two adjacent dimers
+does not contribute due to nilpotency. The matching generating polynomial
+is given by the sum of the coefficients of the non-vanishing products of `\eta`'s.
+
+While the result does not depend on the ordering of the edges in the above
+product, the speed of the algorithm depends on the ordering; in fact doing
+the above product from left to right, one can set `\eta_i=1` as soon
+as `\eta_i` does not appear in the yet unused factors on the right, reducing
+in this way the number of terms.
+A greedy argument to order edges is contained in active_nodes.py
+

 This algorithm is much faster than the Godsil algorithm on large graphs.

@@ -15,7 +30,27 @@
 sage: from sage.graphs.matchpoly import matching_generating_poly
 sage: g = graphs.BuckyBall()
 sage: time sum(matching_generating_poly(g).coefficients())
-CPU times: user 184 ms, sys: 0 ns, total: 184 ms
 Wall time: 184 ms
 1417036634543488

+It can be much faster than the rook algorithm for matching polynomials of bipartite graphs + + +sage: time BipartiteGraph(graphs.LadderGraph(30)).matching_polynomial()[0] +Wall time: 26.9 ms +1346269 + + +It is also much faster than the previous algorithm (now called bipartite) +for computing the sum of the permanental minors, in the case of randomized band matrices + + +sage: from sage.matrix.matrix_misc import permanental_minor_polynomial +sage: n, w = 20, 3 +sage: m = matrix([[i*j + 1 if abs(i-j) <= w else 0 for i in range(n)] for j in range(n)]) +sage: a = list(m); shuffle(a); b = zip(*a); shuffle(b); m1 = matrix(b) +sage: time p1 = permanental_minor_polynomial(m1, algorithm='matching') +Wall time: 174 ms + + +

[735e67df-bbe3-427b-86a6-1ad4def7a38c]

comment:9

I added an explanation on the use of nilpotent elements and two more examples.

[735e67df-bbe3-427b-86a6-1ad4def7a38c]

comment:10

The greedy argument to order edges in active_nodes.py orders edges in such a way that the graph G_i defined by edges[:i] has few nodes which have not the same degree as in G (active nodes). This leads to few \eta_i in the corresponding polynomial, so there are few terms. The greedy algorithm adds edges without increasing the number of active nodes, if possible; otherwise it adds a short path between two active nodes.

This algorithm can be probably improved using some graph algorithm already in Sage, maybe an algorithm to reorder vertices used to visualize graphs.

[rwst]

comment:11

Please try to get 100% coverage (sage -coverage) also in private methods. Please adhere to usual documentation guidelines: esp. start paragraphs with big caps, do not use >>> for examples.

[fchapoton]

comment:12

failing doctests, see patchbot report

plus very bad formatting of the doc

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

923dcc5	small change in ``_add_links_no_incr``; improved the documentation
5b10b9e	``ordered_links`` has now a single parameter;
c87c039	renamed `_obj_free`; used the shorter name "ButPer"

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 2961204 to c87c039

[735e67df-bbe3-427b-86a6-1ad4def7a38c]

Description changed:

--- 
+++ 
@@ -37,7 +37,7 @@

sage: time BipartiteGraph(graphs.LadderGraph(30)).matching_polynomial()[0] -Wall time: 26.9 ms +Wall time: 65.7 ms 1346269

[735e67df-bbe3-427b-86a6-1ad4def7a38c]

Description changed:

--- 
+++ 
@@ -24,16 +24,13 @@
 A greedy argument to order edges is contained in active_nodes.py

-This algorithm is much faster than the Godsil algorithm on large graphs.
+The Godsil algorithm is faster for small graphs, which take little time
+with either algorithms; it becomes progressively slower as the number
+of vertices increases; e.g. on my computer
+for `graphs.KnightGraph([4, n])`, the Godsil algorithm 
+for `n=5` is 25% faster, for `n=6` 5x slower, for `n=9` 4000x slower.

-```
-sage: from sage.graphs.matchpoly import matching_generating_poly
-sage: g = graphs.BuckyBall()
-sage: time sum(matching_generating_poly(g).coefficients())
-Wall time: 184 ms
-1417036634543488
-```
-It can be much faster than the `rook` algorithm for matching polynomials of bipartite graphs
+The new algorithm can be much faster than the `rook` algorithm for matching polynomials of bipartite graphs

sage: time BipartiteGraph(graphs.LadderGraph(30)).matching_polynomial()[0]

[jdemeyer]

comment:17

I personally don't care about this ticket, but let me mention that it doesn't apply and that it should be rebased on top of #23126.