To write legible answers you will need to be familiar with both Markdown and Latex

Before you turn this problem in, make sure everything runs as expected. To do so, restart the kernel and run all cells (in the menubar, select Runtime→→Restart and run all).

Show your work!

Whenever you are asked to find the solution to a problem, be sure to also show how you arrived at your answer.

Resources

[1] Different color plots plotly

[2] Plotting surfaces with plotly

[3] Plotting scatter plots with plotly

[4] Color scales

%matplotlib inline

NAME = ""
STUDENT_ID = ""

Local Search on the Ackley Surface

---

Background

In many optimization problems, the path to the goal is irrelevant; the goal state itself is the solution. Local search is widely used for problems with very large search spaces, the solutions returned are good but usually not optimal or the best. Local search algorithms keep a single "current" state or a small set of states and iteratively try to improve on them. These makes this class of algorithms usually very memory efficient.

In this problem we explore the Ackley function, a non-convex function typically used as a test for optimization algorithms.

The Ackley function has many local optima but only one global optimum that is $f(0,0) = 0$. Interestingly enough it might not be that easy to find this solution.

The Ackley function is defined as:

$$f(x,y):= -20exp[-0.2\sqrt{0.5(x^2 + y^2)}] -exp[0.5(cos(2\pi x) + cos(2\pi y))] + e + 20$$

The figure below illustrates the result of running three optimization methods on the Ackley Function. These solutions were found using methods unknown methods A, B and C.

Learning Objectives

• Ability to implement Stochastic Hill Climbing with Restarts (SHCR), Simulated Annealing (SA), and Local Beam Search (LBS)
• Understand the tradeoffs and differences between the three algorithms
• Develop the ability to make principled judgements in optimization problems

Deliverables

In this question, you will implement three optimization algorithms:

• Stochastic Hill Climbing with Restarts (SHCR)
• Simulated Annealing (SA)
• Local Beam Search (LBS)

1. Complete this Notebook with implementations of SHCR, SA and LBS.
2. For each algorithm visualize the candidate solutions found at each

step using the provided create_ackley_figure() function. Please see the image above for an example of an image that should be included in your submission.

import numpy as np
from numpy import arange
from numpy import exp
from numpy import sqrt
from numpy import cos
from numpy import e
from numpy import pi
from numpy import meshgrid
from numpy.random import randn
from numpy.random import rand
from numpy.random import seed
from numpy import asarray
import plotly.graph_objects as go

# ackley objective function
def f(x, y):
	return -20.0 * exp(-0.2 * sqrt(0.5 * (x**2 + y**2))) - exp(0.5 * (cos(2 * pi * x) + cos(2 * pi * y))) + e + 20

def create_ackley_figure(solution_points=np.array([]), soln_names=['A','B,','C']):
  """
  :param solution_points A numpy array of dimensions [i, n, 2]
  where i is the number of solution points you want to plot and n is the 
  number of steps used in estimating your solution. In each case use the same 
  number of steps if you want to plot all solutions on the same figure.

  :return An interactive Ackley function figure with/without solution points.
  """
  # define range for input
  r_min, r_max = -5.0, 5.0
  # sample input range uniformly at 0.1 increments
  xaxis = arange(r_min, r_max, 0.1) # (100,)
  yaxis = arange(r_min, r_max, 0.1) # (100,)
  # create a mesh from the axis
  x, y = meshgrid(xaxis, yaxis)

  # compute targets
  results = f(x, y)

  figure = go.Figure(data=[go.Surface(z=results, x=x, y=y,colorscale='deep',showscale=False)])
  figure.update_layout(title='Ackley Function', autosize=False,
                  width=500, height=500,
                  title_x=0.5,
                  margin=dict(l=65, r=50, b=65, t=90))
  
  if solution_points.size:
    soln_colors = ['Reds', 'Rainbow', 'Viridis']
    for i in range(len(solution_points)):
      x_sln = solution_points[i][:,0]
      y_sln = solution_points[i][:,1]
      zdata = f(x_sln, y_sln)
      figure.add_scatter3d(name=soln_names[i], x=x_sln, y=y_sln, z=zdata, mode='markers',
      marker=dict(size=1, color=x_sln.flatten(), colorscale=soln_colors[i]))
      
      figure.update_layout(showlegend=True)
      
  return figure

### RUN THIS ###
"""
This is the Ackley Function.
"""
ack_figure = create_ackley_figure()
ack_figure

General Function Definitions

# objective function
def objective(v):
  """
  This function defines the Ackley surface, it has a minimum of 0. This 
  can be used to test if we have found points that yield the minimum value.
  :param v, a tuple representing a 2D point being tested.

  returns a value in the range [-5.0, 5.0]
  """
  x, y = v
  return -20.0 * exp(-0.2 * sqrt(0.5 * (x**2 + y**2))) - exp(0.5 * (cos(2 * pi * x) + cos(2 * pi * y))) + e + 20

  
# check if a point is within the bounds of the search
def in_bounds(point, bounds=asarray([[-5.0, 5.0], [-5.0, 5.0]])):
  """
  It is possible our optimzation method could move us outside of our search 
  space, so it is helpful to check.
  :param point a tuple representing a 2D point
  :param a 2D array, that describes the bounds of each dimension of our point.
  returns a Boolean of whether the point lies in the bounds or not.
  """
	# enumerate all dimensions of the point
  for d in range(len(bounds)):
    # check if out of bounds for this dimension
    if point[d] < bounds[d, 0] or point[d] > bounds[d, 1]:
      return False
  return True

Part 1) Stochastic Hill Climbing with Restarts

1. Implement SHCR in the cell below.
2. Visualize and comment on the candidate points visited.

