vector size specified is too large for 32*32 matrix

I have tried the following code in R in order to calculate the permanent of matrix . But I get this error:

Error in vector("list", gamma(n + 1)) : vector size specified is too large.

I tried many solutions but I could not solve the problem.

    library(combinat)
    perm <- function(A)
    {
      stopifnot(is.matrix(A))
      n <- nrow(A)
      if(n != ncol(A)) stop("Matrix is not square.")
      if(n < 1) stop("Matrix has a dimension of size 0.")
      sum(sapply(combinat::permn(n), function(sigma) prod(sapply(1:n, function(i) A[i, sigma[i]]))))
    }
    
    #We copy our test cases from the Python example.
    
    
    testData <- list("A" = rbind(c(0.5,0.3,0.4,0.4,0.4,0.3,0.3,0.4,0.3,0.4,0.5,0.6,0.5,0,0,0.3,0,0,0,0.6,0.6,0.4,0.1,0.2,0.4,0.7,0.4,0.7,0.8,0.9,0.4,0.4),
                                 c(0.7,0.7,0.7,0.7,0.7,0.5,0.6,0.5,0.6,0.7,0.8,0.8,0.8,0,0,0.7,0,0.6,0.7,0.9,0.8,0.7,0.5,0.5,0.7,0,0.7,0.8,0,0.9,0.6,0.7),
                                 c(0.6,0.3,0.5,0.7,0.7,0.4,0.6,0.6,0.6,0.6,0.7,0.7,0.8,0,0,0.6,0,0.6,0.6,0.8,0.8,0.6,0.4,0.5,0.8,0.9,0.6,0.8,0,1,0.6,0.7),
                                 c(0.6,0.3,0.3,0.5,0.7,0.4,0.5,0.7,0.4,0.6,0.7,0.7,0.8,0.4,0.4,0.5,0,0.5,0.6,0.8,0.7,0.5,0.4,0.4,0.7,0.9,0.6,0.8,0,0.9,0.6,0.7),
                                 c(0.6,0.3,0.3,0.3,0.5,0.3,0.4,0.4,0.4,0.4,0.5,0.6,0.6,0,0,0.3,0,0,0,0.7,0.6,0.4,0.3,0.3,0.4,0,0.4,0.6,0,0.7,0.5,0.4),
                                 c(0.7,0.5,0.6,0.6,0.7,0.6,0.6,0.6,0.6,0.6,0.8,0.8,0.8,0,0,0.6,0,0.6,0.6,0.9,0.7,0.6,0.5,0.6,0.7,0.9,0.6,0.8,0,0.9,0.6,0.7),
                                 c(0.7,0.4,0.4,0.5,0.6,0.4,0.4,0.5,0.4,0.5,0.7,0.7,0.7,0,0,0.4,0,0.4,0.7,0.8,0.8,0.6,0.3,0.5,0.7,0.9,0.6,0.7,0,0.9,0.6,0.7),
                                 c(0.6,0.5,0.4,0.3,0.6,0.4,0.5,0.6,0.4,0.5,0.7,0.7,0.7,0,0,0.6,0,0,0,0.8,0.7,0.6,0.4,0.5,0.7,0,0.6,0.7,0,0,0.6,0.6),
                                 c(0.7,0.4,0.4,0.6,0.6,0.4,0.6,0.6,0.5,0.7,0.7,0.7,0.7,0,0,0.6,0,0.5,0.6,0.8,0.7,0.6,0.4,0.5,0.7,0,0.6,0.7,0.9,0.9,0.7,0.8),
                                 c(0.6,0.3,0.4,0.4,0.6,0.4,0.5,0.5,0.3,0.6,0.7,0.7,0.7,0,0,0.6,0,0,0,0.8,0.6,0.6,0.5,0.5,0.6,0.8,0.6,0.6,0,0,0.6,0.7),
                                 c(0.5,0.2,0.3,0.3,0.5,0.2,0.3,0.3,0.3,0.3,0.4,0.6,0.5,0,0,0.4,0,0.3,0.4,0.7,0.5,0.4,0.2,0.2,0.4,0,0.4,0.5,0.7,0.7,0.4,0.4),
                                 c(0.4,0.2,0.3,0.3,0.4,0.2,0.3,0.3,0.3,0.3,0.4,0.4,0.5,0,0,0.3,0,0.4,0.4,0.7,0.5,0.3,0.3,0.2,0.4,0,0.4,0.5,0.8,0.9,0.4,0.4),
                                 c(0.5,0.2,0.2,0.2,0.4,0.2,0.3,0.3,0.3,0.3,0.5,0.5,0.3,0,0,0.3,0,0,0,0.8,0.4,0.4,0.1,0.1,0.4,0,0.4,0.4,0.8,0.8,0.4,0.4),
                                 c(1,1,1,0.6,1,1,1,1,1,1,1,1,1,0.5,0.7,0.7,0,0.7,0.7,0.9,0.7,0.7,0.5,0.6,0.7,0.9,0.7,0.8,0.9,1,0.7,0.8),
                                 c(1,1,1,0.6,1,1,1,1,0,1,1,0,1,0,0.4,0.5,0,0.5,0.5,0.9,0,0.6,0.4,0.5,0.7,0.9,0.6,0,0,0.9,0.6,0.7),
                                 c(0.7,0.3,0.4,0.5,0.7,0.4,0.6,0.4,0,0.4,0.6,0,0.7,0,0,0.6,0,0.6,0.7,0.8,0.8,0.6,0.4,0.5,0.7,0,0.6,0,0,0.9,0.7,0.7),
                                 c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0.3,0.7,0.7,1,0.7,0.7,0.4,0.6,0.7,1,0.7,0.9,1,1,0.7,0.7),
                                 c(1,0.4,0.4,0.5,1,0.4,0.6,1,0.5,1,0.7,0.6,1,0,0,0.4,0,0.2,0.6,0.8,0.7,0.7,0.4,0.5,0.7,0,0.7,0.8,0,1,0.8,0.7),
                                 c(1,0.3,0.4,0.4,1,0.4,0.3,1,0.4,1,0.6,0,1,0,0,0.3,0,0.4,0.1,0.8,0,0.6,0.4,0.6,0.7,0,0.7,0,0,0.9,0.7,0.7),
                                 c(0.4,0.1,0.2,0.2,0.3,0.1,0.2,0.2,0.2,0.2,0.3,0,0.2,0,0,0.2,0,0.2,0.2,0.5,0,0.3,0.1,0.1,0.3,0.7,0.4,0.4,0.7,0.8,0.4,0.4),
                                 c(0.4,0.2,0.2,0.3,0.4,0.3,0.2,0.3,0.3,0.4,0.5,0.5,0.6,0.3,1,0.2,0,0.3,1,1,0.4,0.4,0.3,0.3,0.4,0.8,0.5,0.7,0.9,0.9,0.4,0.4),
                                 c(0.6,0.3,0.4,0.5,0.6,0.4,0.4,0.4,0,0.4,0.6,0,0.6,0,0,0.4,0,0.3,0.4,0.7,0,0.5,0.3,0.4,0.5,0,0.5,0,0,1,0.6,0.6),
                                 c(0.9,0.5,0.6,0.6,0.7,0.5,0.7,0.6,0.6,0.5,0.8,0.7,0.9,0,0,0.6,0,0.6,0.6,0.9,0.7,0.7,0.5,0.7,0.8,0,0.6,0.8,0,0,0.6,0.7),
                                 c(0.8,0.5,0.5,0.6,0.7,0.4,0.5,0.5,0.5,0.5,0.8,0.8,0.9,0,0,0.5,0,0,0,0.9,0.7,0.6,0.3,0.6,0.8,0,0.7,0.8,0.9,0.9,0.7,0.8),
                                 c(0.6,0.3,0.2,0.3,0.6,0.3,0.3,0.3,0.3,0.4,0.6,0.6,0.6,0.3,0.3,0.3,0,0.3,0.3,0.7,0.6,0.5,0.2,0.2,0.4,0.8,0.5,0.7,0.9,0.9,0.5,0.6),
                                 c(0.3,1,0.1,0,0,0.1,0.1,0,0,0,0,0,0,0,0,0,0,0,0,0.3,0,0,0,0,0.2,0.2,0.3,0.3,0.5,0.5,0,0),
                                 c(0.6,0.3,0.4,0,0,0.4,0.4,0.4,0.4,0.4,0.6,0,0.6,0,0,0.4,0,0.3,0.3,0.6,0,0.5,0.4,0.3,0.5,0,0.2,0,0,0.8,0.4,0.5),
                                 c(0.3,0.2,0.2,0.2,0.4,0.2,0.3,0.3,0.3,0.4,0.5,0.5,0.6,0,0,1,0,0,0,0.6,0.3,1,0.2,0.2,0.3,0.7,1,0.4,0.9,0.9,0.5,0.5),
                                 c(0.2,1,0,0,0,1,0,1,0.1,1,0.3,0.2,0.2,0,0,0,0,0,0,0,0,0,1,0.1,0.1,0.5,0,0.1,0.5,0.5,0,0),
                                 c(0.1,0.1,0,0,0,0.1,0,1,0.1,1,0.3,0.1,0.2,0,0,0,0,0,0,0,0,0,1,0.1,0.1,0.5,0,0.1,0.5,0.4,0,0),
                                 c(0.6,0.4,0.4,0.4,0.5,0.4,0.4,0.4,0.3,0.4,0.6,0.6,0.6,0,0,0.3,0,0,0,0.6,0.6,0.4,0.4,0.3,0.5,0,0.6,0.5,1,1,0.4,0.7),
                                 c(0.6,0.3,0.3,0,0,0.3,0.3,0.4,0.2,0.3,0.6,0,0.6,0,0,0.3,0,0,0.3,0.6,0,0.4,0.3,0.2,0.4,0,0.5,0.5,0,1,0.3,0.3)))
                                       
