说说泰勒展开、lop1p 及 expm1

在数学里,泰勒展开是拉格朗日定理(微分中值定理)和牛顿迭代法在高阶可微函数的推广。在机器学习领域,我们已经见过其威力;通过应用泰勒展开,我们可以大幅降低计算的复杂度。

这里,我们会回顾一下泰勒展开相关数学定义,并说说泰勒展开在计算机科学领域的另一个应用——高精数值计算。

泰勒其人

说起来,现在名气最大的「泰勒」,当属那位甩一个男友就写一首歌的当红女歌手(泰勒·斯威夫特)了。不过此处我们说的泰勒,是英国数学家布鲁克·泰勒(Brook Taylor)。

泰勒 1685 年生于英格兰米德萨斯郡,1731 年逝于英格兰伦敦。尽管后世尊其为数学家,但实际上泰勒在剑桥大学攻读的学位是法学学士和博士学位。泰勒留下的工作不多,一方面由于其从事数学工作的时间比同时期其他数学家要短,另一方面由于它写的文稿往往语焉不详,难以理解。由此可见,好好说话的重要性。

泰勒最著名的工作有两个,但都集中在 1715 年发表的著作《Methodus Incrementorum Directa et Inversa》当中。其一是,从这部著作出发,后世引出了高等数学的一个分支:有限差分方法。其二是,这部著作中,提出了对无穷可微函数的泰勒级数。后者被他的后辈,著名的法国数学家约瑟夫·拉格朗日(Joseph Lagrange)尊为「导数计算的基础」(le principal fondement du calcul différentiel)。

泰勒展开与麦克劳林展开

从原始的微分定义开始

按照微分定义,函数在一点的微分,是它在该点附近的最佳线性近似:

$$f(a+h) = f(a) + h\cdot f’(a) + o(h).$$

其中,$o(h)$ 表示 $h$ 的高阶无穷小。这也就是说,对于可微函数 $f(x)$ 来说,有近似

$$f(a + h) \approx f(a) + h\cdot f’(a).$$

显而易见,$f(x)$ 与 $f(a) + (x - a)\cdot f’(a)$ 在 $a$ 附近是很接近的。以数学的语言表述,这两个表达式在 $a$ 处的零阶导数(函数自身)与一阶导数相等。由此我们可以预计,如果一个函数足够光滑($n + 1$ 阶可微),那么若有一个多项式在 $a$ 处的前 $n$ 阶导数值与函数自身的导数值相等,那么这一多项式就可以是该函数的近似。

泰勒定理

定理的描述如下。

设 $n \in \mathbb{Z}^+$ 是一个正整数;设 $f(x)$ 是定义在包含 $a$ 的区间上的函数,并且在 $a$ 处 $n+1$ 阶可微。那么,对于这个函数在定义区间内的任意 $x$,都有:
$$f(x) = f(a) + \frac{f’(a)}{1!}(x - a) + \frac{f’’(a)}{2!}(x - a)^2 + \cdots + \frac{f^{(n)}}{n!}(x - a)^n + \text{Res}_n(x).$$
此处的多项式,称为 $f(x)$ 在 $a$ 处的泰勒展开式,并且 $\text{Res}_n(x)$ 作为泰勒公式的余项是 $(x - a)^n$ 的高阶无穷小。

泰勒公式的余项有多种形式,在不同意义上描述了泰勒展开与原函数的接近程度。

  • 仅只说余项是高阶无穷小,即 $\text{Res}_n(x) = o[(x - a)^n]$,则是「皮亚诺」余项;
  • 拉格朗日余项可视作微分中值定理的推广,即 $\text{Res}_n(x) = \frac{f^{(n + 1)}(\theta)}{(n + 1)!}(x - a)^{n + 1}$,其中 $\theta \in (a, x)$;
  • 积分余项则课视作积分基本定理的推广,即 $\text{Res}_n(x) = \int_a^x \frac{f^{(n + 1)}(t)}{n!}(x - t)^{n} \mathop{}\!\mathrm{d}t$。

泰勒级数

泰勒定理描述了 $n + 1$ 阶可微函数的 $n$ 阶泰勒展开。相应的,对于在 $a$ 处无穷可微的函数,则可以有在 $a$ 处的无穷阶泰勒展开。与此同时,余项因为趋于 $0$,可以取消。

$$f(x) = \sum_{n = 0}^{\infty} \frac{f^{(n)}(a)}{n!} (x - a)^n.$$

特别地,若 $a = 0$,则泰勒展开特化为麦克劳林展开。

泰勒级数与计算精度

浮点运算的精度损失

