Monthly Archives: April 2018

智能运维

黑盒函数的探索

Leave a comment

黑盒函数的定义

在工程上和实际场景中，黑盒函数（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)=0 的 x 的解。特别的，对于二次多项式而言，也就是 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})<0 和 f(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 足够小。

总结

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

 

Advertisement
Complex Analysis

Riemann Zeta 函数（二）

Leave a comment

在上一篇文章里面，我们已经给出了 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 函数没有零点。

 

Mathematical Analysis

从调和级数到 RIEMANN ZETA 函数（一）

Leave a comment

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 函数的取值。

Computer Science

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

Leave a comment

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

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

v2-58d234f9633710b627fd7199d233b7d1_hd

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 也产生了浓厚的兴趣，撰写过两篇文章“强化学习与泛函分析”，“深度学习与强化学习”。

v2-480b992cb49efc66a880eb4514edaae0_hd

2017 年：再整旗鼓

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

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

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

v2-4fcc9c421892c91f9ce9cdd3589a508e_hd

正式接触到运维项目是 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 时间序列进行异常检测。这一类方法笔者总结过，那就是“时间序列简介（一）”，最终我们做到了“百万条曲线一人挑”，成功去掉了制定阈值的业务效果。

v2-fadf803defdd739088284af0c7ec9039_hd

2018年：走向未来

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

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

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

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

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

v2-b70865721c58f2cb1a53c2ccc2a2f824_hd

张戎

2018年2月

数据挖掘与机器学习

启发式优化算法

Leave a comment

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

启发式优化算法的背景

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年