Class: MixedModels::NelderMead

Inherits:
Object
  • Object
show all
Defined in:
lib/mixed_models/NelderMeadWithConstraints.rb

Overview

Nelder Mead Minimizer with Bound Contraints. A multidimensional minimization methods with the possibility to impose constraints lower_bound <= x <= upper_bound for all i.

Usage

min=MixedModels::NelderMead.new(start_point: [1,2]) {|x| (x[0] - 2)**2 + (x[1] - 5)**2}
while min.converging?
  min.iterate
end
min.x_minimum
min.f_minimum

Instance Attribute Summary collapse

Class Method Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(start_point:, lower_bound: nil, upper_bound: nil, epsilon: 1e-6, max_iterations: 1e6, &f) ⇒ NelderMead

Arguments

  • start_point - an Array specifying the initial point for the minimization

  • lower_bound - an Array of lower bounds for each coordinate of the optimal solution

  • upper_bound - an Array of upper bounds for each coordinate of the optimal solution

  • epsilon - a small number specifying the thresholds for the convergence check: absolute_threshold = epsilon and relative_threshold = 100 * epsilon

  • max_iterations - the maximum number of iterations

  • f - the objective function as a Proc object



75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
 
# File 'lib/mixed_models/NelderMeadWithConstraints.rb', line 75

def initialize(start_point:, lower_bound: nil, upper_bound: nil, 
               epsilon: 1e-6, max_iterations: 1e6, &f)
  @start_point = start_point
  @lower_bound = lower_bound
  @upper_bound = upper_bound

  @rho   = 1.0 # Reflection coefficient
  @khi   = 2.0 # Expansion coefficient
  @gamma = 0.5 # Contraction coefficient
  @sigma = 0.5 # Shrinkage coefficient

  @epsilon             = epsilon
  @max_iterations      = max_iterations
  @relative_threshold  = 100 * @epsilon
  @absolute_threshold  = @epsilon
  @x_minimum           = nil
  @f_minimum           = nil
  @f                   = f

  n = start_point.length
  # create and initialize start configurations
  if @start_configuration == nil
    # sets the start configuration point as unit
    self.start_configuration = Array.new(n) { 1.0 }
  end

  if lower_bound.nil? then
    @lower_bound = Array.new(n) { -Float::INFINITY }
  else
    raise "Lower bound should be of the same length as the start point" unless lower_bound.length == n
    @lower_bound = lower_bound 
  end
  if upper_bound.nil? then
    @upper_bound = Array.new(n) { Float::INFINITY }
  else
    raise "Upper bound should be of the same length as the start point" unless upper_bound.length == n
    @upper_bound = upper_bound 
  end
  0.upto(n-1) do |i|
    raise "Lower bounds should be smaller than upper bounds" unless @lower_bound[i] < @upper_bound[i]
  end

  @iterations  = 0
  @evaluations = 0
  # create the simplex for the first time
  build_simplex(start_point)
  evaluate_simplex
end

Instance Attribute Details

#epsilon ⇒ Object (readonly)

Returns the value of attribute epsilon.



62
63
64
 
# File 'lib/mixed_models/NelderMeadWithConstraints.rb', line 62

def epsilon
  @epsilon
end

#f_minimum ⇒ Object (readonly)

Returns the value of attribute f_minimum.



62
63
64
 
# File 'lib/mixed_models/NelderMeadWithConstraints.rb', line 62

def f_minimum
  @f_minimum
end

#lower_bound ⇒ Object (readonly)

Returns the value of attribute lower_bound.



62
63
64
 
# File 'lib/mixed_models/NelderMeadWithConstraints.rb', line 62

def lower_bound
  @lower_bound
end

#max_iterations ⇒ Object (readonly)

Returns the value of attribute max_iterations.



62
63
64
 
# File 'lib/mixed_models/NelderMeadWithConstraints.rb', line 62

def max_iterations
  @max_iterations
end

