家樂的深度學習筆記「3」 – 線性回歸

• 2020 年 3 月 26 日
• 筆記

線性回歸

線性回歸的輸出是一個連續值，適用於回歸問題，如房價預測、氣溫、銷售額等連續值問題。

基本要素

以最基礎的房價預測為例，假設價格只取決於面積（平方米）和房齡（年），下面來介紹大多數深度學習模型的基本要素和表示方法。

模型

設房屋面積為x，房齡為x，售出價格為y。我們需要建立基於輸入x和x來計算輸出y的表達式，也就是模型。
y = xw+ xw + b
其中w1和w2是權重，b是偏差，且均為標量。它們是線性回歸模型的參數。y是線性回歸對真實價格y的預測或估計。通常允許它們之間有一定誤差。

模型訓練

構建好模型後，需要通過數據來尋找特定的模型參數值，使模型在數據上的誤差儘可能小。這個過程稱為模型訓練。

訓練數據

收集到的一系列真實數據被稱作訓練數據集或訓練集。以房價數據集為例，一棟房屋被稱為一個樣本，其真實出售價格叫作標籤，用來預測標籤的兩個因素叫做特徵，特徵用來表徵樣本的特點。
假設採集的樣本數為n，索引為i的樣本的特徵為x和x，標籤為y。對於索引為i的房屋，線性回歸模型的房屋價格預測表達式為
y(i) = x(i)w+ x(i)w + b

損失函數

平方函數通常用來衡量價格預測值與真實值之間的誤差，其在評估索引為i的樣本誤差的表達式為

給定訓練數據集，這個誤差只與模型參數相關，因此將其記為以模型參數為參數的函數。
在機器學習里，將衡量誤差的函數稱為損失函數。這裡使用的平方誤差函數也稱為平方損失。

通常使用訓練數據集中所有樣本誤差的平均來衡量模型預測的質量，即

在模型訓練中，我們希望找出一組模型參數，記為$$w_1^, w_2^, b^*$$，來使訓練樣本平均損失最小：

優化算法

當模型和損失函數形式較為簡單時，上面的誤差最小化問題的解可以直接用公式表達出來。這類解叫作解析解(analytical solution)。
然而大多數深度學習模型並沒有解析解，只能通過優化算法有限次迭代模型參數來儘可能降低損失函數的值。這類解叫作數值解(numerical solution)。
以下為抄書內容：
在求解數值解的優化算法中，小批量隨機梯度下降(mini-batch stochastic gradient descent)在深度學習中被廣泛使用。其先選取一組模型參數的初始值，如隨機選取；接下來對參數進行多次迭代，使每次迭代都可能降低損失函數的值。在每次迭代中，先隨機均勻採樣一個由固定數目訓練數據樣本所組成的小批量（mini-batch）(mathcal{B})，然後求小批量中數據樣本的平均損失有關模型參數的導數（梯度），最後用此結果與預先設定的一個正數的乘積作為模型參數在本次迭代的減小量。

在訓練本節線性回歸模型過程中，模型的每個參數將作如下迭代：

在上式中，(|mathcal{B}|)代表每個小批量中的樣本個數（批量大小，batch size），(eta)稱作學習率（learning rate）並取正數。需要強調的是，這裡的批量大小和學習率的值是人為設定的，並不是通過模型訓練學出的，因此叫作超參數（hyperparameter）。我們通常所說的「調參」指的正是調節超參數，例如通過反覆試錯來找到超參數合適的值。
在少數情況下，超參數也可以通過模型訓練學出。

模型預測

模型訓練完成後，將模型參數(w_1, w_2, b)在優化算法停止時的值記為(hat{w}_1, hat{w}_2, hat{b})，這裡得到的不一定是最小化損失函數的最優解，而是對最優解的一個近似。然後，就可以用學出的線性回歸模型來估算訓練數據集以外任意一個樣本的標籤了。
這裡的估算也叫作模型預測、模型推斷或模型測試。

表示方法

下面解釋線性回歸與神經網絡的聯繫，以及線性回歸的矢量計算表達式。

