Category Archives: 智能运维

一篇关于时间序列异常检测的论文

近期阅读了一篇论文《Rapid Deployment of Anomaly Detection Models for Large Number of Emerging KPI Streams》,这篇文章基于之前的 ROCKA 系统做了一些额外的工作。ROCKA 系统是用来做时间序列的实时聚类的,而这篇文章是在 ROCKA 系统的基础上增加了时间序列异常检测的功能。通常来说,时间序列异常检测可以使用有监督的方法来解决,参考 Opperentice 系统。而本篇文章使用了半监督学习的思路来解决异常检测的问题,下面来详细分析一下这篇文章的细节,本文的作者把这个系统称为 ADS(Anomaly Detection through Self-training)。

数据集的情况

在论文中,作者使用了两份数据集,分别是已经历史上的 70 条时间序列,另外还有新来的 81 条时间序列。在 ADS 系统中,历史上的 70 条时间序列被划分成 5 类,并且已经可以找出每个类的质心位置,并且每条历史上的时间序列通常来说会大于三个星期(3 weeks)。本篇论文的评价指标是 F-Score,也属于机器学习领域里面比较常用的衡量模型效果的指标。整体来看,这篇文章的数据集大约是 200 条时间序列,时间序列的时间间隔通常来说是五分钟(不过其余的运维场景会有一分钟的数据采集粒度),而一般来说都拥有大半年甚至一年的时间跨度。那么时间点的个数预估是 200 * (1440 / 5) * 365。假设异常的数据:正常的数据 = 1:10000(也就是说平均每条时间序列每周至少发生一次异常),于是这批时间序列数据的异常数据量大约是 200 * (1440 / 5) * 365 / 10000 = 2102,也就是说异常的样本大约是 2102 个左右,剩下的都是正常的样本。PS:当然如果异常的数据:正常的数据的比例大于 1:10000 的话,异常的样本还会更多一些。整体来看,时间序列异常检测是一个样本极其不均衡的场景。

ADS 的系统架构

按照作者之前论文的经验,时间序列异常检测通常都是先做聚类,然后再根据每一个类的特点来做一个异常检测模型,之前的技术架构就是 ROCKA + Opperentice。因为 ROCKA 可以根据时间序列的走势和趋势来进行时间序列的实时分类/聚类,然后 Opperentice 就是做时间序列异常检测的模型。在本文的场景下,作者把 70 条时间序列分成了5 类,因此只需要维护 5 个时间序列的异常检测模型就可以了。当然把时间序列切分成更多的类也是可以的,只是需要维护的时间序列异常检测就变多了,人工成本会加大。

ADS系统架构
ADS 的系统架构
ROCKA系统架构
ROCKA 的系统架构

如果看到上面两幅图,有心的读者一定会发现其实 ADS 就是基于 ROCKA 所做的工作。ADS 先对时间序列进行了分类,然后进行了特征提取的工作,再通过半监督学习模型,最后进行异常检测。也就是说,ADS 会走下面四个步骤:

  1. ADS 先把历史上的时间序列进行聚类;
  2. 通过算法获得每一个类的质心,并且标记出质心曲线的异常点和正常点;
  3. 对新来的时间序列进行实时聚类,划分到合适的类别;
  4. 基于新来的时间序列(没有标记)和历史上的时间序列(有标记)使用无监督算法来重新训练一个新的模型,进行该类别的时间序列异常检测。

ADS 的细节

对于时间序列的聚类框架 ROCKA,之前的一篇 BLOG 里面已经详细介绍过,这里将不会再赘述。而 ADS 的另一个模块就是半监督学习算法 Contrastive Pessimistic Likelihood Estimation(CPLE),详细的论文细节可以参考论文《Contrastive Pessimistic Likelihood Estimation for Semi-Supervised Classification》。CPLE 有几个好处:

  1. CPLE 是半监督学习算法中比较健壮的,因为它并没有过多的假设条件,并且也符合这篇文章的业务场景,同时拥有质心曲线(有标记)和新来的曲线(无标记),使用半监督学习也是符合常理的。除了 CPLE,其实在实战过程中也可以多尝试其他的半监督模型,具体可以参考周志华的《机器学习》。
  2. CPLE 的复杂度比较低,计算快。
  3. CPLE 支持增量学习,因此,当越来越多新的时间序列进入 ADS,这个模型也会随之而调整并提高准确率。

整体来看,ADS = ROCKA + CPLE,而在论文中,它的对比模型就是 ROCKA + Opperentice。而且在 CPLE 中,也使用了与 Opperentice 系统类似的特征,如下图所示。

ADS特征
特征工程

其实,从本质上来看,就是半监督学习与有监督学习在这份数据集合上面的比较。从这篇论文里面所展示的数据来看,CPLE 有一定的优势。

ADS效果对比1
Average best F-scores of ADS, iForest, Donut, Opperentice, ROCKA + Opperentice
ADS效果对比2
The New KPI Streams

整体来看,本篇文章介绍了时间序列异常检测的一种方案,也就是把时间序列先进行聚类的操作,然后根据不同的类来进行异常检测。在异常检测的方法中,不仅可以使用 Random Forest,GBDT,XGBoost 等有监督学习方法,也可以使用 CPLE 等半监督算法。具体在业务中如何使用,其实只能够根据具体的数据来进行合理地选择了。

 

Advertisements

两篇关于时间序列的论文

这次整理的就是清华大学裴丹教授所著的两篇与时间序列相关的论文。一篇是关于时间序列聚类的,《Robust and Rapid Clustering of KPIs for Large-Scale Anomaly Detection》;另外一篇文章是关于时间序列异常检测的,重点检测时间序列上下平移的,《Robust and Rapid Adaption for Concept Drift in Software System Anomaly Detection》。本文将会整理一下这两篇文章的关键技术点。

Robust and Rapid Clustering of KPIs for Large-Scale Anomaly Detection

在互联网公司中,通常会拥有海量的的时间序列,而海量的时间序列就有着各种各样的形状和走势。因此,就有学者提出可以先对时间序列进行分类,然后根据不同的类使用不同的检测模型来进行异常检测。如果要做时间序列的分类,就先需要做聚类的操作,无论从 KMeans,DBSCAN,还是层次聚类来说,都会消耗一定的运算时间。所以,如何在较短的时间内进行聚类或者分类的操作则是这个系统的关键之处。于是,这篇文章提出了一个将时间序列快速聚类的方法。

时间序列 -> 时间序列分类

-> 根据每一类时间序列使用不同的异常检测模型

而在做时间序列聚类的时候,也有着不少的挑战。通常挑战来自于以下几点:

  1. 形状:通常来说,时间序列随着业务的变化,节假日效应,变更的发布,将会随着时间的迁移而造成形状的变化。
  2. 噪声:无论是从数据采集的角度,还是系统处理的角度,甚至服务器的角度,都有可能给时间序列带来一定的噪声数据,而噪声是需要处理掉的。
  3. 平移:在定时任务中,有可能由于系统或者人为的原因,时间序列的走势可能会出现一定程度的左右偏移,有可能每天 5:00 起的定时任务由于前序任务的原因而推迟了。
  4. 振幅:通常时间序列都存在一条基线,而不同的时间序列有着不同的振幅,振幅决定了这条时间序列的振荡程度,而振幅或者基线其实也是会随着时间的迁移而变化的。

从整篇论文来看,ROCKA 系统是为了做实时的时间序列分类判断的。要想做成实时的分类判断,就需要有离线和在线两个模块。其中离线是为了做模型训练或者聚类的,在线是为了使用离线处理好的模块来做曲线分类的。

ROCKA系统架构
ROCKA系统架构

从整个系统来看,离线模块需要做以下几件事情:首先需要收集一批时间序列数据,也就是所谓的 Raw Time Series Data(Raw),通过预处理模块,实施基线提取,再进行聚类的操作,获得相应的聚类结果和质心。在线模块同样也要做类似的事情:首先对于每一条新来的时间序列数据,也就是所谓的 New Time Series Data(Raw),通过预处理模块,实施基线提取,然后使用已经聚类好的离线模块来进行实时的分类。

下面,我们来逐一分析每个模块的作用。

预处理模块(Preprocessing)

