PyTorch Python 曲線擬合和優(yōu)化問題白話文
剛開始深度學習和SLAM卿捎，其中涉及到了很多的優(yōu)化問題配紫，查看一些文章，好多都是對優(yōu)化部分一句帶過午阵。網(wǎng)上的一些講機器學習和梯度下降的博文躺孝，大多用線性擬合講解梯度下降。今天底桂，用代碼實現(xiàn)了一個較為復雜函數(shù)的優(yōu)化問題植袍，份別用了不同的代碼方式，包括使用PyTorch自動求導戚啥，PyTorch搭建神經(jīng)網(wǎng)絡和Scipy.optimize模塊奋单。另外，這里不涉及具體的理論猫十，只提供了代碼示例和必要的注釋览濒。具體代碼請戳：PyTorch SGD,PyTorch NN ,Scipy Optimization (有些明明用的行內公式寫的，預覽的時候也正常拖云，這里看著全部變成行間公式了)
問題描述
該問題描述來自《視覺SLAM14講》贷笛。無償安利一下高翔老師這本，非常適合SLAM初學者入門宙项。以下這個優(yōu)化問題在《視覺SLAM14講》使用了C++代碼乏苦，配合Google的Ceres庫。
假設有一條滿足以下方程的曲線：
其中尤筐，為曲線的參數(shù),
為噪聲汇荐。很明顯，這是一個非線性模型盆繁。設有
個關于
的觀測數(shù)據(jù)掀淘，目的是根據(jù)這些數(shù)據(jù)擬合出曲線的參數(shù)。這個問題可以轉換為以下的優(yōu)化問題：
注意Ｓ桶骸８锫Α！這里待估計的參數(shù)是冕碟，而不是
!!!
開始擼代碼吧：
首先拦惋，我們產(chǎn)生一些觀測數(shù)據(jù)，這些觀測數(shù)據(jù)在三個不同方式的代碼中都是一樣的：
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
給個圖吧安寺，我們生成的數(shù)據(jù)就長這個樣子：
1. 梯度下降實現(xiàn)
這里不在詳細說明什么是梯度下降了厕妖，我認為核心公式就是：
即，不斷更新參數(shù),怎么更新呢我衬，用損失函數(shù)相對于
的偏導叹放，這個在代碼中有體現(xiàn)的饰恕。這個代碼中，我們使用了PyTorch的自動求導的功能井仰，這里只貼出了部分代碼埋嵌，完整代碼戳這里.
# 首先轉換以下numpy的數(shù)據(jù)類型，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)
# 這里的學習率也是多次實驗之后實驗出來的雹嗦，可以在
# 可以改為一個大點的值，然后在我下方標注“打斷點”的
# 地方打上斷點進行調試合是，可以發(fā)現(xiàn)a.grad的值非常大了罪，
# 這個可以根據(jù)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()))
注釋中說了，可以根據(jù)曲線的圖像理解蛾茉，我就把曲線的圖像貼出來：
可見讼呢，改曲線是非常非常陡峭的，
關于學習率的設定，在注釋中已經(jīng)說了吼鳞，可以按照這種方法打個斷電幕帆，我相信就一幕了然，總是聽說學習率很重要赖条。一些其它的博客，使用線性規(guī)劃常熙，好像隨便設置一個學習率纬乍，比如0.01, 0.1之類的，就可以的到很好的結果裸卫。用這個例子你試試仿贬，跑不死你......
下面給出迭代500個epoch 和 1000000 個epoch的截圖：
500：
1000000:
可見，迭代1000000的時候墓贿，得到的參數(shù)
總結：
- 梯度下降中，學習率的設置確實很重要
- 梯度下降不一定是最好的方法队伟，對初值的設置比較敏感
2. PyTorch 神經(jīng)網(wǎng)絡實現(xiàn)
這次穴吹，我搭建了一個有兩個隱藏層的全連接網(wǎng)絡，每個全連接層有10個神經(jīng)元嗜侮，具體看代碼吧:
# 網(wǎng)絡架構
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)
# 均方差損失函數(shù)
loss_func = nn.MSELoss()
# 定義優(yōu)化器
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()
訓練完成后是這個樣子的：
注意左圖中有一些綠色的點港令，這是按照獨同分布產(chǎn)生的隨機數(shù)據(jù)(可能會有和訓練數(shù)據(jù)相同的吧，只是測試以下锈颗，不是重點)顷霹。代碼中使用了Adam優(yōu)化器，可以算是一種自適應的優(yōu)化方法了击吱。一般在實驗的時候淋淀，可以先用Adam，RMSprop等跑一跑覆醇，有了對超參數(shù)的大體估計朵纷，可以換用SGD。因為Adam等雖然速度快叫乌，但不一定能找到最有點柴罐，而SGD雖然速度慢，但執(zhí)著憨奸。
另一方面革屠，我們無法在以上這個模型中找到曲線方程中的參數(shù), 但神經(jīng)網(wǎng)絡就是奇妙地利用簡單的線性方程加
ReLU函數(shù)，給出了我們這么好的預測結果排宰。這難道就是人們把神經(jīng)網(wǎng)絡稱為“黑盒”的原因嗎似芝？
3. Scipy 實現(xiàn)
最后，我們使用專門為優(yōu)化而生的scipy.optimize模塊來解決曲線擬合板甘，或者說參數(shù)估計問題党瓮，簡直是“穩(wěn)，準盐类，狠”寞奸。
這里測試了兩個方法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)
除了生成數(shù)據(jù)的代碼，真正實現(xiàn)連同輸出在跳，總共四行代碼枪萄，更要命的是，參數(shù)估計的結果相當準確猫妙。
在看看最小而成函數(shù)的結果（注意瓷翻，上面的curve_fit實際上也是使用了最小二乘法，詳細的可以看文檔）。下面的代碼借鑒了scipy數(shù)值優(yōu)化與參數(shù)估計齐帚。
from scipy.optimize import leastsq
# https://blog.csdn.net/suzyu12345/article/details/70046826
def func(x, p):
""" 數(shù)據(jù)擬合所用的函數(shù): 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):
""" 實驗數(shù)據(jù)x, y和擬合函數(shù)之間的差妒牙，p為擬合需要找到的系數(shù) """
return 0.5 * (y - func(x, p))**2
p0 = [0.3, 0.1, -0.1]
# 調用leastsq進行數(shù)據(jù)擬合, residuals為計算誤差的函數(shù)
# p0為擬合參數(shù)的初始值
# args為需要擬合的實驗數(shù)據(jù)，也就是对妄，residuals誤差函數(shù)
# 除了P之外的其他參數(shù)都打包到args中
plsq = leastsq(residuals, p0, args=(Y_np, X_np))
print(plsq[0])
結束
現(xiàn)在看的文章湘今，基本上每篇都涉及優(yōu)化，我覺得要學習優(yōu)化饥伊，對我這種數(shù)學功底不好的人來說象浑，絕對不能僅僅看公式(因為看不懂)， 而是要結合代碼去看琅豆。上面的代碼還是太簡單了愉豺，只屬于松離合，點油門的起步作用茫因。希望有高人能指點一二蚪拦，或者推薦一些適合如入門的論文+開源代碼以供學習。
