Package minfx :: Module dogleg
[hide private]
[frames] | no frames]

Source Code for Module minfx.dogleg

  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  """Dogleg trust region 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 numpy import dot, float64, identity, outer, sort, sqrt 
 31  from numpy.linalg import inv 
 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 dogleg(func=None, dfunc=None, d2func=None, args=(), x0=None, min_options=(), func_tol=1e-25, grad_tol=None, maxiter=1e6, delta_max=1e10, delta0=1e5, eta=0.0001, mach_acc=1e-16, full_output=0, print_flag=0, print_prefix=""):
41 """Dogleg trust region algorithm. 42 43 Page 71 from 'Numerical Optimization' by Jorge Nocedal and Stephen J. Wright, 1999, 2nd ed. The dogleg method is defined by the trajectory p(tau):: 44 45 / tau . pU 0 <= tau <= 1, 46 p(tau) = < 47 \ pU + (tau - 1)(pB - pU), 1 <= tau <= 2. 48 49 where: 50 51 - tau is in [0, 2] 52 - pU is the unconstrained minimiser along the steepest descent direction. 53 - pB is the full step. 54 55 pU is defined by the formula:: 56 57 gT.g 58 pU = - ------ g 59 gT.B.g 60 61 and pB by the formula:: 62 63 pB = - B^(-1).g 64 65 If the full step is within the trust region it is taken. Otherwise the point at which the dogleg trajectory intersects the trust region is taken. This point can be found by solving the scalar quadratic equation:: 66 67 ||pU + (tau - 1)(pB - pU)||^2 = delta^2 68 """ 69 70 if print_flag: 71 if print_flag >= 2: 72 print(print_prefix) 73 print(print_prefix) 74 print(print_prefix + "Dogleg minimisation") 75 print(print_prefix + "~~~~~~~~~~~~~~~~~~~") 76 min = Dogleg(func, dfunc, d2func, args, x0, min_options, func_tol, grad_tol, maxiter, delta_max, delta0, eta, mach_acc, full_output, print_flag, print_prefix) 77 if min.init_failure: 78 print(print_prefix + "Initialisation of minimisation has failed.") 79 return None 80 results = min.minimise() 81 return results
82 83
84 -class Dogleg(Trust_region, Bfgs, Newton):
85 - def __init__(self, func, dfunc, d2func, args, x0, min_options, func_tol, grad_tol, maxiter, delta_max, delta0, eta, mach_acc, full_output, print_flag, print_prefix):
86 """Class for Dogleg trust region minimisation specific functions. 87 88 Unless you know what you are doing, you should call the function 'dogleg' rather than using this class. 89 """ 90 91 # Function arguments. 92 self.func = func 93 self.dfunc = dfunc 94 self.d2func = d2func 95 self.args = args 96 self.xk = x0 97 self.func_tol = func_tol 98 self.grad_tol = grad_tol 99 self.maxiter = maxiter 100 self.full_output = full_output 101 self.print_flag = print_flag 102 self.print_prefix = print_prefix 103 self.mach_acc = mach_acc 104 self.delta_max = delta_max 105 self.delta = delta0 106 self.eta = eta 107 108 # Initialisation failure flag. 109 self.init_failure = 0 110 111 # Setup the Hessian modification options and algorithm. 112 self.hessian_type_and_mod(min_options) 113 self.setup_hessian_mod() 114 115 # Initialise the function, gradient, and Hessian evaluation counters. 116 self.f_count = 0 117 self.g_count = 0 118 self.h_count = 0 119 120 # Initialise the warning string. 121 self.warning = None 122 123 # Constants. 124 self.n = len(self.xk) 125 self.I = identity(len(self.xk)) 126 127 # Set the convergence test function. 128 self.setup_conv_tests() 129 130 # Type specific functions. 131 if self.hessian_type and match('[Bb][Ff][Gg][Ss]', self.hessian_type): 132 self.setup_bfgs() 133 self.specific_update = self.update_bfgs 134 self.hessian_update = self.hessian_update_bfgs 135 self.d2fk = inv(self.Hk) 136 elif self.hessian_type and match('[Nn]ewton', self.hessian_type): 137 self.setup_newton() 138 self.specific_update = self.update_newton 139 self.hessian_update = self.hessian_update_newton
140 141
142 - def dogleg(self):
143 """The dogleg algorithm.""" 144 145 # Calculate the full step and its norm. 146 try: 147 pB = -dot(self.Hk, self.dfk) 148 except AttributeError: 149 # Backup the Hessian as the function self.get_pk may modify it. 150 d2fk_backup = 1.0 * self.d2fk 151 152 # The modified Newton step. 153 pB = self.get_pk() 154 155 # Restore the Hessian. 156 self.d2fk = d2fk_backup 157 norm_pB = sqrt(dot(pB, pB)) 158 159 # Test if the full step is within the trust region. 160 if norm_pB <= self.delta: 161 return pB 162 163 # Calculate pU. 164 curv = dot(self.dfk, dot(self.d2fk, self.dfk)) 165 pU = - dot(self.dfk, self.dfk) / curv * self.dfk 166 dot_pU = dot(pU, pU) 167 norm_pU = sqrt(dot_pU) 168 169 # Test if the step pU exits the trust region. 170 if norm_pU >= self.delta: 171 return self.delta * pU / norm_pU 172 173 # Find the solution to the scalar quadratic equation. 174 pB_pU = pB - pU 175 dot_pB_pU = dot(pB_pU, pB_pU) 176 dot_pU_pB_pU = dot(pU, pB_pU) 177 fact = dot_pU_pB_pU**2 - dot_pB_pU * (dot_pU - self.delta**2) 178 tau = (-dot_pU_pB_pU + sqrt(fact)) / dot_pB_pU 179 180 # Decide on which part of the trajectory to take. 181 return pU + tau * pB_pU
182 183
184 - def hessian_update_bfgs(self):
185 """BFGS Hessian update.""" 186 187 # BFGS matrix update. 188 sk = self.xk_new - self.xk 189 yk = self.dfk_new - self.dfk 190 if dot(yk, sk) == 0: 191 self.warning = "The BFGS matrix is indefinite. This should not occur." 192 rk = 1e99 193 else: 194 rk = 1.0 / dot(yk, sk) 195 196 if self.k == 0: 197 self.Hk = dot(yk, sk) / dot(yk, yk) * self.I 198 199 a = self.I - rk*outer(sk, yk) 200 b = self.I - rk*outer(yk, sk) 201 c = rk*outer(sk, sk) 202 matrix = dot(dot(a, self.Hk), b) + c 203 self.Hk = matrix 204 205 # Calculate the Hessian. 206 self.d2fk = inv(self.Hk)
207 208
209 - def hessian_update_newton(self):
210 """Empty function.""" 211 212 pass
213 214
215 - def new_param_func(self):
216 """Find the dogleg minimiser.""" 217 218 self.pk = self.dogleg() 219 self.xk_new = self.xk + self.pk 220 self.fk_new, self.f_count = self.func(*(self.xk_new,)+self.args), self.f_count + 1 221 self.dfk_new, self.g_count = self.dfunc(*(self.xk_new,)+self.args), self.g_count + 1 222 223 if self.print_flag >= 2: 224 print(self.print_prefix + "Fin.") 225 print(self.print_prefix + " pk: " + repr(self.pk)) 226 print(self.print_prefix + " xk: " + repr(self.xk)) 227 print(self.print_prefix + " xk_new: " + repr(self.xk_new)) 228 print(self.print_prefix + " fk: " + repr(self.fk)) 229 print(self.print_prefix + " fk_new: " + repr(self.fk_new))
230 231
232 - def update(self):
233 """Update function. 234 235 Run the trust region update. If this update decides to shift xk+1 to xk, then run the BFGS or Newton updates. 236 """ 237 238 self.trust_region_update() 239 240 if not self.shift_flag: 241 self.hessian_update() 242 else: 243 self.specific_update()
244