Class: Minimization::NonLinearConjugateGradientMinimizer

Inherits:
Object
  • Object
show all
Defined in:
lib/multidim/conjugate_gradient.rb

Overview

Conjugate Gradient minimizer class The beta function may be :fletcher_reeves or :polak_ribiere

Direct Known Subclasses

FletcherReeves, PolakRibiere

Constant Summary collapse

MAX_ITERATIONS_DEFAULT = 
100000
EPSILON_DEFAULT = 
1e-6

Instance Attribute Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(f, fd, start_point, beta_formula) ⇒ NonLinearConjugateGradientMinimizer

Returns a new instance of NonLinearConjugateGradientMinimizer.



46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
 
# File 'lib/multidim/conjugate_gradient.rb', line 46

def initialize(f, fd, start_point, beta_formula)
  @epsilon     = EPSILON_DEFAULT
  @safe_min    = 4.503599e15
  @f           = f
  @fd          = fd
  @start_point = start_point

  @max_iterations     = MAX_ITERATIONS_DEFAULT
  @iterations         = 0
  @update_formula     = beta_formula
  @relative_threshold = 100 * @epsilon
  @absolute_threshold = 100 * @safe_min

  @initial_step = 1.0 # initial step default
  @converging   = true

  # do initial steps
  @point = @start_point.clone
  @n      = @point.length
  @r      = gradient(@point)
  0.upto(@n - 1) do |i|
    @r[i] = -@r[i]
  end
  
  # Initial search direction.
  @steepest_descent = precondition(@point, @r)
  @search_direction = @steepest_descent.clone

  @delta = 0
  0.upto(@n - 1) do |i|
      @delta += @r[i] * @search_direction[i]
  end
  @current = nil
end

Instance Attribute Details

#converging ⇒ Object (readonly) Also known as: converging?

Returns the value of attribute converging.



37
38
39
 
# File 'lib/multidim/conjugate_gradient.rb', line 37

def converging
  @converging
end

#f_minimum ⇒ Object (readonly)

Returns the value of attribute f_minimum.



36
37
38
 
# File 'lib/multidim/conjugate_gradient.rb', line 36

def f_minimum
  @f_minimum
end

#initial_step ⇒ Object

Returns the value of attribute initial_step.



39
40
41
 
# File 'lib/multidim/conjugate_gradient.rb', line 39

def initial_step
  @initial_step
end

#x_minimum ⇒ Object (readonly)

Returns the value of attribute x_minimum.



35
36
37
 
# File 'lib/multidim/conjugate_gradient.rb', line 35

def x_minimum
  @x_minimum
end

Instance Method Details

#converged(previous, current) ⇒ Object 



110
111
112
113
114
115
116
 
# File 'lib/multidim/conjugate_gradient.rb', line 110

def converged(previous, current)
  p          = f(previous)
  c          = f(current)
  difference = (p - c).abs
  size       = [p.abs, c.abs].max
  return ((difference <= size * @relative_threshold) or (difference <= @absolute_threshold))
end

#f(x) ⇒ Object 



81
82
83
 
# File 'lib/multidim/conjugate_gradient.rb', line 81

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

#find_upper_bound(a, h, search_direction) ⇒ Object 



89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
 
# File 'lib/multidim/conjugate_gradient.rb', line 89

def find_upper_bound(a, h, search_direction)
  ya   = line_search_func(a, search_direction).to_f
  yb   = ya
  step = h
  # check step value for float max value exceeds
  while step < Float::MAX
    b  = a + step
    yb = line_search_func(b, search_direction).to_f
    if (ya * yb <= 0)
      return b
    end
    step *= [2, ya / yb].max
  end
  # raise error if bracketing failed
  raise "Unable to bracket minimum in line search."
end

#gradient(x) ⇒ Object 



85
86
87
 
# File 'lib/multidim/conjugate_gradient.rb', line 85

def gradient(x)
  return @fd.call(x)
end

#iterate ⇒ Object

iterate one step of conjugate gradient minimizer

Usage:

f = proc{ |x| (x - 2)**2 + (x - 5)**2 + (x - 100)**2 } fd = proc{ |x| [ 2 * (x[0] - 2) , 2 * (x[1] - 5) , 2 * (x[2] - 100) ] } min = Minimization::FletcherReeves.new(f, fd, [0, 0, 0]) while(min.converging?) min.iterate end min.x_minimum min.f_minimum



157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
 
# File 'lib/multidim/conjugate_gradient.rb', line 157

def iterate
  @iterations  += 1
  @previous     = @current
  @current      = Minimization::PointValuePair.new(@point, f(@point))
  # set converging parameter
  @converging   = !(@previous != nil and converged(@previous.point, @current.point))
  # set results
  @x_minimum    = @current.point
  @f_minimum    = @current.value

  # set search_direction to be used in solve and find_upper_bound methods
  ub   = find_upper_bound(0, @initial_step, @search_direction)
  step = solve(0, ub, 1e-15, @search_direction)

  # Validate new point
  0.upto(@point.length - 1) do |i|
    @point[i] += step * @search_direction[i]
  end

  @r = gradient(@point)
  0.upto(@n - 1) do |i|
    @r[i] = -@r[i]
  end

  # Compute beta
  delta_old            = @delta
  new_steepest_descent = precondition(@point, @r)
  @delta                = 0
  0.upto(@n - 1) do |i|
    @delta += @r[i] * new_steepest_descent[i]
  end

  if (@update_formula == :fletcher_reeves)
    beta = @delta.to_f / delta_old
  elsif(@update_formula == :polak_ribiere)
    deltaMid = 0
    0.upto(@r.length - 1) do |i|
      deltaMid += @r[i] * @steepest_descent[i]
    end
    beta = (@delta - deltaMid).to_f / delta_old
  else
    raise "Unknown beta formula type"
  end
  @steepest_descent = new_steepest_descent

  # Compute conjugate search direction
  if ((@iterations % @n == 0) or (beta < 0))
    # Break conjugation: reset search direction
    @search_direction = @steepest_descent.clone
  else
    # Compute new conjugate search direction
    0.upto(@n - 1) do |i|
      @search_direction[i] = @steepest_descent[i] + beta * @search_direction[i]
    end
  end
end

#line_search_func(x, search_direction) ⇒ Object 



127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
 
# File 'lib/multidim/conjugate_gradient.rb', line 127

def line_search_func(x, search_direction)
  # current point in the search direction
  shifted_point = @point.clone
  0.upto(shifted_point.length - 1) do |i|
    shifted_point[i] += x * search_direction[i]
  end

  # gradient of the objective function
  gradient = gradient(shifted_point)

  # dot product with the search direction
  dot_product = 0
  0.upto(gradient.length - 1) do |i|
    dot_product += gradient[i] * search_direction[i]
  end

  return dot_product
end

#precondition(point, r) ⇒ Object 



106
107
108
 
# File 'lib/multidim/conjugate_gradient.rb', line 106

def precondition(point, r)
  return r.clone # case: identity preconditioner has been used as the default
end

#solve(min, max, start_value, search_direction) ⇒ Object

solver to use during line search



119
120
121
122
123
124
125
 
# File 'lib/multidim/conjugate_gradient.rb', line 119

def solve(min, max, start_value, search_direction)
  # check start_value to eliminate unnessasary calculations ...
  func        = proc{|x| line_search_func(x, search_direction)}
  root_finder = Minimization::BrentRootFinder.new(func)
  root        = root_finder.find_root(min, max, func)
  return root
end