Source Code for Module minimise.exact_trust_region

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2003 Edward d'Auvergne                                        # 
  4  #                                                                             # 
  5  # This file is part of the program relax.                                     # 
  6  #                                                                             # 
  7  # relax is free software; you can redistribute it and/or modify               # 
  8  # it under the terms of the GNU General Public License as published by        # 
  9  # the Free Software Foundation; either version 2 of the License, or           # 
 10  # (at your option) any later version.                                         # 
 11  #                                                                             # 
 12  # relax is distributed in the hope that it will be useful,                    # 
 13  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 14  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 15  # GNU General Public License for more details.                                # 
 16  #                                                                             # 
 17  # You should have received a copy of the GNU General Public License           # 
 18  # along with relax; if not, write to the Free Software                        # 
 19  # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA   # 
 20  #                                                                             # 
 21  ############################################################################### 
 22   
 23   
 24  from LinearAlgebra import cholesky_decomposition, eigenvectors, inverse, solve_linear_equations 
 25  from Numeric import diagonal, dot, identity, matrixmultiply, outerproduct, sort, sqrt, transpose 
 26  from re import match 
 27   
 28  from bfgs import Bfgs 
 29  from newton import Newton 
 30  from base_classes import Hessian_mods, Min, Trust_region 
 31   
 32   
 33 -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=""): 
 34      """Exact trust region algorithm.""" 
 35   
 36      if print_flag: 
 37          if print_flag >= 2: 
 38              print print_prefix 
 39          print print_prefix 
 40          print print_prefix + "Exact trust region minimisation" 
 41          print print_prefix + "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" 
 42      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) 
 43      if min.init_failure: 
 44          print print_prefix + "Initialisation of minimisation has failed." 
 45          return None 
 46      results = min.minimise() 
 47      return results 
 48   
 49   
 50 -class Exact_trust_region(Hessian_mods, Trust_region, Min, Bfgs, Newton): 
 51 -    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): 
 52          """Class for Exact trust region minimisation specific functions. 
 53   
 54          Unless you know what you are doing, you should call the function 'exact_trust_region' rather 
 55          than using this class. 
 56          """ 
 57   
 58          # Function arguments. 
 59          self.func = func 
 60          self.dfunc = dfunc 
 61          self.d2func = d2func 
 62          self.args = args 
 63          self.xk = x0 
 64          self.func_tol = func_tol 
 65          self.grad_tol = grad_tol 
 66          self.maxiter = maxiter 
 67          self.full_output = full_output 
 68          self.print_flag = print_flag 
 69          self.print_prefix = print_prefix 
 70          self.mach_acc = mach_acc 
 71          self.lambda0 = lambda0 
 72          self.delta_max = delta_max 
 73          self.delta = delta0 
 74          self.eta = eta 
 75   
 76          # Initialisation failure flag. 
 77          self.init_failure = 0 
 78   
 79          # Setup the Hessian modification options and algorithm. 
 80          self.hessian_type_and_mod(min_options) 
 81          self.setup_hessian_mod() 
 82   
 83          # Initialise the function, gradient, and Hessian evaluation counters. 
 84          self.f_count = 0 
 85          self.g_count = 0 
 86          self.h_count = 0 
 87   
 88          # Initialise the warning string. 
 89          self.warning = None 
 90   
 91          # Constants. 
 92          self.n = len(self.xk) 
 93          self.I = identity(len(self.xk)) 
 94   
 95          # Set the convergence test function. 
 96          self.setup_conv_tests() 
 97   
 98          # Type specific functions. 
 99          if match('[Bb][Ff][Gg][Ss]', self.hessian_type): 
