ctc/CTC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CTC（ Connectionist Temporal Classification，连接时序分类）是一种用于序列建模的工具，其核心是定义了特殊的**目标函数/优化准则**[1]。\n",
    "\n",
    "> jupyter notebook 版见 [repo](https://github.com/DingKe/ml-tutorial/tree/master/ctc)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# 1. 算法\n",
    "这里大体根据 Alex Grave 的开山之作[1]，讨论 CTC 的算法原理，并基于 numpy 从零实现 CTC 的推理及训练算法。\n",
    "\n",
    "## 1.1 序列问题形式化。\n",
    "序列问题可以形式化为如下函数：\n",
    "\n",
    "$$N_w: (R^m)^T \\rightarrow (R^n)^T$$\n",
    "其中，序列目标为字符串（词表大小为 $n$），即 $N_w$ 输出为 $n$ 维多项概率分布（e.g. 经过 softmax 处理）。\n",
    "\n",
    "网络输出为：$y = N_w(x)$，其中，$y_k^t$ $t$ 表示时刻第 $k$ 项的概率。\n",
    "\n",
    "![](https://distill.pub/2017/ctc/assets/full_collapse_from_audio.svg)\n",
    "**图1. 序列建模【[src](https://distill.pub/2017/ctc/)】**\n",
    "\n",
    "\n",
    "虽然并没为限定 $N_w$ 具体形式，下面为假设其了某种神经网络（e.g. RNN）。\n",
    "下面代码示例 toy $N_w$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.24654511  0.18837589  0.16937668  0.16757465  0.22812766]\n",
      " [ 0.25443629  0.14992236  0.22945293  0.17240658  0.19378184]\n",
      " [ 0.24134404  0.17179604  0.23572466  0.12994237  0.22119288]\n",
      " [ 0.27216255  0.13054313  0.2679252   0.14184499  0.18752413]\n",
      " [ 0.32558002  0.13485564  0.25228604  0.09743785  0.18984045]\n",
      " [ 0.23855586  0.14800386  0.23100255  0.17158135  0.21085638]\n",
      " [ 0.38534786  0.11524603  0.18220093  0.14617864  0.17102655]\n",
      " [ 0.21867406  0.18511892  0.21305488  0.16472572  0.21842642]\n",
      " [ 0.29856607  0.13646801  0.27196606  0.11562552  0.17737434]\n",
      " [ 0.242347    0.14102063  0.21716951  0.2355229   0.16393996]\n",
      " [ 0.26597326  0.10009752  0.23362892  0.24560198  0.15469832]\n",
      " [ 0.23337289  0.11918746  0.28540761  0.20197928  0.16005275]]\n",
      "[[ 1.]\n",
      " [ 1.]\n",
      " [ 1.]\n",
      " [ 1.]\n",
      " [ 1.]\n",
      " [ 1.]\n",
      " [ 1.]\n",
      " [ 1.]\n",
      " [ 1.]\n",
      " [ 1.]\n",
      " [ 1.]\n",
      " [ 1.]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "np.random.seed(1111)\n",
    "\n",
    "T, V = 12, 5\n",
    "m, n = 6, V\n",
    "\n",
    "x = np.random.random([T, m])  # T x m\n",
    "w = np.random.random([m, n])  # weights, m x n\n",
    "\n",
    "def softmax(logits):\n",
    "    max_value = np.max(logits, axis=1, keepdims=True)\n",
    "    exp = np.exp(logits - max_value)\n",
    "    exp_sum = np.sum(exp, axis=1, keepdims=True)\n",
    "    dist = exp / exp_sum\n",
    "    return dist\n",
    "\n",
    "def toy_nw(x):\n",
    "    y = np.matmul(x, w)  # T x n \n",
    "    y = softmax(y)\n",
    "    return y\n",
    "\n",
    "y = toy_nw(x)\n",
    "print(y)\n",
    "print(y.sum(1, keepdims=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.2 align-free 变长映射\n",
    "上面的形式是输入和输出的一对一的映射。序列学习任务一般而言是多对多的映射关系（如语音识别中，上百帧输出可能仅对应若干音节或字符，并且每个输入和输出之间，也没有清楚的对应关系）。CTC 通过引入一个特殊的 blank 字符（用 % 表示），解决多对一映射问题。\n",
    "\n",
    "扩展原始词表 $L$ 为 $L^\\prime = L \\cup \\{\\text{blank}\\}$。对输出字符串，定义操作 $B$：1）合并连续的相同符号；2）去掉 blank 字符。\n",
    "\n",
    "例如，对于 “aa%bb%%cc”，应用 $B$，则实际上代表的是字符串 \"abc\"。同理“%a%b%cc%” 也同样代表 \"abc\"。\n",
    "$$\n",
    "B(aa\\%bb\\%\\%cc) = B(\\%a\\%b\\%cc\\%) = abc\n",
    "$$\n",
    "\n",
    "通过引入blank 及 $B$，可以实现了变长的映射。\n",
    "$$\n",
    "L^{\\prime T} \\rightarrow L^{\\le T}\n",
    "$$\n",
    "\n",
    "\t因为这个原因，CTC 只能建模输出长度小于输入长度的序列问题。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.3 似然计算\n",
    "和大多数有监督学习一样，CTC 使用最大似然标准进行训练。\n",
    "\n",
    "给定输入 $x$，输出 $l$ 的条件概率为：\n",
    "$$\n",
    "p(l|x) =  \\sum_{\\pi \\in B^{-1}(l)} p(\\pi|x)\n",
    "$$\n",
    "\n",
    "其中，$B^{-1}(l)$ 表示了长度为 $T$ 且示经过 $B$ 结果为 $l$ 字符串的集合。\n",
    "\n",
    "**CTC 假设输出的概率是（相对于输入）条件独立的**，因此有：\n",
    "$$p(\\pi|x) = \\prod y^t_{\\pi_t}, \\forall \\pi \\in L^{\\prime T}$$\n",
    "\n",
    "\n",
    "然而，直接按上式我们没有办理有效的计算似然值。下面用动态规划解决似然的计算及梯度计算, 涉及前向算法和后向算法。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.4 前向算法\n",
    "\n",
    "在前向及后向计算中，CTC 需要将输出字符串进行扩展。具体的，$(a_1,\\cdots,a_m)$ 每个字符之间及首尾分别插入 blank，即扩展为 $(\\%, a_1,\\%,a_2, \\%,\\cdots,\\%, a_m,\\%)$。下面的 $l$ 为原始字符串，$l^\\prime$ 指为扩展后的字符串。\n",
    "\n",
    "定义\n",
    "$$\n",
    "\\alpha_t(s) \\stackrel{def}{=} \\sum_{\\pi \\in N^T: B(\\pi_{1:t}) = l_{1:s}} \\prod_{t^\\prime=1}^t y^t_{\\pi^\\prime}\n",
    "$$\n",
    "\n",
    "显然有，\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\alpha_1(1) = y_b^1,\\\\\n",
    "\\alpha_1(2) = y_{l_1}^1,\\\\\n",
    "\\alpha_1(s) = 0, \\forall s > 2\n",
    "\\end{align}\n",
    "$$\n",
    "根据 $\\alpha$ 的定义，有如下递归关系：\n",
    "$$\n",
    "\\alpha_t(s) = \\{  \\begin{array}{l}\n",
    "(\\alpha_{t-1}(s)+\\alpha_{t-1}(s-1)) y^t_{l^\\prime_s},\\  \\  \\    if\\  l^\\prime_s = b \\ or\\  l_{s-2}^\\prime = l_s^{\\prime}  \\\\\n",
    "(\\alpha_{t-1}(s)+\\alpha_{t-1}(s-1) + \\alpha_{t-1}(s-2)) y^t_{l^\\prime_s} \\ \\ otherwise\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4.1 Case 2\n",
    "递归公式中 case 2 是一般的情形。如图所示，$t$ 时刻字符为 $s$ 为 blank 时，它可能由于两种情况扩展而来：1）重复上一字符，即上个字符也是 a，2）字符发生转换，即上个字符是非 a 的字符。第二种情况又分为两种情形，2.1）上一字符是 blank；2.2）a 由非 blank 字符直接跳转而来（$B$） 操作中， blank 最终会被去掉，因此 blank 并不是必须的）。\n",
    "![](https://distill.pub/2017/ctc/assets/cost_regular.svg)\n",
    "**图2. 前向算法 Case 2 示例【[src](https://distill.pub/2017/ctc/)】**\n",
    "\n",
    "### 1.4.2 Case 1\n",
    "递归公式 case 1 是特殊的情形。\n",
    "如图所示，$t$ 时刻字符为 $s$ 为 blank 时，它只能由于两种情况扩展而来：1）重复上一字符，即上个字符也是 blank，2）字符发生转换，即上个字符是非 blank 字符。$t$ 时刻字符为 $s$ 为非 blank 时，类似于 case 2，但是这时两个相同字符之间的 blank 不能省略（否则无法区分\"aa\"和\"a\"），因此，也只有两种跳转情况。\n",
    "\n",
    "![](https://distill.pub/2017/ctc/assets/cost_no_skip.svg)\n",
    "**图3. 前向算法 Case 1 【[src](https://distill.pub/2017/ctc/)】**\n",
    "\n",
    "我们可以利用动态规划计算所有 $\\alpha$ 的值，算法时间和空间复杂度为 $O(T * L)$。\n",
    "\n",
    "似然的计算只涉及乘加运算，因此，CTC 的似然是可导的，可以尝试 tensorflow 或 pytorch 等具有自动求导功能的工具自动进行梯度计算。下面介绍如何手动高效的计算梯度。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  2.46545113e-01   1.67574654e-01   0.00000000e+00   0.00000000e+00\n",
      "    0.00000000e+00   0.00000000e+00   0.00000000e+00]\n",
      " [  6.27300235e-02   7.13969720e-02   4.26370730e-02   0.00000000e+00\n",
      "    0.00000000e+00   0.00000000e+00   0.00000000e+00]\n",
      " [  1.51395174e-02   1.74287803e-02   2.75214373e-02   5.54036251e-03\n",
      "    0.00000000e+00   0.00000000e+00   0.00000000e+00]\n",
      " [  4.12040964e-03   4.61964998e-03   1.22337658e-02   4.68965079e-03\n",
      "    1.50787918e-03   1.03895167e-03   0.00000000e+00]\n",
      " [  1.34152305e-03   8.51612635e-04   5.48713543e-03   1.64898136e-03\n",
      "    2.01779193e-03   1.37377693e-03   3.38261905e-04]\n",
      " [  3.20028190e-04   3.76301179e-04   1.51214552e-03   1.22442454e-03\n",
      "    8.74730268e-04   1.06283215e-03   4.08416903e-04]\n",
      " [  1.23322177e-04   1.01788478e-04   7.27708889e-04   4.00028082e-04\n",
      "    8.08904808e-04   5.40783712e-04   5.66942671e-04]\n",
      " [  2.69673617e-05   3.70815141e-05   1.81389560e-04   1.85767281e-04\n",
      "    2.64362267e-04   3.82184328e-04   2.42231029e-04]\n",
      " [  8.05153930e-06   7.40568461e-06   6.52280509e-05   4.24527009e-05\n",
      "    1.34393412e-04   1.47631121e-04   1.86429242e-04]\n",
      " [  1.95126637e-06   3.64053019e-06   1.76025677e-05   2.53612828e-05\n",
      "    4.28581244e-05   5.31947855e-05   8.09585256e-05]\n",
      " [  5.18984675e-07   1.37335633e-06   5.65009596e-06   1.05520069e-05\n",
      "    1.81445380e-05   1.87825719e-05   3.56811933e-05]\n",
      " [  1.21116956e-07   3.82213679e-07   1.63908339e-06   3.27248912e-06\n",
      "    6.69699576e-06   7.59916314e-06   1.27103665e-05]]\n"
     ]
    }
   ],
   "source": [
    "def forward(y, labels):\n",
    "    T, V = y.shape\n",
    "    L = len(labels)\n",
    "    alpha = np.zeros([T, L])\n",
    "\n",
    "    # init\n",
    "    alpha[0, 0] = y[0, labels[0]]\n",
    "    alpha[0, 1] = y[0, labels[1]]\n",
    "\n",
    "    for t in range(1, T):\n",
    "        for i in range(L):\n",
    "            s = labels[i]\n",
    "            \n",
    "            a = alpha[t - 1, i] \n",
    "            if i - 1 >= 0:\n",
    "                a += alpha[t - 1, i - 1]\n",
    "            if i - 2 >= 0 and s != 0 and s != labels[i - 2]:\n",
    "                a += alpha[t - 1, i - 2]\n",
    "                \n",
    "            alpha[t, i] = a * y[t, s]\n",
    "            \n",
    "    return alpha\n",
    "\n",
    "labels = [0, 3, 0, 3, 0, 4, 0]  # 0 for blank\n",
    "alpha = forward(y, labels)\n",
    "print(alpha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最后可以得到似然 $p(l|x) = \\alpha_T(|l^\\prime|) + \\alpha_T(|l^\\prime|-1)$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.81811271177e-06\n"
     ]
    }
   ],
   "source": [
    "p = alpha[-1, labels[-1]] + alpha[-1, labels[-2]]\n",
    "print(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.5 后向计算\n",
    "类似于前向计算，我们定义后向计算。\n",
    "首先定义\n",
    "$$\n",
    "\\beta_t(s)   \\stackrel{def}{=} \\sum_{\\pi \\in N^T: B(\\pi_{t:T}) = l_{s:|l|}} \\prod_{t^\\prime=t}^T y^t_{\\pi^\\prime}\n",
    "$$\n",
    "\n",
    "显然，\n",
    "$$\n",
    "\\begin{align}\n",
    "\\beta_T(|l^\\prime|) = y_b^T,\\\\\n",
    "\\beta_T(|l^\\prime|-1) = y_{l_{|l|}}^T,\\\\\n",
    "\\beta_T(s) = 0, \\forall s < |l^\\prime| - 1\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "易得如下递归关系：\n",
    "$$\n",
    "\\beta_t(s) = \\{  \\begin{array}{l}\n",
    "(\\beta_{t+1}(s)+\\beta_{t+1}(s+1)) y^t_{l^\\prime_s},\\  \\  \\    if\\  l^\\prime_s = b \\ or\\  l_{s+2}^\\prime = l_s^{\\prime}  \\\\\n",
    "(\\beta_{t+1}(s)+\\beta_{t+1}(s+1) + \\beta_{t+1}(s+2)) y^t_{l^\\prime_s} \n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  1.25636660e-05   7.74586366e-06   8.69559539e-06   3.30990037e-06\n",
      "    2.41325357e-06   4.30516936e-07   1.21116956e-07]\n",
      " [  3.00418145e-05   2.09170784e-05   2.53062822e-05   9.96351200e-06\n",
      "    8.39236521e-06   1.39591874e-06   4.91256769e-07]\n",
      " [  7.14014755e-05   4.66705755e-05   7.46535563e-05   2.48066359e-05\n",
      "    2.77113594e-05   5.27279259e-06   1.93076535e-06]\n",
      " [  1.69926001e-04   1.25923340e-04   2.33240296e-04   7.60839197e-05\n",
      "    9.89830489e-05   1.58379311e-05   8.00005392e-06]\n",
      " [  4.20893778e-04   2.03461048e-04   6.84292101e-04   1.72696845e-04\n",
      "    3.08627225e-04   5.50636993e-05   2.93943967e-05]\n",
      " [  4.81953899e-04   8.10796738e-04   1.27731424e-03   8.24448952e-04\n",
      "    7.48161143e-04   1.99769340e-04   9.02831714e-05]\n",
      " [  9.80428697e-04   1.03986915e-03   3.68556718e-03   1.66879393e-03\n",
      "    2.56724754e-03   5.68961868e-04   3.78457146e-04]\n",
      " [  2.40870506e-04   2.30339872e-03   4.81028886e-03   4.75397134e-03\n",
      "    4.31752827e-03   2.34462771e-03   9.82118206e-04]\n",
      " [  0.00000000e+00   1.10150469e-03   1.28817322e-02   9.11579592e-03\n",
      "    1.35011919e-02   6.24293419e-03   4.49124231e-03]\n",
      " [  0.00000000e+00   0.00000000e+00   9.52648414e-03   3.36188472e-02\n",
      "    2.50664437e-02   2.01536701e-02   1.50427081e-02]\n",
      " [  0.00000000e+00   0.00000000e+00   0.00000000e+00   3.93092725e-02\n",
      "    4.25697510e-02   6.08622868e-02   6.20709492e-02]\n",
      " [  0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00\n",
      "    0.00000000e+00   1.60052748e-01   2.33372894e-01]]\n"
     ]
    }
   ],
   "source": [
    "def backward(y, labels):\n",
    "    T, V = y.shape\n",
    "    L = len(labels)\n",
    "    beta = np.zeros([T, L])\n",
    "\n",
    "    # init\n",
    "    beta[-1, -1] = y[-1, labels[-1]]\n",
    "    beta[-1, -2] = y[-1, labels[-2]]\n",
    "\n",
    "    for t in range(T - 2, -1, -1):\n",
    "        for i in range(L):\n",
    "            s = labels[i]\n",
    "            \n",
    "            a = beta[t + 1, i] \n",
    "            if i + 1 < L:\n",
    "                a += beta[t + 1, i + 1]\n",
    "            if i + 2 < L and s != 0 and s != labels[i + 2]:\n",
    "                a += beta[t + 1, i + 2]\n",
    "                \n",
    "            beta[t, i] = a * y[t, s]\n",
    "            \n",
    "    return beta\n",
    "\n",
    "beta = backward(y, labels)\n",
    "print(beta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.6 梯度计算\n",
    "下面，我们利用前向、后者计算的 $\\alpha$ 和 $\\beta$ 来计算梯度。\n",
    "\n",
    "根据 $\\alpha$、$\\beta$ 的定义，我们有：\n",
    "$$\\alpha_t(s)\\beta_t(s) = \\sum_{\\pi \\in  B^{-1}(l):\\pi_t=l_s^\\prime} y^t_{l_s^\\prime} \\prod_{t=1}^T y^t_{\\pi_t} = y^t_{l_s^\\prime} \\cdot \\sum_{\\pi \\in  B^{-1}(l):\\pi_t=l_s^\\prime} \\prod_{t=1}^T y^t_{\\pi_t}$$\n",
    "则\n",
    "$$\n",
    "\\frac{\\alpha_t(s)\\beta_t(s)}{ y^t_{l_s^\\prime}} = \\sum_{\\pi \\in  B^{-1}(l):\\pi_t=l_s^\\prime} \\prod_{t=1}^T y^t_{\\pi_t}  = \\sum_{\\pi \\in  B^{-1}(l):\\pi_t=l_s^\\prime} p(\\pi|x) \n",
    "$$\n",
    "于是，可得似然\n",
    "$$\n",
    "p(l|x) = \\sum_{s=1}^{|l^\\prime|} \\sum_{\\pi \\in  B^{-1}(l):\\pi_t=l_s^\\prime} p(\\pi|x) =  \\sum_{s=1}^{|l^\\prime|}  \\frac{\\alpha_t(s)\\beta_t(s)}{ y^t_{l_s^\\prime}} \n",
    "$$\n",
    "\n",
    "\n",
    "为计算 $\\frac{\\partial p(l|x)}{\\partial y^t_k}$，观察上式右端求各项，仅有 $s=k$ 的项包含 $y^t_k$，因此，其他项的偏导都为零，不用考虑。于是有：\n",
    "\n",
    "$$\n",
    "\\frac{\\partial p(l|x)}{\\partial y^t_k} = \\frac{\\partial \\frac{\\alpha_t(k)\\beta_t(k)}{ y^t_{k}} }{\\partial y^t_k} \n",
    "$$\n",
    "\n",
    "利用除法的求导准则有：\n",
    "$$\n",
    "\\frac{\\partial p(l|x)}{\\partial y^t_k} =  \\frac{\\frac{2 \\cdot \\alpha_t(k)\\beta_t(k)}{ y^t_{k}} \\cdot y^t_{k} -  \\alpha_t(k)\\beta_t(k) \\cdot 1}{{y^t_k}^2} = \\frac{\\alpha_t(k)\\beta_t(k)}{{y^t_k}^2}\n",
    "$$\n",
    "\n",
    "> 求导中，分子第一项是因为 $\\alpha(k)\\beta(k)$ 中包含为两个 $y^t_k$ 乘积项（即 ${y^t_k}^2$），其他均为与 $y^t_k$ 无关的常数。\n",
    "\n",
    "$l$ 中可能包含多个 $k$ 字符，它们计算的梯度要进行累加，因此，最后的梯度计算结果为：\n",
    "$$\n",
    "\\frac{\\partial p(l|x)}{\\partial y^t_k} = \\frac{1}{{y^t_k}^2} \\sum_{s \\in lab(l, k)} \\alpha_t(s)\\beta_t(s)\n",
    "$$\n",
    "其中，$lab(s)=\\{s: l_s^\\prime = k\\}$。\n",
    "\n",
    "一般我们优化似然函数的对数，因此，梯度计算如下：\n",
    "$$\n",
    "\\frac{\\partial \\ln(p(l|x))}{\\partial y^t_k} =\\frac{1}{p(l|x)} \\frac{\\partial p(l|x)}{\\partial y^t_k}\n",
    "$$\n",
    "其中，似然值在前向计算中已经求得： $p(l|x) = \\alpha_T(|l^\\prime|) + \\alpha_T(|l^\\prime|-1)$。\n",
    "\n",
    "对于给定训练集 $D$，待优化的目标函数为：\n",
    "$$\n",
    "O(D, N_w) = -\\sum_{(x,z)\\in D} \\ln(p(z|x))\n",
    "$$\n",
    "\n",
    "得到梯度后，我们可以利用任意优化方法（e.g. SGD， Adam）进行训练。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 2.50911241  0.          0.          2.27594441  0.        ]\n",
      " [ 2.25397118  0.          0.          2.47384957  0.        ]\n",
      " [ 2.65058465  0.          0.          2.77274592  0.        ]\n",
      " [ 2.46136916  0.          0.          2.29678159  0.02303985]\n",
      " [ 2.300259    0.          0.          2.37548238  0.10334851]\n",
      " [ 2.40271071  0.          0.          2.19860276  0.23513657]\n",
      " [ 1.68914157  0.          0.          1.78214377  0.51794046]\n",
      " [ 2.32536762  0.          0.          1.75750877  0.92477606]\n",
      " [ 1.92883907  0.          0.          1.45529832  1.44239844]\n",
      " [ 2.06219335  0.          0.          0.7568118   1.96405515]\n",
      " [ 2.07914466  0.          0.          0.33858403  2.35197258]\n",
      " [ 2.6816852   0.          0.          0.          2.3377753 ]]\n"
     ]
    }
   ],
   "source": [
    "def gradient(y, labels):\n",
    "    T, V = y.shape\n",
    "    L = len(labels)\n",
    "    \n",
    "    alpha = forward(y, labels)\n",
    "    beta = backward(y, labels)\n",
    "    p = alpha[-1, -1] + alpha[-1, -2]\n",
    "    \n",
    "    grad = np.zeros([T, V])\n",
    "    for t in range(T):\n",
    "        for s in range(V):\n",
    "            lab = [i for i, c in enumerate(labels) if c == s]\n",
    "            for i in lab:\n",
    "                grad[t, s] += alpha[t, i] * beta[t, i] \n",
    "            grad[t, s] /= y[t, s] ** 2\n",
    "                \n",
    "    grad /= p\n",
    "    return grad\n",
    "    \n",
    "grad = gradient(y, labels)\n",
    "print(grad)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "将基于前向-后向算法得到梯度与基于数值的梯度比较，以验证实现的正确性。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1e-05\n",
      "1e-06\n",
      "[0, 3]：3.91e-06\n",
      "[1, 0]：3.61e-06\n",
      "[1, 3]：2.66e-06\n",
      "[2, 0]：2.67e-06\n",
      "[2, 3]：3.88e-06\n",
      "[3, 0]：4.71e-06\n",
      "[3, 3]：3.39e-06\n",
      "[4, 0]：1.24e-06\n",
      "[4, 3]：4.79e-06\n",
      "[5, 0]：1.57e-06\n",
      "[5, 3]：2.98e-06\n",
      "[6, 0]：5.03e-06\n",
      "[6, 3]：4.89e-06\n",
      "[7, 0]：1.05e-06\n",
      "[7, 4]：4.19e-06\n",
      "[8, 4]：5.57e-06\n",
      "[9, 0]：5.95e-06\n",
      "[9, 3]：3.85e-06\n",
      "[10, 0]：1.09e-06\n",
      "[10, 3]：1.53e-06\n",
      "[10, 4]：3.82e-06\n"
     ]
    }
   ],
   "source": [
    "def check_grad(y, labels, w=-1, v=-1, toleration=1e-3):\n",
    "    grad_1 = gradient(y, labels)[w, v]\n",
    "    \n",
    "    delta = 1e-10\n",
    "    original = y[w, v]\n",
    "    \n",
    "    y[w, v] = original + delta\n",
    "    alpha = forward(y, labels)\n",
    "    log_p1 = np.log(alpha[-1, -1] + alpha[-1, -2])\n",
    "    \n",
    "    y[w, v] = original - delta\n",
    "    alpha = forward(y, labels)\n",
    "    log_p2 = np.log(alpha[-1, -1] + alpha[-1, -2])\n",
    "    \n",
    "    y[w, v] = original\n",
    "    \n",
    "    grad_2 = (log_p1 - log_p2) / (2 * delta)\n",
    "    if np.abs(grad_1 - grad_2) > toleration:\n",
    "        print('[%d, %d]：%.2e' % (w, v, np.abs(grad_1 - grad_2)))\n",
    "\n",
    "for toleration in [1e-5, 1e-6]:\n",
    "    print('%.e' % toleration)\n",
    "    for w in range(y.shape[0]):\n",
    "        for v in range(y.shape[1]):\n",
    "            check_grad(y, labels, w, v, toleration)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以看到，前向-后向及数值梯度两种方法计算的梯度差异都在 1e-5 以下，误差最多在 1e-6 的量级。这初步验证了前向-后向梯度计算方法原理和实现的正确性。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.7 logits 梯度\n",
    "在实际训练中，为了计算方便，可以将 CTC 和 softmax 的梯度计算合并，公式如下：\n",
    "$$\n",
    "\\frac{\\partial \\ln(p(l|x))}{\\partial y^t_k} = y^t_k - \\frac{1}{y^t_k \\cdot p(l|x)} \\sum_{s \\in lab(l, k)} \\alpha_t(s)\\beta_t(s)\n",
    "$$\n",
    "\n",
    "这是因为，softmax 的梯度反传公式为：\n",
    "\n",
    "$$\n",
    "\\frac{\\partial \\ln(p(l|x))}{\\partial u^t_k} = y^t_k (\\frac{\\partial \\ln(p(l|x))}{\\partial y^t_k}  - \\sum_{j=1}^{V} \\frac{\\partial \\ln(p(l|x))}{\\partial y^t_j} y^t_j)\n",
    "$$\n",
    "\n",
    "接合上面两式，有:\n",
    "$$\n",
    "\\frac{\\partial \\ln(p(l|x))}{\\partial u^t_k} = \\frac{1}{y^t_k p(l|x)} \\sum_{s \\in lab(l, k)} \\alpha_t(s)\\beta_t(s) - y^t_k\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.59941504485e-15\n"
     ]
    }
   ],
   "source": [
    "def gradient_logits_naive(y, labels):\n",
    "    '''\n",
    "    gradient by back propagation\n",
    "    '''\n",
    "    y_grad = gradient(y, labels)\n",
    "    \n",
    "    sum_y_grad = np.sum(y_grad * y, axis=1, keepdims=True)\n",
    "    u_grad = y * (y_grad - sum_y_grad) \n",
    "    \n",
    "    return u_grad\n",
    "\n",
    "def gradient_logits(y, labels):\n",
    "    '''\n",
    "    '''\n",
    "    T, V = y.shape\n",
    "    L = len(labels)\n",
    "    \n",
    "    alpha = forward(y, labels)\n",
    "    beta = backward(y, labels)\n",
    "    p = alpha[-1, -1] + alpha[-1, -2]\n",
    "    \n",
    "    u_grad = np.zeros([T, V])\n",
    "    for t in range(T):\n",
    "        for s in range(V):\n",
    "            lab = [i for i, c in enumerate(labels) if c == s]\n",
    "            for i in lab:\n",
    "                u_grad[t, s] += alpha[t, i] * beta[t, i] \n",
    "            u_grad[t, s] /= y[t, s] * p\n",
    "                \n",
    "    u_grad -= y\n",
    "    return u_grad\n",
    "    \n",
    "grad_l = gradient_logits_naive(y, labels)\n",
    "grad_2 = gradient_logits(y, labels)\n",
    "\n",
    "print(np.sum(np.abs(grad_l - grad_2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "同上，我们利用数值梯度来初步检验梯度计算的正确性："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1e-05\n",
      "[0, 6]：-8.00e-02, -8.00e-02, 1.03e-05\n",
      "[4, 0]：4.29e-01, 4.29e-01, 1.10e-05\n",
      "[5, 6]：-7.59e-02, -7.59e-02, 1.22e-05\n",
      "[6, 2]：-1.38e-01, -1.38e-01, 1.23e-05\n",
      "[6, 3]：3.33e-01, 3.33e-01, 1.02e-05\n",
      "1e-06\n",
      "[0, 0]：3.88e-01, 3.88e-01, 7.03e-06\n",
      "[0, 1]：-1.59e-01, -1.59e-01, 2.78e-06\n",
      "[0, 2]：-8.89e-02, -8.89e-02, 3.47e-06\n",
      "[0, 3]：4.57e-01, 4.57e-01, 1.64e-06\n",
      "[0, 4]：-6.32e-02, -6.32e-02, 7.19e-06\n",
      "[0, 5]：-7.98e-02, -7.98e-02, 5.46e-06\n",
      "[0, 6]：-8.00e-02, -8.00e-02, 1.03e-05\n",
      "[0, 7]：-1.32e-01, -1.32e-01, 2.21e-06\n",
      "[0, 8]：-1.04e-01, -1.04e-01, 7.75e-06\n",
      "[0, 9]：-1.38e-01, -1.38e-01, 5.95e-06\n",
      "[1, 0]：3.41e-01, 3.41e-01, 2.79e-06\n",
      "[1, 1]：-1.18e-01, -1.18e-01, 6.08e-06\n",
      "[1, 3]：5.04e-01, 5.04e-01, 4.06e-06\n",
      "[1, 4]：-9.96e-02, -9.96e-02, 5.77e-06\n",
      "[1, 5]：-8.22e-02, -8.22e-02, 4.03e-06\n",
      "[1, 6]：-9.46e-02, -9.46e-02, 4.49e-06\n",
      "[1, 7]：-1.49e-01, -1.49e-01, 3.96e-06\n",
      "[1, 8]：-1.24e-01, -1.24e-01, 4.96e-06\n",
      "[1, 9]：-7.48e-02, -7.47e-02, 5.94e-06\n",
      "[2, 0]：3.29e-01, 3.29e-01, 3.47e-06\n",
      "[2, 1]：-9.42e-02, -9.42e-02, 1.63e-06\n",
      "[2, 2]：-9.17e-02, -9.17e-02, 4.47e-06\n",
      "[2, 3]：4.50e-01, 4.50e-01, 2.14e-06\n",
      "[2, 5]：-1.07e-01, -1.07e-01, 6.33e-06\n",
      "[2, 6]：-5.42e-02, -5.42e-02, 1.71e-06\n",
      "[2, 7]：-9.68e-02, -9.68e-02, 7.69e-06\n",
      "[2, 9]：-1.21e-01, -1.21e-01, 9.06e-06\n",
      "[3, 0]：4.42e-01, 4.42e-01, 9.21e-06\n",
      "[3, 1]：-6.71e-02, -6.71e-02, 5.75e-06\n",
      "[3, 2]：-1.16e-01, -1.16e-01, 5.26e-06\n",
      "[3, 3]：4.03e-01, 4.03e-01, 6.39e-06\n",
      "[3, 4]：-1.07e-01, -1.07e-01, 2.42e-06\n",
      "[3, 5]：-1.25e-01, -1.25e-01, 8.90e-06\n",
      "[3, 6]：-1.17e-01, -1.17e-01, 2.08e-06\n",
      "[3, 7]：-1.32e-01, -1.32e-01, 2.21e-06\n",
      "[3, 8]：-6.90e-02, -6.90e-02, 1.72e-06\n",
      "[3, 9]：-1.13e-01, -1.13e-01, 6.68e-06\n",
      "[4, 0]：4.29e-01, 4.29e-01, 1.10e-05\n",
      "[4, 1]：-7.17e-02, -7.17e-02, 7.76e-06\n",
      "[4, 3]：3.25e-01, 3.25e-01, 2.93e-06\n",
      "[4, 4]：-5.91e-02, -5.91e-02, 1.88e-06\n",
      "[4, 5]：-7.78e-02, -7.78e-02, 5.82e-06\n",
      "[4, 6]：-9.08e-02, -9.08e-02, 8.04e-06\n",
      "[4, 7]：-1.12e-01, -1.12e-01, 2.11e-06\n",
      "[4, 8]：-1.26e-01, -1.26e-01, 5.92e-06\n",
      "[4, 9]：-1.40e-01, -1.40e-01, 4.89e-06\n",
      "[5, 0]：1.86e-01, 1.86e-01, 6.60e-06\n",
      "[5, 2]：-9.52e-02, -9.52e-02, 5.01e-06\n",
      "[5, 3]：4.97e-01, 4.97e-01, 6.92e-06\n",
      "[5, 4]：2.50e-02, 2.50e-02, 5.20e-06\n",
      "[5, 5]：-7.93e-02, -7.93e-02, 2.78e-06\n",
      "[5, 6]：-7.59e-02, -7.59e-02, 1.22e-05\n",
      "[5, 7]：-1.05e-01, -1.05e-01, 9.53e-06\n",
      "[5, 8]：-1.25e-01, -1.25e-01, 9.07e-06\n",
      "[5, 9]：-6.74e-02, -6.74e-02, 3.08e-06\n",
      "[6, 0]：1.49e-01, 1.49e-01, 5.17e-06\n",
      "[6, 1]：-6.67e-02, -6.66e-02, 4.40e-06\n",
      "[6, 2]：-1.38e-01, -1.38e-01, 1.23e-05\n",
      "[6, 3]：3.33e-01, 3.33e-01, 1.02e-05\n",
      "[6, 4]：1.72e-01, 1.72e-01, 7.79e-06\n",
      "[6, 6]：-6.73e-02, -6.73e-02, 7.13e-06\n",
      "[6, 7]：-6.08e-02, -6.08e-02, 8.58e-06\n",
      "[6, 8]：-9.98e-02, -9.98e-02, 6.14e-06\n",
      "[6, 9]：-1.31e-01, -1.31e-01, 8.89e-06\n",
      "[7, 0]：2.39e-01, 2.39e-01, 2.53e-06\n",
      "[7, 2]：-1.61e-01, -1.61e-01, 1.78e-06\n",
      "[7, 3]：1.02e-01, 1.02e-01, 4.73e-06\n",
      "[7, 4]：3.63e-01, 3.63e-01, 7.51e-06\n",
      "[7, 5]：-6.02e-02, -6.02e-02, 4.63e-06\n",
      "[7, 6]：-1.05e-01, -1.05e-01, 3.29e-06\n",
      "[7, 8]：-8.83e-02, -8.83e-02, 4.16e-06\n",
      "[7, 9]：-6.05e-02, -6.05e-02, 2.89e-06\n",
      "[8, 0]：2.92e-01, 2.92e-01, 1.38e-06\n",
      "[8, 1]：-9.70e-02, -9.70e-02, 2.27e-06\n",
      "[8, 2]：-9.87e-02, -9.87e-02, 4.36e-06\n",
      "[8, 4]：5.07e-01, 5.07e-01, 2.30e-06\n",
      "[8, 5]：-1.28e-01, -1.28e-01, 7.41e-06\n",
      "[8, 6]：-1.09e-01, -1.09e-01, 6.32e-06\n",
      "[8, 7]：-1.20e-01, -1.20e-01, 6.01e-06\n",
      "[8, 8]：-7.00e-02, -7.00e-02, 1.44e-06\n",
      "[8, 9]：-1.55e-01, -1.55e-01, 2.47e-06\n",
      "[9, 0]：4.90e-01, 4.90e-01, 2.04e-06\n",
      "[9, 1]：-7.69e-02, -7.69e-02, 1.54e-06\n",
      "[9, 2]：-9.72e-02, -9.72e-02, 9.59e-06\n",
      "[9, 4]：2.61e-01, 2.61e-01, 2.16e-06\n",
      "[9, 5]：-9.27e-02, -9.27e-02, 5.07e-06\n",
      "[9, 6]：-7.70e-02, -7.70e-02, 1.03e-06\n",
      "[9, 7]：-8.30e-02, -8.30e-02, 5.42e-06\n",
      "[9, 8]：-1.17e-01, -1.17e-01, 4.26e-06\n",
      "[9, 9]：-1.39e-01, -1.39e-01, 3.53e-06\n"
     ]
    }
   ],
   "source": [
    "def check_grad_logits(x, labels, w=-1, v=-1, toleration=1e-3):\n",
    "    grad_1 = gradient_logits(softmax(x), labels)[w, v]\n",
    "    \n",
    "    delta = 1e-10\n",
    "    original = x[w, v]\n",
    "    \n",
    "    x[w, v] = original + delta\n",
    "    y = softmax(x)\n",
    "    alpha = forward(y, labels)\n",
    "    log_p1 = np.log(alpha[-1, -1] + alpha[-1, -2])\n",
    "    \n",
    "    x[w, v] = original - delta\n",
    "    y = softmax(x)\n",
    "    alpha = forward(y, labels)\n",
    "    log_p2 = np.log(alpha[-1, -1] + alpha[-1, -2])\n",
    "    \n",
    "    x[w, v] = original\n",
    "    \n",
    "    grad_2 = (log_p1 - log_p2) / (2 * delta)\n",
    "    if np.abs(grad_1 - grad_2) > toleration:\n",
    "        print('[%d, %d]：%.2e, %.2e, %.2e' % (w, v, grad_1, grad_2, np.abs(grad_1 - grad_2)))\n",
    "\n",
    "np.random.seed(1111)\n",
    "x = np.random.random([10, 10])\n",
    "for toleration in [1e-5, 1e-6]:\n",
    "    print('%.e' % toleration)\n",
    "    for w in range(x.shape[0]):\n",
    "        for v in range(x.shape[1]):\n",
    "            check_grad_logits(x, labels, w, v, toleration)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. 数值稳定性\n",
    "\n",
    "CTC 的训练过程面临数值下溢的风险，特别是序列较大的情况下。下面介绍两种数值上稳定的工程优化方法：1）log 域（许多 CRF 实现的常用方法）；2）scale 技巧（原始论文 [1] 使用的方法）。\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1 log 域计算\n",
    "\n",
    "log 计算涉及 logsumexp 操作。\n",
    "[经验表明](https://github.com/baidu-research/warp-ctc)，在 log 域计算，即使使用单精度，也表现出良好的数值稳定性，可以有效避免下溢的风险。稳定性的代价是增加了运算的复杂性——原始实现只涉及乘加运算，log 域实现则需要对数和指数运算。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ninf = -np.float('inf')\n",
    "\n",
    "def _logsumexp(a, b):\n",
    "    '''\n",
    "    np.log(np.exp(a) + np.exp(b))\n",
    "\n",
    "    '''\n",
    "    \n",
    "    if a < b:\n",
    "        a, b = b, a\n",
    "        \n",
    "    if b == ninf:\n",
    "        return a\n",
    "    else:\n",
    "        return a + np.log(1 + np.exp(b - a)) \n",
    "    \n",
    "def logsumexp(*args):\n",
    "    '''\n",
    "    from scipy.special import logsumexp\n",
    "    logsumexp(args)\n",
    "    '''\n",
    "    res = args[0]\n",
    "    for e in args[1:]:\n",
    "        res = _logsumexp(res, e)\n",
    "    return res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1.1 log 域前向算法\n",
    "基于 log 的前向算法实现如下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8.60881935942e-17\n"
     ]
    }
   ],
   "source": [
    "def forward_log(log_y, labels):\n",
    "    T, V = log_y.shape\n",
    "    L = len(labels)\n",
    "    log_alpha = np.ones([T, L]) * ninf\n",
    "\n",
    "    # init\n",
    "    log_alpha[0, 0] = log_y[0, labels[0]]\n",
    "    log_alpha[0, 1] = log_y[0, labels[1]]\n",
    "\n",
    "    for t in range(1, T):\n",
    "        for i in range(L):\n",
    "            s = labels[i]\n",
    "            \n",
    "            a = log_alpha[t - 1, i]\n",
    "            if i - 1 >= 0:\n",
    "                a = logsumexp(a, log_alpha[t - 1, i - 1])\n",
    "            if i - 2 >= 0 and s != 0 and s != labels[i - 2]:\n",
    "                a = logsumexp(a, log_alpha[t - 1, i - 2])\n",
    "                \n",
    "            log_alpha[t, i] = a + log_y[t, s]\n",
    "            \n",
    "    return log_alpha\n",
    "\n",
    "log_alpha = forward_log(np.log(y), labels)\n",
    "alpha = forward(y, labels)\n",
    "print(np.sum(np.abs(np.exp(log_alpha) - alpha)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1.2 log 域后向算法\n",
    "基于 log 的后向算法实现如下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.10399945005e-16\n"
     ]
    }
   ],
   "source": [
    "def backward_log(log_y, labels):\n",
    "    T, V = log_y.shape\n",
    "    L = len(labels)\n",
    "    log_beta = np.ones([T, L]) * ninf\n",
    "\n",
    "    # init\n",
    "    log_beta[-1, -1] = log_y[-1, labels[-1]]\n",
    "    log_beta[-1, -2] = log_y[-1, labels[-2]]\n",
    "\n",
    "    for t in range(T - 2, -1, -1):\n",
    "        for i in range(L):\n",
    "            s = labels[i]\n",
    "            \n",
    "            a = log_beta[t + 1, i] \n",
    "            if i + 1 < L:\n",
    "                a = logsumexp(a, log_beta[t + 1, i + 1])\n",
    "            if i + 2 < L and s != 0 and s != labels[i + 2]:\n",
    "                a = logsumexp(a, log_beta[t + 1, i + 2])\n",
    "                \n",
    "            log_beta[t, i] = a + log_y[t, s]\n",
    "            \n",
    "    return log_beta\n",
    "\n",
    "log_beta = backward_log(np.log(y), labels)\n",
    "beta = backward(y, labels)\n",
    "print(np.sum(np.abs(np.exp(log_beta) - beta)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1.3 log 域梯度计算\n",
    "在前向、后向基础上，也可以在 log 域上计算梯度。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.97588081849e-14\n"
     ]
    }
   ],
   "source": [
    "def gradient_log(log_y, labels):\n",
    "    T, V = log_y.shape\n",
    "    L = len(labels)\n",
    "    \n",
    "    log_alpha = forward_log(log_y, labels)\n",
    "    log_beta = backward_log(log_y, labels)\n",
    "    log_p = logsumexp(log_alpha[-1, -1], log_alpha[-1, -2])\n",
    "    \n",
    "    log_grad = np.ones([T, V]) * ninf\n",
    "    for t in range(T):\n",
    "        for s in range(V):\n",
    "            lab = [i for i, c in enumerate(labels) if c == s]\n",
    "            for i in lab:\n",
    "                log_grad[t, s] = logsumexp(log_grad[t, s], log_alpha[t, i] + log_beta[t, i]) \n",
    "            log_grad[t, s] -= 2 * log_y[t, s]\n",
    "                \n",
    "    log_grad -= log_p\n",
    "    return log_grad\n",
    "    \n",
    "log_grad = gradient_log(np.log(y), labels)\n",
    "grad = gradient(y, labels)\n",
    "#print(log_grad)\n",
    "#print(grad)\n",
    "print(np.sum(np.abs(np.exp(log_grad) - grad)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2 scale\n",
    "\n",
    "### 2.2.1 前向算法\n",
    "\n",
    "为了下溢，在前向算法的每个时刻，都对计算出的 $\\alpha$ 的范围进行缩放：\n",
    "$$\n",
    "C_t \\stackrel{def}{=} \\sum_s\\alpha_t(s)\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\hat{\\alpha}_t = \\frac{\\alpha_t(s)}{C_t}\n",
    "$$\n",
    "\n",
    "缩放后的 $\\alpha$，不会随着时刻的积累变得太小。$\\hat{\\alpha}$ 替代 $\\alpha$，进行下一时刻的迭代。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def forward_scale(y, labels):\n",
    "    T, V = y.shape\n",
    "    L = len(labels)\n",
    "    alpha_scale = np.zeros([T, L])\n",
    "\n",
    "    # init\n",
    "    alpha_scale[0, 0] = y[0, labels[0]]\n",
    "    alpha_scale[0, 1] = y[0, labels[1]]\n",
    "    Cs = []\n",
    "    \n",
    "    C = np.sum(alpha_scale[0])\n",
    "    alpha_scale[0] /= C\n",
    "    Cs.append(C)\n",
    "\n",
    "    for t in range(1, T):\n",
    "        for i in range(L):\n",
    "            s = labels[i]\n",
    "            \n",
    "            a = alpha_scale[t - 1, i] \n",
    "            if i - 1 >= 0:\n",
    "                a += alpha_scale[t - 1, i - 1]\n",
    "            if i - 2 >= 0 and s != 0 and s != labels[i - 2]:\n",
    "                a += alpha_scale[t - 1, i - 2]\n",
    "                \n",
    "            alpha_scale[t, i] = a * y[t, s]\n",
    "            \n",
    "        C = np.sum(alpha_scale[t])\n",
    "        alpha_scale[t] /= C\n",
    "        Cs.append(C)\n",
    "            \n",
    "    return alpha_scale, Cs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由于进行了缩放，最后计算概率时要时行补偿：\n",
    "$$\n",
    "p(l|x) = \\alpha_T(|l^\\prime|) + \\alpha_T(|l^\\prime|-1) = (\\hat\\alpha_T(|l^\\prime|) + \\hat\\alpha_T(|l^\\prime|-1) * \\prod_{t=1}^T C_t\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\ln(p(l|x)) = \\sum_t^T\\ln(C_t) + \\ln(\\hat\\alpha_T(|l^\\prime|) + \\hat\\alpha_T(|l^\\prime|-1))\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(-13.202925982240107, -13.202925982240107, 0.0)\n"
     ]
    }
   ],
   "source": [
    "labels = [0, 1, 2, 0]  # 0 for blank\n",
    "\n",
    "alpha_scale, Cs = forward_scale(y, labels)\n",
    "log_p = np.sum(np.log(Cs)) + np.log(alpha_scale[-1][labels[-1]] + alpha_scale[-1][labels[-2]])\n",
    "\n",
    "alpha = forward(y, labels)\n",
    "p = alpha[-1, labels[-1]] + alpha[-1, labels[-2]]\n",
    "\n",
    "print(np.log(p), log_p, np.log(p) - log_p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2.2 后向算法\n",
    "后向算法缩放类似于前向算法，公式如下：\n",
    "\n",
    "$$\n",
    "D_t \\stackrel{def}{=} \\sum_s\\beta_t(s)\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\hat{\\beta}_t = \\frac{\\beta_t(s)}{D_t}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.71362347  0.18910147  0.07964328  0.01763178]\n",
      " [ 0.70165268  0.15859852  0.11849423  0.02125457]\n",
      " [ 0.67689676  0.165374    0.13221504  0.02551419]\n",
      " [ 0.71398181  0.11936432  0.13524265  0.03141122]\n",
      " [ 0.70769657  0.13093688  0.12447135  0.0368952 ]\n",
      " [ 0.63594568  0.1790638   0.14250065  0.04248987]\n",
      " [ 0.63144322  0.1806382   0.13366043  0.05425815]\n",
      " [ 0.33926289  0.35149591  0.24988622  0.05935497]\n",
      " [ 0.30303623  0.26644554  0.33088584  0.0996324 ]\n",
      " [ 0.12510056  0.3297143   0.3956509   0.14953425]\n",
      " [ 0.          0.22078343  0.5153114   0.26390517]\n",
      " [ 0.          0.          0.550151    0.449849  ]]\n"
     ]
    }
   ],
   "source": [
    "def backward_scale(y, labels):\n",
    "    T, V = y.shape\n",
    "    L = len(labels)\n",
    "    beta_scale = np.zeros([T, L])\n",
    "\n",
    "    # init\n",
    "    beta_scale[-1, -1] = y[-1, labels[-1]]\n",
    "    beta_scale[-1, -2] = y[-1, labels[-2]]\n",
    "    \n",
    "    Ds = []\n",
    "    \n",
    "    D = np.sum(beta_scale[-1,:])\n",
    "    beta_scale[-1] /= D\n",
    "    Ds.append(D)\n",
    "\n",
    "    for t in range(T - 2, -1, -1):\n",
    "        for i in range(L):\n",
    "            s = labels[i]\n",
    "            \n",
    "            a = beta_scale[t + 1, i] \n",
    "            if i + 1 < L:\n",
    "                a += beta_scale[t + 1, i + 1]\n",
    "            if i + 2 < L and s != 0 and s != labels[i + 2]:\n",
    "                a += beta_scale[t + 1, i + 2]\n",
    "                \n",
    "            beta_scale[t, i] = a * y[t, s]\n",
    "            \n",
    "        D = np.sum(beta_scale[t])\n",
    "        beta_scale[t] /= D\n",
    "        Ds.append(D)\n",
    "            \n",
    "    return beta_scale, Ds[::-1]\n",
    "\n",
    "beta_scale, Ds = backward_scale(y, labels)\n",
    "print(beta_scale)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2.3 梯度计算\n",
    "\n",
    "$$\n",
    "\\frac{\\partial \\ln(p(l|x))}{\\partial y^t_k} = \\frac{1}{p(l|x)} \\frac{\\partial p(l|x)}{\\partial y^t_k} = \\frac{1}{p(l|x)} \\frac{1}{{y^t_k}^2} \\sum_{s \\in lab(l, k)} \\alpha_t(s)\\beta_t(s) \n",
    "$$\n",
    "\n",
    "考虑到\n",
    "$$\n",
    "p(l|x) =  \\sum_{s=1}^{|l^\\prime|}  \\frac{\\alpha_t(s)\\beta_t(s)}{ y^t_{l_s^\\prime}} \n",
    "$$\n",
    "以及\n",
    "$$\n",
    "\\alpha_t(s) = \\hat\\alpha_t(s) \\cdot \\prod_{k=1}^t C_k\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\beta_t(s) = \\hat\\beta_t(s) \\cdot \\prod_{k=t}^T D_k\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\frac{\\partial \\ln(p(l|x))}{\\partial y^t_k} = \\frac{1}{\\sum_{s=1}^{|l^\\prime|} \\frac{\\hat\\alpha_t(s)  \\hat\\beta_t(s)}{y^t_{l^\\prime_s}}} \\frac{1}{{y^t_k}^2} \\sum_{s \\in lab(l, k)} \\hat\\alpha_t(s)  \\hat\\beta_t(s)\n",
    "$$\n",
    "\n",
    "式中最右项中的各个部分我们都已经求得。梯度计算实现如下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.86256607096e-15\n"
     ]
    }
   ],
   "source": [
    "def gradient_scale(y, labels):\n",
    "    T, V = y.shape\n",
    "    L = len(labels)\n",
    "    \n",
    "    alpha_scale, _ = forward_scale(y, labels)\n",
    "    beta_scale, _ = backward_scale(y, labels)\n",
    "    \n",
    "    grad = np.zeros([T, V])\n",
    "    for t in range(T):\n",
    "        for s in range(V):\n",
    "            lab = [i for i, c in enumerate(labels) if c == s]\n",
    "            for i in lab:\n",
    "                grad[t, s] += alpha_scale[t, i] * beta_scale[t, i]\n",
    "            grad[t, s] /= y[t, s] ** 2\n",
    "         \n",
    "        # normalize factor\n",
    "        z = 0\n",
    "        for i in range(L):\n",
    "            z += alpha_scale[t, i] * beta_scale[t, i] / y[t, labels[i]]\n",
    "        grad[t] /= z\n",
    "        \n",
    "    return grad\n",
    "    \n",
    "labels = [0, 3, 0, 3, 0, 4, 0]  # 0 for blank\n",
    "grad_1 = gradient_scale(y, labels)\n",
    "grad_2 = gradient(y, labels)\n",
    "print(np.sum(np.abs(grad_1 - grad_2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2.4 logits 梯度\n",
    "类似于 y 梯度的推导，logits 梯度计算公式如下：\n",
    "$$\n",
    "\\frac{\\partial \\ln(p(l|x))}{\\partial u^t_k} = \\frac{1}{y^t_k Z_t} \\sum_{s \\in lab(l, k)} \\hat\\alpha_t(s)\\hat\\beta_t(s) - y^t_k\n",
    "$$\n",
    "其中，\n",
    "$$\n",
    "Z_t \\stackrel{def}{=} \\sum_{s=1}^{|l^\\prime|} \\frac{\\hat\\alpha_t(s)\\hat\\beta_t(s)}{y^t_{l^\\prime_s}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. 解码\n",
    "训练和的 $N_w$ 可以用来预测新的样本输入对应的输出字符串，这涉及到解码。\n",
    "按照最大似然准则，最优的解码结果为：\n",
    "$$\n",
    "h(x) = \\underset{l \\in L^{\\le T}}{\\mathrm{argmax}}\\ p(l|x)\n",
    "$$\n",
    "\n",
    "然而，上式不存在已知的高效解法。下面介绍几种实用的近似破解码方法。\n",
    "\n",
    "## 3.1 贪心搜索 （greedy search）\n",
    "虽然 $p(l|x)$ 难以有效的计算，但是由于 CTC 的独立性假设，对于某个具体的字符串 $\\pi$（去 blank 前），确容易计算：\n",
    "$$\n",
    "p(\\pi|x) = \\prod_{k=1}^T p(\\pi_k|x)\n",
    "$$\n",
    "\n",
    "因此，我们放弃寻找使 $p(l|x)$ 最大的字符串，退而寻找一个使 $p(\\pi|x)$ 最大的字符串，即：\n",
    "\n",
    "$$\n",
    "h(x) \\approx B(\\pi^\\star)\n",
    "$$\n",
    "其中，\n",
    "$$\n",
    "\\pi^\\star = \\underset{\\pi \\in N^T}{\\mathrm{argmax}}\\ p(\\pi|x)\n",
    "$$\n",
    "\n",
    "简化后，解码过程（构造 $\\pi^\\star$）变得非常简单（基于独立性假设）： 在每个时刻输出概率最大的字符:\n",
    "$$\n",
    "\\pi^\\star = cat_{t=1}^T(\\underset{s \\in L^\\prime}{\\mathrm{argmax}}\\ y^t_s)\n",
    "$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 3 5 5 5 5 1 5 3 4 4 3 0 4 5 0 3 1 3 3]\n",
      "[1, 3, 5, 1, 5, 3, 4, 3, 4, 5, 3, 1, 3]\n"
     ]
    }
   ],
   "source": [
    "def remove_blank(labels, blank=0):\n",
    "    new_labels = []\n",
    "    \n",
    "    # combine duplicate\n",
    "    previous = None\n",
    "    for l in labels:\n",
    "        if l != previous:\n",
    "            new_labels.append(l)\n",
    "            previous = l\n",
    "            \n",
    "    # remove blank     \n",
    "    new_labels = [l for l in new_labels if l != blank]\n",
    "    \n",
    "    return new_labels\n",
    "\n",
    "def insert_blank(labels, blank=0):\n",
    "    new_labels = [blank]\n",
    "    for l in labels:\n",
    "        new_labels += [l, blank]\n",
    "    return new_labels\n",
    "\n",
    "def greedy_decode(y, blank=0):\n",
    "    raw_rs = np.argmax(y, axis=1)\n",
    "    rs = remove_blank(raw_rs, blank)\n",
    "    return raw_rs, rs\n",
    "\n",
    "np.random.seed(1111)\n",
    "y = softmax(np.random.random([20, 6]))\n",
    "rr, rs = greedy_decode(y)\n",
    "print(rr)\n",
    "print(rs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.2 束搜索（Beam Search）\n",
    "显然，贪心搜索的性能非常受限。例如，它不能给出除最优路径之外的其他其优路径。很多时候，如果我们能拿到 nbest 的路径，后续可以利用其他信息来进一步优化搜索的结果。束搜索能近似找出 top 最优的若干条路径。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "([1, 3, 5, 1, 5, 3, 4, 3, 4, 5, 3, 1, 3], -29.261797539205567)\n",
      "([1, 3, 5, 1, 5, 3, 4, 3, 3, 5, 3, 1, 3], -29.279020152518033)\n",
      "([1, 3, 5, 1, 5, 3, 4, 2, 3, 4, 5, 3, 1, 3], -29.300726142201842)\n",
      "([1, 5, 1, 5, 3, 4, 3, 4, 5, 3, 1, 3], -29.310307014773972)\n",
      "([1, 3, 5, 1, 5, 3, 4, 2, 3, 3, 5, 3, 1, 3], -29.317948755514308)\n",
      "([1, 5, 1, 5, 3, 4, 3, 3, 5, 3, 1, 3], -29.327529628086438)\n",
      "([1, 3, 5, 1, 5, 4, 3, 4, 5, 3, 1, 3], -29.331572723457334)\n",
      "([1, 3, 5, 5, 1, 5, 3, 4, 3, 4, 5, 3, 1, 3], -29.332631809924511)\n",
      "([1, 3, 5, 4, 1, 5, 3, 4, 3, 4, 5, 3, 1, 3], -29.334649090836038)\n",
      "([1, 3, 5, 1, 5, 3, 4, 3, 4, 5, 3, 1, 3], -29.33969505198154)\n",
      "([1, 3, 5, 2, 1, 5, 3, 4, 3, 4, 5, 3, 1, 3], -29.339823066915415)\n",
      "([1, 3, 5, 1, 5, 4, 3, 3, 5, 3, 1, 3], -29.3487953367698)\n",
      "([1, 5, 1, 5, 3, 4, 2, 3, 4, 5, 3, 1, 3], -29.349235617770248)\n",
      "([1, 3, 5, 5, 1, 5, 3, 4, 3, 3, 5, 3, 1, 3], -29.349854423236977)\n",
      "([1, 3, 5, 1, 5, 3, 4, 3, 4, 5, 3, 3], -29.350803198551016)\n",
      "([1, 3, 5, 4, 1, 5, 3, 4, 3, 3, 5, 3, 1, 3], -29.351871704148504)\n",
      "([1, 3, 5, 1, 5, 3, 4, 3, 3, 5, 3, 1, 3], -29.356917665294006)\n",
      "([1, 3, 5, 2, 1, 5, 3, 4, 3, 3, 5, 3, 1, 3], -29.357045680227881)\n",
      "([1, 3, 5, 1, 5, 3, 4, 5, 4, 5, 3, 1, 3], -29.363802591012263)\n",
      "([1, 5, 1, 5, 3, 4, 2, 3, 3, 5, 3, 1, 3], -29.366458231082714)\n"
     ]
    }
   ],
   "source": [
    "def beam_decode(y, beam_size=10):\n",
    "    T, V = y.shape\n",
    "    log_y = np.log(y)\n",
    "    \n",
    "    beam = [([], 0)]\n",
    "    for t in range(T):  # for every timestep\n",
    "        new_beam = []\n",
    "        for prefix, score in beam:\n",
    "            for i in range(V):  # for every state\n",
    "                new_prefix = prefix + [i]\n",
    "                new_score = score + log_y[t, i]\n",
    "                \n",
    "                new_beam.append((new_prefix, new_score))\n",
    "                \n",
    "        # top beam_size\n",
    "        new_beam.sort(key=lambda x: x[1], reverse=True)\n",
    "        beam = new_beam[:beam_size]\n",
    "        \n",
    "    return beam\n",
    "    \n",
    "np.random.seed(1111)\n",
    "y = softmax(np.random.random([20, 6]))\n",
    "beam = beam_decode(y, beam_size=100)\n",
    "for string, score in beam[:20]:\n",
    "    print(remove_blank(string), score)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.3 前缀束搜索（Prefix Beam Search）\n",
    "直接的束搜索的一个问题是，在保存的 top N 条路径中，可能存在多条实际上是同一结果（经过去重复、去 blank 操作）。这减少了搜索结果的多样性。下面介绍的前缀搜索方法，在搜索过程中不断的合并相同的前缀[2]。参考 [gist](https://gist.github.com/awni/56369a90d03953e370f3964c826ed4b0)，前缀束搜索实现如下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "([1, 5, 4, 1, 3, 4, 5, 2, 3], (-18.189863809114193, -17.613677981426175))\n",
      "([1, 5, 4, 5, 3, 4, 5, 2, 3], (-18.19636512622969, -17.621013424585406))\n",
      "([1, 5, 4, 1, 3, 4, 5, 1, 3], (-18.317018960331531, -17.666629973270073))\n",
      "([1, 5, 4, 5, 3, 4, 5, 1, 3], (-18.323388267369936, -17.674125139073176))\n",
      "([1, 5, 4, 1, 3, 4, 3, 2, 3], (-18.415808498759556, -17.862744326248826))\n",
      "([1, 5, 4, 1, 3, 4, 3, 5, 3], (-18.366422766638632, -17.898463479112884))\n",
      "([1, 5, 4, 5, 3, 4, 3, 2, 3], (-18.42224294936932, -17.870025672291458))\n",
      "([1, 5, 4, 5, 3, 4, 3, 5, 3], (-18.372199113900191, -17.905130493229173))\n",
      "([1, 5, 4, 1, 3, 4, 5, 4, 3], (-18.457066311773847, -17.880630315602037))\n",
      "([1, 5, 4, 5, 3, 4, 5, 4, 3], (-18.462614293487096, -17.88759583852546))\n",
      "([1, 5, 4, 1, 3, 4, 5, 3, 2], (-18.458941701567706, -17.951422824358747))\n",
      "([1, 5, 4, 5, 3, 4, 5, 3, 2], (-18.464527031120184, -17.958629487208658))\n",
      "([1, 5, 4, 1, 3, 4, 3, 1, 3], (-18.540857550725587, -17.920589910093689))\n",
      "([1, 5, 4, 5, 3, 4, 3, 1, 3], (-18.547146092248852, -17.928030266681613))\n",
      "([1, 5, 4, 1, 3, 4, 5, 3, 2, 3], (-19.325467801462263, -17.689203224408899))\n",
      "([1, 5, 4, 5, 3, 4, 5, 3, 2, 3], (-19.328748799764973, -17.694105969982637))\n",
      "([1, 5, 4, 1, 3, 4, 5, 3, 4], (-18.79699026165903, -17.945090229238392))\n",
      "([1, 5, 4, 5, 3, 4, 5, 3, 4], (-18.803585534273239, -17.95258394264377))\n",
      "([1, 5, 4, 3, 4, 3, 5, 2, 3], (-19.181531846082809, -17.859420073785095))\n",
      "([1, 5, 4, 1, 3, 4, 5, 2, 3, 2], (-19.439349296385199, -17.884502168470895))\n"
     ]
    }
   ],
   "source": [
    "from collections import defaultdict\n",
    "\n",
    "def prefix_beam_decode(y, beam_size=10, blank=0):\n",
    "    T, V = y.shape\n",
    "    log_y = np.log(y)\n",
    "    \n",
    "    beam = [(tuple(), (0, ninf))]  # blank, non-blank\n",
    "    for t in range(T):  # for every timestep\n",
    "        new_beam = defaultdict(lambda : (ninf, ninf))\n",
    "             \n",
    "        for prefix, (p_b, p_nb) in beam:\n",
    "            for i in range(V):  # for every state\n",
    "                p = log_y[t, i]\n",
    "                \n",
    "                if i == blank:  # propose a blank\n",
    "                    new_p_b, new_p_nb = new_beam[prefix]\n",
    "                    new_p_b = logsumexp(new_p_b, p_b + p, p_nb + p)\n",
    "                    new_beam[prefix] = (new_p_b, new_p_nb)\n",
    "                    continue\n",
    "                else:  # extend with non-blank\n",
    "                    end_t = prefix[-1] if prefix else None\n",
    "                    \n",
    "                    # exntend current prefix\n",
    "                    new_prefix = prefix + (i,)\n",
    "                    new_p_b, new_p_nb = new_beam[new_prefix]\n",
    "                    if i != end_t:\n",
    "                        new_p_nb = logsumexp(new_p_nb, p_b + p, p_nb + p)\n",
    "                    else:\n",
    "                        new_p_nb = logsumexp(new_p_nb, p_b + p)\n",
    "                    new_beam[new_prefix] = (new_p_b, new_p_nb)\n",
    "                    \n",
    "                    # keep current prefix\n",
    "                    if i == end_t:\n",
    "                        new_p_b, new_p_nb = new_beam[prefix]\n",
    "                        new_p_nb = logsumexp(new_p_nb, p_nb + p)\n",
    "                        new_beam[prefix] = (new_p_b, new_p_nb)\n",
    "                \n",
    "        # top beam_size\n",
    "        beam = sorted(new_beam.items(), key=lambda x : logsumexp(*x[1]), reverse=True)\n",
    "        beam = beam[:beam_size]\n",
    "        \n",
    "    return beam\n",
    "\n",
    "np.random.seed(1111)\n",
    "y = softmax(np.random.random([20, 6]))\n",
    "beam = prefix_beam_decode(y, beam_size=100)\n",
    "for string, score in beam[:20]:\n",
    "    print(remove_blank(string), score)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.4 其他解码方法\n",
    "上述介绍了基本解码方法。实际中，搜索过程可以接合额外的信息，提高搜索的准确度（例如在语音识别任务中，加入语言模型得分信息, [前缀搜索+语言模型](https://github.com/PaddlePaddle/DeepSpeech/blob/develop/decoders/decoders_deprecated.py\n",
    ")）。\n",
    "\n",
    "本质上，CTC 只是一个训练准则。训练完成后，$N_w$ 输出一系列概率分布，这点和常规基于交叉熵准则训练的模型完全一致。因此，特定应用领域常规的解码也可以经过一定修改用来 CTC 的解码。例如在语音识别任务中，利用 CTC 训练的声学模型可以无缝融入原来的 WFST 的解码器中[5]（e.g. 参见 [EESEN](https://github.com/srvk/eesen)）。\n",
    "\n",
    "此外，[1] 给出了一种利用 CTC 顶峰特点的启发式搜索方法。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. 工具\n",
    "\n",
    "[warp-ctc](https://github.com/baidu-research/warp-ctc) 是百度开源的基于 CPU 和 GPU 的高效并行实现。warp-ctc 自身提供 C 语言接口，对于流利的机器学习工具（[torch](https://github.com/baidu-research/warp-ctc/tree/master/torch_binding)、 [pytorch](https://github.com/SeanNaren/deepspeech.pytorch) 和 [tensorflow](https://github.com/baidu-research/warp-ctc/tree/master/tensorflow_binding)、[chainer](https://github.com/jheymann85/chainer_ctc)）都有相应的接口绑定。\n",
    "\n",
    "[cudnn 7](https://developer.nvidia.com/cudnn) 以后开始提供 CTC 支持。\n",
    "\n",
    "Tensorflow 也原生支持 [CTC loss](https://www.tensorflow.org/api_docs/python/tf/nn/ctc_loss)，及 greedy 和 beam search 解码器。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 小结\n",
    "1. CTC 可以建模无对齐信息的多对多序列问题（输入长度不小于输出），如语音识别、连续字符识别 [3,4]。\n",
    "2. CTC 不需要输入与输出的对齐信息，可以实现端到端的训练。\n",
    "3. CTC 在 loss 的计算上，利用了整个 labels 序列的全局信息，某种意义上相对逐帧计算损失的方法，\"更加区分性\"。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References\n",
    "1. Graves et al. [Connectionist Temporal Classification : Labelling Unsegmented Sequence Data with Recurrent Neural Networks](ftp://ftp.idsia.ch/pub/juergen/icml2006.pdf).\n",
    "2. Hannun et al. [First-Pass Large Vocabulary Continuous Speech Recognition using Bi-Directional Recurrent DNNs](https://arxiv.org/abs/1408.2873).\n",
    "3. Graves et al. [Towards End-To-End Speech Recognition with Recurrent Neural Networks](http://jmlr.org/proceedings/papers/v32/graves14.pdf).\n",
    "4. Liwicki et al. [A novel approach to on-line handwriting recognition based on bidirectional long short-term memory networks](https://www.cs.toronto.edu/~graves/icdar_2007.pdf).\n",
    "5. Zenkel et al. [Comparison of Decoding Strategies for CTC Acoustic Models](https://arxiv.org/abs/1708.004469).\n",
    "5. Huannun. [Sequence Modeling with CTC](https://distill.pub/2017/ctc/)."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}