1. 多项式拟合问题

  多项式拟合（polynominal curve fitting）是一种线性模型，模型和拟合参数的关系是线性的。多项式拟合的输入是一维的，即$\textbf{x} = x$，这是多项式拟合和线性回归问题的主要区别之一。

  多项式拟合的目标是构造输入$x$的$M$阶多项式函数，使得该多项式能够近似表示输入$x$和输出$y$的关系，虽然实际上$x$和$y$的关系并不一定是多项式，但使用足够多的阶数，总是可以逼近表示输入$x$和输出$y$的关系的。

  多项式拟合问题的输入可以表示如下：

$$\begin{equation*} D = \{(x_1,y_1),(x_2,y_2),...,(x_i,y_i),...,(x_N,y_N)\} \\ x_i \in R \\ y_i \in R \end{equation*}$$

  目标输出是得到一个多项式函数：

$$\begin{equation*} \begin{aligned} f(x) &= w_1 x^1 + w_2 x^2 + w_i x^i + ... + w_Mx^M + b\\ &= (\sum_{i=1}^{M} w_i x^i) + b \end{aligned} \end{equation*}$$

其中$M$表示最高阶数为$M$。

  可见在线性拟合的模型中，共包括了$(M+1)$个参数，而该模型虽然不是输入$x$的线性函数，但却是$(M+1)$个拟合参数的线性函数，所以称多项式拟合为线性模型。对于多项式拟合问题，其实就是要确定这$(M+1)$个参数，这里先假设阶数$M$是固定的（$M$是一个超参数，可以用验证集来确定$M$最优的值，详细的关于$M$值确定的问题，后面再讨论），重点就在于如何求出这$(M+1)$个参数的值。

2.优化目标

  多项式拟合是利用多项式函数逼近输入$x$和输出$y$的函数关系，通过什么指标来衡量某个多项式函数的逼近程度呢？（其实这就是误差/损失函数）。拟合/回归问题常用的评价指标是均方误差（在机器学习中的模型评估与度量博客中，我进行了介绍）。多项式拟合问题也同样采用该评价指标，以均方误差作为误差/损失函数，误差函数越小，模型越好。

$$\begin{equation*} E(\mathbf {w}, b) = \frac{1}{N} \sum_{i=1}^{N} {\lbrack f(x_i) - y_i\rbrack}^2 \end{equation*}$$

  系数$\frac{1}{N}$是一常数，对优化结果无影响，可以去除，即将均方误差替换为平方误差：

$$\begin{equation*} E(\mathbf{w}, b) = \sum_{i=1}^{N} {\lbrack f(x_i) - y_i\rbrack}^2 \end{equation*}$$
   到这里，就成功把多项式拟合问题变成了最优化问题，优化问题可表示为：

$$\begin{equation*} \mathop {\arg \min} _ {\mathbf{w} , b} E(\mathbf{w},b) \end{equation*}$$

即需要求得参数$\{w_1, ...,w_M, b\}$的值，使得$E(\mathbf{w}, b)$最小化。那么如何对该最优化问题求解呢？

3. 优化问题求解

3.1 求偏导，联立方程求解

   直观的想法是，直接对所有参数求偏导，令偏导为0，再联立这$M+1$个方程求解（因为共有$M+1$个参数，故求偏导后也是得到$M+1$个方程）。

$$\begin{equation*} \begin{aligned} E(\mathbf{w}, b) &= \sum_{i=1}^{N} {\lbrack f(x_i) - y_i\rbrack}^2 \\ &= \sum_{i=1}^{N} {\lbrack (w_1 x_i^1 + w_2 x_i^2 + w_i x_i^j + ... + w_Mx_i^M + b) - y_i\rbrack}^2 \end{aligned} \end{equation*}$$

利用$E(\mathbf{w},b)$对各个参数求偏导，如下：

$$\begin{equation*} \begin{aligned} \frac{\partial E(\mathbf{w}, b)}{\partial w_j} &= 2 \sum_{i=1}^{N} {\lbrack (w_1 x_i^1 + w_2 x_i^2 + w_i x_i^j + ... + w_Mx_i^M + b) - y_i\rbrack} x_i^j \\ \frac{\partial E(\mathbf{w}, b)}{\partial b} &= 2 \sum_{i=1}^{N} {\lbrack (w_1 x_i^1 + w_2 x_i^2 + w_i x_i^j + ... + w_Mx_i^M + b) - y_i\rbrack} \end{aligned} \end{equation*}$$

