1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
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
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
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
84 self.f_count = 0
85 self.g_count = 0
86 self.h_count = 0
87
88
89 self.warning = None
90
91
92 self.setup_conv_tests()
93
94
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
101 """Find the Cauchy point."""
102
103
104 curv = dot(self.dfk, dot(self.d2fk, self.dfk))
105 norm_dfk = sqrt(dot(self.dfk, self.dfk))
106
107
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
118 self.pk = - self.tau_k * self.delta * self.dfk / norm_dfk
119
120
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
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