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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4073 by caltinay, Wed Oct 31 01:01:24 2012 UTC revision 4074 by gross, Thu Nov 15 03:30:59 2012 UTC
# Line 26  __all__ = ['AbstractMinimizer', 'Minimiz Line 26  __all__ = ['AbstractMinimizer', 'Minimiz
26    
27  import logging  import logging
28  import numpy as np  import numpy as np
29    
30  try:  try:
31      from esys.escript import Lsup, sqrt, EPSILON      from esys.escript import Lsup, sqrt, EPSILON
32  except:  except:
# Line 92  def line_search(f, x, p, gf, fx, alpha_m Line 93  def line_search(f, x, p, gf, fx, alpha_m
93          return f(x+a*p, *args)          return f(x+a*p, *args)
94      def gradphi(a, *args):      def gradphi(a, *args):
95          gf_new[0]=f.getGradient(x+a*p, *args)          gf_new[0]=f.getGradient(x+a*p, *args)
96          return f.getInner(gf_new[0], p)          return f.getDualProduct(p, gf_new[0])
97          #return f.getDirectionalDerivative(x+a*p, p, *args)          #return f.getDirectionalDerivative(x+a*p, p, *args)
98      def phiargs(a):      def phiargs(a):
99          try:          try:
# Line 110  def line_search(f, x, p, gf, fx, alpha_m Line 111  def line_search(f, x, p, gf, fx, alpha_m
111      else:      else:
112          phi0=fx          phi0=fx
113      lslogger.debug("phi(0)=%e"%(phi0))      lslogger.debug("phi(0)=%e"%(phi0))
114      gphi0=f.getInner(gf, p) #gradphi(0., *args0)      gphi0=f.getDualProduct(p, gf) #gradphi(0., *args0)
115      lslogger.debug("grad phi(0)=%e"%(gphi0))      lslogger.debug("grad phi(0)=%e"%(gphi0))
116      old_phi_a=phi0      old_phi_a=phi0
117      i=1      i=1
# Line 221  class AbstractMinimizer(object): Line 222  class AbstractMinimizer(object):
222          """          """
223          Outputs a summary of the completed minimization process to the logger.          Outputs a summary of the completed minimization process to the logger.
224          """          """
225          self.logger.warning("Number of function evaluations: %d"%self._f.Value_calls)          if hasattr(self._f, "Value_calls"):
226          self.logger.warning("Number of gradient evaluations: %d"%self._f.Gradient_calls)        self.logger.warning("Number of function evaluations: %d"%self._f.Value_calls)
227          self.logger.warning("Number of inner product evaluations: %d"%self._f.Inner_calls)        self.logger.warning("Number of gradient evaluations: %d"%self._f.Gradient_calls)
228          self.logger.warning("Number of argument evaluations: %d"%self._f.Arguments_calls)        self.logger.warning("Number of inner product evaluations: %d"%self._f.DualProduct_calls)
229          self.logger.warning("Number of argument evaluations: %d"%self._f.Arguments_calls)
230    
231    
232  ##############################################################################  ##############################################################################
# Line 252  class MinimizerLBFGS(AbstractMinimizer): Line 254  class MinimizerLBFGS(AbstractMinimizer):
254              else:              else:
255                  raise KeyError("Invalid option '%s'"%o)                  raise KeyError("Invalid option '%s'"%o)
256    
     def getNorm(self, x):  
         return sqrt(self._f.getInner(x, x))  
   
257      def run(self, x):      def run(self, x):
258          args=self._f.getArguments(x)          args=self._f.getArguments(x)
259          gf=self._f.getGradient(x, *args)          gf=self._f.getGradient(x, *args)
260          fx=self._f(x, *args)          fx=self._f(x, *args)
261            if self._f.provides_inverse_Hessian_approximation:
262                self._f.updateHessian()
263                invH_scale = None
264            else:
265            invH_scale = self._initial_H
266          k=0          k=0
         H=self._initial_H  
