minfx.ncg

38 """Line search Newton conjugate gradient algorithm. 39 40 Page 140 from 'Numerical Optimization' by Jorge Nocedal and Stephen J. Wright, 1999, 2nd ed. The algorithm is: 41 42 - Given initial point x0. 43 - while 1: 44 - Compute a search direction pk by applying the CG method to Hk.pk = -gk, starting from x0 = 0. Terminate when ||rk|| <= min(0.5,sqrt(||gk||)), or if negative curvature is encountered. 45 - Set xk+1 = xk + ak.pk, where ak satisfies the Wolfe, Goldstein, or Armijo backtracking conditions. 46 """ 47 48 if print_flag: 49 if print_flag >= 2: 50 print(print_prefix) 51 print(print_prefix) 52 print(print_prefix + "Newton Conjugate Gradient minimisation") 53 print(print_prefix + "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~") 54 min = Ncg(func, dfunc, d2func, args, x0, min_options, func_tol, grad_tol, maxiter, a0, mu, eta, full_output, print_flag, print_prefix) 55 if min.init_failure: 56 print(print_prefix + "Initialisation of minimisation has failed.") 57 return None 58 results = min.minimise() 59 return results

63 - def __init__(self, func, dfunc, d2func, args, x0, min_options, func_tol, grad_tol, maxiter, a0, mu, eta, full_output, print_flag, print_prefix):

64 """Class for newton conjugate gradient minimisation specific functions. 65 66 Unless you know what you are doing, you should call the function 'ncg' rather than using this class. 67 """ 68 69 # Function arguments. 70 self.func = func 71 self.dfunc = dfunc 72 self.d2func = d2func 73 self.args = args 74 self.xk = x0 75 self.func_tol = func_tol 76 self.grad_tol = grad_tol 77 self.maxiter = maxiter 78 self.full_output = full_output 79 self.print_flag = print_flag 80 self.print_prefix = print_prefix 81 82 # Set a0. 83 self.a0 = a0 84 85 # Line search constants for the Wolfe conditions. 86 self.mu = mu 87 self.eta = eta 88 89 # Initialisation failure flag. 90 self.init_failure = 0 91 92 # Setup the line search options and algorithm. 93 self.line_search_options(min_options) 94 self.setup_line_search() 95 96 # Initialise the function, gradient, and Hessian evaluation counters. 97 self.f_count = 0 98 self.g_count = 0 99 self.h_count = 0 100 101 # Initialise the warning string. 102 self.warning = None 103 104 # Set the convergence test function. 105 self.setup_conv_tests() 106 107 # Calculate the initial Newton function value, gradient vector, and Hessian matrix. 108 self.fk, self.f_count = self.func(*(self.xk,)+self.args), self.f_count + 1 109 self.dfk, self.g_count = self.dfunc(*(self.xk,)+self.args), self.g_count + 1 110 self.d2fk, self.h_count = self.d2func(*(self.xk,)+self.args), self.h_count + 1

111 112

113 - def get_pk(self):

114 """The CG algorithm.""" 115 116 # Initial values at i = 0. 117 xi = zeros(len(self.xk), float64) 118 ri = self.dfk * 1.0 119 pi = -ri 120 dot_ri = dot(ri, ri) 121 len_g = sqrt(dot_ri) 122 residual_test = min(0.5 * len_g, dot_ri) 123 #residual_test = min(0.5, sqrt(len_g)) * len_g 124 125 # Debugging. 126 if self.print_flag >= 2: 127 print(self.print_prefix + "Initial data:") 128 print(self.print_prefix + "\tx0: " + repr(xi)) 129 print(self.print_prefix + "\tr0: " + repr(ri)) 130 print(self.print_prefix + "\tp0: " + repr(pi)) 131 132 i = 0 133 while True: 134 # Matrix product and curvature. 135 Api = dot(self.d2fk, pi) 136 curv = dot(pi, Api) 137 138 # Negative curvature test. 139 if curv <= 0.0: 140 if i == 0: 141 if dot_ri == 0.0: 142 return xi 143 else: 144 ai = dot_ri / curv 145 return xi + ai*pi 146 else: 147 return xi 148 if sqrt(dot_ri) <= residual_test: 149 return xi 150 151 # Conjugate gradient part 152 ai = dot_ri / curv 153 xi_new = xi + ai*pi 154 ri_new = ri + ai*Api 155 dot_ri_new = dot(ri_new, ri_new) 156 bi_new = dot_ri_new / dot_ri 157 pi_new = -ri_new + bi_new*pi 158 159 # Debugging. 160 if self.print_flag >= 2: 161 print("") 162 print(self.print_prefix + "Iteration i = " + repr(i)) 163 print(self.print_prefix + "Api: " + repr(Api)) 164 print(self.print_prefix + "Curv: " + repr(curv)) 165 print(self.print_prefix + "ai: " + repr(ai)) 166 print(self.print_prefix + "xi+1: " + repr(xi_new)) 167 print(self.print_prefix + "ri+1: " + repr(ri_new)) 168 print(self.print_prefix + "bi+1: " + repr(bi_new)) 169 print(self.print_prefix + "pi+1: " + repr(pi_new)) 170 171 # Update i+1 to i. 172 xi = xi_new * 1.0 173 ri = ri_new * 1.0 174 pi = pi_new * 1.0 175 dot_ri = dot_ri_new 176 i = i + 1

177 178

179 - def new_param_func(self):

180 """The new parameter function. 181 182 Find the search direction, do a line search, and get xk+1 and fk+1. 183 """ 184 185 # Get the pk vector. 186 self.pk = self.get_pk() 187 188 # Line search. 189 self.line_search() 190 191 # Find the new parameter vector and function value at that point. 192 self.xk_new = self.xk + self.alpha * self.pk 193 self.fk_new, self.f_count = self.func(*(self.xk_new,)+self.args), self.f_count + 1 194 self.dfk_new, self.g_count = self.dfunc(*(self.xk_new,)+self.args), self.g_count + 1

195 196

197 - def update(self):

198 """Function to update the function value, gradient vector, and Hessian matrix.""" 199 200 self.xk = self.xk_new * 1.0 201 self.fk = self.fk_new 202 self.dfk = self.dfk_new * 1.0 203 self.d2fk, self.h_count = self.d2func(*(self.xk,)+self.args), self.h_count + 1

Source Code for Module minfx.ncg