/[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 4040 - (hide annotations)
Sun Oct 28 22:37:45 2012 UTC (7 years ago) by jfenwick
File MIME type: text/x-python
File size: 16560 byte(s)
Some py3 changes that didn't get committed for some reason

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     args_a=phiargs(alpha)
120     phi_a=phi(alpha, *args_a)
121     lslogger.debug("iteration %d, alpha=%e, phi(alpha)=%e"%(i,alpha,phi_a))
122     if (phi_a > phi0+c1*alpha*gphi0) or ((phi_a>=old_phi_a) and (i>1)):
123     alpha, phi_a, gphi_a = _zoom(phi, gradphi, phiargs, old_alpha, alpha, old_phi_a, phi_a, c1, c2, phi0, gphi0)
124     break
125     gphi_a=gradphi(alpha, *args_a)
126     if np.abs(gphi_a) <= -c2*gphi0:
127     break
128     if gphi_a >= 0:
129     alpha, phi_a, gphi_a = _zoom(phi, gradphi, phiargs, alpha, old_alpha, phi_a, old_phi_a, c1, c2, phi0, gphi0)
130     break
131    
132     old_alpha=alpha
133     # the factor is arbitrary as long as there is sufficient increase
134     alpha=2.*alpha
135     old_phi_a=phi_a
136     i+=1
137    
138     return alpha, phi_a, gf_new[0]
139    
140     ##############################################################################
141     class AbstractMinimizer(object):
142     """
143     Base class for function minimization methods.
144     """
145    
146     TOLERANCE_REACHED, MAX_ITERATIONS_REACHED=list(xrange(2))
147    
148     def __init__(self, f, tol=1e-5, imax=300):
149     """
150     Initializes a new minimizer for a given cost function.
151    
152     :param f: the cost function to be minimized
153     :type f: CostFunction
154     """
155     self._f=f
156     self._tol = tol
157     self._imax = imax
158     self._result = None
159     self._callback = None
160     self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)
161    
162     def setTolerance(self, tol):
163     """
164     Sets the tolerance for the stopping criterion. The minimizer stops when
165     an appropriate norm is less than `tol`.
166     """
167     self._tol = tol
168    
169     def setMaxIterations(self, imax):
170     """
171     Sets the maximum number of iterations before the minimizer terminates.
172     """
173     self._imax = imax
174    
175     def setCallback(self, callback):
176     """
177     Sets a callback function to be called after every iteration.
178     The arguments to the function are: (k, x, fx, gfx), where
179     k is the current iteration, x is the current estimate, fx=f(x) and
180     gfx=grad f(x).
181     """
182 caltinay 3948 if callback is not None and not callable(callback):
183 caltinay 3946 raise TypeError("Callback function not callable.")
184     self._callback = callback
185    
186     def _doCallback(self, *args):
187     if self._callback is not None:
188     self._callback(*args)
189    
190     def getResult(self):
191     """
192     Returns the result of the minimization.
193     """
194     return self._result
195    
196     def getOptions(self):
197     """
198     Returns a dictionary of minimizer-specific options.
199     """
200     return {}
201    
202     def setOptions(self, **opts):
203     """
204     Sets minimizer-specific options. For a list of possible options see
205     `getOptions()`.
206     """
207     raise NotImplementedError
208    
209     def run(self, x0):
210     """
211 caltinay 3990 Executes the minimization algorithm for *f* starting with the initial
212     guess ``x0``.
213    
214     :return: `TOLERANCE_REACHED` or `MAX_ITERATIONS_REACHED`
215 caltinay 3946 """
216     raise NotImplementedError
217    
218     def logSummary(self):
219     """
220     Outputs a summary of the completed minimization process to the logger.
221     """
222     self.logger.warning("Number of function evaluations: %d"%self._f.Value_calls)
223     self.logger.warning("Number of gradient evaluations: %d"%self._f.Gradient_calls)
224     self.logger.warning("Number of inner product evaluations: %d"%self._f.Inner_calls)
225     self.logger.warning("Number of argument evaluations: %d"%self._f.Arguments_calls)
226    
227    
228     ##############################################################################
229     class MinimizerLBFGS(AbstractMinimizer):
230     """
231     Minimizer that uses the limited-memory Broyden-Fletcher-Goldfarb-Shanno
232     method.
233     """
234    
235     # History size
236     _m = 15
237    
238     # Initial Hessian multiplier
239     _initial_H = 1
240    
241     def getOptions(self):
242     return {'historySize':self._m,'initialHessian':self._initial_H}
243    
244     def setOptions(self, **opts):
245     for o in opts:
246     if o=='historySize':
247     self._m=opts[o]
248     elif o=='initialHessian':
249     self._initial_H=opts[o]
250     else:
251     raise KeyError("Invalid option '%s'"%o)
252    
253     def getNorm(self, x):
254     return sqrt(self._f.getInner(x, x))
255    
256     def run(self, x):
257     args=self._f.getArguments(x)
258     gf=self._f.getGradient(x, *args)
259     fx=self._f(x, *args)
260     k=0
261     H=self._initial_H
262     error=2*self._tol
263     s_and_y=[]
264    
265     self._doCallback(k, x, fx, gf)
266    
267     while error > self._tol and k < self._imax:
268 caltinay 3948 #self.logger.info("\033[1;31miteration %d\033[1;30m, error=%e"%(k,error))
269     self.logger.info("iteration %d, error=%e"%(k,error))
270 caltinay 3946 # determine search direction
271     p = -self._twoLoop(H, gf, s_and_y)
272    
273     self.logger.debug("H = %s"%H)
274     self.logger.debug("grad f(x) = %s"%gf)
275     self.logger.debug("p = %s"%p)
276     self.logger.debug("x = %s"%x)
277    
278     # determine step length
279     alpha, fx_new, gf_new = line_search(self._f, x, p, gf, fx)
280     self.logger.debug("alpha=%e"%(alpha))
281     # execute the step
282     x_new = x + alpha*p
283     if gf_new is None:
284     gf_new=self._f.getGradient(x_new)
285     delta_x=x_new-x
286     delta_g=gf_new-gf
287     s_and_y.append((delta_x,delta_g))
288 caltinay 3950 if fx_new==0.:
289     error=fx
290     else:
291     error=abs(fx_new-fx)/abs(fx_new)
292 caltinay 3946 x=x_new
293     gf=gf_new
294     fx=fx_new
295     k+=1
296     self._doCallback(k, x, fx, gf)
297     if (error<=self._tol): break
298    
299     # delete oldest vector pair
300     if k>self._m: s_and_y.pop(0)
301    
302     # set the new scaling factor (approximation of inverse Hessian)
303     gnorm=self.getNorm(gf)
304     denom=self.getNorm(delta_g)
305     if denom < EPSILON * gnorm:
306     H=1.
307     self.logger.debug("Break down in H update. Resetting to 1.")
308     else:
309     H=self._f.getInner(delta_x,delta_g)/denom**2
310    
311     if k >= self._imax:
312     reason=self.MAX_ITERATIONS_REACHED
313     self.logger.warning("Maximum number of iterations reached!")
314     else:
315     reason=self.TOLERANCE_REACHED
316     self.logger.warning("Success after %d iterations! Final error=%e"%(k,error))
317     self._result=x
318     return reason
319    
320     def _twoLoop(self, H, gf, s_and_y):
321     """
322     Helper for the L-BFGS method.
323     See 'Numerical Optimization' by J. Nocedal for an explanation.
324     """
325     q=gf
326     alpha=[]
327     for s,y in reversed(s_and_y):
328     rho=1./(self._f.getInner(s, y))
329     a=rho*self._f.getInner(s, q)
330     alpha.append(a)
331     q=q-a*y
332    
333     r=H*q
334     for s,y in s_and_y:
335     rho=1./(self._f.getInner(s, y))
336     beta=rho*self._f.getInner(y,r)
337     a=alpha.pop()
338     r=r+s*(a-beta)
339     return r
340    
341     ##############################################################################
342     class MinimizerBFGS(AbstractMinimizer):
343     """
344     Minimizer that uses the Broyden-Fletcher-Goldfarb-Shanno method.
345     """
346    
347 caltinay 3950 # Initial Hessian multiplier
348     _initial_H = 1
349    
350     def getOptions(self):
351     return {'initialHessian':self._initial_H}
352    
353     def setOptions(self, **opts):
354     for o in opts:
355     if o=='initialHessian':
356     self._initial_H=opts[o]
357     else:
358     raise KeyError("Invalid option '%s'"%o)
359    
360 caltinay 3946 def run(self, x):
361 caltinay 3950 args=self._f.getArguments(x)
362     gf=self._f.getGradient(x, *args)
363     fx=self._f(x, *args)
364 caltinay 3946 k=0
365     try:
366     n=len(x)
367     except:
368     n=x.getNumberOfDataPoints()
369     I=np.eye(n)
370 caltinay 3950 H=self._initial_H*I
371 caltinay 3946 gnorm=Lsup(gf)
372 caltinay 3950 self._doCallback(k, x, fx, gf)
373    
374 caltinay 3946 while gnorm > self._tol and k < self._imax:
375     self.logger.info("iteration %d, gnorm=%e"%(k,gnorm))
376    
377     # determine search direction
378     d=-self._f.getInner(H, gf)
379 caltinay 3950
380     self.logger.debug("H = %s"%H)
381     self.logger.debug("grad f(x) = %s"%gf)
382     self.logger.debug("d = %s"%d)
383     self.logger.debug("x = %s"%x)
384    
385 caltinay 3946 # determine step length
386     alpha, fx, gf_new = line_search(self._f, x, d, gf, fx)
387     self.logger.debug("alpha=%e"%alpha)
388     # execute the step
389     x_new=x+alpha*d
390     delta_x=x_new-x
391     x=x_new
392     if gf_new is None:
393     gf_new=self._f.getGradient(x_new)
394     delta_g=gf_new-gf
395     gf=gf_new
396     k+=1
397 caltinay 3950 self._doCallback(k, x, fx, gf)
398 caltinay 3946 gnorm=Lsup(gf)
399     if (gnorm<=self._tol): break
400    
401     # update Hessian
402 caltinay 3950 denom=self._f.getInner(delta_x, delta_g)
403     if denom < EPSILON * gnorm:
404     denom=1e-5
405     self.logger.debug("Break down in H update. Resetting.")
406     rho=1./denom
407 caltinay 3946 self.logger.debug("rho=%e"%rho)
408 caltinay 3950 A=I-rho*delta_x[:,None]*delta_g[None,:]
409     AT=I-rho*delta_g[:,None]*delta_x[None,:]
410 caltinay 3946 H=self._f.getInner(A, self._f.getInner(H,AT)) + rho*delta_x[:,None]*delta_x[None,:]
411     if k >= self._imax:
412     reason=self.MAX_ITERATIONS_REACHED
413     self.logger.warning("Maximum number of iterations reached!")
414     else:
415     reason=self.TOLERANCE_REACHED
416     self.logger.warning("Success after %d iterations! Final gnorm=%e"%(k,gnorm))
417    
418     self._result=x
419     return reason
420    
421     ##############################################################################
422     class MinimizerNLCG(AbstractMinimizer):
423     """
424     Minimizer that uses the nonlinear conjugate gradient method
425 caltinay 3950 (Fletcher-Reeves variant).
426 caltinay 3946 """
427    
428     def run(self, x):
429     i=0
430     k=0
431 caltinay 3950 args=self._f.getArguments(x)
432     r=-self._f.getGradient(x, *args)
433     fx=self._f(x, *args)
434 caltinay 3946 d=r
435     delta=self._f.getInner(r,r)
436     delta0=delta
437 caltinay 3950 self._doCallback(i, x, fx, -r)
438    
439     while i<self._imax and Lsup(r)>self._tol:
440     self.logger.info("iteration %d"%i)
441     self.logger.debug("grad f(x) = %s"%(-r))
442     self.logger.debug("d = %s"%d)
443     self.logger.debug("x = %s"%x)
444    
445     alpha, fx, gf_new = line_search(self._f, x, d, -r, fx, c2=0.4)
446 caltinay 3946 self.logger.debug("alpha=%e"%(alpha))
447     x=x+alpha*d
448     r=-self._f.getGradient(x) if gf_new is None else -gf_new
449     delta_o=delta
450     delta=self._f.getInner(r,r)
451     beta=delta/delta_o
452     d=r+beta*d
453     k=k+1
454     try:
455     lenx=len(x)
456     except:
457     lenx=x.getNumberOfDataPoints()
458     if k == lenx or self._f.getInner(r,d) <= 0:
459     d=r
460     k=0
461     i+=1
462 caltinay 3950 self._doCallback(i, x, fx, gf_new)
463 caltinay 3946
464     if i >= self._imax:
465     reason=self.MAX_ITERATIONS_REACHED
466     self.logger.warning("Maximum number of iterations reached!")
467     else:
468     reason=self.TOLERANCE_REACHED
469     self.logger.warning("Success after %d iterations! Final delta=%e"%(i,delta))
470    
471     self._result=x
472     return reason
473    
474    
475     if __name__=="__main__":
476     # Example usage with function 'rosen' (minimum=[1,1,...1]):
477     from scipy.optimize import rosen, rosen_der
478 caltinay 3947 from esys.downunder.costfunctions import CostFunction
479 caltinay 3946 import sys
480     N=100
481     x0=np.array([4.]*N) # initial guess
482    
483     class RosenFunc(CostFunction):
484     def _getInner(self, f0, f1):
485     return np.dot(f0, f1)
486     def _getValue(self, x, *args):
487     return rosen(x)
488     def _getGradient(self, x, *args):
489     return rosen_der(x)
490    
491     f=RosenFunc()
492     m=None
493     if len(sys.argv)>1:
494     method=sys.argv[1].lower()
495     if method=='nlcg':
496     m=MinimizerNLCG(f)
497     elif method=='bfgs':
498     m=MinimizerBFGS(f)
499    
500     if m is None:
501     # default
502     m=MinimizerLBFGS(f)
503     #m.setOptions(historySize=10000)
504    
505 caltinay 3950 logging.basicConfig(format='[%(funcName)s] \033[1;30m%(message)s\033[0m', level=logging.INFO)
506 caltinay 3946 m.setTolerance(1e-5)
507 caltinay 3950 m.setMaxIterations(600)
508 caltinay 3946 m.run(x0)
509     m.logSummary()
510     print("\tLsup(result)=%.8f"%np.amax(abs(m.getResult())))
511    
512     #from scipy.optimize import fmin_cg
513     #print("scipy ref=%.8f"%np.amax(abs(fmin_cg(rosen, x0, rosen_der, maxiter=10000))))
514    

  ViewVC Help
Powered by ViewVC 1.1.26