permanent of square matrix
|
01-02-2024, 10:41 PM
(This post was last modified: 01-02-2024 11:22 PM by Albert Chan.)
Post: #1
|
|||
|
|||
permanent of square matrix
Computation of the permanent of a matrix #7151
0-based numpy npperm(M) translated to 1-based HP Prime cas code Code: #cas Code does not generate sub-matrixes, very efficient! > m := ranm(9,9) Code: [[ 86, -97, -82, 7, -27, 26, -89, 63, -49], > permanent(m) /* time = 0.008 sec */ −3500631868051794540 Lets try naive implementation > per(m) /* time = 95.68 sec */ −3500631868051794540 |
|||
01-04-2024, 01:13 PM
Post: #2
|
|||
|
|||
RE: permanent of square matrix
Translated above to HP71B BASIC
Code: 10 DESTROY ALL @ OPTION BASE 1 @ INPUT "N? ";N OP 9×9 random matrix example Code: N? 9 Emulator running at 40x --> Physical HP71B take 3.12*40 ≈ 125 s ≈ 2 min |
|||
01-04-2024, 08:06 PM
Post: #3
|
|||
|
|||
RE: permanent of square matrix
(01-02-2024 10:41 PM)Albert Chan Wrote: d := makelist(2,1,n); I am beginning to understand why d cells = ±2 ±aij - (±2) * aij = ∓aij Next time we process the same row, we want to flip it back, so ±2 --> ∓2 ∓aij - (∓2) * aij = ±aij > permanent([[a11,a12,a13],[a21,a22,a23],[a31,a32,a33]]) Code: ( First term = product of (v = column sums) Other terms, v -= d[j]*M[j], flipped only sign of M j-th row. This make it very efficient, instead of summing them all from scratch. sign flips caused *all* same column products to cancel (nice!) We are left with permutations with different columns (all with plus sign!) 4 terms generated too many copies, thus shrink by 2^(3-1) = 4 correction. (*) > expand(Ans) Code: +a11*a22*a33 (*) Because correction is always pow-of-2. We may eliminate variable s of ±1 t1 - t2 + t3 - t4 = -(t4 - (t3 - (t2 - t1))) Also, M last row sign never flipped, d only need (n-1) cells permanent(M), version 2 Code: #cas |
|||
01-04-2024, 10:59 PM
Post: #4
|
|||
|
|||
RE: permanent of square matrix
(01-04-2024 08:06 PM)Albert Chan Wrote: First term = product of (v = column sums) Selection of what next j-th row is by idea of Gray Code Quote:Gray Code is an alternative binary representation, cleverly devised so that, Lets examine 5×5 matrix sign flip pattern, to confirm indeed it uses Gray Code. Again, last row sign never flipped, thus only 2^(5-1) = 16 sign flip patterns. I reversed bits, to match gray code pattern exactly. lua> n = 5 lua> f = {1,2,3,4,5} lua> b = {0,0,0,0} lua> function bits(k) b[k]=1-b[k]; print(table.concat(b)) end lua> while f[1]<n do j=f[1]; f[1]=1; f[j]=f[j+1]; f[j+1]=j+1; bits(n-j) end 0001 0011 0010 0110 0111 0101 0100 1100 1101 1111 1110 1010 1011 1001 1000 |
|||
01-07-2024, 07:47 PM
(This post was last modified: 01-24-2024 03:46 PM by John Keith.)
Post: #5
|
|||
|
|||
RE: permanent of square matrix
(01-04-2024 08:06 PM)Albert Chan Wrote: I am beginning to understand why d cells = ±2 That is fantastic! Much faster than any other method that I have seen. Even better, it turns out that we can eliminate the arrays d and f entirely, making the program shorter and a bit faster. First of all, we can notice that the number of iterations performed by the program is 2^(n-1)-1, so we can change the WHILE loop to a FOR loop and keep track of which iteration we are in. What f[1], and thus j are actually computing is which bit is flipped to generate the next Gray code. This value is the 2-adic valuation of k, where k is the index of the current iteration. See A001511 and A007814. Determining which values of d[j] are negative is a bit harder but once again the OEIS comes to the rescue. The value of d[j] is negative iff the odd part of k is of the form 4k+3. See A091067. Both of these values can be found by testing the least significant bit of k and if 0, shifting right until the lsb is a 1. If binary numbers are inconvenient, we can divide by 2 and test the result MOD 2. The number of shifts required is the valuation, and the remaining value after shifting is the odd part. If the odd part MOD 4 = 3, d[j] is negative, else it is positive. The following program is a rough translation of your Prime program implementing the changes described above. The program converts the matrix to a nested list to take advantage of list processing commands. Code:
|
|||
01-07-2024, 11:47 PM
Post: #6
|
|||
|
|||
RE: permanent of square matrix
(01-07-2024 07:47 PM)John Keith Wrote: That is fantastic! Much faster than any other method that I have seen. Translated to Lua. Yes, code very compact! Code: function permanent(M) bit.ffs(x) internally call gcc __builtin_ffs(x) int __builtin_ffs (int x) Wrote:Returns one plus the index of the least significant 1-bit of x, or if x is zero, returns zero. Test with OP 9×9 random numbers lua> m = { : { 86, -97, -82, 7, -27, 26, -89, 63, -49}, : {-86, -64, -30, 70, 22, 42, 56, -9, -13}, : { 46, 24, 49, -6, 0, 81, 18, -49, -33}, : {-70, 8, 63, -64, 2, 62, -37, -80, -23}, : { 65, -85, 28, -44, -22, 93, 91, 31, -21}, : { 88, 76, -66, 66, 5, -23, 79, -88, 9}, : { 6, -69, -8, 31, 89, 2, 97, -92, 80}, : {-24, -17, 41, 64, -49, 0, 89, 18, -51}, : {-42, -65, 22, -94, 76, -10, 16, -63, 39}} lua> lua> permanent(m) /* time = 0.00019 s */ -3500631868051796000 |
|||
01-09-2024, 09:49 PM
Post: #7
|
|||
|
|||
RE: permanent of square matrix
(01-04-2024 01:13 PM)Albert Chan Wrote: Translated above to HP71B BASIC This brings up a problem with this algorithm. The variable p is much larger than the final value, up to 3 digits larger for the sizes of matrices we are dealing with. This limits the size of the largest matrix we can use on the 71, 28 and 48 which are limited to 12-digit numbers. Valentin's program is not as fast, but it does not have this problem. Perhaps there exists a happy medium somewhere, but the math behind these algorithms is way above my head. |
|||
01-10-2024, 08:32 AM
Post: #8
|
|||
|
|||
RE: permanent of square matrix
(01-09-2024 09:49 PM)John Keith Wrote: The variable p is much larger than the final value, up to 3 digits larger for the sizes of matrices we are dealing with. This limits the size of the largest matrix we can use on the 71, 28 and 48 which are limited to 12-digit numbers. Valentin's program is not as fast, but it does not have this problem. Perhaps there exists a happy medium somewhere, but the math behind these algorithms is way above my head. Just do v/2, use dj=+/-1 again and double the result. If all vi are odd, this will not help ;-) Cheers, Werner 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE |
|||
01-10-2024, 02:50 PM
(This post was last modified: 01-10-2024 08:20 PM by Albert Chan.)
Post: #9
|
|||
|
|||
RE: permanent of square matrix
(01-10-2024 08:32 AM)Werner Wrote: Just do v/2, use dj=+/-1 again and double the result. If all vi are odd, this will not help ;-) Numerical issue is not 2^(n-1) scaling factor (even for OP 9×9 matrix, factor is only 2^8=256). It is really catastrophic cancellation of p's during its iterations. For illustration purpose, I scaled lua code final result of -p/terms to simply return p For base 2, this scaling served no purpose, except reduce tiny chance of overflow. For base 10, this scaling may make numerical issue worse (some vi may be odd) Code: function permanent(M) Debug with print((-1)^k*p) (★) inside for k loop, tested with OP 9×9 random matrix. lua> permanent(m) -5556546744184800 20429002103295824 55255400982040850 31542151598286920 ... 314811585639751230 -119381930932063740 -81181015636135740 -81169499916475780 -73135689828461120 -241381971380131900 248282227143724100 ... 528128789659814300 524276224665997300 -2888266392290130400 -2753810694416846000 ... -3101014455289180000 -3584102548532368000 -3432517605155301000 -3500631868051796000 -- final p = permanent(m) This is probably a bad example, with final p very close to true result. (error = 3 ULP) But if true result is tiny, we could see why this algorithm have cancellation issues. (★) (-1)^k factor to account for p = t - p sign reversal. In other words, if rest of t's contributed nothing, this is the final p we get. |
|||
01-24-2024, 04:03 PM
(This post was last modified: 01-24-2024 04:05 PM by John Keith.)
Post: #10
|
|||
|
|||
RE: permanent of square matrix
I noticed an inefficiency in the original Python program which has carried through our subsequent programs. This line, v -= 2*d[j]*M[j] results in an array (or list) multiplication in every iteration. If the matrix is multiplied by 2 (after computing the column sums), we can add the currently selected row to v if k MOD 4 = 3, else subtract it. I have updated my RPL program in post #5 to reflect this change, which results in a 33% speedup.
|
|||
01-25-2024, 03:49 PM
Post: #11
|
|||
|
|||
RE: permanent of square matrix
(01-24-2024 04:03 PM)John Keith Wrote: This line, v -= 2*d[j]*M[j] results in an array (or list) multiplication in every iteration. I guess it depends on programming language used. Python was using numpy for fast v list update, rewritten it without multiply does not gain much. Note: 2*d[j] is just a number, not a list For LuaJIT, the cost is smaller still (OP 9×9 matrix, speedup under 5%) Array access cost a lot more than multiply. Bit operations cost even more. (this is just relatively speaking ... LuaJIT is FAST!) Below LuaJIT version is the one I kept. Code: function permanent(M) I gave up on the 5% speedup, for cleaner code. Still, for OP 9×9 matrix, its speed is doubled the bitwise version. (post #9) |
|||
01-25-2024, 09:18 PM
Post: #12
|
|||
|
|||
RE: permanent of square matrix
(01-25-2024 03:49 PM)Albert Chan Wrote: I guess it depends on programming language used. 2*d[j] is just a number, but it is multiplied by row j of the matrix, requiring n multiplies in every iteration for an n X n matrix. Of course, the amount of improvement depends on the language and the hardware that the program is run on. Considering that RPL is an interpreted language running on an 80's era CPU (real or emulated) one would expect a greater degree of improvement than a compiled language running on a modern CPU. I have not tried the program on the Prime, but I expect that the degree of improvement would be somewhere in between. |
|||
01-26-2024, 09:25 AM
Post: #13
|
|||
|
|||
RE: permanent of square matrix
(01-25-2024 03:49 PM)Albert Chan Wrote: Note: 2*d[j] is just a number, not a list Yes, but d[j] is +/-1, and he performs an addition or subtraction accordingly, no list/vector multiply necessary. In the 42S, I decided not to do that as multiplying the matrix by 2 inevitably duplicates it, and memory is relatively tight. Werner 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE |
|||
01-26-2024, 10:08 PM
(This post was last modified: 01-26-2024 10:09 PM by John Keith.)
Post: #14
|
|||
|
|||
RE: permanent of square matrix
(01-26-2024 09:25 AM)Werner Wrote: Yes, but d[j] is +/-1, and he performs an addition or subtraction accordingly, no list/vector multiply necessary. Also worth noting that the column sums are computed from the un-multiplied matrix so p does not get any larger than it does in the original program. I would think that computing the permanent on a 42S would be glacially slow, but should be fine on Free42 with greatly increased speed and 34-digit precision. |
|||
02-14-2024, 03:14 PM
Post: #15
|
|||
|
|||
RE: permanent of square matrix
I've found an old (1978) FREE book, described how permanent algorithm was designed.
Combinatorial Algorithms for computers and calculators, 2nd edition, by Herbert S. Wilf, Albert Nijenhuis Chapter 28, Computation of the Permanent Function Wrote:Observe that a direct application of (3) to an nxn matrix would require about n * n! operations |
|||
02-15-2024, 01:25 PM
(This post was last modified: 02-15-2024 01:26 PM by Namir.)
Post: #16
|
|||
|
|||
RE: permanent of square matrix
How come the permanent of a matrix has been kept a virtual secret in matrix math (especialy numerical analysis) books? We all heard about the determiant and how it can tell if a matrix is invertible, but nothing about the permanent of a matrix.
Namir |
|||
02-15-2024, 06:14 PM
Post: #17
|
|||
|
|||
RE: permanent of square matrix
I don't think that the permanent has been "kept secret", the determinant is just more useful and easier to compute.
|
|||
02-15-2024, 08:08 PM
Post: #18
|
|||
|
|||
RE: permanent of square matrix | |||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 6 Guest(s)