神經網絡圖

在深度學習中，可以使用神經網絡圖直觀地表現模型結構。

神經網絡圖隱去了模型參數權重和偏差。
上述神經網絡輸入分別為x和x，因此輸入層的輸入個數為2。輸入個數也叫特徵數或特徵向量維度。
上述神經網絡中，直接使用神經網絡的輸出o作為線性回歸的輸出，即$$hat{y} = o$$。
由於輸入層並不涉及計算，按照慣例，上述神經網絡的層數為1。所以，線性回歸是一個單層神經網絡，輸出層中負責計算o的單元又叫神經元。在線性回歸中，o的計算依賴於x和x。也就是說，輸出層中的神經元和輸入層中各個輸入完全連接。因此，這裡的輸出層又叫全連接層(fully-connected layer)，或稠密層(dense layer)。

矢量計算表達式

在模型訓練或預測時，常常會同時處理多個數據樣本並用到矢量計算。
矢量加法快於標量加法，應該儘可能使用矢量運算以提升效率。

廣義上講，當數據樣本數為(n)，特徵數為(d)時，線性回歸的矢量計算表達式為
(boldsymbol{hat{y}} = boldsymbol{X} boldsymbol{w} + b)
其中模型輸出(boldsymbol{hat{y}} in mathbb{R}^{n times 1})， 批量數據樣本特徵(boldsymbol{X} in mathbb{R}^{n times d})，權重(boldsymbol{w} in mathbb{R}^{d times 1})， 偏差(b in mathbb{R})。相應地，批量數據樣本標籤(boldsymbol{y} in mathbb{R}^{n times 1})。設模型參數(boldsymbol{theta} = [w_1, w_2, b]^top)，我們可以重寫損失函數為

小批量隨機梯度下降的迭代步驟將相應地改寫為

其中梯度是損失有關3個為標量的模型參數的偏導數組成的向量：

從零實現

這裡不利用高級框架，只用NDArray和autograd實現線性回歸。

from mxnet import gpu, nd, autograd  import matplotlib.pyplot as plt  import random  

這裡可以順帶固定下隨機數種子。

random.seed(233)  

生成數據集

設訓練數據集樣本數為1000，輸入個數（特徵數）為2。給定隨機生成的批量樣本特徵(boldsymbol{X} in mathbb{R}^{1000 times 2})，使用線性回歸模型真實權重(boldsymbol{w} = [2, -3.4]^top)和偏差(b = 4.2)，以及一個隨機噪聲項(epsilon)來生成標籤，其中噪聲項服從均值為0、標準差為0.01的正態分佈。噪聲代表了數據集中無意義的干擾。
因為使用了一台四卡2080Ti的服務器學習，為了盡量不對別人帶來困擾，這裡使用第四張卡，MXNet不會像tf一樣默認佔用整張卡的顯存，而是一點點佔用，這點值得表揚。

ctx = gpu(3)    num_examples = 1000  num_inputs = 2  true_w = [2, -3.4]  true_b = 4.2    features = nd.random.normal(scale=1, shape=(num_examples, num_inputs), ctx=ctx)  labels = true_w[0] * features[:, 0] + true_w[1] * features[:, 1] + true_b  labels += nd.random.normal(scale=0.01, shape=labels.shape, ctx=ctx)  

看一眼生成的訓練數據集部分特徵與標籤的樣子：

features[:5], labels[:5]  

(   [[ 0.34113222  0.76659095]    [ 0.6508758  -1.5585263 ]    [ 1.9786316  -0.8465744 ]    [-0.05428567  0.64376295]    [-1.7093595   0.12594718]]   <NDArray 5x2 @gpu(3)>,   [ 2.2640598  10.818717   11.050849    1.8893591   0.35090753]   <NDArray 5 @gpu(3)>)  

生成一個散點圖可以直觀的觀察兩者之間的線性關係。

plt.scatter(features[:, 1].asnumpy(), labels.asnumpy(), 1);  

改一下生成格式為svg，即矢量圖，這樣可以省下圖片的空間。
如果不保存的話，畫圖語句末尾的分號可以省去打印該圖的類型信息，只顯示圖。

