在新加坡的这五年---助教篇

很多 PHD 在国外攻读博士学位的时候,除了做科研和发表论文,还要负责一定量的教学工作,毕竟这些 PHD 以后是非常有可能去大学或者科研机构去任教的。回想当年(2009 年,i.e. 将近十年前),收到新加坡国立大学 Offer 的时候,除了一些 Congratulations 的话语和愿意提供给我 Offer 的信件之外,最重要的就是被数学系告知每周都需要做一定量的教学工作。在攻读博士学位期间,笔者不仅给本科生修改过作业,还给本科生上过习题课;而每次到了期中期末的时候,还要负责监考。在新加坡读书那五年,当助教的过程中,身边发生了不少有趣的事情,时隔多年也不知道能否一一想起,于是提笔写下此文,希望多年之后还能够想起那么多有趣的人和有趣的事。

KentRidgeRoad

NUS 的 Kent Ridge Road,Block S17 过街天桥。一条走了 3000 次的路。

助教开始

刚进入 NUS 的时候,也就是 2010 年秋季。低年级的 PHD 们自然还处于一个适应期和修课的阶段,于是系里面也没有给低年级的 PHD 们很重的工作量,也没有让低年级的 PHD 们去教本科生的习题课。于是我们拿到的任务就是修改作业,同时授课老师会给我们提供一份作业答案,我们只需要照着那份答案修改即可。其实修改作业还算好,没有给笔者造成很大的压力。当时压力比较大的就是给本科生上 Maple(一种数学软件)。说到压力大并不是说这门课有多难,而是要英语给本科生讲这门课。作为一个英语一直不太好的人,当时的听说读写还是存在很大的困难,何况在刚到新加坡没多久的情况下就要给本科生讲习题课。但是,这些都是数学系给的任务,不做的话就没有办法顺利拿到学位,一切都只能够硬着头皮上。当时教 Maple 的时候不仅需要给学生讲一些基础概念,还要负责学生们的上机操作。通常来说,学生在运行程序出 Bug 的时候除了自己 Debug 之外,更多的则是直接问我们各种各样的问题。所以,在教课的过程中需要不停地给学生们查看 Maple 程序里面的 bug,解答各种奇怪的数学问题。有一次在教课的过程中,有个学生对某个习题提出了自己的意见,不过当时自己没有回答出来,当时也不知道是自己的英语问题还是其他的问题。不过事后还是重新梳理了思路,给了学生一个合理的解释。在教 Maple 的过程中,除了学到了一些新的知识之外,最大的收获应该就是自己的英语口语能力的到了提升,并且能够独立地讲习题课。

当年刚到 NUS 的时候还有一些小插曲,由于当时数学系刚从 Block S14 搬到 Block S17,办公室的座位十分紧张,并不能够保证所有 PHD 都有办公座位。因此就需要通过抽签的方式来决定哪些人有座位,哪些人没有座位。而笔者非常悲剧的通过两次抽签,成为了当年系里面唯一一个没有自己座位的 PHD。所以笔者在第一个学期每次从老师那里拿到作业的时候,都只能够拿到 Science Library 去修改,然后带回 Block 602(当时在外租的房子)。只有到了需要归还给教课老师的时候,再从 Block 602 拿到系里面。所以,当时修改作业给我的最大印象就是把学生们的作业搬来搬去是非常累的一件事情。话说,在当年没有办公室座位的时候,笔者只能去 Science Library 去看书。但是,自从有了办公室之后,除了日常的借书之外,就几乎没有去过 Science Library 了。

ScienceLibrary2

Science Library 的门口

ScienceLibrary1

Science Library 中的 GTM 系列

在修改了一年的作业之后,通过一年多的课程的学习,除了数学知识的增加之外,自己的英语能力也得到了相应的提高,至少用英语讲几节习题课是没有太大的问题了。在 2011 年的时候,应该是在 2011 年的第一个学期,学校组织了一次教学的培训,那就是 CDTL(Centre for Development of Teaching & Learning)那边组织的提供了两天的培训,系里面不少的同学都去参加了那个培训。在培训的第一天都是在听那些教课的老师传授经验,至今能够记住的只有一句话:“每个人能够专注的时间非常短”。第二天则是所有来参加培训的 PHD 们都有一次十分钟的练习机会,然后老师和其余的 PHD 们都会来提供意见与建议。当年 CDTL 培训完之后,给每个学员提供了一个杯子,记得当时我还问培训老师这个我能否带走,老师给出了肯定的答案。于是,这个杯子就被我从 CDTL 拿到了 S17,然后从新加坡带回国,最后放在了我目前的办公桌上。

