最小二乘法的C语言实现

最小二乘法是一种数学优化技术，常用于拟合数据点，使得实际值与预测值之间的误差平方和达到最小。在C语言中，你可以通过迭代或者使用一些线性代数库（如BLAS或LAPACK）来实现。这里是一个简单的线性方程组求解的例子，假设我们有线性模型 y = A * x + b，其中A是矩阵，x是未知向量，y是已知向量，我们想要找到最小误差的x。

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

// 矩阵乘法函数
void matrix_multiply(int a_rows, int a_cols, double* a, int b_cols, double* b, double* result) {
    for (int i = 0; i < a_rows; ++i) {
        for (int j = 0; j < b_cols; ++j) {
            result[i * b_cols + j] = 0.0;
            for (int k = 0; k < a_cols; ++k)
                result[i * b_cols + j] += a[i * a_cols + k] * b[k * b_cols + j];
        }
    }
}

// 求解最小二乘法
double* least_squares(double** a, int m, double* b, int n) {
    // 初始化A'A
    double det = 0.0;
    double** AtA = malloc(sizeof(double*) * n);
    for (int i = 0; i < n; ++i) {
        AtA[i] = malloc(sizeof(double) * n);
        det += pow(a[i], 2); // 假设第一个行是单位向量
    }

    // 计算系数
    for (int i = 1; i < m; ++i) {
        for (int j = 0; j < n; ++j) {
            double sum = 0.0;
            for (int k = 0; k < n; ++k)
                sum += a[i][k] * AtA[j][k];
            AtA[j][n - 1] += a[i][n - 1] * sum;
            for (int k = 0; k < n; ++k)
                AtA[j][k] -= a[i][n - 1] * a[i][k];
        }
        det += pow(a[i][n - 1], 2);
    }

    // 计算系数除以det
    double* x = malloc(sizeof(double) * n);
    for (int i = 0; i < n; ++i) {
        double invDet = 1 / det;
        x[i] = (AtA[i][n - 1] * invDet);
        for (int j = 0; j < n; ++j)
            x[i] -= a[i][j] * x[j] * invDet;
    }

    return x;
}

int main() {
    // 假设输入数据和矩阵尺寸
    int m = 4, n = 3;
    double a[][] = { {1, 2, 3}, {4, 5, 6}, {7, 8, 9}, {1, 0, 1} };
    double b[] = {1, 2, 3, 4};

    // 解决最小二乘问题并打印结果
    double* x = least_squares((double**)a, m, b, n);
    printf("最小二乘解: ");
    for (int i = 0; i < n; ++i)
        printf("%f ", x[i]);
    free(x);
    return 0;
}