如何在GSL非线性最小二乘拟合中限制参数搜索范围？

数据老陈

2026-8-23

为GSL非线性最小二乘拟合添加参数区间约束

我使用GSL的非线性最小二乘拟合例程实现自定义函数的拟合，已完成固定部分参数、估计其余参数的功能，现在需要给参数搜索范围添加区间限制，避免陷入参数空间错误区域的局部最优解。

以下是现有C++封装GSL的实现代码：

template <typename F, size_t... Is>
auto gen_tuple_impl(F func, std::index_sequence<Is...> )
{
    return std::make_tuple(func(Is)...);
}

template <size_t N, typename F>
auto gen_tuple(F func)
{
    return gen_tuple_impl(func, std::make_index_sequence<N>{} );
}

template <class R, class... ARGS>
struct function_ripper {
    static constexpr size_t n_args = sizeof...(ARGS);
};

template <class R, class... ARGS>
auto constexpr n_params(R (ARGS...) )
{
    return function_ripper<R, ARGS...>();
}


auto internal_solve_system(gsl_vector* initial_params, gsl_multifit_nlinear_fdf *fdf,
             gsl_multifit_nlinear_parameters *params) -> std::vector<double>
{
  // This specifies a trust region method
  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;

  auto *work = gsl_multifit_nlinear_alloc(T, params, fdf->n, fdf->p);
  int info;

  // initialize solver
  gsl_multifit_nlinear_init(initial_params, fdf, work);
  //iterate until convergence
  gsl_multifit_nlinear_driver(max_iter, xtol, gtol, ftol, nullptr, nullptr, &info, work);

  // result will be stored here
  gsl_vector * y    = gsl_multifit_nlinear_position(work);
  auto result = std::vector<double>(initial_params->size());

  for(int i = 0; i < result.size(); i++)
  {
    result[i] = gsl_vector_get(y, i);
  }

  auto niter = gsl_multifit_nlinear_niter(work);
  auto nfev  = fdf->nevalf;
  auto njev  = fdf->nevaldf;
  auto naev  = fdf->nevalfvv;

  // nfev - number of function evaluations
  // njev - number of Jacobian evaluations
  // naev - number of f_vv evaluations
  //logger::debug("curve fitted after ", niter, " iterations {nfev = ", nfev, "} {njev = ", njev, "} {naev = ", naev, "}");

  gsl_multifit_nlinear_free(work);
  gsl_vector_free(initial_params);
  return result;
}

template<auto n>
auto internal_make_gsl_vector_ptr(const std::array<double, n>& vec) -> gsl_vector*
{
    auto* result = gsl_vector_alloc(vec.size());
    int i = 0;
    for(const auto e: vec)
    {
        gsl_vector_set(result, i, e);
        i++;
    }
    return result;
}


template<typename C1>
struct fit_data
{
    const std::vector<double>& t;
    const std::vector<double>& y;
    // the actual function to be fitted
    C1 f;
};


template<typename FitData, int n_params>
int internal_f(const gsl_vector* x, void* params, gsl_vector *f)
{
    auto* d  = static_cast<FitData*>(params);
    // Convert the parameter values from gsl_vector (in x) into std::tuple
    auto init_args = [x](int index)
    {
        return gsl_vector_get(x, index);
    };
    auto parameters = gen_tuple<n_params>(init_args);

    // Calculate the error for each...
    for (size_t i = 0; i < d->t.size(); ++i)
    {
        double ti = d->t[i];
        double yi = d->y[i];
        auto func = [ti, &d](auto ...xs)
        {
            // call the actual function to be fitted
            return d->f(ti, xs...);
        };
        auto y = std::apply(func, parameters);
        gsl_vector_set(f, i, yi - y);
    }
    return GSL_SUCCESS;
}

using func_f_type   = int (*) (const gsl_vector*, void*, gsl_vector*);
using func_df_type  = int (*) (const gsl_vector*, void*, gsl_matrix*);
using func_fvv_type = int (*) (const gsl_vector*, const gsl_vector *, void *, gsl_vector *);

template<auto n>
auto internal_make_gsl_vector_ptr(const std::array<double, n>& vec) -> gsl_vector*;


auto internal_solve_system(gsl_vector* initial_params, gsl_multifit_nlinear_fdf *fdf,
             gsl_multifit_nlinear_parameters *params) -> std::vector<double>;

