Category Archives: 智能运维

黑盒函数的探索

黑盒函数的定义

在工程上和实际场景中,黑盒函数(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 足够小。

总结

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

 

Advertisements

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

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

 

时间序列的相似性

在文本挖掘中,可以通过 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} 范数的距离,基于相关性的距离,基于周期图表的计算方法,基于模型的计算方法。

 

 

时序数据与事件的关联分析

文章是:”Correlating Events with Time Series for Incident Diagnosis” 是微软在2014年的工作,并且发表在KDD上。

本文提出了一种无监督和统计判别的算法,可以检测出事件(E)与时间序列(S)的关联关系,并且可以检测出时间序列(S)的单调性(上升或者下降)。在这篇文章中,选择的事件有CPU(Memory, Disk)Intensive Program,Query Alert;选择的时间序列有 CPU(Memory)Usage,Disk Transfer Rate。时间序列的特点是它们的值域范围都是[0,1]。

Table1

案例是:时间序列是CPU的Usage,事件是Disk Intensive task和CPU intensive task。

Figure1

关联关系的挖掘分成三个部分:

(1)是否存在关联性(Existence of Dependency):在事件(E)与时间序列(S)之间是否存在关联关系。

(2)关联关系的因果关系(Temporal Order of Dependency):是事件(E)导致了时间序列(S)的变化还是时间序列(S)导致了事件(E)的发生。

(3)关联关系的单调性影响(Monotonic Effect of Dependency):用于判断时间序列(S)是发生了突增或者是突降。

基本概念:

给定一个事件序列(E),事件发生的时间戳是T_{E}=(t_{1},\cdots,t_{n}),这里n表示有n个事件发生。时间序列(S)表示为S=(s_{1},\cdots,s_{m}),这里的m表示时间序列的长度。时间序列的时间戳可以选择一个等差序列,等差用\tau来表示,并且T_{S}=(t(s_{1}),\cdots,t(s_{n})),and t(s_{i}) =t(s_{i-1})+\tau

e_{i}来表示某个事件,\ell_{k}^{rear}(S,e_{i})表示序列S在事件e_{i}之后的长度为k的子序列,\ell_{k}^{front}(S,e_{i})表示序列S在事件e_{i}之前的长度为k的子序列。如果事件E与时间序列S之间存在关联关系,那么

\Gamma^{front}=\{\ell_{k}^{front}(S,e_{i}), i=1,\cdots,n\}

\Gamma^{rear}=\{\ell_{k}^{rear}(S,e_{i}),i=1,\cdots,n\}应该是不一样的。

定义一:如果事件序列E和时间序列S是相关的,并且S->E,当且仅当\Gamma^{front}=\{\ell_{k}^{front}(S,e_{i}), i=1,\cdots,n\}和随机选择的子序列分布不一致。

定义二:如果事件序列E和时间序列S是相关的,并且E->S,当且仅当\Gamma^{rear}=\{\ell_{k}^{rear}(S,e_{i}),i=1,\cdots,n\}和随机选择的子序列分布不一致,并且\Gamma^{front}=\{\ell_{k}^{front}(S,e_{i}), i=1,\cdots,n\}和随机选择的子序列分布一致。

定义三:如果事件序列E和时间序列S是相关的,那么S->E或者E->S

定义四:如果E->S (or S->E),并且时间序列相比E之前是增加了,那么记为E\stackrel{+}{\longrightarrow} S (or S\stackrel{+}{\longrightarrow} E)。如果E->S (or S->E),并且时间序列相比E之前是减少了,那么记为E\stackrel{-}{\longrightarrow} S (or S\stackrel{-}{\longrightarrow} E)。

方法论:

第一步:最邻近算法(类似kNN)(Nearest Neighbor Method)

在计算时间序列之间距离的时候,使用DTW算法或者DTW-D算法会优于L1或者L2算法。

\Gamma^{front}来做例子,\Gamma^{front}=\{\ell_{k}^{front}(S,e_{i}), i=1,\cdots,n\}\Theta =\{\theta_{1},\cdots,\theta_{\tilde{n}}\} 是随机选择的,Z=\Gamma \cup \Theta,可以标记为Z_{1},\cdots,Z_{p},其中p=n+\tilde{n}Z_{i}=\ell_{k}^{front}(S,e_{i}) when 1\leq i\leq nZ_{i}=\theta_{i-n} when n+1\leq i\leq p。可以使用记号A=A_{1}\cup A_{2},其中A_{1}=\Gamma^{front}A_{2}=\Theta=\{\theta_{1},\cdots,\theta_{\tilde{n}}\}是随机选择的。

