ctc_loss 公式推导与C++实现

CTC简介

本文不介绍CTC的背景与具体应用，有关CTC的基本知识可以通过以下文章了解:
Hannun Awni经典博客.

总体思路

初始时我们得到一个输入矩阵，它的行代表时间步，长度为T,列代表不同的字符在不同时间片下的概率，长度为alphabet_size(包含空字符)。同时我们有一个lebel序列，它代表我们的正确的输出。
我们需要做的是根据label序列利用该矩阵前向计算一遍概率，并把每个字符在每个时间步下的概率保存在矩阵α(t,u)中。然后反向计算一遍概率，存储在矩阵β(t,u)中。最后利用我们推导出的公式计算每一个输入的梯度。
在本文中，我们的实现只以一个样本为例。

准备工作

1. softmax

我们初始得到的输入矩阵并不表示概率，需要对该矩阵的每一行进行softmax。
softmax公式:$S_i=\frac{e^{v_i}}{{\sum }_j{e}^{v_j}}$

2.label的扩展

CTC的算法要求我们对原laebl进行扩展，具体是在laebl的开头，结尾以及每两个字符中间加上空符号。如原label为 {1,3,5,5,7} 则扩展后的符号为:
{0,1,0,3,0,5,0,5,0,7,0} 注:这里的数字对应的符号在符号表中的下标，0默认为空字符。

3.取对数

