This is the mail archive of the gsl-discuss@sourceware.org mailing list for the GSL project.
Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|
Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |
Other format: | [Raw text] |
hi all, currently the nm_simplex minimizer is initialized by giving it a starting point (vector) and a vector of initial stepsizes. from this nm_simplex_set constructs the initial simplex. i had sometimes difficulties with this initializer because it results in a "very regular" initial simplex. the first trial steps often only change one parameter. (just as an example. i tried to optimize the contour of a mirror defined by a spline interpolation of several points. the current minimizer would start by trying to shift only individual points. this, of course were not successful steps. it then contracted for a long time until it finally took off correctly sometimes). so my proposal (see patch) is to change the call to gsl_multimin_fminimizer_set() to use a starting vector and a single value for the initial step size. nm_simplex_set() then contructs a regular n-pod of size initial_step_size with the starting point at its center. you can fine-tune the orientation of this initial simplex by using different values for the environment variable GSL_NM_SIMPLEX_SEED. notes: -make test in multimin fails to build because it needs to link with ../rng/libgslrng.la (you have to edit the makefile to make it work. i just did not find where to configure this thing) -because this is a change in the api it may make sense to change it even more drastic i.e. to only use one generic gsl_multimin_minimizer_set() (the one currently used for the fdf minimizer) and ignore the parameters that do not apply for the given minimizer. this would be similar to how some ode solvers make use of the jacobian and some do not. cheers -- Dr. Ivo Alxneit Laboratory for Solar Technology phone: +41 56 310 4092 Paul Scherrer Institute fax: +41 56 310 2688 CH-5232 Villigen http://solar.web.psi.ch Switzerland gnupg key: 0x515E30C7
diff -Naur /usr/local/src/gsl-1.7/multimin/fminimizer.c multimin/fminimizer.c --- /usr/local/src/gsl-1.7/multimin/fminimizer.c 2005-06-26 15:25:35.000000000 +0200 +++ multimin/fminimizer.c 2006-02-23 16:16:42.000000000 +0100 @@ -74,14 +74,14 @@ gsl_multimin_fminimizer_set (gsl_multimin_fminimizer * s, gsl_multimin_function * f, const gsl_vector * x, - const gsl_vector * step_size) + double step_size) { if (s->x->size != f->n) { GSL_ERROR ("function incompatible with solver size", GSL_EBADLEN); } - if (x->size != f->n || step_size->size != f->n) + if (x->size != f->n) { GSL_ERROR ("vector length not compatible with function", GSL_EBADLEN); } @@ -89,8 +89,10 @@ s->f = f; gsl_vector_memcpy (s->x,x); - - return (s->type->set) (s->state, s->f, s->x, &(s->size), step_size); + + s->size = step_size; + + return (s->type->set) (s->state, s->f, s->x, step_size); } void diff -Naur /usr/local/src/gsl-1.7/multimin/gsl_multimin.h multimin/gsl_multimin.h --- /usr/local/src/gsl-1.7/multimin/gsl_multimin.h 2005-06-26 15:25:35.000000000 +0200 +++ multimin/gsl_multimin.h 2006-02-23 16:16:33.000000000 +0100 @@ -84,8 +84,7 @@ int (*alloc) (void *state, size_t n); int (*set) (void *state, gsl_multimin_function * f, const gsl_vector * x, - double * size, - const gsl_vector * step_size); + double step_size); int (*iterate) (void *state, gsl_multimin_function * f, gsl_vector * x, double * size, @@ -117,7 +116,7 @@ gsl_multimin_fminimizer_set (gsl_multimin_fminimizer * s, gsl_multimin_function * f, const gsl_vector * x, - const gsl_vector * step_size); + double step_size); void gsl_multimin_fminimizer_free(gsl_multimin_fminimizer *s); diff -Naur /usr/local/src/gsl-1.7/multimin/simplex.c multimin/simplex.c --- /usr/local/src/gsl-1.7/multimin/simplex.c 2005-06-26 15:25:35.000000000 +0200 +++ multimin/simplex.c 2006-02-23 16:16:21.000000000 +0100 @@ -33,7 +33,9 @@ #include <config.h> #include <stdlib.h> +#include <math.h> #include <gsl/gsl_blas.h> +#include <gsl/gsl_rng.h> #include <gsl/gsl_multimin.h> typedef struct @@ -45,6 +47,169 @@ } nmsimplex_state_t; +static double a_m(const int m, const double c) +{ + return (c * sqrt(1 - c) / sqrt((1 + (m - 1) * c) * (1 + m * c))); +} + +static double b_m(const size_t m, const double c) +{ + return (sqrt((1 - c) * (1 + m * c) / (1 + (m - 1) * c))); +} + +static gsl_matrix *init_npod(const size_t n) +{ +/* + * given the following n+1 vectors ei: 0<i<=n in n dimensions + * + * e0(a0, 0, 0, ...) + * e1(a0, b1, 0, ..) + * e2(a0, a1, b2, .) + * en-1(a0, ......, bn) + * en(a0, ....., -bn) + * + * that form a regular normalized n+1_pod in n dimensions, i.e. + * a regular triangle in 2D or a regulart tetrahedron in 3D + * + * requesting that + * all (ei, ei) = 1 + * all (ei, ej) = c + */ + + size_t i, j; + const double c = -1.0 / n; + gsl_vector_view v; + gsl_matrix *e; + + + e = gsl_matrix_calloc(n + 1, n); + if (e == NULL) + { + return (NULL); + } + +/* + * init e_0 + */ + gsl_matrix_set(e, 0, 0, 1.0); + +/* + * init e_1 to e_n-1 + */ + for (i = 1; i < n; i++) + { + for (j = 0; j < i; j++) + gsl_matrix_set(e, i, j, a_m((int) j, c)); + + gsl_matrix_set(e, i, j, b_m(j, c)); + + for (j = i + 1; j < n; j++) + gsl_matrix_set(e, i, j, 0.0); + } + +/* + * init e_n + */ + v = gsl_matrix_row(e, n - 1); + gsl_matrix_set_row(e, n, &v.vector); + + gsl_matrix_set(e, n, n - 1, -gsl_matrix_get(e, n - 1, n - 1)); + + return (e); + +} + +static gsl_matrix *init_basis(const size_t n) +{ +/* + * orthonolmal basis, randomly oriented + */ + size_t i, j; + char *val; + unsigned int seed; + gsl_rng *r; + gsl_matrix *a, *basis; + + a = gsl_matrix_alloc(n, n); + if (a == NULL) + { + return (NULL); + } + + basis = gsl_matrix_alloc(n, n); + if (basis == NULL) + { + gsl_matrix_free(a); + return (NULL); + } + + val = getenv("GSL_NM_SIMPLEX_SEED"); + if (val) + seed = strtoul(val, 0, 0); + else + seed = 0; + + r = gsl_rng_alloc(gsl_rng_taus); + if (r == NULL) + { + gsl_matrix_free(a); + gsl_matrix_free(basis); + return (NULL); + } + gsl_rng_set(r, seed); + + +/* + * get n random vectors of length n + */ + for (i = 0; i < n; i++) + for (j = 0; j < n; j++) + gsl_matrix_set(a, i, j, 2.0 * gsl_rng_uniform(r) - 1.0); + + gsl_rng_free(r); + +/* + * orthogonalize + * + * FIXME: we do not check for colinear vectors + */ + gsl_matrix_memcpy(basis, a); + + for (i = 0; i < n; i++) + { + gsl_vector_view v1 = gsl_matrix_row(a, i); + gsl_vector_view v2 = gsl_matrix_row(basis, i); + + for (j = 0; j < i; j++) + { + double t1, t2; + gsl_vector_view v3 = gsl_matrix_row(basis, j); + + gsl_blas_ddot(&v1.vector, &v3.vector, &t1); + gsl_blas_ddot(&v3.vector, &v3.vector, &t2); + + gsl_blas_daxpy(-t1 / t2, &v3.vector, &v2.vector); + + } + } + + gsl_matrix_free(a); + +/* + * normalize basis vectors + */ + for (i = 0; i < n; i++) + { + gsl_vector_view v1 = gsl_matrix_row(basis, i); + double t1 = gsl_blas_dnrm2(&v1.vector); + + gsl_vector_scale(&v1.vector, 1 / t1); + + } + + return (basis); + +} static double nmsimplex_move_corner (const double coeff, const nmsimplex_state_t * state, size_t corner, gsl_vector * xc, @@ -216,47 +381,65 @@ } static int -nmsimplex_set (void *vstate, gsl_multimin_function * f, - const gsl_vector * x, - double *size, const gsl_vector * step_size) +nmsimplex_set(void *vstate, gsl_multimin_function * f, + const gsl_vector * x, double step_size) { int status; size_t i; - double val; + gsl_matrix *n_pod, *basis; nmsimplex_state_t *state = (nmsimplex_state_t *) vstate; - gsl_vector *xtemp = state->ws1; + size_t n = x->size; + +/* + * prepare regular, normalized n_pod + * equilateral triangle 2D + * equilateral tetrahedron 3D + */ + n_pod = init_npod(n); + if (n_pod == NULL) + { + GSL_ERROR("failed to allocate space for n_pod", GSL_ENOMEM); + } - /* first point is the original x0 */ +/* + * construct a randomly oriented, orthonormal basis in n dimensions + */ + basis = init_basis(n); + if (basis == NULL) + { + gsl_matrix_free(n_pod); + GSL_ERROR("failed to allocate space for basis", GSL_ENOMEM); + } - val = GSL_MULTIMIN_FN_EVAL (f, x); - gsl_matrix_set_row (state->x1, 0, x); - gsl_vector_set (state->y1, 0, val); +/* + * scale n_pod by step_size and express it in basis. result in state->x1 + */ + gsl_blas_dgemm(CblasNoTrans, CblasTrans, step_size, n_pod, basis, 0.0, + state->x1); - /* following points are initialized to x0 + step_size */ + gsl_matrix_free(basis); + gsl_matrix_free(n_pod); - for (i = 0; i < x->size; i++) +/* + * now state->x1 contains n+1 vectors forming a n+1 dimensional, randomly + * oriented simplex of size *size centered at the origin. + * we now offset if by the initial vector x and initialize the + * corresponding y values in state->y1. + */ + for (i = 0; i <= n; i++) { - status = gsl_vector_memcpy (xtemp, x); + double val; + gsl_vector_view v = gsl_matrix_row(state->x1, i); - if (status != 0) - { - GSL_ERROR ("vector memcopy failed", GSL_EFAILED); - } + gsl_vector_add(&v.vector, x); - val = gsl_vector_get (xtemp, i) + gsl_vector_get (step_size, i); - gsl_vector_set (xtemp, i, val); - val = GSL_MULTIMIN_FN_EVAL (f, xtemp); - gsl_matrix_set_row (state->x1, i + 1, xtemp); - gsl_vector_set (state->y1, i + 1, val); + val = GSL_MULTIMIN_FN_EVAL(f, &v.vector); + gsl_vector_set(state->y1, i, val); } - /* Initialize simplex size */ - - *size = nmsimplex_size (state); - - return GSL_SUCCESS; + return GSL_SUCCESS; } static void diff -Naur /usr/local/src/gsl-1.7/multimin/test.c multimin/test.c --- /usr/local/src/gsl-1.7/multimin/test.c 2006-02-22 16:27:14.000000000 +0100 +++ multimin/test.c 2006-02-24 08:43:43.000000000 +0100 @@ -58,7 +58,7 @@ test_fdf("Rosenbrock", &rosenbrock, rosenbrock_initpt,*T); T++; } -printf("simplex\n"); + test_f("Roth", &roth_fmin, roth_initpt); test_f("Wood", &wood_fmin, wood_initpt); test_f("Rosenbrock", &rosenbrock_fmin, rosenbrock_initpt); @@ -143,15 +143,12 @@ gsl_vector *x = gsl_vector_alloc (f->n); - gsl_vector *step_size = gsl_vector_alloc (f->n); + double step_size = 1.0; gsl_multimin_fminimizer *s; (*initpt) (x); - for (i = 0; i < f->n; i++) - gsl_vector_set (step_size, i, 1); - s = gsl_multimin_fminimizer_alloc(T, f->n); gsl_multimin_fminimizer_set (s, f, x, step_size); @@ -185,7 +182,6 @@ gsl_multimin_fminimizer_free(s); gsl_vector_free(x); - gsl_vector_free(step_size); return status; }
Attachment:
signature.asc
Description: This is a digitally signed message part
Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|
Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |