CODEWITHSUNDEEP
rmatrixpermanent
Alexandra Bolotskikh
Alexandra Bolotskikh

Reputation: 69

Calculate permanent of a matrix in R

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

Upvotes: 5

Views: 1248

Answers (3)

ThomasIsCoding
ThomasIsCoding

Reputation: 102599

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

Upvotes: 0

IRTFM
IRTFM

Reputation: 263451

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

Upvotes: 3

Ben Bolker
Ben Bolker

Reputation: 226761

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.

Upvotes: 5

Related Questions