🚀 go-pugleaf

In article <MATT.93Jan10152553@physics2.berkeley.edu>, matt@physics2.berkeley.edu (Matt Austern) writes:

   > Does anybody out there have a recommendation for a minimization
   > routine?  I'll be using it to minimize a function of two or three
   > variables; the function is reasonably well behaved, but I am not able
   > to calculate the derivatives very easily.
   > Code in C or C++ would be preferable.  

For small numbers of variables (N<12), I find that the Hooke and Jeeves
algorithm [1963] works quite well.  Its good points are

  *1* Doesn't require the user to supply derivatives

  *2* Will minimize discontinuous functions

  *3* (A consequence of *2*) Will minimize functions whose derivatives
      are discontinuous

  *4* Very short subroutine (2700 bytes), easily transliterated from
      one programming language to another

  *5* (A consequence of *2*) One of the very quickest ways of getting
      a problem up and running.


The algorithm has VERY MUCH fallen into academic disfavour.  Principally
this is because it is a _heuristic_, and because no proofs are known of
its rates of convergence.  The darlings of the I-gotta-see-a-proof
community are called "Huang's class of Variable Metric" algorithms,
and their parentage is quasi-Newton.  Generally these require that
the function have continuous derivatives.  Often they require that
the subroutine user _supply_ the derivatives, leading to disadvantage
*5*.  However, for certain assumptions (including the cheerfully
unrealistic one "for starting points sufficiently close to the
optimum"), Huang's class has provably quadratic convergence.

So, if microseconds are life-or-death, if you must have the very very very
fastest possible algorithm, prepare to use a quasi-Newton method.  In
this case I recommend BFGS or LM, the latter available from netlib
(send index from minpack).

On the other hand, if you can accept a runtime of 8 seconds instead of 5,
and if you want avoid derivatives, and if you are using small numbers
of variables, then it's worth giving Hooke and Jeeves a try.

Pragmatically speaking, quoting Fred Dryer, "It works for me."  I use it
for modeling semiconductor devices, especially short-channel MOSFETs.

   Technical note: because H&J accepts functions with discontinuous derivatives,
                   it mates beautifully with penalty functions to give quick
                   and easy _constrained_ minimization.  And you don't have
                   to sweat the details of slope-matching the penalties or
                   barriers.


#include <stdio.h>
#include <math.h>

#define		VARS		(100)    /* max # of variables */
#define		RBEGIN		(0.6)    /* stepsize geometric shrink   */
#define		DBEGIN		(0.6)    /* delta is a fraction of coord */
#define		EPSMIN		(1E-7)   /* ending value of stepsize */
#define		IMAX		(5000)   /* max # of iterations */


#define		ONETHIRD	(1.0/3.0)    /* used in one of the f() */


/* global variables */
/* global variables */

int	funevals = 0;





/* a few test functions -- to use one, change its name to f() */
/* a few test functions -- to use one, change its name to f() */
/* a few test functions -- to use one, change its name to f() */






/* Powells quartic function */
double	f(x, n)
	double	x[VARS] ;
	int	n ;
{
	double	a, b, c, d, t1, t2 ;

	funevals++;
	a = x[0] ;  b = x[1] ;  c = x[2] ;  d = x[3] ;

	t1  = (a + (10.0*b)) * (a + (10.0*b)) ;
	t1 += 5.0 * (c - d) * (c - d) ;
	t1 += (b - (2.0*c)) * (b - (2.0*c)) * (b - (2.0*c)) * (b - (2.0*c)) ;
	t1 += 10.0 * (a-d) * (a-d) * (a-d) * (a-d) ;

	return( t1 );
}






/* Rosenbrocks parabolic (banana) function */
double	rosen(x, n)
	double	x[VARS] ;
	int	n ;
{
	double	a, b, c ;

	funevals++;
	a = x[0] ;  b = x[1] ;
	c = 100.0 * (b - (a*a)) * (b - (a*a)) ;
	return(c + ( (1.0-a) * (1.0-a) ) ) ;
}







/* Schittowski problem 210: very difficult banana function */
double	prob210(x, n)
	double	x[VARS] ;
	int	n ;
{
	double	a, b, c ;

	funevals++;
	a = x[0] ;  b = x[1] ;
	c = 1.0e6 * (b - (a*a)) * (b - (a*a)) ;
	return(c + ( (1.0-a) * (1.0-a) ) ) ;
}







/* Schittowski problem 281: summation of ten cubics */
double	prob281(x, n)
	double	x[VARS] ;
	int	n ;
{
	double	summ, index ;
	int	i ;

	funevals++;

	summ = 0.0 ;
	for(i=0; i<n; i++)
	{
		index = (double) (1 + i) ;
		summ += (index * index * index) * (x[i] - 1.0) * (x[i] - 1.0) ;
	}

	summ = pow( summ, ONETHIRD ) ;
	return( summ );
}






/* given a point, look for a better one nearby, one coord at a time */

double	best_nearby(delta, point, prevbest, nvars)

	double	delta[VARS], point[VARS] ;
	double	prevbest ;
	int	nvars ;
{
	double	z[VARS] ;
	double	minf, ftmp ;
	int	i ;

	minf = prevbest ;
	for(i=0; i<nvars; i++) z[i] = point[i] ;

	for(i=0; i<nvars; i++)
	{
		z[i] = point[i] + delta[i] ;
		ftmp = f(z, nvars);
		if(ftmp < minf) minf = ftmp ;
		else
		{
			delta[i] = 0.0 - delta[i] ;
			z[i] = point[i] + delta[i] ;
			ftmp = f(z, nvars);
			if(ftmp < minf) minf = ftmp ;
			else z[i] = point[i] ;
		}
	}
	for(i=0; i<nvars; i++) point[i] = z[i] ;
	return(minf);
}







int	hooke(nvars, startpt, endpt, rho, epsilon, itermax)

	double	startpt[VARS], endpt[VARS] ;
	int	nvars, itermax ;
	double	rho, epsilon ;
{

	double	delta[VARS] ;
	double	newf, fbefore, steplength, tmp ;
	double	xbefore[VARS], newx[VARS] ;

	int	i, j, keep ;
	int	iters, iadj ;

	for(i=0; i<nvars; i++)
	{
		newx[i] = xbefore[i] = startpt[i] ;
		delta[i] = fabs(startpt[i] * DBEGIN) ;
		if(delta[i] == 0.0) delta[i] = DBEGIN ;
	}

	iadj = 0 ;
	steplength = rho ;
	iters = 0;
	fbefore = f(newx, nvars);
	newf = fbefore ;

	while( (iters < itermax) && (steplength > epsilon))
	{
		iters++;
		iadj++;


printf("\nAfter %5d funevals, f(x) =  %.4le at\n", funevals, fbefore);
for(j=0; j<nvars; j++) printf("   x[%2d] = %.4le\n", j, xbefore[j]) ;

		/* find best new point, one coord at a time */
		for(i=0; i<nvars; i++) { newx[i] = xbefore[i]; }
		newf = best_nearby(delta, newx, fbefore, nvars);


		/* if we made some improvements, pursue that direction */
		keep = 1 ;
		while( (newf < fbefore)  && (keep == 1) )
		{
			iadj = 0;

			for(i=0; i<nvars; i++)
			{
				/* firstly, arrange the sign of delta[] */
				if(newx[i] <= xbefore[i])
					delta[i] = 0.0 - fabs(delta[i]) ;
				else	delta[i] = fabs(delta[i]);

				/* now, move further in this direction */
				tmp = xbefore[i] ;
				xbefore[i] = newx[i] ;
				newx[i] = newx[i] + newx[i] - tmp ;
			}
			fbefore = newf ;
			newf = best_nearby(delta, newx, fbefore, nvars);


			/* if the further (optimistic) move was bad.... */
			if(newf >= fbefore) break ;


			/* make sure that the differences between the new */
			/* and the old points are due to actual */
			/* displacements; beware of roundoff errors that */
			/* might cause newf < fbefore */
			keep = 0;
			for(i=0; i<nvars; i++)
			{
				keep = 1 ;
				if(fabs(newx[i] - xbefore[i]) >
					(0.5 * fabs(delta[i])) ) break ;
				else keep = 0;
			}
		}

		if((steplength >= epsilon) && (newf >= fbefore) )
		{
			steplength = steplength * rho ;
			for(i=0; i<nvars; i++) { delta[i] *= rho ; }
		}

	}


	for(i=0; i<nvars; i++) endpt[i] = xbefore[i] ;
	return(iters);

}






main()
{

	double	startpt[VARS], endpt[VARS] ;
	int	nvars, itermax ;

	double	rho, epsilon ;
	int	i, jj ;



	/* test function problem 281 */
	/* starting guess */
	nvars = 10;

	startpt[0] =    1.7528323e-02 ;
	startpt[1] =   -2.2503207e+00 ;
	startpt[2] =   -1.0325044e+00 ;
	startpt[3] =    1.1047138e+00 ;
	startpt[4] =    6.7813668e-01 ;
	startpt[5] =    8.6220218e-01 ;
	startpt[6] =    9.6336051e-01 ;
	startpt[7] =    9.3734542e-01 ;
	startpt[8] =    9.9586705e-01 ;
	startpt[9] =    9.5568268e-01 ;


	/* rosenbrock test function   also prob 210 */
	/* starting guess */
	nvars = 2;
	startpt[0] = -1.2 ;
	startpt[1] =  1.0 ;



	/* powell test function */
	/* starting guess */
	nvars = 4;
	startpt[0] =  3.0 ;
	startpt[1] = -1.0 ;
	startpt[2] =  0.0 ;
	startpt[3] =  1.0 ;


	itermax = IMAX ;
	rho = RBEGIN;
	epsilon = EPSMIN ;

	jj = hooke(nvars, startpt, endpt, rho, epsilon, itermax) ;

	printf("\n\n\nHOOKE RETURNED\n");
	for(i=0; i<nvars; i++) printf("x[%3d] = %15.7le \n", i, endpt[i]);

}