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.

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 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)

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)