First off, you are probably better off rethinking your approach. With that said, let's get started attacking your problem.
This is a very difficult problem mainly due to the limitations of resources. Below, I have a solution that will complete in a reasonable amount of time on a home computer given that you have access to a decent amount of storage (at least 7 TB). The algorithm below does not require that much memory and can be tuned to reduce memory usage.
Before we begin, we note that merely generating that many permutations seems impossible at first. However with the help of highly optimized C++ code and parallel computing, the task is brought back into the realm of possibility. This was demonstrated in my answer to the OP's previous question. We utilized RcppAlgos (I am the author) and the parallel package to generate about 36 million permutations per second in chunks of one million using 8 cores.
Now, we are charged with carrying out specific computations on each permutation as fast as possible. The computation is as follows:
(SLS') / 4, where S is a permutation, L is a 41 x 41 matrix
Here are a couple of base R approaches (N.B. m1[x, ] %*% m2 %*% m1[x, ] is the same as m1[x, ] %*% m2 %*% as.matrix(m1[x, ], ncol = 1)):
baseTest1 <- function(m1, m2) {
vapply(1:nrow(m1), function(x) {
m1[x, ] %*% m2 %*% m1[x, ]
}, FUN.VALUE = 1.1111, USE.NAMES = FALSE) / 4
}
baseTest2 <- function(m1, m2) {
temp <- m1 %*% m2
vapply(1:nrow(m1), function(x) {
crossprod(temp[x, ], m1[x, ])
}, FUN.VALUE = 1.1111, USE.NAMES = FALSE) / 4
}
Let's think about this a little bit. We have a bunch of permutations of the numbers one and negative one. When we multiply these permutations by a matrix of real numbers, say M, we end up simply adding and subtracting values from the M. I bet we can speed this up quite a bit using Rcpp and avoid wasteful (and useless) identity multiplications (i.e. multiplications by 1).
#include <Rcpp.h>
//[[Rcpp::export]]
Rcpp::NumericVector makeVecCpp(Rcpp::NumericMatrix A,
Rcpp::NumericMatrix B,
unsigned long int mySize) {
Rcpp::NumericVector result = Rcpp::no_init_vector(mySize);
double temp = 0;
for (std::size_t i = 0; i < mySize; ++i) {
for (std::size_t j = 0; j < 41u; ++j) {
for (std::size_t k = 0; k < 41u; ++k) {
if (A(i, j) + A(i, k)) {
temp += B(j, k);
} else {
temp -= B(j, k);
}
}
}
result[i] = temp / 4;
temp = 0;
}
return result;
}
Now let's see if they give the same results and also benchmark them:
options(scipen = 999)
library(RcppAlgos)
library(microbenchmark)
set.seed(42)
M <- matrix(rnorm(41*41), nrow = 41, ncol = 41)
negOne <- permuteGeneral(c(1L, -1L), freqs = c(21, 20), upper = 100000)
all.equal(baseTest1(negOne, M), baseTest2(negOne, M))
# [1] TRUE
all.equal(baseTest1(negOne, M), makeVecCpp(negOne, M, 100000))
# [1] TRUE
microbenchmark(base1 = baseTest1(negOne, M), base2 = baseTest2(negOne, M),
myRcpp = makeVecCpp(negOne, M, 100000), times = 25)
Unit: milliseconds
expr min lq mean median uq max neval
base1 555.0256 582.2273 597.6447 593.7708 599.1380 690.3882 25
base2 471.0251 494.2367 541.2632 531.1858 586.6774 632.7279 25
myRcpp 202.7637 207.2463 210.0255 209.0399 209.9648 240.6664 25
Our Rcpp implementation is the clear winner!! Moving on, we incorportate this into our final answer:
## WARNING Don't run this unless you have a few DAYS on your hand
library(parallel)
## break up into even intervals of one hundred thousand
firstPart <- mclapply(seq(1, 269128900000, 100000), function(x) {
negOne <- permuteGeneral(c(1L, -1L), freqs = c(21, 20),
lower = x, upper = x + 99999)
vals <- makeVecCpp(negOne, M, 100000)
write.csv(vals, paste0("myFile", x, ".csv", collapse = ""))
x
}, mc.cores = 8)
## get the last few results and complete analysis
lastPart <- permuteGeneral(c(1L, -1L), freqs = c(21, 20),
lower = 269128900001, upper = 269128937220)
vals <- makeVecCpp(lastPart, M, 37220)
write.csv(vals, paste0("myFile", 269128900001, ".csv", collapse = ""))
You will note that we avoid storing everything in memory by writing every one hundred thousand results to main storage hence the need for a huge hard drive. When I tested this, each file was about 2.5 Mb which would total to about 6.5 TB:
a <- 2.5 * (2^20) ### convert to bytes
a * (269128937220 / 1e5) / 2^40 ## get terabytes
[1] 6.416534
To give you an idea of how long this compuation will take, here is a timing for the first one hundred million results:
system.time(firstPart <- mclapply(seq(1, 100000000, 100000), function(x) {
negOne <- permuteGeneral(c(1L, -1L), freqs = c(21, 20),
lower = x, upper = x + 99999)
vals <- makeVecCpp(negOne, M, 100000)
write.csv(vals, paste0("myFile", x, ".csv", collapse = ""))
x
}, mc.cores = 8))
user system elapsed
529.931 9.557 80.690
80 seconds ain't that bad! That means we will only have to wait around for about 2.5 days!!!!!:
(269128937220 / 100000000 / 60 / 60 / 24) * 80
[1] 2.491935
If you really want to reduce this time, you will have to utilize a high performance computing service.
All results were obtained on a MacBook Pro 2.8GHz quad core (with 4 virtual cores.. 8 total).
RcppAlgospackage. In principle by doingpermuteGeneral(c(1,-1), freqs = c(21,20)). But they are too numerous:Error in CombinatoricsRcpp(v, m, repetition, freqs, lower, upper, constraintFun, : The number of rows cannot exceed 2^31 - 1.– Stéphane Laurent Jul 30 '18 at 13:332^31 - 1results you have utilize the argumentslowerandupperin order to generate permutations in chunks. – Joseph Wood Jul 30 '18 at 23:05