一些π的求解方法
Kotori Y 27 Posts
 2022-08-06 20:01 2022-08-10 22:22      

圆周率是一个数学常数，为一个圆的周长和其直径的比率，近似值约等于3.14159265，常用符号 来表示。

二十世纪中期计算机技术的发展、革新再次引发了计算π位数的热潮。美国数学家约翰·伦奇及李维·史密斯在1949年利用桌上型计算机计算到1,120位。同年，乔治·韦斯纳（George Reitwiesner）及约翰·冯·诺伊曼带领的团队利用反三角函数（arctan）的无穷级数，通过ENIAC计算到了小数第2,037位，花了70小时的电脑工作时间。这一纪录后来多次由其他透过arctan级数计算出的结果打破（1957年到7480位小数，1958年到第一万位数，1961年到第十万位小数），直到1973年，人们计算出了小数点后的第一百万位小数。下面将写几个我所知道的的一些求解方法。

1. 蒙特卡洛方法

首先是朴素的统计学方法，蒙特卡洛方法是以概率统计理论为指导的一类非常重要的数值计算方法，通过进行大量重复试验计算事件发生的频率，按照大数定律（即当试验次数充分大时，频率充分地接近于概率）可以求得的近似值。其中一种方法为：随机地往一个与边长为的正方形内接的圆内投置大量小球。最终，落在圆的小球个数与总投球数的比值即为：

代码示例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
import math
import random
import matplotlib.pyplot as plt
import numpy as np


def mc_approx(total: int = 100000) -> float:
    f, ax = plt.subplots(figsize=(8, 8))

    inner = []
    outer = []
    for _ in range(total):
        x, y = random.random(), random.random()
        if (x ** 2 + y ** 2) <= 1:
            inner.append([x, y])
            continue
        outer.append([x, y])

    ax.scatter(*np.array(inner).T, color='red', s=2)
    ax.scatter(*np.array(outer).T, color='green', s=2)
    plt.show()

    return len(inner) / total * 4  # 代码为1/4的内接圆，因此乘以4

循环10,000次所求得的圆周率为3.1536，精度较差。

2. 二分查找算法

在计算机科学中，二分查找算法（binary search algorithm），是一种在有序数组中查找某一特定元素的搜索算法。搜索过程从数组的中间元素开始，如果中间元素正好是要查找的元素，则搜索过程结束；如果某一特定元素大于或者小于中间元素，则在数组大于或小于中间元素的那一半中查找，而且跟开始一样从中间元素开始比较。如果在某一步骤数组为空，则代表找不到。这种搜索算法每一次比较都使搜索范围缩小一半。

而我们知道，在的范围内是单调递减的，并且在内大于零，内小于零。因此，可以通过二分法判断正弦函数的正负来不断逼近。

代码示例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
def bin_search(epoch: int = 100) -> float:
    lower = 3
    upper = 4

    x = (lower + upper) / 2

    for _ in range(epoch):
        if math.sin(x) > 0:
            lower = x
        elif math.sin(x) < 0:
            upper = x
        else:
            return x
        x = (lower + upper) / 2

    return x

通过循环100次可得到结果为：3.141592653589793，精度相较蒙特卡罗方法大幅提升。

3. 梯度下降法

虽然二分查找对于求解的问题已表现得很让人满意了，但对于其他不具有明显分界点或者更高维度的函数就不适用了。这里介绍一种被广泛用于求解极值的方法——梯度下降法。

3.1 导数与微分

导数（derivative）是微积分学中的一个概念。函数在某一点的导数是指这个函数在这一点附近的变化率。导数的本质是通过极限的概念对函数进行局部的线性逼近。当函数的自变量在一点上产生一个改变量时，函数值的增量与自变量的改变量的比值在趋于0时的极限如果存在（即存在），即为在处的导数，记作、或。几何意义上来说，函数在处的导数就是该函数过点的切线的斜率。

而函数的微分（Differential of a function）是指对函数的局部变化的一种线性描述。微分可以近似地描述当函数自变量的取值作足够小的改变时，函数的值是怎样改变的。微分和导数是两个不同的概念。但是，对一元函数来说，可微与可导是完全等价的概念。可微的函数，其微分等于导数乘以自变量的微分，换句话说，函数的微分与自变量的微分之商等于该函数的导数。因此，导数也叫做微商。于是函数的微分又可记作。

e.g. 1

求的近似解。

(1)

设存在函数，那么(1)可以看成在处产生一个微小的该变量

e.g. 2

比较, 与之间的大小

因此有

3.2 偏导数

偏导数（partial derivative）的定义是：一个多变量的函数(或称多元函数)，对其中一个变量(导数)微分，而保持其他变量恒定。

函数关于变量的偏导数写为或。偏导数符号是全导数符号的变体，由阿德里安-马里·勒让德引入，并在雅可比的重新引入后得到普遍接受。

e.g.

注意：对于二元函数，在处的偏导数存在，不能得出函数在处连续的结论。

3.3 方向导数

以二元函数为例，偏导数是垂直于两个坐标轴的函数的变化率，而这对于真实世界的的问题来说是不足够的，因此产生了方向导数的概念。

设是平面上以为起始点的一条射线，是与同方向的单位向量。射线的参数方程为