求导之后，将各个点$(x_i,y_i)$的值带入偏导公式，联立方程求解即可。

  针对该解法，可以举个例子详细说明，比如有两个点$(2,3),(5,8)$,需要利用二阶多项式$f(x) = w_1x + w_2x^2 + b$拟合。求解过程如下：

1. 该二阶多项式对参数求偏导得到
   

   $$\begin{equation*} \begin{aligned} \frac{\partial E(\mathbf{w}, b)}{\partial w_j} &= 2 \sum_{i=1}^{2} {\lbrack (w_1 x_i^1 + w_2 x_i^2 + b) - y_i\rbrack} x_i^j \\ &= [(w_1x_1 + w_2x_1^2 + b) - y_1] x_1^j + [(w_1x_2 + w_2x_2^2 + b) - y_2] x_2^j\\ \frac{\partial E(\mathbf{w}, b)}{\partial b} &= 2 \sum_{i=1}^{2} {\lbrack (w_1 x_i^1 + w_2 x_i^2 + b) - y_i\rbrack} \\ &= [(w_1x_1 + w_2x_1^2 + b) - y_1] + [(w_1x_2 + w_2x_2^2 + b) - y_2]\\ \end{aligned} \end{equation*}$$

2. 将点$(2,3),(5,8)$带入方程，可以得到3个方程，
   

   $$\begin{equation*} \begin{aligned} 2b + 7w_1 + 29w_2= 11 \\ 7b + 29w_1 + 133 w_2 = 46 \\ 29b + 133w_1 + 641w_2 = 212 \end{aligned} \end{equation*}$$

3. 联立这三个方程求解，发现有无穷多的解，只能得到$3w_1 + 21w_2 = 5$，这三个方程是线性相关的，故没有唯一解。

  该方法通过求偏导，再联立方程求解，比较复杂，看着也很不美观。那么有没有更加方便的方法呢？

3.2 最小二乘法

   其实求解该最优化问题（平方和的最小值）一般会采用最小二乘法（其实最小二乘法和求偏导再联立方程求解的方法无本质区别，求偏导也是最小二乘法，只是这里介绍最小二乘的矩阵形式而已）。最小二乘法（least squares），从英文名非常容易想到，该方法就是求解平方和的最小值的方法。

  可以将误差函数以矩阵的表示($N$个点，最高$M$阶)为：

$$\begin{equation*} {\lVert \mathbf{Xw - y} \rVert}_2 \end{equation*}$$

其中，把偏置$b$融合到了参数$\bf w$中，

$$\begin{equation*} \mathbf{w} = \{b, w_1, w_2, ..., w_M\} \end{equation*}$$

$\mathbf X$则表示输入矩阵，

$$\begin{gather*} \begin{bmatrix} 1 & x_1 & x_1^2 &... &x_1^M \\ 1 & x_2 & x_2^2 & ... & x_2^M \\ ... & ... & ... & ... & ... \\ 1 & x_N & x_N^2 & ... & x_N^M \\ \end{bmatrix} \end{gather*}$$

$\mathbf y$则表示标注向量，

$$\begin{equation*} \mathbf{y} = \{y_1,y_2,...,y_N\}^T \end{equation*}$$

因此，最优化问题可以重新表示为

$$\begin{equation*} \min_{\mathbf w} {\lVert \mathbf{Xw - y} \rVert}_2 \end{equation*}$$

对其求导，

$$\begin{equation*} \begin{aligned} \frac {\partial {\lVert \mathbf{Xw - y} \rVert}_2}{\partial \mathbf{w}} &= \frac{\partial (\mathbf{Xw} - \mathbf{y})^T (\mathbf{Xw} - \mathbf{y})}{\partial \mathbf{w}} \\ &= \frac{\partial (\mathbf{w}^T\mathbf{X}^T - \mathbf{y}^T) (\mathbf{Xw} - \mathbf{y})}{\partial \mathbf{w}}\\ &= \frac{\partial (\mathbf{w}^T\mathbf{X}^T\mathbf{Xw} - \mathbf{y}^T\mathbf{Xw} - \mathbf{w}^T\mathbf{X}^T\mathbf{y} + \mathbf{y}^T\mathbf{y}) }{\partial \mathbf{w}} \\ \end{aligned} \end{equation*}$$

