Pytorch学习2020春-1-线性回归

提示：扫一扫查出行
【扫一扫了解最新限行尾号】
复制提示

线性回归

主要内容包括：

1. 线性回归的基本要素
2. 线性回归模型从零开始的实现
3. 线性回归模型使用pytorch的简洁实现

线性回归的基本要素

模型

为了简单起见，这里我们假设价格只取决于房屋状况的两个因素，即面积（平方米）和房龄（年）。接下来我们希望探索价格与这两个因素的具体关系。线性回归假设输出与各个输入之间是线性关系:

数据集

我们通常收集一系列的真实数据，例如多栋房屋的真实售出价格和它们对应的面积和房龄。我们希望在这个数据上面寻找模型参数来使模型的预测价格与真实价格的误差最小。在机器学习术语里，该数据集被称为训练数据集（training data set）或训练集（training set），一栋房屋被称为一个样本（sample），其真实售出价格叫作标签（label），用来预测标签的两个因素叫作特征（feature）。特征用来表征样本的特点。

损失函数

在模型训练中，我们需要衡量价格预测值与真实值之间的误差。通常我们会选取一个非负数作为误差，且数值越小表示误差越小。一个常用的选择是平方函数。 它在评估索引为 的样本误差的表达式为

此处的\(\mathbf{w}^\top\mathbf{x}^{(i)}\)是两个矩阵相乘，向量点积可以用行矩阵乘以列矩阵来实现

优化函数 - 随机梯度下降

当模型和损失函数形式较为简单时，上面的误差最小化问题的解可以直接用公式表达出来。这类解叫作解析解（analytical solution）。本节使用的线性回归和平方误差刚好属于这个范畴。然而，大多数深度学习模型并没有解析解，只能通过优化算法有限次迭代模型参数来尽可能降低损失函数的值。这类解叫作数值解（numerical solution）。

插曲：什么是梯度下降？

大多数机器学习或者深度学习算法都涉及某种形式的优化。 优化指的是改变 （此处x是一个向量）以最小化或最大化某个函数 的任务。 我们通常以最小化 指代大多数最优化问题。 最大化可经由最小化算法最小化 来实现。

我们把要最小化或最大化的函数称为目标函数或准则。 当我们对其进行最小化时，我们也把它称为代价函数、损失函数或误差函数。

下面，我们假设一个损失函数为 ，其中 然后要使得最小化它。

注意：这里只是假设，不用知道这个目标函数就是平方损失函数等等，然后肯定有人问既然要最小化它，那求个导数，然后使得导数等于0求出不就好了吗？Emmmm...是的，有这样的解法，可以去了解正规方程组求解。说下这里不讲的原因，主要是那样的方式太难求解，然后在高维的时候，可能不可解，但机器学习或深度学习中，很多都是超高维的，所以也一般不用那种方法。总之，梯度下降是另一种优化的不错方式，比直接求导好很多。

什么是梯度？

梯度概念是建立在偏导数与方向导数概念基础上的。所谓偏导数，简单来说是对于一个多元函数，选定一个自变量并让其他自变量保持不变，只考察因变量与选定自变量的变化关系。数学上说，是指对于多元函数\(y=f(x_1,x_2,...x_n)\)，假设其偏导数都存在，则该函数共有n个偏导数，可以表示为：

偏导数只能表示多元函数沿某个坐标轴方向的导数，如对于二元函数\(z=x^2+y^2\;,\;{{\partial z} \over {\partial x}} = 2x\)表示函数沿X轴方向的导数，而\({{\partial z} \over {\partial y}} = 2y\)表示函数沿Y轴方向的导数。

除开沿坐标轴方向上的导数，多元函数在非坐标轴方向上也可以求导数，这种导数称为方向导数。很容易发现，多元函数在特定点的方向导数有无穷多个，表示函数值在各个方向上的增长速度。一个很自然的问题是：在这些方向导数中，是否存在一个最大的方向导数，如果有，其值是否唯一？为了回答这个问题，便需要引入梯度的概念。

一般来说，梯度可以定义为一个函数的全部偏导数构成的向量(这一点与偏导数与方向导数不同，两者都为标量)。一般将函数f的梯度记为\(∇f\)

定理：如果函数z=f(x,y)在点p(x,y)是可微分的，函数的增量可以表示为

两边同除以Rho可以得到

同时取rho趋于零的极限

可知，这个关于theta的函数最大值在特殊点取到，这时方向正好是x方向梯度和y方向梯度向量的和向量方向

