Sample code for Procedure 4.1

The procedure in the following page is the C implementation of the local optimization procedure given in Chapter 4. It has a probing vector that adapts its size and direction according to the local properties of the objective function. The constant NDIM is the number of dimensions, THRESHOLD is the minimum step size, INIT_SIZE is the initial step size. The variable count keeps the number of the calls to the objective function. The two arguments to the procedure are the objective function f and the initial starting point x. The procedure prints out successively better solutions and leaves the final solution in the argument vector x.

void local_optimize(double (*f)(), double x[])
{
  double u[NDIM], v[NDIM], xv[NDIM];
  double fx, fxv, vr;
  int vi, vvec, i, iter, maxiter;

  for(i=0; i<NDIM; i++) u[i] = v[i] = 0;
  vi = -1; vvec = 1; vr = -INIT_SIZE;
  fx = f(x); fxv = 1E10;
  printf("%d. %.4f <= ", count, fx);
  for(i=0; i<NDIM; i++) { printf("%.4f ", x[i]); }; printf("\n");
  while(fabs(vr) >= THRESHOLD) {
    maxiter = ((fabs(vr) < 2*THRESHOLD) ? 2*NDIM : 2);
    iter = 0;
    while((fxv >= fx) && (iter < maxiter)) {
      if(iter == 0) { for(i=0; i<NDIM; i++) xv[i] = x[i]; }
      else xv[vi] -= vr;
      if(vvec) vvec = 0;
      vr = -vr;
      if(vr > 0) vi = ((vi+1) % NDIM);
      xv[vi] += vr;
      fxv = f(xv);
      iter++;
    }
    if(fxv >= fx) {
      fxv = 1E10;
      vr /= 2;
    } else {
      fx = fxv; printf("%d. %.4f <= ", count, fx);
      for(i=0; i<NDIM; i++) { x[i] = xv[i]; printf("%.4f ", x[i]); }
      printf("\n");
      if(iter == 0) {
	if(vvec) {
	  for(i=0; i<NDIM; i++) { u[i] += v[i]; v[i] *= 2; xv[i] += v[i]; }
	  vr *= 2;
	} else {
	  u[vi] += vr; vr *= 2; xv[vi] += vr;
	}
	fxv = f(xv);
      } else {
	for(i=0; i<NDIM; i++) xv[i] += u[i];
	xv[vi] += vr;
	fxv = f(xv);
	if(fxv >= fx) {
	  for(i=0; i<NDIM; i++) { u[i] = 0; xv[i] = x[i]; }
	  u[vi] = vr; vr *= 2;
	  xv[vi] += vr; fxv = f(xv);
	} else {
	  for(i=0; i<NDIM; i++) x[i] = xv[i]; fx = fxv;
	  u[vi] += vr;
	  for(i=0; i<NDIM; i++) v[i] = 2*u[i]; vvec = 1;
	  for(i=0; i<NDIM; i++) xv[i] += v[i]; fxv = f(xv);
	  for(vr=0,i=0; i<NDIM; i++) vr += v[i]*v[i];
	  vr = sqrt(vr);
	}}}}}