最小二乘法的实现

---

问题描述

给定一组基函数$f_0(x),f_1(x),\dots,f_{k-1}(x)$和一组数据点$(x_i,y_i)$，求系数$c_1,c_2,\dots,c_k$使得函数$f(x)=\sum c_if_i(x)$对于数据点的拟合程度最高，即最小化：

$$\sum_i(f(x_i)-y(i))^2$$

分析

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

$$\begin{equation} A=\left( \begin{array}{cccc} f_0(x_0) & f_1(x_0) & \dots &f_{k-1}(x_0)\\ f_1(x_1) & f_2(x_1) & \dots &f_{k-1}(x_1)\\ \vdots & \vdots & \dots & \vdots & \\ f_1(x_{n-1}) & f_2(x_{n-1}) & \dots &f_{k-1}(x_{n-1})\\ \end{array} \right) \\ \textbf{c}=\left( \begin{array}{c} c_0\\c_1\\\vdots\\c_{k-1} \end{array} \right), \textbf{y} = \left( \begin{array}{c} y_0\\y_1\\\vdots\\y_{n-1} \end{array} \right) \end{equation}$$

则误差矩阵：

$$\eta = A\textbf{c}-\textbf{y}$$

最小化：

$$\sum_{i}\eta_i^2=\sum_i (\sum_kA_{ki}\textbf{c}_i-\textbf{y}_i)^2$$

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

$$A^T(A\textbf{c}-\textbf{y})=\textbf{0}$$

解得：

$$\textbf c=(A^TA)^{-1}A^T\textbf y$$

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

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;
}