R循环内调用Rcpp函数崩溃，单独调用正常问题排查

问题描述 #

编写了包含StrausProcess_Unc和Birth_Death_Sim_hard_sample_2两个导出函数的Rcpp代码，编译无报错，单独重复调用Birth_Death_Sim_hard_sample_2也正常，但在R循环中多次调用该函数时程序总会崩溃。

R复现代码 #

sz = c()
for(i in 1:20000){
  sz[i] = length( Birth_Death_Sim_hard_sample_2(50, c(0.5, 5.0, 10.0, 250.0, 500.0), 0.5, 0.5,
                                                1e-300, 35, 5, 2400, 0,
                                                242, 2) )
}

Rcpp代码 #

#include <RcppArmadilloExtensions/sample.h>
#include <random>
#include <iostream>
#include <cmath>
#include <vector>
#include <numeric>
#include <algorithm>
#include <stdio.h>    
#include <math.h>
#include<Rmath.h>
using namespace Rcpp;

// [[Rcpp::export]]
double StrausProcess_Unc(NumericVector xi, double MX, double MN, double a, double d, double b){
  if (sum(xi > MX) == 0 & sum(xi < MN) == 0) {
    int K_dot = xi.size();
    double term_2 = 0;
    if (K_dot > 1) {
      for (int i = 0; i < (K_dot - 1); ++i) {
        for (int j = (i + 1); j < K_dot; ++j) {
          term_2 += (abs(xi[i] - xi[j]) <= d) * log(a);
        }
      }
    }
    return K_dot * log(b) + term_2;
  } else {
    for (int i = 0; i < xi.size(); ++i) {
      xi[i] = std::max(std::min(xi[i], MX), MN);
    }
    int K_dot = xi.size();
    double term_2 = 0;
    if (K_dot > 1) {
      for (int i = 0; i < (K_dot - 1); ++i) {
        for (int j = (i + 1); j < K_dot; ++j) {
          term_2 += (abs(xi[i] - xi[j]) <= d) * log(a);
        }
      }
    }
    return K_dot * log(b) + term_2;
  }
}

// [[Rcpp::export]]
NumericVector Birth_Death_Sim_hard_sample_2(int niter, NumericVector mu_init, double q_move, double q_birth,
                                          double a, double d, double beta, double MX, double MN,
                                          double Upper, double Lower){
  double num;
  double den;
  int old_n;
  int n = mu_init.size();
  NumericVector mu_cand(n);
  NumericVector mu_old(n);
  mu_old = clone(mu_init);
  double m_birth;
  NumericVector m_birth_1(2); //vector to be used only when we reach length 1
  int j_star; //position for birth sample
  int j_star_1; //position of component to be killed
  for(int iter=1; iter<niter; ++iter){
    //move change current means
    old_n = mu_old.size();
    if(R::runif(0,1) <= q_move){
      for(int k=0; k<old_n; ++k){
        mu_cand[k] = R::rlnorm(log(mu_old[k]),0.1);
      }
      num = StrausProcess_Unc(mu_cand,  MX,  MN,  a,  d,  beta);

      den = StrausProcess_Unc(mu_old,  MX,  MN,  a,  d,  beta);
      if(R::runif(0,1)<= exp(num-den)){
        mu_old = clone(mu_cand);
      }else{
        mu_old = clone(mu_old);
      }
    }else{
      if(old_n==1){// give birth
        m_birth = R::runif(Lower,Upper);//R::rlnorm(mu_0,sigma_0);
        m_birth_1[0] = mu_old[0];
        m_birth_1[1] = m_birth;
        num = StrausProcess_Unc(m_birth_1,  MX,  MN,  a,  d,  beta);
        den = StrausProcess_Unc(mu_old,  MX,  MN,  a,  d,  beta)  + R::dunif(m_birth,Lower,Upper,1);
        if(R::runif(0,1) <= exp(num-den)){
          mu_old = clone(m_birth_1);
        }else{
          mu_old = clone(mu_old);
        }
      }else{
         if(R::runif(0,1) <= q_birth & old_n<100){//give birth
          j_star = RandInt(old_n+1);
          NumericVector mu_cand_i(old_n+1);
          mu_cand_i[j_star] = R::runif(Lower,Upper);//R::rlnorm(mu_0,sigma_0);
          if(j_star==0){
            Range r(j_star + 1, old_n);
            mu_cand_i[r] = clone(mu_old);
          }else if(j_star== old_n){
            Range r(0, old_n - 1);
            mu_cand_i[r] = clone(mu_old);
          }else{
             Range r(0, j_star-1);
             Range u(j_star+1,old_n);
             Range x(j_star,old_n-1);
             mu_cand_i[r] = mu_old[r];
             mu_cand_i[u] =  mu_old[x];
          }
          num = StrausProcess_Unc(mu_cand_i,  MX,  MN,  a,  d,  beta);
          den = StrausProcess_Unc(mu_old,  MX,  MN,  a,  d,  beta) + R::dunif(mu_cand_i[j_star],Lower,Upper,1);
          if(R::runif(0,1) <= exp(num-den)){
            mu_old = clone(mu_cand_i);
          }else{
            mu_old = clone(mu_old);
          }
        }else{//give death
          j_star_1 = RandInt(old_n);
          NumericVector mu_cand_j(old_n-1);
          if(j_star_1==0){
            Range r(j_star_1+1, old_n-1);
            mu_cand_j = mu_old[r];
          }else if(j_star_1 == (old_n-1) ){
            Range r(0, old_n-2);
            mu_cand_j = mu_old[r];
          }else{
             Range r(0, j_star_1-1);
             Range u(j_star_1+1,old_n-1);
             Range x(j_star_1,old_n-2);
             mu_cand_j[r] = mu_old[r];
             mu_cand_j[x] =  mu_old[u];
          }
          num = StrausProcess_Unc(mu_cand_j,  MX,  MN,  a,  d,  beta) + R::dunif(mu_old[j_star_1],Lower,Upper,1);
          den = StrausProcess_Unc(mu_old,  MX,  MN,  a,  d,  beta);
          if(R::runif(0,1)<=exp(num-den)){
            mu_old = clone(mu_cand_j);
          }else{
            mu_old = clone(mu_old);
          }
        }
      }
    }
  }
  return(mu_old);
}

