Dynamic Time Warping - UCR Suit 优化详解

前言

本文是对 Searching and Mining Trillions of Time Series Subsequences under Dynamic Time Warping[1]读书笔记的总结,内容上比较偏重Paper中提到的几个优化算法。因为我个人在初读此文章时,对文中提到的一些优化的基础知识不是很了解,所以想在此详细展开讨论一下,以便日后可以给阅读原文存在困难的朋友提供一些基础理论层面上的剖析。

在叙述中我选择按顺序先从 DTW 基础概念入手,以便读者可以更直观地理解DTW是怎么工作的,为什么需要优化它?然后剖析在当时比较流行的几种优化方法,从而引出核心部分:四种新的优化 DTW 的算法。由于该文已经被无数案例证实过其实用价值,故在此不会做更多的试验分析结果方向的赘述,而是侧重分析优化算法的理论部分。

为了方便以后能够重现该论文需要了解:原始论文中提到了三个假设,且该论文的优化理论也建立在这三个假设上。具体内容如下:

(1) 所有时间子序列必须被归一化处理

(2) DTW 是已知最好的时间序列搜索方法

(3) 对任意长度的时间子序列的查询没有适当的方法建立索引


1 DTW 核心思想

1.1 关于时序数据的一些基础知识

DTW 是解决对比时序数据中子序列搜索的方法。所谓时序数据[2]是指时间序列数据。时间序列数据是同一统一指标按时间顺序记录的数据列。

Building on the work done for rCharts profiled in the Timely Portfolio piece, the dygraphs R package provides an interface to the dygraphs javascript library. With just a few lines of R code, it is now possible to produce charts that approach the polished look of the professional

如图所示时序数据中的数据点分别代表 IBM 和 Linkedin 公司2010 - 2016 的股票信息。为了更方便读者理解数据这里借助图表更直观地展示出来。通过观察不难看出蓝色和红色线在图表上呈现出不同的形态,换句话说这两个曲线相似度不太高(因为他们表现出来的形状很不一样)。

在现实中我们经常会遇到对比两段时序数据相似度的问题,比如图中的股票曲线,通过计算两个股票走势是否相似,用于辅助投资决策。 计算相似度一个最直观的方法就是用欧式距离 Euclidean Distance (ED)

图中 Q 代表查询 (Query 的缩写),C (候选对象 Candidate)代表被查询的子序列。i 时序C的起始点, kC 的结束点。如图中定义所示,在计算距离时需要对每个 Q 在 i 点上的值 减去 C 在 i 点(与q在同一个时间点)的值,再对差求平方,最后将 n 个时间点的平方和开根号求出距离。这个距离就是所谓的两个时间序列相似度的度量,即欧式距离。

问题是:欧式距离在计算中 Q 和 C 是一一对应的,也就是说 Q 在 i 上的点只能对应 C 在 i 时间上的点计算。这样计算的缺点是忽略了两个时间序列数据上的位移(offset)。比如上图中的曲线,两个波峰并没有重叠在一起,这导致了尽管他们形状(相似度)很接近,但由于这个 offset 的存在,欧式距离的结果有可能会偏大。另一个更经典的例子是 sin 和 cos 函数。其实他们只差了 \pi/2 的距离,如果把其中一个向左或向右移动 \pi/2 那么他们就会完全重合在一起。如果这两条线重合在一起,就代表他们之间的欧氏距离为0,即:形状完全一样的两条曲线。正如上文所讲,因为欧氏距离没有考虑到数据上的偏移量对计算带来的误差,所以在现实中计算距离的结果并不理想,而这个缺点却很好的被 DTW 算法所弥补。

1.2 DTW 简介
由于 DTW 论文[3]出现得很早, 具体细节不是本文重点,所以这里只做关键性的介绍,以便读者在短时间内能对该算法有个初步的了解,方便阅读后面内容。

DTW 中的 W 是 Warping 的缩写,意为“扭曲”。如图

Q 和 C 上的点并不是一一对应的,而是一对多的对应关系。在计算距离时可以用一个 n x n 的矩阵来辅助计算(如 Q 与 C 不等长则为 m x n 的矩阵)。n 为 Q 和 C 的长度,在计算中用 i 和 j 表示矩阵中的元素下标。计算的目标是从矩阵的左下角(0,0)开始到右上角(n,n)结束的众多路径中找到一条 “Warping Path”,且规定该路径必须为对角上两点间的最短路径。矩阵中的值代表连接两点间的开销。

图中红色曲线就是上文谈到的 Warping Path,即在给定开销(红色线相邻两个值的权重)情况下,自左下到右上角的最短距离。在求解中可以把它考虑成一个动态规划在有向无环图中求最短路径的问题。注意:在这个例子中每次只能向上和向右移动,实际中可以调大移动步长。

