lm.cpp

用gsl实现Levenberg-Marquardt算法

#define GSL_DLL
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlin.h>

#include "expfit.h"

#ifdef _DEBUG
#pragma comment(lib, "libgsl_dll_d.lib")
#pragma comment(lib,"libgslcblas_dll_d.lib")
#else
#pragma comment(lib, "libgsl_dll.lib")
#pragma comment(lib,"libgslcblas_dll.lib")
#endif


#define N 40

void print_state (size_t iter, gsl_multifit_fdfsolver * s);

int
main (void)
{
	const gsl_multifit_fdfsolver_type *T;
	gsl_multifit_fdfsolver *s;
	int status;
	unsigned int i, iter = 0;
	const size_t n = N;
	const size_t p = 3;
	
	gsl_matrix *covar = gsl_matrix_alloc (p, p);
	double y[N], sigma[N];
	struct data d = { n, y, sigma};
	gsl_multifit_function_fdf f;
	double x_init[3] = { 1.0, 0.0, 0.0 };//参数初始值
	gsl_vector_view x = gsl_vector_view_array (x_init, p);//这一句有问题

	const gsl_rng_type * type;
	gsl_rng * r;
	
	gsl_rng_env_setup();
	
	type = gsl_rng_default;
	r = gsl_rng_alloc (type);
	
	f.f = &expb_f;
	f.df = &expb_df;
	f.fdf = &expb_fdf;
	f.n = n;
	f.p = p;
	f.params = &d;
	
	/* 生成一组待拟合数据 */
	for (i = 0; i < n; i++)
	{
		double t = i;
		y[i] = 1.0 + 5 * exp (-0.1 * t)//真实的参数为1.0, 5, 0.1 
			+ gsl_ran_gaussian (r, 0.1);
		//y[i]是真实的值再加上一点误差
		sigma[i] = 0.1;
		printf ("data: %u %g %g\n", i, y[i], sigma[i]);

	};
	
	//设置优化器相关环境
	T = gsl_multifit_fdfsolver_lmsder;//这里lmsder中的lm表示levenberg-marquardt
	s = gsl_multifit_fdfsolver_alloc (T, n, p);
	gsl_multifit_fdfsolver_set (s, &f, &x.vector);
	
	print_state (iter, s);//打印初始的迭代信息

	do
	{
		iter++;
		status = gsl_multifit_fdfsolver_iterate (s);
		
		printf ("status = %s\n", gsl_strerror (status));
		
		print_state (iter, s);
		
		if (status)
			break;
		
		//判断是否该停止迭代
		status = gsl_multifit_test_delta (s->dx, s->x,
			1e-4, 1e-4);
	}
	while (status == GSL_CONTINUE && iter < 500);
	
	gsl_multifit_covar (s->J, 0.0, covar);
	
#define FIT(i) gsl_vector_get(s->x, i)
#define ERR(i) sqrt(gsl_matrix_get(covar,i,i))
	
	{ 
		double chi = gsl_blas_dnrm2(s->f);
		double dof = n - p;
		double c = GSL_MAX_DBL(1, chi / sqrt(dof)); 
		
		printf("chisq/dof = %g\n",  pow(chi, 2.0) / dof);
		
		printf ("A      = %.5f +/- %.5f\n", FIT(0), c*ERR(0));
		printf ("lambda = %.5f +/- %.5f\n", FIT(1), c*ERR(1));
		printf ("b      = %.5f +/- %.5f\n", FIT(2), c*ERR(2));
	}
	
	printf ("status = %s\n", gsl_strerror (status));
	

	gsl_multifit_fdfsolver_free (s);
	gsl_matrix_free (covar);
	gsl_rng_free (r);
	return 0;
}

void
print_state (size_t iter, gsl_multifit_fdfsolver * s)
{
	printf ("iter: %3u x = % 15.8f % 15.8f % 15.8f "
		"|f(x)| = %g\n",
		iter,
		gsl_vector_get (s->x, 0), 
		gsl_vector_get (s->x, 1),
		gsl_vector_get (s->x, 2), 
		gsl_blas_dnrm2 (s->f));
}