Package minimise :: Module cauchy_point
[hide private]
[frames] | no frames]

Source Code for Module minimise.cauchy_point

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2003 Edward d'Auvergne                                        # 
  4  #                                                                             # 
  5  # This file is part of the program relax.                                     # 
  6  #                                                                             # 
  7  # relax is free software; you can redistribute it and/or modify               # 
  8  # it under the terms of the GNU General Public License as published by        # 
  9  # the Free Software Foundation; either version 2 of the License, or           # 
 10  # (at your option) any later version.                                         # 
 11  #                                                                             # 
 12  # relax is distributed in the hope that it will be useful,                    # 
 13  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 14  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 15  # GNU General Public License for more details.                                # 
 16  #                                                                             # 
 17  # You should have received a copy of the GNU General Public License           # 
 18  # along with relax; if not, write to the Free Software                        # 
 19  # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA   # 
 20  #                                                                             # 
 21  ############################################################################### 
 22   
 23   
 24  from Numeric import dot, sqrt 
 25   
 26  from base_classes import Trust_region, Min 
 27   
 28   
29 -def cauchy_point(func=None, dfunc=None, d2func=None, args=(), x0=None, func_tol=1e-25, grad_tol=None, maxiter=1e6, delta_max=1e5, delta0=1.0, eta=0.2, full_output=0, print_flag=0, print_prefix=""):
30 """Cauchy Point trust region algorithm. 31 32 Page 69 from 'Numerical Optimization' by Jorge Nocedal and Stephen J. Wright, 1999, 2nd ed. 33 The Cauchy point is defined by: 34 35 delta 36 pCk = - tau_k ------- dfk 37 ||dfk|| 38 39 where: 40 delta_k is the trust region radius, 41 dfk is the gradient vector, 42 43 / 1 if dfk . Bk . dfk <= 0 44 tau_k = < 45 \ min(||dfk||**2/(delta . dfk . Bk . dfk), 1) otherwise. 46 """ 47 48 if print_flag: 49 if print_flag >= 2: 50 print print_prefix 51 print print_prefix 52 print print_prefix + "Cauchy point minimisation" 53 print print_prefix + "~~~~~~~~~~~~~~~~~~~~~~~~~" 54 min = Cauchy_point(func, dfunc, d2func, args, x0, func_tol, grad_tol, maxiter, delta_max, delta0, eta, full_output, print_flag, print_prefix) 55 results = min.minimise() 56 return results
57 58
59 -class Cauchy_point(Trust_region, Min):
60 - def __init__(self, func, dfunc, d2func, args, x0, func_tol, grad_tol, maxiter, delta_max, delta0, eta, full_output, print_flag, print_prefix):
61 """Class for Cauchy Point trust region minimisation specific functions. 62 63 Unless you know what you are doing, you should call the function 'cauchy_point' rather than 64 using this class. 65 """ 66 67 # Function arguments. 68 self.func = func 69 self.dfunc = dfunc 70 self.d2func = d2func 71 self.args = args 72 self.xk = x0 73 self.func_tol = func_tol 74 self.grad_tol = grad_tol 75 self.maxiter = maxiter 76 self.full_output = full_output 77 self.print_flag = print_flag 78 self.print_prefix = print_prefix 79 self.delta_max = delta_max 80 self.delta = delta0 81 self.eta = eta 82 83 # Initialise the function, gradient, and Hessian evaluation counters. 84 self.f_count = 0 85 self.g_count = 0 86 self.h_count = 0 87 88 # Initialise the warning string. 89 self.warning = None 90 91 # Set the convergence test function. 92 self.setup_conv_tests() 93 94 # Initial values before the first iteration. 95 self.fk, self.f_count = apply(self.func, (self.xk,)+self.args), self.f_count + 1 96 self.dfk, self.g_count = apply(self.dfunc, (self.xk,)+self.args), self.g_count + 1 97 self.d2fk, self.h_count = apply(self.d2func, (self.xk,)+self.args), self.h_count + 1
98 99
100 - def new_param_func(self):
101 """Find the Cauchy point.""" 102 103 # Calculate the curvature and norm. 104 curv = dot(self.dfk, dot(self.d2fk, self.dfk)) 105 norm_dfk = sqrt(dot(self.dfk, self.dfk)) 106 107 # tau_k. 108 if curv <= 0.0: 109 self.tau_k = 1.0 110 else: 111 self.tau_k = min(norm_dfk ** 3 / (self.delta * curv), 1.0) 112 113 if self.print_flag >= 2: 114 print self.print_prefix + "dfk . Bk . dfk: " + `curv` 115 print self.print_prefix + "tau_k: " + `self.tau_k` 116 117 # Cauchy point. 118 self.pk = - self.tau_k * self.delta * self.dfk / norm_dfk 119 120 # Find the new parameter vector and function value at that point. 121 self.xk_new = self.xk + self.pk 122 self.fk_new, self.f_count = apply(self.func, (self.xk_new,)+self.args), self.f_count + 1 123 self.dfk_new, self.g_count = apply(self.dfunc, (self.xk_new,)+self.args), self.g_count + 1
124 125
126 - def update(self):
127 """Update function. 128 129 Function to update the function value, gradient vector, and Hessian matrix. 130 """ 131 132 self.xk = self.xk_new * 1.0 133 self.fk = self.fk_new 134 self.dfk = self.dfk_new * 1.0 135 self.d2fk, self.h_count = apply(self.d2func, (self.xk,)+self.args), self.h_count + 1
136