minfx.levenberg

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

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 # 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 = 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 # 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, 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 # 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 range(len(self.errors)): 154 # Loop over all function parameters. 155 for param_j in range(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 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

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" + repr(self.lm_matrix)) 175 #print("\nself.dfk:\n" + repr(self.dfk)) 176 self.pk = solve(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 = 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 * 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: " + 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

216 - def update(self):

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()

Source Code for Module minfx.levenberg_marquardt