目标
掌握最小二乘法求解（无惩罚项的损失函数）、掌握加惩罚项（2范数）的损失函数优化、梯度下降法、共轭梯度法、理解过拟合、克服过拟合的方法(如加惩罚项、增加样本)。
要求
生成数据，加入噪声；
用高阶多项式函数拟合曲线；
用解析解求解两种loss的最优解（无正则项和有正则项）；
优化方法求解最优解（梯度下降，共轭梯度）；
用你得到的实验数据，解释过拟合。
用不同数据量，不同超参数，不同的多项式阶数，比较实验效果。
语言不限，可以用matlab，python。求解解析解时可以利用现成的矩阵求逆。梯度下降，共轭梯度要求自己求梯度，迭代优化自己写。不许用现成的平台，例如pytorch，tensorflow的自动微分工具。
实验环境
Windows 10
PyCharm 2020.3
Python 3.7.0
实验步骤
1. 生成实验数据
运用随机数生成器，生成一组满足在 [0, 1] 区间上均匀分布的（x, y）数据，其中
y = sin ( 2 π x ) + e y = \sin(2\pi x) + e y=sin(2πx)+e
而e即为为满足 (0, 0.01) 正态分布的一个噪声。
def generate_sin2pi_sample(number):
x = np.linspace(0, 1, number).T
y = np.sin(2 * np.pi * x) + np.random.normal(0, 0.1, number).T
return x, y
接下来将利用这组数据进行训练。
接下来为了便于矩阵计算，将会用到x的Vandermende矩阵X。
Vandermende矩阵的计算代码如下：
def vandermonde(x, order):
van = np.zeros(shape=(order + 1, x.size))
van[0] = 1
for i in range(order):
van[i + 1] = van[i] * x.T
return van
2. 利用解析解拟合正弦曲线
先通过最简单的思路解决问题。
当我们想要拟合一条曲线时，通常我们得到的只有几个离散的数据样本 (x, y)，并且与实际的曲线之间还有一些误差。我们对于目标函数的了解仅此而已。
在这种已知条件受限的情况下，我们首先需要假设一个大概的模型，比如本次实验我们采用了多项式函数 Y(x, w) 作为拟合的模型。在此之上，还需要对这种模型的一些待定参数进行调整，使得模型能够在给定样本以外的输入上，也能做到与目标函数尽可能贴近。这种贴近的标准可以用一个与待定参数以及输入样本相关的函数来标定，称为loss函数。
这里我们采用了以下公式作为loss函数：
E ( w ) = 1 2 ∑ i = 1 n ( Y ( w , x i ) − y ( x i ) ) 2 E(w)=\frac{1}{2} \sum_{i=1}^n{(Y(w,x_i)-y(x_i))^2} E(w)=21i=1∑n(Y(w,xi)−y(xi))2
设多项式的次数为n，式中的 w 为 n+1 维的系数列向量，Y(w, x) 即为作为假设模型的n阶多项式函数。利用由 x 引申出的Vandermende矩阵，可以更简洁地表示成矩阵形式：
E ( w ) = 1 2 ( X T w − y ) T ( X T w − y ) E(w) = \frac{1}{2}(X^Tw-y)^T(X^Tw-y) E(w)=21(XTw−y)T(XTw−y)
在这种条件下，我们认为当 w 可以使得 E(w) 达到最小时，预测模型 Y(w, x) 就在给定的样本上尽可能地逼近了目标函数，由此这样的 w 将作为预测模型 Y(x, w) 的参数列向量。故只要 E(w) 对 w 的导数满足：
E ′ ( w ) = X X T w − X y = 0 E'(w) = XX^Tw-Xy = 0 E′(w)=XXTw−Xy=0
⇒ w = ( X X T ) − 1 X y \Rightarrow w = (XX^T)^{-1}Xy ⇒w=(XXT)−1Xy
由此，我们得到了直接通过样本得到结果的算法：
def analysis_solution(x, y, order, rambda=0):
van_x = vandermonde(x, order)
w = np.dot(np.linalg.inv(np.dot(van_x, van_x.T) + rambda * np.eye(order + 1)), np.dot(van_x, y))
return w
我针对不同的多项式次数（order），以及给定样本训练集大小（sample number），利用同分布的10个样本作为验证集，得到多项式拟合函数的 E(w) （error），以此作为拟合偏差的评估标准。
这里同一行的多项式拟合曲线次数一致，同一列的训练集大小一致。从总体上看，多项式次数和训练集个数都不是越高或者越低越好。
这里的图像中，order=7，sample_number=100时拟合效果最好，order=7，sample_number=50则次之。
横向比较：order=3时，sample_number越大，误差反而越大。order=12时，随着sample_number增大，error先增大，后减小。
纵向比较：order=3时，拟合效果都是最好的。order=12时，图像都出现了各种程度的震荡，拟合程度也不佳。
3. 利用有正则项的解析解拟合正弦曲线
在进行拟合的时候，随着拟合曲线训练程度加深，将会逐步贴近目标函数曲线。但是如果只用目前形式 E(w) ，函数在最优化的时候，只会保证尽可能讲过所有的给定训练集上的点，但是在其他点上的拟合效果就未必有那么好了。
从一定程度上，曲线的复杂程度与系数列向量 w 是相关的。w 的“绝对值”越小，曲线越区域平缓，拟合程度一般都会更好。所以，我们将 w 的2-范数也添加到 E(w) 中，使得训练过程中，算法也会将 w 的2-范数纳入最优化的考察范围中。由此改写 E(w) 如下：
E ( w ) = 1 2 ( X T w − y ) T ( X T w − y ) + λ 2 w T w E(w) =\frac{1}{2}(X^Tw-y)^T(X^Tw-y) + \frac{\lambda}{2}w^Tw E(w)=21(XTw−y)T(XTw−y)+2λwTw
代码同无正则项的解析解之算法，只是传入的参数 rambda 将会是一个正数。
def analysis_solution(x, y, order, rambda=0):
van_x = vandermonde(x, order)
w = np.dot(np.linalg.inv(np.dot(van_x, van_x.T) + rambda * np.eye(order + 1)), np.dot(van_x, y))
return w
取 rambda = 0.001，图像如下：
取rambda=0.01，图像如下：
取rambda=0.1，图像如下：
可以看到，随着rambda增加，图像整体曲线将趋于平缓，尤其是有效解决了高次数（order=12）时，曲线震荡的问题。但是，rambda也不是越大越好，如果rambda过大，将会使得曲线整体与训练集的点的偏移增加。拟合效果反而变差了。
4. 利用梯度下降法拟合正弦曲线
梯度下降法是一种迭代求解的技术，针对想要最优化为最小值的多元函数 E(w) ，梯度下降法每次迭代都求解当前 w 的梯度，并让 w 向逆梯度方向发生一定变化，由此不断迭代到想要的解的邻域内，在一定误差范围内即可输出当前 w 作为 min{E(w)} 的解。
这种方法实际上只能找到局部最优解，但是由于 E(w) 是一个凸函数，因此局部最优即为全局最优，适合使用梯度下降法。
具体算法代码如下，其中输入参数rambda=0：
def gradient_descent(x, y, order, acceptable_grad, learning_rate, rambda):
van_x = vandermonde(x, order)
w = np.zeros(order + 1).T
grad = np.dot(van_x, Y(w, van_x) - y) + (rambda * w)
cnt = 0
while np.sum(np.abs(grad)) >= acceptable_grad:
w = w - learning_rate * grad
grad = np.dot(van_x, Y(w, van_x) - y) + (rambda * w)
if cnt % 100000 == 0:
print("cnt[", cnt, "] grad:")
print(grad)
print("w:")
print(w)
cnt += 1
return w
其中的cnt局部变量和标准输出语句是为了确认迭代速度，结束条件是当前 w 所有在梯度的 1-范数 小于一定值acceptable_grad。
以下为输出结果：
当acceptable_grad = 0.01时
当acceptable_grad = 0.1时：
另外，梯度下降法比起其他算法，需要的时间都是没有保证的。
5. 利用有正则项的梯度下降法拟合正弦曲线
同样的道理，正则项也可以放到梯度下降的损失函数的。
以下为acceptable_grad=0.001，rambda=0.01的图像。
6. 利用共轭梯度法拟合正弦曲线
共轭梯度法每次选取一个与之前变化方向共轭的新方向，让 w 向特定的方向进行变化，理论上可以在n步以内达到最优解。（n为目标函数所在空间的维数）
代码如下，其中rambda=0
def conjugate_gradient(x, y, order, acceptable_grad, rambda):
van_x = vandermonde(x, order)
w = np.zeros(order + 1).T
A = np.dot(van_x, van_x.T) + rambda * np.eye(order + 1)
b = np.dot(van_x, y)
p = b - np.dot(A, w)
beta = 0
cnt = 0
while True:
grad = np.sum(np.abs(np.dot(A, w) - b))
if grad < acceptable_grad:
print("cnt: ", cnt)
print(w)
return w
if cnt % 100000 == 0:
print("cnt: ", cnt, " grad: ", grad)
r = b - np.dot(A, w)
beta = beta * np.dot(-1 * r.T, -1 * r)
p = r + beta * p
beta = 1 / np.dot(-1 * r.T, -1 * r)
w = w + (-1) * (np.dot(-1 * r.T, p) / np.dot(np.dot(p.T, A), p)) * p
cnt += 1
print("cnt: ", cnt)
print(w)
return w
这里在n步之后输出 w，或者在n不之内已经达到了预想的误差范围内的话，直接输出结果。
7. 利用有正则项的共轭梯度法拟合正弦曲线
同理使用正则项，应用上一节的代码。其中rambda>0。
rambda=1e-3：
rambda=1e-10：
rambda=1e-100：
这里可以看出正则项的出现，反而让拟合结果恶化了。
理解过拟合
可以看到在上述实验结果当中，当rambda很小，样本数量较少，以及多项式次数较高时，都容易产生类似于图中的现象。也许函数还能比较好的拟合给定的训练集，但是明显偏离了目标函数。
要解决这种问题，一般有3中思路：
一个是降低多项式次数：
一个是增加训练集的大小：
一个是添加正则项：
其中，正则项的参数lambda的选取是有讲究的。lambda不能太小，否则正则项的作用微乎其微；正则项也不能太大，否则拟合曲线将会变成不会经过训练集样本的“平滑”曲线：
思考
- 我们在进行曲线拟合的时候，运用的评估标准实际上并不能够直接反应拟合曲线对于其他点的拟合情况，而只能让曲线尽可能接近给定的训练集。同样，在使用正则项的时候，正则项也并不能直接反应曲线的曲折程度，也只是一个简洁的量化标准。要想得到一组好的参数，还是需要不断和其他参数以及模型相互比较，才能得到比较实在的结果。
- 运用一个最优化算法时，时间复杂度也是一个非常重要的评估标准。比如相对于解析解法和共轭梯度法，梯度下降法的时间复杂度是几乎不可控的，经常会一跑一晚上。在参数优化选择的算法还不明确的情况下，最好的解决办法就是将整个迭代过程的数值变化输出出来，实时监控；当变化非常小，但是离目标值还非常远的时候，就可以停止当前迭代，及时止损，使用更加宽松的参数条件得到结果。
- 在进行训练集以及验证集的选择时，在各种向度上都实现等比例划分是困难的。因此这次我使用了两个满足相同正态分布的样本集合，分别作为训练集和验证集。这是严谨的吗？从分布上看，二者确实做到了等价。那么在现实样本当中，为什么还要对得到的样本进行划分呢？前一天选取训练集，后一天选取验证集不好吗？其实在实际场景下，不会有那么大量的数据。在数据量不大的时候，又要拆分成小的集合时，单个集合本身已经很难体现原有的分布了。因此尽管可能有违我们对于“随机采样”的直觉，但是我们需要对数据进行一些预处理，使得划分处的集合能够体现出相同的分布特性。
- 接受条件的选择应该基于损失函数的导数，而非损失函数本身。即使用梯度下降法时，结束条件应该是 E’(w)<rate，而不是 E(w) < rate，尽管E(w)可能为0。一方面，利用导数更具有普适性；另一方面，E(w)的最小值在一定步长的情况下，可能非常难达到其邻域。
- 更多、更能覆盖样本空间整体的训练集，可以让拟合曲线更贴合现实。我们这次实验选取的样本点均匀分布在考察区间内，实际上是简化了模型训练难度。如果样本都集中在考察区间一个点的附近，那么训练出来的模型在其他的点上拟合效果就很难控制了。
- 更复杂的模型表现会更加灵活，但也需要更多训练样本来支持其模型的训练；如果样本数量不够，也可以加入正则项，降低其模型的复杂度。
- 参数的选取都是一个折衷的选择，选大选小都有各自的优劣。模型的选取也是。
