Category Archives: Computer Science

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

随着 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 足够小。

总结

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

 

如何从零到一地开始机器学习?

导语:作为一个数学系出身,半路出家开始搞机器学习的人,在学习机器学习的过程中自然踩了无数的坑,也走过很多本不该走的弯路。于是很想总结一份如何入门机器学习的资料,也算是为后来人做一点点微小的贡献。

在 2016 年 3 月,随着 AlphaGo 打败了李世乭,人工智能开始大规模的进入人们的视野。不仅是互联网的工程师们很关注人工智能的发展,就连外面的吃瓜群众也开始关注人工智能对日常生活的影响。随着人脸识别能力的日益增强,个性化新闻推荐 app 的横行天下,TensorFlow 等开源工具被更多的人所知晓,于是就有越来越多的人开始逐步的转行到人工智能的领域,无论是计算机出身的后台开发人员,电子通信等工程师,还是数学物理等传统理科人士,都有人逐步开始转行到机器学习的领域。

作为一个数学系出身,半路出家开始搞机器学习的人,在学习机器学习的过程中自然踩了无数的坑,也走过很多本不该走的弯路。于是很想总结一份如何入门机器学习的资料,也算是为后来人做一点点微小的贡献。

作为一个转行的人,自然要介绍一下自己的专业背景。笔者在本科的时候的专业是数学与应用数学,外行人可以理解为基础数学。在博士期间的研究方向是动力系统和分形几何,所做的还是基础数学,和计算机的关系不大。如果有人想了解笔者究竟在做什么科研的话,可以参考知乎文章“复动力系统(1)— Fatou集与Julia集”。至于机器学习的话,在读书期间基本上也没接触过,甚至没听说过还有这种东西。不过在读书期间由于专业需要,C++ 之类的代码还是能够写一些的,在 UVA OJ 上面也留下过自己的足迹。


2015 年:尝试转型

行路难,行路难,多歧路,今安在?

在 2015 年毕业之后机缘巧合,恰好进入腾讯公司从事机器学习的相关工作。不过刚进来的时候压力也不小,现在回想起来的话,当时走了一些不该走的弯路。用李白的《行路难》中的诗词来描述当时的心情就是“行路难,行路难,多歧路,今安在?”在 2015 年 10 月份,第一次接触到一个不大不小的项目,那就是 XX 推荐项目。而这个项目是当时组内所接到的第二个推荐项目,当年的推荐系统还是搭建在大数据集群上的,完全没有任何说明文档和前端页面,当时的整个系统和全部流程复杂而繁琐。不过在接触这个系统的过程中,逐步开始学习了 Linux 操作系统的一些简单命令,SQL 的使用方法。了解 SQL 的话其实不只是通过了这个系统,通过当时的 ADS 值班,帮助业务方提取数据,也把 SQL 的基础知识进一步的加深了。SQL 的学习的话,在2015年读过两本非常不错的入门教材《SQL基础教程》与《HIVE编程指南》。Linux 的相关内容阅读了《Linux 命令行与 Shell 脚本编程大全》之后也就大概有所了解了。于是工作了一段时间之后,为了总结一些常见的 SQL 算法,写过一篇文章 “HIVE基础介绍“。

在做推荐项目的过程中,除了要使用 SQL 来处理数据,要想做机器学习,还需要了解常见的机器学习算法。当年接触到的第一个机器学习算法就是逻辑回归(Logistic Regression),既然提到了机器学习的逻辑回归,无法避免的就是交叉验证的概念,这个是机器学习中的一个基本概念。通过物品的类别属性和用户的基本特征来构造出新的特征,例如特征的内积(inner product)。后来在学习的过程中逐步添加了特征的外积和笛卡尔积,除了特征的交叉之外,还有很多的方法来构造特征,例如把特征标准化,归一化,离散化,二值化等操作。除了构造特征之外,如何判断特征的重要性则是一个非常关键的问题。最常见的方法就是查看训练好的模型的权重,另外还可以使用 Pearson 相关系数和 KL 散度等数学工具来粗糙的判断特征是否有效。在此期间也写过一些文章“交叉验证”,“特征工程简介”,“KL散度”。关于特征工程,除了阅读一些必要的书籍之外,最重要的还是要实践,只有实践才能够让自己的经验更加丰富。

