与Armadillo与`copy_aux_mem`与三角矩阵求解的奇怪/不一致的行为

发布于 2025-02-05 14:54:33 字数 1681 浏览 2 评论 0 原文

考虑以下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（），因为这是我设置示例的最简单方法。使用3D矩阵，我

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＆lt; -4L 一起使用。但是，使用5D矩阵，我得到

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。我发布此信息的原因是由于单位测试升级到新版本的Rcpparmadillo之后，该单元测试失败了。我获得了0.10.8.1.0版本的一致且预期的行为。

如果我删除 arma :: trimatu（）调用，则差异仍然存在。

原文

Consider the following C++ code

// [[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());
}

I am using Rcpp::sourceCpp() as this is the easiest way for me to set the example up. With a 3D matrix I get

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

This exactly what I will expect (the input argument is changed as it is used for the memory for the result). It also works with n <- 4L. However with a 5D matrix I get

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"

This is unlike what I would expect (the input argument is not changed).

Are my expectations wrong?

The above is with RcppArmadillo version 0.11.1.1.0. The reason I am posting this is because of a unit test which failed after I upgraded to the new version of RcppArmadillo. I get a consistent and expected behavior with version 0.10.8.1.0.

The difference persists if I remove the arma::trimatu() calls.