🚀 go-pugleaf

In article <MATT.93Jan10152553@physics2.berkeley.edu> matt@physics2.berkeley.edu (Matt Austern) writes:
>From: matt@physics2.berkeley.edu (Matt Austern)
>Subject: Wanted: minimization routine
>Date: 10 Jan 93 15:25:53
>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.
Try the Nelder's simplex method.
(Computer Journal, vol. 7, pp 308-313, 1964)
Eduardo Lucero

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 a few variables (N<12), I have had good luck with
Hooke and Jeeves.  It will minimize discontinuous functions,
and/or functions whose derivatives are discontinuous.

Hooke and Jeeves does not use derivatives.

It is a heuristic method, from 1963, and has fallen VERY MUCH
into disfavour among academics because no proof exists of
its convergence properties.  (Whereas other methods, notably
those of quasi-Newton parentage like Huangs class of Variable
Metric algos, have provably quadratic convergence given
certain assumptions -- assumptions that include a damn good
starting point for the optimization search).

Pragmatically speaking, to quote Fred Dryer, "It works for me."
And because it doesn't need or use derivatives, it is usually
(again, pragmatically speaking) quicker to get a problem up and
running on H&J.

If it works for you too, great.  If not, you are entitled to
a full refund of your purchase price.


#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]);

}