在做推荐系统的时候,之前都是通过逻辑回归算法(Logistic Regression)离线地把模型的权重算好,然后导入线上系统,再进行实时的计算和打分。除了离线的算法之外,在 2015 年的 12 月份了解到了能够在线学习的 FTRL 算法。调研了之后在 2016 年初在组内进行了分享,同时在 zr9558.com 上面分享了自己的总结,最近把该文章转移到自己的微信公众号上“Follow the Regularized Leader”。

在做 XX 推荐项目的过程中,了解到了数据才是整个机器学习项目的基石,如果数据的质量不佳,那就需要进行数据的预处理,甚至推动开发人员去解决数据上报的问题。通常来说,要想做好一个推荐项目,除了特征工程和算法之外,最重要的就是数据的核对。当时的经验是需要核对多方的数据,那就是算法离线计算出来的结果,线上计算出来的结果,真实产品中所展示的结果这三方的数据必须要完全一致,一旦不一致,就需要复盘核查,而不是继续推进项目。在此期间,踩过无数的数据的坑,因此得到的经验就是一定要反复的核查数据。


2016:从零到一

站在巨人的肩膀上,才能看得更远。-—学习推荐系统

站在巨人的肩膀上,才能看得更远。”到了 2016 年的 2 月份,除了 XX 推荐项目的首页个性化调优算法之外,还开启了另外一个小项目,尝试开启首页的 tab,那就是针对不同的用户推荐不同的物品。这个小项目简单一点的做法就是使用 ItemCF 或者热传导传播的算法,在用户收听过某个节目之后,就给用户推荐相似的节目。这种场景其实在工业界早就有了成功的案例,也不算是一个新的场景。就好比与用户在某电商网站上看中了某本书,然后就被推荐了其他的相关书籍。之前也写过一篇推荐系统的简单算法“物质扩散算法”,推荐给大家参考一下。至于 ItemCF 和热传导算法的相关内容,会在后续的 Blog 中持续完善。

“读书千遍,其义自见。”在使用整个推荐系统的过程中,笔者只是大概知道了整个系统是如何搭建而成的。而要整体的了解机器学习的相关算法,光做项目则是远远不够的。在做推荐业务的这段时间,周志华老师的教材《机器学习》在2016年初上市,于是花了一些时间来阅读这本书籍。但是个人感觉这本书难度不大,只是需要另外一本书结合着看才能够体会其中的精妙之处,那就是《机器学习实战》。在《机器学习实战》中,不仅有机器学习相关算法的原理描述,还有详细的源代码,这足以让每一个初学者从新手到入门了。

路漫漫其修远兮,吾将上下而求索

说到从零到一,其实指的是在这一年体验了如何从零到一地做一个新业务。到了 2016 年的时候,为了把机器学习引入业务安全领域,在部门内部成立了 XX 项目组,这个项目在部门内部其实并没有做过大规模的尝试,也并没有成功的经验,甚至也没有一个合适的系统让人使用,而且安全业务和推荐业务基本上不是一回事。因为对于推荐系统而言,给用户的推荐是否准确决定了 CTR 是否达标,但是对于安全系统而言,要想上线打击黑产的话,准确率则需要 99% 以上才行。之前的推荐系统用得最多的算法就是逻辑回归,而且会存储物品和用户的两类特征,其余的算法主要还是 ItemCF 和热传导算法。这就导致了当时做 XX 项目的时候,之前的技术方案并不可用,需要基于业务安全的实际场景来重新搭建一套框架体系。

但是当时做安全项目的时候并没有实际的业务经验,而且暂定的计划是基于 XX1 和 XX2 两个业务来进行试点机器学习。为了做好这个项目,一开始笔者调研了几家号称做机器学习+安全的初创公司,其中调研的最多的就是 XX 这家公司,因为他们家发表了一篇文章,里面介绍了机器学习如何应用在业务安全上,那就是搭建一套无监督+有监督+人工打标签的对抗体系。笔者还是总结了当时两三个月所学的异常点检测算法,文章的链接如下:“异常点监测算法(一)”,“异常点检测算法(二)”,“异常点检测算法(三)”,“异常点检测算法综述”。

在 2016 年底的时候,说起来也是机缘巧合,有的同事看到了我在 2016 年 11 月份发表的 KM 文章“循环神经网络”,就来找笔者探讨了一下如何构建游戏 AI。当时笔者对游戏AI的应用场景几乎不了解,只知道 DeepMind 做出了 AlphaGo,在 2013 年使用了深度神经网络玩 Atari 游戏。在12月份花费了一定的时间研究了强化学习和深度学习,也搭建过简单的 DQN 网络进行强化学习的训练。通过几次的接触和交流之后总算 2017 年 1 月份做出一个简单的游戏 AI,通过机器学习也能够进行游戏 AI 的自主学习。虽然不在游戏部门,但是通过这件事情,笔者对游戏 AI 也产生了浓厚的兴趣,撰写过两篇文章“强化学习与泛函分析”,“深度学习与强化学习”。


