In an earlier post, I introduced the golden section search method – a modification of the bisection method for numerical optimization that saves computation time by using the golden ratio to set its test points.  This post contains the R function that implements this method, the R functions that contain the 3 functions that were minimized by this method, and the R script that ran the minimization.

I learned some new R functions while learning this new algorithm.

– the curve() function for plotting curves

– the cat() function for concatenating strings and variables and, hence, for printing debugging statements

As usual, I made sure to change the working directory in R to the folder that contains all of these files, and I used the source() function to call the relevant functions in my script for running everything.

Here is the function that I minimized in my earlier blog post; I titled it “f.R”.

f = function(x)
{
   abs(x - 2) + (x - 1)^2
}

Here is the function that implemented the golden section search method; I called it “golden.section.search.R”.

##### Implementing the golden section search method
##### a modification of the bisection method with the golden ratio
##### By Eric Cai - The Chemical Statistician
golden.section.search = function(f, lower.bound, upper.bound, tolerance)
{
   golden.ratio = 2/(sqrt(5) + 1)

   ### Use the golden ratio to set the initial test points
   x1 = upper.bound - golden.ratio*(upper.bound - lower.bound)
   x2 = lower.bound + golden.ratio*(upper.bound - lower.bound)

   ### Evaluate the function at the test points
   f1 = f(x1)
   f2 = f(x2)

   iteration = 0

   while (abs(upper.bound - lower.bound) > tolerance)
   {
      iteration = iteration + 1
      cat('', '\n')
      cat('Iteration #', iteration, '\n')
      cat('f1 =', f1, '\n')
      cat('f2 =', f2, '\n')

      if (f2 > f1)
      # then the minimum is to the left of x2
      # let x2 be the new upper bound
      # let x1 be the new upper test point
      {
         cat('f2 > f1', '\n')
         ### Set the new upper bound
         upper.bound = x2
         cat('New Upper Bound =', upper.bound, '\n')
         cat('New Lower Bound =', lower.bound, '\n')
         ### Set the new upper test point
         ### Use the special result of the golden ratio
         x2 = x1
         cat('New Upper Test Point = ', x2, '\n')
         f2 = f1

         ### Set the new lower test point
         x1 = upper.bound - golden.ratio*(upper.bound - lower.bound)
         cat('New Lower Test Point = ', x1, '\n')
         f1 = f(x1)
      } 
      else 
      {
         cat('f2 < f1', '\n')
         # the minimum is to the right of x1
         # let x1 be the new lower bound
         # let x2 be the new lower test point

         ### Set the new lower bound
         lower.bound = x1
         cat('New Upper Bound =', upper.bound, '\n')
         cat('New Lower Bound =', lower.bound, '\n')

         ### Set the new lower test point
         x1 = x2
         cat('New Lower Test Point = ', x1, '\n')

         f1 = f2

         ### Set the new upper test point
         x2 = lower.bound + golden.ratio*(upper.bound - lower.bound)
         cat('New Upper Test Point = ', x2, '\n')
         f2 = f(x2)
      }
   }

   ### Use the mid-point of the final interval as the estimate of the optimzer
   cat('', '\n')
   cat('Final Lower Bound =', lower.bound, '\n')
   cat('Final Upper Bound =', upper.bound, '\n')
   estimated.minimizer = (lower.bound + upper.bound)/2
   cat('Estimated Minimizer =', estimated.minimizer, '\n')
}

Here is the script that ran everything; I called it “minimization.R”.

##### Finding the minimizers of functions using the bisection method with the golden ratio
##### By Eric Cai - The Chemical Statistician

# Calling the user-defined functions in the working directory
source('golden.section.search.R')
source('f.R')

# printing the PNG images into the working directory
png('INSERT YOUR DIRECTORY PATH HERE/cusped function.png')
# plotting the curve of my user-defined function
curve(f, from = 1, to = 3, main = expression(paste('f(x) = |x - 2| + (x - 1)'^'2')))
dev.off()

# finding the minimizer of my user-defined function using my golden bisection method
golden.section.search(f, 1, 3, 1e-5)