RcppArmadillo中each_slice+Lambda实现Cholesky分解失效问题求助

你的问题出在对Armadillo each_slice 函数的用法理解上，咱们一步步拆解：

为什么chol_array返回原矩阵？ #

1. 传值参数导致操作副本：你写的Lambda表达式里用了arma::mat X作为参数，这是传值传递——也就是说，Lambda拿到的是每个slice的副本，你在里面计算arma::chol(X)并返回，根本不会修改原Sigma_chol里的slice。
2. 未利用each_slice的返回值：Armadillo的each_slice如果传入的Lambda有返回值（比如返回一个arma::mat），它会生成并返回一个新的cube，而不是修改调用它的原cube。你现在只是调用了each_slice，但没把返回值赋值给Sigma_chol，所以原Sigma_chol还是你初始复制的原矩阵。

两种修正方法 #

方法1：用引用参数直接修改现有slice

把Lambda的参数改成引用类型，这样就能直接修改Sigma_chol里的每个slice：

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
// [[Rcpp::export]]
arma::cube chol_array_fixed1(arma::cube Sigma) {
  arma::cube Sigma_chol = Sigma;
  // 使用引用参数& X，直接修改原slice
  Sigma_chol.each_slice([](arma::mat& X) { 
    X = arma::chol(X); 
  });
  return Sigma_chol;
}

方法2：利用each_slice的返回值生成新cube

不需要提前复制原矩阵，直接让each_slice返回处理后的新cube：

// [[Rcpp::export]]
arma::cube chol_array_fixed2(arma::cube Sigma) {
  // each_slice会遍历每个slice，应用Lambda生成新矩阵，最终返回完整的新cube
  return Sigma.each_slice([](const arma::mat& X) { 
    return arma::chol(X); 
  });
}

验证效果 #

在R里测试这两个修正版本，结果会和chol_array2完全一致：

Sigma <- array(crossprod(matrix(rnorm(9), 3, 3)), dim = c(3, 3, 2))
chol_array_fixed1(Sigma)
chol_array_fixed2(Sigma)
chol_array2(Sigma)

关键总结 #

• 如果要修改现有cube的slice，Lambda必须用引用参数（arma::mat& X），否则操作的是副本，不会生效。
• 如果要生成新cube，直接利用each_slice的返回值，Lambda返回处理后的矩阵即可，代码更简洁。

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