Scott Ritchie - 1 year ago 77

R Question

A recently asked question has lead me to believe the syntactic sugar for

`*`

`Rcpp`

Here's what we're trying to achieve in

`Rcpp`

`R`

`> m <- matrix(0:3, 2, 2)`

> m * 3

[,1] [,2]

[1,] 0 6

[2,] 3 9

I've created some minimal examples demonstrating both the problem above, and also some unexpected behaviour along the way. First note that I'm consistently using

`List`

`#include <Rcpp.h>`

using namespace Rcpp;

// [[Rcpp::export]]

List FooMat() {

// Create a fill a 2x2 matrix

NumericMatrix tmp(2,2);

for (int i = 0; i < 4; i++) {

tmp[i] = i;

}

return List::create(tmp);

}

// [[Rcpp::export]]

List FooMat2() {

// Create a fill a 2x2 matrix

NumericMatrix tmp(2,2);

for (int i = 0; i < 4; i++) {

tmp[i] = i;

}

NumericVector x(1);

x[1] = 3;

return List::create(tmp * x);

}

// [[Rcpp::export]]

List FooMat3() {

// Create a fill a 2x2 matrix

NumericMatrix tmp(2,2);

for (int i = 0; i < 4; i++) {

tmp[i] = i;

}

NumericVector x(1);

x[1] = 3;

return List::create(tmp * x[1]);

}

// [[Rcpp::export]]

List FooMat4() {

// Create a fill a 2x2 matrix

NumericMatrix tmp(2,2);

for (int i = 0; i < 4; i++) {

tmp[i] = i;

}

return List::create(tmp * 3);

}

Now if we source the file we get some odd behaviour:

`# Proof that we can return a NumericMatrix in a List:`

> FooMat()

[[1]]

[,1] [,2]

[1,] 0 2

[2,] 1 3

# Multiply the whole NumericMatrix by a whole NumericVector

# whose size is 1. Unsafe behaviour?

> FooMat2()

[[1]]

[1] 0.000000e+00 3.000000e+00 1.388988e-309 2.083483e-309

# Multiply the whole NumericMatrix by the first element of

# The NumericVector. Results are correct, but `*` converts

# the answer to a NumericVector instead of a NumericMatrix

> FooMat3()

[[1]]

[1] 0 3 6 9

# Same as FooMat3() except now we just multiply the NumericMatrix

# by an integer

> FooMat4()

[[1]]

[1] 0 3 6 9

One, the syntactic sugar for

`*`

`Rcpp`

`NumericVector`

`FooMat2()`

Answer Source

As I stated in previous answers, when I need to do actual math on matrices, I use Armadillo objects:

```
R> cppFunction('arma::mat scott(arma::mat x, double z) {
+ return(x*z); }',
+ depends="RcppArmadillo")
R> scott(matrix(1:4,2), 2)
[,1] [,2]
[1,] 2 6
[2,] 4 8
R>
```

Sugar operations are nice, but not complete. We will certainly take patches, though.

And as we said a few times before: rcpp-devel is the proper support channel.

**Edit** (Oct 2016 or 2 1/2 years later): Searching for something else just got me back here. In the Rcpp 0.12.* series, some work when into operations between matrix and vector so the basic 'matrix times scalar' now works as you'd expect:

```
R> cppFunction("NumericMatrix testmat(NumericMatrix m, double multme) {
+ NumericMatrix n = m * multme;
+ return n; }")
R> testmat(matrix(1:4,2), 1)
[,1] [,2]
[1,] 1 3
[2,] 2 4
R> testmat(matrix(1:4,2), 3)
[,1] [,2]
[1,] 3 9
[2,] 6 12
R>
```

I'd probably still use RcppArmadillo for math on matrices though.