IMG_1666

CDTL 的杯子

除了学校提供的助教培训课程之外,数学系主管教学的老师们也会给学生们提供培训课程。于是,到了第二年,系里面主管教学的 VT 老师就给第二年的 PHD 们群发邮件,告诉大家系里面将会组织一次培训,每个 PHD 都收到了一份习题和相应的答案。然后大家将会在一个下午现场操练一下,分别讲一讲那份习题,VT 老师会给大家相应的建议。在这次培训之后,大家就算正式接受过学校和院系的助教培训了,就要给学生们带习题课了。

助教生涯

说到第二年助教,那是自己唯一一次教线性代数的课程,Linear Algebra,MA1102。当时所用的教材是 NUS 的几位老师所写的,作者有 VT,NG 等。虽然这门课的难度并没有国内的难度大,但是教这门课确实也花费了不少的功夫。因为选择这门课的学生不仅有数学系和工程系的学生,甚至还有文科的学生。对于理工科的学生来说,这门课可能难度还行;但是对于文科,例如 Arts 之类的学生来说,这门课的难度估计已经远远超过了他们的想象。记得当时也有一些来自文科院系的学生反映了学这门课比较困难,所以每次这些学生遇到困难的时候,我总是会在下课之后多给他们讲一些内容。一来学生选了这门课就要拿到相应的学分,二来有的学生确实学习态度十分端正。说来说去,线性代数也算是自己第一次带一门课,可能对于老师来说,第一次独自带一门课,对学生总是会好很多。现在回想起来,通过这门课也认识了一些很不错的朋友,虽然毕业之后大家都再也没有见过面,不过当年能够在一起共同学习一门课也是一种缘分。

linear_algebra_i_ma_siu_lin_victor_tan_and_ng_kah_loon_1516770173_954ef566 (1)

Linear Algebra 的教材封面

数学系的老师除了教授本专业的课程之外,最重的课程应该就是 MA1505 和 MA1506。这两门课的名字叫做 Mathematics (I) 和 Mathematics (II),按照国内的教学内容来看就是大学数学。为什么说课程量重呢?因为这两门课是面向全校所有院系的,不仅有数学系的人会听课,还有计算机系(School of Computing)的人,甚至还有整个工程院系(Faculty of Engineering)的本科生。所以每个学期的学生人数都是一千多人,教课的时候老师们都只能够选择那种很大的课堂,然后分成几个很大的班级,分别教授同样的课程。记得每次教这两门课的时候都要用数学系的四个老师去教,而助教的人数就更多了,每次都要花费数十个 PHD 才能够提供足够多的习题课供学生们选择。而数十个 PHD 则几乎是一个年级的全部 PHD,所以,随着后来 PHD 们数量的减少,教学的任务就变得越来越重。在 2011 年的时候,每个 PHD 只需要教课 4 小时可以了,但是到了 2014 年的时候,每个 PHD 都必须有 5 小时以上的工作量才能够把整个习题课解决掉。

IMG_1660

贴一张当年 MA1505 课程安排的工作量

因为要教 4-5 小时的习题课,于是大家都倾向于把习题课安排在同一天,因为这样的话只需要去一天就好了。这样虽然能够减少去的次数,不过教四节课那天实在是非常累,讲完之后基本上都不愿意继续说话。在 UTown 没有修好之前,MA1505 通常都会安排在 Faculty of Engineering;而在 UTown 修好之后,MA1505 就改到了 UTown 的教室。除了习题课之外,MA1505 通常都会给学生们提供答疑的时间,一般来说每个人每周去一个小时就可以了。有一次,我和富贵,何老师同时教一门课,于是三个人一合计就把时间选择在了某一天的早上,从 9:00-12:00,每次一个人连去三小时,一个学期算下来虽然工作的时间量没变化,但是去的次数少了很多。由于 UTown 或者 Engineering 距离 NUS 的游泳池非常近,于是自己在教课之前或者教课之后会选择去游泳池健身。不过整体来说,NUS 的助教的工作量还是属于能够接受的范围内,一般来说不会占用一个 PHD 很多的时间。

IMG_1658

上面这幅图是在 2013 年当助教的时候,在 Faculty of Engineering 买的午饭,也是当年发的第一个微信朋友圈。

IMG_1659

当年 MPSH 的 Swimming Pool