梯度下降：我们知道曲面上方向导数的最大值的方向就代表了梯度的方向，因此我们在做梯度下降的时候，应该是沿着梯度的反方向进行权重的更新，可以有效的找到全局的最优解。这个 的更新过程可以描述为

续：

在求数值解的优化算法中，小批量随机梯度下降（mini-batch stochastic gradient descent）在深度学习中被广泛使用。它的算法很简单：先选取一组模型参数的初始值，如随机选取；接下来对参数进行多次迭代，使每次迭代都可能降低损失函数的值。在每次迭代中，先随机均匀采样一个由固定数目训练数据样本所组成的小批量（mini-batch），然后求小批量中数据样本的平均损失有关模型参数的导数（梯度），最后用此结果与预先设定的一个正数的乘积作为模型参数在本次迭代的减小量。

学习率: η代表在每次优化中，能够学习的步长的大小

批量大小: B是小批量计算中的批量大小batch size

总结一下，优化函数的有以下两个步骤：

• (i)初始化模型参数，一般来说使用随机初始化；
• (ii)我们在数据上迭代多次，通过在负梯度方向移动参数来更新每个参数。

1.基本设置

题外话：这次安装torch的过程很痛苦，失败很多次，只是因为没有更换镜像源。使用一下的pip语句（通过阿里镜像源）安装可以避免从国外镜像源下载过慢而出错的结果,一般在安装出错的时候都要多考虑一下是不是因为网络的问题，即时更换镜像,记得把下面那个xxx改成你要安装的镜像名

pip3 install xxx -i http://mirrors.aliyun.com/pypi/simple/ --trusted-host mirrors.aliyun.com


矢量运算用torch集成的函数会更快一些

import torch
import time# init variable a, b as 1000 dimension vector
n = 1000
a = torch.ones(n)
b = torch.ones(n)


# define a timer class to record time
class Timer(object):
"""Record multiple running times."""
def __init__(self):
self.times = []
self.start() def start(self):
# start the timer
self.start_time = time.time() def stop(self):
# stop the timer and record time into a list
self.times.append(time.time() - self.start_time)
return self.times[-1] def avg(self):
# calculate the average and return
return sum(self.times)/len(self.times) def sum(self):
# return the sum of recorded time
return sum(self.times)


现在我们可以来测试了。首先将两个向量使用for循环按元素逐一做标量加法。

timer = Timer()
c = torch.zeros(n)
for i in range(n):
c[i] = a[i] + b[i]
'%.5f sec' % timer.stop()


另外是使用torch来将两个向量直接做矢量加法：

timer.start()
d = a + b
'%.5f sec' % timer.stop()


结果很明显,后者比前者运算速度更快。因此，我们应该尽可能采用矢量计算，以提升计算效率。

2.线性回归模型从零开始的实现

# import packages and modules
%matplotlib inline
import torch
from IPython import display
from matplotlib import pyplot as plt
import numpy as np
import randomprint(torch.__version__)


2.1 生成数据集

我们构造一个简单的人工训练数据集，它可以使我们能够直观比较学到的参数和真实的模型参数的区别。设训练数据集样本数为1000，输入个数（特征数）为2。给定随机生成的批量样本特征\(\boldsymbol{X} \in \mathbb{R}^{1000 \times 2}\)，我们使用线性回归模型真实权重\(\boldsymbol{w} = [2, -3.4]^\top\)和偏差\(b = 4.2\)，以及一个随机噪声项\(\epsilon\)来生成标签

其中噪声项\(\epsilon\)服从均值为0、标准差为0.01的正态分布。噪声代表了数据集中无意义的干扰。

# 设置 feature（自变量）的数量
num_inputs = 2
# 设置样本（可探测出的已知的自变量和因变量值）的数量
num_examples = 1000# 在这里为了生存伪实验数据，先定好参数（weight）的值
true_w = [2, -3.4]
true_b = 4.2# torch.randn()，正态分布
features = torch.randn(num_examples, num_inputs,
dtype=torch.float32)
labels = true_w[0] * features[:, 0] + true_w[1] * features[:, 1] + true_b
labels += torch.tensor(np.random.normal(0, 0.01, size=labels.size()),
dtype=torch.float32)


randn()的定义如下

def randn(size: _int, generator: Generator, names: Optional[List[Union[str, None]]],

out: Optional[Tensor]=None, dtype: _dtype=None, layout: _layout=strided,

device: Union[_device, str, None]=None, requires_grad:_bool=False) -> Tensor: ...

可以知道这个randn是用来生成随机数的

Returns a tensor filled with random numbers from a normal distribution

with mean 0 and variance 1 (also called the standard normal

distribution).

