Calculate permanent of a matrix in R

Question

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

Preface

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.) – Nonfeasance

en.wikipedia.org/wiki/Permanent – Poetaster

It's basically like determinant but each term comes in with a + sign – Preface

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 ... – Poetaster

Answer 1

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.

Poetaster

Answer 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

Nonfeasance

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? – Poetaster

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....) – Nonfeasance

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.) – Poetaster

Answer 3

You can try the recursion version of "permanent" implementation

permanent <- function(x) {
    n <- nrow(x)
    s <- 0
    if (n > 2) {
        for (j in 1:n) {
            s <- s + x[1, j] * Recall(x[-1, -j])
        }
        return(s)
    } else {
        x[1, 1] * x[2, 2] + x[1, 2] * x[2, 1]
    }
}

such that

> A <- matrix(1:9, 3)

> permanent(A)
[1] 450

and

> A <- matrix(1:16, 4)

> permanent(A)
[1] 55456

Coca