对于集合 Ax\in A 而言,NN_{r}(x,A) 表示A-\{x\}中距离x最近的第r个元素,对于两个不相交的集合A_{1}A_{2},可以定义方程:

I_{r}(x,A_{1},A_{2})=1 when x\in A_{i} \&\& NN_{r}(x,A)\in A_{i},

I_{r}(x,A_{1},A_{2})=0 when otherwise.

该方程I_{r}(x,A_{1},A_{2})表示x与x的第r个最近的邻居是否在同一个子集内。

定义

T_{r,p}=\frac{1}{pr}\sum_{i=1}^{p}\sum_{j=1}^{r}I_{j}(x_{i},A_{1},A_{2}),

在这里p=n+\tilde{n}表示样本的总个数,x_{i}表示集合A的第i个元素。从直觉上讲,如果T_{r,p}小,则说明两类samples A_{1},A_{2}混合得非常好,表示无异常情况;如果T_{r,p}大,则说明两类samples A_{1},A_{2}有区分度,很多元素与它的邻居集中在某个子集里面,说明 A_{1} 这个集合与 A_{2} 有区分度。

根据文献里面的观点,当p足够大的时候,(pr)^{\frac{1}{2}}(T_{r,p}-\mu_{r})/\sigma_{r}遵循标准Gauss分布,其参数是\mu_{r}=(\lambda_{1})^{2}+(\lambda_{2})^{2}, \sigma_{r}^{2}=\lambda_{1}\lambda_{2}+4\lambda_{1}^{2}\lambda_{2}^{2},

\lambda_{1}=n/p=n/(n+\tilde{n}), \lambda_{2}=\tilde{n}/(n+\tilde{n})

根据传统的Gauss分布Test方法,\Gamma^{front}\Theta有显著的不同,当(pr)^{\frac{1}{2}}(T_{r,p}-\mu_{r})/\sigma_{r}^{2}>\alpha,在这里,参数可以按照以下标准设置:

\alpha = 1.96 for P=0.025

\alpha = 2.58 for P=0.001

如果\Gamma^{front}\Theta存在显著性偏差,那么说明\Gamma^{front}应该返回异常的标识。类似的,如果使用\Gamma^{rear}并且它与\Theta存在显著性偏差,那么说明\Gamma^{rear}应该返回异常的标识。

 

第二步:关联顺序的挖掘(Mining Existence and Temporal Order)

Figure3

如果前面的子序列\Gamma^{front}与随机选择的子序列\Theta有显著偏差,那么说明时序的变化导致了事件的发生,S\rightarrow E

如果后面的子序列\Gamma^{rear}与随机选择的子序列\Theta有显著偏差,那么说明事件导致了时序的变化,E\rightarrow S

在Figure 3中,CPU Intensive Program 导致了 CPU Usage,并且 CPU Usage 导致了 SQL Query Alert。

 

第三步:单调性的影响类型(Mining Effect Type)

现在需要判断时间序列是突增还是突降了,需要引入t_{score}的概念。

对于\Gamma^{front}=\{\ell_{k}^{front}(S,e_{i}), i=1,\cdots,n\}\Gamma^{rear}=\{\ell_{k}^{rear}(S,e_{i}), i=1,\cdots,n\}而言,其中n是E中的事件个数。t_{score}就可以定义为:

t_{score}=\frac{\mu_{\Gamma^{front}} - \mu_{\Gamma^{rear}}}{\sqrt{\frac{\sigma_{\Gamma^{front}}^{2}+\sigma_{\Gamma^{rear}}^{2}}{n}}}.

那么,如果t_{score}>\alpha,可以得到 E\stackrel{-}{\longrightarrow}S 或者 S\stackrel{-}{\longrightarrow} E;如果t_{score}<-\alpha,可以得到 E\stackrel{+}{\longrightarrow}S 或者 S\stackrel{+}{\longrightarrow} E

其中参数可以设置为:

\alpha = 1.96 for P=0.025

\alpha = 2.58 for P=0.001

 

算法综述:

algorithm

其中,5,6行是为了计算相关性,D_{r} 是 True 表示 \Gamma^{rear} 有异常,否则表示正常;D_{f} 是 True 表示 \Gamma^{front} 有异常,否则表示正常。

7-13行是 E\rightarrow S 的情形,因为\Gamma^{rear} 异常,同时 \Gamma^{front} 正常,说明事件导致了时间序列的变化。7-13行是为了计算 t_{score} 的范围,判断是显著的提升还是下降。