2017 年:再整旗鼓

在做日常项目的同时,在 2017 年也接触量子计算。在后续几个月的工作中,持续调研了量子计算的基础知识,一些量子机器学习的技术方案,写了两篇文章“量子计算(一)”,“量子计算(二)”介绍了量子计算的基础概念和技巧。

三十功名尘与土,八千里路云和月

提到再整旗鼓,其实指的是在 2017 年再次从零到一的做全新的项目。到了 2017 年 7 月份,随着业务安全的机器学习框架已经逐渐完善,XX 项目也快走到了尾声,于是就又有了新的项目到了自己的手里,那就是智能运维项目。运营中心这边还在探索和起步阶段,业界的智能运维(AIOPS)的提出也是在2017年才逐步开始,那就是从手工运维,自动化运维,逐步走向人工智能运维的阶段,也就是所谓的 AIOPS。只有这样,运营中心才有可能实现真正的咖啡运维阶段。

正式接触到运维项目是 2017 年 8 月份,从跟业务运维同学的沟通情况来看,当时有几个业务的痛点和难点。例如:Monitor 时间序列的异常检测,哈勃的根因分析,ROOT 系统的根源分析,故障排查,成本优化等项目。在 AIOPS 人员短缺,并且学术界并不怎么研究这类技术方案的前提下,如何在运维中开展机器学习那就是一个巨大的难题。就像当年有神盾系统,无论怎么做都可以轻松的接入其余推荐业务,并且也有相对成熟的内部经验,学术界和工业界都有无数成功的案例。但是智能运维这一块,在 2017 年才被推广出来,之前都是手工运维和 DevOps 的一些内容。于是,如何尽快搭建一套能够在部门内使用的智能运维体系就成了一个巨大的挑战。面临的难题基本上有以下几点:

1. 历史包袱沉重

2. AIOPS 人员短缺

3. 没有成熟的系统框架

在这种情况下,外部引进技术是不可能了,只能够靠自研,合作的同事主要是业务运维和运营开发。当时第一个接触的智能运维项目就是哈勃的多维下钻分析,其业务场景就是一旦发现了成功率等指标下跌之后,需要从多维的指标中精准的发现异常,例如从运营商,省份,手机等指标中发现导致成功率下跌的原因,这就是经典的根因分析。这一块在调研之后发现,主要几篇文章可以参考,综合考虑了之后撰写了一份资料,那就是“根因分析的探索”。PS:除了哈勃多维下钻之外,个人感觉在 BI 智能商业分析中,其实也可以是这类方法来智能的发现“为什么DAU下跌?”“为什么收入没有达到预期”等问题。

除了哈勃多维下钻之外,Monitor 的时间序列异常检测算法则是更为棘手的项目。之前的 Monitor 异常检测算法,就是靠开发人员根据曲线的特点设定三个阈值(最大值,最小值,波动率)来进行异常检测。这样的结果就是准确率不准,覆盖率不够,人力成本巨大。在上百万条曲线都需要进行异常检测的时候,每一条曲线都需要人工配置阈值是完全不合理的。于是,导致的结果就是每周都需要有人值班,有了问题还不一定能够及时发现。而对于时间序列算法,大家通常能够想到的就是 ARIMA 算法,深度学习的 RNN 与 LSTM 算法,Facebook 近期开源的 Prophet 工具。这些方法笔者都调研过,并且未来会撰写相关的文章介绍 ARIMA,RNN,Prophet 的使用,欢迎大家交流。

