3 Local search - A trust region algorithm with BFGS update.
9 #include "stogo_config.h"
15 # define IF_NLOPT_CHECK_EVALS stop->nevals++; \
16 if (nlopt_stop_evalstime(stop)) \
19 # define IF_NLOPT_CHECK_EVALS
22 int local(Trial &T, TBox &box, TBox &domain, double eps_cl, double *mgr,
23 Global &glob, int axis, RCRVector x_av
25 , nlopt_stopping *stop
29 int k_max, info, outside ;
30 int k, i, good_enough, iTmp ;
33 double maxgrad, delta, f, f_new, tmp ;
34 double alpha, gamma, beta, d2, s2, nom, den, ro ;
35 double nrm_sd, nrm_hn, snrm_hn, nrm_dl ;
36 RVector x(n), g(n), h_sd(n), h_dl(n), h_n(n), x_new(n), g_new(n) ;
37 RVector s(n),y(n),z(n),w(n) ; // Temporary vectors
38 RMatrix B(n), H(n) ; // Hessian and it's inverse
44 cout << "Local Search, x=" << x << endl;
47 if (box.OutsideBox(x, domain) != 0) {
48 cout << "Starting point is not inside the boundary. Exiting...\n" ;
53 // Check if we are close to a stationary point located previously
54 if (box.CloseToMin(x, &tmp, eps_cl)) {
56 cout << "Close to a previously located stationary point, exiting" << endl;
62 // Initially B and H are equal to the identity matrix
64 for (i=0 ; i<n ; i++) {
69 RVector g_av(x_av.GetLength());
71 f=glob.ObjectiveGradient(x,g,OBJECTIVE_AND_GRADIENT);
75 f=glob.ObjectiveGradient(x_av,g_av,OBJECTIVE_AND_GRADIENT);
84 // Elaborate scheme to initalize delta
85 delta=delta_coef*norm2(g) ;
88 if (!box.InsideBox(z)) {
89 if (box.Intersection(x,g,z)==TRUE) {
91 delta=min(delta,delta_coef*norm2(z)) ;
94 // Algorithm broke down, use INI1
95 delta = (1.0/7)*box.ShortestSide(&iTmp) ;
101 delta = box.ClosestSide(x)*delta_coef ;
102 if (delta<MacEpsilon)
103 // Patch to avoid trust region with radius close to zero
104 delta = (1.0/7)*box.ShortestSide(&iTmp) ;
107 delta = delta_coef*box.ShortestSide(&iTmp) ;
111 // Use a simple scheme for the 1D minimization (INI1)
112 delta = (1.0/7.0)*box.ShortestSide(&iTmp) ;
115 k=0 ; good_enough = 0 ; info=LS_New ; outside=0 ;
117 while (good_enough == 0) {
121 cout << "Maximum number of iterations reached\n" ;
127 // Update maximal gradient value
128 maxgrad=max(maxgrad,normInf(g)) ;
130 // Steepest descent, h_sd = -g
135 if (nrm_sd < epsilon) {
136 // Stop criterion (gradient) fullfilled
138 cout << "Gradient small enough" << endl ;
144 // Compute Newton step, h_n = -H*g
145 gemv('N',-1.0, H, g, 0.0, h_n) ;
146 nrm_hn = norm2(h_n) ;
148 if (nrm_hn < delta) {
152 cout << "[Newton step] " ;
156 gemv('N',1.0,B,g,0.0,z) ;
162 alpha=(nrm_sd*nrm_sd)/tmp ; // Normalization (N38,eq. 3.30)
164 nrm_sd=fabs(alpha)*nrm_sd ;
166 if (nrm_sd >= delta) {
167 gamma = delta/nrm_sd ; // Normalization (N38, eq. 3.33)
171 cout << "[Steepest descent] " ;
175 // Combination of Newton and SD steps
180 snrm_hn=nrm_hn*nrm_hn ;
182 den = tmp-s2 + sqrt((tmp-d2)*(tmp-d2)+(snrm_hn-d2)*(d2-s2)) ;
187 // Normalization (N38, eq. 3.31)
191 axpy((1-beta),h_sd,h_dl) ;
193 cout << "[Mixed step] " ;
201 axpy(1.0,h_dl,x_new) ;
203 // Check if x_new is inside the box
204 iTmp=box.OutsideBox(x_new, domain) ;
207 cout << "x_new is outside the box " << endl ;
210 if (outside>max_outside_steps) {
211 // Previous point was also outside, exit
215 else if (iTmp == 2) {
217 cout << " x_new is outside the domain" << endl ;
228 f_new=glob.ObjectiveGradient(x_new,g_new,OBJECTIVE_AND_GRADIENT);
231 f_new=glob.ObjectiveGradient(x_av,g_av,OBJECTIVE_AND_GRADIENT);
233 IF_NLOPT_CHECK_EVALS;
235 gemv('N',0.5,B,h_dl,0.0,z);
236 ro = (f_new-f) / (dot(g,h_dl) + dot(h_dl,z)); // Quadratic model
244 // Update the Hessian and it's inverse using the BFGS formula
245 #if 0 // changed by SGJ to compute OBJECTIVE_AND_GRADIENT above
247 glob.ObjectiveGradient(x_new,g_new,GRADIENT_ONLY);
250 glob.ObjectiveGradient(x_av,g_av,GRADIENT_ONLY);
254 IF_NLOPT_CHECK_EVALS;
264 // Check curvature condition
266 if (alpha <= sqrt(MacEpsilon)*nrm_dl*norm2(y)) {
268 cout << "Curvature condition violated " ;
273 gemv('N',1.0,B,h_dl,0.0,z) ; // z=Bh_dl
274 beta=-1/dot(h_dl,z) ;
278 // Update Hessian inverse
279 gemv('N',1.0,H,y,0.0,z) ; // z=H*y
280 gemv('T',1.0,H,y,0.0,w) ; // w=y'*H
282 beta=(1+beta/alpha)/alpha ;
284 // It should be possible to do this updating more efficiently, by
285 // exploiting the fact that (h_dl*y'*H) = transpose(H*y*h_dl')
286 ger(beta,h_dl,h_dl,H) ;
287 ger(-1/alpha,z,h_dl,H) ;
288 ger(-1/alpha,h_dl,w,H) ;
291 if (nrm_dl < norm2(x)*epsilon) {
292 // Stop criterion (iteration progress) fullfilled
294 cout << "Progress is marginal" ;
299 // Check if we are close to a stationary point located previously
300 if (box.CloseToMin(x_new, &f_new, eps_cl)) {
301 // Note that x_new and f_new may be overwritten on exit from CloseToMin
303 cout << "Close to a previously located stationary point, exiting" << endl;
310 copy(x_new,x) ; copy(g_new,g) ; f=f_new ;
313 cout << " x=" << x << endl ;
319 cout << "Step is no good, ro=" << ro << " delta=" << delta << endl ;
325 // Make sure the routine returns correctly...
326 // Check if last iterate is outside the boundary
327 if (box.OutsideBox(x, domain) != 0) {
328 info=LS_Out; f=DBL_MAX;
331 if (info == LS_Unstable) {
332 cout << "Local search became unstable. No big deal but exiting anyway\n" ;
336 T.xvals=x ; T.objval=f ;