利用多项式函数拟合数据点，多项式函数形式如下：
令
，
则多项式函数可化为线性代数形式： 为了评价拟合函数的优劣，需要建立损失函数，测量每个样本点目标值与预测值之间的误差，拟合的目标是让误差最小。计算误差时使用均方根误差
最小二乘法
为了取得误差函数的最小值，直接令函数的导数等于0，可以得到唯一解。此解称为解析解。
不带正则项的解析解
误差函数为：
对其求导得：
令导数等于0，得：
可以通过矩阵运算直接目标函数：
W = (X.T * X).I * X.T * T
随着阶数越高，误差越小，但是过拟合越来越严重：
下图为9阶多项式在n = 15， 30， 60， 100时的拟合情况，可以看到曲线随着数据的增加拟合得越接近正弦曲线。
带正则项的解析解
误差函数变成：
这里的m+1为W中元素的个数，，主要是为了平衡高阶多项式项更多带来的影响，不加对下面W的求解也没有影响
，也就是正则系数，是一个超参数，需要设置一个较好的值才能得到较好的结果。这里取的是0.005
求导：
求出W:
同样可以直接由矩阵运算求得
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 |
共轭梯度法
该方法解决的是二次函数极值问题
其中Q为正定矩阵，本问题中
共轭梯度法每次沿着一个方向找到该方向的极小值，后面在沿着其他方向求极小值，由于寻找的方向两两共轭，后面的查找不会影响前面已经查找过的方向的最小值，理论上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
共轭梯度法拟合的情况也很好，几乎与带正则项的解析解相同，并且速度要比梯度下降法快很多。
附：矩阵求导用到的公式有
-
(A为对称矩阵)
