首先，您最好重新考虑您的方法。话虽如此，让我们开始解决您的问题。
这是一个非常困难的问题，主要是由于资源的限制。下面，我有一个解决方案，如果您可以访问相当数量的存储空间（至少7 TB），它将在家用计算机上在合理的时间内完成。下面的算法不需要那么多内存，可以调整以减少内存使用。
在我们开始之前，我们注意到一开始仅仅生成这么多排列似乎是不可能的。然而，在高度优化的C++ 代码和并行计算的帮助下，这项任务又回到了可能的领域。这在我的answer 对 OP 的上一个问题中得到了证明。我们利用RcppAlgos（我是作者）和parallel 包在使用8 个核心的100 万块中每秒生成约3600 万个排列。
现在，我们负责尽可能快地对每个排列进行特定计算。计算如下：
(SLS') / 4, where S is a permutation, L is a 41 x 41 matrix
这里有几个base R 方法（注意m1[x, ] %*% m2 %*% m1[x, ] 与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
}
让我们稍微考虑一下。我们有一堆数字一和负一的排列。当我们将这些排列乘以实数矩阵时，例如 M，我们最终只是简单地从 M 中添加和减去值。我敢打赌，我们可以使用 Rcpp 加快这一速度，并避免浪费（和无用）的身份乘法（即乘以 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;
}
现在让我们看看它们是否给出了相同的结果并对其进行基准测试：
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
我们的Rcpp 实现显然是赢家！！接下来，我们将其纳入我们的最终答案：
## 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 = ""))
您会注意到，我们通过将每十万个结果写入主存储器来避免将所有内容存储在内存中，因此需要一个巨大的硬盘驱动器。当我对此进行测试时，每个文件大约是2.5 Mb，总计大约是6.5 TB：
a <- 2.5 * (2^20) ### convert to bytes
a * (269128937220 / 1e5) / 2^40 ## get terabytes
[1] 6.416534
为了让您了解此计算需要多长时间，以下是前一亿个结果的时间安排：
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 秒还不错！这意味着我们只需要等待大约 2.5 天！！！！！！：
(269128937220 / 100000000 / 60 / 60 / 24) * 80
[1] 2.491935
如果你真的想减少这个时间，你将不得不使用高性能计算服务。
所有结果均在 MacBook Pro 2.8GHz 四核（4 个虚拟核心......总共 8 个）上获得。