基于c语言的非线性最小二乘拟合

实现非线性最小二乘拟合可以使用 Levenberg-Marquardt 算法。该算法采用逐步调整步长的方式，在保证收敛的同时，能够保持算法的速度和精度。以下是一个基于C语言的非线性最小二乘拟合的示例代码：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define MAXN 100
#define MAXM 100

typedef struct {
    int n;    // 数据点个数
    int m;    // 待拟合参数个数
    double x[MAXN];    // 数据点x坐标
    double y[MAXN];    // 数据点y坐标
    double p[MAXM];    // 待拟合参数
    double **jacobi;   // 雅可比矩阵
    double *residuals; // 残差向量
} lm_data_type;

typedef void (*lm_func)(lm_data_type *, double *);

static void lm_mat_print(double **mat, int m, int n)
{
    int i, j;
    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++)
            printf("%lf ", mat[i][j]);
        printf("\n");
    }
}

static double **lm_mat_create(int m, int n)
{
    double **mat = (double **)malloc(m * sizeof(double *));
    int i;
    for (i = 0; i < m; i++)
        mat[i] = (double *)malloc(n * sizeof(double));
    return mat;
}

static void lm_mat_free(double **mat, int m)
{
    int i;
    for (i = 0; i < m; i++)
        free(mat[i]);
    free(mat);
}

static void lm_copy(double *dst, double *src, int n)
{
    int i;
    for (i = 0; i < n; i++)
        dst[i] = src[i];
}

static void lm_mat_copy(double **dst, double **src, int m, int n)
{
    int i, j;
    for (i = 0; i < m; i++)
        for (j = 0; j < n; j++)
            dst[i][j] = src[i][j];
}

static double lm_norm(double *v, int n)
{
    int i;
    double sum = 0.0;
    for (i = 0; i < n; i++)
        sum += (v[i] * v[i]);
    return sqrt(sum);
}

static void lm_func_wrapper(lm_data_type *data, double *fvec, double *p, lm_func f)
{
    int i;
    lm_copy(data->p, p, data->m);
    f(data, fvec);
}

static void lm_numerical_jacobian(lm_data_type *data, double *fvec, lm_func f)
{
    int i, j;
    double eps = 1e-8;
    double **jacobi = data->jacobi;
    double *residuals = data->residuals;
    double tmp;

    for (j = 0; j < data->m; j++) {
        double tmp_p = data->p[j];
        data->p[j] = tmp_p + eps;
        f(data, fvec);
        for (i = 0; i < data->n; i++)
            jacobi[i][j] = (fvec[i] - residuals[i]) / eps;
        data->p[j] = tmp_p;
    }
}

static void lm_levenberg_marquardt(lm_data_type *data, lm_func f)
{
    int i, j, k;
    int n = data->n;
    int m = data->m;
    double eps1 = 1e-8;
    double eps2 = 1e-8;
    double tau = 1e-3;
    double lambda = 0.001;
    double **jacobi = data->jacobi;
    double *residuals = data->residuals;
    double *fvec = (double *)malloc(n * sizeof(double));
    double *p = (double *)malloc(m * sizeof(double));
    double *q = (double *)malloc(m * sizeof(double));
    double *g = (double *)malloc(m * sizeof(double));
    double *s = (double *)malloc(m * sizeof(double));
    double *h = (double *)malloc(m * m * sizeof(double));
    double rho, nu, phi, psi;

    lm_copy(p, data->p, m);
    f(data, fvec);
    lm_copy(residuals, fvec, n);

    while (1) {
        lm_numerical_jacobian(data, fvec, f);
        for (i = 0; i < n; i++)
            residuals[i] = data->y[i] - fvec[i];

        phi = lm_norm(residuals, n);
        if (phi < eps1)
            break;

        for (i = 0; i < m; i++)
            g[i] = 0.0;
        for (i = 0; i < n; i++)
            for (j = 0; j < m; j++)
                g[j] += jacobi[i][j] * residuals[i];

        for (i = 0; i < m; i++) {
            for (j = 0; j < m; j++)
                h[i * m + j] = 0.0;
            h[i * m + i] = 1.0;
        }
        for (i = 0; i < m; i++)
            for (j = 0; j < n; j++)
                for (k = 0; k < m; k++)
                    h[i * m + k] += jacobi[j][i] * jacobi[j][k];

        while (1) {
            for (i = 0; i < m; i++)
                for (j = 0; j < m; j++)
                    s[i] = h[i * m + j] * g[j];
            nu = lm_norm(s, m);
            if (nu < eps2) {
                lm_copy(data->p, p, m);
                lambda *= tau;
                break;
            }

            for (i = 0; i < m; i++)
                q[i] = p[i] + s[i];
            lm_copy(data->p, q, m);
            f(data, fvec);

            lm_copy(residuals, fvec, n);
            psi = lm_norm(residuals, n);
            rho = (phi - psi) / (nu * nu);

            if (rho > 0.0) {
                lm_copy(p, q, m);
                lm_copy(residuals, fvec, n);
                phi = psi;

                if (rho < 0.25)
                    lambda *= 4.0;
                else if (rho > 0.75)
                    lambda *= 0.5;
            } else {
                lambda *= tau;
            }

            if (lambda < 1e-10)
                break;
        }
    }

    lm_mat_free(jacobi, n);
    free(fvec);
    free(p);
    free(q);
    free(g);
    free(s);
    free(h);
}

static void lm_example_func(lm_data_type *data, double *fvec)
{
    int i;
    double a = data->p[0];
    double b = data->p[1];
    double c = data->p[2];

    for (i = 0; i < data->n; i++)
        fvec[i] = a * exp(-b * data->x[i]) + c - data->y[i];
}

int main()
{
    int i;
    lm_data_type data;
    data.n = 10;
    data.m = 3;
    data.x[0] = 0.0;
    data.y[0] = 2.0;
    data.x[1] = 1.0;
    data.y[1] = 2.5;
    data.x[2] = 2.0;
    data.y[2] = 3.0;
    data.x[3] = 3.0;
    data.y[3] = 4.0;
    data.x[4] = 4.0;
    data.y[4] = 5.0;
    data.x[5] = 5.0;
    data.y[5] = 6.0;
    data.x[6] = 6.0;
    data.y[6] = 7.0;
    data.x[7] = 7.0;
    data.y[7] = 8.0;
    data.x[8] = 8.0;
    data.y[8] = 9.0;
    data.x[9] = 9.0;
    data.y[9] = 10.0;
    data.p[0] = 1.0;
    data.p[1] = 1.0;
    data.p[2] = 1.0;
    data.jacobi = lm_mat_create(data.n, data.m);
    data.residuals = (double *)malloc(data.n * sizeof(double));
    lm_levenberg_marquardt(&data, lm_example_func);

    printf("Fitted parameters:\n");
    for (i = 0; i < data.m; i++)
        printf("%lf\n", data.p[i]);

    lm_mat_free(data.jacobi, data.n);
    free(data.residuals);
    return 0;
}

在这个示例代码中，我们使用了一个示例函数 `lm_example_func`，该函数是 $y=a\exp(-bx)+c$，其中 $a$、$b$ 和 $c$ 是待拟合参数。在实际应用中，我们需要将 `lm_example_func` 替换为我们自己的目标函数。同时，我们还需要将数据点的 x 和 y 坐标、待拟合参数的个数和初值作为参数传递给 `lm_data_type` 结构体。最后，我们调用 `lm_levenberg_marquardt` 函数即可完成拟合。