267          error=2*self._tol          error=2*self._tol
268          s_and_y=[]          s_and_y=[]
269    
# Line 270  class MinimizerLBFGS(AbstractMinimizer): Line 273  class MinimizerLBFGS(AbstractMinimizer):
273              #self.logger.info("\033[1;31miteration %d\033[1;30m, error=%e"%(k,error))              #self.logger.info("\033[1;31miteration %d\033[1;30m, error=%e"%(k,error))
274              self.logger.info("iteration %d, error=%e"%(k,error))              self.logger.info("iteration %d, error=%e"%(k,error))
275              # determine search direction              # determine search direction
276              p = -self._twoLoop(H, gf, s_and_y)              p = -self._twoLoop(invH_scale, gf, s_and_y, x, *args)
277    
278              self.logger.debug("H = %s"%H)              if invH_scale: self.logger.debug("H = %s"%invH_scale)
279              self.logger.debug("grad f(x) = %s"%gf)              self.logger.debug("grad f(x) = %s"%gf)
280              self.logger.debug("p = %s"%p)              #self.logger.debug("p = %s"%p)
281              self.logger.debug("x = %s"%x)              #self.logger.debug("x = %s"%x)
282    
283              # determine step length              # determine step length
284              alpha, fx_new, gf_new = line_search(self._f, x, p, gf, fx)              alpha, fx_new, gf_new = line_search(self._f, x, p, gf, fx)
# Line 283  class MinimizerLBFGS(AbstractMinimizer): Line 286  class MinimizerLBFGS(AbstractMinimizer):
286              # execute the step              # execute the step
287              x_new = x + alpha*p              x_new = x + alpha*p
288              if gf_new is None:              if gf_new is None:
289                  gf_new=self._f.getGradient(x_new)              args=self._f.getArguments(x_new)
290                    gf_new=self._f.getGradient(x_new, args)
291    
292              delta_x=x_new-x              delta_x=x_new-x
293              delta_g=gf_new-gf              delta_g=gf_new-gf
294              s_and_y.append((delta_x,delta_g))              s_and_y.append((delta_x,delta_g, self._f.getDualProduct(delta_x, delta_g) ))
295                
296                # this needs to be reviewed
297              if fx_new==0.:              if fx_new==0.:
298                  error=fx                  error=fx
299              else:              else:
300                  error=abs(fx_new-fx)/abs(fx_new)                  error=abs(fx_new-fx)/abs(fx_new)
301                
302                self._f.updateHessian()
303              x=x_new              x=x_new
304              gf=gf_new              gf=gf_new
305              fx=fx_new              fx=fx_new
# Line 301  class MinimizerLBFGS(AbstractMinimizer): Line 310  class MinimizerLBFGS(AbstractMinimizer):
310              # delete oldest vector pair              # delete oldest vector pair
311              if k>self._m: s_and_y.pop(0)              if k>self._m: s_and_y.pop(0)
312    
313              # set the new scaling factor (approximation of inverse Hessian)              if not self._f.provides_inverse_Hessian_approximation:            
314              gnorm=self.getNorm(gf)            # set the new scaling factor (approximation of inverse Hessian)
315              denom=self.getNorm(delta_g)            denom=self._f.getDualProduct(delta_g, delta_g)
316              if denom < EPSILON * gnorm:            if denom > 0:
317                  H=1.                invH_scale=self._f.getDualProduct(delta_x,delta_g)/denom
318                  self.logger.debug("Break down in H update. Resetting to 1.")            else:
319              else:                invH_scale=self._initial_H
320                  H=self._f.getInner(delta_x,delta_g)/denom**2                self.logger.debug("Break down in H update. Resetting to initial value %s."%self._initial_H)
321    
322          if k >= self._imax:          if k >= self._imax:
323              reason=self.MAX_ITERATIONS_REACHED              reason=self.MAX_ITERATIONS_REACHED
# Line 319  class MinimizerLBFGS(AbstractMinimizer): Line 328  class MinimizerLBFGS(AbstractMinimizer):
328          self._result=x          self._result=x
329          return reason          return reason
330    
331      def _twoLoop(self, H, gf, s_and_y):      def _twoLoop(self, invH_scale, gf, s_and_y, x, *args):
332          """          """
333          Helper for the L-BFGS method.          Helper for the L-BFGS method.
334          See 'Numerical Optimization' by J. Nocedal for an explanation.          See 'Numerical Optimization' by J. Nocedal for an explanation.
335          """          """
336          q=gf          q=gf
337          alpha=[]          alpha=[]
338          for s,y in reversed(s_and_y):          for s,y, rho in reversed(s_and_y):
339              rho=1./(self._f.getInner(s, y))              a=self._f.getDualProduct(s, q)/rho
             a=rho*self._f.getInner(s, q)  