14-20行是 S\rightarrow E 的情形,因为\Gamma^{front} 异常,就导致了事件的发生。14-20行是为了计算 t_{score} 的范围,判断是显著的提升还是下降。

参数:

时间序列的长度 k 可以设置为第一次达到顶峰的长度,

最邻近的元素个数 r=\ln(p),其中p是样本的总个数。

Figure5

Figure6

其他算法:

(1)Pearson Correlation

(2)J-Measure Correlation

 

 

Opprentice: Towards Practical and Automatic Anomaly Detection Through Machine Learning

本文是运维系统智能化的一次探索工作,论文的作者是清华大学的裴丹教授,论文的题目是《Opprentice: Towards Practical and Automatic Anomaly Detection Through Machine Learning》。目的是基于机器学习的 KPI(Key Performance Indicator)的自动化异常检测。

标题 Opprentice 来源于(Operator’s Apprentice),意思就是运维人员的学徒。本文通过运维人员的业务经验来进行异常数据的标注工作,使用时间序列的各种算法来提取特征,并且使用有监督学习模型(例如 Random Forest,GBDT,XgBoost 等)模型来实现离线训练和上线预测的功能。本文提到系统 Opprentice 使用了一个多月的历史数据进行分析和预测,已经可以做到准确率>=0.66,覆盖率>=0.66 的效果。

1. Opprentice的介绍

系统遇到的挑战:

Definition Challenges: it is difficult to precisely define anomalies in reality.(在现实环境下很难精确的给出异常的定义)

Detector Challenges: In order to provide a reasonable detection accuracy, selecting the most suitable detector requires both the algorithm expertise and the domain knowledge about the given service KPI (Key Performance Indicators). To address the definition challenge and the detector challenge, we advocate for using supervised machine learning techniques. (使用有监督学习的方法来解决这个问题)

该系统的优势:

(i) Opprentice is the first detection framework to apply machine learning to acquiring realistic anomaly definitions and automatically combining and tuning diverse detectors to satisfy operators’ accuracy preference.

(ii) Opprentice addresses a few challenges in applying machine learning to such a problem: labeling overhead, infrequent anomalies, class imbalance, and irrelevant and redundant features.

(iii) Opprentice can automatically satisfy or approximate a reasonable accuracy preference (recall>=0.66 & precision>=0.66). (准确率和覆盖率的效果)

2. 背景描述:

KPIs and KPI Anomalies:

