CMU 10-414/714: Deep Learning Systems --hw1

hw1

needle库中的两个重要文件

1. autograd.py文件：实现计算图框架、自动微分框架（automatic differentiation framework）
2. ops.py文件：实现各种算子。反向传播的粒度很细，会以一个个算子为单位，而不是以网络层为单位。其实现基础是numpy，后面会自己实现自己的运算库

needle中的一些核心类

1. Value：计算图中的节点，代表计算过程中的数值
2. Op：计算图的算子，即加减乘除等操作。其通过compute()函数进行前向计算，通过gradient()函数进行反向传播
3. Tensor：是Value的子类，表示计算图中的一个多维数组。课程提供了一些功能（如操作符重载）来保证用tensor进行数学运算
4. TensorOp：是Op的子类，表示返回Tensor的算子

实现功能

1. 实现前向传播和反向传播
   • compute()函数实现前向运算，运算的输入都是numpy.ndarray（后面的实验会实现自己的NDArray，因此这里为了对齐用import numpy as array_api），即compute()函数只处理raw data，而不处理自动求导的Tensor
   • gradient()函数实现反向传播，其输入为Tensor，且函数内的运算都是TensorOp算子
   • 需要实现运算的算子如下。这里代码其实没有什么高深的，只要搞清楚每个算子是怎么求导的就好了，比如PowerScalar算子实现$a^{scalar}$，其gradient()函数就是返回$(out\_grad\times self.scalar\times lh{s}^{self.scalar-1})$
     • PowerScalar：乘方
     • EWiseDiv：除法
     • DivScalar：张量除以标量
     • MatMul：矩阵乘法
       • 梯度必须和输入数据的维度一致：损失函数L对输入数据X的梯度，告诉我们在每个特征维度上，改变输入数据时损失函数的变化情况。若想要计算L对每个输入特征的梯度，这些梯度将组成一个与输入数据X相同形状的张量，每个元素表示在相应特征维度上的变化率
     • Summation：在指定的维度上求和
     • BroadcastTo：broadcast成指定的形状
     • Reshape：转变成之后定的形状
     • Negate：取负数
     • Transpose：交换两个维度
2. 拓扑排序：原理很简单，数据结构里学过。就是深度优先遍历node_list，如果一个node没有被visited，就深度优先遍历它的inputs结点，直到一个入度为0的结点，就将其加入topo_order中

def find_topo_sort(node_list: List[Value]) -> List[Value]:
    visited = set()
    topo_order = []
    for node in node_list:
      if node not in visited:
        topo_sort_dfs(node, visited, topo_order)
    return topo_order

def topo_sort_dfs(node, visited, topo_order):
    if node in visited:
      return

    for next in node.inputs:
      topo_sort_dfs(next, visited, topo_order)
    
    visited.add(node)
    topo_order.append(node)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

3. 实现reverse mode differentiation：

def compute_gradient_of_variables(output_tensor, out_grad):
  # 创建一个空字典，用于存储每个节点在计算图中的梯度信息
  node_to_output_grads_list: Dict[Tensor, List[Tensor]] = {
   }
  # 将输出张量的梯度存储在字典中，作为输出节点的梯度
  node_to_output_grads_list[output_tensor] = [out_grad]

  # find_topo_sort([output_tensor])从输出张量开始计算图的拓扑排序（因为输出张量通过node.inputs连着整张图，所以可以实现）
  reverse_topo_order = list(reversed(find_topo_sort([output_tensor])))

  for node in reverse_topo_order:
    # 计算当前节点的梯度：就是把所有偏导数加到一起
    sum_grad = node_to_output_grads_list[node][0]
    for t in node_to_output_grads_list[node][1:]:
      sum_grad = sum_grad + t 
    node.grad = sum_grad
	
	# 叶子节点即没有op的节点，跳过后面的for循环
    if node.is_leaf():
      continue

	# 计算节点的输入节点相对于
    for i, grad in enumerate(node.op.gradient_as_tuple(node.grad, node)):
      input_ = node.inputs[i]
      if input_ not in node_to_output_grads_list:
        node_to_output_grads_list[input_] = []
      node_to_output_grads_list[input_].append(grad)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

举个例子有助于理解自动微分的过程原理：

4. 重新实现softmax loss：使用Tensor且y是one-hot向量组成的$m\times k$矩阵

def softmax_loss(Z, y_one_hot):
    a = ndl.ops.summation(Z*y_one_hot)
    b = ndl.ops.summation(ndl.ops.log(ndl.ops.summation(ndl.ops.exp(Z), axes=(1,))))
    return (b-a) / Z.shape[0]
