符号计算中的深度学习方法

符号计算

符号计算一直是计算数学的重要领域之一。在开源领域,Python 的 SymPy 就可以支持符号计算。在商业化领域,Maple,Matlab,Mathematica 都能够进行符号计算。它们不仅能够做简单的实数和复数加减乘除,还能够支持数学分析,线性代数,甚至各种各样的大学数学课程。

随着人工智能的进一步发展,深度学习不仅在图像识别,自然语言处理方向上发挥着自身的价值,还在各种各样的领域展示着自己的实用性。在 2019 年底,facebook 两位研究员在 arxiv 上挂出了一篇文章《Deep Learning for Symbolic Mathematics》,在符号计算方向上引入了深度学习的工具。

要想了解符号运算,就要先知道在计算机中,是怎么对数学公式进行表示的。较为常见的表达式形如:

  • 2 + 3 * (5 + 2)
  • 3x^{2}+\cos(2x)+1
  • \frac{\partial^{2}\psi}{\partial x^{2}} - \frac{1}{v^{2}}\frac{\partial^{2}\psi}{\partial t^{2}}

在这里,数学表达式通常都会被表示成树的结构,其中树的内部节点是由算子(operator),函数(function)组成的,叶子节点由数字,变量,函数等组成。例如:

dlsm_figure1.png
图 1

图 1 的三幅图分别对应着上面的三个数学表达式。

在 Python 的 SymPy 工具中,同样可以对数学公式进行展示。其表示方法就是用 sympy.srepr

>>> import sympy
>>> x, y = sympy.symbols("x y")
>>> expr = sympy.sin(x+y) + x**2 + 1/y - 10
>>> sympy.srepr(expr)
"Add(Pow(Symbol('x'), Integer(2)), sin(Add(Symbol('x'), Symbol('y'))), Integer(-10), Pow(Symbol('y'), Integer(-1)))"
>>> expr = sympy.sin(x*y)/2 - x**2 + 1/y
>>> sympy.srepr(expr)
"Add(Mul(Integer(-1), Pow(Symbol('x'), Integer(2))), Mul(Rational(1, 2), sin(Mul(Symbol('x'), Symbol('y')))), Pow(Symbol('y'), Integer(-1)))"
dlsm_figure2.png
图 2

SymPy 的 srepr 函数的输出用树状结构来表示就是形如图 2 这种格式。叶子节点要么是 x,y 这种变量,要么是 -1 和 2 这种整数。对于一元函数而言,例如 sin 函数,就是对应唯一的一个叶子。对于二元函数而言,例如 pow,mul,add,则是对应两个叶子节点。

论文方案

在 Deep Learning for Symbolic Mathematics 这篇论文中,作者们的大致思路是分成以下几步的:

  1. 生成数据;
  2. 训练模型;
  3. 预测结果;

第一步生成数据是为了让深度学习模型是大量的已知样本来做训练;第二步训练模型是用第一步的数据把端到端的深度学习模型进行训练;第三步预测结果是给一个函数或者一个微分方程,使用已经训练好的模型来预测结果,对预测出来的多个结果进行排序,选择最有可能的那个结果作为符号计算的值。

众所周知,深度学习的训练是依赖大量的样本数据的,那么要想用深度学习来解决符号计算的问题,就要解决样本少的问题。在这篇论文中,作者们把精力投入了三个领域,分别是:

  1. 函数的积分;
  2. 一阶常微分方程;
  3. 二阶常微分方程。

在生成数据之前,作者们对数据的范围进行了必要的限制:

  1. 数学表达式最多拥有 15 个内部节点;
  2. L = 11 表示叶子节点的值只有 11 个,分别是变量 x\{-5,-4,-3,-2,-1,1,2,3,4,5\}
  3. p_{1} = 15 表示一元计算只有 15 个,分别是 \exp, \log, \sqrt, \sin, \cos, \tan, \arcsin, \arccos, \arctan, sinh, cosh, tanh, arcsinh, arccosh, arctanh
  4. p_{2} = 4 表示二元计算只有四个,分别是 +, -, *, /;

意思就是在这个有限的范围内去生成深度学习所需要的数据集。

积分数据的生成

在微积分里面,积分指的是求导的逆运算,那么形如 (f',f) 的表达式就可以作为深度学习的积分训练数据。生成积分的话其实有多种方法:

第一种方法:前向生成(Forward Generation,简写为 FWD)。主要思路就是在以上的数据范围内随机生成各种各样的方程 f,然后使用 SymPy 或者 Mathematica 等工具来计算函数 f 的积分 F,那么 (f,F) 就可以作为一个训练集。当然,有的时候函数 f 的积分是无法计算出来的,那么这种计算表达式就需要进行放弃,就不能放入训练集。

第二种方法:反向生成(Backward Generation,简写为 BWD)。由于积分是求导的逆运算,可以在以上的数据范围内随机生成各种各样的方程 f,然后计算它们的导数 f',于是 (f',f) 就可以放入积分数据的训练集。

第三种方法:分部积分(Backward generation with integration by parts,简写为 IBP)。根据分部积分的公式,如果 F'=f, G'=g,那么 \int F(x)g(x)dx=F(x)G(x) - \int f(x)G(x)dx。对于两个随机生成的函数 F, G,可以计算出它们的导数 f, g。如果 fG 在训练集合里面,那么就把 Fg 的积分计算出来放入训练集合;反之,如果 Fg 在训练集合里面,那么就把 fG 的积分计算出来放入训练集合。如果 fGFg 都没有在训练集合,并且都无法积分出来,那么就放弃该条数据。

三种方法的比较可以参见表 1,从表 1 可以看出通过分部积分(Integration by parts)所获得的样本表达式相对于前向生成和反向生成都比较长。

dlsm_table1.png
表 1

一阶常微分方程的生成

