Bizarre behavior of Rcpp NumericMatrix when multiplied by scalar on left/right

Question

Can anyone explain the following behavior?

When declaring a new NumericMatrix, y, as the original matrix, x, multiplied by a scalar, c, the order of the scalar/matrix multiplication matters. If I multiply the scalar on the left and the matrix on the right (e.g. NumericMatrix y = c * x;), then I get bizarre behavior. The original matrix, x, is changed!

However if I put the original matrix to the left and multiply the scalar on the right (e.g. NumericMatrix y = x * c;), x remains unchanged.

This doesn't seem to affect other data types. I've tested with int and NumericVector.

Example: Problem appears when using `NumericMatrix`

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix testfun(NumericMatrix x) {

  NumericMatrix y(x.rows(), x.cols());

  y = x * 2;

  std::cout << x; // x is unmodified

  y = 2 * x; 

  std::cout << x; // x is now modified

  return x;
}



/*** R
x <- matrix(2, nrow = 3, ncol = 3)

print(x)

y <- testfun(x = x)

print(y)

print(x)

*/

The output is as follows.

> x <- matrix(2, nrow = 3, ncol = 3)

> print(x)
     [,1] [,2] [,3]
[1,]    2    2    2
[2,]    2    2    2
[3,]    2    2    2

> y <- testfun(x = x)
2.00000 2.00000 2.00000
2.00000 2.00000 2.00000
2.00000 2.00000 2.00000
4.00000 4.00000 4.00000
4.00000 4.00000 4.00000
4.00000 4.00000 4.00000

> print(y)
     [,1] [,2] [,3]
[1,]    4    4    4
[2,]    4    4    4
[3,]    4    4    4

> print(x)
     [,1] [,2] [,3]
[1,]    4    4    4
[2,]    4    4    4
[3,]    4    4    4

Here is my session info

> sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Catalina 10.15.1

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
[1] compiler_3.6.1            tools_3.6.1               RcppArmadillo_0.9.800.1.0
[4] Rcpp_1.0.2                RcppProgress_0.4.1        packrat_0.5.0            
[7] RcppParallel_4.4.4

.

