R Error: function returns NA at 3.07653253930756e+181 distance from x

I am using the R programming language.

Based on the code following website https://www.stat.cmu.edu/~ryantibs/statcomp-F15/lectures/optimization.pdf (note: for some reason, this website does not open in Google Chrome - please try Microsoft Edge Explorer), I am trying to use the "Gradient Descent Optimization Algorithm" to optimize (i.e. find the minimum value) of the function : f(x) = x^3 - 2x - 5

I first defined the function that I want to optimize:

#define function to be optimized
func2 <- function(x) {
 return( x[1]^3 - 2* x[1] - 5)
}

Next, I defined the function for the Gradient Descent Optimization Algorithm:

#load library
library(numDeriv)


#define gradient descent function
grad.descent = function(f, x0, max.iter=200, step.size=0.05,
stopping.deriv=0.01, ...) {
n = length(x0)
xmat = matrix(0,nrow=n,ncol=max.iter)
xmat[,1] = x0
for (k in 2:max.iter) {
# Calculate the gradient
grad.cur = grad(f,xmat[,k-1],...)
# Should we stop?
if (all(abs(grad.cur) < stopping.deriv)) {
k = k-1; break
}
# Move in the opposite direction of the grad
xmat[,k] = xmat[,k-1] - step.size * grad.cur
}
xmat = xmat[,1:k] # Trim
return(list(x=xmat[,k], xmat=xmat, k=k))
}

Finally, I tried to optimize the function :

# I think this serves as an initialization value
 x0 = c(-1.9)

#run gradient descent algorithm
 gd = grad.descent(func2,x0,step.size=1/3)

Problem : But this returns the following error:

 Error in grad.default(f, xmat[, k - 1], ...) : 
  function returns NA at 3.07653253930756e+181 distance from x. 

Can someone please show me what I am doing wrong?

Thanks!

If you want to impose boundaries, you can do it like this:

#load library
library(numDeriv)

#define gradient descent function
grad.descent = function(f, x0, max.iter=200, step.size=0.05,
                        stopping.deriv=0.01, boundaries = NULL, verbose = TRUE, ...) {
  n = length(x0)
  xmat = matrix(0,nrow=n,ncol=max.iter)
  xmat[,1] = x0
  for (k in 2:max.iter) {
    if (verbose) message(paste(xmat[, k-1], collapse = ", "))
    # Calculate the gradient
    grad.cur = grad(f,xmat[,k-1],...)
    # Should we stop?
    if (all(abs(grad.cur) < stopping.deriv)) {
      k = k-1; break
    }
    # Move in the opposite direction of the grad
    xmat[,k] = xmat[,k-1] - step.size * grad.cur
    if (!is.null(boundaries)) {
      xmat[,k] <- ifelse(xmat[,k] < boundaries[1], boundaries[1], xmat[,k])
      xmat[,k] <- ifelse(xmat[,k] > boundaries[2], boundaries[2], xmat[,k])
      if (all(xmat[, k] == xmat[, k-1] | abs(grad.cur) < stopping.deriv))) break #stop if boundaries
    }
  }
  xmat = xmat[,1:k, drop = FALSE] # Trim
  return(list(x=xmat[,k], xmat=xmat, k=k))
}

# starting values
x0 = c(-1.9, 1.9)

#use functions that are actually vectorized 
#if you want to use multiple starting values
f1 <- \(x) x^3 - 2* x - 5
grad.descent(f1,x0,step.size=1/3, boundaries = c(-5, 5))

f2 <- \(x) x^2
grad.descent(f2,x0,step.size=1/3, boundaries = c(-5, 5))