通常预处理模块包含几个关键点:

  1. 缺失值:缺失值指的是在该上报数据的时间戳上并没有相应的数据上报,数据处于缺失的状态。通常的办法就是把数据补齐,而数据补齐的方法有很多种,最简单的就是使用线性插值(Linear Interpolation)的方式来补齐。
  2. 标准化:对于一个时间序列而言,有可能它的均值是10万,有可能只有10,但是它们的走势有可能都是一样的。所以在这个时候需要进行归一化的操作。最常见的有两种归一化方法:一种是标准化,另外一种是最大最小值归一化。如果 [x_{1},x_{2}, \cdots, x_{n}] 表示原始的时间序列的话,标准化指的是 \hat{x}_{i} = (x_{i} -\mu) /\sigma,其中 \mu\sigma 分别表示均值和标准差。最大最小值归一化指的是 \hat{x}_{i} = (x_{i} - \min) / (\max - \min),其中 \max, \min 分别表示这段时间内的最大值与最小值。

基线提取模块(Baseline Extraction)

基线提取指的是把时间序列分成基线和剩余项两个部分,假设时间序列是 [x_{1},x_{2},\cdots,x_{n}],基线提取就是:

x_{i} = baseline_{i} + residual_{i},

其中 baseline_{i}, residual_{i} 分别指的是 x_{i} 的基线和剩余项。

在处理基线的时候,有几件事情需要注意:

  1. 异常值的处理:通常需要移除一些明显异常的值,然后使用线性插值等方法来把这些移除的值补上。
  2. 使用简单的移动平均算法加上一个窗口值 w 来提取基线。假设时间序列是 [x_{1},\cdots,x_{n}]SMA_{i} = \sum_{j=i-w+1}^{i} x_{j} / wR_{i} = x_{i} - SMA_{i}。也就是说 x_{i} = SMA_{i} + R_{i}
  3. 提取基线的方式其实还有很多,使用带权重的移动平均算法(Weighted Moving Average),指数移动平均算法(Exponentially Weighted Moving Average)都可以提取基线,甚至使用深度学习中的 Autoencoder 或者 VAE 算法都能够提取基线。

基于密度的聚类算法(Clustering)

使用预处理和基线提取技术之后,可以得到原始时间序列的基线值,然后根据这些基线值来进行时间序列的聚类操作。一般来说,时间序列的聚类有两种方法:

  1. 通过特征工程的方法,从时间序列中提取出必要的时间序列特征,然后使用 KMeans 等算法来进行聚类。
  2. 通过相似度计算工具,对比两条时间序列之间的相似度,相似的聚成一类,不相似的就分成两类。

对于第一种方法而言,最重要的就是特征工程;对于第二种方法而言,最重要的是相似度函数的选择。在这篇文章中,作者选择了第二种方法来进行时间序列的聚类。对于两条时间序列 X = [x_{1},\cdots,x_{m}]Y = [y_{1},\cdots,y_{m}] 而言,为了解决左右平移的问题,需要考虑一个偏移量 |s|,然后计算它们之间的内积。

CC_{s}(X,Y) = \begin{cases}  \sum_{i=1}^{m-s}x_{s+i}\cdot y_{i}, \text{ where } s\geq 0, \\ \sum_{i=1}^{m+s}x_{i} \cdot y_{i-s} \text{ where } s<0.\end{cases}

通过这个偏移量 |s| 就可以计算出最大的相似度,然后计算出两条时间序列之间的距离。i.e.

NCC(X,Y) = \max_{s}\bigg(\frac{CC_{s}(X,Y)}{||X||_{2}\cdot||Y||_{2}}\bigg),

SBD(X,Y) = 1 - NCC(X,Y).

其中 NCC \in [-1,1] 指的是 Normalized version of Cross-Correlation,SBD \in [0,2] 指的是 Shape-based distance。

而聚类的另一个重要指标就是质心的选择,在这里,每个类的质心可以设置为:

\text{Centroid} = argmin_{X \in \text{cluster}_{i}} \sum_{Y \in \text{cluster}_{i}} SBD(X,Y)^{2}.

实时分类(Assignment)

对于一条新的时间序列(Raw Data),同样需要经过预处理,基线提取等步骤,然后计算它与之前每一个质心的距离。然后进行距离的从小到大排序,最小的那一类可能就是所需要的。当然,当最小距离大于某个 threshold ,就说明这条新的时间序列曲线很可能不属于之前的任何一类。通过人工查看之后,可以考虑新增一类,并且更新之前所做的聚类模型。

聚类与时间序列异常检测

如果要做海量的时间序列异常检测的话,通常有以下两种做法。

  1. 先对时间序列进行聚类或者分类,针对不同的时间序列类型来使用不同的模型;
  2. 在时间序列异常检测中加入分类特征。

对于第一种方法而言,需要针对不同形状的时间序列维护不同的模型,而且如果第一层的聚类/分类错误了,那么使用的模型也会出现错误。对于第二种方法而言,关键在于样本的积累和分类特征的构建。

Robust and Rapid Adaption for Concept Drift in Software System Anomaly Detection

这篇文章是关于时间序列异常检测的,而清华大学的 Netman 实验室做时间序列异常检测的相关工作比较多。从 2015 年开始的 Opperentice(Random Forest 模型),Funnel(SST 模型,变更相关),到后来的 Donut(VAE 模型),都是时间序列异常检测的相关文章。而这篇问题提到的 StepWise 系统针对的场景是关于指标的迁移的,所谓概念漂移(Concept Drift)指的就是时间序列会随着变更,发布,调度或者其他事件的发生而导致上下漂移或者左右漂移。一个比较典型的例子就是关于网络流量的漂移,通常来说在某个特点的时间点,时间序列会出现猛跌或者猛涨的情况,但是下跌或者上涨之后的走势和历史数据是及其相似的,只是绝对值上有所变化。

StepWiseExample.png
概念漂移的例子

正如上图所示,在时间戳 08-02 附近,时间序列出现了一个下跌的情况。但是根据历史数据所计算出来的上界(Upper Bound)和下界(Lower Bound)却会把未来一段时间的序列都判断为异常。但是前后的曲线走势却是一样的,其实只有下跌的那一小段时间应该被判断为异常,其他时间都是正常的。基于以上的业务场景,这篇文章的作者就提出了概念漂移(Concept Drift)的一些方法。

在整个 StepWise 系统背后,有两个比较关键的地方。其中第一个关键之处就是如何判断异常点。在这个业务场景下,需要检测的是上图发生暴跌的点。而暴跌或者暴涨只是业务运维用来描述时间序列的一个词语。在统计学领域,这种点通常称为变点(Change Point)。所以概念漂移的检测可以转化为一个变点检测的问题,正是论文里面写的。

StepWise系统架构
StepWise 的系统架构

Insight 1:Concept drift detection can be converted into a spike detection problem in the change score stream.

如果是进行变点检测的话,其实可以参考 2016 年的论文 Funnel:Assessing Software Changes in Web-based Services。里面使用了 Singular Spectrum Transform 来进行检测,并且使用 Difference in Difference 来判断时间序列是否真正出现了异常。关于 SST(Singular Spectrum Transform)的具体细节可以参考 Yasser Mohammad, Toyoaki Nishida 所撰写的 Robust Singular Spectrum Transform。SST 算法是基于 SVD 算法的,有一定的时间复杂度。而基于 SST 又有学者提出了 Robust SST,具体可以详见论文。而 StepWise 的其中两步就是 Detection 和 Classification,其中的 Detection 使用了 Improved SST,Classification 使用了Difference in Difference。而 Funnel 系统的其中两步也是 Improved SST 和 Difference in Difference,这与 Funnel 有异曲同工之妙,Funnel 系统详情请见下图。

Funnel系统架构
FUNNEL系统架构

而突变前后的时间序列可以分别叫做 Old Concept 和 New Concept,由于前后的走势几乎一致,所以本文的第二个关键点就是相似度的判断。

Insight 2:The relationship between new concept and old concept can be almost perfectly fitted by linear regression.

从数学的角度来讲,Insight 2 的陈述是想判断两条时间序列是否相似。假设 X=[x_{1},\cdots, x_{n}]Y = [y_{1},\cdots, y_{n}],那么以下陈述是等价的。

  1. 存在一个线性关系 y = kx 使得二维点集 \{(x_{i},y_{i}), 1\leq i\leq n\} 能够被很好的拟合好,也就是说此刻的方差较小。
  2. [x_{1},\cdots,x_{n}][y_{1},\cdots,y_{n}] 的 Pearson 系数很高;
  3. [x_{1},\cdots,x_{n}][y_{1},\cdots,y_{n}] 标准化(z-score)之后的曲线几乎一致。
  4. 分别存在两个值 \mu_{1},\mu_{2} 使得 [x_{1}/\mu_{1},\cdots,x_{n}/\mu_{1}][y_{1}/\mu_{2},\cdots,y_{n}/\mu_{2}] 几乎一致。

所以,在 StepWise 中,作者通过线性回归的算法来判断 old concept 与 new concept 是否存在线性关系,也就是说这两者是否只是平移的关系。其实也可以尝试上面所说的其余方法。