from IPython import display  display.set_matplotlib_formats('svg')  point_size = 1  plt.scatter(features[:, 1].asnumpy(), labels.asnumpy(), point_size);  plt.savefig('散點圖.svg')  

書里改變了圖的大小，我認為沒什麼必要，這樣的大小在屏幕正好觀察。

# plt.rcParams['figure.figsize'] = (3.5, 2.5)  # plt.scatter(features[:, 1].asnumpy(), labels.asnumpy(), 1);  

讀取數據集

這裡使用到了python的生成器知識點，詳見我的博客Python高級特性介紹 – 迭代的99種姿勢與協程里有詳細介紹。
按小批讀取訓練數據集，這裡每次都判斷了min，感覺有點多餘但是自己又不知道怎麼更好的修改。

# count = []  def data_iter(batch_size, features, labels, ctx):      num_examples = len(features)      indices = list(range(num_examples))      random.shuffle(indices)      for i in range(0, num_examples, batch_size):          j = nd.array(indices[i:min(i + batch_size, num_examples)], ctx=ctx)  #         count.append(j)          yield features.take(j), labels.take(j)  # print(count[-1])  

這個函數建好之後讓我們試試效果：

batch_size = 10  for X, y in data_iter(batch_size, features, labels, ctx):      print(X, y)      break  

[[ 0.1555565  -0.89841115]   [ 0.19113457 -0.23471509]   [-0.1649106  -1.0402957 ]   [ 0.59347206 -0.04631801]   [-1.473519    0.7591845 ]   [ 1.1765232  -0.4977442 ]   [ 0.35199913 -0.9372813 ]   [ 0.9608316   0.19639738]   [-1.253151   -0.55073833]   [-0.09366591 -0.84877795]]  <NDArray 10x2 @gpu(3)>  [ 7.580011   5.38007    7.406578   5.5393977 -1.3400362  8.245509    8.093137   5.4496117  3.5762053  6.8978252]  <NDArray 10 @gpu(3)>  

初始化模型參數

將模型的權重初始化成均值為0、標準差為0.01的正態隨機數，偏差則初始化成0。

w = nd.random.normal(scale=0.01, shape=(num_inputs, 1), ctx=ctx)  b = nd.zeros(shape=(1, ), ctx=ctx)  

在模型訓練中，需要對這些參數求梯度來迭代參數，所以在這裡需要創建他們的梯度。

w.attach_grad()  b.attach_grad()  

定義模型

模型就是我們期望擬合的出來的函數，在這即線性回歸的矢量計算表達式的實現，下面使用dot做矩陣乘法。

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

定義損失函數

使用前述平方損失來定義線性回歸的損失函數，在這裡真實值y要reshape成預測值y_hat的形狀，實測預測值為二維數組[batch_size, 1]，而真實值為一維數組[batch_size]，這是對書後練習一的解答。

def squared_loss(y_hat, y):  #     print(y_hat, y, y.reshape(y_hat.shape))      return (y_hat - y.reshape(y_hat.shape)) ** 2 / 2  

定義優化算法

以下sgd函數實現了前述小批量隨機梯度下降算法。它通過不斷迭代模型參數來優化損失函數。這裡自動求梯度模塊計算得來的梯度是一個批量樣本的梯度和，在筆記「2」中有所介紹，將其除以批量大小即得出梯度的平均值用於迭代。

def sgd(params, lr, batch_size):      for param in params:          param -= lr * param.grad / batch_size  

訓練模型

訓練需要確定（學習速率、批量大小、迭代輪數、模型、損失函數、優化算法），這樣就可以開始我們的訓練，俗稱跑模型。
前述筆記「2」中提到，由於(l)並不是一個標量，運行(l.backward())將對(l)中元素求和得到新的變量，再求該變量有關模型參數的梯度。
在一個迭代周期(epoch)中，將完整遍歷一遍data_iter函數，即對訓練數據集中所有樣本都使用一次，如果樣本數不能被批量大小整除，最後一個批次的樣本數目會減少，這是對書後練習三的回答。
這裡的迭代周期個數num_epochs和學習速率lr都是超參數，在實踐中，大多數超參數都需要通過反覆試錯來不斷調節，俗稱調參。

