Complete before class on October 13] Assigned material 6

Outline

Root finding and optimization

In this module, we will learn about how to use python to find the root of a function, and solving optimization problems. Both of them are important tools in economic analysis.

Section 1. Root finding with Newton-Raphson method

1.1. Define a function

1.2. Plot the function

1.3. The Newton-Raphson algorithm

Section 2. Optimization

Section 3. The SciPy library

3.1. Introduction to SciPy library

3.2. Root finding using SciPy

3.3. Minimization with SciPy

3.4. Maximization with SciPy

=========================================================================================================

Section 1. Root finding

Mathematically, finding the root of a function $LaTeX: f\left(x\right)$ is to find the value of x such that $LaTeX: f\left(x\right)=0$ . Let's first define an example function to do this. The function we are looking at is $LaTeX: x^4\:-\:3x^2$ . As a starting example, the analytic solution of this function is not difficult. In practice, we often have to deal with much more complex functions that are hard to solve with a pencil. Python can be very helpful in solving this type of problems.

1.1 Define a function

In[1]:

import numpy as np

def f(x):
    return x**4 - 3*x**2

1.2. Plot the function

It is always a good idea to plot the function out and see what it looks like visually. To plot the function, let's use the matplotlib library.

In[2]:

import matplotlib.pyplot as plt

# create 1000 evenly spaced points between -2 and 2, and let it be the domain of the function.
x = np.linspace(-2, 2, 1000)

# evaluate the function at x values
y = f(x)

# plot the function
plt.style.use('seaborn-whitegrid')
fig = plt.figure()
ax = fig.add_subplot(1 ,1, 1)
ax.set.title("example function")
ax.hlines(0, -2, 2, linestyle='--', color='r')  # highlight the horizontal line y = 0
plt.show()

Figure 1. The Example Function

1.3. The Newton-Raphson algorithm

As we can see from the function plot above, the function has more than one roots. Such exploratory analysis is definitely helpful. Let's be rigorous and find the numerical solutions. We can use the Newton-Raphson algorithm for root finding. The figure below illustrates the underlying logic of the Newton-Raphson algorithm.

Newton_Raphson algorithm

Remember, our goal is to find the value of LaTeX: x such that $LaTeX: f\left(x\right)\:=\:0$ . We start with an initial guess LaTeX: x_0 . The Slope of the tangent line through the initial guess can be defined as $LaTeX: f\left(x_0\right)\:-\:0$ over $LaTeX: x_0\:-\:x_1$ , or:

$LaTeX: slope\:=f\:'\left(x\right)\:=\:\:\frac{f\left(x_0\right)\:-\:0}{x_0\:-\:x_1}$

where LaTeX: x_1 is the point where the tangent line cross the x-axis.. This point has coordinates and LaTeX: y_1=0 .

Since this is a tangent line, we can also write the slope of this line as the derivative of function $LaTeX: f\left(x\right)$ , denote as $LaTeX: f\:'\left(x\right)$ . Then, we can use this equation to solve for the value of LaTeX: x_1 :

