Source Code for Module minfx.newton

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2003-2013 Edward d'Auvergne                                   # 
  4  #                                                                             # 
  5  # This file is part of the minfx optimisation library,                        # 
  6  # https://sourceforge.net/projects/minfx                                      # 
  7  #                                                                             # 
  8  # This program is free software: you can redistribute it and/or modify        # 
  9  # it under the terms of the GNU General Public License as published by        # 
 10  # the Free Software Foundation, either version 3 of the License, or           # 
 11  # (at your option) any later version.                                         # 
 12  #                                                                             # 
 13  # This program is distributed in the hope that it will be useful,             # 
 14  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 15  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 16  # GNU General Public License for more details.                                # 
 17  #                                                                             # 
 18  # You should have received a copy of the GNU General Public License           # 
 19  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 20  #                                                                             # 
 21  ############################################################################### 
 22   
 23  # Module docstring. 
 24  """Newton (or Newton-Raphson - NR) optimization. 
 25   
 26  This file is part of the U{minfx optimisation library<https://sourceforge.net/projects/minfx>}. 
 27  """ 
 28   
 29  # Python module imports. 
 30  from math import acos 
 31  from numpy import dot, float64, identity, sqrt, zeros 
 32  from re import match 
 33   
 34  # Minfx module imports. 
 35  from minfx.base_classes import Hessian_mods, Line_search, Min 
 36   
 37   
 38 -def newton(func=None, dfunc=None, d2func=None, args=(), x0=None, min_options=(), func_tol=1e-25, grad_tol=None, maxiter=1e6, a0=1.0, mu=0.0001, eta=0.9, mach_acc=1e-16, full_output=0, print_flag=0, print_prefix=""): 
 39      """Newton minimisation.""" 
 40   
 41      if print_flag: 
 42          if print_flag >= 2: 
 43              print(print_prefix) 
 44          print(print_prefix) 
 45          print(print_prefix + "Newton minimisation") 
 46          print(print_prefix + "~~~~~~~~~~~~~~~~~~~") 
 47      min = Newton(func, dfunc, d2func, args, x0, min_options, func_tol, grad_tol, maxiter, a0, mu, eta, mach_acc, full_output, print_flag, print_prefix) 
 48      if min.init_failure: 
 49          print(print_prefix + "Initialisation of minimisation has failed.") 
 50          return None 
 51      results = min.minimise() 
 52      return results 
 53   
 54   
 55 -class Newton(Hessian_mods, Line_search, Min): 
 56 -    def __init__(self, func, dfunc, d2func, args, x0, min_options, func_tol, grad_tol, maxiter, a0, mu, eta, mach_acc, full_output, print_flag, print_prefix): 
 57          """Class for Newton minimisation specific functions. 
 58   
 59          Unless you know what you are doing, you should call the function 'newton' rather than using this class. 
 60          """ 
 61   
 62          # Function arguments. 
 63          self.func = func 
 64          self.dfunc = dfunc 
 65          self.d2func = d2func 
 66          self.args = args 
 67          self.xk = x0 
 68          self.func_tol = func_tol 
 69          self.grad_tol = grad_tol 
 70          self.maxiter = maxiter 
 71          self.mach_acc = mach_acc 
 72          self.full_output = full_output 
 73          self.print_flag = print_flag 
 74          self.print_prefix = print_prefix 
 75   
 76          # Set a0. 
 77          self.a0 = a0 
 78   
 79          # Line search constants for the Wolfe conditions. 
 80          self.mu = mu 
 81          self.eta = eta 
 82   
 83          # Initialisation failure flag. 
 84          self.init_failure = 0 
 85   
 86          # Setup the line search and Hessian modification algorithms. 
 87          self.line_search_algor = None 
 88          self.hessian_mod = None 
 89   
 90          # Test if the options are a tuple. 
 91          if not isinstance(min_options, tuple): 
 92              print(self.print_prefix + "The minimisation options " + repr(min_options) + " is not a tuple.") 
 93              self.init_failure = 1; return 
 94   
 95          # Test that no more thant 2 options are given. 
 96          if len(min_options) > 2: 
 97              print(self.print_prefix + "A maximum of two minimisation options is allowed (the line search algorithm and the Hessian modification).") 
 98              self.init_failure = 1; return 
 99   
