最小二乘法的实现

---

问题描述

给定一组基函数f0(x),f1(x),…,fk−1(x)和一组数据点(xi,yi)，求系数c1,c2,…,ck使得函数f(x)=∑cifi(x)对于数据点的拟合程度最高，即最小化：

∑i(f(xi)−y(i))2

分析

原问题可以表述成矩阵形式：

A=⎛⎝⎜⎜⎜⎜f0(x0)f1(x1)⋮f1(xn−1)f1(x0)f2(x1)⋮f2(xn−1)…………fk−1(x0)fk−1(x1)⋮fk−1(xn−1)⎞⎠⎟⎟⎟⎟c=⎛⎝⎜⎜⎜⎜c0c1⋮ck−1⎞⎠⎟⎟⎟⎟,y=⎛⎝⎜⎜⎜⎜y0y1⋮yn−1⎞⎠⎟⎟⎟⎟

则误差矩阵：

η=Ac−y

最小化：

∑iη2i=∑i(∑kAkici−yi)2

依次对cp,p=0,1,2求导并令导数为0，得到方程：

AT(Ac−y)=0

解得：

c=(ATA)−1ATy

如果选取幂函数作为基函数，则是用多项式拟合一个函数。

Code

#include <bits/stdc++.h>using namespace std;const int MAXN = 505;double A[MAXN][MAXN], AT[MAXN][MAXN];double B[MAXN][MAXN], M[MAXN][MAXN], R[MAXN][MAXN];void mul(double A[MAXN][MAXN], double B[MAXN][MAXN], int n, int m, int q){    memset(M, 0, sizeof M);    for (int i = 1; i <= n; i++)        for (int j = 1; j <= q; j++)             for (int k = 1; k <= m; k++)                M[i][j] += A[i][k]*B[k][j];}void T(double A[MAXN][MAXN], int n, int m){    for (int i = 1; i <= n; i++)        for (int j = 1; j <= m; j++)            M[j][i] = A[i][j];}void inverse(double A[MAXN][MAXN], int n){    memcpy(M, A, sizeof A), memset(R, 0, sizeof R);    for (int i = 1; i <= n; i++) R[i][i] = 1;    for (int i = 1; i <= n; i++) {        double mx = A[i][i];        int id = i;        for (int j = i+1; j <= n; j++)            if (A[j][i] > mx)                mx = A[j][i], id = i;        swap(A[i], A[id]);        for (int j = 1; j <= n; j++) {            if (j == i) continue;            double v = -M[j][i]/M[i][i];            for (int k = 1; k <= n; k++)                M[j][k] += M[i][k]*v, R[j][k] += R[i][k]*v;        }    }    for (int i = 1; i <= n; i++) {        for (int j = 1; j <= n; j++)            R[i][j] /= M[i][i];    }    memcpy(M, R, sizeof R);}double x[MAXN], y[MAXN], ATy[MAXN];int n, k;int main(){    scanf("%d%d", &n, &k);    for (int i = 1; i <= n; i++)        scanf("%lf%lf", &x[i], &y[i]);    for (int i = 1; i <= n; i++) {        A[i][1] = 1;        for (int j = 2; j <= k; j++)            A[i][j] = A[i][j-1]*x[i];    }    T(A, n, k);    memcpy(AT, M, sizeof M);    for (int i = 1; i <= k; i++) {        ATy[i] = 0;        for (int j = 1; j <= n; j++)            ATy[i] += AT[i][j]*y[j];    }    mul(AT, A, k, n, k);    memcpy(B, M, sizeof M);    inverse(B, k);    memcpy(B, M, sizeof M);    for (int i = 1; i <= k; i++) {        double c = 0;        for (int j = 1; j <= k; j++)            c += B[i][j]*ATy[j];        printf("%.3f\n", c);    }    return 0;}