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 LinearAlgebra import inverse
25 from Numeric import Float64, dot, matrixmultiply, sqrt, zeros
26
27 from base_classes import Line_search, Min
28
29
30 -def ncg(func=None, dfunc=None, d2func=None, args=(), x0=None, min_options=None, func_tol=1e-25, grad_tol=None, maxiter=1e6, a0=1.0, mu=0.0001, eta=0.9, full_output=0, print_flag=0, print_prefix=""):
31 """Line search Newton conjugate gradient algorithm.
32
33 Page 140 from 'Numerical Optimization' by Jorge Nocedal and Stephen J. Wright, 1999, 2nd ed.
34
35 The algorithm is:
36
37 Given initial point x0.
38 while 1:
39 Compute a search direction pk by applying the CG method to Hk.pk = -gk, starting from
40 x0 = 0. Terminate when ||rk|| <= min(0.5,sqrt(||gk||)), or if negative curvature is
41 encountered.
42 Set xk+1 = xk + ak.pk, where ak satisfies the Wolfe, Goldstein, or Armijo backtracking
43 conditions.
44 """
45
46 if print_flag:
47 if print_flag >= 2:
48 print print_prefix
49 print print_prefix
50 print print_prefix + "Newton Conjugate Gradient minimisation"
51 print print_prefix + "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
52 min = Ncg(func, dfunc, d2func, args, x0, min_options, func_tol, grad_tol, maxiter, a0, mu, eta, full_output, print_flag, print_prefix)
53 if min.init_failure:
54 print print_prefix + "Initialisation of minimisation has failed."
55 return None
56 results = min.minimise()
57 return results
58
59
60 -class Ncg(Line_search, Min):
61 - def __init__(self, func, dfunc, d2func, args, x0, min_options, func_tol, grad_tol, maxiter, a0, mu, eta, full_output, print_flag, print_prefix):
62 """Class for newton conjugate gradient minimisation specific functions.
63
64 Unless you know what you are doing, you should call the function 'ncg' rather than using
65 this class.
66 """
67
68
69 self.func = func
70 self.dfunc = dfunc
71 self.d2func = d2func
72 self.args = args
73 self.xk = x0
74 self.func_tol = func_tol
75 self.grad_tol = grad_tol
76 self.maxiter = maxiter
77 self.full_output = full_output
78 self.print_flag = print_flag
79 self.print_prefix = print_prefix
80
81
82 self.a0 = a0
83
84
85 self.mu = mu
86 self.eta = eta
87
88
89 self.init_failure = 0
90
91
92 self.line_search_options(min_options)
93 self.setup_line_search()
94
95
96 self.f_count = 0
97 self.g_count = 0
98 self.h_count = 0
99
100
101 self.warning = None
102
103
104 self.setup_conv_tests()
105
106
107 self.fk, self.f_count = apply(self.func, (self.xk,)+self.args), self.f_count + 1
108 self.dfk, self.g_count = apply(self.dfunc, (self.xk,)+self.args), self.g_count + 1
109 self.d2fk, self.h_count = apply(self.d2func, (self.xk,)+self.args), self.h_count + 1
110
111
113 """The CG algorithm."""
114
115
116 xi = zeros(len(self.xk), Float64)
117 ri = self.dfk * 1.0
118 pi = -ri
119 dot_ri = dot(ri, ri)
120 len_g = sqrt(dot_ri)
121 residual_test = min(0.5 * len_g, dot_ri)
122
123
124
125 if self.print_flag >= 2:
126 print self.print_prefix + "Initial data:"
127 print self.print_prefix + "\tx0: " + `xi`
128 print self.print_prefix + "\tr0: " + `ri`
129 print self.print_prefix + "\tp0: " + `pi`
130
131 i = 0
132 while 1:
133
134 Api = dot(self.d2fk, pi)
135 curv = dot(pi, Api)
136
137
138 if curv <= 0.0:
139 if i == 0:
140 if dot_ri == 0.0:
141 return xi
142 else:
143 ai = dot_ri / curv
144 return xi + ai*pi
145 else:
146 return xi
147 if sqrt(dot_ri) <= residual_test:
148 return xi
149
150
151 ai = dot_ri / curv
152 xi_new = xi + ai*pi
153 ri_new = ri + ai*Api
154 dot_ri_new = dot(ri_new, ri_new)
155 bi_new = dot_ri_new / dot_ri
156 pi_new = -ri_new + bi_new*pi
157
158
159 if self.print_flag >= 2:
160 print ""
161 print self.print_prefix + "Iteration i = " + `i`
162 print self.print_prefix + "Api: " + `Api`
163 print self.print_prefix + "Curv: " + `curv`
164 print self.print_prefix + "ai: " + `ai`
165 print self.print_prefix + "xi+1: " + `xi_new`
166 print self.print_prefix + "ri+1: " + `ri_new`
167 print self.print_prefix + "bi+1: " + `bi_new`
168 print self.print_prefix + "pi+1: " + `pi_new`
169
170
171 xi = xi_new * 1.0
172 ri = ri_new * 1.0
173 pi = pi_new * 1.0
174 dot_ri = dot_ri_new
175 i = i + 1
176
177
179 """The new parameter function.
180
181 Find the search direction, do a line search, and get xk+1 and fk+1.
182 """
183
184
185 self.pk = self.get_pk()
186
187
188 self.line_search()
189
190
191 self.xk_new = self.xk + self.alpha * self.pk
192 self.fk_new, self.f_count = apply(self.func, (self.xk_new,)+self.args), self.f_count + 1
193 self.dfk_new, self.g_count = apply(self.dfunc, (self.xk_new,)+self.args), self.g_count + 1
194
195
197 """Function to update the function value, gradient vector, and Hessian matrix."""
198
199 self.xk = self.xk_new * 1.0
200 self.fk = self.fk_new
201 self.dfk = self.dfk_new * 1.0
202 self.d2fk, self.h_count = apply(self.d2func, (self.xk,)+self.args), self.h_count + 1
203