【手撕版】MXNet應用之線性回歸

近日 Keras 作者 Fran?ois Chollet 在 Twitter 上公布了一項調查：過去三月中 ArXiv 上（截至 3 月 7 日）被提及最多的幾大開源框架。可以看到，Tensorflow穩坐首位，

同期Fran?ois Chollet還發布了github上開源框架的人氣排名，可以看出MXNet的排名顯著上升（從第六到第三）。

MXNet 是亞馬遜（Amazon）[2] 選擇的深度學習庫。它擁有類似於 Theano 和 TensorFlow 的數據流圖，為多 GPU 配置提供了良好的配置，有著類似於 Lasagne 和 Blocks 更高級別的模型構建塊，可以在多數硬體上運行（包括手機）。MXNet 允許混合符號編程和命令式編程，從而最大限度提高效率和生產力，其核心是一個動態的依賴調度，能夠自動並行符號和命令的操作。

科賽社區近日收羅了MXNet系列教程，感謝fksad的貢獻，欲深入學習的同學可登錄網站，Fork後練習。

《mxnet_ch02_linear_regression》

作者：fksad

原文鏈接：科賽 - Kesci.com

參考課程鏈接：Tsinghua Open Source Mirror

儘管強大的深度學習框架可以減少大量重複性工作，但若過於依賴它提供的便利，你就會很難深入理解深度學習是如何工作的。因此，我們的第一個教程是如何只利用ndarray和autograd來實現一個線性回歸的訓練。

線性回歸

給定一個數據點集合 X和對應的目標值y，線性回歸模型的目標就是找到一條由方向向量w和截距b確定的線，使其上的點能近似表示數據集的樣本X[i] 和y[i]。用數學符號來表示就是：

優化的目標是要最小化所有數據點上的平方誤差（樣本點到擬合直線距離的平方和）

線性回歸模型是最簡單、但也是最有用的神經網路。一個神經網路就是一個由節點（神經元）和有向邊組成的集合。我們一般把一些節點組成層，每一層先從下面一層的節點獲取輸入，然後輸出給上面的層使用。要計算一個節點值，我們需要將輸入節點值加權求和（權數值即 w），然後再套上一個激活函數（activation function）。

對於線性回歸而言，它是一個兩層神經網路（或者按照西瓜書的說法，是一個單隱層的神經網路，隱藏層其實就是輸出層），其中第一層是（下圖橙色點）輸入，每個節點對應輸入數據點的一個維度，第二層是單輸出節點（下圖綠色點），它使用身份函數（ f(x)=x）作為激活函數。

創建數據集

這裡我們使用一個數據集來作簡要說明。我們使用如下方法來生成數據；隨機數值 $X[i]$，其相應的標註為 y[i]：

使用數學符號表示：

這裡噪音服從均值為 0 , 標準差為0.01的正態分布。

from mxnet import ndarray as ndfrom mxnet import autograd# 設定x的維數和樣本量num_inputs = 2num_examples = 1000# 設定待估參數真值true_w = [2, -3.4]true_b = 4.2# 生成模擬樣本X = nd.random_normal(shape=(num_examples, num_inputs)) # 2維向量y = true_w[0] * X[:, 0] + true_w[1] * X[:, 1] + true_b # 1維向量（標量）y += .01 * nd.random_normal(shape=y.shape)print(X[0], y[0])

如果有興趣，可以使用安裝包中已包括的Python繪圖包matplotlib，生成第二個特徵值 (X[:, 1]) 和目標值Y的散點圖，更直觀地觀察兩者間的關係。如果線性關係不明顯，可以考慮改變模型設定。

import matplotlib.pyplot as pltplt.scatter(X[:, 1].asnumpy(),y.asnumpy())plt.show()

數據讀取

當我們開始訓練神經網路的時候，我們需要不斷讀入數據塊。這裡我們定義一個函數 —— 它每次返回 batch_size個無放回的隨機的樣本（X）和對應的目標(Y)。我們通過python的 yield來構造一個迭代器。也就是在一輪訓練中，我們用樣本量為batch_size的無重複子樣本進行多次基於隨機梯度下降的優化。如果想用全樣數據一次性訓練模型，不妨把這個變數設置的和樣本量一樣大。

