stogo/local.cc

   1
   2 /*
   3    Local search - A trust region algorithm with BFGS update.
   4 */
   5
   6 #include <iostream>
   7 #include <stdlib.h>
   8
   9 #include "stogo_config.h"
  10 #include "global.h"
  11 #include "local.h"
  12 #include "tools.h"
  13
  14 int local(Trial &T, TBox &box, TBox &domain, double eps_cl, double *mgr,
  15           Global &glob, int axis, RCRVector x_av) {
  16
  17   int k_max, info, outside ;
  18   int k, i, good_enough, iTmp ;
  19   int n=box.GetDim();
  20
  21   double maxgrad, delta, f, f_new, tmp ;
  22   double alpha, gamma, beta, d2, s2, nom, den, ro ;
  23   double nrm_sd, nrm_hn, snrm_hn, nrm_dl ;
  24   RVector x(n), g(n), h_sd(n), h_dl(n), h_n(n), x_new(n), g_new(n) ;
  25   RVector s(n),y(n),z(n),w(n) ; // Temporary vectors
  26   RMatrix B(n), H(n) ;          // Hessian and it's inverse
  27
  28   k_max = max_iter*n ;
  29   x=T.xvals ;
  30
  31 #ifdef LS_DEBUG
  32   cout << "Local Search, x=" << x << endl;
  33 #endif
  34
  35   if (box.OutsideBox(x, domain) != 0) {
  36     cout << "Starting point is not inside the boundary. Exiting...\n" ;
  37     exit(1) ;
  38     return LS_Out ;
  39   }
  40
  41   // Check if we are close to a stationary point located previously
  42   if (box.CloseToMin(x, &tmp, eps_cl)) {
  43 #ifdef LS_DEBUG
  44      cout << "Close to a previously located stationary point, exiting" << endl;
  45 #endif
  46      T.objval=tmp;
  47      return LS_Old ;
  48    }
  49
  50   // Initially B and H are equal to the identity matrix
  51   B=0 ; H=0 ;
  52   for (i=0 ; i<n ; i++) {
  53     B(i,i)=1 ;
  54     H(i,i)=1 ;
  55   }
  56
  57   RVector g_av(x_av.GetLength());
  58   if (axis==-1) {
  59     f=glob.ObjectiveGradient(x,g,OBJECTIVE_AND_GRADIENT);
  60   }
  61   else {
  62     x_av(axis)=x(0);
  63     f=glob.ObjectiveGradient(x_av,g_av,OBJECTIVE_AND_GRADIENT);
  64     g(0)=g_av(axis);
  65   }
  66   FC++;GC++;
  67
  68   if (axis == -1) {
  69     // Skipping AV
  70 #ifdef INI3
  71     // Elaborate scheme to initalize delta
  72     delta=delta_coef*norm2(g) ;
  73     copy(g,z) ;
  74     axpy(1.0,x,z) ;
  75     if (!box.InsideBox(z)) {
  76       if (box.Intersection(x,g,z)==TRUE) {
  77         axpy(-1.0,x,z) ;
  78         delta=min(delta,delta_coef*norm2(z)) ;
  79       }
  80       else {
  81         // Algorithm broke down, use INI1
  82         delta = (1.0/7)*box.ShortestSide(&iTmp) ;
  83       }
  84     }
  85 #endif
  86 #ifdef INI2
  87     // Use INI2 scheme
  88     delta = box.ClosestSide(x)*delta_coef ;
  89     if (delta<MacEpsilon)
  90       // Patch to avoid trust region with radius close to zero
  91       delta = (1.0/7)*box.ShortestSide(&iTmp) ;
  92 #endif
  93 #ifdef INI1
  94     delta = delta_coef*box.ShortestSide(&iTmp) ;
  95 #endif
  96   }
  97   else {
  98     // Use a simple scheme for the 1D minimization (INI1)
  99     delta = (1.0/7.0)*box.ShortestSide(&iTmp) ;
 100   }
 101
 102   k=0 ; good_enough = 0 ; info=LS_New ; outside=0 ;
 103   maxgrad=*mgr ;
 104   while (good_enough == 0) {
 105     k++ ;
 106     if (k>k_max) {
 107 #ifdef LS_DEBUG
 108       cout << "Maximum number of iterations reached\n" ;
 109 #endif
 110       info=LS_MaxIter ;
 111       break ;
 112     }
 113
 114     // Update maximal gradient value
 115     maxgrad=max(maxgrad,normInf(g)) ;
 116
 117     // Steepest descent, h_sd = -g
 118     copy(g,h_sd) ;
 119     scal(-1.0,h_sd) ;
 120     nrm_sd=norm2(h_sd) ;
 121
 122     if (nrm_sd < epsilon) {
 123       // Stop criterion (gradient) fullfilled
 124 #ifdef LS_DEBUG
 125       cout << "Gradient small enough" << endl ;
 126 #endif
 127       good_enough = 1 ;
 128       break ;
 129     }
 130
 131     // Compute Newton step, h_n = -H*g
 132     gemv('N',-1.0, H, g, 0.0, h_n) ;
 133     nrm_hn = norm2(h_n) ;
 134
 135     if (nrm_hn < delta) {
 136       // Pure Newton step
 137       copy(h_n, h_dl) ;
 138 #ifdef LS_DEBUG
 139       cout << "[Newton step]      " ;
 140 #endif
 141     }
 142     else {
 143       gemv('N',1.0,B,g,0.0,z) ;
 144       tmp=dot(g,z) ;
 145       if (tmp==0) {
 146         info = LS_Unstable ;
 147         break ;
 148       }
 149       alpha=(nrm_sd*nrm_sd)/tmp ; // Normalization (N38,eq. 3.30)
 150       scal(alpha,h_sd) ;
 151       nrm_sd=fabs(alpha)*nrm_sd ;
 152
 153       if (nrm_sd >= delta) {
 154         gamma = delta/nrm_sd ; // Normalization (N38, eq. 3.33)
 155         copy(h_sd,h_dl) ;
 156         scal(gamma,h_dl) ;
 157 #ifdef LS_DEBUG
 158         cout << "[Steepest descent]  " ;
 159 #endif
 160       }
 161       else {
 162         // Combination of Newton and SD steps
 163         d2 = delta*delta ;
 164         copy(h_sd,s) ;
 165         s2=nrm_sd*nrm_sd ;
 166         nom = d2 - s2 ;
 167         snrm_hn=nrm_hn*nrm_hn ;
 168         tmp = dot(h_n,s) ;
 169         den = tmp-s2 + sqrt((tmp-d2)*(tmp-d2)+(snrm_hn-d2)*(d2-s2)) ;
 170         if (den==0) {
 171           info = LS_Unstable ;
 172           break ;
 173         }
 174         // Normalization (N38, eq. 3.31)
 175         beta = nom/den ;
 176         copy(h_n,h_dl) ;
 177         scal(beta,h_dl) ;
 178         axpy((1-beta),h_sd,h_dl) ;
 179 #ifdef LS_DEBUG
 180         cout << "[Mixed step]        " ;
 181 #endif
 182       }
 183     }
 184     nrm_dl=norm2(h_dl) ;
 185
 186     //x_new = x+h_dl ;
 187     copy(x,x_new) ;
 188     axpy(1.0,h_dl,x_new) ;
 189
 190     // Check if x_new is inside the box
 191     iTmp=box.OutsideBox(x_new, domain) ;
 192     if (iTmp == 1) {
 193 #ifdef LS_DEBUG
 194       cout << "x_new is outside the box " << endl ;
 195 #endif
 196       outside++ ;
 197       if (outside>max_outside_steps) {
 198         // Previous point was also outside, exit
 199         break ;
 200       }
 201     }
 202     else if (iTmp == 2) {
 203 #ifdef LS_DEBUG
 204       cout << " x_new is outside the domain" << endl ;
 205 #endif
 206       info=LS_Out ;
 207       break ;
 208     }
 209     else {
 210       outside=0 ;
 211     }
 212
 213     // Compute the gain
 214     if (axis==-1)
 215       f_new=glob.ObjectiveGradient(x_new,x_new,OBJECTIVE_ONLY);
 216     else {
 217       x_av(axis)=x_new(0);
 218       f_new=glob.ObjectiveGradient(x_av,x_new,OBJECTIVE_ONLY);
 219     }
 220     FC++;
 221     gemv('N',0.5,B,h_dl,0.0,z);
 222     ro = (f_new-f) / (dot(g,h_dl) + dot(h_dl,z)); // Quadratic model
 223     if (ro > 0.75) {
 224       delta = delta*2;
 225     }
 226     if (ro < 0.25) {
 227       delta = delta/3;
 228     }
 229     if (ro > 0) {
 230       // Update the Hessian and it's inverse using the BFGS formula
 231       if (axis==-1)
 232         glob.ObjectiveGradient(x_new,g_new,GRADIENT_ONLY);
 233       else {
 234         x_av(axis)=x_new(0);
 235         glob.ObjectiveGradient(x_av,g_av,GRADIENT_ONLY);
 236         g_new(0)=g_av(axis);
 237       }
 238       GC++;
 239
 240       // y=g_new-g
 241       copy(g_new,y);
 242       axpy(-1.0,g,y);
 243
 244       // Check curvature condition
 245       alpha=dot(y,h_dl);
 246       if (alpha <= sqrt(MacEpsilon)*nrm_dl*norm2(y)) {
 247 #ifdef LS_DEBUG
 248         cout << "Curvature condition violated " ;
 249 #endif
 250       }
 251       else {
 252         // Update Hessian
 253         gemv('N',1.0,B,h_dl,0.0,z) ; // z=Bh_dl
 254         beta=-1/dot(h_dl,z) ;
 255         ger(1/alpha,y,y,B) ;
 256         ger(beta,z,z,B) ;
 257
 258         // Update Hessian inverse
 259         gemv('N',1.0,H,y,0.0,z) ; // z=H*y
 260         gemv('T',1.0,H,y,0.0,w) ; // w=y'*H
 261         beta=dot(y,z) ;
 262         beta=(1+beta/alpha)/alpha ;
 263
 264         // It should be possible to do this updating more efficiently, by
 265         // exploiting the fact that (h_dl*y'*H) = transpose(H*y*h_dl')
 266         ger(beta,h_dl,h_dl,H) ;
 267         ger(-1/alpha,z,h_dl,H) ;
 268         ger(-1/alpha,h_dl,w,H) ;
 269       }
 270
 271       if (nrm_dl < norm2(x)*epsilon) {
 272         // Stop criterion (iteration progress) fullfilled
 273 #ifdef LS_DEBUG
 274         cout << "Progress is marginal" ;
 275 #endif
 276         good_enough = 1 ;
 277       }
 278
 279       // Check if we are close to a stationary point located previously
 280       if (box.CloseToMin(x_new, &f_new, eps_cl)) {
 281         // Note that x_new and f_new may be overwritten on exit from CloseToMin
 282 #ifdef LS_DEBUG
 283         cout << "Close to a previously located stationary point, exiting" << endl;
 284 #endif
 285         info = LS_Old ;
 286         good_enough = 1 ;
 287       }
 288
 289       // Update x, g and f
 290       copy(x_new,x) ; copy(g_new,g) ; f=f_new ;
 291
 292 #ifdef LS_DEBUG
 293       cout << " x=" << x << endl ;
 294 #endif
 295
 296     }
 297     else {
 298 #ifdef LS_DEBUG
 299       cout << "Step is no good, ro=" << ro << " delta=" << delta << endl ;
 300 #endif
 301     }
 302
 303   } // wend
 304
 305   // Make sure the routine returns correctly...
 306   // Check if last iterate is outside the boundary
 307   if (box.OutsideBox(x, domain) != 0) {
 308     info=LS_Out; f=DBL_MAX;
 309   }
 310
 311   if (info == LS_Unstable) {
 312     cout << "Local search became unstable. No big deal but exiting anyway\n" ;
 313     exit(1);
 314   }
 315
 316   T.xvals=x ; T.objval=f ;
 317   *mgr=maxgrad ;
 318   if (outside>0)
 319     return LS_Out ;
 320   else
 321     return info ;
 322 }