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}{r}u-1\\ u-\end{array}$$