线性回归原理及实现

1. 公式推导

原式

$$\begin{array}{rl}J\left(\theta \right)& =\frac{1}{2}{\left(X\theta -Y\right)}^T\left(X\theta -Y\right)\\ & \\ & =\frac{1}{2}\left({\theta }^T{X}^T-Y^T\right)\left(X\theta -Y\right)\\ & \\ & =\frac{1}{2}\left({\theta }^T{X}^{T}X\theta -{\theta }^T{X}^{T}Y-Y^{T}X\theta +Y^{T}Y\right)\end{array}$$

求导

求导知识
$$\frac{d{X}^{T}AX}{dX}=2AX$$

$$\frac{d{X}^{T}A}{dX}=A$$

$$\frac{dAX}{dX}=A^T$$

---

开始求导

$$\begin{array}{rl}\frac{dJ(\theta )}{d\theta }& =\frac{1}{2}(\frac{d{\theta }^T{X}^{T}X\theta }{d\theta }-\frac{d{\theta }^T{X}^{T}Y}{d\theta }-\frac{d{Y}^{T}X\theta }{d\theta }+\frac{d{Y}^{T}Y}{d\theta })\\ & \\ & =\frac{1}{2}(2X^{T}X\theta -X^{T}Y-X^{T}Y)\\ & \\ & =X^{T}X\theta -X^{T}Y\end{array}$$

令导数为0

$$X^{T}X\theta -X^{T}Y=0$$

$$\theta =(X^{T}X{)}^{-1}{X}^{T}Y$$

2. torch实现

import torch
from torch import nn


class LinearReg(nn.Module):
    def __init__(self, input_dim, out_dim):
        super(LinearReg, self).__init__()
        self.input_dim = input_dim
        self.out_dim = out_dim
        self.linear = nn.Linear(input_dim, out_dim, bias=True)

    def forward(self, x):
        out = self.linear(x)
        return out


def main():
    data = torch.randn(100, 2)*10
    # 假设 y=3*x_1 + 2*x_2 + 5
    weight = torch.tensor([[3.], [2]])
    bias = torch.tensor([[5.]])
    # 构建数据集,增加一些扰动
    y = (data @ weight + bias) + torch.randn(100, 2)*2

    model = LinearReg(input_dim=2, out_dim=1)
    loss_func = nn.MSELoss(reduction="mean")
    optimizer = torch.optim.SGD(lr=5e-3, params=model.parameters())
    epochs = 1000
    for step in range(epochs):
        pred = model(data)
        loss = loss_func(pred, y)
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()

        if (step+1) % 10 == 0:
            print(f"{step}/{epochs} steps, loss: {loss.item():.4f}")
    print("train finished")

    # 打印模型权重
    print("训练后模型权重如下")
    print("weight", model.linear.weight)
    print("bias", model.linear.bias)


if __name__ == "__main__":
    main()

输出如下

手动求梯度实现

import torch


def diff_mse(x, y, w):
    """
    delta_w = x.t@(x@w-y)
    """
    return x.transpose(0, 1)@(x@w-y)/x.shape[0]


def mse_loss(x, y, w):
    return 0.5*torch.mean(torch.square(x@w-y))


def get_batch_data(x, y, batch_size, step):
    data_len = x.shape[0]
    start = step*batch_size
    end = min(start + batch_size, data_len)
    return x[start:end], y[start:end]


def train(epochs, batch_size, lr):
    data = torch.randn(100, 2)*2    # [100, 2]
    weight = torch.tensor([[3.], [2]])   # [2, 1]
    y = data@weight + torch.randn(100, 1)*2

    param_w = torch.randn(2, 1)
    steps = data.shape[0]//batch_size
    for epoch in range(epochs):
        for step in range(steps):
            x, lb = get_batch_data(data, y, batch_size, step)
            loss = mse_loss(x, lb, param_w)
            grad = diff_mse(x, lb, param_w)
            param_w = param_w - lr*grad
            if step % 10 == 0:
                print(f"epoch:{epoch}; step:{step}; loss:{loss.item()}")
    print(f"train finished, param w: {param_w}")


if __name__ == "__main__":
    train(epochs=200, batch_size=8, lr=5e-4)