GSL非线性最小二乘拟合调试报错及NaN/Inf结果问题求助

我基于GSL的gsl_multifit_nlinear接口实现非线性最小二乘拟合，参考官方示例拆分主函数以支持Levenberg-Marquardt和Dogleg两种求解方法。但运行程序返回包含NaN和Inf的异常结果，设置断点调试时也会报错，即使配置了trust region（trs）也无法解决。

问题代码 #

/* parameters to model */
#define N      100    /* number of data points to fit */
#define TMAX   (3.0)  /* time variable in [0,TMAX] */

struct data {
    size_t n;
    double* t;
    double* y;
};

int
expb_f(const gsl_vector* x, void* data,
    gsl_vector* f)
{
    size_t n = ((struct data*)data)->n;
    double* t = ((struct data*)data)->t;
    double* y = ((struct data*)data)->y;

    double A = gsl_vector_get(x, 0);
    double lambda = gsl_vector_get(x, 1);
    double b = gsl_vector_get(x, 2);

    size_t i;

    for (i = 0; i < n; i++)
    {
        /* Model Yi = A * exp(-lambda * t_i) + b */
        double Yi = A * exp(-lambda * t[i]) + b;
        gsl_vector_set(f, i, Yi - y[i]);
    }

    return GSL_SUCCESS;
}

int
expb_df(const gsl_vector* x, void* data,
    gsl_matrix* J)
{
    size_t n = ((struct data*)data)->n;
    double* t = ((struct data*)data)->t;

    double A = gsl_vector_get(x, 0);
    double lambda = gsl_vector_get(x, 1);

    size_t i;

    for (i = 0; i < n; i++)
    {
        /* Jacobian matrix J(i,j) = dfi / dxj, */
        /* where fi = (Yi - yi)/sigma[i],      */
        /*       Yi = A * exp(-lambda * t_i) + b  */
        /* and the xj are the parameters (A,lambda,b) */
        double e = exp(-lambda * t[i]);
        gsl_matrix_set(J, i, 0, e);
        gsl_matrix_set(J, i, 1, -t[i] * A * e);
        gsl_matrix_set(J, i, 2, 1.0);
    }

    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 rcond;

    /* compute reciprocal condition number of J(x) */
    gsl_multifit_nlinear_rcond(&rcond, w);

    fprintf(stderr, "iter %2zu: A = %.4f, lambda = %.4f, b = %.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),
        1.0 / rcond,
        gsl_blas_dnrm2(f));
}

void solve_system(gsl_vector* x, gsl_multifit_nlinear_fdf* fdf, gsl_multifit_nlinear_parameters* fdf_params)
{

    const gsl_multifit_nlinear_type* T = gsl_multifit_nlinear_trust;
    const size_t max_iter = 100;
    const double xtol = 1.0e-8;
    const double gtol = 1.0e-8;
    const double ftol = 0.0;
    const size_t n = fdf->n;
    const size_t p = fdf->p;

    double t[N], y[N], weights[N];
    
  
    gsl_multifit_nlinear_workspace* w = gsl_multifit_nlinear_alloc(T, fdf_params, n, p);


    gsl_vector* f= gsl_multifit_nlinear_residual(w);
    x = gsl_multifit_nlinear_position(w);
    gsl_matrix* J;
    gsl_matrix* covar = gsl_matrix_alloc(p, p);
  
    double x_init[3] = { 1.0, 1.0, 0.0 }; /* starting values */

    gsl_vector_view x0 = gsl_vector_view_array(x_init, p);
    gsl_vector_view wts = gsl_vector_view_array(weights, n);
    
    
    
   
    int status,info;
    double chisq0, chisq, rcond;
    gsl_rng* r;
    size_t i;
    gsl_rng_env_setup();
    r = gsl_rng_alloc(gsl_rng_default);

    /* this is the data to be fitted */
    for (i = 0; i < n; i++)
    {
        double ti = i * TMAX / (n - 1.0);
        double yi = 1.0 + 5 * exp(-1.5 * ti);
        double si = 0.1 * yi;
        double dy = gsl_ran_gaussian(r, si);

        t[i] = ti;
        y[i] = yi + dy;
        weights[i] = 1.0 / (si * si);
        printf("data: %g %g %g\n", ti, y[i], si);
    };


    /* initialize solver */
  
    gsl_multifit_nlinear_winit(&x0.vector, &wts.vector, fdf, w);

    /* compute initial cost function */
    gsl_blas_ddot(f, f, &chisq0);

    /* solve the system with a maximum of 100 iterations */
    status = gsl_multifit_nlinear_driver(max_iter, xtol, gtol, ftol,callback, NULL, &info, w);

    /* compute covariance of best fit parameters */
    J = gsl_multifit_nlinear_jac(w);
    gsl_multifit_nlinear_covar(J, 0.0, covar);

    /* compute final cost */
    gsl_blas_ddot(f, f, &chisq);

  #define FIT(i) gsl_vector_get(w->x, i)
  #define ERR(i) sqrt(gsl_matrix_get(covar,i,i))

    fprintf(stderr, "summary from method '%s/%s'\n",
        gsl_multifit_nlinear_name(w),
        gsl_multifit_nlinear_trs_name(w));
    fprintf(stderr, "number of iterations: %zu\n",
        gsl_multifit_nlinear_niter(w));
    fprintf(stderr, "function evaluations: %zu\n", fdf->nevalf);
    fprintf(stderr, "Jacobian evaluations: %zu\n", fdf->nevaldf);
    fprintf(stderr, "reason for stopping: %s\n",
        (info == 1) ? "small step size" : "small gradient");
    fprintf(stderr, "initial |f(x)| = %f\n", sqrt(chisq0));
    fprintf(stderr, "final   |f(x)| = %f\n", sqrt(chisq));

    {
        double dof = n - p;
        double c = GSL_MAX_DBL(1, sqrt(chisq / dof));

        fprintf(stderr, "chisq/dof = %g\n", chisq / dof);

        fprintf(stderr, "A      = %.5f +/- %.5f\n", FIT(0), c * ERR(0));
        fprintf(stderr, "lambda = %.5f +/- %.5f\n", FIT(1), c * ERR(1));
        fprintf(stderr, "b      = %.5f +/- %.5f\n", FIT(2), c * ERR(2));
    }

    fprintf(stderr, "status = %s\n", gsl_strerror(status));

    gsl_multifit_nlinear_free(w);
    gsl_matrix_free(covar);
 
    gsl_rng_free(r);
   
}


