【Learning Notes】CTC 原理及实现

原文

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

  jupyter notebook 版见 repo.

 

1. 算法

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

 

1.1 序列问题形式化。

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

 

Nw:(Rm)T→(Rn)TNw:(Rm)T→(Rn)T 
其中，序列目标为字符串（词表大小为 nn），即 NwNw 输出为 nn 维多项概率分布（e.g. 经过 softmax 处理）。

网络输出为：y=Nwy=Nw，其中，ytkykt tt 表示时刻第 kk 项的概率。

 
图1. 序列建模【src】

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

 

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

[[ 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.]]
123456789101112131415161718192021222324

 

1.2 align-free 变长映射

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

扩展原始词表 LL 为 L′=L∪{blank}L′=L∪{blank}。对输出字符串，定义操作 BB：1）合并连续的相同符号；2）去掉 blank 字符。

例如，对于 “aa%bb%%cc”，应用 BB，则实际上代表的是字符串 “abc”。同理“%a%b%cc%” 也同样代表 “abc”。 
B(aa%bb%%cc)=B(%a%b%cc%)=abcB(aa%bb%%cc)=B(%a%b%cc%)=abc

通过引入blank 及 BB，可以实现了变长的映射。 
L′T→L≤TL′T→L≤T

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

 

1.3 似然计算

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

给定输入 xx，输出 ll 的条件概率为： 
p(l|x)=∑π∈B−1(l)p(π|x)p(l|x)=∑π∈B−1(l)p(π|x)

其中，B−1(l)B−1(l) 表示了长度为 TT 且示经过 BB 结果为 ll 字符串的集合。

CTC 假设输出的概率是（相对于输入）条件独立的，因此有： 
p(π|x)=∏ytπt,∀π∈L′Tp(π|x)=∏yπtt,∀π∈L′T

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

 

1.4 前向算法

在前向及后向计算中，CTC 需要将输出字符串进行扩展。具体的，(a1,⋯,am)(a1,⋯,am) 每个字符之间及首尾分别插入 blank，即扩展为 (%,a1,%,a2,%,⋯,%,am,%)(%,a1,%,a2,%,⋯,%,am,%)。下面的 ll 为原始字符串，l′l′ 指为扩展后的字符串。

定义 
αt(s)=def∑π∈NT:B(π1:t)=l1:s∏t′=1tytπ′αt(s)=def∑π∈NT:B(π1:t)=l1:s∏t′=1tyπ′t

显然有，

 