#start_point ⇒ Object (readonly)

Returns the value of attribute start_point.



62
63
64
 
# File 'lib/mixed_models/NelderMeadWithConstraints.rb', line 62

def start_point
  @start_point
end

#upper_bound ⇒ Object (readonly)

Returns the value of attribute upper_bound.



62
63
64
 
# File 'lib/mixed_models/NelderMeadWithConstraints.rb', line 62

def upper_bound
  @upper_bound
end

#x_minimum ⇒ Object (readonly)

Returns the value of attribute x_minimum.



62
63
64
 
# File 'lib/mixed_models/NelderMeadWithConstraints.rb', line 62

def x_minimum
  @x_minimum
end

Class Method Details

.minimize(start_point:, lower_bound: nil, upper_bound: nil, epsilon: 1e-6, max_iterations: 1e6, &f) ⇒ Object

Convenience method to minimize

Arguments

  • start_point - an Array specifying the initial point for the minimization

  • lower_bound - an Array of lower bounds for each coordinate of the optimal solution

  • upper_bound - an Array of upper bounds for each coordinate of the optimal solution

  • epsilon - a small number specifying the thresholds for the convergence check: absolute_threshold = epsilon and relative_threshold = 100 * epsilon

  • max_iterations - the maximum number of iterations

  • f - the objective function as a Proc object

Usage

minimizer=MixedModels::NelderMead.minimize(start_point: [0,0]) {|x| (x[0] - 1) ** 2 + (x[1] - 5) ** 2}



373
374
375
376
377
378
379
380
381
382
 
# File 'lib/mixed_models/NelderMeadWithConstraints.rb', line 373

def self.minimize(start_point:, lower_bound: nil, upper_bound: nil, 
                  epsilon: 1e-6, max_iterations: 1e6, &f)
  min=MixedModels::NelderMead.new(start_point: start_point, lower_bound: lower_bound, 
                                  upper_bound: upper_bound, epsilon: epsilon, 
                                  max_iterations: max_iterations, &f)
  while min.converging?
    min.iterate
  end
  return min
end

Instance Method Details

#build_simplex(start_point) ⇒ Object

Build an initial simplex

Arguments

  • start_point - starting point of the minimization search



164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
 
# File 'lib/mixed_models/NelderMeadWithConstraints.rb', line 164

def build_simplex(start_point)
  n = start_point.length
  raise "dimension mismatch" if n != @start_configuration.length
  # set first vertex
  @simplex = Array.new(n+1)
  @simplex[0] = PointValuePair.new(move_into_bounds(start_point), Float::NAN)

  # set remaining vertices
  0.upto(n - 1) do |i|
    conf_i   = @start_configuration[i]
    vertex_i = Array.new(n)
    0.upto(n - 1) do |k|
      vertex_i[k] = start_point[k] + conf_i[k]
    end
    @simplex[i + 1] = PointValuePair.new(move_into_bounds(vertex_i), Float::NAN)
  end
end

#compare(v1, v2) ⇒ Object

compares 2 PointValuePair points



226
227
228
229
230
231
232
233
234
 
# File 'lib/mixed_models/NelderMeadWithConstraints.rb', line 226

def compare(v1, v2)
  if v1.value == v2.value
    return 0
  elsif v1.value > v2.value
    return 1
  else
    return -1
  end
end

#converging? ⇒ Boolean

checks whether the function is converging. Returns true if not converged yet, false when converged.

Returns:

  • (Boolean) 


195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
 
# File 'lib/mixed_models/NelderMeadWithConstraints.rb', line 195

def converging?
  # check the convergence in a given direction comparing the previous and current values
  def point_converged?(previous, current)
    pre        = previous.value
    curr       = current.value
    diff       = (pre - curr).abs
    size       = [pre.abs, curr.abs].max
    return ((diff <= (size * @relative_threshold)) and (diff <= @absolute_threshold))
  end

  # returns true if converging is possible atleast in one direction
  if @iterations > 0
    # given direction is converged
    converged = true
    0.upto(@simplex.length - 1) do |i|
      converged &= point_converged?(@previous[i], @simplex[i])
    end
    return !converged
  end

  # if no iterations were done, convergence undefined
  return true
