Opened 7 years ago
Closed 6 years ago
#16603 closed enhancement (fixed)
permanental_minor_vector, matching polynomial
Reported by:  pernici  Owned by:  

Priority:  major  Milestone:  sage6.4 
Component:  linear algebra  Keywords:  
Cc:  ncohen, jsp  Merged in:  
Authors:  Mario Pernici  Reviewers:  Vincent Delecroix 
Report Upstream:  N/A  Work issues:  
Branch:  901d36c (Commits)  Commit:  901d36c702cfab6c13e1d14083cef7e8ee931c37 
Dependencies:  Stopgaps: 
Description (last modified by )
Added function permanental_minor_polynomial
, which computes the polynomial of the sums of the permanental minors in an efficient way, see arXiv:1406.5337.
As a benchmark, consider dancing.sage. On my computer dance(10)
took 3 hours using sage6.2; it takes 0.08 seconds with this patch, using the
ButeraPernici algorithm, 13 seconds using the Godsil algorithm.
In the case of banded matrices, the speedup is greater, since for fixed bandwidth the new algorithm is polynomial; here is a permanent which cannot be computed with sage6.2
sage: from sage.matrix.matrix_misc import permanental_minor_polynomial sage: n, w = 100, 5 sage: m = matrix([[i*j + 1 if abs(ij) <= w else 0 ....: for i in range(n)] for j in range(n)]) sage: time p = permanental_minor_polynomial(m, permanent_only=True) Wall time: 773 ms
This can be used to compute the matching polynomial of bipartite graphs much faster than matching_polynomial for certain graphs. For example, using the attached file to compute the reduced adjacency matrix of a 2d grid (n1, n2):
sage: from grid import sq_mat sage: from sage.matrix.matrix_misc import permanental_minor_polynomial sage: m = matrix(sq_mat(6, 6)) sage: time a = permanental_minor_polynomial(m) Wall time: 10.6 ms sage: g = BipartiteGraph(m) sage: time p = g.matching_polynomial() Wall time: 46.9 s
For n1,n2=10,10 permanental_minor_vector(m) takes 0.3s For n1,n2=50,6 it takes 0.5s
Attachments (1)
Change History (40)
comment:1 Changed 7 years ago by
 Cc ncohen jsp added
comment:2 Changed 7 years ago by
 Branch set to u/pernici/ticket/16603
 Created changed from 07/01/14 18:24:37 to 07/01/14 18:24:37
 Modified changed from 07/01/14 18:29:25 to 07/01/14 18:29:25
comment:3 Changed 7 years ago by
 Commit set to b0d714572f15dd915b1caabbb11f5f96f4020c85
comment:4 Changed 7 years ago by
 Description modified (diff)
Changed 7 years ago by
comment:5 Changed 7 years ago by
 Description modified (diff)
comment:6 Changed 7 years ago by
 Description modified (diff)
 Summary changed from permanental_minor_vector to permanental_minor_vector, matching polynomial
comment:7 Changed 6 years ago by
 Milestone changed from sage6.3 to sage6.4
comment:8 Changed 6 years ago by
 Branch changed from u/pernici/ticket/16603 to public/ticket/16603
 Commit changed from b0d714572f15dd915b1caabbb11f5f96f4020c85 to f65acdedf20554ed8cbd88fde232da4229d96f85
comment:9 Changed 6 years ago by
 Branch changed from public/ticket/16603 to u/pernici/ticket/16603
 Modified changed from 10/31/14 19:35:58 to 10/31/14 19:35:58
comment:10 Changed 6 years ago by
 Commit changed from f65acdedf20554ed8cbd88fde232da4229d96f85 to 947ddd16da731019f25a951e43bf551ad997a700
I made a commit in which I documented the auxiliary function, but by mistake I committed it to u/pernici/ticket/16603 instead of public/ticket/16603
New commits:
947ddd1  Added missing documentation in ``_prm_mul``;

comment:11 followup: ↓ 12 Changed 6 years ago by
Do you want this to be reviewed ? If yes, please turn it into "needs review".
comment:12 in reply to: ↑ 11 Changed 6 years ago by
Replying to chapoton:
Do you want this to be reviewed ? If yes, please turn it into "needs review".
How do I turn it into "needs review" ?
comment:13 Changed 6 years ago by
 Status changed from new to needs_review