由于在CTC_loss的计算过程中有大量的概率相乘，这些浮点数往往较小，极有可能发生underflow，为了数值的稳定，我们往往先对这些概率取对数,将概率的加法与乘法转化为对数的计算:
概率乘法转化:$\log (a\cdot b)=\log (a)+\log (b)$
概率加法转化:$\log (a+b)=\log (1+e^{(\log (b)-\log (a))}）$

需要特别注意，由于log(0)无法表示，
而我们在计算中会涉及大量的log(0)的运算，
所以我们需要定义一个-INFINITY,
在初始化时需要将所有初始化为-INFINITY,它表示log(0)。
double log_sum(double a, double b)
{
	//log(a+b)=log(a)+log(1+exp(logb-loga));
	if (a == -INFINITY)
		return b;
	if (b == -INFINITY)
		return a;

	return a + log(1 + exp(b - a));
}

double log_prod(double a,double b)
{
	if (a == -INFINITY && b == -INFINITY)
		return -INFINITY;
	if (a == -INFINITY)
		return b;
	if (b == -INFINITY)
		return a;

	return a + b;
}

4.相关代码

在具体实现中，我主要参考了百度的工作，他们在16年的时候开源了warp_ctc的代码: warp_ctc.

数据结构(类)设计:

class cpu_ctc {
public:
	cpu_ctc(vector<int> labels, vector<vector<double>> Prob);
	void compute();
	double compute_alphas();
	double compute_betas();
	void compute_gradients();
	
private:
	int* s_inc;
	int* e_inc;
	int blank_label = 0;   //空字符在字符表里对应的序号 
						 //默认为0; 
	double prob_;        //前向传播计算出的概率(取对数);
	int L,S,T;
	//L表示标签序列的长度，S表示扩展后的序列长度，T表示输入的时间步数量。
	int repeats;  //表示标签序列中紧接着的重复的标签数量。
	
	vector<vector<double>> alphas;  //计算前向概率; 
	vector<vector<double>> betas;   //计算后向概率
	vector<vector<double>> gradients;   //计算梯度; 
	vector<vector<double>> Prob;
	vector<vector<double>> tmp_alphas;
	vector<vector<double>> tmp_betas;
	vector<vector<double>> tmp_gradients;

	int* labels_with_blanks;
	//注:出现repeats后相当于原label需要扩充一个空! 
	int prepare(vector<int> labels, vector<vector<double>> Prob);
};

未加注释的设计的具体意义将在后文中详细介绍。

初始化工作:

//在构造函数中进行一些数据的初始化工作;
//写的较为朴素;
cpu_ctc::cpu_ctc(vector<int> labels, vector<vector<double>> Prob)
{
	for (int i = 0; i < Prob.size(); i++)  //初始化输入概率
	{
		vector<double> single;
		for (int j = 0; j < Prob[i].size(); j++)
			single.push_back(Prob[i][j]);
		this->Prob.push_back(single);
	}

	for (int i = 0; i < Prob.size(); i++) //alpha数组和beta数组
	//为正常的概率计算数组，全初始化为0.
	{
		vector<double> tmp;
		for (int j = 0; j < 2 * labels.size() + 1; j++)
			tmp.push_back(0);

		this->alphas.push_back(tmp);
		this->betas.push_back(tmp);
		//this->gradients.push_back(tmp);
	}
	
	//tmp_alphas和tmp_betas为log域下概率数组，全部初始化为-INFINITY
	for (int i = 0; i < Prob.size(); i++)
	{
		vector<double> tmp;
		for (int j = 0; j < 2 * labels.size() + 1; j++)
			tmp.push_back(-INFINITY);

		this->tmp_alphas.push_back(tmp);
		this->tmp_betas.push_back(tmp);
		//this->gradients.push_back(tmp);
	}
//初始化梯度数组;
	for (int i = 0; i < Prob.size(); i++)
	{
		vector<double> tmp;
		vector<double> tmp_;
		for (int j = 0; j < Prob[0].size(); j++)
		{
			tmp.push_back(0);
			tmp_.push_back(-INFINITY);
		}
			
		this->tmp_gradients.push_back(tmp_);
		this->gradients.push_back(tmp);
	}
	L = labels.size();
	T = Prob.size();
	repeats = this->prepare(labels, Prob);
}

prepare函数:
prepare函数的目的主要有两个：
1.生成s_int和e_int数组; 这两个数组在计算前向和反向概率中有用。
2.用blank符号扩充labels;

//主要模仿的是百度warp-ctc的开源代码；
int cpu_ctc::prepare(vector<int> labels, vector<vector<double>> Prob)
{
	S = 2 * labels.size() + 1;
	s_inc = (int*)malloc(sizeof(int) * (labels.size() * 2 + 1));
	e_inc = (int*)malloc(sizeof(int) * (labels.size() * 2 + 1));
	labels_with_blanks = (int*)malloc(sizeof(int) * (labels.size() * 2 + 1));

	int s_counter = 0;
	int e_counter = 0;

	this->s_inc[s_counter++] = 1;  //remain=0,start+=1;   
	int repeats = 0;

	for (int i = 1; i < labels.size(); i++)
	{
		if (labels[i - 1] == labels[i])  //repeat; 
		{
			repeats++;
			s_inc[s_counter++] = 1;
			s_inc[s_counter++] = 1;
			e_inc[e_counter++] = 1;
			e_inc[e_counter++] = 1;
		}
		else
		{
			this->s_inc[s_counter++] = 2;
			this->e_inc[e_counter++] = 2;
		}
	}
	this->e_inc[e_counter++] = 1;

	for (int i = 0; i < labels.size(); i++)
	{
		this->labels_with_blanks[2 * i] = this->blank_label;
		this->labels_with_blanks[2 * i + 1] = labels[i];
	}

	labels_with_blanks[2 * labels.size()] = this->blank_label;
	return repeats;
}

具体理解可以参考这两篇文章:
CTC实现——compute ctc loss（1）.
CTC实现——compute ctc loss（2）.

前向计算

原理

前向传播主要有两个任务，一是在给定输入x的情况下，计算得到label l的概率，即$p(l|x)$。而一个序列$l$可能通过多条路径映射得到，随着$l$长度的增加，相应路径的数量成指数增加，所以我们需要一种高效的算法帮助我们计算。
设原label的长度为U,经过扩展后长度变为2U+1,设扩展后的序列为$l^{\prime }$。对于一个特定的序列$l$,我们定义前向变量$$\alpha (t,u)=\sum _{\pi \in V(t,u)}\prod _{i=1}^t{y}_{\pi }^i$$
其中$V(t,u)$代表所有经过映射之后为序列$l$,且长度为t的序列，且在第t步输出为$l_u^{{}^{\prime }}$的集合。
由于$\alpha (t,u)$的每一个后续状态一定依赖于前面的状态转移，可以借助动态规划算法来求解,且$\alpha$数组大小为$(2U+1)\cdot T$。
初始状态:
所有label的正确开头都只能是空或者第一个字符，所以初始化状态为:
$$\begin{gathered} \alpha (1,1)=y_0^1 \\ \alpha (1,2)=y_1^1 \\ \alpha (1,u)=0u>2 \end{gathered}$$
递推关系:
$$\alpha (t,u)=y_{l_u^{{}^{\prime }}}^t\sum _{i=f(u)}^u\alpha (t-1,i)$$
其中:
$$f(u)=\{\begin{array}{rc}u-1& if{l}_u^{{}^{\prime }}=blank|l_{u-2}^{{}^{\prime }}=l_u{}^{\prime }\\ u-2& otherwise\end{array}$$
即当当前字符是空或者与上一个字符相同(这样必须通过一个空字符过渡),当前状态依赖于上一个时间步的前两个状态，即通过当前字符和当前的前一个字符转移而来。
而其他情况都可以由三个状态转移而来。
最终$p(l|x)={\sum }_{i=0}^{2u+1}\alpha (T-1,i)$

具体实现

特别注意，这里既有正常的概率计算，即alphas数组所对应的计算，也有log域下的计算，即tmp_alphas所对应的计算。
依旧参考了百度的逻辑。

double cpu_ctc::compute_alphas()
{
	int start = ((L + repeats - T) < 0) ? 0 : 1,
		end = (L * 2 + 1) > 1 ? 2 : 1;
	
	for (int i = start; i < end; i++)  
	{
		alphas[0][i] = Prob[0][labels_with_blanks[i]];      
		tmp_alphas[0][i] = log(Prob[0][labels_with_blanks[i]]);                        
	}

	for (int t = 1; t < T; ++t)   //模拟时间  
	{
		int remain = L + repeats - (T - t);   //label数量和剩余时间的差; 
		if (remain >= 0)                     
			start += this->s_inc[remain];      

		if (t <= L + repeats)	//若t>L,则已到末尾;       
			end += e_inc[t - 1];          
		                         
		int startloop = start;            
		if (start == 0)                              
		{
			alphas[t][0] = alphas[t-1][0]*Prob[t][blank_label];  //概率相乘; 
			tmp_alphas[t][0] = tmp_alphas[t - 1][0]+log(Prob[t][blank_label]);

			startloop += 1;
		}

		for (int i = startloop; i < end; i++)
		{
			double prev_sum = alphas[t - 1][i] + alphas[t - 1][i - 1];
			double tmp_sum = log_sum(tmp_alphas[t - 1][i], tmp_alphas[t - 1][i - 1]);

			if (labels_with_blanks[i] != blank_label && i != 1 && labels_with_blanks[i - 2] != labels_with_blanks[i])
			{
				prev_sum += alphas[t - 1][i - 2];
				tmp_sum = log_sum(tmp_sum,tmp_alphas[t-1][i-2]);
			}
			alphas[t][i] = prev_sum*Prob[t][labels_with_blanks[i]];
			tmp_alphas[t][i] = tmp_sum + log(Prob[t][labels_with_blanks[i]]);
		}
	}
	double score = 0;
	for (int i = start; i < end; i++)
	{
		score += alphas[T-1][i];
	}
	double tmp_score = -INFINITY;
	for (int i = start; i < end; i++)
	{
		tmp_score = log_sum(tmp_score, tmp_alphas[T - 1][i]);
	}
	prob_ = log(score);
	return score;
}

反向计算

原理

反向计算的主要任务是计算存储$\beta (t,u)$数组，并将最终求得的概率与前向求得的概率相验证。
其中:$$\beta (t,u)=\sum _{\pi \in W(t,u)}\prod _{i=1}^{T-t-1}{y}_{\pi }^{i+t}$$
$W(t,u)$的含义是所有在t时刻输出为$l_u^{{}^{\prime }}$的合法路径的从t开始至结尾的部分。
而$\beta (t,u)$就是从t+1时刻起一直到序列结尾的所有合法路径（映射后为$l$）的概率之和(不包含时间t)。
这样$\alpha (t,u)\cdot \beta (t,u)$的值就是，在t时间步输出$l_u^{{}^{\prime }}$的所有合法路径的概率之和。
反向数组的求解可以理解为前向的逆过程，其总体思路与前向相似。
初始化:
$$\begin{gathered} \beta (T-1,U^{\prime })=\beta (T-1,U^{\prime }-1)=1 \\ \beta (T-1,u)=0u<U^{\prime }-1 \end{gathered}$$
其递推公式为:$$\beta (t,u)=y_{l_u^{{}^{\prime }}}^t\sum _{i=u}^{g(u)}\beta (t+1,i)$$
其中:
$$g(u)=\{\begin{array}{rc}u+1& if{l}_u^{{}^{\prime }}=blank|l_{u+2}^{{}^{\prime }}=l_u{}^{\prime }\\ u+2& otherwise\end{array}$$

具体实现

反向计算的实现与前向计算极为相似，相当于前向计算的逆过程。
实现了正常概率与log域下的计算：

double cpu_ctc::compute_betas()
{
	int start = S > 1 ? (S - 2) : 0,
		end = (T > (S / 2) + repeats) ? S : S - 1;

	for (int i = start; i < end; i++)
	{
		betas[T - 1][i] = Prob[T - 1][labels_with_blanks[i]];
		tmp_betas[T - 1][i] = log(Prob[T - 1][labels_with_blanks[i]]);
	}

	for (int t = T - 2; t >= 0; t--)
	{
		int remain = (S / 2) + repeats - (T - t);
		if (remain >= -1)
			start -= s_inc[remain + 1];
		if (t < (S / 2) + repeats)
			end -= e_inc[t];
		int endloop = end;
		for (int i = start; i < endloop; ++i) {
			double next_sum=0;
			double tmp_sum = -INFINITY;
			if (i == S - 1)
			{
				next_sum = betas[t + 1][i];
				tmp_sum = log(betas[t + 1][i]);
			}
				
			else
			{
				next_sum = betas[t + 1][i] + betas[t + 1][i + 1];
				tmp_sum = log_sum(tmp_betas[t + 1][i], tmp_betas[t + 1][i + 1]);
			}

			if (labels_with_blanks[i] != blank_label && i != (S - 2) && labels_with_blanks[i] != labels_with_blanks[i + 2]) {
				next_sum += betas[t + 1][i + 2];
				tmp_sum = log_sum(tmp_sum, tmp_betas[t + 1][i + 2]);
			}
	
			betas[t][i] = next_sum * Prob[t][labels_with_blanks[i]];
			tmp_betas[t][i] = tmp_sum + log(Prob[t][labels_with_blanks[i]]);
		}
	}

	double score = 0;
	double tmp_score = -INFINITY;
	for (int i = start; i < end; i++)
	{
		score += betas[0][i];
		tmp_score = log_sum(tmp_score, tmp_betas[0][i]);
	}

	for (int i = 0; i < tmp_betas.size(); i++)
		for (int j = 0; j < tmp_betas[0].size(); j++)
			if (tmp_betas[i][j] != -INFINITY)
				tmp_betas[i][j] = (tmp_betas[i][j]- log(Prob[i][labels_with_blanks[j]]));
	for (int i = 0; i < betas.size(); i++)
		for (int j = 0; j < betas[i].size(); j++)
			betas[i][j] = betas[i][j] / Prob[i][labels_with_blanks[j]];
	
	return score;
}

需要特别注意，在计算出beta数组后，我们将其除以Prob[i][labels_with_blanks[j]]，这是因为alpha数组和beta数组相乘时多计算了一次$y_k^t$自身的概率。

损失函数

CTC的损失函数定义如下:
$$L(S)=-ln\prod _{(x,z)\in S}p(z|x)=-\sum _{(x,z)\in S}lnp(z|x)$$
其中S为训练集。损失函数可以解释为:给定标签序列和输入，最终输出正确序列的概率，为了方便计算，我们将这个概率取负对数。而我们将取负对数之后的loss最小化就是将输出概率最大化。
而我们在前文提到过，对于给定的$t$和$u$,%=$\alpha (t,u)\beta (t,u)$有具体的意义。具体写成公式为:
$$\alpha (t,u)\beta (t,u)=\sum _{\pi \in X(t,u)}\prod _{t=1}^T{y}_{{\pi }_t}^t$$
$X(t,u)$表示在t时刻经过$l_u^{{}^{\prime }}$所有路径的集合。
进一步可以转换为
$$\alpha (t,u)\beta (t,u)=\sum _{\pi \in X(t,u)}p(\pi |x)$$
所以给定任意时刻的t，我们的概率可以这么计算:
$$p(z|x)=\sum _{u=0}^{|z^{{}^{\prime }}-1|}\alpha (t,u)\beta (t,u)$$
其意义是任意时刻，从扩展后的序列的第一个位置直到最后一个位置的$\alpha \beta$乘积之和就是$p(z|x)$的概率，只要理解了上述所说的所有公式的含义，这个式子也是可以直观理解的。在任意t时刻，如果合法的路径有经过$l_u^{{}^{\prime }}$的,乘积一定不为0，而不经过一定为0.所有乘积加和一定是最终概率。这也很容易通过程序验证。

梯度计算

公式推导

相信大家在看前向过程和反向过程的时候都是一知半解，搞懂了这些计算的过程，但并不明白计算这些中间结果到底有什么用。其实，前面所做的都是为最终求输入的梯度做准备的。
假设经softmax后的网络输出$y_k^t$代表第t个时间片第k个字符的输出概率。则损失函数对网络输出$y_k^t$的偏导数为:
$$\frac{\partial L(x,y)}{\partial y_k^t}=-\frac{\partial \ln p(x|z)}{\partial y_k^t}=-\frac{1}{p(x|z)}\frac{\partial p(x|z)}{\partial y_k^t}$$
关于p(z|x)的意义已经在前文中解释过。
在这里，我们再引入一个集合 B ( z , k ) = { u : z u ′ = k } B(z,k)=\{u: z_{u}^{'}=k\} B(z,k)={u:zu′​=k},它的意义是，k出现在序列z’中的所有路径的集合(注:k可以在路径中任意位置出现)。这样，我们很容易得到:
$$\frac{\partial \alpha (t,u)\beta (t,u)}{\partial y_k^t}=\{\begin{array}{rc}\frac{\alpha (t,u)\beta (t,u)}{y_k^t}& ifkoccursin{z}^{\prime }\\ 0& otherwise\end{array}$$
因此，损失函数关于输出的偏导数可以写为:
$$\frac{\partial L(x,y)}{\partial y_k^t}=-\frac{1}{p(x|z)}\frac{\partial p(x|z)}{\partial y_k^t}=-\frac{1}{p(x|z)y_k^t}\sum _{u\in B(z,k)}\alpha (u,t)\beta (u,t)$$
最后，由于该输出是经过softmax之后的结果，我们需要得到损失函数对于softmax之前的网络输出的${\alpha }_k^t$的偏导数。
$$\frac{\partial L(x,y)}{\partial a_k^t}=-\sum _{k^{\prime }}\frac{\partial L(x,y)}{\partial y_{k^{\prime }}^t}\frac{\partial y_{k^{\prime }}^t}{a_k^t}$$
其中，$$y_k^t=\frac{e^{a_k^t}}{{\sum }_{k^{\prime }}{e}^{a_{k^{\prime }}^t}}$$
针对未经softmax的网络层的梯度推导较为复杂，具体推导如下:
$$\frac{\partial y_{k^{\prime }}^t}{\partial a_k^t}=\frac{\partial \frac{e^{a_k^t}}{{\sum }_{k^{\prime }}{e}^{a_{k^{\prime }}^t}}}{\partial a_k^t}$$
这里需要特别注意，当$y_{k^{\prime }}^t$中的k’与$a_k^t$中的k相同与不相同时求导结果是不同的,当$k^{\prime }=k$时:
$$\frac{\partial \frac{e^{a_k^t}}{{\sum }_{k^{\prime }}{e}^{a_{k^{\prime }}^t}}}{\partial a_k^t}=\frac{e^{a_k^t}({\sum }_{k^{\prime }}{e}^{a_k^t}-e^{a_k^t})}{({\sum }_{k^{\prime }}{e}^{a_k^t}{)}^2}$$
注意:这里的$\frac{e^{a_k^t}}{{\sum }_{k^{\prime }}{e}^{a_{k^{\prime }}^t}}$可以提取为$y_k^t$,所以原式可化简为:
$$\begin{gathered} =y_k^t\frac{{\sum }_{k^{\prime }}{e}^{a_k^t}-e^{a_k^t}}{{\sum }_{k^{\prime }}{e}^{a_{k^{\prime }}^t}} \\ =y_k^t(1-y_k^t) \end{gathered}$$
同理，当$k^{\prime }!=k$时:
$$\begin{gathered} \frac{\partial y_{k^{\prime }}^t}{\partial a_k^t}=\frac{-e_{k^{\prime }}^t{e}_k^t}{({\sum }_{k^{\prime }}{e}^{a_k^t}{)}^2} \\ =-y_k^t{y}_k^t \end{gathered}$$
当$k^{\prime }=k$时(所有)，有:
$$\frac{\partial L(x,y)}{\partial a_k^t}=\frac{y_k^t-1}{p(z|x)}\sum _{u\in B(z,k)}\alpha (u,t)\beta (u,t)$$
当$k^{\prime }!=k$时(k’!=k的所有情况)，有:
$$\frac{\partial L(x,y)}{\partial a_k^t}=\frac{y_k^t}{p(z|x)}\sum _{k^{\prime }}\sum _{u\in B(z,k^{\prime })}\alpha (u,t)\beta (u,t)$$
相加得:
$$\frac{\partial L(x,y)}{\partial a_k^t}=\frac{y_k^t}{p(z|x)}\sum _k\sum _{u\in B(z,k)}\alpha (u,t)\beta (u,t)-\frac{1}{p(z|x)}\sum _{u\in B(z,k)}\alpha (u,t)\beta (u,t)$$
特别注意:
${\sum }_k{\sum }_{u\in B(z,k)}\alpha (u,t)\beta (u,t)$就是$p(z|x)$(可以从公式的实际意义出发，该式意义就是所有可以映射为z序列的序列之和)
于是,最终的式子可以化简为:
$$y_k^t-\frac{1}{p(z|x)}\sum _{u\in B(z,k)}\alpha (u,t)\beta (u,t)$$
这就是我们最终推导出的相对于未经softmax的网络层的梯度。

具体实现

梯度的计算主要难在公式的推到上，将最终的公式推导出后具体实现较为简单，就是对公式的简单模拟。

void cpu_ctc::compute_gradients()   //将梯度改为log域!!!
{
 	for (int t = 0; t < T; t++)
	{
		for (int k = 0; k < Prob[0].size(); k++)  //a{k,t} 表示在t时刻k字符的梯度; 
		{
			double prob_sum = 0;
			double tmp_sum = -INFINITY;
			
			for (int i = 0; i < 2 * L + 1; i++)
			{
				if (labels_with_blanks[i] == k)
				{
					prob_sum += alphas[t][i] * betas[t][i];
					tmp_sum = log_sum(tmp_sum, log_prod(tmp_alphas[t][i],tmp_betas[t][i]));
				}
			}

			gradients[t][k] = Prob[t][k] - prob_sum / exp(prob_);
			if(tmp_sum!=-INFINITY)
				tmp_gradients[t][k] = Prob[t][k] -  exp(tmp_sum - prob_);
			else
				tmp_gradients[t][k] = Prob[t][k];
		}
	}
}