OK, 'sweep' is a fair bit (15x) slower than the Rcpp code Giovanni
wrote, but you still should be using the return value of the C++
function to get the result back to R, rather than writing directly into
function arguments.
Here is a simple benchmark, comparing Giovanni's 'f_cpp' function with
my solution using 'sweep', with N=10000, M=500, K=10:
> XX <- array(x, dim=c(N, M, K))
> benchmark(f_cpp=f_cpp(x, mu), sweep=sweep(XX, 2:3, mu),
columns=c('test', 'relative', 'elapsed'),
order='relative')
test relative elapsed
1 f_cpp 1.000 69.306
2 sweep 15.269 1058.241
I would have also used f_R, but it runs out of memory with these dimensions.
Davor
On 12-10-18 08:02 PM, Davor Cubranic wrote:
Giovanni,
Anything you pass from R into your C++ routine should be considered
read-only. It's typically a bad idea to do in-place updates of the R
object you passed in as an argument.
In addition, I doubt that the code you wrote (addition/subtraction of
matrix columns) is going to be that much faster in C++ over R simply
because R already does those operations with optimized C-level code.
Maybe not as you wrote it, but you can use available R's functionality a
lot more effectively with something like this:
sweep(array(x, dim=c(N, M, K)), 2:3, mu)
This should do exactly what you're trying to do with your C++ function,
except storing it as a 3D array instead of a list of 2D matrices. (It's
more efficient and works better with everything else in R.)
If it's the large number of iterations in your Gibbs sampler that you're
worried about, then write the entire algorithm in C++ and return the
final result.
Davor
On 12-10-18 07:28 PM, Giovanni Petris wrote:
Thanks Dirk!
By the way, I just figured out a working solution for what I wanted to
achieve - 'code_f' below.
Best,
Giovanni
=====
### Small reproducible example
K <- 2; N <- 5; M <- 3
x <- matrix(rnorm(N * M), N, M)
mu <- matrix(rnorm(M * K, mean = 1 : K), M, K, TRUE)
### What I would do in R
f_R <- function() {
## center observations in 'x' using columns of 'mu'
w <- lapply(1L:K, function(i) x - rep(mu[, i], each = N))
## eventually more stuff here...
w
}
### For the Rcpp part, set aside working space in the global environment
workspace <- lapply(1:K, function(k) matrix(0.0, N, M))
code_f <- '
Rcpp::NumericMatrix cpp_x(x);
Rcpp::NumericMatrix cpp_mu(mu);
int K = cpp_mu.ncol();
int M = cpp_mu.nrow();
int N = cpp_x.nrow();
Environment glob = Environment::global_env();
List cpp_work(glob.get("workspace"));
NumericMatrix tmpMatrix;
for (int k=0; k < K; k++) {
tmpMatrix = as<SEXP>(cpp_work[k]);
for (int i=0; i < M; i++)
tmpMatrix(_, i) = cpp_x(_, i) - cpp_mu(i, k);
}
return R_NilValue;
'
## library(Rcpp); library(inline)
f_cpp <- cxxfunction(signature(x = "numeric", mu = "numeric"), code_f,
plugin = "Rcpp")
## try it
f_cpp(x, mu)
all.equal(workspace, f_R())
________________________________________
From: Dirk Eddelbuettel [e...@debian.org]
Sent: Thursday, October 18, 2012 9:15 PM
To: Giovanni Petris
Cc: Dirk Eddelbuettel; rcpp-devel@lists.r-forge.r-project.org
Subject: RE: [Rcpp-devel] How to modifying the elements of a list in
the global environment?
On 19 October 2012 at 01:32, Giovanni Petris wrote:
| Hi Dirk,
|
| Thank you for the quick reply. I will look for more examples on the
net.
|
| About your suggestion of allocating scrap space inside the C++
routine, am I wrong to think that when the matrices are large and the
function is called repeatedly within a Gibbs sampler loop, this is not
a very efficient approach?
Yes, if you want to use it several times you can pass it around for
reuse.
As I mentioned, a reference is just a pointer and hence unrelated to
the size
of the object.
Dirk
| Thanks again!
|
| Best,
| Giovanni
|
| ________________________________________
| From: Dirk Eddelbuettel [e...@debian.org]
| Sent: Thursday, October 18, 2012 8:06 PM
| To: Dirk Eddelbuettel
| Cc: Giovanni Petris; rcpp-devel@lists.r-forge.r-project.org
| Subject: Re: [Rcpp-devel] How to modifying the elements of a list in
the global environment?
|
| Giovanni,
|
| On 18 October 2012 at 17:46, Dirk Eddelbuettel wrote:
| | | > code_f <- '
| | | + Rcpp::NumericMatrix cpp_x(x);
| | | + Rcpp::NumericMatrix cpp_mu(mu);
| | | + int K = cpp_mu.ncol();
| | | + int M = cpp_mu.nrow();
| | | + int N = cpp_x.nrow();
| | | + Environment glob = Environment::global_env();
| | | + List cpp_work(glob.get("workspace"));
| | | +
| | | + for (int k=0; k < K; k++) {
| | | + cpp_work[k] = clone(cpp_x);
| | | + for (int i=0; i < M; i++)
| | | + cpp_work(_, i) = cpp_work(_, i) - cpp_mu(i, k);
|
| That can't work. You have the types confused.
|
| You could just pass the list object down, and then pick from it.
Working
| examples for that are in the list archives, in the RcppExamples
package, and
| in other places on the net.
|
| Or, as I suggested, just allocate the scrap space inside the C++
routine.
|
| Hth, Dirk
|
| --
| Dirk Eddelbuettel | e...@debian.org | http://dirk.eddelbuettel.com
--
Dirk Eddelbuettel | e...@debian.org | http://dirk.eddelbuettel.com
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