340              alpha.append(a)              alpha.append(a)
341              q=q-a*y              q=q-a*y
342    
343          r=H*q          if self._f.provides_inverse_Hessian_approximation:    
344          for s,y in s_and_y:               r= self._f.getInverseHessianApproximation(x, q, *args)
345              rho=1./(self._f.getInner(s, y))          else:
346              beta=rho*self._f.getInner(y,r)           r= invH_scale * q
347            
348            for s,y,rho in s_and_y:
349                beta=self._f.getDualProduct(r, y)/rho
350              a=alpha.pop()              a=alpha.pop()
351              r=r+s*(a-beta)              r=r+s*(a-beta)
352          return r          return r
# Line 377  class MinimizerBFGS(AbstractMinimizer): Line 388  class MinimizerBFGS(AbstractMinimizer):
388              self.logger.info("iteration %d, gnorm=%e"%(k,gnorm))              self.logger.info("iteration %d, gnorm=%e"%(k,gnorm))
389    
390              # determine search direction              # determine search direction
391              d=-self._f.getInner(H, gf)              d=-self._f.getDualProduct(H, gf)
392    
393              self.logger.debug("H = %s"%H)              self.logger.debug("H = %s"%H)
394              self.logger.debug("grad f(x) = %s"%gf)              self.logger.debug("grad f(x) = %s"%gf)
# Line 401  class MinimizerBFGS(AbstractMinimizer): Line 412  class MinimizerBFGS(AbstractMinimizer):
412              if (gnorm<=self._tol): break              if (gnorm<=self._tol): break
413    
414              # update Hessian              # update Hessian
415              denom=self._f.getInner(delta_x, delta_g)              denom=self._f.getDualProduct(delta_x, delta_g)
416              if denom < EPSILON * gnorm:              if denom < EPSILON * gnorm:
417                  denom=1e-5                  denom=1e-5
418                  self.logger.debug("Break down in H update. Resetting.")                  self.logger.debug("Break down in H update. Resetting.")
# Line 409  class MinimizerBFGS(AbstractMinimizer): Line 420  class MinimizerBFGS(AbstractMinimizer):
420              self.logger.debug("rho=%e"%rho)              self.logger.debug("rho=%e"%rho)
421              A=I-rho*delta_x[:,None]*delta_g[None,:]              A=I-rho*delta_x[:,None]*delta_g[None,:]
422              AT=I-rho*delta_g[:,None]*delta_x[None,:]              AT=I-rho*delta_g[:,None]*delta_x[None,:]
423              H=self._f.getInner(A, self._f.getInner(H,AT)) + rho*delta_x[:,None]*delta_x[None,:]              H=self._f.getDualProduct(A, self._f.getDualProduct(H,AT)) + rho*delta_x[:,None]*delta_x[None,:]
424          if k >= self._imax:          if k >= self._imax:
425              reason=self.MAX_ITERATIONS_REACHED              reason=self.MAX_ITERATIONS_REACHED
426              self.logger.warning("Maximum number of iterations reached!")              self.logger.warning("Maximum number of iterations reached!")
# Line 434  class MinimizerNLCG(AbstractMinimizer): Line 445  class MinimizerNLCG(AbstractMinimizer):
445          r=-self._f.getGradient(x, *args)          r=-self._f.getGradient(x, *args)
446          fx=self._f(x, *args)          fx=self._f(x, *args)
447          d=r          d=r
448          delta=self._f.getInner(r,r)          delta=self._f.getDualProduct(r,r)
449          delta0=delta          delta0=delta
450          self._doCallback(i, x, fx, -r)          self._doCallback(i, x, fx, -r)
451    
# Line 449  class MinimizerNLCG(AbstractMinimizer): Line 460  class MinimizerNLCG(AbstractMinimizer):
460              x=x+alpha*d              x=x+alpha*d
461              r=-self._f.getGradient(x) if gf_new is None else -gf_new              r=-self._f.getGradient(x) if gf_new is None else -gf_new
462              delta_o=delta              delta_o=delta
463              delta=self._f.getInner(r,r)              delta=self._f.getDualProduct(r,r)
464              beta=delta/delta_o              beta=delta/delta_o
465              d=r+beta*d              d=r+beta*d
466              k=k+1              k=k+1
# Line 457  class MinimizerNLCG(AbstractMinimizer): Line 468  class MinimizerNLCG(AbstractMinimizer):
468                  lenx=len(x)                  lenx=len(x)
469              except:              except:
470                  lenx=x.getNumberOfDataPoints()                  lenx=x.getNumberOfDataPoints()
471              if k == lenx or self._f.getInner(r,d) <= 0:              if k == lenx or self._f.getDualProduct(r,d) <= 0:
472                  d=r                  d=r
473                  k=0                  k=0
474              i+=1              i+=1
# Line 477  class MinimizerNLCG(AbstractMinimizer): Line 488  class MinimizerNLCG(AbstractMinimizer):
488  if __name__=="__main__":  if __name__=="__main__":
489      # Example usage with function 'rosen' (minimum=[1,1,...1]):      # Example usage with function 'rosen' (minimum=[1,1,...1]):
490      from scipy.optimize import rosen, rosen_der      from scipy.optimize import rosen, rosen_der
491      from esys.downunder.costfunctions import CostFunction      from esys.downunder import  MeteredCostFunction
492      import sys      import sys
493      N=100      N=10
494      x0=np.array([4.]*N) # initial guess      x0=np.array([4.]*N) # initial guess
495    
496      class RosenFunc(CostFunction):      class RosenFunc(MeteredCostFunction):
497          def _getInner(self, f0, f1):          def __init__(self):
498           super(RosenFunc, self).__init__()
499               self.provides_inverse_Hessian_approximation=False
500            def _getDualProduct(self, f0, f1):
501              return np.dot(f0, f1)              return np.dot(f0, f1)
502          def _getValue(self, x, *args):          def _getValue(self, x, *args):
503              return rosen(x)              return rosen(x)
# Line 504  if __name__=="__main__": Line 518  if __name__=="__main__":
518          m=MinimizerLBFGS(f)          m=MinimizerLBFGS(f)
519          #m.setOptions(historySize=10000)          #m.setOptions(historySize=10000)
520    
521      logging.basicConfig(format='[%(funcName)s] \033[1;30m%(message)s\033[0m', level=logging.INFO)      logging.basicConfig(format='[%(funcName)s] \033[1;30m%(message)s\033[0m', level=logging.DEBUG)
522      m.setTolerance(1e-5)      m.setTolerance(1e-5)
523      m.setMaxIterations(600)      m.setMaxIterations(600)
524      m.run(x0)      m.run(x0)

Legend:
Removed from v.4073  
changed lines
  Added in v.4074

  ViewVC Help
Powered by ViewVC 1.1.26