在继续对其求导之前，需要先补充一些矩阵求导的先验知识（常见的一些矩阵求导公式可以参见转载的博客https://blog.csdn.net/lipengcn/article/details/52815429），如下：

$$\begin{equation*} \frac{\partial \mathbf{x}^T\mathbf{a}}{\partial \mathbf{x}} = \mathbf{a} \\ \frac{\partial \mathbf{ax}}{\partial \mathbf{x}} = \mathbf{a}^T \\ \frac{\partial \mathbf{x}^T\mathbf{A}}{\partial \mathbf{x}} = \mathbf{Ax} + \mathbf{A}^T\mathbf{x} \end{equation*}$$

根据上面的矩阵求导规则，继续进行损失函数的求导

$$\begin{equation*} \begin{aligned} \frac {\partial {\lVert \mathbf{Xw - y} \rVert}_2}{\partial \mathbf{w}} &= \frac{\partial (\mathbf{w}^T\mathbf{X}^T\mathbf{Xw} - \mathbf{y}^T\mathbf{Xw} - \mathbf{w}^T\mathbf{X}^T\mathbf{y} + \mathbf{y}^T\mathbf{y}) }{\partial \mathbf{w}} \\ &= \mathbf{X}^T\mathbf{Xw} + (\mathbf{X}^T\mathbf{X})^T\mathbf{w} - (\mathbf{y}^T\mathbf{X})^T - \mathbf{X}^T\mathbf{y} \\ &= 2 \mathbf{X}^T\mathbf{Xw} - 2\mathbf{X}^T\mathbf{y} \end{aligned} \end{equation*}$$

其中$\mathbf{X}^T\mathbf{Xw} = (\mathbf{X}^T\mathbf{X})^T\mathbf{w}$.令求导结果等于0，即可以求导问题的最小值。

$$\begin{equation*} \begin{aligned} 2 \mathbf{X}^T\mathbf{Xw} - 2\mathbf{X}^T\mathbf{y} = 0 \\ \mathbf{w} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y} \end{aligned} \end{equation*}$$

  再利用最小二乘法的矩阵形式对前面的例子进行求解，用二阶多项式拟合即两个点$(2,3),(5,8)$。

1. 表示输入矩阵 $\mathbf{X}$和标签向量$\mathbf{y}$
   

   $$\begin{gather*} \mathbf{X} = \begin{bmatrix} 1 & 2& 4 \\ 1 & 5 & 25\\ \end{bmatrix} \\ \mathbf{y} = \begin{bmatrix} 3 & 8 \\ \end{bmatrix} ^T\\ \end{gather*}$$

2. 计算$\mathbf{X}^T\mathbf{X}$
   

   $$\begin{gather*} \mathbf{X}^T\mathbf{X} = \begin{bmatrix} 2 & 7& 29 \\ 7 & 29 & 133\\ 29 & 133 & 641\\ \end{bmatrix} \\ \end{gather*}$$

3. 矩阵求逆，再做矩阵乘法运算
   但 $\mathbf{X}^T\mathbf{X}$不可逆，故无唯一解。

  关于矩阵的逆是否存在，可以通过判断矩阵的行列式是否为0（$det(\mathbf{A}) \stackrel{?}{=} 0$ 来判断，也可以通过初等行变换，观察矩阵的行向量是否线性相关，在这个例子下，矩阵不可逆，故有无穷多解。但如果新增一个点$(4,7)$，则就可以解了。

  其实这和数据集的点数和选择的阶数有关，如果点数小于阶数则会出现无穷解的情况，如果点数等于阶数，那么刚好有解可以完全拟合所有数据点，如果点数大于阶数，则会求的近似解。

  那么对于点数小于阶数的情况，如何求解？在python的多项式拟合函数中是可以拟合的，而且效果不错，具体算法不是很了解，可以想办法参考python的ployfit()函数的实现。

4. 拟合阶数的选择

   在前面的推导中，多项式的阶数被固定了，那么实际场景下应该如何选择合适的阶数$M$呢？

1. 一般会选择阶数$M$小于点数$N$
2. 把训练数据分为训练集合验证集，在训练集上，同时用不同的$M$值训练多个模型，然后选择在验证集误差最小的阶数 M M <script type="math/tex" id="MathJax-Element-5573">M</script>

   5. 后续

     如果后续还想写的话，可以考虑正则化问题。

