Armadillo使用copy_aux_mem与三角矩阵求解时的不一致行为问题

开发者阿泽

2026-9-1

问题复现

考虑如下C++代码：

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

// [[Rcpp::export(rng = false)]]
void possible_bug(arma::vec &x, arma::mat const &sig_chol){
  if(x.n_elem != sig_chol.n_rows * sig_chol.n_cols)
    throw std::runtime_error("boh");
  
  arma::mat x_mat(x.begin(), sig_chol.n_rows, sig_chol.n_rows, false);
  x_mat = arma::solve
    (arma::trimatu(sig_chol),
     arma::solve(arma::trimatu(sig_chol), x_mat).t());
}

使用Rcpp::sourceCpp()可以最简便地搭建复现示例。当矩阵阶数n=3时，运行如下R测试代码得到的结果为：

set.seed(1)
n <- 3L
x <- rnorm(n * n)
sig_chol <- rWishart(1, n, diag(n)) |> drop() |> chol()

x_cp <- x + 0.
possible_bug(x = x_cp, sig_chol)
all.equal(x_cp, x)
#R> [1] "Mean relative difference: 1.000271"
all.equal(x_cp, 
          solve(sig_chol, t(solve(sig_chol, matrix(x, n)))) |> c())
#R> [1] TRUE

该结果符合预期：输入参数被修改，其内存被复用来存储计算结果，n=4时行为也与预期一致。但当矩阵阶数n=5时，运行如下测试代码得到的结果为：

set.seed(1)
n <- 5L
x <- rnorm(n * n)
sig_chol <- rWishart(1, n, diag(n)) |> drop() |> chol()

x_cp <- x + 0.
possible_bug(x = x_cp, sig_chol)
all.equal(x_cp, x)
#R> [1] TRUE
all.equal(x_cp, 
          solve(sig_chol, t(solve(sig_chol, matrix(x, n)))) |> c())
#R> [1] "Mean relative difference: 2.041895"

该结果不符合预期：输入参数未被修改。上述测试基于RcppArmadillo 0.11.1.1.0版本，升级至该新版本后单元测试失败，旧版本0.10.8.1.0中可得到一致且符合预期的行为。
移除arma::trimatu()调用后，上述不一致行为仍然存在。

原因说明

对该行为的预期本身不符合Armadillo的API约定，旧版本的正常表现只是实现层面的巧合，新版本的行为变化直接暴露了写法的问题：

构造x_mat时传入false作为copy_aux_mem参数，仅代表Armadillo可以将你提供的外部内存作为矩阵的初始存储，但从来没有承诺对该矩阵的赋值操作一定会写入这块外部内存。这个构造方式的设计初衷是让Armadillo直接读取外部内存避免拷贝，只有在矩阵尺寸不变、赋值操作未触发别名安全保护的特殊情况下，才会直接写入外部内存。
赋值语句右侧是嵌套solve+转置的表达式模板，本身存在读写重叠问题：需要先读取x_mat做第一次solve，转置后做第二次solve，最后把结果写回x_mat。Armadillo的别名安全机制会在检测到这种读写重叠时，先把整个右侧表达式求值到内部临时矩阵，避免计算过程中覆盖还未读取的原始数据导致结果错误。
0.10.x旧版本的别名检测阈值较高，n=3、n=4的小矩阵走了直接写入目标内存的计算路径，刚好没有触发临时对象逻辑，所以能看到x被修改；0.11.x版本调整了别名检测逻辑，n≥5的矩阵会触发临时对象求值，计算结果先存在Armadillo内部分配的临时内存中。此时对于copy_aux_mem=false构造的矩阵，新版本的赋值逻辑会直接将x_mat的内部内存指针指向临时矩阵的内存，而不是把结果拷贝到你提供的外部内存，等函数退出时x_mat和临时对象一同销毁，传入的x内存从头到尾没有被写入，自然表现为x_cp和原始x完全一致。
移除arma::trimatu()后问题仍然存在，也能印证问题和三角矩阵标记无关，核心是别名检测逻辑+外部内存矩阵的赋值规则变化。
额外说明：哪怕在旧版本，原有写法也存在计算错误风险——如果第一次solve计算时直接写入x_mat内存，会在计算未完成时覆盖原始x_mat值，导致第二次solve的输入是半计算半原始的错误数据，只是n较小时计算顺序刚好没有触发这个问题。

修复方案

不要依赖copy_aux_mem=false构造的矩阵做隐式内存覆写，两种稳妥的修改方式：

显式求值后手动拷贝到外部内存，保留原有修改入参的逻辑：

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

// [[Rcpp::export(rng = false)]]
void possible_bug_fixed(arma::vec &x, arma::mat const &sig_chol){
  if(x.n_elem != sig_chol.n_rows * sig_chol.n_cols)
    throw std::runtime_error("boh");
  const unsigned int n = sig_chol.n_rows;
  arma::mat x_mat(x.begin(), n, n, false);
  // 先把右侧表达式显式求值为临时矩阵
  arma::mat res = arma::solve(
    arma::trimatu(sig_chol),
    arma::solve(arma::trimatu(sig_chol), x_mat).t()
  );
  // 手动将结果拷贝到x_mat绑定的外部内存，避免指针重绑定
  std::memcpy(x_mat.memptr(), res.memptr(), sizeof(double) * n * n);
}

按值传参、直接返回结果，从根源上避免内存别名问题，这也是Rcpp更推荐的写法：

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

// [[Rcpp::export(rng = false)]]
arma::vec possible_bug_fixed(arma::vec x, arma::mat const &sig_chol){
  if(x.n_elem != sig_chol.n_rows * sig_chol.n_cols)
    throw std::runtime_error("boh");
  const unsigned int n = sig_chol.n_rows;
  arma::mat x_mat(x.begin(), n, n, false);
  x_mat = arma::solve(
    arma::trimatu(sig_chol),
    arma::solve(arma::trimatu(sig_chol), x_mat).t()
  );
  return x;
}

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