Rosenbrock Function

rosenbrock.r

rosen_height <- function(x,y) {
  100 * (y-x^2)^2 + (1-x)^2
}

gradient <- function(x,y) {
  val_x <- -400*x*(y-x^2)-2*(1-x)
  val_y <- 200*(y-x^2)
  return( c(val_x, val_y) )
}

# Check function looking for insufficient decrease in height
insufficient_decrease <- function(x,y,step_len, dir_x, dir_y, const_1){
  new_height <- rosen_height( x + step_len * dir_x, y + step_len * dir_y )

  # something proportional to directional derivative.
  #   result will be 1x1 matrix, and we need a number.
  directional_prop <- (gradient(x,y) %*% c(dir_x, dir_y))[1,1]

  hurdle <- rosen_height(x, y) + const_1 * step_len * directional_prop
  return( new_height > hurdle )
}

# Check function looking for when there is no curvature
no_curvature <- function(x,y,step_len, dir_x, dir_y, const_2){
  # new directional derivative
  lhs <- gradient(x + step_len*dir_x, y + step_len*dir_y) %*% c(dir_x,dir_y)
  # old directional derivative * stringency constant
  rhs <- const_2 * gradient(x,y) %*% c(dir_x,dir_y)

  # inner products cast to matrix, taking [1,1] gets us back to simple boolean
  result <- (lhs < rhs)[1,1]
  return(result)
}

# Main function
rosenbrock <- function(init_x, init_y, step_len, dir_x, dir_y, const_1, const_2){
  result <- NA
  while(is.na(result)){
    if( insufficient_decrease(init_x, init_y, step_len, dir_x, dir_y, const_1) ){
      # Update step lengths
      step_max <- step_len
      step_len <- (step_min + step_max) / 2
    } else if( no_curvature(init_x, init_y, step_len, dir_x, dir_y, const_2) ){
      step_min <- step_len
      if(step_max == Inf){
        step_len <- 2 * step_len
      } else {
        step_len <- (step_min + step_max) / 2
      }
    } else {
      result <- step_len
    }
  }
  return(result)
}

#################
# Smoke Testing #
#################
# Starting point
init_x <- 1 #formerly x0[1]
init_y <- 3 #formerly x0[2]

step_len <- 1 # formerly a0
step_min <- 0 # formerly amin
step_max <- Inf # formerly amax

# Constants controlling the stringency of the optimization
const_1 <- 0.1 # formerly c1
const_2 <- 0.5 # formerly c2

# direction, p
dir_x <- 2 # formerly p[1]
dir_y <- 1 # formerly p[2]

# Do the thing! Expect 0.5 result
rosenbrock(init_x, init_y, step_len, dir_x, dir_y, const_1, const_2)