线性回归(Linear Regression)原理及手工实现_解析解法、梯度下降法求解最优解

一、什么是线性回归

直观来讲，就是$y=kx+b$的形式，可以用一条直线来进行拟合。
一元线性回归：只有一个自变量和一个因变量，$y=k_1{x}_1+b$的形式表示。
多元线性回归：两个及以上的自变量和一个因变量，$y=k_1{x}_1+k_2{x}_2+...+k_n{x}_n+b$的形式表示。
在这里可以再与多项式回归进行区分，$y=k_1{x}_1+k_2{x}_2^2+...+k_n{x}_n^n+b$的形式表示。

二、最小二乘法与解析解

“基于均方误差（mean-square error, MSE）最小化来进行模型求解的方法称为“最小二乘法””–《机器学习》周志华
最小二乘法(又称最小平方法)可采用均方误差来定义损失函数：$$L(w)=\sum _{i=1}^n(w^{\mathrm{T}}{x}_i-y_i{)}^2$$若要最小化损失函数，即可转化为对其求导，令导数为0，即可得到$w$的解析解。
因此可以定义目标函数$J(w)=argminL(w)$,找到一个$w$,使得损失函数$L(w)$的值最小。
最终得到的解析解：$$\hat{w}=(X^{\mathrm{T}}X{)}^{-1}{X}^{\mathrm{T}}y$$详细推导过程可参考:https://www.bilibili.com/video/BV1hW41167iL
“像线性回归这样的简单问题存在解析解，但并不是所有的问题都存在解析解。解析解可以进行很好的数学分析，但解析解的限制很严格，导致它无法应用在深度学习里。” --https://www.bilibili.com/video/BV1hW41167iL
因此引入梯度下降(gradient descent)方法,是深度学习模型中常见的优化方法，可以用于求解最小二乘问题(线性非线性都可以)。

三、梯度下降法

梯度下降法的计算过程就是沿梯度下降的方向求解极小值,进行多次迭代:$$w=w-\eta \frac{\partial }{\partial w}J(w)$$其中$w$为模型的可学习参数，$J(w)$为目标函数，$\eta$为学习率。
批量梯度下降法(Batch Gradient Descent, BGD)：使用整个训练集去计算目标函数的梯度。由于每次使用全部数据来计算梯度去更新参数，速度会很慢，而且数据很难一次性全部载入内存当中进行计算。首先定义单个样本的损失函数：$$l^{(i)}(w)=\frac{1}{2}(w^{\mathrm{T}}{x}^{(i)}-y^{(i)}{)}^2$$采用全部训练集数据更新的批量梯度下降法如下所示，其中$N$为训练集样本的总个数，$j$为迭代的次数：$$w_{j+1}=w_j-\frac{\eta }{N}\sum _{i=1}^N\frac{\partial }{\partial w_j}{l}^{(i)}(w_j)$$ 随机梯度下降法(Stochastic Gradient Descent, SGD)：每次只采用一个样本来计算梯度。速度快但是每次迭代方向变化大，导致不能快速收敛到最优解。$$w_{j+1}=w_j-\eta \frac{\partial }{\partial w_j}{l}^{(i)}(w_j)$$小批量梯度下降法(Mini-Batch Gradient Descent, MBGD)：对批量梯度下降和随机梯度下降进行了折中，随机抽取固定数量的训练样本$\beta$$(1<|\beta |<N)$，与批量梯度下降法不同的是将$N$改为了$|\beta |$：$$w_{j+1}=w_j-\frac{\eta }{|\beta |}\sum _{i=1}^{|\beta |}\frac{\partial }{\partial w_j}{l}^{(i)}(w_j)$$

四、线性回归从零实现

假设我们有1000个样本的训练集，每个样本$x^{(i)}$是3维的，也就是有三个特征，$x_1^{(i)}$，$x_2^{(i)}$，$x_3^{(i)}$。有要拟合的直线:$$Y=w_1{x}_1+w_2{x}_2+w_3{x}_3+b=W^{\mathrm{T}}X+b$$假设$Y=2x_1-3x_2+4x_3-5$，如何使用线性回归模型和梯度下降法来不断迭代求得参数呢？