KPIs: The KPI data are the time series data with the format of (time stamp, value). In this paper, Opprentice pays attention to three kinds of KPIs: the search page view (PV), which is the number of successfully served queries; The number of slow responses of search data centers (#SR); The 80th percentile of search response time (SRT).

Anomalies: KPI time series data can also present several unexpected patterns (e.g. jitters, slow ramp ups, sudden spikes and dips) in different severity levels, such as a sudden drop by 20% or 50%.

OpprenticeFigure1

问题和目标:

覆盖率(recall):# of true anomalous points detected / # of the anomalous points

准确率(precision):# of true anomalous points detected / # of anomalous points detected

1-FDR(false discovery rate):# of false anomalous points detected / # of anomalous points detected = 1 – precision

The quantitative goal of opprentice is precision>=0.66 and recall>=0.66.

The qualitative goal of opprentice is automatic enough so that the operators would not be involved in selecting and combining suitable detectors, or tuning them.

3. Opprentice Overview: (Opprentice系统的概况)

OpprenticeFigure2

(i) Opprentice approaches the above problem through supervised machine learning.

(ii) Features of the data are the results of the detectors.(Basic Detectors 来计算出特征)

(iii) The labels of the data are from operators’ experience.(人工打标签)

(iv) Addressing Challenges in Machine Learning: (机器学习遇到的挑战)

(1) Label Overhead: Opprentice has a dedicated labeling tool with a simple and convenient interaction interface. (标签的获取)

(2) Incomplete Anomaly Cases:(异常情况的不完全信息)

(3) Class Imbalance Problem: (正负样本比例不均衡)

(4) Irrelevant and Redundant Features:(无关和多余的特征)

4. Opprentice’s Design:

Architecture: Operators label the data and numerous detectors functions are feature extractors for the data.

OpprenticeFigure3

Label Tool:

人工使用鼠标和软件进行标注工作

OpprenticeFigure4

Detectors:

(i) Detectors As Feature Extractors: (Detector用来提取特征)

Here for each parameter detector, we sample their parameters so that we can obtain several fixed detectors, and a detector with specific sampled parameters a (detector) configuration. Thus a configuration acts as a feature extractor:

data point + configuration (detector + sample parameters) -> feature,

(ii) Choosing Detectors: (Detector的选择,目前有14种较为常见的)

Opprentice can find suitable ones from broadly selected detectors, and achieve a relatively high accuracy. Here, we implement 14 widely-used detectors in Opprentice.

Opprentice has 14 widely-used detectors:

OpprenticeTable3

Diff“: it simply measures anomaly severity using the differences between the current point and the point of last slot, the point of last day, and the point of last week.

MA of diff“: it measures severity using the moving average of the difference between current point and the point of last slot.

The other 12 detectors come from previous literature. Among these detectors, there are two variants of detectors using MAD (Median Absolute Deviation) around the median, instead of the standard deviation around the mean, to measure anomaly severity.

(iii) Sampling Parameters: (Detector的参数选择方法,一种是扫描参数空间,另外一种是选择最佳的参数)

Two methods to sample the parameters of detectors.

(1) The first one is to sweep the parameter space. For example, in EWMA, we can choose \alpha \in \{0.1,0.3,0.5,0.7,0.9\} to obtain 5 typical features from EWMA; Holt-Winters has three [0,1] valued parameters \alpha,\beta,\gamma. To choose \alpha,\beta,\gamma \in \{0.2,0.4,0.6,0.8\}, we have 4^3 features; In ARIMA, we can estimate their “best” parameters from the data, and generate only one set of parameters, or one configuration for each detector.

Supervised Machine Learning Models:

Decision Trees, logistic regression, linear support vector machines (SVMs), and naive Bayes. 下面是决策树(Decision Tree)的一个简单例子。

OpprenticeFigure5

Random Forest is an ensemble classifier using many decision trees. It main principle is that a group of weak learners (e.g. individual decision trees) can together form a strong learner. To grow different trees, a random forest adds some elements or randomness. First, each tree is trained on subsets sampled from the original training set. Second, instead of evaluating all the features at each level, the trees only consider a random subset of the features each time. The random forest combines those trees by majority vote. The above properties of randomness and ensemble make random forest more robust to noises and perform better when faced with irrelevant and redundant features than decisions trees.

Configuring cThlds: (阈值的计算和预估)

(i) methods to select proper cThlds: offline part

OpprenticeFigure6

We need to figure cThlds rather than using the default one (e.g. 0.5) for two reasons.

(1) First, when faced with imbalanced data (anomalous data points are much less frequent than normal ones in data sets), machine learning algorithems typically fail to identify the anomalies (low recall) if using the default cThlds (e.g. 0.5).

(2) Second, operators have their own preference regarding the precision and recall of anomaly detection.

The metric to evaluate the precision and recall are:

(1) F-Score: F-Score = 2*precision*recall/(precision+recall).

(2) SD(1,1): it selects the point with the shortest Euclidean distance to the upper right corner where the precision and the recall are both perfect.

(3) PC-Score: (本文中采用这种评估指标来选择合适的阈值)

If r>=R and p>=P, then PC-Score(r,p)=2*r*p/(r+p) + 1; else PC-Score(r,p)=2*r*p/(r+p). Here, R and P are from the operators’ preference “recall>=R and precision>=P”. Since the F-Score is no more than 1, then we can choose the cThld corresponding to the point with the largest PC-Score.

(ii) EWMA Based cThld Prediction: (基于EWMA方法的阈值预估算法)

OpprenticeFigure7

In online detection, we need to predict cThlds for detecting future data.

Use EWMA to predict the cThld of the i-th week ( or the i-th test set) based on the historical best cThlds. Specially, EWMA works as follows:

If i=1, then cThld_{i}^{p}=cThld_{1}^{p}= 5-fold prediction

Else i>1, then cThld_{i}^{p}=\alpha\cdot cThld_{i-1}^{b}+(1-\alpha)\cdot cThld_{i-1}^{p}, where cThld_{i-1}^{b} is the best cThld of the (i-1)-th week. cThld_{i}^{p} is the predicted cThld of the i-th week, and also the one used for detecting the i-th week data. \alpha\in [0,1] is the smoothing constant.

For the first week, we use 5-fold cross-validation to initialize cThld_{1}^{p}. As \alpha increases, EWMA gives the recent best cThlds more influences in the prediction. We use \alpha=0.8 in this paper.

5. Evaluation(系统评估)

在 Opprentice 系统中,红色表示 Opprentice 系统的方法,黑色表示其他额外的方法。

OpprenticeFigure8

Opprentice has 14 detectors with about 9500 lines of Python, R and C++ code. The machine learning block is based on the scikit-learn library.

Random Forest is better than decision trees, logistic regression, linear support vector machines (SVMs), and naive Bayes.