PyTorch执行矩阵求导详细过程

本文用最小二乘法来简述PyTorch中的梯度下降具体执行过程

之前已经推到过最小二乘法的矩阵求导过程：链接

一、数学形式

之前虽然知道矩阵求导是怎么来的，但是具体PyTorch是怎么实现的还未知，所以这里进行了求解和验证，采用了SGD梯度下降，为方便计算，batch_size=3，并且采用一元线性回归来实现，所以这里w只有一维，带有bias偏置

数学表达式为
$$Y=Xw+b$$

$$Y={\left(\begin{array}{c}{y}_1\\ y_2\\ ...\\ y_n\end{array}\right)}_{n\times 1}X={\left(\begin{array}{c}{x}_1^T\\ x_2^T\\ ...\\ x_n^T\end{array}\right)}_{n\times p}w={\left(\begin{array}{c}{w}_1\\ w_2\\ ...\\ w_n\end{array}\right)}_{p\times 1}$$

注意为了方便下面验证，取一元线性回归，即参数
$$Y={\left(\begin{array}{c}{y}_1\\ y_2\\ ...\\ y_n\end{array}\right)}_{n\times 1}X={\left(\begin{array}{c}{x}_1^T\\ x_2^T\\ ...\\ x_n^T\end{array}\right)}_{n\times 1}\to {\left(\begin{array}{c}{x}_1\\ x_2\\ ...\\ x_n\end{array}\right)}_{n\times 1}w={\left(\begin{array}{c}w\end{array}\right)}_{1\times 1}$$
本来$x^T$为$1\times p$的矩阵，即n元，现在改为一元，转置相当于没变化了

接下来给出矩阵求导以及非矩阵求导的公式，这里我们取batch_size=3，即$n=3$

损失函数为：
$$L=\frac{1}{n}\sum _{i=1}^n(y_i-x_i^{T}w{)}^2$$

1. 矩阵求导

$$\frac{dL}{dw}=2\frac{1}{n}{X}^T(Xw-Y)$$

矩阵求导的时候没有算上偏置$b$，下面非矩阵求导会算上，而且主要是用非矩阵求导的过程来验证PyTorch运行过程，所以不影响

将$n=3$带入得
$$\begin{array}{rl}\frac{dL}{dw}& =2\frac{1}{n}{X}^T(Xw-Y)\\ & =\frac{2}{3}\left(\begin{array}{ccc}{x}_1& x_2& x_3\end{array}\right)(\left(\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right)\left(\begin{array}{c}w\end{array}\right)-\left(\begin{array}{c}{y}_1\\ y_2\\ y_3\end{array}\right))\\ & =\frac{2}{3}\left(\begin{array}{ccc}{x}_1& x_2& x_3\end{array}\right)\left(\begin{array}{c}{x}_{1}w-y_1\\ x_{2}w-y_2\\ x_{3}w-y_3\end{array}\right)\\ & =\frac{2}{3}\left(\begin{array}{c}{x}_1^{2}w-x_1{y}_1\\ x_2^{2}w-x_2{y}_2\\ x_3^{2}w-x_3{y}_3\end{array}\right)\end{array}$$

2. 非矩阵求导

$$\frac{dL}{dw}=2\frac{1}{n}(w\sum _{i=1}^n{x}_i^2-\sum _{i=1}^n(y_i-b)x_i)$$

将$n=3$带入得
$$\begin{array}{rl}\frac{dL}{dw}& =2\frac{1}{n}(w\sum _{i=1}^n{x}_i^2-\sum _{i=1}^n(y_i-b)x_i)\\ & =\frac{2}{3}(w(x_1^2+x_2^2+x_3^2)-[((y_1-b)x_1)+((y_2-b)x_2)+((y_3-b)x_3)])\\ & =\frac{2}{3}(w(x_1^2+x_2^2+x_3^2)-[(x_1{y}_1-x_{1}b)+(x_2{y}_2-x_{2}b)+(x_3{y}_3-x_{3}b)])\\ & =\frac{2}{3}(w(x_1^2+x_2^2+x_3^2)-(x_1{y}_1+x_2{y}_2+x_3{y}_3)+(x_{1}b+x_{2}b+x_{3}b))\end{array}$$
上述就是对$w$求导的结果，由于$b$相同，直接给出结果
$$\frac{dL}{db}=\frac{2}{n}(nb-\sum _{i=1}^n(y_i-w{x}_i))$$

