Inspired Algorithms

Greedy randomized adaptive search procedure: Step-wise construction with local search

2009-02-02T21:30:00.002+11:00

The Iterated Local Search algorithm introduced the combination of a biased global search algorithm followed by a local search to refine the results. This general problem solving pattern is also the basis of the Greedy Randomized Adaptive Search Procedure (GRASP).

The GRASP involves alternating a greedy adaptive solution construction procedure with a local search algorithm to refine the coarser search result. The approach was designed for combinatorial problems where the individual components of a candidate solution can be discriminated in the context of the goal (such as cost). This greedy step-wise construction involves the preparation of Restricted Candidate List (RCL) of components from which to select from, and the random selection of a component from the RCL to add to the growing solution. This process continues until a viable candidate solution is constructed and can be passed onto a local search procedure for evaluation.

The power of GRASP is in the preparation of the RCL. Typically the value of individual components to a final solution is used to restrict entry to the list (a threshold). Alternative approaches involve a fixed RCL size (of the top components), and an adaptive threshold or size over time. Management of the RCL requires careful attention, ensuring sufficient flexibility in the selection of components for a solution, whilst biasing the search towards a global optimum supplemented with a local search technique. It is interesting to note that the use of a strict RCL that selects the best component during each step of construction will result in a greedy search.

Greedy randomized adaptive search procedure
This guide provide a demonstration of the GRASP applied to a standard instance of the Traveling Salesman Problem. This classical combinatorial optimization provides a good example for the step-wise construction technique used by the GRASP.

The Berlin52TSP class provides an implementation of a small sized TSP with 52 cities. The problem data is stored as constants within the class. The COORDINATES constant provides a 2D array of city coordinates and the OPTIMAL_TOUR constant stores a representation of the optimal tour (permutation) of cities that used in the algorithms termination criterion. The initialize() constructor in turn calls the build_distance_matrix method to pre-calculate a matrix of the Euclidean distance between all combinations of cities. The city-to-city distance calculation represents the most expensive aspect of candidate solution evaluations which is eased via the pre-calculation requiring only a lookup of city ID's to determine the intervening distance. The evaluate(permutation) makes use of this pre-calculation via the dist(c1, c2) method, evaluating candidate solutions in the form a 51-city permutation where each city ID is specified between 1 and 52 inclusively.

class Berlin52TSP
  OPTIMAL_TOUR = [1,49,32,45,19,41,8,9,10,43,33,51,11,52,14,13,47,26,
    27,28,12,25,4,6,15,5,24,48,38,37,40,39,36,35,34,44,46,16,29,50,20,
    23,30,2,7,42,21,17,3,18,31,22]
        
  COORDINATES = [[565, 575],[25, 185],[345, 750],[945, 685],[845, 655],
  [880, 660],[25, 230],[525, 1000],[580, 1175],[650, 1130],[1605, 620], 
  [1220, 580],[1465, 200],[1530, 5],[845, 680],[725, 370],[145, 665],
  [415, 635],[510, 875], [560, 365],[300, 465],[520, 585],[480, 415],
  [835, 625],[975, 580],[1215, 245],[1320, 315],[1250, 400],[660, 180],
  [410, 250],[420, 555],[575, 665],[1150, 1160],[700, 580],[685, 595],
  [685, 610],[770, 610],[795, 645],[720, 635],[760, 650],[475, 960],
  [95, 260],[875, 920],[700, 500],[555, 815],[830, 485],[1170, 65],
  [830, 610],[605, 625],[595, 360],[1340, 725],[1740, 245]]
  
  attr_reader :num_cities

  def initialize()
    @num_cities = COORDINATES.length        
    @distance_matrix = Array.new(@num_cities) {Array.new(@num_cities)}
    build_distance_matrix
    @optimal_tour_length = evaluate(OPTIMAL_TOUR) # calculate
  end
  
  def build_distance_matrix
    # symmetrical matrix along the diag
    puts "pre-calculating the TSP distance matrix"
    @distance_matrix.each_with_index do |row, i|
      row.each_index do |j|
        row[j] = euc_2d(COORDINATES[i], COORDINATES[j])
      end
    end
  end

  def evaluate(permutation)
    dist = 0    
    permutation.each_with_index do |c1, i|
      c2 = (i==@num_cities-1) ? permutation[0] : permutation[i+1] 
      dist += dist(c1,c2)
    end
    return dist
  end
  
  # expects city numbers [1,52]
  def dist(c1, c2)
     @distance_matrix[c1-1][c2-1]
  end
  
  def euc_2d(c1, c2)
    # As defined in TSPLIB'95 (EUC_2D)
    Math::sqrt((c1[0] - c2[0])**2 + (c1[1] - c2[1])**2).round
  end

  def is_optimal?(scoring)
    scoring == optimal_score
  end

  def optimal_score
    @optimal_tour_length
  end
  
  # true if s1 is better score than s2
  def is_better?(s1, s2)
    s1 < s2 # minimizing
  end
end

The Solution class provides a generic container for candidate solutions. Problem specific solution data (in this case permutations) are stored in the immutable @data instance variable, and a mutable solution costing is stored in the @score instance variable.

class Solution
  attr_reader :data
  attr_accessor :score
  
  def initialize(data)
    @data = data
    @score = 0.0/0.0 # NaN
  end
  
  def to_s
    "[#{@data.inspect}] (#{@score})"
  end    
end

The GRASP class provides an implementation of the Greedy Randomized Adaptive Search Procedure. The initialize(max_iterations) constructor for the algorithm sets up the maximum number of iterations for the algorithms termination criteria, as well as the algorithm specific @max_rcl_size that defines the maximum size of the Restricted Candidate List each step during construction, defaulted to 5 components. The search(problem) method provides the entry point to the search that executes the main loop until the termination criterion should_stop?(curr_it, problem) returns true. The main loop involves first a call to greedy_randomized_solution(current, problem) to construct a candidate solution, followed by a call to local_search_solution(candidate, problem) to refine the result.

The greedy_randomized_solution(solution, problem) provides a simplistic realization of the GRASP with a fixed size RCL. The construction starts with a random city and proceeds by first preparing a set of viable (unused) cities to which to connect. The candidates are ordered by their distance from the 'current' city, the list is trimmed to be no longer than the @max_rcl_size instance variable, and finally a random component from the list is selected and appended to the permutation.

The local search algorithm defined in the local_search_solution(solution, problem) method involves the application of the classical TSP local search called 2-opt (the two_opt_solution(solution) method) where two cities are selected and reconnected. This procedure is repeated 15 times as an approximation of locating the local optimum for the provided candidate solution.

class GRASP
  attr_accessor :max_iterations, :max_rcl_size
  attr_reader :best_solution
  
  def initialize(max_iterations)
    @max_iterations = max_iterations  
    @max_rcl_size = 5 # default  
  end
  
  # execute the algorithm
  def search(problem)
    # initialize the search
    @best_solution = nil
    current = generate_initial_solution(problem)
    evaluate_candidate_solution(current, problem)
    curr_it = 0
    begin
      # greedy randomized solution
      candidate = greedy_randomized_solution(current, problem)  
      evaluate_candidate_solution(candidate, problem) # eval    
      # local search
      candidate = local_search_solution(candidate, problem)
      # greedy acceptance
      current = candidate if problem.is_better?(candidate.score, current.score)
      curr_it += 1
    end until should_stop?(curr_it, problem)
    return @best_solution
  end
  
  def should_stop?(curr_it, problem)
    (curr_it >= @max_iterations) or problem.is_optimal?(best_solution.score)
  end
  
  def generate_initial_solution(problem)
    all = Array.new(problem.num_cities) {|i| (i+1)}
    permutation = Array.new(all.length) {|i| all.delete_at(rand(all.length))}
    return Solution.new(permutation)
  end

  def greedy_randomized_solution(solution, problem)    
    # construct a valid perm
    perm = []
    perm << rand(problem.num_cities) + 1 # random starting point
    last = perm[0]
    while perm.length < problem.num_cities
      # create restricted candidate list for components
      all = Array.new(problem.num_cities) {|i| (i+1)}
      # ensure we can only add unused components
      rcl = all - perm 
      # order by distance, asc (best have smallest distance)
      rcl.sort! {|i,j| problem.dist(last, i) <=> problem.dist(last, j)}
      # trim to fixed size
      rcl.pop until rcl.length <= @max_rcl_size
      # select random
      last = rcl[rand(rcl.length)]
      perm << last
    end    
    return Solution.new(perm)
  end

  def local_search_solution(solution, problem)
    # greedy iterated 2-opt
    15.times do
      candidate = two_opt_solution(solution)
      evaluate_candidate_solution(candidate, problem)
      if problem.is_better?(candidate.score, solution.score)
        solution = candidate 
      end
    end
    return solution
  end

  def two_opt_solution(solution)
    perm = Array.new(solution.data) # copy
    # select a sub-sequence
    c1, c2 = rand(perm.length), rand(perm.length)
    c2 = rand(perm.length) while c1 == c2
    # ensure c1 is low and c2 is high
    c1, c2 = c2, c1 if c2 < c1
    # reverse sub-sequence
    perm[c1...c2] = perm[c1...c2].reverse
    return Solution.new(perm)
  end

  def evaluate_candidate_solution(solution, problem)
    solution.score = problem.evaluate(solution.data)
    # keep track of the best solution found
    if @best_solution.nil? or
      problem.is_better?(solution.score, @best_solution.score)
      @best_solution = solution
      puts "> new best: #{solution.score}"               
    end
  end