import randomdef data_iter(batch_size = 10): batch_size: int; 每次返回的樣本來那個 # 產生一個隨機索引 idx = list(range(num_examples)) random.shuffle(idx) for i in range(0, num_examples, batch_size): j = nd.array(idx[i:min(i+batch_size, num_examples)]) # 將索引切分成每 batch_size 個一份 yield nd.take(X, j), nd.take(y, j) # 與 X[j], y[j] 等價# 隨機讀取一個子集for data, label in data_iter(): print(data, label) break

初始化模型參數

同第一張第二部分介紹的步驟一樣，我們先隨機初始化模型參數，之後訓練時我們需要對這些參數求導來更新它們的值，使損失盡量減小；因此我們需要創建它們的梯度（分配內存）。

# 有 num_inputs 個權重和 1 個截距w = nd.random_normal(shape=(num_inputs, 1))b = nd.zeros((1,))params = [w, b]# 給係數列表的每個元素分配存放梯度的內存for param in params: param.attach_grad()

定義模型

之後我們需要定義模型形式，等同於第一部分的顯示定義待估係數的函數形式。線性模型中，就是輸入的數據矩陣乘以權重再加上截距

def net(X): return nd.dot(X, w) + b

損失函數

損失函數是我們需要最小化的目標，我們常用平方誤差（均方誤差）來衡量預測值（擬合的直線）和真實值之間的差距。

def square_loss(yhat, y): # 注意這裡我們把y變形成yhat的形狀來避免矩陣形狀的自動轉換 return (yhat - y.reshape(yhat.shape)) ** 2

優化

之後的工作就是帶入已有數據並利用梯度下降法求得使得損失函數最小的係數組合了。

當然，線性回歸的係數其實有顯式解

但絕大部分模型並沒有。所以我們這裡通過隨機梯度下降來求解，保證模型係數估計過程的一致性。

每一步，我們將模型參數沿著梯度的反方向走一段距離，這個距離一般叫學習率（learning rate） lr。（我們會之後一直使用這個函數，我們將其保存在utils.py）。訓練通常需要迭代數據數次，在這裡使用epochs表示迭代總次數；一次迭代中，我們每次隨機讀取固定數個數據點，計算梯度並更新模型參數。看看這個優化演算法能否實現我們的目的。

def SGD(params, lr): for param in params: param[:] = param - lr * param.grad # param[:] 代表對被遍歷的 params 做原地修改

# 模型函數def real_fn(X): # w = [1, -3.4], b = 4.2 return 2 * X[:, 0] - 3.4 * X[:, 1] + 4.2# 繪製損失隨訓練次數降低的折線圖，以及預測值和真實值的散點圖def plot(losses, X, sample_size=100): losses: list; 誤差序列 X: nd.array; 訓練集數據 sample_size: int; 訓練集樣本量 xs = list(range(len(losses))) f, (fg1, fg2) = plt.subplots(1, 2) # 子圖1是訓練誤差和輪數的折線圖 fg1.set_title(Loss during training) fg1.plot(xs, losses, -r) # 子圖2是估計的樣本和真值的對比散點圖 fg2.set_title(Estimated vs real function) fg2.plot(X[:sample_size, 1].asnumpy(), net(X[:sample_size, :]).asnumpy(), or, label=Estimated) fg2.plot(X[:sample_size, 1].asnumpy(), real_fn(X[:sample_size, :]).asnumpy(), *g, label=Real) fg2.legend() plt.show()

epochs = 5learning_rate = .001niter = 0losses = []moving_loss = 0smoothing_constant = .01# 訓練5輪for e in range(epochs): total_loss = 0 # 利用 data_iter 函數抽取選聯數據 for data, label in data_iter(batch_size=10): # 顯示定義這批數據的損失函數 with autograd.record(): output = net(data) loss = square_loss(output, label) # 求導並利用梯度下降做優化 loss.backward() SGD(params, learning_rate) # 累加得到這一輪訓練的總誤差 total_loss += nd.sum(loss).asscalar() # 記錄每讀取一批數據點後，損失的移動平均值的變化； niter += 1 curr_loss = nd.mean(loss).asscalar() # 計算加權移動平均後的誤差 moving_loss = (1 - smoothing_constant) * moving_loss + (smoothing_constant) * curr_loss # 利用移動平均修正誤差 est_loss = moving_loss/(1-(1-smoothing_constant)**niter) # 利用整個數據集完成一輪訓練後繪圖 losses.append(est_loss) print("Epoch %s, batch %s. Moving avg of loss: %s. Average loss: %f" % (e, niter, est_loss, total_loss/num_examples)) plot(losses, X)

