CTC（ Connectionist Temporal Classification，连接时序分类）是一种用于序列建模的工具，其核心是定义了特殊的目标函数/优化准则[1]。

jupyter notebook 版见 repo.

1. 算法

这里大体根据 Alex Graves 的开山之作[1]，讨论 CTC 的算法原理，并基于 numpy 从零实现 CTC 的推理及训练算法。

1.1 序列问题形式化。

序列问题可以形式化为如下函数：

$$mathcal{N}_w: (mathcal{R}^m)^T rightarrow (mathcal{R}^n)^T$$
其中，序列目标为字符串（词表大小为 $n$），即 $mathcal{N}_w$ 输出为 $n$ 维多项概率分布（e.g. 经过 softmax 处理）。

网络输出为：$y = mathcal{N}_w$，其中， $y_k^t$ $t$ 表示时刻第 $k$ 项的概率。

图1. 序列建模【src】

虽然并没为限定 $mathcal{N}_w$ 具体形式，下面为假设其了某种神经网络（e.g. RNN）。
下面代码示例 toy $mathcal{N}_w$：

import numpy as np

np.random.seed(1111)

T, V = 12, 5
m, n = 6, V

x = np.random.random([T, m])  # T x m
w = np.random.random([m, n])  # weights, m x n

def softmax(logits):
    max_value = np.max(logits, axis=1, keepdims=True)
    exp = np.exp(logits - max_value)
    exp_sum = np.sum(exp, axis=1, keepdims=True)
    dist = exp / exp_sum
    return dist

def toy_nw(x):
    y = np.matmul(x, w)  # T x n 
    y = softmax(y)
    return y

y = toy_nw(x)
print(y)
print(y.sum(1, keepdims=True))

[[ 0.24654511  0.18837589  0.16937668  0.16757465  0.22812766]
 [ 0.25443629  0.14992236  0.22945293  0.17240658  0.19378184]
 [ 0.24134404  0.17179604  0.23572466  0.12994237  0.22119288]
 [ 0.27216255  0.13054313  0.2679252   0.14184499  0.18752413]
 [ 0.32558002  0.13485564  0.25228604  0.09743785  0.18984045]
 [ 0.23855586  0.14800386  0.23100255  0.17158135  0.21085638]
 [ 0.38534786  0.11524603  0.18220093  0.14617864  0.17102655]
 [ 0.21867406  0.18511892  0.21305488  0.16472572  0.21842642]
 [ 0.29856607  0.13646801  0.27196606  0.11562552  0.17737434]
 [ 0.242347    0.14102063  0.21716951  0.2355229   0.16393996]
 [ 0.26597326  0.10009752  0.23362892  0.24560198  0.15469832]
 [ 0.23337289  0.11918746  0.28540761  0.20197928  0.16005275]]
[[ 1.]
 [ 1.]
 [ 1.]
 [ 1.]
 [ 1.]
 [ 1.]
 [ 1.]
 [ 1.]
 [ 1.]
 [ 1.]
 [ 1.]
 [ 1.]]

1.2 align-free 变长映射

上面的形式是输入和输出的一对一的映射。序列学习任务一般而言是多对多的映射关系（如语音识别中，上百帧输出可能仅对应若干音节或字符，并且每个输入和输出之间，也没有清楚的对应关系）。CTC 通过引入一个特殊的 blank 字符（用 % 表示），解决多对一映射问题。

扩展原始词表 $L$ 为 $L^prime = L cup {text{blank}}$。对输出字符串，定义操作 $mathcal{B}$：1）合并连续的相同符号；2）去掉 blank 字符。

例如，对于 “aa%bb%%cc”，应用 $mathcal{B}$，则实际上代表的是字符串 “abc”。同理“%a%b%cc%” 也同样代表 “abc”。

$$mathcal{B}(aa%bb%%cc) = mathcal{B}(%a%b%cc%) = abc$$

通过引入blank 及 $mathcal{B}$，可以实现了变长的映射。

$$L^{prime T} rightarrow L^{le T}$$

因为这个原因，CTC 只能建模输出长度小于输入长度的序列问题。

1.3 似然计算

和大多数有监督学习一样，CTC 使用最大似然标准进行训练。

给定输入 $x$，输出 $l$ 的条件概率为：

$$p(l|x) = sum_{pi in mathcal{B}^{-1}(l)} p(pi|x)$$

