📖 Root finding and optimization

Two classic algorithms for solving $f(x)=0$, and what they become when the equation we are solving is the first order condition of a likelihood. Every estimator in this course sits on top of one of these.

Bisection method¶

The first of two very important classic algorithms for equation solving.

Solve equations of the form (we focus on the scalar case today)

The latter condition requires that the function $f(x)$ takes different signs at the endpoints $a$ and $b$.

The algorithm is similar to binary search, but in continuous space.

Input: function f(x)
       brackets [a,b] such that f(a)f(b)<0
       convergence tolerance epsilon
       maximum number of iterations max_iter

Algorithm:
  step 0: ensure all conditions are satisfied
  step 1: compute the sign of the function at (a+b)/2
  step 2: replace a with (a+b)/2 if f(a)f((a+b)/2)>0, otherwise replace b with (a+b)/2
  step 3: repeat steps 1-2 until |a-b| < epsilon, or max_iter number of iterations is reached
  step 4: return (a+b)/2

def bisection(f,a=0,b=1,tol=1e-6,maxiter=100,callback=None):
  '''Bisection method for solving equation f(x)=0
  on the interval [a,b], with given tolerance and number of iterations.
  Callback function is invoked at each iteration if given.
  '''
  if f(a)*f(b)>0:
    raise ValueError('Function has the same sign at the bounds')
  for i in range(maxiter):
    err = abs(b-a)
    if err<tol: break
    x = (a+b)/2
    a,b = (x,b) if f(a)*f(x)>0 else (a,x)
    if callback != None: callback(err=err,x=x,iter=i)
  else:
    raise RuntimeError('Failed to converge in %d iterations'%maxiter)
  return x

Source

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [9, 6]

f = lambda x: -4*x**3+5*x+1
a,b = -3,-.5  # upper and lower limits
xd = np.linspace(a,b,1000)  # x grid
plt.plot(xd,f(xd),c='red')  # plot the function
plt.plot([a,b],[0,0],c='black')  # plot zero line
ylim=[f(a),min(f(b),0)]
plt.plot([a,a],ylim,c='grey')  # plot lower bound
plt.plot([b,b],ylim,c='grey')  # plot upper bound
def plot_step(x,**kwargs):
    plot_step.counter += 1
    plt.plot([x,x],ylim,c='grey')
plot_step.counter = 0  # new public attribute
bisection(f,a,b,callback=plot_step)
print('Converged in %d steps'%plot_step.counter)
plt.show()

Converged in 22 steps

Bisection is slow but bulletproof: given a valid bracket it always converges, at a linear rate, halving the interval each step. Remember this when we get to solving Bellman equations — a robust but slow method makes an excellent fallback inside a poly-algorithm.

Newton–Raphson method¶

The second of the two classic methods for solving an equation $f(x)=0$; gradient based.

General form

• Equation solving

• Finding a maximum/minimum based on the FOC, in which case $f(x)=Q'(x)$

Derivation using a Taylor series expansion¶

Take the first two terms, assume $f(x)=0$ is the solution, and let $x_0=x_i$ and $x=x_{i+1}$

The main idea of the Newton–Raphson method is to iterate on this equation starting from some $x_0$

It is applicable to systems of equations, in which case $x\in\mathbb{R}^n$ and $f: \mathbb{R}^n \to \mathbb{R}^n$, with the derivative replaced by the Jacobian.

Input: function f(x)
       gradient function f'(x)
Algorithm:
1. Start with some good initial value
2. Update x using the Newton step above
3. Iterate until convergence

def newton(fun,grad,x0,tol=1e-6,maxiter=100,callback=None):
    '''Newton method for solving equation f(x)=0
    with given tolerance and number of iterations.
    Callback function is invoked at each iteration if given.
    '''
    for i in range(maxiter):
        x1 = x0 - fun(x0)/grad(x0)
        err = abs(x1-x0)
        if callback != None: callback(err=err,x0=x0,x1=x1,iter=i)
        if err<tol: break
        x0 = x1
    else:
        raise RuntimeError('Failed to converge in %d iterations'%maxiter)
    return (x0+x1)/2

Source

