汪群超 Chun-Chao Wang

Dept. of Statistics, National Taipei University, Taiwan

Lesson 7: Root-finding and minimization of unit-variate function

Objective:

  1. Compute f(x) = 0 in both numerical and symbolic ways.
  2. Find min f(x) in both numerical and symbolic ways.
  3. Get familiar with the graphing techniques in locating the real roots and the local extrema.

Prerequisite:


範例 1:Numerical root-finding for polynomial functions

  1. Solve f(x) = x^2-3x+2=0 for x
  2. Solve f(x) = x^4-8x^3+16x^2-8x+1=0 for x
  3. Draw the graph for each function and mark the locations of the roots.

Note:
  1. Polynomial function usually receives a special treatment in locating roots. In this example, the numpy.polynomial.polynomial.polyroots is employed.
  2. In this example, we focus on the real roots and need to take out imaginary roots.
  3. For polynomial functions, all roots are found in a single command (This is not the case for the general funcion).
import numpy as np
import matplotlib.pyplot as plt 
import numpy.polynomial.polynomial as poly

coef = [2, -3, 1] # from low to high
r = poly.polyroots(coef)
print('The roots are {}'.format(r))

coef = [8, -2, 16, -8, 1]
r = poly.polyroots(coef)
r_x =np.real(r[np.isreal(r)]) # try to print inside out to get the idea 
x = np.linspace(-1, 5, 100)
y = poly.polyval(x, coef)

plt.style.use('ggplot')
plt.plot(x, y, color = 'g', linewidth = 2)
plt.scatter(r_x, [0, 0], s = 20, c = 'r')
plt.grid(True)
plt.title('Root finding of a polynomial function')
plt.show()


範例 2:Numerical root-finding for the non-polynomial functions

  1. Solve f(x) = x - e^{-x/2} for x
  2. Draw the graph for each function and mark the locations of the roots.

Note:
  1. For non-Polynomial functions, it is usual to compute a root at a time in a bounded range. In this example, the sicpy.optimize.root_scalar is employed. Check with result of sicpy.optimize.root_scalar to see how to get the required root, e.g. print(dir(sol))
  2. The following codes use “plt.text” to mark the root position. To put a marker exactly on the position of the root, horizontal and vertical alignments are necessary.
import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt

def f(x) :
    return x - np.exp(-x/2)

# ---- First, draw the graph ------
x = np.linspace(0, 3, 1000)
plt.plot(x, f(x), color = 'g', linewidth = 1)
# ---- Second, compute the root ------
sol = opt.root_scalar(f, bracket=[0.5, 1], method='brentq')  
print('The root is at x = {:.4f}'.format(sol.root))  
# mark the root on top of function line
plt.text(sol.root, 0, 'X', color = 'r',
    horizontalalignment='center',
    verticalalignment='center')  
plt.title('$f(x) = x - e^{-x/2}$ and its root')
plt.xlabel('X'), plt.grid(True)
plt.show()
Another way to compute a root for the general function is as follows.
import scipy.optimize as opt

f = lambda x : x - np.exp(-x/2)
sol = opt.root(f, x0 = 1) 
# print(dir(sol))
print('The root is x = {:.4f}'.format(sol.x[0]))
print('The function value is f = {:.6f}'.format(sol.fun[0]))


範例 3:Numerical root-finding

solve the following equation
F(x) = \int_{-\infty}^{x} \frac{1}{\sqrt{2\pi}}e^{-t^2/2} \;dt = 0.95

Note:

  1. This is a root-finding problem by explicitly rewriting as f(x) = 0, where f(x) = F(x) - 0.95
  2. F(x) can be readily as the cumulative density function of the standard normal distribution.
  3. Before computing the root(s) of a complicated function, it is necessary to draw the graph of f(x) so as to get a rough idea where the root(s) are located.
import numpy as np
import scipy.optimize as opt
from scipy.stats import norm
import matplotlib.pyplot as plt
import scipy.integrate as integral


def f(x, prob) :
    g = lambda x : np.exp(-x**2/2) / np.sqrt(2*np.pi)
    return integral.quad(g, -np.inf, x)[0] - prob

# Draw function graph
prob = 0.95
x = np.linspace(-5, 5, 100)
vec_f = np.vectorize(f)
F = vec_f(x, prob)
plt.plot(x, F, linewidth = 2)

sol = opt.root_scalar(f, args = prob, bracket = [0, 3])   # pass prob as an extra argument
print('The root is x = {:.3f}'.format(sol.root)) 

plt.hlines(0, x.min(), x.max(), color = 'g')
plt.vlines(sol.root, F.min(), F.max(), color = 'g')
plt.grid(True)
plt.title('$F(x) - {:.2f}$'.format(prob))
plt.show()

# Draw the area in pdf 
y = norm.pdf(x)
plt.plot(x, y)
x_area = np.linspace(x.min(), sol.root, 100 )
y_area = norm.pdf(x_area)
plt.fill_between(x_area, y_area, 0, color = 'b', alpha = 0.3)
plt.text(0, 0.1, prob)
plt.title('$F(x)={:.2f}, x = {:.3f}$'.format(prob, sol.root))
plt.show()



練習:Numerical root-finding

Let the \beta pdf function be \beta(x;a,b).

  1. Solve \beta(x;2, 6) = \beta(x;4, 2).
  2. Draw a graph to indicate the position of the root.

範例 4:Symbolic root-finding

