C++基于GSL实现带参数固定的双自变量非线性最小二乘拟合

原框架核心逻辑无需大幅调整，仅需修改数据存储、残差计算、对外接口三个模块即可支持双独立变量，原有的仿函数固定参数能力完全保留，无需改动原有求解器核心逻辑。

具体修改点 #

• 修改fit_data结构体：原结构仅存储1组自变量和观测值，新增第二个独立变量的常引用成员，同时调整长度校验逻辑，保证两组自变量长度与观测值长度完全一致。
• 修改残差计算函数internal_f：遍历样本时同时读取两个独立变量的对应索引值，和自由参数一起传入待拟合函数计算预测值，再和观测值做差得到残差。
• 调整对外curve_fit接口：新增第二个独立变量向量的入参，初始化fit_data时同步传入两组自变量。

修改后完整核心代码 #

// 原gen_tuple、n_params、internal_solve_system、internal_make_gsl_vector_ptr逻辑完全保留，无需修改

// 调整fit_data结构体，支持双自变量存储
template<typename C1>
struct fit_data
{
    const std::vector<double>& x1; // 第一个独立变量
    const std::vector<double>& x2; // 第二个独立变量
    const std::vector<double>& y;  // 观测因变量
    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);
    // 从GSL向量提取待优化的自由参数
    auto init_args = [x](int index)
    {
        return gsl_vector_get(x, index);
    };
    auto parameters = gen_tuple<n_params>(init_args);

    // 逐样本计算残差
    for (size_t i = 0; i < d->y.size(); ++i)
    {
        double x1i = d->x1[i];
        double x2i = d->x2[i];
        double yi = d->y[i];
        auto func = [x1i, x2i, &d](auto ...xs)
        {
            // 传入两个独立变量+自由参数调用拟合函数
            return d->f(x1i, x2i, xs...);
        };
        double y_pred = std::apply(func, parameters);
        gsl_vector_set(f, i, yi - y_pred);
    }
    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<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.x1.size() == fd.y.size() && fd.x2.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.y.size();
    fdf.p   = initial_params->size;
    fdf.params = &fd;

    // 保留原Levenberg-Marquardt+几何加速度配置
    fdf_params.trs = gsl_multifit_nlinear_trs_lmaccel;
    return internal_solve_system(initial_params, &fdf, &fdf_params);
}

// 调整对外curve_fit接口，新增第二个自变量入参
template <typename Callable, auto n>
auto curve_fit(Callable f, const std::array<double, n>& initial_params, const std::vector<double>& x1, const std::vector<double>& x2, const std::vector<double>& y) -> std::vector<double>
{
    assert(initial_params.size() == n);
    auto params = internal_make_gsl_vector_ptr(initial_params);
    auto fd = fit_data<Callable>{x1, x2, y, f};
    return curve_fit_impl(internal_f<decltype(fd), n>, nullptr, nullptr, params,  fd);
}

固定参数仿函数适配方法 #

原仿函数固定参数的逻辑完全兼容，仅需按照双自变量签名定义仿函数即可。以固定参数a的高斯函数为例，仿函数写法调整为：

struct gaussian_fixed_a
{
    double a;
    gaussian_fixed_a(double a) : a{a} {}
    // 签名顺序：第一个自变量x，第二个自变量b，待优化自由参数c
    double operator()(double x, double b, double c) const { 
        const double z = (x - b) / c;
        return a * std::exp(-0.5 * z * z); 
    }
};

原有测试代码无需调整逻辑，直接调用即可，返回结果就是自由参数c的最优拟合值，拟合精度和单自变量场景一致。

如果后续需要支持3个及以上自变量，按照相同逻辑在fit_data中新增对应变量的引用、残差计算时读取对应值传入拟合函数、对外接口新增对应入参即可，框架扩展性不受影响。

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