Understand and avoid infinite recursion R

Question

I cannot find any advice on how to deal with infinite recursion in R. I would like to illustrate my problem in the most general way so that others can benefit from it. Feel free to edit it.

I used space for loop

for (k in 1:n){ for (i in 1:m){ f[i,,k] <- an expression that depends on g[i-1,,k] g[i,,k] <- a mixture of g[i-1,,k] and f[i,,k]}} 

This worked fine, but now I was hoping to find the k that best met my criteria. Therefore, I decided to turn it into a function, in order to subsequently optimize or unify it. I wrote something like this:

 f <- function(i,k){an expression that depends on g(i-1,k)} g <- function(i,k){an expression that depends on g(i-1,k) and f(i,k)} 

I thought the two problems were similar, but to my great surprise, I get endless recursion errors.

I read about maximum memory, but I'm sure there is a more aesthetic way to do this.

My reproducible example:

 library(memoise) gradient <- function(x,y,tau){if (xy > 0) {- tau} else {(1-tau)}} aj <- c(-3,-4,-2,-3,-5,-6,-4,-5,-1,rep(-1,15)) f <- function(x,vec){sum(x^vec)-1} root <- uniroot(f, interval=c(0,200), vec=aj)$root memloss<-function(i,k){if (i==1) {c(rep(0,24))} else if (i <= 0 | k < -5) {0} else {gradient(dailyreturn[i-1],weight(i-1,k)%*%A[i-1,],0.0025)*A[i-1,]}} memweight <- function(i,k){if (i==1) {c(rep(root,24)^aj)} else if (i <= 0 | k < -5) {0} else {(exp(- (2^(k)/(sqrt(1415))) * loss(i,k))) / (weight(i-1,k) %*% exp(- 2^(k)/(sqrt(1415)) * loss(i,k)) ) * weight(i-1,k)}} loss <- memoize(memloss) weight <- memoize(memweight) 

where dailyreturn is a vector (2080 long)

A - matrix 1414 x 24

I hope this helps.

Answer 1

There are three problems.

First you need an initial case for your recursion. The following leads to infinite recursion (the value of i constantly decreasing, but never stops).

 f <- function(i) g(i-1) g <- function(i) g(i-1) + f(i) f(5) 

The following will stop.

 f <- function(i) g(i-1) g <- function(i) if( i <= 0 ) 0 else g(i-1) + f(i) f(5) 

The second problem is that some of these values ​​will be recounted exponentially.

 f(500) # Too long 

In more abstract terms, consider a graph whose vertices f(i) and g(i) , for all values ​​of i , with edges corresponding to function calls. Recursion allows you to explore this graph as if it were a tree. But in this case, it is not a tree, and you end up evaluating the same function (exploring the same node) many times. The following code draws this graph.

 library(igraph) n <- 5 g <- graph.empty() g <- g + vertices( paste0("f(", 1:n, ")" ) ) g <- g + vertices( paste0("g(", 0:n, ")" ) ) for( i in 1:n) { g <- g + edge( paste0("f(", i ,")"), paste0( "g(", i-1, ")" ) ) g <- g + edge( paste0("g(", i ,")"), paste0( "f(", i, ")" ) ) g <- g + edge( paste0("g(", i ,")"), paste0( "g(", i-1, ")" ) ) } plot(g) 

One workaround is to save the values ​​that you have already calculated to avoid recalculating them: this is called memoization .

 library(memoise) f <- function(i) G(i-1) g <- function(i) if( i <= 0 ) 1 else G(i-1) + F(i) F <- memoize(f) G <- memoize(g) f(500) 

When you memoise a function, the number of recursive calls becomes linear, but it can still be too large. You can increase the limit as suggested in the original error message:

 options( expressions = 5e5 ) 

If this is not enough, you can pre-populate the table using ever larger i values. In your example:

 options( expressions = 5e5 ) loss(1000,10) # Does not work: Error: protect(): protection stack overflow loss(500,10) # Automatically stores the values of loss(i,100) for i=1:500 loss(1000,10) # Works 

Thirdly, there may be function calls that unnecessarily increase the size of the call stack. In your example, if you type traceback() after an error, you will see that many intermediate functions are in the call stack because weight(i,k) and loss(i,k) are used inside the function arguments. If you move these calls outside the function arguments, the call stack is smaller and it seems to work.

 library(memoise) gradient <- function(x,y,tau){ if (xy > 0) { - tau } else { (1-tau) } } aj <- c(-3,-4,-2,-3,-5,-6,-4,-5,-1,rep(-1,15)) f <- function(x,vec){sum(x^vec)-1} root <- uniroot(f, interval=c(0,200), vec=aj)$root memloss<-function(i,k){ cat( "loss(", i, ",", k, ")\n", sep="" ) if (i==1) { c(rep(0,24)) } else if (i <= 0 | k < -5) { 0 } else { w <- weight(i-1,k) # Changed gradient(dailyreturn[i-1],w%*%A[i-1,],0.0025)*A[i-1,] } } memweight <- function(i,k){ cat( "weight(", i, ",", k, ")\n", sep="" ) if (i==1) { c(rep(root,24)^aj) } else if (i <= 0 | k < -5) { 0 } else { w <- weight(i-1,k) # Changed l <- loss(i,k) # Changed (exp(- (2^(k)/(sqrt(1415))) * l)) / (w %*% exp(- 2^(k)/(sqrt(1415)) * l) ) * w } } loss <- memoize(memloss) weight <- memoize(memweight) A <- matrix(1, 1414, 24) dailyreturn <- rep(1,2080) options( expressions = 1e5 ) loss(1400,10)