通过当年的助教,确实把之前所学的很多内容重新学习了一遍,虽然之前在本科的时候学过数学分析和高等代数,但是多年过去难免忘记很多细节。借助助教的机会重温了当年的许多课程,也算是把自己多年的所学进一步巩固了一下。时至今日,一旦有一些空闲时间,笔者就会翻开之前所学的数学书,看一看当年的内容是否已经生疏,是否需要重新学习一遍。

助教故事

提到当 NUS 数学系的助教,就不能够不提到大峰哥。当年在 VT 给大家培训的时候,我们就发现只有大峰哥把习题全部自己做了一遍,然后逐字逐句给大家解释得很清楚,同时 VT 老师也对大风歌给出了非常高的评价。于是,在未来几年里面,大峰哥都主要在给 VT 老师当助教。大峰哥在数学系做助教的时候,尽职尽责,不仅能够做到认真备课,还把课堂上的所有 PPT 和资料自行整理了一遍。这里的整理并不只是把资料整理好给学生,而是把所有的课堂要点和考试关键点用 LaTex 整理出来,并整理成 PDF 发给学生当学习资料。而且每次大峰哥在考试之前都开设专门的答疑时间,给全体学生们的考试答疑解惑。最终的结果就是大峰哥连续 N 次拿到 Best Tutor 的职位,每次的学生评分都在 4.5 以上,据说有一次的评分达到了 4.8(满分是5)。记得之前笔者讲得最好的时候貌似也就是 3.7 左右的得分,从来没有上过 4 分,与大峰哥的得分形成了鲜明的对比。

BestTutor1

NUS 数学系的 Best Tutor 们

虽然说笔者从来没有拿到过 Best Tutor,但是自认为教课也算认真负责。那几年只有一次讲线性代数,一次讲数学物理方法。剩下的时间全部是讲 MA1505 和 MA1506,每一学期几乎都是面对相同的教学内容,只是学生换了一波又一波。因此,笔者在 2013 年秋整理了 MA1505 的讲义和资料,不仅使用 WordPress 整理在个人主页上(zr9558.com),还在习题课的时候打开自己的主页,使用所做的资料给学生们授课。除此之外,在每次考试之前,无论是期中考试还是期末考试,都会给学生们划出考试重点,甚至在考试之前还会写出考题预测(虽然也没有命中某个题目,但是足以命中某类题目)。即使笔者那么做,也并没有提高学生的评分。不过,在教 MA1505 和 MA1506 的时候,还是认识了不少优秀的学生。当时来新加坡就读的大部分都是 SM2,SM3 这些政府项目,也就是说这些大陆来的学生在毕业之后是需要留在新加坡工作几年的,需要履行当年所签下的合同。

TutorialsMA1505

以上图片是当年在 zr9558.com 上面给学生所整理的资料

除此之外,还有许多有趣的人和有趣的事情尚未写下来,后续笔者会持续更新当年求学路上的点点滴滴。

 

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

总结

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

 

Riemann Zeta 函数(二)

在上一篇文章里面,我们已经给出了 Riemann Zeta 函数的定义,

\zeta(s) = \sum_{n=1}^{\infty} \frac{1}{n^s}.

其定义域是 [1,\infty)\subseteq\mathbb{R}. 根据级数与定积分的等价关系可以得到:

  1. s = 1 时,\zeta(1) = \infty;
  2. s>1 时,\zeta(s)<\infty.

本文将会重点讲两个内容:

  1. 如何把 Riemann Zeta 函数从 [1,\infty)\subseteq \mathbb{R} 上延拓到 \{s\in \mathbb{C}: \Re(s)>0\} 上;
  2. Riemann Zeta 函数在 \{s\in\mathbb{C}: \Re(s)\geq 1\} 上没有零点。

Riemann Zeta 函数定义域的延拓

如果想把 Riemann Zeta 函数的定义域从 [1,\infty)\subseteq \mathbb{R} 延拓到更大的区域 \{s\in\mathbb{C}:\Re(s)>0\} 上,就需要给出 Riemann Zeta 函数在 \{s\in \mathbb{C}: \Re(s)>0\} 上的定义。而且在原始的定义域 [1,\infty)\subseteq\mathbb{R} 上面,新的函数的取值必须与原函数的取值保持一致。

首先,我们将会在 [1,\infty)\subseteq \mathbb{R} 上面证明如下恒等式:

\zeta(s) = \frac{s}{s-1} - s\int_{1}^{\infty}\frac{\{x\}}{x^{s+1}}dx.

