最小二乘拟合问题求解算法（含python代码）

1. 最小二乘拟合原理

        参考西安交大数值分析教材

        当f(x)是定义在点集X上的列表函数时，内积为

 2. 最小二乘拟合问题的求解过程

最小二乘拟合问题的求解过程有三种方法

法一：

       （1）根据问题的特点选择一组线性无关的基函数，

（2）解方正规方程组

（3）求得拟合函数

法二：

（1）利用三项递推关系式构造正交多项式

（2）由式5.2.11得到最小二乘拟合多项式

法三：

（1）做变量替换，利用已有的正交基函数构造最优平方逼近函数

（2）转到法二求解

3. python程序以及说明

        本程序实现了上述的法一和法二，对于法一，程序选择了xk，(k=1，2，3 ……)作为基函数，然后使用改进平方根法求解正规方程组得到系数，最后利用系数构造拟合函数。对于法二，程序中使用最高项系数为1的正交多项式的三项递推关系构造正交基函数，然后通过式5.2.11求解系数。正交多项式的三项递推关系如下

        python代码如下

import numpy as np
from sympy.abc import t
from sympy import simplify

class LeastSquareFit():
    '''
    用于计算列表函数的最小二乘拟合
    输入参数：
    X：自变量列表
    Y：函数值列表
    n：多项式阶数
    JiFunction：基函数形式
    '''
    
    def __init__(self,X,Y,n,JiFunction='正交多项式'):
        self.x=np.array(X)
        self.y=np.array(Y)
        self.n=n
        self.w=self.myset_w(f=1) #设置权函数
        self.JiFunction=self.generate_JiFuction(n,JiFunction) #生成基函数
        self.c=[0]*(n+1)  #初始化拟合多项式的系数
    
    #计算权函数
    def myset_w(self,f):
        if type(f) is int or float:
            y=[f]*len(self.x)
        else:
            y=[f.subs(t,x) for x in self.x]
        return y

    #生成正交基函数
    def generate_JiFuction(self,m,function_name):
        if function_name == '正交多项式':
            value_num=m
            g={}
            g[0]=1
            beta0=np.dot([(t*g[0]).subs(t,x) for x in self.x] , np.array(self.w)*([g[0]]*len(self.x)))
            gama0=np.dot([g[0]]*len(self.x) , np.array(self.w)*([g[0]]*len(self.x)))
            g[1]=t-beta0/gama0
            for k in range(1,value_num):
                beta=np.dot([(t*g[k]).subs(t,x) for x in self.x] , np.array(self.w)*[g[k].subs(t,x) for x in self.x])
                gama=np.dot([g[k].subs(t,x) for x in self.x] , np.array(self.w)*[g[k].subs(t,x) for x in self.x])
                if k==1:
                    gama_pre=np.dot([g[k-1]]*len(self.x) , np.array(self.w)*([g[k-1]]*len(self.x)))
                else:
                    gama_pre=np.dot([g[k-1].subs(t,x) for x in self.x] , np.array(self.w)*[g[k-1].subs(t,x) for x in self.x])
                b=beta/gama
                c=gama/gama_pre
                g[k+1]=simplify((t-b)*g[k]-c*g[k-1])
            self.g=g
        #生成线性无关的x的k次方函数
        elif function_name == 'x^k': 
            value_num=m
            g={}
            g[0]=1
            g[1]=t
            for k in range(1,value_num):
                g[k+1]=t**(k+1)
            self.g=g
 
    #正交基函数正规方程组求解
    def ZhengJiao_solve_ck(self):
        for k in range(self.n+1):
            if k==0:
                c1=np.dot([self.g[k]]*len(self.x) , np.array(self.w)*self.y)
                c2=np.dot([self.g[k]]*len(self.x) , np.array(self.w)*([self.g[k]]*len(self.x)))
            else:
                c1=np.dot([self.g[k].subs(t,x) for x in self.x] , np.array(self.w)*self.y)
                c2=np.dot([self.g[k].subs(t,x) for x in self.x] , np.array(self.w)*[self.g[k].subs(t,x) for x in self.x])
            self.c[k]=c1/c2

    #最小二乘多项式正规方程组求解
    def ZXEC_solve_ck(self):
        #计算A和b
        G_pre=[]
        for key in self.g.keys():
            if key==0:
                pass
            else:
                G_pre.append(self.g[key])
        G1=[[g.subs(t,x) for g in G_pre] for x in self.x]
        one_temp=np.ones(len(self.x))
        one_temp=one_temp.reshape(len(self.x),1)
        G2=np.c_[one_temp,G1]
        G=np.array(G2)
        W=np.diag(self.w)
        A_temp=np.dot(G.T,W)
        A=np.dot(A_temp,G)
        b=np.dot(np.dot(G.T,W),self.y)
        #LDLT分解
        L,U=self.my_LU(A)
        D=np.diag(U)
        D=np.diag(D)
        #解y
        y={}
        y[0]=b[0]
        for i in range(1,self.n+1):
            Ly=[L[i][k]*y[k] for k in range(0,i)]
            sum_Ly=sum(Ly)
            y[i]=b[i]-sum_Ly
        #解x
        x={}
        x[self.n]=y[self.n]/U[self.n][self.n]
        for j in range(self.n-1,-1,-1):
            Ux=[U[j][k]*x[k] for k in range(j+1,self.n+1)]
            sum_Ux=sum(Ux)
            x[j]=(y[j]-sum_Ux)/U[j][j]
        #将x变为c
        for i in range(len(x)):
            self.c[i]=x[i]

    #LU分解
    def my_LU(self,A):
        A=np.asarray(A)
        U=np.zeros_like(A)
        L=np.zeros_like(A)
        U[0,:]=A[0,:]
        L[1:A.shape[0],0]=A[1:A.shape[0],0]/U[0][0]
        for i in range(1,A.shape[0]-1):
            U[i][i]=A[i][i]-sum([L[i][k]*U[k][i] for k in range(0,i)])
            for j in range(i+1,A.shape[1]):
                U[i][j]=A[i][j]-sum([L[i][k]*U[k][j] for k in range(0,i)])
                L[j][i]=(A[j][i]-sum([L[j][k]*U[k][i] for k in range(0,i)]))/U[i][i]
        U[A.shape[0]-1][A.shape[1]-1]=A[A.shape[0]-1][A.shape[1]-1]-\
                                      sum([L[A.shape[0]-1][k]*U[k][A.shape[1]-1] for k in range(0,A.shape[1]-1)])
        for i in range(L.shape[0]):
            L[i][i]=1
        return L,U   
    
    #计算逼近多项式
    def solve_pn(self):
        pn=[]
        for k in range(self.n+1):
            pn.append(simplify(self.c[k]*self.g[k]))
        result=0
        for p in pn:
            result=result+p
        return simplify(result)

    #计算误差
    def LeastSquareFit_error(self):
        y_pn=self.solve_pn()
        e=0
        for index,value in enumerate(self.x):
            temp=abs(y_pn.subs(t,value)-self.y[index])
            e=e+temp**2
        return pow(e,0.5)