此处的size是一个可以读入列表的参数，如果输入一些逗号隔开的整数，*号给函数一个讯息，让它把整数整合为一个列表，如果输入两个，第一个就是行数，第二个是列数

lables的初始设置利用了python的广播机制

列向量，行向量，标量，会自动复制，变成和最大的矩阵相同的size，然后相加减

此处调用tensor[:k]就是调用第k列向量，展现出了的形式就是一个列向量，那么被赋值的就是一个列向量，因此labels就是我们的因变量的合集

最后给labels加上一些偏差，torch.tensor()是给定已知数据然后帮你集成出一个tensor的函数，数据生成则需要在其参数内部另外嵌入其他函数来生成，最好生成个列表，labels.size()会给出一个一维列表，分别告知numpy.random.normal()函数目标tensor的每个维度是多大

2.2 生成散点图

通过生成第二个特征features[:, 1]和标签 labels 的散点图，可以更直观地观察两者间的线性关系。

def use_svg_display():
# 用矢量图显示
display.set_matplotlib_formats('svg')def set_figsize(figsize=(3.5, 2.5)):
use_svg_display()
# 设置图的尺寸
plt.rcParams['figure.figsize'] = figsizeset_figsize()
plt.scatter(features[:, 1].numpy(), labels.numpy(), 1); # 加分号只显示图


我们将上面的plt作图函数以及use_svg_display函数和set_figsize函数定义在d2lzh包里。以后在作图时，我们将直接调用d2lzh.plt。由于plt在d2lzh包中是一个全局变量，我们在作图前只需要调用d2lzh.set_figsize()即可打印矢量图并设置图的尺寸。

plt.scatter()必须封入数组，所以只用.numpy()的attribute。将该tensor以NumPy的形式返回ndarray，两者共享相同的底层内存。原tensor改变后会相应的在ndarray有反映，反之也一样。

2.3 读取数据

def data_iter(batch_size, features, labels):
num_examples = len(features)
indices = list(range(num_examples))
random.shuffle(indices) # random read 10 samples
for i in range(0, num_examples, batch_size):
j = torch.LongTensor(indices[i: min(i + batch_size, num_examples)]) # the last time may be not enough for a whole batch
yield features.index_select(0, j), labels.index_select(0, j)


代码是这样运行的：先生成一个index，0~1000-1的整数顺序列表，用random.shuffle(*list)的方式打乱（洗牌），接着用循环，十个十个跳着，把每10个单位长度的区间作为一个Batch，用yield生成占用最小内存的表格式返回给使用这个迭代函数的流程

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


2.4 初始化模型参数

w = torch.tensor(np.random.normal(0, 0.01, (num_inputs, 1)), dtype=torch.float32)
b = torch.zeros(1, dtype=torch.float32)w.requires_grad_(requires_grad=True)
b.requires_grad_(requires_grad=True)


这里，注意区别tensor.requires_grad和tensor.requires_grad_()两个东西，前面是调用变量的属性值，后者是调用内置的函数，来改变属性。

2.5 求导

Reference: https://www.cnblogs.com/timlong/p/11742535.html
https://www.cnblogs.com/timlong/p/11746848.html
https://en.wikipedia.org/wiki/Automatic_differentiation


2.5.1 手动微分法

研究者针对问题的模型自己求解出梯度表达式

2.5.2 数值微分法

根据梯度的定义式 $$\lim\limits_{\rho\rightarrow 0}\frac{f(x+\Delta x,y+\Delta y)-f(x,y)}{\rho}$$ 取一个很小的\(\rho\)值，计算出梯度的近似解

（其实和下面自动求导的二元数原理没什么差别）

2.5.3 符号微分法

符号微分技术在数学软件如MATLAB、Maple及Mathematica等均有应用。其原理是基于下面的求导法则

例如要求符号表达式\(f(x)=2x+x^2\)的导数

就需要先将其展开为一棵符号树（利用我Vincent项目里面的iGraph module绘制的）