证明:当 s=1 时,上述等式显然成立,两侧都是 \infty.

\frac{s}{s-1}-s\int_{1}^{\infty}\frac{\{x\}}{x^{s+1}}dx

= \frac{s}{s-1} - s\sum_{n=1}^{\infty}\int_{n}^{n+1}\frac{\{x\}}{x^{s+1}}dx 

= \frac{s}{s-1} - s\sum_{n=1}^{\infty}\int_{n}^{n+1}\frac{x-n}{x^{s+1}}dx 

= \frac{s}{s-1} - s\sum_{n=1}^{\infty}\bigg(\int_{n}^{n+1}\frac{1}{x^{s}}dx - \int_{n}^{n+1}\frac{n}{x^{s+1}}dx\bigg)

= \frac{s}{s-1} - s\int_{1}^{\infty}\frac{1}{x^{s}}dx + \sum_{n=1}^{\infty}n\cdot\int_{n}^{n+1}\frac{s}{x^{s+1}}dx

= \sum_{n=1}^{\infty}n\cdot\bigg(\frac{1}{n^{s}}-\frac{1}{(n+1)^{s}}\bigg)

= \sum_{n=1}^{\infty}\bigg(\frac{1}{n^{s-1}}-\frac{1}{(n+1)^{s-1}} + \frac{1}{(n+1)^{s}}\bigg)

= \sum_{n=1}^{\infty}\frac{1}{n^{s}}.

从右式的表达式

\frac{s}{s-1} - s \int_{1}^{\infty}\frac{\{x\}}{x^{s+1}}dx

可以看出 \zeta(s) 可以延拓到 \{s \in\mathbb{C}:\Re(s)>0\} 上。而且右侧的函数在 \{s\in\mathbb{C}:\Re(s)>0,s\neq 1\} 是解析的,并且 s=1 是该函数的一个极点。进一步的分析可以得到,我们得到一个关于 (s-1)\zeta(s) 的解析函数,而且 \lim_{s\rightarrow 1}(s-1)\zeta(s)=1. 综上所述:

  1. Riemann Zeta 函数可以延拓到 \{s\in\mathbb{C}:\Re(s)>0\} 上;
  2. Riemann Zeta 函数在 \{s\in\mathbb{C}:\Re(s)>0, s\neq 1\} 上是解析的;s=1 是 Riemann Zeta 函数的极点。

 

Riemann Zeta 函数的非零区域

著名的 Riemann 猜想说的是 \zeta(s) 函数的所有非平凡零点都在直线 \{s\in\mathbb{C}:\Re(s)=1/2\} 上。因此,数学家首先要找出的就是 Riemann Zeta 函数的非零区域。而本篇文章将会证明 Riemann Zeta 函数在 \{s\in\mathbb{C}:\Re(s)\geq 1\} 上面没有零点。

\Re(s)>1 区域

首先,我们要证明当 \Re(s)>1 时,\zeta(s)\neq 0.

在这里,就需要使用一个重要的恒等式:当 \Re(s)>1 时,

\zeta(s) =\sum_{n=1}^{\infty}\frac{1}{n^{s}}

= \prod_{p}\bigg(1+\frac{1}{p^{s}}+\frac{1}{p^{2s}}+\cdots\bigg)

= \prod_{n=1}^{\infty}\bigg(1-\frac{1}{p_{n}^{s}}\bigg)^{-1},

其中这里的 p 表示所有的素数相乘,而 p_{n} 表示第 n 个素数。

下面我们证明:

\bigg|1-\frac{1}{p_{n}^{s}}\bigg|^{-1}\geq 1-\frac{1}{p_{n}^{\sigma}-1} .

事实上,令 s = \sigma + i t,,当 \sigma=\Re(s)>1 时,我们有

\bigg|1-\frac{1}{p_{n}^{s}}\bigg|^{-1} = \bigg(1+\frac{1}{p_{n}^{s}}+\frac{1}{p_{n}^{2s}}+\cdots\bigg)

\geq 1-\frac{1}{|p_{n}^{s}|}- \frac{1}{|p_{n}^{2s}|} -\cdots

= 1- \frac{1}{p_{n}^{\sigma}} - \frac{1}{p_{n}^{2\sigma}} -\cdots

= 1- \frac{1}{p_{n}^{\sigma}-1}.

因此,