设函数在点的某个领域内有定义，为上的另一点，且。如果函数增量与到的距离的比值在时的极限存在，那么称此极限为函数沿着方向时的方向导数，记做，即。

3.4 梯度

上节中的公式，通过证明等于，写成向量数量积的形式为。其中为函数在处的梯度，记做或.显然，梯度是个向量。

方向导数可进一步转换为。当：

  • ：，方向导数有最大值，即函数增长最快的方向，值为梯度的模；

  • ：，方向导数有最小值，即函数减小最快的方向，值为梯度的模；

  • ：，方向导数为0，函数值不变。

那么对于某个函数极值，只需要沿着梯度的正方向或者反方向即可快速找到。

3.5 梯度下降法的应用

求

设存在损失函数,显然当时，函数有最小值

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
from typing import Callable, List

from scipy.misc import derivative


def dx(func, n):
    return derivative(func, n, dx=1e-6)


def sqrt_2(x):
    # loss = 1/2 * (x^2 - 2)^2
    return 1 / 2 * (x ** 2 - 2) ** 2


def grad_desc(func: Callable, init_solve: float = 0.0, epochs: int = 100, eta: float = 1e-1) -> List[float]:
    """
    
    Args:
        func: function, object function
        init_solve: float, initial solution
        epochs: int, number of iteration
        eta: float, learning rate

    Returns: list of float, solution of each epoch

    """
    s = init_solve
    solves = [s]

    for epoch in range(epochs):
        s += -dx(func, s) * eta
        solves.append(s)

    return solves
  
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
# 可视化
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.animation import FuncAnimation


def draw_anime(func, X, solves, outfile):
    y = [func(x) for x in X]

    fig, ax = plt.subplots(figsize=(8, 4.5))

    ax.plot(X, y)
    scat = ax.scatter(solves[0], func(solves[0]), color='m')
    label = ax.text(0.7, 0.9, f'x={solves[0]}', va='top', transform=ax.transAxes)

    def animate(i): 
        scat.set_offsets(np.c_[solves[i], func(solves[i])])
        label.set_text(f'epoch={i}\nx={solves[i]}')

    anim = FuncAnimation(fig, animate, interval=100, frames=len(solves))

    plt.draw()
    plt.show()

    anim.save(outfile)

solves = grad_desc(sqrt_2, init_solve=3.5, epochs=100, eta=1e-2)
draw_anime(sqrt_2, X = np.arange(1, 4, 0.01), solves=solves, outfile='sqrt2.gif')

求解

设存在损失函数，显然当时，函数有最小值

1
2
3
4
5
6
7
8
9
import math


def solve_pi(x):
    # loss = 1/2 * sin(x)^2
    return 1/2 * math.sin(x) ** 2
    
solves = grad_desc(solve_pi, init_solve=2, epochs=100, eta=1e-1)
draw_anime(solve_pi, X = np.arange(0, 4, 0.01), solves=solves, outfile='solve_pi.gif')

陷入局部最优

梯度下降法最大的一个问题就是容易陷入局部最优，假设存在这样一个函数，用上述方法会得到这样的结果

1
2
3
4
5
def spec_func(x):
    return (x ** 2 - 2 * np.cos(np.pi * x)) / 8
  
solves = grad_desc(spec_func, init_solve=5, epochs=100, eta=1e-1)
draw_anime(spec_func, X = np.arange(-6, 6, 0.01), solves=solves, outfile='grad1.gif')

对于这个问题，可以参考现实中小球的滚动：由于惯性的存在，小球的运动不会因为到达一个最低点而停止。因此我们可以把梯度看成点受的力：时刻实际受到的力。其中为衰减系数，为时刻的梯度。

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
def moment_grad_desc(
        func: Callable, init_solve: float = 0.0, 
        epochs: int = 100, eta: float = 1e-1, beta: float = 0.9) -> List[float]:
    """
    
    Args:
        func: function, object function
        init_solve: float, initial solution
        epochs: int, number of iteration
        eta: float, learning rate
        beta: float, decreasing rate

    Returns: list of float, solution of each epoch

    """
    v = 0
    s = init_solve
    solves = [s]

    for epoch in range(epochs):
        v = v * beta - dx(func, s) * eta
        s += v
        solves.append(s)

    return solves
  
solves = moment_grad_desc(spec_func, init_solve=5, epochs=100, eta=1e-1, beta=0.8)
draw_anime(spec_func, X = np.arange(-6, 6, 0.01), solves=solves, outfile='grad2.gif')

相比之下，动量梯度下降法的运动轨迹更加“真实”，也更加快速地找到了全局最优。

此外，还存在很多优化的方法，例如AdaGrad、RMSProp与Adam来解决不同的问题。这个就后面另开一篇来写吧。

参考文献

  1. 高等数学（第七版） - 同济大学数学系
  2. 导数 - wikipedia
  3. 微分 - wikipedia
  4. 方向导数 - wikipedia
  5. 【4K算法详解】【梯度】从二分法到梯度下降，带你领你的思维走过六十年长征！ - 哔哩哔哩
  • Post title：一些π的求解方法
  • Post author：Kotori Y
  • Create time：2022-08-06 20:01
  • Update time：2022-08-10 22:22
  • Post link：https://blog.iamkotori.com/2022/08/06/π的一些求解方法/
  • Copyright Notice：All articles in this blog are licensed under BY-NC-SA unless stating additionally.
 Comments