Mar 12

## How to Vectorize Nested Loop in R?

Posted by abiao at 14:24 | Code » R/Splus | Comments(9) | Reads(19994)
Could any R expert here help me to vectorize my for loop? Thanks in advance for your favor. The reason I am in trouble is the variable inside my "for" function are updated after each loop, which makes me feel difficult to use lapply, sapply or whatever.

Simplifed codes are listed below:
for (i in 1:N) { #N could be a large number, AdjS and VarS are initially given and updated for each i
PredY <- H %*% PredS
PredError <- (Y[i,] - t(PredY))
VarY <- (H %*% VarS) %*% t(H)
InvVarY <- solve(VarY)
KG <- (VarS %*% t(H)) %*% InvVarY
VarS <- (diag(3) - KG %*% H) %*% VarS
ll[i] <- PredError %*% InvVarY %*% t(PredError)
}

The point is how to vectorize the for loop while allowing AdjS and VarS to be updated. I appreciate your help.

Tags: ,
You ought to be using a cholesky decomposition and backsolve/forwardsolve for KG and II, rather than inverting a matrix.Besides this, if matrices are not sparse, you could try RcppArmadillo to apply a matrix templating library that will try to optimize some of this, and limit memory allocation. It will also speed up the looping.
btw, if this , perchance, is a kalman filter implementation, (the names of the variables are quite suggestive... KG, Pred*, etc.) there is an extensive lit on stable computations, incl cholesky and svd based, as well as several fast C implementations for R. why reinvent the wheel for an application (it looks like you've hard-coded in a 3 by 3 matrix, so you aren't looking for a general solution).
thanks a lot, BC for your comments, it is indeed for kalman filter,  I found a very good package FKF which is written in C and relies fully on linear algebra subroutines contained in BLAS and LAPACK, so quite efficient for my purpose. Back to the question, is it possible to use laaply or other "apply" function for this loop case? let's forget kalman filter. thank you.
you could write a function and work with global variables, i.e. your current code looks like:

tmp<-100
start<-100
for (j in 1:4) {
i[j]<-start*tmp
tmp<-tmp/2
}

which is the same as:

tmp<-100
test<-function(x) {
tmp2<-x*tmp
tmp<<-tmp/2
return(tmp2)
}
sapply(1:4,function(x) i[x]<-test(100))

just to add some explantory stuff:
the difference between <- and <<- is that the first one is locally defined (i.e. within a function) and the second one is globally defined (i.e. within and outside of a function) and thus you can change the global variables and can work with some apply functions.
thank you so much, Chris, I am testing it, have a nice weekend.
I'd like to see timings on this compared to the original though -- as one of the major reasons for vectorization is to avoid lookups/indirection, and facilitate pre-allocation of memory. Once you are repeatedly storing and overwriting something in a global variable, you are adding back lookups (to find the global variable) and memory allocation (for each array modification of II, which in R copies the entire array). Plus, the original code was much easier to understand.You can find an example on the Rcpp devel list, that I put in for my "Site URI", that uses the inline package and matrix templating, in this discussion:If you explore this, be sure to read the rest of the thread where Romain Francois, the creator of Rcpp sugar, suggests how to improve the code.
Sure, BC, there is always a trade-off between time spent on programming and on code running, I can't agree more on your perception of the importance of easy-to-understood codes. In the next few days I will test and compare the performance of both methods suggested by you and Chris.A single post teaches me at least 2 things: <<- for global variable and Rcpp library, what a fruitful weekend. Thanks guys.
patmmccann
You should ask this on stackexchange, those guys are experts.
instead of  "<<-" I sugest to create a new environment and than use assign. It is explicit and safer.
Pages: 1/1 1
 Add a comment Emots Enable HTML Enable UBB Enable Emots Hidden Remember Nickname   Password   Optional Site URI   Email   [Register]