end

#evaluate_simplex ⇒ Object

Evaluate all the non-evaluated points of the simplex, and sort the points in the simplex from best to worst



184
185
186
187
188
189
190
191
 
# File 'lib/mixed_models/NelderMeadWithConstraints.rb', line 184

def evaluate_simplex
  # evaluate the objective function at all non-evaluated simplex points
  @simplex.each_with_index do |v,i|
    @simplex[i].value = f(v.point) if v.value.nan?
  end
  # sort the simplex from best to worst
  @simplex.sort!{ |x1, x2| x1.value <=> x2.value }
end

#f(x) ⇒ Object 



124
125
126
 
# File 'lib/mixed_models/NelderMeadWithConstraints.rb', line 124

def f(x)
  return @f.call(x)
end

#increment_iterations_counter ⇒ Object

increment iteration counter by 1



220
221
222
223
 
# File 'lib/mixed_models/NelderMeadWithConstraints.rb', line 220

def increment_iterations_counter
  @iterations += 1
  raise "iteration limit reached" if @iterations > @max_iterations
end

#iterate ⇒ Object

Iterate the simplex one step. Use this when iteration needs to be done manually

Usage

minimizer=MixedModels::NelderMead.new(start_point: [0,0]) {|x| (x[0] - 1) ** 2 + (x[1] - 5) ** 2}
while minimizer.converging?
  minimizer.iterate
end
minimizer.x_minimum
minimizer.f_minimum



347
348
349
350
351
352
353
354
355
356
 
# File 'lib/mixed_models/NelderMeadWithConstraints.rb', line 347

def iterate
  # set previous simplex as the current simplex
  @previous = Array.new(@simplex.length)
  @simplex.each_with_index { |v,i| @previous[i] = PointValuePair.new(v.point, v.value) }
  # iterate simplex
  iterate_simplex
  # set results
  @x_minimum = @simplex[0].point
  @f_minimum = @simplex[0].value
end

#iterate_simplex ⇒ Object 



251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
 
# File 'lib/mixed_models/NelderMeadWithConstraints.rb', line 251

def iterate_simplex
  increment_iterations_counter
  
  # the simplex has n+1 point if dimension is n
  n           = @simplex.length - 1
  best        = @simplex[0]
  secondWorst = @simplex[n - 1]
  worst       = @simplex[n]
  x_worst     = worst.point
  
  # compute the centroid of the best vertices
  # (dismissing the worst point at index n)
  centroid   = Array.new(n, 0)
  0.upto(n - 1) do |i|
    x = @simplex[i].point
    0.upto(n - 1) { |j| centroid[j] += x[j] }
  end
  scaling = 1.0 / n
  0.upto(n - 1) { |j| centroid[j] *= scaling }

  # compute the reflection point
  xr = Array.new(n)
  0.upto(n - 1) do |j|
    xr[j] = centroid[j] + @rho * (centroid[j] - x_worst[j])
  end
  xr = move_into_bounds(xr)
  reflected = PointValuePair.new(xr, f(xr))
  if ((compare(best, reflected) <= 0) && (compare(reflected, secondWorst) < 0))
    # accept the reflected point
    replace_worst_point(reflected)
  elsif (compare(reflected, best) < 0)
    # compute the expansion point
    xe = Array.new(n)
    0.upto(n - 1) do |j|
      xe[j] = centroid[j] + @khi * (xr[j] - centroid[j])
    end
    xe = move_into_bounds(xe)
    expanded = PointValuePair.new(xe, f(xe))
    if (compare(expanded, reflected) < 0)
      # accept the expansion point
      replace_worst_point(expanded)
    else
      # accept the reflected point
      replace_worst_point(reflected)
    end
  else
    if (compare(reflected, worst) < 0)
      # perform an outside contraction
      xc = Array.new(n)
      0.upto(n - 1) do |j|
        xc[j] = centroid[j] + @gamma * (xr[j] - centroid[j])
      end
      xc = move_into_bounds(xc)
      out_contracted = PointValuePair.new(xc, f(xc))
      if (compare(out_contracted, reflected) <= 0)
        # accept the contraction point
        replace_worst_point(out_contracted)
        return
      end
    else
      # perform an inside contraction
      xc = Array.new(n)
      0.upto(n - 1) do |j|
        xc[j] = centroid[j] + @gamma * (x_worst[j] - centroid[j])
      end
      xc = move_into_bounds(xc)
      in_contracted = PointValuePair.new(xc, f(xc))
      if (compare(in_contracted, worst) < 0)
        # accept the contraction point
        replace_worst_point(in_contracted)
        return
      end
    end
    # if contraction failed, perform a shrink
    x_smallest = @simplex[0].point
    0.upto(n) do |i|
      x = @simplex[i].get_point_clone
      0.upto(n - 1) do |j|
        x[j] = x_smallest[j] + @sigma * (x[j] - x_smallest[j])
      end
      @simplex[i] = PointValuePair.new(x, Float::NAN)
    end
    evaluate_simplex
  end
