为什么我的 Rcpp 函数在使用大型输入矩阵时会崩溃?

Why is my Rcpp function crashing when using a large input matrix?

我做了一个简单的 Rcpp 函数来计算可以从输入矩阵的所有行组合计算出的所有皮尔逊相关系数 E。结果以 4 位小数精度(整数格式)存储在向量 v 中。如果 E 的维度不是太大,则该函数可以正常工作,但当我使用与我要使用该函数处理的真实数据相似的数据大小时进行测试时,该函数就会崩溃。

这是 Rccp 代码:

#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
void pearson(NumericMatrix E, IntegerVector v){
    int rows = E.nrow();
    int cols = E.ncol();
    int j, irow, jrow;
    double rowsum;
    NumericVector means(rows);
    int k = 0;
    double cov, varx, vary;
    double pearson;

    for(irow = 0; irow < rows; irow++){
        rowsum = 0;
        for(j = 0; j < cols; j++){
            rowsum += E(irow, j);
        }
        means[irow] = rowsum / cols;
    }
    
    for(irow = 0; irow < rows - 1; irow++){
        for(jrow = irow + 1; jrow < rows; jrow++){
            cov = 0;
            varx = 0;
            vary = 0;
            for(j = 0; j < cols; j++) {
                cov += (E(irow, j) - means[irow]) * (E(jrow, j) - means[jrow]);
                varx += std::pow(E(irow, j) - means[irow], 2);
                vary += std::pow(E(jrow, j) - means[jrow], 2);
            }
            pearson = cov / std::sqrt(varx * vary);
            v[k] = (int) (pearson * 10000);
            k++;
        }
    }

}

然后为了在 R 中测试它,我从以下开始:

library(Rcpp)
sourceCpp("pearson.cpp")
testin <- matrix(rnorm(1000 * 1100), nrow = 1000, ncol = 1100)
testout <- integer( (nrow(testin) * (nrow(testin) - 1)) / 2 )
pearson(testin, testout) # success!

然而,当增加输入大小时,R 会话在执行以下脚本的最后一行后崩溃:

library(Rcpp)
sourceCpp("pearson.cpp")
testin <- matrix(rnorm(16000 * 17000), nrow = 16000, ncol = 17000)
testout <- integer( (nrow(testin) * (nrow(testin) - 1)) / 2 )
pearson(testin, testout) # sad

我觉得这很奇怪,因为我可以在执行函数之前很好地分配输入和输出。在函数内部,输出向量通过引用进行修改。无法弄清楚出了什么问题。目前我正在使用 16GB RAM 的机器。

编辑:sessionInfo()

的输出
R version 4.0.4 (2021-02-15)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

locale:
[1] LC_COLLATE=Spanish_Mexico.1252 
[2] LC_CTYPE=Spanish_Mexico.1252   
[3] LC_MONETARY=Spanish_Mexico.1252
[4] LC_NUMERIC=C                   
[5] LC_TIME=Spanish_Mexico.1252    

attached base packages:
[1] stats     graphics  grDevices
[4] utils     datasets  methods  
[7] base     

other attached packages:
[1] Rcpp_1.0.5

loaded via a namespace (and not attached):
[1] compiler_4.0.4

只是为了结束这个问题,我尝试 运行 函数只是分配输入,而不是 运行 评论中建议的实际算法 return很好。我认为在 Windows 中,由于某种原因,当输入达到一定大小时,window 会变暗并在 R 控制台的 window 名称旁边显示“无响应”。然而,该功能仍然是 运行ning,因为如果有足够的时间它最终会完成并且 R 控制台的 window 将 return 恢复正常。这个过程花了这么长时间,而且 window 看起来像 Rcpp 崩溃时的样子,这让我认为这个过程不是 运行ning,而是某种崩溃。

我最终做的是借助 RcppParallel 的一些创建者非常有用的 tutorial 编写该算法的并行版本。由于内存限制,我无法使用基本 R cor() 函数,因此并行版本非常适合我的需要。