PyTorch Python 曲线拟合和优化问题白话文
刚开始深度学习和SLAM，其中涉及到了很多的优化问题，查看一些文章，好多都是对优化部分一句带过。网上的一些讲机器学习和梯度下降的博文，大多用线性拟合讲解梯度下降。今天，用代码实现了一个较为复杂函数的优化问题，份别用了不同的代码方式，包括使用PyTorch自动求导，PyTorch搭建神经网络和Scipy.optimize模块。另外，这里不涉及具体的理论，只提供了代码示例和必要的注释。具体代码请戳：PyTorch SGD,PyTorch NN ,Scipy Optimization (有些明明用的行内公式写的，预览的时候也正常，这里看着全部变成行间公式了)
问题描述
该问题描述来自《视觉SLAM14讲》。无偿安利一下高翔老师这本，非常适合SLAM初学者入门。以下这个优化问题在《视觉SLAM14讲》使用了C++代码，配合Google的Ceres库。
假设有一条满足以下方程的曲线：
其中，为曲线的参数,
为噪声。很明显，这是一个非线性模型。设有
个关于
的观测数据，目的是根据这些数据拟合出曲线的参数。这个问题可以转换为以下的优化问题：
注意！！！这里待估计的参数是，而不是
!!!
开始撸代码吧：
首先，我们产生一些观测数据，这些观测数据在三个不同方式的代码中都是一样的：
np.random.seed(1000)
torch.random.manual_seed(1000)
a = 2.0
b = 3.0
c = 1.0
N = 300
X_np = np.random.rand(N) * 1
noise = np.random.normal(loc=0.0, scale=1.0, size=N)*20
Y_np = np.exp(a*X_np**2 + b*X_np + c)
Y_observed = Y_np + noise
给个图吧，我们生成的数据就长这个样子：
1. 梯度下降实现
这里不在详细说明什么是梯度下降了，我认为核心公式就是：
即，不断更新参数,怎么更新呢，用损失函数相对于
的偏导，这个在代码中有体现的。这个代码中，我们使用了PyTorch的自动求导的功能，这里只贴出了部分代码，完整代码戳这里.
# 首先转换以下numpy的数据类型，PyTorch的求导好像支持float32，在其它地方改也可以
X_np = X_np.astype('float32')
noise = noise.astype(np.float32)
Y_np = Y_np.astype('float32')
Y_observed = Y_observed.astype('float32')
# 将Numpy 转为 torch.Tensor
X = torch.from_numpy(X_np).float().view(-1,1)
Y = torch.from_numpy(Y_observed).float().view(-1,1)
定义曲线方程：
def equation(a,b,c,X):
return torch.exp(a*X**2 + b*X + c)
开始训练：(注意代码中的注释)
# y_pred = equation(a,b,c,X)
# 这里的学习率也是多次实验之后实验出来的，可以在
# 可以改为一个大点的值，然后在我下方标注“打断点”的
# 地方打上断点进行调试，可以发现a.grad的值非常大，
# 这个可以根据y = np.exp(a*x*x + b*x +c)的图像理解
lr = 0.00000001
losses = []
for e in range(1000000):
pred = equation(a, b, c, X)
loss = torch.mean((Y-pred)**2)
if a.grad is not None and b.grad is not None and c.grad is not None:
# PyTorch计算的梯度是会累加的，所以每次计算前要清零，
# 但是，在没有调用backward()之前，a.grad是None，直接调用会报错，所以加了条件判断
a.grad.zero_()
b.grad.zero_()
c.grad.zero_()
loss.backward() # 计算梯度
a.data = a.data - lr*a.grad # 打断点
b.data = b.data - lr*b.grad
c.data = c.data - lr*c.grad
losses.append(loss.item())
print("Epoch: {} Loss: {}".format(e, loss.item()))
注释中说了，可以根据曲线的图像理解，我就把曲线的图像贴出来：
可见，改曲线是非常非常陡峭的，
关于学习率的设定，在注释中已经说了，可以按照这种方法打个断电，我相信就一幕了然，总是听说学习率很重要。一些其它的博客，使用线性规划，好像随便设置一个学习率，比如0.01, 0.1之类的，就可以的到很好的结果。用这个例子你试试，跑不死你......
下面给出迭代500个epoch 和 1000000 个epoch的截图：
500：
1000000:
可见，迭代1000000的时候，得到的参数
总结：
- 梯度下降中，学习率的设置确实很重要
- 梯度下降不一定是最好的方法，对初值的设置比较敏感
2. PyTorch 神经网络实现
这次，我搭建了一个有两个隐藏层的全连接网络，每个全连接层有10个神经元，具体看代码吧:
# 网络架构
class LS(nn.Module):
def __init__(self):
super(LS, self).__init__()
self.hidden1 = nn.Linear(1, 10)
self.hidden2 = nn.Linear(10, 10)
self.pred = nn.Linear(10, 1)
def forward(self, input):
out_hidden1 = nn.ReLU()(self.hidden1(input))
out_hidden2 = nn.ReLU()(self.hidden2(out_hidden1))
result = (self.pred(out_hidden2))
return result
net = LS()
print(net)
# 均方差损失函数
loss_func = nn.MSELoss()
# 定义优化器
optimizer = torch.optim.Adam(net.parameters())
# optimizer = torch.optim.SGD(net.parameters(), lr=0.1)
# 开始训练
losses = []
start_time = time.time()
for e in range(10000):
# for i in range(N):
# x = X[i]
# y = Y[i]
# pred = net(x)
# loss = loss_func(y, pred)
# optimizer.zero_grad()
# loss.backward()
# optimizer.step()
# # print(loss.item())
pred = net(X)
loss = loss_func(Y, pred)
print("Epoch: ", e, " ", loss.item())
optimizer.zero_grad()
loss.backward()
optimizer.step()
losses.append(loss.item())
end_time = time.time()
训练完成后是这个样子的：
注意左图中有一些绿色的点，这是按照独同分布产生的随机数据(可能会有和训练数据相同的吧，只是测试以下，不是重点)。代码中使用了Adam优化器，可以算是一种自适应的优化方法了。一般在实验的时候，可以先用Adam，RMSprop等跑一跑，有了对超参数的大体估计，可以换用SGD。因为Adam等虽然速度快，但不一定能找到最有点，而SGD虽然速度慢，但执着。
另一方面，我们无法在以上这个模型中找到曲线方程中的参数, 但神经网络就是奇妙地利用简单的线性方程加
ReLU函数，给出了我们这么好的预测结果。这难道就是人们把神经网络称为“黑盒”的原因吗？
3. Scipy 实现
最后，我们使用专门为优化而生的scipy.optimize模块来解决曲线拟合，或者说参数估计问题，简直是“稳，准，狠”。
这里测试了两个方法curve_fit和leastsq, 先贴上curve_fit的代码：
def ls_func(x_data, a, b, c):
return np.exp(a*x_data**2 + b*x_data + c)
popt, pcov = optimize.curve_fit(ls_func, X_np, Y_observed)
print(popt)
除了生成数据的代码，真正实现连同输出，总共四行代码，更要命的是，参数估计的结果相当准确。
在看看最小而成函数的结果（注意，上面的curve_fit实际上也是使用了最小二乘法，详细的可以看文档）。下面的代码借鉴了scipy数值优化与参数估计。
from scipy.optimize import leastsq
# https://blog.csdn.net/suzyu12345/article/details/70046826
def func(x, p):
""" 数据拟合所用的函数: A*sin(2*pi*k*x + theta) """
# A, k, theta = p
a, b, c = p
# return A*np.sin(2*np.pi*k*x+theta)
return np.exp(a*x**2 + b*x + c)
def residuals(p, y, x):
""" 实验数据x, y和拟合函数之间的差，p为拟合需要找到的系数 """
return 0.5 * (y - func(x, p))**2
p0 = [0.3, 0.1, -0.1]
# 调用leastsq进行数据拟合, residuals为计算误差的函数
# p0为拟合参数的初始值
# args为需要拟合的实验数据，也就是，residuals误差函数
# 除了P之外的其他参数都打包到args中
plsq = leastsq(residuals, p0, args=(Y_np, X_np))
print(plsq[0])
结束
现在看的文章，基本上每篇都涉及优化，我觉得要学习优化，对我这种数学功底不好的人来说，绝对不能仅仅看公式(因为看不懂)， 而是要结合代码去看。上面的代码还是太简单了，只属于松离合，点油门的起步作用。希望有高人能指点一二，或者推荐一些适合如入门的论文+开源代码以供学习。
