maliasadi commented 3 years ago

Context: Currently, The Walrus uses the Ryser's formula to calculate permanents of arbitrary matrices that has time-complexity of O(n 2^n) if it's implemented with Gray code ordering. Another useful permanent algorithm with the same complexity is through BBFG's formula that's introduced by David G. Glynn in"The permanent of a square matrix". This pull request adds both C++ and Python support for permanent calculation using the latter formula.

Description of the Change:

Implemented a templated version of the BBFG formula in C++. It's a sequential algorithm (nthreads=1) from "The permanent of a square matrix" by David G. Glynn with using Gray code ordering. This is added in the ./include/permanent.hpp.
Added two wrappers for Python integration (perm_BBFG_qp and perm_BBFG_cmplx) and one for the future parallel vs. serial integration (perm_BBFG).
Added tests to ./tests/libwalrus_unittest.cpp
Integrated algorithms in ./thewalrus/libwalrus.pyx
Added tests to ./thewalrus/tests/test_permanent.py for integrated python methods naming perm_BBFG, perm_BBFG_real and perm_BBFG_complex.
...

Benefits: As @mlxd suggested, having both methods implemented will allow for performance and accuracy comparisons and trade-offs.

Possible Drawbacks:

Related GitHub Issues: Issue #263 created by @mlxd

maliasadi commented 3 years ago

A brief discussion on new algorithms in ./include/permanent.hpp:

My first implementation of BBFG algorithm is perm_BBFG_serial that is in fact my first understanding of this formula presented in "The permanent of a square matrix" with adding Gray code ordering support. In my second try to parallel the code, I realized that there are still some embarrassingly obvious optimizations there, so, I optimized my first implementation and name the new one perm_BBFG_serial2. I compare the performance of both BBFG algorithms with the Ryser's implementation (the permanent function) for n*n double random matrices with ntrials=10, nthreads=1 and g++ -O3:

[mat=20*20] BBFG_serial: 17ms
[mat=20*20] BBFG_serial2: 12ms
[mat=20*20] Ryser_permanent: 21ms

[mat=25*25] BBFG_serial: 544ms
[mat=25*25] BBFG_serial2: 442ms
[mat=25*25] Ryser_permanent: 818ms

[mat=30*30] BBFG_serial: 21182ms
[mat=30*30] BBFG_serial2: 18218ms
[mat=30*30] Ryser_permanent: 33242ms

[mat=31*31] BBFG_serial: 46229ms
[mat=31*31] BBFG_serial2: 38278ms
[mat=31*31] Ryser_permanent: 67662ms

[mat=32*32] BBFG_serial: 88446ms
[mat=32*32] BBFG_serial2: 74259ms
[mat=32*32] Ryser_permanent: 141167ms

Please, let me know what you think about further improvements.

maliasadi commented 3 years ago

Thank you @mlxd for the suggestions!

New benchmark comparing permanent algorithms:

[mat=20*20] BBFG_serial: 12ms
[mat=20*20] BBFG_serial2: 10ms
[mat=20*20] Ryser_permanent: 21ms

[mat=25*25] BBFG_serial: 471ms
[mat=25*25] BBFG_serial2: 420ms
[mat=25*25] Ryser_permanent: 821ms

[mat=30*30] BBFG_serial: 19510ms
[mat=30*30] BBFG_serial2: 16157ms
[mat=30*30] Ryser_permanent: 31629ms

[mat=31*31] BBFG_serial: 41823ms
[mat=31*31] BBFG_serial2: 35887ms
[mat=31*31] Ryser_permanent: 65591ms

[mat=32*32] BBFG_serial: 88514ms
[mat=32*32] BBFG_serial2: 71134ms
[mat=32*32] Ryser_permanent: 135052ms

[mat=35*35] BBFG_serial: 751871ms
[mat=35*35] BBFG_serial2: 732640ms
[mat=35*35] Ryser_permanent: 1391576ms

jakeffbulmer commented 3 years ago

Hi @maliasadi, I'd be interested to see if there are any numerical accuracy differences between the Ryser and BBFG algorithms. I have a feeling that BBFG may perform better. The simplest test for this I can think of is that the permant of an n x n all-ones matrix, M, should give n!, so computing |perm(M) - n!| / n! for different n should be interesting.