α1(1)=y1b,α1(2)=y1l1,α1(s)=0,∀s>2α1(1)=yb1,α1(2)=yl11,α1(s)=0,∀s>2 
根据 αα 的定义，有如下递归关系： 
αt(s)={(αt−1(s)+αt−1(s−1))ytl′s,   if l′s=b or l′s−2=l′s(αt−1(s)+αt−1(s−1)+αt−1(s−2))ytl′s  otherwiseαt(s)={(αt−1(s)+αt−1(s−1))yls′t,   if ls′=b or ls−2′=ls′(αt−1(s)+αt−1(s−1)+αt−1(s−2))yls′t  otherwise

 

1.4.1 Case 2

递归公式中 case 2 是一般的情形。如图所示，tt 时刻字符为 ss 为 blank 时，它可能由于两种情况扩展而来：1）重复上一字符，即上个字符也是 a，2）字符发生转换，即上个字符是非 a 的字符。第二种情况又分为两种情形，2.1）上一字符是 blank；2.2）a 由非 blank 字符直接跳转而来（BB） 操作中， blank 最终会被去掉，因此 blank 并不是必须的）。 
 
图2. 前向算法 Case 2 示例【src】

 

1.4.2 Case 1

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

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

我们可以利用动态规划计算所有 αα 的值，算法时间和空间复杂度为 O(T∗L)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)1234567891011121314151617181920212223242526

最后可以得到似然 p(l|x)=αT(|l′|)+αT(|l′|−1)p(l|x)=αT(|l′|)+αT(|l′|−1)。

 

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

6.81811271177e-06
1

 

1.5 后向计算

类似于前向计算，我们定义后向计算。 
首先定义 
βt(s)=def∑π∈NT:B(πt:T)=ls:|l|∏t′=tTytπ′βt(s)=def∑π∈NT:B(πt:T)=ls:|l|∏t′=tTyπ′t

显然， 
βT(|l′|)=yTb,βT(|l′|−1)=yTl|l|,βT(s)=0,∀s<|l′|−1βT(|l′|)=ybT,βT(|l′|−1)=yl|l|T,βT(s)=0,∀s<|l′|−1

易得如下递归关系： 
βt(s)={(βt+1(s)+βt+1(s+1))ytl′s,   if l′s=b or l′s+2=l′s(βt+1(s)+βt+1(s+1)+βt+1(s+2))ytl′sβt(s)={(βt+1(s)+βt+1(s+1))yls′t,   if ls′=b or ls+2′=ls′(βt+1(s)+βt+1(s+1)+βt+1(s+2))yls′t

 

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

 

1.6 梯度计算

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

根据 αα、ββ 的定义，我们有： 
αt(s)βt(s)=∑π∈B−1(l):πt=l′sytl′s∏t=1Tytπt=ytl′s⋅∑π∈B−1(l):πt=l′s∏t=1Tytπtαt(s)βt(s)=∑π∈B−1(l):πt=ls′yls′t∏t=1Tyπtt=yls′t⋅∑π∈B−1(l):πt=ls′∏t=1Tyπtt 
则 
αt(s)βt(s)ytl′s=∑π∈B−1(l):πt=l′s∏t=1Tytπt=∑π∈B−1(l):πt=l′sp(π|x)αt(s)βt(s)yls′t=∑π∈B−1(l):πt=ls′∏t=1Tyπtt=∑π∈B−1(l):πt=ls′p(π|x) 
于是，可得似然 
p(l|x)=∑s=1|l′|∑π∈B−1(l):πt=l′sp(π|x)=∑s=1|l′|αt(s)βt(s)ytl′sp(l|x)=∑s=1|l′|∑π∈B−1(l):πt=ls′p(π|x)=∑s=1|l′|αt(s)βt(s)yls′t

为计算 ∂p(l|x)∂ytk∂p(l|x)∂ykt，观察上式右端求各项，仅有 s=ks=k 的项包含 ytkykt，因此，其他项的偏导都为零，不用考虑。于是有：

 

∂p(l|x)∂ytk=∂αt(k)βt(k)ytk∂ytk∂p(l|x)∂ykt=∂αt(k)βt(k)ykt∂ykt

利用除法的求导准则有： 
∂p(l|x)∂ytk=2⋅αt(k)βt(k)ytk⋅ytk−αt(k)βt(k)⋅1ytk2=αt(k)βt(k)ytk2∂p(l|x)∂ykt=2⋅αt(k)βt(k)ykt⋅ykt−αt(k)βt(k)⋅1ykt2=αt(k)βt(k)ykt2

  求导中，分子第一项是因为 α(k)β(k)α(k)β(k) 中包含为两个 ytkykt 乘积项（即 ytk2ykt2），其他均为与 ytkykt 无关的常数。

ll 中可能包含多个 kk 字符，它们计算的梯度要进行累加，因此，最后的梯度计算结果为： 
∂p(l|x)∂ytk=1ytk2∑s∈lab(l,k)αt(s)βt(s)∂p(l|x)∂ykt=1ykt2∑s∈lab(l,k)αt(s)βt(s) 
其中，lab(s)={s:l′s=k}lab(s)={s:ls′=k}。

一般我们优化似然函数的对数，因此，梯度计算如下： 
∂ln(p(l|x))∂ytk=1p(l|x)∂p(l|x)∂ytk∂ln⁡(p(l|x))∂ykt=1p(l|x)∂p(l|x)∂ykt 
其中，似然值在前向计算中已经求得： p(l|x)=αT(|l′|)+αT(|l′|−1)p(l|x)=αT(|l′|)+αT(|l′|−1)。

对于给定训练集 DD，待优化的目标函数为： 
O(D,Nw)=−∑(x,z)∈Dln(p(z|x))O(D,Nw)=−∑(x,z)∈Dln⁡(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)123456789101112131415161718192021

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

 

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

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
1234567891011121314151617181920212223

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

 

1.7 logits 梯度

在实际训练中，为了计算方便，可以将 CTC 和 softmax 的梯度计算合并，公式如下： 
∂ln(p(l|x))∂ytk=ytk−1ytk⋅p(l|x)∑s∈lab(l,k)αt(s)βt(s)∂ln⁡(p(l|x))∂ykt=ykt−1ykt⋅p(l|x)∑s∈lab(l,k)αt(s)βt(s)

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

 

∂ln(p(l|x))∂utk=ytk(∂ln(p(l|x))∂ytk−∑j=1V∂ln(p(l|x))∂ytjytj)∂ln⁡(p(l|x))∂ukt=ykt(∂ln⁡(p(l|x))∂ykt−∑j=1V∂ln⁡(p(l|x))∂yjtyjt)

接合上面两式，有: 
∂ln(p(l|x))∂utk=1ytkp(l|x)∑s∈lab(l,k)αt(s)βt(s)−ytk∂ln⁡(p(l|x))∂ukt=1yktp(l|x)∑s∈lab(l,k)αt(s)βt(s)−ykt

 

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

1.34961486431e-15
1

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

 

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

 

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 res12345678910111213141516171819202122232425

 

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

8.60881935942e-17
1

 

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

1.10399945005e-16
1

 

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

4.97588081849e-14
1

 

2.2 scale

 

2.2.1 前向算法

为了避免下溢，在前向算法的每个时刻，都对计算出的 αα 的范围进行缩放： 
Ct=def∑sαt(s)Ct=def∑sαt(s)

 

α^t=αt(s)Ctα^t=αt(s)Ct

缩放后的 αα，不会随着时刻的积累变得太小。α^α^ 替代 αα，进行下一时刻的迭代。

 

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, Cs12345678910111213141516171819202122232425262728293031

由于进行了缩放，最后计算概率时要时行补偿： 
p(l|x)=αT(|l′|)+αT(|l′|−1)=(α^T(|l′|)+α^T(|l′|−1)∗∏t=1TCtp(l|x)=αT(|l′|)+αT(|l′|−1)=(α^T(|l′|)+α^T(|l′|−1)∗∏t=1TCt

 

ln(p(l|x))=∑tTln(Ct)+ln(α^T(|l′|)+α^T(|l′|−1))ln⁡(p(l|x))=∑tTln⁡(Ct)+ln⁡(α^T(|l′|)+α^T(|l′|−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)123456789

(-13.202925982240107, -13.202925982240107, 0.0)
1

 

2.2.2 后向算法

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

 

Dt=def∑sβt(s)Dt=def∑sβt(s)

 

β^t=βt(s)Dtβ^t=βt(s)Dt

 

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

 

2.2.3 梯度计算

∂ln(p(l|x))∂ytk=1p(l|x)∂p(l|x)∂ytk=1p(l|x)1ytk2∑s∈lab(l,k)αt(s)βt(s)∂ln⁡(p(l|x))∂ykt=1p(l|x)∂p(l|x)∂ykt=1p(l|x)1ykt2∑s∈lab(l,k)αt(s)βt(s)

考虑到 
p(l|x)=∑s=1|l′|αt(s)βt(s)ytl′sp(l|x)=∑s=1|l′|αt(s)βt(s)yls′t 
以及 
αt(s)=α^t(s)⋅∏k=1tCkαt(s)=α^t(s)⋅∏k=1tCk

 

βt(s)=β^t(s)⋅∏k=tTDkβt(s)=β^t(s)⋅∏k=tTDk

 

∂ln(p(l|x))∂ytk=1∑|l′|s=1α^t(s)β^t(s)ytl′s1ytk2∑s∈lab(l,k)α^t(s)β^t(s)∂ln⁡(p(l|x))∂ykt=1∑s=1|l′|α^t(s)β^t(s)yls′t1ykt2∑s∈lab(l,k)α^t(s)β^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)))123456789101112131415161718192021222324252627

6.86256607096e-15
1

 

2.2.4 logits 梯度

类似于 y 梯度的推导，logits 梯度计算公式如下： 
∂ln(p(l|x))∂utk=1ytkZt∑s∈lab(l,k)α^t(s)β^t(s)−ytk∂ln⁡(p(l|x))∂ukt=1yktZt∑s∈lab(l,k)α^t(s)β^t(s)−ykt 
其中， 
Zt=def∑s=1|l′|α^t(s)β^t(s)ytl′sZt=def∑s=1|l′|α^t(s)β^t(s)yls′t

 

3. 解码

训练和的 NwNw 可以用来预测新的样本输入对应的输出字符串，这涉及到解码。 
按照最大似然准则，最优的解码结果为： 
h(x)=argmaxl∈L≤T p(l|x)h(x)=argmaxl∈L≤T p(l|x)

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

 

3.1 贪心搜索 （greedy search）

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

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

 

h(x)≈B(π⋆)h(x)≈B(π⋆) 
其中， 
π⋆=argmaxπ∈NT p(π|x)π⋆=argmaxπ∈NT p(π|x)

简化后，解码过程（构造 π⋆π⋆）变得非常简单（基于独立性假设）： 在每个时刻输出概率最大的字符: 
π⋆=catTt=1(argmaxs∈L′ yts)π⋆=catt=1T(argmaxs∈L′ yst)

 

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

[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]
12

 

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

 

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

 

3.4 其他解码方法

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

本质上，CTC 只是一个训练准则。训练完成后，NwNw 输出一系列概率分布，这点和常规基于交叉熵准则训练的模型完全一致。因此，特定应用领域常规的解码也可以经过一定修改用来 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 解码器。

 

小结

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

 

References

Graves et al. Connectionist Temporal Classification : Labelling Unsegmented Sequence Data with Recurrent Neural Networks.
Hannun et al. First-Pass Large Vocabulary Continuous Speech Recognition using Bi-Directional Recurrent DNNs.
Graves et al. Towards End-To-End Speech Recognition with Recurrent Neural Networks.
Liwicki et al. A novel approach to on-line handwriting recognition based on bidirectional long short-term memory networks.
Zenkel et al. Comparison of Decoding Strategies for CTC Acoustic Models.
Huannun. Sequence Modeling with CTC.
————————————————
版权声明：本文为CSDN博主「MoussaTintin」的原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/JackyTintin/article/details/79425866