论文中还提到了 Sakoe-Chiba Band ,它硬性限制了“扭曲” 不能超过 R 个长度。因为本文重点介绍基础知识,在此不对更多赘述。
注意:R (步长最大值)会直接影响 DTW 的计算复杂度,R 越大 DTW 复杂度越高。

2 时间序列已知优化算法(2012年以前)

论文中还提到了利用多核并行化优化。因为在这里不涉及到具体数学优化算法,所以不做介绍。

2.1 去平方根操作

在 DWT 和 ED 计算中都有带有开平方根的操作,为了方便理解后边优化在这里就去掉了平方根,因为这两种累加计算是单调递增的,所以去掉平方根对排序相似度结果上没有任何影响。

2.2 使用下界

通过使用下界可以提前筛选掉部分不符合条件的候选结果。我们可以这样理解下界这个概念:如果想找到班里最高的同学,常规做法是准确测量全班同学的身高并按身高排序,再择测得最高身高数据的那个。这样做当然可行,但想得到最终结果我们必须测量每一个同学的身高,这个工作量是相对比较大的。但也可以换个方法考虑:比如通过目测估计,选出一个中等偏上高的同学作为标杆,所有目测比他矮的同学肯定都不是要找的最高的那个,所以就不用花时间准确测量这些目测都不达标同学的身高。这样就通过这个中等偏上身高的同学(下界)这届排除了一大批是最高那个同学的可能。最后只需要测量小部分比这个下界高的同学,就能很快找到最终结果。

计算下界有很多种不同的方法,有些虽然复杂度小,但筛选效果未必很理想,有些效果好,但是计算量却很高。本文中提到了下面两种被工业广泛应用的计算下界方法。

2.2.1 LB_Kim_FL

该算法其实是 LB_Kim 算法的一个在归一化处理后的简化版本。原始的 LB_Kim 会一共取了4个观测点,即:初始点,终止点,时序数据最大值,时序数据最小值。我们在这里提到的 LB_Kim_FL 只用了起始和终止值,因为在归一化后的极值在这里几乎不起任何作用,但是在省略它们后整个算法的时间复杂度却从 O(n) 降为 O(1)。复杂度的下降是省略了遍历整个数组寻找全局最大和最小值的过程而做到的。经实践证明,这种方法在大多数情况下十分地粗暴有效,但在起始值和终止值差别较大时会对筛选结果造成一定的误差。

计算公式为:

LB_Kim_FL (Q, C) = max { | First (Q) - First (C) |, | Last (Q) - Last (C)| }

2.2.2 LB_Keogh

与简单粗暴的 LB_Kim 不同,LB_Keogh 计算时不只关注四个关键值,对查询子序列 Q 引入了 Upper Bound 和 Lower Bound 的概念。这种方法很好的解决了 LB_Kim 只关注起始值与终止值的问题,从而达到了更加精确地定义下界的目的,但由于它需要遍历整个样本才能完成计算,故复杂度为O(n)。

具体计算步骤如下:

(1) 定义一个时间序列子序列 Q : { q_{1} ... q_{n} } 以及窗口值 r 且在其范围内寻找 max 和 min 分别作为 Upper Bound,Lower Bound。

(2) 确定 Upper Bound : U_{i} = max (q_{ i-r},  q_{i + r})

(3) 确定 Lower Bound: L_{i} = min ( q_{i - r}, q_{i + r} )

(4) 利用公式计算 LB_Keogh:

LB_Keogh (Q, C) = \sqrt{\sum_{i = 1}^{n}{\left( c_{i} - U_{i} \right)^{2} I_{ \{c_{i} > U_{i} \} } + \sum_{i = 1}^{n}{  (c_{i} - L_{i})^{2}  I_{ \{c_{i} < L_{i} \}}    }}}

引用一位大牛对这个算法的理解: LB_Keogh 方法首先求出 Q 序列的上下包络线,然后对data序列 ( C ) 与上下包络线进行比较,如果不在上下包络线的范围内,就对该点与对应的包络线上的点求欧几里得距离(此处的欧几里得距离主要是指Y轴上的距离,并非二维),最终求和得到误差,与之前得到的误差进行比较,整条data中误差最小的就是目标序列。

2.3 Early Abandoning of ED and LB_Keogh

我们在之前讨论的 Lower Bound 剪枝其实就是一种 Early Abandoning。在计算中涉及到一个 best-so-far 的概念,在初始化时这个值被设置为无穷大。本着一旦找到更好的结果就替换的原则,然后用下界(Lower Bound)与这个值对比,如果下界大于它,那么就可以在此终止后面计算且筛选掉这条不符合条件的数据,反之则取代(更新)当前 best-so-far。

2.4 Early Abandoning of DTW