其中，$mathcal{B}^{-1}(l)$ 表示了长度为 $T$ 且示经过 $mathcal{B}$ 结果为 $l$ 字符串的集合。

CTC 假设输出的概率是（相对于输入）条件独立的，因此有：

$$p(pi|x) = prod y^t_{pi_t}, forall pi in L^{prime T}$$

然而，直接按上式我们没有办理有效的计算似然值。下面用动态规划解决似然的计算及梯度计算, 涉及前向算法和后向算法。

1.4 前向算法

在前向及后向计算中，CTC 需要将输出字符串进行扩展。具体的，$(a_1,cdots,a_m)$ 每个字符之间及首尾分别插入 blank，即扩展为 $(%, a_1,%,a_2, %,cdots,%, a_m,%)$。下面的 $l$ 为原始字符串，$l^prime$ 指为扩展后的字符串。

定义

$$alpha_t(s) stackrel{def}{=} sum_{pi in N^T: mathcal{B}(pi_{1:t}) = l_{1:s}} prod_{t^prime=1}^t y^t_{pi^prime}$$

显然有，

$$begin{align} alpha_1(1) = y_b^1,\ alpha_1(2) = y_{l_1}^1,\ alpha_1(s) = 0, forall s > 2 end{align}$$
根据 $alpha$ 的定义，有如下递归关系：
$$alpha_t(s) = { begin{array}{l} (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} \ (alpha_{t-1}(s)+alpha_{t-1}(s-1) + alpha_{t-1}(s-2)) y^t_{l^prime_s} otherwise end{array}$$

1.4.1 Case 2

递归公式中 case 2 是一般的情形。如图所示，$t$ 时刻字符为 $s$ 为 blank 时，它可能由于两种情况扩展而来：1）重复上一字符，即上个字符也是 a，2）字符发生转换，即上个字符是非 a 的字符。第二种情况又分为两种情形，2.1）上一字符是 blank；2.2）a 由非 blank 字符直接跳转而来（$mathcal{B}$） 操作中， blank 最终会被去掉，因此 blank 并不是必须的）。

图2. 前向算法 Case 2 示例【src】

1.4.2 Case 1

递归公式 case 1 是特殊的情形。
如图所示，$t$ 时刻字符为 $s$ 为 blank 时，它只能由于两种情况扩展而来：1）重复上一字符，即上个字符也是 blank，2）字符发生转换，即上个字符是非 blank 字符。$t$ 时刻字符为 $s$ 为非 blank 时，类似于 case 2，但是这时两个相同字符之间的 blank 不能省略（否则无法区分”aa”和”a”），因此，也只有两种跳转情况。

图3. 前向算法 Case 1 示例【src】

我们可以利用动态规划计算所有 $alpha$ 的值，算法时间和空间复杂度为 $O(T * L)$。

似然的计算只涉及乘加运算，因此，CTC 的似然是可导的，可以尝试 tensorflow 或 pytorch 等具有自动求导功能的工具自动进行梯度计算。下面介绍如何手动高效的计算梯度。

def forward(y, labels):
    T, V = y.shape
    L = len(labels)
    alpha = np.zeros([T, L])

    # init
    alpha[0, 0] = y[0, labels[0]]
    alpha[0, 1] = y[0, labels[1]]

    for t in range(1, T):
        for i in range(L):
            s = labels[i]

            a = alpha[t - 1, i] 
            if i - 1 >= 0:
                a += alpha[t - 1, i - 1]
            if i - 2 >= 0 and s != 0 and s != labels[i - 2]:
                a += alpha[t - 1, i - 2]

            alpha[t, i] = a * y[t, s]

    return alpha

labels = [0, 3, 0, 3, 0, 4, 0]  # 0 for blank
alpha = forward(y, labels)
print(alpha)

最后可以得到似然 $p(l|x) = alpha_T(|l^prime|) + alpha_T(|l^prime|-1)$。

p = alpha[-1, labels[-1]] + alpha[-1, labels[-2]]
print(p)

6.81811271177e-06

1.5 后向计算

类似于前向计算，我们定义后向计算。
首先定义

$$beta_t(s) stackrel{def}{=} sum_{pi in N^T: mathcal{B}(pi_{t:T}) = l_{s:|l|}} prod_{t^prime=t}^T y^t_{pi^prime}$$

显然，

$$begin{align} beta_T(|l^prime|) = y_b^T,\ beta_T(|l^prime|-1) = y_{l_{|l|}}^T,\ beta_T(s) = 0, forall s < |l^prime| - 1 end{align}$$

易得如下递归关系：

$$beta_t(s) = { begin{array}{l} (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} \ (beta_{t+1}(s)+beta_{t+1}(s+1) + beta_{t+1}(s+2)) y^t_{l^prime_s} end{array}$$

def backward(y, labels):
    T, V = y.shape
    L = len(labels)
    beta = np.zeros([T, L])

    # init
    beta[-1, -1] = y[-1, labels[-1]]
    beta[-1, -2] = y[-1, labels[-2]]

    for t in range(T - 2, -1, -1):
        for i in range(L):
            s = labels[i]

            a = beta[t + 1, i] 
            if i + 1 < L:
                a += beta[t + 1, i + 1]
            if i + 2 < L and s != 0 and s != labels[i + 2]:
                a += beta[t + 1, i + 2]

            beta[t, i] = a * y[t, s]

    return beta

beta = backward(y, labels)
print(beta)

1.6 梯度计算

下面，我们利用前向、后者计算的 $alpha$ 和 $beta$ 来计算梯度。

根据 $alpha$、$beta$ 的定义，我们有：

$$alpha_t(s)beta_t(s) = sum_{pi in mathcal{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 mathcal{B}^{-1}(l):pi_t=l_s^prime} prod_{t=1}^T y^t_{pi_t}$$
则
$$frac{alpha_t(s)beta_t(s)}{ y^t_{l_s^prime}} = sum_{pi in mathcal{B}^{-1}(l):pi_t=l_s^prime} prod_{t=1}^T y^t_{pi_t} = sum_{pi in mathcal{B}^{-1}(l):pi_t=l_s^prime} p(pi|x)$$
于是，可得似然
$$p(l|x) = sum_{s=1}^{|l^prime|} sum_{pi in mathcal{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}}$$

为计算 $frac{partial p(l|x)}{partial y^t_k}$，观察上式右端求各项，仅有 $s=k$ 的项包含 $y^t_k$，因此，其他项的偏导都为零，不用考虑。于是有：

$$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}$$

利用除法的求导准则有：

$$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}$$

求导中，分子第一项是因为 $alpha(k)beta(k)$ 中包含为两个 $y^t_k$ 乘积项（即 ${y^t_k}^2$），其他均为与 $y^t_k$ 无关的常数。

$l$ 中可能包含多个 $k$ 字符，它们计算的梯度要进行累加，因此，最后的梯度计算结果为：

$$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)$$
其中， $lab(s)={s: l_s^prime = k}$。