end

#move_into_bounds(point) ⇒ Object

Check if a given point is within the bounds given by @lower_bound and @upper_bound, and if that’s not the case then move the point inside the bounded region. The returned value is the shifted point if it was necessary to move it (otherwise the originally supplied point).

Arguments

  • point - an array with the coordinates of the point



149
150
151
152
153
154
155
156
157
 
# File 'lib/mixed_models/NelderMeadWithConstraints.rb', line 149

def move_into_bounds(point)
  n = point.length
  raise "dimension mismatch" if n != @start_configuration.length
  0.upto(n-1) do |i|
    point[i] = @lower_bound[i] if @lower_bound[i] > point[i]
    point[i] = @upper_bound[i] if @upper_bound[i] < point[i]
  end
  return point
end

#point_converged?(previous, current) ⇒ Boolean

check the convergence in a given direction comparing the previous and current values

Returns:

  • (Boolean) 


197
198
199
200
201
202
203
 
# File 'lib/mixed_models/NelderMeadWithConstraints.rb', line 197

def point_converged?(previous, current)
  pre        = previous.value
  curr       = current.value
  diff       = (pre - curr).abs
  size       = [pre.abs, curr.abs].max
  return ((diff <= (size * @relative_threshold)) and (diff <= @absolute_threshold))
end

#replace_worst_point(point_value_pair) ⇒ Object

Replace the worst point of the simplex by a new point

Arguments

  • point_value_pair - point to insert



241
242
243
244
245
246
247
248
249
 
# File 'lib/mixed_models/NelderMeadWithConstraints.rb', line 241

def replace_worst_point(point_value_pair)
  n = @simplex.length - 1
  0.upto(n - 1) do |i|
    if (compare(@simplex[i], point_value_pair) > 0)
      point_value_pair, @simplex[i] = @simplex[i], point_value_pair
    end
  end
  @simplex[n] = point_value_pair
end

#start_configuration=(steps) ⇒ Object

only the relative position of the n vertices with respect to the first one are stored



130
131
132
133
134
135
136
137
138
 
# File 'lib/mixed_models/NelderMeadWithConstraints.rb', line 130

def start_configuration=(steps)
  n = steps.length
  @start_configuration = Array.new(n) { Array.new(n, 0) }
  0.upto(n - 1) do |i|
    vertex_i = @start_configuration[i]
    raise "equals vertices #{i-1} and #{i} in simplex configuration" if steps[i] == 0.0
    0.upto(i) { |j| vertex_i[j] = steps[j] }
  end
end