Back Propagation

• 算法特征
  ①. 统一看待线性运算与非线性运算; ②. 确定求导变量loss影响链路; ③. loss影响链路梯度逐级反向传播.
• 算法推导
  Part Ⅰ
  以如下简单正向传播链为例, 引入线性运算与非线性运算符号,
  相关运算流程如下,
  $$
  egin{equation*}
  egin{split}
  & ext{linear operation } & I^{(l+1)} = W^{(l)}cdot O^{(l)} + b^{(l)} \
  & ext{non-linear operation }quad & O^{(l+1)} = f^{(l+1)}(I^{(l+1)})
  end{split}
  end{equation*}
  $$
  其中, $O^{(l)}$、$I^{(l+1)}$、$O^{(l+1)}$分别为第$l$层输出(output)、第$l+1$层输入(input)、第$l+1$层输出, $W^{(l)}$、$b^{(l)}$、$f^{(l+1)}$分别为相关weight、bias及activation function. 对于线性运算, bias可合并至weight, 则
  $$
  egin{equation*}
  ext{linear operation }quad I^{(l+1)} = ilde{W}^{(l)}cdot ilde{O}^{(l)} = g^{(l)}( ilde{O}^{(l)})
  end{equation*}
  $$
  其中, $ ilde{W}^{(l)} = [W^{(l)mathrm{T}}, b^{(l)}]^mathrm{T}$, $ ilde{O}^{(l)}=[O^{(l)mathrm{T}},1]^mathrm{T}$.
  对于上述简单正向传播链, 基础影响链路如下,
  $$
  egin{equation*}
  egin{split}
  & ilde{O}^{(l)} &quad ightarrowquad I^{(l+1)} \
  & ilde{W}^{(l)} &quad ightarrowquad I^{(l+1)} \
  &I^{(l+1)} &quad ightarrowquad O^{(l+1)}
  end{split}
  end{equation*}
  $$
  根据链式法则, 影响链路起点由求导变量决定, 终点位于Loss处.
  Part Ⅱ
  现以如下神经网络为例, 加以阐述,
  

  其中, $(x_1, x_2)$为网络输入, $(y_1, y_2)$为网络输出. 统一处理其中线性运算与非线性运算, 并将该网络完全展开,
  

  现以$w_1^{(0)}$为例, 作为求导变量确定loss影响链路. $w_1^{(0)}$之loss影响链路(红色+绿色箭头)如下,
  

  拆解为基础影响链路, 如下,
  $$
  egin{equation*}
  egin{split}
  & w_1^{(0)} &quad ightarrowquad I_1^{(1)} \
  & I_1^{(1)} &quad ightarrowquad O_1^{(1)} \
  & O_1^{(1)} &quad ightarrowquad I_1^{(2)} \
  & O_1^{(1)} &quad ightarrowquad I_2^{(2)} \
  & I_1^{(2)} &quad ightarrowquad y_1 \
  & I_2^{(2)} &quad ightarrowquad y_2 \
  & y_1 &quad ightarrowquad L \
  & y_2 &quad ightarrowquad L
  end{split}
  end{equation*}
  $$
  定义源(source)为存在多个下游分支的节点(如: $O_1^{(1)}$), 定义汇(sink)为存在多个上游分支的节点(如: $L$). 根据链式求导法则, 在影响链路上, 下游变量对源变量求导需要在源变量处求和, 汇变量对上游变量求导无需特殊处理, 即,
  $$
  egin{equation*}
  egin{split}
  & ext{source variable $x$: }quad &frac{partial J(u(x), v(x))}{partial x} &= frac{partial J(u(x), v(x))}{partial u}cdotfrac{partial u}{partial x} + frac{partial J(u(x), v(x))}{partial v}cdotfrac{partial v}{partial x} \
  & ext{sink variable $J$: }quad &frac{partial J(u(x), v(y))}{partial x} &= frac{partial J(u(x), v(y))}{partial u}cdotfrac{partial u}{partial x}
  end{split}
  end{equation*}
  $$
  由此, 根据梯度沿影响链路的反向传播, 可确定loss对$w_1^{(0)}$之偏导,
  $$
  egin{equation*}
  frac{partial L}{partial w_1^{(0)}} = left(frac{partial L}{partial y_1}cdotfrac{partial y_1}{partial I_1^{(2)}}cdotfrac{partial I_1^{(2)}}{partial O_1^{(1)}} + frac{partial L}{partial y_2}cdotfrac{partial y_2}{partial I_2^{(2)}}cdotfrac{partial I_2^{(2)}}{partial O_1^{(1)}} ight)cdotfrac{partial O_1^{(1)}}{partial I_1^{(1)}}cdotfrac{partial I_1^{(1)}}{partial w_1^{(0)}}
  end{equation*}
  $$
  再以$b_2^{(1)}$为例, 作为求导变量确定loss影响链路. $b_2^{(1)}$之loss影响链路(红色+绿色箭头)如下,
  拆解为基础影响链路, 如下,
  $$
  egin{equation*}
  egin{split}
  & b_2^{(1)} &quad ightarrowquad I_2^{(2)} \
  & I_2^{(2)} &quad ightarrowquad y_2 \
  & y_2 &quad ightarrowquad L
  end{split}
  end{equation*}
  $$
  同样, 根据梯度沿影响链的反向传播, 可确定loss对$b_2^{(1)}$之偏导,
  $$
  egin{equation*}
  frac{partial L}{partial b_2^{(1)}} = frac{partial L}{partial y_2}cdotfrac{partial y_2}{partial I_2^{(2)}}cdotfrac{partial I_2^{(2)}}{partial b_2^{(1)}}
  end{equation*}
  $$
