Root finding: bisection and Newton

Solving f(x) = 0

Many quant questions come down to finding the x where a function crosses zero: the yield that prices a bond, the volatility that matches an option, the rate that makes a project's value zero. A root of f is an x with f(x) = 0. Usually there is no formula, so a program closes in on it. Two methods do almost all the work. Bisection is slow but certain; Newton's method is fast but needs care. Both receive the function itself as an argument, which Python allows because functions are values.

def bisection(f, a, b, tol=1e-10):
    if f(a) * f(b) > 0:
        raise ValueError("f(a) and f(b) must have opposite signs")
    while b - a > tol:
        mid = (a + b) / 2
        if f(a) * f(mid) <= 0:       # the root is in the left half
            b = mid
        else:
            a = mid
    return (a + b) / 2

root = bisection(lambda x: x * x - 2, 0, 2)
print(root)                          # 1.414213562355144
print(root * root)                   # 1.9999999999492266

Bisection: slow and sure

If f is continuous and f(a) and f(b) have opposite signs, a root lies between them. Evaluate the midpoint: one of the two halves still has opposite signs at its ends, so the root is trapped there. Each step halves the interval, whatever f looks like, so the error after n steps is (b - a) / 2**n and you can say beforehand how many steps you need: ceil(log2((b - a) / tol)). That guarantee, with no derivative required, is why bisection is the safe default.

# How many halvings does a given tolerance cost? About 3.3 per decimal digit
import math

def bisection_steps(a, b, tol):
    if b - a <= tol:
        return 0
    return math.ceil(math.log2((b - a) / tol))

for tol in (1e-2, 1e-4, 1e-8, 1e-12):
    print(tol, bisection_steps(0, 2, tol))       # 1e-2: 8, 1e-4: 15, 1e-8: 28, 1e-12: 41

Newton: follow the tangent

Newton's method replaces f by its tangent line at the current guess and jumps to where that line crosses zero: x_new = x - f(x) / df(x). Near a simple root the error is roughly *squared* each step, so the number of correct digits doubles: from 1 to 2, 4, 8, 16 in four steps. It needs the derivative df (for implied volatility that is vega, which you already have), and a good start.

# Newton for sqrt(2): print every iterate and the error that remains
import math

def newton_step(f, df, x):
    return x - f(x) / df(x)

f = lambda x: x * x - 2
df = lambda x: 2 * x
x = 1.0
for step in range(6):
    print(step, x, abs(x - math.sqrt(2)))
    x = newton_step(f, df, x)

# Newton on the implied volatility of a call, using vega as the derivative
import math

def norm_cdf(x):
    return 0.5 * (1 + math.erf(x / math.sqrt(2)))

def norm_pdf(x):
    return math.exp(-0.5 * x * x) / math.sqrt(2 * math.pi)

def bs_call(spot, strike, rate, sigma, t):
    d1 = (math.log(spot / strike) + (rate + 0.5 * sigma ** 2) * t) / (sigma * math.sqrt(t))
    return spot * norm_cdf(d1) - strike * math.exp(-rate * t) * norm_cdf(d1 - sigma * math.sqrt(t))

def vega(spot, strike, rate, sigma, t):
    d1 = (math.log(spot / strike) + (rate + 0.5 * sigma ** 2) * t) / (sigma * math.sqrt(t))
    return spot * norm_pdf(d1) * math.sqrt(t)

market = 8.0
sigma = 0.5                                              # a rough start
for step in range(1, 6):
    error = bs_call(100, 100, 0.05, sigma, 0.5) - market
    sigma -= error / vega(100, 100, 0.05, sigma, 0.5)
    print(step, round(sigma, 10))

When Newton fails

The tangent is a local picture. If the derivative is zero or tiny, the jump is huge or undefined; if the start is far from the root, the iterates can leap away or cycle (for atan(x) from x = 2 they grow without bound). So a robust solver guards the division, limits the number of iterations, and raises an error rather than return a wrong number. Many libraries start with a bracket and use Newton inside it, falling back to bisection when a step leaves it.

• A root of f is an x with f(x) = 0; solvers take f (and df) as function arguments.
• Bisection: keep the half with a sign change; ceil(log2((b - a) / tol)) steps, guaranteed, no derivative.
• Newton: x - f(x) / df(x); digits double each step near a simple root, but it needs df, a good start, and guards for a zero derivative and for no convergence.

Next up: binomial trees, pricing an option by stepping the stock up and down and working backwards.

Lab · what's graded

• Trap a root by halving with bisection25 pts
• Count the halvings with bisection_steps25 pts
• Follow the tangent once with newton_step25 pts
• Converge on a root with newton25 pts