【论文解读】【多元时间序列异常检测】 Detection and Characterization of Anomalies in Multivariate Time Series

这是一篇比较经典的多元时间序列数据异常检测算法的论文阅读解析，包括了算法代码的实现，论文原文已经上传到CSDN方便下载阅读，点击这里下载原始论文

一、这篇论文说了什么？

针对多元时间序列的异常数据检测，提出了一种图遍历算法，作者认为，多元数据序列数据都可以看成n个数据点，然后每个数据点上对应了多个变量，比如每天的天气数据，有以下变量：温度、光照、湿度、空气质量等。这是常见的多元时间序列数据。那么每个数据点都有一个温度、光照、湿度和空气质量4个变量，那么就认为整个数据都可以构成一个图，图有节点和连接每两个节点的边，每条边都有自己的权值，我们认为这是一个带权图，作者认为，n个数据点，可以构成一个n x n的图，每两个节点都进行连接，那么每两个节点之间连接的边权就认为是这两个数据节点的相似度数值，于是我们得到了一个n x n大小的相似矩阵，第(i,j)个数值就代表了节点i和节点j的数据相似度，于是我们把这个矩阵和马尔科夫链结合起来，如果对马尔科夫链了解不多的可以先去了解下，这里简单说一下，马尔科夫链是一种概率估计预测算法，通过已经存在的原始数据进行概率计算，然后得到一个概率矩阵，然后利用这个矩阵进行迭代运算，如果一开始有x个状态，那么最后也会有x个状态对应的出现的概率，在第step次就可以知道状态xi发生的概率，那么对应到我们的论文中来，作者巧妙结合了马尔科夫链的特性，将我们上述得到的相似矩阵进行这样的迭代运算，最后得到n个数值，也就是对应的n个数据节点的发生概率，如果某个数据节点发生概率很低，我们认识这个数据节点和其他的n-1个数据节点是不太相似的，然后认为这个数据节点是异常。这个比较好理解。

上述方法提出来后就已经可以去进行异常检测了，但是检测效果不是很理想，因为这样的计算方法太过笼统，没法专注于某个具体的变量的异常中（比如温度或空气质量），每个变量可以单独作为一条时间序列数据，那么也可以得到单独的一个相似矩阵K，总共P个变量，就可以得到P个相似矩阵，我们通过某些方法，让这个相似矩阵可以专注于某个变量的异常检测，于是提出了核矩阵对齐的思想。

作者的想法是，把P个变量（P条时间序列）划分为P-1个变量X，把我们要检测的目标变量为Y，然后得到了两个矩阵Kx和Ky，然后再提出一个转换的相似矩阵Ka，Ka矩阵是从Kx过渡到ky的，也就是说ka应该在包含Kx的基本特征信息的基础上，同时偏向于Ky的特征信息，为啥要这样处理？可以这样理解，一开始非常笼统的计算得到的相似矩阵K没有偏向性，因为变量之间本身就是相互影响的，这样直接计算得到的相似矩阵K可能表达效果就不是很好，通过Kx和Ky的这样转变方法，我们可以理解为在拥有了其他的P-1个变量信息后，再向目标变量进行一定偏向，那么整体得到的相似矩阵Ka就可以很好的在保证了接受所有变量的数据信息的情况下把我们的目标变量的异常信息表达出来。这里也叫特定目标变量检测。

那么通过上面的方法，我们得到了很具体的针对某个变量i的目标异常检测方法，那么现在怎么去计算整体的表现呢？这里作者提出的方法很简单粗暴，因为Kai可以很好的表达特定的目标变量i的异常信息情况，那么针对整体的数据，也就是针对所有的P个变量，我们可以得到P个特定的目标相似矩阵Ka1,Ka2,Ka3,…,Ka(p-1),Ka§，
那么此时整体的相似矩阵就直接为Ka1xKa2xKa3…xKa(p-1)xKa§，注意这里是矩阵乘法相乘。

好了，说到这里整个论文的内容就已经说完了，现在梳理下我们要干什么？计算相似矩阵，然后根据作者提出的迭代公式进行迭代计算，最后得出的数值最小的数据点就是异常。

二、具体公式

2.1迭代公式

作者根据马尔科夫链提出的迭代公式