一阶常微分方程只有该函数的一阶导数,因此在构造函数的时候,作者们用了一个技巧。在随机生成表达式的时候,叶子节点的元素只从 \{x, -5,-4,-3,-2,-1,1,2,3,4,5\} 选择一个,于是随机把其中的一个整数换成变量 c。例如:在 x\log(2/x) 中就把 2 换成 c,于是得到了一个二元函数 f(x,c) = x\log(c/x)。那么就执行以下步骤:

  1. 生成二元函数 f(x,c) = x\log(c/x)
  2. 求解 c,得到 c = xe^{f(x)/x}
  3. x 进行求导得到 0 = e^{f(x)/x}(1+f'(x)-f(x)/x) = 0
  4. 简化后得到 x+xf'(x) - f(x) =0,也就是 x+xy'-y=0

此时,一阶常微分方程的训练数据就是 (x+xy'-y=0, x\log(c/x))

二阶常微分方程不仅可能由该函数的一阶导数,还必须有二阶导数,那么此时就要求导两次:

  1. 生成三元函数 f(x,c_{1},c_{2}) = c_{1}e^{x}+c_{2}e^{-x}
  2. 求解 c_{2} 得到 c_{2} = f(x,c_{1},c_{2}) e^{x}- c_{1}e^{2x}
  3. x 进行求导得到 0 = e^{x}(\partial f(x,c_{1},c_{2})/\partial x + f(x,c_{1},c_{2})) - 2c_{1}e^{2x} = 0
  4. 求解 c_{1} 得到 c_{1} = e^{-x}(\partial f(x,c_{1},c_{2})/\partial x+f(x))/2
  5. x 进行求导得到 0 = e^{-x}(\partial^{2} f(x,c_{1},c_{2})/\partial x^{2} - f(x,c_{1},c_{2}))=0
  6. 简化后得到 \partial^{2} f(x,c_{1},c_{2})/\partial x^{2} - f(x,c_{1},c_{2})=0,也就是 y''-y=0

那么此时的二阶微分方程的训练数据就是 (y''-y=0, c_{1}e^{x}+c_{2}e^{-x})

需要注意的事情就是,在生成数据的时候,一旦无法求出 c_{1}, c_{2} 的表达式,该数据就需要放弃,重新生成新的数据。

数据处理

  1. 数学表达式的简化(expression simplification):例如 x+1+1 可以简化成 x+2\sin^{2}(x)+\cos^{2}(x) 可以简化成 1。
  2. 参数的简化(coefficient simplification):例如 \log(x^{2}) + c\log(x) 可以简化成 c\log(x)
  3. 无效表达式的过滤(invalid expression filter):例如 \sqrt{2}, \log(0) 等。

树状结构的表达式,是使用前缀表达式来写成一个句子的。例如 2+3 就写成 + 2 3,2 + x 就写成 + 2 x。

模型训练

在这里,作者们用了 Transformer 模型,8 attention heads,6 layers,512 个维度,并且发现用更复杂的模型并没有提升其效果。在预测的时候,使用不同的 Beam Size,其准确率是不一样的,在 Beam Size = 50 的时候,效果比较好。参见表 2。

dlsm_table2.png
表 2

在与其他数学软件的比较中, 作者限制了 Mathematica 的运行时间为 30s。

dlsm_table3.png

并且举出了几个现有模型能够计算出来,Mathematica 无法计算的例子。

dlsm_table4.png

在求解的等价性方面,可以根据 Score 逆序排列,然后这个例子的 Top10 都是等价的。

dlsm_table5.png

在使用 SymPy 的过程中,可以获得各种各样的积分表达式如下:

dlsm_table7.png

结论

符号计算在 1960 年代末就已经在研究了,有诸多的符号计算软件,例如 Matlab,Mathematica,Maple,PARI,SAGE 等。在这篇论文中,作者们使用标准的 seq2seq 模型来对生成的数据进行训练,然后获得积分,一阶常微分方程,二阶常微分方程的解。在传统符号计算的基础上,开拓了一种新的思路。

用 Python 来研究数学 — SymPy 符号工具包介绍

SymPy 的简单介绍

SymPy 是一个符号计算的 Python 库,完全由 Python 写成,为许多数值分析,符号计算提供了重要的工具。SymPy 的第一个版本于 2007 年开源,并且经历了十几个版本的迭代,在 2019 年已经基于修正的 BSD 许可证开源了 1.4 版本。SymPy 的开源地址和官方网站分别是:

  1. GitHub 链接:https://github.com/sympy/sympy
  2. SymPy 官方网站:https://www.sympy.org/en/index.html
sympy_logo
SymPy 的 logo

SymPy 的 1.4 版本文档中,可以看出,SymPy 可以支持很多初等数学,高等数学,甚至研究生数学的符号计算。在初等数学和高等数学中,SymPy 可以支持的内容包括但不限于:

  1. 基础计算(Basic Operations);
  2. 公式简化(Simplification);
  3. 微积分(Calculus);
  4. 解方程(Solver);
  5. 矩阵(Matrices);
  6. 几何(geometry);
  7. 级数(Series);

在更多的数学领域中,SymPy 可以支持的内容包括但不限于:

  1. 范畴论(Category Theory);
  2. 微分几何(Differential Geometry);
  3. 常微分方程(ODE);
  4. 偏微分方程(PDE);
  5. 傅立叶变换(Fourier Transform);
  6. 集合论(Set Theory);
  7. 逻辑计算(Logic Theory)。
sympy_tutorial
SymPy 的教学目录

SymPy 的工具库介绍

SymPy 的基础计算

在数学中,基础的计算包括实数和复数的加减乘除,那么就需要在程序中描述出实数与复数。著名的欧拉公式

e^{i\pi}+1 = 0

正好用到了数学中最常见的五个实数。在 SymPy 里面,e, i, \pi, \infty 是用以下符号来表示的:其中 sympy.exp() 表示以 e 为底的函数。

sympy.exp(1), sympy.I, sympy.pi, sympy.oo

而想要计算欧拉公式的话,只需要输入下面的公式即可:

>>> sympy.exp(sympy.I * sympy.pi) + 1
0

如果需要看 e, \pi 的小数值,可以使用 evalf() 函数,其中 evalf() 函数里面的值表示有效数字的位数。例如下面就是精确到 10 位有效数字。当然,也可以不输入。

>>> sympy.E.evalf(10)
2.718281828
>>> sympy.E.evalf()
2.71828182845905
>>> sympy.pi.evalf(10)
3.141592654
>>> sympy.pi.evalf()
3.14159265358979

除此之外,如果需要查看某个实数的有效数字,也是类似操作的:

>>> expr = sympy.sqrt(8)
>>> expr.evalf()
2.82842712474619

而对于实数的加减乘除,则可以如下操作:

>>> x, y= sympy.symbols("x y")
>>> x + y
x + y
>>> x - y
x - y
>>> x * y
x*y
>>> x / y
x/y

而对于复数的加减乘除,则是类似的操作,令两个复数分别是 z_{1} = x_{1} + i y_{1}z_{2} = x_{2} + i y_{2}

>>> x1, y1, x2, y2 = sympy.symbols("x1 y1 x2 y2")
>>> z1 = x1 + y1 * sympy.I
x1 + I*y1
>>>  z2 = x2 + y2 * sympy.I
x2 + I*y2
>>> z1 + z2
x1 + x2 + I*y1 + I*y2
>>> z1 - z2
x1 - x2 + I*y1 - I*y2
>>> z1 * z2
(x1 + I*y1)*(x2 + I*y2)
>>> z1 / z2
(x1 + I*y1)/(x2 + I*y2)

对于多项式而言,有的时候我们希望将其展开,有的时候则需要将其合并,最终将其简化成最简单的形式。

>>> sympy.expand((x+1)**2)
x**2 + 2*x + 1
>>> sympy.expand((x+1)**5)
x**5 + 5*x**4 + 10*x**3 + 10*x**2 + 5*x + 1
>>> sympy.factor(x**3+1)
(x + 1)*(x**2 - x + 1)
>>> sympy.factor(x**2+3*x+2)
(x + 1)*(x + 2)
>>> sympy.simplify(x**2 + x + 1 - x)
x**2 + 1
>>> sympy.simplify(sympy.sin(x)**2 + sympy.cos(x)**2)
1

在多变量的场景下,SymPy 也可以对其中的某个变量合并同类项,同时还可以计算某个变量的某个次数所对应的系数是多少,例如:

>>> expr = x*y + x - 3 + 2*x**2 - x**2 + x**3 + y**2 + x**2*y**2
>>> sympy.collect(expr,x)
x**3 + x**2*(y**2 + 1) + x*(y + 1) + y**2 - 3
>>> sympy.collect(expr,y)
x**3 + x**2 + x*y + x + y**2*(x**2 + 1) - 3
>>> expr.coeff(x, 2)
y**2 + 1
>>> expr.coeff(y, 1)
x

有理函数形如 f(x) = p(x)/q(x),其中 p(x)q(x) 都是多项式。一般情况下,我们希望对有理函数进行简化,合并或者分解的数学计算。

在需要合并的情形下,如果想把有理函数处理成标准格式 p(x)/q(x) 并且去除公因子,那么可以使用 cancel 函数。另一个类似的就是 together 函数,但是不同之处在于 cancel 会消除公因子,together 不会消除公因子。例如:

expr = \frac{x^{2}+3x+2}{x^{2}+x}

>>> expr = (x**2 + 3*x + 2)/(x**2 + x)
>>> sympy.cancel(expr)
(x + 2)/x
>>> sympy.together(expr)
(x**2 + 3*x + 2)/(x*(x + 1))

除了合并和消除公因子之外,有的时候还希望对分子和分母进行因式分解,例如:

expr = (x**2 + 3*x + 2)/(x**2 + x)
>>> sympy.factor(expr)
(x + 2)/x
>>> expr = (x**3 + 3*x**2 + 2*x)/(x**5+x)
>>> sympy.factor(expr)
(x + 1)*(x + 2)/(x**4 + 1)
>>> expr = x**2 + (2*x+1)/(x**3+1)
>>> sympy.factor(expr)
(x**5 + x**2 + 2*x + 1)/((x + 1)*(x**2 - x + 1))

合并的反面就是部分分式展开(Partial Fraction Decomposition),它是把有理函数分解成多个次数较低的有理函数和的形式。这里需要用 apart 函数:

>>> expr = (x**4 + 3*x**2 + 2*x)/(x**2+x)
>>> sympy.apart(expr)
x**2 - x + 4 - 2/(x + 1)
>>> expr = (x**5 + 1)/(x**3+1)
>>> sympy.apart(expr)
x**2 - (x - 1)/(x**2 - x + 1)

在 SymPy 里面,同样支持各种各样的三角函数,例如:三角函数的简化函数 trigsimp,三角函数的展开 expand_trig,

>>> expr = sympy.sin(x)**2 + sympy.cos(x)**2
>>> sympy.trigsimp(expr)
1
>>> sympy.expand_trig(sympy.sin(x+y))
sin(x)*cos(y) + sin(y)*cos(x)
>>> sympy.expand_trig(sympy.cos(x+y))
-sin(x)*sin(y) + cos(x)*cos(y)
>>> sympy.trigsimp(sympy.sin(x)*sympy.cos(y) + sympy.sin(y)*sympy.cos(x))
sin(x + y)
>>> sympy.trigsimp(-sympy.sin(x)*sympy.sin(y) + sympy.cos(x)*sympy.cos(y))
cos(x + y)

同样的,在乘幂上面,同样有简化函数 powsimp,效果与之前提到的 simplify 一样。除此之外,还可以根据底数来做合并,即分别使用 expand_power_exp 函数与 expand_power_base 函数。

>>> sympy.powsimp(x**z*y**z*x**z)
x**(2*z)*y**z
>>> sympy.simplify(x**z*y**z*x**z)
x**(2*z)*y**z
>>> sympy.expand_power_exp(x**(y + z))
x**y*x**z
>>> sympy.expand_power_base(x**(y + z))
x**(y + z)

作为指数的反函数对数,sympy 也是有着类似的展开合并函数,expand_log,logcombine 承担了这样的角色。

\ln(xy) = \ln(x) + \ln(y)

\ln(x/y) = \ln(x) - \ln(y)

>>> sympy.expand_log(sympy.log(x*y), force=True)
log(x) + log(y)
>>> sympy.expand_log(sympy.log(x/y), force=True)
log(x) - log(y)

 

SymPy 的微积分工具

下面,我们会从一个最基本的函数 f(x) = 1/x 出发,来介绍 SymPy 的各种函数使用方法。如果想进行变量替换,例如把 x 变成 y,那么可以使用 substitution 方法。除此之外,有的时候也希望能够得到函数 f 在某个点的取值,例如 f(1) = 1/1 = 1,那么可以把参数换成 1 即可得到函数的取值。例如,

>>> import sympy
>>> x = sympy.Symbol("x")
>>> f = 1 / x
1/x
>>> y = sympy.Symbol("y")
>>> f = f.subs(x,y)
1/y
>>> f = f.subs(y,1)
1

在微积分里面,最常见的概念就是极限,SymPy 里面的极限函数就是 limit。使用方法如下:

>>> f = 1/x
>>> sympy.limit(f,x,0)
oo
>>> sympy.limit(f,x,2)
1/2
>>> sympy.limit(f,x,sympy.oo)
0
>>> g = x * sympy.log(x)
>>> sympy.limit(g,x,0)
0

对于函数 f(x) = 1/x 而言,它的导数计算函数是 diff,n 阶导数也可以用这个函数算。

>>> f = 1/x
>>> sympy.diff(f,x)
-1/x**2
>>> sympy.diff(f,x,2)
2/x**3
>>> sympy.diff(f,x,3)
-6/x**4
>>> sympy.diff(f,x,4)
24/x**5

提到 n 阶导数,就必须要提一下 Taylor Series 了。对于常见函数的 Taylor Series,SymPy 也是有非常简便的方法,那就是 series 函数。其参数包括 expr, x, x0, n, dir,分别对应着表达式,函数的自变量,Taylor Series 的中心点,n 表示阶数,dir 表示方向,包括”+-“,”-“,”+”,分别表示 x\rightarrow x0, x\rightarrow x0^{-}, x\rightarrow x0^{+}

sympy.series.series.series(exprx=Nonex0=0n=6dir='+') >>> g = sympy.cos(x) >>> sympy.series(g, x) 1 - x**2/2 + x**4/24 + O(x**6) >>> sympy.series(g, x, x0=1, n=10) cos(1) - (x - 1)*sin(1) - (x - 1)**2*cos(1)/2 + (x - 1)**3*sin(1)/6 + (x - 1)**4*cos(1)/24 - (x - 1)**5*sin(1)/120 - (x - 1)**6*cos(1)/720 + (x - 1)**7*sin(1)/5040 + (x - 1)**8*cos(1)/40320 - (x - 1)**9*sin(1)/362880 + O((x - 1)**10, (x, 1))

积分的计算函数是 integrate,包括定积分与不定积分:

\int\frac{1}{x}dx = \ln(x)+C

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

>>> f = 1/x
>>> sympy.integrate(f,x)
log(x)
>>> sympy.integrate(f, (x,1,2))
log(2)

对于广义积分而言,就需要用到 \infty 这个概念了,但是在 SymPy 里面的写法还是一样的。

\int_{-\infty}^{0}e^{-x^{2}}dx=\frac{\sqrt{\pi}}{2}

\int_{0}^{+\infty}e^{-x}dx = 1

\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}e^{-x^{2}-y^{2}}dxdy = \pi

>>> g = sympy.exp(-x**2)
>>> sympy.integrate(g, (x,-sympy.oo,0))
sqrt(pi)/2
>>> g = sympy.exp(-x)
>>> sympy.integrate(g, (x, 0, sympy.oo))
1
>>> h = sympy.exp(-x**2 - y**2)
>>> sympy.integrate(h, (x,-sympy.oo, sympy.oo), (y, -sympy.oo, sympy.oo))
pi

 

SymPy 的方程工具

在初等数学中,通常都存在一元一次方程,一元二次方程等,并且在不同的域上有着不同的解。SymPy 里面的相应函数就是 solveset,根据定义域的不同,可以获得完全不同的解。

\{x\in\mathbb{R}: x^{3}-1=0\}

\{x\in\mathbb{C}:x^{3}-1=0\}

\{x\in\mathbb{R}:e^{x}-x=0\}

\{x\in\mathbb{R}:e^{x}-1=0\}

\{x\in\mathbb{C}:e^{x}-1=0\}

>>> sympy.solveset(sympy.Eq(x**3,1), x, domain=sympy.S.Reals)
{1}
>>> sympy.solveset(sympy.Eq(x**3,1), x, domain=sympy.S.Complexes)
{1, -1/2 - sqrt(3)*I/2, -1/2 + sqrt(3)*I/2}
>>> sympy.solveset(sympy.Eq(x**3 - 1,0), x, domain=sympy.S.Reals)
{1}
>>> sympy.solveset(sympy.Eq(x**3 - 1,0), x, domain=sympy.S.Complexes)
{1, -1/2 - sqrt(3)*I/2, -1/2 + sqrt(3)*I/2}
>>> sympy.solveset(sympy.exp(x),x)
EmptySet()
>>> sympy.solveset(sympy.exp(x)-1,x,domain=sympy.S.Reals)
{0}
>>> sympy.solveset(sympy.exp(x)-1,x,domain=sympy.S.Complexes)
ImageSet(Lambda(_n, 2*_n*I*pi), Integers)

在这里,Lambda 表示的是数学公式,第一个是自变量,第二个是函数,最后是自变量的定义域。

在线性代数中,最常见的还是多元一次方程组,那么解法是一样的:

\begin{cases}x+y-10=0 \\ x-y-2=0\end{cases}

>>> sympy.solve([x+y-10, x-y-2], [x,y])
{x: 6, y: 4}

对于三角函数,也是类似的写法:

\begin{cases} \sin(x-y)=0 \\ \cos(x+y)=0 \end{cases}

>>> sympy.solve([sympy.sin(x-y), sympy.cos(x+y)], [x,y])
[(-pi/4, 3*pi/4), (pi/4, pi/4), (3*pi/4, 3*pi/4), (5*pi/4, pi/4)]

 

SymPy 的矩阵工具

在矩阵论中,最常见的就是单位矩阵了,而单位矩阵只与一个参数有关,那就是矩阵的大小。下面就是 3*3,3*2,2*3 大小的矩阵。

>>> sympy.eye(3)
Matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]])
>>> sympy.eye(3,2)
Matrix([
[1, 0],
[0, 1],
[0, 0]])
>>> sympy.eye(2,3)
Matrix([
[1, 0, 0],
[0, 1, 0]])

另外还有全零和全一矩阵,意思就是矩阵中的所有值全部是零和一。

>>> sympy.ones(2,3)
Matrix([
[1, 1, 1],
[1, 1, 1]])
>>> sympy.zeros(3,2)
Matrix([
[0, 0],
[0, 0],
[0, 0]])

而对角矩阵也可以使用 diag 轻松获得:

>>> sympy.diag(1,1,2)
Matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 2]])

而矩阵的加法,减法,乘法,逆运算,转置,行列式,SymPy 都是可以支持的:

A = \left(\begin{array}{cc} 1 & 1 \\ 0 & 2 \end{array}\right)

B = \left(\begin{array}{cc} 1 & 0 \\ 1 & 1 \end{array}\right)

>>> A = sympy.Matrix([[1,1],[0,2]])
>>> B = sympy.Matrix([[1,0],[1,1]])
>>> A
Matrix([
[1, 1],
[0, 2]])
>>> B
Matrix([
[1, 0],
[1, 1]])
>>> A + B
Matrix([
[2, 1],
[1, 3]])
>>> A - B
Matrix([
[ 0, 1],
[-1, 1]])
>>> A * B
Matrix([
[2, 1],
[2, 2]])
>>> A * B**-1
Matrix([
[ 0, 1],
[-2, 2]])
>>> B**-1
Matrix([
[ 1, 0],
[-1, 1]])
>>> A.T
Matrix([
[1, 0],
[1, 2]])
>>> A.det()
2

在某些情况下,需要对矩阵进行加上一行或者加上一列的操作,在这里就可以使用 row_insert 或者 col_insert 函数:第一个参数表示插入的位置,第二个参数就是相应的行向量或者列向量。而在删除上就很简单了,直接使用 row_del 或者 col_del 即可。

>>> A
Matrix([
[1, 1],
[0, 2]])
>>> A.row_insert(1, sympy.Matrix([[10,10]]))
Matrix([
[ 1, 1],
[10, 10],
[ 0, 2]])
>>> A.col_insert(0, sympy.Matrix([3,3]))
Matrix([
[3, 1, 1],
[3, 0, 2]])
>>> A.row_del(0)
>>> A
Matrix([[0, 2]])
>>> A.col_del(1)
>>> A
Matrix([[0]])

在对角化方面,同样可以使用 eigenvals(),eigenvecs(), diagonalize() 函数:

>>> A
Matrix([
[1, 1],
[0, 2]])
>>> A.eigenvals()
{2: 1, 1: 1}
>>> A.eigenvects()
[(1, 1, [Matrix([
[1],
[0]])]), (2, 1, [Matrix([
[1],
[1]])])]
>>> A.diagonalize()
(Matrix([
[1, 1],
[0, 1]]), Matrix([
[1, 0],
[0, 2]]))

在 eigenvals() 返回的结果中,第一个表示特征值,第二个表示该特征值的重数。在特征向量 eigenvecs() 中,第一个表示特征值,第二个表示特征值的重数,第三个表示特征向量。在对角化 diagonalize() 中,第一个矩阵表示 P,第二个矩阵表示 DA = P*D*P^{-1}

在矩阵中,最常见的还是多元一次方程的解。如果要求 Ax =b 的解,可以有以下方案:

>>> A = sympy.Matrix([[1,1],[0,2]])
>>> A
Matrix([
[1, 1],
[0, 2]])
>>> b = sympy.Matrix([3,5])
>>> b
Matrix([
[3],
[5]])
>>> A**-1*b
Matrix([
[1/2],
[5/2]])
>>> sympy.linsolve((A,b))
{(1/2, 5/2)}
>>> sympy.linsolve((A,b),[x,y])
{(1/2, 5/2)}

 

SymPy 的集合论工具

集合论可以说是数学的基础,在任何数学的方向上都能够看到集合论的身影。在 SymPy 里面,有一个类叫做 sympy.sets.sets.set。在集合论里面,常见的就是边界,补集,包含,并集,交集等常见的操作。但是感觉 SymPy 中的集合论操作主要集中在实数域或者复数域上。

对于闭区间 I=[0,1] 和开区间 J = (0,1) 而言,在 SymPy 中使用以下方法来表示:

I = sympy.Interval(0,1)
J = sympy.Interval.open(0,1)
K = sympy.Interval(0.5,2)

其开始和结束的点可以分别使用 start 和 end 来表示:

>>> I.start
0
>>> I.end
1

其长度用 measure 来表示:

>>> I.measure
1

其闭包用 closure 来表示:

>>> I.closure
Interval(0, 1)

其内点用 interior 来表示:

>>> I.interior
Interval.open(0, 1)

判断其边界条件可以使用 left_open 或者 right_open 来做:

>>> I.left_open
False
>>> I.right_open
False

对于两个集合之间的操作,可以参考以下方法:

I = sympy.Interval(0,1)
K = sympy.Interval(0.5,2)
>>> I.intersect(K)
Interval(0.500000000000000, 1)
>>> I.union(K)
Interval(0, 2)
>>> I-K
Interval.Ropen(0, 0.500000000000000)
>>> K-I
Interval.Lopen(1, 2)
>>> I.symmetric_difference(K)
Union(Interval.Ropen(0, 0.500000000000000), Interval.Lopen(1, 2))

实数集 \mathbb{R} 在 SymPy 中用 sympy.S.Reals 来表示,自然数使用 sympy.S.Naturals,非负整数用 sympy.S.Naturals0,整数用 sympy.S.Integers 来表示。补集的计算可以用减号,也可以使用 complement 函数。

>>> sympy.S.Reals
Reals
>>> sympy.S.Reals-I
Union(Interval.open(-oo, 0), Interval.open(1, oo))
>>> I.complement(sympy.S.Reals)
Union(Interval.open(-oo, 0), Interval.open(1, oo))
>>> sympy.S.Reals.complement(I)
EmptySet()
>>> I.complement(K)
Interval.Lopen(1, 2)
>>> I.complement(sympy.S.Reals)
Union(Interval.open(-oo, 0), Interval.open(1, oo))

 

SymPy 的逻辑工具

在逻辑运算中,我们可以使用 A, B, C 来代表元素。&, |, ~, >> 分别表示 AND,OR,NOT,imply。而逻辑运算同样可以使用 sympy.simplify_logic 简化。

A, B, C = sympy.symbols("A B C")
>>> sympy.simplify_logic(A | (A & B))
A
>>> sympy.simplify_logic((A>>B) & (B>>A))
(A & B) | (~A & ~B)
>>> A>>B
Implies(A, B)

 

SymPy 的级数工具

SymPy 的级数工具有一部分放在具体数学(Concrete Mathematics)章节了。有的时候,我们希望计算某个级数是发散的,还是收敛的,就可以使用 is_convergence 函数。考虑最常见的级数:

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

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

>>> n = sympy.Symbol("n", integer=True)
>>> sympy.Sum(1/n, (n,1,sympy.oo)).is_convergent()
False
>>> sympy.Sum(1/n**2, (n,1,sympy.oo)).is_convergent()
True

如果想计算出收敛级数的值,加上 doit() 函数即可;如果想计算有效数字,加上 evalf() 函数即可。

>>> sympy.Sum(1/n**2, (n,1,sympy.oo)).evalf()
1.64493406684823
>>> sympy.Sum(1/n**2, (n,1,sympy.oo)).doit()
pi**2/6
>>> sympy.Sum(1/n**3, (n,1,sympy.oo)).evalf()
1.20205690315959
>>> sympy.Sum(1/n**3, (n,1,sympy.oo)).doit()
zeta(3)

除了加法之外,SymPy 也支持连乘,其符号是 sympy.Product,考虑

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

\prod_{n=1}^{+\infty}\cos\left(\frac{\pi}{n}\right)

>>> sympy.Product(n/(n+1), (n,1,sympy.oo)).is_convergent()
False
>>> sympy.Product(sympy.cos(sympy.pi/n), (n, 1, sympy.oo)).is_convergent()
True

 

SymPy 的 ODE 工具

在常微分方程(Ordinary Differential Equation)中,最常见的就是解方程,而解方程主要是靠 dsolve 函数。例如想求解以下的常微分方程:

df/dx + f(x) = 0,

d^{2}f/dx^{2} + f(x) = 0

d^{3}f/dx^{3} + f(x) = 0

可以使用 dsolve 函数:

>>> f = sympy.Function('f')
>>> sympy.dsolve(sympy.Derivative(f(x),x) + f(x), f(x))
Eq(f(x), C1*exp(-x))
>>> sympy.dsolve(sympy.Derivative(f(x),x,2) + f(x), f(x))
Eq(f(x), C1*sin(x) + C2*cos(x))
>>> sympy.dsolve(sympy.Derivative(f(x),x,3) + f(x), f(x))
Eq(f(x), C3*exp(-x) + (C1*sin(sqrt(3)*x/2) + C2*cos(sqrt(3)*x/2))*sqrt(exp(x)))

而常微分方程对于不同的方程类型也有着不同的解法,可以使用 classify_ode 来判断常微分方程的类型:

>>> sympy.classify_ode(sympy.Derivative(f(x),x) + f(x), f(x))
('separable', '1st_exact', '1st_linear', 'almost_linear', '1st_power_series', 'lie_group', 'nth_linear_constant_coeff_homogeneous', 'separable_Integral', '1st_exact_Integral', '1st_linear_Integral', 'almost_linear_Integral')
>>> sympy.classify_ode(sympy.Derivative(f(x),x,2) + f(x), f(x))
('nth_linear_constant_coeff_homogeneous', '2nd_power_series_ordinary')
>>> sympy.classify_ode(sympy.Derivative(f(x),x,3) + f(x), f(x))
('nth_linear_constant_coeff_homogeneous',)

 

SymPy 的 PDE 工具

在偏微分方程(Partitial Differential Equation)中,同样可以直接求解和判断偏微分方程的类型,分别使用函数 pdsolve() 和 classify_pde()。假设 f = f(x,y) 是一个二元函数,分别满足以下偏微分方程:

\partial f/\partial x + \partial f/\partial y =0

\partial f/\partial x + \partial f/\partial y + f = 0

\partial f/\partial x + \partial f/\partial y + f + 10 = 0

>>> f = sympy.Function("f")(x,y)
>>> sympy.pdsolve(sympy.Derivative(f,x)+sympy.Derivative(f,y),f)
Eq(f(x, y), F(x - y))
>>> sympy.pdsolve(f.diff(x)+f.diff(y)+f,f)
Eq(f(x, y), F(x - y)*exp(-x/2 - y/2))
>>> sympy.pdsolve(f.diff(x)+f.diff(y)+f+10,f)
Eq(f(x, y), F(x - y)*exp(-x/2 - y/2) - 10)

查看类型就用 classify_pde() 函数:

>>> sympy.classify_pde(f.diff(x)+f.diff(y)+f)
('1st_linear_constant_coeff_homogeneous',)
>>> sympy.classify_pde(f.diff(x)+f.diff(y)+f+10,f)
('1st_linear_constant_coeff', '1st_linear_constant_coeff_Integral')
>>> sympy.classify_pde(f.diff(x)+f.diff(y)+f+10,f)
('1st_linear_constant_coeff', '1st_linear_constant_coeff_Integral')

不过目前的 PDE 解法貌似只支持一阶偏导数,二阶或者以上的偏导数就不支持了。

 

SymPy 的数论工具

在数论中,素数就是一个最基本的概念之一。而素数的批量计算,比较快的方法就是筛法(sieve method)。在 sympy 中,同样有 sympy.sieve 这个工具,用于计算素数。如果想输出前100个素数,那么

>>> sympy.sieve._reset()
>>> sympy.sieve.extend_to_no(100)
>>> sympy.sieve._list
array('l', [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631])

如果想输出一个区间内的所有素数,可以使用 primerange(a,b) 函数:

>>> [i for i in sympy.sieve.primerange(10,100)]
[11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

search() 函数是为了计算某个数附近是第几个素数:

>>> sympy.sieve.search(10)
(4, 5)
>>> sympy.sieve.search(11)
(5, 5)

如果只想获得第 n 个素数,则使用函数 sympy.ntheory.generate.prime(n) 即可。如果是希望计算 x 后面的下一个素数,使用 sympy.ntheory.generate.nextprime(x) 即可。判断 x 是否是素数,可以使用 sympy.ntheory.generate.isprime(x)。

>>> sympy.ntheory.generate.prime(10)
29
>>> sympy.ntheory.generate.nextprime(10)
11
>>> sympy.ntheory.generate.nextprime(11)
13
>>> sympy.ntheory.generate.isprime(11)
True
>>> sympy.ntheory.generate.isprime(12)
False

除此之外,SymPy 的数论方法还有很多,需要读者根据 SymPy 的官方文档自行探索。

 

SymPy 的范畴论工具

SymPy 还支持范畴论(Category Theory)的一些计算方法,在这里简要地列举一下。

>>> A = sympy.categories.Object("A")
>>> B = sympy.categories.Object("B")
>>> f = sympy.categories.NamedMorphism(A,B,"f")
>>> f.domain
Object("A")
>>> f.codomain
Object("B")

由于范畴论是数学的“黑话”,因此其余方法留给范畴论的科研人员自行翻阅。

总结:

整体来看,SymPy 是一个非常卓越的 Python 开源符号计算库。在符号计算领域,不仅支持常见的微积分,线性代数,几何运算,还支持集合论,微分方程,数论等诸多数学方向。后续笔者将会持续跟进并研究这一卓越的开源工具库。

 

参考文献:

  1. Meurer A, Smith C P, Paprocki M, et al. SymPy: symbolic computing in Python[J]. PeerJ Computer Science, 2017, 3: e103.
  2. GitHub:https://github.com/sympy/sympy
  3. SymPy:https://www.sympy.org/en/index.html
  4. Sympy 维基百科:https://en.wikipedia.org/wiki/SymPy
  5. GreatX’s Blog:数值 Python:符号计算:https://vlight.me/2018/04/01/Numerical-Python-Symbolic-Computing/
  6. SymPy 符号计算-让Python帮我们推公式:https://zhuanlan.zhihu.com/p/83822118
  7. Python 科学计算利器—SymPy库:https://www.jianshu.com/p/339c91ae9f41

 

新加坡国立大学的数据科学与机器学习项目介绍

新加坡国立大学(National University of Singapore)是一所综合性的大学,根据泰晤士报和世界大学排名来看,NUS 在整个亚洲的排名是非常靠前的。同时,数学系(Department of Mathematics)则是在理学院(Science)下的一个院系。

NUS数学系的介绍

新加坡国立大学数学系的前身可以追溯到 1929 年的 Raffles College。当时理学院开设了数学,化学,物理三门课程,不过总共也就只有十个学生和三位教师,其中有一位是数学教师。第一届数学系的领导(从 1931 年到 1959 年)是Alexander Oppenheim教授,他是在美国芝加哥大学获得的博士学位。从 1929 年开始,在新加坡的教育系统中,数学教育事业得到了巨大的发展,对现在的新加坡国立大学和南洋理工大学的建立起到了至关重要的作用。

随着 NUS 的建立,数学系就进入了一个新的时代。新的校区在 Kent Ridge,1986 年理学院和数学系就在这里成立。这个时候,数学系就有了巨大的发展,不仅在本科生的招生规模方面有了巨大提高,在研究生项目规模上也有了一定的深度的提升。

NUS数学系首页

新加坡国立大学的数学系与国内的数学系有所不同。一般情况下,国内的数学系能够提供的专业包括数学与应用数学(Mathematics and Applied Mathematics),信息与计算科学(Information and Computing Science)与统计学(Statistics),有的时候会加上金融数学(Financial Mathematics)这一方向。而新加坡国立大学的数学系(Department of Mathematics)与统计系(Statistics)是分开的两个院系,虽然学生可以互相之间选择对方的课程,但是两者却是分属不同的院系。

NUS数学系本科专业

从数学系的首页来看,对于本科生而言,通常都有机会进行双学位的选择,例如:

  1. 计算机科学与数学;
  2. 经济学与数学;

因此,有不少的本科生都会有两个专业的学位证书。目标是为了让学生能够在未来从事数学研究和工业界的工作都打下坚实的基础。

NUS数学系硕士专业

NUS 的项目分成两块,coursework 项目和 Research 项目。第一个主要是以授课型的研究生为主,后者主要是为博士生或者准备攻读博士的人而准备的。对于硕士生而言,一般一两年就可以硕士毕业,只要修课的学分满了即可。

对于众多硕士生而言,进入 NUS 之前就要根据自身情况来选择一个合适的方向进行申请。对于硕士生的项目,数学系可以提供数学,金融数学等方向的课程,并且近期也提供了数据科学与机器学习专业(Data Science and Machine Learning)的课程。

项目介绍

申请条件

数据科学与机器学习项目是数学系,统计系,计算机系联合举办的为期一至两年的硕士生项目。期望学生的本科背景是数学方向,应用数学方向,统计与物理方向等。同时对学生的英语能力有一定的要求,希望对于母语不是英语的学生能够达到托福 85 分以上或者雅思 6.0 分以上的成绩。

项目课程要求

在课程安排方面,这个项目会有 20 个学分的课程,5 门核心课程,包括:

  1. Introduction to Big Data for Industry;
  2. Optimisation for Large-Scale Data-Driven Inference;
  3. Foundations of Machine Learning/Theory and Algorithms for Machine Learning;
  4. Cloud Computing;
  5. DSML Industry Consulting and Applications Project.

而选修课方面包括机器学习,数据挖掘,大数据,计算机视觉,金融数学等方向的专业课程。同样也会要求选择 20 个学分的选修课程即可。

学费

在学费方面,该项目是没有奖学金的,而一年的学费是 45000 新币,大约为 23 万人民币。

NUS申请流程

如果对这个项目感兴趣的同学,需要在 2020 年的 3 月 15 日之前提交申请,才有机会在 2020 年 8 月入学。

NUS数学系联系方式

如果有任何问题的话,可以考虑发邮件给 askmathpg@nus.edu.sg 或者拨打上面的电话号码。

 

读博的意义

读博的意义

在本科数学系,身边的学生大致分为几种情况:第一种学生是高中数学学的比较好,也没有什么特别爱好的专业,于是通过高考或者保送进入了数学系;第二种学生是为了转行计算机或者金融,就来数学系过渡一下,未来从事其他行业的工作;第三种学生是想读其他的专业,但是由于高考分数不足,被分配到了数学系;第四种学生就是为了攀登科学的高峰,立志要在数学界从事科研工作,在填报志愿的时候就果断进入数学系。一般情况下,学生在拿到本科学位之后,如果还想继续从事数学科研的话,就会选择继续攻读数学方向的博士学位,无论在国内还是国外继续学业。640-6

在数学系本科或者硕士的时候,学生们通常都是按部就班的上课,做作业,上讨论班,或者每周集体分享一些论文或者学术资料。到了最后一个学期需要写论文的时候其实也不需要有太多的创新工作(更直接一点的话就是做不出来什么创新的工作),学生们大概就是整理一篇前沿的学术论文或者学术书籍就算完成了论文。在这两个阶段,导师基本上都会告诉学生,什么课题可以做得出来,什么课题比较难被做出来。或者就是导师其实已经知道某个问题的解决步骤和方法了,只需要学生执行完成就可以了。虽然在数学系的本科或者硕士都是按部就班的学习,但是这种学习是非常有必要的,毕竟不积跬步无以至千里,不积小流无以成江海。

QQ20191027-154338@2x

但是到了博士阶段就不一样了,导师给博士生的题目基本上就是学术界未知的问题。如果这个问题被解决了,那就不能够作为一个博士论文的题目。在某些极端情况下,导师连这个题目最终能不能够被学生做出来都不知道,只能够大致判断出这个问题的难易程度。数学系的纯数学方向和其他专业还不一样,某些专业跑跑数据改改模型就能够发表所谓的优秀论文了。而想要做纯数学的研究不仅需要通过天马行空的想象能力来寻找灵感,还需要有强大的逻辑推理能力把之前的想象整理成严谨的学术论文。而这种能力,仅仅靠本科或者硕士的训练是远远不够的,按部就班的读书考试也许能够应付小学,中学,大学,甚至研究生阶段,但是却无法完成博士阶段的任务目标。

QQ20191027-154259@2x

当年笔者在博士期间的研究方向是“复平面上多项式的 Julia 集合正测度”(PS:估计不是相关方向的人都不太读得懂这个题目的意思)。首先要知道复平面,其次要知道 Julia 集合,再次要知道正测度是什么意思。即使这些名词和概念都知道了,也不足以撬动这个题目。要想解决这个题目,除了导师的必要指导之外,还要自行去查阅各种论文资料来详细阅读。记得当年读过好几篇三五十页的 Annals of Mathematics(数学年刊)的论文,还读过一篇上百页的预印稿(关键是这篇论文有一个核心步骤还是错误的)。在此情况下,导师也只能够给博士生圈定一个论文的范围,告诉学生可以去参考其中的思路和方法,至于学生能不能够读懂这些文献,是导师无法保证的。而纯数学论文与其他专业最明显的不同就是满篇都是数学公式的推导,有很多地方充斥着“显然”,“易得”,“显而易见”诸多词汇。有一种情况是写论文的作者没有想明白,然后糊弄了一把;另外一种情况是这个地方真的是显而易见的,只是读者没有明白。纯数学论文的其中一页读上一两天并不是一件罕见的事情,一篇论文读一个学期能读明白也算完成了一件还不错的任务。PS:数学系有的教材就可以让学生一页读上好几天,例如 GTM 52。某些专业做不出来实验还可以把原因归结为材料不足,经费不足,设备不够,但是纯数学专业看不懂论文只能够把绝大部分原因归结为自己,因为草稿纸,笔,打印机,网络都是买得起的,客观上并不存在任何阻碍。

Julia_0.4_0.6Julia-set_z2+c_-0.742_0.1_5000

在整个科研的过程中,一般都会经历:高峰->低谷->逐步爬坡回升->重新达到高峰的这样一个流程状态。正如:邓宁-克鲁格心理效应。从博士生的五年生涯来看,通常都会经历这几个不同的阶段:

1. 愚昧山峰:刚入学第一年的时候;

2. 低谷:第二年,第三年的时候;

3. 回升:第四年,第五年的时候。

而这种状态的经历是在读本科,硕士阶段都无法获得的,只有在经历了博士的求学之路才能够逐渐体会。当然,也有的导师也会帮学生写论文,整理论文,甚至送论文给学生,这样的学生可以在学生圈谋求一份教职,但是也永远体会不到这些心路历程了。这份心路历程在人生中是一场宝贵的财富,教会了如何在绝望之谷的时候迅速反弹回升,如何在逆境中持续成长,如何在被打击之后迅速恢复。

640-7

后来由于种种原因笔者没有继续从事博士期间的科研项目,到了企业从事普通的研发工作。在工作中,通常领导们也不会直接告诉下属项目的细节和执行方案,只会大致说出项目方向,该达成的目标是什么,下属们需要自行思考这个项目的背景是什么,为什么要启动这个项目,该怎么实现这个项目,怎么样做才能够把项目做好做完。这个过程跟读博期间的过程有点类似,但是又有所不同。相同的是都只能够拿到一个大致的方向,同时需要发挥自主能动性去做完细节,争取身边的各种资源,并且保证论文或者项目的落地。不同的是,工作期间并不会给下属很长的空窗期,并不会像博士期间给学生大几个月甚至一年的空闲时间去研究某个课题或者方向。基本上都是一旦制定了方向,就要立刻给出大致的方案和解决思路,然后给出项目的排期和执行计划,并且给出项目的风险点是什么,最终开始逐步落实方案和具体的工作内容。不过,在工作中,在大多数情况下也不会要求一定要做一个行业第一的产品,而且很多项目由于种种原因也无法做到行业第一。有的时候,行业第一的量化标准还不太好做,如果有明确规则的,例如挣钱的多少,日活月活多少,比赛的第一,这种是可以制定出来的。但是大多数情况下其实也不太好制定出来。不过在工作中还是可以精益求精,以行业的第一为目标来逐步提升的。可以在现有的情况和条件下把项目做得尽量完善,在提升效率方面比之前能够更加好,在节省人力方面能够做到较少人力就维护较多的事情。

在工作之后,虽然也很难体会到一篇论文读半年的时候,但是在博士期间学会的资料整理,收集,汇总,甚至日常的笔记和资料总结,都是在工作中非常常用的技能。甚至一些时间管理和项目管理的经验,也是在工作中必备的技能之一。

640-2

640-3
PHD期间的时间安排

640-4

640-5
工作的时间安排

很难说读博这件事情对所有的人都有意义,因为每个人都有着不同的生长环境和家庭条件,在如此复杂的情况下,给出一个固定的答案其实没有任何意义。何况每个人的想法也会随着时间的迁移,环境的变化而做出改变。之前想做的一些事情可能后面就不感兴趣了,一直想追求的一些目标可能也不再继续追求了。只能够根据具体的情况,在不同的阶段来定制化地给每个人相应的建议和支持(个性化推荐)。无论做什么样的事情,最好都要保持一种乐观向上的心态,不能够一直处于一种低估和不能自拔的阶段。如果在本科,硕士,甚至博士期间都没有找好自己的定位,也可以尝试在工作之后逐步把自己拉上正规,在期间调整自己的状态和节奏,让自己的幸福感逐步提升。无论在哪里,无论最终是否做科研,最终都是为了过上更好的生活。

张戎

2019年10月27日

深度学习在时间序列分类中的应用

本篇博客将会分享几篇文章,其内容主要集中在深度学习算法在时间序列分类中的应用。

无论是图像分类,文本分类,还是推荐系统的物品分类,都是机器学习中的常见问题和应用场景。同样的,时间序列的分类问题也是研究时间序列领域的重要问题之一。近期,神经网络算法被用于物体识别,人脸识别,语音分类等方向中,于是有学者用深度学习来做时间序列的分类

假设

X=\{x_{1},\cdots,x_{n}\}

是一个长度为 n 的时间序列,高维时间序列

X = \{X^{1},\cdots,X^{M}\}

则是由 M 个不同的单维时间序列而组成的,对于每一个 1\leq i\leq M 而言,时间序列 X^{i} 的长度都是 n. 而时间序列的分类数据通常来说都是这种格式:数据集

D=\{(X_{1},Y_{1}),(X_{2},Y_{2}),\cdots,(X_{N},Y_{N})\}

表示时间序列与之相应的标签,而 Y_{i} 是 one hot 编码,长度为 K(表示有 K 个类别)。

整体来看,时间序列分类的深度学习方案大体是这个样子的:输入的是时间序列,通过某个神经网络算法进行端到端的训练,最后输出相应的分类概率。

ReviewFigure1.png

而做时间序列分类的深度学习算法分成生成式(Generative)判别式(Discriminative)两种方法。在生成式里面包括 Auto Encoder 和 Echo State Networks 等算法,在判别式里面,包括时间序列的特征工程和各种有监督算法,还有端到端的深度学习方法。在端到端的深度学习方法里面,包括前馈神经网络,卷积神经网络,或者其余混合模型等常见算法。

ReviewFigure5.png

深度学习算法在时间序列分类中的应用:Baseline

这一部分将会介绍用神经网络算法来做时间序列分类的 Baseline,其中包括三种算法,分别是多层感知机(MLP),FCN(Fully Convolutional Network)和 ResNet。其论文的全名是《Time Series Classification from Scratch with Deep Neural Networks: A Strong Baseline》。这篇论文中使用的神经网络框架如下图所示:

DNN_Baseline_结构1.png

多层感知机(MLP)模型使用了全连接层,每个隐藏层大约 500 个神经元,然后使用 ReLU 作为激活函数,同时使用 Dropout 来防止过拟合,最后一层是 Softmax 层。MLP 中一个基础的块包括:

\tilde{x} = f_{dropout, p}(x),

y = W\cdot \tilde{x} + b,

h = ReLU(y).

除了前馈神经网络之外,全卷积网络(FCN)同样可以作为时间序列的特征提取工具,一个卷积块包括:

y = W \otimes x + b,

s = BN(y),

h = ReLU(s),

在这里,\otimes 指的是卷积算子,BN 指的是 Batch Normalization,ReLU 则是激活函数。

Residual Network 是在 FCN 的基础上进行的改造。令 Block_{k} 来表示第 k 个卷积块,而 Residual 块就定义为:

h_{1} = Block_{k_{1}}(x),

h_{2} = Block_{k_{2}}(h_{1}),

h_{3} = Block_{k_{3}}(h_{2}),

y=h_{3} + x,

\hat{h}=ReLU(y).

其中,k_{1} = 64, k_{2} = 128, k_{3} = 128

ReviewFigure6_ResNet.png

评价指标

Mean Per Class Error (in Multi-class Classification only) is the average of the errors of each class in your multi-class data set. This metric speaks toward misclassification of the data across the classes. The lower this metric, the better.

模型的评价指标使用的是 Mean Per-Class Error,指的是在多分类场景下,每一类(Class)错误率的平均值。换句话说,一个数据集 D=\{d_{k}\}_{1\leq k\leq K} 是由 K 个类的元素构成的,每个类的标签是 C=\{c_{k}\}_{1\leq k\leq K},通过模型其实可以计算出模型对每一个类的错误率 e_{k},那么模型的 MPCE 就是:MPCE= \sum_{1\leq k\leq K} e_{k}/K.

其实验结论是:

DNN_Baseline_实验数据1.png

MSCNN

MSCNN 的全称是 Multi-Scale Convolutional Neural Networks,相应的论文是《Multi-Scale Convolutional Neural Networks for Time Series Classification》

在时间序列的分类算法里面,通常来说,可以分成以下几种:

  1. 基于距离的方法(distance-based methods):kNN,SVM(相似核),DTW;
  2. 基于特征的方法(feature-based methods):SVM,逻辑回归等;
  3. 基于神经网络的方法(neural network-based methods):CNN 等;

正如前文所提到的,一条时间序列通常可以写作 T=\{t_{1},\cdots,t_{n}\},其中 t_{i} 表示在时间戳 i 下的取值,并且时间序列 T 的长度是 n。在时间序列分类的场景下,每一条时间序列对应着唯一的一个标签(label),也就是说 D=\{(T_{i},y_{i})\}_{i=1}^{N}。其中 D 集合里面包含 N 条时间序列,每条时间序列 T_{i} 对应着一个标签 y_{i}y_{i} 表示分类值集合 \mathcal{C} = \{1,\cdots,C\} 中的元素,C\in \mathbb{Z}^{+}

MSCNN 的整体结构:

在 Multi-Scale Convolutional Neural Network(MSCNN)中,包括几个串行的阶段,

  1. 变换阶段(Transformation Stage):包括恒等变换,下采样,谱变换等变换方式,每一种方式都是一个分支,并且也是卷积神经网络的输入;
  2. 局部卷积(Local Convolution Stage):使用卷积层来对不同的输入提取特征,不同的输入分支之间是相互独立的,输出的时候都会经过一个最大值池化(max pooling)的过程;
  3. 整体卷积(Full Convolution Stage):把上一步提取到的特征进行拼接(concatenate),然后使用全连接层并且加上一个 softmax 层来做多分类。

如下图所示,MSCNN 是一个端到端的训练网络结构,所有参数都是通过后向传播算法得到的。

MSCNN结构1.png

首先来看神经网络的第一步,变换阶段(Transformation Stage),也就是神经网络的多尺度的输入。在不同的尺度下,神经网络能够提炼到不同类型的特征。长期的特征(long-term features)反映了时间序列的整体趋势,短期的特征(short-term features)反映了时间序列的局部的微妙变化。要想判断时间序列的形状,不仅要参考整体的特征,也要参考局部的特征,这两者对于判断时间序列的形状都具有一定的辅助作用。

在 Transformation Stage,identity map 指的是恒等变换,也就是说时间序列是原封不动的作为神经网络的输入数据。对于 Smoothing Transformation,指的就是对时间序列进行必要的平滑操作,将新的时间序列作为神经网络的输入数据。在这种情况下,我们可以对时间序列 T=\{t_{1},\cdots,t_{n}\} 进行移动平滑,i.e.

T^{\ell}=(x_{i}+x_{i+1}+\cdots+x_{i+\ell-1})/\ell, 0\leq i\leq n-\ell+1,

其中的 \ell\in \mathbb{Z}^{+} 表示窗口长度。对于不同的窗口长度 \ell,我们可以的到不同的时间序列平滑序列,但是它们的长度都是一样的,都是原始的时间序列长度 n

而下采样(down sampling)指的则是对时间序列的间隔进行抽样操作。假设时间序列 T=\{t_{1},\cdots,t_{n}\},下采样的比例是 k,也就是说我们每隔 k 个点保留时间序列的取值,i.e.

T^{k} = \{t_{1+k\cdot i}\}, 0\leq i \leq [(n-1)/k].

用这种方法,我们可以对 k=1,2,3,\cdots 来进行下采样的时间序列提取。在进行了恒等变换,平滑变换,下采样之后,时间序列就可以变成多种形式,作为神经网络的输入。

其次,在神经网络部分,本文使用了一维(1-D)的卷积层和最大值池化的方法来提取特征,并且在局部卷积阶段之后把提炼到的抽象特征进行拼接(concatenate)。拼接完了之后,持续使用卷积层和池化层进行特征的提取,然后使用全连接层(fully connected layers)和 softmax 层来进行时间序列类别的预测。

数据增强

在深度学习里面,由于是端到端的训练网络,因此是需要相对多的样本数据的,于是有的时候需要进行数据增强(data augmentation)。也就是在现有的基础上获得更多的训练数据。对于时间序列 T=\{t_{1},\cdots,t_{n}\},可以定义一个子序列:

S_{i:j} = \{t_{i}, t_{i+1},\cdots,t_{j}\}, 1\leq i,\leq j\leq n

对于正整数 s\in \mathbb{Z}^{+},可以生成 n-s+1 个子序列如下所示:

Slicing(T,s) = \{S_{1:s},S_{2:s+1},\cdots,S_{n-s+1:n}\},

这些子序列的标签与原始的时间序列 T 是一样的。

本文用到的数据集情况如下表所示:

MSCNN数据集1.png

实验数据如下图所示:

MSCNN实验1MSCNN数据集2

结论

本文使用了 MCNN 来对变换之后的时间序列进行特征提取,并且进行了端到端的模型训练。并且也讨论了卷积神经网络使用在 shapelet learning 上的一些逻辑和方法,然后解释了 MSCNN 在时间序列分类上能够有不错表现的原因。但是所有的 TSC 数据集都不算特别大,对端到端的训练模式有一定的限制。

GASF 和 GADF 方法

这篇文章《Imaging Time Series to Improve Classification and Imputation》介绍了如何把时间序列转换成图像,包括 GASF 方法和 GADF 方法。

假设时间序列是 X = \{x_{1},\cdots, x_{n}\},长度是 n,我们可以使用归一化方法把时间序列压缩到 [0,1] 或者 [-1,1]

\tilde{x}_{0}^{i} = (x_{i} - \min(X))/(\max(X) -\min(X)),

\tilde{x}_{-1}^{i} = ((x_{i} - \max(X)) + (x_{i}-\min(X)))/(\max(X) - \min(X)),

此时的 \tilde{x}_{0}^{i}\in[0,1], \forall 1\leq i\leq n\tilde{x}_{-1}^{i} \in [-1,1],\forall 1\leq i\leq n。于是可以使用三角函数来代替归一化之后的值。下面通用 \tilde{x}_{i} 来表示归一化之后的时间序列,令 \phi_{i} = \arccos(\tilde{x}_{i})\tilde{x}_{i} \in [-1,1]1\leq i\leq n。因此,\phi_{i}\in[0,\pi],于是,\sin(\phi_{i}) \geq 0

定义矩阵 GASF(Gramian Angular Summation Field)

GASF = [\cos(\phi_{i}+\phi_{j})]_{1\leq i,j\leq n}

于是,

GASF = [\cos(\phi_{i})\cdot \cos(\phi_{j}) - \sin(\phi_{i})\cdot \sin(\phi_{j})]_{n\times n}

\tilde{X}=(\cos(\phi_{1}),\cdots,\cos(\phi_{n}))^{T},可以得到

GASF = \tilde{X} \cdot \tilde{X}^{T} - \sqrt{I - \tilde{X}^{2}} \cdot \sqrt{I - \tilde{X^{T}}^{2}}

以上的都是 element 乘法和加法,I 表示单位矩阵。它的对角矩阵是

diag(GASF) = \{\cos(2\phi_{1}),\cdots, \cos(2\phi_{n})\}

= \{2\cos^{2}(\phi_{1})-1,\cdots,2\cos^{2}(\phi_{n})-1)\} = \{GASF_{ii}\}_{1\leq i\leq n}.

如果是使用 min-max normalization 的话,是可以从 diag(GASF) 反推出 \tilde{x}_{i} 的。因为,2 \tilde{x}_{i}^{2} - 1 = 2\cos^{2}(\phi_{i}) - 1 = GASF_{ii},可以得到 y_{i} = \sqrt{(GASF_{ii}+1)/2}

定义 GADF(Gramian Angular Difference Field)如下:

GADF = [\sin(\phi_{i}-\phi_{j})]_{1\leq i,j\leq n}

= [\sin(\phi_{i})\cdot cos(\phi_{j}) - \cos(\phi_{i})\cdot\sin(\phi_{j})]_{1\leq i,j\leq n}

= \sqrt{1-X^{2}}\cdot X^{T} - X \cdot \sqrt{1-(X^{T})^{2}}.

GASFGADF1.png

Markov Transition Field(MTF)

除了 GSAF 和 GSDF 之外,《Imaging Time Series to Improve Classification and Imputation》,《Encoding Time Series as Images for Visual Inspection and Classification Using Tiled Convolutional Neural Networks》,《Encoding Temporal Markov Dynamics in Graph for Time Series Visualization》也提到了把时间序列转换成矩阵 Image 的算法 MTF。在 pyts 开源工具库里面,也提到了 MTF 算法的源码。

假设时间序列是 X = \{x_{1},\cdots,x_{n}\},我们把它们的值域分成 Q 个桶,那么每一个 x_{i} 都可以被映射到一个相应的 q_{j} 上。于是我们可以建立一个 Q\times Q 的矩阵 Ww_{ij} 表示在桶 j 中的元素被在桶 i 中的元素跟随的概率,也就是说 w_{ij} = P(x_{t}\in q_{i}|x_{t-1}\in q_{j}),同时,它也满足 \sum_{j=1}^{Q}w_{ij} =1。于是,得到矩阵 W = (w_{ij})_{1\leq i,j\leq Q}

MTF1.png

除此之外,我们也能够计算一个迁移概率矩阵 M。其中 $m_{ij}$ 表示桶 i 中的元素迁移至桶 j 中的概率 P(q_{i}\rightarrow q_{j}),同样有 \sum_{1\leq j\leq Q} m_{ij} =1。因此,我们同样可以构造出一个 Q\times Q 的矩阵将时间序列可视化。

时间序列的降维方法有两种:

  1. 分段聚合(PAA):使用局部平均等方法,把时间序列进行降维;
  2. 核变换(Kernel):使用 Bivariate Gaussian 核或者均值核来把时间序列进行降维。

在把时间序列进行可视化之后,对于时间序列分类的场景,就可以使用 CNN 的技术方案来做了。如下图所示:

tildeCNN.png

其实验数据效果如下:

tildeCNN实验数据1.png

Time Le-Net

在本篇文章《Data Augmentation for Time Series Classification using Convolutional Neural Networks》中,主要用到了卷积神经网络来做时间序列的分类。

time_lenet_1.png

除此之外,也使用了不少数据增强(Data Augmentation)的技术。包括前面提到的 Window Slicing(WS)方法。也考虑了 Warping 的变换技巧,例如 Warping Ratio = 1/2 或者 2。这种时间扭曲指标比率可以通过交叉验证来选择。该方法叫做 Window Warping(WW)技术。

另外也有其余论文使用卷积神经网络做时间序列分类,例如《Convolutional neural networks for time series classification》,如下图所示:

CNN_figure1.png

Multi-Channels Deep Convolutional Neural Networks

在高维时间序列的分类中,有人提出用多通道的卷积神经网络来进行建模。

multichannelcnn_figure3.png

整体来看,分成四个部分。前三个部分作为特征提取工具,最后一层作为分类工具。

  1. Filter Layer:
  2. Activation Layer:
  3. Pooling Layer:
  4. Fully-Connected Layer:

实验对比数据如下:

multichannelcnn_table1.png

结论:

在本篇博客中,列举了一些深度学习算法在时间序列分类中的应用,也介绍了部分数据增强的方法和时间序列数据变换的方法。从以上各篇文章的介绍来看,深度学习在时间序列分类领域上应该是大有可为的。

参考资料:

  1. Wang Z, Yan W, Oates T. Time series classification from scratch with deep neural networks: A strong baseline[C]//2017 international joint conference on neural networks (IJCNN). IEEE, 2017: 1578-1585.
  2. Cui Z, Chen W, Chen Y. Multi-scale convolutional neural networks for time series classification[J]. arXiv preprint arXiv:1603.06995, 2016.
  3. Wang Z, Oates T. Imaging time-series to improve classification and imputation[C]//Twenty-Fourth International Joint Conference on Artificial Intelligence. 2015.
  4. Wang Z, Oates T. Encoding time series as images for visual inspection and classification using tiled convolutional neural networks[C]//Workshops at the Twenty-Ninth AAAI Conference on Artificial Intelligence. 2015.
  5. Liu L. Encoding Temporal Markov Dynamics in Graph for Time Series Visualization, Arxiv, 2016.
  6. Fawaz H I, Forestier G, Weber J, et al. Deep learning for time series classification: a review[J]. Data Mining and Knowledge Discovery, 2019, 33(4): 917-963.
  7. Zhao B, Lu H, Chen S, et al. Convolutional neural networks for time series classification[J]. Journal of Systems Engineering and Electronics, 2017, 28(1): 162-169.
  8. Le Guennec A, Malinowski S, Tavenard R. Data augmentation for time series classification using convolutional neural networks[C]. 2016.

时间序列的标签

本篇文章是为了介绍一种基于少量样本标记而获得更多样本的方法,论文的原文是《Label-Less: A Semi-Automatic Labeling Tool for KPI Anomalies》,是清华大学与多家公司(必示科技,中国建设银行等)的合作论文。

在时间序列异常检测中,因为标注的成本比较大,于是需要寻找一种较少而高效地标注时间序列异常点的方法。在该论文中,Alibaba,Tencent,Baidu,eBay,Sogou提供了上千条时间序列(每条时间序列大约是2-6个月的时间跨度),作者们进行了 30 条 KPIs 的标注工作。但是其标注成本依旧是很大的,于是作者们想到了一种异常相似搜索(anomaly similarity search)的算法,目标是对已经标注好的时间序列异常模式进行模版搜索。目的就是达到 label-less,也就是较少的标注而获得更多的标注数据。

在本篇论文中,在异常检测的过程中,作者们使用了时间序列的预测模型(time series prediction models)来获得时间序列的特征,使用了孤立森林(Isolation Forest) 来对时间序列的特征来做无监督的异常检测。并且其效果由于 one class svm 算法和 local outlier factor 算法。在搜索的部分,作者使用了加速版的 DTW 算法(accelerated dynamic time warping approach)来做相似度的搜索和模式的匹配。其中也尝试了各种技巧和方法,包括 constrained DTW,LB Keogh 方法,early stopping 算法等工具。

整个 Label-Less 的架构图如下表示:

Fig2.png

其中的 Operators 指的是业务运维人员,面对着无标记的多条时间序列曲线。系统首先会进行无监督的异常检测算法啊,包括时间序列的预处理(归一化等)操作,然后使用差分(Difference),移动平均算法(moving average),带权重的移动平均算法(weighted moving average),指数移动平均(ewma),holt winters,ARIMA 等算法来做特征的提取。此时,对于不同的时间序列预测工具,我们可以得到不同的预测值,然后把预测值减去实际值并且取绝对值,就得到时间序列的误差序列。i.e. |p_{i} - x_{i}| 就作为数据点 x_{i} 的特征。

在这种情况下,由于用了六个时间序列预测算法,因此原始的时间序列 X (n\times 1) 就可以变成特征矩阵 X' (n\times 1)。对于特征矩阵 X' 可以使用 isolation forest 来做无监督的异常检测并且做阈值的设定;如下图所示:

Fig3.png

