M次多项式拟合

　

多项式拟合

M次多项式拟合问题实际就是一个最小二乘法的问题，作者在《统计学习方法》中并没有给出具体的推导公式，下文给出具体的推导公式。

设M次多项式为

fM(xI,w)=w0+w1xi+w2xi2+⋯+wMxiM=∑j=0Mwjxj=XiTw" role="presentation" style="text-align: center; position: relative;">fM(xI,w)=w0+w1xi+w2x2i+⋯+wMxMi=∑j=0Mwjxj=XTiw$$f_M(x_I,\mathbf{w})=w_0+w_1{x}_i+w_2{x}_i^2+\cdots +w_M{x}_i^M=\sum _{j=0}^M{w}_j{x}^j=X_i^T\mathbf{w}$$

其中w=(w0,w1,⋯,wM)" role="presentation" style="position: relative;">w=(w0,w1,⋯,wM)$w=(w_0,w_1,\cdots ,w_M)$,Xi为矩阵

X=[1x1x12⋯x1M⋮⋮1xnxn2⋯xnM]" role="presentation" style="text-align: center; position: relative;">X=⎡⎣⎢⎢⎢1⋮1x1xnx21x2n⋯⋯xM1⋮xMn⎤⎦⎥⎥⎥$$X=\left[\begin{array}{ccccc}1& x_1& x_1^2& \cdots & x_1^M\\ ⋮& & & & ⋮\\ 1& x_n& x_n^2& \cdots & x_n^M\end{array}\right]$$

的第 i 行。于是我们可以得到损失函数

L(w)=12∑i=1N(f(xi,w)−yi)2=12∑i=1N(∑j=0Mwjxij−yi)2=12||Xw−y||22" role="presentation" style="text-align: center; position: relative;">L(w)=12∑i=1N(f(xi,w)−yi)2=12∑i=1N(∑j=0Mwjxji−yi)2=12||Xw−y||22$$\begin{array}{rl}L(\mathbf{w})& =\frac{1}{2}\sum _{i=1}^N(f(x_i,\mathbf{w})-y_i{)}^2\\ & =\frac{1}{2}\sum _{i=1}^N(\sum _{j=0}^M{w}_j{x}_i^j-y_i{)}^2\\ & =\frac{1}{2}||X\mathbf{w}-\mathbf{y}|{|}_2^2\end{array}$$

这里的损失函数前面加了12" role="presentation" style="position: relative;">12$\frac{1}{2}$，是为了方便计算。此外，这里转换为向量更加方便求解，书中是直接公式推导，没有使用向量方法，我们等会再看那种方法。

为了使损失函数最小化，只需要对w" role="presentation" style="position: relative;">w$\mathbf{w}$进行求导，然后令导数为0就可以获得结果。

∂L(w)∂w=XT(Xw−y)=XTXw−XTy" role="presentation" style="text-align: center; position: relative;">∂L(w)∂w=XT(Xw−y)=XTXw−XTy$$\begin{array}{rl}\frac{\mathrm{\partial }L(\mathbf{w})}{\mathrm{\partial }\mathbf{w}}& =X^T(X\mathbf{w}-\mathbf{y})\\ & =X^{T}X\mathbf{w}-X^T\mathbf{y}\end{array}$$

从而可以求出来w=(XTX)−1XTy" role="presentation" style="position: relative;">w=(XTX)−1XTy$\mathbf{w}=(X^{T}X{)}^{-1}{X}^T\mathbf{y}$

" role="presentation" style="position: relative;">$$

接下来不使用向量进行求导：

L(w)=12∑i=1N(∑j=0Mwjxij−yi)2∂L(w)∂wk=12∑i=1N[2(∑j=0Mwjxij−yi)xik]=∑i=1N∑j=0Mwjxij+k−∑i=0Nyixik=0" role="presentation" style="position: relative;">L(w)=12∑i=1N(∑j=0Mwjxji−yi)2∂L(w)∂wk=12∑i=1N[2(∑j=0Mwjxji−yi)xki]=∑i=1N∑j=0Mwjxj+ki−∑i=0Nyixki=0$$\begin{gathered} L(\mathbf{w})=\frac{1}{2}\sum _{i=1}^N(\sum _{j=0}^M{w}_j{x}_i^j-y_i{)}^2 \\ \begin{array}{rl}\frac{\mathrm{\partial }L(\mathbf{w})}{\mathrm{\partial }{w}_k}& =\frac{1}{2}\sum _{i=1}^N[2(\sum _{j=0}^M{w}_j{x}_i^j-y_i)x_i^k]\\ & =\sum _{i=1}^N\sum _{j=0}^M{w}_j{x}_i^{j+k}-\sum _{i=0}^N{y}_i{x}_i^k=0\end{array} \end{gathered}$$