下面验证时仅针对有偏置$b$的一元线性回归，并且只验证$w$

二、PyTorch

这里我们batch_size=3，随机取输入$x$和$y$分别为

输出为

loss为

$w$的梯度为-150.28

更新前后$w$与$b$的参数为

• $w:-0.34487\to -0.19458$
• $b:0.86449\to 0.88858$

而学习率为0.001，恰好$-0.34487\times -[0.001\ast (-150.28)]=-0.19458$

下面我们手动计算进行验证看看是否与上述输出一致

import numpy as np

x = np.array([4.6667, 8.2222, 2.8889])
y = np.array([9.7349, 17.2363, 6.3046])
y_ = np.array([-0.7449, -1.9711, -0.1318])

# 计算loss
f = y_-y
l = np.power(f, 2)
print(l)
l_ = np.sum(l)
print(l_/3)
'''
[109.82620804 368.92421476  41.42724496]
173.39255591999998
'''


# 计算相应参数
print(np.sum(x))
print(np.power(x, 2))
print(x*y)
print(np.sum(np.power(x, 2)))
print(np.sum(x*y))
'''
15.7778
[21.77808889 67.60457284  8.34574321]
[ 45.42985783 141.72030586  18.21335894]
97.72840494
205.36352263000003
'''


# 计算w梯度
gw = (2/3)*((-0.3449 * 97.7284)-(205.3635-15.7778*0.8645))
print(gw)
'''
-150.28674470666664
'''

可以发现我们的计算是一致的

三、执行代码

'''
    torch.optim.SGD
'''

import os
import torch
import torch.nn as nn
import matplotlib.pyplot as plt
from torch.utils.data import Dataset, DataLoader


# 1. 创建数据集
def create_linear_data(nums_data, if_plot=False):
    """
    Create data for linear model
    Args:
        nums_data: how many data points that wanted
    Returns:
        x with shape (nums_data, 1)
    """
    x = torch.linspace(2, 10, nums_data)
    x = torch.unsqueeze(x,dim=1)
    k = 2
    y = k * x + torch.rand(x.size())

    if if_plot:
        plt.scatter(x.numpy(),y.numpy(),c=x.numpy())
        plt.show()
    # data = torch.cat([x, y], dim=1)
    datax = x
    datay = y
    return datax, datay


datax, datay = create_linear_data(10, if_plot=False)
length = len(datax)


# 2. Dummy DataSet
class LinearDataset(Dataset):

    def __init__(self, length, datax, datay):
        self.len = length
        self.datax = datax
        self.datay = datay

    def __getitem__(self, index):
        datax = self.datax[index]
        datay = self.datay[index]
        return {'x': datax, 'y': datay}

    def __len__(self):
        return self.len


# 3. Parameters and DataLoaders
batch_size = 3
device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")
linear_loader = DataLoader(dataset=LinearDataset(length, datax, datay),
                         batch_size=batch_size, shuffle=True)


# 4. model
class Model(torch.nn.Module):
    """
    Linear Regressoin Module, the input features and output
    features are defaults both 1
    """
    def __init__(self):
        super().__init__()
        self.linear = torch.nn.Linear(1, 1)

    def forward(self, input):
        output = self.linear(input)
        print("In Model: #############################################")
        print("input: ", input.size(), "output: ", output.size())
        return output

model = Model()
# model = torch.nn.DataParallel(model, device_ids=device_ids)
model.cuda()


# 6. loss && optimizer
loss = torch.nn.MSELoss()
optimizer = torch.optim.SGD(model.parameters(), lr=0.001)


# 7. train
Linear_loss = []
for i in range(10):
    # print("epoch: %d\n", i)
    for data in linear_loader:
        x, y = data['x'], data['y']
        x = x.cuda()
        y = y.cuda()
        y_pred = model(x)
        lloss = loss(y_pred, y)
        optimizer.zero_grad()
        lloss.backward()
        for param in model.parameters():  # before optimize
            print(param.item())
        optimizer.step()
        for param in model.parameters():  # after optimize
            print(param.item())
        Linear_loss.append(lloss.item())
        print("Out--->input size:", x.size(), "output_size:", y_pred.size())


# 8. plot
loss_len = len(Linear_loss)
axis_x = range(loss_len)
plt.plot(axis_x, Linear_loss)
plt.show()