其实以上的几种时间序列预测和异常检测算法,主要还是基于单条时间序列来做的,而且基本上是针对那些比较平稳,具有历史规律的时间序列来进行操作的。如果针对每一条曲线都单独搭建一个时间序列模型的话,那和阈值检测没有任何的区别,人力成本依旧巨大。而且在 Monitor 的实际场景下,这些时间序列异常检测模型都有着自身的缺陷,无法做到“百万条KPI曲线一人挑”的效果。于是在经历了很多调研之后,我们创新性的提出了一个技术方案,成功的做到了“百万条曲线”的异常检测就用几个模型搞定。那就是无监督学习的方案加上有监督学习的方案,第一层我们使用无监督算法过滤掉大部分的异常,第二层我们使用了有监督的算法来提升准确率和召回率。在时间序列异常检测的各类算法中,通常的论文里面都是针对某一类时间序列,使用某一类模型,效果可以达到最优。但是在我们的应用场景下,见过的曲线千奇百怪,笔者都说不清楚有多少曲线的形状,因此只用某一类时间序列的模型是绝对不可取的。但是,在学习机器学习的过程中,有一种集成学习的办法,那就是把多个模型的结果作为特征,使用这些特征来训练一个较为通用的模型,从而对所有的 Monitor 时间序列进行异常检测。这一类方法笔者总结过,那就是“时间序列简介(一)”,最终我们做到了“百万条曲线一人挑”,成功去掉了制定阈值的业务效果。


2018年:走向未来

亦余心之所善兮,虽九死其犹未悔。

在转行的过程中,笔者也走过弯路,体会过排查数据问题所带来的痛苦,经历过业务指标达成所带来的喜悦,感受过如何从零到一搭建一套系统。在此撰写一篇文章来记录笔者这两年多的成长经历,希望能够尽微薄之力帮助到那些有志向转行来做机器学习的人。从这两年做项目的经历来看,要想从零到一地做好项目,在一开始就必须要有一个好的规划,然后一步一步的根据项目的进展调整前进的方向。但是如果没有一个足够的知识积累,就很难找到合适的前进方向。

亦余心之所善兮,虽九死其犹未悔。”在某些时候会有人为了短期的利益而放弃了一个长远的目标,但是如果要让自己走得更远,最佳的方案是让自己和团队一起成长,最好的是大家都拥有一个长远的目标,不能因为一些微小的波动而放任自己。同时,如果团队或个人急于求成,往往会导致失败,而坚持不懈的学习则是做科研和开展工作的不二法门。

诗人陆游曾经教育过他的后辈:“汝果欲学诗,功夫在诗外”。意思是说,如果你想真正地写出好的诗词,就要在生活上下功夫,去体验生活的酸甜苦辣,而不是抱着一本诗词歌赋来反复阅读。如果看过天龙八部的人就知道,鸠摩智当时上少林寺去挑战,在少林高僧面前展示出自己所学的少林七十二绝技,诸多少林高僧无不大惊失色。而当时的虚竹在旁边观战,就对少林高僧们说:“鸠摩智所耍的招数虽然是少林绝技,但是本质上却是使用小无相功催动出来的。虽然招数相同,但是却用的道家的内力。”为什么少林的高僧们没有看出来鸠摩智武功的关键之处呢,那是因为少林高僧们在练功的时候,一直抱着武学秘籍在修炼,一辈子练到头了也就13门绝技。其实从鸠摩智的个人修炼来看,修练武学的关键并不在武学秘籍里。没有找到关键的佛经,没有找到运功的法门,无论抱着武学秘籍修炼多少年,终究与别人有着本质上的差距。

笔者在 SNG 社交网络运营部的这两年多,用过推荐项目,做过安全项目,正在做运维项目,也算是部门内唯一一个(不知道是否准确)做过三种项目的人,使用过推荐系统,从零到一搭建过两个系统。目前笔者的个人兴趣集中在 AIOPS 这个场景下,因为笔者相信在业务运维这个传统领域,机器学习一定有着自己的用武之地。相信在不久的将来,AIOPS 将会在运维上面的各个场景落地,真正的走向咖啡运维。

张戎

2018年2月

启发式优化算法

启发式优化算法的一些调研

启发式优化算法的背景

1. 首先,我们需要有一个目标函数 f(\bold{x}), 其中的 \bold{x} \in \mathbb{R}^{n}n 是一个正整数。然后,目标函数 f(\bold{x}) 不一定需要是显式的(所谓显式指的是能够精确写出 f(\bold{x}) 的表达式,例如逻辑回归的目标函数,神经网络的目标函数,XGBoost 的目标函数等),隐式的目标函数使用启发式优化算法也是可以研究的(即无法写出 f(\bold{x}) 表达式)。

2. 其次,需要计算目标函数 f(\bold{x}) 的最大值 \max f(\bold{x}) 或者最小值 \min f(\bold{x})

3. 再次,如果存在多个函数 f_{i}(\bold{x}), 1\leq i\leq m, 并且每一个 f_{i}(\bold{x}) 的值域都在 [0,1] 内。如果目标函数是

