Dear Sir/Madam,
I am currently working on a gsl least squares method. In this specific case , i have it set to the levenberg-marquardt with geodesic acceleration. i used the same example from the gsl website and added a few changes to it like so : #include <stdlib.h>#include <stdio.h>#include <gsl/gsl_vector.h>#include <gsl/gsl_matrix.h>#include <gsl/gsl_blas.h>#include <gsl/gsl_multifit_nlinear.h>#include <gsl/gsl_rng.h>#include <gsl/gsl_errno.h>#include <gsl/gsl_randist.h> #define max_iterations 100 struct data { double* t; double* y; size_t n; }; static double scaled_norm(const gsl_vector* D, const gsl_vector* a) { const size_t n = a->size; double e2 = 0.0; size_t i; for (i = 0; i < n; ++i) { double Di = gsl_vector_get(D, i); double ai = gsl_vector_get(a, i); double u = Di * ai; e2 += u * u; } return sqrt(e2); } /* model function: a * exp( -1/2 * [ (t - b) / c ]^2 ) */doublegaussian(const double a, const double b, const double c, const double t) { const double z = (t - b) / c; return (a * exp(-0.5 * z * z)); } intfunc_f(const gsl_vector* x, void* params, gsl_vector* f) { struct data* d = (struct data*)params; double a = gsl_vector_get(x, 0); double b = gsl_vector_get(x, 1); double c = gsl_vector_get(x, 2); size_t i; for (i = 0; i < d->n; ++i) { double ti = d->t[i]; double yi = d->y[i]; double y = gaussian(a, b, c, ti); gsl_vector_set(f, i, yi - y); } return GSL_SUCCESS; } intfunc_df(const gsl_vector* x, void* params, gsl_matrix* J) { struct data* d = (struct data*)params; double a = gsl_vector_get(x, 0); double b = gsl_vector_get(x, 1); double c = gsl_vector_get(x, 2); size_t i; for (i = 0; i < d->n; ++i) { double ti = d->t[i]; double zi = (ti - b) / c; double ei = exp(-0.5 * zi * zi); gsl_matrix_set(J, i, 0, -ei); gsl_matrix_set(J, i, 1, -(a / c) * ei * zi); gsl_matrix_set(J, i, 2, -(a / c) * ei * zi * zi); } return GSL_SUCCESS; } intfunc_fvv(const gsl_vector* x, const gsl_vector* v, void* params, gsl_vector* fvv) { struct data* d = (struct data*)params; double a = gsl_vector_get(x, 0); double b = gsl_vector_get(x, 1); double c = gsl_vector_get(x, 2); double va = gsl_vector_get(v, 0); double vb = gsl_vector_get(v, 1); double vc = gsl_vector_get(v, 2); size_t i; for (i = 0; i < d->n; ++i) { double ti = d->t[i]; double zi = (ti - b) / c; double ei = exp(-0.5 * zi * zi); double Dab = -zi * ei / c; double Dac = -zi * zi * ei / c; double Dbb = a * ei / (c * c) * (1.0 - zi * zi); double Dbc = a * zi * ei / (c * c) * (2.0 - zi * zi); double Dcc = a * zi * zi * ei / (c * c) * (3.0 - zi * zi); double sum; sum = 2.0 * va * vb * Dab + 2.0 * va * vc * Dac + vb * vb * Dbb + 2.0 * vb * vc * Dbc + vc * vc * Dcc; gsl_vector_set(fvv, i, sum); } return GSL_SUCCESS; } void callback(const size_t iter, void* params, const gsl_multifit_nlinear_workspace* w) { gsl_vector* f = gsl_multifit_nlinear_residual(w); gsl_vector* x = gsl_multifit_nlinear_position(w); double avratio = gsl_multifit_nlinear_avratio(w); double rcond; (void)params; /* not used */ gsl_multifit_nlinear_trust_state* state = (gsl_multifit_nlinear_trust_state*)w->state; /* compute reciprocal condition number of J(x) */ gsl_multifit_nlinear_rcond(&rcond, w); fprintf(stderr, "iter %2zu: a = %.4f, b = %.4f, c = %.4f, |a|/|v| = %.4f cond(J) = %8.4f, |f(x)| = %.4f\n", iter, gsl_vector_get(x, 0), gsl_vector_get(x, 1), gsl_vector_get(x, 2), avratio, 1.0 / rcond, gsl_blas_dnrm2(f)); } voidsolve_system(gsl_vector* x, gsl_multifit_nlinear_fdf* fdf, gsl_multifit_nlinear_parameters* params) { const gsl_multifit_nlinear_type* T = gsl_multifit_nlinear_trust; const size_t max_iter = 200; const double xtol = 1.0e-8; const double gtol = 1.0e-8; const double ftol = 1.0e-8; const size_t n = fdf->n; const size_t p = fdf->p; size_t iter = 0; int status, info; gsl_multifit_nlinear_workspace* work = gsl_multifit_nlinear_alloc(T, params, n, p); gsl_vector* f = gsl_multifit_nlinear_residual(work); gsl_vector* y = gsl_multifit_nlinear_position(work); double chisq0, chisq, rcond; double trust_radius_list[max_iterations]; double step_size_list[max_iterations]; /* initialize solver */ gsl_multifit_nlinear_init(x, fdf, work); /* store initial cost */ gsl_blas_ddot(f, f, &chisq0); /* iterate until convergence */ if (callback) callback(iter, NULL, work); /*initialize scale and trust region*/ gsl_vector* diag = gsl_vector_alloc(p); work->params.scale->init(work->J, diag); double scale = scaled_norm(diag, work->x); double trust_region_radius = 0.3 * GSL_MAX(1.0, scale); trust_radius_list[0] = trust_region_radius; // initialize the trust region radius step_size_list[0] = gsl_blas_dnrm2(work->dx); do { status = gsl_multifit_nlinear_iterate(work); /*update the trust region information*/ work->params.scale->update(work->J, diag); scale = scaled_norm(diag, work->dx); trust_region_radius = scale; trust_radius_list[iter + 1] = trust_region_radius; step_size_list[iter + 1] = gsl_blas_dnrm2(work->dx); if (status == GSL_ENOPROG && iter == 0) { info = status; } ++iter; if (callback) callback(iter,NULL, work); /* test for convergence */ status = gsl_multifit_nlinear_test(xtol, gtol, ftol, &info, work); } while (status == GSL_CONTINUE && iter < max_iter); if (status == GSL_ETOLF || status == GSL_ETOLX || status == GSL_ETOLG) { info = status; status = GSL_SUCCESS; } /* check if max iterations reached */ if (iter >= max_iter && status != GSL_SUCCESS) status = GSL_EMAXITER; int i; for (i = 0; i < iter + 1; i++) { fprintf(stdout, "trust-radius: %g ", trust_radius_list[i]); } /* store final cost */ gsl_blas_ddot(f, f, &chisq); /* store cond(J(x)) */ gsl_multifit_nlinear_rcond(&rcond, work); gsl_vector_memcpy(x, y); /* print summary */ fprintf(stderr, "NITER = %zu\n", gsl_multifit_nlinear_niter(work)); fprintf(stderr, "NFEV = %zu\n", fdf->nevalf); fprintf(stderr, "NJEV = %zu\n", fdf->nevaldf); fprintf(stderr, "NAEV = %zu\n", fdf->nevalfvv); fprintf(stderr, "initial cost = %.12e\n", chisq0); fprintf(stderr, "final cost = %.12e\n", chisq); fprintf(stderr, "final x = (%.12e, %.12e, %12e)\n", gsl_vector_get(x, 0), gsl_vector_get(x, 1), gsl_vector_get(x, 2)); fprintf(stderr, "final cond(J) = %.12e\n", 1.0 / rcond); gsl_multifit_nlinear_free(work); } intmain(void) { const size_t n = 300; /* number of data points to fit */ const size_t p = 3; /* number of model parameters */ const double a = 5.0; /* amplitude */ const double b = 0.4; /* center */ const double c = 0.15; /* width */ const gsl_rng_type* T = gsl_rng_default; gsl_vector* f = gsl_vector_alloc(n); gsl_vector* x = gsl_vector_alloc(p); gsl_multifit_nlinear_fdf fdf; gsl_multifit_nlinear_parameters fdf_params = gsl_multifit_nlinear_default_parameters(); struct data fit_data; gsl_rng* r; size_t i; gsl_rng_env_setup(); r = gsl_rng_alloc(T); fit_data.t = (double *)malloc(n * sizeof(double)); fit_data.y = (double *)malloc(n * sizeof(double)); fit_data.n = n; /* generate synthetic data with noise */ for (i = 0; i < n; ++i) { double t = (double)i / (double)n; double y0 = gaussian(a, b, c, t); double dy = gsl_ran_gaussian(r, 0.1 * y0); fit_data.t[i] = t; fit_data.y[i] = y0 + dy; } /* define function to be minimized */ fdf.f = func_f; fdf.df = func_df; fdf.fvv = func_fvv; fdf.n = n; fdf.p = p; fdf.params = &fit_data; /* starting point */ gsl_vector_set(x, 0, 1.0); gsl_vector_set(x, 1, 0.0); gsl_vector_set(x, 2, 1.0); fdf_params.trs = gsl_multifit_nlinear_trs_lmaccel; solve_system(x, &fdf, &fdf_params); /* print data and model */ { double A = gsl_vector_get(x, 0); double B = gsl_vector_get(x, 1); double C = gsl_vector_get(x, 2); for (i = 0; i < n; ++i) { double ti = fit_data.t[i]; double yi = fit_data.y[i]; double fi = gaussian(A, B, C, ti); printf("%f %f %f\n", ti, yi, fi); } } gsl_vector_free(f); gsl_vector_free(x); gsl_rng_free(r); return 0; } This code initializes the trust radius and updates it and stores it in a list. I started with the same initialization method as trust.c which a source file in the gsl library used to calculate the trust region information. work->params.scale->init(work->J, diag); double scale = scaled_norm(diag, work->x); double trust_region_radius = 0.3 * GSL_MAX(1.0, scale); but for the update , i used the step instead: work->params.scale->update(work->J, diag); scale = scaled_norm(diag, work->dx); trust_region_radius = scale; i did not do factor up and down to make the trust region bigger or smaller since i assumed that that was already done if i multiply the diagonal of the scaling matrix with the step. I tested the values by debugging both from the source file and from my method and they were different. the factor up and factor down parameters are accessed like so : work->params.factor_up; work->params.factor_down; and those increase or decrease the trust region. is the reason why my code isn't giving similar results because i didnt use factor_up and factor_down. even though they should already be taken into account by the step? Keep in mind that my trust region calculation is done the other way around , meaning that the trust region was already calculated internally but i am trying to recalculate it so i can plot it. i also realized that the predicted reduction is calculated differently for every method. I would really appreciate your help with this. thank you so much, Sincerely, Aleja Carmento.