######################测试程序####################
if __name__ == "__main__":
    #------------------使用正交多项式---------------------
    # X=[0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
    # Y=[5.1234, 5.3057, 5.5687, 5.9375, 6.4370, 7.0978, 7.9493, 9.0253, 10.3627]
    # test=LeastSquareFit(X,Y,4,JiFunction='正交多项式')
    # test.myset_w(f=1) #设置权函数
    # test.ZhengJiao_solve_ck() #求解系数c
    # print("多项式系数", test.c)
    # result=test.solve_pn() #求解拟合的多项式
    # print("拟合多项式", result)
    # my_error=test.LeastSquareFit_error() #计算误差
    # print("误差为：", my_error)
    #------------------使用x^k---------------------------
    X = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
    Y = [5.1234, 5.3057, 5.5687, 5.9375, 6.4370, 7.0978, 7.9493, 9.0253, 10.3627]
    test=LeastSquareFit(X,Y,4,JiFunction='x^k')
    test.myset_w(f=1) #设置权函数
    test.ZXEC_solve_ck() #求解多项式系数c
    print("多项式系数", test.c)
    result=test.solve_pn() #求解拟合的多项式
    print("拟合多项式", result)
    my_error=test.LeastSquareFit_error() #计算误差
    print("误差为：", my_error)

4. 程序使用说明

        首先给定输入数据X和Y，然后确定多项式阶数比如4，最后选择基函数，可选“正交基函数”或“x^k”。确定完参数后，具体使用步骤为：（1）类实例化；（2）调用myset_w()设置权函数；（3）调用ZXEC_solve_ck()（或ZhengJiao_solve_ck()）求解多项式系数；（4）调用solve_pn()生成拟合多项式；（5）调用LeastSquareFit_error()计算误差。