walker's code blog

coder, reader

二分法、牛顿法和梯度下降法开根号和解方程

二分法

二分法的思路是每次排除一半样本的试错方法，把样本一分为二（A和B），那么目标值不在A就在B里，不断缩小范围。

就像在玩一个猜价格的游戏，通过告诉你猜高了还是低了，你总能猜到正确价格一样，设计好一个计算差值的函数能大体上判断出你下一轮尝试的值应该在前一半还是后一半，总能迭代到足够接近的结果。

对于求平方根来说，我们没什么过多的设计，直接对中值取平方，高了就取小的一半，低了就取大的一半，实测小的数字是没问题的，这里仅仅用来演示思路。

import math

def binary_sqrt(n):
    epsilon = 1e-10         # quit flag
    start = 0
    end = n
    mid = start + (end - start) / 2
    diff = mid ** 2 - n
    while abs(diff) >= epsilon:
        # 值过大则尝试小的一半，否则就尝试大的一半，修改边界值即可
        if diff > 0:
            end = mid
        else:
            start = mid
        mid = start + (end - start) / 2
        diff = mid ** 2 - n
    return mid

for i in range(1,11):
    print(f'estimated:\t{binary_sqrt(i)}, \t sqrt({i}): \t {math.sqrt(i)}')

output:

estimated:	0.9999999999708962, 	 sqrt(1): 	 1.0
estimated:	1.4142135623842478, 	 sqrt(2): 	 1.4142135623730951
estimated:	1.7320508075645193, 	 sqrt(3): 	 1.7320508075688772
estimated:	2.0000000000000000, 	 sqrt(4): 	 2.0
estimated:	2.2360679775010794, 	 sqrt(5): 	 2.23606797749979
estimated:	2.4494897427794060, 	 sqrt(6): 	 2.449489742783178
estimated:	2.6457513110653963, 	 sqrt(7): 	 2.6457513110645907
estimated:	2.8284271247393917, 	 sqrt(8): 	 2.8284271247461903
estimated:	2.9999999999890860, 	 sqrt(9): 	 3.0
estimated:	3.1622776601579970, 	 sqrt(10): 	 3.1622776601683795

牛顿法

我就算不画图也能把这事说清楚。

