Understand and avoid infinite recursion R

Question

I can't 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 to run a double 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 was running fine but now I hoped to find the k that would best fit my criterion. So I decided to turn it into a function so that i could later optimize or uniroot it. I wrote something like that :

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 infinite recursion errors.

I read about max memory but I am sure there is a more aesthetic way of doing it.

My reproducible example :

library(memoise)

gradient <- function(x,y,tau){if (x-y > 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 (of length 2080)

A is a 1414 x 24 matrix

I hope that helps.

Vincent Zoonekynd · Accepted Answer · 2013-05-30T19:16:08.420

There are three problems.

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

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

The following would 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 those values will be re-computed an exponential number of times.

f(500) # Too long

In more abstract terms, consider the graph whose vertices are 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) a very large number of times. The following code draw 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)

Function call graph

One workaround is to store the values you have already computed to avoid recomputing 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 the function, the number of recursive calls becomes linear, but it can still be too large. You can increase the limit as suggested by the initial error message:

options( expressions = 5e5 )

If this is not sufficient, you can prepopulate the table, by using increasingly larger values of i. With 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

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

library(memoise)
gradient <- function(x,y,tau){
  if (x-y > 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)

Hi and thank you, I applied what you said but I still get : Error: evaluation nested too deeply: infinite recursion / options(expressions=)? — user1627466, May 30 '13 at 10:19
Please provide a [reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). — Vincent Zoonekynd, May 30 '13 at 10:40
@user1627466 try calculating the first N terms (times thru loop) and plotting the results to see whether there's any sort of converging trend. The part V.Zoonekynd said about picking an appropriate initial condition is very important as well. Try plotting the data from your loop-based setup to see what range of `k` is likely to enclose the optimal value. — Carl Witthoft, May 30 '13 at 11:36
@Carl Thank I'll try that. Vincent : It works. I don't really understand the reasoning behind it but it works. I'm surprised R does not memorize past results on its own. — user1627466, May 30 '13 at 11:44

Understand and avoid infinite recursion R

1 Answers1

Linked