基于一系列运算就可得到\(f'(x)=2+2x\)

有一点要注意的是，符号微分不一定会得到简化的导数，因为计算机可能并不能进行智能的简化。所以，如果表达式树结构较复杂时，得到的导数表达式会相当复杂，也许出现表达式爆炸现象。

2.5.4 自动求导法

感觉自动求导本质上是一种优化了的数值方法，就是用了一种专门的算法来解决问题，经过优化后能够handle非常复杂的多元函数

最基本的求导还是数值方法，是基于二元数（一个实数加上一个实数与小量的乘积）\(a+b\epsilon\)

\(a,b\in\mathbb{R}\),\(\epsilon\)是无穷小量，满足\(\epsilon^2=0\)

二元数的基本性质有

回忆一下微分的定义，我们知道

类似的，对于二元数

因此通过这一的求值可以计算出导数值

而自动求导和粗暴数值解的区别在于它结合了符号和数值方法的精华，通过多步计算和算法优化，能高效且误差小地算出导数

它的基础是链式法则，用一些新变量整体替换

把这些处理出的部分的、整体的变量当做节点，根据运算逻辑处理出一张DAG（有向无环图）

有两种模式

1. forward accumulation computes the recursive relation \(\displaystyle {\frac {dw_{i}}{dx}}={\frac {dw_{i}}{dw_{i-1}}}{\frac {dw_{i-1}}{dx}}\) with \(\displaystyle w_{3}=y\)
2. reverse accumulation computes the recursive relation \(\displaystyle {\frac {dy}{dw_{i}}}={\frac {dy}{dw_{i+1}}}{\frac {dw_{i+1}}{dw_{i}}}\) with \(\displaystyle w_{0}=x\)

2.5.3.1 Forward Mode（前向模式）

$\displaystyle {\frac {\partial y}{\partial x}}={\frac {\partial y}{\partial w_{n-1}}}{\frac {\partial w_{n-1}}{\partial x}}={\frac {\partial y}{\partial w_{n-1}}}\left({\frac {\partial w_{n-1}}{\partial w_{n-2}}}{\frac {\partial w_{n-2}}{\partial x}}\right)={\frac {\partial y}{\partial w_{n-1}}}\left({\frac {\partial w_{n-1}}{\partial w_{n-2}}}\left({\frac {\partial w_{n-2}}{\partial w_{n-3}}}{\frac {\partial w_{n-3}}{\partial x}}\right)\right)=\cdots $

举个例子（来自wikipedia）

先求出w1和w2的导数

接着一步步地向上走，如图中箭头所示，最后得到

“Forward accumulation is more efficient than reverse accumulation for functions f : ℝn → ℝm with m ≫ n as only n sweeps are necessary, compared to m sweeps for reverse accumulation.”

个人认为前向模式是一种便于人理解的方法，另一种方式反馈模式才是计算机能够操作的，把前向模式转化为一个递归流程

2.5.3.2 Reverse Mode（反馈模式）

$\displaystyle {\frac {\partial y}{\partial x}}={\frac {\partial y}{\partial w_{1}}}{\frac {\partial w_{1}}{\partial x}}=\left({\frac {\partial y}{\partial w_{2}}}{\frac {\partial w_{2}}{\partial w_{1}}}\right){\frac {\partial w_{1}}{\partial x}}=\left(\left({\frac {\partial y}{\partial w_{3}}}{\frac {\partial w_{3}}{\partial w_{2}}}\right){\frac {\partial w_{2}}{\partial w_{1}}}\right){\frac {\partial w_{1}}{\partial x}}=\cdots $

2.6 定义模型

def linreg(X, w, b):
return torch.mm(X, w) + b


2.7 定义损失函数

单元素损失原始公式

批量损失公式

其中\(|\mathcal{B}|\)是批量取样的大小

def squared_loss(y_hat, y):
return (y_hat - y.view(y_hat.size())) ** 2 / 2


2.8 定义优化函数

此处使用小批量随机下降

def sgd(params, lr, batch_size):
for param in params:
param.data -= lr * param.grad / batch_size # ues .data to operate param without gradient track


此处的\(\eta\)是学习率，剩下的部分是对批量损失函数求导的结果表达式

2.9 训练

当数据集、模型、损失函数和优化函数定义完了之后就可来准备进行模型的训练了。

# super parameters init
lr = 0.03
num_epochs = 5net = linreg
loss = squared_loss# training
for epoch in range(num_epochs): # training repeats num_epochs times
# in each epoch, all the samples in dataset will be used once # X is the feature and y is the label of a batch sample
for X, y in data_iter(batch_size, features, labels):
l = loss(net(X, w, b), y).sum()
# calculate the gradient of batch sample loss
l.backward()
# using small batch random gradient descent to iter model parameters
sgd([w, b], lr, batch_size)
# reset parameter gradient
w.grad.data.zero_()
b.grad.data.zero_()
train_l = loss(net(features, w, b), labels)
print('epoch %d, loss %f' % (epoch + 1, train_l.mean().item()))


w, true_w, b, true_b


2.9.1 Tensor的AutoGrad

Reference: https://www.jianshu.com/p/a105858567df

y是一个tensor，那么(y ** 2)/2就是y的对应项相乘再每个项都除以2形成的新tensor，tensor.sum()方法是将tensor的每个项相加，形成一个只有一个元素的tensor，可以看做标量(scalar)

经过之前的语句的控制，定义的tensor和经过运算得到的衍生的tensor都具有自动求导性质，每次运算的路径会被记录下来，形成一个计算图，在调用自动求导的时候，会利用计算图计算出梯度

被调用过的计算图部分默认被销毁，那么计算图就只剩下z、q了，已经不完整，无法求导了。

所以这个时候，无论你是想再次运行p.backward()还是q.backward()，都无法进行，报错如下：

RuntimeError: Trying to backward through the graph a second time, but the buffers have already been freed. Specify retain_graph=True when calling backward the first time.

此时有两种解决办法

1.更改你的backward函数，添加参数retain_graph=True，重新进行backward，这个时候你的计算图就被保留了，不会报错。

但是这样会吃内存！，尤其是，你在大量迭代进行参数更新的时候，很快就会内存不足，memory out了。

2.你实际根本没必要对一个计算图backward多次，而你不小心多跑了一次backward函数。

通常，你要是在IPython里面联系PyTorch的时候，因为你会反复运行一个单元格的代码，所以很容易一不小心把backward运行了多次，就会报错。这个时候，你就检查一下代码，防止backward运行多次即可。

3. 线性回归用pytorch简洁实现

import torch
from torch import nn
import numpy as np
torch.manual_seed(1)print(torch.__version__)
torch.set_default_tensor_type('torch.FloatTensor')


1.4.0


# 生成数据集，这里和从零实现的代码是一样的
num_inputs = 2
num_examples = 1000true_w = [2, -3.4]
true_b = 4.2features = torch.tensor(np.random.normal(0, 1, (num_examples, num_inputs)), dtype=torch.float)
labels = true_w[0] * features[:, 0] + true_w[1] * features[:, 1] + true_b
labels += torch.tensor(np.random.normal(0, 0.01, size=labels.size()), dtype=torch.float)


# 读取数据集
import torch.utils.data as Databatch_size = 10# combine featues and labels of dataset
dataset = Data.TensorDataset(features, labels)# put dataset into DataLoader
data_iter = Data.DataLoader(
dataset=dataset, # torch TensorDataset format
batch_size=batch_size, # mini batch size
shuffle=True, # whether shuffle the data or not
num_workers=2, # read data in multithreading
)


# 定义模型
class LinearNet(nn.Module):
def __init__(self, n_feature):
super(LinearNet, self).__init__() # call father function to init
self.linear = nn.Linear(n_feature, 1) # function prototype: `torch.nn.Linear(in_features, out_features, bias=True)` def forward(self, x):
y = self.linear(x)
return ynet = LinearNet(num_inputs)
print(net)


LinearNet(
(linear): Linear(in_features=2, out_features=1, bias=True)
)


# 生成多层网络
# ways to init a multilayer network
# method one
net = nn.Sequential(
nn.Linear(num_inputs, 1)
# other layers can be added here
)# method two
net = nn.Sequential()
net.add_module('linear', nn.Linear(num_inputs, 1))
# net.add_module ......# method three
from collections import OrderedDict
net = nn.Sequential(OrderedDict([
('linear', nn.Linear(num_inputs, 1))
# ......
]))print(net)
print(net[0])


# 初始化模型参数
from torch.nn import initinit.normal_(net[0].weight, mean=0.0, std=0.01)
init.constant_(net[0].bias, val=0.0) # or you can use `net[0].bias.data.fill_(0)` to modify it directly


# 定义损失函数
loss = nn.MSELoss() # nn built-in squared loss function
# function prototype: `torch.nn.MSELoss(size_average=None, reduce=None, reduction='mean')`


# 定义优化函数
import torch.optim as optimoptimizer = optim.SGD(net.parameters(), lr=0.03) # built-in random gradient descent function
print(optimizer) # function prototype: `torch.optim.SGD(params, lr=, momentum=0, dampening=0, weight_decay=0, nesterov=False)`


# 训练
num_epochs = 3
for epoch in range(1, num_epochs + 1):
for X, y in data_iter:
output = net(X)
l = loss(output, y.view(-1, 1)) # 讲输出的tensor从左到右，从上到下顺序转化成一个列向量
optimizer.zero_grad() # reset gradient, equal to net.zero_grad()
l.backward()
optimizer.step()
print('epoch %d, loss: %f' % (epoch, l.item()))


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