C语言实现最小二乘多项式曲线拟合（附带源码）

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

一、项目背景与目标说明

尽管多项式拟合的数学表达形式较为简洁，但在实际编程实现过程中仍需关注多个关键问题：数值稳定性、多项式阶数的选择、过拟合现象（高阶多项式容易在数据点之间产生剧烈振荡）、以及条件数恶化等问题。特别是构造 Vandermonde 矩阵时，该矩阵通常具有病态特性，容易导致求解结果不稳定。

本项目旨在开发一个基于 C 语言的多项式拟合程序，具备以下特点：

• 代码结构清晰，易于理解与运行，适合教学使用
• 支持任意阶数 m 的拟合（建议 m 远小于数据点数量 n）
• 采用正规方程法，并结合带部分主元的高斯消去法进行求解，兼顾通用性与可读性
• 能够计算出拟合系数，并对新的输入值进行预测
• 提供拟合优度（R）和残差相关统计信息
• 支持从文件或标准输入读取数据，附带示例程序以便直接测试验证

二、核心功能需求解析

三、关键技术点分析（数学原理与数值方法）

实现过程中涉及的主要数学与数值计算要点包括：

• 构建设计矩阵并形成正规方程组
• 利用最小二乘法求解超定线性系统
• 使用带部分主元的高斯消去法提升求解精度与稳定性
• 评估模型性能，计算 R 值、SSE（误差平方和）、SST（总平方和）等指标

四、模块化实现方案设计

整个程序采用模块化架构，各功能独立封装，便于维护与理解。主要模块如下：

1. 数据读取模块

负责从指定文件或标准输入中读取成对的 double 类型数据 (x, y)，存储于动态分配的数组中，并返回实际读取的数据点数量 n。

read_data()

2. 设计矩阵与正规方程构造模块

根据输入数据和设定的多项式阶数 m，构建用于最小二乘求解的设计矩阵 A 和观测向量 b。具体计算：

A = X^TX，其中 A 为 (m+1)×(m+1) 的对称矩阵

b = X^Ty，b 为长度为 m+1 的列向量

最终返回矩阵 A 和向量 b。

build_normal_equations()

A = V^T V

b = V^T y

3. 线性方程组求解模块

采用带部分主元的高斯消去法求解正规方程 A·coef = b，使用双精度浮点运算以提高数值稳定性。求解完成后通过回代得到拟合系数向量 coef。

gauss_solve()

4. 拟合评估与统计分析模块

包含以下功能：

• 根据拟合系数计算任意 x 对应的预测值 p(x)
• 计算 SSE（残差平方和）、SST（总偏差平方和）
• 计算决定系数 R = 1 - SSE/SST，用于衡量拟合优度

evaluate_poly()

compute_statistics()

5. 用户交互模块

解析命令行参数，控制程序流程；输出拟合结果；接收用户输入的 x 值并返回对应的预测 y 值。

6. 错误处理机制

在检测到矩阵接近奇异、条件数过大或设置过高阶数时，给出提示信息并安全退出，或建议降低拟合阶数以避免不稳定的计算结果。

性能优化与实现细节

• 使用 malloc 和 free 动态管理矩阵与向量内存
• 计算 x^k 时采用递推方式（即逐步乘积）而非幂函数调用，提升效率
• 利用 A 矩阵的对称性，在构造时仅计算上三角部分后镜像填充，节省约一半计算量（当前版本为增强可读性仍采用全填充方式）
• 高斯消去过程中，推荐使用行指针交换代替元素逐个交换以提高效率，但出于教学目的采用直观的元素交换方式，并在注释中提供指针优化建议

pow

A

五、完整代码实现

整合上述所有模块，形成可编译运行的完整 C 程序，包含详细注释，支持命令行调用与数据测试。