错误分析与修正 #

1. 核心崩溃原因：内存越界 #

mu_cand初始化为与mu_init相同的大小，但循环中mu_old的长度会随birth/death操作动态变化。当mu_old长度超过初始n时，for(int k=0; k<old_n; ++k)会访问mu_cand[k]超出其长度的位置，直接触发内存越界，多次调用后必然导致程序崩溃。

修正：将mu_cand的定义移到move操作块内部，每次根据当前old_n动态创建对应大小的向量：

if(R::runif(0,1) <= q_move){
  NumericVector mu_cand(old_n); // 动态匹配当前mu_old的长度
  for(int k=0; k<old_n; ++k){
    mu_cand[k] = R::rlnorm(log(mu_old[k]),0.1);
  }
  num = StrausProcess_Unc(mu_cand,  MX,  MN,  a,  d,  beta);
  den = StrausProcess_Unc(mu_old,  MX,  MN,  a,  d,  beta);
  if(R::runif(0,1)<= exp(num-den)){
    mu_old = clone(mu_cand);
  }
}

2. 意外修改输入向量 #

StrausProcess_Unc函数中直接修改传入的xi向量（xi[i] = std::max(std::min(xi[i], MX), MN);），而Rcpp的NumericVector默认按引用传递，这会意外修改调用者的向量（如mu_cand、mu_old），导致后续计算逻辑混乱。

修正：在函数内部克隆输入向量，避免修改原数据：

double StrausProcess_Unc(NumericVector xi_in, double MX, double MN, double a, double d, double b){
  NumericVector xi = clone(xi_in); // 克隆输入，隔离修改
  // 后续逻辑不变
}

3. 不必要的内存操作 #

代码中多处出现mu_old = clone(mu_old);，完全是冗余操作，只会浪费内存和计算资源，直接删除即可。

4. 数值稳定性优化 #

当a=1e-300时，log(a)会得到极小的负数（约-690.7755），重复累加可能导致数值下溢。可以提前计算一次log(a)并复用，减少重复计算：

double log_a = log(a);
for (int i = 0; i < (K_dot - 1); ++i) {
  for (int j = (i + 1); j < K_dot; ++j) {
    if (abs(xi[i] - xi[j]) <= d) {
      term_2 += log_a;
    }
  }
}

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