+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: HP Prime Software Library (/forum-15.html)
+--- Thread: permanent of square matrix (/thread-21097.html)
permanent of square matrix - Albert Chan - 01-02-2024 10:41 PM
Computation of the permanent of a matrix #7151
0-based numpy npperm(M) translated to 1-based HP Prime cas code
Code:
#cas
permanent(M) :=
BEGIN
LOCAL n,d,j,s,f,v,p;
n := dim(M)[j:=s:=1];
d := makelist(2,1,n);
f := range(1,n+1);
v := map(sum, transpose(M));
p := product(v);
WHILE j < n DO
v -= d[j]*M[j]
d[j] := -d[j]
s := -s
p += s*product(v);
f[1] := 1;
f[j] := f[j+1];
f[j+1] := j+1;
j := f[1];
END;
RETURN p/2**(n-1)
END;
#endCode does not generate sub-matrixes, very efficient!
> m := ranm(9,9)
Code:
[[ 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]]> permanent(m) /* time = 0.008 sec */
−3500631868051794540
Lets try naive implementation
> per(m) /* time = 95.68 sec */
−3500631868051794540
RE: permanent of square matrix - Albert Chan - 01-04-2024 01:13 PM
Translated above to HP71B BASIC
Code:
10 DESTROY ALL @ OPTION BASE 1 @ INPUT "N? ";N
20 DIM M(N,N),D(N),F(N),V(N) @ MAT INPUT M @ T0=TIME
30 MAT D=(2) @ MAT V=CSUM(M) @ P=1 @ S=1
40 FOR I=1 TO N @ F(I)=I @ P=P*V(I) @ NEXT I
50 WHILE F(1)<N @ I=F(1) @ T=1
60 FOR J=1 TO N @ V(J)=V(J)-D(I)*M(I,J) @ T=T*V(J) @ NEXT J
70 P=P-S*T @ S=-S @ D(I)=-D(I) @ F(1)=1 @ F(I)=F(I+1) @ F(I+1)=I+1
80 END WHILE
90 P=P/2^(N-1) @ DISP "P=";P,TIME-T0OP 9×9 random matrix example
Code:
N? 9
M(1,1)? 86, -97, -82, 7, -27, 26, -89, 63, -49
M(2,1)? -86, -64, -30, 70, 22, 42, 56, -9, -13
M(3,1)? 46, 24, 49, -6, 0, 81, 18, -49, -33
M(4,1)? -70, 8, 63, -64, 2, 62, -37, -80, -23
M(5,1)? 65, -85, 28, -44, -22, 93, 91, 31, -21
M(6,1)? 88, 76, -66, 66, 5, -23, 79, -88, 9
M(7,1)? 6, -69, -8, 31, 89, 2, 97, -92, 80
M(8,1)? -24, -17, 41, 64, -49, 0, 89, 18, -51
M(9,1)? -42, -65, 22, -94, 76, -10, 16, -63, 39
P=-3.50063186806E18 3.12Emulator running at 40x --> Physical HP71B take 3.12*40 ≈ 125 s ≈ 2 min
RE: permanent of square matrix - Albert Chan - 01-04-2024 08:06 PM
(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:
(
+( a11+a21+a31)*( a12+a22+a32)*( a13+a23+a33)
−( a11-a21+a31)*( a12-a22+a32)*( a13-a23+a33)
+(-a11-a21+a31)*(-a12-a22+a32)*(-a13-a23+a33)
−(-a11+a21+a31)*(-a12+a22+a32)*(-a13+a23+a33)
) / 4First 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
+a11*a23*a32
+a12*a21*a33
+a12*a23*a31
+a13*a21*a32
+a13*a22*a31(*) 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
permanent(M) :=
BEGIN
LOCAL n,d,j,f,v,p;
if (n:=len(M))<2 THEN RETURN M[1][1] END
d := makelist(2,2,n);
f := range(1,n+1);
v := map(sum,transpose(M));
p := product(v);
WHILE (j:=f[1]) < n DO
v -= d[j]*M[j]
p := product(v) - p;
d[j] := -d[j]
f[1] := 1;
f[j] := f[j+1];
f[j+1] := j+1;
END;
RETURN -p/2**(n-1)
END;
#endRE: permanent of square matrix - Albert Chan - 01-04-2024 10:59 PM
(01-04-2024 08:06 PM)Albert Chan Wrote: 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.
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,
between any two adjacent numbers, only one bit changes at a time.
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
RE: permanent of square matrix - John Keith - 01-07-2024 07:47 PM
(01-04-2024 08:06 PM)Albert Chan Wrote: I am beginning to understand why d cells = ±2
...
permanent(M), version 2
Code:
#cas
permanent(M) :=
BEGIN
LOCAL n,d,j,f,v,p;
if (n:=len(M))<2 THEN RETURN M[1][1] END
d := makelist(2,2,n);
f := range(1,n+1);
v := map(sum,transpose(M));
p := product(v);
WHILE (j:=f[1]) < n DO
v -= d[j]*M[j]
p := product(v) - p;
d[j] := -d[j]
f[1] := 1;
f[j] := f[j+1];
f[j+1] := j+1;
END;
RETURN -p/2**(n-1)
END;
#end
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:
\<< DUP 2 * @ Multiply matrix by 2
AXL DUP SIZE 1. @ Convert matrix to nested list
IF > @ Size > 1?
THEN SWAP TRN AXL @ List of columns of original matrix
:: \GSLIST DOSUBS @ Column sums
DUP \PILIST UNROT @ Product of column sums
DUP SIZE 2. OVER 1. - ^ \-> n w @ n = size, w = 2^(n-1)
\<< 1. w 1. - @ 2^(n-1) - 1 iterations
FOR k 1. k @ Row index = 1
WHILE DUP 2. MOD NOT @ While even
REPEAT 2. / SWAP 1. + SWAP @ Divide by 2 and increment index
END 4. PICK ROT GET @ Get matrix row
SWAP 4. MOD 3. SAME @ Odd part MOD 4 = 3?
:: ADD :: - IFTE @ Add or subtract sums
DUP \PILIST 4. ROLL - UNROT @ Update product
NEXT DROP2 @ Drop lists
NEG w R\->I / @ Return negative product / w
\>>
ELSE DROP 1. GET @ Else drop list, get single element
END
\>>RE: permanent of square matrix - Albert Chan - 01-07-2024 11:47 PM
(01-07-2024 07:47 PM)John Keith Wrote: 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.
Translated to Lua. Yes, code very compact!
Code:
function permanent(M)
local p,n,v = 1,#M,{}
if n==1 then return M[1][1] end
for i = 1,n do
local t = 0
for j = 1,n do t = t+M[j][i] end
v[i], p = t, p*t
end
local terms = 2^(n-1)
for k = 1,terms-1 do
local j, t = bit.ffs(k), 1
local d = bit.test(k,j) and -2 or 2
for i = 1,n do v[i] = v[i]-d*M[j][i]; t = t*v[i] end
p = t - p
end
return -p/terms
endbit.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
RE: permanent of square matrix - John Keith - 01-09-2024 09:49 PM
(01-04-2024 01:13 PM)Albert Chan Wrote: Translated above to HP71B BASIC
Code:
10 DESTROY ALL @ OPTION BASE 1 @ INPUT "N? ";N
20 DIM M(N,N),D(N),F(N),V(N) @ MAT INPUT M @ T0=TIME
30 MAT D=(2) @ MAT V=CSUM(M) @ P=1 @ S=1
40 FOR I=1 TO N @ F(I)=I @ P=P*V(I) @ NEXT I
50 WHILE F(1)<N @ I=F(1) @ T=1
60 FOR J=1 TO N @ V(J)=V(J)-D(I)*M(I,J) @ T=T*V(J) @ NEXT J
70 P=P-S*T @ S=-S @ D(I)=-D(I) @ F(1)=1 @ F(I)=F(I+1) @ F(I+1)=I+1
80 END WHILE
90 P=P/2^(N-1) @ DISP "P=";P,TIME-T0
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.
RE: permanent of square matrix - Werner - 01-10-2024 08:32 AM
(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
RE: permanent of square matrix - Albert Chan - 01-10-2024 02:50 PM
(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)
local p,n,v = -2,#M,{}
if n==1 then return M[1][1] end
for i = 1,n do
local t = 0
for j = 1,n do t = t+M[j][i] end
v[i] = t/2
p = p*v[i]
end
for k = 1,2^(n-1)-1 do
local j, t = bit.ffs(k), -2
local d = bit.test(k,j) and -1 or 1
for i = 1,n do v[i] = v[i]-d*M[j][i]; t = t*v[i] end
p = t - p
end
return p
endDebug 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.
RE: permanent of square matrix - John Keith - 01-24-2024 04:03 PM
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.
RE: permanent of square matrix - Albert Chan - 01-25-2024 03:49 PM
(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)
local p,n, d,f,v = 1,#M, {},{},{}
if n==1 then return M[1][1] end
for j = 1,n do
d[j] = 2
f[j] = j
local t = 0
for i = 1,n do t = t+M[i][j] end -- sum col
v[j], p = t, p*t
end
while true do
local i = f[1]
if i >= n then return -p/2^(n-1) end
local di, r, t = d[i], M[i], 1
for j=1,n do v[j] = v[j]-di*r[j]; t = t*v[j] end
p = t - p
d[i] = -d[i]
f[1] = 1
f[i] = f[i+1]
f[i+1] = i+1
end
endI gave up on the 5% speedup, for cleaner code.
Still, for OP 9×9 matrix, its speed is doubled the bitwise version. (post #9)
RE: permanent of square matrix - John Keith - 01-25-2024 09:18 PM
(01-25-2024 03:49 PM)Albert Chan Wrote: I guess it depends on programming language used.
...
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!)
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.
RE: permanent of square matrix - Werner - 01-26-2024 09:25 AM
(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
RE: permanent of square matrix - John Keith - 01-26-2024 10:08 PM
(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.
In the 42S, I decided not to do that as multiplying the matrix by 2 inevitably duplicates it, and memory is relatively tight.
Werner
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.
RE: permanent of square matrix - Albert Chan - 02-14-2024 03:14 PM
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
for the calculation of per(A). We reduce this labor in three steps.
(1) A method due to Ryser evaluates per(A) in about (n^2/2) * 2^n operations.
It requires an average of n^2/2 calculations for each of the 2^n subsets of {1,2,...,n}.
Ryser’s method is derived below and appears in Eqs. (17) and (18).
(2) We reduce the above by a factor of 2 by a method which makes it necessary to process
only the subsets of {1,2,...,n-1}. Thus in Eq. (24) we do about n^2/2 calculations for
each of the 2^(n-1) subsets of {1,2,... , n-1]}, or (n^2/4) * 2^n operations altogether.
(3) A further reduction by a factor of n/2 is accomplished by arranging the sequence of subsets
so as to follow a Hamilton walk on an (n-1)-cube (see Chapter 1), If this is done, only a slight
change in the calculation for a set S will yield the result of the calculation for the successor of S.
In this way, our final algorithm PERMAN does about n calculations for each of the 2^(n-1)
subsets of {1,2,...,n-1}, for a total of n * 2^(n-1) operations (multiplication and addition)
to compute the permanent of an nxn matrix.
RE: permanent of square matrix - Namir - 02-15-2024 01:25 PM
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
RE: permanent of square matrix - John Keith - 02-15-2024 06:14 PM
I don't think that the permanent has been "kept secret", the determinant is just more useful and easier to compute.
RE: permanent of square matrix - Namir - 02-15-2024 08:08 PM
(02-15-2024 06:14 PM)John Keith Wrote: I don't think that the permanent has been "kept secret", the determinant is just more useful and easier to compute.
None of my math teachers ever discussed the permanent of a matrix -- unlike the detreminant.
Namir