1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24 """Newton (or Newton-Raphson - NR) optimization.
25
26 This file is part of the U{minfx optimisation library<https://sourceforge.net/projects/minfx>}.
27 """
28
29
30 from math import acos
31 from numpy import dot, float64, identity, sqrt, zeros
32 from re import match
33
34
35 from minfx.base_classes import Hessian_mods, Line_search, Min
36
37
38 -def newton(func=None, dfunc=None, d2func=None, args=(), x0=None, min_options=(), func_tol=1e-25, grad_tol=None, maxiter=1e6, a0=1.0, mu=0.0001, eta=0.9, mach_acc=1e-16, full_output=0, print_flag=0, print_prefix=""):
39 """Newton minimisation."""
40
41 if print_flag:
42 if print_flag >= 2:
43 print(print_prefix)
44 print(print_prefix)
45 print(print_prefix + "Newton minimisation")
46 print(print_prefix + "~~~~~~~~~~~~~~~~~~~")
47 min = Newton(func, dfunc, d2func, args, x0, min_options, func_tol, grad_tol, maxiter, a0, mu, eta, mach_acc, full_output, print_flag, print_prefix)
48 if min.init_failure:
49 print(print_prefix + "Initialisation of minimisation has failed.")
50 return None
51 results = min.minimise()
52 return results
53
54
55 -class Newton(Hessian_mods, Line_search, Min):
56 - def __init__(self, func, dfunc, d2func, args, x0, min_options, func_tol, grad_tol, maxiter, a0, mu, eta, mach_acc, full_output, print_flag, print_prefix):
57 """Class for Newton minimisation specific functions.
58
59 Unless you know what you are doing, you should call the function 'newton' rather than using this class.
60 """
61
62
63 self.func = func
64 self.dfunc = dfunc
65 self.d2func = d2func
66 self.args = args
67 self.xk = x0
68 self.func_tol = func_tol
69 self.grad_tol = grad_tol
70 self.maxiter = maxiter
71 self.mach_acc = mach_acc
72 self.full_output = full_output
73 self.print_flag = print_flag
74 self.print_prefix = print_prefix
75
76
77 self.a0 = a0
78
79
80 self.mu = mu
81 self.eta = eta
82
83
84 self.init_failure = 0
85
86
87 self.line_search_algor = None
88 self.hessian_mod = None
89
90
91 if not isinstance(min_options, tuple):
92 print(self.print_prefix + "The minimisation options " + repr(min_options) + " is not a tuple.")
93 self.init_failure = 1; return
94
95
96 if len(min_options) > 2:
97 print(self.print_prefix + "A maximum of two minimisation options is allowed (the line search algorithm and the Hessian modification).")
98 self.init_failure = 1; return
99
100
101 none_option = None
102 for opt in min_options:
103 if match('[Nn]one', opt):
104 none_option = opt
105 continue
106 if self.line_search_algor == None and self.valid_line_search(opt):
107 self.line_search_algor = opt
108 elif self.hessian_mod == None and self.valid_hessian_mod(opt):
109 self.hessian_mod = opt
110 else:
111 if self.line_search_algor:
112 print(self.print_prefix + "The minimisation option " + repr(opt) + " from " + repr(min_options) + " is not a valid Hessian modification.")
113 elif self.hessian_mod:
114 print(self.print_prefix + "The minimisation option " + repr(opt) + " from " + repr(min_options) + " is not a valid line search algorithm.")
115 else:
116 print(self.print_prefix + "The minimisation option " + repr(opt) + " from " + repr(min_options) + " is neither a valid line search algorithm or Hessian modification.")
117 self.init_failure = 1; return
118 if none_option and not self.hessian_mod:
119 if self.valid_hessian_mod(none_option):
120 self.hessian_mod = none_option
121 elif none_option and not self.line_search_algor:
122 if self.valid_line_search(none_option):
123 self.line_search_algor = none_option
124 elif none_option:
125 print(self.print_prefix + "Invalid combination of options " + repr(min_options) + ".")
126 self.init_failure = 1; return
127
128
129 if self.line_search_algor == None:
130 self.line_search_algor = 'Back'
131
132
133 if self.hessian_mod == None:
134 self.hessian_mod = 'GMW'
135
136
137 self.setup_line_search()
138 self.setup_hessian_mod()
139
140
141 self.f_count = 0
142 self.g_count = 0
143 self.h_count = 0
144
145
146 self.warning = None
147
148
149 self.n = len(self.xk)
150 self.I = identity(len(self.xk))
151
152
153 self.setup_conv_tests()
154
155
156 self.setup_newton()
157
158
159 self.update = self.update_newton
160
161
163 """The new parameter function.
164
165 Find the search direction, do a line search, and get xk+1 and fk+1.
166 """
167
168
169 self.pk = self.get_pk()
170
171
172 self.line_search()
173
174
175 self.xk_new = self.xk + self.alpha * self.pk
176 self.fk_new, self.f_count = self.func(*(self.xk_new,)+self.args), self.f_count + 1
177 self.dfk_new, self.g_count = self.dfunc(*(self.xk_new,)+self.args), self.g_count + 1
178
179
180 if self.print_flag >= 2:
181 print("\n" + self.print_prefix + "New param function.")
182 print(self.print_prefix + "pk: " + repr(self.pk))
183 print(self.print_prefix + "alpha: " + repr(self.alpha))
184 print(self.print_prefix + "xk: " + repr(self.xk))
185 print(self.print_prefix + "xk+1: " + repr(self.xk_new))
186 print(self.print_prefix + "fk: " + repr(self.fk))
187 print(self.print_prefix + "fk+1: " + repr(self.fk_new))
188 print(self.print_prefix + "dfk: " + repr(self.dfk))
189 print(self.print_prefix + "dfk+1: " + repr(self.dfk_new))
190 print(self.print_prefix + "d2fk:\n" + repr(self.d2fk))
191
192
193
194
195 print(self.print_prefix + "Angle to the unit vector pointing along the first dimension.")
196 unit_vect = zeros(self.n, float64)
197 unit_vect[0] = 1.0
198 dfk_angle = acos(dot(self.dfk, unit_vect) / sqrt(dot(self.dfk, self.dfk)))
199 pk_angle = acos(dot(self.pk, unit_vect) / sqrt(dot(self.pk, self.pk)))
200 print(self.print_prefix + "steepest descent: " + repr(dfk_angle))
201 print(self.print_prefix + "Newton step: " + repr(pk_angle))
202
203
205 """Setup function.
206
207 The initial Newton function value, gradient vector, and Hessian matrix are calculated.
208 """
209
210 self.fk, self.f_count = self.func(*(self.xk,)+self.args), self.f_count + 1
211 self.dfk, self.g_count = self.dfunc(*(self.xk,)+self.args), self.g_count + 1
212 self.d2fk, self.h_count = self.d2func(*(self.xk,)+self.args), self.h_count + 1
213
214
216 """Function to update the function value, gradient vector, and Hessian matrix."""
217
218 self.xk = self.xk_new * 1.0
219 self.fk = self.fk_new
220 self.dfk = self.dfk_new * 1.0
221 self.d2fk, self.h_count = self.d2func(*(self.xk,)+self.args), self.h_count + 1
222