end

The algorithm can be executed by first creating an instance of the problem, the algorithm, and passing the problem to the algorithms search method. The algorithm will displaying its progress during the search and return the best result found. The global random number generated is seeded with a fixed number to ensure a consistent sequence of random numbers are used which can prove useful during testing.

srand(1) # set the random number seed to 1
algorithm = GRASP.new(500) # limit to 500 iterations 
problem = Berlin52TSP.new # create a problem
best = algorithm.search(problem) # execute the search
puts "Best Solution: #{best}" # display the best solution

The demonstrated GRASP provides much opportunity for extension. One may modify the implementation of the TSP problem instance to provide a class that can load problems form definition files provided by the TSPLIB. The size of the RCL list may be tuned to vary the greediness of solution construction. Further, the preparation of the RCL may be modified to make use of the current working solution passed into the function, perhaps using existing selected components to bias the selection of those components during construction. Finally, adaptive RCL management procedures may be adopted as well as alternative local search procedures for the TSP.

Source Code
Download: GRASP.rb

Further Reading

Papers

T. A. Feo and M. G. C. Resende. Greedy randomized adaptive search procedures. (pre-print) Journal of Global Optimization, 6, 109-133, 1995.
M. G. C. Resende and C. C. Ribeiro Greedy randomized adaptive search procedures (pre-print). In F. Glover and G. Kochenberger, editors, Handbook of Metaheuristics, Springer, pages. 219-249, 2003.
L. Pitsoulis and M. G. C. Resende Greedy randomized adaptive search procedures (pre-print). In P. M. Pardalos and M. G. C. Resende, editors, Handbook of Applied Optimization, Oxford University Press, pages. 168-181, 2002

Sites

Greedy randomized adaptive search procedure. Wikipedia
TSPLIB

Have you found a bug? Know of a great reference? Got an idea for a guide? Please send an email to jasonb@InspiredAlgorithms.com.

Need Help? Want to discuss this post? Please visit the Inspired Algorithms User Group.

Tabu Search: Using explicit memory

2009-01-19T11:00:00.053+11:00

A problem with local search techniques is that after locating a local optima, they will continue to re-sample the same optima and neighboring solutions in the search space over and over again until a stopping condition is triggered. Tabu Search is an algorithm to address this problem by adding memory to a search algorithm, preventing or penalizing the algorithm for re-sampling (cycling) the same solutions and in turn encouraging further exploration of the search space.

The core search algorithm embedded within a tabu search generates candidate solutions in the neighborhood of a current working solution. A tabu list is maintained that holds a set of recently visited neighboring solutions (or components thereof) that may not be visited again (they are taboo). Moves are held in the tabu list for a defined tabu tenure that is typically fixed, proportional to the problem size, or generated randomly within a pre-defined range.

The tabu list is referred to as the short term memory of a search algorithm. A searches intermediate term memory refer to the factors that direct the process to promising areas (expected payoff) of the search space inferred from past experience, referred to as intensification. Examples of such intermediate memory processes include the inclusion of common components from recently visited good quality solutions, adapting the tabu tenures when locally optimal solutions are located, and restarting the search (clearing the tabu list) with a locally optimal solution. Long term memory processes are used to counter the intensification of intermediate memory by encouraging exploration in the search space, referred to as diversification. Examples include the use of penalties per component that are increased with the frequency of inclusion, and changing the acceptance criteria to allow non-improving moves.

Finally, so-called aspiration criteria may be used to permit neighbourhood moves on the tabu list (overcome the tabu). A standard aspiration criteria is to accept a tabu move if is better than a the current working position or better than the best solution found by the search to date.

Tabu Search for the Traveling Salesman Problem
This guide applies the tabu search to a classical combinatorial optimization problem called the Traveling Salesman Problem (TSP). The general class of problem involves finding an order to visit a set of cities to minimize the overall distance traveled.

The Berlin52TSP class provides a definition of a TSP instance called Berlin 52 (with 52 nodes or cities). The class is large given that it has the city coordinates hard coded in the COORDINATES class constant, as well as the optimal permutation in the OPTIMAL_TOUR constant. The evaluate(permutation) expects a valid permutation of 52 non-repeating cities to visit as an array of integers (between 1 and 52 inclusively) and is evaluated as the rounded Euclidean distance of the cycle.

class Berlin52TSP
  OPTIMAL_TOUR = [1,49,32,45,19,41,8,9,10,43,33,51,11,52,14,13,47,26,
    27,28,12,25,4,6,15,5,24,48,38,37,40,39,36,35,34,44,46,16,29,50,20,
    23,30,2,7,42,21,17,3,18,31,22]
        
  COORDINATES = [[565, 575],[25, 185],[345, 750],[945, 685],[845, 655],
  [880, 660],[25, 230],[525, 1000],[580, 1175],[650, 1130],[1605, 620], 
  [1220, 580],[1465, 200],[1530, 5],[845, 680],[725, 370],[145, 665],
  [415, 635],[510, 875], [560, 365],[300, 465],[520, 585],[480, 415],
  [835, 625],[975, 580],[1215, 245],[1320, 315],[1250, 400],[660, 180],
  [410, 250],[420, 555],[575, 665],[1150, 1160],[700, 580],[685, 595],
  [685, 610],[770, 610],[795, 645],[720, 635],[760, 650],[475, 960],
  [95, 260],[875, 920],[700, 500],[555, 815],[830, 485],[1170, 65],
  [830, 610],[605, 625],[595, 360],[1340, 725],[1740, 245]]
  
  attr_reader :num_cities

  def initialize()
    @num_cities = COORDINATES.length        
    @optimal_tour_length = evaluate(OPTIMAL_TOUR) # calculate
  end
  
  def evaluate(permutation)
    dist = 0    
    permutation.each_with_index do |c1, i|
      c2 = (i==@num_cities-1) ? permutation[0] : permutation[i+1] 
      dist += euc_2d(COORDINATES[c1-1], COORDINATES[c2-1])
    end
    return dist
  end
  
  def euc_2d(c1, c2)
    # As defined in TSPLIB'95 (EUC_2D)
    Math::sqrt((c1[0] - c2[0])**2 + (c1[1] - c2[1])**2).round
  end

  def is_optimal?(scoring)
    scoring == optimal_score
  end

  def optimal_score
    @optimal_tour_length
  end
  
  # true if s1 is better score than s2
  def is_better?(s1, s2)
    s1 < s2 # minimizing
  end
end

The Solution class provides a container for candidate solutions, with an immutable @data instance variable for problem specific solution representation and a mutable @score instance variable to hold its costing. The overridden to_s method provides a pretty string representation of the solution for display on the command line. Important to this example is the overridden ==(other) method that is used to test equality between solution objects. Here, the method is specialized assuming the @data holds an array, and compares two Solution objects based on the equality of the permutations they hold. Permutation comparison is not dependent on the order the permutation is recorded, so equality is checked in both directly and with one of the permutations reversed. This equality method is used to assess whether a solution has already been added to the tabu list (later on). Finally, a mutable @count instance variable is maintained to use in the frequency a tabu solution is desired (later).

class Solution
  attr_reader :data
  attr_accessor :score, :count
  
  def initialize(data)
    @data = data
    @score = 0.0/0.0 # NaN
    @count = 0
  end
  
  def ==(other)
    # permutation is direction independant
    @data == other.data or @data.reverse == other.data
  end
  
  def to_s
    "[#{@data.inspect}] (#{@score})"
  end    
end

The TabuSearchAlgorithm class offers an implementation of a iterative 2-opt search algorithm with a tabu list and aspiration criteria for a given TSP instance. The initialize(max_iterations) constructor records the maximum number of iterations the algorithm will execute and sets additional tabu search specific configuration. The @tabu_tenure instance variable defines the the size of the tabu list under the assumption that one candidate is added to the list each iteration. The @aspiration_frequency instance variable defines the maximum number of requests that may be made for a tabu solution on the tabu list before it may be reused (prematurely liberated from the tabu list).

