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

Source Code for Module minimise.newton

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2003, 2005 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 LinearAlgebra import cholesky_decomposition, eigenvectors, inverse, solve_linear_equations 
 25  from math import acos 
 26  from Numeric import Float64, array, dot, identity, matrixmultiply, sort, sqrt, trace, transpose, zeros 
 27  from re import match 
 28   
 29  from base_classes import Hessian_mods, Line_search, Min 
 30   
 31   
32 -def newton(func=None, dfunc=None, d2func=None, args=(), x0=None, min_options=(), func_tol=1e-25, grad_tol=None, maxiter=1e6, a0=1.0, mu=0.0001, eta=0.9, mach_acc=1e-16, full_output=0, print_flag=0, print_prefix=""):
33 """Newton minimisation.""" 34 35 if print_flag: 36 if print_flag >= 2: 37 print print_prefix 38 print print_prefix 39 print print_prefix + "Newton minimisation" 40 print print_prefix + "~~~~~~~~~~~~~~~~~~~" 41 min = Newton(func, dfunc, d2func, args, x0, min_options, func_tol, grad_tol, maxiter, a0, mu, eta, mach_acc, full_output, print_flag, print_prefix) 42 if min.init_failure: 43 print print_prefix + "Initialisation of minimisation has failed." 44 return None 45 results = min.minimise() 46 return results
47 48
49 -class Newton(Hessian_mods, Line_search, Min):
50 - def __init__(self, func, dfunc, d2func, args, x0, min_options, func_tol, grad_tol, maxiter, a0, mu, eta, mach_acc, full_output, print_flag, print_prefix):
51 """Class for Newton minimisation specific functions. 52 53 Unless you know what you are doing, you should call the function 'newton' rather than using 54 this class. 55 """ 56 57 # Function arguments. 58 self.func = func 59 self.dfunc = dfunc 60 self.d2func = d2func 61 self.args = args 62 self.xk = x0 63 self.func_tol = func_tol 64 self.grad_tol = grad_tol 65 self.maxiter = maxiter 66 self.mach_acc = mach_acc 67 self.full_output = full_output 68 self.print_flag = print_flag 69 self.print_prefix = print_prefix 70 71 # Set a0. 72 self.a0 = a0 73 74 # Line search constants for the Wolfe conditions. 75 self.mu = mu 76 self.eta = eta 77 78 # Initialisation failure flag. 79 self.init_failure = 0 80 81 # Setup the line search and Hessian modification algorithms. 82 self.line_search_algor = None 83 self.hessian_mod = None 84 85 # Test if the options are a tuple. 86 if type(min_options) != tuple: 87 print self.print_prefix + "The minimisation options " + `min_options` + " is not a tuple." 88 self.init_failure = 1; return 89 90 # Test that no more thant 2 options are given. 91 if len(min_options) > 2: 92 print self.print_prefix + "A maximum of two minimisation options is allowed (the line search algorithm and the Hessian modification)." 93 self.init_failure = 1; return 94 95 # Sort out the minimisation options. 96 none_option = None 97 for opt in min_options: 98 if match('[Nn]one', opt): 99 none_option = opt 100 continue 101 if self.line_search_algor == None and self.valid_line_search(opt): 102 self.line_search_algor = opt 103 elif self.hessian_mod == None and self.valid_hessian_mod(opt): 104 self.hessian_mod = opt 105 else: 106 if self.line_search_algor: 107 print self.print_prefix + "The minimisation option " + `opt` + " from " + `min_options` + " is not a valid Hessian modification." 108 elif self.hessian_mod: 109 print self.print_prefix + "The minimisation option " + `opt` + " from " + `min_options` + " is not a valid line search algorithm." 110 else: 111 print self.print_prefix + "The minimisation option " + `opt` + " from " + `min_options` + " is neither a valid line search algorithm or Hessian modification." 112 self.init_failure = 1; return 113 if none_option and not self.hessian_mod: 114 if self.valid_hessian_mod(none_option): 115 self.hessian_mod = none_option 116 elif none_option and not self.line_search_algor: 117 if self.valid_line_search(none_option): 118 self.line_search_algor = none_option 119 elif none_option: 120 print self.print_prefix + "Invalid combination of options " + `min_options` + "." 121 self.init_failure = 1; return 122 123 # Default line search algorithm. 124 if self.line_search_algor == None: 125 self.line_search_algor = 'Back' 126 127 # Default Hessian modification. 128 if self.hessian_mod == None: 129 self.hessian_mod = 'GMW' 130 131 # Line search and Hessian modification initialisation. 132 self.setup_line_search() 133 self.setup_hessian_mod() 134 135 # Initialise the function, gradient, and Hessian evaluation counters. 136 self.f_count = 0 137 self.g_count = 0 138 self.h_count = 0 139 140 # Initialise the warning string. 141 self.warning = None 142 143 # Constants. 144 self.n = len(self.xk) 145 self.I = identity(len(self.xk)) 146 147 # Set the convergence test function. 148 self.setup_conv_tests() 149 150 # Newton setup function. 151 self.setup_newton() 152 153 # Set the update function. 154 self.update = self.update_newton
155 156
157 - def new_param_func(self):
158 """The new parameter function. 159 160 Find the search direction, do a line search, and get xk+1 and fk+1. 161 """ 162 163 # Calculate the Newton direction. 164 self.pk = self.get_pk() 165 166 # Line search. 167 self.line_search() 168 169 # Find the new parameter vector and function value at that point. 170 self.xk_new = self.xk + self.alpha * self.pk 171 self.fk_new, self.f_count = apply(self.func, (self.xk_new,)+self.args), self.f_count + 1 172 self.dfk_new, self.g_count = apply(self.dfunc, (self.xk_new,)+self.args), self.g_count + 1 173 174 # Debugging. 175 if self.print_flag >= 2: 176 print "\n" + self.print_prefix + "New param function." 177 print self.print_prefix + "pk: " + `self.pk` 178 print self.print_prefix + "alpha: " + `self.alpha` 179 print self.print_prefix + "xk: " + `self.xk` 180 print self.print_prefix + "xk+1: " + `self.xk_new` 181 print self.print_prefix + "fk: " + `self.fk` 182 print self.print_prefix + "fk+1: " + `self.fk_new` 183 print self.print_prefix + "dfk: " + `self.dfk` 184 print self.print_prefix + "dfk+1: " + `self.dfk_new` 185 print self.print_prefix + "d2fk:\n" + `self.d2fk` 186 187 eigen = eigenvectors(self.d2fk) 188 print self.print_prefix + "Eigenvalues: " + `eigen[0]` 189 190 print self.print_prefix + "Angle to the unit vector pointing along the first dimension." 191 unit_vect = zeros(self.n, Float64) 192 unit_vect[0] = 1.0 193 dfk_angle = acos(dot(self.dfk, unit_vect) / sqrt(dot(self.dfk, self.dfk))) 194 pk_angle = acos(dot(self.pk, unit_vect) / sqrt(dot(self.pk, self.pk))) 195 print self.print_prefix + "steepest descent: " + `dfk_angle` 196 print self.print_prefix + "Newton step: " + `pk_angle`
197 198
199 - def setup_newton(self):
200 """Setup function. 201 202 The initial Newton function value, gradient vector, and Hessian matrix are calculated. 203 """ 204 205 self.fk, self.f_count = apply(self.func, (self.xk,)+self.args), self.f_count + 1 206 self.dfk, self.g_count = apply(self.dfunc, (self.xk,)+self.args), self.g_count + 1 207 self.d2fk, self.h_count = apply(self.d2func, (self.xk,)+self.args), self.h_count + 1
208 209
210 - def update_newton(self):
211 """Function to update the function value, gradient vector, and Hessian matrix.""" 212 213 self.xk = self.xk_new * 1.0 214 self.fk = self.fk_new 215 self.dfk = self.dfk_new * 1.0 216 self.d2fk, self.h_count = apply(self.d2func, (self.xk,)+self.args), self.h_count + 1
217