对于计算完 LB_Keogh 再计算 DTW 的复杂度也有相应的优化方法,在原文中作者也介绍了他们自己的方法。文章指出在 LB_Keogh 已经计算完一次后,这个结果可以被重复使用在计算 DTW 中,以省略重复计算带来的开销。

如上图所示,在计算 DTW 时从 K = 0 位置开始,计算至 K = 10 结束,这时我们得到了部分 DTW 距离累加的结果。在此可以用该部分结果与之前 LB_Keogh 的 K = 11 后的结果拼接,并把这个拼接的结果当做 Lower Bound 对比 best-so-far 决定是否要提前终止计算。

这种混合 Lower Bound 计算可以表示为:

LB = DTW(Q_{1: K}, C_{1: K}) + LB_Keogh( Q_{K+1, n}, C_{K+1, n} )

3 本文提出的四种创新优化算法

3.1 Early Abandoning Z-Normalization

作者在这里吐槽了一下之前人研究成果的疏忽,他指出尽管计算归一化的开销远大于计算欧氏距离的开销,但迄今为止(2012年)还没人提出过对这个归一化过程做优化,并在此引出他们的创意:在计算欧式距离或者 LB_Keogh 的过程中同时进行 online Z-Normalization,而不是分别做。而且一旦发现 C 子序列已经可以排除在最优结果之外了也没有必要继续对该序列做剩余归一化的操作。这个方案可以节省的对时间序列做完整归一化操作带来的巨大开销,因为这里的归一化操作只是针对部分子序列的计算而不是全部。

常规来讲想要对某个子序列进行 Z-Normalization 需要先计算其整体的均值和标准差。

即: \mu = \frac{1}{m} \sum_{}^{}{x_{i}}  , \sigma^{2} = \frac{1}{m} \sum_{}^{}{ x_{i}^{2} - \mu^{2} } (这里的标准差计算和公式略有不同,因为作者展开了平方差项)

然后用遍历样本后求得的值带入归一化公式: Normalized-C_{i} = \frac{x - \mu}{\sigma}

用这种常规方法计算均值和标准差是需要遍历整个子序列的,这个造成了一定的开销。在此原文作者提出了他改良过的方法,专门计算子序列归一化涉及到的均值和标准差的计算。

即: \mu = \frac{1}{m}  \left(  \sum_{i = 1}^{k}{x_{i}}  -  \sum_{i=1}^{k-m}{x_{i}}  \right), \sigma^{2} = \frac{1}{m}  \left(  \sum_{i = 1}^{k}{x_{i}^{2}}  -  \sum_{i=1}^{k-m}{x_{i}^{2}}  \right) - \mu^{2} 且在计算过程中储存中间输出的部分子序列加和结果加和平方结果,这样做可以避免重复计算。

这里的 m 指的是时间子序列的长度, k 为整个时间序列中从起点到截取的第 k 个时间段序列。

此处需要用图解释 m 和 k 在整个时序数据中的位置 TBD

论文中,作者提到此处存在浮点计算在上亿条数据下回出现累加误差的问题,并提出这里每比对100万子序列后,就要进行一次完整的Z-normalization。

3.2 Reordering Early Abandoning

前面 Early Abandoning 的计算方式都是从子序列的起点自左向右进行计算的。作者在这里提出一种先快速找到子序列 Q 和 C 之间差值之和最大的部分,然后以它为中心计算 LB 来判断这个子序列是否大于best-so-far值,从而可达到降低计算成本的目的。

引用下原文的图:

左图中计算的顺序是从左向右的,在计算到第9个时间点时发现 Lower Bound 已经超过了 best-so-far 值,然后终止了计算并筛选掉 C。 右图中将顶点放在离 Q 峰值很近的位置,且交替的从初始点两侧分别取值累加,这样只运算了 5 次就超过了 best-so-far 值,这里就可以提前终止运算并排除当前 C是最优解的可能性。紧接着作者引出他的优化问题,如何选择合适的初始点来做优化,并给出了他的经验:在已经做完归一化的时序子序列中,由于本文讨论的是 Z-Normalization 也就是说归一化操作后的序列均值为0,方差为1 的高斯分布。使用离均值 0 最远距离的点一般应该是对计算距离贡献最大的点,也是起始的最佳位置,并在文章后边的实验中给予了证实。(此处脑补下高斯概率分布)

3.3 Reversing the Query/Data Role in LB_Keogh

在讨论 LB_Keogh 算法时我们是围绕 Q 进行展开的,而不是全部的时间序列 C。这样只需要对 Q 进行一次计算,而不是针对每一个可能的结果进行逐次计算。但有时候对 Q 的计算有可能效果并不理想,这里作者提出了可以 just-in-time 的互换 Q 和 C 的角色来加速计算。