The search(problem) provides the entry point for the search for a given TSP instance. A random starting point is used to initialize the search provided by the generate_initial_solution(problem) method. The algorithm loop generally involves generating neighbouring candidate solutions one at a time, evaluating a solution, and using it as the current working position greedily if it has an improved problem specific scoring. Neighboring candidate solutions are generated by calling the two_opt_solution(solution) method for a given current working point. This standard local search for the TSP was chosen for its speed, effectiveness, and simplicity and involves selecting two break points in a permutation and reversing the sub-sequence between the points.

The generation of neighboring candidate solutions loops until a candidate solution is generated that is not already on the tabu list. Importantly, this loop increments the number of times each tabu solution is requested, liberating a specific candidate solution if it's request frequency exceeds @aspiration_frequency, providing a basic aspiration criteria. The tabu list is maintained at the bottom of the main algorithm loop, first adding the current candidate solution regardless of whether it is accepted or not, and then whittling down the tabu_list to @tabu_tenure when required (treating the tabu_list array like a first-in-last-out FILO queue).

class TabuSearchAlgorithm
  attr_accessor :max_iterations
  attr_reader :best_solution
  
  def initialize(max_iterations)
    @max_iterations = max_iterations
    @tabu_tenure = 100 # maximum number of solutions on the list
    @aspiration_frequency = 5 # max tries to visit a taboo before allowing
  end
  
  # execute a search on the provided problem
  def search(problem)
    tabu_list = [] # prep the tabu list
    # random starting point
    current = generate_initial_solution(problem)
    tabu_list << current # make candidate taboo
    evaluate_candidate(current, problem)
    curr_it = 0       
    begin
      # generate new neighbour that is not tabo
      selection_ok = true
      begin
        selection_ok = true # optimism
        candidate = two_opt_solution(current) # generate
        # check tabu
        if tabu_list.include?(candidate)
          other = tabu_list.find {|t| t==candidate}
          if (other.count+=1) >= @aspiration_frequency
            puts "> we have an aspiration!!!"
            tabu_list.delete_if {|t| t==candidate} # ensure no duplicates
          else
            puts "> rejected tabu"
            selection_ok = false 
          end
        end
      end until selection_ok
      evaluate_candidate(candidate, problem) # evaluate
      # manage short term memory
      tabu_list << candidate # make candidate taboo
      tabu_list.delete_at(0) while tabu_list.length > @tabu_tenure
      # greedy acceptance
      current = candidate if problem.is_better?(candidate.score, current.score)
      curr_it += 1
    end until should_stop?(curr_it, problem)
    return @best_solution
  end
  
  def should_stop?(curr_it, problem)
    (curr_it >= @max_iterations) or problem.is_optimal?(best_solution.score)
  end
  
  def generate_initial_solution(problem)
    all = Array.new(problem.num_cities) {|i| (i+1)}
    permutation = Array.new(all.length) {|i| all.delete_at(rand(all.length))}
    return Solution.new(permutation)
  end

  def two_opt_solution(solution)
    perm = Array.new(solution.data) # copy
    # select a sub-sequence    
    c1, c2 = rand(perm.length), rand(perm.length)
    c2 = rand(perm.length) while c1 == c2
    # ensure c1 is low and c2 is high
    c1, c2 = c2, c1 if c2 < c1
    # reverse sub-sequence
    perm[c1...c2] = perm[c1...c2].reverse
    return Solution.new(perm)
  end

  def evaluate_candidate(solution, problem)
    solution.score = problem.evaluate(solution.data)
    # keep track of the best solution found
    if @best_solution.nil? or 
      problem.is_better?(solution.score, @best_solution.score)
      @best_solution = solution
      puts " > new best: #{solution.score}"               
    end
  end
end

The algorithm is demonstrated by first initializing the global random number generator to 1 to ensure the same sequence of random numbers on each execution (useful for testing). The problem is created and algorithm is initialized with a maximum of 2000 iterations. The search is executed by passing the problem instance to the algorithms search function, and after the execution of the run the best solution is displayed.

srand(1) # set the random number seed to 1
algorithm = TabuSearchAlgorithm.new(2000) # limit to 2000 iterations 
problem = Berlin52TSP.new # create a problem
best = algorithm.search(problem) # execute the search
puts "Best Solution: #{best}" # display the best solution

During the execution of the algorithm, the scoring of the best solution is displayed each time a new best is located. The loop of the main algorithm displays a message each time a tabu solutions aspiration criteria is triggered or when a generated candidate solution is denied based on its existence in the tabu list. You will observe that as the search gets closer to the global optimum, the number of candidates denied increases and a number of liberation events occur as the aspiration criteria are triggered.

There is plenty of room for extension on this example. One may change the configuration of the tabu search, increasing or decreasing both the tabu tenure as well as the aspiration frequency. Also experimentation of the search algorithm with and without a tabu list may be interesting to see if and when the short term memory of this tabu search example are useful. Alternative representations may be used for the tabu list such as the selection of cities to swap in the local search. Alternatively, the local search may be adjusted to construct tours in a step-wise (per-city) manner and maintain tabu information per-edge in an adjacency matrix. This matrix would then discourage or prevent the use of specific edges during solution construction. Finally, additional memory processes may be introduced such as search-restarts with the best found solution for intermediate memory, and the solution acceptance criteria may be made less greedy to allow non-improving moves to be taken.

Source Code
Download: TabuSearch.rb

Further Reading

Papers:

Fred Glover and Manuel Laguna, Tabu Search. Springer, 1997
F. Glover, Tabu Search: Part I. (pre-print) INFORMS Journal on Computing 1(3), pages 190-206, 1989
F. Glover, Tabu Search: Part II. (pre-print) INFORMS Journal on Computing 2(1), pages 4-32, 1990
Fred Glover and Eric Taillard, A user's guide to tabu search. Annals of Operations Research 41(1), pages 1-28, 1993
F. Glover, Tabu Search: A Tutorial. (pre-print) INTERFACES 20(4), pages 74-94, 1990
Tabu Search Applied to TSP. In Johannes Josef Schneider and Scott Kirkpatrick, Stochastic Optimization, Springer, pages 441-447. 2006

Sites:

Tabu Search. Wikipedia
Tabu Search. Wolfram MathWorld
TSPLIB

Multiple Restart Search: Multiple runs to address convergence

2009-01-17T11:00:00.040+11:00

While searching for the global optimum of a given problem domain the may question arise as to whether it is worth continuing with the search or starting the search again. This strategy of repeatedly executing a search on a problem domain with new initial conditions is called Random, Early, or Multiple Restart Search.

The Multiple Restart Search algorithm involves the application of another search technique until a predefined condition is met, after which the embedded search technique is executed again with different starting conditions. The conditions that govern when a restart will occur are refereed to as the restart schedule and examples include: a fixed number of function evaluations, a fixed time, or the convergence of the search process (no improvements for a fixed number of iterations or function evaluations).

Restarting or multiple executions of a technique is a common practice with stochastic optimization algorithms as they are influenced by their initial starting conditions and convergent behavior that may lead to local rather than the desired global optimum of a problem's search space. Restart schedules are also used with local search and hill climbing searching algorithms as they allow potentially multiple local optima to be explored improving the algorithms approximation of the global optima. Such algorithms are commonly referred to a Random-Restart Hill Climbing algorithms.

The multiple executions of an embedded search algorithm (the restarts) are traditionally executed sequentially as the name suggests. Alternatively, the algorithm may exploit parallel or distributed hardware and execute each algorithm run in parallel, offering a linear (and in some cases a super-linear) speedup of locating or approximating a globally optimal solution.

Multiple Restart Search for Function Optimization
This guide provides an example of a the multiple restart search algorithm applied to a continuous function optimization problem. The embedded search algorithm that is repeatedly restarted is a hill climbing algorithm, and as such the generalized approach may be referred to as a Random-Restart Hill Climber.

The SquaringFunction class defines an instance of a function optimization problem. The problem has a configurable number of parameters (@num_dimensions) the squares of which are summed together (hence the name). The evaluate(vector) method expects a vector of real valued numbers, the in_bounds?(vector) convenience method determines whether or not a given vector is valid within the bounds of the problems search space [-5.12, 5.12], and the is_optimal?(scoring) method determines whether or not the given scoring is optimal (equal to zero).

