最小二乘法的正交多项拟合算法

最小二乘法的正交多项拟合算法

数学表达式

参考清华大学 李庆杨《数值分析》第五版 P69

构造正交多项式 ${p_k(x)}_{k=0}^n$
$$\begin{gathered} p_0(x)=1 \\ {p}_1(x)=(x-A_1)p_0(x) \\ ... \\ {p}_{k+1}(x)=(x-A_{k+1})p_k(x)-B_k{p}_{k-1}(x) \\ k=1,2,3,...,n-1 \end{gathered}$$
其中
$$\begin{gathered} A_{k+1}=\frac{(x{p}_k,p_k)}{(p_k,p_k)}=\frac{{\sum }_{i=0}^m{x}_i{p}_k^2(x_i)}{{\sum }_{i=0}^m{p}_k^2(x_i)},k=0,1,2,...,n-1 \\ {B}_k=\frac{(p_k,p_k)}{(p_{k-1},p_{k-1})}=\frac{{\sum }_{i=0}^m{p}_k^2(x_i)}{{\sum }_{i=0}^m{p}_{k-1}^2(x_i)},k=1,2,...,n-1 \end{gathered}$$
利用正交多项式构建 $a^{\ast }$
$$a_k^{\ast }=\frac{(f,p_k)}{(p_k,p_k)}=\frac{{\sum }_{i=0}^{m}f(x_i)p_k(x_i)}{{\sum }_{i=0}^m{p}_k^2(x_i)}$$
则待拟合曲线的正交多项式为
$$S^{\ast }=a_0^{\ast }{p}_0(x)+a_1^{\ast }{p}_1(x)+...+a_n^{\ast }{p}_n(x)$$

计算示例

给定数据表，求二次最小二乘拟合多项式

$x_i$	0	0.5	0.6	0.7	0.8	0.9	1.0
$y_i$	1.0	1.75	1.96	2.19	2.44	2.71	3.0

1. 求第一项 $a_0^{\ast }$
   $$a_0^{\ast }=\frac{(f,p_0)}{(p_0,p_0)}=\frac{{\sum }_{i=0}^{6}f(x_i)p_0(x_i)}{{\sum }_{i=0}^6{p}_0^2(x_i)}$$

$p_0(x_i)=1;$

${\sum }_0^{6}f(x_i)=15.05;$

求得 $a_0^{\ast }=15.05/7=2.15$

2. 求第二项 $a_1^{\ast }$
   $$a_1^{\ast }=\frac{(f,p_1)}{(p_1,p_1)}=\frac{{\sum }_{i=0}^{6}f(x_i)p_1(x_i)}{{\sum }_{i=0}^6{p}_1^2(x_i)}$$

先确定 $p_1(x_i)=x-A_1$， 求 $A_1$
$$A_1=\frac{(x{p}_0,p_0)}{(p_0,p_0)}=\frac{{\sum }_{i=0}^6{x}_i}{{\sum }_{i=0}^6{p}_0\ast p_0}=\frac{4.5}{7}=0.64$$
则 $p_1(x_i)=x-0.64$，求得
$$a_1^{\ast }=\frac{1.30}{0.657}\approx 1.98$$

3. 求第三项 $a_2^{\ast }$
   $$a_2^{\ast }=\frac{(f,p_2)}{(p_2,p_2)}=\frac{{\sum }_{i=0}^{6}f(x_i)p_2(x_i)}{{\sum }_{i=0}^6{p}_2^2(x_i)}$$
   $p_2(x)=(x-A_2)p_1(x)-B_1{p}_0(x)=(x-A_2)(x-A_1)-B_1$
   $$A_2=\frac{(x{p}_1,p_1)}{(p_1,p_1)}=\frac{{\sum }_{i=0}^6{x}_i\ast (x_i-0.64{)}^2}{{\sum }_{i=0}^6(x_i-0.64{)}^2}=\frac{0.2242}{0.6572}=0.34$$

$$B_1=\frac{(p_1,p_1)}{(p_0,p_0)}=\frac{0.6572}{7}=0.094$$

$p_2(x)=x^2-0.98x+0.12$

求得【此处需要高精度计算】
$$a_2^{\ast }\approx \frac{0.0686}{0.0686}\approx 1$$
最终
$$S^{\ast }=a_0^{\ast }{p}_0(x)+a_1^{\ast }{p}_1(x)+a_2^{\ast }{p}_2(x)=x^2+x+1$$

正交多项式是目前为止多项式拟合最好的方法，不用解线性方程组，根据次数 n，只用递推公式方便地计算正交多项式；MATLAB 正交多项式最小二乘拟合函数 p = polyfit (x,y,n)，输入参数 $x=[x_0,x_1,...x_m];y=[y_0,y_1,...,y_m]$，返回次数为 n 的多项式 p(x) 的系数，p 中的系数按降幂排列，p 的长度为 n + 1。

Python 代码

以拟合线性或二次曲线为例，由于1-2次曲线，最多只有 $a_0^{\ast }..a_2^{\ast },p_0..p_2$ 需要计算，代码量比较小

def linearfit(xl: list, yl: list)->list:
    if len(xl) != len(yl) or len(xl) < 2:
        return []

    # 生成必要数据
    A = 0.0
    asl = [0, 0]

    # 定义p函数
    p = [lambda n:1, lambda n: n-A]
    # 定义x，y映射函数
    def f(i): return yl[xl.index(i)]

    # 计算第一项
    asl[0] = sum(map(f, xl))/sum(map(p[0], xl))

    # 计算第二项
    A = sum(map(lambda n: n, xl)) / sum(map(p[0], xl)) #确定p1函数的项
    asl[1] = sum(map(lambda n: f(n)*p[1](n), xl)) / sum(map(lambda n:p[1](n)**2, xl))

    # 计算返回系数
    cof = [asl[0]-A*asl[1],asl[1]]
    return cof


def QuadFit(xl: list, yl: list) -> list:
    if len(xl) != len(yl) or len(xl) < 3:
        return []

    # 定义x，y映射函数
    def f(i): return yl[xl.index(i)]

    # 生成必要数据
    A = [0.0, 0.0, 0.0]
    B = [0.0, 0.0, 0.0]
    asl = [0, 0, 0]
    # 定义p函数
    p = [lambda n:1, lambda n: n-A[1], lambda n:((n-A[2])*(n-A[1])-B[1])]

    # 计算第一项
    asl[0] = sum(map(f, xl))/sum(map(p[0], xl))

    # 计算第二项
    A[1] = sum(map(lambda n: n, xl)) / sum(map(p[0], xl))
    asl[1] = sum(map(lambda n: f(n)*p[1](n), xl)) / \
        sum(map(lambda n: p[1](n)**2, xl))

    # 计算第三项
    A[2] = sum(map(lambda n: n*p[1](n)**2, xl)) / \
        sum(map(lambda n: p[1](n)**2, xl))
    B[1] = sum(map(lambda n: p[1](n)**2, xl)) / \
        sum(map(lambda n: p[0](n)**2, xl))
    asl[2] = sum(map(lambda n: f(n)*p[2](n), xl)) / \
        sum(map(lambda n: p[2](n)**2, xl))

    def rst(n): return asl[0] + asl[1] * (n-A[1]) + \
        asl[2]*((n-A[2])*(n-A[1])-B[1])

    # 计算返回系数
    cof = [rst(0), asl[1]-asl[2]*(A[1]+A[2]), asl[2]]
    return cof