整体来看,StepWise 大致分成两个关键步骤。首先使用类似 Funnel 系统的思想,先进行异常检测(SST),再使用判断一下是不是因为软件的变更引起的(DID),最后使用一层过滤逻辑(线性拟合)来判断是否出现了概念漂移的情况。在实际使用过程中,无论是异常检测还是过滤逻辑,都需要根据具体的业务来做,很难找到一个固定的模型来解决所有的难题。

 

浅谈机器学习的一些小众方向

随着 DeepMind 的 AlphaGo 在 2016 年战胜了李世石,“人工智能”这个词开始进入大众的视野。从那时起,不管是大型互联网公司还是初创企业都开始大规模招聘机器学习的相关从业者,无论社招的求职者还是校招的应聘学生都出现了大规模的增长。由于机器学习的人才短缺并且大量应届生涌入,以至于现在某些公司的校园招聘出现了算法工程师简历太多,并且移动端岗位,web 开发岗位的简历略有不足的情况,导致这些互联网公司甚至通过邮件的方式来呼吁应届生尽量修改投递职位。

字节跳动

就这几年的人工智能发展情况和笔者的个人经验而言,人工智能可以大致分成以下几个方向:

  1. 计算机视觉方向
  2. 自然语言处理方向
  3. 语音识别方向
  4. 机器学习方向

CV,NLP & Speech Recognition

计算机视觉方向(Computer Vision)无论是在学校还是在公司,都有着大量的从业者,并且 ImageNet 项目可以提供上千万的标注图片供大家使用。既然 ImageNet 是开源的数据集,那么无论是学校的教授还是学生,不管是大型互联网公司还是初创企业,都可以轻易地获取到这些数据集,不仅可以进行 CV 算法的研究工作,还可以进行相关的工程实践。由于计算机视觉方向的历史悠久,不管是计算机系,工程系,甚至数学系,都有着大量的老师和相应的学生从事该方向的研究工作,因此,学校或者研究所能够对工业界输出的计算机视觉人才数量也是可观的。

与计算机视觉方向相比,自然语言处理方向(Natural Language Processing)在学校里面也有不少的教授从事相关研究。不过要想让计算机理解人类的语言可不是一件容易的事情。尤其是中文还拥有多音字,语义双关等情形,而且理解中文很可能还要基于上下文来前后推敲。如果和聊天机器人聊过就会发现,其实聊天机器人和人类的聊天给用户的感觉是完全不一样的。语音方向笔者不是很了解,也只是道听途说而已,在这里就不在赘述了。

ImageNet

机器学习

除了以上三个方向,人工智能的另外一个研究方向自然就是机器学习了。在周志华老师的教材《机器学习》中,无监督学习,有监督学习,半监督学习,强化学习等方向都已经在该教材中进行了详细的解释。貌似几年前强化学习这个方向也是不温不火,但是在 AlphaGo 崛起之后,深度学习和强化学习就已经开始进入了大多数人的视野。随着围棋被攻克之后,德州扑克AI,或者其他的游戏 AI 也被很多学者和大型游戏公司所关注。DeepMind 也在 2017 年开放了星际争霸的研究平台,今年无论是在 Dota2 还是星际争霸上,游戏 AI 相比之前都有了巨大的突破。

starcraft2

除了强化学习之下的游戏 AI 之外,其实机器学习一直在一个领域发挥着巨大的用处,那就是推荐系统。无论是广告推荐,YouTube 视频推荐,甚至今年非常火的抖音 APP,推荐系统在其中的作用都不容忽视。关于推荐系统的书其实有很多,笔者也没有一一读过,不过就近些年的发展状况来看,无论是在学术界还是工业界,从零到一搭建一套推荐系统已经不是壁垒,如何搭建一套结合业务场景的优秀推荐系统才是难题。而推荐系统中常用的各种模型,例如逻辑回归(logistic regression),SVD,ItemCF & UserCF,甚至深度神经网络,在各种开源框架之下(Spark,Tensorflow等),只要拥有足够的计算资源,训练出一个可以使用的模型已经没有太大的难度。难度在于算法工程师如何贴近业务并且理解业务,在此基础上如何使用机器学习算法将内容库里面的优质内容推荐给用户,而不引起用户的反感,点击率如何在合理的范围内进一步提升。搭建一套推荐系统已经不是难题,如何结合多种多样的推荐场景才是关键,怎么结合业务来使用推荐系统则是算法工程师需要思考的问题。

Tensorflow

机器学习+安全业务

就笔者的个人经验来看,推荐系统或者游戏 AI 其实只是机器学习的一个应用领域。既然机器学习能够应用在推荐系统或者游戏 AI 上,那么为何不能够应用在别的领域上呢?

对于一些大型互联网公司而言,推荐系统能够给用户们带来足够优质的体验,游戏 AI 能够帮助玩家提升自己的技艺。但是在给用户带来优质体验的时候,总有一些黑产用户在伺机而动,通过 APP 的各种 bug 来寻找赚钱的机会,给正常用户带来各种各样的骚扰。在游戏中,有一些人使用了外挂等技术,破坏了游戏中的平衡。在金融行业中,一直都有黑产用户正在进行各种各样违法犯罪的事情,例如信用卡欺诈等,给正常用户带来了不少的损失。在社交网络中,有一些用户通过社交网络传播着各种各样的不良信息,无论是谣言,虚假广告还是各种各样的假冒伪劣产品宣传,都给正常用户带来了不好的体验。因此,安全业务一直是互联网公司和金融公司的重点业务,安全业务一直是保护着互联网公司能够正常运行的基石。各种各样的安全实验室在大型互联网公司里面并不罕见,也是必须要配备的力量。对于业务安全上,无论是盗号,刷帖,传播虚假消息等都是需要关注的对象。在黑产力量日益壮大的情况下,打击黑产的人力也越来越多。随着人力的增多,如何使用机器学习算法来进行人类经验的传承,或者说随着黑产技术的升级如何才能够尽快的提升互联网公司的黑产对抗能力,这些都是值得做的工作。除了互联网公司之外,银行等金融机构也需要进行信用卡的风控评级,打击信用卡盗刷,黑色产业的资金链条挖掘等。因此,银行等金融机构对于业务安全上面的要求有的时候可能比互联网公司还要严格。

黑客图片1

能够用在安全领域上的机器学习算法有很多,最容易想到的当然就是异常检测。无论是高维异常检测,还是图(Graph)上的异常检测,都在业务安全领域有着巨大的应用场景。异常检测算法可以从众多的数据中发现数据中的异常点,然后通过人工审核等方式进行数据的标注,并且可以使用有监督学习模型进行训练和上线预测。整体来说,就是使用无监督算法,有监督算法,图挖掘算法等机器学习常见技术来进行恶意黑产的打击工作。对于从事业务安全+机器学习方向的算法工程师来说有一些潜在的优势,那就是业务安全方向是工业界的刚需。但是学术界并不完全有能力培养相关的人才,因为互联网或者金融公司的数据都具有保密性,很难把数据像 ImageNet 一样开放给全世界,共同享受数据带来的巨大优势。如果没有基础的数据,那么学校的教授或者学生就无法接触到这个领域,也就无法在学校提升相关的技术。虽然异常检测等其他机器学习算法会在学术中有所突破,但是安全的业务经验只有在做过相关业务之后,真正地打击过黑产用户之后才能够有更深层次的体会和理解。一个没有接触过安全业务的人,即使他的学术造诣再高,在短时间内也是很难提出一些靠谱想法或者技术方案的。

机器学习+运维业务

在这里做一个不恰当的比喻来方便大家理解。

如果把 APP 比喻成一栋楼房的话,那么后台开发就是搭建钢筋水泥的人,前台开发就是负责刷墙贴砖的人,设计师是负责把这栋楼设计得更加美观的人,安全人员就好比楼房的保卫人员,那么运维人员就是这栋大楼的检修人员。

在一些互联网公司,运维人员也被称为技术运营人员,整体来说就是保障APP或者业务稳定运营的。例如:网络抖动了该怎么办,交换机何时宕机,大量用户无法登陆APP了该怎么办,APP的某个页面无法打开了该怎么办等诸如此类的问题。为了保障业务的稳定运营,就需要有一定数量的技术运营同事来维护整个业务的正常运行。正所谓“天有不测风云,人有旦夕祸福”,公司拥有安全人员和运维人员好比买保险,在没有黑客攻击或者业务正常运行的时候,通常存在感略低。但是一旦业务出了问题,第一个要召集的人肯定就是安全和运维人员。因此,无论是安全工作还是运维工作,都是大型互联网公司和金融机构必不可少的力量。