所以要拟合多项式系数w0,w1,⋯,wM" role="presentation" style="position: relative;">w0,w1,⋯,wM$w_0,w_1,\cdots ,w_M$,需要求解下方的方程组。为了方便，省略∑" role="presentation" style="position: relative;">∑$\sum$的上下标记

[N∑xi∑xi2⋯∑xiM∑xi∑xi2∑xi3⋯∑xIM+1⋮⋱⋮∑xiM∑xiM+1∑xiM+2⋯∑xi2M]{w0w1⋮wM}=[∑yi∑xiyi∑xi2yi⋮∑xiMyi]" role="presentation" style="text-align: center; position: relative;">⎡⎣⎢⎢⎢⎢⎢N∑xi⋮∑xMi∑xi∑x2i∑xM+1i∑x2i∑x3i∑xM+2i⋯⋯⋱⋯∑xMi∑xM+1I⋮∑x2Mi⎤⎦⎥⎥⎥⎥⎥⎧⎩⎨⎪⎪⎪⎪w0w1⋮wM⎫⎭⎬⎪⎪⎪⎪=⎡⎣⎢⎢⎢⎢⎢⎢∑yi∑xiyi∑x2iyi⋮∑xMiyi⎤⎦⎥⎥⎥⎥⎥⎥$$\left[\begin{array}{ccccc}N& \sum x_i& \sum x_i^2& \cdots & \sum x_i^M\\ \sum x_i& \sum x_i^2& \sum x_i^3& \cdots & \sum x_I^{M+1}\\ ⋮& & & \ddots & ⋮\\ \sum x_i^M& \sum x_i^{M+1}& \sum x_i^{M+2}& \cdots & \sum x_i^{2M}\end{array}\right]\left\{\begin{array}{c}{w}_0\\ w_1\\ ⋮\\ w_M\end{array}\right\}=\left[\begin{array}{c}\sum y_i\\ \sum x_i{y}_i\\ \sum x_i^2{y}_i\\ ⋮\\ \sum x_i^M{y}_i\end{array}\right]$$

之后就可以计算出各个w的值了。

下文中是一个Python多项式拟合的代码，参考别人的博客：HTTPs://blog.csdn.net/xiaolewennofollow/article/details/46757657。对于上面给出的线性方程组的公式，我们只需要求出左右矩阵和向量就可以获得w向量的结果。

import matplotlib.pyplot as plt
import math
import numpy
import random

fig=plt.figure()
ax=fig.add_subplot(111)

#生成数据点
x=numpy.arange(-1,1,0.1)
y = [((a*a-1)*(a*a-1)*(a*a-1)+0.5)*numpy.sin(a*2) for a in x]
plt.plot(x,y)
i=0
x_offset=[]
y_offset=[]
#生成的曲线上的点进行偏移，相当于加上噪声
for xx in x:
    yy=y[i]
    d=float(random.randint(60,140))/100
    i+=1
    x_offset.APPend(xx*d)
    y_offset.append(yy*d)
ax.plot(x_offset,y_offset,color='m',linestyle='',marker='.')
#这个函数是用来求系数的，order是多项式的幂次数,x_offset，y_offset是生成的数据
def get_w(order,x_offset,y_offset):
    #存储从0到m次的幂方和
    saveMat=[]
    for j in range(0,2*order+1):
        sum=0
        for i in range(0,len(x_offset)):
            sum+=(x_offset[i]**j)
        saveMat.append(sum)

    #求左边的矩阵
    matLeft=[]
    for row in range(0,order+1):
        rowvector=saveMat[row:row+order+1]
        matLeft.append(rowvector)

    matLeft=numpy.array(matLeft)

    #求右边的向量
    matRight=[]
    for i in range(0,order+1):
        y=0.0
        for k in range(0,len(x_offset)):
            y+=y_offset[k]*(x_offset[k]**i)
        matRight.append(y)

    matRight=numpy.array(matRight)

    W=numpy.linalg.solve(matLeft,matRight)
    return W

'''
order=3
W=get_w(order,x_offset,y_offset)
    #进行曲线的拟合
xxa= numpy.arange(-1,1,0.1)
yya=[]
for i in range(0,len(xxa)):
    yy=0.0
    for j in range(0,order+1):
        dy=(xxa[i]**j)
        dy*=W[j]
        yy+=dy
    yya.append(yy)
ax.plot(xxa,yya,color='g',linestyle='-',marker='')
ax.legend()
plt.show()
'''

看下图order为3和9时候的曲线拟合情况(绿色的线是进行拟合的曲线）

我们可以看出随着M（多项式次数）的增大，模型越来越复杂，拟合的程度越来越好，训练误差也越来越小，训练误差不断逼近0。但是测试误差确并不如此，它会随着模型的复杂度的增大先减小再增大。下图就是训练误差和测试误差与模型复杂度的关系图。

相关阅读