Solve each equation in a symbolic approach provided by sympy. Also, practice to draw the symbolic function and put marker on the root location.

  1. f(x) = x^2-3x+2 = 0
  2. f(x) = x - e^{-x/2} = 0

Note:
  1. The codes below use sympy.solve and sympy.solveset as an alternative approach.
  2. The result by sympy.solve is also in symbolic form. It needs to be transferred to floats for further use. It is a good technique to learn.
  3. The last part of the codes are commented for your reference. These codes demonstrate how to print symbolic contents.
import sympy as sym
from sympy import plot

x = sym.Symbol('x')
coef = [1, -3, 2] # from high to low
f = np.polyval(coef, x)
# f = x**2 - 3*x + 2
r = sym.solve(f)
# convert r to float numbers
r_float = [float(x) for x in r]

fig1 = plot(f, xlim = [-1, 4], ylim = [-1, 8], \
    show = False, title = 'Symbolic root-finding {}'.format(r_float))
fig2 = plot(0*x, xlim = [-1, 4], line_color = 'g', show = False)
fig1.append(fig2[0])
fig1.show()

# Alternative: sympy claims that solveset will replace solve in the future
rset = sym.solveset(f, x, domain = sym.S.Reals) #sym.Interval(0, 1)
rset_float = [float(x) for x in rset] # convert to float
print('The roots are {}'.format(rset_float))
#-------------------------------------------
# sym.pprint(f)
# exp = sym.sin(x) / sym.cos(x)
# sym.pprint(exp)
# sym.pprint(sym.simplify(exp))
#--------------------------------------------

範例 5:Numerical minimization

Find the minimum of f(x) = x^4 - 8x^3+16x^2-2x+8.

Note:

  1. Use scipy.optimize.minimize_scalar to solve for the local minimum of a function.
  2. The command minimize_scala can only compute a local minimum at a time. The location of the local minimum is determined by the attribute “bracket” which can contain two or three numbers. Refer to the document for the details. For example, for this function, bracket = [2, 3, 5], which means f(3)< f(2), and f(3) < f(5).
  3. The default algorithm used by scipy.optimize.minimize_scalar is “Brent”. For its details and other selections, please refer to the document.
  4. There is another local minimum except for the one shown on the graph. Please find it out.
import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt

# draw the graph
coef = [1, -8, 16, -2, 8] # high to low
f= lambda x : np.polyval(coef, x)
x = np.linspace(-1, 5, 1000)
plt.plot(x, f(x))
plt.xlabel('X'), plt.grid(True)

# compute the local minimum
res = opt.minimize_scalar(f, bracket=[2, 3, 5])
print(res) # find out the return values
print('The function has a local minimum at x = {:.4f}'.format(res.x))
print('The corresponding function value is {:.4f}'.format(res.fun))

plt.text(res.x, res.fun, 'X', color = 'r',
    horizontalalignment='center',
    verticalalignment='center')
plt.title('The local minimum of $f(x)$')
plt.show()


練習:Numerical minimization

\displaystyle \min_{x \in (0,3)} \tan^{-1}\frac{5}{x}+\tan^{-1}\frac{2}{3-x}

Note:

  1. You mat want to try another algorithm for scipy.optimize.minimize_scalar, e.g. res = opt.minimize_scalar(f, bounds = [0.1, 2], method = ‘bounded’), where the method “bounded” is employed and the attribute “bounds” is used to set up the searching boundary.

範例 6:Numerical minimization

Solve the following maximization problem.

\displaystyle \max_{0 \leq \theta \leq 3\pi} V(\theta),
where V(\theta)=\frac{\pi}{3}(\frac{10\pi - 5\theta}{2\pi})^2\sqrt{25 - (\frac{10\pi - 5\theta}{2\pi})^2}

Note:

  1. The maximization problem is considered as the minimization problem with a minus sign on the function, i.e., \max_x f(x) \equiv \min_x -f(x).
  2. This example demonstrates how to deal with a complicated function when a single line anonymous function is not appropriate.
import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt

# First, draw function graph
# define a function for a complicated funcion
def V(t) :
    tmp = ((10*np.pi - 5*t) / (2*np.pi))**2
    V = np.pi / 3 * tmp * np.sqrt(25 - tmp)
    return V

t = np.linspace(0, 3*np.pi, 1000)
plt.plot(t, V(t))

f = lambda x : -V(x)
res = opt.minimize_scalar(f, bounds = [0, 2])

plt.text(res.x, -res.fun, '+', color ='r', fontsize = 16,
    horizontalalignment='center',
    verticalalignment='center')
plt.xlabel('X'), plt.grid(True)
plt.title('The fnction maximum is {:.4f} at x = {:.4f}'.format(-res.fun, res.x))
plt.show()

習題:

  1. \displaystyle \min_x \sqrt{\frac{x^2+1}{x+1}}

  2. \displaystyle \min_{-4 \leq x \leq 3} (x+1)^5 sin(x-3)

  3. solve L(x) = 10 for x, where
    L(x) = \int_a^x \sqrt{1 + (f'(t))^2}\; dt and a = 0.

  4. \displaystyle \max_{\lambda} \ln \Pi_{i=1}^N f(x_i; \lambda) , where f(x_i; \lambda) denotes the likelihood function of the exponential distribution with parameter \lambda. Let the sample size N be 10, 20 ,30 50, 100, 300, 500, 1000. Generate the random sample x_i for each sample size and estimate \lambda by maximizing the log likelihood function \ln \Pi_{i=1}^N f(x_i; \lambda). Let the true \lambda = 0.5.