/[escript]/trunk/downunder/py_src/minimizers.py
ViewVC logotype

Annotation of /trunk/downunder/py_src/minimizers.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4047 - (hide annotations)
Tue Oct 30 09:10:07 2012 UTC (7 years ago) by gross
File MIME type: text/x-python
File size: 16607 byte(s)
More on magnetic inversion
1 caltinay 3946
2 jfenwick 3981 ##############################################################################
3 caltinay 3946 #
4     # Copyright (c) 2003-2012 by University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 caltinay 3946 #
7     # Primary Business: Queensland, Australia
8     # Licensed under the Open Software License version 3.0
9     # http://www.opensource.org/licenses/osl-3.0.php
10     #
11 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     # Development since 2012 by School of Earth Sciences
13     #
14     ##############################################################################
15 caltinay 3946
16 caltinay 4007 """Generic minimization algorithms"""
17    
18 caltinay 3946 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
19 jfenwick 3981 http://www.uq.edu.au
20 caltinay 3946 Primary Business: Queensland, Australia"""
21     __license__="""Licensed under the Open Software License version 3.0
22     http://www.opensource.org/licenses/osl-3.0.php"""
23     __url__="https://launchpad.net/escript-finley"
24    
25 caltinay 3947 __all__ = ['AbstractMinimizer', 'MinimizerLBFGS', 'MinimizerBFGS', 'MinimizerNLCG']
26    
27 caltinay 3946 import logging
28     import numpy as np
29     try:
30     from esys.escript import Lsup, sqrt, EPSILON
31     except:
32     Lsup=lambda x: np.amax(abs(x))
33     sqrt=np.sqrt
34     EPSILON=1e-18
35    
36 jfenwick 4040 import sys
37     if sys.version_info[0]>2:
38     xrange=range
39    
40 caltinay 3946 lslogger=logging.getLogger('inv.minimizer.linesearch')
41     zoomlogger=logging.getLogger('inv.minimizer.linesearch.zoom')
42    
43     def _zoom(phi, gradphi, phiargs, alpha_lo, alpha_hi, phi_lo, phi_hi, c1, c2, phi0, gphi0, IMAX=25):
44     """
45     Helper function for `line_search` below which tries to tighten the range
46     alpha_lo...alpha_hi. See Chapter 3 of 'Numerical Optimization' by
47     J. Nocedal for an explanation.
48     """
49     i=0
50     while True:
51     alpha=alpha_lo+.5*(alpha_hi-alpha_lo) # should use interpolation...
52     args_a=phiargs(alpha)
53     phi_a=phi(alpha, *args_a)
54     zoomlogger.debug("iteration %d, alpha=%e, phi(alpha)=%e"%(i,alpha,phi_a))
55     if phi_a > phi0+c1*alpha*gphi0 or phi_a >= phi_lo:
56     alpha_hi=alpha
57     else:
58     gphi_a=gradphi(alpha, *args_a)
59     zoomlogger.debug("grad(phi(alpha))=%e"%(gphi_a))
60     if np.abs(gphi_a) <= -c2*gphi0:
61     break
62     if gphi_a*(alpha_hi-alpha_lo) >= 0:
63     alpha_hi = alpha_lo
64     alpha_lo=alpha
65     phi_lo=phi_a
66     i+=1
67     if i>IMAX:
68     gphi_a=None
69     break
70     return alpha, phi_a, gphi_a
71    
72     def line_search(f, x, p, gf, fx, alpha_max=50.0, c1=1e-4, c2=0.9, IMAX=15):
73     """
74     Line search method that satisfies the strong Wolfe conditions.
75     See Chapter 3 of 'Numerical Optimization' by J. Nocedal for an explanation.
76    
77     :param f: callable objective function f(x)
78     :param x: start value for the line search
79     :param p: search direction
80     :param gf: value for the gradient of f at x
81     :param fx: value of f(x)
82     :param alpha_max: algorithm terminates if alpha reaches this value
83     :param c1: value for Armijo condition (see reference)
84     :param c2: value for curvature condition (see reference)
85     :param IMAX: maximum number of iterations to perform
86     """
87     # this stores the latest gradf(x+a*p) which is returned
88     gf_new=[gf]
89    
90     def phi(a, *args):
91     """ phi(a):=f(x+a*p) """
92     return f(x+a*p, *args)
93     def gradphi(a, *args):
94     gf_new[0]=f.getGradient(x+a*p, *args)
95     return f.getInner(gf_new[0], p)
96     #return f.getDirectionalDerivative(x+a*p, p, *args)
97     def phiargs(a):
98     try:
99     args=f.getArguments(x+a*p)
100     except:
101     args=()
102     return args
103    
104     old_alpha=0.
105     # we assume gf is properly scaled so alpha=1 is a reasonable starting value
106     alpha=1.
107     if fx is None:
108     args0=phiargs(0.)
109     phi0=phi(0., *args0)
110     else:
111     phi0=fx
112     lslogger.debug("phi(0)=%e"%(phi0))
113     gphi0=f.getInner(gf, p) #gradphi(0., *args0)
114     lslogger.debug("grad phi(0)=%e"%(gphi0))
115     old_phi_a=phi0
116     i=1
117    
118     while i<IMAX and alpha>0. and alpha<alpha_max:
119 gross 4047
120 caltinay 3946 args_a=phiargs(alpha)
121     phi_a=phi(alpha, *args_a)
122     lslogger.debug("iteration %d, alpha=%e, phi(alpha)=%e"%(i,alpha,phi_a))
123     if (phi_a > phi0+c1*alpha*gphi0) or ((phi_a>=old_phi_a) and (i>1)):
124     alpha, phi_a, gphi_a = _zoom(phi, gradphi, phiargs, old_alpha, alpha, old_phi_a, phi_a, c1, c2, phi0, gphi0)
125     break
126 gross 4047
127 caltinay 3946 gphi_a=gradphi(alpha, *args_a)
128     if np.abs(gphi_a) <= -c2*gphi0:
129     break
130     if gphi_a >= 0:
131     alpha, phi_a, gphi_a = _zoom(phi, gradphi, phiargs, alpha, old_alpha, phi_a, old_phi_a, c1, c2, phi0, gphi0)
132     break
133    
134     old_alpha=alpha
135     # the factor is arbitrary as long as there is sufficient increase
136     alpha=2.*alpha
137     old_phi_a=phi_a
138     i+=1
139    
140     return alpha, phi_a, gf_new[0]
141    
142     ##############################################################################
143     class AbstractMinimizer(object):
144     """
145     Base class for function minimization methods.
146     """
147    
148     TOLERANCE_REACHED, MAX_ITERATIONS_REACHED=list(xrange(2))
149    
150     def __init__(self, f, tol=1e-5, imax=300):
151     """
152     Initializes a new minimizer for a given cost function.
153    
154     :param f: the cost function to be minimized
155     :type f: CostFunction
156     """
157     self._f=f
158     self._tol = tol
159     self._imax = imax
160     self._result = None
161     self._callback = None
162     self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)
163    
164     def setTolerance(self, tol):
165     """
166     Sets the tolerance for the stopping criterion. The minimizer stops when
167     an appropriate norm is less than `tol`.
168     """
169     self._tol = tol
170    
171     def setMaxIterations(self, imax):
172     """
173     Sets the maximum number of iterations before the minimizer terminates.
174     """
175     self._imax = imax
176    
177     def setCallback(self, callback):
178     """
179     Sets a callback function to be called after every iteration.
180     The arguments to the function are: (k, x, fx, gfx), where
181     k is the current iteration, x is the current estimate, fx=f(x) and
182     gfx=grad f(x).
183     """
184 caltinay 3948 if callback is not None and not callable(callback):
185 caltinay 3946 raise TypeError("Callback function not callable.")
186     self._callback = callback
187    
188     def _doCallback(self, *args):
189     if self._callback is not None:
190     self._callback(*args)
191    
192     def getResult(self):
193     """
194     Returns the result of the minimization.
195     """
196     return self._result
197    
198     def getOptions(self):
199     """
200     Returns a dictionary of minimizer-specific options.
201     """
202     return {}
203    
204     def setOptions(self, **opts):
205     """
206     Sets minimizer-specific options. For a list of possible options see
207     `getOptions()`.
208     """
209     raise NotImplementedError
210    
211     def run(self, x0):
212     """
213 caltinay 3990 Executes the minimization algorithm for *f* starting with the initial
214     guess ``x0``.
215    
216     :return: `TOLERANCE_REACHED` or `MAX_ITERATIONS_REACHED`
217 caltinay 3946 """
218     raise NotImplementedError
219    
220     def logSummary(self):
221     """
222     Outputs a summary of the completed minimization process to the logger.
223     """
224     self.logger.warning("Number of function evaluations: %d"%self._f.Value_calls)
225     self.logger.warning("Number of gradient evaluations: %d"%self._f.Gradient_calls)
226     self.logger.warning("Number of inner product evaluations: %d"%self._f.Inner_calls)
227     self.logger.warning("Number of argument evaluations: %d"%self._f.Arguments_calls)
228    
229    
230     ##############################################################################
231     class MinimizerLBFGS(AbstractMinimizer):
232     """
233     Minimizer that uses the limited-memory Broyden-Fletcher-Goldfarb-Shanno
234     method.
235     """
236    
237     # History size
238     _m = 15
239    
240     # Initial Hessian multiplier
241     _initial_H = 1
242    
243     def getOptions(self):
244     return {'historySize':self._m,'initialHessian':self._initial_H}
245    
246     def setOptions(self, **opts):
247     for o in opts:
248     if o=='historySize':
249     self._m=opts[o]
250     elif o=='initialHessian':
251     self._initial_H=opts[o]
252     else:
253     raise KeyError("Invalid option '%s'"%o)
254    
255     def getNorm(self, x):
256     return sqrt(self._f.getInner(x, x))
257    
258     def run(self, x):
259     args=self._f.getArguments(x)
260     gf=self._f.getGradient(x, *args)
261     fx=self._f(x, *args)
262     k=0
263     H=self._initial_H
264     error=2*self._tol
265     s_and_y=[]
266    
267     self._doCallback(k, x, fx, gf)
268    
269     while error > self._tol and k < self._imax:
270 caltinay 3948 #self.logger.info("\033[1;31miteration %d\033[1;30m, error=%e"%(k,error))
271     self.logger.info("iteration %d, error=%e"%(k,error))
272 caltinay 3946 # determine search direction
273     p = -self._twoLoop(H, gf, s_and_y)
274    
275     self.logger.debug("H = %s"%H)
276     self.logger.debug("grad f(x) = %s"%gf)
277     self.logger.debug("p = %s"%p)
278     self.logger.debug("x = %s"%x)
279    
280     # determine step length
281     alpha, fx_new, gf_new = line_search(self._f, x, p, gf, fx)
282     self.logger.debug("alpha=%e"%(alpha))
283     # execute the step
284     x_new = x + alpha*p
285     if gf_new is None:
286     gf_new=self._f.getGradient(x_new)
287     delta_x=x_new-x
288     delta_g=gf_new-gf
289     s_and_y.append((delta_x,delta_g))
290 caltinay 3950 if fx_new==0.:
291     error=fx
292     else:
293     error=abs(fx_new-fx)/abs(fx_new)
294 caltinay 3946 x=x_new
295     gf=gf_new
296     fx=fx_new
297     k+=1
298     self._doCallback(k, x, fx, gf)
299     if (error<=self._tol): break
300    
301     # delete oldest vector pair
302     if k>self._m: s_and_y.pop(0)
303    
304     # set the new scaling factor (approximation of inverse Hessian)
305     gnorm=self.getNorm(gf)
306     denom=self.getNorm(delta_g)
307     if denom < EPSILON * gnorm:
308     H=1.
309     self.logger.debug("Break down in H update. Resetting to 1.")
310     else:
311     H=self._f.getInner(delta_x,delta_g)/denom**2
312    
313     if k >= self._imax:
314     reason=self.MAX_ITERATIONS_REACHED
315     self.logger.warning("Maximum number of iterations reached!")
316     else:
317     reason=self.TOLERANCE_REACHED
318     self.logger.warning("Success after %d iterations! Final error=%e"%(k,error))
319     self._result=x
320     return reason
321    
322     def _twoLoop(self, H, gf, s_and_y):
323     """
324     Helper for the L-BFGS method.
325     See 'Numerical Optimization' by J. Nocedal for an explanation.
326     """
327     q=gf
328     alpha=[]
329     for s,y in reversed(s_and_y):
330 gross 4047 print "TT",s, y, self._f.getInner(s, y)
331 caltinay 3946 rho=1./(self._f.getInner(s, y))
332     a=rho*self._f.getInner(s, q)
333     alpha.append(a)
334     q=q-a*y
335    
336     r=H*q
337     for s,y in s_and_y:
338     rho=1./(self._f.getInner(s, y))
339     beta=rho*self._f.getInner(y,r)
340     a=alpha.pop()
341     r=r+s*(a-beta)
342     return r
343    
344     ##############################################################################
345     class MinimizerBFGS(AbstractMinimizer):
346     """
347     Minimizer that uses the Broyden-Fletcher-Goldfarb-Shanno method.
348     """
349    
350 caltinay 3950 # Initial Hessian multiplier
351     _initial_H = 1
352    
353     def getOptions(self):
354     return {'initialHessian':self._initial_H}
355    
356     def setOptions(self, **opts):
357     for o in opts:
358     if o=='initialHessian':
359     self._initial_H=opts[o]
360     else:
361     raise KeyError("Invalid option '%s'"%o)
362    
363 caltinay 3946 def run(self, x):
364 caltinay 3950 args=self._f.getArguments(x)
365     gf=self._f.getGradient(x, *args)
366     fx=self._f(x, *args)
367 caltinay 3946 k=0
368     try:
369     n=len(x)
370     except:
371     n=x.getNumberOfDataPoints()
372     I=np.eye(n)
373 caltinay 3950 H=self._initial_H*I
374 caltinay 3946 gnorm=Lsup(gf)
375 caltinay 3950 self._doCallback(k, x, fx, gf)
376    
377 caltinay 3946 while gnorm > self._tol and k < self._imax:
378     self.logger.info("iteration %d, gnorm=%e"%(k,gnorm))
379    
380     # determine search direction
381     d=-self._f.getInner(H, gf)
382 caltinay 3950
383     self.logger.debug("H = %s"%H)
384     self.logger.debug("grad f(x) = %s"%gf)
385     self.logger.debug("d = %s"%d)
386     self.logger.debug("x = %s"%x)
387    
388 caltinay 3946 # determine step length
389     alpha, fx, gf_new = line_search(self._f, x, d, gf, fx)
390     self.logger.debug("alpha=%e"%alpha)
391     # execute the step
392     x_new=x+alpha*d
393     delta_x=x_new-x
394     x=x_new
395     if gf_new is None:
396     gf_new=self._f.getGradient(x_new)
397     delta_g=gf_new-gf
398     gf=gf_new
399     k+=1
400 caltinay 3950 self._doCallback(k, x, fx, gf)
401 caltinay 3946 gnorm=Lsup(gf)
402     if (gnorm<=self._tol): break
403    
404     # update Hessian
405 caltinay 3950 denom=self._f.getInner(delta_x, delta_g)
406     if denom < EPSILON * gnorm:
407     denom=1e-5
408     self.logger.debug("Break down in H update. Resetting.")
409     rho=1./denom
410 caltinay 3946 self.logger.debug("rho=%e"%rho)
411 caltinay 3950 A=I-rho*delta_x[:,None]*delta_g[None,:]
412     AT=I-rho*delta_g[:,None]*delta_x[None,:]
413 caltinay 3946 H=self._f.getInner(A, self._f.getInner(H,AT)) + rho*delta_x[:,None]*delta_x[None,:]
414     if k >= self._imax:
415     reason=self.MAX_ITERATIONS_REACHED
416     self.logger.warning("Maximum number of iterations reached!")
417     else:
418     reason=self.TOLERANCE_REACHED
419     self.logger.warning("Success after %d iterations! Final gnorm=%e"%(k,gnorm))
420    
421     self._result=x
422     return reason
423    
424     ##############################################################################
425     class MinimizerNLCG(AbstractMinimizer):
426     """
427     Minimizer that uses the nonlinear conjugate gradient method
428 caltinay 3950 (Fletcher-Reeves variant).
429 caltinay 3946 """
430    
431     def run(self, x):
432     i=0
433     k=0
434 caltinay 3950 args=self._f.getArguments(x)
435     r=-self._f.getGradient(x, *args)
436     fx=self._f(x, *args)
437 caltinay 3946 d=r
438     delta=self._f.getInner(r,r)
439     delta0=delta
440 caltinay 3950 self._doCallback(i, x, fx, -r)
441    
442     while i<self._imax and Lsup(r)>self._tol:
443     self.logger.info("iteration %d"%i)
444     self.logger.debug("grad f(x) = %s"%(-r))
445     self.logger.debug("d = %s"%d)
446     self.logger.debug("x = %s"%x)
447    
448     alpha, fx, gf_new = line_search(self._f, x, d, -r, fx, c2=0.4)
449 caltinay 3946 self.logger.debug("alpha=%e"%(alpha))
450     x=x+alpha*d
451     r=-self._f.getGradient(x) if gf_new is None else -gf_new
452     delta_o=delta
453     delta=self._f.getInner(r,r)
454     beta=delta/delta_o
455     d=r+beta*d
456     k=k+1
457     try:
458     lenx=len(x)
459     except:
460     lenx=x.getNumberOfDataPoints()
461     if k == lenx or self._f.getInner(r,d) <= 0:
462     d=r
463     k=0
464     i+=1
465 caltinay 3950 self._doCallback(i, x, fx, gf_new)
466 caltinay 3946
467     if i >= self._imax:
468     reason=self.MAX_ITERATIONS_REACHED
469     self.logger.warning("Maximum number of iterations reached!")
470     else:
471     reason=self.TOLERANCE_REACHED
472     self.logger.warning("Success after %d iterations! Final delta=%e"%(i,delta))
473    
474     self._result=x
475     return reason
476    
477    
478     if __name__=="__main__":
479     # Example usage with function 'rosen' (minimum=[1,1,...1]):
480     from scipy.optimize import rosen, rosen_der
481 caltinay 3947 from esys.downunder.costfunctions import CostFunction
482 caltinay 3946 import sys
483     N=100
484     x0=np.array([4.]*N) # initial guess
485    
486     class RosenFunc(CostFunction):
487     def _getInner(self, f0, f1):
488     return np.dot(f0, f1)
489     def _getValue(self, x, *args):
490     return rosen(x)
491     def _getGradient(self, x, *args):
492     return rosen_der(x)
493    
494     f=RosenFunc()
495     m=None
496     if len(sys.argv)>1:
497     method=sys.argv[1].lower()
498     if method=='nlcg':
499     m=MinimizerNLCG(f)
500     elif method=='bfgs':
501     m=MinimizerBFGS(f)
502    
503     if m is None:
504     # default
505     m=MinimizerLBFGS(f)
506     #m.setOptions(historySize=10000)
507    
508 caltinay 3950 logging.basicConfig(format='[%(funcName)s] \033[1;30m%(message)s\033[0m', level=logging.INFO)
509 caltinay 3946 m.setTolerance(1e-5)
510 caltinay 3950 m.setMaxIterations(600)
511 caltinay 3946 m.run(x0)
512     m.logSummary()
513     print("\tLsup(result)=%.8f"%np.amax(abs(m.getResult())))
514    
515     #from scipy.optimize import fmin_cg
516     #print("scipy ref=%.8f"%np.amax(abs(fmin_cg(rosen, x0, rosen_der, maxiter=10000))))
517    

  ViewVC Help
Powered by ViewVC 1.1.26