lib/multidim/powell.rb
# = powell.rb -
# Minimization- Minimization algorithms on pure Ruby
# Copyright (C) 2010 Claudio Bustos
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
#
# This algorith was adopted and ported into Ruby from Apache-commons
# Math library's PowellOptimizer.java and
# BaseAbstractMultivariateVectorOptimizer.java
# files. Therefore this file is under Apache License Version 2.
#
# Powell's Algorithm for Multidimensional minimization
require "#{File.expand_path(File.dirname(__FILE__))}/point_value_pair.rb"
require "#{File.expand_path(File.dirname(__FILE__))}/../minimization.rb"
module Minimization
class ConjugateDirectionMinimizer
attr_accessor :max_iterations
attr_accessor :max_brent_iterations
attr_accessor :x_minimum
attr_accessor :f_minimum
# default maximum Powell's iteration value
Max_Iterations_Default = 100
# default Brent iteration value
MAX_BRENT_ITERATION_DEFAULT = 10 # give a suitable value
def initialize(f, initial_guess, lower_bound, upper_bound)
@iterations = 0
@max_iterations = Max_Iterations_Default
@evaluations = 0
@max_brent_iterations = MAX_BRENT_ITERATION_DEFAULT
@converging = true
# set minimizing function
@f = f
@start = initial_guess
@lower_bound = lower_bound
@upper_bound = upper_bound
# set maximum and minimum coordinate value a point can have
# while minimization process
@min_coordinate_val = lower_bound.min
@max_coordinate_val = upper_bound.max
# validate input parameters
check_parameters
end
# return the convergence of the search
def converging?
return @converging
end
# set minimization function
def f(x)
@f.call(x)
end
# validate input parameters
def check_parameters
if (!@start.nil?)
dim = @start.length
if (!@lower_bound.nil?)
# check for dimension mismatches
raise "dimension mismatching #{@lower_bound.length} and #{dim}" if @lower_bound.length != dim
# check whether start point exeeds the lower bound
0.upto(dim - 1) do |i|
v = @start[i]
lo = @lower_bound[i]
raise "start point is lower than lower bound" if v < lo
end
end
if (!@upper_bound.nil?)
# check for dimension mismatches
raise "dimension mismatching #{@upper_bound.length} and #{dim}" if @upper_bound.length != dim
# check whether strating point exceeds the upper bound
0.upto(dim - 1) do |i|
v = @start[i]
hi = @upper_bound[i]
raise "start point is higher than the upper bound" if v > hi
end
end
if (@lower_bound.nil?)
@lower_bound = Array.new(dim)
0.upto(dim - 1) do |i|
@lower_bound[i] = Float::INFINITY # eventually this will occur an error
end
end
if (@upper_bound.nil?)
@upper_bound = Array.new(dim)
0.upto(dim - 1) do |i|
@upper_bound[i] = -Float::INFINITY # eventually this will occur an error
end
end
end
end
# line minimization using Brent's minimization
# == Parameters:
# * <tt>point</tt>: Starting point
# * <tt>direction</tt>: Search direction
#
def brent_search(point, direction)
n = point.length
# Create a proc to minimize using brent search
# Function value varies with alpha value and represent a point
# of the minimizing function which is on the given plane
func = proc{ |alpha|
x = Array.new(n)
0.upto(n - 1) do |i|
# create a point according to the given alpha value
x[i] = point[i] + alpha * direction[i]
end
# return the function value of the obtained point
f(x)
}
# create Brent minimizer
line_minimizer = Minimization::Brent.new(@min_coordinate_val, @max_coordinate_val, func)
# iterate Brent minimizer for given number of iteration value
0.upto(@max_brent_iterations) do
line_minimizer.iterate
end
# return the minimum point
return {:alpha_min => line_minimizer.x_minimum, :f_val => line_minimizer.f_minimum}
end
end
# = Powell's Minimizer.
# A multidimensional minimization methods
# == Usage.
# require 'minimization'
# f = proc{ |x| (x[0] - 1)**2 + (2*x[1] - 5)**2 + (x[2]-3.3)**2}
# min = Minimization::Powell.minimize(f, [1, 2, 3], [0, 0, 0], [5, 5, 5])
# min.f_minimum
# min.x_minimum
#
class Powell < ConjugateDirectionMinimizer
attr_accessor :relative_threshold
attr_accessor :absolute_threshold
# default of relative threshold
RELATIVE_THRESHOLD_DEFAULT = 0.1
# default of absolute threshold
ABSOLUTE_THRESHOLD_DEFAULT =0.1
# == Parameters:
# * <tt>f</tt>: Minimization function
# * <tt>initial_guess</tt>: Initial position of Minimization
# * <tt>lower_bound</tt>: Lower bound of the minimization
# * <tt>upper_bound</tt>: Upper bound of the minimization
#
def initialize(f, initial_guess, lower_bound, upper_bound)
super(f, initial_guess.clone, lower_bound, upper_bound)
@relative_threshold = RELATIVE_THRESHOLD_DEFAULT
@absolute_threshold = ABSOLUTE_THRESHOLD_DEFAULT
end
# Obtain new point and direction from the previous point,
# previous direction and a parameter value
# == Parameters:
# * <tt>point</tt>: Previous point
# * <tt>direction</tt>: Previous direction
# * <tt>minimum</tt>: parameter value
#
def new_point_and_direction(point, direction, minimum)
n = point.length
new_point = Array.new(n)
new_dir = Array.new(n)
0.upto(n - 1) do |i|
new_dir[i] = direction[i] * minimum
new_point[i] = point[i] + new_dir[i]
end
return {:point => new_point, :dir => new_dir}
end
# Iterate Powell's minimizer one step
# == Parameters:
# * <tt>f</tt>: Function to minimize
# * <tt>starting_point</tt>: starting point
# * <tt>lower_bound</tt>: Lowest possible values of each direction
# * <tt>upper_bound</tt>: Highest possible values of each direction
# == Usage:
# minimizer = Minimization::Powell.new(proc{|x| (x[0] - 1)**2 + (x[1] -1)**2},
# [0, 0, 0], [-5, -5, -5], [5, 5, 5])
# while minimizer.converging?
# minimizer.iterate
# end
# minimizer.x_minimum
# minimizer.f_minimum
#
def iterate
@iterations += 1
# set initial configurations
if(@iterations <= 1)
guess = @start
@n = guess.length
# initialize all to 0
@direc = Array.new(@n) { Array.new(@n) {0} }
0.upto(@n - 1) do |i|
# set diagonal values to 1
@direc[i][i] = 1
end
@x = guess
@f_val = f(@x)
@x1 = @x.clone
end
fx = @f_val
fx2 = 0
delta = 0
big_ind = 0
alpha_min = 0
0.upto(@n - 1) do |i|
direction = @direc[i].clone
fx2 = @f_val
# Find line minimum
minimum = brent_search(@x, direction)
@f_val = minimum[:f_val]
alpha_min = minimum[:alpha_min]
# Obtain new point and direction
new_pnd = new_point_and_direction(@x, direction, alpha_min)
new_point = new_pnd[:point]
new_dir = new_pnd[:dir]
@x = new_point
if ((fx2 - @f_val) > delta)
delta = fx2 - @f_val
big_ind = i
end
end
# convergence check
@converging = !(2 * (fx - @f_val) <= (@relative_threshold * (fx.abs + @f_val.abs) + @absolute_threshold))
# storing results
if((@f_val < fx))
@x_minimum = @x
@f_minimum = @f_val
else
@x_minimum = @x1
@f_minimum = fx
end
direction = Array.new(@n)
x2 = Array.new(@n)
0.upto(@n -1) do |i|
direction[i] = @x[i] - @x1[i]
x2[i] = 2 * @x[i] - @x1[i]
end
@x1 = @x.clone
fx2 = f(x2)
if (fx > fx2)
t = 2 * (fx + fx2 - 2 * @f_val)
temp = fx - @f_val - delta
t *= temp * temp
temp = fx - fx2
t -= delta * temp * temp
if (t < 0.0)
minimum = brent_search(@x, direction)
@f_val = minimum[:f_val]
alpha_min = minimum[:alpha_min]
# Obtain new point and direction
new_pnd = new_point_and_direction(@x, direction, alpha_min)
new_point = new_pnd[:point]
new_dir = new_pnd[:dir]
@x = new_point
last_ind = @n - 1
@direc[big_ind] = @direc[last_ind]
@direc[last_ind] = new_dir
end
end
end
# Convenience method to minimize
# == Parameters:
# * <tt>f</tt>: Function to minimize
# * <tt>starting_point</tt>: starting point
# * <tt>lower_bound</tt>: Lowest possible values of each direction
# * <tt>upper_bound</tt>: Highest possible values of each direction
# == Usage:
# minimizer = Minimization::Powell.minimize(proc{|x| (x[0] - 1)**2 + (x[1] -1)**2},
# [0, 0, 0], [-5, -5, -5], [5, 5, 5])
# minimizer.x_minimum
# minimizer.f_minimum
#
def self.minimize(f, starting_point, lower_bound, upper_bound)
min = Minimization::Powell.new(f, starting_point, lower_bound, upper_bound)
while min.converging?
min.iterate
end
return min
end
end
end