一种高效稳定的多项式拟合算法的C语言实现

创作时间:

作者:

@小白创作中心

一种高效稳定的多项式拟合算法的C语言实现

引用

CSDN

https://blog.csdn.net/jz_ddk/article/details/145943315

多项式拟合是数据拟合和数值计算中的基础方法之一，广泛应用于工程、科研等领域。本文介绍了一种基于Cholesky分解的高效稳定多项式拟合算法，并提供了完整的C语言实现代码。该算法通过引入正则化项和Cholesky分解，有效提升了计算效率和数值稳定性，能够更好地处理病态问题。

关键实现原理

正则化处理：在法方程矩阵的对角线元素上添加λ（A[i][j] += lambda），通过吉洪诺夫正则化改善矩阵条件数
Cholesky分解优势：
时间复杂度，比高斯消元快一倍
数值稳定性更好，避免主元选取问题
天然保持矩阵的对称性

实现步骤

构造带正则项的法方程矩阵
进行Cholesky分解得到下三角矩阵L
通过前代和回代求解线性方程组

数值稳定性改进

正则项确保矩阵正定
使用平方根分解避免平方误差累积
严格的下三角计算顺序保证分解稳定性

使用建议

多项式次数不宜过高（一般≤10）
λ取值通过交叉验证选择
数据建议进行归一化处理
可通过计算残差平方和评估拟合效果

性能对比

该实现相比传统高斯消元法：

计算效率提高约50%
对病态问题的容忍度提升1-2个数量级
内存访问模式更规则，有利于缓存优化

完整实现代码

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

// 动态创建二维数组
double** create_matrix(int rows, int cols) {
    double** mat = (double**)malloc(rows * sizeof(double*));
    for (int i = 0; i < rows; i++) {
        mat[i] = (double*)calloc(cols, sizeof(double));
    }
    return mat;
}

// Cholesky分解
int cholesky(double** A, double** L, int n) {
    for (int i = 0; i < n; i++) {
        for (int j = 0; j <= i; j++) {
            double sum = 0.0;
            
            if (j == i) { // 对角线元素
                for (int k = 0; k < j; k++)
                    sum += L[j][k] * L[j][k];
                L[j][j] = sqrt(A[j][j] - sum);
                if (L[j][j] <= 0) return 0; // 矩阵不正定
            } else { // 非对角线元素
                for (int k = 0; k < j; k++)
                    sum += L[i][k] * L[j][k];
                L[i][j] = (A[i][j] - sum) / L[j][j];
            }
        }
    }
    return 1;
}

// 解下三角方程组 Ly = b
void forward_substitution(double** L, double* y, double* b, int n) {
    for (int i = 0; i < n; i++) {
        y[i] = b[i];
        for (int j = 0; j < i; j++) {
            y[i] -= L[i][j] * y[j];
        }
        y[i] /= L[i][i];
    }
}

// 解上三角方程组 L^T x = y
void backward_substitution(double** L, double* x, double* y, int n) {
    for (int i = n-1; i >= 0; i--) {
        x[i] = y[i];
        for (int j = i+1; j < n; j++) {
            x[i] -= L[j][i] * x[j]; // 注意L转置的索引
        }
        x[i] /= L[i][i];
    }
}

// 多项式拟合主函数
void poly_fit(double* x, double* y, int num_points, 
             int degree, double lambda, double* coeffs) {
    
    int n = degree + 1;
    double** A = create_matrix(n, n);
    double* b = (double*)calloc(n, sizeof(double));
    
    // 构建法方程矩阵（带正则项）
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            for (int k = 0; k < num_points; k++) {
                A[i][j] += pow(x[k], i + j);
            }
            if (i == j) A[i][j] += lambda; // 添加正则项
        }
    }
    
    // 构建右端向量
    for (int i = 0; i < n; i++) {
        for (int k = 0; k < num_points; k++) {
            b[i] += y[k] * pow(x[k], i);
        }
    }
    
    // Cholesky分解
    double** L = create_matrix(n, n);
    if (!cholesky(A, L, n)) {
        fprintf(stderr, "矩阵不正定，无法分解\n");
        exit(1);
    }
    
    // 解方程组
    double* y_vec = (double*)malloc(n * sizeof(double));
    forward_substitution(L, y_vec, b, n);
    backward_substitution(L, coeffs, y_vec, n);
    
    // 释放内存
    for (int i = 0; i < n; i++) free(A[i]);
    free(A);
    free(b);
    for (int i = 0; i < n; i++) free(L[i]);
    free(L);
    free(y_vec);
}

void polynomial_predict(double *x, double *y, unsigned int length , double *coeffs, int degree) {
    double result = 0;
    for (int j =0;j<length;j++) {
        result = 0;
        for (int i = 0; i <= degree; i++){
            result += coeffs[i] * pow(x[j], i);
        }
        y[j]=result;
    }
}

int main() {
    // 示例数据：带噪声的正弦曲线
    int num_points = 10;
    double x[num_points], y[num_points];
    for (int i = 0; i < num_points; i++) {
        x[i] = i;
        y[i] = sin(2*M_PI*0.5*x[i]/10) + 0.1 * ((double)rand()/RAND_MAX - 0.5);
    }
    
    int degree = 6;   // 多项式次数
    double lambda = 0.1; // 正则化参数
    double coeffs[degree+1];
    
    poly_fit(x, y, num_points, degree, lambda, coeffs);
    
    printf("拟合多项式系数：\n");
    for (int i = 0; i <= degree; i++) {
        printf("a%d = %f\n", i, coeffs[i]);
    }
    unsigned int predict_num = 15;
    double x1[predict_num],y1[predict_num];
    for (int i=0;i<predict_num;i++){
        x1[i] = i;
    }
    polynomial_predict(x1, y1, predict_num , coeffs, degree);
    //存储原始数据和外推数据结果
    FILE *src = fopen("source.txt", "w");
    for (int i=0;i<num_points;i++){
        fprintf(src,"%.16lf\n",y[i]);
    }
    fclose(src);
    FILE *dest = fopen("dest.txt", "w");
    for (int i=0;i<predict_num;i++){
        fprintf(dest,"%.16lf\n",y1[i]);
    }
    fclose(dest);
    
    return 0;
}