int
main(void)
{
    
    size_t i;
    const size_t n = N;
    const size_t p = 3;

    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();
    
    double t[N], y[N];
    struct data d = { n, t, y };


    

    /* define the function to be minimized */
    fdf.f = expb_f;
    fdf.df = expb_df;   /* set to NULL for finite-difference Jacobian */
    fdf.fvv = NULL;     /* not using geodesic acceleration */
    fdf.n = n;
    fdf.p = p;
    fdf.params = &d;

   
   

    fprintf(stderr, "=== Solving system with lm ===\n");
    fdf_params.trs = gsl_multifit_nlinear_trs_lm;
    solve_system(x, &fdf, &fdf_params);

    fprintf(stderr, "=== Solving system with dogleg ===\n");
    fdf_params.trs = gsl_multifit_nlinear_trs_dogleg;
    solve_system(x, &fdf, &fdf_params);

 
    
    gsl_vector_free(f);
    gsl_vector_free(x);

    return 0;
}

核心问题与修复方案 #

• 数据传递错误：原代码中solve_system内部定义了局部t、y数组，但fdf结构体关联的是main中的数组，导致拟合时使用未初始化的垃圾数据，直接引发NaN/Inf。修复时将数据生成逻辑移到main函数中，确保拟合使用正确的数据集。
• 参数指针覆盖：solve_system中的x = gsl_multifit_nlinear_position(w);完全多余，且会覆盖传入的指针，直接删除该行即可，求解器内部会维护参数向量。
• 重复数据生成：原代码每次调用solve_system都会重新生成随机数据，导致两次拟合的数据集不一致。修复后仅在main中生成一次数据，两次拟合复用同一组数据。

修复后代码 #

/* parameters to model */
#define N      100    /* number of data points to fit */
#define TMAX   (3.0)  /* time variable in [0,TMAX] */

struct data {
    size_t n;
    double* t;
    double* y;
};

int
expb_f(const gsl_vector* x, void* data,
    gsl_vector* f)
{
    size_t n = ((struct data*)data)->n;
    double* t = ((struct data*)data)->t;
    double* y = ((struct data*)data)->y;

    double A = gsl_vector_get(x, 0);
    double lambda = gsl_vector_get(x, 1);
    double b = gsl_vector_get(x, 2);

    size_t i;

    for (i = 0; i < n; i++)
    {
        /* Model Yi = A * exp(-lambda * t_i) + b */
        double Yi = A * exp(-lambda * t[i]) + b;
        gsl_vector_set(f, i, Yi - y[i]);
    }

    return GSL_SUCCESS;
}

