linear_operators.iterative.algorithms

1 """ 2 Implements iterative algorithm class. 3 """ 4 import numpy as np 5 from copy import copy 6 7 from linesearch import * 8 from criterions import * 9 from norms import norm2 10 11 # defaults 12 TOL = 1e-6 13 GTOL = 1e-6 14 MAXITER = None 15 16 # stop conditions

17 -class StopCondition(object):

18 - def _test_maxiter(self, algo):

19 return algo.iter_ >= self.maxiter

20 - def _test_tol(self, algo):

21 self.resid = np.abs(algo.last_criterion - algo.current_criterion) 22 self.resid /= algo.first_criterion 23 return self.resid < self.tol

24 - def _test_gtol(self, algo):

25 return algo.current_gradient_norm < self.gtol

26 _all_tests = [_test_maxiter, _test_tol, _test_gtol]

27 - def __init__(self, maxiter=None, tol=None, gtol=None, cond=np.any):

28 self.cond = cond 29 self.maxiter = maxiter 30 self.tol = tol 31 self.gtol = gtol 32 self.all_val = [self.maxiter, self.tol, self.gtol] 33 # filter out tests with None values 34 self.tests_val = [val for val in self.all_val 35 if val is not None] 36 self.tests = [test 37 for test, val in zip(self._all_tests, self.all_val) 38 if val is not None] 39 # store values for printing 40 self.resid = None

41 - def __call__(self, algo):

42 return self.cond([test(self, algo) for test in self.tests])

43 - def str(self, algo):

44 """ 45 Returns a string with current condition values. 46 """ 47 if self.resid is not None and self.tol is not None: 48 return "\t %1.2e / %1.2e" % (self.resid, self.tol) 49 else: 50 return "\t Residual"

51 52 default_stop = StopCondition(maxiter=MAXITER, tol=TOL, gtol=GTOL) 53 54 # update types 55

56 -def fletcher_reeves(algo):

57 """ 58 Fletcher-Reeves descent direction update method. 59 """ 60 return algo.current_gradient_norm / algo.last_gradient_norm

61

62 -def polak_ribiere(algo):

63 """ 64 Polak-Ribiere descent direction update method. 65 """ 66 b = np.dot(algo.current_gradient.T, 67 (algo.current_gradient - algo.last_gradient)) 68 b /= np.norm(algo.last_gradient) 69 return b

70 71 # callback function 72

73 -class Callback(object):

74 - def __init__(self, verbose=False, savefile=None, shape=()):

75 self.verbose = verbose 76 self.savefile = savefile 77 self.shape = shape 78 self.im = None

79 - def print_status(self, algo):

80 if self.verbose: 81 if algo.iter_ == 1: 82 print('Iteration \t Criterion') 83 print_str = "\t%i \t %e" % (algo.iter_, algo.current_criterion) 84 print_str += algo.stop_condition.str(algo) 85 print(print_str)

86 - def save(self, algo):

87 if self.savefile is not None: 88 var_dict = { 89 "iter":algo.iter_, 90 "criterion":algo.current_criterion, 91 "solution":algo.current_solution, 92 } 93 np.savez(self.savefile, **var_dict)

94 - def imshow(self, algo):

95 import pylab 96 if algo.iter_ == 1: 97 self.im = pylab.imshow(algo.current_solution.reshape(self.shape)) 98 else: 99 self.im.set_data(algo.current_solution.reshape(self.shape)) 100 pylab.draw() 101 pylab.show()

102 - def plot(self, algo):

103 import pylab 104 if algo.iter_ == 1: 105 self.im = pylab.plot(algo.current_solution)[0] 106 else: 107 y = algo.current_solution 108 self.im.set_ydata(y) 109 pylab.ylim((y.min(), y.max())) 110 pylab.draw() 111 pylab.show()

112 - def __call__(self, algo):

113 if self.verbose: 114 self.print_status(algo) 115 if self.savefile is not None: 116 self.save(algo) 117 if self.shape is not None: 118 if len(self.shape) == 1: 119 self.plot(algo) 120 elif len(self.shape) == 2: 121 self.imshow(algo)

122 123 default_callback = Callback() 124 125 # algorithms 126

127 -class Algorithm(object):

128 """ 129 Abstract class to define iterative algorithms. 130 131 Attributes 132 ---------- 133 134 iter_ : int 135 Current iteration number. 136 137 Methods 138 ------- 139 140 initialize : Set variables to initial state 141 iterate : perform one iteration and return current solution 142 callback : user-defined function to print status or save variables 143 cont : continue the optimization skipping initialiaztion 144 __call__ : perform the optimization unt stop_condition is reached 145 """

146 - def initialize(self):

147 self.iter_ = 0 148 self.current_solution = None

149 - def callback(self):

150 pass

151 - def iterate(self):

152 """ 153 Perform one iteration and returns current solution. 154 """ 155 self.iter_ += 1 156 self.callback(self) 157 # return value not used in loop but usefull in "interactive mode" 158 return self.current_solution

159 - def at_exit(self):

160 """ 161 Perform some task at exit. 162 Does nothing by default. 163 """ 164 pass

165 - def __call__(self):

166 """ 167 Perform the optimization. 168 """ 169 self.initialize() 170 self.iterate() # at least 1 iteration 171 self.cont() 172 self.at_exit() 173 return self.current_solution

174 - def cont(self):

175 """ 176 Continue an interrupted estimation (like call but avoid 177 initialization). 178 """ 179 while not self.stop_condition(self): 180 self.iterate() 181 return self.current_solution