maliasadi commented 3 years ago

Hi @maliasadi, I'd be interested to see if there are any numerical accuracy differences between the Ryser and BBFG algorithms. I have a feeling that BBFG may perform better. The simplest test for this I can think of is that the permant of an n x n all-ones matrix, M, should give n!, so computing |perm(M) - n!| / n! for different n should be interesting.

Hi @jakeffbulmer, although both Ryser and BBFG are among "exact" algorithms and not statistical, like Markov Chain Monte Carlo algorithm for calculating permanent, that's a very good concern. Thanks for the proposed test formula! I gave it a try with these algorithms and the results are as follows. It seems that your feeling is right :)

n: 10   n!: 3628800
bbfg: 3.6288e+06    Ryser: 3.6288e+06
bbfg-diff: 0        Ryser-diff: 0
n: 11   n!: 39916800
bbfg: 3.99168e+07   Ryser: 3.99168e+07
bbfg-diff: 0        Ryser-diff: 0
n: 12   n!: 479001600
bbfg: 4.79002e+08   Ryser: 4.79002e+08
bbfg-diff: 0        Ryser-diff: 0
n: 13   n!: 6227020800
bbfg: 6.22702e+09   Ryser: 6.22702e+09
bbfg-diff: 0        Ryser-diff: 0
n: 14   n!: 87178291200
bbfg: 8.71783e+10   Ryser: 8.71783e+10
bbfg-diff: 0        Ryser-diff: 0
n: 15   n!: 1307674368000
bbfg: 1.30767e+12   Ryser: 1.30767e+12
bbfg-diff: 0        Ryser-diff: 0
n: 16   n!: 20922789888000
bbfg: 2.09228e+13   Ryser: 2.09228e+13
bbfg-diff: 0        Ryser-diff: 0
n: 17   n!: 355687428096000
bbfg: 3.55687e+14   Ryser: 3.55687e+14
bbfg-diff: 3.99136e-16  Ryser-diff: 2.8958e-13
n: 18   n!: 6402373705728000
bbfg: 6.40237e+15   Ryser: 6.40237e+15
bbfg-diff: 0        Ryser-diff: 3.56696e-12
n: 19   n!: 121645100408832000
bbfg: 1.21645e+17   Ryser: 1.21645e+17
bbfg-diff: 1.4011e-15   Ryser-diff: 1.34466e-11
n: 20   n!: 2432902008176640000
bbfg: 2.4329e+18    Ryser: 2.4329e+18
bbfg-diff: 1.31777e-15  Ryser-diff: 1.55759e-10

codecov[bot] commented 3 years ago

Codecov Report

Merging #267 (c03fabf) into master (ecb5aaa) will not change coverage. The diff coverage is 100.00%.

@@            Coverage Diff            @@
##            master      #267   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files           21        21           
  Lines         1418      1419    +1     
=========================================
+ Hits          1418      1419    +1

Impacted Files	Coverage Δ
thewalrus/__init__.py	`100.00% <100.00%> (ø)`
thewalrus/_permanent.py	`100.00% <100.00%> (ø)`

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update ecb5aaa...c03fabf. Read the comment docs.

jakeffbulmer commented 3 years ago

Hi @jakeffbulmer, although both Ryser and BBFG are among "exact" algorithms and not statistical, like Markov Chain Monte Carlo algorithm for calculating permanent, that's a very good concern. Thanks for the proposed test formula! I gave it a try with these algorithms and the results are as follows. It seems that your feeling is right :)

Thanks @maliasadi - that's exactly what I was hoping to see!

maliasadi commented 3 years ago

@mlxd @nquesada: Could you please check and re-run ci/circleci: osx-wheels for me? This failed test may not be related to my changes.

nquesada commented 3 years ago

Hey @maliasadi : could you also add your name here: https://github.com/XanaduAI/thewalrus/blob/master/.github/ACKNOWLEDGMENTS.md

co9olguy commented 3 years ago

Thanks @maliasadi!

nquesada commented 3 years ago

indeed, really nice addition @maliasadi , aka commander of threads :)

maliasadi commented 3 years ago

Thanks @maliasadi!

Couldn't done it without @mlxd and @nquesada suggestions!

XanaduAI / thewalrus

Add support for calculating the permanent using the BBFG algorithm #267

Codecov Report