Quantitative finance collector
C++ Matlab VBA/Excel Java Mathematica R/Splus Net Code Site Other
Mar 12

How to Vectorize Nested Loop in R?

Posted by abiao at 14:24 | Code » R/Splus | Comments(9) | Reads(21397)
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
        PredS <- F %*% AdjS                                                                              
        PredY <- H %*% PredS
        PredError <- (Y[i,] - t(PredY))
        VarY <- (H %*% VarS) %*% t(H)
        InvVarY <- solve(VarY)
        KG <- (VarS %*% t(H)) %*% InvVarY
        AdjS <- PredS + PredError
        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 First page 1 Final page
Add a comment
Emots
Enable HTML
Enable UBB
Enable Emots
Hidden
Remember
Nickname   Password   Optional
Site URI   Email   [Register]