为何RRHO分析的Rcpp实现比原生R代码运行速度慢？

我在实现RRHO分析时，分别编写了原生R版本的calculate_hyper_overlap函数和Rcpp版本的calculate_hyper_overlap_cpp函数，但benchmark测试显示Rcpp版本运行速度显著慢于原生R代码，现分析该性能差异的产生原因。

原生R代码 #

## Compute the overlaps between two *character* atomic vector:
hyper_test <- function(sample1, sample2, n) {
    count <- length(intersect(sample1, sample2))
    m <- length(sample1)
    k <- length(sample2)
    # under-enrichment
    if (count <= m * k / n) {
        sign <- -1L
        pvalue <- stats::phyper(
            q = count, m = m, n = n - m,
            k = k, lower.tail = TRUE, log.p = FALSE
        )
    } else {
        # over-enrichment
        sign <- 1L
        pvalue <- stats::phyper(
            q = count, m = m, n = n - m,
            k = k, lower.tail = FALSE, log.p = FALSE
        )
    }

    c(count = count, pvalue = pvalue, sign = sign)
}

calculate_hyper_overlap <- function(sample1, sample2, n, stepsize) {
    row_ids <- seq.int(stepsize, length(sample1), by = stepsize)
    col_ids <- seq.int(stepsize, length(sample2), by = stepsize)
    indexes <- expand.grid(
        row_ids = row_ids,
        col_ids = col_ids
    )
    overlaps <- apply(as.matrix(indexes), 1L, function(x) {
        hyper_test(
            sample1[seq_len(x[["row_ids"]])],
            sample2[seq_len(x[["col_ids"]])],
            n = n
        )
    }, simplify = FALSE)
    overlaps <- data.table::transpose(overlaps)
    number_of_obj <- length(row_ids)
    matrix_counts <- matrix(
        overlaps[[1L]],
        nrow = number_of_obj
    )
    matrix_pvals <- matrix(
        overlaps[[2L]],
        nrow = number_of_obj
    )
    matrix_signs <- matrix(
        overlaps[[3L]],
        nrow = number_of_obj
    )
    list(
        counts = matrix_counts,
        pvalue = matrix_pvals,
        signs = matrix_signs
    )
}

Rcpp代码 #

// [[Rcpp::export]]
List calculate_hyper_overlap_cpp(CharacterVector sample1, CharacterVector sample2, int n, int stepsize)
{
    int list1_len = floor((sample1.size() - stepsize) / stepsize) + 1;
    int list2_len = floor((sample2.size() - stepsize) / stepsize) + 1;
    IntegerMatrix counts(list1_len, list2_len);
    NumericMatrix pvalue(list1_len, list2_len);
    IntegerMatrix signs(list1_len, list2_len);
    for (int i = 0; i < list1_len; i++)
    {
        for (int j = 0; j < list2_len; j++)
        {
            CharacterVector list1 = sample1[Range(0, (i + 1) * stepsize - 1)];
            CharacterVector list2 = sample2[Range(0, (j + 1) * stepsize - 1)];
            int count = intersect(list1, list2).size();
            counts(i, j) = count;
            int m = list1.size(), k = list2.size();
            if (count <= m * k / n)
            // under-enrichment
            {
                pvalue(i, j) = R::phyper(count, m, n - m, k, true, false);
                signs(i, j) = -1;
            }
            else
            // over-enrichment
            {
                pvalue(i, j) = R::phyper(count, m, n - m, k, false, false);
                signs(i, j) = 1;
            }
        }
    }
    return List::create(Named("counts") = counts,
                        Named("pvalue") = pvalue,
                        Named("signs") = signs);
}

测试代码 #

n <- 200
sample1 <- rnorm(n)
sample2 <- rnorm(n)
names(sample1) <- names(sample2) <- paste0("gene", seq_len(n))
bench_res <- bench::mark(
    res1 <- calculate_hyper_overlap_cpp(
        names(sample1), names(sample2),
        n = n, stepsize = 3L
    ),
    res2 <- calculate_hyper_overlap(
        names(sample1), names(sample2),
        n = n, stepsize = 3L
    ),
    check = FALSE
)
dplyr::select(bench_res, where(~ !is.list(.x)))

性能差异原因分析 #

• 频繁临时对象创建与内存拷贝：Rcpp代码的嵌套循环中，每次迭代都会通过Range创建新的CharacterVector子向量，这会触发大量内存分配与数据拷贝操作。而原生R的子集操作依赖内部优化（如延迟拷贝、内存共享），实际拷贝开销远低于显式的Rcpp对象创建。
• intersect函数的实现效率差异：Rcpp中调用的intersect本质是调用R的基础函数，并非原生C++的高效集合实现；而原生R的intersect经过高度优化，底层采用哈希表或排序后快速查找的方式，处理字符向量交集的效率远高于Rcpp版本。
• 循环交互的额外开销：Rcpp的嵌套循环中，每次迭代都要在C环境与R对象系统之间切换，频繁的对象交互会累积大量开销；而原生R的apply循环经过内部优化，批量处理的效率反而优于显式的C循环（尤其是当循环次数较多时）。

内容的提问来源于stack exchange，提问作者Yun