Class: Minimization::Brent

Inherits:
Unidimensional show all
Defined in:
lib/minimization.rb

Overview

Direct port of Brent algorithm found on GSL. See Unidimensional for methods.

Usage

min=Minimization::Brent.new(-1000,20000  , proc {|x| (x+1)**2}
min.expected=1.5  # Expected value
min.iterate
min.x_minimum
min.f_minimum
min.log

Constant Summary collapse

GSL_SQRT_DBL_EPSILON =
1.4901161193847656e-08

Constants inherited from Unidimensional

Unidimensional::EPSILON, Unidimensional::MAX_ITERATIONS

Instance Attribute Summary

Attributes inherited from Unidimensional

#epsilon, #expected, #f_minimum, #iterations, #log, #log_header, #x_minimum

Instance Method Summary collapse

Methods inherited from Unidimensional

#f, #log_summary, minimize

Constructor Details

#initialize(lower, upper, proc) ⇒ Brent


217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
# File 'lib/minimization.rb', line 217

def initialize(lower,upper, proc)
  super

  @do_bracketing=true

  # Init

  golden = 0.3819660;      #golden = (3 - sqrt(5))/2

  v = @lower + golden * (@upper - @lower);
  w = v;

  @x_minimum = v ;
  @f_minimum = f(v) ;
  @x_lower=@lower
  @x_upper=@upper
  @f_lower = f(@lower) ;
  @f_upper = f(@lower) ;

  @v = v;
  @w = w;

  @d = 0;
  @e = 0;
  @f_v=f(v)
  @f_w=@f_v
end

Instance Method Details

#bracketingObject


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
# File 'lib/minimization.rb', line 251

def bracketing
  eval_max=10
  f_left = @f_lower;
  f_right = @f_upper;
  x_left = @x_lower;
  x_right= @x_upper;
  golden = 0.3819660;      # golden = (3 - sqrt(5))/2 */
  nb_eval=0

  if (f_right >= f_left)
    x_center = (x_right - x_left) * golden + x_left;
    nb_eval+=1;
    f_center=f(x_center)
  else
    x_center = x_right ;
    f_center = f_right ;
    x_right = (x_center - x_left).quo(golden) + x_left;
    nb_eval+=1;
    f_right=f(x_right);
  end


  begin
    @log << ["B#{nb_eval}", x_left, x_right, f_left, f_right, (x_left-x_right).abs, (f_left-f_right).abs]
    if (f_center < f_left )
      if (f_center < f_right)
        @x_lower = x_left;
        @x_upper = x_right;
        @x_minimum = x_center;
        @f_lower = f_left;
        @f_upper = f_right;
        @f_minimum = f_center;
        return true;
      elsif (f_center > f_right)
        x_left = x_center;
        f_left = f_center;
        x_center = x_right;
        f_center = f_right;
        x_right = (x_center - x_left).quo(golden) + x_left;
        nb_eval+=1;
        f_right=f(x_right);
      else # f_center == f_right */
        x_right = x_center;
        f_right = f_center;
        x_center = (x_right - x_left).quo(golden) + x_left;
        nb_eval+=1;
        f_center=f(x_center);
      end
    else # f_center >= f_left */
      x_right = x_center;
      f_right = f_center;
      x_center = (x_right - x_left) * golden + x_left;
      nb_eval+=1;
      f_center=f(x_center);
    end
  end while ((nb_eval < eval_max) and
  ((x_right - x_left) > GSL_SQRT_DBL_EPSILON * ( (x_right + x_left) * 0.5 ) + GSL_SQRT_DBL_EPSILON))
  @x_lower = x_left;
  @x_upper = x_right;
  @x_minimum = x_center;
  @f_lower = f_left;
  @f_upper = f_right;
  @f_minimum = f_center;
  return false;

end

#brent_iterateObject

Generate one iteration.


336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
# File 'lib/minimization.rb', line 336

def brent_iterate
  x_left = @x_lower;
  x_right = @x_upper;

  z = @x_minimum;
  d = @e;
  e = @d;
  v = @v;
  w = @w;
  f_v = @f_v;
  f_w = @f_w;
  f_z = @f_minimum;

  golden = 0.3819660;      # golden = (3 - sqrt(5))/2 */

  w_lower = (z - x_left)
  w_upper = (x_right - z)

  tolerance =  GSL_SQRT_DBL_EPSILON * z.abs

  midpoint = 0.5 * (x_left + x_right)
  _p,q,r=0,0,0
  if (e.abs > tolerance)

    # fit parabola */

    r = (z - w) * (f_z - f_v);
    q = (z - v) * (f_z - f_w);
    _p = (z - v) * q - (z - w) * r;
    q = 2 * (q - r);

    if (q > 0)
      _p = -_p
    else
      q = -q;
    end
    r = e;
    e = d;
  end

  if (_p.abs < (0.5 * q * r).abs and _p < q * w_lower and _p < q * w_upper)
    t2 = 2 * tolerance ;

    d = _p.quo(q);
    u = z + d;

    if ((u - x_left) < t2 or (x_right - u) < t2)
      d = (z < midpoint) ? tolerance : -tolerance ;
    end
  else

    e = (z < midpoint) ? x_right - z : -(z - x_left) ;
    d = golden * e;
  end

  if ( d.abs >= tolerance)
    u = z + d;
  else
    u = z + ((d > 0) ? tolerance : -tolerance) ;
  end

  @e = e;
  @d = d;

  f_u=f(u)

  if (f_u <= f_z)
    if (u < z)
      @x_upper = z;
      @f_upper = f_z;
    else
      @x_lower = z;
      @f_lower = f_z;
    end
    @v = w;
    @f_v = f_w;
    @w = z;
    @f_w = f_z;
    @x_minimum = u;
    @f_minimum = f_u;
    return true;
  else
    if (u < z)
      @x_lower = u;
      @f_lower = f_u;
      return true;
    else
      @x_upper = u;
      @f_upper = f_u;
      return true;
    end

    if (f_u <= f_w or w == z)
      @v = w;
      @f_v = f_w;
      @w = u;
      @f_w = f_u;
      return true;
    elsif f_u <= f_v or v == z or v == w
      @v = u;
      @f_v = f_u;
      return true;
    end

  end
  return false

end

#expected=(v) ⇒ Object


245
246
247
248
249
# File 'lib/minimization.rb', line 245

def expected=(v)
  @x_minimum=v
  @f_minimum=f(v)
  @do_bracketing=false
end

#iterateObject

Start the minimization process If you want to control manually the process, use brent_iterate


319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
# File 'lib/minimization.rb', line 319

def iterate
  k=0
  bracketing if @do_bracketing
  while k<@max_iteration and (@x_lower-@x_upper).abs>@epsilon
    k+=1
    result=brent_iterate
    raise FailedIteration,"Error on iteration" if !result
    begin
      @log << [k, @x_lower, @x_upper, @f_lower, @f_upper, (@x_lower-@x_upper).abs, (@f_lower-@f_upper).abs]
    rescue =>@e
      @log << [k, @e.to_s,nil,nil,nil,nil,nil]
    end
  end
  @iterations=k
  return true
end