|\zeta(s)| \geq \prod_{n=1}^{\infty}\bigg|1-\frac{1}{p_{n}^{s}}\bigg|^{-1} \geq\prod_{n=1}^{\infty}\bigg(1-\frac{1}{p_{n}^{\sigma}-1}\bigg).

同时,

\lim_{n\rightarrow \infty} \bigg(1- \frac{1}{p_{n}^{\sigma}-1}\bigg) = 1 ,

1-\frac{1}{p_{n+1}^{\sigma}-1} \geq 1- \frac{1}{p_{n}^{\sigma}-1} ,

\sum_{n=1}^{\infty}\frac{1}{p_{n}^{\sigma}}\leq \sum_{n=1}^{\infty}\frac{1}{n^{\sigma}}<\infty when \sigma>1.

所以,当 \Re(s)>1 时,\zeta(s) \neq 0.

\Re(s) =1 直线

Claim 1. 下面我们将会证明恒等式:对于 \sigma >1, \text{ } t\in\mathbb{R},

\Re(\ln\zeta(\sigma + it)) = \sum_{n=2}^{\infty}\frac{\Lambda(n)}{n^{\sigma}\ln(n)}\cos(t\ln(n)) ,

其中当 n 形如 p^{\alpha}, p 是素数,\alpha \geq 1. \Lambda(n) = \ln(p). 而对于其余的 n, \Lambda(n)=0.

事实上,根据 Euler 公式,

\zeta(s) = \prod_{p}\bigg(1-\frac{1}{p^{s}}\bigg)^{-1}.

s = \sigma + it, 可以得到

\ln\zeta(s) = -\sum_{p}\ln\bigg(1-\frac{1}{p^{s}}\bigg)

= \sum_{p}\sum_{\alpha=1}^{\infty}\frac{1}{\alpha p^{\alpha s}}

= \sum_{p}\sum_{\alpha=1}^{\infty}\frac{1}{\alpha p^{\alpha\sigma}}\cdot p^{-i\alpha t}

= \sum_{p}\sum_{\alpha = 1}^{\infty}\frac{1}{\alpha p^{\alpha\sigma}}\cdot e^{-i\alpha t \ln p}

进一步,

\Re(\ln\zeta(s)) = \sum_{p}\sum_{\alpha =1}^{\infty}\frac{1}{\alpha p^{\alpha\sigma}}\cos(\alpha t \ln p)

并且右侧等于

RHS = \sum_{n=2}^{\infty}\frac{\Lambda(n)}{n^{\sigma}\ln(n)}\cos(t\ln(n))

= \sum_{p}\sum_{\alpha = 1}^{\infty} \frac{\ln(p)}{p^{\alpha\sigma}\ln(p^{\alpha})}\cos(t\ln(p^{\alpha}))

= \sum_{p}\sum_{\alpha = 1}^{\infty}\frac{1}{\alpha p^{\alpha\sigma}}\cos(\alpha t\ln p).

所以,恒等式成立,Claim 1 证明完毕。

Claim 2.

\Re(3\ln\zeta(\sigma) + 4\ln\zeta(\sigma+it) + \ln\zeta(\sigma+2it))\geq 0,

其中 \sigma>1, t\in\mathbb{R}. 换句话说

|\zeta(\sigma)^{3}\zeta(\sigma+it)^{4}\zeta(\sigma+2it)|\geq 1.

事实上,

从三角函数的性质可以得到:

3+4\cos(\theta)+\cos(2\theta) = 3 + 4\cos(\theta)+2\cos^{2}(\theta)-1

= 2(\cos(\theta)-1)^{2}\geq 0,

所以,从 Claim 1 可以得到

\Re(3\ln\zeta(\sigma) + 4\ln\zeta(\sigma+it) + \ln\zeta(\sigma+2it))

= \sum_{n=2}^{\infty} \frac{\Lambda(n)}{n^{\sigma}\ln(n)} \cdot ( 3 + 4\cos(t\ln(n)) + \cos(2t\ln(n))) \geq 0.

进一步地,使用 \Re(\ln(z)) = \ln(|z|) 可以得到

0\leq 3\ln|\zeta(\sigma)| + 4\ln|\zeta(\sigma+it)| + \ln|\zeta(\sigma+2it)|

= \ln|\zeta(\sigma)^{3}\zeta(\sigma+it)^{4}\zeta(\sigma+2it)|,

可以推导出 |\zeta(\sigma)^{3}\zeta(\sigma+it)^{4}\zeta(\sigma+2it)|\geq 1. 因此 Claim 2 证明完毕。