众所周知,计算机的浮点运算,由于内部表示的问题(小数的二进制表示),总是有限的。通常情况,这无关紧要,但是某些情况仍会导致不少问题。例如说,有经验的程序员都知道,不能直接用相等去比较两个不知经过何种运算的来的浮点数。因为,它们在逻辑上可能相等,但是由于多次计算的精度损失,它们在内部表示上却并不相等。

为了说明精度损失,我们看一个在机器学习任务中常见的函数。

$$f(x) = \ln(1 - \exp(-x)).$$

这个函数的定义域是 $(0, +\infty)$。

在 $x$ 取值适中的时候,我们可以直接计算函数值。

但是,当 $x \to 0^+$ 时,$\exp(-x)$ 十分接近 $1$;因此 $1 - \exp(-x)$ 损失精度之后,可能直接变成 $0$。此时,再进行对数运算,就不合适了。

又或者,当 $x$ 较大时,$\exp(-x)$ 十分接近 $0$;因此 $1 - \exp(-x)$ 损失精度之后,可能直接变成 $1$。此时,再进行对数运算,就只能得到 $0$ 了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
In [1]: import math

In [2]: def f(x):
...: return math.log(1 - math.exp(-x))
...:

In [3]: f(1)
Out[3]: -0.45867514538708193

In [4]: f(10)
Out[4]: -4.54009603704951e-05

In [5]: f(100)
Out[5]: 0.0

In [6]: f(1000)
Out[6]: 0.0

In [7]: f(0.1)
Out[7]: -2.3521684610440903

In [8]: f(0.01)
Out[8]: -4.610166019324902

In [9]: f(1e-8)
Out[9]: -18.420680750029838

In [10]: f(1e-15)
Out[10]: -34.53957599234088

In [11]: f(1e-20)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-11-f08323787079> in <module>()
----> 1 f(1e-20)

<ipython-input-2-80da3292b612> in f(x)
1 def f(x):
----> 2 return math.log(1 - math.exp(-x))

ValueError: math domain error

泰勒展开拯救世界

如果是简单的计算任务,损失一些精度可能无关紧要。但是,在一些不容出错的领域——比如利息计算,少许的精度损失,可能导致成百上千的差距。然而,精度损失是计算机系统的固有问题,我们无法轻易解决。此时,就轮到泰勒展开出场,去拯救世界了。

我们观察两个函数的泰勒展开,

$$\begin{aligned}
g(x) ={}& \ln(1 + x) = x - \frac{x^2}{2} + \frac{x^3}{3} - \cdots \\
h(x) ={}& \exp(x) - 1 = x + \frac{x^2}{2!} + \frac{x^3}{3!} + \cdots
\end{aligned}$$

并且注意到以下事实,

$$\begin{aligned}
f(x) ={}& g(-\exp(-x)) \\
f(x) ={}& \ln(-h(-x)).
\end{aligned}$$

因此,计算机在计算 $g(x)$ 和 $h(x)$ 时,可以借助泰勒展开,将对数函数、指数函数,转换成精度更高的幂函数计算,以此提高精度。另一方面,利用上述变形,我们就可以借助 $g(x)$ 和 $h(x)$ 去计算 $f(x)$ 了。特别地,在 $x \to 0^+$ 时,使用 $\ln(-h(-x))$ 的效果好;在 $x \to \infty$ 时,使用 $g(-\exp(-x))$ 的效果好。

实际实现看看

$g(x)$ 和 $h(x)$ 在计算领域是两个重要的函数,多数语言的标准库中,都提供了这两个函数的泰勒展开实现。通常 $g(x)$ 的名字叫做 log1p 而 $h(x)$ 的名字叫做 expm1。我们用 log1pexpm1 重新实现目标函数看看。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
In [12]: def f(x, boundary = math.log(2)):
...: if x <= boundary:
...: return math.log(-math.expm1(-x))
...: else:
...: return math.log1p(-math.exp(-x))
...:

In [13]: f(1)
Out[13]: -0.45867514538708193

In [14]: f(10)
Out[14]: -4.5400960370489214e-05

In [15]: f(100)
Out[15]: -3.720075976020836e-44

In [16]: f(1000)
Out[16]: -0.0

In [17]: f(0.1)
Out[17]: -2.3521684610440907

In [18]: f(1e-8)
Out[18]: -18.420680748952364

In [19]: f(1e-15)
Out[19]: -34.538776394910684

In [20]: f(1e-20)
Out[20]: -46.051701859880914

可以看到,新版本的目标函数,能够区分出自变量取值 $x = 100$ 和 $x = 1000$ 之间的差异了;另外,原本提示定义域错误的 $x = 1e-20$ 现在也能计算出正确的值了。