# stochastic hill climbing with restarts algorithm
def stochastic_hillclimbing(objective, bounds, n_iterations, step_size, start_pt):
  """
  :param objective, the function we are trying to optimize
  :param bounds, the boundaries of the problem
  :param n_iterations number of times to repeat stochastic hill climbing
  :param step_size how much we should move
  :param start_pt the point we start optimizing from.
  returns [solution, solution_value, candidates]
  """
  ### YOUR CODE HERE ###
  return

def random_restarts(objective, bounds, n_iter, step_size, n_restarts):
  """
  :param objective, the function we are trying to optimize
  :param bounds, the boundaries of the problem
  :param n_iter number of times to repeat stochastic hill climbing 
  :param step_size how much we should move
  :param n_restarts, the number of times we restart.
  returns [best solution, best solution value, visited points]
  """
  ### YOUR CODE HERE ###
  return -1, -1, -1

# seed the pseudorandom number generator
seed(240)
# define range for input
bounds = asarray([[-5.0, 5.0], [-5.0, 5.0]])
# define the total iterations
n_iter = 1000
# define the maximum step size
step_size = 0.05
# total number of random restarts
n_restarts = 30
# perform the hill climbing search
best, score, shcr_candidates = random_restarts(objective, bounds, n_iter, step_size, n_restarts)
print('Done!')
print('f(%s) = %f' % (best, score))

# An example of plotting solutions from 3 optimziation methods.
# all_method_candidates has shape [3, 1000, 2], the first dim is the number of solutions
# the second dim is the number of steps taken in each method.
shcr_candidates = np.random.rand(100,2) # EXAMPLE PURPOSES ONLY replace with your actual solution.
shcr_candidates = shcr_candidates.reshape(-1,2)
all_method_candidates = np.array([shcr_candidates])
ack_figure = create_ackley_figure(all_method_candidates)
ack_figure

[YOUR ANSWER HERE]

Part 2) Simulated Annealing

1. Implement SA in the cell below.
2. Visualize and comment on the candidate points visited.

def get_temperature_schedule(decayRate, temp, diff):
  """
  :param decayRate A float in range [0,1) used to decrease the temperature
  :param temp Temperature
  :param diff Difference between value of candidate solution and value of best solution so far.
  """
  # calculate temperature for current epoch
  t = temp*decayRate
  # calculate acceptance criterion
  criterion = exp(-diff / t)
  return criterion, t

# simulated annealing algorithm
def simulated_annealing(objective, bounds, n_iterations, step_size, temp):
  """ 
  :param objective, the function we are trying to optimize
  :param bounds, the boundaries of the problem
  :param n_iterations number of times to repeat stochastic hill climbing
  :param step_size how much we should move
  :param temp temperature for each epoch
  returns [solution, solution_value, candidates] 
  """
  ### YOUR CODE HERE ###
  return -1, -1, -1

# seed the pseudorandom number generator
seed(240)
# define range for input
bounds = asarray([[-5.0, 5.0], [-5.0, 5.0]])
# define the total iterations
n_iterations = 1000
# define the maximum step size
step_size = 0.1
# initial temperature
temp = 10
# perform the simulated annealing search
best, score, sa_candidates = simulated_annealing(objective, bounds, n_iterations, step_size, temp)
print('Done!')
print('f(%s) = %f' % (best, score))

sa_candidates = np.random.rand(40,2)
sa_candidates = sa_candidates.reshape(-1,2)
all_method_candidates = np.array([shcr_candidates,sa_candidates])
ack_figure = create_ackley_figure(all_method_candidates)
ack_figure

[YOUR ANSWER HERE]

Part 3) Local Beam Search

1. Implement LBS in the cell below.
2. Visualize and comment on the candidate points visited.
3. Compare all 3 implementations by commenting on the distribution of points on the Ackley surface and the empirical run time of each method.

# generate the neighbors based on step size
def generate_neighbors(point, step_size=0.1):
    point = point.reshape(-1,1)
    possible_steps = possible_steps = np.array([[step_size*i for i in range(-10,11)],
                           [step_size*i for i in range(-10,11)]])
    neighbors = point + possible_steps
    return neighbors

def local_beam_search(objective, bounds, step_size, k, n_iterations):
  """ 
  :param objective, the function we are trying to optimize
  :param bounds, the boundaries of the problem
  :param k, how many candidates to consider
  :param n_iterations, how long to search for.
  returns [solution, solution_value, candidates] 
  """
  ### YOUR CODE HERE ###
  return -1, -1, -1

# seed the pseudorandom number generator
seed(240)
# define range for input
bounds = asarray([[-5.0, 5.0], [-5.0, 5.0]])
# define the total iterations
n_iterations = 10
# define the maximum step size
step_size = 0.1
# candidates to consider
k = 1
# perform the local beam search
sequences = local_beam_search(objective, bounds, step_size, k, n_iterations)

[YOUR ANSWER HERE]

Part 4) Create a function for which SHCR tends to outperform SA and LBS.

• The function should accept two variables def f(x, y): <your code here>.
• Evaluate the performance of each algorithm by 10 trials for each search algorithm and report the distributions of solutions found.

def schr_best(x, y):
  ## YOUR CODE HERE ##
  return

Part 5) Repeat and create a function on which SA outperforms the other two.

def sa_best(x, y):
  ## YOUR CODE HERE ##
  return

Part 6) Repeat and create a function on which LBS outperforms the other two.

def lbs_best(x, y):
  ## YOUR CODE HERE ##
  return

Part 7) Explain your approach to parts 4-6

Describe how you constructed each function and your hypothesis as to why they lead to the relative behaviors of the optimization methods.