$LaTeX: x_1\:=\:x_0\:-\:\frac{f\left(x_0\right)}{f\:'\:\left(x_0\right)}$

Next, we put a new tangent line through the value of $LaTeX: f\left(x_1\right)$ and search where this tangent line crosses the a-axis (i.e., LaTeX: x_2 )

Similar to the previous steps, we can solve for the value of LaTeX: x_2 :

$LaTeX: x_2\:=\:x_1\:-\:\frac{f\left(x_1\right)}{f\:'\:\left(x_1\right)}$

Next, we will repeat the above steps using the iterative procedure:

$LaTeX: x_{n+1}\:=\:x_n\:-\:\frac{f\left(x_n\right)}{f\:'\:\left(x_n\right)}$

In order for this function to work well, we need to be able to calculate the derivative of the function around the root point and preferably over the entire domain of the function. Our target is to find the LaTeX: x that letting $LaTeX: f\left(x\right)$ is very, very, very close to 0. To apply the Newton-Raphson algorithm, we need to first find the derivative of the function. In this example, it should be

$LaTeX: f\:'\left(x\right)\:=\:4x^3\:-\:6x$

In[3]:

def fp(x):
    return 4*x**3 - 6*x

Let's then define the Newton-Raphson algorithm

In[4]:

def newton_raphson(f,fp, initial_guess, tolerance=1e-9 , max_iteration=100):
    """
    This function will apply the Newton-Raphson method to find the root of a given function
    Parameters:
    ------------
    f : function
       The original function we wanted to find roots
    fp : function
        The first order derivative of the original function
    initial_guess: list
        A list of starting points. 
    tolerance : float, optional
        defines how close to zero needs to be
    max_iteration : int
        defines maximum iterations if not converge
    Return:
    roots : list
        a list of roots found
    """
    roots = [ ]

    for x0 in initial_guess:
        x = x0
        fx = f(x)
        fpx = fp(x)
        iteration = 0
        # continue the iteration until stopping conditions are met
        while (abs(f(x)) > tolerance) and (iteration < max_iteration):
            x = x - fx/fpx
            fx = f(x)
            iteration += 1
    
        if abs(f(x)) < tolerance:
            roots.append(np.round(x,3))
            
    return roots

There are two caveats of using the Newton-Raphson method for root finding.

First, it requires that the function is differentiable. For non-differential functions, this algorithm won't work.
Second, this method is very sensitive to the starting value. Some function has more than one roots (such as the example above). In order to find all of the roots, we will need to choose several different starting points. In addition, if the starting point is very far from the root, the algorithm will converge slowly.

Based on the function plotting, we know the approximate positions of the roots. Let's use the function we defined above to find the roots.

In[5]:

results = newton_raphson(f,fp,[-2, -1, 0, 1, 2])
print(set(results)) # return a set of unique elements

Out[5]:

{0.0, 1.732, -1.732}

=========================================================================================================

2. Optimization

Optimization means find the maximum or minimum values of a function, or finding when these max/min occurs. It is an important tool in economic analysis. Think about some of the economic problems we have learned: from utility maximization to profit maximization, to cost minimization; mathematically, all of these are optimization problems.

In order to solve an optimization problem, a basic method is to look for the root of the first-order derivative so that $LaTeX: f\:'\left(x\right)\:=\:0$ .

Still using the example function above: from the previous section, we have figured out the first derivative of $LaTeX: f\left(x\right)$ . Again, it is always a good idea to plot the function first. When we have both $LaTeX: f\left(x\right)$ and $LaTeX: f\:'\left(x\right)$ on the figure, the relationship becomes more clear.

function and first order derivative

In[6]:

fig = plt.figure()
ax1 = fig.add_subplot(2 ,1, 1)
ax1.plot(x,f(x))
ax1.hlines(0, -2, 2, linestyle='--', color='r')
ax1.set_title("f (x) ")

ax2 = fig.add_subplot(2, 1, 2)
ax2.plot(x,fp(x))
ax2.hlines(0, -2, 2, linestyle='--', color='r')
ax2.set_title("f'(x)")
plt.show()

To find the roots of the first derivative, we would need to figure out the second-order derivative of $LaTeX: f\left(x\right)$ . In this example, it should be

$LaTeX: f\:''\left(x\right)\:=\:12x^2\:-\:6$

Then, we can apply the Newton-Raphson method to find the roots.

In[7]:

# find the roots of the first-order derivative
def fpp(x):
    return 12*x**2 - 6

opt = newton_raphson(fp,fpp,[-2, -1.5, -1, -0.5, -0.05, 0.05, -.5, 1, 1.5, 2])
print(set(opt))

Out[7]:

{0.0, 1.225, -1.225}

=========================================================================================================

3. The SciPy library

3.1. Introduction to SciPy library

As an open-source language, there are many high-quality libraries that can help us do our job faster and easier. We have used numpy to perform scientific calculations, random to generate random numbers, matplotlib to visualize data. Today, we will learn about another very helpful package: SciPy (pronounced "Sigh Pie"). SciPy library can help us performing jobs in optimization, linear algebra, signal and image processing, and many other tasks common in science and engineering.

3.2 Root finding using SciPy

From section 1 and section 2, we have learned about how to do root finding and optimizations. The procedure was not insanely complex, however, it was not insanely straightforward either. It would be nice if we don't have to define the root-finding algorithms ourselves, or to calculate the first order and second order derivatives manually. The scipy.optimize module can help us achieve that.

As all other libraries and modules, to use the SciPy library, the first thing is to import it to python. Instead of importing the whole library, we can only import the necessary modules from the library. In this practice, we will need the two modules:

In[8]:

from scipy.optimize import root
from scipy.optimize import minimize

The root function takes two required arguments: the first one is the function. In our example, it should be the function f we defined above. The second one is the starting point. Since the searching process will terminated once it finds a root. we still need to specify several different starting points if multiple roots exist.

In[9]:

root1 = root(f,-2)
root2 = root(f, -0.05)
root3 = root(f, 2)

Let's take root1 as an example, and look at how to read the result panel.

In[10]:

print(root1)

Out[10]:

    fjac: array([[-1.]])
     fun: array([0.])
 message: 'The solution converged.'
    nfev: 9
     qtf: array([-9.39692768e-11])
       r: array([10.39230675])
  status: 1
 success: True
       x: array([-1.73205081])

In the result panel, we should first read the message on the third row. Here, it says 'The solution converged'. It means we find a root! Next, we can jump to the last row, which says x: array([-1.73205081]). This is the root we just found. For now, we can skip all other fields in the panel. As a reminder, this is the result panel of root1. We need to check the other two result panels (root2 and root3) to see the other roots.

3.3. Minimization with SciPy

The process of find the minimum of a function is similar to root finding. the minimize function also takes two required arguments: the function need to be minimized, and the starting point. In our example function $LaTeX: x^4\:-\:3x^2$ , we know that there are two local minimum points within the domain [-2, 2]. To find both of them, we will need to run the minimization process twice, with different starting points.

In[11]:

neg_min = minimize(f,  -1.5)
pos_min = minimize(f, 1)

Let's print out neg_min and see what's in there.

In[12]:

print(neg_min)

Out[12]:

      fun: -2.2499999999991216
 hess_inv: array([[0.08374484]])
      jac: array([4.73856926e-06])
  message: 'Optimization terminated successfully.'
     nfev: 12
      nit: 4
     njev: 6
   status: 0
  success: True
        x: array([-1.22474449])

From the result panel, the message we have is 'Optimization terminated successfully'. It means we have found the minimum of the function on the domain. If we look at the top line of the result (fun: -2.2499999999991216), that is the minimum value we just found. The last row of the result panel x: array([-1.22474449]) tells that when x is at -1.22474449, the function is minimized.

3.4. Maximization with SciPy

Maximization problem is very common in economic analysis. For example, in consumer theory, we need to figure out how the consumers maximize utility; In firm and production theory, we need to calculate how a firm maximize profits.

Essentially the process of maximization are the same as minimization. However, the SciPy library does not have a "maximization" module. To solve the maximization problem with SciPy, we can simply take the negative value of the objective function, and then use the minimize module to perform the task.

Still using our example function, function here. Let's define a new function that takes the negative of the original function.

In[13]:

def neg_f(x):
    return -f(x)

Let's plot this function out and see what it looks like.

Negative function

In[14]:

fig = plt.figure()
ax = fig.add_subplot(1 ,1, 1)
ax.plot(x,neg_f(x))
ax.hlines(0, -2, 2, linestyle='--', color='r')
ax.set_title("negative function")
plt.show()

As you can see, finding the local maximum of the original function also means finding the minimum of the negative function. We can then use the SciPy to do the job.

In[15]:

max_point = minimize(neg_f,-0.05)
print (max_point)

out[15]:

      fun: 5.823264206482406e-12
 hess_inv: array([[0.16694958]])
      jac: array([8.40407969e-06])
  message: 'Optimization terminated successfully.'
     nfev: 8
      nit: 2
     njev: 4
   status: 0
  success: True
        x: array([1.39322937e-06])

=========================================================================================================