Source Code for Module minfx.exact_trust_region

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2003-2013 Edward d'Auvergne                                   # 
  4  #                                                                             # 
  5  # This file is part of the minfx optimisation library,                        # 
  6  # https://gna.org/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  """Exact trust region optimization. 
 25   
 26  This file is part of the U{minfx optimisation library<https://gna.org/projects/minfx>}. 
 27  """ 
 28   
 29  # Python module imports. 
 30  from numpy import dot, identity, sort, sqrt, transpose 
 31  from numpy.linalg import cholesky, eig, inv, solve 
 32  from re import match 
 33   
 34  # Minfx module imports. 
 35  from minfx.base_classes import Hessian_mods, Min, Trust_region 
 36  from minfx.bfgs import Bfgs 
 37  from minfx.newton import Newton 
 38   
 39   
 40 -def exact_trust_region(func=None, dfunc=None, d2func=None, args=(), x0=None, min_options=(), func_tol=1e-25, grad_tol=None, maxiter=1e6, lambda0=0.0, delta_max=1e5, delta0=1.0, eta=0.2, mach_acc=1e-16, full_output=0, print_flag=0, print_prefix=""): 
 41      """Exact trust region algorithm.""" 
 42   
 43      if print_flag: 
 44          if print_flag >= 2: 
 45              print(print_prefix) 
 46          print(print_prefix) 
 47          print(print_prefix + "Exact trust region minimisation") 
 48          print(print_prefix + "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~") 
 49      min = Exact_trust_region(func, dfunc, d2func, args, x0, min_options, func_tol, grad_tol, maxiter, lambda0, delta_max, delta0, eta, mach_acc, full_output, print_flag, print_prefix) 
 50      if min.init_failure: 
 51          print(print_prefix + "Initialisation of minimisation has failed.") 
 52          return None 
 53      results = min.minimise() 
 54      return results 
 55   
 56   
 57 -class Exact_trust_region(Trust_region, Bfgs, Newton): 
 58 -    def __init__(self, func, dfunc, d2func, args, x0, min_options, func_tol, grad_tol, maxiter, lambda0, delta_max, delta0, eta, mach_acc, full_output, print_flag, print_prefix): 
 59          """Class for Exact trust region minimisation specific functions. 
 60   
 61          Unless you know what you are doing, you should call the function 'exact_trust_region' rather than using this class. 
 62          """ 
 63   
 64          # Function arguments. 
 65          self.func = func 
 66          self.dfunc = dfunc 
 67          self.d2func = d2func 
 68          self.args = args 
 69          self.xk = x0 
 70          self.func_tol = func_tol 
 71          self.grad_tol = grad_tol 
 72          self.maxiter = maxiter 
 73          self.full_output = full_output 
 74          self.print_flag = print_flag 
 75          self.print_prefix = print_prefix 
 76          self.mach_acc = mach_acc 
 77          self.lambda0 = lambda0 
 78          self.delta_max = delta_max 
 79          self.delta = delta0 
 80          self.eta = eta 
 81   
 82          # Initialisation failure flag. 
 83          self.init_failure = 0 
 84   
 85          # Setup the Hessian modification options and algorithm. 
 86          self.hessian_type_and_mod(min_options) 
 87          self.setup_hessian_mod() 
 88   
 89          # Initialise the function, gradient, and Hessian evaluation counters. 
 90          self.f_count = 0 
 91          self.g_count = 0 
 92          self.h_count = 0 
 93   
 94          # Initialise the warning string. 
 95          self.warning = None 
 96   
 97          # Constants. 
 98          self.n = len(self.xk) 
 99          self.I = identity(len(self.xk)) 
100   
101          # Set the convergence test function. 
102          self.setup_conv_tests() 
103   
104          # Type specific functions. 
105          if match('[Bb][Ff][Gg][Ss]', self.hessian_type): 
106              self.setup_bfgs() 
107              self.specific_update = self.update_bfgs 
108              self.d2fk = inv(self.Hk) 
109              self.warning = "Incomplete code." 
110          elif match('[Nn]ewton', self.hessian_type): 
111              self.setup_newton() 
112              self.specific_update = self.update_newton 
113   
114   
115 -    def new_param_func(self): 
116          """Find the exact trust region solution. 
117   
118          Algorithm 4.4 from page 81 of 'Numerical Optimization' by Jorge Nocedal and Stephen J. Wright, 1999, 2nd ed.  This is only implemented for positive definite matrices. 
119          """ 
120   
121          # Matrix modification. 
122          pB, matrix = self.get_pk(return_matrix=1) 
123   
124          # The exact trust region algorithm. 
125          l = 0 
126          lambda_l = self.lambda0 
127   
128          while l < 3: 
129              # Calculate the matrix B + lambda(l).I 
130              B = matrix + lambda_l * self.I 
131   
132              # Factor B + lambda(l).I = RT.R 
133              R = cholesky(B) 
134   
135              # Solve RT.R.pl = -g 
136              y = solve(R, self.dfk) 
137              y = -solve(transpose(R), y) 
138              dot_pl = dot(y, y) 
139   
140              # Solve RT.ql = pl 
141              y = solve(transpose(R), y) 
142   
143              # lambda(l+1) update. 
144              lambda_l = lambda_l + (dot_pl / dot(y, y)) * ((sqrt(dot_pl) - self.delta) / self.delta) 
145   
146              # Safeguard. 
147              lambda_l = max(0.0, lambda_l) 
148   
149              if self.print_flag >= 2: 
150                  print(self.print_prefix + "\tl: " + repr(l) + ", lambda(l) fin: " + repr(lambda_l)) 
151   
152              l = l + 1 
153   
154          # Calculate the step. 
155          R = cholesky(matrix + lambda_l * self.I) 
156          y = solve(R, self.dfk) 
157          self.pk = -solve(transpose(R), y) 
158   
159          if self.print_flag >= 2: 
160              print(self.print_prefix + "Step: " + repr(self.pk)) 
161   
162          # Find the new parameter vector and function value at that point. 
163          self.xk_new = self.xk + self.pk 
164          self.fk_new, self.f_count = self.func(*(self.xk_new,)+self.args), self.f_count + 1 
165          self.dfk_new, self.g_count = self.dfunc(*(self.xk_new,)+self.args), self.g_count + 1 
166   
167   
168 -    def old_param_func(self): 
169          """Find the exact trust region solution. 
170   
171          More, J. J., and Sorensen D. C. 1983, Computing a trust region step.  SIAM J. Sci. Stat. Comput. 4, 553-572.  This function is incomplete. 
172          """ 
173   
174          self.warning = "Incomplete code, minimisation bypassed." 
175          print(self.print_prefix + "Incomplete code, minimisation bypassed.") 
176          return 
177   
178          # Initialisation. 
179          iter = 0 
180          self.l = self.lambda0 
181          self.I = identity(len(self.dfk)) 
182   
183          # Initialise lL, lU, lS. 
184          self.lS = -1e99 
185          b = 0.0 
186          for j in range(len(self.d2fk)): 
187              self.lS = max(self.lS, -self.d2fk[j, j]) 
188              sum = 0.0 
189              for i in range(len(self.d2fk[j])): 
190                  sum = sum + abs(self.d2fk[i, j]) 
191              b = max(b, sum) 
192          a = sqrt(dot(self.dfk, self.dfk)) / self.delta 
193          self.lL = max(0.0, self.lS, a - b) 
194          self.lU = a + b 
195   
196          # Debugging. 
197          if self.print_flag >= 2: 
198              print(self.print_prefix + "Initialisation.") 
199              eigen = eig(self.d2fk) 
200              eigenvals = sort(eigen[0]) 
201              for i in range(len(self.d2fk)): 
202                  print(self.print_prefix + "\tB[" + repr(i) + ", " + repr(i) + "] = " + repr(self.d2fk[i, i])) 
203              print(self.print_prefix + "\tEigenvalues: " + repr(eigenvals)) 
204              print(self.print_prefix + "\t||g||/delta: " + repr(a)) 
205              print(self.print_prefix + "\t||B||1: " + repr(b)) 
206              print(self.print_prefix + "\tl:  " + repr(self.l)) 
207              print(self.print_prefix + "\tlL: " + repr(self.lL)) 
208              print(self.print_prefix + "\tlU: " + repr(self.lU)) 
209              print(self.print_prefix + "\tlS: " + repr(self.lS)) 
210   
211          # Iterative loop. 
212          return 
213          while True: 
214              # Safeguard lambda. 
215              if self.print_flag >= 2: 
216                  print(self.print_prefix + "\n< Iteration " + repr(iter) + " >") 
217                  print(self.print_prefix + "Safeguarding lambda.") 
218                  print(self.print_prefix + "\tInit l: " + repr(self.l)) 
219                  print(self.print_prefix + "\tlL: " + repr(self.lL)) 
220                  print(self.print_prefix + "\tlU: " + repr(self.lU)) 
221                  print(self.print_prefix + "\tlS: " + repr(self.lS)) 
222              self.l = max(self.l, self.lL) 
223              self.l = min(self.l, self.lU) 
224              if self.l <= self.lS: 
225                  if self.print_flag >= 2: 
226                      print(self.print_prefix + "\tself.l <= self.lS") 
227                  self.l = max(0.001*self.lU, sqrt(self.lL*self.lU)) 
228              if self.print_flag >= 2: 
229                  print(self.print_prefix + "\tFinal l: " + repr(self.l)) 
230   
231              # Calculate the matrix 'B + lambda.I' and factor 'B + lambda(l).I = RT.R' 
232              matrix = self.d2fk + self.l * self.I 
233              pos_def = 1 
234              if self.print_flag >= 2: 
235                  print(self.print_prefix + "Cholesky decomp.") 
236                  print(self.print_prefix + "\tB + lambda.I: " + repr(matrix)) 
237                  eigen = eig(matrix) 
238                  eigenvals = sort(eigen[0]) 
239                  print(self.print_prefix + "\tEigenvalues: " + repr(eigenvals)) 
240              try: 
241                  func = cholesky 
242                  R = func(matrix) 
243                  if self.print_flag >= 2: 
244                      print(self.print_prefix + "\tCholesky matrix R: " + repr(R)) 
245              except "LinearAlgebraError": 
246                  if self.print_flag >= 2: 
247                      print(self.print_prefix + "\tLinearAlgebraError, matrix is not positive definite.") 
248                  pos_def = 0 
249              if self.print_flag >= 2: 
250                  print(self.print_prefix + "\tPos def: " + repr(pos_def)) 
251   
252              if pos_def: 
253                  # Solve p = -inv(RT.R).g 
254                  p = -dot(inv(matrix), self.dfk) 
255                  if self.print_flag >= 2: 
256                      print(self.print_prefix + "Solve p = -inv(RT.R).g") 
257                      print(self.print_prefix + "\tp: " + repr(p)) 
258   
259                  # Compute tau and z if ||p|| < delta. 
260                  dot_p = dot(p, p) 
261                  len_p = sqrt(dot_p) 
262                  if len_p < self.delta: 
263   
264                      # Calculate z. 
265   
266                      # Calculate tau. 
267                      delta2_len_p2 = self.delta**2 - dot_p 
268                      dot_p_z = dot(p, z) 
269                      tau = delta2_len_p2 / (dot_p_z + sign(dot_p_z) * sqrt(dot_p_z**2 + delta2_len_p2**2)) 
270   
271                      if self.print_flag >= 2: 
272                          print(self.print_prefix + "||p|| < delta") 
273                          print(self.print_prefix + "\tz: " + repr(z)) 
274                          print(self.print_prefix + "\ttau: " + repr(tau)) 
275                  else: 
276                      if self.print_flag >= 2: 
277                          print(self.print_prefix + "||p|| >= delta") 
278                          print(self.print_prefix + "\tNo doing anything???") 
279   
280                  # Solve q = inv(RT).p 
281                  q = dot(inv(transpose(R)), p) 
282   
283                  # Update lL, lU. 
284                  if len_p < self.delta: 
285                      self.lU = min(self.lU, self.l) 
286                  else: 
287                      self.lL = max(self.lL, self.l) 
288   
289                  # lambda update. 
290                  self.l_corr = dot_p / dot(q, q) * ((len_p - self.delta) / self.delta) 
291   
292              else: 
293                  # Update lambda via lambda = lS. 
294                  self.lS = max(self.lS, self.l) 
295   
296                  self.l_corr = -self.l 
297   
298              # Update lL 
299              self.lL = max(self.lL, self.lS) 
300              if self.print_flag >= 2: 
301                  print(self.print_prefix + "Update lL: " + repr(self.lL)) 
302   
303              # Check the convergence criteria. 
304   
305              # lambda update. 
306              self.l = self.l + self.l_corr 
307   
308              iter = iter + 1 
309   
310          # Find the new parameter vector and function value at that point. 
311          self.xk_new = self.xk + self.pk 
312          self.fk_new, self.f_count = self.func(*(self.xk_new,)+self.args), self.f_count + 1 
313          self.dfk_new, self.g_count = self.dfunc(*(self.xk_new,)+self.args), self.g_count + 1 
314   
315   
316 -    def safeguard(self, eigenvals): 
317          """Safeguarding function.""" 
318   
319          if self.lambda_l < -eigenvals[0]: 
320              self.lambda_l = -eigenvals[0] + 1.0 
321              if self.print_flag >= 2: 
322                  print(self.print_prefix + "\tSafeguarding. lambda(l) = " + repr(self.lambda_l)) 
323          elif self.lambda_l <= 0.0: 
324              if self.print_flag >= 2: 
325                  print(self.print_prefix + "\tSafeguarding. lambda(l)=0") 
326              self.lambda_l = 0.0 
327   
328   
329 -    def update(self): 
330          """Update function. 
331   
332          Run the trust region update.  If this update decides to shift xk+1 to xk, then run the Newton update. 
333          """ 
334   
335          self.trust_region_update() 
336          if self.shift_flag: 
337              self.specific_update() 
338