Understanding strange behavior in Rcpp

Question

I have run into something I cannot wrap my head around. It's part of a larger coding effort, but a minimal example is here:

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
Rcpp::List foo(arma::vec & tau2, const arma::vec & nu) {

  arma::vec bet = Rcpp::rnorm(3);
  tau2 = R::rgamma(1, arma::as_scalar(sum(pow(bet, 2)/nu)));

  return Rcpp::List::create(Rcpp::Named("nu") = nu,
                            Rcpp::Named("tau2") = tau2);
}

(tau2, although a scalar, is a vector here because I want to pass by reference: function pass by reference in RcppArmadillo)

What is puzzling me is that if I now run the following R code:

n <- 3
m <- matrix(0, n, 1)
for (r in 1:1000) {
  tau2 <- 1.0
  nu <- matrix(1, n, 1)
  upd <- foo(tau2, nu)
}

I get:

error: element-wise division: incompatible matrix dimensions: 3x1 and 18x1
Error in foo(tau2, nu) : 
  element-wise division: incompatible matrix dimensions: 3x1 and 18x1

where the 18x1 varies; mostly it's 0x1 but it's always a multiple of 3.

Looking at the output:

> nu
         [,1]     [,2]     [,3]     [,4]
[1,] 4.165242 4.165242 4.165242 4.165242
[2,] 4.165242 4.165242 4.165242 4.165242
[3,] 4.165242 4.165242 4.165242 4.165242
> upd
$nu
     [,1]
[1,]    1
[2,]    1
[3,]    1

$tau2
         [,1]
[1,] 4.165242

That is, despite declaring nu as a constant reference (which I do because I do not want it changed), it is altered. The value it is filled with is upd$tau2 (but why?).

Strangely, I can make the behavior go away by seemingly meaningless changes by:

• putting tau2 <- 1.0 or nu <- matrix(1, n, 1) (or both) outside of the loop
• removing the reference in the first argument (i.e. using arma::vec tau2)
• not dividing pow(bet, 2) by nu
• changing to nu <- rep(1, n)

Perhaps the most confusing part is that if I select the code chunk inside of the loop and repeatedly run it, it works(!). However, if I run the R code using the loop it crashes on the second iteration.

Because I seem to be able to fix the problem, I'm mostly interested in learning what is going on here. I suspect it's just a consequence of my lack of expertise in C++ and recklessness with various variable types, so knowing what is causing all of this would be very valuable.

Answer 1

Two fixes:

1. tau2 is a double (mimic @dirkeddelbuettel here)
2. Temporary variable for the NumericVector of length n being generated prior to saving into bet

Code:

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

// [[Rcpp::export]]
Rcpp::List foo(double tau2, const arma::vec & nu) {

  int n = nu.n_elem;

  Rcpp::NumericVector x = Rcpp::rnorm(n);
  arma::vec bet = arma::vec(x.begin(), n, true, false);

  tau2 = R::rgamma(1, arma::as_scalar(sum(pow(bet, 2) / nu)));

  return Rcpp::List::create(Rcpp::Named("nu") = nu,
                            Rcpp::Named("tau2") = tau2);
}

Test case:

n <- 3
m <- matrix(0, n, 1)
for (r in 1:1000) {
  tau2 <- 1.0
  nu <- matrix(1, n, 1)
  upd <- foo(tau2, nu)
}
upd
#> $nu
#>      [,1]
#> [1,]    1
#> [2,]    1
#> [3,]    1
#> 
#> $tau2
#> [1] 3.292889

Answer 2

If I change the interface to using a double it all works:

Code

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
Rcpp::List foo(double & tau2, const arma::vec & nu) {

  arma::vec bet = Rcpp::rnorm(3);
  tau2 = R::rgamma(1, arma::as_scalar(sum(pow(bet, 2)/nu)));

  return Rcpp::List::create(Rcpp::Named("nu") = nu,
                            Rcpp::Named("tau2") = tau2);
}

/*** R
n <- 3
m <- matrix(0, n, 1)
for (r in 1:1000) {
  tau2 <- 1.0
  nu <- matrix(1, n, 1)
  upd <- foo(tau2, nu)
}
*/

Demo

R> sourceCpp("/tmp/hejseb.cpp")

R> n <- 3

R> m <- matrix(0, n, 1)

R> for (r in 1:1000) {
+   tau2 <- 1.0
+   nu <- matrix(1, n, 1)
+   upd <- foo(tau2, nu)
+ }
R> upd
$nu
     [,1]
[1,]    1
[2,]    1
[3,]    1

$tau2
[1] 1.77314

R> 

I am not sure if those are the numbers you expected. I don't really have time to work through what you are trying to do.