minimise.dogleg

34 """Dogleg trust region algorithm. 35 36 Page 71 from 'Numerical Optimization' by Jorge Nocedal and Stephen J. Wright, 1999, 2nd ed. 37 The dogleg method is defined by the trajectory p(tau): 38 39 / tau . pU 0 <= tau <= 1, 40 p(tau) = < 41 \ pU + (tau - 1)(pB - pU), 1 <= tau <= 2. 42 43 where: 44 tau is in [0, 2] 45 pU is the unconstrained minimiser along the steepest descent direction. 46 pB is the full step. 47 48 pU is defined by the formula: 49 50 gT.g 51 pU = - ------ g 52 gT.B.g 53 54 and pB by the formula: 55 56 pB = - B^(-1).g 57 58 If the full step is within the trust region it is taken. Otherwise the point at which the 59 dogleg trajectory intersects the trust region is taken. This point can be found by solving the 60 scalar quadratic equation: 61 ||pU + (tau - 1)(pB - pU)||^2 = delta^2 62 """ 63 64 if print_flag: 65 if print_flag >= 2: 66 print print_prefix 67 print print_prefix 68 print print_prefix + "Dogleg minimisation" 69 print print_prefix + "~~~~~~~~~~~~~~~~~~~" 70 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) 71 if min.init_failure: 72 print print_prefix + "Initialisation of minimisation has failed." 73 return None 74 results = min.minimise() 75 return results

79 - 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):

80 """Class for Dogleg trust region minimisation specific functions. 81 82 Unless you know what you are doing, you should call the function 'dogleg' rather than using 83 this class. 84 """ 85 86 # Function arguments. 87 self.func = func 88 self.dfunc = dfunc 89 self.d2func = d2func 90 self.args = args 91 self.xk = x0 92 self.func_tol = func_tol 93 self.grad_tol = grad_tol 94 self.maxiter = maxiter 95 self.full_output = full_output 96 self.print_flag = print_flag 97 self.print_prefix = print_prefix 98 self.mach_acc = mach_acc 99 self.delta_max = delta_max 100 self.delta = delta0 101 self.eta = eta 102 103 # Initialisation failure flag. 104 self.init_failure = 0 105 106 # Setup the Hessian modification options and algorithm. 107 self.hessian_type_and_mod(min_options) 108 self.setup_hessian_mod() 109 110 # Initialise the function, gradient, and Hessian evaluation counters. 111 self.f_count = 0 112 self.g_count = 0 113 self.h_count = 0 114 115 # Initialise the warning string. 116 self.warning = None 117 118 # Constants. 119 self.n = len(self.xk) 120 self.I = identity(len(self.xk)) 121 122 # Set the convergence test function. 123 self.setup_conv_tests() 124 125 # Type specific functions. 126 if self.hessian_type and match('[Bb][Ff][Gg][Ss]', self.hessian_type): 127 self.setup_bfgs() 128 self.specific_update = self.update_bfgs 129 self.hessian_update = self.hessian_update_bfgs 130 self.d2fk = inverse(self.Hk) 131 elif self.hessian_type and match('[Nn]ewton', self.hessian_type): 132 self.setup_newton() 133 self.specific_update = self.update_newton 134 self.hessian_update = self.hessian_update_newton

135 136

137 - def dogleg(self):

138 """The dogleg algorithm.""" 139 140 # Calculate the full step and its norm. 141 try: 142 pB = -matrixmultiply(self.Hk, self.dfk) 143 except AttributeError: 144 # Backup the Hessian as the function self.get_pk may modify it. 145 d2fk_backup = 1.0 * self.d2fk 146 147 # The modified Newton step. 148 pB = self.get_pk() 149 150 # Restore the Hessian. 151 self.d2fk = d2fk_backup 152 norm_pB = sqrt(dot(pB, pB)) 153 154 # Test if the full step is within the trust region. 155 if norm_pB <= self.delta: 156 return pB 157 158 # Calculate pU. 159 curv = dot(self.dfk, dot(self.d2fk, self.dfk)) 160 pU = - dot(self.dfk, self.dfk) / curv * self.dfk 161 dot_pU = dot(pU, pU) 162 norm_pU = sqrt(dot_pU) 163 164 # Test if the step pU exits the trust region. 165 if norm_pU >= self.delta: 166 return self.delta * pU / norm_pU 167 168 # Find the solution to the scalar quadratic equation. 169 pB_pU = pB - pU 170 dot_pB_pU = dot(pB_pU, pB_pU) 171 dot_pU_pB_pU = dot(pU, pB_pU) 172 fact = dot_pU_pB_pU**2 - dot_pB_pU * (dot_pU - self.delta**2) 173 tau = (-dot_pU_pB_pU + sqrt(fact)) / dot_pB_pU 174 175 # Decide on which part of the trajectory to take. 176 return pU + tau * pB_pU

177 178

179 - def hessian_update_bfgs(self):

180 """BFGS Hessian update.""" 181 182 # BFGS matrix update. 183 sk = self.xk_new - self.xk 184 yk = self.dfk_new - self.dfk 185 if dot(yk, sk) == 0: 186 self.warning = "The BFGS matrix is indefinite. This should not occur." 187 rk = 1e99 188 else: 189 rk = 1.0 / dot(yk, sk) 190 191 if self.k == 0: 192 self.Hk = dot(yk, sk) / dot(yk, yk) * self.I 193 194 a = self.I - rk*outerproduct(sk, yk) 195 b = self.I - rk*outerproduct(yk, sk) 196 c = rk*outerproduct(sk, sk) 197 matrix = matrixmultiply(matrixmultiply(a, self.Hk), b) + c 198 self.Hk = matrix 199 200 # Calculate the Hessian. 201 self.d2fk = inverse(self.Hk)

202 203

204 - def hessian_update_newton(self):

205 """Empty function.""" 206 207 pass

208 209

210 - def new_param_func(self):

211 """Find the dogleg minimiser.""" 212 213 self.pk = self.dogleg() 214 self.xk_new = self.xk + self.pk 215 self.fk_new, self.f_count = apply(self.func, (self.xk_new,)+self.args), self.f_count + 1 216 self.dfk_new, self.g_count = apply(self.dfunc, (self.xk_new,)+self.args), self.g_count + 1 217 218 if self.print_flag >= 2: 219 print self.print_prefix + "Fin." 220 print self.print_prefix + " pk: " + `self.pk` 221 print self.print_prefix + " xk: " + `self.xk` 222 print self.print_prefix + " xk_new: " + `self.xk_new` 223 print self.print_prefix + " fk: " + `self.fk` 224 print self.print_prefix + " fk_new: " + `self.fk_new`

225 226

227 - def update(self):

228 """Update function. 229 230 Run the trust region update. If this update decides to shift xk+1 to xk, then run the BFGS 231 or Newton updates. 232 """ 233 234 self.trust_region_update() 235 236 if not self.shift_flag: 237 self.hessian_update() 238 else: 239 self.specific_update()

Source Code for Module minimise.dogleg