而另外的一部分的异常相似搜索(anomaly similarity search)是在第一部分的基础上在做的,Unsupervised Anomaly Detection 会输出疑似异常或者候选异常,并且基于已知的异常模板(Anomaly Template)进行相似度的匹配,此时可以使用 accelerated DTW 算法,选择出最相似的 Top-K 异常,然后运维人员进行标注,得到更多的样本。

由于,对于两条长度分别是 mn 的时间序列,DTW 相似度算法的时间复杂度是 O(mn),因此在搜索的时候需要必要的加速工作。在这种地方,作者们使用了 LB-Kim,LB-Keogh,LB-Keogh-Reverse 算法来做搜索的加速工作。而这些的时间复杂度是 O(m+n)。整体的思路是,如果两条时间序列 qc 的 LB-Kim,LB-Keogh,LB-Keogh-Reverse 的下界大于某个阈值,则不计算它们之间的 DTW 距离。否则就开始计算 DTW。并且在计算 DTW 的时候,如果大于下界,则会提前终止(early stopping),不会继续计算下去。如果都没有大于阈值,则把这个候选曲线和 dist 距离放入列表,最后根据列表中的 dist 来做距离的逆序排列。

整体流程如下:

AnomalySimilaritySearch

其运行速度也比直接使用 DTW 快不少:

Table3.png

Label-Less 的交互页面如下所示:

图(a)表示使用无监督算法获得的疑似异常;

图(b)表示使用异常搜索算法获得的异常结果。

Fig10.png

下图则表示模板,m 表示模板的长度,c 表示相似的异常候选集个数;

Fig11

总结:

整体来看,本文提供了一种通过少量人工标注无监督算法相似度算法来获得更多样本的方法。在候选的时间序列条数足够多的时候,是可以进行时间序列的相似度匹配的。这给未来在运维领域提供海量的时间序列标注数据给予了一定的技术支持。

 

怀念 Jett

​之前笔者写过不少文章,其中不仅包括技术和数学文章,还有一些日常学习工作的心得,甚至还写过几篇关于旅游的流水账,但是却很少写关于人物的文章。最近我打算写一篇关于人物的文章,用于纪念在工作中一起战斗过的伙伴 Jett。

许多年前,我还在学校读书的时候,就有前辈告诉我,交朋友的时候基本上就只能够在学校,离开了学校之后在职场中想交到合适的朋友是比较困难的。当时作为一枚职场菜鸟的我,也只能够表示对这句话的认同。直到在工作后遇到了 Jett,我对这个观点就有了自己的看法。

屏幕快照 2019-08-10 下午10.02.49

我是2015年7月作为校招生进入公司,Jett 则是2015年10月通过社招进公司的。当年虽然两人身处同一个部门,但是由于不在一个中心,也不是同一条业务线,并且这两个中心平时也没啥日常工作的交集,因此我和 Jett 同处一个部门许久却互相不认识对方。就好比在学校里面,虽然身处同一个院系,甚至同一栋宿舍楼,由于种种原因而不认识对方。

随着工作时间的增加,到 2016 年的时候笔者就已经过了新手保护期。所谓新手保护期,就是每一个新人其实都有一段适应公司的过程,公司也会给人一定的保护时间,但是一旦过了这个时间就不算新人了。大约到了2016年5月份的时候,部门内部启动了一个用机器学习的方法来做安全对抗的项目。当年的安全中心并没有机器学习的从业者,因此需要从其他中心抽调人手去支持他们的项目。在一些机缘巧合之下,我作为数据团队的代表,就和安全团队的同事们同在安全大数据项目组。在公司里面,机器学习的人员通常都是和开发人员通力合作,各司其职,共建一个项目。在项目初期,其实 Jett 也不负责与我对接,后续随着项目的进一步发展和迭代,Jett 就成为了安全中心与我对接的人。

虽然 Jett 是和我对接的人,但是 Jett 主要也是负责平台开发,之前 Jett 在进入公司之前好像没有做过前端,但是在项目初期好像做了一些前端方面的工作,在邮件里面得到了领导的赞扬。这个应该是我对 Jett 的第一印象了,那就是一个靠谱的开发人员。当时虽然是成立了项目组,但是还是作为一个试点的工作来做,无论是算法的调研和平台的工作,都由我们两个来完成。这些项目在初期都是一大堆人在里面,但是随着项目的进展,人数都会越来越少,然后领奖的时候又出现了一大堆人。在项目的初期,我当时应该也只会 SQL,连 Python 都不算很熟悉,更没有开发过大型项目了。即使是作为一个试点的项目,其实个人也有着巨大的压力,毕竟这个方向在公司当年还处于一个相对不那么清晰的状态,也不知道该怎么做,甚至也不知道有什么资料可以查。不过作为一个 PHD,虽然写代码的能力不行,但是读论文和搜索资料的能力还是有点的。当时记得花了很多时间在网上找来找去,总算找到一两家号称做这个方向的创业公司,运气够好的是当时它们居然都发表了相关的论文。后来把它们的论文整理成资料和PPT,跟相关的同事沟通了一下,就打算复现这几篇论文。

屏幕快照 2019-08-10 下午10.01.49

当然,做项目总是有着各种各样的风险,况且当年本人的技术实力确实也不怎么样,写过不少的 SQL,做过很多数据分析的工作,查过一些数据方面的问题,至于其他方面也就那么回事。在项目的调研阶段和初期,大约是5月份-8月份,其实也没有什么特别大的产出,主要也是在积攒技术经验,把 Python 的各个工具库熟悉了一下。不过在公司里面总是有考核的,到了9月份的时候压力就比较大了,毕竟最终还是要对上线的效果负责的。在9月底的时候,我和 Jett 被拉去星巴克喝了一杯咖啡,然后被下达了任务和命令,也就是在十月份完成机器学习的效果指标。

屏幕快照 2019-08-10 下午9.59.46

既然被拉到星巴克去喝咖啡,也就是说明在9月份的时候,模型的效果都不算太好,估计再这样干下去也没有办法完成既定的目标。不过目标这种事情制定的时候都是根据当时的情况来制定的,很有可能会随着时间的变化而产生变化。虽然说目标总是在发生变化,但是整体的大思路是没有改变的。当时,作为开发人员的 Jett 问我要不要修改 XGBoost 的模型源码,我连忙说不用不用,其实就算让我改,我也没有本事修改别人的源码。既然不会修改模型的源码,那就只能够从数据和特征来入手了。在机器学习领域,只要把数据和特征处理好了,基本上就能够保证模型的效果和质量。于是 Jett 和我当时把成千上万条数据一起看了一遍,结果是发现正负样本有一部分混在了一起,于是模型的效果无论如何都做不好。

屏幕快照 2019-08-10 下午10.05.45
XGBoost 的 Github 主页

对于在校学生而言,通常来说都是 Python 跑完模型,然后得到一个模型文件,用它继续离线预测就可以得到最终的结果。但是在工业界,很多机器学习项目都需要进行上线的工作。这种时候只靠一个机器学习人员的战斗力是无法解决问题的,不可避免地需要有开发人力的介入。此时 Jett 发挥了作为一位开发人员的强大战斗力,一个人就能把线上的代码全部完成,没有让我撰写任何一行 C++ 的代码。在码农界有一种“结伴编程”的说法,也就是两个人共同搞一份代码,共同搞一个项目。我是负责数据处理和离线模块,Jett 是负责平台开发。有一次在核对数据的时候,把 Python 的预测结果和 C++ 的预测结果进行核对。也就是为了交流方便,那次我把 Mac Air 搬到了 Jett 的桌子上,两个人并排坐在一起,一起核对数据的准确性和可靠性。

不过核对数据只是机器学习项目的第一步,并且也是长期需要做的一步,因为数据总是会出现各种各样的问题。即使数据没有问题,也不代表最终的模型效果达标。到了2016年10月中旬的时候,我俩的压力也已经很大,每天晚上都在万利达五楼加班到十点以后。每天我都在训练模型和数据分析,Jett 每天也在做平台开发的工作,虽然我不懂他在做什么,但是总是感觉很忙的样子。而且项目经理每两天都会催一次进度,顺带着会有各种大大小小的会议。其实当时 Jett 是作为业务部门的人,压力更大的是压在他那一边。在思考了大半个月之后,在一个星期日的夜晚,我去华润万家买了一堆小本子,把整个项目的细节仔仔细细地思考了一遍,突然灵光一现,发现其实有一个优化点,于是就顺手写在了小本子上。到了周一的时候,我把小本子上面的想法去实现了一遍,结果效果瞬间提升,业务指标瞬间完成。当我把这个好消息告诉 Jett 的时候,Jett 还在 RTX 上说了一句:“等我去哭一会”。最终我俩猛加班两天,把模型往线上一扔,效果直线上升,而此刻的报表系统早已齐备,瞬间就看到了效果的提升,甩掉前面的模型五条大马路。

6EC4948295154B525FD67EAC9B309320.jpg
当年撰写思路的小本子

在2017年初的时候,模型的整体效果已经很好了,整体的路线其实已经基本走通。因此后续更多的是平台建设方面的工作,这一块做得更多一点的是 Jett,最终他还真是完成了无算法人员参与就能够自动接入各种业务的全自动流程。不过随着时间的流逝,我俩都各自都被抽调去做其它的项目,合作的时间和机会在项目成功之后就变得越来越少了。万幸的是,在2017年9月5日-2017年9月9日,我俩共同去了台湾进行了一次团建。参见《五日台湾行》。

640

640-2
2017年的台湾

后来,由于种种原因,我也不再继续从事安全方面的工作,与 Jett 的交集也越来越少。虽然后面也有吃饭聚餐等时候,但是却没有了在一个项目中一起为了一个目标奋斗的机会,心中不免有少许遗憾。随着去年的大调整,Jett 和我已经不在同一个部门。今年 Jett 由于家庭原因离开了公司,恐怕短期内已经再也没有通力协作,互帮互助的机会了。也许这个项目是我俩共同做的最后一个项目,但是当年在万利达五楼共同奋斗的时光,却是一段珍贵的回忆。

zr9558's Blog