牛顿法用的是斜率的思想，对$f(x)=0$，选一个足够接近目标值($x$)的点($x_0$)，计算其切线与X轴的交点($x_1$），这个交点$x_1$往往比$x_0$更接近$x$，数次迭代后肯定越来越接近目标值$x$。

  1. 问题转化成一个求函数上任一点($x_n$)的切线与X轴的交点($x_{n+1}$)的问题(我们假设n+1在n的左边，即向左来逼近$x_0$)
  2. 令$\Delta x = x_n - x_{n+1}, \Delta y = f(x_n) - f(x_{n+1})$，则$f'(x_n) = 该点斜率 = \frac{\Delta y}{\Delta x}$
  3. 展开:$f'(x_n) = \frac{f(x_n) - f(x_{n +1})}{x_n - x_{n+1}}$
  4. $\because f(x_{n+1})=0\ \Rightarrow x_{n +1} = x_n - \frac{f(x_n)}{f'(x_n)}$
  5. 至此，我们用$x_n$推出了一个离$x_0$更近的点：$x_{n+1}$
  6. 如此往复则可以推出足够精度的解。

而求任意正整数$a$的平方根，

  1. 函数就变成了 $g(x) = a, \Rightarrow g(x) = x^2$，
  2. 那么我们有: $f(x) = g(x) - a = 0 = x^2 - a = 0$
  3. $f'(x) = 2x$
  4. $f(x),f'(x)$都有了，就能代入上面得到的公式：$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$了
  5. 得到$x_{n+1} = x_n - \frac{x_n^2-a}{2x_n}$

现在可以写代码了，不断去迭代，求下一个$x_{n+1}$:

def newton_sqrt(n):
    x_n = n / 2
    epsilon = 1e-10             # quit flag

    f_x = lambda a : a**2 - n   # f(x)=x^2 - a
    df_x = lambda a : 2*a       # derivative of f(x)
    x_next = lambda a : a - f_x(a) / df_x(a)

    while abs(x_n ** 2 - n) > epsilon:
        x_n = x_next(x_n)
    return x_n

for i in range(1, 10):
    print(f'sqrt({i})\t{newton_sqrt(i)}')

output

sqrt(1)	1.000000000000001
sqrt(2)	1.4142135623746899
sqrt(3)	1.7320508075688772
sqrt(4)	2.0
sqrt(5)	2.23606797749979
sqrt(6)	2.4494897427831788
sqrt(7)	2.6457513110646933
sqrt(8)	2.8284271247493797
sqrt(9)	3.0

梯度下降法

梯度下降法的数学原理是$f(x_1,x_2,\dots$)的gradient（$\nabla f$）就是其最陡爬升方向（steepest ascent）。

可以拿这个当成结论，也可以去感性认识，而要去证明的话，网上有数不清的教程，在花书(《Deep Learning深度学习》)和可汗学院里，都是用的directional derivate来解释（证明）的，即”自定义方向上的瞬时变化率“，也是我认为在如果有多变量微积分的基础下，比较容易让人接受且简单直白的解释：

  • $\overset{\rightarrow}{v} \cdot \nabla f = |\overset{\rightarrow}{v}|\cdot|\nabla f|\cdot cos\theta$
  • $\overset{\rightarrow}{v}$ 就是指的任意方向，如果是x, y等方向，那就是普通的偏导了。
  • 显然上式当$\theta = 0$时拥有最大值，即$\overset{\rightarrow}{v}$指向的是$\nabla f$的方向，那就是梯度的方向了
  • 所以梯度方向就是爬升最陡峭的方向

在一元方程里，”梯度“就是过某点的斜率（slope)，或者说函数的导数（derivative）。

我们要到局部最小值，显然就应该向相向的方向走。并且由于越接近目标值（谷底），斜率越小，所以即使我们选择一个固定的步长（learning rate），也是会有一个越来越小的步进值去逼近极值，而无需刻意去调整步长。

以上是思路，只是要注意它$\color{red}{并不是作用到要求的函数本身}$上去的，而是精心设计的loss，或者说diff、error函数。

其实它跟前面的二分法很类似，就是不断指导里应该在哪个区间里去尝试下一个x值，再来结果与真值的差异（而牛顿法则是直接朝着直值去逼近，并不是在“尝试“）。

二分法里我随便设计了一个判断loss的函数（即中值的平方减真值），而梯度下降里不能那么随便了，它需要是一个连续的函数（即可微分），还要至少拥有局部极小值：

我们令$e(x)$表示不同的x取值下与目标值$Y$的差的平方（损失函数loss），既然是一个简单二次函数，就能求极值，且它的最小值意味着当x值为该值时估算原函数$f(x)=Y$的误差最小，有：

$e(x) = \frac{1}{2}(f(x) - Y)^2$ (1/2的作用仅仅是为了取导数时消除常数项，简化计算)
$e'(x) = (f(x) - Y) \cdot f'(x) = \Delta y \cdot f'(x)\quad \color{green}{\Leftarrow Chain\ Rule}$
$\Delta x = e'(x) \cdot lr = \Delta y \cdot f'(x) \cdot lr\ \color{red}{\Leftarrow这里得到了更新x的依据}$
$x_{n+1} = x_n - \Delta x = x_n - \Delta y \cdot f'(x) \cdot lr \Leftarrow 公式有了$

这时可以写代码了

def gradient_sqrt(n):
    x       = n / 2       # first try
    lr      = 0.01        # learning rate
    epsilon = 1e-10       # quit flag

    f_x     = lambda a : a**2
    df_dx   = lambda a : 2*a
    delta_y = lambda a : f_x(a) -n
    e_x     = lambda a : delta_y(a)**2 * 0.5     # funcon of loss
    de_dx   = lambda a : delta_y(a) * df_dx(a)   # derivative of loss
    delta_x = lambda a : de_dx(a) * lr

    count   = 0
    while abs(x**2 - n) > epsilon:
        count += 1
        x = x - delta_x(x)
    return x, count

for i in range(1, 10):
    print(f'sqrt({i}): {gradient_sqrt(i)}次')

output

sqrt(1): (0.9999999999519603, 593)次
sqrt(2): (1.4142135623377403, 285)次
sqrt(3): (1.7320508075423036, 181)次
sqrt(4): (2.0, 0)次
sqrt(5): (2.236067977522142, 103)次
sqrt(6): (2.449489742798969, 87)次
sqrt(7): (2.645751311082885, 73)次
sqrt(8): (2.828427124761154, 63)次
sqrt(9): (3.00000000001166, 55)次

Bonus

梯度下降解二次方程

  • 求解方程：$(x_1 - 3)^2 + (x_2 + 4)^2 = 0$的根

$f(x) = (x_1 - 3)^2 + (x_2 + 4)^2 = 0$

$e(x) = \frac{1}{2}(f(x)-Y)^2$

$\frac{\partial}{\partial x_1}e(x)=(f(x)-Y)\cdot(f(x) -Y)' = (f(x)-Y)\cdot\frac{\partial}{\partial x_1}((x_1 - 3)^2 + (x_2 + 4)^2-Y)$

$\therefore \begin{cases} \frac{\partial}{\partial x_1}e(x)=\Delta y \cdot 2(x_1 - 3) \ \frac{\partial}{\partial x_2}e(x)=\Delta y \cdot 2(x_2 + 4) \end{cases} $

def gradient_f(n):
    x1, x2  = 1, 1        # first try
    lr      = 0.01        # learning rate
    epsilon = 1e-4        # quit flag

    f_x     = lambda x1, x2 : (x1-3)**2 + (x2+4)**2
    dfx1    = lambda x : 2 * (x - 3)
    dfx2    = lambda x : 2 * (x + 4)
    delta_y = lambda x1, x2 : f_x(x1, x2) - n
    e_x     = lambda x1, x2 : delta_y(x1, x2)**2 * 0.5     # cost function
    dedx1   = lambda x1, x2 : delta_y(x1, x2) * dfx1(x1)   # partial derivative of loss \
    dedx2   = lambda x1, x2 : delta_y(x1, x2) * dfx2(x2)   # with Chain Rule
    delt_x1 = lambda x1, x2 : dedx1(x1, x2) * lr
    delt_x2 = lambda x1, x2 : dedx2(x1, x2) * lr

    count   = 0
    while abs(f_x(x1, x2) - n) > epsilon:
        count += 1
        x1 -= delt_x1(x1, x2)
        x2 -= delt_x2(x1, x2)
    return x1, x2, count

a, b, c = gradient_f(0)
print(f'''
a \t= {a}
b \t= {b} 
f(a, b) = {(a-3)**2 + (b+4)**2}
count \t= {c}''')

output

a 	= 2.9967765158140387
b 	= -3.9905337923806563 
f(a, b) = 9.999993698966316e-05
count 	= 249990

之所以做两个练习， 就是因为第一个是故意把过程写得非常详细，如果直接套公式的话，而不是把每一步推导都写成代码，可以简单很多（其实就是最后一步的结果）:$\Delta x = \Delta y \cdot f'(x) \cdot lr$

梯度下降解反三角函数

  • 求解arcsin(x)，在$x = 0.5$和$x = \frac{\sqrt{3}}{2}$的值

即估算两个x值，令$f(x)=sin(x)=0.5$和$f(x)=sin(x)=\frac{\sqrt{3}}{2}$
这次不推导了，套一次公式吧$\Delta x = \Delta y \cdot f'(x) \cdot lr$

import math

def arcsin(n):
    x       = 1           # first try
    lr      = 0.1        # learning rate
    epsilon = 1e-8        # quit flag

    f_x     = lambda x : math.sin(x)
    delta_y = lambda x : f_x(x) - n
    delta_x = lambda x : delta_y(x) * math.cos(x) * lr

    while abs(f_x(x) - n) > epsilon:
        x -= delta_x(x)

    return math.degrees(x)

print(f'''sin({arcsin(0.5)}) ≈ 0.5
sin({arcsin(math.sqrt(3)/2)}) ≈ sqrt(3)/2
''')

output

sin(30.000000638736502) ≈ 0.5
sin(59.999998857570986) ≈ sqrt(3)/2
walker / posts / #二分法 #牛顿法 #梯度下降 / Views