結論

綜上，利用NDArray和autograd就可以很容易實現的一個模型。對比係數估計值和真值，我們發現差別很小，可見演算法整體很合理。

print(true_w, w)print(true_b, b)

線性回歸 - gluon版

使用 mxnet 的高層抽象包可以很方便地訓練一個線性模型。

數據導入

此處我們沿用手撕版的數據集，利用 gluon 的內置函數完成建模。

from mxnet import gluon# 利用 data 模塊的 ArrayDataset 函數完成數據導入dataset = gluon.data.ArrayDataset(X, y)# 同上，利用 data 模塊的 Dataloader 函數完成上文 data_iter 的工作batch_size = 10data_iter = gluon.data.DataLoader(dataset, batch_size, shuffle=True)# 二者一毛一樣for data, label in data_iter: print(data, label) break

定義模型

手撕版中，我們需要先聲明模型參數，然後再使用它們來構建模型。在gluon中，我們把線性回歸模型視作神經網路的特例，利用gluon提供大量預定義的層，我們只需要關注使用哪些層來構建模型。

線性模型使用對應的Dense層；之所以稱為dense層，是因為輸入的所有節點都與後續的節點相連。在這個例子中僅有一個輸出，但在大多數後續章節中，我們會用到具有多個輸出的網路（比如分類器）。

我們之後還會介紹如何構造任意結構的神經網路，但對於初學者來說，構建模型最簡單的辦法是利用 Sequential 串聯所有層。輸入數據之後，Sequential 會依次執行每一層，並將前一層的輸出，作為輸入提供給後面的層。

我們要做的是：

1. 定義一個空的模型
2. 加入一個Dense層，同時必須定義輸出的維度，線性模型中這個參數是 1
3. 此處我們不需要顯示的指定輸入維度，因為線性回歸是只有輸入層和輸出層的神經網路，輸入層維度就是訓練集的維度
4. 初始化參數 - 此處不傳入參數，利用默認設置初始化
5. 設定損失函數 - gluon 內置平方損失
6. 進行優化 - 創建一個Trainer的實例，將學習率等參數參數傳遞給它就行。
7. 訓練模型 - 基本同手撕版

使用gluon使模型訓練過程更為簡潔。我們不需要挨個定義相關參數、損失函數，也不需使用隨機梯度下降。gluon的抽象和便利的優勢將隨著我們著手處理更多複雜模型的愈發顯現。不過在完成初始設置後，訓練過程本身和前面沒有太多區別，唯一的不同在於我們不再是調用SGD，而是trainer.step來更新模型（此處一併省略之前繪製損失變化的折線圖和散點圖的過程，有興趣的同學可以自行嘗試）。

# 定義空模型net = gluon.nn.Sequential()# 添加 Dense 層 net.add(gluon.nn.Dense(1))# 初始化參數net.initialize()# 指定損失函數square_loss = gluon.loss.L2Loss()# 指定優化演算法和學習率trainer = gluon.Trainer( net.collect_params(), sgd, {learning_rate: 0.1}) # 學習率大了不少# 進行 epochs 輪訓練epochs = 5batch_size = 10for e in range(epochs): total_loss = 0 for data, label in data_iter: with autograd.record(): output = net(data) loss = square_loss(output, label) loss.backward() trainer.step(batch_size) total_loss += nd.sum(loss).asscalar() print("Epoch %d, average loss: %f" % (e, total_loss/num_examples))# 訓練好的參數結果就是dense層的權重值dense = net[0]est_w = dense.weight.data()[0]est_b = dense.bias.data()[0]# 和真實值對比，沒毛病print(true_b, est_b)print(true_w, est_w)

思考題

• 在訓練的時候，為什麼我們用了比前面要大10倍的學習率呢？（提示：可以嘗試運行 help(trainer.step)來尋找答案。）
• 如何拿到weight的梯度呢？（提示：嘗試 help(dense.weight)）

推薦閱讀：