1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
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
30 from numpy import dot, float64, identity, outer, sort, sqrt
31 from numpy.linalg import inv
32 from re import match
33
34
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
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
109 self.init_failure = 0
110
111
112 self.hessian_type_and_mod(min_options)
113 self.setup_hessian_mod()
114
115
116 self.f_count = 0
117 self.g_count = 0
118 self.h_count = 0
119
120
121 self.warning = None
122
123
124 self.n = len(self.xk)
125 self.I = identity(len(self.xk))
126
127
128 self.setup_conv_tests()
129
130
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
143 """The dogleg algorithm."""
144
145
146 try:
147 pB = -dot(self.Hk, self.dfk)
148 except AttributeError:
149
150 d2fk_backup = 1.0 * self.d2fk
151
152
153 pB = self.get_pk()
154
155
156 self.d2fk = d2fk_backup
157 norm_pB = sqrt(dot(pB, pB))
158
159
160 if norm_pB <= self.delta:
161 return pB
162
163
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
170 if norm_pU >= self.delta:
171 return self.delta * pU / norm_pU
172
173
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
181 return pU + tau * pB_pU
182
183
185 """BFGS Hessian update."""
186
187
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
206 self.d2fk = inv(self.Hk)
207
208
210 """Empty function."""
211
212 pass
213
214
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
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()