Hello all,
I've been finding that the new simplex algorithm seems to fail
to converge much more frequently than the old one. It took me
awhile to generate a sufficiently portable case which demonstrates
this, but attached is a program which shows that the simplex
minimizer succeeds, but the simplex2 minimizer stalls and gets
stuck (the true minimum is at (1,0,0)).
Unfortunately I haven't yet figured out why it's stalling, but
in the meantime I thought I'd post this code and see if anyone
else has been having similar difficulties. It's a bit of an unusual
function, but it's continuous everywhere and all of the other
minimizers I try happily succeed without complaint.
Take care,
Andrew
( http://o2scl.sourceforge.net )
#include <stdlib.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_errno.h>
double pi;
double spring(const gsl_vector *x, void *pa) {
double x0=gsl_vector_get(x,0);
double x1=gsl_vector_get(x,1);
double x2=gsl_vector_get(x,2);
double theta=atan2(x1,x0);
double r=sqrt(x0*x0+x1*x1);
double z=x2;
while (z>pi) z-=2.0*pi;
while (z<-pi) z+=2.0*pi;
double tmz=theta-z;
double rm1=r-1.0;
double ret=exp(tmz*tmz+rm1*rm1)+fabs(x2/10.0);
if (!finite(ret)) exit(-1);
return ret;
}
int main(int argv, char *argc[]) {
pi=acos(-1.0);
int nvar=3;
{
const gsl_multimin_fminimizer_type *T=gsl_multimin_fminimizer_nmsimplex;
gsl_multimin_fminimizer *s=0;
gsl_multimin_function minex_func;
gsl_vector *gx=gsl_vector_alloc(nvar);
gsl_vector_set(gx,0,1.0);
gsl_vector_set(gx,1,0.0);
gsl_vector_set(gx,2,7.0*pi);
int iter=0;
int status;
double size;
gsl_vector *ss;
ss=gsl_vector_alloc(nvar);
gsl_vector_set_all(ss,1.0);
minex_func.n=nvar;
minex_func.f=&spring;
minex_func.params=0;
s=gsl_multimin_fminimizer_alloc(T,nvar);
gsl_multimin_fminimizer_set(s,&minex_func,gx,ss);
do {
iter++;
status=gsl_multimin_fminimizer_iterate(s);
if (status) break;
size=gsl_multimin_fminimizer_size(s);
status=gsl_multimin_test_size(size,1.0e-4);
if (iter%20==0) {
printf("%d %6.5e %6.5e %6.5e\n",iter,
gsl_vector_get(s->x,0),
gsl_vector_get(s->x,1),
gsl_vector_get(s->x,2));
}
} while (status == GSL_CONTINUE && iter < 600);
gsl_vector_free(ss);
gsl_vector_free(gx);
gsl_multimin_fminimizer_free(s);
}
printf("\n");
{
const gsl_multimin_fminimizer_type *T=gsl_multimin_fminimizer_nmsimplex2;
gsl_multimin_fminimizer *s=0;
gsl_multimin_function minex_func;
gsl_vector *gx=gsl_vector_alloc(nvar);
gsl_vector_set(gx,0,1.0);
gsl_vector_set(gx,1,0.0);
gsl_vector_set(gx,2,7.0*pi);
int iter=0;
int status;
double size;
gsl_vector *ss;
ss=gsl_vector_alloc(nvar);
gsl_vector_set_all(ss,1.0);
minex_func.n=nvar;
minex_func.f=&spring;
minex_func.params=0;
s=gsl_multimin_fminimizer_alloc(T,nvar);
gsl_multimin_fminimizer_set(s,&minex_func,gx,ss);
do {
iter++;
status=gsl_multimin_fminimizer_iterate(s);
if (status) break;
size=gsl_multimin_fminimizer_size(s);
status=gsl_multimin_test_size(size,1.0e-4);
if (iter%20==0) {
printf("%d %6.5e %6.5e %6.5e\n",iter,
gsl_vector_get(s->x,0),
gsl_vector_get(s->x,1),
gsl_vector_get(s->x,2));
}
} while (status == GSL_CONTINUE && iter < 600);
gsl_vector_free(ss);
gsl_vector_free(gx);
gsl_multimin_fminimizer_free(s);
}
return 0;
}
_______________________________________________
Bug-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/bug-gsl
