4

How can I find a permanent of a square matrix (for a general dimension nxn) in R? In particular, I'm trying to find a pdf of order statistics for independent but not identically distributed populations, which includes calculating a permanent of a matrix whose elements are pdfs and cdfs of the original populations

thanks

Alexandra Bolotskikh
  • 59
  • 3
  • 2
    For those of us with insufficient background, would you like to explain what a "permanent" is. (I'm wondering if you meant to write "determinant", in whcihc case you would probably want to delete this question and just do a search using the standard `??`-facilities in R.) – IRTFM Jun 11 '14 at 23:00
  • 1
    http://en.wikipedia.org/wiki/Permanent – Ben Bolker Jun 11 '14 at 23:02
  • It's basically like determinant but each term comes in with a + sign – Alexandra Bolotskikh Jun 11 '14 at 23:13
  • whether @BondedDust's brute-force solution is adequate, or whether you will need to implement Ryser's algorithm in R, or wrap a C implementation (and whether you need to modify the C implementation pointed to below) will depend on how big your matrices are and whether they are integer or floating-point ... – Ben Bolker Jun 12 '14 at 00:08

2 Answers2

4

tl;dr this is a non-trivial computational problem that doesn't seem to have been implemented in R, and is computationally hard enough that a compiled solution may be necessary. Your best bet would be to write R code wrapping this open source C implementation.

Based on the relevant Wikipedia article, "Ryser" looks like a good search term for finding implementations of this computation. library("sos"); findFn("Ryser") finds only the help for Spearman's rank correlation, which says

Calculation of the exact null distribution of Spearman's rank correlation statistics is exponentially hard in n. This package uses precomputed exact distribution for n <= 22 obtained using Ryser's formula applied to an appropriate monomial permanent.

This isn't even a general implementation, but a special case. Googling "permanent Ryser" doesn't come up with any implementations until we get down to here, which is MATLAB code. Googling "permanent Ryser implementation" comes up with this CodeProject page, which gives fairly straightforward C code licensed under the fairly permissive Code Project Open License.

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
2

The permanent of matrix(1:9, 3) would then just be:

 install.packages("permute"); library(permute)
  A<-matrix(1:9, 3)

  # Error: sum( apply( allPerms(1:3), 1, function(r) prod( A[1:3, r]) )  )

The allPerms function seems to leave out the original vector, hence the need for one of Ben Bolker's corrections and I should have used cbind to construct the indices for the items of A:

sum( apply( rbind(1:3,allPerms(1:3)), 1,
                               function(r) prod( A[cbind(1:3, r)]) ) )

The fact that these values are all positive and there is no subtraction suggests a reason why this "naive" implementation of that definition is not recommended.

 A <- matrix(1:16,4)    
 sum( apply( rbind(1:4,allPerms(1:4)), 1, 
                       function(r) prod( A[cbind(1:4, r)]) ) )
#[1] 55456
Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
IRTFM
  • 258,963
  • 21
  • 364
  • 487
  • based on the wikipedia page I think the right answer is `1*5*9+3*4*8+2*6*7+3*5*7+6*8*1+2*4*9`=450? – Ben Bolker Jun 11 '14 at 23:31
  • Then I guess I misunderstood the subscripting? Or Gavin's code doesn't deliver the right number of permutations of three items. You have 6 and I got 5. (And now I'm wondering why there are not 6....) – IRTFM Jun 11 '14 at 23:36
  • Hmm ... corrected: `sum( apply( rbind(1:3,allPerms(1:3)), 1, function(r) prod( A[cbind(1:3, r)]) ) )` gives 450. (`allPerms` doesn't include the original vector, *and* the subscripting is different.) – Ben Bolker Jun 11 '14 at 23:36