多项式拟合C#

Class1.cs-拟合类

using System;
using System.Collections.Generic;
using System.Text;

namespace 最小二乘法拟合多项式
{
    class Class1
    {
        ///<summary>
        ///用最小二乘法拟合二元多次曲线
        ///</summary>
        ///<param name="arrX">已知点的x坐标集合</param>
        ///<param name="arrY">已知点的y坐标集合</param>
        ///<param name="length">已知点的个数</param>
        ///<param name="dimension">方程的最高次数</param>
        public double[] MultiLine(double[] arrX, double[] arrY, int length, int dimension)//二元多次线性方程拟合曲线
        {
            int n = dimension + 1;                  //dimension次方程需要求 dimension+1个 系数
            double[,] Guass = new double[n, n + 1];      //高斯矩阵 例如：y=a0+a1*x+a2*x*x
            for (int i = 0; i < n; i++)
            {
                int j;
                for (j = 0; j < n; j++)
                {
                    Guass[i, j] = SumArr(arrX, j + i, length);
                }
                Guass[i, j] = SumArr(arrX, i, arrY, 1, length);
            }
            return ComputGauss(Guass, n);
        }

        public double SumArr(double[] arr, int n, int length) //求数组的元素的n次方的和
        {
            double s = 0;
            for (int i = 0; i < length; i++)
            {
                if (arr[i] != 0 || n != 0)
                    s = s + Math.Pow(arr[i], n);
                else
                    s = s + 1;
            }
            return s;
        }

        public double SumArr(double[] arr1, int n1, double[] arr2, int n2, int length)
        {
            double s = 0;
            for (int i = 0; i < length; i++)
            {
                if ((arr1[i] != 0 || n1 != 0) && (arr2[i] != 0 || n2 != 0))
                    s = s + Math.Pow(arr1[i], n1) * Math.Pow(arr2[i], n2);
                else
                    s = s + 1;
            }
            return s;

        }

        public double[] ComputGauss(double[,] Guass, int n)
        {
            int i, j;
            int k, m;
            double temp;
            double max;
            double s;
            double[] x = new double[n];
            for (i = 0; i < n; i++) x[i] = 0.0;//初始化

            for (j = 0; j < n; j++)
            {
                max = 0;
                k = j;
                for (i = j; i < n; i++)
                {
                    if (Math.Abs(Guass[i, j]) > max)
                    {
                        max = Guass[i, j];
                        k = i;
                    }
                }


                if (k != j)
                {
                    for (m = j; m < n + 1; m++)
                    {
                        temp = Guass[j, m];
                        Guass[j, m] = Guass[k, m];
                        Guass[k, m] = temp;
                    }
                }
                if (0 == max)
                {
                    // "此线性方程为奇异线性方程" 
                    return x;
                }

                for (i = j + 1; i < n; i++)
                {
                    s = Guass[i, j];
                    for (m = j; m < n + 1; m++)
                    {
                        Guass[i, m] = Guass[i, m] - Guass[j, m] * s / (Guass[j, j]);
                    }
                }

            }//结束for (j=0;j<n;j++)

            for (i = n - 1; i >= 0; i--)
            {
                s = 0;
                for (j = i + 1; j < n; j++)
                {
                    s = s + Guass[i, j] * x[j];
                }
                x[i] = (Guass[i, n] - s) / Guass[i, i];
            }
            return x;
        }//返回值是函数的系数
        //例如：y=a0+a1*x 返回值则为a0 a1
        //例如：y=a0+a1*x+a2*x*x 返回值则为a0 a1 a2
    }
}

Progran.cs - 测试类

using System;
using System.Collections.Generic;
using System.Text;

