minimise.levenberg

31 """Levenberg-Marquardt minimimisation. 32 33 Function options 34 ~~~~~~~~~~~~~~~~ 35 36 chi2_func: User supplied chi-squared function which is run with the function parameters and 37 args as options. 38 39 dchi2_func: User supplied chi-squared gradient function which is run with the function 40 parameters and args as options. 41 42 dfunc: User supplied function which should return a vector of partial derivatives of the 43 function which back calculates values for the chi-squared function. 44 45 params: The initial function parameter values. 46 47 errors: The experimental errors. 48 49 args: A tuple containing the arguments to send to chi2_func and dchi2_func. 50 51 maxiter: The maximum number of iterations. 52 53 full_output: A flag specifying what should be returned. 54 55 56 Output 57 ~~~~~~ 58 59 If full_output = 0, the parameter values and chi-squared value are returned as a tuple. 60 61 If full_output = 1, the parameter values, chi-squared value, number of iterations, and the 62 warning flag are returned as a tuple. 63 """ 64 65 if print_flag: 66 if print_flag >= 2: 67 print print_prefix 68 print print_prefix 69 print print_prefix + "Levenberg-Marquardt minimisation" 70 print print_prefix + "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" 71 min = Levenberg_marquardt(chi2_func, dchi2_func, dfunc, errors, args, x0, func_tol, grad_tol, maxiter, full_output, print_flag, print_prefix) 72 results = min.minimise() 73 return results

77 - def __init__(self, chi2_func, dchi2_func, dfunc, errors, args, x0, func_tol, grad_tol, maxiter, full_output, print_flag, print_prefix):

78 """Class for Levenberg-Marquardt minimisation specific functions. 79 80 Unless you know what you are doing, you should call the function 81 'levenberg_marquardt' rather than using this class. 82 """ 83 84 # Function arguments. 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 # Initialise the function, gradient, and Hessian evaluation counters. 99 self.f_count = 0 100 self.g_count = 0 101 self.h_count = 0 102 103 # Initialise the warning string. 104 self.warning = None 105 106 # Set the convergence test function. 107 self.setup_conv_tests() 108 self.orig_conv_test = self.conv_test 109 self.conv_test = self.test_mod 110 111 # The initial chi-squared value, chi-squared gradient vector, and derivative function matrix. 112 self.fk, self.f_count = apply(self.chi2_func, (self.xk,)+self.args), self.f_count + 1 113 self.dfk, self.g_count = -0.5 * apply(self.dchi2_func, (self.xk,)+self.args), self.g_count + 1 114 self.df = self.dfunc() 115 self.dfk_new = None 116 117 # Initial value of lambda (the Levenberg-Marquardt fudge factor). 118 self.l = 0.001 119 120 # N. 121 self.n = len(self.xk)

122 123

124 - def create_lm_matrix(self):

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. 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 # Create the Levenberg-Marquardt matrix with elements equal to zero. 147 self.lm_matrix = zeros((self.n, self.n), Float64) 148 149 # Calculate the inverse of the variance to minimise calculations. 150 i_variance = 1.0 / self.errors**2 151 152 # Loop over the error points from i=1 to n. 153 for i in xrange(len(self.errors)): 154 # Loop over all function parameters. 155 for param_j in xrange(self.n): 156 # Loop over the function parameters from the first to 'param_j' to create the Levenberg-Marquardt matrix. 157 for param_k in xrange(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

167 - def new_param_func(self):

168 """Find the new parameter vector self.xk_new.""" 169 170 # Create the Levenberg-Marquardt matrix. 171 self.create_lm_matrix() 172 173 # Solve the Levenberg-Marquardt equation to get the vector of function parameter changes. 174 #print "\nself.lm_matrix:\n" + `self.lm_matrix` 175 #print "\nself.dfk:\n" + `self.dfk` 176 self.pk = solve_linear_equations(self.lm_matrix, self.dfk) 177 178 # Find the new parameter vector and function value at that point. 179 self.xk_new = self.xk + self.pk 180 self.fk_new, self.f_count = apply(self.chi2_func, (self.xk_new,)+self.args), self.f_count + 1 181 182 if self.fk_new <= self.fk: 183 # Move to xk+1 and shrink lambda. 184 self.dfk_new, self.g_count = -0.5 * apply(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 # Don't move and increase lambda. 190 if self.l <= 1e99: 191 self.l = self.l * 10.0 192 self.move_flag = 0 193 194 # Print out. 195 if self.print_flag >= 2: 196 print self.print_prefix + "xk_new: " + `self.xk_new` 197 print self.print_prefix + "lm_matrix: " + `self.lm_matrix` 198 print self.print_prefix + "df: " + `self.df` 199 print self.print_prefix + "l: " + `self.l` 200 print self.print_prefix + "fk+1: " + `self.fk_new` 201 print self.print_prefix + "fk: " + `self.fk` 202 print self.print_prefix + "move_flag: " + `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 209 movement during an iteration due to an uphill step being encountered. 210 """ 211 212 if self.move_flag: 213 if self.orig_conv_test(fk_new, fk, dfk_new): 214 return 1

215 216

217 - def update(self):

218 """Update function 219 220 Update the chi-squared value, chi-squared gradient vector, and derivative function matrix. 221 """ 222 223 if self.move_flag: 224 self.xk = self.xk_new * 1.0 225 self.fk = self.fk_new 226 self.dfk = self.dfk_new * 1.0 227 self.df = self.dfunc()

Source Code for Module minimise.levenberg_marquardt