int
expb_df(const gsl_vector* x, void* data,
    gsl_matrix* J)
{
    size_t n = ((struct data*)data)->n;
    double* t = ((struct data*)data)->t;

    double A = gsl_vector_get(x, 0);
    double lambda = gsl_vector_get(x, 1);

    size_t i;

    for (i = 0; i < n; i++)
    {
        /* Jacobian matrix J(i,j) = dfi / dxj, */
        /* where fi = (Yi - yi)/sigma[i],      */
        /*       Yi = A * exp(-lambda * t_i) + b  */
        /* and the xj are the parameters (A,lambda,b) */
        double e = exp(-lambda * t[i]);
        gsl_matrix_set(J, i, 0, e);
        gsl_matrix_set(J, i, 1, -t[i] * A * e);
        gsl_matrix_set(J, i, 2, 1.0);
    }

    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 rcond;

    /* compute reciprocal condition number of J(x) */
    gsl_multifit_nlinear_rcond(&rcond, w);

    fprintf(stderr, "iter %2zu: A = %.4f, lambda = %.4f, b = %.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),
        1.0 / rcond,
        gsl_blas_dnrm2(f));
}

void solve_system(gsl_vector* x, gsl_multifit_nlinear_fdf* fdf, gsl_multifit_nlinear_parameters* fdf_params, double weights[])
{
    const gsl_multifit_nlinear_type* T = gsl_multifit_nlinear_trust;
    const size_t max_iter = 100;
    const double xtol = 1.0e-8;
    const double gtol = 1.0e-8;
    const double ftol = 0.0;
    const size_t n = fdf->n;
    const size_t p = fdf->p;

    gsl_multifit_nlinear_workspace* w = gsl_multifit_nlinear_alloc(T, fdf_params, n, p);
    gsl_vector* f = gsl_multifit_nlinear_residual(w);
    gsl_matrix* J;
    gsl_matrix* covar = gsl_matrix_alloc(p, p);
  
    double x_init[3] = { 1.0, 1.0, 0.0 }; /* starting values */
    gsl_vector_view x0 = gsl_vector_view_array(x_init, p);
    gsl_vector_view wts = gsl_vector_view_array(weights, n);
    
    int status, info;
    double chisq0, chisq, rcond;

    /* initialize solver */
    gsl_multifit_nlinear_winit(&x0.vector, &wts.vector, fdf, w);

    /* compute initial cost function */
    gsl_blas_ddot(f, f, &chisq0);

    /* solve the system with a maximum of 100 iterations */
    status = gsl_multifit_nlinear_driver(max_iter, xtol, gtol, ftol, callback, NULL, &info, w);

    /* compute covariance of best fit parameters */
    J = gsl_multifit_nlinear_jac(w);
    gsl_multifit_nlinear_covar(J, 0.0, covar);

    /* compute final cost */
    gsl_blas_ddot(f, f, &chisq);

#define FIT(i) gsl_vector_get(w->x, i)
#define ERR(i) sqrt(gsl_matrix_get(covar,i,i))

    fprintf(stderr, "summary from method '%s/%s'\n",
        gsl_multifit_nlinear_name(w),
        gsl_multifit_nlinear_trs_name(w));
    fprintf(stderr, "number of iterations: %zu\n",
        gsl_multifit_nlinear_niter(w));
    fprintf(stderr, "function evaluations: %zu\n", fdf->nevalf);
    fprintf(stderr, "Jacobian evaluations: %zu\n", fdf->nevaldf);
    fprintf(stderr, "reason for stopping: %s\n",
        (info == 1) ? "small step size" : "small gradient");
    fprintf(stderr, "initial |f(x)| = %f\n", sqrt(chisq0));
    fprintf(stderr, "final   |f(x)| = %f\n", sqrt(chisq));

    {
        double dof = n - p;
        double c = GSL_MAX_DBL(1, sqrt(chisq / dof));

        fprintf(stderr, "chisq/dof = %g\n", chisq / dof);
        fprintf(stderr, "A      = %.5f +/- %.5f\n", FIT(0), c * ERR(0));
        fprintf(stderr, "lambda = %.5f +/- %.5f\n", FIT(1), c * ERR(1));
        fprintf(stderr, "b      = %.5f +/- %.5f\n", FIT(2), c * ERR(2));
    }

    fprintf(stderr, "status = %s\n", gsl_strerror(status));

    gsl_multifit