2. 多项式回归详解：从原理推导到实战预测

1 前言

在上一章节，我们学习了如何使用线性回归对数据进行预测，但发现一元线性回归并不能很好的反应某些数据。

因此本章将要继续学习更高级的内容，让我们拥有回归预测不那么线性的数据，也就是多项式回归。

所谓多项式，指的是由变量和系数常量通过有限次加减乘除以及自然幂次的乘方运算得到的表达式，是整式的一种。

当未知数只有一个时，就被成为一元多项式，也就是上一章节所介绍的内容，比如$x^2-x+4$。

而未知数不止一个的多项式就称为多元多项式，例如$x^3-3xy^2z-yz+1$就是一个三元多项式。

2 认识多项式

为了更加直观的演示，首先这里提供一组数据，然后将其画出来：

from matplotlib import pyplot as plt

x = [4, 7, 13, 26, 31, 41, 59, 62, 71, 78]
y = [21, 35, 49, 54, 43, 30, 31, 48, 66, 74]

plt.scatter(x, y)
plt.show()

效果如下：

很明显可以看出来，这组数据具有明显的波动性，如果用前面学过的内容、用直线去拟合这组数据肯定效果不好。

因此下面我们来尝试使用多项式来拟合数据，首先尝试用二次多项式。

一个标准的一元高阶多项式函数如下：

$y(x, \mathbf{w}) = w_0 + w_1 x + w_2 x^2 + \cdots + w_m x^m = \sum_{j=0}^{m} w_j x^j$

其中m代表多项式的阶数，$x^j$表示x的j次幂，w代表x的系数。

因此当我们尝试用上面的一元高阶多项式去拟合数据时，就需要确定多项式系数w与阶数m，这是多项式的两个基本要素。

由于阶数我们可以实现尝试确定、比如这里设定阶数为2，那么就只需要求解系数w即可：

$y(x, \mathbf{w}) = w_0 + w_1 x + w_1 x^2$

此时，我们又回到了前面的内容，依旧是自定义损失函数：

def func(p, x):
    # 根据公式，定义 2 次多项式函数
    w0, w1, w2 = p
    x = np.asarray(x)
    f = w0 + w1 * x + w2 * x**2
    return f


def err_func(p, x, y):
    # 残差函数（观测值与拟合值之间的差距）
    f=func(p,x)
    return f-y

然后调用求解最佳拟合参数的函数即可：

from scipy.optimize import least_squares
#...
res=least_squares(err_func,[0,0,0],args=(x,y))

这里使用了scipy库内部的least_squares函数，通过最小二乘法求解最佳拟合参数，而无需自己去手动实现。

如果不存在该库可能需要自己安装一下：

uv pip install scipy

除非你需要更多自定义行为，否则只需要上面三个参数即可：损失函数、初始系数值、训练数据。

这里为了简单，直接设置了三个0作为初始值，该值并不影响结果，但它会影响最终最终多项式的次数，所以这里涉及三个，代表求解2次多项式（0、1、2次共三个系数）。

下面是一个完整的实例程序：

import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import least_squares
x = [4, 7, 13, 26, 31, 41, 59, 62, 71, 78]
y = [21, 35, 49, 54, 43, 30, 31, 48, 66, 74]

def func(p, x):
    w0, w1, w2 = p
    x = np.asarray(x)
    f = w0 + w1 * x + w2 * x**2
    return f
def err_func(p,x,y):
    f=func(p,x)
    return f-y

res=least_squares(err_func,[0,0,0],args=(x,y))

tmp_x=np.linspace(0,100,1000)
tmp_y=func(res.x,tmp_x)
plt.plot(tmp_x, tmp_y, 'r-', label='Fitted Curve')

plt.scatter(x, y)
plt.show()

得到结果的res上的x字段便是求解得到了拟合系数，然后这里生成了一组0-100共1000个数据作为横坐标、再通过该系数、横坐标计算得到纵坐标，将其画出来。

效果如下：

