1. 时间序列问题
首先，我们还是从问题开始入手，"时间序列"就是一系列事件(events)随着时间的变化而改变，比方说一个人来了个后空翻，而旁边放了个相机将这个动作按每秒拍10张照片的频率给记录了下来，那么如果动作一共是10s，那就会有100张图片，这每一张图片都是在不同的时间下的发生的事件，而这些事件组成了一个时间序列。同理，如果你说了一段话，一个录音笔按照某个频率(如1s记录10个音)把你说的话给记录了下来，那么记录的每一个音也组成了一个时间序列。所以你会发现生活中很多东西都是时间序列：一天内的室内温度、一个月以来的某支股票价格等等。
那么对于这些时间序列，我们往往会想了解两个时间序列之间的相似度，比方说你家安装了个语音识别开门系统，那么你会先记录你“我要开门”的声音，系统记录下来了之后，你以后每次说“我要开门”，他都会对这两个时间序列进行匹配，如果相似度高就给你开门，相似度不高就不给你开门。而DTW就是计算时间序列相似度的一个非常经典和常用的方法。
2. 相似度计算
2.1 Euclidean Distance (ED)
在介绍DTW的相似度计算之前，我们可以先说说一个更基础也更简单的时间序列相似度计算方法: Euclidean Distance。通过下图，我们可以看到一个时间序列Q和一个时间序列C，这两个时间序列的长度都是N。那么ED的计算方法就是在每个时间段将Q的值减去C的值，再把这些值加起来。比方说在第i个时间点，用Qi - Ci，我们将会得到每个时间段对应的Di，然后我们把从i=1到i=N的Di都算出来然后相加，得到的就是这两个时间序列之间的总距离D，总距离越大，这两个时间序列的相似度就越小。当然，我们往往不是单纯相减，因为会有正负的问题会相互抵消，所以我们会用平方距离或者用绝对值，公式如下：
然而ED存在一个问题，就是你会发现其实Q和C长得很像，形状基本上一样，只不过有一个time shift，就可能是同一个人说的“我要开门”，但是一次比另一次完全慢两拍。而在这种情况下，你会发现ED计算出来的相似度并不大，这是为什么呢？原因就是出现在时间序列的匹配上，换句话说Qi减错Ci了，可能是应该减下一个时间段的C，也就是C(i+1)，而DTW就是为了解决这种匹配的问题。
2.2 Dynamic Time Warping (DTW)
接下来我们来介绍DTW，在介绍DTW之前我们需要说明一下，其实我们在对比两个时间序列的时候，这两个时间序列的长度不一定是相同的，比方说你喝醉酒之后说“我要开门”就很有可能比清醒的时候说语速要慢很多，所以一个时间序列的长度是n，一个时间序列的长度可能是m。因此，DTW里我们也不需要将Q和C完全一一对应，DTW求的是最匹配的Q和C之间的距离，而定义最匹配也就到这个时间段为止哪个C和Q对应的总距离最短，那么就是这个Q和这个C对应。也就是说某个时间段的Qi可能完全没有与之匹配的Ci(如"喝醉的时候发出奇怪的噪音，不需要匹配“），而某个时间段的Qi可能匹配3个Ci（如“喝醉的时候发‘我’的音拖了3个时间段“）。以上只是让你大概了解一下背后的intuition，我们可以正式的介绍了，请看下图。
我们可以看到，时间序列Q和时间序列C之间不再是一一对应的了，也不是像ED算法一样单纯的将一个时间段上的Q(i)-C(i)， 而是会寻找到最匹配的点，如Q(i)-C(i+1)，比方说现在Q时间序列的最高点减去C时间序列的最高点，Q时间序列的最低点减去C时间序列的最低点，那么你可想而知最后的distance肯定相对ED算法较小，最后的相似度也会更大。因此通过这个办法，你应该理解我们可以有效解决之前ED算法出现的问题。那么这个算法是如何实现的呢？为了高效的实现，我们可以通过一个机器学习中非常经典且常用的算法：Dynamic Programming(动态规划)。
2.2.1 Dynamic Programming Implementation
其实我们到现在应该都了解求两个时间序列之间的相似度，就是需要求这两个时间序列之间的距离之和，而DTW求的距离是根据最匹配的点来求，而最匹配的定义就是让总距离的值最小。因此动态规划需要解决某个时间点上的Q和哪个时间点上的C匹配，可以让总距离最小。
首先，我们如上图右建立一个Matrix(我们称该Matrix为"DA")，y轴是时间序列Q，长度是N，x轴是时间序列C，长度是M，因此是个NxM的matrix。这个matrix上的左下角就是Q的第一个时间段的值，Q(1)，和C的第一个时间段的值，C(1)，的绝对差。而右上角的点是Q的最后一个时间段的值，Q(N)，和C的最后一个时间段的值，C(M)，的绝对差。DTW就是需要求从左下角的起点到右上角的终点，如何匹配才能让总距离最短(注意：Q的起点需要匹配C的起点，Q的终点需要匹配C的终点，其余的点由匹配程度确定)，也就是找到一条从起点到终点的使总距离最短的‘路径‘。
如果你觉得有点难理解的话，没有关系，我们现在先一起把Matrix填满：
【注：该Matrix里的Y轴-“R”时间序列就是前一个图中的“Q”时间序列，X轴-“V”时间序列就是“C”时间序列】
从上图我们可以看出，Matrix左边第一列填满的方式为：|R(i) - V(j)| + DA[i-1, 0]，这是什么意思呢？举个例子，如果我们看标记出来的红色框内的一个值20，这个值代表的是，如果“R”时间序列(Y轴)的第5个时间段的值R(4) 和“V“时间序列(X轴)的第一个时间段的值V(0)匹配的话【注：index从0开始】，那么截止到目前的最短总距离该会是多少？也就是说不仅仅要求这两个点之间的距离（8-1=7），还要加上之前的匹配所计算出来的最短总距离，而我们很显然之前匹配出来的最短总距离DA[3, 0]是13(也就是下方的格子)，那就是7+13 =20。
同理，如果你看红色框选出来的另一个值11，这个值对应的是“R”时间序列(Y轴)的第4个时间段的值R(3)和“V“时间序列(X轴)的第4个时间段的值V(3)匹配的话，这事这两个点之间的距离为（9-3=6），而之前的匹配出来的最短总距离是多少呢，这时我们看与这个格子相近的左下方的3个格子，看看哪个数最小，我们就选哪个，很明显最小值是5，所以这个格子的值为（6+5=11）。
你可能会纳闷为什么是左下方这3个格子呢？其实你要回想，这个matrix的每个格子都代表的是：如果匹配这个格子，到目前为止的最短距离是多少，所以你只要不断的把目前的距离加上之前的最短距离，直到最后，你就会找到一个最短总距离的解了。而另外一个红框算出来是9，具体细节我就不解释了，其实计算方法和之前两个红框都是异曲同工。
【注：为了简化，我们的序列是不会有配对交叉的情况出现的，如R(2)配对C(3), R(3)配对C(2)，这样就出现了交叉，你也可以想象在现实生活中确实也不太常见(语序颠倒？)，这也是为什么我们的表格计算只需要考虑这个格子的左下方的三个格子，而不是他附近的全部格子， 比方说如果你考虑右上方的格子就必然会产生交叉了】
最后，当你把整个matrix填满之后，你只需要从右上角的终点开始往起点追溯，比方说终点的值是15，他附近的3个格子最小值是15，因此前一个匹配就是15的这两个点，也就是“R”时间序列的第9个点和“C”时间序列的第9个点。同理，你可以一直追溯回起点，这样你就找到了，这两个时间序列的所有应该匹配的点，以及最短的总距离15。
2.2.2 Code Snapshot
如果感兴趣的话，以上就是DTW基于Dynamic Programming的代码实现，如果你理解了上述填满Matrix的过程，那么这个代码也一定非常好理解。其中，与上面解释的过程稍稍有一点不一样的地方在于代码第19行，也就是你会看到设置了一个boundary。比方说对于R(1)， 你不再需要去将这个点和全部C的点尝试配对和计算最短总距离，因为我们压根不希望一个R点配对很多很多的C点，因为如果你套入现实情况考虑会发现不会有什么意义，因此对于R(1)你只需要计算与前后10个时间段的C的配对情况就好。如果在图上表示这个Boundary，请参照下图的绿色线，而R值的意思就是你希望一个Q最多能连接多少个时间点的C。
而“euc_dist"只是一个很简单的function，就是求两个值相减后的绝对值。最后一点，就是代码第24到26行，在之前我们只考虑格子左下方的三个格子，而在这里我们会考虑更多的格子，因为会有跳过的情况(有些时间点可能根本不需要匹配)。
DTW的原理以及算法实现基本上就是这样，但其实DTW也存在一些缺点，比如时间复杂度太高: O(NM)，我们的前辈也尝试了非常多提高DTW效率的优化方法，接下来我们一起来探讨。
3. DTW 的效率优化方法
从上述的过程中你可以直观感受到，基本上每个时间序列上的每个点都要做计算，所以时间复杂度为O(NM)，然而对于现实生活中的问题，这两个序列可能是非常大的，如万亿级别，因此我们很明显需要做大量的效率优化。不过为了更好地理解，我们先考虑一个情形，比方说你现在有一个很长很长的时间序列“T"，然后你有一个较短的时间序列“Q”，你需要在这个很长的时间序列“T”里看看能不能找到和这个“Q”非常相近的时间序列。于是乎，你把“T”分解成一个个小一点的时间序列“C”(如分成了m个“C”)，然后你会使用我们以上学到的办法来尝试计算每个“C”和“Q”的相似度，从而找到最接近的“C”有多接近。好，接下来就让我们开始吧。
3.1 Early Abandoning Z-Normalisation
首先，什么是Early Abandoning？我们先通过下图来了解一下。
Early Abandoning就是比方说当你在计算距离时，你突然发现目前这两个时间序列还没算出最终的距离，就已经发现比之前某个“C”和“Q”的距离要长了，也就是说对于目前的这个“C”已经没有计算下去的必要了，因为已经有一个“C”(best-so-far)一定比目前这个“C”更接近“Q”了，所以你可以直接提前abandon(放弃)计算这个“C”和“Q”的相似度。
好的，那接下来我们说Normalisation，首先我们也需要说明为什么要normalisation呢？其实normalisation是非常重要的一步，而且其实我们每次对比一个"C"和"Q"的时候，都需要确保这两个时间序列已经normalised过了(当然Q只需要normalised一次，因为都是同一个Q)。为什么有这个必要呢，我们看下图。
上图的中间就是一段录像("T")中的截图(“C”)，而左边和右边的照片是分别左右scale了110%和90%的情况，尽管这个改变看上去并不明显，但是他们可以将检测错误率增加一倍以上。因此，你想现实生活中，可能用户离镜头近一点远一点，或者高度差别，都会带来很大的影响，由此可见normalisation的重要性。
因此，Early Abandoning Z-Normalisation背后的想法就是，如果你可以使用online Z-normalisation(即计算两个序列的距离同时，同步normalised时间序列)，那么每当你Early Abandoning时，你不仅仅可以节省掉继续计算该序列的距离的时间复杂度，你还可以减掉normalisation的时间复杂度。实现的办法则是通过针对每次要对比的“C”，如果我们把“C”设置成为长度m，也就是每个“C”的长度都是m，那么我们可以使用下方公式先求出这个“C”时间序列的平均值和标准差。
然后先不把整个C序列里的每个时间点的值都normalised，只在每次匹配这个时间点的距离的时候顺便normalised，这样的话，只要early abandoning，就会节省normalised的复杂度。以下是《Searching and Mining Trillions of Time Series Subsequences under Dynamic Time Warping》的作者给出的online normlisation的代码思想，尤其注意第11行就是online normalisation的实现，感兴趣的朋友也可以去Reference查看完整源代码。
3.2 Reordering Early Abandoning
了解了Early Abandoning后，其实Reordering也比较好理解，我们一起来看下图。
你可以看到上图的左边是传统的计算距离的顺序，从左到右，可是这不一定是最优秀的顺序。比方说我们现在已经有一个best-so-far的距离为50，那么左边这个顺序可能要计算9次距离才发现大于best-so-far，于是乎才停止。可是右边这个图，通过使用更优秀的顺序，我们计算5次就发现已经大于best-so-far，所以节省了时间。
那么这个最优秀顺序要如何确定呢？其实你想求的是要看看Q的每个时间段对应C的哪个时间段的距离最大，就从哪个开始，可是这样好像运算量也有点大。而一个比较好的办法就是，我们知道通过normalised，许多C序列都会是高斯分布，因此，我们如果从离C序列的平均值最远的点开始，往往就是对距离的增长贡献最大的点。
3.3 Reversing the Query/Data Role
在介绍这个方法之前，需要简单介绍一下什么是Lower Bounding。
上图左边是一个“Q”和一个“C”时间序列，我们现在知道我们要计算“Q”和“C”之间的距离，但是不是每个时间点都需要去做计算的，而去掉这些点我们可以大幅度提高我们的效率。有人就提出了Lower Bounding的概念，你可以看到右图上，Q(红线)上下出现了一个U线和一个L线，分别表示upper bound和lowe bound，现在你不再是需要用Q来和C计算距离了，你只需要用L和U来计算，具体用U还是用L，就是看U和L哪个离C比较近，所以你可以看到右图左边是用L来计算距离，右边则是用U来计算距离。
U和L的制定有非常多的方法，包括：LBKeogh、LB_kim。。。，在这里我们就不一一介绍了，我在下面给出了LB_Keogh的公式，相信结合上面的图片都能很容易理解。而这些不同的LB方法的区别主要在upper bound和lower bound离原来的Q线有多接近，过于接近的话就没有起到太多效率优化的效果，离原来的Q线太远的话，就会少算太多的点，需要权衡。
说回Reversing the Query/Data Role，其实通常来说不需要reversing，一般采取的办法就是在Q线上设置U线和L线，因为这样Q只用算一次，所以U和L也只需要算一次。但有的时候，如果所有的Lower Bounds都没有办法起到剪枝效果，那么我们也可以reverse，就是在C线上设置一个U线和L线，如下图所示，具体在什么时候需要Reverse，你继续往下看就知道了。
3.4 Cascading Lower Bounds
之前提到了Lower Bounds，确实这是一个提高时间序列计算效率的最有效的方法之一了，因此也有不少不同的Lower Bounds，也很难定义哪个比别的就优秀，像我们前面说的，这是一个权衡的问题。我们的前辈就试过把一些主要的Lower Bounds在50个不同的数据集上做了测试，并制作了以下的图，Y轴是U线与L线离Q线的紧凑度(Tightness)，X轴是时间复杂度，而作者指出在一般情况下，应该选择虚线上的Lower Bounds方法，因为根据这个图，如果不在这个虚线上，意味着在同一个紧凑程度下，存在更有效率的Lower Bounds。
而Cascading Lower Bounds的意思是，我们可以在计算DTW的时候，通过Cascading的方法，把所有的LB方法都用上，也就是说我们可以先尝试时间复杂度为O(1)的Lower Bounds，如LB_KimFL, 可如果C没有被剪枝成功的话(也就是并没有超过best-so-far)，我们可以再使用LB_KeoghEQ，也就是时间复杂度为O(N)了，如果还是没有用，那就可以使用我们之前介绍的Reversing方法，也就是LB_KeoghEC，以此类推直至没有任何Lower bounds的情况。通过以上办法，我们可以剪枝掉99.9999%的DTW时间序列距离计算。
4. 总结
以上就是我对DTW时间序列相似度计算以及效率优化的方法汇总，大部分的内容(尤其效率优化方面)来自《Searching and Mining Trillions of Time Series Subsequences under Dynamic Time Warping》，原作者实现的源代码我也放在了Reference上，有兴趣的朋友可以去拿来啃一啃。总的来说，经过将DTW通过一些列的剪枝优化，这已经成为了一个计算时间序列相似度简单又高效的办法，因此也可以应用在各类的时间序列相似度计算上。
Reference
- 《Searching and Mining Trillions of Time Series Subsequences under Dynamic Time Warping》
- Sakoe, H., Chiba, S. (1978). Dynamic programming algorithm optimization for spoken word recognition. IEEE Trans. Acoustic Speech and Signal Processing, v26, pp. 43-49.
- Souza, C.F.S., Pantoja, C.E.P, Souza, F.C.M. Verificação de assinaturas offline utilizando Dynamic Time Warping. Proceedings of IX Brazilian Congress on Neural Networks, v1, pp. 25-28. 2009.
- Mueen, A., Keogh. E. Extracting Optimal Performance from Dynamic Time Warping
- 源代码：http://www.cs.ucr.edu/~eamonn/UCRsuite.html
794
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