上文讲过 LB_Keogh 的计算复杂度是 O(n) 而 DTW 是动规问题,高于 O(n) 。比起算所有的 C 的 DTW值,这里计算所有 C 的 LB_Keogh 还是比计算 DTW 节省开销的。具体看后边逐层优化,这个优化算法是在比较后面才被用到的。越往后剩下的 C 子序列越少,因为大部分不符合最优解的 C 已经被简单的 LB 方法排除出去了。


3.4 Cascading Lower Bounds

原文作者在50个数据集上使用18种不同的计算 Lower Bound 的方式测试了筛选效果,并绘制了下图。

从图中可以看出每个 LB 算法都可以从计算复杂度和筛选宽松度来评价。在此基础上作者引出了“层级筛选”的概念,即:先使用复杂度较低的算法快速过滤掉一批不符合最优结果的候选 C,再逐渐使用复杂度高且比复杂度低算法更严谨的筛选算法过滤掉剩余数据。准确的说:先使用 O(1) 复杂度的 LB_Kim_FL,然后对 Q 进行 LB_Keogh_EQ 剪枝,此处的算法复杂度可能介于O(1) 与 O(n) 之间,如果该剪枝算法执行完后 LB 仍没有达到 best-so-far 再用 3.3 章节提到的互换 Q 和 C 的 LB_Keogh_EC 算法辅助剪枝效果,最后使用 提前终止 DTW 的算法处理剩余数据。利用这个顺序剪枝可以节省99.9999%的DTW算法的时间开销。

总的来说,这个层级筛选的方法是整个 Paper 的精华所在,之前作者在考虑到计算复杂度和筛选精准的前提下分别介绍了几种高效计算 LB 的方法,最后从18种中选出四个组合效果最好的,并按复杂度从小到大依次就情况使用。

4 总结

通过使用URC优化的 DTW 可以大大提升搜索速度,这方便了大文本时间序列相似度对比在工业中的落地。在写完这个总结时我脑子里从日常的需求中浮现出如下应用场景:

(1)实时系统监控: 当今的分布式系统例如 Hadoop, Kubernetes 等每天都在工业界产生大量时序日志数据。如果可以将这些数据导入分布式时序数据库中的话,可以通过对特定指标的搜索来更好地做日志分析,甚至由于改善后的速度提升做系统的实时监控。

(2)股票 K 线分析:假设可以观察股票 K 线的走势预测未来股票的涨跌(完全非专业想法,不知是否正确)。那么可以通过查找历史数据中某个盈利股票在暴涨前的一段时序数据,将这个数据作为 Q 去搜索当前所有股票 C 与这个股票最相近的一个,如果他们的相似度高,也许可以预测这个股票在未来的一段时间会涨。。。

(3)客服中心 Trend 预测: 以本人所在保险行业为例,客服中心(Call Center)的工作量与时间周期,突发性事件(如:自然灾害,政治变动,传染病)有极大的相关性。一旦某事件发生,会引起 Hotline 在未来的一段时间段打入电话数量暴增,导致需要同时在线的客服人员数量翻倍增长。为了节约成本,优化安排客服人员工作量,可以使用 DTW 这类算法对历史数据进行搜索,如发现某个 Trend 会在未来的一段时间造成工作量非正常的变化,则可以提前组织相应工作以便做到节约成本(如:通过优化客服人员数量、工时)和提高客服质量的目的。

(4)销售 Trend 预测: 方法同 (3),可以预测下个阶段销售数量,以便优化供应链管理。

(5)工业中可以结合分布式流式处理中间件 Kafka 以及 Kafka Streams, Spark Streaming, Flink 等分布式框架,结合时序 NoSQL 数据库实现海量数据高并发实时监控 / 预测。

阅读论文后的思考

DTW 是以计算时序数据的相似度方法进行搜索,来达到对下一阶段时序预测的目的。单从目的上讲,这和机器学习中的生成模型效果类似。比如 HMM,甚至深度学习中的 LSTM 模型也是通过考虑数据在过去时间段上的特征来生成下一时间段的预测。不知用两种方法哪一个会在实战中效果更佳?以及如何将二者(距离计算 VS 概率计算)对比优缺点、适用场景?


由于时间问题本文尚未完成参考资料和图片的引用。。。TBD

由于本人长期生活在德国,中文写作的能力退化了很多。。。我会努力改进的!

参考

  1. ^dtw urc 优化原论文 https://www.cs.ucr.edu/~eamonn/SIGKDD_trillion.pdf
  2. ^时序数据 https://baike.baidu.com/item/%E6%97%B6%E5%BA%8F%E6%95%B0%E6%8D%AE
  3. ^DTW 论文 https://www.aaai.org/Papers/Workshops/1994/WS-94-03/WS94-03-031.pdf
编辑于 2019-10-29

文章被以下专栏收录