I don't see how "pass by reference vs. pass by value" could be the issue. If that was the case, I would expect to see the same behavior with `int` and `NumericVector`. Why would this only affect `NumericMatrix`? — Tommy Jones, Dec 28 '19 at 23:55
To make the point more succinctly, using `Rcpp::clone` as suggested in the linked 'answer' still appears to modify `x` as a side-effect when `clone` is passed an expression: `auto z = Rcpp::clone(2 * x);` modifies `x` when `x` is a `Rcpp::NumericMatrix`. So I agree that answer is _not_ satisfactory at face value. — smarchese, Dec 29 '19 at 00:48
@TommyJones While you're not right re: `int`, I would also expect `Rcpp::NumericVector`; I don't think the existing discussions (at least not the ones I've seen) on `Rcpp` classes and pass by reference vs. value are sufficient to explain the behavior observed with a `NumericVector` and how it contrasts with the observed behavior of a `NumericMatrix`. Maybe I'm missing something, but I think your question needs a new answer. It may be helpful if you can edit your question with a very succinct blurb at the top explaining why it's different than the marked duplicate (which should be easy for you) — duckmayr, Dec 29 '19 at 02:25
It may also help if you can get your point across in a third or quarter of the space. Make it more compact and concise -- you are desiring to get other people's attention for free, so make it _as easy and obvious_ as possible for them to give you the help. — Dirk Eddelbuettel, Dec 29 '19 at 14:15
I'm sorry that I wasn't more succinct. In my defense, the issue was stated right at the top and the first example illustrates it. Also _my_ issue is solved. (I reversed order of the multiplication.) I made this post to help others who may come across this issue. I spent a good chunk of yesterday searching the issues on Rcpp's GitHub page, mailing list, and stack overflow to ensure that it _hadn't_ been answered before posting a question here. And I appreciate that you're doing this for free, as am I. I spent hours yesterday answering user questions about my own package. :( — Tommy Jones, Dec 29 '19 at 16:08

Dirk Eddelbuettel · Accepted Answer · 2019-12-30T22:02:54.527

4

The question is way too long as posted and hides its point. We do not need the second and third examples. All we need is this code:

Code

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix testfun(NumericMatrix x) {
  NumericMatrix y(x.rows(), x.cols());
  y = x * 2;
  std::cout << x; // x is unmodified
  y = 2 * x;      // contrast with x * 2
  std::cout << x; // x is now modified
  return x;
}

/*** R
print(x <- matrix(2, nrow = 2, ncol = 2))
print(y <- testfun(x = x))
print(x)
*/

Output

R> Rcpp::sourceCpp("~/git/stackoverflow/59515517/question.cpp")

R> print(x <- matrix(2, nrow = 2, ncol = 2))
     [,1] [,2]
[1,]    2    2
[2,]    2    2

R> print(y <- testfun(x = x))
2.00000 2.00000
2.00000 2.00000
4.00000 4.00000
4.00000 4.00000
     [,1] [,2]
[1,]    4    4
[2,]    4    4

R> print(x)
     [,1] [,2]
[1,]    4    4
[2,]    4    4
R>

Issue

Behaviour of x * 2 is different from 2 * x. The latter has a side effect. This is likely a bug,

Bigger issue

There is not really a point to do matrix algebra with Rcpp. The implementation is rudimentary and incomplete. If you want to do "math", use RcppArmadillo or RcppEigen.

Arma implementation

#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
arma::mat testfun(arma::mat x) {
  arma::mat y(x.n_rows, x.n_cols);
  y = x * 2;
  std::cout << x; // x is unmodified
  y = 2 * x;      // contrast with x * 2
  std::cout << x; // x is now modified
  return x;
}

/*** R
print(x <- matrix(2, nrow = 2, ncol = 2))
print(y <- testfun(x = x))
print(x)
*/

Arma output

R> Rcpp::sourceCpp("~/git/stackoverflow/59515517/answer.cpp")

R> print(x <- matrix(2, nrow = 2, ncol = 2))
     [,1] [,2]
[1,]    2    2
[2,]    2    2

R> print(y <- testfun(x = x))
   2.0000   2.0000
   2.0000   2.0000
   2.0000   2.0000
   2.0000   2.0000
     [,1] [,2]
[1,]    2    2
[2,]    2    2

R> print(x)
     [,1] [,2]
[1,]    2    2
[2,]    2    2
R>

I'll see if I can get to the Rcpp bug you found here.

Edit: It is a bug, and I don't quite understand it yet, but I added something to the issue where the code originated: https://github.com/RcppCore/Rcpp/issues/365

Edit 2: The fix is now in master. Thanks to KK for the PR, to everbody else hinting in the comments what may be behind this, and especially to Ralf for getting us all to try the fix.

edited Dec 30 '19 at 22:02

answered Dec 29 '19 at 14:27

Dirk Eddelbuettel

360,940
56
644
725

3

Yeah, when I saw this, I looked at `Matrix.h` also and couldn't figure out why the order (which was `lhs` and which was `rhs`) mattered either. Good luck... it is bizarre indeed! Solid advice on using Armadillo instead – duckmayr Dec 29 '19 at 19:16
2

Ralf suggested a fix which I had also contemplated but not tried. So this may be squashed "soon". – Dirk Eddelbuettel Dec 29 '19 at 20:21
2

I saw that same asymmetry, but didn't try anything with it because I couldn't reason out why it'd matter. Awesome that it seems to work though! – duckmayr Dec 29 '19 at 20:22
3

Ditto! I have no idea why that change matters. Isn't C++ fun? ;-) – Dirk Eddelbuettel Dec 29 '19 at 20:24
3

I have no idea why it should matter either. It was more the physicist in me striving for symmetry ;-) – Ralf Stubner Dec 29 '19 at 20:28
1

Thank you for looking into this, all. And thank you for making Rcpp an awesome tool! – Tommy Jones Dec 29 '19 at 21:10
2

Feel free to accept the answer again. It did lead us to a fix :) – Dirk Eddelbuettel Dec 29 '19 at 22:24
@DirkEddelbuettel in the Arma implementation, this comment is wrong: `std::cout << x; // x is now modified`. The Arma output shows that `x` is not modified. Armadillo doesn't have the Rcpp bug. Copy paste error? – mtall Jan 04 '20 at 04:03
What file and line? – Dirk Eddelbuettel Jan 04 '20 at 04:14
@DirkEddelbuettel In your answer above. Code under the title **Arma implementation**. Line 11. – mtall Jan 04 '20 at 12:52
I see. Copy and past effect from the original answer. The two go together, and I just changed the typed. The result matters, and it is correct. – Dirk Eddelbuettel Jan 04 '20 at 12:56
Even though the result is correct, the comment in the code is wrong and misleading. `x` is not modified. – mtall Jan 04 '20 at 13:25