Claim 3. \zeta(1+it)\neq 0 对于所有的 \{t\in\mathbb{R}: t\neq 0\} 成立。

反证法:假设 \zeta(s)s=\sigma + it (t\neq 0) 存在阶数为 m 的零点。也就是说:

\lim_{\sigma\rightarrow 1^{+}} \frac{\zeta(\sigma+it)}{(\sigma+it-1)^{m}}=c\neq 0, 其中 m\geq 1.

从 Riemann Zeta 函数的延拓可以知道,\lim_{\sigma\rightarrow 1^{+}}(\sigma -1)\zeta(\sigma) = 1. 并且 \zeta(s)\{s\in\mathbb{C}:\Re(s)>0, s\neq 1\} 上是解析函数。

从 Claim 2 可以得到:

|(\sigma-1)^{3}\zeta(\sigma)^{3}(\sigma+it-1)^{-4m}\zeta(\sigma+it)^{4}\zeta(\sigma+2it)|

\geq |\sigma-1|^{3}|\sigma-1+it|^{-4m}

\geq |\sigma-1|^{3}\cdot |\sigma-1|^{-4m}

= \frac{1}{|\sigma-1|^{4m-3}}.

\sigma\rightarrow 1^{+}, 可以得到左侧趋近于一个有限的值,但是右侧趋近于无穷,所以得到矛盾。也就是说当 t\neq 0 时, \zeta(1+it)\neq 0 成立。

根据之前的知识,s= 1\zeta(s) 的极点,所以我们得到了本篇文章的主要结论:\zeta(s)\{s\in\mathbb{C}:\Re(s)\geq 1\} 上面没有零点。

 

总结

本篇文章从 Riemann Zeta 函数的延拓开始,证明了 Riemann Zeta 函数在 \{s\in\mathbb{C}:\Re(s)\geq 1\} 上没有零点。在下一篇文章中,笔者将会证明在 \Re(s)=1 附近一个“狭长”的区域上,Riemann Zeta 函数没有零点。

 

从调和级数到 RIEMANN ZETA 函数(一)

Riemann Zeta 函数

Riemann Zeta 函数(Riemann zeta function),\zeta(s),是一个关于复数 s 的方程。在复平面上,当复数 s 的实数部分 \sigma=\Re s >1 时,\zeta(s) 就是如下的级数形式:

\zeta(s) = \sum_{n=1}^{\infty}\frac{1}{n^{s}}.

调和级数的概念与性质

既然提到了级数,首先让我们来回顾一下级数的定义是什么?

级数的定义:在数学中,一个有穷或者无穷的序列 (x_{0},x_{1},x_{2},...) 的形式和 S = x_{0}+x_{1}+x_{2}+... 称为级数,里面的每一项都称为级数的通项。

级数收敛的定义:令 S_{n}=x_{0}+...+x_{n},如果存在有限的 S 使得 \lim_{n\rightarrow \infty}S_{n}=S,那么就称该级数收敛。否则,该级数就称为发散级数。

然后下面我们来研究一下调和级数的基本性质。调和级数的表达式写出来十分简单,那就是 Riemann Zeta 函数在 s=1 的取值,i.e.

\zeta(1) = \sum_{n=1}^{+\infty}\frac{1}{n}.

提到级数的收敛或发散,就必须要提到关于级数收敛的等价定理(Cauchy 判别法),那就是:级数 S_{n} 收敛当且仅当对任意的 \epsilon>0,存在 N 使得对于任意的 m, n>N 都有 |S_{m}-S_{n}|<\epsilon.

既然是等价定理,那么就可以使用 Cauchy 判别法来判断调和级数是否收敛。

Method 1.

S_{n}=\sum_{k=1}^{n}\frac{1}{k},

直接通过计算得到

|S_{2n}-S_{n}|=\frac{1}{n+1}+...+\frac{1}{2n}>\frac{1}{2n}+...+\frac{1}{2n}=\frac{1}{2},

说明该级数是不收敛的,也就是调和级数是发散的。

除了基于 Cauchy 收敛准则的证明之外,能否写出判断调和级数发散的其他方法呢?答案是肯定的。以下有一种使用初等数学方法就能够解释调和级数发散的方法。

Method 2.

\sum_{n=1}^{+\infty}\frac{1}{n}

=1+\frac{1}{2}+(\frac{1}{3}+\frac{1}{4})+(\frac{1}{5}+\frac{1}{6}+\frac{1}{7}+\frac{1}{8})+...