随着机器学习的发展,智能运维(Artificial Intelligence Operations),也就是所谓的 AIOps,也开始被众多技术公司所关注。提到技术运营工作,根据 2018 年的《企业级AIOps实施建议白皮书V0.6》 的观点,可以大致分成以下三个方向:

  1. 质量保障;
  2. 效率提升;
  3. 成本管理。

其中质量保障就是为了保证业务的正常,高效,稳定地运转。在质量保障的过程中,无法避免的就需要进行异常检测。在运维领域,异常检测的范围非常广,不仅包括大家耳熟能详的时间序列异常检测,还包括多维数据下钻分析,甚至还包括日志模板提取和异常挖掘。除了质量保障之外,效率提升也是一个方面,无论是自动化运维领域还是使用 NLP 的技术来构建智能聊天机器人,甚至使用机器学习等技术来进行智能扩缩容,机器学习技术在运维领域都有着巨大的发挥空间。

AIOps场景

在智能运维领域,最重要的任务之一就是时间序列异常检测,这里的时间序列不仅包括服务器的各种各样的指标(CPU,进程,PKG等),还有网络出入流量等交换机数据,甚至包括各种各样的业务指标(在线用户数,失败数,请求量等)。各种各样的时间序列组合在一起就形成了一个时间序列数据库,而且这些时间序列通常来说都是按照分钟量级来收集数据,因此,时间序列项目完全符合机器学习项目的各种条件。在时间序列异常检测或者趋势预测中,时间序列和机器学习,甚至深度学习结合的各种技术都可以在这里有着一定的用武之地。

timeseries

除了时间序列之外,服务器的异常挖掘,多维度数据分析都是智能运维中非常有挑战的项目。除了质量保障之外,效率提升中的智能聊天机器人将有希望把运维人员从繁重的客服任务中解放出来,智能扩缩容技术将有机会取代原来很多“拍脑袋”所做出来的容量估计。对于一家正常经营的公司而言,质量保障和效率提升只是其中的两个方面,如何有效地进行成本的管理则是非常重要的项目。如果成本预算过少,那么明年的项目发展将会受到限制;如果成本预算过多,那么明年的资源势必造成各种浪费。因此,无论是质量保障,效率提升,还是成本管理,都是技术运营领域的核心问题。

成本

机器学习+其他领域

除了以上笔者接触过或者略微了解过的领域之外,其实机器学习在其他的领域应该都是有着自己的用武之地。在量化分析方向,据说有的团队已经开始用机器学习的方法进行股票交易。在化学或者生物学领域,也有学者使用机器学习的方法来挖掘数据之间的信息。总之,除了人工智能在那几个经典领域的应用之外,机器学习的方法应该有希望应用到各行各业中,改变原来的工作方式,提升原有学科的效率。机器学习本身并不是一个新的东西,只要运用得当,机器学习在各行各业都有着强大的创造力和生命力。

 

基于自编码器的时间序列异常检测算法

随着深度学习的发展,word2vec 等技术的兴起,无论是 NLP 中的词语,句子还是段落,都有着各种各样的嵌入形式,也就是把词语,句子,段落等内容转换成一个欧氏空间中的向量。然后使用机器学习的方法来进行文本的聚类和相似度的提取,甚至进行情感分类等操作。那么在表示学习(Representation Learning)方向上,除了刚刚提到的自然语言之外,语音,图像,甚至图论中的Graph都可以进行嵌入的操作,于是就有了各种各样的表示算法。既然提到了表示学习,或者特征提取的方法,而且在标注较少的情况下,各种无监督的特征提取算法就有着自己的用武之地。除了 NLP 中的 word2vec 之外,自编码器(Auto Encoder)也是一种无监督的数据压缩算法,或者说特征提取算法。本文将会从自编码器的基础内容出发,在时间序列的业务场景下,逐步展开基于自编码器的时间序列表示方法,并且最终如何应用与时间序列异常检测上。

自编码器

AutoEncoder3

提到自编码器(Auto Encoder),其实它就是一种数据压缩算法或者特征提取算法。自编码器包含两个部分,分别是编码层(encoder)和解码层(decoder),分别可以使用 \phi\psi 来表示,也就是说:

\phi: X\rightarrow F,

\psi: F\rightarrow X,

\phi,\psi = argmin_{\phi,\psi}||X-(\psi\circ\phi)X||^{2},

其目标函数就是为了拟合一个恒等函数。对于最简单的情况,可以令 X = \mathbb{R}^{n}, F=\mathbb{R}^{m},并且编码器和解码器都是前馈神经网络,也就是说:

z = f(Ax+c),

x'=g(Bx+d),