1
2
3
4

5. 重新为一个两次神经网络实现SGD
   与hw0的区别：W1和W2变成Tensor，X和y还是numpy array；hw0中是手动计算反向传播公式然后用代码实现，这次则是使用自动微分机制，只需算出loss值，然后调用.backward()即可自动对每个参数求导，再对W1和W2更新即可

def nn_epoch(X, y, W1, W2, lr = 0.1, batch=100):
    for i in range(X.shape[0]//batch+1):
      # 首先不变的，采样batch
      start, end = i * batch, min((i+1)*batch, X.shape[0])
      m = end - start
      if m == 0:
          break
      Xb, yb = X[start:end], y[start:end]

      # hw0这里是计算梯度，是通过手动计算得出数学公式然后实现的
      # 现在则是使用自动微分机制，算出loss后直接.backward()即可
      # 下面并没有对X和y设置自动微分，因为参数更新过程中只需要对参数矩阵W1和W2自动微分即可，不需要对输入数据和得出来的结果微分（node_to_grad第一个值是loss）
      Xb, yb = X[start:end], y[start:end]
      Xb = ndl.Tensor(Xb, requires_grad=False)
      Z1 = Xb @ W1
      A1 = ndl.ops.relu(Z1)
      Z2 = A1 @ W2
      y_one_hot = np.zeros(Z2.shape, dtype=np.float32)
      y_one_hot[np.arange(Z2.shape[0]), yb] = 1
      y_one_hot = ndl.Tensor(y_one_hot, requires_grad=False)
      loss = softmax_loss(Z2, y_one_hot)
      loss.backward()
      
      # 更新参数矩阵
      W1.data = W1.data - lr * ndl.Tensor(W1.grad.numpy().astype(np.float32)) 
      W2.data = W2.data - lr * ndl.Tensor(W2.grad.numpy().astype(np.float32))
    return W1, W2
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

相关知识补充

一、pytorch中的正向传播、反向传播和计算图

1. Tensor
   pytorch提供的autograd包能够根据输入和前向传播过程自动构建计算图，并执行反向传播。
   
   Tensor是autograd包的核心类，若将其属性.requires_grad设为True，它将开始追踪在其上的所有操作（这样就可以利用链式法则进行梯度传播了）。完成计算后，可调用.backward()来完成所有的梯度计算。此Tensor的梯度将累积到.grad属性中。
   
   Function是另一个核心类，Tensor和Function结合就可以构建一个记录有整个计算过程的有向无环图。每个Tensor的.grad_fn属性记录着这个Tensor最近是通过哪一个运算得到的，若其不是通过某个运算得到的，.grad_fn属性值为None，像grad_fn属性值为None的称为叶子节点

x = torch.ones(2, 2, requires_grad=True)
y = x + 2
z = y * y * 3
out = z.mean()
print(z, out)
1
2
3
4
5

输出为：

tensor([[27., 27.],
        [27., 27.]], grad_fn=<MulBackward>) tensor(27., grad_fn=<MeanBackward1>)
1
2

2. 梯度
   
   $$out=\frac{1}{4}\sum _{i=1}^4{z}_i=\frac{1}{4}\sum _{i=1}^{4}3(x_i+2{)}^2$$ 所以
   $$\frac{\partial out}{\partial x_i}{|}_{x_i=1}=\frac{9}{2}=4.5$$
   
   数学上，若有一个函数值和自变量都为向量的函数$\vec{y}=f(\vec{x})$，那么$\vec{y}$关于$\vec{x}$的梯度就是一个雅可比矩阵：
   J = ( ∂ y 1 ∂ x 1 ⋯ ∂ y 1 ∂ x n ⋮ ⋱ ⋮ ∂ y m ∂ x 1 ⋯ ∂ y m ∂ x n ) J=\left(
   ∂y1∂x1⋯∂y1∂xn⋮⋱⋮∂ym∂x1⋯∂ym∂xn" role="presentation" style="position: relative;">∂y1∂x1⋮∂ym∂x1⋯⋱⋯∂y1∂xn⋮∂ym∂xn$$\begin{array}{ccc}\frac{\mathrm{\partial }{y}_1}{\mathrm{\partial }{x}_1}& \cdots & \frac{\mathrm{\partial }{y}_1}{\mathrm{\partial }{x}_n}\\ ⋮& \ddots & ⋮\\ \frac{\mathrm{\partial }{y}_m}{\mathrm{\partial }{x}_1}& \cdots & \frac{\mathrm{\partial }{y}_m}{\mathrm{\partial }{x}_n}\end{array}$$
   \right) J=