一般我们优化似然函数的对数，因此，梯度计算如下：

$$frac{partial ln(p(l|x))}{partial y^t_k} =frac{1}{p(l|x)} frac{partial p(l|x)}{partial y^t_k}$$
其中，似然值在前向计算中已经求得： $p(l|x) = alpha_T(|l^prime|) + alpha_T(|l^prime|-1)$。

对于给定训练集 $D$，待优化的目标函数为：

$$O(D, N_w) = -sum_{(x,z)in D} ln(p(z|x))$$

得到梯度后，我们可以利用任意优化方法（e.g. SGD， Adam）进行训练。

def gradient(y, labels):
    T, V = y.shape
    L = len(labels)

    alpha = forward(y, labels)
    beta = backward(y, labels)
    p = alpha[-1, -1] + alpha[-1, -2]

    grad = np.zeros([T, V])
    for t in range(T):
        for s in range(V):
            lab = [i for i, c in enumerate(labels) if c == s]
            for i in lab:
                grad[t, s] += alpha[t, i] * beta[t, i] 
            grad[t, s] /= y[t, s] ** 2

    grad /= p
    return grad

grad = gradient(y, labels)
print(grad)

将基于前向-后向算法得到梯度与基于数值的梯度比较，以验证实现的正确性。

def check_grad(y, labels, w=-1, v=-1, toleration=1e-3):
    grad_1 = gradient(y, labels)[w, v]

    delta = 1e-10
    original = y[w, v]

    y[w, v] = original + delta
    alpha = forward(y, labels)
    log_p1 = np.log(alpha[-1, -1] + alpha[-1, -2])

    y[w, v] = original - delta
    alpha = forward(y, labels)
    log_p2 = np.log(alpha[-1, -1] + alpha[-1, -2])

    y[w, v] = original

    grad_2 = (log_p1 - log_p2) / (2 * delta)
    if np.abs(grad_1 - grad_2) > toleration:
        print('[%d, %d]：%.2e' % (w, v, np.abs(grad_1 - grad_2)))