损失函数就是 L(x,x')=||x-x'||^{2} = ||x-g(Bf(Ax+c)+d)||^{2}, 其中 x\in X=\mathbb{R}^{n}, z\in F =\mathbb{R}^{m}. fg 分别是编码层和解码层的激活函数,A,cB,d 分别是编码层和解码层的矩阵和相应的向量。具体来说它们的矩阵大小分别是 A_{m\times n}, c_{m\times 1}, B_{n\times m}, d_{n\times 1}.

AutoEncoder2

对于自编码器而言,它的输入层的维度等于输出层的维度,隐藏层的维度是需要小于输入层的维度的。只有这样,自编码器才可以学习到数据分布的最显著特征。如果隐藏层的维度大于或者等于输入层的维度,其实是没有任何意义的,具体的解释可以参考下面这个Claim。

Claim. 对于自编码器而言,其中隐藏层的维度 m 一定是要小于输入层的维度 n 的。

Proof. 如果 n=m,那么令 A=B=I_{n}, c=d=0, f=g=id 就可以得到一个自编码器,而这个自编码器对于提取特征没有任何的意义。同理,当 m>n 时,A 是一个 m\times n 矩阵,B 是一个 n\times m 矩阵。从线性代数的角度来看,有无数个矩阵 A, B 满足 BA=I_{n}。这种情况下对于提取特征也是没有意义的。而当 m<n 时,其实无法找到矩阵 A,B 使得 BA=I_{n}. 如果存在 BA=I_{n}, 那么

n = rank(I_{n})=rank(BA) \leq \min\{rank(A),rank(B)\} \leq \min\{m,m\}=m.

这就导致了矛盾。因此,只有在 m<n 的情况下提取特征才是有意义的。

对于自编码器而言,其本质上也是一个神经网络,那么它的激活函数其实不仅可以选择 sigmoid, 还可以使用 tanh,ReLU,LeakyReLU 等其余激活函数,其本质上都是为了拟合一个恒等变换,中间层则作为一个特征提取的工具。在训练的时候,同样是使用反向传播算法,可以使用不同的优化函数,例如 SGD,Momentum,AdaGrad,RMSProp,Adam 等。

在图像领域,有学者尝试使用自编码器来进行图像的重构工作,图像的特征提取等内容,整体来看也能达到不错的效果,请看下图:

AutoEncoder1

从上图来看,基于均方误差的自编码器是无法重构出乒乓球的。由于该自编码器的容量有限,目标函数是均方误差,因此自编码器并没有意识到乒乓球是图片中的一个重要物品。

 

时间序列异常检测:

时间序列异常检测一直是学术界和工业界都关注的问题,无论使用传统的 Holt-Winters,ARIMA,还是有监督算法进行异常检测,都是统计学和传统机器学习的范畴。那么随着深度学习的兴起,是否存在某种深度学习算法来进行异常检测呢?其实是存在的。请看上图,左边一幅图有一个白色的小乒乓球,但是随着自编码器进行重构了之后,白色的小乒乓球已经在重构的图像中消失了。那么根据异常检测的观点来看,小乒乓球其实就可以作为图片中的异常点。只要在图片的局部,重构出来的图片和之前的图片存在着巨大的误差,那么原始图片上的点就有理由认为是异常点。

在这个思想下,针对时间序列异常检测而言,异常对于正常来说其实是少数。如果我们使用自编码器重构出来的时间序列跟之前有所差异的话,其实我们就有理由认为当前的时间序列存在了异常。其实,简单来看,基于自编码器的时间序列异常检测算法就是这样的:

原始时间序列

-> Auto Encoder(Encoder 和 Decoder)

-> 重构后的时间序列

-> 通过重构后的时间序列与原始时间序列的整体误差和局部误差来判断异常点

简单来说,只要输出的时间序列在局部的信息跟原始的时间序列不太一致,就有理由认为原始的时间序列存在着异常。

那么,首先我们需要提取时间序列中的一些子序列,例如我们可以提取今天(today),昨天(yesterday),一周前(week)的数据,基于同样的时间戳把它们重叠在一起,也就是下图这个形式。其中,蓝线表示一周前的数据,黑线表示昨天的数据,红色表示今天的数据。

AutoEncoder4

基于一条很长的时间序列,我们可以提取它的很多子序列,从而构造出很多的片段序列。这些片段序列就可以形成自编码器的输入数据,而自编码器是模拟一个恒等变换,因此它会把有异常的点尽量磨平,而正常的点则保持原样。所以,通过大量子片段来进行训练数据的输入,自编码器就能够得到一个较为合理的权重。得到了一个训练好的自编码器之后,对于任何一个子片段,都可以重构出一个新的片段。例如上面的子片段就可以重构成下图:对于今天的数据(today),那个凸起被直接抹平;对于昨天的数据(yesterday)而言,那个凹下去的部分也被磨平。基于时间序列重构前和重构后的数据差异,可以获得时间序列的异常点。

AutoEncoder5

除此之外,还有很多时间序列的异常点可以被自编码器(AutoEncoder)发现,例如下面四幅图,无论是上涨,还是下跌,其实都可以被自编码器(AutoEncoder)发现异常。

总结

通常来说,在时间序列异常检测场景中,异常的比例相对于正常的比例而言都是非常稀少的。因此,除了有监督算法(分类,回归)之外,基于无监督算法的异常检测算法也是必不可少的。除了 HoltWinters,ARIMA 等算法之外,本文尝试了一种新的异常检测算法,基于深度学习模型,利用自编码器的重构误差和局部误差,针对时间序列的异常检测的场景,初步达到了一个还不错的效果。这种方法可以用来提供部分异常样本,加大异常检测召回率的作用。但是这种方法也有一定的弊端:

  1. 从理论上说,它只能对一个时间序列单独训练一个模型,不同类型的时间序列需要使用不同的模型。这样的话,其实维护模型的成本比较高,不太适用于大规模的时间序列异常检测场景;
  2. 对周期型的曲线效果比较好,如果是毛刺型的数据,有可能就不太适用;因为长期的毛刺型数据就可以看成正常的数据了。
  3. 每次调参需要人为设置一定的阈值,不同的时间序列所需要的阈值是不一样的。

参考文献

  1. Unsupervised Anomaly Detection via Variational Auto-Encoder for Seasonal KPIs in Web Applications, Haowen XU, etc., 2018
  2. Deep Learning, Ian Goodfellow, etc., 2016
  3. https://zr9558.com/2016/06/12/replicator-neural-networks/

 

黑盒函数的探索

黑盒函数的定义

在工程上和实际场景中,黑盒函数(black box function)指的是只知道输入和输出对应关系,而不知道它的内部结构,同时也不能写出具体表达式的一类函数。正如下图所示,每次给定一组输入,通过黑盒函数的计算,都能够得到一组输出的值,但是却无法写出 Black box 函数的精确表达式。

Blackbox3D-withGraphs

与之相反的是函数或者系统称之为白盒函数(open system),它不仅能够根据具体的输入获得相应的输出,还能够知道该函数的具体表达式和细节描述。例如 f(x) = \sin(x)f(x) = \ln(x) 等都是白盒函数。

黑盒函数的研究对象

无论是白盒函数还是黑盒函数,都有很多的学术界人士和工业界人士去研究。通常来说,对于一个函数 f(x) 而言,我们都可以研究该函数的以下性质:

  1. 最大值与最小值,i.e. \max f(x)\min f(x).
  2. 根,i.e. \{x:f(x) = 0\}.
  3. 函数的单调性与凹凸性等。

对于具有明显表达式的函数,例如 f(x) = \sin(x) 等,我们能够使用的方法和技巧都很多,其方法包括但不限于导数,积分,Taylor 级数等等。但是对于黑盒函数,我们能够使用的方法和技巧就会一定的限制。本文将从如何研究一个函数的根,最大值和最小值等方向入手,逐步向大家展示黑盒函数研究中所遇到的数学与机器学习方法。

黑盒函数的根

对于多项式 p(x) = a_{n}x^{n}+\cdots + a_{0} 而言,多项式的根指的是使得 p(x)=0x 的解。特别的,对于二次多项式而言,也就是 ax^{2} +bx+c=0, 它的根可以表示为:

x = \frac{-b\pm\sqrt{b^{2}-4ac}}{2a}.

对于一般函数 f(x) 而言,它的根指的是 \{x:f(x)=0\} 这个集合。下面我们来介绍一下如何计算一个函数的根。

二分法

在数学分析中,介值定理(Intermediate value theorem)描述了连续函数在两点之间的连续性,具体描述是:

[介值定理] 如果连续函数 f(x) 的定义域包含 [a,b], 而且通过 (a,f(a))(b,f(b)) 两点,它也必定通过区间 [a,b] 内的任意一点 (c,f(c)), 其中 a<c<b.

从介值定理可以得到,如果我们知道 f(x_{1})<0f(x_{2})>0, 那么必定存在 x_{0} \in (x_{1},x_{2}) 使得 f(x_{0})=0。根据这个定理,我们可以提出二分法来计算函数的根。

如果要计算 f(x) = 0 的解,其一般步骤是:

  1. 先找到一个区间 [a,b], 使得 f(a)f(b)<0;
  2. 求这个区间的中点 m=(a+b)/2, 并求出 f(m) 的取值;
  3. 如果 f(m)=0, 那么 m 就是函数的根;如果 f(m)f(a)>0, 就选择 [m,b] 为新的区间,否则选择 [a,m] 为新的区间;
  4. 重复第 2 步和第 3 步直到达到最大迭代次数或者最理想的精度为止。

 

Bisection_method

 

牛顿法(Newton’s Method)

牛顿法的大致思想是:选择一个接近 f(x) 零点的初始点 x_{0}, 计算这个点相应的函数取值 f(x_{0}) 与导数值 f'(x_{0}), 然后写出通过点 (x_{0},f(x_{0})) 的切线方程,并且计算出该切线与横轴的交点 x_{1}. i.e.

x_{1} = x_{0} - f(x_{0})/f'(x_{0}).

我们可以不停地重复以上过程,就得到一个递推公式:

x_{n+1}= x_{n}-f(x_{n})/f'(x_{n}).

在一定的条件下,牛顿法必定收敛。也就是说 x_{n} 随着 n 趋近于无穷,将会趋近于 f(x)=0 的解。

Ganzhi001

割线法

根据导数的定义:

f'(x_{0}) = \lim_{x\rightarrow x_{0}}\frac{f(x)-f(x_{0})}{x-x_{0}}

可以得到,当 x 靠近 x_{0} 的时候,可以用右侧的式子来估计导数值,i.e.

f'(x_{0}) \approx \frac{f(x)-f(x_{0})}{x-x_{0}}.

当我们不能够计算 f(x) 的导数的时候,就可以用上式来代替。

于是,割线法与牛顿法的迭代公式非常相似,写出来就是:

x_{n+1} = x_{n} - \frac{x_{n}-x_{n-1}}{f(x_{n})-f(x_{n-1})}f(x_{n}).

在这里,割线法需要两个初始值 x_{0}x_{1}, 并且它们距离函数 f(x)=0 的根越近越好。

割线法1

备注

对于黑盒函数而言,我们是不知道它们的表达式的,因此以上的方法和技巧在黑盒函数的使用上就有限制。例如牛顿法是需要计算函数的导数值的,因此不适用在这个场景下。但是对于二分法与割线法,只需要计算函数在某个点的取值即可,因此可以用来寻找黑盒函数的根。

黑盒函数的最大值与最小值

对于能够写出表达式的函数而言,如果要寻找 f(x) 的最大值与最小值,可以计算 f(x) 的导数 f'(x), 然后令 f'(x) =0, 就可以得到函数的临界点(critical point),再根据周围的点导数的性质即可判断这个点是否是局部最大值或者局部最小值。

Weierstrass 逼近定理

对于黑盒函数而言,通常来说我们只知道一组输入和相应的输出值。如果只考虑一维的情况而言,那就是 \{(x_{i},y_{i})\in \mathbb{R}^{2},0\leq i\leq n\}n+1 个样本。根据 Weierstrass 逼近定理可以知道:

  1. 闭区间上的连续函数可以用多项式级数一致逼近;
  2. 闭区间上的周期为 2\pi 的连续函数可以用三角函数级数一致逼近。

用数学符号来描述就是:

[Weierstrass 逼近定理] 假设 f(x) 是闭区间 [a,b] 连续的实函数。对于任意的 \epsilon>0,存在一个多项式 p(x) 使得对于任意的 x\in[a,b],|f(x)-p(x)|<\epsilon.

因此,如果要计算黑盒函数的最大值和最小值,可以使用一个多项式去拟合这 n+1 个点,于是计算这个多项式的最大值与最小值即可。

Lagrange 插值公式

按照之前的符号,如果我们拥有 n+1 个样本 \{(x_{i},y_{i}), 0\leq i\leq n\}, 那么我们可以找到一个多项式 p(x) 使得 p(x_{i}) = y_{i} 对每一个 0\leq i\leq n 都成立。根据计算,可以得到该多项式是:

p(x) = \sum_{i=0}^{n}\bigg(\prod_{0\leq j\leq n, j\neq i}\frac{x-x_{j}}{x_{i}-x_{j}}\bigg)y_{i}.

于是,要计算黑盒函数的最大值与最小值,就可以转化成计算 p(x) 的最大值与最小值。

除了数学方法之外,机器学习中有一种算法叫做启发式优化算法,也是用来计算黑盒函数的最大值和最小值的,例如粒子群算法与模拟退火算法等。

粒子群算法(Particle Swarm Optimization)

PSO 最初是为了模拟鸟群等动物的群体运动而形成的一种优化算法。PSO 算法是假设有一个粒子群,根据群体粒子和单个粒子的最优效果,来调整每一个粒子的下一步行动方向。假设粒子的个数是 N_{p},每一个粒子 \bold{x}_{i}\in \mathbb{R}^{n} 都是 n 维欧几里德空间里面的点,同时需要假设粒子的速度 \bold{v}_{i}\in\mathbb{R}^{n}。在每一轮迭代中,需要更新两个最值,分别是每一个粒子在历史上的最优值和所有粒子在历史上的最优值,分别记为 \bold{x}_{i}^{*}1\leq i \leq N_{p})和 \bold{x}^{g}。在第 t+1 次迭代的时候,