namespace 最小二乘法拟合多项式
{
    class Program
    {
        static void Main(string[] args)
        {
            double[] xCoe;

            double[] xLeastValue = {
0.005,
1.01,
2.019,
3.025,
4.027,
5.025,
6.016,
7.007,
7.999,
8.998,
10.006};
            double[] zLeastValue = {
0,
0.25,
0.5,
0.75,
1,
1.25,
1.5,
1.75,
2,
2.25,
2.5};
            //例如：y=a0+a1*x+a2*x*x 返回值则为a0 a1 a2
            Class1 sb = new Class1();
            xCoe = sb.MultiLine(xLeastValue, zLeastValue, 11, 5);
            Console.WriteLine("5阶多项式因子：");
            for(int i = 0;i < xCoe.Length; i++)
            {
                Console.WriteLine("a"+i+":"+xCoe[i]);
            }

            Console.WriteLine("多项式计算出来的位移值 标准位移值 两者差值 线性度：");
            for (int i = 0; i < xLeastValue.Length; i++)
            {
                double Y = xCoe[0] + xCoe[1] * Math.Pow(xLeastValue[i], 1) + xCoe[2] * Math.Pow(xLeastValue[i], 2) + xCoe[3] * Math.Pow(xLeastValue[i], 3)
                    + xCoe[4] * Math.Pow(xLeastValue[i], 4) + xCoe[5] * Math.Pow(xLeastValue[i], 5);
                double dev = Y - zLeastValue[i];
                double error = dev / 2.5;
                Console.WriteLine(Y + "   " + zLeastValue[i] + "  " + dev + "    " + error);
            }

             Console.ReadKey(true);
        }


    }
}

多项式拟合之一般方程法

之前经常会用到多项式拟合来拟合函数关系，十分好用，得出的表达式使用起来超级方便。在此对多项式拟合原理进行分析。

1 什么是多项式拟合？

$$\begin{array}{rl}y& =a_0+a_{1}x+a_2{x}^2+\cdots +a_n{x}^n\\ & =\sum _{i=0}^n{a}_i{x}^i\end{array}$$

这是多项式拟合的输出表达式。

2 怎样拟合？

考虑输入样本数据如下：
$$\begin{gathered} x=[x_1,x_2,...,x_m{]}^T \\ y=[y_1,y_2,...,y_m{]}^T \end{gathered}$$
构建方程组：
$$\begin{array}{rl}{y}_1& =a_0+a_1{x}_1+a_2{x}_1^2+\cdots +a_n{x}_1^n\\ y_2& =a_0+a_1{x}_2+a_2{x}_2^2+\cdots +a_n{x}_2^n\\ ⋮& \\ y_m& =a_0+a_1{x}_m+a_2{x}_m^2+\cdots +a_n{x}_m^n\end{array}$$

令：

则上面的方程组可以写为：
$$Ab=f$$
拟合过程就是求$b$的过程

3 求$b$

现在一个很直观的想法，就是：
$$\begin{array}{r}Ab=f\\ b=A^{-1}f\end{array}$$
但是，没这么简单。考虑一下矩阵维度
$$A_{[m,n+1]}{b}_{[n+1,1]}=f_{[m,1]}$$
强行规定输入的$x=[x_1,x_2,...,x_m]$不能重复，则$A$的各个列向量都线性无关

• 当$m=n+1$时，$A$是可逆矩阵，可以直接计算

• 当$m<n+1$时，方程数少于未知数个数，方程不可解，拟合失败

• 当$m>n+1$时，怎么算呢？

  正交矩阵

  对于一个矩阵A，如果有$A{A}^T=E$，则称之为正交矩阵

  QR分解

  将一个矩阵A，分解为一个正交矩阵和一个非奇异上三角阵的过程
  $$A_{[m,n+1]}=Q_{[m,n+1]}{R}_{[n+1,n+1]}$$
  其中$Q$是正交矩阵，$R$是非奇异上三角阵

  那么，如何进行QR分解呢？

  施密特正交化

  这个请大家自行学习。

  完成QR分解后
  $$\begin{array}{rl}Ab& =f\\ QRb& =f\\ Rb& =Q^{T}f\\ b& =R^{-1}{Q}^{t}f\end{array}$$

matlab测试函数

测试数据

%% 准备测试数据
x=1:0.1:50;
y=0.001*x.^4+100*rand(size(x));

plot(x,y,'b')
title('测试数据')

拟合与测试

%% 拟合
n=4;
c = mypolyfit(x',y',n);
%% 测试效果
y1 = mypolyval(c,x');
plot(x,y,'b',x,y1','r')
title('测试拟合效果')

mypolyfit函数

function c = mypolyfit(x,f,n)
A = x.^(0:n);
c = A\f;

mypolyval函数

function y = mypolyval(c,s)
n = length(c)-1;
B = s.^(0:n);
y = B*c;