Category Archives: 智能运维

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

随着 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

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

成本

机器学习+其他领域

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

 

Advertisements

黑盒函数的探索

黑盒函数的定义

在工程上和实际场景中,黑盒函数(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 等,以及时间序列数据库的搜索算法。

 

如何理解时间序列?— 从Riemann积分和Lebesgue积分谈起

Riemann 积分和 Lebesgue 积分是数学中两个非常重要的概念。本文将会从 Riemann 积分和 Lebesgue 积分的定义出发,介绍它们各自的性质和联系。

积分

Riemann 积分

Riemann 积分虽然被称为 Riemann 积分,但是在 Riemann 之前就有学者对这类积分进行了详细的研究。早在阿基米德时代,阿基米德为了计算曲线 x^{2} 在 [0,1] 区间上与 X 坐标轴所夹的图形面积,就使用了 Riemann 积分的思想。 他把 [0,1] 区间等长地切割成 n 段,每一段使用一个长方形去逼近 x^{2} 这条曲线的分段面积,再把 n 取得很大,所以得到当 n 趋近于无穷的时候,就知道该面积其实是 1/3。

下面来看一下 Riemann 积分的详细定义。

Riemann Integral1

考虑定义在闭区间 [a,b] 上的函数 f(x)

取一个有限的点列 a=x_{0}<x_{1}<\cdots<x_{n}=b\lambda = \max(x_{i+1}-x_{i}) 表示这些区间长度的最大值,在这里 0\leq i\leq n-1。在每一个子区间上[x_{i},x_{i+1}] 上取出一个点 t_{i}\in[x_{i},x_{i+1}]。而函数 f(x) 关于以上取样分割的 Riemann 和就是以下公式:

\sum_{i=0}^{n-1}f(t_{i})(x_{i+1}-x_{i}).

当我们说该函数 f(x) 在闭区间 [a,b] 上的取值是 s 的意思是:

对于任意 \epsilon>0,存在 \delta>0 使得对于任意取样分割,当 \lambda<\delta 时,就有

|\sum_{i=0}^{n-1}f(t_{i})(x_{i+1}-x_{i}) - s|<\epsilon.

通常来说,用符号来表示就是:\int_{a}^{b}f(x)=s.

用几幅图来描述 Riemann 积分的思想就是:

 

Lebesgue 积分

Riemann 积分是为了计算曲线与 X 轴所围成的面积,而 Lebesgue 积分也是做同样的事情,但是计算面积的方法略有不同。要想直观的解释两种积分的原理,可以参见下图:

RiemannANDLebesgue
Riemann 积分(上)与 Lebesgue 积分(下)

Riemann 积分是把一条曲线的底部分成等长的区间,测量每一个区间上的曲线高度,所以总面积就是这些区间与高度所围成的面积和。

Lebesgue 积分是把曲线化成等高线图,每两根相邻等高线的差值是一样的。每根等高线之内含有它所圈着的长度,因此总面积就是这些等高线内的面积之和。

用再形象一点的语言来描述就是:吃一块汉堡有多种方式

  1. Riemann 积分:从一个角落开始一口一口吃,每口都包含所有的配料;
  2. Lebesgue 积分:从最上层开始吃,按照“面包-配菜-肉-蛋-面包”的节奏,一层一层来吃。

再看一幅图的表示就是:Riemann 积分是按照蓝色的数字顺序相加的,Lebesgue 积分是按照红色的数字来顺序相加的。

Riemann and Lebesgue1

 

基于这些基本的思想,就可以给出 Lebesgue 积分的定义:

简单函数指的是对指示函数的有限线性组合,i.e. \sum_{k}a_{k}I_{S_{k}},这里的 a_{k} 是系数,S_{k} 是可测集合,I 表示指示函数。当 a_{k} 非负时,令

\int(\sum_{k}a_{k}1_{S_{k}})d\mu = \sum_{k}a_{k}\int 1_{S_{k}}d\mu = \sum_{k}a_{k}\mu(S_{k}).

如果 f 是一个非负可测函数时,可以定义函数 f 在可测集合 E 上的 Lebesgue 积分是:

\int_{E}f d\mu = \sup\{\int_{E}sd\mu: \bold{0}\leq s\leq f\},

这里的 s 指的是非负简单函数,\bold{0} 表示零函数,这里的大小关系表示对定义域内的每个点都要成立。

而对于可测函数时,可以把可测函数 f 转换成 f= f^{+}-f^{-},而这里的 f^{+}f^{-} 都是非负可测函数。所以可以定义任意可测函数的 Lebesgue 积分如下:

\int fd\mu = \int f^{+}d\mu - \int f^{-}d\mu..

Riemann 积分与Lebesgue 积分的关系

定义了两种积分之后,也许有人会问它们之间是否存在矛盾?其实,它们之间是不矛盾的,因为有学者证明了这样的定理:

如果有界函数 f 在闭区间 [a,b] 是 Riemann 可积的,则它也是 Lebesgue 可积的,并且它们的积分值相等:

(R)\int_{a}^{b}f(x)dx = (L)\int_{[a,b]}f(x)dx.

左侧是表示 Riemann 积分,右侧表示 Lebesgue 积分。

形象化一点的语言描述就是:无论从角落一口一口地吃汉堡,还是从顶至下一层一层吃,所吃的汉堡都是同一个。

但是 Lebesgue 积分比 Riemann 积分有着更大的优势,例如 Dirichlet 函数,

  1. x 是有理数时,D(x) = 1
  2. x 是无理数时,D(x) = 0.

Dirichlet 函数是定义在实数轴的函数,并且值域是 \{0,1\},无法画出函数图像,它不是 Riemann 可积的,但是它 Lebesgue 可积。

 

时间序列

提到时间序列,也就是把以上所讨论的连续函数换成离散函数而已,把定义域从一个闭区间 [a,b] 换成 \{1,2,3,\cdots\} 这样的定义域而已。所以,之前所讨论的很多连续函数的想法都可以应用在时间序列上。

时间序列的表示 — 基于 Riemann 积分

现在我们可以按照 Riemann 积分的计算方法来表示一个时间序列的特征,于是就有学者把时间序列按照横轴切分成很多段,每一段使用某个简单函数(线性函数等)来表示,于是就有了以下的方法:

  1. 分段线性逼近(Piecewise Linear Approximation)
  2. 分段聚合逼近(Piecewise Aggregate Approximation)
  3. 分段常数逼近(Piecewise Constant Approximation)

说到这几种算法,其实最本质的思想就是进行数据降维的工作,用少数的数据来进行原始时间序列的表示(Representation)。用数学化的语言来描述时间序列的数据降维(Data Reduction)就是:把原始的时间序列 \{x_{1},\cdots,x_{N}\}\{x_{1}^{'},\cdots, x_{D}^{'}\} 来表示,其中 D<N。那么后者就是原始序列的一种表示(representation)。

分段聚合逼近(Piecewise Aggregate Approximation)— 类似 Riemann 积分

在这种算法中,分段聚合逼近(Piecewise Aggregate Approximation)是一种非常经典的算法。假设原始的时间序列是 C = \{x_{1},\cdots, x_{N}\},定义 PAA 的序列是:\overline{C} = \{\overline{x}_{1},\cdots,\overline{x}_{w}\}

其中

\overline{x}_{i} = \frac{w}{N} \cdot \sum_{j=\frac{N}{w}(i-1)+1}^{\frac{N}{w}i} x_{j}.

在这里 1\leq i\leq w。用图像来表示那就是:

PAA

至于分段线性逼近(Piecewise Linear Approximation)和分段常数逼近(Piecewise Constant Approximation),只需要在 \overline{x}_{i} 的定义上稍作修改即可。

符号特征(Symbolic Approximation)— 类似用简单函数来计算 Lebesgue 积分

在推荐系统的特征工程里面,特征通常来说可以做归一化,二值化,离散化等操作。例如,用户的年龄特征,一般不会直接使用具体的年月日,而是划分为某个区间段,例如 0~6(婴幼儿时期),7~12(小学),13~17(中学),18~22(大学)等阶段。

其实在得到分段特征之后,分段特征在某种程度上来说依旧是某些连续值,能否把连续值划分为一些离散的值呢?于是就有学者使用一些符号来表示时间序列的关键特征,也就是所谓的符号表示法(Symbolic Representation)。下面来介绍经典的 SAX Representation。

如果我们希望使用 \alpha 个符号来表示时间序列,那么我们其实可以考虑正态分布 N(0,1),用\{z_{1/\alpha},\cdots,z_{(\alpha-1)/\alpha}\} 来表示 Gauss 曲线下方的一些点,而这些点把 Gauss 曲线下方的面积等分成了 \alpha 段。用 \{l_{1},\cdots,l_{\alpha}\} 表示 \alpha 个字母。

SAX 方法的流程如下:

  1. 正规化(normalization):把原始的时间序列映射到一个新的时间序列,新的时间序列满足均值为零,方差为一的条件。
  2. 分段表示(PAA):\{x_{1},\cdots, x_{N}\} \Rightarrow  \{\overline{x}_{1},\cdots,\overline{x}_{w}\}
  3. 符号表示(SAX):如果 \overline{x}_{i}<z_{1/\alpha},那么 \hat{X}_{i}=l_{1};如果 z_{(j-1)/\alpha}\leq \overline{x}_{i}<z_{j/\alpha},那么 \hat{X}_{i} = l_{j},在这里 2\leq j\leq \alpha;如果 \overline{x}_{i}\geq z_{(\alpha-1)/\alpha},那么 \hat{X}_{i} = l_{\alpha}

于是,我们就可以用 \{l_{1},\cdots,l_{\alpha}\}\alpha 个字母来表示原始的时间序列了。

SAX

时间序列的表示 — 基于 Lebesgue 积分

要想考虑一个时间序列的值分布情况,其实就类似于 Lebesgue 积分的计算方法,考虑它们的分布情况,然后使用某些函数去逼近时间序列。要考虑时间序列的值分布情况,可以考虑熵的概念。

熵(Entropy)

通常来说,要想描述一种确定性与不确定性,熵(entropy)是一种不错的指标。对于离散空间而言,一个系统的熵(entropy)可以这样来表示:

\text{entropy}(X) = -\sum_{i=1}^{\infty}P\{x=x_{i}\}\ln(P\{x=x_{i}\}).

如果一个系统的熵(entropy)越大,说明这个系统就越混乱;如果一个系统的熵越小,那么说明这个系统就更加确定。

提到时间序列的熵特征,一般来说有几个经典的熵指标,其中有一个就是 binned entropy。

分桶熵(Binned Entropy)

从熵的定义出发,可以考虑把时间序列的值进行分桶的操作,例如,可以把 [min, max] 这个区间等分为十个小区间,那么时间序列的取值就会分散在这十个桶中。根据这个等距分桶的情况,就可以计算出这个概率分布的熵(entropy)。i.e. Binned Entropy 就可以定义为:

\text{binned entropy}(X) = -\sum_{k=0}^{\min(maxbin, len(X))} p_{k}\ln(p_{k})\cdot 1_{(p_{k}>0)},

其中 p_{k} 表示时间序列 X 的取值落在第 k 个桶的比例(概率),maxbin 表示桶的个数,len(X) 表示时间序列 X 的长度。

如果一个时间序列的 Binned Entropy 较大,说明这一段时间序列的取值是较为均匀的分布在 [min, max] 之间的;如果一个时间序列的 Binned Entropy 较小,说明这一段时间序列的取值是集中在某一段上的。

总结

在本篇文章中,笔者从 Riemann 积分和 Lebesgue 积分出发,介绍了它们的基本概念,性质和联系。然后从两种积分出发,探讨了时间序列的分段特征,时间序列的熵特征。在未来的 Blog 中,笔者将会介绍时间序列的更多相关内容。

 

时间序列的相似性

在文本挖掘中,可以通过 Word2Vec 生成的向量以及向量的内积,或者根据语义和词性来判断两个词语是否是近义词。在时间序列的挖掘中,同样可以找到一些方法来描述两条时间序列是否相似。

在介绍时间序列的距离之前,笔者感觉需要回顾一下数学中度量空间内积空间的定义。

度量空间

在数学里面,集合 M 上的距离函数定义为 d: M\times M\rightarrow \mathbb{R},其中 \mathbb{R} 表示实数集合,并且函数 d 满足以下几个条件:

  1. d(x,y)\geq 0,并且 d(x,y)=0 当且仅当 x=y;
  2. d(x,y)=d(y,x),也就是满足对称性;
  3. d(x,z)\leq d(x,y)+d(y,z),也就是三角不等式。

满足这三个条件的距离函数称为度量,具有某种度量的集合则叫做度量空间

Remark.

在本文下面的时间序列距离的各种定义中,这些距离不一定满足三角不等式。例如动态时间算法(DTW)就不满足三角不等式。

内积空间

一个内积空间是域 F(其中 F=\mathbb{R} 或者 F=\mathbb{C})上的向量空间 V 与一个内积(映射)所构成,V 上的一个内积定义为正定,非退化的共轭双线性形式,记成 <\cdot,\cdot>:V\times V\rightarrow F,它满足以下设定:

1. 对于任意的 x,y\in V,有 <x,y> =\overline{<y,x>}.

2. 共轭双线性形式指的是:

\forall a\in F, \forall x,y\in V, <ax,y>=a<x,y>,

\forall x,y,z\in F, <x+y,z> = <x,z> + <y,z>.

\forall b\in F, \forall x,y\in V, <x,by> = \overline{b}<x,y>,

\forall x,y,z\in F,<x,y+z> = <x,y>+<x,z>.

3. 非负性:\forall x\in V, <x,x>\geq 0.

4. 非退化:从 V 到对偶空间 V^{*} 的映射:x\mapsto<x,\cdot> 是同构映射。

 

基于欧几里德距离的相似度计算

假设两条时间序列曲线为 X_{T} = \{x_{1},\cdots,x_{T}\}Y_{T} = \{y_{1},\cdots,y_{T}\},于是可以使用欧几里德空间里面的 L^{1}, L^{2}, L^{p}, L^{\infty} 范数来表示两个时间序列之间的距离。用公式来描述就是:

d_{L^{1}}(X_{T},Y_{T}) = \sum_{t=1}^{T}|x_{t}-y_{t}|,

d_{L^{p}}(X_{T}, Y_{T}) = (\sum_{t=1}^{T}|x_{t}-y_{t}|^{p})^{1/p},

d_{L^{2}}(X_{T}, Y_{T}) = (\sum_{t=1}^{T}|x_{t}-y_{t}|^{2})^{1/2},

d_{L^{\infty}}(X_{T},Y_{T}) = \max_{1\leq t\leq T}|x_{t}-y_{t}|.

 

基于相关性的相似度计算方法

Pearson 系数(Pearson Coefficient)

\text{COR}(X_{T},Y_{T}) = \frac{\sum_{t=1}^{T}(x_{t}-\overline{X}_{T})\cdot(y_{t}-\overline{Y}_{T})}{\sqrt{\sum_{t=1}^{T}(x_{t}-\overline{X}_{T})^{2}}\cdot\sqrt{\sum_{t=1}^{T}(y_{t}-\overline{Y}_{T})^{2}}},

其中,

\overline{X}_{T} = \sum_{t=1}^{T}x_{t}/T, \overline{Y}_{T} = \sum_{t=1}^{T}y_{t}/T

Pearson 系数的性质如下:

  1. 如果两条时间序列 X_{T} = Y_{T},则 \text{COR}(X_{T},Y_{T}) =1 表是它们是完全一致的,如果两条时间序列 X_{T} = -Y_{T},则 \text{COR}(X_{T},Y_{T}) = -1 表示它们之间是负相关的。
  2. -1\leq \text{COR}(X_{T},Y_{T})\leq 1.

可以基于 Pearson 系数来制定两条时间序列之间的距离:

d_{COR,1}(X_{T},Y_{T}) = \sqrt{2\cdot(1-COR(X_{T},Y_{T}))},

d_{COR,2}(X_{T},Y_{T}) = \sqrt{\big(\frac{1-COR(X_{T},Y_{T})}{1+COR(X_{T},Y_{T})}\big)^{\beta}},

其中\beta\geq 0.

The First Order Temporal Correlation Coefficient

这个相关性系数与 Pearson 系数类似,但是略有不同,其定义为:

\text{CORT}(X_{T},Y_{T}) = \frac{\sum_{t=1}^{T-1}(x_{t+1}-x_{t})\cdot(y_{t+1}-y_{t})}{\sqrt{\sum_{t=1}^{T-1}(x_{t+1}-x_{t})^{2}}\cdot\sqrt{\sum_{t=1}^{T-1}(y_{t+1}-y_{t})^{2}}},

\text{CORT}(X_{T},Y_{T}) 的性质:

  1. -1\leq \text{CORT}(X_{T},Y_{T}) \leq 1
  2. \text{CORT}(X_{T},Y_{T}) =1 表示两条时间序列持有类似的趋势, 它们会同时上涨或者下跌,并且涨幅或者跌幅也是类似的。
  3. \text{CORT}(X_{T},Y_{T})=-1 表示两条时间序列的上涨和下跌趋势恰好相反。
  4. \text{CORT}(X_{T},Y_{T})=0 表示两条时间序列在单调性方面没有相关性。

基于 CORT,同样可以定义时间序列的距离,用公式描述如下:

d_{CORT}(X_{T},Y_{T}) = \phi_{k}[CORT(X_{T},Y_{T})]\cdot d(X_{T},Y_{T}),

其中,d(X_{T},Y_{T}) 可以用 d_{L^{1}}, d_{L^{p}}, d_{L^{2}}, d_{L^{\infty}}, d_{DTW}, d_{Frechet} 来计算,而

\phi_{k}(u) = 2/(1+\exp(ku)), k\geq 0

是一个递减函数。

基于自相关系数的距离(Autocorrelation-based distance)

假设时间序列是 X_{T} = \{x_{1},\cdots,x_{T}\},对于任意的 k<T,可以定义自相关系数为:

\hat{\rho}_{k} = \frac{1}{(T-k)\sigma^{2}}\sum_{t=1}^{T-k}(x_{t}-\mu)\cdot(x_{t+k}-\mu),

其中 \mu, \sigma^{2} 分别表示该时间序列的均值和方差。该公式相当于是比较整个时间序列 X_{T}=\{x_{1},\cdots,x_{T}\} 的两个子序列的相似度(Pearson 系数),这两个子序列分别是 \{x_{1},\cdots,x_{T-k}\}\{x_{k+1},\cdots,x_{T}\}

于是,通过给定一个正整数 L<T,可以对每一个时间序列得到一组自相关系数的向量,用公式描述如下:

\hat{\rho}_{X_{T}} = (\hat{\rho}_{1,X_{T}},\cdots,\hat{\rho}_{L,X_{T}})^{T}\in \mathbb{R}^{L},

\hat{\rho}_{Y_{T}} = (\hat{\rho}_{1,Y_{T}},\cdots,\hat{\rho}_{L,Y_{T}})^{T}\in\mathbb{R}^{L}.

对于 i>L 的情况,可以假定 \hat{\rho}_{i,X_{T}} = 0\hat{\rho}_{i,Y_{T}} = 0。于是,可以定义时间序列之间的距离如下:

d_{ACF}(X_{T},Y_{T}) = \sqrt{(\hat{\rho}_{X_{T}}-\hat{\rho}_{Y_{T}})^{T}\Omega(\hat{\rho}_{X_{T}}-\hat{\rho}_{Y_{T}})}.

其中的 \Omega 表示一个 L\times L 的矩阵。它有着很多种选择,例如:

(1)\Omega = I_{L} 表示单位矩阵。用公式表示就是

d_{ACFU}(X_{T},Y_{T}) =\sqrt{\sum_{i=1}^{L}(\hat{\rho}_{i,X_{T}}-\hat{\rho}_{i,Y_{T}})^{2}}.

(2)\Omega = diag\{p(1-p),p(1-p)^{2},\cdots,p(1-p)^{L}\} 表示一个 L\times L 的对角矩阵,其中 0<p<1。此时相当于一个带权重的求和公式。

d_{ACFU}(X_{T},Y_{T}) =\sqrt{\sum_{i=1}^{L}p(1-p)^{i}(\hat{\rho}_{i,X_{T}}-\hat{\rho}_{i,Y_{T}})^{2}}.

除了自相关系数(Autocorrelation Coefficients)之外,也可以考虑偏自相关系数(Partial Autocorrelation Coefficients),使用 PACFs 来取代 ACFs。这样,使用同样的定义方式就可以得到 d_{PACFU}d_{PACFG} 两个距离公式。

 

基于周期性的相似度计算方法

这里会介绍基于周期图表(Periodogram-based)的距离计算方法。其大体思想就是通过 Fourier 变换得到一组参数,然后通过这组参数来反映原始的两个时间序列时间的距离。用数学公式来描述就是:

I_{X_{T}}(\lambda_{k}) = T^{-1}|\sum_{t=1}^{T}x_{t}e^{-i\lambda_{k}t}|^{2},

I_{Y_{T}}(\lambda_{k}) = T^{-1}|\sum_{t=1}^{T}y_{t}e^{-i\lambda_{k}t}|^{2}.

其中 \lambda_{k} = 2\pi k/Tk = 1,\cdots,nn=[(T-1)/2]。这里的 [\cdot] 表示 Gauss 取整函数。

(1)用原始的特征来表示距离:

d_{P}(X_{T},Y_{T}) = \frac{1}{n}\sqrt{\sum_{k=1}^{n}(I_{X_{T}}(\lambda_{k})-I_{Y_{T}}(\lambda_{k}))^{2}}.

(2)用正则化之后的特征来描述就是:

d_{P}(X_{T},Y_{T}) = \frac{1}{n}\sqrt{\sum_{k=1}^{n}(NI_{X_{T}}(\lambda_{k})-NI_{Y_{T}}(\lambda_{k}))^{2}},

其中 NI_{X_{T}}(\lambda_{k})=I_{X_{T}}(\lambda_{k})/\hat{\gamma}_{0,X_{T}}NI_{Y_{T}}(\lambda_{k})=I_{Y_{T}}(\lambda_{k})/\hat{\gamma}_{0,Y_{T}}\hat{\gamma}_{0,X_{T}}\hat{\gamma}_{0,Y_{T}} 表示 X_{T}, Y_{T} 的标准差(sample variance)。

(3)用取对数之后的特征表示:

d_{LNP}(X_{T},Y_{T}) = \frac{1}{n}\sqrt{\sum_{k=1}^{n}(\ln NI_{X_{T}}(\lambda_{k})-\ln NI_{Y_{T}}(\lambda_{k}))^{2}}.

 

 

基于模型的相似度计算

Piccolo 距离

基于模型的相似度判断本质上是用一个模型和相应的一组参数去拟合某条时间序列,然后得到最优的一组参数,计算两个时间序列所得到的最优参数的欧几里德距离即可。

ARMA(p,q) 模型有自己的 AR 表示,因此可以得到相应的一组参数 (\pi_{1},\pi_{2},\cdots),所以,对于每一条时间序列,都可以用一组最优的参数去逼近。如果

\hat{\prod}_{X_{T}}=(\hat{\pi}_{1,X_{T}},\cdots,\hat{\pi}_{k_{1},X_{T}})^{T},

\hat{\prod}_{X_{T}}=(\hat{\pi}_{1,X_{T}},\cdots,\hat{\pi}_{k_{1},X_{T}})^{T}

分别表示 AR(k_{1})AR(k_{2}) 对于时间序列 X_{T}Y_{T} 的参数估计,则 Piccolo 距离如下:

d_{PIC}(X_{T},Y_{T}) =\sqrt{\sum_{j=1}^{k}(\hat{\pi}_{j,X_{T}}'-\hat{\pi}_{j,Y_{T}}')^{2}},

其中 k=\max(k_{1},k_{2})\hat{\pi}_{j,X_{T}}'=\hat{\pi}_{j,X_{T}}j\leq k_{1},并且 \hat{\pi}_{j,X_{T}}' = 0k_{1}<j\leq k\hat{\pi}_{j,Y_{T}}'=\hat{\pi}_{j,Y_{T}}j\leq k_{2},并且 \hat{\pi}_{j,Y_{T}}' = 0k_{2}<j\leq k

Maharaj 距离

按照之前的描述,可以增加一个矩阵来修改 Piccolo 距离:

d_{MAH}(X_{T},Y_{T}) =\sqrt{T}(\hat{\prod}'_{X_{T}}-\hat{\prod}'_{Y_{T}})^{T}\hat{V}^{-1}(\hat{\prod}'_{X_{T}}-\hat{\prod}'_{Y_{T}}).

其中 \hat{\prod}'_{X_{T}}\hat{\prod}'_{Y_{T}} 表示 AR(k) 模型对于 X_{T}Y_{T} 的参数估计,和 Piccolo 距离一样。\hat{V} = \sigma_{X_{T}}^{2}R_{X_{T}}^{-1}(k) + \sigma_{Y_{T}}^{2}R_{Y_{T}}^{-1}(k)\sigma_{X_{T}}^{2}\sigma_{Y_{T}}^{2} 表示时间序列的方差,R_{X_{T}}R_{Y_{T}} 表示时间序列的 sample covariance 矩阵。

基于 Cepstral 的距离

考虑时间序列 X_{T} 满足 AR(p) 的结构,i.e. X_{t}=\sum_{r=1}^{p}\phi_{r}X_{t-r}+\epsilon_{t},这里的 \phi_{r} 表示 AR 模型的参数,\epsilon_{t} 表示白噪声(均值为 0,方差为 1 的 Gauss 正态分布)。于是可以从这些参数定义 LPC 系数如下:

\psi_{1}=\phi_{1}

\psi_{h}=\phi_{h}+\sum_{m=1}^{h-1}(\phi_{m}-\psi_{h-m})1<h\leq p

\psi_{h}=\sum_{m=1}^{p}(1-\frac{m}{h})\phi_{m}\psi_{h-m}p<h

所以,LPC 的距离定义是:

d_{LPC, Cep}(X_{T},Y_{T}) =\sqrt{\sum_{i=1}^{T}(\psi_{i,X_{T}}-\psi_{i,Y_{T}})^{2}}.

 

总结

在本文中,介绍了时间序列之间距离的计算方法,包括基于 L^{p} 范数的距离,基于相关性的距离,基于周期图表的计算方法,基于模型的计算方法。

 

 

时间序列的表示与信息提取

提到时间序列,大家能够想到的就是一串按时间排序的数据,但是在这串数字背后有着它特殊的含义,那么如何进行时间序列的表示(Representation),如何进行时间序列的信息提取(Information Extraction)就成为了时间序列研究的关键问题。

就笔者的个人经验而言,其实时间序列的一些想法和文本挖掘是非常类似的。通常来说句子都是由各种各样的词语组成的,并且一般情况下都是“主谓宾”的句子结构。于是就有人希望把词语用一个数学上的向量描述出来,那么最经典的做法就是使用 one – hot 的编码格式。i.e. 也就是对字典里面的每一个词语进行编码,一个词语对应着一个唯一的数字,例如 0,1,2 这种形式。one hot 的编码格式是这行向量的长度是词典中词语的个数,只有一个值是1,其余的取值是0,也就是 (0,…,0,1,0,…,0) 这种样子。但是在一般情况下,词语的个数都是非常多的,如何使用一个维度较小的向量来表示一个词语就成为了一个关键的问题。几年前,GOOGLE 公司开源了 Word2vec 开源框架,把每一个词语用一串向量来进行描述,向量的长度可以自行调整,大约是100~1000 不等,就把原始的 one-hot 编码转换为了一个低维空间的向量。在这种情况下,机器学习的很多经典算法,包括分类,回归,聚类等都可以在文本上得到巨大的使用。Word2vec 是采用神经网络的思想来提取每个词语与周边词语的关系,从而把每个词语用一个低维向量来表示。在这里,时间序列的特征提取方法与 word2vec 略有不同,后面会一一展示这些技巧。

 

时间序列的统计特征

提到时间序列的统计特征,一般都能够想到最大值(max),最小值(min),均值(mean),中位数(median),方差(variance),标准差(standard variance)等指标,不过一般的统计书上还会介绍两个指标,那就是偏度(skewness)和峰度(kuriosis)。如果使用时间序列 X_{T} = \{x_{1},\cdots,x_{T}\} 来表示长度为 T 的时间序列,那么这些统计特征用数学公式来表示就是:

\mu=\frac{1}{T}\sum_{i=1}^{T}x_{i},

\sigma^{2} = \sum_{i=1}^{T}\frac{1}{T}(x_{i}-\mu)^{2},

\text{skewness}(X) = E[(\frac{X-\mu}{\sigma})^{3}]=\frac{1}{T}\sum_{i=1}^{T}\frac{(x_{i}-\mu)^{3}}{\sigma^{3}},

\text{kurtosis}(X) = E[(\frac{X-\mu}{\sigma})^{4}]=\frac{1}{T}\sum_{i=1}^{T}\frac{(x_{i}-\mu)^{4}}{\sigma^{4}} .

其中 \mu\sigma 分别表示时间序列 X_{T} 的均值和方差。

 

时间序列的熵特征

为什么要研究时间序列的熵呢?请看下面两个时间序列:

时间序列(1):(1,2,1,2,1,2,1,2,1,2,…)

时间序列(2):(1,1,2,1,2,2,2,2,1,1,…)

在时间序列(1)中,1 和 2 是交替出现的,而在时间序列(2)中,1 和 2 是随机出现的。在这种情况下,时间序列(1)则更加确定,时间序列(2)则更加随机。并且在这种情况下,两个时间序列的统计特征,例如均值,方差,中位数等等则是几乎一致的,说明用之前的统计特征并不足以精准的区分这两种时间序列。

通常来说,要想描述一种确定性与不确定性,熵(entropy)是一种不错的指标。对于离散空间而言,一个系统的熵(entropy)可以这样来表示:

\text{entropy}(X) = -\sum_{i=1}^{\infty}P\{x=x_{i}\}\ln(P\{x=x_{i}\}).

如果一个系统的熵(entropy)越大,说明这个系统就越混乱;如果一个系统的熵越小,那么说明这个系统就更加确定。

提到时间序列的熵特征,一般来说有几个经典的例子,那就是 binned entropyapproximate entropysample entropy。下面来一一介绍时间序列中这几个经典的熵。

Binned Entropy

从熵的定义出发,可以考虑把时间序列 X_{T} 的取值进行分桶的操作。例如,可以把 [\min(X_{T}), \max(X_{T})] 这个区间等分为十个小区间,那么时间序列的取值就会分散在这十个桶中。根据这个等距分桶的情况,就可以计算出这个概率分布的熵(entropy)。i.e. Binned Entropy 就可以定义为:

\text{binned entropy}(X) = -\sum_{k=0}^{\min(maxbin, len(X))} p_{k}\ln(p_{k})\cdot 1_{(p_{k}>0)},

其中 p_{k} 表示时间序列 X_{T} 的取值落在第 k 个桶的比例(概率),maxbin 表示桶的个数,len(X_{T}) = T 表示时间序列 X_{T} 的长度。

如果一个时间序列的 Binned Entropy 较大,说明这一段时间序列的取值是较为均匀的分布在 [\min(X_{T}), \max(X_{T})] 之间的;如果一个时间序列的 Binned Entropy 较小,说明这一段时间序列的取值是集中在某一段上的。

Approximate Entropy

回到本节的问题,如何判断一个时间序列是否具备某种趋势还是随机出现呢?这就需要介绍 Approximate Entropy 的概念了,Approximate Entropy 的思想就是把一维空间的时间序列提升到高维空间中,通过高维空间的向量之间的距离或者相似度的判断,来推导出一维空间的时间序列是否存在某种趋势或者确定性。那么,我们现在可以假设时间序列 X_{N}: \{x_{1},\cdots, x_{N}\} 的长度是 N,同时 Approximate Entropy 函数拥有两个参数,mr,下面来详细介绍 Approximate Entropy 的算法细节。

Step 1. 固定两个参数,正整数 m 和正数 r,正整数 m 是为了把时间序列进行一个片段的提取,正数 r 是表示时间序列距离的某个参数。i.e. 需要构造新的 m 维向量如下:

X_{1}(m) = (x_{1},\cdots, x_{m})\in\mathbb{R}^{m},

X_{i}(m) = (x_{i},\cdots, x_{m+i-1})\in\mathbb{R}^{m},

X_{N-m+1}(m) = (x_{N-m+1},\cdots, x_{N})\in\mathbb{R}^{m}.

Step 2. 通过新的向量 X_{1}(m),\cdots, X_{N-m+1}(m),可以计算出哪些向量与 X_{i} 较为相似。i.e.

C_{i}^{m}(r) = (\text{number of }X_{j}(m)\text{ such that } d(X_{i}(m), X_{j}(m))\leq r)/(N-m+1),

在这里,距离 d 可以选择 L^{1}, L^{2}, L^{p}, L^{\infty} 范数。在这个场景下,距离 d 通常选择为 L^{\infty} 范数。

Step 3. 考虑函数

\Phi^{m}(r) = (N-m+1)^{-1}\cdot \sum_{i=1}^{N-m+1} \ln(C_{i}^{m}(r)),

Step 4. Approximate Entropy 可以定义为:

\text{ApEn}(m,r) = \Phi^{m}(r)-\Phi^{m+1}(r).

Remark.

  1. 正整数 m 一般可以取值为 2 或者 3,r>0 会基于具体的时间序列具体调整;
  2. 如果某条时间序列具有很多重复的片段(repetitive pattern)或者自相似性(self-similarity pattern),那么它的 Approximate Entropy 就会相对小;反之,如果某条时间序列几乎是随机出现的,那么它的 Approximate Entropy 就会相对较大。

Sample Entropy

除了 Approximate Entropy,还有另外一个熵的指标可以衡量时间序列,那就是 Sample Entropy,通过自然对数的计算来表示时间序列是否具备某种自相似性。

按照以上 Approximate Entropy 的定义,可以基于 mr 定义两个指标 AB,分别是

A = \#\{\text{vector pairs having } d(X_{i}(m+1),X_{j}(m+1))<r \text{ of length } m+1  \},

B = \#\{ \text{vector pairs having } d(X_{i}(m), X_{j}(m))<r \text{ of length } m\}.

其中,\# 表示集合的元素个数。根据度量 d(无论是 L^{1}, L^{2}, L^{\infty})的定义可以知道A\leq B,因此 Sample Entropy 总是非负数,i.e.

\text{SampEn} = -\ln(A/B) \geq 0.

Remark.

  1. Sample Entropy 总是非负数;
  2. Sample Entropy 越小表示该时间序列具有越强的自相似性(self similarity)。
  3. 通常来说,在 Sample Entropy 的参数选择中,可以选择 m = 2, r = 0.2 \cdot std.

 

时间序列的分段特征

即使时间序列有一定的自相似性(self-similarity),能否说明这两条时间序列就完全相似呢?其实答案是否定的,例如:两个长度都是 1000 的时间序列,

时间序列(1): [1,2] * 500

时间序列(2): [1,2,3,4,5,6,7,8,9,10] * 100

其中,时间序列(1)是 1 和 2 循环的,时间序列(2)是 1~10 这样循环的,它们从图像上看完全是不一样的曲线,并且它们的 Approximate Entropy  和 Sample Entropy 都是非常小的。那么问题来了,有没有办法提炼出信息,从而表示它们的不同点呢?答案是肯定的。

首先,我们可以回顾一下 Riemann 积分和 Lebesgue 积分的定义和不同之处。按照下面两幅图所示,Riemann 积分是为了算曲线下面所围成的面积,因此把横轴划分成一个又一个的小区间,按照长方形累加的算法来计算面积。而 Lebesgue 积分的算法恰好相反,它是把纵轴切分成一个又一个的小区间,然后也是按照长方形累加的算法来计算面积。

RiemannANDLebesgue

之前的 Binned Entropy 方案是根据值域来进行切分的,好比 Lebesgue 积分的计算方法。现在我们可以按照 Riemann 积分的计算方法来表示一个时间序列的特征,于是就有学者把时间序列按照横轴切分成很多段,每一段使用某个简单函数(线性函数等)来表示,于是就有了以下的方法:

  1. 分段线性逼近(Piecewise Linear Approximation)
  2. 分段聚合逼近(Piecewise Aggregate Approximation)
  3. 分段常数逼近(Piecewise Constant Approximation)

说到这几种算法,其实最本质的思想就是进行数据降维的工作,用少数的数据来进行原始时间序列的表示(Representation)。用数学化的语言来描述时间序列的数据降维(Data Reduction)就是:把原始的时间序列 \{x_{1},\cdots,x_{N}\}\{x_{1}^{'},\cdots, x_{D}^{'}\} 来表示,其中 D<N。那么后者就是原始序列的一种表示(representation)。

分段聚合逼近(Piecewise Aggregate Approximation)— 类似 Riemann 积分

在这种算法中,分段聚合逼近(Piecewise Aggregate Approximation)是一种非常经典的算法。假设原始的时间序列是 C = \{x_{1},\cdots, x_{N}\},定义 PAA 的序列是:\overline{C} = \{\overline{x}_{1},\cdots,\overline{x}_{w}\}

其中

\overline{x}_{i} = \frac{w}{N} \cdot \sum_{j=\frac{N}{w}(i-1)+1}^{\frac{N}{w}i} x_{j}.

在这里 1\leq i\leq w。用图像来表示那就是:

PAA

至于分段线性逼近(Piecewise Linear Approximation)和分段常数逼近(Piecewise Constant Approximation),只需要在 \overline{x}_{i} 的定义上稍作修改即可。

符号逼近(Symbolic Approximation)— 类似 Riemann 积分

在推荐系统的特征工程里面,特征通常来说可以做归一化,二值化,离散化等操作。例如,用户的年龄特征,一般不会直接使用具体的年月日,而是划分为某个区间段,例如 0~6(婴幼儿时期),7~12(小学),13~17(中学),18~22(大学)等阶段。

其实在得到分段特征之后,分段特征在某种程度上来说依旧是某些连续值,能否把连续值划分为一些离散的值呢?于是就有学者使用一些符号来表示时间序列的关键特征,也就是所谓的符号表示法(Symbolic Representation)。下面来介绍经典的 SAX Representation。

如果我们希望使用 \alpha 个符号,例如 \{l_{1},\cdots,l_{\alpha}\} 来表示时间序列。同时考虑正态分布 N(0,1),用\{z_{1/\alpha},\cdots,z_{(\alpha-1)/\alpha}\} 来表示 Gauss 曲线下方的一些点,而这些点把 Gauss 曲线下方的面积等分成了 \alpha 段。

SAX 方法的流程如下:

Step 1. 正规化(normalization):也就是该时间序列被映射到均值为零,方差为一的区间内。

Step 2. 分段表示(PAA):\{x_{1},\cdots, x_{N}\} \Rightarrow  \{\overline{x}_{1},\cdots,\overline{x}_{w}\}

Step 3. 符号表示(SAX):如果 \overline{x}_{i}<z_{1/\alpha},那么 \hat{X}_{i}=l_{1};如果 z_{(j-1)/\alpha}\leq \overline{x}_{i}<z_{j/\alpha},那么 \hat{X}_{i} = l_{j};如果 \overline{x}_{i}\geq z_{(\alpha-1)/\alpha},那么 \hat{X}_{i} = l_{\alpha}

于是,我们就可以用 \{l_{1},\cdots,l_{\alpha}\}\alpha 个字母来表示原始的时间序列了。

SAX

 

总结

在本篇文章中,我们介绍了时间序列的一些表示方法(Representation),其中包括时间序列统计特征,时间序列的熵特征,时间序列的分段特征。在下一篇文章中,我们将会介绍时间序列的相似度计算方法。