for toleration in [1e-5, 1e-6]:
    print('%.e' % toleration)
    for w in range(y.shape[0]):
        for v in range(y.shape[1]):
            check_grad(y, labels, w, v, toleration)

1e-05
1e-06
[0, 3]：3.91e-06
[1, 0]：3.61e-06
[1, 3]：2.66e-06
[2, 0]：2.67e-06
[2, 3]：3.88e-06
[3, 0]：4.71e-06
[3, 3]：3.39e-06
[4, 0]：1.24e-06
[4, 3]：4.79e-06
[5, 0]：1.57e-06
[5, 3]：2.98e-06
[6, 0]：5.03e-06
[6, 3]：4.89e-06
[7, 0]：1.05e-06
[7, 4]：4.19e-06
[8, 4]：5.57e-06
[9, 0]：5.95e-06
[9, 3]：3.85e-06
[10, 0]：1.09e-06
[10, 3]：1.53e-06
[10, 4]：3.82e-06

可以看到，前向-后向及数值梯度两种方法计算的梯度差异都在 1e-5 以下，误差最多在 1e-6 的量级。这初步验证了前向-后向梯度计算方法原理和实现的正确性。

1.7 logits 梯度

在实际训练中，为了计算方便，可以将 CTC 和 softmax 的梯度计算合并，公式如下：

$$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)$$

这是因为，softmax 的梯度反传公式为：

$$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)$$

接合上面两式，有:

$$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$$

def gradient_logits_naive(y, labels):
    '''
    gradient by back propagation
    '''
    y_grad = gradient(y, labels)

    sum_y_grad = np.sum(y_grad * y, axis=1, keepdims=True)
    u_grad = y * (y_grad - sum_y_grad) 

    return u_grad

def gradient_logits(y, labels):
    '''
    '''
    T, V = y.shape
    L = len(labels)

    alpha = forward(y, labels)
    beta = backward(y, labels)
    p = alpha[-1, -1] + alpha[-1, -2]

    u_grad = np.zeros([T, V])
    for t in range(T):
        for s in range(V):
            lab = [i for i, c in enumerate(labels) if c == s]
            for i in lab:
                u_grad[t, s] += alpha[t, i] * beta[t, i] 
            u_grad[t, s] /= y[t, s] * p

    u_grad -= y
    return u_grad

grad_l = gradient_logits_naive(y, labels)
grad_2 = gradient_logits(y, labels)

print(np.sum(np.abs(grad_l - grad_2)))

1.34961486431e-15

同上，我们利用数值梯度来初步检验梯度计算的正确性：

def check_grad_logits(x, labels, w=-1, v=-1, toleration=1e-3):
    grad_1 = gradient_logits(softmax(x), labels)[w, v]

    delta = 1e-10
    original = x[w, v]

    x[w, v] = original + delta
    y = softmax(x)
    alpha = forward(y, labels)
    log_p1 = np.log(alpha[-1, -1] + alpha[-1, -2])

    x[w, v] = original - delta
    y = softmax(x)
    alpha = forward(y, labels)
    log_p2 = np.log(alpha[-1, -1] + alpha[-1, -2])

    x[w, v] = original

    grad_2 = (log_p1 - log_p2) / (2 * delta)
    if np.abs(grad_1 - grad_2) > toleration:
        print('[%d, %d]：%.2e, %.2e, %.2e' % (w, v, grad_1, grad_2, np.abs(grad_1 - grad_2)))

np.random.seed(1111)
x = np.random.random([10, 10])
for toleration in [1e-5, 1e-6]:
    print('%.e' % toleration)
    for w in range(x.shape[0]):
        for v in range(x.shape[1]):
            check_grad_logits(x, labels, w, v, toleration)

2. 数值稳定性

CTC 的训练过程面临数值下溢的风险，特别是序列较大的情况下。下面介绍两种数值上稳定的工程优化方法：1）log 域（许多 CRF 实现的常用方法）；2）scale 技巧（原始论文 [1] 使用的方法）。

2.1 log 域计算