class SquaringFunction
  attr_reader :dimensions, :min, :max

  def initialize(dimensions=2)
    @dimensions = dimensions
    @min, @max = -5.12, +5.12
  end

  def evaluate(vector)
    vector.inject(0) {|sum, x| sum + (x ** 2.0)}
  end  
  
  def in_bounds?(vector)
    vector.each {|x| return false if x>@max or x<@min}
    return true    
  end

  def is_optimal?(scoring)
    scoring == optimal_score
  end

  def optimal_score
    0.0
  end
  
  # true if s1 has a better score than s2
  def is_better?(s1, s2)
    s1 < s2 # minimizing
  end
end

The Solution class provides a container for candidate solutions. The @data instance variable maintains an immutable solution representation (in this case real valued vectors), whereas the @score instance variable maintains a mutable problem specific scoring for the maintained data.

class Solution
  attr_reader :data
  attr_accessor :score
  
  def initialize(data)
    @data = data
    @score = 0.0/0.0 # NaN
  end
  
  def to_s
    "[#{@data.inspect}] (#{@score})"
  end    
end

The MultipleRestartSearchAlgorithm provides an implementation of the random-restart hill climbing algorithm. The initialize(max_restarts) constructor stores some algorithm configuration parameters in instance variables, specifically the maximum number of restarts provided as an argument, and the number of non-improving steps taken per run of the embedded hill climbing algorithm (set to 100). The search(problem) method manages the restart schedule generating a random starting point which is used to seed a run of the embedded hill climbing algorithm on each restart.

The problem independent hill_climb(problem, current) method does the hard work for the search, performing a fixed step hill climbing procedure on a provided starting point. The step size is set to 10% of the possible search space, and is unchanged throughout the run. The main loop of this procedure involves the generation of two steps (a positive and negative application of the step size) for a single problem component (function parameter). If one of the two steps results in an improvement the candidate solution is taken as the current position otherwise the number of non-improving iterations is incremented. This counter is rest each time an improving move is made, and is used to measure the convergence of a run of the hill climbing algorithm. The procedure loops through all components of the problem (one per iteration) seeking improvements. A run completes on convergence, which is occurs when the no_improve_count variable exceeds the @max_no_improvements parameter.

The take_step(problem, current, step_size, index, add) method creates a single new candidate solution from a given solution by mutating a specific component of the solution in a specific direction. Interestingly, this is the only problem-specific part of the algorithm, that in turn calls the next_bfloat(min, max) utility method to generates random values that are used to perturb a given vector.

class MultipleRestartSearchAlgorithm
  attr_accessor :max_iterations
  attr_reader :best_solution
  
  def initialize(max_restarts)
    @max_restarts = max_restarts
    @max_no_improvements = 100 # max steps for hill climber
  end
  
  # execute a random search on the provided problem
  def search(problem)    
    curr_it = 0
    begin
      # generate a random solution 
      current = generate_random_solution(problem)
      evaluate_candidate_solution(current, problem)
      # hill climb until convergence
      run_best = hill_climb(problem, current)
      curr_it += 1      
      puts "> finished run #{curr_it} with: #{run_best.score}"   
    end until should_stop?(curr_it, problem)
    return @best_solution
  end
  
  def should_stop?(curr_it, problem)
    (curr_it >= @max_restarts) or problem.is_optimal?(best_solution.score)
  end
  
  def generate_random_solution(problem)
    real_vector = Array.new(problem.dimensions) do
      next_bfloat(problem.min, problem.max)
    end
    return Solution.new(real_vector)
  end
  
  def hill_climb(problem, current)
    step_size = (problem.max-problem.min)*0.10
    no_improve_count = 0
    index = 0
    begin
      # take a step
      step1 = take_step(problem, current, step_size, index, true)
      evaluate_candidate_solution(step1, problem)
      step2 = take_step(problem, current, step_size, index, false)
      evaluate_candidate_solution(step2, problem)      
      # check for improvement
      if problem.is_better?(step1.score, current.score) or 
        problem.is_better?(step2.score, current.score)
        # store the best
        if problem.is_better?(step1.score, current.score)
          current = step1
        else
          current = step2
        end
        no_improve_count = 0 # reset
      else
        no_improve_count += 1 # count consecutative no improvements
      end
      index = (index==problem.dimensions-1) ? 0 : index+1
    end until no_improve_count >= @max_no_improvements
    return current
  end
  
  def take_step(problem, current, step_size, index, add)
    vector = nil
    begin # keep stepping until a valid point is generated
      offset = next_bfloat(0, step_size)
      vector = Array.new(current.data)
      vector[index] += (add) ? offset : -offset
    end until problem.in_bounds?(vector)
    return Solution.new(vector)
  end

  def next_bfloat(min, max)
    min + ((max - min) * rand)
  end
  
  def evaluate_candidate_solution(solution, problem)
    solution.score = problem.evaluate(solution.data)
    # keep track of the best solution found
    if @best_solution.nil? or
      problem.is_better?(solution.score, @best_solution.score)
      @best_solution = solution                  
    end
  end  
end

The algorithm is demonstrated by first seeding the global random number generator to ensure a consistent sequence of random numbers (useful for testing) and creating instance of the problem and algorithm. The search is executed by passing the problem instance to the algorithm's search method that displays the status of the algorithms as the best result found after each restart. The demonstration ends by displaying an approximation of the problem's global optima (which is 0.0 in each dimension) as the best solution encountered by the algorithm.

srand(1) # set the random number seed to 1
algorithm = MultipleRestartSearchAlgorithm.new(10) # limit to 10 restarts
problem = SquaringFunction.new(5) # create a problem with 5 dimensions
best = algorithm.search(problem) # execute the search
puts "Best Solution: #{best}" # display the best solution

This general procedure is called a meta-search algorithm as it can generally be applied to any other direct search technique. The demonstration may be extended by elaborating upon the hill climbing algorithm such as using an adaptive rather than a fixed step size. Further, the embedded hill climbing (local search) algorithm may be replaced with a population-based global stochastic search algorithm. Additionally, restarted algorithm runs may exploit the results found by previous runs, and alternate restart schedules may be used such as a maximum number of function evaluations or relative improvement over previous runs.

Source Code
Download: MultipleRestartSearch.rb

Further Reading

Papers:

X. Hu, R. Shonkwiler, and M. Spruill. Random restart in global optimization. Technical Report 110592-015, Georgia Tech School of Mathematics, 1994.
M. Magdon-Ismail and A. Atiya. The early restart algorithm. Neural Computation, 12(6), pages 1303–1312, 2000.
A. Fukunaga. Restart scheduling for genetic algorithms. (pre-print) In Parallel Problem Solving from Nature - PPSN V, pages 357-366. 1998.
Marco Muselli, A Theoretical Approach to Restart in Global Optimization. (pre-print) In Journal of Global Optimization, 10(1), pages 1-16. 1997

Sites:

Hill climbing. Wikipedia

Adaptive Random Search: Varying step size based on performance

2009-01-15T11:00:00.003+11:00

The Random Search algorithm involves the blind generation and evaluation of candidate solutions to a problem, where as the Localized Random Search algorithm makes the simple addition of successively varying the best solution found by the search so far. The Adaptive Random Search algorithm, like Localized Random Search is concerned with searching the neighborhood of the best solution found by the search so far, although unlike that algorithm the definition of neighborhood is not fixed. Instead it is changed dynamically over the course of the search.

When exploring the neighborhood of a candidate solution, the size of the steps taken within the neighborhood may be fixed in what is referred to as Fixed Step Size Random Search. The Adaptive Random Search algorithm involves varying the step size during the search in response to the relative improvements or non-improving steps as defined by the objective function evaluation. A small number of samples are taken each iteration with varied small and large step sizes in an attempt to approximate an optimal step size to converge the search to the local optimum. The step size that generate an improved candidate solution is adopted as the current step size for the search, and improved solutions are taken as the current working position. When a given step size stops producing improvements, its size is reduced to encourage the search algorithm to more carefully explore the smaller localized neighborhood. Larger step sizes are frequently explored in an effort to potentially escape local optima, and very large step sizes are sampled with a low frequency in an effort to explore distant regions of the search space.

Adaptive Random Search for Continuous Function Optimization
This guide provides an example of the Adaptive Random Search algorithm applied to an instance of the continuous function optimization problem.

The ExponentFunction class defines an instance of a function optimization problem with a variable number of function parameters or dimensions. The evaluate(vector) method expects a vector of real values in the problem bounds and evaluates it as the sum of the fourth power of each function parameter. The global optimum of this minimization function is located at 0.0 which is accessible via the optimal_score method and is used by is_optimal?(scoring) to test whether the algorithm has located to the global optimum. Finally, the problem provides a in_bounds?(vector) method for assessing whether a given vector falls within the hypercube of the valid problem space, useful for algorithms like Adaptive Random Search that generate new samples from existing samples and may stray out of the problem's bounds.