100          # Sort out the minimisation options. 
101          none_option = None 
102          for opt in min_options: 
103              if match('[Nn]one', opt): 
104                  none_option = opt 
105                  continue 
106              if self.line_search_algor == None and self.valid_line_search(opt): 
107                  self.line_search_algor = opt 
108              elif self.hessian_mod == None and self.valid_hessian_mod(opt): 
109                  self.hessian_mod = opt 
110              else: 
111                  if self.line_search_algor: 
112                      print(self.print_prefix + "The minimisation option " + repr(opt) + " from " + repr(min_options) + " is not a valid Hessian modification.") 
113                  elif self.hessian_mod: 
114                      print(self.print_prefix + "The minimisation option " + repr(opt) + " from " + repr(min_options) + " is not a valid line search algorithm.") 
115                  else: 
116                      print(self.print_prefix + "The minimisation option " + repr(opt) + " from " + repr(min_options) + " is neither a valid line search algorithm or Hessian modification.") 
117                  self.init_failure = 1; return 
118          if none_option and not self.hessian_mod: 
119              if self.valid_hessian_mod(none_option): 
120                  self.hessian_mod = none_option 
121          elif none_option and not self.line_search_algor: 
122              if self.valid_line_search(none_option): 
123                  self.line_search_algor = none_option 
124          elif none_option: 
125              print(self.print_prefix + "Invalid combination of options " + repr(min_options) + ".") 
126              self.init_failure = 1; return 
127   
128          # Default line search algorithm. 
129          if self.line_search_algor == None: 
130              self.line_search_algor = 'Back' 
131   
132          # Default Hessian modification. 
133          if self.hessian_mod == None: 
134              self.hessian_mod = 'GMW' 
135   
136          # Line search and Hessian modification initialisation. 
137          self.setup_line_search() 
138          self.setup_hessian_mod() 
139   
140          # Initialise the function, gradient, and Hessian evaluation counters. 
141          self.f_count = 0 
142          self.g_count = 0 
143          self.h_count = 0 
144   
145          # Initialise the warning string. 
146          self.warning = None 
147   
148          # Constants. 
149          self.n = len(self.xk) 
150          self.I = identity(len(self.xk)) 
151   
152          # Set the convergence test function. 
153          self.setup_conv_tests() 
154   
155          # Newton setup function. 
156          self.setup_newton() 
157   
158          # Set the update function. 
159          self.update = self.update_newton 
160   
161   
162 -    def new_param_func(self): 
163          """The new parameter function. 
164   
165          Find the search direction, do a line search, and get xk+1 and fk+1. 
166          """ 
167   
168          # Calculate the Newton direction. 
169          self.pk = self.get_pk() 
170   
171          # Line search. 
172          self.line_search() 
173   
174          # Find the new parameter vector and function value at that point. 
175          self.xk_new = self.xk + self.alpha * self.pk 
176          self.fk_new, self.f_count = self.func(*(self.xk_new,)+self.args), self.f_count + 1 
177          self.dfk_new, self.g_count = self.dfunc(*(self.xk_new,)+self.args), self.g_count + 1 
178   
179          # Debugging. 
180          if self.print_flag >= 2: 
181              print("\n" + self.print_prefix + "New param function.") 
182              print(self.print_prefix + "pk:    " + repr(self.pk)) 
183              print(self.print_prefix + "alpha: " + repr(self.alpha)) 
184              print(self.print_prefix + "xk:    " + repr(self.xk)) 
185              print(self.print_prefix + "xk+1:  " + repr(self.xk_new)) 
186              print(self.print_prefix + "fk:    " + repr(self.fk)) 
187              print(self.print_prefix + "fk+1:  " + repr(self.fk_new)) 
188              print(self.print_prefix + "dfk:    " + repr(self.dfk)) 
189              print(self.print_prefix + "dfk+1:  " + repr(self.dfk_new)) 
190              print(self.print_prefix + "d2fk:\n" + repr(self.d2fk)) 
191   
192              #eigen = eigenvectors(self.d2fk) 
193              #print(self.print_prefix + "Eigenvalues: " + repr(eigen[0])) 
194   
195              print(self.print_prefix + "Angle to the unit vector pointing along the first dimension.") 
196              unit_vect = zeros(self.n, float64) 
197              unit_vect[0] = 1.0 
198              dfk_angle = acos(dot(self.dfk, unit_vect) / sqrt(dot(self.dfk, self.dfk))) 
199              pk_angle = acos(dot(self.pk, unit_vect) / sqrt(dot(self.pk, self.pk))) 
200              print(self.print_prefix + "steepest descent: " + repr(dfk_angle)) 
201              print(self.print_prefix + "Newton step:      " + repr(pk_angle)) 
202   
203   
204 -    def setup_newton(self): 
205          """Setup function. 
206   
207          The initial Newton function value, gradient vector, and Hessian matrix are calculated. 
208          """ 
209   
210          self.fk, self.f_count = self.func(*(self.xk,)+self.args), self.f_count + 1 
211          self.dfk, self.g_count = self.dfunc(*(self.xk,)+self.args), self.g_count + 1 
212          self.d2fk, self.h_count = self.d2func(*(self.xk,)+self.args), self.h_count + 1 
213   
214   
215 -    def update_newton(self): 
216          """Function to update the function value, gradient vector, and Hessian matrix.""" 
217   
218          self.xk = self.xk_new * 1.0 
219          self.fk = self.fk_new 
220          self.dfk = self.dfk_new * 1.0 
221          self.d2fk, self.h_count = self.d2func(*(self.xk,)+self.args), self.h_count + 1 
222