print(sapply(testData, function(x) list( Permanent = perm(x))))

Solution

Calculating matrix permanents is a computationally hard problem (e.g. "Computing the permanent" has its own dedicated Wikipedia article separate from the mathematical definition of the operation ...)

The naive method you've implemented (and that is discussed here) will require computing a vector of length 32! = 2.6e35 (which will take more than 10^24 terabytes of memory, let alone the time taken to do all of the computations ...)

Lucky for you, since the last time this was discussed on SO (in 2014), someone has posted the BosonSampling package to CRAN, which uses compiled code (and an efficient algorithm!) to do the computation. tl;dr BosonSampling::rePerm(testData$A) completes in less than 2 minutes on my laptop.

setup

## repeated definition from above, for convenience
library(BosonSampling)
library(combinat)
perm <- function(A)  {
    stopifnot(is.matrix(A))
    n <- nrow(A)
    if(n != ncol(A)) stop("Matrix is not square.")
    if(n < 1) stop("Matrix has a dimension of size 0.")
    sum(sapply(combinat::permn(n), function(sigma) prod(sapply(1:n, function(i) A[i, sigma[i]]))))
}

check equivalence of naive and sophisticated solutions on small data

set.seed(101)
t1 <- matrix(rnorm(16), 4, 4)
t2 <- matrix(rnorm(36), 6, 6)
stopifnot(all.equal(perm(t1), rePerm(t1)))
stopifnot(all.equal(perm(t2), rePerm(t2)))

