目錄

• 線性迴歸
  • 基本要素
  • 表示方法
  • 從零實現
  • 簡潔實現

線性迴歸

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

基本要素

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

模型

設房屋面積為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的樣本誤差的表示式為
\[\ell^{(i)}(w_1, w_2, b) = \frac{1}{2} \left(\hat{y}^{(i)} - y^{(i)}\right)^2\]

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

通常使用訓練資料集中所有樣本誤差的平均來衡量模型預測的質量，即
\[\ell(w_1, w_2, b) =\frac{1}{n} \sum_{i=1}^n \ell^{(i)}(w_1, w_2, b) =\frac{1}{n} \sum_{i=1}^n \frac{1}{2}\left(x_1^{(i)} w_1 + x_2^{(i)} w_2 + b - y^{(i)}\right)^2\]
在模型訓練中，我們希望找出一組模型引數，記為\[w_1^*, w_2^*, b^*\]，來使訓練樣本平均損失最小：

\[w_1^*, w_2^*, b^* = \operatorname*{argmin}_{w_1, w_2, b}\ \ell(w_1, w_2, b)\]

優化演算法

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

在訓練本節線性迴歸模型過程中，模型的每個引數將作如下迭代：
\[\begin{aligned} w_1 &\leftarrow w_1 - \frac{\eta}{|\mathcal{B}|} \sum_{i \in \mathcal{B}} \frac{ \partial \ell^{(i)}(w_1, w_2, b) }{\partial w_1} = w_1 - \frac{\eta}{|\mathcal{B}|} \sum_{i \in \mathcal{B}}x_1^{(i)} \left(x_1^{(i)} w_1 + x_2^{(i)} w_2 + b - y^{(i)}\right),\\ w_2 &\leftarrow w_2 - \frac{\eta}{|\mathcal{B}|} \sum_{i \in \mathcal{B}} \frac{ \partial \ell^{(i)}(w_1, w_2, b) }{\partial w_2} = w_2 - \frac{\eta}{|\mathcal{B}|} \sum_{i \in \mathcal{B}}x_2^{(i)} \left(x_1^{(i)} w_1 + x_2^{(i)} w_2 + b - y^{(i)}\right),\\ b &\leftarrow b - \frac{\eta}{|\mathcal{B}|} \sum_{i \in \mathcal{B}} \frac{ \partial \ell^{(i)}(w_1, w_2, b) }{\partial b} = b - \frac{\eta}{|\mathcal{B}|} \sum_{i \in \mathcal{B}}\left(x_1^{(i)} w_1 + x_2^{(i)} w_2 + b - y^{(i)}\right). \end{aligned}\]
在上式中，\(|\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\)，我們可以重寫損失函式為
\[\ell(\boldsymbol{\theta})=\frac{1}{2n}(\boldsymbol{\hat{y}}-\boldsymbol{y})^\top(\boldsymbol{\hat{y}}-\boldsymbol{y})\]
小批量隨機梯度下降的迭代步驟將相應地改寫為
\[\boldsymbol{\theta} \leftarrow \boldsymbol{\theta} - \frac{\eta}{|\mathcal{B}|} \sum_{i \in \mathcal{B}} \nabla_{\boldsymbol{\theta}} \ell^{(i)}(\boldsymbol{\theta})\]
其中梯度是損失有關3個為標量的模型引數的偏導陣列成的向量：
\[\nabla_{\boldsymbol{\theta}} \ell^{(i)}(\boldsymbol{\theta})= \begin{bmatrix} \frac{ \partial \ell^{(i)}(w_1, w_2, b) }{\partial w_1} \\ \frac{ \partial \ell^{(i)}(w_1, w_2, b) }{\partial w_2} \\ \frac{ \partial \ell^{(i)}(w_1, w_2, b) }{\partial b} \end{bmatrix} = \begin{bmatrix} x_1^{(i)} (x_1^{(i)} w_1 + x_2^{(i)} w_2 + b - y^{(i)}) \\ x_2^{(i)} (x_1^{(i)} w_1 + x_2^{(i)} w_2 + b - y^{(i)}) \\ x_1^{(i)} w_1 + x_2^{(i)} w_2 + b - y^{(i)} \end{bmatrix} = \begin{bmatrix} x_1^{(i)} \\ x_2^{(i)} \\ 1 \end{bmatrix} (\hat{y}^{(i)} - y^{(i)})\]

從零實現

這裡不利用高階框架，只用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的全部程式碼。
線性迴歸簡潔實現.