template<typename C1>
auto curve_fit_impl(func_f_type f, func_df_type df, func_fvv_type fvv, gsl_vector* initial_params, fit_data<C1>& fd) -> std::vector<double>
{
    assert(fd.t.size() == fd.y.size());

    auto fdf = gsl_multifit_nlinear_fdf();
    auto fdf_params = gsl_multifit_nlinear_default_parameters();

    fdf.f   = f;
    fdf.df  = df;
    fdf.fvv = fvv;
    fdf.n   = fd.t.size();
    fdf.p   = initial_params->size;
    fdf.params = &fd;

    // "This selects the Levenberg-Marquardt algorithm with geodesic acceleration."
    fdf_params.trs = gsl_multifit_nlinear_trs_lmaccel;
    return internal_solve_system(initial_params, &fdf, &fdf_params);
}


template <typename Callable, auto n>
auto curve_fit(Callable f, const std::array<double, n>& initial_params, const std::vector<double>& x, const std::vector<double>& y) -> std::vector<double>
{
    // We can't pass lambdas without convert to std::function.
    //constexpr auto n = 3;//decltype(n_params(f))::n_args - 5;
    //constexpr auto n = 2;
    assert(initial_params.size() == n);

    auto params = internal_make_gsl_vector_ptr(initial_params);
    auto fd = fit_data<Callable>{x, y, f};
    return curve_fit_impl(internal_f<decltype(fd), n>, nullptr, nullptr, params,  fd);
}

待拟合的自定义函数gaussian：

double gaussian(double x, double b, double a, double c)
{
    const double z = (x - b) / c;
    return a * std::exp(-0.5 * z * z);
}

struct gaussian_fixed_a
{
    double a;
    gaussian_fixed_a(double a) : a{a} {}
    double operator()(double x, double b, double c) const { return gaussian(x, b, a, c); }
};

测试用的数据集生成与拟合代码：

int main()
{
    auto device = std::random_device();
    auto gen    = std::mt19937(device());

    auto xs = linspace<std::vector<double>>(0.0, 1.0, 300);
    auto bs = linspace<std::vector<double>>(0.4, 1.4, 300);
    auto ys = std::vector<double>(xs.size());

    double a = 5.0, c = 0.15;

    for(size_t i = 0; i < xs.size(); i++)
    {
        auto y =  gaussian(xs[i], a, bs[i], c);
        auto dist  = std::normal_distribution(0.0, 0.1 * y);
        ys[i] = y + dist(gen);
    }
    gaussian_fixed_a g(a);
    auto r = curve_fit(g, std::array{0.11}, xs, bs, ys);

    std::cout << "result: " << r[0] << ' ' << '\n';
    std::cout << "error : " << r[0] - c << '\n';
}

我需要在数值优化的信赖域中定义参数边界，请问该如何实现？

解决方案

GSL的非线性最小二乘模块支持通过约束信赖域方法实现参数区间限制，核心是利用gsl_multifit_nlinear_constraint结构体定义参数上下边界，并将其关联到求解器工作区。以下是具体修改步骤：

1. 扩展拟合数据结构，添加边界存储

修改fit_data结构体，新增参数上下限的引用：

template<typename C1>
struct fit_data
{
    const std::vector<double>& t;
    const std::vector<double>& y;
    C1 f;
    // 新增：待估计参数的上下边界，长度与参数数量一致
    const std::vector<double>& param_lower;
    const std::vector<double>& param_upper;
};

2. 更新求解函数，添加约束逻辑

修改internal_solve_system函数，创建并初始化约束对象，使用带约束的工作区分配函数：

auto internal_solve_system(gsl_vector* initial_params, gsl_multifit_nlinear_fdf *fdf,
             gsl_multifit_nlinear_parameters *params, const std::vector<double>& lower, const std::vector<double>& upper) -> std::vector<double>
{
  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;

  // 创建约束结构体，为每个参数设置上下限
  gsl_multifit_nlinear_constraint *constr = gsl_multifit_nlinear_constraint_alloc(fdf->p);
  for (size_t i = 0; i < fdf->p; ++i) {
      gsl_multifit_nlinear_constraint_set(constr, i, lower[i], upper[i]);
  }

  // 分配带约束的工作区
  auto *work = gsl_multifit_nlinear_alloc_with_constraint(T, params, constr, fdf->n, fdf->p);
  int info;

  gsl_multifit_nlinear_init(initial_params, fdf, work);
  gsl_multifit_nlinear_driver(max_iter, xtol, gtol, ftol, nullptr, nullptr, &info, work);

  // 提取拟合结果
  gsl_vector * y    = gsl_multifit_nlinear_position(work);
  auto result = std::vector<double>(initial_params->size());
  for(int i = 0; i < result.size(); i++)
  {
    result[i] = gsl_vector_get(y, i);
  }

  auto niter = gsl_multifit_nlinear_niter(work);
  auto nfev  = fdf->nevalf;
  auto njev  = fdf->nevaldf;
  auto naev  = fdf->nevalfvv;

  // 释放资源
  gsl_multifit_nlinear_free(work);
  gsl_multifit_nlinear_constraint_free(constr);
  gsl_vector_free(initial_params);
  return result;
}

