Bierlaire, M. (2015). Optimization: Principles and Algorithms. EPFL Press
def newtonFinDiffOneVariable(obj,x0,eps,tau,maxiter = 100):
xk = x0
# The function returns both f and g. Here, we need only f, so the [0] index.
f = obj(xk)[0]
k = 0
iters = list()
iters.append([k,xk,f])
while abs(f) > eps and k < maxiter:
if (abs(xk)>= 1):
s = tau * xk
else:
s = tau
fs = obj(xk+s)[0]
xk = xk - s * f / (fs - f)
f = obj(xk)[0]
k += 1
iters.append([k,xk,f])
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 $\tau=10^{-7}$
tau = 1.0e-7
(root,iters) = newtonFinDiffOneVariable(ex0703,2.0,1.0e-15,tau)
print("x= {} F(x)={}".format(root,ex0703(root)))
Table 8.1, page 204
print("k\txk\t\tF(xk)")
for k in iters:
print("{0}\t{1:+E}\t{2:+E}".format(*k))
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()
Run the example with $x_0=2$ and $\tau=0.1$
tau = 0.1
(root,iters) = newtonFinDiffOneVariable(ex0703,2.0,1.0e-15,tau)
print("x= {} F(x)={}".format(root,ex0703(root)))
Table 8.2, page 205
print("k\txk\t\tF(xk)")
for k in iters:
print("{0}\t{1:+E}\t{2:+E}".format(*k))
table= np.array(iters)
plt.xlabel("Iteration")
plt.ylabel("F(x)")
plt.plot(table[1:,2])
plt.show()
def secantOneVariable(obj,x0,a0,eps,maxiter=100):
xk = x0
ak = a0
f = obj(xk)[0]
k = 0
iters = list()
iters.append([k,xk,f,ak])
while abs(f) > eps and k < maxiter:
xold = xk
xk = xk - f / ak
fold = f
f = obj(xk)[0]
ak = (fold - f) / (xold - xk)
k += 1
iters.append([k,xk,f,ak])
return(xk,iters)
Example: $F(x)=x^2-2$, with $x_0=2$ and $a_0=1$
x0 = 2
a0 = 1
(root,iters) = secantOneVariable(ex0703,x0,a0,1.0e-15)
print("x= {} F(x)={}".format(root,ex0703(root)))
Table 8.3, page 207
print("k\txk\t\tF(xk)\t\tak")
for k in iters:
print("{0}\t{1:+E}\t{2:+E}\t{3:+E}".format(*k))
table= np.array(iters)
plt.xlabel("Iteration")
plt.ylabel("F(x)")
plt.plot(table[1:,2])
plt.show()
import numpy as np
from scipy import linalg
def newtonFinDiffNVariables(obj,x0,tau,eps,maxiter=100):
n = len(x0)
xk = x0
f = obj(xk)[0]
k = 0
iters = list()
iters.append([k,xk,f])
I = np.eye(n,n)
while linalg.norm(f) > eps and k < maxiter:
A = np.empty([n,n])
for col in range(n):
xi = xk.item(col)
if (abs(xi) >= 1):
s = tau * xi
elif xi >= 0:
s = tau
else:
s = -tau
ei = I[col].reshape(n,1)
fp = obj(xk + s * ei)[0]
A[:,col:col+1] = (fp - f) / s
# We need to enforce the column vector shape
d = linalg.solve(A,-f).reshape(2,1)
xk = xk + d
f = obj(xk)[0]
k += 1
iters.append([k,xk,f])
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 $\tau=10^{-7}$
# 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
x0 = np.array([[1,1]]).T
tau = 1.0e-7
(root,iters) = newtonFinDiffNVariables(ex0711,x0,tau,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 8.4, page 209
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)))
norm = [linalg.norm(iters[k][2]) for k in range(len(iters))]
plt.xlabel("Iteration")
plt.ylabel("||F(x)||")
plt.plot(norm)
plt.show()
Run the example with $x_0= \left(\begin{array}{c} 1 \\ 1 \end{array}\right)$ and $\tau=0.1$
# 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
x0 = np.array([[1,1]]).T
tau = 0.1
(root,iters) = newtonFinDiffNVariables(ex0711,x0,tau,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 8.5, page 210
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)))
norm = [linalg.norm(iters[k][2]) for k in range(len(iters))]
plt.xlabel("Iteration")
plt.ylabel("||F(x)||")
plt.plot(norm)
plt.show()
def secantNVariables(obj,x0,eps,maxiter=100):
n = len(x0)
xk = x0
f = obj(xk)[0]
k = 0
A = np.eye(n,n)
iters = list()
iters.append([k,xk,f,A])
while linalg.norm(f) > eps and k < maxiter:
xold = xk
d = linalg.solve(A,-f).reshape(2,1)
xk = xk + d
fold = f
f = obj(xk)[0]
y = f-fold
update = np.outer(y - A.dot(d), d ) / np.inner(d[:,0],d[:,0])
A = np.add(A,update)
k += 1
iters.append([k,xk,f,A])
return (xk,iters)
Example 7.11: $F(x)=\left(\begin{array}{c}(x_1+1)^2+ x_2^2 - 2 \\ e^{x_1} + x_2^3 - 2 \end{array}\right)$, $x_0=\left(\begin{array}{c} 1 \\ 1 \end{array}\right)$
# 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
x0 = np.array([[1,1]]).T
(root,iters) = secantNVariables(ex0711,x0,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 8.6, page 215
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)))
norm = [linalg.norm(iters[k][2]) for k in range(len(iters))]
plt.xlabel("Iteration")
plt.ylabel("||F(x)||")
plt.plot(norm)
plt.show()
Table 8.7, page 215.
for k in range(len(iters)):
x = iters[k][1]
A = iters[k][3]
J = ex0711(x)[1]
print("{}\t{:+E}\t{:+E}\t{:+E}\t{:+E}".format(k,J[0,0],J[0,1],A[0,0],A[0,1]))
print(" \t{:+E}\t{:+E}\t{:+E}\t{:+E}".format(J[1,0],J[1,1],A[1,0],A[1,1]))