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

给个图吧，我们生成的数据就长这个样子：

Real Data

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()))

注释中说了，可以根据曲线的图像理解，我就把曲线的图像贴出来：

Curve
可见，改曲线是非常非常陡峭的， 500
1000000:
1000000
可见，迭代1000000的时候，得到的参数 NN

注意左图中有一些绿色的点，这是按照独同分布产生的随机数据(可能会有和训练数据相同的吧，只是测试以下，不是重点)。代码中使用了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])

结束

现在看的文章，基本上每篇都涉及优化，我觉得要学习优化，对我这种数学功底不好的人来说，绝对不能仅仅看公式(因为看不懂)， 而是要结合代码去看。上面的代码还是太简单了，只属于松离合，点油门的起步作用。希望有高人能指点一二，或者推荐一些适合如入门的论文+开源代码以供学习。