182

183 -class ConjugateGradient(Algorithm):

184 """ 185 Apply the conjugate gradient algorithm to a Criterion instance. 186 187 Parameters 188 ---------- 189 190 criterion : Criterion 191 A Criterion instance. It should have following methods and attributes: 192 __call__ : returns criterion values at given point 193 gradient : returns gradient (1st derivative) of criterion at given point 194 n_variable: the size of the input vector of criterion 195 196 x0 : ndarray (None) 197 The first guess of the algorithm. 198 199 callback : function (default_callback) 200 Perform some printing / saving operations at each iteration. 201 202 stop_condition : function (default_stop) 203 Defines when the iterations should stop 204 205 update_type : function (fletcher_reeves) 206 Type of descent direction update : e.g. fletcher_reeves, polak_ribiere 207 208 line_search : function (optimal step) 209 Line search method to find the minimum along each direction at each 210 iteration. 211 212 Returns 213 ------- 214 215 Returns an algorithm instance. Optimization is performed by 216 calling the this instance. 217 218 """

219 - def __init__(self, criterion, x0=None, 220 callback=default_callback, 221 stop_condition=default_stop, 222 update_type=fletcher_reeves, 223 line_search=optimal_step, **kwargs):

224 self.criterion = criterion 225 self.gradient = criterion.gradient 226 self.n_variables = self.criterion.n_variables 227 # functions 228 self.callback = callback 229 self.stop_condition = stop_condition 230 self.update_type = update_type 231 self.line_search = line_search 232 self.kwargs = kwargs 233 # to store values 234 self.current_criterion = np.inf 235 self.current_solution = None 236 self.current_gradient = None 237 self.current_gradient_norm = None 238 self.current_descent = None 239 self.last_criterion = np.inf 240 self.last_solution = None 241 self.last_gradient = None 242 self.last_gradient_norm = None 243 self.last_descent = None

244 - def initialize(self):

245 """ 246 Initialize required values. 247 """ 248 Algorithm.initialize(self) 249 self.first_guess() 250 self.first_criterion = self.criterion(self.current_solution) 251 self.current_criterion = self.first_criterion

252 - def first_guess(self, x0=None):

253 """ 254 Sets current_solution attribute to initial value. 255 """ 256 if x0 is None: 257 self.current_solution = np.zeros(self.n_variables) 258 else: 259 self.current_solution = copy(x0)

260 # update_* functions encode the actual algorithm

261 - def update_gradient(self):

262 self.last_gradient = copy(self.current_gradient) 263 self.current_gradient = self.gradient(self.current_solution)

264 - def update_gradient_norm(self):

265 self.last_gradient_norm = copy(self.current_gradient_norm) 266 self.current_gradient_norm = norm2(self.current_gradient)

267 - def update_descent(self):

268 if self.iter_ == 0: 269 self.current_descent = - self.current_gradient 270 else: 271 self.last_descent = copy(self.current_descent) 272 b = self.update_type(self) 273 self.current_descent = - self.current_gradient + b * self.last_descent

274 - def update_solution(self):

275 self.last_solution = copy(self.current_solution) 276 a = self.line_search(self) 277 self.current_solution += a * self.current_descent

278 - def update_criterion(self):

279 self.last_criterion = copy(self.current_criterion) 280 self.current_criterion = self.criterion(self.current_solution)

281 - def iterate(self):

282 """ 283 Update all values. 284 """ 285 self.update_gradient() 286 self.update_gradient_norm() 287 self.update_descent() 288 self.update_solution() 289 self.update_criterion() 290 Algorithm.iterate(self)

291

292 -class QuadraticConjugateGradient(ConjugateGradient):

293 """ 294 A subclass of ConjugateGradient using a QuadraticCriterion. 295 """

296 - def __init__(self, model, data, priors=[], hypers=[], **kwargs):

297 store = kwargs.pop("store", True) 298 criterion = QuadraticCriterion(model, data, hypers=hypers, 299 priors=priors, store=store) 300 ConjugateGradient.__init__(self, criterion, **kwargs)

301

302 -class HuberConjugateGradient(ConjugateGradient):

303 """ 304 A subclass of ConjugateGradient using an HuberCriterion. 305 """

306 - def __init__(self, model, data, priors=[], hypers=[], deltas=None, **kwargs):

307 store = kwargs.pop("store", True) 308 criterion = HuberCriterion(model, data, hypers=hypers, priors=priors, 309 deltas=deltas, store=store) 310 ConjugateGradient.__init__(self, criterion, **kwargs)

311 312 # for backward compatibility 313

314 -def acg(model, data, priors=[], hypers=[], **kwargs):

315 algorithm = QuadraticConjugateGradient(model, data, priors=priors, 316 hypers=hypers, **kwargs) 317 sol = algorithm() 318 return sol

319

320 -def hacg(model, data, priors=[], hypers=[], deltas=None, **kwargs):

321 algorithm = HuberConjugateGradient(model, data, priors=priors, 322 hypers=hypers, deltas=deltas, **kwargs) 323 sol = algorithm() 324 return sol

325 326 # other 327

328 -def normalize_hyper(hyper, y, x):

329 """ 330 Normalize hyperparamaters so that they are independent of pb size 331 """ 332 nx = float(x.size) 333 ny = float(y.size) 334 return np.asarray(hyper) * ny / nx

335

Source Code for Module linear_operators.iterative.algorithms