/******************************************************************************
 * 文件：poly_fit_least_squares.c
 * 功能：最小二乘多项式曲线拟合（Least Squares Polynomial Fitting）
 * 编译：gcc -std=c99 poly_fit_least_squares.c -o polyfit -lm
 *
 * 使用说明：
 *   ./polyfit data.txt m
 *   其中 data.txt 包含若干行 "x y"（空格或 tab 分隔）
 *   m 为拟合多项式的阶数（非负整数）
 *
 *   若没有提供文件，则从标准输入读取，输入结束可用 Ctrl+D (Linux) / Ctrl+Z (Windows)。
 *
 * 输出：
 *   - 拟合系数 a0..am
 *   - SSE、SST、R^2
 *   - 可交互输入 x 值得到预测 p(x)
 *
 * 备注：
 *   本实现采用正规方程 A a = b，A = V^T V，使用带部分主元的高斯消去求解。
 *   对高阶拟合或病态数据，建议使用 QR 分解或 SVD（文章最后给出扩展建议）。
 ******************************************************************************/

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

/* 判定小主元为零的阈值（相对） */
#define EPS_ZERO 1e-14

/* 动态数组装载数据点 */
typedef struct {
    int n;
    double *x;
    double *y;
} DataSet;

/* 读取数据：从文件或 stdin 读取，格式每行 "x y" */
DataSet read_data_from_file(const char *filename) {
    DataSet ds = {0, NULL, NULL};
    FILE *fp = NULL;
    if (filename) {
        fp = fopen(filename, "r");
        if (!fp) {
            fprintf(stderr, "无法打开文件：%s\n", filename);
            return ds;
        }
    } else {
        fp = stdin;
        printf("从标准输入读取数据，输入每行: x y，输入结束请输入 EOF（Ctrl+D/Ctrl+Z）\n");
    }

    int capacity = 128;
    ds.x = (double*)malloc(sizeof(double) * capacity);
    ds.y = (double*)malloc(sizeof(double) * capacity);
    if (!ds.x || !ds.y) {
        fprintf(stderr, "内存分配失败\n");
        if (fp != stdin) fclose(fp);
        free(ds.x); free(ds.y);
        ds.x = ds.y = NULL;
        ds.n = 0;
        return ds;
    }

    ds.n = 0;
    double xv, yv;
    while (fscanf(fp, "%lf %lf", &xv, &yv) == 2) {
        if (ds.n >= capacity) {
            capacity *= 2;
            ds.x = (double*)realloc(ds.x, sizeof(double) * capacity);
            ds.y = (double*)realloc(ds.y, sizeof(double) * capacity);
            if (!ds.x || !ds.y) {
                fprintf(stderr, "内存扩展失败\n");
                break;
            }
        }
        ds.x[ds.n] = xv;
        ds.y[ds.n] = yv;
        ds.n++;
    }

    if (fp != stdin) fclose(fp);
    return ds;
}

/* 构造正规方程 A = V^T V 和 b = V^T y
 * size = m+1
 * A is allocated as size*size (row-major)
 * b is allocated as size
 * 返回指针，调用者负责 free
 */