1、首先生成数据集

对于$Y=2x_1-3x_2+4x_3-5$，训练集中$X\in R^{1000\times 3}$，$Y\in R^{1000\times 1}$。

import torch
def generate_data(w, b, num_examples):  
    """ Y = 2a - 3b + 4c - 5 + noise """
    X = torch.normal(0, 1, (num_examples, len(w)))
    y = torch.matmul(X, w) + b
    y += torch.normal(0, 0.01, y.shape)
    return X, y.reshape((-1, 1))

true_w = torch.tensor([2., -3., 4.])
true_b = 5.
features, labels = generate_data(true_w, true_b, 1000)

2、小批量读取数据集

1000个数据中，batch_size=8，整除后共125组。batch_X[125,8,3]，batch_y[125,8,1]。
这里先将1000个样本数据打乱，然后按照顺序、无放回地8个8个取样本组成一个batch。因此，在后面训练模型的时候，每一个epoch，都会按组(125组)来遍历，每次会遍历完所有组。

def data_iterater(X, y, batch_size):
    num = len(X)
    indices = list(range(num))
    random.shuffle(indices)  # 将顺序打乱
    batch_X = torch.zeros([num//batch_size, batch_size, len(true_w)])
    batch_y = torch.zeros([num//batch_size, batch_size, 1])
    for id, i in enumerate(range(0, num, batch_size)):
        batch_indices = torch.tensor(indices[i: min(i + batch_size, num)])
        batch_X[id,:,:] = features[batch_indices]
        batch_y[id,:,:] = labels[batch_indices]
    return batch_X, batch_y

batch_size = 8
batch_X, batch_y = data_iterater(features, labels, batch_size)

3、定义模型

# 定义线性模型
def Linear_Model(X, w, b):
    return torch.matmul(X, w) + b

4、定义损失函数、优化方法

# 定义均方损失函数
def loss_sq(predict, y):
    return (predict - y) ** 2 / 2
    
# 定义小批量梯度下降优化方法
def sgd(params, lr, batch_size):
    with torch.no_grad():
        for param in params:
            param -= lr * param.grad / batch_size
            param.grad.zero_()

5、训练

训练前初始化模型参数$w,b$。以及相应超参数学习率$lr$，模型迭代次数$epochs$等等。

w = torch.normal(0, 0.01, size=(3,1), requires_grad=True)
b = torch.zeros(1, requires_grad=True)
lr = 1e-3
epochs = 100
for epoch in range(epochs):
    """在每个epoch中，小批量梯度遍历整个数据集，也就是每次batch_size大小取样本
    并将训练数据集中所有样本都使用一次（假设样本数能够被批量大小整除）。"""
    for (X, y) in zip(batch_X, batch_y):
        predict = Linear_Model(X, w, b)
        loss = loss_sq(predict, y)
        loss.sum().backward()  # 计算梯度
        sgd([w, b], lr, batch_size)  # 更新参数值
    with torch.no_grad():
        train_l = loss_sq(Linear_Model(features, w, b), labels)
        print(f'epoch {epoch + 1}, loss {float(train_l.mean()):f}')
print(f'true_w {true_w}, true_b {true_b}')
print(f'pred_w {w}, pred_b {b}')

6、结果

epoch 1, loss 21.021652
epoch 2, loss 16.366686
epoch 3, loss 12.744938
epoch 4, loss 9.926495
epoch 5, loss 7.732749
.......................
epoch 99, loss 0.000052
epoch 100, loss 0.000052
true_w tensor([ 2., -3.,  4.]), true_b 5.0
pred_w tensor([[ 2.0000],
        [-2.9998],
        [ 3.9995]], requires_grad=True), pred_b tensor([5.0002], requires_grad=True)

五、参考

[1]https://zh-v2.d2l.ai/chapter_linear-networks/linear-regression-scratch.html
[2]https://www.bilibili.com/video/BV1hW41167iL