SICP Series (1) - Square Roots with Newton's Method

Newton’s Iteration and Square Roots

Principle excerpted from: Implementing the square root function with Newton’s iteration

Newton’s iteration:

$x_{n+1}= x_n-f(x_n)/f’(x_n)$

Suppose the input number is S and the square root sought is x, satisfying $S=x^2$. Then we can define the function $f(x)=x^2−S$, and the problem finally transforms into finding the root of the equation $f(x)=0=x^2−S$.

$x_{n+1}=x_n-(x_n^2-S)/ 2x_n$

After simplification it becomes:

$x_{n+1}=(x_n^2+S)/2x_n$

The Function Derivative Problem

But we can also skip the simplification and let the program compute $f’(x_n)$ for us. Before encountering SICP I had never even thought of such an operation. Now the problem transforms into: given $f(x)$, how to obtain $f’(x)$. We know mathematically $f’(x)$ is defined as:

$f’(x) = \frac{f(x+\Delta x) - f(x))}{\Delta x}$

Then based on $f(x)$ we can define $f’(x)$ — a function that accepts a value and returns the derivative value of the original function at that value. This is the definition in the scheme language:

(define dx 0.00001)
(define DERIV
  (lambda (f)
    (lambda (X)
      (/
        (- (f (+ X dx)) (f X))
        dx
      )
    )
  )
)

Fixed Points

The formula we got here is:

$x_{n+1}=x_n-(x_n^2-S)/ 2x_n$

If we treat $x_{n+1}$ as y and $x_n$ as x, the expression above is finding the fixed point of the formula below:

$f(x) = x - (x^2-S)/2x$

The result we want is that after many iterations, $x_{n+1}$ approximately equals $x_n$; therefore what we seek is exactly this iteration function’s fixed point $f(x)=x$. We can extract a tool for finding an iteration function’s fixed point, defined as follows:

(define (fixed-point f first-guess)
  (define tolerance 0.00001)
  (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)
)

Final Code

(define (newton f guess)
  (define DF (DERIV F))
  (fixed-point
    (lambda (x) (- x (/ (f x) (DF x))))
    guess
  )
)
(define (sqrt A)
  (newton
    (lambda (Y) (- A (square Y)))
    1
  )
)
(sqrt 1234)

Summary

I think the approach worth borrowing is: after giving $f(x)$, being able to use (DERIV f) to produce $f’(x)$ — DERIV accepts a function and returns that function’s derivative. Also the multiple abstractions of the execution process, e.g. the fixed-point and newton functions: they abstract the concrete process and give it a name, meaning we can use these functions to find any function’s fixed point, or use Newton’s iteration to find zeros.

Appendix: JS Implementation

It feels like many return statements must be written — maybe not as elegant as the original.

const dx = 0.00001;
function deriv(f) {
  return (x) => {
    return (f(x+dx)-f(x))/dx;
  }
}

function fixed_point(f, first_guess) {
  const tolerance = 0.00001;
  const close_enough = function(v1, v2) {
    return Math.abs(v1-v2) < tolerance;
  };
  const try_guess = function(guess) {
    const next_guess = f(guess);
    if (close_enough(guess, next_guess)) {
      return next_guess;
    }
    return try_guess(next_guess);
  };
  return try_guess(first_guess);
}

function newton(f, guess) {
  const DF = deriv(f);

  return fixed_point(
    (x) => {return x - f(x)/DF(x)},
    guess,
  )
}

function sqrt(A) {
  return newton(
    (x) => (A - x*x),
    1
  )
}
sqrt(1234);

// 35.12833614050059