class ExponentFunction
  attr_reader :dimensions, :min, :max

  def initialize(dimensions=2)
    @dimensions = dimensions
    @min, @max = -5.12, +5.12
  end

  def evaluate(vector)
    vector.inject(0) {|sum, x| sum + (x ** 4.0)}
  end  
  
  def in_bounds?(vector)
    vector.each {|x| return false if x>@max or x<@min}
    return true    
  end

  def is_optimal?(scoring)
    scoring == optimal_score
  end

  def optimal_score
    0.0
  end
  
  # true if s1 has a better score than s2
  def is_better?(s1, s2)
    s1 < s2 # minimizing
  end
end

The Solution class provides a convenience container for holding immutable real-valued vectors in the @data instance variable, and a mutable @score instance variable for holing the costing of the vector assigned by a problems objective function.

class Solution
  attr_reader :data
  attr_accessor :score
  
  def initialize(data)
    @data = data
    @score = 0.0/0.0 # NaN
  end
  
  def to_s
    "[#{@data.inspect}] (#{@score})"
  end    
end

The AdaptiveRandomSearchAlgorithm class provides a reasonably generic implementation of the algorithm. The initialize(max_iterations) constructor expects a maximum number of algorithm iterations as an argument and goes on to define a set of algorithm configuration parameters as class instance variables. The search(problem) method is large and contains all the problem agnostic conditional logic of the algorithm. The search starts by generating a random starting point for the search (current) and setting an initial step size of 10% of the size of the search domain. The main loop of the algorithm is concerned with first generating a new sample using the step_size and a second sample using a larger step size. The larger step size is fixed at @small_factor * step_size although periodically (once every @large_factor_multiple iterations) a much larger factor (@large_factor) of the current step size is used.

If one of the generated steps results in a candidate solution better than the current solution it is replaced. If the larger step size resulted in the improvement, it is used the current step size. If neither generated step sizes resulted in an improved candidate solution, a count of non-improving steps is incremented. If this count exceeds @maximum_no_improvements then the current step size is reduced by a constant factor (step_size / @small_factor).

The generate_random_solution(problem) method used for the starting point of the search generates a random vector in the bounds of the problem space. The take_step(problem, current, step_size) is responsible for stochastically generating a new solution from a provided existing solution within a neighborhood defined by the provided step_size. This works by generating a random vector between +/- step_size which is then added to the current sample. This means the neighbourhood is a step_size sized hypercube around the origin. Steps are taken by this method until a valid (within the problem bounds) candidate solution is generated.

class AdaptiveRandomSearchAlgorithm
  attr_accessor :max_iterations
  attr_reader :best_solution
  
  def initialize(max_iterations)
    @max_iterations = max_iterations
    @small_factor = 1.3 # normal step size incrase(*) and decrease(/) rate
    @large_factor = 10 # step size factor increase for larger jump
    @large_factor_multiple = 100 # max steps before larger jump is tried
    @maximum_no_improvements = 50 # max non-improvement steps before size decreased
  end
  
  # execute a random search on the provided problem
  def search(problem)    
    # random starting point
    current = generate_random_solution(problem)
    evaluate_candidate_solution(current, problem)
    curr_it = 0
    step_size = (problem.max-problem.min)*0.1 # 10% of the domain
    no_change_cnt = 0
    begin
      # current step
      step = take_step(problem, current, step_size)
      evaluate_candidate_solution(step, problem)
      # take second larger step
      factor = @small_factor
      factor = @large_factor if curr_it.modulo(@large_factor_multiple) == 0
      larger_step_size = step_size * factor
      larger_step = take_step(problem, current, larger_step_size)
      evaluate_candidate_solution(larger_step, problem)
      # check for improvement
      if problem.is_better?(step.score, current.score) or 
        problem.is_better?(larger_step.score, current.score)        
        # select new step size if larger step is better
        if problem.is_better?(larger_step.score, step.score)
          step_size = larger_step_size # increase step size
          current = larger_step # new best solution
        else
          current = step # new best solution
        end
        no_change_cnt = 0 # reset counter
      # check for step size decrease
      elsif (no_change_cnt+=1) >= @maximum_no_improvements
        no_change_cnt = 0 # reset counter
        step_size /= @small_factor # decrease step size
      end
      curr_it += 1      
    end until should_stop?(curr_it, problem)
    return @best_solution
  end
  
  def should_stop?(curr_it, problem)
    (curr_it >= @max_iterations) or problem.is_optimal?(best_solution.score)
  end
  
  def generate_random_solution(problem)
    real_vector = Array.new(problem.dimensions) do
      next_bfloat(problem.min, problem.max)
    end
    return Solution.new(real_vector)
  end
  
  def take_step(problem, current, step_size)
    vector = nil
    begin # keep stepping until a valid point is generated
      step = Array.new(problem.dimensions) do
        v = next_bfloat(-step_size, +step_size) 
      end
      data = current.data
      vector = Array.new(data.length) {|i| data[i] + step[i]}
    end while !problem.in_bounds?(vector)
    return Solution.new(vector)
  end

  def next_bfloat(min, max)
    min + ((max - min) * rand)
  end
  
  def evaluate_candidate_solution(solution, problem)
    solution.score = problem.evaluate(solution.data)
    # keep track of the best solution found
    if @best_solution.nil? or
      problem.is_better?(solution.score, @best_solution.score)
      @best_solution = solution
      puts "> new best: #{solution.score}"               
    end
  end  
end

The algorithm is demonstrated by first seeding the global random number number generator to 1 to ensure the same sequence of random numbers is ued for each run. A new instance of the algorithm is created allowing a maximum of 5000 iterations, and an instance of the ExponentFunction problem is created in 5 dimensions. The search is executed by passing the problem to the algorithm, and the best result achieved by the algorithm is printed once the search is completed.

srand(1) # set the random number seed to 1
algorithm = AdaptiveRandomSearchAlgorithm.new(5000) # limit to 5000 iterations 
problem = ExponentFunction.new(5) # create a problem with 5 dimensions
best = algorithm.search(problem) # execute the search
puts "Best Solution: #{best}" # display the best solution

Once the general approach is understood, one may begin playing with the parameters defined in the algorithms constructor that govern the dynamic adaptation of the step size. For example, smaller or larger factors may be used, and the frequency of large or small step size increase and decreases may be varied. The acceptance criteria for when the current working solution is replaced may be made probabilistic rather than deterministically greedy. Finally, one may experiment with hyper-spherical sampling neighbourhoods in the next_step function and with sampling distributions other than uniform such as Gaussian with the origin as the mean and the step_size as the standard deviation.

Source Code
Download: AdaptiveRandomSearch.rb

Further Reading

Papers:

M. Schumer and Steiglitz, Adaptive step size random search. IEEE Transactions on Automatic Control 13(3), pages 270- 276. 1968
S. F. Masri and G. Bekey, Global Optimization Algorithm Using Adaptive Random Search. Applied Mathematics and Computation 7, pages 353-375. 1980
Rajeeva Kumar, Pierre T. Kabamba, David C. Hyland, Analysis and parameter selection for an adaptive random search algorithm. Mathematics and Computers in Simulation 68(2), pages 95-103. 2005
Andrei A. Prudius, Adaptive Random Search Methods for Simulation Optimization. Ph.D. Thesis, Georgia Institute of Technology. 2007

Sites:

Random optimization. Wikipedia
Stochastic optimization. Wikipedia

Iterated Local Search: Cycling global and local search

2009-01-13T11:08:00.022+11:00

A Hill Climbing algorithm such as the one demonstrated in Localized Random Search is an example of a Local Search Algorithm. A local search involves moving through a search space from one neighboring candidate solution to the next, such as from vertex to vertex along the edges of a graph. The Iterated Local Search algorithm involves the repeated application of a local search algorithm applied to the candidate solutions found by a broader search process that involves a biased random walk through the search space.

The algorithm works by first selecting a starting point for the search either randomly or via a domain specific construction heuristic. The starting position is then optimized to an approximation of the local optimum using a local search strategy. The algorithm loop involves three steps: a perturbation of the current position, the application of the local search to the perturbation, and an acceptance decision of whether or not the locally optimizing candidate solution should replace the current working position for the search.

A good starting point is important for shorter algorithm runs, whereas this importance degrades as the length of the run is increased because the algorithm is given more opportunity to improve. The Iterative Local Search algorithm resembles random restart of a local search technique, although out-performs this technique given the dependence between each new starting position of the search. The perturbation of the current working position provides long jumps in the search space that are further refined by the local search. If the jumps are too large, performance degrades as the behaviors resembles random restart. If the jumps are too small in the search space, the broader search technique resembles the local search algorithm.