可以看到，此时拟合曲线已经出现了弯曲，但和数据点的拟合效果却并不好，这说明一元二次项依旧无法正确的反应数据点的变化趋势。

3 n次多项式

此时的方式自然就是增加拟合方程的次数了，二次不行，那就三次、四次，直到能正确拟合数据。

为了简化这个步骤，这里封装一个n次多项式的代码：

import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import least_squares
x = [4, 7, 13, 26, 31, 41, 59, 62, 71, 78]
y = [21, 35, 49, 54, 43, 30, 31, 48, 66, 74]

def fit_func(p, x):
    """根据公式，定义 n 次多项式函数"""
    f = np.poly1d(p)
    return f(x)


def err_func(p, x, y):
    """残差函数（观测值与拟合值之间的差距）"""
    ret = fit_func(p, x) - y
    return ret


def n_poly(n):
    """n 次多项式拟合"""
    p_init = np.random.randn(n)  # 生成 n 个随机数
    parameters = least_squares(err_func, p_init, args=(np.array(x), np.array(y)))
    return parameters.x


# 绘制拟合图像时需要的临时点
x_temp = np.linspace(0, 100, 1000)

# 绘制子图
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

axes[0, 0].plot(x_temp, fit_func(n_poly(4), x_temp), "r")
axes[0, 0].scatter(x, y)
axes[0, 0].set_title("m = 3")

axes[0, 1].plot(x_temp, fit_func(n_poly(5), x_temp), "r")
axes[0, 1].scatter(x, y)
axes[0, 1].set_title("m = 4")

axes[0, 2].plot(x_temp, fit_func(n_poly(6), x_temp), "r")
axes[0, 2].scatter(x, y)
axes[0, 2].set_title("m = 5")

axes[1, 0].plot(x_temp, fit_func(n_poly(7), x_temp), "r")
axes[1, 0].scatter(x, y)
axes[1, 0].set_title("m = 6")

axes[1, 1].plot(x_temp, fit_func(n_poly(8), x_temp), "r")
axes[1, 1].scatter(x, y)
axes[1, 1].set_title("m = 7")

axes[1, 2].plot(x_temp, fit_func(n_poly(9), x_temp), "r")
axes[1, 2].scatter(x, y)
axes[1, 2].set_title("m = 8")

plt.scatter(x, y)
plt.show()

上面的代码中，直接调用np库的poly1d函数，通过传入的系数、生成一个多项式函数对象，避免我们手动去定义。

此外，这里调用的least_squares函数初始化系统也直接更改为了n个随机数作为初始化值、或者你更改为n个0的一维数组也可以。

此时只需要调用n_poly函数、传入数字n、就能自动拟合n-1次多项式的系数，并将其作为子图画了出来：

可以看到，但m=4，也就是4次多项式时，效果就已经明显优于m=3时了，但随着m次数的增长，其拟合的效果却并没有越来越好。

当m大于5之后，整个曲线便出现了非常明显的震荡，其预测的曲线明显偏离原数据，出现快速、大范围的增加、或减少。

而这就是线性回归中的过拟合现象。

4 scikit-learn多项式拟合

事实上，对于多项式回归拟合来说，除了像上面那样自定定义多项式来拟合外，也可以使用scikit-learn这个库的sklearn.preprocessing.PolynomialFeatures类直接进行多项式回归拟合结果。

这个类的主要作用便是产生多项式特征矩阵。

比如一个二次多项式，我们一般的写法是：$y(x,w)=w_0+w_1x+w_2x^2$。

但我们可以将其转化一下，令：$x=x_1,x^2=x_2$。

那么此时原方程就变成了：$y(x,w)=w_0+w_1x_1+w_2x_2$。

看到这里是不是突然就感觉熟悉了，这不就是上一章节用到的多元线性回归嘛！甚至上一章的实例中，我们用的还是三元进行预测数据。

这类似于高中学过的换元法，将一元高次方程组转换为了多元一次多项式。

转换为矩阵的形式如下，如果训练数据X的值为：

$\begin{bmatrix} 1\\ -1\\ 2 \end{bmatrix}$