lr = 0.03  num_epochs = 5  net = linreg  loss = squared_loss    for epoch in range(num_epochs):      for X, y in data_iter(batch_size, features, labels, ctx):          with autograd.record():              l = loss(net(X, w, b), y)          l.backward()          sgd([w, b], lr, batch_size)      train_l = loss(net(features, w, b), labels)      print('epoch %d, loss %f' % (epoch, train_l.mean().asnumpy()))  

epoch 0, loss 0.027449  epoch 1, loss 0.000099  epoch 2, loss 0.000052  epoch 3, loss 0.000051  epoch 4, loss 0.000051  

訓練完後就能看到我們訓練的數據與真實數據的距離了。

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

[2, -3.4]  [[ 2.0001543]   [-3.4000373]]  <NDArray 2x1 @gpu(3)>  4.2  [4.1994386]  <NDArray 1 @gpu(3)>  

全部代碼

這裡給出我的notebook的全部代碼。
線性回歸從零實現.html

簡潔實現

上面手寫的適合理解原理，但是人和動物的區別就是人會使用工具，可以在前人的基礎上進行自己的工作，MXNet就提供了Gluon接口更方便的實現線性回歸的訓練。

from mxnet.gluon import nn, loss as gloss, data as gdata  from mxnet import gpu, nd, autograd, init, gluon  import matplotlib.pyplot as plt  import random  

設一下隨機數種子，方便復現。

random.seed(233)  

生成數據集

和上面一樣，就不多解釋了。

ctx = gpu(3)    num_examples = 1000  num_inputs = 2  true_w = [2, -3.4]  true_b = 4.2    features = nd.random.normal(scale=1, shape=(num_examples, num_inputs), ctx=ctx)  labels = true_w[0] * features[:, 0] + true_w[1] * features[:, 1] + true_b  labels += nd.random.normal(scale=0.01, shape=labels.shape, ctx=ctx)  

查看一下生成的訓練數據集的樣子。

features[:5], labels[:5]  

(   [[ 0.34113222  0.76659095]    [ 0.6508758  -1.5585263 ]    [ 1.9786316  -0.8465744 ]    [-0.05428567  0.64376295]    [-1.7093595   0.12594718]]   <NDArray 5x2 @gpu(3)>,   [ 2.2640598  10.818717   11.050849    1.8893591   0.35090753]   <NDArray 5 @gpu(3)>)  

圖像方式查看：

from IPython import display  display.set_matplotlib_formats('svg')  point_size = 1  plt.scatter(features[:, 1].asnumpy(), labels.asnumpy(), point_size);  plt.savefig('散點圖.svg')  

讀取數據集

這裡就可以使用Gluon提供的data包來讀取數據了。

batch_size = 10  dataset = gdata.ArrayDataset(features, labels)  data_iter = gdata.DataLoader(dataset, batch_size, shuffle=True)  

需要注意的是，這裡X是來自我們指定的GPU，而y來自了CPU，雖然不明白這是feature還是bug，但是官方論壇上有相應討論：論壇鏈接，比較可惜的是，在發文時間2020-03-26，官網論壇炸了，很奇怪為什麼亞馬遜背書的MXNet的aws服務器會掛掉，從谷歌快照上得到的解決方案還是把數據重新拷貝到GPU上，well。

for X, y in data_iter:      print(X, y)      break  

[[ 0.5275261  -1.0045179 ]   [ 0.44837257 -0.43323204]   [-0.01466356 -0.48579317]   [ 1.1003594   0.7395018 ]   [-0.22922567 -1.9637256 ]   [ 0.4390832  -0.60523504]   [-0.75654405 -1.4318156 ]   [-1.1557755   0.4982119 ]   [ 0.5148314   0.9012168 ]   [ 0.42361584 -1.2854003 ]]  <NDArray 10x2 @gpu(3)>  [ 8.658156    6.5668592   5.829563    3.9018939  10.418489    7.1303806    7.569174    0.18902431  2.1697466   9.418344  ]  <NDArray 10 @cpu(0)>  