100              self.setup_bfgs() 
101              self.specific_update = self.update_bfgs 
102              self.d2fk = inverse(self.Hk) 
103              self.warning = "Incomplete code." 
104          elif match('[Nn]ewton', self.hessian_type): 
105              self.setup_newton() 
106              self.specific_update = self.update_newton 
107   
108   
109 -    def new_param_func(self): 
110          """Find the exact trust region solution. 
111   
112          Algorithm 4.4 from page 81 of 'Numerical Optimization' by Jorge Nocedal and Stephen J. 
113          Wright, 1999, 2nd ed. 
114   
115          This is only implemented for positive definite matrices. 
116          """ 
117   
118          # Matrix modification. 
119          pB, matrix = self.get_pk(return_matrix=1) 
120   
121          # The exact trust region algorithm. 
122          l = 0 
123          lambda_l = self.lambda0 
124   
125          while l < 3: 
126              # Calculate the matrix B + lambda(l).I 
127              B = matrix + lambda_l * self.I 
128   
129              # Factor B + lambda(l).I = RT.R 
130              R = cholesky_decomposition(B) 
131   
132              # Solve RT.R.pl = -g 
133              y = solve_linear_equations(R, self.dfk) 
134              y = -solve_linear_equations(transpose(R), y) 
135              dot_pl = dot(y, y) 
136   
137              # Solve RT.ql = pl 
138              y = solve_linear_equations(transpose(R), y) 
139   
140              # lambda(l+1) update. 
141              lambda_l = lambda_l + (dot_pl / dot(y, y)) * ((sqrt(dot_pl) - self.delta) / self.delta) 
142   
143              # Safeguard. 
144              lambda_l = max(0.0, lambda_l) 
145   
146              if self.print_flag >= 2: 
147                  print self.print_prefix + "\tl: " + `l` + ", lambda(l) fin: " + `lambda_l` 
148   
149              l = l + 1 
150   
151          # Calculate the step. 
152          R = cholesky_decomposition(matrix + lambda_l * self.I) 
153          y = solve_linear_equations(R, self.dfk) 
154          self.pk = -solve_linear_equations(transpose(R), y) 
155   
156          if self.print_flag >= 2: 
157              print self.print_prefix + "Step: " + `self.pk` 
158   
159          # Find the new parameter vector and function value at that point. 
160          self.xk_new = self.xk + self.pk 
161          self.fk_new, self.f_count = apply(self.func, (self.xk_new,)+self.args), self.f_count + 1 
162          self.dfk_new, self.g_count = apply(self.dfunc, (self.xk_new,)+self.args), self.g_count + 1 
163   
164   
165 -    def old_param_func(self): 
166          """Find the exact trust region solution. 
167   
168          More, J. J., and Sorensen D. C. 1983, Computing a trust region step.  SIAM J. Sci. Stat. 
169          Comput. 4, 553-572. 
170   
171          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 xrange(len(self.d2fk)): 
187              self.lS = max(self.lS, -self.d2fk[j, j]) 
188              sum = 0.0 
189              for i in xrange(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 = eigenvectors(self.d2fk) 
200              eigenvals = sort(eigen[0]) 
201              for i in xrange(len(self.d2fk)): 
202                  print self.print_prefix + "\tB[" + `i` + ", " + `i` + "] = " + `self.d2fk[i, i]` 
203              print self.print_prefix + "\tEigenvalues: " + `eigenvals` 
204              print self.print_prefix + "\t||g||/delta: " + `a` 
205              print self.print_prefix + "\t||B||1: " + `b` 
206              print self.print_prefix + "\tl:  " + `self.l` 
207              print self.print_prefix + "\tlL: " + `self.lL` 
208              print self.print_prefix + "\tlU: " + `self.lU` 
209              print self.print_prefix + "\tlS: " + `self.lS` 
210   
211          # Iterative loop. 
212          return 
213          while 1: 
214              # Safeguard lambda. 
215              if self.print_flag >= 2: 
216                  print self.print_prefix + "\n< Iteration " + `iter` + " >" 
217                  print self.print_prefix + "Safeguarding lambda." 
218                  print self.print_prefix + "\tInit l: " + `self.l` 
219                  print self.print_prefix + "\tlL: " + `self.lL` 
220                  print self.print_prefix + "\tlU: " + `self.lU` 
221                  print self.print_prefix + "\tlS: " + `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: " + `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: " + `matrix` 
237                  eigen = eigenvectors(matrix) 
238                  eigenvals = sort(eigen[0]) 
239                  print self.print_prefix + "\tEigenvalues: " + `eigenvals` 
240              try: 
241                  func = cholesky_decomposition 
242                  R = func(matrix) 
243                  if self.print_flag >= 2: 
244                      print self.print_prefix + "\tCholesky matrix R: " + `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: " + `pos_def` 
251   
252              if pos_def: 
253                  # Solve p = -inverse(RT.R).g 
254                  p = -matrixmultiply(inverse(matrix), self.dfk) 
255                  if self.print_flag >= 2: 
256                      print self.print_prefix + "Solve p = -inverse(RT.R).g" 
257                      print self.print_prefix + "\tp: " + `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: " + `z` 
274                          print self.print_prefix + "\ttau: " + `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 = inverse(RT).p 
281                  q = matrixmultiply(inverse(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: " + `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 = apply(self.func, (self.xk_new,)+self.args), self.f_count + 1 
313          self.dfk_new, self.g_count = apply(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) = " + `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 
333          Newton update. 
334          """ 
335   
336          self.trust_region_update() 
337          if self.shift_flag: 
338              self.specific_update() 
339