c++armadillorcpparmadillo

Strange/inconsistent behavior with armadillo with `copy_aux_mem` and solving with a triangular matrix

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.

Solution

  • This is an intended change in version 11.1.1. See https://gitlab.com/conradsnicta/armadillo-code/-/issues/210#note_974299524

    The way to go is to set strict to true in the constructor of arma::mat and other objects.