\bold{v}_{i}(t+1) = \bold{v}_{i}(t) + c r_{1}[\bold{x}_{i}^{*}(t) - \bold{x}_{i}(t)] + c r_{2}[\bold{x}^{g}(t) - \bold{x}_{i}(t)],

\bold{x}_{i}(t+1) = \bold{x}_{i}(t)+\bold{v}_{i}(t+1), \text{ } 1\leq i\leq N_{p}.

在这里,c>0,并且 r_{1}, r_{2}[0,1] 中间的随机数。

模拟退火(Simulated Annealing)

SA1SA2SA3

模拟退火算法是为了模拟金属退火现象。其核心思想是构造了温度 T 这个维度,当温度 T 高的时候,分子运动速度快,能够探索更多的区域;温度 T 低的时候,分子运动速度慢,能够探索的区域有限。随着时间的迁移,温度 T 会从一个较大的值逐渐降低到 0。

模拟退火的大体思路是这样的:先设置一个较大的温度值 T,随机初始化 \bold{x}(0)。假设目标函数是 E(\bold{x}), 需要寻找 argmin_{\bold{x}}E(\bold{x}),然后执行以下的程序:

Repeat:

a. Repeat:

i. 进行一个随机扰动 \bold{x} = \bold{x} + \Delta \bold{x}

ii. 计算 \Delta E(\bold{x}) = E(\bold{x}+\Delta\bold{x}) - E(\bold{x})

如果 \Delta E(\bold{x}) <0,也就是 E(\bold{x} + \Delta\bold{x})<E(\bold{x}),选择 \bold{x}+\Delta\bold{x}

否则,按照 P = e^{-\Delta E/T} 的概率来接受 \bold{x} +\Delta\bold{x}

b. 令 T = T-\Delta T

直到 T 足够小。

总结

本文从数学和机器学习的角度,简要介绍了部分计算黑盒函数的最大值,最小值和根的方法,后续将会介绍更多的类似方法。

 

时间序列的自回归模型—从线性代数的角度来看

Fibonacci 序列

在开始讲解时间序列的自回归模型(AutoRegression Model)之前,我们需要先介绍一下线性代数的基础知识。为了介绍线性代数的基础知识,我们可以先看一个简单的例子。考虑整数序列 Fibonacci 序列,它的通项公式是一个递归函数,每项的取值是由前两项所生成的,其具体的公式就是

F_{n}=F_{n-1}+F_{n-2}

其初始值是 F_{0}=0, F_{1} = 1。按照其递归公式来计算,我们可以详细写出前面的几项,那就是:

0,1,1,2,3,5,8,13,21,34,55,89,144,…

但是,计算 Fibonacci 的通项公式则比计算等差数列和等比数列的通项公式复杂的多,因为这里需要使用线性代数的技巧才能解决这个问题。

求解 Fibonacci 序列的通项公式 --- 矩阵对角化

根据 Fibonacci 数列的递归公式,基于矩阵乘法的定义,Fibonacci 序列可以写成如下形式:

\left( \begin{array}{c} F_{n+2}  \\ F_{n+1}  \\ \end{array}\right)= \left( \begin{array}{cc} 1 & 1 \\ 1 & 0 \\ \end{array}\right) \left( \begin{array}{c} F_{n+1}  \\ F_{n}  \\ \end{array}\right) = A \left( \begin{array}{c} F_{n+1}  \\ F_{n}  \\ \end{array}\right),

\left( \begin{array}{c} F_{n}  \\ F_{n-1}  \\ \end{array}\right)= A \left( \begin{array}{c} F_{n-1}  \\ F_{n-2}  \\ \end{array}\right) = \cdots = A^{n-1} \left( \begin{array}{c} F_{1}  \\ F_{0}  \\ \end{array}\right).

这就意味着我们需要计算出矩阵 A 的幂。在线性代数里面,为了计算矩阵的 n 次方,除了通过矩阵乘法的公式直接计算之外,还有一个经典的技巧,那就是将矩阵对角化。详细来说,如果 m\times m 的矩阵 A 能够对角化,那就是存在可逆矩阵 P 使得

P^{-1}AP = diag(\lambda_{1},\cdots,\lambda_{m})

\implies AP = P diag(\lambda_{1},\cdots,\lambda_{m})

\implies A = P diag(\lambda_{1},\cdots,\lambda_{m})P^{-1}.

其中 D = diag(\lambda_{1},\cdots,\lambda_{m}) 表示一个 m\times m 的对角矩阵,其对角的元素从左上角到右下角依次是 \lambda_{1},\cdots,\lambda_{m}。如果把矩阵 P 写成列向量的形式,i.e. P =(\vec{\alpha}_{1},\cdots,\vec{\alpha}_{m}),那么以上的矩阵方程就可以转换为 A\vec{\alpha}_{i} = \lambda_{i}\vec{\alpha}_{i}, 1\leq i\leq m。进一步来说,如果要计算矩阵 A 的幂,就可以得到:

A^{k} = (PDP^{-1})\cdots(PDP^{-1}) = P D^{k} P^{-1}= P diag(\lambda_{1}^{k},\cdots,\lambda_{m}^{k})P^{-1}.

另外,如果想知道一个矩阵 A 的特征值和特征向量,则需要计算以下多项式的根,i.e. 计算关于 \lambda 的多项式的解,

det(\lambda I - A) = 0,

其中 I 是一个单位矩阵(identity matrix)。

按照以上的思路,如果令

A = \left( \begin{array}{cc} 1 & 1 \\ 1 & 0 \\ \end{array}\right),

可以计算出 A 的两个特征值分别是 \phi=(1+\sqrt{5})/2- \phi^{-1} = (1-\sqrt{5})/2,它们所对应的特征向量分别是:

\vec{\alpha}_{1} = (\phi,1)^{T}, \vec{\alpha}_{2} = (-\phi^{-1},1).

因此直接计算可以得到

F_{k} = \frac{1}{\sqrt{5}}\bigg(\frac{1+\sqrt{5}}{2}\bigg)^{k} - \frac{1}{\sqrt{5}}\bigg(\frac{1-\sqrt{5}}{2}\bigg)^{k}=\frac{\phi^{k}-(-\phi)^{-k}}{\sqrt{5}}.

通过上面的计算方法,为了计算 Fibonacci 数列的通项公式,我们可以先把它转换成一个矩阵求幂的问题,于是我们就能矩阵对角化的方法把 Fibonacci 数列的通项公式求出来。

时间序列的弱平稳性

要讲解自回归模型,就必须提到时间序列的弱平稳性。一个时间序列 \{x_{t}\}_{t\geq 0} 具有弱平稳性(Weak Stationary)指的是:

  1. E(x_{t}) 对于所有的 t\geq 0 都是恒定的;
  2. Var(x_{t}) 对于所有的 t\geq 0 都是恒定的;
  3. x_{t}x_{t-h} 的协方差对于所有的 t\geq 0 都是恒定的。