void build_normal_equations(const DataSet *ds, int m, double **A_out, double **b_out) {
    int n = ds->n;
    int sz = m + 1;

    double *A = (double*)calloc(sz * sz, sizeof(double));
    double *b = (double*)calloc(sz, sizeof(double));
    if (!A || !b) {
        fprintf(stderr, "内存分配失败（正规方程）\n");
        free(A); free(b);
        *A_out = NULL; *b_out = NULL;
        return;
    }

    /* 计算 sums: for powers of x up to 2m we may reuse them */
    /* pow_x[k] = sum_i x_i^k */
    int maxpow = 2 * m;
    double *pow_x = (double*)calloc((maxpow + 1), sizeof(double));
    if (!pow_x) {
        fprintf(stderr, "内存分配失败（pow_x）\n");
        free(A); free(b);
        *A_out = NULL; *b_out = NULL;
        return;
    }
    for (int k = 0; k <= maxpow; k++) pow_x[k] = 0.0;

    for (int i = 0; i < n; i++) {
        double xi = ds->x[i];
        double xp = 1.0;
        for (int k = 0; k <= maxpow; k++) {
            pow_x[k] += xp;    /* accumulate xi^k */
            xp *= xi;
        }
    }

    /* 构造 A: A[p,q] = sum_i x_i^{p+q} = pow_x[p+q] */
    for (int p = 0; p < sz; p++) {
        for (int q = 0; q < sz; q++) {
            A[p * sz + q] = pow_x[p + q];
        }
    }

    /* 构造 b: b[p] = sum_i y_i * x_i^p */
    for (int p = 0; p < sz; p++) b[p] = 0.0;
    for (int i = 0; i < n; i++) {
        double xi = ds->x[i];
        double yi = ds->y[i];
        double xp = 1.0;
        for (int p = 0; p < sz; p++) {
            b[p] += yi * xp;
            xp *= xi;
        }
    }

    free(pow_x);
    *A_out = A;
    *b_out = b;
}

/* 带部分主元的高斯消去 + 回代
 * 求解 Ax = b，A 为 n x n，b 长度 n
 * A 和 b 会被破坏
 * 返回 0 成功, 非 0 表示失败（奇异）
 */
int gaussian_solve(double *A, double *b, double *x, int n) {
    /* A: row-major n*n */
    for (int k = 0; k < n; k++) {
        /* 寻找第 k 列中从 k..n-1 的最大主元行 */
        int maxRow = k;
        double maxVal = fabs(A[k * n + k]);
        for (int i = k + 1; i < n; i++) {
            double v = fabs(A[i * n + k]);
            if (v > maxVal) {
                maxVal = v; maxRow = i;
            }
        }
        if (maxVal < EPS_ZERO) {
            /* 主元过小，矩阵可能奇异或病态 */
            return -1;
        }
        /* 交换行 k 与 maxRow（包括 b） */
        if (maxRow != k) {
            for (int j = k; j < n; j++) {
                double tmp = A[k * n + j];
                A[k * n + j] = A[maxRow * n + j];
                A[maxRow * n + j] = tmp;
            }
            double tmpb = b[k];
            b[k] = b[maxRow];
            b[maxRow] = tmpb;
        }

        /* 消元 */
        for (int i = k + 1; i < n; i++) {
            double factor = A[i * n + k] / A[k * n + k];
            /* A[i][k] = 0 */
            A[i * n + k] = 0.0;
            for (int j = k + 1; j < n; j++) {
                A[i * n + j] -= factor * A[k * n + j];
            }
            b[i] -= factor * b[k];
        }
    }

    /* 回代求解 */
    for (int i = n - 1; i >= 0; i--) {
        double sum = b[i];
        for (int j = i + 1; j < n; j++) sum -= A[i * n + j] * x[j];
        if (fabs(A[i * n + i]) < EPS_ZERO) return -2;
        x[i] = sum / A[i * n + i];
    }
    return 0;
}

/* 评估多项式 p(x) 给定系数 coef[0..m] */
double evaluate_poly(const double *coef, int m, double x) {
    /* Horner 方法 */
    double res = coef[m];
    for (int k = m - 1; k >= 0; k--) res = res * x + coef[k];
    return res;
}

/* 计算 SSE、SST、R^2 */
void compute_statistics(const DataSet *ds, const double *coef, int m, double *SSE, double *SST, double *R2) {
    int n = ds->n;
    double mean_y = 0.0;
    for (int i = 0; i < n; i++) mean_y += ds->y[i];
    mean_y /= n;

    double sse = 0.0, sst = 0.0;
    for (int i = 0; i < n; i++) {
        double yi = ds->y[i];
        double yhat = evaluate_poly(coef, m, ds->x[i]);
        double e = yi - yhat;
        sse += e * e;
        double d = yi - mean_y;
        sst += d * d;
    }
    *SSE = sse; *SST = sst;
    if (sst <= 0.0) *R2 = 1.0; /* 常数序列 */
    else *R2 = 1.0 - sse / sst;
}

