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