The goal is to find a balance between the perturbation and local search mechanisms within the algorithm such that the local search technique cannot easily undo the longer jumps made by the perturbations. There is a relationship between the perturbations and the acceptance criteria for replacing the current working position. Always accepting improved points may result in a greedy algorithm that only accepts small steps and cannot explore distant regions of the search space. Longer jumps made by the perturbation mechanism may require a more flexible (less greedy) acceptance criteria for the search to progress.

Iterated Local Search of the Traveling Salesman Problem
The Iterated Local Search algorithm is suited to problem domains where the good solutions in the search space are clustered. The approach is more commonly applied to combinatorial optimization problems like the Traveling Salesman Problem (TSP) where this clustering of good solutions is manifest as the sharing of common edges. This guide demonstrates the application of the Iterated Local Search algorithm applied to a standardized instance of the TSP from the TSP library (TSPLIB).

The Berlin52TSP class contains a the city coordinate data and optimal tour data as COORDINATES and OPTIMAL_TOUR constants respectively, taken from the TSPLIB definition files for the Berlin52 problem instance. The initialize() constructor sets the accessible @num_cities instance variable and calculates the @optimal_tour_length used by the is_optimal?(scoring) method for assessing scores to see if an optimal tour has been discovered. The evaluate(permutation)method calculates the distance of a given tour permutation which is an array of non-repeating integers that specify cities between 1 and 52 inclusively. Evaluation involves calculating the Euclidean distance for each of the city pairs using the problem definition coordinates. The euc_2d(c1, c2) method matches the Euclidean calculation specified by the TSPLIB documentation, most notably involving the rounding of calculated distances to integers for this symmetrical TSP instance.

class Berlin52TSP
  OPTIMAL_TOUR = [1,49,32,45,19,41,8,9,10,43,33,51,11,52,14,13,47,26,
    27,28,12,25,4,6,15,5,24,48,38,37,40,39,36,35,34,44,46,16,29,50,20,
    23,30,2,7,42,21,17,3,18,31,22]
        
  COORDINATES = [[565, 575],[25, 185],[345, 750],[945, 685],[845, 655],
  [880, 660],[25, 230],[525, 1000],[580, 1175],[650, 1130],[1605, 620], 
  [1220, 580],[1465, 200],[1530, 5],[845, 680],[725, 370],[145, 665],
  [415, 635],[510, 875], [560, 365],[300, 465],[520, 585],[480, 415],
  [835, 625],[975, 580],[1215, 245],[1320, 315],[1250, 400],[660, 180],
  [410, 250],[420, 555],[575, 665],[1150, 1160],[700, 580],[685, 595],
  [685, 610],[770, 610],[795, 645],[720, 635],[760, 650],[475, 960],
  [95, 260],[875, 920],[700, 500],[555, 815],[830, 485],[1170, 65],
  [830, 610],[605, 625],[595, 360],[1340, 725],[1740, 245]]
  
  attr_reader :num_cities

  def initialize()
    @num_cities = COORDINATES.length        
    @optimal_tour_length = evaluate(OPTIMAL_TOUR) # calculate
  end
  
  def evaluate(permutation)
    dist = 0    
    permutation.each_with_index do |c1, i|
      c2 = (i==@num_cities-1) ? permutation[0] : permutation[i+1] 
      dist += euc_2d(COORDINATES[c1-1], COORDINATES[c2-1])
    end
    return dist
  end
  
  def euc_2d(c1, c2)
    # As defined in TSPLIB'95 (EUC_2D)
    Math::sqrt((c1[0] - c2[0])**2 + (c1[1] - c2[1])**2).round
  end

  def is_optimal?(scoring)
    scoring == optimal_score
  end

  def optimal_score
    @optimal_tour_length
  end
  
  # true if s1 is better score than s2
  def is_better?(s1, s2)
    s1 < s2 # minimizing
  end
end

The Solution class is a classical container for candidate solutions, holding a TSP permutation in the immutable @data instance variable and managing a mutable permutation scoring in the @score instance variable.

class Solution
  attr_reader :data
  attr_accessor :score
  
  def initialize(data)
    @data = data
    @score = 0.0/0.0 # NaN
  end
  
  def to_s
    "[#{@data.inspect}] (#{@score})"
  end    
end

The IteratedLocalAlgorithm class is quite large. The search(problem) method is generalized and executes the iterated local search algorithm on the provided problem definition. It starts by calling generate_initial_solution(problem) to prepare a starting point for the search and refining it with a call to local_search_solution(current, problem). The algorithms main loop is a text book implementation involving repeated rounds of generating a perturbed versions of the current working solution, refining it with a local search procedure, and in this case using a greedy (candidate must be better) acceptance criteria.

The meat of the algorithm are specific to TSP permutations. The generate_initial_solution(problem) generates a random permutation as the starting point of the search. The perturb_solution(solution) method generates a 4-opt variation of a given permutation, also known as a double-bridge move. Basically the permutation is split into 4 pieces and reordered. This function may be called repeatedly if more drastic perturbations are desired, perhaps on large problem instances. The local_search_solution(solution, problem) method is an iterative application of the 2-opt procedure, that evaluates generated solutions and greedily maintains the best candidate solution discovered. The 2-opt procedure involves breaking the tour permutation into two parts (two break points) and reconnecting the tour back together by switching the end points (reversing the sequence between the break points). This is a classical and fast local search procedure for TSP and here it is repeated 30 times.

class IteratedLocalAlgorithm
  attr_accessor :max_iterations
  attr_reader :best_solution
  
  def initialize(max_iterations)
    @max_iterations = max_iterations
  end
  
  # execute a iterated local search on the provided problem
  def search(problem)
    # random starting point
    current = generate_initial_solution(problem)
    @best_solution = current
    evaluate_candidate_solution(current, problem)
    # local search
    local_best = local_search_solution(current, problem)
    current = local_best if problem.is_better?(local_best.score, current.score)
    curr_it = 0
    begin
      # generate perturbation
      pert_solution = perturb_solution(current)
      evaluate_candidate_solution(pert_solution, problem)
      # local search
      local_best = local_search_solution(pert_solution, problem)
      # greedy acceptance
      current = local_best if problem.is_better?(local_best.score, current.score)
      curr_it += 1
    end until should_stop?(curr_it, problem)
    return @best_solution
  end
  
  def should_stop?(curr_it, problem)
    (curr_it >= @max_iterations) or problem.is_optimal?(best_solution.score)
  end
  
  def generate_initial_solution(problem)
    all = Array.new(problem.num_cities) {|i| (i+1)}
    permutation = Array.new(all.length) {|i| all.delete_at(rand(all.length))}
    return Solution.new(permutation)
  end

  def perturb_solution(solution)    
    data =  solution.data
    length = data.length
    # double-bridge move (4-opt), break into 4 parts (a,b,c,d)
    pos1 = 1 + rand(length / 4)
    pos2 = pos1 + 1 + rand(length / 4)
    pos3 = pos2 + 1 + rand(length / 4)
    # put it back together (a,d,c,b)
    perm = data[0...pos1] + data[pos3...length] + 
      data[pos2...pos3] + data[pos1...pos2]
    return Solution.new(perm)
  end

  def local_search_solution(solution, problem)
    # greedy iterated 2-opt
    30.times do
      candidate = two_opt_solution(solution)
      evaluate_candidate_solution(candidate, problem)
      if problem.is_better?(candidate.score, solution.score)
        solution = candidate 
      end
    end
    return solution
  end

  def two_opt_solution(solution)
    perm = Array.new(solution.data) # copy
    # select a sub-sequence
    c1, c2 = rand(perm.length), rand(perm.length)
    c2 = rand(perm.length) while c1 == c2
    # ensure c1 is low and c2 is high
    c1, c2 = c2, c1 if c2 < c1
    # reverse sub-sequence
    perm[c1...c2] = perm[c1...c2].reverse
    return Solution.new(perm)
  end

  def evaluate_candidate_solution(solution, problem)
    solution.score = problem.evaluate(solution.data)
    # keep track of the best solution found
    if problem.is_better?(solution.score, @best_solution.score)
      @best_solution = solution
      puts " > new best: #{solution.score}"               
    end
  end
end

The demonstration of the algorithm involves first seeding the global random number generator to one, and constructing both an instance of the problem and the algorithm. The algorithm is executed on the problem instance and the algorithm stop condition is triggered if the optimal solution is discovered or a maximum of 2000 iterations of the main algorithm loop are completed.