/* 打印多项式系数为可读形式 */
void print_polynomial(const double *coef, int m) {
    printf("拟合多项式（阶数 %d）：\n  p(x) = ", m);
    for (int k = 0; k <= m; k++) {
        printf("%.10g", coef[k]);
        if (k >= 1) printf(" * x^%d", k);
        if (k < m) printf(" + ");
    }
    printf("\n");
}

/* 主程序 */
int main(int argc, char *argv[]) {
    const char *filename = NULL;
    int m = 2; /* 默认二次拟合 */
    if (argc >= 2) {
        filename = argv[1];
    }
    if (argc >= 3) {
        m = atoi(argv[2]);
        if (m < 0) m = 0;
    }

    DataSet ds = read_data_from_file(filename);
    if (ds.n <= 0) {
        fprintf(stderr, "没有读取到数据点，程序退出\n");
        return 1;
    }
    if (ds.n <= m) {
        fprintf(stderr, "数据点数量 n=%d 必须大于拟合阶数 m=%d（需 n > m）\n", ds.n, m);
        free(ds.x); free(ds.y);
        return 1;
    }

    /* 构造正规方程 */
    double *A = NULL;
    double *b = NULL;
    build_normal_equations(&ds, m, &A, &b);
    if (!A || !b) {
        fprintf(stderr, "构造正规方程失败\n");
        free(ds.x); free(ds.y);
        return 1;
    }

    int sz = m + 1;
    double *coef = (double*)calloc(sz, sizeof(double));
    if (!coef) {
        fprintf(stderr, "系数数组分配失败\n");
        free(A); free(b); free(ds.x); free(ds.y);
        return 1;
    }

    /* 复制 A 和 b，因为高斯消去会破坏它们；若需要保留可复制 */
    double *A_copy = (double*)malloc(sizeof(double) * sz * sz);
    double *b_copy = (double*)malloc(sizeof(double) * sz);
    if (!A_copy || !b_copy) {
        fprintf(stderr, "内存分配失败\n");
        free(A); free(b); free(coef); free(ds.x); free(ds.y);
        free(A_copy); free(b_copy);
        return 1;
    }
    memcpy(A_copy, A, sizeof(double) * sz * sz);
    memcpy(b_copy, b, sizeof(double) * sz);

    int status = gaussian_solve(A_copy, b_copy, coef, sz);
    if (status != 0) {
        fprintf(stderr, "高斯求解失败（矩阵可能奇异或病态），status=%d\n", status);
        fprintf(stderr, "建议降低拟合阶数或使用 QR/SVD 方法。\n");
        free(A); free(b); free(A_copy); free(b_copy); free(coef);
        free(ds.x); free(ds.y);
        return 1;
    }

    /* 输出系数与统计量 */
    print_polynomial(coef, m);
    double SSE, SST, R2;
    compute_statistics(&ds, coef, m, &SSE, &SST, &R2);
    printf("样本数 n = %d, 多项式阶数 m = %d\n", ds.n, m);
    printf("残差平方和 SSE = %.12g\n", SSE);
    printf("总平方和 SST = %.12g\n", SST);
    printf("决定系数 R^2 = %.12g\n", R2);

    /* 交互预测 */
    printf("\n请输入若干 x 值（按回车），输入非数字或 EOF 结束预测：\n");
    double xt;
    while (scanf("%lf", &xt) == 1) {
        double ypred = evaluate_poly(coef, m, xt);
        printf("p(%.8g) = %.12g\n", xt, ypred);
    }

    /* 清理 */
    free(A); free(b); free(A_copy); free(b_copy); free(coef);
    free(ds.x); free(ds.y);
    return 0;
}

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：曲线拟合 最小二乘 多项式 C语言 Statistics

相关内容：C语言源码实现