Click on "Modify ticket", then on the button "needs review". I did it for you.
comment:14 Changed 6 years ago by
 Status changed from needs_review to needs_work
Hello,
This code is very cool! You can have a look at my branch u/vdelecroix/16603
where I provide two commits that I describe below.
 I edited a bit the code, I hope it is cleaner (and actually also a bit faster)
 Please, could you fill the section
ALGORITHM
ofpermanental_minor_vector
. It would be nice to have here a rough description of the algorithm.
 I edited a lot the documentation. First of all the file
matrix_misc.py
was not part of the matrix documentation! After my commits it is. If you compile the documentation (with eithermake doc
orsage docbuild reference html
) then you will see the result in the sectionMatrices and Spaces of Matrices
>Miscellaneous matrix functions
. I changed the name_prm_mul
toprm_mul
in order to make it appear in the documentation.
 Currently the situtation for benchmarks is unfair:
 the code of Jaap Spies to compute the permanent and the permanent minors can be optimized a lot!
 your code is in Python while Jaap's one is in Cython
 As you mentioned this algorithm can be used in several contexts. It would be nice to have some methods (of graphs and matrices) that call it. But if you prefer, we can do it in a further ticket.
Vincent
comment:15 followup: ↓ 16 Changed 6 years ago by
Hello,
how can I access your branch u/vdelecroix/16603 ?
comment:16 in reply to: ↑ 15 Changed 6 years ago by
Replying to pernici:
Hello,
how can I access your branch u/vdelecroix/16603 ?
How do you do with your branch u/pernici/ticket/16603
? All users have write access to public/*
and u/username/*
and read access everywhere. Just do in the Sage root:
$ git fetch trac u/vdelecroix/16603 $ git checkout b 16603_vincent_branch FETCH_HEAD
The first command download my two commits. The special temporary branch FETCH_HEAD
is a pointer to the last commit. The second command creates a new branch called 16603_vincent_branch
and move to it.
You can then merge these commits (if you like them) in your branch with
$ git checkout my_original_branch $ git merge 16603_vincent_branch
Is that clear enough?
Vincent
comment:17 followup: ↓ 18 Changed 6 years ago by
Previously I used ./sage dev checkout ticket 16603 Following your instructions I got your branch; thanks.
comment:18 in reply to: ↑ 17 Changed 6 years ago by
Replying to pernici:
Previously I used ./sage dev checkout ticket 16603 Following your instructions I got your branch; thanks.
I see. The sage dev
commands are a bit outdated (and give you much less freedom). You might better use git
or the enhanced git trac
that Volker Braun wrote. Have a look at the documentation.
comment:19 Changed 6 years ago by
 Description modified (diff)
I rewrote the description (formatting the code).
comment:20 Changed 6 years ago by
 Commit changed from 947ddd16da731019f25a951e43bf551ad997a700 to f511623cce7bde7f022fa212a505526b3b62f0db
comment:21 Changed 6 years ago by
Hello,
I will push several commits in a minute on my branch u/vdelecroix/16603
. Again you are free to use/discuss/reject them. I discuss the most important changes below.
Change in your original code and documentation:
 I changed your initial function into
permanental_minor_polynomial
. It looks more natural as it is what you compute... I moved the part of the code that makes it a list inmatrix.rook_vector
 I added a keyword
var
(as it is the case formatrix.charpoly
for example) and change the argument name fromm
toA
(since A is used in the docstring)
 I optimized a bit the method
prm_mul
(it was possible to remove the loop computingexp
) and introduced the listvars_to_do
inplace of the setdone_vars
.
 Sometimes it might happen that some coefficient of the partial product become zero. I added a test in
prm_mul
. If you think this is not worth it, I can revert that change.
 Again, I edited a lot the documentation of
permanental_minor_polynomial
. Please reread carefully, I tried to make it more understandable but I might have made many mistakes.
 Do you know an adaptation of your method to compute only a specific permanental minor? (I guess that to compute the
k
th permanental minor it is just enough to ignore the termst^l^
withl > k
)
Methods of matrices:
 does your algorithm have a name? I used ButeraPernici (see the item just below, it has some importance).
 To allow the computation using your algorith, I made changes in
matrix.permanental_minor
andmatrix.permanent
. After my commits you can dosage: m = matrix(ZZ,4,4,range(16)) sage: m.permanent(algorithm="Ryser") 22432 sage: m.permanent(algorithm="ButeraPernici") 22432
it will be easier to make benchmarks (but I repeat that right now it is not fair for both, see also #17457).
 Right now, in all methods except
matrix.rook_vector
, the default choice for algorithm is Ryser (what was before). Where do you think we should put yours in first?
 why
matrix.rook_vector
is only defined for {0,1} matrices... it is well defined for any matrix, isn't it?
 possibly we can add a method
matrix.permanental_minor_polynomial
but it will overlap a lot with rook vector. I am not sure what to do as these methods are very specific to matrix combinatorics. Perhaps you have a better overview than me on that.
Despite this big list of remarks, I am sure this is not far for being ready for integration! This is a really nice piece of work.
Vincent
comment:22 Changed 6 years ago by
 Status changed from needs_work to needs_info
changes pushed on u/vdelecroix/16603
.
comment:23 Changed 6 years ago by
 Commit changed from f511623cce7bde7f022fa212a505526b3b62f0db to 19350a4574f4caa2e0b5ad08f66e51930b647a42
Branch pushed to git repo; I updated commit sha1. New commits:
97f2062  trac #16603: documentation + argument modifs

63e37cb  trac #16603: various optimizations

8a228d8  trac #16603: modify methods of matrices

c76177e  Added `prec` parameter to `permanental_minor_polynomial`;

6f2c1c7  used default "ButeraPernici" algorithm in `rook_vector`

19350a4  Added "Godsil" algorithm to `rook_vector`.

comment:24 Changed 6 years ago by
Hello Vincent,
Do you know an adaptation of your method to compute only a specific permanental minor? (I guess that to compute the kth permanental minor it is just enough to ignore the terms t^{l} with l > k)
I adapted the algorithm to compute fastly kth permanental minor for small k
;
testing on A = matrix(n, k, range(n)
using this algorithm
A.permanent_minor(k)
is faster for `2 <= k <= n2
does your algorithm have a name? I used ButeraPernici (see the item just below, it has some importance).
No, it did not have a name.
Right now, in all methods except matrix.rook_vector, the default choice for algorithm is Ryser (what was before). Where do you think we should put yours in first?
In the example above, for n x n matrices in permanental_minor
the Ryser algorithm is generally faster for k=n
(permanent), slightly faster for
k=1,n1
; maybe one could make a default option which uses the Ryser
algorithm in those cases, the new algorithm in the other cases.
why matrix.rook_vector is only defined for {0,1} matrices... it is well defined for any matrix, isn't it?
Usually rook vector
is defined only on {0,1} matrices; I added the "Godsil" algorithm which
is faster than the "Ryser" algorithm.
comment:25 Changed 6 years ago by
 Description modified (diff)
comment:26 Changed 6 years ago by
 Description modified (diff)
comment:27 Changed 6 years ago by
 Description modified (diff)
comment:28 Changed 6 years ago by
 Status changed from needs_info to needs_review
comment:29 Changed 6 years ago by
 Branch changed from u/pernici/ticket/16603 to u/vdelecroix/16603
 Commit changed from 19350a4574f4caa2e0b5ad08f66e51930b647a42 to 3e7297b35f94986b8e0f37f30a3b3d6ee3d3b699
 Reviewers set to Vincent Delecroix
Hello,
I really think that the documentation is half of a software: it provides specification of functions and help the user to orient himself/herself. Next time you implement a function in Sage pay attention that your new functionality is referenced from other part of the code and that it also make references to relevant functions. There is a lot of possibilities for organization of the documentation (see the section "Docstring" of the developer guide).
I was not happy with the current behavior and documentation of rook_vector
and permanental_minor
(these methods should work for rectangular matrices whatever their shape). I guess that this is my last commit for that ticket. Just have a look and if you agree we can set the ticket to positive review.
Next step would be optimisation of Ryser algorithm: #17457.
Vincent
New commits:
3e7297b  trac #16603: documentation commit

comment:30 Changed 6 years ago by
 Branch changed from u/vdelecroix/16603 to u/pernici/ticket/16603
 Modified changed from 12/18/14 23:23:23 to 12/18/14 23:23:23
comment:31 Changed 6 years ago by
 Commit changed from 3e7297b35f94986b8e0f37f30a3b3d6ee3d3b699 to e4359bb3d2444375eed077c72274e4f6a4ed87bc
Hello,
I added the matching_polynomial
method to the
BipartiteGraph?
class, so one can optionally
use the "rook" algorithm.
Now
rook_vector
is much faster when the matrix is made mostly of ones.
I fixed a couple of bugs and added documentation.
A simple optimization I did not include here is the one in case of matrices which can be put in block form permuting rows and columns. I did not include it because it does not seem to be a commonly occurring case. The same optimization can be used in the computation of determinants of large matrices of this kind. Do you think that it is worth opening a ticket on it, or is it too obvious and useless?
Sooner or later I will open a ticket for a faster algorithm for the matching polynomial of graphs (non bipartite).
New commits:
dd4636a  fixed bug in ``permanental_minor_poly``;

a60cc74  added argument `complement`, eliminated argument `check` in ``rook_vector``;

e4359bb  added ``matching_polynomial`` method to ``BipartiteGraph`` class

comment:32 Changed 6 years ago by
 Branch changed from u/pernici/ticket/16603 to u/vdelecroix/16603
 Commit changed from e4359bb3d2444375eed077c72274e4f6a4ed87bc to e03f93ae3016870383b50bf2a3abea4b5234e51b
Hello,
First of all, there is already a lot of material on this ticket. It would be good to stop here and open tickets for new features (please put me on cc on them).
I push a commit as I do not like the complement
keyword that sounds like you want the rook vector of the complement matrix. I replaced it with use_complement
. I know that it is different from the specification for BipartiteGraphs
but it fits with IndependentSets
sage: from sage.graphs.independent_sets import IndependentSets sage: I = IndependentSets(graphs.GridGraph((5,5))) sage: I.cardinality() 55447 sage: I = IndependentSets(graphs.GridGraph((5,5)), complement=True) sage: I.cardinality() 66
Nevertheless, I kept the keyword complement
to do the following as it is very useful
sage: identity_matrix(6).rook_vector(complement=True) [1, 30, 315, 1420, 2715, 1854, 265]
Vincent
New commits:
e03f93a  trac #16603: change keywords for `rook_vector` + doc

comment:33 Changed 6 years ago by
 Branch changed from u/vdelecroix/16603 to u/pernici/ticket/16603
 Modified changed from 01/09/15 09:22:05 to 01/09/15 09:22:05
comment:34 Changed 6 years ago by
 Commit changed from e03f93ae3016870383b50bf2a3abea4b5234e51b to 127418fc436cad3561c298a8f5c685c1e4b4fb81
Hello, I made one more small change. Thank you for reviewing and improving the code and documentation.
New commits:
127418f  fixed a bug in ``matching_polynomial``; added test.

comment:35 Changed 6 years ago by
 Status changed from needs_review to positive_review
Hello,
Let us stop there with positive review! If you want to go any further open a new ticket (and please put me in cc).
Vincent
comment:36 Changed 6 years ago by
 Status changed from positive_review to needs_work
Doctest failure in src/sage/misc/sagedoc.py
comment:37 Changed 6 years ago by
 Branch changed from u/pernici/ticket/16603 to u/vdelecroix/16603
 Commit changed from 127418fc436cad3561c298a8f5c685c1e4b4fb81 to 901d36c702cfab6c13e1d14083cef7e8ee931c37
 Status changed from needs_work to needs_review
comment:38 Changed 6 years ago by
 Status changed from needs_review to positive_review
comment:39 Changed 6 years ago by
 Branch changed from u/vdelecroix/16603 to 901d36c702cfab6c13e1d14083cef7e8ee931c37
 Resolution set to fixed
 Status changed from positive_review to closed
Branch pushed to git repo; I updated commit sha1. New commits:
fixed bug in `permanental_minor_vector` in case of vanishing permanent