\min_{x}(\max_{i}f_{i}(x)-\min_{i}f_{i}(x)),

希望求该目标函数的最小值,那么意思就是使得每一个 1\leq i\leq m, f_{i}(\bold{x})\bold{x} 的定义域内都相差得不多,也就是这 m 条曲线几乎重合。

4. 拥有目标函数的目的是寻找最优的 \bold{x} 使得它满足 argmin_{\bold{x}}f(x) 或者 argmax_{\bold{x}}f(x)

5. 在这里的 \bold{x}\in \mathbb{R}^{n} 一般都是连续的特征,而不是那种类别特征。如果是类别的特征,并不是所有的启发式优化算法都适用,但是遗传算法之类的算法在这种情况下有着一定的用武之地。

启发式优化算法的实验方法

1. 论文中的方法:找到一些形状怪异的函数作为目标函数,需要寻找这些函数的全局最大值或者全局最小值。除了全局的最大最小值点之外,这些函数也有很多局部极值点。例如 Rastrigin 函数(Chapter 6,Example 6.1,Search and Optimization by Metaheuristics)

\min_{\bold{x}}f(\bold{x}) = 10\cdot n +\sum_{i=1}^{n}(x_{i}^{2}-10\cos(2\pi x_{i})),

其中,\bold{x}\in[-5.12,5.12]^{n}. 或者 Easom 函数(Chapter 2,Example 2.1,Search and Optimization by Metaheuristics)

\min_{\bold{x}}f(\bold{x}) = - \cos(x_{1})\cos(x_{2})\exp(-(x_{1}-\pi)^{2}-(x_{2}-\pi)^{2}),

其中,\bold{x} \in [-100,100]^{2}.

2. 最朴素的算法:随机搜索法(Random Search)。也就是把 \bold{x} 的定义域切分成 m^{n} 份,其中 \bold{x}\in\mathbb{R}^{n}. 计算每一个交点的函数取值,然后统计出其最大值或者最小值就可以了。实际中不适用,因为当 n\geq 100 的时候,这个是指数级别的计算复杂度。

3. 启发式优化算法:针对某一个目标函数,使用某种启发式优化算法,然后与随机搜索法做对比,或者与其余的启发式优化算法做对比。

4. 在线调优:启发式优化算法大部分都可以做成实时调优的算法,调优的速度除了算法本身的收敛速度之外,还有给定一个 \bold{x} 之后,获取到 f(\vec{x}) 的时间。

注:可以尝试使用启发式优化算法在线调优推荐系统中的 AUC 或者 CTR。不过 CTR 一般具有波动性,需要几个小时甚至一天的数据才能够统计出一个相对靠谱的取值。同时,如果用户量不够多的话,如果切流量的话,可能会造成一些流量的 CTR 偏高,另外一些流量的 CTR 偏低。并且流量的条数会有上限,一般来说不会超过 100 左右,甚至只有 10 个不同的流量。因此,在其余场景下使用这类方法的时候,一定要注意获取目标函数 value 的时间长度。时间长度越短,能够进行的尝试(探索次数)就会越多,收敛的希望就会越大;反之,时间长度越长,能够进行的尝试(探索次数)就会受到限制,收敛的时间长度就会变长。

启发式优化算法的一些常见算法

以下的内容暂时只涉及最基本的算法,有很多算法的变种暂时不做过多的描述。

粒子群算法(Particle Swarm Optimization)

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] 中间的随机数。

注:需要解决的问题是如何设置这 N_{p} 个粒子,也就是构造粒子群的问题。在每次只能够设置调整一次 \bold{x} 的时候,可以把时间窗口按照连续的 N_{p} 段来进行切分,在一个时间段内的这 N_{p} 个时间点可以当成 N_{p} 个粒子,然后进行下一轮的迭代即可。

模拟退火(Simulated Annealing)

其核心思想是构造了温度 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 足够小。

遗传算法的个人理解:

  1. 如何选择合适的 \Delta \bold{x},可以选择 Gaussian 正态分布,并且可以同时修改 \bold{x}n 个维度,也可以每次只修改一个维度或者部分维度的值;
  2. 一开始在温度 T 很大的时候,即使 E(\bold{x}+\Delta\bold{x})\geq E(\bold{x}) 的时候,都会以极大概率接受 E(\bold{x}+\Delta\bold{x}),因为此时 \lim_{T\rightarrow\infty}e^{-\Delta E/T}=1。然后一开始都在四处游荡;
  3. 当温度 T 逐渐降低的时候,当时 \bold{x} 不一定在一个合适的位置,因此,需要记录历史上探索过的满足 E(\bold{x}) 最小的值,从历史上的值中选择出一个合适的值继续迭代。因此,在温度下降的时候,如何选择一个合适的值继续迭代可能是一个关键问题。
  4. 从遗传算法和 PSO 的对比效果来看,因为 PSO 算法会随时记录当前粒子和所有粒子的最优状态,然后在整个种群中,总会有少数粒子处于当前的最优状态,所以 PSO 的效果从实验效果上来看貌似比不记录状态的遗传算法会好一些。