log 计算涉及 logsumexp 操作。
经验表明，在 log 域计算，即使使用单精度，也表现出良好的数值稳定性，可以有效避免下溢的风险。稳定性的代价是增加了运算的复杂性——原始实现只涉及乘加运算，log 域实现则需要对数和指数运算。

ninf = -np.float('inf')

def _logsumexp(a, b):
    '''
    np.log(np.exp(a) + np.exp(b))

    '''

    if a < b:
        a, b = b, a

    if b == ninf:
        return a
    else:
        return a + np.log(1 + np.exp(b - a)) 

def logsumexp(*args):
    '''
    from scipy.special import logsumexp
    logsumexp(args)
    '''
    res = args[0]
    for e in args[1:]:
        res = _logsumexp(res, e)
    return res

2.1.1 log 域前向算法

基于 log 的前向算法实现如下：

def forward_log(log_y, labels):
    T, V = log_y.shape
    L = len(labels)
    log_alpha = np.ones([T, L]) * ninf

    # init
    log_alpha[0, 0] = log_y[0, labels[0]]
    log_alpha[0, 1] = log_y[0, labels[1]]

    for t in range(1, T):
        for i in range(L):
            s = labels[i]

            a = log_alpha[t - 1, i]
            if i - 1 >= 0:
                a = logsumexp(a, log_alpha[t - 1, i - 1])
            if i - 2 >= 0 and s != 0 and s != labels[i - 2]:
                a = logsumexp(a, log_alpha[t - 1, i - 2])

            log_alpha[t, i] = a + log_y[t, s]

    return log_alpha

log_alpha = forward_log(np.log(y), labels)
alpha = forward(y, labels)
print(np.sum(np.abs(np.exp(log_alpha) - alpha)))

8.60881935942e-17

2.1.2 log 域后向算法

基于 log 的后向算法实现如下：

def backward_log(log_y, labels):
    T, V = log_y.shape
    L = len(labels)
    log_beta = np.ones([T, L]) * ninf

    # init
    log_beta[-1, -1] = log_y[-1, labels[-1]]
    log_beta[-1, -2] = log_y[-1, labels[-2]]

    for t in range(T - 2, -1, -1):
        for i in range(L):
            s = labels[i]

            a = log_beta[t + 1, i] 
            if i + 1 < L:
                a = logsumexp(a, log_beta[t + 1, i + 1])
            if i + 2 < L and s != 0 and s != labels[i + 2]:
                a = logsumexp(a, log_beta[t + 1, i + 2])

            log_beta[t, i] = a + log_y[t, s]

    return log_beta

log_beta = backward_log(np.log(y), labels)
beta = backward(y, labels)
print(np.sum(np.abs(np.exp(log_beta) - beta)))

1.10399945005e-16

2.1.3 log 域梯度计算

在前向、后向基础上，也可以在 log 域上计算梯度。

def gradient_log(log_y, labels):
    T, V = log_y.shape
    L = len(labels)

    log_alpha = forward_log(log_y, labels)
    log_beta = backward_log(log_y, labels)
    log_p = logsumexp(log_alpha[-1, -1], log_alpha[-1, -2])

    log_grad = np.ones([T, V]) * ninf
    for t in range(T):
        for s in range(V):
            lab = [i for i, c in enumerate(labels) if c == s]
            for i in lab:
                log_grad[t, s] = logsumexp(log_grad[t, s], log_alpha[t, i] + log_beta[t, i]) 
            log_grad[t, s] -= 2 * log_y[t, s]

    log_grad -= log_p
    return log_grad

log_grad = gradient_log(np.log(y), labels)
grad = gradient(y, labels)
#print(log_grad)
#print(grad)
print(np.sum(np.abs(np.exp(log_grad) - grad)))

4.97588081849e-14

2.2 scale

2.2.1 前向算法

为了避免下溢，在前向算法的每个时刻，都对计算出的 $alpha$ 的范围进行缩放：

$$C_t stackrel{def}{=} sum_salpha_t(s)$$

$$hat{alpha}_t = frac{alpha_t(s)}{C_t}$$

缩放后的 $alpha$，不会随着时刻的积累变得太小。$hat{alpha}$ 替代 $alpha$，进行下一时刻的迭代。

