Is it possible to vectorise the sequential update

问题:

Is it possible to vectorise code like the following?

length(x) <- 100;
x[1]      <- 1;
y         <- rnorm(100);

for(i in 2:100) {
    x[i] <- 2 * y[i] * x[i-1];
}

I appreciate that this is a trivial example, but it serves to illustrate the idea.

I often need to write code where the i-th value in a vector depends on the (i-1)-th value and if possible, I'd like to write this without needing a for loop, as profiling suggests the functions with this type of operation are the major bottlenecks in my code.

Is this operation vectorizable so I do not need to use a for() loop in the calculation?

回答1:

In general, if you want a vectorised solution you need to solve the recurrence relation.

回答2:

In the example you have you could work out the formula for x[i] and see if it can be vectorized. In this case I think cumprod might work.

x <- c(1, cumprod(2*y)[1:99])

For some cases case you can also use the filter command in convolution or recursive mode. See ?filter

However if it is isn't possible to work out a formula for the n'th value that fits one of the molds above, you could try using a package like inline or Rcpp to write this in loop in C/C++.

回答3:

The interior of this plot command is equivalent. Rather interesting to repeatedly run it:

plot(c(1, 2^(2:length(x)-1)*cumprod(rnorm(99) )) )

回答4:

I don't have full details on this yet, but it looks the function filter() is going to be useful to do what I need.

回答5:

You can write the non-vertorized code in C++:

library(inline)
myfun <- cxxfunction(signature(y="numeric"), body='
Rcpp::NumericVector yvec(y);
int ysize = yvec.size();
Rcpp::NumericVector result(ysize);
if (ysize > 0) {
    result[0] = 1;
    for (int i = 1; i < ysize; i++) {
        result[i] = 2 * yvec[i] * result[i-1];
    }
}
return result;
', plugin="Rcpp")

Then call this function from R:

y <- rnorm(100);
x <- myfun(y);