• 代码实现
  现以如下简单feed-forward网络为例进行算法实施,
  输入层为$(r, g, b)$, 输出层为$(x,y,lv)$且不取激活函数, 中间隐藏层取激活函数为双曲正切函数$ anh$. 采用如下损失函数,
  $$
  egin{equation*}
  L = sum_ifrac{1}{2}(ar{x}^{(i)}-x^{(i)})^2 + frac{1}{2}(ar{y}^{(i)} - y^{(i)})^2 + frac{1}{2}(ar{lv}^{(i)} - lv^{(i)})^2
  end{equation*}
  $$
  其中, $i$为data序号, $(ar{x}, ar{y}, ar{lv})$为相应观测值. 相关training data采用如下策略生成,
  $$
  egin{equation*}
  left{
  egin{split}
  x &= r + 2g + 3b \
  y &= r^2 + 2g^2 + 3b^2 \
  lv &= -3r - 4g - 5b
  end{split}
  ight.
  end{equation*}
  $$
  具体实现如下,
  

  # Back Propagation之实现
  # 优化器使用Adam
  
  import numpy
  from matplotlib import pyplot as plt
  
  
  numpy.random.seed(1)
  
  
  # 生成training data
  def getData(n=100):
      rgbRange = (-1, 1)
      r = numpy.random.uniform(*rgbRange, (n, 1))
      g = numpy.random.uniform(*rgbRange, (n, 1))
      b = numpy.random.uniform(*rgbRange, (n, 1))
      x_ = r + 2 * g + 3 * b
      y_ = r ** 2 + 2 * g ** 2 + 3 * b ** 2
      lv_ = -3 * r - 4 * g - 5 * b
      RGB = numpy.hstack((r, g, b))
      XYLv_ = numpy.hstack((x_, y_, lv_))
      return RGB, XYLv_
      
      
  class BPEx(object):
  
      def __init__(self, RGB, XYLv_):
          self.__RGB = RGB
          self.__XYLv_ = XYLv_
          
          self.__rgb = None           # (1, 3)
          self.__O0 = None            # (1, 4)
          self.__W0 = None            # (4, 5)
          self.__I1 = None            # (1, 5)
          self.__O1 = None            # (1, 6)
          self.__W1 = None            # (6, 3)
          self.__I2 = None            # (1, 3)
          self.__xylv_ = None         # (1, 3)
          self.__L = None             # scalar
          
          self.__grad_W0_I1 = None    # 基础影响链路(W0->I1)之梯度
          self.__grad_I1_O1 = None    # 基础影响链路(I1->O1)之梯度
          self.__grad_O1_I2 = None    # 基础影响链路(O1->I2)之梯度
          self.__grad_W1_I2 = None    # 基础影响链路(W1->I2)之梯度
          self.__grad_I2_L = None     # 基础影响链路(I2->Loss)之梯度
          
          self.__gradList_W0_I1 = list()
          self.__gradList_I1_O1 = list()
          self.__gradList_O1_I2 = list()
          self.__gradList_W1_I2 = list()
          self.__gradList_I2_L = list()
          
          self.__init_weights()
          self.__bpTag = False        # 是否反向传播之标志
          
          
      def calc_xylv(self, rgb, W0=None, W1=None):
          '''
          默认在当前W0、W1之基础上进行预测, 实际使用需要配合优化器
          '''
          if W0 is not None:
              self.__W0 = W0
          if W1 is not None:
              self.__W1 = W1
          
          self.__calc_xylv(numpy.array(rgb).reshape((1, 3)))
          xylv = self.__I2
          return xylv
          
          
      def get_W0W1(self):
          return self.__W0, self.__W1
          
          
      def calc_JVal(self, W):
          self.__bpTag = True
          self.__clr_gradList()
      
          self.__W0 = W[:20, 0].reshape((4, 5))
          self.__W1 = W[20:, 0].reshape((6, 3))
          
          JVal = 0
          for rgb, xylv_ in zip(self.__RGB, self.__XYLv_):
              self.__calc_xylv(rgb.reshape((1, 3)))
              self.__calc_loss(xylv_.reshape((1, 3)))
              JVal += self.__L
              self.__add_gradList()
              
          self.__bpTag = False
          return JVal
          
          
      def calc_grad(self, W):
          '''
          此处W仅为统一调用接口
          '''
          grad_W0 = numpy.zeros(self.__W0.shape)
          grad_W1 = numpy.zeros(self.__W1.shape)
          for grad_W0_I1, grad_I1_O1, grad_O1_I2, grad_W1_I2, grad_I2_L in zip(self.__gradList_W0_I1,
              self.__gradList_I1_O1, self.__gradList_O1_I2, self.__gradList_W1_I2, self.__gradList_I2_L):
              grad_W1_curr = grad_I2_L * grad_W1_I2
              grad_W1 += grad_W1_curr
              
              term0 = numpy.sum(grad_I2_L.T * grad_O1_I2, axis=0)   # (1, 5) source
              term1 = term0 * grad_I1_O1
              grad_W0_curr = term1 * grad_W0_I1
              grad_W0 += grad_W0_curr
              
          grad = numpy.vstack((grad_W0.reshape((-1, 1)), grad_W1.reshape((-1, 1))))
          return grad
          
          
      def init_seed(self):
          n = self.__W0.size + self.__W1.size
          seed = numpy.random.uniform(-1, 1, (n, 1))
          return seed
          
          
      def __add_gradList(self):
          self.__gradList_W0_I1.append(self.__grad_W0_I1)
          self.__gradList_I1_O1.append(self.__grad_I1_O1)
          self.__gradList_O1_I2.append(self.__grad_O1_I2)
          self.__gradList_W1_I2.append(self.__grad_W1_I2)
          self.__gradList_I2_L.append(self.__grad_I2_L)
          
          
      def __clr_gradList(self):
          self.__gradList_W0_I1.clear()
          self.__gradList_I1_O1.clear()
          self.__gradList_O1_I2.clear()
          self.__gradList_W1_I2.clear()
          self.__gradList_I2_L.clear()
          
          
      def __init_weights(self):
          self.__W0 = numpy.zeros((4, 5))
          self.__W1 = numpy.zeros((6, 3))
          
          
      def __update_grad_by_calc_xylv(self):
          self.__grad_W0_I1 = numpy.tile(self.__O0.T, (1, 5))
          self.__grad_I1_O1 = 1 - self.__I1_tanh ** 2
          self.__grad_O1_I2 = self.__W1.T[:, :-1]
          self.__grad_W1_I2 = numpy.tile(self.__O1.T, (1, 3))
          
          
      def __calc_xylv(self, rgb):
          '''
          rgb: shape=(1, 3)之numpy array
          '''
          self.__rgb = rgb
          self.__O0 = numpy.hstack((self.__rgb, [[1]]))
          self.__I1 = numpy.matmul(self.__O0, self.__W0)
          self.__I1_tanh = numpy.tanh(self.__I1)
          self.__O1 = numpy.hstack((self.__I1_tanh, [[1]]))
          self.__I2 = numpy.matmul(self.__O1, self.__W1)
          
          if self.__bpTag:
              self.__update_grad_by_calc_xylv()
              
              
      def __update_grad_by_calc_loss(self):
          self.__grad_I2_L = self.__I2 - self.__xylv_
              
              
      def __calc_loss(self, xylv_):
          '''
          xylv_: shape=(1, 3)之numpy array
          '''
          self.__xylv_ = xylv_
          self.__L = numpy.sum((self.__xylv_ - self.__I2) ** 2) / 2
          
          if self.__bpTag:
              self.__update_grad_by_calc_loss()
  
  
  class Adam(object):
  
      def __init__(self, _func, _grad, _seed):
          '''
          _func: 待优化目标函数
          _grad: 待优化目标函数之梯度
          _seed: 迭代起始点
          '''
          self.__func = _func
          self.__grad = _grad
          self.__seed = _seed
  
          self.__xPath = list()
          self.__JPath = list()
  
  
      def get_solu(self, alpha=0.001, beta1=0.9, beta2=0.999, epsilon=1.e-8, zeta=1.e-6, maxIter=3000000):
          '''
          获取数值解,
          alpha: 步长参数
          beta1: 一阶矩指数衰减率
          beta2: 二阶矩指数衰减率
          epsilon: 足够小正数
          zeta: 收敛判据
          maxIter: 最大迭代次数
          '''
          self.__init_path()
  
          x = self.__init_x()
          JVal = self.__calc_JVal(x)
          self.__add_path(x, JVal)
          grad = self.__calc_grad(x)
          m, v = numpy.zeros(x.shape), numpy.zeros(x.shape)
          for k in range(1, maxIter + 1):
              print("iterCnt: {:3d},   JVal: {}".format(k, JVal))
              if self.__converged1(grad, zeta):
                  self.__print_MSG(x, JVal, k)
                  return x, JVal, True
  
              m = beta1 * m + (1 - beta1) * grad
              v = beta2 * v + (1 - beta2) * grad * grad
              m_ = m / (1 - beta1 ** k)
              v_ = v / (1 - beta2 ** k)
  
              alpha_ = alpha / (numpy.sqrt(v_) + epsilon)
              d = -m_
              xNew = x + alpha_ * d
              JNew = self.__calc_JVal(xNew)
              self.__add_path(xNew, JNew)
              if self.__converged2(xNew - x, JNew - JVal, zeta ** 2):
                  self.__print_MSG(xNew, JNew, k + 1)
                  return xNew, JNew, True
  
              gNew = self.__calc_grad(xNew)
              x, JVal, grad = xNew, JNew, gNew
          else:
              if self.__converged1(grad, zeta):
                  self.__print_MSG(x, JVal, maxIter)
                  return x, JVal, True
  
          print("Adam not converged after {} steps!".format(maxIter))
          return x, JVal, False
  
  
      def get_path(self):
          return self.__xPath, self.__JPath
  
  
      def __converged1(self, grad, epsilon):
          if numpy.linalg.norm(grad, ord=numpy.inf) < epsilon:
              return True
          return False
  
  
      def __converged2(self, xDelta, JDelta, epsilon):
          val1 = numpy.linalg.norm(xDelta, ord=numpy.inf)
          val2 = numpy.abs(JDelta)
          if val1 < epsilon or val2 < epsilon:
              return True
          return False
  
  
      def __print_MSG(self, x, JVal, iterCnt):
          print("Iteration steps: {}".format(iterCnt))
          print("Solution:
  {}".format(x.flatten()))
          print("JVal: {}".format(JVal))
  
  
      def __calc_JVal(self, x):
          return self.__func(x)
  
  
      def __calc_grad(self, x):
          return self.__grad(x)
  
  
      def __init_x(self):
          return self.__seed
  
  
      def __init_path(self):
          self.__xPath.clear()
          self.__JPath.clear()
  
  
      def __add_path(self, x, JVal):
          self.__xPath.append(x)
          self.__JPath.append(JVal)
          
          
  class BPExPlot(object):
  
      @staticmethod
      def plot_fig(adamObj):
          alpha = 0.001
          epoch = 50000
          x, JVal, tab = adamObj.get_solu(alpha=alpha, maxIter=epoch)
          xPath, JPath = adamObj.get_path()
          
          fig = plt.figure(figsize=(6, 4))
          ax1 = fig.add_subplot(1, 1, 1)
          
          ax1.plot(numpy.arange(len(JPath)), JPath, "k.", markersize=1)
          ax1.plot(0, JPath[0], "go", label="seed")
          ax1.plot(len(JPath)-1, JPath[-1], "r*", label="solution")
          
          ax1.legend()
          ax1.set(xlabel="$epoch$", ylabel="$JVal$", title="solution-JVal = {:.5f}".format(JPath[-1]))
          
          fig.tight_layout()
          # plt.show()
          fig.savefig("plot_fig.png", dpi=100)
          
          
          
  if __name__ == "__main__":
      RGB, XYLv_ = getData(1000)
      bpObj = BPEx(RGB, XYLv_)
      
      # rgb = (0.5, 0.6, 0.7)
      # xylv = bpObj.calc_xylv(rgb)
      # print(rgb)
      # print(xylv)
      # func = bpObj.calc_JVal
      # grad = bpObj.calc_grad
      # seed = bpObj.init_seed()
      # adamObj = Adam(func, grad, seed)
      # alpha = 0.1
      # epoch = 1000
      # adamObj.get_solu(alpha=alpha, maxIter=epoch)
      # xylv = bpObj.calc_xylv(rgb)
      # print(rgb)
      # print(xylv)
      # W0, W1 = bpObj.get_W0W1()
      # xylv = bpObj.calc_xylv(rgb, W0, W1)
      # print(rgb)
      # print(xylv)
      
      func = bpObj.calc_JVal
      grad = bpObj.calc_grad
      seed = bpObj.init_seed()
      adamObj = Adam(func, grad, seed)
      BPExPlot.plot_fig(adamObj)

• 结果展示
  

  可以看到, 在training data上总体loss随epoch增加逐渐降低.

• 使用建议
  ①. 某层之output节点可能受多个该层之input节点影响(如: softmax激活函数), 此时input节点具备source特性;
  ②. 正向计算过程可确定基础影响链路之梯度, 反向传播过程可串联基础影响链路之梯度;
  ③. 初值的选取对非凸问题的优化比较重要, 权重初值尽量不取全0.
• 参考文档
  ①. Rumelhart D E, Hinton G E, Williams R J. Learning internal representations by error propagation[R]. California Univ San Diego La Jolla Inst for Cognitive Science, 1985.
  ②. Adam (1) - Python实现