3. 修改拟合实现函数，传递边界参数

更新curve_fit_impl，将边界信息传递给求解函数：

template<typename C1>
auto curve_fit_impl(func_f_type f, func_df_type df, func_fvv_type fvv, gsl_vector* initial_params, fit_data<C1>& fd) -> std::vector<double>
{
    assert(fd.t.size() == fd.y.size());
    assert(fd.param_lower.size() == initial_params->size());
    assert(fd.param_upper.size() == initial_params->size());

    auto fdf = gsl_multifit_nlinear_fdf();
    auto fdf_params = gsl_multifit_nlinear_default_parameters();

    fdf.f   = f;
    fdf.df  = df;
    fdf.fvv = fvv;
    fdf.n   = fd.t.size();
    fdf.p   = initial_params->size;
    fdf.params = &fd;

    fdf_params.trs = gsl_multifit_nlinear_trs_lmaccel;
    // 传递边界给求解函数
    return internal_solve_system(initial_params, &fdf, &fdf_params, fd.param_lower, fd.param_upper);
}

4. 更新外部调用接口，支持传入边界

修改curve_fit模板函数，新增边界参数：

template <typename Callable, auto n>
auto curve_fit(Callable f, const std::array<double, n>& initial_params, 
               const std::vector<double>& x, const std::vector<double>& y,
               const std::vector<double>& param_lower, const std::vector<double>& param_upper) -> std::vector<double>
{
    assert(initial_params.size() == n);
    assert(param_lower.size() == n);
    assert(param_upper.size() == n);

    auto params = internal_make_gsl_vector_ptr(initial_params);
    auto fd = fit_data<Callable>{x, y, f, param_lower, param_upper};
    return curve_fit_impl(internal_f<decltype(fd), n>, nullptr, nullptr, params,  fd);
}

5. 在主函数中使用约束拟合

以高斯函数拟合为例，指定参数c的区间（如0.1 <= c <= 0.2）：

int main()
{
    auto device = std::random_device();
    auto gen    = std::mt19937(device());

    auto xs = linspace<std::vector<double>>(0.0, 1.0, 300);
    auto bs = linspace<std::vector<double>>(0.4, 1.4, 300);
    auto ys = std::vector<double>(xs.size());

    double a = 5.0, c = 0.15;

    for(size_t i = 0; i < xs.size(); i++)
    {
        auto y =  gaussian(xs[i], a, bs[i], c);
        auto dist  = std::normal_distribution(0.0, 0.1 * y);
        ys[i] = y + dist(gen);
    }
    gaussian_fixed_a g(a);
    
    // 定义参数c的上下边界
    std::vector<double> lower = {0.1};
    std::vector<double> upper = {0.2};
    auto r = curve_fit(g, std::array{0.11}, xs, bs, ys, lower, upper);

    std::cout << "result: " << r[0] << ' ' << '\n';
    std::cout << "error : " << r[0] - c << '\n';
}

关键注意事项

约束仅对待估计参数生效，固定参数无需设置边界
初始参数值必须落在指定区间内，否则求解器会自动将初始值投影到约束区间
约束信赖域方法会在迭代过程中强制参数保持在设定区间，避免进入无效参数区域

内容的提问来源于stack exchange，提问作者CaffèSospeso

如何在GSL非线性最小二乘拟合中限制参数搜索范围？

1. 扩展拟合数据结构，添加边界存储 #

2. 更新求解函数，添加约束逻辑 #

3. 修改拟合实现函数，传递边界参数 #

4. 更新外部调用接口，支持传入边界 #

5. 在主函数中使用约束拟合 #