定義模型

Gluon提供的大量預定義層讓我們可以只關注使用哪些層來構造模型，是不是和Keras看起來很像（
有一種在搭積木的感覺。

net = nn.Sequential()  

線性回歸作為單層神經網絡，輸出層的神經元和輸入層中各個輸入完全連接。因此，線性回歸的輸出層又叫全連接層，這裡添加一層全連接層，並定義該層的輸出個數為1。

net.add(nn.Dense(1))  

初始化模型參數

init模塊是MXNet提供的，通過其指定權重參數每個元素的初始化值。

net.initialize(init.Normal(sigma=0.01), ctx=ctx)  

定義損失函數

Gluon里提供了loss模塊，在這裡我們直接使用L2範數損失，其實L2就是前面用的平方損失，範數的定義就是這樣，比如L3範數就是真實值與預測值差的立方的立方根。

loss = gloss.L2Loss()  

定義優化算法

優化算法依然選擇sgd，在這裡直接使用字符串指定即可。該算法用來迭代net實例所有通過add函數嵌套的層所包含的全部參數。這些參數使用collect_params()函數獲取。

trainer = gluon.Trainer(net.collect_params(), 'sgd', {'learning_rate': 0.03})  

訓練模型

使用Gluon訓練時，使用Trainer實例的step函數來迭代訓練模型。

num_epochs = 3  for epoch in range(1, num_epochs + 1):      for X, y in data_iter:          with autograd.record():              l = loss(net(X), y.as_in_context(ctx))          l.backward()          trainer.step(batch_size)      l = loss(net(features), labels)      print('epoch %d, loss %f' % (epoch, l.mean().asnumpy()))  

epoch 1, loss 0.027621  epoch 2, loss 0.000095  epoch 3, loss 0.000051  

從net中獲得需要的層後即可訪問其權重和偏差，另外也可以獲取其梯度，這是對練習三的回答。

dense = net[0]  print(true_w, dense.weight.data())  print(true_b, dense.bias.data())  dense.weight.grad()  

[2, -3.4]  [[ 2.000109  -3.3995216]]  <NDArray 1x2 @gpu(3)>  4.2  [4.1992593]  <NDArray 1 @gpu(3)>  Out[14]:    [[0.02693137 0.00585834]]  <NDArray 1x2 @gpu(3)>  

練習

練習一說，如果(l)替換為mean()，trainer.step(batch_size)相應的需要調整為trainer.step(1)，這是因為迭代求得的梯度會被(1/batch_size)正則化（Gradient will be normalized by 1 / batch_size）。
並且如果需要獲取其中的梯度值，也可以自己調用allreduce_grads() + update()實現，step()函數本身就相當於這兩個函數的組合，並且需要在autograde.backward()之後、record()作用域之外調用。

net.initialize(init.Normal(sigma=0.01), ctx=ctx, force_reinit=True)  

num_epochs = 3  for epoch in range(1, num_epochs + 1):      for X, y in data_iter:          with autograd.record():              l = loss(net(X), y.as_in_context(ctx)).mean()          l.backward()          trainer.step(1)      l = loss(net(features), labels)      print('epoch %d, loss %f' % (epoch, l.mean().asnumpy()))  

epoch 1, loss 0.027576  epoch 2, loss 0.000097  epoch 3, loss 0.000051  

這裡查看的結果和剛才是幾乎相同的。

dense = net[0]  print(true_w, dense.weight.data())  print(true_b, dense.bias.data())  

[2, -3.4]  [[ 2.0000014 -3.3993788]]  <NDArray 1x2 @gpu(3)>  4.2  [4.1991076]  <NDArray 1 @gpu(3)>  

另外關於練習二，官方文檔的鏈接在這gluon.loss提供的損失函數、init提供的初始化方法，筆記用到了再另外介紹。

全部代碼

這裡給出我的notebook的全部代碼。
線性回歸簡潔實現.html