多项式曲线拟合

2023-08-04 08:51| 来源: 网络整理| 查看: 265

利用多项式函数拟合数据点，多项式函数形式如下： $y(x, W) = w_0 + w_1 x + \cdots + w_m x^m = \sum_{i=0}^{m} {w_i x^i}$ 令 $W = \left[ \begin{matrix} w_0 \\ w_1 \\ \vdots \\ w_m \end{matrix} \right]$ ， $X = \left[ \begin{matrix} 1 & x_1 & \cdots & x_1^m \\ 1 & x_2 & \cdots & x_2^m \\ \vdots & \vdots & \ddots & \vdots\\ 1 & x_n & \cdots & x_n^m \end{matrix} \right]$ 则多项式函数可化为线性代数形式： $y(x, W) = XW$ 为了评价拟合函数的优劣，需要建立损失函数，测量每个样本点目标值与预测值之间的误差，拟合的目标是让误差最小。计算误差时使用均方根误差 $E_{RMS} = \sqrt {\frac {E(W^*)}{N}}$

最小二乘法

为了取得误差函数的最小值，直接令函数的导数等于0，可以得到唯一解。此解称为解析解。不带正则项的解析解误差函数为： $E(W) = \frac 1 2 \sum_{n=1}^{N} {(y(x_n, W) - t_n)^2} = \frac 1 2 (XW - T)^{\top}(XW - T)$ 对其求导得： $\frac {\partial E(W)}{\partial W} = {X^{\top}XW - X^{\top}T}$ 令导数等于0，得： $W = (X^{\top}X)^{-1}X^{\top}T$ 可以通过矩阵运算直接目标函数：W = (X.T * X).I * X.T * T

随着阶数越高，误差越小，但是过拟合越来越严重：

观察计算出来的W，会看到，阶数越高，多项式的系数越大的离谱

为了解决过拟合问题，可以增加样本数据，也可以加入正则项（惩罚项）下图为9阶多项式在n = 15， 30， 60， 100时的拟合情况，可以看到曲线随着数据的增加拟合得越接近正弦曲线。

带正则项的解析解误差函数变成：这里的m+1为W中元素的个数，，主要是为了平衡高阶多项式项更多带来的影响，不加对下面W的求解也没有影响

$\lambda$ ，也就是正则系数，是一个超参数，需要设置一个较好的值才能得到较好的结果。这里取的是0.005

求导： $\frac {\partial \widetilde E (W)}{\partial W} = \frac 1 {m+1}({(X^{\top}X + \lambda E_{m+1})W - X^{\top}T})$ 求出W: $W = (X^{\top}X + \lambda E_{m+1})^{-1}X^{\top}T$ 同样可以直接由矩阵运算求得W = (X.T * X + lambda * np.eye(m+1)).I * X.T * T

由于加上了正则项，高阶多项式没有出现明显的过拟合

多项式系数的大小也较为正常：

最小二乘法看似简单，但是涉及到矩阵求逆运算，当数据规模较大时，矩阵求逆速度较慢，因此需要其他优化方法，如梯度下降法，共轭梯度法，牛顿法等。梯度下降法

梯度下降法的原理可以看这篇文章：深入浅出--梯度下降法及其实现

拟合的目标时让误差函数值最小，也就是找函数图像的最低点，沿着曲线梯度的反方向一直走，一定会走到一个局部最低点，也就是局部最优解，在二次型，即为全局最优解。

def gradient(W, X, T): '''在带正则项的最小二乘函数中计算W处的梯度''' return (1 / X.shape[1]) * (X.T * X * W - X.T * T + lambd * W) def gradient_descent(X, T: '目标值', learning_rate): '''进行梯度下降的迭代''' # 初始化W为全为0的列向量 W = np.mat(np.zeros(X.shape[1])).T grad = gradient(W, X, T) while np.all(np.absolute(grad) > 1e-9): W = W - learning_rate * grad grad = gradient(W, X, T) return W

拟合的曲线与带正则项的解析解几乎一致

但是阶数越高，迭代的次数越多，速度越慢，进行十几万次迭代的时间约为十几秒

阶数 0 1 2 3 4 5 6 迭代次数 28 352 8337 93643 133747 160849 181179 共轭梯度法

该方法解决的是二次函数极值问题 $\min \limits_{x \in R^n} \frac 1 2 x^{\top}Qx - b^{\top}x$ 其中Q为正定矩阵，本问题中 $Q = \frac 1 {m+1}(X^{\top}X + \lambda E_{m+1})$

共轭梯度法每次沿着一个方向找到该方向的极小值，后面在沿着其他方向求极小值，由于寻找的方向两两共轭，后面的查找不会影响前面已经查找过的方向的最小值，理论上n次求极小值就可以得到n为问题的极小值。

共轭梯度法理解起来比较难，不过编程实现按步骤来就可以了：

def conjugate_gradient(X, T): '''进行共轭梯度下降的迭代''' Q = (1 / X.shape[1]) * (X.T * X + lambd * np.mat(np.eye(X.shape[1]))) W = np.mat(np.zeros(X.shape[1])).T r = -gradient(W, X, T) p = r for i in range(1, X.shape[1]): a = float((r.T * r) / (p.T * Q * p)) r_prev = r W = W + a * p r = r - a * Q * p p = r + float((r.T * r) / (r_prev.T * r_prev)) * p return W

共轭梯度法拟合的情况也很好，几乎与带正则项的解析解相同，并且速度要比梯度下降法快很多。附：矩阵求导用到的公式有

$\frac {\partial {\mathbf {a}}^{\top }{\mathbf {x}}}{\partial {\mathbf {x}}} = \frac {\partial {\mathbf {x}}^{\top }{\mathbf {a}}}{\partial {\mathbf {x}}} = {\mathbf {a}}$ $\frac {\partial {\mathbf {x}}^{\top}{\mathbf {x}}}{\partial {\mathbf {x}}} = 2\mathbf{x}$ $\frac {\partial {\mathbf {x^{\top}Ax}}}{\partial {\mathbf {x}}} = 2\mathbf{Ax}$ (A为对称矩阵)

参考：深入浅出--梯度下降法及其实现 -简书梯度下降法和共轭梯度法有何异同？-知乎

【本文地址】

公司简介

联系我们