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: sage-6.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 pernici)

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 sage-6.2; it takes 0.08 seconds with this patch, using the Butera-Pernici 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 sage-6.2

sage: from sage.matrix.matrix_misc import permanental_minor_polynomial
sage: n, w = 100, 5
sage: m = matrix([[i*j + 1 if abs(i-j) <= 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)

grid.py (985 bytes) - added by pernici 7 years ago.

Download all attachments as: .zip

Change History (40)

comment:1 Changed 7 years ago by pernici

  • Authors set to Mario Pernici
  • Cc ncohen jsp added

comment:2 Changed 7 years ago by pernici

  • 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 git

  • Commit set to b0d714572f15dd915b1caabbb11f5f96f4020c85

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

b0d7145fixed bug in `permanental_minor_vector` in case of vanishing permanent

comment:4 Changed 7 years ago by pernici

  • Description modified (diff)

Changed 7 years ago by pernici

comment:5 Changed 7 years ago by pernici

  • Description modified (diff)

comment:6 Changed 7 years ago by pernici

  • Description modified (diff)
  • Summary changed from permanental_minor_vector to permanental_minor_vector, matching polynomial

comment:7 Changed 6 years ago by vbraun_spam

  • Milestone changed from sage-6.3 to sage-6.4

comment:8 Changed 6 years ago by chapoton

  • Branch changed from u/pernici/ticket/16603 to public/ticket/16603
  • Commit changed from b0d714572f15dd915b1caabbb11f5f96f4020c85 to f65acdedf20554ed8cbd88fde232da4229d96f85

I have had a quick look and made a review patch. You need to document the auxiliary function.


New commits:

55d328fMerge branch 'u/pernici/ticket/16603' in 6.4.rc0
f65acdetrac #16603 first review commit

comment:9 Changed 6 years ago by pernici

  • 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 pernici

  • 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:

947ddd1Added missing documentation in ``_prm_mul``;

comment:11 follow-up: Changed 6 years ago by chapoton

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 pernici

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 chapoton

  • 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 vdelecroix

  • 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 of permanental_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 either make doc or sage -docbuild reference html) then you will see the result in the section Matrices and Spaces of Matrices -> Miscellaneous matrix functions. I changed the name _prm_mul to prm_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 follow-up: Changed 6 years ago by pernici

Hello,

how can I access your branch u/vdelecroix/16603 ?

comment:16 in reply to: ↑ 15 Changed 6 years ago by vdelecroix

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 follow-up: Changed 6 years ago by pernici

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 vdelecroix

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 vdelecroix

  • Description modified (diff)

I rewrote the description (formatting the code).

comment:20 Changed 6 years ago by git

  • Commit changed from 947ddd16da731019f25a951e43bf551ad997a700 to f511623cce7bde7f022fa212a505526b3b62f0db

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

a6fd53btrac #16603: speedup + doc
f611f58trac #16603: add matrix_misc.py to the doc
f511623Added documentation to `permanental_minor_vector`

comment:21 Changed 6 years ago by vdelecroix

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 in matrix.rook_vector
  • I added a keyword var (as it is the case for matrix.charpoly for example) and change the argument name from m to A (since A is used in the docstring)
  • I optimized a bit the method prm_mul (it was possible to remove the loop computing exp) and introduced the list vars_to_do inplace of the set done_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 terms t^l^ with l > k)

Methods of matrices:

  • does your algorithm have a name? I used Butera-Pernici (see the item just below, it has some importance).
  • To allow the computation using your algorith, I made changes in matrix.permanental_minor and matrix.permanent. After my commits you can do
    sage: 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 vdelecroix

  • Status changed from needs_work to needs_info

changes pushed on u/vdelecroix/16603.

comment:23 Changed 6 years ago by git

  • Commit changed from f511623cce7bde7f022fa212a505526b3b62f0db to 19350a4574f4caa2e0b5ad08f66e51930b647a42

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

97f2062trac #16603: documentation + argument modifs
63e37cbtrac #16603: various optimizations
8a228d8trac #16603: modify methods of matrices
c76177eAdded `prec` parameter to `permanental_minor_polynomial`;
6f2c1c7used default "ButeraPernici" algorithm in `rook_vector`
19350a4Added "Godsil" algorithm to `rook_vector`.

comment:24 Changed 6 years ago by pernici

Hello Vincent,

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 terms tl with l > k)

I adapted the algorithm to compute fastly k-th 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 <= n-2

does your algorithm have a name? I used Butera-Pernici (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,n-1; 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 pernici

  • Description modified (diff)

comment:26 Changed 6 years ago by pernici

  • Description modified (diff)

comment:27 Changed 6 years ago by pernici

  • Description modified (diff)

comment:28 Changed 6 years ago by pernici

  • Status changed from needs_info to needs_review

comment:29 Changed 6 years ago by vdelecroix

  • 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:

3e7297btrac #16603: documentation commit

comment:30 Changed 6 years ago by pernici

  • 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 pernici

  • 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:

dd4636afixed bug in ``permanental_minor_poly``;
a60cc74added argument `complement`, eliminated argument `check` in ``rook_vector``;
e4359bbadded ``matching_polynomial`` method to ``BipartiteGraph`` class

comment:32 Changed 6 years ago by vdelecroix

  • 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:

e03f93atrac #16603: change keywords for `rook_vector` + doc

comment:33 Changed 6 years ago by pernici

  • 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 pernici

  • 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:

127418ffixed a bug in ``matching_polynomial``; added test.

comment:35 Changed 6 years ago by vdelecroix

  • 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 vbraun

  • Status changed from positive_review to needs_work

Doctest failure in src/sage/misc/sagedoc.py

comment:37 Changed 6 years ago by vdelecroix

  • 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

Funny! Done...


New commits:

901d36cfix doctest in misc/sagedoc.py

comment:38 Changed 6 years ago by vdelecroix

  • Status changed from needs_review to positive_review

comment:39 Changed 6 years ago by vbraun

  • Branch changed from u/vdelecroix/16603 to 901d36c702cfab6c13e1d14083cef7e8ee931c37
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.