f = lambda x: -4*x**3+5*x+1
g = lambda x: -12*x**2+5
a,b = -3,-.5  # upper and lower limits
xd = np.linspace(a,b,1000)  # x grid
def plot_step(x0,x1,iter,**kwargs):
    plot_step.counter += 1
    if iter<5:
        plt.plot(xd,f(xd),c='red')  # plot the function
        plt.plot([a,b],[0,0],c='black')  # plot zero line
        ylim = [min(f(b),0),f(a)]
        plt.plot([x0,x0],ylim,c='grey') # plot x0
        l = lambda z: g(x0)*(z - x1)
        plt.plot(xd,l(xd),c='green')  # plot the tangent line
        plt.ylim(bottom=10*f(b))
        plt.title('Iteration %d'%(iter+1))
        plt.show()
plot_step.counter = 0  # new public attribute
newton(f,g,x0=-2.5,callback=plot_step)
print('Converged in %d steps'%plot_step.counter)

Converged in 7 steps

Newton–Raphson is fast but fragile: quadratic convergence near the root, but it can diverge, cycle, or run off to a different root if started badly or if $f'(x)$ is near zero.

Measuring the complexity of Newton and bisection methods¶

• What is the size of the input $n$?

• The desired precision of the solution!

• Thus, attention to the errors in the solution as the algorithm proceeds

• The rate of convergence is part of the computational complexity of these algorithms

Computational complexity: calculating a root of a function $f(x)$ with $n$-digit precision, provided that a good initial approximation is known, is $O\big(\log(n) F(n)\big)$, where $F(n)$ is the cost of calculating $f(x)/f'(x)$ with $n$-digit precision.

From root finding to optimization¶

Estimation problems in this course are of the form

where $q$ is the negative log-likelihood contribution, a moment condition, or some other criterion. In nonlinear models the optimum rarely has a closed form and must be found numerically — by solving the first order conditions, which is exactly the root finding problem above.

The relevant algorithms — Newton–Raphson, BHHH, BFGS — all belong to the class of quasi-Newton methods with the form

where

• $g$ is the index of the iteration

• $\lambda$ is the step size

• $\sum_i H_i(\theta^{(g)})$ is the sample sum of the Hessians of $q(w_i;\theta)$ evaluated at $\theta^{(g)}$, reflecting the curvature of $Q_N=\sum_i q(w_i;\theta)$

• $\sum_i s_i(\theta^{(g)})$ is the sample sum of the scores evaluated at $\theta^{(g)}$, measuring the slope of $Q_N$

Newton–Raphson for optimization¶

Second order Taylor expansion of the objective function around $\theta^{(g)}$:

First order condition:

Properties

• Uses slope and curvature

• Moves downhill if the Hessian is positive definite

Step size $\lambda$ and global convergence¶

In addition to the plain Newton–Raphson update, it is typical to include a step size search. This ensures global convergence of the algorithm, at a slower rate.

1. Start with the full Newton step, $\lambda=1$

2. If the step decreases the objective, proceed with $\lambda=1$

3. If the step does not decrease the objective, reduce $\lambda$ to half its value

4. Repeat steps 2–3 until an improving step is found

The step sizes are therefore in the sequence $\lambda=1,1/2,1/4,1/8,\ldots$, so the approach is often called step-halving line search.

This addition is computationally cheap because the Hessian and the scores do not have to be recomputed for different $\lambda$ values — only the function values.

Newton–Raphson works best when the objective is close to quadratic and convex. Yet the Hessian may fail to be positive definite far from the optimum, leading to uphill moves. The algorithm can be improved by ensuring that the curvature matrix is always positive definite — which is exactly what BHHH does.

BHHH¶

BHHH (Berndt, Hall, Hall and Hausman) approximates the Hessian by the outer product of the scores:

Motivation

• Information identity: $\E(s_i s_i') \to \E(H_i)$

• Derived for the ML estimator, but works more generally for M-estimators

Advantages

• Always positive definite

• No need to compute the Hessian — only first derivatives

Disadvantages

• The approximation is valid only

  • at the true parameters

  • for large $N$

  • for well-specified models

BFGS and other quasi-Newton methods¶

BFGS builds up an approximation to the inverse Hessian iteratively, ensuring positive definiteness. Other approaches to approximating the Hessian, such as the DFP (Davidon–Fletcher–Powell) method, were developed for general optimization problems.

BHHH, being based on the statistical properties of the scores, is more specific to econometric M-estimation problems — and for that reason it is often not found in standard optimization packages.