符号
意义
$m$ 样本个数。
$k$ 样本标签向量的维度。
$T$ 循环神经网络的时间步数，在矩阵右上角时表示矩阵的转置。
$d^l$ 神经网络第 l 层神经元的个数。
$\sigma(·)$ 是某一种常用的激活函数。比如 sigmoid、tanh 和 ReLu 函数等。
$\rho(·)$ 神经网络输出层的变换函数。比如 logistic 和 Softmax 函数等。
$\Omega(·)$ 循环神经网络不同时间步同一层节点梯度距离函数。
$\Psi(·)$ 循环神经网络同一时间步相邻层节点梯度距离函数。
$\mathcal{O}(·)$ 时间复杂度。
$C$ 代价函数。
$Loss$ 损失函数，意义同 $C$ 。
$sum(·)$ 求和函数。可以根据矩阵行或列求和，比如 $sum(D, axis=0)$ 就是按照矩阵 $D$ 的行求和。
$where(·)$ 条件函数。可以根据矩阵元素满足的条件操作，比如 ${Y}=where(A>=0.5, 1, 0)$ 就是若矩阵 $A$ 中元素 $a_{ij}$ 大于等于 0.5 则 $Y$ 矩阵中对应位置元素 $y_{ij}$ 的值取 1 否则取 0。
$*$ 按矩阵元素乘积。比如 $A * B$ 。
$\odot$ 按矩阵元素乘积。比如 $A \odot B$ 。
符号
向量形式
维度
意义
$x$ $[x_0,x_1,...,x_n]^T$ $(n, 1)$
$x$ 表示一个有 $n$ 维特征样本实例， $x_i$ 表示向量的第 $i$ 维特征。
$w_{jk}^{l}$ 标量
$(1, 1)$ 是从 $l-1$ 层的第 $k$ 个神经元到 $l$ 层的第 $j$ 个神经元的权重。
$s_{jk}^{l}$ 标量
$(1, 1)$ 是第 $l$ 层的记忆向量的第 $k$ 个元素到 $l$ 层的第 $j$ 个神经元的权重。
$b_j^l$ 标量
$(1, 1)$ 是第 $l$ 层的第 $j$ 个神经元的偏置。
$z_j^l$ 标量
$(1, 1)$ 是第 $l$ 层的第 $j$ 个神经元的带权输入。
$a_j^l$ 标量
$(1, 1)$ 是第 $l$ 层的第 $j$ 个神经元的激活值。
$y_j$ 标量
$(1, 1)$ 样本标签向量的第 $j$ 个值。
$\hat y_j$ 标量
$(1, 1)$ 模型预测的样本标签向量的第 $j$ 个值。
$\delta_j^l$ 标量
$(1, 1)$ 是第 $l$ 层的第 $j$ 个神经元的误差值。
$z^l$ $[z_0,z_1,...,z_{d^l}]^T$ $(d^l, 1)$ 表示第第 $l$ 层的带权输入。第 $j$ 个元素是 $z_j^l$ 。
$a^{l}$ $[a_0,a_1,...,a_{d^l}]^T$ $(d^l, 1)$ 是激活向量。第 $j$ 行的元素是 $a_j^l$ 。
$b^l$ $[b_0,b_1,...,b_{d^l}]^T$ $(d^l, 1)$ 偏置向量，第 $j$ 行的元素是 $b_j^l$
$h^l$ $[h_0,h_1,...,h_{d^l}]^T$ $(d^l, 1)$ 记忆向量，第 $j$ 行的元素是 $h_j^l$
$W^{l}$ $\left[\begin{matrix} w_{i1}^l & w_{i2}^l &···& w_{ik}^l & \ \w_{j1}^l & w_{j2}^l &···& w_{jk}^l\end{matrix} \right]$
$(d^l, d^{l-1})$ 权重矩阵，第 $j$ 行 $k$ 列的元素是 $w_{jk}^{l}$ 。
$S^{l}$ $\left[\begin{matrix} s_{i1}^l &···& s_{ik}^l & \ \s_{k1}^l &···& s_{kk}^l\end{matrix} \right]$
$(d^l, d^{l})$ 记忆向量的权重方阵，第 $j$ 行 $k$ 列的元素是 $s_{jk}^{l}$ 。
$\delta^l$ $[\delta_0,\delta_1,...,\delta_{d^l}]^T$ $(d^l, 1)$ 是误差向量，$\delta^l$ 的第 $j$ 个元素是 $\delta_j^l$ 。
$z^{v,l}$ $[z_0,z_1,...,z_{d^l}]^T$ $(d^l, 1)$ 表示神经网络第 $l$ 层的小批量样本中的第 $v$ 个样本的带权输入向量。
$a^{v,l}$ $[a_0,a_1,...,a_{d^l}]^T$ $(d^l, 1)$ 表示神经网络第 $l$ 层的小批量样本中的第 $v$ 个样本的激活向量。
$\delta^{v,l}$ $[\delta_0,\delta_1,...,\delta_{d^l}]^T$ $(d^l, 1)$ 表示神经网络第 $l$ 层的小批量样本中的第 $v$ 个样本的误差。
$y$ $[y_0,y_1,...,y_{k}]^T$ $(k, 1)$ 标签向量，第 $j$ 个值是 $y_j$ 。
$\hat y$ $[\hat y_0,\hat y_1,...,\hat y_{k}]^T$ $(k, 1)$ 模型预测的标签向量，第 $j$ 个值是 $\hat y_j$ 。
$X$ $[{(x^{(1)})}^T,{(x^{(2)})}^T,...,{(x^{(m)})}^T]^T$ $(m, n)$ 样本矩阵，第 $i$ 行是第 $i$ 个样本向量 $x^{(i)}$ 的转置。
$Y$ $[{(y^{(1)})}^T,{(y^{(2)})}^T,...,{(y^{(m)})}^T]^T$ $(m, k)$ 标签矩阵，第 $i$ 行是第 $i$ 个标签向量 $y^{(i)}$ 的转置。
$\hat Y$ $[{(\hat y^{(1)})}^T,{(\hat y^{(2)})}^T,...,{(\hat y^{(m)})}^T]^T$ $(m, k)$ 估计的标签矩阵，第 $i$ 行是第 $i$ 个标签向量 $\hat y^{(i)}$ 的转置。
$Z^l$ $[{(z^{1,l})}^T,{(z^{2,l})}^T,...,{(z^{m,l})}^T]^T$ $(m, d^l)$ 带权输入矩阵，第 $v$ 行是第 $v$ 个样本在第 $l$ 层的带权输入向量 $z^{v,l}$ 的转置。
$A^l$ $[{(a^{1,l})}^T,{(a^{2,l})}^T,...,{(a^{m,l})}^T]^T$ $(m, d^l)$ 激活矩阵，第 $v$ 行是第 $v$ 个样本在第 $l$ 层的激活向量 $a^{v,l}$ 的转置。
$\Delta^l$ $[{(\delta^{1,l})}^T,{(\delta^{2,l})}^T,...,{(\delta^{m,l})}^T]^T$ $(m, d^l)$ 误差矩阵，第 $v$ 行是第 $v$ 个样本在第 $l$ 层的误差向量 $\delta^{v,l}$ 的转置。
$\Omega(j, i)$ $(m, d^l)$ 循环神经网络不同时间步同一层节点梯度距离函数, 它的值是$\dfrac{\partial Z_{t_j}^{(l)}}{\partial Z_{t_i}^{(l)}}$。
$\Omega$ $(T, T)$ 梯度距离函数 $\Omega(·)$ 的状态转移矩阵，它的 $j$ 行 $i$ 列位置对应的元素是 $\Omega(j, i)$ 。
$\widetilde \Omega$ $(T, T+1)$
$\Omega$ 矩阵的增广矩阵。
无数的深度学习任务需要处理顺序数据。图像字幕，语音合成和音乐生成都要求模型产生作为序列的输出。在其他领域，例如时间序列预测，视频分析和音乐信息检索，模型必须从作为序列的输入中学习。交互式任务，例如翻译自然语言，参与对话和控制机器人，往往需要两种能力。递归神经网络（RNN）是连接模型，其通过节点网络中的循环捕获序列的动态。与标准前馈神经网络不同，循环网络保留可以表示来自任意长的上下文窗口的信息的状态。
循环神经网络引入了状态变量。在一个序列中，循环神经网络当前时刻的状态不仅保存了过去时刻的信息，还与当前时刻的输入共同决定当前时刻的输出。但随着时间步数 $T$ （循环次数）增加，循环神经网络越来越难训练。
循环层是 RNN 难以训练的主要原因，本文分析了计算循环神经网络反向传播算法中循环层的时间复杂度为 $\mathcal{O}(T^3)$ ， 为了降低计算复杂度，本文定义了两种梯度距离函数，然后由距离函数推导出 $\Omega$ 梯度距离矩阵，最后使用动态规划算法求解 $\Omega$ 梯度距离矩阵，在保证梯度沿时间步依赖关系反向传播的前提下使循环层的时间复杂度降低为 $\mathcal{O}(T^2)$ 。
如今训练 RNN 的主流方法是使用基于时间的反向传播算法 BPTT（Back Propagation Trough Time）。多年来对于 BPTT 算法的改进和应用领域研究也在不断地进行着。对于 BPTT 的算法改进主要有两种方式：一种是基于启发式进行的改进,一种是基于数值优化方式的改进。而基于数值优化的方式主要有共扼梯度法、LM 算法等。但是它们要么计算准确度有缺陷（启发式方法），要么计算量太大（共轭梯度法）。本文提出的动态规划算法来加速基于时间的反向传播算法，大大提高了计算速度，且保证了原有数值准确度。不仅如此，本文的动态规划方法同样可以应用到改进后的 BPTT 算法中，与其它改进算法有良好的相容性。本文研究的动态动态规划算法加速循环神经网络梯度的反向传播的方法已经代码实现，并且开源提供给其他研究者使用和进一步开发。
|计算公式|维度变换|
|-|-|-|
|$\Delta_{t=i}^L = A^L_{t=i} - Y_{t=i}$|$(m , d^L) = (m , d^L) - (m , d^L)$|
|$\Delta^{l}{t=T}=\Delta^{(l+1)} {t=T}W^{(l+1)} \odot \sigma'(Z^{(l)}{t=T}) $|$(m, d^{l})=(m, d^{l+1})(d^{l+1}, d^{l}) * (m, d^{l})$|
|$\Delta^{l} {t=T-i}=\Delta^{(l+1)}{t=T-i}W^{(l+1)} \odot \sigma'(Z^{(l)} {t=T-i}) + \sum_{j=0}^{i-1} \left{ \Delta^{(l)}{t=T-j}[(S^{(l)})^T]^{i-j} \prod {k=T-i}^{T-j-1}\sigma'(Z^{(l)}{t=k}) \right}$|$(m, d^{l})=(m, d^{l+1})(d^{l+1}, d^{l}) * (m, d^{l}) + (m, d^l)(d^l, d^l)^{i-j} * (m, d^l)^{i-j}$|
|$\dfrac{\partial C}{\partial b^{l}} =\dfrac{1}{m}[{\sum {t=i}^T \Delta_{t=i}^{(l)}}]^T $|$ (d^l, 1) = {(1, d^l)}^T$|
|$\dfrac{\partial C}{\partial W^{l}} = \dfrac{1}{m} \sum_{t=i}^T [\Delta {t}^{(l)}]^T A_t^{(l-1)}$|$(d^l, d^{(l-1)}) = {(m, d^l)}^T (m, d^{l-1})$|
|$\dfrac{\partial C}{\partial S^{l}} = \dfrac{1}{m} \sum {t=i}^T [\Delta_{t}^{(l)}]^T H_t^{(l)}$|$(d^l, d^{l}) = {(m, d^l)}^T (m, d^{l})$|
$\sigma'(Z_{t=i}^{(l)})$ 表示要用激活函数的导函数对矩阵 $Z$ 中的每一个元素进行计算，而矩阵 $Z$ 的维度为 $(m, d^l)$ 含有 $m * d^l$ 个元素且激活函数的导函数一般很复杂（可能含有指数运行），所以该操作非常耗时，是本文主要优化的对象。
对于梯度的递推计算式子，有几点需要说明，
$$\Delta^{l}{t=T-i}=\Delta^{(l+1)} {t=T-i}W^{(l+1)} \odot \sigma'(Z^{(l)}{t=T-i}) + \sum {j=0}^{i-1} \left{ \Delta^{(l)}{t=T-j}[(S^{(l)})^T]^{i-j} \prod {k=T-i}^{T-j-1}\sigma'(Z^{(l)}_{t=k}) \right}$$
$[S^{(l)}]^{i-j}$ 表示 $i-j$ 次方个 $S^{(l)}$ ，$\prod$ 使用的是 $\odot$ 连乘，即矩阵对应位置元素相乘。
为了简化上面的式子，同时为了看清计算式的本质，也是后面使用动态规划方法求解梯度递推计算式子的基础。
定义循环神经网络不同时间步同一层节点梯度距离函数,
$$\Omega(Z_{t_j}^{(l)}, Z_{t_i}^{(l)})=\dfrac{\partial Z_{t_j}^{(l)}}{\partial Z_{t_i}^{(l)}}=\prod_{t=t_i}^{t_{j}-1}(S^T)^{(l)} \odot \sigma'(Z_{t}^{(l)}) $$
定义循环神经网络同一时间步相邻层节点梯度距离函数,
$$\Psi(Z_{t_i}^{(l+1)}, Z_{t_i}^{(l)})=\dfrac{\partial Z_{t_i}^{(l+1)}}{\partial Z_{t_i}^{(l)}}=W^{(l+1)} \odot \sigma'(Z^{(l)}_{t_i})$$
梯度的递推计算式子就可以简化为，
$$\Delta^{(l)}{t=T-i}=\Delta^{(l+1)} {t=T-i}\Psi(Z_{t=T-i}^{(l+1)}, Z_{t=T-i}^{(l)}) + \sum_{j=0}^{i-1} \Delta^{(l)}{t=T-j}\Omega(Z {T-j}^{(l)}, Z_{T-i}^{(l)}) $$
该式子的意义就是，循环神经网络 $l$ 层时刻 $t$ 的反向传播的误差 $\Delta^{(l)}{t=T-i}$ 来自同一时刻上一层的误差 $\Delta^{(l+1)} {t=T-i}$ 和该时刻之后沿时间步反向传播的误差和 $\sum_{j=0}^{i-1} \Delta^{(l)}_{t=T-j}$ ， 当然这些误差传递过程中会有一些“损耗”也就是要乘上对应的距离函数。
下面重新给出矩阵形式反向传播算法的计算公式，
|计算公式|维度变换|
|-|-|-|
|$\Delta_{t=i}^L = A^L_{t=i} - Y_{t=i}$|$(m , d^L) = (m , d^L) - (m , d^L)$|
|$\Delta^{l}{t=T}=\Delta^{(l+1)} {t=T}\Psi(Z_{t=T}^{(l+1)}, Z_{t=T}^{(l)})$|$(m, d^{l})=(m, d^{l+1})(d^{l+1}, d^{l}) * (m, d^{l})$|
|$\Delta^{(l)}{t=T-i}=\Delta^{(l+1)} {t=T-i}\Psi(Z_{t=T-i}^{(l+1)}, Z_{t=T-i}^{(l)}) + \sum_{j=0}^{i-1} \Delta^{(l)}{t=T-j}\Omega(Z {T-j}^{(l)}, Z_{T-i}^{(l)}) $|$ (m, d^{l})=(m, d^{l+1})(d^{l+1}, d^{l}) * (m, d^{l}) + (m, d^l)(d^l, d^l)^{i-j} * (m, d^l)^{i-j}$|
|$\dfrac{\partial C}{\partial b^{l}} =\dfrac{1}{m}[{\sum_{t=i}^T \Delta_{t=i}^{(l)}}]^T $|$ (d^l, 1) = {(1, d^l)}^T$|
|$\dfrac{\partial C}{\partial W^{l}} = \dfrac{1}{m} \sum_{t=i}^T [\Delta {t}^{(l)}]^T A_t^{(l-1)}$|$(d^l, d^{(l-1)}) = {(m, d^l)}^T (m, d^{l-1})$|
|$\dfrac{\partial C}{\partial S^{l}} = \dfrac{1}{m} \sum {t=i}^T [\Delta_{t}^{(l)}]^T H_t^{(l)}$|$(d^l, d^{(l)}) = {(m, d^l)}^T (m, d^{l})$|
直接使用 RNN 梯度反向传播算法计算公式的时间复杂度 循环神经网络与前馈神经网络的关键区别就是，含有记忆向量的循环层，循环层虽然赋予了循环神经网络“记忆”能力，但在梯度反向传播带来了极大的复杂性。循环层的梯度公式如下，
$$\Delta^{(l)}{t=T-i}=\Delta^{(l+1)} {t=T-i}\Psi(Z_{t=T-i}^{(l+1)}, Z_{t=T-i}^{(l)}) + \sum_{j=0}^{i-1} \Delta^{(l)}{t=T-j}\Omega(Z {T-j}^{(l)}, Z_{T-i}^{(l)}) $$
上面式子中关键就是要求出距离函数 $\Omega$ 和 $\Psi$ 的值。
注意到 $\Omega(·)$ 距离函数中会重复的计算许多 $\sigma'(Z_{t})$ ，而 $Z_t$ 是 (小批量样本数，样本特征维度) 的大型矩阵，$\sigma'(·)$ 是激活函数的导数，是要对 $Z_t$ 矩阵中的每个元素使用该激活函数的导数，所以计算一次比较耗时。下面以一个实例来展示直接使用该公式计算的效果。
$\Delta^{l}{t=4}=\Delta^{(l+1)} {t=4}\Psi(Z_{t=4}^{(l+1)}, Z_{t=4}^{(l)})$
$\Delta^{l}{t=3}=\Delta^{(l+1)} {t=3}\Psi(Z_{t=3}^{(l+1)}, Z_{t=3}^{(l)}) + \Delta^{l}{t=4}[S^T] \sigma'(Z {t=3}^{(l)})$
$\Delta^{l}{t=2}=\Delta^{(l+1)} {t=2}\Psi(Z_{t=2}^{(l+1)}, Z_{t=2}^{(l)}) + \Delta^{l}{t=3}[S^T] \sigma'(Z {t=2}^{(l)}) + \Delta^{l}{t=4}[S^T]^2 \sigma'(Z {t=3}^{(l)})\sigma'(Z_{t=2}^{(l)})$
$\Delta^{l}{t=1}=\Delta^{(l+1)} {t=1}\Psi(Z_{t=1}^{(l+1)}, Z_{t=1}^{(l)}) + \Delta^{l}{t=2}[S^T] \sigma'(Z {t=1}^{(l)}) + \Delta^{l}{t=3}[S^T]^2 \sigma'(Z {t=2}^{(l)})\sigma'(Z_{t=1}^{(l)}) + \Delta^{l}{t=4}[S^T]^3 \sigma'(Z {t=3}^{(l)})\sigma'(Z_{t=2}^{(l)})\sigma'(Z_{t=1}^{(l)})$
$\Delta^{l}{t=0}=\Delta^{(l+1)} {t=0}\Psi(Z_{t=0}^{(l+1)}, Z_{t=0}^{(l)}) + \Delta^{l}{t=1}[S^T] \sigma'(Z {t=0}^{(l)}) + \Delta^{l}{t=2}[S^T]^2 \sigma'(Z {t=1}^{(l)})\sigma'(Z_{t=0}^{(l)}) + \Delta^{l}{t=3}[S^T]^3 \sigma'(Z {t=2}^{(l)})\sigma'(Z_{t=1}^{(l)})\sigma'(Z_{t=0}^{(l)}) + \Delta^{l}{t=4}[S^T]^4 \sigma'(Z {t=3}^{(l)})\sigma'(Z_{t=2}^{(l)})\sigma'(Z_{t=1}^{(l)})\sigma'(Z_{t=1}^{(0)})$
可以发现，随着有很多 $\sigma'(Z_{t=i}^{(l)})\sigma'(Z_{t=j}^{(l)})$ 乘法是在重复计算，而矩阵乘法操作是非常消耗时间的，加上 $\sigma'(Z_{t=i}^{(l)})$ 表示要用激活函数的导函数对矩阵 $Z$ 中的每一个元素进行计算，所以一次 $\sigma'(Z_{t=i}^{(l)})$ 矩阵相乘的计算量是非常大的。下面定量的计算直接使用 RNN 梯度反向传播算法计算公式的乘法次数。
假设 RNN 的循环时间步数为 $T(T\ge3)$ ，只统计占主要计算量的同一层节点梯度距离函数矩阵乘法的次数和记忆向量的权重矩阵 $S$ 的矩阵乘法次数。
先统计 $\sigma'(Z_i)$ 乘 $\sigma'(Z_j)$ 的矩阵乘法次数。
对于 $\Delta_{t=0}$ 矩阵乘法次数为 $1+2+3+···+(T-2)+(T-1)=\dfrac{T(T-1)}{2}$
对于 $\Delta_{t=1}$ 矩阵乘法次数为 $1+2+3+···+(T-3)+(T-2)=\dfrac{(T-1)(T-2)}{2}$
对于 $\Delta_{t=k}$ 矩阵乘法次数为 $1+2+3+···+(T-3)+(T-2)=\dfrac{(T-k)(T-(k+1))}{2}$
所以$\sigma'(Z_i)$ 乘 $\sigma'(Z_j)$ 的矩阵乘法次数为，
$$
\begin{aligned}
\sum_{t=2}^T \dfrac{t(t-1)}{2} &= \dfrac{1}{2}(\sum_{t=2}^T t^2-\sum_{t=2}^T t) \\
&= \dfrac{1}{2}(\dfrac{T(T+1)(2T+1)}{6}-\dfrac{T(T+1)}{2})\\
&= \dfrac{1}{6}(T^3-T) \\
\end{aligned}
$$
因为 $S$ 乘 $S$ 的矩阵乘法次数与$\sigma'(Z_i)$ 乘 $\sigma'(Z_j)$ 的矩阵乘法次数相同，所以RNN 的循环时间步数为 $T(T\ge3)$ 时循环层的距离函数矩阵乘法次数为 $\dfrac{1}{3}(T^3-T)$ ，时间复杂度为 $\mathcal{O}(T^3)$ 。
如上图所示，RNN 梯度反向传播主要就是要计算 $\Delta^l_{t=0},\Delta^l_{t=1}..\Delta^l_{t=T}$ ，但是梯度反向传播时前一时间步的梯度 $\Delta^l_{t=T-1}$ 依赖后一时间步 $\Delta^l_{t=T-1}$ ，而计算不同时间步 $\Delta^l_{t}$ 需要计算的 $\sigma'(Z_{t}^{(l)})$ 不同，也就是说 $\sigma'(Z_{t}^{(l)})$ 的计算根时间有关系。而动态规划擅长用来解决时间步依赖求最优化的问题，这也是本文选用动态规划算法的出发点。
$\sigma'(Z_{t=i}^{(l)})$ 表示要用激活函数的导函数对矩阵 $Z$ 中的每一个元素进行计算，而矩阵 $Z$ 的维度为 $(m, d^l)$ 含有 $m * d^l$ 个元素且激活函数的导函数一般很复杂（可能含有指数运行），所以该操作非常耗时，是本文主要优化的对象。
循环神经网络不同时间步同一层节点梯度距离函数,
$$\Omega(Z_{t_j}^{(l)}, Z_{t_i}^{(l)})=\dfrac{\partial Z_{t_j}^{(l)}}{\partial Z_{t_i}^{(l)}}=\prod_{t=t_i}^{t_{j}-1}(S^T)^{(l)} \odot \sigma'(Z_{t}^{(l)}) $$
定义 $\Omega$ 动态规划梯度矩阵（又称之为 $\Omega$ 梯度距离矩阵），用 $\Omega(j, i)$ 代表 $\Omega(Z_{t_j}^{(l)}, Z_{t_i}^{(l)})$ 。
$\Omega(j, i)$ $0$ $1$ $···$ $T-2$ $T-1$ $T$
$0$
$1$ $(1,0)$
$2$ $(2, 1)$
$···$ $(···, ···)$
$T-1$ $(T-1, T-2)$
$T$ $(T, T-1)$
首先，主对角线下面的第 $1$ 条对角线上的元素可以直接计算，坐标满足$(j-i=1)$的元素。
$$\Omega(T, T-1)=\Omega(Z_{T}^{(l)}, Z_{T-1}^{(l)})=\dfrac{\partial Z_{T}^{(l)}}{\partial Z_{T-1}^{(l)}}=\sigma'(Z_{T-1}^{(l)})$$
然后计算主对角线下面的第 $2$ 条对角线上的元素，坐标满足$(j-i=2)$的元素。
$\Omega(2, 0)=\Omega(2, 1) * \Omega(1, 0)$
$\Omega(3, 1)=\Omega(3, 2) * \Omega(2, 1)$
$\Omega(4, 2)=\Omega(4, 3) * \Omega(3, 2)$
然后计算主对角线下面的第 $k$ 条对角线上的元素，坐标满足$(j-i=k)$的元素。
计算使用的动态规划状态转移方程如下，
$$\Omega(j, i)=\Omega(j, i+1) * \Omega(i+1, i)$$
该方程表示的意义是，在 $\Omega$ 函数空间中 $j$ 和 $i$ 两点距离 $\Omega(j, i)$ 等于 $j$ 点
到 $i+1$ 点的距离 $\Omega(j, i+1)$ 乘以 $i+1$ 点到 $i$ 点的距离 $\Omega(i+1, i)$ 。该方程的其实就是链式求导法则在时间约束下的计算方程，证明和理解是容易的，如下，
方程左边等于 $\Omega(j, i)=\Omega(Z_{j}^{(l)}, Z_{i}^{(l)})=\dfrac{\partial Z_{j}^{(l)}}{\partial Z_{i}^{(l)}}$
方程右边等于 $\Omega(j, i+1) * \Omega(i+1, i) =\Omega(Z_{j}^{(l)}, Z_{i+1}^{(l)}) * \Omega(Z_{i+1}^{(l)}, Z_{i}^{(l)})=\dfrac{\partial Z_{j}^{(l)}}{\partial Z_{i+1}^{(l)}}*\dfrac{\partial Z_{i+1}^{(l)}}{\partial Z_{i}^{(l)}}=\dfrac{\partial Z_{j}^{(l)}}{\partial Z_{i}^{(l)}}$
方程左边等于右边，所以方程 $\Omega(j, i)=\Omega(j, i+1) * \Omega(i+1, i)$ 成立。
$\Omega$ 动态规划梯度矩阵计算顺序图如下，
RNN 循环层的递推梯度公式如下，
$$\Delta^{(l)}{t=T-i}=\Delta^{(l+1)} {t=T-i}\Psi(Z_{t=T-i}^{(l+1)}, Z_{t=T-i}^{(l)}) + \sum_{j=0}^{i-1} \Delta^{(l)}{t=T-j}\Omega(Z {T-j}^{(l)}, Z_{T-i}^{(l)}) $$
为了尽可能减少循环，加速运算，可以把上面介绍的 $\Omega$ 动态规划梯度矩阵再做增强，变成 $\Omega$ 动态规划梯度矩阵的增广矩阵 $\widetilde \Omega$ 。
增广矩阵 $\widetilde \Omega$ 的原来的第 $T$ 列（更名为 $\Delta^{(l)}$ 列）用来储存 RNN 第 $l$ 层各个时间步反向传播的误差，最后一列 $\Psi$ 用来储存 RNN 循环层的递推梯度公式第一项的计算结果。
与用 $\Omega(j, i)$ 代表 $\Omega(Z_{t_j}^{(l)}, Z_{t_i}^{(l)})$ 一样，这里用$\Psi_{t=k}$ 代表 $\Psi(Z_{t=k}^{(l+1)}, Z_{t=k}^{(l)})$ 。
$$\Psi_{t=k}=\Psi(Z_{t=k}^{(l+1)}, Z_{t=k}^{(l)})=\dfrac{\partial Z_{t=k}^{(l+1)}}{\partial Z_{t=k}^{(l)}}$$
$\widetilde \Omega(j, i)$ $0$ $1$ $···$ $T-1$ $\Delta^{(l)}$ $\Psi$
$0$ $\Delta^{(l)}_{t=0}$ $\Delta^{(l+1)}{t=0}\Psi {t=0}$
$1$ $(1,0)$ $\Delta^{(l)}_{t=1}$ $\Delta^{(l+1)}{t=1}\Psi {t=1}$
$2$ $(2, 0)$ $(2, 1)$ $\Delta^{(l)}_{t=2}$ $\Delta^{(l+1)}{t=2}\Psi {t=2}$
$···$ $(···, ···)$ $···$ $···$
$T$ $(T, 0)$ $(T, 1)$ $(T, T-1)$ $\Delta^{(l)}_{t=T}$ $\Delta^{(l+1)}{t=T}\Psi {t=T}$
增广矩阵 $\widetilde \Omega$ 的计算方法是，首先按照时间顺序 $(t=0, t=1···)$ 依次计算出最后一列 $\Psi$ 的值，然后用 $\Psi$ 列最后一个元素 $\Delta^{(l+1)}{t=T}\Psi {t=T}$ 初始化 $\Delta^{(l)}$ 列最后一个元素 $\Delta^{(l)}_{t=T}$ 。
然后 $i$ 从 $T-1$ 列到第 $0$ 列遍历该增广矩阵 $(i=T-1, i=T-2,...,i=1, i=0)$ ，动态规划的转移方程是，
$$\Delta^{(l)}[i]=\Psi[i]+(\Delta^{(l)})^T\widetilde \Omega[:,i]$$
该方程表示，$\Delta^{(l)}$ 列第 $i$ 个元素的值等于 $\widetilde \Omega$ 列第 $i$ 个元素的值加上 $\Delta^{(l)}$ 列的转置乘以 $\widetilde \Omega$ 矩阵第 $i$ 列的值。
该方程就是 RNN 循环层梯度向量分量计算式 $\Delta^{(l)}{t=T-i}=\Delta^{(l+1)} {t=T-i}\Psi(Z_{t=T-i}^{(l+1)}, Z_{t=T-i}^{(l)}) + \sum_{j=0}^{i-1} \Delta^{(l)}{t=T-j}\Omega(Z {T-j}^{(l)}, Z_{T-i}^{(l)}) $ 的矩阵形式表达。
$\widetilde \Omega$ 动态规划梯度矩阵计算顺序图如下，
使用动态规划算法求解 RNN 梯度反向传播的时间复杂度 与直接使用 RNN 梯度反向传播算法计算公式的时间复杂度假设一致，假设 RNN 的循环时间步数为 $T(T\ge3)$ ，只统计占主要计算量的同一层节点梯度距离函数矩阵乘法的次数和记忆向量的权重矩阵 $S$ 的矩阵乘法次数。
很明显用动态规划算法求解时，只需要计算 $\Omega$ 动态规划梯度矩阵的左下的小三角中元素的值，而每求一个元素的值，根据 $\Omega$ 动态规划梯度矩阵的状态转移方程 $\Omega(j, i)=\Omega(j, i+1) * \Omega(i+1, i)$ 可知求解每一个元素只需做一次矩阵乘法，所以一共需要 $\sigma'(Z_i)$ 乘 $\sigma'(Z_j)$ 的矩阵乘法次数为，
$$
\sum_{i=1}^T 1+2+···+(T-3)+(T-2) = \dfrac{(T-1)(T-2)}{2}
$$
因为 $S$ 乘 $S$ 的矩阵乘法次数与$\sigma'(Z_i)$ 乘 $\sigma'(Z_j)$ 的矩阵乘法次数相同，所以RNN 的循环时间步数为 $T(T\ge3)$ 时循环层的距离函数矩阵乘法次数为 $(T-1)(T-2)$ ，时间复杂度为 $\mathcal{O}(T^2)$ 。
其实，因为 $S$ 矩阵的值是固定的，可以先把它与自身的乘积 $S^2,S^3,···,S^{T-1}$ 计算出来，只需要 $T-2$ 次矩阵乘法，所以实际上使用动态规划算法加上这个优化只需要 $\dfrac{T(T-2)}{2}$ 次矩阵乘法运算，当然时间复杂度还是 $\mathcal{O}(T^2)$ 。
随着时间步数 $T$ 增加，循环神经网络越来越难训练。本文分析了计算循环神经网络反向传播算法中循环层的时间复杂度为 $\mathcal{O}(T^3)$ ， 为了降低计算复杂度，本文定义了两种梯度距离函数，然后由距离函数推导出 $\Omega$ 梯度距离矩阵，最后使用动态规划算法求解 $\Omega$ 梯度距离矩阵，使循环层的时间复杂度降低为 $\mathcal{O}(T^2)$ 。
循环神经网络不同时间步同一层节点梯度距离函数,
$$\Omega(Z_{t_j}^{(l)}, Z_{t_i}^{(l)})=\dfrac{\partial Z_{t_j}^{(l)}}{\partial Z_{t_i}^{(l)}}=\prod_{t=t_i}^{t_{j}-1}(S^T)^{(l)} \odot \sigma'(Z_{t}^{(l)}) $$
和同一时间步相邻层节点梯度距离函数,
$$\Psi(Z_{t_i}^{(l+1)}, Z_{t_i}^{(l)})=\dfrac{\partial Z_{t_i}^{(l+1)}}{\partial Z_{t_i}^{(l)}}=W^{(l+1)} \odot \sigma'(Z^{(l)}_{t_i})$$
然后使用动态规划算法优化了 $\Omega$ 的计算，主要用到了 $\Omega$ 梯度距离矩阵的转移方程，
$$\Omega(j, i)=\Omega(j, i+1) * \Omega(i+1, i)$$
接着提出 $\Omega$ 梯度距离矩阵的增广矩阵 $\widetilde \Omega$ ，使循环层的所有计算都矩阵化，进一步简化计算公式和加速计算速度，主要用到了 $\widetilde \Omega$ 矩阵的转移方程，
$$\Delta^{(l)}[i]=\Psi[i]+(\Delta^{(l)})^T\widetilde \Omega[:,i]$$
至此就可以计算出限制循环神经网络计算效率的循环层梯度反向传播速度的 $\Delta$ 。
根据以上研究提出全矩阵形式的循环神经网络反向传播计算公式，
计算公式
维度变化
$\Delta^L = A^L - Y$ $(T, m, d^L)=(T, m, d^L)-(T, m, d^L)$
$\Delta^{l}{t=T}=\Delta^{(l+1)} {t=T}\Psi(Z_{t=T}^{(l+1)}, Z_{t=T}^{(l)})$
$(m, d^{l})=(m, d^{l+1})(d^{l+1}, d^{l}) * (m, d^{l})$
$\Delta^{(l)}[i]=\Psi[i]+(\Delta^{(l)})^T\widetilde \Omega[:,i]$ $(m, d^{l})=(m, d^{l+1})(d^{l+1}, d^{l}) * (m, d^{l}) + (m, d^l)(d^l, d^l)^{i-j} * (m, d^l)^{i-j}$
$\dfrac{\partial C}{\partial b^{l}} =\dfrac{1}{m}({sum(\Delta^{(l)}},axis=0,1))^T $ $(d^l, 1) = {(sum(T,m, d^l)),axis=0,1)}^T$
$\dfrac{\partial C}{\partial W^{l}} = \dfrac{1}{m} sum((\Delta^{(l)}).swapaxes(1,2) A^{(l-1)},axis=0)$ $(d^l,d^{l-1})=sum((T,d^l,m)(T,m,d^{l-1}),axis=0)$
$\dfrac{\partial C}{\partial S^{l}} = \dfrac{1}{m} sum((\Delta^{(l)}).swapaxes(1,2) H^{(l)},axis=0)$ $(d^l,d^{l})=sum((T,d^l,m)(T,m,d^{l}),axis=0)$
说明：
${sum(\Delta^{(l)}},axis=0,1))^T$ 表示对维度是 $(T, m, d^l)$ 的三维张量 $\Delta^{(l)}$ 的第一和第二维度上的元素进行求和操作，最后得到一个维度是 $(1,d^l)$ 向量再进行转置运算；
$sum((\Delta^{(l)}).swapaxes(1,2) A^{(l-1)},axis=0)$ 表示对维度是 $(T, m, d^l)$ 的三维张量 $\Delta^{(l)}$ 的后面两个维度表示的矩阵 $\Delta^{(l)}[i]$ 进行转置变成 $(T,d^l,m)$ ，再与维度是 $(T, m, d^l)$ 的三维张量 $A^{(l-1)}$ 后两个维度表示的矩阵 $A^{(l-1)}[i]$ 进行矩阵乘积，结果为 $(T, d^l,d^{l-1})$ ，最后在对该结果的第一维度上的元素进行求和操作，得到矩阵 $(d^l,d^{l-1})$ 。
