Bierlaire, M. (2015). Optimization: Principles and Algorithms. EPFL Press.
The numbering of the algorithms, tables and page refer to the book.
Objective: find the machine epsilon
def machineEpsilon():
eps = 1.0
while eps + 1.0 != 1.0:
eps /= 2.0
return eps
eps = machineEpsilon()
print("Machine epsilon: {}".format(eps))
# The function obj return a tuple containing the value of the function and its derivative
def newtonOneVariable(obj,x0,eps,maxiter=100):
xk = x0
(f,g) = obj(xk)
k = 0
iters = list()
iters.append([k,xk,f,g])
while abs(f) > eps and k < maxiter:
xk = xk - f / g
(f,g) = obj(xk)
k += 1
iters.append([k,xk,f,g])
return (xk,iters)
Example: $F(x)=x^2-2$
def ex0703(x):
f = x * x - 2.0
g = 2.0 * x
return (f,g)
Run the example with $x_0=2$ and $\varepsilon=10^{-15}$.
(root,iters) = newtonOneVariable(ex0703,2.0,1.0e-15)
print("x= {} F(x)={}".format(root,ex0703(root)[0]))
Table 7.1 page 186
print("k\txk\t\tF(xk)\t\tF'(xk)")
for k in iters:
print("{0}\t{1:+E}\t{2:+E}\t{3:+E}".format(*k))
We plot the value of the function as the iterations progress, excluding the starting point
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
table= np.array(iters)
plt.xlabel("Iteration")
plt.ylabel("F(x)")
plt.plot(table[1:,2])
plt.show()
Another example: $F(x)=x-\sin(x)$
import numpy as np
def ex0704(x):
f = x - np.sin(x)
g = 1 - np.cos(x)
return (f,g)
Run the example with $x_0=1$ and $\varepsilon=10^{-15}$.
(root,iters) = newtonOneVariable(ex0704,1.0,1.0e-15)
print("x= {} F(x)={}".format(root,ex0704(root)[0]))
Table 7.2, p. 187.
print("k\txk\t\tF(xk)\t\tF'(xk)")
for k in iters:
print("{0}\t{1:+E}\t{2:+E}\t{3:+E}".format(*k))
We plot the value of the function as the iterations progress, excluding the starting point
table = np.array(iters)
plt.xlabel("Iteration")
plt.ylabel("F(x)")
plt.plot(table[1:,2])
plt.show()
Another example: $F(x)=\arctan(x)$
def ex0705(x):
f = np.arctan(x)
g = 1.0 / (1.0 + x * x)
return (f,g)
Run the example with $x_0=1.5$ and $\varepsilon=10^{-15}$. We set maxiter
to 10, as the algorithm is not converging.
(root,iters) = newtonOneVariable(ex0705,1.5,1.0e-15,maxiter=10)
print("x= {} F(x)={}".format(root,ex0705(root)[0]))
Table 7.3, page 188.
print("k\txk\t\tF(xk)\t\tF'(xk)")
for k in iters:
print("{0}\t{1:+E}\t{2:+E}\t{3:+E}".format(*k))
Plot the iterations, excluding the starting point. We observe the oscillations.
table = np.array(iters)
plt.xlabel("Iteration")
plt.ylabel("F(x)")
plt.plot(table[1:,2])
plt.show()
# The function f return a tuple containing the vector with the values of the functions
# and its Jacobian in numpy array format
from scipy import linalg
def newtonNVariables(obj,x0,eps,maxiter=100):
xk = x0
(f,J) = obj(xk)
k = 0
iters = list()
iters.append([k,xk,f,J])
while linalg.norm(f) > eps and k < maxiter:
# We need to enforce the column vector shape
d = linalg.solve(J,-f).reshape(2,1)
xk = xk + d
(f,J) = obj(xk)
k += 1
iters.append([k,xk,f,J])
return (xk,iters)
Example: $F(x)=\left(\begin{array}{c}(x_1+1)^2+ x_2^2 - 2 \\ e^{x_1} + x_2^3 - 2 \end{array}\right)$
# Note that the indices of array start at 0, not 1.
def ex0711(x):
f1 = (x.item(0) + 1) * (x.item(0) + 1) + x.item(1) * x.item(1) - 2
f2 = np.exp(x.item(0)) + x.item(1) * x.item(1) * x.item(1) - 2
# Note the double brackets, as the vector must be a matrix with one column
f = np.array([[f1],[f2]])
J11 = 2 * (x.item(0) + 1)
J12 = 2 * x.item(1)
J21 = np.exp(x.item(0))
J22 = 3 * x.item(1) * x.item(1)
J = np.array([ [ J11, J12 ],[J21 ,J22 ]])
return (f,J)
Run the example with $x_0= \left(\begin{array}{c} 1 \\ 1 \end{array}\right)$ and $\varepsilon = 10^{-15}$
# The .T represents transpose, to obtain a column vector. Note the double brackets syntax to emphasize
# that it is a matrix with only one column
x = np.array([[1,1]]).T
(root,iters) = newtonNVariables(ex0711,x,1.0e-15)
# Value at the solution
f = ex0711(root)[0]
print("x= ({},{}) F(x)=({},{})".format(root.item(0),root.item(1),f.item(0),f.item(1)))
Table 7.4, page 195.
print("k\txk\t\tF(xk)\t\t||F(xk)||")
for k in range(len(iters)):
print("{}\t{:+E}\t{:+E}\t{:+E}".format(iters[k][0],iters[k][1].item(0),iters[k][2].item(0),linalg.norm(iters[k][2])))
print(" \t{:+E}\t{:+E}".format(iters[k][1].item(1),iters[k][2].item(1)))
We plot the evolution of the norm of the function with the iterations.
norm = [linalg.norm(iters[k][2]) for k in range(len(iters))]
plt.xlabel("Iteration")
plt.ylabel("||F(x)||")
plt.plot(norm)
plt.show()
Example: $F(x)=\left(\begin{array}{c}x_1^3 - 3 x_1 x_2^2 -1 \\ x_2^3 - 3x_1^2 x_2\end{array}\right)$
# Note that the indices of array start at 0, not 1.
def ex0712(x):
# Note the double brackets, as the vector must be a matrix with one column
f = np.array([[x.item(0)**3 - 3 * x.item(0) * x.item(1)**2 - 1.0],[x.item(1)**3 - 3 * x.item(0)**2 * x.item(1)]])
J = np.array([ [ 3 * x.item(0)**2 - 3 * x.item(1)**2, -6 * x.item(0) * x.item(1) ],[-6 *x.item(0) * x.item(1) ,3* x.item(1)**2 - 3 * x.item(0)**2]])
return (f,J)
Run the example with $x_0= \left(\begin{array}{c} 1 \\ 1 \end{array}\right)$ and $\varepsilon = 10^{-15}$
# The .T represents transpose, to obtain a column vector. Note the double brackets syntax to emphasize
# that it is a matrix with only one column
x = np.array([[1,1]]).T
(root,iters) = newtonNVariables(ex0712,x,1.0e-15)
# Valur at the solution
f = ex0712(root)[0]
print("x= ({},{}) F(x)=({},{})".format(root.item(0),root.item(1),f.item(0),f.item(1)))
Run the example with $x_0= \left(\begin{array}{c} -1 \\ -1 \end{array}\right)$ and $\varepsilon = 10^{-15}$
# The .T represents transpose, to obtain a column vector. Note the double brackets syntax to emphasize
# that it is a matrix with only one column
x = np.array([[-1,-1]]).T
(root,iters) = newtonNVariables(ex0712,x,1.0e-15)
# Valur at the solution
f = ex0712(root)[0]
print("x= ({},{}) F(x)=({},{})".format(root.item(0),root.item(1),f.item(0),f.item(1)))
Run the example with $x_0= \left(\begin{array}{c} 0 \\ 1 \end{array}\right)$ and $\varepsilon = 10^{-15}$
# The .T represents transpose, to obtain a column vector. Note the double brackets syntax to emphasize
# that it is a matrix with only one column
x = np.array([[0,1]]).T
(root,iters) = newtonNVariables(ex0712,x,1.0e-15)
# Valur at the solution
f = ex0712(root)[0]
print("x= ({},{}) F(x)=({},{})".format(root.item(0),root.item(1),f.item(0),f.item(1)))