SICP开坑——使用Scheme搜索不动点

最近去读了SICP的第一章，初次领悟到了函数式编程的魅力。曾经尝试过好几次入坑，但是都是跟着网上奇奇怪怪的教程入的，不久就放弃了。这次要一鼓作气读完这本大部头！flag立起来了

对于求一些式子的值，比如$\sqrt{x}$，曾经在课上学过牛顿法等去计算，但是一直又都没有真正理解为什么这样去计算，在SICP中，介绍了一种寻找不动点来计算这类值的方法。所谓不动点，就是当函数的输入等于输出时的一个点，这个时候有$f(x) = x$，以刚才的$y = \sqrt{x}$举例子，可以有以下推导：

$y = \sqrt{x}\\ y^2 = x \\ y = \frac{x}{y}$

如果设$f(y) = \frac{x}{y}$，就可发现，计算$\sqrt{x}$的值就等价于寻找$f(y)$这个函数的不动点，只要找到一个不动点，就可以确定这个开2次方的结果。

使用普通的迭代算法，即首先随意选择一个猜测值，然后带入函数中，如果不是不动点就将函数的输出作为输入再代入一次，直到找到一个不动点。

但是新的问题又来了，怎么保证这样做下去一定会找到一个不动点呢？关于这一点，SICP没有告诉我们，因为涉及到了大量的数学知识。但是我发现了这个定理

设$y_0$为$f$的一个不定点，当且仅当$|f'(y_0)| < 1$时，在$y_{0}$附近不动点有局部收敛性，也就是说，要保证这个方法收敛，就必须保证$|f'(y_0)| < 1$。如果将之前的函数求导，会发现是等于1的，并不小于1，这种方法不能收敛，SICP 给的方法是对原方程进行一些变形，每次得到的不是$f(y)$，而是$f(y)$和$y$的平均值，这种方法叫做平均阻尼，也就是变成下面这种样子：

$y = \frac{1}{2}(\frac{x}{y} + y)$

这样的话得到了一个不一样的函数，计算出导数就是小于1的了。

SICP 在章节的最后留了一个问题，这样的变形对于求2次方是可以的，但是求4次方时就无能为力了，但是可以做两次平均阻尼来完成，也就是变成：

$y = \frac{1}{2}(\frac{1}{2}(\frac{x}{y} + y) + y)$

可以证明这个方程仍然和原方程相等，但是是在求n次方根的时候，要做多少次平均阻尼才能收敛呢，虽然 SICP 让做实验得出结论，但是有了上面的定理，我尝试使用公式推导出来了这个结果，以下是推导过程：

设求的是n次方根，需要做a次平均阻尼才能让这个函数收敛，不妨先观察 a=1,2,3 时候这个公式的变形情况

a = 1时

$y = \frac{1}{2}(\frac{x}{y^{n-1}}+y)\\\\ y = \frac{x}{2y^{n-1}}+\frac{y}{2}$

a = 2时

$y = \frac{1}{2}((\frac{x}{2y^{n-1}} + \frac{y}{2}) + y)\\\\ y = \frac{x}{2^2y^{n-1}} + \frac{y(1+2)}{2^2}$

a = 3时

$y = \frac{1}{2}((\frac{x}{2^2y^{n-1}} + \frac{y(1+2)}{2^2}) + y)\\\\ y = \frac{x}{2^3y^{n-1}} + \frac{y(1+2+2^2)}{2^3}$

通过对上面几种情况的归纳，可以得出一般的规律：

$y = \frac{x}{2^ay^{n-1}} + \frac{y(1+2^1+2^2+...+2^{a-1})}{2^a}\\\\ y = \frac{x}{2^ay^{n-1}} + \frac{y(2^a-1)}{2^a}$

这个时候即可设

$f(y) = \frac{x}{2^ay^{n-1}} + \frac{y(2^a-1)}{2^a}$

对这个函数求导，可以得到

$\frac{df}{dy} = (1-n)\frac{x}{2^ay^n} + \frac{2^a-1}{2^a}$

可以知道，在$f(y)$的不动点处，$y^n = x$，将这个结论带入，并且已知导数的绝对值小于1，可以得到

$|\frac{1-n}{2^a} +\frac{2^a-1}{2^a}| < 1\\\\ |\frac{2^a-n}{2^a}| < 1\\\\ |1-\frac{n}{2^a}| < 1\\\\ |\frac{n}{2^a}-1| <1\\\\ -1 < \frac{n}{2^a}-1 < 1\\\\ 0 < \frac{n}{2^a} < 2$

由于我们知道n一定大于0，所以可以省略大于0的条件，得到下面的式子

$\frac{n}{2^a} < 2\\\\ n < 2^{a+1}\\\\ a+1 > \log_2n\\\\ a > \log_2n-1$

真相大白了，原来做平均阻尼的次数只需要$\log_2n$次就够了，最终的练习是使用Scheme编写一个能够计算n次方根的过程，下面是代码

(define (repeated f n)
  (if (= n 1)
    f
    (lambda (x)
      (f ((repeated f (- n 1)) x)))))

(define tolerance 0.00001)

(define (fixed-point f first-guess)
  (define (close-enough? v1 v2)
    (< (abs (- v1 v2)) tolerance))
  (define (try guess)
    (let ((next (f guess)))
      (if (close-enough? guess next)
        next
        (try next))))
  (try first-guess))

(define (average a b)
  (/ (+ a
        b)
     2.0))

(define (average-damp f)
  (lambda (x)
    (average
      (f x)
      x)))

(define (ilog2 x)
  (define (iter it res)
    (if (= it 1)
      res
      (iter
        (arithmetic-shift it -1)
        (+ res 1))))
  (iter x 0))



(define (pow a b)
  (cond
    ((= a 0) 0.0)
    ((= b 0) 1.0)
    (else
      (let ((res (pow a (arithmetic-shift b -1))))
        (* res
           res
           (if (= (remainder b 2) 0.0)
             1.0
             a))))))

(define (sqrtn x n)
  (fixed-point ((repeated
                  average-damp
                  (ilog2 n))
                (lambda (y)
                  (/ x
                     (pow y (- n 1)))))
               (average x
                        (pow x n))))

(sqrtn 2 12)

我在最后最终做了一个$\sqrt[12]{2}$的计算，结果是1.059465254018959，代入计算器中基本符合计算要求，可以观察得到之前的推演是正确的

总结：

之前对偏数学上理解程序很弱，需要进一步提高