另外,时间序列的自相关方程(AutoCorrelation Function)指的是对于 h = 1,2,3,\cdots,可以定义 ACF 为

ACF(x_{t},x_{t-h}) = \frac{Covariance(x_{t},x_{t-h})}{\sqrt{Var(x_{t})\cdot Var(x_{t-h})}}.

如果时间序列 \{x_{t}\}_{t\geq 0} 在弱平稳性的假定下,ACF 将会简化为

ACF(x_{t},x_{t-h}) = \frac{Covariance(x_{t},x_{t-h})}{Var(x_{t})}.

时间序列的自回归模型(AutoRegression Model)

AR(1) 模型

AR(1) 模型指的是时间序列 \{x_{t}\}_{t\geq 0} 在时间戳 t 时刻的取值 x_{t} 与时间戳 t - 1 时刻的取值 x_{t-1} 相关,其公式就是:

x_{t}=\delta+\phi_{1}x_{t-1}+w_{t}

这个时间序列 \{x_{t}\}_{t\geq 0} 满足如下条件:

  1. w_{t}\sim N(0,\sigma_{w}^{2}),并且 w_{t} 满足 iid 条件。其中 N(0,\sigma_{w}^{2}) 表示 Gauss 正态分布,它的均值是0,方差是 \sigma_{w}^{2}
  2. w_{t}x_{t} 是相互独立的(independent)。
  3. x_{0},x_{1},\cdots弱平稳的,i.e. 必须满足 |\phi_{1}|<1

如果选择初始条件 x_{0}=1,则可以得到一些 AR(1) 模型的例子如下图所示:

AR Models

从 AR(1) 以上的定义出发,我们可以得到:

  1. E(x_{t}) = \delta/(1-\phi_{1}).
  2. Var(x_{t}) = \sigma_{w}^{2}/(1-\phi_{1}^{2}).
  3. Covariance(x_{t},x_{t-h}) = \phi_{1}^{h}.

Proof of 1. 从 AR(1) 的模型出发,可以得到

E(x_{t}) = E(\delta + \phi_{1}x_{t-1}+w_{t})  = \delta + \phi_{1}E(x_{t-1}) = \delta + \phi_{1}E(x_{t}),

从而,E(x_{t}) = \delta/(1-\phi_{1}).

Proof of 2. 从 AR(1) 的模型出发,可以得到

Var(x_{t}) = Var(\delta + \phi_{1}x_{t-1}+w_{t})

= \phi_{1}^{2}Var(x_{t-1}) +Var(w_{t}) = \phi_{1}^{2}Var(x_{t}) + \sigma_{w}^{2},

从而,Var(x_{t}) =\sigma_{w}^{2}/(1-\phi_{1}^{2}).

Proof of 3.\mu = E(x_{t}), \forall t\geq 0. 从 x_{t} 的定义出发,可以得到:

x_{t}-\mu = \phi_{1}(x_{t-1}-\mu)+w_{t}

= \phi_{1}^{h}(x_{t-h}-\mu) + \phi_{1}^{h-1}w_{t-h+1}+\cdots+\phi_{1}w_{t-1}+w_{t},

从而,

\rho_{h} = Covariance(x_{t},x_{t-h}) = \frac{E((x_{t}-\mu)\cdot(x_{t-h}-\mu))}{Var(x_{t})}=\phi_{1}^{h}.

AR(1) 模型与一维动力系统

特别的,如果假设 w_{t} 恒等于零,就可以得到 x_{t} =\delta + \phi_{1}x_{t-1} 对于所有的 t\geq 1 都成立。也就是可以写成一个一维函数的迭代公式:

f(x) = \phi_{1}x + \delta,

下面我们要计算 f^{n}(x) 的收敛性,这里的 f^{n}(x) = f\circ\cdots\circ f(x) 表示函数 fn 次迭代。

Method 1. 

通过 f(x) 的定义直接计算可以得到:

f^{n}(x) = \phi_{1}^{n}x+ \frac{1-\phi_{1}^{n}}{1-\phi_{1}}\delta,

n\rightarrow \infty,可以得到 f^{n}(x)\rightarrow \delta/(1-\phi_{1})。这与 E(x_{t}) = \delta/(1-\phi_{1}) 其实是保持一致的。

另外,如果 |\phi_{1}|>1,可以从公式上得到 f^{n}(x) \rightarrow \inftyn\rightarrow \infty

Method 2.

将原函数转换成 Lipschitz 函数的形式,i.e. 如果 |\phi_{1}|<1,那么

f(x)-\frac{\delta}{1-\phi_{1}} = \phi_{1}(x-\frac{\delta}{1-\phi_{1}})

\implies |f(x)-\frac{\delta}{1-\phi_{1}}| <\frac{1+|\phi_{1}|}{2}\cdot|x-\frac{\delta}{1-\phi_{1}}|

\implies |f^{n}(x)-\frac{\delta}{1-\phi_{1}}|<\bigg(\frac{1+|\phi_{1}|}{2}\bigg)^{n}\cdot|x-\frac{\delta}{1-\phi_{1}}|.

因为 (1+|\phi_{1}|)/2<1,我们可以得到 \lim_{n\rightarrow\infty}f^{n}(x)=\delta/(1-\phi_{1}). i.e. f^{n}(x) 趋近于 \delta/(1-\phi_{1}).

反之,如果 |\phi_{1}|>1,很容易得到

f(x)-\frac{\delta}{1-\phi_{1}} = \phi_{1}(x-\frac{\delta}{1-\phi_{1}})

\implies |f(x)-\frac{\delta}{1-\phi_{1}}| >\frac{1+|\phi_{1}|}{2}\cdot|x-\frac{\delta}{1-\phi_{1}}|

\implies |f^{n}(x)-\frac{\delta}{1-\phi_{1}}|>\bigg(\frac{1+|\phi_{1}|}{2}\bigg)^{n}\cdot|x-\frac{\delta}{1-\phi_{1}}|.

因此,在 |\phi_{1}|>1 这种条件下,f^{n}(x)\rightarrow \infty as n\rightarrow \infty. 特别地,对于一阶差分方程 x_{t} =\delta + \phi_{1}x_{t-1} 而言,如果 |\phi_{1}|>1,那么 x_{t} 的取值会越来越大,这与现实的状况不相符,所以在时间序列的研究中,一般都会假设 |\phi_{1}|<1

AR(p) 模型

按照之前类似的定义,可以把 AR(1) 模型扩充到 AR(p) 模型,也就是说:

1. AR(1) 模型形如:

x_{t}=\delta+\phi_{1}x_{t-1}+w_{t}.

2. AR(2) 模型形如:

x_{t} = \delta + \phi_{1}x_{t-1}+\phi_{2}x_{t-2}+w_{t}.

3. AR(p) 模型形如:

x_{t} = \delta + \phi_{1}x_{t-1}+\phi_{2}x_{t-2}+\cdots+\phi_{p}x_{t-p}+w_{t}.

AR(p) 模型的稳定性 --- 基于线性代数

对于 AR(2) 模型,可以假定 \delta = 0 并且忽略误差项,因此可以得到简化版的模型形如:

x_{t}= \phi_{1}x_{t-1} + \phi_{2}x_{t-2}.

写成矩阵的形式形如:

\left( \begin{array}{c} x_{t+2}  \\ x_{t+1}  \\ \end{array}\right)= \left( \begin{array}{cc} \phi_{1} & \phi_{2} \\ 1 & 0 \\ \end{array}\right) \left( \begin{array}{c} x_{t+1}  \\ x_{t}  \\ \end{array}\right) = A \left( \begin{array}{c} x_{t+1}  \\ x_{t}  \\ \end{array}\right).

求解其特征多项式则是基于 det(\lambda I - A) = 0,求解可以得到 \lambda^{2}-\phi_{1}\lambda - \phi_{2} =0,i.e. A^{k} = P diag(\lambda_{1}^{k}, \lambda_{2}^{k})P^{-1}。当 \lambda_{1}, \lambda_{2} 都在单位圆内部的时候,也就是该模型 x_{t+2} = \phi_{1}x_{t+1}+\phi_{2}x_{t} 满足稳定性的条件。

对于更加一般的 AR(p) 模型,也就是考虑 p 阶差分方程

x_{t} = \phi_{1}x_{t-1}+\phi_{2}x_{t-2}+\cdots+\phi_{p}x_{t-p}.

可以用同样的方法将其转换成矩阵的形式,那就是:

\left(\begin{array}{c} x_{t+p} \\ x_{t+p-1} \\ \vdots \\ x_{t+1}\\ \end{array}\right) = \left(\begin{array}{ccccc} \phi_{1} & \phi_{2} &\cdots & \phi_{p-1} & \phi_{p} \\ 1 & 0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & 1 & 0 \\ \end{array}\right) \left(\begin{array}{c} x_{t+p-1} \\ x_{t+p-2} \\ \vdots \\ x_{t} \\ \end{array}\right) = A \left(\begin{array}{c} x_{t+p-1} \\ x_{t+p-2} \\ \vdots \\ x_{t} \\ \end{array}\right)

计算 det(\lambda I - A) = 0,可以得到其特征多项式为:

\lambda^{p}-\phi_{1}\lambda^{p-1}-\phi_{2}\lambda^{p-2}-\cdots-\phi_{p}=0.

当每个特征值都在单位圆盘内部的时候,i.e. |\lambda_{i}|<1, \forall 1\leq i\leq p,该 p 阶差分方程

x_{t} = \phi_{1}x_{t-1}+\phi_{2}x_{t-2}+\cdots+\phi_{p}x_{t-p}

存在稳定性的解。

 

时间序列的搜索

在之前的时间序列相似度算法中,时间戳都是一一对应的,但是在实际的场景中,时间戳有可能出现一定的偏移,但是两条时间序列却又是十分相似的。例如正弦函数 \sin(x) 和余弦函数 \cos(x),只是平移了 \pi/2 个长度而已。本文将会介绍一些基于形状的时间序列的距离算法,并且介绍如何在给定时间序列的情况下,在时间序列数据库中寻找相似的时间序列。

 

基于动态规划的相似度计算方法

DTW 的经典算法

基于动态规划的相似度计算的典型代表就是 DTW(Dynamical Time Warping )和 Frechet 距离。在这里将会主要介绍 DTW 算法。详细来说,DTW 算法是为了计算语音信号处理中,由于两个人说话的时长不一样,但是却是类似的一段话,欧几里德算法不完全能够解决这类问题。在这种情况下,DTW 算法就被发展出来。DTW 算法是为了计算两条时间序列的最佳匹配点,假设我们有两条时间序列 Q 和 C,长度都是 n,并且 Q:\{q_{1},q_{2},\cdots,q_{n}\}C:\{c_{1},c_{2},\cdots,c_{n}\}。首先我们可以建立一个 n\times n 的矩阵,(i,j) 位置的元素是 dist(q_{i},c_{j}),这里的 dist 可以使用 L^{1}, L^{2}, L^{p}, L^{\infty} 范数。其次,我们想找到一条路径,使得这个矩阵的累积距离最小,而这条路则是两条时间序列之间的最佳匹配。在这里,我们可以假设这条路径是 W: \{w_{1},\cdots,w_{K}\},其中 W 的每个元素表示时间序列 Q 中的第 i 个元素和时间序列 C 中的第 j 个元素之间的距离. i.e. w_{k}=(q_{i}-c_{j})^{2}

DTW1

现在我们需要找到一条路径使得

W^{*} = argmin_{W}(\sqrt{\sum_{k=1}^{K}w_{k}}).

这条路径就是动态规划的解,它满足一个动态规划方程:对于 1\leq i\leq n, 1\leq j\leq n,有

DTW(i,j)

= dist(q_{i},c_{j}) + \min(DTW(i-1,j-1), DTW(i-1,j), DTW(i,j-1)).

其初始状态是 DTW(0,0)=0, DTW(i,0)=\infty, DTW(0,j)=\infty, \forall 1\leq i\leq n, 1\leq j\leq n. 最终的取值 DTW(n,n) 就是我们需要的解,也就是两条时间序列的 DTW 距离。按照上面的算法,DTW 算法的时间复杂度是 O(n^{2})。特别地,

  1. 如果 dist(q_{i},c_{j}) = (q_{i}-c_{j})^{2} 时,则 \sqrt{DTW(n,n)} 表示最后的距离;
  2. 如果 dist(q_{i},c_{j}) = |q_{i}-c_{j}| 时,则 DTW(n,n) 表示最后的距离。
  3. 如果 dist(q_{i},c_{j}) = |q_{i}-c_{j}|^{p} 时,则 (DTW(n,n))^{1/p} 表示最后的距离。

Remark. 

DTW 算法不满足三角不等式。例如:x={0}, y={1,2}, z={1,2,2},则

DTW(x,z)=5>DTW(x,y)+DTW(y,z) = 3 + 0 =3.

DTW 的加速算法

某些时候,我们可以添加一个窗口长度的限制,换言之,如果要比较 q[i]c[j] 的话,i 与 j 需要满足 |i-j|\leq w,这里的 w 表示窗口长度。因此算法的描述如下:

初始条件和之前一样。

DTW(i,j) = dist(q_{i},c_{j}) + \min(DTW(i-1,j-1), DTW(i-1,j), DTW(i,j-1)),

这里的 j 取值范围是:对每一个 1\leq i\leq n,需要 j 满足\max(1,i-w) \leq j\leq \min(m,i+w)

 

相似时间序列的搜索

相似的时间序列的搜索问题一般是这样的:

Question 1. 给定一条时间序列 Q 和一个时间序列的数据库 D=\{C_{i}\}_{i=1}^{\infty}。通过某种相似度或者距离计算方法,计算出给定的时间序列 Q 与时间序列数据库中 D 中最相似的时间序列。

Question 2. 给定一条时间序列 Q 和一个时间序列的数据库 D=\{C_{i}\}_{i=1}^{\infty},以及正整数 K。从数据库 D 中寻找与给定的时间序列 Q 最相似的 K 条时间序列。

DTW 算法的下界 LB_Kim

对于两条时间序列 Q 和 C 而言,可以分别提取它们的四个特征,那就是最大值,最小值,第一个元素的值,最后一个元素的值。这样可以计算出 LB_Kim

LBKim(Q,C)

= \max\{|\max(Q)-\max(C)|,|\min(Q)-\min(C)|,|First(Q)-First(C)|,|Last(Q)-Last(C)| \}.

可以证明 LBKim(Q,C)\leq DTW(Q,C).

DTW 算法的下界 LB_Yi

有学者在 LB_Kim 的基础上提出了一种下界的计算方法,那就是 LB_Yi,有兴趣的读者可以参见:efficient retrieval of similar time sequences under time warping, 1998.

DTW 算法的下界 LB_Keogh

但是,如果是基于 DTW 的距离计算方法,每次都要计算两条时间序列的 DTW 距离,时间复杂度是 O(n^{2})。于是就有学者研究是否存在 DTW 距离的下界表示,也就是说找到一个合适的下界,Lower Bound of DTW。每次判断 Lower Bound of DTW 是否小于当前的最小距离,如果下界高于最小距离,就不需要进行 DTW 的计算;否则开始计算 DTW 的值。如果下界的计算速度足够快,并且下界足够精准的话,可以大量的压缩搜索的时间。于是,Keogh 提出了下界的计算方法。

(1)首先定义时间序列 Q 的上下界。i.e. Q:\{q_{1},\cdots,q_{n}\},给定一个窗口的取值 r,得到 U_{i}=\max(q_{i-r},q_{i+r})L_{i} = \min(q_{i-r},q_{i+r})

(2)定义公式:

LBKeogh(Q,C)

= \sqrt{\sum_{i=1}^{n}(c_{i}-U_{i})^{2}I_{\{c_{i}>U_{i}\}} + \sum_{i=1}^{n}(c_{i}-L_{i})^{2}I_{\{c_{i}<L_{i}\}}}.

其中,LBKeogh 满足性质:

对于任意两条长度为 n 的时间序列 Q 和 C,对于任意的窗口 r\geq 1,有不等式 LBKeogh(Q,C)\leq DTW(Q,C) 成立。

所以,可以把搜索的算法进行加速:

Lower Bounding Sequential Scan(Q)
    best_so_far = infinity 
    for all sequences in database
        LB_dist = lower_bound_distance(C_{i},Q)
        if LB_dist < best_so_far
            true_dist = DTW(C_{i},Q)
            if true_dist < best_so_far
                best_so_far = true_dist
                index_of_best_match = i
            endif
        endif
    endfor

在论文 Exact Indexing of Dynamic Time Warping 里面,作者还尝试将 Piecewise Constant Approximation 与 LB_Keogh 结合起来,提出了 LB_PAA 的下界。它满足条件 LBPAA(Q,C)\leq LBKeogh(Q,C)\leq DTW(Q,C).

 

总结

本文初步介绍了 DTW 算法以及它的下界算法,包括 LB_Kim, LB_Keogh 等,以及时间序列数据库的搜索算法。