>1+\frac{1}{2}+(\frac{1}{4}+\frac{1}{4})+(\frac{1}{8}+\frac{1}{8}+\frac{1}{8}+\frac{1}{8})+...

=1+\frac{1}{2}+\frac{1}{2}+\frac{1}{2}+...=+\infty.

既然都提到了高等数学,那么当然不能仅仅局限于使用初等数学的技巧来解决问题。而且如果只是用初等数学的方法,在拓展性方面就会受到极大的限制。

Method 3. 调和级数的发散可以通过定积分的技巧来进行解决。

HarmonicSeries

1+\frac{1}{2}+...+\frac{1}{n}

>\int_{1}^{2}\frac{1}{x}dx + \int_{2}^{3}\frac{1}{x}dx+...+\int_{n}^{n+1}\frac{1}{x}dx

=\int_{1}^{n+1}\frac{1}{x}dx=\ln(n+1)

因此,\sum_{n=1}^{\infty}\frac{1}{n}=+\infty.

从上面的定积分的方法可以预计出调和级数的量级大约是对数的量级,那么能否精确的估计出来呢?例如下面这个问题:

问题:\lim_{n\rightarrow +\infty}\frac{\sum_{k=1}^{n}\frac{1}{k}}{ln(n)}=?

通过 L’Hospital 法则可知:\lim_{x\rightarrow 0}x/\ln(1+x)=1.

通过 Stolz 定理可知:

\lim_{n\rightarrow +\infty}\frac{\sum_{k=1}^{n}\frac{1}{k}}{ln(n)}

= \lim_{n\rightarrow +\infty}\frac{\frac{1}{n}}{\ln(n/(n-1))}

= \lim_{x\rightarrow 0}\frac{x}{\ln(1+x)}=1

除此之外,我们同样可以证明

\lim_{n\rightarrow+\infty}(1+\frac{1}{2}+...+\frac{1}{n}-\ln(n))

这个极限是存在并且有限的。

调和级数的推广

那么,如果在考虑 \zeta(2) 也就是级数

\zeta(2) = \sum_{n=1}^{\infty}\frac{1}{n^{2}}

是否收敛的时候,能否用到以上类似的技巧呢?首先,确实也存在各种各样的初等数学技巧,例如:

Method 1.

\sum_{n=1}^{+\infty}\frac{1}{n^{2}}<1+\sum_{n=2}^{+\infty}\frac{1}{n(n-1)}=1+\sum_{n=2}^{+\infty}(\frac{1}{n-1}-\frac{1}{n})=2.

Method 2. 使用数学归纳法。也就是要证明:

\sum_{k=1}^{n}1/k^{2}\leq 2-\frac{1}{n}.

n=1 的时候,公式是正确的。假设 n 的时候是正确的,那么我们有\sum_{k=1}^{n}1/k^{2}\leq 2-\frac{1}{n}。计算可得:

\sum_{k=1}^{n+1}\frac{1}{k^{2}}

<2-\frac{1}{n}+\frac{1}{(n+1)^{2}}

= 2- \frac{1}{n+1}-\frac{1}{n(n+1)^{2}}

\leq 2-\frac{1}{n+1}.

因此,不等式正确,所以 \sum_{n=1}^{+\infty}1/n^{2} 收敛。

其次,在判断调和级数发散的时候,使用的定积分的方法同样可以应用在这个场景下。

Method 3.

1+\frac{1}{2^{2}}+...+\frac{1}{n^{2}}

<1+\int_{1}^{2}\frac{1}{x^{2}}dx+...+\int_{n-1}^{n}\frac{1}{x^{2}}dx

=1+\int_{1}^{n}\frac{1}{x^{2}}dx=1+1-\frac{1}{n}<2.

那么这个是针对次数等于2的情况,对于一般的情形,

\zeta(s)=\sum_{n=1}^{+\infty}\frac{1}{n^{s}},\sigma = \Re(s)>1.

使用定积分的技术,同样可以证明对于任意的 \sigma = \Re(s)>1,都有 \zeta(s) 是收敛的。但是 \zeta(1) 是发散的。

Riemann Zeta 函数中某些点的取值

除此之外,既然 \zeta(s)\sigma = \Re(s)>1 的时候收敛,能否计算出某些函数的特殊值呢?答案是肯定的,例如,我们可以使用 Fourier 级数来计算出 \zeta(2), \zeta(4), \zeta(6),... 的取值。首先,我们回顾一下 Fourier 级数的一些性质:

假设 f(x) 是一个关于 2\pi 的周期函数, i.e. f(x)=f(x+2\pi) 对于所有的 x \in \mathbb{R} 都成立。那么函数 f(x) 的 Fourier 级数就定义为

a_{0}+\sum_{n=1}^{\infty} (a_{n} \cos(nx) +b_{n} \sin(nx)),

其中,a_{0}= \frac{1}{2\pi} \int_{-\pi}^{\pi} f(x) dx,

a_{n}= \frac{1}{\pi} \int_{-\pi}^{\pi} f(x) \cos(nx) dx n\geq 1,

b_{n}= \frac{1}{\pi} \int_{-\pi}^{\pi} f(x) \sin(nx) dx n\geq 1,

定理 1. 如果 f(x) 在区间 (-\pi, \pi) 上满足 Lipschitz 条件,那么

f(x) =a_{0}+\sum_{n=1}^{\infty} (a_{n} \cos(nx) +b_{n} \sin(nx)).

定理 2. Parseval’s 恒等式.

\frac{1}{\pi} \int_{-\pi}^{\pi} |f(x)|^{2} dx= 2a_{0}^{2}+ \sum_{n=1}^{\infty} (a_{n}^{2}+b_{n}^{2}).

下面我们就来证明下列恒等式:

\sum_{n=1}^{\infty} \frac{1}{(2n-1)^{2}}=\frac{\pi^{2}}{8}

\sum_{n=1}^{\infty} \frac{1}{n^{2}}=\frac{\pi^{2}}{6}

\sum_{n=1}^{\infty} \frac{1}{(2n-1)^{4}}=\frac{\pi^{4}}{96}

\sum_{n=1}^{\infty} \frac{1}{n^{4}}=\frac{\pi^{4}}{90}

证明:

选择在区间 (-\pi, \pi) 上的函数 f(x)=|x|,并且该函数是关于 2\pi 的周期函数。

使用 a_{n}b_{n} 的公式,我们可以得到函数 f(x)=|x| 的 Fourier 级数是

\frac{\pi}{2} + \sum_{n=1}^{\infty} \frac{2((-1)^{n}-1)}{\pi}  \cdot \frac{cos(nx)}{n^{2}}

从定理1, 令 x=0, 可以得到

0= \frac{\pi}{2} + \sum_{n=1}^{\infty} \frac{2((-1)^{n}-1)}{n^{2} \pi}  = \frac{\pi}{2} + \sum_{m=1}^{\infty} \frac{-4}{(2m-1)^{2}\pi}  = \frac{\pi}{2} - \frac{4}{\pi} \sum_{m=1}^{\infty} \frac{1}{(2m-1)^{2}}

因此,\sum_{n=1}^{\infty} \frac{1}{(2n-1)^{2}}=\frac{\pi^{2}}{8} .

假设 S=\sum_{n=1}^{\infty} \frac{1}{n^{2}} , 可以得到

S=\sum_{odd} \frac{1}{n^{2}} + \sum_{even} \frac{1}{n^{2}}  = \frac{\pi^{2}}{8} + \frac{1}{4} S .

因此 S=\frac{\pi^{2}}{6} .

从 Parserval’s 恒等式,我们知道

\frac{2\pi^{2}}{3}= \frac{1}{\pi} \int_{-\pi}^{\pi} x^{2}dx  = 2\cdot (\frac{\pi}{2})^{2} + \sum_{n=1}^{\infty} \frac{4((-1)^{n}-1)^{2}}{\pi^{2}\cdot n^{4}}  = \frac{\pi^{2}}{2} + \sum_{m=1}^{\infty} \frac{16}{\pi^{2} (2m-1)^{4}}

因此 \sum_{n=1}^{\infty} \frac{1}{(2n-1)^{4}} = \frac{\pi^{4}}{96} .

假设 S=\sum_{n=1}^{\infty} \frac{1}{n^{4}} , 得到

S=\sum_{odd} \frac{1}{n^{4}} + \sum_{even} \frac{1}{n^{4}}  = \frac{\pi^{4}}{96} + \frac{1}{16} S

因此, S=\frac{\pi^{4}}{90} .

总结

本篇文章从调和级数的发散性开始,介绍了判断调和级数是否收敛的几种方法。进一步考虑了其他级数的收敛性,并通过 Fourier 级数的方法计算出了部分 Riemann Zeta 函数的取值。

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

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

在 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月

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

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}

存在稳定性的解。

 

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

 

zr9558's Blog