1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24 """Levenberg-Marquardt 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 float64, zeros
31 from numpy.linalg import solve
32
33
34 from minfx.base_classes import Min
35
36
37 -def levenberg_marquardt(chi2_func=None, dchi2_func=None, dfunc=None, errors=None, args=(), x0=None, func_tol=1e-25, grad_tol=None, maxiter=1e6, print_flag=0, print_prefix="", full_output=False):
38 """Levenberg-Marquardt minimisation.
39
40 @keyword chi2_func: User supplied chi-squared function which is run with the function parameters and args as options.
41 @keyword chi2_func: func
42 @keyword dchi2_func: User supplied chi-squared gradient function which is run with the function parameters and args as options.
43 @keyword dchi2_func: func
44 @keyword dfunc: User supplied function which should return a vector of partial derivatives of the function which back calculates values for the chi-squared function.
45 @keyword dfunc: func
46 @keyword errors: The experimental errors.
47 @keyword errors: numpy rank-1 array
48 @keyword args: A tuple containing the arguments to send to chi2_func and dchi2_func.
49 @keyword args: tuple
50 @keyword x0: The initial function parameter values.
51 @keyword x0: numpy rank-1 array
52 @keyword func_tol: The function tolerance value. Once the function value between iterations decreases below this value, minimisation is terminated.
53 @type func_tol: float
54 @keyword grad_tol: The gradient tolerance value.
55 @type grad_tol: float
56 @keyword maxiter: The maximum number of iterations.
57 @keyword maxiter: int
58 @keyword print_flag: A flag specifying how much information should be printed to standard output during minimisation. 0 means no output, 1 means minimal output, and values above 1 increase the amount of output printed.
59 @type print_flag: int
60 @keyword print_prefix: The text to add out to the front of all printouts.
61 @type print_prefix: str
62 @keyword full_output: A flag specifying what should be returned. If full_output is False, the parameter values and chi-squared value are returned as a tuple. If full_output is True, the parameter values, chi-squared value, number of iterations, and the warning flag are returned as a tuple.
63 @keyword full_output: bool
64 """
65
66 if print_flag:
67 if print_flag >= 2:
68 print(print_prefix)
69 print(print_prefix)
70 print(print_prefix + "Levenberg-Marquardt minimisation")
71 print(print_prefix + "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")
72 min = Levenberg_marquardt(chi2_func, dchi2_func, dfunc, errors, args, x0, func_tol, grad_tol, maxiter, full_output, print_flag, print_prefix)
73 results = min.minimise()
74 return results
75
76
78 - def __init__(self, chi2_func, dchi2_func, dfunc, errors, args, x0, func_tol, grad_tol, maxiter, full_output, print_flag, print_prefix):
79 """Class for Levenberg-Marquardt minimisation specific functions.
80
81 Unless you know what you are doing, you should call the function 'levenberg_marquardt' rather than using this class.
82 """
83
84
85 self.chi2_func = chi2_func
86 self.dchi2_func = dchi2_func
87 self.dfunc = dfunc
88 self.errors = errors
89 self.args = args
90 self.xk = x0
91 self.func_tol = func_tol
92 self.grad_tol = grad_tol
93 self.maxiter = maxiter
94 self.full_output = full_output
95 self.print_flag = print_flag
96 self.print_prefix = print_prefix
97
98
99 self.f_count = 0
100 self.g_count = 0
101 self.h_count = 0
102
103
104 self.warning = None
105
106
107 self.setup_conv_tests()
108 self.orig_conv_test = self.conv_test
109 self.conv_test = self.test_mod
110
111
112 self.fk, self.f_count = self.chi2_func(*(self.xk,)+self.args), self.f_count + 1
113 self.dfk, self.g_count = -0.5 * self.dchi2_func(*(self.xk,)+self.args), self.g_count + 1
114 self.df = self.dfunc()
115 self.dfk_new = None
116
117
118 self.l = 0.001
119
120
121 self.n = len(self.xk)
122
123
125 """Function to create the Levenberg-Marquardt matrix.
126
127 The matrix is::
128
129 _n_
130 \ / 1 d y(xi) d y(xi) \
131 LM_matrix_jk = > | ---------- . ------- . ------- . (1 + lambda) |
132 /__ \ sigma_i**2 dj dk /
133 i=1
134
135 where j == k is one of the function parameters, and::
136
137 _n_
138 \ / 1 d y(xi) d y(xi) \
139 LM_matrix_jk = > | ---------- . ------- . ------- |
140 /__ \ sigma_i**2 dj dk /
141 i=1
142
143 where j != k are function parameters.
144 """
145
146
147 self.lm_matrix = zeros((self.n, self.n), float64)
148
149
150 i_variance = 1.0 / self.errors**2
151
152
153 for i in range(len(self.errors)):
154
155 for param_j in range(self.n):
156
157 for param_k in range(param_j + 1):
158 if param_j == param_k:
159 matrix_jk = i_variance[i] * self.df[i, param_j] * self.df[i, param_k] * (1.0 + self.l)
160 self.lm_matrix[param_j, param_k] = self.lm_matrix[param_j, param_k] + matrix_jk
161 else:
162 matrix_jk = i_variance[i] * self.df[i, param_j] * self.df[i, param_k]
163 self.lm_matrix[param_j, param_k] = self.lm_matrix[param_j, param_k] + matrix_jk
164 self.lm_matrix[param_k, param_j] = self.lm_matrix[param_k, param_j] + matrix_jk
165
166
168 """Find the new parameter vector self.xk_new."""
169
170
171 self.create_lm_matrix()
172
173
174
175
176 self.pk = solve(self.lm_matrix, self.dfk)
177
178
179 self.xk_new = self.xk + self.pk
180 self.fk_new, self.f_count = self.chi2_func(*(self.xk_new,)+self.args), self.f_count + 1
181
182 if self.fk_new <= self.fk:
183
184 self.dfk_new, self.g_count = -0.5 * self.dchi2_func(*(self.xk_new,)+self.args), self.g_count + 1
185 if self.l >= 1e-99:
186 self.l = self.l * 0.1
187 self.move_flag = 1
188 else:
189
190 if self.l <= 1e99:
191 self.l = self.l * 10.0
192 self.move_flag = 0
193
194
195 if self.print_flag >= 2:
196 print(self.print_prefix + "xk_new: " + repr(self.xk_new))
197 print(self.print_prefix + "lm_matrix: " + repr(self.lm_matrix))
198 print(self.print_prefix + "df: " + repr(self.df))
199 print(self.print_prefix + "l: " + repr(self.l))
200 print(self.print_prefix + "fk+1: " + repr(self.fk_new))
201 print(self.print_prefix + "fk: " + repr(self.fk))
202 print(self.print_prefix + "move_flag: " + repr(self.move_flag))
203
204
205 - def test_mod(self, fk_new, fk, dfk_new):
206 """Modified convergence test.
207
208 This is needed to prevent the Levenberg-Marquardt minimiser from terminating if there is no movement during an iteration due to an uphill step being encountered.
209 """
210
211 if self.move_flag:
212 if self.orig_conv_test(fk_new, fk, dfk_new):
213 return 1
214
215
217 """Update function
218
219 Update the chi-squared value, chi-squared gradient vector, and derivative function matrix.
220 """
221
222 if self.move_flag:
223 self.xk = self.xk_new * 1.0
224 self.fk = self.fk_new
225 self.dfk = self.dfk_new * 1.0
226 self.df = self.dfunc()
227