def forward_scale(y, labels):
    T, V = y.shape
    L = len(labels)
    alpha_scale = np.zeros([T, L])

    # init
    alpha_scale[0, 0] = y[0, labels[0]]
    alpha_scale[0, 1] = y[0, labels[1]]
    Cs = []

    C = np.sum(alpha_scale[0])
    alpha_scale[0] /= C
    Cs.append(C)

    for t in range(1, T):
        for i in range(L):
            s = labels[i]

            a = alpha_scale[t - 1, i] 
            if i - 1 >= 0:
                a += alpha_scale[t - 1, i - 1]
            if i - 2 >= 0 and s != 0 and s != labels[i - 2]:
                a += alpha_scale[t - 1, i - 2]

            alpha_scale[t, i] = a * y[t, s]

        C = np.sum(alpha_scale[t])
        alpha_scale[t] /= C
        Cs.append(C)

    return alpha_scale, Cs

由于进行了缩放，最后计算概率时要时行补偿：

$$p(l|x) = alpha_T(|l^prime|) + alpha_T(|l^prime|-1) = (hatalpha_T(|l^prime|) + hatalpha_T(|l^prime|-1) * prod_{t=1}^T C_t$$

$$ln(p(l|x)) = sum_t^Tln(C_t) + ln(hatalpha_T(|l^prime|) + hatalpha_T(|l^prime|-1))$$

labels = [0, 1, 2, 0]  # 0 for blank

alpha_scale, Cs = forward_scale(y, labels)
log_p = np.sum(np.log(Cs)) + np.log(alpha_scale[-1][labels[-1]] + alpha_scale[-1][labels[-2]])

alpha = forward(y, labels)
p = alpha[-1, labels[-1]] + alpha[-1, labels[-2]]

print(np.log(p), log_p, np.log(p) - log_p)

(-13.202925982240107, -13.202925982240107, 0.0)

2.2.2 后向算法

后向算法缩放类似于前向算法，公式如下：

$$D_t stackrel{def}{=} sum_sbeta_t(s)$$

$$hat{beta}_t = frac{beta_t(s)}{D_t}$$

def backward_scale(y, labels):
    T, V = y.shape
    L = len(labels)
    beta_scale = np.zeros([T, L])

    # init
    beta_scale[-1, -1] = y[-1, labels[-1]]
    beta_scale[-1, -2] = y[-1, labels[-2]]

    Ds = []

    D = np.sum(beta_scale[-1,:])
    beta_scale[-1] /= D
    Ds.append(D)

    for t in range(T - 2, -1, -1):
        for i in range(L):
            s = labels[i]

            a = beta_scale[t + 1, i] 
            if i + 1 < L:
                a += beta_scale[t + 1, i + 1]
            if i + 2 < L and s != 0 and s != labels[i + 2]:
                a += beta_scale[t + 1, i + 2]

            beta_scale[t, i] = a * y[t, s]

        D = np.sum(beta_scale[t])
        beta_scale[t] /= D
        Ds.append(D)

    return beta_scale, Ds[::-1]

beta_scale, Ds = backward_scale(y, labels)
print(beta_scale)

2.2.3 梯度计算

$$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)$$

考虑到

$$p(l|x) = sum_{s=1}^{|l^prime|} frac{alpha_t(s)beta_t(s)}{ y^t_{l_s^prime}}$$
以及
$$alpha_t(s) = hatalpha_t(s) cdot prod_{k=1}^t C_k$$

$$beta_t(s) = hatbeta_t(s) cdot prod_{k=t}^T D_k$$

$$frac{partial ln(p(l|x))}{partial y^t_k} = frac{1}{sum_{s=1}^{|l^prime|} frac{hatalpha_t(s) hatbeta_t(s)}{y^t_{l^prime_s}}} frac{1}{{y^t_k}^2} sum_{s in lab(l, k)} hatalpha_t(s) hatbeta_t(s)$$

式中最右项中的各个部分我们都已经求得。梯度计算实现如下：

def gradient_scale(y, labels):
    T, V = y.shape
    L = len(labels)

    alpha_scale, _ = forward_scale(y, labels)
    beta_scale, _ = backward_scale(y, labels)

    grad = np.zeros([T, V])
    for t in range(T):
        for s in range(V):
            lab = [i for i, c in enumerate(labels) if c == s]
            for i in lab:
                grad[t, s] += alpha_scale[t, i] * beta_scale[t, i]
            grad[t, s] /= y[t, s] ** 2

        # normalize factor
        z = 0
        for i in range(L):
            z += alpha_scale[t, i] * beta_scale[t, i] / y[t, labels[i]]
        grad[t] /= z

    return grad

