如何在Rcpp/RcppArmadillo中高效置换矩阵行以加速平稳自助法？

安全阿岳

2026-7-11

优化平稳自助法的RcppArmadillo实现方案

核心优化思路

跳过显式构造置换后矩阵的步骤，直接通过随机索引对原矩阵行进行视图访问，逐次计算列和，大幅减少内存拷贝与冗余计算；同时结合预分配内存、并行计算进一步压缩耗时。

基础高效实现代码

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
arma::mat bootstrap_colsum_fast(const arma::mat& X, int B) {
    int N = X.n_rows;
    int M = X.n_cols;
    
    // 预分配结果矩阵：M行×B列，存储B次置换后的列和
    arma::mat col_sums(M, B, arma::fill::zeros);
    
    for (int b = 0; b < B; ++b) {
        // 生成0-based随机置换索引
        arma::uvec idx = arma::shuffle(arma::linspace<arma::uvec>(0, N-1, N));
        // 通过行索引视图直接计算列和，无数据拷贝
        col_sums.col(b) = arma::sum(X.rows(idx), 0).t();
    }
    
    return col_sums;
}

多线程加速版本

利用OpenMP实现并行循环，进一步提升大B值场景下的效率：

#include <RcppArmadillo.h>
#include <omp.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::plugins(openmp)]]

// [[Rcpp::export]]
arma::mat bootstrap_colsum_omp(const arma::mat& X, int B, int n_threads = 4) {
    int N = X.n_rows;
    int M = X.n_cols;
    
    arma::mat col_sums(M, B, arma::fill::zeros);
    
    omp_set_num_threads(n_threads);
    
    #pragma omp parallel for
    for (int b = 0; b < B; ++b) {
        arma::uvec idx = arma::shuffle(arma::linspace<arma::uvec>(0, N-1, N));
        col_sums.col(b) = arma::sum(X.rows(idx), 0).t();
    }
    
    return col_sums;
}