进化策略(Evolutionary Strategies)

ES 和 GA 有一些方法论上面的不同,这个在书上有论述,等撰写了 GA 在一起看这一块的内容。

进化策略(Evolutionary Strategies)的大体流程与遗传算法(Genetic Algorithm) 相似,不同的是 ES 的步骤是交叉/变异(Crossover/Mutation)在前,选择(Selection)在后。

对于经典的进化策略而言,交叉算子(Crossover Operator)可以按照如下定义:两个父代(parents)\bold{x}_{1}\bold{x}_{2},那么子代(children)的交叉可以定义为:

\bold{x}'=(\bold{x}_{1}+\bold{x}_{2})/2,

这里的相加就是向量之间的相加。

对于经典的进化策略而言,变异算子(Mutation Operator)可以按照如下定义:对于一个向量 \bold{x} = (x_{1},\cdots,x_{n}),使用 Gaussian 正态分布可以构造出一个后代如下:

x_{i}'=x_{i}+N(0,\sigma_{i}),\text{ } i=1,\cdots,n,

这里的 \sigma_{i}, \text{ } i=1,\cdots,n 是基于具体问题的,很难给出一个通用的值。

就经验而谈,ES 算法一般来说比 GA 更加 Robust。在交叉(Crossover)或者变异(Mutation)之后,ES 算法就需要进行选择,通常来说,ES 算法有两种方式,分别是 (\lambda+\mu)(\lambda,\mu) 两种策略。这里的 \mu 表示人口数量(population size),\lambda 表示从人口中产生的后代的数量。

(1)在 (\lambda+\mu) 策略中,从 \mu 人口中通过变异或者交叉得到 \lambda 个人,然后从这 (\lambda+\mu) 个人中选择出 \mu 个最优的人。(\lambda+\mu) 是精英策略,每一次迭代都比上一次要好。

(2)在 (\lambda, \mu) 策略中,\mu 个最优的人是从 \lambda(\lambda\geq\mu) 中选择出来的。

传统的梯度搜索法,例如 Newton-Raphson 方法,需要能够计算目标函数 f(x) 的导数,也就是说目标函数必须可导。但是在 Evolutionary Gradient Search 中,不需要 f(x) 是可导的,而是使用导数的近似值来计算。

如果是 (1,\lambda)-ES 的方案,它整体来说其实只有一个人口。从当前的点 \bold{x} 开始,通过 Gaussian Mutation 的方法可以生成 \lambda 个后代,记为 \bold{t}_{1},\cdots,\bold{t}_{\lambda}。计算它们的函数取值 f(\bold{t}_{1}),\cdots,f(\bold{t}_{\lambda})。那么估算的梯度就可以计算出来:

\bold{g}=\sum_{i=1}^{\lambda}(f(\bold{t}_{i})-f(\bold{x}))(\bold{t}_{i}-\bold{x}),

归一化之后得到 \bold{e}=\bold{g}/||\bold{g}||.

Evolutionary gradient search 生成了两个点,分别是:

\bold{x}_{1}=\bold{x}+(\sigma \psi)\bold{e}, \bold{x}_{2}=\bold{x}+(\sigma/\psi)\bold{e},

这里的 \psi>1. 新的个体定义为: \bold{x}'=\bold{x}+\sigma'\bold{e}. 如果 f(\bold{x}_{1})>f(\bold{x}_{2}), 那么 \sigma' = \sigma\psi; 如果 f(\bold{x}_{1})\leq f(\bold{x}_{2}), 那么 \sigma'=\sigma/\psi.

进化策略(ES)包括 gradient search 和 gradient evolution,并且 covariance matrix adaptation(CMA)ES 在一定的限制条件下加速了搜索的效率。

参考资料:

  1. Search and Optimization by Metaheuristics,Kelin DU, M.N.S. SWAMY,2016年

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

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 中,笔者将会介绍时间序列的更多相关内容。