labels = [0, 3, 0, 3, 0, 4, 0]  # 0 for blank
grad_1 = gradient_scale(y, labels)
grad_2 = gradient(y, labels)
print(np.sum(np.abs(grad_1 - grad_2)))

6.86256607096e-15

2.2.4 logits 梯度

类似于 y 梯度的推导，logits 梯度计算公式如下：

$$frac{partial ln(p(l|x))}{partial u^t_k} = frac{1}{y^t_k Z_t} sum_{s in lab(l, k)} hatalpha_t(s)hatbeta_t(s) - y^t_k$$
其中，
$$Z_t stackrel{def}{=} sum_{s=1}^{|l^prime|} frac{hatalpha_t(s)hatbeta_t(s)}{y^t_{l^prime_s}}$$

3. 解码

训练和的 $N_w$ 可以用来预测新的样本输入对应的输出字符串，这涉及到解码。
按照最大似然准则，最优的解码结果为：

$$h(x) = underset{l in L^{le T}}{mathrm{argmax}} p(l|x)$$

然而，上式不存在已知的高效解法。下面介绍几种实用的近似破解码方法。

3.1 贪心搜索 （greedy search）

虽然 $p(l|x)$ 难以有效的计算，但是由于 CTC 的独立性假设，对于某个具体的字符串 $pi$（去 blank 前），确容易计算：

$$p(pi|x) = prod_{k=1}^T p(pi_k|x)$$

因此，我们放弃寻找使 $p(l|x)$ 最大的字符串，退而寻找一个使 $p(pi|x)$ 最大的字符串，即：

$$h(x) approx B(pi^star)$$
其中，
$$pi^star = underset{pi in N^T}{mathrm{argmax}} p(pi|x)$$

简化后，解码过程（构造 $pi^star$）变得非常简单（基于独立性假设）： 在每个时刻输出概率最大的字符:

$$pi^star = cat_{t=1}^T(underset{s in L^prime}{mathrm{argmax}} y^t_s)$$

def remove_blank(labels, blank=0):
    new_labels = []

    # combine duplicate
    previous = None
    for l in labels:
        if l != previous:
            new_labels.append(l)
            previous = l

    # remove blank     
    new_labels = [l for l in new_labels if l != blank]

    return new_labels

def insert_blank(labels, blank=0):
    new_labels = [blank]
    for l in labels:
        new_labels += [l, blank]
    return new_labels

def greedy_decode(y, blank=0):
    raw_rs = np.argmax(y, axis=1)
    rs = remove_blank(raw_rs, blank)
    return raw_rs, rs

np.random.seed(1111)
y = softmax(np.random.random([20, 6]))
rr, rs = greedy_decode(y)
print(rr)
print(rs)

[1 3 5 5 5 5 1 5 3 4 4 3 0 4 5 0 3 1 3 3]
[1, 3, 5, 1, 5, 3, 4, 3, 4, 5, 3, 1, 3]

3.2 束搜索（Beam Search）

显然，贪心搜索的性能非常受限。例如，它不能给出除最优路径之外的其他其优路径。很多时候，如果我们能拿到 nbest 的路径，后续可以利用其他信息来进一步优化搜索的结果。束搜索能近似找出 top 最优的若干条路径。

def beam_decode(y, beam_size=10):
    T, V = y.shape
    log_y = np.log(y)

    beam = [([], 0)]
    for t in range(T):  # for every timestep
        new_beam = []
        for prefix, score in beam:
            for i in range(V):  # for every state
                new_prefix = prefix + [i]
                new_score = score + log_y[t, i]

                new_beam.append((new_prefix, new_score))

        # top beam_size
        new_beam.sort(key=lambda x: x[1], reverse=True)
        beam = new_beam[:beam_size]

    return beam

np.random.seed(1111)
y = softmax(np.random.random([20, 6]))
beam = beam_decode(y, beam_size=100)
for string, score in beam[:20]:
    print(remove_blank(string), score)

3.3 前缀束搜索（Prefix Beam Search）

直接的束搜索的一个问题是，在保存的 top N 条路径中，可能存在多条实际上是同一结果（经过去重复、去 blank 操作）。这减少了搜索结果的多样性。下面介绍的前缀搜索方法，在搜索过程中不断的合并相同的前缀[2]。参考 gist，前缀束搜索实现如下：