evaluate time scaling

(this takes a few minutes on my laptop)

f <- function(n, type = "naive") {
    set.seed(101)
    m <- matrix(rnorm(n^2), n, n)
    t <- system.time(if (type == "naive") perm(m) else rePerm(m))
    return(t[["elapsed"]])
}
n0 <- 3; n1 <- 10; n2 <- 32
ntimes <- sapply(n0:n1, f)
ptimes <- sapply(n0:n2, f, type = "rePerm")

dd <- data.frame(n = c(n0:n1, n0:n2),
                 time = c(ntimes, ptimes),
                 type = rep(c("naive", "rePerm"), c(n1-n0+1, n2-n0+1)))
saveRDS(dd, file = "SO72004136.rds")
library(ggplot2); theme_set(theme_bw())
(ggplot(dd) +
 aes(n, time, colour = type) +
 geom_point() + geom_line() +
 scale_y_log10() +
 scale_x_log10() +
 expand_limits(x=32) +
 labs(y = "elapsed time (sec)") +
 geom_smooth(method = "lm",
             lty = 2,
             colour = "black",
             data = subset(dd, type == "rePerm" & n>16))
)

lm(log(time) ~ log(n),
   data = subset(dd, type == "rePerm" & n>16))

Running the naive approach for n=10 takes 190 seconds; rePerm for n=29 takes 18 seconds. The computational time scales roughly as time \propto n^17 (i.e. the slope of the log-log regression is approximately 17), which is frighteningly bad ...

Once I had figured out that the problem was likely to take minutes rather than hours, I computed rePerm(testData$A) = 1.271369e+20 in 115 seconds. (I sincerely hope you don't have to do computations much bigger than this!)