srand(1) # set the random number seed to 1
algorithm = IteratedLocalAlgorithm.new(1000) # limit to 1000 iterations 
problem = Berlin52TSP.new # create a problem
best = algorithm.search(problem) # execute the search
puts "Best Solution: #{best}" # display the best solution

There is plenty of room for extension of this interesting guide. The problem definition may be made generic such that it loads problem instance descriptions and optimal tours from TSPLIB files. The problem class may be made computationally more efficient by pre-calculating the city distance matrix on construction and looking up these distance evaluations in the evaluate method. The use of a randomly generated starting point for the search suggests that the algorithm may need to be executed for a long time to approximate the global optima (longer than 2000 iterations). The algorithm may be adjusted to use a deterministically generated nearest-neighbor solution as the starting point of the search that expected to be a much better quality solution. Finally, more advanced local search procedures exist for the TSP (such as Lin-Kernighan), which may be integrated into the algorithm.

Source Code
Download: IteratedLocalSearch.rb

Further Reading

Papers

Helena Ramalhinho-Lourenço, Olivier C. Martin, Thomas Stützle, Iterated Local Search. (pre-print) Handbook of Metaheuristics, pages 320-353. Springer. 2003
Thomas Stützle, Iterated local search for the quadratic assignment problem. (pre-print) European Journal of Operational Research. 174(3), pages 1519-1539. 2006
David S. Johnson, Local optimization and the Traveling Salesman Problem. Automata, Languages and Programming, pages 446-461. 1990

Sites

Have you found a bug? Know of a great reference? Got an idea for a guide? Please send an email to jasonb@InspiredAlgorithms.com

Need Help? Want to discuss this post? Please visit the Inspired Algorithms User Group.

Localized Random Search: Sampling around known solutions

2009-01-11T10:00:00.005+11:00

This guide extends Random Search by generating new candidate solutions as variations of the the best known solution. This approach is generally referred to as Localized Random Search or as Mutation Hill Climbing by the evolutionary computation community. This approach works by firstly selecting a starting point for the search by generating a random solution or selecting the best solution from a set of randomly generated solutions. The algorithm then proceeds to generate new candidate solutions as variations of the current best, replacing the current best if and only if the new mutant candidate solution is better or has the same objective scoring.

The search is localized to a known 'good' area of the problem space, and can rapidly locate the locally optimal solution. Naturally, the starting point for the localized search is critically important as once a starting point for the search is selected, the algorithm can at best return the locally optimal solution. As such, the approach is commonly used to refine the results of another more broadly focused search algorithm. Additionally, Localized Random Search can operate upon a set of candidate starting points in parallel or be restarted after it has converged to a local optimum, both strategies of which may increase the chances better approximating the global optimal solution.

Mutation Hill Climbing Algorithm on the OneMax Problem
The evolutionary computation community have produced some interesting work on Localized Random Search. Specifically regarding the binary string optimization under difficult and deceptive cost functions, referring to the algorithm as a stochastic hill climber or mutation hill climbing algorithm. This guide demonstrates a mutation hill climbing algorithm on a binary string optimization problem called OneMax.

The OneMax problem involves the optimization of a string of binary digits (1's and 0's) of an arbitrary fixed length where the cost assigned to a given string reflects the number of 1's in the string. The optimal scoring in this maximization problem is equal to the length of the string, sometimes referred to as unity. It is an interesting black box function because although gradient information is provided (more or less 1's), the algorithm is not informed which positions in the string are right or wrong.

The OneMaxFunction class encapsulates the problem definition. It expects a string of characters comprised of 1's and 0's which are assessed via the evaluate(bitstring) method. The complexity of the problem is managed by the @length instance variable, which defaults to 16 bits.

class OneMaxFunction
  attr_reader :length

  def initialize(length=16)
    @length = length
  end

  def evaluate(bitstring)
    bitstring.inject(0) {|sum, x| sum + ((x=='1') ? 1 : 0)}
  end  

  def is_optimal?(scoring)
    scoring == optimal_score
  end

  def optimal_score
    @length
  end
  
  # true if s1 has the same or better score than s2
  def is_same_or_better?(s1, s2)
    s1 >= s2 # maximizing
  end
end

The Solution class is a simple convenience class, providing an immutable (problem non-specific) representation of candidate solutions with the data instance variable, and a mutable score instance variable for assigned costing. The to_s method enumerates data to ensure that the binary strings used in this example print nicely (for example: 101100101).

class Solution
  attr_reader :data
  attr_accessor :score
  
  def initialize(data)
    @data = data
    @score = 0.0/0.0 # NaN
  end
  
  def to_s
    "[#{@data.collect{|x|x}}] (#{@score})"
  end    
end

The LocalizedRandomSearchAlgorithm class is initialized with a maximum number of iterations to complete. The important method of this class is search(problem) which executes a search on the provided problem definition. The method calls generate_random_solution(problem) to create a random candidate solution as a starting point for the search. The algorithm then loops until the stop condition is triggered should_stop?(curr_it, problem) creating mutated variations of the current best known solution (@best_solution instance variable), and replacing the current best known solution of the mutated candidate has the same or equal scoring. This acceptance criteria (of not just accepting improved scores) is important as it allows the search to continue across flat regions of the search space, potentially escaping local optima.

The generate_mutate_solution(solution) method creates a new candidate solution from a provided solution. It works by enumerating the provided solutions data one bit (character) at a time, probabilistically mutating the string by inverting a give bit position ('0'=>'1' or '1'=>'0'). This is achieved by calling the should_mutate?(length) method for each bit in the string which permits a bit-flip mutation with the probability of 1/L, where L is the length of the bitstring being optimized. This has be demonstrated to be an optimal mutation rate for the OneMax function and related linear binary optimization problems (see references), can provide a good heuristic for problems with unknown response surfaces.

class LocalizedRandomSearchAlgorithm
  attr_accessor :max_iterations
  attr_reader :best_solution
  
  def initialize(max_iterations)
    @max_iterations = max_iterations
  end
  
  # execute a random search on the provided problem
  def search(problem)
    @best_solution = generate_random_solution(problem) # starting point
    @best_solution.score = problem.evaluate(@best_solution.data)
    curr_it = 0
    begin
      # generate mutant
      candidate = generate_mutate_solution(@best_solution)
      candidate.score = problem.evaluate(candidate.data)
      # compare to current best
      if problem.is_same_or_better?(candidate.score, @best_solution.score)
        @best_solution = candidate
        puts " > new best: #{@best_solution.score}"    
      end
      curr_it += 1
    end until should_stop?(curr_it, problem)
    return @best_solution
  end
  
  def should_stop?(curr_it, problem)
    (curr_it >= @max_iterations) or problem.is_optimal?(best_solution.score)
  end
  
  def generate_random_solution(problem)
    bitstring = Array.new(problem.length) {|i| (rand<0.5) ? "1" : "0"}
    return Solution.new(bitstring)
  end

  def generate_mutate_solution(solution)
    data = solution.data
    bitstring = Array.new(data.length) do |i| 
      if should_mutate?(data.length)
         (data[i]=='0') ? "1" : "0" # invert
      else
        data[i]
      end
    end
    return Solution.new(bitstring)
  end
  
  def should_mutate?(length)
    rand < 1.0/length
  end
end

The algorithm is demonstrated by first seeding the global random number generator to 1, creating an instance of the problem initialized with 32 bits, and an instance of the algorithm with a maximum of 1000 algorithm iterations. The algorithm is executed by calling and passed the problem instance. The algorithms current progress is displayed as improved or equally best scored solutions are found and the best solution is displayed at the end of the run.

srand(1) # set the random number seed to 1
algorithm = LocalizedRandomSearchAlgorithm.new(1000) # limit to 1000 iterations
problem = OneMaxFunction.new(32) # create a problem with 32 bits
best = algorithm.search(problem) # execute the search
puts "Best Solution: #{best}" # display the best solution

This implementation can be extended upon by modifying the LocalizedRandomSearchAlgorithm class to manage a population of 'best' solutions in parallel, or to generate a small set of mutate candidate solutions from the current best solution each iteration. The algorithm implementation is generalized and can be applied to other problems that align to the OneMaxFunction interface by specializing the generate_random_solution(problem) and generate_mutate_solution(solution) methods.

Source Code
Download: LocalizedRandomSearch.rb

Further Reading

Papers

Thomas Bäck, "Optimal Mutation Rates in Genetic Search". Proceedings of the Fifth International Conference on Genetic Algorithms, pages 2-9, 1993.
Heinz Mühlenbein, "How Genetic Algorithms Really Work: I. Mutation and Hillclimbing". Parallel Problem Solving from Nature 2, PPSN-II, Belgium, pp. 15-26, 1992.
James C. Spall, Stochastic Optimization. In James E. Gentle, Wolfgang Härdle, Yuichi Mori, Handbook of Computational Statistics. 2004

