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

Question

我做了一个简单的 Rcpp 函数来计算所有 pearson 相关系数，这些系数可以从输入矩阵 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

Answer 1

只是为了结束这个问题，我尝试运行仅分配输入的函数，而不是按照 cmets 中的建议运行实际算法，它返回就好了。我认为在 Windows 中，当输入达到一定大小时，窗口会变暗并在 R 控制台的窗口名称旁边显示“无响应”。但是，该功能仍在运行，因为如果留出足够的时间，它将最终完成，并且 R 控制台的窗口将恢复正常。进程花费了这么长时间，而且窗口看起来像 Rcpp 崩溃时的样子，这一事实让我认为进程没有运行，并且它是某种崩溃。

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