from collections import defaultdict

def prefix_beam_decode(y, beam_size=10, blank=0):
    T, V = y.shape
    log_y = np.log(y)

    beam = [(tuple(), (0, ninf))]  # blank, non-blank
    for t in range(T):  # for every timestep
        new_beam = defaultdict(lambda : (ninf, ninf))

        for prefix, (p_b, p_nb) in beam:
            for i in range(V):  # for every state
                p = log_y[t, i]

                if i == blank:  # propose a blank
                    new_p_b, new_p_nb = new_beam[prefix]
                    new_p_b = logsumexp(new_p_b, p_b + p, p_nb + p)
                    new_beam[prefix] = (new_p_b, new_p_nb)
                    continue
                else:  # extend with non-blank
                    end_t = prefix[-1] if prefix else None

                    # exntend current prefix
                    new_prefix = prefix + (i,)
                    new_p_b, new_p_nb = new_beam[new_prefix]
                    if i != end_t:
                        new_p_nb = logsumexp(new_p_nb, p_b + p, p_nb + p)
                    else:
                        new_p_nb = logsumexp(new_p_nb, p_b + p)
                    new_beam[new_prefix] = (new_p_b, new_p_nb)

                    # keep current prefix
                    if i == end_t:
                        new_p_b, new_p_nb = new_beam[prefix]
                        new_p_nb = logsumexp(new_p_nb, p_nb + p)
                        new_beam[prefix] = (new_p_b, new_p_nb)

        # top beam_size
        beam = sorted(new_beam.items(), key=lambda x : logsumexp(*x[1]), reverse=True)
        beam = beam[:beam_size]

    return beam

np.random.seed(1111)
y = softmax(np.random.random([20, 6]))
beam = prefix_beam_decode(y, beam_size=100)
for string, score in beam[:20]:
    print(remove_blank(string), score)

3.4 其他解码方法

上述介绍了基本解码方法。实际中，搜索过程可以接合额外的信息，提高搜索的准确度（例如在语音识别任务中，加入语言模型得分信息, 前缀束搜索+语言模型）。

本质上，CTC 只是一个训练准则。训练完成后，$N_w$ 输出一系列概率分布，这点和常规基于交叉熵准则训练的模型完全一致。因此，特定应用领域常规的解码也可以经过一定修改用来 CTC 的解码。例如在语音识别任务中，利用 CTC 训练的声学模型可以无缝融入原来的 WFST 的解码器中[5]（e.g. 参见 EESEN）。

此外，[1] 给出了一种利用 CTC 顶峰特点的启发式搜索方法。

4. 工具

warp-ctc 是百度开源的基于 CPU 和 GPU 的高效并行实现。warp-ctc 自身提供 C 语言接口，对于流利的机器学习工具（torch、 pytorch 和 tensorflow、chainer）都有相应的接口绑定。

cudnn 7 以后开始提供 CTC 支持。

Tensorflow 也原生支持 CTC loss，及 greedy 和 beam search 解码器。

小结

1. CTC 可以建模无对齐信息的多对多序列问题（输入长度不小于输出），如语音识别、连续字符识别 [3,4]。
2. CTC 不需要输入与输出的对齐信息，可以实现端到端的训练。
3. CTC 在 loss 的计算上，利用了整个 labels 序列的全局信息，某种意义上相对逐帧计算损失的方法，”更加区分性”。

References

1. Graves et al. Connectionist Temporal Classification : Labelling Unsegmented Sequence Data with Recurrent Neural Networks.
2. Hannun et al. First-Pass Large Vocabulary Continuous Speech Recognition using Bi-Directional Recurrent DNNs.
3. Graves et al. Towards End-To-End Speech Recognition with Recurrent Neural Networks.
4. Liwicki et al. A novel approach to on-line handwriting recognition based on bidirectional long short-term memory networks.
5. Zenkel et al. Comparison of Decoding Strategies for CTC Acoustic Models.
6. Huannun. Sequence Modeling with CTC.

最后

以上就是调皮雪碧为你收集整理的【Learning Notes】CTC 原理及实现1. 算法2. 数值稳定性3. 解码4. 工具小结References的全部内容，希望文章能够帮你解决【Learning Notes】CTC 原理及实现1. 算法2. 数值稳定性3. 解码4. 工具小结References所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。