Sites

Hill Climbing. Wikipedia

To get help or discuss this post, please visit the Inspired Algorithms User Group.

Random Search: A baseline technique

2009-01-10T10:00:00.009+11:00

The classical baseline for stochastic optimization is the random search algorithm, also referred to as blind search and stochastic generate and test. The strategy of this algorithm involves randomly generating viable solutions from the problems search space and evaluating them. It is used as a baseline for comparison because it is general enough to be applied to any problem for which candidate solutions can be generated and tested.

Each new randomly generated candidate solution is independent of the generated candidate solutions that came before it. Importantly, the algorithm keeps track of the candidate solution found so far using a simple ratcheting function (current best is only replaced by something as good as or better). Given that the random search algorithm samples the problem search space in an unbiased (uniform) manner it may discover areas of interest not considered by conventional wisdom.

The penalty of the algorithms simplicity and lack of memory is its inefficiency. The algorithm does not make use of any problem specific information such as gradient of the response surface. It may perform worse than a complete enumeration of the search space given that there are no mechanisms to bias or prevent the procedure from re-investigating already visited regions of the search space.

Random Search for Continuous Function Optimization
This guide demonstrates the random search algorithm applied to a simple summed squaring function. We start by defining the problem as a class with methods for assessing and comparing candidate solutions. The function is defined as the sum of the squared parameters, where the number of function parameters defines the difficulty or number of dimensions of the problem. Each function parameter is constrained between -5.12 and +5.12.

Candidate solutions are assessed via the evaluate method and are expected as arrays of floating point values. The global solution to the problem is known at 0.0 in each dimension which evaluates to the optimal score of 0.0. In one-dimension the function's response surface can be graphed which looks like the letter 'u' with the optima in the middle of the function. In two-dimensions the response surface looks like a basin. The known optimal scoring is made available in the is_optimal?(scoring) function, useful for an algorithm for accessing whether the optimal solution has been achieved. The problem also defines a solution score comparison, enforcing the minimization constrain of the problem.

class SquaringFunction
  attr_reader :dimensions, :min, :max

  def initialize(dimensions=2)
    @dimensions = dimensions
    @min, @max = -5.12, +5.12
  end

  def evaluate(vector)
    vector.inject(0) {|sum, x| sum +  (x ** 2.0)}
  end  

  def is_optimal?(scoring)
    scoring == optimal_score
  end

  def optimal_score
    0.0
  end
  
  # true if s1 has a better score than s2
  def is_better?(s1, s2)
    s1 < s2 # minimizing
  end
end

In order to address this problem using the random search algorithm, a candidate solution class and algorithm class are defined. The Solution class is not problem specific, proving a immutable container for solution data in an instance variable which must be set on construction, as well as a score which can be assigned at any time. A default score of NaN (not a number) is assigned defensively to ensure that an error occurs if the solution is not evaluated and is compared to another solution that is evaluated.

class Solution
  attr_reader :data
  attr_accessor :score
  
  def initialize(data)
    @data = data
    @score = 0.0/0.0 # NaN
  end
  
  def to_s
    "[#{@data.inspect}] (#{@score})"
  end    
end

The RandomSearchAlgorithm class requires a maximum number of iterations as parameters on construction which is stored as an instance variable. The core method of the class is search(problem) that executes a random search on the provided problem definition. A single iteration of the algorithm involves the generation of a random candidate solution and its evaluation. The algorithm loops until the stop condition should_stop?(curr_it, problem) is triggered which occurs if the maximum number of iterations have elapsed or if the optimal solution is discovered. The generation of random solutions in the generate_random_solution(problem) method is the only problem-specific logic in the algorithm, creating a random floating point vector within there bounds of the problem space. Finally, the evaluate_candidate_solution(solution, problem) method evaluates a given candidate solution against the problem definition and keeps track of the best solution found.

class RandomSearchAlgorithm
  attr_accessor :max_iterations
  attr_reader :best_solution
  
  def initialize(max_iterations)
    @max_iterations = max_iterations
  end
  
  # execute a random search on the provided problem
  def search(problem)    
    @best_solution = nil
    curr_it = 0
    begin
      candidate = generate_random_solution(problem)
      evaluate_candidate_solution(candidate, problem)
      curr_it += 1
    end until should_stop?(curr_it, problem)
    return @best_solution
  end
  
  def should_stop?(curr_it, problem)
    (curr_it >= @max_iterations) or problem.is_optimal?(best_solution.score)
  end
  
  def generate_random_solution(problem)
    real_vector = Array.new(problem.dimensions) do
      next_bfloat(problem.min, problem.max)
    end
    return Solution.new(real_vector)
  end

  def next_bfloat(min, max)
    min + ((max - min) * rand)
  end
  
  def evaluate_candidate_solution(solution, problem)
    solution.score = problem.evaluate(solution.data)
    # keep track of the best solution found
    if @best_solution.nil? or
      problem.is_better?(solution.score, @best_solution.score)
      @best_solution = solution
      puts " > new best: #{solution.score}"               
    end
  end  
end

Now we can execute the algorithm on the problem. This involves initializing the global random number generator to 1 and initializing the algorithm with the maximum number of 1000 iterations to execute. The problem instance is created with 5 decision variables (dimensions) and the algorithm is executed on the instance.

srand(1) # set the random number seed to 1
algorithm = RandomSearchAlgorithm.new(1000) # limit to 1000 iterations
problem = SquaringFunction.new(5) # create a problem with 5 dimensions
best = algorithm.search(problem) # execute the search
puts "Best Solution: #{best}" # display the best solution

This example provides a demonstration of a standard separation of problem, algorithm, and solution, as well as modularized responsibility within each class. For example, the RandomSearchAlgorithm class can be directly applied to any function optimization problem, or with a replacement of the generate_random_solution(problem) method it can be applied to any problem definition that aligns to the SquaringFunction classes simple interface.

Source Code
Download: RandomSearch.rb

Further Reading
Given the simplicity of the approach there are few resources dedicated to describing it. Instead the following provide some general resources for the broader field of stochastic optimization.

Books

James C. Spall, Introduction to Stochastic Search and Optimization. Wiley-IEEE, 2005
Anatoly Zhigljavsky and Antanas Žilinskas, Stochastic Global Optimization (Springer Optimization and Its Applications). Springer, 2008
James C. Spall, Stochastic Optimization. In James E. Gentle, Wolfgang Härdle, Yuichi Mori, Handbook of Computational Statistics. 2004

Sites

Stochastic Optimization. Wikipedia

To get help with or discuss this post, please visit the Inspired Algorithms User Group.

Welcome to Inspired Algorithms

2009-01-10T08:50:00.002+11:00

Hi, my name is Jason Brownlee and I am really excited to present Inspired Algorithms. The goal of this project is to write and offer the largest set of explained and executable computational intelligence algorithms on the planet.

Toward this goal I (at this stage, alone) will endeavor to write and publish a significant number of independent guides that provide a functional introduction and working program code for a wide array of computational intelligence algorithms. The best of these guides will be collected and re-released as a compendium in the form of (1) a free e-book, and (2) a self-published printed book offered at a low cost.

The phrase "inspired algorithms" is used in this project as an umbrella term to refer to machine learning algorithms drawn from across many fields of artificial intelligence not limited to computational intelligence, biologically inspired computing, metaheuristics, natural computation, and biomemetics.

An intention of this project is openness. The guides will be readable for a general audience with some programming, computer science, or related background. The programs will be understandable for the average programmer. A diverse set of popular and interesting algorithms will be covered, and a useful set of books, papers, and websites will be listed with each guide for further reading and extension.

An Inspired Algorithms Open Source Project has been created to maintain all source code presented in the guides. Source code is presented in the Ruby Programming Language and is released under the GNU General Public License v3.

An Inspired Algorithms User Group has been created where we can discuss algorithms, guides and help each other out with what can sometimes be a complex but always rewarding field. There are no constraints on who can participate, and I encourage interested amateurs, professional developers, and career academics alike to ask and answer questions. I hope together we can refine this set of guides and algorithm implementation to make them both accurate, complete, and easy to understand.

You may subscribe to the guides on this site and read them as they are released by pointing your news reader to the Inspired Algorithms Feed.

If you're interested in the author, you can find out more about me on my LinkedIn profile as well as my personal blog.

I look forward to this adventure of learning and discovery, and hope something here sparks your interest or curiosity in what is a truly exciting and powerful frontier of computation.