其中d是一个参数，需要自己调整，n为节点个数，也就是矩阵S的大小，S为相似矩阵K进行归一化处理（把每个数据映射在0到1之间）得到的结果，c为一个nx1的一维向量，表示每个节点和其他n-1个节点总的一个连接的相似度权重，观察上面的公式，S为nxn，c为nx1，那么Sxc得到的还是一个nx1的一维向量，于是可以这样反复迭代运算下去，那么S我们知道是相似矩阵K归一化得到的，相似矩阵K是我们接下来要计算的，那么c初始化赋值怎么弄呢？这里作者在文中没提到，所以这里我认为c可以初始化赋值为1，因为表达每个节点之间相似关系的信息都在S里，经过不断迭代我们可以知道最后的关系，所以c初始化为1没有问题。

2.2相似矩阵计算

一开始的相似矩阵K计算公式如下

其中p为变量的数目，如果p=1就是单变量时间序列异常检测了，上面这个计算函数叫RFB径向基函数，xki和xkj为对应的不同时间节点上的数据点。

这个是针对单个数据节点，经过变换也可以去计算一段连续的时间片段上的数据点，如下

其中d为时间片段的宽度。

当然上面这个方法计算得到的相似矩阵效果不是很好，然后作者又想进行核矩阵对齐，上面我们说了，将P个变量划分为P-1个变量X和目标变量Y，分别得到相似矩阵Kx和Ky，然后再求中间的一个转换的相似矩阵Ka可以包括Kx的信息的同时偏向于Ky的表达，也就是Ka和Ky的相似度要尽可能大，于是提出计算公式如下：

其中

上面这个公式不难理解辣，就是让矩阵Ka和Ky按上面的计算方式进行计算，最后会得到一个数值，这个数值我们叫相似值，然后希望这个数值尽可能大，我们根据变量Y的时间序列数据可以算出相似矩阵Ky，但是现在我们不知道Ka是多少，然后再进行分析，因为Ka是Kx演变过来的，所以我们认为可以对Kx进行特征变量拆分，把Kx拆分成一些基本的信息变量vi，因为x是由P-1个变量，注意，论文里的P是针对了X的变量数目，下面我描述的时候为了符合论文内容以P表示变量内容，读者注意区分。因为Kx是由集合X里面的P个变量计算得到的，为了拆分作者认为可以拆分成P个基础信息向量v1,v2,v3,…,v(p-1),vp，那么转换为Ka就需要再乘一个转换的参数ai，也是P个为a1,a2,a3,…,a(p-1),ap。
于是Ka和Kx的基础信息向量vi的转换公式为：

上面这个公式，vi是一个nx1的向量，那么vi的转置vi’是一个1xn的向量，那么vixvi’得到的是一个nxn的向量，这个时候再乘转换参数ai，再把这P个基础的矩阵信息加起来，得到我们最后的转换的相似矩阵Ka。这里应该比较好理解的。

再往下是计算vi，因为vi是矩阵Kx的基础信息向量，怎么把kx分解为P个vi呢，作者提出了两种分解方式。

2.2.1求解特征方程

求解特征方程比较好理解了吧，线性代数的知识，任何矩阵都可以通过求解其特征方程计算得到该矩阵对应的特征值和特征向量，计算得到特征值和特征向量后，我们的公式就变换如下：

注意上面特征向量乘特征向量的转置为单位向量，所以得到了上面的变换过程。现在我们知道了基础信息向量vi为求解特征方程后的特征向量，现在的目标是求解ai，将上述的公式再进行转换得到下面的

然后这种情况下得到了答案如下

但是作者又提出，认为这种方式太过于贴近了Ky而忽略了Kx原本的一些特征信息，也就是对Ky而言过拟合了，于是作者又将公式进行一定的转换，用来避免过拟合，如下

通过引入Kx的特征值来修正ai，这样求解得到的答案为

至此，通过求解特征方程来计算得到ai从而得到Ka的方法结束。

2.2.2替换法

这种方法比较简单粗暴，直接用原始变量Xi替换掉vi用来作为特征向量计算，作者后面的实验也用的第二种方法，这种方法在实现上也比较容易。

再通过求解得到答案为

首先看左边，Xi和Ky都是已经知道的，左边计算得到一个数值，再看右边Xi和Xj都是知道的，右边会得到一个含有未知数aj的方程，整体构成一个齐次线性方程组。
通过求解齐次线性方程组便可以算出具体的每个aj，然后结束。

三、复现的代码
复现的代码不一定完全正确，只是供参考学习使用。
因为写的代码使用了相应的数据集，这里我直接把代码和数据集一起打包上传到Github。
点击这里下载复现代码

希望对你的学习有所帮助，如果有问题请及时指出，谢谢~